Numerical methods to compute a minimal realization of a port-Hamiltonian system

TL;DR

This paper introduces three numerical algorithms for deriving port-Hamiltonian models from time-domain input-output data, enhancing the tools available for modeling complex multi-physics systems.

Contribution

It presents novel indirect numerical methods for realizing port-Hamiltonian systems from data, with implementations and numerical examples demonstrating their effectiveness.

Findings

Algorithms successfully derive pH models from data

Methods are implemented in MATLAB

Numerical examples validate the approaches

Abstract

Port-Hamiltonian (pH) systems are a very important modeling tool in almost all areas of systems and control, in particular in network based model of multi-physics multi-scale systems. They lead to remarkably robust models that can be easily interconnected. This paper discusses the derivation of pH models from time-domain input-output data. While a direct construction of pH models is still an open problem, we present three different indirect numerical methods for the realization of pH systems. The algorithms are implemented in MATLAB and their performance is illustrated via several numerical examples.

Figures

Tables

Table 1. Table 1: Computation time of the algorithms for the presented examples

Examples	Simple method	PRBT	Nearest
Illustrative example	0.82441	3.45088	2.56567
Ladder network	Failed (time$>3600s$)	719.74701	10.26562

Equations

$$\dot{x}$$

$$y$$

$$W=\left[\begin{array}[]{lc}R&P\\[2.84526pt] P^{T}&S\end{array}\right]$$

$$E\dot{x}$$

$$y$$

$$W=\left[\begin{array}[]{lc}Q^{T}RQ&Q^{T}P\\[2.84526pt] P^{T}Q&S\end{array}\right]\geq 0.$$

$$\displaystyle E=\begin{bmatrix}M_{1}&0&0\\ 0&M_{2}&0\\ 0&0&0\end{bmatrix},\ A=\begin{bmatrix}0&G&0\\ -G^{\top}&D&-N^{\top}\\ 0&N&0\end{bmatrix},\ B=\begin{bmatrix}0\\ B_{2}\\ 0\end{bmatrix},\$$

$$\begin{array}[]{rcl}E\dot{\hat{x}}&=&A\hat{x}+B\hat{u},\quad\hat{x}(0)=0,\\ \hat{y}&=&C\hat{x}+D\hat{u}.\end{array}$$

$$E_1\dot{x}$$

$$0$$

$$y$$

$$A_{2}x(0)+B_{2}u(0)=0.$$

$$\dot{z}$$

$$y$$

$$\displaystyle\left[\begin{array}[]{cc}A_{11}&A_{12}\end{array}\right]$$

$${\tilde{B}}_1$$

$$\tilde{A}$$

$$\tilde{C}$$

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

$$\dot{\hat{z}}$$

$$y$$

$$\left[\begin{array}{cc}{\hat{A}}^T\hat{Q}+\hat{Q}\hat{A}& \hat{Q}\hat{B}-{\hat{C}}^T\\ {\hat{B}}^T\hat{Q}-\hat{C}& -(\hat{D}+{\hat{D}}^T)\end{array}\right]\le 0.$$

$$U_{2}^{T}(D+D^{T})U_{2}=\left[\begin{array}[]{cc}\hat{D}_{1}+\hat{D}_{1}^{T}&0\\ 0&0\end{array}\right],\ \hat{B}U_{2}=\left[\begin{array}[]{cc}\hat{B}_{1}&\hat{B}_{2}\end{array}\right],\ U_{2}^{T}\hat{C}=\left[\begin{array}[]{c}\hat{C}_{1}\\ \hat{C}_{2}\end{array}\right],$$

$$\left[\begin{array}{cc}{\hat{A}}^T\hat{Q}+\hat{Q}\hat{A}& \hat{Q}{\hat{B}}_1-{{\hat{C}}_1}^T\\ {{\hat{B}}_1}^T\hat{Q}-{\hat{C}}_1& -({\hat{D}}_1+{{\hat{D}}_1}^T)\end{array}\right]\le 0$$

$$Q$$

$$F$$

$$S$$

$$AX{E}^T+EX{A}^T$$

$$EX{C}^T-P_{l}B$$

$$A^{T}YE+E^{T}YA$$

