Document [original]

Port-Hamiltonian Realizations

of Linear Time Invariant Systems

Christopher Beattie∗and Volker Mehrmann†and Hongguo Xu‡

January 7, 2016

Abstract

The question when a general linear time invariant control system is equivalent to a

port-Hamiltonian systems is answered. Several equivalent characterizations are derived

which extend the characterizations of [38] to the general non-minimal case. An explicit

construction of the transformation matrices is presented. The methods are applied in the

stability analysis of disc brake squeal.

Keywords: port-Hamiltonian system, passivity, stability, system transformation, linear ma-

trix inequality, Lyapunov inequality, even pencil, quadratic eigenvalue problem.

AMS subject classification. 93A30, 93B17, 93B11.

1 Introduction

The synthesis of system models that describe complex physical phenomena often involves

the coupling of independently developed subsystems originating within different disciplines.

Systematic approaches to coupling such diversely generated subsystems prudently follows a

system-theoretic network paradigm that focusses on the transfer of energy, mass, and other

conserved quantities among the subsystems. When the subsystem models themselves arise

from variational principles, then the aggregate system typically has structural features that

reflects underlying conservation laws and often it may be characterized as a port-Hamiltonian

(PH) system, see [11, 15, 27, 29, 30, 33, 32, 34, 35, 36] for some major references. Although

PH systems may be formulated in a more general framework, we will restrict ourselves to

input-state-output PH systems, which have the form

x= (J−R)∇xH(x) + (F−P)u(t),(1)

y(t) = (F+P)T∇xH(x) + (S+N)u(t),

∗Department of Mathematics, Virginia Tech, Blacksburg, VA 24061, USA. [email protected]. Supported by

Einstein Foundation Berlin, through an Einstein Visiting Fellowship.

2Institut f¨ur Mathematik MA 4-5, TU Berlin, Str. des 17. Juni 136, D-10623 Berlin, FRG.

[email protected]. The authors gratefully acknowledge the support by the Deutsche Forschungs-

gemeinschaft (DFG) as part of the collaborative research center SFB 1029 Sub- stantial efficiency increase in

gas turbines through direct use of coupled unsteady combustion and flow dynamics, project A02 Development

of a reduced order model of pulsed detonation combuster and by Einstein Center ECMath in Berlin.

3Department of Mathematics, University of Kansas, Lawrence, KS 66045, USA. [email protected]. Partially

supported by Alexander von Humboldt Foundation and by Deutsche Forschungsgemeinschaft, through the DFG

Research Center Matheon Mathematics for Key Technologies

where x∈Rnis the n-dimensional state vector;H:Rn→[0,∞) is the Hamiltonian, a

continuously differentiable scalar-valued vector function describing the distribution of internal

energy among the energy storage elements of the system; J=−JT∈Rn×nis the structure

matrix describing the energy flux among energy storage elements within the system; R=

RT∈Rn×nis the dissipation matrix describing energy dissipation/loss in the system; F±P∈

Rn×mare the port matrices, describing the manner in which energy enters and exits the

system, and S+N, with S=ST∈Rm,m and N=−NT∈Rm,m, describing the direct

feed-through of input to output. The matrices, R,P, and Smust satisfy

K=[R P

PTS]≥0; (2)

that is, Kis symmetric positive-semidefinite. This implies, in particular, that Rand Sare

also positive semidefinite, R≥0 and S≥0.

Port-Hamiltonian systems generalize the classical notion of Hamiltonian systems expressed

in our notation as ˙

x=J∇xH(x). The analog of the conservation of energy for Hamiltonian

systems is for PH systems (1) the dissipation inequality:

H(x(t1)) −H(x(t0)) ≤∫t1

y(t)Tu(t)dt, (3)

which has a natural interpretation as an assertion that the increase in internal energy of the

system, as measured by H, cannot exceed the total work done on the system. H(x) is a

storage function associated with the supply rate,y(t)Tu(t). In the language of system theory,

(3) constitutes the property that (1) is a passive system [10].

One may verify with elementary manipulations that the inequality in (3) is an immediate

consequence of the inequality in (2), and holds even when the coefficient matrices J,R,F,

P,S, and Ndepend on xor explicitly on time t(see, [27]) or, indeed (with care taken to

define suitable operator domains), when they are defined as linear operators acting on infinite

dimensional spaces [21, 35]. Notice that with a null input, u(t) = 0, the the dissipation

inequality asserts that H(x) is non-increasing along any unforced system trajectory. Thus,

H(x) defines a Lyapunov function for the unforced system, so PH systems are implicitly

Lyapunov stable [19]. Similarly, H(x) is non-increasing along any system trajectory that

produces a null output, y(t) = 0, so PH systems also have Lyapunov stable zero dynamics

[9].

PH systems constitute a class of systems that is closed under power-conserving intercon-

nection. This means that port-connected PH systems produce an aggregate system that must

also be port-Hamiltonian. This aggregate system hence, will be guaranteed to be both stable

and passive. Modeling with PH systems, thus, represents physical properties in such a way

as to facilitate automated modeling [23] and at same time encoding physical properties in the

structure of equations. This framework also provides compelling motivation to identify and

preserve PH structure when producing reduced order surrogate models, see [1, 18, 31].

Now consider a general linear time-invariant (LTI) system:

x=A x +B u

y=C x +D u,(4)

with A∈Rn×n,B∈Rn×m,C∈Rm×n, and D∈Rm×m. The system (4) is passive if there

exists a continuously differentiable storage function H:Rn→[0,∞) such that (3) holds for

all admissible inputs u.

So the natural question arises: When is (4) equivalent to a port-Hamiltonian system?

Here specifically, we consider equivalence to LTI port-Hamiltonian systems taking the

form ˙

ξ= (J−R)Qξ+ (F−P)φ,

η= (F+P)TQξ+ (S+N)φ,(5)

with J=−JT,R=RT≥0, Q=QT>0, S=ST≥0, N=−NT, where J,R,Q∈Rn×n,

F,P∈Rn×m,S,N∈Rm×m, and Kas defined in (2) is positive semidefinite.

The notion of system equivalence that we consider here engages three invertible transfor-

mations connecting (4) and (5), one on each of the input, output, and state space:

u(t) = ˜

Vφ(t),η(t) = VTy(t),and x(t) = Tξ(t) (with ˜

V,V,Tinvertible).

Within this context, the supply rate associated with (4) is transformed as

y(t)Tu(t) = η(t)TV−1˜

Vφ(t).

We wish to constrain the permissible transformations characterizing system equivalence so as

to be power conserving; that is, so that supply rates remain invariant, y(t)Tu(t) = η(t)Tφ(t).

Thus, we assume that ˜

V=Vand we say that (4) is equivalent to (5) if there exist invertible

matrices, Vand T, such that

u(t) = Vφ(t),η(t) = VTy(t),and x(t) = Tξ(t).(6)

Since PH systems are structurally passive, our starting point is the following characteriza-

tion of passivity introduced in [39] for minimal linear time invariant systems. The system (4)

is minimal if it is both controllable and observable. The system (4) (and more specifically, the

pair of matrices (A,B) with A∈Rn×n,B∈Rn×m) is controllable if rank [sI−A B ]=n

for all s∈C. Similarly, the system (4) (and the pair (A,C) with A∈Rn×n,C∈Rm×nis

observable if rank [sI−A

C]=nfor all s∈C.

Theorem 1 ([39]) Assume that the LTI system (4) is minimal. The matrix inequality

[ATQ+QA QB −CT

BTQ−C−(D+DT)]≤0 (7)

has a solution Q=QT>0if and only if (4) is a passive system, in which case:

i) H(x) = 1

2xTQx defines a storage function for (4) associated with the supply rate yTu,

satisfying (3).

ii) There exist maximum and minimum symmetric solutions to (7): 0<Q−≤Q+such

that for all symmetric solutions Qto (7), Q−≤Q≤Q+.

This result yields an immediate consequence for PH realizations.

Corollary 2 Assume that the LTI system (4) is minimal. Then (4) has a PH realization if

and only if it is passive. Moreover, if (4) is passive then every equivalent system to (4) (as

generated by transformations as in (6)) may also be directly expressed as a PH system.

Proof. If (4) has a PH realization then a fortiori it is passive. Conversely, if (4) is passive

then (7) has a positive definite solution ˆ

Q=ˆ

QT=TTT, written in terms of a Cholesky or

square root factor T. Then we can define directly

Q=I,J=1

2(TAT−1−(TAT−1)T),R=−1

2(TAT−1+ (TAT−1)T)

F=1

2(TB + (CT−1)T),P=1

2(TB + (CT−1)T),

S=1

2(D+DT)N=1

2(D−DT)

(8)

and (7) can be written in terms of these defined quantities as

−2[TT0

0I][R P

PTS][T0

0I]≤0,

which verifies (2). Thus, J,R,Q,F,P,S,and Nas defined in (8) indeed determine a PH

system.

The matrix inequality (7) actually includes two properties, the Lyapunov inequality in the

upper left block,

ATQ+QA ≤0,(9)

via Lyapunov’s theorem [26], guarantees that the uncontrolled system ˙x=Axis stable, i.e., A

has all eigenvalues in the closed left half plane and those on the imaginary axis are semisimple,

while passivity is encoded in the solvability of the full matrix inequality (7). If the inequality

(9) is strict, then the system is asymptotically stable, i.e., all eigenvalues of Aare in the open

left half plane, and strictly passive, i.e., also the dissipation inequality (3) is strict, see [24]

for a detailed analysis.

We will discuss the solution of these two types matrix inequalities in the more general

situation of non-minimal systems in Section 2. Based on these results we will then construct

explicit transformations from of a general LTI system (minimal or not) to port-Hamiltonian

form (5) in Section 3.

2 Lyapunov and Riccati inequalities

The solution of Lyapunov and Riccati inequalities as they arise in the characterizations (9)

and (7) is usually addressed via the solution of semidefinite programming, see [4]. We discuss

