[basicstyle=]CpH represents a circuit by a [basicstyle=]pHedge object and its Kron matrix \(F\)  plus tree \(\texttt {T}\)  and cotree \({\texttt {T}^*}\) .
The constructor makes a [basicstyle=]CpH object [basicstyle=]P from a [basicstyle=]pHcircuit object [basicstyle=]P0 and optionally a proposed tree T.
If T is not given, this constructor computes one. If a normal tree \(\texttt {T}\)  and cotree \({\texttt {T}^*}\)  have been found, the Kron matrix \(F\)  is computed as in Subsection REF .
[basicstyle=]P stores sufficient data to reverse the construction, i.e., recover [basicstyle=]P0 exactly from [basicstyle=]P.
There are two methods to construct from a CpH object [basicstyle=]P a function \(f(t,x,\dot{x})\)  that implements CpH Model 2 and is accepted by [basicstyle=]ode15i [1]}: one that “interprets” [basicstyle=]P as this \(f\) , and one that for more efficiency generates a Matlab [basicstyle=]daefcn function for this \(f\) .
We actually solve the DAE by a slight adaptation of [basicstyle=]ode15i which at each \(t_i\)  returns numerical approximations of \(\dot{x}_i\)  as well as of \(x_i\) .
This supports computing the control-output \(y\) , which is a function of \(t\) , \(x\)  and \(\dot{x}\) .
Computation of a normal tree and its Kron matrix
From Theorem REF , a normal tree and its Kron matrix \(F\)  are a byproduct of computing a suitable RREF of the incidence matrix \(A\) .
Matlab provides the function [basicstyle=][R,T] = rref(A) to do this using Gauss-Jordan elimination with partial pivoting.
The resulting [basicstyle=]R is the echelon form; [basicstyle=]T, the integer vector of pivot column indices, is the tree; and \(F\)  can be formed from the non-pivot columns, see Corollary REF .
A faster version of Matlab's RREF, [basicstyle=]frref, that uses sparse QR algorithm to compute RREF for sparse matrices has been provided by [2]}.
Since floating point is used in the computations, the obtained \(F\)  can differ from integer values at roundoff level, so we round it to the nearest integer.
Both [basicstyle=]rref and [basicstyle=]frref are very inefficient for finding a tree.
A much faster approach is to use Kruskal's algorithm (KA) [3]} for finding a minimum spanning tree; see below and Figure REF . In this algorithm, we do not sort the edges but provide them in the order \({\hspace{0.2pt}{\hspace{0.2pt}\texttt {V}}}, {\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}, {\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}, {\hspace{0.2pt}{\hspace{0.2pt}\texttt {L}}}, {\hspace{0.2pt}{\hspace{0.2pt}\texttt {I}}}\) .
The ordering in each of these sets does not matter.
Starting with \(\texttt {T}= \emptyset \) , KA adds an edge \(e\)  to \(\texttt {T}\)  if \(\texttt {T}\cup \lbrace e\rbrace \)  is acyclic.
This results in a normal tree as long as (A1) and (A2) hold.
If a normal tree \(\texttt {T}\)  and cotree \({\texttt {T}^*}\)  are found, the Kron matrix \(F\)  is then computed from (REF ) by
[basicstyle=]F = -inv(A(1:n-1,T))*A(1:n-1,coT),
using Matlab notation.
Surprisingly, using the [basicstyle=]inv function is more efficient than finding \(F\)  through Matlab's backslash operator (see Figure REF ).
<FIGURE>The timing results are produced with Matlab version R2021a on Mac OSX with 2.2 GHz Quad-Core Intel Core i7
and 16 GB 1600 MHz DDR3.
group-minimum-digits = 3,group-separator = ,
We use sparse incidence matrices obtained from random connected graphs with \(n\)  vertices and \(b = 3n\)  edges.
In Figure REF , we plot the CPU time
of [basicstyle=]rref, [basicstyle=]frref, and KA for
\(n = 20, 40, \ldots , 1280\) . For each set of data points
(number of edges, CPU time), we also plot \(c_1b^c\) , where the constants \(c_1, c\)  are determined from a least-squares fit of these data.
In Figure REF , we plot the CPU time for KA
and computing \(F\)  using [basicstyle=]inv and \,
where \(n = 2000, 4000, \ldots , 128000\) , and also plot the corresponding \(c_1b^c\)  fits.
The time for the overall computation of \(F\)  is dominated by [basicstyle=]inv.
This matrix can be directly derived from the incidence matrix of the tree by an algorithm that is purely based on logical operations [4]}. We are planning to compare the efficiency of this approach to our current method.

