電気回路:一次遅れ系の過渡現象

はじめに

本記事では、RL回路とRC回路といった一次遅れ系についての過渡現象特性をPythonでシミュレーションしたのちに、時定数がそれらの特性把握のカギとなることを確認する。
ただし、Pythonで用いたライブラリは、微分方程式を解くために使用されるsympyを用いた。

RC直列回路について

微分方程式

この微分方程式は右辺が０でないので、特殊解を用いて工夫する必要がある。

この場合、解の１つである特殊解は以下の様になる。

次に、右辺を０にした微分方程式である

について解くと以下の様になる。

 

ゆえに、コンデンサにたまる電荷は以下の様になる。

ここで、q(t=0)=0という境界条件（初期条件）を用いると、

プログラミングコード

これを、pythonで書くと以下の様になる。

%reset -f
import sympy
import numpy as np
import matplotlib.pyplot as plt

# シンボル（変数記号）の定義
q,i = sympy.symbols('q i', cls=sympy.Function)
t,E,R,C,C1 = sympy.symbols('t E R C C1')

# 方程式 eq1 を立てる
eq1 = sympy.Eq( E, R * sympy.Derivative(q(t), t) + q(t)/C )
print("微分方程式の表示")
display(eq1)

eq2 = sympy.dsolve(eq1, q(t))
print("Qについて解く")
display(eq2)

#初期条件の代入
# eq2の左辺を0、右辺にt=0を代入
eq3 = sympy.Eq(0, eq2.rhs.subs(t,0),)  
C1a = sympy.solve(eq3,C1)[0]
print("初期条件の代入")
display(C1a)
#Qの解を表示する
# 初期条件を与えて求めた積分定数C1aを代入
eq4 = eq2.subs(C1,C1a)
print("Qの解の表示")
display(eq4) 

# 式を整理(展開)
eq4 = eq4.expand()
print("Qの解の展開表示")
display(eq4) 

# eq4 の右辺を t 微分
eq5 = sympy.Eq( i(t), sympy.diff(eq4.rhs,t))
print("電荷を微分して電流値を求める")
display( eq5 )

i = eq5.rhs

# 抵抗両端の電圧 VR 
VR = R*i
print("抵抗Rの両端の電圧")
display(VR)

# コンデンサ両端の電圧 VC
VC = E - VR

print("コンデンサの電圧")
display(VC)
display(VC.collect(E)) # 電圧Eで括る 

# パラメータ:E=5[V] (t>0)でのステップ応答　
prm = { E:sympy.Piecewise((0,t=0)), R:50e3, C:40e-6 }

# パラメータを代入
i = i.subs(prm)
VR=VR.subs(prm)
VC=VC.subs(prm)

# sympy オブジェクトを numpyユニバーサル関数に変換
f_VR = sympy.lambdify(t, VR, 'numpy')
f_VC = sympy.lambdify(t, VC, 'numpy')

tt = np.linspace(-2, 15, 200)

# 
plt.figure(dpi=96)
plt.plot(tt, f_VC(tt),label='VC')
plt.plot(tt, f_VR(tt),label='VR')
plt.legend() 
plt.show()



#時定数を求める
t_1=C.subs(prm)*R.subs(prm)
print("時定数CR:{}秒".format(t_1))

以下、コードの実行結果である。

RL直列回路について

 

微分方程式

電流に関する微分方程式は以下の様になる。

この微分方程式は右辺が０でないので、特殊解を用いて工夫する必要がある。

この場合、解の１つである特殊解は以下の様になる。

次に、右辺を０にした微分方程式である

を解く方法を考える。

と変形すれば、

ゆえに、求める電流値は以下の様になる。

ここで、i(t=0)=0という境界条件（初期条件）を用いると、

 

プログラミングコード

これをPythonで表すと以下のコードになる。

import numpy as np
import matplotlib.pyplot as plt

# シンボル（変数記号）の定義
q,i = sympy.symbols('q i', cls=sympy.Function)
t,E,R,L,C1 = sympy.symbols('t E R L C1')

