Optimal robustness of port-Hamiltonian systems

TL;DR

This paper develops methods to construct port-Hamiltonian systems with maximal robustness and passivity, providing algorithms for optimal and near-optimal solutions to ensure system passivity and stability.

Contribution

It introduces a systematic approach to realize the most robust port-Hamiltonian systems from a given transfer function, including algorithms for optimal and suboptimal passive system approximations.

Findings

Maximal passivity radius realized by normalized port-Hamiltonian systems.

Algorithm for constructing optimal port-Hamiltonian realization via linear matrix inequalities.

Suboptimal method for approximating the nearest passive system to a non-passive one.

Abstract

We construct optimally robust port-Hamiltonian realizations of a given rational transfer function that represents a passive system. We show that the realization with a maximal passivity radius is a normalized port-Hamiltonian one. Its computation is linked to a particular solution of a linear matrix inequality that defines passivity of the transfer function, and we provide an algorithm to construct this optimal solution. We also consider the problem of finding the nearest passive system to a given non-passive one and provide a simple but suboptimal solution.

Figures

Equations

$$\begin{array}[]{rcl}\dot{x}&=&Ax+Bu,\ x(0)=0,\\ y&=&Cx+Du,\end{array}$$

$$\mathcal{T}(s)=D+C(sI-A{)}^{-1}B.$$

$$\Phi (s):={\mathcal{T}}^{\mathsf{H}}(-s)+\mathcal{T}(s)$$

$$W(X,\mathcal{M}):=\left[\begin{array}[]{cc}-A^{\mathsf{H}}X-X\,A&C^{\mathsf{H}}-X\,B\\ C-B^{\mathsf{H}}X&D^{\mathsf{H}}+D\end{array}\right],$$

$$W(X)\ge 0$$

$${\mathbb{X}}^{>}:=\left\{X\in {\mathbb{H}}_n|W(X)\ge 0,X>0\right\},$$

$${\mathbb{X}}^{\gg }:=\left\{X\in {\mathbb{H}}_n|W(X)>0,X>0\right\}.$$

$$\mathcal{H}(x(t_1))-\mathcal{H}(x(t_0))\le {\int }_{t_0}^{t_1}\Re (y(t{)}^{\mathsf{H}}u(t))dt$$

$$\mathsf{Ricc}(X):=-XA-A^{\mathsf{H}}X-(C^{\mathsf{H}}-XB)(D^{\mathsf{H}}+D{)}^{-1}(C-B^{\mathsf{H}}X)=0.$$

$$H:=\left[\begin{array}[]{cc}A-B(D^{\mathsf{H}}+D)^{-1}C&-B(D^{\mathsf{H}}+D)^{-1}B^{\mathsf{H}}\\ C^{\mathsf{H}}(D^{\mathsf{H}}+D)^{-1}C&-(A-B(D^{\mathsf{H}}+D)^{-1}C)^{\mathsf{H}}\end{array}\right],$$

$$S(s):=\left[\begin{array}[]{cc|c}0&A-sI_{n}&B\\ A^{\mathsf{H}}+sI_{n}&0&C^{\mathsf{H}}\\ \hline\cr B^{\mathsf{H}}&C&D^{\mathsf{H}}+D\end{array}\right].$$

$$\begin{array}[]{rcl}\dot{x}&=&(J-R)Qx+(G-K)u,\\ y&=&(G+K)^{\mathsf{H}}Qx+(S+N)u,\end{array}$$

$$\mathcal{V}:=\left[\begin{array}[]{cccc}J&G\\ -G^{\mathsf{H}}&N\end{array}\right]=-\mathcal{V}^{\mathsf{H}},\ \mathcal{W}:=\left[\begin{array}[]{cccc}R&K\\ K^{\mathsf{H}}&S\end{array}\right]=\mathcal{W}^{\mathsf{H}}\geq 0,\ Q=Q^{\mathsf{H}}>0.$$

$$\{A_T,B_T,C_T,D_T\}:=\{TA{T}^{-1},TB,C{T}^{-1},D\}$$

$$\displaystyle\left[\begin{array}[]{cccc}T^{-\mathsf{H}}&0\\ 0&I_{m}\end{array}\right]\left[\begin{array}[]{cccc}-A^{\mathsf{H}}X-XA&C^{\mathsf{H}}-XB\\ C-B^{\mathsf{H}}X&D^{\mathsf{H}}+D\end{array}\right]\left[\begin{array}[]{cccc}T^{-1}&0\\ 0&I_{m}\end{array}\right]$$

$$\displaystyle\qquad=\left[\begin{array}[]{cccc}-A_{T}&-B_{T}\\ C_{T}&D_{T}\end{array}\right]+\left[\begin{array}[]{cccc}-A_{T}^{\mathsf{H}}&C^{\mathsf{H}}_{T}\\ -B^{\mathsf{H}}_{T}&D^{\mathsf{H}}_{T}\end{array}\right]\geq 0.$$

$$\mathcal{S}:=\left[\begin{array}[]{cccc}-A_{T}&-B_{T}\\ C_{T}&D_{T}\end{array}\right]$$

$$\left[\begin{array}[]{cccc}R&K\\ K^{\mathsf{H}}&S\end{array}\right]:=\frac{{\mathcal{S}}+{\mathcal{S}}^{\mathsf{H}}}{2}\geq 0,\quad\left[\begin{array}[]{cccc}J&G\\ -G^{\mathsf{H}}&N\end{array}\right]:=\frac{{\mathcal{S}}-{\mathcal{S}}^{\mathsf{H}}}{2}.$$

$$W(X,\mathcal{M})=\left[\begin{array}[]{cc}0&C^{\mathsf{H}}\\ C&D^{\mathsf{H}}+D\end{array}\right]-\left[\begin{array}[]{cc}A^{\mathsf{H}}X+X\,A&X\,B\\ B^{\mathsf{H}}X&0\end{array}\right]$$

$$W(X,\mathcal{M}+{\Delta }_{\mathcal{M}})$$

$$\|\Delta_{\mathcal{M}}\|_{F}:=\left\|\left[\begin{array}[]{ccc}\Delta_{A}&\Delta_{B}\\ \Delta_{C}&\Delta_{D}\end{array}\right]\right\|_{F},\quad\|\Delta_{\mathcal{M}}\|_{2}:=\left\|\left[\begin{array}[]{ccc}\Delta_{A}&\Delta_{B}\\ \Delta_{C}&\Delta_{D}\end{array}\right]\right\|_{2},$$

$${\rho }_{\mathcal{M}}(X):=\underset{{\Delta }_{\mathcal{M}}\in {\mathbb{C}}^{n+m,n+m}}{\inf }\left\{\Vert {\Delta }_{\mathcal{M}}\Vert |\det W(X,\mathcal{M}+{\Delta }_{\mathcal{M}})=0\right\}.$$

$${\rho }_{\mathcal{M}}:=\underset{X\in {\mathbb{X}}^{\gg }}{\sup }\underset{{\Delta }_{\mathcal{M}}\in {\mathbb{C}}^{n+m,n+m}}{\inf }\{\Vert {\Delta }_{\mathcal{M}}\Vert |\det W(X,\mathcal{M}+{\Delta }_{\mathcal{M}})=0\}=\underset{X\in {\mathbb{X}}^{\gg }}{\sup }{\rho }_{\mathcal{M}}(X).$$

$$\displaystyle\left[\begin{array}[]{cc}-X&0\\ 0&I_{m}\end{array}\right]\left[\begin{array}[]{cc}A+\Delta_{A}&B+\Delta_{B}\\ C+\Delta_{C}&D+\Delta_{D}\end{array}\right]$$

