The external port \(\mathcal {P}\)  comprises only independent current and voltage sources.

To obtain a well-defined DAE system we append the implicit relation
\(r(\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}},\nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}) = 0\) 

as the constitutive relations for the dissipative elements to the system (REF ).
Due to M1-2, we can cast the DAE to input-output form as in
(REF , ).
For a given input \(u\) , i.e., given source functions, the state equation (REF ) is a well-defined DAE system that can be solved for \((t,x,\dot{x})\)  with a numerical integration scheme (cf. Section ).
Once a solution for \((t,x,\dot{x})\)  is obtained, the output \(y\)  can be computed using () if required.
Therefore, in the following we restrict our analysis to the CpH DAE
\(0=f(t,x,\dot{x},u)\) 

that is given by the state equation (REF ) of the control input-output form with a given input \(u\) .
In (REF ) the vector \(x\)  consists of the energy state variables \(z=(q,\varphi )\)  plus the \(\widehat{x}\)  components \(\nu _{\hspace{0.2pt}\texttt {d}}\) , \(\nu _{\hspace{0.2pt}\texttt {D}}\) , \(\iota _{\hspace{0.2pt}\texttt {d}}\) , \(\iota _{\hspace{0.2pt}\texttt {D}}\) .
This results in one DAE variable for each capacitor and inductor, and
two variables for each dissipative element in the circuit, for a total DAE size of
\(N = n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}+n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {L}}}+ 2 n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\) .
Theorem 4.1 
If the assumptions M1-1 and M1-2 hold, then the CpH DAE (REF ) is SA-amenable.

The proof follows standard lines of the \(\Sigma \) -method summarised in Appendix .
We construct the signature matrix of (REF ).
We choose suitable offsets, and show in Lemma REF —deploying (A3), (A4), (A5) and the normality of the tree—that the system Jacobian defined by these offsets is nonsingular.
The Theorem then follows.
Using (REF ) in (REF ) we get
\(0 = \left\lbrace \begin{array}{lc}\begin{bmatrix} f_{\hspace{0.2pt}\texttt {v}}\\f_{\hspace{0.2pt}\texttt {c}}\\ f_{\hspace{0.2pt}\texttt {d}}\\ f_{\hspace{0.2pt}\texttt {l}}\end{bmatrix} =\begin{bmatrix}\iota _{\hspace{0.2pt}\texttt {v}}\\ \iota _{\hspace{0.2pt}\texttt {c}}\\ \iota _{\hspace{0.2pt}\texttt {d}}\\ \iota _{\hspace{0.2pt}\texttt {l}}\end{bmatrix}- \begin{bmatrix}F_{{\hspace{0.2pt}\texttt {C}}{\hspace{0.2pt}\texttt {v}}}^F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {v}}}^F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {v}}}^F_{{\hspace{0.2pt}\texttt {I}}{\hspace{0.2pt}\texttt {v}}}^F_{{\hspace{0.2pt}\texttt {C}}{\hspace{0.2pt}\texttt {c}}}^F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {c}}}^F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {c}}}^F_{{\hspace{0.2pt}\texttt {I}}{\hspace{0.2pt}\texttt {c}}}^0&F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {d}}}^F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {d}}}^F_{{\hspace{0.2pt}\texttt {I}}{\hspace{0.2pt}\texttt {d}}}^0&0&F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {l}}}^F_{{\hspace{0.2pt}\texttt {I}}{\hspace{0.2pt}\texttt {l}}}^{bmatrix}\begin{bmatrix}\iota _{\hspace{0.2pt}\texttt {C}}\\ \iota _{\hspace{0.2pt}\texttt {D}}\\ \iota _{\hspace{0.2pt}\texttt {L}}\\ \iota _{\hspace{0.2pt}\texttt {I}}\end{bmatrix} \\[6ex]\begin{bmatrix}f_{\hspace{0.2pt}\texttt {C}}\\ f_{\hspace{0.2pt}\texttt {D}}\\f_{\hspace{0.2pt}\texttt {L}}\\ f_{\hspace{0.2pt}\texttt {I}}\end{bmatrix} =\begin{bmatrix}\nu _{\hspace{0.2pt}\texttt {C}}\\ \nu _{\hspace{0.2pt}\texttt {D}}\\ \nu _{\hspace{0.2pt}\texttt {L}}\\ \nu _{\hspace{0.2pt}\texttt {I}}\end{bmatrix}+ \begin{bmatrix}F_{{\hspace{0.2pt}\texttt {C}}{\hspace{0.2pt}\texttt {v}}} &F_{{\hspace{0.2pt}\texttt {C}}{\hspace{0.2pt}\texttt {c}}} &0 &0 \\F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {v}}} &F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {c}}} &F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {d}}} &0 \\F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {v}}} &F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {c}}} &F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {d}}} &F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {l}}}\\F_{{\hspace{0.2pt}\texttt {I}}{\hspace{0.2pt}\texttt {v}}} &F_{{\hspace{0.2pt}\texttt {I}}{\hspace{0.2pt}\texttt {c}}} &F_{{\hspace{0.2pt}\texttt {I}}{\hspace{0.2pt}\texttt {d}}} &F_{{\hspace{0.2pt}\texttt {I}}{\hspace{0.2pt}\texttt {l}}}\end{bmatrix}\begin{bmatrix}\nu _{\hspace{0.2pt}\texttt {v}}\\ \nu _{\hspace{0.2pt}\texttt {c}}\\ \nu _{\hspace{0.2pt}\texttt {d}}\\ \nu _{\hspace{0.2pt}\texttt {l}}\end{bmatrix} \\[6ex]r(\iota _{\hspace{0.2pt}\texttt {d}}, \nu _{\hspace{0.2pt}\texttt {D}}, \nu _{\hspace{0.2pt}\texttt {d}}, \iota _{\hspace{0.2pt}\texttt {D}})\end{bmatrix}\end{array}\right.We remove the f_{\hspace{0.2pt}\texttt {v}} and f_{\hspace{0.2pt}\texttt {I}} rows in (\ref {eq:CpH_2}), as they give(\ref {eq:ctrlIOform4b}),set \nu _{\hspace{0.2pt}\texttt {v}}= V(t),  \iota _{\hspace{0.2pt}\texttt {I}}= I(t), and using(\ref {eq:zconstiteqns}) write theresulting equations constituting (\ref {eq:CpH_Model1}) as\) 

\(\begin{array}{l@{\;=\;}|c|c|c|c}\multicolumn{1}{c}{} &\multicolumn{1}{c}{\mbox{\footnotesize $q=(q_{\hspace{0.2pt}\texttt {C}}, q_{\hspace{0.2pt}\texttt {c}})$}} & \multicolumn{1}{c}{\mbox{\footnotesize $\varphi =(\varphi _{\hspace{0.2pt}\texttt {l}}, \varphi _{\hspace{0.2pt}\texttt {L}})$}} & \multicolumn{1}{c}{\mbox{\footnotesize $\widehat{x}=(\iota _{\hspace{0.2pt}\texttt {d}}, \nu _{\hspace{0.2pt}\texttt {D}}, \nu _{\hspace{0.2pt}\texttt {d}}, \iota _{\hspace{0.2pt}\texttt {D}})$}} & \\\cline {2-4}f_{\hspace{0.2pt}\texttt {C}}&\frac{\partial \mathcal {H}}{\partial q_{\hspace{0.2pt}\texttt {C}}}(q_{\hspace{0.2pt}\texttt {c}},q_{\hspace{0.2pt}\texttt {C}}) + F_{{\hspace{0.2pt}\texttt {C}}{\hspace{0.2pt}\texttt {c}}} \frac{\partial \mathcal {H}}{\partial q_{\hspace{0.2pt}\texttt {c}}}(q_{\hspace{0.2pt}\texttt {c}},q_{\hspace{0.2pt}\texttt {C}}) &&&+ F_{{\hspace{0.2pt}\texttt {C}}{\hspace{0.2pt}\texttt {v}}} V(t) \\f_{\hspace{0.2pt}\texttt {c}}&- F_{{\hspace{0.2pt}\texttt {C}}{\hspace{0.2pt}\texttt {c}}}^_{\hspace{0.2pt}\texttt {C}}+ \dot{q}_{\hspace{0.2pt}\texttt {c}}&- F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {c}}}^{\mathcal {H}}{\varphi _{\hspace{0.2pt}\texttt {L}}}(\varphi _{\hspace{0.2pt}\texttt {l}},\varphi _{\hspace{0.2pt}\texttt {L}}) &- F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {c}}}^_{\hspace{0.2pt}\texttt {D}}& - F_{{\hspace{0.2pt}\texttt {I}}{\hspace{0.2pt}\texttt {c}}}^Ḭ(t) \\ \cline {2-4}f_{\hspace{0.2pt}\texttt {l}}&&\frac{\partial \mathcal {H}}{\partial \varphi _{\hspace{0.2pt}\texttt {l}}}(\varphi _{\hspace{0.2pt}\texttt {l}},\varphi _{\hspace{0.2pt}\texttt {L}}) - F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {l}}}^{\mathcal {H}}{\varphi _{\hspace{0.2pt}\texttt {L}}}(\varphi _{\hspace{0.2pt}\texttt {l}},\varphi _{\hspace{0.2pt}\texttt {L}}) &&- F_{{\hspace{0.2pt}\texttt {I}}{\hspace{0.2pt}\texttt {l}}}^Ḭ(t) \\f_{\hspace{0.2pt}\texttt {L}}& F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {c}}}\frac{\partial \mathcal {H}}{\partial q_{\hspace{0.2pt}\texttt {c}}}(q_{\hspace{0.2pt}\texttt {c}},q_{\hspace{0.2pt}\texttt {C}}) &+ F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {l}}} \dot{\varphi }_{\hspace{0.2pt}\texttt {l}}+ \dot{\varphi }_{\hspace{0.2pt}\texttt {L}}&+ F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {d}}} \nu _{\hspace{0.2pt}\texttt {d}}&+ F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {v}}} V(t) \\\cline {2-4}f_{\hspace{0.2pt}\texttt {d}}&&- F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {d}}}^{\mathcal {H}}{\varphi _{\hspace{0.2pt}\texttt {L}}}(\varphi _{\hspace{0.2pt}\texttt {l}},\varphi _{\hspace{0.2pt}\texttt {L}}) &+ \iota _{\hspace{0.2pt}\texttt {d}}- F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {d}}}^_{\hspace{0.2pt}\texttt {D}}&- F_{{\hspace{0.2pt}\texttt {I}}{\hspace{0.2pt}\texttt {d}}}^Ḭ(t) \\f_{\hspace{0.2pt}\texttt {D}}& F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {c}}} \frac{\partial \mathcal {H}}{\partial q_{\hspace{0.2pt}\texttt {c}}}(q_{\hspace{0.2pt}\texttt {c}},q_{\hspace{0.2pt}\texttt {C}}) &&+\nu _{\hspace{0.2pt}\texttt {D}}+ F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {d}}} \nu _{\hspace{0.2pt}\texttt {d}}&+ F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {v}}} V(t) \\r &&&{r(\iota _{\hspace{0.2pt}\texttt {d}}, \nu _{\hspace{0.2pt}\texttt {D}}, \nu _{\hspace{0.2pt}\texttt {d}}, \iota _{\hspace{0.2pt}\texttt {D}})} & \\ \cline {2-4}\end{array}\) 


