AN EFFICIENT ALGORITHM FOR FINDING
MULTIPLE DC SOLUTIONS BASED ON SPICEORIENTED NEWTON HOMOTOPY METHODyyyyyDepartment of Electrical and Electronic Engineering,

A.Ushida , Y. Yamagami , I. Kinouchi , Y.Nishio , Y.Inouez

Tokushima University, Tokushima, 770-8506 JAPANzGraduate school of Integrated Science and Art, University of East Asia,
Simonoseki, Yamaguchi, 751-8503 JAPANABSTRACTFor circuit designing, it is very important to calculate thedc operating points. It is known that if the circuit contains
positive feedback loops such as ip- op and negative resistance circuits, it may happened to have many dc solutions.
It is very dicult to nd all of the solutions for these circuits. In this paper, we show a very simple Spice oriented
Newton homotopy method which can eciently nd out themultiple dc solutions, and perhaps all of the solutions.
1. INTRODUCTION

Many circuits such as ip- op, Schmitt-trigger, Hop eld
neural, negative resistance circuits [1,2], and so on have multiple dc solutions. There are also many important large scale
circuits which contain some of the above sub-circuits [3,11],
where a number of the dc solutions may be rapidly increasedby the combinations of the solutions from the sub-circuits.
In these circuits, it becomes very dicult to nd all of thesolutions. Therefore, many papers already have been published [3-14]. In the references [3,10,11], nonlinear resistive
elements are assumed to be piecewise-linear characteristic.
Especially, Yamamura's algorithm [11] can eciently nd
all of the solutions. The interval method [12-14] can nd allof the solutions of nonlinear algebraic equations in a given
accuracy. Although the method is elegant in the mathematical point of view, the computational eciency is not goodfor the large scale systems. Note that all of the above methods need to formulate the nonlinear algebraic equations for agiven circuit which is a dicult task for large scale circuits inthe practice. The others [4-9] are based on the continuation
method, which trace the solution curves [4,5], and nd thesolutions laying on the curves. The homotopy method [8]can trace all of the solutions if we can nd out all of the disjoint branches. However, the algorithm to nd branches is
not so easy when the variables of system are increased. Trajkovic et al.[9] have developed a Spice-like simulator having
a globally convergent property which is based on an arti cialparameter homotopy method. Inoue [6,7] also has developed
a simple Spice oriented simulator which eciently traces thedriving point characteristic curve, where the curve is traced
with the transient analysis tool of Spice, eciently. The
method can nd out all solutions only if they are laying on
the driving point characteristic curve. Thus, there is not
yet any practical Spice-like algorithm to nd all of the dc
solutions.

In this paper, we propose a Newton homotopy method
such that
F(x; ) = f (x) 0 (1 0 ) g(x0 ) = 0(1)0f (x) : Rn 7! Rn ; g(x0 ) : Rn 7! Rp
where n is the number of nodes, and p is the number of

positive feedback loops in the circuit 1 . In our homotopy
method, the initial point (1 0 )g(x0 ) with  = 0 is corresponding to the input current sources to the speci ed p
nodes, whose sources are continuously changed by a variable
, and the dc solutions are found at  = 1 on the solution
curve. By changing the initial guesses g(x0 ), we can trace
the another solution curves. In this way, it is possible tond the all of the solutions. We developed a user friendly
simulator using Spice to implement our Newton homotopy
method. From many benchmark examples, we found that
our simulator can nd all of the solutions on one branch if
we choose a suitable initial guess. It has been said that thecontinuation method is a time-consuming for large scale systems because we need to trace the solution curves in a large
dimensional space. However, we found from the theorem in
section 2 that we are enough only to consider the solution
curves in the p-dimensional space. p is usually very small
compared to n because it is equal to the number of positive
feedback loops. Therefore, the simulator can be usefully applied to the dc analysis of large scale ICs.
2. NEWTON HOMOTOPY METHOD

In this section, we consider the ideas of our Newton homotopy method. Let us assume the circuit equation of anonlinear resistive circuit as follows;f (v) = 0;(2)