$$\displaystyle\qquad+\left[\begin{array}[]{cc}A^{\mathsf{H}}+\Delta_{A}^{\mathsf{H}}&C^{\mathsf{H}}+\Delta_{C}^{\mathsf{H}}\\ B^{\mathsf{H}}+\Delta_{B}^{\mathsf{H}}&D^{\mathsf{H}}+\Delta_{D}^{\mathsf{H}}\end{array}\right]\left[\begin{array}[]{cc}-X&0\\ 0&I_{m}\end{array}\right]>0.$$

$$\hat{W}:=W(X),\;\hat{X}:=\left[\begin{array}[]{cc}X&0\\ 0&I_{m}\end{array}\right],\;\Delta_{\mathcal{S}}:=\left[\begin{array}[]{cc}-\Delta_{A}&-\Delta_{B}\\ \Delta_{C}&\Delta_{D}\end{array}\right],$$

$$W(X,\mathcal{M}+{\Delta }_{\mathcal{M}})=\hat{W}+\hat{X}{\Delta }_{\mathcal{S}}+{\Delta }_{\mathcal{S}}^{\mathsf{H}}\hat{X}>0$$

$$M(\gamma):=\left[\begin{array}[]{cc}\gamma\hat{X}\hat{W}^{-\frac{1}{2}}\\ \hat{W}^{-\frac{1}{2}}/\gamma\end{array}\right]\left[\begin{array}[]{cc}\gamma\hat{W}^{-\frac{1}{2}}\hat{X}&\hat{W}^{-\frac{1}{2}}/\gamma\end{array}\right]$$

$$M(\underline{\gamma})z=\underline{\lambda}_{\max}z,\quad z:=\left[\begin{array}[]{cc}u\\ v\end{array}\right],$$

$${\gamma }^2{\hat{W}}^{-\frac{1}{2}}{\hat{X}}^2{\hat{W}}^{-\frac{1}{2}}+{\hat{W}}^{-1}/{\gamma }^2.$$

$${\rho }_{\mathcal{M}}(X)=\underset{\gamma }{\min }\Vert {\Delta }_{\mathcal{S}}{\Vert }_2=\underset{\gamma }{\min }\Vert {\Delta }_{\mathcal{S}}{\Vert }_F\ge \alpha \beta /2.$$

$$\alpha \beta /2\le {\rho }_{\mathcal{M}}(X)\le \alpha \beta /(1+|v^{\mathsf{H}}w|),$$

Full text

{AMS}

93D09, 93C05, 49M15, 37J25

1 Introduction

We consider realizations of linear dynamical systems that are denoted as positive-real or passive and their associated transfer functions. In particular, we study positive real transfer functions which play a fundamental role in systems and control theory: they represent e. g., spectral density functions of stochastic processes and arise in spectral factorizations. Positive-real transfer functions form a convex set, and this property has lead to the extensive use of convex optimization techniques in this area [8]. Passive systems and their relationships with positive-real transfer functions are well studied, starting with the works [18, 22, 29, 30] and the topic has recently received considerable attention in the context of port-Hamiltonian (pH) system models, [24, 27].

In this paper we show that in the set of continuous-time pH realizations of positive-real transfer functions, there is a subset that achieves optimal robustness, in the sense that their passivity radius is maximal. Considering the Laplace transform of the linear time-invariant system

[TABLE]

denoted as $\mathcal{M}:=\{A,B,C,D\}$, the transfer function is given by

[TABLE]

Here $u:\mathbb{R}\to\mathbb{C}^{m}$, $x:\mathbb{R}\to\mathbb{C}^{n}$, and $y:\mathbb{R}\to\mathbb{C}^{m}$ are vector-valued functions denoting, respectively, the input, state, and output of the system. Denoting real and complex $n$-vectors ($n\times m$ matrices) by $\mathbb{R}^{n}$, $\mathbb{C}^{n}$ ($\mathbb{R}^{n\times m}$, $\mathbb{C}^{n\times m}$), respectively, the coefficient matrices satisfy $A\in\mathbb{C}^{n\times n}$, $B\in\mathbb{C}^{n\times m}$, $C\in\mathbb{C}^{m\times n}$, and $D\in\mathbb{C}^{m\times m}$.

We restrict ourselves to systems which are minimal, i. e., the pair $(A,B)$ is controllable (for all $s\in\mathbb{C}$, \operatorname{rank}\mbox{\small[,sI-A;|;B,]}=n), and the pair $(A,C)$ is observable (i. e., $(A^{\mathsf{H}},C^{\mathsf{H}})$ is controllable). Here, the Hermitian transpose and the transpose of a vector or matrix $V$ is denoted by $V^{\mathsf{H}}$ and $V^{\mathsf{T}}$, respectively, and the identity matrix is denoted by $I_{n}$ or $I$ if the dimension is clear. Throughout this article we will use the following notation. We denote the set of Hermitian matrices in $\mathbb{C}^{n\times n}$ by ${\mathbb{H}_{n}}$. Positive definiteness (semidefiniteness) of $A\in{\mathbb{H}_{n}}$ is denoted by $A>0$ ($A\geq 0$). The real and imaginary parts of a complex matrix $Z$ are written as $\Re(Z)$ and $\Im(Z)$, respectively, and $\imath$ is the imaginary unit. We consider functions over ${\mathbb{H}_{n}}$, which is a vector space if considered as a real subspace of $\mathbb{R}^{n\times n}+\imath\mathbb{R}^{n\times n}$.

The paper is organized as follows. In Section 2 we characterize the classes of passive systems and of port-Hamiltonian (pH) systems. We then show in Section 3 the relevance of pH systems in estimating the passivity radius of passive systems and construct in Section 4 realizations with optimal robustness margin for passivity. In Section 5 we describe an algorithm to compute this optimal robustness margin. In Sections 6 and 7 we show how to use these ideas to estimate the distance to the set of passive systems and to the set of stable systems.

2 Passive systems and port-Hamiltonian realizations

The concepts of positive-realness and passivity are well studied. We briefly recall some important properties following [28], and refer to the literature for a more detailed survey. Consider a continuous-time system as in (1) and its transfer function $\mathcal{T}(s)$ as in (2). The transfer function $\mathcal{T}(s)$ is called positive-real if the matrix-valued rational function

[TABLE]

is positive semidefinite for $s$ on the imaginary axis, i. e., $\Phi(\imath\omega)\geq 0$ for all $\omega\in\mathbb{R}$ and it is called strictly positive-real if $\Phi(\imath\omega)>0$ for all $\ \omega\in\mathbb{R}$.

For $X\in{\mathbb{H}_{n}}$ and a system $\mathcal{M}=\{A,B,C,D\}$, we consider the matrix function

[TABLE]

which we also denote by $W(X)$, when it is clear from the context which model we refer to. If $\mathcal{T}(s)$ is positive-real, then there exists $X\in{\mathbb{H}_{n}}$ such that the linear matrix inequality (LMI)

[TABLE]

holds. In this context we will make frequent use of the sets

[TABLE]

A system $\mathcal{M}:=\left\{A,B,C,D\right\}$ is called passive if there exists a state-dependent storage function, $\mathcal{H}(x)\geq 0$, such that for any $t_{1},t_{0}\in\mathbb{R}$ with $t_{1}>t_{0}$, the dissipation inequality

[TABLE]

holds. If for all $t_{1}>t_{0}$, the inequality in (6) is strict then the system is called strictly passive.

