はじめに
RLC回路の過渡応答解析は一般的には微分方程式を用いて考察する。しかし、その微分方程式を解析的に解くのは難易度が高く、ましては応答をグラフ化するのは手計算では限界がある。そこで、微分方程式をプログラムを用いて解くことでRLC回路の過渡応答の様子を伺うことを本記事では試みる。具体的には、微分方程式を状態方程式に変換した上で、簡易な差分法を用いて解析し、代数的な解と比較する。以下は解析結果のグラフである。
問題設定
以下のようなRLC直列回路に正弦波電圧を加えたときの電流の応答を調べる。ただし、R=2[Ω],L=1[mH],C=10[μF]として、電源周波数を50[Hz]、電源電圧を1[V]とした。
微分方程式と状態方程式
RLC直列回路の微分方程式は以下のようになる。ただし、電源電圧とコンデンサに蓄積された電荷の瞬時値をそれぞれ$e(t),q(t)$とした。
e(t)=R\frac{dq(t)}{dt}+L\frac{d^2 q(t)}{dt^2}+\frac{q(t)}{C}
ゆえに、
\frac{d^2 q(t)}{dt^2}=-\frac{q(t)}{LC}-\frac{R \frac{dq(t)}{dt}}{L}+\frac{e(t)}{L}
これを状態方程式にすると以下のようになる。
\begin{pmatrix}
\frac{dq(t)}{dt} \\
\frac{d^2 q(t)}{dt^2} \\
\end{pmatrix}
=
\begin{pmatrix}
0 & 1 \\
-\frac{R}{L} & -\frac{1}{CL} \\
\end{pmatrix}
\begin{pmatrix}
q(t) \\
\frac{d q(t)}{dt} \\
\end{pmatrix}
+
\begin{pmatrix}
0 \\
\frac{1}{L} \\
\end{pmatrix}
e(t)
したがって、以下のように表す。
\dot{\textbf{x}}(t)
=
\begin{pmatrix}
0 & 1 \\
-\frac{R}{L} & -\frac{1}{CL} \\
\end{pmatrix}
\textbf{x(t)}
+
\begin{pmatrix}
0 \\
\frac{1}{L} \\
\end{pmatrix}
e(t)
ただし、
\textbf{x}(t)=\begin{pmatrix}
q(t) \\
\frac{d q(t)}{dt} \\
\end{pmatrix}
\dot{\textbf{x}}(t)=\textbf{v}(t)=
\begin{pmatrix}
\frac{dq(t)}{dt} \\
\frac{d^2 q(t)}{dt^2} \\
\end{pmatrix}
とする。
ここで、位置ベクトル、速度ベクトルをそれぞれ$\textbf{x}(t),\dot{\textbf{x}}(t)$として、差分法で微分方程式を解くことを考える。
この条件のもとで、
I(t)=\frac{dq(t)}{dt}
を求める。
アルゴリズム
微分方程式(状態方程式)を解くアルゴリズムは以下の通りである。
ここでは、ある物体がx軸に対して直線運動していると考える。
変位と速度(差分法の原理)
微小区間($t=t$から$t=t+\Delta t$)に速度$v(t)$（変位の微分）と変位$x(t)$が以下のよう与えられたとする。
$x(t+\Delta t)$は微分の定義式から以下のように考えることができる。
これは、$t=t$のときの情報のみを頼りに$t=t+\Delta t$のときの情報を推定することができることを意味している。
\textbf{x}(t+\Delta t)=\textbf{x}(t)+\textbf{v}(t)\Delta t
これを$t=0$から$t=t$まで加算することで$t=t$での位置$\textbf{x}(t)$を求めることができる。
プログラム
さて、以上の考察をもとにして以下のようなプログラムを作成した。
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
# 微分方程式:V =R Q'*L Q''+Q/C
#パラメータ
R = 2 # 抵抗値
L = 1e-3 # インダクタンス
C = 10e-6 # 静電容量
#初期条件 （電荷と電荷の微分量つまり電流について）
x=np.array([0,0])
#初期時間
t=0
#時間と状態量（ここでは電荷と電荷の微分量）
time_ary=[]
X=[]
#微小時間の時間間隔(小さくするほど精度は上がるが計算時間が比例して増大する)
delta_t=0.00001
while(t<0.02):
E=np.sin(2*math.pi*50*t)
v=np.array([[0, 1], [-1/(L*C),-R/L]])@x+np.array([0,1/L])*E
#位置の微小変化
x=x+v*delta_t
#print(v)
#演算結果を記録
time_ary.append(t)
X.append(x[1])
#時間の更新
t=t+delta_t
plt.xlabel("応答時間[s]")
plt.ylabel("電流[A]")
plt.plot(time_ary,X)
#グラフを保存する
plt.savefig("RLC_状態方程式_電流.png")
#グラフの表示
plt.show()
これを実行すると以下のようなグラフが出力される。
まず、コンデンサの初期電荷量は0なので、電流は急激に流れて高速に振動をした。次に、十分に時間が経過すると電源周波数で緩やかな振動をした。
代数的解法
上記の計算の精度がどれほどなものかを検討するために、代数的な解法を以下に載せる。
一般的に、線形微分方程式は過渡解（特殊解）と定常解の線形和で表すことができる。なので、まずは定常解について考える。
定常解
本記事の微分方程式
Esin (\omega t)=R\frac{dq(t)}{dt}+L\frac{d^2 q(t)}{dt^2}+\frac{q(t)}{C}
の定常解を$q_0(t)=A sin(\omega t+\alpha)$とおく。定常解を微分方程式に代入することで$A,\alpha$を求める。三角関数の合成の法則より、以下のようになる。
A=\sqrt{\frac{1}{(-\omega^2 L-\frac{1}{C})^2+(\omega R)^2}},\alpha=-tan^{-1}(\frac{\omega R}{-L \omega^2 +\frac{1}{C}})
これを用いることで、
I_0(t)=A\omega cos(\omega t+ \alpha)
過渡解
次に、過渡解について考える。
微分方程式
0=R\frac{dq(t)}{dt}+L\frac{d^2 q(t)}{dt^2}+\frac{q(t)}{C}
について考えればよい。
従って、
q_1(t)=A_1 e^{\lambda_1 t}+B_1 e^{\lambda_2 t}
と表すことができる。
ただし、$A_1,B_1$は任意定数で初期条件によって決定される。
また、$\lambda_1,\lambda_2$は、特性方程式
L\lambda^2+R\lambda + \frac{1}{C}=0
の解である。
したがって、電流は
I_1(t)=\lambda_1 A_1 e^{\lambda_1 t}+\lambda_2 B_1 e^{\lambda_2 t}
となる。
一般解
以上の考察から、コンデンサに蓄えられる電荷と流れる電流は以下のように表すことができる。
q(t)=q_0(t)+q_1(t)=A sin(\omega t+\alpha)+A_1 e^{\lambda_1 t}+B_1 e^{\lambda_2 t}
I(t)=I_0(t)+I_1(t)=A\omega cos(\omega t+ \alpha)+\lambda_1 A_1 e^{\lambda_1 t}+\lambda_2 B_1 e^{\lambda_2 t}
初期条件
ここで、$q(t=0)=0,I(t=0)=0$という初期条件を用いることで$A_1,A_2$は以下のように表すことができる。
A_1=\frac{A(\lambda_2 sin\alpha -\omega cos\alpha)}{\lambda_1-\lambda_2}
B_1=\frac{-A(\lambda_1 sin\alpha -\omega cos\alpha)}{\lambda_1-\lambda_2}
描写
上記の考察による代数解をグラフ化するために以下のようなプログラムを作成した。
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
# 微分方程式:V =R Q'*L Q''+Q/C
#パラメータ
R = 2 # 抵抗値
L = 1e-3 # インダクタンス
C = 10e-6 # 静電容量
E=1.0
f=50
t=np.linspace(0,0.02,1000)
omega = 2*math.pi*f
# 定常解
A=(1/((-L*omega**2+1/C)**2+(omega*R)**2))**0.5
alpha = -np.arctan(omega*R/(-L*omega**2+1/C))
q0=A*np.sin(omega*t+alpha)
I0=A*np.cos(omega*t+alpha)*omega
#plt.plot(t,q0)
# 過渡解
lambda1= (-R +(R**2-4*L*1/C)**0.5)/(2*L)
lambda2= (-R -(R**2-4*L*1/C)**0.5)/(2*L)
# 初期条件
A1=A*(lambda2*np.sin(alpha)-omega*np.cos(alpha))/(lambda1-lambda2)
B1=-A*(lambda1*np.sin(alpha)-omega*np.cos(alpha))/(lambda1-lambda2)
q1=A1*np.exp(lambda1*t)+B1*np.exp(lambda1*t)
I1=lambda1*A1*np.exp(lambda1*t)+lambda2*B1*np.exp(lambda1*t)
# 一般解
q=q0+q1
I=I0+I1
plt.xlabel("応答時間[s]")
plt.ylabel("電流[A]")
plt.plot(t,I)
#グラフを保存する
plt.savefig("RLC_状態方程式_電流_代数解.png")
#グラフの表示
plt.show()
これを実行すると以下のようなグラフが出力される。
このグラフは、差分法により解析した解とほぼ一致している。したがって、差分法による解法の精度は解析時間内においては十分であると言うことができる。
まとめ
今回は、RLC直列回路の過渡応答について、状態方程式を用いることでプログラミングによりシミュレーションした。結果、解析解と代数解はほぼ一致した。また、振動解を実用計算時間内において十分に再現することが可能であるということが分かった。具体的には、電流値の応答自体は、初期は共振周波数付近の振動をするが、時間が経つに従って電源周波数での緩やかな振動に落ち着いた。