On approximating the nearest \Omega-stable matrix

TL;DR

This paper introduces a convex optimization approach to approximate matrices with eigenvalues in specific regions of the complex plane, relevant for system stability, using dissipative Hamiltonian matrices and linear matrix inequalities.

Contribution

It formulates the nearest -stable matrix problem within a convex framework using LMI and dissipative Hamiltonian parametrization, enabling efficient solutions.

Findings

Convex reformulation of the approximation problem.

Application of block coordinate descent for solutions.

Illustrative examples demonstrating effectiveness.

Abstract

In this paper, we consider the problem of approximating a given matrix with a matrix whose eigenvalues lie in some specific region \Omega, within the complex plane. More precisely, we consider three types of regions and their intersections: conic sectors, vertical strips and disks. We refer to this problem as the nearest \Omega-stable matrix problem. This includes as special cases the stable matrices for continuous and discrete time linear time-invariant systems. In order to achieve this goal, we parametrize this problem using dissipative Hamiltonian matrices and linear matrix inequalities. This leads to a reformulation of the problem with a convex feasible set. By applying a block coordinate descent method on this reformulation, we are able to compute solutions to the approximation problem, which is illustrated on some examples.

Tables

Table 1. Table 4.1: Comparison of the algorithms for randomly generated matrices A 𝐴 A . The table displays the average relative error in percent with the standard deviation obtained by each initialization and, in brackets, the number of times the algorithm found the best solution out of the 20 runs (up to 0.01%). The best results are highlighted in bold.

	Identity	LMI	True
$ϵ=0.01$	0.53 $\pm$ 0.70 (18)	0.00 $\pm$ 0.01 (20)	0.21 $\pm$ 0.12 (0)
$ϵ=0.05$	0.99 $\pm$ 0.53 (14)	0.25 $\pm$ 0.40 (19)	1.37 $\pm$ 0.74 (1)
$ϵ=0.10$	2.06 $\pm$ 1.64 (18)	2.09 $\pm$ 2.24 (12)	2.76 $\pm$ 1.63 (0)
$ϵ=0.20$	4.93 $\pm$ 2.01 (16)	4.79 $\pm$ 2.35 (11)	6.29 $\pm$ 2.63 (0)

Equations

$$\dot{x}(t)=Ax(t)+Bu(t),$$

$$\Omega_{C}(a,\theta):=\left\{x+iy\in{\mathbb{C}}\ \big{|}\ \sin(\theta)(x-a)<\cos(\theta)y<-\sin(\theta)(x-a),\,x\leq a\right\}.$$

$$\Omega_{V}(h,k):=\left\{x+iy\in{\mathbb{C}}\ \big{|}\ -k<x<-h\right\}.$$

$$\Omega_{D}(-q,r):=\left\{z\in{\mathbb{C}}\ \big{|}\ |z+q|<r\right\}.$$

$$\underset{X\in {\mathbb{S}}_{\Omega }^{n,n}}{\inf }{\Vert A-X\Vert }_F^2,$$

$$\mathrm{Re}(\lambda )=-\frac{v^{\ast }Rv}{v^{\ast }{Q}^{-1}v}\text{and}\mathrm{Im}(\lambda )=-i\frac{v^{\ast }Jv}{v^{\ast }{Q}^{-1}v}.$$

$$v^{\ast }(J-R)Q=\lambda v^{\ast }⟹v^{\ast }(J-R)v=\lambda v^{\ast }{Q}^{-1}v$$

$$v^{\ast }(-J-R)v=\bar{\lambda }{v}^{\ast }{Q}^{-1}v.$$

$$2v^{\ast }Jv$$

$$-2v^{\ast }Rv$$

$$\left[\begin{array}[]{cc}\alpha(R+aQ^{-1})&-\beta J\\ \beta J&\alpha(R+aQ^{-1})\end{array}\right]\succ 0.$$

$$\displaystyle-2\left[\begin{array}[]{cc}v^{*}&0\\ 0&v^{*}\end{array}\right]\left[\begin{array}[]{cc}\alpha\,(R+aQ^{-1})&-\beta\,J\\ \beta\,J&\alpha\,(R+aQ^{-1})\end{array}\right]\left[\begin{array}[]{cc}v&0\\ 0&v\end{array}\right]$$

$$\displaystyle\Longrightarrow\quad 2\left[\begin{array}[]{cc}-\alpha v^{*}(R+aQ^{-1})v&\beta v^{*}Jv\\ -\beta v^{*}Jv&-\alpha v^{*}(R+aQ^{-1})v\end{array}\right]$$

$$\displaystyle\Longrightarrow\quad 2\left[\begin{array}[]{cc}-\alpha v^{*}Rv&\beta v^{*}Jv\\ -\beta v^{*}Jv&-\alpha v^{*}Rv\end{array}\right]-\alpha\left[\begin{array}[]{cc}av^{*}Q^{-1}v&0\\ 0&av^{*}Q^{-1}v\end{array}\right]$$

