Speed-up hyperspheres homotopic path tracking algorithm for PWL circuits simulations

In the present work, we introduce an improved version of the hyperspheres path tracking method adapted for piecewise linear (PWL) circuits. This enhanced version takes advantage of the PWL characteristics from the homotopic curve, achieving faster path tracking and improving the performance of the homotopy continuation method (HCM). Faster computing time allows the study of complex circuits with higher complexity; the proposed method also decrease, significantly, the probability of having a diverging problem when using the Newton–Raphson method because it is applied just twice per linear region on the homotopic path. Equilibrium equations of the studied circuits are obtained applying the modified nodal analysis; this method allows to propose an algorithm for nonlinear circuit analysis. Besides, a starting point criteria is proposed to obtain better performance of the HCM and a technique for avoiding the reversion phenomenon is also proposed. To prove the efficiency of the path tracking method, several cases study with bipolar (BJT) and CMOS transistors are provided. Simulation results show that the proposed approach can be up to twelve times faster than the original path tracking method and also helps to avoid several reversion cases that appears when original hyperspheres path tracking scheme was employed.

a method capable to avoid the reversion phenomenon (Yamamura 1993b). Additionally, a selection criteria for the path tracking starting point is addressed. This paper is organized as follows. "Original HCM scheme for studying PWL circuits" section provides a brief description on PWL modelling, a short introduction to the HCM method, and a summary of the proposed method in Vazquez-Leal et al. (2014). The suggested technique for avoiding the reversion phenomenon, the proposed path tracking method, and selection criteria for the starting point are provided in "Proposed homotopy scheme" section. In "Cases study" section, five cases study of nonlinear circuits are presented and solved using an HCM and the proposed path tracking method. Numerical simulations and discussion about results are provided in "Numerical simulation and discussion" section. Finally, our conclusions about this work are given in "Conclusions" section.

Original HCM scheme for studying PWL circuits
Equilibrium system of equations The equilibrium system of equations are obtained applying the modified nodal analysis (MNA) (Ho et al. 1975), which is a method that allows the systematic study of circuits containing devices incompatible with the classic nodal analysis like voltage sources, voltage-dependent voltage sources, among others.
As a result of the MNA a set of equations of the form will be obtained; where x represents the vector of unknowns (electrical variables) of the circuit.

Piecewise linear (PWL) model
The Piecewise Linear Model is an approximation of a nonlinear equation to a set of linear equations which, altogether, exhibit the same behavior as the original system. In this work, Chua's model serves as base for the proposed homotopy scheme; this model is described as follows model parameters are computed by where β represents the breakpoints, represents the number of breakpoints, and J i represents the slope of the i-th straight line segment in the PWL model.

Homotopy continuation method (HCM)
To solve a system of equations using the HCM, first, the actual solution is introduced in a set of solutions described by (1) where is the homotopic parameter and f (x) is the system of equations to be solved. To find the solution of the original system, the HCM starts from a known solution of the homotopic system, which is commonly given at = 0. Afterwards, a path tracking method is employed to calculate subsequent points within the homotopic curve. Each time the homotopic path intercepts the solution line (generally placed at = 1) a solution to the original system is found (Vazquez-Leal et al. 2011a).
The HCM is commonly used to find multiple solutions as it does not stop the calculation of solutions once it has found one, unlike the NRM which is designed to find just one solution per simulation.
An important issue of the HCM is the possibility that path tracking fails by following an incorrect path, this will cause losing solutions or not finding any solution at all, even if the Homotopy path exist and several, or all, solutions may be located.

Homotopy formulation
The homotopy formulation used in this work is Newton's homotopy given by As reported in Vazquez-Leal et al. (2014), it is possible to model devices using PWL techniques by applying Newton's Homotopy during DC analysis. Nevertheless, it is important to notice that experiments in Vazquez-Leal et al. (2014) proved that homotopic curves are also PWL as long as Newton's Homotopy is employed. This important characteristic will be applied in this work to propose a novel initial point selection scheme and a new scheme to accelerate the trace for the homotopic path.

Modified spheres algorithm (MSA)
In Vazquez-Leal et al. (2014) the modified spheres algorithm path tracking method was used to find the homotopic curve. This method consists in including the equation of a sphere within the original homotopic system (Vazquez-Leal et al. 2011b;Torres-Munoz et al. 2014;Yamamura 1993b;Oliveros-Munoz and Jimenez-Islas 2013;Jimenez-Islas 1996). It is expressed as where c is the center of the sphere, r is its radius, and n is the number of variables from the equilibrium system of equations. By incorporating (6) into the homotopic system (4) the system of equations become H n (f n (x), ) = 0, S(x 1 , x 2 , · · · , x n , ) = 0, the system contains n + 1 equations and n + 1 variables. Figure 1 shows the application of the MSA. When the center of the sphere is located at c 1 , the NRM is applied with predictor vector k 1 ; corrector steps are applied to achieve the intersection between the sphere and the homotopic path which will be used as center of the next sphere S 2 .