Applying the Newton homotopy method to (2), we haveF(v; ) = f (v) + ( 0 1)f (v0 ) = 0(3)

where  is an additional variable and v0 is an initial guess.
Thus, (3) has n + 1 variables and n equations, so that the
1 We will prove in section 2 that the minimum number of variables to trace the solution curve is equal to the positive feedback
loops p[15,16].

V-447
0-7803-6685-9/01/$10.00©2001 IEEEv 2 Rn ; f (1) : Rn 7! Rnsolutions of (3) are given by the solution curves in (n + 1)space. For simplicity, put vn+1 = .
Then, we have from (3)F(v) = 0; F(1) : Rn+1 7! Rn(4)Let us trace the curve from the initial point v0 with vn+1 =0 (or  = 0). Then, the solutions of (2) are found at the
points vn+1 = 1 (or  = 1) on the solution curves. In our arclength method, a point v on the solution curve is describedby a function of the distance s from the initial guess v0 . It
satis esn+1
X2(ds) = (dvi )2(5)should be chosen in our solution curve tracing algorithm. We
call the indispensable variables to trace the solution curve
essential variables.
Theorem In our Newton homotopy method, the minimum
number of the variables to be chosen is equal to a numberof positive feedback loops contained in the circuit, whichmeans that one node voltage in each positive feedback loop
should be chosen as the essential variable.
proof The above curve tracing algorithm can trace the solution curve only if the rank of J is n +1. In this case, when
the solution curve traces the multiple solutions, it must passthrough the turning points p1 and/or p2 as shown in Fig.1[1,15,16,18]. Hence, it satis es the relationin the Euclidean space. Combining (5) and (4), we have thefollowing set of algebraic-di erential equations:
F1 (v1 ; v2 ; . . . ; vn ; vn+1 ) = 0
F2 (v1 ; v2 ; . . . ; vn ; vn+1 ) = 0(6)::::::::::::::::::::::::::Fn (v1 ; v2 ; . . . ; vn ; ; vn+1 ) = 0(11)d = 0; or dvn+1 = 0
dvi
dvi

i=1

n+1 
X
i=1

for some i at the turning points.
p1
ρ =1
or
vn+1=1dvi 2 = 1
ds

It can be solved by the backward-di erence method [17],where the derivative of kth order formular at s = sj+1 isgiven bydvi
= hk0 vij+1 + Qk;i (vij ; . . . ; vij0k+1 ) (7)ds s=sj +1
Substituting (7) into (6) at vj+1 , the algebraic-di erential

equations are transformed into the nonlinear algebraic equation as follows:
9
F1 (v1j +1 ; v2j+1 ; . . . ; vnj +1+1 ) = 0>>>F2 (v1j +1 ; v2j+1 ; . . . ; vnj +1=+1 ) = 0::::::::::::::::::::::::::::::::::::(8)>
j +1 j +1
j +1
>
F
(
v
;
v
;
.
.
.
;
v
)
=
0
n
>
1
2
n
+1
;
Pn+1
j +1
+ Qk;i (:::)g2 0 h2 = 0
i=1 f k0 vi
It is solved by the Newton method
vj +1 = vj 0 [J]01 v=vj H(vj )(9)where H(vj ) in (9) is a set of the functions given by (8) with
vj instead of vj +1 , and the Jacobian matrix is given by
2
3
@F1
@F1
@F1
...
@v1
@v2
@vn+1
@F2
@F2
@F2
6
7
...
@v1
@v2
@vn+1
6
7
6
7
j
.
.
J(v ) = 6 ..
7
.
.
6
7
@Fn
@Fn
@Fn
4
5
...
@v1
@v2
@vn+1
P1 (v1 ) P2 (v2 ) . . . Pn+1 (vn+1 ) v=vj
(10)
and
Pi (vi ) = 2 k0 [ k0 vi + Qk;i (vij ; . . . ; vij 0k+1 )]The solution curve can be continuously traced by the above
backward-di erence method [17] only if the matrix J(vj ) is
nonsingular on the solution curve. Now, we have the following important theorem about how many and which variables

