Ellipsoidal Set Based Command Governors for Constrained Linear Systems with Bounded Disturbances

This paper provides a solution for a linear command governor (CG) that employs invariant and constraint-admissible ellipsoid. The motivation is to substitute the typical polyhedral set used in almost all CG schemes with the ellipsoidal one, which is much easier to construct. However the price for this ofﬂine computational efﬁciency is that the size of the feasible set can be relatively small, and the online computational burden is heavier than that of polyhedral set based CGs. The proposed solution overcomes these two weaknesses and offers a very attractive alternative to polyhedral set based CG. Two numerical examples with comparison to earlier solutions from the literature illustrate the effectiveness of the proposed algorithm.


I. INTRODUCTION AND PROBLEM FORMULATION
Consider the following linear discrete-time system subject to bounded disturbances x(k + 1) = Ax(k) + Bu(k) + Dw(k) (1) where x(k) ∈ R n is the measured state, u(k) ∈ R m is the command, and w(k) ∈ R l is the disturbance.
In CG, (1) generally represents the closed-loop system and reflects the combined closed-loop dynamics of the plant and of a stabilizing controller. Thus, it is assumed that all eigenvalues of A are strictly in the interior of the unit circle.
The state x(k) and command u(k) are subject to linear constraints, (x, u) ∈ L xu , with where g ∈ R q , g > 0. The matrices F x , F u are of appropriate dimension. The inequalities (2) are to be interpreted element-wise. The disturbance w(k) is known only to the extend that it belongs to the bounded polytopic set W W = conv{w j , j = 1, p} where the verticesw j , j = 1, p are known, and conv(·) denotes convex hull. The description (3) implies that w(k) satisfies, ∀k ≥ 0 where γ j (k), j = 1, p are unknown parameters that satisfy p j=1 γ j (k) = 1, γ j (k) ≥ 0 Define r(k) as an arbitrary reference sequence for u(k). In CG, the objective is to design a command u(k) = U (x(k), r(k)) such that, for a given reference r(k), the controlled system x(k + 1) = Ax(k) + BU (x(k), r(k)) + Dw(k) fulfills the command and state constraints (2) despite the uncertainties. Furthermore, if the reference r(k) is reachable, then u(k) should be steered to r(k); if r(k) is not reachable, then the controlled system should be steered to the closest possible steady controlled variable.
A typical procedure in CG consists of two phases: offline and online.
• Offline Phase: Under the condition that u(k + 1) = u(k), construct a robustly invariant and constraint-admissible set Ω of the closed-loop system, subject to the constraints (2), and ∀w ∈ W   x(k + 1) with .
• Online Phase: At each time instant k, calculate the command u(k) by solving the following optimization problem, Q 0 min u(k) (u(k) − r(k)) T Q(u(k) − r(k)) , The main idea of CG is the following. If there is no constraint violation, u(k) = r(k) so that the CG does not interfere with the desired operation of the system. If u(k) = r(k) would cause a constraint violation, the CG seeks the closest admissible command to the current reference r(k).
In the extreme case, thanks to the robustly positive invariance of Ω, u(k) = u(k − 1), which implies that the CG temporarily isolates the system from further variations of the reference to guarantee safety, in terms of constraint satisfaction.
In set invariance theory [1], ellipsoidal set and polyhedral set are widely considered. The popularity of ellipsoidal sets is due to computational efficiency via the use of linear matrix inequality (LMI) formulation. It is well know [1] that constructing a polyhedral Ω is generally harder than the computation of an ellipsoidal one.
For high dimensional systems, the computation of polyhedral invariant constraint-admissible set might be prohibitive, as it might takes many hours and may require more memory than available.
However to the best of the author's knowledge, in the CG literature [2], [3], polyhedral invariant constraint-admissible sets are uniquely employed. This is because (a) If Ω is a polyhedral set, then (8) is a quadratic programming (QP) problem. In the other case, (8) is a second order cone programming problem [4]. Hence the online computational burden with a polyhedral set is smaller than that with an ellipsoidal one.
(b) With linear state and input constraints (2), polyhedral invariant constraint-admissible sets are preferred to the ellipsoidal ones, as they offer a larger approximation of the domain of attraction.
The objective of this paper is to present a solution to overcome the two disadvantages (a), (b) of ellipsoidal set based CG. For (a), it is shown that the online optimization problem can be boiled down to a problem of finding a unique positive root of a univariate polynomial. Hence the online computational burden is dramatically reduced. For problem (b), we adopt the strategy that is proposed in [5] for enlarging the domain of attraction. The main idea is to introduce an extra degree of freedom through the use of a vanishing part, which is the output of a fictitious dynamics. As written in [6], there is a lack of design methods able to systematically design the fictitious dynamics. We aim at covering this lack by providing a convex formulation for the optimization of fictitious dynamics. Despite a significant improvements in the performance, the proposed approach does not induce any extra online computational load.
This paper is structured as follows. Section II contains results on ellipsoidal set based CG for the nominal case. Then in Section III the results are extended to the robust case with bounded disturbances. Two simulated examples with comparison to earlier solutions from the literature are evaluated in Section IV before drawing the conclusions in Section V.
Notation: A positive definite (semi-definite) matrix P is denoted by P 0 (P 0). 0 n×m is the zero matrix of dimension n × m, 0 n /I n is an n × n zero/identity matrix. We denote by R the set of real numbers, by R n×m the set of real n × m matrices , and by S n the set of positive semi-definite n × n matrices. For a given P 0, E(P ) represents the following ellipsoid For symmetric matrices, the symbol ( * ) denotes each of its symmetric block.

