An Improved Control Vector Iteration Approach for Nonlinear Dynamic Optimization (I) Problems Without Path Constraints

PROCESS SYSTEMS ENGINEERING Chinese Journal of Chemical Engineering, 20(6) 1053—1058 (2012)

An Improved Control Vector Iteration Approach for Nonlinear Dynamic Optimization (I) Problems Without Path Constraints* HU Yunqing (胡云卿)1, LIU Xinggao (刘兴高)1,** and XUE Anke (薛安克)2 1 State 2

Key Laboratory of Industry Control Technology, Zhejiang University, Hangzhou 310027, China Institute of Information and Control, Hangzhou Dianzi University, Hangzhou 310018, China

Abstract This study proposes an efficient indirect approach for general nonlinear dynamic optimization problems without path constraints. The approach incorporates the virtues both from indirect and direct methods: it solves the optimality conditions like the traditional indirect methods do, but uses a discretization technique inspired from direct methods. Compared with other indirect approaches, the proposed approach has two main advantages: (1) the discretized optimization problem only employs unconstrained nonlinear programming (NLP) algorithms such as BFGS (Broyden-Fletcher-Goldfarb-Shanno), rather than constrained NLP algorithms, therefore the computational efficiency is increased; (2) the relationship between the number of the discretized time intervals and the integration error of the four-step Adams predictor-corrector algorithm is established, thus the minimal number of time intervals that under desired integration tolerance can be estimated. The classic batch reactor problem is tested and compared in detail with literature reports, and the results reveal the effectiveness of the proposed approach. Dealing with path constraints requires extra techniques, and will be studied in the second paper. Keywords nonlinear dynamic optimization, control vector iteration, discretization

1

INTRODUCTION

Many chemical engineering processes are essentially dynamic. However, traditional optimizations of these processes are generally based on their steady-state models, which can not satisfy the increasing competition and tighter operation in modern chemical plants. As a result, dynamic optimization methods plays a growing important role for providing the optimal operation strategies [1-4]. The mathematical model of a dynamic process is generally described by a set of ordinary differential equations (ODEs). Dynamic optimization aims to minimize or maximize some specified performance indexes of a dynamic process by determining the optimal operation strategy, while ensuring product quality, operating safety, environmental requirements and other necessary constraints. Numerical methods for solving nonlinear dynamic optimization problems can be classified into two categories: direct and indirect methods. Direct methods discretize the original infinite-dimensional dynamic optimization problem, resulting in a finite-dimensional static optimization problem that can be solved by many available nonlinear programming (NLP) solvers [5]. Indirect methods achieve the solution by solving the optimality necessary conditions, which are taken as a two-point boundary value problem (TPBVP). Most indirect methods need a rather precise initial guess for the control, state and adjoint variables, some indirect methods even require priori knowledge about the structure of the optimal solution. Unfortunately, none of these requirements is easy to satisfy in practice [6].

Nevertheless, this paper proposes a novel indirect approach that can conquer all the above drawbacks. The approach solves the optimality conditions like the traditional indirect methods do, but uses a discretization technique which is inspired and introduced from direct methods. After discretizing the control variables, the original dynamic optimization problem is reduced to an unconstrained NLP problem, and then iteratively solved by the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm. The function evaluation required by BFGS algorithm is performed by solving an initial value problem (IVP) of ODEs, and the gradient information (the gradients of the objective function with respect to the control variables) can be proved to be equal to the gradients of Hamiltonian with respect to the control variables. In order to improve the computational efficiency, we use the four-step Adams predictor-corrector algorithm to solve the IVPs. Furthermore, the relationship between the minimal number of discretized time intervals and the integration error of the four-step Adams predictor-corrector method is established, so that the control and state profiles can be refined when the result is not satisfactory. 2

PROBLEM FORMULATION

The mathematical statement of the nonlinear dynamic optimization problem (P1) considered in this study is formulated as tf

min J [u(t )] = Φ [ x (tf ), tf ] + ∫ L[u(t ), x (t ), t ]dt (1a) t0

