Model Predictive Control with Integral Action: a Simple Mpc Algorithm

A simple Model Predictive Control (MPC) algorithm of velocity (incremental) form is presented. The proposed MPC controller is insensitive to slowly varying system and measurement trends and therefore has integral action. The presented algorithm is illustrated by both simulations and practical experiments on a quadruple tank MIMO process.

MPC has its origin in the theory of optimal control problems, Pontryagin et al. (1956) and Bellmann (1957), and in lecture notes on optimal control theory at MIT in the early 1960-ties which resulted in the book Athans and Falb (1966).
Parts of the optimal control theory are believed developed during the Apollo program in the early 60-ties where there was a great focus on minimum time (time optimal control) and minimum fuel (fuel optimal control) problems, see e.g.Athans and Falb (1966) Ch. 7 where the minimum principle by Pontryagin et al. (1956) is used to solve the optimal control problems.This also successful resulted in the state feedback optimal control law in the Apollo Lunar-Module Autopilot Widnall (1970).
The above mentioned optimal control problems are most often based on a fixed time optimization interval from say a constant initial time t 0 = 0 to a final time, say L, i.e. a time optimal control problem where a cost functional J(u) = L 0 dt = L and L is free, are optimized with respect to the control action u subject to a model and constraints.MPC problems are however based on a moving horizon (i.e. a receding horizon) strategy Garca et al. (1989) were we consider the initial time of the optimization interval to be equal to the present time t, and the final time of the optimization interval to be equal to the present time plus a prediction horizon, say t+L where L is the prediction horizon.There is no theoretical difference between standard optimal control problems and MPC problems and notice that the functional J(u) = t+L t dt = L.An early description of the moving horizon or receding horizon optimal control strategy where the initial integration time is set equal to the present time, t, and the final horizon is set equal to t + L where L > 0 is a constant time interval is as presented in Propoi (1963) and Balchen et al. (1970) p. 208 (in Norwegian).
An early survey paper on MPC is Garca et al. (1989).A survey of both linear and nonlinear MPC is given in, e.g.Maeder et al. (2009), Qin and Badgwell (2003).Some early work on non-linear MPC is presented in Balchen et al. (1992) and the references therein.
The Generalized Predictive Control (GPC) algorithm by Clarke et al. (1987) is an algorithm based on an input and output Controlled Auto Regressive and Integrated Moving Average (CARIMA) model.The GPC algorithm gives integral action of the closed loop system and the present state is estimated using the model and some old known input and output data.
Standard MPC algorithms usually do not achieve integral action and there is one main reason for this.The answer is that integral action is not necessarily optimal.However, if integral action is wanted there is some common methods to achieve this.
One commonly used method is to augment an integrator at the input, i.e. augment the plant model with an integrator u k = u k−1 + ∆u k .MPC algorithms is of state feedback type and in this case the MPC algorithm is a function of the plant state estimate xk and the estimate ûk−1 .The plant control is then u k = u k−1 + ∆u * k where ∆u * k is the MPC calculated optimal control deviation.One should notice that in this case and due to unknown disturbances the actual previous control u k−1 is not necessarily equal to the estimate ûk−1 .This method is among others described in Mayne et al. (2000).
Another commonly used strategy to incorporate integral action is to augment an integrator model d k = d k−1 of a constant disturbance d influencing the plant model state equation and the output equation.This strategy may be viewed as putting the integrator at the output.This strategy is also among others described in Mayne et al. (2000).
Recently in Ruscio (2012) a simple Linear Quadratic (LQ) optimal controller algorithm with integral action is presented.In this method constant or slowly varying disturbances in the state and output equations are removed from the problem by working on deviation models.
The contributions of this paper may be itemized as follows: • In this paper an MPC algorithm with integral action, along the same lines as used in the LQ optimal controller with integral action method, in Ruscio ( 2012) is presented.

