Simulating circuits with impasse points

In this paper circuits with impasse points, i.e. with jumps in their configuration space will be analyzed. These non-regularized circuits exhibit a fold in their configuration space, which can lead to difficulties during the simulation with standard circuit simulators like SPICE. The former developed geometric approach to simulate these circuits without regularization will be extended by a detailed discussion of which coordinate system has to be chosen. Furthermore, two new approaches for a numerically efficient calculation of the hit points will be shown.


Introduction
There is a class of electronic circuits whose configuration space manifold S is folded within the embedding space E of currents and voltages.This fold can be related to so-called impasse points and leads under certain conditions to jumps from one stable part of S to another.A classical transient solution of such non-regularized circuits exhibiting impasse points is not possible (see e.g.Chua (1980), Reissig (1996), Chua and Deng (1989a)).However, a common method to overcome these simulation problems is to regularize the circuit by adding suitable located parasitic inductors L's or capacitors C's considering Tikhonov's Theorem (for further literature see Reissig (1996)).In previous works, a geometric concept was developed to simulate those circuits without regularization (see e.g.Thiessen and Mathis (2011a), Thiessen et al. (2013)).There, several electronic circuits exhibiting these behavior were studied and different approaches to calculate S, jump and hit points were analyzed (e.g.Thiessen et al. (2012b), Thiessen et al. (2012a), Thiessen et al. (2011), Thiessen and Mathis (2011b), Thiessen and Mathis (2011a)).In this work improvements in the topic of hit point calculation will be shown.Especially the difficulty of multiple hit points will be studied in detail.Furthermore, a detailed description of the geometric system analysis will be given.There, the question of which coordinate system, i.e. system of equations, is best placed to deal with those non generic circuits will be answered.

Geometric system analysis
Common circuit simulators (e.g.SPICE) are based on MNA, which leads in the description of electronic circuits to a quasilinear differential-algebraic system of equations (DAE) (cf.Riaza (2008)): This system of equations can also be formulated as In the sense of Reich Reich (1990) eq. ( 2) is a triple (R k , B, h), where B : R k → R k×k is a linear operator and h : R k → R k a diffeomorphism.Reich had proven that the DAE (2) is regular, if there is a differentiable manifold S ⊂ R k and a vector field v : S → T S, such that a differential mapping w : I → S (I ⊂ R) is a solution of the vector field for all t ∈ I , if and only if the mapping c : j • w : I → R k is a solution of the DAE, where j : S → R k is the natural injection Reich (1990).
This work focuses on the investigation of non-regular DAEs where impasse points exists.These impasse points will be called "jump points", because the transients will be continued by jumps in E from one point on S to another.
The systems of equations of the considered nonlinear dynamical circuits can be characterized by a semi explicit DAE: The vector x ∈ R n corresponds to the capacitor voltages and inductor currents and the vector y ∈ R m to additional voltages and currents in an electronic circuit.
Taken into account that S can exhibit a fold e.g.respectively an input voltage v s (t) and not respectively t, the input sources have to be treated differently.For describing the behavior of the circuit for any input values, the independent and time dependent input sources will be replaced by norators and treated as further vector z fo unknowns (cf.Fig. 1).Therefore, an additional vector z ∈ R η has to be considered in the describing system of equations.The resulting autonomous, semi explicit system of equations then is described by: These system of Eq. ( 4) forms the basis of the further investigations.The distinction in x, y and z is essential for calculating S and the jump and hit points.S can be defined as a subspace of E = R k (where k = n+m+η) and is represented by the solution set of the independent algebraic Eq. (4b).The calculation of the jump and hit points will be shown in Section 3 and 4.
Remark: As can be seen, the system of equations (1) do not distinct between x, y and z.One method to modify the system of equations (1) to (4) is shown in Thiessen et al. (2012a) and Thiessen et al. (2013).Another method for deriving a system of equations of type ( 4) is by applying the Augmented Nodal Analysis (ANA) shown in Riaza (2008).The resulting system of equations derived by the Augmented Nodal Analysis (ANA) is Here, the vectors v c and i l are counted among the vector x.
The vector y is composed of i c , i u and e and the input vectors i s and v s are counted among the vector z.