s1

s2

s3

p

2

v

i

Fig.1 The schematic diagram of a solution curve,

s1 ; s2 ; s3 are solution points, and p1 ; p2 are turning points
Now, di erentiating (4) by the arc-length s, and applying

Cramer's formula, we haved = detjJ0(n+1);0(n+1) j dvi(12)ds
detjJ0(n+1);0i j ds
Thus, we nd from eqs.(11) and (12) that the cofactor
J0(n+1);0i must be nonsingular, but J0(n+1);0(n+1) is singular. Therefore, if the element Pi (vij ) at (n + 1; i) positionin J(vij ) is set to a nonzero value, the (n + 1) 2 (n + 1)
Jacobian matrix (10) will become nonsingular, and our algorithm can trace the curve at the turning points p1 andp2 . Therefore, we found that the node voltage vi must be
chosen as the essential variable.
In the practical Spice simulation, we need to set a probe at
the node i to get the voltage vi .
Q.E.D.
Corollary A node in positive feedback loop is found in sucha manner that when the node voltage at that point is xed,
the characteristic does not have a point satisfying the property given by (11), any more.
From the theorem, our homotopy method can be carriedout by solving the following set of algebraic-di erential equations:9f1 (v1 ; v2 ; . . . ; vn ) 0 (1 0 )I01 = 0 >
>

and

>
>
::::::::::::::::::::::::::::>
fp (v1 ; v2 ; . . . ; vn ) 0 (1 0 )I0p = 0 =fp+1 (v1 ; v2 ; . . . ; vn ) = 0>>>>::::::::::::::::::::::::::::>;
fn (v1 ; v2 ; . . . ; vn ) = 0

p 
X
i=1

dvi 2 +  d 2 0 1 = 0
ds
ds

(13)

(14)

where fv1 ; v2 ; . . . ; vp g are the essential variables, and fI01 ;. . . ; I0p g are the initial guess current sources.
V-448

3. REALIZATION OF NEWTON
HOMOTOPY METHOD

We have developed our Newton homotopy simulator using
Spice as shown in Fig.2.
ρ

ρ -1

Σ

4.1 Transistor-diode circuit

VCCS

-1
STC

I 0 i , i=1,2,...pInitial guessvi

i=1,2,...p( ρ -1 ) I 0 i

Nonlinearresistive circuit

i=1,2,...pFig.2 Circuit diagram of our Newton homotopy method
We rst set the initial guesses fI01 ; . . . ; I0p g, and get
the input current sources ( 0 1)I0i ; i = 1; . . . ; p with
VCCSs(voltage-controlled current source) and multipliers.
Secondly, the voltage  satisfying the relation (14) can beproduced by STC(solution curve tracing circuit) [6] as shownin Fig.3.
.
ρSchmitt-trigger circuit, whose sub-circuits have 3 dc solutions, two of them are stable and the other is unstable [18].
Negative resistance circuits [2] have also multiple solutions.
If these sub-circuits are coupled in each other, then, they
may have many dc solutions. In this section, we solve some
important benchmark problems.
Fig.4 is a well-known benchmark circuit [3-11] having 11
dc solutions. The circuit consists of 3 Schmitt-trigger circuits. The circuit is constructed in such manner that theoutput voltages from 2 sub-circuits in the right hand side
are inputted into the left hand side circuit through diodes.
In this example, we used 3 probes at 3 positive feedback
loops as shown in Fig.4. With the initial guess I01 = 1[mA]
I02 = 1[mA] I03 = 5[mA], we can trace all of the solutionsat  = 1 on the solution curve as shown in Fig.5. We also
found that our Newton homotopy method can trace all of
the solutions for any initial guess in relatively large region.
The solutions are shown in Table 1.
ρ -1( ρ -1 ) I 03
VCCS

I 01

I s = kI v

VCCS

I 02

( ρ -1 ) I 01

2.7k

1k
15k

15k

3k

2.7k

.

I ρ=k I ρ 2

ρv= Σ vi2

R

1k

3k

