RLC直列回路の伝達関数・状態空間モデルとPythonによるシミュレーション

# 概要

制御工学の勉強メモ。第2弾。

[前回](https://qiita.com/code0327/items/10fb56090a1e56046fa4)は、バネ・マス・ダンパ系の質点の運動方程式から、伝達関数／状態空間モデルを求めて、制御系のPythonライブラリ「[Python Control Systems Library](https://python-control.readthedocs.io/en/0.8.3/index.html)」を使ってシミュレーションをしました。

今回は、**RLC直列回路**を対象に、回路方程式から伝達関数／状態空間モデルを求めて、Pythonを使ってシミュレーションをしてみます。<br>![2020-04-28_09h00_39.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/152805/f10d0f88-eba1-d2ca-9ebc-75858504a7d3.png)

# 回路方程式

通常は、次のように **電流** $i(t)$ を使って回路方程式を立てます。

$$ v_{i}(t) = Ri(t) + L\frac{\mathrm{d}}{\mathrm{dt}}i(t) + \frac{1}{C}\int^{t}_{0} i(\tau)\,\mathrm{d} \tau$$

しかし、ここでは [前回](https://qiita.com/code0327/items/10fb56090a1e56046fa4)の力学系との対応をとるために、次の関係を使って **電荷** $q(t)$ を使って回路方程式を書いていきます。

$$ i(t) = \frac{\mathrm{d}}{\mathrm{dt}}q(t)$$

電荷 $q(t)$ で回路方程式を書くと次のようになります。

$$ v_{i}(t) = R\frac{\mathrm{d}}{\mathrm{dt}}q(t) + L\frac{\mathrm{d}^2}{\mathrm{dt}^2}q(t) + \frac{1}{C}q(t) $$

微分をドット表記に変えて、式を整理します。

$$ L \ddot{q}(t) + R \dot{q}(t) + \frac{1}{C}q(t) = v_{i}(t) $$

これは、[前回](https://qiita.com/code0327/items/10fb56090a1e56046fa4)で示したバネ・マス・ダンパ系における1自由度の質点の運動方程式（下記）と、同じ形になっていることが分かります。

$$ m\ddot{y}(t) + c\dot{y}(t) + ky(t) = f(t)$$

# 伝達関数でシステムをモデル化

力学のモデルでは、力 $f(t)$ をシステムに対する**入力**、変位 $y(t)$ をシステムからの**出力**として伝達関数を求めました。

仮に、それに対応させれば、入力電圧 $v\_i(t)$ をシステムに対する**入力**、電荷 $q(t)$ をシステムからの**出力**として伝達関数を求める形になるのですが、**電荷** $q(t)$ **を直接的に観測（計測）することは簡単ではありません**。また、実際のところ、関心があるのは各素子の電圧になります。よって、**電気回路の場合は、電圧を出力として伝達関数を求めます**。

ここでは、入力電圧 $v\_i(t)$ をシステムに対する**入力** $u(t)$、コンデンサの両端電圧 $v_{C}(t)$ を、システムの**出力** $y(t)$ とした伝達関数を求めていきます。

$$ 入力\,\,\,u(t) = v_{i}(t) $$

$$ 出力\,\,\,y(t) = v_{C}(t) = \frac{1}{C}q(t) $$

上記の出力を、$q(t)$ についての式にします。

$$ q(t) = Cy(t) $$

これを、先の回路方程式 $ L \ddot{q}(t) + R \dot{q}(t) + \frac{1}{C}q(t) = v_{i}(t)$ に代入すれば、次を得ます。

$$ LC \ddot{y}(t) + RC \dot{y}(t) + y(t) = u(t) $$

これをラプラス変換して、システムの**伝達関数** $G(s)=Y(s)/U(s)$ を求めると次のようになります（この際、$\dot{y}(0)=0$ と $y(0)=0$ を仮定します）。

$$ LCs^2Y(s) + RCsY(s) + Y(s) = U(s)$$

$$ \big(LCs^2+ RC s + 1\big) Y(s) = U(s)$$

$$ G(s) = \frac{Y(s)}{U(s)} = \frac{1}{LC s^2 + RC s + 1}$$ 

なお、力学モデルでは、力を入力、変位を出力として、質量を $m$、減衰係数を $c$、バネ定数を $k$ とした伝達関数は次のようになりました。力学系も回路系も非常に似た形になっていることが分かります。

$$ G(s) = \frac{Y(s)}{F(s)} = \frac{1}{m s^2 + c s + k}$$ 

# Pythonでモデル化・シミュレーション

Python の **Python Control Systems Library** を利用して、シミュレーションをします。伝達関数を設定してシステムをモデル化して、**ステップ応答**をシミュレーションします。

力学系では、質点に一瞬だけ力を加えたときの変位を観測するために**インパルス応答**のシミュレーションをしましたが、電気回路では、入力電圧 $v_{i}(t)=1\, (t\ge0)$ としたときのコンデンサ電圧 $v_{C}$ を観測するという設定のほうが自然なので**ステップ応答**をみていきます。

```py:準備（GoogleColab.環境）
!pip install --upgrade sympy
!pip install slycot
!pip install control
```

```py:ステップ応答のシミュレーション（伝達関数）
import numpy as np
import control.matlab as ctrl
import matplotlib.pyplot as plt

R = 2       # 抵抗値
L = 1e-3    # インダクタンス
C = 10e-6   # 静電容量

sys = ctrl.tf((1),(L*C,R*C,1))
print(sys)

t_range = (0,0.01)
y, t = ctrl.step(sys, T=np.linspace(*t_range, 400))

plt.figure(figsize=(7,4),dpi=120,facecolor='white')
plt.hlines(1,*t_range,colors='gray',ls=':')
plt.plot(t,y)
plt.xlim(*t_range)
plt.ylim(bottom=0)
plt.show()
```

```:実行結果
      1
------------
s^2 + s + 10
```

![f1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/152805/23110198-7a19-ac8d-863c-a3328745ddbc.png)

入力電圧を $v_{i}(t)=1\, (t\ge0)$ とすると、コンデンサ電圧 $v_{C}$ は、最初に振動して、最終的には $1\,\[\mathrm{V}\]$ に落ち着くことが分かります。

# 状態空間モデルでモデル化

システムを、次のような状態空間モデル $\mathcal{P}$ として記述します。

```math:
\mathcal{P}\,:\,\left\{
\begin{array}{l}
\dot{\boldsymbol{x}} = \mathbf{A}\boldsymbol{x} + \mathbf{B}\boldsymbol{u}& 状態方程式\\
\boldsymbol{y} = \mathbf{C}\boldsymbol{x} + \mathbf{D}\boldsymbol{u}& 出力方程式（観測方程式）
\end{array}
\right.
```

準備として、回路方程式 $ LC \ddot{y}(t) + RC \dot{y}(t) + y(t) = u(t) $ を、 $\ddot{y}(t)=...$ に変形しておきます。

$$ \ddot{y}(t) =  -\frac{1}{LC} y(t)- \frac{R}{L} \dot{y}(t) + \frac{1}{LC} u(t)$$

$x\_1(t)=y(t)$、$x\_2(t)=\dot{y}(t)$ として、**状態**　$\boldsymbol{x}=\[x\_1\,\,x\_2]^T$ を定めます。

これより、$\dot{x}\_1(t)=\dot{y}(t)=x\_2$ となります。また、$\dot{x}\_2(t)=-\frac{1}{LC} y(t)- \frac{R}{L} \dot{y}(t) + \frac{1}{LC} u(t)$ となります。

以上より、次の**状態方程式**を得ることができます。

```math:
\left[
    \begin{array}{c}
      \dot{x}_1 \\
      \dot{x}_2
    \end{array}
  \right]
=\left[
    \begin{array}{cc}
      0 & 1  \\
      -\frac{1}{LC} & -\frac{R}{L} 
    \end{array}
  \right]
  \left[
    \begin{array}{c}
      x_1 \\
      x_2
    \end{array}
  \right]
+ \left[
    \begin{array}{c}
      0 \\
      \frac{1}{LC}
    \end{array}
  \right] u
```

```math:
\dot{\boldsymbol{x}} 
=\left[
    \begin{array}{cc}
      0 & 1  \\
      -\frac{1}{LC} & -\frac{R}{L} 
    \end{array}
  \right]
\boldsymbol{x} + \left[
    \begin{array}{cc}
      0 \\
      \frac{1}{LC}
    \end{array}
  \right] u
```

$$\dot{\boldsymbol{x}} 
=\mathbf{A}\boldsymbol{x} + \mathbf{B} u$$


また、次の**出力方程式**（観測方程式）を得る（**コンデンサ電圧** $v_{C}$ に相当する状態 $x_1$ と、**電流**に相当する $Cx\_2 (= C\dot{y}(t) = C\dot{v}_{C}(t) = C\big(\frac{1}{C}\dot{q}(t)\big)= \dot{q}(t) = i(t)\,) $ を観測する形）。

```math:
y = \left[
    \begin{array}{cc}
      1 & 0\\
      0 & C\\
    \end{array}
  \right]   \left[
    \begin{array}{c}
      x_1 \\
      x_2
    \end{array}
  \right]
```

$$y
=\mathbf{C}\boldsymbol{x} + \mathbf{D} u$$

ただし、$\mathbf{D}=0$。

# Pythonでモデル化・シミュレーション

状態空間モデルで、**ステップ応答**をシミュレーションします。

```py:グラフで日本語利用するための準備（GoogleColab.環境）
!pip install japanize-matplotlib
```

```py:ステップ応答のシミュレーション（状態空間モデル）
import numpy as np
import control.matlab as ctrl
import matplotlib.pyplot as plt
import japanize_matplotlib

R = 2       # 抵抗値
L = 1e-3    # インダクタンス
C = 10e-6   # 静電容量

A = [[0,1],[-1/(L*C), -R/L]]
B = [[0],[1/(L*C)]]
C = [[1,0],[0,C]]
D = 0

sys = ctrl.ss(A,B,C,D)
#print(sys)

t_range = (0,0.01)
y, t = ctrl.step(sys, T=np.linspace(*t_range, 400))

fig, ax = plt.subplots(nrows=2, figsize=(7,5),dpi=120,sharex='col')
fig.patch.set_facecolor('white') 
fig.subplots_adjust(hspace=0.1)  
for i,s in enumerate(['コンデンサ電圧','電流']) :
  ax[i].hlines(0,*t_range,colors='gray',ls=':')
  ax[i].plot(t,y[:,i])
  ax[i].set_ylabel(s)
  ax[i].set_xlim(*t_range)
plt.show()
```

![f2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/152805/29f01399-8970-a198-8580-ec5bbb865919.png)