If $D^{\mathsf{H}}+D$ is invertible, then the minimum rank solutions of (4) in {\mathbb{X}^{\raisebox{1.39998pt}{{>}}}} are those for which $\operatorname{rank}W(X)=\operatorname{rank}(D^{\mathsf{H}}+D)=m$, which in turn is the case if and only if the Schur complement of $D^{\mathsf{H}}+D$ in $W(X)$ is zero. This Schur complement is associated with the continuous-time algebraic Riccati equation (ARE)

[TABLE]

Solutions $X$ to (7) yield a spectral factorization of $\Phi(s)$, and each solution corresponds to an invariant subspace spanned by the columns of U:=\left[\begin{array}[]{cc}I_{n}&-X^{\mathsf{T}}\end{array}\right]^{\mathsf{T}} that remains invariant under multiplication with the Hamiltonian matrix

[TABLE]

i. e., $U$ satisfies $HU=UA_{F}$ for a closed loop matrix $A_{F}=A-BF$ with $F:=(D^{\mathsf{H}}+D)^{-1}(C-B^{\mathsf{H}}X)$, see e.g., [11].

We can also associate with $\Phi$ a system pencil

[TABLE]

Then the Schur complement of $S(s)$ is the transfer function $\Phi(s)$ and the generalized eigenvalues of $S(s)$ are the zeros of $\Phi(s)$. The properties of the system can, however, be checked in a much more numerically robust way using the pencil $S(s)$ rather than the matrix $H$. This is, in particular, true if $D^{\mathsf{H}}+D$ is singular or ill-conditioned, see [6] for a detailed analysis and appropriate algorithms.

A special class of realizations of passive systems is that of port-Hamiltonian systems.

Definition 2.1**.**

A linear time-invariant port-Hamiltonian (pH) system has the state-space form

[TABLE]

and the system matrices satisfy the symmetry conditions

[TABLE]

Port-Hamiltonian systems were introduced from a different point of view in [24], but they also have a storage function and satisfy a dissipation inequality, and hence they are passive. Thus, there must be a coordinate transformation between a passive system and a representation (10) as a pH system. We briefly recall the construction of such a possible transformation.

Consider a minimal state-space model $\mathcal{M}:=\{A,B,C,D\}$ of a passive linear time-invariant system and let X\in{\mathbb{X}^{\raisebox{1.39998pt}{{>}}}} be a solution of the LMI (4). We then use a symmetric factorization $X=T^{\mathsf{H}}T$, which implies the invertibility of $T$, and define a new realization

[TABLE]

so that

[TABLE]

We can then use the Hermitian and skew-Hermitian part of the matrix

[TABLE]

to define the coefficients of a pH representation via

[TABLE]

This construction yields $Q=I_{n}$ because of the chosen factorization $X=T^{\mathsf{H}}T$. Note that the factor $T$ is unique up to a unitary factor $U$, since $T^{\mathsf{H}}U^{\mathsf{H}}UT=T^{\mathsf{H}}T$, but this factor $U$ will not affect the results described in this paper.

There is a lot of freedom in the representation of the system, since we could have used any matrix $X$ from the set {\mathbb{X}^{\raisebox{1.39998pt}{{>}}}}, or we could have chosen a representation where $Q$ was not the identity matrix. In the remainder of this paper, we will restrict ourselves to pH models where $Q=I_{n}$. The freedom remaining is thus the choice of the matrix $X$ from {\mathbb{X}^{\raisebox{1.39998pt}{{>}}}}, which we will use to make the representation ’maximally’ robust or well-conditioned to perturbations.

Remark 2.2**.**

*We stress that when the model $\mathcal{M}$ is real, then all the definitions and properties discussed above still hold. Moreover, the sets {\mathbb{X}^{\raisebox{1.39998pt}{{>}}}} and {\mathbb{X}^{\raisebox{1.39998pt}{{\gg}}}} can be constrained to be real without altering any of the results, since the real part $X_{\Re}$ of a Hermitian matrix $X$ is symmetric and positive (semi-)definite whenever $X$ is Hermitian and (semi-)definite. When the model $\mathcal{M}$ is real it therefore follows that whenever $W(X)\geq 0$ (or $W(X)>0$) then we also have that $W(X_{\Re})\geq 0$ (or $W(X_{\Re})>0$), and it then suffices to verify these conditions over the real symmetric matrices only. When doing that, the corresponding pH realizations will also be real. Finally, we point out that the extremal solutions $X_{-}$ and $X_{+}$ of the Riccati equations are also real when the model $\mathcal{M}$ is real. *

3 The passivity radius

Our goal to achieve robust pH representations of a passive system can be realized in different ways. A natural measure for this optimality is a large passivity radius $\rho_{\mathcal{M}}$, which is the smallest perturbation (in an appropriate norm) to the coefficients of a model $\mathcal{M}$ that makes the system non-passive. In this section we recall and extend a few results on passivity radii from [5] that will be employed in the next section.

Once we have determined a solution X\in{\mathbb{X}^{\raisebox{1.39998pt}{{>}}}} to the LMI (4), we can determine the representations (10) as discussed in Section 2 and the system is automatically passive (but not necessarily strictly passive). For each such representation we can determine the passivity radius and then choose the solution X\in{\mathbb{X}^{\raisebox{1.39998pt}{{>}}}} which is most robust under perturbations ${\Delta_{\mathcal{M}}}$ of the model parameters $\mathcal{M}:=\{A,B,C,D\}$. This is a very suitable approach for perturbation analysis, since as soon as we fix $X$, the matrix

[TABLE]

is linear in the unknowns $A,B,C,D$ and when we perturb the coefficients, then we preserve strict passivity as long as

[TABLE]

Hence, given X\in{\mathbb{X}^{\raisebox{1.39998pt}{{\gg}}}}, we can look for the smallest perturbation $\Delta_{\mathcal{M}}$ to our model $\mathcal{M}$ that makes $\det W(X,\mathcal{M}+\Delta_{\mathcal{M}})=0$. To measure the size of the perturbation $\Delta_{\mathcal{M}}$ of a state space model $\mathcal{M}$ we will use the Frobenius norm or the spectral norm

[TABLE]

and we use also the following $X$-passivity radius, which was introduced in [5] and gives a bound for the usual passivity radius.

Definition 3.1**.**

For X\in{\mathbb{X}^{\raisebox{1.39998pt}{{\gg}}}} the $X$-passivity radius is defined as

[TABLE]

Note that in order to compute $\rho_{\mathcal{M}}(X)$ for the model $\mathcal{M}$, we must have a matrix X\in{\mathbb{X}^{\raisebox{1.39998pt}{{\gg}}}}, since $W(X,\mathcal{M})$ must be positive definite to start with and also $X$ should be positive definite to obtain a state-space transformation from it. The following relation between the $X$-passivity radius and the usual passivity radius was also given in [5].

Lemma 3.2**.**

The passivity radius for a given model $\mathcal{M}$ satisfies

[TABLE]

We now provide an exact formula for the $X$-passivity radius based on a one-parameter optimization problem. For this, we rewrite the condition $W(X,\mathcal{M}+\Delta_{\mathcal{M}})>0$ as

[TABLE]

Setting

[TABLE]

inequality (35) can be written as the LMI

[TABLE]

