Computation of the analytic center of the solution set of the linear matrix inequality arising in continuous- and discrete-time passivity analysis

In this paper formulas are derived for the analytic center of the solution set of linear matrix inequalities (LMIs) defining passive transfer functions. The algebraic Riccati equations that are usually associated with such systems are related to boundary points of the convex set defined by the solution set of the LMI. It is shown that the analytic center is described by closely related matrix equations, and their properties are analyzed for continuous- and discrete-time systems. Numerical methods are derived to solve these equations via steepest ascent and Newton-like methods. It is also shown that the analytic center has nice robustness properties when it is used to represent passive systems. The results are illustrated by numerical examples.


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 transfer functions which play a fundamental role in systems and control theory: they represent e. g., spectral density functions of stochastic processes, show up in spectral factorizations, are the Hermitian part of a positive real transfer function, characterize port-Hamiltonian systems, and are also related to algebraic Riccati equations.
Positive transfer functions form a convex set, and this property has lead to the extensive use of convex optimization techniques in this area (especially for so-called linear matrix inequalities [5]). In order to optimize a certain scalar function f (X) over a convex set, one often defines a barrier function b(X) that becomes infinite near the boundary of the set, and then finds the minimum of c · f (X) + b(X), c 0, as c → +∞. These minima (which are functions of the parameter c) are called the points of the central path. The starting point of this path (c = 0) is called the analytic center of the set.
In this paper we present an explicit set of equations that define the analytic center of the solution set of the linear matrix inequality defining a passive transfer function. We also show how these equations relate to the algebraic Riccati equations that typically arise in the spectral factorization of transfer functions. We discuss transfer functions both on the imaginary axis (i. e. the continuous-time case), as well as on the unit circle (i. e. the discrete-time case). In the continuous-time setting the transfer function arises from the Laplace transform of the systemẋ = Ax + Bu, x(0) = 0, where u : R → C m , x : R → C n , and y : R → C m are vector-valued functions denoting, respectively, the input, state, and output of the system. Denoting real and complex n-vectors (n × m matrices) by R n , C n (R n×m , C n×m ), respectively, the coefficient matrices satisfy A ∈ C n×n , B ∈ C n×m , C ∈ C m×n , and D ∈ C m×m . In the discrete-time setting the transfer function arises from the z-transform applied to the system x k+1 = Ax k + Bu k , x 0 = 0, y k = Cx k + Du k , respectively. We restrict ourselves to systems which are minimal, i. e. the pair (A, B) is controllable (for all λ ∈ C, rank [ λI − A B ] = n), and the pair (A, C) is reconstructable (i. e. (A H , C H ) is controllable). Here, the conjugate transpose (transpose) of a vector or matrix V is denoted by V H (V T ) and the identity matrix is denoted by I n or I if the dimension is clear. We furthermore require that input and output port dimensions are equal to m and assume that rank B = rank C = m.
Passive systems and their relationships with positive-real transfer functions are well studied, starting with the works [13,17,20,21,22,23] and the topic has recently received a revival in the work on port-Hamiltonian (pH) systems, [18,19]. For a summary of the relationships see [2,20], where also the characterization of passivity via the solution set of an associated linear matrix inequality (LMI) is highlighted.
The paper is organized as follows. After some preliminaries in Section 2, in Section 3 we study the analytic centers of the solution sets of LMIs associated with the continuous-and discrete-time case. In Section 4 we discuss numerical methods to compute the analytic centers using steepest ascent as well as Newton-like methods and show that the analytic centers can be computed efficiently. In Section 5 lower bounds for the distance to non-passivity (the passivity radius) are derived using smallest eigenvalues of the Hermitian matrices associated with the linear matrix inequalities evaluated at the analytic center. The results are illustrated with some simple examples where the analytic center can be calculated analytically. In Appendix A we derive formulas for the computation of the gradients and the Hessian of the functions that we optimize and in Appendix B we clarify some of the differences that arise between the continuous-and the discrete-time case.

