2021-02-23

matlabでの伝達関数の定義法

matlabでの伝達関数の定義は単純で、例えば次のような伝達関数 $G(s)$
\begin{align}
G(s)=\dfrac{1}{s^{2}+2s+3}
\end{align}
であれば

Np = [0, 1] 
Dp = [1, 2, 3]
P = tf(Np, Dp)

とすればいい。
matlabのLisenceを所持していない場合、Python_Controlパッケージを用いることでほぼ同様に取り扱うことができる。
Python_Controlがインストールされている環境で

from control.matlab import *
Np = [0, 1]      
Dp = [1, 2, 3]  
P = tf(Np, Dp)
print('P(s)=', P)

とすれば同様の結果を得られる。

2021-01-26

浮体に働く復元力

数学物理

浮体に働く浮力ベクトルとその大きさは
\begin{align}
B = mg \hspace{5mm} W = - \rho g V
\end{align}
で表すことができる。浮体座標系から見れば
\begin{align}
\boldsymbol{f}_{g}^{n} = {}^{t}
\begin{pmatrix}
0 & 0 & B
\end{pmatrix} \hspace{5mm}
\boldsymbol{f}_{b}^{n} = - {}^{t}
\begin{pmatrix}
0 & 0 & W
\end{pmatrix}
\end{align}
重心と浮心の座標を
\begin{align}
\boldsymbol{r}_{g}= {}^{t} \begin{pmatrix}
x_{g} & y_{g} & z_{g}
\end{pmatrix}
\hspace{5mm} \boldsymbol{r}_{b} = {}^{t} \begin{pmatrix}
x_{b} & y_{b} & z_{b}
\end{pmatrix}
\end{align}
とすれば、グローバル座標系への変換は回転行列を用いれば次のようになる。

\begin{align}
\boldsymbol{f}_{g}^{b} = \boldsymbol{R}_{b}^{n} (\boldsymbol{\textit{η}})^{-1} \boldsymbol{f}_{g}^{n} \\
\boldsymbol{f}_{b}^{b} = \boldsymbol{R}_{b}^{n} (\boldsymbol{\textit{η}})^{-1} \boldsymbol{f}_{g}^{n}
\end{align}

これより，浮体に働く復元力ベクトルを $\boldsymbol{g} (\boldsymbol{\eta})$ は
\begin{align}
\boldsymbol{g}(\boldsymbol{\textit{η}}) &= -
\begin{pmatrix}
\boldsymbol{f}_{g}^{b} + \boldsymbol{f}_{b}^{b} \\
\boldsymbol{r}_{g}^{b} \boldsymbol{S} (\boldsymbol{f}_{g}^{b}) + \boldsymbol{r}_{g}^{b} \boldsymbol{S} (\boldsymbol{f}_{g}^{b})
\end{pmatrix} \\
&= -
\begin{pmatrix}
\boldsymbol{R}_{b}^{n} (\boldsymbol{\textit{η}})^{-1} \left ( \boldsymbol{f}_{g}^{b} + \boldsymbol{f}_{b}^{b} \right ) \\
\boldsymbol{r}_{g}^{b} \boldsymbol{S} \left( \boldsymbol{R}_{b}^{n} (\boldsymbol{\textit{η}})^{-1} \boldsymbol{f}_{n}^{b} \right ) + \boldsymbol{r}_{b}^{b} \boldsymbol{S} \left( \boldsymbol{R}_{b}^{n} (\boldsymbol{\textit{η}})^{-1} \boldsymbol{f}_{n}^{b} \right )
\end{pmatrix}
\end{align}
したがって次を得る。
\begin{align}
\boldsymbol{g}(\boldsymbol{\textit{η}}) =
\begin{pmatrix}
(m-\rho V) g \ \sin{\theta} \\
-(m-\rho V) g\ \sin{\phi} \cos{\theta} \\
-(m-\rho V) g\ \cos{\phi} \cos{\theta} \\
-(my_{g}-\rho V y_{b} ) g\ \cos{\phi} \cos{\theta} + (mz_{g}-\rho V z_{b} ) g\ \sin{\phi} \cos{\theta} \\
(mz_{g}-\rho V z_{b} ) g\ \sin{\theta} + (mx_{g}-\rho V x_{b} ) g\ \cos \phi \cos{\theta} \\
-(mx_{g}-\rho V x_{b} ) g\ \sin{\phi} \cos{\theta} - (my_{g}-\rho V y_{b} ) g\ \sin{\theta}
\end{pmatrix}
\end{align}
ここで、浮体の浮力が中性浮力であれば $W=B$ が成り立つので
\begin{align}
\begin{split}
\boldsymbol{g}(\boldsymbol{\textit{η}}) = \begin{pmatrix}
0\\
0\\
0\\
-(y_{g} - y_{b} ) mg \ \cos{\phi} \cos{\theta} + (z_{g} - z_{b} ) mg \ \sin{\phi} \cos{\theta} \\
(z_{g} - z_{b} ) mg\ \sin{\theta} + (x_{g} - x_{b} ) mg\ \cos \phi \cos{\theta} \\
-(x_{g} - x_{b} ) mg\ \sin{\phi} \cos{\theta} - (y_{g} - y_{b} ) mg\ \sin{\theta}
\end{pmatrix}
\end{split}
\end{align}

2021-01-19

回転行列から回転角を求める

数学

回転行列 $\boldsymbol{R}$ について