$$\displaystyle v^{*}Q^{-1}v\left[\begin{array}[]{cc}\alpha\,(\lambda_{1}-a)&\beta\,i\lambda_{2}\\ -\beta\,i\lambda_{2}&\alpha\,(\lambda_{1}-a)\end{array}\right]\prec 0.$$

$$\left[\begin{array}[]{cc}\alpha(AX+XA^{T}-aX)&\beta(AX-XA^{T})\\ \beta(XA^{T}-AX)&\alpha(AX+XA^{T}-aX)\end{array}\right]\prec 0.$$

$$R=-\frac{AX+(AX{)}^T}{2},J=\frac{AX-(AX{)}^T}{2},\text{and}Q=X^{-1}.$$

$$\left[\begin{array}[]{cc}\alpha\,(R+aQ^{-1})&-\beta\,J\\ \beta\,J&\alpha\,(R+aQ^{-1})\end{array}\right]=-1/2\left[\begin{array}[]{cc}\alpha(AX+XA^{T}-aX)&\beta(AX-XA^{T})\\ \beta(XA^{T}-AX)&\alpha(AX+XA^{T}-aX)\end{array}\right]\succ 0.$$

$$k{Q}^{-1}\succ R\succ h{Q}^{-1}.$$

$$x^{\ast }Q(J-R)Qx=\lambda x^{\ast }Qx⟹\mathrm{Re}(\lambda )=-\frac{x^{\ast }QRQx}{x^{\ast }Qx}.$$

$$k>\frac{x^{\ast }QRQx}{x^{\ast }Qx}>h.$$

$$AX+X{A}^T+2hX\prec 0\text{and}AX+X{A}^T+2kX\succ 0.$$

$$Q=X^{-1},R=-\frac{AX+X{A}^T}{2},J=\frac{AX-X{A}^T}{2}.$$

$$\left[\begin{array}[]{cc}rQ^{-1}&-qQ^{-1}\\ -qQ^{-1}&rQ^{-1}\end{array}\right]\succ\left[\begin{array}[]{cc}0&J-R\\ (J-R)^{T}&0\end{array}\right].$$

$$\displaystyle\left[\begin{array}[]{cc}v^{*}&0\\ 0&v^{*}\end{array}\right]\left[\begin{array}[]{cc}rQ^{-1}&-qQ^{-1}\\ -qQ^{-1}&rQ^{-1}\end{array}\right]\left[\begin{array}[]{cc}v&0\\ 0&v\end{array}\right]\succ\left[\begin{array}[]{cc}v^{*}&0\\ 0&v^{*}\end{array}\right]\left[\begin{array}[]{cc}0&J-R\\ (J-R)^{T}&0\end{array}\right]\left[\begin{array}[]{cc}v&0\\ 0&v\end{array}\right].$$

$$\displaystyle\left[\begin{array}[]{cc}r&-q\\ -q&r\end{array}\right]v^{*}Q^{-1}v\succ\left[\begin{array}[]{cc}0&v^{*}(J-R)v\\ v^{*}(J-R)^{T}v&0\end{array}\right].$$

$$\displaystyle\left[\begin{array}[]{cc}r&-q\\ -q&r\end{array}\right]\succ\left[\begin{array}[]{cc}0&\frac{v^{*}(J-R)v}{v^{*}Q^{-1}v}\\ \frac{v^{*}(J-R)^{T}v}{v^{*}Q^{-1}v}&0\end{array}\right].$$

$$\displaystyle\left[\begin{array}[]{cc}r&-q\\ -q&r\end{array}\right]\succ\left[\begin{array}[]{cc}0&i\lambda_{2}+\lambda_{1}\\ -i\lambda_{2}+\lambda_{1}&0\end{array}\right]\quad\Longrightarrow\quad\left[\begin{array}[]{cc}r&-q-\lambda\\ -q-\overline{\lambda}&r\end{array}\right]\succ 0.$$

$$\left[\begin{array}[]{cc}-rX&qX+AX\\ qX+XA^{T}&-rX\end{array}\right]\prec 0.$$

$$Q=X^{-1},R=-\frac{AX+X{A}^T}{2},J=\frac{AX-X{A}^T}{2}.$$

$$\displaystyle 0\succ\left[\begin{array}[]{cc}-rX&qX+AX\\ qX+XA^{T}&-rX\end{array}\right]$$

Full text

Keywords: stability radius, linear time-invariant systems, $\Omega$-stability, convex optimization

1 Introduction

Let us consider the following linear time-invariant (LTI) systems of the form

[TABLE]