Consistent initial conditions
Computing consistent initial conditions (CICs),
as in Step 6 of the -method in §, is challenging for general DAEs. Usually, this involves index reduction and finding hidden constraints [5]}, [6]}.
A solver such as [basicstyle=]ode15i requires consistent values for derivatives as well. Matlab provides the function [basicstyle=]decic that seeks CICs for an index-1 DAE prior to solving by [basicstyle=]ode15i. The [basicstyle=]decic function works well for a semi-explicit DAE but a CpH DAE is usually not so, and [basicstyle=]decic may give wrong \(\dot{x}\)  values.
An example is \(x_1+x_2=0,\;x_1-\dot{x}_2=0\) , which comes from a circuit comprising just a resistor \(R_1=1\)  and capacitor \(C_2=1\) .
In our CpH approach, having correct approximations for derivatives is crucial, as we insert them in () to obtain the output \(y\) .
The DAE being SA-amenable index-1, as here, eases the CICs task.
From results of §, for both models 1 and 2, one can fix \(n_{\hspace{0.2pt}\texttt {c}}\)  charges \(q_j\)  and \(n_{\hspace{0.2pt}\texttt {L}}\)  fluxes \(\varphi _j\) .
In case of independent elements, they can be exactly the tree capacitor charges and cotree inductor fluxes. In general, a correct choice comes from rank properties of the \(f_{\hspace{0.2pt}\texttt {C}}\)  and \(f_{\hspace{0.2pt}\texttt {l}}\)  rows of \(\mathbf {J}\) , see (REF ).
This ensures a unique solution—locally in the nonlinear case, globally in the linear.
We use the \(\Sigma \) -method's SSS, cf. (REF ). At stage \(k=-1\) , for a given initial guess \(\widetilde{x}_0\)  at \(t_0\)  solve \(0=f_{\hspace{0.2pt}\texttt {C}},f_{\hspace{0.2pt}\texttt {l}},f_{\hspace{0.2pt}\texttt {d}},f_{\hspace{0.2pt}\texttt {D}},r\)  in (REF ) for the whole vector \(x\) : a full-rank, generally underdetermined system.
Let \(x_0\)  be the result.
Now let \(x_1(h)\)  denote the result of one step of implicit Euler from \((t_0,x_0)\)  to \(t_1=t_0+h\) .
For a DAE with analytic coefficients, this has a series expansion \(x_1(h) = x_0 + \dot{x}(t_0)h + c_1h^2 + c_2h^3 + \cdots \)  where \(x(t)\)  is the exact solution through \((t_0,x_0)\) .
Hence the central difference \(\frac{x_1(h)-x_1(-h)}{2h}\)  is an \(O(h^2)\)  accurate estimate of \(\dot{x}(t_0)\) .
Error estimates, and higher-order approximations, can be found by Richardson extrapolation.
We aim to detail this in a future paper.

Numerical Examples
First we consider the circuit from Example REF  with parameter values
\(C_1=C_2=5\cdot 10^{-6}\) , \(G=R=1\) , \(L_1=L_2=0.1\)  and source functions \(V(t)=10t\sin (200\pi t)\) , \(I(t)=10\sin (10t)\) .
We integrated
the Model 2 formulation over the time interval \([0,0.2]\)  using the initial guess \(\widetilde{x}_0=[1,\dots ,1]^.The numerical solution for \) x=[qC2,qC1,L1,L2,R,G]\( is depicted in Figure~\ref {fg:outputRUNEX} (a). The output \) y=[V,I]\( is given in Figure~\ref {fg:outputRUNEX} (b).The analytical solution can be found for this problem, and the numerical one agrees with it (see red dashed line in Figure~\ref {fg:outputRUNEX} (b)).\)
<FIGURE>The second example is the Diode clipper circuit in Figure REF  (a) from [7]} (it can be found in audio-distortion devices).
<FIGURE>It contains one linear resistor (\(R=1000\) ) and two diodes that are treated as nonlinear conductors (voltage-controlled dissipators) modeled by Shockley's equation, i.e., the diode current is given as a function of its branch voltage \(\iota (\nu )=I_s\cdot (e^{\nu /\nu _T}-1)\) , with typical parameters \(I_s= 10^{-13}\)  (the saturation current), and \(\nu _T= 0.025\)  (the thermal voltage).Standard notations for this diode model: the \(T\)  in \(\nu _T\)  is about temperature, not a tree.
The circuit contains a voltage source with \(V(t)=(\widehat{V}\,t/t_{max})\,\sin (\omega t)\) , with \(\omega = 2\pi f\) , a linearly increasing sinusoidal excitation with frequency \(f={1}{}\)  and maximum amplitude \(\widehat{V}=2\)  at \(t_{\max }=0.03\) , and an artificial current source \(I(t)=0\)  that acts as voltmeter to measure the output voltage.
The Model 2 formulation is integrated over the interval \([0, t_{\max }]\)  using an initial guess \(\widetilde{x}_0=[1,\dots ,1]^.The control-output \) y=[V,I]\( is plotted in Figure~\ref {fg:DiodeClipper} (b).\)Matlab code for the two examples is given in Appendix .
To show larger circuits being handled, we created an example where a graph is generated and populated randomly with circuit elements.
If it satisfies (A1), (A2) (checked during construction) its CpH Model 2 form is constructed and solved with random initial values.
<FIGURE>An example with \(b=15\)  branches and \(n=9\)  nodes is in Figure REF  (a); Figure REF  (b) displays the \(F\)  matrix, annotated with the specifications of the elements.
The \(x\)  solution components and the \(y\)  output components
are plotted in Figure REF .
<FIGURE>
Conclusions
Old and new.
The port-Hamiltonian concept is a powerful way to bring theory nearer physics, and suggest new computational approaches.
SA-amenability is a powerful concept in another way, making possible, or cheap to implement, numerical methods that without it are impossible or expensive.
Our contribution is to show that a particular mix of—mostly well known—ingredients gives a DAE that is both port-Hamiltonian and SA-amenable.
That it has index \(\le 1\)  was an unexpected bonus, which we suspect comes from obeying the pH injunction that the Hamiltonian be an algebraic function of the state variables.
Writing Kirchhoff's laws in terms of the \(F\)  matrix goes back at least to [8]}, [9]}, where also the zero blocks of \(F\)  are exploited though less transparently than here.
The existence of a normal tree is shown in [9]}.
Our main proof is a distillation of that in [11]} that the BOM circuit model's DAE is SA-amenable.
However, the BOM DAE is not pH and can have index 2.
Applying pH theory to electrical circuits is a topic of current research, see e.g. [12]}, [13]}.
In [12]}, selection of a ground node—not needed by CpH—plays a crucial role in defining the Dirac structure; and there is no discussion of the DAE index.
The aim of [13]} is to put the MNA equations in port-Hamiltonian form; this can lead to a DAE of index 2.
Neither work discusses SA-amenability.
To the best of our knowledge the CpH method construction, and SA-amenability and related proofs, are new.
So also appears to be our proof of Theorem REF , and the resulting direct counts of the number of linearly independent \({\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}{\hspace{0.2pt}{\hspace{0.2pt}\texttt {V}}}\) -cycles (including \({\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}\) -cycles) \(n_{\hspace{0.2pt}\texttt {C}}=n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}-r_{{\hspace{0.2pt}{\hspace{0.2pt}\texttt {V}}}{\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}} + r_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {V}}}\) , and similarly of \({\hspace{0.2pt}{\hspace{0.2pt}\texttt {L}}}{\hspace{0.2pt}{\hspace{0.2pt}\texttt {I}}}\) -cutsets \(n_{\hspace{0.2pt}\texttt {l}}=r_{{\hspace{0.2pt}{\hspace{0.2pt}\texttt {V}}}{\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}{\hspace{0.2pt}{\hspace{0.2pt}\texttt {L}}}} - r_{{\hspace{0.2pt}{\hspace{0.2pt}\texttt {V}}}{\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}}\) .
Our CICs method in §REF  is new, but related to one in [16]}.
Our method is easy to automate: as said in §REF , from a [basicstyle=]CpH object representing a circuit we generate code simulating that circuit, in the form of a [basicstyle=]daefcn function accepted by an index-1 solver such as [basicstyle=]ode15i.
Other advantages of being SA-amenable is the systematic reduction to an ODE Subsection REF , the easy eigenvalue analyis for LTI systems Subsection REF , or the use of Mattsson-Söderlind dummy derivatives [17]}.

