---
abstract: |
    A novel optimization method is proposed to minimize a convex function subject to bilinear matrix inequality (BMI) constraints.  The key idea is to decompose the bilinear mapping as a difference between two positive semidefinite convex mappings.  At each iteration of the algorithm the concave part is linearized, leading to a convex subproblem. Applications to various output feedback controller synthesis problems are presented.  In these applications the subproblem in each iteration step can be turned into a convex optimization problem with linear matrix inequality (LMI) constraints.  The performance of the algorithm has been benchmarked on the data from COMPl$_\textrm{e}$ib library. 

    0.1cm **Keywords:** [Static feedback controller design, linear time-invariant system, bilinear matrix inequality, semidefinite programming, convex-concave decomposition.]{} 
author:
- 'Tran Dinh Quoc[^1] [ ]{} Suat Gumussoy[^2] [ $\cdot$]{} Wim Michiels Moritz Diehl'
date: 'Technical Report, July 2011'
title: 'Combining Convex-Concave Decompositions and Linearization Approaches for solving BMIs, with application to Static Output Feedback'
---

Introduction
============

Optimization involving matrix constraints have broad interest and applications in static state/output feedback controller design, robust stability of systems, topology optimization (see, e.g. [@BenTal2004; @Boyd1994; @Leibfritz2004; @Kocvara2005]).  Many problems in these fields can be reformulated as an optimization problem of linear matrix inequality (LMI) constraints [@Boyd1994; @Leibfritz2004] which can be solved efficiently and reliably by means of interior point methods for semidefinite programming (SDP) [@BenTal2004; @Nesterov1994] and efficient open-source software tools such as Sedumi [@Sturm1999], SDPT3 [@Tutunku2001].  However, solving optimization problems involving nonlinear matrix inequality constraints is still a big challenge in practice.  The methods and algorithms for nonlinear matrix constrained optimization problems are still limited [@Correa2004; @Freund2007; @Kocvara2005]. 

In control theory, many problems related to the design of a reduced-order controller can be conveniently reformulated as a feasibility problem or an optimization problem with bilinear matrix inequality (BMI) constraints by means of, for instance, Lyapunov’s theory.  The BMI constraints make the problems much more difficult than the LMI ones due to their nonconvexity and possible nonsmoothness.  It has been shown in [@Blondel1997] that the optimization problems involving BMI are NP-hard.  Several approaches to solve optimization problems with BMI constraints have been proposed.  For instance, Goh *et al* [@Goh1995] considered problems in robust control by means of BMI optimization using global optimization methods.  Hol *et al* in [@Hol2005] proposed to used a sum-of-squares approach to fixed order $H$-infinity synthesis.  Apkarian and Tuan [@Apkarian1998] proposed local and global methods for solving BMIs also based on techniques of global optimization.  These authors further considered this problem by proposing parametric formulations and difference of two convex functions (DC) programming approaches.  A similar approach can be found in [@Alamo2008].  However, finding a global optimum is in general impractical while global optimization methods are usually recommended to a low dimensional problem.  Our method developed in this paper is classified as a local optimization method which aims as finding a local optimum based on solving a sequence of convex semidefinite programming problems.  Sequential semidefinite programming method for nonlinear SDP and its application to robust control was considered by Fares *et al* in [@Fares2002].  Thevenet *et al* [@Thevenet2006] studied spectral SDP methods for solving problems involving BMI arising in controller design.  Another approach is based on the fact that problems with BMI constraints can be reformulated as problems with LMI constraints with additional rank constraints.  In [@Orsi2006] Orsi *et al* developed a Newton-like method for solving problems of this type. 

In this paper, we are interested in optimization problems arising in static output feedback controller design for a linear, time-invariant system of the form: $$\label{eq:LTI}
\left\{\begin{array}{cl}
&\dot{x} = Ax + B_1w + Bu,\\
&z = C_1x + D_{11}w + D_{12}u,\\
&y = Cx + D_{21}w,
\end{array}\right.$$ where $x\in\mathbb{R}^{n}$ is state vector, $w\in\mathbb{R}^{n_w}$ is the performance input, $u\in\mathbb{R}^{n_u}$ is input vector, $z\in\mathbb{R}^{n_z}$ is the performance output, $y\in\mathbb{R}^{n_y}$ is physical output vector, $A\in\mathbb{R}^{n\times n}$ is state matrix, $B\in\mathbb{R}^{n\times n_u}$ is input matrix and $C\in\mathbb{R}^{n_y\times n}$ is the output matrix.  Using a static feedback controller of the form $u = Fy$ with $F\in\mathbb{R}^{n_u\times n_y}$, we can write the closed-loop system as follows: $$\label{eq:ss_LTI}
\left\{\begin{array}{cl}
\dot{x}_F = A_Fx_F + B_Fw,\\
z = C_Fx_F + D_Fw. 
\end{array}\right.$$ The stabilization, ${\mathcal{H}}_2$, ${\mathcal{H}}_\infty$ optimization and other control problems for this closed-loop system will be considered. 

0.1cm **Contribution.**  Many control problems can be expressed as optimization problems of BMI constraints and these optimization problems can conveniently be reformulated as optimization problems of difference of two *positive semidefinite convex* (psd-convex) mappings (or convex-concave decomposition) constraints (see Definition \[de:psd\_convex\] below).  In this paper, we propose to use this reformulation leading to a new local optimization method for solving some classes of optimization problems involving BMI constraints.  We provide a practical algorithm and prove the convergence of the algorithm under certain standard assumptions. 

The algorithm proposed in this paper is very simple to implement by using available SDP software tools.  Moreover, it does not require any globalization strategy such as line-search procedures to guarantee global convergence to a local minimum.  The method still works in practice for nonsmooth optimization problems, where the objective function and the concave parts are only subdifferentiable, but not necessarily differentiable.  Note that our method is different from the standard DCA approach in [@Pham1998] since we work directly with positive semidefinite matrix inequality constraints instead of transforming into DC representations as in [@Apkarian1998]. 

We show that our method is applicable to many control problems in static state/output feedback controller design.  The numerical results are benchmarked using the data from COMPl$_\textrm{e}$ib library.  Note, however, that this method is also applicable to other nonconvex optimization problems with matrix inequality constraints which can be written as a convex-concave decomposition. 

0.1cm **Outline of the paper.**  The remainder of the paper is organized as follows.  Section \[sec:premi\] provides some preliminary results which will be used in what follows.  Section \[sec:MDCP\] presents the formulation of optimization problems involving convex-concave matrix inequality constraints and a fundamental assumption, Assumption A\[as:A1\].  The algorithm and its convergence results are presented in Section \[sec:SSDP\].  Applications to control problems on static feedback controller design and numerical benchmarking are given in Section \[sec:applications\].  The last section contains some concluding remarks. 

Preliminaries {#sec:premi}
=============

Let $\mathcal{S}^p$ be the set of symmetric matrices of size $p\times p$, $\mathcal{S}^p_{+}$, and resp.,  $\mathcal{S}^p_{++}$ be the set of symmetric positive semidefinite, resp.,  symmetric positive definite matrices.  For given matrices $X$ and $Y$ in $\mathcal{S}^p$, the relation $X\succeq Y$ (resp.,  $X\preceq Y$) means that $X-Y\in\mathcal{S}^p_{+}$ (resp.,  $Y-X\in\mathcal{S}^p_{+}$) and $X\succ Y$ (resp.,  $X\prec
Y$) is $X-Y\in\mathcal{S}^p_{++}$ (resp.,  $Y-X\in\mathcal{S}^p_{++}$).  The quantity $X\circ Y := \textrm{trace}(X^TY)$ is an inner product of two matrices $X$ and $Y$ defined on $\mathcal{S}^p$, where $\textrm{trace}(Z)$ is the trace of matrix $Z$.

\[de:psd\_convex\][@Shapiro1997] A matrix-valued mapping $G : \mathbb{R}^n\to\mathcal{S}^p$ is said to be positive semidefinite convex (*psd-convex*) on a convex subset $C\subseteq\mathbb{R}^n$ if for all $t\in[0,1]$ and $x,y\in C$, one has $$\label{eq:psd_convex}
G(tx+(1-t)y) \preceq tG(x) + (1-t)G(y).$$ If holds true for $\preceq$ instead of $\prec$ and $t\in (0,1)$ then $G$ is said to be *strictly psd-convex* on $C$.  Alternatively, if we replace $\preceq$ in by $\succeq$ then $G$ is said to be psd-concave on $C$.  It is obvious that any convex function $f : \mathbb{R}^n\to\mathbb{R}$ is psd-convex ($p=1$). 

A function $f :\mathbb{R}^n\to\mathbb{R}$ is said to be *strongly convex* with the parameter $\rho > 0$ if $f(\cdot) - \frac{1}{2}\rho{\|\cdot\|}^2$ is convex. 

The derivative of a matrix-valued mapping $G$ at $x$ is a linear mapping $DG$ from $\mathbb{R}^n$ to $\mathbb{R}^{p\times p}$ which is defined by $$DG(x)h := \sum_{i=1}^nh_i\frac{\partial{G}}{\partial x_i}(x), ~\forall h\in \mathbb{R}^n.$$ -0.1cm For a given convex set $X\in\mathbb{R}^n$, the matrix-valued mapping $G$ is said to be differentiable on a subset $X$ if its derivative $DG(x)$ exists at every $x\in X$. The definitions of the second order derivatives of matrix-valued mappings can be found, e.g., in [@Shapiro1997].  Let $A : \mathbb{R}^n\to \mathcal{S}^p$ be a linear mapping defined as $Ax = \sum_{i=1}^nx_iA_i$, where $A_i \in \mathcal{S}^p$ for $i=1,\dots, n$. The adjoint operator of $A$, $A^{*}$, is defined as $A^{*}Z = (A_1\circ Z, A_2\circ Z, \dots, A_n\circ Z)^T$ for any $Z\in\mathcal{S}^p$.

\[le:psd\_convex\_properties\]

-   A matrix-valued mapping $G$ is psd-convex on $X$ if and only if for any $v\in\mathbb{R}^p$ the function $\varphi(x) := v^TG(x)v$ is convex on $X$.

-   A mapping $G$ is psd-convex on $X$ if and only if for all $x$ and $y$ in $X$, one has $$\label{eq:psd_convex_pro1}
    G(y) - G(x) \succeq DG(x)(y-x).$$

The proof of the statement a) can be found in [@Shapiro1997].  We prove b).  Let $\varphi(x) = v^TG(x)v$ for any $v\in\mathbb{R}^p$. If $G$ is psd-convex then $\varphi$ is convex.  We have $\varphi(y) - \varphi(x) \geq \nabla\varphi(x)^T(y-x)$. Now, $\nabla\varphi(x)^T(y-x) =
\sum_{i=1}^n(y_i-x_i)v^T\frac{\partial{G}}{\partial{x}_i}(x)v = v^T[DG(x)(y-x)]v$. Hence, $v^T\left[G(y)-G(x)-DG(x)\right.$ $\left.(y-x)\right]v \geq 0$ for all $v$. We conclude that holds.  Conversely, if holds then, for any $v$, we have $v^T\left[G(y) - G(x) - DG(x)(y-x)\right]v \geq 0$, which is equivalent to $\varphi(y) - \varphi(x) \geq
\nabla\varphi(x)^T(y-x)$. Thus $\varphi$ is convex.  By virtue of a), the mapping $G$ is psd-convex. 

For simplicity of discussion, throughout this paper, we assume that all the functions and matrix-valued mappings are *twice differentiable* on their domain [@Shapiro1997; @Thevenet2006].  However, this assumption can be reduced to the *subdifferentiability* of the objective function and the concave parts of the matrix-valued mappings. 

\[de:DC\_decomp\] A matrix-valued mapping $F :\mathbb{R}^n\to\mathcal{S}^p$ is said to be a *psd-convex-concave* mapping if $F$ can be represented as a difference of two psd-convex mappings, i.e. $F(x) = G(x) - H(x)$, where $G$ and $H$ are psd-convex.  The pair $(H,G)$ is called a psd-DC (or psd-convex-concave) decomposition of $F$.

Note that each given psd-convex-concave mapping possesses many psd-convex-concave decompositions. 

Optimization of convex-concave matrix inequality constraints {#sec:MDCP}
============================================================

Psd-convex-concave decomposition of BMIs
----------------------------------------

Instead of using the vector $x$ as a decision variable, we use from now on the matrix $X$ as a matrix variable in $\mathbb{R}^{m\times n}$. Note that any matrix $X$ can be considered as an $m\times n$-column vector by vectorizing with respect to its columns, i.e. $x = \textrm{vec}(X) := (X_{11},
X_{21}, \dots, X_{mn})^T$. The inverse mapping of $\textrm{vec}$ is called $\textrm{mat}$. Since $\textrm{vec}$ and $\textrm{mat}$ are linear operators, the psd-convexity is still preserved under these operators. 

A mapping $F : \mathbb{R}^{p\times q}\times\mathcal{S}^{p} \to \mathcal{S}^p$ given by $F(X,Y) := XQ^{-1}X^T - Y$, where $Q\in\mathcal{S}^q$ is symmetric positive definite, is called a *Schur psd-convex*[^3] mapping. 

Consider a bilinear matrix form $$\label{eq:bilinearform}
F(X, Y) := X^TY + Y^TX.$$ By using the Kronecker product, we can write $F$ as $\textrm{vec}(F(X,Y)) = (I_n\otimes X^T)\textrm{vec}(Y) + (I_y\otimes Y^T)\textrm{vec}(X) =
(\sum_{i,j}x_iy_j)$, where $I_n$, $I_y$ are appropriate identity matrices, $\otimes$ denotes the Kronecker product.  Hence, the vectorization of $F(X,Y)$ is indeed a bilinear form of two vectors $x := \textrm{vec}(X)$ and $y:=\textrm{vec}(Y)$.

The following lemma shows that the bilinear matrix form can be decomposed as a difference of two psd-convex mappings. 

\[le:key\_lemma\_a\]

-   The mapping $f(X) := X^TX$, $g(X) := XX^T$) are psd-convex on $\mathbb{R}^{m\times n}$. The mapping $f(X) := X^{-1}$ is psd-convex on $\mathcal{S}^{p}_{++}$.

-   The bilinear matrix form $X^TY + Y^TX$ can be represented as a psd-convex-concave mapping of at least three forms: $$\begin{aligned}
    \label{eq:DC_form}
    X^TY + Y^TX & =  (X+Y)^T(X+Y) - (X^TX + Y^TY)\nonumber\\
    & = X^TX + Y^TY - (X-Y)^T(X-Y)\\ 
    & = \frac{1}{2}[(X+Y)^T(X+Y) - (X-Y)^T(X-Y)].\nonumber\end{aligned}$$

The statement b) provides at least three different explicit psd-convex-concave decompositions of the bilinear form $X^TY+Y^TX$. Intuitively, we can see that the first decomposition has a “strong curvature” on the second term, while the second and the third decompositions have “less curvature” on the second term due in case of a compensation between $X$ and $Y$.

The following result will be used to transform Shur psd-convex constraints to LMI constraints. 

\[le:key\_lemma\]

-   Suppose that $A\in\mathcal{S}^n$. Then the matrix inequality $BB^T - A \prec (\preceq)~ 0$ is equivalent to $$\label{eq:Schur_comp}
    \begin{bmatrix}A & B\\ B^T & I\end{bmatrix} \succ (\succeq)~ 0.$$

-   Suppose that $A\in\mathcal{S}^n$, $D \succ 0$, then we have: $$\begin{bmatrix}
    A - BB^T & C\\ C^T & D
    \end{bmatrix} \succ (\succeq) ~ 0 ~~
    \Longleftrightarrow~~
    \begin{bmatrix}
    A & B & C\\ B^T & I & O\\ C^T & O & D
    \end{bmatrix} \succ (\succeq)~0.$$

The proof of this lemma immediately follows by applying Schur’s complement and Lemma \[le:psd\_convex\_properties\] [@Boyd2004].  We omit the proof here. 

Optimization involving convex-concave matrix inequality constraints
-------------------------------------------------------------------