where $A\in\mathbb{R}^{n,n}$, $B\in\mathbb{R}^{n,m}$, $x$ is the state vector, $u$ is the input vector. Let also $\Omega$ be a subset of the complex plane. An LTI system is called $\Omega$-stable if all eigenvalues of $A$ lie inside $\Omega$. The most well-known example is when $\Omega$ is the open left-half of the complex plane which characterizes stability of continuous LTI systems. For discrete LTI systems of the form $x(l+1)=Ax(l)$ for $l\in\mathbb{N}$, stability requires $\Omega$ to be the unit disk. Another motivation for enforcing eigenvalues in specific regions of the complex plane is that the transient response is related to the location of its eigenvalues [12, 1]. For example, the step response of a system with eigenvalues $\lambda=-\tau\omega_{n}\pm i\omega_{d}$ is fully characterized in terms of the undamped natural frequency $\omega_{n}=|\lambda|$, the damping ratio $\tau$, and the damped natural frequency $\omega_{d}$ [12]. By constraining $\lambda$ to lie in a prescribed region, specific bounds can be put on these quantities to ensure a satisfactory transient response; see also [17] and the references therein. In this paper, we focus on three regions of the complex plane and their intersections, namely:

• •

Conic sector: the conic sector region of parameters $a,\,\theta\in{\mathbb{R}}$ with $0\leq\theta\leq\pi/2$, denoted by $\Omega_{C}(a,\theta)$, is defined as

[TABLE]

• •

Vertical strip: the vertical strip region of parameters $h<k$, denoted by $\Omega_{V}(h,k)$, is defined as

[TABLE]

Note that $h$ (resp. $k$) can possibly be equal to $-\infty$ (resp. $+\infty$) in which case $\Omega_{V}$ is a half space. In particular, $\Omega_{V}(0,+\infty)$ is the open left half of the complex plane, corresponding to stable matrices for continuous LTI systems.

• •

Disks centred on the real line: the disk centred at $(-q,0)$ with radius $r>0$, denoted by $\Omega_{D}(-q,r)$, is defined as

[TABLE]

In particular, $\Omega_{D}(0,1)$ is the unit disk, corresponding to stable matrices for discrete LTI systems.

The sets $\Omega$ considered in this paper can be either any of $\Omega_{C}$, $\Omega_{V}$, $\Omega_{D}$, or the intersection of such sets; see Figure 1.1 for an illustration. Note that $\Omega$ is symmetric with respect to the real line. Note also that it is useless to consider more than one set of the type $\Omega_{V}$ since the intersection of such sets can be simplified to a single set $\Omega_{V}$. However, it makes sense to consider the intersection of several sets of the types $\Omega_{C}$ and $\Omega_{D}$, that is, the intersection of several conic sectors and several disks.

In this paper, we consider the problem of computing the nearest $\Omega$-stable matrix to a given matrix. More precisely, for a given matrix $A\in{\mathbb{R}}^{n,n}$ we are interested in solving the following optimization problem

[TABLE]

where ${\|\cdot\|}_{F}$ denotes the Frobenius norm of a matrix and $\mathbb{S}_{\Omega}^{n,n}$ is the set of all $\Omega$-stable matrices of size $n\times n$. This problem is important for example in system identification where one needs to identify a stable system from observations. For example consider the $\Omega$-stability region as in Figure 1.1 which is the intersection of the regions $\Omega_{C}$, $\Omega_{V}$ and $\Omega_{D}$. A solution of the nearest $\Omega$-stable matrix problem for this region is useful to identify a stable system with a minimum decay rate $h$, a minimum damping ratio $\tau=\cos\theta$, and a maximum undamped natural frequency $\omega_{d}=r\sin\theta$. Such a system bounds the maximum overshoot, the frequency of oscillatory modes, the delay time, the rise time, and the setting time [3, 12]. Such nearness problems for LTI systems have recently been studied; see for example [15, 7, 13, 9] for the continuous LTI systems, and [15, 14, 10, 5] for discrete LTI systems. The converse of (1.1) is the distance to $\Omega$-instability, that is, for a given $\Omega$-stable matrix $A$, find the smallest perturbation $\Delta_{A}$ with respect to some norm such that the perturbed matrix $A+\Delta_{A}$ has at least one eigenvalue outside $\Omega$. It is the more constrained version of the widely studied distance to instability [2, 11].

Notation

Throughout the paper, $X^{T}$ and $\|X\|$ stand for the transpose and the spectral norm of a real square matrix $X$, respectively. We write $X\succ 0$ $(X\preceq 0)$ and $X\succeq 0$ $(X\preceq 0)$ if $X$ is symmetric and positive definite (negative definite) or positive semidefinite (negative semidefinite), respectively. By $I_{m}$ we denote the identity matrix of size $m\times m$.

Outline and contribution

In order to tackle (1.1), we extend the approach proposed in [7] that tackles the nearest stable matrix problem for continuous LTI systems. In [7], stable matrices for continuous LTI systems are parametrized using dissipative Hamiltonian (DH) matrices of the form $A=(J-R)Q$ where $J^{T}=-J$, $R\succeq 0$ and $Q\succ 0$. In fact, it turns out that a matrix is stable if and only if it is a DH matrix. In Section 2, we show how to impose the eigenvalues of a matrix of the form $A=(J-R)Q$ to lie in the sets $\Omega_{C}$, $\Omega_{V}$ and $\Omega_{D}$ using linear matrix inequalities (LMIs) derived in [3]. This allows us to reformulate in Section 3 the problem (1.1) into an equivalent optimization problem with a convex feasible set onto which it is easy to project. We then propose in Section 4 a block coordinate descent algorithm to tackle this problem and illustrate the effectiveness of this algorithm on several examples.