the more explicit characterization via invariant subspaces.

2.1 Solution of Lyapunov inequalities.

The stability of Ais a necessary condition for (9) and (7) which require that TAT−1+

T−TATTT≤0, or equivalently, that the Lyapunov inequality (9) has a positive definite

solution Q=TTT.

It is well known that the equality case in (9) always has a positive definite solution if A

is stable, see [26]. In the following we recall, see e.g. [4], a characterization of the complete

set of solutions of the inequality case.

If Ais stable, but not asymptotically stable, then due to the fact that the eigenvalues on

the imaginary axis are all semi-simple, using the real Jordan form of A, see e.g. [20], there

exist a nonsingular matrix M∈Rn×nsuch that

MAM−1= diag (A1, α2J2, . . . , αrJr),(10)

where A1∈Rn1×n1is asymptotically stable, α2, . . . , αr≥0 are real and distinct, and Jj=

[0Inj

−Inj0],j= 2, . . . , r. To characterize the solution set of (9), we make the ansatz that

Q=MTdiag (Q1,ˆ

Q2, . . . , ˆ

Qr)M,(11)

and separately consider the determination of the block Q1and the other blocks. Let W(n1)

be the set of symmetric positive semidefinite matrices Θ1∈Rn1×n1with the property that

Θ1x= 0 for any eigenvector xof A1. Then for any Θ1∈ W(n1) we define Q1to be the

unique symmetric positive definite solution of the Lyapunov equation AT

1Q1+Q1A1=−Θ1,

see [26]. The other matrices ˆ

Qj,j= 2, . . . , r are chosen of the form

Qj=[YjZj

−ZjYj]>0,

with Zj=−ZT

j, when αj>0 or an arbitrary ˆ

Qj>0 when αj= 0.

We have the following characterization of the solution set of (9).

Lemma 3 Let A∈Rn×n. Then the Lyapunov inequality (9) has a symmetric positive definite

solution Q∈Rn×nif and only if Ais stable.

If Ais asymptotically stable, then the solution set is given by the set of all symmetric

positive definite solutions of the Lyapunov equation ATQ+QA =−Θ, where Θis any

symmetric positive semidefinite matrix with the property that Θx= 0 for any eigenvector x

of A, or in other words (A, Θ) is observable.

If Ais stable, but not asymptotically stable, then with the transformation (10), any solu-

tion of (9) has the form (11) solving the Lyapunov equation

ATQ+QA =−MTdiag(Θ1,0, . . . , 0)M(12)

with Θ1∈ W(n1).

Proof. Consider a transformation to the form (10) and set, for a symmetric matrix Q,

M−TQM−1= [Qij] partitioned accordingly. Form

M−T(ATQ+QA)M−1

=⎡

⎢

⎣

1Q11 +Q11A1AT

1Q12 +α2Q12J2. . . AT

1Q1,r +αrQ1,rJr

α2JT

2QT

12 +QT

12A1α2(JT

2Q22 +Q22J2). . . α2JT

2Q2,r +αrQ2,rJr

.....

αrJT

rQT

1,r +QT

1,rA1αrJT

rQT

2,r +α2QT

2,rJ2. . . αr(JT

rQr,r +Qr,rJr)

⎤

⎥

⎦,

again partitioned accordingly. For any j= 2, . . . , r with αr= 0 partition

Qjj =[Qj,11 Qj,12

j,12 Qj,22 ]

according to the block structure in Jj, so that

jQjj +QjjJj=[−Qj,12 −QT

j,12 Qj,11 −Qj,22

Qj,11 −Qj,22 Qj,12 +QT

j,12 ].

Then, αj(JT

jQjj +QjjJj)≤0 if and only if Qj,12+QT

j,12 = 0 and it follows that Qj,11−Qj,22 =

0. Thus any matrix of the form Qjj =[YjZj

−ZjYj]>0 with Zj=−ZT

j, is a positive definite

solution. If αj= 0, then clearly any Qjj >0 is a solution.

In both cases the resulting diagonal blocks Qjj >0 satisfy

αj(JT

jQjj +QjjJj) = 0.

Thus, for i=j, the blocks Qij have to satisfy one of the two matrix equations

1Q1j+αjQ1jJj= 0, αiJT

iQij +αjQijJj= 0,

which implies that Qij = 0 for all i=j, since the spectra of the respective pairs A1, αjJj

and αiJi, αjJjare different and thus the corresponding Sylvester equations only have 0 as

solution [26].

For any positive semidefinite matrix Θ1satisfying Θ1x= 0 for any eigenvector xof A1,

by Lyapunov’s Theorem [26], there exists a symmetric positive definite matrix Q1satisfying

1Q1+Q1A1=−Θ1. Thus, the Lyapunov inequality (9) has a solution Q>0 and it has

the form (11) satisfying (12). The results for an asymptotically stable Acan be easily derived

as a special case with all ˆ

Q2, . . . , ˆ

Qrvoid in (11).

Conversely, if the Lyapunov inequality has a solution Q>0, and if x= 0 is an eigenvector

of A, i.e., Ax =λx, then

(¯

λ+λ)x∗Qx ≤0.

Since x∗Qx >0, it follows that ¯

λ+λ≤0, i.e., the real part of λis non-positive. Suppose

that Ahas a purely imaginary eigenvalue iα. Since a similarity transformation of Adoes not

change the Lyapunov inequality (9), we may assume that Ais in (complex) Jordan canonical

form, and we may consider each Jordan block separately. If Akis a single real Jordan block

of size k×k,k > 1 with eigenvalue iα, then for any Hermitian matrix ˆ

Qk= [qij] we have

that

A∗

kQk+QkAk=⎡

⎢

⎣

0q11 . . . q1,k−1

q11 q12 + ¯q12 . . . q1k+q2,k−1

.....

¯q1,k−1¯q2,k−1+ ¯q1,k . . . qk−1,k + ¯qk−1,k

⎤

⎥

⎦.

This matrix cannot be negative semidefinite, since otherwise q11 = 0, which would contradict

that ˆ

Qkis positive definite.

Therefore, Ahas to be stable, i.e., if it has purely imaginary eigenvalues, these eigenvalues

must be semi-simple.

Remark 1 The construction in Lemma 3 relies on the computation of the real Jordan form

of A, which usually cannot be computed in finite precision arithmetic. For the numerical

computation of the solution it is better to use the real Schur form, see [16].

2.2 Solution of Riccati inequalities

By Corollary 2 we can characterize (at least in the minimal case) the existence of a transfor-

mation to PH form via the existence of a symmetric positive definite matrix Qsolving the

linear matrix inequality [−ATQ−QA CT−QB

C−BTQ D +DT]≥0,(13)

In this subsection we will discuss the explicit solution of this matrix inequality in the general

(non-minimal) case. It is clear that for the inequality to hold, S:= D+DThas to positive

semidefinite. However, we will see below that we can restrict ourselves to the case that S:=

D+DTis positive definite, so we only consider this case here and, using Schur complements,

(13) is equivalent to the Riccati inequality

(A−BS−1C)TQ+Q(A−BS−1C) + QBS−1BTQ+CTS−1C≤0.(14)

We first investigate the influence of the purely imaginary eigenvalues of A. Clearly, a necessary

condition for (13) to be solvable is that Qsatisfies ATQ+QA ≤0. Following Lemma 3, if

Ahas the form (10) then Qmust have the form (11), and ATQ+QA has the form (12).

Written in compact form, we get

MAM−1=[A10

0A2],M−TQM−1=[Q10

0Q2],

where

A2= diag(α2J2, . . . , αrJr) = −AT

2,Q2= diag( ˆ

Q2, . . . , ˆ

Qr),

satisfying AT

2Q2+Q2A2= 0, and

M−T(ATQ+QA)M−1=[AT

1Q1+Q1A10

0 0 ].

Setting

MB =[B1

B2],CM−1=[C1C2].

and premultiplying M−Tand post-multiplying M−1to the first block row and column of

(13), respectively, one has

⎡

⎣−AT

1Q1−Q1A10CT

1−Q1B1

0 0 CT

2−Q2B2

C1−BT

1Q1C2−BT

2Q2S⎤

⎦≥0.

Therefore, Q2must be positive definite satisfying

2Q2=C2,AT

2Q2+Q2A2= 0,(15)

and [−AT

1Q1−Q1A1CT

1−Q1B1

C1−BT

1Q1S]≥0,(16)

or equivalently it has to satisfy the reduced Riccati inequality

Ψ(Q1) := (A1−B1S−1C1)TQ1+Q1(A1−B1S−1C1)+Q1B1S−1BT

1Q1+CT

1S−1C1≤0.(17)

To study the solvability of (17), we note that A1is asymptotically stable. We claim that

A1−B1S−1C1is necessarily asymptotically stable as well. To show this, suppose that the

inequality (17) has a solution Q1>0. Since Q1B1S−1BT

1Q1+CT

1S−1C1≥0, it follows

from Lemma 3 that A1−B1S−1C1is stable, and there exits an invertible matrix M1such

that M1(A1−B1S−1C1)M−1

1has real Jordan form (10), M−T

1Q1M−1

1has the form (11) and

following (12), we have

M−T

1(Q1B1S−1BT

1Q1+CT

1S−1C1)M−1

1= diag(Θ,0, . . . , 0) ≥0.

Then, due to the positive definiteness of Sit follows that C1M−1

1=[C11 0. . . 0],

and by making use of the block diagonal structure of M−T

1Q1M−1

1, we also have M1B1=

[BT

11 0. . . 0]T. Thus it follows that

M1A1M−1

1= diag(A11, α2J2, . . . , αrJr).

Since A1is asymptotically stable, all αjJjmust be void, which implies that A1−B1S−1C1

must be asymptotically stable as well.

In order to characterize the solution of the reduced Riccati inequality (17), we use the

following lemmas. The first is the orthogonal version of the Kalman decomposition, [22, 37].