as long as the system is still passive. In order to violate this condition, we need to find the smallest $\Delta_{\mathcal{S}}$ such that $\det W(X,\mathcal{M}+\Delta_{\mathcal{M}})=0$. The following theorem, based on results of [21] and [5], gives the minimal perturbation $\Delta_{\mathcal{S}}$ in both the Frobenius norm and the spectral norm. We point out that the definition of $\Delta_{\mathcal{S}}$ yields that $\|\Delta_{\mathcal{S}}\|_{2}=\|\Delta_{\mathcal{M}}\|_{2}$ and $\|\Delta_{\mathcal{S}}\|_{F}=\|\Delta_{\mathcal{M}}\|_{F}$.

Theorem 3.3**.**

Consider the matrices $\hat{X},\hat{W}$ in (36) and the pointwise positive semidefinite matrix function

[TABLE]

in the real parameter $\gamma$. Then the largest eigenvalue $\lambda_{\max}(M(\gamma))$ is a unimodal function of $\gamma$, (i.e. it is first monotonically decreasing and then monotonically increasing with growing $\gamma$). At the minimizing value $\underline{\gamma}$, $M(\underline{\gamma})$ has an eigenvector $z$, i.e.

[TABLE]

*where $\|u\|_{2}^{2}=\|v\|_{2}^{2}=1$. The minimum norm perturbation $\Delta_{\mathcal{S}}$ is of rank $1$ and is given by $\Delta_{\mathcal{S}}=uv^{\mathsf{H}}/\underline{\lambda}_{\max}$. It has norm $1/\underline{\lambda}_{\max}$ both in the spectral norm and in the Frobenius norm. *

In [5] also the following simple bound for $\underline{\lambda}_{\max}$ was derived.

Corollary 3.4**.**

Consider the matrices $\hat{X},\hat{W}$ in (36) and the pointwise positive semidefinite matrix function $M(\gamma)$ as in (38). The largest eigenvalue of $M(\gamma)$ is also the largest eigenvalue of

[TABLE]

An upper bound for $\underline{\lambda}_{\max}$ is given by $\underline{\lambda}_{\max}\leq\frac{2}{\alpha\beta}$ where $\alpha^{2}:=\lambda_{\min}(\hat{W})$ and $\beta^{2}=\lambda_{\min}(\hat{X}^{-1}\hat{W}\hat{X}^{-1})$. The corresponding lower bound for $\|\Delta_{\mathcal{S}}\|_{2}$ and $\|\Delta_{\mathcal{S}}\|_{F}$ is given by

[TABLE]

The following theorem, also proven in [5], constructs a rank one perturbation which makes the matrix $W_{\Delta_{\mathcal{M}}}$ singular and therefore gives an upper bound for $\rho_{M}(X)$.

Theorem 3.5**.**

Let $\mathcal{M}=\{A,B,C,D\}$ be a given model and assume that we are given a matrix X\in{\mathbb{X}^{\raisebox{1.39998pt}{{\gg}}}}, then the $X$-passivity radius $\rho_{\mathcal{M}}(X)$ is bounded by

[TABLE]

where $u$, $v$ and $w$ are vectors of norm $1$, satisfying

[TABLE]

*Moreover, if $v$ and $w$ are linearly dependent, then $\rho_{\mathcal{M}}(X)=\alpha\beta/2$. *

The following corollary shows how these results can be applied to pH systems.

Corollary 3.6**.**

If for a given system $\mathcal{M}$ we have that $X=I_{n}$ then the corresponding representation of the system is port-Hamiltonian, i.e., it has the representation $\mathcal{M}:=\{J-R,G-K,G^{\mathsf{H}}+K^{\mathsf{H}},S+N\}$ and the X-passivity radius of this model is given by

[TABLE]

Proof 3.7**.**

*The proof follows directly from Theorem 3.5, since under the given assumption we have $\alpha=\beta$ and we can choose $u=w$. *

Finally, we also show that pH realizations always have a better $X$-passivity radius than models that are not in pH form.

Theorem 3.8**.**

*Let $\mathcal{M}=\{A,B,C,D\}$ be a given model and let X\in{\mathbb{X}^{\raisebox{1.39998pt}{{\gg}}}}, then the port-Hamiltonian model $\mathcal{M}_{T}=\{TAT^{-1},TB,CT^{-1},D\}$ constructed from any matrix $T$ such that $X=T^{\mathsf{H}}T$, has an $I$-passivity radius $\rho_{\mathcal{M}_{T}}(I)$ which is at least as large as $\rho_{\mathcal{M}}(X)$. *

Proof 3.9**.**

It follows from Corollary 3.4 that $\rho_{\mathcal{M}}(X)$ satisfies

[TABLE]

while $\rho_{\mathcal{M}_{T}}(I)$ satisfies

[TABLE]

But the matrix inequality

[TABLE]

implies that

[TABLE]

*for all values of $\gamma$, and therefore $\rho_{\mathcal{M}}^{-1}(X)\geq\rho_{\mathcal{M}_{T}}^{-1}(I)$ or $\rho_{\mathcal{M}}(X)\leq\rho_{\mathcal{M}_{T}}(I)$. Note also that any other factorization $X=(UT)^{\mathsf{H}}(UT)$ yields the same result, since $\rho_{\mathcal{M}_{(UT)}}(I)=\rho_{\mathcal{M}_{T}}(I)$. *

4 Maximizing the passivity radius

The main goal of our paper is the maximization of the passivity radius over all pH representations of a passive system. For this we now have a closer look at the constrained LMI

[TABLE]

and obtain the following theorem.

Theorem 4.1**.**

*Let $\mathcal{M}:=\{A,B,C,D\}$ be a minimal realization of a passive system, and let $X$ be any matrix in {\mathbb{X}^{\raisebox{1.39998pt}{{>}}}}. Then there is a unique $\xi^{*}(X)$ which is maximal for the matrix inequality (41). Moreover, this value of $\xi^{*}(X)$ is also the passivity radius of the pH system $\mathcal{M}_{T}=\{TAT^{-1},TB,CT^{-1},D\}$, where $X=T^{\mathsf{H}}T$, in both the spectral and Frobenius norm. *

Proof 4.2**.**

Every X\in{\mathbb{X}^{\raisebox{1.39998pt}{{>}}}} is strictly positive definite, and can thus be factorized as $X=T^{\mathsf{H}}T$ with $\det T\neq 0$. Thus, we can define the transformed system $\mathcal{M}_{T}=\{TAT^{-1},TB,CT^{-1},D\}$. It is obvious that the passivity LMI $W(X,\mathcal{M}_{T})\geq 0$ of the transformed system $\mathcal{M}_{T}$ is satisfied with $X=I_{n}$ and that it is related to the passivity LMI $W(X,\mathcal{M})\geq 0$ of $\mathcal{M}$ via

[TABLE]

It also follows that (41) is satisfied if and only if $W(I,\mathcal{M}_{T})\geq\xi I_{n+m}$ is satisfied. But the largest value $\xi^{*}(X)$ of $\xi$ for which this holds is clearly equal to

[TABLE]

Since state-space transformations do not change the transfer function, it follows that $\mathcal{M}_{T}$ is a particular pH realization of the transfer function of $\mathcal{M}$ and that

[TABLE]

*as in Definition 2.1. Then, it follows from Corollary 3.6 that $\xi^{*}(X)$ is also equal to the passivity radius $\rho_{\mathcal{M}_{T}}(I)$ of the port-Hamiltonian system $\mathcal{M}_{T}:=\{J-R,G-K,G^{\mathsf{H}}+K^{\mathsf{H}},S+N\}$, and this for both the spectral and Frobenius norm. *