2 DH matrices and $\Omega$-stability

Let us formally define a DH matrix.

Definition 1**.**

A matrix $A\in{\mathbb{R}}^{n,n}$ is said to be a DH matrix if $A=(J-R)Q$ for some $J,R,Q\in{\mathbb{R}}^{n,n}$ such that $J^{T}=-J$, $R\succeq 0$ and $Q\succ 0$.**

It was shown in [7] that a matrix is stable (that is, all its eigenvalues are in the left half of the complex plane) if and only if it is a DH matrix. In this section, we obtain a parametrization of the sets $\Omega_{C}$, $\Omega_{V}$ and $\Omega_{D}$ in terms of DH matrices with extra constraint on $J$, $R$ and $Q$; see Section 2.1, 2.2 and 2.3, respectively. Note that for $\Omega$, the constraint $R\succeq 0$ can be removed to allow the region to intersect with the right half of the complex plane. This will allow us in particular to model $\Omega$-stability for discrete time system that correspond to $\Omega_{D}(0,1)$ or its intersection with $\Omega_{C}$ and $\Omega_{V}$. The following lemma will be frequently used in the following subsections.

Lemma 1**.**

Let $A=(J-R)Q$, where $J,R,Q\in{\mathbb{R}}^{n,n}$ such that $J^{T}=-J$, $R^{T}=R$, and $Q^{T}=Q$ is invertible. Let $\lambda\in{\mathbb{C}}$, and $v\in{\mathbb{C}}^{n}\setminus\{0\}$ be such that $v^{*}A=\lambda v^{*}$. Then

[TABLE]

Proof.

Given $v$ a left eigenvector of $A$ corresponding to eigenvalue $\lambda$, i.e., $v^{*}A=\lambda v^{*}$. This implies that

[TABLE]

and by taking the conjugate of (2.1), we get

[TABLE]

From (2.1) and (2.2) we have

[TABLE]

and

[TABLE]

2.1 Parametrization for conic sectors $\Omega_{C}$

Consider the $\Omega_{C}(a,\theta)$ region with parameters $a\in{\mathbb{R}}$ and $0\leq\theta\leq\pi/2$ and set $\alpha:=\sin(\theta)$ and $\beta:=\cos(\theta)$. To parametrize $\Omega_{C}(a,\theta)$ in terms of DH matrices, let us first prove the following elementary lemma.

Lemma 2**.**

Let $\lambda=\lambda_{1}+i\lambda_{2}$, where $\lambda_{1}$, $\lambda_{2}\in{\mathbb{R}}$. Then \left[\begin{array}[]{cc}\alpha\,(\lambda_{1}-a)&\beta\,i\lambda_{2}\\ -\beta\,i\lambda_{2}&\alpha\,(\lambda_{1}-a)\end{array}\right]\prec 0 if and only if $\lambda\in\Omega_{C}(a,\theta)$.

Proof.

The proof follows using the fact that \left[\begin{array}[]{cc}\alpha\,(\lambda_{1}-a)&\beta\,i\lambda_{2}\\ -\beta\,i\lambda_{2}&\alpha\,(\lambda_{1}-a)\end{array}\right] is Hermitian and therefore it is negative definite if and only if both eigenvalues $\mu_{1}=\alpha(\lambda_{1}-a)+\beta\lambda_{2}$ and $\mu_{2}=\alpha(\lambda_{1}-a)-\beta\lambda_{2}$ are negative which is true if and only if $\alpha(\lambda_{1}-a)<\beta\lambda_{2}<-\alpha(\lambda_{1}-a)$, i.e., $\lambda\in\Omega_{C}(a,\theta)$.

Theorem 1**.**

Let $A\in{\mathbb{R}}^{n,n}$. Then $A$ is $\Omega_{C}(a,\theta)$-stable if and only if $A=(J-R)Q$ for some $J,R,Q\in R^{n,n}$ such that $Q\succ 0$ and

[TABLE]

Proof.

First suppose that $A=(J-R)Q$ for some $J,R,Q$ satisfying $Q\succ 0$ and (2.3). Let $\lambda=\lambda_{1}+i\lambda_{2}$ be an eigenvalue of $A$ and let $v\in{\mathbb{C}}^{n}\setminus\{0\}$ be a left eigenvector of $A$ corresponding to eigenvalue $\lambda$. Since \left[\begin{array}[]{cc}\alpha\,(R+aQ^{-1})&-\beta\,J\\ \beta\,J&\alpha\,(R+aQ^{-1})\end{array}\right]\succ 0 and $v\neq 0$, we have that

[TABLE]

Thus by using Lemma 1 in (2.10), we obtain

