Ind. Eng. Chem. Res. 1994,33, 2123-2133
2123
Solution of a Class of Multistage Dynamic Optimization Problems. 2. Problems with Path Constraints V. S. Vassiliadis, R. W. H. Sargent, and C. C. Pantelides' Centre for Process S y s t e m Engineering, Imperial College of Science, Technology and Medicine, London SW7 2BY, United Kingdom
This paper considers the treatment of general equality and inequality path constraints in the context of the control vector parametrization approach t o the optimization of dynamic systems described by mixed sets of differential and algebraic equations (DAEs) of index not exceeding 1. Equality path constraints are handled by incorporating as many of them as possible within the DAE system itself without increasing its index. This allows a subset of the control variables to be determined from the solution of the augmented DAE system. The issues involved in establishing an appropriate partitioning of the control variable vector are examined. Inequality path constraints are handled through the combined application of the discretization of these constraints at a finite number of points, and forcing an integral measure of their violation to zero. Numerical experiments demonstrating the advantages of this hybrid technique over its individual components are presented. 1. Introduction
The first part of this paper (Vassiliadis et al., 1994) presented an algorithm for the optimization of systems, the transient behavior of which is described by multistage sets of differential and algebraic equations (DAEs) of the form
f&,
t ) = 0 t E B k - 1 , tk), k = 1, ..., NS- 1; t E ItNs-1, t J , k NS (1) where NS is the number of stages and t is time. The sets f , Y , 4 u,
of differential and algebraic variables are denoted by x ( t ) E X C Rn and y ( t ) E Y C Rm,respectively, with f being the time derivatives of x ( t ) . u ( t ) E U C R* and u E V C Rq are sets of time-varying control variables and timeinvariant parameters, respectively, and f k : X X Rn X Y X U X V X [tk-l, tk) R"+m. A variety of constraints enforced at a finite set of distinct time points (corresponding to the stage boundaries) during the time horizon of interest were considered. These included stage end-point constraints, as well as general initial and junction conditions. In this second part of the paper, we turn our attention to the handling of equality and inequality constraints that need to be enforced over finite time intervals. We assume that the latter coincide entirely with one or more stages, defining additional artificial stages for this purpose, if necessary. We are therefore interested in stage path constraints of the form
-
c F ( x ( t ) , fW, Y W ,
w ,t ) = 0 0,
k = 1,...,NS - 1; where
cp: X
cF(x(t),
X
Rn
X
Y
X
t
E It,-,, [tk-l, t k )
We begin this section by considering established approaches to dealing with equality path constraints and pointing out some numerical problems that may arise from their application. We then proceed to consider a uniform manner of treating all path equalities (including the original DAE system (1))in the optimization problem, analyze the complications that may arise, and propose a technique that circumvents them. A simple example illustrating the application of this technique is presented. 2.1. Established Techniques for Handling Equality Path Constraints. A widely used method involves the introduction of integral penalty terms in the objective function to be minimized (Bryson and Ho, 1975). In the context of the multistage problem considered here, the penalty term would be of the form NP
(4)
-
RXh,
and
w ,Y W , w , u, t ) 5 0 k = 1,...,NS - 1;
2. Equality Path Constraints
tk),
t E [tNgl, tr], k = NS (2)
UX VX
Section 2 of this paper considers equality path constraints (2), attempting to establish a unified manner of treating all path equalities (including both (1)and (2)) in the problem, while paying special attention to issues relating to the index of the resulting DAE systems. Section 3 considers inequality path constraints (3), proposing a hybrid technique which combines elements of earlier approaches for handling such constraints. Numerical results demonstrating the effectiveness of this technique are presented in section 4.
t E ltk-1, tk), t E [t,,,, t f ] , k = NS (3)
-
where c$ X X Rn X Y X U X V X [&-I, t k ) Rxk. In particular, we are concerned with the effective and efficient handling of such path constraints in the context of control vector parametrization algorithms.
* Author to whom correspondence should be addressed. Electronic mail address:
[email protected].
where K is a large positive constant. An alternative approach that seeks to avoid the numerical difficulties which may be caused by the use of the penalty term involves the conversion of the path constraint to an equivalent end-point constraint (Sargent and Sullivan, 1977). In the context of the multistage problem, this would involve the introduction of a new differential variable Z ( t ) defined through the stage differential equations
3 = Ilcpcxct,, f ( t ) ,Y ( t ) ,u ( 0 ,u, t)l1,2, t E [ t k - , , t k ) , k = 1,..., NS - 1; t E [ t N g 1 , tfl, k = NS (5a) and the initial and junction conditions 1994 American Chemical Society
2124
Ind. Eng. Chem. Res., Vol. 33, No. 9, 1994
R(tk)
= Z(t,),
k = 1, ...,NS - 1
(5c)
The path constraints (2) are then equivalent to the single end-point constraint Z(tf) = 0
(6)
In some cases, if there are significant differences in the magnitude of the different elements of the path constraint vector cp(.), and/or of the same element in different stages, better numerical behavior could be ensured by disaggregating 2 into measures of the square violation of individual path constraints in individual stages, and imposing a separate end-point constraint on each one of them. A common characteristic of all the techniques discussed so far is that the use of square norms implies that the penalty terms or end-point constraints introduced have zero gradients with respect to the optimization parameters at the solution. This, in turn, may result in a reduced convergence rate near the solution, and as noted by Goh and Teo (19881, the success of such techniques depends very strongly on the line-search merit function used by the optimization algorithm. Interestingly, the approach described later in this paper for dealing with inequality path constraints can also be used to alleviate these problems (see concluding remarks in section 5). 2.2. A Unified Treatment of Path Equalities. An alternative way of handling equality path constraints (2) is based on the observation that there is, in principle, no mathematical difference between these and the system equations (1) even if, in practice, their origins may be different. For instance, eqs 1 would normally describe the underlying physics of the system under consideration, while constraints 2 could express desired properties to be attained under the influence of the controls u(t). We can therefore view the combined set of equations (1)and (2) as a DAE system to be solved simultaneously. The addition of hk new equations to the system removes an equal number of degrees of freedom from the problem, effectively converting hk elements of the control vector u ( t ) into algebraic variables, which, like y ( t ) , will be determined by the solution of the combined DAE system. The remaining elements of u(t) continue to be treated as control variables, the time variation of which will be established by the optimization. We denote the two types of control variables as &(t) and &(t) respectively. It should be noted that this partitioning of the u ( t ) may differ from stage to stage. However, the above simple approach suffers from a number of complications: (a) The combined DAE system may be of index exceeding unity (see Brenan et al. (1989) for a discussion of DAE index). In fact, this will always be the case if one or more of the constraints (2) is a function of the differential variables x ( t )only. Unfortunately, the numerical solution of most high index problems is much more difficult than that of problems of index 0 or 1. (b) Unless hk = T , there mayhe more than one possible partitioning of u into and i i k . Different choices may lead to a different index for the combined DAE system. (c) General solvability conditions for high index problems have not yet been established. Therefore, it may be difficult to guarantee that, for a given partitioning of the u ( t )vector, solutions ( x ( t ) , Y(t)Liik(t)) exist for d l possible variations of the elements of
(d) In many applications, the space of practically implementable control functions may be restricted (e.g., to be piecewise constant or linear) for certain elements of u(t). As the solution of DAE systems do not, in general, belong to such function spaces, the cqrresponding control variables must be delegated to the i i k partition of u. In view of the above complications, it may be inevitable in some cases that some of the constraints (2) must be handled by the optimization using one of the already established techniques of section 2.1. Nevertheless, we seek to establish, for every stage k, a maximal subset of the constraints that can be treated as part of the DAE system without giving rise to the above problems. An appropriate partitioning of u ( t ) (7)
also needs to be established. Let EL@) be the subset of the control variables u which can-take any functional form, and denote its cardinality by T . Clearly, the partition will be a subset of EL, while i i k will be a superset of its complement. We now form the matrix
and denote its rank by r k . Through the application of a set of row and column pivoting operatings, we effect the transformation
where Lk is a lower triangular matrix with unit diagonal, k is a rk X rk upper triangular matrix. and u Assuming that the DAE system (1)is of index 1, with the matrix
being nonsingular, we can guarantee that r,Ln+m (11) and it is then possible to select the permutation matrices Pk and Qk to be of the special structure
where Pk' and Qk' are ( n + m ) X (n m! permutation matrices, and Pk" is a hk X hk and Qk" a a X T permutation matrix. This has the effect of ensuring that row pivoting takes place separately among the rows corresponding to f k ( ' ) and cp(.),and also that column pivoting takes place separately among the columns corresponding to {x,y ) on one hand and to EL on the other. Thus the rows and columns in the null submatrix 0 in (8)always correspond to subsets of cp(.) and EL, respectively. Assuming that Pk and Qk are such that it is everywhere possible to determine a matrix Lk that will transform Ak in the manner indicated by (91, then the DAE system formed from (1) and the first ( r k - n - m) elements of Pk"cp(') has a global index of 1with respect to a variable vector comprising x , y , and the first (rk - n - m ) elements of Qk"G. Therefore, this system can be solved to establish
Ind. Eng. Chem. Res., Vol. 33, No. 9, 1994 2125 a unique time variation of these variables for any given variation of the rest of the control variables u. In the context of the control vector_parametrization algorithm considered here, the last (a - rk + n + n) elements of &''a together with any elements of u not in Ei will be manipulated directly by the optimization procedure. Also the last ( x k - rk + n + m) elements of cF(.)will have to be treated using one of the techniques of section 2.1. 2.3. A Simple Example. We consider the following linear DAE system comprising three equations, and two differential, one algebraic, and three control variables:
kl
-kl + x1 + x 2 + 2y1 + u1= 0
(13a)
+ 3X1 + 2X2- y1+ 3U1+ 5U2 = 0
(13~)
We now impose the following two equality path constraints on the solution of the above system:
+ x2-
x1
1= 0
(14a)
2x1 + x 2 - 2u1 - u3 = 0
A=
(: 0'
2
1
0
J'
1 !
0
0
0
0
2
0
0
- 2 0
0
")
2 0
-1
0 -1
0 1
3.1. Established Techniques for Handling Inequality Path Constraints. Most of the established approaches for dealing with inequality path constraints are similar to the techniques for handling equality path constraints presented in section 2.1. Thus, they rely on defining a measure of the constraint violation over the entire horizon, and then penalizing it in the objective function, or forcing it directly to zero through an endpoint constraint. One possible measure of the violation of a scalar constraint (3) is
where the exponent f is typically either 1or 2. Although f = 1is sufficient, using f = 2 has the advantage that the function (max(0,w))r has first-order continuity at w = 0. Vectors of path constraints can be handled by defining either a separate violation for each element of the vector, or a combined violation for the entire vector. A possible measure for the latter case is
(15)
-1
1 0
0
0 0 using the permutation matrices
1 l 70
3. Inequality Path Constraints
0
which can be transformed to - 1 0
latter will be treated as an algebraic variable to be determined by the solution of the augmented DAE system. On the other hand, the optimization procedure will manipulate directly the remaining control variables, ul(t) and u3(t),with constraint 14a being handled using one of the techniques of section 2.1.
(14b)
Furthermore, we assume that no restrictions exist on the functional form that any of the three control variables can take. Thus, in this case, Ei 1 u. To apply the procedure of section 2.2, we form the Jacobian matrix of the combined system (13) and (14) with respect to {kl,k2, y1, u1,u2,u3):
- 1 0
P and Q suggests that this augmentation will involve path constraint 14b and the second control variable u&). The
0 0
where [crlj denotes the j t h element of cp. It is well recognized that a major common disadvantage of all the above techniques arises from the use of the max operator in the violation measures (18) and (19). More specifically, when the constraint is inactive, both the measures and their gradients with respect to all optimization parameters are zero. Thus no useful information can be conveyed to the optimizer regarding the proximity of the current point to the boundary of the feasible region. Overall, this may result in inefficient behavior involving excessive oscillations between feasible and infeasible choices of the optimization parameters in successive optimization steps or during the line-search procedure. An alternative way of handling inequality path constraints (Jacobson and Lele, 1969) involves the introduction of squared slack variables s which convert the inequalities to equalities. Thus a scalar constraint (3) is modified to c F ( x ( t ) , k ( t ) , y(t), u ( t ) , u, t )
0 and the lower triangular matrix
L=
0
)
1
1 (17c) -4 1 0 0 0 1 Clearly the Rank of A is 4, and hence the original DAE system (13)can be augmented by one more equation and one more variable. The form of the permutation matrices
(! -2
0 0
+ s2 = 0
(20)
In principle, s ( t ) could be treated as a control variable to be determined by the optimization while (20) could be appended to the original DAE system (11, effectively determining the time variation of one of the control variables u ( t ) . In practice, this technique suffers from one major disadvantage. In the very common situation of (3) involving only differential variables x ( t ) , the augmentation of (1) by (20) results in a high-index DAE system which must then be reduced to an index 1problem by a procedure involving repeated differentiations of (20) (and, possibly, some of the equations in (l)), and algebraic manipulations.
2126 Ind. Eng. Chem. Res., Vol. 33, No. 9, 1994
If the index of the system is K , then ( K - 1)differentiations are required, introducing time derivatives $, s, ..., s+l). The latter is then treated as the control variable to be manipulated by the optimization, while the lower order derivatives and s ( t ) itself are treated as distinct variables, related through equations of the form
Zj(0) = 0, j = 1,
e..,
(24b)
Xk
junction conditions Zj(tk)
= Rj(ti), j = 1,
..e,
X,, k = 1, NS - 1 ( 2 4 ~ ) e..,
and end-point constraints Zj(tf) Ie, j = 1, ...,i;,
3.2. A Hybrid Approach to Inequality Path Constraints. Integral expressions of the type shown in eq 18 provide a single measure of the violation of an inequality path constraint over the entire time horizon of interest, permitting, for instance, the replacement of each such constraint by a single end-point constraint (cf. eq 6). The techniques discussed in section 3.1 are particularly convenient for handling inequality path constraints in the context of the control parametrization approach to the solution of dynamic optimization problems. On the other hand, the complete discretization approach, of the type advocated, for instance, by Cuthrell and Biegler (1987) and Logsdon and Biegler (19891, allows a much simpler treatment of such constraints. Thus, each inequality path constraint is enforced as a point inequality constraint at all the points considered by the discretization, in exactly the same manner as the system equations (1). It should be noted that, even within control vector parametrization algorithms, it is in principle possible to adopt an approach to inequality path constraints similar to that of complete discretization algorithms. Thus, the path constraints (3) could be replaced by
i = 1,2, ..., k = 1, ...,NS (22) where Tik E [ t k - l , tk] are a given set of points in stage k . The point constraints (22) and their gradients with respect to the optimization parameters can readily be evaluated through the integration of (1)and the sensitivity equations (see Vassiliadis et al., 1994). The disadvantage of the approach described above is that a rather large number of points Tik might be necessary to ensure that the path constraints (22) are actually not violated between consecutive Tik. We therefore adopt a hybrid approach that combines the advantages of the compactness of representation of (3) afforded by the use of the integral violation measures together with the increased information provided to the optimizer by the point constraints (22). More specifically, we replace (3) by a set of point constraints at the stage boundaries only c F ( x ( T i k ) , *(Ti,),
Y ( T i k ) , u ( T i k ) , u, Ti,)
CF(x(tk), i ( t k ) , Y(tk), U(tk),
50
u, t k ) 5 0
k = 0, NS (23a)
in addition to an end-point constraint forcing the integral violation of each element of the constraint vector (3) to zero. This is achieved by introducing a new differential variable R for each inequality path constraint, defined through
with initial conditions
(244 The use of a small nonnegative constant e on the righthand side of (24d) to relax the end-point constraint (6) has been suggested by Walsh (1993), and has been found to lead to significant performance improvements in some cases. In any case, it is assumed that the number of inequality path constraints remains the same from stage to stage although their functional form may change. The slight complication arising if this is not the case can readily be dealt with at the implementation level. In summary, each inequality path constraint requires the introduction of an extra simple ordinary differential equation (ODE) (24a), and is converted to- 2NS point constraints (23) and one end-point constraint (24b). The motivation for this is that (23) will provide some guiclance to the optimization during the search for the solution, while (24b) will ensure that the final solution satisfies the constraint everywhere within the DAE integration tolerance and the constant e. 4. Numerical Experiments with Inequality Path Constraints The hybrid technique presented in section 3.2 has been implemented within the DAEOPT dynamic optimization code (Vassiliadis, 1993a). In this section, we demonstrate its effectiveness through results obtained by applying DAEOPT to three dynamic optimization problems involving inequality path constraints. The first example is derived from a well-known optimal control test problem, while the other two are chemical engineering problems concerning the optimal operation and design of chemical reactors. 4.1. Constrained van der Pol Oscillator. We use the form of this problem employed by Gritsis (1990).This can be written as min x3(5)
(25)
subject to i1 = (1- x?)xl
2, = x12
- x2 + u
+ x; + u2
@a)
(26~)
with the initial conditions x ( 0 ) = [O,1,OIT,and the control variable u ( t ) restricted between lower and upper bounds of -0.3 and 1.0, respectively. To the above, we add the inequality path constraint x l ( t ) L - 0.4, t E [O,51 (27) For the purposes of the solution using DAEOPT, the time domain [0, 51 was divided into 10 control intervals, the lengths of which were bounded between 0.05 and 4.5. A piecewise linear approximation to the control variable u ( t ) was employed. The optimization was started with a uniform control profile u ( t ) = 0.7, t E [O, 51. The above problem was solved both without and with constraint 27. For the unconstrained problem, an optimal value of 2.868 75 was obtained for the objective function,
Ind. Eng. Chem. Res., Vol. 33, No. 9,1994 2127 Van der Pol Oscillator
:I
0.90 1.00
0.70
am -
am -
t I
050 0.40 0.60
0.30
-
0.20
-
0.10
-
4.00
-
4.10
-
-
Van der Pol constr.
4.mL I
4.m -
-
1
\
i
4.30
4.30
I
I
I
I
I
0.00
1.M
2.00
3.00
4.00
I 1
0.00
5.00
Figure 1. Control variable profile for unconstrained van der Pol problem.
I 1.00
I 3.03
I
ZW
I
I
4.00
5.00
t
tirns
Figure 3. Control variable profile for constrained van der Pol problem.
-
Van der Pol Oscillator
Van der Pol constr.
0.w 1.00
0.80
-
0.70
-
0.60
-
0.M 0.40
-
0.30
-
0.10 0.20
4.00
-
4.10 4.20
-
4.30 4.40-
4.500.00
1.00
2.00
3 .00
4.00
5.00
-.-
1 1 0.00
I 1.00
I
I
zoo
3.00
I 4.00
5.00
Figure 2. Differential variable profiles for unconstrained van der Pol problem.
Figure 4. Differential variable profiles for constrained van der Pol problem.
which compares well with the value of 2.868 reported by Gritsis (1990). The corresponding profiles of the control and differential variables are shown in Figures 1 and 2, respectively. We now consider the solution of the constrained problem. For comparison purposes, we solve this problem with three different ways of handling the inequality path constraint (27). In particular, we employ (a) the hybrid technique of section 3.2,(b) the integral violation measure (cf. eqs 24)only, and (c) the discretized constraint at the 11 stage (i.e., control interval) boundaries (cf. eq 23)only. A value of e = lo-' (cf. eq 24b) was used for (a) and (b). Run c fails to produce a solution that satisfies the constraint within the control intervals themselves. On the other hand, the other two runs produce very similar results, with the path constraint being satisfied everywhere.
The optimal values of the objective function (2.95474and 2.954 36 for runs a and b, respectively) again compare well with the value of 2.960 reported by Gritsis (1990). The corresponding profiles of the control and differential variables are shown in Figures 3 and 4,respectively. Detailed computational satistics for the unconstrained and the three constrained runs are shown in Table 1. By comparing the results for the two constrained runs a and b that produce valid solutions to the problem, it can be seen that the use of the hybrid technique greatly improves the performance of the algorithm-in fact, the computational cost of (a) is comparable to that of (c) while ensuring solution feasibility. A somewhat less expected result is that the cost of solving the constrained problem using the hybrid technique is less than that of solving the (theoretically simpler)
2128 Ind. Eng. Chem. Res., Vol. 33, No. 9, 1994 Table 1. Computational Statistics for van der Pol Oscillator Problem. run QP LS 20 30 unconstrained constrained (a) 16 25 constrained (b) 34 51 constrained (c)b 15 20
Table 2. Kinetic and Thermodynamic Data for Batch Reaction Problem CPU 105.0 99.8 217.4 94.4
QP, number of quadratic programming subproblems solved by the optimization; LS, number of function evaluations during line search; CPU, CPU seconds on SUN SPARCStation 2 workstation. b Solution fails to satisfy path constraint.
unconstrained problem. This effect has, in fact, been observed in a large part of the extensive test problem set examined by Vassiliadis (1993b). One possible explanation is that the incorporation of the discretized constraints (23) restricts the feasible region, thus limiting the optimization search for some problems. 4.2. Optimization of Batch Reactor Operation. Here we consider the optimal operation of a batch reactor used to carry out the chemical reaction A+B-C in the presence of the side reaction B+C-D Both reactions are strongly exothermic, and therefore direct mixing of the entire necessary amounts of the reactants A and B must be avoided. Instead, the operation is performed by initially charging an amount of A in the reactor vessel, with B being added slowly during the reaction. The reactor vessel is fitted with a cooling water jacket which is used to remove the generated heat of reaction. Under a suitable set of assumptions, the transient behavior of the above system is described by the following mathematical model: r;lA= -Vr,
(28a)
reactionj 1 2
Reaction Data Ej (K) 0.008 3000.0 0.002 2400.0
Aj (m3/(mol-s))
A H j (kJ/mol)
-100.0 -75.0
Component Data (mol/m3) ai (kJ/(moEK)) @i (kJ/(mol.K2)) 11250 0.1723 0.474 X 10-9 B 16 OOO 0.2000 0.500X 10-9 C 10 400 0.1600 0.550X 10-9 D 10 OOO 0.1550 0.323 X 10-9 T,r = 298 K hp = 20 kJ/mol
component i A
pio
of components i (=A, B, C or D) in the system, while MT denotes the total molar holdup. The total energy holdup is denoted by H and the temperature by T. The rates of the two reactions are denoted by rj, j = 1, 2, and the corresponding rate constants by kj. The volume of liquid in the system is denoted by V , and an ideal mixture assumption is used to calculate the liquid density and enthalpy (cf. eqs 28i and 281). The values of the necessary kinetic and thermodynamic constants are given in Table 2. At the start of processing, the reactor vessel contains 9000 mol of pure A at a temperature of 350 K. It is desired to operate the above system so as to maximize the production of component C by manipulating the rate of addition of B and the external cooling load, i.e., FBand Q. Upper bounds of 10 mol/s and 1000 kW are imposed on these two control variables, respectively. In order to ensure that the reactions are sufficiently quenched without excessive subcooling, the final temperature must lie in the narrow range [295 K, 300 Kl. However, the overall duration of the reaction is left free to be determined by the solution of the optimization problem which can be written as m u
tf&dt),Q(t),tE[o, td
MC(tf)
(29a)
subject to the DAE constraints (28) and also 0 IFB(t) I10, t E [o, tf]
(29b)
0 IQ(t)5 1000,t E [O, tfl
(29c)
295 IT(tf) 5 300 (29d) For safety reasons, we also impose an upper bound of 520 K on the temperature attained during the reaction. This introduces the path constraint
V =MAa +MBT +MCa +MDT PA
PB
PC
(28i)
PD
In the above model, Mi, Xi, and Ci denote respectively the molar holdup, mole fraction, and molar concentration
T(t) 5 520, t E [0, tf] (29e) The above problem was solved in DAEOPT without and with the path constraint (29e). In both cases, five control intervals were employed, with piecewise constant controlvariables FB(t) and Q(t). The length of the control intervals was bounded between 60 and 1800 s. The optimization was initialized with each of the five control intervals having a length of 200 8, with &(t) = 5 mol/s and Q ( t )= 0 for all t in the time horizon. The solution of the unconstrained problem results in an optimal processing time of 2690 8, with the final amount of C being 3804 mol. The control variable profiles are shown in Figure 5. The optimal operating policy essentially involves adding a constant flow rate of B for approximately the first half of the processing time. The supply of B is then switched off, and a practically constant cooling load
Ind. Eng. Chem. Res., Vol. 33, No. 9,1994 2129 Batch Rcoctor I10 3
9m 8m
im 6m
sm 4m 3m 100
im
......._...... J. _ ......._... J.__. A.._.__. i.....-...-.i _..-... b.10 i
.I
i
om S
om
1.m
100
~ Q I o l m o o u o z m u o Batch Reactor
(a) Feed rate of component B Batch Reactor 6ylm
_..........................
.... :"'
.
aom
~
p-'
nom
p-'
smm
F-.
I
I
a m
! .-a !
p-.
-7I
L...
mm
i ! '-7 j
I-.
..-.
i
-7
C".
...d
som
xom
i I
j
i
:...A
................i....._...............-..._................. 1.......-......i!L--i ..........!hn 10 am
as
I .
! ;
IAO
is
am
3
YS. O
(b) Cooling load Figure 5. Control variable profiles for unconstrained batch reaction problem.
is applied to bring the temperature of the reactor contents down to the required level at the end of processing. Figure 6 shows the evolution of the molar holdups and the temperature. It can be seen clearly from the former that the reactions effectively terminate as soon as the addition of B is stopped. Also the reactor temperature satisfies the end-point constraint (29d), but exceeds the safety limit of 520 K over part of the horizon. We now re-solve the above problem, this time imposing the path constraint (29e) handled (a) by the hybrid technique of section 3.2, and (b) by the simple integral violation measure on ita own without the use of any constraint discretization. The objective function values obtained by the two techniques are identical to four significant digits, corresponding to the production of 3630 mol of C. Thus the imposition of an upper limit on the
om 1.00 am Figure 6. Molar holdup and temperature profiles for unconstrained batch reaction problem.
operating temperature causes a decrease of about 4.6% in the amount of C produced for a constant initial charge of A. The optimal processing time is also reduced to 1936 s. Figure 7 shows the optimal control variable profiles, and Figure 8 shows the reactor molar holdups and temperature. It can be seen that, in this case, there is a period with simultaneous addition of B and a nonzero cooling load, which practically coincides with the period during which the inequality path constraint is active. Some computational statistics for the three runs performed are shown in Table 3. The significant computational advantage of the hybrid technique (a) over its simpler counterpart (b) is again evident. 4.3. Steady-State Optimization of Tubular Reactor Design. This problem, which was also considered by Vasantharajan and Biegler (1990),involves two reactants, A and B, fed into a tubular reactor of length L in the ratio 1:3 to undergo an exothermic reaction. The reactor is jacketed, and the heat removed is used to produce steam. The feed stream is preheated through heat exchange with the exit stream, as shown in Figure 9.
2130 Ind. Eng. Chem.
Res.,Vol. 33, No.9,1994 Batch R u d o r
BStChRW .I0
i
I
1' r I
s
9 _i
Ec
4-
-.ii
I
!i
r+-
A
i -4I
I
b.10
om
OJO
im
1
zm
1-
Bakb R u d o r
am
Lo
, I
I=
u
@) cooling load Figure 7. Control variable profilea for constrained batch reaction problem. Table 3. Compntational Statistics tor Batch m i o n Problem. Nn unconstrained constrained (a) constrained (b)
QP
Is
34
35
954.8
42 53
48 60
1266.3
om OJO im 1W Figure 8. Molar holdup and temperature profdea for constrained bateh reaction problem.
hl
CPU 1664.8
QP. number of quadratic programming subproblems solved by the optimization; IS,number of function evaluations during line search;CPU,CPUsefondsonSUN SPARCStation 10141workstation. a
The variation of the conversion q along the length of the reactor under steady-state conditions is described by the equation q = 0.3(1- q)em1-'lh (30~) where q denotes the derivative of q with respect to the
--
10.
..c
Figure 9. Exothermic tubular reactor.
axial position z, and T is the normalized temperature, defined as the ratio of the temperature to the inlet temperature TR. T i s determined by the energy balance: T = 1.5(T- TJT,)
+ (2/3)q
(30b)
Finally, the rate of steam production is proportional to
Ind. Eng. Chem. Res., Vol. 33, No. 9,1994 2131 Table 4. Design Parameters for Reaction Problem init est 250.0 462.23 1.0 425.25
oarameter
TP TR L TS
lower bound 120.0 400.0 0.5 400.0
unconstrained problem this Vasantharajan work and Biegler 184.0 188.4 500.0 500.0 1.25 1.25 474.0 470.1 -171.4879 -171.2340
upper bound 1Ooo.o 500.0 1.25 500.0
obj function value
constrained problem this Vasantharajan work and Biegler 213.6 232.1 500.0 500.0 1.25 1.25 451.5 450.5 -152.8518 -148.0847
the temperature difference between the reactor and the steam jacket, and the associated profit is assumed to be given by
P = JL(TRT- T,) dz which can be written in differential form as
P = TRT- T, (30~) The initial conditions associated with equations (30) are q(0) = 0; T(0)= 1; P(0) = 0
(31)
The operation of the feed preheater is governed by a simple steady-state energy balance of the form
TR - TF = TRT(L)- Tp (32a) where the feed temperature, T F is , assumed to be 110 O C . Note that T(L) corresponds to the exit (normalized) temperature from the reactor. A minimum approach temperature of 10 "C is also specified through the constraint TRT(L)I TR + 10 (32b) The objective of the optimization is to determine the values of the four design parameters L, TR,Tp, and Ts that minimize the cost of the reactor. The latter takes account of the capital cost, assumed to be proportional to the length L, and the revenue P, leading to the objective function: min L - P(L)
L,TRTP,Ts
(33)
This constitutes a dynamic optimization problem with four time-invariant parameters and no control variables. The axial position z and the length of the reactor L correspond to "time" t and "final time" tf. The optimization was carried out using the starting point and bounds on the optimization parameters employed in the earlier study by Vasantharajan and Biegler (1990). These and the optimal values obtained at optimization respectively and integration tolerances of 10-5 and are shown in Table 4. The corresponding objective function is -171.4879. Vasantharajan and Biegler report an objective function of -171.4600; however, inserting the optimal values of the parameters that they report into eqs 30and integrating them a t an integration tolerance of yields a slightly worse objective of -171.2340. The variation of the conversion and the temperature along the lengthof reactor for the problem described above is shown in Figure loa, which clearly demonstrates the existence of a hot spot. In order to eliminate this hot spot, we consider the same problem with the additional path constraint
T(z)I 1.45,
E [O, Ll
(34) The optimal values of the design parameters for the constrained case are also listed in Table 4, while Figure 10b shows the axial variation of the conversion and z
(a) Unconstrained problem
(b) Constrained problem Figure 10. Axial conversion and temperature variation in tubular reactor.
temperature. Again the objective function actually corresponding to the optimization parameter values of Vasantharajan and Biegler is slightly different to that reported by those authors (-148.4730). Computational statistics for this problem are shown in Table 5. In comparison to the simple integral violation technique (b), the use of the hybrid technique (a) results
2132 Ind. Eng. Chem. Res., Vol. 33, No. 9, 1994 Table 5. Computational Statistics for Tubular Reactor Problem. run unconstrained constrained (a) constrained (b)
QP
LS
CPU
7 9 10
6 8 10
29.1 76.8 37.2
QP, number of quadratic programming subproblems solved by the optimization; LS, number of function evaluations during line search; CPU, CPU seconds on SUN SPARCStation 2 workstation. a
in adecrease in both the number of quadratic programming subproblems solved and the number of function evaluations during the line search. However, in this case, the decrease is outweighed by the increased cost of solving the quadratic programs, which now involve more constraints. This results in an increase in the overall CPU time. 5. Concluding Remarks
This paper has considered general methods for handling equality and inequality path constraints in the context of controlvector parametrization algorithms for the dynamic optimization of multistage DAE systems. The idea of treating the equality path constraints together with the underlying system equations is not in itself new. For instance, Bryson and Ho (1975, Chapter 3) consider means of imposing equality constraints of various types to systems described by ODES. Although the concept of DAE index was not known at that time, they propose a method based on repeated differentiations and algebraic manipulations that is effectively identical to reducing the index of the resulting DAE system to unity. However, the treatment presented in this paper is more general covering in a unified manner all types of equality constraints. It also differs in that it incorporates within the DAE system only those constraints that will not significantly increase the difficulty of its solution by increasing its index. The rest of the constraints are left to be handled by the optimization through conversion to end-point constraints. Of course, any end-point constraints derived from equality path constraints will still suffer from the problem of their gradients with respect to the optimization parameters being zero at the solution. However, it is interesting to note that these problems can be alleviated by adopting a hybrid approach similar to that for inequality constraints. Thus, in addition to the end-point constraints, we enforce each equality constraint as a point constraint at the stage boundaries. The proposed approach for handling inequality path constraints is also a hybrid of techniques proposed earlier in the literature, namely those based on constraint discretization at a finite number of points on one hand and forcing the integral constraint violation to zero on the other. As already explained, in their combined application, the former technique plays a key role in guiding the optimization search, while the latter ensures true feasibility of the solution obtained. In addition to the three examples presented here, further numerical evidence of the relative reliability and efficiency of the hybrid approach applied to both equality and inequality path constraints is presented by Vassiliadis (1993b).
Acknowledgment This research was supported by a joint SERC/AFRC grant.
Nomenclature A, = preexponential Arrhenius constant for reaction j Ak = matrix defined by eq 8 Bk = matrix defined by eq 9 c p = equality path constraints over stage k = inequality path constraints over stage k Ci = concentration of component i E, = activation energy for reaction j (expressed in K) f k = equations for stage k FB = feed flow rate of component B h~ = specific molar enthalpy of reactor feed stream hio = specific molar enthalpy of pure component i H = total energy holdup in reactor i = standard index for components (=A, B, C, D) I = unit matrix j = standard index for inequality path constraints k = standard stage index k, = rate constant for reaction j K = penalty parameter Lk = lower triangular matrix defined by eq 9 rn = number of algebraic variables y Mi = molar holdup of component i n = number of differential variables x NS = number of stages Pk = row permutation matrix defined by eq 9 q = number of time-invariant parameters u Q = reactor cooling load Q k = column permutation matrix defined by eq 9 rj = rate of reaction j rk = rank of matrix Ak s = slack variable t = time t o = initial time tf = final time t k = end time of stage k T = temperature T,,f = reference temperature for enthalpy calculations u = control variables a = control variables u,the time variation of which can take any functional form i i k = control variables u that will be determined by solution - of augmented DAE system in stage k C k = control variables u that will be determined by optimization in stage k U = domain of u uk = upper triangular matrix defined by eq 9 u = time-invariant parameters V = domain of u V = volume of reactor contents x = differential variables f = measure of integral constraint violation x i = mole fraction of component i f = time derivatives of differential variables x X = domain of x y = algebraic variables Y = domain of y
cr
Greek Letters
coefficient of linear term in pure component specific enthalpy expression (28m) p; = coefficient of quadratic term in pure component specific enthalpy expression (28m) AHj = enthalpy of reaction j e = end-point constraint relaxation parameter (cf. eq 24d) { = exponent of constraint violation in integral violation measure (18) x k = number of equality path constraints k k = number of inequality path constraints y = number of control variables u ?r = number of control variables 1 2 pio = molar density of pure component i a, =
Ind. Eng. Chem. Res., Vol. 33, No. 9,1994 2133
= points at which inequality path constraints are enforced in stage k (cf. eq 22)
7ik
Literature Cited Brenan, K. E.; Campbell, S. L.; Petzold, L. R. Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations; North-Holland, Elsevier Scientific Publishing Co.: New York, 1989. Bryson, A. E.; Ho, Y-C. Applied Optimal Control; Hemisphere Publishing Corporation: Washington, DC, 1975. Cuthrell, J. E.; Biegler, L. T. On the Optimization of DifferentialAlgebraic Process Systems. AZChE J. 1987,8, 1257-1270. Goh, C. J.;Teo, K. L. Control Parameterization: AUnified Approach to Optimal Control Problems with General Constraints. Automatica 1988, 24, 3-18. Gritsis, D. M. The Dynamic Simulation and Optimal Control of SystemsDescribedby Index Two Differential-AlgebraicEquations. Ph.D. Thesis, University of London, 1990. Jacobson, D. H.; Lele, M. M. A Transformation Technique for Optimal Control Problems with a State Variable Inequality Constraint. IEEE Trans. Autom. Control 1969, AC-14,457-464. Logsdon, J. S.; Biegler, L. T. Accurate Solution of DifferentialAlgebraic Optimization Problems. Ind. Eng. Chem. Res. 1989, 28, 1628-1639.
Sargent, R. W. H.; Sullivan, G. R. The Development of an Efficient Optimal Control Package. R o c . ZFZP Conf. Optim. Tech. 1977, 8th, part 2. Vasantharajan, S.; Biegler, L. T. Simultaneous Strategies for Optimization of Differential-AlgebraicSystems with Enforcement of Error Criteria. Comput. Chem. Eng. 1990,14, 1083-1100. Vassiliadis, V. S. DAEOPT-A Differential-Algebraic Dynamic Optimization Code, us. 2.0; Technical Report, Centre for Process Systems Engineering, Imperial College: London, 1993a. Vassiliadis, V. S.Computational Solution of Dynamic Optimization Problems with General Differential-AlgebraicConstraints. Ph.D. Thesis, University of London, 1993b. Vassiliadis, V. S.;Sargent, R. W. H.; Pantelides, C. C. Solution of a Class of Multistage DynamicOptimization Problems. 1. Problems without Path Constraints. Ind. Eng. Chem. Res. 1994, preceding paper in this issue. Walsh, S. P. Integrated Design of Chemical Waste Water Treatment Systems. Ph.D. Thesis, University of London, 1993. Received for review October 11, 1993 Revised manuscript received May 26, 1994 Accepted June 15,1994.
* Abstract published in Advance ACS Abstracts, August 1, 1994.