Starting point criteria
Chua's model possess the characteristic of having upper and lower limits, this means that its linear behavior will not suffer any change because there are no breakpoints affecting its linearity. In addition, when using an HCM, the starting point is selected far from the interest area as a means to induce the homotopy to run through an entire region where a solution may be located, this is known as feasible region. Therefore, two tests are proposed for the selection of the starting point.
• Test 1 Every value of the starting point vector must not have a value located within the feasible region. In electrical terms, the values for voltage variables must be higher (or equal) than the highest positive voltage supply or lower (or equal) than the lowest negative voltage supply; the values for current variables should be in a higher range than the possible normal operating currents for the circuit, for instance, flipflops typically work in the range of milliamperes so setting the values in the range of amperes would be suitable for current variables. • Test 2 This work proposes circuits modeled by devices of type i = y(u). Therefore, the test will focus on this kind of elements. Nevertheless, a simple extrapolation of the explanation in this section may be extended to elements of type u = y(i). As it has been explained in Test 1, it is important that initial point should be located above or below the maximum and minimum values of the power supply, respectively. By doing this, the chance of the homotopic path to cross all the feasible region of solutions is increased. However, it is important to understand that initial point consist in a set of electrical values (nodal currents and nodal voltages). Therefore, given the nature of the MNA formulation, most of the electrical variables for the circuits under study are nodal voltages, being the only variable of type current the unknown current from the power source. Besides, we know that every PWL device has an specific number of breakpoints ( ) and within them are breakpoints that may cause multiple operating points. Therefore, it is important to assure that proposed initial point in Homotopy trajectory terms of nodal voltages (v k and v n for each pair of terminals of the PWL devices), produce a voltage drop u outside the bounded region by the lower breakpoint (B L ) and upper breakpoint (B U ) as it can be seen in Fig. 2. The shaded region is the feasible region for solutions although solutions, in fact, could be located out of this region for some devices. Therefore, Test 1 and Test 2 are complementary; helping to propose initial points capable to provide the most number of solutions for simulation or path tracking.

Avoiding the reversion phenomenon
The MSA path tracking method consists in calculating the point where the circumference of a sphere crosses the homotopic path; nevertheless, this sphere always intersects the homotopic path in two points, as for our purposes we are only interested in just one. It is possible that the NRM calculations for the new point in the path converge to a point already found, this situation cause a backward path tracking, thus causing the method to fail. This situation is known as the reversion phenomenon (Yamamura 1993b). This work introduces a technique capable to avoid the reversion phenomenon. It consists in perturbing the hypersphere equation (6) to avoid one of the interceptions between the hypersphere and the homotopic path.
To understand the proposed technique, it is necessary to study the concept of an inverted sphere which is described by where x n+1 represents the homotopic variable ( ), c inv represents the inverted sphere center and r inv represents its radius. This equation is the same as obtaining the square root of (6). When substituting any value of x located inside the sphere in (8), it will generate a negative number inside the square root and create an empty region in the domain of real numbers as shown in Fig. 3. In order to take advantage of the properties from the inverted sphere we add (6) and (8), placing the center of the inverted sphere at the unwanted interception and the center of the original sphere at the point obtained in the previous iteration, the result is here j represents the j-th iteration and K inv is an arbitrary constant used to reduce the contribution of the inverted sphere to the equation in order to deform the shape of the sphere as little as possible.
As shown in Fig. 4, the resulting system has a cavity in its circumference. In fact, this empty spot covers the unwanted interception of the original sphere and the homotopic path. This factor allows the NRM to find the next point in the homotopic path as it is the only solution for which is the new system of equations to trace the homotopic path.
Homotopic path

Fig. 4 MSA with an inverted hypersphere
In order to avoid a possible oscillation of the NRM or iterations with complex numbers, when it approaches the empty region, a limit in the number of iterations is established. When the limit is reached, the NRM is executed with a different starting point as described in Torres-Munoz et al. (2014). Empirically, we found that the radius of the inverted sphere should be smaller than the radius of the original one, a range between 1,000 to 10,000 times smaller is proposed in this work. If the inverted sphere radius is too small, the NRM might find the unwanted root even if it does not exist, because the radius from the inverted sphere is smaller than the error tolerance from the NRM as seen in Fig. 5a. If it is too big, the homotopic path tracking may fail because next point in the path could be inside the empty region of the original sphere as shown in Fig. 5b.

