Ind. Eng. Chem. Res. 2000, 39, 2305-2314
2305
Numerical Solution of Dynamic Optimization Problems with Flexible Inequality Constraints by Adaptive Stochastic Algorithm Cheng-Liang Chen* and Daim-Yuang Sun Department of Chemical Engineering, National Taiwan University, Taipei 10617, Taiwan, R.O.C.
An algorithm for solving nonlinear dynamic optimization problems subjected to flexible path constraints is proposed. By using the concept of fuzzy set to define the degree of acceptability for path constraints and degree of satisfaction for objective function, the original problem is transformed into a fuzzy decision making problem. To determine the optimum of the transformed system, fuzzy inference is incorporated into the computational platform, the integrated controlled random search for dynamic system (ICRS/DS). Because of the intrinsic structure, such an algorithm is very straightforward and efficient. Three examples are employed to demonstrate the simplicity and capability for the proposed algorithm. 1 Introduction 1960s,1,2
the dynamic optimization of chemiSince the cal processes has received considerable attention. Several solution methods have been proposed in this field. Those methods can be classified into three major categories, including the following: (1) variational formulation (maximum principle of Pontryagin); (2) transformation techniques, which may be further grouped into complete parametrization3,4 and control vector parametrization;5-7 (3) iterative dynamic programming (IDP).8 Moreover, all these researches assume the system being operated under deterministic environment. However, real situations confronted in operating chemical processes are sometimes not so rigid. For example, engineers may subjectively allow the temperature of a batch reactor to rise some degrees to attain higher conversion, provided there appear to be no imminent risks. In optimizing such a system, one must take the higher temperature into account. Despite successful progress in solving most deterministic dynamic optimization problems, all the aforementioned methods fail to be implementable, as flexible factors appear in the system. In general, the recognition for the flexible factors in operating chemical processes is subjective and qualitative. Therefore, the language of fuzzy sets developed by Zadeh9 is a suitable tool on describing subjective tendency. By using [0, 1] multivalued logics, which is in contrast to traditional {0, 1} binary logics, subjective and qualitative properties can be quantitatively depicted. Bellman and Zadeh further extended fuzzy set theory to the field of decision making.10 In the past, research and applications for fuzzy steady-state optimization have been fruitful.11 In this study, the authors attempt to seek an integration of fuzzy set theory and optimization technique to resolve the dynamic optimization problems for nonlinear chemical processes restricted by a flexible environment. We consider a flexible dynamic optimization problem as a dynamic fuzzy decision making problem. After the objective function and flexible constraints are graded by the language of fuzzy set and the continuous control * Corresponding author: E-mail:
[email protected]. Telephone: 886-2-23636194. Fax: 886-2-23623040.
input is discretized, the fuzzy inference can be easily incorporated into a computational platform, the socalled integrated controlled random search for dynamic system (ICRS/DS),12 which is an efficient and powerful dynamic optimization technique, to solve the problem. The rest of this paper is organized as follows: Section 2 introduces the formulation of the problems including the fuzzification of flexible constraints and objective function. Section 2.3 reviews fuzzy optimization theory. Section 3 presents the solution algorithm, a modification of ICRS/DS. Section 4 provides the computational illustrations of the proposed algorithm. Finally, conclusions are made in section 5. 2. Problem Formulation Consider the following dynamic optimization problem
max J[x(tf)]
u(t)∈Ω ˜
(1)
where J[x(tf)] denotes the objective function to be maximized and Ω ˜ ≡ {u(t) | h(x(t), u(t)) ) x3 (t) - f(x(t), ˜ ; u(t) e u(t) e u j (t)} u(t)) ) 0, x(0) ) x0; g(x(t), u(t)) e b represents the feasible control input space, in which all control policies u(t) satisfy system dynamics, x3 (t) ) f(x(t), u(t)), and flexible inequality constraints, g(x(t), u(t)) e b ˜ . Here, x(t) is an (n × 1) state vector with x0 as the initial value; u(t) is an (m × 1) control vector ˜ ) bounded by u(t) and u j (t); g ≡ [g1, ..., gK]′, and b [b˜ 1, ..., b˜ K]′. The flexible inequality constraints imply that, as the reasonable limiting value, bk, and its acceptable maximal tolerance, pk, for each constraint can be preliminarily assigned, the values of gk(x(t), u(t)) located in (bk, bk + pk) all can be regarded as satisfying the constraint. In this study, the final time end tf is specified in advance. Notably, the type of objective function in our problem is very general. Other types of objective functions can always be cast into this form by defining an additional state variable. Moreover, there is no loss in generality in maximizing, since an expression can be minimized by simply changing the sign of the objective function. 2.1. Fuzzification of Flexible Constraints. For each inequality level gk(x(t), u(t)) between bk and bk + pk, the extent of relaxation accepted by the decision
10.1021/ie990638n CCC: $19.00 © 2000 American Chemical Society Published on Web 06/16/2000
2306
Ind. Eng. Chem. Res., Vol. 39, No. 7, 2000
with the following degree of acceptability µC
µC ) T{µC1, ..., µCk}
Figure 1. Linear type monotonic decreasing membership function for flexible inequality constraints.
Figure 2. Exponential type monotonic decreasing membership function for flexible inequality constraints.
maker (DM) can be variable. It takes for granted that all values for gk(x(t), u(t)) less than bk can be thoroughly accepted. gk values larger than bk + pk, on the other hand, must be rejected. However, for a value of gk(x(t), u(t)) within bk and bk + pk, the extent of acceptance by the DM will decrease with an increase in its value. To quantitatively demonstrate the above linguistic depiction, a fuzzy set Ck with µCk [gk(x(t), u(t))] denoting the degree of acceptability for the inequality constraint gk(x(t), u(t)) is defined as follows:
{
if gk(x(t), u(t)) < bk 1 FCk[gk(x(t),u(t)),bk,pk] if bk e gk(x(t),u(t)) e bk + pk if gk(x(t), u(t)) > bk + pk 0 (2)
where FCk denotes some monotonic decreasing function used to describe the change of degree of acceptability, as the value of gk(x, u) changes from bk to bk + pk. The general criteria for choosing function FCk are mostly subjective and indefinite and depend on individual preference. Two types of monotonic decreasing membership functions commonly used in describing the degree of acceptability are listed in the following (also shown in Figure 1 and Figure 2)
{
FCk [gk(x(t), u(t)), bk, pk] ) (bk + pk) - gk(x, u) pk
where µC represents the degree of acceptability for aggregating fuzzy inequality constraints and T, the socalled T-norm operator, denotes the operation used in the fuzzy intersection.13 Appendix A.1 provides the related definitions and operators used in the fuzzy set theory. 2.2. Fuzzification of Objective Function. From the previous statements, the flexible dynamic optimization problem is actually a parameter sensitivity problem related to the right-hand side (rhs) bounded on the path constraint. In general, the way to determine the effect of change in parameter is to solve a series of new problems once for each of the change made. This is very inefficient, however, from a computational point of view. To promote the computational efficiency, we attempted to bound the acceptable objective values of eq 1 and to grade them by the language of the fuzzy set, the fuzzification of objective values. Then, using fuzzy inference, the DM can select the resultant optimal value. To bound the objective values, we redefine eq 1 as two subproblems.
J0 ) max J[x(tf)] u(t)∈Ω0
for linear type
1 - e-R([(bk+pk)-gk(x,u)]/pk) for exponential type 1 - e-R
(3)
where R ∈ R1 is the shape-adjusting factor. Now, a new fuzzy set for the whole flexible constraints can be aggregated as
C ) C1 ∩ ... ∩ CK
(4)
(6)
where Ω0 ≡ {u(t) | h(x, u) ) x3 (t) - f(x(t), u(t)) ) 0, x(0) j (t)}, where b ≡ ) x0; g(x(t), u(t)) e b; u(t) e u(t) e u [b1, ..., bK]′, and
J1 ) max J[x(tf)]
µCk[gk(x(t), u(t))] )
(5)
u(t)∈Ω1
(7)
where Ω1 ≡ {u(t) | h(x, u) ) x3 (t) -f(x(t), u(t)) ) 0, x(0) j (t)}, where bp ) ) x0; g(x(t), u(t)) e bp; u(t) e u(t) e u [b1 + p1, ..., bK + pK]′. Herein, J0 and J1 represent the optima for the two subproblems, respectively. Assume at least one inequality constraint in eq 6 is active. From the analysis of Euler-Lagrange, we find that once the value on the rhs of gk(x, u) is increased, the optimal objective value will be correspondingly augmented. The detailed mathematical analysis is provided in Appendix A.2. By the above finding, we can conclude that the resultant optimum of eq 1 is bounded by the optimum of eq 6 and that of eq 7, that is J1 > J > J0. Moreover, owing to the variant degrees of acceptability for the value of gk(x(t), u(t)) from bk to bk + pk, the values from J0 to J1 differ as well in terms of the satisfaction of the DM. The values of J larger than J1 have the highest satisfaction. However, the values less than J0, are not acceptable. For other intermediate values between J0 and J1, the degree of satisfaction will increase with a higher value. To numerically demonstrate the change of satisfaction, the DM can define a fuzzy satisfaction, J, with the degree of satisfaction, µJ[J(x(tf)]
{
0 if J < J0 0 1 µJ[J(x(tf)] ) FJ(J, J , J ) if J0 e J e J1 1 if J > J1
(8)
Ind. Eng. Chem. Res., Vol. 39, No. 7, 2000 2307
highest degree of feasibility, to µ/D(x0|u(t)).
µ/D(x0|u*(t)) ) max µD(x0|u(t)) u(t)∈Ω ˜
) max T{µJ(x(tf)), µC(x(t))} u(t)∈Ω ˜
(13)
) max T{µJ(x(t)), µC1(x(t)), ..., µCK(x(t))} u(t)∈Ω ˜
Figure 3. Linear type monotonic increasing membership function for the objective.
Figure 4. Exponential type monotonic increasing membership function for the objective.
where FJ denote some monotonic increasing functions used to describe the change of degree of satisfaction, as the value of J changes from J0 to J1. Similar to the use of FCk, linear-type function and exponential-type function are two functions commonly used in FJ, such as shown in Figure 3 and Figure 4
{
FJ(J, J0, J1) ) J - J0 for linear type J1 - J0 (9) 0 1 0 1 - e-R(J-J /J -J ) for exponential type 1 - e-R
where R ∈ R1 is the shape-adjusting factor. Like the selection for FCk, the choice for FJ is also subjective and dependent on the DM’s preference. 2.3. Fuzzy Optimization. As objective and constraints are respectively fuzzified, the optimization problem can be linguistically described as14
Thus, eq 1 can now be interpreted as finding an optimal continuous control trajectory u*(t) over t ∈ [t0, tf] which simultaneously maintains the degree of satisfaction of objective and that of acceptability of fuzzy inequality constraints as large as possible. Comparing the above-mentioned procedure with penalty function method, we can find the analogy between them. By adding penalty terms, which consist of the constraints, to the objective function, the original problem is converted into an unconstrained optimization problem. To achieve a reliable solution, there are a sequence of problems with updating penalty factor(s) to be solved. Such procedure is tedious. In the proposed method, the DM is allowed to subjectively grade the validity for the values of path constraints between upper and lower bounds, which is similar to penalize the violation of constraints. And, after the possible objective values are bounded, the subjective preference of the DM for these values can be directly integrated into the algorithm. All these can be done by easily adjusting the type and/or shape of membership functions. Therefore, such a feature makes the interaction between the DM and the algorithm possible. 3. Solution Method 3.1. Multistage Fuzzy Control. Bellman and Zadeh10 formulated multistage fuzzy optimization problems for a time-invariant discrete system. To extend their formulation to our case, eq 1 must be rewritten as follows
max
u0,...,uP-1(t)∈Ω ˜d
J(xP)
(14)
(11)
where Ω ˜ d ) {u0, ..., uP-1| h(x, u) ) x3 (t) - f(x(t), ui) ) 0, x(0) ) x0; g(x(t), ui) e b ˜ ; u(t) e ui < u j (t) ∀t ∈[iT, (i + 1)T], i ) 0, ..., P - 1}, which represents the feasible space of discretized control input; T ) tf/P denotes time duration. With the discretized control input, [u0, u1, ..., uP-1], the degree of feasibility defined in eq 12 should be recast as
and can be measured by the degree of feasibility, µD(x0|u(t))
µD(x0|u0, ..., uP-1) ) T[µC(x0), ..., µC(xP-1), µJ(xP)] (15)
attain fuzzy objective J and satisfy fuzzy constraints C (10) The above statement implies that control input can be viewed as the aggregation (or intersection) of the fuzzy objective J and the fuzzy inequality constraint C
D)C∩J
µD(x0|u(t)) ) T{µJ(x(tf)), µC(x(t))} ) T{µJ(x(tf)), µC1(x(t)), ..., µCK, (x(t))} (12) where µD(x0|u(t)) denotes the degree of feasibility for the continuous control input u(t) starting from the initial state x0. For a collection of control actions, the most feasible control input u*(t) is the one having the
where [µC(xi) ≡ T{µC1(xi), ..., µCK(xi)} denotes the degree of acceptability for aggregating fuzzy inequality constraints at the ith stage, and µJ(xP) represents the degree of satisfaction for the fuzzy objective at final time end. The problem now becomes finding a consecutive optimal control policies u0, u1, ..., uP-1 with the highest degree of feasibility. By this, the degree of satisfaction for fuzzy objective and the degree of acceptability for
2308
Ind. Eng. Chem. Res., Vol. 39, No. 7, 2000
fuzzy constraints can be attained as much as possible, that is / µ/D(x0|u/0, ..., uP-1 ))
)
max
u0,...,uP-1∈Ω ˜d
max
u0,...,uP-1∈Ω ˜d
µD(x0|u0, ..., uP-1) T{µC(x0), ..., µC (xP-1), µJ(xP)} (16)
3.2. Fuzzy ICRS. Under the framework of eq 16, a computational technique which can determine the degree of acceptability at each grid point and the degree of satisfaction at final time end is needed. However, because the use of operator T will cause nondifferentiability, and nonconvexity is introduced by general smoothing of such an operator, the solution of the resulting problem might be impossible with a standard gradient-based NLP procedure. In recent work, ICRS/ DS12 was proposed to solve general dynamic optimization problems. By using control parametrization technique, a dynamic optimization problem is transferred into a NLP problem and then is solved by an adaptive stochastic algorithm. Because ICRS/DS does not depend on gradient information to search and the incorporation of fuzzy inference into it is very straightforward, the algorithm is used as our computation platform for solving the nonlinear dynamic optimization problems with flexible inequality constraints. These works are summarized as follows 1. Preliminary Preparations: (a) Select the type of monotonic functions for FCk, k ) 1, ... K, and FJ. In this study, a linear-type function is adopted to describe the change of decreasing membership for the path constraints and that of increasing membership for objective function. That is
{
µCk[gk(x, u)] ) if gk < bk 1 (bk + pk) - gk(x,u) if bk e gk e bk + pk (17) pk if gk > bk + pk 0
and
{
0 if J < J0 J - J0 µJ[J(x(tf)] ) 1 if J0 e J e J1 J - J0 1 if J > J1
(18)
(b) Obtain the parametrization of the control input, u(t), t ∈ [t0, tf]. The control input in ICRS/DS is parametrized using a grid of P + 1 points and approximated by some specified function. For most problems, the grid numbers from 5 to 10 are adequate. One of the following functions, such as piecewise constant, piecewise linear, B-splines and Lagrange polynomial, ..., etc., can be used in approximating the control profile. Since a variable-length piecewise linear function can cause superior solution,12 we also adopt such a function in the approximation. Furthermore, because the extension of multiple controls is straightforward, we will focus only on the case of a single input system.
Thus, the control for the iteration l during t ∈ [θil, θi+1l] is given as
ul(t) )
l ωi+1 - ωli l θi+1 - θli
(t - θli) + ωli
i ) 0, ..., P - 1
where ωi and θi denote control and time coordinates of the control parametrization grid, respectively. Notably, since θ0 ) t0 and θP ) tf are fixed, there are P - 1 time grids. The original dynamic optimization problem with single control input will result in a NLP problem of dimension 2P, where the goal is to find (P + 1) ω/i s and (P - 1) θ/i s, which optimize the objective function. Therefore, the decision vector for any iteration l in ICRS/ DS consists of
ξli )
{
θli for i ) 1, ..., P - 1 l ωi-P for i ) P, ..., 2P
(19)
Also, these decision variables are restricted by their physical bounds
ξli ∈ [ξi,ξhi] )
{
[θi,θ h i] for i ) 1, ..., P - 1 j i-P] for i ) P, ..., 2P [ωi-P,ω
(20)
where θi (or ωi) and θ h i (or ω j i) denote the upper/lower bounds on time (or control), respectively. (c) Select heuristic parameters. Basic ICRS/DS contains three heuristic parameters: K1 controls the initial step size of search (with default value 1/3); K2 is an stepsize reducing factor (with default value 1/2); ne governs the rate of convergence (with default value 4). Notably, for problem without local optima, these settings are adequate. For highly multimodal problems, K1 can be adjusted from 1.0 to 2.0.12 2. Solution algorithm. Step 0. Solve eq 6 and eq 7 by basic ICRS/DS algorithm to determine J0 and J1, and record optimal control input for each subproblem. Step 1. Read heuristic parameters, K1, K2, ne, and the upper/lower bounds for the decision variables. Set maximal loop value to Lmax, (with defaults value 5). Set l ) 0 and L ) 1. Put J0, J1, bk, and bk + pk into selected membership functions. Use the decision vector obtained from the previous subproblems as the initial guess. Step 2. Evaluate the degree of feasibility, µlD for the initial control input. The procedures are stated as follows. Step 2a. Use the initial decision variable presented above, integrating the differential equation(s). Record the state’s value at each time grid and the value of objective function at tf. To avoid the violation of path constraints occurring between grid points, additional data recording between grids can be used to monitor the state path. Step 2b. Calculate the degree of acceptability of each constraint at each time grid, µC(xi), i ) 0, ..., P- 1, and the degree of satisfaction of objective function, µJ[xP]. Step 2c. Aggregate µC(xi), i ) 0, ..., P - 1 and µJ[xP] by a selected T-norm to determine µ0D. Step 2d. As µ0D is not equal to zero, the decision vector is a candidate, so go to Step 3. Otherwise, return to Step 2a.
Ind. Eng. Chem. Res., Vol. 39, No. 7, 2000 2309
Step 3. Generate a new decision vector ξl+1 using i the old one, ξli
ξl+1 ) ξli + γiK1σli, i ) 1, ..., 2P i
(21)
where σli, which is defined as {min {(ξhi - ξi), (ξi - ξi}}, denotes the standard deviation and γi is a pseudorandom number obtained from normal distribution of the zero mean and the unity standard deviation. Go to Step 4 if ξl+1 ∈ (ξi,ξi). Otherwise redo Step 3. i Step 4. Evaluate the new degree of feasibility µl+1 D with the new decision vector, following the same procedure as in Steps 2a-d. Go to Step 6 if the new decision is feasible and the degree of feasibility is improved l (µl+1 D > µD). Otherwise, go to Step 5. Step 5. Increase the failure counter, F r F + 1. Return to Step 3 if F e (ne × 2P). If F > (ne × 2P), then add the loop counter, L r L + 1. If L e Lmax, then F r 1, σli r K2σli, i ) 1, ..., 2P, and return to Step 3. Otherwise, save the data and stop. Step 6. Check the convergence criterion by degree of feasibility according to the specified tolerance e: l |µl+1 D - µD |
|µlD|
ξi )
ξl+1 i
(
dx1 ) (1 - x1)k1 - x1k2 dt dx2 ) 300[(1 - x1)k1 - x1k2] - u(x2 - 290) (23) dt where
k1 ) 1.7536 × 105 exp
( (
k2 ) 2.4885 × 1010 exp
) )
-1.1374 × 104 1.9872x2
-2.2748 × 104 1.9872x2
(24)
and with the initial condition
x(0) ) [0 380]T
(25)
where x1 and x2 denote the normalized concentration of the desired product and the temperature, respectively. The control variable is normalized coolant flow rate and is bounded by
470 0
To determine the band of objective value, the original problem is divided into two subproblems, one is restricted by
x2(t) e 450
(29)
and the other is restricted by
x2(t) e 470
(30)
2310
Ind. Eng. Chem. Res., Vol. 39, No. 7, 2000
Figure 5. Concentration profile in example 1.
Figure 6. Temperature trajectory in example 1.
To avoid the violation of path constraint, eqs 29 and 30, in determining J0 and J1, we introduce an additional state such as
dx3 ) dt x2 - 450 (or x2 - 470) if x2 > 450 (or x2 > 470) if x2 e 450 (or x2 e 470) 0 (31)
{
with initial conditions x3(0) ) 0. The objective function in the subproblems is redefined by the augmented performance index to be minimized as
-x1(tf) + ηx3(tf) where η is a positive penalty factor. As Luus reported,15 η is not very important in the case, and is set as 1000 in our calculation. The optimal values for J0 and J1 are 0.66904 and 0.67984, respectively. By using both values, we also employ the linear-type function to depict the change of degree of feasibility for the fuzzy objective:
{
µJ[J(x1)] ) 0 if J < 0.66904 J - 0.66904 if 0.66904 e J e 0.67984 (32) 0.67984-0.66904 1 if J > 0.67984 By using P ) 10 and a Zadeh-min operator, the resulting concentration for the desired product x2(tf) is 0.67698 and the degree of feasibility for the control input (coolant flow rate) is 0.7350, such as shown in Figure 5. The resulting optimum is between those considering conservative constraints and loose constraints, respectively. The temperature trajectory and coolant flow rate are also exhibited in Figure 6 and Figure 7. Figure 6 shows that the temperature trajectory is kept in the allowable range. However, during the period from t ) 1.87 min to t ) 2.91 min, the temperature is allowed to move from 451.93 to 455.30 K for increasing the concentration of the desired product. It is worthwhile to note the dip in the control input around t ) 2.63 min and then the rise to the second peak. Luus reported a similar phenomenon15 and concluded that it is due to the state constraint. We conjecture that the dip cause the temperature not extremely exceeding 450 K. Therefore, the degree of acceptability can be maintained at an acceptable level. For the next rising peak, it keeps the temperature for increasing the concentration of the desire product.
Figure 7. Resulting control input in example 1.
Example 2: A Mathematical System with Nonlinear Inequality Constraint.15 The system is governed by three differential equations as follows
dx1 ) x2 dt dx2 ) -x2 + u dt
(33)
dx3 ) x12 + x22 + 0.005u2 dt with the initial condition
x(0) ) [0 - 1 0]T
(34)
the control input is bounded by
-20 e u e 20
(35)
The flexible constraint is
g(x) ) x2 + 0.5 - 8(t - 0.5)2 e0˜
(36)
The objective function to be minimized is
x3(tf)
(37)
where tf ) 1. Now, if the limiting value and maximal tolerance for the constraint g(x) are set as 0 and 0.15, and the linear
Ind. Eng. Chem. Res., Vol. 39, No. 7, 2000 2311
Figure 8. Trajectory for x3 in example 2.
Figure 9. Resulting control input in example 2.
function is adopted to depict the degree of acceptability, we obtain
Table 1. Data for Constraint Trajectory in Example 2
{
if g(x) < 0 1 0.15 - g(x) if 0 e g(x) e 0.15 µC[g(x)] ) 0.15 - 0 if g(x) > 0.15 0
times grid
g(x)
times grid
g(x)
0.00 0.20 0.22 0.37
-2.50000 -0.23295 -0.16851 0.05199
0.60 0.65 0.91 1.00
0.05199 0.03821 -0.84322 -1.43279
(38)
As in the previous statement, we introduce an additional state to detect the violation in computing J0 and J1
{
dx4 g(x) if g(x) > 0 (or g(x) > 0.15) ) if g(x) e 0 (or g(x) e 0.15) 0 dt
(39)
with initial conditions x4(0) ) 0. The objective function in the subproblems is redefined by the augmented performance index to be maximized as
-x3(tf) - ηx4(tf) where the positive penalty factor η is set as 1000 in this case. By using ICRS/DS, the optimal values of the two subproblems restricted by
g(x) ) x2 + 0.5 - 8(t - 0.5)2 e 0
(40)
and
g(x) ) x2 + 0.5 - 8(t - 0.5)2 e 0.15
(41)
respectively, are J0 ) -0.17659 and J1 ) -0.12596. To comply with eq 18, the linear function characterizes the degree of satisfaction µJ[J] as follows
{
µJ [J] ) 0 if J < -0.17659 J - (-0.17659) if -0.17659 e J e -0.12596 -0.12596 - (-0.17659) 1 if J > -0.12596
(42)
In this example, P is set as 7 and the Zadeh-min is used as the operator of fuzzy intersection. The resulting optimal value for x3 is 0.14351, and the degree of satisfaction for control input is 0.6534. The trajectory of x3 and the control profile for this example are shown in Figure 8 and Figure 9. As shown in Table 1 and Figure 10, the value of the inequality state constraint g(x) is smaller than 0 during t ∈ [0, 0.37]. To improve the degree of feasibility, it is allowable to slightly violate the constraint during t ∈ [0.37, 0.65].
Figure 10. Trajectory for g(x) in example 2.
Example 3: Fed-Batch Fermentor Problem.4,15-18 The system is governed by four differential equations:
( )
dx1 x1 ) h1x1 u dt 500x4
( )
x2 dx2 ) h2x1 - 0.01x2 u dt 500x4 h1x1 h2x1 0.029x3x1 dx3 + )dt 0.47 1.2 0.0001 + x3
(
1-
dx4 u ) dt 500 where
h1 )
0.11x3 0.006x1 + x3
and
h2 )
0.0055x3 0.0001 + x3(1 + 10x3)
)
x3 u (43) 500 x4
2312
Ind. Eng. Chem. Res., Vol. 39, No. 7, 2000
with the initial conditions
x(0) ) [1.5 0 0 7]T
(44)
the constraint on feed rate (control input) is
0 e u e 50
(45)
and the flexible path constraints are
0 e x1 e 40 ˜
(46)
˜ 0 e x3 e 25
(47)
0 e x4 e 10 ˜
(48)
The limiting values for x1, x3, and x4 are 40, 25, and 10, respectively. The maximal tolerance for the same states are 5, 5, and 2. The objective to be maximized is the amount of penicillin at the final time, tf ) 126, given by
x2(tf)x4(tf)
(49)
As the linear functions are used, the mathematical description of the degree of acceptability can be defined as
{ { {
if x1 < 40 1 45 - x1 if 40 e x1 e 45 µC[x1] ) 45-40 if x1 > 45 0
(50)
if x3 < 25 1 30 - x3 if 25 e x3 e 30 µC[x3] ) 30-25 if x3 > 30 0
(51)
Acknowledgment The authors would like to thank the anonymous reviewers for their constructive criticism and suggestions. In particular, the authors would give our appreciation to Dr. Julio R. Banga of Instituto Investigacions Marinas (CSIC) in Spain for his sincere helps on this research. Nomenclature
if x4 < 10
1 12 - x4 if 10 e x4 e 12 µC[x4] ) 12-10 if x4 > 12 0
such tedious procedure, we first determine the upper and lower bounds on the resulting optimum by solving the problems subjected to the loosest constraints and that subjected to tightest constraints, respectively. From mathematical analysis, we prove that the resulting optimal value is within the lower/upper optimal bounds. Because of the variant degree of acceptability for the bounding value of inequality constraints, the objective values can be considered as ones having different degrees of satisfaction. By using the bounding values on inequality constraints and objective value, the depiction of the degree of acceptability for each flexible constraint and that of the degree of satisfaction can be built. Such a procedure is the so-called fuzzification of objective function and inequality constraints. Also, the original problem has been formulated as a fuzzy nonlinear dynamic optimization problem. To seek the resultant optimal value, the control input in the model is discretized. Then, the fuzzy inference is incorporated into ICRS/DS to solve the problem. Three examples are used to exhibit the facility in this algorithm. Results of these examples demonstrate that the proposed modification is very simple and efficient on solving the dynamic optimization problems with flexible inequality constraints.
(52)
By using basic ICRS/DS, the bounds on objective value can be determined as J0 ) 86.93 and J1 ) 137.87. Therefore, the change of degree of satisfaction can be characterized as follows
{
µJ[J(x2x4)] ) 0 if J < 86.93 x2(tf)x4(tf) - 86.93 if 86.93 e J e 137.87 (53) 137.87-86.93 1 if J > 137.87 The resultant optimal for the amount of penicillin at final time end is 123.49 with 0.73 of degree of feasibility for the feed rate. 5. Conclusions This work proposes an algorithm to solve the dynamic optimization problem with flexible inequality constraints. Such kind of problem is actually a sensitivity problem. To solve this problem, the DM needs to solve a series of optimization problems with relaxing the bounding value of inequality constraints. To overcome
FJ ) membership function for fuzzified objective function FCk ) membership function for flexible path constraints T ) T-norm operator for fuzzy intersection g ) path inequality constraints of the state h ) dynamic portion of the system u(t) ) m × 1 control vector x(t) ) n × 1 state vector bk ) reasonable imiting value for path constrain, gk F ) failure counter in ICRS/DS algorithm J ) Objective function K1 ) heuristic parameter in ICRS/DS algorithm K2 ) heuristic parameter in ICRS/DS algorithm L ) loop counter in ICRS/DS algorithm ne ) heuristic parameter in ICRS/DS algorithm P ) numbers of time stage pk ) acceptable maximal tolerance for path constrain, gk Greek Letters R ) the shape-adjusting factor µD ) the degree of feasibility for the continuous control policy, u(t) µCk ) the degree of acceptability for the inequality constraint gk (x(t), u(t)) µJ ) the degree of satisfaction for fuzzy satisfaction, J ω ) parameter grid for control input θ ) parameter grid for time Ω ˜ ) feasible control input space ξ ) decision variables
Ind. Eng. Chem. Res., Vol. 39, No. 7, 2000 2313
constraints are also zero. Thus, eq 59 in its stationary point, (ξ*, λ*, µ*), must satisfy the following condition:
Superscripts * ) value of variable at optimal point
A Appendix
∂L
A.1. Fuzzy sets and Fuzzy Intersection. Definition 1 (Fuzzy Set). A fuzzy set A in a universe of discourse X ) {x}, is defined as a set of pairs.
∂ξl
A ) {(x, µA(x))|x ∈ X}
(54)
where µA: X f [0, 1] is a membership function of x in A. It is worthwhile to note that the membership function µA(x) represents the extent that x belongs to A and is normally limited to values from 0 to 1. The higher the value of µA(x), the higher the extent that x belongs to A. Definition 2 (Fuzzy Intersection). Let A and B be two fuzzy sets with membership functions µA(x) and µB(x), respectively. The membership functions of the intersection D ) A ∩ B, is defined by
µD(x) ) T{µA(x), µB(x)}, x ∈ X
{
min{µA(x), µB(x)} Zadeh min µA(x)‚µB(x) algebraic product
max J(ξ) (57)
g(ξ) e b where ξ ∈ R2P denotes the decision vector, and J: R2P f R, h: R2P f Rn, g: R2P f RK. To proceed to the analyses, we need the following definitions. Definition 3. An inequality gk e bk is said to be active at ξ if gk(ξ) ) bk. It is inactive at ξ if gk(ξ) < bk. Definition 4. Suppose that ξ* satisfies h(ξ*) ) 0, g(ξ*) e b, and let A be the index set of active inequality constraints, that is
A(ξ*) } {k: gk(ξ*) ) bk, 1 e k e K}
(58)
Definition 5. For some function f of ξ, a point ξ* at which f ′(ξ) ) 0 is called a stationary point. By convention, equality constraints are always considered to be active. Let us introduce a Lagrange function to solve eq 57
L(ξ, λ, µ,) ) J(ξ) +
∂ξl
+
∂hi
λi ∑ ∂ξ ∀i
λihi(ξ) + ∑µk[gk(ξ) - bk + yk2] ∑ ∀i ∀k
(59)
where Y ) [y1, ..., yK]′, the vector of the slack variable, is used to transformed the inequality constraints in eq 57 into equalities; λ′is are the Lagrange multipliers associated with h′is, and µ′ks, are the Karush-KuhnTucker (KKT) multipliers associated with g′ks, respectively.20 From mathematical programming theory,21 the corresponding KKT multipliers for inactive inequality constraints are zero, and the slack variables for active
∂gk
µk ∑ ∂ξ ∀k
+
l
)0
l ) 1, ..., 2P
l
(60)
Notably, the perturbation of ξ from the optimal point, ξ*, will cause J(ξ) < J(ξ*) for the maximization problem, therefore
0 > J(ξ) - J(ξ*) ≡ dJ(ξ*) )
[
[
)-
∂J
dξl ∑ ∀l ∂ξ l
∂hi
∑ ∑ λi ∀l ∀i ∂ξ
)-
∂hi
λi∑ dξl ∑ ∀i ∀l ∂ξ l
+
µk ∑ ∂ξ ∀k
] [
)-
l
]
∂gk
-
ξ ) ξ*
dξl
l ξ ) ξ*
∂gk
]
µk∑ dξl ∑ ∀k ∀l ∂ξ l
ξ ) ξ*
λidhi(ξ*) - ∑µkdgk(ξ*) ∑ ∀i ∀k
(61)
In this equation, dhi(ξ*) ) 0, ∀i, because the equality constraints must be satisfied, and µk ) 0, ∀k ∉A. Thus we can conclude that
0 > dJ(ξ*) ) -
(56)
A.2. The Effect of Relaxing Inequality Constraints on Optimal Value. Consider the following constrained maximization problem
s.t. h(ξ) ) 0
∂J
(55)
where T denotes a general class of intersection operators for the fuzzy set. Two frequently used T-norm operators19 are
µD(x) )
)
∑ µkdgk(ξ*)
(62)
k∈A
Furthermore, since dg(ξ*) e 0, ∀k ∈ A, for any feasible perturbation. Thus
µk < 0 ∀k ∈ A
(63)
Next, let us consider the sensitivity of relaxation of the bounding value for the active inequality constraints at ξ*. It is clear that
dgk ) dbk ∀k ∈ A
(64)
Multiplying the above equation by the associated KKT multipliers, µk, and summing them up, we obtain
∑ µkdbk ) k∈A ∑ µkdgk ) -dJ > 0
(65)
k∈A
Thus
dJ ) -µk > 0 dbk
∀k ∈ A
(66)
It implies that once the value on the rhs of an active inequality constraint gk(ξ*) e bk, k ∈ A is relaxed, the optimal objective value will be correspondingly augmented. Literature Cited (1) Denbigh, K. G. Optimal temperature sequences in reactors. Chem. Eng. Sci. 1958, 8, 125. (2) Aris, R. On Denbigh’s optimum temperature sequence. Chem. Eng. Sci. 1960, 12, 56. (3) Biegler, L. T. Solution of dynamic optimization problems by successive quadratic programming and orthogonal collocation. Comput. Chem. Eng. 1984, 8, 243. (4) Cuthrell, J. J.; Biegler, L. T. Simultaneous optimization and solution methods for batch reactor profiles. Comput. Chem. Eng. 1989, 13, 49.
2314
Ind. Eng. Chem. Res., Vol. 39, No. 7, 2000
(5) Goh, C. K.; Teo, K. L. Control parametrization: a unified approach to optimal problems with general constraints. Automatica 1988, 24, 3. (6) Vassiladis, V. S.; Sargent, R. W. H.; Pantelide, C. C. Solution of a class of multistage dynamic optimization problems. 1. Problems without path constraints. Ind. Eng. Chem. Res. 1994, 33, 2111. (7) Vassiladis, V. S.; Sargent, R. W. H.; Pantelide, C. C. Solution of a class of multistage dynamic optimization problems. 2. Problems with path constraints. Ind. Eng. Chem. Res. 1994, 33, 2123. (8) Luus, R. Optimal control of batch reactor by iterative dynamic programming. J. Process. Control 1994, 4, 218. (9) Zadeh, L. A. Fuzzy sets. Inf. Control 1965, 8, 338. (10) Bellman, R. E.; Zadeh, L. A. Decision making in a fuzzy environment. Manage. Sci. 1970, 17, 141. (11) Sakawa, M. Fuzzy sets and interative multiobjective optimization; Plenum Press: New York and London, 1993. (12) Carrasco, E. F.; Banga, J. R. Dynamic optimization of batch reactors using adaptive stochastic algorithms. Ind. Eng. Chem. Res. 1997, 36, 2252. (13) Klir, G. J.; Yuan, B. Fuzzy sets andfuzzy logics-theory and application; Prentice Hall: London, 1995. (14) Kacprzyk, J. Multistage fuzzy control; John Wiley and Sons: New York, 1997.
(15) Wichaya, M.; Luus, R. Optimal control of inequality state constrained systems. Ind. Eng. Chem. Res. 1997, 36, 1686. (16) Luus, R. Optimization of Fed-Batch Fermentors by Iterative Dynamic Programming. Biotechnol. Bioeng. 1993, 41, 599. (17) Luus, R. Piecewise Linear Continuoue Optimal Control by Iterative Dynamic Programming. Ind. Eng. Chem. Res. 1993, 32, 859. (18) Dadebo, S. A.; McAuley, K. B. Dynamic Optimization of Constrained Chemical Engineering Problems Using Dynamic Programming. Comput. Chem. Eng. 1995, 19, 513. (19) Zimmermann, H. J. Fuzzy set theory and its applications; Kluwer Academic Publisher: London, 1991. (20) Chong, E. K. P.; Zak, S. H. An Introduction to Optimization; John Wiley and Sons: New York, 1996. (21) Rao, S. S. Engineering Optimization; John Wiley and Sons: New York, 1996.
Received for review August 23, 1999 Revised manuscript received April 24, 2000 Accepted April 25, 2000 IE990638N
