Multidomain modeling of nonlinear electromagnetic circuits using wave digital filters

Accurate models of electromagnetic systems can be derived by coupling electric and magnetic equivalent circuits together. The different nature of such physical domains constitutes a big challenge that puts a continuous strain on software simulators. To face this problem, a simulation approach based on Wave Digital Filters (WDFs) is proposed in this manuscript. The method is employed to solve nonlinear electromagnetic systems containing complex magnetic equivalent circuits while maintaining the modularity of the electric and magnetic subsystems. The nonlinearities can be locally handled, enabling the possibility to choose a dedicated model for each one of them. The proposed algorithm is a hierarchical generalization of the Scattering Iterative Method, which has shown, over the past few years, promising performance for the simulation of large nonlinear circuits. In addition, the method constitutes a further step towards the development of novel general‐purpose circuit simulators based on WDF principles. In a comparison with mainstream circuit simulation software, the proposed approach turns out to be considerably faster and thus particularly promising for parametric analyses of electromagnetic systems.


| INTRODUCTION
Simulation and performance forecasting of large electrical and electronic systems continue to be one of the challenges in the design of new devices. The progress in electronics, in fact, is strictly dependent on how fast, accurate, and reliable the simulation algorithms are. Reference systems to be simulated, however, are often characterized by multiple interdependent physical domains. In particular, in a large number of applications, magnetic components such as inductors or transformers are integrated together with electrical or electronic circuits. Suffice to think of the fundamental role played by such devices in the context of power converters, where increasingly complex operation conditions are taken into account. 1,2 Equivalent circuits of electromagnetic problems can be obtained exploiting, for example, the Partial Element Equivalent Circuit (PEEC) 3,4 technique or the Magnetic Equivalent Circuit (MEC) method. [5][6][7] While the first leads to circuits in the electrical domain, the latter allows us to embed detailed physics-based models of magnetic devices and to simulate them together with electrical subsections in a unified (homogeneous) fashion. Hence, MECs are suitable for multidomain simulations. Nonetheless, electromagnetic systems share some features that make systemlevel simulations extremely computationally intensive: (i) magnetic devices are inherently characterized by distributed phenomena, and thus, aiming at an accurate modeling, MECs usually involve a large number of elements; (ii) MECs feature, in addition, a high number of nonlinear elements since magnetic materials exhibit strongly nonlinear characteristics; (iii) magnetic and electric subsections need to be simultaneously simulated.
The most widespread approach to multidomain simulation (e.g., used in general-purpose solvers based on the modified nodal analysis [MNA] such as LTspice 8,9 ) involves the solution of one large nonlinear system, encompassing equations from both the magnetic and the electric domains, by means of multivariate iterative solvers, for example, Newton-Raphson (NR). These iterative methods are usually based on large Jacobian matrices, 10 whose size is typically in the order of the number of nodes or meshes of the reference circuit, and require thus to perform computationally demanding operations that might impair both the efficiency and the robustness of the simulation. Moreover, such an approach, even though general, does not preserve the modularity of the two physical domains, which can, instead, be exploited to effectively solve the circuit. A way to enhance the efficiency of the simulation is to apply reduced-order modeling techniques in order to extract the essential system dynamics. 11 Such techniques, however, can introduce inaccuracies due to the truly nonlinear nature of MECs that may undermine the reliability of the overall results.
In this paper, with the purpose of providing efficient simulations of the aforementioned multidomain systems, we propose a technique based on a WD relaxation method, called Scattering Iterative Method (SIM), which, recently, has shown promising performance for solving circuits with an arbitrary number of nonlinear elements. 12 The method is based on the Wave Digital Filter (WDF) theory, introduced by Fettweis in the 1970s. 13 The WDF approach finds its roots in the traveling-wave formulation of lumped electrical elements introduced by Belevitch 14 and shares some features with the transmission-line modeling (TLM) method. 15 Circuit elements and topological interconnections (called connection networks) are described by means of scattering equations and modeled as input/output blocks. Wave Digital (WD) structures can be then derived in a block-wise fashion following, for example, one of the procedures explained in Werner et al. 16 or in Bernardini et al. 17 The signal variables into play are no longer voltages and currents (Kirchhoff variables) but incident and reflected waves. This peculiar description of the reference circuit leads to several numerical advantages, first and foremost the possibility to separately model topology and circuit elements, enhancing the modularity of the circuit representation.
The WD approach has already shown to be very effective for the solution of linear electromagnetic circuits obtained through the PEEC technique, 18 raising interest in the solution of nonlinear electromagnetic systems. As far as the simulation of nonlinear circuits in the WD domain is concerned, several iterative methods have been developed. For example, in Olsen et al., 19 the nonlinearities are grouped and solved with a multivariate NR solver, whereas in Schwerdtfeger and Kummert 20 the modularity is preserved using the multidimensional WDF formalism. Among those, SIM is a fixed-point method able to solve circuits with an arbitrary number of one-port nonlinear elements using independent scalar solvers (e.g., scalar NR solvers). 12,21 This implies several powerful features that differentiate SIM from other WD techniques, 19,20 such as greater efficiency, greater robustness, and the possibility of parallelizing many operations. As discussed in Bernardini et al., 12,22 SIM turned out to be more efficient than standard Spice-like simulators when exploited for the analysis of photovoltaic arrays characterized by thousands of nonlinear elements. In other studies, 21,23-25 the method is applied, instead, to emulate diode-based audio circuits for virtual analog applications, showing performance comparable to state-of-the-art techniques. However, in order to employ SIM for solving magnetic and electric subcircuits coupled together, several adjustments and extensions are required.
In this manuscript, we thus provide a generalization of SIM, referred to as Hierarchical SIM (HSIM), which is able to solve structures consisting of an arbitrary number of WD junctions, preserving accordingly the modularity of the physical domains. Hierarchical generalizations of relaxation methods have been already presented in the literature, for example, in the context of the Waveform Relaxation (WR) method for circuit simulation, 26,27 proving to perform, in many cases, better than the original method. It is worth pointing out, however, that the proposed HSIM and hierarchical WR are two different kinds of relaxation methods, since the first operates in the WD domain, while the second in the Kirchhoff domain. Moreover, while in the hierarchical WR method, each level of the tree is analyzed independently of the others; according to HSIM, all levels are analyzed in the same iterative procedure. In particular, the WD structures used in the proposed HSIM present themselves as generalizations of the WD binary connection trees discussed in Sarti and De Sanctis 28 and are scanned following a forward-backward computational flow. The multidomain implementation is achieved by introducing a WD block able to couple the magnetic and the electric subsystems together.
Several nonlinear and non-ideal magnetic effects are taken into account, such as core saturation, eddy currents, and flux leakage. Moreover, the magnetic nonlinearities are modeled in the WD domain following a data-driven approach based on Canonical Piecewise-Linear (CPWL) functions. 29,30 It is worth highlighting that, being based on MECs, the methodology discussed in the present manuscript can be employed to solve the subset of electromagnetic systems that can be modeled by means of equivalent circuits.
The paper is organized as follows. Section 2 describes the magnetic equivalent model employed for the simulation of a reference nonlinear transformer. Section 3 provides background knowledge on traditional WDFs and explains a particular modeling technique, based on CPWL functions, used to implement nonlinear reluctances in the WD domain. Section 4 shows the general WD implementation of a two-port junction able to couple electric domain and magnetic domain, whereas Section 5 presents the HSIM algorithm. In Section 6, the proposed methodology is exploited to efficiently simulate a full-wave rectifier, whose results are then compared in terms of accuracy and efficiency to those of two general-purpose solvers based on MNA, 8 such as MathWorks Simscape 31 and LTspice. 9 Conclusions are finally drawn in Section 7.