11k

15k

5k

12V

V1

STC

Int

11k

V2

Σ

a

I 03

( ρ -1 ) I 02

-1

ρ

p .

VCCS

15k

15k

i=1

1k

50

V3

1[A]

1k
50

k I =1k I =1VCCS

N VCCS

Fig.4 Transistor-diode circuit

Fig.3 Solution curve tracing circuit (STC)
The input voltages to STC are the essential node voltages
v1 ; . . . ; vp in the positive feedback loops. At rst, each voltage is transfered to the corresponding current source by
VCCS, so that the STC does not give any disturbance to thenonlinear resistive circuit. The transfered current sources
are di erentiated by the time t instead of the arc-length s by
the use of di erentiators or inductor circuits. Thus, we havev_ 1 ; . . . ; v_ p , which are again squared to produce v_ 12 ; . . . ; v_ p2
with multipliers, and they are summed by an adder again.
Thus, we have
p
X
v=
v_ i2
i=1

The voltage v is further transformed to the corresponding
current source Is as shown in Fig.3. Note that if we assume
the node voltage at node "a" is _ , I = _2 is easily produced
by a multiplier and a VCCS. Thus, we can realize the circuit
satisfying the relation (14). Furthermore, the node voltage
_ at "a" is integrated, and we can get . Note that R in
Fig.3 is a very large resistance only to avoid the L-J cutset.
4. ILLUSTRATIVE EXAMPLES

Large scale circuits usually consist of many sub-circuits
having positive feedback loops, such as ip- op circuit and

ρ=1

V3

4
2
1

S6

S2 S3 S4 S 5

S1

S9

S8

0

S11
S10

S7
V2

ρ

V1

-2
-4
0

20

40

60

80

s

100

120

140

160

180

Fig.5 A set of solution curves with our Newton homotopy
method, where we get the solutions at  = 1.
4.2 Flip- op and Schmitt-trigger circuits

Let us apply our Newton homotopy method to a circuit
shown in Fig.6, where ip- op and Schmitt-trigger subcircuits are coupled with a resistor R0 . We found that, for alarge R0 , the two sub-circuits are weakly coupled each other,and it has 9 dc solutions. If we choose the smaller R0 , thenumber of solutions will be decreased because of their strong
coupling. We need 2 probes at each positive feedback loop
as shown in Fig.6. We can trace 9 dc solutions with theV-449

initial guess I01 = 3[mA] I02 = 5[mA] and R0 = 2[k ]. In
the strongly coupled circuit with R0 = 10[ ], however, thecircuit only has 3 dc solutions.
ρ -1
ρ

STC

VCCS

-1
Σ

I 01

4k

VCCS

( ρ -1 ) I 01

4k

( ρ -1 ) I 02

I 02

1k

3k

11k

2k
V1

8k

8k

15k

R0

V2

12V
T3

T

T1

2

0.1k

T4

0.1k

1k

50

Fig.6 Flip-Schmitt circuit
6
V2

V [V]

V1

4
ρ
2
S1

S3
S5

0

-2

S6 S7

S4

ρ=1

S2

0

20

40

s

60

S9
S8

80

100

Fig.7 A set of solution curves with our Newton homotopy
method, where we get the solutions at  = 1.
Table Dc solutions of transistor-diode circuit
Solution
V1
V2
V3
S1
S2
S3
S4
S5
S6
S7
S8
S9
S10
S11

0.719
0.691
0.036
0.028
0.027
0.023
0.024
0.028
0.027
0.023
0.021

we need to try some curve tracing algorithms from the different initial guesses. Thus, it is possible to nd all of thesolutions. However, how to nd the all of the solutions is
the future problem.

10.222
10.221
10.216
1.134
0.271
0.271
1.135
10.21810.2231.134
0.271

10.222
10.22010.222
10.22310.223
1.135
1.135
1.135
0.271
0.271
0.271