The corresponding schematic \(\Sigma \)  matrix is
\(\Sigma \ &=\ \begin{blockarray}{rcccccccclc}& q_{\hspace{0.2pt}\texttt {C}}& q_{\hspace{0.2pt}\texttt {c}}& \varphi _{\hspace{0.2pt}\texttt {l}}& \varphi _{\hspace{0.2pt}\texttt {L}}& \iota _{\hspace{0.2pt}\texttt {d}}& \nu _{\hspace{0.2pt}\texttt {D}}& \nu _{\hspace{0.2pt}\texttt {d}}& \iota _{\hspace{0.2pt}\texttt {D}}& \mbox{\footnotesize $c_i$} & \\\begin{block}{r[cc|cc|cccc]ll}f_{\hspace{0.2pt}\texttt {C}}& \colorbox {lightgray!110}{\le 0} & \colorbox {lightgray!110}{\le 0} & - & - & - & - & - & - & \mbox{\footnotesize $1$} & n_{\hspace{0.2pt}\texttt {C}}\\f_{\hspace{0.2pt}\texttt {c}}& \colorbox {lightgray!110}{\le 1} & \colorbox {lightgray!110}{\widehat{I}{+}1} & {\le 0} & {\le 0} & - & - & - & {\le 0} & \mbox{\footnotesize $0$} & n_{\hspace{0.2pt}\texttt {c}}\\ \cline {2-9}f_{\hspace{0.2pt}\texttt {l}}& - & - & \colorbox {lightgray!110}{\le 0} & \colorbox {lightgray!110}{\le 0} & - & - & - & - & \mbox{\footnotesize $1$} & n_{\hspace{0.2pt}\texttt {l}}\\f_{\hspace{0.2pt}\texttt {L}}& {\le 0} & {\le 0} & \colorbox {lightgray!110}{\le 1} & \colorbox {lightgray!110}{\widehat{I}{+}1} & - & - & {\le 0} & - & \mbox{\footnotesize $0$} & n_{\hspace{0.2pt}\texttt {L}}\\ \cline {2-9}f_{\hspace{0.2pt}\texttt {d}}& - & - & \le 0 & \le 0 &\colorbox {lightgray!110}{\widehat{I}} & - & - & \colorbox {lightgray!110}{\le 0} & \mbox{\footnotesize $1$} & n_{\hspace{0.2pt}\texttt {d}}\\f_{\hspace{0.2pt}\texttt {D}}& \le 0 & \le 0 & - & - &- & \colorbox {lightgray!110}{\widehat{I}} & \colorbox {lightgray!110}{\le 0} & - & \mbox{\footnotesize $1$} & n_{\hspace{0.2pt}\texttt {D}}\\r & - & - & - & - &\colorbox {lightgray!110}{\le 0} & \colorbox {lightgray!110}{\le 0} & \colorbox {lightgray!110}{\le 0} & \colorbox {lightgray!110}{\le 0} & \mbox{\footnotesize $1$} & n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}{=}n_{\hspace{0.2pt}\texttt {D}}{+}n_{\hspace{0.2pt}\texttt {d}}\\\end{block}\mbox{\footnotesize $d_j$}& \mbox{\footnotesize $1$} & \mbox{\footnotesize $1$} & \mbox{\footnotesize $1$} & \mbox{\footnotesize $1$} & \mbox{\footnotesize $1$} & \mbox{\footnotesize $1$} &\mbox{\footnotesize $1$} & \mbox{\footnotesize $1$} & & & \\& n_{\hspace{0.2pt}\texttt {C}}& n_{\hspace{0.2pt}\texttt {c}}& n_{\hspace{0.2pt}\texttt {l}}& n_{\hspace{0.2pt}\texttt {L}}& n_{\hspace{0.2pt}\texttt {d}}& n_{\hspace{0.2pt}\texttt {D}}& n_{\hspace{0.2pt}\texttt {d}}& n_{\hspace{0.2pt}\texttt {D}}& &\end{blockarray}\) 

All entries \(\sigma _{ij}\)  are \({-\!\infty }\) , 0 or 1 because at most first derivatives occur in the DAE.
Here \(\le 0\)  denotes a block with 0 or \({-\!\infty }\)  entries; \(\le 1\)  denotes a block with 1 or \({-\!\infty }\)  entries; \(-\)  is a block of \({-\!\infty }\)  entries; \(\widehat{I}\)  is a square block with 0's on the main diagonal and \({-\!\infty }\)  elsewhere; \(\widehat{I}+1\)  (elementwise addition) is 1 on the main diagonal and \({-\!\infty }\)  elsewhere.
We assign provisional offset values \(c_i, d_j\) .
Their meaning is blockwise, e.g. \(c_i=0\)  everywhere in the \(f_{\hspace{0.2pt}\texttt {c}}\)  and \(f_{\hspace{0.2pt}\texttt {L}}\)  rows and \(c_i=1\)  elsewhere.
Lemma 4.2 
The system Jacobian \(\mathbf {J}\)  defined by the provisional offsets is nonsingular.

\(\mathbf {J}\)  is defined by (REF ) in Appendix .
The blocks with light grey background in (REF ) are those where \(d_j-c_i>\sigma _{ij}\) , and therefore
\(\mathbf {J}_{ij}=0\)  in these blocks. Hence \(\mathbf {J}\)  has block-triangular form
\(\mathbf {J} \ &=\ \begin{blockarray}{rcccl}& q_{\hspace{0.2pt}\texttt {C}},q_{\hspace{0.2pt}\texttt {c}}& \varphi _{\hspace{0.2pt}\texttt {l}},\varphi _{\hspace{0.2pt}\texttt {L}}& \iota _{\hspace{0.2pt}\texttt {d}},\nu _{\hspace{0.2pt}\texttt {D}},\nu _{\hspace{0.2pt}\texttt {d}},\iota _{\hspace{0.2pt}\texttt {D}}& \\\begin{block}{r[c|c|c]l}f_{\hspace{0.2pt}\texttt {C}},f_{\hspace{0.2pt}\texttt {c}}& \colorbox {lightgray}{\mathbf {J}_{\mathcal {C}}} & 0 & 0 &\mbox{\footnotesize $n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}=n_{\hspace{0.2pt}\texttt {C}}{+}n_{\hspace{0.2pt}\texttt {c}}$} \\ \cline {2-4}f_{\hspace{0.2pt}\texttt {l}},f_{\hspace{0.2pt}\texttt {L}}& 0 & \colorbox {lightgray}{\mathbf {J}_{\mathcal {L}}} & 0 &\mbox{\footnotesize $n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {L}}}=n_{\hspace{0.2pt}\texttt {l}}{+}n_{\hspace{0.2pt}\texttt {L}}$} \\\cline {2-4}f_{\hspace{0.2pt}\texttt {d}},f_{\hspace{0.2pt}\texttt {D}},r & \ast & \ast &\colorbox {lightgray}{\mathbf {J}_{\mathcal {D}}} &\mbox{\footnotesize $2n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}=n_{\hspace{0.2pt}\texttt {d}}{+}n_{\hspace{0.2pt}\texttt {D}}{+}n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}$} \\\end{block}& \mbox{\footnotesize $n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}$} & \mbox{\footnotesize $n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {L}}}$} & \mbox{\footnotesize $2n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}$}\end{blockarray}\) 

where \(\ast \)  marks a possibly nonzero block.
It suffices to show each of \(\mathbf {J}_{\mathcal {C}},\mathbf {J}_{\mathcal {L}},\mathbf {J}_{\mathcal {D}}\)  is nonsingular.
The \(\mathbf {J}_{\mathcal {C}}\)  block.
Since \(d_j-c_i\)  equals 0 in the \(f_{\hspace{0.2pt}\texttt {C}}\)  rows and 1 in the \(f_{\hspace{0.2pt}\texttt {c}}\)  rows, by (REF )
\(\mathbf {J}_{\mathcal {C}} &= \begin{bmatrix}\displaystyle \frac{\partial f_{\hspace{0.2pt}\texttt {C}}}{\partial q_{\hspace{0.2pt}\texttt {C}}}& \displaystyle \frac{\partial f_{\hspace{0.2pt}\texttt {C}}}{\partial q_{\hspace{0.2pt}\texttt {c}}}\\[2ex]\displaystyle \frac{\partial f_{\hspace{0.2pt}\texttt {c}}}{\partial \dot{q}_{\hspace{0.2pt}\texttt {C}}}& \displaystyle \frac{\partial f_{\hspace{0.2pt}\texttt {c}}}{\partial \dot{q}_{\hspace{0.2pt}\texttt {c}}}\end{bmatrix}= \begin{bmatrix}\displaystyle \frac{\partial ^2 \mathcal {H}}{\partial q_{\hspace{0.2pt}\texttt {C}}^2} + F_{{\hspace{0.2pt}\texttt {C}}{\hspace{0.2pt}\texttt {c}}} \frac{\partial ^2\mathcal {H}}{\partial q_{\hspace{0.2pt}\texttt {C}}\partial q_{\hspace{0.2pt}\texttt {c}}} & \displaystyle \frac{\partial ^2\mathcal {H}}{\partial q_{\hspace{0.2pt}\texttt {c}}\partial q_{\hspace{0.2pt}\texttt {C}}} + F_{{\hspace{0.2pt}\texttt {C}}{\hspace{0.2pt}\texttt {c}}} \frac{\partial ^2 \mathcal {H}}{\partial q_{\hspace{0.2pt}\texttt {c}}^2} \\[3ex]-F_{{\hspace{0.2pt}\texttt {C}}{\hspace{0.2pt}\texttt {c}}}^ I\end{bmatrix} \\\nonumber \) 

Here, \(I\)  denotes an identity matrix of suitable size.
By (A4), the matrix
\(\mathbf {C}& = \begin{bmatrix}\displaystyle \frac{\partial ^2 \mathcal {H}}{\partial q_{\hspace{0.2pt}\texttt {C}}^2} & \displaystyle \frac{\partial ^2\mathcal {H}}{\partial q_{\hspace{0.2pt}\texttt {C}}\partial q_{\hspace{0.2pt}\texttt {c}}} \\[3ex]\displaystyle \frac{\partial ^2\mathcal {H}}{\partial q_{\hspace{0.2pt}\texttt {c}}\partial q_{\hspace{0.2pt}\texttt {C}}}& \displaystyle \frac{\partial ^2 \mathcal {H}}{\partial q_{\hspace{0.2pt}\texttt {c}}^2} \end{bmatrix}\) 

is PD. Setting \(M = \mathbf {C}\)  and \(N=-F_{{\hspace{0.2pt}\texttt {C}}{\hspace{0.2pt}\texttt {c}}}^ in Lemma~\ref {lm:PD2}, \) JC \( is of the form (\ref {eq:PD2}) and hence nonsingular.\)
The \(\mathbf {J}_{\cal L}\)  block.
Similarly,
\(\mathbf {J}_{\cal L} &=\begin{bmatrix}\displaystyle \frac{\partial f_{\hspace{0.2pt}\texttt {l}}}{\partial \varphi _{\hspace{0.2pt}\texttt {l}}}& \displaystyle \frac{\partial f_{\hspace{0.2pt}\texttt {l}}}{\partial \varphi _{\hspace{0.2pt}\texttt {L}}}\\[2ex] \displaystyle \frac{\partial f_{\hspace{0.2pt}\texttt {L}}}{\partial \dot{\varphi }_{\hspace{0.2pt}\texttt {l}}}& \displaystyle \frac{\partial f_{\hspace{0.2pt}\texttt {L}}}{\partial \dot{\varphi }_{\hspace{0.2pt}\texttt {L}}}\end{bmatrix}=\begin{bmatrix}\displaystyle \frac{\partial ^2 \mathcal {H}}{\partial \varphi _{\hspace{0.2pt}\texttt {l}}^2} -F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {l}}}^{\mathcal {H}}{\varphi _{\hspace{0.2pt}\texttt {L}}}{\varphi _{\hspace{0.2pt}\texttt {l}}} &\displaystyle \frac{\partial ^2\mathcal {H}}{\partial \varphi _{\hspace{0.2pt}\texttt {l}}\partial \varphi _{\hspace{0.2pt}\texttt {L}}} -F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {l}}}^{^2 \mathcal {H}}{\varphi _{\hspace{0.2pt}\texttt {L}}^2}\\[3ex]F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {l}}} & I\end{bmatrix}.\) 

By (A5) the matrix
\(\mathbf {L}& = \begin{bmatrix}\displaystyle \frac{\partial ^2 \mathcal {H}}{\partial \varphi _{\hspace{0.2pt}\texttt {l}}^2} & \displaystyle \frac{\partial ^2\mathcal {H}}{\partial \varphi _{\hspace{0.2pt}\texttt {l}}\partial \varphi _{\hspace{0.2pt}\texttt {L}}} \\[3ex]\displaystyle \frac{\partial ^2\mathcal {H}}{\partial \varphi _{\hspace{0.2pt}\texttt {L}}\partial \varphi _{\hspace{0.2pt}\texttt {l}}}& \displaystyle \frac{\partial ^2 \mathcal {H}}{\partial \varphi _{\hspace{0.2pt}\texttt {L}}^2} \end{bmatrix}\) 

