An extension of Newton-type algorithms for nonlinear process control

Pergamon

1995 Auromatica, Vol. 31, No. 2, pp. 281436, Science Ltd Copyright 0 1995 Elswier Printed in Great Britain. All rights reserved ccm-109m5 $9.50 + 0.00

ooos1098(94)oooa6-7

Brief Paper

An Extension

NUN0

for Nonlinear

of Newton-type Algorithms Process Control* M. C. DE OLIVEIRA

t#t and LORENZ T. BIEGLERt

Key Words-Process control; chemical industry; model reference control; nonlinear control systems; nonlinear programming; quadratic programming.

properties under appropriate conditions. An additional advantage of the method is that it provides a consistent treatment of both linear and nonlinear models. The next two sections describe the computational formulation with a review of the stability properties in Section 3. In particular, we show that in the absence of active constraints, the method becomes a nonlinear extension of the linear quadratic regulator, which guarantees closed-loop stability for a large range of tuning parameters. Finally, the application of this algorithm to various process examples is considered in Section 4.

Abstract-This study extends the multistep, Newton-type formulation [Li and Biegler (1989) C/tern. Eng. Res. Des., 67, 562-5771 for nonlinear constrained process control. The computation of the control law is based on an augmented performance index, which improves the performance of the method over the original formulation. A number of modifications are also introduced in the computational algorithm to extend the range of problems, eliminate steady-state offsets, and extend the output predictive horizon to infinity. We show that the global stability properties of the original Newton controller are preserved by the present modifications. Moreover in the absence of active constraints, the method behaves as an extension of the linear quadratic regulator for, exhibiting therefore excellent stability properties for a large range of tuning parameters. The capabilities of the method are illustrated by application to various process examples.

2. Computational framework The use of operator theory for the synthesis of nonlinear control systems was first suggested by Economou (1985). This proposed Newton approach was later generalized by Li and Biegler (1989) to an SQP algorithm, in order to allow the inclusion of process constraints, and the consideration of a multistep predictive horizon. Here, we consider in this section a modified multistep Newton-type formulation based on a more general quadratic performance index, for a process model of the form

1. Introduction In the last decade, a substantial increase in the use of optimization based predictive algorithms for process control was observed. Model predictive control methods based on linear models have been tested extensively, and are presently well-accepted in industrial practice; they are seen to provide control systems with relatively high performance and reliability (Garcia et al, 1989). Alternative formulations for predictive control using nonlinear process models have also been considered. These extensions have been recently reviewed by Bequette (1991). Nonlinear model predictive control (MPC) algorithms are advantageous in situations where substantial incentives for improved control quality exist, due to the strong nonlinear nature of the system itself or when large changes in the operating conditions can be anticipated during routine operation, such as in batch processes, or during the start-up and shutdown of continuous processes. This paper describes a nonlinear control algorithm that extends the basic Newton formalism of Li and Biegler (1989). A number of important characteristics demonstrated by this approach make the algorithm particularly attractive for process control. These include a computationally efficient method for the determination of the model sensitivities and the successive quadratic programming (SQP) nature of the algorithm, together with the existence of global convergence

2 =&(x, u, d),

(la)

Y = g,(x).

(Lb)

where u E R”e represents the vector of system inputs (manipulations), x E Iw”s is the vector of state variables, y E lR”o is the vector of system outputs, and d E LW is a vector of disturbances. It is assumed that n, s ni, that f, and g, are continuously differentiable functions. In order to simplify the determination of the optimal control laws, a discrete-time framework is used. Denoting the present sampling instant by tk, the correspondent input and output predictive horizons are designated by %$ = [tk, tk + rik] and %& = [tk, cl, + t,J, respectively. This allows the correspondent predictive problem at tk to be formulated as

minJZi("7Y) =

u* I,%

2 (Yk+i- Y,p,k+i)TQyi(Yk+i-Ysp,k+i)

i=l

+ hk+i-I - u,,k+,-dTQui(Uk+r-, - u,.k+i--l), (24

* Received 25 April 1993; revised 10 January 1994; revised 20 March 1994; received in final form 14 April 1994. This paper was not presented at any IFAC meeting. This paper was recommended for publication in revised form by Associate Editor B. Wayne Bequette under the direction of Editor Yaman Arkun. Corresponding author Professor L. T. Biegler. Tel. +I 412 268 2232; Fax +l 412 268 7139. t Department of Chemical Engineering, Carnegie Mellon University, Pittsburgh, PA 15213, U.S.A. $ Present address: Department of Chemical Engineering, Universidade de Coimbra, Large Marques de Pombal, 3000 Coimbra, Portugal.

i= fp(x,u, 4,

PI

Y = &Ax).

(2c)

(24

Y/SY 'Y"? U,5U~U,

and

u(t) = u(rk + r,, - T); t E [zk + tih, rk + tohI. (2e)

a reference input trajectory, Here, u, E R”a represents accomplishing the same role as y,, for the outputs, and compatible to these set-points. The notation can be simplified 281

282

Brief Papers

by de&ring the augmented predictive trajectory, as

input and output vectors for a

been satisfied. If so, implement uk = V,, and update ii for the next interval as V=[ii;+,

and a nominal input trajectory

u:+,,,~.,]‘r.

ii E lR”P’= [UZ i;+]

At each iteration, the difference between the current and reference input trajectories is represented by AU = V - V. Since the Newton framework requires a linearized approximation of the system equations, we consider the expansion _of Y in Taylor series around the nominal trajectory Y Y = Y + Y”,AV,

(3)

with Y, = aY/X//,=, E 5?(n~p)x(n~m),The dynamic matrix y, represents the sensitivity of the process to input changes around the nominal trajectory considered. It can be expressed in terms of the derivatives of the state transition function of (la), leading to a lower block triangular matrix Ym defined by 0

zk+,=zk+r~-,

S,, otherwise,

i g if i = j, S = Ck+,Ik+, 1 1, i C,+,(II$L, Qk+,-,,,)Ik +, ~I otherwise.

(5)

Here C corresponds to the derivative of the output function, while @ and F represent the sensitivities of the nonlinear model (la) to changes in the initial conditions and inputs, respectively. These two quantities can be computed exactly by integration of a set of linear equations together with the process model (Li and Biegler, 1989). Using the linearized model given by (3) and defining V’ = V, - V, with E = ySp- Y, this can be expressed as the following quadratic program (QP):

Liz+,-,

U:+,_,]r.

+k,(y,+-~s~,~+r),

i = 1,.

,p,

with zk E R”o, k, E Iw”~““~, is added to the discrete model. Using the deviation variables denoted in a generic form by the ’ superscript, with u’ = u - ri, x’ =x -X, the first-order approximation of this state equation is given by ik+r=fk+r-,

+kl(~k+r+Ck+i(4k+r~,X;+r-l -

k1.vsp.k

+rk+,-&--1))

t,.

Similarly, defining z;+, = ik +, - Z,,+,, leads to the following augmented LTV approximation for the system equations:

(74

-L&AU)

+ (AV - V’)TQ,(AV V,d c AU 5 Cl,,<,, ylds

li:+,mz

Zk = zo.

(4)

~tn,j =

...

Go to Step 1, with I = r,,,. Step 6. If the convergence criteria for V are still not satisfied, check the QP iteration counter i. If i
ifj > i, ifj
I-$LI~=(E-L&AV)~Q,(E

li:+z

Y,AV 5 Y,d,

~ V’),

(64

0)

(6b)

(6~)

with lJld = V, - 0, VU, = V,, - 6, V’ = V7 - i?, Y,d = Y, - y, and Y, = Y, - Y. Here Q, = diag {Q,,} E Rg(n~p)x(n~p),and Q2 = diag {Q,} E lRtn~m)xtnlm?In comparison to (6) the objective function considered by Li and Biegler (1989) corresponds to Q, = Inop and Qz = 0. However, the use of the augmented performance index (2a) gives the method improved stability properties, as discussed in Section 3. The algorithm for the determination of the nonlinear control law is stated below. It assumes that each nonlinear optimization problem for a given predictive horizon is iterated until convergence, and that the value of the state vector xk is also available to the control algorithm. Conrrol algorithm (converged form). Using the initial conditions tk = to, X* = xg and U = V,: Step 1. At t= I~, set i = 0 (iteration count for the QP subproblem). Step 2. Starting with the initial condition x~, integrate the state and sensitivity equations (la) in the output predictive horizon %&, for the nominal trajectory V = U. This gives x, and Ik+,. Compute also the output sensitivities Ck+l, t_;

The correspondent augmented predictive model equations can therefore be written as Y = Y + y,,,AV, where p = [yT z;lT, y”=_ky ?3=, (i = k, , p). The augmented dynamic matrix Y, is obtained from (5), by replacing the sensitivity coefficients C, @ and I by -their augmented counterparts in (7). Defining also Qy = diag {Q,, I,,$ Ysr = rY,& OIT,allows the control algorithm described m Section 2 to be used directly. 2.2. Control law computation for an infinite output horizon. It is well-known that large predictive horizons make the stability properties of the resultant control law less dependent on the tuning parameters of the algorithm. The previous algorithm can be extended to use an infinite output horizon, by introduction of a number of rearrangements that deal with the fact that some of the previous quantities would have an infinite-dimension in this case. We assume here, that the open-loop system to be controlled is already largely asymptotically stable. Since the formulation requires a set-point for the infinite horizon, we assume also that the specified set-point trajectory converges to a finite and feasible limit y*. In this case, it can be shown that the linear and quadratic terms of (6a) are given by

>...>P). _ Step 3. From X compute Y and E, update the deviation bounds Vid, VU,,, &,, Yud and solve the corresponding QP (6)

for AU. Step 4. Using a line search, determine an appropriate step sjze h_ for the current input correction dV = AAV. Update VcVtdV. Step 5. Check whether the convergence

criteria for V have

W

b, = -2 2

ST,Q,,E,

- 2QY,V,‘,

i, j = 1,.

*=, where bi E l??~ and H, E Iw”~““~represent

, m,

(8b)

the linear and

Brief Papers Hessian block elements of (6a), corresponds to the Kronecker delta

respectively.

Here

while Sjj is given by (5). A more detailed treatment issue is given in Oliveira (1994).

S,,

of this

3. Stability properties Results in this section comprise a generalization of the stability results of Li and Biegler (1989) for OL-stable systems and a general unconstrained treatment. The effects of the presence of active constraints on the closed-loop stability, including some soft-constraint handling strategies are presented in Oliveira and Biegler (1994). For open-loop asymptotically stable systems, Li and Biegler (1989) demonstrate that with Q, = Z, Qz = 0 the controller can be made globally stable, by choosing a sufficiently long sampling time, and using a line-search algorithm to enforce a descent property. The next theorem shows that similar properties exist for the formulation presented in the last section. Theorem 1. Assume that the open-loop nonlinear system described by (1) is largely asymptotically stable, and that a constant set-point y,, = y* is used. For (1y, - y*llo, 2 6 > 0, the application of the multistep Newton-type controller correspondent to the previous algorithm with m,p > 1 converges monotonically for sampling times T > T*, where T* is finite. Therefore we have II Ys+I -Y*l12a,
proof for this theorem can be found in the Appendix. It relies on the fact that for open-loop asymptotically stable systems, a sufficiently large but finite T can be found such that 119115l, for any small predetermined E. Therefore choosing T > T*, where

guarantees the stability of the Newton controller with open-loop asymptotically stable systems. This situation is representative of the characteristics of the majority of chemical processes. However, since the above results rely on the boundedness of the sensitivity coefficients, they do not extend to the OL-unstable case. In fact, a strong statement about the non-existence of global stability properties of unstable systems can be made (Sontag, 1984). Also, the requirement of a sufficiently large sampling time might not always be readily met in practice. This motivates us to look closer at the stabilization properties of the underlying state-feedback structure. In the absence of constraints, the predictive control problem at r = rk is given by (6a). The optimal solution for this case can be obtained in analytical form, leading to the control law f~ =

(CQ,%

The first one is a Continuous Stirred-tank Reactor (CSTR), and the second is a Fluid Catalytic Cracking (FCC) unit. Both of these examples exhibit interesting characteristics in terms of nonlinear behavior over an extended range of operating conditions. 4.1. CSTR example. The system considered here consists of an ideal CSTR with an isothermal pseudo first-order reaction A + B +P, in the presence of an excess concentration of A (Oliveira, 1994). The corresponding process model is

~=(c,,-**)~+(c”‘-“~)~-~.

We assume that both state variables are measured and that the process output is given by {x1, x2}. The vector of manipulated variables is {u,, uZ}, while Car and Ca, constitute potential input disturbances. Under nominal operating conditions this model can be shown to have three possible equilibrium points, with the middle one being OL-unstable. For simulation, we used a sampling time T = 1, predictive horizons with 10 sampling intervals and Q, = diag{l, l}, Q, = diag (0.01, 0.01).Figures 1 and 2 illustrate the response of the closed-loop system to a +lO% step change in C,,, at t = 1. Despite being unmeasured and unpredicted by the control algorithm, this disturbance is well rejected; for comparison the correspondent open-loop response goes to a steady state of x~,+~= 9.4. Since the disturbance is not explicitly accounted for by the controller, the final steady states of the output variables exhibit a small offset when no integral action is used (about 0.4%). This offset is, however, eliminated when integral action is added, as can be observed from the corresponding output profiles. Using a tolerance Ed” = 0.5 X 10w3, the number of iterations required reached a maximum of five or six during the sampling times where the major change (disturbance or set-point change) occurred, and was typically smaller after that (around two). This resulted in a total computational time of about 1 minute of CPU time for each case, using a mixture of ‘Mathematics’ (Wolfram, 1988) and Fortran code, on DecStation 5ooO/240. 4.2. FCC example. Next we consider the control of a FCC

I.

I04 102 ;

100 .U

4. Examples We consider, in this section, the application of the previous algorithm to the control of two nonlinear process examples.

. .. .

-*-h-*-c.+-.++-

98 96 J

+ Q&‘(%Q,E* + QzU,),

with E* = Y,, - Y*, where Y* corresponds to the system response for a zero input. Hence, the extended Newton formulation behaves in the unconstrained case as a generalization of the discrete Linear Quadratic Regulator (LQR) control for nonlinear models, possessing therefore similar properties to this last approach, with both infinite horizons and finite horizons with additional stabilization conditions (see Oliveira (1994) for additional details). These modifications allow the unconstrained Newton-type controller to achieve closed-loop stability with a perfect model, without requiring an extensive search for the appropriate tuning parameters. As mentioned, a detailed treatment of the effects of the presence of active constraints can be found in Oliveira and Biegler (1994) and Oliveira (1994).

(W

0

5

10

13

20

15

20

t : ,.i 0.68 0.66 2

0.64 0.62 0.60

-I

”

5

10

t Fig. 1. Output response to a 10% step change in CB1 (. . . open-loop; --- k, = 0; - k, = {l, 1)). The dots indicate the discrete sampling points used.

284

Brief Papers 1.05

Table

0.85

I. Steady state values for the FCC model: OL-unstable point; (b) stable steady state

(4

(b)

Steady state values

Unstable

Stable

C,, (lb coke/lb cat) r,, (OF) C,, (lb coke/lb cat) T’p(“F)

7.985 x 10-s 957.6 2.347 x lo- ’ II63

7.167 x 10-a 1149 3.589 x lo-“ 1396

(a)

0.80

1.1

2

1.05

1.00

0.95

I

0

5

10

15

20

Fig. 2. Input profiles for a 10% step change in CHI (. open-loop; --- k, = 0; - k, = {l, 1)).

process. These units are commonly used to convert heavy petroleum feedstocks into lighter hydrocarbon products; a key step in petroleum refining. The present study is based on the semi-empirical model of Lee and Kugelman (1973). which consists of balance equations for the mass of coke (carbonaceous material) and energy, both in the reactor and regenerator vessels (Oliveira, 1994). The state variables are C, (coke content in the spent catalyst), TX (reactor bed temperature), C, (coke content in the regenerated catalyst) and Trs (regenerator bed temperature). The model considers also five main input variables, which are F, (air flow rate), T, (air temperature), F; (feed rate), G (feed temperature) and 6 (catalyst recirculation rate).

The selection of control objectives was based both on the state variables with more direct effects on the reaction environment and potential difficulties on their measurements. We adopted in this study T,, and C, as controlled variables, while using F, and F, as manipulated variables. For simulation, a sampling time of 0.5 hours and predictive horizons with 5 sampling intervals were used. The tuning parameters in the objective were chosen as QY= diag {l, O.OOl}, and Q, = diag {O.l, O.l}. Figure 3 illustrates the behavior of the process during an unanticipated change in the set-point, from point (a) to (b) of Table 1. In the absence of saturation of the input constraints, the new set-point is reached at the end of the second interval, with just a small overshoot in C,,. This indicates that a more cautious response can also be achieved, by imposing additional limits on the rate of change of the manipulations, if desired. The corresponding input profiles are presented in Fig. 4. More iterations were required in the present case to converge the input profile to optimality, in comparison to the previous example. The maximum number of iterations during a given predictive horizon was 15, over the intervals where the largest changes occurred. This number decreases to 2 iterations at the end of each profile, when the final value of the response is approached. Each run required around 5 minutes of CPU time, in a DecStation 5000/240, for a control horizon of 10 hours. 5. Conclusions

An extended formulation of a Newton-type algorithm for nonlinear model predictive control was presented. The

1200 722500 1150 722000 1100 ,K +

d

721500

1050 721000 1000

720500

950 J 0

I

2

3

4

I

5

2

t

3

4

5

3

4

5

t

0.0025 9.70 106 0.0020 9.65 IO6 p v

0.0015 z 0.0010

9.60 lo6 9.55 106 9.50 106

0.0005 L

0

I

2

3

4

5

Fig. 3. System response to an unanticipated set-point change from the operating point (a) to (b).

0

I

2

Fig. 4. Input profiles for an unanticipated from (a) to (b).

set-point change

285

Brief Papers proposed approach is based on the use of an augmented quadratic performance index and leads to a computationally efficient framework for nonlinear model predictive control. In addition, the inclusion of an integral term to eliminate steady-state offsets, and the possibility of extending the predictive horizon to infinity are considered. The stability properties of the method are considerably improved by the extended formulation. In the unconstrained case, the controller becomes a nonlinear extension of LQR, which guarantees good stability properties with a large range of tuning parameters. The algorithm was successfully applied with two nonlinear models. In particular, the behavior demonstrated with the FCC model illustrates the efficacy of the method for challenging nonlinear systems. Simultaneously, the small computational times required with both examples demonstrate the feasibility of the on-line implementation of this nonlinear strategy. Acknowledgements-Financial support from the Department of Energy, in the form of grant DE-FGOZ85ER13396 is gratefully acknowledged. The first author is also grateful for additional financial support provided by the University of Coimbra and FLAD. A more complete description of the control alaorithm is available as Technical Renort 06-157. Engineering Research Design Center, Carnegie Mellon University, Pittsburgh, PA. This report is also available on the Internet by anonymous ftp from ftp.deq.uc.pt (IP address 192.84.13.35). References

Bequette, B. W. (1991). Nonlinear control of chemical processes-a review. Ind. Eng. Chem. Res., 30,1391-1413. Economou, C. G. (1985). An Operator Theory Approach to Nonlinear Controller Desian. Ph.D. Thesis. California Institute of Technology, Pasadena, CA. Garcia, C. E., D. M. Prett and M. Morari (1989). Model predictive control: theory and practice-a survey. Automatica, 25,335-348. Lee, W. and A. M. Kugelman (1973). Number of steady-state

operating points and local stability of open-loop fluid catalytic cracker. Ind. Eng. Chem. Process Des. Develop.,

From (A.l) and (A.2) J2W

II%+, -

=

Y*IIzp~

IIYk -y*112py+ O(a).

(A.3)

Assuming now that O< 6
(A.4)

Now substituting (A.3) into (A.4) and for E sufficiently small that O(c) < 612, we have J2(A~~)~Ily~-~*ll1p,+~~~)-~A~ll~~-~*ll~~-~.

(A.5) Also llyk+i

-y*ll’p+

A2 IlAu~ll&~Jz(AW

+ WA*).

64.6)

Choosing A so that O(A’) < SA/4, we have from (A.5) and (A.6) the desired result for m = p = 1: llY*+1 -y*v

q 5 llyk.+1-

y*ll$~+

A2 IIWl&,

Note, that if the control problem is solved to optimality in each horizon then Jf 5 J2(AAuk) < J2(0) it Jz follows from the NLP, and the desired result is also obtained immediately. In the multistep case, optimality conditions of (6) can be written as

2(CQ, % + Q&W - ‘WtQ, E -

2Q2U’ - Ay,) = 0,

(A.7a)

A:u(AU - U,) = 0, A:,(-AU + U,J = 0, AT&&AU - Y&) = 0, A;,(-&AU + Y,‘,)= 0,

(A.7c)

+ (Au,, - Au/) + %(A,,

A,ZO, AU E ]I&, U,],

(A.7b) (A.7d)

X,AU

l

[Y/d, Y&l.

(A.7e)

Also for large but finite T such that ]]ek ]I = O(e), the matrix 9, becomes

l2,197-204.

Li, W. C. and L. T. Biegler (1989). Multistep, Newton-type control strategies for constrained nonlinear processes. Chem. Eng. R&

Des., 67, 562-577.

Oliveira. N. M. C. and L. T. Biealer (1994). Constraint handling and stability propertie: of ‘model predictive control. AIChE. J., 40, 1138-1155. Oliveira, N. M. C. (1994). Newton-type algorithms for Nonlinear Constrained Chemical Process Control. Ph.D. Thesis, Carnegie Mellon University, Pittsburgh, PA. Sontag, E. D. (1984). An algebraic approach to bounded controllability of linear systems. ht. 1. Control, 39, 181-188. Wolfram, S. (1988). Mathematics. A System for Doing Mathematics by Computer. Addison-Wesley, Redwood City, CA. Appendix Proof of Theorem

1. Here we start by looking at the single controller (with m = p = l), and examine a single iteration of the Newton algorithm. The solution from (6) is Auf, and for a stepsize A E [0, l] and a constant set-point ysp =y*, the objective is defined by J2Wuk)

=

lIPk+l- y*ll'a+ A2 IIA~,cl12ak.

(A.11

(‘4.8)

For m > 1, taking into account the block diagonal structure of Qi = diag {Q,}, Q2 = diag{Q,}, the first equation in (A.7b) becomes

2[(C,+,r,)TQ,(C,+,r,) + Q, + o(~~)lAu, + O(E)bUk+,+. . . + O(E~-~)AU~+,,-, - 2](C,+,I,)TQ,e,+, + O(e)] - 2Quu; + (L - L/h + [(~k+lwTQy(h.u - bh + wi = 0, with ek+, = y* -yk+,. Note that AU remains always bounded as l --) 0, since (ymQ, ym + QJ is positive definite. Assuming that (6) has a feasible solution, the previous equation can therefore be written as AU, =

We assume now that a sufficiently large sampling time but finite T is used, so that ]]@,I]= O(E), for a small predefined positive quantity E. From (3) we have yk+, =&+, + AC, 4,TkAu,, and

[(ck+lrk)TQy(ck+lrk)

+

QY]-’

(Ayu +

=

Q

”

&

_

(huu

Au* + O(c),

-

2

*“h

1 +

2 - Ay’)‘)

Ote)

(A.9)

Brief Papers where Au* represents the solution of the single step controller. A large T, therefore, has a decoupling effect in the computation of the control law (6), and thus yk+, =yz+, + O(e). Clearly J2(Auk) =&(A@) + O(8), and from (AS)

J,(hAu,) < jlyk -y*l&+

O(e) - Sh + O(e*).

Hence, for sufficiently small E, such that O(r) < 6A/2, the remainder of the previous proof applies, and we have the desired result from (A.6). 0