We point out that Theorem 4.1 applies to all matrices in {\mathbb{X}^{\raisebox{1.39998pt}{{>}}}}, and therefore also to all matrices in {\mathbb{X}^{\raisebox{1.39998pt}{{\gg}}}}, which can be distinguished as follows.

Corollary 4.3**.**

*The maximal value $\xi^{*}(X)$ of a matrix X\in{\mathbb{X}^{\raisebox{1.39998pt}{{>}}}} for a given model $\mathcal{M}$ equals zero if $X$ is a boundary point of {\mathbb{X}^{\raisebox{1.39998pt}{{>}}}} and is strictly positive if and only if $X$ is in {\mathbb{X}^{\raisebox{1.39998pt}{{\gg}}}}. *

Proof 4.4**.**

*If $X$ is a boundary point of {\mathbb{X}^{\raisebox{1.39998pt}{{>}}}} then $\det W(X,\mathcal{M})=0$ and for those $X$, we have $\xi^{*}(X)=0$. If $X$ belongs to {\mathbb{X}^{\raisebox{1.39998pt}{{\gg}}}}, then $W(X,\mathcal{M})>0$ and $\operatorname{diag}(X,I_{m})>0$. Therefore there exists a $\xi>0$ such that $W(X,\mathcal{M})>\xi\operatorname{diag}(X,I_{m})$, and hence $\xi^{*}(X)>0$. Conversely, if $\xi^{*}(X)>0$ then $W(X,\mathcal{M})>0$ and X\in{\mathbb{X}^{\raisebox{1.39998pt}{{\gg}}}}. *

In order to maximize the passivity radius, it is clear that we need to look at {\mathbb{X}^{\raisebox{1.39998pt}{{\gg}}}}. For a given $X$ in {\mathbb{X}^{\raisebox{1.39998pt}{{\gg}}}}, we therefore consider the passivity LMI $W(X,\mathcal{M}_{\xi})$ for the modified model $\mathcal{M}_{\xi}:=\{A+\frac{\xi}{2}I_{n},B,C,D-\frac{\xi}{2}I_{m}\}$ with a $\xi$ chosen such that

[TABLE]

We have the following Lemma.

Lemma 4.5**.**

*For every $X>0$ in {\mathbb{X}^{\raisebox{1.39998pt}{{\gg}}}} and any $0\leq\xi_{-}<\xi_{+}\leq\xi^{*}(X)$, the systems $\mathcal{M}_{\xi_{-}}$ and $\mathcal{M}_{\xi_{+}}$ are passive. Moreover, the whole solution set of $W(X,\mathcal{M}_{\xi_{+}})\geq 0$ is included in the solution set of $W(X,\mathcal{M}_{\xi_{-}})>0$. *

Proof 4.6**.**

The LMIs for two different values $\xi_{-}<\xi_{+}$ are related as

[TABLE]

Since X\in{\mathbb{X}^{\raisebox{1.39998pt}{{\gg}}}}, $\xi^{*}(X)>0$ and $\operatorname{diag}(X,I_{m})>0$, it follows that

[TABLE]

*Hence, the systems $\mathcal{M}_{\xi_{-}}$ and $\mathcal{M}_{\xi_{+}}$ are passive, since the associated LMIs have a nonempty solution set. Now consider any $X$ for which $W(X,\mathcal{M}_{\xi_{+}})\geq 0$. Since $\xi_{+}$ is strictly positive, so is $\xi^{*}(X)$ and hence X\in{\mathbb{X}^{\raisebox{1.39998pt}{{\gg}}}}. It then follows from (43) that $W(X,\xi_{-})>0$. Hence, the solution set of $W(X,\mathcal{M}_{\xi_{+}})\geq 0$ is included in the solution set of $W(X,\mathcal{M}_{\xi_{-}})>0$. *

Lemma 4.5 implies that for a given X\in{\mathbb{X}^{\raisebox{1.39998pt}{{\gg}}}}, the solution sets of $W(X,\mathcal{M}_{\xi})\geq 0$ are shrinking with increasing $\xi$. It remains to find the $\xi^{*}(X)$ that corresponds to the largest possible passivity radius. We can obtain this value by relating it to the passivity of the transfer function

[TABLE]

of the modified system $\mathcal{M}_{\xi}$. Note that we have assumed that the associated system is minimal, a property which is not changed by the shift. It follows from the discussion of Section 2 that this transfer function corresponds to a strictly passive system if and only if (i) the transfer function $\mathcal{T}_{\xi}(s)$ is asymptotically stable and (ii) the matrix function $\Phi_{\xi}(s):=\mathcal{T}^{\mathsf{H}}_{\xi}(-s)+\mathcal{T}_{\xi}(s)$ is strictly positive on the $\imath\omega$ axis, with $\omega=\infty$ included. It has been presented in Section 2 that the zeros of $\Phi_{\xi}(s)$ are also the eigenvalues of the Hamiltonian matrix

[TABLE]

provided that $D^{\mathsf{H}}+D-\xi I_{m}>0$ and the realization of $\mathcal{M}_{\xi}$ is minimal.

The three algebraic conditions corresponding to strict passivity of $\mathcal{T}_{\xi}(s)$ are, therefore, given by

1. A1.

$A+\xi I_{n}/2$ has all its eigenvalues in the open left half plane (stability). 2. A2.

$D^{\mathsf{H}}+D-\xi I_{m}$ has strictly positive eigenvalues (positive-realness at $\omega=\infty$). 3. A3.

$H_{\xi}$ has no eigenvalues on the $\imath\omega$ axis (positive-realness at finite $\omega$).

All of these conditions are phrased in terms of eigenvalues of certain matrices that depend on the parameter $\xi$. Since eigenvalues are continuous functions of the matrix elements, one can consider limiting cases for the above conditions. As explained in Section 2, the passive transfer functions are limiting cases of the strictly passive ones. These limiting cases correspond to the first value of $\xi$ where one of the three algebraic conditions fails. Note that condition A3. is more robustly expressed in terms of the eigenvalues of the matrix pencil

[TABLE]

It is obvious that the conditions A1.-A3. are not satisfied anymore for large enough $\xi$. For instance, for $\xi>\lambda_{\min}(D^{\mathsf{H}}+D)$ the second condition fails and $\lambda_{\min}(D^{\mathsf{H}}+D)$ is thus a simple upper bound for $\xi^{*}(X)$ for any $X$.

Theorem 4.7**.**

Let $\mathcal{M}$ be a strictly passive system. Then there is a bounded supremum $\Xi:=\sup_{\xi}\{\xi\;|\;\mathcal{T}_{\xi}(s)\mathrm{\;is\;strictly\;passive}\}$, for which the following properties hold

$\mathcal{T}_{\Xi}(s)$* is passive,* 2. 2.

the solution set of $W(X,\mathcal{M}_{\Xi})\geq 0$ is not empty, 3. 3.

the solution of $W(X,\mathcal{M}_{\Xi})>0$ is empty, 4. 4.

for any $\xi<\Xi$ the solution set of $W(X,\mathcal{M}_{\xi})>0$ is non-empty, 5. 5.

$\Xi:=\sup_{X}\xi^{*}(X)$* for all X\in{\mathbb{X}^{\raisebox{1.39998pt}{{>}}}}.*

Proof 4.8**.**

