Polyhedral Feasible Set Computation of MPC-Based Optimal Control Problems

Lantao Xie,Lei Xie,Hongye Su,Senior Member,IEEE,and Jingdai Wang

Abstract—Feasible sets play an important role in model predictive control(MPC)optimal control problems(OCPs).This paper proposes a multi-parametric programming-based algorithm to compute the feasible set for OCP derived from MPC-based algorithms involving both spectrahedron(represented by linear matrix inequalities)and polyhedral(represented by a set of inequalities)constraints.According to the geometrical meaning of the inner product of vectors,the maximum length of the projection vector from the feasible set to a unit spherical coordinates vector is computed and the optimal solution has been proved to be one of the vertices of the feasible set.After computing the vertices,the convex hull of these vertices is determined which equals the feasible set.The simulation results show that the proposed method is especially effcient for low dimensional feasible set computation and avoids the non-unicity problem of optimizers as well as the memory consumption problem that encountered by projection algorithms.

I.INTRODUCTION

M ODEL predictive control(MPC)is a widely used method in process industries and fields such as automobile,energy,environment,aerospace and medical treatment,etc[1],[2]because of its ability to effectively handle the complex dynamics of systems with multiple inputs and outputs,system constraints,and conflicting control objectives.A key role played in a MPC framework is the so-called feasible set,i.e.the largest subset of the state space for which there exists a control action satisfying all the constraints[3].MPC-based controllers will finally result in solving a standard optimal control problems(OCPs)such as linear programming(LP),quadratic programming(QP),semidefinite programming(SDP),dynamic programming(DP),etc.If the corresponding OCP is designed to be stable or recursively feasible,once the initial state is in the feasible set,there always exists an optimal control action at each following sample time.

As stated in[4],a larger feasible set usually means that the corresponding MPC algorithm is less conservative.Thus,when comparing conservativeness of different MPC algorithms,feasible sets are often compared,such as in[5],[6].Feasible sets can also be used to modify initial work states to ensure feasibility,and its shape description plays an important role in the effectiveness of many of the approaches for approximate explicit MPC[3].By restricting the state of the second step in the prediction horizon in the feasible set and computing each set iteratively,the algorithm proposed in[7]obtains the recursively feasible set and guarantees recursive feasibility.Even with the consideration of techniques used for stability,feasible sets are closely related with the prediction horizon and system constraints;generally longer horizons and looser system constraints result in larger feasible sets.

Oval-shaped and polyhedral feasible sets are two commonly used convex sets in MPC algorithms.Ellipsoid can be described as a weighted 2-norm and its size is related with some measure of the weighted matrix.By performing a maximum volume optimization in terms of the weighted matrix,one can obtain the largest ellipsoid,i.e.,the oval-shaped feasible set,such as that shown in[8].A polyhedron is the intersection of finitely many halfspaces.Computation of a polyhedron is much more complicated than of a ellipsoid.The standard approach to compute polyhedral feasible set is orthogonal projection[9].However,the orthogonal projection needs the polyhedron to be described by its vertices(V-representation)or the intersection of halfspaces(H-representation),which cannot be applied when the polyhedron is written in linear matrix inequality(LMI)form.Moreover,orthogonalprojection often turns out to be computationally intractable in high spatial dimensions.This often happens when the feasible set is computed for MPC-based OCPs with long prediction horizons or with large number of slack variables such as that in some robust MPC[10]and stochastic MPC algorithms[7].Alternatives to compute feasible set involve parametric linear programming(PLP)methods[11]and set relation methods[3].However,PLP methods may lead to numerical difficulties due to the non-unicity of the optimizers[12],while set relation methods are only efficient for OCPs without slack variables.

Based on the final OCP derived from MPC-based algorithms,this paper proposes a novel multi-parametric programming problem where the final OCP involves both spectrahedron(represented by LMI)and polyhedral(represented by a set of inequalities)constraints.According to the geometrical meaning of an inner product of vectors,the maximum length of the projection vector from the feasible set to a unit spherical coordinates vector is computed and the optimal solution is proven to be one of the vertices of the feasible set.After computing all the vertices of the feasible set,the convex hull of these vertices is computed which is equivalent to the feasible set.The simulation results show that the proposed method is especially efficient for low dimensional feasible set computation and avoids the non-unicity problem of optimizers,as well as the memory consumption problem encountered by projection algorithms.

