Control of Integrating Processes Using Dynamic Matrix Control

0263±8762/98/$10.00+0.00 q Institution of Chemical Engineers Trans IChemE, Vol 76, Part A, May 1998

CONTROL OF INTEGRATING PROCESSES USING DYNAMIC MATRIX CONTROL* Y. P. GUPTA Department of Chemical Engineering, Dalhousie University, Halifax, Canada

D

ynamic Matrix Control (DMC) has been popular for the control of chemical and petroleum processes. These processes commonly include integrating process units, which produce a ramp change in the output for a step change in input. When DMC is used for control of integrating process units, a steady-state offset occurs for sustained load changes. This offset is not acceptable for certain applications. A simple modi® cation in the DMC algorithm that eliminates this offset is presented in this paper. The performance of the proposed algorithm is presented on SISO and MIMO example problems through simulations and experimental results. Keywords: level control; model predictive control; dynamic matrix control

INTRODUCTION

However, when the DMC algorithm is used for control of processes containing integrating process units, steady-state offsets occur for sustained load changes. For applications where the performance objective is to ensure that the controlled variables remain very close to their set points, these steady-state offsets are not acceptable. Lee et al.7 and Lundstrom et al.8 have presented a state-space based MPC algorithm that facilitates the handling of integrating process units, and addresses some other issues. The state-space model is expressed in terms of step-response parameters for systems of stable and/or integrating dynamics. A Riccati equation needs to be solved in case of integrating process units. In this paper, an alternative approach to eliminate the steady-state offsets that are encountered when dealing with integrating process units is presented. The proposed approach does not require the formulation of the MPC problem in state-space form and is much simpler. Because of this advantage, it can be implemented directly in the stepresponse formulation of the DMC algorithm, which is currently one of the most popular and widely used MPC algorithms in the process industry. Many industrialists have preferred the DMC formulation because it is simple, intuitive and allows the formation of the prediction vector in a natural way. The elimination of the steady-state offsets is demonstrated by presenting the performance of the proposed algorithm on SISO and MIMO example problems through simulations and experimental results.

The integrating process units are commonly present in chemical and petroleum processes that produce a ramp change in the output for a step change in the input. This results from the process unit’ s material or energy imbalance. Non-self regulatory level processes are typical examples of integrating process units. Control of these process units is important for the control of the overall process. Levels are usually controlled by adjusting an inlet or an outlet ¯ ow and two distinct control objectives are involved. For some applications the objective is to control the level tightly and the variations in the manipulated variable are not important. This objective is encountered in the control of level in chemical reactors and ¯ otation columns. In other applications the control objective is to minimize the variations in the manipulated variable while keeping the level within speci® ed limits. This objective is referred to as averaging level control and is encountered in the control of level in surge tanks and distillation columns. Several ef® cient control algorithms (Marlin1 , Campo and Morari2 , Cheung and Luyben3 , Shunta and Feherari4 ) are available to meet the above two objectives whenever these units can be considered as SISO systems. However, real-time process optimization usually requires the operation of the process at the intersection of constraints. To accomplish this, a constrained multivariable control scheme needs to be implemented. The control of the outputs of integrating process units needs to be considered along with the control of other process outputs and all of the constraints need to be considered simultaneously. Decomposition of the problem, where control of the integrating process units is considered separately, does not provide the optimal solution that is otherwise possible. Model Predictive Control (MPC) is best suited for control of constrained MIMO processes (Prett and Garcia5 ). Cutler6 suggested the use of Dynamic Matrix Control (DMC) for the control of integrating process units.

THE STEADY-STATE OFFSET The steady-state offset that results from the use of the DMC algorithm for sustained load changes can be determined from its transfer function. A transfer function of the DMC algorithm for open-loop stable processes is given by Gupta9 . For integrating process units, the predicted response due to past inputs is a straight line (ramp) and a different expression applies for the change in predicted output due to past inputs. This changes the transfer function.

* Based on a paper presented at the 46th Canadian Chemical Engineering Conference, Kingston, Ontario, Canada, October 1996.

465

466

GUPTA

Therefore, the applicable transfer function of the DMC algorithm is ® rst derived for a SISO case. The control move for the DMC algorithm is given by Cutler and Ramaker1 0 : D u(k) = Cec (k)

(1)

where the C matrix is calculated from the dynamic matrix (A) of the process. The size of the A matrix is Np ´ Nm . To suppress the control moves, the diagonal elements of the AT A matrix are increased by multiplying them by a move suppression factor f. Let us say that the matrix C in equation (1) includes the effect of the move suppression used. The ® rst control move that is implemented can be expressed as: D u(k) =