*The existence of a bounded supremum follows from the fact that $\mathcal{T}_{\xi}(s)$ is strictly passive only if $\xi$ is smaller than $\lambda_{\min}(D^{\mathsf{H}}+D)$. Property 1. holds because $\mathcal{T}_{\Xi}(s)$ is the limit of $\mathcal{T}_{\xi}(s)$ for $\xi\rightarrow\Xi$. Property 2. is a direct consequence of the previous property. Property 3. follows by contradiction, since if $W(X,\mathcal{M}_{\Xi})>0$ would not be empty, then $\xi^{*}(X)$ for $X$ in the solution set of $W(X,\mathcal{M}_{\Xi})>0$, would be larger than $\Xi$. Property 4. follows from Lemma 4.5, where we use any $X$ in the solution set of $W(X,\mathcal{M}_{\Xi})\geq 0$ and choose $\xi_{+}=(\Xi+\xi)/2$ and $\xi_{-}=\xi$ to show that $X$ also lies in the solution set of $W(X,\mathcal{M}_{\xi})>0$. Property 5. follows from $\xi^{*}(X)=\max\{\xi\;|\;W(X,\mathcal{M}_{\xi})\geq 0\}$, which expresses that $\mathcal{T}_{\xi}(s)$ is passive. *

We now link the value of $\Xi$ in Theorem 4.7 to the passivity radius of an optimally robust pH realization.

Theorem 4.9**.**

Let $\mathcal{M}:=\{A,B,C,D\}$ be a minimal realization of a strictly passive transfer function $\mathcal{T}(s):=C(sI-A)^{-1}B+D$. Then

[TABLE]

*is the largest possible passivity radius out of all realizations of this transfer function. The models with such an optimal passivity radius correspond to a solution $X$ of $W(X,\mathcal{M}_{\Xi})\geq 0$ and a pH realization is given by $\mathcal{M}_{T}:=\{T^{-1}AT,BT,T^{-1}C,D\}$ where $X:=T^{\mathsf{H}}T$. *

Proof 4.10**.**

Consider realizations $\mathcal{M}_{T}:=\{T^{-1}AT,BT,T^{-1}C,D\}$ with $X:=T^{\mathsf{H}}T$ and $X\in W(X,\mathcal{M})\geq 0$. It was shown in Theorem 4.1 that such realizations have passivity radius equal to $\xi^{*}(X)$ and Theorem 4.7 shows that the supremum of all $\xi^{*}(X)$ is precisely $\Xi$.

*Let us now consider an arbitrary model $\mathcal{M}$. Then its passivity radius is $\rho_{\mathcal{M}}=\rho_{\mathcal{M}}(X)$ for some X\in{\mathbb{X}^{\raisebox{1.39998pt}{{\gg}}}}. It follows that the passivity radius of the pH realization $\mathcal{M}_{T}$ derived from $X=T^{\mathsf{H}}T$ is larger than or equal to that of $\mathcal{M}$. Moreover, the corresponding passivity radius is $\xi^{*}(X)$. To complete the proof we point out that the matrices $X$ that maximize $\xi^{*}(X)$ are in $W(X,\mathcal{M}_{\Xi})\geq 0$. *

In this section we have derived a characterization of the passivity radius of a strictly passive system. In the next section we show how this can be computed numerically.

5 Computing the optimal passivity radius

In this section we describe algorithms that compute, within a given tolerance $\tau$, an approximation of the optimal $\Xi$ as in Theorem 4.9 for a given minimal realization $\mathcal{M}:=\{A,B,C,D\}$ of a passive system.

First of all, if $\mathcal{M}$ is passive but not strictly passive then $\Xi=0$. If $\mathcal{M}$ is strictly passive, then simple upper bounds for $\Xi$ are given by the conditions A1. and A2. in Section 4, i.e.,

[TABLE]

The procedure to check passivity is then to verify condition A3. for $0\leq\xi\leq\Xi_{up}$ , namely that $S_{\xi}$, (or $H_{\xi}$) has no generalized eigenvalues on the imaginary axis. We therefore first recall some basic properties of the scalar function

[TABLE]

Theorem 5.1**.**

*The real function $\gamma(\xi,\omega):=\lambda_{\min}\Phi_{\xi}(\imath\omega)$ in (45) is continuous in the real variables $\xi$ and $\omega$ and it has the following properties. It is positive for all $\omega$ in the interval $\xi\in[0,\Xi)$, it is zero for at least one value of $\omega$ at $\xi=\Xi$, and it is negative in some open interval(s) of $\omega$ whenever $\xi\in(\Xi,\Xi_{up}]$, provided that $\Xi<\Xi_{up}$. *

Proof 5.2**.**

*The continuity follows trivially from the fact that eigenvalues of a (Hermitian) matrix are continuous functions of the parameters of the matrix. It is clear that if $\xi<\Xi$ then the transfer function $\Phi(\imath\omega)$ is positive definite for all $\omega$ and so is its smallest eigenvalue $\gamma(\xi,\omega)$. When we increase $\xi$ and reach the limiting value $\Xi$, then the transfer function is passive but not strictly passive anymore, and hence $\gamma(\xi,\omega)$ must loose strict positivity in at least one point $\imath\omega$, i. e., $S_{\xi}$ must must have at least one eigenvalue on the imaginary axis. When we further increase $\xi$, $S_{\xi}$ will have more purely imaginary eigenvalues, and there will be purely imaginary eigenvalues all the way up to $\Xi_{up}$ because at $\omega=\pm\infty$ the function eigenvalues of $\Phi_{\xi}(\pm\infty)$ are those of $D^{\mathsf{H}}+D-\xi I_{m}$. By continuity, $\gamma(\xi,\omega)=0$ must therefore intersect the zero-level for all $\xi\in[\Xi,\Xi_{up}]$. Notice also that when $D^{\mathsf{H}}+D-\xi I_{m}>0$, then the number of generalized eigenvalues on the imaginary axis is bounded by $2n$, since the pencil $S_{\xi}$ is then regular and has $m$ infinite generalized eigenvalues. These three different cases are also depicted in Figure 1. *

As a consequence of Theorem 5.1, the smallest value of $\xi$ in the interval $[0,\Xi_{up}]$, where condition A3. fails is equal to $\Xi$. (Note that this could be equal to $\Xi_{up}$.) One can then apply bisection to this interval and check the presence of purely imaginary eigenvalues in the above interval. Putting $\Xi_{lo}=0$, we have the following procedure.

Bisection procedure for computing $\Xi$:

[TABLE]

Since the interval for $\Xi$ shrinks by a factor $2$ at each step of the iteration, then in $k=\lceil\log_{2}(\Xi_{up}/\tau)\rceil$ steps we will have $\Xi_{up}-\Xi_{lo}\leq\tau$.

One can also make use of the computed eigenvalue decompositions of

[TABLE]

to construct a method with faster convergence properties. This approach relies on the fact that $\gamma(\xi,\omega)$ is a continuous function with smoothness properties at the local minima. It is inspired from a method in [7] for the computation of the $L_{\infty}$-norm of a transfer function.

Let $\hat{\xi}\in[\Xi,\Xi_{up}]$, let $\Omega$ be the set of purely imaginary eigenvalues of $S_{\hat{\xi}}$, let $\partial\Omega$ be the corresponding set of derivatives of these eigenvalues with respect to $\omega$ , and let $\gamma(\hat{\xi},\omega)$ be the smallest eigenvalue of $\Phi_{\hat{\xi}}(\imath\omega)$. Then we exploit the following ideas.

The real roots $\omega_{i}$ of $\gamma(\hat{\xi},\omega)$ are a subset of $\Omega$. They can be identified by looking at the derivatives in $\partial\Omega$. A detailed description of this selection procedure is described in [26]. 2. 2.

