Optimal Control of Distributed Processes Using Restricted Second Order Information

Copyright © IFAC Advanced Control of Chemical Processes, Pisa. Italy, 2000

OPTIMAL CONTROL OF DISTRIBUTED PROCESSES USING RESTRICTED SECOND ORDER INFORMATION. Eva Balsa-Canto' , Julio R. Banga' Antonio A. Alonso" and Vassilios S. Vassiliadis •••

• Chemical Engineering Lab. Instituto de Investigaciones Marinas (CSIC) , Vigo, Spain. e-mail: [email protected]. julio @iim.csic.es •• Department of Chemical Engineering, Universidad de Vigo, Vigo , Spain. e-mail:[email protected] ••• Department of Chemical Engineer'ing, Cambridge Univer'sity, Cambridge, UK. e-mail:[email protected]

Abstract: The open loop optimal control of combined lumped and distributed process systems is considered here. A method based on the Control Vector Parameterization approach, which uses exact second order information, is presented for the efficient computation of the control policy that optimizes a specific performance criterion. This approach makes use of first and second order sensitivities to obtain exact gradients and restricted second order information by the solution of an extended Initial Value Problem (IVP) . We illustrate the use of this method considering the solution of several case studies of increasing complexity. All these case studies were successfully solved with a modest computational effort, proving its capability to solve large scale, challenging problems in an efficient way. Copyright © 2000 IFAC Keywords: Distributed parameter systems, optimal control, sensitivity analysis, process control , optimization problems.

1. INTRODUCTION.