Let us consider the following optimization problem: $$\label{eq:MDC_program}
\left\{\begin{array}{cl}
\displaystyle\min_{x}&f(x)\\
\textrm{s.t.}&  G_i(x) - H_i(x) \preceq 0, ~i=1,\dots, l,\\
&x\in\Omega,
\end{array}\right.$$ where $f :\mathbb{R}^n\to\mathbb{R}$ is convex, $\Omega\subseteq\mathbb{R}^n$ is a nonempty, closed convex set, and $G_i$ and $H_i$ ($i=1,\dots,l$) are psd-convex.  Problem is referred to as a convex optimization with psd-convex-concave matrix inequality constraints. 

Let $\Omega$ be a polyhedral in $\mathbb{R}^n$. Then, if $f$ is nonlinear or one of the mappings $G_i$ or $H_i$ ($i=1,\dots,l)$ is nonlinear then is a nonlinear semidefinite program.  If $H_i$ ($i=1,\dots, l$) are linear then is a convex nonlinear SDP problem.  Otherwise, it is a nonconvex nonlinear SDP problem. 

Let us define $L(x,\Lambda) := f(x) + \sum_{i=1}^l\Lambda_i\circ [G_i(x)-H_i(x)]$ as the Lagrange function of , where $\Lambda_i\in\mathcal{S}^p$ ($i=1,\dots, l$) considered as Lagrange multipliers.  The generalized KKT condition of is presented as: $$\label{eq:DC_KKT}
\left\{\begin{array}{cl}
&0 \in \nabla f(x) + \sum_{i=1}^lD[G_i(x)-H_i(x)]^{*}\Lambda_i + N_{\Omega}(x),\\
&G_i(x)-H_i(x)\preceq 0,~ \Lambda_i\succeq 0,\\
&[G_i(x)-H_i(x)]\circ\Lambda_i=0, ~i=1,\dots, l.
\end{array}\right.$$ Here, $N_{\Omega}(x)$ is the normal cone of $\Omega$ at $x$ defined as $$N_{\Omega}(x) := \left\{\begin{array}{cl}
\left\{w\in\mathbb{R}^n ~|~ w^T(y-x)\geq 0, ~\forall y\in \Omega\right\}, &~\textrm{if}~ x\in\Omega,\\
\emptyset, &~\textrm{otherwise}. 
\end{array}\right.$$ A pair $(x^{*}, \Lambda^{*})$ satisfying is called a KKT point, $x^{*}$ is called a stationary point and $\Lambda^{*}$ is the corresponding multiplier of .  The generalized optimality condition for nonlinear semidefinite programming can be found in the literature (e.g., [@Shapiro1997; @Sun2006]). 

Let us denote by $$\label{eq:feasible_set}
\mathcal{D} := \left\{x\in\Omega ~|~ G_i(x)-H_i(x)\preceq 0, ~l=1,\dots,l\right\},$$ the feasible set of and $\textrm{ri}(\mathcal{D})$ is the relative interior of $\mathcal{D}$ which is defined by $$\textrm{ri}(\mathcal{D}) := \left\{
x\in\textrm{ri}(\Omega) ~|~ G_i(x) - H_i(x) \prec 0, ~i=1,\dots, m\right\},$$ where $\textrm{ri}(\Omega)$ is the set of classical relative interiors of $\Omega$ [@Boyd2004; @Horst1996].  The following condition is a fundamental assumption in this paper. 

\[as:A1\] $\textrm{ri}(\mathcal{D})$ is nonempty. 

Note that this assumption is crucial for our method, because, as we shall see, it requires a strictly feasible starting point $x^0\in\textrm{ri}(\mathcal{D})$. Finding such a point is in principle not an easy task.  However, in many problems, this assumption is always satisfied.  In Section \[sec:applications\] we will propose techniques to determine a starting point for the control problems under consideration. 

The algorithm and its convergence {#sec:SSDP}
=================================

In this section, a local optimization method for finding a stationary point of problem is proposed.  Motivated from the DC programming algorithm developed in [@Pham1998] and the convex-concave procedure in [@Sriperumbudur2009] for scalar functions, we develop an iterative procedure for finding a stationary point of .  The main idea is to linearize the nonconvex part of the psd-convex-concave matrix inequality constraints and then transform the linearized subproblem into a quadratic semidefinite programming problem.  The subproblem can be either directly solved by means of interior point methods or transformed into a quadratic problem with LMI constraints.  In the latter case, the resulting problem can be solved by available software tools such as Sedumi [@Sturm1999] and SDPT3 [@Tutunku2001]. 

The algorithm
-------------

Suppose that $x^k\in\Omega$ is a given point, the linearized problem of around $x^k$ is written as $$\label{eq:psd_convex_subprob}
\left\{ \begin{array}{cl}
\displaystyle\min_{x}  & \left\{ f_k(x) := f(x) + \frac{\rho_k}{2}\|Q_k(x-x_k)\|^2_2\right\}\\
\textrm{s.t.}   &\!\!\!  G_i(x) \!-\!  H_i(x^k) \!-\!  DH_i(x^k)(x \!-\!  x^k) \!  \preceq \!0, \!~i=1,\dots, l,\\
&\!\!\!  x\in\Omega. 
\end{array}\right.$$ Here, we add a regularization term into the objective function of the original problem, where $Q_k$ is a given matrix that projects $x-x_k$ in a certain subspace of $\mathbb{R}^n$ and $\rho_k \geq 0$ is a regularization parameter.  Since $G_i$ ($i=1,\dots, l$) are psd-convex and the objective function is convex, problem is convex.  The *linearized convex-concave SDP algorithm* for solving is described as follows. 

\[alg:A1\][ ]{}\
**Initialization:** Choose a positive number $\rho_0$ and a matrix $Q_0\in\mathcal{S}^n_{+}$. Find an initial point $x^0\in\textrm{ri}(\mathcal{D})$. Set $k := 0$.\
**Iteration $k$:** For $k=0, 1,\dots$. Perform the following steps:

-   *Step 1:* Solve the convex semidefinite program to obtain a solution $x^{k+1}$.

-   *Step 2:* If ${\|x^{k+1} - x^{k}\|}\leq \varepsilon$ for a given tolerance $\varepsilon>0$ then terminate.  Otherwise, update $\rho_k$ and $Q_k$ (if necessary), set $k:=k+1$ and go back to *Step 1*. 

The following main property of the method makes an implementation very easy.  If the initial point $x^0$ belongs to the relative interior of the feasible set $\mathcal{D}$, i.e. $x^0\in\textrm{ri}(\mathcal{D})$, then Algorithm \[alg:A1\] generates a sequence ${x^k}$ which still belongs to $\mathcal{D}$. Consequently, no line-search procedure is needed to ensure the global convergence. 

This property follows from the fact that the linearization of the concave part $-H_i$ is its an upper approximation of this mapping (in the sense of positive semidefinite cone), i.e. $$-H_i(x) \preceq -H_i(x^k) - DH_i(x^k)(x-x^k), \forall x\in\Omega,$$ which is equivalent to $$G_i(x) - H_i(x) \preceq G_i(x) - H_i(x^k) - DH_i(x^k)(x-x^k), \forall x\in\Omega.$$ Hence, if the subproblem has a solution $x^{k+1}$ then it is feasible to .  Geometrically, Algorithm \[alg:A1\] can be seen as an inner approximate method. 

The main tasks of an implementation of Algorithm \[alg:A1\] consist of:

-   determining an initial point $x^0\in\textrm{ri}(\mathcal{D})$, and

-   solving the convex semidefinite program repeatedly. 

As mentioned before, since $\mathcal{D}$ is nonconvex, finding an initial point $x^0$ in $\textrm{ri}(\mathcal{D})$ is, in principle, not an easy task.  However, in some practical problems, this can be done by exploiting the special structure of the problem (see the examples in Section \[sec:applications\]). 

To solve the convex subproblem , we can either implement an interior point method and exploit the structure of the problem or transform it into a standard SDP problem and then make use of available software tools for SDP.  The regularization parameter $\rho_k$ and the projection matrix $Q_k$ can be fixed at appropriate choices for all iterations, or adaptively updated. 

\[le:KKT\_point\] If $x^{k}$ is a solution of linearized at $x^k$ then it is a stationary point of . 

Suppose that $\Lambda^{k+1}$ is a multiplier associated with $x^k$, substituting $x^k$ into the generalized KKT condition of we obtain .  Thus $x^k$ is a stationary point of . 

Convergence analysis
--------------------

-0.3cm In this subsection, we restrict our discussion to the following special case. 

\[as:A2\] The mappings $G_i$ ($i=1,\dots, l$) are Schur psd-convex and $\Omega$ is formed by a finite number of LMIs.  In addition, $f$ is convex quadratic on $\mathbb{R}^n$ with a convexity parameter $\rho_f\geq 0$.

This assumption is only technical for our implementation.  If the mapping $G_i$ is Schur psd-convex then the linearized constraints of problem can directly be transformed into LMI constraints (see Lemma \[le:key\_lemma\]).  In practice, $G_i$ ($i=1,\dots,l$) can be a general psd-convex mappings and $f$ can be a general convex function. 

Under Assumption A\[as:A2\], the convex subproblem can be transformed equivalently into a quadratic semidefinite program of the form: $$\label{eq:LSDP}
\left\{\begin{array}{cl}
\displaystyle\min_{z\in\mathbb{R}^{n_z}} & \frac{1}{2}z^TBz + h^Tz\\
\textrm{s.t.} & A(z) + C \preceq 0,
\end{array}\right.$$ where $A$ is a linear mapping from $\mathbb{R}^{n_z}$ to $\mathcal{S}^{p_z}$, $C\in\mathcal{S}^{p_z}$ and $B$ is a symmetric matrix, by means of Lemma \[le:key\_lemma\]. 

A vector $\hat{z}$ is said to satisfy the *Slater condition* of if $A(\hat{z}) + C \prec 0$. Suppose that the triple $(\bar{z},\bar{V},\bar{S})$ satisfies the KKT condition of (see [@Freund2007]), where $\bar{z}$ is a primal stationary point, $\bar{V}$ is a Lagrange multiplier and $\bar{S}$ is a slack variable associated with $\bar{z}$ and $\bar{V}$. Then, problem is said to satisfy the *strict complementarity* condition at $(\bar{z},\bar{V},\bar{S})$ if $\bar{V} + \bar{S} \succ 0$.

Let $\bar{z}$ be a stationary point of .  We say that $0\neq p\in\mathbb{R}^{n_z}$ is a *feasible direction* to if $\bar{z}
+ \varepsilon p$ is a feasible point of for all $\varepsilon>0$ sufficiently small.  As in [@Freund2007], we assume that the *second order sufficient condition* holds for at $\bar{z}$ with modulus $\mu > 0$ if for all feasible directions $p$ at $\bar{z}$ with $p^T(h+B\bar{z}) = 0$, one has $p^TBp \geq \mu{\|p\|}^2$. We say that the convex problem is solvable and satisfies the *strong second order sufficient condition* if there exists a KKT point $(\bar{z},\bar{V},\bar{S})$ of the KKT system of satisfies the second order sufficient condition and the strict complementary condition. 

\[as:A3\] The convex subproblem is *solvable* and satisfies the *strong second order sufficient condition*. 

Assumption A\[as:A3\] is standard in optimization and is usually used to investigate the convergence of the algorithms [@Fares2002; @Freund2007; @Shapiro1997]. 

The following lemma shows that $\Delta x^k := x^{k+1}-x^k$ is a descent direction of problem whose proof is given in the Appendix.  \[app:A1\]. 

\[le:descent\_dir\] Suppose that $\{(x^k,\Lambda^k)\}_{k\geq 0}$ is a sequence generated by Algorithm \[alg:A1\].  Then:

-   $\mathrm{a)}$ The following inequality holds for $k\geq 0$: $$\begin{aligned}
    \label{eq:descent_dir}
    f(x^{k + 1}) - f(x^k) \leq -\frac{\rho_f}{2}{\|x^{k + 1} - x^k\|}_2^2 - \rho_k{\|Q_k(x^{k+1} - x^k)\|}^2_2,\end{aligned}$$ where $\rho_f$ is the convexity parameter of $f$.

-   $\mathrm{b)}$ If there exists at least one constraint $i_0$, $i_0\in\{1,2,\dots,l\}$, to be strictly feasible at $x^k$, i.e. $G_{i_0}(x^k) - H_{i_0}(x^k)
    \prec 0$, then $f(x^{k+1}) <
    f(x^k)$ provided that $\Lambda^{k+1}_{i_0} \succ 0$.

-   $\mathrm{c)}$ If $\rho_k > 0$ and $Q_k$ is full-row-rank then $\Delta x^k$ is a sufficiently descent direction of . 

The following theorem shows the convergence of Algorithm \[alg:A1\] in a particular case. 

\[th:convergence\] Under Assumptions A\[as:A1\], A\[as:A2\] and A\[as:A3\], suppose that $f$ is bounded from below on $\mathcal{D}$, where $\mathcal{D}$ is assumed to be bounded in $\mathbb{R}^n$. Let $\{(x^k,\Lambda^k)\}$ be a sequence generated by Algorithm \[alg:A1\] starting from $x^0\in\textrm{ri}(\mathcal{D})$. Then if either $f$ is strongly convex or $\rho_k \equiv \rho > 0$ and $Q_k \equiv Q$ is full-row-rank for all $k\geq 0$ then every accumulation point $(x^{*},\Lambda^{*})$ of $\{(x^k, \Lambda^k)\}$ is a KKT point of .  Moreover, if the set of the KKT points of is finite then the whole sequence $\{ (x^k, \Lambda^k) \}$ converges to a KKT point of . 

Let $M(x^0) := \left\{x^k\right\}$ be a sequence of the sample points generated by Algorithm \[alg:A1\] starting from $x^0$. For a given $x\in\Omega$, let us define the following mapping: -0.3cm $$\begin{aligned}
A_{\textrm{sol}}(x) := \text{arg}\!\min \Big\{\!\!\!\!\! &&  f(y) + \frac{\rho}{2}{\|Q(y - x)\|}_2^2 ~~|~~ y\in\Omega, \\
&& G_i(y) -  H_i(x) - DH_i(x)(y - x) \preceq 0, ~i=1,\dots,m
\Big\}.\notag\end{aligned}$$ Then, $A_{\textrm{sol}}$ is a multivalued mapping and it can be considered as the solution mapping of the convex subproblem .  Note that the sequence $\{x^k\}$ generated by Algorithm \[alg:A1\] satisfies $x^{k+1} \in A_{\textrm{sol}}(x^k)$.\
We first prove that $A_{\textrm{sol}}$ is a closed mapping.  Indeed, since the convex subproblem satisfies Slater’s condition and has a solution that satisfies the strict complementarity and the second order sufficient condition, applying Theorem 1 in [@Freund2007] we conclude that the mapping $A_{\textrm{sol}}$ is differentiable in a neighborhood of the solution.  Consequently, it is closed due to the compactness of $\mathcal{D}$.

On the other hand, since $f$ is either strongly convex or $\rho_k \equiv \rho >0$ for all $k\geq 0$ and $Q_k\equiv Q$ is full-row-rank, it follows from Lemma \[le:descent\_dir\] that the objective function $f$ is strictly monotone on $M(x^0)$. Since $M(x^0)\subseteq\mathcal{D}$ and $\mathcal{D}$ is compact, $M(x^0)$ is also compact.  Applying Theorem 2 in [@Meyer1976] we conclude that every limit points of the sequence $\{x^k\}$ belongs to the set of stationary points $S^{*}$. Moreover, $S^{*}$ is connected and if $S^{*}$ is finite then the whole sequence $\{x^k\}$ converges to $x^{*}$ in $S^{*}$.

\[re:quadratic\] The condition that $f$ is quadratic in Assumption \[as:A2\] can be relaxed to $f$ being twice continuously differentiable.  However, in this case, we need a direct proof for Theorem \[th:convergence\] instead of applying Theorem 1 in [@Freund2007]. 

Applications to robust controller design {#sec:applications}
========================================

In this section, we apply the method developed in the previous section to the following static state/output feedback controller design problems:

1.   Sparse linear static output feedback controller design;

2.   Spectral abscissa and pseudo-spectral abscissa optimization;

3.   $\mathcal{H}_2$ optimization;

4.   $\mathcal{H}_{\infty}$ optimization;

5.   and mixed $\mathcal{H}_2/\mathcal{H}_{\infty}$ synthesis. 

We used the system data from [@Hassibi1999; @Ostertag2008] and the COMP$\textrm{l}_\textrm{e}$ib library [@Leibfritz2003].  All the implementations are done in Matlab 7.11.0 (R2010b) running on a PC Desktop Intel(R) Core(TM)2 Quad CPU Q6600 with 2.4GHz and 3Gb RAM.  We use the YALMIP package [@Lofberg2004] with the SeDuMi 1.1 solver [@Sturm1999] to solve the LMI optimization problems arising in Algorithm \[alg:A1\] at the initial phase (Phase 1) and subproblem .  The Matlab codes can be downloaded at <http://www.kuleuven.be/optec/software/BMIsolver>.  We also benchmarked our method with various examples and compared our results with HIFOO [@Gumussoy09] and PENBMI [@Henrion2005] for all control problems.  HIFOO is an open-source Matlab package for fixed-order controller design.  It computes a fixed-order controller using a hybrid algorithm for nonsmooth, nonconvex optimization based on quasi-Newton updating and gradient sampling.  PENBMI [@Henrion2005] is a commercial software for solving optimization problems with quadratic objective and BMI constraints.  PENBMI is free licensed for academic purposes.  We initialized the initial controller for HIFOO and the BMI parameters for PENBMI to the initial values of our method.  As we shall see, we can reformulate the spectral abscissa optimization problem as a rank constrained LMI problem.  Therefore, we also compared our results with LMIRank [@Orsi2006], a MATLAB toolbox for solving rank constrained LMI problems, for the spectral abscissa optimization. 

Note that all problems addressed here lead to at least one BMI constraint.  To apply the method developed in the previous section, we propose a unified scheme to treat these problems. 

\[scheme:A1\][ ]{}

-   *Step 1.*  Find a convex-concave decomposition of the BMI constraints as $G(x)-H(x) \preceq 0$.

-   *Step 2.*  Find a starting point $x^0\in\textrm{ri}(\mathcal{D})$.

-   *Step 3.*  For a given $x^k$, linearize the concave part to obtain the convex constraint $G(x) - H_k(x) \preceq 0$, where $H_k$ is the linearization of $H$ at $x^k$.

-   *Step 4.*  Reformulate the convex constraint as LMI constraint by means of Lemma \[le:key\_lemma\]. 

-   *Step 5.*  Apply Algorithm \[alg:A1\] with SDP solver to solve the problem. 

Sparse linear constant output-feedback design {#subsec:sparse_SOF}
---------------------------------------------

Let us consider a BMI optimization problem of sparse linear constant output-feedback design given as: $$\label{eq:BMI_p1}
\begin{array}{cl}
\displaystyle\min_{\alpha, P, F} & -\sigma\alpha+\sum_{i=1}^{n_u}\sum_{j=1}^{n_y}|F_{ij}|\\
\text{s.t} ~& (A+BFC)^TP + P(A+BFC) + 2\alpha P \prec 0, ~~P = P^T,~ P\succ 0. 
\end{array}$$ Here, matrices $A$, $B$, $C$ are given with appropriate dimensions, $P$ and $F$ are referred as variables and $\sigma > 0$ is a weighting parameter.  The objective function consists of two terms: the first term $\sigma\alpha$ is to stabilize the system (or to maximize the decay rate) and the second one is to ensure the sparsity of the gain matrix $F$. This problem is a modification of the first example in [@Hassibi1999].  Let us illustrate <span style="font-variant:small-caps;">Scheme A.\[scheme:A1\]</span> for solving this problem. 

*Step 1:* Let $B_F := A + BFC + \alpha I$, where $I$ is the identity matrix.  Then, applying Lemma \[le:key\_lemma\_a\] we can write $$\begin{aligned}
&(A \!+\!  BFC)^T\!\!P \!+\!  P(A\!+\!BFC) + 2\alpha P = B_F^TP + PB_F\nonumber\\
& = B_F^TB_F + P^TP - (B_F-P)^T(B_F-P),\label{eq:BMI_p1_dc1}\\
& = \frac{1}{2}\left[(B_F+P)^T(B_F+P) - (B_F-P)^T(B_F-P)\right].\label{eq:BMI_p1_dc2}\end{aligned}$$ In our implementation, we use the decomposition . 

If we denote by $$\begin{aligned}
\label{eq:G_H_mapping}
G(\alpha,P,F) :=\frac{1}{2}(B_F+P)^T(B_F+P), ~~ \textrm{and} ~~ H(\alpha,P,F) :=\frac{1}{2}(B_F-P)^T(B_F-P),$$ then the BMI constraint in can be written equivalently as a psd-convex-concave matrix inequality constraint (of a variable $x$ formed from $(\alpha, P, F)$ as $x := (\alpha,
\text{vec}(P)^T, \textrm{vec}(F)^T)^T$) as follows: $$\begin{aligned}
\label{eq:DC_const_p1}
G(\alpha,P,F) - H(\alpha,P,F) \prec 0.\end{aligned}$$ Note that the objective function of is convex but nonsmooth which is not directly suitable for the SSDP approach in [@Correa2004], but, the nonconvex problem can be reformulated in the form of using slack variables. 

*Steps 2-5:* The implementation is carried out as follows:

-   **Phase 1.** (* Determine a starting point $x^0\in\textrm{ri}(\mathcal{D})$*).  Set $F^0 := 0$, $\alpha^0 := -\alpha_0(A^T+A)/2$ where $\alpha_0(\cdot)$ is the maximum real part of the eigenvalues of the matrix, and compute $P=P^0$ as the solution of the LMI feasibility problem $$\label{eq:LMI_phase1}
    (A+BF^0 C)^TP + P(A+BF^0 C) + 2\alpha^0 P \prec 0.$$ The above choice for $(\alpha^0,F^0)$ originates from the property that $P^0 = I$ renders the left hand size of (\[eq:LMI\_phase1\]) negative semi-definite (but not negative definite). 

-   **Phase 2.**  Perform Algorithm \[alg:A1\] with a starting point $x^0$ found at Phase 1. 

Let us now illustrate *Step 4* of <span style="font-variant:small-caps;">Scheme A.\[scheme:A1\]</span>.  After linearizing the concave part of the convex-concave reformulation of the last BMI constraint in at $(F^k, P^k, \alpha^k)$ we obtain the linearization: $$\label{eq:Shur_form2}
(A \!+\!  BFC \!+\!  \alpha I \!+\!  P)^T(A \!+\!  BFC \!+\!  \alpha I \!+\!  P) - H_k(F,P,\alpha) \prec 0,$$ where $H_k(F,P,\alpha)$ is a linear mapping of $F$, $P$ and $\alpha$. Now, applying Lemma \[le:key\_lemma\], can be transformed to an LMI constraint: $$\begin{bmatrix}
H_k(F, P, \alpha) & (A \!+\!  BFC \!+\!  \alpha I \!+\!  P)^T\\
(A \!+\!  BFC \!+\!  \alpha I \!+\!  P) & I
\end{bmatrix}
\succ 0.$$ With the above approach we solved problem for the same system data as in [@Hassibi1999].  Here, matrices $A$, $B$ and $C$ are given as: $$\begin{aligned}
A \!\!=\!\! { \scriptsize\left[\!\begin{array}{rrrrr}
   \!\!\!\!  2.45  \!\!\!\!&\!\!  -0.90 \!\!\!&\!\!   1.53 \!\!\!&\!\!  -1.26 \!\!\!&\!\!   1.76\!\!\\
   \!\!\!\!  -0.12 \!\!\!\!&\!\!  -0.44 \!\!\!&\!\!  -0.01 \!\!\!&\!\!   0.69 \!\!\!&\!\!   0.90\!\!\\
   \!\!\!\!  2.07  \!\!\!\!&\!\!  -1.20 \!\!\!&\!\!  -1.14 \!\!\!&\!\!   2.04 \!\!\!&\!\!  -0.76\!\!\\
   \!\!\!\!  -0.59 \!\!\!\!&\!\!   0.07 \!\!\!&\!\!   2.91 \!\!\!&\!\!  -4.63 \!\!\!&\!\!  -1.15\!\!\\
   \!\!\!\!  -0.74 \!\!\!\!&\!\!  -0.23 \!\!\!&\!\!  -1.19 \!\!\!&\!\!  -0.06 \!\!\!&\!\!  -2.52\!\! 
\end{array}\!\right]}\!, \ 
B \!\!=\!\! { \scriptsize\left[\!\begin{array}{rrrrr}
   \!\!\!  0.81 \!\!&\!\!  -0.79 \!\!&\!\!   0.00 \!\!&\!\!   0.00 \!\!&\!\!\!  -0.95\!\!\\
   \!\!\!-0.34 \!\!&\!\!  -0.50 \!\!&\!\!   0.06 \!\!&\!\!   0.22 \!\!&\!\!\!   0.92\!\!\\
   \!\!\!-1.32 \!\!&\!\!   1.55 \!\!&\!\!  -1.22 \!\!&\!\!  -0.77 \!\!&\!\!\!  -1.14\!\!\\
   \!\!\!-2.11 \!\!&\!\!   0.32 \!\!&\!\!   0.00 \!\!&\!\!  -0.83 \!\!&\!\!\!   0.59\!\!\\
   \!\!\!  0.31 \!\!&\!\!  -0.19 \!\!&\!\!  -1.09 \!\!&\!\!   0.00 \!\!&\!\!\!   0.00\!\! 
\end{array}\!\right]}, \
C \!\!=\!\! { \scriptsize\begin{bmatrix}
    0.00 \!\!&\!\!   0.00 \!\!&\!\!   0.16 \!\!&\!   0.00 \!\!\!&\!\!  -1.78\!\\
    1.23 \!\!&\!\!  -0.38 \!\!&\!\!   0.75 \!\!&\!  -0.38 \!\!\!&\!\!  -0.00\!\\
    0.46 \!\!&\!\!   0.00 \!\!&\!\!  -0.05 \!\!&\!   0.00 \!\!\!&\!\!   0.00\!\\
    0.00 \!\!&\!\!  -0.12 \!\!&\!\!  0.23 \!\!&\!  -0.12 \!\!\!&\!\!   1.14\! 
\end{bmatrix}}\!.\end{aligned}$$ The weighting parameter $\sigma$ is chosen by $\sigma = 3$. Algorithm \[alg:A1\] is terminated if one of the following conditions is satisfied:

-   subproblem encounters a numerical problem;

-   ${\|\Delta x^k\|}_{\infty}/({\|x^k\|}_{\infty}+1)\leq 10^{-3}$;

-   the maximum number of iterations, $K_{\max}$, reaches;

-   or the objective function is not significantly improved after two successive iterations (i.e. $|f^{k+1}-f^k| \leq 10^{-4}(1+|f^k|)$, for some $k=\bar{k}$ and $k=\bar{k}+1$, where $f^k :=
    f(x^k)$). 

In this example, Algorithm \[alg:A1\] is terminated after $15$ iterations, whereas the objective function is not significantly improved.  However, after the $2^{\textrm{nd}}$ iteration, matrix $F$ only has $3$ nonzero elements, while the decay rate $\alpha$ is $1.17316$. This value is much higher than the one reported in [@Hassibi1999], $\alpha = 0.3543$ after $6$ iterations.  We obtain the gain matrix $F$ as $$F = {\footnotesize\begin{bmatrix}
 0.6540 &  0.0000 & 0.0000 & 0.0000 \\
 0.0000 & -0.4872 & 0.0000 & 0.0000 \\
 0.0000 &  0.0000 & 0.0000 & 0.0000 \\
 0.0000 &  0.0000 & 0.0000 & 0.0000 \\
 0.0000 &  1.1280 & 0.0000 & 0.0000
\end{bmatrix}.}$$ With this matrix, the maximum real part of the eigenvalues of the closed-loop matrix in (\[eq:ss\_LTI\]), $A_F := A + BFC$, is $\alpha_0(A_F) := -1.40706$. Simultaneously, $\alpha_0(A_F^TP+PA_F+2\alpha P) = -0.327258 < 0$ and $\alpha_0(P) = 0.587574 > 0$. Note that $\alpha_0(A_F) \neq -\alpha$ due to the in-activeness of the BMI constraint in at the 2nd iteration step. 

Spectral abscissa and pseudo-spectral abscissa optimization {#subsec:SOF_design}
-----------------------------------------------------------

One popular problem in control theory is to optimize the spectral abscissa of the closed-loop system $\dot{x} = (A+BFC)x$. Briefly, this problem is presented as an unconstrained optimization problem of the form: $$\label{eq:pure_spectral_abscissa_opt}
\min_{F\in\mathbb{R}^{n_u\times n_y}} \alpha_0(A+BFC),$$ where $\alpha_0(A+BFC) := \sup\left\{ \textrm{Re}(\lambda) ~|~ \lambda\in\Lambda(A+BFC)\right\}$ is the spectral abscissa of $A+BFC$, $\textrm{Re}(\lambda)$ denotes the real part of $\lambda\in\mathbb{C}$ and $\Lambda(A+BFC)$ is the spectrum of $A+BFC$. Problem has many drawbacks in terms of numerical solution due to the nonsmoothness and non-Lipschitz continuity of the objective function $\alpha_0$ [@Burke2002]. 

In order to apply the method developed in this paper, we reformulate problem as an optimization problem with BMI constraints [@Burke2002; @Leibfritz2004]: $$\label{eq:BMI_p2}
\left\{\begin{array}{cl}
\displaystyle
\max_{P,F,\beta} & \beta\\
\textrm{s.t.}~ &(A+BFC)^TP + P(A+BFC) + 2\beta P \prec 0, ~~P = P^T, ~P \succ 0. 
\end{array}\right.$$ Here, matrices $A\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times n_u}$, $C\in\mathbb{R}^{n_y\times n}$ are given.  Matrices $P\in\mathbb{R}^{n\times n}$ and $F\in\mathbb{R}^{n_u\times n_y}$ and the scalar $\beta$ are considered as variables.  If the optimal value of is strictly positive then the closed-loop feedback controller $u = Fy$ stabilizes the linear system $\dot{x} =
(A+BFC)x$.

Problem is very similar to .  Therefore, using the same trick as in , we can reformulate in the form of .  More precisely, if we define $B_F := A + BFC + \beta I$ then the bilinear matrix mapping $A_F^TP + PA_F$ can be represented as a psd-convex-concave decomposition of the form and problem can be rewritten in the form of .  We implement Algorithm \[alg:A1\] for solving this resulting problem using the same parameters and the stopping criterions as in Subsection \[subsec:sparse\_SOF\].  In addition, we regularize the objective function by adding the term $\frac{\rho_F}{2}{\|F-F^k\|}_F^2 + \frac{\rho_P}{2}{\|P-P^k\|}_F^2$, with $\rho_F
= \rho_P = 10^{-2}$. The maximum number of iterations $K_{\max}$ is set to $150$.

We test for several problems in COMP$\textrm{l}_\textrm{e}$ib and compare our results with the ones reported by HIFOO, PENBMI and LMIRank.  For LMIRank, we implement the algorithm proposed in [@Orsi2006].  We initialize the value of the decay rate $\alpha^0$ at $10^{-4}$ and perform an iterative loop to increase $\alpha$ as $\alpha^{k+1} := \alpha^k + 0.1$. The algorithm is terminated if either the problems (12) or (21) in [@Orsi2006] with a correspondence $\alpha$ can not be solved or the maximum number of iterations $K_{\max} = 100$ is reached. 

  ---------------------------------------- --------------- ------------------ ------------------ ------------------ ----------------- --------- ------------ -- --
  `Name`                                     $\alpha_0(A)$              HIFOO            LMIRANK             PENBMI   $\alpha_0(A_F)$    `Iter`    `time[s]`    
  [AC1]{}                                        [0.000]{}    [**-0.2061**]{}    [**-8.4766**]{}    [**-7.0758**]{}   [**-0.8535**]{}   [ 41]{}    [12.44]{}    
  81 [AC4]{}                                     [2.579]{}    [**-0.0500**]{}    [**-0.0500**]{}    [**-0.0500**]{}   [**-0.0500**]{}   [ 14]{}     [4.60]{}    
  [**AC5**$^{\textbf{a}}$]{}                     [0.999]{}    [**-0.7746**]{}    [**-1.8001**]{}    [**-2.0438**]{}   [**-0.7389**]{}   [ 28]{}    [63.33]{}    
  [AC7]{}                                        [0.172]{}    [**-0.0322**]{}    [**-0.0204**]{}     [**0.0896**]{}   [**-0.0673**]{}   [150]{}   [111.46]{}    
  2 [AC8]{}                                      [0.012]{}    [**-0.1968**]{}    [**-0.4447**]{}     [**0.4447**]{}   [**-0.0755**]{}   [ 24]{}    [21.95]{}    
  6 [AC9]{}                                      [0.012]{}    [**-0.3389**]{}    [**-0.5230**]{}    [**-0.4450**]{}   [**-0.3256**]{}   [ 78]{}    [74.57]{}    
  7 [AC11]{}                                     [5.451]{}    [**-0.0003**]{}    [**-5.0577**]{}          [**-**]{}   [**-3.0244**]{}   [ 61]{}    [38.44]{}    
  52 [AC12]{}                                    [0.580]{}   [**-10.8645**]{}    [**-9.9658**]{}    [**-1.8757**]{}   [**-0.3414**]{}   [150]{}    [86.72]{}    
  100 [HE1]{}                                    [0.276]{}    [**-0.2457**]{}    [**-0.2071**]{}    [**-0.2468**]{}   [**-0.2202**]{}   [150]{}    [87.64]{}    
  32.8000, 4 [HE3]{}                             [0.087]{}    [**-0.4621**]{}    [**-2.3009**]{}    [**-0.4063**]{}   [**-0.8702**]{}   [ 47]{}    [48.80]{}    
  700.8300, 25 [HE4]{}                           [0.234]{}    [**-0.7446**]{}    [**-1.9221**]{}    [**-0.0909**]{}   [**-0.8647**]{}   [ 63]{}    [71.66]{}    
  21 [HE5]{}                                     [0.234]{}    [**-0.1823**]{}          [**-**]{}    [**-0.2932**]{}   [**-0.0587**]{}   [150]{}   [178.71]{}    
  [HE6]{}                                        [0.234]{}    [**-0.0050**]{}    [**-0.0050**]{}    [**-0.0050**]{}   [**-0.0050**]{}   [ 12]{}    [41.00]{}    
  2 [REA1]{}                                     [1.991]{}   [**-16.3918**]{}    [**-5.9736**]{}    [**-1.7984**]{}   [**-3.8599**]{}   [ 77]{}    [79.23]{}    
  61 [REA2]{}                                    [2.011]{}    [**-7.0152**]{}   [**-10.0292**]{}    [**-3.5928**]{}   [**-2.1778**]{}   [ 40]{}    [39.15]{}    
  [REA3]{}                                       [0.000]{}    [**-0.0207**]{}    [**-0.0207**]{}    [**-0.0207**]{}   [**-0.0207**]{}   [150]{}   [362.21]{}    
  [DIS2]{}                                       [1.675]{}    [**-6.8510**]{}   [**-10.1207**]{}    [**-8.3289**]{}   [**-8.4540**]{}   [ 28]{}    [37.28]{}    
  [DIS4]{}                                       [1.442]{}   [**-36.7203**]{}    [**-0.5420**]{}   [**-92.2842**]{}   [**-8.0989**]{}   [ 72]{}   [124.23]{}    
  [WEC1]{}                                       [0.008]{}    [**-8.9927**]{}    [**-8.7350**]{}    [**-0.9657**]{}   [**-0.8779**]{}   [150]{}   [305.36]{}    
  [IH]{}                                         [0.000]{}    [**-0.5000**]{}    [**-0.5000**]{}    [**-0.5000**]{}   [**-0.5000**]{}    [ 7]{}    [39.41]{}    
  [CSE1]{}                                       [0.000]{}    [**-0.4509**]{}    [**-0.4844**]{}    [**-0.4490**]{}   [**-0.2360**]{}   [ 38]{}   [158.67]{}    
  [TF1]{}                                        [0.000]{}          [**-**]{}          [**-**]{}    [**-0.0618**]{}   [**-0.1544**]{}   [ 56]{}   [137.98]{}    
  [TF2]{}                                        [0.000]{}          [**-**]{}          [**-**]{}    [**-1.0e-5**]{}   [**-1.0e-5**]{}    [ 8]{}    [20.41]{}    
  [TF3]{}                                        [0.000]{}          [**-**]{}          [**-**]{}    [**-0.0032**]{}   [**-0.0031**]{}   [ 93]{}   [237.93]{}    
  [NN1]{}                                        [3.606]{}    [**-3.0458**]{}    [**-4.4021**]{}    [**-4.3358**]{}   [**-0.8746**]{}   [ 12]{}    [37.53]{}    
  53.1700, 65 [**NN5**$^{\textbf{a}}$]{}         [0.420]{}    [**-0.0942**]{}    [**-0.0057**]{}    [**-0.0942**]{}   [**-0.0913**]{}   [ 11]{}    [42.62]{}    
  [NN9]{}                                        [3.281]{}    [**-2.0789**]{}    [**-0.7048**]{}          [**-**]{}   [**-0.0279**]{}   [ 33]{}   [111.41]{}    
  [NN13]{}                                       [1.945]{}    [**-3.2513**]{}    [**-4.5310**]{}    [**-9.0741**]{}   [**-3.4318**]{}   [150]{}   [572.50]{}    
  [NN15]{}                                       [0.000]{}    [**-6.9983**]{}   [**-11.0743**]{}    [**-0.0278**]{}   [**-0.8353**]{}   [150]{}   [524.80]{}    
  [NN17]{}                                       [1.170]{}    [**-0.6110**]{}    [**-0.5130**]{}          [**-**]{}   [**-0.6008**]{}   [ 99]{}   [342.67]{}    
  ---------------------------------------- --------------- ------------------ ------------------ ------------------ ----------------- --------- ------------ -- --

  : Computational results for in COMP$\textrm{l}_{\textrm{e}}$ib[]{data-label="tb:exam2"}

-0.4cm

The numerical results of four algorithms are reported in Table \[tb:exam2\].  Here, we initialize the algorithm in HIFOO with the same initial guess $F^0 = 0$. Since PENBMI and our methods solve the same BMI problems, they are initialized by the same initial values for $P$, $F$ and $\beta$.

The notation in Table \[tb:exam2\] consists of: `Name` is the name of problems, $\alpha_0(A)$, $\alpha_0(A_F)$ are the maximum real part of the eigenvalues of the open-loop and closed-loop matrices $A$, $A_F$, respectively; `iter` is the number of iterations, `time[s]` is the CPU time in second.  The columns titled HIFOO, LMIRank and PENBMI give the maximum real part of the eigenvalues of the closed-loop system for a static output feedback controller computed by available software HIFOO [@Gumussoy09], LMIRank [@Orsi2006] and PENBMI [@Henrion2005], respectively.  Our results can be found in the sixth column.  The entries with a dash sign indicate that there is no feasible solution found.  Algorithm $1$ fails or makes only slow progress toward to a local solution with $6$ problems: AC18, DIS5, PAS, NN6, NN7, NN12 in COMP$\textrm{l}_\textrm{e}$ib.  Problems AC5 and NN5 are initialized with a different matrix $F^0$ to avoid numerical problems. 

Note that Algorithm \[alg:A1\] as well as the algorithms implemented in HIFOO, LMIRank and PENBMI are local optimization methods, which only report a local minimizer and these solutions may not be the same.  To apply the LMIRank package for solving problem , we have used a direct search procedure for finding $\alpha$. The computational time of this procedure is very high compared with the other methods. 

To conclude this subsection, we show that our method can be applied to solve the optimization problem of pseudo-spectral abscissa in static feedback controller designs.  This problem is described as follows (see [@Burke2002; @Leibfritz2004]): $$\label{eq:BMI_p3b}
\begin{array}{cl}
&\displaystyle\max_{\beta,\mu,\omega,F,P}\beta\\
& \textrm{s.t.}~ \begin{bmatrix} 2\beta P \!+\!  A_F^T\!\!P \!+\!  PA_F \!+\!  \mu I \!-\!  \omega I_z \!\!\! & \!\!\!  \varepsilon P\\ \varepsilon P
&\omega I \end{bmatrix} \!\preceq \!  0, ~ P\succ 0, ~ P = P^T, ~ \mu > 0,
\end{array}$$ where $A_F = A + BFC$ as before and $\omega\leq 0$.

Using the same notation $B_F = A + BFC + \beta I$ as in and applying the statement b) of Lemma \[le:key\_lemma\], the BMI constraint in this problem can be transformed into a psd-convex-concave one: $$\begin{aligned}
\label{eq:BMI_reform}
&\begin{bmatrix}
B_F^T\!\!B_F \!+\!  P^T\!\!P \!+\! ( \mu \!-\!  \omega)I \!\!&\!\!  \epsilon P\\ \epsilon P \!\!&\!\!  \omega I
\end{bmatrix}
\!-\!  \begin{bmatrix}(B_F\!-\!P)^T\!\!(B_F\!-\!P) \!\!&\!\!  0\\ 0 \!\!&\!\!  0\end{bmatrix} \!\preceq \!0.\end{aligned}$$ If we denote the linearization of $(B_F-P)^T(B_F-P)$ at the iteration $k$ by $H_k$, i.e. $H_k = (B_F-P)^T(B_F^k-P^k) + (B_F^k -P^k)^T(B_F-P) -
(B_F^k-P^k)^T(B_F^k-P^k)$, then the linearized constraint in the subproblem can be represented as an LMI thanks to Lemma \[le:key\_lemma\]: $$\begin{bmatrix}
H_k + (\omega-\mu)I & B_F^T & P & -\varepsilon P\\ B_F & I & 0 & 0\\
P & 0 & I & 0\\ -\varepsilon P & 0 & 0 & -\omega I
\end{bmatrix} \succeq 0.$$ Hence, Algorithm \[alg:A1\] can be applied to solve problem . 

If we define $\bar{F} := BFC$ then the bilinear matrix mapping $A_F^TP + PA_F$ can be rewritten as $$\begin{aligned}
A_F^TP + PA_F &= \frac{1}{2}\left[(P+\bar{F})^T(P+\bar{F}) - (P-\bar{F})^T(P-\bar{F})\right] - A^TP - PA.\end{aligned}$$ Using this decomposition, one can avoid the contribution of matrix $A$ on the bilinear term.  Consequently, Algorithm \[alg:A1\] may work better in some specific problems. 

$\mathcal{H}_2$ optimization: BMI formulation {#subsec:H2_prob}
---------------------------------------------

In this subsection, we consider an optimization problem arising in the $\mathcal{H}_2$ synthesis of the linear system .  Let us assume that $D_{12} = 0$ and $D_{21} = 0$, then this problem is formulated as the following optimization problem with BMI constraints [@Leibfritz2003].  $$\label{eq:BMI_H2_prob}
\begin{array}{cl}
\displaystyle\min_{F, Q, X} &\textrm{trace}(X)\\
\textrm{s.t.}\!\!\!\!\!\!\!\! &  ~ (A+BFC)Q + Q(A+BFC)^T + B_1B_1^T \preceq 0, ~~\begin{bmatrix}X & C_1Q\\ QC_1^T & Q \end{bmatrix}\succeq 0, ~~ Q \succ 0. 
\end{array}$$ Here, we also assume that $B_1B_1^T$ is positive definite.  Otherwise, we use $B_1B_1^T + \varepsilon I$ instead of $B_1B_1^T$ with $\varepsilon = 10^{-5}$ in . 

In order to apply Algorithm \[alg:A1\] for solving problem , a starting point $x^0\in\textrm{ri}(\mathcal{D})$ is required.  This task can be done by performing some extra steps called *Phase 1*.  The algorithm is now split in two phases as follows. 

**Phase 1:** (*Determine a starting point $x^0$*). 

-   *Step 1.*  If $\alpha_0(A + A^T) < 0$ then we set $F^0 := 0$. Otherwise, go to Step 3. 

-   *Step 2.*  Solve the following optimization problem with LMI constraints: $$\label{eq:LMI_H2_prob}
    \begin{array}{cl}
    \displaystyle\min_{Q, X} \textrm{trace}(X) ~~~~
    \textrm{s.t.}  ~~ A_{F^0}Q + QA_{F^0}^T + B_1B_1^T \prec 0, ~~\begin{bmatrix}X & C_1Q\\ QC_1^T & Q \end{bmatrix}\succ 0, ~~Q \succ 0,
    \end{array}$$ where $A_{F^0} := A + BF^0C$. If this problem has a solution $Q^0$ and $X^0$ then terminate Phase 1 and using $F^0$ together with $Q^0, X^0$ as a starting point $x^0$ for Phase 2.  Otherwise, go to Step 3. 

-   *Step 3.*  Solve the following feasibility problem with LMI constraints: $$\label{eq:LMI_H2_phase1}
    \begin{array}{cl}
    \textrm{Find}& \textrm{$P\succ 0$ and $K$ such that:}\\
    &\begin{bmatrix}PA+A^TP+KC+C^TK^T & PB_1\\ B_1^TP &-\sigma_0^2I_{w}\end{bmatrix} \preceq 0, ~~\begin{bmatrix}X & C_1\\ C_1^T & P\end{bmatrix} \succeq 0,
    \end{array}$$ to obtain $K^{*}$ and $P^{*}$, where $\sigma_0$ is a given regularization factor.  Compute $F^{*} := B^{+}(P^{*})^{-1}K^{*}$, where $B^{+}$ is a pseudo-inverse of $B$, and resolve problem with $F^0 := F^{*}$. If problem has a solution $Q^0$ and $X^0$ then terminate Phase 1 and set $x^0 := (F^0, Q^0, X^0)$. Otherwise, perform Step 4. 

-   *Step 4.*  Apply the method in Subsection \[subsec:SOF\_design\] to solve the following BMI feasibility problem: $$\begin{aligned}
    \label{eq:BMI_H2_phase1}
    \textrm{Find }& \textrm{$F$ and $Q\succ 0$ such that}: ~~ (A+BFC)Q + Q(A+BFC)^T + B_1B_1^T \prec 0.\end{aligned}$$ If this problem has a solution $F^0$ then go back to Step 2.  Otherwise, declare that no strictly feasible point is found. 

**Phase 2:** (*Solve problem* ).  Perform Algorithm \[alg:A1\] with the starting point $x^0$ found at Phase 1. [ $\square$]{}

Note that Step 3 of Phase 1 corresponds to determining a full state feedback controller and approximating it subsequently with an output feedback controller.  Step 4 of Phase 1 is usually expensive.  Therefore, in our numerical implementation, we terminate Step 4 after finding a point such that $\alpha_0((A+BFC)Q + Q(A+BFC)^T + B_1B_1^T) \leq -0.1$.

\[re:stop\_phase1\] The algorithm described in Phase 1 is finite.  It terminates either at Step 4 if no feasible point is found or at Step 2 if a feasible point is found.  Indeed, if a feasible matrix $F^0$ is found at Step 4, the first BMI constraint of is feasible with some $Q \succ 0$. Thus we can find an appropriate matrix $X$ such that $X -
CQC^T \succ 0$, which implies the second LMI constraint of is satisfied.  Consequently, problem has a solution. 

The method used in Phase 1 is heuristic.  It can be improved when we apply to a certain problem.  However, as we can see in the numerical results, it performs quite acceptable for majority of the test problems. 

In the following numerical examples, we implement Phase 1 and Phase 2 of the algorithm using the decomposition $$\begin{aligned}
A_FQ + QA_F^T + B_1B_1^T &= \frac{1}{2}(A_F+Q)(A_F+Q)^T  + B_1B_1^T - \frac{1}{2}(A_F-Q)(A_F-Q)^T\end{aligned}$$ for the BMI form at the left-top corner of the first constraint in .  The regularization parameters and the stopping criterion for Algorithm \[alg:A1\] are chosen as in Subsection \[subsec:sparse\_SOF\] and $K_{\max} = 300$.

  ---------------- -------- ------- ------- -------- -------- ----------------- ----------------- -------------------- --------- ------------- --
  `Name`            $n_x$    $n_y$   $n_u$   $n_z$    $n_w$       HIFOO                 PENBMI      $\mathcal{H}_2$       `iter`     `time`    
  [AC1$^{b}$]{}     [5]{}    [3]{}   [3]{}   [2]{}    [3]{}    [**0.0250**]{}      [**0.0061**]{}       [**0.0540**]{}     [3]{}     [3.380]{} 
  [AC2$^{b}$]{}     [5]{}    [3]{}   [3]{}   [5]{}    [3]{}    [**0.0257**]{}      [**0.0075**]{}       [**0.0540**]{}     [3]{}     [3.390]{} 
  [AC3]{}           [5]{}    [4]{}   [2]{}   [5]{}    [5]{}    [**2.0964**]{}      [**2.0823**]{}       [**2.1117**]{}   [210]{}    [73.380]{} 
  [AC4]{}           [4]{}    [2]{}   [1]{}   [2]{}    [2]{}    [**11.0269**]{}          [**-**]{}      [**11.0269**]{}     [2]{}     [0.990]{} 
  [AC6]{}           [7]{}    [4]{}   [2]{}   [7]{}    [7]{}    [**2.8648**]{}      [**2.8648**]{}       [**2.8664**]{}   [153]{}   [124.230]{} 
  [AC7]{}           [9]{}    [2]{}   [1]{}   [1]{}    [4]{}    [**0.0172**]{}      [**0.0162**]{}       [**0.0176**]{}     [1]{}     [0.670]{} 
  [AC8]{}           [9]{}    [5]{}   [1]{}   [2]{}    [10]{}   [**0.6330**]{}      [**0.7403**]{}       [**0.6395**]{}   [300]{}   [282.760]{} 
  [AC12$^{b}$]{}    [4]{}    [4]{}   [3]{}   [1]{}    [3]{}    [**0.0022**]{}      [**0.0106**]{}       [**0.0992**]{}    [36]{}    [28.540]{} 
  [AC15$^{b}$]{}    [4]{}    [3]{}   [2]{}   [6]{}    [4]{}    [**1.5458**]{}      [**1.4811**]{}       [**1.5181**]{}   [264]{}    [85.390]{} 
  [AC16$^{b}$]{}    [4]{}    [4]{}   [2]{}   [6]{}    [4]{}    [**1.4769**]{}      [**1.4016**]{}       [**1.4427**]{}   [300]{}    [99.790]{} 
  [AC17]{}          [4]{}    [2]{}   [1]{}   [4]{}    [4]{}    [**1.5364**]{}      [**1.5347**]{}       [**1.5507**]{}   [171]{}    [49.350]{} 
  [HE2]{}           [4]{}    [2]{}   [2]{}   [4]{}    [4]{}    [**3.4362**]{}      [**3.4362**]{}       [**4.7406**]{}   [263]{}    [97.310]{} 
  [HE3$^{b}$]{}     [8]{}    [6]{}   [4]{}   [10]{}   [1]{}    [**0.0197**]{}      [**0.0071**]{}       [**0.1596**]{}   [249]{}   [217.360]{} 
  [HE4$^{b}$]{}     [8]{}    [6]{}   [4]{}   [12]{}   [8]{}    [**6.6436**]{}      [**6.5785**]{}       [**7.1242**]{}   [300]{}   [412.830]{} 
  [REA1]{}          [4]{}    [3]{}   [2]{}   [4]{}    [4]{}    [**0.9442**]{}      [**0.9422**]{}       [**1.0622**]{}   [249]{}    [80.810]{} 
  [REA2$^{b}$]{}    [4]{}    [2]{}   [2]{}   [4]{}    [4]{}    [**1.0339**]{}      [**1.0229**]{}       [**1.1989**]{}   [300]{}   [101.730]{} 
  [DIS1]{}          [8]{}    [4]{}   [4]{}   [8]{}    [1]{}    [**0.6705**]{}      [**0.1174**]{}       [**0.7427**]{}   [300]{}   [255.810]{} 
  [DIS2]{}          [3]{}    [2]{}   [2]{}   [3]{}    [3]{}    [**0.4013**]{}      [**0.3700**]{}       [**0.3819**]{}     [4]{}     [1.370]{} 
  [DIS3]{}          [6]{}    [4]{}   [4]{}   [6]{}    [6]{}    [**0.9527**]{}      [**0.9434**]{}       [**1.0322**]{}   [300]{}   [210.470]{} 
  [DIS4]{}          [6]{}    [6]{}   [4]{}   [6]{}    [6]{}    [**1.0117**]{}      [**0.9696**]{}       [**1.0276**]{}   [300]{}   [210.690]{} 
  [WEC1$^{b}$]{}    [10]{}   [4]{}   [3]{}   [10]{}   [10]{}   [**7.3940**]{}      [**8.1032**]{}      [**12.9093**]{}   [119]{}   [190.150]{} 
  [WEC2$^{b}$]{}    [10]{}   [4]{}   [3]{}   [10]{}   [10]{}   [**6.7908**]{}      [**7.6502**]{}      [**12.2102**]{}   [261]{}   [407.470]{} 
  [AGS]{}           [12]{}   [2]{}   [2]{}   [12]{}   [12]{}   [**6.9737**]{}      [**6.9737**]{}       [**6.9838**]{}    [18]{}    [28.830]{} 
  [BDT1]{}          [11]{}   [3]{}   [3]{}   [6]{}    [1]{}    [**0.0024**]{}           [**-**]{}       [**0.0017**]{}    [51]{}    [64.410]{} 
  [MFP]{}           [4]{}    [2]{}   [3]{}   [4]{}    [4]{}    [**6.9724**]{}      [**6.9724**]{}       [**7.0354**]{}   [300]{}   [114.560]{} 
  [PSM]{}           [7]{}    [3]{}   [2]{}   [5]{}    [2]{}    [**0.0330**]{}      [**0.0007**]{}       [**0.1753**]{}   [300]{}   [217.250]{} 
  [EB2$^{b}$]{}     [10]{}   [1]{}   [1]{}   [2]{}    [2]{}    [**0.0640**]{}      [**0.0084**]{}       [**0.1604**]{}   [114]{}   [131.380]{} 
  [EB3]{}           [10]{}   [1]{}   [1]{}   [2]{}    [2]{}    [**0.0732**]{}      [**0.0072**]{}       [**0.0079**]{}     [7]{}     [6.240]{} 
  [TF1$^{b}$]{}     [7]{}    [4]{}   [2]{}   [4]{}    [1]{}    [**0.0945**]{}           [**-**]{}       [**0.1599**]{}   [192]{}   [166.810]{} 
  [TF2]{}           [7]{}    [3]{}   [2]{}   [4]{}    [1]{}    [**11.1803**]{}          [**-**]{}      [**11.1803**]{}     [3]{}     [2.810]{} 
  [TF3$^{b}$]{}     [7]{}    [3]{}   [2]{}   [4]{}    [1]{}    [**0.1943**]{}      [**0.1424**]{}       [**0.2565**]{}   [138]{}   [128.250]{} 
  [NN2]{}           [2]{}    [1]{}   [1]{}   [2]{}    [2]{}    [**1.1892**]{}      [**1.1892**]{}       [**1.1892**]{}     [4]{}     [1.090]{} 
  [NN4]{}           [4]{}    [3]{}   [2]{}   [4]{}    [4]{}    [**1.8341**]{}      [**1.8335**]{}       [**1.8590**]{}   [222]{}    [67.260]{} 
  [NN8]{}           [3]{}    [2]{}   [2]{}   [3]{}    [3]{}    [**1.5152**]{}      [**1.5117**]{}       [**1.5725**]{}   [183]{}    [50.170]{} 
  [NN11]{}          [16]{}   [5]{}   [3]{}   [3]{}    [3]{}    [**0.1178**]{}      [**0.0790**]{}       [**0.1263**]{}    [39]{}    [91.070]{} 
  [NN13$^{b}$]{}    [6]{}    [2]{}   [2]{}   [3]{}    [3]{}    [**26.1012**]{}    [**26.1314**]{}      [**62.3995**]{}   [138]{}   [112.750]{} 
  [NN14$^{b}$]{}    [6]{}    [2]{}   [2]{}   [3]{}    [3]{}    [**26.1448**]{}    [**26.1314**]{}      [**62.3995**]{}   [138]{}   [110.020]{} 
  [NN15]{}          [3]{}    [2]{}   [2]{}   [4]{}    [1]{}    [**0.0245**]{}           [**-**]{}       [**0.0210**]{}     [6]{}     [2.060]{} 
  [NN16$^{b}$]{}    [8]{}    [4]{}   [4]{}   [4]{}    [8]{}    [**0.1195**]{}      [**0.1195**]{}       [**0.1195**]{}     [3]{}    [23.030]{} 
  [NN17]{}          [3]{}    [1]{}   [2]{}   [2]{}    [1]{}    [**3.2530**]{}      [**3.2404**]{}       [**3.3329**]{}   [300]{}    [88.770]{} 
  ---------------- -------- ------- ------- -------- -------- ----------------- ----------------- -------------------- --------- ------------- --

  : $\mathcal{H}_2$ synthesis benchmarks on COMP$\textrm{l}_{\textrm{e}}$ib plants[]{data-label="tb:H2_prob"}

-0.4cm

We test the algorithm for many problems in COMP$\textrm{l}_{\textrm{e}}$ib and the computational results are reported in Table \[tb:H2\_prob\].  For the comparison purpose, we also carry out the test with HIFOO [@Gumussoy09] and PENBMI [@Henrion2005], and the results are put in the columns marked by HIFOO and PENBMI in Table \[tb:H2\_prob\], respectively.  The initial controller for HIFOO is set to $F^0$ and the BMI parameters for PENBMI are initialized with $(F, Q, X)=(F^0, Q^0,
X^0)$. Here, $n, n_y, n_z, n_w, n_u$ are the dimensions of problems, the columns titled HIFOO and PENBMI give the $\mathcal{H}_2$ norm of the closed-loop system for the static output feedback controller computed by HIFOO and PENBMI; `iter` and `time[s]` are the number of iterations and CPU time in second of Algorithm \[alg:A1\], respectively, included Phase 1 and Phase 2.  Problems marked by “b” mean that Step 4 in Phase 1 is performed.  In Table \[tb:H2\_prob\], we only report the problems that were solved by Algorithm \[alg:A1\].  The numerical results allow us to conclude that Algorithm \[alg:A1\], PENBMI and HIFOO report similar values for majority of the test problems in COMP$\textrm{l}_{\textrm{e}}$ib. 

If $D_{12} \neq 0$ then the second LMI constraint of becomes a BMI constraint: $$\label{eq:BMI_H2_const2}
\begin{bmatrix} X & (C_1 + D_{12}FC)Q \\ Q(C_1+D_{12}FC)^T & Q\end{bmatrix} \succeq 0,$$ which is equivalent to $X - C_FQC_F^T \succeq 0$, where $C_F := C_1 + D_{12}FC$. Since $f(Q) := Q^{-1}$ is convex on $\mathbf{S}^{n_x}_{++}$ (see Lemma \[le:key\_lemma\_a\] a)), this BMI constraint can be reformulated as a convex-concave matrix inequality constraint of the form: $$\label{eq:BMI_H2_const2_DC}
\begin{bmatrix} X & C_F \\ C_F^T & O\end{bmatrix} + \begin{bmatrix}O & O\\ O & Q^{-1}\end{bmatrix}\succeq 0.$$ By linearizing the concave term $-f(Q)$ at $Q = Q^k$ as $(Q^k)^{-1} - (Q^{k})^{-1}(Q-Q^k)(Q^k)^{-1}$ (see [@Boyd2004]), the resulting constraint can be written as an LMI constraint.  Therefore, Algorithm \[alg:A1\] can be applied to solve problem in the case $D_{12}\neq 0$.

$\mathcal{H}_{\infty}$ optimization: BMI formulation {#subsec:Hinf_prob}
----------------------------------------------------

Alternatively, we can also apply Algorithm \[alg:A1\] to solve the optimization with BMI constraints arising in $\mathcal{H}_{\infty}$ optimization of the linear system .  Let us assume that $D_{21} = 0$, then this problem is reformulated as the following optimization problem with BMI constraints [@Leibfritz2003]: $$\label{eq:BMI_Hinf_prob}
\begin{array}{cl}
\displaystyle\min_{F, X, \gamma} &\gamma \\
\textrm{s.t.}\!\!\!\!\!\!\!\! &  ~\begin{bmatrix}A_F^TX + XA_F &  XB_1 & C_F^T\\ B_1^TX & -\gamma I_w & D_{11}^T\\ C_F & D_{11} & -\gamma I_z\end{bmatrix} \prec
0, ~~ X \succ 0, ~ \gamma > 0. 
\end{array}$$ Here, as before, $A_F = A + BFC$ and $C_F = C_1 + D_{12}FC$. The bilinear matrix term $A_F^TX + XA_F$ at the top-corner of the last constraint can be decomposed as or .  Therefore, we can use these decompositions to transform problem into .  After linearization, the resulting subproblem is also rewritten as a standard SDP problem by applying Lemma \[le:key\_lemma\].  We omit this specification here. 

To determine a starting point, we perform Phase 1 which is similar to the one carried out in the ${\mathcal{H}}_2$-optimization subsection. 

**Phase 1.** (* Determine a starting point $x^0\in\textrm{ri}(\mathcal{D})$*). 

-   *Step 1.*  If $\alpha_0(A^T+A) < 0$ then set $F^0 = 0$. Otherwise, go to Step 3. 

-   *Step 2.*  Solve the following optimization with LMI constraints $$\label{eq:LMI_Hinf_prob}
    \begin{array}{cl}
    \displaystyle\min_{\gamma, X} &\gamma ~~~ \textrm{s.t.}  ~~~ 
    \begin{bmatrix}A_{F^0}^TX + XA_{F^0} & XB_1 & C_{F^0}^T\\ B_1^TX & -\gamma I_w & D_{11}^T\\ C_{F^0} & D_{11} & -\gamma I_z \end{bmatrix} \!\prec 0, ~~X \succ 0,
    ~\gamma > 0,
    \end{array}$$ where $A_{F^0} := A + BF^0C$ and $C_{F^0} := C_1 + D_{12}F^0C$. If this problem has a solution $\gamma^0$ and $X^0$ then terminate Phase 1 and using $F^0$ together with $\gamma^0, X^0$ as a starting point $x^0$ for Phase 2.  Otherwise, go to Step 3. 

-   *Step 3.*  Solve the following feasibility problem of LMI constraints: $$\label{eq:LMI_Hinf_phase1}
    \begin{array}{cl}
    \textrm{Find $P\succ 0$, $\gamma>0$ and $K$ such that:} \begin{bmatrix}\!  PA^T\!\!  + \!AP \!+\!  K^TB^T\!\!  +\!  BK \!\!&  B_1& \!\!  PC_1 \!+\!  K^T\!D_{12}^T\!\! 
    \\ B_1^T \!&\!  -\gamma
    I_{w} \!&\!  D_{11}^T\!\!  \\C_1P \!+\!  D_{12}K \!&\!  D_{11} \!&\! 
    -\gamma I_z\!\!\end{bmatrix} \!\prec\!  0,\\
    \end{array}$$ to obtain $K^{*}$, $\gamma^{*}$ and $P^{*}$. Compute $F^{*} := K^{*}(P^{*})^{-1}C^{+}$, where $C^{+}$ is a pseudo-inverse of $C$, and resolve problem with $F^0 := F^{*}$. If problem has a solution $X^0$ and $\gamma^0$ then terminate Phase 1.  Set $x^0 := (F^0, X^0,\gamma^0)$. Otherwise, perform Step 4. 

-   *Step 4.*  Apply the method in Subsection \[subsec:SOF\_design\] to solve the following BMI feasibility problem: $$\begin{aligned}
    \label{eq:BMI_Hinf_phase1}
    \textrm{Find } ~\textrm{$F$ and $P\succ 0$ such that}: ~~ (A+BFC)^TP + P(A+BFC) \prec 0.\end{aligned}$$ If this problem has a solution $F^0$ then go back to Step 2.  Otherwise, declare that no strictly feasible point for is found. [  ]{}[$\square$]{}

As in the $\mathcal{H}_2$ problem, Phase 1 of the $\mathcal{H}_{\infty}$ is also terminated after finitely many iterations.  In this subsection, we also test this algorithm for several problems in COMP$\textrm{l}_{\textrm{e}}$ib using the same parameters and the stopping criterion as in the previous subsection.  The computational results are shown in Table \[tb:H\_inf\_problems\]. 

  ---------------- -------- -------- -------- -------- -------- ----------------- ----------------------- ------------------------ --------- -------------- --
  `Name`            $n_x$    $n_y$    $n_u$    $n_z$    $n_w$               HIFOO              PENBMI       $\mathcal{H}_\infty$      `iter`      `time`    
  [AC1]{}           [5]{}    [3]{}    [3]{}    [2]{}    [3]{}      [**0.0000**]{}               [**-**]{}           [**0.0177**]{}    [93]{}     [28.050]{} 
  [AC2]{}           [5]{}    [3]{}    [3]{}    [5]{}    [3]{}      [**0.1115**]{}               [**-**]{}           [**0.1140**]{}    [99]{}     [32.540]{} 
  [AC3]{}           [5]{}    [4]{}    [2]{}    [5]{}    [5]{}      [**4.7021**]{}               [**-**]{}           [**3.4859**]{}   [210]{}     [76.170]{} 
  [AC4]{}           [4]{}    [2]{}    [1]{}    [2]{}    [2]{}      [**0.9355**]{}               [**-**]{}          [**69.9900**]{}     [2]{}      [2.620]{} 
  [AC6]{}           [7]{}    [4]{}    [2]{}    [7]{}    [7]{}      [**4.1140**]{}               [**-**]{}           [**4.1954**]{}   [167]{}    [138.570]{} 
  [AC7]{}           [9]{}    [2]{}    [1]{}    [1]{}    [4]{}      [**0.0651**]{}          [**0.3810**]{}           [**0.0548**]{}   [300]{}    [278.330]{} 
  [AC8]{}           [9]{}    [5]{}    [1]{}    [2]{}    [10]{}     [**2.0050**]{}               [**-**]{}           [**3.0520**]{}   [247]{}    [298.070]{} 
  [AC9$^{b}$]{}     [10]{}   [5]{}    [4]{}    [2]{}    [10]{}     [**1.0048**]{}               [**-**]{}           [**0.9237**]{}   [300]{}    [470.910]{} 
  [AC11$^{b}$]{}    [5]{}    [4]{}    [2]{}    [5]{}    [5]{}      [**3.5603**]{}               [**-**]{}           [**3.0104**]{}    [68]{}     [60.260]{} 
  [AC12]{}          [4]{}    [4]{}    [3]{}    [1]{}    [3]{}      [**0.3160**]{}               [**-**]{}           [**2.3025**]{}   [300]{}    [181.870]{} 
  [AC15]{}          [4]{}    [3]{}    [2]{}    [6]{}    [4]{}     [**15.2074**]{}        [**427.4106**]{}          [**15.1995**]{}   [105]{}     [36.700]{} 
  [AC16]{}          [4]{}    [4]{}    [2]{}    [6]{}    [4]{}     [**15.4969**]{}               [**-**]{}          [**14.9881**]{}   [186]{}     [68.820]{} 
  [AC17]{}          [4]{}    [2]{}    [1]{}    [4]{}    [4]{}      [**6.6124**]{}               [**-**]{}           [**6.6373**]{}   [129]{}     [42.400]{} 
  [HE1$^{b}$]{}     [4]{}    [1]{}    [2]{}    [2]{}    [2]{}      [**0.1540**]{}          [**1.5258**]{}           [**0.1807**]{}   [300]{}    [143.520]{} 
  [HE2]{}           [4]{}    [2]{}    [2]{}    [4]{}    [4]{}      [**4.4931**]{}               [**-**]{}           [**6.7846**]{}   [177]{}     [67.470]{} 
  [HE3]{}           [8]{}    [6]{}    [4]{}    [10]{}   [1]{}      [**0.8545**]{}          [**1.6843**]{}           [**0.9243**]{}   [105]{}     [95.000]{} 
  [HE4$^{b}$]{}     [8]{}    [6]{}    [4]{}    [12]{}   [8]{}     [**23.3448**]{}               [**-**]{}          [**22.8713**]{}   [252]{}    [325.580]{} 
  [HE5$^{b}$]{}     [8]{}    [2]{}    [4]{}    [4]{}    [3]{}      [**8.8952**]{}               [**-**]{}          [**37.3906**]{}   [300]{}    [430.820]{} 
  [REA1]{}          [4]{}    [3]{}    [2]{}    [4]{}    [4]{}      [**0.8975**]{}               [**-**]{}           [**0.8815**]{}    [96]{}     [34.430]{} 
  [REA2$^{b}$]{}    [4]{}    [2]{}    [2]{}    [4]{}    [4]{}      [**1.1881**]{}               [**-**]{}           [**1.4188**]{}   [300]{}    [118.320]{} 
  [REA3]{}          [12]{}   [3]{}    [1]{}    [12]{}   [12]{}    [**74.2513**]{}         [**74.4460**]{}          [**74.5478**]{}     [2]{}    [234.800]{} 
  [DIS1]{}          [8]{}    [4]{}    [4]{}    [8]{}    [1]{}      [**4.1716**]{}               [**-**]{}           [**4.1943**]{}    [93]{}     [66.130]{} 
  [DIS2]{}          [3]{}    [2]{}    [2]{}    [3]{}    [3]{}      [**1.0548**]{}          [**1.7423**]{}           [**1.1546**]{}    [54]{}     [17.120]{} 
  [DIS3]{}          [6]{}    [4]{}    [4]{}    [6]{}    [6]{}      [**1.0816**]{}               [**-**]{}           [**1.1382**]{}   [285]{}    [195.960]{} 
  [DIS4]{}          [6]{}    [6]{}    [4]{}    [6]{}    [6]{}      [**0.7465**]{}               [**-**]{}           [**0.7498**]{}   [126]{}     [93.220]{} 
  [TG1$^{b}$]{}     [10]{}   [2]{}    [2]{}    [10]{}   [10]{}    [**12.8462**]{}               [**-**]{}          [**12.9336**]{}    [45]{}     [84.380]{} 
  [AGS]{}           [12]{}   [2]{}    [2]{}    [12]{}   [12]{}     [**8.1732**]{}        [**188.0315**]{}           [**8.1732**]{}    [24]{}     [55.290]{} 
  [WEC2]{}          [10]{}   [4]{}    [3]{}    [10]{}   [10]{}     [**4.2726**]{}         [**32.9935**]{}           [**6.6082**]{}   [300]{}    [476.010]{} 
  [WEC3]{}          [10]{}   [4]{}    [3]{}    [10]{}   [10]{}     [**4.4497**]{}        [**200.1467**]{}           [**6.8402**]{}   [300]{}    [425.330]{} 
  [BDT1]{}          [11]{}   [3]{}    [3]{}    [6]{}    [1]{}      [**0.2664**]{}               [**-**]{}           [**0.8562**]{}    [29]{}     [40.910]{} 
  [MFP]{}           [4]{}    [2]{}    [3]{}    [4]{}    [4]{}     [**31.5899**]{}               [**-**]{}          [**31.6079**]{}   [171]{}     [57.430]{} 
  [IH]{}            [21]{}   [10]{}   [11]{}   [11]{}   [21]{}     [**1.9797**]{}               [**-**]{}           [**1.1858**]{}   [114]{}   [1206.340]{} 
  [CSE1]{}          [20]{}   [10]{}   [2]{}    [12]{}   [1]{}      [**0.0201**]{}               [**-**]{}           [**0.0220**]{}     [3]{}     [20.250]{} 
  [PSM]{}           [7]{}    [3]{}    [2]{}    [5]{}    [2]{}      [**0.9202**]{}               [**-**]{}           [**0.9227**]{}    [87]{}     [67.470]{} 
  [EB1]{}           [10]{}   [1]{}    [1]{}    [2]{}    [2]{}      [**3.1225**]{}         [**39.9526**]{}           [**2.0276**]{}   [300]{}    [295.420]{} 
  [EB2]{}           [10]{}   [1]{}    [1]{}    [2]{}    [2]{}      [**2.0201**]{}         [**39.9547**]{}           [**0.8148**]{}    [84]{}     [73.770]{} 
  [EB3]{}           [10]{}   [1]{}    [1]{}    [2]{}    [2]{}      [**2.0575**]{}    [**3995311.0743**]{}           [**0.8153**]{}    [84]{}     [75.820]{} 
  [NN1]{}           [3]{}    [2]{}    [1]{}    [3]{}    [3]{}     [**13.9782**]{}         [**14.6882**]{}          [**18.4813**]{}   [300]{}    [144.940]{} 
  [NN2]{}           [2]{}    [1]{}    [1]{}    [2]{}    [2]{}      [**2.2216**]{}               [**-**]{}           [**2.2216**]{}     [9]{}      [4.060]{} 
  [NN4]{}           [4]{}    [3]{}    [2]{}    [4]{}    [4]{}      [**1.3627**]{}               [**-**]{}           [**1.3802**]{}   [156]{}     [51.480]{} 
  [NN8]{}           [3]{}    [2]{}    [2]{}    [3]{}    [3]{}      [**2.8871**]{}   [**78281181.1490**]{}           [**2.9345**]{}   [180]{}     [51.830]{} 
  [NN9$^{b}$]{}     [5]{}    [2]{}    [3]{}    [4]{}    [2]{}     [**28.9083**]{}               [**-**]{}          [**32.1222**]{}   [300]{}    [129.920]{} 
  [NN11$^{b}$]{}    [16]{}   [5]{}    [3]{}    [3]{}    [3]{}      [**0.1037**]{}               [**-**]{}           [**0.1566**]{}     [9]{}    [366.890]{} 
  [NN15]{}          [3]{}    [2]{}    [2]{}    [4]{}    [1]{}      [**0.1039**]{}               [**-**]{}           [**0.1194**]{}     [6]{}      [3.740]{} 
  [NN16]{}          [8]{}    [4]{}    [4]{}    [4]{}    [8]{}      [**0.9557**]{}               [**-**]{}           [**0.9656**]{}    [48]{}     [37.950]{} 
  [NN17]{}          [3]{}    [1]{}    [2]{}    [2]{}    [1]{}     [**11.2182**]{}               [**-**]{}          [**11.2381**]{}   [117]{}     [32.160]{} 
  ---------------- -------- -------- -------- -------- -------- ----------------- ----------------------- ------------------------ --------- -------------- --

  : $\mathcal{H}_{\infty}$ synthesis benchmarks on COMP$\textrm{l}_{\textrm{e}}$ib plants[]{data-label="tb:H_inf_problems"}

-0.4cm

The numerical results computed by HIFOO and PENBMI are also included in Table \[tb:H\_inf\_problems\].  Here, the notation is the same as in Table \[tb:H2\_prob\], except that $\mathcal{H}_{\infty}$ denotes the $\mathcal{H}_{\infty}$-norm of the closed-loop system for the static output feedback controller.  We can see from Table \[tb:H\_inf\_problems\] that the optimal values reported by Algorithm \[alg:A1\] and HIFOO are almost similar for many problems whereas in general PENBMI has difficulties in finding a feasible solution. 

$\mathcal{H}_2/\mathcal{H}_{\infty}$ optimization: BMI formulation {#subsec:H2Hinf_prob}
------------------------------------------------------------------

Motivated from the $\mathcal{H}_2$ and $\mathcal{H}_{\infty}$ optimization problems, in this subsection we consider the mixed $\mathcal{H}_2/\mathcal{H}_{\infty}$ synthesis problem.  Let us assume that $D_{11} = 0$, $D_{21} = 0$ and the performance output $z$ is divided in two components, $z_1$ and $z_2$. Then the linear system becomes: $$\label{eq:mixed_LTI_z1z2}
\left\{\begin{array}{cl}
&\dot{x} = Ax + B_1w + Bu,\\
&z_1 = C_1^{z_1}x + D_{12}^{z_1}u,\\
&z_2 = C_1^{z_2}x + D_{12}^{z_2}u,\\
&y   = Cx. 
\end{array}\right.$$ The mixed $\mathcal{H}_2/\mathcal{H}_{\infty}$ control problem is to find a static output feedback gain $F$ such that, for $u = Fy$, the $\mathcal{H}_2$-norm of the closed loop from $w$ to $z_2$ is minimized, while the $\mathcal{H}_{\infty}$-norm from $w$ to $z_1$ is less than some imposed level $\gamma$ [@Boyd1994; @Leibfritz2004; @Ostertag2008]. 

This problem leads to the following optimization problem with BMI constraints [@Ostertag2008]: $$\label{eq:H2Hinf_prob}
\begin{array}{cl}
\displaystyle\min_{F, P_1,P_2,Z} \!\!\!\!\!\!\!\! &  \textrm{trace}(Z)\\
\textrm{s.t.}\!\!\!\!\!\!\!\!   &\begin{bmatrix}A_F^TP_1 + P_1A_F + (C_{F}^{z_1})^TC_{F}^{z_1} & P_1B_1\\ B_1^TP_1 & -\gamma^2I\end{bmatrix}\prec 0,\\
&\begin{bmatrix}A_F^TP_2 + P_2A_F & P_2B_1\\B_1^TP_2 & -I\end{bmatrix}\prec 0, ~~
\begin{bmatrix}P_2 & (C_{F}^{z_2})^T\\C_{F}^{z_2} & Z \end{bmatrix}\succ 0, ~ P_1 \succ 0, ~P_2\succ 0,
\end{array}$$ where $A_F := A + BFC$, $C_{F}^{z_1} := C_1^{z_1} + D_{12}^{z_1}FC$ and $C_{F}^{z_2} := C_1^{z_2} + D_{12}^{z_2}FC$. Note that if $C = I_{n_x}$, the identity matrix, then this problem becomes a mixed $\mathcal{H}_2/\mathcal{H}_{\infty}$ of static state feedback design problem considered in [@Ostertag2008].  In this subsection, we test Algorithm \[alg:A1\] for the static state feedback and output feedback cases. 

0.2cm **Case 1.** * The static state feedback case ($C = I_{n_x}$)*.  First, we apply the method in [@Ostertag2008] to find an initial point via solving two optimization problems with LMI constraints.  Then, we use the same approach as in the previous subsections to transform problem into an optimization problem with psd-convex-concave matrix inequality constraints.  Finally, Algorithm \[alg:A1\] is implemented to solve the resulting problem.  For convenience of implementation, we introduce a slack variable $\eta$ and then replace the objective function in by $f(x) = \eta^2$ with an additional constraint $\textrm{trace}(Z) \leq \eta^2$.

In the first case, we test Algorithm \[alg:A1\] with three problems.  The first problem was also considered in [@Hassibi1999] with $$\begin{aligned}
&A   = {\footnotesize\begin{bmatrix}-1.40 & -0.49 & -1.93 \\ -1.73 & -1.69 & -1.25 \\ 0.99 & 2.08 & -2.49\end{bmatrix}}, ~
B_1  = {\footnotesize\begin{bmatrix}-0.16 & -1.29 \\ 0.81 & 0.96 \\ 0.41 & 0.65 \end{bmatrix}}, ~ B = {\footnotesize\begin{bmatrix} 0.25 \\ 0.41 \\ 0.65
\end{bmatrix}},\\
&C_1^{z_1} = {\footnotesize\begin{bmatrix} -0.41 & 0.44 & 0.68 \end{bmatrix}}, ~ C_1^{z_2} = {\footnotesize\begin{bmatrix} -1.77 & 0.50 & -0.40
\end{bmatrix}}, ~ D_{12}^{z_1} = D_{12}^{z_2} = 1,
~\textrm{and}~ \gamma = 2.\end{aligned}$$ If the tolerance $\varepsilon = 10^{-3}$ is chosen then Algorithm \[alg:A1\] converges after $17$ iterations and reports the value $\eta = 0.7489$ with $F = {\footnotesize \begin{bmatrix}
1.9485  & 0.3990  & -0.2119
\end{bmatrix}.}$ This result is similar to the one shown in [@Ostertag2008].  If we regularize the subproblem with $\rho = 0.5\times 10^{-3}$ and $Q=I_{PF}$ then the number of iterations is reduced to $10$ iterations. 

The second problem is `DIS4` in COMP$\textrm{l}_{\textrm{e}}$ib [@Leibfritz2003].  In this problem, we set $C_1^{z_1} = C_1^{z_2}$ and $D_{12}^{z_1}
= D_{12}^{z_2}$ as in [@Ostertag2008].  Algorithm \[alg:A1\] converges after $24$ iterations with the same tolerance $\varepsilon = 10^{-3}$. It reports $\eta = 1.6925$ and $\gamma = 1.1996$ with $$F = {\footnotesize\begin{bmatrix}
  -0.8663 \!&\!   -0.6504 \!&\!   -1.1115 \!&\!   -0.1951 \!&\!   -0.6099 \!&\!    0.2065 \\
   0.1591 \!&\!   -0.4941 \!&\!   -0.6322 \!&\!   -0.5409 \!&\!   -1.2895 \!&\!    0.2774 \\
  -0.7017 \!&\!   -0.0785 \!&\!    0.6121 \!&\!   -0.8919 \!&\!    0.2518 \!&\!   -0.2354 \\
  -0.0522 \!&\!   -0.5556 \!&\!   -0.5838 \!&\!    0.4497 \!&\!   -1.4279 \!&\!   -0.6677
\end{bmatrix}}.$$ If we regularize the subproblem with $\rho = 0.5\times 10^{-3}$ and $Q= I_{PF}$ then the number of iterations is $18$ iterations. 

The third problem is `AC16` in COMP$\textrm{l}_{\textrm{e}}$ib [@Leibfritz2003].  In this example we also choose $C_1^{z_1} = C_1^{z_2}$ and $D_{12}^{z_1} = D_{12}^{z_2}$ as in the previous problem.  As mentioned in [@Ostertag2008], if we choose a starting value $\gamma_0 = 100$, then the LMI problem can not be solved by the SDP solvers (e.g., Sedumi, SDPT3) due to numerical problems.  Thus, we rescale the LMI constraints using the same trick as in [@Ostertag2008].  After doing this, Algorithm \[alg:A1\] converges after $298$ iterations with the same tolerance $\varepsilon = 10^{-3}$. The value of $\eta$ reported in this case is $\eta = 12.3131$ and $\gamma = 20.1433$ with $$F = {\footnotesize\begin{bmatrix}
   -1.8533 &  0.1737 &  0.6980 &  6.4208 \\
    4.2672 & -0.9668 & -1.5952 & -2.9240
\end{bmatrix}}.$$ The results obtained by Algorithm \[alg:A1\] for solving problems `DIS4` and `AC16` in this paper confirm the results reported in [@Ostertag2008]. 

**Case 2.** * The static output feedback case.*  Similarly to the previous subsections, we first propose a technique to determine a starting point for Algorithm \[alg:A1\].  We described this phase algorithmically as follows. 

**Phase 1.** (* Determine a starting point $x^0$*). 

*Step 1.*  If $\alpha_0(A^T+A) < 0$ then set $F^0 = 0$. Otherwise, go to Step 3. 

*Step 2.*  Solve the following linear SDP problem: $$\label{eq:H2Hinf_LMI_prob}
\begin{array}{cl}
\displaystyle\min_{P_1,P_2, Z} \!\!\!\!\!\!\!\! &  \textrm{trace}(Z)\\
\textrm{s.t.}\!\!\!\!\!\!\!\!   &\begin{bmatrix}A_{F^0}^TP_1 + P_1A_{F^0} + (C_{F^0}^{z_1})^TC_{F^0}^{z_1} & P_1B_1\\ B_1^TP_1 & -\gamma^2I\end{bmatrix}\prec
0,\\
&\begin{bmatrix}A_{F^0}^TP_2 + P_2A_{F^0} & P_2B_1\\B_1^TP_2 & -I\end{bmatrix}\prec 0, ~\begin{bmatrix}P_2 & (C_{F^0}^{z_2})^T\\C_{F^0}^{z_2} & Z
\end{bmatrix}\succ 0, ~~ P_1 \succ 0, ~P_2\succ 0,
\end{array}$$ where $A_{F^0} = A + BF^0C$, $C_{F^0}^{z_1} = C_1^{z_1}+D_{12}^{z_1}F^0C$ and $C_{F^0}^{z_2} = C_1^{z_2}+D_{12}^{z_2}F^0C$. If this problem has an optimal solution $P_1^0, P_2^0$ and $Z^0$ then terminate Phase 1.  Set $x^0 := (F^0, P_1^0, P_2^0, Z^0)$ for a starting point of Algorithm \[alg:A1\].  Otherwise go to Step 3. 

*Step 3.*  Solve the following LMI feasibility problem: $$\label{eq:H2Hinf_LMI_phase1}
\begin{array}{cl}
&\textrm{Find $Q\succ 0$, $W$ and $Z$ such that:}\\
&~~\begin{bmatrix}
AQ \!+\!  QA^T\!\!  +\!  BW \!+\!  W^T\!B^T \!\!&\!\!  B_1 \!\!&\!\! ( C_1 \!+\!  D_{12}W)^T\!\! \\\!\!  B_1^T \!\!&\!\!  -I_w \!\!&\!\!  O \!\!\\ \!\!C_1 \!+\!  D_{12}W
\!\!&\!\!  O &\!\!  -\gamma^2I_z\!\! 
\end{bmatrix} \!\prec\!  0,\\
&~~\begin{bmatrix}AQ + QA^T + BW + W^TB^T & B_1 \\ B_1^T & -I_w\end{bmatrix}\prec 0,~
\begin{bmatrix}Q & (C_1Q + D_{12}W)^T \\C_1Q + D_{12}W & Z \end{bmatrix}\succ 0,\\
\end{array}$$ to obtain a solution $Q^{*}$, $W^{*}$ and $Z^{*}$. Set $F^{*} := W^{*}(Q^{*})^{-1}C^{+}$, where $C^{+}$ is the pseudo-inverse of $C$.  Solve again problem with $F^0 := F^{*}$. If problem has solution then terminate Phase 1.  Otherwise, perform Step 4. 

*Step 4.*  Solve the following optimization with BMI constraints: $$\label{eq:H2Hinf_BMI_prob}
\begin{array}{cl}
\displaystyle\max_{\beta, F, P_1,P_2} \!\!\!\!\!\!\!\! &\!\!\!  \beta\\
\!\!\!\!\!\!\!\!\textrm{s.t.}\!\!\!\!\!\!\!\! &  \!\!\!P_1 \succ 0, ~P_2\succ 0,\\
\!\!\!\!\!\!\!\! & \!\!  A_F^T\!P_1 \!+\!  P_1A_F^T \!  +\! ( C_{F}^{z_1})^T\!C_{F}^{z_1} + \frac{1}{\gamma^2}P_1B_1B_1^TP_1\preceq -2\beta P_1,\\
\!\!\!\!\!\!\!\! & \!\!  A_F^TP_2 + P_2A_F + P_2B_1B_1^TP_2 \preceq -2\beta P_2
\end{array}$$ to obtain an optimal solution $F^{*}$ corresponding to the optimal value $\beta^{*}$. If $\beta^{*} > 0$ then set $F^0 := F^{*}$ and go back to Step 2 to determine $P_1^0$, $P_2^0$ and $Z^0$. Otherwise, declare that no strictly feasible point of problem is found. [ $\square$]{}

Since at Step 4 of Phase 1, it requires to solve an optimization problem with two BMI constrains.  This task is usually expensive.  In our implementation, we only terminate this step after find a strictly feasible point with a feasible gap $0.1$. If matrix $C$ is invertible then the matrix $F^{*}$ at Step 3 is $F^{*} = W^{*}(Q^{*})^{-1}C^{-1}$. Hence, we can ignore Step 4 of Phase 1. 

To avoid the numerical problem in Step 3, we can reformulate problem equivalently to the following one:$$\label{eq:H2Hinf_LMI2_phase1}
\begin{array}{cl}
&\textrm{Find $Q\succ 0$, $W$ and $Z$ such that:}\\
&~\begin{bmatrix}
AQ \!+\!  QA^T\!\!  +\!  BW \!+\!  W^T\!B^T \!\!\!\!&\!\!\!\!  B_1 \!\!&\!\! ( C_1 \!+\!  D_{12}W)^T\!\! \\\!\!  B_1^T \!\!\!\!&\!\!\!\!  -\gamma I_w \!\!&\!\!  O \!\!\\
\!\!C_1 \!+\!  D_{12}W \!\!&\!\!  O &\!\! 
-\gamma I_z\!\! 
\end{bmatrix} \!\prec\!  0,~
\begin{bmatrix}Q \!\!\! &\!\!\! ( C_1Q + D_{12}W)^T \\C_1Q + D_{12}W \!\!\!&\!\!\!  Z \end{bmatrix} \!  \succ\!  0.\\
\end{array}$$ We test the algorithm described above for several problems in COMP$\textrm{l}_{\textrm{e}}$ib with the level values $\gamma = 4$ and $\gamma = 10$. In these examples, we assume that the output signals $z_1\equiv z_2$. Thus we have $C_1^{z_1} = C_1^{z_2} = C_1$ and $D_{12}^{z_1} = D_{12}^{z_2} = D_{12}$. The parameters and the stopping criterion of the algorithm are chosen as in Subsection \[subsec:H2\_prob\].  The computational results are reported in Table \[tb:H2Hinf\_problems\] with $\gamma = 4$ and $\gamma = 10$.

  --------------- -------------------------------------------- --------- ------------ ------------------------------------------- --------- ------------ --
  `Name`            $\mathcal{H}_2/\mathcal{H}_{\infty}$                       `time`   $\mathcal{H}_2/\mathcal{H}_{\infty}$         `iter`       `time` 
  [[ ]{}AC1]{}                         [**0.0585 / 0.0990**]{}     [3]{}     [4.22]{}                     [**0.0585 / 0.0990**]{}     [3]{}     [4.27]{} 
  [[ ]{}AC2]{}                         [**0.1067 / 0.1723**]{}     [6]{}     [7.31]{}                     [**0.1070 / 0.1727**]{}     [3]{}     [7.15]{} 
  [[ ]{}AC3]{}                         [**5.2770 / 3.9999**]{}    [51]{}   [281.53]{}                     [**4.5713 / 5.1298**]{}    [18]{}    [19.18]{} 
  [[ ]{}AC6]{}                                   [**- / -**]{}    [- ]{}        [-]{}                     [**4.0297 / 4.8753**]{}   [283]{}   [330.64]{} 
  [[ ]{}AC7]{}                         [**0.0415 / 0.0961**]{}     [1]{}     [3.39]{}                     [**0.0420 / 0.1286**]{}     [2]{}     [3.91]{} 
  [[ ]{}AC8]{}                         [**1.2784 / 2.2288**]{}    [43]{}    [60.78]{}                     [**1.3020 / 2.5719**]{}    [23]{}    [31.59]{} 
  [[ ]{}AC11]{}                        [**4.0704 / 4.0000**]{}    [76]{}   [175.75]{}                     [**4.0021 / 5.1949**]{}   [117]{}   [122.86]{} 
  [[ ]{}AC12]{}                        [**0.0924 / 0.3486**]{}    [18]{}    [73.46]{}                     [**1.4454 / 1.6444**]{}   [300]{}   [234.13]{} 
  [[ ]{}AC17]{}                                [**- /  -** ]{}     [-]{}        [-]{}                     [**4.1228 / 6.6472**]{}     [2]{}   [11.620]{} 
  [[ ]{}HE1]{}                         [**0.1123 / 0.2257**]{}     [2]{}   [131.18]{}                     [**0.0973 / 0.2080**]{}     [1]{}    [30.97]{} 
  [[ ]{}HE2]{}                                 [**-  /  -**]{}     [-]{}        [-]{}                     [**4.7302 / 9.8931**]{}    [75]{}    [55.48]{} 
  [[ ]{}REA1]{}                        [**1.8214 / 1.4740**]{}    [30]{}    [25.64]{}                     [**1.8213 / 1.4730**]{}    [30]{}    [26.65]{} 
  [[ ]{}REA2]{}                        [**3.5014 / 3.5180**]{}    [42]{}    [22.09]{}                     [**3.5015 / 3.5209**]{}    [45]{}    [23.26]{} 
  [[ ]{}DIS1]{}                                [**-  /  -**]{}     [-]{}        [-]{}                     [**2.8505 / 4.7904**]{}    [15]{}    [30.51]{} 
  [[ ]{}DIS2]{}                        [**1.5079 / 1.8500**]{}    [18]{}     [7.92]{}                     [**1.5079 / 1.8520**]{}    [21]{}     [7.92]{} 
  [[ ]{}DIS3]{}                        [**2.0577 / 1.7934**]{}    [27]{}    [25.03]{}                     [**2.0577 / 1.7766**]{}    [30]{}    [24.54]{} 
  [[ ]{}DIS4]{}                        [**1.6926 / 1.1952**]{}    [21]{}    [18.62]{}                     [**1.6926 / 1.2009**]{}    [24]{}    [21.55]{} 
  [[ ]{}AGS]{}                               [**-   /   -**]{}     [-]{}        [-]{}                     [**7.0332 / 8.2035**]{}     [8]{}   [196.73]{} 
  [[ ]{}PSM]{}                         [**1.5115 / 0.9248**]{}   [177]{}   [160.41]{}                     [**1.5115 / 0.9248**]{}   [180]{}   [167.31]{} 
  [[ ]{}EB2]{}                         [**0.7765 / 1.0828**]{}     [7]{}     [9.70]{}                     [**0.7768 / 1.0867**]{}    [10]{}    [13.16]{} 
  [[ ]{}EB3]{}                         [**0.8406 / 0.9249**]{}     [1]{}     [3.21]{}                     [**0.8383 / 0.9418**]{}     [1]{}     [2.93]{} 
  [[ ]{}EB4]{}                         [**1.0147 / 1.0707**]{}     [6]{}    [59.55]{}                     [**0.9981 / 1.2146**]{}    [12]{}   [111.26]{} 
  [[ ]{}NN2]{}                         [**1.5651 / 2.4834**]{}    [12]{}     [5.37]{}                     [**1.5651 / 2.4876**]{}    [12]{}     [5.49]{} 
  [[ ]{}NN4]{}                         [**1.8778 / 2.0501**]{}   [202]{}   [154.49]{}                     [**1.8779 / 2.0519**]{}   [213]{}   [161.00]{} 
  [[ ]{}NN8]{}                         [**2.3609 / 3.9999**]{}    [21]{}    [15.71]{}                     [**2.3376 / 4.6514**]{}    [15]{}     [6.57]{} 
  [[ ]{}NN15]{}                        [**0.0820 / 0.1010**]{}    [42]{}    [18.75]{}                       [**0.0771/0.1012**]{}    [24]{}    [10.47]{} 
  [[ ]{}NN16]{}                        [**0.3187 / 0.9574**]{}    [90]{}    [96.44]{}                     [**0.3319 / 0.9572**]{}   [258]{}   [303.87]{} 
  --------------- -------------------------------------------- --------- ------------ ------------------------------------------- --------- ------------ --

  : $\mathcal{H}_2/\mathcal{H}_{\infty}$ synthesis benchmarks on COMP$\textrm{l}_{\textrm{e}}$ib plants[]{data-label="tb:H2Hinf_problems"}

-0.4cm

Here, $\mathcal{H}_2/\mathcal{H}_{\infty}$ are the $\mathcal{H}_2$ and $\mathcal{H}_{\infty}$ norms of the closed-loop systems for the static output feedback controller, respectively.  With $\gamma = 10$, the computational results show that Algorithm \[alg:A1\] satisfies the condition ${\|P_{\infty}(s)\|}_{\infty} \leq \gamma = 10$ for all the test problems.  While, with $\gamma = 4$, there are $5$ problems reported infeasible, which are denoted by “-”.  The $\mathcal{H}_{\infty}$-constraint of three problems: AC3, AC11 and NN8 is active with respect to $\gamma = 4$.

Concluding remarks
==================

We have proposed a new algorithm for solving many classes of optimization problems involving BMI constraints arising in static feedback controller design.  The convergence of the algorithm has been proved under standard assumptions.  Then, we have applied our method to design static feedback controllers for various problems in robust control design.  The algorithm is easy to implement using the current SDP software tools.  The numerical results are also reported for the benchmark collection in COMP$\textrm{l}_{\textrm{e}}$ib.  Note, however, that our method depends crucially on the psd-convex-concave decomposition of the BMI constraints.  In practice, it is important to look at the specific structure of the problems and find an appropriate psd-convex-concave decomposition for Algorithm \[alg:A1\].  The method proposed can be extended to general nonlinear semidefinite programming, where the psd-convex-concave decomposition of the nonconvex mappings are available.  From a control design point of view, the application to more general reduced order controller synthesis problems and the extension towards linear parameter varying or time-varying systems are future research directions. 

Acknowledgments
===============

[Research supported by Research Council KUL: CoE EF/05/006 Optimization in Engineering(OPTEC), IOF-SCORES4CHEM, GOA/10/009 (MaNet), GOA/10/11, several PhD/postdoc and fellow grants; Flemish Government: FWO: PhD / postdoc grants, projects G.0452.04, G.0499.04, G.0211.05, G.0226.06, G.0321.06, G.0302.07, G.0320.08, G.0558.08, G.0557.08, G.0588.09, G.0377.09, G.0712.11, research communities (ICCoS, ANMMM, MLDM); IWT: PhD Grants, Belgian Federal Science Policy Office: IUAP P6/04; EU: ERNSI; FP7-HDMPC, FP7-EMBOCON, Contract Research: AMINAL.  Other: Helmholtz-viCERP, COMET-ACCM. ]{} 

Proof of Lemma \[le:descent\_dir\]. {# app:A1 .unnumbered}
-----------------------------------

For any matrices $A, B\in\mathcal{S}^p_{+}$, we have $A\circ B \geq 0$. From *Step 1* of Algorithm \[alg:A1\], we have $x^{k+1}$ is the solution of the convex subproblem and $\Lambda^{k+1}$ is the corresponding multiplier, under Assumption \[as:A3\], they must satisfy the following generalized Kuhn-Tucker condition: $$\label{eq:KT_cond}
\left\{\begin{array}{cl}
&\!\!  0\in \nabla f(x^{k\!+\!1}) \!+\!  \rho_kQ^T_kQ_k(x^{k\!+\!1} \!-\!  x^k) \!+\!  \left\{\sum_{i=1}^lD[G_i(x) \!-\!  H_i(x^k)\right.\\
&\!\!  \left.- DH_i(x^k)(x-x^k)]|_{x^{k+1}}\right\}^{*}\Lambda^{k+1}_i + N_{\Omega}(x^{k+1}),\\
&\!\!G_i(x^{k+1}) - H_i(x^k) - DH(x^k)(x^{k+1}-x^k) \preceq 0, ~~ \Lambda_i \succeq 0,\\
&\!\!  \left[\!G_i(x^{k\!+\!1}) \!-\!  H_i(x^k) \!-\!  DH(x^k)(x^{k\!+\!1}\!-\!x^k)\!\right]\circ \Lambda_i^{k\!+\!1} = 0. 
\end{array}\right.$$ Noting that $D\left[G_i(x) - H_i(x^k) - DH_i(x^k)(y-x^{k})\right]|_{x=x^{k+1}} = DG_i(x^{k+1}) - DH_i(x^k)$ for $i=1,\dots l$, it follows from the first line of and the convexity of $f$ that $$\begin{aligned}
\label{eq:term1}
&f(y) \!-\!  f(x^{k\!+\!1}) \!+\!  \left\{\sum_{i=1}^l[DG_i(x^{k\!+\!1})\!-\!DH_i(x^k)]^{*}\Lambda_i^{k+1}\right\}^T\!\!\!(y-x^{k+1})\notag\\
&\geq \left\{\nabla f(x^{k+1}) + \sum_{i=1}^l[DG_i(x^{k+1})-DH_i(x^k)]^{*}\Lambda_i^{k+1}\right\}^T\!\!\!(y-x^{k+1})\notag\\
&+\!  \frac{\rho_f}{2}{\|y \!-\!  x^{k\!+\!1}\|}_2^2 \!  \geq\!  \frac{\rho_f}{2}{\|y \!-\!  x^{k\!+\!1}\|}^2_2 \!+\! 
\rho_k(y\!-\!x^{k\!+\!1})^TQ_k^TQ_k(x^k\!-\!x^{k\!+\!1}),  ~~\forall y\in\Omega.\end{aligned}$$ On the other hand, we have $$\begin{aligned}
\label{eq:term2}
&\left\{ [DG_i(x^{k + 1})  - DH_i(x^k)]^{*}\Lambda_i^{k + 1}\right\}^T\!(y - x^{k + 1}) \nonumber\\
& =  \Lambda_i^{k + 1} \circ
[DG_i(x^{k + 1})(y  - x^{k + 1}) - DH_i(x^k)(y - x^{k + 1})].\end{aligned}$$ Since $G_i$ and $H_i$ are psd-convex, applying Lemma \[le:psd\_convex\_properties\] we have $$\begin{aligned}
& G_i(x^k) \!-\!  G_i(x^{k\!+\!1}) \!\succeq\!  DG_i(x^{k\!+\!1}\!)(x^k \!-\!  x^{k\!+\!1}), \nonumber\\
\textrm{and}~& H_i(x^{k\!+\!1}) \!-\!  H_i(x^k) \!\succeq\!  DH_i(x^k)(x^{k\!+\!1} \!-\!  x^k),
~i=1,\dots, l.\end{aligned}$$ Summing up these inequalities we obtain $$\begin{aligned}
G_i(x^k) \!-\!  H_i(x^k) \!-\! [ G_i(x^{k\!+\!1}) \!-\!  H_i(x^{k\!+\!1})]  \!\succeq\! [ DG_i(x^{k \!+\!  1})(x^k \!-\!  x^{k+1}) \!-\! 
DH_i(x^k)(x^k-x^{k+1})].\end{aligned}$$ Using the fact that $\Lambda^{k+1}_i\succeq 0$, this inequality implies that $$\begin{aligned}
\label{eq:term3}
& \Lambda^{k + 1}_i  \circ  \left\{ G_i(x^k) - H_i(x^k)  - [G_i(x^{k + 1} ) -  H_i(x^{k\!+\!1})] \right\} \nonumber\\
& \succeq  \Lambda^{k\!+\!1}_i  \circ [DG_i(x^{k\!+\!1})(x^k - x^{k\!+\!1})  -  DH_i(x^k)(x^k  -  x^{k\!+\!1})].\end{aligned}$$ Substituting $y=x^k$ into and then combining the consequence, , and the last line of to obtain $$\begin{aligned}
\label{eq:term4}
f(x^k) \!-\!  f(x^{k\!+\!1}) \!+\!  \sum_{i=1}^l\Lambda_i^{k\!+\!1}\circ[G_i(x^k)-H_i(x^k)] \geq \frac{\rho_f}{2}{\|x^{k \!+ \!  1}-x^k\|}_2^2 \!+\! 
\rho_k{\|Q_k(x^{k\!+\!  1} \!-\!  x^k)\|}_2^2.\end{aligned}$$ Now, since $x^k$ is the solution of the convex subproblem linearized at $x^{k-1}$. One has $G_i(x^k) - H_i(x^k) \preceq 0$. Moreover, since $\Lambda^{k+1}_i\succeq 0$, we have $\Lambda^{k+1}_i\circ\left[G_i(x^k) - H_i(x^k)\right] \leq 0$. Substituting this inequality into , we obtain $$f(x^k) - f(x^{k+1}) \geq \frac{\rho_f}{2}{\|x^k-x^{k+1}\|}_2^2 + \rho_k{\|Q_k(x^{k+1}-x^k)\|}_2^2.$$ This inequality is indeed which proves the item a).  If there exists at least one $i_0\in \left\{1,\dots, l\right\}$ such that $G_{i_0}(x^k) - H_{i_0}(x^k) \prec 0$ and $\Lambda^{k+1}_{i_0} \succ 0$ then $\Lambda_{i_0}^{k+1}\circ\left[G_{i_0}(x^k) -
H_{i_0}(x^k)\right] < 0$. Substituting this inequality into we conclude that $f(x^{k+1}) < f(x^k)$ which proves item b).  The last statement c) follows directly from the inequality . [ $\square$]{}

[99]{}

T. Alamo, J.M. Bravo, M.J. Redondo and E.F. Camacho, *A set-membership state estimation algorithm based on DC programming*, Automatica, vol.  44, no.  1, pp.  216–224, 2008. 

P. Apkarian, and H.D. Tuan, “Robust control via concave optimization: local and global algorithms”, *Proc of CDC*, 1998. 

A. Ben-Tal and A.K. Nemirovski, *Lectures on Modern Convex Optimization: Analysis, Algorithms, and Engineering Applications*, SIAM, Philadelphia, 2001. 

V.D. Blondel and J.N. Tsitsiklis, “NP-hardness of some linear control design problems” *SIAM J. on Control, Signals and Systems*, vol.  35, no.  21, pp.  18–27, 1997. 

S.P. Boyd, L.E. Ghaoui, E. Feron and V. Balakrishnan, *Linear matrix inequalities in system and control theory*, Vol.  15, SIAM studies in applied mathematics, SIAM, Philadelphia, 1994. 

S.P. Boyd and L. Vandenberghe, *Convex Optimization*, Cambridge University Press, 2004.  J.V. Burke, A.S. Lewis and M.L. Overton, “Two numerical methods for optimizing matrix stability”, *Linear Algebra and Its Applications*, vol.  351/352, pp.  117–145, 2002. 

J.V. Burke, D. Henrion, A.S. Lewis and M.L. Overton, “Stabilization via Nonsmooth, Nonconvex Optimization”, *IEEE Trans.  Automatic Control*, vol.  51, pp.  1760–1769, 2006. 

R. Correa and H. Ramirez, “A global algorithm for nonlinear semidefinite programming”, *SIAM J. Optim.*,  vol.  15, no.  1, pp.  303–318, 2004. 

B. Fares, D. Noll and P. Apkarian, “Robust control via sequential semidefinite programming”, *SIAM J. Control Optim.*,  vol.  40, no.  6, pp.  1791–1820, 2002. 

R.W. Freund, F. Jarre and C.H. Vogelbusch, “Nonlinear semidefinite programming: sensitivity, convergence, and an application in passive reduced-order modeling”, *Math.  Program.*,  vol.  109, Ser B, pp.  581–611, 2007. 

K.C. Goh, “Robust control synthesis via bilinear matrix inequalities”, Ph.D. thesis, University of Southern California, Los Angeles, CA, 1995. 

S. Gumussoy, D. Henrion, M. Millstone and M.L. Overton, “Multiobjective Robust Control with HIFOO $2.0$”, *Proceedings of the IFAC Symposium on Robust Control Design*, Haifa, 2009. 

A. Hassibi, J. How and S. Boyd, “A path following method for solving BMI problems in control”, *Proceedings of American Control Conference*, vol.  2, pp.  1385–1389, 1999. 

D. Henrion, J. Loefberg, M. Kocvara and M. Stingl, “Solving polynomial static output feedback problems with PENBMI”, *Proc.  joint IEEE Conf.  Decision Control and Europ.  Control Conf.*,  Sevilla, Spain, 2005. 

C.W.J. Hol and C.W. Scherer, *A sum-of-squares approach to fixed-order H-infinity synthesis*, in D. Henrion and A. Garulli (eds), *Positive Polynomials in Control*, Springer, pp.  45–71, 2005. 

R.A. Horn and C.R. Johnson, *Topics in Matrix Analysis*, Cambridge University Press, Cambridge, 1991. 

R. Horst and H. Tuy, *Global Optimization: Deterministic Approaches*, 3rd ed.  Springer-Verlag, Berlin, 1996. 

Software, vol.18, no.  3, pp.  317–333, 2003. 

M. Kočvara, F. Leibfritz, M. Stingl and D. Henrion, *A nonlinear SDP algorithm for static output feedback problems in $\mathrm{COMPL}_{\mathrm{e}}\mathrm{ib}$*, Proc.  IFAC World Congress, Prague, Czech Rep., 2005. 

F. Leibfritz and W. Lipinski, “Description of the benchmark examples in COMPleib 1.0”, Dept.  Math.,  Univ.  Trier, Trier, Germany, Tech.  Rep., 2003. 

F. Leibfritz, “COMPleib: Constraint matrix optimization problem library - a collection of test examples for nonlinear semidefinite programs, control system design and related problems”.  Dept.  Math.,  Univ.  Trier, Trier, Germany, Tech.  Rep., 2004. 

J. Löfberg., “ YALMIP : A Toolbox for Modeling and Optimization in MATLAB”. * Proceedings of the CACSD Conference*, Taipei, Taiwan, 2004. 

A.W. Marshall and I. Olkin, *Inequalities: Theory of Majorization and Its Applications*, Academic Press, New York, 1979. 

R. R. Meyer, “Sufficient conditions for the convergence of monotonic mathematical programming algorithms”, *Journal of Computer and System Sciences*, vol.  12, pp.  108–121, 1976. 

Mathematics, Courant Institute of Mathematical Y. Nesterov and A.K. Nemirovski, *Interior-Point Polynomial Methods in Convex Programming*, SIAM Series in Applied Mathematics, Philadelphia, PA, 1994. 

R. Orsi, U. Helmke and J.B. Moore, “A Newton-like method for solving rank constrained linear matrix inequalities”, Automatica, vol.  42, no.  11, pp.  1875–1882, 2006. 

E. Ostertag, “An improved path-following method for mixed $H_2/H_{\infty}$ controller design”, *IEEE Trans.  Automatic Control*, vol.  53, no.  8, pp.  1967–1971, 2008. 

D.T. Pham and H.A. Le Thi, “A DC optimization algorithms for solving the trust region subproblem”, *SIAM J. Optim.*,  vol.  8, pp.  476–507, 1998. 

A. Shapiro, “First and second order analysis of nonlinear semidefinite programs”, *Math.  Program.*  vol.  77, no.  1, pp.  301–320, 1997. 

B.K. Sriperumbudur and G.R.G. Lanckriet, “On the convergence of the concave-convex procedure”, *Neural Information Processing Systems*, NIPS, 2009. 

J.F. Sturm, “Using SeDuMi 1.02: A Matlab toolbox for optimization over symmetric cones”, *Optim.  Methods Software*, vol.  11-12, pp.  625–653, 1999. 

D. Sun, “The strong second order sufficient condition and constraint non-degeneracy in nonlinear semidefinite programming and their implications”, *Mathematics of Operation Research*, vol.  31, no.  4, pp.  761–776, 2006. 

R.H. Tütünkü, K.C. Toh and M.J. Todd, “SDTP3 - a Matlab software package for semidefinite-quadratic-linear programming version 3.0”, 2001. 

J.B. Thevenet, D. Noll and P. Apkarian, “Nonlinear spectral SDP method for BMI-constrained problems: applications to control design”, *Informatics in Control, Automation and Robotics*, vol.  1, pp.  61–72, 2006. 

L.N. Trefethen, “Pseudo-spectra of linear operators”, *SIAM Review*, vol.  39, pp.  383–406, 1997. 

F. Zhang, *The Schur complement and its applications*, Springer, New-York, 2005. 

[^1]: Department of Electrical Engineering (ESAT-SCD) and Optimization in Engineering Center (OPTEC), K.U. Leuven, Kasteelpark Arenberg 10, B-3001 Leuven, Belgium.  Email: {quoc.trandinh,moritz.diehl}@ esat.kuleuven.be

[^2]: Department of Computer Science, Katholieke Universiteit Leuven, Leuven, Belgium.  Email: {suat.gumussoy, wim.michiels}@cs.kuleuven.be

[^3]: Due to Schur’s complement form