is PD.
Setting \(M=\mathbf {L}\)  and \(N=F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {l}}}\)  in Lemma REF , \(\mathbf {J}_{\cal L}\)  is of the form (REF ) and hence nonsingular.
The \(\mathbf {J}_\mathcal {D}\)  block.
We have
\(\frac{\partial (f_{\hspace{0.2pt}\texttt {d}}, f_{\hspace{0.2pt}\texttt {D}})}{\partial (\iota _{\hspace{0.2pt}\texttt {d}}, \nu _{\hspace{0.2pt}\texttt {D}}, \nu _{\hspace{0.2pt}\texttt {d}}, \iota _{\hspace{0.2pt}\texttt {D}})} =\begin{bmatrix}I& 0 & 0 &-F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {d}}}^0 & I& F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {d}}} & 0 \end{bmatrix} =: \begin{bmatrix}I&S\end{bmatrix} \in \mathbb {R}^{n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\times 2n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}},\) 
where \(S\)  is skew-symmetric.
By (A6), \((H,H^{\prime }) {=} \bigl (\partial r/\partial (\nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}), \partial r/\partial (\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}})\bigr ) {=} \bigl (\partial r/\partial (\nu _{\hspace{0.2pt}\texttt {d}}, \nu _{\hspace{0.2pt}\texttt {D}}), \partial r/\partial (\iota _{\hspace{0.2pt}\texttt {d}}, \iota _{\hspace{0.2pt}\texttt {D}})\bigr )\)  is a positive definite matrix pair (PDMP).
By Lemma REF  (iv), \((\mathcal {D},\mathcal {D}^{\prime }) = \bigl (\partial r/\partial (\iota _{\hspace{0.2pt}\texttt {d}}, \nu _{\hspace{0.2pt}\texttt {D}}),\linebreak \partial r/\partial (\nu _{\hspace{0.2pt}\texttt {d}}, \iota _{\hspace{0.2pt}\texttt {D}})\bigr )\) , which is obtained by swapping some columns of \(H\)  with corresponding columns of \(H^{\prime }\) , is also a PDMP.
Thus Jacobian \(\mathbf {J}_\mathcal {D}\)  has the form
\(\mathbf {J}_\mathcal {D}&={\begin{bmatrix}\displaystyle \frac{\partial f_{\hspace{0.2pt}\texttt {d}}}{\partial \iota _{\hspace{0.2pt}\texttt {d}}}&\displaystyle \frac{\partial f_{\hspace{0.2pt}\texttt {d}}}{\partial \nu _{\hspace{0.2pt}\texttt {D}}}& \displaystyle \frac{\partial f_{\hspace{0.2pt}\texttt {d}}}{\partial \iota _{\hspace{0.2pt}\texttt {D}}}& \displaystyle \frac{\partial f_{\hspace{0.2pt}\texttt {d}}}{\partial \nu _{\hspace{0.2pt}\texttt {d}}}\\\displaystyle \frac{\partial f_{\hspace{0.2pt}\texttt {D}}}{\partial \iota _{\hspace{0.2pt}\texttt {d}}}&\displaystyle \frac{\partial f_{\hspace{0.2pt}\texttt {D}}}{\partial \nu _{\hspace{0.2pt}\texttt {D}}}&\displaystyle \frac{\partial f_{\hspace{0.2pt}\texttt {D}}}{\partial \iota _{\hspace{0.2pt}\texttt {D}}}&\displaystyle \frac{\partial f_{\hspace{0.2pt}\texttt {D}}}{\partial \nu _{\hspace{0.2pt}\texttt {d}}}\\[1ex]\cline {1-4}\displaystyle \frac{\partial r}{\partial \iota _{\hspace{0.2pt}\texttt {d}}}&\displaystyle \frac{\partial r}{\partial \nu _{\hspace{0.2pt}\texttt {D}}}& \displaystyle \frac{\partial r}{\partial \iota _{\hspace{0.2pt}\texttt {D}}}& \displaystyle \frac{\partial r}{\partial \nu _{\hspace{0.2pt}\texttt {d}}}\end{bmatrix}\begin{matrix}n_{\hspace{0.2pt}\texttt {d}}\\n_{\hspace{0.2pt}\texttt {D}}\\n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\end{matrix}= \begin{bmatrix}I&S \\\cline {1-2} \mathcal {D}&\mathcal {D}^{\prime }\end{bmatrix}\begin{matrix}n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\\n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\end{matrix}}\) 

satisfying the conditions of Lemma REF  and hence is nonsingular.
Proof of Theorem REF .
Lemma REF  has shown that the Jacobian \(\mathbf {J}\)  defined by the provisional offsets \(c_i\) , \(d_j\)  given in (REF ) is nonsingular.
Thus, there exists a transversal \(\mathcal {T}\)  such that \(\mathbf {J}_{ij}\ne 0\)  for all \((i,j)\in \mathcal {T}\) , which means the DAE (REF ) is structurally well-posed (has a finite transversal).
Namely \(d_j-c_i = \sigma _{ij}\)  at all such \((i,j)\) , see (REF ), and by [1]}, \(\mathcal {T}\)  is an HVT.
The value of \(\mathcal {T}\)  is
\(\mathop {\textup {Val}}\mathcal {T}= \sum _{(i,j)\in \mathcal {T}}\sigma _{ij}= \sum _j d_j - \sum _i c_i = n_{\hspace{0.2pt}\texttt {c}}+n_{\hspace{0.2pt}\texttt {L}},\) 

which is also the number of degrees of freedom.
Since the offsets \(c_i,d_j\)  satisfy (REF ) in Appendix  they are valid offsets.
Hence, the DAE (REF ) is SA-amenable.
\(\Box \)

Theorem 4.3 
Assume that M1-1 and M1-2 hold.
If \(n_{\hspace{0.2pt}\texttt {C}}=n_{\hspace{0.2pt}\texttt {l}}=n_{\hspace{0.2pt}\texttt {d}}=n_{\hspace{0.2pt}\texttt {D}}=0\)  the system (REF ) is an ODE.
Otherwise, (REF ) is a DAE of index 1.

If \(n_{\hspace{0.2pt}\texttt {C}}=n_{\hspace{0.2pt}\texttt {l}}=n_{\hspace{0.2pt}\texttt {d}}=n_{\hspace{0.2pt}\texttt {D}}=0\)  it is immediately clear that the system (REF ) simplifies to an ODE,
as (REF ) does not contain algebraic equations.
In the other case—at least one of \(n_{\hspace{0.2pt}\texttt {C}}, n_{\hspace{0.2pt}\texttt {l}}, n_{\hspace{0.2pt}\texttt {d}}, n_{\hspace{0.2pt}\texttt {D}}\)  does not vanish—since the DAE is SA-amenable by Theorem REF  we can apply the standard solution scheme (SSS), see (REF ) in § for solving (REF ).
By Theorem REF  the Jacobian \(\mathbf {J}_k\)  is of full row rank for \(k=-1,0\)  such that the generally underdetermined system in stage \(k=-1\)  can be solved for \(x\)  if \(\textsc {dof} = n_{\hspace{0.2pt}\texttt {c}}+n_{\hspace{0.2pt}\texttt {L}}\)  initial conditions are provided.
In stage \(k=0\)  we get a nonsingular system for all of \(\dot{x}\) .
Thus, one differentiation of equations is needed to convert the original DAE to an implicit ODE, and the DAE has index 1.
Corollary 4.4 The following statements are equivalent:
[(i)]
system (REF ) is an ODE;

the circuit contains neither dissipators \({\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\) , \({{\hspace{0.2pt}{\hspace{0.2pt}\texttt {L}}}{\hspace{0.2pt}{\hspace{0.2pt}\texttt {I}}}}\) -cutsets or \({{\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}{\hspace{0.2pt}{\hspace{0.2pt}\texttt {V}}}}\) -cycles;

\(A_{{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}}=\emptyset \)  and \(\mathop {\rm rank}\nolimits [A_{{\hspace{0.2pt}{\hspace{0.2pt}\texttt {V}}}}\;A_{{\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}}]={{n}{-}1}\)  as well as \(\ker [A_{{\hspace{0.2pt}{\hspace{0.2pt}\texttt {V}}}}\;A_{{\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}}]=\lbrace 0\rbrace \) .
Under (A3) the circuit contains \({{\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}{\hspace{0.2pt}{\hspace{0.2pt}\texttt {V}}}}\) -cycles if and only if \(n_{\hspace{0.2pt}\texttt {C}}>0\)  (since otherwise components from \({\hspace{0.2pt}\texttt {C}}\)  could be swapped with other component from the tree contradicting the optimality of (A3)). The same argument holds for \({{\hspace{0.2pt}{\hspace{0.2pt}\texttt {L}}}{\hspace{0.2pt}{\hspace{0.2pt}\texttt {I}}}}\) -cutsets.
The corresponding rank conditions on \(A\)  follow from Theorem REF .
Remark 4.5 If one is interested in controlling the circuit, i.e., the input \(u\)  is treated as an unknown itself, then the whole system (REF ) has to be considered.
The same holds if the interconnection of different port-Hamiltonian systems via inputs and outputs is taken into account.
The result in Theorem REF  also holds for this more general setting, although with different offsets and index results.
Since we are interested in the simulation (and not control) of the circuit equations, we do not present these more general results here.

Remark 4.6 The offsets provided in (REF ) are shown to be valid offsets, however they are non-canonical, i.e., not the smallest ones possible (see §).
Using these offsets from (REF ) we get a structural index \(\nu _S\)  which corresponds to the results of Theorem REF : \(\nu _S=0\)  if \(n_{\hspace{0.2pt}\texttt {C}}=n_{\hspace{0.2pt}\texttt {l}}=n_{\hspace{0.2pt}\texttt {d}}=n_{\hspace{0.2pt}\texttt {D}}=0\) , and \(\nu _S=1\)  otherwise.
If we use the canonical offsets, i.e., \(c_i=0\)  and \(d_j=0\)  in the rows/columns belonging to the dissipative structure, we get a structural index which overestimates the index by one.

CpH Circuit Model 2
For the CpH Circuit Model 2 we make the following assumptions.
[M2-1.]
The assumptions (A1)–(A6) hold.

The external port \(\mathcal {P}\)  comprises only independent current and voltage sources.

The dissipative structure is described by a mixed form as in §REF , (REF ).

As in §REF , we cast the DAE into an input-output form and restrict to the CpH DAE
\(0=f(t,x,\dot{x},u),\) 

where \(u\)  is given input.
Port \(\mathcal {P}\)  contributes no DAE variables, so Algorithm REF  has \(x=(z,\widehat{x})\)  with energy state variables \(z=(q,\varphi )\) , dissipative variables \(\widehat{x}\) —one variable for each dissipative edge, for a total DAE size of
\(N = n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}+n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {L}}}+ n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}= {b}- n_\mathcal {P},\) 

i.e., the total number of branches that are not external ports.
Theorem 4.7 If the assumptions M2-1–M2-3 hold, then the CpH DAE (REF ) is SA-amenable.

The proof is similar to that of Theorem REF .
The version of (REF ) for Model 2 is

\(\begin{array}{l@{\;=\;}|c|c|c|c}\multicolumn{1}{c}{} &\multicolumn{1}{c}{\mbox{\footnotesize $q=(q_{\hspace{0.2pt}\texttt {C}}, q_{\hspace{0.2pt}\texttt {c}})$}} & \multicolumn{1}{c}{\mbox{\footnotesize $\varphi =(\varphi _{\hspace{0.2pt}\texttt {l}}, \varphi _{\hspace{0.2pt}\texttt {L}})$}} & \multicolumn{1}{c}{\mbox{\footnotesize $\widehat{x}=(\iota _{\hspace{0.2pt}\texttt {d}}, \nu _{\hspace{0.2pt}\texttt {D}})$}} & \\\cline {2-4}f_{\hspace{0.2pt}\texttt {C}}&\frac{\partial \mathcal {H}}{\partial q_{\hspace{0.2pt}\texttt {C}}}(q_{\hspace{0.2pt}\texttt {c}},q_{\hspace{0.2pt}\texttt {C}}) {+} F_{{\hspace{0.2pt}\texttt {C}}{\hspace{0.2pt}\texttt {c}}} \frac{\partial \mathcal {H}}{\partial q_{\hspace{0.2pt}\texttt {c}}}(q_{\hspace{0.2pt}\texttt {c}},q_{\hspace{0.2pt}\texttt {C}}) &&&+ F_{{\hspace{0.2pt}\texttt {C}}{\hspace{0.2pt}\texttt {v}}} V(t) \\f_{\hspace{0.2pt}\texttt {c}}&- F_{{\hspace{0.2pt}\texttt {C}}{\hspace{0.2pt}\texttt {c}}}^_{\hspace{0.2pt}\texttt {C}}+ \dot{q}_{\hspace{0.2pt}\texttt {c}}&- F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {c}}}^{\mathcal {H}}{\varphi _{\hspace{0.2pt}\texttt {L}}}(\varphi _{\hspace{0.2pt}\texttt {l}},\varphi _{\hspace{0.2pt}\texttt {L}}) &- F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {c}}}^_{\hspace{0.2pt}\texttt {D}}& - F_{{\hspace{0.2pt}\texttt {I}}{\hspace{0.2pt}\texttt {c}}}^Ḭ(t) \\ \cline {2-4}f_{\hspace{0.2pt}\texttt {l}}&&\frac{\partial \mathcal {H}}{\partial \varphi _{\hspace{0.2pt}\texttt {l}}}(\varphi _{\hspace{0.2pt}\texttt {l}},\varphi _{\hspace{0.2pt}\texttt {L}}) {-} F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {l}}}^{\mathcal {H}}{\varphi _{\hspace{0.2pt}\texttt {L}}}(\varphi _{\hspace{0.2pt}\texttt {l}},\varphi _{\hspace{0.2pt}\texttt {L}}) &&- F_{{\hspace{0.2pt}\texttt {I}}{\hspace{0.2pt}\texttt {l}}}^Ḭ(t) \\f_{\hspace{0.2pt}\texttt {L}}& F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {c}}}\frac{\partial \mathcal {H}}{\partial q_{\hspace{0.2pt}\texttt {c}}}(q_{\hspace{0.2pt}\texttt {c}},q_{\hspace{0.2pt}\texttt {C}}) &+ F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {l}}} \dot{\varphi }_{\hspace{0.2pt}\texttt {l}}+ \dot{\varphi }_{\hspace{0.2pt}\texttt {L}}&+ F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {d}}} \nu _{\hspace{0.2pt}\texttt {d}}&+ F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {v}}} V(t) \\\cline {2-4}f_{\hspace{0.2pt}\texttt {d}}&&- F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {r}}}^{\mathcal {H}}{\varphi _{\hspace{0.2pt}\texttt {L}}}(\varphi _{\hspace{0.2pt}\texttt {l}},\varphi _{\hspace{0.2pt}\texttt {L}}) &+ \iota _{\hspace{0.2pt}\texttt {d}}{-} F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {d}}}^_{\hspace{0.2pt}\texttt {D}}&- F_{{\hspace{0.2pt}\texttt {I}}{\hspace{0.2pt}\texttt {d}}}^Ḭ(t) \\f_{\hspace{0.2pt}\texttt {D}}& F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {c}}} \frac{\partial \mathcal {H}}{\partial q_{\hspace{0.2pt}\texttt {c}}}(q_{\hspace{0.2pt}\texttt {c}},q_{\hspace{0.2pt}\texttt {C}}) &&+\nu _{\hspace{0.2pt}\texttt {D}}{+} F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {d}}} \nu _{\hspace{0.2pt}\texttt {d}}&+ F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {v}}} V(t) \\ \cline {2-4}\end{array}\) 