| THE MAGNETIC MODEL
The modeling procedure considered in this paper is based on the well-known MEC method, 7,32 widely studied in the literature, here briefly re-introduced to show the peculiar aspects of the proposed approach.
Let us consider, without any loss of generality, a two-winding transformer having an EI geometry ( Figure 1A). By applying the traditional voltage/magneto-motive-force analogy, 33 we then derive the equivalent magnetic circuit shown in Figure 1B. The model contains several reluctances that can be divided into three sets: nonlinear core reluctances (ℛ 1 ,…, ℛ 10 ), which model the nonlinear B À H curve of the magnetic material, where B is the flux density and H is the magnetic field, air-gap reluctances (ℛ 11 , …, ℛ 13 ), and leakage reluctances (ℛ 14 ,…, ℛ 20 ). ℛ a1 and ℛ a2 are the auxiliary leakage reluctances and constitute a special case since the potential at their terminals is known. In addition, we suppose the magnetic core to be laminated by means of n lam ¼ 20 magnetic sheets.
Due to eddy currents, the magneto-motive force (m.m.f.) is not uniform along the cross-section of the core. The main approach to take into account such a phenomenon, followed for example in Bottauscio et al., 7 requires to divide its cross-section into different "concentric" areas where both electric and magnetic quantities are assumed constant. Considering for illustration purpose the case of four areas, the cross-section of each lamination can be modeled by means of the ladder network shown in Figure 2. Elements ℛ lad1 , ℛ lad2 , ℛ lad3 , and ℛ lad4 are the nonlinear reluctances of the four areas, whereas conductances ℒ 1 , ℒ 2 , and ℒ 3 model the losses due to the eddy currents. In this scenario, each nonlinear reluctance ℛ 1 , …, ℛ 10 of the MEC model is substituted with the ladder network in Figure 2, whose components should be tuned in order to describe the geometry of the corresponding core piece. The formulas employed to compute ℛ 11 , …,ℛ 20 , ℛ a1 , and ℛ a2 can be found, instead, in Cale et al. 34 The final result is a MEC composed of 20 laminations that in turn are composed of two m.m.f. sources, 12 linear reluctances (air-gap and leakage reluctances), and 10 ladder networks of seven elements each. The total number of one-port elements of the MEC model would thus be 84 Â 20 ¼ 1680, that represents the complexity of the model. In Section 4, we will show how the proposed approach can be used to reduce this complexity. Nevertheless, the amount of nonlinear one-ports remains quite high and an efficient simulation tool is thus required, especially whether the transformer model is used when simulating large nonlinear circuits.

| MODELING CIRCUITS IN THE WAVE DIGITAL DOMAIN
The peculiarity of the Wave Digital approach lies in the choice of waves as circuit variables. In fact, each pair of port variables in the Kirchhoff domain, that is, a port voltage v and a port current i, is substituted with a pair of waves, whose most used definition is the following: 13,35 where a and b are the incident and the reflected waves, respectively, while Z is a free parameter called port resistance.
The inverse mapping of Equation (1) is which holds true if and only if Z ≠ 0. The free parameter Z constitutes a powerful degree of freedom in the port description since it can be arbitrarily selected. A proper choice of Z is fundamental for enhancing the speed of numerical solutions of WD structures. 12 Finally, WDF principles allow us to derive a port-wise description of the reference circuit and to separately model circuit elements and topological connection networks using one-port or multi-port WD blocks, whose input-output scattering relations can be derived as follows.

