1230
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 8, AUGUST 2006
Impulse Control: Boolean Programming and Numerical Algorithms K. H. Kyung
Abstract—A numerical algorithmic approach to the impulse control problem is considered. Impulse controls are modelled by Boolean binary variables. The impulse Gâteaux derivatives for impulse times, impulse volumes and Boolean variables are derived, and these are applied to the numerical algorithms. These algorithms require significantly less computation time and memory storage than the quasi-variational inequalities by Bensoussan– Lions. By using our algorithms, complicated models of hybrid or constrained systems can be more easily treated numerically than by using Pontryagin’s Minimum Principle. Numerical experiments are performed for models on capacity expansion in a manufacturing plant, and on impulse control of Verhulst systems and Lotka–Volterra systems; the results confirm the effectiveness of the proposed method. Index Terms—Boolean programming, hybrid system, impulse control, impulse variational inequality.
I. INTRODUCTION
W
E consider the impulse control problem as follows:
(1) (2) (3) (4) (5) (6) have the binary value 1 or 0; where the Boolean variables is the prejump time; is the postjump time; for small ; . At each impulse time , the system is controlled impulsively with the impulse scale with its effect if . We suppose ; a fixed setup cost
is incurred at each impulse; the
Manuscript received September 9, 2003; revised September 21, 2004 and January 9, 2006. Recommended by Associate Editor Y. Wardi. The author is with the Department of Business Administration, Yonsei University, Seoul 120-749, Korea (e-mail:
[email protected]). Digital Object Identifier 10.1109/TAC.2006.879913
presence of prevents the system from being continuously conis the impulse control. trolled. The triplet There have been many approaches to solving impulse control problems. First, the problems have been treated using the celebrated quasi-variational inequalities by Bensoussan–Lions [3] and many others [18], [27]; this approach gives the global solution. Second, Pontryagin’s Minimum Principle has been applied to solve impulse control problems by [23]–[26], and others; this approach provides good insight into the necessary conditions for optimality. Third, impulse control problems are transformed into equivalent non-impulse control problems by introducing an artificial time to which Pontryagin’s Minimum Principle can be applied directly; this approach has often been applied in econometrics and management science since being developed by Vind [29], e.g., in [1], [12], [15]. Recently, there have been many important researches on related problems, [7]–[11], [19], [28], [31]. In the impulse control theory, there remain some unresolved difficulties. First, an enormous amount of computation time is required to solve the quasi-variational inequalities by Bensoussan–Lions [3] when there are many state variables, i.e., the problem of the “rise of dimension” [2] appears. The method uses unnecessary amount of memory to find a whole field of extremals on the value function, if we want only one optimal path from a known initial point. Second, the constrained or complicated models with impulses cannot easily be treated numerically by Pontryagin’s Minimum Principle, as seen in [12], [15], and [25]. In this paper, a numerical algorithmic approach is proposed: impulse decisions are modelled by Boolean binary variables , and binarities are transformed into equivalent constraints including continuous variables ; these conditions are treated using augmented Lagrangian methods, although they have been treated in references using the Lagrange function [17], the penalty approach [30] or proximal relaxations [16]. The Gâteaux derivatives for impulse times , impulse voland Boolean variables are derived, and these are umes applied to numerical algorithms using variational inequalities; the convergence of algorithms is proved. The algorithms are extended to hybrid systems as well as constrained impulse systems. The algorithms require significantly less computation time and memory storage than the quasi-variational inequalities by Bensoussan–Lions [3]. By using the proposed algorithms, complicated models with impulses can be more easily treated numerically than by using Pontryagin’s Minimum Principle [12], [23], [24], [29]. Numerical experiments are performed for models on capacity expansion in a manufacturing plant, and on impulse control
0018-9286/$20.00 © 2006 IEEE
KYUNG: IMPULSE CONTROL: BOOLEAN PROGRAMMING AND NUMERICAL ALGORITHMS
1231
of nonlinear Verhulst systems and Lotka–Volterra systems; the computational results confirm the effectiveness of the proposed method.
(14) (15) with vector
II. ASSUMPTIONS AND DEFINITIONS
. We define the adjoint , Hamiltonian functions for
belongs to the Sobolev space (16) (17) (18) and , defined later, also belong to . Assumptions follow. is continuous. A1): A2): . is continuous. A3): . A4): is continuous. A5): A6): . is continuous. A7): . A8): A9): is continuous. A10): . . A11): We assume that are strictly positive. We with assume the continuous differentiability of their arguments. Under assumptions A1)–A4), the existence and uniqueness of the solution of (2)-(4) follow [2], [19], [24]. We . use notations We assume the Gâteaux differentiability of and with re. The models with quadratic and linear spect to have many advantages in numerical algorithms. Variation Equations: We define the variables
(19) (20) with : Transpose of tonian. We use notations or ; but We assume that lutions .
: impulse Hamil. for so-
III. IMPULSIVE VARIATIONAL INEQUALITIES AND NUMERICAL ALGORITHMS A. Impulsive Variational Inequalities For problem (P0), we prove the following theorem. Theorem 1 (Impulse Gâteaux Derivatives): Assuming the with regard to in Gâteaux differentiability of problem (P0), we have the impulse Gâteaux derivatives i (21) ii (22) iii
satisfying the following differential equations for (the details are treated in Appendix A):
(23)
1
Proof: i) Gâteaux Differential : We can express the cost function (1) as a function of the impulse time for given (7) (8) (9)
(10) (11)
(24)
(12)
(13)
and [same results for for small similar to (134)]. In the case of the delayed impulse occurring at with delayed effects appearing at
1232
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 8, AUGUST 2006
, the function (24) is denoted as I): ref.(16), (19)
II): ref.(7)
II) : ref.(7)
(25) For small
, we compute now, from (24) and (25)
III): ref.(9),(17) III): ref.(9),(17)
By inserting the last results
with the definition and with notafor simplicity. The known tions is written as . Assuming the continuous differentiability of with regard to and with regard to , we have, denoting
into (26), we obtain
IV): ref.(8)
IV): ref.(8)
ref.(18)
(26)
Replacing by in (16), (19)] in the part of
V): ref.(17)
[the definition , we obtain V): ref.(17)
I): ref.(16), (19)
KYUNG: IMPULSE CONTROL: BOOLEAN PROGRAMMING AND NUMERICAL ALGORITHMS
The integration
ref.(19)
1233
can be transformed, using (16), as follows:
ref.(19)
I ref.(10) ref.(20) I ref.(10) From this, we have the Gâteaux differential with the gradient , and we obtain the result i). : ii). Gâteaux Differential Given , we can express the objective function (1) as a function of as follows:
II
ref(12)
II
ref(12)
From the last results of
, one obtains
ref.(18)
III ref.(11)
(27) We compute the and as ,
-derivative of
at a point
, writing
as III ref.(11)
From this, we have
Assuming the continuous differentiability of to and with regard to , we have
with regard
and we obtain the result ii). : iii). Gâteaux Differential , we express the objective function (1) as a Given function of the impulse variable
(28)
1234
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 8, AUGUST 2006
We compute the -derivative of at a point , writing as and as , assuming the continuous differentiability of with regard to
From this, we have the Gâteaux differential and the gradient . With these, we obtain the result iii). Lemma 1 (Impulse Variational Inequalities): Assuming that we have Gâteaux differentials i) and ii) in Theorem 1, we have the impulse variational inequalities i ii
The integration
Proof: i). Variational Inequality for : For optimal
, we have
For , we have ii). Variational Inequality for : For optimal
. , we have
can be transformed, using (16), as follows:
I ref.(13)
I ref.(13)
II II From the last results of
ref.(15) ref.(15) , one obtains the following:
ref.(18)
For , we have . Remark 1: (Representation of Impulse Phenomena and the Impulse-Time Derivative): Impulses and delayed impulses are represented in objective functions as (24) and (25), and in system equations as (135) and (136). These representations are in (21). helpful in deriving the impulse-time derivative Remark 2: (Necessary Conditions for Optimality): Condiand for are necessary for tions the optimality of the problem (P0); they can be derived also by Pontryagin’s Minimum Principle [6], [23], [25], [26]. Such necessary conditions, however, generate two-point boundary value problems, including impulses which are difficult to solve numerically. Remark 3 (Numerical Algorithms): The explicit forms of Gâteaux derivatives (21)–(23) in Theorem 1 are powerful in numerical implementation, especially in cases of complicated models or constrained problems; the algorithms require significantly less computation time and memory storage than the quasi-variational inequalities by Bensoussan–Lions [3]. in and Remark 4 (Models With ): If the state variables with multiple in (1)–(3), then Theorem jumps are included in impulses 1 and its proof become somewhat complicated. If one impulse is optimal, then Theorem 1 can be proved with some extensions of definitions and assumptions, e.g.,
III ref.(14) III ref.(14)
; then, similar results to those of Theorem 1 can be derived in practically the same way. In [3], [2], [18], and [27] on quasi-variational in , i.e., those models inequality, no models have . Such models have the form with no in are clear in their meaning and make the
KYUNG: IMPULSE CONTROL: BOOLEAN PROGRAMMING AND NUMERICAL ALGORITHMS
computation easy; many practical problems can be modelled as in . (1)–(3) with no
converges to the local solution metric convergence rates
1235
, for
, with geo-
B. Boolean Programming (35)
The Boolean binary conditions in (5) can be transformed into the equivalent quadratic equality conditions infor efficient algorithms: cluding continuous variables
(36) (37)
(29)
(38)
For the problem (1)–(4) and (29), we consider the dual problem subject to (2)–(4)
Proof: We examine the augmented Lagrangian function (32). The function (32) can be written by the second-order Taylor series as the following:
(30) with . In 0–1 programming, quadratic Boolean models were treated in [17], using a in (P); in [30], 0-1 problems were Lagrange function with treated using a penalty approach; in [16], set covering problems with binary variables were treated using proximal relaxations. We consider the problem (P) with the augmented Lagrangian for the efficient numerical solution. function By denoting , the saddle point of the problem (P) is such that
(31) (39) where denote optimal solutions. If is the saddle point satisfying (31), then is the solution of the original problem (P0). For the approximation of the saddle point , we use the coordinate gradient projection method (or point relaxation method, coordinate descent method [5], [13], [20]–[22]). For the problem (P), we consider the coordinate augmented Lagrangian function
where the relation is used from (29) for the local , and we use the relations solution . is the error term with orders The last term higher than second. We may find the next iterate near to for by (33) such that
(32) where the variables in (P) are assumed to be temin numerical iterations; we denote porarily known as as for simplicity. Theorem 2 (Boolean Programming): We assume that at the local solution of the . If we start from problem (P) and (32) for the known near to , then the successive sequence generated by the algorithm i).
(33)
ii).
(34)
(40) where
. Then, we have
(41)
1236
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 8, AUGUST 2006
We can write
From (P) and (32), one obtains (42)
(53)
(43) where Assuming first-order Taylor series, one obtains
by (34) and using the
for the solution in (29). The convergence condition for (51), (46) is, from (47), (53)
(44) where for the local solution and (42), one obtains
(54)
in (29). From (44)
for the solution in (29). We . With the , we have for , and
where (45) After subtracting one obtains
from use the definition: assumption
from both sides of the previous equation,
therefore, the inequality (54) is satisfied. and . Therefore, from (51), (46), (54), Therefore, algorithm i) and ii) converges. Remark 5: (Meaning of : Impulse Gain): We consider the capacity expansion problem (Example-1 in Section VI) with convex cost and linear system (46) (47)
From (42) and (44), one obtains (48)
(49) We assume the solution From (42) and (47)–(49), one obtains
has one or no impulse: . The system is
We can represent the objective function
(50) From (50), one obtains (51) (52)
as
KYUNG: IMPULSE CONTROL: BOOLEAN PROGRAMMING AND NUMERICAL ALGORITHMS
This can be rearranged for temporarily known
as
1237
, one obtains . We have two cases. and , then
optimality condition i) If
where
.
are defined as
, e.g., then
If . ii) If (we use
and instead of
, then for
): . If , e.g., then . . From i) and ii), we have In a different approach to the problem (P), we consider the following problems subject to (2)–(4): (P11)
We build an augmented Lagrangian function
(P12)
(P13)
We define the notion of “Impulse Gain”
as (P1)
where . In (P1), the term is added for concavifying with and for avoiding “ill-conditioning” in numerical regard to solutions. For (P1), we consider the coordinate Lagrangian with : known where the first term means costs caused by nonexpansion and the second term denotes the capacity expansion cost plus costs , one has after expansion. From (55) For the problem (P1), we prove the following theorem. Theorem 3: (Proximal Augmented Lagrangian): If the assumptions in Theorem 2 are satisfied, then the successive sequence, with (56) If
, then one obtains the condition . If , then one has the inequality . If , then it means that nonimpulse costs are greater than impulse costs and, therefore, capacity expansion (impulse) is the optimal decision, i.e., , and we denote as . If , , and we denote as . From the then
(57) (58) (59) converges to the local solution
of (P1).
1238
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 8, AUGUST 2006
Proof: iii) and iv): (Transformation): as follows: We define
(60)
If we minimize :
with respect to
, then we obtain iii) with
(61) Inserting
into
, we have
(62) Fig. 1. Saddle-points and contour-lines.
From (60), (62), and by introducing
, we obtain
(63) If we maximize
with respect to
, then we obtain iv) with
Remark 7: (Choose of and ): In order to fulfill the convergence condition in (37), one must know the local solution . At the beginning of the iterations, however, they are not known. help to satisfy One can guess that large enough values of creates the convergence condition. However, an overly large the ill-conditioning problem in solving (33). We consider an example (Fig. 1)
(64) Using
in (64), we have from (60) and (63) with
(65) Therefore, problem (P1) or (55) is transformed into the problem (P) or (32). i), ii): (Convergence): From (65) and (32), the proof of i), ii) is the same as that of Theorem 2. Remark 6: (Concavifying and avoiding ill-posedness; relations between , and ): From in (63), one , obtains the following relations: If , i.e., if we use , then it means that we use then . Therefore, using for concavifying in (P1) small has the same effect as choosing small in iteration (33),(34) for avoiding the ill-conditioning.
; 2-cases are depicted for has the saddle point ; anis non-optimal. If other saddle point the large is selected at the beginning of the iterafor approximates tion, then the iterative solution in this example. Therepossibly to the false solution fore, it is good to start with a small value of . After ill-conditioning problems are avoided in iterative solutions (i.e., for in the last figure), the large helps the iterain tive solution to approximate quickly to the solution with each iteration or keep this example. One may increase as constant. For example, in CASE-3 described later in numerical examples, was chosen; it was increased with until . C. Numerical Algorithms and Convergence We apply Theorems 1 and 2 to the solution of the impulse control problems (P0) and (P). For the approximation of the saddle point satisfying (31), we use the coordinate gradient projecusing impulse Gâteaux derivatives detion algorithm for rived in Theorem 1 and the augmented Lagrangian method for discussed in Theorem 2.
KYUNG: IMPULSE CONTROL: BOOLEAN PROGRAMMING AND NUMERICAL ALGORITHMS
Algorithm I: The following procedures i), ii), and iii) are repeated at each iteration step . i) One determines the impulse variables
1239
Algorithm II-2: Instead of Algorithm II-1, one may use the following iterations: (74) (75)
for with temporarily known the iterative algorithm
(66)
(67) (68) (69) ii) Find the state for the given i.e., integrate (2) with (3) and (4):
(76)
by
,
where and the iterations are equal to Algorithm I. Other iterations are similar to for Algorithm I. Algorithm II-2 uses the alternating iteration from to , which helps the convergence in some cases. We consider the coordinate gradient projection method (66), ; we use the notation (67): instead of for simplicity; or denotes the onto for with projection of . The projection operator is a contraction operator: . is assumed to be -convex (77) Also, we assume that
satisfies the Lipschitz condition
iii) Integrate the adjoint system (16) with (17) and (18)
(78)
In i), the derivatives are defined in Theorem 1; we use from (P); means the projection on the in ii) and iii) boundary for . The integration for can be approximately obtained using numerical methods, e.g., Euler Method or Runge–Kutta method. Here, we apply the Arrow–Hurwicz-type algorithm (68) for preventing ill-conditioning in iterative solutions, although the algorithms (33) in can be a small Theorem 2 are Uzawa-type. The step size constant or selected by step size rules [5], [22]. Algorithm II-1: We use Theorem 3 with some revisions
Then, we have the following lemma. Lemma 2 (Convergence Condition of the Coordinate ): If of (P0) is Gradient Projection Method for -convex and if satisfies the Lipschitz condition with , then the gradient method with projection (66) converges of for coordinate by coordinate to the local minimum . Similarly, if is -convex and satisfies the Lipschitz condition with , then the if gradient method with projection (67) converges coordinate-wise of for . to the local minimum Proof: The local solution is characterized by the following variational inequality [ref. i) in Lemma 1]: (79) which can be written by
(70)
(80)
(71) (72)
with
, and by the definition of projection, we have
(73) where the iterations for are equal to Algorithm I. In each iteration, we find the variables in (2)–(4), (16)–(18) in (70), (72) similar to Algorithm I. The iterations of have the Arrow–Hurwicz type for preventing ill-conditioning. By using Algorithm II-1 with adequate and in in (55), one obtains solutions more easily in some cases than by using Algorithm I. In cases of sensitive models, one often obtains many near-optimal solutions, the cost function values of which are not much different from that of the optimal solution; one can compare such solutions with the results obtained by Algorithm I.
(81) We can therefore write
(82) Since the projection operator is a contraction, we have
(83)
1240
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 8, AUGUST 2006
Inequality (83) can be written
, calculating constructed in terms of the vector in its components in turn through the minimization of in the same way as (66)–(68): one variable (93) (84)
with
Because is -convex (77), and because satisfies the Lipschitz condition with (78), we have (85) Therefore, for
, we have (86) (87)
which proves the convergence of the algorithm (66) for . using Similar to , we have the convergence results for the algorithm: ; here, denotes the coordinate-wise projection of the onto for with :
The projection means the same as (66)–(68); we use the noinstead of for known for simplicity. tation , the minimization of in (30) is equivWith respect to alent to the minimization of in (P0). as follows: We define for , where are used as the same -convex in (87), (89), (37), and (38). We assume that is -convex, and and satisfy the Lipschitz and conditions for all and as in Lemma 2. We assume that near to the solution as in Theorem 2. From the convergence of in (86), satisfying (87), (89), (88), (35) with step-sizes . and (37), we have : We introduce the intermediate vector
(88) (89) which proves the convergence of (67). We prove the convergence of the algorithms (66)–(69), i.e., the coordinate gradient projection algorithm (66), (67) and the augmented Lagrangian method (68), (69). Theorem 4 (Convergence of the Coordinate Gradient Projection Method and the Augmented Lagrangian Method): If the assumptions in Lemma 2 and Theorem 2 are satisfied, then the successive sequence , generated by the algorithm converges to the local solu(66)–(69) starting from tion of the problem (30) or (31) with geometric convergence rates
(94) From the definitions (93), (94), (86), (88), and (35), we have, , the following result: with
(95) (90) (91) where near to . Proof: The control variable
From the inequality (95), one obtains the following result:
with can be written in detail as
(92) with . in (30) with respect to For the minimization of for given , the successive descent directions are chosen in the directions of the coordinate axes, taken in cyclic order. Starting is with an initial vector , each vector
with , which proves the convergence of (90). In a similar way, we consider the convergence (91) for with as (36) and (37). The variable is calculated in the same way as (69) or (34) (96) for
. We introduce the intermediate vector (97)
KYUNG: IMPULSE CONTROL: BOOLEAN PROGRAMMING AND NUMERICAL ALGORITHMS
with . From the convergence of in (36) with satisfying (37), we have the following result:
(98) From (98), we have the following inequality:
with , which proves the convergence of (91). converging to the solution Then, one obtains satisfying the saddle point condition (31). Remark 8: (Convexity and Proximity): If in (1) is convex for , and in (2) is linear for , then the coordinate gradient projection algorithm (66)–(69) converges efficiently to the sois not convex for , then the proximal terms for lution. If can be added to (refer to [4]): . The convexified helps for the iterative algorithm to converge quickly to the global solution in some cases. The convexity for was treated with an example in Remark 5. Remark 9: (Good Convergence in Using Boolean Variables ): If one omits the Boolean variables in the model (P0) (1)–(6) (e.g., if we set in the Example-1 CASE 1-1), then the model becomes somewhat simple, and, in some cases, ; one can use gradient algorithms using Theorem 1 with but, if we omit in models, then the convergence of algorithms in becomes unstable compared to our algorithms containing models, and the optimality of the impulse decision (1-impulse or 0-impulse) is difficult to confirm. IV. HYBRID CONTROL
1241
satisfying the adjoint (16)–(18), but with different Hamiltonian functions including continuous control , as the following:
For numerical solution, we use the coordinate gradient methods. in iterations, For temporarily known impulse variables . We prove we try to find the optimal continuous control the following lemma. and Variational Lemma 3: (Gâteaux Differential Inequality): Assuming the Gâteaux differentiability of with respect to , we have the following Gâteaux derivative and vari: ational inequality, for given i). ii). Proof: See Appendix B. , we see that is linear For and is the continuous with regard to gradient of at the point . Algorithm: Procedures of the coordinate gradient method for hybrid control follow. using I)-i). One determines the impulse variables the same algorithms in the previous section, for temporarily at each iteration step . known ii). Find the trajectory , for known , i.e., integrate the state equation , . iii) Integrate the adjoint equation, for known
. and II). For known impulse controls from step I), one determines the continuous control the gradient method as follows. is computed by i). The gradient of at
by
We now consider the control problem of hybrid systems in addition to the which contain the continuous control in (P0): impulse control ii).
is obtained in the direction
as follows:
(99)
(103)
(100) (101) (102)
The iterations are repeated starting with the new control . The integration for can be approximately obtained by using numerical methods. In case of coordinate-wise constraints, we use the gradient projection method, for
where . In the previous sections, the did not appear in the model. In this seccontinuous control tion, we consider the hybrid systems containing discrete impulse and continuous ordinary control control at the same time. We introduce the adjoint variable ,
(104) where boundary
denotes the projection of .
onto the
1242
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 8, AUGUST 2006
V. CONSTRAINED IMPULSE CONTROL We consider the constrained impulse control problem as the following:
, similar to the problem P0),(1)–(5) or to the hybrid control (99)–(102). We define the Hamiltonian functions and adjoint systems (117) (118)
(105) (106)
(119) (120)
(107) (108)
In order to find the saddle-point of the problem , we use the coordinate gradient projection algorithm similar to Algorithm I in the Section III-C
(109) (110) (111) (112)
(121)
The inequality constraints (111), (112) can be transformed into equality constraints as
(123)
(122)
(124) (125)
(113) (114)
(126) For (105)–(110) and (113)–(114), we formulate the following augmented Lagrangian: (127)
(128)
(115)
(129) where defined as the following:
, and
are where projection, and and jectories state equation
Corresponding to the constrained impulse control problem (115), (106)–(110), we consider the following saddle-point problem:
, and
for new
means the coordinate . At each step , find tra, i.e., integrate the
integrate the adjoint equations
(116) To the problem , (115), (106)–(110), we apply the coordinate-wise gradient projection methods with
The convergence of the algorithm can be analyzed by methods similar to those used in previous sections.
KYUNG: IMPULSE CONTROL: BOOLEAN PROGRAMMING AND NUMERICAL ALGORITHMS
1243
VI. NUMERICAL EXAMPLES Example-1: Impulsive Capacity Expansion: (CASE 1-1): A capacity expansion problem in a manufacturing plant is formulated as:
Fig. 2. Impulse control Example 1-1.
where denotes the production capacity, : The capacity expansion time in the planning period : The capacity expansion scale, : The decision to expand or not to expand, : The forecasted sales of the product, : The depreciation rate of the production capacity, : The discount rate, : the cost coefficients. The problem is to find the capacity expansion time and scale that will satisfy the increasing sales. The real cost data are approximated by a quadratic function (see, e.g., [14]). The following data are given. . The optimal solution is One impulse is optimal. (CASE 1-2): The jump cost is equal to that of (CASE 1-1). The optimal solution is
,
. . Another set of data is
. Two jumps are optimal. The results of (CASE 1-1),(CASE 1-2) are depicted in Figs. 2 and 3. Example-2: Linear-System Quadratic-Costs With Constraints: We consider the constrained impulse control problem as the following:
(CASE 2-1). No constraints.
Fig. 3. Impulse control Example 1-2.
The following example for the capacity expansion in a plant is assumed to be given:
where denotes the production capacity of the equipment, : The product-inventory level, : The pro: The forecasted sales of the product, : duction rate, : The capacity expansion The capacity expansion time, : The scale, : The decision to expand capacity or not,
1244
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 8, AUGUST 2006
Fig. 4. Impulse control Example 2-1.
Fig. 5. Impulse control Example 2-2.
depreciation rate of the production capacity, : The depreciation rate of the product inventory, : The discount rate, : The cost coefficients. In contrast to (CASE 1-1), the product inventory and production rate are included in the model in order to compensate for the shortage in or excess of capacity. The problem is to find the capacity expansion time and scale that will satisfy the increasing demand while taking into account the inventory and production are rate. There are no constraints, i.e., . The following equal to 0, and data are assumed to be given:
Example-3: Harvesting of Nonlinear Verhulst-System: (CASE 3): A harvesting problem of a resource population (e.g., forest, fish, rabbit, etc.) is formulated as: Minimize
(130)
The optimal solution is . One impulse is optimal. (CASE 2-2). Constraints: In addition to formulas in (CASE 2-1), the following constraints are assumed to be given: . For constraints, we define vectors and matrices
There are constraints on , i.e., . The is half of that found in (CASE 2-1). Another jump cost set of data is equal to that of (CASE 2-1). The optimal solution is . Two jumps are optimal. The results of (CASE 2-1) and (CASE 2-2) are depicted in Figs. 4 and 5.
where denotes the population density, : The harvesting : The harvesting scale, : The time during the period decision to harvest or not to harvest, : The objective density of the population, : The system parameter, the carrying : The cost coefficients. capacity, : The discount rate, The problem describes what needs to be solved for the system to preserve the objective state of the population density with minimum harvesting costs. The nonlinear Verhulst model (130) describes the population dynamics under environmental limitations. The serves as the growth rate. The term population leads to an exponential growth but, as time goes on and the population reaches larger values, the quadratic becomes important and produces a term decrease in the rate of growth. The exact solution of the nonlinear differential equation (130) for one jump with is: for , and at the impulse appears: . After impulse
, the variable evolves as: for
, and so on until . The real cost data are approximated by a quadratic function as in CASE-1. The following data are assumed to be given: .
KYUNG: IMPULSE CONTROL: BOOLEAN PROGRAMMING AND NUMERICAL ALGORITHMS
1245
Fig. 7. Impulse control Example 4. Fig. 6. Impulse control Example 3.
The optimal solution is: . Two impulses are optimal. The results of (CASE 3) are depicted in Fig. 6. For large seems to be linear. Example-4: Lotka–Volterra Prey–Predator Systems: (CASE 4): An impulse control problem of a Lotka–Volterra system (e.g., management of a tourist aquarium with a prey-predator system of small fish and big fish) is
(131)
from tourists, or he plans to fish the excess fishes after considering the different increasing rates of the prey–predators. The real cost data are approximated by a quadratic function as in CASE-1. The following data are assumed to be given:
. The optimal solution is . One impulse is optimal. One solves approximately the nonlinear differential equations (131) and (132) and adjoint equations as
From many test results, it can be seen that the previous Taylor series method is much better for solving the previous nonlinear chaotic systems than either the Runge–Kutta method or the Euler method. Some detailed techniques are needed in numerical implementation. The results of (CASE 4) are depicted in Fig. 7.
(132) VII. CONCLUDING REMARKS where denote the population density of small fish and big fish, : The fishing time during the period : The fishing scale, : The decision to fish or : not to fish, : The objective density of the population, The objective scales of the impulses, : The discount rate, : The cost coefficients. The nonlinear differential equations (131) and (132) have been introduced by Volterra to model the dynamics of a biological system composed of two populations: prey and predator. If the system is in a different state from the stable or objective , then the system evolves very sensitively to the pastates , and the system often reaches unstable rameters chaotic states. The aquarium manager waits for the system to reach the objective state in order to obtain maximum profits
In this paper, first, impulse decisions are modelled by Boolean , and binarities are transformed binary variables , including contininto equivalent constraints uous variables ; these conditions are treated using augmented Lagrangian methods. Second, the Gâteaux derivatives for imand Boolean variables pulse times , impulse volumes are derived. Third, the Gâteaux derivatives are used in the numerical algorithm “coordinate gradient projection method,” and the convergence of the algorithms is proved. The algorithms are extended to hybrid systems and constrained impulse systems. Numerical experiments are performed for models on capacity expansion in a manufacturing plant, and on impulse control of nonlinear Verhulst systems and Lotka–Volterra systems. In cases with few jumps, i.e., one or two jumps, the solutions to impulse control problems can be found very efficiently by
1246
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 8, AUGUST 2006
our algorithms, assuming the convexity of the objective functions and the linearity of the systems. Using our algorithms, accuracy is greatly increased, and computational requirements are significantly reduced compared to the quasi-variational inequality method by Bensoussan–Lions [3], which uses unnecessary amount of memory to find a whole field of extremals on the value function. Also, by using our algorithms, the complicated impulse control models with constraints can be more easily treated numerically than by using Pontryagin’s Minimum Principle, which yields two-point boundary value problems including impulses. If the optimal solution has multiple jumps or if one treats hybrid or constrained systems, then one often obtains multiple near-optimal solutions. For checking the uniqueness of the solution, the near-optimal solutions can be compared with the approximate global solution which may be obtained by the quasivariational inequalities by Bensoussan–Lions [3]. The global convexity for all impulse variables is difficult to analyze. Further research is necessary to investigate such a subject. APPENDIX A. Variations in Variation in th impulse for small as
: We denote the times related to the , and
and we have . Hence,
satisfies
(137) It follows that 1 system for
satisfies the differential
(138) (139) (140) Variation in : The mapping which associates with the solution of the state equation is assumed to be Gâteaux-differentiable with regard to appearing at , i.e., for exists, and . satisfies
(133) if
(141) (134)
and its effect apWe assume that the th impulse occurs at pears at ; the delayed impulse occurs at and its effect for . (For , the appears at and the delayed effect apdelayed impulse occurs at as (134)). pears at If the th impulse occurs at and its effect appears at , then for as ; we denote the state , it assuming the definition of impulse times (133) for satisfies
It follows that 1 with given
satisfies the differential system for )
(142) (143) (144) Variation in as
: We can define
similar to
(135) can be treated using in (134) with small (The case and its changes). If the delayed th impulse occurs at , then we denote the state as effect appears at
(136) If we assume the Gâteaux-differentiability of with regard to and the differentiability of with regard to , then we define
for . It follows that satisfies the difwith given ferential system for
(145) (146) (147)
KYUNG: IMPULSE CONTROL: BOOLEAN PROGRAMMING AND NUMERICAL ALGORITHMS
B. Proof of Lemma 3 in Hybrid Control Variation Equations: The mapping, which associates the solution with of the state (100), is assumed to be continuous from to , and Gâteaux-differentiable with regard to , i.e., , exists and is denoted by , with . satisfies the equation
(148) If we assume to be continuously differentiable with regard to and , then satisfies for
(149) (150) where and are computed for . (Proof of Lemma 3): : Proof: i). Gâteaux Differential Assuming that is continuously differentiable with regard to and , we have
where is used for . The rest of the proof is quite similar to the proof in the ordinary control theory, for given . Replacing with [ref. (16)], we obtain
with the relations:
; we used (149), (150), and
. We obtain result i).
1247
ii). Variational Inequality for : From i) for optimal , we have
For result ii).
, we have
. We obtain
REFERENCES [1] K. J. Arrow and M. Kurz, Public Investment, The Rate of Return, and Optimal Fiscal Policy. Baltimore, MD: Johns Hopkins Univ. Press, 1970. [2] M. Bardi and I. Capuzzo-Dolcetta, Optimal Control and Viscosity Solutions of Hamilton-Jacobi-Bellman Equations. Boston, MA: Birkäuser, 1997. [3] A. Bensoussan and J.-L. Lions, Impulse Control and Quasi-Variational Inequalities. Paris, France: Bordas, 1984. [4] D. P. Bertsekas, “Convexification procedures and decomposition algorithms for large-scale nonconvex optimization problems,” J. Optim. Theory Appl., vol. 29, pp. 169–197, 1979. [5] D. P. Bertsekas, Constrained Optimization and Lagrange Multiplier Methods. San Diego, CA: Academic, 1982. [6] A. Blaquière, “Impulsive optimal control with finite or infinite time horizon,” J. Optim. Theory Appl., vol. 46, no. 4, pp. 431–439, 1985. [7] M. S. Branicky, V. S. Borkar, and S. K. Mitter, “A unified framework for hybrid control: Model and optimal control theory,” IEEE Trans. Autom. Control, vol. 43, no. 1, pp. 31–45, Jan. 1998. [8] A. Bressan and F. Rampazzo, “On differential systems with quadratic impulses and their applications to Lagrangian mechanics,” SIAM J. Control Optim., vol. 31, pp. 1205–1220, 1993. [9] F. Camilli and M. Falcone, “Analysis and approximation of the infinitehorizon problem with impulsive controls,” Automat. Remote Control, vol. 58, no. 7, pp. 1203–1215, 1997. [10] C. G. Cassandras, D. L. Pepyne, and Y. Wardi, “Optimal control of a class of hybrid systems,” IEEE Trans. Autom. Control, vol. 46, no. 3, pp. 398–415, Mar. 2001. [11] B. De Schutter, “Optimal control of a class of linear hybrid systems with saturation,” SIAM J. Control Optim., vol. 39, no. 3, pp. 835–851, 2000. [12] J. R. Dorroh and G. Ferreyra, “A multistate, multicontrol problem with unbounded controls,” SIAM J. Control Optim., vol. 32, pp. 1322–1331, 1994. [13] R. Glowinski, J.-L. Lions, and R. Trémolières, Numerical Analysis of Variational Inequalities. Amsterdam, The Netherlands: North-Holland, 1981. [14] C. C. Holt, F. Modigliani, J. Muth, and H. A. Simon, Planning Production, Inventories and Workforces. Englewood Cliffs, NJ: Prentice-Hall, 1960. [15] M. I. Kamien and N. L. Schwartz, Dynamic Optimization: The Calculus of Variations and Optimal Control in Economics and Management. New York: North-Holland, 1981. [16] K. C. Kiwiel, P. O. Lindberg, and A. Nõu, “Bregman proximal relaxation of large-scale 0-1 problems,” Comput. Optim. Appl., vol. 15, pp. 33–44, 2000. [17] P. F. Körner, “An efficient method for obtaining sharp bounds for nonlinear Boolean programming problems,” RAIRO Recherche Opérationnelle, vol. 26, pp. 113–120, 1992. [18] P. L. Lions and B. Perthame, “Quasi-variational inequalities and ergodic impulse control,” SIAM J. Control Optim., vol. 24, pp. 604–615, 1986. [19] B. M. Miller, “The generalized solutions of nonlinear optimization problems with impulse control,” SIAM J. Control Optim., vol. 34, pp. 1420–1440, 1996. [20] S. Mitter, “Successive approximation methods for the solution of optimal control problems,” Automatica, vol. 3, pp. 135–149, 1966. [21] E. Polak, Optimization: Algorithms and Consistent Approximations. New York: Springer-Verlag, 1997. [22] B. T. Polyak, Introduction to Optimization. New York: Optimization Software, Inc., 1987. [23] R. Rempala and J. Zabczyk, “On the maximum principle for deterministic impulse control problems,” J. Optim. Theory Appl., vol. 59, pp. 281–288, 1988. [24] A. M. Samoilenko and N. A. Perestyuk, Impulsive Differential Equations. Singapore: World Scientific, 1997.
1248
[25] A. Seierstad and K. Sydsaeter, Optimal Control Theory with Economic Applications. Amsterdam, The Netherlands: North-Holland, 1987. [26] S. P. Sethi and G. L. Thompson, Optimal Control Theory: Applications to Management Science. Boston, MA: Martinus Nijhoff, 1981. [27] A. Sulem, “Quasi-variational inequalities and impulse control problems,” in Optimal Pricing, Inflation, and The Cost of Price Adjustment, E. Sheshinski and Y. Weiss, Eds. Cambridge, MA: MIT Press, 1993, pp. 57–96. [28] A. J. van der Schaft and J. M. Schumacher, “Complementarity modeling of hybrid systems,” IEEE Trans. Autom. Control, vol. 43, no. 3, pp. 483–490, Mar. 1998. [29] K. Vind, “Control systems with jumps in the state variables,” Econometrica, vol. 35, no. 2, pp. 273–277, 1967. [30] S. Wang, K. L. Teo, H. W. J. Lee, and L. Caccetta, “Solving 0-1 programming problems by a penalty approach,” Opsearch (Oper. Res. Soc. India), vol. 34, no. 3, pp. 196–206, 1997. [31] S. T. Zavalishchin and A. N. Sesekin, Dynamic Impulse Systems. Dordrecht, The Netherlands: Kluwer, 1997.
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 8, AUGUST 2006
K. H. Kyung received the Bachelor of Commerce and M.B.A. degress from Yonsei University, Seoul, Korea, and the Dr.oec. degree from Hochschule St. Gallen für Wirtschats und Sozialwissenschaften (now Universität St. Gallen), Switzerland. He is currently a Professor of Operations Research and Applied Mathematics in the Department of Business Administration at Yonsei University. His research experiences include: Visiting Scholar at ENSMP, Centre Automatique et Systeme (France); Humboldt-Universität zu Berlin, Institut für Physik (Germany); Department of Electrical Engineering and Computer Sciences, University of California at Berkeley; and the Laboratory for Manufacturing and Productivity, Massachusetts Institute of Technology, Cambridge. His research interests include impulse control, hybrid systems, optimal control theory, dynamic programming, stochastic control with applications to telecommunication networks, neural networks, transportation, and production systems. Dr. Kyung was a recipient of the Swiss Federal Scholarship.