社団法人 電子情報通信学会
THE INSTITUTE OF ELECTRONICS,INFORMATION AND COMMUNICATION ENGINEERS

信学技報
TECHNICAL REPORT OF IEICE.

SPICE を用いた平均化法による解析アルゴリズムの提案
片岡大茂†

上手 洋子†

山上 喜廣†

西尾芳文†† 徳島大学〒 770-8506 徳島県徳島市南常三島町 2-1
E-mail: †{hiroshige,uwate,yamagami,nishio}@ee.tokushima-u.ac.jp
あらまし

非線形システムの解析において、さまざまな近似解法が用いられる。特に弱非線型システムにおいては平均化法が用いられることが一般的である。我々は平均化法を満たす回路を SPICE 上で実現し、さらにニュートンホモトピー方を用いて解析することにより、平均化法で得られる複数の平衡点を 1 度のシミュレーションで探索するアルゴリズムを提案する。例題として 2 つの van der Pol 発振器をインダクタで結合した結合発振器の解析を行う。
キーワードSPICE, 平均化法, ニュートンホモトピー法

SPICE-Oriented Algorithm for Coupled Oscillators
by Using Averaging Method
Hiroshige KATAOKA† , Yoko UWATE† , Yoshihiro YAMAGAMI† , and Yoshifumi NISHIO†† Dept. of E.E. Eng., Tokushima University, 2-1 Minami-Josanjima, Tokushima 770-8506, JAPANE-mail: †{hiroshige,uwate,yamagami,nishio}@ee.tokushima-u.ac.jp
Abstract For designing oscillatory circuit, it is important to simulate and search the characteristics. If we analyze
a weakly nonlinear circuit, we often apply averaging method which is one of the approximate solving method forthe weakly nonlinear systems. In this study, we combine the averaging method to Newton homotopy method andanalyze a circuit by using SPICE. By solving Newton homotopy method, we obtain multiple equilibrium points in asingle SPICE simulation. As an example, we analyze the property of two coupled van der Pol oscillators. By using
our proposed SPICE-oriented algorithm, we obtain three equilibrium points. The result show us, our proposed
method is convenient for search the equilibrium points in averaging method.
Key words SPICE, averaging method, newton homotopy method

1. ま え が き

複数平衡点を一度の解析で探索するアルゴリズムを提案する。
平均化法は、発振器の接点方程式を正規化した連立方程式に対結合発振器の同期現象の解析は、蛍の群集発光や生物の鼓動、して定常発振状態を仮定した変数を与え、さらに微小変化項をレーザー発振など様々な自然現象のメカニズムを解明するため平均近似することで得られる連立方程式を解くことで発振器のに非常に重要であり、多くの研究成果が報告されている。また、定常状態の振幅・位相差を求める手法である。それに対して本同様に電気電子の分野でも同期現象に関する様々な研究が行わ研究では正規化後の連立方程式に定常発振状態を仮定した変数れている [1]-[3]。弱非線形システムにおいては、平均化法 [4]、を与え、0〜2π の平均化は定積分ではなく台形公式を用いて近摂動法や漸近法など、いくつかの近似解法を用いて平衡点が求似計算を行い回路で実現する。さらにニュートンホモトピー法

められる。遠藤らは、van der Pol 発振器を梯子上、環状、二を用いて、回路を解析することにより複数平衡点を 1 度のシ次元配列状に結合したシステムに対して平均化法を用いた理論ミュレーションで探索することが可能である。ニュートンホモ解析と回路実験結果を報告している [5]-[7]。しかし、平均法はトピー法は階曲線追跡回路を用いることにより、簡単に SPICE結合する発振機の数が大きくなるほど、平衡点の数は大きくな上で実現することが可能である。

なり解析は容易ではなくなる。
本研究では、平均化法を SPICE に適用し、ニュートンホモトピー法を用いて解析することにより、発振器の解析における2. 平 均 化 法
ここでは平均化法について説明する。平均化法は次のような—1—である。定常解を求める場合、振幅・位相は一定になるので2 階の非線型微分方程式に用いられることが多い。
ẍ + x = εf (t, x, ẋ)(1)

上式において右辺は非線形関数で 0 < ε ≪ 1 は微小な定数とする。式 (1) は次のような 1 階の連立方程式で表すことができる。

