Structure-preserving Interpolatory Model Reduction for Port-Hamiltonian Differential-Algebraic Systems

Structure-pr eserving Interpolatory Model R eduction f or P ort-Hamiltonian Differential- Algebraic Sy stems C. A. Beattie, S. Gugercin, and V . Mehrmann Dedicated to Thanos Ant oulas on the occasion of his 70th birthday Abstract W e e xamine interpolator y model reduction methods that are w ell-suited f or treating larg e scale por t-Hamiltonian differential-alg ebraic sys tems in a w a y that is able to preserve and tak e advantag e of the underlying structural f eatures of the sy stem. W e introduce approaches that incor porate regular ization tog ether with a prudent selection of interpolation data. W e f ocus on linear time-inv ariant sy stems and present a sy stematic treatment of a variety of model classes that include combinations of inde x- 1 and index- 2 systems, descr ibing in par ticular ho w constraints ma y be represented in the transfer function so that the polynomial par t can be preser v ed with interpolator y methods. W e propose an algorithm to generate effectiv e inter polation data and illustrate its effectiveness on a numer ical e xample. 1 Introduction Port-Hamiltonian ( pH ) sys tems are network -based models that ar ise within a mod- eling frame w ork in which a phy sical sy stem model is decomposed into hierarchies of submodels that are interconnected principally through the e x chang e of energy . At some lev el the submodels typically reflect one of a variety of core modeling paradigms that describe phenomenological aspects of the dynamics ha ving differ - ent phy sical character , such as e.g. electrical, thermodynamic, or mechanical. The pH frame w ork is able to knit tog ether submodels f eaturing dramatically different ph ysics through a disciplined f ocus on ener gy flux as a principal mode of sy stem interconnection; pH str ucture is inher ited via po w er conser ving interconnection, and a variety of phy sical proper ties and functional constraints (e.g., passivity and energy C. A. Beattie Department of Mathematics, Virginia T ech, Blacksburg, V A 24061-0123, USA, e-mail: beattie@vt.edu S. Gugercin Department of Mathematics, Virginia T ech, Blacksburg, V A 24061-0123, USA, e-mail: gugercin@vt.edu V . Mehrmann Institut für Mathematik, T echnische Univ ersität Berlin, D-10829, Berlin, Ger many e-mail: mehr mann@math.tu-berlin.de 1 2 C. A. Beattie, S. Gugercin, and V . Mehr mann and momentum conservation) are encoded directl y into the structure of the model equations [2, 30]. Interconnection of submodels often creates further cons traints on sys tem behavior and ev olution, or iginating as conser vation la ws (e.g., Kirchoff’ s la ws or mass balance), or as position and v elocity limitations in mechanical systems. As a result sys tem models are often naturally posed as combinations of dynami- cal sy stem equations and algebraic constraint equations, i.e. as P or t-Hamiltonian descriptor syst ems or port-Hamiltonian differ ential-alg ebraic equations ( pHDAE ). When a pHD AE sy stem is linear ized around a stationary solution, one obtains a linear time-in v ariant pHDAE with specially s tructured coefficient matr ices, see [2]. Although the approach we dev elop here can be extended easily to more general settings, we nar ro w our f ocus to a particular formulation as one of the simpler among alternative formulations laid out in [2, 23]. Definition 1. A linear time-inv ar iant DAE sys tem of the f orm E ¤ x = ( J − R ) x + ( B − P ) u , x ( 𝑡 0 ) = 0 , y = ( B + P ) 𝑇 x + ( S + N ) u , (1) with E , J , R ∈ I R 𝑛 × 𝑛 , B , P ∈ I R 𝑛 × 𝑚 , S = S 𝑇 , N = − N 𝑇 ∈ I R 𝑚 × 𝑚 , on a compact interval I ⊂ I R is a pHD AE sys tem if the f ollo wing proper ties are satisfied: 1. The differential-algebraic operator E 𝑑 𝑑 𝑡 − J : C 1 ( I , I R 𝑛 ) → C 0 ( I , R 𝑛 ) is skew-adjoint , i.e., w e ha v e that J 𝑇 = − J and E = E 𝑇 , 2. E is positiv e semidefinite, i.e., E ≥ 0 , and 3. the passivity matr ix W =  R P P 𝑇 S  ∈ R ( 𝑛 + 𝑚 ) × ( 𝑛 + 𝑚 ) is symmetric positive semi-definite, i.e., W = W 𝑇 ≥ 0 . The associated quadratic Hamiltonian function H : R 𝑛 → R of the system is H ( x ) = 1 2 x 𝑇 Ex . In Definition 1, the Hessian matr ix of H ( x ) is the matr ix E which is also referred to as energy matrix ; R is the dissipation matrix ; J is the structure matrix descr ibing the energy flux among inter nal energy storag e elements; and B ± P are por t matrices describing how energy enters and leav es the sys tem. S and N are matr ices associated with a direct f eed-thr ough from input u to output y . If the sys tem has a solution x ∈ C 1 ( I , R 𝑛 ) in I when f or a giv en input function u , then the dissipation inequality 𝑑 𝑑 𝑡 H ( x ) = u 𝑇 y −  x u  𝑇 W  x u  Structure-preserving Inter polatory Model Reduction for pHDAE s 3 must hold, i.e., the sys tem is both passive and Ly apunov stable , as H defines a L yapuno v function , see [2]. In practice, netw ork -based automated modeling tools such as modelica , Ma tlab/Simulink , or 20-sim 1 each may produce sy s tem models that are o verde- termined as a consequence of redundant modeling. Thus, the automated generation of a system model usually must be f ollo w ed by a ref ormulation and regularization procedure (that preserves the pH structure and makes the sy s tem model compatible with standard simulation, control and optimization tools. pHD AE models ma y be very larg e and comple x, e.g., models that ar ise from semi-discretized (in space) continuum models in hydrodynamics [11, 19, 29], or mechanics [25]. In such cases, model reduction techniques are necessary to apply control and optimization methods. Preser vation of the pH , the constraint, and the interconnection str ucture is necessar y to maintain model integr ity . For liner time -in variant pH sys tems with positiv e-definite E , such methods are well-de v eloped, including tangential inter polation approaches [15, 16], moment matching [26, 28, 31] as well as effor t and flo w constraint reduction methods [27]. Inter polator y approaches ha v e been e xtended to nonlinear pH sy stems as well, in [4] and [9]. In the case of singular E , only recently in [11, 18] first model reduction techniques ha v e been introduced that preserv e the pH structure, while for unstructured D AEs such constraint preserving methods hav e been introduced in [17, 19, 25, 29]. W e discuss str ucture preser ving model reduction methods incor porating both regularization and inter polation f or linear pHD AE systems ha ving the f orm (1). F or different model classes we describe how cons traints are represented in the transf er function, and how the pol ynomial par t can be preser v ed with inter polatory model reduction methods. The ke y step identifies, as in [2, 24], all the redundancies as w ell as both e xplicit and implicit sy s tem constraints, partitioning the sys tem equations into redundant, algebraic, and dynamic par ts. Only the dynamic part is then reduced in a wa y that preser v es the structure. In §3, w e consider pHDAE s of the form (1) and note impor tant simplifications of the general regular ization procedure. Structure preser ving model reduction of pHD AE s using interpolator y methods is discussed first in g eneral ter ms in §4, and then in more detail for se v eral specific model structures. W e propose an algor ithm to g enerate effectiv e interpolation data in §5 f ollow ed by a numer ical e xample that demonstrates its effectiveness. 2 General Differential- Algebraic Sy stems In this section w e recall some basic properties of g eneral linear constant coefficient D AE sy stems E ¤ x = Ax + Bu , x ( 𝑡 0 ) = 0 , y = Cx + Du , (2) 1 https://www.modelica.org/ , http://www.mathworks.com , https://www.20sim.com 4 C. A. Beattie, S. Gugercin, and V . Mehr mann where E , A ∈ I R 𝑛 × 𝑛 , B ∈ I R 𝑛 × 𝑚 , C ∈ I R 𝑝 × 𝑛 , D ∈ I R 𝑝 × 𝑚 ; f or details, see, e.g., [7]. The matr ix pencil 𝜆 E − A is said to be regular if det ( 𝜆 E − A ) ≠ 0 f or some 𝜆 ∈ I C . For regular pencils, the finite eig envalues are the values 𝜆 ∈ I C for which det ( 𝜆 E − A ) = 0 . If the rev er sed pencil 𝜆 A − E has the eigen value 0 then this is the infinite eig envalue of 𝜆 E − A . With zero initial conditions x ( 𝑡 0 ) = 0 as in (2) and a regular pencil 𝜆 E − A , w e obtain, in the frequency domain, the rational transf er function H ( 𝑠 ) = C ( 𝑠 E − A ) − 1 B + D , which maps the Laplace transf orm of the input functions u to the Laplace transf orm of the output function y . If the transfer function is wr itten as H ( 𝑠 ) = G ( 𝑠 ) + P ( 𝑠 ) , where G ( 𝑠 ) is a proper rational matr ix function and P ( 𝑠 ) is a pol ynomial matr ix function, then the finite eigen v alues are the poles of the proper rational part, while infinite eigen values belong to the polynomial part. In the f ollo wing, a matr ix with or thonormal columns spanning the r ight nullspace of the matrix M is denoted b y 𝑆 ∞ ( M ) and a matr ix with or thonormal columns spanning the left nullspace of M b y 𝑇 ∞ ( M ) . Regular pencils can be anal yzed via the W eierstraß Canonical Form; see, e.g., [13]. The index 𝜈 of the pencil 𝜆 E − A is the size of the larg est bloc k associated with the eigen v alue ∞ in the W eierstraß canonical form and if E is nonsingular , then the pencil has inde x 0 . T o analyze general D AE systems and also to understand the proper ties of the transf er function, some controllability and observability conditions are needed, see e.g. [7, 10]. A system that satisfies C1 : rank [ 𝜆 E − A , B ] = 𝑛 f or all 𝜆 ∈ I C and C2 : rank [ E , A 𝑆 ∞ ( E ) , B ] = 𝑛 is called strong ly controllable [7]. If a sy stem satisfies analogous observability conditions O1 : rank  𝜆 E 𝑇 − A 𝑇 C 𝑇  = 𝑛 for all 𝜆 ∈ I C , and O2 : rank  E 𝑇 A 𝑇 𝑇 ∞ ( E ) C 𝑇  = 𝑛 , then it is called str ong ly obser vable . If a sys tem is strongly controllable and strongly obser vable then it is called minimal . Conditions (C1) and (C2) are preserved under non-singular equivalence transfor - mations as w ell as under state and output f eedback. N ote, how e v er , that regularity or non-regularity of the pencil, the inde x, and the polynomial par t of the transf er function are in general not preserv ed under state or output f eedback. For sy stems that satisfy (C2), there exis ts a suitable linear state feedbac k matr ix F 1 that such that 𝜆 E − ( A + BF 1 ) ) is regular and of inde x at most one. Also if conditions (C2) and (O2) hold, then there e xists a linear output feedbac k matr ix F 2 so that the pencil 𝜆 E − ( A + BF 2 C ) has this proper ty , see [7]. 3 Regularization of PHD AE Systems In general, it cannot be guaranteed that a system, g enerated either from a realization procedure or an automated modeling procedure, has a regular pencil 𝜆 E − A . There- f ore, typically a regularization procedure has to be applied, see [6, 8, 20]. For pHD AE s Structure-preserving Inter polatory Model Reduction for pHDAE s 5 the structure helps simplify theses g eneral regularization procedures significantl y , since the pencil 𝜆 E − ( J − R ) associated with a free dissipativ e Hamiltonian DAE sys tem (i.e., u = 0 ) has many nice properties [23]: The inde x of the system is at most tw o; all eigen values are in the closed left half plane; and all eig en values on the imaginary axis are semi-simple (e x cept for possibly the eigen v alue 0 , which may ha v e Jordan blocks of size at most two). Further more, a singular pencil can only occur when E , J , R ha v e a common nullspace; see [22]. Theref ore, if one is able to efficiently compute this common nullspace, it is possible to remov e the singular part. Lemma 1. F or the pHD AE in (1) , ther e exis ts an orthogonal basis tr ansf ormation matrix V ∈ I R 𝑛 × 𝑛 suc h that in the new variable e x =  e x 𝑇 1 e x 𝑇 2 e x 𝑇 3  𝑇 = V 𝑇 x , the sys tem has the form       E 1 0 0 0 0 0 0 0 0             ¤ e x 1 ¤ e x 2 ¤ e x 3       =       J 1 − R 1 0 0 0 0 0 0 0 0             e x 1 e x 2 e x 3       +       B 1 − P 1 B 2 0       u , (3) y =  ( B 1 + P 1 ) 𝑇 B 𝑇 2 0        e x 1 e x 2 e x 3       + ( S + N ) u , wher e 𝜆 E 1 − ( J 1 − R 1 ) is regular and B 2 has full ro w rank. Also, the subsyst em  E 1 0 0 0   ¤ e x 1 ¤ e x 2  =  J 1 − R 1 0 0 0   e x 1 e x 2  +  B 1 − P 1 B 2  u , (4) y =  ( B 1 + P 1 ) 𝑇 B 𝑇 2   e x 1 e x 2  + ( S + N ) u , that is obtained by remo ving the thir d equation and the v ariable x 3 is still a pHD AE . Proof. F irst deter mine an or thogonal matrix V 1 (via SVD) such that V 𝑇 1 ( 𝜆 E − ( J − R ) ) V 1 = 𝜆  E 1 0 0 0  −  J 1 − R 1 0 0 0  , V 𝑇 1 ( B − P ) =  B 1 − P 1 ˜ B 2 − ˜ P 2  . Such V 1 e xists, since E , J , R hav e a common nullspace when the pencil 𝜆 E − ( J − R ) is singular [22]. Then a ro w compression of ˜ B 2 − ˜ P 2 via an orthogonal matrix ˜ V 2 and a congr uence transformation with V 2 = diag ( I , ˜ V 2 ) is per f ormed, so that with V = diag ( V 1 , V 2 ) , w e obtain the zero patter n in (3). U pdating the output eq uation accordingly and using the fact that the transformed passivity matrix ˜ W =  V 𝑇 R V V 𝑇 P P 𝑇 V 𝑇 S  ∈ R ( 𝑛 + 𝑚 ) × ( 𝑛 + 𝑚 ) is still semidefinite; it follo ws that P 2 = 0 and P 3 = 0 , giving the desired form. The ne xt result presents a condensed f orm and sho ws that conditions (C2) and (O2) are equiv alent and hold f or a subsystem of (4). 6 C. A. Beattie, S. Gugercin, and V . Mehr mann Lemma 2. F or the pHD AE in (4) there exists an orthogonal basis transf ormation b V in the state space and U in the control space such that in the new variables b x =  b x 𝑇 1 b x 𝑇 2 b x 𝑇 3 b x 𝑇 4 b x 𝑇 5 b x 𝑇 6  𝑇 = b V 𝑇  e x 𝑇 1 e x 𝑇 2  𝑇 , and  u 𝑇 1 u 𝑇 2 u 𝑇 3  𝑇 = U 𝑇 u the sys tem has the form             E 11 E 12 0 0 0 0 E 21 E 22 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0                        ¤ b x 1 ¤ b x 2 ¤ b x 3 ¤ b x 4 ¤ b x 5            =             J 11 − R 11 J 12 − R 12 J 13 − R 13 J 14 J 15 0 J 21 − R 21 J 22 − R 22 J 23 − R 23 J 24 0 0 J 31 − R 21 J 32 − R 32 J 33 − R 33 0 0 0 J 41 J 42 0 0 0 0 J 51 0 0 0 0 0 0 0 0 0 0 0                         b x 1 b x 2 b x 3 b x 4 b x 5 b x 6             +             B 11 − P 11 B 12 − P 12 B 13 − P 13 B 21 − P 21 B 22 − P 22 B 23 − P 23 B 31 − P 31 B 32 − P 32 B 33 − P 33 0 B 42 B 43 0 0 B 53 0 0 B 63                   u 1 u 2 u 3       , (5)       y 1 y 2 y 3       =       ( B 11 + P 11 ) 𝑇 ( B 21 + P 21 ) 𝑇 B 31 + P 31 ) 𝑇 0 0 0 ( B 12 + P 12 ) 𝑇 ( B 22 + P 22 ) 𝑇 B 32 + P 32 ) 𝑇 B 𝑇 42 0 0 ( B 13 + P 13 ) 𝑇 ( B 23 + P 23 ) 𝑇 B 33 + P 33 ) 𝑇 B 𝑇 43 B 𝑇 53 B 𝑇 63                   b x 1 b x 2 b x 3 b x 4 b x 5 b x 6             +       D 11 D 12 D 13 D 21 D 22 D 23 D 31 D 32 D 33             u 1 u 2 u 3       , wher e E 22 , J 33 − R 33 , J 15 , and B 42 and B 63 ar e invertible. F urthermore, the subsy stem obtained by dele ting the fir st, fifth, and sixt h bloc k r ow and column satisfies (C2) and equiv alently (O2) . Proof. The proof follo ws again by a sequence of or thogonal transf or mations. Star ting from (4) in the first step one deter mines an or thogonal matrix V 1 (via a spectral decomposition of E 1 ) such e V 𝑇 1 E 1 e V 1 =  E 11 0 0 0  with E 11 > 0 , and then f orms a congr uence transf ormation with b V 1 = diag ( e V 1 , I ) yielding b V 𝑇 1 ( J − R ) b V 1 =  ˜ J 11 − ˜ R 11 ˜ J 12 − ˜ R 12 ˜ J 21 − ˜ R 21 ˜ J 22 − ˜ R 22  , b V 𝑇 1 ( B − P ) =  ˜ B 1 − ˜ P 1 ˜ B 2 − ˜ P 2  . Ne xt compute a full rank decomposition e V 𝑇 2 ( e J 22 − e R 22 ) e V 2 =  b J 22 − b R 22 0 0 0  , where b J 22 − b R 22 is inv er tible and b R 22 ≥ 0 . This exis ts, since e J 22 − e R 22 has a positiv e semidefinite symmetr ic par t. Then an appropr iate cong ruence transf or mation with b V 2 = diag ( I , e V 2 , I ) yields Structure-preserving Inter polatory Model Reduction for pHDAE s 7 b V 𝑇 2 b V 𝑇 1 E b V 1 b V 2 =         E 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0         , b V 𝑇 2 b V 𝑇 1 ( J − R ) b V 1 b V 2 =         b J 11 − b R 11 b J 12 − b R 12 b J 13 0 b J 21 − b R 21 b J 22 − b R 22 0 0 b J 31 0 0 0 0 0 0 0         , b V 𝑇 2 b V 𝑇 1 ( B − P ) =          b B 1 − b P 1 b B 2 − b P 2 b B 3 b B 4          , where b J 22 − b R 22 is in v ertible and b B 4 has full ro w rank. Then one performs an orthogonal decomposition e V 𝑇 3  b B 3 b B 4  U =       0 B 42 B 43 0 0 B 53 0 0 B 63       with B 42 and B 63 square nonsingular, where the number of row s in B 63 is that of B 4 and applies an appropr iate congruence transformation with b V 3 = diag ( I , I , b V 3 , I ) so that one obtains block matr ices           E 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0           ,           ¯ J 11 − ¯ R 11 ¯ J 12 − ¯ R 12 ¯ J 13 − ¯ R 13 ¯ J 14 0 ¯ J 21 − ¯ R 21 ¯ J 22 − ¯ R 22 0 0 0 ¯ J 31 0 0 0 0 ¯ J 41 0 0 0 0 0 0 0 0 0           , As final step one computes a column compression of the full row rank matrix ¯ J 41 and applies a an appropriate cong ruence transformation. This yields the desired f orm. The f act that the se v eral P blocks do not occur f ollow s ag ain from the semidefi- niteness of the transf or med passivity matr ix. Since J 51 , B 42 and B 63 are in v ertible, it f ollow s immediately that u 3 = 0 and ˆ x 1 = 0 and that ˆ x 5 is uniq uely determined b y all the other variables. Consider ing the subsys tem obtained by remo ving the first, fifth, and sixth bloc k row and column, it f ollow s by the symmetry str ucture that the condition (C2) holds if and only (O2) holds. The procedure to compute the condensed f orm (3) immediately separates the dynamical part (given by the second block ro w), the algebraic inde x- 1 condition (giv en by the third bloc k row s), and the inde x- 2 conditions (giv en by the fifth and first bloc k ro w). In the follo wing, f or ease of notation, we denote the sy stem that is obtained b y remo ving the variables b x 1 , b x 4 and the cor responding first and fifth equation by E 𝑟 ¤ x 𝑟 = ( J 𝑟 − R 𝑟 ) x 𝑟 + ( B 𝑟 − P 𝑟 ) u , x 𝑟 ( 𝑡 0 ) = 0 , y 𝑟 = ( B 𝑟 + P 𝑟 ) 𝑇 x + ( S 𝑟 + N 𝑟 ) u . If w e apply an output feedbac k u = − K 𝑟 y 𝑟 with K 𝑟 = K 𝑇 𝑟 > 0 , that makes the (2,2) block in the closed loop system in v ertible, then this corresponds to a f eedback f or 8 C. A. Beattie, S. Gugercin, and V . Mehr mann the state-to-output map in (4) of the form y 𝑟 = ( I + ( S 𝑟 + N 𝑟 ) K 𝑟 ) − 1 ( ( B 𝑟 + P 𝑟 ) 𝑇 x 𝑟 , which we can inser t into the first equation to obtain E 𝑟 ¤ x 𝑟 = ( J 𝑟 − R 𝑟 ) x 𝑟 + ( B 𝑟 − P 𝑟 ) K 𝑟 y 𝑟 =  ( J 𝑟 − R 𝑟 ) − ( B 𝑟 − P 𝑟 ) ( K − 1 𝑟 + ( S 𝑟 + N 𝑟 ) ) − 1 ( ( B 𝑟 + P ) 𝑇  x 𝑟 . The negativ e of the symmetr ic par t of the system matr ix is the Schur complement of ˜ W 𝑟 =  R 𝑟 P 𝑟 P 𝑇 𝑟 K − 1 𝑟 + S 𝑟 + N 𝑟  ≥ 0 . Hence the closed loop sys tem is regular and the closed loop system is still pH . The procedures described abov e are computationally demanding, since they typ- ically require larg e-scale singular value decompositions or spectral decompositions. Fortunately , in many practical cases the condensed f orm is already a v ailable directly from the modeling procedure, so that the transf er function can be f ormed and the model reduction method can be directly applied. 4 Interpolatory model reduction of pHD AE s Giv en an order - 𝑛 pHD AE as in (1), we want to construct an order- 𝑟 reduced pHD AE , with 𝑟  𝑛 , ha ving the same structured f orm b E ¤ x 𝑟 = ( b J − b R ) x 𝑟 + ( b B − b P ) u , x 𝑟 ( 𝑡 0 ) = 0 , y 𝑟 = ( b B + b P ) 𝑇 x 𝑟 + ( b S + b N ) u , (6) such that b E , b J , b R ∈ I R 𝑟 × 𝑟 , b B , b P ∈ I R 𝑟 × 𝑚 , b S = b S 𝑇 , b N = − b N 𝑇 ∈ I R 𝑚 × 𝑚 satisfy the same requirements as in Definition 1 and that the output y 𝑟 ( 𝑡 ) of (6) is an accurate appro ximation to the or iginal output y ( 𝑡 ) ov er a wide range of admissible inputs u ( 𝑡 ) . W e will enf orce accuracy b y constructing the reduced model (6) via inter polation. Let H ( 𝑠 ) = C ( 𝑠 E − A ) − 1 B + D and b H ( 𝑠 ) = b C ( 𝑠 b E − b A ) − 1 b B + b D denote the transf er functions of (1) and (6), where C = ( B + P ) 𝑇 , A = J − R , B = B − P , D = S + N , and similar l y f or the reduced-order (“hat”) quantities. Given the right-inter polation points { 𝜎 1 , 𝜎 2 , . . . , 𝜎 𝑟 } ∈ I C tog ether with the corresponding left-tang ent directions { 𝒃 1 , 𝒃 2 , . . . , 𝒃 𝑟 } ∈ I C 𝑚 and the left-inter polation points { 𝜇 1 , 𝜇 2 , . . . , 𝜇 𝑟 } ∈ I C together with the cor responding right-tangent directions { ℓ 1 , ℓ 2 , . . . , ℓ 𝑟 } ∈ I C 𝑚 , we would like to construct b H ( 𝑠 ) suc h that it tang entially inter polates H ( 𝑠 ) , i.e., H ( 𝜎 𝑖 ) 𝒃 𝑖 = b H ( 𝜎 𝑖 ) 𝒃 𝑖 and ℓ 𝑇 𝑖 H ( 𝜇 𝑖 ) = ℓ 𝑇 𝑖 b H ( 𝜇 𝑖 ) , f or 𝑖 = 1 , 2 , . . . , 𝑟 . (7) These tangential inter polation conditions can be easil y enf orced via a Petro v -Galerkin projection [1, 5, 12]. Construct V ∈ I C 𝑛 × 𝑟 and Z ∈ I C 𝑛 × 𝑟 using Structure-preserving Inter polatory Model Reduction for pHDAE s 9 V =  ( 𝜎 1 E − A ) − 1 B 𝒃 1 , ( 𝜎 2 E − A ) − 1 B 𝒃 2 , · · · ( 𝜎 𝑟 E − A ) − 1 B 𝒃 𝑟  and (8) Z =  ( 𝜇 1 E − A ) − 𝑇 C 𝑇 ℓ 1 , ( 𝜇 2 E − A ) − 𝑇 C 𝑇 ℓ 2 , · · · ( 𝜇 𝑟 E − A ) − 1 C 𝑇 ℓ 𝑟  . (9) Then the interpolator y reduced model can be obtained via projection: b E = Z 𝑇 EV , b A = Z 𝑇 A V , b B = Z 𝑇 B , b C = C V , and b D = D . (10) In the setting of pHD AE s two fundamental issues ar ise: Firs t, the reduced quantities in (10) are no long er guaranteed to hav e the pH structure. This is easiest to see in the reduced quantity b A = Z 𝑇 A V = Z 𝑇 JV − Z 𝑇 R V . If w e decompose b A into its symmetric and sk ew -symmetr ic par ts, we can no long er guarantee that the symmetric part is positive semi-definite. This could be resol v ed by using a Galerkin projection, i.e., with Z = V . In this case one only satisfies the inter polation conditions associated with right interpolation data. Ho w e v er , this does not resolv e the second issue since in the g eneric case when 𝑟 < rank ( E ) , the reduced q uantity b E is e xpected to be a nonsingular matrix; thus the reduced sys tem will be an ODE . This means that the polynomial par ts of H ( 𝑠 ) and b H ( 𝑠 ) do not match, leading to unbounded errors. Structure-preservation inter polator y reduction of pH sys tems in the most general setting of tangential inter polation has been studied in [15, 16]. Ho w e v er , this w ork f o- cused on the ODE case. On the other hand, [17] de v eloped the tang ential inter polation frame w ork f or reducing unstructur ed D AE s with guaranteed pol ynomial matching. Only recently in [11, 18], the combined problem has been inv es tigated. W e now dev elop a treatment of structure-preserving inter polatory model reduction problem f or index- 1 and index- 2 pHD AE s in the general setting of tang ential inter polation. 4.1 Semi-explicit inde x- 1 pHD AE systems The simplest class of pHDAE s are semi-explicit index- 1 pHDAE s of the form  E 11 0 0 0  ¤ x ( 𝑡 ) =  J 11 − R 11 J 12 − R 12 − J 𝑇 12 − R 𝑇 12 J 22 − R 22  x ( 𝑡 ) +  B 1 − P 1 B 2 − P 2  𝑢 ( 𝑡 ) , y ( 𝑡 ) =  B 𝑇 1 + P 𝑇 1 B 𝑇 2 + P 𝑇 2  x ( 𝑡 ) + ( S + N ) u ( 𝑡 ) . (11) where E 11 and J 22 − R 22 are nonsingular . W e ha v e the follo wing inter polation result. Theorem 1. Consider the pHD AE sys tem in (11) . Let the interpolation points { 𝜎 1 , 𝜎 2 , . . . , 𝜎 𝑟 } ∈ I C and the corresponding tang ent dir ections { 𝒃 1 , 𝒃 2 , . . . , 𝒃 𝑟 } ∈ I C 𝑚 be giv en. Construct the interpolator y model reduction basis V as V =  V 1 V 2  =  ( 𝜎 1 E − A ) − 1 B 𝒃 1 , · · · , ( 𝜎 𝑟 E − A ) − 1 B 𝒃 𝑟  ∈ I C 𝑛 × 𝑟 , (12) wher e V is partitioned conformably with the syst em, and define the matrices 10 C. A. Beattie, S. Gugercin, and V . Mehr mann B =  𝒃 1 𝒃 2 · · · 𝒃 𝑟  ∈ I C 𝑚 × 𝑟 and D = D − ( B 𝑇 2 + P 𝑇 2 ) ( J 22 − R 22 ) − 1 ( B 2 − P 2 ) ∈ I C 𝑚 × 𝑚 . Let A =  A 11 A 12 A 21 A 22  = J − R , partitioned according ly to (11) . Then the transf er function b H ( 𝑠 ) of the r educed model b E ¤ x 𝑟 ( 𝑡 ) = ( b J − b R ) x 𝑟 ( 𝑡 ) + b B u ( 𝑡 ) , y 𝑟 ( 𝑡 ) = b C x 𝑟 ( 𝑡 ) + b Du ( 𝑡 ) (13) with b E = V 𝑇 1 E 11 V 1 b C = C V + ( B 𝑇 2 + P 𝑇 2 ) A − 1 22 ( B 2 − P 2 ) B , b A = V 𝑇 A V + B 𝑇 ( D − D ) B , b D = D = D − ( B 𝑇 2 + P 𝑇 2 ) A − 1 22 ( B 2 − P 2 ) b B = V 𝑇 B + B 𝑇 ( B 𝑇 2 + P 𝑇 2 ) A − 1 22 ( B 2 − P 2 ) , matc hes the polynomial par t of H ( 𝑠 ) and tang entially interpolates it, i.e., H ( 𝜎 𝑖 ) 𝒃 𝑖 = b H ( 𝜎 𝑖 ) 𝒃 𝑖 , for 𝑖 = 1 , 2 , . . . , 𝑟 . Define b P = 1 2 ( − b G + b C ) , and decompose b D = b S + b N , b A = b J − b R into their symmetric and skew-symmtric part. Then, the reduced model (13) is a pHD AE syst em if the r educed passivity matrix c W =  b R b P b P 𝑇 b S  is positiv e semidefinite. Proof. W e emplo y a Galerkin projection using the inter polator y model reduction basis V to obtain the intermediate reduced model e E = V 𝑇 1 E 11 V 1 , e A = e J − e R = V 𝑇 JV − V 𝑇 R V , e B = V 𝑇 B , e C = C V and e D = D . Ev en though this reduced model is a pHD AE sy stem due to the one-sided Galerkin model reduction and it satisfies the tangential interpolation conditions, it will not match the transfer H ( 𝑠 ) function at 𝑠 = ∞ , i..e, its pol ynomial par t, giv en b y lim 𝑠 →∞ H ( 𝑠 ) = D = D − ( B 𝑇 2 + P 𝑇 2 ) A − 1 22 ( B 2 − P 2 ) . A remedy to this problem, proposed in [3, 21] and emplo y ed in the general DAE setting in [17], is to modify the D -term in the reduced model to matc h the pol ynomial par t and at the same time to shift the other reduced quantities appropr iately so as to keep the tang ential inter polation property . Using this, we obtain a modified reduced model Structure-preserving Inter polatory Model Reduction for pHDAE s 11 b E = e E = V 𝑇 1 E 11 V 1 , b A = e A + B 𝑇 ( D − D ) B = V 𝑇 JV − V 𝑇 R V − B 𝑇 ( B 𝑇 2 + P 𝑇 2 ) 𝐴 − 1 22 ( B 2 − P 2 ) B , b B = e B + B 𝑇 ( B 𝑇 2 + P 𝑇 2 ) 𝐴 − 1 22 ( B 2 − P 2 ) = V 𝑇 1 ( B 1 − P 1 ) + V 𝑇 2 ( B 2 − P 2 ) + B 𝑇 ( B 𝑇 2 + P 𝑇 2 ) 𝐴 − 1 22 ( B 2 − P 2 ) , (14) b C = e C + ( B 𝑇 2 + P 𝑇 2 ) 𝐴 − 1 22 ( B 2 − P 2 ) B = ( B 𝑇 1 + P 𝑇 1 ) V 𝑇 1 + ( B 𝑇 2 + P 𝑇 2 ) V 𝑇 2 + ( B 𝑇 2 + P 𝑇 2 ) 𝐴 − 1 22 ( B 2 − P 2 ) B , and b D = D = D − ( B 𝑇 2 + P 𝑇 2 ) 𝐴 − 1 22 ( B 2 − P 2 ) , which satisfies the or iginal tangential inter polation conditions, and matches the polynomial par t due to the modified b D -term. Ho we ver , f or this system to be pH we need to chec k that the associated passivity matrix is still positive semidefinite, which after rewriting the in put and output matrix in the usual w a y is e xactly the condition on c W in the assertion. W e then hav e that the reduced model not onl y satisfies the interpolation conditions and matches the polynomial par t at 𝑠 = ∞ , but also is pH . Remar k 1. Note that if the input does not influence the algebraic equations, i.e., if B 2 = P 2 = 0 , then the shift of the constant ter m is not necessary , and the f ormulas simplify significantly , i.e., b A = V 𝑇 A V , b B = e B = V 𝑇 B , b C = C V , and b D = D = D . Another solution to preser ving pH structure via inter polation can be obtained through the f ollo wing theorem. Theorem 2. Consider a full-order pHD AE sys tem of the form (11) . Let inter - polation points { 𝜎 1 , 𝜎 2 , . . . , 𝜎 𝑟 } ∈ I C and the corresponding tang ent dir ections { 𝒃 1 , 𝒃 2 , . . . , 𝒃 𝑟 } ∈ I C 𝑚 be given. Construct the interpolatory model reduction ba- sis V as in (12) . Then the reduced model  b V 𝑇 1 E 11 b V 1 0 0 0  ¤ x ( 𝑡 ) =  V 𝑇 1 ( J 11 − R 11 ) V 1 V 𝑇 1 ( J 12 − R 12 ) ( − J 𝑇 12 − R 𝑇 12 ) V 1 J 22 − R 22  x ( 𝑡 ) +  V 𝑇 1 ( B 1 − P 1 ) B 2 − P 2  u ( 𝑡 ) y 𝑟 ( 𝑡 ) =  ( B 𝑇 1 + P 1 ) 𝑇 V 1 B 2 + P 𝑇 2  x ( 𝑡 ) + Du ( 𝑡 ) (15) r etains the pH structur e, tang entially interpolates t he original model, and matches the polynomial par t. Proof. W e first note that the subspace spanned b y the columns of V =  V 𝑇 1 V 𝑇 2  𝑇 is contained in the subspace spanned by the columns of b V : = diag ( V 1 , I ) . Then the sys tem in (15) results from reducing the or iginal system in (11) via b V . Since span ( V ) ⊆ span ( b V ) , this reduced DAE automatically satisfies the inter polation conditions and since b V does not alter the matr ix J 22 − R 22 and the matrices B 2 , P 2 , the polynomial par t of the transfer function is D − ( B 𝑇 2 + P 𝑇 2 ) ( J 22 − R 22 ) − 1 ( B 2 − P 2 ) , matching that of the or iginal model. The reduced sys tem in (15) is pH as this one- sided projection retains the original pH structure. Remar k 2. Theorem 2 presents a seemingly easy solution to structure-preser ving interpolator y model reduction of inde x-1 pHD AE s compared to Theorem 1. How ev er , 12 C. A. Beattie, S. Gugercin, and V . Mehr mann this may not be the maximal reduction that is possible because redundant algebraic conditions cannot be remo v ed; see [25]. 4.2 Semi-explicit pHD AE systems with index- 2 constraints As next class w e study semi-explicit index- 2 systems. W e first consider the case that the input does not affect the algebraic equations. Theorem 3. Consider an index- 2 pHD AE sys tem of the form  E 11 0 0 0   ¤ x 1 ( 𝑡 ) ¤ x 2 ( 𝑡 )  =  J 11 − R 11 J 12 − J 𝑇 12 0   x 1 ( 𝑡 ) x 2 ( 𝑡 )  +  B 1 − P 1 0  u ( 𝑡 ) y ( 𝑡 ) =  B 𝑇 1 + P 𝑇 1 0   x 1 ( 𝑡 ) x 2 ( 𝑡 )  + Du ( 𝑡 ) , (16) with E 11 > 0 and se t A 11 = J 11 − R 11 . Given interpolation points { 𝜎 1 , 𝜎 2 , . . . , 𝜎 𝑟 } and associated tang ent dir ections { 𝒃 1 , 𝒃 2 , . . . , 𝒃 𝑟 } , let the v ector s v 𝑖 , for 𝑖 = 1 , 2 , . . . , 𝑟 , be the firs t bloc k of the solution of  A 11 − 𝜎 𝑖 E 11 J 12 − J 𝑇 12 0   v 𝑖 z  =  ( B 1 − P 1 ) 𝒃 𝑖 0  . (17) Define V = [ v 1 , v 2 , . . . , v 𝑟 ] . Then the reduced model b E ¤ x 𝑟 = ( b J − b R ) x 𝑟 + b B u ( 𝑡 ) , y 𝑟 = b C x 𝑟 + b D , (18) with b E = V 𝑇 E 11 V , b J = V 𝑇 J 11 V , b R = V 𝑇 R 11 V , b B = V 𝑇 B 1 − V 𝑇 P 1 , b C = B 𝑇 1 V 𝑇 + P 𝑇 V 𝑇 1 , and b D = D , (19) is still pH , matches the polynomial part of the original transf er function, and satisfies the tangential interpolation conditions, i.e., H ( 𝜎 𝑖 ) 𝒃 𝑖 = b H ( 𝜎 𝑖 ) 𝒃 𝑖 , for 𝑖 = 1 , 2 , . . . , 𝑟 . Proof. N ote first that the regular ity of 𝜆 E − ( J − R ) and the index- 2 condition imply that − J 𝑇 12 E − 1 11 J 12 is in v ertible, see [2, 20]. Follo wing [17], we wr ite (16) as 𝚷 E 11 𝚷 ¤ x 1 ( 𝑡 ) = 𝚷 A 11 𝚷 x 1 ( 𝑡 ) + 𝚷 B 1 u ( 𝑡 ) , y ( 𝑡 ) = C 1 𝚷 x 1 ( 𝑡 ) + Du ( 𝑡 ) , (20) together with the algebraic equation x 2 ( 𝑡 ) = − ( J 𝑇 12 E − 1 11 J 12 ) − 1 J 𝑇 12 E − 1 11 A 11 x 1 ( 𝑡 ) − ( J 𝑇 12 E − 1 11 J 12 ) − 1 J 𝑇 12 E − 1 11 B 1 u ( 𝑡 ) , Structure-preserving Inter polatory Model Reduction for pHDAE s 13 with the projector 𝚷 = I − E − 1 11 J 12 ( J 𝑇 12 E − 1 11 J 12 ) − 1 J 𝑇 21 The equiv alent sys tem (20) is now an implicit ODE pH sys tem that can be reduced with standard model reduction techniques. How e v er , this would require computing the projector 𝚷 e xplicitl y , see [25]. F or g eneral index- 2 D AE sy stems one can av oid this computational step in inter polator y model reduction, see [17]. T o adapt this idea to pHD AE sys tems, w e construct V using (17) and then compute the reduced-order quantities via one-sided projection as in (19). This constr uction of V , as in [17], guarantees that the reduced model in (18) tang entially interpolates the original pHD AE sys tem in (20) and the polynomial par t of the transf er function in (16) is given by D = S + N par titioned in its symmetr ic and ske w-symmetric part. Since (18) is an implicit ODE pH system with the ex act D -term, it matches the polynomial par t of the original transfer function H ( 𝑠 ) . It remains to show that (18) is pH . By constr uction in (19), b J is sk e w-symmetric, b R is symmetric positive semidefinite, and E 11 is symmetric positive definite. Moreov er ,  V 𝑇 R 11 V V 𝑇 P 1 P 𝑇 1 V S  =  V 𝑇 0 0 I   R 11 P 1 P 𝑇 1 S   V 0 0 I  ≥ 0 , since the original model is PH, and therefore the pH-structure is retained. The situation becomes more complicated when the second block in B is nonzero, i.e., the sys tem has the form  E 11 0 0 0   ¤ x 1 ( 𝑡 ) ¤ x 2 ( 𝑡 )  =  J 11 − R 11 J 12 − J 𝑇 12 0   x 1 ( 𝑡 ) x 2 ( 𝑡 )  +  B 1 − P 1 B 2 − P 2  u ( 𝑡 ) , y ( 𝑡 ) =  B 𝑇 1 + P 𝑇 1 B 2 + P 𝑇 2   x 1 ( 𝑡 ) x 2 ( 𝑡 )  + Du ( 𝑡 ) . (21) W e ha v e the f ollo wing theorem. Theorem 4. Consider a pHD AE sys tem of the form (21) and define the matrices Z = − ( A 21 E − 1 11 A 12 ) − 1 =  J 𝑇 12 E − 1 11 J 12  − 1 C = ( B 𝑇 1 − P 𝑇 1 ) − ( B 𝑇 2 − P 𝑇 2 ) ZJ 𝑇 12 E − 1 11 ( J 11 − R 11 ) , B = ( B 1 − P 1 ) + ( J 11 − R 11 ) E − 1 11 J 12 Z ( B 2 − P 2 ) , D 0 = D − ( B 𝑇 2 + P 𝑇 2 ) ZJ 𝑇 12 E − 1 11 ( B 1 − P 1 ) , and D 1 = − ( B 𝑇 2 + P 𝑇 2 ) Z ( B 2 − P 2 ) . Giv en interpolation points { 𝜎 1 , 𝜎 2 , . . . , 𝜎 𝑟 } and associated tang ent directions { 𝒃 1 , 𝒃 2 , . . . , 𝒃 𝑟 } , let the v ectors v 𝑖 be the first bloc ks of the solutions of  J 11 − R 11 − 𝜎 𝑖 E 11 J 12 − J 𝑇 12 0   v 𝑖 z  = B 𝒃 𝑖 , for 𝑖 = 1 , 2 , . . . , 𝑟 , (22) 14 C. A. Beattie, S. Gugercin, and V . Mehr mann and set V : = [ v 1 , v 2 , . . . , v 𝑟 ] , u 1 : = u , u 2 : = ¤ u , and b D : =  D 0 D 1  = b S + b N . Then, the reduced model b E ¤ x 𝑟 = ( b J − b R ) x 𝑟 +  b B − b P 0   u 1 u 2  , y 𝑟 = ( b B + b P ) 𝑇 x 𝑟 + b D  u 1 u 2  , (23) with b E = V 𝑇 E 11 V , b J = V 𝑇 J 11 V , b R = V 𝑇 R 11 V , b B = 1 2  V 𝑇 B + V 𝑇 C 𝑇  , and b P = 1 2  V 𝑇 C 𝑇 − V 𝑇 B  , (24) satisfies the interpolation conditions, matches the polynomial part of the transf er function, and preserves the pH structur e, pr ovided that the reduced passivity matrix c W =  b R b P b P 𝑇 b S  is positiv e semidefinite. Proof. The proof follo ws similar to the proof of Theorem 3. Follo wing [17], the state x 1 can be decomposed as x 1 = x 𝑐 + x 𝑔 , where x 𝑔 = E − 1 11 A 12 ( J 𝑇 12 E − 1 11 J 12 ) − 1 ( B 2 − P 2 ) u ( 𝑡 ) and x 𝑐 ( 𝑡 ) satisfies J 𝑇 12 x 𝑐 = 0 . Then, one can rewrite (21) as 𝚷 E 11 𝚷 ¤ x 0 ( 𝑡 ) = 𝚷 A 11 𝚷 x 0 ( 𝑡 ) + 𝚷 B 1 u ( 𝑡 ) y ( 𝑡 ) = C 𝚷 x 1 ( 𝑡 ) + D 0 u ( 𝑡 ) + D 1 ¤ u ( 𝑡 ) . (25) As bef ore, the ODE part can be reduced with usual model reduction techniques. Follo wing [17], ho w e v er , we achiev e this without computing the projectors 𝚷 𝑙 and 𝚷 𝑙 e xplicitl y , instead b y constructing V using (22) and then applying one-sided model reduction with V to obtain the reduced model b E ¤ x 𝑟 = ( b J − b R ) x 𝑟 + e Bu ( 𝑡 ) , y 𝑟 = e C 𝑇 x 𝑟 + D 0 u ( 𝑡 ) + D 1 ¤ u ( 𝑡 ) , (26) where b E = V 𝑇 E 11 V , b J = V 𝑇 J 11 V , b R = V 𝑇 R 11 V , e B = V 𝑇 B , and e C = C V . This reduced model, by constr uction, satisfies the tangential inter polation conditions. Note that the reduced model in (26) has e xactl y the same realization as the reduced model in (23) e x cept for the reduced e B and e C ter ms. The reduced terms in (26) already hav e the desired pH structure. T o recov er the symmetry in (25), w e deternmine matrices b B and b P such that e B = V 𝑇 B = b B − b P and e C = C V = ( b B + b P ) 𝑇 via b B = 1 2  e B + e C 𝑇  = 1 2 V 𝑇  B + C 𝑇  and b P = 1 2  e C 𝑇 − e B  = 1 2 V 𝑇  C 𝑇 − B  , reco v ering (24). The final requirement to retain the pH structure is, then, again that  b R b P b P 𝑇 1 2  D + D 𝑇   ≥ 0 , which is the fina l condition in the statement of the theorem. This approach has the disadv antag e that one has to introduce the derivativ e of u as an extra input, which ma y lead to difficulties when applying standard control and optimization methods. For this reason it is usually pref erable to first perform an Structure-preserving Inter polatory Model Reduction for pHDAE s 15 inde x reduction via an appropr iate output f eedback, see Section 3, and then appl y the results from Section 4.1. But note that this chang es the polynomial par t of the transf er function. 4.3 Semi-explicit pHD AE systems with index- 1 and index- 2 constraints Finall y we consider semi-explicit index- 2 sy stems which also ha v e an inde x- 1 par t, [11, 18, 25]. In this case we only consider the special case with state       E 11 E 22 0 E 21 E 22 0 0 0 0             ¤ x 1 ¤ x 2 ¤ x 3       =       J 11 − R 11 J 12 − R 12 J 13 J 21 − R 21 J 22 − R 22 0 J 31 0 0             x 1 x 2 x 3       +       B 1 − P 1 B 2 − P 2 0       u , (27) output y =  ( B 1 + P 1 ) 𝑇 ( B 2 + P 2 ) 𝑇 0        x 1 x 2 x 3       + ( S + N ) u ,  E 11 E 22 E 21 E 22  > 0 , J 22 − R 22 , and J 31 are nonsingular . Theorem 5. Consider an index- 2 pHD AE sys tem of the form (27). Construct an interpolatory model reduction basis V =  V 𝑇 1 V 𝑇 2 V 𝑇 3  𝑇 =  ( 𝜎 1 E − A ) − 1 B 𝒃 1 , · · · ( 𝜎 𝑟 E − A ) − 1 B 𝒃 𝑟  ∈ I C 𝑛 × 𝑟 , (28) partitioned as the sys tem. Define b V : = diag ( I , V 2 , I ) . Then the r educed sys tem b E ¤ x 𝑟 = ( b J − b R ) x 𝑟 + b B u ( 𝑡 ) , y 𝑟 = b C x 𝑟 + b D , (29) with b E = b V 𝑇 E b V , b J = b V 𝑇 J b V , b R = b V 𝑇 R b V , b B = b V 𝑇 B = b V 𝑇 B − b V 𝑇 P , b C = C b V = B 𝑇 b V 𝑇 + P 𝑇 b V 𝑇 , and b D = D , (30) is still pH , matc hes the polynomial par t of the tr ansf er function, and satisfies the tang ential interpolation conditions, i.e., H ( 𝜎 𝑖 ) 𝒃 𝑖 = b H ( 𝜎 𝑖 ) 𝒃 𝑖 , f or 𝑖 = 1 , 2 , . . . , 𝑟 . Proof. It f ollo ws from the definitions of V and b V that span ( V ) ⊆ span ( b V ) . There- f ore, the resulting reduced system automatically satisfies the inter polation conditions. Since b V does not alter the alg ebraic constraints, the pol ynomial par t of its transfer function is still D , matching that of the original model. The reduced system in (30) is a pHD AE as this one-sided projection retains the original pH structure. 16 C. A. Beattie, S. Gugercin, and V . Mehr mann 5 Numerical Eperiments The preceding anal ysis presumed that inter polation points and tang ent directions w ere specified beforehand. W e consider no w how one might make choices that generall y produce effective approximations with respect to the H 2 sys tem measure. H 2 -inspired structure-preserving interpolation: between the full model H ( 𝑠 ) and the reduced model H 𝑟 ( 𝑠 ) :   H − H 𝑟   H 2 =  1 2 𝜋 ∫ ∞ −∞   H ( 𝚤 𝜔 ) − H 𝑟 ( 𝚤 𝜔 )   2 𝐹 𝑑 𝜔  1 / 2 , where 𝚤 2 = − 1 and k M k 𝐹 denote the Frobenius norm of a matrix M . T o ha v e a finite H 2 error nor m the polynomial parts of H 𝑟 ( 𝑠 ) and H ( 𝑠 ) need to match. T o make this precise, wr ite H ( 𝑠 ) = G ( 𝑠 ) + P ( 𝑠 ) and H 𝑟 ( 𝑠 ) = G 𝑟 ( 𝑠 ) + P 𝑟 ( 𝑠 ) as where G ( 𝑠 ) and G 𝑟 ( 𝑠 ) are strictly proper transf er functions, and P ( 𝑠 ) and P 𝑟 ( 𝑠 ) are the polynomial par ts. Theref ore, if H 𝑟 ( 𝑠 ) is the H 2 -optimal approximation to H ( 𝑠 ) , then P 𝑟 ( 𝑠 ) = P ( 𝑠 ) and G 𝑟 ( 𝑠 ) is the H 2 -optimal appro ximation to G ( 𝑠 ) . This sugg ests that one decomposes H ( 𝑠 ) into its rational and polynomial par ts H ( 𝑠 ) = G ( 𝑠 ) + P ( 𝑠 ) and then applies H 2 optimal reduction to G ( 𝑠 ) . How ev er , this requires the e xplicit construction of G ( 𝑠 ) . This problem was resolv ed in [17] for unstructured inde x- 1 and index- 2 D AE s. On the other hand, for the ODE case, [16] proposed a pH structure preserving algor ithm f or minimizing the H 2 norm. In this section, we aim to unify these tw o approaches. Firs t, w e br iefly revisit the interpolator y H 2 optimality conditions; for details w e ref er the reader to [1, 5, 14] and the references therein. Since H 2 optimality f or the D AE case boils do wn to optimality f or the ODE par t, w e f ocus on the latter . Let G 𝑟 ( 𝑠 ) = Í 𝑟 𝑖 = 1 c 𝑖 b 𝑇 𝑖 𝑠 − 𝜆 𝑖 be the pole-residue decomposition of G 𝑟 ( 𝑠 ) . For simplicity w e assume simple poles. If G 𝑟 ( 𝑠 ) is an H 2 optimal approximation to G ( 𝑠 ) , then G ( − 𝜆 𝑖 ) b 𝑖 = G 𝑟 ( − 𝜆 𝑖 ) b 𝑖 , c 𝑇 𝑖 G ( − 𝜆 𝑖 ) = c 𝑇 𝑖 G 𝑟 ( − 𝜆 𝑖 ) , and c 𝑇 𝑖 G 0 ( − 𝜆 𝑖 ) b 𝑖 = c 𝑇 𝑖 G 0 𝑟 ( − 𝜆 𝑖 ) b 𝑖 . for 𝑖 = 1 , 2 , . . . , 𝑟 where 0 denotes the derivate with respect to 𝑠 . Theref ore an H 2 optimal reduced model is a bitang ential Hermite inter polant where the inter polation points are the mir ror images of its poles and the tangent directions are the residue directions. Since the optimality conditions depend on the reduced model to be computed, this requires an iterative algor ithm, which the Iterativ e Rational Krylo v Algor ithm [14]. For other approaches in H 2 optimal appro ximation, see, e.g., [1, 5]. Follo wing [16], in order to preserv e the structure, we will satisfy only a subset of these conditions. W e will make sure that the inter polation points will be the mirror images of the reduced order poles and enf orce either left or r ight-tang ential interpolation conditions without the derivate conditions. Ho w ev er, intuitivel y , one might e xpect that in the pH setting this may not cause too much de viation from true optimality , since this a Galerkin projection, i.e., Z = V . Consider the semi-explicit inde x- 2 case. Starting with an initial selection of interpolation points { 𝜎 1 , 𝜎 2 , . . . , 𝜎 𝑟 } and the tangent directions { 𝒃 1 , 𝒃 2 , . . . , 𝒃 𝑟 } , construct the inter polator y reduced pHDAE b H ( 𝑠 ) as in Theorem 3. Let b H ( 𝑠 ) = Structure-preserving Inter polatory Model Reduction for pHDAE s 17 Í 𝑟 𝑖 = 1 b ℓ 𝑖 b 𝒃 𝑇 𝑖 𝑠 − 𝜆 𝑖 + D be the pole-residue decomposition. Since initially the optimality condition 𝜎 𝑖 = − 𝜆 𝑖 ( b J − b R , b E ) is not (generall y) satisfied, choose − 𝜆 𝑖 ( b J − b R , b E ) as the next set of inter polation points and b 𝒃 𝑖 as the next set of tang ent directions. This is repeated until conv er gence upon which the reduced model H 𝑟 ( 𝑠 ) is not only a structure-preserving pHDAE , but also satisfies 𝜎 𝑖 = − 𝜆 𝑖 ( b J − b R , b E ) . Numerical Example: W e illustrate the discussed procedure with the incompress- ible fluid flo w model of the Oseen equations, from [18, §4.1] 𝜕 𝑡 𝑣 = − ( 𝑎 · ∇) 𝑣 + 𝜇 Δ 𝑣 − ∇ 𝑝 + 𝑓 in Ω × ( 0 , 𝑇 ] , 𝑣 = 0 , on 𝜕 Ω × ( 0 , 𝑇 ] , 0 = − div 𝑣 , in Ω × ( 0 , 𝑇 ] , 𝑣 = 𝑣 0 , in Ω × 0 , where 𝑣 and 𝑝 are the v elocity and pressure variables, 𝜇 > 0 is the viscosity , and Ω = ( 0 , 1 ) 2 with boundar y 𝜕 Ω . 𝑓 is an externally imposed body f orce that f or simplicity is assumed to be separable: 𝑓 ( 𝑥 , 𝑡 ) = 𝑏 ( 𝑥 ) 𝑢 ( 𝑡 ) . A finite-difference discretization on a stagg ered rectangular gr id leads to a single-input/single-output inde x- 2 pHD AE of the form (16), see [18]. In our model, we used a uniform grid with 50 gr id points yielding a descriptor sys tem pHD AE of order 𝑛 = 7399 , of which 𝑛 1 = 4900 degrees of freedom are for velocity and 𝑛 2 = 2499 f or pressure. W e initialized our approach to reduce the order of the tranf er function to 𝑟 = 1 , 2 , . . . , 10 with logarithmically spaced inter polation points in the inter val [ 10 − 2 , 10 4 ] . For ev er y 𝑟 v alue, w e compute k H − b H k ∞ / k H k ∞ where k H k ∞ = sup 𝜔 ∈ I R | H ( 𝚤 𝜔 ) | . The results in Figure 1 sho w that the reduced transf er function accuratel y approximates the or iginal one; f or the reduction to order 𝑟 = 10 , the relative er ror is around 10 − 5 , illustrating the success of the proposed framew ork. 1 2 3 4 5 6 7 8 9 10 r 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 || H - H r || / || H || Relative H error vs reduced order Fig. 1 Evolution of the model reduction error for Oseen e xample as 𝑟 varies 18 C. A. Beattie, S. Gugercin, and V . Mehr mann 6 Conclusion W e ha v e presented inter polator y model reduction methods f or se v eral classes of larg e scale linear time-inv ar iant por t-Hamiltonian differential-algebraic sy stems. W e ha v e sho wn ho w constraints are represented in the transfer function so that the polynomial part can be preserved with inter polatory methods. W e ha v e illustrated the results with numerical example from flo w control. A ckno wledgment The w orks of Beattie and Gug ercin w ere supported in par ts b y NSF through Grant DMS-1819110. The w ork of Mehrmann was suppor ted b y Deutsche Forschungsg emeinschaft through CR C 910 within project A02. Parts of this work were completed while Beattie and Gugercin visited TU Berlin, f or which Gugercin ackno wledg es suppor t through CR C 910 and Beattie ackno w ledg es support through CR C TRR 154 as well as DFG R esearch T raining Group D AED ALUS. Ref erences 1. A.C. Antoulas, , C.A. Beattie, and S. Gugercin. Interpolatory Methods for Model Reduction . Society for Industrial and Applied Mathematics, to appear , 2020. 2. C. Beattie, V . Mehr mann, H. Xu, and H. Zwart. Linear por t-Hamiltonian descr iptor systems. Mathematics of Control, Signals, and Syst ems , 30(4):1–27, 2018. 3. C.A. Beattie and S. Gugercin. Interpolator y projection methods for str ucture-preserving model reduction. Syst ems & Control Lett ers , 58(3):225–232, 2009. 4. Christopher Beattie and Serkan Gugercin. Structure-preser ving model reduction f or nonlinear port-Hamiltonian sys tems. In 2011 50th IEEE Confer ence on Decision and Contr ol and European Control Confer ence , pages 6564–6569. IEEE, 2011. 5. Christopher Beattie and Serkan Gug ercin. Chapter 7: Model Reduction by Rational Interpo- lation , pages 297–334. 2017. 6. A. Binder , V . Mehrmann, A. Miedlar, and P . Schulze. A Matlab toolbox for the regularization of descr iptor sys tems ar ising from generalized realization procedures. Preprint 24–2015, Institut für Mathematik, TU Berlin, 2015. 7. A. Bunse-Gerstner, R. By ers, V . Mehrmann, and N. K. Nichols. F eedback design for regular- izing descr iptor sys tems. Linear Algebr a Appl. , 299(1–3):119–151, 1999. 8. S. L. Campbell, P . Kunkel, and V . Mehr mann. Regularization of linear and nonlinear descriptor sys tems. In L. T . Biegler, S. L. Campbell, and V . Mehr mann, editors, Contr ol and Op timization with Differential-Alg ebraic Constr aints , number 23 in Adv ances in Design and Control, pages 17–36, 2012. 9. S. Chaturantabut, C. Beattie, and S. Gugercin. Structure-preserving model reduction f or nonlinear port-Hamiltonian sy stems. SIAM Journal on Scientific Computing , 38(5):B837– B865, 2016. 10. L. Dai. Singular Contr ol Syst ems , volume 118 of Lecture No tes in Control and Information Sciences . Spr inger - V erlag, Berlin, 1989. 11. H. Egger , T . Kugler , B. Liljegren-Sailer, N. Marheinek e, and V . Mehrmann. On s tructure- preserving model reduction f or damped w av e propag ation in transpor t networks. SIAM Jour nal on Scientific Computing , 40(1):A331–A365, 2018. 12. K. Gallivan, A. V andendor pe, and P . V an Dooren. Model reduction of MIMO sys tems via tangential inter polation. SIAM Jour nal on Matrix Analysis and Applications , 26(2):328–349, 2005. 13. F . R. Gantmacher . The Theor y of Matrices , volume II. Chelsea Publishing Company , Ne w Y ork, N.Y ., 1959. Structure-preserving Inter polatory Model Reduction for pHDAE s 19 14. S. Gugercin, A.C. Antoulas, and C.A. Beattie. H 2 Model Reduction f or Larg e-Scale Linear Dynamical Systems. SIAM J. Matrix Anal. Appl. , 30(2):609–638, 2008. 15. S. Gugercin, R. V . Pol yug a, C. Beattie, and A. v an der Schaft. Interpolation-based H2 model reduction f or por t-Hamiltonian systems. In Proceedings 48th IEEE Confer ence on Decision and Control , pages 5362–5369. IEEE, 2009. 16. S. Gug ercin, R.V . P olyug a, C.A. Beattie, and A.J. van der Sc haft. Structure-preserving tangen- tial inter polation for model reduction of por t-Hamiltonian systems. Automatica , 48:1963–1974, 2012. 17. S. Gugercin, T . Stykel, and S. W y att. Model reduction of descriptor sys tems by inter polator y projection methods. SIAM Jour nal on Scientific Computing , 35(5):B1010–B1033, 2013. 18. S.-A. Hauschild, N. Marheineke, and V . Mehrmann. Model reduction tech- niques f or linear constant coefficient port-hamiltonian differential-alg ebraic sys tems. Preprint 02-2019, Ins titut für Mathematik, TU Berlin, D-10623 Berlin, FR G, 2019. url: http://www.math.tu-berlin.de/preprints/, https://arxiv.org/abs/1901.10242 . 19. M. Heinkensc hloss, D. C. Sorensen, and K. Sun. Balanced truncation model reduction for a class of descr iptor sy stems with application to the Oseen equations. SIAM Jour nal on Scientific Computing , 30(2):1038–1063, 2008. 20. P . Kunk el and V . Mehrmann. Differ ential-Alg ebr aic Equations . European Mathematical Society (EMS), Zür ich, 2006. 21. A.J. Ma y o and A.C. Antoulas. A framew ork for the solution of the g eneralized realization problem. Linear Alg ebra Appl. , 425(2-3):634–662, 2007. 22. C. Mehl, V . Mehr mann, and M. W ojtylak. Distance problems for dissipative hamiltonian and port-hamiltonian descr iptor sys tems. T echnical repor t. In preparation. 23. C. Mehl, V . Mehrmann, and M. W ojtylak. Linear alg ebra proper ties of dissipative por t- Hamiltonian descr iptor sy stems. SIAM Journal on Matrix Analysis and Applications , 39(3):1489–1519, 2018. 24. V . Mehr mann and R. Morandin. Structure-preser ving discretization for por t-hamiltonian descriptor sys tems. Preprint 05–2019, Inst. f. Mathematik, TU Berlin, 2019. 25. V . Mehrmann and T . Stykel. Balanced truncation model reduction for larg e-scale sys tems in descriptor form. In Dimension Reduction of Larg e-Scale Syst ems , pages 83–115. Springer , Berlin, 2005. 26. R. V . P ol yuga and A. van der Schaft. Structure-preserving moment matching for por t- Hamiltonian sys tems: Arnoldi and Lanczos. IEEE T ransactions on Automatic Contr ol , 56(6):1458–1462, 2011. 27. R. V . P olyug a and A. van der Schaft. Effor t- and flow -constraint reduction methods for structure-preserving model reduction of por t-Hamiltonian systems. Systems and Control Let- ter s , 61(3):412–421, 2012. 28. R.V . Polyug a and A.J. van der Schaft. Structure preserving model reduction of por t- Hamiltonian systems by moment matching at infinity . Automatica , 46:665–672, 2010. 29. T . Styk el. Balanced truncation model reduction f or semidiscretized Stok es equation. Linear Alg ebr a and its Applications , 415(2-3):262–289, 2006. 30. A. van der Schaft. Port-Hamiltonian differential-algebraic systems. In Sur vey s in Differential- Alg ebr aic Equations. I , pages 173–226. Spr inger , Heidelberg, 2013. 31. T . W olf, B. Lohmann, R. Eid, and P . Koty czka. Passivity and structure preser ving order reduction of linear port-Hamiltonian sy stems using Kr ylo v subspaces. European Journal of Control , 16(4):401–406, 2010.