II. NOMINAL CASE
In this section, linear discrete-time systems of the form are considered. This case is addressed first in order to introduce the main ideas.

A. Reducing the Online Computation Burden
Here we show how to convert (8) to a problem of finding a unique positive root of a univariate polynomial with Ω being an ellipsoid, i.e., The result is a generalization of [7] where the case r = 0 was considered. For this purpose, P a is decomposed as with P xx ∈ S n , P ss ∈ S m and P ux ∈ R m×n . For simplicity, the time index k is dropped. Rewrite It is well known [4] that (12) is a convex second order cone program (SOCP). Note that if P uu = 0, then (12) becomes a QP problem. In general, QP can be solved much more efficiently than SOCP.
Using the Karush-Kuhn-Tucker condition [4], the optimal solution (12) satisfies the following with R uu = Q −1 P uu , R ux = Q −1 P ux . Note that by using the eigenvalue/vector decomposition of R uu , one can find the inverse of (I m + λR uu ) analytically.
For a given x, the constraint in (12) defines an ellipsoid, denoted as E u , within which u must lie. If r ∈ E u then u = r is the solution of (11). If r / ∈ E u , the minimizing u should lie on the boundary of E u , i.e., Substituting (13) in (14), one obtains a univariate polynomial of λ, say P λ , of degree 2m. This equation can have only two real roots, one positive corresponding to the shortest distance from the boundary of E u to r, the other negative corresponding to the largest distance. There is a number of well known techniques, e.g., Newton-Raphson method [8] that can be used to find the unique positive root of P λ with a quadratic rate of convergence to the solution. Once λ is found, the command u(k) is computed using (13). As consequence, the optimization problem (12) can be solved efficiently.

