Model predictive control for nonlinear batch processes with asymptotically perfect tracking

~

Computers chem. Engng, Vol. 21, Suppl., pp. $873-$879, 1997

Pergamon

© 1997 Elsevier Science Ltd All rights reserved Printed in Great Britain

PII:S0098-1354(97)00159-2

0098-1354/97 $17.00+0.00

Model Predictive Control for Nonlinear Batch Processes with Asymptotically Perfect Tracking Kwang Soon Lee*

Jay H. Lee

Dept. Chem. Eng., Sogang Univ. Shinsoo-1, Mapogu, Seoul 121-742, Korea

Dept. Chem. Eng., Auburn Univ. Auburn, AL, 36849-5127, USA

A b s t r a c t - We propose a novel model predictive control(MPC) algorithm designed specifically for nonlinear batch or other repetitive processes. Unlike existing MPCs which duplicate the same control error in the repeated batches, the proposed algorithm can achieve perfect tracking (for square systems) despite model uncertainty as the number of batch runs increases. As a special case, we also propose a maximum horizon algorithm which ensures convergence just as infinite horizon MPC does for continuous processes. Numerical examples are provided to demonstrate the performance of the proposed batch MPC algorithm.

1

Introduction

From the viewpoint of control system design, batch processes present a unique challenge seldom found in continuous processes. In a batch process, the control problem is usually given as a tracking problem for time-varying reference trajectories defined over a finite interval. During the course of a typical batch, process characteristics very often change with time (for example by reaction) and at the same time process variables swing over wide ranges. All these may cause the process to significantly expose its nonlinearity and prevents us from finding a reliable process model. Because of this fact, the linear control techniques have very often proved unsuccessful in batch process control. Through the evolution for the last fifteen years or so, model predictive control(MPC) is now regarded by many as one of the standard control paradigms for industrial processes (Lee, 1996). MPC has many attractive features, which well suit to industrial demands, such as easy accommodation of conflicting control requirements of a multivariable system in an optimization criterion, constraint handling, decision based on prediction, etc. Despite the reputation in continuous process control, however, MPC may not fit batch process control with its current form. By the model-based nature, MPC under tracking control cannot effectively cope with various uncertainties although the imbedded integral action can remove offset during regulatory control. In addition, the tracking error will be duplicated through batches since only the current measurements are fed back to the control law. One of the aspects of batch operation often left unexplored is that it is repetitive. Exploring this aspect (for example, feeding back the measurements in the past batches) may give us a new opportunity to renovate batch process control. Based on the above considerations, the objective of

the present study is placed on developing a novel MPC algorithm for nonlinear batch and other repetitive processes, where both precise tracking and disturbance rejection are important. The key aspect of the developed algorithm is that measurements not only at the current time but also from the previous batch are fed back for control input calculation. With this idea, the proposed batch MPC algorithm has several outstanding traits that are not found in existing MPCs. The most important consequence among them is that minimum possible tracking error can be attained against unknown repeated disturbances as well as model error as the batch number increases. As a special case of the algorithm, we also propose and investigate the so-called maximum horizon batch MPC algorithm. This algorithm guarantees convergence for any process model as infinite horizon MPC in continuous process control. The rationale of proposed batch MPC is indeed borrowed from quadratic criterion-based iterative learning control(Q-ILC) the authors have developed recently as an input trajectory refinement technique for batch or other repetitive processes (Lee et al., 1997). We will briefly review Q-ILC in a subsequent section.

2 2.1

General Framework Model Development

We consider an nu-input/nu-output discrete-time nonlinear batch process. Since a discrete-time batch operation is defined over finite time indices, it is more convenient to consider the process as an algebraic system which relates a finite input sequence to a finite output sequence defined over the underlying discrete-time domain. If we use u, y, and d to denote the input, output, and disturbance vec-

*To whom all correspondence should be addressed: tel (82-2)705-8477, fax (82-2)3272-0319, kslee~ccs.sogang.ac.kr $873

PSE '97-ESCAPE-7 Joint Conference

$874

tors (including initial conditions), respectively, and assume additive random disturbance to the output, the algebraic relationship at the k th batch of the process can be written as Yk = Y k + V k = A / ( u k , d k ) + V k

