はじめに
RLC回路がもしも直列接続だけで構成されていたら、過渡解析は比較的扱いやすいといえる。というのは、分岐がなく電流はどこを計測しても同じ値を示すからである。しかし、並列の場合は分岐が存在しているので、それぞれの枝に分担される電流は異なる。ただ、それらの和は最初の分岐前のときのものと等しくなるというキルヒホッフの第一法則の情報しかわからない。そこで、今回は並列回路を含む直列回路RLC回路の過渡応答を扱う。そのためにはまず回路条件より微分方程式を導出する。次に、現代制御理論によりそれを状態方程式に直す。最後に、数値解析の一種であるオイラー法よりも精度が高いルンゲクッタを用いることで微分方程式、つまり状態方程式を解きシミュレーションを行う。
回路設定
回路設定としては、以下のような回路モデルで考える。
まず、コンデンサ$C_1$には、初期電荷として$Q$だけ蓄えられているが$C_2$には電荷は蓄えられていないものとする。また、$t=0$時点において、スイッチを投入し、コンデンサ1,2の電荷である$Q_1,Q_2$の過渡応答をみるものとする。
微分方程式
$C_1,C_2,R$の回路構成から、以下のような微分方程式をキルヒホッフの第二法則（電圧則）より導出することができる。
\frac{Q_1}{C_1}=-R\dot{Q_1}+\frac{Q_2}{C_2}
一方で、$L,C_2$の回路構成から以下のような微分方程式をキルヒホッフの第二法則（電圧則）より導出することができる。
\frac{Q_2}{C_2}=-L(\ddot{Q_1}+\ddot{Q_2})
状態方程式
今回の場合は、未知の状態変数が$Q_1,Q_2$と２つ存在しているので、状態方程式を用いて微分方程式を解くことを考える。
\frac{Q_1}{C_1}=-R\dot{Q_1}+\frac{Q_2}{C_2}
の両辺を微分すると、
\frac{\dot{Q_1}}{C_1}=-R\ddot{Q_1}+\frac{\dot{Q_2}}{C_2}
となるので、
\begin{pmatrix}
\dot{Q_1} \\
\dot{Q_2} \\
\ddot{ Q_1} \\
\ddot{ Q_2} \\
\end{pmatrix}
=
\begin{pmatrix}
\frac{1}{RC_1} & -\frac{1}{RC_2} &0&0 \\
0 & 0 &0&1\\
0&-\frac{1}{C_2L}&\frac{1}{RC_1}&-\frac{1}{RC_2}\\
0&0&-\frac{1}{RC_1}&\frac{1}{RC_2}\\
\end{pmatrix}
\begin{pmatrix}
Q_1 \\
Q_2\\
\dot{Q_1}\\
\dot{Q_2}\\
\end{pmatrix}
数値解析
微分方程式を数値解析しシミュレーションを行う場合、一般的にまずは微小区間に時間を区切りその中で漸化式のように微小時間経過後の状態と現在の状態を比較するという方法が取られる。しかし、その場合のデメリットとして微小時間を無限小にできない以上どうしてもそれらを足し交わせる（積分）するときに誤差が生じてしまうことである。そこで、今回はそのような誤差をなるべく少なくするために、ルンゲクッタ法を用いるものとする。ルンゲクッタ法のPythonプログラムの例は以下の記事の『例.ロジスティック方程式』を参照されたい。
プログラム
初期条件
まず、電荷の条件は、$Q_1(t=0)=Q,Q_2(t=0)=0$である。
したがって、
\frac{Q_1}{C_1}=-R\dot{Q_1}+\frac{Q_2}{C_2}
より、
\dot{Q_1}(t=0)=-\frac{Q}{C_1R}
また、
0=\frac{Q_2(t=0)}{C_2}=-L(\ddot{Q_1}(t=0)+\ddot{Q_2}(t=0))
\ddot{Q_1}(t=0)+\ddot{Q_2}(t=0)=0
となる。
一方で、スイッチ投入前において、リアクトルには電流は流れていなかったので、
\dot{Q_1}(t=0)+\dot{Q_2}(t=0)=0
となる。
プログラム
それでは、プログラムを示す。
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import japanize_matplotlib
t=0
delta_t=1.0*10**-6
R=10
L=1.0*10**-3
C_1=10*10**-5
C_2=10*10**-5
#C_!の初期電圧V
V=1.0*10**2
Q=C_1*V
Q_1=Q
Q_2=0
#データを入れる箱を作る
t_ary=[]
Q_1_ary=[]
Q_2_ary=[]
dQ_1_ary=[]
dQ_2_ary=[]
#状態方程式
def v(x,t):
return np.array([[-1/(R*C_1),1/(R*C_2),0,0],[0,0,0,1],[0,0,-1/(R*C_1),1/(R*C_2)],[0,-1/(C_2*L),1/(C_1*R),-1/(C_2*R)]])@x
#メインプログラム
while t<1.0*10**-1:
if 'dQ_1'and'dQ_2' in locals():
#Runge-Kutta法
k1 = delta_t*v(x, t)
k2 = delta_t*v(x+0.5*k1, t+0.5*delta_t)
k3 = delta_t*v(x+0.5*k2, t+0.5*delta_t)
k4 = delta_t*v(x+k3, t+delta_t)
x = x + (k1 + 2*k2 + 2*k3 + k4)/6
t=t+delta_t
t_ary.append(t)
#箱に値を現在の値を入れる
Q_1_ary.append(x[0])
Q_2_ary.append(x[1])
dQ_1_ary.append(x[2])
dQ_2_ary.append(x[3])
else:
#初期条件
dQ_1=-(Q_1)/(R*C_1)
dQ_2=-dQ_1
x=np.array([Q_1,Q_2,dQ_1,dQ_2])
#グラフの箱を作る
fig = plt.figure(figsize=(8, 5))
ax1 = fig.add_subplot(221)
ax2 = fig.add_subplot(222)
ax3 = fig.add_subplot(223)
ax4 = fig.add_subplot(224)
#グラフの図示
plt.subplots_adjust(wspace=0.5, hspace=0.5)
ax1.set_title('Q_1')
ax1.set_xlabel('time[s]')
ax1.set_ylabel('Q_1[C]')
ax1.plot(t_ary,Q_1_ary)
ax2.set_title('Q_2')
ax2.set_xlabel('time[s]')
ax2.set_ylabel('Q_2[C]')
ax2.plot(t_ary,Q_2_ary)
ax3.set_title('dQ_1/dt')
ax3.set_xlabel('time[s]')
ax3.set_ylabel('dQ_1[A]')
ax3.plot(t_ary,dQ_1_ary)
ax4.set_title('dQ_2/dt')
ax4.set_xlabel('time[s]')
ax4.set_ylabel('dQ_2[A]')
ax4.plot(t_ary,dQ_2_ary)
plt.savefig("parallel_RLC.png")
plt.show()
結果
このプログラムを実行した結果は以下のようになった。
このように、最初は$C_1$がとてつもない勢いで放電されて一部は$C_2$に充電される。しかし、その間にリアクトルがあるので、電荷は振動を行い減少してしまう。
まとめ
今回は、複雑なRLC回路の過渡解析をPythonを用いた数値解析を用いて実行することができるかを調査した。そこで、解法をより汎用的かつ多変数でも扱えるようにするため、微分方程式を一旦現代制御の状態方程式にすることを考えた。また、そのような式を4次のランゲクッタ法を用いて数値解析を行った。結果、振動減衰を再現できていることからかなり高い精度で過渡解析を行うことが可能であるということが分かった。