伝達関数から周波数解析

古典制御

今回は伝達関数から周波数解析する方法を書く。
パソコンでプログラムからメソッドやライブラリを使うと簡単にできる。
しかし、手書きでできるところは手で計算し、プログラムもfor文を使って書いて計算させるとまた理解が深まる。

システム

以下のような線形時不変系システムを考える。

ここで入力から出力までの伝達関数は以下のように示せる。
\begin{align}
G(s)=\frac{X(s)}{F(s)}=\frac{1}{ms^2+cs+k}\\
G(s)=\frac{Kw_{n}^2}{s^2+2\zeta w_{n} s+w_{n}^2} \tag{1}
\end{align}
ここで$w_{n}$は固有角振動数を表し、$\zeta$は減衰比、$K$はゲインである。
それぞれ次のような関係である。

\begin{align}
w_{n}=\sqrt{\frac{k}{m}}, \quad \zeta=\frac{c}{2\sqrt{mk}},\quad K=\frac{1}{k}
\end{align}

周波数特性

(1)式の$s$に$jw$を代入すると周波数伝達関数になる。
\begin{gather}
G(jw)=\frac{X(jw)}{F(jw)}=\frac{Kw_{n}^2}{(jw)^2+2\zeta w_{n} jw+w_{n}^2}\\
=\frac{Kw_{n}^2}{w_{n}^2-w^2+2j\zeta w_{n}w}\tag{2}
\end{gather}
ここで$G(jw)=a+bj$の形にするために、
$(w_{n}^2-w^2)=A$と置きます。そして次のように整理していきます。
\begin{gather}
G(jw)=\frac{K w_{n}^2}{A+2\zeta w_{n} w j}\\
=\frac{K w_{n}^2 (A-2 \zeta w_{n}wj)}{(A+2 \zeta w_{n}wj)(A-2 \zeta w_{n}wj)}\\
=\frac{Kw_{n}^2A-2Kw_{n}^3 \zeta w j}{A^2+4 \zeta^2 w_{n}^2w^2}\\
=\frac{K w_{n}^2 A}{A^2+4 \zeta^2 w_{n}^2 w^2}+(-\frac{2Kw_{n}^3 w}{A^2+4 \zeta^2 w_{n}^2 w^2})j
\end{gather}
したがって
\begin{gather}
a=\frac{K w_{n}^2 A}{A^2+4 \zeta^2 w_{n}^2 w^2}\\
b=-\frac{2Kw_{n}^3 w}{A^2+4 \zeta^2 w_{n}^2 w^2}
\end{gather}

ゲイン特性

ゲイン特性は$\sqrt{a^2+b^2}$より求まる。
\begin{gather}
\sqrt{a^2+b^2}=\sqrt{(\frac{K w_{n}^2 A}{A^2+4 \zeta^2 w_{n}^2 w^2})^2+(-\frac{2Kw_{n}^3 w}{A^2+4 \zeta^2 w_{n}^2 w^2})^2}\\
=\sqrt{\frac{K^2 w_{n}^4A^2+4K^2w_{n}^6w^2 \zeta^2}{(A^2+4 \zeta w_{n}^2w^2)^2}}\\
=\sqrt{\frac{K^2w_{n}^2(A^2+4w_{n}^2w^2\zeta^2)}{(A^2+4 \zeta w_{n}^2w^2)^2}}\\
=\frac{Kw_{n}^2\sqrt{A^2+4 \zeta w_{n}^2w^2}}{\sqrt{A^2+4 \zeta w_{n}^2w^2} \sqrt{A^2+4 \zeta w_{n}^2w^2}}\\
(w_{n}^2-w^2)=Aより、
=\frac{Kw_{n}^2}{\sqrt{(w_{n}^2-w^2)^2+4 \zeta^2 w_{n}^2w^2}}
\end{gather}

ゲイン特性をdBで表示するためには次式を用いる
\begin{align}
G_{log}(w)=20\log_{10}|G(w)| \rm{[dB]}
\end{align}

もっと簡単に

