RAKUS Developers Blog | ラクス エンジニアブログ

株式会社ラクスのITエンジニアによる技術ブログです。

Python初心者がSymPyで電気回路を解いてみた

はじめに

こんにちは、crowd_kです。

プログラムを本格的に触り始めて、1年弱が経ちました。

というのも、私は「電気電子学科」出身であるためプログラミングは授業の一環でほんの少し触った程度しかありませんでした。
しかし、そんな授業で"Python"を少しだけ学んだことがあります。
内容は、"Python"を使って確率統計を解くというものでした。
紙を使って手計算すると、1問とくのに何十分もかけ、A4 の紙を1枚使い切ったり と大変だったものが"Python"を使って解くと簡単にあっさりと解けてしまい、とても驚いた記憶があります。

当時、「確率統計よりも厄介な(と感じる)、電気回路や電磁気を"Python"を使って解くことができたら便利ではないか」と思ったことがあります。
しかし、就活や卒業研究・論文 等で忙しかったため、実際にトライすることはありませんでした。

ふとそんなことを思い出し、"Python"の勉強になるのではと思ったので
簡単な電気回路を"Python"を使って解いて記事にまとめてみました。

Pythonの環境

まずは、Pythonのインストールです。
「これからPythonを勉強するなら3系が良い」という記事を多く見たので、3系を選択しました。
また、当時Pythonjupyter notebookで動かしていたので、そちらもあわせてインストール、
また同時に、sympyもインストールしました。

pip install jupyter --user
pip install sympy --user

コマンドプロンプトjupyter notebookと入力すれば、ブラウザでjupyter notebookが立ち上がります。
これでひとまず、準備は完了です。

実際に解いてみる

まずはじめに、先ほどインストールしたsympyをインポートしておきます。

import sympy as sym

今回はRL回路とRC回路の過渡応答について解いてみたいと思います。
RL回路とは、抵抗(R)とコイル(L)でできた電気回路のことで、RC回路は抵抗とコンデンサー(C)でできた回路のことです。

RL回路

まずは、変数を定義します。
変数の宣言は、sympy.symbolsで下のように定義できるようです。

R, L, t, V = sym.symbols("R L t V")

ここで、R=抵抗[Ω]、L=コイル[H]、t=時間[s]、V=電圧[V] を表します。

また、sympy.Function関数も定義することができるようです。

i = sym.Function('i')

i=電流[A]です。
上記は、電流iは時間tによって変化することを表します。
f:id:crowd_k:20200121143625p:plain

上記の回路を解いていこうと思います。
一般的に以下のような式で表されます。

f:id:crowd_k:20200121144638p:plain

これを"python"を使って表していきたいと思います。
抵抗の電圧を表す式は単純で、以下のようになります。(すごく直感的に書くことができるので、楽です)

R*i(t) 

次に、コイルでの電圧ですが、Lと電流i(t)を時間tで微分したものの積で表されています。
sympyにおいて微分は、sympy.diff()と表すことができるので、

L*i(t).diff(t,1)

となります。
これらを使うと、RL回路の回路方程式は以下で表すことができることがわかります。

eq1 = sym.Eq(R*i(t) + L*i(t).diff(t,1), V)

eq1について出力してみたら、下のようにうまくいきました。

f:id:crowd_k:20200121145953p:plain

今回は、R = 100、L = 100, V = 100000 として、解きたいと思います。(値は適当です)

eq1 = sym.Eq(100*i(t) + 100*i(t).diff(t,1), 100000)

方程式を解くには、sympy.solveを使います。
また、微分方程式を解く上で必要な初期値を引数として渡すことで、それを考慮した計算をしてくれるみたいです。(非常に便利)
今回は、t = 0のとき、(回路の電源をONにする直前)の電流を0とします。

ans1 = sym.dsolve(eq1, ics = {i(0):0})

結果を出力すると、以下の式が得られました

f:id:crowd_k:20200121152636p:plain

これで、解き終わりました。
これが何を表しているのか、グラフに書いてわかりやすくしようと思います。

from sympy.plotting import plot

で、plotを使って、表示してみます。

plot(1000 - 1000*exp(-t),(t,0,10))

f:id:crowd_k:20200121153054p:plain