(1)

where yT

=

[yT(1) yT(2)...yT(N)]

uT

=

[uT(0) uT(1) ' ' . u T ( N -

1)]

(2)

and likewise for Y, d, and v. Through Taylor expansion, we derive a lineaxized output error transition model between two successive batches. Yk+l = Yk + GkAUk+I + DkAdk+l + . . .

(3)

where

Auk+l

= u k + l - uk

(4)

and likewise for Adk+l. Let Yd denote the desired output sequence yT = [yT(1 ) yT(2). ' .yT(N)]

(5)

and define

2.2

R e v i e w of Q-ILC

Q-ILC is a version of iterative learning control(ILC)(Moore, 1993), where the open-loop input trajectory is successively refined using previous batch data. It can be formulated under several different settings(Lee et al., 1997). In this subsection, we briefly introduce observer-based Q-ILC on which the present research is based. The objective of QILC is written in the deterministic context as

Given the model in (7), design a learning algorithm with the property of

-~ min II~ll~

II~kll~

as k -~ ~

(8)

when Awk = 0 and vk = O. In order to achieve the above objective, Q-ILC solves the following minimization at the end of the k th batch: min 21 {ll~k+xll~ + II~nk+l I1~ Izk } ~u~+,

(91

where 77k represents the information at the end of the k th batch; Q and R are PSD (positive semidefinite) and PD (pesitive definite) matrices. Note the cost function has a penalty term on the input change between two adjacent batches, which provides an integral mode in the resulting control law. Unconstrained Q-ILC Algorithm

ek

:

Yd -- Yk

ek

=

Yd -- Yk

Awk+ 1

=

DkAdk+l

(6) +...

then the following locally linearized model is obtained in a state space form: ek+l

=

ek -- G k A u k + l

ek

:

ek q- Vk

Remarks

By the certainty-equivalence principle, the solution can be obtained by solving min 51 {ll~k+11kllQ + IIAuk+~IIR Iz~ } Au~+,

(I0)

where the state estimate ek+llk is obtained from the following state estimator:

q- A W k + l (7)

:

1. dk which is persistent over batches is eliminated in the above model. Implication of this is that such a disturbance plays as a bias term and an integral control over the batch index, k, will be able to remove its effects as in continuous process control. 2. By causality, G is represented by a lower block triangular matrix with ny × n~ time-varying pulse response coefficient matrix g(i, j) as its (i, j)th block element. 3. vk can be assumed to be zero-mean independent random since zero-mean output disturbance vectors from two different batches can be reasonably assumed to be independent no matter how the disturbance is correlated in a batch sequence. Likewise, if we consider the case (for example, batch reactors) where the same d persists at least for several batches with some additive random perturbation and Gk is found through an unbiased model identification technique, we may presume zero-mean independence of Awk, too.

eklk

:

ek]k-1 q- K (ek - eklk-1 )

ek+llk

=

~kl k -- G k A U k + l

(11) (12)

The resulting algorithm is uk+l = uk + HQk~klk

(13)

The observer gain K can be considered as a de facto tuning factor. If statistical informations of Awk and Vk are available, it can be replaced by an optimal gain. For instance, if Vk and Awk are zero-mean i.i.d (independent, identically distributed) sequences with covariances of V and W , respectively, the following steady state Kalman gain can be used:

K

=

Ppd(Ppd + V ) -1

(14)

P

=

W+WP-1V

(15)

where Ppd is the positive definite solution of the algebraic Riccati equation in (15). The above Q-ILC algorithm has the following important properties: convergence: Under the assumption that (7) represents the true process and Awk = vk = 0, (8) is achieved regardless of the choice of Q and It. This

$875

PSE '97-ESCAPE-7 Joint Conference is true for constrained algorithms, too under some legitimate assumptions on the constraints. robustness: (8) can be shown to be fulfilled even when there is model error within certain limit. The limit can be expanded by increasing t t at the expense of convergence rate as long as sign reversal of the gains of G (in terms of singular values) doesn't occur. disturbance sensitivity: "Learning" is a procedure to find the model inverse in a progressive way using input/output data. Therefore, ILC is apt to sensitively respond to high frequency components in the outputs. The high frequency sensitivity can be adjusted at will without losing the desired asymptotic convergence through adjusting R, but at the expense of convergence rate. The same effect can also be obtained by adjusting the observer gain. 2.3

Incorporation trol

of Feedback

Con-

Although Q-ILC has many attractive features, it provides only an open-loop input signal and lacks real-time correction capability. To use Q-ILC in real processes, feedback control needs to be supplemented. In general, feedback control can be described by the following two-degrees-of-freedom equation: u(t) = Hl(q-1)yd(t) -- H2(q-1)y(t)

(16)

by inaccurate control input. To avoid this problem, we have to redesign the predictor such that correct output prediction can be obtained in spite of model error by feeding back previous batch information as in Q-ILC.

3

Batch MPC

3.1 3.1.1

Algorithms

Unconstrained

Batch

MPC

A l g o r i t h m Derivation

Assume that we are at t in the k + 1th batch, that is, we just measured Yk+l (t) and have to determine Uk+l(t). Let the prediction and control horizons be chosen as p and m, respectively. As was reviewed in subsection 2.2, the output error estimate in Q-ILC converges to the true output error notwithstanding model error within a certain bound. This property comes from the integrating process structure (with respect to the batch index) of (7) where previous batch information is fed back. In order to make the most of this property, we derive the MPC predictor as a partition of the observer equation in (11) and (12). For this purpose, we forward-shift the batch index of the a posteriori update equation, (11), by one. Now, to simplify the presentation, we first define the following vector and matrix partitions:

Addition of feedback control to (13) will yield the following combined algorithm: Uk+,

~--" U k "~

H?fiklk q" Sl,kYd -- H2,kYk+l

(17)

Au(t)

=

: L,X,,(t-

From this equation, we can see that two-degreesof-freedom control results in output offset. When Awk = vk = 0, convergence (to some limit) of the controlled system necessarily implies convergence of Uk. This in turn tells ( H O + H 1 ) y d = (HQ+H2)yoo, which leads to Yk ~ Yd unless H I = H2. From this consideration, we note that feedback control should be of an error-driven type for zero tracking offset. Error-driven PID control can be one of possible candidates for the feedback algorithm(Lee et al., 1994). Blind C o m b i n a t i o n o f M P C w i t h Q - I L C Among a lot of different control techniques, MPC is, of course, one of the most attractive candidates for the feedback algorithm. Owing to the modelbased predictive nature, however, regular MPC has some prospective problems. In the standard MPC design procedure, prediction is made by incorporating the feedforward signal together and the control input is calculated to compensate for the prediction error. The consequence is that an open-loop feedforward bias signal does not change the resulting control input. In any case whether the learning input is given or not, MPC generates virtually the same control signal according to the performance criterion. When the model has uncertainty, prediction will be inaccurate and tracking error will be produced

E(t) =

e't)

,

AU(t)= I

1)

, E(t)=

:

LA,,(t

-

L

e(t +:P)

1)

(IS)

(11 (t), G2 (t), Ga (t) = partitions of G associated with

(E(t), AU(t)), (g(t), AU(t)), and (g(t), AU(t))

pairs,

respectively Also, we use ^ to denote an estimate, - and + to represent the a priori and a posteriori estimates replacing the conditional notations k+l]k and kWllk+l, respectively. For further simplicity, we omit the argument t where there is no lconfusion. According to the above definitions, the output error prediction over the prediction horizon is denoted by E++I" Using (11) with the batch index shifted by one, the output error prediction can be described by

where E is the partition of the observer gain K corresponding to (gk+l, Ek+l) pair. K: can be regarded as a tuning parameter. When K is the Kalman gain, however, K: becomes the optimal predictor gain under the condition that measurements up to t at k + I th batch are available.

$876

PSE '97-ESCAPE-7 Joint Conference

~k+l a n d Ek+l in (19) can be further expanded using (12) as functions of future control movement and other available information. ~k+l = ~ : -- G2,kAUk+ 1 -- G3,kAUk+ 1 Ek-t-1 = Ek+ -- G1, k m u k + l

(20) (21)

In standard MPC (for example, QDMC), input rate is penalized in the performance criterion to provide an integral action for elimination of zero-frequency offset. Likewise, elimination of tracking error in batch process control is possible by providing the contol law with integral action over the batch index (penalizing input change between two adjacent batches). From this stipulation, the following quadratic criterion is considered for input calculation:

min 1 {[ig++ll[22 +HAHk+11122}

AL4~+I

3.2

MPC

(22)

• -i {IIG3kAHk+ 1 - '~"~k+l1A2 2 +,, ,.IlAU~+IlI-,j~'~ mm .Au~+t 2 (23) where

(24)

In the above, ~ - 1 represents the output error prediction when no control is implemented. For unconstrained case, the following closed-form solution is readily obtained. AUk+ 1 = (G3,kA2G3,k T "-rF 2)~ - l f'~J3,k~X2Vk-{-1 , T ~ ~uc (25)

From AHk+I, the first element is taken and implemented. The same calculation is repeated with the time index. 3.1.2

Batch

, , F JAU,~+~I [GTkA2(G~ k -- EGI,k) G3T,,kA2G3,k "1- 2][AUk+I J

=

^

Horizon

As a special case of the proposed batch MPC algorithm, we consider an algorithm with maximum horizons. In this algorithm, the control as well as prediction horizons are constantly changing to have the maximum length up to the batch terminal time, i.e., p = m = N for Au(0); p = m = N - 1 for Au(1), and so on. We will call this algorithm maximum horizon batch MPC(MHB-MPC) and show it has an exquisite convergence property similar to infinite horizon MPC(Muske and Rawlings, 1993) for continuous processes. MHB-MPC can be rewritten in a rather compact form. After some straightforward manipulation, (25) can be rearranged as

Substitution of (19) and (20) into (22) yields

^uc

Maximum

Properties

Offset From (25) and the related equations, we can note AHk+I is represented by a linear function of control error and its estimate, i.e., an error-driven form. This is also true for AUk+l which is stacking of the first element of AHk+l(t) for t E [0, N - 1]. Therefore, we can see no tracking offset will occur if the algorithm converges.

GLA2][ E+

E,+, + [ GLA2 o] [ E +I ]

(26)

Owing to the maximum horizon policy, we have e T -- [FT vcT] and A u T = [AU T AuT]. Also we note that (26) is reduced to (13) at t = 0 since A u = AH and e = C at the initial time. This implies Auk+l (0) is determined by the Q-ILC equation alone. Now, it can be easily verified that the three coefficient matrices in (26) are the [ (2, 1) (2, 2) ] partitions of GTA(I - / C ( t ) ) G k + F, G T A ( I - / C ( t ) ) , and GTA/C(t), respectively where

A=IAx 0 ] 0

A2

, F=

It, 0] 0

0]

F2

/C 0

with K:(0) = 0. In the above, we use the argument t only for )C to explicitly show it is redefined with t while other matrices remain same. From this, we can see that AUk+l(t) is the t + 1th element of the following formula: AUk+l

= (GTA(I - E(t))Gk + F) -1 x

GTA ((I -/C(t))~kl* + ]~(t)ek+l)

(27)

MHB-MPC inputs over the whole batch horizon can be summarized in a vector form as AUk+l = H~eklk + Hk/ek+l

Convergence As demonstrated through an example in section 4, batch MPC with arbitrary p and m does not always converge even when the process is linear and the exact model is given. Through extensive numerical simulations with linear time-invariant processes, we found that stable linear processes always converge while a part of unstable processes show divergence. This behavior is very similar to that of standard finite horizon MPC.

[-GLA

(28)

where N-I

N-1

H~ = E Hk(t)(I --/C(t)), H I = E Hk(tlKZ(t) t=0 t=0 Hk(t) = N ( t ) ( G T A ( I - ]~(t))Gk + F ) - I G T A (29) and N(t) is a block square matrix with zero elements except I at (t + 1, t + 1) position.

$877

PSE '97-ESCAPE-7 Joint Conference 3.2.1

Properties of MHB-MPC

Offset

As a special form of batch MPC, MHB-MPC will not result in tracking offset if converges.

constrained minimization for constrained case. Implement the first element of AUk+I. Set k = k + 1. step 2 - batch u p d a t e Update @klk using (11) and (12). Update Gk if necessary. Go to step 1.

4 Convergence

Dynamics of the controlled system can be described with respect to the state, ek, and observer error, eklk = ek -- eklk, by substituting (28) into (7), (11), and (12). From this result, convergence of the controlled system is found to be determined by I - K and (I + G k H k / ) - l ( I - GkHq). Since (7) is observable, we can assume I - K is exponentially stable without loss of generality. What remains to show is (I + GkHk/) -1 (I - GkH~) is stable. Through extensive numerical test, we found that (I + G ~ H I ) - I ( I - GkH~) is independent of K and can be represented by the following simpler formula: (I + GkHf,k) - 1 ( I -- GkHq,k) = I - G k L k

(30)

where Lk = ( G / A G k + r ) - I G T A

(31)

This tells the closed-loop dynamics is independent of the observer dynamics, which is a well-known property in linear control systems. Using this property, the following theorem is derived: T h e o r e m 1 Under the assumption that Awk = vk = 0 and the range space of Gk is invariant for all k, the process of (7) under M H B - M P C achieves the following:

II~kll~-+ minll~ll~

exponentially

as k ~ oo. (32)

< P r o o f > The above directly follows the convergence of observer-based unconstrained Q-ILC.(Lee et al., 1997) Q.E.D

3.3

Constrained Batch MPC

In industrial processes, process variables axe very often constrained in a way or others. The constraints may be imposed on the inputs, input rate (in terms of discrete-time index), and the outputs. Besides, in order for the linearized model in (7) to be a valid representation of the original nonlinear system, we may consider an additional constraint imposed on input change. The minimization criterion of (23) together with linear inequality constraints constitutes a quadratic programming just as in regular constrained MPC.

3.4

Implementation

The proposed batch MPC algorithm is implemented as follows: Starting with u0,/~010, and Go, let k = 0. At k + 1, s t e p 1 - M P C Advancing with t from 0 to N - 1, calculate AL/k+I according to (25) or by solving a

Numerical

Illustrations

We provide three numerical examples; two for linear SISO systems and one for a nonlinear batch reactor. In all the examples, process models (G) are derived from the zero-order hold equivalents of the continuous-time models with sampling period of h=l. E x a m p l e 1. C o n t i n u o u s vs. B a t c h M P C We assumed that the true and nominal processes are derived from G(s) = 1/15s 2 + 8s + 1 and G(s) = 0.8/12s 2 + 7s + 1, respectively. Figure 1 shows output responses of the blind combination of regular MPC with Q-ILC. As was discussed in subsection 2.3, the same tracking error results over the repeated batch runs. In Fig. 2, output responses of batch MPC w i t h p = m = 5 are shown. This time, the same large square-wave disturbance was assumed for w(t) in the repeated batches. In addition, time-correlated random noise with covarinace 0.22 was assumed for v(t) and also superimposed on w(t). The corresponding partition of the steady state Kalman filter was used for/C. From the figure, we can see that proposed batch MPC attains desired asymptotic zero tracking offset while completely rejecting the repeated disturbance and effectively attenuating the random noise through Kalman filtering. E x a m p l e 2. M H B - M P C o f U n s t a b l e P r o c e s s In this example, control of a linear unstable process was attempted. In the first case, we assumed that both the true and nominal process models are derived from G(s) = 1/15s 2 + 2s - 1 with no model error. Figure 3 shows results of batch MPC with p = m = 3. We can see that the output diverges as batch number increases. Next, we applied MHBMPC to the process. This time, the nominal model was derived from G(s) = 0.8/20s 2 + 3 s - 1 to demonstrate the robustness of MHB-MPC together. This model has the same unstable pole as the true process but different values for the steady state gain as well as stable pole. As can be seen from Fig. 4, MHB-MPC steers the unstable process to converge even when there is model error. Example 3. M H B - M P C of a Batch Reactor with Second-Order Exothermic Reaction As the lastexample, we considered temperature control of a nonlinear batch reactor where a secondorder exothermic reaction, A --+ B, takes place. It is assumed that the reactor has a cooling jacket whose temperatuer can be directlymanipulated. Dynamics of the reactor system is described by

dT UA Mcp(T-Tj)dt dCA s 2 dt = - koe- r r CA

AHV ~k°e

_ ~__a_ 2 RrCA(33 ) (34)

$878

PSE

'97-ESCAPE-7 Joint Conference

The following values are used for parameters:

UA/MCp = 0.0288(1/min) AHV/MCp = -5.79(°K • liter/mol) E/R = - 1 3 5 5 0 ( ° K ) k0 = 0.47(l/tool. min) T(0) = 298.15(°K) CA(O) = 0.9(mol/liter) To obtain a linear model for MHB-MPC design, we conducted two open-loop tests where the difference between the two inputs is kept constant. Figure 5 shows the test inputs and resulting outputs. Difference of the output responses was taken as the step response from which G was constructed. This rough model was used throughout the successive batch runs without further updating. In Fig. 6, results of MHB-MPC are shown. It can be observed that the reactor temperature starts to closely approach to the reference trajectory after first several transient runs while coping with nonlinearity as well as disturbances.

5

Conclusion

As a way to overcome the limitations of regular model predictive control(MPC) in batch process control, we have proposed a novel batch MPC algorithm utilizing the concept of quadratic criterion-based iterative learning control(Q-ILC). As a special case of batch MPC, we also proposed maximum horizon batch MPC(MHB-MPC) as a convergence guaranteeing algorithm. The proposed batch MPC algorithm carries out prediction based on not only the current measurement but also the previous batch information. Thanks to this feature, batch MPC retains the advantages of not only Q-ILC such as zero tracking offset even in case of model error and complete rejection of repeated disturbances but also regular MPCs such as real-time disturbance rejection and constraint handling as so forth. Through

a series of numerical study, the proposed algorithm was found to exhibit outstanding performance in batch process control. Through the present study, it is believed that we have set a cornerstone for batch process control study which is only in its infant stage. Further extension of the proposed algorithm to a more general, unified formulation for feedback-learning combination is being under investigation. Acknowledgement The authors would like to acknowledge LG Yonam Foundation(Korea), the Automation Research Center (Korea) at POSTECH, Korea Science and Engineering Foundation, and NSF (USA) for financial support.

References

[1]

Moore, K. L., Itertive Learning Control for Deterministic System, Springer-Verlager, NY, 1993.

[2] Muske, K. R. and Rawlings, J. B. ,"Stability of Constrained Receding Horizon Control", IEEE Trans. A.C., 38, 1512-1510, 1993.

[3] Lee, J. H., "Recent Advances in Model Predictive Control and Other Related Areas", Proc. CPC-V, Tahoe City, 1996.

[4] Lee, K. S., Bang, S. H, and Chang, K. S., "Feedback- assisted Iterative Learning Control based on Inverse Process Model", J. Process Control, 4, 77-89, 1994.

[5] Lee, K. S., Lee, J. H., and Kim, Won Cheol, "Model-based Learning Control with a Quadratic Criterion for Time-Varying Linear Systems" submitted to A UTOMATICA, 1997. OUTPUT

OUTPUT

9

9 8 7 6 5 4 3 2 1 0

-1

1'0 2b ~'o 4'o s'0 6'0 T'0 80 ]1ME

Fig. 1. Output response of the blind combination of MPC with Q-ILC (Example 1).

'JO

20

30

40 TIME

50

60

70

80

Fig. 2. Output response of batch MPC. (Example 1)

$879

PSE '97-ESCAPE-7 Joint Conference Output

OUTPUT

10

12 I£

10

20

30

40 50 ]]me

60

70

80

10

Fig. 3. B a t c h M P C o f an unstable process

20

30

Fig. 4. M H B - M P C

in Example 2.

40 50 TIME

60

70

o f an unstable process

in Example 2.

TEST INPUT

TE S T OUTPUT

31

32

30

31 TEST-INPUT #2

30

29

2g 28 TEST-INPUT #1

28

TEST-OUTPUT#1

27 27 26

26

25 24

25 0

50

100 TIME

150

200

24

0

50

100 TIME

150

200

60

80

Fig. 5. Step test for batch reactor identification(Example 3).

O UTP UT

INPUT

38 ull 35

36

u7 u4

yl 34

30

32

25

30 20 28 15

26 24

10 0

20

40 TIME

60

80

0

Fig. 6. MHB-MPC of the batch reactor in Example 3.

20

80

40 TIME