Projection onto the Set of Rank-Constrained Structured Matrices for Reduced-Order Controller Design

Abstract

In this paper, we propose an efficient numerical computation method of reduced-order controller design for linear time-invariant systems. The design problem is described by linear matrix inequalities (LMIs) with a rank constraint on a structured matrix, due to which the problem is non-convex. Instead of the heuristic method that approximates the matrix rank by the nuclear norm, we propose a numerical projection onto the rank-constrained set based on the alternating direction method of multipliers (ADMM). Then the controller is obtained by alternating projection between the rank-constrained set and the LMI set. We show the effectiveness of the proposed method compared with existing heuristic methods, by using 95 benchmark models from the COMPL${}_e$ib library.

1. Introduction

It is well known that a stabilizing output-feedback controller and an $H^{\infty }$ controller of a linear time-invariant system can be obtained by solving linear matrix inequalities (LMIs), assuming that the order of the controller is more than or equal to that of the controlled plant model [1,2]. Since the set of optimization variables described by LMIs is convex, the problem can be efficiently solved by convex optimization solvers such as Sedumi [3], SDPT3 [4], and MOSEK [5]. Also, LMIs are easily coded with YALMIP [6] and CVX [7] on MATLAB, and CVXPY [8] on Python.

Practically, it is preferred for implementation to use a low-order controller, especially a static controller of a high-order plant, which we call a reduced-order controller. To obtain a reduced-order controller that has a lower order than the plant is however known to be NP-hard [9] due to a rank constraint [10]. Therefore, we need to employ a heuristic method to efficiently obtain an approximated reduced-order controller. Actually, a couple of heuristic methods have been proposed; the $XY$-centring algorithm [11], the cone complementarity linearization algorithm [12], and alternating projection methods [13,14], to name a few.

More recently, the nuclear norm minimization with LMIs has been proposed to cope with this hard problem [15,16,17,18,19]. This is based on the fact that the nuclear norm of a matrix well approximates the matrix rank [20]. Since the nuclear norm is a convex function and the set described by LMIs is also a convex set, the problem boils down to a convex optimization problem that can be solved very efficiently. The nuclear norm heuristic has been recently applied to, e.g., principal component analysis [21], image denoising [22,23], and system identification [24,25].

Although the nuclear norm heuristic is widely used for rank-constrained problems, we show by numerical examples in this paper that this is not necessarily efficient for reduced-order controller design. Instead, we propose a new method to solve the reduced-order controller design problem by extending the alternating projection method proposed in [13]. The idea is to compute a more precise projection onto the set of rank-constrained structured matrices by the alternating direction method of multipliers (ADMM) [26]. By numerical examples in Section 4, we show that the proposed method significantly improves the precision of the solution compared to the nuclear norm minimization [15] and the original alternating projection method [13].

The organization of this paper is as follows: In Section 2, we show two reduced-order control problems that are described as rank-constrained LMI problems. In Section 3, we propose the alternating projection algorithm to solve the rank-constrained LMI problem. Numerical examples are shown in Section 4 to illustrate the effectiveness of the proposed method with 95 benchmark models from the COMPL${}_e$ib library [27]. A summary is given in Section 5.

We note that the MATLAB programs to check the numerical examples and the results of stability tests for 95 benchmark models shown in Section 4 are available at the web page of [28].

Notation

Let A be a matrix. The transpose of A is denoted by $A^{\top }$, the trace by $\mathrm{tr}\left(A\right)$, and the rank by $\mathrm{rank}\left(A\right)$. The i-th singular value of A is denoted by ${\sigma }_i\left(A\right)$. In this paper, we use two kinds of matrix norms: one is the Frobenius norm $\parallel A\parallel$ of A is defined by

$$\parallel A\parallel \triangleq \sqrt{\mathrm{tr}\left(A^{\top }A\right)}=\sqrt{\sum _{i=1}^n{\sigma }_i^2\left(A\right)}.$$

(1)