Received 2012-04-12, accepted 2012-07-21. * Supported by the National Natural Science Foundation of China (U1162130), the National High Technology Research and Development Program of China (2006AA05Z226), and the Outstanding Youth Science Foundation, Zhejiang Province (R4100133). ** To whom correspondence should be addressed. E-mail: [email protected]

1054

Chin. J. Chem. Eng., Vol. 20, No. 6, December 2012

s.t.

dx = F [u(t ), x (t ), t ] dt

(1b)

x (t0 ) = x0 L

(1c)

u ≤ u(t ) ≤ u

U

(1d)

t0 ≤ t ≤ t f

(1e)

where x (t ) ∈ R n is the vector of state variables,

differential- algebraic system which is very difficult to solve. If Eq. (1d) is assumed to be absent here, Eq. (2) is then reduced to Eq. (4), H [u(t ), x (t ), λ (t ), t ] = L + λ (t )T F

and Eqs. (3a)-(3f) are reduced to much simpler Eqs. (5a)-(5c): dx ∂H = = F [u(t ), x (t ), t ] x (t0 ) = x0 (5a) dt ∂λ (t )

u(t ) ∈ R r is the vector of control variables; uL ∈ R r U

3

APPROACH DESCRIPTION

3.1

dλ ∂H =− dt ∂x ( t )

r

and u ∈ R are the lower and upper bounds on u(t) respectively; the initial time t0 and the final time tf of the process are both fixed; the nonlinear dynamic process is described by the vector function F, which is differentiable with respect to u(t ) and x (t ) under the local Lipschitz condition; the consistent initial conditions are given by x0 .

Optimality necessary conditions for P1

The Hamiltonian of P1 is defined as H [u(t ), x (t ), λ (t ), μ (t ),ν (t ), t ] = L + λ (t ) F + μ (t )T [uL − u(t )] + ν (t )T [u U − u(t )] (2) T

where λ (t ) ∈ R n is the vector of adjoint variables for F, μ (t ) ∈ R r is the vector of adjoint variables for the

path constraints uL ≤ u(t ) , and ν (t ) ∈ R r is the vector of adjoint variables for u(t ) ≥ u U . The optimality necessary conditions for P1 can be derived from Biegler [1], and are expressed by Eqs. (3a)-(3f): dx ∂H = dt ∂λ (t ) dλ ∂H =− dt ∂x ( t )

x (t0 ) = x0

λ(tf ) =

(3a)

∂Φ[ x(t), t] ∂x(t) t =t

f

∂H =0 ∂u(t )

μ (t ) ⊥ [uL − u(t )] = 0 uL − u(t ) ≤ 0

(3b)

ν (t ) ⊥ [u U − u(t )] = 0 u U − u(t ) ≤ 0

μ (t ) ≥ 0 ν (t ) ≥ 0

(3c) (3d) (3e) (3f)

where H, L, and F represent H [u(t ), x (t ), λ (t ), μ (t ), ν (t ), t ] , L[u(t ), x (t ), t ] , and F [u(t ), x (t ), t ] respectively. Clearly, Eqs. (3a)-(3f) constitute a complicated

(4)

λ(tf ) =

∂Φ[ x(t), t] ∂x(t) t =t

f

∂H =0 ∂u(t ) 3.2

(5b) (5c)

Control vector iteration and discretization

Control vector iteration (CVI) is a classic indirect approach for dynamic optimization problems without path constraints. Its idea is to keep Eqs. (4), (5a), (5b) exactly satisfied during the solution process, while the satisfaction of Eq. (5c) is achieved gradually in an iterative manner. When Eq. (5c) is finally satisfied, the optimal operation strategy for P1 is obtained. For more details, please refer to [7]. In this study, the structure of the traditional CVI is improved in the following aspects: first, just unconstrained NLP algorithms (e.g. BFGS algorithm) are enough to get the solution; second, a control variable discretization technique is introduced to make the solution process numerical as expressed by Eqs. (6) and (7): tk =