Here, we use that by (REF ) \(\iota _{\hspace{0.2pt}\texttt {d}}, \iota _{\hspace{0.2pt}\texttt {D}}, \nu _{\hspace{0.2pt}\texttt {d}}, \nu _{\hspace{0.2pt}\texttt {D}}\)  are explicit functions of the \(\widehat{x}=\widehat{x}^\mathcal {R}\)  variables.
This gives the schematic form (REF ) of the signature matrix:
\(\Sigma \ &=\ \begin{blockarray}{rccccclc}& q_{\hspace{0.2pt}\texttt {C}}& q_{\hspace{0.2pt}\texttt {c}}& \varphi _{\hspace{0.2pt}\texttt {l}}& \varphi _{\hspace{0.2pt}\texttt {L}}& \widehat{x} & \mbox{\footnotesize $c_i$} & \\ \begin{block}{r[cc|cc|c]lc}f_{\hspace{0.2pt}\texttt {C}}&\colorbox {lightgray!110}{\le 0} &\colorbox {lightgray!110}{\le 0} & - & - & - & \mbox{\footnotesize $1$} & n_{\hspace{0.2pt}\texttt {C}}\\f_{\hspace{0.2pt}\texttt {c}}&\colorbox {lightgray!110}{\le 1} &\colorbox {lightgray!110}{\widehat{I}{+}1}&{\le 0} &{\le 0} & {\le 0} & \mbox{\footnotesize $0$} & n_{\hspace{0.2pt}\texttt {c}}\\ \cline {2-6}f_{\hspace{0.2pt}\texttt {l}}& - & - &\colorbox {lightgray!110}{\le 0} &\colorbox {lightgray!110}{\le 0} & - & \mbox{\footnotesize $1$} & n_{\hspace{0.2pt}\texttt {l}}\\f_{\hspace{0.2pt}\texttt {L}}&\le 0 &{\le 0} &\colorbox {lightgray!110}{\le 1} &\colorbox {lightgray!110}{\widehat{I}{+}1}& {\le 0} & \mbox{\footnotesize $0$} & n_{\hspace{0.2pt}\texttt {L}}\\ \cline {2-6}f_{\hspace{0.2pt}\texttt {d}}& - & - &\le 0 &\le 0 &\colorbox {lightgray!110}{\le 0} & {\mbox{\footnotesize $1$}} & n_{\hspace{0.2pt}\texttt {d}}\\f_{\hspace{0.2pt}\texttt {D}}&\le 0 &\le 0 & - & - &\colorbox {lightgray!110}{\le 0} & {\mbox{\footnotesize $1$}} & n_{\hspace{0.2pt}\texttt {D}}\\ \end{block}\mbox{\footnotesize $d_j$}& \mbox{\footnotesize $1$} & \mbox{\footnotesize $1$} & \mbox{\footnotesize $1$} & \mbox{\footnotesize $1$} & {\mbox{\footnotesize $1$}} & & \\& n_{\hspace{0.2pt}\texttt {C}}& n_{\hspace{0.2pt}\texttt {c}}& n_{\hspace{0.2pt}\texttt {l}}& n_{\hspace{0.2pt}\texttt {L}}& n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}& &\end{blockarray}\) 


As in (REF ), provisional offsets are added, and the blocks with light-grey background are those that cannot contribute nonzeros to the associated system Jacobian.
The \(\mathbf {J}_{\mathcal {C}}\) , \(\mathbf {J}_{\cal L}\)  blocks are the same as before.
We need to show that the Jacobian
\(\mathbf {J}_{\mathcal {D}}=\begin{bmatrix}\partial f_{\hspace{0.2pt}\texttt {d}}/\partial \widehat{x}\\\partial f_{\hspace{0.2pt}\texttt {D}}/\partial \widehat{x}\end{bmatrix}\) 

is nonsingular.
Let \(P_1\)  be the permutation matrix that swaps \(\iota _{\hspace{0.2pt}\texttt {D}}\)  and \(\nu _{\hspace{0.2pt}\texttt {D}}\)  in \([\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}};\nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}]=[\iota _{\hspace{0.2pt}\texttt {d}};\iota _{\hspace{0.2pt}\texttt {D}};\nu _{\hspace{0.2pt}\texttt {d}};\nu _{\hspace{0.2pt}\texttt {D}}]\) , and let \(P=P_1P_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {G}}}\)  where \(P_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {G}}}\)  is as in (REF ). Denote \(S=\left[{\begin{matrix}0& -F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {d}}}^ F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {d}}} & 0\end{matrix}}\right]\) .
The part of the \(f_{\hspace{0.2pt}\texttt {d}}\) , \(f_{\hspace{0.2pt}\texttt {D}}\)  equations that depends on \(\widehat{x}\)  is
\(\begin{bmatrix}\iota _{\hspace{0.2pt}\texttt {d}}\\ \nu _{\hspace{0.2pt}\texttt {D}}\end{bmatrix} + \begin{bmatrix}0& -F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {d}}}^ F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {d}}} & 0\end{bmatrix}\begin{bmatrix}\nu _{\hspace{0.2pt}\texttt {d}}\\ \iota _{\hspace{0.2pt}\texttt {D}}\end{bmatrix}=[I,\;S]P_1\, \begin{bmatrix}\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\\ \nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\end{bmatrix}= [I,\;S]\, P \begin{bmatrix}\widehat{x}\\ \rho (\widehat{x})\end{bmatrix}.\) 

Denote \(\rho ^{\prime }=\partial \rho /\partial \widehat{x}\) . Then
\(\mathbf {J}_{\mathcal {D}}&= \frac{\partial }{\partial \widehat{x}} \begin{bmatrix}f_{\hspace{0.2pt}\texttt {d}}\\f_{\hspace{0.2pt}\texttt {D}}\end{bmatrix}= [I,\;S]\; P\; \begin{bmatrix}I\\ \rho ^{\prime }\end{bmatrix},\) 

where \(\rho ^{\prime }\)  is PD by assumption.
By Lemma REF  \((I,\rho ^{\prime ) is a NDMP.Let [C_1, C_2] = [I,\rho ^{\prime ]P. By Lemma~\ref {lm:PDMP1}, (C_1, C_2) is a NDMP, C_1 and C_2 are nonsingular, andC^{-1}_2 C_1 is PD.Then\mathbf {J}_{\mathcal {D}} = [I, S]\begin{bmatrix}C^1\\C_2^{bmatrix}= C_1^ SC_2^ (C_1 - C_2S)^ [C_2(C^{-1}_2 C_1 - S)]^Since S is skew-symmetric the nonsingularity of \mathbf {J}_{\mathcal {D}} follows fromLemma~\ref {lm:PDone}.\end{bmatrix}}In a similar way as before we get the following results.\begin{theorem}Assume that {\bf M2-1} to {\bf M2-3} hold. If n_{\hspace{0.2pt}\texttt {C}}=n_{\hspace{0.2pt}\texttt {l}}=n_{\hspace{0.2pt}\texttt {d}}=n_{\hspace{0.2pt}\texttt {D}}=0, thesystem (\ref {eq:CpH_Model2}) is an ODE.Otherwise, (\ref {eq:CpH_Model2}) is a DAE of index 1.\end{theorem}}{\begin{examp}\small \sf The running example can be described in the form of Model 2.Below is shown the signature matrix, alongside the actual equations given by (\ref {eq:ex4_1}) with (\ref {eq:ex4_2}) inserted.{{\begin{@align}{1}{-1}\Sigma &=\raisebox {-1.2ex}{\begin{blockarray}{rccccccl}& q_{{\hspace{0.2pt}\texttt {C}}_2} & q_{{\hspace{0.2pt}\texttt {C}}_1} & \varphi _{{\hspace{0.2pt}\texttt {L}}_1} & \varphi _{{\hspace{0.2pt}\texttt {L}}_2} & \iota _{{\hspace{0.2pt}\texttt {R}}} & \nu _{{\hspace{0.2pt}\texttt {G}}} &\mbox{\footnotesize $c_i$} \\\begin{block}{c[cc|cc|cc]l}f_{{\hspace{0.2pt}\texttt {C}}_2} &\colorbox {lightgray}{ 0 } &\colorbox {lightgray}{ 0 } & - & - & - & - &\mbox{\footnotesize $1$} \\f_{{\hspace{0.2pt}\texttt {C}}_1} &\colorbox {lightgray}{ 1 } &\colorbox {lightgray}{ 1 } & - & 0 & - & - &\mbox{\footnotesize $0$} \\ \cline {2-7}f_{{\hspace{0.2pt}\texttt {L}}_1} & - & - &\colorbox {lightgray}{ 0 } &\colorbox {lightgray}{ 0 } & - & - &\mbox{\footnotesize $1$} \\f_{{\hspace{0.2pt}\texttt {L}}_2} & - & 0 &\colorbox {lightgray}{ 1 } &\colorbox {lightgray}{ 1 } & 0 & - &\mbox{\footnotesize $0$} \\ \cline {2-7}f_{{\hspace{0.2pt}\texttt {R}}} & - & - & - & - &\colorbox {lightgray}{ 0 } &\colorbox {lightgray}{ 0 } &\mbox{\footnotesize $1$} \\f_{{\hspace{0.2pt}\texttt {G}}} & - & - & - & 0 &\colorbox {lightgray}{ 0 } &\colorbox {lightgray}{ 0 } &\mbox{\footnotesize $1$} \\\end{block}\mbox{\footnotesize $d_j$}& \mbox{\footnotesize $1$} & \mbox{\footnotesize $1$} & \mbox{\footnotesize $1$} & \mbox{\footnotesize $1$} & \mbox{\footnotesize $1$} & \mbox{\footnotesize $1$} & \end{blockarray}},\quad {\left\lbrace \begin{array}{r@{\;=\;}l@{}c}0& q_{{\hspace{0.2pt}\texttt {C}}_2}/C_2-q_{{\hspace{0.2pt}\texttt {C}}_1}/C_1 & +V(t) \\0& \dot{q}_{{\hspace{0.2pt}\texttt {C}}_1} +\dot{q}_{{\hspace{0.2pt}\texttt {C}}_2}-\varphi _{{\hspace{0.2pt}\texttt {L}}_2}/L_2 & + I(t)\\0& \varphi _{{\hspace{0.2pt}\texttt {L}}_1}/L_1 -\varphi _{{\hspace{0.2pt}\texttt {L}}_2}/L_2 & +I(t) \\0& \dot{\varphi }_{{\hspace{0.2pt}\texttt {L}}_2} +q_{{\hspace{0.2pt}\texttt {C}}_1}/C_1+R\iota _{{\hspace{0.2pt}\texttt {R}}} +\dot{\varphi }_{{\hspace{0.2pt}\texttt {L}}_1} &-V(t) \\0& \nu _{{\hspace{0.2pt}\texttt {G}}} +R\iota _{{\hspace{0.2pt}\texttt {R}}} & -V(t) \\0& \iota _{{\hspace{0.2pt}\texttt {R}}}- G\nu _{{\hspace{0.2pt}\texttt {G}}}-\varphi _{{\hspace{0.2pt}\texttt {L}}_2}/L_2 &\end{array}\right.}\end{@align}}}From Figure~\ref {fg:CpH8by5ex_graph0} one can verify each equation is KCL applied to a cutset, or KVL applied to a cycle, of tree \texttt {T}.\end{examp}}\)
CpH Circuit Model 3
For the CpH Circuit Model 3 we make the following assumptions.
[M3-1.]
The assumptions (A1)–(A6) hold.

The behavior of the external ports \(\mathcal {P}\)  is described by
an implicit relation \(p(\iota _\mathcal {P},\nu _\mathcal {P})=0\) .

The dissipative structure is described by a user-chosen mixed form (REF ) (cf. §REF ).