Matrix inequalities $A\succ 0$, $A⪰0$, $A\prec 0$, and $A⪯0$ respectively mean A is positive definite, positive semidefinite, negative definite, and negative semidefinite. For $A\in {\mathbb{R}}^{n\times m}$ with $r=\mathrm{rank}\left(A\right)<n$, $A^{\perp }$ is a matrix that satisfies

$$A^{\perp }\in {\mathbb{R}}^{(n-r)\times n},A^{\perp }A=0,A^{\perp }{A}^{\perp \top }\succ 0.$$

(2)

By ${\mathcal{S}}_n$, we denote the set of $n\times n$ real symmetric matrices.

For a closed subset $\Omega$ of ${\mathcal{S}}_n$, the projection operator of $X\in {\mathcal{S}}_n$ onto $\Omega$ is denoted by ${\Pi }_{\Omega }$, that is,

$${\Pi }_{\Omega }\left(X\right)\in \underset{Z\in \Omega }{\arg \min }\parallel Z-X\parallel ,$$

(3)

and the distance from X to $\Omega$ is defined by

$$\mathrm{dist}(X,\Omega )\triangleq \underset{Z\in \Omega }{min}\parallel Z-X\parallel .$$

(4)

2. Reduced-Order Controller Design Problems

In this section, we show two examples of reduced-order controller design.

2.1. Reduced-Order Stabilizing Controllers

Let us consider the following linear time-invariant system:

$$\dot{x}\left(t\right)=Ax\left(t\right)+Bu\left(t\right),y\left(t\right)=Cx\left(t\right),t\ge 0,$$

(5)

where $x\left(t\right)\in {\mathbb{R}}^n$, $u\left(t\right)\in {\mathbb{R}}^m$, $y\left(t\right)\in {\mathbb{R}}^p$, $A\in {\mathbb{R}}^{n\times n}$, $B\in {\mathbb{R}}^{n\times m}$, and $C\in {\mathbb{R}}^{p\times n}$. For this system, we consider an output-feedback controller $u=Ky$, whose order is assumed to be $n_c<n$. Then, the reduced-order output-feedback controller design is described as the following feasibility problem [10].

Problem 1

(Stabilizing controller). Find $X_1,X_2\in {\mathcal{S}}_n$ such that the rank constraint

$$\mathrm{rank}\left[\begin{array}{cc}{X}_1& I\\ I& X_2\end{array}\right]\le n+n_c,$$

(6)

and LMIs

$$\begin{array}{c}-\left[\begin{array}{cc}{X}_1& I\\ I& X_2\end{array}\right]⪯0,\end{array}$$

(7)

$$\begin{array}{c}{B}^{\perp }(A{X}_1+X_1{A}^{\top })B^{\perp \top }⪯-ϵI,\end{array}$$

(8)

$$\begin{array}{c}{C}^{\top \perp }(X_{2}A+A^{\top }{X}_2)C^{\top \perp \top }⪯-ϵI\end{array}$$

(9)

hold for some $ϵ>0$

We note that the inequality “$⪯-ϵI$” can be “$\prec 0$,” however for the projection-based algorithm described in Section 3, we introduce small $ϵ>0$ to make the subsets closed.

2.2. Reduced-Order $H^{\infty }$ Controllers

Let us consider the following generalized plant:

$$G:\left[\begin{array}{c}\dot{x}\left(t\right)\\ z\left(t\right)\\ y\left(t\right)\end{array}\right]=\left[\begin{array}{ccc}A& B_1& B_2\\ C_1& D_{11}& D_{12}\\ C_2& D_{21}& 0\end{array}\right]\left[\begin{array}{c}x\left(t\right)\\ w\left(t\right)\\ u\left(t\right)\end{array}\right],t\ge 0,$$

(10)

where we assume that $(A,B_2)$ is stabilizable and $(C_2,A)$ is detectable. Let $T_{zw}$ denote the feedback connection (or the linear fractional transformation) [29] of G and a controller K such that $u=Ky$. Then, the problem is to seek a controller K of order $n_c\le n$ such that the $H^{\infty }$ norm of $T_{zw}$ satisfies $\parallel T_{zw}{\parallel }_{\infty }<\gamma$ with a given $\gamma >0$. This problem is described as LMIs with a rank constraint [10,30].