# 方程式 eq1 を立てる
eq1 = sympy.Eq( E, L * sympy.Derivative(i(t), t) + R*i(t) )
print("微分方程式の表示")
display(eq1)

eq2 = sympy.dsolve(eq1, i(t))
print("Qについて解く")
display(eq2)

#初期条件の代入
# eq2の左辺を0、右辺にt=0を代入
eq3 = sympy.Eq(0, eq2.rhs.subs(t,0),)  
C1a = sympy.solve(eq3,C1)[0]
print("初期条件の代入")
display(C1a)
#iの解を表示する
# 初期条件を与えて求めた積分定数C1aを代入
eq4 = eq2.subs(C1,C1a)
print("iの解の表示")
display(eq4) 

# 式を整理(展開)
eq4 = eq4.expand()
print("iの解の展開表示")
display(eq4) 


i = eq4.rhs
print("電流値の表示")
display(i)
# 抵抗両端の電圧 VR 
VR = R*i
print("抵抗Rの両端の電圧")
display(VR)

# コンデンサ両端の電圧 VL
VL = E - VR

print("コイルの電圧")
display(VL)
display(VL.collect(E)) # 電圧Eで括る 

# パラメータ
prm = { E:sympy.Piecewise((0,t<0),(5,t>=0)), R:5, L:40}

# パラメータを代入
i = i.subs(prm)
VR=VR.subs(prm)
VL=VL.subs(prm)

# sympy オブジェクトを numpyユニバーサル関数に変換
f_VR = sympy.lambdify(t, VR, 'numpy')
f_VL = sympy.lambdify(t, VL, 'numpy')

tt = np.linspace(-2, 15, 200)

# 
plt.figure(dpi=96)
plt.plot(tt, f_VL(tt),label='VL')
plt.plot(tt, f_VR(tt),label='VR')
plt.legend() 
plt.show()

#時定数を求める
t_1=L.subs(prm)/R.subs(prm)
print("時定数L/R:{}秒".format(t_1))

以下コードの実行結果である。

時定数を考慮したグラフの作成

RL回路

上記のコードの場合、時定数が変わるごとに手動で観測期間を変化させなければならない。なので、それを自動化したものを以下に述べる。

%reset -f
import sympy
import numpy as np
import matplotlib.pyplot as plt

# シンボル（変数記号）の定義
q,i = sympy.symbols('q i', cls=sympy.Function)
t,E,R,L,C1 = sympy.symbols('t E R L C1')

# 方程式 eq1 を立てる
eq1 = sympy.Eq( E, L * sympy.Derivative(i(t), t) + R*i(t) )
print("微分方程式の表示")
display(eq1)

eq2 = sympy.dsolve(eq1, i(t))
print("Qについて解く")
display(eq2)

#初期条件の代入
# eq2の左辺を0、右辺にt=0を代入
eq3 = sympy.Eq(0, eq2.rhs.subs(t,0),)  
C1a = sympy.solve(eq3,C1)[0]
print("初期条件の代入")
display(C1a)
#iの解を表示する
# 初期条件を与えて求めた積分定数C1aを代入
eq4 = eq2.subs(C1,C1a)
print("iの解の表示")
display(eq4) 

# 式を整理(展開)
eq4 = eq4.expand()
print("iの解の展開表示")
display(eq4) 


i = eq4.rhs
print("電流値の表示")
display(i)
# 抵抗両端の電圧 VR 
VR = R*i
print("抵抗Rの両端の電圧")
display(VR)

# コンデンサ両端の電圧 VL
VL = E - VR

print("コイルの電圧")
display(VL)
display(VL.collect(E)) # 電圧Eで括る 

#パラメーターの入力

r=int(input("抵抗値を入力してください"))
l=int(input("インダクタンスを入力してください"))

# パラメータ
prm = { E:sympy.Piecewise((0,t=0)), R:r, L:l*10e3}

