Optimal control of a two-stage reactor system

Optimal Control of a Hybrid Dynamical System: Two-stage Reactor System T. Hirmajer and M. Fikar Institute of Information Engineering, Automation, and Mathematics Department of Information Engineering and Process Control FCFT STU, Radlinsk´eho 9, 812 37 Bratislava, Slovakia fax: +421 2 52 49 64 69 e-mail: {tomas.hirmajer, miroslav.fikar}@stuba.sk Technical Report TH0002 http://www.kirp.chtf.stuba.sk/~hirmajer/

30.3.2006

Version 1.10

Abstract This report presents determination of optimal operation policies for a batch reactor system described by two sets of differential equations. In each stage behavior of the reactor is described by one set of differential equations. When certain conditions are satisfied the transition from the first to the second stage occurs. The maximum conversion problem is investigated subject to inequality constraints. Dynamic optimization based on control vector parameterization (CVP) is used to find optimal control profile with the regard of a general stage-transition procedure. Gradients of the resulting nonlinear programming problem (NLP) are obtained by adjoint method based on the optimal control theory. Keywords: Batch reactor, dynamic optimization, hybrid systems, control vector parameterization.

Contents 1 Introduction

4

2 A Two-Stage Reactor System 2.1 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Optimization Problem Definition . . . . . . . . . . . . . . . . . . . . . . .

6 6 8

3 General Problem Formulation

9

4 Gradient of the Cost Function

10

5 Procedure 5.1 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Implementation of Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Integration of Adjoint Equations . . . . . . . . . . . . . . . . . . . . . . . .

12 12 13 13

6 Results and Discussion 6.1 Jacobians . . . . . . . . 6.2 Discontinuities . . . . . . 6.3 Gradients Computations 6.4 Results . . . . . . . . . .

14 14 14 14 15

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

7 Conclusions

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

19

1

List of Figures 1 2 3 4 5 6

Two-stage reactor system . . . . . . . . . . . . . . . Optimal state trajectories of components A, B, C . . Optimal state trajectories of components D, E, F . . Optimal control profile for 6 discretization intervals . Optimal control profile for 10 discretization intervals Optimal control profile for 20 discretization intervals

2

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

6 17 17 18 18 19

List of Tables 1

Effect of the number of discretization intervals P on the performance index J0 16

3

1 Introduction

1

4

Introduction

Batch and semi-batch processes are of considerable importance in the chemical industry. They are widely used in the production of fine chemicals, specialties (pharmaceuticals), polymers, biotechnological, etc. Many authors have been trying to find a way how to increase productivity of batch production (Srinivasan et al. 2003, Bonvin 1998). Inherent dynamic nature is a characteristic sign of such batch process. These processes consist of the most interesting and challenging problems. Therefore, the modeling of the batch reactors and the optimization of this type of problems requires the use of dynamic optimization techniques where the total effectiveness of such process can be improved in comparison to the original operation. Systems with a switching structure are typical for batch processes, which are described with different reactions stages. For example wastewater treatment (Lin and Cheng 2001, Zhao et al. 1995, 1994), where nitrification/denitrification process takes place (Coelho et al. 2000, Isaacs 1997), optimal operation policies in batch reactors with switching structure (Barton et al. 2000, Schlegel and Marquardt 2006), optimal control of a product supply process (Manon et al. 2002), etc. The above mentioned group of processes can be viewed as systems exhibiting hybrid behavior, combining periods of continuous operation with discrete state changes. Hence many units in process industries are described by multiple sets of differential and algebraic equations. As such they are difficult to control and optimize in transient regimes if switching between the sets is to be taken into account. The switches can involve different regimes of operation (occurrence multiphase phenomena, explosive areas in mixtures of gasses, etc) or external actions (addition of second unit when production increases, etc). There are several approaches to solution of such dynamic optimization problems. If the process to be optimized can be described accurately enough by piece-wise linear and logic formulation, powerful algorithms in the area of explicit model predictive control exists (Bemporad and Morari 1999). If fully nonlinear processes are concerned, original dynamic optimization problem has to be approximated by some simplified formulation. One approach is complete discretization of state and control variables – orthogonal collocation. Such formulation can be found in Avraam et al. (1998), Biegler et al. (2002), Cuthrell and Biegler (1987), Logsdon and Biegler (1989). This approach is currently the most versatile and applicable also to mixed continuous/discrete cases (Avraam et al. 1998) but the size of resulting nonlinear programming problem (NLP) is very large. Other possibility is to leave the states intact and approximate only the control variables as piece-wise constant, or with some higher order approximations. This approach is known as control vector parameterization (CVP). Although there are some issues with CVP, especially with truly hybrid systems and unknown sequences of stage changes (Feehery 1998), in general it is considered as a very reliable and accepted method. With CVP, different formulations can be found, depending on how gradients of the resulting NLP are calculated (Rosen and Luus 1991). In Caracotsios and Stewart (1985), Vassiliadis et al. (1994), system of sensitivity equations is formed and the gradients are http://www.kirp.chtf.stuba.sk/~ hirmajer