Lemma 4 Consider a general LTI system of the form (4). Then there exists a real orthogonal

matrix Usuch that

UTAU =⎡

⎣

A11 0A13

A21 A22 A23

0 0 A33 ⎤

⎦=: [˜

A11 ˜

A12

0˜

A22 ]

(18)

UTB=⎡

⎣

0⎤

⎦=: [˜

0],CU =[C10C3]=: [˜

C1˜

C2],

where the pair (˜

A11,˜

B1)is controllable and the pair (A11,C1)is observable.

Proof. Applying first the controllability and then the observability staircase forms of [37] to

the system, there exists a real orthogonal matrix Usuch that UTAU is in the form (18) with

(˜

A11,˜

B1) controllable and (A11,C1) observable.

Note that by using the block structure in Lemma 4, we have that (A11,B1) is controllable

as well.

The next lemma considers the Riccati inequality in the form (refstcf).

Lemma 5 Consider the reduced Riccati inequality (17) and suppose that both A1and A1−

B1S−1C1are asymptotically stable, with S>0. Then there exists a transformation to the

condensed form (18) of A1,B1,C1, such that

UT(A1−B1S−1C1)U=[˜

A11 −˜

B1S−1˜

C1˜

A12 −˜

B1S−1˜

0˜

A22 ]

=⎡

⎣

A11 −B1S−1C10A13 −B1S−1C3

A21 −B2S−1C1A22 A23 −B2S−1C3

0 0 A33 ⎤

⎦,(19)

where A11 −B1S−1C1,A22, and A33 are all asymptotically stable. In addition, also A11 is

asymptotically stable, (˜

A11,˜

B1)and (A11,B1)are controllable, and (A11,C1)is observable.

The proof is straightforward.

The following Lemma is also well-known, see e.g. [25, 28], but we present a proof in our

notation.

Lemma 6 Suppose that (A,B)is controllable, (A,C)is observable, Ais asymptotically

stable, and S>0. Then the Riccati equation

(A−BS−1C)TQ+Q(A−S−1C) + QBS−1BTQ+CTS−1C= 0 (20)

has a solution Q>0if and only if the Hamiltonian matrix

H:= [A−BS−1C BS−1BT

−CTS−1C−(A−BS−1C)T](21)

has a Lagrangian invariant subspace, i.e., there exist square matrices W1,W2,Esuch that

1W2=WT

2W1,[WT

1WT

2]Thas full column rank, and

H[W1

W2]=[W1

W2]E.(22)

Moreover, if (22) holds, then Q:= W2W−1

1>0solves (20).

Proof. Suppose that Q>0 solves (20). It is straightforward to prove that [WT

1WT

2]T=

[I Q ]Tsatisfies the required conditions and (22) with E=A+BS−1BTQ.

Conversely, suppose there exist W1and W2satisfying all the described conditions. With-

out loss of generality we may assume [WT

1WT

2]Thas orthogonal columns. Then it is well

known, [6, 7] that there exists a CS decomposition, i.e., there exist orthogonal matrices U,V

and diagonal matrices ∆,Γ such that

UTW1V=[∆ 0

0 0 ],UTW2V=[Γ 0

0I],

