Robust Model Predictive Control Using Linear Matrix Inequalities for the Treatment of Asymmetric Output Constraints

Abstract

One of the main advantages of predictive control approaches is the capability of dealing explicitly with constraints on the manipulated and output variables. However, if the predictive control formulation does not consider model uncertainties, then the constraint satisfaction may be compromised. A solution for this inconvenience is to use robust model predictive control (RMPC) strategies based on linear matrix inequalities (LMIs). However, LMI-based RMPC formulations typically consider only symmetric constraints. This paper proposes a method based on pseudoreferences to treat asymmetric output constraints in integrating SISO systems. Such technique guarantees robust constraint satisfaction and convergence of the state to the desired equilibrium point. A case study using numerical simulation indicates that satisfactory results can be achieved.

1. Introduction

Model-based predictive control (MPC) is a strategy in which a sequence of control actions is obtained by minimizing a cost function considering the predictions of a process model within a certain prediction horizon. At each sample time, only the first value of this sequence is applied to the plant, and the optimization is repeated in order to use feedback information [1, 2]. One of the main advantages of MPC is the possibility to consider explicitly the physical and operational constraints of the system during the design of the control loop [1–3]. However, if there are mismatches between the nominal model and the actual behavior of the process, then the performance of the control loop can be degraded and the optimization problem may even become unfeasible. Thus, the study of new strategies for design of robust MPC (RMPC) with guaranteed stability and constraint satisfaction properties even in the presence of uncertainties is an area with great potential for research [4–7].

In this context, Kothare et al. [7] proposed an RMPC strategy with infinite horizon employing linear matrix inequalities (LMIs) for dealing with model uncertainty and symmetric constraints on the manipulated and output variables. This approach was later extended to encompass multimodel representations [8, 9], setpoint management [10], integrator resetting [11], and offline solutions [12, 13].

Within this scope, Cavalca et al. [4] proposed a heuristic procedure that allows the inclusion of asymmetric constraints on the plant output, but without stability or constraint satisfaction guarantees. The present paper presents a formal strategy to handle asymmetric output constraints in the control of integrating single input, single output (SISO) systems, for the case where the output is linear in states. It is shown that robust constraint satisfaction is achieved, as well as convergence of the state trajectory to the desired equilibrium point. The effectiveness of the proposed method is illustrated by means of numerical simulations.

The remainder of this paper is organized as follows. Section 2 reviews the design of RMPC based on LMI. Section 3 formalizes the proposed technique for the inclusion of asymmetric constraints. Section 4 presents a case study consisting of a discrete time model of a double integrator. The results are evaluated through numerical simulations as presented in Section 5. Concluding remarks are presented in Section 6.

Throughout the text, 𝐼 represents an identity matrix, the notation (⋅∣𝑘) is used in predictions with respect to time 𝑘, and the ∗ superscript indicates an optimal solution. For brevity of notation, only the upper triangular part of symmetric matrices is explicitly presented.

2. LMI-Based RMPC

Consider a linear time-invariant system described by an uncertain model of the following form: 𝑥(𝑘+𝑦1)(𝑘=)𝐴=𝑥𝐶(𝑥𝑘)(𝑘+)𝐵,𝑢(𝑘),(1) where 𝑥(𝑘)∈ℝ^𝑛, 𝑢(𝑘)∈ℝ^𝑝, and 𝑦(𝑘)∈ℝ^𝑞 are the states, the manipulated and the output variables, respectively, for each time 𝑘, and 𝐴, 𝐵, and 𝐶 are constant matrices of appropriate dimensions. The uncertainty is represented in polytopic form; that is, matrices 𝐴 and 𝐵 are unknown to the designer, but they are assumed to belong to a convex polytope Ω with 𝐿 known vertices (𝐴_𝑖,𝐵_𝑖), 𝑖=1,2,…,𝐿, so that [2, 7]: (𝐴,𝐵)=^𝐿_𝑖=1𝜆_𝑖𝐴_𝑖,𝐵_𝑖(2) for some unknown set of coefficients 𝜆₁,𝜆₂,…,𝜆_𝐿 that satisfy ^𝐿_𝑖=1𝜆_𝑖=1𝜆_𝑖≥0,𝑖=1,…,𝐿.(3)