\begin{align}
\begin{pmatrix}
\cos \psi \cos \theta& \cos \psi \sin \phi \sin \theta - \cos \phi \sin \psi & \sin \phi \sin \psi + \cos \phi \cos \psi \sin \theta \\
\cos \theta \sin \psi& \cos \phi \cos \psi + \sin \phi \sin \psi \sin \theta & \cos \phi \sin \psi \sin \theta - \cos \psi \sin \phi \\
-\sin \theta & \cos \theta \sin \phi & \cos \phi \cos \theta
\end{pmatrix}=
\begin{pmatrix}
R_{11} & R_{12} & R_{13} \\
R_{21} & R_{22} & R_{23} \\
R_{31} & R_{32} & R_{33}
\end{pmatrix}
\end{align}
と定義する。これより回転行列から回転角は
\begin{align}
\phi &= \tan^{-1} \dfrac{R_{32} }{ R_{33}} \\
\theta &= - \sin^{-1} R_{31} =- \tan^{-1} \dfrac{R_{31} }{\sqrt{1- R_{31}^{2} }} \hspace{10mm} (\theta \neq \pm \frac{\pi}{2} )\\
\psi &= \tan^{-1} \dfrac{R_{21} }{ R_{11}}
\end{align}
と求まる。

2021-01-15

クォータニオンを定義する

数学

クォータニオンは１つの実部 $\textit{η}$ と３つの要素を持つ虚部 $\boldsymbol{\varepsilon}$ からなる。 $\boldsymbol{\varepsilon}$ は
\begin{align}
\boldsymbol{\varepsilon}= {}^{t}
\begin{pmatrix}
\varepsilon_{1} & \varepsilon_{2} & \varepsilon_{3}
\end{pmatrix}
\end{align}
である。ここで、 $\boldsymbol{q}$ は
\begin{align}
\boldsymbol{q} {}^{t} \boldsymbol{q} = 1
\end{align}
を満たしている。これを考慮して単位クォータニオンを定義すれば
\begin{align}
Q=\{ \boldsymbol{q} | \boldsymbol{q} {}^{t} \boldsymbol{q} = 1 , \boldsymbol{q} = {}^{t} ( \textit{η} \ {}^{t} \boldsymbol{\varepsilon} ) , \textit{η} \in \mathbb{R} \land \ \boldsymbol{\varepsilon} \in \mathbb{R}^3 \}
\end{align}
流儀により多少異なる定義になる場合があるが、その場合であってもクォータニオンの持つ性質に影響はない。
ここで $\boldsymbol{\varepsilon},\textit{η}$ はそれぞれ

\begin{align}
\boldsymbol{\varepsilon} = \boldsymbol{\lambda} \sin \frac{\delta}{2} \hspace{10mm} \textit{η}= \cos \frac{\delta}{2}
\end{align}

であるから $\boldsymbol{q}$ は

\begin{align}
\boldsymbol{q}= \textit{η}+ \varepsilon_{1} \boldsymbol{i} + \varepsilon_{2} \boldsymbol{j} + \varepsilon_{3} \boldsymbol{k} =\cos \frac{\delta}{2} + \sin \frac{\delta}{2} \boldsymbol{i} + \sin \frac{\delta}{2} \boldsymbol{j} + \sin \frac{\delta}{2} \boldsymbol{k}
\end{align}
となる。

次回以降、クォータニオンと回転行列の関係などを調べていく

2021-01-04

ボード線図を作図する

制御

前回、二次遅れ系の伝達関数 $G(s)$
\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}
のボード線図は
f:id:sironekoblog:20210104171542j:plain
となる。

同様に二次遅れ系
\begin{align}
G=\dfrac{1}{s^2+5s+6}
\end{align}
のボード線図は
f:id:sironekoblog:20210104171539j:plain
となる。これらは片対数グラフで描画刺されている。

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

2020-12-26

因数分解と零点

数学

中学生のころに学ぶ因数分解とは
\begin{align}
(x+a)(x+b)=x^2+(a+b)x+ab
\end{align}
のように式を積の形で表す操作である。上の例では $x+a,x+b$ が因数に当たる。このような式の段階では分かりづらいが、方程式が出現してからはその有用さは明らかである。
次のような一次、二次方程式を解くことを考える。
\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}
となり
f:id:sironekoblog:20201226032403j:plain
のようなグラフとなる。非常にわかりにくいが、 $x=-2,x=-3$ が $x$ 軸との共有点になっている。 $x$ 軸との共有点こそが零点である。

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

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

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

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

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

2020-12-23

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

制御

二次遅れ系の伝達関数 $G(s)$
\begin{align}
G(s)=\dfrac{K}{Ts+1}
\end{align}
を例に、ゲイン線図と位相線図からなるボード線図を作図する。
はじめに $s=j \omega$ としてフーリエ変換すると
\begin{align}
G(j \omega)&=\dfrac{K}{1 + j \omega T }
\end{align}
また、このシステムのゲイン $|G(j \omega)|$ は
\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}
を求めればよい。

2020-12-22

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

制御

ばねマスダンパ系の状態方程式を導出する。
まず、ばねマスダンパ系の運動方程式は
\begin{align}
M \ddot{x}(t) + C \dot{x}(t) + k x(t) = f(t)
\end{align}
ここで $M$ は質量、 $C$ は減衰係数、 $k$ はばね定数、 $x(t)$ は変位である。
はじめに
\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}
を得る。この式を状態方程式という。
観測できる量を $x$ とすれば
\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}
を得る。この式を出力方程式という。