Speed-up hyperspheres path tracking method (SHPT)
This work proposes a modification to the path tracking method presented in Vazquez-Leal et al. (2014). It is capable to reduce the computing time for tracking homotopic paths having PWL characteristics. Taking advantage of the local linearity of PWL models, a parameterized straight line equation is deduced from the first two points (x 0 and x 1 ) obtained using the MSA (see Fig. 6). After that, this linear equation is used to calculate the next point (x 2 ) in the straight line and the values obtained are substituted in the homotopic system of equations, if these values satisfy the system of equations, the next iteration will be performed in the same way to find (x 3 , x 4 , . . .) as depicted in Fig. 7a; otherwise, it means that we found a break point and the MSA must be used again two times in order to obtain a new straight line equation as depicted in Fig. 7b.

Homotopy formulation
The Newton Homotopy should be formulated based on (1). As shown in Vazquez-Leal et al. (2014), the homotopic curves obtained with this formulation have a PWL nature, this will hold as long as all of the nonlinear devices in the circuit are PWL modelled. This property allows to execute the following steps for the proposed path tracking method.

Starting point criteria
Once the homotopic system is defined, the SP should be arranged in a way to accomplish both Test 1 and Test 2 (see Fig. 2).

Modified hypersphere equation (hypersphere iterations)
Once the starting point (x 0 ) is defined, the sphere equation (6) should be formulated with a center located at c = x 0 . With this formulation the NRM must be applied to (10).
If the found solution is within the region of < 0 (Fig. 6a), it will be used as the center of the inverted sphere for the next iteration and the center of the normal sphere will be located at x 0 in order to induce the NRM to find the interception on the positive side of and follow the path heading to = 1. If the found solution is within the region of > 0 (Fig. 6b), the next sphere will have a center located at x 1 and the inverted sphere will have a center at x 0 .

Straight line equation formulation
Once the direction of the path tracking has been set to the positive region of , it is possible to formulate the straight line equation given by where j represents the j-th iteration and m represents the slope of the straight line that crosses from x j to x j−1 . Next iterations will be predicted substituting j+1 = 2 j − j−1 in (11), as shown in Fig. 7. If the values calculated using (11) satisfy (7) it can be assured that they belong to the homotopic path. This procedure shall be repeated for next predictions until the obtained values no longer satisfies (7). It means that iterations with straight lines found a breakpoint on the homotopic path just like the one in Fig. 7b. Therefore, two new iterations using the hypersphere tracing point should be performed to create a new straight line and predict points for the new segment.
The straight line path tracking does not need any correcting steps like those needed when using the NRM, this results in greatly reducing the computing resources and time required. It also avoids the calculation of different starting points when NRM fails. Furthermore, the diverging issues present sometimes in the NRM are moderated for these iterations. Finally, no reversion phenomenon will appear as (11) has only one solution and it leads to the forward path tracking.

Path tracking technique near breakpoints
When the homotopic path crosses a breakpoint of the PWL model (see Fig. 7b), (11) will not satisfy (7). At this point, (10) has to be solved, placing the center of the noninverted sphere at the last calculated point that was part of the homotopic path and the center of the inverted sphere at the second to the last calculated point; NRM is applied to calculate the next point in the path. Afterwards, another point should be calculated by solving (10) in order to repeat the hypersphere iterations procedure explained in "Straight line equation formulation" section. When a path has a high density of straight line segments, the SHPT will tend to slow down, although not as slow like the method proposed in Vazquez-Leal et al. (2014). This characteristic shows that SHPT, compared to MSA, requires lower computation time (CP) or could perform almost identical if the homotopic curve exhibits a high number of break points.

Find zero strategy
As reported in Vazquez-Leal et al. (2014), if the homotopic path crosses the solution line ( = 1), the exact solution to the PWL system can be obtained by calculating (11) using the last point before the homotopic path crossed the solution line and the next point after it, as shown in Fig. 8, and substituting j+1 = 1. In other words, a linear interpolation at = 1 is performed among the two iterations crossing = 1.