At each time 𝑘, the sequence of future controls 𝑈_∞={𝑢(𝑘∣𝑘),𝑢(𝑘+1∣𝑘),𝑢(𝑘+2∣𝑘),…} is obtained as solution of the following min-max optimization problem: min_𝑈∞max_{(𝐴,𝐵)∈Ω}𝐽_∞𝐴,𝐵,𝑈_∞,(4) where 𝐽_∞𝐴,𝐵,𝑈_∞=^∞_𝑗=0‖𝑥(𝑘+𝑗∣𝑘)‖^2𝑊_𝑥+‖𝑢(𝑘+𝑗∣𝑘)‖^2𝑊_𝑢(5) in which 𝑊_𝑥>0 and 𝑊_𝑢>0 are symmetric weighting matrices. By assumption, all states are available for feedback so that 𝑥(𝑘∣𝑘)=𝑥(𝑘).

This min-max problem can be replaced with the following convex optimization problem with variables 𝛾∈ℝ, 𝑄∈ℝ^𝑛×𝑛, Σ∈ℝ^𝑝×𝑛, and LMI constraints [2, 7]: subjecmtitno_𝛾,𝑄,Σ𝛾⎡⎢⎢⎣⎡⎢⎢⎢⎢⎢⎢⎣𝑄𝑄⋅𝑥01(𝑘)0⎤⎥⎥⎦>𝐴0_𝑖,((67))⋅𝛾𝐼0𝑊𝑄+_𝑥1/2𝐵_𝑖Σ⋅⋅𝛾𝐼𝑊_𝑢1/2𝑄⋅⋅⋅𝑄Σ⎤⎥⎥⎥⎥⎥⎥⎦>0,𝑖=1,2,…,𝐿.(8)

If the problem (6)–(8) has a solution 𝛾^∗𝑘,𝑄^∗𝑘,Σ^∗𝑘, then the optimal control sequence is given by 𝑢(𝑘+𝑗∣𝑘)=𝐾^∗𝑘𝑥(𝑘+𝑗∣𝑘),(9) where 𝐾^∗𝑘=Σ^∗𝑘𝑄^∗𝑘⁻¹.(10)

Symmetric constraints on the manipulated variables of the form |𝑢_𝑙(𝑘+𝑗∣𝑘)|<𝑢_𝑙, 𝑙=1,2,…,𝑝,𝑗≥0 and on the output variables of the form |𝑦_𝑚(𝑘+𝑗+1∣𝑘)|<𝑦_𝑚, 𝑚=1,2,…,𝑞,𝑗≥0 can be imposed by including additional LMIs [7] ⎡⎢⎢⎣𝑋⋅𝑄Σ⎤⎥⎥⎦>0(11) with ⎡⎢⎢⎣𝑄𝐴𝑋_𝑙𝑙<𝑢_𝑙2,𝑙=1,2,…,𝑝,_𝑖𝑄+𝐵_𝑖Σ^𝑇𝐶^𝑇𝑚⋅𝑦^2𝑚⎤⎥⎥⎦>0,𝑖=1,2,…,𝐿(12) for 𝑚=1,2,…,𝑞, where 𝐶_𝑚 denotes the 𝑚th row of 𝐶.

Let ℙ(𝑥(𝑘)) denote the optimization problem (6) with constraints (7), (8), (11), and (12). Suppose that the control law uses the concept of receding horizon; that is, 𝑢(𝑘)=𝐾^∗𝑘𝑥(𝑘) with 𝐾^∗𝑘 recalculated at each sampling time 𝑘. The following lemma is concerned with convergence of the state trajectory to the origin under the constraints imposed on the operation of the plant.

Lemma 1. If ℙ(𝑥(𝑘₀)) is feasible at some time 𝑘=𝑘₀, then ‖𝑥(𝑘)‖→0 as 𝑘→∞ with satisfaction of the input and output constraints.

The proof of Lemma 1 follows directly from the recursive feasibility and asymptotic stability properties demonstrated in [7].