[TABLE]

This implies that \left[\begin{array}[]{cc}\alpha\,(\lambda_{1}-a)&\beta\,i\lambda_{2}\\ -\beta\,i\lambda_{2}&\alpha\,(\lambda_{1}-a)\end{array}\right]\prec 0 since $Q$ is positive definite. Thus Lemma 2 implies that $A$ is $\Omega_{C}(a,\theta)$-stable.

For the ‘only if’ part, since $A$ is $\Omega_{C}(a,\theta)$-stable, by [3, Theorem 2.2], there exists $X\succ 0$ such that

[TABLE]

Let

[TABLE]

Then $(J-R)Q=A$ and it follows from (2.22) that

[TABLE]

As a consequence of (2.3) in Theorem 1, the matrix $J$ is skew-symmetric. However, the matrix $R$ may not be positive definite (when $a>0$) and therefore the $\Omega_{C}(a,\theta)$-stable matrix $A$ need not be a DH matrix. But when $a\leq 0$, then (2.3) implies that $R+aQ^{-1}\succ 0$, or equivalently, $R\succ-aQ^{-1}$ since $a\leq 0$ and $Q\succ 0$. As a result $R$ is positive semidefinite. Therefore in this case $A$ is $\Omega_{C}(a,\theta)$-stable if and only if $A$ is a DH matrix satisfying (2.3).

2.2 Parametrization for vertical strips $\Omega_{V}$

We can characterize $\Omega_{V}$-stability as follows.

Theorem 2**.**

Let $A\in{\mathbb{R}}^{n,n}$ and $h<k$. Then $A$ is $\Omega_{V}(h,k)$-stable if and only if $A=(J-R)Q$ for some $J,R,Q\in{\mathbb{R}}^{n,n}$ such that $J^{T}=-J$, $R^{T}=R$, $Q\succ 0$, and

[TABLE]

Proof.

First suppose that $A=(J-R)Q$, where $J^{T}=-J$, $R^{T}=R$, $Q\succ 0$ such that $kQ^{-1}\succ R\succ hQ^{-1}$. Let $\lambda$ be an eigenvalue of $A$ and $x\in{\mathbb{C}}^{n}\setminus\{0\}$ such that $Ax=\lambda x$ or $(J-R)Qx=\lambda x$. Since $Q$ is invertible, this implies that

[TABLE]

Since $x^{*}Qx>0$ as $Q\succ 0$ and $R$ satisfies $kQ^{-1}\succ R\succ hQ^{-1}$, we have $kx^{*}Qx>x^{*}QRQx>hx^{*}Qx$. This implies that

[TABLE]

From (2.24) and (2.25), we have that $-k<\mathop{\mathrm{Re}}{(\lambda)}<-h$.

Conversely, let $A$ be $\Omega_{V}(h,k)$-stable. Then from [3], there exists a single $X\succ 0$ such that

[TABLE]

Define

[TABLE]

Then clearly $A=(J-R)Q$. Also in view of (2.26) $R$ satisfies $kQ^{-1}\succ R\succ hQ^{-1}$.

It is easy to see that in Theorem 2 when $h\geq 0$ $A$ is $\Omega_{V}$-stable if and only if $A$ is a DH matrix since $R\succ hQ^{-1}$.

2.3 Parametrization for disks $\Omega_{D}$

We first recall a result from [3] and state it as a lemma in the following.

Lemma 3**.**

[3]* Consider the $\Omega_{D}(-q,r)$ region where $q\in{\mathbb{R}}$ and $r>0$, and let $\lambda\in{\mathbb{C}}$. Then $\lambda\in\Omega_{D}(-q,r)$ if and only if \left[\begin{array}[]{cc}-r&q+\lambda\\ q+\overline{\lambda}&-r\end{array}\right]\prec 0.*

We can characterize $\Omega_{D}$-stability as follows.

Theorem 3**.**

Let $A\in{\mathbb{R}}^{n,n}$, $q\in{\mathbb{R}}$ and $r>0$. Then $A$ is $\Omega_{D}(-q,r)$-stable if and only if $A=(J-R)Q$ for some $J,R,Q\in{\mathbb{R}}^{n,n}$ such that $J^{T}=-J$, $R$ is symmetric, $Q\succ 0$ and

[TABLE]

Proof.

First suppose that $A=(J-R)Q$ with $J^{T}=-J$, $R^{T}=R$ and $Q\succ 0$ such that (2.27) holds. Let $\lambda=\lambda_{1}+i\lambda_{2}$, where $\lambda_{1},\lambda_{2}\in{\mathbb{R}}$ be an eigenvalue of $A$ and $v\in{\mathbb{C}}^{n}\setminus\{0\}$ be a corresponding left eigenvector. Since (2.27) holds, we have that

[TABLE]

This implies that

[TABLE]

Since $v^{*}Q^{-1}v>0$ as $Q\succ 0$, we obtain

[TABLE]

Thus in view of Lemma 1, we have