where ∆ is invertible, and ∆2+ Γ2=I. (Note that the zero block in UTW1Vcould be void.

Partition

UTAU =[A11 A12

A21 A22 ],UTBS−1BTU=[G11 G12

12 G22 ],

UTCTS−1CU =[H11 H12

12 H22 ],VTEV =[E11 E12

E21 E22 ].

Then (22) takes the form

⎡

⎢

⎣

A11 A12 G11 G12

A21 A22 GT

12 G22

−Y11 −Y12 −AT

11 −AT

−YT

12 −Y22 −AT

12 −A22

⎤

⎥

⎦⎡

⎢

⎣

∆ 0

0 0

Γ 0

⎤

⎥

⎦=⎡

⎢

⎣

∆ 0

0 0

Γ 0

⎤

⎥

⎦[E11 E12

E21 E22 ].

By comparing the (2,2) block on both sides one has G22 = 0. Due to the positive semi-

definiteness of S, then also G12 = 0. By comparing the (1,2) blocks on both sides of the above

equation and using the nonsingularity of ∆, it follows that E12 = 0. Then by comparing the

(3,2) block on both sides one has A21 = 0. Hence,

UTAU =[A11 A12

0A22 ],UTBS−1BTU=[G11 0

0 0 ],

implying that (A,B) is not controllable if W1is not invertible. Therefore, W1must be

invertible. Similarly, the observability of (A,C) implies that W2is invertible. Define Q=

W2W−1

1, which is then invertible and also symmetric because of the relation WT

1W2=

2W1. From (22), we then obtain

H[I

Q]=[I

Q](W1RW−1

1),

which implies that Qsolves (20). The positive definiteness of Qfollows from the fact that Q

is invertible, Ais asymptotically stable, and QBS−1BTQ+CTS−1C≥0.

Remark 2 The previous lemmas still hold if Ais replaced by A−BS−1C, provided that

this matrix is also asymptotically stable. Note that (A,B) is controllable if and only if

(A−BS−1C,B) is controllable and (A,C) is observable if and only if (A−BS−1C,C) is

observable.

Lemma 6 shows that the Riccati equation (20) has a solution Q>0 whenever the Hamil-

tonian matrix Hhas a Lagrangian invariant subspace. The existence of such an invariant

subspace depends only on the purely imaginary eigenvalues of H, e.g., [14]. When such

an invariant subspace exists, then there are many such invariant subspaces. Note that the

eigenvalues of Hare symmetric with respect to the imaginary axis in the complex plane. If

(22) holds, then the union of the eigenvalues of Eand −ETform the spectrum of H. One

particular choice is that the spectrum of Eis in the closed left half complex plane, another

choice is that it is in the closed right half complex plane. The two corresponding solutions

of the Riccati equation (20) are the minimal solution Q−and the maximal solution Q+and

all other solution of the Riccati equation lie (in the Loewner ordering of symmetric matrices)

between these extremal solutions.

Example 1 Consider the example B=S= 1, C=−1, A=−1−α, where α > 0. So

Aand A−BS−1C=−αare both asymptotically stable. The Riccati equation (20) is

q2−2αq + 1 = 0, which does not have a real positive semidefinite solution when α∈(0,1).

If α= 1, then it has a unique solution q= 1 associated with E= 0. If α > 1, then it has two

solutions α±√α2−1>0 with E=±√α2−1.

Suppose that A1,B1,C1from (17) have the condensed form (18) with an orthogonal

matrix U. Partition

UTQ1U=[˜

Q11 ˜

Q12

12 ˜

Q22 ].

The reduced Riccati inequality (17) then is equivalent to

UTΨ(Q1)U=[˜

Ψ11(Q1)˜

Ψ12(Q1)

Ψ12(Q1)T˜

Ψ22(Q1)]≤0,(23)

where (leaving off arguments)

Ψ11 = ( ˜

A11 −˜

B1S−1˜

C1)T˜

Q11 +˜

Q11(˜

A11 −˜

B1S−1˜

C1) + ˜

Q11 ˜

B1S−1˜

1˜

Q11 +˜

1S−1˜

C1,

Ψ12 = ( ˜

A11 −˜

B1S−1(˜

C1−˜

1˜

Q11))T˜

Q12 +˜

Q12 ˜

A22 +˜

Q11(˜

A12 −˜

B1S−1˜

C2) + ˜

1S−1˜

Ψ22 =˜

22 ˜

Q22 +˜

Q22 ˜

A22 + ( ˜

A12 −˜

B1S−1˜

C2)T˜

Q12 +˜

12(˜

A12 −˜

B1S−1˜

C2),

+˜

12 ˜

B1S−1˜

1˜

Q12 +˜

2S−1˜

C2.

For (17) to have a solution Q1>0, it is necessary that

Ψ11(Q1) = ˜

Ψ11(˜

Q11)≤0

has a positive definite solution ˜

Q11 or equivalently, the dual Riccati inequality

(˜

A11 −˜

B1S−1˜

C1)˜

Y+˜

Y(˜

A11 −˜

B1S−1˜

C1)T+˜

Y˜

1S−1˜

C1˜

Y+˜

B1S−1˜

1≤0

has a solution ˜

Y=˜

Q−1

11 >0. Using the partitioning in (18) for ˜

Y=[Y11 Y12

12 Y22 ], then one

can write this as [Φ11(˜

Y) Φ12(˜

Φ12(˜

Y)TΦ22(˜

Y)]≤0,

where

Φ11(˜

Y) = (A11 −B1S−1C1)Y11 +Y11(A11 −B1S−1C1)T+Y11CT

1S−1C1Y11 +B1S−1BT

Again, it is necessary that Φ11(˜

Y) = Φ11(Y11)≤0 has a solution Y11 >0, or equivalently

that the dual inequality

(A11 −B1S−1C1)TQ11 +Q11(A11 −B1S−1C1) + Q11B1S−1BT

1Q11 +CT

1S−1C1≤0 (24)

has a positive definite solution. Since (A11,B1) is controllable and (A11,C1) is observable,

the inequality (24) is solvable if and only if the corresponding Riccati equation is solvable,

see [4, 25]. From Lemma 6 we have as necessary condition for the solvability of (13) or (14)

that A11 −B1S−1C1is asymptotically stable and that the Hamiltonian matrix

H=[A11 −B1S−1C1B1S−1BT

−CT

1S−1C1−(A11 −B1S−1C1)T].(25)

has a Lagrangian invariant subspace.

We now show these two conditions as well as (15) are also sufficient for the existence of a

positive definite solution by explicit constructions. Clearly, we only need to consider (16) or

(17).

We first show that ˜

Ψ11 ≤0 has a solution ˜

Q11 >0, where ˜

Ψ11 is defined in (23). Suppose

the Hamiltonian matrix Hhas a Lagrangian invariant subspace. By Lemma 6, we can deter-

mine a matrix Q0

11 >0 that solves the equation (24) with the eigenvalues of A11−B1S−1(C1−

1Q0

11) in the closed left half complex plane. Then ˜

11 =: [Q0

11 0

0 0 ]is a solution of

Ψ11 = 0, but ˜

11 is only positive semidefinite if A22 is present. If A11 −B1S−1(C1−BT

1Q0

11)

has no purely imaginary eigenvalues, or equivalently Hdoes not have purely imaginary eigen-

values, then the Hamiltonian matrix

H=[˜

A11 −˜

B1S−1˜

C1˜

B1S−1˜

−˜

1S−1˜

C1−(˜

A11 −˜

B1S−1˜

C1)T]

does not have purely imaginary eigenvalues either, since the spectrum of ˜

His the union of

the spectra of H,A22 and −A22, and A22 is asymptotically stable. Consider the Riccati

equation ˜

Ψ11 + Ξ = 0,(26)

with Ξ ≥0 being chosen such that ( ˜

A11 −˜

B1S−1˜

C1,˜

1S−1˜

C1+ Ξ) is observable; such Ξ

always exists. Recall that ( ˜

A11 −˜

B1S−1˜

C1,˜

B1) is controllable. The Riccati equation (26)

corresponds to the Hamiltonian matrix

H−[0 0

Ξ 0 ].

For a sufficiently small (in norm) Ξ, by continuity, this new Hamiltonian matrix has a La-

grangian invariant subspace (when its eigenvalues are away from the imaginary axis), and

by Lemma 6, ˜

Ψ = −Ξ≤0 has a positive definite solution ˜

Q11 with all the eigenvalues of

A11 −˜

B1S−1(˜

C1−˜

1˜

Q11) in the closed left half complex plane.

If A11 −B1S−1(C1−BT

1Q0

11) has purely imaginary eigenvalues (including 0), then there

exists an invertible matrix Lsuch that

L(A11 −B1S−1(C1−BT

1Q0

11)L−1=[Σ10

0 Σ2],

where Σ1has only purely imaginary eigenvalues and Σ2is asymptotically stable. So there

exists an invertible matrix ˜

Lsuch that

L˜

11 ˜

L−1=˜

L[A11 −B1S−1(C1−BT

1Q0

11) 0

A21 −B2S−1(C1−BT

1Q0

11)A22 ]˜

L−1=⎡

⎣

Σ10 0

0 Σ20

0 Σ32 A22 ⎤

⎦=: [Σ10

0 Σ0

2],

where ˜

11 =˜

A11 −˜

B1S−1(˜

C1−˜

1˜

11)

and Σ0

2is asymptotically stable. Since similarity transformations will not affect our analysis,

without loss of generality we may assume that both Land ˜

Lare identity matrices.

Subtracting

(˜

A11 −˜

B1S−1˜

C1)T˜

11 +˜

11(˜

A11 −˜

B1S−1˜

C1) + ˜

11 ˜

B1S−1˜

1˜

11 +˜

1S−1˜

C1= 0

from (26) yields the Riccati equation

(˜

11)T˜

Y+˜

Y˜

11 +˜

Y˜

B1S−1˜

1˜

Y+ Ξ = 0,(27)

where ˜

Y=˜

Q11 −˜

11. Repartition

B1=[B0

according to ˜

11, and set

Ξ = [0 0

0 Ξ22 ],˜

Y=[0 0

0Y2].

Then the Riccati equation (27) is reduced to

(Σ0

2)TY2+Y2Σ0

2+Y2B0

2S−1(B0

2)TY2+ Ξ22 = 0.

Since ( ˜

11,˜

B1) is controllable, so is (Σ0

2,B0

2). Recall also that Σ0

2is asymptotically stable.

Analogous to the previous case one can choose (sufficiently small) Ξ22 ≥0 with (Σ0

2,Ξ22)

observable and then the reduced Riccati equation has a positive definite solution Y2with the

eigenvalues of Σ0

2+B0

2S−1(B0

2)TY2in the closed left half complex plane. Then

Q11 =˜

11 +˜

Y=[Q0

11 0

0 0 ]+[0 0

0Y2]=⎡

⎣

11 Q0

12 0

(Q0

12)TQ0

22 0

0 0 0 ⎤

⎦+⎡

⎣

0 0 0

0Y11 Y12

0YT

12 Y22 ⎤

⎦>0,

where Q0

11 has the same size of Σ1and Y22 has the same size of A22, and is a solution to

Ψ11 =−[0 0

0 Ξ22 ]≤0,

since

A11 −˜

B1S−1(˜

C1−˜

1˜

Q11) = ˜

11 +˜

B1S−1˜

1˜

Y=[Σ1B0

1S−1(B0

2)TY2

0 Σ0

2+B0

2S−1(B0

2)TY2],

has its eigenvalues is in the closed left half complex plane.

Once we have such a solution ˜

Q11, we can solve ˜

Ψ12 = 0 for ˜

Q12. This Sylvester equation

has a unique solution because ˜

A22 =A33 is asymptotically stable and the eigenvalues of

A11 −˜

B1S−1(˜

C1−˜

1˜

Q11) are in the closed left half complex plane.

Having (leaving off arguments) solved the equality ˜

Ψ12 = 0, we finally need to solve the

inequality ˜

Ψ22 ≤0. Let

Θ = [I˜

Q−1

11 ˜

Q12

0I].

From the equations ˜

Ψ11 =−Ξ and ˜

Ψ12 = 0, we obtain

Θ−T(UTΨ(Q1)U)Θ−1= Θ−T[−Ξ 0

0˜

Ψ22 ]Θ−1=[−Ξ Ξ ˜

Q−1

11 ˜

Q12

12 ˜

Q−1

11 Ξ˜

Ψ22 −˜

12 ˜

Q−1

11 Ξ˜

Q−1

11 ˜

Q12 ].

On the other hand

Θ−T(UTQ1U)Θ−1=[˜

Q11 0

0 Π ],Π = ˜

Q22 −˜

12 ˜

Q−1

11 ˜

Q12,

and

Θ(UTA1U)Θ−1=[˜

A11 ˜

A12 −˜

A11 ˜

Q−1

11 ˜

Q12 +˜

Q−1

11 ˜

Q12 ˜

A22

0˜

A22 ]

Θ(UTB1) = UTB1=[˜

0],(C1U)Θ−1=[˜

C1˜

C2−˜

C1˜

Q−1

11 ˜

Q12 ].

The (2,2) block of Θ−T(UTΨ(Q1)U)Θ−1can be rewritten as

22Π+Π˜

A22 + ( ˜

C2−˜

C1˜

Q−1

11 ˜

Q12)TS−1(˜

C2−˜

C1˜

Q−1

11 ˜

Q12),

so that

Ψ22 =˜

22Π+Π˜

A22 + ( ˜

C2−˜

C1˜

Q−1

11 ˜

Q12)TS−1(˜

C2−˜

C1˜

Q−1

11 ˜

Q12) + ˜

12 ˜

Q−1

11 F˜

Q−1

11 ˜

Q12.

Since ˜

A22 is asymptotically stable and

(˜

C2−˜

C1˜

Q−1

11 ˜

Q12)TS−1(˜

C2−˜

C1˜

Q−1

11 ˜

Q12) + ˜

12 ˜

Q−1

11 F˜

Q−1

11 ˜

Q12 ≥0,

one can always select ˜

Ξ≥0 such that the Lyapunov equation

Ψ22 =−˜

has a unique solution Π >0. Clearly, for ˜

Q22 = Π + ˜

12 ˜

Q−1

11 ˜

Q12, the matrix

Q1=U[˜

Q11 ˜

Q12

12 ˜

Q22 ]UT>0

satisfies

Ψ(Q1) = −U[ξ0

0˜

Ξ]UT≤0,

i.e., Q1>0 solves (16) and (17).

We summarize the conditions for the existence of a positive definite solution of the matrix

inequality (13) in the following theorem.

Theorem 7 Consider a general LTI system of the form (4) with Astable and S=D+DT>

0. Let Mbe invertible such that

MAM−1=[A10

0A2],MB =[B1

B2],CM−1=[C1C2],

where A2is diagonalizable and contains all the purely imaginary eigenvalues of A. Let

A1,B1,C1have the condensed form (19) with an orthogonal matrix U. Then the LMI (13)

has a positive definite solution Qif and only if the following conditions hold.

(a) There exists a positive definite matrix Q2satisfying

2Q2=C2,A2Q2=Q2A2.

(b) The block A11 −B1S−1C1is asymptotically stable.

If these conditions hold, then the LMI (13) has a positive definite solution of the form

Q=MT[Q10

0Q2]M>0, where Q1>0solves (17) and Q2is determined from condition

(a), and

(A−BS−1C)TQ+Q(A−BS−1C) + QBS−1BTQ+CTS−1C

=−MT⎡

⎣U[Ξ 0

0˜

Ξ]UT0

0 0 ⎤

⎦M≤0,

where Ξ,˜

Ξas in above construction.

Proof. The proof follows from the given construction.

Remark 3 Relation (22) is equivalent to the property that the even matrix pencil

λ⎡

⎣

0I0

−I0 0

0 0 0 ⎤

⎦−⎡

⎣

0A B

AT0CT

BTC S ⎤

⎦(28)

is regular and has index at most one, i.e., the eigenvalues at ∞are all semi-simple, (this

follows because S>0) and that it has a deflating subspace spanned by the columns of

[Q−I YT]Twith Y=S−1(C−BTQ), i.e.,

⎡

⎣

0I0

−I0 0

0 0 0 ⎤

⎦⎡

⎣

−I

Y⎤

⎦(A−BY) = ⎡

⎣

0A B

AT0CT

BTC S ⎤

⎦⎡

⎣

−I

Y⎤

⎦.(29)

A pencil λN −M is called even if N=−NTand M=MT. Since even pencils have the

Hamiltonian spectral symmetry in the finite eigenvalues, this means that there are equally

many eigenvalues in the open left and in the open right half plane.

Numerically, to compute Qit is preferable to work with the even pencil (29) over the

Hamiltonian matrix (21), since explicit inversion of Scan be avoided. Numerically stable

structure preserving methods for this task are available, see [2, 3].

Note that this procedure can be directly applied to the original data in (4) if S=D+DT

is positive definite and we attempt to solve the Riccati equation of (14). The situation is more

complicated if Sis singular, because then it is well-known that the pencil (28) has Jordan

blocks of size larger than 1 (the index of (28) is larger than 1) associated with the eigenvalue

∞, which means that the number of finite eigenvalues is less than 2n. In this case one has to

first deflate the high index and possibly singular part of the system. This can be done with

the help of the staircase form [8] which has been implemented in a production code [5]. The

situation is also more difficult if Sis indefinite, since in this case both BS−1BTand CTS−1C

cannot be guaranteed to be positive semidefinite and Lemma 6 may no longer be valid.

Remark 4 To show that the system (4) is passive it is sufficient that the LMI (13) has a

positive semidefinite solution Q. In this case the conditions can be relaxed. First of all,

the condition BT

2Q2=C2can be relaxed to Ker B2⊆Ker CT

2, rank C2B2= rank C2, and

C2B2≥0. Also, A2may have purely imaginary eigenvalues with Jordan blocks. For example

in the extreme case when C2= 0, Q2= 0 always satisfies the conditions for any A2.

Secondly, for (16) or (17), we still require that A11 −B1ˆ

S−1C1is asymptotically stable

and that the Hamiltonian matrix ˆ

Hhas a Lagrangian invariant subspace. Since this only

requires that Q1≥0, a solution can be determined in a simpler way. We may simply set

Q11 =˜

11 and ˜

Ψ11 = 0. Then with the block structure the solution of ˜

Ψ12 = 0 has a form

Q12 =[Q12

0]. To solve ˜

Ψ22 = 0 for ˜

Q22, one can show Q1=U[˜

Q11 ˜

Q12

12 ˜

Q22 ]UT≥0 and

solves the Riccati equation Ψ(Q1) = 0. Also, under certain conditions ˜

A22 =A33 does not

have to be asymptotically stable.

This shows that for non-minimal systems passivity is more general than the PH property.

3 Construction of port-Hamiltonian realizations

In the last section we have seen that the existence of a port-Hamiltonian realization for (4)

reduces to the existence of a nonsingular matrix T∈Rn×nor a positive definite matrix

Q=TTTsuch that the matrix inequalities (9) or (7) hold. In this section we develop a

constructive procedure to check these conditions. In our previous considerations the matrix

Vthat is used for basis change in the input space suffices to be invertible, but to implement

the transformations in a numerical stable algorithms, we require Vto be real orthogonal.

Consider V0= [V0,1,V0,2], where V0,1is chosen so that the columns of V0,1form an

orthonormal basis of the kernel of D+DT. To construct such a V0we can use a singular value

or rank-revealing QR decomposition, [16]. Then for the symmetric part Sof the feedthrough

matrix we have

S=1

2(VT

0DV0+VT

0DV0) = 1

2VT

0(D+DT)V0=[0 0

0S2],(30)

where S2is symmetric and invertible. Set

[B1B2]=BV0,[CT

1CT

2]:= CTV0,(31)

each partitioned compatibly with V0as in (30).

Scaling the second block row and column of the matrix inequality (7) with VT

0and V0

respectively, the matrix inequality

⎡

⎣

ATQ+QA QB1−CT

1QB2−CT

1Q−C10 0

2Q−C20−S2⎤

⎦≤0 (32)

has a solution Q=QT>0 if and only if the matrix inequality

[ATQ+QA QB2−CT

2Q−C2−S2]≤0 (33)

has a positive definite solution Qsatisfying the constraint BT

1Q−C1. We will characterize

conditions when this constraint is satisfied in the following subsections. For the solution of

the matrix inequality, then we can restrict the input and output space to the invertible part

of D+DT, and once the transformation matrices have been constructed, this transformation

can be extended to the full system, making sure that the constraint is satisfied.

3.1 Construction of the transformation in the case D =−DT

To explicitly construct the transformation to port Hamiltonian form let us first discuss the

special case that D=−DT, i.e., S= 0. Considering the matrix Kin (2), to have K≥0, we

must have P= 0, and the block V0,2in the transformation of the feedthrough term is void,

while V0=V0,1is any orthogonal matrix. Then we can express the conditions (32) and (33)

in terms of the factor Tin a factorization Q=TTT.

Corollary 8 For a state-space system of the form (4) with D+DT= 0 the following two

statements are equivalent:

1. There exists a change of basis x=Tz with an invertible matrix T∈Rn×nsuch that

the resulting realization in the new basis has PH structure as in (5).

2. There exists an invertible matrix Tsuch that

a) (TB)T=CT−1and b) (TAT−1)T+TAT−1≤0.(34)

Note that without the constraint (34a), if Ais stable, then by Lemma 3 the second condition

(34b) can always be satisfied. Adding the constraint (34a), however, makes the question

nontrivial.

We have the following characterization of the transformation matrices Tthat satisfy (34a).

Lemma 9 Consider B,CT∈Rn×m, and assume that rank B=r.

a) There exists an invertible transformation Tsatisfying condition (34a) if and only if

Ker CT=Ker B,rank CB =rand CB ≥0, or equivalently, there exists an invertible

(orthogonal) matrix Vsuch that

BV =[B10],CTV=[CT

10],C1B1=YYT>0,

where B1,CT

1∈Rn×rhave full column rank and Y∈Rr×ris invertible.

b) Let NB∈Rn×(n−r)have columns that form a basis of Ker BT. If Condition a) is