Future work.
A key feature of port-Hamiltonian theory is that a “Dirac join” of several pH systems—joining their Dirac structures in a way that makes a Dirac structure on the union of the systems—is again a pH system, whose Hamiltonian is the sum of the Hamiltonians of the parts.
Our Matlab implementation does not yet have this property and it is a priority to provide it.
This would follow in the steps of the port-Hamiltonian Modelica system recently reported in [18]}.
At a more modest level, our implementation currently only supports independent, possibly nonlinear, circuit elements.
We are working on upgrading it to handle the full generality of models in this paper, i.e. circuits containing dependent elements forming multi-terminal components (e.g. transistors) and other more complex circuit configurations (e.g. controlled sources).
In a future paper we aim to give proofs and algorithm details for the method in §REF .

Acknowledgments
We thank Prof Volker Mehrmann and Prof Caren Tischendorf for valuable advice.
We thank Jacek Kierzenka of the MathWorks for writing us an adapted version of the solver [basicstyle=]ode15i that supports writing a DAE in control input-output form (REF ).

Terminology of circuit graphs vs other kinds of graph
We summarise where it is necessary or convenient for circuit graph terms to differ from what is “often” the case, i.e. common in other graph theory areas.
Let \(G=(V,{\hspace{0.2pt}\texttt {E}},\alpha )\)  be a circuit graph.
Riaza [19]} reserves the terms branch, node for elements of a circuit, and edge, vertex for the corresponding elements of the circuit graph. We keep this distinction.

A circuit graph allows several edges between a given pair of vertices and forbids edges from a vertex to itself.
Often the opposite of one or both of these is true.

Often the induced subgraph of a subset of edges discards any vertex that is not the end of some edge in the subset.
Here it is more useful to keep all vertices by default.

Often a spanning subgraph of \(G=(V,{\hspace{0.2pt}\texttt {E}},\alpha )\)  simply means one whose vertex set is all of \(V\) .
Here, being connected is essential.

A link is often called a chord; but both link and chord have multiple meanings in graph literature.

Giving \(G\)  an orientation is just a housekeeping aid for using Kirchoff's laws.
It doesn't change the topology and is thus quite different from making \(G\)  a directed graph.
The circuit physics is the same whatever orientation is chosen.
The Signature Method (\(\Sigma \) -method) for DAEs
Our preferred SA method for DAEs is the Signature method (\(\Sigma \) -method) [20]}.
The \(\Sigma \) -method can be applied to regular nonlinear DAEs of arbitrary high order \(k\)  of the form
\(f(t,x,\dot{x},\dots ,x^{(k)})=0,\) 

with \(f:\mathbb {T}\times \mathbb {R}^n\times \dots \times \mathbb {R}^n\rightarrow \mathbb {R}^n\)  sufficiently smooth.
The \(\Sigma \) -method consists of the following steps:
Building the signature matrix \(\Sigma =[\sigma _{ij}]_{i,j=1,\dots ,n}\)  with
\(\sigma _{ij}={\left\lbrace \begin{array}{ll}\text{highest order of derivative of $x_j$ in $f_i$,} \\ -\infty \text{ if $x_j$ does not occur in $f_i$.}\end{array}\right.}\) 

Finding a highest value transversal (HVT) of \(\Sigma \) , i.e.,
a set
\({\cal T}=\lbrace (1,j_1),(2,j_2),\dots ,(n,j_n)\rbrace \) ,
where \((j_1,\dots ,j_n)\)  is a permutation of \((1,\dots ,n)\) ,
with maximal value \(\text{Val}({\cal T}):=\sum _{(i,j)\in {\cal T}}\sigma _{ij}.\)

Computing offsets \(c=[c_i]_{i=1,\dots ,n}\)  and \(d=[d_j]_{j=1,\dots ,n}\)  with \(c_i \ge 0\) , \(d_j \ge 0\)  such that
\(\text{$d_j-c_i\ge \sigma _{ij}$ for all $i,j=1,\dots ,n$, with equality for $(i,j)\in \mathcal {T}$}.\) 