Problem 2

($H^{\infty }$ controller). Find $X_1,X_2\in {\mathcal{S}}_n$ such that the rank constraint (6), the LMI (7), and the following LMIs:

$$\begin{array}{c}{\left[\begin{array}{c}{B}_2\\ D_{12}\\ 0\end{array}\right]}^{\perp }\left[\begin{array}{ccc}A{X}_1+X_1{A}^{\top }& X_1{C}_1^{\top }& B_1\\ C_1{X}_1& -\gamma I& D_{11}\\ B_1^{\top }& D_{11}^{\top }& -\gamma I\end{array}\right]{\left[\begin{array}{c}{B}_2\\ D_{12}\\ 0\end{array}\right]}^{\perp \top }⪯-ϵI,\end{array}$$

(11)

$$\begin{array}{c}{\left[\begin{array}{c}{C}_2^{\top }\\ D_{21}^{\top }\\ 0\end{array}\right]}^{\perp }\left[\begin{array}{ccc}{X}_{2}A+A^{\top }{X}_2& X_2{B}_1& C_1^{\top }\\ B_1^{\top }{X}_2& -\gamma I& D_{11}^{\top }\\ C_1& D_{11}& -\gamma I\end{array}\right]{\left[\begin{array}{c}{C}_2^{\top }\\ D_{21}^{\top }\\ 0\end{array}\right]}^{\perp \top }⪯-ϵI\end{array}$$

(12)

hold for some $ϵ>0$.

3. Algorithms

In this section, we propose a new algorithm based on the projection onto the set of rank-constrained structured matrices.

3.1. Proposed Algorithm

First, we define function $F(X_1,X_2)$ such that $F(X_1,X_2)⪯0$ is equivalent to the LMIs to be solved. For Problem 1 for example, $F(X_1,X_2)⪯0$ means that the LMIs (7)–(9) hold. Namely, we consider the following problem including Problems 1 and 2.

Problem 3.

Find a pair of matrices $(X_1,X_2)\in {\mathcal{S}}_n^2$ such that

$$\mathrm{rank}\left[\begin{array}{cc}{X}_1& I\\ I& X_2\end{array}\right]\le r,F(X_1,X_2)⪯0.$$

(13)

Then, we propose alternating projection [31] to solve the rank-constrained LMI problems. For this, we define the two closed subsets of ${\mathcal{S}}_n^2$:

$$\begin{array}{cc}\hfill {\Omega }_r& \triangleq \left\{(X_1,X_2)\in {\mathcal{S}}_n^2:\mathrm{rank}\left[\begin{array}{cc}{X}_1& I\\ I& X_2\end{array}\right]\le r\right\},\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill \Lambda & \triangleq \left\{(X_1,X_2)\in {\mathcal{S}}_n^2:F(X_1,X_2)⪯0\right\}.\hfill \end{array}$$

(15)

Then, the problem is that of finding a pair of matrices $(X_1,X_2)$ in ${\Omega }_r\cap \Lambda$ with $r=n+n_c$. For this, we adapt alternating projection between ${\Omega }_r$ and $\Lambda$. The iterative algorithm is given by

$$\begin{array}{cc}\hfill Z\left[k\right]& ={\Pi }_{{\Omega }_r}\left(X\left[k\right]\right),\hfill \\ \hfill X[k+1]& ={\Pi }_{\Lambda }\left(Z\left[k\right]\right),k=0,1,2,\dots ,\hfill \end{array}$$

(16)

where $X\left[0\right]=(X_1\left[0\right],X_2\left[0\right])\in {\mathcal{S}}_n^2$ is a given initial guess of $(X_1,X_2)$. The computation of the projection operators ${\Pi }_{{\Omega }_r}$ and ${\Pi }_{\Lambda }$ are shown in the next following subsections.

We show a stability result for this algorithm as follows:

Lemma 1.

