システムとモデリング

modelica, Julia, Design Structure Matrix, SysML, 他モデリング全般について。

トップ > Julia言語 > 【プログラミング言語Julia】DifferentialEquations.jlでRLC回路の微分代数方程式(DAE)を解く

【プログラミング言語Julia】DifferentialEquations.jlでRLC回路の微分代数方程式(DAE)を解く

Julia言語 微分方程式

今回はJuliaの強力な微分方程式数値計算パッケージDifferentialEquations.jlを使用してRLC回路を解いていきます。

パッケージの詳細は以下です。

diffeq.sciml.ai

また関連する記事は以下になります。

otepipi.hatenablog.com

otepipi.hatenablog.com

RLC回路

今回は以下サイトでmodelicaの例題として使用されているRLC回路を計算します f:id:Otepipi:20210603181236p:plain

mbe.modelica.university

こちらの方程式は次のようになります。 \begin{align} V = i_RR \end{align}

\begin{equation} C\frac{dV}{dt} = i_C \end{equation}

\begin{equation} L\frac{di_L}{dt} = V_b - V \end{equation}

\begin{equation} i_L = i_R + i_C \end{equation}

数式の整理

この式をDifferentialEquations.jlで読み込める形式に変換します。以下のように残差方程式に変形します。

\begin{align} 0 = -V + i_RR \end{align}

\begin{equation} 0 = -C\frac{dV}{dt} + i_C \end{equation}

\begin{equation} 0 = -L\frac{di_L}{dt} + V_b - V \end{equation}

\begin{equation} 0 = -i_L + i_R + i_C \end{equation}

さらに、このうち変数 V, i_R, i_C, i_Lをそれぞれu[1], u[2], u[3], u[4]に置き換え、それら変数の時間微分をdu[1], du[2], du[3], du[4]に置換します。さらにパラメーター R, C, L, V_bをそれぞれp[1],p[2],p[3],p[4]に置き換えます。 最後に、各数式行の左辺0をout[1],out[2],out[3],out[4]に変換すれば完成です。

\begin{align} out[1] = -u[1] + p[1]*u[2] \end{align}

\begin{equation} out[2] = -p[2] * du[1] + u[3] \end{equation}

\begin{equation} out[3] = -p[3]*du[4] + p[4] - u[1] \end{equation}

\begin{equation} out[4] = u[2] + u[3] - u[4] \end{equation}

では、Juliaでプログラミングしていきます。

Juliaプログラム

まず使用するパッケージを宣言します。

using DifferentialEquations, Plots

次に、先程作成した数式を関数として定義してやります。

function f(out,du,u,p,t)
  out[1] = -u[1] + p[1]*u[2]
  out[2] = -p[2] * du[1] + u[3]
  out[3] = -p[3]*du[4] + p[4] - u[1]
  out[4] = u[2] + u[3] - u[4]
end

f (generic function with 1 method)

初期条件、パラメーター、タイムスパンを定義します。

u₀ = [0.0, 0.0, 0.0, 0.0]
du₀ = [0.0, 0.0, p[4]/p[3], p[4]/p[3]]
p = [100, 1e-3, 1, 24]
tspan = (0.0,100000.0)

differential_varsはu[1]~u[4]のうち微分方程式で定義される変数をtrueにします。代数方程式で定義されるu[4]を除いてtrueにしてやります。

differential_vars = [true,true,true,false]
prob = DAEProblem(f,du₀,u₀,tspan,p,differential_vars=differential_vars)

今回はDAEソルバーとしてSundials.jlを使用します。

using Sundials
sol = solve(prob,IDA())

retcode: Success
Interpolation: 3rd order Hermite
t: 263-element Vector{Float64}:
      0.0
      2.9462782549439476e-8
      5.892556509887895e-8
      1.178511301977579e-7
      2.357022603955158e-7
      4.714045207910316e-7
      9.428090415820632e-7
      1.8856180831641265e-6
      3.771236166328253e-6
      7.542472332656506e-6
      1.5084944665313012e-5
      3.0169889330626023e-5
      6.0339778661252046e-5
      ⋮
     58.33849213536743
    114.02472507418413
    225.39719095181755
    448.14212270708435
    893.631986217618
   1784.6117132386853
   3566.5711672808197
   7130.490075365089
  14258.327891533627
  28514.003523870702
  57025.35478854485
 100000.0
u: 263-element Vector{Vector{Float64}}:
 [0.0, 0.0, 0.0, 0.0]
 [2.0833327195237355e-11, 2.0833327195237358e-13, 7.071065728526617e-7, 7.071067811859337e-7]
 [6.249997544759982e-11, 6.249997544759982e-13, 1.4142129373708854e-6, 1.4142135623706398e-6]
 [1.964284626964042e-10, 1.9642846269640422e-12, 2.82842516045069e-6, 2.8284271247353165e-6]
 [7.186140704521984e-10, 7.186140704521984e-12, 5.6568470632868595e-6, 5.656854249427564e-6]
 [2.7314763464374074e-9, 2.7314763464374075e-11, 1.1313681183707797e-5, 1.1313708498471261e-5]
 [1.0754021372318532e-8, 1.075402137231853e-10, 2.2627309453714978e-5, 2.26274169939287e-5]
 [4.277019566458438e-8, 4.2770195664584375e-10, 4.5254406261876154e-5, 4.5254833963832796e-5]
 [1.7079168266682338e-7, 1.707916826668234e-9, 9.05079598188744e-5, 9.050966773570107e-5]
 [6.827936630492731e-7, 6.827936630492732e-9, 0.00018101250599955734, 0.00018101933393618782]
 [2.7306758133720213e-6, 2.7306758133720214e-8, 0.0003620113488338514, 0.0003620386555919852]
 [1.0921551727595293e-5, 1.0921551727595293e-7, 0.0007239679974318806, 0.0007240772129491565]
 [4.368039055157781e-5, 4.368039055157781e-7, 0.0014477168362187045, 0.00144815364012422]
 ⋮
 [24.0, 0.24000000000000002, -3.239828388946411e-20, 0.24]
 [24.0, 0.24, 3.1396372559514076e-20, 0.24]
 [24.0, 0.24, -4.962917911371169e-22, 0.24]
 [24.0, 0.24, -4.972552283062986e-22, 0.24]
 [24.0, 0.24000000000000002, -4.976132462749414e-22, 0.24000000000000002]
 [24.0, 0.24, -4.975237407815362e-22, 0.24]
 [24.0, 0.24, -4.9752374077754535e-22, 0.24]
 [24.0, 0.24, -4.975237407775453e-22, 0.24]
 [24.0, 0.24, 9.043745455248383e-25, 0.24]
 [24.0, 0.24000000000000002, 9.040249044989204e-25, 0.24000000000000002]
 [24.0, 0.24, 9.041123139850977e-25, 0.24]
 [24.0, 0.24, 9.040856702059545e-25, 0.24]

無事方程式が解けたようなので、結果をプロットします。

plot(sol.t,sol[1,:],xlim=[0,2],label="V")

f:id:Otepipi:20210603183827p:plain

問題なくプロットできました。

今回はここまでにします。