Jump points
It was shown by Chua and Deng (1989b), that almost all singular points of an autonomous semi explicit DAE are in fact impasse points, i.e. jump points.However, e.g. in Chua and Deng (1989a), Chua and Deng (1989b) and Andronov et al. (1966), the fact that S could contain a fold respectively z was neglected.Taken this into account and knowing that singular points are points, where the local solvability to y is not guaranteed (e.g.Thiessen et al. (2011), Chua and Deng (1989a)), the necessary jump condition can be specified by (cf.Chua and Deng (1989a), Andronov et al. (1966)): A point that is specified by eq. ( 6) and whose neighborhood includes each a Lyapunov-stable and -unstable point, is defined as proper jump point P j .This sufficient jump condition can be verified by calculating the eigenvalues λ i of the characteristic equation det ∂ y f(x, y, z) − λ • E = 0, where E is the identity matrix (cf.theory of discontinuous oscillators e.g.Andronov et al. (1966), Mishchenko and Rozov (1980)).
The set of all points fulfilling these two conditions is called jump-set , which represents a l − 1-dimensional subset of S. Of course, the calculation of the zero set of all points fulfilling the m + 1 algebraic equations specified by Eq. ( 6) is difficult.However, not all zeros of Eq. ( 6) are of interest, but only in the actual chosen jump point during a simulation.Hence, the dynamics on S will be traced till reaching a stopping point P s .This stopping point is defined as a point, where the step size of the numerical solver reaches a lower boundary (which is related to the machine constant of the simulating computer).In the next step, the "nearest" point on will be calculated by choosing a suitable norm and defined as the actual jump point P j (cf.Fig. 2 (a)).

Hit points
The dynamics on S is specified by Eq. ( 3a) and (3b).Since the dynamics is not defined at points on , the trajectory increases very fast nearby while tracing it.Therefore the numerical integration is stopped at P s when the step size reaches the lower boundary.From there, the jump point P j is calculated as described in Sect.3. To trace the dynamics, we use a variable order solver based on the numerical differentiation formulas (NDFs) Shampine and Reichelt (1997).
Considering that the voltage across a capacitance and the current through an inductance is preserved, both have inertia through a jump process and do not change (i.e.x j = x h ).Assuming that a jump happens instantaneous, i.e. t j = t h , a further restriction is the fixed value of z (i.e.v s (t j ) = v s (t h ), i s (t j ) = i s (t h )) during a jump.Consequently, a jump takes place in a tangential space of R m , which corresponds to the coordinate space of y.In the following, the jump space will be denoted by J S. This corresponds to the jump postulate of Chua and Alexander (1971).
Because we introduced E, a hit point P h can be calculated by the intersection of J S defined in P j and S excluding the jump point itself (P h ∈ (J S P j ∩ S) ) (see Fig. 2 (b)).Thus, the problem of the hit point calculation can be defined as with the constraint In Thiessen and Mathis (2011a), Sarangapani et al. and Thiessen et al. (2012a) the hit point calculation was done in two steps: First a point P j outside S was chosen so that P h ∈ (J S P j ∩ S).In the next step, the actual hit point P h ∈ (J S P j ∩ S) with the initial condition P h was calculated.There, the difficulty was to chose a suitable point P j outside S, so that the numerical solver is able to find P h .
Another approach (cf.Thiessen et al. (2011)) to calculate a hit point was to use a bisection method.This provided a set of possible points from which the closest to the corresponding jump point was chosen.Now, an efficient approach to calculate the hit points of circuits with only one fold in S based on the penalty function will be shown.

Penalty function
The basic idea of a penalty function is to convert the zero problem of Eq. ( 7) with a constraint (8) in a zero problem without constraint (cf.Florian Jarre (2003)).
Thus, the optimization problem h(y) → min with y = y j , where h : R m → R m can be transformed to The penalty function p(y, r) consists of the weighted sum of the objective function h(y) and the penalty function l(y).The penalty function itself can be weighted by the vector r = (r 1 , r 2 , . . ., r m ) T which is chosen to be r = r •(1, 1, . . ., 1) T in this work.
Here, an easy penalty function for the constraint of eq. ( 8) was chosen: It follows that l(y) → 0 for y = y j and l(y) → ∞ for y ≈ y j .
To achieve a more robust solution, suitable initial conditions y h,0 outside S has to be chosen.It became apparent, that the addition of the inverse of the last step y to y j , yields a more robust numerical solution than by choosing e.g.y h,0 = 0.This can be explained by the switching process of the considered circuits.
From P h the dynamics can be traced, till reaching again.

Multiple hit points
Another problem appears if the configuration space S of the electronic circuit is multiple folded, so that there are multiple possible hit points (cf.Fig. 3).In these cases, the penalty function approach is not suitable and a new approach based on the homotopy method is needed.