Preliminaries
Throughout this article we will use the following notation. We denote the set of Hermitian matrices in C n×n by H n . Positive definiteness (semidefiniteness) of A ∈ H n is denoted by A 0 (A 0). The real and imaginary parts of a complex matrix Z are written as (Z) and (Z), respectively, and ı is the imaginary unit. We consider functions over H n , which is a vector space if considered as a real subspace of R n×n + ıR n×n . We will identify C m×n with R m×n +ıR m×n , but we note that this has implications when one is carrying out differentiations, see Appendix A. The Frobenius scalar product for matrices X, Y ∈ R n×n + ıR n×n is given by where we have partitioned X, Y as X = X r + ıX i , Y = Y r + ıY i with real and imaginary parts in R n×n . As we are mainly concerned with this scalar product, we will drop the subscript R. We will make frequent use of the following properties of this inner product given by The concepts of positive-realness and passivity are well studied. In the following subsections we briefly recall some important properties following [10,20], where we repeat a few observations from [2]. See also [20] for a more detailed survey.

Positive-realness and passivity, continuous-time
Consider a continuous-time system M as in (1) and the transfer function T c as in (2). The transfer function T c (s) is called positive real if the matrix-valued rational function is positive semidefinite for s on the imaginary axis, i. e. Φ c (ıω) 0 for all ω ∈ R and it is called strictly positive real if Φ c (ıω) 0 for all ω ∈ R.
We associate with Φ c a system pencil where R := D + D H . Here (3) has a Schur complement which is the transfer function Φ c (s) and the generalized eigenvalues of S c (s) are the zeros of Φ c (s).
For X ∈ H n we introduce the matrix function If T c (s) is positive real, then the linear matrix inequality (LMI) has a solution X ∈ H n and we have the sets An important subset of X c are those solutions to (4) for which the rank r of W c (X) is minimal (i. e. for which r = rank Φ c (s)). If R is invertible, then the minimum rank solutions in X c are those for which rank W c (X) = rank(R) = m, which in turn is the case if and only if the Schur complement of R in W c (X) is zero. This Schur complement is associated with the continuous-time algebraic Riccati equation (ARE) Solutions X to (6) produce a spectral factorization of Φ c (s), and each solution corresponds to a Lagrangian invariant subspace spanned by the columns of U c := I n −X T T that remains invariant under the action of the Hamiltonian matrix i. e. U c satisfies H c U c = U c A Fc for a closed loop matrix A Fc = A − BF c with F c := R −1 (C − B H X) (see e.g., [8]). Each solution X of (6) can also be associated with an extended Lagrangian invariant subspace for the pencil S c (s) (see [4]), spanned by the columns of U c := −X T I n −F T c T . In particular, U c satisfies The sets X c , X Ï c are related to the concepts of passivity and strict passivity see [20]. If for the system M := {A, B, C, D} of (3) the LMI (4) has a solution X ∈ X c then M is (Lyapunov) stable (i.e. all eigenvalues are in the closed left half plane with any eigenvalues occurring on the imaginary axis being semisimple), and passive, and if there exists a solution X ∈ X Ï c then M is asymptocially stable, (i.e. all eigenvalues are the open left half plane) and strictly passive. Furthermore, if M is passive, then there exist maximal and minimal solutions X − X + of (4) in X c such that all solutions X of W c (X) 0 satisfy 0 ≺ X − X X + , which implies that X c is bounded. For more details on the different concepts discussed in this section, see [2].

Positive-realness and passivity, discrete-time
For each of the results of the previous subsection there are discrete-time versions which we briefly recall in this section, see [12,17]. Note, that these results can be obtained by applying a bilinear transform (see Appendix B) to the continuous-time counterparts.
The transfer function T d (s) in (2) is called positive real if the matrix-valued rational function We consider an associated the matrix function where again R = D + D H , the sets