1 Introduction

5

calculated from its solution. The advantage of this method is easy formulation of the problem and forward integration of both states and sensitivity equations. The drawback of this method is generation of a large system of differential equations as each optimized parameter corresponds to a set of differential equations with the same dimension as the number of states of the optimized process. Another possibility that is pursued in this work is to calculate the gradients of NLP via optimal control theory using the so-called co-state, or adjoint equations (Ruban 1997, Goh and Teo 1988). The advantage is that the number of differential equations is not proportional to the number of optimized parameters, but to the number of constraints. On the other side, adjoint equations have to be solved in opposite direction of time which makes the implementation more difficult. When dealing with processes comprised of a large number of state equations and only a small number of state-dependent constraints, this approach has favorable properties compared to calculation of sensitivities. The main aim of this work is to show feasibility of dynamic optimization based on adjoint CVP for a chemical reactor exhibiting hybrid dynamics. The reactor is in initial stage described by one set of differential equations and after an addition of one reactant in the middle of the operation by another set of differential equations. The chemical reactor and the dynamic optimization problem is taken from Vassiliadis et al. (1994) where it has been solved using the system of sensitivity equations. We will show attractiveness of the proposed approach and its advantages over Vassiliadis et al. (1994).

http://www.kirp.chtf.stuba.sk/~ hirmajer

2 A Two-Stage Reactor System

2

6

A Two-Stage Reactor System

The presented reactor system is based on the assumption of perfect mixing and ideal liquid mixture behavior. The kinetic model allows for the calculation of the concentration of all components every time.

2.1

Model

Figure 1: Two-stage reactor system We assume a non-linear system of two batch reactors shown in Figure 1. The first reactor is initially loaded with volume V1 of an aqueous solution of component A of concentration x1 (t0 ) and with a solid catalyst. The heating coil is used as a control variable for the first stage – manipulates the reactor temperature over the time. When the first reaction chain 2A → B → C

(1)

is finished at time tn , an amount of dilute aqueous solution of component B of concentration xs2 is added to the product of the first reactor stage and loaded to the second reactor where the three parallel reactions B → D,

B → E,

2B → F

(2)

take place with fixed temperature – isothermal conditions (second stage). Detailed information about this process is described in Vassiliadis et al. (1994). Presented two-stage reactor system is described by two sets of differential equations f1

http://www.kirp.chtf.stuba.sk/~ hirmajer

2.1

Model

7

and f2 :  dx  1

 dt   dx2      −2k1 (u)x21  dt   dx  k1 (u)x21 − k2 (u)x2   3       k (u)x 2 2 dt   , f1 =   dx4  =   0      dt    0  dx  5   0    dt  dx6 dt

(3)

 dx  1

 dt   dx2      0  dt   dx  −0.02x2 − 0.05x2 − 0.00008x22    3      0 , dt   f2 =    dx4  =  0.02x2      dt    0.05x 2  dx  2 5   0.00004x2    dt  dx6 dt

