今回は伝達関数から周波数解析する方法を書く。
パソコンでプログラムからメソッドやライブラリを使うと簡単にできる。
しかし、手書きでできるところは手で計算し、プログラムも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であることに注意。
ボード線図からわかる情報については以下の記事で書いた。
コメント