やりたいこと

JuliaでSPICEみたいな回路シミュレーションをしてみます。今回はタイトルの通り、RLC直列回路を修正節点法で求めて、一次の後退差分法で演算します。また、修正節点法との比較として閉路解析法でも同じ回路を後退差分法で演算します。
今回の特徴としては、SymPyを使って代数をそのまま計算したところで、Julia+SymPyだと式通りの計算ができるのでなかなか面白いかなと思います。

参考図書

[1] 書籍：電子回路シミュレーション
[2] 論文：パワーエレクトロニクスを対象としたリアルタイム・シミュレータに関する研究

上記の２つの参考図書は、SPICEを理解するためには非常に良いと思います。修正節点法に関しては、基本的には[2]の48ページ目くらいに記載の通り計算をしました。また、[1]も非常にわかりやすい書籍だと思います。閉路解析法は、適当に直感で作りました。

閉路解析法の場合

V=RI+L\frac{dI}{dt}+\frac{1}{C}\int Idt\\
V=RI(n)+L\frac{I(n)-I(n-1)}{\Delta t}+V_C(n-1)+\frac{I(n)}{C} \Delta t \\
I(n) = \frac{V+L\frac{I(n-1)}{\Delta t}-V_C(n-1)}{R+\frac{L}{\Delta t}+\frac{\Delta t}{C}} \\
V_C(n) = V_C(n-1) + \frac{I(n)}{C} \Delta t

最初の行が立式。2行目が差分法への展開。最後の2行が計算で使用する式です。今回の例だと、閉路は一つになるので結構簡単な式でいけます。また、（個人的には？）直感的な表現なので、式を考えるときに便利だと感じています。一方で、SPICEのようにネットリストから機械的に式を生成するためには、どのように式を作ればよいのか、判断が難しくなるそうです。

using Plots
gr()

# パラメータ
ts = 1e-3
tl = 0.0:ts:1.0
Rset, C, L = 10.0, 1e-1, 1e-1
G1 = 1/Rset
fin = 10 # Hz

ix, Vc = 0.0, 0.0
vina, vca, ixa = [], [], []
for t in tl
    Vin = 100*sin(2pi*fin*t)
    global ix = (Vin - Vc + L/ts*ix)/(L/ts + Rset)
    global Vc = Vc + ix/C * ts    
    push!(vina, Vin);    push!(vca, Vc);    push!(ixa, ix)
end

# 描画
pl1 = plot(tl, Vector[vina, vca], label=["Vin" "Vc"])
pl2 = plot(ixa, label=["i"])
plot(pl1, pl2, layout=(2,1))

修正節点法の場合

修正節点法は、名前の通り節点法を修正した方式になります。節点法と考え方が同じですが、節点法の場合は電流源しか扱えないのに対して電圧源を扱えるのが特徴です。節点法ですべてやりきるためには、電圧源を全て電流源に置き換える必要があり、いろいろと大変なので、SPICEでは修正節点法が使用されているそうです。ちなみに、修正節点法では、LとCは離散化モデルとしてコンダクタンス（抵抗）と電流源と電圧源に置き換えることになります。式の導出は参考図書[2]の48ページ目にお任せして、SymPyを使って計算をすると以下のようになります。

using SymPy, Plots
gr()

@syms G_1 G_C R_L E_in I_C E_L 

# SymPyで式を作る
A = [G_1 -G_1    0     0  1
    -G_1  G_1    0     1  0
       0    0  G_C    -1  0
       0    1   -1  -R_L  0
       1    0    0     0  0]

funcx = inv(A) * [0; 0; -I_C; E_L; E_in]
Ein_a, V3_a, I3_a = [], [], []

# 計算を実施
I3, V3 = 0.0, 0.0
for t in tl
    global I3, V3
    Gc, JC = C/ts, -C/ts * V3
    RL, EL = L/ts, -L/ts * I3
    Ein = 100*sin(2pi*fin*t)

    xxxx = subs(funcx, G_1=>G1, G_C=>Gc, R_L=>RL, E_in=>Ein, E_L=>EL, I_C=>JC)
    V3 = convert(Float64, xxxx[3])
    I3 = convert(Float64, xxxx[4])
        
    push!(Ein_a, Ein);     push!(V3_a, V3);    push!(I3_a, I3)
end
pl1 = plot([Ein_a, V3_a], label=["Vin" "Vc"])
pl2 = plot(I3_a, label="i")
plot(pl1, pl2, layout=(2,1))

jupyter notebookを使えば、SymPyはきれいな数式が表示されてとても気持ち良いです。計算は遅いので、lambdifyを使えば、若干速くなるかもしれません。jupyter notebookは後ほどリンクします。

最後に

時間が全然足りなくて色々と力及ばずな内容になりましたが、参考図書は非常に素晴らしいものなので、その紹介ができた点は良かったかなと思います。