Problem formulation
Given a process model where x k ∈ R n is the state vector, u k ∈ R r is the control input vector, y k ∈ R m is the output (measurement) vector, A, B and D are system matrices of appropriate dimensions, and x 0 is the initial state.
The disturbances v and w may both be unknown, i.e., v may be an unknown constant or slowly varying process disturbance, and w may be an unknown constant or slowly varying measurements noise vector.v and w may represent trends or drifts.
Note that the variables u k and y k in the model eqs.
(1) and ( 2) are the actual input and output variable, respectively.Furthermore, note that the model eqs.
(1) and (2) may arise from linearizing non-linear models around some nominal steady state variables.The model may also rise from system identification based on trended input and output data.Hence, in these cases, the external noise variables v and w are known, but the resulting control algorithm to be presented in this paper is insensitive to these noise variables v and w.Furthermore the system and the measurements may be influenced by drifts and in these cases the noise variables v and w may be unknown and slowly varying.Hence, the model eqs.( 1) and ( 2) is a realistic model.
We will study the MPC controller subject to the following performance index, where is the control rate of change (or control increment), r k is a reference signal and Q i and P i are symmetric positive semi-definite weighting matrices of appropriate dimensions.For finite prediction horizon L, then Q L may be chosen as the solution to Riccati equation of the problem to ensure closed loop nominal stability.
The above MPC objective criterion may be written in more compact form as where Q ∈ R Lm×Lm is a block diagonal matrix with Q i ∀ i = 1, . . ., L on the block diagonal.P ∈ R Lr×Lr is defined similar with P i ∀ i = 1, . . ., L on the block diagonal.The notation used to define the vectors in eq. ( 4) is defined in Appendix A.
In this paper we consider input rate of change and amplitude constraints.These constraints may be formulated as a linear inequality, where the matrix A and the vector b k are defined later in Sec.3.2.The MPC problem is now equivalent to a Quadratic Programming (QP) problem, i.e., the objective function eq. ( 4) with the process model, eqs (1) and (2), is minimized with respect to the unknown vector of future control increments, subject to the process constraints eq. ( 5), i.e., ∆u * k|L = arg min The simplified MPC strategy of including a control horizon, 1 ≤ L u ≤ L, and instead calculating a reduced number of future controls ∆u k|Lu , will be discussed and solved later in the paper.

Model discussion
Note that the model eqs.
(1) and (2), when v and w are constant vectors, is not unique.The constant trends v and w may be incorporated in the model by including one additional state.We find that the following model is equivalent with the initial state vector as and with z 0 = 1 in order to take care of the constant trends.

Problem solution
In order to solve the MPC optimal control problem eq. ( 6) we need a model which is independent of the unknown disturbances.For the sake of generality we are focusing on state space modeling.
From the state eq.( 1) we have where ∆x k = x k − x k−1 .From the measurement eq.
(2) we have Augmenting eqs.( 10) with ( 11) gives the state space model Hence, we have a strictly proper state space model of the form The state space model eqs ( 12), ( 13) (or equivalently ( 14), ( 15) ) may be used to define a Prediction Model (PM) of the form where and where O L is the extended observability matrix of the pair Ã, D and H d L ∈ R mL×(L−1)r the Toepliz matrix of impulse response matrices D Ãi−1 B ∈ R m×r .See Appendix A for definitions.
The performance index eq.( 4) with the PM eq. ( 16) can be written as a quadratic function on standard form, i.e., where The constant term J 0 in eq. ( 19) is not a function of the unknown ∆u k|L and then not needed and therefore not presented.
Notice that when the constraints in the MPC problem eq. ( 5) is inactive, then the unconstrained MPC controls is given by ∆u Usually we have process constraints and this will be discussed in the next Sec.3.2.

Constraints
It make sense to specify input rate of change constraints, i.e.
Using the relationship where S ∈ R rL×rL and c ∈ R rL×r are matrices with ones and zeroes as defined in the Appendix A.
We find that the constraints may be written as the linear matrix inequality where and The MPC algorithm is then to minimize the objective eq. ( 19) with respect to the constraints given by the linear inequality eq. ( 26).This is a standard Quadratic Programming (QP) problem in terms of the unknown future control increments and the optimal solution is ∆u * k|L , as defined in eq. ( 6).A receding horizon strategy is used and only the first control increment ∆u * k in the calculated ∆u * k|L , is used for control.The actual control action to the process is then

Reducing the number of unknowns future control actions
Usually when presenting MPC algorithms a control horizon is defined, and this control horizon is usually less than the prediction horizon L in order to reduce the number of unknown and then reducing the computation time.

Computing only the present control action
Consider now the extreme case in which the future control actions are equal to the present control action, i.e. such that u k+i−1 = u k ∀ i = 1, 2, . ... In this case the only unknown control action is u k , and equivalently ∆u k .
In this case we have that ∆u k+i = 0 ∀ i = 1, 2, . ... This gives and hence we have the control objective In this case we furthermore have a simple PM of the form where and where O L is the extended observability matrix of the pair ( D, Ã).The term p L is unchanged and given by eq. ( 17).This gives the control objective as a function of the increment ∆u k only, i.e. as follows where In this case we find that the constraints may be written as the linear matrix inequality where The MPC algorithm is then to minimize the objective eq. ( 33) with respect to the constraints given by the linear inequality eq. ( 36).This is a standard quadratic programming problem in terms of the unknown and the optimal minimizing solution is ∆u * k .A receding horizon strategy is used.The actual control action to the process is, The strategy presented in this section is considerably reducing the computational time of the MPC algorithm.This strategy is demonstrated to work considerably well for the control of the quadruple tank process.