and the system pencil
If the system is positive real then, see [20], there exists X ∈ H n such that W d (X) 0. If that is the case, a transfer function T d (z) := C(zI n − A) −1 B + D is called passive and strictly passive if even W d (X) 0. We again have an associated discrete-time Riccati equation defined as from which one directly obtains a spectral factorization of Φ d (z). The solutions of the discretetime Riccati equation can be obtained by computing a Lagrangian invariant subspace spanned by the columns of U d := I n −X T T of the symplectic matrix Each solution X of (9) can also be associated with an extended Lagrangian invariant subspace for the pencil S d (z) (see [4]), spanned by the columns of U d := −X T I n −F T d T .
In particular, U d satisfies Again, if the system is passive, then there exist maximal and minimal solutions X − X + in X Ï d , such that all solutions X of W d (X) 0 satisfy 0 ≺ X − X X + , which implies that X d is bounded.

The analytic center
If the sets X Ï c , X Ï d in (5), respectively (8), are non-empty, then we can define their respective analytic center. Following the discussion in [10], we first consider the continuous-time case, the discrete-time case is derived in an analogous way. We choose a scalar barrier function b(X) := − ln det W c (X), which is bounded from below but becomes infinitely large when W c (X) becomes singular. We define the analytic center of the domain X Ï c as the minimizer of this barrier function.

The continuous-time case
The solutions X + and X − of the Riccati equation Ricc c (X) = 0 in (6), are both on the boundary of X c , and hence are not in X Ï c . Since we assume that X Ï c is non-empty, the analytic center is well defined, see, e. g., , Section 4.2 in [16].
To characterize the analytic center, we first need to find a variation of the gradient b X of the barrier function b at point X along direction ∆ X , which is equal to where b X = W c (X) −1 and ∆W c (X)[∆ X ] is the incremental step in the direction ∆ X , for details see Appendix A. It appears that X is an extremal point of the barrier function if and only if The equation for the extremal point then becomes Defining For a point X ∈ X Ï c it is obvious that we also have P c = Ricc c (X) 0, and hence (10) is equivalent to and this is equivalent to where we have set We emphasize that P c is nothing but the Riccati operator Ricc c (X) defined in (6), and that A Fc is the corresponding closed loop matrix. For the classical Riccati solutions we have P c = Ricc c (X) = 0 and the corresponding closed loop matrix is well-known to have its eigenvalues equal to a subset of the eigenvalues of the corresponding Hamiltonian matrix (7).
Since P c = Ricc c (X) 0, it follows that P c has a Hermitian square root T c satisfying P c = T 2 c . Transforming (11) with the invertible matrix T c , we obtain is skew-Hermitian and has all its eigenvalues on the imaginary axis, and so does A Fc . Therefore, the closed loop matrix A Fc of the analytic center has a spectrum that is also central.
It is important to also note that which implies that we are also finding a stationary point of det Ricc c (X), since det R is constant and non-zero.
Since the matrix P c is positive definite and invertible, we can rewrite the equations defining the analytic center as where X = X H and P c = P H c 0. We can compute the analytic center by solving these three equations which actually form a cubic equation in X.
Note that even though the eigenvalues of the closed loop matrix F c associated with the analytic center are all purely imaginary, the eigenvalues of the original system and the poles of the transfer function stay invariant under the state space transformation T c .

The discrete-time case
For discrete-time systems, the increment of W d (X) equals and the equation for the extremal point becomes This is equivalent to which is a non-homogenous discrete-time Lyapunov equation. Since (A, B) is controllable (by assumption), so is (A Fc , B) and it follows then from (13) that the eigenvalues of A F d are now strictly inside the unit circle. This is clearly different from the continuous-time case, where the spectrum of A Fc was on the boundary of the stability region. The equations defining the discrete-time analytic center then become

Numerical computation of the analytic center
In this section we present methods for the numerical computation of the analytic center. Suppose that we are at a point X 0 ∈ X Ï c (X Ï d ) and want to perform the next step using an increment ∆ X . We discuss a steepest ascent and a Newton-like method to obtain that increment.