The sequences $\left\{X\right[k\left]\right\}$ and $\left\{Z\right[k\left]\right\}$ generated by (16) satisfy the following inequalities:

$$\begin{array}{cc}\hfill \parallel X[k+1]-Z[k+1]\parallel & \le \parallel X\left[k\right]-Z\left[k\right]\parallel ,\hfill \\ \hfill \mathrm{dist}(X[k+1],{\Omega }_r)& \le \mathrm{dist}(X\left[k\right],{\Omega }_r),\hfill \\ \hfill \mathrm{dist}\left(Z\right[k+1],\Lambda )& \le \mathrm{dist}\left(Z\right[k],\Lambda ),k=0,1,2,\dots \hfill \end{array}$$

(17)

where dist is the distance function defined in (4).

The proof is given in Appendix A.

Remark 1.

We can also adopt the Dykstra algorithm [32] that gives an element in ${\Omega }_r\cap \Lambda$ that is a projection (i.e., one of the nearest points) on ${\Omega }_r\cap \Lambda$ from the initial guess $X\left[0\right]\in {\mathcal{S}}_n^2$. The algorithm is described as follows:

$$\begin{array}{cc}\hfill Z\left[k\right]& ={\Pi }_{{\Omega }_r}(X\left[k\right]+P\left[k\right]),\hfill \\ \hfill P[k+1]& =X\left[k\right]+P\left[k\right]-Z\left[k\right],\hfill \\ \hfill X[k+1]& ={\Pi }_{\Lambda }(Z\left[k\right]+Q\left[k\right]),\hfill \\ \hfill Q[k+1]& =Z\left[k\right]+Q\left[k\right]-X[k+1],k=0,1,2,\dots ,\hfill \end{array}$$

(18)

where we set $P\left[0\right]=Q\left[0\right]=0$.

3.2. Projection onto the Set ${\Omega }_r$ of Rank-Constrained Structured Matrices

Here we consider the projection of $(X_1,X_2)$ onto the set ${\Omega }_r$ of rank-constrained structured matrices in (14). This projection can be written by definition as

$${\Pi }_{{\Omega }_r}(X_1,X_2)\in \underset{(Z_1,Z_2)\in {\Omega }_r}{\arg \min }\parallel Z_1-X_1{\parallel }^2+{\parallel Z_2-X_2\parallel }^2.$$

(19)

We note that, since the set ${\Omega }_r$ is closed but non-convex, multiple solutions may exist for the minimization in (19) may exist.

For this, we propose a precise projection based on alternating direction method of multipliers (ADMM). For the minimization problem (19), we introduce the indicator function ${\mathcal{I}}_r$ defined by

