Port-Hamiltonian ( pH ) systems are network-based models that arise within a modeling framework in which a physical system model is decomposed into submodels that are interconnected principally through the exchange of energy. Submodels may themselves be decomposed into submodels and at some level of branching they typically reflect one of a variety of core modeling paradigms that describe phenomenological aspects of the dynamics having different physical character, such as e.g. electrical, thermodynamic, or mechanical. The pH framework is able to knit together submodels featuring dramatically different physics through disciplined focus on energy flux as a principal mode of system interconnection; pH structure is inherited via power conserving interconnection, and a variety of physical properties and functional constraints (e.g., passivity and energy and momentum conservation) are encoded directly into the structure of the model equations BeaMXZ18 ; MehM19_ppt ; Sch13 . Interconnection of submodels often creates further constraints on system behavior and evolution, originating as conservation laws (e.g., Kirchoff’s laws or mass balance), or as position and velocity limitations in mechanical systems. As a result system models are often naturally posed as combinations of dynamical system equations and algebraic constraint equations, see e.g., BreCP96 ; KunM06 ; Ria08 . Port-Hamiltonian descriptor systems and port-Hamiltonian differential-algebraic equations (pHDAE) are alternate terms that refer to the same class of dynamical system models that aggregate pH dynamical system models together with algebraic constraints. These are among the most effective system representations for network-based models of multi-physics problems.
When a pHDAE system is linearized around a stationary solution, one obtains a linear time-invariant pHDAE with specially structured coefficient matrices, see BeaMXZ18 . Our focus in this work is on the construction of high fidelity reduced models for such pHDAE systems. Model reduction is important for such systems since they are often components in larger pHDAE systems and preserving the associated interconnection structure for energy flux through subsystems in such a way that retains in the reduced system pH properties as well as algebraic side constraints may be critical in order to retain overall simulation quality.
Although the approach we develop here can be extended easily to more general settings, we narrow our focus to a particular formulation as one of the simpler among alternative formulations laid out in BeaMXZ18 ; MehMW18 .
A linear time-invariant DAE system of the form
with , , , , , , on a compact interval is a pHDAE system if the following properties are satisfied:
The differential-algebraic operator
is skew-adjoint, i.e., we have that and ,
is positive semidefinite, i.e., , and
the passivity matrix
is symmetric positive semi-definite, i.e., .
The associated quadratic Hamiltonian function of the system is
In Definition 1, the Hessian matrix of is the energy matrix, . is the dissipation matrix; is the structure matrix describing the energy flux among internal energy storage elements; and are port matrices describing how energy enters and leaves the system. and are matrices associated with a direct feed-through from input to output . If the system has a (classical) solution in when presented with a given input function , then the dissipation inequality
must hold, i.e., the system is both passive and Lyapunov stable, as defines a Lyapunov function, see BeaMXZ18 .
In practice, network-based automated modeling tools such as modelica111https://www.modelica.org/, Matlab/Simulink222http://www.mathworks.com, or 20-sim333https://www.20sim.com each may produce system models that are very often overdetermined as a consequence of redundant modeling, a situation that is often difficult to avoid. Thus, they contain either explicit or hidden algebraic constraints and so, the automated generation of a system model usually must be followed by a reformulation and regularization procedure that makes the system model compatible with standard simulation, control and optimization tools. When the system is pH, this regularization step should additionally respect the pH structure.
Note that system models may be very large and very complex, e.g., models that arise from semi-discretized (in space) continuum models in hydrodynamics EggKLMM18 ; HeiSS08 ; Sty06 , or mechanics GraMQSW16 ; MehS05 . In such cases, model reduction techniques are critical to enable application of control and optimization methods. Preservation of both pH and constraint structure is necessary to maintain model integrity.
For linear time-invariant pH systems with strictly positive-definite Hamiltonians (i.e., positive-definite ), structure-preserving model reduction methods are well-developed, including tangential interpolation approaches GugPBS12 ; GugPBS09
, moment matchingPolS10 ; PolS11 ; WolLEK10 as well as effort and flow constraint reduction methods PolS12 . Interpolatory approaches have been extended to nonlinear pH systems as well, in BG2011pHstructPresMOR and ChaBG16 .
We have special interest in the case of singular , for which only recently in EggKLMM18 ; HauMM19_ppt have approaches been proposed extending model reduction techniques that preserve pH structure to this case. In order to obtain physically meaningful results it is essential to preserve both explicit and hidden constraints; see GugSW13 ; HeiRA11 ; MehS05 ; Sty06 for the case of general unstructured DAE systems.
We discuss model reduction methods incorporating both regularization and interpolation, for linear pHDAE systems having the form (1), taking advantage and respecting the structure. We present a systematic treatment of different model classes and describe how in these different formulations constraints are represented in the transfer function, and so how they can be preserved with interpolatory model reduction methods. The key step identifies, as in BeaMXZ18 ; MehM19_ppt , all the redundancies as well as both explicit and implicit system constraints, partitioning the system equations into redundant, algebraic, and dynamic parts. Only the dynamic part is then reduced in a way that preserves structure.
We proceed in §2 to recall basic properties of general DAE systems which will be important for later developments. In §3, we restrict consideration to pHDAEs and note important simplifications of the general regularization procedure to pHDAEs of the form (1). Structure preserving model reduction of pHDAEs using interpolatory methods is discussed first in general terms in §4, and then in more detail for specific model structures that include semi-explicit index- pHDAE systems (§4.1), semi-explicit pHDAE systems with index- constraints (§4.2), and finally, semi-explicit pHDAE systems with a combination of index- and index- constraints (4.3). We propose an algorithm to generate effective interpolation data in §5.1 followed by a few numerical examples that demonstrate its effectiveness.
2 General Differential-Algebraic Systems
In this section we recall some basic properties of general linear constant coefficient DAE systems
where , , , ; for details, see, e.g., BunBMN99 .
The matrix pencil is said to be regular if for some . For regular pencils, the finite eigenvalues are the values for which . If the reversed pencil
has the eigenvaluethen this is the infinite eigenvalue of .
With zero initial conditions as in (3) and a regular pencil , we obtain, in the frequency domain, the rational transfer function
which maps Laplace transforms of the input functions to the Laplace transforms of the corresponding output functions . If the transfer function is written as
where is a proper rational matrix function and is a polynomial matrix function, then the finite eigenvalues are the poles of the proper rational part, while infinite eigenvalues belong to the polynomial part.
In the following, a matrix with orthonormal columns spanning the right nullspace of the matrix is denoted by and a matrix with orthonormal columns spanning the left nullspace of by . These matrices are not uniquely determined although the spaces are. For ease of notation, we speak of these matrices as the corresponding spaces.
Regular pencils can be analyzed via the Weierstraß Canonical Form; see, e.g., Gan59b .
If is a regular pencil, then there exist nonsingular matrices and for which
where is a matrix in Jordan canonical form whose diagonal elements are the finite eigenvalues of the pencil and is a nilpotent matrix, also in Jordan form. and are unique up to permutation of Jordan blocks.
The index of the pencil is the index of nilpotency of the matrix in (5), and if is nonsingular, then the pencil has index zero. A simple characterization of a pencil to have index less than or equal to one is that is square and nonsingular.
2.1 Controllability and observability
To analyze general DAE systems and also to understand the properties of the transfer function, some controllability and observability conditions are needed, see e.g. BunBMN99 ; Dai89 . Consider the conditions
(C1) characterizes controllability of the dynamical system, while systems that satisfy condition (C2) are controllable at infinity or impulse controllable Dai89 . If both (C1) and (C2) hold, the system is strongly controllableBunBMN99 . Analogous observability conditions dual to (C1) and (C2) appear as
Condition (O1) characterizes observability of the dynamical system, and systems that satisfy both (O1) and (O2) are called strongly observable. Condition (O2) is observability at infinity or impulse observability. If a system is strongly controllable and strongly observable then the underlying dynamical system is minimal.
(C1) and (C2) are preserved under non-singular equivalence transformations as well as under state and output feedback, i.e., if the system satisfies (C1) or (C2), then for any non-singular , , and any and , the system satisfies the same condition for all of the following three choices:
Analogous properties hold for (O1) and (O2).
Note, however, that regularity or non-regularity of the pencil, the index, and the polynomial part of the transfer function are in general not preserved under state or output feedback. For systems that satisfy (C2), there exists a suitable linear state feedback matrix that such that is regular and of index at most one. Also if conditions (C2) and (O2) hold, then there exists a linear output feedback matrix so that the pencil has this property, see BunBMN99 . In practical control design it is therefore common, see KunM06 , to first perform a regularization procedure that we will briefly discuss in the next subsection.
In general, it cannot be guaranteed that a system, generated either from a realization procedure or an automated modeling procedure, has a regular pencil . Therefore, typically a regularization procedure has to be applied, see BinMMS15_ppt ; CamKM12 ; KunM06 . To briefly describe this procedure, we write the system in behavior form combining input and state to a descriptor vector
descriptor vector, i. e.,
with block matrices (having block rows/columns) of the form
It follows from the general theory of DAE (see KunM06 ) that for the system (12) there exists an integer such that the coefficients of the derivative array (13) associated with satisfy the following three properties.
, i. e., there exists a matrix of size and maximal rank satisfying .
, i. e., can be partitioned as , with of size and of size , such that has full row rank and . Furthermore, there exists a matrix of size and maximal rank satisfying .
, i. e., there exists a matrix of size and maximal rank satisfying with .
Furthermore, system (12) has the same solution set as the system
where , and .
System (14) is a reformulation of (12) (using the original model and its derivatives) and it does not alter the solution set. Since the constructed submatrices and are obtained from the block matrix by transformations from the left, it follows (see KunMR01 ) that derivatives of the input function are not needed, and the partitioning into parts associated with original states and original controls is not mixed. Including the output equation, we thus obtain a reformulated system of the form
where with , and with .
Next we focus on pHDAEs and show how structure helps simplify the analysis.
3 Linear port-Hamiltonian systems
In this section we analyze and simplify the regularization procedure of §2.2 for pHDAEs of the form (1) by judiciously exploiting structure. The pencil associated with a free dissipative Hamiltonian DAE system (i.e., ) has many nice properties MehMW18 : The index of the system is at most two; all eigenvalues are in the closed left half plane; and all eigenvalues on the imaginary axis are semi-simple (except for possibly the eigenvalue , which may have Jordan blocks of size at most two). Furthermore, a singular pencil can only occur when have a common nullspace; see MehMW19_ppt . Therefore, if one is able to efficiently compute this common nullspace, it is possible to remove the singular part via a congruence transformation.
For the pHDAE in (1), there exists an orthogonal basis transformation matrix such that in the new variable , the system has the form
where is a regular pencil and has full row rank. Also, the subsystem
that is obtained by removing the third equation and the variable is still a pHDAE.
First determine an orthogonal matrix(via SVD) such that
Such exists, since have a common nullspace when the pencil is singular MehMW19_ppt . Then a row compression of via an orthogonal matrix and a congruence transformation with is performed, so that with , we obtain the zero pattern in (16). Updating the output equation accordingly and using the fact that the transformed passivity matrix
is still semidefinite; it follows that and , giving the desired form.
Lemma 1 shows that for pHDAEs the general regularization procedure can be significantly simplified. The next result presents a condensed form and shows that the controllability conditions (C2) and (O2) are equivalent and hold for (18)-(19).
Next compute a spectral decomposition , where . Then an appropriate congruence transformation with yields
Finally performing a spectral decomposition , with nonsingular, and an appropriate congruence transformation with yields the desired form together with the updating of the output equation accordingly. The fact that the blocks , , blocks do not appear follows again from the semidefiniteness of the transformed passivity matrix. By the symmetry structure it follows directly that the condition (C2) holds if and only (O2) holds. Since , , and is invertible, this holds if and only if the matrix has full rank, which is the case because (which is the block in (18)-(19)) has full row rank and the matrices , , were assumed not to have a left nullspace anymore; thus has full row rank.
Note that even though the procedure to compute the condensed form (16)-(17) was used to prove that the controllability conditions (C2) and (O2) hold, it also immediately separates the dynamical part (given by the first block row) the algebraic index- conditions (with and without dissipation, given by the second and third block rows), and the index- conditions (given by the fourth block row). The last row is the singular part of the free system. However, since the conditions (C2) and (O2) hold, the system can be made regular by output feedback. Note that this is already displayed in the subsystem (18)-(19), so the regularization can be made immediately after creating (18)-(19). Let us denote, for ease of notation, this subsystem by
If we apply an output feedback with , that makes the (2,2) block in the closed loop system invertible, i.e., it turns the singular part into an index- equation. This corresponds to a feedback for the state-to-output map in (18) of the form , which we can insert into the first equation to obtain
The negative of the symmetric part of the system matrix is the Schur complement of the matrix
Hence the closed loop system is regular and the closed loop system is still pH.
The regularization procedures described above are computationally demanding, since they typically require large-scale singular value decompositions or spectral decompositions. Fortunately, in many practical cases the condensed form is already available directly from the modeling procedure, so that the transfer function can be formed and the model reduction method can be directly applied. In the following sections we assume that this is the case and extend the tangential interpolation model reduction procedure ofGugSW13 to three classes of semi-explicit systems, namely systems with only index- equations, systems with only index- equations, and systems with both index- and index- equations.
4 Interpolatory model reduction of pHDAEs
Given an order- pHDAE as in (1), we want to construct an order- reduced pHDAE, with , having the same structured form
such that , , , , , satisfy the same requirements as in Definition 1 and that the output of (22) is an accurate approximation to the original output over a wide range of admissible inputs . We will enforce accuracy by constructing the reduced model (22) via interpolation.
Let and denote the transfer functions of (1) and (22) where , , , , and similarly for the reduced-order (“hat”) quantities. Given the right-interpolation points together with the corresponding left-tangent directions and the left-interpolation points together with the corresponding right-tangent directions , we would like to construct such that it tangentially interpolates , i.e.,
Then the interpolatory reduced model can be obtained via projection:
In the setting of pHDAEs two fundamental issues arise: First, the reduced quantities in (26) are no longer guaranteed to have the pH structure. This is easiest to see in the reduced quantity . If we decompose into its symmetric and skew-symmetric parts, we can no longer guarantee that the symmetric part is positive semi-definite. This could be resolved by using a Galerkin projection, i.e., with . In this case one only satisfies the interpolation conditions associated with right interpolation data. However, this does not resolve the second issue since in the generic case when , the reduced quantity is expected to be a nonsingular matrix; thus the reduced system will be an ODE. This means that the polynomial parts of and do not match, leading to unbounded errors.
Structure-preservation interpolatory reduction of pH systems in the most general setting of tangential interpolation has been studied in GugPBS12 ; GugPBS09 . However, this work focused on the ODE case. On the other hand, GugSW13 developed the tangential interpolation framework for reducing unstructured DAEs with guaranteed polynomial matching. Only recently in EggKLMM18 ; HauMM19_ppt , the combined problem has been investigated. Our goal in the rest of the paper is to develop a complete treatment of structure-preserving interpolatory model reduction problem for index- and index- pHDAEs in the general setting of tangential interpolation.
In the rest of this section, in some instances we will use the generic transfer function notations and for full and reduced pHDAEs. It will be clear from the context what the matrices , , , , , and their reduced counter-parts are in terms of the pH structure.
4.1 Semi-explicit index- pHDAE systems
The simplest class of pHDAEs are semi-explicit index- pHDAEs of the form
where and are nonsingular. We have the following interpolation result.
Consider the full-order pHDAE system in (27). Let the interpolation points and the corresponding tangent directions be given. Construct the interpolatory model reduction basis as
where is partitioned conformably with the system, and define the matrices
Let partitioned accordingly to (27). Then the transfer function of the reduced model
matches the polynomial part of and tangentially interpolates it, i.e.,
Define , and decompose , into their symmetric and skew-symmtric part. Then, the reduced model (29) is a pHDAE system if the reduced passivity matrix is positive seimdefinite.
We employ a Galerkin projection using the interpolatory model reduction basis to obtain the intermediate reduced model
Even though this reduced model is a pHDAE system due to the one-sided Galerkin model reduction and it satisfies the tangential interpolation conditions, it will not match the transfer function at