A steepest ascent method
In order to formulate an optimization scheme to compute the analytic center, we can use the gradient of the barrier function b(X) with respect to X to obtain a steepest ascent method.
In the continuous-time case, we then need to take a step ∆ X for which The maximum is obtained by choosing ∆ X proportional to the gradient The corresponding optimal stepsize α for the increment ∆ X can be obtained from the determinant of the incremented LMI W c (X 0 + α∆ X ) 0.
In the discrete-time case, we obtain the increment from The maximum is obtained by choosing ∆ X proportional to and the stepsize α for the increment ∆ X can again be obtained from the determinant of the incremented LMI W d (X 0 + α∆ X ) 0.
Remark 4.1. The detailed explanation how to compute the stepsize α will be done later as a special case of the derivation of the Newton step, see subsection 4.2. The idea is to find the second order Taylor expansion of the function f (X 0 + α∆ X ) = − ln det G(X 0 + α∆ X ) and then to maximize this quadratic function in the scalar α.

A Newton-like method
For the computation of a Newton-like increment ∆ X we also need the Hessian of the barrier function b. In order to simplify the derivation we first equivalently reformulate the barrier function into a more suitable form.

The continuous-time case
In the continuous-time case, we have that Carrying out an additional congruence transformation with It is clear that the determinant of the congruence transformation introduces a factor det(P 0 ). Finally, the determinant of the transformed matrix is, up to a constant det(R 0 ), equal to This is the multiplying factor of the current value of det W c (X 0 ) and we can make it larger than 1 ifÂ Fc is not skew-Hermitian yet. Introduce f (X) := − ln det(G(X)), In the set of Hermitian matrices (over the reals), the gradient of f (X) then is given by and the Hessian is given by A second order approximation of f (at X = 0) is given by and we want the gradient of f to be 0. For the Newton step we want to determine ∆ = ∆ H such that

∂T
(2) f ∂∆ (∆) = 0, i. e. we require that Using the properties of the scalar product, we obtain that this is equivalent to If we fix a direction ∆ and look for α such that f (α∆) is maximal, then the Newton step can be computed in an analogous way. With g(α) = f (α∆), we then have and thus the Newton correction in α is given by

The discrete-time case
For the discrete-time case, we have that transforming with Z from the left and Z r from the right, and substitutingB = P The determinant of the transformed matrix is equal to We introduce f (X) := − ln det(G(X)), and compute the gradient and the Hessian of f (X). The computation of the gradient is not as straight-forward as in the continuous-time case, since we consider non-Hermitian matrices. It is given by see Appendix A for more details. Revisiting the steps for the derivation of G(X), we notice that det(G(X)) is still real and the solution ∆ is still unique and Hermitian. Thus, the Hessian is given by , and a second order approximation of f (at X = 0) is given by We want the gradient of f to be 0, so for the Newton step we determine ∆ = ∆ H such that

∂T
(2) f ∂∆ (∆) = 0, or equivalently Using the properties of the scalar product, we obtain that If we fix a direction ∆ and look for α such that f (α∆) is maximal, then the Newton correction in α is given by Remark 4.2. To carry out the Newton step, we have to solve equation (15) in the continuoustime case or (16) in the discrete-time case. This can be done via Kronecker products (for the cost of increasing the system dimension to n 2 ), i. e. via in the continuous-time case, or in the discrete-time case.