Forming the system Jacobian \(\mathbf {J} = \begin{bmatrix}\mathbf {J}_{ij}\end{bmatrix}_{i,j=1}^n\)  where (the two formulae are equivalent)
\({\mathbf {J}}_{ij}:={\left\lbrace \begin{array}{ll}\displaystyle \partial f_i/\partial x_j^{(\sigma _{ij})} &\text{ if } d_j-c_i=\sigma _{ij}, \\0 &\text{ otherwise}\end{array}\right.}\quad = \quad {\left\lbrace \begin{array}{ll}\displaystyle \partial f_i/\partial x_j^{(d_j-c_i)} &\text{if $d_j-c_i\ge 0$}, \\0 &\text{otherwise}\end{array}\right.}\) 


Form the reduced derivative array \({\cal F}(t,{\cal X})=0\)  where
\({\cal X}=[x_1,\dot{x}_1,\dots ,x_1^{(d_1)}, \dots , x_n,\dot{x}_n,\dots ,x_n^{(d_n)}]^\ {\cal F}=[f_1,\dot{f}_1,\dots ,f_1^{(c_1)}, \dots , f_n,\dot{f}_n,\dots ,f_n^{(c_n)}]^\) 


Success check: if \({\cal F}(t,{\cal X})=0\) , considered locally as an algebraic system of \(M=n+\sum _{i=1}^n c_i\)  equations in \(N=n+\sum _{j=1}^n d_j\)  unknowns, has a solution \((t^*,{\cal X}^*)\in \mathbb {T}\times \mathbb {R}^N\)  where \(\mathbf {J}\)  is nonsingular, then \((t^*,{\cal X}^*)\)  is a consistent point, the method “succeeds”, and the DAE is SA-amenable.
In case of success the structural index
\(\nu _s:=\max _i c_i + (\text{$1$ if some $d_j=0$, $0$ otherwise})\) 

is \(\ge \)  the differentiation index \(\nu _d\) .
Also \(\text{Val}(\Sigma )\) , defined as the value of the highest value transversal \(\mathcal {T}\)  and also equal to \(\sum _j d_j - \sum _i c_i\) , equals the number of degrees of freedom of the system.
We have the standard solution scheme (SSS): starting with \(k = -\max d_j\) , in order of increasing stage number \(k\) ,
\(\begin{aligned}& \text{solve the equations}& &f_i^{(k+c_i)}=0&&\text{for those $i$ such that $k+c_i\ge 0$} \\& \text{for the unknowns}&&x_j^{(k+d_j)}&&\text{for those $j$ such that $k+d_j\ge 0$}\end{aligned}\biggr \rbrace \;.\) 

One interpretation of this is as a numerical scheme that first finds a consistent point and then expands the DAE solution through this point as a Taylor series; then, \(k\)  goes up to some Taylor order \(k=k_{\max }>0\) .
However the meaning used in this paper is that we go up to \(k=0\) : stages up to \(k=-1\)  are “numerical”, and a consistent point is a set of values that solve these equations; while \(k=0\)  is “symbolic”—a set of equations by which (using the Implicit Function Theorem) the highest derivatives \(x_j^{(d_j)}\)  are proved to be functions of lower derivatives, thus defining an implicit ODE.
Let \(\mathbf {J}_k\)  be the Jacobian of the equations in (REF ) at stage \(k\) .
Theorem B.1 
([20]})
For an SA-amenable DAE, each \(\mathbf {J}_k\)  is of full row rank, and for all \(k\ge 0\)  is equal to the system Jacobian \(\mathbf {J}\)  (hence is nonsingular).

Port-Hamiltonian dynamics
Port-Hamiltonian system and Dirac structure
A physical system is modelled by flows of energy within it, these flows being represented by ports: (flow, effort) pairs of variables whose product has the dimension of power.
The following paradigm is taken from [22]}.
The \(\mathpzc {F}_\bullet \)  are linear spaces and \(\mathpzc {E}_\bullet =\mathpzc {F}_\bullet ^*\)  their duals.
In our context they are some \(\mathbb {R}^d\) , so the duality inner product \(\langle \,\mathpzc {f}\hspace{-1.0pt},\,\mathpzc {e}\hspace{0.5pt}\,\rangle =\mathpzc {f}\hspace{-1.0pt}\mathop {\cdot }\mathpzc {e}\hspace{0.5pt}\) .
In general \(\mathcal {Z}\)  is a manifold, with \(T_z\mathcal {Z}\)  its tangent space at the state \(z\in \mathcal {Z}\) , and \(T_z^*\mathcal {Z}\)  the corresponding co-tangent space.
In our context \(\mathcal {Z}\)  is some \(\mathbb {R}^d\) , so \(T_z\mathcal {Z}\)  and \(T_z^*\mathcal {Z}\)  can also be regarded as \(\mathbb {R}^d\) , again with duality inner product.
A port-Hamiltonian system gives mathematical form to Figure REF . It is defined by
[(pH1)]
A state space \(\mathcal {Z}\)  and a Hamiltonian \(\mathcal {H}: \mathcal {Z}\rightarrow \mathbb {R}\) , defining energy-storage, corresponding to an energy-storing port \((\mathpzc {f}\hspace{-1.0pt}_\mathcal {S},\mathpzc {e}\hspace{0.5pt}_\mathcal {S})\in T_z\mathcal {Z}\times T_z^*\mathcal {Z}\) .

A dissipative structure
\(\mathcal {R}&\subset \mathpzc {F}_\mathcal {R}\times \mathpzc {E}_\mathcal {R}&\text{having $\langle \,\mathpzc {f}\hspace{-1.0pt},\,\mathpzc {e}\hspace{0.5pt}\,\rangle \le 0$ for all $(\mathpzc {f}\hspace{-1.0pt},\mathpzc {e}\hspace{0.5pt})\in \mathcal {R}$},\) 