t f − t0 k, N

k = 1, 2," , N

u(t ) ≈ u = [u1 , u2 ,", uN ]

(6) (7)

where N represents the number of discretized time intervals and uk represents the value of u(t ) over the kth time interval; third, the control and state variable profiles can be refined automatically when the result is beyond a pre-specified tolerance; last, problems with path constrains can also be solved in its framework, and will be discussed in the second paper. Note that tk can be taken as the integration points directly if N is large enough. Moreover, if the integration of Eq. (5a) forwards and the integration of Eq. (5b) backwards share these same integration points, the result of the integration phase will be more reliable and stable. Although the path constraints expressed by Eq. (1d) do not appear in Eqs. (5a)-(5c), they should be satisfied in the solution. We apply the clipping technique to handle Eq. (1d) independently in this paper, as was used by previous researchers [8, 9]: if any of the control variable is beyond the constraints as imposed by Eq. (1d) , then set the control variable as the boundary value of the violated constraint in that time interval.

1055

Chin. J. Chem. Eng., Vol. 20, No. 6, December 2012

3.3 Function evaluation and gradient information calculation

3.4 Adams predictor-corrector algorithm and profiles refining

Suppose u(t ) is given, then x (t ) can be obtained by solving the state equations Eq. (5a) and the value of the objective function J [u(t )] can also be evaluated. Furthermore, λ (t ) can be obtained by solving the adjoint equations Eq. (5b) and the value of H can be obtained from Eq. (4). The gradients of the objective function J [u(t )] with respect to the control variable u(t ) is equal to the gradients of Hamiltonian with respect to u(t ) , as expressed in Eq. (8): ∂H ∇J [u(t )] = (8) ∂u(t ) In order to prove this formula, we just need to prove Eq. (9):

The high-precision four-step Adams predictorcorrector algorithm, which is a combination of the explicit Adams-Bashforth formula and the implicit Adams-Moulton formula, is chosen to solve Eqs. (5a), (5b). The time points tk in Eq. (6) are set as the integration points. The explicit formula Eq. (10): h xˆ (tk +1 ) = x(tk ) + {55F [ x(tk ), uk , tk ] − 24 59F [ x(tk −1 ), uk −1 , tk −1 ] +

tf 0

T

⎛ ∂H ⎞ ⎟ δ u dt ⎝ ∂u ⎠

δ J ( u) = ∫ ⎜ t

(9)

Proof Since the objective function J [u(t )] in Eq. (1a) can be equivalently expressed by J [u(t )] = tf

Φ [ x (tf ), tf ] + ∫ {H − λ (t )T F }dt the variation of J[u(t)] t0

is

δ J [u(t )] T

tf ⎡ ⎡ ∂Φ ⎤ ⎛ dx ⎞ ⎤ T =⎢ ⎥ δ x (tf ) + ∫t0 δ ⎢⎣ H − λ (t ) ⎜⎝ dt ⎟⎠ ⎥⎦ dt ∂ x ( t ) f ⎦ ⎣ tf tf ⎡ ⎛ dx ⎞ ⎤ = ∫ d[λ Tδ x (t )] + ∫ δ ⎢ H − λ T ⎜ ⎟ ⎥ dt t0 t0 ⎝ dt ⎠ ⎦ ⎣ T ⎤ tf ⎡⎛ dλ ⎞ T ⎛ dx ⎞ = ∫ ⎢⎜ ⎟ δ x (t ) + λ ⎜ ⎟ ⎥ d t + t0 ⎣⎝ dt ⎠ ⎝ dt ⎠ ⎦ T T ⎫⎪ tf ⎧ ⎪ ⎡ ∂H ⎤ ⎛ ∂H ⎞ + x δ ( t ) ⎜ ∂u(t ) ⎟ δ u ⎬ dt + ∫t0 ⎨⎩⎪⎢⎣ ∂x (t ) ⎥⎦ ⎝ ⎠ ⎭⎪ tf

∫t

0

=∫

tf

t0 tf

∫t

0

⎧⎪ ⎡ ∂H ⎤ ⎫ ⎛ dx ⎞ ⎛ dx ⎞ δλ (t ) − λ (t )T δ ⎜ ⎟ − ⎜ ⎟ δλ (t ) ⎬ dt ⎨⎢ ⎥ ⎝ dt ⎠ ⎝ dt ⎠ ⎭ ⎩⎪ ⎣ ∂λ (t ) ⎦ T

T

T ⎧⎪ ⎡⎛ dλ ⎞ ∂H ⎤ T ⎫⎪ ⎡ ∂H ⎤ x ( t ) δ δ u(t ) ⎬ dt + + + ⎨ ⎢⎜ ⎟ ⎥ ⎢ ⎥ ⎣ ∂u(t ) ⎦ ⎩⎪ ⎣⎝ dt ⎠ ∂x (t ) ⎦ ⎭⎪

⎧⎪ ⎡ ∂H ⎛ dx ⎞T ⎤ ⎫⎪ − ⎜ ⎟ ⎥ δλ (t ) ⎬ dt ⎨⎢ ⎩⎪ ⎣ ∂λ (t ) ⎝ dt ⎠ ⎦ ⎭⎪ T

tf ⎡ ∂H ⎤ δ u(t )dt =∫ ⎢ 0 ∂u(t ) ⎥ ⎦ ⎣

Eq. (8) is thus proved. Consequently, the gradients of the objective function J [u(t )] with respect to the control variable u(t ) can be obtained by calculating the gradients of Hamiltonian with respect to u(t ) .

37F [ x(tk −2 ), uk −2 , tk −2 ] −

(10)

9F [ x(tk −3 ), uk −3 , tk −3 ]}

predicts an approximation of x (tk +1 ) [denoted by xˆ (tk +1 ) ], then the implicit formula Eq. (11) corrects this prediction: h {9F [ xˆ (tk +1 ), uk +1 , tk +1 ] + 24 19 F [ x (tk ), uk , tk ] − 5F [ x (tk −1 ), uk −1 , tk −1 ] +

x (tk +1 ) = x (tk ) +

F [ x (tk − 2 ), uk − 2 , tk − 2 ]} (11) The algorithm is not self-starting so the fourth-order Runge-Kutta algorithm is used for initialization. From Eq. (6), the time step size h can be expressed as t −t h= f 0 (12) N The local truncation error of Eq. (10) and Eq. (11) are both fifth-order, as provided by Eq. (13) and Eq. (14) respectively [10]: ⎛ 251 ⎞ ⎛ d x OAdams-Bashforth (h5 ) = ⎜ ⎟⎜ ⎝ 720 ⎠ ⎜⎝ dt 5 t =ξ

⎞ 5 ⎟⎟ h ⎠

(13)

5 ⎛ 19 ⎞ ⎛ d x OAdams-Moulton (h5 ) = ⎜ − ⎟⎜ 5 ⎝ 720 ⎠ ⎜ dt t =ξ ⎝

⎞ 5 ⎟h ⎟ ⎠

(14)

5

where ζ is a constant belongs to the range of [tk −3 , tk +1 ] , ξ is a constant belongs to the range of [tk − 2 , tk +1 ] . The global truncation error of the four-step Adams predictor-corrector algorithm is of one order lower than its local truncation error [10], that means, it is fourth-order precision, and can be written as Eq. (15): OThe four-step Adams (h 4 ) = η h 4

(15)

where η is a constant depends on the fifth derivative of x(t) with respect to t in [tk −3 , tk +1 ] . The global error goes to zero as the step size h goes to zero, in contrast, the number of time intervals N will be very large and produce an extremely highdimensional NLP problem. By taking the expression

1056

Chin. J. Chem. Eng., Vol. 20, No. 6, December 2012

of Eqs. (12) and (15) into account, the minimum of N can be estimated: 4

⎛ t f − t0 ⎞ (16) ⎟ ≤ ε1 ⎝ N ⎠ where ε1 is the integration tolerance. In the proposed approach, η does not have to be the accurate value; the value of η is assigned previously, if the result is not satisfactory, a larger value will be assigned to η again.

η⎜

3.5

Figure 2

The batch reactor problem

Algorithm outline

The structure of the improved CVI approach is sketched in Fig. 1, the corresponding numerical algorithm is stated as follows: Step 1 Set the tolerances ε1 for Adams algorithm, ε 2 for BFGS algorithm, and η, choose the time interval number N according to Eq. (16). Step 2 Assign the initial guess for control variables u(0) according to experience or priori knowledge about the problem. Set the number of iteration k = 0 . Step 3 Integrate the state equations Eq. (5a) forwards to get the initial state trajectory x ( k ) (t ) , and then integrate the adjoint equations Eq. (5b) backwards to get the initial adjoint trajectory λ ( k ) (t ) , thus J [u( k ) ] and ∂H ( k ) / ∂u( k ) are obtained. Step 4 Optimize the resulting NLP problem, Eq. (1a), by BFGS algorithm where the gradient trajectory is ∂H ( k ) / ∂u( k ) , and handle the new control variables

by the clipping technique, in this way a set of new control variables unew are obtained. Step 5 Repeat the integration phase for Eqs. (5a) and (5b) once again with the input unew , then x new (t ) ,

λ new (t ) , J [unew ] and ∂H ( k ) / ∂u( k ) are obtained. Step 6

If J [unew ] − J [u(k ) ] < ε 2 is satisfied, set

u* := unew , then take u* as the approximate solution of P1. Otherwise, set u( k ) := unew , x ( k ) (t ) := x new (t ) ,

λ ( k ) (t ) := λ new (t ) , J [u( k ) ] := J [unew ] and (∂H ( k ) / ∂u( k ) ) := (∂H new / ∂unew ) , then go back to Step 3. 4

CASE TESING

The classic batch reactor problem in chemical engineering is tested in this section. The problem is about a batch reactor operating over a one hour period and producing two products according to the parallel reaction mechanism: A → B , A → C . Both reactions are irreversible and first order in A, and have rate constants given by ⎛ E ⎞ ki = ki 0 exp ⎜ − i ⎟ i = 1, 2 (17) ⎝ RT ⎠ where k10 = 106 s−1, k20 = 5.1×1011 s−1, E1 = 41.84

kJ·mol−1, E2 = 83.69 kJ·mol−1. The objective is to find the optimal temperature operating policy below 412 K that maximizes the yield of B. The problem is formulated as: max M B (1.0) (18a) subject to

dM A = −[k1 + k2 ]M A (t ) dt dM B = k1M A (t ) dt M A (0) = M A0 M B (0) = M B0

Figure 1

The structure of improved CVI approach

(18b)

(18c)

T ≤ 412 K

(18d)

0 ≤ t ≤1

(18e)

1057

Chin. J. Chem. Eng., Vol. 20, No. 6, December 2012

(a) k1(t);

(b)

k2(t)

MA(t);

MB(t)

Figure 3 Optimal control profiles (a) and optimal state profiles (b) Table 1

this work

①

①

Results for batch reactor problem

u(0)

N

J [u* ]

BFGS iterations

CPU time/s

2

100

0.57325

21

1.11

2

200

0.57326

25

2.28

1

100

0.57330

21

1.09

3

100

0.57325

21

1.02

5

100

0.57352

50

Vassiliadis [11]

0.573530

Renfro [12]

0.57000

Logsdon [13]

0.57349

2.27 93.2

15.7

The best literature result.

The computing platform consists of Pentium E6550/2.93GHz CPU and DDRII/1000MHz system memory, both of the desired integration tolerance and optimization tolerance are set to 10−6, μ is set to 100. The number of BFGS iterations and CPU time can be used to measure the computational efficiency. In this problem N should be greater than 100 according to Eq. (16). A0 is set to 1, B0 is set to 0. The algorithm in Section 3 is applied and an optimal objective function value of 0.57325 is obtained, with 21 BFGS iterations. Figs. 3 (a) and 3 (b) show the optimal control profiles and corresponding state profiles respectively. In Fig. 3 (a), the satisfaction of the constraint Eq. (18d) is illustrated: T reaches 412 K when k1 equals to 5 and k2 equals to 12.5, these values are exactly consistent with formula Eq. (17). In Fig. 3 (b), the yield of B increases during the whole operation process, and finally reaches the maximum as the time hits 1 h. Table 1 shows the results under different conditions: when u(0) = 2 and N = 100 , the objective function value is 0.57325; when N is doubled, the objective function value is increased to 0.57326. It is revealed that as N increases, the solution turns more precise, while the computing process costs doubled CPU time. In order to observe the effect of different u(0) , the problem is resolved with u(0) = 1 , u(0) = 3 ,

and u(0) = 5 . The best objective function value of 0.57352 is obtained when u(0) = 5 , which is very close to the best literature value of 0.57353 with N = 10 by Vassiliadis (on SUN SPARC II workstation) [11], the relative error between them is less than 1.75×10−5. The uniform convergence for various initial control variables also illustrates the effectiveness of the approach. 5

CONCLUSIONS

The main contributions of this paper are to clarify an efficient indirect approach for nonlinear dynamic optimization problems without path constraints, and to apply the corresponding numerical algorithm on the benchmark case. The proposed approach can overcome most of the drawbacks in the traditional indirect approaches. REFERENCES 1 2

Biegler, L.T., Nonlinear Programming: Concepts, Algorithms, and Applications to Chemical Processes, SIAM, Philadelphia (2010). Hong, W.R., Tan, P.C., Wang, S.Q., Li, P., “A comparison of arithmetic operations for dynamic process optimization approach”, Chin. J.

1058

3

4

5

6 7 8

Chin. J. Chem. Eng., Vol. 20, No. 6, December 2012

Chem. Eng., 18 (1), 80-85 (2010). Woranee, P., Paisan, K., Amornchal, A., “Batch-to-batch optimization of batch crystallization processes”, Chin. J. Chem. Eng., 16 (1), 26-29 (2008). Zhang, B., Chen, D.Z., Wu, X.H., “Graded optimization strategy and its application to chemical dynamic optimization with fixed boundary”, Journal of Chemical Industry and Engineering (China), 56 (7), 1276-1280 (2005). (in Chinese) Gill, P.E., Murray, W., Saunders, M.A., “SNOPT: An SQP algorithm for large-scale constrained optimization”, SIAM. J. Optimiz., 12 (4), 979-1006 (2002). Bell, M.L., Sargent, R.W.H., “Optimal control of inequality constrained DAE systems”, Comput. Chem. Eng., 24 (11), 2385-2404 (2000). Kirk, D.E., Optimal Control Theory: An Introduction, Dover Publications, New York (2004). Biegler, L.T., “Solution of dynamic optimization problems by

9

10 11

12

13

successive quadratic programming and orthogonal collocation”, Comput. Chem. Eng., 8 (3-4), 243-247 (1984). Luus, R., “Further developments in the new approach to boundary condition iteration in optimal control”, Can. J. Chem. Eng., 79 (6), 968-976 (2001). Burde n, R.L., Faires, J.D., Numerical Analysis, Brooks/Cole, New York (2011). Vassiliadis, V.S., “Computational solution of dynamic optimization problems with general differential-algebraic constraints”, Ph.D. Thesis, University of London, UK (1993). Renfro, J.G., Morshedi, A.M., Asbjornsen, O.A., “Simultaneous optimization and solution of systems described by differential/algebraic equations”, Comput. Chem. Eng., 11 (5), 503-517 (1987). Logsdon, J.S., Biegler, L.T., “Decomposition strategies for large-scale dynamic optimization problems”, Chem. Eng. Sci., 47 (4), 851-864 (1992).