Notations:The sets of reals is denoted by R.Given two integers a and b where a＜b,N[a,b]=[a,a+1,a+2,...,b]and given two sets X⊆Rnand Y⊆Rn,the Minkowski set addition is defined by X⊕Y={x+y|x∈X,y∈Y}and the Pontryagin set difference is then defined by X⊖Y={x|x⊕Y⊆X}.Given set Z⊆Rn+m,its projection onto Rnis defined by Projx(Z)={x ∈ Rn|∃y∈ Rmsuch that(x,y)∈Z}.Given matrix Qi,diag(Q1,...,Qn)denotes the block diagonal matrix with matrix Qion the main diagonal.E(·)denotes expectation.|→a|denotes the length of vector →a.[x]iis the i-th element of x.‖x‖ denotes the Euclidean norm.The semidefinite matrix A is denoted by A≿0.Given two vectors →a,→b,the angle between →a and →b is expressed as∠(→a,→b).conv{v1,v2,...,vn}denotes the convex hull of points(v1,v2,...,vn).

II.PROBLEM SETUP

Consider the controlled system to be described as

where xkis system state and ukis system input,△f is model uncertainty,dkis deterministic noise and wkrepresents stochastic noise with proper dimensions.Constraints to the systems include expectation constraints which can be expressed as HE(x)+GE(u)≥0[13],probability constraints(or chance constraints)(Pro{Hx+Gu≥ 0}≥ 1-ϵ)[14],and hard constraints(Hx+Gu≥0)[15].By using linear matrix inequalities(LMI)or constraints tightening methods,the optimal control problems of different algorithms can befinally converted to a uniform uncertainty-free form:

where≿denotes partial order[16](i.e.A≿B if and only if A-B is positive semidefinite)and H(xk,θk)is a square matrix.The symbol“≥”represents the vector inequality or componentwise inequality[17]when G(xk,θk)is a vector or a matrix,i.e.,A ≥ B means that aij≥ bij,for∀i,j.xkis the initial state and θkrepresents the inputs as well as the necessary slack variables depending on different systems and different MPC-based algorithms.

The feasible set of OCP(2)is the largest subset of the state space for which there exists a θksatisfying all the constraints that can be defined as

Def nition 1(Polytope[18]):A bounded convex polyhedron is a polytope or,equivalently,a set P is called a polytope if it can be expressed as the convex hull of finitely many points,i.e.,P=conv{x1,x2,...,xp}.A k-dimensional polytope is called a k-polytope.This means that some(k+1)-subfamily of(x1,x2,...,xp)is affinely independent,but no such(k+2)-subfamily is affinely independent.The convex hull of k+1 affinely independent points(i.e.,k-polytope)is a k-simplex.A simplex is a k-simplex if and only if it has k+1 vertices.

From this definition we can see that a 1-simplex is a closed segment,a 2-simplex is a triangle and a 3-simplex is a tetrahedron.We adopt the convention that the empty set is a polytope of dimension-1.

Assumption 1:The feasible set of OCP(2)is a(p-1)-simplex.

Thus,F can be represented as[19]:

where(v1,v2,...,vp)are the p vertices of F.Note that the description of a feasible set by a set of inequalities(H-representation)in(3)and the convex linear combinations(V-representation)in(4)are equivalent.The goal of this paper is to find all the vertices of F.

III.FEASIBLE SET COMPUTATION

As the value of inner product between two vectors →a,→b can be treated as the length of the projection vector from →b to →a if→a is a unit vector,the multi-parametric programming problem(mPPP)for computing the vertices of feasible set then can be expressed as:

where ‖c‖ =1.In fact,the restriction to the length of vector→c is not necessary as long as the cost function is to compute the maximum length of the projection fromto a vector in some direction.However,if we chooseto be unit,the cost function of mPPP(5)would be more meaningful as cTx==where →p is the projection fromto