上のような結果が得られました。
i(t)は時間経過とともに徐々に1000[A]に近づいている(すぐには立ち上がらない)のがわかります。
抵抗のみなら、電源をONにしたと同時に電気が流れますが、コイルを挟むことで、ゆっくり立ち上がるようになります。
これは、コイルに流れる電流の変化したとき、その向きに対して逆向きの力を加えようとする性質からで、正しい答えが導き出せたと思います。

RC回路

f:id:crowd_k:20200121174633p:plain

同様に、RC回路を解いてみようと思います。
まずは、変数と関数の宣言です。

R, C, t, V = sym.symbols("R C t V")

i = sym.Function('i')

R=抵抗[Ω]、C=コンデンサー[F]、t=時間[s]、V=電圧[V] としています。

f:id:crowd_k:20200121154213p:plain

RC回路の場合、上記のような式になります。
抵抗の電圧を示す式に関しては、RL回路と同様の考え方ができます。
しかし、コンデンサーはコイルとは異となります。
コンデンサーの電圧は、Cの逆数と、電流i(t)をtで積分したものの積で表すことができます。
sympyでは、積分sympy.integrateと表すことができるので、

eq2 = sym.Eq(R*i(t) + sym.integrate(1/C * i(t), t), V)

で、RC回路の回路方程式を表すことができます。
積分がいると面倒なので、(この状態の方程式を解く方法がわからなかった)ひと工夫して解こうと思います。
それは、i(t)をtで積分したものを新たにq(t)として、解く方法です。

f:id:crowd_k:20200121160633p:plain

裏を返すと、以下のようにもなります。

f:id:crowd_k:20200121160844p:plain

これらを使うと、回路方程式は以下のようになりRL回路と同様に微分を含めた回路方程式になります。

f:id:crowd_k:20200121173515p:plain

これを、コードで書き直したものが下のものになります。

eq2 = sym.Eq(R*q(t).diff(t,1) + 1/C*q(t), V)

これなら、RL回路と同様に方程式を解いてくれそうです。
R = 100、C = 0.01, V = 100 とし、初期条件にt=0のときのqを0に(スイッチを入れる前は電流が流れていないこととし、総電荷量も0になる)を当てはめると

eq2 = sym.Eq(100*q(t).diff(t,1) + 1/0.01*q(t), 100)

ans2 = sym.dsolve(eq2, ics = {q(0):0})

となります。
ans2の実際の出力結果がこちら。

f:id:crowd_k:20200121162407p:plain

この答えは、i(t)をtで積分したq(t)を表しています。
i(t)に戻すために、両辺をtで微分します。

ans2_i = (1.0 - 1.0*exp(-1.0*t)).diff(t, 1))

f:id:crowd_k:20200121163136p:plain

同様に、グラフに表してみたいと思います。

from sympy.plotting import plot

plot(ans2_i, (t,0,10))

f:id:crowd_k:20200121163429p:plain

RC回路における電流i(t)の時間変化を表すグラフです。
立ち上がりは抵抗のみのときのように一瞬で立ち上がりますが、徐々に衰退しているのがわかります。
RL回路とは逆に時間経過とともに電流は流れなくなります。
コンデンサーは電荷をため込む性質(充電バッテリーのように)があります。
時間経過とともに電気をため込み、満タンになると電気が通らなくなります。
それがグラフで読み取ることができたので、こちらも正しい結果であると思います。

さいごに

本当ならこの流れでRLC回路も解いてみたかったのですが、解くことができませんでした。

f:id:crowd_k:20200121165244p:plain

積分微分が両方混ざった回路方程式になります。
なので、「積分を含んだ方程式を解く」または「二回微分を含む方程式を解く」
ことができる必要があります。
また、ベクトルの内積外積、線積分や面積分といった計算が中心になる電磁気も、解けるようになったら面白そうです。

まだまだ勉強不足だと思うので、試行錯誤して今後チャレンジしていきたいです。


  • エンジニア中途採用サイト
    ラクスでは、エンジニア・デザイナーの中途採用を積極的に行っております!
    ご興味ありましたら是非ご確認をお願いします。
    20210916153018
    https://career-recruit.rakus.co.jp/career_engineer/

  • カジュアル面談お申込みフォーム
    どの職種に応募すれば良いかわからないという方は、カジュアル面談も随時行っております。
    以下フォームよりお申込みください。
    forms.gle

  • イベント情報
    会社の雰囲気を知りたい方は、毎週開催しているイベントにご参加ください! rakus.connpass.com

Copyright © RAKUS Co., Ltd. All rights reserved.