この記事はJulia Advent Calendar 2017の21日目の記事です。
やりたいこと
SymPyを使って、回路インピーダンスの周波数特性のボード線図を描いてみます。
JuliaにはSymPy.jlというパッケージがあり、Pythonの代数計算ライブラリSymPyを簡単に使うことができます。SymPyをPythonではなくJuliaで実行するメリットは今のところわかりませんが、Juliaで計算した結果と一緒に扱いたい場合には、SymPy.jlを使うと便利だと思っています。
今回の例は簡単な内容ですが、もうちょっと式が複雑になると、SymPyのLambdifyが非常に便利なので、Lambdifyは便利だよ、という意味もあります。
また、実行結果はGitHubのJupyter notebookを参照ください。
RLC直列共振回路
今回考える回路は、RLC直列共振回路です。回路図は以下のURLにあるような回路です。
http://wakariyasui.sakura.ne.jp/p/elec/koukairo/kyousinn.html
RLCに適当にパラメータを与えたものをPlotするのは、以下のようなソースになります。また、回路の並列接続の計算は関数化すると楽なので、とりあえず関数にしています。
using SymPy, Plots
parallel(z) = 1/sum(map(x->1/x, z))
serial(z) = sum(z)
@syms L1 R1 C1 omega s z_inv z_inv2 T
ZR1 = R1
ZL1 = s * L1
ZC1 = 1/(s * C1)
r1, l1, c1 = 1, 1e-3, 1e-6
eq1 = serial([ZR1, ZL1, ZC1])
eq2 = simplify(subs(eq1, (R1, r1), (L1, l1), (C1, c1), (s, im * omega)))
eq3 = lambdify(eq2)
x = logspace(-2,10,1000)
p1 = plot(x, 20*log.(abs.(eq3(2pi*x))), xscale=:log10, label="gain")
p2 = plot(x, rad2deg.(angle.(eq3(2pi*x))), xscale=:log10, label="angle", ylims=(-90, 90))
plot(p1, p2, layout=(2,1))
こんな感じで、Lambdifyを使うとSymPyで導出した式を関数化することができ、計算を高速化することができます（どれくらい早いかはちゃんと調べてませんが）。
RLC直列共振＋一部並列共振回路
先の、RLC直列共振回路のコンデンサ部分に並列で、更にRLC直列回路が接続された場合を考えます。
その場合、先ほどのコードのeq1が以下のようになります。
eq1 = serial([ZR1, ZL1, parallel([ZC1, serial([ZR1, ZL1, ZC1])])])
eq1以外同じコードにして、実行すればOKです。
共振周波数が２つになります。こんな感じで直並列の関係が増えてくると、自分で計算するのはミスが増えたり面倒だったりするので、普通は回路シミュレータの周波数特性解析機能を使ったりするんじゃないかなと思います。しかし、SymPyを使うと簡単に計算ができるので、ちょっとした計算ならばSymPyで十分じゃないかと思っています（回路シミュレータの方が楽な気もしますが）。
更に、入力信号電圧波形をFFTして周波数特性に変換して求めたインピーダンスに突っ込んでやれば、流れる電流の周波数特性がわかります。この辺りをやる場合は、わざわざ回路シミュレータを使用するよりも、Juliaで計算した方が楽なような気もします。
ちなみに、JuliaのControlSystemsパッケージがあるため、ボード線図はそのあたりを使えば簡単に書けますが、個人的には、ControlSystemsの描画機能を気に入っていないので、logspaceで自分で描く方が自由に調整できて良いと思っています。
自己紹介
メーカーで、電気関係の仕事をやっています。
この業界では有償ツールを使用して、目的達成のために近道をするのが一般的なのですが、近道をし過ぎると原理が疎かになるため、できる限りは自力で計算式を入力してみるのが良いように思っています。そっちの方が、長中期的にみると目的達成の近道のようにも感じています。
というのは言い訳で、勤務時間中にJuliaで遊びたいだけだという気もしています。