| Topological connection networks
The interconnections among elements are handled by a multi-port WD junction that is described through a set of scattering relations. Let us consider an N-port topological connection network. 17,36 Once v t is defined as the vector of independent port voltages of size χ Â 1 and i l as the vector of independent port currents of size ψ Â 1, where χ þ ψ ¼ N, we can express the vector of port voltages v ¼ ½v 1 , …, v N T and the vector of port currents i ¼ ½i 1 F I G U R E 2 Ladder network modeling the cross-section of each lamination B is the fundamental loop matrix, whereas Q is the fundamental cut-set matrix. Given that topological connection networks are reciprocal, the orthogonality property QB T ¼ 0 (or BQ T ¼ 0) holds true. 17,37,38 Matrix Q of size χ Â N and matrix B of size ψ Â N can be derived by performing a tree-cotree decomposition of a directed graph representation of the reference circuit. 39 The WD realization of the connection network is an N-port scattering junction characterized by a vector of incident waves a ¼ ½a 1 , …, a N T and a vector of reflected waves b ¼ ½b 1 where Z ¼ diag½Z 1 ,…, Z N is a diagonal matrix having port resistances as diagonal entries. The scattering relation between a and b is where S is a scattering matrix of size N Â N. According to the method discussed in Bernardini et al. 17 and Martens and Meerkötter, 40 matrix S can be computed choosing one of the following two equivalent formulas: where I is the N Â N identity matrix. If χ > ψ, Equation (7) is computationally cheaper than Equation (6), since it requires the inversion of a ψ Â ψ matrix rather than a χ Â χ matrix. If ψ > χ, the opposite holds true.

| Linear circuit elements
The discrete-time model of a large class of linear one-port circuit elements, including resistors, resistive voltage/current sources, and dynamic elements, such as capacitors and inductors discretized using stable linear multi-step methods, 21 is characterized by the following equation where k is the sampling index, v[k] is the port voltage, i[k] is the port current, V g ½k is a voltage parameter, and R g ½k is a resistive parameter. Substituting in Equation (8) the linear transformation of variables shown in Equation (2), a WD realization of such a generic linear one-port can be written as If we set Z½k ¼ R g ½k, the instantaneous dependence between b[k] and a[k] is eliminated and Equation (9) reduces to b½k ¼ V g ½k; in this case, the linear one-port element is said to be adapted. 13,35

| Nonlinear circuit elements
Contrary to linear one-ports, in general, nonlinear circuit elements cannot be adapted; hence, the dependence of the reflected wave b[k] from the incident wave a[k] cannot be eliminated. However, in order to minimize the instantaneous dependence of b[k] on a[k], the free parameter Z can be dynamically updated during the simulation in such a way that it is set as close as possible to the slope of the line tangent to the operating point on the v À i characteristic of the element, as already discussed in other studies. 12,21,[23][24][25] The possibility given by the WD formalism to separately model topology and circuit elements allows us to implement each nonlinear element with a dedicated technique, as confirmed by the examples described in the following subsections.

| Diodes
The diodes considered in this work are described with the extended Shockley model, 12 whose implicit function is where, for the sake of clarity, the sampling index k is omitted and I s is the saturation current, η is the ideality factor, V t is the thermal voltage, whereas R s and R p are the series and parallel resistances of the p-n junction. Since Equation (10) is a transcendental function, finding a scattering relation is not straightforward. Possible approaches for the WD implementation of Equation (10) are provided in other studies. 12,21,23,24,41 Here, instead, we provide an explicit WD scattering relation based on the Wright ω function. [42][43][44] Substituting Equation (1) in Equation (10) and solving for b, we get

| Nonlinear reluctances
With the purpose of having WD models of nonlinear reluctances able to describe off-the-shelf magnetic components, we assume the considered EI core to be built with a Voestalpine isovac electrical steel (235-35 A), 45 whose nonlinear B À H characteristic is depicted in Figure 3A. The curve exhibits saturation for high-amplitude induction fields. We fit the original isovac characteristic selecting eight points on the B À H curve and employing a first-order interpolator. We do this by using a highly-flexible modeling methodology based on CPWL functions. 29,30 CPWL representations are explicit and global, and they can be implemented with low storage requirements as well as low computational cost. In addition, they allow us to implement single-valued functions with arbitrary accuracy, without resorting to look-up tables or iterative solvers. The accuracy, in fact, can be improved by increasing the number of linear segments employed in the description. As an example, Figure 3B shows the result achieved with an interpolation of just eight points: the maximum value of the absolute error ε between the reference nonlinearity and the piecewise linear curve is ε max ¼ 0:28 T, whereas the mean value is ε ¼ 94:4 mT. A CPWL representation of a generic scalar WD mapping b = h(a) (without jump discontinuities) can be written as where J + 1 is the number of segments and J is the number of points with coordinates (a j , b j ). All the other parameters are, instead, defined as where m j is the slope of the jth segment and it is computed as Without any loss of generality, let us assume that the slopes of the segments with index j ¼ 0 and j ¼ J, corresponding to the intervals [À ∞, a 1 ] and [a J , + ∞], to be m 0 ¼ m 1 and m J ¼ m JÀ1 , respectively. Let us define p B ¼ ½p B1 ,…, p BJ T as the vector of the B-coordinates and p H ¼ ½p H1 , …, p HJ T as the vector of the H-coordinates, both sorted in increasing order since we are dealing with monotonically increasing characteristics. For each nonlinear reluctance, we can compute the vectors of coordinates ℱ p ¼ ½ℱ 1 , …, ℱ J T and ϕ p ¼ ½ϕ 1 , …, ϕ J T referred to the m.m.f. and the magnetic flux, imposing where Γ is the length of the magnetic path in meters and Λ is its cross-section in square meters. We then collect the slopes (i.e., reluctance values) of the different segments belonging to the same ℱ À ϕ curve in the vector ℛ p ¼ ½ℛ 1 , …, ℛ JÀ1 T . Each parameter ℛ j is computed as where δℱ j ¼ ðℱ jþ1 Àℱ j Þ and δϕ j ¼ ðϕ jþ1 À ϕ j Þ. Remembering that ℱ is the magnetic equivalent of the electric voltage and ϕ is the magnetic equivalent of the electric current, the corresponding points in the WD domain will have coordinates where Z p is a free parameter properly chosen among the ℛ j values, which are all precomputed. If Z p > 0, it can be verified that the set of points in the WD domain is still sorted in increasing order. 30 Once we have derived the vectors a p and b p , the remaining parameters required by the CPWL representation can be computed according to Equation (13).

