概要¶
ターゲット
- 制御工学を一通りかじった程度の人間が分かるような理論要点集がほしい
- 一次遅れ・二次遅れについて上記理論を適用した要点集がほしい
- 理論は周波数領域の話ばかりなので時間領域でどうなるか具体的に知りたい
- 周波数領域での制御設計をプログラムベースに落とし込んだらどうなるかまとめたい
参考書 (非常にいい本だった)
- 現場で役立つ制御工学の基本 (涌井信二 他)
from numpy import *
from matplotlib.pyplot import *
from control.matlab import *
import itertools as it
import pandas as pd
import matplotlib as mpl
zoom=1.3
mpl.rcParams['figure.figsize'] = [6.4 *zoom , 4.8*zoom]
mpl.rcParams['figure.dpi'] = 80*zoom
mpl.rcParams['lines.linewidth'] = 2.0
ion()
基本的なこと¶
機械と電気の対応¶
| 並進 | 回転 | 電気 |
|---|---|---|
| 力 $f$ | トルク $\tau$ | 電圧 $v$ |
| 位置 $x$ | 角度 $\theta$ | 電荷 $q$ |
| バネ係数 $f=kx$ | バネ係数 $\tau=k\theta$ | 静電容量 $v=\frac{1}{C}q$ |
| 速度 $v=\frac{dx}{dt}$ | 角速度 $\omega=\frac{d\theta}{dt}$ | 電流 $i=\frac{dq}{dt}$ |
| ダンパ係数 $f=dv$ | ダンパ係数 $f=d\omega$ | 抵抗値 $v=Ri$ |
| 加速度 $a=\frac{dv}{dt}$ | 角加速度 $\alpha=\frac{d\omega}{dt}$ | ??? $ \frac{di}{dt}$ |
| 質量 $m$ | 慣性モーメント $J$ | インダクタンス $L$ |
周波数伝達関数¶
伝達関数$G(s)$に正弦波を入力したときの出力は、定常状態にて振幅が入力の$\alpha$倍、位相が入力から$\beta$[rad]遅れた正弦波になる $\alpha,\beta$は周波数伝達関数$G(j\omega)$を用いて以下で表す事ができる(付録で証明)
$$\alpha=|G(j \omega)| = \sqrt{(Re(G(j\omega))^2 + Im(G(j\omega))^2}$$$$\beta=\angle G(j\omega)=tan^{-1}\left(\frac{Im(G(j\omega))}{Re(G(j\omega))}\right)$$ボード線図はゲイン$=20log(|G(j \omega)|)$と、位相遅れ$\angle G(j\omega)$をそれぞれ横軸=周波数でプロットしたもの。 実験的には、正弦波の周波数を振りながら入力し、定常出力の振幅と位相遅れをプロットしていけば良い
ゲインの考え方の例は以下
- ゲイン20db = 入力の10倍
- ゲイン-20db = 入力の0.1倍
- ゲイン-3db = -20db+17db = 0.1倍*7.07倍 = 0.7倍
- ゲイン35db = 20db+15db = 10倍*5.62倍 = 56.2倍
- ゲイン2dbと5dbがどれだけ違う？ = (5db-2db) = 3db = 5dbは2dbの1.41倍
- ゲインをgとおくと倍率=$10^\frac{g}{20}=1.122^g$倍
- ゲインg1,g2の倍率比=$10^\frac{g1}{20} / 10^\frac{g2}{20}=1.122^{(g1-g2)}$倍
ゲインと$\alpha$倍の関係は以下の対応表参照。
x=np.linspace(-20,20,41)
y=pow(10,x/20)
figure()
plot(x,y,'.-')
grid()
xlabel("gain")
ylabel("ratio")
p = pd.Series(data=y,index=x)
print(p)
極配置の見方¶
- 根と実軸のなす角$\theta$が大きいと，減衰係数$\zeta$は小さい($cos\theta = \zeta$)
- 原点に近い極は周波数が低い
- 右半平面の極は発散する
安定性判別¶
─┬─C(s)─G(s)─┬─>
└───H(s)────┘- CL(s) : closed-loop 伝達関数
- F(s) : Forward 方向伝達関数
- OL(s) : open-loop 伝達関数 (一巡伝達関数)
- C(s) : コントローラ
- G(s) : 制御対象システム
- H(s) : 検出器
安定条件
- closed-loop伝達関数のすべての極がs平面の左半平面にあること
- フィードバック含めた系全体の極配置=応答を考えている
- ラウス・フルビッツの方法で解析的に確認できる
- 一巡伝達関数の位相余裕・ゲイン余裕が確保されていること
- 入力が一周して戻ってきたときに正帰還になっていないかを考えている
- ボード線図・ナイキスト線図で実験的に確認できる
- だいたいゲイン余裕10-20db，位相余裕40-60degくらいが良いらしい
外部安定
- LTIの安定と言われるのは大体BIBO(Boundary Input Boundary Output)安定
- ざっくりインパルス入力の応答が収束すること
- インパルス入力の応答の積分値が収束=有界=Boundaryって言葉が出てくる
- 系全体の入出力関係がBIBO安定 = 外部安定
- 系のサブシステム全ての入出力関係が安定 = 内部安定
- 極零相殺が起きてるようなシステムは外部安定に見えても，内部安定でない場合があるから注意
一次遅れ¶
一般系 $$G(s) = \frac{K}{Ts+1}$$ $$T\frac{dy}{dt}+y=Kx$$
インパルス応答解 $$\mathcal{L}^{-1}\left[\frac{K}{Ts+1}\right] = \frac{K}{T}e^{-\frac{t}{T}}$$ ステップ応答解 $$\mathcal{L}^{-1}\left[\frac{K}{Ts+1}\frac{1}{s}\right] = K \left( 1-e^{-\frac{t}{T}} \right)$$
周波数特性 $$|G(j\omega)| = K \frac{1}{\sqrt{1+(\omega T)^2}}$$ $$\angle G(j\omega) = tan^{-1}(T\omega)$$
- 折点周波数$f_0 =\frac{1}{2\pi T}$ Hz
- 低域ゲイン$g = 20 log_{10}(K)$ db
- 折点周波数にてゲイン=-3db,位相遅れ45deg
- ゲイン-3dbになる周波数を「バンド幅」という人もいる．一次遅れだと折点周波数に等しい
代表的な例¶
LR回路 (電流が電圧に対して一次遅れ) $$ \begin{eqnarray} L\frac{di}{dt} + Ri &=&v \\ \iff \frac{L}{R}\frac{di}{dt} + i &=&\frac{1}{R}v \end{eqnarray} $$ より， $$T=\frac{L}{R}$$ $$K=\frac{1}{R}$$
RC回路 (電荷∝充電電圧が電圧に対して一次遅れ) $$Ri+\frac{1}{C} \int qdt = v$$ ここで，$i=\frac{dq}{dt}$より， $$ \begin{eqnarray} R\frac{dq}{dt} + \frac{1}{C}q &=&v \\ \iff RC\frac{dq}{dt} + q &=& Cv \end{eqnarray} $$ より， $$T=RC$$ $$K=C$$
時間応答¶
Tが小さいと応答性が良い Kは応答性には影響せず，Kが大きいと出力スケールが大きくなる
def sys_1st(K,T):
return tf([K],[T,1])
def sys_2nd(K,wn,zeta):
return tf([K*wn*wn],[1,2*zeta*wn,wn*wn])
def get_responces(sysparam,sysfunc,resfunc,label_prefix=""):
systems = [sysfunc(s) for s in sysparam]
out = [resfunc(sys) for sys in systems]
data = [(tup[1],tup[0],s,label_prefix) for s,tup in zip(sysparam,out)]
return data
def draw_alldata(data,linspec='-'):
figure()
for t,y,s,pre in data:
plot(t,y,linspec,label=str(pre)+str(s))
legend()
grid()
seed = [0.5,1.0,2.0]
data_all = [
("impulse,var=T", get_responces(seed,lambda x : sys_1st(1,x),impulse,"T=")),
("impulse,var=K", get_responces(seed,lambda x : sys_1st(x,1),impulse,"K=")),
("step,var=T", get_responces(seed,lambda x : sys_1st(1,x),step,"T=")),
("step,var=K", get_responces(seed,lambda x : sys_1st(x,1),step,"K=")),
]
for data in data_all:
draw_alldata(data[1])
title(data[0])
xlabel("time")
ylabel("out")
5倍ルール¶
s平面上で，原点に近い極は応答が遅い = 過渡応答が最後まで残る = 支配的な応答になる. 一般的に2つの極の原点からの距離が5倍離れていたら，遠い極は無視して良い．
図は２つの一次遅れ系のカスケード結合の応答． 支配極時定数 = 10 に対して，被支配極を1,5,2(s平面上での原点からの距離10倍，5倍，2倍)
def get_dom_responces(pole_dom,pole_non,resfunc):
sys_dominant = sys_1st(1,pole_dom)
sys_nondom = sys_1st(1,pole_non)
sys_series = series(sys_dominant,sys_nondom)
systems = [sys_dominant,sys_nondom,sys_series]
out = [resfunc(sys) for sys in systems]
data = [(tup[1],tup[0],"",label) for tup,label
in zip(out,["dominant T="+str(pole_dom),"non dominant T="+str(pole_non),"series responce"])]
return data
data_all =[
("pole ratio = 10 , impulse",get_dom_responces(10,1,impulse)),
("pole ratio = 10 , step",get_dom_responces(10,1,step)),
("pole ratio = 5 , impulse",get_dom_responces(5,1,impulse)),
("pole ratio = 5 , step",get_dom_responces(5,1,step)),
("pole ratio = 2 , impulse",get_dom_responces(2,1,impulse)),
("pole ratio = 2 , step",get_dom_responces(2,1,step)),
]
for data in data_all:
draw_alldata(data[1])
title(data[0])
xlabel("time")
ylabel("out")
ボード線図¶
図はボード線図と，折点周波数の[0.1倍,1倍,10倍]でそれぞれ振幅1の正弦波入力したときの時間応答． シミュレーション点でのゲインはそれぞれ[-23db,-3db,ほぼ0db] 接点周波数での出力振幅=-3db=0.707 (入力の70% or 入力を30%減衰) 接点周波数*10での出力振幅=-23db=00.7 (入力の7%) 接点周波数/10での出力振幅≒0db=1 (入力の100%)
def draw_bode(sys):
bode(sys)
ax1 = gcf().axes[0]
lin1 = ax1.lines[-1]
ax2 = gcf().axes[1]
lin2 = ax2.lines[-1]
return ax1,ax2,lin1,lin2
T = 5
syst = sys_1st(1,T)
fn = 1/(2*pi*T)
figure()
ax1,ax2,lin1,_ = draw_bode(syst,)
"1/(Ts+1)"
ax1.plot([fn,fn*10],[0,-20],"--",color=lin1._color,label="-20 db/dec")
ax1.plot(fn,-3,'or',label="-3db")
ax2.plot(fn,-45,'og',label="-45deg")
gcf().legend()
# draw servo analyzer sim point
freq = [fn/10,fn,fn*10]
ttl = ["(1/(2*pi*T))/10","1/(2*pi*T)","(1/(2*pi*T))*10"]
for i in range(3):
ax1.plot([freq[i],freq[i]],[0,-20],'--c')
ax2.plot([freq[i],freq[i]],[-90,0],'--c')
# servo analyzer sim
N=5
resol = 20
pts = linspace(0,N,resol*N)
theta = pts*2*pi
u = sin(theta)
figure()
for i,f in enumerate(freq):
t = theta / (2*pi*f)
data = lsim(syst,U=u,T=t)
subplot(3,1,i+1)
plot(t,u,'k--',label="input")
plot(data[1],data[0],label="responce")
title("f = " +ttl[i])
gspec = gcf().axes[i].get_yticks()
tck = linspace(-1,1,11)
gcf().axes[i].set_yticks(tck,minor=True)
grid(which="both")
tight_layout()
拡張1次遅れ¶
1次遅れと同様に，時定数TとスケールKで決まるシステムは実は4パターンある
- 右肩下がり : $\frac{1}{Ts+1}$
- 右肩上がり : $Ts+1$
- 左肩下がり : $\frac{Ts}{Ts+1}$
- 左肩上がり : $\frac{Ts+1}{Ts}$
どれも折点周波数は同じ. PIDを考えるときにちょくちょく出てくるから知っとくと役に立つ．
figure()
syss = [
(tf([1],[T,1]) , "1/(Ts+1)"),
(tf([T,1],[1]) , "Ts+1"),
(tf([T,0],[T,1]), "Ts/(Ts+1)"),
(tf([T,1],[T,0]), "(Ts+1)/Ts"),
]
# ,tf([-T,0],[T,1])
for sys in syss:
ax1,ax2,lin1,lin2 = draw_bode(sys[0])
lin1.set_label(sys[1])
gcf().axes[0].plot(fn,0,"kx")
gcf().axes[0].plot([0.1*fn,10*fn],[20,-20],"k:")
gcf().axes[0].plot([0.1*fn,10*fn],[-20,20],"k:")
gcf().legend()
二次遅れ¶
一般系 $$G(s) = \frac{\omega_n^2}{s^2 + 2\zeta\omega_ns + \omega_n^2}$$ $$\frac{d^2x}{dt^2}+2\zeta \omega_n\frac{dx}{dt}+\omega_n^2x = \frac{1}{k}\omega_n^2f$$
極 $$s=-\zeta\omega_n \pm \omega_n\sqrt{\zeta^2-1}$$
応答特性は，極が実根・重根・虚根となるようなそれぞれの$\zeta$の値によって大分類できる. | $\zeta$ | 種類 | |---|---| | $\zeta=0$ | 単振動 | | $0<\zeta<1$ | 不足制動 | | $\zeta=1$ | 臨界振動 | | $1<\zeta$ | 過振動 |
周波数特性 $$|G(j\omega)| = \frac{\omega_n^2}{\sqrt{(\omega^2 - \omega_n^2)^2+(2\zeta \omega_n\omega)^2}}$$ $$ \begin{eqnarray} \angle G(j\omega) = \begin{cases} -tan^{-1}\frac{2\zeta\omega_n}{\omega_n^2 - \omega^2} \\ -90 \\ -180-tan^{-1}\frac{2\zeta\omega_n}{\omega_n^2 - \omega^2} \end{cases} \end{eqnarray} $$
ゲイン=-3dbとなる周波数(バンド幅)は $$\omega_b = \omega_n\sqrt{(1-2\zeta^2)^2+\sqrt{(1-2\zeta^2)^2+1}}$$
代表例¶
機械系 $$ \begin{eqnarray} m\frac{d^2x}{dt^2}+d\frac{dx}{dt}+kx &=& f \\ \iff \frac{d^2y}{dt^2}+2\zeta \omega_n\frac{dy}{dt}+\omega_n^2y &=& \frac{1}{k}\omega_n^2x \end{eqnarray} $$ ただし$\zeta=\frac{d}{2\sqrt{km}}$,$\omega_n=\sqrt{\frac{k}{m}}$
時間応答と極配置¶
図の説明
- 1: $\omega_n=5$に固定して$\zeta=0,0.5,1,2$の4パターンを同時プロットしたステップ応答と極配置
- 2~4: $\zeta=0,0.5,1,2$のそれぞれに対して，$\omega_n=1,3,5$の3パターンを同時プロットしたステップ応答と極配置
zeta = [0,0.5,1,2]
omega = [1,3,5]
def get_poles(sysparam,sysfunc,label_prefix=""):
systems = [sysfunc(s) for s in sysparam]
polelist = [pole(sys) for sys in systems]
data = [(p.real,p.imag,s,label_prefix) for s,p in zip(sysparam,polelist)]
return data
# basic step responce
data_all_step = [
("var=zeta, omega=5",get_responces(zeta,lambda x : sys_2nd(1,5,x),step,"zeta=")),
]
# pole placement
data_all_pole = [
("var=zeta, omega=5",get_poles(zeta,lambda x : sys_2nd(1,5,x),"zeta=")),
]
for z in zeta:
data_all_step.append((
"var=omega , zeta="+str(z),
get_responces(omega,lambda x : sys_2nd(1,x,z),step,"omega=")))
data_all_pole.append((
"var=omega , zeta="+str(z),
get_poles(omega,lambda x : sys_2nd(1,x,z),"omega=")))
for datas,datap in zip(data_all_step,data_all_pole):
draw_alldata(datas[1])
title(datas[0])
xlabel("time")
ylabel("out")
draw_alldata(datap[1],'o')
title(datap[0])
xlabel("real")
ylabel("imag")
axis("equal")
ボード線図¶
wn = 5
zeta = [0,0.2,0.5,1/sqrt(2),1,1.5,2]
syst = [sys_2nd(1,wn,z) for z in zeta]
f=figure()
for sy,z in zip(syst,zeta):
bode(sy)
ax1 = f.axes[0]
ax2 = f.axes[1]
lin1 = ax1.lines[-1]
lin2 = ax2.lines[-1]
lin1.set_label("zeta = {0:.3f}".format(z))
lin2.set_label("zeta = {0:.3f}".format(z))
if 0< z and z<1/sqrt(2):
wp = wn*sqrt(1-2*z*z)
gp = 1/(2*z*sqrt(1-z*z))
ax1.plot(wp/(2*pi),20*log10(gp),'x',label='resonant freq')
ax1 = f.axes[0]
ax1.legend()
ax2 = f.axes[1]
ax2.legend()
PID制御¶
┌─P─┐
──┼─I─┼─>
└─D─┘P,I,Dそれぞれのコントローラを並列に接続したとき，コントローラの伝達関数は以下． $$C(s)=K_p+\frac{K_i}{s}+K_ds=K_p \left(1 + \frac{1}{T_is}+ T_ds \right)$$ ただし，以下の変数変換を行った． $$T_i = \frac{K_p}{K_i} \\ T_d = \frac{K_d}{K_p}$$
なお，$T_i,T_d$はそれぞれ積分時間，微分時間と呼ばれる． 厳密には時間ではなくて，Kpの次元(または逆数の次元)に対して，時間次元を乗じたものになる． ここで，$T_i >> T_d$を仮定するとコントローラの伝達関数は以下のようになる． $$C(s) = K_p\frac{T_is+1}{T_is}(T_ds+1) = C_p(s)C_i(s)C_d(s)$$ ただし，以下の変数変換を行った． $$ \begin{eqnarray} C_p(s) &=& K_p \\ C_i(s) &=& \frac{T_is+1}{T_is} \\ C_d(s) &=& T_ds+1 \end{eqnarray} $$
つまり，PIDコントローラは$C_p,C_i,C_d$の直列接続で表すことができる．以降それぞれを比例コントローラ，微分コントローラ，積分コントローラと呼ぶ． また，$C_i,C_d$は"拡張1次遅れボード線図"で見たように，それぞれ1次遅れ系の仲間であり，$T_i,T_d$が時定数になっている．(そういう意味で微分時間，積分時間と言っているのかも)
─Cp─Ci─Cd─>ところで，上記のカスケードPIDコントローラは$T_i >> T_d$を仮定していた．これをボード線図上で考えると，折点周波数=1/2\pi T$だから，「積分コントローラ(左肩上がりの一次遅れ)の折点周波数」 << 「微分コントローラ(右肩上がりの一次遅れ)の折点周波数」になる． つまり，ボード線図上でPIDコントローラ全体のゲインは＼＿／のような凹型になるべきである． (参考書によれば折点周波数の大小関係を逆転するような熟練技もあるらしい)
また，微分ゲイン積分ゲインをそれぞれ定性的に解釈すると
- 微分を効かせる=Kdを上げる=Tdを下げる=折点周波数が低域に来る=ゲイン＿／の斜め部分が高域に対して大きくなる・効く範囲が低域に広がる
- 積分を効かせる=Kiを上げる=Tiを上げる=折点周波数が高域に来る=ゲイン＼＿の斜め部分が低域に対して大きくなる・効く範囲が高域に広がる
PD制御¶
$$C_{pi}(s) = Kp(1+Tds)$$右肩上がりの一次遅れ系と同じ形
Kp_list=[1,3]
Td_list=[0.1,0.5]
Ti_list=[10,50]
ctl = [(tf([Kp*Td,Kp],[1]),Kp,Td) for Kp,Td in it.product(Kp_list,Td_list)]
figure()
for sys,Kp,Td in ctl:
ax1,ax2,line,_ = draw_bode(sys)
line.set_label("Kp={0},Td={1}".format(Kp,Td))
fcf = 1/(2*pi*Td)
ax1.plot(fcf,mag2db(Kp),'x',color=line._color)
ax1.plot([fcf,fcf*10],[mag2db(Kp),mag2db(Kp*10)],'--',color=line._color)
ax2.plot(fcf,45,'x',color=line._color)
gcf().axes[0].legend()
不完全微分を用いたPD制御¶
微分コントローラ$C_d(s)$は入力周波数→∞に対して出力→∞になるため，例えばステップ入力などに対して極大な出力をしてしまう． 実装上は出力に対して最大値制限をかけることで回避できるが，それだとボード線図上で扱いにくい． そこで，$C_d(s)$に追加の微分抑制コントローラ$C_{d0}$をカスケード結合することを考える．
─Cd─Cd0─>微分抑制コントローラ$C_{d0}$の要件は以下．
- 低域ゲイン = $0 db$
- 高域ゲイン = $-20db / dec$
- → 右肩下がりの一次遅れ系の形をしていればよい．
- 微分コントローラよりも高域に接点周波数を持つ
- 微分コントローラのゲインが最大制限値に等しくなるような接点周波数を持つ
具体的に微分抑制コントローラの形を考える．
- 微分抑制コントローラの接点周波数$f_0$を微分コントローラの接点周波数$f_d$の$N$倍とする
- $f_0 = Nf_d = \frac{N}{2\pi T_d}$
- このとき，$N$は出力の最大制限値[db]から決める．
- $f_0$における微分コントローラゲイン = $20log_{10}(N)$
- なぜならボード線図上で20db/decだから
- さらに，$f_0$におけるPIコントローラ全体のゲイン = $20log_{10}(K_pN)$
- ↑が出力最大制限値とイコールになるように$N$を決める
- 微分抑制コントローラの時定数$T_n$を接点周波数の式から求める
- $T_0 = \frac{T_d}{N}$
- よって，微分抑制コントローラの伝達関数$C_{d0}$は以下 $$C_{d0}(s) = \frac{1}{Td/N +1}$$
なので，不完全微分を用いたPDコントローラの伝達関数は $$ \begin{eqnarray} C'_{pd}(s) &=& K_p \left(1 + \frac{Tds}{ \frac{Td}{N}s+1} \right) \\ &=& K_p \left( \frac{(1+N) \frac{T_d}{N} s +1 }{\frac{T_d}{N}s+1} \right) \\ &=& K_p \frac{\beta T'_d s +1}{T'_ds+1} = K_p \frac{ T'_d s +1}{ \alpha T'_ds+1} \\ \end{eqnarray} $$ ただし， $$T'_d = \frac{T_d}{N} \\ \alpha = 1+N \\ \beta = \frac{1}{1+N}$$
この形は$\alpha < 1 , \beta > 1$であるから，位相進み補償器に等しい. ボード線図を見ると，確かに全領域で位相>0になっている． 位相進め要素である微分コントローラに細工をしたら位相進み補償器になる，と考えるとなんとなく位相進み補償器の得体が知れてくるように思う．
Kp = Kp_list[-1]
Td = Td_list[-1]
N=5
fcf = N/(2*pi*Td)
fc = 1/(2*pi*Td)
fmax = sqrt(fc*fcf)
imcompctl = tf([Kp*Td,Kp],[Td/N,1])
ctl = tf([Kp*Td,Kp],[1])
sub = tf([Kp],[Td/N,1])
figure()
ax1,_,line,_ = draw_bode(ctl)
line.set_label("comprete diff")
ax1.plot(fc,mag2db(Kp),"x",color=line._color)
ax1.plot([fc,fc*10],[mag2db(Kp),mag2db(Kp*10)],"--",color=line._color)
ax1,_,line,line2 = draw_bode(sub)
line.set_label("O(1) sys")
line.set_linestyle(":")
line2.set_linestyle(":")
ax1.plot([fcf,fcf*10],[mag2db(Kp),mag2db(Kp/10)],"--",color=line._color)
ax1.plot(fcf,mag2db(Kp),"x",color=line._color)
ax1.plot([fcf,fcf],[mag2db(Kp),mag2db(Kp*N)],"k:")
ax1,ax2,line,_ = draw_bode(imcompctl)
line.set_label("imcomprete diff N={0}".format(N))
ax1.plot(fcf,mag2db(N*Kp),'x',color=line._color)
ax1.plot([fcf,fcf*N],[mag2db(N*Kp),mag2db(N*Kp)],"--",color=line._color)
gcf().legend()
位相進み補償器としての不完全微分¶
伝達関数は以下で与えられる $$C_{plead}(s) = K \frac{ T s +1}{ \alpha Ts+1} \ \ \ \mathrm{where} \ \alpha<1$$ この形は，上記の不完全微分を用いたPDコントローラと等価であることはここまでの議論でわかっている．
2つの折点周波数は低い方から以下で与えられる $$f_L = \frac{1}{2\pi T} \\ f_H \frac{1}{2\pi \alpha T}$$
位相を最大にすすめる周波数,角速度は以下で与えられる $$f_{max}=\sqrt{f_L f_H} = \frac{1}{2 \pi T \sqrt{\alpha}} \\ \omega_{max}=2\pi f_{max}$$
最大位相進みにおけるゲイン$G_{max}$と位相$\phi_{max}$は以下で与えられる $$ G_{max} = |C_{plead}(j\omega)|_{\omega=\omega_{max}} = \frac{K}{\sqrt{alpha}} \\ \phi_{max} = \angle C_{plead}(j\omega)_{\omega=\omega_{max}} = tan^{-1}\frac{1-\alpha}{2\sqrt{\alpha}} $$
したがって，最大位相進みから位相進み補償器を設計する場合，係数$\alpha$は以下で与えられる． $$\alpha = \frac{1-sin\phi_{max}}{1+sin\phi_{max}}$$
位相進み補償器は，特にゲインがゼロクロスする周波数領域にて位相を進める事によって位相余裕確保に利用する．
PI制御¶
$$C_pi(s) = Kp\frac{T_is+1}{T_is}$$左肩上がりの一次遅れ系と同じ形
# pi
ctl = [(tf([Kp*Ti,Kp],[Ti,0]),Kp,Ti) for Kp,Ti in it.product(Kp_list,Ti_list)]
figure()
for sys,Kp,Ti in ctl:
ax1,ax2,line,_ = draw_bode(sys)
line.set_label("Kp={0},Ti={1}".format(Kp,Ti))
fcf = 1/(2*pi*Ti)
ax1.plot(fcf,mag2db(Kp),'x',color=line._color)
ax1.plot([fcf,fcf/10],[mag2db(Kp),mag2db(Kp*10)],'--',color=line._color)
ax2.plot(fcf,-45,'x',color=line._color)
gcf().axes[0].legend()
不完全積分を用いたPI制御¶
入力周波数→0(直流)のとき，積分コントローラ出力→∞になってしまう
不完全微分と同様に抑制コントローラ$C_{i0}$をカスケード結合する．
積分抑制コントローラの要件は以下．
- 低域ゲイン = $20 db/dec$
- 高域ゲイン = $0 db$
- → 左肩下がりの一次遅れ系の形をしていればよい．
- 積分コントローラよりも低域に接点周波数を持つ
- 積分コントローラのゲインが最大制限値に等しくなるような接点周波数を持つ
具体的に積分抑制コントローラの形を考える．
- 積分抑制コントローラの接点周波数$f_0$を積分コントローラの接点周波数$f_i$の$1/N$倍とする
- $f_0 = f_i/N = \frac{1}{2\pi T_d N}$
- このとき，$N$は出力の最大制限値[db]から決める．
- $f_0$における積分コントローラゲイン = $-20log_{10}(1/N)$
- なぜならボード線図上で-20db/decだから
- さらに，$f_0$におけるPIコントローラ全体のゲイン = $20log_{10}(K_pN)$
- ↑が出力最大制限値とイコールになるように$N$を決める
- 積分抑制コントローラの時定数$T_0$を接点周波数の式から求める
- $T_0 = T_iN$
- よって，積分抑制コントローラの伝達関数$C_{i0}$は以下 $$C_{i0}(s) = N\frac{T_is+1}{NT_is+1}$$
なので，不完全積分を用いたPIコントローラの伝達関数は $$C'_{pi}(s) = K_p N\frac{T_is+1}{NT_is+1}$$
この形はN>1より，後述の位相遅れ補償器に等しい
Kp = Kp_list[-1]
Ti = Ti_list[-1]
N=5
fcf = 1/(2*pi*Ti*N)
fc = 1/(2*pi*Ti)
imcompctl = tf([Kp*N*Ti,Kp*N],[N*Ti,1])
ctl = tf([Kp*Ti,Kp],[Ti,0])
sub = tf([Kp*N*Ti,0],[N*Ti,1])
figure()
ax1,_,line,_ = draw_bode(ctl)
line.set_label("comprete integ")
ax1.plot(fc,mag2db(Kp),"x",color=line._color,label="fc = 1/(2pi*Ti)")
ax1.plot([fc,fc/10],[mag2db(Kp),mag2db(Kp*10)],"--",color=line._color)
ax1,_,line,line2 = draw_bode(sub)
line.set_label("cascaded O(1) system")
line.set_linestyle(":")
line2.set_linestyle(":")
ax1.plot(fcf,mag2db(Kp),"x",color=line._color,label="fcf = fc/N")
ax1.plot([fcf,fcf/10],[mag2db(Kp),mag2db(Kp/10)],"--",color=line._color)
ax1.plot([fcf,fcf],[mag2db(Kp),mag2db(Kp*N)],"k:")
ax1,_,line,_ = draw_bode(imcompctl)
line.set_label("imcomprete integ N={0}".format(N))
ax1.plot(fcf,mag2db(N*Kp),'x',color=line._color)
ax1.plot([fcf,fcf/N],[mag2db(N*Kp),mag2db(N*Kp)],"--",color=line._color)
gcf().legend()
位相遅れ補償器としての不完全積分¶
伝達関数は以下で与えられる $$C_{plag}(s) = K \beta \frac{Ts+1}{\beta Ts+1} \ \ \ \mathrm{where} \ \beta>1$$ この形は，上記の不完全積分を用いたPIコントローラと等価であることはここまでの議論でわかっている．
2つの折点周波数は低い方から以下で与えられる $$f_L = \frac{1}{2\pi \beta T} \\ f_H \frac{1}{2\pi T}$$
位相を最大に遅らせる周波数,角速度は以下で与えられる $$f_{min}=\sqrt{f_L f_H} = \frac{1}{2 \pi T \sqrt{\beta}} \\ \omega_{min}=2\pi f_{min}$$
最大位相遅れにおけるゲイン$G_{min}$と位相$\phi_{min}$は以下で与えられる $$ G_{min} = |C_{plag}(j\omega)|_{\omega=\omega_{max}} = \frac{K}{\sqrt{\beta}} \\ \phi_{min} = \angle C_{plag}(j\omega)_{\omega=\omega_{max}} = tan^{-1}\frac{\beta-1}{2\sqrt{\beta}} $$
位相遅れ補償器は低周波領域でのゲイン改善に用いる(位相を遅らせる代わりにゲインを上げる使い方)． 基本的に位相を遅らせることは特性改善に繋がらないので，位相を補償する用途では使わない． 低周波ゲインを改善する = 直流成分の入力に対して出力が大きくなる = 直流成分の偏差に対して出力が大きくなる = 定常偏差に対して出力が大きくなる = 定常偏差が小さくなる，と考えればOK
PID制御¶
$$C(s) = K_p\frac{T_is+1}{T_is}(T_ds+1)$$今まで出てきた要素の重ね合わせ
# pid
Kp=Kp_list[0]
Td=Td_list[0]
Ti=Td*20
cp = tf([Kp],[1])
ci = tf([Ti,1],[Ti,0])
cd = tf([Td,1],[1])
figure()
_,_,l1,_ = draw_bode(series(cp,ci,cd))
l1.set_label("pid")
_,_,l1,l2 = draw_bode(series(cp,ci))
l1.set_linestyle(":")
l2.set_linestyle(":")
l1.set_label("pi")
_,_,l1,l2 = draw_bode(series(cp,cd))
l1.set_linestyle(":")
l2.set_linestyle(":")
l1.set_label("pd")
gcf().axes[0].legend()
ノッチフィルタ¶
基本的にPID制御で振動しないように制御するためには応答性=低域ゲインを犠牲にする必要がある． (共振点でのゲインを下げる=比例ゲインを下げる=全体のゲインが下がる)
ノッチフィルタを使えば，PID特性をそのままに，特定周波数のゲインだけを下げることができる． $$C(s) = \frac{s^2+2R\zeta_N\omega_Ns+\omega_N^2}{s^2+2\zeta_N\omega_Ns+\omega_N^2}$$
- $\omega_N$ : ゲインを下げたい周波数
- $R$ : ノッチの深さ．$\omega=\omega_N$におけるゲインが$20log_{10}(R)$になる
- $\zeta_N$ : ノッチの幅． $\frac{1}{10}\zeta_N\omega_N < \omega_N < 10\zeta_N\omega_N$の範囲内のゲインが0db以下になる
def draw_bode2(sys,omega):
bode(sys,omega)
f= gcf()
ax1=f.axes[0]
ax2=f.axes[1]
line1 = ax1.lines[-1]
line2 = ax2.lines[-1]
return ax1,ax2,line1,line2
def Cn(R,z,w):
return tf([1,2*R*z*w,w**2],[1,2*z*w,w**2])
w=0.3
osamp =linspace(w/100,w*100,10000)
figure()
_,_,l,_=draw_bode2(Cn(0.01,1,w),osamp)
l.set_label("R=0.01")
_,_,l,_=draw_bode2(Cn(0.02,1,w),osamp)
l.set_label("R=0.02")
_,_,l,_=draw_bode2(Cn(0.1,1,w),osamp)
l.set_label("R=0.1")
gcf().axes[0].legend()
figure()
_,_,l,_=draw_bode2(Cn(0.1,1,w),osamp)
l.set_label("z=1")
_,_,l,_=draw_bode2(Cn(0.1,2,w),osamp)
l.set_label("z=2")
_,_,l,_=draw_bode2(Cn(0.1,10,w),osamp)
l.set_label("z=10")
gcf().axes[0].legend()
PID制御によるドローン高さ制御¶
普通の二次系とかモータ制御系はきっとググれば分かる． 前々から疑問だった，減衰項のない単振動の制御を考えてみる． 例としてドローンのZ制御を扱う.
ちなみに調べたところによると，市販ドローンは速度制御マイナーループを作り，位置制御メジャーループを構成しているようだ．
対象システム¶
要件
- 入力：z方向推力f [N]
- 出力：z方向位置z [m]
- 初期条件：z=z0 離陸済みの状態（空中からスタート）
- 制約：常に初期位置から位置偏差3m以内で制御すること
システムパラメータ
- 重量m=3[kg]
- 最大推力=6.2[kgf]
制御対象モデリング $$f-mg = f' =m\ddot{z} $$ ここで$f'$は重力と釣り合う推力
伝達関数 $$G(s)=\frac{1}{ms^2}$$
伝達関数のボード線図から次のことが分かる
- 直流入力に対してゲイン->∞なので，入力fが釣り合い推力から少しでもずれると発散する
- 常に位相=-180degなので，位相余裕を確保するために位相進め要素必須
def comp_bode(Cs,name):
sysp_ol = series(G,Cs)
sysp_cl = feedback(sysp_ol)
omega = list( linspace(2*pi*0.01 , 2*pi*10,10000) )
figure()
ax1,ax2,line1,line2 = draw_bode2(G,omega)
line1.set_label("G(s)")
ax1,ax2,line1,line2 = draw_bode2(Cs,omega)
line1.set_label("C(s)")
ax1,ax2,line1,line2 = draw_bode2(sysp_cl,omega)
line1.set_label("CL")
ax1,ax2,line1,line2 = draw_bode2(sysp_ol,omega)
line1.set_label("OL")
gcf().legend()
gcf().suptitle(name)
return Cs,sysp_ol,sysp_cl
m = 3
G = tf([1],[m,0,0])
figure()
ax1,ax2,line1,line2 = draw_bode(G)
line1.set_label("G(s)")
コントローラ要件
- 出力は最大推力$f'_{max}$=6.4-3=3.4kgf以下
- 偏差=入力は3m以下
-> 偏差3mのときに最大推力を出力するようなコントローラにする つまり，コントローラゲインについて以下が成り立つようにしたい． $$|C(s)|_{max} < 20log \left( \frac{f'_{max} g}{\Delta z} \right)$$ ※g:重力加速度
ratio_Cmax = 3.2*9.8/3
gain_Cmax = 20*log10(ratio_Cmax)
print("max gain of controller = {0:.4f} N/m = {1:.4f} db".format(ratio_Cmax,gain_Cmax))
P制御¶
いったんP制御だけで試してみる． PID制御において，DとI要素の接点周波数の間のゲイン（ボード線図の＼＿／の底の部分）はKpによって決まるので，まずはKpを最大コントローラゲインにして設計していく.
ボード線図を見ると閉ループに0.3Hz程度の共振点が現れている． また，CL伝達関数を見ても減衰項のない二次系になっているから，単振動を続けると思われる．
# P control
def tfp(Kp):
return tf([Kp],[1])
Kp = ratio_Cmax
Cp = tfp(Kp)
C,ol,cl = comp_bode(Cp,"P")
print("Kp={0:.4f}".format(Kp))
print("C(s)=")
C
print("OL(s)=")
ol
print("CL(s)=")
cl
実際にステップ入力の結果を見てみる． (Xcosを利用)
応答波形から以下が分かる
- 周期3.5sec=0.29Hzで振動する．
- CL周波数特性のゲインピーク部分に相当する．
- 入力推力は偏差と同相(比例制御だから)
- 偏差=0になったときに推力=0にしても，既に速度を持っているから通り過ぎてしまう．
- 偏差=0になる前に減速に転じる必要がある
- 例えばz=0.5(約1.5sec)の時点で加速から減速に転じたい
- つまり偏差に対して推力の位相を時間軸左側に進めたい
- 推力位相を進めるためには位相進め要素である微分器を入れる
PD制御¶
微分コントローラの接点周波数はいくつにするか？ -> CL特性の位相が落ち込むところに効いてほしいから，0.29Hz近辺を狙ってみる．
CL伝達関数の形は減衰ありの二次系になっている． このことからD要素はダンパとして働くことが分かる． 未知のシステムに対してはD要素を入れておけば事故が減って良さそう．
# PD control
def tfd(fcorner):
Td = 1/(2*pi*fcorner)
return tf([Td,1],[1]),Td
fpeek = 0.29
Cd,Td = tfd(fpeek)
C,ol,cl=comp_bode(series(Cp,Cd),"PD")
print("Td={0:.4f} -> Kd={1:.4f}".format(Td,Td*Kp))
print("C(s)=")
C
print("OL(s)=")
ol
print("CL(s)=")
cl
ステップ応答は以下のようになった
ステップ入力を見ると，無事に目標値に収束している． また，偏差と推力の関係を見ても，狙い通り推力の位相が進んでいる．
進み量は，時間応答で見るとt=2~6の4secで1周期(0.25Hzに対して，0.7sec進み=(360x0.7/4)=63deg程度． ボード線図で0.25HzのC特性を見ると，位相=40deg程度． 目測誤差もあるとして大体ボード線図通りになっている．
一方で，ステップ入力時の出力が無限大に発散してしまっている． そこで，微分の代わりに不完全微分を使うことにする
不完全微分PD制御¶
微分コントローラの折点周波数をそのままに推力の抑制を考える． 設計パラメータは，推力抑制の最大値をいくつにするか？である．
今，比例ゲインの設定は偏差=3mのときに最大推力となるように設計し，約20dbという設定値が得られている．ステップ入力の1m程度の偏差が最大(これ以上の偏差になる前に目標値収束させる)とすると，約60db程度のコントローラゲインにしても最大推力以内に収まると考えられる． そこで，(適当だけど) 不完全微分含めたPDコントローラ全体の高域ゲインを40dbとなるように設計してみる．
# PDn PDimcomplete control
def tfdn(fcorner,ratio):
# ratio : 20log(1) at low freq, 20log(ratio) at high freq
Td = 1/(2*pi*fcorner)
N = ratio
return tf([Td,1],[Td/N,1]),Td,N
outratio = 2
Cdn,_,_ = tfdn(fpeek,outratio)
C,ol,cl = comp_bode(series(Cp,Cdn),"PDn")
print("C(s)=")
C
print("OL(s)=")
ol
print("CL(s)=")
cl
ステップ入力の結果を示す．
ギリギリだけど推力が最大値内に収まってくれた． トレードオフとして，位相進め量が完全微分よりも小さくなっている(45deg程度)ボード線図見れば分かる． その結果，収束までの時間が完全微分よりも長くなっている
PID制御¶
せっかくなので積分コントローラも使うことにする．
ol
今の制御系の一巡伝達関数は上記の形． 以下のように変形できるので，積分器を2つ持っていることが分かる(2型) $$\frac{1}{s^2} \frac{b_1s+b_0}{a_1s+a_0}$$ なので，偏差の二階微分要素に対しては定常偏差が生じる つまり，加速度・推力外乱に対して定常偏差が生じる． 例として-5Nの直流外乱を付与したシミュレーションを示す．
コントローラ設計値は積分時間=接点周波数の位置
- 積分入れると位相が遅れるから，0.3Hzよりも十分低域に設定してゼロクロスに影響しないようにしたい
- 狙いたいのは直流外乱だから，低域でOK
- 接点周波数を高周波にすれば定常偏差消えるまでの時間が早くなるはずだけどどのくらいにしたら？
いったんは0.05Hzで設計してみる．
def tfi(fcorner):
Ti = 1/(2*pi*fcorner)
return tf([Ti,1],[Ti,0]),Ti
fi = 0.05
Ci,Ti = tfi(fi)
print("Ti={0:.4f} -> Ki={1:.4f}".format(Ti,Kp/Ti))
C,ol,cl = comp_bode(series(Cp,Cdn,Ci),"PDnI")
print("C(s)=")
C
print("OL(s)=")
ol
print("CL(s)=")
cl
ステップ入力結果は以下のようになった
ノッチフィルタ¶
ここまでのPID制御ではステップ入力に対して振動しながら収束している． 現状の応答性をそのままに振動抑制を考える．
PID制御系のCL特性を見ると，0.34Hzにて15dbのピークを持っている． また，ゲインピークの裾は0.34Hzを中心に±0.2Hz程度の幅を持っている.
# notch
def tfnotch(gain_db,delta_freq,freq):
# gain_db : gain = gain_db at freq
# delta_freq : gain = 0 at freq+-delta_freq
R= 10**(gain_db/20)
z=(freq + delta_freq) / freq
w=2*pi*freq
return tf([1,2*R*z*w,w**2],[1,2*z*w,w**2]),R,z,w
fpeek = 0.34
gain_db = -15
delta_f = 0.2
Cnc,R,z,w = tfnotch(gain_db,delta_f,fpeek)
_ = comp_bode(series(Cp,Cdn,Ci,Cnc),"notch")
print("R={0:.4f},z={1:.4f},w={2:.4f}".format(R,z,w))
print("C(s)=")
C
print("OL(s)=")
ol
print("CL(s)=")
cl
ステップ入力は以下のようになった
0.34Hz振動は発生していないが，0.2Hz程度の不安定振動が生じた．
コントローラ位相が0.2sec=15deg遅れている．
なんでだろう？
より高度な制御¶
フィードフォワード制御¶
フィードバック制御では一つのコントローラによって系の安定性と応答性を両立しなくてはいけない． 安定性と応答性を分けて考えられるようにするのがフィードフォワード制御．
図のようなコントローラ$C(s)$と，外乱$d$を受けるプラント$P(s)$からなる系を考える．
ブロック線図を変形すると以下のようになる．
したがって，入力$r$と外乱$d$の閉ループ伝達関数は以下の式になる． $$Gr(s) = \frac{Y}{R} = \frac{PC}{1+PC}$$ $$Gd(s) = \frac{Y}{D} = \frac{P}{1+PC}$$
ここで$Gr(s)$は目標値への追従性を表す伝達関数，$Gd(s)$は外乱への応答=安定性を表す伝達関数になっている． この式からも，コントローラ$C(s)$の設計が安定性と応答性の両方に効く事がわかる．
安定性と応答性を分離して設計するためにフィードフォワードコントローラ$Cf(s)$を図のように追加する．
入力$r$と外乱$d$の閉ループ伝達関数は以下の式になる． $$Gr(s) = \frac{Y}{R} = \frac{P(C+Cf)}{1+PC}$$ $$Gd(s) = \frac{Y}{D} = \frac{P}{1+PC}$$
この構成ならコントローラ$C(s)$で安定性設計を行い，$Cf(s)$で応答性設計が独立してできるようになっている． 特に，$Cf(s) = P(s)^{-1}$であれば，フィードフォワードコントローラへの入力(目標値) = 出力値になるため，外乱の無い系に対して完全な制御(偏差なし，遅れなし)を実現することができる．
一般的な自然界のシステムは遅れ系であることから，$P(s)$はプロパー(伝達関数の分母の次数のほうが大きい)になる．そのため$Cf(s) = P(s)^{-1}$を構成しようとすると$Cf(s)$は非プロパーなコントローラになる． 連続系における議論では非プロパーなコントローラは実現できない． なぜなら，非プロパー分子の次数が大きい=微分器が含まれる=ステップ入力に対して無限大を出力するようなコントローラを作らなきゃいけないから． そこで$Cf(s)$がプロパーになるような補償器$F(s)$を$C(s),Cf(s)$にカスケード結合することでフィードフォワード制御を実現する工夫が必要．
参考書には難しいことが書いてあるが，現代はデジタル制御ばっかりなので，離散系で差分器を使えば$Cf(s)=P(s)^{-1}$なるコントローラが近似できてしまうからあまり気にしなくて良さそう．
外乱オブザーバ¶
制御対象$P(s)$のモデル$Pn(s)$を持っているとき，モデルによる予測と，実際の観測を用いて入力を補正する方法．
図のような系を構成する． これもFF制御と同様に，$Pn(s)^{-1}$が非プロパーになってしまうが，離散系なら気にしなくていい． 本来は$Pn(s)^{-1}F(s)$がプロパーとなるような補償器$F(s)$を構成する必要がある．
カスケード制御¶
モータ制御でよく使われている． 電流制御コントローラ，速度制御コントローラ，位置制御コントローラと複数のコントローラをカスケード結合して対象を制御する方法．
設計目標
- インナーループの高帯域化 = 外乱やモデル化誤差の補償
メリット
- アウターループから見ると既知のモデル通りの動きに見えて設計が簡単になる．
- インナーループ内の外乱の影響をアウターループが受けにくい = アウターループを高帯域化できる
アンチワインドアップ¶
制御入力値は現実的には飽和を起こす． 飽和した場合，基本的に応答性は劣化する．
応答性劣化の要因の一つは積分コントローラの出力（らしい） 原理はさておき，制御量飽和時には積分コントローラへの偏差入力をゼロにして偏差が溜まらないようにしておくと良いらしい．
単純なPIDならIf文で実現してもいい． URL参考
その他細かい話¶
周波数伝達関数とはなんぞや¶
$G(s)$を$G(j\omega)$にしたものはどういう経緯で出てきたのか？ 一般のLTI(Linear Time Invariant : 時不変)システムは以下の式．
$$ a_n \frac{d^ny}{dt^n} + \cdots + a_1\frac{dy}{dt} + a_0 = x \\ \iff Y(s) (a_n s^n + \cdots + a_1s + a_0) = X(s) \\ \iff Y(s) H(s) = X(s)$$ここで$H(s)$は伝達関数$G(s)$の逆数
このシステムの周波数特性を知りたいというのは言い換えれば， 入力$x=Asin(\omega t)=Ae^{j\omega t}$のときのLTIシステム方程式の特殊解(=定常解)が知りたいということ． 特特解を$y_p$とすると，その解の形は以下のようになる． $$y_p = \alpha sin(\omega t+\beta) = \alpha e^{j \omega t + \beta}$$ LTI方程式に代入すると以下のようになる． $$ \alpha [ a_n (j\omega)^n + \cdots + a_1(j\omega) + a_0] e^{j\omega t + \beta} = Ae^{j\omega t} \\ \iff \alpha ( H(j\omega) ) e^{j\omega t + \beta} = Ae^{j\omega t} $$
※指数関数の微分 $$\frac{d}{dt} e^{j\omega t} = \frac{d}{dt}(j\omega t)e^{j\omega t} = j\omega e^{j\omega t}$$
ここでたまたま$H(j\omega)$の形が出てくる．つまり以下のような対応関係になっているということ． $$ \begin{eqnarray} \mathcal{L} \left[ \frac{d^n}{dt^n}y \right] &=& s^n Y \\ \frac{d^n}{dt^n} y &=& (j \omega)^n y \ \ \ where \ y= e^{j\omega t} \end{eqnarray} $$ なので，物理的な意味というよりも数学的な意味で$s \rightarrow j\omega$にすると便利みたい．
実際にゲインの計算までやってみる．システム方程式を更に解き進めると以下のようになる．
$$ \begin{eqnarray} \alpha ( H(j\omega) ) e^{j\omega t + \beta} &=& Ae^{j\omega t} \\ \iff \alpha (Xj+Y)e^{\beta}e^{j\omega t} &=& Ae^{j\omega t} \\ \iff \alpha [ (Xj+Y)(cos\beta + j sin\beta )]e^{j\omega t} &=& Ae^{j\omega t} \\ \iff \alpha [(Ycos\beta -X sin\beta ) + (Xcos\beta + Ysin\beta )j ]e^{j\omega t} &=& Ae^{j\omega t} \\ \iff \alpha \sqrt{X^2+Y^2} [cos(\beta-\phi)+jsin(\beta-\phi)]e^{j\omega t} &=& Ae^{j\omega t} \\ \iff \alpha \sqrt{X^2+Y^2} e^{j\omega t +(\beta-\phi)} &=& Ae^{j\omega t} \\ \end{eqnarray} $$ここで以下の変数変換を用いた $$ X=a_1\omega-a_3\omega^3+a_5\omega^5-\cdots \\ Y=a_0-a_2\omega^2+a_4\omega^4-\cdots \\ \phi=arctan\frac{Y}{-X}=-arctan\frac{Y}{X}=arctan\frac{X}{Y} $$
よって $$\alpha = \frac{1}{A}\sqrt{X^2+Y^2} \\ \beta=arctan\frac{Y}{X} $$
一般に伝達関数でよく使われるのは$G(s)=1/H(s)$だから，$\alpha,\beta$が$G(j\omega)$を使って表せないか考える．
$$G(j\omega)=\frac{1}{Xj+Y}=\frac{Y}{X^2+Y^2}+\frac{X}{X^2+Y^2}j$$だから，以下の式が成り立つ． $$\alpha = \frac{1}{A}\sqrt{Re[G(j\omega)]^2+Im[G(j\omega)]^2} \\ \beta = arctan\frac{Im[G(j \omega)]}{Re[G(j\omega)]} $$
したがって，ゲインと位相遅れは周波数伝達関数$G(j\omega)$を使って表すことができる．
モータの伝達関数¶
モータ回路方程式¶
LRC直列回路と等価としてモデリング
$$ \begin{eqnarray} v(t) &=& L \frac{di}{dt} + R i(t) + \frac{1}{C} \int_0^t i(t) dt \\ &=& L \frac{di}{dt} + R i(t) + e(t) \end{eqnarray} $$ここで，コンデンサにおける電圧降下$e$は逆起電力と呼ばれる．（なぜC成分が逆起になるのかは調査したが理由不明） ラプラス変換すると以下のようになる $$ v(s) = (sL + R) i(s) + e(s) $$
運動方程式¶
$$ \begin{eqnarray} \tau(t) &=& J\frac{d^2\theta}{dt^2}+D\frac{d\theta}{dt} \\ &=& J\frac{d\omega(t)}{dt}+D\omega(t) \end{eqnarray} $$$$ T(s) = (sJ+D)\omega(s) = \frac{1}{s}(sJ+D)\theta(s)$$電気-機械接続部の方程式¶
逆起はモータ回転速度に比例する(逆起電力定数$K_e$) $$ e(t) = K_e \frac{d\theta}{dt}$$ $$ \begin{eqnarray} e(s) &=& sK_e \theta(s) \\ &=& K_e \omega(s) \end{eqnarray} $$
モータトルクは電流に比例する(トルク定数$K_t$) $$ \tau(t) = K_t i(t)$$ $$ \tau(s) = K_t i(s)$$
式整理¶
$$ \begin{eqnarray} v(s) &=& (sL + R) i(s) + K_e \omega(s) \\ \iff i(s) &=& \frac{1}{sL+R}(v(s)-K_e\omega(s)) \\ &=& \frac{1}{sL+R}v_e(s) \end{eqnarray} $$ただし実電圧を$v_e(s)$と置いた．
ブロック線図を書くと以下のようになる