For a given value of $\hat{\xi}$, $\gamma(\hat{\xi},\omega)$ we can locally approximate $\gamma$ in the neighborhood of a local minimum $\omega^{*}$ by a quadratic approximation $\gamma^{*}+(\omega-\omega^{*})^{2}$. A detailed description of this result is described in [7]. 3. 3.

If $\omega_{1}<\omega_{2}$ are two “consecutive” zeros of $\gamma(\hat{\xi},\omega)$ then the smallest real root $\tilde{\xi}$ of $S_{\xi}(\hat{\omega})$ lies between [math] and $\hat{\xi}$ in that interval. This follows from the fact that the smallest eigenvalue $\gamma(\hat{\xi},\omega)$ for the current value $\hat{\xi}$ is negative in that interval of $\omega$. Moreover, at the midpoint $\hat{\omega}:=(\omega_{1}+\omega_{2})/2$ the smallest real root $\tilde{\xi}$ is an improved upper bound for $\Xi$, in the sense that close to $\Xi$, the corresponding error improves quadratically, i.e., $|\Xi-\tilde{\xi}|\approx|\Xi-\hat{\xi}|^{2}$. A detailed description of this result follows from [7] and the continuity of the two-variable function $\gamma(\xi,\omega)$.

Since we can compute the zeros $\omega_{i}$ of $\gamma(\hat{\xi},\omega)$) for a given value of $\hat{\xi}$, and we can find the smallest real root $\tilde{\xi}$ of $S_{\xi}(\imath\hat{\omega})$ for a given value of $\hat{\omega}$, we can then obtain an algorithm with improved local convergence for computing $\Xi$ based on the above ideas. A detailed analysis of its convergence requires more investigation and will be considered in forthcoming work.

An improved algorithm:

$\hat{\xi}:=\Xi_{up}-\tau$. 2. 2.

Compute the eigenvalues of $S_{\hat{\xi}}(\imath\omega)$ and select those corresponding to real zeros of $\gamma(\hat{\xi},\omega)$. 3. 3.

if $\gamma(\hat{\xi},\omega)$ has no real zeros, then $\Xi_{lo}=\hat{\xi}$, stop, 4. 4.

else let $\hat{\omega}:=(\omega_{1}+\omega_{2})/2$ for the largest interval $[\omega_{1},\omega_{2}]$ of these roots.

       Compute the real eigenvalues $\xi_{i}$ of $S_{\xi}(\imath\hat{\omega})$ and update $\Xi_{up}:=\min_{i}\xi_{i}$.

       Define the next guess $\hat{\xi}:=\Xi_{up}-\tau$ and go to step 2.

This algorithm will typically require only a few iterations and stops with an interval $[\Xi_{lo},\Xi_{up}]$ of size $\tau$. Each step of both, the bisection and the faster converging algorithms has a complexity that is cubic in the matrix dimensions, but they are guaranteed to provide a required accuracy. For large scale problems, this complexity may become a problem, but one can combine it with iterative structured Krylov space techniques like those discussed in [19] or model reduction based techniques as in [2, 3].

6 The distance to passivity

In real world applications the system data are usually inaccurate and therefore only an approximate model is known. Often the real physical problem is nonlinear or a partial differential equation and it is approximated by a finite element or finite difference model [17], it may be obtained by a realization or system identification procedure [10, 16, 23], or the result of a model reduction procedure [4]. Then it is in general not clear that in the approximation process passivity is preserved. In this situation one approximates the non-passive system by a (hopefully) nearby passive system, by introducing small or minimal perturbations to the model $\mathcal{M}:=\{A,B,C,D\}$, see [1, 9, 10, 12, 15, 23, 25]. Therefore, in this section we start with a system $\mathcal{M}$ that is not passive and ask the question about the smallest perturbation $\Delta_{\mathcal{M}}$ of the model that makes the system $\mathcal{M}+\Delta_{\mathcal{M}}$ passive.

It is clear from our previous analysis that to find the smallest perturbation that makes the system passive is equivalent to finding the smallest perturbation $\Delta_{\mathcal{M}}$ such that the LMI $W(X,\mathcal{M}+\Delta M)\geq 0$ has a Hermitian and positive semidefinite solution $X$. Moreover, if the perturbed system remains minimal, then we expect $X>0$. We recall the notations

[TABLE]

Definition 6.1**.**

The distance to passivity of a minimal system $\mathcal{M}\!=\!\{A,B,C,D\}$ is the minimum norm $\|\Delta_{\mathcal{S}}\|_{2}$ ($\|\Delta_{\mathcal{S}}\|_{F}$) such that there exists a matrix $X>0$ satisfying

[TABLE]

In the following we need an extension of Lemma 4.5, for which we consider the LMI for the modified model $\mathcal{M}_{-\xi}:=\{A-\frac{\xi}{2}I_{n},B,C,D+\frac{\xi}{2}I_{m}\}$ with the corresponding transfer function

[TABLE]

and LMI

[TABLE]

Lemma 6.2**.**

*Let $\mathcal{M}:=\{A,B,C,D\}$ be a non-passive system. Then for every $X>0$ in $\mathbb{H}_{n}$ there exists a $\xi^{*}(X)>0$ such that the LMI (51) for the system $\mathcal{M}_{-\xi^{*}(X)}$ holds. Moreover, for every value $\xi>\xi^{*}(X)$, the system $\mathcal{M}_{-\xi}$ is passive. *

Proof 6.3**.**

*Clearly we have that $W(X,\mathcal{M}_{-\xi})=W(X,\mathcal{M})+\xi\hat{X}$. Since $W(X,\mathcal{M})$ is bounded from below and $\hat{X}>0$, the inequality $W(X,\mathcal{M}_{-\xi})\geq 0$ holds for a sufficiently large value of $\xi$. Let $\xi^{*}(X)$ be the smallest $\xi$ value for which $W(X,\mathcal{M}_{-\xi})\geq 0$ holds, then $W(X,\mathcal{M}_{-\xi})=W(X,\mathcal{M}_{-\xi^{*}(X)})+(\xi-\xi^{*}(X))\hat{X}$, which implies that the passivity condition holds for all $\xi>\xi^{*}(X)$. *

To find the optimal $\xi$, let us consider first perturbations $\Delta_{\mathcal{S}}$ that are a multiple of the identity.

Theorem 6.4**.**

The minimum norm perturbation of the type

[TABLE]

*that makes system $\mathcal{M}$ passive, has spectral norm $\Xi/2$ and Frobenius norm $\Xi\sqrt{n+m}/2$, where $\Xi$ is the minimum value of $\xi$ such that the model $\mathcal{M}_{-\xi}:=\{A-\xi I_{n}/2,B,C,D+\xi I_{m}/2\}$ with transfer function $T_{-\xi}(s)$ is passive. *

Proof 6.5**.**

It follows from (46) that $\xi$ must satisfy

[TABLE]

*for some $X>0$. By Lemma 6.2 there exists a bounded minimal solution, which we call $\Xi$. The state-space model corresponding to $\mathcal{S}+\frac{1}{2}\xi I_{n+m}$ is $\mathcal{M}_{-\xi}$ with transfer function (47). Therefore, $\Xi$ is the smallest value of $\xi$ that makes the model $\mathcal{M}_{-\xi}$ with transfer function $T_{-\xi}(s)$ become passive. We can then choose $X>0$ from the solution set of $W(X,\mathcal{M}_{-\Xi})\geq 0$ to satisfy (46). *

The minimal value $\Xi$ of $\xi$ for the restricted class of perturbations being a multiple of the identity can then be computed with the algorithms described in the previous section.