We restrict to \(\mathcal {P}={\hspace{0.2pt}{\hspace{0.2pt}\texttt {V}}}\stackrel{\cdot }{\smash{\cup }}{\hspace{0.2pt}{\hspace{0.2pt}\texttt {I}}}\)  and
further split these sets into controlled and independent sources as
\({\hspace{0.2pt}{\hspace{0.2pt}\texttt {V}}}={\hspace{0.2pt}{\hspace{0.2pt}\texttt {V}}}^\text{c}\stackrel{\cdot }{\smash{\cup }}{\hspace{0.2pt}{\hspace{0.2pt}\texttt {V}}}^\text{i},\quad {\hspace{0.2pt}{\hspace{0.2pt}\texttt {I}}}={\hspace{0.2pt}{\hspace{0.2pt}\texttt {I}}}^\text{c}\stackrel{\cdot }{\smash{\cup }}{\hspace{0.2pt}{\hspace{0.2pt}\texttt {I}}}^\text{i}.\) 
Independent sources are given by explicit source functions \(\iota _{{\hspace{0.2pt}\texttt {I}}^\text{i}}=I(t)\)  and \(\nu _{{\hspace{0.2pt}\texttt {v}}^\text{i}}=V(t)\)  and can again be seen as input \(u=\begin{bmatrix}\nu _{{\hspace{0.2pt}\texttt {v}}^\text{i}}\\\iota _{{\hspace{0.2pt}\texttt {I}}^\text{i}}\end{bmatrix}\) , with corresponding output \(y=\begin{bmatrix}\iota _{{\hspace{0.2pt}\texttt {v}}^\text{i}}\\\nu _{{\hspace{0.2pt}\texttt {I}}^\text{i}}\end{bmatrix}\) .
Moreover, we pose the following assumptions for
controlled voltage source (CVS) and controlled current source (CCS):
[VS-1.]

No CVS is part of a \({\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}{\hspace{0.2pt}{\hspace{0.2pt}\texttt {V}}}\) -cycle.

Controlling voltages of CVS are voltages of capacitances or independent voltage sources.

Controlling currents of CVS are currents of inductances or independent current sources.

[CS-1.]

No CCS is part of an \({\hspace{0.2pt}{\hspace{0.2pt}\texttt {L}}}{\hspace{0.2pt}{\hspace{0.2pt}\texttt {I}}}\) -cutset.

Controlling voltages of CCS are voltages of capacitances, independent voltage sources or resistances.

Controlling currents of CCS are currents of inductances, independent current sources or resistances.

They are adapted from the conditions for controlled sources given in [2]}.
Under these assumptions and
using the constitutive relations \(\nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}= \partial \mathcal {H}/\partial q\) , \(\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {L}}}= \partial \mathcal {H}/\partial \varphi \) , the source functions of independent sources, as well as
(REF ) according to M3-3, the controlled sources are described by
\(\nu _{{\hspace{0.2pt}\texttt {v}}^\text{c}}&=\gamma _{{\hspace{0.2pt}\texttt {v}}}\left(q_{\hspace{0.2pt}\texttt {C}},q_{\hspace{0.2pt}\texttt {c}},\varphi _{\hspace{0.2pt}\texttt {L}},\varphi _{\hspace{0.2pt}\texttt {l}},t\right),\\\iota _{{\hspace{0.2pt}\texttt {I}}^\text{c}}&=\gamma _{{\hspace{0.2pt}\texttt {I}}}\left(q_{\hspace{0.2pt}\texttt {C}},q_{\hspace{0.2pt}\texttt {c}},\varphi _{\hspace{0.2pt}\texttt {L}},\varphi _{\hspace{0.2pt}\texttt {l}},\nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {G}}},\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {R}}},t\right).\) 

Casting the DAE to input-output form as in §REF  and restricting
to the state equation given by (REF ) with a given input
\(u=\begin{bmatrix}\nu _{{\hspace{0.2pt}\texttt {v}}^\text{i}}\\\iota _{{\hspace{0.2pt}\texttt {I}}^\text{i}}\end{bmatrix}\)  leads to the CpH DAE
\(0=f(t,x,\dot{x},u)\) 

of the form
\(0= &\begin{bmatrix} f_{{\hspace{0.2pt}\texttt {v}}^\text{c}}\\f_{\hspace{0.2pt}\texttt {c}}\\ f_{\hspace{0.2pt}\texttt {d}}\\ f_{\hspace{0.2pt}\texttt {l}}\end{bmatrix} = \begin{bmatrix}\iota _{{\hspace{0.2pt}\texttt {v}}^\text{c}}\\\dot{q}_{\hspace{0.2pt}\texttt {c}}\\ \iota _{\hspace{0.2pt}\texttt {d}}\\ \partial \mathcal {H}/\partial \varphi _{\hspace{0.2pt}\texttt {l}}\end{bmatrix}- \begin{bmatrix}0 &F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {v}}^\text{c}}^F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {v}}^\text{c}}^F_{{\hspace{0.2pt}\texttt {I}}^\text{c} {\hspace{0.2pt}\texttt {v}}^\text{c}}^F_{{\hspace{0.2pt}\texttt {C}}{\hspace{0.2pt}\texttt {c}}}^F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {c}}}^F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {c}}}^F_{{\hspace{0.2pt}\texttt {I}}^\text{c} {\hspace{0.2pt}\texttt {c}}}^0&F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {d}}}^F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {d}}}^F_{{\hspace{0.2pt}\texttt {I}}^\text{c} {\hspace{0.2pt}\texttt {d}}}^0&0&F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {l}}}^0 \end{bmatrix}\begin{bmatrix}\dot{q}_{\hspace{0.2pt}\texttt {C}}\\ \iota _{\hspace{0.2pt}\texttt {D}}\\ \partial \mathcal {H}/\partial \varphi _{\hspace{0.2pt}\texttt {L}}\\ \iota _{{\hspace{0.2pt}\texttt {I}}^\text{c}}\end{bmatrix}- \begin{bmatrix}F_{{\hspace{0.2pt}\texttt {I}}^\text{i} {\hspace{0.2pt}\texttt {v}}^\text{c}}^F_{{\hspace{0.2pt}\texttt {I}}^\text{i} {\hspace{0.2pt}\texttt {c}}}^ F_{{\hspace{0.2pt}\texttt {I}}^\text{i} {\hspace{0.2pt}\texttt {d}}}^F_{{\hspace{0.2pt}\texttt {I}}^\text{i} {\hspace{0.2pt}\texttt {l}}}^{bmatrix}I(t),\\[.5ex]0= &\begin{bmatrix}f_{\hspace{0.2pt}\texttt {C}}\\ f_{\hspace{0.2pt}\texttt {D}}\\f_{\hspace{0.2pt}\texttt {L}}\\ f_{{\hspace{0.2pt}\texttt {I}}^\text{c}}\end{bmatrix} = \begin{bmatrix}\partial \mathcal {H}/\partial q_{\hspace{0.2pt}\texttt {C}}\\ \nu _{\hspace{0.2pt}\texttt {D}}\\ \dot{\varphi }_{\hspace{0.2pt}\texttt {L}}\\ \nu _{{\hspace{0.2pt}\texttt {I}}^\text{c}}\end{bmatrix}+ \begin{bmatrix}0 &F_{{\hspace{0.2pt}\texttt {C}}{\hspace{0.2pt}\texttt {c}}} &0 &0 \\F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {v}}^\text{c}} &F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {c}}} &F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {d}}} &0 \\F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {v}}^\text{c}} &F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {c}}} &F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {d}}} &F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {l}}}\\F_{{\hspace{0.2pt}\texttt {I}}^\text{c} {\hspace{0.2pt}\texttt {v}}^\text{c}} &F_{{\hspace{0.2pt}\texttt {I}}^\text{c} {\hspace{0.2pt}\texttt {c}}} &F_{{\hspace{0.2pt}\texttt {I}}^\text{c} {\hspace{0.2pt}\texttt {d}}} &0\end{bmatrix}\begin{bmatrix}\nu _{{\hspace{0.2pt}\texttt {v}}^\text{c}}\\ \partial \mathcal {H}/\partial q_{\hspace{0.2pt}\texttt {c}}\\ \nu _{\hspace{0.2pt}\texttt {d}}\\ \dot{\varphi }_{\hspace{0.2pt}\texttt {l}}\end{bmatrix}+\begin{bmatrix}F_{{\hspace{0.2pt}\texttt {C}}{\hspace{0.2pt}\texttt {v}}^\text{i}}\\F_{{\hspace{0.2pt}\texttt {D}}{\hspace{0.2pt}\texttt {v}}^\text{i}}\\F_{{\hspace{0.2pt}\texttt {L}}{\hspace{0.2pt}\texttt {v}}^\text{i}}\\ F_{{\hspace{0.2pt}\texttt {I}}^\text{c} {\hspace{0.2pt}\texttt {v}}^\text{i}}\end{bmatrix}V(t),\\0=& \begin{bmatrix}p_{{\hspace{0.2pt}\texttt {I}}^\text{c}}\\p_{{\hspace{0.2pt}\texttt {v}}^\text{c}}\end{bmatrix} = \begin{bmatrix}\iota _{{\hspace{0.2pt}\texttt {I}}^\text{c}}\\\nu _{{\hspace{0.2pt}\texttt {v}}^\text{c}}\end{bmatrix}-\begin{bmatrix}\gamma _{{\hspace{0.2pt}\texttt {I}}}\left(q_{\hspace{0.2pt}\texttt {C}},q_{\hspace{0.2pt}\texttt {c}},\varphi _{\hspace{0.2pt}\texttt {L}},\varphi _{\hspace{0.2pt}\texttt {l}},\nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {G}}},\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {R}}},t\right)\\\gamma _{{\hspace{0.2pt}\texttt {v}}}\left(q_{\hspace{0.2pt}\texttt {C}},q_{\hspace{0.2pt}\texttt {c}},\varphi _{\hspace{0.2pt}\texttt {L}},\varphi _{\hspace{0.2pt}\texttt {l}},t\right)\end{bmatrix},\end{bmatrix}where, similar as before we use that \iota _{\hspace{0.2pt}\texttt {d}}, \iota _{\hspace{0.2pt}\texttt {D}}, \nu _{\hspace{0.2pt}\texttt {d}}, \nu _{\hspace{0.2pt}\texttt {D}} are explicit functions of \nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {G}}} and \iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {R}}}.The assumptions {\bf VS-1} and {\bf CS-1} imply that F_{{\hspace{0.2pt}\texttt {C}}{\hspace{0.2pt}\texttt {v}}^\text{i}}=0 and F_{{\hspace{0.2pt}\texttt {I}}^\text{i}{\hspace{0.2pt}\texttt {l}}}=0.In (\ref {eq:CpH_Model3}) the vector x consists of the energy state variables z=(q,\varphi ) plus the \widehat{x} components \nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {G}}}, \iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {R}}},\nu _{{\hspace{0.2pt}\texttt {v}}^\text{c}}, \iota _{{\hspace{0.2pt}\texttt {I}}^\text{c}}, \iota _{{\hspace{0.2pt}\texttt {v}}^\text{c}}, and \nu _{{\hspace{0.2pt}\texttt {I}}^\text{c}} resulting in one DAE variable for each capacitor, inductor and dissipator andtwo DAE variables for each controlled source in the circuit, for a total DAE size of N = n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}+n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {L}}}+ n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}+ 2(n_{{\hspace{0.2pt}\texttt {v}}^\text{c}}+n_{{\hspace{0.2pt}\texttt {I}}^\text{c}}).\) 

Theorem 4.8 If the assumptions M3-1–M3-3, VS-1–VS-3, and CS-1–CS-3 hold,
and in addition \(F_{{\hspace{0.2pt}\texttt {I}}^\text{c}{\hspace{0.2pt}\texttt {d}}}=0\)  (i.e., there are no cycles containing twig resistors and controlled current sources),
then the CpH DAE (REF ) is SA-amenable.

For Model 3 the signature matrix takes a schematic form as in Figure REF , where \(\Sigma \)  has been ordered such that the corresponding Jacobian is in block triangular form.
<FIGURE>Again, a transversal can be found in the diagonal shaded blocks.
Blocks with light-grey background cannot contribute nonzero entries to the system Jacobian \(\mathbf {J}\) .
Thus, \(\mathbf {J}\)  takes the form
\(\mathbf {J} \ &=\ \begin{blockarray}{rccccccl}& q_{\hspace{0.2pt}\texttt {C}},q_{\hspace{0.2pt}\texttt {c}}& \varphi _{\hspace{0.2pt}\texttt {l}},\varphi _{\hspace{0.2pt}\texttt {L}}& \nu _{{\hspace{0.2pt}\texttt {v}}^\text{c}}&\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {R}}},\nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {G}}}&\iota _{{\hspace{0.2pt}\texttt {I}}^\text{c}}&\iota _{{\hspace{0.2pt}\texttt {v}}^\text{c}},\nu _{{\hspace{0.2pt}\texttt {I}}^\text{c}} \\\begin{block}{r[c|c|c|c|c|c]l}f_{\hspace{0.2pt}\texttt {C}},f_{\hspace{0.2pt}\texttt {c}}& \colorbox {lightgray}{\mathbf {J}_{\mathcal {C}}} & 0 &0&0&0 & 0 &\mbox{\footnotesize $n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}=n_{\hspace{0.2pt}\texttt {C}}{+}n_{\hspace{0.2pt}\texttt {c}}$} \\ \cline {2-7}f_{\hspace{0.2pt}\texttt {l}},f_{\hspace{0.2pt}\texttt {L}}& 0 & \colorbox {lightgray}{\mathbf {J}_{\mathcal {L}}} & 0 &0&0&0 &\mbox{\footnotesize $n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {L}}}=n_{\hspace{0.2pt}\texttt {l}}{+}n_{\hspace{0.2pt}\texttt {L}}$} \\\cline {2-7}p_{{\hspace{0.2pt}\texttt {v}}^\text{c}}& \ast &\ast &\colorbox {lightgray}{I}&0&0&0&\mbox{\footnotesize $n_{{\hspace{0.2pt}\texttt {v}}^\text{c}}$} \\\cline {2-7}f_{\hspace{0.2pt}\texttt {d}},f_{\hspace{0.2pt}\texttt {D}}& \ast & \ast & \ast &\colorbox {lightgray}{\mathbf {J}_{\mathcal {D}}}&0&0 &\mbox{\footnotesize $n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}=n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {R}}}{+}n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {G}}}$} \\\cline {2-7}p_{{\hspace{0.2pt}\texttt {I}}^\text{c}}& \ast &\ast &0&\ast &\colorbox {lightgray}{I}&0&\mbox{\footnotesize $n_{{\hspace{0.2pt}\texttt {I}}^\text{c}}$}\\\cline {2-7}f_{{\hspace{0.2pt}\texttt {v}}^\text{c}},f_{{\hspace{0.2pt}\texttt {I}}^\text{c}}& \ast & \ast &\ast &\ast &\ast &\colorbox {lightgray}{I} &\mbox{\footnotesize $n_{{\hspace{0.2pt}\texttt {v}}^\text{c}}{+}n_{{\hspace{0.2pt}\texttt {I}}^\text{c}}$} \\\end{block}& \mbox{\footnotesize $n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}$} & \mbox{\footnotesize $n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {L}}}$} &\mbox{\footnotesize $n_{{\hspace{0.2pt}\texttt {v}}^\text{c}}$}& \mbox{\footnotesize $n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}$} &\mbox{\footnotesize $n_{{\hspace{0.2pt}\texttt {I}}^\text{c}}$}&\mbox{\footnotesize $n_{{\hspace{0.2pt}\texttt {v}}^\text{c}}$}+\mbox{\footnotesize $n_{{\hspace{0.2pt}\texttt {I}}^\text{c}}$}\end{blockarray}\) 