Including a control horizon
In order to reduce the number of unknown input variables it is common to include a control horizon, say L u where 1 ≤ L u ≤ L. In this section we study the objective Hence, the compact form of this objective can be written as where P Lu = P (1 : rL u , 1 : rL u ).
In this case we find a PM of the form as in eq. ( 31) with where F L is defined in eq. ( 18).We are using MAT-LAB notation in eq. ( 41) and hence F is defined from all rows in F L and the first rL u columns.Substituting this PM into the objective eq. ( 40) gives where The constraints are then formulated as the linear inequality where in this case, for 1 ≤ L u ≤ L we have Hence we have a quadratic programming problem in the unknown vector ∆u k|Lu ∈ R rLu of future control actions and only the first vector, ∆u k is used.

Constant references and large prediction horizon
We will here discuss a special case which leads to a particular simple solution to the optimal control problem.
Consider the case where the references are constant, i.e., and that the prediction horizon is large or infinite.
Then we may use the LQ index where we have defined ỹk = y k − r.
If r is a non-zero constant reference then the measurements eq. ( 11) can be written as The state and output eqs.( 12) and ( 13) can then be rewritten as Hence, we have a strictly proper state space model of the form The state space model ( 50) and (51) (or equivalently ( 52) and ( 53)) with the performance index (48) defines a standard LQ optimal control problem.The optimal control is of the form which can be rewritten as The LQ optimal controller (55) gives y = r in steady state since the closed loop system is stable due to the properties of the LQ optimal controller.The states are seldom measured in practice.In this case we can use a state observer to define the deviation state ∆x k .However, another solution is to define ∆x k in terms of some past and known outputs . . ., y k−1 , y k and some known inputs . . ., u k−1 and the model matrices A, B and D.

Numerical examples
Example 6.1 (4 tank process (simulations)) Consider the quadruple tank process, Johansson (2000), with the non-linear state space model derived from mass balances and Bernulli's/Torricelli's law.By equating the potential energy and kinetic energy, i.e. mgh = 1 2 mv 2 and solving for the velocity we obtain v = √ 2gh.Multiplying with the area, a, of the outlet hole of the tank we obtain the volumetric flow-rate, q, out of the tank as q = av = a √ 2gh.
Hence, a mass balance of the four tank process gives the state space model where A i ∀ i = 1, . . ., 4 is the cross-section area of tank i, a i ∀ i = 1, . . ., 4 is the cross-section area of the outlet pipe of tank i.
The flow k 1 u 1 from pump 1 may be divided into a flow γ 1 k 1 u 1 into tank 1 and a flow (1 − γ 1 )k 1 u 1 to tank 4, i.e. such that the flow from pump number 1 is k 1 u 1 = γ 1 k 1 u 1 + (1 − γ 1 )k 1 u 1 .Similarly, the flow k 2 u 2 from the second pump may be divided into a flow γ 2 k 2 u 2 into tank 2 and a flow (1 − γ 2 )k 2 u 2 into tank 3.Here γ 1 and γ 2 are fixed parameters.The system is non-minimum phase when choosing these parameters such that, 0 < γ 1 + γ 2 < 1, and the system is minimum phase when, 1 < γ 1 + γ 2 < 2. The numerical values for the above parameters, as well as nominal values for the states and control inputs, are chosen as presented in Johansson (2000).
The 4 tank process is studied in a number of papers, see e.g.Gatzke et al. (2000) where Internal Model Control (IMC) and Dynamic Matrix Control (DMC) were used.Here we use the proposed MPC controller with integral action as presented in Sec.3.1.
The results after using the MPC controller in Sec. 3 in order to control the non-linear model eqs.( 56)-( 59) are presented in Figures 1 and 2. The MATLAB quadprog.mfunction is used to solve the QP problem as described in Sec.3.2.The weighting matrices were chosen simply as P = I 2 and Q = 0.0001I 2 .Only the minimum phase case is illustrated.