Convergence
In this subsection, we show that the functions that we consider here actually have a globally converging Newton method. For this we have to analyze some more properties of our functions and refer to [6,16] for more details. Recall that a smooth function f : R n → R is selfconcordant if it is a closed and convex function with open domain and 3 2 in the case n = 1, and if n > 1, then f is self-concordant if it is self-concordant along every direction in its domain. In particular, if n = 1 then f (x) = − ln(x) is self-concordant and in general, if f is self-concordant and in addition A ∈ C n×m , b ∈ R n , then f (Ax + b) is also self-concordant. These results can be easily extended to the real space of complex matrices showing that the function b(X) = − ln det(W (X)) is self-concordant. Let the discrete-time case respectively. In both cases λ(X) can be easily computed during the Newton step and gives an estimate of the residual of the current approximation of the solution.
Furthermore, for every X ∈ X Ï c (X Ï d ) the quadratic form of the Hessian in the original coordinates can be expressed as Using the Courant-Fischer theorem twice, see e. g. [3], this implies that Note that ∆W [∆ X ] F = 0 for controllable (A, B) and ∆ X = 0. Minimizing the left-hand side over all ∆ X with ∆ X 2 F = 1 yields uniform positivity of the Hessian, since the spectrum of W (X) is bounded.
Hence, it follows, see e. g. [6], that the Newton method is quadratically convergent, whenever λ(X) < .25 in some intermediate step. Once this level is reached, the methods stays in the quadratically converging regime. If the condition does not hold, then one has to take a smaller stepsize (1 + λ(X)) −1 ∆ X in order to obtain convergence.

Initialization
Note that for the reformulations of the Newton step we have to assume that the starting value X 0 is in the interior of the domain. In this section, we show how to compute an initial point X 0 ∈ X We start the optimization from a model M that is minimal and strictly passive. It then follows that the solution set of W c (X 0 , M) 0 has an interior point X 0 0 such that To construct such an X 0 , let α := λ min W c (X 0 ) > 0 and β := max( X 0 2 , 1) > 0. Then, for 0 < 2ξ ≤ α/β, we have the inequality In order to compute a solution X 0 for this LMI, we rewrite it, up to a scaling factor, as for the modified model M ξ := {A + ξI n , B, C, D − ξI m }. The solution set of this shifted LMI can be obtained from the extremal solutions X − (ξ) and X + (ξ) of the Riccati equations for M ξ . It therefore follows that The reasoning for the discrete-time case is very similar. Starting from a strictly passive and minimal model M, we have the inequality In order to compute a solution X 0 for this LMI, we rewrite it as an LMI The solution set X − (ξ) X 0 X + (ξ) of this scaled LMI is again strictly included in the original solution set.
The procedure to find an inner point is thus to choose one of the Riccati solutions X − (ξ) or X + (ξ) of shifted or scaled problems, respectively, or some kind of average of both, since they are then guaranteed to be an interior point of the original problem.
Another possibility to compute an initial point is to take the geometric mean of the minimal and maximal solution of the Riccati equations (6), respectively (9), denoted by X − and X + , which is defined by X 0 = X − (X −1 − X + ) 1 2 , see [15]. However, e. g., if X − and X + are multiples of the identity matrix, then the geometric mean is a convex combination of X − and X + and will not necessarily be in the interior.

Numerical results
We have implemented the steepest ascent method of Subsection 4.1 and the Newton method introduced in Subsection 4.2. The software package is written in python 3.6. The code and all the examples can be downloaded under [1].
We have performed several experiments to test convergence for the different methods developed in this paper. (b) Convergence of the relative error between the current value of the objective function det(W c (X k )) and the value det(W c (X c )) at the analytic center Example 4.1. As a prototypical example consider a randomly generated continuous-time example with n = 30 and m = 10, i. e. the overall dimension of the matrix W c (X) is 40 × 40 and we have a total of 465 unknowns. As one would expect, the steepest ascent method shows linear convergence behavior, whereas the Newton method has quadratic convergence as soon as one is close enough to the analytic center. Figure 1 shows the convergence behavior using the Newton method. Note, that the barrier function det(W (X)) increases monotonously, whereas the distance of the argument X to the analytic center X c slightly increases in the linearly convergent phase. The number of steps required in the steepest ascent approach, however, is much higher than in the Newton approach.
Also, the initial point computed by the geometric mean approach turns out to be much better in all the practical examples, even though one cannot guarantee positivity in some extreme cases.
Note that one has to be extremely careful with the implementation of the algorithm. Without explicitly forcing the intermediate solutions X k to be Hermitian in finite precision arithmetic, the intermediate Riccati residuals P k may diverge from the Hermitian subspace.