Np i=1

ci ec (k + i)

(2)

th

where ci represents the i element in the ® rst row of C matrix and ec (k + i) represents the error at the it h step that needs to be compensated. The error, ec (k + i), is given by, ec (k + i) = Desired value of Adjusted predicted value output at point ± of output at point (k + i) (k + i) due to past control moves = ySP ± yp (k + i) (3) where, yp (k + i) = y(k) + D y(k + i)

i = 1, 2, ..., Np

(4)

y(k) represents the measured value of output at the current control instant and D y(k+i ) represents the expected change in output at i steps ahead due to past control moves. D y(k+i ) is the residual effect of past control moves that has not already shown up in the measured value of output. To account for unmeasured disturbances and model mismatch, the predicted response in the DMC algorithm is adjusted by an amount equal to the difference between the actual output and the model prediction. Note that the expression for the predicted response given above includes this adjustment. Since the predicted response due to past inputs is a straight line, its slope depends only upon the input implemented at the previous control instant. The other past inputs have no effect. Therefore, D y(k+i ) is give by, (5) D y(k + i) = ia1 u(k - 1) i = 1, 2, ..., Np where u(k- 1) represents the input at the previous control instant as a deviation from the initial steady state, and a1 is the ® rst step response coef® cient. If the input has not changed over the last N control intervals, D y(k + i) equals zero in the case of open-loop stable processes. Note that this is not necessarily true for integrating units. By combining equations (2) through (5), D u(k) =

=

Np i=1

ci [ySP

- y(k) - ia1 u(k - 1)]

(6)

where b1 b2

= =

Np i =1

ci

(10)

ici

(11)

Np i =1

The control move can also be expressed as, (12) D u(k) ; u(k) - u(k - 1) Substituting for D u(k) in equation (9) and by taking the ztransform of both sides, gives the transfer function of the DMC controller, u(z)/ e(z) = b1 / [1 + (a1 b2 - 1)z- 1 ]

D(z) ;

(13)

Now the steady-state offset can be calculated that results in the case of sustained disturbances. Consider a liquid-level control system where the inlet ¯ ow (Fi ) is the load variable and the outlet ¯ ow (F0 ) is the manipulated variable. The open-loop transfer functions for level (h) are given by h(s)/ F 0 (s) ;

Gp (s) = - 1/ (As)

(14)

h(s)/ F 1 (s) ;

Gd (s) = 1/ (As)

(15)

where A is the cross-sectional area of the tank. Assuming the dynamics of the measuring sensor and ® nal control element to be negligible, for a zero-order hold, the pulse transfer function, HGp (z), and the discrete-time closed-loop response of the system to load changes, h(z), are as follows (Stephanopoulos1 1 ): HGp (z) = -(1/ A)Tz- 1 / (1 - z- 1 )

(16)

h(z) = Z[Gd (s)F i (s)]/ [1 + HGp (z)D(z)]

(17)

For a unit step change in load, Z[Gd (s)F i (s)] = Z[(1/ A)t ] = (1/ A)Tz- 1 / (1 - z- 1 )2

(18)

By substituting equations (13), (16) and (18) into equation (17), then, h(z) =

(1/ A)Tz- 1 [1 + (a1 b2 - 1)z- 1 ] (1 - z )[(1 - z- 1 )(1 + (a1 b2 - 1)z- 1 ) - (1/ A)Tb1 z- 1 ] (19) Using ® nal-value theorem, the steady-state offset for a unit step change in load is given by, -1

offset

=a b /b 1 2

1

(20)

where b1 and b2 are de® ned by equations (10) and (11). This expression for the steady-state offset was checked against closed-loop simulation results obtained on the system in Example Problem 3. Several runs were conducted using different values of the move suppression factor. In each case, the observed offset was the same as given by equation (20).

Np i =1

ci [e(k) - ia1 u(k - 1)]

= e(k)

Np i=1

(7) PROPOSED CONTROL ALGORITHM

Np

ci - a1 u(k - 1)

= b1 e(k) - a1 b2 u(k - 1)

i =1

ici

(8) (9)

The adjusted predicted response of outputs due to past control moves is needed for calculating the errors to be compensated at any control instant. For a SISO case, this is shown in equation (3). In the proposed algorithm, the Trans IChemE, Vol 76, Part A, May 1998

CONTROL OF INTEGRATING PROCESSES USING DYNAMIC MATRIX CONTROL