The \(\ast \)  denote some (zero or nonzero) blocks that play no role for the further analysis.
Under the given assumptions the diagonal blocks, and thus \(\mathbf {J}\) , are nonsingular as discussed before.
Theorem 4.9 Assume that M3-1–M3-3, VS-1–VS-3, and CS-1–CS-3 hold, and \(F_{{\hspace{0.2pt}\texttt {I}}^\text{c}{\hspace{0.2pt}\texttt {d}}}=0\) .
If \(n_{\hspace{0.2pt}\texttt {C}}=n_{\hspace{0.2pt}\texttt {l}}=n_{\hspace{0.2pt}\texttt {d}}=n_{\hspace{0.2pt}\texttt {D}}=n_{{\hspace{0.2pt}\texttt {v}}^\text{c}}=n_{{\hspace{0.2pt}\texttt {I}}^\text{c}}=0\)  the system (REF ) is an ODE,
otherwise a DAE of index 1.

Along the same lines as the proof of Theorem REF .
Reduction to ODE
An advantage of knowing that the DAE is SA-amenable is that the \(\Sigma \) -method's standard solution scheme (see (REF )) tells us which equations need to be differentiated in order to find a numerical solution.
This can be used in several ways, e.g. to construct a scheme for reducing the DAE to an explicit ODE.
It is explained in [3]} why this might be desirable, e.g. for storage reasons, when the number of DOF is small compared to the total number of variables.
We show this process for a Model 2 formulation; with small changes it also works for the other model formulations.

Algorithm 4.1 (Reduction of Model 2 form to explicit ODE)  

Consider the Model 2 equations (REF ).
To simplify notation we do not distinguish between column and row vectors.
[]
I. Corresponding to stage \(k{=}-1\)  of the \(\Sigma \) -method:
[I.1.]

As \(\mathbf {J}_\mathcal {C}\)  is nonsingular, its \(f_{\hspace{0.2pt}\texttt {C}}\)  rows \(\partial f_{\hspace{0.2pt}\texttt {C}}/\partial q\)  (with \(q=(q_{\hspace{0.2pt}\texttt {C}},q_{\hspace{0.2pt}\texttt {c}})\) ) have full row rank. Choose columns to form a nonsingular \(n_{\hspace{0.2pt}\texttt {C}}\times n_{\hspace{0.2pt}\texttt {C}}\)  matrix.
Split \(q\)  as \((\overline{q},\widehat{q})\)  with \(\widehat{q}\)  belonging to the chosen columns.

By the IFT we can solve \(f_{\hspace{0.2pt}\texttt {C}}=0\)  for \(\widehat{q}\)  in terms of \(\overline{q}\)  and the right-hand-side item \(V(t)\) .
Let \(\mathop {\textrm {\sl fun}}\)  denote some function that does not need a specific name, then:
\( \widehat{q} = \mathop {\textrm {\sl fun}}(\overline{q},V(t)) \) 

Similarly choose columns of \(\partial f_{\hspace{0.2pt}\texttt {l}}/\partial \varphi \)  (with \(\varphi =(\varphi _{\hspace{0.2pt}\texttt {l}},\varphi _{\hspace{0.2pt}\texttt {L}})\) ) to form a nonsingular \(n_{\hspace{0.2pt}\texttt {l}}\times n_{\hspace{0.2pt}\texttt {l}}\)  matrix.
Split \(\varphi \)  as \((\overline{\varphi },\widehat{\varphi })\)  with \(\widehat{\varphi }\)  belonging to the chosen columns.

By the IFT we can solve \(f_{\hspace{0.2pt}\texttt {l}}=0\)  for \(\widehat{\varphi }\)  in terms of \(\overline{\varphi }\)  and the right-hand-side item \(I(t)\) :
\( \widehat{\varphi }= \mathop {\textrm {\sl fun}}(\overline{\varphi },I(t)). \) 


Recall \(\iota _{\hspace{0.2pt}\texttt {d}}, \iota _{\hspace{0.2pt}\texttt {D}}, \nu _{\hspace{0.2pt}\texttt {d}}, \nu _{\hspace{0.2pt}\texttt {D}}\)  are explicit functions of \(\widehat{x}\) .
As \(\mathbf {J}_\mathcal {D}= \partial (f_{\hspace{0.2pt}\texttt {d}},f_{\hspace{0.2pt}\texttt {D}})/\partial \widehat{x}\)  is nonsingular, by the IFT we can solve \((f_{\hspace{0.2pt}\texttt {d}},f_{\hspace{0.2pt}\texttt {D}})=0\)  for \(\widehat{x}\)  in terms of \(V(t)\) , \(I(t)\) , as well as of \(q_{\hspace{0.2pt}\texttt {C}},q_{\hspace{0.2pt}\texttt {c}},\varphi _{\hspace{0.2pt}\texttt {l}},\varphi _{\hspace{0.2pt}\texttt {L}}\)  which reordered are the same as \(\overline{q},\widehat{q},\overline{\varphi },\widehat{\varphi }\) .
That is,
\( \widehat{x} = \mathop {\textrm {\sl fun}}(\overline{q},\widehat{q},\overline{\varphi },\widehat{\varphi },V(t),I(t)). \) 


Now \(V(t),I(t)\)  comprise the control input \(u(t)\)  so write the results of steps I.2, I.4, I.5 as
\(\widehat{q} &= \mathop {\textrm {\sl fun}}(\overline{q},u(t)), &\widehat{\varphi }&= \mathop {\textrm {\sl fun}}(\overline{\varphi },u(t)),&\widehat{x} &= \mathop {\textrm {\sl fun}}(\overline{q},\widehat{q},\overline{\varphi },\widehat{\varphi },u(t)).\) 


Eliminate \(\widehat{q},\widehat{\varphi }\)  by substituting the first two equations of (REF ) into the third to get
\(\widehat{x} = \mathop {\textrm {\sl fun}}(\overline{q},\overline{\varphi },u(t)).\) 


II. Corresponding to stage \(k{=}0\)  of the \(\Sigma \) -method, so we differentiate \(f_{\hspace{0.2pt}\texttt {C}},f_{\hspace{0.2pt}\texttt {l}}\) :
[II.1.]

Solve \((\dot{f}_{\hspace{0.2pt}\texttt {C}},f_{\hspace{0.2pt}\texttt {c}})=0\)  as a nonsingular linear system with matrix \(\mathbf {J}_\mathcal {C}\)  for \(\dot{q}= (\dot{\overline{q}},\dot{\widehat{q}})\)  as a function of \((\overline{q},\widehat{q},\overline{\varphi },\widehat{\varphi },\widehat{x},\dot{V}(t),I(t))\) .
Discard \(\dot{\widehat{q}}\)  which will not be used.

Solve \((\dot{f}_{\hspace{0.2pt}\texttt {l}},f_{\hspace{0.2pt}\texttt {L}})=0\)  as a nonsingular linear system with matrix \(\mathbf {J}_\mathcal {L}\)  for \(\dot{\varphi }= (\dot{\overline{\varphi }},\dot{\widehat{\varphi }})\)  as a function of \((\overline{q},\widehat{q},\overline{\varphi },\widehat{\varphi },\widehat{x},V(t),\dot{I}(t))\) .
Discard \(\dot{\widehat{\varphi }}\)  which will not be used.

Regard \(V(t),I(t)\)  as comprising \(u(t)\)  as above.
Eliminate \(\widehat{q},\widehat{\varphi }\)  and \(\widehat{x}\)  using (REF , REF ) to get \(\dot{\overline{q}}\)  and \(\dot{\overline{\varphi }}\)  as functions of \((\overline{q},\overline{\varphi },u(t),\dot{u}(t))\)  yielding a first-order ODE of size \(n_{\hspace{0.2pt}\texttt {c}}{+}n_{\hspace{0.2pt}\texttt {L}}=\textsc {dof} \)
\(\left.\begin{aligned}\dot{\overline{q}} &= \mathop {\textrm {\sl fun}}(\overline{q},\overline{\varphi },u(t),\dot{u}(t)), \\\dot{\overline{\varphi }} &= \mathop {\textrm {\sl fun}}(\overline{q},\overline{\varphi },u(t),\dot{u}(t)),\end{aligned}\right\rbrace .\) 


Thus, by three (in general nonlinear) solves corresponding to an \(n_{\hspace{0.2pt}\texttt {C}}\times n_{\hspace{0.2pt}\texttt {C}}\)  part of the \(\mathcal {C}\)  block, an \(n_{\hspace{0.2pt}\texttt {l}}\times n_{\hspace{0.2pt}\texttt {l}}\)  part of the \(\mathcal {L}\)  block, and all of the \(\mathcal {D}\)  block,
followed by two linear solves corresponding to the whole \(\mathcal {C}\)  and \(\mathcal {L}\)  blocks, we convert the DAE to an ODE in \((\overline{q},\overline{\varphi })\) , with \(u(t), \dot{u}(t)\)  as control input.
The remaining DAE variables, which comprise \(\widehat{q},\widehat{\varphi }\)  and \(\widehat{x}\) , can be retrieved from a solution of the ODE by the purely algebraic equations (REF ).
Remark 4.10 Algorithm REF  is the \(\Sigma \) -method's version of the dummy derivatives process [4]}, specialised to the present case. The discarded variables \(\dot{\widehat{q}},\dot{\widehat{\varphi }}\)  correspond to the “dummy derivatives” in the sense of [4]}.
See [3]} for more details.
In the linear time-invariant case the columns to make nonsingular matrices can be chosen once for all.
But in general a matrix might become singular as solution proceeds, and one needs a new choice of columns (dummy pivoting).
A cheap way to find the needed \(n_{\hspace{0.2pt}\texttt {C}}\)  columns of the \(f_{\hspace{0.2pt}\texttt {C}}\)  rows, and \(n_{\hspace{0.2pt}\texttt {l}}\)  columns of the \(f_{\hspace{0.2pt}\texttt {l}}\)  rows, is the Scholz–Steinbrecher method [7]}: using some HVT of \(\Sigma \) , choose the columns where it intersects the \(f_{\hspace{0.2pt}\texttt {C}}\)  rows and \(f_{\hspace{0.2pt}\texttt {l}}\)  rows respectively.
In general, this gives nonsingular matrices with high probability, but is not guaranteed.
For CpH DAEs, it always provides a right choice if each inductor and each capacitor is independent. Namely, it chooses \(q_{\hspace{0.2pt}\texttt {C}}\)  and \(\varphi _{\hspace{0.2pt}\texttt {l}}\)  for \({\widehat{q}}\) , \({\widehat{\varphi }}\)  and the resulting nonsingular matrices are the inverses of the corresponding conductance or inductance matrix, respectively.
In code, packing up the repeated nonlinear and linear solves within an “[basicstyle=]odefun” routine to give to a standard ODE solver makes for opaque programming.
It is more convenient to use a “reverse communication” interface so that the solves can be done inline in the calling program.
The NAG Library in 2015 provided a reverse communication Runge–Kutta solver for just such applications, replaced since 2018 by a variable-order variable-step Adams code [basicstyle=]d02qgf [8]}.