Remark 2. As shown in Appendix A of Kothare et al. [7], ℙ(𝑥(𝑘)) is equivalent to minimizing a function 𝑉(𝑥(𝑘∣𝑘))=𝑥(𝑘∣𝑘)^𝑇𝑃_𝑘𝑥(𝑘∣𝑘) with 𝑃_𝑘=𝛾^∗𝑘(𝑄^∗𝑘)⁻¹>0. Such 𝑉(𝑥(𝑘∣𝑘)) function is found to be an upper bound for the cost function 𝐽_∞(𝐴,𝐵,𝑈_∞) in (5) and can be used as a candidate Lyapunov function in the proof of asymptotic stability.

Remark 3. For a regulation problem around a point different from the origin, a change of variables can be used, so that the new origin corresponds to the desired equilibrium point [7]. In the present work, the new value for the reference signal will be termed a pseudoreference. It is assumed that the process is integrating, and, therefore, the control value in steady state (𝑢_ss) is zero for any value of the pseudoreference. Otherwise, the determination of 𝑢_ss would not be trivial since the system matrices 𝐴 and 𝐵 are subjected to model uncertainties.

Remark 4. The RMPC problem formulation presented in this section considers that the matrices 𝐴 and 𝐵 are unknown but dot not vary with time. This is a particular case of the general framework introduced in [7], which was concerned with time-varying matrices 𝐴(𝑘) and 𝐵(𝑘).

3. Treatment of Asymmetric Constraints

The present work is concerned with regulation problems around the origin involving a SISO system with output variable 𝑦(𝑘)=𝐶𝑥(𝑘). This is a particular case of the problem described in Section 2, with 𝑝=𝑞=1. Therefore, the indexes 𝑙 and 𝑚 in (12) can be omitted. Moreover, the system is assumed to be integrating, so that 𝑢_ss=0 regardless of the pseudoreference, as discussed in Remark 3. It is also considered that the manipulated variable 𝑢(𝑘) is subjected to a symmetric constraint 𝑢 as in Section 2. Suppose that the constraints on 𝑦 are of the following form: 𝑦_min<𝑦(𝑘)<𝑦_max(13) with 𝑦_min<0 and 𝑦_max>0. If |𝑦_min|=𝑦_max, then the symmetric constraint formulation presented in Section 2 can be applied directly by making 𝑦=𝑦_max. If |𝑦_min|≠𝑦_max, a different approach is required. One alternative is to adopt a more conservative constraint, that is, 𝑦=𝑦^max−𝑦,_minif𝑦_max<||𝑦_min||,if𝑦_max>||𝑦_min||.(14)

However, this procedure may not be convenient in the following cases:(i)𝑦_max<|𝑦_min| with 𝑦_min≤𝑦(𝑘)<−𝑦_max,(ii)𝑦_max>|𝑦_min| with −𝑦_min<𝑦(𝑘)≤𝑦_max.

In these cases, the initial value of the output is admissible under the original asymmetric constraints, but not under the more conservative constraint (14). Figure 1(a) provides an illustration for case (i). As can be seen, the imposition of the more conservative constraint −𝑦_max<𝑦(𝑘)<𝑦_max makes the output variable 𝑦(𝑘) be located outside of the range of admissible values.

(a)

(b)

Figure 1 

(a) Illustration of the more conservative constraint (14) in case (i). (b) Heuristic solution proposed by Cavalca et al. [4].

Cavalca et al. [4] proposed a heuristic solution, based on a time-varying pseudoreference 𝑟(𝑘)=min{(𝑦_max+𝑦(𝑘))/2,0} as shown in Figure 1(b). The problem of asymmetric output constraints is then redefined in terms of new symmetric constraints (𝑎/𝑎) around 𝑟(𝑘). However, it should be noted that this technique does not lead to guaranteed stability and constraint satisfaction.