| MAGNETIC/ELECTRIC JUNCTION
It is desirable to connect the magnetic subcircuits to the electric subcircuits in a modular fashion, such that a flexible multidomain model of the reference circuit can be developed. This is accomplished by designing a two-port element with memory, that we call Magnetic/Electric (ME) junction, able to convert magnetic variables into electrical variables and vice versa. 33 As discussed in Section 2, since we are considering a laminated core, we should in principle separately model each lamination. However, considering the approximation that all of them have the same MEC model and that there is no flux exchange among them, we can simulate just one lamination and take into account the whole stack using a multiplication factor n lam in the electric domain. The constitutive equations of the ME junction in the continuous-time domain therefore become vðtÞ ¼ Àn t n lam dϕðtÞ dt ℱðtÞ ¼ n t iðtÞ where n t is the number of winding turns.

| WD implementation of the ME junction
The ME junction is characterized by two ports: the port facing the electrical subcircuit called "electric (or electrical) port" and the port facing the magnetic subcircuit called "magnetic port." Let us mark the relative signals with the subscript "e" and "m," respectively. Therefore, we define the electrical and the magnetic variables in Equation (18) as In order to implement the ME junction in the discrete-time domain, the differential relation in Equation (18) is numerically approximated by using the Backward Euler discretization method as where T s is the sampling period (or time step) and k is the discrete-time sampling index. In order to derive a WD realization of the ME junction, we write the Kirchhoff variables as a function of the wave variables a e , b e , a m , and b m v e ½k ¼ a e ½k þ b e ½k 2 , where Z e and Z m are the free parameters referred to the electric port and the magnetic port, respectively. By replacing Equation (20) into Equation (19) and solving for the reflected waves, we obtain the following system of equations: where , and β½k ¼ Z e ½k n t þ n t n lam T s Z m ½k . It is worth noticing that, in case the port resistances Z e and Z m are varied, the two matrices S ME and M ME need to be recomputed.
In order to speed up the solution of WD structures, it is desirable to remove as many implicit relations as possible. Since the two diagonal entries of matrix S ME are equal, we can make both the electric and the magnetic ports reflection-free at each sampling step k by properly setting the free parameters Z e and Z m . This is done by satisfying the following constraint: either setting Z e as or, equivalently, setting Z m as For example, Figure 4 shows an adapted ME junction, where the two T-shaped symbols indicate that the two ports are made reflection-free.
Backward Euler (BE) is the simplest implicit method that can be used for the discretization of the time derivative in Equation (18). Although BE shows good properties as far as stability is concerned (it is an L-stable method 46 ), its Local Truncation Error (LTE) is OðT 2 s Þ, and thus, in certain application scenarios, T s must be very small to match the desired accuracy. In these scenarios, linear multi-step discretization methods, such as Backward Discretization Formulas (BDFs), 21,46 can be taken into account for the modeling of the ME junction, since, at the cost of worse stability properties, a lower LTE (and thus better accuracy) can be accomplished. Starting from Equation (18) and following the same design procedure applied with the BE method, we can derive WD models of the ME junction by means of BDFs of order Q. Considering the fixed-step scenario, Equation (21) becomes where S ME is the same defined in Equation (22) F I G U R E 4 ME junction adapted employing Equation (25). The adaptation process eliminates at once the instantaneous reflections at both the ME ports, and it is symbolically represented by a T-shaped stub The coefficients φ q weighing past samples of a m and b m in Equation (26) can be computed, instead, starting from the coefficients of the considered BDF. 46 Table 1 presents the φ q values obtained taking into account zero-stable BDF methods. 46 It is worth pointing out that, since BE can be treated as a BDF of order Q ¼ 1 (see Table 1), Equation (21) can be seen as a special case of Equation (26). Then, the adaptation conditions, that is, Equation (24) and Equation (25), required to remove the instantaneous reflections at both the electrical and the magnetic ports of the ME junction, are updated as It follows that linear multi-step discretization methods with variable step size could also be employed for the WD modeling of ME junctions by adopting an approach similar to the one applied in Bernardini et al. 21 to nonlinear circuits with linear capacitors and inductors.