Find u(t) to minimize: J[x , y, u , v , ~ , tJ

Many chemical and biochemical processes must be modeled using a combination of both lumped and distributed parameter systems. The resulting mathematical models consist of sets of Partial Differential Equations (PDEs) and/or Ordinary Differential and Algebraic Equations (DAEs) : expressing the physical laws of conservation of mass, energy and momentum.

(1)

Subject to: lJI(x , x~ , x~~ , Xt , y , y ,

.. ., v,~,t) = 0

(2)

with ~ E n c R3 the space variables, x(~ , t) E X c RV the state variables depending on the position and time, y(t) EYe R the subset of state variables which are position invariant, x~ = 8x/8( x~~ = 8 2 x/8e , Xt = 8x/8t , Y = dy / dt , u E U c R" the control variables and v EVe RP the time invariant parameters. In order to solve the set of PDAE (2) , appropriate initial and boundary conditions must be imposed.

The efficient optimal control of such class of processes is considered in this work. A recently developed deterministic method is presented as an efficient and robust alternative. This method is based on the Control Vector Parameterization (CVP) approach and makes use of first and second order sensitivities to compute exact gradients and restricted second order information (as Hessianvector dot product). The capabilities of this new method are shown in this work considering the solution of several case studies of increasing complexity.

Path and point constraints may be imposed to the system to express conditions related to safe or proper operation: g(x(~ , t), y(t),

u(t), v,~, t) :::; 0

(3)

The control variables u and the parameters v may be subject to upper and lower bounds: uL(t) :::; u(t) :::; uU(t) and v L :::; v:::; vU (4)

The general dynamic optimization (open loop optimal control) problem considered here can be stated as follows: 881

2. SOLUTION APPROACH. Basically, existing methods for the solution of Optimal Control Problems (OCPs) can be classified in two classes (see review in VassiJiadis, 1993): indirect and direct approaches. Indirect approaches intend to solve the original problem through the application of the maximum principle of Pontryagin, either the distributed version or the discrete distributed version (Neittaanmaki and Tiba, 1994), for the problems considered here. These methods result in a two or multi-point boundary value problem which is usually very difficult to solve. Direct approaches transform the original problem into a Nonlinear Programming Problem (NLP), either using control parameterization (VassiJiadis, 1993) or complete (control and states) parameterization (Cuthrell and Biegler, 1989). The resulting NLP may be solved using a standard optimization method. In this work the CVP method will be used due to its ability to handle large dynamic optimization problems without solving very large NLPs. It proceeds dividing the duration of the process into a number of control intervals and approximating the control function using a low order polynomial. E.g. the controls may be approximated by some type of discretization which involves some (or all) of the parameters in set v, U = u(v) . Here piecewise constant approximations will be used, so the decision variables in the optimization process are the control levels. The solution of the master NLP involves the solution of the initial value problem (IVP) (system of PDAES, in our case), for each function evaluation.

This property allows on-line implementations in a model predictive control (MPC) framework (Sanchez et al., 1998) . After transformation or discretization of the original PDAEs, Eqn.2, the resulting set of DAEs describing the system can be written as follows:

= On,t E [to , t,]

f(z,z , u , v)

(5)

with the following set of initial conditions:

z(to)

= zo(v)

(6)

where: z, z E Rn are the new state (output) variables and their time derivatives respectively, u E R" are control (input) variables, and v E RP are time-invariant parameters. The functions fare such that: f : Rn X Rn X R" X RP H Rn. In Vassiliadis (1993) , the use of the first order sensitivity equations within the CVP approach was proposed to obtain gradients with respect to the optimization parameters. These can be obtained by a chain rule differentiation applied to the system in equations (5) with respect to the time-invariant parameters: 8r 8z + 8f 8z + 8f OU + or 8z OV 8z OV OU ov ov

= Onxp

(7)

with the initial conditions: 8z (to) 8v

= 8zo (v) 8v

(8)

The system Jacobian in Eqn. 7 is the partitioned matrix:

or]

Of [ oz ' 8z

Several approaches for the solution of PDAEs have been proposed in the literature (see reviews in Schiesser, 1994 and Oh, 1995). Most of them apply standard spatial discretization techniques to transform the original PDEs into ODEs or DAEs. The finite difference, finite elements and finite volume approaches proceed dividing the domain of interest into many smaller sub domains, and applying an approximation method within each subdomain. This approach has been followed for the optimization by several authors (Petzold et al., 1996; Blatt and Schittkowski, 1997; Banga et al.,1994).

(9)

which is the same as the one in the resulting equations used in the integration of Eqn. (5) , as noted in Leis and Kramer (1985) . A second differentiation (as presented by VassiJiadis et al ., 1999) of Eqn. (7) and (8) with respect to v yields the following equations: {[ a! a70 @

{[

Alternative approaches involve projection of the PDEs on functions globally defined on the spatial domain. These may correspond to the eigenfunctions of the spatial operator or more general sets such as those employed in the so-called proper orthogonal decomposition techniques (Holmes et al., 1996). The main advantage of these procedures is the reduction on computation time in objective function evaluation.

I

I PJ {jVI a2,;, + [af az

z} e IJa P {jVI + 2

a:l)T] [aWfa:l e ( av Qv 2

n

a2 f av az + ouo., a2 f au+ a2 f] + azaz ov 8va,;,

+ [I n

e

)T] [ az8,;, a2 f av a';; + "lfi.'I a2 f ov az + OUOZ a2 f au a2 f ] (az 8v ov + 8V8z

+[af au

n.

1PJ {jVI a u+

2

'
[I

n

e (au)T] av

(10)

with initial conditions given by:

02 Z 8v 2 (to) 882

=

82 z0

ov 2 (v)

(11)

where ~ is the Kronecker matrix-product operator, and the second order derivatives result in tensors of third order. To clarify the similarity between the first and second order sensitivities: we present the product of a time-invariant vector p E RP with the Hessians for each state. The result of Eqn. (10) is post-multiplied by p and comparing terms the equivalent form is derived:

8f . T 8z s

8f T A(' + 8zs + Z,Z,tl:V ) =

0n.p " l

(

12

)

where s, s E RP x n are the matrices whose n columns are respectively given by:

(2) Derives symbolically the set of first order sensitivities and the exact product equations, as presented in Eqns. (7) and (12). (3) Generates the Jacobian as given in (9) and its sparse structure in order to guarantee efficiency in the integration. (4) Transforms the resulting equations into a suitable format ready to be used by a standard solver. Note that a proper solver must be used. In particular, a method which allows the solution of sparse and possibly stiff systems , as for example: LSODES (Hindmarsh , 1983) . In this paper we present results obtained using the numerical method of lines (NMOL) (Schiesser, 1991) in step (1) above. 3. CASE STUDIES.

3.1 Temperature control oJ a slab .

with their correspondent initial conditions: Si(tO)

= ( pT [ 82z0 8v 2 '' (v) ])

We first consider the simple problem of adjusting the temperature distribution of a slab of length L=l , to a desired temperature T set = 0.2 using the external temperature (To tl as the control variable. The system is described in terms of the heat conduction equation subject to boundary conditions expressing the flux of heat to the slab. The statement is (Ramirez, 1994):

(14)

The Jacobian for the equation (12) is the same as in equation (9). This allows efficient integration of Eqns. (7) and (12) in a decoupled direct way, simultaneously with that of Eqn. (5). An interesting feature of the derivation in equation (12) is that for Truncated Newton type algorithms (e.g. Nash, 1984): only the product of the Hessian matrix with a search direction vector is required being very efficient in the optimization of large scale problems. The TN approach uses an iterative scheme (conjugate-gradient based) to solve the Newton equations of the optimization problem: H(x) p = -g(x), where H(x) is the Hessian matrix, p is the search direction and g(x) is the gradient vector. An iterative solution of the set of these equations requires the product , = H(x) p. In the case of the dynamic optimization problems considered here, this product can be obtained exactly with comparable computational cost to that of the first order sensitivity evaluation using Eqn. (12) . This approach was tested with the TN (Truncated Newton) optimization solver of Nash (1984) after the necessary modifications in order to introduce these exact products.

Find To t(t) to minimize: L

]=

J

(Tset - T(x, tf ))2 dx

(15)

o subject to: (16) with the initial and boundary conditions: T(x,O)

Tx(O , t)

= K(T -

To

=0

tl, Tx(L, t)

=0

(17)

Results are reported in Table 1, where p is the control discretization level, OB] is the objective function value, T max and T min are the maximum and minimum temperature values at t f , NI/ N F / CC are the iteration, function evaluation and H.p evaluation numbers respectively, and Tcpu is the computation time, in seconds, on a PC Pentium III/450 MHz.

In order to fully automate the application of this procedure, a Mathematica™ 3.0 notebook has been developed to derive symbolically the extended initial value problem (IVP). It proceeds as follows :

Table 1 p

10 80 160 320

(1) Derives symbolically the set of DAEs, from the original PDAE (Eqn. 2), through the application of the discretization or transformation procedure.

883

OB] 6.5 7.9 6.6 4.4

x x x x

10 5 10- 5 10- 5 10- 5

Tmax 0.203 0.202 0.202 0.203

T min NI/NF/CG Tcpu 1.3 0.196 3/4/7 29 0.195 2/3/8 87 0.195 2/3/8 767 0.196 4/5/13

The results achieved with this new methodology are in good agreement with the ones reported by Ramirez (1994) . The optimal control profile obtained for the discretization level p = 320 is shown in Figure 1. 0 .50

UA5

This method is based on the use of the evp approach and a Sequential Quadratic Programming method for the solution of the master NLP. Results obtained with the new method compare satisfactorily with the ones obtained with this standard SQP based strategy. However the reduction in number of iterations resulted in speed ups of almost 5 times with the new approach.

U.-IO

Table 2. Results for M = 100IN and (P=O.OI)

IUS 0 .)0

~

.....

O.2S

OBJ 4.00 x 10 4 3.84 x 104 3.78 x 104 3.73 x 104

p

U .2U

IUS

5 10 20 40

0 . 10 0 .05

0 .00 0 .0

0.2 TUn€'

0.1

(U

0 .4

Fig. 1. Optimal control profile for the heating of a slab (p = 320).

Tmax

T min

100.532 100.841 100.654 100.640

99.037 98.760 98.904 98.948

NI/NF/CG Tc pu 22 3/4/8 47 3/4/9 164 3/4/9 835 3/4/9

ISO

1l'0

140

G ~ :;

3.2 Temperature control of a cylinder.

Cl

.....

I~()

100 SO

This system is described by the heat conduction equation in 2 dimensions (center plane of the cylinder due to symmetry) . The objective is to maximize the uniformity on the final temperature around a desired value {Tsed restricting the movements of the external temperature (control variable).

60 40 ()

250

500

7501000 Time Is)

1250

1500

Fig. 2. Optimal control profile for the heating of a cylinder (p = 40).

Find Tout (t) to minimize:

'0'

subject to:

....

(19)

2

where Z is the vector of differences between the temperatures inside the cylinder and T set , M is a weight matrix, rand z are the radial and axial coordinates respectively and P is weight factor which restricts the control variable movement. The initial and boundary conditions are:

T(r, z, 0) = 25.0 Tr{O, z, t) = Tz(r, 0, t) = 0.0 Tr{R , z, t) h{Tout - T{R, z, t)) Tz{r , L, t) = h(Tout - T{r, L, t))

=

° ° Fig. 3. Final temperature profile in the center plane of the cylinder.

(20)

3.3 Diffusion througn. a membrane. This example is taken from Neittaanmaki and Tiba (1994). A membrane M separates two compartments I and n. M consists of an inactive protein coreticulated with an enzyme E. At time t=O the membrane does not contain any substrate or product. Then , by diffusion, some substrate S comes from I and II intro the membrane, where the enzyme E is a catalyst of the reaction to produce a certain product P.

Results obtained for Tset = 100 0 e are reported in Table 2 and the optimal control profile and final temperature distribution are shown in Figures 2 and 3 respectively. For the sake of comparison, a standard dynamic optimization method , was also considered for the solution of this case study. 884

Assume that the concentration of an activator (a) in I and II is at our disposal. The objective is to regulate the concentration of substrate (Cl) at a certain point (x=0.5), that is, to make it close to a desired concentration Zd = a x arctan(30.0t) . The statement is: Find al (t) and a2(t) to minimize:

=

J

(Cl (0.5, t) - Zd(t))2dt

(21)

(24)

Results obtained with the TN approach are presented in table 3 and the optimal control profile and the final concentration are shown in figures 4 and 5 respectively. Table 3 p

OBJ

NI/NF/CG

10 15 25

0.94079 0.93660 0.93307

10/15/22 8/11/17 12/15/24

Tcpu 38.8 48.5 130.0

........... .Actiuator 1

45 .0

30 .0

!"

20 .0

'E

flc:

"

u

_.ActhQtot 1

25 .0

15 .0 10.0 5 .0 0 .0 0 .0

0 .2

3 .0

< .0

2.0

0 .0 0 .<

0 .6

Reliable models are needed in order to obtain meaningful solutions for the optimal control problem. However, the solution of very refined models, based on Maxwell's equations, requires huge computational resources. In this work, a model of intermediate complexity was considered . This model is based on the application of Lambert's law, taking into account several internal reflections (Van Remmen et ai., 1996), to express the field decay inside the load. These problems were successfully solved, and typical results (optimal controls) for a cylindrical food (potato) load are presented in Figure 6.

50 .0

...

i

Find P(t) and Toven(t) over [to, tf] to maximize the uniformity of the final product quality, as measured by the "cook" C-value, while assuring a minimum reached C-value and a constraint on the maximum temperature at final time, subject to the heat conduction equation in 3 dimensions with a heat generation term

These results are better, in terms of objective value, number of iterations and computation time, than the ones obtained by Schittkowski (1997) using PDECON, a recently developed software for the solution of optimal control problems described by systems of DAEs or one-dimensional, time-dependent PDEs, which is based on the use of a SQP algorithm within the CVP approach.

tl c:

50

Despite its well-known advantages, microwave heating of bioproducts (e.g. foods) also has a number of drawbacks associated with its non-uniform heating patterns. Thus, it would be highly desirable to obtain optimal operating policies for microwave combination ovens in order to achieve better uniformity of temperature and/or quality. Different types of optimal control problems were stated considering the microwave power (P(t)) and the oven (air) temperature (Toven(t)) as the control variables (Saa et ai. , 1998). These problems are very challenging due to the distributed and highly non linear nature of the system 's dynamics. As an example, we present here the statement of one of these optimal control problems:

= 0,

c2 (x , 0) = 0 Cl (0, t) = 10, Cl (1, t) = 10 c2 (0, t) = ai , c2 (1, t) = a 2

35 .0

Oestred

3.4 Optimal control of microwave heating.

with the initial and boundary conditions:

16

_

• - - · ....chle..ed

6 .0

a

j

... ~..

Fig. 5. The desired concentration of substrate vs. the concentration in the optimum.

(23)

.~

j

.

Tin,.

(22)

40 .0

7 .0

0 .2

subject to:

0

0 .0

, .0

o

cl(x,O)

9 .0

~

tf

J

10.0

0 .4

0 .6

0 .6

1 .0

Time

It was found that, using these optimal operating policies, it is possible to achieve a quite uniform final quality or temperature distribution inside cylindrical or spherical loads.

Fig. 4. Optimal control profile for the diffusion through a membrane.

885

~r---------------------------~ '40 T ou t [DEG Cl

'20

P o IwJ!

'00

~

80

I

I'----.J-60

j

40 50

o

~--

o

20

__- -L -__- -_ 300

600

- - - -_____ _ _- 4

!lOO

'200

'500

0 '800

time (a)

Fig. 6. Optimal controls for maximum quality uniformity inside a cylindrical food load.

4. CONCLUSIONS. The results achieved are comparable to or better than the ones obtained with methods based on SQP solvers. Speed-ups of around 5 times and important reductions in number of iterations and evaluations are achieved with the method proposed here, proving that efficiency and robustness are considerably enhanced using exact second order information. This strategy allows the solution of the large scale dynamic systems, resulting from the application of the numerical method of lines to the original combined lumped/distributed parameter system, in reasonable computation times even for high discretization levels of the control variable. The use of reduced order models for the approximation of distributed processes (Sanchez et al. , 1998) will allow the reduction in computation times and thus the use of the TN approach with exact second order information within a model predictive control framework. AKNOWLEDGMENTS: This work was supported in part by the EU (project FAIR CT96 - 1192) and the Spanish Government (CICyT project ALI971939-CE). Author Balsa-Canto thanks the Diputacion Provincial de Pontevedra, Spain, for a predoctoral fellowship.

5. REFERENCES Banga, J . R. , A. A. Alonso, R. 1. P. Matin and R. P. Singh (1994). Optimal control of heat and mass transfer in food and bioproducts processing. Comp. fj Chem. Eng. 18, S699S705. Blatt, M. and K. Schittkowski (1997). PDECON: A FORTRAN code for solving control problems based on ordinary, algebraic and partial differential equations. Technical report. Dep. of Mathematics, University of Bayreuth. Cuthrell, J .E. and L.T. Biegler (1989) . Simultaneous optimization and solution methods 886

for batch reactor control profiles.. Comput. Chem. Eng. 13, 49. Hindmarsh, A.C. (1983) . ODEPACK, a Syst ematized Collection of ODE Solvers. pp. 55- 64. North Holland , Amsterdam. Holmes, P., J. L. Lumley and G. Berkooz (1996) . Turbulence , Coherent Structures, Dynamical Systems and Symmetry .. Cambridge University Press. Leis , J. R. and M. A. Kramer (1985). Sensitivity analysis of systems of differential and algebraic equations. Comp. fj Chem. Eng. 9(3), 93-96. Nash, S. G . (1984) . Newton-type minimization via the lanczos method. SIA M .1. Num. Anal. 21 , 770-778. Neittaanmaki, P. and D. Tiba (1994) . Optimal Control of Nonlinear Parabolic Systems. Marcel Dekker, Inc. Oh, M. (1995). Modelling and Simulation of Combined Lumped and Distributed Processes .. PhD thesis. Imperial College. London. Petzold , L. , J . B . Rosen, P.E . Gill , L.O. J ay and K. Park (1996) . Numerical optimal control of parabolic PDEs using DASOPT. IMA workshop on large scale optimization. Ramirez, W. F . (1994). Process Control and Identification .. Academic Press, Inc., San Diego, USA . Saa, J. , A. A. Alonso and J .R. Banga (1998) . Optimal control of microwave heating using mathematical models of medium complexity. In : ACoFoP (Automatic Control of Good and Biological Processes). Goteborg, Sweden. Sanchez, 1. , J .R. Banga and A.A . Alonso (1998) . Microwave heating of foods : Model predictive control using reduced order models. In: ACoFoP IV ( Automatic Control of Food and Biological Processes). Goteborg, Sweden. Schiesser, W. E. (1994). Computational Mathematics in Engineering and Applied Science: ODEs, DAEs and PDEs. CRC Press, Inc. Schiesser, W.E. (1991). The Numer'ical Method of Lines. Academic Press, New York. Van Remmen , H.H.J. , C.T. Ponne, H.H . Nijuis, P.V. Bartels and P.J .A.M Kerkhof (1996). Microwave heating distributions in slabs, spheres and cylinders with relation to food processing .. Journal of Food Science 61(6),1105- 1117. Vassiliadis , V. S. (1993) . Computational Solution of Dynamic Optimization Problems with General Differential-Algebraic Constraints. PhD thesis. Imperial College, University of London, London , U.K. Vassiliadis, V. S., E . Balsa-Canto and J . R. Banga (1999). Second order sensitivities of general dynamic systems with application to optimal control problems.. Chem. Eng. Sci. 54(17),3851 - 3860.