SHPT algorithm
In this section we introduce an algorithm for the SHPT method. (11) The general procedure work as follows: first, the system of equations is generated as shown in (14). Second, the user provides a starting point. This starting point has to fulfill Test 1 and Test 2 already introduced in the previous section. For the case when starting point does not fulfill any of the tests, the user is requested to provide another starting point and the verification process is performed until a valid point is achieved. Then, the process continues by generating the Homotopic formulation for the system. Afterwards, the hypersphere formulation is given. It is important to mention that, at this point, the inverse hypersphere formulation is also generated to avoid the reversion phenomenon. Notice that the hypersphere and inverse hypersphere formulation will be updated after every iteration by adjusting its center. Once all necessary equations are already provided and a valid starting point is given, the main loop is started. It begins with a hypersphere iteration, the result of the iteration and the starting point allows the calculation of another point that allows calculation of a predictive straight line in the next block. Once this step is done, the Straight Line Iteration block is performed as explained in the corresponding section. As iteration continues, detection for crossing at = 1 is performed; for the case that it is not detected, break point detection is applied. If a break point is not found, the process returns to the Straight Line Iteration block. For the case that a break point is detected, two hypersphere iterations are performed to correct the homotopic path (see Fig. 7). When = 1 is located, linear interpolation (see Fig. 8) is applied to provide a very accurate approximation that is stored in a solution vector. Once the store block is executed, the process returns to the straight line iteration block. The algorithm will repeat until maximum number of iterations limit is reached. Figure 9 shows the flow diagram for the entire process.

Cases study
In this section we provide five circuits to be analysed using the proposed path tracking method. All the proposed SP for the circuits fulfils the Test 1 and Test 2 starting point criteria. Also, we compare the proposed path tacking method against the proposed in Vazquez- , this comparison includes the results using the inverted sphere and without it. The value for K inv was set to 50000 for all the cases.

Example 1
The circuit in Fig. 10a (Tadeusiewicz and Kuczynski 2013) has two BJT Transistors which were modelled using the Ebers-Moll model shown in Fig. 10b. The PWL equations for these devices are   Figure 10c shows the equivalent circuit once the Ebers-Moll has been substituted.
Applying the MNA method to the circuit, the equilibrium system of equations are obtained as which is used as the base to create a Newton homotopy, and results in the following homotopic system of equations The proposed starting point for the homotopic path tracking is shown in Table 1. Table 2 shows that the SP fulfils both Test 1 and Test 2 criteria. The first column contains the variables of the NAES and the PWL model dependent variables u = v k − v n , where k and n represent the diodes nodes; the second column shows that every variable is greater than or equal to the value of the voltage source (V S1 = 5), accomplishing Test 1; the third column shows the fulfilment of Test 2 as every PWL model dependent variable u = v k − v n is less than or equal to the lowest breakpoint (B L = 0) in (12).
The resulting homotopic path when using the proposed SP is depicted in Fig. 11 for v 4 and v 2 .
In Tadeusiewicz and Kuczynski (2013) three operating points were reported for this circuit, in Fig. 11 can be seen that all of them were found using the proposed methodology. The solutions found are listed in Table 3.
The errors of the obtained solutions are calculated by where N represents the number of equations in the original system and f (S i ) represents the substitution of the solution S i in the vector of equations (1).

Example 2
The circuit shown in Fig. 12 is the classical Chua's circuit, this circuit has nine solutions. The values of the resistors are taken from Reyes (1994) and BJT Transistors were modelled using a simplified version of the Ebers-Moll model, see Fig. 13. The V-I characteristics of the model in this model are given by The equilibrium system of equations obtained using the MNA has sixteen variables, thirteen are voltage variables and three are current variables from the three independent voltage sources.  In order to be able to calculate all the solutions of this circuit three SP were needed. These SPs are listed in Table 4 and all of them fulfil the starting point selection criteria proposed in this work.
The resulting homotopic paths for variable v 5 are shown in Fig. 14    solutions because S 5 was found using the SP3 and SP2, and S 9 was found using SP1 and SP3. Numerical solutions are listed in Tables 5 and 6.

Example 3
This circuit has been proposed for this work (see Fig. 15). Contains three flip-flops connected in a cascade configuration. The transistors are modelled using the Ebers-Moll model. The V-I characteristics for diodes D 1 , D 2 , D 3 , and D 4 are given by where v D represents the voltage drop in diode D and i D represents its current. This circuit contains twenty variables; calculating one SP (listed in Table 7) found five different operating points. Table 8 shows the operating points and Fig. 16 presents the homotopy path for v 12 .