adjusted predicted response of outputs that exhibit a ramp for a step change in input is calculated in a different manner. The rest of the procedure for control calculation remains the same. Advantage is taken of the fact that the predicted response due to past inputs is a straight line passing through the output at current control instant. The slope of the predicted response is determined from the slope of the output trajectory between the current and the previous control instants. Note that this slope includes the effect of unmeasured disturbances and any model mismatch that may be present. The adjusted predicted response yp (k + i) at i steps ahead is thus calculated by using the following equation. yp (k + i) = y(k) + i[y(k) - y(k - 1)]

i = 1, 2, ..., Np

(21) Note that the above calculation is affected by the noise on the output signal. The effect of noise and its reduction is considered in Example Problem 1. The adjusted predicted response of outputs that are open-loop stable is determined in the usual way. Let us consider a three-input±three-output case where the ® rst output exhibits a ramp for a step change in input and the other two outputs are open-loop stable. In this case, the adjusted predicted response of the ® rst output due to past control moves will be calculated using equation (21). For the last two outputs, the adjusted predicted response will be calculated from the following equation, yp (k + i) = [y(k) + D y(k + i)]

i = 1, 2, ..., Np (22)

where, D y(k + i = D Ai+ 1 D u(k - 1) + D Ai+ 2 D u(k - 2) + ...

+ D Ai+ N - 1 D u(k - N + 1)

(23)

and a2,1,i+ j - a2,1, j a2,2,i+ k - a2,2, j a2,3,i+ j - a2,3, j a3,1,i+ j - a3,1, j a3,2,i+ j - a3,2, j a3,3,i+ j - a3,3, j

_i j ; D A +

(24) In equation (24), a2 ,1 ,j refers to the change in the 2nd output due to a unit step change in 1st input occurring j steps ahead. a3 ,2 , j refers to the change in the 3rd output due to a unit step change in 2nd input occurring j steps ahead, and so on. Note that whenever j > N, a2 ,1 , j equals a2 ,1 ,N , a3 ,2 , j equals a3 ,2 ,N and so on, because the open-loop response is assumed to settle in N control intervals. Now, let us derive the transfer function for the proposed algorithm for a SISO case. By substituting equation (21) into equation (3), the error ec (k + i) is given by, ec (k + i) = ySP - y(k) - i[y(k) - y(k - 1)]

= 1, 2, ..., N = e(k) - i[y(k) - y(k - 1)] i = 1, 2, ..., N i

p

(25)

p

(26)

By adding and subtracting yS P , the above equation can be written as, ec (k + i) = e(k) - i[ySP

- y(k - 1) - (ySP - y(k))](27) (28) = e(k) - i[e(k - 1) - e(k)]

Trans IChemE, Vol 76, Part A, May 1998

467

By substituting equation (25) into equation (2), D u(k) =

Np i=1

ci [e(k) - i(e(k - 1) - e(k))]

(29)

Combining equations (10), (11), (12) and (29) yields, u(k) - u(k - 1) = b1 e(k) - b2 [e(k - 1) - e(k)]

(30)

By taking z - transform of both sides, the transfer function of the proposed algorithm is: (b1 + b2 ) - b2 z- 1 (31) (1 - z- 1 Therefore, it can be seen that in SISO cases, the proposed modi® cation reduces the DMC controller to a PI controller and thus eliminates the offset. For MIMO situations, the proposed modi® cation is made within the MPC framework, as shown above for the three-input±three-output case. Therefore, all of the advantages that the MPC framework offers in dealing with a constrained multivariable optimization problem, e.g., considering all constraints collectively, the ability to use the desired weights on different outputs in the objective function, the handling of process interactions, are preserved.

D(z) ;

u(z)/ e(z) =

CONTRO L PERFORMANC E The performance of the proposed control algorithm is presented on three example problems. The load changes were considered to be `unmeasured,’ that is, their effect was not fed forward in the control calculations. In all of the ® gures, the changes considered occur at time equal to zero. The negative time region shows the initial steady-state values of the controlled and manipulated variables. The results obtained are compared with the results from the regular DMC algorithm. The values of the tuning parameters (move suppression, control horizon, prediction horizon, weights on outputs) were the same in the two algorithms. These values and other details are given in the Appendix. Example Problem 1 This example considers an isothermal CSTR where a ® rst-order chemical reaction A ! B takes place. The total mass balance and the component mass balance are described by the following differential equations.

= (F - F)/ (Area) / dt = F (C - C )/ V - kC

dh/ dt dC A

(32)

i

i

Af

A

A

(33)

A two-input-two-output system is considered where the level in the reactor and concentration of component A are controlled by manipulating outlet ¯ ow and inlet feed composition. The simulation results obtained for step changes in set points are shown in Figure 1. Since the responses of the two algorithms are essentially the same, the response curves lie on top of each other and therefore cannot be distinguished in the ® gure. This comparison shows that the proposed modi® cation does not change the performance of the DMC algorithm for set point changes. The simulation results obtained for a step increase of 0.2 litre/second in feed ¯ ow to the reactor are shown in Figure 2. In this case, the DMC algorithm results in a steady-state offset not only in

468

GUPTA

Figure 1. Comparison of control performance for step changes in set points.

Figure 3. Effect of noise on control performance.

level but in concentration also. The proposed algorithm brings both of the outputs back to their set points. The effect of noise on the output signal can now be considered. It is assumed that the noise is normally distributed and that approximately 99.75% of the measurements are within 6 1% of the true value of the output. To simulate the

effect of the noise, normally distributed random numbers with a mean of zero and standard deviation (r ) of 0.33% were generated and added to the outputs. These outputs were then used in the control calculations. The results obtained for a step increase of 0.2 litre/second in feed ¯ ow are shown in Figure 3. As expected, the proposed algorithm was affected more than the DMC algorithm by the noise. The effect of the noise can be reduced by using a analog ® lter and/or by reading the outputs a number of times and using the average value. The second approach is very easy to implement through the data acquisition software. At any control instant, the output channels can be read a number of times in a very short time and an average value for each of the outputs can be calculated. The averaging reduces the standard deviation in the output. To study the effect of this approach, the outputs were read in the following order, y1 , y2 , y1 , y2 , and so on. Each of the two outputs was read 50 times and the averaged values for the outputs were used in the control algorithms. The resulting responses for a step increase of 0.2 litre/second in feed ¯ ow are shown in Figure 4 and are fairly smooth. This simulation shows that the adverse effect of the noise can be reduced.

Example Problem 2 This example considers a nonisothermal CSTR where an exothermic ® rst-order chemical reaction A ! B takes place. The total mass balance, the component mass balance and the energy balance are described by the following differential equations.

= (F - F)/ (Area) / dt = F (C - C )/ V - kC

dh/ dt Figure 2. Comparison of control performance for a step increase in feed ¯ ow.

dC A

(34)

i

i

Af

A

A

(35)

Trans IChemE, Vol 76, Part A, May 1998

CONTROL OF INTEGRATING PROCESSES USING DYNAMIC MATRIX CONTROL

469

Figure 4. Reduction in the effect of noise on control performance.

= F (T - T )/ V + JkC - Q/ (q where J = (-D H)/ (q c ) k = k edT / dt

i

i

A

cp V )

p

E/ RT

0

(36) (37) (38)

A three-input-three-output system is considered where the level, concentration of component A and reactor temperature are controlled by manipulating outlet ¯ ow, inlet feed composition and the amount of heat removed by the coolant. Again, normally distributed random numbers with a mean of zero and standard deviation of 0.33% were generated and added to the outputs as noise. Each of the three outputs was read 100 times and the averaged values were used in the control algorithms. The simulation results obtained for a step increase of 0.1 litre/second in feed ¯ ow to the reactor are shown in Figure 5. The DMC algorithm again results in a steady-state offset not only in the level but also in the temperature. The offsets resulted even when noise was not added. The proposed algorithm brings all of the outputs back to their set points. The simulation results were also checked for step changes in set points. Again, the two algorithms performed in essentially the same way.

Figure 5. Comparison of control performance for a step increase in feed ¯ ow.

for a step decrease of about 2 litres/minute in feed to the column are shown in Figure 6. With the DMC algorithm, the level continuously drops for this sustained disturbance and results in a steady-state offset of 11.5 cm. With the proposed algorithm, the level is controlled and brought back to its set point. Step changes in set point were also considered. Both of the algorithms took the level to the new set points and their transient responses were close to each other.

Example Problem 3

OTHER COMMENTS

This example problem considers the control of water level in an actual column, 10 cm in diameter and 2.5 m in height. The outlet ¯ ow from the bottom of the column was manipulated by changing the speed of a pump. The column was interfaced with an IBM type personal computer where the control algorithms were executed. To reduce the effect of noise, the output was read 50 times at each control calculation and the averaged value of the output was used in the control algorithms. The experimental results obtained

(a) To avoid the steady-state offsets that result from the use of the DMC algorithm, the control of the outputs of integrating process units can be considered separately from the control of other outputs. However, this decomposition would restrict the determination of the true optimum of the optimization problem that needs to be solved at every control instant. The proposed algorithm allows the consideration of all inputs, outputs and constraints in one optimization problem.

Trans IChemE, Vol 76, Part A, May 1998

470

GUPTA b1, b2 ci e(k) ec(k+i) f h k (k) N

Figure 6. Comparison of control performance for a step decrease in feed ¯ ow.

(b) The constraints are not considered in this study because the proposed modi® cation only affects the prediction part of the algorithm. The optimization part remains the same. (c) The comparison made in the paper shows that when the regular DMC algorithm is used in MIMO situations, the offset occurs not only in the liquid-level but also in other outputs. The comparison also shows that the proposed modi® cation does not alter the performance of the DMC for set point changes. CONCLUSIONS The proposed modi® cation is very simple and can be implemented directly in the step-response formulation used in the DMC algorithm. For processes containing integrating process units, the modi® cation eliminates the steady-state offsets that are otherwise encountered for sustained load changes. APPENDIX: PARAMETERS Example 1: Area = 0.05 m2 N = 100 f = 1.00001 w1 = 1 Example 2: Area = 0.05 m2

k = 0.25 s ± 1 Np = 100 T=6 s w2 = 1

Nm = 2

3 D H = ±50 000 000 q = 1 000 000 g/m cal/kmole cp = 1 cal/(g 8 K) E/R = 1000 8 K k0 =4.872997584s ± 1 N = 30 Np = 30 Nm = 2 f = 1.005 T=3 s w1 = 1 w2 = 2.5 w3 = 1/400 Heat removal is plotted as Q/(q cp ) Example 3: Area = 0.0081 m2 N = 25 Np = 25 Nm = 2 f = 1.3 T=2 s

NOMENCLATUR E A a1

area of cross-section change in the output, one step ahead, due to a unit step change in the input

de® ned by equations (10) and (11) ith element in the ® rst row of C matrix difference in the set point and measured value of output error in the output at i steps ahead that needs to be compensated move suppression factor level rate constant denotes current control instant number of control intervals in which the open-loop response settles Np number of points on each output trajectory at which the error is minimized Nm number of control moves into the future including the current one T control interval u(k) manipulated variable wi weight on the ith controlled variable y(k) measured value of the controlled (output) variable at current control instant yp(k + i) adjusted predicted value of output at point (k+i ) due to past inputs ySP set point of the controlled variable D y(k+i) change in output at i steps ahead from current value due to past control moves boldface indicates a vector or matrix

REFERENCES 1. Marlin, T. E., 1995, Process Control, (McGraw-Hill, New York), p. 581±599. 2. Campo, P. and Morari, M.,1989, Model predictive optimal averaging level control, AIChE J, 35: 579±591. 3. Cheung, T. F. and Luyben, W., 1979, Liquid-level control in single tanks and cascades of tanks with proportional-only and proportionalintegral feedback controllers, IEC Fund, 18: 15±21. 4. Shunta, J. and Feherari, W., 1976, Non-linear control of liquid level, Instr Techn, 43±48. 5. Prett, D. M. and Garcia, C. E., 1988, Fundamental Process Control, (Butterworths-Heinemann, Boston, MA), p. 3. 6. Cutler, C. R., 1982, Dynamic matrix control of imbalanced systems, ISA Trans, 21: 1±6. 7. Lee, J. H., Morari, M., Garcia, C. E., 1994, State-space interpretation of model predictive control, Automatica, 30: 707±717. 8. Lundstrom, P., Lee, J. H., Morari, M. and Skogestad, S., 1995, Limitations of dynamic matrix control, Comput Chem Eng, 19: 409± 421. 9. Gupta, Y. P., 1993, Characteristic equations and robust stability of a simpli® ed predictive control algorithm, Can J Chem Eng, 71: 617± 624. 10. Cutler, C. R. and Ramaker, B. L., 1980, Dynamic matrix controlÐ A computer control algorithm, Proc JACC, Paper WP5-B. 11. Stephanopoulos, G., 1984, Chemical Process Control, (Prentice-Hall, Englewood, NJ), p. 617± 623.

ACKNOWLEDGEMENTS The ® nancial support provided by the Natural Sciences and Engineering Research Council of Canada is gratefully acknowledged The help of Ms M. Bhagi in obtaining the experimental results is appreciated.

ADDRESS Correspondence concerning this paper should be addressed to Professor Y. P. Gupta, Department of Chemical Engineering, DalTech, Dalhousie University, Halifax, Nova Scotia, Canada B3J 2X4 (E-mail: yash.gupta@ dal.ca). The manuscript was received 3 November 1997 and accepted for publication after revision 17 April 1998.

Trans IChemE, Vol 76, Part A, May 1998