Computation of bounds for the passivity radius
Once we have found a solution X ∈ X Ï c , respectively X ∈ X Ï d , we can use this solution to find an estimate of the passivity radius of our system, i. e. the smallest perturbation ∆ M to the system coefficients M = {A, B, C, D} that puts the system on the boundary of the set of passive systems, so that an arbitrary small further perturbation makes the system nonpassive. In this section we derive lower bounds for the passivity radius in terms of the smallest eigenvalue of a scaled version of the matrices W c (X, M) or W d (X, M), respectively. Since the analytic center is central to the solution set of the LMI, we choose it for the realization of the transfer function, since then we expect to maximize a very good lower bound for the passivity radius.

The continuous-time case
As soon as we fix X ∈ X Ï c , the matrix is linear as a function of the coefficients A, B, C, D. When perturbing the coefficients, we thus preserve strict passivity, as long as We thus suppose that W c (X, M) 0 and look for the smallest perturbation ∆ M to our model M that makes det W c (X, M + ∆ M ) = 0. To measure the model perturbation, we propose to use the norm of the perturbation of the system pencil We have the following lower bound in terms of the smallest eigenvalue λ min of a scaled version of W c (X, M).
Lemma 5.1. The X-passivity radius, defined for a given X ∈ X Ï c as for Proof. We first note that If we introduce the n × (n + m) matrix Z c := −X 0 , then (18) is equivalent to If we replace the matrix Z c I m+n by the matrix U c = Z c I m+n Y c with orthonormal columns, which we can e. g. obtain from a QR decomposition [11], then we obtain Therefore, the smallest perturbation of the matrix Y c W c (X, M)Y c to make Y c W c (X, M + ∆ M )Y c singular must have a 2-norm which is at least as large as λ min (Y c W c (X)Y c ), and since the perturbation is a contraction of the proposed one, the lower bound in (17) follows.

The discrete-time case
In the discrete-time case, for a fixed X the LMI takes the form and its perturbed version is where again R := D + D H and ∆ R := ∆ D + ∆ H D . Note that, in contrast to the continuous-time case, for given X ∈ X Ï d , W d (X, M + ∆ M ) is not linear in the perturbations. Nevertheless, we have an analogous bound as in Lemma 5.1 also in the discrete-time case.
Lemma 5.2. The X-passivity radius, defined for a given X ∈ X Ï d as where Proof. We first observe that since again W d (X, M+∆ M ) is just the Schur complement with respect to the leading 2n×2n matrix. Note that this matrix (20) is linear in the perturbation parameters, since X is fixed.
Using the definition of the matrix Z d , then from (20), it follows that we can consider If we replace the matrix Z d I m+n by the matrix with orthonormal columns Again, since the perturbation is a contraction of the proposed one, the (approximate) lower bound in (19) follows.

Examples with analytic solution
In this subsection, to illustrate the results, we present simple examples of scalar transfer functions (m = 1) of first degree (n = 1). Consider first an asymptotically stable continuous-time system and transfer function T (s) = d + cb s−a i. e. with a < 0. Then and its determinant is det(W c (x)) = −4adx − (c − bx) 2 , which is maximal at the central point We then get For the transfer function to be strictly passive, it must be asymptotically stable and positive on the imaginary axis and hence also at 0 and ∞. Thus, we have the conditions The function Φ c (ıω) = 2d − 2acb a 2 +ω 2 is a unimodal function, which reaches its minimum either at 0 (namely Φ c (0) = p b 2 a 2 ) or at ∞ (namely Φ c (∞) = 2d) and hence the conditions in (21) are sufficient to check passivity. Thus, for the model M, strict passivity gets lost when either one of the following happens Therefore, it follows that At the analytic center x a we have and the smallest perturbation of the parameters that makes this determinant go to 0, yields exactly the same conditions as (21). This illustrates that the X-passivity radius at the analytic center yields a very good condition for strict passivity of the model. In the discrete-time case the transfer function is T (z) = d+ cb z−a and for it be asymptotically stable we need a 2 < 1, when we assume the coefficients to be real. Then The function Φ d (z) = bc 1 z −a + bc z−a + 2d will be minimal on the unit circle at z = 1 or z = −1. Thus positivity will be lost, when either a reaches 1 or −1, or bc − (a − 1)d = 0 or bc − (a + 1)d = 0. This is exactly the condition also reflected in the determinant of W (x c ) at the analytic center x a . This again illustrates that the X-passivity radius at the analytic center gives a good bound the passivity radius of the system.