$$E^{T}YB-{P_r}^T{C}^T$$

$$\left[\hat{E},\hat{A},\hat{B},\hat{C},\hat{D}\right]=\left[W_{b}E{T}_b,W_{b}A{T}_b,W_{b}B,C{T}_b,D\right]$$

$$A^{T}YE+E^{T}YA+{(B^{T}YE-C{P}_r)}^T{(D+D^T)}^{-1}(B^{T}YE-C{P}_r)=0,$$

$$Y=P_l^{T}Y{P}_l.$$

$$\underset{(\tilde{E},\tilde{A},\tilde{B},\tilde{C},\tilde{D})\in \mathcal{D}}{\inf }\mathcal{F}(\tilde{E},\tilde{A},\tilde{B},\tilde{C},\tilde{D})$$

$$\begin{gathered} \mathcal{F}(\tilde{E},\tilde{A},\tilde{B},\tilde{C},\tilde{D})={\Vert A-\tilde{A}\Vert }_F^2+{\Vert B-\tilde{B}\Vert }_F^2+{\Vert C-\tilde{C}\Vert }_F^2 \\ +{\Vert D-\tilde{D}\Vert }_F^2+{\Vert E-\tilde{E}\Vert }_F^2. \end{gathered}$$

$$\underset{M,J,R,F,P,S}{\mathop{\inf}}\,\mathcal{G}(M,J,R,F,P,S)\text{ subject to }{{J}^{T}}=-J,\ M\geq 0,\ \left[\begin{matrix}R&P\\ {{P}^{T}}&S\\ \end{matrix}\right]\geq 0,\$$

Full text

Keywords: port-Hamiltonian system, positive-real system, realization, dissipation inequality.

AMS Subject Classification: 93B15, 93B20, 93C05

1 Introduction