(4)

with kinetic constants defined as k1 (u) = 0.0444e−2500/u ,

(5)

k2 (u) = 6889.0e−5000/u ,

(6)

where x1−6 [mol m−3 ] – state values – concentration of the components; u [K] – control value – profile of the reactor temperature for the first stage; k1 (u) [m3 mol−1 min−1 ], k2 (u) [min−1 ] – rate constants. The mixing operation at the switching time tn is described by following equations: − V2 x1 (t+ n ) = V1 x1 (tn ), − s V2 x2 (t+ n ) = V1 x2 (tn ) + Sx2 , − V2 x3 (t+ n ) = V1 x3 (tn ),

(7) (8) (9)

where V1 [m3 ] – volume of the material loaded into the first reactor; S [m3 ] – amount of the solution of the component B added after first stage with concentration xs2 [mol m−3 ]; V2 [m3 ] – volume of the material loaded in the second reactor.

http://www.kirp.chtf.stuba.sk/~ hirmajer

2.2

2.2

Optimization Problem Definition

8

Optimization Problem Definition

The aim of the optimization is to reach a maximum concentration of the amount of component D not later than at time tP subject to a constraint that the final concentration of D should be greater than some desired values xw 4. The decision variables are the temperature profile of the first reaction stage, the durations of the two stages, and the amount S of the component B added at tn . Thus, the cost functional is defined as max

S,∆ti ,u[0,tn ]

J0 = V2 x4 (tP ),

(10)

subject to the constraints J1 = x4 (tP ) − xw 4 ≥ 0, J2 = tP −

P X

∆ti ≥ 0,

(11) (12)

i=1

−3 where xw 4 [mol m ] – required minimum value of the concentration of the component D at the final time, tP [min] – processing time for the both reaction stages. The concrete values of the parameters are defined as: volume V1 = 0.1 m3 ; initial states x1 (t0 ) = 2000 mol m−3 , x2−6 (t0 ) = 0 mol m−3 ; concentration of added component −3 B xs2 = 600 mol m−3 ; desired minimum state xw 4 = 150 mol m ; maximum final time tP = 180 min.

http://www.kirp.chtf.stuba.sk/~ hirmajer

3 General Problem Formulation

3

9

General Problem Formulation

Presented process of the optimization which is described by sets of the differential equations (3), (4) can be written in a general form x˙ = f i (t, x(t), u), ti−1 ≤ t ≤ ti , i = 1, n, x˙ = f n+1 (t, x(t), u), tn ≤ t ≤ tP ,

(13) (14)

with initial condition x(t0 ) = x0 (t0 , u). Next, we consider here the nx – dimensional process state vector x(t) and a constant nu – dimensional vector of optimized parameters u. Note that u includes not only control trajectory, but also other time independent parameters p. It is assumed that system states are governed by different state equations f i that are continuously differentiable. The switching instant ti (u) that can be a function of u, at which one set of equations (3) is replaced by another (4) is determined by the switching condition g i (ti , x(t− i ), u) = 0,

i = 1, n

(15)

We assume that the functions g i are continuously differentiable with respect to all variables. At switching instants, the vector x(t) can have breaks defined by the equations (8) − − x(t+ i ) = x(ti ) + ∆i (ti , x(ti ), u),

i = 1, n,

(16)

+ where ∆i are also continuously differentiable vector functions and x(t− i ), x(ti ) are the values of the vector x(t) before and after the switching instant, respectively. In (16) we consider additive jumps that are superimposed on the continuous trajectory at points ti (e.g. if ∆i = 0, then x(t) is continuous at the point ti ). In the next step we define the cost J (10) which is to be optimized in a general Bolza form (Bryson and Ho 1975)

J(u) = G(tP , x(tP ), u) +

ZtP

F (t, x(t), u)dt,

(17)

t0