We will start with the computation of a 2-dimensional feasible set.

A.2-D Case

When x∈R2,the feasible set can be shown in a plane.For convenience,we assume the p vertices(v1,v2,...,vp)of F are ordered.Here,“ordered”means viis adjacent to vi+1,vpis adjacent to v1,and the index increases anticlockwise.Let c=cα=(cosα,sinα),where c is the unit vector of the polar coordinator.Letdenote the boundary of F and let int(F)denote the interior of F,where F=⊕int(F).

Lemma 1:Given a point x ∈ F and a vectorif the angle between-andis equal or less than 90°for∀i∈ N[1,p],then x is on the boundary of polytope F.

Proof:If p≤2,it is obvious that x∈as=F.When p≥3,assume x is a interior point of F and xqis a point on the boundary of F withand λ ＞ 0.Assume xqis on the edge ofThen,there exists a i such that∠vjxvj+1+∠vj+1xvj+2+···+∠vj+i-1xvj+i≤ 180°and∠vjxvj+1+∠vj+1xvj+2+···+∠vj+ixvj+i+1＞ 180°where the subscript of v represents the remainder when dividing by p if it is larger than p.As[0°,90°],which means that∠vj+ixq,∠vj+i+1xq ∈ [0°,90°],∠vj+ixvj+i+1=360°-∠vj+ixq-∠vj+i+1xq≥ 180°.Since∠vj+ixvj+i+1is a interior angle of△vj+ixvj+i+1and less than 180°and conflicts with∠vj+ixvj+i+1≥ 180°,x cannot be a interior point of F.Thus,x∈ F.

Lemma 2:The solution x∗to problem(5)is on the boundary of F,i.e.,x∗∈

Proof:As x∗is the solution to problem(5),we have cx∗≥ cvi,∀i ∈ N[1,p].That isThus,the angle betweenis equal to or less than 90°for∀i∈ N[1,p].According to Lemma 1,x∗would be on the boundary of F.

Lemma 3:There exists a certain αi∈ [0°,360°]such that the solution x∗to problem(5)would be the vertice viof F for∀i∈ N[1,p].

Proof:From Lemma 2 we know that x∗∈Assume∃xα∈such that xα/=viand cαxα≥cαvifor∀α ∈[0°,360°].That is∃xα(/=vi)∈such thatfor∀α ∈ [0°,360°].Let α = α1such that90°andThen,let α2=180°+α1+ϵ where ϵ is a extremely small angle such thatThus,we cannot find a xα2∈such thatas Fig.1 shows a contradiction.

Fig.1. Feasible set in 2-D.

From Lemma 3 we know that,if we can find all the αi,we will obtain all the vertices of F.Thus,if we traverse α in[0°,360°]with a proper interval,the feasible set can be computed.This procedure can be described as

1)Let k=1.

2)For α =0°to 360°with step length ε (i.e.,α = α+ε for the next step).

3)Compute xαby solving J(x,cα)from OCP(5).

4)Let xk=xαand k=k+1.

5)Endfor.

6)Compute the convex hull of points set P =(x1,x2,...,xend)and let S=conv{P}.

Th eorem 1:There exists an appropriate ε such that S=F.

Proof:From Lemma 3 we get that there exists αi∈[0°,360°]such that the solution x∗to problem(5)would be the vertice viof F.ε is chosen such that αi=Niε,where Niis a proper integer for∀i∈ N[1,p].Thus,vi,∀i∈ N[1,p]would be in the points set P and S=conv{P}=F.

From the proof of theorem 1,we know that the optimal ε would be the maximum to make αi/ε an integer so that the total computation time would be reduced to a minimum.

B.3-D Case

If x ∈ R3,let c=cα,β=(sinβ,cosβ sinα,cosβ cosα)which is a unit vector in spherical coordinates.For vertice vi,we assume its adjacent ni vertices can be expressed as(vi+1,...,vi+ni).A similar lemma for the existence of these vertices can be described as:

