2730
Ind. Eng. Chem. Res. 2004, 43, 2730-2735
Real-Time Application of Scheduling Quasi-Min-Max Model Predictive Control to a Bench-Scale Neutralization Reactor Yaohui Lu,† Yaman Arkun,*,‡ and Ahmet Palazogˇ lu§ School of Chemical Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332, College of Engineering, Koc¸ University, Rumelifeneri, Sariyer, Istanbul, 34450 Turkey, and Department of Chemical Engineering and Materials Science, University of California, One Shields Avenue, Davis, California 95616-5294
Scheduling quasi-min-max model predictive control is an MPC algorithm to control nonlinear systems. In this paper, real-time application of the scheduling quasi-min-max MPC algorithm to a bench-scale pH neutralization reactor is presented. The control performance is compared with multilinear model-based MPC and a scheduling IMC-PID controller. 1. Introduction
2. Scheduling Quasi-Min-Max MPC
Model predictive control (MPC), also known as moving or receding-horizon control, has been a very effective technology in chemical process industries where the presence of operational constraints creates a challenging control task.1 Basically, MPC solves online an open-loop constrained optimization problem at each time instant and implements only the first computed control move. The optimization is repeated at the next time step by updating the initial condition with the new state. In most of the applications of MPC, finite-horizon implementation of the algorithm has been reported. Yet, the finite-horizon formulation is difficult to analyze theoretically.2 Imposing terminal equality and inequality constraints as done by Chen and Allgo¨wer3 and Michalska and Mayne,4 respectively, is one common way to guarantee stability. For a nice summary of the state of the art in nonlinear MPC theory, the reader is referred to work by Mayne et al.5 Scheduling quasi-min-max MPC is an infinitehorizon MPC algorithm developed by Lu and Arkun6 for nonlinear systems. In this algorithm, the nonlinear system is expressed as a combination of a dynamic linear model with a linear parameter-varying (LPV) model. The linear dynamic model is used to express the current “local” dynamic behavior of the nonlinear system, and the LPV model is used to approximate the future nonlinear behavior. The closed-loop stability is guaranteed by the feasibility of the associated recedinghorizon optimization problem. The next section introduces briefly the quasi-minmax MPC formulation and is followed by the description of the pH neutralization problem. For the first time, the results of real-time application of quasi-min-max MPC are presented and the comparison of its performance to various benchmark control strategies is discussed. Finally, conclusions and future directions are presented.
It is assumed that the nonlinear first-principle model is known:
* To whom correspondence should be addressed. E-mail:
[email protected]. † Georgia Institute of Technology. Presently at Air Products & Chemicals Inc., Allentown, PA. ‡ Koc¸ University. § University of California.
x˜˘ ) f(x˜ ,u˜ ) y˜ ) Cx˜
(1)
where x˜ is the state variable, u˜ is the control variable, and y˜ is the output variable. By using a Taylor series expansion around the current point k, we will have
|
∂f x˜˘ ≈ f[x˜ (k|k),u(k-1)] + [x˜ - x˜ (k|k)] + ∂x˜ x˜ (k|k),u˜ (k-1) ∂f [u˜ - u˜ (k - 1)] (2) ∂u˜ x˜ (k|k),u˜ (k-1)
|
where x˜ (k|k) is the estimated state variable and u˜ (k-1) is the control action calculated from the previous point. If the deviation variables are defined around the desired setpoint (x˜ sp, u˜ sp, y˜ sp): x ) x˜ - x˜ sp, u ) u˜ - u˜ sp, and y ) y˜ - y˜ sp, then the continuous state-space model can be written as6
x˘ ≈ Ac(k) x + Bc(k) u + π(k) y ) Cx
(3)
where
π(k) ) f[x˜ (k|k),u˜ (k-1)] - Ac(k) [x˜ (k|k) - x˜ sp] Bc(k) [u˜ (k-1) - u˜ sp] and
| |
Ac(k) )
∂f ∂x˜ x˜ (k|k),u˜ (k-1)
Bc(k) )
∂f ∂u˜ x˜ (k|k),u˜ (k-1)
The discrete-time state-space model can be obtained using the sampling time ∆t
x(k+1|k) ≈ A(k) x(k|k) + B(k) u(k) + θ(k) where
10.1021/ie030311t CCC: $27.50 © 2004 American Chemical Society Published on Web 04/22/2004
(4)
Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004 2731
A(k) ) eA (k) ∆t
umin(k) e u(k) e umax(k)
B(k) ) Ac(k)-1(eA (k) ∆t - I)Bc(k)
ymin(k+1) e y(k+1|k) e ymax(k+1)
θ(k) ) -Ac(k)-1(I - eA (k) ∆t)π(k)
(ii) The objective function starting from the next step and extending to infinity J∞1 (k) is bounded from above by a worst-case value
c
c
c
While eq 3 explains the local dynamics, we assume that the future nonlinear behavior can be approximated by a LPV model (see work by Johansen and Foss7 and Banarjee et al.8). First of all, the plant operating space is partitioned into several local descriptions by linear models that are valid in some regimes. Then, a global model is constructed among the regions by using a parameter vector as an interpolating or a model validity function. This defines an LPV model of the form
x(k+i+1|k) ) A[F(k+i|k)] x(k+i|k) + B[F(k+i|k)] u(k+i) y(k+i|k) ) Cx(k+i|k)
J∞1 (k) e xT(k+1|k) P(k) x(k+1|k)
Φ(k) ) J10(k) + xT(k+1|k) P(k) x(k+1|k)
where
A[F(k+i|k)] )
Fj(k+i|k) Aj ∑ j)1 N
B[F(k+i|k)] )
Fj(k+i|k) Bj ∑ j)1
(6)
Aj and Bj are local models that can be obtained around different operating points. Here, N is the number of the local models included in the LPV model, and Fj(k+i|k) is the scheduling parameter reflecting the validity of the local linear models. The technique of global linearization suggests that a LPV approximation of a nonlinear system is possible.9 In summary, a linear model (eq 4) is obtained to express the current linear dynamics, and the future nonlinear behavior is assumed to vary within a region constructed by a LPV model (eq 5). Scheduling quasimin-max MPC minimizes an infinite-horizon objective function based on the combination of the linear model and the LPV model. The formulation of the algorithm can be expressed as
min U∞ 0
J∞0
∞
)
[x(k+i|k)TQx(k+i|k) + u(k+i)TRu(k+i)] ∑ i)0
) xT(k|k) Qx(k|k) + uT(k) Ru(k) + ∞
[x(k+i|k)TQx(k+i|k) + uT(k+i) Ru(k+i)] ∑ i)1 ) J10(k) + J∞1 (k)
The optimization can be solved by semidefinite programming.6 Lyapunov stability is guaranteed when the feasible solutions of the algorithm are implemented in a receding-horizon fashion. It should be noted that the above algorithm should have feasible solutions at all sampling times (i.e., solvable at all times), and feasibility at time k does not necessarily imply feasibility at time k + 1. To guard against such feasibility problems, the algorithm can be modified to include constraints on all of the future inputs and outputs (as was done by Lu and Arkun10). Feasibility problems were not encountered in the pH experimentation. 3. pH Neutralization Reactor, Experimental Setup, and Modeling The real-time application of the scheduling quasimin-max algorithm was conducted at the University of California at Davis by using a bench-scale pH neutralization experiment. An acid stream (HCl solution) and an alkaline stream (NaOH and NaHCO3 solution) are pumped to a well-mixed tank. The pH value is measured through a sensor located in the tank. The goal of the controller is to drive the system to different pH conditions. More details about the experimental apparatus can be found in ref 11. The process model based on component balances can be written as
1 1 x˜˘ ) (x˜ 1ini - x˜ 1) - x˜ 1u˜ τ τ
(7)
1 1 x˜˘ ) - x˜ 2 + (x˜ 2ini - x˜ 2)u˜ τ τ
where Q and R are positive-definite and U∞0 stands for all of the control actions from the current time to infinity.
U∞0
) {u(k+i), i ) 0, 1, 2, ...}
(8)
The optimization is solved subject to the following constraints: (i) Constraints on the control action that will be implemented on the plant u(k) and the resulting output y(k+1|k)
(11)
where
(5)
N
(10)
where P(k) is a positive-definite matrix that will be obtained from optimization. (iii) The Lyapunov stability constraint forces the objective function of quasi-min-max to decrease monotonically
Φ(k) < Φ(k-1) ig1
(9)
1 1 x˜˘ ) - x˜ 3 + (x˜ 3ini - x˜ 3)u˜ τ τ
(12)
where
τ)
qB V , u˜ ) qA qA
(13)
x˜ 1 is the concentration of HCl, x˜ 2 is the concentration of
2732 Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004
uncertainty are associated with the acid, base, and buffer concentrations in the storage tanks. As a result, experimental and model titration curves differ as noted in Figure 1. The first step is to build a state-space model of the form of eq 3. If we apply the first-order Taylor expansion to eq 12 around the current point [x˜ (k|k), u˜ (k1)], we will have
x˘ ) A(k) x + B(k) u + Π(k)
(16)
where x and u are deviation variables defined at the desired setpoint, x ) x˜ - x˜ sp(k), u ) u˜ - u˜ sp(k),
A(k) ) -[1 + u(k-1)]/τ 0 -[1 + u(k-1)]/τ 0 0 -[1 + u(k-1)]/τ 0 0
[
Figure 1. Steady-state curve of the pH value versus the input u˜ ) qB/qA. The continuous line is the model prediction, and the + symbols denote the experimental values. The circles indicate the points around which the local models are obtained.
NaOH, and x˜ 3 is the concentration of NaHCO3. In the experiment, these concentrations are not measured. V is the volume of the reactor, qA is the flow rate of the acid, and qB is the flow rate of the base. In the experiment, the acid flow is kept constant and the manipulated variable is the base flow rate. The values of the parameters are as follows:
and
[
][]
-(1/τ)x˜ 1(k|k) B1 B(k) ) (1/τ)[x˜ 2ini - x˜ 2(k|k)] ≡ B2 B3 (1/τ)[x˜ 3ini - x˜ 3(k|k)] and
Π(k) ) f(x˜ (k|k),u˜ (k-1)) - A(k) [x˜ (k|k) - x˜ sp(k)] B(k) [u˜ (k-1) - u˜ sp(k)] (17)
[ ]
Π1(k) ≡ Π2(k) Π3(k)
x˜ 1ini ) 0.0012 mol of HCl/L x˜ 2ini ) 0.002 mol of NaOH/L
From the knowledge of h(x˜ ,y˜ ) ) 0, we can have a specific function of η that
x˜ 3ini ) 0.0025 mol of NaHCO3/L qA ) 1 L/min
y˜ ) η(x˜ )
V ) 2.500 L The pH value is obtained through the following nonlinear relationship:
Kw h(x˜ ,y˜ ) ) ξ + x˜ 2 + x˜ 3 - x˜ 1 ξ
x˜ 3 )0 Kxξ 1+ Kw
(18)
and the first-order Taylor expansion can be used again to linearize the function η
][
x˜ 1 - x˜ 1sp ∂η ∂η ∂η y˜ - y˜ sp ) ∂x˜ ∂x˜ ∂x˜ x˜ 2 - x˜ 2sp 1 2 3 x˜ 3 - x˜ 3sp
[
(14)
y ) Cx ξ ) 10
(15)
Kx ) 10-7 mol/L Kw ) 10
2
(19)
2
mol /L
From the model, it is observed that the pH value exhibits a strong nonlinearity with respect to the input (u˜ ) qB/qA). The steady-state curve is depicted in Figure 1. The actual pH system deviates from the model because of the presence of unknown dissolved salts in the water used in the experiment. Other sources of
(20)
where
and y˜ is the pH value with
-14
]
or
where -y˜
]
∂η ) ∂x˜ i
∂h ∂x˜ i ξ ln(10)
∂h ∂ξ
i ) 1-3
(21)
The linearized state-space model (eqs 16 and 20) is not an observable system. The maximally unobservable subspace is two-dimensional. However, the system is detectable because the unobservable subspace is stable. Next, we construct an observable and controllable reduced-order system whose output is the pH value that we want to control. To accomplish this, the output
Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004 2733
equation (eq 20) is differentiated to get
y˘ ) )
∂η ∂η ∂η x˘ + x˘ + x˘ ∂x˜ 1 1 ∂x˜ 2 2 ∂x˜ 3 3 ∂η 1 1 - [1 + u˜ (k-1)]x1 - x˜ 1(k|k) u + Π1(k) + ∂x˜ 1 τ τ ∂η 1 1 - [1 + u˜ (k-1)]x2 + [x˜ 2ini - x˜ 2(k|k)]u + ∂x˜ 2 τ τ ∂η 1 1 - [1 + u˜ (k-1)]x3 + [x˜ 3ini Π2(k) + ∂x˜ 3 τ τ
[ [
]
]
[
Figure 2. Strategy to update the current linear model.
]
x˜ 3(k|k)]u + Π3(k)
[
]
∂η 1 ∂η ∂η ) - [1 + u˜ (k - 1)] x1 + x2 + x3 + τ ∂x˜ 1 ∂x˜ 2 ∂x3 3
3
∂η
∂η
Bi u + ∑ Πi(k) ∑ ˜ ˜ i)1 ∂x i)1∂x i
i
3 3 ∂η ∂η 1 ) - [1 + u˜ (k-1)]y + Bi u+ Πi(k) (22) τ ∂x˜ i ˜i i)1 i)1∂x
∑
∑
or
y˘ ) ap(k) y + bp(k) u + Ω(k)
(23)
where
1 ap ) - (1 + u˜ ) τ
3
bp )
∂η
Bi ∑ ∂x˜ i)1
u
i 3
Ω)
∂η
∑ ˜ i)1∂x
tion and the actual measurement. This updating of the current dynamics is critical to the success of the scheduling quasi-min-max algorithm as shown in section 4 and is reflected in the robustness and superior performance of the proposed control strategy.
Πi(k) (24)
i
Equation 23 is a state-space model, and the state variable is the pH value itself. The state-space model is time-varying and depends on the operating conditions. ap(k) is a function of u˜ (k-1|k-1), which is available. However, both Bi and ∂η/∂x˜ i are functions of the concentrations x˜ (k|k), which are not available from the plant. An updating strategy then is needed to get these concentrations by using the first-principle model as well as the plant measurements. We use a nonlinear observer to estimate the states. An open-loop model-based observer updates the two unobservable (but stable) states x˜ 2(t) and x˜ 3(t). Specifically, the states x˜ 2(t) and x˜ 3(t) are computed from the original nonlinear model (eq 12) in an open-loop fashion. The third observable state x˜ 1(t) is estimated using the nonlinear output equation (14) in the following way. First, ξ is computed from eq 15 using the pH measurement. Next, the states x˜ 2(t) and x˜ 3(t) are calculated from eq 12 and substituted into the output equation (14), and finally x˜ 1(t) is computed to satisfy eq 14. This procedure is depicted in Figure 2, in which we show the plant, the observer, and the controller. The updating strategy shown in Figure 2 addresses the mismatch between the theoretical model and the actual experimental system. As shown in eqs 23 and 24, state-space model parameters ap(k) and bp(k) are updated in real-time and cover the mismatch between the model gains and the experiment gains. On the other hand, updating of Ω(k) covers the measurement errors and any other discrepancies between the model predic-
4. Discussion of the Experimental Results In the experiments, three pH values are selected, 5, 7, and 9 (see Figure 1). Local models are obtained around these three operating points by using the strategy discussed in the previous section. The goal of the control algorithm is to track the variations in the pH value, as shown in Figures 3 and 4. The experiments may start at any initial pH value, while the setpoint is pH ) 7. Then the setpoint changes from pH ) 7 to 9, from pH ) 9 back to 7, from pH ) 7 to 5, and finally from pH ) 5 to 9. In addition to scheduling quasi-minmax MPC, two other scheduling control algorithms are also tested. These two controllers are the scheduling IMC-PID controller and the multilinear model-based MPC controller. Because of the presence of noise in the system, we did not use the derivative action. The algorithm of the scheduling IMC-PI controller can be expressed as
∫
u˜ ) uss + K1(y˜ - y˜ ss) + K2 [y˜ (t) - y˜ ss] dt
(25)
where 3
K1 )
∑ j)1
3
φjK1j
K2 )
φjK2j ∑ j)1
(26)
and the controller gains K1j and K2j are obtained from the local linear models. For the chosen three steadystate operating points, the state-space model (23) can be easily converted into the first-order model
y˜ j(s) )
k1j u˜ (s) k2js + 1 j
j ) 1-3
(27)
where k1j ) -bpj/apj and k2j ) -1/apj. Note that at steady state Ω(k) is equal to zero. On the basis of the wellknown tuning rules of IMC,12 the gain and the integral constant can be found as
K1j )
k2j λjk1j
K2j )
1 λjk1j
(28)
where λj is the tuning parameter. At high-sensitivity regions, pH ) 5 and 9, the best tuned values of λ are 100 s, and at the low-sensitivity region, pH ) 7, the best value is found to be 10 s. φj is the normalized Gaussian
2734 Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004
Figure 3. Comparison between the scheduling IMC-PID controller and the scheduling quasi-min-max MPC controller.
Figure 4. Comparison between the scheduling MPC controller and the scheduling quasi-min-max MPC controller.
function (see ref 13) and can be calculated from
[( ∑ [ ( exp -
φj )
3
exp -
i)1
)] )]
y˜ j(k|k) - y˜ (k) 2σj
2
y˜ i(k|k) - y˜ (k) 2σi
2
j ) 1-3 (29)
σj is the covariance of the measured signal and was selected as 0.25 in the controller tuning. The multilinear model-based MPC algorithm is modified from the algorithm by Kwon and Pearson,14 in which it was based on a single-linear time-varying model, and in this paper, it is designed based on the multilinear model. The formulation is as follows: N
min
J(k) )
N u(k+i|k)|i)1
[x(k+i|k)TQx(k+i|k) + ∑ i)1 u(k+i)TRu(k+i)] (30)
subject to umin(k) e u(k+i) e umax(k)
i ) 1, 2, ..., N
(31)
with the terminal constraint
x(k+N|k) ) 0
(32)
and 3
x(k+i+1|k) )
φj(k) [Ajx(k+i|k) + Bju(k+i)] ∑ j)1
i ) 1, 2, ..., N (33)
where the normalized weights are also obtained from eq 29. In the experimental runs, the control and prediction horizons were selected as N ) 10, which is considered to be long enough for the system dynamics. The state variable and input variable weights are the same as the weights used in the scheduling quasi-minmax MPC algorithm. A comparison of scheduling quasi-min-max MPC with scheduling IMC-PID is shown in Figure 3. The upper plot shows the setpoint tracking of the pH value, while the lower plot displays the manipulated variable moves. One can observe that scheduling quasi-min-max
MPC has a much better control performance; i.e., the pH value reaches the setpoint in a shorter time while the calculated control action is larger and quicker. It is notable that the control action is large when needed (during the initial transient). The response of the scheduling IMC-PI controller is sluggish especially around pH ) 5 and 9, where the pH value is highly sensitive to the changes in the manipulated variable. Tuning efforts were made, but the PI controller could not be tuned to give a better performance consistently across all of the transitions. The main reason is that the controller gains are obtained based on the theoretical model and there is no robustness feature to cover the mismatch between the model and the actual experiment. However, because of the more accurate estimation of current nonlinearities and a better prediction of the future dynamics, scheduling quasi-min-max MPC generates more aggressive control actions resulting in faster response times. The updating strategy shown in Figure 2 addresses the model/experiment mismatch directly and, consequently, improves the robustness and performance of the controller. It is also notable that the IMC-PI controller saturates at the lower limit when the setpoint is moved down to pH ) 5 (no antireset feature was used). This is not an issue for the MPC algorithm because constraints are optimally handled by the predictive control law. Remark. Scheduling quasi-min-max MPC does not have any tuning parameters except for the number of models used to describe the LPV representation of the nonlinear system. In all of the runs, the LPV model was constructed from the current updated model [A(k|k), b(k|k)] and the local model that corresponds to the new setpoint (final condition) [Af, Bf].6 Thus, the LPV model changes as the current model does, thereby accommodating not only noise and disturbance effects but also nonlinearities. A comparison of scheduling quasi-min-max MPC versus multilinear model-based MPC is depicted in Figure 4. Faster response and better tracking of the setpoint can be observed again for scheduling quasimin-max MPC. For the multilinear model-based MPC algorithm, even though multilinear models are considered and a model is obtained by interpolating among those models, this one model is used to predict all of the next N steps. The prediction based on this linear time-invariant model cannot account for the future
Ind. Eng. Chem. Res., Vol. 43, No. 11, 2004 2735
system dynamic changes even with a very long prediction horizon. Furthermore, any mismatch between the model and the experiment affects the performance. However, in the scheduling quasi-min-max MPC algorithm, the model contains two parts: the current linear model, which expresses the current behavior and addresses the model/experiment mismatch, and the LPV model, which explains the possible future nonlinear behaviors. Because of the better prediction and accurate expression of the current behavior, scheduling quasimin-max MPC can generate faster and larger control actions. This is especially apparent when the setpoint is moved from pH ) 7 to 5 and back; the control action generated by quasi-min-max MPC is appropriately aggressive and moves the valve to the limit, yet multilinear MPC is overly conservative and sends a rather cautious control signal. 5. Conclusions In this paper, a real-time application of the scheduling quasi-min-max MPC algorithm on a bench-scale pH neutralization reactor is discussed. Two other control algorithms are also tested for comparison: a scheduling IMC-PI controller and a multilinear model-based MPC with a terminal constraint. From the analysis of the experimental results, we observe that scheduling quasimin-max MPC has a better control performance because of its unique model-handling approach: a current linear model, which is updated online to capture the current dynamics and any mismatch between the model and the actual system, and a LPV model, which is used to cover the possible future nonlinear behaviors. Scheduling quasi-min-max MPC is shown to be a useful
method to control nonlinear processes that go through significant transitions. Literature Cited (1) Qin, J.; Badgwell, T. An overview of industrial model predictive control technology. In Fifth International Conference on Chemical Process Control CPC V; J. C. Kantor, C. E. G.; Carnahan, B., Eds.; 1996. (2) Bitmead, R. R.; Gevers, M.; Wertz, V. Adaptive optimal control: the thinking man’s GPC; Prentice Hall: Englewood Cliffs, NJ, 1990. (3) Chen, H.; Allg¨ower, F. Automatica 1998, 34, 1205-1217. (4) Michalska, H.; Mayne, D. Q. IEEE Trans. Autom. Control 1993, 38, 1623-1633. (5) Mayne, D. Q.; J. B.; Rawlings, C. V. R.; Scokaert, P. O. M. Automatica 2000, 36, 789-814. (6) Lu, Y.; Arkun, Y. J. Process Control 2002, 12, 589-604. (7) Johansen, T. A.; Foss, B. A. 12th IFAC World Congr., Sydney, Australia 1993, 1, 431-434. (8) Banerjee, A.; Arkun, Y.; Ogunnaike, B.; Pearson, R. AIChE J. 1997, 43, 1204-1226. (9) Boyd, S.; Ghaoui, L. E.; Feron, E.; Balakrishnan, V. Linear Matrix Inequalities in System and Control Theory; SIAM: Philadelphia, PA, 1994. (10) Lu, Y.; Arkun, Y. Automatica 2000, 36, 527-540. (11) Ga´lan, O.; Romagnoli, J. A.; Palazogˇlu, A. Chem. Eng. Sci. 2000, 55, 4435-4450. (12) Morari, M.; Zafiriou, E. Robust Process Control; Prentice Hall: Englewood Cliffs, NJ, 1989. (13) Brown, M. D.; Lightbody, D.; Irwin, G. W. IEE Proc.: Control Theory Appl. 1997, 144, 505-514. (14) Kwon, W. H.; Pearson, A. E. IEEE Trans. Autom. Control 1978, AC-23, 479-481.
Received for review April 14, 2003 Revised manuscript received March 5, 2004 Accepted March 18, 2004 IE030311T
