ボード線図を作図する

前回、二次遅れ系の伝達関数
\begin{align}
G(s)=\dfrac{K}{Ts+1}
\end{align}
を例に、ゲイン線図と位相線図からなるボード線図を作図する過程を示した。
これによれば、ゲイン線図は
\begin{align}
\left | G(j\omega) \right | = \dfrac{K\sqrt{1+T^{2}\omega^{2}}}{1+T^{2}\omega^{2}}
\end{align}
を、位相線図は
\begin{align}
\angle G (j \omega) =\tan^{-1} \dfrac{\Im G(j \omega) }{\Re G(j \omega) }
\end{align}
を求める必要がある。
定義に従い数値計算ソフトを用いて描画すれば一次遅れ系
\begin{align}
G=\dfrac{5}{s+3}
\end{align}
のボード線図は

となる。

同様に二次遅れ系
\begin{align}
G=\dfrac{1}{s^2+5s+6}
\end{align}
のボード線図は

となる。これらは片対数グラフで描画刺されている。

次回はボード線図からゲイン余裕、位相余裕を求める。

因数分解と零点

中学生のころに学ぶ因数分解とは
\begin{align}
(x+a)(x+b)=x^2+(a+b)x+ab
\end{align}
のように式を積の形で表す操作である。上の例ではが因数に当たる。このような式の段階では分かりづらいが、方程式が出現してからはその有用さは明らかである。
次のような一次、二次方程式を解くことを考える。
\begin{align}
ax+b = 0 (a \neq 0) \hspace{10mm} cx^2+dx+e =0 (c \neq 0)
\end{align}
これらの方程式の解はそれぞれ
\begin{align}
x=\dfrac{b}{a} (a \neq 0) \hspace{10mm} x=\dfrac{-b \pm \sqrt{b^2 -4ac }}{2a}
\end{align}
で与えられる。このように解が求まるということは
\begin{align}
&x+\dfrac{b}{a} =0 \hspace{10mm} \left ( x+\dfrac{-b \pm \sqrt{b^2 -4ac}}{2a} \right ) \left ( x-\dfrac{-b \pm \sqrt{b^2 -4ac }}{2a} \right) =0
\end{align}
と因数分解できるということである。
例えば
\begin{align}
(x+2)(x+3)=0
\end{align}
について、展開して整理すれば
\begin{align}
x^2+5x+6=0
\end{align}
となり

のようなグラフとなる。非常にわかりにくいが、が軸との共有点になっている。軸との共有点こそが零点である。

零点の重要性は上の結果から次のように考えることができる。

• 1 関数を零点を用いて分解
• 2 零点から関数を作る

この２つの問題は非常に重要である。

ここでは２について、次のような零点
\begin{align}
x = 0, \pm \pi , \pm 2 \pi , \pm 3 \pi \cdots
\end{align}
を持つ関数を考える。先ほどの例のように考えればこのような零点を持つ関数は
\begin{align}
f(x) &= (x-\pi )(x+\pi)(x-2\pi)(x+2\pi)(x-3\pi)(x+3\pi) \cdots \\
&=\prod_{a=零点} (x-a)
\end{align}
と表示できるはずである。今回の関数は明らかにであるから
\begin{align}
\sin x =\prod_{a=零点} (x-a)
\end{align}
の表示を得る。

数学で１つの式に別の表示を与えることは重要である。実際この表示はオイラーがバーゼル問題を解くために使用している。
バーゼル問題については今回の記事の趣旨と異なるので次回以降投稿する。

ボード線図作図のためのゲイン線図と位相線図

二次遅れ系の伝達関数
\begin{align}
G(s)=\dfrac{K}{Ts+1}
\end{align}
を例に、ゲイン線図と位相線図からなるボード線図を作図する。
はじめにとしてフーリエ変換すると
\begin{align}
G(j \omega)&=\dfrac{K}{1 + j \omega T }
\end{align}
また、このシステムのゲインは
\begin{align}
\left | G(j\omega) \right |=\sqrt{\Re (j\omega)^{2}+\Im (j\omega)^{2}}
\end{align}
これらより、ゲインは
\begin{align}
\left | G(j\omega) \right | = \dfrac{K\sqrt{1+T^{2}\omega^{2}}}{1+T^{2}\omega^{2}}
\end{align}
最後に、ゲイン線図は対数グラフで表現するので
\begin{align}
20\log_{10}|G(j\omega)| = 20\log_{10} \left (\dfrac{K\sqrt{1+T^{2}\omega^{2}}}{1+T^{2}\omega^{2}} \right )
\end{align}
となる。
位相線図も同様にフーリエ変換し
\begin{align}
G(j \omega)&=\dfrac{K}{1 + j \omega T }
\end{align}
について、その角度
\begin{align}
\angle G (j \omega) =\tan^{-1} \dfrac{\Im G(j \omega) }{\Re G(j \omega) }
\end{align}
を求めればよい。