Lemma 4:There exists a certain αi∈ [0°,360°]and βj∈[0°,360°]such that the solution x∗to problem(5)would be vertice viof F for∀i∈ N[1,p].

Proof:Assume ∃xα,β∈ F such that xα,β/= viand cα,βxα,β≥ cα,βvifor ∀α ∈ [0°,360°]and ∀β ∈ [0°,360°].That is ∃xα,β(/=vi) ∈ F such thatfor∀α ∈ [0°,360°]and ∀β ∈ [0°,360°].Let α = α1,β = β1such thatand90°,j=2,3,...,ni.Then,let α2=180°+ α1+ ϵ1,β2=360°- β1+ ϵ2where ϵ1and ϵ2are extremely small angles such thatfor some j(such j always exists;otherwise,F will be concave or vicannot be a vertice).Thus,we cannot find a xα2,β2∈ F such thatas Fig.2 shows this contradiction.

Fig.2.Feasible set in 3-D.

Again,if we traverse α and β through 0°to 360°,we will find all the needed vertices.This procedure can be described as:

1)Let k1=1,k2=1.

2)For β =0°to 360°with step length ε1

3) For α =0°to 360°with step length ε2

4) Compute xα,βby solving J(x,cα,β)from OCP(5).

5) Let xk1,k2=xα,βand k2=k2+1.

6) Endfor.

7) k1=k1+1.

8)Endfor.

9)Compute the convex hull of points set P=(x1,1,x1,2,...,xend,end)and let S=conv{P}.

Theorem 2:There exists an appropriate ε =(ε1,ε2)such that S=F.

Proof:Similar to Theroem 1,omitted here.

Parameters ε1and ε2can be chosen individually according to prior knowledge of the structure of the feasible set.If the vertices of the feasible set are smooth in some axis,the according ε can be chosen a little larger.The best ε would still be the maximum(ε1,ε2)to make αi/ε1and βj/ε2integers so that the number of circulation can be reduced to minimum.

C.High er Dimensional Case

When x∈Rnwhere n≥4 and considering the n-dimensional spherical coordinates,we specify c=(c1,c2,...,cn)where:

Theoretically,if we traverse (φ1,φ2,...,φn-1) from 0°to 360°with appropriate interval(ε1,ε2,...,εn-1),we will end up with a convex hull conv{P} =conv{x1,...,1,...,xend,...,end} with all the vertices vi∈conv{P},∀i∈ N[1,p].Thus,S=conv{P}just like in the 2-D and 3-D cases.However,if the computation time for every loop is T0with the majority comes from solving the mPPP(5),the total computation time Tt=(360°/ε0)n-1T0where we treat ε0= ε1= ···= εn-1will increase so rapidly as the dimension of x increases that the proposed method may not be tractable.

An alternative is to consider the orthogonal projection if the feasible set is well structured and the projection is computationally tractable.From the definition of orthogonal projection,we know that given P={x,θ|H(x,θ) ≿ 0,G(x,θ) ≥ 0},F=Projx(P).However,LMI H(x,θk)≿ 0 cannot be handled in projection algorithms in which only V-representation and H-representation are allowed.In fact,the region of LMI H={x,θk|H(x,θk) ≿ 0}is called a spectrahedron[20].A lot of work on the connections between polyhedrons and spectrahedrons has been done such as[20],[21].Recently,[22]proposed a method to convert the spectrahedron to an identical polyhedron.In particular,we have to find a transform matrix M such that MTH(x,θk)M=diag(Q(x, θk),D(x, θk))where D(x,θk))is a diagonal matrix map where[22]proved that H={x,θk|D(x,θk)≿ 0}.Thus,H can be expressed in H-representation H={x,θk|Dii(x,θk) ≥ 0}.Though[22]proposed an approach to compute M by seeking the joint invariant subspace and its orthonormal basis,this procedure is rather complicated.

Remark 1:The proposed method is particularly efficient for low dimensional feasible set computation.When the total computation time Ttbecomes unacceptable,the transformed orthogonal projection then should be considered.