| HIERARCHICAL SIM
As a reference circuit, let us consider the full-wave rectifier shown in Figure 5, where the transformer is modeled as in Section 2, and R M and L M represent the internal resistance and inductance of an electric motor. Applying the proposed multidomain modeling approach to the circuit in Figure 5 gives rise to the model shown in Figure 6. The magnetic domain is represented in red, whereas the electrical one is in black. We can notice the presence of two ME junctions, since the transformer has two windings. The rectangular one-ports represent, instead, the ladder networks used for modeling the eddy currents (see Figure 2).
The WD realization of the network in Figure 6 is shown in Figure 7. Considering that the reference circuit is composed of many nonlinear one-ports, such a WD structure is characterized by multiple implicit relations between the port variables, called delay-free-loops (DFLs), that cannot be avoided. In order to find the solution of the reference circuit at each sampling step, an iterative method is thus required.
In this Section, we propose a hierarchical version of the so-called SIM 12 suitable for solving nonlinear electromagnetic circuits like the one in Figure 6. In its original conception, SIM has been formalized to solve networks composed of one single topological junction and, for this reason, it cannot be directly employed to solve the WD structures T A B L E 1 Coefficients related to WD models of ME junctions based on BDFs with fixed step size Backward Euler Reference circuit showing a diode-based full wave rectifier. R M and L M represent an electric motor considered in this manuscript (see Figure 7). We thus provide a generalization of SIM, referred to as HSIM, able to solve WD structures consisting of an arbitrary number of ME junctions and WD topological junctions. The arrangement of the WD structure can be seen as a non-binary generalization of the WD Binary Connection Tree (BCT) discussed in F I G U R E 6 Diode-based rectifier employing a multidomain transformer model, where N T1 ¼ n lam n t1 , N T2 ¼ n lam n t2 , and the generic rectangular one-ports represent the different ladder networks shown in Figure 2. The electric (in black) and the magnetic (in red) domains are coupled by means of ME junctions F I G U R E 7 WD multidomain implementation of the diode-based rectifier. The blocks S lad1 ,…, S lad10 represent the different ladder networks, whose WD realization is shown in Figure 8. The magnetic blocks are depicted in red, whereas the electric blocks in black. The magnetic domain is coupled to the electric domain by means of the ME junctions. The T-shaped stubs indicate port adaptation Sarti and De Sanctis. 28 As a further difference, the BCT structure discussed in Sarti and De Sanctis 28 is characterized by only one non-adaptable nonlinear element, while in this manuscript, the WD structure contains multiple nonlinearities. The considered WD structure has many levels that go from level l ¼ 1 (the root of the tree) to level l ¼ L consisting only of leaves (one-port circuit elements). The largest WD topological junction, that is, the one with the biggest scattering matrix, is selected to be the root of the tree. This is done because the root will be the only non-adapted junction of the structure and the computational cost to make one port of a junction reflection-free is proportional to its size. Hence, the WD connection tree considered in this paper consists of a root (the biggest topological junction), a number of nodes (all the other topological junctions and ME junctions), and a number of leaves(the linear or nonlinear one-port elements). Further ahead, WD topological junctions will be also referred to as topological nodes, while ME junctions will be also referred to as ME nodes.
As far as notation is concerned, we use different subscripts to uniquely identify the leaves (one-ports) and the nodes (junctions) of the connection tree. Waves referred to leaves, as well as wave vectors and scattering matrices referred to nodes, do have a subscript in the form l, u to denote the uth leaf/node of the lth level. A variable referred to the nth port of the uth node in level l is indicated with the subscript l, u, n. Vectors with no subscript are a concatenation of the vectors of port variables referred to the whole circuit. It is important to stress that the waves incident to a certain level are the waves reflected by an adjacent level and vice versa.
HSIM is characterized by a forward-backward computational flow: during the forward scan, we compute the waves scattered from the leaves, passing through the nodes, and reaching the root; during the backward scan, we compute all the waves scattered by the root and reaching the leaves, passing through the nodes. In the forward paths from the leaves to the root, only one wave is reflected from the uth node in level l towards level l À 1. Conversely, in the backward paths from the root to the leaves, N l, u À 1 waves are propagated from the uth node in level l towards level l + 1, where N l, u is the number of ports of the uth WD block in level l.
The proposed HSIM algorithm comprises six stages that are performed at each sampling step k: 1. Initialization and Update: the port resistance Z l, u of each one-port element is set as close as possible to the tangent slope at the current working point on its v À i (or ℱ À ϕ) characteristic. For linear one-ports, this is straightforward. In fact, the optimal value of the free parameter is known a priori and can be set according to the adaptation conditions known in the traditional WDF theory. 13,35 For nonlinear elements, instead, it can only be estimated exploiting the solution of the circuit at previous sampling steps. 12,21,23 As far as topological nodes are concerned, we make the port facing the previous level reflection-free by means of the corresponding adaptation condition. For ME nodes, the F I G U R E 8 WD implementation of the ladder network shown in Figure 2. The waves a and b connect each ladder block to the main topological junction S m1 , as represented in Figure 7. The T-shaped stubs indicate port adaptation same action is carried out using Equation (24) or Equation (25). The matrices of the WD junctions need to be updated whenever their port resistances are changed. 2. Leaves Scattering Stage: waves b l, u [k] reflected by one-port elements (leaves) in level l are computed according to their WD scattering equations. In particular, for adapted linear one-ports we use Equation (9), while for nonlinear elements the reflected waves are computed employing a nonlinear mapping in the form (see Equation 11 and where γ is the HSIM iteration index. It is worth noticing that level L of the WD structure is entirely constituted of leaves, but other leaves can be present also in lower levels. where s l, u, n [k] is the row vector of the scattering matrix S l, u [k] corresponding to the nth port facing level l À 1. The nth entry of vector s l, u, n [k] is set to zero, in compliance with the junction adaptation condition for port n. On the other hand, the wave reflected by a ME node is computed using one of the two scalar scattering relations in Equation (21), depending on which port (electric or magnetic) is connected to level l À 1. whereas if the WD block is a ME node, b l, u [k] is scalar, and it is computed according to one of the two scattering relations in Equation (21), depending on which port (magnetic or electric) is connected to level l + 1.

Convergence Check: convergence is met when
where v (γ ) [k] is the vector collecting all the port voltages of the circuit at iteration γ, while ε HSIM is a small tolerance (e.g., ε HSIM ¼ 10 À5 ). Figure 9 shows a possible flow chart describing the HSIM algorithm at each sampling step k, where the WD structure is scanned following a level-by-level approach. As a consequence, the computational operations of the Leaves Scattering Stage and the Forward Scattering Stage can be interleaved, resulting in a more compact representation. Figure 10 shows, instead, the HSIM computational flow-graph related to the WD structure in Figure 7. As an example, let us consider the leaf-root path of diode D 1 . The computation starts at Level 4 and by means of the Leaves Scattering Stage (green arrows) we compute the wave reflected by diode D 1 and incident to node S e2 in Level 3. In the Forward Scattering Stage (orange arrows), the wave reflected by the port of the junction S e2 facing the electrical port of the junction M 1 /E 2 in Level 2 is computed. Then, the wave reflected by the magnetic port of the same ME junction and incident to the root S m1 is evaluated. Once all the other waves incident to the root are computed, the Root Scattering Stage (red arrows) can take place and all the waves incident to each leaf/node in Level 2 are assessed at a time. Finally, the waves incident to node S e2 in Level 3 and to diode D 1 in Level 4 are computed by means of the Backward Scattering Stage (blue arrows); in particular, all the waves reflected by the topological node S e2 and incident to Level 4 are computed at a time. The same analysis can be carried out for all the other leaf-root paths. As a final remark, it is clear from Figure 10 that a lot of operations are parallelizable (e.g., each subtree of the root can be solved in an independent fashion), and this can be exploited to enhance the efficiency of the HSIM algorithm.
F I G U R E 9 Flow-chart of the HSIM algorithm at each sampling step k; γ is the iteration index, U l is the total number of leaves/nodes in level l, while L is the total number of levels

| On the convergence of HSIM
In other studies, 12,21,23,25 it is shown that the convergence of SIM is guaranteed under the following hypotheses: (i) the topological junction is lossless and reciprocal; (ii) each circuit element is described by a continuously differentiable monotonically increasing v À i (or ℱ À ϕ) characteristic (hereafter indicated with v l, u (i l, u ) or ℱ l,u ðϕ l,u Þ); (iii) each port resistance is positive. However, it is not straightforward to extend the convergence analysis of SIM to HSIM. In particular, the WD structures considered in this manuscript are characterized by some nonreciprocal junctions (ME junctions), and in addition, some of their elements (the nonlinear reluctances) are modeled by means of CPWL functions that are not continuously differentiable 8 ϕ l,u ℝ (i.e., ℱ l,u ðϕ l,u Þ = 2 C 1 ). Nonetheless, in the following, we provide some considerations on the HSIM convergence properties.
HSIM convergence is reached when all the DFLs of the WD structure are iteratively solved. According to the setting of free parameters indicated in the "Initialization and Update" stage of HSIM, some ports of the topological junctions F I G U R E 1 0 Computational flow graph of the HSIM algorithm related to the WD structure shown in Figure 7. The circuit elements (leaves) are represented with green circles, topological and ME junctions with blue circles, whereas the root with a red circle. Incident/ reflected waves to/from a level are computed by means of the Leaves Scattering Stage (green arrows), Forward Scattering Stage (orange arrows), Root Scattering Stage (red arrows), and Backward Scattering Stage (blue arrows) and ME junctions are made reflection-free. Therefore, the only DFLs left to be solved are those involving all the pairs of nonlinear elements that cannot be adapted. For the SIM algorithm, it has been shown that by setting the free parameters Z l, u equal to the slope of the tangent at the current operating point on the element v À i characteristic, the number of iterations needed to reach convergence is minimized. 25,47 Although we cannot ensure HSIM convergence for each value of Z l, u , we can suppose that, by following the same approach, the speed of converge can be increased.
In other studies, 25,47 it has also been shown that, by expressing f l, u as a function of a l, u , the reflection coefficient can be rewritten as the derivative of f l, u with respect to a l, u , i.e., f 0 l,u ða l,u Þ, which depends also on Z l, u . It is possible then to minimize the instantaneous dependence of b l, u on a l, u by minimizing jf 0 l,u ða l,u Þj. As an example, when f l, u is linear, according to Equation (9), we have that f 0 l,u ða l,u Þ ¼ ðR gl,u À Z l,u Þ=ðR gl,u þ Z l,u Þ; therefore, jf 0 l,u ða l,u Þj can be easily minimized by setting Z l,u ¼ R gl,u . For nonlinear circuit elements characterized by v l, u (i l, u ) C 1 , instead, the minimization can be accomplished by setting the free parameter Z l, u equal to the tangent slope of the v À i characteristic at the current operating point. Unfortunately, when dealing with CPWL functions such an analysis cannot be directly applied because CPWL functions are not continuously differentiable, since their derivative exhibits discontinuities at the segment endpoints. However, we can similarly suppose that by setting Z l,u ¼ ℛ jl,u , where ℛ jl,u is the slope of the segment of the function ℱ l,u ðϕ l,u Þ evaluated at the current operating point, the instantaneous dependence between b l, u and a l, u is minimized. Finally, it is worth pointing out that, as already mentioned in Section 3.3.2, a CPWL function is characterized by a finite set of slopes. As a consequence, there exists a finite set of optimal free parameters Z l, u that can be precomputed. In Section 6, as an experimental validation of what discussed above, we show that setting the free parameters of the nonlinear elements as close as possible to the tangent slope at the current operating point on the v À i (or ℱ À ϕ) characteristic is fundamental for increasing the HSIM convergence speed.

| NUMERICAL RESULTS
In this section, we provide some numerical results concerning the HSIM algorithm and we carry out a performance comparison with two general-purpose solvers, that is, MathWorks Simscape 31 and LTspice. 9 We consider the circuit in Figure 5, which represents an AC/DC conversion stage with a transformer and a rectifier, widely used in industrial applications. Through simulations, we then investigate how harmonic distortion and losses are affected by the nonlinearities as core materials are varied. Table 2  in Equation (7), it is possible to compute the scattering matrix S lad related to the ladder network shown in Figure 2; Port 1 of this WD block is connected to the largest topological junction (root), whereas elements ℛ lad1 , ℒ 1 , ℛ lad2 , ℒ 2 , ℛ lad3 , ℒ 3 , and ℛ lad4 in this order are connected to ports ranging from 2 to 8. The scattering matrix S e2 of the junction connected to the transformer secondary side, instead, can be computed employing  Elements R F , D 1 , D 2 , D 3 , D 4 , C F , L M , and R M are connected to ports ranging from 2 to 9. Table 3 summarizes, instead, the port connections of the biggest topological junction S m1 (i.e., the root), whose fundamental loop matrix is : ð36Þ As input signal, we consider the sinusoidal line voltage of the electrical distribution system at 50 Hz. We thus set V in ½k ¼ Asinð2πkf 0 =f s Þ, where A ¼ 220V rms is the amplitude, k is the sampling index, f 0 ¼ 50 Hz is the fundamental frequency, and f s ¼ 96 kHz is the sampling frequency. Moreover, we take into account the typical harmonic spectrum of the distribution system that consists mainly of odd harmonics. In particular, we consider odd harmonics ranging from the 3 rd to the 15 th, whose magnitudes are taken from Cherian et al. 48 Figure 11 shows voltage v e1 and current i e1 at the transformer primary side and voltage v e2 and current i e2 at the transformer secondary side. The output current, that is, the current feeding the electrical motor, is, instead, shown in Figure 12. The blue and green curves representing the Simscape (SSC) and the LTspice outcomes, respectively, and the dashed orange curves representing the output of the WD implementation are overlapped. Referring to Figure 11A, for example, the maximum discrepancy between the WD and Simscape curves is around 0.03 V, which is negligible compared to the peak value of 253.4 V, demonstrating the accuracy of the proposed method. Figure 13 shows, instead, the norm of the error jjv (γ ) À v (γ À 1) jj 2 averaged over samples as a function of the number of HSIM iterations γ. We see that after four iterations, the mean error is already in the order of 10 À5 , that is, the order of the threshold ε HSIM (see Equation 33). This is compliant with the value of the mean number of iterations per sample required to reach convergence that we computed to be 3.7. As a further test on the HSIM speed of convergence, we computed the same mean number of iterations but without updating the free parameters, that is, ignoring what suggested in the "Initialization and Update" stage of the HSIM algorithm (see Section 5). We set the free parameters of WD nonlinear reluctances equal to the mean slopes of the segments of the CPWL functions, whereas the free parameters of WD diodes equal to one. The result is an average number iterations per sample equal to 169.6. It is thus evident that properly updating the free parameters during the simulation, by setting them as close as possible to the tangent slope at the current operating point on the v À i (or ℱ À ϕ) characteristics, is the key for generally increasing the HSIM convergence speed.
As far as the efficiency is concerned, we compared HSIM (implemented in a Matlab script) to both Simscape and LTspice, arranging the same circuit in the two environments by means of standard components, apart from the nonlinear reluctances that are implemented with CPWL functions using ad hoc blocks. Moreover, Backward Euler is chosen as discretization method. With the purpose of making a fair comparison with general-purpose simulators, we have not performed any optimization related to the particular circuit under consideration. Moreover, 100 identical  F I G U R E 1 1 Voltage (A) and current (B) at the primary side and voltage (C) and current (D) at the secondary side of the transformer. The dashed orange curve refers to the WD implementation, the blue curve refers to the Simscape (SSC) implementation, while the green curve refers to the LTspice implementation. In each plot, the three curves are overlapped, confirming the accuracy of the proposed method F I G U R E 1 2 Output current i out (i.e., the current feeding the electrical motor) vs. time. The dashed orange curve refers to the WD implementation, the blue curve refers to the Simscape (SSC) implementation, while the green curve refers to the LTspice implementation. The current waveform underlines the correct behavior of the rectifier, which, starting from a sinusoid, provides at steady state an almost constant signal to the load simulations of 0.35 s are performed (on a Intel i5 dual-core processor), and the mean simulation time (indicated with t sim ) is assessed. Table 4 shows the results of such a comparison, pointing out the effectiveness of the WD approach. The Simscape simulation is run considering the same fixed time-step of HSIM, whereas the LTspice simulation is run considering a variable step size, since it is the only possibility available in this environment. Nevertheless, even considering a maximum step size of 100 Á T s , LTspice results to be less efficient than HSIM. As regards Simscape, instead, it is known that it relies on efficient and optimized functions written in C++. 31 However, although the HSIM algorithm is currently implemented in a prototype MATLAB script, it has already shown to be faster with respect to the Simulink toolbox and, since HSIM is a generalization of the SIM algorithm, 12 we are confident that the speed difference will be even more pronounced as the size of the circuit under test becomes larger. Figure 14 shows a comparison between the magnitude of the input signal harmonics (represented in orange) and the magnitude of the harmonics of the voltage at the transformer secondary side (represented in blue). The total harmonic distortion (THD) of the input voltage is where V h and V 1 are the rms values of the hth harmonic and of the fundamental frequency, respectively. Voltage v e2 features, instead, a richer spectrum due to the harmonics introduced by the nonlinearity of the magnetic material. In this case, we have that THD % ¼ 20:4%. Thanks to the performed magnetic modeling, we are able to fully determine the magnetic flux distribution inside the transformer core. For example, Figure 15 shows the distribution of the magnetic flux ϕ inside the four areas of Ladder 1 (named "Lad 1 " in Figure 6). In the bottom right corner, the cross-section of the transformer is depicted: the thicknesses of the first three areas are chosen to be δ 1 < δ 2 < δ 3 , leading to Λ 1 < Λ 2 < Λ 3 , where Λ 1 , Λ 2 , and Λ 3 are the cross-sections of ℛ lad1 , ℛ lad2 , and ℛ lad3 , respectively (see Figure 2). The cross-section of the central parallelepiped, that is, of ℛ lad4 , results to be Λ 4 < Λ 1 . Given that the reluctance is inversely proportional to the cross-section, flux ϕ 4 (represented in red) flowing through ℛ lad4 is the smallest and flux ϕ 3 (represented in purple) flowing through ℛ lad3 is the largest, as we can appreciate looking at Figure 15.
Comparison between the magnitude of the harmonics concerning input signal V in (represented in orange) and the magnitude of the harmonics concerning voltage v e2 at the transformer secondary side (represented in blue). The input signal is characterized by odd harmonics ranging from the 3rd to the 15th and by a THD % ¼ 8:25%. Voltage v e2 presents, instead, a THD % ¼ 20:4% F I G U R E 1 5 Magnetic flux distribution over time among the nonlinear reluctances of Ladder 1. With reference to Figure 2, flux ϕ 1 (represented in green) flows through ℛ lad1 , flux ϕ 2 (represented in blue) flows through ℛ lad2 , flux ϕ 3 (represented in purple) flows through ℛ lad3 , and flux ϕ 4 (represented in red) flows through ℛ lad4 . In the bottom right corner, it is shown the cross-section of the transformer and the subdivisions into "concentric" areas: the smallest area is associated to ℛ lad4 , whereas the biggest area to ℛ lad3 The speed up provided by HSIM allows the designer to run several simulations in a reduced amount of time w.r.t. the mentioned mainstream MNA-based methods. HSIM reveals thus to be a valid and promising candidate for parametric analyses requiring repeated simulations, since it allows the designer to tackle the choice of specific circuit parameters in less span of time. In this regard, we run four simulations varying the material of the magnetic core. Besides isovac, we considered other three materials with different values of resistivity and saturation induction fields: iron, M15 (steel), and Hiperco50 (soft magnetic alloy). Table 5 summarizes the results of such an analysis, showing the THD% at the transformer secondary side and the losses for all the materials taken into account.

| CONCLUSIONS
The flexible WD methodology proposed in this manuscript can be employed to simulate nonlinear electromagnetic circuits in a multidomain fashion. We have shown how to implement a transformer with an arbitrary number of laminations and an EI geometry taking into account flux leakage and eddy currents. However, the described method can be applied to whatever MEC model of transformers or inductors. We have provided a general WD description of a two-port junction, referred to as magnetic/electric (ME) junction, that couples magnetic domain and electric domain, converting magnetic variables into electric variables and vice versa. The method allows us to locally handle nonlinear elements, enabling the possibility of choosing the most suitable modeling technique for each nonlinearity. In this work, we have shown how to model diodes in the WD domain by means of the Wright ω function and we have proposed a methodology, based on CPWL representations, to model nonlinear reluctances. The choice of CPWL models allows us to employ specific B À H characteristics, contrary to common models that approximate a generic magnetic nonlinearity. The presented HSIM algorithm is a generalization of the SIM 12,21 and turns out to be characterized by a high degree of modularity, being at the same time accurate, robust, and efficient. We have compared its performance to the one of MathWorks Simscape and LTspice obtaining promising results and we expect that such a performance gain would be more pronounced dealing with larger nonlinear circuits. 12 The results discussed in this manuscript encourage the development of new general-purpose simulators based on WD principles that take advantage of the latest findings in the WDF literature. 21,47 As far as future works are concerned, it is worth noticing that most of the HSIM operations can be easily parallelized, since each nonlinearity is locally handled; therefore, the investigation of a multithreaded implementation of HSIM is promising. Another possible development is to model reluctances with alternative piecewise representation of nonlinear functions that are differentiable. 49 A further interesting extension of the proposed approach would be searching for strategies to effectively accommodate nonlinear reluctances with hysteresis. Moreover, it would be interesting to test HSIM for the solution of electromagnetic circuits modeled using the PEEC method. 4 Finally, it is worth extending the proposed multidomain simulation methodology to other physical domains in addition to the electric and the magnetic ones, aiming at implementing more complex multiphysics systems.

DATA AVAILABILITY STATEMENT
Data sharing is not applicable to this article as no datasets were generated or analyzed during the current study.