Experimental results on a quadruple tank process
A practical experiment with the MPC algorithm were performed on a quadruple tank process.The quadruple tank process is further discussed in Example 6.1.The sampling rate in all experiments is one second.We started with empty tanks in all experiments.Hence, this may be viewed as a test for robustness for unknown non-linearities when using the proposed MPC controller.The quadruple tank process setup results in a non-minimum phase behavior.
The experiments are described in the following items.
1.The non-linear model, Eqs. ( 56)-( 59) with measurements of the levels h 1 and h 2 in the lower tanks, were used in en Extended Kalman Filter (EKF) to estimate the four states in vector x needed in the MPC algorithm in order to formulate the term p L in the PM.
The experiment is illustrated in Fig. (3).In order to reduce the computational time using LabView the simple MPC algorithm in Sec 4 were used.From Fig.
(3) we see that this MPC algorithm performs well.
2. An open loop input experiment is designed as illustrated in Fig. ( 4) and the corresponding outputs, i.e. the levels in the two lower tanks, also illustrated in Fig. (5).
3. The input and output data are collected into data matrices U ∈ R N ×2 , and Y ∈ R N ×2 where the number of samples is N = 5459.The first N ID = 4000 first samples were used for identification.Hence, the last 1459 samples may be used for validation of the identified state space models.The data was also centered before use in the identification methods.5).The PE criterion evaluated for the validation data was V PEM = 3.38.
The best DSR model with n = 4 states was found with parameters L = 2 and J = 29.The simulated outputs are illustrated in Fig. ( 5).The PE criterion evaluated for the validation data was V DSR = 3.07.
7. Two SISO PI controllers were tuned by using the model based tuning method in Ruscio (2010).The model used was the DSR model.The experimental results using this decentralized control strategy is illustrated in Figs. ( 7) and ( 8).
8. The LQ optimal control strategy eq. ( 19) was implemented.The Kalman filter identified by the DSR method were used to identify the present state deviation ∆x k = x k − x k−1 needed in the controller.The experimental results using this LQ optimal controller with integral action strategy is illustrated in Figs. ( 7) and ( 8).
The conclusions drawn from these experimental results are discussed in the following.
Interestingly the identified state space models, both from PEM and DSR, fits the real data better than the FP model.Here the simulated output, i.e. the behavior from the input u, to the output y, is used in order to calculate the PE criterion.The results using the FP model, the PEM model and the DSR model are V fp = 7.57, V PEM = 3.38 and V DSR = 3.07 , respectively.Interestingly the DSR model is slightly better to fit the validation data compared to the PEM model.
Based on this conclusion we are using the identified DSR model for both tuning the PI controllers and for use in the LQ optimal controller with integral action strategy eq. ( 20).The deterministic part of the model, i.e. x k+1 = Ax k + Bu k and y k = Dx k , were used to tune the PI controller strategy (by first using the RGA pairing strategy, Bristol (1966)), as well as for the calculation of the feedback matrices G 1 and G 2 .Furthermore the DSR identified Kalman filter gain matrix K were used in the Kalman filter on deviation form as presented in Ruscio (2012), for estimating the deviation states ∆x k needed in eq. ( 20).
As we see from Figs. ( 7) and ( 8) the LQ strategy works very well compared to the PI controller strategy.This is justified by comparing the Integrated Absolute (IAE) indices.The DSR model gave IAE indices 1.6849 and 1.3290 for level one and two, respectively, and for the PI controllers 2.2723 and 2.5141 for level one and two, respectively.It is also worth mentioning that it is very difficult to tune PI controllers for this process due to the non-minimum phase behavior of the process.

Discussion
The presented MPC algorithm is based on a state space model of the plant and is therefore flexible to be used for MIMO systems.The algorithm may be combined with any state observer for estimating the present state, e.g. the Kalman filter, Jazwinski (1989) or simply with a state observer based on past inputs and outputs as described in Ruscio and Foss (1998).
The algorithm as presented in this paper is believed to work very similar as the GPC algorithm in Clarke et al. (1987).However, as mentioned in the introduction Sec. 1 the GPC algorithm is based on the use of an input and output CARIMA model.Such models are practical only for SISO systems.CARIMA models are capable of removing the influence of constant disturbances as in the state space model description in Eqs. ( 1) and (2) (SISO systems assumed).
The PM used in the GPC algorithm may be written on the form as in Eqs. ( 16) -( 18).Se e.g.Bitmead et al. (1990).One difference between the presented MPC algorithm and the GPC algorithm, is that the state estimate in the GPC algorithm, e.g. as in Eq. ( 17) is calculated based on the smallest number of past inputs and outputs.The presented MPC algorithm is more flexible with respect to state observers to be used.Se e.g.Ruscio and Foss (1998) for a state observer along these lines.