In this paper we consider the construction of minimal realizations of linear time-invariant (LTI) port-Hamiltonian (pH) systems from input-output data in the time domain. Such pH systems are structured representations of dynamical systems [1, 45, 29, 35, 46, 51, 52, 56] that typically arise from energy based modeling via bond graphs [12, 20]. For several application domains this modeling framework has been implemented in the package, 20-sim (http://www.20sim.com/).

The pH framework is extremely powerful, since it encodes underlying physical principles such as conservation laws directly into the structure of the system model. Standard LTI pH systems take the form

[TABLE]

Here $J=-J^{T}\in\mathbb{R}^{n,n}$, $R=R^{T}\in\mathbb{R}^{n,n}$ is positive semidefinite; $F\pm P\in\mathbb{R}^{n,m}$ are port matrices; $S+N$, with $S=S^{T}\in\mathbb{R}^{m,m}$ and $N=-N^{T}\in\mathbb{R}^{m,m}$ is the feed-through from input to output, and $Q=Q^{T}\in\mathbb{R}^{n,n}$ positive semidefinite is associated with a quadratic energy functional (Hamiltonian) ${\mathcal{H}}(x)=\frac{1}{2}x^{T}Qx$. There is an intricate and subtle relationship between pH systems, positive real systems and passive systems, see e.g. [4, 15]. To assure that an LTI system of the form (1) is Lyapunov stable and passive one further requires that the block matrix

[TABLE]

associated with the passivity of the system, is also positive semidefinite (which we denote as $W\geq 0$ or $W>0$ in the positive definite case). The pH structure has many further spectral properties [43], it is very robust under structured perturbations [41, 42], the family of pH systems is closed under power-conserving interconnection [37], and model reduction of pH systems via Galerkin projection yields (smaller) pH systems [3, 31, 48]. It also has recently been shown how to robustly transform a passive and Lyapunov stable standard LTI system into a robust pH system [4].

Another important topic that has recently been addressed is the incorporation of constraints in the system models. Then one obtains pH differential-algebraic equations (pHDAEs), or pH descriptor systems, see [6, 53], which in the LTI case take the form

[TABLE]

In the descriptor case it is allowed that $E,Q\in\mathbb{R}^{n,\ell}$ are rectangular, provided that $Q^{T}E=E^{T}Q\geq 0$ and it is required that $Q^{T}JQ=-Q^{T}J^{T}Q\in\mathbb{R}^{n,n}$, $Q^{T}RQ=Q^{T}R^{T}Q\in\mathbb{R}^{n,n}$, that the quadratic Hamiltonian is ${\mathcal{H}}(x)=\frac{1}{2}x^{T}E^{T}Qx$ and the passivity matrix satisfies

[TABLE]

It has been shown that the discussed properties of standard pH systems carry over to pHDAEs in [6, 32, 53, 54].

Considering the remarkable system properties and the fact that in many applications a detailed first principle physical modeling of a complicated technical system may be too demanding, it is natural to study the question of how to construct realizations of pH (or pHDAE) systems from input-output data in the time domain. In particular, it is of great importance to develop numerical methods for the construction of such realizations of minimal order, usually called minimal realizations. Most current pH modeling approaches work directly on the physical network realization, see e.g. [56, 60], while direct data based approaches so far are rare, and direct pH modeling from time-domain data is an open problem.

An example, where a data based approach is urgently required, is the modeling of gas transport networks as pHDAE systems. A general hierarchy of first principles based models has been derived in [23]. For the flow in the pipe network, after discontinuous Galerkin discretization, see [24], the model takes the pHDAE form (3) with $Q=I$, $P=0$, $S+N=0$,

[TABLE]

where $M_{1}=M_{1}^{T},M_{2}=M_{2}^{T}>0$ are mass matrices and $x_{1}$, $x_{2}$ are coefficient vectors associated with pressure and mass flux, respectively, and $x_{3}$ represents Lagrange multipliers associated with conservation laws for mass and momentum at network nodes. The matrix $B_{2}$ is associated with the transport in and out of the pipe network and can be used for coupling with other network components, in particular the compressor station for which no first principles model, but huge real world data sets of time domain input-output data are available. In the context of rewriting the whole gas transport system as a network of hierarchies of pH submodels, it would be ideal to generate a model component in pH form.

Another example, where PH modeling from data would be important, is in the modeling of a cable-driven parallel robot (CDPR) [55]. PH modeling in this case is used to design physics-shaping controllers, and the PH model was derived analytically. The analysis in [55] shows the limitations and the complexity of trying to derive large and complex nonlinear PH models analytically. Similar requirements arise e.g. in complex nonlinear pH modeling for power electronics [21] and continuous stirred tank reactors [33].

In this paper we present three different numerical methods for generating such structured (minimal) pH system realizations. All these require first an interpolation step where the input-output sequence is interpolated to obtain a descriptor (DAE) system. The resulting realizations from the interpolation step are also usually not well formulated for systems and control applications. For this reason they are followed by a regularization procedure, see [11, 16], which also guarantees that the structural properties are satisfied.

In the first method, a straightforward approach based on a combination of well-known classical computational control techniques is presented. A model is generated from classical unstructured realization techniques; this model is converted to a standard minimal state space system and then a linear matrix inequality (LMI) is solved. Using the solution of the LMI, a state space transformation is determined that takes the system to standard pH form.

The second method also uses a combination of standard techniques but replaces the transformation to a minimal system and the solution of the LMI by a positive real balanced truncation.

The third method follows a suggestion by [28] that was introduced for dissipative Hamiltonian systems and tackles the minimal pH realization problem as an optimization problem by computing the nearest minimal pH system to a given unstructured realization.

All three approaches are implemented in MATLAB. We compare the three approaches, and illustrate their performance at hand of several numerical examples.

2 Realization and regularization of descriptor systems

In this section we recall the basis for realizing and regularizing a linear time-invariant descriptor (DAE) system

[TABLE]

In a realization process, if the given input-output data are sampled frequency response data, then the Loewner approach [2, 40] is currently one of the most promising and best analyzed method. If the given data is in terms of input-output time-domain data, the case that we study here, then a modified Loewner approach has recently been proposed in [47]. This method generates an LTI descriptor system of the form (5) and is particularly useful if the system is given as a black box routine that simulates the model from chosen inputs or if the input-output data are just given as a time series.

To construct a realization, we assume that inputs $u$ and the resulting outputs $y$ of the system are evaluated or measured at a sequence of time points $t_{i}$ and recorded as $k$ pairs $(u_{i},y_{i})$ of vectors $u_{i},y_{i}\in\mathbb{R}^{m}$, $i=1,\ldots,k$, where $k$ is selected to ensure that the outputs have typically entered a steady state. Note that for a pH realization we would need the same number of inputs and outputs. The choice of the input sequence $u_{i}$, the time interpolation points $t_{i}$, and the number of samples $k$ are very flexible and can be adapted to best correspond to the real system at hand.

The resulting model realization is typically a linear descriptor system of the form (5), see [2, 40, 47]. However, these models are often not so easy to use with current numerical methods for analysis and control design. Often the resulting systems are even close to singular systems, i.e. systems that are not uniquely solvable for a given input $u$. A descriptor system of the form (5) is called non-regular (singular) if $E,A$ are not square matrices or in the square case if $\det(sE-A)\equiv 0$. To check regularity and to obtain a robustly regular system, a regularization procedure as those introduced in [10, 14, 16, 34, 38], and implemented for LTI descriptor systems in [11], needs to be applied. The regularization procedure of [11, 16] combines a behavioral approach and appropriate feedbacks. We will not present this procedure here, see [11] for a detailed description and a MATLAB implementation. As final output one obtains a reformulated system of the form

[TABLE]

where $E_{1},A_{1}\in\mathbb{R}^{d,n}$, $E_{2},A_{2}\in\mathbb{R}^{a,n}$, with $n=a+d$, $B_{1}\in\mathbb{R}^{d,m}$, $B_{2}\in\mathbb{R}^{a,m}$, $C\in\mathbb{R}^{m,n}$, $D\in\mathbb{R}^{m,m}$ and the matrix \left[\begin{array}[]{c}E_{1}\\ A_{2}\end{array}\right] is invertible. This property of (6) is usually referred to as the reformulated system being regular index at most one [13] or strangeness-free [38] as a free system for $u=0$.

In the form (6) consistency of initial conditions can be easily checked because these have to satisfy

[TABLE]

Note, however, that usually the resulting regularized DAE system will not have a pH structure, so in the next section we will present three procedures to obtain pHDAE realizations.

3 PH realization via an unstructured approach

In this section we present a straightforward approach to create a pH realization from a regularized unstructured realization of the form (6). It transforms (6) to a standard state-space form, then to a minimal realization and then this to a pH formulation.

Given a system of the form (6), the first step computes matrices $\tilde{A},\tilde{B},\tilde{C},\tilde{D}$ representing the same input-output map as (6) but in standard state space form

[TABLE]

Such a procedure has recently been suggested in [25]. However, since we have already regularized the descriptor system to the form (6), this procedure can be simplified and made more numerically robust. As a first step (since by construction $E_{1}$ has full row rank, we perform a singular value decomposition, see e.g. [30], E_{1}=U_{1}\left[\begin{array}[]{cc}\Sigma_{E}&0\end{array}\right]V^{T}, with orthogonal matrices $U_{1},V$, and $\Sigma_{E}\in\mathbb{R}^{d,d}$. We form

[TABLE]

partitioned accordingly. As a direct consequence of the regularization procedure it follows that $A_{22}$ is invertible. Multiplying the first block row by $\Sigma_{E}^{-1}$ and performing a Schur complement with respect to the matrix $A_{22}$, we obtain a standard state-space system of the form (8) with

[TABLE]

together with a constraint equation $x_{2}=-A_{22}^{-1}(A_{21}x_{1}+B_{2}u)$ which fixes $x_{2}$ in terms of $x_{1}$ and $u$. For further construction we only consider the system in $x_{1}$, since it represents the transfer function in frequency domain given by

[TABLE]

But we should keep in mind that the constraint equation cannot be ignored completely, since it restricts the initial conditions for $x,u$ in the behavior setting. Note that even if the original system did not have a feed-through term, the algebraic constraints will lead to such a feed-through term.

The presented reduction procedure is computationally acceptable if the condition numbers, see e.g. [30], of $\Sigma_{E}$ and $A_{22}$ are small, in relation to the modeling, interpolation, and round-off errors in the regularized system (6). If these condition numbers are too large, then this system is close to a singular or high index system and should be treated with great care.

The resulting system in standard state space form (8) may not be of minimal order. Using the well-known result, see e.g. [17], that a standard state space system of the form (8) with a proper frequency domain transfer function $\tilde{G}(s)$ is minimal if and only if $(\tilde{A},\tilde{B})$ is controllable and $(\tilde{A},\tilde{C})$ is observable, one can make the system minimal by removing the uncontrollable and observable parts, e.g. by the techniques in [17, 18, 22, 36, 50], in particular, using the numerically stable staircase form of [57] and its implementation in the SLICOT library [7, 44].

Then, keeping only the part that is controllable and observable, results in a minimal realization

[TABLE]

with a minimal proper transfer function $\hat{G}(s)=\hat{D}+\hat{C}(sI-\hat{A})^{-1}\hat{B}$. Note that by this minimalization procedure the state-space dimension may have been further reduced.

Each step of this four-stage process of realization, regularization, transformation to standard state-space form, and making the system minimal is well understood from its numerical properties, and also concerning its difficulties and robustness under perturbations, thus we have refrained here from presenting more details than necessary.

Using the well-known relation between pH systems and systems that are (asymptotically) stable and (strictly) passive, [58, 59], see also [4] for a discussion of some limiting cases, the transformation to pH form can then be determined by computing a positive definite solution $\hat{Q}$ to the linear matrix inequality (LMI)

[TABLE]

See also [5] for a more detailed discussion on numerical properties and a direct computation of this solution via a structured eigenvalue problem.

Note the LMI (12) in general may not have a positive definite solution if $D+D^{T}$ is singular. A commonly used practice is to use a perturbation $\hat{D}+\hat{D}^{T}+\varepsilon I$ to make this term positive definite. However, there is no perturbation analysis of the effect of this perturbation on the solution procedure as a whole. Another, more reliable approach is to first perform a projection of the problem, by computing an orthogonal matrix $U_{2}$ such that

[TABLE]

where $\hat{D}_{1}+\hat{D}_{1}^{T}>0$. Then one has to compute a positive definite solution of the LMI

[TABLE]

satisfying the constraint $\hat{Q}\hat{B}_{2}=\hat{C}_{2}^{T}$. According to the results of [58, 59], if the system is stable and passive then such a $Q=Q^{T}>0$ exists. If not, then a further regularization procedure is needed that can be based on [5], however, in the following we assume that such a $\hat{Q}>0$ can be computed. Then one can determine a Cholesky factorization $\hat{Q}=T{T}^{T}$, with a nonsingular matrix $T$. This Cholesky factor can be used to transform (11) to pH form as in [4, 5] leading to the system matrices

[TABLE]

We have implemented this procedure in MATLAB using the LMI approach as well as a Riccati equation approach [4], see Figure 1 for a block diagram of the procedure. To avoid unnecessary computation, we allow three different types of starting data for the resulting MATLAB function, either input/output data, or a general descriptor system or a standard state space system. Depending on these data some blocks of the procedure may be skipped.

Numerical results obtained with this procedure are presented in Section 6.

In this section we have discussed a five step procedure realization, regularization, transformation to standard state-space form, making the system minimal, transformation to PH form to generate a standard state-space pH realization from given time domain input-output data. It should be noted that although every single step in the discussed procedure is well understood from the point of view of error analysis and perturbation theory, the procedure may face a significant number of difficulties. These include a singularity or close to singularity of the system (5), an ill-conditioning of the matrices $\Sigma_{E}$, $A_{22}$, a close to non-controllability or non-observability of the minimal system (11), or that this system is not stable or not passive. If the system is not stable, then this can be fixed by an appropriate feedback, since the system is controllable. If the system is not passive, then a perturbation to the nearest passive system can be performed. But even if the system is passive it may happen that the Cholesky factor $T$ is ill-conditioned.

It is almost impossible to carry out a detailed error and perturbation analysis for the complete procedure, since there are too many factors that can go wrong in the procedure. Furthermore, since each of the steps builds upon the next, this simple straightforward approach will not always yield satisfactory results. Finally, as we will see below, it is computationally demanding and does not scale well to large scale problems.

4 Transformation based on positive real balanced truncation

Another straightforward approach to generate a minimal pH realization from an unstructured standard realization is to exploit the relationship between pH systems and positive real systems. Consider a DAE system of the form (5) with a transfer function $G(s)=D+C(sE-A)^{-1}B=p(s)+r(s)$, where $p=\sum_{i=0}^{\nu}M_{i}s^{i}$ with ${{M}_{i}}\in{{\mathbb{R}}^{m\times m}}$, is the polynomial and $r$ the proper rational part of $G$. Let $P_{r}$ and $P_{l}$ be the spectral projectors onto the right and left deflating subspaces of $\lambda E-A$ corresponding to the finite eigenvalues. It has been shown in [9] that if the system (5) is $R$-minimal, i.e. the dynamical part of the system is controllable and observable and if the transfer function $G$ is passive, then the projected Lur’e equation

[TABLE]

is solvable for ${{K}_{C}}\in{{\mathbb{R}}^{n,m}},{{J}_{c}}\in{{\mathbb{R}}^{m,m}}$ and $X\in{{\mathbb{R}}^{n,n}}$. Furthermore, if ${M}_{k}=0$ for $k>1$, ${M}_{1}={M}_{1}^{T}\geq 0$ and the projected Lur’e equations have a solution, then $G$ is positive real.

When, starting from a system of the form (5), we have performed the regularization procedure and obtained the system (6), then we have ${{M}_{0}}=D$ and $M_{\ell}=0$ for $\ell>0$. Then, if the resulting system is passive, the solution of (15) can be computed, and also instead of (15), its dual form can be solved equivalently.

[TABLE]

for ${{K}_{O}}\in{{\mathbb{R}}^{m,n}},{{J}_{O}}\in{{\mathbb{R}}^{m,m}}$ and $Y\in{{\mathbb{R}}^{n,n}}$ .

The R-minimality condition in the first part of Theorem 3.1 can be weakened. For standard systems, less restrictive conditions on the involved matrices are given in [19, 39].

The solution of the Lur’e equations can then be used to compute the balanced truncation, where the minimal solutions of the projected Lur’e equations, ${{X}_{\min}}$ and ${{Y}_{\min}}$, respectively, are the positive real controllability and observability Gramians. The system is called positive real balanced if ${{X}_{\min}}={{Y}_{\min}}=\mbox{\rm diag}(\Pi,0)$ with $\Pi=\mbox{\rm diag}({{\pi}_{1}},...,{{\pi}_{{{n}_{f}}}})$, where the ${{\pi}_{j}}$ are the positive real characteristic values of the system, i.e. the square roots of the non-zero eigenvalues of the matrix ${{X}_{\min}}{{E}^{T}}{{Y}_{\min}}E$.

After computing the balancing transformation, one can compute a transformation

[TABLE]

such that a balanced truncation is obtained and states corresponding to insignificant ${{\pi}_{j}}$ are truncated, see [49]. If $D+{{D}^{T}}$ is nonsingular, then the positive real Lur’e system is equivalent to the projected Riccati equation

[TABLE]

If $D+{{D}^{T}}$ is rank deficient, then actually the Riccati equation is not well defined and as in the previous section no positive definite solution exists that can be used for the transformation. Again commonly a perturbation ($\varepsilon I$ is used to make $D+{{D}^{T}}+\varepsilon I>0$. This method has been implemented in the MORLAB toolbox [8], where the truncation method is applied using Newton-type methods, where Lyapunov equations are solved in each step using the matrix sign function method. For illustration, the Positive Real Balanced truncation (PRBT) method is presented in Figure 2.

Again, although there are much fewer steps in this procedure, it is very difficult to carry out a detailed error and perturbation analysis for the complete procedure, since still many steps can go wrong in this procedure. However, the procedure scales much better than that of the previous section for large scale problems.

5 Nearest Port Hamiltonian system

Having discussed in the previous section two procedures for the construction of a pHODE realization of a descriptor system arising from a given input-output sequence and having also remarked on the potential numerical difficulties in computing the desired pHODE realization, in this section we proceed in a completely different way. Many of the difficulties discussed in the previous two approaches arise from the fact that the interpolation procedure may lead to a system being a (close to) singular, unstable or non-passive descriptor system even though the system from which the data were generated was originally stable and passive. Thus we may just ask the question whether we can just solve an optimization problem to compute the nearest stable, passive and regular and index at most one (as a free system) phDAE system. Let us denote the set of such systems by $\mathcal{D}$. For a system with given realization $(E,A,B,C,D)$ we would then have to solve the minimization problem

[TABLE]

where

[TABLE]

This a very difficult optimization problem, since the feasible set $\mathcal{D}$ is unbounded, highly non-convex, and neither open nor closed, see [27, 28] for a detailed analysis. However, requiring the structure of pHDAE form comes to help, because we can reformulate the problem as the optimization problem

[TABLE]

where

[TABLE]

For a slightly different version of this optimization problem, see [27, 28], where an algorithm to solve this optimization problem is based on a fast projected gradient method (FGM) which has the advantage of being in general much faster than the standard projected gradient method even in the nonconvex case. The algorithm presented in [28] also includes a restarting procedure which is necessary in our case, since the objective function is not convex and hence FGM without restart is not guaranteed to converge [26]. It is also possible to add weights to each term of the objective function.

Choosing good initial points is crucial to obtain good solutions for nonconvex optimization problems and this is still an open problem. In [28], an LMI based initialization procedure as discussed in our first approach is proposed. This LMI-based initialization works well when the initial system is close to being passive. However, solving the LMI may require a large computational effort if the system is large scale.

6 Numerical examples

In this section we present some numerical examples to illustrate the performance of the different methods. All tests were performed on a laptop Intel Core i7-5600U CPU @2.6GHz, 16GB RAM. The computation times of the algorithms are measured in seconds.

Consider a descriptor system (5) of order $5$ with

[TABLE]

Applying the first method results in a minimal pHODE system with matrices $E=Q=I$, $S=N=0$,

[TABLE]

of order $4$.

Applying the PRBT method results in a pHODE system with $E=Q=I$, $N=0$, $S=0.2204$,

[TABLE]

Finally, using the nearest pHDAE algorithm results in a pHODE system with $S=0.2204$, $N=0$,

[TABLE]

Note that $Q=E^{-1}$ is not the identity for the nearest pHDAE algorithm. This constitutes an alternative to the previous methods, but the resulting system has order equal to the original system, and it is necessary to employ a structure preserving minimalization procedure if one wants a minimal pHDAE system.

A comparison of the transfer functions of the three pH realizations with the original system is shown in the Bode plots in Figure 3.

As second example we consider an electrical ladder network consisting of resistors, capacitors and inductors. This example has order $2000$. It is tested for all the presented algorithms. The computation time for the straightforward approach exceeds the allowed time of $3600s$, this is why it not included in the bode plot comparison in Figure 4. The resulting pHODE system is obtained without losing the frequency response of the system. Also, in the case of the PRBT approach, the resulting system has order $22$ which is the numerically minimal order. This is a huge improvement compared to the given system which has order $2000$.

In terms of computation time, the difference between the different methods are large. Table 1 shows the computation times for the previous examples in seconds $(s)$.

The straightforward simple method presented very good results for the simple illustrative example of order $5$ but does not scale well for larger systems. The nearest pH algorithm gave the best results in terms of computation time but suffers from the problems mentioned before as it depends on the initialization and it can get stuck in a local minimum. Also it does not lead to a minimal system. The positive real balanced truncation method, although requiring more time than the nearest pH is very robust, and also results in a reduced model.

7 Conclusion

In this paper, we presented three methods to obtain a port-Hamiltonian system from time domain input-output data that are realized as a descriptor system or a standard state space representation. The first approach is a five step procedure that may have many numerical difficulties, is hard to analyze, and does not scale well in computation time. The second method is based on positive real balanced truncation and finally the third method solves an optimization problem to compute the nearest port-Hamiltonian system. The three procedures were implemented in MATLAB. This implementation and testing indicates that the method based on the positive real balanced truncation is the most robust especially for large scale systems. However, when initiated well, the nearest port-Hamiltonian system algorithm is the most promising but needs further implementation improvements. It should in particular be combined with a model reduction step as in [32] and a minimalization step.

Ideally one would like to construct a pH representation directly from the given input-output data, but this is still an open problem.