B. Enlarging the Domain of Attraction
This section provides the main contribution of the paper. We show how to enlarge the domain of attraction of the ellipsoidal set based CG. For this purpose, the command u(k) is decomposed as where v(k) is used to follow r(k), andȳ(k) is a vanishing part, which is the output of the whereĀ ∈ Rn ×n ,C ∈ R m×n are design parameters.
Combining (10), (15), (16) one gets, under the condition v(k + 1) = v(k), Using (2), (15), the following constraints are obtained DefineΩ as any invariant constraint-admissible set for (17), (18). At each time instant, the CG calculates v(k) andx(k) by solving the following optimization problem whereQ 0 is any matrix satisfyingĀ Denote v * (k) andx * (k) as the optimal solution of (19). The command u(k) is computed as It is shown [5] that the control scheme (19), (21) retains the same theoretical properties as the standard CG (8) such as constraint satisfaction, finite settling time, but has a larger feasible set and hence faster response. The price to be paid is the increased offline and online computational burdens whenΩ is a polyhedral set [5]. It will be shown later in this section that withΩ being an ellipsoid, the online computational burden remains almost the same as the standard CG (8).
For the moment, let us focus on the design parametersĀ,C.
Note that the idea of decomposing u(k) as (15) is well known [9], [5]. It is first introduced in [9], whereĀ = γIn with 0 < γ < 1,n = m were chosen. In [5], another choice consists in the use of shift registers. It was shown that with this choice, CG can be considered as a particular MPC scheme. Another popular choice uses Laguerre sequences [10] which possess orthogonality properties. As written in Introduction, the problem of designingĀ,C for the CG is still open. The aim of this section is to provide a convex formulation for optimizingĀ,C.
Rewrite (17) as We want to findĀ,C that give the largest possible domain of attraction on the x a −subspace.
This is done by looking for an invariant constraint-admissible ellipsoid E(P e ) = {z ∈ R (na+n) : , that the volume of the projection from E(P e ) onto the x a −subspace is maximized, It is well-known [11], [12] that invariance and constraint-admissibility of E(P e ) for (22) under with The problem of maximizing the volume of E xa can be recast as The optimization problem (26) is non-convex in P e ,Ā,C. It will be shown that (26) can be reformulated as a convex semi-definite (SDP) optimization problem by using a nonlinear parameter transformation technique. For this purpose, first define the matrices P e , P −1 e partitioned into blocks as where X ∈ S na , Y ∈ S na , U ∈ R na×n , V ∈ R na×n , and (?) denotes blocks in P e , P −1 e with no importance for the derivation to be presented in the sequel. Note that the condition P e P −1 e = I na+n implies Ifn < n a , then (28) imposes a non-convex rank constraint, since it implies that rank(X −Y ) =n.
Hence in the rest of the paper it is assumed thatn ≥ n a .
The following theorem holds.
where M =CV T , N = UĀV T .
It is well known [12] that (25) is satisfied if and only if ∃Π 0 such that Note that (34) is (31). Pre and post multiplying (33) respectively by yields (30). The proof is complete.
Since T T P −1 e T = Y , the projection of E(P e ) onto the x a −subspace is given as The maximization of the volume of E(Y −1 ) is performed by solving the following problem Denote X * , Y * , Π * , M * , N * as a solution of (35). Since (35) is a convex SDP optimization problem ∀n ≥ n a , it follows that Y * is independent ofn. In other words, there is no advantage to be gained by usingn > n a . Therefore,n = n a is considered in the rest of the paper.
Note that the cut of E(P e ) throughx = 0 is given by In the interest of the size of the domain of attraction the set E(Y −1 ) should be maximized, sine . On the other hand E(X −1 ) should also be maximized because this is the region on which the tracking performance is primary purpose.
To address both objectives, the following SDP optimization problem can be used for some tuning parameter µ ≥ 0.
To sum up, the method to design the fictitious dynamics is summarized in Algorithm 1. Since can be chosen arbitrary close to 1, the obtained invariant ellipsoid E(P e ) can be used as a reliable approximation of the domain of attraction.

Remark 2:
There are several U, V satisfying (28). In this paper the following solution Now we show how to extend the result in Section II-A in order to solve (19). Decompose P e as with P xx ∈ S n , P ss ∈ S m+na , P sx ∈ R (m+na)×n . Define Note that (41) has the same form as (12). Hence it can be solved effectively.

III. ROBUST CASE
In this section, we extend the preceding development to the robust case, namely we show how to enlarge the domain of attraction.
The following result is recalled. It is a slight modification of Theorem 2 in [13]. It concerns an LMI condition for invariance of an ellipsoid E(P ) = {x ∈ R n : x T P x ≤ 1} . The condition will be used to optimize the fictitious dynamics.
Lemma 1: If there exist P 0 and 0 < α < 1 such that, ∀j = 1, p hold, then E(P ) is robustly invariant for the following autonomous system, ∀w(k) ∈ W x(k + 1) = Ax(k) + Dw(k) (43) Remark 3: In Theorem 2 of [13], the disturbance w is assumed to be constrained in an ellipsoid as w T w ≤ 1. Hence Lemma 1 exploits the polytopic nature of w and avoids the bounding of the set W via an ellipsoid.
As in Section II, to enlarge the domain of attraction, the command u(k) is decomposed as where v(k) is used to follow r(k), andȳ(k) is the output of the following system whereĀ ∈ R na×na ,C ∈ R m×na ,D ∈ R na×l are variables to be optimized. Recall that n a = n+m.

Remark 4:
It should be noted that the fictitious dynamics of the CG in [5] are a special case of our scheme withD = 0 na×l .
Combining (1) x a (k + 1) x(k + 1) where x a (k) = x(k) T v(k) T T . Or equivalently, with Using [12] and Lemma 1, robust invariance and constraint admissibility of E(P e ) = {z ∈ R 2na : where Similarly to Section II, P e , P −1 e are partitioned into blocks as The following theorem holds.
Pre and post multiplying (49) by respectively yields (52). The proof is complete.
The cut of E(P e ) throughx = 0 is given by E(X −1 ) = {x a ∈ R na : x T a X −1 x a ≤ 1}. The projection of E(P e ) onto the x a − state space is given as E(Y −1 ) = {x a ∈ R na : x T a Y −1 x a ≤ 1}. As in Section II, one would like to maximize the size of E(X −1 ), and of E(Y −1 ) for the performance and for the domain of attraction. This can be done by solving the following optimization problem max for some tuning parameter µ ≥ 0.
Since α is a variable in the design of the fictitious dynamics, the constraints (52) are bilinear matrix inequalities (BMI) due to the multiplication terms αX, αY . One could solve directly the BMI problem (56). In this paper, an alternative procedure is considered. The idea is to perform a univariate search on 0 < α < 1. At each iteration, an SDP problem is solved for a fixed α.
The procedure stops when no more progress is made with respect to a given tolerance.

IV. EXAMPLES
This section demonstrates the potential benefit of the new methods by simulations of two examples system. The CVX toolbox [14] was used to solve the SDP optimization problems. For comparison purposes, the ellipsoidal set based CG without the fictitious dynamics, i.e.,Ā = 0 andC = 0, will be denoted as the standard CG.
A. Example 1 Consider the system (1)  The state and command constraints are The disturbance set is W = {w : −0.1 ≤ w ≤ 0.1}.
The command and state constraints are  The size of the feasible set for our approach is logdet(Y ) = 21.3821. As comparison, for the standard CG, the size of of the feasible set is logdet(P a ) = 3.1909. Note that procedures in [15] to calculate an invariant and constraint-admissible polyhedral set could not be terminated after one hour.
For the initial condition x(0) = [0 0 0 0] T and for the reference r = 15.8540, Fig. 3 shows, respectively, the state and command trajectories of the closed-loop system as functions of time using our approach (solid blue), and the standard CG (dashed red).
Using the TIC/TOC function of MATLAB 2018b, we found that the online computation times is converted into a convex LMI problem. Then, the robust case with bounded disturbances is addressed. It is shown that our fictitious dynamics scheme provides a generalization of the one in [5]. The optimization of the fictitious dynamics requires a solution of an BMI problem, which can be solved by performing a univariate search on a scalar parameter.
The method offers a good control performance over larger feasible sets with small online computational burden as compared to the polyhedral set based CG scheme. Simulation results verifying the performance of the proposed approach have been presented.
We believe that the proposed approaches offer a very attractive alternative to polyhedral set based CG.