(2)式はより丁寧に書くと
\begin{align}
\frac{Kw_{n}^2+0j}{(w_{n}^2-w^2)+2\zeta w_{n}wj}
\end{align}
と、分母分子、それぞれ虚部と実部を分けられる。
ここからゲインは、
\begin{align}
|G(w)| = \frac{\sqrt{(Kw_n)^2+0^2}}{\sqrt{(w_n^2-w^2)^2+(2\zeta w_n w)^2}}
\end{align}
と求めることもできる。

位相特性

位相特性は次式で求まる。
\begin{gather}
\large
\angle G(w) = \frac{X(w)}{F(w)}=\frac{b}{a}\\\\
=\tan^{-1}\frac{-\frac{2Kw_{n}^3w\zeta}{A^2+4\zeta^2 w_{n}^2w^2}}{\frac{Kw_{n}^2A}{A^2+4\zeta^2w_{n}^2w^2}}\\\\
=\tan^{-1}\frac{-2Kw_{n}^3w\zeta}{Kw_{n}^2A}\\\\
=\tan^{-1}\frac{-2\zeta w_{n}w}{w_{n}^2-w^2}
\end{gather}

もっと簡単に

(2)式はより丁寧に書くと
\begin{align}
\frac{Kw_{n}^2+0j}{(w_{n}^2-w^2)+2\zeta w_{n}wj}
\end{align}
と、分母分子、それぞれ虚部と実部を分けられる。

したがって位相特性は、分母と分子で位相を分けて考えると
\begin{align}
&\tan^{-1}\frac{0}{Kw_n^2}-\tan^{-1}\frac{2\zeta w_n w}{w_n^2-w^2}\\
&=-\tan^{-1}\frac{2\zeta w_n w}{w_n^2-w^2}
\end{align}

プログラム

では上の章で求めたゲイン特性と位相特性をプログラムで記述していきます。
ライブラリやメソッドを使うのではなく、for文で書くとより理解が深ります。
ここではmatlabファイルとpythonで記述します。
質量や粘性係数、ばね定数、外力は以下のように与えた。
\begin{align}
m = 5, ~ k = 20,~ c = 1,~ f=1
\end{align}

matlab(octave)


close all
clear
Tsim = 20;
sampleT = 0.01;
%--- -- パラメータ --------
m = 5;
k = 20;
c = 1;
omega = sqrt (k/m )
zeta = c /2/ sqrt (m *k)
K = 1/ k
omegad = omega * sqrt (1 -2* zeta ^2)

%--- -- 時間応答 ------
A = [0 1; - omega ^2 -2* zeta * omega ];
b = [0; K * omega ^2];
c = [1 0];
x0 = [0; 0];
Tn = round ( Tsim / sampleT );

f = ones (1 , Tn +1);
x_s = x0 ;
t_s = 0;

for i = 1: Tn + 1
t_s(i+1)=i*sampleT ;
% ルンゲクッタ法==
k1 = A* x_s (: , i) + b *f(i );
k2 = A*(x_s(:,i)+sampleT/2*k1)+b*f(i);
k3 = A*(x_s(:,i)+sampleT/2*k2)+b*f(i);
k4 = A*(x_s(:,i)+sampleT*k3)+b*f(i);
x_s(:,i+1)=x_s(:,i)+sampleT/6*(k1+2*k2+2*k3+k4);
%===========
t_o(i)=t_s(i);
x_o(i)=c*x_s(:,i);
end

%--- ---- 周波数応答 ----------
w = 0.1:0.01:10;
for i = 1 : length (w)
G(i) = K* omega ^2/ sqrt ( omega ^4+ w(i )^4+2*(2* zeta ^2 -1)* omega ^2* w (i )^2);
Gl ( i) = 20* log10 ( G(i ));
Ph ( i) = atan2 ( -2* zeta * omega *w (i) , omega ^2 - w(i )^2);
end


%--- -- Figure --------
fss = 20;
figure

subplot (111)
a1 = plot ( t_o , x_o ); grid ;
set ( gca ,'linewidth',1 , 'fontsize', fss )
set (a1 ,'linewidth',2);
xlabel ('time(s)','fontsize', fss );
ylabel   ('x[m]','fontsize', fss );
%ylim([0 0.06]);

