- Email: [email protected]
Sensitivity analysis of LP-MPC cascade control systems
Journal of Process Control 19 (2009) 16–24
Contents lists available at ScienceDirect
Journal of Process Control journal homepage: www.elsevier.com/locate/jprocont
Sensitivity analysis of LP-MPC cascade control systems Alexei Nikandrov, Christopher L.E. Swartz * Department of Chemical Engineering, McMaster University, 1280 Main Street West, Hamilton, Ontario, Canada L8S 4L7
a r t i c l e
i n f o
Article history: Received 26 June 2007 Received in revised form 13 February 2008 Accepted 15 February 2008
Keywords: Linear programming Model predictive control LP-MPC Sensitivity
a b s t r a c t Model predictive control (MPC) has found wide application in the chemical process industry as well as other industrial sectors. Commercial MPC systems are typically implemented in conjunction with a steady-state linear or quadratic programming optimizer, whose key functions are to track the economic optimum and to provide feasible set-points to the model predictive controller. The two-level system is complementary to real-time optimization which typically utilizes more complex models and is executed less frequently. Despite the widespread adoption of LP-MPC systems, occurrences of poor performance have been reported, where large variations in the computed set-points were observed. In this paper, we analyze the sensitivity of the LP solution to variation in the LP model bias, through which feedback to the LP layer occurs. We consider both multi-input, single-output (MISO) and multi-input, multi-output (MIMO) systems. Principles are illustrated through graphical representation as well as case studies. The performance of the two-level LP-MPC closed-loop system is evaluated and explained using results of the LP sensitivity analysis. Ó 2008 Elsevier Ltd. All rights reserved.
1. Introduction Model predictive control (MPC) is arguably the advanced control algorithm of choice in the chemical process industry, and has made inroads into other industrial sectors as well [10,13]. A dynamic model is utilized within the control algorithm to predict the effect of future plant inputs on the controlled outputs. Future inputs are computed in accordance with a performance objective, typically as the solution of an optimization problem. The inputs corresponding to the first control interval are implemented, and the calculation process repeated at the end of the sampling period, with the model predictions adjusted using the difference between the measured and predicted outputs. Details of the algorithm may be found, inter alia, in [4,9,13]. Industrial MPC systems are generally implemented in conjunction with a linear programming (LP) or quadratic programming (QP) steady-state optimizer [5,13,17,19,20]. The LP (or QP) typically uses a static model consistent with the dynamic MPC model, and is executed at the same frequency as the model predictive controller. The plant economic optimum may shift due to disturbances; thus the steady-state LP (QP) provides a bridge between a higher-level and less frequently executed real-time optimization (RTO) layer and the model predictive controller by making setpoint adjustments in response to changing conditions between RTO executions. The LP formulation may involve minimization of the deviation between the set-points and target values determined * Corresponding author. Tel.: +1 905 525 9140; fax: +1 905 521 1350. E-mail address: [email protected] (C.L.E. Swartz). 0959-1524/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.jprocont.2008.02.005
by the real-time optimizer, or optimization of an economic criterion directly. Fig. 1 illustrates the location of the steady-state LP (QP) within a plant automation hierarchy, as given in [19]. We note that the process would typically include local PID-type controllers and that several variants are possible, such as the presence of a plant unit optimization layer between the LP and plant-wide RTO layers [5,13]. Despite the apparent success of two-level LP-MPC systems, instances of poor performance have been reported [6,16]. Shah et al.[16], in the context of control performance monitoring, describe an industrial application in which the variation in the set-points exceeds that of the corresponding controlled variables. Kozub[6] also reports set-points being noisy relative to their controlled variables, and excessive variation in the set of inputs which are at their constraints. This motivates an investigation into the potential causes of such behavior. Ying and Joseph [19] provide stability theorems for LP-MPC and QP-MPC cascade control systems with no plant/model mismatch. Consideration of model uncertainty is included in a case study based on the Shell Standard Control Problem [12]. Kassmann et al.[5] present a formulation for robust steady-state target calculation. For elliptic uncertainty on the model parameters, the target calculation problem takes the form of a second order cone program (SOCP) which the authors solve using software based on a primaldual interior point algorithm. Steady-state targets are also computed in [11,14]. However, this is driven primarily not by economics, but rather to provide steady-state values for the plant inputs and states for inclusion in an MPC formulation in which
A. Nikandrov, C.L.E. Swartz / Journal of Process Control 19 (2009) 16–24
17
^ j kÞ ¼ yðkÞ C x ^ ðk j k 1Þ dðk ^ ^ dðk þ i j kÞ ¼ dðk j kÞ; i ¼ 1; . . . ; P where ^ ðk j kÞ ¼ x ^ðk j k 1Þ ¼ Ax ^ðk 1 j k 1Þ þ Buðk 1Þ: x The above estimation scheme is based on an assumption of steplike disturbances, and may result in poor control performance in the presence of measurement noise and certain other disturbance types. A more general state estimation framework that addresses these limitations is described in [7–9,15]. In the application studies that follow, we consider only the original DMC disturbance estimation scheme. 2.2. LP The set-points are determined through the solution of an LP of the form Fig. 1. LP (QP) within process automation hierarchy [19].
deviations of the states and inputs from corresponding steadystates are weighted by positive definite matrices. Feedback in the two-level LP-MPC configuration occurs through a bias term in the steady-state LP model. In this paper, we analyze the sensitivity of the LP solution to variations in the LP model bias. We consider first the LP level separately for multi-input, singleoutput (MISO) and multi-input, multi-output (MIMO) systems. This analysis is coupled with graphical representations to provide insights into the sensitivity effects. Thereafter, the performance of a two-level LP-MPC cascade system is evaluated, with observed performance related to earlier sensitivity results. 2. Problem definition The control structure we consider in this work is a two-layer LPMPC system in which the LP provides output set-points and/or input target values to the model predictive controller. The MPC-controlled system feeds back information to the LP layer, from which model adjustments may be determined. We present in this section the MPC and LP formulations used, and also describe the LP model bias update method.
minu ;y s:t:
p X
i þ ai y
m X
i bi u
i¼1
i¼1 ss
þd y¼G u y
min
6y6y
ð2Þ
max
umin 6 u 6 umax where Gss is the plant model gain matrix and d is a model bias that is updated in accordance with the actual plant outputs. We use the overbar to identify variables involved in the LP optimization. Com and y become the MPC set-points (includponents of the optimal u ing input target values). We typically assign a number of set-points equal to the number of manipulated inputs; specification of an excess number of set-points would in general result in offset. In real-time optimization, the bias update, d, is typically calculated as the difference between the measured outputs and those resulting from applying the plant inputs to the steady-state model [2]. In the present framework, this would correspond to d ¼ y Gss u; where y and u correspond to the steady-state plant outputs and inputs respectively. For higher LP execution frequencies than transitions between steady-states, the MPC disturbance estimate may be used [19], ^ðk j k 1Þ: d¼yy
2.1. MPC The model predictive control implementation we use here is quadratic dynamic matrix control (QDMC) [3], using a state-space formulation of the internal dynamic model [9]. The model predictions take the form, ^ ðk þ i 1 j kÞ; i ¼ 1; . . . ; M ð1aÞ ^ðk þ i j kÞ ¼ Ax ^ðk þ i 1 j kÞ þ Bu x ^ ðk þ M 1 j kÞ; i ¼ M þ 1; . . . ; P ^ðk þ i j kÞ ¼ Ax ^ðk þ i 1 j kÞ þ Bu x ð1bÞ ^ ^ ^ yðk þ i j kÞ ¼ C xðk þ i j kÞ þ dðk þ i j kÞ; i ¼ 1; . . . P; ð1cÞ
3. Sensitivity analysis: MISO systems For simplicity of analysis, a multi-input, single-output (MISO) system is considered first; such a system has only one bias term in its formulation. The LP then has the form: þ min ay ;y u
s:t:
¼ y
m X i¼1 m X
i bi u
g ss i ui þ d
ð3Þ
i¼1
where P is the prediction horizon, M is the control move horizon, ^ ðk þ i j kÞ 2 Rp represents the predicted values of the outputs at y time step k þ i based on information available at time step k, and ^ ðk þ i j kÞ 2 Rm , ^ ðk þ i j kÞ 2 Rn , and inputs, u the predicted states, x are similarly defined. A 2 Rnn ; B 2 Rnm and C 2 Rpn are model coefficient matrices. The disturbance estimate in the originally proposed Dynamic Matrix Control (DMC) and QDMC algorithms [1,3] is computed as the difference between the measured and predicted outputs, and assumed constant over the prediction horizon. In the present framework, this becomes:
6 ymax ymin 6 y i 6 umax umin 6u ; i i
i ¼ 1; . . . ; m
where a and bi are price coefficients and g ss i are steady-state gains of the process. According to the properties of an LP, a finite solution to problem (3) lies at the boundary of the feasible region, which is determined by the inequality constraints. For fixed cost function coefficients and optimization variable bounds, the particular location of the solution depends only on the value of the bias d. From this viewpoint, the LP problem can be considered as a mapping from R1 ! Rmþ1 , where m þ 1 is the number of optimization variables.
A. Nikandrov, C.L.E. Swartz / Journal of Process Control 19 (2009) 16–24
We first consider the range of d for which the LP is feasible. This can be done by reformulating the optimization problem (3), posing d as an optimization variable. The upper bound for d can be found as 8 max d 8 > m > P > > max y > subject to : g ss > > > > i ui > > u ; y > i > i¼1 < < m P g ss d¼y subject to : ð4Þ ) i ui > > i¼1 > > > > 6 ymax > ymin 6 y > min max > > 6y : min > y 6y > > i 6 umax ui 6 u : min i i 6 umax ui 6 u i
1.5 y u1 u2
1
0.5
y, u1, u2
18
—0.5
The lower bound can be found analogously to (4), formulated as a minimization problem. If all the steady-state gains g ss i are either nonnegative or nonpositive, the solution of the optimization problems can be written analytically as
dmax
g ss i P 0; 8i m P min ¼ ymax g ss i ui
dmax
i¼1
dmin ¼ ymin
m P i¼1
max g ss i ui
g ss i 6 0; 8i m P max ¼ ymax g ss i ui i¼1
dmin ¼ ymin
m P i¼1
1.5
function can be found using simple calculations. For d < 1:0, ¼ ymin ¼ 1:0 and u 2 ¼ umax y ¼ 0:5. Then, the relation between 2 1 and d can be expressed as follows: u ss ¼ g ss y 1 u1 þ g 2 u2 þ d ) u1 ¼
a¼
u 1 maxy;u1 ;u2 f ¼ y ¼ 0:4u 1 þ 0:4u 2 þ d s:t: y ð6Þ
2 6 0:5 0:5 6 u First, we calculate the range of the bias term d for which a feasible solution of the LP formulation (6) can be calculated. Since the steady-state gains are both positive, formulas (5) can be used:
dmin ¼ ymin
1
1 1 1 ss g ss 2 ss d y u g ss g1 2 g1 1 ð7Þ
where
The objective of maximizing f ¼ y u1 gives an LP of the following form:
i¼1 m X
0.5
) u1 jd<1:0 ¼ ad þ b;
0:5 6 u1 ; u2 6 0:5:
m X
0
Fig. 2. LP solution as a function of bias for case study 1.
min g ss i ui
6 1:0 1:0 6 y 1 6 0:5 0:5 6 u
—0.5
ð5Þ
with constraints:
dmax ¼ ymax
—1
d
0:4 0:4 u1 ðsÞ þ u2 ðsÞ 3s þ 1 5s þ 1
1:0 6 y 6 1:0;
—1
—1.5 —1.5
In all other cases the solution depends on the sign of each steadystate gain. The effect of the bias term on the LP solution is illustrated through the following case study. Case study 1. Consider a MISO process with two inputs described by following model: yðsÞ ¼
0
min g ss ¼ 1:0 ð0:4Þð0:5Þ ð0:4Þð0:5Þ ¼ 1:4 i ui
max g ss ¼ 0:1 ð0:4Þð0:5Þ ð0:4Þð0:5Þ ¼ 1:4 i ui
i¼1
If the bias term is within this range, the solution of the LP exists. The dependance of this solution on the bias is shown graphically in Fig. 2. Each line in the figure represents an optimization variable. For any feasible value of d, at least two variables lie at their constraints. There are two particular values of the bias (d ¼ 1:0 and d ¼ 1:0) where all three variables lie at their constraints. Moreover, at these points the constraints change their activity. When is active while 0:5 6 u 1 is not active, d < 1:0, constraint 1 6 y is 1 is active while 1 6 y and when d > 1:0, constraint 0:5 6 u not active. This means that for any given bias value, some constraints are active while others may not be, and there are particular values of the bias around which small deviations cause changes in the constraint activity. When an optimization variable in the above problem is not at a constraint, its value is a linear function of the bias. The slope of this
1 ; g ss 1
b¼
1 min 1 max y ss g ss 2 u2 g ss g 1 1
From (7) it follows that the slope of this linear function is the negative of the reciprocal of the corresponding steady-state gain. Since ss the steady-state gains for u1 and u2 are the same ðg ss 1 ¼ g 2 ¼ 0:4Þ, their solution lines are collinear. Therefore, if a steady-state gain is small then the linear function which represents the solution for the corresponding input has a large slope and vice versa. It is important to take this observation into account because deviations in the bias value could result in larger or smaller variations in the LP solution, depending on the values of the steady-state gains and steadystate bias value. is not at a constraint, it varies From Fig. 2 we observe that if y linearly with d with slope of unity. This is true in general for MISO systems for which ag ss i þ bi –0, and can be seen as follows. , By using the equality constraint in problem (3) to eliminate y we may formulate the LP equivalently as ss ss minðag ss 1 þ b1 Þu1 þ ðag 2 þ b2 Þu2 þ þ ðag m þ bm Þum þ ad u
s:t:
ss max d ymin d 6 g ss 1 u1 þ þ g m um 6 y
umin i
i 6 6u
umax ; i
ð8Þ
i ¼ 1; . . . ; m
If neither of the first inequality constraints is active, then by inspection the optimal solution is given by max ui ; ðag ss i þ bi Þ < 0 i ¼ u ss ; ðag umin i i þ bi Þ > 0 and the optimal output, ¼ y
m X
g ss i ui þ d
i¼1
is linear in d with unity gain. We now consider the response of the LP solution to noise in the bias term, which could be a result of high-frequency disturbances
A. Nikandrov, C.L.E. Swartz / Journal of Process Control 19 (2009) 16–24
19
1 constraints, the optimization problem can be cast in terms of u 2 alone using direct substitution. Optimization problem (9) and u can thus be reformulated as: ss ss ss min f ¼ a1 g ss 11 þ a2 g 21 þ b1 u1 þ a1 g 12 þ a2 g 22 þ b2 u2 1 ;u 2 u
þ a1 d1 þ a2 d2 s:t: ss max ymin 6 g ss 1 11 u1 þ g 12 u2 þ d1 6 y1 ss max ymin 6 g ss 2 21 u1 þ g 22 u2 þ d2 6 y2
1 6 umax umin 6u 1 1 2 6 umax umin 6u 2 2 This could also be written as ss ss ss min f ¼ a1 g ss 11 þ a2 g 21 þ b1 u1 þ a1 g 12 þ a2 g 22 þ b2 u2 1 ;u 2 u
s:t:
8 ss max > 1 þ y1 gssd1 2 6 gg11 1 :u 2 P gg11 þ u ss 1 g ss 12 12 8 ss max 6 gg21 1 þ y2 gssd2 > u ss u > < 2 22 22 g ss ymin d2 pair2 21 2 P þ u u > 2 1 g ss g ss > 22 22 :
Fig. 3. Effect of bias noise on LP solution for case study 1.
ð10Þ
1 6 umax 6u umin 1 1 or measurement noise. The effect of noise applied to the LP bias with a steady-state value 1.0 on the LP solution is shown in Fig. 3. Since the bias steady-state value is 1.0 and the noise is mod 1 is constant, as expected from Fig. 2. From erate, the solution for u Fig. 2 we deduce that if the noise is negative (steady-state bias plus 2 is at its upper bound and the solution for noise is less than 1.0), u the output replicates the noise. If the noise is positive (noisy bias is greater than 1.0), the output is at its maximum value, while the 2 is amplified by 2:5 in the opposite direction to that solution for u of the noise. These results are verified by the simulation results in Fig. 3. We observe in this section that for MISO systems, the output set-point calculated by the LP is either constant or varies linearly with the bias parameter with gain of unity, assuming unique LP solutions for specified bias parameter values. However, the variation in the input target values may be larger than that of the bias. We show in the next section that for MIMO systems, the variation in output set-points may exceed that of the bias term due to interaction effects. 4. Sensitivity analysis: MIMO systems We focus in this section two-input, two-output systems. Key characteristics extend readily to the more general MIMO case. The LP optimization for the general two-input, two-output system takes the form: min
1 ;u 2 ;y 1 ; y 2 u
s:t:
1 þ a2 y 2 þ b1 u 1 þ b2 u 2 a1 y
ss 1 ¼ g ss y 11 u1 þ g 12 u2 þ d1 ss ss y2 ¼ g 21 u1 þ g 22 u2 þ d2
1 6 ymax ymin 6y 1 1
ð9Þ
2 6 ymax ymin 6y 2 2 1 6 umax umin 6u 1 1 2 6 umax umin 6u 2 2 A key difference between this formulation and the formulations considered previously is that two bias terms are now involved, and they are both updated at every iteration. Since optimization 2 relate to u 1 and u 2 through the model equality 1 and y variables y
2 6 umax umin 6u 2 2 The term a1 d1 þ a2 d2 in the objective function is omitted since it is constant for any given biases d1 and d2 , and does not affect the solution. From formulation (10), several inferences can be made regarding the feasible region of the problem and possible location of its 2 plane, the slope of the line which rep1 u solution. First, on the u resents the objective function (so-called ‘‘isocost line” [18]) does not depend on the bias terms. Its slope depends only on the cost coefficients and the steady-state gains of the process. Second, the output inequality constraints are presented in pairs, where each pair comprises two collinear constraints. The slope of each pair of constraints depends on the ratio of the steady-state gains, ss ss ss g ss i1 =g i2 ; i ¼ 1; 2. If the steady-state gains, g i1 and g i2 , are the same sign then the slope of the constraint is negative, and vice versa. The distance between the constraints in each pair does not depend on the bias term, and is only determined by the steady-state gains and the original output constraint limits. The distances between the constraints in each pair are D1 ¼
ymax d1 ymin d1 ymax ymin 1 1 ss ¼ 1 ss 1 g ss g g 12 12 12
D2 ¼
ymax d2 ymin d2 ymax ymin 2 2 ss ¼ 2 ss 2 : g ss g g 22 22 22
If d1 is an arithmetic mean of ymax and ymin (or d2 for the pair ymax 1 1 2 and ymin 2 ) then this pair of constraints is symmetric around the ori) or down (if di is gin. Otherwise it is shifted up (if di is closer to ymin i 2 versus u 1 plot (see Fig. 4). The feasible area ), on a u closer to ymax i for each pair of constraints lies between the constraint lines. Third, the minimum and maximum constraints for the input variables form a rectangle with the feasible area inside. The resulting feasible region of the problem is the intersection of all the constraint feasible regions and is schematically presented in Fig. 4. The solution of the optimization problem can be found by shifting the isocost line within the feasible region in the direction of decreasing cost until it cannot be improved any further. The resulting position of the isocost line will determine the optimal values of the optimization variables. At every iteration, the biases d1 and d2 are updated. If these values differ from their previous values, the shape of the feasible re-
20
A. Nikandrov, C.L.E. Swartz / Journal of Process Control 19 (2009) 16–24
Fig. 4. Graphical representation of the LP optimization problem for a 2 2 system.
gion changes, which may impact the solution of the optimization problem. As mentioned above, none of the slopes depends on the bias terms and the feasible region changes only because the constraint lines shift up and down. We consider here three scenarios in which the effect of bias variations on the optimal LP solution, hence MPC set-points, is explored. 4.1. Scenario 1 In this scenario, illustrated in Fig. 5, the optimal solution is deand ymax constraints. If the fined by the intersection of the ymax 1 2 combination of the steady-state gains is such that constraints are almost collinear, then the LP solution is expected to be sensitive to small changes in the biases. If one of the biases changes slightly, the new intersection of the shifted constraint lines (and, therefore, new solution) may appear some distance from the previous solution as shown in Figs. 5(a) and (b). However, the output set-points determined by the solution of the LP remain constant at the upper output constraint limits.
This situation is illustrated in Fig. 6. Analogously to the previous case, the solution of the LP can become sensitive to the changes in the bias term if the constraint lines are almost collinear. However, unlike the previous scenario, only one output set-point will be constant, and the other may fluctuate in response to variations in the bias value. In fact, we show in the case study below that it is possible for bias variations to be amplified in the set-point of the second output. Case study 2 We consider here a plant described by the following transfer function model, 1:2 1:4 u1 ðsÞ u2 ðsÞ 3s þ 1 3s þ 1 0:5 5:5 y2 ðsÞ ¼ u1 ðsÞ þ u2 ðsÞ 3s þ 1 3s þ 1 y1 ðsÞ ¼
ð11Þ
with constraints 1:0 6 y1 6 1:0 1:0 6 y2 6 1:0 0:5 6 u1 6 0:5
4.2. Scenario 2
0:5 6 u2 6 0:5
Here, the optimal solution is defined by the intersection of one 2 . of the output constraint limits and the maximum constraint on u
Maximization of f ¼ 0:3y1 þ 0:5y2 þ 0:6u1 þ 8:3u2 yields the corresponding LP:
21
A. Nikandrov, C.L.E. Swartz / Journal of Process Control 19 (2009) 16–24
Fig. 5. Effect of the bias on the LP solution: Scenario 1.
Fig. 6. Effect of the bias on the LP solution: Scenario 2.
1 þ 0:5y 2 þ 0:6u 1 þ 8:3u 2 max f ¼ 0:3y
1 ;y 2 ;u 1 ;u 2 y
s:t:
1 ¼ 1:2u 1 1:4u 2 þ d1 y 1 þ 5:5u 2 þ d2 2 ¼ 0:5u y 1 6 1:0 1:0 6 y 2 6 1:0 1:0 6 y 1 6 0:5 0:5 6 u 2 6 0:5 0:5 6 u
ð12Þ
Fluctuations around the steady-state bias values, d1 ¼ 0:4, d2 ¼ 1:8, results in the scenario depicted in Fig. 6. Simulations were conducted in which bias d2 was subjected to white noise, while bias d1 was kept at its steady-state value. The results of the simulation for manipulated variable (MV) targets and controlled variable (CV) set-points are presented in Fig. 7. Since the solution of the LP lies at the intersection of the constraints for y2 and u2 , the corresponding set-points and target
values are constant throughout the simulation period. Fig. 7(a) 1 with larger variashows that noise in d2 causes fluctuations in u 2 ¼ umax 2 is always equal to ymax ¼ 1 and u , tion. Indeed, since y 2 2 then the LP solution for any value of bias d2 must satisfy the equality: ss max 1 ¼ ymax ¼ g ss þ d2 ) u 2 21 u1 þ g 22 u2
1 d2 þ const: g ss 21
g ss 21 ¼ 0:5 and therefore, any changes in d2 cause changes in u1 that are twice as large and in the opposite direction. At the same time, 1 through the steady-state model equality 1 is related to u since y 1 2 and d1 are constants, any changes in u constraints in (12), and u 1 with coefficient g ss result in changes in y (which is 1:2 in this 11 case study). Summarizing both effects, it can be concluded that ss 1 such that: y 1 ¼ ðg ss noise in d2 results in changes in y 11 =g 21 Þd2 . In ss this case study g ss 11 =g 21 ¼ 2:4, and the effect on y1 is confirmed
22
A. Nikandrov, C.L.E. Swartz / Journal of Process Control 19 (2009) 16–24
Fig. 7. 2 2 system sensitivity: Case study 2. (a) Biases and input target values. (b) Biases and output set points.
1 are ampliby the results presented in Fig. 7. Thus, variations in y 2 remains constant. fied, while y 4.3. Scenario 3 In this scenario, we consider the solution to lie at the intersection of an output and input constraint as shown in Fig. 8. Variation in the bias can again cause amplification in one of the output setpoints, depending on the values of the process gains. Consider the response of the optimal solution to variations in 1 and y 1 at their maximum constraints. the bias term, d1 , with u From (10), we would have 2 ¼ u
g ss ymax d1 1 11 max u1 þ 1 ss ¼ ss d1 þ const g ss g 12 g 12 12
2 depends on the value Therefore, the sensitivity of the solution for u 2 also result in non. The fluctuations in u of the steady-state gain g ss 12 constant set-points for y2 . If g ss 22 is larger than unity, then any per 2 will be amplified, inducing even larger variations turbations in u 2 , and vice versa. in y
5. Performance of two-level LP-MPC system We consider here the performance of a two-level LP-MPC control system with noise applied to the bias. Graphically, the system is shown in Fig. 9. Case study 3 This case study is an extension of case study 2. It comprises the same plant, constraints and LP objective function, but we consider the closed-loop LP-MPC cascade system. The controlled variables are y1 and y2 , and the MPC parameters are as follows: Q ¼ I, R ¼ 0, S ¼ I, P ¼ 50, M ¼ 2, with sampling period, T s ¼ 0:3. The steady-state bias values were achieved by introducing output step-like disturbances passing through first-order filters with unity steady-state gains. Also, white noise was applied to output y2 , while output y1 was maintained uncorrupted. A perfect plant model was assumed. The results of the two-level control system operation after the transient effects of the step disturbance had dissipated are presented in Fig. 10. White noise with variance 0.0013 caused fluctu1 with variance 0.0074. 1 with variance 0.0051 and in y ations in u
Fig. 8. Effect of the bias on the LP solution: Scenario 3.
A. Nikandrov, C.L.E. Swartz / Journal of Process Control 19 (2009) 16–24
23
6. Conclusions
Fig. 9. Two-level LP-MPC control system.
This confirms that in this case study, white noise in the second output channel may cause fluctuations in the first output set-point with much larger variation. Even though the set-points for first output are corrupted by the noise, the control performance for the corresponding output is quite good. This could likely be ascribed to the relatively short control move horizon, as control move severity tends to attenuate with decreasing M. The variance of the output (0.000493) is much smaller than that of its set-point. However the performance of output y2 is quite poor (its variance is 0.0018) considering that its set point is constant at steady-state. This can be explained by fact that fast changes in yset 1 cause frequent changes in u1 which in turn affects y2 . u2 operates near its upper constraint and becomes saturated for periods of time during which it cannot be used for control.
Industrial implementations of MPC typically include a steadystate optimization layer to track the economic optimum and to provide feasible set-points to the model predictive controller. While two-level LP-MPC systems are relatively widely applied, reports of excessive set-point variation motivate the investigation of potential causes of such behavior. The present study focused on the sensitivity of the LP solution to variations in the LP model bias term, since the bias is updated prior to each LP execution based on the difference between measured and predicted plant outputs. A sensitivity analysis was conducted first for multi-input, single output systems. The steady-state LP model has only one bias term in this case, and its effect on the LP solution is conveniently summarized graphically. A key observation is that if the output is not at one of its limits, it varies linearly with the bias term with gain of unity, provided that the LP solution is unique. However, changes in the inputs may exceed those of the bias term if the gain relating the output to the corresponding input is less than unity. MIMO systems were explored through a detailed analysis of a two-input, two-output system. A graphical analysis was developed by expressing the problem in terms of the inputs only. Scenarios in which the solution is defined by two output constraints, and an output and an input constraint, were considered. In the former case, the set-points remain constant for small perturbations in the bias term. However, in the latter case, variations in the bias may result in amplified variations in the output that is not at its bound, depending on the magnitudes of the process gains. Simulation of the two-level LP-MPC system was conducted to evaluate the response of the two-level closed-loop system to a step disturbance combined with measurement noise. The results of the sensitivity analysis were used to explain the observed performance of the two-level system.
Fig. 10. Two-level control system response: Case study 3.
24
A. Nikandrov, C.L.E. Swartz / Journal of Process Control 19 (2009) 16–24
The focus of this paper has been on analysis of the sensitivity of LP-generated set-points to variations in the bias term induced, for example, by measurement and/or process noise. A question that naturally arises is how adverse amplification effects can be avoided or mitigated. One possibility would be the use of an alternative control structure. Another option might be to filter the bias update; however, the impact on the overall performance of the two-level LP-MPC system would need to be evaluated. A systematic investigation of these and other factors within a design framework would be a potentially useful avenue for future research. References [1] C.R. Cutler, B.L. Ramaker, Dynamic matrix control – a computer control algorithm, in: AIChE 86th National Meeting, Houston, TX, 1979. [2] J.F. Forbes, T.E. Marlin, Model accuracy for economic optimizing controllers: the bias update case, Ind. Eng. Chem. Res. 33 (1994) 1919–1929. [3] C.E. Garcia, A.M. Morshedi, Quadratic programming solution of dynamic matrix control (QDMC), Chem. Eng. Commun. 46 (1986) 73–87. [4] C.E. Garcia, D.M. Prett, M. Morari, Model predictive control: theory and practice – a survey, Automatica 25 (3) (1989) 335–348. [5] D.E. Kassmann, T.A. Badgwell, R.B. Hawkins, Robust steady-state target calculation for model predictive control, AIChE J. 46 (5) (2000) 1007–1024. [6] D.J. Kozub, Controller performance monitoring and diagnosis. Industrial perspective, in: Fifteenth Triennial World Congress of the International Federation of Automatic Control, Barcelona, Spain, 2002.
[7] J.H. Lee, M. Morari, C.E. Garcia, State-space interpretation of model predictive control, Automatica 30 (4) (1994) 707–717. [8] S. Li, K.Y. Lim, D.G. Fisher, A state space formulation for model predictive control, AIChE J. 35 (2) (1989) 241–249. [9] J.M. Maciejowski, Predictive Control with Constraints, Springer, 2002. [10] M. Morari, J.H. Lee, Model predictive control: past, present and future, Comput. Chem. Eng. 23 (1999) 667–682. [11] K.R. Muske, Steady-state target optimization in linear model predictive control, in: Proceedings of ACC, Albuquerque, NM, 1997. [12] D.M. Prett, M. Morari, The Shell Process Control Workshop, Butterworths, Stoneham, MA, 1987. [13] S.J. Qin, T.A. Badgwell, A survey of industrial model predictive control technology, Control Eng. Practice 11 (2003) 733–764. [14] C.V. Rao, J.B. Rawlings, Steady states and constraints in model predictive control, AIChE J. 45 (6) (1999) 1266–1278. [15] N.L. Ricker, Model predictive control with state estimation, Ind. Eng. Chem. Res. 29 (1990) 374–382. [16] S.L. Shah, R. Patwardhan, B. Huang, Multivariate controller performance analysis: Methods, applications and challenges, in: Chemical Process Control VI, Tuscon, Arizona, 2001. [17] R.C. Sorensen, C.R. Cutler, LP integrates economics into dynamic matrix control, Hydocarbon Process. (September) (1998) 57–65. [18] W.L. Winston, Operations Research: Applications and Algorithms, fourth ed., Duxbury Press, 2004. [19] C.-M. Ying, B. Joseph, Performance and stability analysis of LP-MPC and QPMPC cascade control systems, AIChE J. 45 (7) (1999) 1521–1534. [20] C. Yousfi, R. Tournier, Steady state optimization inside model predictive control, in: Proceedings of ACC, Boston, 1991.
Contents lists available at ScienceDirect
Journal of Process Control journal homepage: www.elsevier.com/locate/jprocont
Sensitivity analysis of LP-MPC cascade control systems Alexei Nikandrov, Christopher L.E. Swartz * Department of Chemical Engineering, McMaster University, 1280 Main Street West, Hamilton, Ontario, Canada L8S 4L7
a r t i c l e
i n f o
Article history: Received 26 June 2007 Received in revised form 13 February 2008 Accepted 15 February 2008
Keywords: Linear programming Model predictive control LP-MPC Sensitivity
a b s t r a c t Model predictive control (MPC) has found wide application in the chemical process industry as well as other industrial sectors. Commercial MPC systems are typically implemented in conjunction with a steady-state linear or quadratic programming optimizer, whose key functions are to track the economic optimum and to provide feasible set-points to the model predictive controller. The two-level system is complementary to real-time optimization which typically utilizes more complex models and is executed less frequently. Despite the widespread adoption of LP-MPC systems, occurrences of poor performance have been reported, where large variations in the computed set-points were observed. In this paper, we analyze the sensitivity of the LP solution to variation in the LP model bias, through which feedback to the LP layer occurs. We consider both multi-input, single-output (MISO) and multi-input, multi-output (MIMO) systems. Principles are illustrated through graphical representation as well as case studies. The performance of the two-level LP-MPC closed-loop system is evaluated and explained using results of the LP sensitivity analysis. Ó 2008 Elsevier Ltd. All rights reserved.
1. Introduction Model predictive control (MPC) is arguably the advanced control algorithm of choice in the chemical process industry, and has made inroads into other industrial sectors as well [10,13]. A dynamic model is utilized within the control algorithm to predict the effect of future plant inputs on the controlled outputs. Future inputs are computed in accordance with a performance objective, typically as the solution of an optimization problem. The inputs corresponding to the first control interval are implemented, and the calculation process repeated at the end of the sampling period, with the model predictions adjusted using the difference between the measured and predicted outputs. Details of the algorithm may be found, inter alia, in [4,9,13]. Industrial MPC systems are generally implemented in conjunction with a linear programming (LP) or quadratic programming (QP) steady-state optimizer [5,13,17,19,20]. The LP (or QP) typically uses a static model consistent with the dynamic MPC model, and is executed at the same frequency as the model predictive controller. The plant economic optimum may shift due to disturbances; thus the steady-state LP (QP) provides a bridge between a higher-level and less frequently executed real-time optimization (RTO) layer and the model predictive controller by making setpoint adjustments in response to changing conditions between RTO executions. The LP formulation may involve minimization of the deviation between the set-points and target values determined * Corresponding author. Tel.: +1 905 525 9140; fax: +1 905 521 1350. E-mail address: [email protected] (C.L.E. Swartz). 0959-1524/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.jprocont.2008.02.005
by the real-time optimizer, or optimization of an economic criterion directly. Fig. 1 illustrates the location of the steady-state LP (QP) within a plant automation hierarchy, as given in [19]. We note that the process would typically include local PID-type controllers and that several variants are possible, such as the presence of a plant unit optimization layer between the LP and plant-wide RTO layers [5,13]. Despite the apparent success of two-level LP-MPC systems, instances of poor performance have been reported [6,16]. Shah et al.[16], in the context of control performance monitoring, describe an industrial application in which the variation in the set-points exceeds that of the corresponding controlled variables. Kozub[6] also reports set-points being noisy relative to their controlled variables, and excessive variation in the set of inputs which are at their constraints. This motivates an investigation into the potential causes of such behavior. Ying and Joseph [19] provide stability theorems for LP-MPC and QP-MPC cascade control systems with no plant/model mismatch. Consideration of model uncertainty is included in a case study based on the Shell Standard Control Problem [12]. Kassmann et al.[5] present a formulation for robust steady-state target calculation. For elliptic uncertainty on the model parameters, the target calculation problem takes the form of a second order cone program (SOCP) which the authors solve using software based on a primaldual interior point algorithm. Steady-state targets are also computed in [11,14]. However, this is driven primarily not by economics, but rather to provide steady-state values for the plant inputs and states for inclusion in an MPC formulation in which
A. Nikandrov, C.L.E. Swartz / Journal of Process Control 19 (2009) 16–24
17
^ j kÞ ¼ yðkÞ C x ^ ðk j k 1Þ dðk ^ ^ dðk þ i j kÞ ¼ dðk j kÞ; i ¼ 1; . . . ; P where ^ ðk j kÞ ¼ x ^ðk j k 1Þ ¼ Ax ^ðk 1 j k 1Þ þ Buðk 1Þ: x The above estimation scheme is based on an assumption of steplike disturbances, and may result in poor control performance in the presence of measurement noise and certain other disturbance types. A more general state estimation framework that addresses these limitations is described in [7–9,15]. In the application studies that follow, we consider only the original DMC disturbance estimation scheme. 2.2. LP The set-points are determined through the solution of an LP of the form Fig. 1. LP (QP) within process automation hierarchy [19].
deviations of the states and inputs from corresponding steadystates are weighted by positive definite matrices. Feedback in the two-level LP-MPC configuration occurs through a bias term in the steady-state LP model. In this paper, we analyze the sensitivity of the LP solution to variations in the LP model bias. We consider first the LP level separately for multi-input, singleoutput (MISO) and multi-input, multi-output (MIMO) systems. This analysis is coupled with graphical representations to provide insights into the sensitivity effects. Thereafter, the performance of a two-level LP-MPC cascade system is evaluated, with observed performance related to earlier sensitivity results. 2. Problem definition The control structure we consider in this work is a two-layer LPMPC system in which the LP provides output set-points and/or input target values to the model predictive controller. The MPC-controlled system feeds back information to the LP layer, from which model adjustments may be determined. We present in this section the MPC and LP formulations used, and also describe the LP model bias update method.
minu ;y s:t:
p X
i þ ai y
m X
i bi u
i¼1
i¼1 ss
þd y¼G u y
min
6y6y
ð2Þ
max
umin 6 u 6 umax where Gss is the plant model gain matrix and d is a model bias that is updated in accordance with the actual plant outputs. We use the overbar to identify variables involved in the LP optimization. Com and y become the MPC set-points (includponents of the optimal u ing input target values). We typically assign a number of set-points equal to the number of manipulated inputs; specification of an excess number of set-points would in general result in offset. In real-time optimization, the bias update, d, is typically calculated as the difference between the measured outputs and those resulting from applying the plant inputs to the steady-state model [2]. In the present framework, this would correspond to d ¼ y Gss u; where y and u correspond to the steady-state plant outputs and inputs respectively. For higher LP execution frequencies than transitions between steady-states, the MPC disturbance estimate may be used [19], ^ðk j k 1Þ: d¼yy
2.1. MPC The model predictive control implementation we use here is quadratic dynamic matrix control (QDMC) [3], using a state-space formulation of the internal dynamic model [9]. The model predictions take the form, ^ ðk þ i 1 j kÞ; i ¼ 1; . . . ; M ð1aÞ ^ðk þ i j kÞ ¼ Ax ^ðk þ i 1 j kÞ þ Bu x ^ ðk þ M 1 j kÞ; i ¼ M þ 1; . . . ; P ^ðk þ i j kÞ ¼ Ax ^ðk þ i 1 j kÞ þ Bu x ð1bÞ ^ ^ ^ yðk þ i j kÞ ¼ C xðk þ i j kÞ þ dðk þ i j kÞ; i ¼ 1; . . . P; ð1cÞ
3. Sensitivity analysis: MISO systems For simplicity of analysis, a multi-input, single-output (MISO) system is considered first; such a system has only one bias term in its formulation. The LP then has the form: þ min ay ;y u
s:t:
¼ y
m X i¼1 m X
i bi u
g ss i ui þ d
ð3Þ
i¼1
where P is the prediction horizon, M is the control move horizon, ^ ðk þ i j kÞ 2 Rp represents the predicted values of the outputs at y time step k þ i based on information available at time step k, and ^ ðk þ i j kÞ 2 Rm , ^ ðk þ i j kÞ 2 Rn , and inputs, u the predicted states, x are similarly defined. A 2 Rnn ; B 2 Rnm and C 2 Rpn are model coefficient matrices. The disturbance estimate in the originally proposed Dynamic Matrix Control (DMC) and QDMC algorithms [1,3] is computed as the difference between the measured and predicted outputs, and assumed constant over the prediction horizon. In the present framework, this becomes:
6 ymax ymin 6 y i 6 umax umin 6u ; i i
i ¼ 1; . . . ; m
where a and bi are price coefficients and g ss i are steady-state gains of the process. According to the properties of an LP, a finite solution to problem (3) lies at the boundary of the feasible region, which is determined by the inequality constraints. For fixed cost function coefficients and optimization variable bounds, the particular location of the solution depends only on the value of the bias d. From this viewpoint, the LP problem can be considered as a mapping from R1 ! Rmþ1 , where m þ 1 is the number of optimization variables.
A. Nikandrov, C.L.E. Swartz / Journal of Process Control 19 (2009) 16–24
We first consider the range of d for which the LP is feasible. This can be done by reformulating the optimization problem (3), posing d as an optimization variable. The upper bound for d can be found as 8 max d 8 > m > P > > max y > subject to : g ss > > > > i ui > > u ; y > i > i¼1 < < m P g ss d¼y subject to : ð4Þ ) i ui > > i¼1 > > > > 6 ymax > ymin 6 y > min max > > 6y : min > y 6y > > i 6 umax ui 6 u : min i i 6 umax ui 6 u i
1.5 y u1 u2
1
0.5
y, u1, u2
18
—0.5
The lower bound can be found analogously to (4), formulated as a minimization problem. If all the steady-state gains g ss i are either nonnegative or nonpositive, the solution of the optimization problems can be written analytically as
dmax
g ss i P 0; 8i m P min ¼ ymax g ss i ui
dmax
i¼1
dmin ¼ ymin
m P i¼1
max g ss i ui
g ss i 6 0; 8i m P max ¼ ymax g ss i ui i¼1
dmin ¼ ymin
m P i¼1
1.5
function can be found using simple calculations. For d < 1:0, ¼ ymin ¼ 1:0 and u 2 ¼ umax y ¼ 0:5. Then, the relation between 2 1 and d can be expressed as follows: u ss ¼ g ss y 1 u1 þ g 2 u2 þ d ) u1 ¼
a¼
u 1 maxy;u1 ;u2 f ¼ y ¼ 0:4u 1 þ 0:4u 2 þ d s:t: y ð6Þ
2 6 0:5 0:5 6 u First, we calculate the range of the bias term d for which a feasible solution of the LP formulation (6) can be calculated. Since the steady-state gains are both positive, formulas (5) can be used:
dmin ¼ ymin
1
1 1 1 ss g ss 2 ss d y u g ss g1 2 g1 1 ð7Þ
where
The objective of maximizing f ¼ y u1 gives an LP of the following form:
i¼1 m X
0.5
) u1 jd<1:0 ¼ ad þ b;
0:5 6 u1 ; u2 6 0:5:
m X
0
Fig. 2. LP solution as a function of bias for case study 1.
min g ss i ui
6 1:0 1:0 6 y 1 6 0:5 0:5 6 u
—0.5
ð5Þ
with constraints:
dmax ¼ ymax
—1
d
0:4 0:4 u1 ðsÞ þ u2 ðsÞ 3s þ 1 5s þ 1
1:0 6 y 6 1:0;
—1
—1.5 —1.5
In all other cases the solution depends on the sign of each steadystate gain. The effect of the bias term on the LP solution is illustrated through the following case study. Case study 1. Consider a MISO process with two inputs described by following model: yðsÞ ¼
0
min g ss ¼ 1:0 ð0:4Þð0:5Þ ð0:4Þð0:5Þ ¼ 1:4 i ui
max g ss ¼ 0:1 ð0:4Þð0:5Þ ð0:4Þð0:5Þ ¼ 1:4 i ui
i¼1
If the bias term is within this range, the solution of the LP exists. The dependance of this solution on the bias is shown graphically in Fig. 2. Each line in the figure represents an optimization variable. For any feasible value of d, at least two variables lie at their constraints. There are two particular values of the bias (d ¼ 1:0 and d ¼ 1:0) where all three variables lie at their constraints. Moreover, at these points the constraints change their activity. When is active while 0:5 6 u 1 is not active, d < 1:0, constraint 1 6 y is 1 is active while 1 6 y and when d > 1:0, constraint 0:5 6 u not active. This means that for any given bias value, some constraints are active while others may not be, and there are particular values of the bias around which small deviations cause changes in the constraint activity. When an optimization variable in the above problem is not at a constraint, its value is a linear function of the bias. The slope of this
1 ; g ss 1
b¼
1 min 1 max y ss g ss 2 u2 g ss g 1 1
From (7) it follows that the slope of this linear function is the negative of the reciprocal of the corresponding steady-state gain. Since ss the steady-state gains for u1 and u2 are the same ðg ss 1 ¼ g 2 ¼ 0:4Þ, their solution lines are collinear. Therefore, if a steady-state gain is small then the linear function which represents the solution for the corresponding input has a large slope and vice versa. It is important to take this observation into account because deviations in the bias value could result in larger or smaller variations in the LP solution, depending on the values of the steady-state gains and steadystate bias value. is not at a constraint, it varies From Fig. 2 we observe that if y linearly with d with slope of unity. This is true in general for MISO systems for which ag ss i þ bi –0, and can be seen as follows. , By using the equality constraint in problem (3) to eliminate y we may formulate the LP equivalently as ss ss minðag ss 1 þ b1 Þu1 þ ðag 2 þ b2 Þu2 þ þ ðag m þ bm Þum þ ad u
s:t:
ss max d ymin d 6 g ss 1 u1 þ þ g m um 6 y
umin i
i 6 6u
umax ; i
ð8Þ
i ¼ 1; . . . ; m
If neither of the first inequality constraints is active, then by inspection the optimal solution is given by max ui ; ðag ss i þ bi Þ < 0 i ¼ u ss ; ðag umin i i þ bi Þ > 0 and the optimal output, ¼ y
m X
g ss i ui þ d
i¼1
is linear in d with unity gain. We now consider the response of the LP solution to noise in the bias term, which could be a result of high-frequency disturbances
A. Nikandrov, C.L.E. Swartz / Journal of Process Control 19 (2009) 16–24
19
1 constraints, the optimization problem can be cast in terms of u 2 alone using direct substitution. Optimization problem (9) and u can thus be reformulated as: ss ss ss min f ¼ a1 g ss 11 þ a2 g 21 þ b1 u1 þ a1 g 12 þ a2 g 22 þ b2 u2 1 ;u 2 u
þ a1 d1 þ a2 d2 s:t: ss max ymin 6 g ss 1 11 u1 þ g 12 u2 þ d1 6 y1 ss max ymin 6 g ss 2 21 u1 þ g 22 u2 þ d2 6 y2
1 6 umax umin 6u 1 1 2 6 umax umin 6u 2 2 This could also be written as ss ss ss min f ¼ a1 g ss 11 þ a2 g 21 þ b1 u1 þ a1 g 12 þ a2 g 22 þ b2 u2 1 ;u 2 u
s:t:
8 ss max > 1 þ y1 gssd1 2 6 gg11 1 :u 2 P gg11 þ u ss 1 g ss 12 12 8 ss max 6 gg21 1 þ y2 gssd2 > u ss u > < 2 22 22 g ss ymin d2 pair2 21 2 P þ u u > 2 1 g ss g ss > 22 22 :
Fig. 3. Effect of bias noise on LP solution for case study 1.
ð10Þ
1 6 umax 6u umin 1 1 or measurement noise. The effect of noise applied to the LP bias with a steady-state value 1.0 on the LP solution is shown in Fig. 3. Since the bias steady-state value is 1.0 and the noise is mod 1 is constant, as expected from Fig. 2. From erate, the solution for u Fig. 2 we deduce that if the noise is negative (steady-state bias plus 2 is at its upper bound and the solution for noise is less than 1.0), u the output replicates the noise. If the noise is positive (noisy bias is greater than 1.0), the output is at its maximum value, while the 2 is amplified by 2:5 in the opposite direction to that solution for u of the noise. These results are verified by the simulation results in Fig. 3. We observe in this section that for MISO systems, the output set-point calculated by the LP is either constant or varies linearly with the bias parameter with gain of unity, assuming unique LP solutions for specified bias parameter values. However, the variation in the input target values may be larger than that of the bias. We show in the next section that for MIMO systems, the variation in output set-points may exceed that of the bias term due to interaction effects. 4. Sensitivity analysis: MIMO systems We focus in this section two-input, two-output systems. Key characteristics extend readily to the more general MIMO case. The LP optimization for the general two-input, two-output system takes the form: min
1 ;u 2 ;y 1 ; y 2 u
s:t:
1 þ a2 y 2 þ b1 u 1 þ b2 u 2 a1 y
ss 1 ¼ g ss y 11 u1 þ g 12 u2 þ d1 ss ss y2 ¼ g 21 u1 þ g 22 u2 þ d2
1 6 ymax ymin 6y 1 1
ð9Þ
2 6 ymax ymin 6y 2 2 1 6 umax umin 6u 1 1 2 6 umax umin 6u 2 2 A key difference between this formulation and the formulations considered previously is that two bias terms are now involved, and they are both updated at every iteration. Since optimization 2 relate to u 1 and u 2 through the model equality 1 and y variables y
2 6 umax umin 6u 2 2 The term a1 d1 þ a2 d2 in the objective function is omitted since it is constant for any given biases d1 and d2 , and does not affect the solution. From formulation (10), several inferences can be made regarding the feasible region of the problem and possible location of its 2 plane, the slope of the line which rep1 u solution. First, on the u resents the objective function (so-called ‘‘isocost line” [18]) does not depend on the bias terms. Its slope depends only on the cost coefficients and the steady-state gains of the process. Second, the output inequality constraints are presented in pairs, where each pair comprises two collinear constraints. The slope of each pair of constraints depends on the ratio of the steady-state gains, ss ss ss g ss i1 =g i2 ; i ¼ 1; 2. If the steady-state gains, g i1 and g i2 , are the same sign then the slope of the constraint is negative, and vice versa. The distance between the constraints in each pair does not depend on the bias term, and is only determined by the steady-state gains and the original output constraint limits. The distances between the constraints in each pair are D1 ¼
ymax d1 ymin d1 ymax ymin 1 1 ss ¼ 1 ss 1 g ss g g 12 12 12
D2 ¼
ymax d2 ymin d2 ymax ymin 2 2 ss ¼ 2 ss 2 : g ss g g 22 22 22
If d1 is an arithmetic mean of ymax and ymin (or d2 for the pair ymax 1 1 2 and ymin 2 ) then this pair of constraints is symmetric around the ori) or down (if di is gin. Otherwise it is shifted up (if di is closer to ymin i 2 versus u 1 plot (see Fig. 4). The feasible area ), on a u closer to ymax i for each pair of constraints lies between the constraint lines. Third, the minimum and maximum constraints for the input variables form a rectangle with the feasible area inside. The resulting feasible region of the problem is the intersection of all the constraint feasible regions and is schematically presented in Fig. 4. The solution of the optimization problem can be found by shifting the isocost line within the feasible region in the direction of decreasing cost until it cannot be improved any further. The resulting position of the isocost line will determine the optimal values of the optimization variables. At every iteration, the biases d1 and d2 are updated. If these values differ from their previous values, the shape of the feasible re-
20
A. Nikandrov, C.L.E. Swartz / Journal of Process Control 19 (2009) 16–24
Fig. 4. Graphical representation of the LP optimization problem for a 2 2 system.
gion changes, which may impact the solution of the optimization problem. As mentioned above, none of the slopes depends on the bias terms and the feasible region changes only because the constraint lines shift up and down. We consider here three scenarios in which the effect of bias variations on the optimal LP solution, hence MPC set-points, is explored. 4.1. Scenario 1 In this scenario, illustrated in Fig. 5, the optimal solution is deand ymax constraints. If the fined by the intersection of the ymax 1 2 combination of the steady-state gains is such that constraints are almost collinear, then the LP solution is expected to be sensitive to small changes in the biases. If one of the biases changes slightly, the new intersection of the shifted constraint lines (and, therefore, new solution) may appear some distance from the previous solution as shown in Figs. 5(a) and (b). However, the output set-points determined by the solution of the LP remain constant at the upper output constraint limits.
This situation is illustrated in Fig. 6. Analogously to the previous case, the solution of the LP can become sensitive to the changes in the bias term if the constraint lines are almost collinear. However, unlike the previous scenario, only one output set-point will be constant, and the other may fluctuate in response to variations in the bias value. In fact, we show in the case study below that it is possible for bias variations to be amplified in the set-point of the second output. Case study 2 We consider here a plant described by the following transfer function model, 1:2 1:4 u1 ðsÞ u2 ðsÞ 3s þ 1 3s þ 1 0:5 5:5 y2 ðsÞ ¼ u1 ðsÞ þ u2 ðsÞ 3s þ 1 3s þ 1 y1 ðsÞ ¼
ð11Þ
with constraints 1:0 6 y1 6 1:0 1:0 6 y2 6 1:0 0:5 6 u1 6 0:5
4.2. Scenario 2
0:5 6 u2 6 0:5
Here, the optimal solution is defined by the intersection of one 2 . of the output constraint limits and the maximum constraint on u
Maximization of f ¼ 0:3y1 þ 0:5y2 þ 0:6u1 þ 8:3u2 yields the corresponding LP:
21
A. Nikandrov, C.L.E. Swartz / Journal of Process Control 19 (2009) 16–24
Fig. 5. Effect of the bias on the LP solution: Scenario 1.
Fig. 6. Effect of the bias on the LP solution: Scenario 2.
1 þ 0:5y 2 þ 0:6u 1 þ 8:3u 2 max f ¼ 0:3y
1 ;y 2 ;u 1 ;u 2 y
s:t:
1 ¼ 1:2u 1 1:4u 2 þ d1 y 1 þ 5:5u 2 þ d2 2 ¼ 0:5u y 1 6 1:0 1:0 6 y 2 6 1:0 1:0 6 y 1 6 0:5 0:5 6 u 2 6 0:5 0:5 6 u
ð12Þ
Fluctuations around the steady-state bias values, d1 ¼ 0:4, d2 ¼ 1:8, results in the scenario depicted in Fig. 6. Simulations were conducted in which bias d2 was subjected to white noise, while bias d1 was kept at its steady-state value. The results of the simulation for manipulated variable (MV) targets and controlled variable (CV) set-points are presented in Fig. 7. Since the solution of the LP lies at the intersection of the constraints for y2 and u2 , the corresponding set-points and target
values are constant throughout the simulation period. Fig. 7(a) 1 with larger variashows that noise in d2 causes fluctuations in u 2 ¼ umax 2 is always equal to ymax ¼ 1 and u , tion. Indeed, since y 2 2 then the LP solution for any value of bias d2 must satisfy the equality: ss max 1 ¼ ymax ¼ g ss þ d2 ) u 2 21 u1 þ g 22 u2
1 d2 þ const: g ss 21
g ss 21 ¼ 0:5 and therefore, any changes in d2 cause changes in u1 that are twice as large and in the opposite direction. At the same time, 1 through the steady-state model equality 1 is related to u since y 1 2 and d1 are constants, any changes in u constraints in (12), and u 1 with coefficient g ss result in changes in y (which is 1:2 in this 11 case study). Summarizing both effects, it can be concluded that ss 1 such that: y 1 ¼ ðg ss noise in d2 results in changes in y 11 =g 21 Þd2 . In ss this case study g ss 11 =g 21 ¼ 2:4, and the effect on y1 is confirmed
22
A. Nikandrov, C.L.E. Swartz / Journal of Process Control 19 (2009) 16–24
Fig. 7. 2 2 system sensitivity: Case study 2. (a) Biases and input target values. (b) Biases and output set points.
1 are ampliby the results presented in Fig. 7. Thus, variations in y 2 remains constant. fied, while y 4.3. Scenario 3 In this scenario, we consider the solution to lie at the intersection of an output and input constraint as shown in Fig. 8. Variation in the bias can again cause amplification in one of the output setpoints, depending on the values of the process gains. Consider the response of the optimal solution to variations in 1 and y 1 at their maximum constraints. the bias term, d1 , with u From (10), we would have 2 ¼ u
g ss ymax d1 1 11 max u1 þ 1 ss ¼ ss d1 þ const g ss g 12 g 12 12
2 depends on the value Therefore, the sensitivity of the solution for u 2 also result in non. The fluctuations in u of the steady-state gain g ss 12 constant set-points for y2 . If g ss 22 is larger than unity, then any per 2 will be amplified, inducing even larger variations turbations in u 2 , and vice versa. in y
5. Performance of two-level LP-MPC system We consider here the performance of a two-level LP-MPC control system with noise applied to the bias. Graphically, the system is shown in Fig. 9. Case study 3 This case study is an extension of case study 2. It comprises the same plant, constraints and LP objective function, but we consider the closed-loop LP-MPC cascade system. The controlled variables are y1 and y2 , and the MPC parameters are as follows: Q ¼ I, R ¼ 0, S ¼ I, P ¼ 50, M ¼ 2, with sampling period, T s ¼ 0:3. The steady-state bias values were achieved by introducing output step-like disturbances passing through first-order filters with unity steady-state gains. Also, white noise was applied to output y2 , while output y1 was maintained uncorrupted. A perfect plant model was assumed. The results of the two-level control system operation after the transient effects of the step disturbance had dissipated are presented in Fig. 10. White noise with variance 0.0013 caused fluctu1 with variance 0.0074. 1 with variance 0.0051 and in y ations in u
Fig. 8. Effect of the bias on the LP solution: Scenario 3.
A. Nikandrov, C.L.E. Swartz / Journal of Process Control 19 (2009) 16–24
23
6. Conclusions
Fig. 9. Two-level LP-MPC control system.
This confirms that in this case study, white noise in the second output channel may cause fluctuations in the first output set-point with much larger variation. Even though the set-points for first output are corrupted by the noise, the control performance for the corresponding output is quite good. This could likely be ascribed to the relatively short control move horizon, as control move severity tends to attenuate with decreasing M. The variance of the output (0.000493) is much smaller than that of its set-point. However the performance of output y2 is quite poor (its variance is 0.0018) considering that its set point is constant at steady-state. This can be explained by fact that fast changes in yset 1 cause frequent changes in u1 which in turn affects y2 . u2 operates near its upper constraint and becomes saturated for periods of time during which it cannot be used for control.
Industrial implementations of MPC typically include a steadystate optimization layer to track the economic optimum and to provide feasible set-points to the model predictive controller. While two-level LP-MPC systems are relatively widely applied, reports of excessive set-point variation motivate the investigation of potential causes of such behavior. The present study focused on the sensitivity of the LP solution to variations in the LP model bias term, since the bias is updated prior to each LP execution based on the difference between measured and predicted plant outputs. A sensitivity analysis was conducted first for multi-input, single output systems. The steady-state LP model has only one bias term in this case, and its effect on the LP solution is conveniently summarized graphically. A key observation is that if the output is not at one of its limits, it varies linearly with the bias term with gain of unity, provided that the LP solution is unique. However, changes in the inputs may exceed those of the bias term if the gain relating the output to the corresponding input is less than unity. MIMO systems were explored through a detailed analysis of a two-input, two-output system. A graphical analysis was developed by expressing the problem in terms of the inputs only. Scenarios in which the solution is defined by two output constraints, and an output and an input constraint, were considered. In the former case, the set-points remain constant for small perturbations in the bias term. However, in the latter case, variations in the bias may result in amplified variations in the output that is not at its bound, depending on the magnitudes of the process gains. Simulation of the two-level LP-MPC system was conducted to evaluate the response of the two-level closed-loop system to a step disturbance combined with measurement noise. The results of the sensitivity analysis were used to explain the observed performance of the two-level system.
Fig. 10. Two-level control system response: Case study 3.
24
A. Nikandrov, C.L.E. Swartz / Journal of Process Control 19 (2009) 16–24
The focus of this paper has been on analysis of the sensitivity of LP-generated set-points to variations in the bias term induced, for example, by measurement and/or process noise. A question that naturally arises is how adverse amplification effects can be avoided or mitigated. One possibility would be the use of an alternative control structure. Another option might be to filter the bias update; however, the impact on the overall performance of the two-level LP-MPC system would need to be evaluated. A systematic investigation of these and other factors within a design framework would be a potentially useful avenue for future research. References [1] C.R. Cutler, B.L. Ramaker, Dynamic matrix control – a computer control algorithm, in: AIChE 86th National Meeting, Houston, TX, 1979. [2] J.F. Forbes, T.E. Marlin, Model accuracy for economic optimizing controllers: the bias update case, Ind. Eng. Chem. Res. 33 (1994) 1919–1929. [3] C.E. Garcia, A.M. Morshedi, Quadratic programming solution of dynamic matrix control (QDMC), Chem. Eng. Commun. 46 (1986) 73–87. [4] C.E. Garcia, D.M. Prett, M. Morari, Model predictive control: theory and practice – a survey, Automatica 25 (3) (1989) 335–348. [5] D.E. Kassmann, T.A. Badgwell, R.B. Hawkins, Robust steady-state target calculation for model predictive control, AIChE J. 46 (5) (2000) 1007–1024. [6] D.J. Kozub, Controller performance monitoring and diagnosis. Industrial perspective, in: Fifteenth Triennial World Congress of the International Federation of Automatic Control, Barcelona, Spain, 2002.
[7] J.H. Lee, M. Morari, C.E. Garcia, State-space interpretation of model predictive control, Automatica 30 (4) (1994) 707–717. [8] S. Li, K.Y. Lim, D.G. Fisher, A state space formulation for model predictive control, AIChE J. 35 (2) (1989) 241–249. [9] J.M. Maciejowski, Predictive Control with Constraints, Springer, 2002. [10] M. Morari, J.H. Lee, Model predictive control: past, present and future, Comput. Chem. Eng. 23 (1999) 667–682. [11] K.R. Muske, Steady-state target optimization in linear model predictive control, in: Proceedings of ACC, Albuquerque, NM, 1997. [12] D.M. Prett, M. Morari, The Shell Process Control Workshop, Butterworths, Stoneham, MA, 1987. [13] S.J. Qin, T.A. Badgwell, A survey of industrial model predictive control technology, Control Eng. Practice 11 (2003) 733–764. [14] C.V. Rao, J.B. Rawlings, Steady states and constraints in model predictive control, AIChE J. 45 (6) (1999) 1266–1278. [15] N.L. Ricker, Model predictive control with state estimation, Ind. Eng. Chem. Res. 29 (1990) 374–382. [16] S.L. Shah, R. Patwardhan, B. Huang, Multivariate controller performance analysis: Methods, applications and challenges, in: Chemical Process Control VI, Tuscon, Arizona, 2001. [17] R.C. Sorensen, C.R. Cutler, LP integrates economics into dynamic matrix control, Hydocarbon Process. (September) (1998) 57–65. [18] W.L. Winston, Operations Research: Applications and Algorithms, fourth ed., Duxbury Press, 2004. [19] C.-M. Ying, B. Joseph, Performance and stability analysis of LP-MPC and QPMPC cascade control systems, AIChE J. 45 (7) (1999) 1521–1534. [20] C. Yousfi, R. Tournier, Steady state optimization inside model predictive control, in: Proceedings of ACC, Boston, 1991.