## 3 Generation of Equations
### Realizability Analysis
#### 3.2.2 Algorithm
\mathbf{u}_{k}\) for\(n=0\) to \(n_{NR}-1\)do \(\mathbf{r}_{n}\leftarrow\mathbf{w}_{k,n}-\mathbf{B_{w}}.\mathbf{z}(\mathbf{w}_{ k,n})-\mathbf{p}_{k}\) \(\mathcal{J}_{n}\leftarrow\mathbf{I_{d}}-\mathbf{B_{w}}.\mathcal{J}_{\mathbf{z}}( \mathbf{w}_{k,n})\) \(\mathbf{w}_{k,n+1}\leftarrow\mathbf{w}_{k,n}-\mathcal{J}_{n}^{-1}\cdot\mathbf{r}_ {n}\) end for \(\mathbf{w}_{k}\leftarrow\mathbf{w}_{k,n_{NR}}\) \(\mathbf{y}_{k}\leftarrow\mathbf{A_{y}}\cdot\mathbf{x}_{k}+\mathbf{B_{y}} \cdot\mathbf{z}(\mathbf{w}_{k})+\mathbf{C_{y}}\cdot\mathbf{u}_{k}\) \(\delta\mathbf{x}_{k}\leftarrow\mathbf{A_{x}}\cdot\mathbf{x}_{k}+\mathbf{B_{x}} \cdot\mathbf{z}(\mathbf{w}_{k})+\mathbf{C_{x}}\cdot\mathbf{u}_{k}\) \(\mathbf{x}_{k+1}\leftarrow\mathbf{x}_{k}+\delta\mathbf{x}_{k}\)  end for
```

**Algorithm 2**Simulation, with \(n_{t}\) the number of time-steps and \(n_{NR}\) the (fixed) number of Newton-Raphson iterations. ### Comparison with Standard Methods

In this section, the proposed approach (PHS structure combined with the discrete gradient method) is compared with two standard methods: the _trapezoidal rule_ (average of the vector field at \(\mathbf{x}(k)\) and \(\mathbf{x}(k+1)\), used in the WDF approach [7]) and the _midpoint rule_ (evaluation of the vector field at \(\frac{\mathbf{x}(k)+\mathbf{x}(k+1)}{2}\), suitable for any differential-algebraic system of equations). Both are known to preserve the passivity of linear undamped systems (see [35] for a detailed analysis). The updates associated with these three methods are given in Table 3. These methods are applied on the same conservative system \(\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}=\mathbf{J}_{\mathbf{x}}\cdot\
abla \mathcal{H}(\mathbf{x})\), the power balance of which is given by \(\frac{\mathrm{d}\mathbf{E}}{\mathrm{d}t}=0\) (with \(\mathtt{D}=\mathtt{S}=0\)). The comparison measure is then the relative error on energy \(\varepsilon(k)=\frac{\left|\mathcal{H}\left(\mathbf{x}(k+1)\right)-\mathcal{ H}\left(\mathbf{x}(k)\right)\right|}{\mathcal{H}\left(\mathbf{x}(0)\right)}\) for \(k\geq 0\). First, notice that for quadratic Hamiltonian \(\mathcal{H}(\mathbf{x})=\frac{\mathbf{x}^{\intercal}\cdot\mathbf{Q}\cdot \mathbf{x}}{2}\) with linear gradient \(\
abla\mathcal{H}(\mathbf{x})=\mathbf{Q}\cdot\mathbf{x}\), the three methods yield the same update:

\[\mathbf{x}(k+1)=\mathbf{x}(k)+T\cdot\mathbf{J}_{\mathbf{x}}\cdot\mathbf{Q} \left(\mathbf{x}(k)+\frac{\delta\mathbf{x}(k)}{2}\right). \tag{29}\]

As a consequence, these three methods induce the same frequency warping (see [36] for the analysis of the bilinear transform derived from the trapezoidal rule). We focus on the nonlinear case. For comparison, we choose a simple nonlinear conservative system with state \(\mathbf{x}=(x_{1},x_{2})^{\intercal}\), non-quadratic Hamiltonian

\[\mathcal{H}(\mathbf{x})=10\,\log\Big{(}\cosh(x_{1})\Big{)}+\Big{(}\cosh(x_{2} )-1\Big{)}, \tag{30}\]

and canonical skew-symmetric matrix \(\mathbf{J}_{\mathbf{x}}=\left( \begin{array}{cc}0&-1\\ 1&0\end{array} \right)\). Figure 4: _Cont._ \begin{table}
\begin{tabular}{c c} \hline \hline
**Method** & **Update** \\ \hline Trapezoidal rule & \(\mathbf{x}(k+1)=\mathbf{x}(k)+T\cdot\mathbf{J}_{\mathbf{x}}\cdot\frac{\
abla \mathcal{H}\left(\mathbf{x}(k)\right)+\
abla\mathcal{H}\left(\mathbf{x}(k+1) \right)}{2}\) \\ \hline Midpoint rule & \(\mathbf{x}(k+1)=\mathbf{x}(k)+T\cdot\mathbf{J}_{\mathbf{x}}\cdot\
abla \mathcal{H}\left(\frac{\mathbf{x}(k)+\mathbf{x}(k+1)}{2}\right)\) \\ \hline PHS with discrete gradient & \(\mathbf{x}(k+1)=\mathbf{x}(k)+T\cdot\mathbf{J}_{\mathbf{x}}\cdot\
abla_{ \mathcal{A}}\mathcal{H}\left(\mathbf{x}(k),\mathbf{x}(k+1)\right)\) \\ \hline \hline \end{tabular}
\end{table} 
Table 3: Updates for the three methods considered in §4.3. PHS stands for port-Hamiltonian system. In each case, the resulting implicit equations are solved by Python iterative solver (see [37]). The Python code is available at the url given in [38].