satisfied, then any Tsatisfying condition (34a) has the form T=UTZT0with

T0=[NT

Y−1C1],TZ=[Z0

0I],(35)

where U∈Rn×nis an arbitrary orthogonal matrix and Z∈R(n−r)×(n−r)is arbitrary

nonsingular.

Proof. Condition (34a) is equivalent to C=BTTTT. Necessary for the existence of an

invertible Twith this property are the conditions

Ker CT=Ker TTTB =Ker B,

0≤CB =BTTTTB,(36)

rank CB = rank BTTTTB = rank B=r.

Conversely, Ker CT=Ker Bis equivalent to the existence of an orthogonal matrix V∈

Rm×msuch that

BV =[B10],CTV=[CT

10],(37)

with B1,CT

1∈Rn×rof full column rank. The conditions rank CB =rand CB ≥0 together

are equivalent to ˜

C1˜

B1>0. So there exists an invertible matrix Y(e.g. a Cholesky factor

or a positive square root) such that C1B1=YYT.

The matrix T0as in (35) is then well defined. Furthermore, T0is invertible, since if

T0Q= 0, then C1Q= 0 and NT

BQ= 0. The latter statement implies that Q∈Ran(B),

so Q=B1zfor some zand, furthermore, C1B1z= 0. This in turn implies that z= 0 and

Q= 0; so T0is injective, and hence invertible.

The invertibility of T0implies that

BTTT

0T0=V−T[BT

0][NBCT

1Y−T][NT

Y−1C1]

=V−T[0 (Y−1C1B1)T

0 0 ][ NT

Y−1C1]

=V−T[(C1B1)(C1B1)−1C1

0]=V−T[C1

0]=C.

Hence, (34a) holds with T=T0.

Now suppose that Tis any invertible transformation satisfying (34a). Then,

BTTT(TT−1

0) = CT−1

0=BTTT

which is equivalent to

1TT

0((TT−1

0)T(TT−1

0)−I)= [0 Y]((TT−1

0)T(TT−1

0)−I)= 0.

From this, it follows that

(TT−1

0)T(TT−1

0) = [ZTZ0

0I]

for some invertible matrix Z, and that TT−1

0[Z−10

0I]=Umust be real orthogonal.

In order to explicitly construct a transformation matrix T0as in (35), it will be useful

to construct bi-orthogonal bases for the two subspaces Ker BTand Ker C. Toward this end,

let NCcontain columns that form a basis of Ker C, so that Ran NC=Ker C. Such a matrix

is easily constructed in a numerically stable way via the singular value decomposition or a

rank-revealing QR decomposition of C, see [16]. Since Band CTare assumed to satisfy

Ker B=Ker CT, we have singular value decompositions

B= [UB,1UB,2][ΣB0

0 0 ]VBT,C=UC[ΣC0

0 0 ][VT

C,1

C,2](38)

with ΣB,ΣC∈Rr×rboth invertible. We then obtain

NC=VC,2,NB=UB,2.(39)

Observe that NT

BNCis nonsingular, since if NT

BNCQ= 0 and z=NCQ, then NT

Bz= 0

implies that z∈Ran B. But then, z=B1Y=NCQimplies C1B1Y= 0, where B1and C1

are defined in (37), and so, Y= 0 and hence z= 0. Then, since NChas full column rank,

we have Q= 0. Thus, NT

BNCis injective, hence invertible.

Performing another singular value decomposition, NT

BNC=˜

U˜

∆˜

VT, with ˜

∆ positive

diagonal, and ˜

U,˜

Vreal orthogonal, we can perform a change of basis ˜

NB=NB˜

U˜

∆−1/2and

NC=NC˜

V˜

∆−1/2and obtain that the columns of ˜

NBform a basis for Ker BT, the columns

of ˜

NCform a basis for Ker Cand these two bases are bi-orthogonal, i.e., ˜

B˜

NC=I, and we

have

T0=[˜

Y−1C1],T−1

0=[˜

NCB1Y−T].(40)

Using the formula (35), we can express the conditions for a transformation to PH form

that we have obtained so far in a more concrete way.

Corollary 10 Consider system (4) with D=−DTand rank B=r. Let the columns of ˜

and ˜

NCspan the kernels of BT,Cand satisfy ˜

B˜

NC=I. Then system (4) is equivalent to

a PH system if and only if