[TABLE]

This implies by using Lemma 3 that $\lambda\in\Omega_{D}(-q,r)$ and therefore $A$ is $\Omega_{D}(-q,r)$-stable.

Conversely, suppose $A$ is $\Omega_{D}(-q,r)$-stable then by [3, Theorem 2.2] there exists $X\succ 0$ satisfying

[TABLE]

Define

[TABLE]

Then $A=(J-R)Q$ with $J^{T}=-J$, $R$ symmetric and $Q\succ 0$. Moreover, by (2.32), we have

[TABLE]

This implies that

[TABLE]

This completes the proof.

We note that in the above theorem, the matrix $R$ need not be positive semidefinite and thus a $\Omega_{D}$-stable matrix need not be a DH matrix. However, if the disc $\Omega_{D}$ completely lies in the left half of the complex plane, then $A$ is a DH matrix. More precisely we have the following result.

Theorem 4**.**

*Let $A\in{\mathbb{R}}^{n,n}$ and $q\geq r>0$. Then $A$ is $\Omega_{D}(-q,r)$-stable if and only if $A$ is a DH matrix such that (2.27) is satisfied. *

Proof.

The proof is immediate from Theorem 3 if we show that $A$ is $\Omega_{D}(-q,r)$-stable implies $R=-\frac{AX+XA^{T}}{2}$ is positive semidefinite. Let $x\in{\mathbb{C}}^{n}\setminus\{0\}$ so that $2x^{*}Rx=-x^{*}(AX)x-x^{*}(AX)^{T}x$. Using (2.32), we also have

[TABLE]

This implies that $-x^{*}(AX)x-x^{*}(AX)^{T}x~{}>~{}2(q-r)x^{*}Xx$ and hence $x^{*}Rx>2(q-r)x^{*}Xx\geq 0,$ since $X\succ 0$ and $q\geq r$ which completes the proof.

3 Reformulation of the nearest $\Omega$-stable matrix problem

In this section, we reformulate the problem (1.1) of finding the nearest $\Omega$-stable matrix to a given unstable matrix $A$ into an equivalent optimization problem with a convex feasible set. In the reformulations of $\Omega_{V}$ and $\Omega_{D}$, we used linear matrix inequalities (LMIs) involving the inverse of $Q^{-1}$, which is not convex; see (2.3), (2.23) and (2.27). Let us introduce the auxiliary variable $P=Q^{-1}\succ 0$. In view of Theorems 1, 2 and 3, this allows us to parametrize the sets $\mathbb{S}_{\Omega_{C}(a,\theta)}^{n,n}$, $\mathbb{S}_{\Omega_{V}(h,k)}^{n,n}$, and $\mathbb{S}_{\Omega_{D}(-q,r)}^{n,n}$ as follows:

[TABLE]

[TABLE]

and

[TABLE]

where $0\leq\theta<\pi/2$, $h<k$ and $r>0$. Note that these sets are non-convex and open. For example, the set of stable matrices (or, equivalently $\Omega_{C}(0,\theta)$ matrices with $\theta=\pi/2$) is non-convex [15]. It can be checked that for any $\theta\in(0,\pi/2)$, the set $\Omega_{C}(0,\theta)$ is also non-convex. Indeed, consider A=\left[\begin{array}[]{cc}\alpha\,(-\beta)&\beta\,i(\alpha)\\ 0&\alpha\,(-\beta)\end{array}\right] and B=\left[\begin{array}[]{cc}\alpha\,(-\beta)&0\\ -\beta\,i(\alpha)&\alpha\,(-\beta)\end{array}\right] with $A,\,B\in\mathbb{S}_{\Omega_{\theta}}^{n,n}$. For $\gamma=1/2$, the matrix $\gamma A+(1-\gamma)B\notin\mathbb{S}_{\Omega_{\theta}}^{n,n}$ as it has an eigenvalue at zero. Note also that the above sets are not closed (because of the constraints of positive definiteness). From an optimization point of view, it does not make much sense to optimize on such sets since the optimal solution(s) may not be attained. Therefore, we will consider the closure of these sets: this amounts to replace all constraints involving a positive definite constraint with a positive semidefinite constraint, that is, replace $\succ 0$ with $\succeq 0$, in the definition of the sets (3.1), (3.2) and (3.3). We will denote the corresponding sets as $\bar{\mathbb{S}}_{\Omega_{C}(a,\theta)}^{n,n}$, $\bar{\mathbb{S}}_{\Omega_{V}(h,k)}^{n,n}$, and $\bar{\mathbb{S}}_{\Omega_{D}(-q,r)}^{n,n}$, respectively. Note that by considering the closure of these sets, as done in [7], we do no change the value of the infimum of (1.1).

Finally, given $0<\theta\leq\pi/2$, $h<k$ and several disks of parameters $(q_{i},r_{i})$ $1\leq i\leq k$, we tackle (1.1) by solving

[TABLE]

where

[TABLE]