Concluding remarks
We have derived conditions for the analytic center of the linear matrix inequalities (LMIs) associated with the passivity of linear continuous-time or discrete-time systems. We have presented numerical methods to compute these analytic centers with steepest ascent and Newton-like methods and we have presented lower bounds for the passivity radii associated with the LMIs evaluated at the respective analytic center.

A Derivatives of functions of complex matrices
In this appendix we present a precise derivation of the formulas for the differentiation of a matrix function with respect to a complex matrix. Here we distinguish between complex vector spaces C n and the corresponding real vector space R n + ıR n . Both spaces can be identified by c : C n → R n + ıR n , c(v) = (v) + ı (v). For matrix spaces of dimension m × n we use the usual identification with the vector spaces C n and R n + ıR n . The space C n is equipped with the standard scalar product x, y C := x H y. By ∂ ∂X we denote the differentiation in a real vector space, whereas the differentiation of a holomorphic function g is denoted by g . Note that if we write c • g(x) = u(x r + ıx i ) + ıv(x r + ıx i ), then by the Cauchy-Riemann equations, see e. g. [9], we have c • g (x) = ∂ ∂xr u(x r + ıx i ) − ı ∂ ∂x i u(x r + ıx i ). Then we have the following result: Lemma A.1. Assume that g : C n×n → C is holomorphic. Then f : R n×n + ıR n×n → R defined by f (X r + ıX i ) := g(X) is differentiable over R with and ∂ ∂X f (X r + ıX i ), ∆ R = g (X), c −1 (∆) C , ∆ = ∆ r + ı∆ i .
For the holomorphic function g(X) = det(X) the following fact is well-known, see e. g. [14] for a proof in the real case, that easily extends to the complex case.
Applying the chain-rule we finally obtain the differentiation formula, which is used throughout this paper.
Corollary A.1. Let f : R n×n + ıR n×n → R with f (X r + ıX i ) = ln det(X) and X ∈ C n×n with det(X) > 0. Then

B Differences between continuous-time and discrete-time systems
Usually, statements for a continuous linear time-invariant system can be transformed back and forth to discrete-time systems using some bilinear transform. However, the equations determining the analytic center in both cases are cubic in X, which suggests that there might not be a one-to-one correspondence. We have shown that the eigenvalues of the feedback system matrix A Fc at the analytic center lie on the imaginary axis in the continuous-time case, whereas they lie inside the unit disk in the discrete-time setting. In this appendix we show that it is indeed necessary to consider the continuous-time and discrete-time case separately by showing that the three equations determining the analytic center are not preserved under the usual bilinear transformations.
where A F d denotes the bilinear transform of the matrix A Fc = A c − B c F c . Thus, the transformed feedback matrix F d can be defined by It needs to be analyzed how the Riccati operator P c (X) is transformed for a fixed X. Clearly, then X fulfills the Riccati equation Ricc c (X) = −Q c with Q c := −P c (X). From the equations Ricc c (X) = −Q c and Ricc d (X) = −Q d one would then expect, that P d (X) = −Q d . However, we have the relation where we compute We thus obtain that which, by considering that P c 0 and equation (23), only coincides with −Q d if F c = 0. Thus we have shown, that if we enforce a feedback, that keeps the feedback system matrix A F d on the unit circle, then the transformed residual of the Riccati operator P c does not correspond to the discrete-time residual P d . In other words, since relation (25) has to hold, the transformation of the feedback (24) cannot be true, and thus the discrete-time feedback system matrix A F d does not lie on the unit circle. Indeed, as mentioned before, the eigenvalues lie strictly inside the unit circle.