corresponding to an energy-dissipating port \((\mathpzc {f}\hspace{-1.0pt}_\mathcal {R}, \mathpzc {e}\hspace{0.5pt}_\mathcal {R}) \in \mathpzc {F}_\mathcal {R}\times \mathpzc {E}_\mathcal {R}\) .

An external port \((\mathpzc {f}\hspace{-1.0pt}_\mathcal {P}, \mathpzc {e}\hspace{0.5pt}_\mathcal {P}) \in \mathpzc {F}_\mathcal {P}\times \mathpzc {E}_\mathcal {P}\) .

A Dirac structure linking (pH1)–(pH3),
\(\mathcal {D}\;&\subset \; T_z\mathcal {Z}\times T_z^*\mathcal {Z}\,\times \, \mathpzc {F}_\mathcal {R}\times \mathpzc {E}_\mathcal {R}\,\times \, \mathpzc {F}_\mathcal {P}\times \mathpzc {E}_\mathcal {P}.\) 

Definition C.1 [22]}
Consider a finite-dimensional linear space \(\mathpzc {F}\)  with dual \(\mathpzc {E}=\mathpzc {F}^*\) ,
and duality inner product \(\langle \,.,\,.\,\rangle :\mathpzc {F}\times \mathpzc {E}\rightarrow \mathbb {R}\) .
A linear subspace \(\mathcal {D}\subset \mathpzc {F}\times \mathpzc {E}\)  is a Dirac structure if \(\langle \,\mathpzc {f}\hspace{-1.0pt},\,\mathpzc {e}\hspace{0.5pt}\,\rangle = 0\) , for all \((\mathpzc {f}\hspace{-1.0pt},\mathpzc {e}\hspace{0.5pt}) \in \mathcal {D}\) , and \({\rm dim}\, \mathcal {D}= {\rm dim}\, \mathpzc {F}\) .

Here, using the natural identification \((V\times V^{\prime })^* = V^*\times V^{\prime *}\)  for linear spaces \({\hspace{0.2pt}\texttt {V}},V^{\prime }\) , the space on the right of (REF ) is viewed as the product \(\mathpzc {F}\times \mathpzc {E}\)  of
\(\mathpzc {F}&= T_z\mathcal {Z}\times \mathpzc {F}_\mathcal {R}\times \mathpzc {F}_\mathcal {P}\quad \text{and}\quad \mathpzc {E}= \mathpzc {F}^* = T_z^*\mathcal {Z}\times \mathpzc {E}_\mathcal {R}\times \mathpzc {E}_\mathcal {P},{with the inner product}\langle \,\mathpzc {f}\hspace{-1.0pt},\,\mathpzc {e}\hspace{0.5pt}\,\rangle  &= \langle \,(\mathpzc {f}\hspace{-1.0pt}_\mathcal {S},\mathpzc {f}\hspace{-1.0pt}_\mathcal {R},\mathpzc {f}\hspace{-1.0pt}_\mathcal {P}),\,(\mathpzc {e}\hspace{0.5pt}_\mathcal {S},\mathpzc {e}\hspace{0.5pt}_\mathcal {R},\mathpzc {e}\hspace{0.5pt}_\mathcal {P})\,\rangle  =\mathpzc {f}\hspace{-1.0pt}_\mathcal {S}\mathop {\cdot }\mathpzc {e}\hspace{0.5pt}_\mathcal {S} + \mathpzc {f}\hspace{-1.0pt}_\mathcal {R}\mathop {\cdot }\mathpzc {e}\hspace{0.5pt}_\mathcal {R} + \mathpzc {f}\hspace{-1.0pt}_\mathcal {P}\mathop {\cdot }\mathpzc {e}\hspace{0.5pt}_\mathcal {P} .\nonumber \) 

The \(\langle \,\mathpzc {f}\hspace{-1.0pt},\,\mathpzc {e}\hspace{0.5pt}\,\rangle = 0\)  in the definition represents the lossless routing of energy between parts of the system.
In the circuit context the Dirac structure is essentially Kirchhoff's laws for a given circuit topology.
That they define lossless routing, is Tellegen's theorem.
Port-Hamiltonian system dynamics
The dynamics of a system is specified by
\(\left(-\dot{z}(t), \frac{\partial \mathcal {H}}{\partial z}(z(t)),\; \mathpzc {f}\hspace{-1.0pt}_\mathcal {R}(t), \mathpzc {e}\hspace{0.5pt}_\mathcal {R}(t),\; \mathpzc {f}\hspace{-1.0pt}_\mathcal {P}(t), \mathpzc {e}\hspace{0.5pt}_\mathcal {P}(t) \right) &\in \mathcal {D}(z(t)), \\(\mathpzc {f}\hspace{-1.0pt}_\mathcal {R}(t), \mathpzc {e}\hspace{0.5pt}_\mathcal {R}(t)) &\in \mathcal {R}(z(t)), \nonumber \) 

for \(t\)  in some real interval.
A DAE that fits this structure—it might actually be an ODE—will be called a port-Hamiltonian DAE (pHDAE).
A key point from the physical modelling angle is that \(\mathcal {H}\)  is an algebraic function of \(z\) —the time derivative \(\dot{z}\)  must not enter—and must include the whole set of variables on which system energy depends (and preferably nothing but them).
Hence we can use the chain rule in the form
\(\dot{\mathcal {H}}(z(t)) = \frac{\partial \mathcal {H}}{\partial z}^(z(t))\;\dot{z}(t) = -\mathpzc {f}\hspace{-1.0pt}_\mathcal {S}\mathop {\cdot }\mathpzc {e}\hspace{0.5pt}_\mathcal {S},\) 

