面倒な双一次変換の計算はPythonにやってもらおう

上の写真の系について，運動方程式は，

$$ M \ddot{x} + B \dot{x} + K x = f $$

ですね．ではこの物体を人間が握るものとして，アドミタンス提示型のインピーダンス制御をかけるましょう．つまり，アクチュエータを用いて任意の$M,B,K$を表現したいのです．入力を$F(s)$，出力を$\dot{X}(s)$とすると，伝達関数$G$は

$$ G = \frac{s}{Ms^2 + Bs + K} $$

という2型になります．(本当は剛性いらないので1型にできるが...)
そろそろ本題に入りますが，目標速度$\dot{x}$はどうやって計算しますか? 僕は双一次変換を使います．今回の例の正解は，

$$
\dot{x}[k] = \frac{1}{2BT_s+KT_s^2+4M} { 2T_s(f[k-2]+f[k]) + (8M-2KT_s^2)\dot{x}[k-1] \

• (2BT_s-KT_s^2-4M)\dot{x}[k-2] }
  $$

ですが，こんな面倒な計算を手でしたくはありません．そこでPythonにやってもらおうということです．ipynbファイルを以下のリポジトリに置きました．2型までですが，伝達関数を入力するだけで双一次変換を一瞬で計算してくれます．SympyはSymbolicな演算をしてくれるだけでなく，簡単にTex表記もしてくれるので非常に便利ですね．