{ẋ = yẏ = εf (t, x, ẋ) − x(2){3. ニュートンホモトピー法
ホモトピー法はニュートン法に比べて大域的な収束性を示し、は代数方程式に補助変数 τ を導入し、初期値 (V(0) , τ = 0) から連続的に方程式を変化し、指定した値 (τ = 1) に達したとき、

ẋ = yẏ = −x(3)式 (3) は線形微分方程式の計算により解は次のようになる。

{られる。

複数解を求めることができるという特徴を持つ。ホモトピー法

ここで ε は微小定数なので ε = 0 と仮定すると次のような解を得ることができる。

ρ̇(t) = 0、 θ̇(t) = 0 とおくことができる。また、非線形関数f (θ, ρ, ϕK ) は SPICE のアナログビへービアモデルを用いて得元の代数方程式に帰着するように変形して解析する方法である。
またホモトピー法の中でもニュートンホモトピー法は SPICEで容易に実現が可能であるため、今回はこれを用いる。二ュートンホモトピー法は次の式で表される。

x = ρ(t) sin(t + θ(t))
y = ρ(t) cos(t + θ(t))(4)

G(V , τ ) = g(V ) − (1 − τ )g(V(0) ) = 0(10)式 (4) を式 (2) に代入し整理することにより、次の式が得ら

ここで、G(V , τ ) は補助方程式と呼ばれ、(V , τ ) = (V(0) , 0)れる。のとき G(V(0) , 0) = 0 を満足し、τ = 1 に達したとき解が求ま ρ̇(t) = εf (t, ρ(t) sin(t + θ(t)), ρ(t) cos(t + θ(t)))
· cos(t + θ(t))εθ̇(t)= − f (t, ρ(t) sin(t + θ(t)), ρ(t) cos(t + θ(t)))
ρ
· sin(t + θ(t))(5)

る。図 1 にニュートンホモトピー法の回路モデルを示す。

v0(τ- 1) g 0(0)VCCS
g 0(V,τ )v1
g 1(V,τ )・・・それぞれ定数 ρ と θ として t = 0, 2π の区間で平均化した次の(τ- 1) g 1(0)VCCS

g 1(0)

1[Ω]平均化法の理論では、ε = εr（小さい定数）が存在し 0 < ε < εrを満たす時、ρ と θ に関する上式の解は右辺の ρ(t) と θ(t) をg 0(0)1[Ω]・・・vM(τ- 1) gM(0)VCCS方程式で近似することができる。
∫ 2περ̇(t)=f (ϕ − θ, ρ sin ϕ, ρ cos ϕ) · cos ϕdϕ


2π 0

∫ 2πεf (ϕ − θ, ρ sin ϕ, ρ cos ϕ) · sin ϕdϕθ̇(t) = −

2πρ
0


ϕ≡t+θ

gM(V,τ )gM (0)

1[Ω]

τSTC

Σ
-1

図 1 ニュートンホモトピー法の回路モデル(6)
本研究では、式 (6) を SPICE 上で直流回路として実現するために、定積分項を台形公式を用いて近似する。

ρ̇(t) =

θ̇(t) =ε(f (θ, ρ, ϕ0 ) + f (θ, ρ, ϕK ))2πε+ (f (θ, ρ, ϕ1 ) cos ϕ1 + f (θ, ρ, ϕ2 ) cos ϕ22π+ . . . + f (θ, ρ, ϕK−1 ) cos ϕK−1 )ε(f (θ, ρ, ϕ0 ) + f (θ, ρ, ϕK ))2πρε(f (θ, ρ, ϕ1 ) sin ϕ1 + f (θ, ρ, ϕ2 ) sin ϕ2+2πρ
+ . . . + f (θ, ρ, ϕK−1 ) sin ϕK−1 )解曲線追跡回路

ニュートンホモトピー法を SPICE で実現する際、補助変数 τを導入するために図 2 に示す解曲線追跡回路が用いられる。解曲線追跡回路とは弧長法を元に考えられた回路である。(n+1)(8)

次元のユークリッド空間 x において、微小変化は式 (11) で表すことができる。

ここで21· 2π, ϕ2 =· 2π,ϕ0 = 0, ϕ1 =KK
K −1. . . , ϕK−1 =· 2π, ϕK = 2π
K

図2(7)

ds =(9)√(dx1 )2 + (dx2 )2 + (dx3 )2 + · · · + (dxn+1 )2(11)ここで弧長 s を時間変数 t、各 x 成分を電圧と仮定すると式—2—

(12) が得られ、図 3 の回路で表すことができる。
)2M (∑dvi
i=1ds(+dvτds

式 (16) はベクトル微分方程式で表すと次のようになる。

)2=1ẍ + Bx = εẋ −(12)1εẋc3(17)このとき、[式 (12) において vi (i = 1, 2, . . . , M ), vτ は、それぞれ式 (12)の vx (x = 1, 2, . . . , n), vx (x = n + 1) に対応する。図 2 においT

33 Tx = [xa , xb ] , xc = [xa , xb ] , B =て、入力は電圧制御電流源 (VCCS) によって決定される。

4. シミュレーション例
図 3 に示す 2 つの van der Pol 発振器がインダクタで結合されたモデルに対し、第 2 章、3 章で説明した手法を用いて解析を行う。

1−γ−γ1](18)である。式 (17) に線形変換を行うと次のようになる。
 ÿ + P −1 BP y = εẏ − 1 εP −1 εẋc3(19) x = Py
上式においてB の固有値は λ1 = 1 − γ 、λ2 = 1 + γ となる。し√√たがって、λ1 に対する固有ベクトルを p1 = [1/ 2, 1/ 2]T 、√√
λ2 に対する固有ベクトルを p2 = [1/ 2, −1/ 2]T とすると、P −1 BP は対角化され次のようになる。
[λ1−1TP BP = P BP =0
図 3 2 つの van der Pol 発振器を{(21)このときωa 2 = λa , ωb 2 = λb f1 (ya , yb , ẏa , ẏb ) = ẏa − 1 g1 (ya , yb , ẏa , ẏb )31 f2 (ya , yb , ẏa , ẏb ) = ẏb − 3 g2 (ya , yb , ẏa , ẏb )(22)である。式 (21) において g1 , g2 は次式で定義される。

[(14)g1]T= P −1 xc =g2d(P T xc )dt)3  ( 1111√√√√
y+yabd 22 ) 2 = 123 1  ( 1dt1√−√√ ya − √ yb
2222


3 232
 2 ya ẏa + 2 ẏa yb + 3ya yb ẏb = (23)
3 2
32yb ẏb + ẏb ya + 3yb ya ẏa22

ここで、(vi , t) を次のように定義する。

本例題では γa = γb = γ 、εa = εb = ε とする。

ÿa − ωa 2 ya = εf1 (ya , yb , ẏa , ẏb )ÿb − ωb 2 yb = εf2 (ya , yb , ẏa , ẏb )式 (13) を変形すると次のようになる。

式 (13) は式 (16) のように正規化することができる。
2 x¨a − εa (1 − xa )x˙a + xa − γa xbx¨ − εb (1 − xb 2 )x˙b + xb − γb xa
 bLiα, γi =εi = √Li + L011+
Ci LiCi L0
 i = a, b(20)λ2表すことができる。

ここで図 3 における a, b 点の接点方程式は次のようになる。
∫∫

1dva
1
vdt+C+i=(va − vb )dt

aara

La ∫dtL0∫

dvb1
 1
vb dt + Cb+ irb =(vb − va )dt(13)
Ldt
L0b


3

ira = −αva + βva
irb = −αvb + βvb 3√αxi , i = a, bvi =
3βt′t= √11
+Ci LiCi L0]

このような対角化により、式 (19) は次のような微分方程式でインダクタで結合した結合発振器

 2()
d va
α3β dva−1−
2C(aαdt ) dt
111++va −
vb = 0

Ca La
Ca L0
Ca L0()
α3β dvb d2 vb
−1−
dt2C(αdt )b



111


++vb −
va = 0

Cb Lb
Cb L0
Cb L 00(15)すなわち[f1f2(16)]11 22 ẏ1 − 2 ya ẏa − 2 ẏa yb − ya yb ẏb =
11 22ẏ2 − yb ẏb − ẏb ya − yb ya ẏa22

(24)である ε が十分に小さいとき (≃ 0)、 ya , ẏa , yb , ẏb に関して次の関係式が得られる。

{ya = ρa sin(ωa t + θa ), ẏa = ρa ωa cos(ωa t + θa )yb = ρb sin(ωb t + θb ), ẏb = ρb ωb cos(ωb t + θb )

(25)

—3—

式 (24) に式 (25) を代入し、次の式が得えられる。

∫ 2π
f (ya , yb , ẏa , ẏb ) · cos(ωi t + θi )ερ̇dϕi (t) =
2πωi
0∫ 2πf (ya , yb , ẏa , ẏb ) · sin(ωi t + θi )εθ̇i (t) =dϕ2πω i ρi0
i = a, b(26)xa = −xb ≃ 2.05 sin(ωb t + θb )(29)√となり、ここで ωb = 1 + γ = 1.00049988 である。式 (28) は同相解を示し、式 (29) は逆相解を示す。.
例題回路では、3 つの平衡点を SPICE を用いて 1 度の解析で求めることができた。われわれの手法を拡大することで、複数個結合した発振器の平衡点をより簡単に解析することができると考えられる。

我々の提案するアルゴリズム中では式 (26) における定積分は台形公式を用いて表され、アナログ回路で実現される。また、台6. ま と め形公式を用いる際に行う積分範囲の分割数は K = 20 とする。

本研究では平均化法にニュートンホモトピー法を導入し、この回路に STC を接続し、ニュートンホモトピー法を用いてSPICE を用いて複数平衡点を簡単に解析するアルゴリズムを解析する。

提案した。接点方程式に平均化法を用いる際、定積分項を台形公式を用いて近似し直流回路で実現した。これにより、SPICE

5. シミュレーション結果上で容易に実現可能なニュートンホモトピー法を用いて解析することができた。例題として 2 つの van der Pol 発振器をインダクタで結合した回路を解析した。我々の提案アルゴリズムにより、3 つの平衡点を 1 度のシミュレーションで解析することができた。
本研究では正規化後の式を SPICE 上で回路として実現し、安定判別は MATLAB を用いた。今後の課題として、解析回路の拡大や安定判別を含めた SPICE による解析があげられる。
図4

謝

解析結果 (γ = 0.001).

図 4 は γ = 0.001 かつ ε = 0.1 における解析結果を示す。図4 において τ = 1 を満たす点を ×, ◦, △ とし、それぞれの値を表 1 に示す。 表 2 は 3 つの平衡点における固有値を示す。
表 1 γ = 0.001 における解time[s]ρaρb1.7359

2.8977741.899µ4.32081.4057

1.57196.63212.6675m

2.8953

表 2 3 つの平衡点における固有値

time[s]λ1

λ21.7359 -17.2000 -8.7900
4.3208-11.9584

5.72846.6321 -17.1001 -8.7699

これらは MATLAB による計算で得られたものである。すべての固有値 λ の実部が 0 より小さいとき、その平衡点は漸近安定となるので、time= 1.7359, 6.6321 のとき漸近安定となる.
今回 x, y は 次のように定義した。

1√

x = P y, P =  12√2
1√2 1 −√
2(27)すなわち、time= 1.7359 のとき定常解は

xa = xb ≃ 2.05 sin(ωa t + θa )となり、ここで ωa =√(28)

1 − γ = 0.99949987 である。また、

time= 6.6321 のとき定常解は辞

本 研 究 の 一 部 は 、日 本 学 術 振 興 会 の 科 学 研 究 費 補 助 金(22500203) により得られた成果である。
文

献[1] T. Kousaka, H. Kawakami, and T. Ueta, “Synchronizationof Electric Fireﬂies by Using Square Wave Generators,” IEICE Trans. Fundamentals, vol. E81-A, pp. 656-663, 1998.
[2] H. Kitajima and H. Kawakami, “Synchronized State Observed in Coupled Four Oscillators,” Proc. NOLTA’ 02,Xi’an, PRC, Oct. 7-11, 2002.
[3] X. Ge, M. Arcak, and K. N. Salama, “Nonlinear Analysisof Cross-Coupled Oscillator Circuits,” IEEE Conference onDecision and Control, Cancun, Mexico, Dec. 9-11, 2008.
[4] A.H. Nayfeh and D.T. Mook, “Nonlinear Oscillations,”
Wiley-Inter Science, 1998.
[5] T. Endo and S. Mori, “Mode Analysis of a Multimode Ladder Oscillator,” IEEE Trans. Circuits and Syst., vol. 23, pp.
100-113, Feb. 1976.
[6] T. Endo and S. Mori, “Mode Analysis of Two -Dimensional
Low-Pass Multimode Oscillator,” IEEE Trans. Circuits andSyst., vol. 23, pp. 517-530, Sep. 1976.
[7] T. Endo and S. Mori, “Mode Analysis of a Ring of a
Large Number of Mutually Coupled van der Pol Oscillators,” IEEE Trans. Circuits and Syst., vol. 25, pp. 7-18,
Jan. 1978.
[8] E. Ikeno and A. Ushida, “The Arc-length Method for theComputation of Characteristic Curves,” IEEE Trans. Circuits Syst., vol. 23, pp. 181-183, 1976.
[9] A. Ushida, Y. Yamagami, I. Kinouchi, Y. Nishio and Y. Inoue, “An Eﬃcient Algorithm for Finding Multiple DC Solutions Based on the SPICE-Oriented Newton Homotopy
Method,” IEEE Trans. Computer-Aided Design of Integrated Circuits Syst., vol. 21, no. 3, pp. 337-348, Mar. 2002.
[10] Y. Inoue, “DC Analysis of Nonlinear Circuits Using
Solution-Tracing Circuits,” IEICE Trans. vol. J74-A, pp.
1647-1655, 1991.

—4—