the second equality because from (REF ), \(-\dot{z}\)  is identified with \(\mathpzc {f}\hspace{-1.0pt}_\mathcal {S}\) , and \(\partial \mathcal {H}/\partial z\)  with \(\mathpzc {e}\hspace{0.5pt}_\mathcal {S}\) .
The Dirac structure enforces the power balance equation
\(\mathpzc {f}\hspace{-1.0pt}_\mathcal {S}\mathop {\cdot }\mathpzc {e}\hspace{0.5pt}_\mathcal {S} + \mathpzc {f}\hspace{-1.0pt}_\mathcal {R}\mathop {\cdot }\mathpzc {e}\hspace{0.5pt}_\mathcal {R} + \mathpzc {f}\hspace{-1.0pt}_\mathcal {P}\mathop {\cdot }\mathpzc {e}\hspace{0.5pt}_\mathcal {P} = 0\) 
and thus the DAE records energy flow: any solution of (REF ) satisfies
\(\dot{\mathcal {H}}= \mathpzc {f}\hspace{-1.0pt}_\mathcal {R}\mathop {\cdot }\mathpzc {e}\hspace{0.5pt}_\mathcal {R} + \mathpzc {f}\hspace{-1.0pt}_\mathcal {P}\mathop {\cdot }\mathpzc {e}\hspace{0.5pt}_\mathcal {P}.\) 

In particular, energy is conserved if there is no dissipation or external power flow.
Nothing is said about what variables occur in a pHDAE.
It will have the first order form \(f(t,x(t),\dot{x}(t)) = 0\) , where vector \(x(t)\) , the DAE variables, contains the state variables \(z(t)\in \mathcal {Z}\)  as a subvector.
What the remaining \(x_i\)  are, depends on the scenarios chosen to specify the dissipative structure (REF ) and the external port behaviour.

Positive definiteness and passivity results
Matrix positive definiteness
We cite without proof some results on positive definite matrices or matrix pairs that are standard, e.g. [24]}, [25]}, and prove some that are special to this paper. The facts about matrix pairs are present, but not so directly stated, in [19]}.
Definition D.1  
A matrix \(A\in \mathbb {R}^{n\times n}\)  is said to be positive definite (PD) if \(x^A̰x>0\)  for all \(x\in \mathbb {R}^n\setminus \lbrace 0\rbrace \) .
A pair of matrices \((B,B^{\prime })\) , \(B,B^{\prime }\in \mathbb {R}^{n\times n}\)  is called a positive definite matrix pair (PDMP) (resp. negative definite matrix pair (NDMP)) if for all \(y,y^{\prime }\in \mathbb {R}^n\)  not both 0, and \(By+B^{\prime }y^{\prime }=0\) , it holds that \(y^y̰^{\prime }>0\)  (resp. \(y^y̰^{\prime }<0\) ).

Note that positive definiteness of a matrix \(A\)  does not imply that \(A\)  is symmetric.
Lemma D.2  
Let \(A\in \mathbb {R}^{n\times n}\) .
[(i)]
If \(A\)  is PD, then \(A\)  is nonsingular, and \(A^{-1}\)  is PD.