Example 4 
The running example's \(\Sigma \) -matrix and equations are in ().
The steps of Algorithm REF  are:
I.1,I.3
Each circuit element is independent, so the Scholz–Steinbrecher method is guaranteed to succeed.
The main diagonal is an HVT. With this HVT, Scholz–Steinbrecher selects \(\widehat{q}=q_{{\hspace{0.2pt}\texttt {C}}_2}\) , \(\widehat{\varphi }=\varphi _{{\hspace{0.2pt}\texttt {L}}_1}\) .
The ODE variables will then be \((\overline{q},\overline{\varphi }) = (q_{{\hspace{0.2pt}\texttt {C}}_1},\varphi _{{\hspace{0.2pt}\texttt {L}}_2})\) .
I.2
Solve \(q_{{\hspace{0.2pt}\texttt {C}}_2} = \mathop {\textrm {\sl fun}}(q_{{\hspace{0.2pt}\texttt {C}}_1},V(t))\) .
I.4
Solve \(\varphi _{{\hspace{0.2pt}\texttt {L}}_1} = \mathop {\textrm {\sl fun}}(\varphi _{{\hspace{0.2pt}\texttt {L}}_2},I(t))\) .
I.5
Solve \((\iota _{{\hspace{0.2pt}\texttt {R}}},\nu _{{\hspace{0.2pt}\texttt {G}}}) = \widehat{x} = \mathop {\textrm {\sl fun}}(\varphi _{{\hspace{0.2pt}\texttt {L}}_2},V(t))\) .
In fact \(\nu _{{\hspace{0.2pt}\texttt {G}}}\)  is not used.
I.6,I.7
\(\widehat{x}\)  is already in the form \(\widehat{x} = \mathop {\textrm {\sl fun}}(\overline{q},\overline{\varphi },u(t))\) .
II.1
Solve \((\dot{f}_{\hspace{0.2pt}\texttt {C}},f_{\hspace{0.2pt}\texttt {c}})=0\)  for \((\dot{q}_{{\hspace{0.2pt}\texttt {C}}_2},\dot{q}_{{\hspace{0.2pt}\texttt {C}}_1}) = \mathop {\textrm {\sl fun}}(\varphi _{{\hspace{0.2pt}\texttt {L}}_2},\dot{V}(t),I(t))\) . Discard \(\dot{q}_{{\hspace{0.2pt}\texttt {C}}_2}\) .
II.2,II.3
Solve \((\dot{f}_{\hspace{0.2pt}\texttt {l}},f_{\hspace{0.2pt}\texttt {L}})=0\)  for \((\dot{\varphi }_{{\hspace{0.2pt}\texttt {L}}_1},\dot{\varphi }_{{\hspace{0.2pt}\texttt {L}}_2}) = \mathop {\textrm {\sl fun}}(q_{{\hspace{0.2pt}\texttt {C}}_1},\iota _{{\hspace{0.2pt}\texttt {R}}},\dot{I}(t),V(t))\) . Discard \(\dot{\varphi }_{{\hspace{0.2pt}\texttt {L}}_1}\) .
Eliminate \(\iota _{{\hspace{0.2pt}\texttt {R}}}\)  by I.5.

This results in an ODE
\((\dot{q}_{{\hspace{0.2pt}\texttt {C}}_1},\dot{\varphi }_{{\hspace{0.2pt}\texttt {L}}_2}) = \mathop {\textrm {\sl fun}}(q_{{\hspace{0.2pt}\texttt {C}}_1},\varphi _{{\hspace{0.2pt}\texttt {L}}_2}, u(t),\dot{u}(t))\) 
as follows.
\(\dot{q}_{{\hspace{0.2pt}\texttt {C}}_1} &= \left(\left(\frac{\varphi _{{\hspace{0.2pt}\texttt {L}}_2}}{L_2}-I(t)\right)/C_2 + \dot{V}(t)\right)\left\bad.\left(\frac{1}{C_1}+\frac{1}{C_2}\right)\right., \\\dot{\varphi }_{{\hspace{0.2pt}\texttt {L}}_2} &= \left(-\frac{q_{{\hspace{0.2pt}\texttt {C}}_1}}{C_1L_1}- \frac{R}{L_1}\left(\frac{GV(t) + \varphi _{{\hspace{0.2pt}\texttt {L}}_2}/L_2}{1+GR}\right) + \frac{V(t)}{L_1} + \dot{I}(t) \right)\left\bad.\left(\frac{1}{L_1}+\frac{1}{L_2}\right)\right..\) 

In this symbolic form it is somewhat elaborate. Program code to perform the algorithm steps numerically would be simpler.
Spectral analysis for linear time-invariant systems
When the non-source circuit elements are constant and linear, the system becomes a linear time-invariant (LTI) system
\(0 = Ax - B\dot{x}- u(t),\) 

with constant matrices \(A\)  and \(B\) .
One may view (REF ) as linearization of a nonlinear DAE system around an equilibrium point.
Eigenvalues and eigenvectors of electrical networks can be used to to examine stability of DC operating points; to determine the exponential stability of equilibrium points; or to study bifurcations [9]}, [10]}.
Several qualitative properties of equilibria can be characterized in terms of the spectrum of the corresponding matrix pencil \(A-\lambda B\) , see e.g. [11]}, [12]}, [13]}, [14]}.
The condition \(\Re \lambda < 0\)  guarantees asymptotic stability of the equilibrium in the sense of Lyapunov.
For several reasons, it can be of interest to ensure there exists a unique equilibrium, which is the case if and only if \(A\)  is invertible (or all eigenvalues are nonzero).
Non-trivial, purely imaginary eigenvalues are important in linear circuit applications since they characterize the existence of periodic solutions describing oscillations.
Also linear stability analysis of a power system via eigenvalues can be used to predict its degree of stability [15]}.
Sometimes a system starts out as passive but discretisation or model reduction may destroy passivity—whether this has happened can be observed by checking the eigenvalues.
As an example of finding basic LTI system data we show here how to compute eigenfunctions for a Model 2 formulation,
it works in the same way for Models 1 or 3.
One might adapt Algorithm REF —which in the LTI case yields a linear constant-coefficient ODE \(\dot{x}=Cx+g(t)\) —and compute eigenvalues and eigenvectors of the matrix \(C\) .
However, it is simpler to derive (REF ) directly and solve the corresponding generalised eigenvalue problem (GEVP).
In the Model 2 formulation, the DAE variables are \(x = (q,\varphi ,\widehat{x})\) , with
\(\widehat{x}=(\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {R}}},\nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {G}}})\)  according to the splitting of \({\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}={\hspace{0.2pt}{\hspace{0.2pt}\texttt {R}}}\stackrel{\cdot }{\smash{\cup }}{\hspace{0.2pt}{\hspace{0.2pt}\texttt {G}}}\)  into current-controlled (resistor) and voltage-controlled (conductor) edges.
The Hamiltonian Hessians in (A4) and (A5) are now constant and denoted by \(\mathbf {C}\in \mathbb {R}^{n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}\times n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}}\)  and \(\mathbf {L}\in \mathbb {R}^{n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {L}}}\times n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {L}}}}\) , respectively. Similarly, we denoted the constant Jacobian of \(\rho (\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {R}}},\nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {G}}})\)  in (REF ) by
\(R\in \mathbb {R}^{n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\times n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}}\) .
The relation (REF ) becomes
\(\begin{bmatrix}\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\\ \nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\end{bmatrix} &= P_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {G}}}\begin{bmatrix}\widehat{x}\\ R\widehat{x}\end{bmatrix} = P_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {G}}}\begin{bmatrix}I\\R\end{bmatrix}\widehat{x}= \begin{bmatrix}\mathbf {D}_\iota \\\mathbf {D}_\nu \end{bmatrix}\widehat{x}, \quad \text{say},\) 

with \(\mathbf {D}_\iota , \mathbf {D}_\nu \in \mathbb {R}^{n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\times n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}}\) , where \(P_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {G}}}\)  is the flip-matrix of the split \({\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}={\hspace{0.2pt}{\hspace{0.2pt}\texttt {R}}}\stackrel{\cdot }{\smash{\cup }}{\hspace{0.2pt}{\hspace{0.2pt}\texttt {G}}}\) .
Let the edge set \({\hspace{0.2pt}{\hspace{0.2pt}\texttt {X}}}= {\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}\stackrel{\cdot }{\smash{\cup }}{\hspace{0.2pt}{\hspace{0.2pt}\texttt {L}}}\stackrel{\cdot }{\smash{\cup }}{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\)  be split into \({\hspace{0.2pt}\texttt {x}}\stackrel{\cdot }{\smash{\cup }}{\hspace{0.2pt}\texttt {X}}\)  as in §REF .
The constitutive relations, being linear, give the current and voltage vectors on \({\hspace{0.2pt}{\hspace{0.2pt}\texttt {X}}}\)  as
\(\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {X}}}&=\begin{bmatrix}\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}\\ \iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {L}}}\\ \iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\end{bmatrix} =\begin{bmatrix}\dot{q}\\ \mathbf {L}\varphi \\ \mathbf {D}_\iota \widehat{x} \end{bmatrix}= \begin{bmatrix}0&&\\ &\mathbf {L}&\\ &&\mathbf {D}_\iota \end{bmatrix}x - \begin{bmatrix}-I&&\\ &0&\\ &&0\end{bmatrix}\dot{x}&= A_\iota x - B_\iota \dot{x},\text{ say},  \\\nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {X}}}&=\begin{bmatrix}\nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}\\ \nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {L}}}\\ \nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\end{bmatrix} =\begin{bmatrix} \mathbf {C}q\\ \dot{\varphi }\\ \mathbf {D}_\nu \widehat{x}\end{bmatrix}= \begin{bmatrix}\mathbf {C}&&\\ &0&\\ &&\mathbf {D}_\nu \end{bmatrix}x - \begin{bmatrix}0&&\\ &-I&\\ &&0\end{bmatrix}\dot{x}&= A_\nu x - B_\nu \dot{x},\text{ say}. \) 

Here \(A_\iota , B_\iota , A_\nu , B_\nu \)  are \(N\times N\)  matrices, where \(N=n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}+n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {L}}}+n_{\hspace{0.2pt}{\hspace{0.2pt}\texttt {D}}}\)  as in (REF ) is the DAE size.
Define the coordinate projection matrices \(\Pi _{\hspace{0.2pt}\texttt {x}}\in \mathbb {R}^{n_{\hspace{0.2pt}\texttt {x}}\times N}\) , \(\Pi _{\hspace{0.2pt}\texttt {X}}\in \mathbb {R}^{n_{\hspace{0.2pt}\texttt {X}}\times N}\)  that extract the \({\hspace{0.2pt}\texttt {x}}\)  and \({\hspace{0.2pt}\texttt {X}}\)  coordinates of a column \(N\) -vector.
They consist of the \({\hspace{0.2pt}\texttt {x}}\) -rows and \({\hspace{0.2pt}\texttt {X}}\) -rows respectively, of the \(N\times N\)  identity matrix. Then \(\iota _{\hspace{0.2pt}\texttt {x}}= \Pi _{\hspace{0.2pt}\texttt {x}}\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {X}}}\) , \(\iota _{\hspace{0.2pt}\texttt {X}}= \Pi _{\hspace{0.2pt}\texttt {X}}\iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {X}}}\)  and so on.
Thus, (REF ) may be written as
\(0 &=\begin{bmatrix}f_{\hspace{0.2pt}\texttt {x}}\\f_{\hspace{0.2pt}\texttt {X}}\end{bmatrix}=\begin{bmatrix}\Pi _{\hspace{0.2pt}\texttt {x}}\, \iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {X}}}- (F_{{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {X}}}}{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {x}}}}}^_{\hspace{0.2pt}\texttt {X}}\, \iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {X}}}+ F_{{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {I}}}}{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {x}}}}}^Ḭ(t)) \\\Pi _{\hspace{0.2pt}\texttt {X}}\, \nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {X}}}+ (F_{{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {X}}}}{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {x}}}}}\Pi _{\hspace{0.2pt}\texttt {x}}\, \nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {X}}}+ F_{{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {X}}}}{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {v}}}}}V(t)) \end{bmatrix}=\begin{bmatrix}F_\iota \, \iota _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {X}}}\\ F_\nu \, \nu _{\hspace{0.2pt}{\hspace{0.2pt}\texttt {X}}}\end{bmatrix}+\begin{bmatrix}-F_{{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {I}}}}{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {x}}}}}^Ḭ(t) \\ \quad F_{{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {X}}}}{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {v}}}}}V(t)\end{bmatrix}, {where F_\iota = \Pi _{\hspace{0.2pt}\texttt {x}}- F_{{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {X}}}}{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {x}}}}}^_{\hspace{0.2pt}\texttt {X}}, F_\nu = \Pi _{\hspace{0.2pt}\texttt {X}}+ F_{{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {X}}}}{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {x}}}}}\Pi _{\hspace{0.2pt}\texttt {x}}. With (\ref {eq:LTI5i}) and (\ref {eq:LTI5v}) this gives}0 &=\begin{bmatrix}F_\iota A_\iota \\ F_\nu A_\nu \end{bmatrix}x - \begin{bmatrix}F_\iota B_\iota \\ F_\nu B_\nu \end{bmatrix}\dot{x}+\begin{bmatrix}-F_{{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {I}}}}{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {x}}}}}^Ḭ(t) \\ \quad F_{{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {X}}}}{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {v}}}}}V(t)\end{bmatrix} \\&= Ax - B\dot{x}- u(t), \) 


