1 Introduction↩︎

In the simulation and analysis of electric circuits, structure-preserving formulations have gained increasing attention [1]–[3] due to their ability to reflect the physical characteristics at both the continuous and discrete levels. Within this context, the port-Hamiltonian (pH) framework facilitates the modeling and the simulation of dynamical systems [1] with an energy balance.

To solve such systems efficiently while maintaining their structural properties, operator splitting methods are a promising approach [4]. In the prior work [5], the introduced JR-decomposition separates the pH system into conserved and dissipative subsystems. However, for electric circuit models, the applicability of the standard JR-decomposition remains restrictive, for instance due to the presence of voltage sources.

In this work, we introduce an enhanced JR-decomposition for electric circuit differential-algebraic equations (DAEs) that extends the applicability of the method to general lumped circuit topologies. We consider the MNA [6] formulation and demonstrate the applicability of the proposed approach.

This paper is structured as follows: we introduce the pH framework and JR-decomposition in Section 2. In Section 3, this decomposition is transferred to MNA models and enhanced to enable usage for a broader class of circuits. Finally, numerical results are presented (Section 4).

2 PH formulation and JR-decomposition↩︎

Port-Hamiltonian (pH) systems are network descriptions, which include energy routing, energy dissipation and energy input with collocated output. This class is closed under the pH coupling of these systems. Therefore, it serves as a coupling framework also for systems from various physical domains, e.g. electric circuit models or multibody systems. Thus, the prime property of pH systems is the energy conservation/dissipation. Very naturally, this yields differential-algebraic equations (DAEs) in time domain.

2.1 Linear port-Hamiltonian systems↩︎

For linear models, this amounts to the following description: given a time interval \(I\subseteq \mathbb{R}\) and a state variable \(x: I \to \mathbb{R}^n\), we consider linear pH-DAEs of the form \[\begin{align} \label{eq:pHDAE} E \dot{x} &= (J-R)x + Bu(t), \quad x(0)=x_0 \\ y &= B^\top x \nonumber \end{align}\tag{1}\] with matrices \(J=-J^\top\), \(R=R^\top \geq 0\), and \(E=E^\top \geq 0\) from \(\mathbb{R}^{n\times n}\), as well as \(B\in\mathbb{R}^{n\times m}\), input \(u(t)\in \mathbb{R}^m\), collocated output \(y(t)\in \mathbb{R}^m\) and (consistent) initial value \(x_0\in\mathbb{R}^n\). Moreover, the solution \(x\) of DAE 1 satisfies the dissipation inequality \[\frac{d}{dt}\mathcal{H}(x(t)) \leq y(t)^\top u(t), \quad \mathcal{H}(x) = \frac{1}{2}x^\top E x.\] For constant coefficient case of 1 , the DAE index is at most 2, see [7].

2.2 JR-decomposition for ODEs and DAEs↩︎

The JR-decomposition is a structural splitting approach for pH-DAEs. It was introduced and discussed in [5]. This decomposition enables tailored, structure-preserving numerical schemes. In fact, the port-Hamiltonian system 1 can be decomposed as follows: \[\label{eq:pHS-splitting-basis} E \dot{x} = f \quad \text{as} \quad \left\{ \begin{array}{l} E \dot{x} = f_1 \\ E \dot{x} = f_2 \end{array} \right\} \quad \text{with} \quad \begin{array}{l} f_1 = Jx, \\ f_2 = -Rx+Bu. \end{array}\tag{2}\] Thus \(f_1\) states the energy-conserving part, whereas \(f_2\) gives the dissipative contribution with sources. This is the so-called JR-decomposition. It enables the use of energy conserving methods for the first part with \(f_1\) and dissipative/general techniques for the second part.
Index-0. If \(E>0\), system 1 is an ODE system and standard theory for splitting applies. In this case, the JR-decomposition is well-defined.
Index-1. Here, \(E\) is not regular. Thus we have two subsystems with at least one DAE subsystem. Standard DAE decoupling techniques double algebraic constraints, see [3], [5]. A well-posed splitting strategy can be constructed using the dissection approach from [8], which separates differential and algebraic components. Depending on the structure of the pH-DAE, the JR-decomposition can occur in two cases:

1. \(R, B\) contribute only to the differential part, or

2. \(J\) contributes only to the differential part.