$${\mathcal{I}}_r\left(Z\right)\triangleq \left\{\begin{array}{cc}0,\hfill & \mathrm{if}\mathrm{rank}\left(Z\right)\le r,\hfill \\ +\infty ,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(20)

Then the minimization problem in (19) is equivalently described as

$$\begin{array}{cc}\hfill \underset{(Z_1,Z_2)\in {\mathcal{S}}_n^2,\tilde{Z}\in {\mathcal{S}}_{2n}}{minimize}& \parallel Z_1-X_1{\parallel }^2+{\parallel Z_2-X_2\parallel }^2+{\mathcal{I}}_r\left(\tilde{Z}\right)\hfill \\ \hfill \mathrm{subject}\mathrm{to}& \tilde{Z}=\left[\begin{array}{cc}{Z}_1& I\\ I& Z_2\end{array}\right]\hfill \end{array}$$

(21)

The convergence of the ADMM algorithm of (21) with nonconvex constraint (20) is discussed in [33].

To solve this optimization problem, we first consider the projection ${\Pi }_{{\mathcal{C}}_r}\left(Z\right)$ onto the set of rank-r matrices

$${\mathcal{C}}_r\triangleq \{Z\in {\mathbb{R}}^{2n\times 2n}:\mathrm{rank}\left(Z\right)\le r\}.$$

(22)

The projection ${\Pi }_{{\mathcal{C}}_r}$ is easily computed via the singular value decomposition $Z=U\Sigma V^{\top }$. Define ${\Sigma }_r$ by setting all but r largest (in magnitude) diagonal entries of $\Sigma$ to 0. Then, the projection ${\Pi }_{{\mathcal{C}}_r}\left(Z\right)$ is given by

$${\Pi }_{{\mathcal{C}}_r}\left(Z\right)=U{\Sigma }_r{V}^{\top }.$$

(23)

Now, the optimization problem in (21) can be efficiently solved by adapting the alternating direction method of multipliers (ADMM) algorithm [26]. The iterative algorithm is given by

$$\begin{array}{cc}{Z}_1[i+1]\hfill & ={\left(1+\frac{\rho }{2}\right)}^{-1}\left(X_1+\frac{\rho }{2}{M}_{11}\left[i\right]\right),\hfill \end{array}$$

(24)

$$\begin{array}{cc}{Z}_2[i+1]\hfill & ={\left(1+\frac{\rho }{2}\right)}^{-1}\left(X_2+\frac{\rho }{2}{M}_{22}\left[i\right]\right),\hfill \end{array}$$

(25)

$$\begin{array}{cc}\tilde{Z}[i+1]\hfill & ={\Pi }_{{\mathcal{C}}_r}\left(\left[\begin{array}{cc}{Z}_1[i+1]& I\\ I& X_2[i+1]\end{array}\right]-W\left[i\right]\right),\hfill \end{array}$$

(26)

$$\begin{array}{cc}W[i+1]\hfill & =W\left[i\right]+\tilde{Z}[i+1]-\left[\begin{array}{cc}{Z}_1[i+1]& I\\ I& Z_2[i+1]\end{array}\right],i=0,1,2,\dots ,\hfill \end{array}$$

(27)

where $\rho >0$ is the step size, and $M_{11}\left[i\right],M_{22}\left[i\right]\in {\mathbb{R}}^{n\times n}$ are defined as

$$\left[\begin{array}{cc}{M}_{11}\left[i\right]& M_{12}\left[i\right]\\ M_{21}\left[i\right]& M_{22}\left[i\right]\end{array}\right]\triangleq \tilde{Z}\left[i\right]+W\left[i\right].$$

(28)

We show in Appendix B how to obtain this iteration algorithm for solving (21).

3.3. Projection onto the Set $\Lambda$ Described by LMIs

The projection of $(X_1,X_2)$ onto the set $\Lambda$ in (15) can be described as convex optimization with LMIs [34,35].

$$\begin{array}{cc}\hfill \underset{Z_1,Z_2,W\in {\mathcal{S}}_n}{minimize}& \mathrm{tr}\left(W\right)\hfill \\ \hfill \mathrm{subject}\mathrm{to}& F(Z_1,Z_2)⪯0,\hfill \\ \hfill & \left[\begin{array}{cc}W& {(Z-X)}^{\top }\\ Z-X& I\end{array}\right]⪰ϵI,\hfill \end{array}$$

(29)

where $Z={\left[Z_1{Z}_2\right]}^{\top }$ and $X={\left[X_1{X}_2\right]}^{\top }$.

Remark 2.

Grigoriadis and Skelton [13] have also proposed an alternating projection method. This method is described as

$$\begin{array}{cc}\hfill Z\left[k\right]& ={\Pi }_{\mathcal{R}}\left({\Pi }_{\mathcal{D}}\left(X\left[k\right]\right)\right)\hfill \\ \hfill X[k+1]& ={\Pi }_{\Lambda }\left(Z\left[k\right]\right),k=0,1,2,\dots ,\hfill \end{array}$$

(30)

where

$$\begin{array}{cc}\hfill \mathcal{D}& \triangleq \left\{Z\in {\mathcal{S}}_{2n}:Z=\left[\begin{array}{cc}{X}_1& 0\\ 0& X_2\end{array}\right],X_1,X_2\in {\mathcal{S}}_n\right\},\hfill \\ \hfill \mathcal{R}& \triangleq \left\{Z\in {\mathcal{S}}_{2n}:\mathrm{rank}(Z+J)\le r\right\},\hfill \\ \hfill J& \triangleq \left[\begin{array}{cc}0& I\\ I& 0\end{array}\right].\hfill \end{array}$$

(31)

We note that${\Omega }_r$ in (14) is the intersection of $\mathcal{D}$ and $\mathcal{R}$. The difference between the proposed algorithm (16) and the above algorithm (30) is the projection onto the set ${\Omega }_r$ of rank-constrained structured matrices. Namely, the algorithm in [13] approximates the projection ${\Pi }_{{\Omega }_r}$ by a composite projection ${\Pi }_{\mathcal{R}}{\Pi }_{\mathcal{D}}$.

Finally, we summarize the proposed algorithm to solve Problem 3 in Algorithm 1.

Algorithm 1 Algorithm to solve Problem 3

Require: Initial guess $(X_1\left[0\right],X_2\left[0\right])\in {\mathcal{S}}_n$ Ensure:$(X_1,X_2)\in {\Omega }_r\cap {\Omega }_{\Lambda }$    for $k=0,1,2,\dots ,N-1$do        $X_1\leftarrow X_1\left[k\right]$                                                                             ▹ Projection onto ${\Omega }_r$        $X_2\leftarrow X_2\left[k\right]$        $\tilde{Z}\left[0\right]\leftarrow 0$        $W\left[0\right]\leftarrow 0$        for $i=0,1,2,\dots ,M-1$ do           $Z_1[i+1]\leftarrow {\left(1+\frac{\rho }{2}\right)}^{-1}\left(X_1+\frac{\rho }{2}{M}_{11}\left[i\right]\right)$           $Z_2[i+1]\leftarrow {\left(1+\frac{\rho }{2}\right)}^{-1}\left(X_2+\frac{\rho }{2}{M}_{22}\left[i\right]\right)$           $\tilde{Z}[i+1]\leftarrow {\Pi }_{{\mathcal{C}}_r}\left(\left[\begin{array}{ccc}{Z}_1[i+1]& II& X_2[i+1]\end{array}\right]-W\left[i\right]\right)$           $[i+1]\leftarrow W\left[i\right]+\tilde{Z}[i+1]-\left[\begin{array}{ccc}{Z}_1[i+1]& II& Z_2[i+1]\end{array}\right]$        end for        $X\leftarrow {[Z_1\left[M\right],Z_2\left[M\right]]}^{\top }$                                                            ▹ Projection onto ${\Omega }_{\Lambda }$        $(Z_1,Z_2,W)\leftarrow \underset{Z_1,Z_2,W\in {\mathcal{S}}_n}{\arg \min }\mathrm{tr}\left(W\right)$ subject to $F(Z_1,Z_2)⪯0$, $\left[\begin{array}{cc}W& {(Z-X)}^{\top }\\ Z-X& I\end{array}\right]⪰ϵI$,        $X_1[k+1]\leftarrow Z_1$        $X_2[k+1]\leftarrow Z_2$ end for $X_1\leftarrow X_1\left[N\right]$ $X_2\leftarrow X_2\left[N\right]$

4. Numerical Examples

In this section, we show some control examples to illustrate the effectiveness of the proposed algorithm. We use benchmark models listed in the COMPL${}_e$ib library [27]. MATLAB programs for the numerical computation in this section can be downloaded from [28]. For numerical optimization in the examples, we use SDPT3 [4] on MATLAB.

4.1. Stabilizing Static Controllers

We consider all benchmark linear time-invariant models whose order is less than 1000 from COMPL${}_e$ib. There are 95 benchmark models to be checked. For these models, we solve Problem 1 by five methods:

1

Nonsmooth $H^{\infty }$ synthesis [36] with hinfstruct funciton in MATLAB (NS$H^{\infty }$)

2

Cone complementarity linearization algorithm [12] (CCL)

3

Nuclear norm minimization [15] (NNM)

4

Alternating projection with approximate projection onto ${\Omega }_r$ [13] (GS96)

5

Alternating projection with the proposed precise projection in Section 3 (Proposed)

Table 1. Stabilizing static controller design results.

	NS${\mathit{H}}^{\mathit{\infty }}$ [36]	CCL [12]	NNM [15]	GS96 [13]	Proposed
# success	75	61	53	56	59
CPU time [s]	0.365	56.5	4.58	663	666

summarizes the number of successful results (i.e., a static stabilizing controller is obtained), and the average computational time.

From this table, we can say that the Nonsmooth $H^{\infty }$ method [36] is the best among the five methods. We note that the alternating methods by [13] and our method show a very long average CPU time since there are a few models (e.g., EB5 and JE1) for which it takes a very long CPU time by the se methods. For example, GS96 takes 53,450 [s] and the proposed method takes 52,315 [s] for EB5. However, there are models for which the proposed method successfully gives a stabilizing static controller while some of the other methods fail. We summarize the results in

Table 2. Stabilizing static controller results (the full list available at [28]).

	NS${\mathit{H}}^{\mathit{\infty }}$ [36]	CCL [12]	NNM [15]	GS96 [13]	Proposed
HE6	unstable	stable	unstable	stable	stable
HE7	stable	stable	unstable	stable	stable
REA3	unstable	stable	unstable	unstable	stable
DIS1	stable	stable	unstable	unstable	stable
PAS	unstable	stable	stable	stable	stable
TF1	unstable	unstable	unstable	unstable	stable
TF2	unstable	stable	stable	stable	stable
NN1	stable	unstable	unstable	unstable	stable
NN11	stable	stable	unstable	stable	stable
NN12	unstable	unstable	unstable	stable	stable
FS	unstable	stable	stable	stable	stable

.

An advantage of the proposed method is found in this table. Since stabilizing static controller design is in general a non-convex problem, and hence there may be no unified approach that gives a solution for any plant models. For example, no methods but the proposed method can compute a stabilizing static controller for TF1. We emphasize our method is an effective method that may provide a solution to a reduced-order controller design problem that cannot be solved by standard methods as nonsmooth $H^{\infty }$, cone complementarity linearization, and nuclear norm minimization.

4.2. Stabilizing Low-Order Controllers

In this section, we focus on the model TF1. In this model, the state-space matrices are given as follows:

$$\begin{array}{cc}\hfill A& =\left[\begin{array}{ccccccc}-1& 0& 0& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& -1& 0& 0\\ -0.088& 0.0345& 0& 0& 1& -0.0032& 0\\ 0& 0& 0.05& 0& 0& 0& -0.00001\end{array}\right],B=\left[\begin{array}{cc}1& 0\\ 0& 0\\ 0& 0\\ 0& 0.09\\ 0& 0\\ 0& 0\\ 0& 0\end{array}\right],\hfill \\ \hfill C& =\left[\begin{array}{ccccccc}0& 0& 1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 0& 1& 0\end{array}\right],D=\left[\begin{array}{cc}0& 0\\ 0& 0\\ 0& 0\\ 0& 0\end{array}\right].\hfill \end{array}$$

(32)

As explained in the previous subsection, a static controller is obtained only with the proposed method. Here, we consider low-order controllers with order $n_c=1$ and $n_c=2$. It is observed that all methods but the nonsmooth $H^{\infty }$ method return stabilizing controllers with order 1 and 2. The results are summarised in

Table 3. Results of stabilizing low-order controllers with $n_c=1$ and 2 for TF1.

	NS${\mathit{H}}^{\mathit{\infty }}$ [36]	CCL [12]	NNM [15]	GS96 [13]	Proposed
TF1 ($n_c=1$)	unstable	stable	stable	stable	stable
CPU time [s]	0.016910	0.214947	0.177307	0.194252	1.836873
TF1 ($n_c=2$)	unstable	stable	stable	stable	stable
CPU time [s]	0.016012	0.263448	0.205307	0.206729	4.250059

.

The state-space matrices $A_K,B_K,C_K,D_K$ of the second order controller ($n_c=2$) with the proposed method is given by

$$\begin{array}{cc}\hfill A_K& =\left[\begin{array}{cc}-1.804\times {10}^4& 1.356\times {10}^4\\ 9512& -7151\end{array}\right],B_K=\left[\begin{array}{cccc}-9193& -2019& -240.4& -1.77\times {10}^4\\ 4848& 1066& 126.8& 9333\end{array}\right],\hfill \\ \hfill C_K& =\left[\begin{array}{cc}-4412& 3317\\ -1.075\times {10}^5& 8.08\times {10}^4\end{array}\right],\hfill \\ \hfill D_K& =\left[\begin{array}{cccc}-2249& -496& -58.91& -4329\\ -5.478\times {10}^4& -1.203\times {10}^4& -1433& -1.055\times {10}^5\end{array}\right].\hfill \end{array}$$

(33)

4.3. $H^{\infty }$ Static Controllers

Then, we consider the reduced-order $H^{\infty }$ control problem formulated in Problem 2. Here we choose AC4, NN1, NN12, and HE6 from COMPL${}_e$ib. AC4 is from a autopilot control problem for an air-to-air missile discussed in [37]. NN1 and NN12 are academic test problems proposed in [38] and [39], respectively. HE6 is a helicopter model that has four inputs, 20 states, and six outputs [40]. For these plant models, we seek the $H^{\infty }$ static controller by using the bisection method on $\gamma$. Namely, we first give a sufficiently large upper bound $\overline{\gamma }$, for example $\overline{\gamma }=100$, and a sufficiently small lower bound $\underline{\gamma }$ (e.g., $\underline{\gamma }=0$). Then we set $\gamma =(\underline{\gamma }+\overline{\gamma })/2=50$ and solve Problem 2. If there is a feasible solution, then we update the upper bound to $\overline{\gamma }=\gamma =50$, otherwise we set the lower bound to $\underline{\gamma }=\gamma =50$. We note that the problem is assumed to be infeasible if a solution of Problem 2 is not found after 15 iterations of (16). We repeat this process until sufficient accuracy is achieved.

Table 4. Upper bounds $\gamma$ of the $H^{\infty }$ norm by static controller.

Model	AC4	NN1	NN12	HE6
$\gamma$	1.000	74.72	28.09	520.0

shows the obtained upper bounds of $\gamma$ for the chosen models.

The obtained static controllers are given as follows:

$$\begin{array}{cc}\hfill K_{\mathrm{AC}4}& =\left[\begin{array}{cc}-0.3228& -0.07534\end{array}\right],\hfill \\ \hfill K_{\mathrm{NN}1}& =\left[\begin{array}{cc}3.5& 60.22\end{array}\right],\hfill \\ \hfill K_{\mathrm{NN}12}& =\left[\begin{array}{cc}22.81& -35.61\\ -5.179& 8.137\end{array}\right],\hfill \\ \hfill K_{\mathrm{HE}6}& =\left[\begin{array}{cccccc}83.25& -0.5581& -0.5931& 0.1238& 0.1546& -0.02533\\ -24.27& 7.876& 0.7263& 0.0309& 0.3272& -0.6042\\ -15.85& 0.1609& -6.578& -2.062& 1.677& 0.1478\\ 65.77& -1.413& 9.07& -16.61& -0.9596& -0.04191\end{array}\right].\hfill \end{array}$$

(34)

It is easy to check that the obtained static controllers really achieve the $H^{\infty }$ norm listed in Table 2. These numerical examples demonstrate the effectiveness of the proposed method.

5. Conclusions

In this paper, we have proposed a novel design method of reduced-order controllers based on projection onto the set of rank-constrained structured matrices. We compared the proposed method with existing methods by numerical examples. We have shown that the proposed algorithm successfully solved benchmark problems that other methods could not. This is thanks to the precise computation of the projection onto the set of rank-constrained structured matrices. Future work includes reduced-order controller design with sparsity constraints on the controller realization, which is a challenging problem that should take two non-convex rank and sparsity constraints.