5. CONCLUSIONS AND REMARKSWe developed a user friendly simulator to nd the multiple dc solutions based on our Newton homotopy method.
The simulator is consisted of a simple solution curve tracing
circuits(STC), an adder(6) and voltage controlled current
sources(VCCS). The solution curves are traced by the use oftransient analysis tool of Spice, and the solutions are foundon the curves at  = 1. In order to trace the multiple solutions, we need to set one probe at each positive feedback
loop to get the essential node voltage. If we choose a suitable
initial guess, it may trace all of the dc solutions. Otherwise,

References[1] L.Jiang, Y.Nishio and A.Ushida, \Stability of characteristic curves of nonlinear resistive circuits," IEEE Trans. on
Circuits and Systems-I, vol.45, pp.634-643, 1998.
[2] L.O.Chua, J.Yu and Y.Yu, \ Bipolar-JFET-MOSFET negative resistance device," IEEE Trans. on Circuits and Systems, vol.CAS-32, pp.46-61, 1985.
[3] M.Tadeussiewicz,\DC analysis of circuits with idealized
diodes con sidering reverse bias breakdown phenomenon,"IEEE Trans. on Circuits and Systems-I, vol.44, pp.312326,1997.
[4] A.Ushida and L.O.Chua, \Tracing solution curves of nonlinear equations with sharp turning points," Int. Jour. of
Circuit Theory and Applications,pp.1-21, 1984.
[5] C.W.Yun and K.S.Chao, \Simple solution curves of nonlinear resistive networks," Int. Jour. of Circuit Theory andApplications, vol.11, pp.47-55, 1983.
[6] Y. Inoue \DC analysis of nonlinear circuits using solutiontracing circuits," Trans. of IEICE (A), vo.J74-A, pp.16471655, 1991.
[7] Y.Inoue \A practical algorithm for DC operating-point
analysis of large scale circuits," Trans. of IEICE (A), vo.J77A, pp.388-398, 1994.
[8] J.Lee and H.-D. Chiang, \ Constructive homotopy methods
for nding all or multiple dc operating points of nonlinear
circuits and systems," IEEE Int. Symp. Circuits and Systems, pp.IV-525-528, 2000.
[9] R.C.Melville, L.Trajkovic, S.-C.Fang
and L.T.Watson, \Arti cial parameter convergent homotopy methods for the dc operating problems," IEEE Trans.
on Computer-Aided Design of Integrated Circuits and Signal Processing, vol.6, pp.861-877, 1993.

[10] T.Ohtsuki, T.Fujisawa and S.Kumagai, \Existence theorems
and a solution algorithm for piecewise-linear resistive networks," SIAM J. Math. Anal., vol.8, no.1, pp.69-99, 1977.
[11] K.Yamamura and T.Ohshima, \Finding all solutions ofpiecewise-linear resistive circuits using linear programing,"
IEEE Trans. on Circuits and Systems-I, vol.45, pp.434-445,
1998.
[12] G.Alefeld and J.Herzberger, Introduction to Interval computations, New York,NY* Academic, 1983.
[13] L.Kolev and V.Mladenov, \ An interval method for nding all operating points of nonlinear resistive circuits,"
Int.J.Circuit Theory and Applications,vol.18, pp.257-267,
1990.
[14] L.Kolev,\An interval method for global nonlinear analysis,"IEEE Trans. on Circuits and Systems-I, vol.47, pp.675-683,
2000.
[15] R.O.Nielsen and A.N.Willson, \A fundamental results concerning the topology of transistor circuits with multiple equilibria," Proc. IEEE vol.68, no.2, pp.196-208, 1980.
[16] T.Nishi and L.O.Chua, \Topological criteria for nonlinear
resistive circuits containing controlled sources to have aunique solution," IEEE Trans. on Circuits and Systems,
vol.CAS-31, pp.722-741, 1984.
[17] L.O.Chua and P-M Lin, Computer Aided Analysis of Electronic Circuits: Algorithms and Computational Techniques,Prentice-Hall, 1975.
[18] M.M.Green and A.N.Willson,\How to identify unstable operating points,"IEEE Trans. on Circuits and Systemsi-I,
vol.CAS-39, pp.820-832, 1992.

V-450