# パラメータを代入
i = i.subs(prm)
VR=VR.subs(prm)
VL=VL.subs(prm)

# sympy オブジェクトを numpyユニバーサル関数に変換
f_VR = sympy.lambdify(t, VR, 'numpy')
f_VL = sympy.lambdify(t, VL, 'numpy')
t_1=int(L.subs(prm)/R.subs(prm))
tt = np.linspace(-0.5*t_1,2*t_1 , 200)

# 
plt.figure(dpi=96)
plt.plot(tt, f_VL(tt),label='VL')
plt.plot(tt, f_VR(tt),label='VR')
plt.legend() 
plt.show()

#時定数を求める
t_1=L.subs(prm)/R.subs(prm)
print("時定数L/R:{}秒".format(t_1))

RC回路

上記のコードの場合、時定数が変わるごとに手動で観測期間を変化させなければならない。なので、それを自動化したものを以下に述べる。

%reset -f
import sympy
import numpy as np
import matplotlib.pyplot as plt

# シンボル（変数記号）の定義
q,i = sympy.symbols('q i', cls=sympy.Function)
t,E,R,C,C1 = sympy.symbols('t E R C C1')

# 方程式 eq1 を立てる
eq1 = sympy.Eq( E, R * sympy.Derivative(q(t), t) + q(t)/C )
print("微分方程式の表示")
display(eq1)

eq2 = sympy.dsolve(eq1, q(t))
print("Qについて解く")
display(eq2)

#初期条件の代入
# eq2の左辺を0、右辺にt=0を代入
eq3 = sympy.Eq(0, eq2.rhs.subs(t,0),)  
C1a = sympy.solve(eq3,C1)[0]
print("初期条件の代入")
display(C1a)
#Qの解を表示する
# 初期条件を与えて求めた積分定数C1aを代入
eq4 = eq2.subs(C1,C1a)
print("Qの解の表示")
display(eq4) 

# 式を整理(展開)
eq4 = eq4.expand()
print("Qの解の展開表示")
display(eq4) 

# eq4 の右辺を t 微分
eq5 = sympy.Eq( i(t), sympy.diff(eq4.rhs,t))
print("電荷を微分して電流値を求める")
display( eq5 )

i = eq5.rhs

# 抵抗両端の電圧 VR 
VR = R*i
print("抵抗Rの両端の電圧")
display(VR)

# コンデンサ両端の電圧 VC
VC = E - VR

print("コンデンサの電圧")
display(VC)
display(VC.collect(E)) # 電圧Eで括る 
#パラメーターの入力

r=int(input("抵抗値を入力してください"))
c=int(input("静電容量を入力してください"))
# パラメータ:E=5[V] (t>0)でのステップ応答　
prm = { E:sympy.Piecewise((0,t=0)), R:r*10e3, C:c*10e-6 }

# パラメータを代入
i = i.subs(prm)
VR=VR.subs(prm)
VC=VC.subs(prm)

# sympy オブジェクトを numpyユニバーサル関数に変換
f_VR = sympy.lambdify(t, VR, 'numpy')
f_VC = sympy.lambdify(t, VC, 'numpy')

#時定数を求める
t_1=C.subs(prm)*R.subs(prm)
print("時定数CR:{}秒".format(t_1))


tt = np.linspace(int(-0.5*t_1), int(2*t_1), 200)



# 
plt.figure(dpi=96)
plt.plot(tt, f_VC(tt),label='VC')
plt.plot(tt, f_VR(tt),label='VR')
plt.legend() 
plt.show()

これらは、各時定数に合わせた範囲において自動的に最適な測定領域でプロットすることが出来る（-0.5T<t<2T:Tは時定数）

まとめ

RL回路やRC回路といった一次遅れ系は、一階級の微分方程式を用いて解析することが出来る。今回は、sympyというライブラリを用いて、一階級の微分方程式を解くことで、各パラメーターの過渡応答を評価した。

参考文献

PythonとRC回路の過渡現象