Homotopy method
The homotopy method is used to solve nonlinear algebraic equations of the form (7). The advantage is that the convergence region is much larger than the one by applying the penalty approach.Starting from an easy zero problem, the system of equations will be deformed by λ till reaching the original zero problem.This continuous deformation process is achieved by solving H(y, λ) = 0.A general homotopy can be given as follows: This homotopy consists of a linear combination of two real functions: h(y), whose zeros are sought and B(y), a function for which a zero is known.The difficulty is to choose the proper function B(y) for the corresponding zero problem Trajkovic and Mathis (1995).
T. Thiessen and W. Mathis: Simulating circuits with impasse points The zero curve of the homotopy map H(y, λ) (homotopy path) can be tracked by different techniques.One technique, which is used in this work, is to use an ODE-based algorithm (cf.Watson (1990)).For this, the parametrization is not based on λ but on the arc length s of the homotopy path.Therefor, the equation H(y(s), λ(s)) = 0 has to be differentiated with respect to s and solved for λ and y.The main advantage for using the arc length method is because therewith regressive homotopy paths can be traced.A disadvantage of such a multiple-step method is, that the error of the approximation of y increases with every step.By using a predictor corrector method, the error accumulation can be counteracted.
In the following the problem of multiple hit points will be displayed by the example of two series connected resonance tunneling diodes (cf.Fig. 4 (a)) Thiessen et al. (2012b).The chosen V -I characteristics of the tunnel diodes is shown in Fig. 4 (b).A good approximation of these V -I characteristics can be achieved by using the following equation (cf.Chang et al. (1993)): where c i are constants related to the peak and valley currents and voltages, e is a scaling factor and m and l are integer fitting factors.By using a norator, as explained in Section 2, the configuration space S of the series connection can be determined by solving the system of equations It is noteworthy that the configuration space of the series connection consists of two separated manifolds: one main part proceeding through the origin and a second separated single loop which can only be reached by choosing suitable initial conditions (cf.Thiessen et al. (2012b)).If the current I is increased from zero, the first jump point of the nonregularized circuit of Fig. 4 (a) will be reached at the point P j .From there, there are five possible hit points as shown in Fig. 5).
For the sake of completeness, the transients of the corresponding regularized circuit are also shown in Fig. 5) by the green lines with arrows.The regularized circuit can be achieved by adding small capacitances parallel to the diodes (cf.Thiessen et al. (2012b)).To calculate the hit points, the FPN homotopy shown in Rahimian et al. (2011) is used.The FPN homotopy consists of a combination of the fixed-point (FP) and Newton (N) function approach and can be formulated in two steps: (1) The function h(y) to be solved is multiplied by the fixedpoint function giving (2) The function B(y) is formed by a combination of the fixed-point and Newton function (cf.Rahimian et al. (2011)) yielding Inserting Eq. ( 14) and eq.( 15) in Eq. ( 12) gives After some rearrangements, eq.( 16) simplifies to With this homotopy, the corresponding homotopy path for calculating the hit points of the series connection can be seen in Fig. 6.The jump point P j was chosen as starting point y 0 , which is also a solution of h(y).To calculate further zeros except P j , the homotopy path was progressed beyond λ start = 1.Therefore, a bifurcation point (BP) was inserted at P j (cf.Dhooge et al. (2006)).
The numerical calculation of the homotopy path was done by a numerical continuation method, which includes a predictor corrector method.Therefore, the MATLAB toolbox CL MATCONT was used Dhooge et al. (2006).
As can be seen from Fig. 5 and Fig. 6 not all hit points can be calculated with this method.The hit points on the separated single loop cannot be calculated, but all hit points on the main part of S. This results due to the fact, that S consists of two separated manifolds.
The hit points on the separated single loop are in fact instable points and could be calculated by choosing suitable initial conditions.But, for tracking the transients, only the stable hit points are of interest, which are P h2 and P h4 .Now the question is, which of these both points is the right hit point ?
In the work of Chua and Alexander (1971) there is an inertia postulate which says: trajectories on a stable branch will continue on a stable part till reaching a jump point and than jump to the "nearest" stable part of S. But a complete proof of this postulate is still pending.One possibility to verify this postulate is to analyze the catchment area of the dynamic of the corresponding regularized circuit, but this will be the focus on further studies.

Conclusions
In this work, a detailed description of the geometric system analysis was given.Thereby it was explained why the system of equations yielding from the ANA is the most suitable for applying the geometric approach.Furthermore, the penalty and the homotopy approach for a numerically efficient calculation of the hit points were shown.The penalty function method turned out to be suitable for circuits with only one fold in S, but not for several folded configuration spaces.By using a homotopy method with a proper homotopy function, all hit points can be calculated assuming S consists only of one manifold.In this work a non generic case, where S consists of two separated manifolds was shown.In those cases, not all hit points can be calculated without further arrangements.Furthermore, the question of choosing the right hit points appeared, but shall be studied in further investigations.

Fig. 2 .
Fig. 2. Calculation principle of (a) P j and (b) P h

Fig. 5 .V
Fig. 5. S of the series connection of D1 and D2