High-performance Embedded Model Predictive Control using Step Response Models*

16th IFAC workshop on Control Applications of Optimization 16th IFAC6-9, workshop on Control ApplicationsGermany of Optimization October 2015. Garmisch-Partenkirchen, 16th IFAC workshop on Control Applications of Optimization October 6-9, 2015. Garmisch-Partenkirchen, Germany Available online at www.sciencedirect.com October 6-9, 2015. Garmisch-Partenkirchen, Germany

ScienceDirect IFAC-PapersOnLine 48-25 (2015) 138–143

High-performance Embedded Model Predictive High-performance Embedded Model Predictive  High-performance Embedded Model Predictive Control using Step Response Models  Control Control using using Step Step Response Response Models Models 

D. K. M. Kufoalor ∗∗ L. Imsland ∗∗ T. A. Johansen ∗∗ ∗∗ D. D. K. K. M. M. Kufoalor Kufoalor ∗ L. L. Imsland Imsland ∗ T. T. A. A. Johansen Johansen ∗∗ ∗ Dept. of Eng. Cybernetics, Norwegian University of Science and Technology ∗ Dept. of Eng. Cybernetics, Norwegian University of Science and Technology ∗ Dept. (NTNU), O.S. Bragstads plass 2DUniversity N-7491 Trondheim, Norway. of Eng. Cybernetics, Norwegian of Science and Technology (NTNU), O.S. Bragstads plass 2D N-7491 Trondheim, (e-mails: kwame.kufoalor, [email protected]) (NTNU), O.S. Bragstads plass 2D N-7491 Trondheim, Norway. Norway. ∗∗ Dept. of Eng. (e-mails: kwame.kufoalor, [email protected]) Cybernetics, and Center for Autonomous Marine Operations (e-mails: kwame.kufoalor, [email protected]) ∗∗ of Eng. Cybernetics, and Center Autonomous ∗∗ Dept. and Systems NTNU. [email protected]) Dept. of Eng. (AMOS), Cybernetics, and (e-mail: Center for for Autonomous Marine Marine Operations Operations and Systems (AMOS), NTNU. (e-mail: [email protected]) and Systems (AMOS), NTNU. (e-mail: [email protected]) Abstract: This paper presents an efficient approach used to obtain high-performance control targets Abstract: This efficient used obtain control for an embedded MPCpresents scheme,an has approach a relatively lowto complexity compared to a Abstract: This paper paper presents anwhich efficient approach used tocomputational obtain high-performance high-performance control targets targets for an embedded MPC scheme, which has a relatively low computational complexity compared to field industrial framework. Thea results show prediction quality is an essential for anproven embedded MPCMPC scheme, which has relatively lowthat computational complexity comparedfactor to aa field proven industrial MPC framework. The results show that prediction quality is an essential factor to consider usingMPC the embedded MPC challenging industrial applications wherefactor step field provenwhen industrial framework. Thescheme results for show that prediction quality is an essential to when using the embedded scheme for industrial where response models used. recursiveMPC prediction implementation and applications a more straightforward to consider consider when are using the A embedded MPC schememodel for challenging challenging industrial applications where step step response are prediction model more straightforward implementation, the standard implementation are examinedand in aathis paper. The results response models models referred are used. used.toA Aasrecursive recursive prediction model implementation implementation and more straightforward implementation, referred are in paper. show that both implementations havestandard inherentimplementation properties that directly affect control in the implementation, referred to to as as the the standard implementation are examined examined in this this performance paper. The The results results show that both implementations have inherent properties that directly affect control performance in presence disturbances. The recursive implementation out as the bestcontrol choiceperformance for the embedded show thatof both implementations have inherent properties turns that directly affect in the the presence of disturbances. The recursive implementation turns out as the best choice for the embedded MPC due to its superior robustness and computational properties. The impact of the prediction model presence of disturbances. The recursive implementation turns out as the best choice for the embedded MPC robustness and properties. The of implementation on the control of a subsea compact separation process also studied. MPC due due to to its its superior superior robustness and computational computational properties. Theisimpact impact of the the prediction prediction model model implementation on the control of a subsea compact separation process is also studied. implementation on the control of a subsea compact separation processbyisElsevier also studied. © 2015, IFAC (International Federation of Automatic Control) Hosting Ltd. All rights reserved. Keywords: Embedded model predictive control, step response models, robustness to model errors. Keywords: Keywords: Embedded Embedded model model predictive predictive control, control, step step response response models, models, robustness robustness to to model model errors. errors. 1. INTRODUCTION 1. 1. INTRODUCTION INTRODUCTION Model Predictive Control (MPC) is well established as the preModel Predictive (MPC) is as preferred method for Control challenging constrained multivariable control Model Predictive Control (MPC) is well well established established as the the preferred method for constrained multivariable control problems, especially in the process industries. Most industrial ferred method for challenging challenging constrained multivariable control problems, especially in industries. Most MPC solutions incorporate several features, ad hoc problems, especially in the the process process industries.including Most industrial industrial MPC solutions incorporate several features, including ad fixes, solutions that contribute to achieve control perforMPC incorporate severalsatisfactory features, including ad hoc hoc fixes, contribute to satisfactory control performancethat for different applications (Qin and Badgwell, 2000; Mafixes, that contribute to achieve achieve satisfactory control performance for applications (Qin 2000; ciejowski, 2002). Due to the possibility of Badgwell, operator intervention mance for different different applications (Qin and and Badgwell, 2000; MaMaciejowski, 2002). Due to the possibility of operator intervention and the reliance conventional low-level controllers, critical ciejowski, 2002).on Due to the possibility of operator intervention and the reliance on low-level controllers, critical issues, as closed-loop stability guaranty and robustness and the such reliance on conventional conventional low-level controllers, critical issues, such as closed-loop stability guaranty and robustness to unknown disturbances, are not the main focus in many issues, such as closed-loop stability guaranty and robustness to unknown disturbances, are main focus industrial applications (Maciejowski, when to unknown disturbances, are not not the the2002). main However, focus in in many many industrial (Maciejowski, However, when consideringapplications safety-critical applications,2002). Johansen (2015) points industrial applications (Maciejowski, 2002). However, when considering safety-critical applications, (2015) out the fact that although safe operationJohansen may be accounted for considering safety-critical applications, Johansen (2015) points points out the fact that although safe operation may be accounted for by building resilient architecture around an MPC scheme, the out the fact that although safe operation may be accounted for by building resilient architecture around MPC effects of model uncertainty and unknown on the by building resilient architecture around an andisturbances MPC scheme, scheme, the effects of model uncertainty and unknown disturbances on control of performance are still fundamental effects model uncertainty and unknownlimitations. disturbances on the the control performance are fundamental limitations. control performance are still still (2002), fundamental limitations. According to Maciejowski Dynamic Matrix Control According to Maciejowski (2002), Dynamic Matrix (DMC) Cutler and Ramaker(2002), (1979) Dynamic and its derivative, DMCAccording to Maciejowski Matrix Control Control (DMC) and Ramaker (1979) DMCPlus, areCutler the most known and and usedits of derivative, the commercially (DMC) Cutler andwidely Ramaker (1979) and its derivative, DMCPlus, are most widely and the available MPC products, with particularly usage in the Plus, are the the most widely known known and used used of ofhigh the commercially commercially available MPC products, with particularly high usage in petrochemicals These products high use finite reavailable MPC sector. products, withMPC particularly usagestep in the the petrochemicals sector. These MPC products use finite step response prediction models thatMPC rely on simple use disturbance modpetrochemicals sector. These products finite step response models that simple els and prediction assumptions on future planton behaviormod(Qin sponse prediction models that rely rely ondisturbance simple disturbance disturbance models and assumptions on future plant disturbance behavior (Qin and Badgwell, 2000; Maciejowski, 2002). Another example of els and assumptions on future plant disturbance behavior (Qin and Badgwell, Maciejowski, 2002). Another example of industrial MPC2000; software packages is SEPTIC (Statoil Estimaand Badgwell, 2000; Maciejowski, 2002). Another example of industrial MPC packages is tion and Prediction Tool for Identification and (Statoil Control)EstimaStrand industrial MPC software software packages is SEPTIC SEPTIC (Statoil Estimation and Prediction Tool for Identification and Control) and Sagli (2003), which achieves good performance forStrand many tion and Prediction Tool for Identification and Control) Strand and Sagli which achieves good performance for applications by using step response models. Nevertheless, the and Sagli (2003), (2003), which achieves good performance for many many applications applications by by using using step step response response models. models. Nevertheless, Nevertheless, the the

 This work is funded by the Research Council of Norway (NFR) and Statoil  This work is funded by the Research Council of Norway (NFR) and Statoil through the PETROMAKS project 215684, and also by NFR,(NFR) Statoil,and andStatoil DNV  This work is funded by the Research Council of Norway through the PETROMAKS through AMOS project project 223254.215684, through the PETROMAKS project 215684, and and also also by by NFR, NFR, Statoil, Statoil, and and DNV DNV through the AMOS project 223254. through the AMOS project 223254.

use of alternative model representations (e.g. state-space moduse alternative representations (e.g. els) is possible formodel applications where a better prediction moderror use of of alternative model representations (e.g. state-space state-space models) is possible for applications where a better prediction propagation will lead to improved control performance. The els) is possible for applications where a better prediction error error propagation will to control performance. The above MPC software incorporate features propagation will lead leadexamples to improved improved control numerous performance. The above MPC software examples incorporate numerous features into a PC/server-based softwareincorporate package, innumerous order to facilitate above MPC software examples features into aa PC/server-based package, order to the of high software control performance different appliintoachievement PC/server-based software package, in in for order to facilitate facilitate the achievement of high control performance for different cations. The features include strategies for correcting inaccurate the achievement of high control performance for different appliapplications. The features include strategies correcting inaccurate predictions in some cases, the plantfor is able to tune cations. Theand, features include strategies foroperator correcting inaccurate predictions in the operator is tune parameters and, for compensating errors (Strand andto predictions and, in some some cases, cases,model the plant plant operator is able able toSagli, tune parameters for compensating model errors (Strand and Sagli, 2003; Qin and Badgwell, 2000; Maciejowski, 2002). parameters for compensating model errors (Strand and Sagli, 2003; Qin and Badgwell, 2000; 2002). 2003; Badgwell, 2000; Maciejowski, Maciejowski, 2002). on typical Due toQin the and limited computational resources available Due to the limited computational resources available typical ultra-reliable industrial hardwareresources (e.g. PLCs), all theon Due to the limited computational available onfeatures typical ultra-reliable industrial hardware (e.g. PLCs), all the features that are foundindustrial in industrial MPC(e.g. software cannot ultra-reliable hardware PLCs),packages all the features that are in software packages cannot always expected to fit intoMPC an embedded framework. that arebefound found in industrial industrial MPC software MPC packages cannot always be expected to fit into an embedded MPC framework. For this reason, the enhancement of MPC performance in a always be expected to fit into an embedded MPC framework. For this reason, the enhancement of MPC performance relatively simple framework is examined in this paper. In For this reason, the enhancement of MPC performance in intheaa relatively simple framework is examined in this paper. In existing MPC literature, two isdifferent approaches are used to relatively simple framework examined in this paper. In the the existing MPC literature, two different approaches are used to derive step response prediction models.approaches The approaches are reexisting MPC literature, two different are used to derive step response prediction models. The approaches are referred to as the standard and recursive implementations in this derive step response prediction models. The approaches are referred as and implementations in paper. Fundamental properties that distinguish the prediction ferred to to as the the standard standard and recursive recursive implementations in this this paper. Fundamental properties that distinguish the prediction strategies used in theproperties standard and implementations paper. Fundamental that recursive distinguish the prediction strategies used the and are analyzed. hardware-in-the-loop simulation restrategies usedInin inaddition, the standard standard and recursive recursive implementations implementations are analyzed. In addition, hardware-in-the-loop results of a subsea compact separation process aresimulation used to show are analyzed. In addition, hardware-in-the-loop simulation results of a subsea compact separation process are used to show the capabilities of using a simplified translation of a highsults of a subsea compact separation process are used to show the capabilities of using translation of performance MPC obtained from SEPTIC. the capabilities of configuration using aa simplified simplified translation of aa highhighperformance MPC configuration obtained from SEPTIC. performance MPC configuration obtained from SEPTIC. 2. PREDICTION MODEL FORMULATIONS 2. 2. PREDICTION PREDICTION MODEL MODEL FORMULATIONS FORMULATIONS The step response model: The N−1 The step step response response model: model: − i) + s(N)u(k − N) (1) y(k) = N−1 N−1 ∑ s(i)∆u(k s(i)∆u(k − i) + s(N)u(k − N) (1) y(k) ∑ s(i)∆u(k − i) + s(N)u(k − N) (1) y(k) = = i=1 ∑ i=1 can be used to describe the dynamic and static interactions i=1 can to dynamic and between the input u (alsothe known as manipulated variable or can be be used used to describe describe the dynamic and static static interactions interactions between the input u (also known as manipulated between the input u (also known as manipulated variable variable or or

Copyright IFAC 2015 138 Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © 2015, IFAC (International Federation of Automatic Control) Copyright IFAC responsibility 2015 138Control. Peer review© of International Federation of Automatic Copyright ©under IFAC 2015 138 10.1016/j.ifacol.2015.11.073

138 138 138

IFAC CAO 2015 October 6-9, 2015. Garmisch-Partenkirchen, D. K.Germany M. Kufoalor et al. / IFAC-PapersOnLine 48-25 (2015) 138–143

MV) and the output y (known as the controlled variable or CV) of an asymptotically stable single-input-single-output (SISO) system, and ∆u(k) = u(k) − u(k − 1). The coefficients s(i) reach constant values after N sampling periods, i.e. s(N + 1) ≈ s(N). 2.1 Standard implementation Based on (1), the future output trajectory can be predicted using y(k ˆ + j|k) =

N−1

∑

i= j+1

s(i)∆u(k ˜ + j − i) + s(N)u(k ˜ + j − N) (2)

i=1

where the prediction horizon is specified by { j = 1, . . . , Hp } (see e.g. Garcia et al. (1989)). The present and future input moves ∆u(·) are unknown, while the past input u(·) ˜ and past input moves ∆u(·) ˜ are known. The constant disturbance model v(k|k) = ym (k) −

N−1

˜ − i) − s(N)u(k ˜ − N), ∑ s(i)∆u(k

(3)

i=1

also known as a bias term, is a usual choice for integral action in industrial MPC schemes, and it provides output feedback using the current measurement ym (k). A measured disturbance variable d (or DV) can be incorporated in a similar way as the MVs. Since the future changes in disturbance are not always known at time k, it is usually assumed that ∆d(k + j) = 0, j ≥ 1. In (2), the known terms of y(k ˆ + j|k) define the free response: yˆ f (k + j|k) =

N−1

˜ + j − i) + s(N)u(k ˜ + j − N) + v(k|k) ∑ s(i)∆u(k

i= j+1

which represents the response at each point along the prediction horizon, if the future inputs remain the same as u(k ˜ − 1). The above development leads to a prediction model of the form j

y(k ˆ + j|k) = ∑ s(i)∆u(k + j − i) + yˆ f (k + j|k)

(4)

i=1

The model (4) can be extended to include measured disturbances and rewritten in a matrix-vector form ¯ Yˆ (k + 1) = Θ∆U(k) + Yˆ f (k + 1), (5) ˆ + 1|k) y(k ˆ + 2|k) . . . y(k ˆ + Hp |k)]T , where Yˆ (k + 1) = [y(k T ∆U(k) = [∆u(k) ∆u(k + 1) . . . ∆u(k + Hu − 1)] ,  

s(1)

(6)

˜ − N), ¯ d ∆D(k) ¯ d ∆d(k) + Ψ ˜ ¯ + ϒ¯ d d(k where D(k|k) =Θ ˜ ˜ − 1) ∆u(k ˜ − 2) . . . ∆u(k ˜ − N + 2)]T , ∆U(k) = [∆u(k ˜ − N + 1) . . . u(k ˜ − N + Hp )]T , u(k ˜ − N) = [u(k ¯ d = [sd (1) sd (2) . . . sd (H p )]T , Θ ϒ¯ = s(N), ϒ¯ d = sd (Nd ),   s(2)

can be computed recursively using ¯ Y˜ (k + 1) = A¯Y˜ (k) + B∆u(k) + B¯ d ∆d(k), (7a) ¯ ˜ y(k) ˜ = C Y (k), (7b) ¯ ˜ where Y (k) is the state vector, C = [1 0 . . . 0]1×N ,       0 1 0 ··· 0 s(1) sd (1) 0 0 1 . . . 0  s(2)   sd (2)      ¯  .. .. ¯ = .. .. .. . . ..  A¯ =  , B = , B    . d . . . . . . .      0 0 0 · · · 1 sd (Nd − 1) s(N − 1) 0 0 0 · · · 1 N×N

s(N)

sd (Nd )

A free response model equivalent to (6) can be obtained by considering the initial state Y˜ (k − 1): ¯ u(k ˜ − 1), Y˜ (k) = A¯Y˜ (k − 1) + B∆ ˜ − 1) + B¯ d ∆d(k (8a) ¯ ˜ y(k) ˜ = C Y (k), (8b)

and the free response (including bias correction) becomes ¯ d ∆d(k) + A¯pY˜ (k) + 1[ym (k) − y(k)], Yˆ f (k + 1) = Θ ˜ (9) ¯ ¯ where A p contains the first H p rows of A . Using (9), the recursive prediction model can also be written as (5). Similar to the standard formulation, the recursive model can be extended to MIMO systems using block matrices/vectors.

The bias calculations reveal fundamental differences between the recursive and standard implementations. Consider the current output prediction in the recursive implementation using (8), without DVs: y(k|k) ˜ = y(k|k ˜ − 1) + s(1)∆u(k ˜ − 1) (10) Note that only s(1) is used directly in (10) at time k, assuming y(k|k ˜ − 1) is known perfectly. If the standard model is used, y(k|k) ˜ =

N−1

∑ s(i)∆u(k − i) + s(N)u(k − N),

(11)

i=1

. . . s(N − 2) s(N − 1) . . . s(N − 1) 0  .. .. .. , . . . . s(H p + 1) s(H p + 2) . . . 0 0

s(3) ¯ = Ψ  ..

a new estimate of the output trajectory ˜ + 1|k + 1) . . . y(k ˜ + N|k + 1)]T , Y˜ (k + 1) = [y(k

3.1 Robustness to small truncation errors

and

¯ U(k) ˜ ¯ Yˆ f (k + 1) = Ψ∆ + ϒ¯ u(k ˜ − N) + D(k|k) + v(k|k)1,

A convenient choice of state variables that leads to a state-space representation can be found in Li et al. (1989); Lee et al. (1994). Given the current output trajectory estimate ˜ y(k ˜ + 1|k) . . . y(k ˜ + N − 1|k)]T , Y˜ (k) = [y(k|k)

3. INHERENT PROPERTIES AND DIFFERENCES

0

... ..  s(2) s(1)  . , ¯ = Θ ..  ..  .. . . . s(H p ) . . . s(H p − Hu + 1)

˜ − N), ∆D(k), ¯ d , d(k ˜ and ∆d(k), respectively. The step for Ψ response coefficients for the measured disturbance models are denoted sd (·), and Nd represents the number of coefficients used. Note that Yˆ f (k + 1) in (5) is the response when no further control moves are applied, and depends only on the state of the plant. A straightforward generalization to multiple-inputmultiple-output (MIMO) systems can be achieved using the superposition principle (see e.g. Maciejowski (2002)). 2.2 Recursive implementation

j

+ ∑ s(i)∆u(k + j − i) + v(k|k),

139

s(3) s(4) .. .

1 is an H p -dimensional vector with all elements equal to one. ¯ u(k ˜ Similar definitions for Ψ, ˜ − N), ∆U(k), and ∆u(k) apply 139

and the entire step response sequence, including s(N), is directly involved in computing y(k|k). An error in s(N) will enter the bias calculation (3) directly at sampling time k, even though all past inputs available at time k describe y(k|k ˜ − 1) perfectly. The above observation implies that an error in s(N) will start affecting (11), N − 1 samples before affecting (10). Since a

IFAC CAO 2015 140 D. K.Germany M. Kufoalor et al. / IFAC-PapersOnLine 48-25 (2015) 138–143 October 6-9, 2015. Garmisch-Partenkirchen,

small truncation error, which is such that s(N + 1) ≈ s(N), will affect mainly the predicted steady-state output value, the error introduced will be matched by the constant model v(k|k) for unknown disturbances, and therefore has no significant effect on prediction quality. However, if the approximated s(N) is used to extend the sequence of step response coefficients, different error effects may be incurred by the two prediction model implementations. In a dense/compact variant of (6), commonly found in the MPC literature (see e.g. Maciejowski (2002); Camacho and Bordons ¯ are replaced by s(N), meaning that (2007)), the zeros in Ψ cumulative errors, due to truncation, are introduced into the prediction model. This may result in poor predictions if not properly compensated for. On the other hand, extending the step response sequence with an approximation s(N) for the recursive ¯ Due to the structure of implementation implies extending B. ¯ only step disturbances on the predicted output values will be B, incurred in this case, and such disturbances can be compensated for using the simple bias update.

”memory”. It rather updates (or corrects) its initial state estimate, which consists of predicted current and future outputs. Therefore, an error in the initial estimate, once captured by the bias term, remains as part of the bias update throughout the time the controller is active. A simple example that clearly illustrates the difference discussed above is treated in Section 5. 3.3 Computational efficiency The differences in computational efficiencies were identified by Li et al. (1989), where the main focus was the difference in computational speed. It should be, however, emphasized that the difference in memory usage is also a key factor to consider for embedded hardware with limited memory/storage resources. Based on the direct comparison of the number of floating point operations involved (cf. (6) and (9)) and the resulting code size in a code generation framework, the recursive implementation is identified as the most promising way to implement step response models on embedded hardware. 4. MPC PROBLEM FORMULATION

3.2 Effect of errors in initial estimates At start-up, an MPC based on step response models will need a means of constructing (or estimating) the state of the plant based on available input history. Since the system is stable, an ideal start-up will involve keeping the inputs constant (i.e. ∆u˜ = 0) for at least N time steps, and assuming that measured disturbances also remain constant (i.e. ∆d˜ = 0) in the same duration. The system will then be at steady state when the controller is activated at k = 0, and the initial state variables can be set to T Y˜ (k − 1) = [ym (0) ym (0) . . . ym (0)] , for the recursive model, and T T ˜ ˜ ∆U(k) = [0 0 . . . 0] , ∆D(k) = [0 0 . . . 0] , T

u(k ˜ − N) = [u(−1) u(−1) . . . u(−1)] , ˜ − N) = [d(−1) d(−1) . . . d(−1)]T , d(k for the standard model. The above ideal initialization has consequences, if changes in inputs or disturbances did occur N sampling intervals prior to time k = 0. In that case, the actions of the controller may not be smooth at start-up or when a switch from manual to automatic control mode is made. If a wrong initial steady state is estimated, the error will be matched by the bias compensation used in both the standard and recursive formulations. However, the controller’s response will not always be the same for the two output prediction implementations. A difference is identified when there is an (unknown) error in the measured disturbance at k = 0. Since the measured disturbance is an external input which has an effect that propagates through the plant, a ”ramp-like” effect of the error will be observed on the output. This effect will be handled in the same way by the bias update of the two prediction model implementations. However, Hp samples after the error was first detected, the first element of the vector u(k ˜ − N) in (6), i.e. u(k ˜ − N + 1), will experience a step change, since the erroneous starting history will have been resolved/cleared. The step change will cause a new ”bumpy” disturbance for the standard model, whereas the recursive model will continue to produce ”smooth” predictions due to the following reason. The recursive model does not reconstruct the state of the plant at each sampling time k based on a finite input history or 140

4.1 Embedded MPC scheme The MPC problem is formulated as: Hp

min

∑

j=Hw

y(k + j|k) − ry (k + j)2Q¯ y +

Hu −1

∑

j=0

∆u(k + j)2P¯

+ ρ¯ h ε¯h + ρ¯ l ε¯l subject to ¯ ∆u ≤ ∆u(k + j) ≤ ∆u, u ≤ u(k + j) ≤ u, y − ε¯l ≤ y(k + j|k) ≤ y¯ + ε¯h , ε¯h ≥ 0, ε¯l ≥ 0, u(k + j) = u(k + j − 1) + ∆u(k + j),

(12a) (12b) (12c) (12d)

j

y(k + j|k) = ∑ s(i)∆u(k + j − i) + yˆ f (k + j|k),

(12e)

i=1

where j ∈ {Hw , . . . , H p }, Hw > 1, for the output constraints, and j ∈ {0, . . . , Hu − 1} for the input constraints. The reference is denoted by ry (k + j). The slack variables ε¯h , ε¯l , weighted by ρ¯ h , ρ¯ l > 0, relax the constraints of the CV to avoid infeasibility in case of large disturbances or prediction model errors. Nominal closed-loop stability can be achieved by an adequate choice ¯ and the horizon lengths H p and Hu . of the weights Q¯ y , P, Problem (12) can be converted into a sparse quadratic programming (QP) problem (see e.g. Kufoalor et al. (2014)), where structure exploiting solvers can be used to compute the solution. An alternative is to use the equality constraints for elimination to produce a more dense QP as presented in Kufoalor et al. (2015). Since the choice of QP solver is not the main focus of this paper, only the sparse formulation is used: min Y T QyY + ∆U T P∆U − 2T T QyY + ρhT εh + ρlT εl (13a) subject to E∆U ≤ e, FU ≤ f , (13b) GY − Mh εh − Ml εl ≤ g, εh ≥ 0, εl ≥ 0, (13c) KU −∆U = Λu(k ˜ − 1), (13d) Y −Θ∆U = Y f , (13e) where the vectors and matrices involved now cover MIMO systems according to the definitions in Kufoalor et al. (2014, 2015). T represents the reference vector, and it is straightforward to deduce the constraint matrices.

IFAC CAO 2015 October 6-9, 2015. Garmisch-Partenkirchen, D. K.Germany M. Kufoalor et al. / IFAC-PapersOnLine 48-25 (2015) 138–143

The paper focuses on Y f in (13e), which represents the free response prediction model. It is the only variable/vector that is different for a given prediction model implementation in the MPC problem, when either (6) or (9) is used. All the remaining matrices and vectors in (13) are the same for both prediction model implementations. The free response Y f is computed only once at each sampling time k, and the predictions enter (13e) as parameters. In order to reduce the number of decision variables, move blocking is implemented (i.e. MVs are fixed to be constant over several time-steps), and CVs are evaluated on only specified evaluation points on the prediction horizon.

outputs

0.2

y1,

0.1

y1,

y2 (recursive);

y2 (standard)

0 -0.1 0

5

10

15

20

25

30

20

25

30

20

25

30

25

30

bias

0.3 0.2 0.1 0 -0.1

0

0.4

4.2 SEPTIC MPC scheme and embedded MPC design

141

5 u (recursive),

10

15

input

u (standard)

0.2 0

The main features of SEPTIC are presented in Strand and Sagli (2003). The intention of this section is to highlight the features that are directly translated into the embedded MPC scheme in subsection 4.1, and also outline the significant differences. Limits can be defined for MVs, MV moves, and CVs in SEPTIC MPC. However, a quadratic slack variable implementation is used to soften constraints on CVs. Note that an ”exact penalty” slack variable implementation is used in (12), where only one slack variable per constraint is needed. In a quadratic or 1-norm formulation, a separate slack variable for every constraint at every CV evaluation point is normally used, and this may result in a larger QP problem. In addition to weights on control targets, the priority of each control target (including priorities on constraints) can be assigned explicitly in SEPTIC. SEPTIC MPC solves a sequence of steady-state quadratic programs to respect the specified control targets with as many of the high priority control targets as possible, and the steady-state targets are used as references for the dynamic optimization problem. This feature is not implemented in the embedded MPC scheme. In this work, SEPTIC is used to obtain the step response models and MPC configurations, which serve as a high-performance control target specification for the embedded controller. The CV evaluation points are automatically generated by SEPTIC based on information from the step response models, prediction horizon, and MV blocks specified by the user. The priority on each constraint specified in SEPTIC is also considered as a tuning parameter when calculating appropriate weights, ρh and ρl , for the embedded MPC’s slack variable implementation. 5. A SIMPLE EXAMPLE Consider the simple mass-damper-spring system:           1 4 0.4 x1 (k + 1) 8 6 x1 (k) = + u(k) + d(k) , x2 (k + 1) 6 0.6 10 −3 2 x2 (k) y1 (k) = x1 (k), y2 (k) = x2 (k). (14) The variable d represents a known (or measured) disturbance. The sampling time is 1s, and the system (14) is asymptotically stable. The state-space model (14) is used to simulate the plant, and the step response data obtained from (14) is used in the prediction models, where N = Nd = 13. The MPC problem (12), excluding slack variables, is used with H p = Hu = 12, u = −2.5, u¯ = 10, y = [0, −10]T , y¯ = [20, 10]T , Q¯ y = diag(1, 0), and P¯ = 1, and changes in ry are made for only y1 . The effect of errors in the initial estimates and the inherent differences discussed in Section 3.2 are illustrated in Fig. 1, using the MPC example described above. A unit error (unknown to 141

-0.2

0

5

10

15

measured disturbance d(k) 4 2 0 -2 0

5

10

15 time [s]

20

Fig. 1. Illustration of the effect of an error in the initial measured disturbance for the recursive implementation (blue) and the standard implementation (red). the controller) was simulated in the disturbance measurement d(k) at time k = 0, and the ideal initialization discussed in Section 3.2 was used. The output reference was set to ry1 = 0, and no other errors were simulated in the entire duration of the simulation. Note that both implementations (recursive in blue, and the standard implementation in red) produce the same control performance when no errors are present. It can be seen in Fig. 1 that the plots of both implementations are almost inseparable until H p − 1 = 11 time steps have elapsed. The controller with the recursive model continues to operate smoothly, while the standard model experiences a new ”bumpy” disturbance effect, as discussed in Section 3.2. The situation reaches its worst point for the standard implementation when the erroneous starting history is fully resolved, 12 time steps after the effect of the error was first detected (see bias in Fig. 1). 6. INDUSTRIAL CASE STUDY The case study in this section is used to verify the practical importance and benefits of the discussions made in sections 3, and the outcome of the simple example above. 6.1 Subsea compact separator and MPC problem setup The subsea compact separation process is well described in Statoil’s patent documentation (Høydal et al., 2013). The setup used in this paper is the same as that used in Kufoalor et al. (2014, 2015). Therefore, only relevant details for this paper are repeated here. The main objective of the MPC is to control the quality of fluid (i.e. gas volume fraction) in the gas (Gasout ) and liquid (Liqout ) outlets of the separator. It is also necessary to control two pressure variables (P1 and P2 ) around their operational points and within their safety limits, while respecting the physical limits of three control valves. The MPC problem setup includes 4 CVs, each with 10 evaluation points. 1 CV (Liqout ) has a high bound, 1 CV (Gasout ) has a low bound, and 2 CVs (P1 and P2 ) have both high and

Inlet flow

0.06 0.04 0.02 0

200 400 time [s]

50 0 0

600

Table 2. Control performance of the standard and recursive implementations compared to the highperformance control targets of SEPTIC MPC

Gas in inlet flow

100 DV2 [%]

DV1 · 100[%]

IFAC CAO 2015 142 D. K.Germany M. Kufoalor et al. / IFAC-PapersOnLine 48-25 (2015) 138–143 October 6-9, 2015. Garmisch-Partenkirchen,

200 400 time [s]

Mean Square Error of each CV

600

MPC test setup SEPTIC MPC Embedded MPC recursive Embedded MPC standard

Fig. 2. Measured disturbances used to simulate hydrodynamic slugging in the subsea compact separation process. Table 1. Number of step response model coefficients (Ni, j ) from the jth MV or DV to the ith CV

P1

P2

0.01 0.01 0.04

0.004 0.003 0.008

Liqout Gasout 3.51 2.96 2.38

0.23 0.22 0.29

2

CV1

Ni, j for each MV and DV CVi CV1 (P1 ) CV2 (P2 ) CV3 (Liqout ) CV4 (Gasout )

MV1 44 94 93 93

MV2 MV3 DV1 DV2 94 94 18 16

93 90 81 85

56 57 53 51

0

−2 0 1

55 56 44 49

CV2

−1 0

Step response models of lengths shown in table 1 were obtained using SEPTIC and a nonlinear simulator, where all the models are such that s(Ni, j + 1) ≈ s(Ni, j ). In order to obtain appropriate matrix-vector sizes in the MIMO step response prediction models, simple extensions were made based on s(Ni, j + 1) ≈ ¯ The resulting lengths are ¯ and B. s(Ni, j ) and the structure of Ψ Ni = 94 for the recursive and N j , defined as the largest N for MV j , for the standard implementation. The simple choice of Ni and N j clearly favors the standard implementation, computationally. However, all matrix-vector operations are prepared offline, where the recursive implementation benefits the most. The significant extensions made through the choice of Ni for the recursive model are also intended to illustrate the model’s robustness to small truncation errors, since an exact steady value, s(Ni, j + 1) = s(Ni, j ), is not used. Using N j also leads to a slight extension of the standard model. Moreover, the standard model is implemented dense in order to verify the robustness properties discussed in Section 3.1. The standard implementation is the same as that used in Kufoalor et al. (2014, 2015). 6.2 Hardware-in-the-loop simulation setup The hardware-in-the-loop simulation setup consists of a nonlinear process simulator (i.e. the plant) running on a PC and the embedded MPC application running on an ABB AC500 PM592-ETH PLC. The PLC has a 400 MHz MPC603e microprocessor, 4MB RAM for user program memory, and 4MB integrated user data memory. The PLC communicates with the plant through an OPC server and an Ethernet connection. Two embedded MPC PLC programs were developed, one based on the standard free response model formulation (6), and the 142

CV3

200

300

400

500

600

100

200

300

400

500

600

100

200

300

400

500

600

100

200

300

400

500

600

0

0.2

low bounds. There are 6 slack variables (1 slack on each CV bound), 3 MVs (control valves), each with 6 move blocking indices, an upper limit, a lower limit, and a rate of change limit. The evaluation points and move blocks are chosen within a horizon of 80. Two measured process disturbances (DVs) are also present, where a challenging hydrodynamic slugging scenario is simulated by the DV sequences shown in Fig. 2. The above MPC problem description is the same for both the embedded MPC and SEPTIC MPC. In addition, all the controlled variables and constraints have priorities that are implemented explicitly in SEPTIC MPC. The required sampling time is 1 s.

100

0

−0.2 0

0.05

CV4

0

−0.05 0

time [s]

Fig. 3. Bias updates for the compact separator’s CVs for both the recursive (——) and standard (——) implementations. other based on the recursive implementation (9). The embedded MPC program code is written in C, and it incorporates the QP solver code, which is the same for both implementations. The tailored interior-point algorithm of CVXGEN (Mattingley and Boyd, 2012), was used to obtain high solution accuracies, in double precision, and the code is compiled with gcc 4.7.0, -mcpu=603e, and optimization level set to -O1. 6.3 Control performance results The impact of the two prediction model implementations on separation quality and pressure control in the subsea separator is shown in table 2. The Mean Square Error values are used as the closed-loop control performance measure. Using the recursive implementation in the embedded MPC scheme produced the same (or slightly better) high-performance control targets as those obtained by SEPTIC MPC on a PC. The standard implementation, on the other hand, produced results that are slightly poorer, but close to the SEPTIC MPC performance. Due to the characteristics of the slugging scenario (i.e. the fast changing ramps of DVs in Fig 2), it is not straightforward to relate the differences observed in control inputs and outputs to the differences discussed in Section 3. However, a closer look at the bias updates, shown in Fig. 3, reveals the underlying differences that can be traced to the properties discussed in Section 3. Similar to the results of the simple example in Fig 1, the obvious differences in prediction quality (in table 2) are related to the differences seen in the bias update plots of Fig. 3. First, it is observed in Fig. 3 that the initial response is not perfectly steady for all CVs. This situation should be expected in practice, since the measured DV, MV, and CV values at startup only produce an approximation of the true steady-state of

IFAC CAO 2015 October 6-9, 2015. Garmisch-Partenkirchen, D. K.Germany M. Kufoalor et al. / IFAC-PapersOnLine 48-25 (2015) 138–143

7. CONCLUSIONS

Table 3. Memory usage and free response Y f prediction time for the subsea compact separation process on an ABB AC500 PLC. Memory and time Data memory [MB] PLC program size [MB] Y f prediction time [ms]

Recursive

Standard

0.20 1.08 0.18

0.20 2.16 4.02

143

the plant. As discussed in Section 3.2, the slight deviation from the true steady-state values results in a similar initial response in both prediction models. However, the aftermath of the initial response is not smooth for the standard implementation, and a ”jump” can be seen near the 100 time step, when the slightly erroneous starting history is resolved for all CVs. The second observation is the effect of small truncation errors. It can be seen from the first two bias plots of Fig. 3 that both prediction models continue in a similar fashion until about halfway between the 100 and 200 time steps, where the standard model starts to have problems in reconstructing the state of the plant based on the past control actions and disturbance history (cf. N j for the DVs). Recall that since the disturbance sequence is unknown to the controller, the future behavior of DVs are assumed constant. The same assumption applies to the constant bias update model used. When this assumption is no longer valid (i.e. from k = 100, and in a number of sections of the DV sequences in Fig. 2), the effects of the small truncation errors incurred in the extended step response data start to show up in the standard implementation. The recursive model continues more smoothly since the effect of the small truncation errors only lead to step disturbances on the predicted output values (see discussions in Section 3.1). The difference in prediction quality is consistent with the control performance results. 6.4 Memory usage and computational time Table 3 shows the memory usage and the computational times obtained on the PLC. The results of the standard implementation are the same as those in Kufoalor et al. (2014, 2015) when CVXGEN is used, and are therefore old results. The new results, using the recursive model, illustrate the benefits of switching to a more computationally efficient implementation. Although the same solver is used in both implementations, the recursive implementation does not depend on the data structures of CVXGEN when computing Y f . Rather, the tailoring procedure exploits the sparse structure of the recursive model, and uses sparse matrix-vector operations to prepare the computations offline. 50% reduction in code size is achieved, and a speed-up of at least ×20 is obtained for the Y f computations.

The speed-up achieved in the recursive implementation does not necessarily lead to a similar gain in the controller’s total computational time. In this case study, CVXGEN’s solution time is the dominating factor. CVXGEN requires 18 iterations in the worst case, leading to a total computational time of 85.5ms for the standard implementation and 81.8ms for the recursive implementation. However, it is shown in Kufoalor et al. (2015) that when the qpOASES solver is used, a total time below 4ms can be obtained for the same MPC setup, where the standard implementation of Y f is used in a condensed version of (13). In that case, the time for computing Y f dominates, and the speed-up in the recursive implementation will become the most significant part of the controller’s total computational time.

143

It is shown in this work that it is possible to achieve high control performance by using an embedded MPC scheme that has low computational complexity. The quality of predictions is identified as an important factor which affects the control performance. In particular, it is shown that the choice of step response prediction model directly affects both control performance and the efficient use of computational resources of an embedded hardware target. The results show that MPC based on the recursive model is more robust to initialization errors and small truncation errors, compared to the standard implementation. The recursive implementation is also more computationally efficient, and the sparse structure can be exploited to obtain MPC code sizes suitable for embedded hardware. ACKNOWLEDGEMENTS We thank G. O. Eikrem, M. Fredriksen, and S. Strand at Statoil, for useful discussions on step response models and implementation aspects of SEPTIC. We also thank Statoil for making SEPTIC MPC tools and models available for this study. REFERENCES Camacho, E.F. and Bordons, C. (2007). Model Predictive Control. Springer. Cutler, C.R. and Ramaker, B.L. (1979). Dynamic matrix control–A computer control algorithm. In AIChE 86th National Meeting. Houston, TX. Garcia, C., Prett, D., and Morari, M. (1989). Model predictive control: theory and practice – a survey. Automatica, 25(3), 335–348. Høydal, J., Kristiansen, O., Eikrem, G.O., and Fjalestad, K. (2013). Method and system for fluid separation with an integrated control system. Patent nr. WO2013091719 A1. Johansen, T.A. (2015). Toward Dependable Embedded Model Predictive Control. IEEE Systems Journal. Kufoalor, D.K.M., Binder, B.J.T., Ferreau, H.J., Imsland, L., Johansen, T.A., and Diehl, M. (2015). Automatic Deployment of Industrial Embedded Model Predictive Control Using qpOASES. In European Control Conference. Linz, Austria. Kufoalor, D.K.M., Richter, S., Imsland, L., Johansen, T.A., Morari, M., and Eikrem, G.O. (2014). Embedded model predictive control on a plc using a primal-dual first-order method for a subsea separation process. In 22nd IEEE Mediterranean Conference on Control and Automation. Palermo, Italy. Lee, J.H., Morari, M., and Garcia, C.E. (1994). State-space interpretation of model predictive control. Automatica, 30(4), 707–717. Li, S., Lim, K.Y., and Fisher, D.G. (1989). A state space formulation for model predictive control. AIChE Journal, 35(2), 241–249. Maciejowski, J.M. (2002). Predictive Control: with constraints. Pearson and Prentice Hall. Mattingley, J. and Boyd, S. (2012). CVXGEN: A Code Generator for Embedded Convex Optimization. Optimization and Engineering, 13(1), 1–27. Qin, S.J. and Badgwell, T.A. (2000). An overview of nonlinear linear model predictive control applications. In Nonlinear Predictive Control, 369–392. Springer-Verlag. Strand, S. and Sagli, J. (2003). MPC in Statoil – advantages with in-house technology. In Int. Symposium on Adv. Control of Chemical Processes (ADCHEM), 97–103. Hong Kong.