\(A\)  is PD iff its symmetric part \((A+A^/2\)  is PD.

If \(A\)  is PD and \(\tilde{A}=-\tilde{A}^ is skew-symmetric, then \) A+A\( is PD.\)
Lemma D.3  
Let \((B,B^{\prime })\)  with \(B,B^{\prime }\in \mathbb {R}^{n\times n}\)  be a PDMP (resp. NDMP). Then
[(i)]
\((B^{\prime },B)\)  is a PDMP (resp. NDMP);

\(B\)  and \(B^{\prime }\)  are nonsingular and \(-B^{-1}B^{\prime }\) , \(-{B^{\prime }}^{-1}B\)  are PD (resp. \(B^{-1}B^{\prime }\) , \({B^{\prime }}^{-1}B\)  are PD);

\((PBQ,PB^{\prime }Q)\)  is a PDMP (resp. NDMP) for \(P,Q\in \mathbb {R}^{n\times n}\)  with \(P\)  nonsingular, \(Q\)  orthogonal.

\([C,C^{\prime }] = [B,B^{\prime }]\;P\)  is a PDMP (resp. NDMP), where \(P\)  is a permutation matrix as in (REF ), i.e., \((C,C^{\prime })\)  results from swapping any given set of columns of \(B\)  with the corresponding columns of \(B^{\prime }\) .
We prove the results for PDMP \((B,B^{\prime })\) .
[(i)]

Follows immediately from \({y^{\prime }}^y̰=y^y̰^{\prime }\) .

Suppose \(B\)  is singular, so there exists \(y\ne 0\)  with \(By=0\) . Setting \(y^{\prime }=0\)  gives \(By+B^{\prime }y^{\prime }=0\)  with \(y,y^{\prime }\)  not both 0, but \(y^y̰^{\prime }=0\) , contradicting the definition of a PDMP.
For the second part let \(A=-B^{-1}B^{\prime }\) . For any \(y^{\prime }\ne 0\) , set \(y=-B^{-1}B^{\prime }y^{\prime }\) .
Then \(y,y^{\prime }\)  are not both 0 and \(By+B^{\prime }y^{\prime }=0\)  so \(0<{y^{\prime }}^y̰ = -{y^{\prime }}^B̰^{-1}B^{\prime }y^{\prime } = {y^{\prime }}^A̰y^{\prime }\) . Hence \(A\)  is PD.

Let \((C,C^{\prime })=(PBQ,PB^{\prime }Q)\)  with \(P,Q\)  be as stated.
For any \(y,y^{\prime }\) , not both 0, suppose that \(0 = Cy+C^{\prime }y^{\prime } = P(Bz + B^{\prime }z^{\prime })\)  where \(z=Qy\) , \(z^{\prime }=Qy^{\prime }\) . Since \(P\)  and \(Q\)  are nonsingular, \(z,z^{\prime }\)  are not both 0 and \(Bz + B^{\prime }z^{\prime }=0\) .
As \((B,B^{\prime })\)  is a PDMP, we have \(0<z^z̰^{\prime } = y^Q̰^Q̰y^{\prime }=y^y̰^{\prime }\)  since \(Q\)  is orthogonal.
Since \(y,y^{\prime }\)  were arbitrary this shows \((C,C^{\prime })\)  that is a PDMP.

Let \(y,y^{\prime }\in \mathbb {R}^n\) , not both 0, satisfy \(Cy+Cy^{\prime }=0\) .
Defining \([w;w^{\prime }]=P_\beta [y;y^{\prime }]\) , we have \(Bw+B^{\prime }w^{\prime } = Cy+C^{\prime }y^{\prime }=0\) , with \(w,w^{\prime }\)  not both 0. So \(w^w̰^{\prime }>0\)  as \((B,B^{\prime })\)  is a PDMP.
By (REF ) it holds that \(w^w̰^{\prime }=y^y̰^{\prime }\)  (as \(w_jw_j^{\prime }\)  is replaced by \(w_j^{\prime }w_j\)  for each \(j\)  where a swap was done, but the sum stays the same).
Hence \(y^y̰^{\prime }>0\) , which proves \((C,C^{\prime })\)  is a PDMP.
Lemma D.4 
Let \(B,B^{\prime }\in \mathbb {R}^{n\times n}\) .
If \(B\)  is nonsingular and \(B^{-1}B^{\prime }\)  is PD then \((B,B^{\prime })\)  is a NDMP. In particular \(A\in \mathbb {R}^{n\times n}\)  is PD if and only if \((I,A)\)  is a NDMP.

Follows directly from the previous lemma.
Lemma D.5 
If \(M = \begin{bmatrix}M_{11} &M_{12} \\M_{21} &M_{22}\end{bmatrix}\begin{matrix}\mbox{\footnotesize $n$}_1 \\\mbox{\footnotesize $n$}_2\end{matrix}\)  is PD, then
\( P = \begin{bmatrix}M_{11} - N^M̰_{21}, &M_{12} - N^M̰_{22} \\N, &I\end{bmatrix} \) 

is nonsingular for an arbitrary \(n_2\times n_1\)  matrix \(N\) .

Take any \(w = [u;v]\)  and suppose \(Pw=0\) , so
\(0 &= (M_{11} - N^M̰_{21})u + (M_{12} - N^M̰_{22})v,\qquad 0 = Nu + v{Premultiply the first equation by u^ and substitute -Nu = v given by the second equation:}0 &= u^M̰_{11}u - u^N̰^M̰_{21}u + u^M̰_{12}v - u^N̰^M̰_{22}v \\&= u^M̰_{11}u + v^M̰_{21}u + u^M̰_{12}v + v^M̰_{22}v \qquad = w^M̰ w,\) 

so \(w=0\)  by the PD assumption, proving \(P\)  is nonsingular.
Lemma D.6 
Let \(C,C^{\prime },S\in \mathbb {R}^{n\times n}\) . If \((C,C^{\prime })\)  is a PDMP and \(S=-S^, then \) M =I S
C C'\( is nonsingular.\) 

Suppose \(z = \begin{bmatrix}x\\y\end{bmatrix}\)  with \(Mz = 0\) .
So \(x + Sy = 0\) , hence \(y^x̰ = y^-Sy)) = -y^S̰y = 0\)  since \(S\)  is skew-symmetric.
Also \(Cx + C^{\prime }y = 0\) , so as \((C,C^{\prime })\)  is a PDMP we have \(x = y = 0\) , i.e. \(z = 0\) . Hence \(M\)  is nonsingular.
Passivity results
There are various notions of passivity used in the literature. For a discussion of different definitions we refer to [27]}.
Most generally, a system is passive if it can either absorb or store energy, but cannot create energy within itself.
Passivity can be characterized by a dissipation inequality [28]}.
For port-Hamiltonian systems, passivity directly follows from the power-balance.
A port-Hamiltonian system is called passive (with respect to the port variables \(\mathpzc {f}\hspace{-1.0pt}_\mathcal {P},\mathpzc {e}\hspace{0.5pt}_\mathcal {P}\)  and storage function \(\mathcal {H}\) ) if
\(\mathcal {H}(z(t_1))-\mathcal {H}(z(t_0))\le \int _{t_0}^{t_1} \mathpzc {e}\hspace{0.5pt}_\mathcal {P}^_\mathcal {P}\, dt\) 
for all time instants \(t_0\le t_1\) .
Hence, a passive system cannot store more energy in its storage function \(\mathcal {H}\)  than it is supplied with.
We also need the concept of passivity for the dissipative structure \(\mathcal {R}\)  described by a relation \(r(\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}},\nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}})=0\) .
Definition D.7 [29]}
The resistive structure \(\mathcal {R}\)  is called passive [strictly passive] if \(\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\mathop {\cdot }\nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}} \ge 0\)  [resp. \(\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\mathop {\cdot }\nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}} > 0\) , except at the origin] for all consistent \((\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}, \nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}})\) .

Definition D.8 [30]}, [19]}
Assume that there exist a disjoint union \(\alpha \stackrel{\cdot }{\smash{\cup }}\beta ={\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\)  such that
\(\left.\left(\left[\frac{\partial r}{\partial \iota _{\alpha }},\frac{\partial r}{\partial \nu _{\beta }}\right]\right)\right|_{(\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}^*,\nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}^*)}\qquad \text{is nonsingular}\) 