where J is the scalar performance index to be minimized (functional of the quality of the dynamic system), G defines final time conditions and F requirements along the time axis. In a similar manner we define constraints to be satisfied (11), (12).

http://www.kirp.chtf.stuba.sk/~ hirmajer

4 Gradient of the Cost Function

4

10

Gradient of the Cost Function

As the original dynamic optimization problem will be transformed into a static one – NLP, we need to derive gradients of the cost function and constraints with respect to optimized parameters – u. These will be obtained from the optimal control theory by the variational method. The gradients computations consider to the optimized parameters – p are similar. The detailed derivations of the gradients are beyond the scope of this paper, the reader is referred to Hirmajer and Fikar (2005). Optimality conditions for our problem differ slightly from the original approach by Pontryagin et al. (1962) and have been derived by Ruban (1997). First, the Hamiltonian H is defined as H(t, x, u, λ) = F + λT f ,

(18)

where the vector of the Lagrange multipliers λ is defined by the differential equation ∂H ˙ λ(t) =− , ∂x

(19)

with the boundary value at the point tP and the switching condition of the Lagrange multipliers at switching instants ∂G , ∂x(t )  P h λT (t− λT (t+ i ) = i ) I +

(20)

λ(tP ) =

 i  ∂∆ ∂∆i i + − f i+1 (ti ) ai + ∂ti ∂xT (t− i ) 
−1 − + + (F (ti ) − F (ti ))ai I − f i (t− )a i i

(21)

The gradient of the cost function is then defined as   t− i   Z P ∂G ∂G ∂H X T + ∂x0 (t0 , u0 )  − + H(tP ) δtP δJ =  dt + + λ (t0 )  δui + ∂tP ∂uTi ∂uTi ∂uTi i=1 t+ i−1

+

n X i=1

"

λT (t+ i )



∂∆i ∂uTi



#  ∂∆ i h ∂G i + + λT (t+ bi δti , δui + H(t− i ) i ) − H(ti ) + ∂ti ∂ti (22)

where the coefficients ai , bi are defined as  −1   ∂g i (ti , x(t− ∂g i (ti , x(t− i ), u) i ), u) ai = − , ∂ti ∂xT (t− i ) bi = −



∂g i (ti , x(t− i ), u) ∂ti

−1 

∂g i (ti , x(t− i ), u) T ∂u

http://www.kirp.chtf.stuba.sk/~ hirmajer



(23)

(24)

4 Gradient of the Cost Function

11

If ai = 0 the switching conditions are simplified considerably (if the value of the jump of state variables of the object ∆i does not depend on x(t− i ), then the multipliers are continuous): λ

T

(t− i )

=λ

T

(t+ i )

h  ∂∆ i i I+ , T ∂x (t− i )

i = 1, n

(25)

If state variables of the process are continuous at switching points (∆i = 0), then the integrand function F of the quality index is continuous, but the Lagrange multipliers contain discontinuity at the moment of switch  
−1 T + + − λT (t− , i = 1, n (26) i ) = λ (ti ) I − f i+1 (ti )ai I − f i (ti )ai Finally, if the times of switch depend explicitly on control parameter (ti = ti (u), ai = 0, i ), then the Lagrange multipliers are continuous. and bi = dt du

http://www.kirp.chtf.stuba.sk/~ hirmajer

5 Procedure

5

12

Procedure

5.1

Algorithm

We assume that the continuous control trajectory is considered to be piece-wise constant over P intervals. This makes it possible to convert the original problem of dynamic optimization into nonlinear programming problem (Fikar and Latifi 2001). In this algorithm we assume that we have the functional J0 and k constraints Jj , where j = 1, k. We further separate the optimized variables into times ti , constant controls u and parameters p. For simplicity assume that the initial time t0 and state x(t0 ) are given and constant. Then we can write the following algorithm: 1. Integrate the system (13) and integral terms Fj together from t = t0 to t = tP . Restart integration with switching conditions (15), states can be discontinuous following the equation (16). 2. For j = 0, k repeat (a) Initialize adjoint variables λj (tP ), according to equations (20). (b) Initialize the intermediate variables J D,j , J P,j as zero. These represent integral part of the gradients. (c) Integrate backwards from t = tP to t = t0 the adjoint system (19) and intermediate variables. Allow for discontinuities of the adjoint equations as given in (25), restart integration at these points, and at the points of changes of dynamics ∂Hj , λ˙ j = − ∂xi ∂Hj , J˙D,j = ∂ui ∂Hj J˙P,j = ∂p

