11482
Ind. Eng. Chem. Res. 2010, 49, 11482–11490
Design of a Robust Nonlinear Model Predictive Controller Based on a Hybrid Model and Comparison to Other Approaches Rosendo Dı´az-Mendoza and Hector Budman* Department of Chemical Engineering, UniVersity of Waterloo, Waterloo, Ontario, Canada N2L 3G1
A methodology to systematically design a model-based nonlinear model predictive controller is presented. The controller is referred to as hybrid since it uses the first-principles model to calculate the value of the controlled variables along the prediction and control horizons whereas uses the empirical model to ensure a terminal condition that accounts for model errors. The empirical Volterra series model was split into nominal and uncertain parts that were then used to formulate a structured singular value based robustness test. The proposed hybrid controller was compared against a robust empirical that uses solely an empirical model and to a nonrobust first principles model based nonlinear model predictive controller. To show the benefits of considering robustness in the controller formulation, extensive simulation studies were conducted that considered mismatch between the real process parameters and the model parameters. It is shown that in some case the performance of the hybrid controller can be superior to the purely empirical and to the first principles based controllers. 1. Introduction One of the most widely used control strategies in the process industry is model predictive control (MPC). An MPC controller minimizes at each sampling instant a cost function comprised of a main term that penalizes the deviations of the predicted outputs from its predefined set-points with an additional term that penalizes manipulated variables movements. At the end of every sampling instant, a new set of plant measurements is obtained and the optimization is repeated considering the new information. The optimization problem is solved, subject to constraints on input and output variables. To calculate the values of the predicted outputs that will be used in the cost function of the controller, it is required to use a model that describes the interrelationship between the manipulated and controlled variables. The original MPC work and most current industrial implementations use linear models for prediction. However, because most processes are nonlinear, the use of nonlinear models is expected that will result in a significant improvement of closed-loop performance.1 Accordingly, the current research is focusing on predictive algorithms based on nonlinear models. MPC algorithms that use nonlinear models to calculate the output predictions are referred to as nonlinear model predictive control (NMPC) algorithms. There are two types of models that have been used for calculating the predicted outputs in NMPC implementations: first principles and empirical. First-principles models are based on the mass, energy, and momentum balances of the process, whereas empirical models are directly calibrated from input-output data. Empirical models are easier to obtain than first-principles models since they can be calculated from a straightforward regression of input-output data. On the other hand, firstprinciples models are expected to provide better extrapolation accuracy than empirical models, implying that they are generally superior for predicting the process behavior at conditions that are different than the ones used for model calibration. Both types of models have been used in previous NMPC algorithms. First principles model based NMPC has been addressed by refs 2-4. * To whom correspondence should be addressed. Tel.: +1 519 888 4567 ext 36980. Fax: +1 519 746 4979. E-mail: hbudman@ cape.uwaterloo.ca.
Some types of empirical models that have been considered for NMPC design are polynomial autoregressive moving average models,5 Hammerstein models,6 Volterra series models,7-9 Volterra-Laguerre models,10 and Wiener models.11-13 Several studies have been conducted on nominal stability and on the use of different types of models for predictions. For example, the use of a terminal condition has been proposed to address the nominal stability of NMPC based on first-principle models.1 However, to the authors’ knowledge, the robustness of these algorithms to model error has not been thoroughly studied as yet and it has been identified as one of the key issues to be addressed for successful implementation.1 For any mathematical model a discrepancy between the process output and the model output will always exists. Ignoring this discrepancy could result in poor closed-loop performance, and in extreme cases closed-loop instability may occur. Thus, it is imperative to consider the effect of mismatch between the model and the process by designing a controller that will be robust with respect to this mismatch. The robustness of NMPC algorithms based on first-principles models have been only addressed through simulation studies14-18 since it is impossible to formulate a general robustness test. First principles model based controllers will be referred to in this work as FP-NMPC. On the other hand it was shown in ref 19 that certain nonlinear empirical models can be used to formulate a NMPC algorithm that can be systematically tuned for robustness with respect to model error. The robust controller proposed in ref 19 by the authors of the current study uses the structure of the Volterra series model to formulate a robust NMPC test as a µ-structured singular value20 (SSV) problem that can be used to calculate online the optimal control actions. Due to the use of an empirical model this controller is referred heretofore as robust empirical nonlinear model predictive controller (RENMPC). For this controller, at each sampling instant, the SSV-based test is calculated and used for bounding the worst case predictions and constraint violations along predefined prediction and control horizons in the presence of disturbances, set-point changes, and uncertainty in the parameters of the Volterra series model. The algorithm was shown to exhibit good robustness properties when
10.1021/ie1016283 2010 American Chemical Society Published on Web 10/08/2010
Ind. Eng. Chem. Res., Vol. 49, No. 22, 2010
operated around a fixed operating condition. On the other hand this previously reported RENMPC algorithm has not been compared as yet to FP-NMPC algorithms in terms of performance around different operating conditions and in terms of computation times. Following the above, two different goals are pursued in the current work. A first goal of the current study was to conduct a comparison to assess the relative advantages and disadvantages of the RENMPC with first principles model based (FP-NMPC). Then, motivated by the fact that the RENMPC requires long computational times, a second goal of the current research is the development of a new robust algorithm that can be computed faster. To that purpose, a third controller is proposed in this work that uses a first-principles model for short-term prediction while it uses an empirical model, a Volterra series, to enforce a robust terminal condition. Due to the combined use of the empirical and first-principles models, this third controller is referred heretofore as a hybrid nonlinear model predictive controller (HNMPC). This new HNMPC algorithm uses the first-principles model for calculating the values of the predicted outputs and the uncertain empirical model to propose a SSV-based robust terminal condition calculation. Since the predictions are calculated from the first-principles model, it is expected that the performance of the new controller will be superior to that of the RENMPC controller proposed in ref 19, while, at the same time, the algorithm ensures that the output is bounded at steady state by enforcing a terminal condition that accounts for model uncertainty. The additional obvious benefit of the hybrid controller is that it requires significantly shorter computation time than the RENMPC since it accounts for model uncertainty only in the terminal condition but does not consider uncertainty along the entire prediction horizon as RENMPC. To assess the benefits and trade-offs of the proposed methodology, the hybrid controller (HNMPC) is compared against the RENMPC controller proposed in ref 19 and FPNMPC in the presence of disturbances, set-point changes, and process parameter uncertainty. This comparison is intended to address the relative performance trade-offs between the expected superior extrapolating accuracy of the FP-NMPC and HNMPC versus the robustness properties of the RENMPC. At the same time the three controllers are compared with respect to computation times. It will be shown that for certain levels of model errors and in the presence of disturbances the HNMPC controller is able to achieve equivalent or better performance to the RENMPC and FP-NMPC controllers with a significant reduction in computation time as compared to the RENMPC. To assess the relative merits of the three controllers (HNMPC, RENMPC, and FP-NMPC), a case study involving a multipleinput-multiple-output (MIMO) chemical process was considered. This manuscript is organized as follows. In sections 2 and 3 a review of the Volterra series model based (RENMPC) and the first principles model based NMPC algorithms (FP-NMPC), respectively, is presented. Section 4 describes the methodology for designing a hybrid robust NMPC controller (HNMPC). Finally section 5 presents the case study, and conclusions are presented in section 6. 2. Robust Nonlinear Model Predictive Control Based on Volterra Series Models (RENMPC) This section briefly describes the robust nonlinear model predictive control algorithm proposed in ref 19. The algorithm is based on Volterra series models that have been shown to effectively capture nonlinear behavior.22 The specific motiva-
11483
tion to use this type of model for designing the robust NMPC is that it has a structure that can be split into two parts where the first part contains the nominal values of the Volterra series coefficients, whereas the second part contains the values of the uncertainty associated to each parameter. This feature is used to formulate an online mathematical robustness test. With the purpose of decreasing the number of parameters that must be identified the robust NMPC uses an autoregressive Volterra series model for which the relationship between the outputs and the inputs for a MIMO process is as follows: nARX
yˆχ(k) )
∑
M-1
∑h
hqχyˆχ(k - q) +
n(χ,1)u1(k
q)1 M-1 M-1
∑ ∑h
i,j(χ,1)u1(k
i)0
n)0
- i) u1(k - j) + ... +
j)i
M-1
∑h
- n) +
n)0
M-1 M-1 n(χ,1)unu(k
- n) +
∑ ∑h i)0
j)i
i,j(χ,nu)unu(k
- i) unu(k - j)
(1) where yˆ is the predicted output, ny is the number of outputs, u is the manipulated variable, nu is the number of inputs, nARX is the number of autoregressive terms, M is referred to as the memory of the system and it is generally chosen to correspond to the settling time of the process being modeled, and χ ) 1, ..., ny. The model parameters, i.e., the hn’s, n ) 0, ..., M - 1, for the linear terms and hi,j’s, i ) 0, ..., M - 1, j ) i, ..., M - 1, for the nonlinear terms, are obtained by least-squares regression or nonlinear optimization using process input-output data by imposing an appropriate input sequence. In ref 23 it was shown that to identify the coefficients for a system with polynomial degree N, it is necessary to use a N + 1 level pseudorandom multilevel sequence (PRMS). Autoregressive Volterra series models are simpler to use and have smaller sensitivity to noisy data since they require a smaller number of parameters as compared to their non-autoregressive counterpart. To provide for robustness when the process is expected to work around several operating regions, the process identification procedure must be conducted around these different operating regions. This can be accomplished by using several PRMS around different operating regions and identifying the Volterra series coefficients from all the collected data. Regardless of the quality of the data collected at different operating conditions and the values of M and nARX, the output of a Volterra series model will always be different from the actual process output due to model truncation and round-off errors. Thus, when using a Volterra series model for modelbased control, it is expected that performance will deteriorate due to the presence of this model error. The effect that the mismatch has on the controller’s performance can be considered by combining the nominal and uncertain parts of the model resulting in a family of models that represents the system to be controlled. The family of models is analyzed by the robust NMPC algorithm to obtain the model that results in the worst closed-loop performance. The basic premise used for designing the controller is that if the worst model in the set satisfies a specific performance index, then it can be assumed that all of the other models, within the family of models considered for robust design, will also satisfy that performance level. The analysis of the worst model that is associated with the worst value of the elements of the Y vector for the robust
11484
Ind. Eng. Chem. Res., Vol. 49, No. 22, 2010
[ ] [
os1 M11 M12 ) os2 M21 M22
][ ] is1 is2
is1 ) ∆os1
(5) (6)
The relationship between the exogenous signals os2 and is2 is given by os2 ) [M21∆[I - M11∆]-1M12 + M22]is2
Figure 1. Interconnection matrix representation.
[ ]
NMPC controller can be formulated in terms of a SSV test,24-26 where the structure of Y is as follows: ysp ˆ 1(k) 1 (k) - y l (k + p) yˆ1(k + p) ysp 1 l Y) ˆ ny(k) ysp ny (k) - y ysp ny (k
∑h
yˆχ(k) )
max ||Y||∞ g kssv S µ∆(M) g kssv
ˆ χ(k qχy
- q) +
where the maximization is carried out with respect to the uncertainty in the Volterra series coefficients defined by the following vectors:
l + p) - yˆny(k + p)
n(χ,1)
i, j(χ,1)
( δhn(χ,1))u1(k - n) +
( δhi, j(χ,1))u1(k - i) u1(k - j) + ... +
The robustness test proposed in (8) is used to obtain a bound on the worst deviation of |Y|∞ for the family of models defined by combining (3) and (9) in the presence of external disturbances. The test in (8) can be reformulated by the following constrained optimization problem:24
n(χ,nu)
( δhn(χ,nu))unu(k - n) +
n)0
M-1 M-1
i, j(χ,nu)
i)0
j)i
with
(10)
M-1
∑ (h
(9)
HiL ) [δh0(i,1), ..., δh0(i,nu), ..., δhM-1(i,1), ..., δhM-1(i,nu)] HiNL ) [δh0,0(i,1), ..., δh0,0(i,nu), ..., δhM-1,M-1(i,1), ..., δhM-1,M-1(i,nu)]
j)i
∑ ∑ (h
L L L HL ) [H1 , H2 , ..., Hny] NL NL NL HNL ) [H1 , H2 , ..., Hny ]
n)0
∑ ∑ (h i)0
∑ (h
(8)
(2)
M-1
q)1 M-1 M-1
where os2 is related to Y and is2 is related to wχ(k) and dχ(k). Details on the construction of this interconnection are given in ref 19. On the basis of these interconnected matrices, the worst value of the elements of the Y vector, i.e., the worst |Y|∞, can be calculated by the following SSV test:24 wrt HL,HNL
where p is the prediction horizon. Following the above, the first step for designing a robust NMPC consists of modifying (1) to include the coefficient uncertainty as follows: nARX
(7)
( δhi, j(χ,nu))unu(k - i) unu(k - j) + icχ(k) + wχ(k) + dχ(k)
(3)
where δhn, n ) 0, ..., M - 1, and δhi,j, i ) 0, ..., M - 1, j ) i, ..., M - 1, are the uncertainties in the coefficients of the linear and nonlinear Volterra series terms, respectively. icχ(k), χ ) 1, ..., ny, is the effect of the initial conditions, wχ(k), χ ) 1, ..., ny, is a feedback term that considers the effect of unmeasured disturbances, and dχ(k), χ ) 1, ..., ny, is the effect of measured disturbances. The feedback term wχ(k) is calculated according to the following equation: wχ(k) ) yreal ˆ χ(k - 1) χ (k - 1) - y
(4)
For the purpose of prediction, the feedback term in (4) is assumed to be equal for all time intervals along the prediction horizon, which is a common assumption for linear MPC implementations. The formulation of the structured singular value based robustness test requires an appropriate interconnection matrix M and accompanying uncertainty description ∆. The interconnection between M and ∆ is shown schematically in Figure 1; the inputs to M are is1 and is2, and the outputs are os1 and os2. If M is partitioned into a structure that is compatible with the structure of the uncertainty matrix ∆, the equations that describe the system of Figure 1 are
max ||Y||∞ )
wrt HL,HNL
max
wrt kssv s.t. µ∆(M)gkssv
(kssv)
(11)
The constrained optimization problem proposed in (11), also known as a skew µ problem, has been shown to be convex.24 As shown in ref 19, the skew µ problem can be extended to account for input and output constraints by appending the values that should be kept within bounds to the prediction vector Y. Following this rationale it is possible to include the following: (i) constraints on manipulated variable changes ∆U to prevent an excessive movement of the manipulated variables; (ii) constraints on the values of the manipulated variables ulimits to account for actuator limits, and (iii) a terminal value condition tc to ensure convergence at steady state. The vector with the appended variables related to manipulated variable changes’ values, manipulated variables value constraints, and the steadystate output prediction can then be used within the following optimization problem:
max
wrt HL,HNL
Y ∆U ulimits tc
||
g kssv S µ∆(M) g kssv
(12)
∞
Then, a skew µ problem is solved where a bound for the augmented vector of variables is solved as follows:
Ind. Eng. Chem. Res., Vol. 49, No. 22, 2010
max
wrt HL,HNL
Y ∆U ulimits tc
||
)
max
wrt kssv s.t. µ∆(M)gkssv
∞
(13)
(kssv)
where the interconnection matrix M has a structure similar to (5); M11 and M22 are zero matrices; M12 is a function of the scalar kssv and the elements of the U vector, i.e., M12 ) f(kssv,U); and M21 is a function of kssv and the elements of the vectors Y, ∆U, ulimits, and tc i.e., M21 ) f(kssv,Y,∆U,ulimits,tc). Complete details on how to formulate the interconnection matrix M, the prediction vector Y, and the vectors quantifying ∆U, ulimits, and tc to be used in eq 13 are provided in ref 19. The RENMPC control law is based on the evaluation of the worst cost function when uncertainty is considered in the Volterra series coefficients by means of the SSV-based test proposed in (13). The objective of the proposed controller is to minimize the value of the worst-case cost function with respect to the manipulated variables vector U at each sampling instant. Accordingly, the control calculation at each sampling instant consists of the solution of the following problem:
J1 ) min
( | |) max
wrt U wrt HL,HNL
Y ∆U ulimits tc
∞
) min wrt U
(
max
wrt kssv s.t. µ∆(M)gkssv
)
(kssv)
ˆ J2 Y ∆UJ2 J2 ) min wrt U ulimits J2 tcJ2
| |
(18) ∞
where ∆UJ2 accounts for manipulated variables movement weighting, ulimitsJ2 accounts for manipulated variables constraints and tcJ2 accounts for terminal condition constraints. The vector ˆ J2 that considers the deviation of the predicted outputs from Y its predefined set-points can be written as
[
]
ysp - yˆ(k) ˆ J2 ) l Y ) sp y - yˆ(k + p - 1) ysp - (g(xˆ(k), u(k)) + fb(k)) l sp y - (g(xˆ(k + p - 1), u(k + m - 1)) + fb(k))
[
]
(19)
where fb is a feedback term that is calculated as follows: fb(k) ) yprocess(k - 1) - yˆ(k - 1)
(20)
The equation used to observe the system states has the following form:
U ) [u1(k), ..., u1(k + m), ..., unu(k), ..., unu(k + m)]
(15) and m is the control horizon. The controller formulation given by (14) has the disadvantage that it does not have integral action in the presence of model error as shown in ref 19. To correct for this situation when the error is within a region specified by ε′, i.e., when |E′|∞ < ε′, then the prediction for the calculation of the worst case is done without considering the uncertainty resulting in zero offset.19 This proposed correction results in a two-mode controller where, for |E′|∞ g ε′, the Volterra model with uncertainty is used to calculate the worst case whereas for |E′|∞ < ε′ the Volterra model without uncertainty is used for that calculation. 3. Nonlinear Model Predictive Control Based on First-Principles Model (FP-NMPC) This section describes the FP-NMPC algorithm used for comparison with the RENMPC and with the HNMPC algorithms. In this case the equations that describe the interrelationship between the inputs, outputs, and states can be written in discrete form as follows:
(16)
where f and g are nonlinear functions obtained from the firstprinciples model, x are the system states, u are the system inputs, and yˆ are the predicted system outputs. If some or all system states are not measured, an observer can be used to estimate their values. The corresponding observer equations are as follows: xˆ(k + 1) ) f(xˆ(k), u(k)) yˆ(k) ) g(xˆ(k), u(k))
where xˆ denotes and observed state. To ensure a fair comparison between the controllers the structure of the cost functions is similar to (13). Accordingly, the cost function can be written as
(14)
where the vector U has the following structure:
x(k + 1) ) f(x(k), u(k)) yˆ(k) ) g(x(k), u(k))
11485
(17)
ˆ ) X
[
][
]
xˆ(k + 1) f(xˆ(k), u(k)) l l ) xˆ(k + p - 1) f(xˆ(k + p - 2), u(k + m - 1)) (21)
In this work, to isolate the effect of the controller from the effect of the observer, the initial conditions of the estimated states were assumed to be equal to the actual states. To ensure a fair comparison, the structure of ∆UJ2 is similar to the one used in ref 19, and the terms ulimitsJ2 and tcJ2 were formulated using a SSV problem as follows: max
|
ulimits J2 tcJ2
|
) ∞
max
wrt kres s.t. µ∆r(Mr)gkres
(kres)
(22)
The structures of ∆r and Mr are provided as Supporting Information. Therefore the cost function of the controller in (18) can be rewritten as ˆ J2 Y J2 ) min ∆UJ2 wrt U kres
| |
(23) ∞
where kres is obtained from (22).
4. Hybrid First Principles Empirical Model Based Robust NMPC (HNMPC) This section presents the methodology to design a HNMPC. The algorithm uses the first-principles model to calculate the output predictions along the predefined prediction horizon, whereas the uncertain empirical model is used to propose a robust stability test. Accordingly the cost function of the controller is
11486
Ind. Eng. Chem. Res., Vol. 49, No. 22, 2010
ˆ J3 Y ∆UJ3 J3 ) min wrt U ulimits J3 tcunc
| |
A
(24)
dh ) q 1 + q2 + q3 - q 4 dt
(29)
where q4 is calculated as
∞
ˆ J3 ) Y ˆ J2, ∆UJ3 ) ∆UJ2, and ulimitsJ3 ) ulimitsJ2. The term where Y tcunc is quantified using the uncertain Volterra series model of (3) to calculate the value of yˆ(k+p), i.e., the value of the predicted outputs at the end of the prediction horizon, as follows:
q4 ) CV(h)n
(30)
The equations that describe the dynamics of the effluent reaction invariants (Wa4,Wb4) can be obtained from a mass balance on each of the ionic species as follows:
nARX
∑h
yˆχ(k + p) )
+ p - q) +
ˆ χ(k qχy
q)1 M-1
∑ (h
n(χ,1)
i)0
i, j(χ,1)
( δhi, j(χ,1))u12(k + m) + ... +
j)i M-1
∑ (h
n(χ,nu)
n)0
dWa4 ) q1(Wa1 - Wa4) + q2(Wa2 - Wa4) + dt q3(Wa3 - Wa4) (31)
Ah
dW b4 ) q1(Wb1 - Wb4) + q2(Wb2 - Wb4) + dt q3(Wb3 - Wb4) (32)
( δhn(χ,1))u1(k + m) +
n)0 M-1 M-1
∑ ∑ (h
Ah
Finally, the pH can be obtained from Wa and Wb by using the following nonlinear relationship:
( δhn(χ,nu))unu(k + m) +
M-1 M-1
∑ ∑ (h i)0
j)i
i, j(χ,nu)
( δhi, j(χ,nu))unu2(k + m) + icχ(k + p) + wχ(k + p) + dχ(k + p)
(25)
To quantify the terms ulimitsJ3 and tcunc, a SSV test similar to (22) was formulated as follows: max
|
ulimits J3 tcunc
|
) ∞
max
wrt krun s.t. µ∆runc(Mrunc)gkrun
(krun)
(26)
The structures of ∆runc and Mrunc are provided as Supporting Information. Thus, the only difference between the FP-NMPC described in the previous section and the HNMPC described in this section is that the constraint in (22) does not consider model uncertainty, whereas in (26) uncertainty is accounted for through the uncertainties in the Volterra model coefficients. Because the FP-NMPC controller in (18) does not consider uncertainty for calculating tc, it is expected that its performance will be inferior to the performance obtained with the HNMPC in those cases where there is a large mismatch between the first-principles model and the real process. Since tcunc in (26) is obtained from the uncertain empirical model, the robustness test will account ˆ J3, and the control for the effect of the mismatch on the vector Y law will minimize the mismatch effect. Therefore, it is expected that the performance of the HNMPC controller will surpass that of the FP-NMPC controller for the same terminal region ε.
Wa + 10pH-14 + Wb
1 + 2 × 10pH-pK2 - 10-pH ) 0 1 + 10pK1-pH + 10pH-pK2 (33)
The control problem consists of maintaining the tank height h and the outlet stream pH at their predefined set-points by manipulating the flows of the acid q1 and base q3 streams, in the presence of perturbations in the flow of the buffer stream q2. The values of the operating conditions for the process are shown in Table 1. Both the RENMPC and HNMPC controllers required the identification of a Volterra series model that describes the interrelationship between the manipulated variables (q1 and q3) and the controlled variables (h and pH). An autoregressive Volterra series model with M ) 2 and nARX ) 1 was used for the RENMPC. These values were found by trial and error to result in a reasonable fit between the process and the model. The sampling time used for identification and simulation purposes was 25 s. To prevent numerical difficulties during the
5. Case Study The system studied involves a pH neutralization process27 schematically shown in Figure 2, where an acid (q1), a base (q3), and a buffer (q2) are mixed in a tank. Assuming perfect mixing, constant density, and complete solubility of the ions involved, the chemical equilibrium can be modeled according to the following two reaction invariants for each stream i ∈ [1,4]: Wai ) [H+]i - [OH-]i - [HCO3-]i - 2[CO32-]i
(27)
Wbi ) [H2CO3]i + [HCO3-]i + [CO32-]i
(28)
The overall mass balance on the tank yields
Figure 2. pH neutralization system.
Ind. Eng. Chem. Res., Vol. 49, No. 22, 2010 nseq
Table 1. Process Operating Conditions A pK1 pK2 q1
) ) ) )
VCL )
Wa1 ) 3 × 10-3 Wa2 ) -3 × 10-2 Wa3 ) -3.05 × 10-3
2
207cm 6.35 10.25 3 × 10-3 M HNO3
q3 ) 3 × 10-3 M NaOH +
NL
VC
Vd
y1 (height) y2 (pH) q1 (acid flow) q3 (base flow)
14 7 q1ss q3ss
5 3 0.2q1ss 0.225q3ss
HL ) 2
identification process, the manipulated and controlled variables were normalized as follows: Vn )
V - Vp Vd
(34)
where Vn is the normalized value of the variable V, the values of Vp and Vd can be found in Table 2. To obtain the nominal and uncertain values of the autoregressive Volterra series coefficients rseq, different PRMS sequences were applied to the system around different operating conditions and a set of Volterra series coefficients were obtained by regression. This set of coefficients was divided into two subsets corresponding to the coefficients of the linear and nonlinear terms. Correspondingly, VC defines the set of Volterra series coefficients as follows: NL VC ) [VCL1 , ..., VCLny, VCNL 1 , ..., VCny ]
(35)
VCLi ) [h0(i,1), ..., h0(i,nu), ..., hM-1(i,1), ..., hM-1(i,nu)] VCNL ) [h0,0(i,1), ..., h0,0(i,nu), ..., hM-1,M-1(i,1), ..., hM-1,M-1(i,nu)] i (36) for the jth PRMS sequence the identified Volterra series coefficients are L L VCidL ) [VCid1 , ..., VCidnseq] NL NL VCidNL ) [VCid1 , ..., VCidnseq]
(37)
where the elements in eq 37 are also vectors defined as follows: L L VCidLj ) [VCidj,1, ..., VCidj,ny] NL NL VCidNL ) [VCidj,1 , ..., VCidj,ny] j
)
(40)
where nseq is the total number of PRMS sequences used (nseq ) rseq × number of operating regions). The elements of the vectors VCL and VCNL are equal to the means of [VCidiL, ..., VCidnseqL] and [VCidiNL, ..., VCidnseqNL], respectively. The associated uncertainty matrix used for the construction of the interconnection matrix M was calculated on the basis of 2 standard deviations as follows:
Table 2. Values for normalization Vp
∑
VCidNL j nseq
j)1
5 × 10-5 M NaHCO3
variable
VCidLj nseq
j)1 nseq
Wb2 ) 3 × 10-2 Wb3 ) 5 × 10-5
q2 ) 3 × 10-2 M NaHCO3
∑
11487
(38)
L VCidj,i ) [h0(i,1), j, ..., h0(i,nu), j, ..., hM-1(i,1), j, ..., hM-1(i,nu), j] NL VCidj,i ) [h0,0(i,1), j, ..., h0,0(i,nu), j, ..., hM-1,M-1(i,1), j, ..., hM-1,M-1(i,nu), j]
(39) The nominal values of the coefficients that define the nominal model and that were used for the construction of the interconnection matrix M were calculated as follows:
NL
H
) 2
1 nseq 1 nseq
nseq
∑ (VCid
L j
- VCL)2
j)1
(41)
nseq
∑
(VCidNL j
- VC )
NL 2
j)1
Two Volterra series were identified to model the process. The first describes the relationship between the tank height and q1 and q3, i.e., h ) h(q1,q3); the second set describes the relationship between the pH and q1 and q3, i.e., pH ) pH(q1,q3). The FP-NMPC and HNMPC controllers use the ordinary differential equations (ODE’s) and algebraic equations that describe the behavior of the process (29-33). Since it was assumed that Wa and Wb could not be measured, an observer21 was used to estimate its values. Simulations were conducted to test the controllers’ abilities to reject disturbances as well as to track set-point changes. A disturbance profile consisting of a superposition of random changes in q2, shown schematically in Figure 3, was used for comparing the controllers, whereas the sequence of set-point changes used for the comparisons is shown schematically in Figure 4. The three controllers were compared with and without model error. The current comparison study considers uncertainty in the two valve parameters: CV and n. For the purpose of discussion, the values of these parameters that are used for the first-principles model will be referred heretofore as CVmodel and nmodel, whereas the actual values occurring in the process will be referred to as CVprocess and nprocess. Although the empirical Volterra series model does not explicitly use CVprocess and nprocess because they are assumed to be unknown, the nominal coefficient values of this model together with the uncertainties associated with these coefficients are identified from input-output data collected from experiments performed around different combinations of the values of CV and n. Thus, the family of models described by the Volterra series model with uncertainty is expected to capture the behavior of the process in the presence of variations in the values of CV and n. A key obstacle for conducting a fair comparison between the controllers is that the robust controllers (RENMPC and HNMPC) are designed for a worst case scenario, whereas the system is not necessarily operated always at worst case conditions. Since such a worst case cannot be found analytically, the approach adopted for conducting the comparative study is to test the closed-loop operation at different combinations of CV and n and then assess the performance on an average basis for all of these combinations. In that way an average performance for different realizations of the uncertain process parameters can be quantified. The key difference between FP-NMPC and RENMPC is that while the FP-NMPC used the same values of CVmodel and
11488
Ind. Eng. Chem. Res., Vol. 49, No. 22, 2010 Table 3. IAE Results for Case 1 [IAEFP-NMPC, IAEHNMPC, IAERENMPC] [0.15, 0.10]
n ) 0.45
n ) 0.50
n ) 0.55
CV ) 9.25 CV ) 8.75 CV ) 8.25
[2.47, 2.47, 4.90] [10.06, 4.15, 5.02] [3.41, 3.77, 5.48]
[3.22, 3.22, 4.05] [3.77, 3.77, 4.40] [4.00, 2.97, 4.41]
[4.33, 3.47, 5.57] [3.10, 2.23, 3.68] [2.37, 2.53, 4.26]
Table 4. IAE Results for Case 2 [IAEFP-NMPC, IAEHNMPC, IAERENMPC] [0.25, 0.10]
n ) 0.45
n ) 0.50
n ) 0.55
CV ) 9.25 CV ) 8.75 CV ) 8.25
[5.16, 4.72, 5.50] [4.64, 3.91, 6.01] [3.94, 4.65, 6.47]
[3.14, 3.14, 4.62] [3.44, 3.44, 4.97] [2.41, 3.01, 5.35]
[3.02, 2.85, 4.62] [3.77, 4.13, 4.45] [4.61, 4.75, 4.67]
Table 5. IAE Results for Case 3 [IAEFP-NMPC, IAEHNMPC, IAERENMPC]
Figure 3. Disturbance profile for case study.
[0.30, 0.1]
n ) 0.45
n ) 0.50
n ) 0.55
CV ) 9.25 CV ) 8.75 CV ) 8.25
[3.56, 3.24, 5.75] [4.35, 4.21, 5.86] [5.34, 5.18, 6.96]
[5.97, 3.68, 5.12] [4.87, 3.58, 6.10] [4.00, 3.54, 6.16]
[3.40, 3.62, 5.14] [3.47, 4.07, 4.49] [2.60, 3.67, 5.28]
Table 6. IAE Results for Case 4 [IAEFP-NMPC, IAEHNMPC, IAERENMPC] [0.10, 0.15]
n ) 0.45
n ) 0.50
n ) 0.55
CV ) 9.25 CV ) 8.75 CV ) 8.25
[2.73, 5.39, 5.07] [2.94, 2.84, 5.18] [6.18, 3.22, 6.10]
[3.59, 2.61, 4.18] [3.68, 4.17, 4.56] [2.40, 4.04, 5.01]
[4.08, 4.08, 4.29] [3.94, 3.94, 3.61] [3.05, 3.07, 4.32]
Table 7. IAE Results for Case 5 [IAEFP-NMPC, IAEHNMPC, IAERENMPC] [0.10, 0.25]
n ) 0.45
n ) 0.50
n ) 0.55
CV ) 9.25 CV ) 8.75 CV ) 8.25
[3.86, 3.04, 5.60] [3.85, 3.14, 6.28] [4.02, 4.21, 7.00]
[2.87, 5.69, 5.13] [3.80, 2.67, 5.64] [3.58, 2.53, 5.57]
[3.61, 3.42, 4.75] [2.65, 2.62, 4.30] [3.29, 3.32, 5.13]
Table 8. IAE Results for Case 6 [IAEFP-NMPC, IAEHNMPC, IAERENMPC] [0.10, 0.30]
n ) 0.45
n ) 0.50
n ) 0.55
CV ) 9.25 CV ) 8.75 CV ) 8.25
[5.23, 5.13, 5.96] [4.01, 3.52, 6.67] [5.30, 4.64, 7.29]
[3.29, 3.85, 5.76] [3.14, 3.14, 5.92] [3.66, 3.36, 5.71]
[5.20, 3.88, 4.22] [3.95, 3.06, 5.16] [4.44, 5.90, 5.39]
Table 9. Average IAE for Case Studies 1-6 Figure 4. pH set-point changes profile for case study.
nmodel for all of the combinations of CVprocess and nprocess tested, RENMPC used a model which was identified on the basis of the data obtained around different combinations of these parameters. On the other hand the hybrid controller (HNMPC) uses the values of CVmodel and nmodel for short-term prediction with the ODE’s while it uses the uncertain Volterra model identified around different combinations of CVprocess and nprocess for calculating the terminal condition given by (26). The objective of the comparison of the controllers is to test whether the variations in CV and n justify the use of the robust controllers (HNMPC and RENMPC) as compared to the first-principles controller (FP-NMPC), where the latter does not take into account these variations but it is expected to provide better prediction accuracy in the absence of these parameters’ variations. The performance of the controllers was assessed by calculating the integral of the absolute error (IAE). The comparisons were conducted for different values of the manipulated variables weights which determine the speed of the closed-loop response. Six different combinations of weights [Wi∆u1,Wi∆u2], i ) 1, ..., m, were used as follows: case 1 [0.15,0.1], case 2 [0.25,0.1], case 3 [0.3,0.1], case 4 [0.1,0.15], case 5 [0.1,0.25], and case 6 [0.1,0.3]. For all of the cases, the prediction horizon p was set to 10 and the control horizon m was set to 2. FP-NMPC and HNMPC were designed with the following parameters’ values: CVmodel ) 8.75 and nmodel ) 0.5. The three controllers were designed with a terminal condition region of tc ) 0.5. For each
case 1 case 2 case 3 case 4 case 5 case 6 average
IAEFP-NMPC
IAEHNMPC
IAERENMPC
4.08 3.79 4.17 3.62 3.50 4.24 3.90
3.17 3.84 3.86 3.70 3.40 4.05 3.67
4.64 5.18 5.65 4.70 5.49 5.79 5.24
case study nine different combinations of CVprocess and nprocess were simulated. Table 3-Table 8 show the values of the IAE for all nine combinations for FP-NMPC, HNMPC, and RENMPC. Table 9 shows the average IAE for the three controllers, and in the last row of that table an overall average of the six case studies is calculated for each controller. The first observation that can be made from Table 3 to Table 8 is that, for all combinations of weights, whenever the simulations are conducted without model error (center value in each table), FP-NMPC and HNMPC outperform RENMPC in terms of IAE. As mentioned above, this result is expected since, in the absence of model error, the predictions used in FP-NMPC and HNMPC are more accurate. It should be mentioned that for the nominal case the IAE for FP-NMPC and HNMPC are not always identical since HNMPC uses a terminal condition that is calculated with the Volterra model which is in slight error with respect to the ODE’s model. The improvements in IAE for FP-NMPC over RENMPC range from 23% for case study 4 (IAEFP-NMPC ) 3.68 versus IAERENMPC ) 4.56 in Table 6) to 88% for case study 6 (IAEFP-NMPC ) 3.14 versus IAERENMPC
Ind. Eng. Chem. Res., Vol. 49, No. 22, 2010
Figure 5. Controlled variables profiles for case study 4: CV ) 8.75; n ) 0.5 (solid line, FP-NMPC; dashed line, RENMPC).
Figure 6. Controlled variables profiles for case study 6: CV ) 8.75; n ) 0.5 (solid line, FP-NMPC; dashed line, RENMPC).
) 5.92 in Table 8). The evolution of the controlled variables with respect to time for these two cases is shown in Figure 5 and Figure 6, respectively. The improvements in IAE for HNMPC over RENMPC range from 9% for case study 4 (IAEHNMPC ) 4.17 versus IAERENMPC ) 4.56 in Table 6) to over 95% for case study 5 (IAEHNMPC ) 2.67 versus IAERENMPC ) 5.64 in Table 7). When model error is present, the IAE results also show some clear trends. An average IAE for each case study was calculated to enable a more systematic comparison between the three controllers. These calculated average values are presented in Table 9. The results from Table 9 indicate that for four out of the six cases that were considered (case studies 1, 3, 5, and 6) the IAE values of HNMPC are lower, than those of FP-NMPC, whereas the other two remaining cases (case studies 2 and 4) result in very similar IAE for both controllers. On the basis of the average IAE, the biggest improvement is of 28% corresponding to case study 1. The improvement obtained for case study 1 can be explained by the fact that the manipulated variable weights Wi∆u1 and Wi∆u2 used for this case were both very small, indicating that more aggressive control results in higher sensitivity to model error resulting in lower performance of FP-NMPC that was not designed for robustness. For case study 1 that used a combination of small manipulated variables’
11489
weights it was found that there is a relatively high number of CV and n combinations for which HNMPC is similar or it significantly outperforms FP-NMPC, e.g., seven cases out of nine in Table 3 (Wi∆u1 ) 0.15 and Wi∆u2 ) 0.1). Thus, the robust hybrid controller can be a better alternative than the nonrobust controller, especially in those cases when the weights penalizing the control actions are small. Table 9 also indicates that when the controller weights increase, the performance of the RENMPC further degrades. For example as Wi∆u1 increases from 0.1 to 0.3, the IAE consistently increases for RENMPC from 4.64 to 5.65 (Table 9). A similar trend is observed when Wi∆u2 increases, resulting in a corresponding increase of the IAE for RNMPC from 4.70 to 5.79 (Table 9). The same trend can be seen for HNMPC: when Wi∆u1 increases from 0.1 to 0.3, the HNMPC IAE’s consistently increase from 3.17 to 3.86 (Table 9), whereas when Wi∆u2 increases, the HNMPC IAE’s increase from 3.70 to 4.05 (Table 9). This trend can be explained by the fact that when the weights are larger, slower control actions result. Because the HNMPC uses a first-principles model for prediction, it is expected that its performance will be similar to the FP-NMPC. When the comparison also considers the hybrid controller, it can be seen that, on the basis of the IAE, HNMPC outperforms FP-NMPC and RENMPC. This trend can be explained by the fact that HNMPC only considers that the prediction at sampling instant k + p has uncertainty; therefore it is a less conservative controller than RENMPC that considers uncertainty in all the predictions from sampling instant k to k + p but it is more robust than FP-NMPC that does not consider uncertainty. The largest improvement of HNMPC over FPNMPC corresponds to a case in Table 3 (CV ) 8.75 and n ) 0.45) indicating that it is very important to account for uncertainty in the calculation of the terminal condition as is done with the HNMPC controller. Also, the overall averages given in the last row of Table 9 confirm that HNMPC can be a much better alternative than either the nonrobust first principles model based controller or the robust controller that considers uncertainty along all of the entire prediction horizon. However, the most important advantage of HNMPC over RENMPC is the significant reduction in computation time resulting from the fact that the SSV calculation is done for a terminal condition in HNMPC, whereas it is done for the entire prediction horizon for RENMPC. In general, to improve the speed of FP-NMPC and HNMPC, the set of ODE’s was discretized, coded in C, and used as a MEX file in Matlab. It was found that the computation time of HNMPC is 3 times longer than for FP-NMPC, but it is 166 times faster than that of RENMPC. Thus, the proposed HNMPC results in good performance, as shown by the results in Table 9, and at the same time it can be calculated much faster than with RENMPC and at comparable speed to FP-NMPC. Therefore, it can be concluded that the hybrid controller is an efficient nonlinear predictive control algorithm from the points of view of both performance and computational effort. The time required could be further decreased by coding the SSV calculation in C or Fortran, but this has been left for future studies. 6. Conclusions A hybrid first principles empirical model based robust NMPC algorithm (HNMPC) was presented and compared against a first principles model based nonrobust NMPC (FP-NMPC) and an empirical model based robust NMPC (RENMPC). The HNMPC and RENMPC controllers exploited the structure of the Volterra series model to account for model error by SSV calculations.
11490
Ind. Eng. Chem. Res., Vol. 49, No. 22, 2010
On the basis of the case studies, it was found that the hybrid controller generally outperforms the other two controllers in terms of closed-loop performance. It is both less conservative than RENMPC because it does not consider uncertainty along the entire prediction horizon as in RENMPC and it is more robust than the FP-NMPC because it enforces a robust terminal condition whereas the FP-NMPC does not account for model error. Apart from the performance improvements, it was shown that a key advantage of the HNMPC over the RENMPC technique reported in a previous work is a significant reduction in calculations by almost 2 orders of magnitude. Also, this drastic reduction in computational effort provides the freedom to increase the prediction horizon if necessary without increasing the dimensions of the interconnection matrix, thus maintaining the same computational demand. It may be that for other processes and for very small manipulated variables’ weights, where the short-term prediction may be more sensitive to model error than the case presented in this study, it could be important to consider uncertainty along the entire prediction horizon as done in the RENMPC algorithm. However, it should be noticed that the enforcement of the terminal condition in the three algorithms inherently leads to detuned controllers, thus making it less critical to account for model uncertainty along the entire horizon as done in the RENMPC method. Although not shown in this paper for brevity, when the terminal condition is removed, the RENMPC controller significantly outperformed the nonrobust FP-NMPC. However, the terminal condition is essential in order to ensure steady-state convergence. Also, a possible disadvantage of the hybrid controller is that it requires a firstprinciples model for short-term prediction. However, if such a model is available, HNMPC was shown to be an efficient and robust nonlinear predictive control algorithm. Acknowledgment We acknowledge NSERC (National Science and Engineering Research Council of Canada) for financial support. Supporting Information Available: Structure of the matrices Mr, ∆r, Mrunc, and ∆runc. This material is available free of charge via the Internet at http://pubs.acs.org. Literature Cited (1) Findeisen, R.; Allgo¨wer, F. An introduction to nonlinear model predictive control. Proceedings of the 21st Benelux Meeting on Systems and Control, Veldhoven, The Netherlands; 2002. (2) Wright, G. T.; Edgar, T. F. Nonlinear model predictive control of a fixed-bed water-gas shift reactor: An experimental study. Comput. Chem. Eng. 1994, 18, 83. (3) Santos, L. O.; Afonso, P. A. F. N. A.; Castro, J. A. A. M.; Oliveira, N. M. C.; Biegler, L. T. On-line implementation of nonlinear MPC: An experimental case study. Control Eng. Pract. 2001, 9, 847. (4) Nygaard, G.; Nævdal, G. Nonlinear model predictive control scheme for stabilizing annulus pressure during oil well drilling. J. Process Control. 2006, 16, 719. (5) He´rnandez, E.; Arkun, Y. Control of nonlinear systems using polynomial ARMA models. AIChE J. 1993, 39, 446.
(6) Fruzzetti, K. P.; Palazog˘lu, A.; McDonald, K. A. Nonlinear model predictive control using Hammerstein models. J. Process Control. 1997, 7, 31. (7) Doyle III, F. J.; Ogunnaike, B. A.; Pearson, R. K. Nonlinear model based control using second-order Volterra models. Automatica 1995, 31, 697. (8) Maner, B. R.; Doyle, F. J., III; Ogunnaike, B. A.; Pearson, R. K. Nonlinear model predictive control of a simulated multivariable polymerization reactor using second-order Volterra models. Automatica 1996, 32, 1285. (9) Maner, B. R.; Doyle III, F. J. Polymerization reactor control using autoregressive-plus Volterra-based MPC. AIChE J. 1997, 43, 1763. (10) Parker, R. S.; Doyle, F. J., III. Optimal control for a continuous bioreactor using an empirical nonlinear model. Ind. Eng. Chem. Res. 2001, 40, 1939. (11) Norquay, S. J.; Palazoglu, A.; Romagnoli, J. A. Application of Wiener model predictive control (WMPC) to a pH neutralization experiment. IEEE Trans. Control Syst. Technol. 1999, 7, 437. (12) Juricˇic´, G. D.; Strmcˇnik, S.; Matko, D. Wiener model based nonlinear model predictive control. Int. J. Syst. Sci. 2000, 31, 189. (13) Jeong, B.-G.; Yoo, K.-Y.; Rhee, H.-K. Nonlinear model predictive control using a Wiener model of a continuous methyl methacrylate polymerization reactor. Ind. Eng. Chem. Res. 2001, 40, 5968. (14) Nagy, Z. K.; Braatz, R. D. Robust nonlinear model predictive control of batch processes. AIChE J. 2003, 49, 1776. (15) Kawohl, M.; Heine, T.; King, R. A new approach for robust model predictive control of biological production processes. Chem. Eng. Sci. 2007, 62, 5212. (16) MagniL. ScattoliniR. Robustness and robust design of MPC for nonlinear discrete-time systems. In Assessment and Future Directions of Nonlinear Model PredictiVe Control, 1st ed.; Findeisen, R., Allgo¨wer, F., Biegler, L. T., Eds.; Springer: Berlin, Heildelberg, Germany, 2007; pp 239254. (17) Diehl, M.; Gerhard, J.; Marquardt, W.; Mo¨nnigmann, M. Numerical solution approaches for robust nonlinear optimal control problems. Comput. Chem. Eng. 2008, 32, 1279. (18) Zavala, V. M.; Biegler, L. T. The advanced-step NMPC controller: Optimality, stability and robustness. Automatica 2009, 45, 86. (19) Dı´az-Mendoza, R.; Budman, H. Structured singular valued based robust nonlinear model predictive controller using Volterra series models. J. Process Control. 2010, 5, 653. (20) Doyle, J. Analysis of feedback systems with structured uncertainties. IEE Proc.sD: Control Theory Appl. 1982, 129, 242. (21) Henson, M. A.; Seborg, D. E. Adaptive nonlinear control of a pH neutralization process. IEEE Trans. Control Syst. Technol. 1994, 3, 169. (22) Parker, R. S.; Heemstra, D.; Doyle, F. J., III; Pearson, R.; Ogunnaike, B. A. The identification of nonlinear models for process control using tailored “plant-friendly” input sequences. J. Process Control. 2001, 11, 237. (23) Nowak, R. D.; Van Been, B. D. Random and pseudorandom inputs for Volterra identification. IEEE Trans. Signal Process. 1994, 42, 2124. (24) Braatz, R. D.; Young, P. M.; Doyle, J. C.; Morari, M. Computational complexity of µ calculation. IEEE Trans. Autom. Control. 1994, 39, 1000. (25) Ma, D. L.; Braatz, R. D. Worst-case analysis of finite-time control policies. IEEE Trans. Control Syst. Technol. 2001, 9, 766. (26) Nagy, Z. K.; Braatz, R. D. Worst-case and distribution analysis of finite-time control trajectories for nonlinear distributed parameter systems. IEEE Trans. Control Syst. Technol. 2003, 11, 697. (27) Nahas, E. P.; Henson, M. A.; Seborg, D. E. Nonlinear internal model control strategy for neural network models. Comput. Chem. Eng. 1992, 16, 1039.
ReceiVed for reView May 14, 2010 Accepted September 14, 2010 IE1016283