at a given consistent \((\iota ^*_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}},\nu ^*_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}})\) .
The resistive structure is called strictly locally passive at \((\iota ^*_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}, \nu ^*_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}})\) if
\(D_{\cal R}:=-\begin{bmatrix}\frac{\partial r}{\partial \nu _{\alpha }}&\frac{\partial r}{\partial \iota _{\beta }}\end{bmatrix}^{-1}\cdot \begin{bmatrix}\frac{\partial r}{\partial \iota _{\alpha }}&\frac{\partial r}{\partial \nu _{\beta }}\end{bmatrix}.\) 

is positive definite at \((\iota ^*_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}, \nu ^*_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}})\) .
The resistive structure \(\mathcal {R}\)  is called strictly
locally passive if \(D_{\cal R}\)  is positive definite for all consistent \((\iota ^*_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}, \nu ^*_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}})\) .

Strict local passivity implies local passivity, the converse is not true.
Also passive components need not to be strictly locally passive as the following example shows.
Example 7 
The ideal transformer can be described by the relation
\( \begin{bmatrix}\nu _1\\\iota _2\end{bmatrix}=\begin{bmatrix}0&k\\-k&0\end{bmatrix}\begin{bmatrix}\iota _1\\\nu _2\end{bmatrix},\quad \text{where $k\in \mathbb {R}$ is the turns ratio}, \) 
see e.g. [32]}.
This resistive structure is passive, since
\(\nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\mathop {\cdot }\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}=\nu _1\iota _1+\nu _2\iota _2=k\nu _2\iota _1-k\nu _2\iota _1=0,\) 
but not strictly locally passive as (A5) is not satisfied.

From [19]} it follows that the positive definiteness of \(D_{\cal R}\)  in Definition REF  does not depend on the specific subsets \(\alpha ,\beta \) . In particular, in a strictly passive circuit one can choose arbitrarily the voltages or currents in every dissipative edge as a variable for a mixed description (REF ).
Defining
\(\bar{E}:=\frac{\partial r}{\partial \iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}}(\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}^*, \nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}^*), \text{ and } \bar{D}:=\frac{\partial r}{\partial \nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}}(\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}^*, \nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}^*)\) 
strict local passivity implies that the incremental conductance and resistance matrices
\(\bar{G}:=-\bar{E}^{-1}\bar{D} \quad \text{and}\quad \bar{R}:= -\bar{D}^{-1}\bar{E}\) 

are well-defined and positive definite (cf. [19]}).
In the linear setting (A5) reads
\(D \nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}+ E \iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}= 0\qquad \text{for some }D,E\in \mathbb {R}^{n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\times n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}}\) 

with positive definite pair \((D,E)\)  so that \(G=-E^{-1}D\)  (or equivalently \(R:= -D^{-1}E\) ) are positive definite.
Note that in the linear setting these properties are globally defined.
Within the port-Hamiltonian setting, the dissipative structure is defined by \(\mathcal {R}\subset \mathpzc {F}_\mathcal {R}\times \mathpzc {E}_\mathcal {R}\)  with \(\mathpzc {e}\hspace{0.5pt}_\mathcal {R}\mathop {\cdot }\mathpzc {f}\hspace{-1.0pt}_\mathcal {R}\le 0\) .
Non-negativity of the dissipative power flow \(\mathpzc {e}\hspace{0.5pt}_\mathcal {R}\mathop {\cdot }\mathpzc {f}\hspace{-1.0pt}_\mathcal {R}\)  respects the common convention that power flows from the “boundary” ports (i.e. ports related to \(\mathcal {P}\)  and \(\mathcal {R}\) ) into the system and from the internal network (the Dirac structure) into the energy storing elements \(\mathcal {S}\) , such that in the power balance equation all flows sum up to zero.
Here, it means the power flows from the energy-dissipating elements into the Dirac structure \(\mathcal {D}\)  (cf. [35]}).
By defining \(\nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}=\mathpzc {e}\hspace{0.5pt}_\mathcal {R}\)  and \(\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}=-\mathpzc {f}\hspace{-1.0pt}_\mathcal {R}\)  this setting corresponds to the sign convention in circuit analysis and in the definition of strict local passivity (e.g. the positive definiteness of conductance and resistance matrices) and
\(\mathpzc {e}\hspace{0.5pt}_\mathcal {R}\mathop {\cdot } \mathpzc {f}\hspace{-1.0pt}_\mathcal {R}\le 0 \quad \Longleftrightarrow \quad \iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\mathop {\cdot }\nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\ge 0. \) 
Using e.g. the voltage-controlled representation of resistors \(\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}= \bar{G} \nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\)  this gives
\(\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\mathop {\cdot }\nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}} =\nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\mathop {\cdot } \bar{G} \nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\ge 0\quad \Longleftrightarrow \quad \bar{G}\ge 0,\) 
i.e., the resistor is passive iff \(\bar{G}\)  is positive definite.

Matlab code
Below are the Matlab functions for the circuit in Example REF  and for the Diode Clipper example generated from the CpH class, both in Model 2 form. They can be used by Matlab's [basicstyle=]ode15i solver.
flexiblefalse0truesolarized@cyansolarized@bluesolarized@redsolarized@base01solarized@base2falsefalse%@trueutf8trueφ1
numbers=left]code/daefcnExample1mm.m
flexiblefalse0truesolarized@cyansolarized@bluesolarized@redsolarized@base01solarized@base2falsefalse%@trueutf8trueφ1
numbers=left]code/daefcnDiodeClippermm.m
