SPICE を用いた複数個の動作点解析- 全解探索への挑戦 SPICE Oriented Analysis of Multiple Operationg Points- Challenge to find out all the solutions 高田 幸司 †，木内 一慶 ‡，山上 喜廣 †，西尾 芳文 †，牛田 明夫 †Kouji TAKATA†, Ikkei KINOUCHI‡, Yoshihiro YAMAGAMI†, Yoshifumi NISHIO† and Akio USHIDA†† 徳島大学 工学部‡ 富士通(Tokushima University)(Fujitsu)

Email ushida@ee.tokushima-u.ac.jp
ロップ回路やシュミット・トリガー回路や負性抵抗回路

1がある．そこで，このような複数個の解を持つ回路にあらまし

対する全解探索アルゴリズムが幾つか提案されている[5-7]．その一つに区分線形化法がある．すなわち，非線直流解析は回路設計における重要なステップである

形特性を区分線形近似を行い領域ごとに解析すればよが，近年，IC の大規模化と複雑化に伴い，特に，大規

いが，あらかじめ解が存在しない領域を削除することに模な回路の複数個の動作点を求めることが必要になってよって効率的に全解を探索する方法が文献 [7] に提案さきた．そこで，ここでは，ニュートン・ホモトピー法にれている．また，区間解析法 [5-6] は連続系にも適用で基づいた SPICE の過渡解析ツールを用いた解析手法を

き理論的にはエレガントであるが，方程式の規模が大き提案する．この場合，回路に含まれている正帰還ルー

くなるにしたがって演算効率が悪くなる．

プの各１点からの節点電圧を変数として取り出すことによって複数個の動作点を探索できる手法を提案する．

一方，複数個の動作点を追跡するための連続法にホ

この手法は回路方程式やヤコビー行列を求めることな

モトピー法がある [8-12]．ここで，ホモトピー・パスをく動作点解析ができるために汎用性の点で優れている．

追跡する各種のアルゴリズムが [10-12] に報告されていまた，全解を探索するために解曲線が孤立した枝を持つ

る．特に，文献 [11] では非線形代数方程式の複数根をような場合にも適用できるアルゴリズムを提案する．

追跡するための汎用シミュレータ (HOMPACK) が報告されている．また，文献 [12] では電子回路のようなス

2

ティフな系にも適用可能な後退差分公式に基づいた解

はじめに曲線追跡アルゴリズムが提案されている．HOMPACK

直流解析は電子回路を設計するために最も重要なスを直流解析に応用したものに文献 [13,14] がある．そこ

テップの一つであり, このことは非線形代数方程式を解

では，各ステップ毎に SPICE からヤコビー行列を参照

くことに帰着される．従来から，SPICE などではこのするようなシステムになっており，このヤコビー行列を目的のためにニュートン法が用いられているが，一般呼び出す回数が多くなるに従って，このシミュレータのに，この手法は複数個の動作点を持つような回路の解析

演算効率は悪くなる．駆動点特性曲線を追跡するためには適用できない．一方，正帰還ループを持つような回に文献 [12] を応用した解曲線追跡回路が井上によって路は複数個の動作点を持つ可能性のあることが報告さ提案されている [15-16]．このアルゴリズムは SPICE のれている [1-4]．その代表的な回路としてフリップ・フ

過渡解析ツールを直接用いるため汎用シミュレータと1

して非常に適している．さらに，これをニュートン・ホ

代数方程式をニュートン法によって解くことによって決

モトピー法を適用できるようにしたものを文献 [17,19]