(27) (28) (29)

(d) Calculate the gradients of Jj with respect to times ti , controls u and parameters p, with help of (22). ∂Jj ∂tP ∂Jj ∂ti ∂Jj ∂ui ∂Jj ∂pT

∂Gj = Hj (t− , P) + ∂tP   ∂∆  i T + − + = Hj (ti ) − Hj (ti ) + λj (ti ) bi , ∂ti = J D,j (ti−1 ) − J D,j (ti ), =

J TP,j (t0 )

+λ

T

i = 1, n,

(31)

i = 1, P ,

∂x0 (t0 , u0 ) (t+ 0) T ∂p

http://www.kirp.chtf.stuba.sk/~ hirmajer

(30)

(32) n

∂Gj X T + + T + λ (ti ) ∂p i=1



∂∆i ∂pT



(33)

5.2

Implementation of Algorithm

13

In this manner, the values of Jj are obtained in the step 1 and the values of gradients in the step 2d. This is all what is needed in NLP algorithm. For numerical reasons, time increments ∆ti will be optimized, rather than absolute time values ti . Therefore, the gradients have to be modified correspondingly. The relations between times and their increments are given as tP =

P X

∆ti

(34)

i=1

Therefore, the following holds for the derivatives P

X ∂Jj ∂tr ∂Jj = ∂∆ti ∂tr ∂∆ti r=1

5.2

(35)

Implementation of Algorithm

The algorithm was implemented FORTRAN 77. The first part of this program contains a module for forward and backward integrations LSODAR (Petzold and Hindmarsh 1997) that is able to handle state events. The second module of the program computes gradients and calls the NLP solver NLPQL (Schittkowski 1981). The user specifies initial conditions: • u0 , x0 , lower and upper control bounds, number of time intervals and subintervals • cost function, constraints, • differential equations of the process.

5.3

Integration of Adjoint Equations

When the adjoint equations are integrated backwards in time, the knowledge of states x(t) is needed. In our case the program stores at first in the forward pass the states at a predefined grid points and interpolates them when adjoint equations are solved. Two interpolations were implemented into the program: linear, and the approximations having continuous states and continuous first order derivatives across boundaries. The second one was used in simulation.

http://www.kirp.chtf.stuba.sk/~ hirmajer

6 Results and Discussion

6 6.1

14

Results and Discussion Jacobians

The Jacobians of the reactor system with of the form  −4k1 (u)x1 0 0 0  2k1 (u)x1 −k2 (u) 0 0  ∂f 1  0 k2 (u) 0 0 =  0 0 0 0 ∂x   0 0 0 0 0 0 0 0