The feasible set of the above optimization problem only involves convex LMI constraints. Of course, the objective function is non-convex and the problem remains difficult, but it is easier to handle non-convex objective rather than a non-convex feasible set.

4 Algorithm and numerical experiments

In this section, we describe a block coordinate descent method to tackle (3.4) and then apply it on some selected examples. To solve the convex optimization subproblems involving LMIs, we use the interior point method SDPT3 (version 4.0) [16, 18] with CVX as a modeling system [4, 8]. Our code is available from https://sites.google.com/site/nicolasgillis/ and the numerical examples presented below can be directly run from this online code. All tests are preformed using Matlab R2015a on a laptop Intel CORE i7-7500U CPU @2.7GHz 24Go RAM.

4.1 Block coordinate descent method for (3.4)

In this section, we propose to use a block coordinate descent (BCD) method to tackle (3.4). The BCD method will alternatively optimize the block of variable $(J,R)$ and $P$:

• •

For $P$ fixed, the subproblem in $(J,R)$ is convex since ${\|A-(J-R)P^{-1}\|}_{F}^{2}$ is a quadratic function of $(J,R)$, and the feasible set is convex in variables $(J,R,P)$. More precisely, this is a linear least squares problem under LMI constraints.

• •

For $(J,R)$ fixed, the problem in variable $P$ is non-convex, because of the objective function ${\|A-(J-R)P^{-1}\|}_{F}^{2}$. In a standard BCD method, one would optimize solely on variable $P$ for $(J,R)$ fixed. However, we have observed that performing a gradient step on all the variables performs better. Denoting $f(J,R,P)={\|A-(J-R)P^{-1}\|}_{F}^{2}$ and $E=(J-R)P^{-1}-A$, the gradient of $f$ is given by

[TABLE]

For the computation of the gradient with respect to $P$, we refer to [6, Appendix] for a similar derivation. For the step length used in the gradient method, we use a backtracking line search, that is, we reduce the step length as long as the objective function increases, while allowing it to increase at the next step.

Note that we have also implemented a projected fast gradient method, as done in [7], but it does not perform as well as BCD. The reason is that, in [7], the feasible set is simply $R\succeq 0$ and $Q=P^{-1}\succeq 0$ which can be projected onto very efficiently with an eigenvalue decomposition. For general $\Omega$ stability, projecting $(J,R,P)$ onto the feasible set roughly has the same computational cost as optimizing exactly over variables $(J,R)$ for $P$ fixed. For all the numerical experiments presented below, we will use 100 iterations of the above scheme.

4.2 Initialization

As expected, our algorithm will be sensitive to the initial choice of the matrices $(J,R,P)$ since we are tackling the difficult non-convex optimization problem (3.4). In this paper, we propose two initializations.

The first initialization, which we will refer to as identity initialization chooses $P=I_{n}$ and then, for $P$ fixed, chooses $J$ and $R$ that minimize (3.4). We will denote this initialization $(J_{i},R_{i},P_{i})$ and the solution obtained by BCD with this initialization as $(J_{i}^{*},R_{i}^{*},P_{i}^{*})$.

The second initialization, which we will refer to as the LMI initialization, will use a solution to a relaxation of the LMIs (2.22), (2.26) and (2.32). In fact, since the input matrix $A$ is in general not $\Omega$-stable, these LMIs do not admit a feasible solution $X\succ 0$. Hence we will replace all inequalities $\succ 0$ with $\succ-\delta I$ and minimize $\delta$ over variables $X\succ I$ and $\delta$. Denoting $(X^{*},\delta^{*})$ an optimal solution of these relaxed LMIs, we take

[TABLE]

as specified in Theorems 1, 2 and 3. Note that the optimal solutions $(X^{*},\delta^{*})$ will be such that $\delta^{*}=0$ if and only if $A$ is $\Omega$-stable in which case the problem would be solved. This is an advantage of this initialization: it will detect whether the input matrix is $\Omega$-stable. We will denote this initialization $(J_{x},R_{x},P_{x})$ and the solution obtained by BCD with this initialization as $(J_{x}^{*},R_{x}^{*},P_{x}^{*})$.

4.3 Synthetic data sets

Let us consider

[TABLE]

so that

[TABLE]

The set $\Omega$ is illustrated on Figure 4.1. To generate a matrix that is not too far from being $\Omega$-stable, we proceed as follows. First, we generate each entry of $J_{0}$, $R_{0}$ and $Q_{0}$ at random using the Gaussian distribution of mean 0 and standard deviation 1 (randn(n) in Matlab). Then we perform the following projection

[TABLE]

to obtain an $\Omega$-stable matrix $A_{t}=(J_{t}-R_{t})P_{t}^{-1}$. Finally, we generate a matrix $N$ where each entry is generated at random using the Gaussian distribution of mean 0 and standard deviation 1, and set

[TABLE]

so that ${\|A_{t}-{A}\|}_{F}/{\|A_{t}\|}_{F}=\epsilon$. The matrix $A$ will in general not be $\Omega$-stable for $\epsilon$ sufficiently large.

