伝達関数の差分化の話（その３）

　そろそろ終わりにしたい気もするけど、やっぱりグラフの比較くらいほしいのと、自分流+=戦法の話をもう少し書いておく。もうちっとだけ続くんじゃ。

　さて、(1)の実装として、(2)と(3)を比較してみよう。

$\displaystyle Y(s) = \frac{1}{1+sT}U(s)$ ...... (1)

$\displaystyle y_n=(1-a)y_{n-1}+a u_{n-1}$ ...... (2)

$\displaystyle y_n=(1-2b)y_{n-1}+b u_n+b u_{n-1}$ ...... (3)

（ $a=\frac{ΔT}{T}$ 、 $b=\frac{ΔT}{2T+ΔT}$ とおいた）

　(2)のいい所は $u_n$ が出てこないところ。

タイマ「CPUさん、時間です！」

CPU「よし、 $y_n$ を計算するか。ADCよ、 $u_n$ はまだか？」

ADC「すいません！ただ今取ってきますので少々お待ちください！」

CPU「ちっ、これだからのろまなやつと組むのは嫌いなんだ。」

みたいな事にならずにCPUの手元に残ってる一個前のyとuの値で済む。一方で(3)のいい所は何にも考えずに実装できるわりにはそれなりに使えるところ。え？茶番よりグラフ？はい。

ΔＴが細かい場合と粗い場合。双一次変換は一次遅れ系ならほぼ誤差は気にしなくてよさげ。+=戦法はΔT/Tが0.1くらい以下ならいいかなって感じ。

　以上で良いかと思ったんだけどエクセルで遊んでたら色々思い出してきたのでもうちょっと書く。以下の二次遅れ系を実装してみる。

$\displaystyle Y(s) = \frac{ω^ 2}{s^ 2 + 2ξωs + ω^ 2}U(s)$ ...... (4)

双一次変換は何も考えずに(4)に(5)をぶち込んでぐつぐつ煮て(6)を得ましょう。

$\displaystyle s ≒ \frac{2}{ΔT}\frac{z-1}{z+1}$ ...... (5)

$y_n = \frac{2(1-c)}{1+c+d}y_{n-1} +\frac{d-c-1}{1+c+d}y_{n-2}+ \frac{c}{1+c+d}(u_{n}+2u_{n-1}+u_{n-2})$ ...... (6)

（ $c=\frac{ΔT^ 2ω^ 2}{4}$ 、 $d=ξωΔT$ とおいた）

頭は使わず手だけ動かせば良いのがいい所。

　一方で+=戦法はそのまま代入しちゃいけない。なぜか？ $s$ が２個以上出てくる式にそのまま入れると $z$ の多項式になって $y_{n-2}$ とかが出てきてしまう。一個前のだけで書ける事だけが利点なのでちょっと頭を使おう。一言で書くと１元高次の方程式を多元１次の方程式にすればいい。

$\displaystyle \ddot{y} + 2ξω\dot{y} + ω^ 2 y = ω^ 2 u$ ...... (4')

とりあえずyの時間微分をxとおく。

$\displaystyle \begin{cases} \; \dot{x}= -2ξωx -ω^ 2 y + ω^ 2 u \\ \; \dot{y}=x\end{cases}$ ...... (4'')

$\displaystyle \begin{cases} \; x_n=x_{n-1} -2ξωΔTx_{n-1} -ΔTω^ 2 y_{n-1} + ΔTω^ 2 u_{n-1} \\ \; y_n=y_{n-1}+ΔTx_n \end{cases}$ ...... (5)

できた。(5)の２行目の最後が $x_{n-1}$ じゃなくで $x_{n}$ なのはずるい気もするけど、２行目みたいに微分だけの時は１行目の結果が使えるからできる。当然お互いにお互いの行の結果を使うようなときは全部 $_{n-1}$ にするしかないけど。

どっちもΔTω=0.3くらいが限界ぽい。周波数の補正とかすれば使えなくもないけど。

添え字が全部 $_{n-1}$ のずる無し版+=戦法だとΔTω=0.03くらいが限界ぽい。そして周波数じゃなくて振幅がずれるのか。なるほど。

（ちなみに勘の良い読者は「(4'')って状態方程式ぽくね？」とか思ったかもしれない。実はその通りで伝達関数を状態方程式にして一階微分方程式にしちゃえば+=戦法以外もやりたい放題(4次ルンゲクッタとか)なんだけど余白が足りないので今回はここまで。もう飽きたけど需要があれば続くかも。）