In both cases, the resulting decomposition yields two subsystems associated with the conservative and dissipative parts, respectively; for details see [5]. The corresponding splitting can then be integrated using Strang splitting, which preserves second-order convergence.
Index-2. This case is more challenging, since there are hidden algebraic constraints. In the present work, we restrict ourselves to index-1 formulations.

Next, we discuss the application of JR-decomposition to circuits.

3 PH modeled electric circuit and JR-decomposition↩︎

We consider linear electric networks with two-terminal elements: capacitances (\(C\)), inductances (\(L\)), resistances (\(G\)), independent current sources \(\imath_S(t)\), and independent voltage sources \(v_S(t)\), for short linear RLC-networks with inputs. Using modified nodal analysis (MNA), [6] these circuits can be represented as port-Hamiltonian system.

3.1 MNA circuit model and pH system↩︎

In what follows, the matrix \(A\) denotes the incidence matrix and is partitioned according to the circuit elements of our RLC-network as \(A=[A_C, A_L , {A_G}, A_I, A_V]\). Using the standard MNA formulation, the circuit equations can be written as the linear pH-DAE \[\begin{align} \label{eq:mna-model} &\underbrace{\begin{pmatrix} A_C C A_C^\top & 0 & 0 \\ 0 & L & 0 \\ 0& 0 & 0 \end{pmatrix}}_{\displaystyle E:=} \begin{pmatrix} \dot{e} \\ \dot{\imath}_L \\ \dot{\imath}_V \end{pmatrix} = \Biggl( \underbrace{\begin{pmatrix} 0 & - A_L & -A_V \\ A_L^ \top & 0 & 0 \\ A_V^\top & 0 &0 \end{pmatrix}}_{\displaystyle J:=} - \underbrace{\begin{pmatrix} {A_G} G {A_G^\top} & 0 & 0\\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}}_{\displaystyle R:=}\Biggl) \begin{pmatrix} e \\ \imath_L \\ \imath_V \end{pmatrix} + \underbrace{\begin{pmatrix} -A_I & \phantom{-}0 \\ \phantom{-}0 & \phantom{-}0 \\ \phantom{-}0 & -I \end{pmatrix}}_{\displaystyle B:=} \underbrace{\begin{pmatrix} \imath_S(t) \\ v_S(t) \end{pmatrix}}_{u(t):=} \end{align}\tag{3}\] equipped with the output relation \[y = \begin{pmatrix} -A_I^\top & \phantom{-}0 & \phantom{-}0 \\ \phantom{-}0 & \phantom{-}0 &-I \end{pmatrix} x = \begin{pmatrix} - A_I^\top e \\ - \imath_V \end{pmatrix}.\] The state variables collected in \(x^\top =\bigl(e^\top,\, \jmath_L^\top,\, \jmath_V^\top\bigr)\) are: node potentials \(e\), currents through inductances \(\imath_L\) and currents through voltage sources \(\imath_V\). Notice, we use \(\dot{\imath}(t)=\frac{\text{d}}{\text{d}t}\imath (t)\). As above, the linear system 3 has maximal index-2; see also [1].

In fact, we have index-0 (an ODE) if \(A_C^\top\) has full rank and there are no voltage sources. The system is of index-1 if there exists neither a \(CV^{+}\)-loop (i.e., a loop of capacitors and voltage sources with at least one voltage source) nor an \(LI\)-cutset (i.e., a cutset composed by inductors and currents sources only). Otherwise, the system is of index-2; see [9]. We note that index-2 circuit systems can be treated in the loop-cutset framework where index-2 variables are merely computed in a post-processing, see [10].

3.2 JR-decomposition for index-1 MNA models↩︎

Now, we analyze the index-1 case of MNA equations in the context of the JR-decomposition. We denote by \(K_E\) the projector onto \(\ker(E)\), which reads for the MNA equations: \(K_E=\text{blkdiag}(K_C, 0, I)\), where \(K_C\) denotes a projector onto \(\ker (A_C^\top)\).

3.3 Enhanced JR-decomposition↩︎

