A comparison of collocation methods for solving dynamic optimization problems

I)

Computers chem. Engng Vol. 19, No.4, pp. ' 7 5 - 3~ 1 1995 Copyright © 1995 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0098-1354(94)00064-6 (X)9~-1354/95 $9.50 + O.oo

Pergamon

A COMPARISON OF COLLOCATION METHODS FOR SOLVING DYNAMIC OPTIMIZATION PROBLEMS D. TIEU,1 W. R. CLUETT 1 and A. PENLlDIS: I

Department of Chemical Engineering and Applied Chemistry, University of Toronto, Toronto, Ontario, Canada M5S IA4 2 Department of Chemical Engineering, University of Waterloo, Waterloo. Ontario, Canada N2L 3G I

(Received 12 February 1993; final revision received 23 May 1994; received for publication 7 June 1994)

Abstract-This paper describes and compares a global end-point collocation method with the SOCOLL method. The effects of the number of interior collocation points (N) and the collocation methods themselves are investigated using an example. A second example shows the effectiveness of the endpoint collocation method with a state constrained problem. For the first example, the overall deviation between the approximate state profile [generated via the nonlinear program (NLP)j and the true state profile (generated via simulating the original model with the optimal control policy) was smaller for the end-point collocation method. The results also showed a nonuniform improvement with respect to increasing N. For the second example, the end-point collocation method (N = 5) gave better results than those found in the literature.

1. INTRODUCTION

This paper describes and compares an end-point collocation method with the SOCOLL (Simultaneous Optimization and COLLocation) method (Biegler, 1984) for solving dynamic optimization problems. The end-point coIlocation method is different from the SOCOLL method in that the endpoint at If is also a coIlocation point. The approximate polynomials in the end-point collocation method are one degree higher than those used in the SOCOLL method. The use of the end-point in collocation formulae has been well-studied for the solution of ordinary differential equations (Villadsen and Michelsen, 1978). For these problems, the end-point collocation methods are preferable to collocation methods without end-point collocation due to improved stability. In the case of solving dynamic optimization problems, collocation formulae are used only to transform the ordinary differential equations into algebraic equations. If the end-point is not used, as in the case of the SOCOLL method, then extrapolation equations for the state variables must be included in the problem formulation. In the case of collocation on finite-element methods (Cuthrell and Biegler, 1987, 1989; Logsdon and Biegler. 1989), additional extrapolation equations have to be added onto the right-hand side of each element for control variables in order to obtain continuous control profiles. The primary objective of this paper is to compare a global end-point coIlocation method with the SOCOLL method using an example in order to

examine whether or not including the end-point provides any improvement to the solution. In particular, the effects of the number of interior collocation points (N) and the choice of collocation method are investigated. A second example, which has been solved by others [control parameterization method (Goh and Teo, 1988) and iterative dynamic programming method (Luus, 1991)1, is included in this paper to demonstrate the applicability of the endpoint collocation method to a state constrained optimization problem. Both examples, which have smooth and continuous profiles and therefore are amenable to global collocation methods, are solved using an optimization package (OPTPAC) developed by Tieu (1991). OPTPAC uses GAMS/MINOS (5.2) to solve the resulting discretizcd nonlinear programming (NLP) problem, and SIMNON (3.0) to test the NLP solution by integrating the original dynamic model with the "optimal" policy suggested by GAMS/MINOS. The default convergence tolerances in both packages were used . Additional applications of OPTPAC for solving minimum end-time problems involving batch polymerization reactors may be found in Tieu et al (1994).

2. END.POINT COl.LOCATION METHOD

A general differential-algebraic problem can be written as follows: max or min x,

375

U, 'i. /1

optimization

J P[x(t), u(t), s, ttl.

(I)

376

D. TIEV et al.

s.t. g[x(t) , u(t), sJ~o,

(2)

8fl=Xo,i=1,2, . . . , ( N + I ).

h[x(t), u(t), sJ = 0 ,

(3)

To work with a dimensionless time scale , one can multiply equation (15) by t, to obtain ;

(4)

r(T;)=

dx/dt = f[x(t) , u(t) , s],

N+ I

x(O) = x, (given),

(6)

~x~x· ,

(7)

x.

In order to transform the differential equations into algebraic equations , the state and control variable vectors are approximated by vectors of polynomials written in Lagrange form, of degrees ~(N + 1) and ~N, respectively: N+I

XN+I(t) =

L ~¢;(t),

(8)

j ~(J

N+I

UN(t) =

L b;tp;(t) , j~

(9)

1

with IIfl=Xfl,

i=I,2, .. . ,(N+l ),

where (17)

q)j(T) =

L~t rj~TJ [p~lj (jJ.P (r-T/»)

l

rj

JP[a;, b., s, tr],

max or min (10)

s.t.

(11)

g[aj, b., s] ~o ,

(20)

h[a;, b., s] =

(21)

0,

L ajq)lT;)-t,f[a;, b., s] =0 ,

UN(t), s]

(22)

j = (J

with llt,=lC(h

N+l

f[XN+ I(t),

(19)

N+ I

r(T;) =

r'(r) = 2>jq) ;(t) -

(18)

Since the approximate polynomials XN+ 1 and UN [equations (8) and (9)], when evaluated at tj (or if written in the dimensionless time variable), become the corresponding coefficients aj and bj , the general differential-algebraic optimization problem can be rewritten as:

where

and k = 0, j denotes k starting from 0 and k =I:- j . Substitution of equations (8) and (9) into (4) yields residual equations:

(16)

j~n

(5)

u. ~u~u* ,

L ajq)lT;) - t,f[XN+I(T;), uN(rj) , 5]=0,

(23)

(12)

(24)

j~fl

with

(25) and

where

i=I,2, . .. ,(N+I) .

and 1=0, i . p denotes 1 starting from 0 and t=l:-j, I:/:.p. To discretize the residual equations, the method of collocation is used which requires:

f

The locations of rio i = 1,2, . .. , N in equation (22) are chosen to correspond to the shifted zeros of an orthogonal Legendre polynomial of degree N as used by Biegler ( 1984). The value of TN+ I in equation (22) is unity. 3. ILLUSTRATIVE EXAMPLES AND DISCUSSION

ll

r'(t)<5(t-t;) dt=O,

i= 1,2, ... . (N+ 1),

(14)

(J

where <5 is the Dirac delta function . Equation (14) becomes: N+I

r ' (t;) =

2: ajq)j(t;) -

f[XN+ I(t;),

UN(t;), s] = 0 (15)

To compare the end-point colloc at ion method with the SOCOLL method, an example is taken from lang and Yang (1989) . In the SOCOLL method , the following extrapolation equation for the state variables is required : N

j~(J

x(tr) = with

2: a,j¢,j(I), j~(J

(26)

:177

Comparison of collocation methods Plot 1 : State Profiles

s

N=3 Solid=Approx. Doued=Simu.

0.8

'"' S 0

Plot 2 : State Profiles

S

~ ......

0.6

80

N=5 Solid =Approx. Dotted=Simu.

0.8 0.6

u 0.4

0.4

6 o

o

0.5

1

0.5 Dimensionless Time

Dimensionless Time Plot 3 : State Profiles

1

S0

1000

o :N=3 x :N=5

N=1 Solid=Approx. Dotted=Simu.

g

0.6

~

0.4

Plot 4 : Control Profiles

'"' ~

......

!

U

700 0.5

0

0.5 Dimensionless Time

Dimensionless Time

Fig. 1. State profiles (N = 3.5,7) and control profiles (N = 3,5) obtained by the end-point collocation method for Example 1. The discrete points in the plots are the NLP solution.

whe re

(27)

Example 1 Th e problem is formulat ed as follows: max

CII(t f) ,

T(r)

whe re

k;=kiilexp(- E;lRT) , i = 1.2, k lll == 65.6 s- I ,

k7!.'= 1970 s- l , E, == 1.0 X 104 cal/ mol, E2 = 1.6 x 104 cal/mol, tf = 12.5 s.

Th is example was first solved using the end-point co llocation method describ ed in Section 2. The

approxi mate and simula tion state profiles for N = 3,5 and 7 (plots 1, 2 and 3) and the optimal contro l profiles for N = 3 and 5 (plot 4) a re show n in Fig. I . Th e co rrespo nding de viati on between the approximate profi les and the simulatio n pro files (plots I , 2 and 3) a nd th e optima l cont rol profi les for N = 5 and 7 (plot 4) are given in Fig . 2. In Fig. I an d plot 4 of Fig. 2, the discrete points are the NLP so lutio ns. The se points are used in generating the Lagrange polynomials. These polynomials are used to plot the solid lines (approximate profiles) which connect all of the discrete point s on the gra phs. The do tte d lines (plots 1. 2 and :I of Fig. l ) are the simulation pro files which result from simulating th e origi na l mod el using the o ptimal co ntro l policies (plots 4 of Figs 1 and 2). The purpose of co mpa ring the ap proxi mate and simulatio n profiles is to check the va lidity of the NLP results. It can be seen that the fit betwee n the approximate and the simulation profiles improves as N increa ses from 3 to 5 to 7 (Fi g. 1) . Th e impro veme nt in the fit is mor e evident as sho wn in plot s 1. 2 and 3 of Fig. 2. The inte grals of the absolu te erro r betwe en the approxi mate and simulated profiles (IA E) for plot s I , 2 a nd 3 of Fig. 2 are tabul ated in Table 1. It can be seen that the IA E de creases as N increases fro m 3 to 5 to 7. However . a larger N did not necessar ily mean better agree me nt between the approximat e and the simulation pro files because

378

D. T IEu et al.

N=5

.... 0.01

S ~ .....,

>

8 -0.01 : -0.04' - - - - - - - ' - - - - - - - - - ' o 0.5 Dimensionless Time

1000

A

Plot 4 : Control Profiles

N=7

x:N=5 +:N=7

....

....0.005

S0

~ .....,

,.' - < ,

0

>

Dimensionless Time

Plot 3 : Dev. A

0.01

5

-0.02 '------~'---------I a 0.5

c:i.

E ~

~-0.005

700

-noi

0

0.5

0

0.5 Dimensionless Time

Dimensionless Time

Fig. 2. Deviations of the approxi mate from the simulation sta te profiles in Fig. 1. T he o ptimal contro l profiles (N =5.7) are shown in plot 4.

these profiles did not match wel1 for N = 4 and 6 (not shown). For comparison , the above proble m was also solved using the SOCOLL method for N =3 , 5 and 7. For each N , the initial starting values and bound s of the variable s in the NLP problems were the same as for the end-point col1ocation method . The approximate and simulation state profiles for N = 3, 5 and 7 (plots 1, 2 and 3) and the optimal control profiles for N = 3 and 5 are shown in Fig. 3, while the corresponding deviation between the approximate and simulation profiles (plots 1, 2 and 3) and the optim al control policies for N= 5 and 7 (plot 4) are given in Fig. 4. The corresponding IAE are also tabulated in Table 1. It can be seen that for both methods the best approximation is N = 7. For the same N, the IAE is

smal1er for the end-point collocation meth od in all three cases. Table 1 also tabulates CB(tt) obtained by the two methods. To decide which method provides the best result for this example , one may make the final selection based on the following thr ee criteria: (a) highest simulation value of CB(tt), (b) lowest percentage difference of Ca(tt) (NLP) relative to the simulation result ; (c) best overall fit of the approximate to the simulation state profiles (lowest IAE, Table 1). The simulation value of CB(tt) is suggested as an indicator in criterio n (a) instead of the approximate value of CB(tr) since the simulation value corresponds to the actual Ca(tr) that would be obtained . Criteria (b) and (c) are suggested to help screen out

Table I. IAE values and performance index for Example I End-point collocation

SOCOLL

Integral of absolute error N

A

B

Total

3

7.112e-3

1.015c-2

1.797e-2

5

2.6lIc-3

3. 12e·3

5.lIlle-3

7

I. l3e-3

1.32e-3

2.45e-3

lntcgral of absolute error CH(lfl"

A

B

Total

CH(lr)'

U.4110036 (1l.471197I1l U.4791172 (0.479862) 0.479996 (0.48(X147)

1.573e-2

1.1I42e-2

3.4 I5e-2

3.IIUe-3

4.57e-3

1I.37e-3

J.54e-3

1.85e-3

3.3ge-3

0.4112446 (U.475693) 1l.4110109 (0.479604) 0.4110048 (0.479993)

" # = NLP or approximate results: ( # ) = simulation results.

1,79

Comparison of collocation methods

Plot I : State Profiles 0.8

S

~ '-'

0.6

~ u

0.4

0

Plot 2 : State Profiles

1

N=3 Solid=Approx. Dotted=Simu.

N=5 Solid=Approx. Dotted=Simu.

0.8

S ~

0.6

~

0.4

'-'

_.... -

u

0.2

B

0.2 0 0

0.5

0.5 Dimensionless Time

Dimensionless Time Plot 3 : State Profiles

1 A

0.8

0:

0.6

~

U

0.4

eu

N=~

x:N=

......

S0 e

c:i.

'-'

I:: 0

Plot 4 : Control Profiles

1000

N=7 Solid=Approx. Dotted=Simu.

!-'

u

0.5

0

0.5

Dimensionless Time

Dimensionless Time

Fig. 3. State profiles (N= 3.5,7) and control profiles (N= 3,5) obtained by the SOCOLL method for Example I. The discrete points in the plots are the NLP solution.

A

i

N=5

...... 0.01

Or---t't------::X-:---?'(,J

.

'-'

B-0.02 ""

N=3

-0.04 \_,/ B -0.06 '---

o

---'0.5

S

-g Ot---r---:,i------"-<-------:~_"I B-0.01, ,, ...

......

-,

~

",/ B

----J

-o.02'---------'~--­

o

Dimensionless Time

0.5 Dimensionless Time Plot 4 : Control Profiles

A

N=7

0.005

x:N=5 +: N=7

g

..

1

, ,

! ", "

-0.005,

-0.01 '----

o

700 --'-0.5

Dimensionless Time

---J

o

0.5 Dimensionless Time

Fig. 4. Deviations of the approximate from the simulation state profiles in Fig. control profiles (N = 5,7) are also shown in plot 4.

1,.

The optimal

380

D. TIEu et

at.

State Profiles

0.8 N

0.6

.....I<

0.4

><

0.2 0

0.1

0.8

0.9

-0.5 L..-_----'--_ _-'--_----'--_ _-'--_----'-_ _-'--_----'-_ _ 0.4 0.5 0.6 0.7 0.8 0.1 o 0.2 0.3

0.9

0

0.2

0.3

0.4

0.5

0.6

0.7

TIME

ControlProfile

0.5

:l

0

..I...-_----L.._--J

TIME Fig. 5. State and control profiles obtained by the end-point collocat ion method (N = 5) for Example 2. The discrete points in the plots are the NLP solution.

large discrepancies between the approximate and simulated results. Criterion (b) is concerned with the NLP solution at t.; while criterion (c) is concerned with the overall fit. Using the above three criteria, the best result is obtained by the end-point collocation method (N = 7). The simulation value of the performance index CB(tt) = 0.480047 is also slightly higher than the literature result (0.47991) obtained by the control vector iteration method in lang and Yang (1989).

Example 2 The problem is formulated as follows: min X2(tt). U( /)

dX2 dt

-=X

2+U 2 I

,

xo.(O) =0.

The final state constraint at tt = 1 is:

Xt(1) = 1. Using the end-point collocation method , the approximate and simulation state profiles for N = 5 and the corresponding optimal control profile are shown in Fig. 5. The performance index values obtained by the end-point collocation method

(N = 5) were 0.92423 (NLP) and 0.92425 (simulation). These values are lower (better) than those obtained by Goh and Teo (0.92518) and Luus (0.92441).

4. CONCLUDING REMARKS

This paper has described and compared an endpoint collocation method with the SOCOLL method (Biegler , 1984) using an example. In this example, the overall deviation of the approximate (NLP) state profiles from the simulation state profiles was smaller for the end-point collocation method than for the SOCOLL method with the same number of interior collocation points N. Within the same collocation method , the overall deviation decreased as N increased from 3 to 5 to 7 (Table 1). However , the impro vement in performance was nonuniform with respect to N. It should be noted that these observations have been derived from an example which has smooth and continuous profiles. In a second example with a terminal state constraint , the end-point collocation method gave a bett er solution than those reported in the literature (Goh and Teo, 1988; Luus , 1991). The performance index obtained by the end-point collocation method (N = 5) for the minimization problem was smaller, and there was no constraint violation in the solution.

3ill

Comparison of collocation methods NOMENCLATURE 8i = jth

column vector of coefficients of state variables hi = jth column vector of coefficients of control variables C A , C B = Concentration of species A and B, respectively f= n X I vector of the derivative of x g = Inequality constraint vector, p x I h = Equality constraint vector, q x I J P = Objective function k" k 2 = Reaction rate constants N = Number of interior collocation points r= Vector of residual equations, n x I R = Gas constant (cal fmol-K) s = Decision variable vector t=Time T=Temperature (K ) u = Control variable vector, m x l x = State variable vector, II X 1 Greek

b = Dirac Delta functi on r=Dimensionless time (1I I r) Subscripts • = Lower bounds f = Final time N, N + 1 = Denote Nand (N + 1) degree polynomials. respectively 0= Initial time (except k",) Superscripts , = Denotes the dim ensional time versions of r • = Upper bounds

REFERENCES

Biegler L. T. Solution of dynamic optimization problems by successive quadratic programm ing a nd orthogona l collocation. Com puters chem . Engng 8, 243-247 (19il4) . Cuthrell J . E. and L. T . Biegle r, On the optimiza tio n of differential-algebraic process systems . AIChE JI 33. 1257-1270 (1987). Cuthrell J. E . and L. T. Biegler, Simultaneou s optimization and so lutio n methods for batch reactor control profiles . Comp uters chem. Engng 13. 49-62 ( 19il9). Goh C. J . and K. L. Teo. Control parameterization: a unified approach to optimal control problems with general constraints. Automatica 24, 3-18 (1988). Jan g S. S. and W. L. Yang , D ynamic optimization of batch emulsion polymerization of vinyl acetate-an orthogonal polynomial initiator policy . Chem. Engng Sci . 44, 515-528 (1989). Logsdon J . S. and L. T. Biegler , Accurate solution of differential-algebraic optimization problems . Ind. Engng Chem . Res. 28, 1628-1639 (1989). Luu s R . Application of iterative dynamic programming to sta te constrained optimal control problems, July 1991. Presented at the Int. Workshop on Chemical Eng ineering Mathematics, Gottingen, Germany (1991). Tieu D ., A dynami c optimization package and the appli cat ion o f an e nd -po int collocation method to batch pol ymer reactors. Master's Thesis, Department of Che mical Enginee ring and Applied Chemistry . University of Toronto ( 199 1). Ti eu D ., W. R . Cluett and A . Penlidis, Optimization of pol ymerization reactor operation: review and case stu dies with the end-point collocation method . Poly . React . Engng J, 2. 275-313 (1994). Villadsen J . a nd L. M. Michelsen , Solution of Differential Equation Model s By Polynomial Approximation . Prentice-Hall , E nglewood Cliffs. NJ (1978) .