1. Ker CT=Ker B,rank CB =r,CB ≥0, and

2. there exists an invertible matrix Zsuch that

[Z0

0I]T0AT−1

0[Z−10

0I]+([Z0

0I]T0AT−1

0[Z−10

0I])T

≤0,(41)

and T0,T−1

0are defined in (40).

Proof. The condition follows from Corollary 8 and the representation (35) by setting U=I

and T0as in (40).

3.2 Construction of the transformation in the case of general D

For the case that Dis general we will present a recursive procedure. The first step is to

perform the transformations (30), (31), and to obtain the following characterization when a

transformation to PH form exists.

Lemma 11 Consider system (4) transformed as in (30) and (31). Let the columns of ˜

NB1

and ˜

NC1form the kernels of BT

1,C1respectively, and satisfy ˜

NC1=I.

Then the system is equivalent to a PH system if and only if

1. Ker CT

1=Ker B1,rank C1B1= rank B1, and C1B1≥0, and

2. there exists an invertible matrix Zsuch that

Y+˜

YT≥0,(42)

with

Y:= [Z0

0I][−T0AT−1

0−T0B2

C2T−1

0VT

0,2DV0,2][Z−10

0I](43)

and

T0=[˜

Y−1C1],T−1

0=[˜

NC1B1Y−T],C1B1=YYT,(44)

and B1,C1are defined similar as in (37) just with B,Cbeing replaced by B1,C1.

Proof. Condition (34a) in this case has the form (TB1)T=C1T−1. If this condition holds,

then applying the result of Lemma 9 to B1and C1one has the form (35). The result is

proved by applying this formula to (34b).

Note that ˜

Yis obtained from the one given in (34b) via a congruence (and similarity)

transformation with diag(U,I), where Uis orthogonal.

We can repartition the middle factor of ˜

Yin (43) as

[−T0AT−1

0−T0B2

C2T−1

0VT

2DV2]=⎡

⎣−˜

B1A˜

NC1−˜

B1AB1Y−T−˜

B1B2

−Y−1C1A˜

NC1−Y−1C1AB1Y−T−Y−1C1B2

C2˜

NC1C2B1Y−TVT

0,2DV0,2⎤

⎦

=: [−˜

A1−˜

C1˜

D1].(45)

Thus, we have that

Y=[Z0

0I][−˜

A1−˜

C1˜

D1][Z−10

0I],

and condition (42) is analogous to condition (9), just replacing T,A,B,C,Dwith Z,˜

A1,

B1,˜

C1,˜

D1. Hence, the existence of Zcan be checked again by using Lemma 9.

This implies that the procedure of checking the existence of a transformation from (4)

to a PH system can be performed in a recursive way. One first performs the transformation

(30) and checks whether a condition as in Part 1 of Lemma 11 does not hold, in which

case a transformation does not exist, or one has that ˜

D+˜

DTin (45) is invertible. In the

latter case, if the associated matrix inequality does not have a positive definite solution, then

a transformation does not exist. Otherwise, a transformation exists and a corresponding

transformation matrix Tcan be constructed by computing Zsatisfying ˜

Y(Z) + ˜

YT(Z)≥0

and the matrix T0is formed accordingly. After this the process is repeated in a recursive

manner.

To formalize the recursive procedure, let

G0=[−A−B

C D ],

and suppose that (34a) is satisfied. Then form

G1:= ˜

0˜

T0G0˜

T−1

0˜

V0,

where

T0=[T00

0I],˜

V0=[I0

0V0],

V0is the matrix in the decomposition (30) times a permutation that interchanges the last

block columns, and T0is obtained from (44). In this way we obtain

G1=⎡

⎢

⎣

−˜

B1A˜

NC1−˜

B1AB1Y−T−˜

B1B20

−Y−1C1A˜

NC1−Y−1C1AB1Y−T−Y−1C1B2−YT

C2˜

NC1C2B1Y−TVT

0,2DV0,2VT

0,2DV0,1

0Y VT

0,1DV0,2VT

0,1DV0,1

⎤

⎥

⎦

by using the fact that

C1T−1

0= (T0B1)T=V[0Y

0 0 ],

where Vis defined in (37) for B1,C1. By (30), we have that VT

0,2DV0,1=−VT

0,1DV0,2,

0,1DV0,1=−VT

0,1DTV0,1, so that we can express G1as

G1=⎡

⎣−˜

A1−˜

B10

C1˜

D1−Γ1

0 ΓT

1Φ1⎤

⎦.

If [−˜

A1−˜

C1˜

D1]satisfies condition (34a), then in an analogous way we construct

T1= diag(T1,I,I),˜

V1= diag(I,V1,I)

such that

G2=˜

1˜

T1G1˜

T−1

1˜

V1=⎡

⎢

⎣

−˜

A2−˜

B20 0

C2˜

D2−Γ2−Γ11

0 ΓT

2Φ2−Γ12

0 ΓT

11 ΓT

12 Φ1

⎤

⎥

⎦

and we continue. Since the size of the matrices Tjis monotonically decreasing, this procedure

will terminate after a finite number of ksteps with

Gk=˜

k−1˜

Tk−1. . . ˜

0˜

T0G0˜

T−1

0˜

V0. . . ˜

T−1

k−1˜

Vk−1=⎡

⎣−˜

Ak−˜

Bk0

Ck˜

Dk−˜

Γk

0˜

ΓT

k˜

Φk

⎤

⎦,

where ˜

Dk+˜

kis positive definite, ˜

Φk=−˜

ΦT

Tj= diag(Tj,I,I),˜

Vj= diag(Iℓj,Vj,I),

and ℓjis the size of Tjfor j= 0, . . . , k −1.

If there exists an invertible matrix Tksuch that

[Tk0

0I][−˜

Ak−˜

Ck˜

Dk][T−1

0I]+([Tk0

0I][−˜

Ak−˜

Ck˜

Dk][T−1

0I])T

≥0,

see the solvability conditions in the previous section, then with ˜

Tk= diag(Tk,I,I), we have

TkGk˜

T−1

k+ ( ˜

TkGk˜

T−1

k)T≥0.(46)

Observe that for each i, j with j > i the matrices ˜

Viand ˜

icommute with ˜

Tjand ˜

T−1

and thus setting ˜

T=˜

Tk. . . ˜

T0,˜

V=˜

V0. . . ˜

Vk−1,

then ˜

TkGk˜

T−1

k=˜

VT˜

TG0˜

T−1˜

and inequality (46) implies that

TG0˜

T−1+ ( ˜

TG0˜

T−1)≥0.

Then the desired transformation matrix Tis positioned in the top diagonal block of ˜

T, and

the matrix Vis positioned in the bottom diagonal block of ˜

Remark 5 The recursive procedure described above requires at each step the computation

of three singular value decompositions, to check the ranks of the matrices Bj,Cjand to

construct the bi-orthogonal bases so that (44) holds.

While each step of this procedure can be implemented in a numerically stable way, the

consecutive rank decisions and the consequence of performing these subtle decisions are hard

to analyze, similar to the case of staircase algorithms [8, 12, 13]. In general the strategy should

be adapted to the goal of the computation, which in our case is to obtain a representation in

PH form that is robust to small perturbations.

3.3 Explicit solution of linear matrix inequalities via even pencils

We have seen that to check the existence of the transformation to PH form and to explicitly

construct the transformation matrices T,Vwe can consider the solution of the linear matrix

inequality (7). As we have discussed before the best way to do this is via the transformation

of the even pencil (28). In this subsection we combine the recursive procedure with the

construction of a staircase like form for this even pencil..

Let for a given real symmetric matrix Qdenote by

Ψ0(Q) := [−ATQ−QA CT−QB

C−BTQ D +DT]

the corresponding block matrix which is supposed to be positive semidefinite.

Let V0,1,B1,C1be defined as in (30), (31). If B1,C1satisfy Part 1. of Lemma 11 then,

since Ker CT

1=Ker B1, there exist real orthogonal matrices ˜

U1,˜

V1(which can be obtained

by performing a permuted singular value decomposition of B1) such that

1B1˜

V0,1=[0 0

ΣB0],˜

0,1C1˜

U1=[C1,1C1,2

0 0 ](47)

where ΣBis invertible and C1,2ΣBis real symmetric and positive definite. Transforming the

desired Qcorrespondingly as

Q:= ˜

1Q˜

U1=[Q11 Q12

12 Q22 ],

then, since the linear matrix inequality (7) implies ˜

QB1=CT

1, it follows in the transformed

variables that

Q22 := CT

1,2Σ−1

B>0,Q12 := CT

1,1Σ−1

B,(48)

and it remains to determine Q11 so that ˜

Q>0. To achieve this, we set

Q0:= [Q12Q−1

22 QT

12 Q12

12 Q22 ]=T−T

0[0 0

0Q22 ]T−1

0,T0=[I0

−Q−1

22 QT

12 I]

and we clearly have that Q0≥0. Then we can rewrite ˜

Qas

Q=[Q10

0 0 ]+Q0=T−T

0([Q10

0 0 ]+[0 0

0Q22 ])T−1

0=T−T

0[Q10

0Q22 ]T−1

partitioned analogously and we obtain that ˜

Q>0 if and only if Q1>0. By performing a

congruence transformation on Ψ0(Q) with

Z0:= [T00

0V0],T0=˜

U1T0,V0=V0[˜

V0,10

0I]

and using the fact that ˜

Q(˜

1B1˜

V0,1)=(˜

0,1C1˜

U1)Tfor any real symmetric ˜

Q1, it follows

that

0Ψ0(Q)Z0=⎡

⎣−(T−1

0AT0)TTT

0QT0−TT

0QT0(T−1

0AT0) 0 TT

0CT

2−TT

0QT0(T−1

0B2)

0 0 0

C2T0−(T−1

0B2)TTT

0QT00 2S2⎤

⎦.

Partitioning

T−1

0AT0=: [A11 A12

A21 A22 ],T−1

0B2=[B13

B23 ],C2T0=[C31 C32 ],