定される．
j+1
)=0
f2 (v1j+1 , v2j+1 , . . . , vn+1ての解を探索できることもあるが，異った枝や閉曲線が存在する場合には探索されない．このような場合には，....................................
j+1
fn (v1j+1 , v2j+1 , . . . , vn+1
)=0
n+1
j+1
} + Qki (・) − h2 = 0
i=1 {αk0 vi

孤立したホモトピー・パス上の１点であるリミット・ポ
イント [18] を見付けることができれば，孤立した解曲線が追跡されることになり全解探索の手がかりとなる．
ここでは，さらに，孤立した駆動点特性曲線の追跡方法

(5)










このヤコビー行列はは次のように与えられるjJ(v ) =

についても議論する．

3











j+1
f1 (v1j+1 , v2j+1 , . . . , vn+1)=0で報告した．この解曲線追跡法は初期値に依存して全

弧長法による解曲線追跡法
駆動点特性曲線やホモトピー法などによって得られた非線形代数方程式は回路変数の他に補助変数を一個多く含んでおり，次のような非線形代数方程式で記述さ











∂f1
∂v1
∂f2
∂v1

∂f1
∂v2
∂f2
∂v2

...

∂fn
∂v1

∂fn
∂v2

..
.

..
.

P1 (v1 ) P2 (vi2 )れる．

f1 (v1 , v2 , . . . , vn , vn+1 ) = 0
f2 (v1 , v2 , . . . , vn , vn+1 ) = 0..........................

...

...

∂f1
∂vk
∂f2
∂vk

...

∂f1
∂vn+1
∂f2
∂vn+1

...

∂fn
∂vk

...

∂fn
∂vn+1

..
.












. . . Pk (vk ) . . . Pn+1 (vn+1 )ただし，P1 (v1 ), . . . , Pn+1 (vn+1 ) は後退差分公式より決(1)

定され，次式の形に与えられる．

Pi (vi ) = 2αk0 [αk0 vi + Qki (vij . . . vij−k+1 )]fn (v1 , v2 , . . . , vn , vn+1 ) = 0
ニュートン・ホモトピー法の場合にはこの方程式は次のすなわち，Pi (vi ) は i 番目の節点電圧を微分することによって得られるが，SPICE 上ではプローブをこの点ように書かれる．

F(v) = f (v)+(ρ−1)f (v0 ) = 0, f (・) : Rn → Rn (2)に設定し，電圧を解曲線追跡回路に取り込むことによっこの方程式の解は (n + 1) 次元空間における解曲線としめには J(v ) の階数が n + 1 であればよい．このとき，て与えられる．そこで，解曲線上の任意の一点を初期値

プローブの数に関して次のような定理が成立する．

v0 に選び，その点からの曲線の弧長を s とするならば，[定理] ニュートン・ホモトピー法を用いて解曲線を追跡

ユークリッド距離としてするために必要とする変数の最少数は回路に含まれて(ds)2 =n+1(dvi )2(3)i=1が成立する．式 (1) と (3) 連立させることによって次のような微分代数方程式が得られる．

..........................

jいる正帰還の数に等しい. すなわち, 各正帰還ループから１つの節点を選び, その節点の電位を変数として選べ
ばよい [4].

[証明] 解曲線が複数個の解を追跡するものと仮定しよう. このとき, 解曲線はρ = 1 の点を何回か横切る.

f1 (v1 , v2 , . . . , vn , vn+1 ) = 0
f2 (v1 , v2 , . . . , vn , vn+1 ) = 0て実現される [15]．この解曲線が連続的に追跡されるた

p1

(4)

ρ =1
or
vn+1=1fn (v1 , v2 , . . . , vn , , vn+1 ) = 0
n+1

i=1

dvi
ds

s2

s3

p

2

2

=1

v

式 (4) の微分代数方程式の求解には後退差分公式が有効∗

s1である ．この場合，s = sj+1

における解は次の非線形∗ SPICE では本求解に対し，予測子ー修正子と可変ステップサイ

ズ方式の数値積分法が採用されており効率よく数値積分を行なっている．

i

図１ 複数個の解を追跡する解曲線の概略図

この概略図が図 1 に示されている. すなわち, s1 , s2 , s3は解であって, p1 , p2 は変数 vi に関する変曲点であり次式を満足する.

v=vj

dvi
dρ
dρ
= 0,
= 0
= 0, ⇒
dvi
ds
ds

(6)

4

ニュートン・ホモトピー法による
動作点解析

このとき, i 番目のノード電圧を変数に選ぶならば解曲線は連続的に追跡されることを示そう. (1) 式を弧長 sで微分し,∂fi dv1∂fi dv2∂fi dvn+1++ ... += 0 (7)∂v1 ds∂v2 ds∂vn+1 dsi = 1, 2, . . . , n

2

上で実現した回路が図２に示されている [15,17]．すな
わち，非線形抵抗回路中の p 個の正帰還ループの節点電圧を選び，これらが弧長の関係式 (4) の最後の式を満足するように解曲線追跡回路 (STC) を構成する．

(7) 式の係数行列を用いて, 次の行列式を考える.
∂fi
 ∂fi ∂fi . . .
 ∂v1 ∂v2
∂vn+1 
 ∂f

∂f1
∂f1
1


 ∂v1 ∂v2 . . . ∂vn+1  ... ... ...
. . . Ai =  ∂f∂f∂fi
 ∂v1i ∂v2i . . . ∂vn+1 ... ... ...
. . . 
 ∂fn ∂fn
 ∂v
. . . ∂v∂fn ∂v
1

前節で述べた弧長法に基づいた解曲線追跡法を SPICE

ρ

vi

n+1

電圧制御
電流源

-1

解曲線
追跡
回路

(8)

ρ−1

Σ

V0,i=I 0 ,i
初期値設定

i=1,2,...,p

この行列式は常に零であり, これを第 1 行目で展開する(ρ−1)Ι 0,i

ことにより
∂fi
∂fi
∂fi
D1 (s) +D2 (s) + . . . +Dn+1 (s) = 0
∂v1
∂v2
∂vn+1

i = 1, 2, . . . , ni=1,2,...p

(9)

ここで,

図２ SPICE 上での回路構成

ここで，I0i は初期条件である．これを適用して図３ (a)Di (s) =


 ∂f1
 ∂v1

 ∂f2
i+1  ∂v1
(−1)

 ...

 ∂fn
 ∂v1

非線形回路

∂f1
∂vi+1
∂f2
∂vi+1

...

...

∂f1
∂vi−1
∂f2
∂vi−1

...

...

...

...

...

∂fn
∂vi−1

∂fn
∂vi+1

...

...

...



∂f1
∂vn+1 ∂f2
∂vn+1 
. . . 

∂fn

∂vn+1

の電圧制御回路 [20] を解析してみよう．
VSCC
20V

Q1

start

STC

Q3

Q

5X

(11)

v

Q
out

start

Q

start

3.54k

5

5k

1k

1k
Startup circuitry5k

図３ (a) 電圧制御回路(12)このようにして, (6) 式と (10) 式から Di (s) は非零であることが分かる. 従って, ヤコビー行列 J(vj ) の (n+1, i)の要素に非零要素の pi (vi ) を設定すると, J(vj ) は正則となり，この変曲点も連続的に追跡できることが分かる. このようにして, 各正帰還回路中に 1 つの節点を選
び, その節点にプローブを設定することにより解曲線が安定に連続的に追跡される. [証明終]

1X

-1

100k

Q2

Q4

(10)

Σ
ρ

あるから, (7) 式と (9) 式から次の関係が得られる.

k が次のように求まる.
2(s)k = ± D12 (s) + D22 (s) + . . . + Dn+1

I0

5k

5k

(1) 式における n × (n + 1) のヤコビ行列の階数は n でdvi= kDi (s), i = 1, 2, . . . , n + 1
ds(10) 式を次式に代入することにより


2 
2
2
dv2
dvn+1
dv1
+
+ ... +
=1
ds
ds
ds

ρ−1

(ρ−1)I 0

20

[V]
15

vout
10
ρ
5
ρ=1

S1

S2 S3

0
0

10

20

30

40

S

50

(b) ニュートン・ホモトピー法による解曲線

この回路は温度変化に強いことが知られている．まず，f1 (v1 , v2 ) + (ρ − 1)I01 = 0

Q4 のベースにプローブを設定して解曲線追跡回路を構

f2 (v1 , v2 ) + (ρ − 1)I02 = 0

成する．I0 = 5[mA] を初期値に選んで追跡を行ったところ，ρ = 1 で３個の全解を追跡している．解 S2 から(14)

ここで, STC は解曲線追跡回路である.

S3 の範囲ではρの値はほぼ１となっており，この領域で応用例 5.1は近似的に関数値が零に近いことを示している．従っ

次式は化学プラントの反応方程式から得られたものであて，ニュートン法などでは解析は困難である．

り, 7 つの解が存在していることが知られている．文献

5

る範囲に選べば全解が追跡されることが示されている．

[21] ではこれに定点ホモトピー法を適用して初期値をあ解曲線が孤立した閉曲線を持つ場
合の解析

f1 (v1 , v2 ) = 2v13 + 2v1 v2 − 22v1 + v22 + 13 = 0

ここでは，回路方程式が式 (1) のように非線形代数方程式系で記述された系を SPICE を用いて解くことを考えよう．SPICE には高機能なアナログビヘービア・モデルが組み込まれており，これらを用いて与えられた代

f2 (v1 , v2 ) = v12 + 2v1 v2 + 2v23 − 14v2 + 9 = 0

上 式 に お い てΓ1 :I01 =40,I02 =50 な る 初 期 値 を 選 び,

SPICE を用いてニュートン・ホモトピー法を適用した
結果が図６に示されている.

数方程式を記述することは容易である．n 変数に対する3

式 (1) の連立非線形代数方程式を電圧制御電流源を用い2.52て実現したのが図４である．ここで，

1.5

Fi (v1 , v2 , · · · , vn ) = fi (v1 , v2 , · · · , vn ) + vi

(13)

1

ρ

1 [Ω ]

-

+

+

v2

vn

1 [Ω ]

-

5

10

s3 s4

s5 s6

s7

-0.5
1 [Ω ]

0

15
20
s[sec]

2530

図６ ニュートン·ホモトピー法による解曲線

F ( v1 , v 2 , .... , v )2nF1 ( v1 , v 2 , .... , v )
n

s2

0

-1
v1

Γ1

s1

0.5

i = 1, 2, . . . , n+(15)Fn ( v1 , v 2 , .... , v )nΓ1 を (v1 ,v2 ) 平面で描いたのが図７であり, 1 本の解曲線で全ての解 (s1 , s2 , · · · , s7 ) を追跡できている. しか

図４ 非線形代数方程式の実現しながら，Γ2 :I01 =80,I02 =-90, Γ3 :I01 =-4,I02 =4 なる初

この非線形代数方程式を２変数の場合にニュートン・ホ

期値を選ぶと, それぞれの解曲線は一部の解のみ追跡しモトピー法を適用して解こうとすると図５のような回ていることが分かる．

路構成となる．
v1

86(ρ -1 ) I01
I01

1 [ Ω]vv2
1 [ Ω]

Γ1
Γ3

2

20

VCCS
F1(v1 , v 2)

Γ2

4

VCCS

I02

(ρ -1 ) I02

ρ -1

-2
-4

Σ

STC
ρ

-6
-6

Γ1
-4

-2

0

2

4

6

v1
-1

F2 (v 1, v 2)
図５ ニュートン・ホモトピー法の実現回路

図７ 異なる初期値による解曲線

そこで，Γ2 の場合ついて考える．ホモトピー・パス
として同じパラメータに対して Γ2 の他に独立した解曲線が存在することが予想される．そこで，その解曲線のi in

リミット·ポイントを求め，その点を初期値にとること+ v4 −

により, 独立した解曲線を追跡することを考える．
すなわち，ニュートン・ホモトピー法の関係式

F1 (v1 , v2 , ρ) = f1 (v1 , v2 ) + (ρ − 1)I01 = 0

i4


(16)

F2 (v1 , v2 , ρ) = f2 (v1 , v2 ) + (ρ − 1)I02 = 0

vin

からρを消去すると次式が得られる．

ここで，v2 から見たリミット・ポイントにおいては式(17) より∂F dv1∂F+=0∂v1 dv2∂v2(18)を満足する必要がある．このようにして

vin

STC

+ v4 −
R3= 1

+
v2
−
+
v1
−

(a)

F (v1 , v2 ) = I20 f1 (v1 , v2 ) − I10 f2 (v1 , v2 ) = 0 (17)F3 (v1 , v2 , I01 , I02 ) =i+ 2
v3 R3= 1
−
i1

+
v2
−
+
v1
−

(b)

図９ (a) 直並列抵抗回路，(b) 駆動点特性追跡回路

ここで，電圧制御型非線形抵抗の特性を次のようにおく
i1 (v1 ) = 7.778v13 − 25.085v12 + 20.54v1



3
2
i2 i(v2 ) = 4.313v2 − 17.466v2 + 18.11v2 
(20)

i4 (v4 ) = 1.0406v45 − 10.731v44 + 39.335v43 −60.373v 42 + 34.138v4
ここで，駆動点特性曲線群のリミット・ポイントを求め

∂F
=0
∂v2

(19)

るために次の回路方程式

f2 (v1 , v2 , v3 , v4 , vin ) = v3 − v2 − v1 = 0

を同時に満足する解をニュートン·ホモトピー法で探索

f3 (v1 , v2 , v3 , v4 , vin ) = i4 (v4 ) − v3 − i2 (v2 ) = 0 f4 (v1 , v2 , v3 , v4 , vin ) = i2 (v2 ) − i1 (v1 ) = 0

することにより図８におけるリミット・ポイントが４個求まった．これらの初期値からニュートン・ホモトピー法を実行することにより全解が図８のように探索された．これは一種の２次元ホモトピー法であることに注意
しよう．








f1 (v1 , v2 , v3 , v4 , vin ) = vin − v4 − v3 = 0

なる補助方程式がえられる．従って，式 (16) と式 (19)と，補助方程式として，リミット・ポイントの条件式

f5 (v1 , v2 , v3 , v4 , vin ) =

dvin
=0
dv4

(22)

を追加したものにニュートン・ホモトピー法を適用する．
8

このときρ = 1 のところで各特性曲線に対応するリミッ

6
4

ト・ポイント１８個が次のように得られた．

Γ2
S4

S2

v2 0

解曲線 (Γ1 )[s1 , s2 , s5 , s6 , s13 , s15 , s17 ]

Γ2

S3

2

解曲線 (Γ2 )[s8 ]，解曲線 (Γ3 )[s3 , s4 , s9 , s11 ]S5-2

解曲線 (Γ4 )[s10 ], 解曲線 (Γ5 )[s16 ], 解曲線 (Γ6 ) [s7 ]
S7

S6

-4

S1

-6
-6

-4

-2

v1

0

2

解曲線 (Γ7 )[s18 ], 解曲線 (Γ8 )[s12 , s14 ]

4

6

20
15

図８ ニュートン·ホモトピー法による

10

リミット·ポイントの求解

ρ

5
s1s3
0 s2

例題 5.2-5駆動点特性曲線が複数個の枝や閉曲線から構成されて-10

いるような場合に，孤立した全てのリミット・ポイント-15を見付けることができるならば全ての閉曲線が探索さ-20れることになる．このような例として図９ (a) の回路を考えることにしよう [12]．

0

s17

s18

ρ=1

50

100

150

200
s[sec]

250

300

350

図１０ ニュートン・ホモトピー法による解曲線(21)[7] K.Yamamura and T.Ohshima, “Finding all solutionsof piecewise-linear resistive circuits using linear programming,” IEEE Trans. on Circuits and Systems-I,
vol.45, pp.434-445, 1998.

10
8
Γ2

Γ1

6

Γ3

Γ7

Γ4

Γ8

Γ6

Iin

[8] W.I.Zangwill and C.B.Garcia, Pathways to Solutions,Fixed Points, and Equilibria, Prentice-Hall, 1981.

4
Γ5

2
0
0

1

2

3

4

5

6

7

8

9

10

Vin

図１１ 駆動点特性曲線

この結果は全てのリミット・ポイントを求めてはいないが，少なくとも閉曲線上の各１つのリミット・ポイントが求まっており，これらの点から出発して全ての孤立した閉曲線を求めることができた．

6

結論

ここで提案した解曲線追跡法に基づいたホモトピー
法は SPICE の過渡解析ツールが利用されるために複数個の解を容易に追跡することが可能である．特に，我々のシミュレータは多くのトランジスタや MOSFET からなる回路に対して回路方程式を誘導したり，そのヤコビー行列を求める必要がないために汎用性の点で優れている．しかしながら，ホモトピー法で全解を探索する
ことは保証されていない．そこで，本稿では孤立した解曲線を探索するためにリミット・ポイントを見付ける方
法について議論した．この方法も完全ではないが，全解探索ために有効な一手法となるであろう．

参考文献[1] R.O.Nielsen and A.N.Willson, ”A fundamental results
concerning the topology of transistor circuits withmultiple equilibra,” Proc. IEEE vol.68, no.2,pp.196208, 1980.
[2] L.O.Chua, J.Yu and Y.Yu, ”Negative resistance devices,” Int. Jour. of Circuit Theory and Applications,vol.11, pp.161-186, 1983.
[3] 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.
[4] T.Nishi and L.O.Chua, ”Topological criteria for nonlinear resistive circuits containing controlled sources
to have a unique solution,” IEEE Trans. on Circuits
and Systems, vol.CAS-31,pp.722-741, 1984
[5] 奥村浩士, 佐伯秀一, 木嶋昭, ”区間解析による非線形回路方程式の求解アルゴリズムに関する一考察,” 電子通信
学会論文誌, vol.J69-A, pp.489-496, 1986.
[6] L.Kolev,“An interval method for global nonlinear
analysis,” IEEE Trans. on Circuits and Systems-I,
vol.47, pp.675-683, 2000.

[9] K.S.Chao, D.K.Liu and C.T.Pan, “A systematic
search method for obtaining multiple solutions of simultaneous nonlinear equations,” IEEE Trans. on
Circuits and Systems, vol.CAS-22, no.9, pp.748-753,
1975.
[10] L.T.Watson, “Globally convergent homotopy algorithms for nonlinear systems of equations,”Nonlinear
Dynamics, vol.1, pp.143-191, 1990.
[11] L.T.Watson, S.C.Billups and A.P.Morgan, “ALGORITHM 652 HOMPACK: A suite of codes for globally convergent homotopy algorithms,” ACM Trans.
on Math. Software, vol.13, no.3, pp.281-310, 1987.
[12] A.Ushida and L.O.Chua, “Tracing solution curves ofnonlinear equations with sharp turning points,” Int.
Jour. of Circuit Theory and Applications,vol.12, pp.121, 1984.
[13] R.C.Melville,
Lj.Trajkovic,
S.C.Fang and L.T.Watson, “Artiﬁcial parameter homotopy methods for the dc operating point problem,”IEEE Trans. on Computer-Aided Design of Integrated
Circuits and Signal Processing, vol.12, no.6, pp.861877, 1993.
[14] S.C.Fang, R.C.Melville and Lj.Trajkovic, “Artiﬁcial
parameter homotopy methods for the dc operating
point problem,” Patent No.5,181,179, Jan. 19, 1993.
[15] 井上靖秋, ”解曲線追跡回路を用いた非線形回路の直流解
析,” 電子情報通信学会論文誌 (A), vol.J74-A, pp.16471655, 1991.
[16] 井上靖秋, ”大規模回路の直流動作点解析法,” 電子情報通信学会論文誌 (A), vol.J77-A, pp.388-398, 1994.
[17] A.Ushida, Y.Yamagami, Y.Nishio, I.Kinouchi and
Y.Inoue, “An eﬃcient algorithm for ﬁnding multiple
DC solutions based on SPICE-oriented Newton homotopy method,”IEEE Trans. CAD, to appear.
[18] M.Kubicek and M.Marek, Computational Methods
in Bifurcation Theory and Dissipative Structures,
Springer-Vrlag, New York, 1983.
[19] 木内，山上，西尾，牛田，“ニュートン・ホモトピー法を用いた直流回路の SPICE シミュレーションについて，”
信学技法 NLP200-7-20, pp.35-42, 2000.
[20] L.B.Goldgeisser and M.M.Green, “A novel algorithm
that ﬁnds multiple operating points of nonlinear circuits automatically,” Proc. IEEE Int. Symp. Circuits
and Systems, Monterey, III, pp.558-561, 1998.
[21] M.Kuno and J.D.Seader, “Computing all real solutions to systems of nonlinear equations with a global
ﬁxed-point homotopy,” Ind. Eng. Chem. Res., vol.27,
pp.1320-1329, 1988.