IV.SIMULATION EXAMPLE

Example 1:Consider the widely used system[10],[23],[24]:

with system constraints x∈X={x|-10≤[x]1≤2,-10≤[x]2≤10},u∈U={u|-1≤u≤1}and w∈W={w|-0.01≤[w]1≤0.01,-0.01≤[w]2≤0.01}.Assume the applied algorithm is the min-max MPC proposed in[8]with open-loop control law.When the prediction horizon N=1,constraints to the mPPP would be

where S={x|0.5[x]2≤0.99,-0.66[x]1-1.33[x]2≤0.96,0.66[x]1+1.32[x]2≤0.96,0.43[x]1+0.21[x]2≤0.95,-0.43[x]1-0.21[x]2≤0.95}is the maximum robust control invariant set which can be computed by methods introduced in[4]and θ=u.Thus,the feasible set computed by the proposed method with ε=5°is shown in Fig.3 in dark green and set in light green is the set P={x,u|Ax+Bu∈X⊖S,x∈X,u∈U}.The result is identical to that computed by projection which means ε=5°works well.

Fig.3.Feasible set of Example 1 with N=1.

Example 2:Consider the continuously stirred tank reactor(CSTR)system in[7]:

where

Constraints for this CSTR system are Pr{-10≤[x]1≤10}≥90%,Pr{-5≤[x]2≤5}≥90%,Pr{-2.8≤[x]3≤2.8}≥90%,-2.156×10-2≤u≤0.2.Deterministic disturbances range-1≤[p]2≤1,-2≤[p]3≤2.[p]1is a persistent deterministic disturbance and[p]1=-1×10-3when t≥30 mins,otherwise,[p]1=0.Stochastic disturbance q∼N(0,10)and 10≤q≤10.

Using the proposed algorithms in[7],the final OCP involves both spectrahedral and polyhedral constraints.Choose ε0= ε1= ε2=3°and the feasible set can be computed accurately when N=12 as Fig.4 shows.The dimension of θ increases when the prediction horizon N increases which in turn increases,the total time to compute the feasible set.We record the total time for the proposed method as well as the orthogonal projection approach based on Fourier method[25]to compute the feasible set.From Fig.5 we can see that the computation time for the proposed method increases slowly as it depends on the computation time of the basic mPPP once the dimension of the feasible set is decided.However,the computation time as well as memory consumption of orthogonal projection increases rapidly when the dimension of θ increases,which ultimately leads to the memory exhaustion when N=10.Thus,when the dimension of θ is large and dimension of x is small,the proposed method would be a good choice to compute the feasible set.

The implementation was performed on laptop with Intel Core i7-4600U@2.10GHz 2.70GHz and 8 GB RAM.mPPP is solved by cvx tool box[26]in MATLAB and the polyhedron operation and orthogonal projection are realized by the multiparametric toolbox(MPT)[27].

Fig.4. Feasible set for example 2 with N=12.

Fig.5.Computation time.

V.CONCLUSION

In this paper,a multi-parametric programming problem has been proposed to compute the feasible set of the final OCP derived from MPC-based algorithms and involving both spectrahedral(represented by LMI)and polyhedral(represented by a set of inequalities)constraints.In fact,the forms or even the convexity and linearity of constraints are not restricted as long as assumption 1 holds.The simulation results showed that the proposed method is especially efficient for low dimensional feasible set computation and avoided the non-unicity problem of optimizers.Because the computation procedure just involves solving the basic mPPP repeatedly,memory consumption is not large.Thus,the proposed algorithm avoids the memory depletion problem encountered by projection algorithms.If we choose a large ε,a glimpse of the feasible set is available as the computed S would be a subset of the feasible set.This can be used to find a rough initial state or check the existence of a nonempty feasible set.

However,if the dimension of the feasible set is very large,the computational burden of the proposed method would be horrific.One may consider transforming the LMI-representation constraints to H-representation using transformation techniques from spectrahedron to polyhedron and then applying the orthogonal projection if the feasible set is well structured and the projection is computationally tractable.