and using that

0QT0=[Q10

0Q22 ],

we obtain that

0Ψ0(Q)Z0=⎡

⎢

⎣

−AT

11Q1−Q1A11 −Q1A12 −AT

21Q22 0CT

31 −Q1B13

−AT

12Q1−Q22A21 −AT

22Q22 −Q22A22 0CT

32 −Q22B23

0 0 0 0

C31 −BT

13Q1C32 −BT

23Q22 0 2S2

⎤

⎥

⎦.

In this way, we have that (7) holds for some Q>0 if and only if

Ψ1(Q1) := ⎡

⎣−AT

11Q1−Q1A11 −Q1A12 −AT

21Q22 CT

31 −Q1B13

−AT

12Q1−Q22A21 −AT

22Q22 −Q22A22 CT

32 −Q22B23

C31 −BT

13Q1C32 −BT

23Q22 2S2⎤

⎦

=: [−AT

1Q1−Q1A1CT

1−Q1B1

C1−BT

1Q1D1+DT

1]≥0

holds for some real symmetric positive definite Q1, where

A1=A11,C1=[−Q22A21

C3,1],B1=[A12 B1,3],

and

D1+DT

1=[−AT

22Q22 −Q22A22 CT

32 −Q22B23

C32 −BT

23Q22 2S2].

This construction has reduced the solution of the linear matrix inequality (7) to the solution

of a smaller linear matrix inequality of the same form. Thus, we can again proceed in a

recursive manner with the same reduction process until either the condition in Part 1. of

Lemma 9 no longer holds (in which case no solution exists) or Dk+DT

kis positive definite

for some k.

This reduction process can be considered as the construction of a structured staircase

form for the even pencil (28). By applying a congruence transformation to the pencil (28)

with the matrix

Y0=⎡

⎣

T−T

00 0

0T00

0 0 V0⎤

⎦,

it follows that

⎡

⎢

⎣

0 0 I0 0 0 0

0 0 0 I000

−I0 0 0 0 0 0

0−I00000

0 0 0 0 0 0 0

⎤

⎥

⎦

−

⎡

⎢

⎣

0 0 A11 A12 0 0 B13

0 0 A21 A22 ΣB0B23

11 AT

21 0 0 0 0 CT

12 AT

22 0 0 CT

12 0CT

0 ΣT

B0C12 0 0 0

0 0 0 0 0 0 0

13 BT

23 C31 C32 0 0 ˜

⎤

⎥

⎦

where ˜

S1:= 2S2. By performing another congruence transformation with the matrix

Y0=

⎡

⎢

⎣

I0 0 0 0 0 0

0I0−Q22 0 0 0

0 0 I0 0 0 0

0 0 0 I0 0 0

0 0 0 0 I0 0

0 0 0 0 0 I0

0 0 0 0 0 0 I

⎤

⎥

⎦

the pencil becomes

⎡

⎢

⎣

0 0 I0 0 0 0

0 0 0 I000

−I0 0 0 0 0 0

0−I00000

0 0 0 0 0 0 0

⎤

⎥

⎦

−

⎡

⎢

⎣

0 0 A11 A12 0 0 B13

0 0 A21 A22 ΣB0B23

11 AT

21 0−AT

21Q22 0 0 CT

12 AT

22 −Q22A21 −AT

22Q22 −Q22A22 0 0 CT

32 −Q22B23

0 ΣT

B0 0 0 0 0

0 0 0 0 0 0 0

13 BT

23 C31 C32 −BT

23Q22 0 0 ˜

⎤

⎥

⎦

By further moving the last block row and column to the fifth position and then the 2nd block

row and column to the fifth position, i.e., by performing a congruence permutation with

⎡

⎢

⎣

I0 0 0 0 0 0

0 0 0 0 I0 0

0I00000

0 0 I0 0 0 0

0 0 0 0 0 I0

0 0 0 0 0 0 I

0 0 0 I000

⎤

⎥

⎦

one obtains

⎡

⎢

⎣

0I0 0 0 0

−I0 0 0 0 0

0 0 0 −Γ10 0

0 0 ΓT

10 0 0

0 0 0 0 0 0

⎤

⎥

⎦−

⎡

⎢

⎣

0A1B10 0 0

10CT

1AT

21 0 0

1C1D1+DT

1∆10 0

0A21 ∆T

10 Σ10

0 0 0 ΣT

10 0

0 0 0 0 0 0

⎤

⎥

⎦

where

Γ1=[I

0],∆1=[AT

23 ],Σ1:= ΣB

and A1,B1,C1,D1+DT

1are as defined before. In this way, we may repeat the reduction

process on the (1,1) block, which corresponds to Ψ1. In order to exploit the block structures

of the pencil we use a slightly different compression technique for D1+DT

1. Note that we

may write

D1+DT

1=[D11 D12

12 S1],

with S1symmetric positive definite. Then we have

D1+DT

1=[I0

S−1

1DT

12 I]T[D11 −D12S−1

1DT

12 0

0S1][ I0

S−1

1DT

12 I].

Let

D11 −D12S−1

1DT

12 =Z1[0 0

0˜

S2]ZT

where ˜

S2is invertible and Z1is orthogonal. Then

[Z10

−S−1

1DT

12Z1I]T

(D1+DT

1)[Z10

−S−1

1DT

12Z1I]=⎡

⎣

0 0 0

0˜

S20

00S1⎤

⎦=: [0 0

0S2].

A necessary condition for the existence of a transformation to PH form is that S2>0 or

equivalently ˜

S2>0. If this holds, then using the fact

[Z10

−S−1

1DT

12Z1I]T

Γ1=[ZT

0],

by performing a congruence transformation on the 3rd block rows and columns with

[Z10

−S−1

1DT

12Z1I]

and another congruence transformation on the fourth block row and column with Z1we

obtain the pencil

⎡

⎢

⎣

0I0 0 0 0 0

−I0 0 0 0 0 0

0 0 0 0 −Γ11 0 0

0 0 0 0 −Γ12 0 0

0 0 ΓT

11 ΓT

12 0 0 0

0 0 0 0 0 0 0

⎤

⎥

⎦

−

⎡

⎢

⎣

0A1B11 B12 0 0 0

10CT

11 CT

21 ∆11 0 0

11 C11 0 0 ∆21 0 0

12 C21 0S2∆31 0 0

0 ∆T

11 ∆T

21 ∆T

31 0 Σ10

0 0 0 0 ΣT

10 0

0 0 0 0 0 0 0

⎤

⎥

⎦

where [Γ11

Γ12 ]= Γ = ⎡

⎣

0 0 ⎤

⎦.

In order to proceed, B11,CT

11 must satisfy the same conditions as B1,CT

1. If these conditions

hold, then we can perform a second set of congruence transformation and transform the pencil

⎡

⎢

⎣

0I0 0 0 0 0 0 0

−I0 0 0 0 0 0 0 0

0 0 0 −Γ20 0 −Γ11 0 0

0 0 ΓT

20 0 0 0 0 0

0 0 0 0 0 0 −Γ12 0 0

0 0 0 0 0 0 −Γ13 0 0

0 0 ΓT

11 0 ΓT

12 ΓT

13 0 0 0

0 0 0 0 0 0 0 0 0

⎤

⎥

⎦

−

⎡

⎢

⎣

0A2B20 0 0 0 0 0

20CT

2ΘT

20 0 ∆11 0 0

2C2D2+DT

2∆20 0 ∆21 0 0

0 Θ2∆T

20 Σ20 0 0 0

0 0 0 ΣT

20 0 ∆41 0 0

0 0 0 0 0 0 ∆51 0 0

0 ∆T

11 ∆T

21 0 ∆T

41 ∆T

51 0 Σ10

0 0 0 0 0 0 ΣT

10 0

0 0 0 0 0 0 0 0 0

⎤

⎥

⎦

where ⎡

⎣

Γ11

Γ12

Γ13 ⎤

⎦=⎡

⎣

Γ12

Γ11 ⎤

⎦,Γ2=[I

0].

This reduction process continues as long as all the required conditions hold, until for some k,

Dk+DT

k>0. If this is the case, then the pencil (28) is reduced to even pencil that has the

eigenvalue ∞with equal algebraic and geometric multiplicity (it is of index one)

λ⎡

⎣

0I0

−I0 0

0 0 0 ⎤

⎦−⎡

⎣

0AkBk

l0CT

kCkDk+DT

⎤

⎦,

and if it has a deflating subspace spanned by a matrix as in (29), then the passive system (4)

can be transformed to a PH system.

Note the above process is actually a special staircase form reduction process that deflates

the singular part and higher index of the eigenvalue infinity of the even pencil (28), [5, 8].

Remark 6 Note that to check the passivity of (4) it is only necessary to have a positive

semidefinite solution to (7). Thus, if one only wants to check passivity, then Part 1. in

Lemma 9,1 can be relaxed to Ker B1⊆Ker CT

1. In this case the transformation to the

form (47) can still be made, but Q22 in (48) is only positive semidefinite, and Ker Q22 ⊆

Ker Q12. Then Q0can still be defined but instead of Q−1

11 one needs to use the Moore-

Penrose pseudoinverse, see [16], of Q11. However, in this case Q0and the resulting solution

Qcannot be positive definite. Thus, in this situation, (4) may be a passive system that

cannot be transformed to a PH system. A simple example is the scalar system

z=−z+ 2u,y= 0z+ 0u.

It is passive (but not strictly passive) with S=zT[0]z, since yTu= 0. But it cannot be

transformed to a PH system, since for any Q=αzand y=βy,u=β−1uwith α, β = 0, one

has ˙

Q=−Q+ 2αβu,y= 0Q+ 0u.

So we would obtain J= 0,R= 1,B=−P=αβ and S= 0. However,

K=[1−αβ

−αβ 0]

is indefinite for any αβ = 0.

What goes wrong in this case in the transformation to PH form is that the resulting