Concluding remarks
A simple state space MPC controller with integral action on velocity (incremental) form for MIMO systems is presented.

Acknowledgment
The author acknowledges the assistance of Mr. Danuskha who did the practical implementation on the quadruple tank process.

A. Notations used
The special structure of a Hankel matrix as well as some matching notations, which are frequently used throughout the paper, are defined in the following.Given a vector where nr is the number of rows in s k .Define integer numbers j and i and define the vector s j|i ∈ R inr as follows , which is defined as an extended vector.Illustrating the reference and the outputs from the process controlled by two single loop PI controllers, and the proposed LQ optimal controller with integral action.The LQ controller were constructed by using the DSR method for system identification.The DSR model was used to identify a Kalman filter for the system.The states were estimated with this Kalman filter and the deterministic part of the model were used to design the controller.Illustrating the reference and the outputs from the process controlled by two single loop PI controllers, and the proposed LQ optimal controller with integral action.The LQ controller where constructed by using the DSR method for system identification.The DSR model was used to identify a Kalman filter for the system.The states where estimated with this Kalman filter and the deterministic part of the model were used to design the controller.
The integer numbers j and i have the following interpretations: • j start index or initial time in the sequence used to form s j|i , i.e., s j , is the upper vector element in the extended vector s j|i .
• i is the number of nr-rows in s j|i .
Examples of such vector processes, s k , to be used in the above definition, are the measured process outputs, y k ∈ R m , inputs, u k ∈ R r and references, r k ∈ R m .
The extended observability matrix, O i , for the pair (D, A) is defined as where the subscript i denotes the number of block rows.
The reversed extended controllability matrix, C d i , for the pair (A, B) is defined as where the subscript i denotes the number of block columns.
The lower block triangular Toeplitz matrix, H d i ∈ R im×(i+g−1)r , for the quadruple matrices (D, A, B, E).
where I r is the r × r identity matrix and 0 r is the r × r matrix of zeroes.

Figure 1 :Figure 2 :
Figure 1: Simulation of the quadruple tank process and the minimum phase case in Example 6.1 with MPC control with integral action.Level references and actual levels illustrated.Prediction horizon L = 100 and control horizon L u = L.
4. A First Principles (FP) model, very similar to the one presented in Example 6.1, were fitted to the process as well as believed possible.Using the input experiment as illustrated in Fig. (4)) gave the simulated outputs as illustrated in Fig. (5).The Prediction Error (PE) criterion evaluated for the validation data was V fp = 7.57. 5.The MATLAB IDENT Toolbox system identification function pem.m were used to identify a n = 4 order state space model.The simulated outputs are illustrated in Fig. (

Figure 3 :
Figure 3: Practical run of the quadruple tank process as in Sec. 7, and the non-minimum phase case in Example 6.1 with MPC control with integral action.A control horizon L u = 1 and the algorithm in Sec. 4.

Figure 4 :
Figure 4: Open loop system identification input experiment, i.e. the volt input to the pumps.

Figure 5 :
Figure 5: This figure illustrates the real measurementsof the level in the two lower tanks as well as the corresponding simulated outputs of the system identification models, from DSR and PEM, as well as the simulated outputs from the first principles model.

Figure 6 :Figure 7 :
Figure 6: Quadruple tank process.Level in tank one upper and tank two lower.Illustrating the reference and the outputs from the process controlled by the proposed MPC controller, for three different cases.Case 1: FP model used in EKF to estimate the present state and linearized FP model in the MPC.Case 2: FP model used in EKF to estimate the present state and identified DSR model in the MPC.Case 3: DSR model and the corresponding Kalman filter used to estimate the present state and DSR model in the MPC.

Figure 8 :
Figure 8: Quadruple tank process.Level in tank one.Illustrating the reference and the outputs from the process controlled by two single loop PI controllers, and the proposed LQ optimal controller with integral action.The LQ controller where constructed by using the DSR method for system identification.The DSR model was used to identify a Kalman filter for the system.The states where estimated with this Kalman filter and the deterministic part of the model were used to design the controller.

I
i denotes the number of block rows and i+g −1 is the number of block columns, and where 0 m×r denotes the m × r matrix with zeroes.Define ∆uk = u k − u k−1 .Using that u k = ∆u k + u k−1 , u k+1 = ∆u k+1 + u k = ∆u k+1 + ∆u k + u k−1and so on, we haveu k|L = S∆u k|L + cu k−1 (63)where S ∈ R Lr×Lr and c ∈ R Lr×r are given by r 0 r • • • 0 r I r I r • • • 0 r .