Unlike the approach described above, which involves a pseudoreference 𝑟(𝑘) that may change at each sampling time 𝑘, the solution proposed in the present paper employs a sequence of pseudoreferences 𝑟_𝑖(𝑖=1,2,…,𝑖_max) which are defined at 𝑘=0 on the basis of the initial output value 𝑦(0). As illustrated in Figure 2, symmetric constraints (𝑎/𝑎,𝑏/𝑏,𝑐/𝑐,…) are established around each pseudoreference. It will be shown that the use of such pseudoreferences, together with a convenient commutation rule, provides robust constraint satisfaction and ensures that the state trajectory converges to the origin. For sake of brevity of presentation, only case (i) will be treated. Case (ii) can be recast into case (i) by defining 𝑦(̆𝑘)=−𝐶𝑥(𝑘) and replacing constraints 𝑦_min≤𝑦(𝑘)≤𝑦_max with −𝑦_max≤𝑦(̆𝑘)≤−𝑦_min.

Figure 2 

Illustration of the proposed pseudoreference scheme.

Given an initial state 𝑥(0) such that 𝑦(0)=𝐶𝑥(0) falls within the scope of case (i), the following algorithm defines the pseudoreferences 𝑟_𝑖, as well as a sequence of associated matrices 𝑄^∗𝑖. These matrices will be subsequently employed in the control law to establish a rule of commutation from one pseudoreference to the next.

It is assumed that the set 𝑋_𝑠 of possible equilibrium values 𝑥_𝑠 for the state of the plant is known from the physics of the process to be controlled.

Algorithm for Determination of the Pseudoreferences and Associated Ellipsoids (PR Algorithm)
Step 1. Define the pseudoreferences 𝑟_𝑖 as follows:
Step  1.1. Let 𝑟₀=[𝑦_max+𝑦(0)]/2Step  1.2. Let 𝑖=1Step  1.3. While |𝑟_𝑖−1|>𝑦_max do        𝑟_𝑖=(𝑦_max+𝑟_𝑖−1)/2        𝑖=𝑖+1       End While Step  1.4. Let 𝑖_max=𝑖Step  1.5. Let 𝑟_𝑖max=0.Step 2. For each 𝑟_𝑖 determine 𝑥_𝑠,𝑖 such that: 𝐶𝑥𝑥_𝑠,𝑖=𝑟_𝑖,(15)_𝑠,𝑖∈𝑋_𝑠.(16)Step 3. Let𝜉₀=𝑥(0)−𝑥_𝑠,0𝜉_𝑖=𝑥_{𝑠,𝑖−1}−𝑥_𝑠,𝑖,1≤𝑖≤𝑖_max𝑦_𝑖=𝑦_max−𝑟_𝑖,0≤𝑖≤𝑖_max.Step 4. Solve the problem ℙ(𝜉_𝑖) with the constraints 𝑦_𝑖, 𝑢 for each 𝑖 (𝑖=0,1,…,𝑖_max) and denote the resulting solution by (𝛾̃^∗𝑖,𝑄^∗𝑖,Σ^∗𝑖).

The matrices 𝑄^∗𝑖 obtained in Step 4 define ellipsoids of the form (𝑥−𝑥_𝑠,𝑖)^𝑇(𝑄^∗𝑖)⁻¹(𝑥−𝑥_𝑠,𝑖)<1, as illustrated in Figure 3 for the case of a second-order system. It is noteworthy that, by construction (in view of LMI (7) with 𝑥(𝑘)=𝜉_𝑖), the ellipsoid 𝑖 contains the center of the ellipsoid (𝑖−1).

Figure 3 

Result of the PR algorithm.

The PR algorithm is said to be feasible if (15) and (16) in Step 2 and the optimization problem ℙ(𝜉_𝑖) in Step 4 are feasible for every 𝑖=1,2,…,𝑖_max. In this case, the resulting 𝑥_𝑠,𝑖,𝑦_𝑖,𝑄^∗𝑖,𝑖=0,1,…,𝑖_max are used in the control algorithm proposed below.