transformation matrix Tthat is needed to satisfy the condition on the input and output

matrices becomes singular.

To illustrate the analysis procedures, consider the following example.

Example 2 In the finite element analysis of disc brake squeal [17] large scale (approx. 1

million degrees of freedom) second order differential equations arise which depend on param-

eters, e.g., the disc speed ω. If no further constraints are incorporated, then in the stationary

case they take the form

M¨

q+ (C1+ωr

ωCR+ω

ωr

CG)˙

q+ (K1+KR+ ( ω

ωr

)2KG)q= 0,

where M=MT>0 is the mass matrix, C1=CT

1≥0 models material damping, CG=−CT

models gyroscopic effects, CR=CT

R≥0 models friction induced damping, K1=KT

1>0

is the stiffness matrix, KR=K2+Nwith K2=KT

2and N=−NT, is a nonsymmetric

matrix modeling circulatory effects, KG=KT

G≥0 is the geometric stiffness matrix, and

ωis the rotational speed of the disc with reference velocity ωr. In industrial brake models,

the matrices D:= C1+ωr

ωCR, and Nare sparse and have very low rank (approx 2000)

corresponding to finite element nodes associated with the brake pad. Setting G:= ω

ωrCG),

K=K1+K2+ ( ω

ωr)2KG, we may assume that K>0

Thus, introducing p=M˙

q, we can write the system in first order form

[˙

q]=([ −G−(I+1

2NK−1)

(I+1

2K−1N) 0 ]−[D1

2NK−1

2K−1N0])[M−10

0K][p

:= (J−R)Q[p

q].

Since

−R:= [D1

2NK−1

2NK−10]

is indefinite, as soon as N= 0 it is clear that this system is not automatically stable and it

is definitely unstable if x∗Rx <0 for any eigenvector xof

J:= [−G−(I+1

2NK−1)

(I+1

2NK−1) 0 ].

4 Conclusions

In this paper we have extended the characterization of [38] when a system is equivalent to

a port Hamiltonian system to the case of general non-minimal systems. We have presented

an explicit construction procedure for the construction of the transformation matrices, which

can be implemented as a numerical method. This implementation is currently under con-

struction.

References

[1] C. Beattie and S. Gugercin, Structure-preserving model reduction for nonlinear port-

Hamiltonian systems, 50th IEEE Conference on Decision and Control and European

Control Conference (CDC-ECC), 2011, IEEE, 2011, pp. 6564–6569.

[2] P. Benner, R. Byers, V. Mehrmann, and H. Xu, Numerical methods for linear-quadratic

and H∞control problems, Dynamical Systems, Control, Coding, Computer Vision: New

Trends, Interfaces, and Interplay (G. Picci and D.S. Gilliam, eds.), Progress in Systems

and Control Theory, vol. 25, Birkh¨auser, Basel, 1999, pp. 203–222.

[3] P. Benner, P. Losse, V. Mehrmann, and M. Voigt, Numerical linear algebra methods

for linear differential-algebraic equations: A survey, DAE Forum, Springer, 2015, p. To

appear.

[4] S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan, Linear matrix inequalities in

systems and control theory, SIAM, Philadelphia, 1994.

[5] T. Br¨ull and V. Mehrmann, STCSSP: A FORTRAN 77 routine to compute a structured

staircase form for a (skew-)symmetric/(skew-)symmetric matrix pencil, Preprint 31-2007,

Institut f¨ur Mathematik, TU Berlin, 2007.

[6] R. Byers, Hamiltonian and symplectic algorithms for the algebraic Riccati equation, Ph.D.

thesis, Cornell University, Dept. Comp. Sci., Ithaca, NY, 1983.

[7] , A Hamiltonian QR-algorithm, SIAM J. Sci. Statist. Comput. 7(1986), 212–229.

[8] R. Byers, V. Mehrmann, and H. Xu, A structured staircase algorithm for skew-

symmetric/symmetric pencils, Electron. Trans. Numer. Anal. 26 (2007), 1–13.

[9] C.I. Byrnes and A. Isidori, Limit sets, zero dynamics, and internal models in the problem

of nonlinear output regulation, IEEE Transactions on Automatic Control 48 (2003),

no. 10, 1712–1723.

[10] C.I. Byrnes, A. Isidori, and J.C. Willems, Passivity, feedback equivalence, and the global

stabilization of minimum phase nonlinear systems, IEEE Trans. Autom. Control 36

(1991), 1228–1240.

[11] M. Dalsmo and A.J. van der Schaft, On representations and integrability of mathematical

structures in energy-conserving physical systems, SIAM J. Control Optim. 37 (1999), 54–

91.

[12] K. Demmel and B. K˚agstr¨om, The generalized Schur decomposition of an arbitrary pen-

cil A−zB: Robust software with error bounds and applications. Part I: Theory and

algorithms, ACM Trans. Math. Softw. 19 (1993), no. 2, 160–174.

[13] , The generalized Schur decomposition of an arbitrary pencil A−zB: Robust

software with error bounds and applications. Part II: Software and applications, ACM

Trans. Math. Softw. 19 (1993), no. 2, 175–201.

[14] G. Freiling, V. Mehrmann, and H. Xu, Existence, uniqueness and parametrization of

Lagrangian invariant subspaces, SIAM J. Matrix Anal. Appl. 23 (2002), 1045–1069.

[15] G. Golo, A.J. van der Schaft, P.C. Breedveld, and B.M. Maschke, Hamiltonian formula-

tion of bond graphs, Nonlinear and Hybrid Systems in Automotive Control (A. Rantzer

R. Johansson, ed.), Springer, Heidelberg, 2003, pp. 351–372.

[16] G. H. Golub and C. F. Van Loan, Matrix computations, 3rd ed., Johns Hopkins Univ.

Press, Baltimore, 1996.

[17] N. Gr¨abner, V. Mehrmann, S. Quraishi, and C. Schr¨odera nd U. von Wagner, Numeri-

cal methods for parametric model reduction in the simulation of disc brake squeal, Tech.

Report Preprint 16/2015, Institut f¨ur Mathematik, TU Berlin, 2015, Submitted for pub-

lication.

[18] S. Gugercin, R.V. Polyuga, C. Beattie, and A.J. van der Schaft, Structure-preserving

tangential interpolation for model reduction of port-Hamiltonian systems, Automatica

48 (2012), 1963–1974.

[19] D. Hinrichsen and A. J. Pritchard, Mathematical systems theory I. modelling, state space

analysis, stability and robustness, Springer-Verlag, New York, NY, 2005.

[20] R. A. Horn and C. R. Johnson, Matrix analysis, Cambridge University Press, Cambridge,

1985.

[21] B. Jacob and Zwart H, Linear port-Hamiltonian systems on infinite-dimensional spaces,

Operator Theory: Advances and Applications, 223, Birkh¨auser/Springer Basel AG, Basel

CH, 2012.

[22] T. Kailath, Linear systems, vol. 1, Prentice-Hall Englewood Cliffs, NJ, 1980.

[23] C. Kleijn, 20-sim 4c 2.1 reference manual, Controlab Products B.V., 2013.

[24] N. Kottensette and P.J. Antsaklis, Relationships between positive real, passive, dissipative

and positive systems, American Control Conference (Baltimore, MD, USA), 2010.

[25] P. Lancaster and L. Rodman, The algebraic Riccati equation, Oxford University Press,

Oxford, 1995.

[26] P. Lancaster and M. Tismenetsky, The theory of matrices, 2nd ed., Academic Press,

Orlando, 1985.

[27] B.M. Maschke, A.J. van der Schaft, and P.C. Breedveld, An intrinsic Hamiltonian for-

mulation of network dynamics: non-standard poisson structures and gyrators, J. Franklin

Inst. 329 (1992), 923–966.

[28] V. Mehrmann, The autonomous linear quadratic control problem, theory and numerical

solution, Lecture Notes in Control and Inform. Sci., vol. 163, Springer-Verlag, Heidelberg,

July 1991.

[29] R. Ortega, A.J. van der Schaft, Y. Mareels, and B.M. Maschke, Putting energy back in

control, Control Syst. Mag. 21 (2001), 18–33.

[30] R. Ortega, A.J. van der Schaft, B.M. Maschke, and G. Escobar, Interconnection and

damping assignment passivity-based control of port-controlled Hamiltonian systems, Au-

tomatica 38 (2002), 585–596.

[31] R. V. Polyuga and A.J. van der Schaft, Structure preserving model reduction of port-

Hamiltonian systems by moment matching at infinity, Automatica 46 (2010), 665–672.

[32] A.J. van der Schaft, Port-Hamiltonian systems: an introductory survey, Proc. of the

International Congress of Mathematicians, vol. III, Invited Lectures (Madrid, Spain)

(J.L. Verona M. Sanz-Sole and J. Verdura, eds.), pp. 1339–1365.

[33] , Port-Hamiltonian systems: network modeling and control of nonlinear physical

systems, Advanced Dynamics and Control of Structures and Machines, CISM Courses

and Lectures, Vol. 444, Springer Verlag, New York, N.Y., 2004.

[34] A.J. van der Schaft and B.M. Maschke, The Hamiltonian formulation of energy conserv-

ing physical systems with external ports, Arch. Elektron. ¨

Ubertragungstech. 45 (1995),

362–371.

[35] , Hamiltonian formulation of distributed-parameter systems with boundary energy

flow, J. Geom. Phys. 42 (2002), 166–194.

[36] , Port-Hamiltonian systems on graphs, SIAM J. Control Optim. (2013).

[37] P. Van Dooren, The computation of Kronecker’s canonical form of a singular pencil,

Linear Algebra Appl. 27 (1979), 103–121.

[38] J. C. Willems, Dissipative dynamical systems – Part I: General theory, Arch. Ration.

Mech. Anal. 45 (1972), 321–351.

[39] , Dissipative dynamical systems – Part II: Linear systems with quadratic supply

rates, Arch. Ration. Mech. Anal. 45 (1972), 352–393.