which has the form (REF ), with
\(A &= \begin{bmatrix}F_\iota \,A_\iota \\F_\nu A_\nu \end{bmatrix},&B &= \begin{bmatrix}F_\iota \,B_\iota \\F_\nu B_\nu \end{bmatrix},&u(t) &= \begin{bmatrix}-F_{{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {I}}}}{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {x}}}}}^Ḭ(t) \\ \quad F_{{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {X}}}}{\hspace{0.2pt}\texttt {{\hspace{0.2pt}\texttt {v}}}}}V(t)\end{bmatrix}.\) 

Matrix rows, i.e. equations, have undergone reordering during this process, but this does not change the DAE's solutions so we may regard (REF , REF ) as a way, for the LTI case, to write (REF ).
Eigenfunctions are solutions of the corresponding homogeneous system \(Ax - B\dot{x}= 0\)  that have the form \(x(t) = ve^{\lambda t}\)  for some constant vector \(v\ne 0\)  and complex scalar \(\lambda \) .
Then \(\lambda ,v\)  solve the GEVP \(Av = \lambda Bv\) , so that \(\lambda \)  is a zero of the characteristic polynomial \(p(\lambda ) = \det (A-\lambda B)\) .
As a result of the tree being normal according to (A3), derivatives occur only in \(f_{\hspace{0.2pt}\texttt {c}}\)  and \(f_{\hspace{0.2pt}\texttt {L}}\)  of (REF ).
Only the corresponding rows of \(B\)  can have nonzeros, a total of \(n_{\hspace{0.2pt}\texttt {c}}+n_{\hspace{0.2pt}\texttt {L}}=\textsc {dof} \)  rows.
Hence by the formula for determinant as a sum of products over matrix transversals, the degree \(\deg p\)  of \(p(\lambda )\)  cannot exceed dof.
Non-singularity of the system Jacobian (a result of (A4)–(A6)) implies that \((A,B)\)  is a regular matrix pencil that has exactly \(n_{\hspace{0.2pt}\texttt {c}}+n_{\hspace{0.2pt}\texttt {L}}\)  eigenvalues or,
equivalently, that \(\det (A-\lambda B)\)  is a polynomial in \(\lambda \)  with degree \(n_{\hspace{0.2pt}\texttt {c}}+n_{\hspace{0.2pt}\texttt {L}}\)This results from the structure of \(B\)  and the Kronecker canonical form, see e.g. [12]}..
These are the finite eigenvalues of the GEVP; its remaining \((N{-}\textsc {dof})\)  eigenvalues are \(\infty \)  and are discarded (see [17]}).
For simplicity assume the problem has dof linearly independent eigenvectors \(v_i\)  (true e.g. if the finite eigenvalues are distinct).
Then each solution of \(Ax - B\dot{x}= 0\)  has a unique expansion
\(x(t) = \sum _{i=1}^{\textsc {dof}} c_i e^{\lambda _i t} v_i.\) 

Example 5 
The running example equations in () can be written in the LTI form (REF ) as
\(0 &=\begin{bmatrix}1/C_2 &-1/C_1 & 0 & 0 & 0 & 0 \\0 & 0 & 0 &-1/L_2 & 0 & 0 \\0 & 0 &1/L_1 &-1/L_2 & 0 & 0 \\0 &1/C_1 & 0 & 0 &R & 0 \\0 & 0 & 0 &-1/L_2 & 1 &-G\\0 & 0 & 0 & 0 &R & 1\end{bmatrix}\begin{bmatrix}q_{C_2}\\ q_{C_1}\\ \varphi _{L_1}\\ \varphi _{L_2}\\ \iota _{R}\\ \nu _{G}\end{bmatrix}-\begin{bmatrix}0 & 0 & 0 & 0 & 0 & 0 \\-1 &-1 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & -1&-1 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 \\\end{bmatrix}\begin{bmatrix}\dot{q}_{C_2}\\ \dot{q}_{C_1}\\ \dot{\varphi }_{L_1}\\ \dot{\varphi }_{L_2}\\ \dot{\iota }_{R}\\ \dot{\nu }_{G}\end{bmatrix}- \begin{bmatrix}-V(t)\\ -I(t)\\ -I(t)\\ V(t)\\ V(t)\\ 0\end{bmatrix} \\&= Ax - B\dot{x}- u(t).\) 
The characteristic polynomial is \(p = \det (A-\lambda B) = c(\alpha \lambda ^2 + \beta \,\lambda + \gamma )\)  where
\(c &=-R/(C_1C_2L_1L_2),&\alpha &=(G+1/R)(C_1+C_2)(L_1+L_2),&\beta &=C_1+C_2,&\gamma &=G+1/R.\) 

The finite eigenvalues are
\(\lambda _1,\lambda _2 &=-\delta \pm \sqrt{\delta ^2 -\eta }&\text{where } \delta &= \frac{1}{2(G+1/R)(L_1+L_2)} = \frac{\beta }{2\alpha },&\eta &= \frac{1}{(C_1+C_2)(L_1+L_2) }= \frac{\gamma }{\alpha }.\) 



Relation to other circuit formulations
Commonly used methods for modelling electrical circuits are the Modified Nodal Analysis (MNA) [18]}, the Modified Loop Analysis (MLA) [19]}, or the analysis via Branch-oriented Models (BOM) [10]} (or tree-based formulations) also known as hybrid analysis in the circuit literature [21]}, [22]}, [23]}, [24]}.
It has been shown in [25]}, [26]}, [27]}, [28]} that the index of the DAE resulting from the hybrid analysis is at most one for strictly locally passive components, in contrast to MNA and other nodal techniques, for which certain configurations (as \({\hspace{0.2pt}{\hspace{0.2pt}\texttt {C}}}{\hspace{0.2pt}{\hspace{0.2pt}\texttt {V}}}\) -cycles or \({\hspace{0.2pt}{\hspace{0.2pt}\texttt {L}}}{\hspace{0.2pt}{\hspace{0.2pt}\texttt {I}}}\) -cutsets) yield index two systems .
While MNA (and also MLA) naturally gives a DAE of port-Hamiltonian structure, see e.g. [29]}, [30]} for related results, simple examples show that the resulting DAE need not be SA-amenable, see [31]}.
In [31]} a branch-oriented model in port-Hamiltonian form based on a normal tree is presented that is closely related to the results presented here, but uses different variables.
A different definition for port-Hamiltonian descriptor systems is given in [33]}, [34]}.
With a slight increase in the dimension of the system
– by adding the variables \(\iota _{\hspace{0.2pt}\texttt {C}}\) , \(\nu _{\hspace{0.2pt}\texttt {l}}\)  and the corresponding relations \(\nu _{\hspace{0.2pt}\texttt {l}}=\dot{\varphi }_{\hspace{0.2pt}\texttt {l}}\) , \(\iota _{\hspace{0.2pt}\texttt {C}}=\dot{q}_{\hspace{0.2pt}\texttt {C}}\)  – our formulation can be transformed into this form.
Due to (A3) \(n_{\hspace{0.2pt}\texttt {l}}\)  and \(n_{\hspace{0.2pt}\texttt {C}}\)  are usually small.
The MNA formulation is based on Kirchhoff's laws written as
\(\nu = \widetilde{A}^\quad \text{ and}\quad \widetilde{A}\iota = 0,\) 

where \(\widetilde{A}\in \mathbb {R}^{({{n}{-}1})\times {b}}\)  is a reduced incidence matrix, i.e., \(A\)  with one row removed, and
\(\eta \)  denotes the \(({n}-1)\) -dimensional vector of node potentials.
By [35]} \(\widetilde{A}_\texttt {T}\)  is nonsingular, and
\( F &= -(\widetilde{A}_\texttt {T}^{-1} \widetilde{A}_{\texttt {T}^*})^ \quad \text{independent of which row was removed}.\) 

A transformation of (REF ) with \(\widetilde{A}_\texttt {T}^{-1}\)  gives
\(\nu = (\widetilde{A}_\texttt {T}^{-1}\widetilde{A})^\)  and \(\widetilde{A}_\texttt {T}^{-1}\widetilde{A}\iota = 0\) ,
where \(\mu =\widetilde{A}_\texttt {T}^\) .
Ordering \(\widetilde{A}\)  into tree and cotree parts as \([\widetilde{A}_\texttt {T}, \widetilde{A}_{\texttt {T}^*}]\)  , similarly the vectors \(\nu \)  and \(\iota \) , then
with (REF ) we get (REF ) and \(\mu \)  becomes the tree voltages \(\nu _\texttt {T}\) .
Similar relations also holds for the MLA equations. The BOM formulation is similar to \((\ref {eq:KL})\)  without using the properties of a normal tree.
The MNA equations contain the variable \(\iota _{\hspace{0.2pt}\texttt {V}}\) , while in the CpH approach \(\iota _{\hspace{0.2pt}\texttt {V}}\)  can be either separated out (as e.g. in Model 2) or considered together with the complementary variable \(\nu _{\hspace{0.2pt}\texttt {I}}\)  as output of the port-Hamiltonian system.
Also note that the CpH method does not require a ground node as do nodal-based methods.
Example 6 
The MNA formulation of the running example, with vertex 2 chosen as ground, is
\(C_1(\dot{\eta }_1 - \dot{\eta }_4)+ G \eta _1+\iota _{{\hspace{0.2pt}\texttt {L}}_1}+\iota _{\hspace{0.2pt}\texttt {V}}&=0\\C_2(\dot{\eta }_3 -\dot{\eta }_4)+ (1/R)\eta _3-\iota _{\hspace{0.2pt}\texttt {V}}&=0\\C_1(\dot{\eta }_4-\dot{\eta }_1)+C_2(\dot{\eta }_4-\dot{\eta }_3)+\iota _{{\hspace{0.2pt}\texttt {L}}_2}&=0\\-\iota _{{\hspace{0.2pt}\texttt {L}}_1}-\iota _{{\hspace{0.2pt}\texttt {L}}_2}+I(t)&=0\\L_1\mathrm {d}\iota _{{\hspace{0.2pt}\texttt {L}}_1}/\mathrm {d}t - \eta _1+\eta _5&=0\\L_2\mathrm {d}\iota _{{\hspace{0.2pt}\texttt {L}}_2}/\mathrm {d}t - \eta _4+\eta _5&=0\\\eta _1-\eta _3-V(t)&=0.\) 


It is SA-unamenable. However, choosing another vertex as ground leads to an SA-amenable system.


Numerical Implementation and Examples
Since the CpH DAEs in this paper have index \(\le 1\) , they can be solved by a DASSL-type code. Here we use Matlab's index-1 solver [basicstyle=]ode15i.
Subsection REF  outlines the main ideas of our Matlab implementation of the CpH formalism. The computation of a normal tree and it's Kron matrix is discussed in Subsection REF .
Subsection REF  discusses computing consistent initial conditions (CICs)
and
Subsection REF  presents numerical results.
Matlab implementation
We use three Matlab classes:
[\(\bullet \) ]
[basicstyle=]pHedge implements independent, possibly nonlinear, circuit elements.
An object is a vector of edge elements each having three fields:
[basicstyle=]type, indicating the kind of edge, one of V, C, R, G, L, I;
[basicstyle=]name, a text name of the edge;
[basicstyle=]value, a number or a function (handle).
If [basicstyle=]value is a number \(\xi \) , it defines a constant current or voltage for a V or I edge, and linear behaviour for a C, R, G or L edge, e.g. \(\nu =\xi \iota \)  for an R edge and \(\iota =\xi \nu \)  for a G edge.
If [basicstyle=]value is a function \(f\) , it defines \(V(t)\)  or \(I(t)\)  for a V or I edge, and the relation \(\nu =f(q)\) , \(\nu =f(\iota )\) , \(\iota =f(\nu )\) , \(\iota =f(\varphi )\)  for a C, R, G or L edge respectively.

A [basicstyle=]pHcircuit object is mainly a [basicstyle=]pHedge object together with an incidence matrix [basicstyle=]A stored in sparse form.
An important method is [basicstyle=]pHjoin which “solders together” specified nodes of one or more [basicstyle=]pHcircuit objects to create a new one combining all the edges of the input objects.
One can specify only a subset of a circuit's nodes as “visible”, i.e., usable by [basicstyle=]pHjoin.

[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 [36]}: 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 [37]}.
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) [38]} 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 [22]}. 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 [40]}, [41]}.
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 [42]} (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 [43]}, [44]}, where also the zero blocks of \(F\)  are exploited though less transparently than here.
The existence of a normal tree is shown in [44]}.
Our main proof is a distillation of that in [31]} 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. [29]}, [30]}.
In [29]}, 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 [30]} 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 [51]}.
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 [4]}.

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 [53]}.
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 [10]} 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) [1]}.
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 
([1]})
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 [57]}.
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 [57]}
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. [59]}, [60]}, and prove some that are special to this paper. The facts about matrix pairs are present, but not so directly stated, in [10]}.
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 [62]}.
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 [63]}.
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 [64]}
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 [65]}, [10]}
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. [67]}.
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 [10]} 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. [10]}).
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. [70]}).
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