Algorithm for Control Using the Pseudoreferences (CPR Algorithm)
Initialization
Let 𝑘=0 and 𝑖=0.Step 1. Read the state 𝑥(𝑘).Step 2. If 𝑖<𝑖_max then  Let 𝑥̃_𝑖+1(𝑘)=𝑥(𝑘)−𝑥_𝑖+1,𝑠  If 𝑥̃^𝑇𝑖+1(𝑘)(𝑄^∗𝑖+1)⁻¹𝑥̃_𝑖+1(𝑘)<1 (condition of transition) then   Let 𝑖=𝑖+1  End If End If.Step 3. Let 𝑥̃_𝑖(𝑘)=𝑥(𝑘)−𝑥_𝑠,𝑖 and solve the problem ℙ(𝑥̃_𝑖(𝑘)) with the constraints 𝑦_𝑖 and 𝑢 in order to obtain (𝛾^∗𝑘,𝑄^∗𝑘,Σ^∗𝑘).Step 4. Calculate 𝐾^∗𝑘=Σ^∗𝑘(𝑄^∗𝑘)⁻¹ and 𝑢(𝑘)=𝐾^∗𝑘𝑥̃_𝑖(𝑘).Step 5. Apply 𝑢(𝑘) to the plant.Step 6. Let 𝑘=𝑘+1, wait a sample time and return to Step 1.

The main result of this work is stated in the following theorem, which is concerned with the satisfaction of constraints and convergence of the state trajectory to the origin under the control law given by the CPR algorithm.

Theorem 5. If the PR algorithm is feasible and the CPR algorithm is applied to control the plant, then ‖𝑥(𝑘)‖→0 as 𝑘→∞, with satisfaction of the input and output constraints.

Proof. For 𝑘=0, the state 𝑥(𝑘) lies in the ellipsoid associated with 𝑄^∗0, which was obtained by solving problem ℙ(𝑥(0)−𝑥_𝑠,0) in the PR algorithm. Therefore, problem ℙ(𝑥̃₀(0)) is feasible by hypothesis. Lemma 1 then guarantees that ‖𝑥̃₀(𝑘)‖→0 (i.e., ‖𝑥(𝑘)−𝑥_𝑠,0‖→0) as 𝑘→∞, with satisfaction of the input and output constraints, under the control law stated in Steps 3 and 4 of the CPR algorithm with 𝑖=0. Since the ellipsoid associated with 𝑄^∗1 contains 𝑥_𝑠,0, by construction, it can be concluded that the condition of transition stated in Step 2 of the CPR algorithm with 𝑖=1 will be satisfied in finite time. Let 𝑘₁ be the first time when this condition is satisfied, that is, 𝑥̃^𝑇1𝑘₁𝑄^∗1⁻¹𝑥̃₁𝑘₁<1.(17) This condition ensures that the optimization problem ℙ(𝑥̃₁(𝑘₁)) is feasible, since, by construction, the solution (𝛾̃^∗1,𝑄^∗1,Σ^∗1) of ℙ(𝜉₁) satisfies the constraints of ℙ(𝑥̃₁(𝑘₁)). In fact, condition (17) ensures that LMI (7) is satisfied with 𝑥(𝑘) and 𝑄 replaced with 𝑥̃₁(𝑘₁) and 𝑄^∗1, respectively. Moreover, the remaining LMIs (8),(11)-(12) are satisfied by 𝑄^∗1 and Σ^∗1 by hypothesis. Therefore, after switching from 𝑖=0 to 𝑖=1, Lemma 1 ensures that ‖𝑥̃₁(𝑘)‖→0 (i.e., ‖𝑥(𝑘)−𝑥_𝑠,1‖→0) as 𝑘→∞ with satisfaction of the input and output constraints. As a result, the condition of transition with 𝑖=1 will be satisfied in finite time. A similar reasoning can be applied to show that the condition of transition will be satisfied for all 𝑖=0,1,…,𝑖_max−1 in finite time. Finally, when 𝑖=𝑖_max, the state 𝑥(𝑘) will be inside the last ellipsoid, centered at 𝑥_𝑠,𝑖max=[00⋯0]^𝑇, and then Lemma 1 will ensure that ‖𝑥(𝑘)‖→0 as 𝑘→∞, again with satisfaction of the constraints.

4. Numerical Example