figure
subplot (211)
a11 = semilogx (w , Gl ); grid ;
set ( gca ,'linewidth',1 , 'fontsize' , fss )
set ( a11 ,'linewidth',2);
ylabel ('Gain[db]', 'fontsize', fss );


subplot (212)
a12 = semilogx (w ,Ph ); grid ;
set ( gca ,'linewidth',1 , 'fontsize' , fss )
set ( a12 ,'linewidth',2);
ylabel ('Phase[rad]','fontsize', fss );
xlabel ( 'omega[rad/s]' ,'fontsize', fss );
時間応答

周波数応答

Python



import numpy as np
import matplotlib.pyplot as plt

# パラメータ
Tsim = 20
sampleT = 0.01
m = 5
k = 20
c = 1
omega = np.sqrt(k / m)
zeta = c / (2 * np.sqrt(m * k))
K = 1 / k
omegad = omega * np.sqrt(1 - 2 * zeta**2)

# 時間応答
A = np.array([[0, 1], [-omega**2, -2 * zeta * omega]])
b = np.array([[0], [K * omega**2]])
c = np.array([1, 0])
x0 = np.array([[0], [0]])
Tn = int(round(Tsim / sampleT))

f = np.ones(Tn + 1)
x_s = np.zeros((2, Tn + 1))
x_s[:, 0] = x0[:, 0]
t_s = np.zeros(Tn )
t_o = np.zeros(Tn  )
x_o = np.zeros(Tn  )

for i in range(Tn):
#     pythonは0から使えるため、i+1にする必要はない。
    t_s[i] = i* sampleT
    # ルンゲクッタ法
    k1 = A @ x_s[:, i] + b[:, 0] * f[i]
    k2 = A @ (x_s[:, i] + sampleT / 2 * k1) + b[:, 0] * f[i]
    k3 = A @ (x_s[:, i] + sampleT / 2 * k2) + b[:, 0] * f[i]
    k4 = A @ (x_s[:, i] + sampleT * k3) + b[:, 0] * f[i]
    x_s[:, i + 1] = x_s[:, i] + sampleT / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
    t_o[i] = t_s[i]
    x_o[i] = c @ x_s[:, i]

# 周波数応答
w = np.arange(0.1, 10, 0.01)
G = np.zeros(len(w))
Gl = np.zeros(len(w))
Ph = np.zeros(len(w))

for i in range(len(w)):
    G[i] = K * omega**2 / np.sqrt(omega**4 + w[i]**4 + 2 * (2 * zeta**2 - 1) * omega**2 * w[i]**2)
    Gl[i] = 20 * np.log10(G[i])
    Ph[i] = np.arctan2(-2 * zeta * omega * w[i], omega**2 - w[i]**2)

# Figure
fss = 20

FONTSIZE=12
fig,axes = plt.subplots(1,1,figsize=(10,10))
axes.plot(t_o,x_o)
axes.grid()
axes.set_xlabel('time, s',fontsize=FONTSIZE)
axes.set_ylabel('x, m',fontsize=FONTSIZE)

fig,axes = plt.subplots(2,1,figsize=(10,10))
axes = axes.flatten()
axes[0].semilogx(w,Gl)
axes[0].grid()
axes[0].set_ylabel('gain, dB',fontsize=FONTSIZE)

axes[1].semilogx(w,Ph)
axes[1].grid()
axes[1].set_xlabel('w, rad/s',fontsize=FONTSIZE)
axes[1].set_ylabel('phase rad',fontsize=FONTSIZE)
時間応答
周波数応答

周波数特性は対数で表示する。このときlog10であることに注意。
ボード線図からわかる情報については以下の記事で書いた。

bode線図からわかる情報
bode線図を初めて大学で勉強した際には周波数によってゲインや位相が変わることくらいしかわかりませんでした。しかし勉強していくにつれbode線図が持つリッチな情報量に気付き始めました。今回はその内容を書いていきます。 想定する読者 ・bod

コメント

タイトルとURLをコピーしました