Since, most likely, the perturbation will not lead to (46) being an equality, we can reduce the Frobenius norm of the perturbation $\Delta_{\mathcal{S}}$ by considering more general perturbations than (52). In order to do that, we use a matrix $X$ from the set $W(X,\mathcal{M}_{-\Xi})\geq 0$, where $\Xi$ was obtained from perturbations as in (52). We use the factorization $X=T^{\mathsf{H}}T$, to transform the system $\mathcal{M}$ to the pH form $\mathcal{M}_{T}=\{TAT^{-1},TB,CT^{-1},D\}:=\{J-R,G-K,G^{\mathsf{H}}+K^{\mathsf{H}},S+N\}$. If we denote the transformed matrix $\mathcal{S}$ as $\hat{\mathcal{S}}:=\hat{T}\mathcal{S}\hat{T}^{-1}$, where $\hat{T}=\operatorname{diag}(T,I_{m})$, then it follows that

[TABLE]

since $\mathcal{M}_{-\Xi}$ is passive. The smallest eigenvalue of $\frac{1}{2}(\hat{\mathcal{S}}^{\mathsf{H}}+\hat{\mathcal{S}})$ is greater or equal than $-\Xi$, but the other eigenvalues may be larger or even positive.

To obtain a solution for the minimum Frobenius norm problem

[TABLE]

let $\hat{T}=\hat{U}^{T}\hat{\Sigma}\hat{V}$ be a singular value decomposition with $\hat{\Sigma}$ diagonal and $\hat{U},\hat{V}$ unitary. Then we can alternatively study the problem

[TABLE]

where $\Delta_{\tilde{\mathcal{S}}}:=\hat{V}\Delta_{\mathcal{S}}\hat{V}^{H}$ and $\tilde{\mathcal{R}}:=\hat{U}\hat{\mathcal{R}}\hat{U}^{\mathsf{H}}$.

For the given $\tilde{\mathcal{R}}$ we have a unitary spectral decomposition

[TABLE]

where $D_{1}$ is diagonal with positive diagonal elements and $D_{2}$ is diagonal with nonnegative diagonal elements. Then $\Delta\tilde{\mathcal{R}}=\tilde{U}^{\mathsf{H}}\operatorname{diag}(0,D_{2})\tilde{U}$ is the minimum Frobenius norm perturbation so that $\tilde{\mathcal{R}}+\Delta\tilde{\mathcal{R}}$ is positive semidefinite. We can then replace the optimization problem (54) by the simpler problem

[TABLE]

Indeed, let $\Delta F$ be the strictly lower triangular part plus $\frac{1}{2}$ the diagonal part of $\Delta\tilde{\mathcal{R}}$. Then clearly $\hat{\Sigma}^{-1}\Delta F\hat{\Sigma}$ is the optimal Frobenius norm solution to (55).

It should be pointed out that we could have used another matrix $X$ from the domain of $W(X,\mathcal{M}_{-\Xi})\geq 0$. The matrices $R$, $K$ and $S$ will then change and the minimum Frobenius norm solution will be affected as well. Note that if we choose a matrix $X$ that does not belong to the solution set of $W(X,\mathcal{M}_{-\Xi})\geq 0$, then because of Lemma 6.2, the value of the most negative eigenvalue in $\Lambda$ will be equal to $-\xi^{*}(X)$, but it follows that $\xi^{*}(X)>\Xi$. It is therefore unlikely that the Frobenius norm of the constructed minimum norm perturbation will be smaller.

7 The distance to stability

We can employ similar arguments as in the previous section for the computation of the smallest perturbation that makes a given system stable, see [13, 14, 20] for other approaches. When searching for the smallest perturbation we only have to study the $(1,1)$ block of the LMI (4), i.e. the case where the matrices $B$, $C$ and $D$ are void, or where $m=0$. In this case we write the LMI (4) as

[TABLE]

while the shifted LMI takes the form

[TABLE]

For any $X$ in the solution set of $W(X,A)\geq 0$, the shifted LMI (57) has a nonempty solution set as long as $\xi$ is smaller than $\xi^{*}(X)$. Instead of the three conditions A1.–A3. for passivity, we only have one condition, namely that the Hamiltonian matrix

[TABLE]

has no purely imaginary eigenvalue for all $0\leq\xi<\xi^{*}(X)$. So if $\Xi$ is the largest value of $\xi^{*}(X)$ over all $X$ in the solution set of $W(X,A)\geq 0$, then the solution set of $W_{\Xi}(X,A)\geq 0$ is not empty, but it has an empty interior.

A well-known formula for the stability radius $\rho_{A}$ of a given matrix $A$ is given by

[TABLE]

The perturbation that achieves the minimum is constructed from the smallest singular value $\sigma_{\min}$ and corresponding singular vectors $u$ and $v$ of $A-\imath\omega I_{n}$ at the minimizing value of $\omega$ :

[TABLE]

A similar result is obtained using the analysis of Section 3. This encourages us to also use the distance to passivity results to derive bounds for the distance to stability.

Let us first look at the solution of the problem for a specific diagonal form of the perturbation. Given an unstable matrix $A$ we study the question what is the smallest perturbation of the type $\Delta_{A}:=-\frac{1}{2}\xi I_{n}$ that makes the matrix stable. Following a similar analysis as in Section 6 we can prove the following theorem.

Theorem 7.1**.**

The minimum norm solution of the form

[TABLE]

*for the stabilization problem of a matrix $A$, has spectral norm $\Xi/2$ and Frobenius norm $\Xi\sqrt{n}/2$, where $\Xi$ is the minimum value of $\xi$ such that the matrix $A_{-\xi}:=A-\xi I_{n}/2$, is stable. *

Proof 7.2**.**

Clearly $\xi$ must satisfy

[TABLE]

for some $X>0$. It follows from Lemma 6.2 that there exists a bounded minimal solution, which we call $\Xi$. The matrix corresponding to $A_{-\xi}:=A-\frac{1}{2}\xi I_{n}$ has the associated Hamiltonian matrix

[TABLE]

*Let $\Xi$ be the smallest value such that $H_{-\xi}$ has no purely imaginary eigenvalues for $\xi<\Xi$. Then the matrix $A_{-\Xi}$ is stable (but not asymptotically stable) and we can then choose $X>0$ from the solution set of $W(X,A_{-\Xi})\geq 0$ to satisfy (60). *

Since we made some of the eigenvalues of the LMI in (60) strictly positive, rather than non-negative, we can further reduce the Frobenius norm of the perturbation $\Delta_{A}$ when removing the imposed restriction of a diagonal perturbation. In order to do that, we use a matrix $X$ from the set $W(X,A_{-\Xi})\geq 0$, where $\Xi$ was obtained from the diagonal perturbation. Using its factorization $X=T^{\mathsf{H}}T$, we transform $A$ to $A_{T}=TAT^{-1}:=J-R$, which satisfies

[TABLE]

since $A_{-\Xi}$ is stable. We know that the smallest eigenvalue equals $-\Xi$, but the others may be larger or even positive. In order to construct a nearly optimal solution to this, we can again use the two stage procedure and the perturbation result derived in the previous section.

8 Conclusion

We have presented analytic formulas and numerical methods to construct optimally robust port-Hamiltonian realizations of a given transfer function of a linear time invariant passive system. We have shown how to use shifted linear matrix inequalities to achieve this goal. The techniques can also be applied to compute a nearby passive system to a given non-passive one or a nearby stable system to a given unstable one.