ばねマスダンパ系の状態方程式

ばねマスダンパ系の状態方程式を導出する。
まず、ばねマスダンパ系の運動方程式は
\begin{align}
M \ddot{x}(t) + C \dot{x}(t) + k x(t) = f(t)
\end{align}
ここでは質量、は減衰係数、はばね定数、は変位である。
はじめに
\begin{align}
x_{1}(t) = x(t) \hspace{5mm} x_{2}(t) = \dot{x}_{1} (t)
\end{align}
とおくと
\begin{align}
\begin{cases}
\dot{x}_{1}(t) = x_{2} (t)\\
\dot{x}_{2}(t) = -C x_{2} -k x_{1} = f
\end{cases}
\end{align}
行列では
\begin{align}
\begin{bmatrix}
\dot{x}_{1}(t)\\
\dot{x}_{2}(t)
\end{bmatrix}=
\begin{bmatrix}
0 & 1\\
-C & -k
\end{bmatrix}
\begin{bmatrix}
x_{1}(t)\\
x_{2}(t)
\end{bmatrix}=
\begin{bmatrix}
0\\ 1
\end{bmatrix}f
\end{align}
と表すことができる。ここで行列を
\begin{align}
\boldsymbol{x}=
\begin{bmatrix}
x_{1}(t)\\
x_{2}(t)
\end{bmatrix} \hspace{5mm}
\boldsymbol{A}=
\begin{bmatrix}
0 & 1\\
-C & -k
\end{bmatrix}\hspace{5mm}
\boldsymbol{B}=
\begin{bmatrix}
0\\ 1
\end{bmatrix}
\end{align}
とおけば
\begin{align}
\dot{\boldsymbol{x}}=\boldsymbol{A}\boldsymbol{x} + \boldsymbol{B} f
\end{align}
を得る。この式を状態方程式という。
観測できる量をとすれば
\begin{align}
y=
\begin{bmatrix}
1 &0
\end{bmatrix}
\begin{bmatrix}
x_{1}\\
x_{2}
\end{bmatrix}
\end{align}
同様に行列を置き換えれば
\begin{align}
\boldsymbol{y}=\boldsymbol{C}\boldsymbol{x}
\end{align}
を得る。この式を出力方程式という。

制御偏差と定常偏差

フィードバック制御系において目標値が変化したり外乱が加わると制御量が変化する。制御量は出力とも呼ばれ、制御量と目標値の差は（制御）偏差と呼ばれる。
適切に設計されたフィードバック制御系であれば徐々に偏差は減少するが、定常状態となった後にも残っている場合がある。このような偏差を定常偏差という。
フィードバック制御系の偏差は
\begin{align}
E(s)=\dfrac{R(s)}{1 + G(s)}
\end{align}
となる。定常偏差は
\begin{align}
e(t \to \infty ) = \lim_{t \to \infty} e(t) = \lim_{t \to 0} s E(s)
\end{align}
で与えられる。

pandasを使ったcsvファイルの読み込み

pandasを使ったcsvファイルの操作が思いのほか使いやすかったので残しておく。
pandasがインストールされている環境で

import pandas as pd
df = pd.read_csv('FILEPASS/FILENAME.csv',encoding = 'cp932',usecols=[2])
print(df)

とすればcsvファイルの2列目を読み込むことができる。
encoding = 'cp932'はエンコード方式の指定、usecols=[2]は読み込む列番号の指定をする。

print(df)は確認
コンソール画面に結果が表示される（はず）

二次遅れ系の過渡応答の極大極小値

前回二次遅れ系の伝達関数
\begin{align}
G(s)=\dfrac{\omega_{n}^2}{s^2 + 2 \zeta \omega_{n} s+ \omega_{n}^2}
\end{align}
についての行過ぎ量を求める式
\begin{align}
f(t) = 1 - (-1)^n e^{ -\zeta \dfrac{n \pi }{\sqrt{1-\zeta^2}} } \hspace{5mm} (n=0,1,2 \cdots)
\end{align}
を求めた。
まず二次遅れ系のインディシャル応答は
\begin{align}
f(t)=1-\dfrac{e^{- \zeta \omega_{n} t} }{\sqrt{1-\zeta^2}} \sin \left (\sqrt{1-\zeta^2} \omega_{n} t + \tan^{-1} \dfrac{\sqrt{1-\zeta^2}}{\zeta} \right )
\end{align}
ここで
\begin{align}
\omega = \sqrt{1-\zeta^2} \omega_{n} \hspace{5mm} \phi=\tan^{-1} \dfrac{\sqrt{1-\zeta^2}}{\zeta}
\end{align}
とすれば
\begin{align}
f(t)=1-\dfrac{e^{- \zeta \omega_{n} t} }{\sqrt{1-\zeta^2}} \sin \left ( \omega t + \phi \right )
\end{align}
と簡潔に表すことができた。また前回の結果から応答の極大極小値の座標は
\begin{align}
\left ( \dfrac{n \pi }{\sqrt{1-\zeta^2} \omega_{n}},1 - (-1)^n e^{ -\zeta \dfrac{n \pi }{\sqrt{1-\zeta^2}} } \right )
\end{align}
となることからグラフを書けば