respect to the states, control and parameter are 0 0 0 0 0 0

 0 0  0 , 0  0 0

 ) −2k1 (u)x21 ( 2500 u2 k1 (u)x21 ( 2500 ) − k2 (u)x2 ( 5000 ) u2 u2   5000  ∂f 1  k2 (u)x2 ( u2 ) , =   0 ∂u     0 0

(36)



 0 0 0 −0.07 − 0.00016x2  ∂f 2  0 0 =  0.02 ∂x 0 0 0.05 0 0.00008x2

and other partial derivatives

6.2

∂f 1 ∂f 2 , ∂pT ∂uT

0 0 0 0 0 0 ,

0 0 0 0 0 0

∂f 2 ∂pT

0 0 0 0 0 0

(37)

 0 0  0  0  0 0

(38)

are equal to 0.

Discontinuities

We consider at the switching time discontinuities ∆i following the equation (16). These discontinuities are described according to the adjusted equations (7),(8),(9):   −x1 (t− ) n S  S  x2 (tn ) − x2 (t− ∆n (tn , x(t− (39) n) n ), u) = V1 + S − −x3 (tn )

6.3

Gradients Computations

Detailed information about the gradients of the process are needed in each NLP iteration. These were derived from the optimal control theory. Detailed gradients formula of the cost and constrains are presented as follows. http://www.kirp.chtf.stuba.sk/~ hirmajer

6.4

Results

15

Cost J0 :
λT0 (tP ) = 0, 0, 0, V2, 0, 0 , ∂J0 − + T + = λT0 (t− i )f (ti ) − λ0 (ti )f (ti ), ∂ti ∂J0 − = λT0 (t− P )f 2 (tP ), ∂tP ∂J0 = J D,0 (ti−1 ) − J D,0 (ti ), ∂ui   ∂J0 ∂G0 ∂∆i T + = + λ (ti ) ∂pT ∂pT ∂pT

(40) (41) (42) (43) (44)

Constraint (J1 ):
λT1 (tP ) = 0, 0, 0, 1, 0, 0 , ∂J1 − T + + = λT1 (t− i )f (ti ) − λ1 (ti )f (ti ), ∂ti ∂J1 − = λT1 (t− P )f 2 (tP ), ∂tP ∂J1 = J D,1 (ti−1 ) − J D,1 (ti ), ∂ui   ∂J1 ∂∆i T + = λ (ti ) ∂pT ∂pT

(45) (46) (47) (48) (49)

Constraint (J2 ): λ2 (tP ) = 0,

6.4

∂J2 = 1, ∂ti

∂J2 = 1, ∂tP

∂J2 ∂J2 = =0 ∂ui ∂pT

(50)

Results

We have considered the following scenario. The first stage contains n time intervals and the second stage only one interval. Therefore we can write the whole number of the intervals: P = n + 1. Several simulations were performed with different number of intervals displayed in Table 1. The initial values of the optimized parameters are: ∆ti = 15 min, ui = 350 K, S = 0.1 m3 with bounds defined as ui ∈ [298, 398] for t ∈ [t0 , tn ], S ∈ [0, 0.1] and ∆ti ∈ [10, 100] for P = 6 otherwise ∆ti ∈ [1, 100]. Optimization and integration tolerances were set to 10−4 and 10−10 , respectively. The optimal concentrations of components A, B, C can be found in Figure 2 and concentrations of components D, E, F in Figure 3 with corresponding control trajectory shown in Figure 4. The first and the second constraint ensures that the final value of the component D is 150 mol m−3 and the optimal duration of the process is 180 min, respectively. The time of http://www.kirp.chtf.stuba.sk/~ hirmajer

6.4

Results

Number of Intervals, P 6 10 20

16 Number of Time of Iterations Switch, tn [min] 39 106.04 51 104.98 64 106.82

Amount of Component B, S [m3 ] 0.0702 0.0705 0.0705

Performance Index, J0 [mol] 25.5365 25.5681 25.5755

Table 1: Effect of the number of discretization intervals P on the performance index J0 the switch and the corresponding injected volume of the component B after the first stage are shown in Table 1. Finally, the optimal value of the cost function is 25.54 mol for the first instance that is slightly less than the one in Vassiliadis et al. (1994), where optimal value of the cost 25.55 mol was found. This can be caused by integration and optimization tolerances. This example illustrates the fact that the adjoint approach can have computational advantages over the sensitivity equations approach. Indeed, this is clear, especially when a large number of variables is optimized with a relative small number of constraints, and the system is described by a large number of equations. As we can see in the example for 6 discretization intervals the number of differential equations to be integrated was 27 times per one NLP iteration (6 differential equations forward and 3 times 7 equations backward (6 adjoint equations, 1 integral part of Hamiltonian) every constraint and cost). The sensitivity equations approach needs to integrate 72 differential equations (number of states multiplied by the number of optimized variables). Thus, the adjoint approach is faster and the difference is larger clear if either the number of states or number of optimized parameters increases. Figures 4, 5, and 6 show optimal control trajectory with increasing number of time intervals. This smoothens the control considerably. However, due to the adjoint approach to CVP, the number of integrations in gradients has not changed as it would using the sensitivity approach. Although the cost function changes only a little, this could result in large profit increase if the desired product would be expensive.

http://www.kirp.chtf.stuba.sk/~ hirmajer

6.4

Results

17

2000 cA

1750

cB cC

c [mol/m3]

1500 1250 1000 750 500 250 0 0

30

60

90 120 time [min]

150

180

Figure 2: Optimal state trajectories of components A, B, C

400 cD

350

cE cF

c [mol/m3]

300 250 200 150 100 50 0 0

30

60

90 time [min]

120

150

180

Figure 3: Optimal state trajectories of components D, E, F

http://www.kirp.chtf.stuba.sk/~ hirmajer

6.4

Results

18

345 340 335

T [K]

330 325 320 315 310 305

0

30

60

90 time [min]

120

150

180

Figure 4: Optimal control profile for 6 discretization intervals

345 340 335

T [K]

330 325 320 315 310 305

0

30

60

90 time [min]

120

150

180

Figure 5: Optimal control profile for 10 discretization intervals

http://www.kirp.chtf.stuba.sk/~ hirmajer

7 Conclusions

19

345 340 335

T [K]

330 325 320 315 310 305

0

30

60

90 time [min]

120

150

180

Figure 6: Optimal control profile for 20 discretization intervals

7

Conclusions

This report has considered the formulation and solution of dynamic optimization problem for a two-stage reactor system described by general stage-transition procedure. Each stage was described by sets of differential equations. Control vector parameterization has been used and gradients for the nonlinear programming have been calculated based on the optimal control theory. This is in contrast to usual approaches where sensitivity equations are preferred due to simplicity of the implementation. On the other hand, the adjoint variable approach has its advantage in systems described by a larger set of differential equations and it can reduce the computational time considerably because evaluate of the gradients is usually much faster than sensitivity approach. Simulation with a two-stage reactor system first presented in Vassiliadis et al. (1994) confirmed the ability of this method to solve a complex problem. As we can see in example with 6 discretization intervals the number of differential equations to be integrated is 27 per one NLP iteration opposite to the sensitivity equations approach where 72 differential equations are to be considered. This difference in the number of integrations between adjoint approach and sensitivity equations rises further up with the increasing number of intervals. It follows that the computational time can be considerably smaller compared to the sensitivity equations.

Acknowledgments The authors are pleased to acknowledge the financial support of the Scientific Grant Agency of the Slovak Republic under grants No. 1/1046/04 and 1/3081/06.

REFERENCES

20

References M. P. Avraam, N. Shah, and C. C. Pantelides. Modelling and optimisation of general hybrid systems in the continuous time domain. Computers and Chemical Engineering, 22(Suppl.):S221–S228, 1998. 4 P.I. Barton, S. Gal´an, and J.R. Banga. Optimization of Hybrid Discrete/Continuous Dynamic Systems. Computers and Chemical Engineering, 24(4):851–864, 2000. 4 A. Bemporad and M. Morari. Control of systems integrating logic, dynamics, and constraints. Automatica, 35(3):407 – 427, 1999. 4 L. T. Biegler, A. M. Cervantes, and A. W¨achter. Advances in simultaneous strategies for dynamic process optimization. Chemical Engineering Science, 57(4):575–593, 2002. 4 D. Bonvin. Optimal operation of batch reactors - a personal view. Journal of Process Control, 8(5-6):355–368, 1998. 4 A. E. Bryson and Yu-Chi Ho. Applied Optimal Control - Optimization, Estimation and Control. Hemisphere publishing corporation, Taylor & Francis, Inc., 1975. 9 M. Caracotsios and W. E. Stewart. Sensitivity analysis of initial value problems with mixed ODEs and algebraic equations. Computers and Chemical Engineering, 9:359–365, 1985. 4 M. A. Z. Coelho, C. Russo, and O. Q. F. Arajo. Optimization of a sequencing batch reactor for biological nitrogen removal. Water Research, 34(10):2809–2817, 2000. 4 J. E. Cuthrell and L. T. Biegler. On the optimization of differential-algebraic process systems. AIChE J., 33(8):1257–1270, 1987. 4 W. F. Feehery. Dynamic Optimization with Path Constraints. PhD thesis, MIT, 1998. 4 M. Fikar and M. A. Latifi. User’s guide for FORTRAN dynamic optimisation code DYNO. Technical Report mf0201, LSGC CNRS, Nancy, France; STU Bratislava, Slovak Republic, 2001. 12 C. J. Goh and K. L. Teo. Control parameterization : a unified approach to optimal control problems with general constraints. Automatica, 24(1):3–18, 1988. 5 T. Hirmajer and M. Fikar. Optimal control of a hybrid dynamical system: Two coupled tanks. Technical report, Department of Information Engineering and Process Control, FCFT STU, Radlinskho 9, 812 37 Bratislava, Slovakia, 2005. 10 S. Isaacs. Automatic adjustment of cycle length and aeration time for improved nitrogen removal in an alternating activated sludge process. 35(1):225–232, 1997. 4

REFERENCES

21

S. H. Lin and K. W. Cheng. A new sequencing batch reactor for treatment of municipal sewage wastewater for agricultural reuse. Desalination, 133(41-51), 2001. 4 J. S. Logsdon and L. T. Biegler. Accurate solution of differential-algebraic optimization problem. IEC Res., 28:1628–1639, 1989. 4 P. Manon, C. V. Roubinet, and G. Gilles. Optimal control of hybrid dynamical systems:application in process engineering. Control Engineering Practice, 10:133–149, 2002. 4 L. R. Petzold and A. C. Hindmarsh. LSODAR. Technical report, Computing and mathematics research division, lawrence livermore national laboratory livermore, ca 94550, 1997. 13 L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, and E. F. Mishchenko. The Mathematical Theory of Optimal Processes. Wiley, New York, 1962. 10 O. Rosen and R. Luus. Evaluation of gradients for piecewise constant optimal control. Computers and Chemical Engineering, 15(4):273–281, 1991. 4 A. I. Ruban. Sensitivity coefficients for discontinuous dynamic systems. Computer and Systems Sciences International, 36:536–542, 1997. 5, 10 K. Schittkowski. NLPQL: A new fortran implementation of a sequential quadratic programming algorithm for parallel computing. Technical report, Department of Mathematics, University of Bayreuth, D - 95440 Bayreuth, 1981. 13 M. Schlegel and W. Marquardt. Detection and exploitation of the control switching structure in the solution of dynamic optimization problems. Journal of Process Control, 16: 275–290, 2006. 4 B. Srinivasan, S. Palanki, and D. Bonvin. Dynamic optimization of batch processes I. Characterization of the nominal solution. Computers and Chemical Engineering, 27: 1–26, 2003. 4 V. S. Vassiliadis, R. W.H. Sargent, and C. C. Pantelides. Solution of a class of multistage dynamic optimization problems. 1. Problems without path constraints, 2. Problems with path constraints. IEC Res., 33:2111–2122, 2123–2133, 1994. 4, 5, 6, 16, 19 H. Zhao, S. H. Isaacs, A. Søeberg, and M. K¨ ummel. A novel control strategy for improved nitrogen removal in an alternating activated sludge process - Part II. Control development. Water Research, 28(3):535–542, 1994. 4 H. Zhao, S. H. Isaacs, H. Søeberg, and M. K¨ ummel. An analysis of nitrogen removal and control strategies in an alternating activated sludge process. Water Research, 29(2): 535–544, 1995. 4

View publication stats