A discrete state-space model of a double integrator will be employed to illustrate the proposed strategy. The matrices of the model are given by𝐴=⎡⎢⎢⎣10𝑇1⎤⎥⎥⎦,𝐵=⎡⎢⎢⎣𝑇²𝑇/2⎤⎥⎥⎦𝜅,𝐶=10.2,(18) where 𝑇 is the sampling time, and 𝜅 is an uncertain gain parameter.

The initial condition is set to 𝑥(0)=[−100^𝑇], the sampling time is 𝑇=0.5 s, and the constraints are defined as 𝑢=𝑢_max=𝑢_min=5, 𝑦_min=−10 and 𝑦_max=0.1. The uncertain parameter 𝜅 is assumed to be in the range 0.9−1.1. The weight matrices of the controller are defined as 𝑊_𝑥=⎡⎢⎢⎣10000.001⎤⎥⎥⎦,𝑊_𝑢=100.(19)

In this case, the set 𝑋_𝑠 of possible equilibrium points is given by 𝑋_𝑠=𝑥₁,𝑥₂∣𝑥₂=0.(20)

In fact, the state variables 𝑥₁ and 𝑥₂ can be regarded as position and velocity, respectively, and thus the equilibrium can only be achieved if the velocity 𝑥₂ is zero.

The pseudoreferences 𝑟_𝑖 and associated symmetric output constraints 𝑦_𝑖, which define the controllers 𝑖=0,…,6, are presented in Table 1. The simulations were performed in the Matlab environment.

Table 1 

Pseudoreferences and associated output constraints.

5. Results and Discussions

As discussed in Section 3, a possible approach to address asymmetric constraints consists of adopting the conservative bounds defined in (14). In the present example, such procedure amounts to setting 𝑦=𝑦_max=0.1. However, the resulting optimization problem ℙ(𝑥(0)) becomes infeasible. In fact, given the control constraint −5<𝑢(𝑘)<5, it is not possible to steer the output from 𝑦(0)=−10 to the range [−0.1,0.1] in a single sampling period.

On the other hand, if the constraints are relaxed by setting 𝑦=−𝑦_min=10, there is no guarantee that the resulting output trajectory will remain within the original (𝑦_min,𝑦_max) bounds. In fact, the inset in Figure 4 shows that such a relaxation of the output constraints does lead to a violation of the original upper bound for 𝜅=0.9.

Figure 4 

Simulation results using relaxed output constraints.

These findings motivate the adoption of the proposed strategy for handling the asymmetric output constraints. Figures 5(a) and 5(b) present the simulation results obtained by using the CPR algorithm. As can be seen, both the output and control constraints were properly enforced, even by using the extreme values of 𝜅 in the simulation.

(a)

(b)

Figure 5 

Simulation results using the proposed strategy: (a) output and (b) control signals.

The commutation between the successive pseudoreferences is illustrated in Figure 6. This graph indicates that the commutation from one pseudoreference to the next occurs in finite time, as expected. The final pseudoreference (𝑖=6) corresponds to the origin, which is the desired equilibrium point for the regulation problem.

Figure 6 

Commutation between pseudoreferences.

6. Conclusion

This paper presented a strategy for handling asymmetric output constraints within the scope of an LMI-based RMPC scheme. For this purpose, a procedure for defining a sequence of pseudoreferences was devised, along with a rule for commutation from one pseudoreference to the next. The proposed approach guarantees constraint satisfaction and convergence of the state trajectory to the origin, provided that the algorithm for determination of the pseudoreferences is feasible. The results of a numerical simulation study indicated that the proposed procedure may be a suitable alternative to the use of either more conservative constraints (which may lead to infeasibility issues) or more relaxed constraints (which do not guarantee satisfaction of the original restrictions). Future research could be concerned with the extension of the proposed approach to multiple input-multiple output (MIMO) systems. In this case, it may be necessary to define different pseudoreferences for each constrained output under consideration.

Acknowledgments

The authors gratefully acknowledge the support of FAPESP (scholarship 2008/54708-6 and grant 2006/58850-6), CAPES (Pró-Engenharias), and CNPq (research fellowships).