In the following, we compare the solutions obtained by our BCD algorithm with the two initializations described in the previous section. We also consider the initialization made of the true $(J_{t},R_{t},Q_{t})$ used to generate $A$, which we will refer to as the true initialization. This will be useful to validate the other two initializations and see whether they are able to obtain good nearby $\Omega$-stable matrices.

Figure 4.1 displays the result of a particular example with $\epsilon=0.1$. We observe that the three initializations lead to different solutions. The best solution (with relative error $\frac{{\|A-(J-R)P^{-1}\|}_{F}}{{\|A\|}_{F}}=2.74\%$) is found with the identity initialization although it has, as expected, the highest initial relative error (41.33%). The LMI initialization has the lowest initial error (5.04%), lower than the true initialization111The true initialization $(J_{t}-R_{t})P_{t}^{-1}=A$ has relative error $\epsilon=10\%$ by construction. However, in BCD, we only use $P$ for the initialization and set $(J,R)$ as the optimal solution for $P$ fixed. (5.64%) but BCD is not able to improve it much (up to 3.9%).

Let us now perform more extensive numerical experiments with these randomly generated matrices. We consider $\epsilon=0.01,0.05,0.1,0.2$, and for each value of $\epsilon$, we generate 10 matrices $A$. Table 4.1 reports the average error and standard deviation for each noise level for these 10 randomly generated matrices. It also reports in bracket the number of times each initialization obtained the best solution (up to 0.01%).

We observe the following

• •

For low noise levels ($\epsilon\leq 0.05$), the LMI initialization performs in average the best, although the identity initialization finds in most of the cases the best solution. This means that, in several cases, the identity initialization is not able to find a good solution. Hence the LMI initialization is more reliable (the standard deviation is smaller for $\epsilon\leq 0.05$). This was expected as the LMI initialization provides an exact solution for $\Omega$-stable matrices and, for small noise levels, it is expected that the relaxed LMIs are close to the original LMIs. Also, recall that the LMI initialization obtains a low error prior to BCD performing any iteration while the identity initialization requires several iterations before achieving a low error; see Figure 4.1 for an illustration.

• •

For larger noise levels ($\epsilon\geq 0.1$), the identity and LMI initializations perform similarly, although the identity initialization finds in more cases the best solution.

• •

In all cases, rather surprisingly, the true initialization does not perform well. For $\epsilon=0.01$, it performs better than the identity initialization while, in all other cases, it is worse than the other two initializations. This is rather surprising, and shows that our two proposed initializations are performing well.

4.4 Discrete-time stability: $\Omega=\Omega_{D}(0,1)$

Stability of discrete-time LTI systems requires the eigenvalues of the matrix $A$ defining the system to belong to the unit disk, that is, to $\Omega_{D}(0,1)$; see, e.g., [15, 10, 5] and the references therein. Therefore, our algorithm can be used to find a nearby stable system for discrete LTI systems. We illustrate this on the example from [10, Section 4.4]:

[TABLE]

with $\rho(A)=\max_{i}|\lambda_{i}(A)|=2.4$. The nonnegative solution provided by the authors with their algorithm is

[TABLE]

with error ${\|A-A_{+}\|}_{F}=1.10$. In [5], the best reported solution is

[TABLE]

with error ${\|A-A_{b}\|}_{F}=0.76$. The solution provided by BCD with the identity initialization is

[TABLE]

with error ${\|A-(J-R)Q\|}_{F}=0.90$, while with the LMI initialization, BCD gives

[TABLE]

with error ${\|A-(J_{x}^{*}-R_{x}^{*})Q_{x}^{*}\|}_{F}=1.40$. Figure 4.2 displays the eigenvalues of the different solutions: again we see that depending on the initialization and the algorithm used, we obtained rather different solutions. Surprisingly, although the LMI initialization achieves the highest approximation error of $A$, it provides an optimal approximation of its eigenvalues: four are perfectly recovered while the last one is approximated by its projection on the unit disk. Hence, although this locally optimal solution does not have a low Frobenius norm error, it has an interesting structure.

5 Conclusion and further work

In this paper, we have proposed a new parametrization of $\Omega$-stable matrices, where $\Omega$ is the intersection of several regions of the complex plane, namely conic sectors, vertical strips and disks centred on the real line. This allowed us to propose an algorithm to tackle the nearest $\Omega$-stable matrix problem where one is given a matrix $A$ and is looking for the nearest $\Omega$-stable matrix. We illustrated the effectiveness of this approach on several examples.

Further work include the design of faster algorithms to tackle (3.4). In fact, our BCD algorithm currently relies on an interior point method to solve the LMIs which does not scale well. For example, on a standard laptop, $n$ can be up to about 50 for which one iterations takes about 1 minute. Another direction of research is to extended our approach to other regions of the complex plane. For example, LMI regions, that is, subsets of the complex plane that is representable by LMIs [3], would be of particular interest.