As we have seen in Section 3.2, the assumptions required for the JR-decomposition for circuits in MNA are somewhat restrictive, e.g. no voltage sources can be present. To relax these assumptions and extend the method to a broader class of problems, we introduce for the MNA equations 3 \[\label{eq:enJR} \begin{equation} J_{en} = \begin{pmatrix} 0 & -A_L & 0 \\ A_L^\top & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \qquad R_{en} = \begin{pmatrix} A_G G A_G^\top & 0 & A_V \\ 0 & 0 & 0 \\ -A_V^\top & 0 & 0 \end{pmatrix}. \end{equation} This way, the voltage sources are moved to the R-part and it gives the decomposition \begin{equation} f_1(x)= J_{en}(x), \qquad f_2(x)=-R_{en}x+Bu, \end{equation}\tag{4}\] which we call enhanced JR-decomposition. It satisfies similar properties: \[J_{en} \;\text{ is skew-symmetric}, \qquad R_{en} \;\text{ has a positive symmetric part, i.e., } \;R_{en}^\top +R_{en} \geq 0.\]

Again, for consistency, we accommodate the following two cases:

1. \(R_{en}, B\) contribute only to the differential part;

2. \(J_{en}\) contributes only to the differential part.

3.4 Convergence Theory↩︎

For the enhance JR-decomposition, we have the following

Next, we apply the proposed method to a numerical example, demonstrating its convergence and practical performance.

4 Numerical Results↩︎

We consider an electric circuit, which is roughly build of two short and coupled transmission lines with an input \(v(t)\) and a load resistance \(R_L\), see Fig. 2. For all parameters, see caption of Fig. 2. The unknown state variables are the node potentials \(e_i,\; i=1,\, \dotsc,\, 6\), the currents (through inductors) \(\imath_1, \, \imath_2\) and the current through the voltage source \(\imath_V\). Modeling the electric network via modified nodal analysis yields a pH-DAE of index-1 for \(x^\top = (e_1,e_2,e_3,e_4,e_5,e_6,\imath_1, \imath_2, \imath_v)\).

In the top row of Fig. 3, we observe the convergence behavior of the Strang splitting for the enhanced JR. The differential variables, in both the 3-mid and iE-mid-eE variants, display consistent second-order convergence. The 2-stage Radau IIA method (3-Radau) shows second-order convergence for both differential and algebraic variables, which aligns with its theoretical order (for index-1 DAEs) when combined with Strang splitting. The non-symmetric variant iE-mid-iE shows only first order as expected. For comparison, we also solve the original (non-split) system using the classical backward differentiation formula of order two (BDF2), a widely used method in circuit simulation. The results indicate that BDF2 exhibits the same second-order convergence behavior as the 3-mid and 3-Radau variants, although with a larger error constant for the differentiable variables. In contrast, the iE-mid-iE method exhibits a clear first-order convergence, as expected from the sole use of the implicit Euler scheme in solving the \(J_{en}\)-subsystem. For the algebraic variables, 3-mid and iE-mid-eE show second order, which is in line with theory. In the bottom row of Figure 3, we observe the satisfaction of the dissipation inequality for Strang splitting with 3-mid variant (the red line giving the dissipation restriction to non-negative values), as well as the 3-Radau variant with smaller time step sizes \(\lesssim 10^{-12}\). However, we observe a violation of the dissipation inequality when using the 3-Radau variant with step sizes \(\geq 10^{-10}\). This is due to the structure-preserving property of the implicit midpoint rule for the conserved part; which is not the case for the Radau method. Moreover, we investigate the dissipation inequality for the solution variants mid-iE-mid and mid-Radau-mid with time steps \(10^{-12}\) and \(10^{-11}\), respectively. The results indicate that both approaches capture the dissipative behavior of the system and remain close to the reference dissipation. In particular, the mid-Radau-mid scheme shows an improved approximation of the dissipation compared to the 3-Radau method with a fixed step size. This can be explained by the fact that, each subsystem in the splitting framework is discretized using a method adapted to its structural properties. This structure-aware choice allows for a more accurate reproduction of the energy balance, for larger step sizes.

5 Conclusions↩︎

We considered the applicability of JR-decomposition-based operator splitting to electric circuit systems in the MNA modeling framework. We highlighted the limitations of the classical JR-splitting when applied to certain circuit topologies. To overcome these limitations, we proposed an enhanced version of the JR-decomposition, which extends the applicability of the splitting technique to a wider range of circuits, including those with voltage sources. A numerical experiment has been conducted to assess the convergence orders, which confirm the theoretical findings.

References↩︎