となる。グラフから明らかのようにのとき極大、のとき極小値をとる。

二次遅れ系のインディシャル応答と行過ぎ量

次のような二次遅れ系の伝達関数
\begin{align}
G(s)=\dfrac{\omega_{n}^2}{s^2 + 2 \zeta \omega_{n} s+ \omega_{n}^2}
\end{align}
について、不足振動となるの行き過ぎ量を考える。
この場合の応答は
\begin{align}
y(t)=1-\dfrac{e^{- \zeta \omega_{n} t} }{\sqrt{1-\zeta^2}} \sin \left (\sqrt{1-\zeta^2} \omega_{n} t + \tan^{-1} \dfrac{\sqrt{1-\zeta^2}}{\zeta} \right )
\end{align}
ここで
\begin{align}
\omega = \sqrt{1-\zeta^2} \omega_{n} \hspace{5mm} \phi=\tan^{-1} \dfrac{\sqrt{1-\zeta^2}}{\zeta}
\end{align}
とすれば
\begin{align}
y(t)=1-\dfrac{e^{- \zeta \omega_{n} t} }{\sqrt{1-\zeta^2}} \sin \left ( \omega t + \phi \right )
\end{align}
と簡潔に表すことができる。
ここで行き過ぎ量は
\begin{align}
\frac{dy(t)}{dt}=0
\end{align}
で求めることができるから
\begin{align}
\frac{dy(t)}{dt}&=\dfrac{d}{dt} \left \{1-\dfrac{e^{- \zeta \omega_{n} t} }{\sqrt{1-\zeta^2}} \sin \left ( \omega t + \phi \right ) \right \} \\
&=- \left (- \zeta \omega_{n} \right ) \dfrac{e^{- \zeta \omega_{n} t} }{\sqrt{1-\zeta^2}} \sin \left ( \omega t + \phi \right ) -\dfrac{e^{- \zeta \omega_{n} t} }{\sqrt{1-\zeta^2}} \omega \cos \left ( \omega t + \phi \right )
\end{align}
ここで条件より
\begin{align}
\zeta \omega_{n} \dfrac{e^{- \zeta \omega_{n} t} }{\sqrt{1-\zeta^2}} \sin \left ( \omega t + \phi \right ) &=\dfrac{e^{- \zeta \omega_{n} t} }{\sqrt{1-\zeta^2}} \omega \cos \left ( \omega t + \phi \right ) \\
\zeta \omega_{n} \sin \left ( \omega t + \phi \right ) &= \omega \cos \left ( \omega t + \phi \right ) \\
\tan \left ( \omega t + \phi \right ) &= \dfrac{\omega}{\zeta \omega_{n} }
\end{align}
これより
\begin{align}
\tan \left ( \omega t + \tan^{-1} \dfrac{\sqrt{1-\zeta^2}}{\zeta} \right ) &= \dfrac{\sqrt{1-\zeta^{2}} }{\zeta}
\end{align}
ここで
\begin{align}
\omega t = n \pi \hspace{5mm} (n=0,1,2 \cdots)
\end{align}
であれば上式は成り立つ。定義より
\begin{align}
\sqrt{1 -\zeta^2}\omega_{n} t= n \pi \hspace{5mm} t = \dfrac{n \pi }{\sqrt{1-\zeta^2} \omega_{n}}
\end{align}
となるから
\begin{align}
f(t) = 1 - \dfrac{e^{-\zeta \dfrac{n \pi }{\sqrt{1-\zeta^2} }}} {\sqrt{1-\zeta^2}} \sin \left ( n \pi + \tan^{-1} \dfrac{\sqrt{1-\zeta^2}}{\zeta} \right )
\end{align}
加法定理
\begin{align}
\sin(a \pm b) = \sin a \cos b \pm \cos a \sin b
\end{align}
より
\begin{align}
f(t) &= 1 - \dfrac{e^{-\zeta \dfrac{n \pi }{\sqrt{1-\zeta^2} }}} {\sqrt{1-\zeta^2}} \left \{ \sin n \pi \cos \left( \tan^{-1} \dfrac{\sqrt{1-\zeta^2}}{\zeta} \right ) + \cos n \pi \sin \left( \tan^{-1} \dfrac{\sqrt{1-\zeta^2}}{\zeta} \right ) \right\}\\
&= 1 - (-1)^n e^{ -\zeta \dfrac{n \pi }{\sqrt{1-\zeta^2}} }
\end{align}
より行き過ぎ量を求める式を得る。