Ind. Eng. Chem. Res. 2000, 39, 3835-3843
3835
Response Surface Tuning Methods in Dynamic Matrix Control of a Pressure Tank System Aiping Jiang† and Arthur Jutan* Department of Chemical and Biochemical Engineering, University of Western Ontario, London, Ontario, Canada N6A 5B9
Dynamic matrix control (DMC) is implemented on a nonlinear, multivariable, experimental pressure tank. Linearized models for the pressure tank are identified using subspace methods, which prove to be suitable over a limited operating range. When the DMC algorithm is implemented, some tuning/design parameters are easy to select but others require trial and error. Here, a novel approach to the selection of these other parameters using response surface methodology is used. Finally, the DMC algorithm is implemented on-line and leads to good control. 1. Introduction Model predictive controllers (MPCs) refer to a class of algorithms that compute a sequence of manipulated variable adjustments to optimize the future behavior of a plant. They are different from some other advanced controllers in that a dynamic optimization problem is solved on-line at each control execution. This reduces the uncertainty effects caused by model mismatch, nonlinearity, and time variation and increases the control efficiency. The first description of MPC applications was presented by Richalet et al., at a 1976 Conference.1 Since then, MPC technology has developed considerably. During the past decade, there has been an increasing trend in the industry toward the use of MPCs especially for multivariable systems. A survey of MPC technology products by Qin and Badgwell2 has shown that there have been more than 2000 applications to date, and it was anticipated that the future of MPC technology is bright. This paper focuses on the application of DMC (dynamic matrix control, one subset algorithm of MPCs) to a multivariable pressure tank system. Conventional PID controllers were not able to satisfactorily control this process due to strong nonlinearities and strong interactions between control loops.3 This process is thus a good candidate for a multivariable DMC application. In the application of DMC, the one important step is the parameter tuning. It is a crucial and sometimes difficult step. Several tuning guidelines proposed for DMC have been suggested, such as those by Clarke et al.,4 Georgiou et al.,5 Maurath et al.,6 and Hinde and Cooper.7 But they are not universal guidelines, and for different processes, determining all the parameters requires a lot of time spent in trial and error. In this paper, we use response surface methodology (RSM) to determine the parameters that are difficult to chose. This leads to a more systematic approach to tuning and helps ensure that the tuning parameters lie within an optimum range. * To whom correspondence should be addressed. E-mail:
[email protected]. Tel.: (519)661-2111, ext. 88322. † Present address: Institute of Petrochemical Engineering, East China University of Science and Technology, Shanghai, China.
Model identification is another key factor necessary for DMC application. This paper uses subspace methods to obtain a satisfactory model of the multivariable pressure tank system. 2. DMC Algorithm DMC is a popular choice within MPC technology. It was proposed by Cutler and Ramaker8 based on a linear step response process model. The quadratic performance objective over a finite prediction horizon is used, and the future plant output behavior is specified by trying to follow the setpoint as closely as possible subject to a move suppression term. For a linear time invariant SISO system, if the step dynamic response coefficients are given by
{0, a1, a2, ..., aN} then the output of the system can be computed from its step response model {ai} as follows: N-1
y(k) )
ai∆u(k - i) + y0 ∑ i)1
(1)
where
∆u(k) ) u(k) - u(k - 1) y0 ) aNu(k - N) and y(k) is the value of the output y at sampling time k; ai is the ith sampled step response coefficient; u(k i) is the value of the input variable at sampling time k - i; y0 is the value of y when k ) 0; N is the truncation point of the step response model. A p-step predictor with m future control actions can be defined as8
Yp(k + j) ) A∆U(k + i) + Y0(k + j) + e(k)
(2)
where Yp is a vector representing the output of the predictor for (j ) 1, 2, ..., p); ∆U is a vector representing the incremental changes in the manipulated variable, for (i ) 0, 1, ..., m - 1); A is a p × m matrix consisting
10.1021/ie990311m CCC: $19.00 © 2000 American Chemical Society Published on Web 08/31/2000
3836
Ind. Eng. Chem. Res., Vol. 39, No. 10, 2000
Figure 1. Pressure tank setup.
of the step response coefficients; and e(k) is the difference between the real system output and model output. The DMC algorithm uses the quadratic performance objective to minimize the square of the deviation between the predicted output trajectory and the setpoint values at p future sampling periods, that is,
Jp ) [(Yp - Yr)TQTQ(Yp - Yr) + ∆UTRTR∆U] T
3. Process Description
T
) [(A∆U - E) Q Q(A∆U(A∆U - E) + ∆UTRTR∆U] (3) where
E ) Yr - Y0 - e and Yr is the setpoint trajectory; Q is the diagonal weighting matrix for the output components; R is the suppression factor that penalizes the objective function for changes in the input ∆U. The solution of the problem is
∆U ) (ATQTQA + R2I)-1ATQTQE ) KE
Parameter design is somewhat difficult for the DMC algorithm because of many parameters and the interaction between parameter choices. In this paper, we use response surface methodology (RSM)9 to help with the choice of the more difficult parameter selections, to obtain a more optimal overall parameter design.
(4)
where I is the identity matrix and K is the gain matrix (m × p). The extension of DMC to a multivariable system is obtained by partitioning the matrices. From above we see that the DMC algorithm uses a linear step response model, which is easy to obtain from transformation of other models, such as an impulse response or state space model. Although m control moves (∆U) can be calculated from eq 4, normally only the first one is applied, and then eq 4 is used again to update the control policy, to keep the predictions close to the actual values of the output. So only the first row of the controller gain matrix K needs to be computed and stored; also, it can be computed off-line to save on-line computation time. In the objective of the DMC controller algorithm, a penalty term R on the control move is used. This results in smaller input moves and a less aggressive output response, so that a degree of robustness to model error is provided. In controller design, beside step response coefficients, a number of design (or tuning) parameters have to be specified. The design parameters include four scalars (N, p, m, and sampling interval Ts) and two weighting matrices (Q and R). In the present study we choose Q ) I, which implies that all p predicted errors are weighted equally.
The system used here is a pressure tank (Figure 1) through which air flows from a regulated supply. Control valves are installed on both the inlet and the outlet of the tank. The pressure in the tank and the outlet flow rate are measured and transmitted to a Data Translation interface board, which connects to a computer. The controlled variables are pressure and flow rate. This is a two-input two-output system and the variation in both inputs will influence both controlled variables. From the process design, it can be seen that this is a strongly interactive system. In fact, a relative gain analysis (RGA) of the process shows that the RGA elements are all close to 0.5. This implies that a control scheme based on multiloop PID would be very difficult.3 This is also a nonlinear system. Not only is the relation between the controlled variables and valve opening nonlinear, but there is also hysteresis in the valves themselves. Identification under these conditions can be challenging. In Figure 2, we show the real-time system open-loop responses in which step signals are applied to each valve, respectively. We found that for three of the four responses, the time constants were 26 s (from outlet valve to pressure), 35 s (from inlet valve to flow rate), and 40 s (from inlet valve to pressure), respectively. However, for the outlet valve to flow response, the dynamics is more complicated and consists of an immediate response, followed by a slower, first-order type decrease to an equilibrium flow. Any identification scheme would need to seek some compromise to accommodate this variety of dynamics. 4. Model Identification The DMC algorithm requires a suitable step response model. It can be obtained from an impulse response model or a state-space model. If an impulse response model is used, the correlation analysis (or time series) method can be used to estimate the model. While this is common practice for single output systems, it is not always efficient for multiinput/ multioutput systems (MIMO).
Ind. Eng. Chem. Res., Vol. 39, No. 10, 2000 3837
Figure 2. Step response of real plant and identified system (solid line, real plant; dashed line, identified system).
For multivariable systems, an alternative approach is to identify directly a state-space model using the socalled “subspace” method. The algorithm, N4SID, proposed by Overschee and De Moor,10 is based on the subspace method. One advantage of this method is that instead of having to guess the orders of the transfer functions in each of the channels, as required by the time series methods, the N4SID algorithm requires the user to select a single parametersthe order of the system. For additional details see Franklin et al.11 A successful identification also requires that the collected measurement data contain significant information about the system.. We set up the data collection system using SIMULINK and a Data Translation interface board. This was accomplished with the help of a Real Time Toolbox available from Humusoft. A number of practical issues needed to be decided prior to the identification experiments. In particular, the selection of a suitable input signal must be chosen. We tried both PRBS (pseudo-random binary signal) and band-limited white noise. The second choice gave somewhat better results. The sample time selection is not trivial for this system because not all the channels for this multivariable system have similar dynamics. Ultimately, a compromise is necessary, with the sample time chosen to capture most of the dynamics of the quickest response. After a number of trial and error stages, we selected a sample time τ ) 1 s. Figure 2 shows one of our identification results where the models are compared to the collected data. The range of operation depicted by the step tests is quite large and exaggerates the model mismatch somewhat. Nevertheless, it is clear that while some channels are modeled quite well, others are less well modeled, over a full operating range. Ultimately, the DMC algorithm will need to correct for this model mismatch if good control is to be realized.
4.1. Identified State Model. The N4SID subspace algorithm was used to identify a state-space model of the system based on the collected data. A fourth-order (n ) 4) model was chosen because it gave the best result compared with the third- and fifthorder model. The model is identified in the process innovations form, where K is the model noise gain matrix and e(k) is a white noise sequence. The identified model is given by
x(k + 1) ) Ax(k) + Bu(k) + Ke(k) y(k) ) Cx(k) + Du(k) + e(k) where
[
(5)
0.7303 -0.3368 -0.4584 0.1903 0.1266 0.8455 -0.5272 -0.3315 A) 0.3119 -0.2158 -0.0821 -0.1729 0.0021 0.0194 0.0080 0.7523
C)
[
0.0179 0.7748 -0.1335 -0.0915 B) 0.0392 -0.6176 0.0619 -0.0002
[
]
0.5922 0.3634 0.7171 0.0101 0.0631 -0.3181 0.1185 -0.6039
[
D)
[ ] 0 0 0 0
-0.2980 -0.2838 0.3440 -0.2012 K) -0.0206 0.6319 -0.0462 -0.0405
]
]
]
This state-space model can also be easily transformed to a step-response model.
3838
Ind. Eng. Chem. Res., Vol. 39, No. 10, 2000
Figure 3. Control results with different values of N (model coefficients).
Model validation steps12 were also carried out and showed that, over a modest operating range (less that those shown in the step response plots), the models were adequate. The ultimate test is, of course, when the models are used in the feedback control scheme. 5. Design Parameters for DMC Control The available DMC publications often describe the control algorithm itself, but very few offer easily implemented and well-accepted guidelines for the selection of the design parameters, and furthermore most of the guidelines given in the open literature concern SISO systems. From the first section, we know there are four design/ tuning parameter scalars (N, p, m, Ts) and one matrix (R) that needs to be selected. Usually, the scalars N, p, and Ts are chosen as design parameters, and the scalar m and matrix R are left as tuning parameters.13,6 Because the value of the sampling interval affects all the choices of the parameters N, p, and m, it is the most important parameter and needs to be well chosen; otherwise, all further design choices may lead to poor control. In our system, the sampling interval was chosen as 1 s. This choice was based on the time constants of the valves, Ts, as well as the rule of thumb that requires Ts to be chosen as approximately 1/10 of the smallest time constant in the process. Parameter N is the number of coefficients that are used in the step response model. The model truncation time is expressed as the multiplication of Ts and N and will affect the value of Y0. To obtain efficient control, most of the dynamics of the system have to be included within the truncation time. Usually, this time should be approximately equal to 90% of the settling time. For considerably smaller N, insufficient information about the process is obtained, and the error becomes
larger (see Figure 3). For considerably larger N, more computational time and storage space will be required. For MIMO systems, this can cause the control quality to decline significantly. In Figure 3, we see the trend of the decline. In our system, when N is in the range of 80-100, the control quality remains approximately constant. Parameter p is the step length of the predictive horizon. A smaller p will speed up the response, but the stability and robustness will suffer. For larger p, a slower dynamic response and better stability will be obtained, but the computational time will increase, so there should be a tradeoff between the speed and stability. Usually, a larger p value is recommended to guarantee the stability, such as p ) N,13 p ) t60/Ts + t65/Ts - 1,5 and p ) (5τp + td)/Ts.7 But this will increase the computational time, and also the dynamics of the feedback system may become sluggish. We propose that the value of p should be chosen approximately equal to the system time constant because this ensures that the main dynamics of the system will have been included. In our system, when the value of p is located in the range of 20-40, there are no large differences in control quality. However, when p ) 10 or p ) 100, the control quality becomes worse, especially in the former case (see Figure 4). We are now left with the more difficult parameters m and R, which need to be determined. R is used to limit the excessive change in the input variable, and m is the number of steps for the control horizon. When the value of the R matrix elements are larger, stable control is more easily achieved, but the dynamic response becomes sluggish and control quality suffers. When the value of m is larger, the controlling scope of the system will increase, but the stability and robustness will decrease. Usually, we choose matrix R ) fI, where I is an identity matrix and f is a constant. So our tuning parameters
Ind. Eng. Chem. Res., Vol. 39, No. 10, 2000 3839
Figure 4. Control results with different values of P.
are m and f. Determining these tuning parameters is often a trial and error procedure as no general methodology exists to fine-tune a process. Furthermore, tuning constants are often left at “safe” values because it is not known how far the current controller settings are from an optimum, or perhaps an unstable condition.14 In this paper, we use a statistical technique (RSM) to determine the tuning parameters. 6. Response Surface Methodology (RSM) The most extensive applications of RSM are in the situations where several input, or design, variables potentially influence some performance measure or quality characteristic of the product or process. The philosophy behind RSM is that, within a given region, a response η can be described as a function of a set of input, or design, variables, x1, x2, ..., xk. It is assumed that these variables can be freely changed and they affect the response in the region of interest. If there is an optimum of the true response, η, in the chosen area, it can be approximated as a local function of the k design variables in a second-order form, k
η ) β0 +
∑ i)1
βixi +
∑ ∑ i