Example 4
This circuit (see Fig. 17) was studied in Tadeusiewicz and Kuczynski (2013); it contains five NMOS transistors and five PMOS transistors. The CMOS transistors are represented by using the Shichmann-Hodges model (Shichman and Hodges 1968) (Fig. 18) and simulated in SPICE setting the following parameters: LEVEL = 1, VT0 = 0.5705, RD = RS = 0, LAMBDA = 0. For NMOS transistors K p = 79.173u, W = 51u, L = 4u, and for PMOS K p = 19.485u, W = 102u, L = 2u. The V-I characteristics for NMOS and (17) PMOS transistors are given by Adby (1980), Tadeusiewicz (2001) and Tadeusiewicz and Kuczynski (2013) where v G represents the voltage at the gate node, v S the voltage at the source node, v D the voltage at the drain node, i S and i D represent the branch currents in source and drain, respectively. The value of k is calculated (Tadeusiewicz and Kuczynski 2013) as    Using the SP shown in Table 9, three operating points were found (Table 10). Figure 19a shows the resulting homotopic path for v 6 and Fig. 19b a zoom to the homotopic path, it is possible to notice the location of the operating points more clearly. (20)

Example 5
The circuit shown in Fig. 20 was introduced in Tadeusiewicz and Kuczynski (2013), it has three operating points. The transistors in this circuit were modelled using (20) and (21) and setting k = 0.5 mA/V 2 for every transistor, except for T12 and T4, for both transistors it was set to k = 1 mA/V 2 . The model was simulated in SPICE using the parameter values given in Example 4, except for T12 and T4, their parameters were set as follows: W = 51u, L = 2u, and W = 102u, L = 1u, respectively. The SP for the circuit is shown in Table 11. It was possible to find three solutions (see Table 12). The homotpic path is shown in Fig. 21. The obtained solutions are shown in Fig. 21d.

Numerical simulation and discussion
This section presents a performance comparison between the SHPT method proposed in this work and the MSA proposed in Vazquez-Leal et al. (2014). Also, the MSA method from Vazquez-Leal et al. (2014) was used without any modification. Besides, the inverted sphere technique was also applied to the MSA method and its performance evaluated. As will be seen, the SHPT path tracking method reduced, significantly, the computing time and the inverted sphere technique allowed, efficiently, to avoid the reversion phenomenon. Besides, an starting point criteria was employed that, in fact, eased the process to find multiple operating points using SHPT and MSA.   Tables 13, 14, and 15 show the comparison between the performance of three different simulations for each example studied in the previous section. SHPT stands for the proposed path tracking method, MSA2 is the method introduced in Vazquez-Leal et al. Example 1 and example 2 did not show any sign of the reversion phenomenon in any of the simulations performed. Nevertheless, the rest of the simulations using the MSA2 method showed reversion. No simulation using the MSA1 method showed reversion. It is important to make notice of the fact that the homotopic path is the same for the three types of simulation; this is because they are based on the same homotopy formulation, and the difference lies in the applied path tracking method. The results prove that the efficiency is improved using the inverted sphere technique proposed in this work, as it helped to avoid the revision phenomenon and allowed to perform the path tracking without issues. The cases where reversion phenomenon were noticed are marked with an asterisk (*), the computing time for these cases was not possible to be calculated because the method locked and no further calculations were possible.
As shown in Tables 13, 14, and 15, the time spent in straight line iterations is minimum compared to the time spent in iterations using hyperspheres. This characteristic allows the acceleration of the path tracking.
Computing time increases as the density of breakpoints grows. Nevertheless, for the worst case scenario, a curve with high density of breakpoints or a non-PWL curve, the proposed method would only spend the same computing time as the method proposed in Vazquez-Leal et al. (2014). The computing time spent by the MSA1 simulations and the SHPT method is noticeable different; the SHPT method performed up to twelve times faster. For the worst case scenario, the difference was 1.89 times faster than MSA1. The algorithm was implemented in Maple 15. Future work will focus on implementing the technique in Fortran, the goal is to improve the computation time and analyse larger circuits.

Conclusions
This work introduced a path algorithm for analysis of PWL circuits using the HCM, this algorithm exhibited improvements in the computing time compared to the algorithm proposed in Vazquez-Leal et al. (2014). Also, a starting point criteria was proposed in order to achieve better performance of the HCM, this criterion was proved useful but it does not assure an homotopic path that travels through all the root of the system. It just increase the probability of finding them. Furthermore, a technique for avoiding the reversion method was suggested and proved to be effective. By experimentation was possible to avoid reversion and allowed the continuation of the path tracking, nevertheless, for systems with a high density of variables (around one hundred) some instabilities were detected. Further work on this technique aims to improve the number of variables that is possible to work on. As it can be seen in Tables 13, 14, and 15, most of the simulation time is spent in NRM iterations, this leads to focus future work on the reduction of time spent on them. As a final comment, the SHPT path tracking method does not have a stop criterion, thus, the development of a stop criterion would be an important improvement.