MPC制御などでは、システムを離散時間で表した状態空間モデルを使います．その一方で立式の段階では微分方程式を立てるので、
連続時間のシステムを離散化する必要があります．
| 連続時間 |
離散時間 |
| \(\dot{x} = A x + B u\) |
\(x_{k+1} = A x_{k} + B u_{k} \) |
システムの離散化にはPythonのモジュールのScipyにscipy.signal.cont2discreteという大変便利なメソッドがありますが、v1.63の時点で離散化方法が下記のように7つもあります．
流石に多すぎてどう違うのかわからないので、ちょっと試してみましょう．
gbtだけは他の手法の一般化という位置づけなので、比較の対象は6種類です．
- gbt: generalized bilinear transformation
- bilinear: Tustin’s approximation (“gbt” with alpha=0.5)
- euler: Euler (or forward differencing) method (“gbt” with alpha=0)
- backward_diff: Backwards differencing (“gbt” with alpha=1.0)
- zoh: zero-order hold (default)
- foh: first-order hold (versionadded: 1.3.0)
- impulse: equivalent impulse response (versionadded: 1.3.0)
最初に結果だけ見たい
はい、まず実験台となるシステム、ドン．
\[
G(s) = \frac{10}{s^2 + 3s + 10}
\]
システムの出典: Convert model from continuous to discrete time - MATLAB c2d - MathWorksより、むだ時間を削除
伝達関数で表現された2次遅れ系ですね．こちらの伝達関数を状態空間の形にします．
\[
\begin{array}{rl}
\dot{x} &= \begin{bmatrix} 0 & 1 \\ -10 & -3 \end{bmatrix} x
+ \begin{bmatrix} 0 \\ 10 \end{bmatrix} u \\
y &= \begin{bmatrix} 1 & 0 \end{bmatrix} x
\end{array}
\]
従いまして、Scipyのシステム定義 \(A, B, C, D\)とする行列はこちらです．
\[
\begin{array}{rl}
A &= \begin{bmatrix} 0 & 1 \\ -10 & -3 \end{bmatrix} \\
B &= \begin{bmatrix} 0 & 10 \end{bmatrix}^{\intercal} \\
C &= \begin{bmatrix} 1 & 0 \end{bmatrix} \\
D &= \begin{bmatrix} 0 \end{bmatrix}
\end{array}
\]
離散化を dt = 0.1 [sec] として比較します．
import numpy as np
from scipy.signal import cont2discrete, lti, dlti, dstep
import matplotlib.pyplot as plt
# システム係数
A = np.array([[0, 1],[-10, -3]])
B = np.array([[0],[10]])
C = np.array([[1, 0]])
D = np.array([[0]])
# 離散時間のステップ値[sec]
dt = 0.1
# 連続時間システムの作成
l_system = lti(A, B, C, D)
## 正解データ向けステップ応答取得
t, x = system.step(T=np.linspace(0,5,100))
# 制御入力のプロット(ステップ応答)
plt.hlines(1, t[0], t[-1], linestyles='dotted')
# 正解データのプロット
plt.plot(t, x, label='Ground Truth')
# 6つの方法を試す．デフォルトのZOHが最初．
for method in ['zoh', 'bilinear', 'euler', 'backward_diff', 'foh', 'impulse']:
# システムの離散化
d_system = cont2discrete((A, B, C, D), dt, method=method)
## 離散データのステップ応答取得
s, x_d = dstep(d_system)
x_d = np.squeeze(x_d)
plt.step(s, x_d, label=method, where='post')
plt.xlim(t[0], t[-1])
plt.ylim(x[0], 1.4)
plt.ylabel('Amplitude')
plt.xlabel('Time [s]')
plt.title('Step Response')
plt.legend()
plt.savefig('000.svg')
結果は下記です．Eulerさん2名が残念な以外は、ちょっと多すぎて何言ってるかわからないですね．
Root Mean Square Errorをとってみると順位はこちら．
| 手法 |
RMSE |
| zoh |
3.80e-15 |
| impulse |
1.92e-01 |
| bilinear |
1.92e-01 |
| foh |
2.04e-01 |
| backward_diff |
3.74e-01 |
| euler |
4.50e-01 |
上位3手法を取り出すとこうなってます．ZOHが優秀です．
離散化の理論
連続時間のシステムモデル
\[
\dot{x} = A x + B u \\
\]
に対しては、Wikipedia JP
離散化で導出が詳しく書かれている通り、
\[
x_{k+1} = e^{A \Delta \text{t}} x_k + A^{-1} (e^{A \Delta \text{t}} - I) B u_k \\
\]
という厳密解が求まる．これにより
\[
\begin{array}{rl}
A_d &= e^{A \Delta \text{t}} \\
B_d &= A^{-1} (e^{A \Delta \text{t}} - I) B \\
C_d &= C \\
D_d &= D \\
\end{array}
\]
Euler
理論に対して非線形システムの局所線形化など線形時不変システム以外の条件だと \(e^{A \Delta \text{t}}\)の計算ちょっと大変すぎるよね．
\(A^{-1}\)だってしょっちゅう計算したくないし．
というわけで、
\[
e^{A \Delta \text{t}} \approx I + A \Delta \text{t}
\]
としてしまって
\[
\begin{array}{rl}
A_d &= I + A \Delta \text{t} \\
B_d &= B \Delta \text{t} \\
C_d &= C \\
D_d &= D \\
\end{array}
\]
と簡略したのがこちらの手法．
Backward
同様に\(e^{A \Delta \text{t}} \approx \left(I - A \Delta \text{t} \right)^{-1} \) と近似したもの．式を整理して、こちら．
\[
\begin{array}{rl}
A_d &= (I - A \Delta \text{t})^{-1} \\
B_d &= (I - A \Delta \text{t})^{-1} B \Delta \text{t} \\
C_d &= C (I - A \Delta \text{t})^{-1} \\
D_d &= D + C \cdot B_d \\
\end{array}
\]
離散化した状態方程式・観測方程式は添字が少しずれるので注意．
\[
\begin{array}{rl}
x_{k} &= A_d \color{red}{x_{k-1}} + B_d u_{k} \\
y_{k} &= C_d \color{red}{x_{k-1}} + D_d u_{k} \\
&= C (A_d x_{k-1} + B_d u_{k}) + D u_{k} \\
&= C x_{k} + D u_{k} \\
\end{array}
\]
Bilinear / Tustin
同様に
\[
\begin{array}{rl}
e^{A \Delta \text{t}} &= \frac{e^{A \Delta \text{t}/2}}{e^{-A \Delta \text{t}/2}} \\
& \approx \left(I + \frac{1}{2} A \Delta \text{t} \right) \left(I - \frac{1}{2} A \Delta \text{t} \right)^{-1}
\end{array}
\]
と近似したもの
Impulse
前記の厳密解を導出する過程では
\[
x_{k+1} = e^{A \Delta \text{t}} x_k + \int_0^{\Delta \text{t}} e^{Av} B u_k \text{d} v \\
\]
を計算している．この際、\([0, \Delta \text{t} ]\) で、\(u_k\) 一定としているが、
ここを一定ではなくインパルス \( \delta(v) u_k \) とする．
すると \(B_d \) の計算結果が変わり、
\[
\begin{array}{rl}
A_d &= e^{A \Delta \text{t}} \\
B_d &= e^{A \Delta \text{t}} \Delta \text{t} B \\
C_d &= C \\
D_d &= D \\
\end{array}
\]
ZOH: Zero Order Hold
デジタル制御：最初の最初 ①で触れたZOHそのもの．制御入力uは、通常、制御周期間では定数(ZOH)であって \( \dot{u} = 0 \) といえる．
この条件を使ってシステム方程式を書き直すと
\[
\begin{bmatrix} \dot{x} \\ \dot{u} \end{bmatrix} = \begin{bmatrix} A & B \\ 0 & 0 \end{bmatrix} \begin{bmatrix} x \\ u \end{bmatrix}
\]
改めてベクトルrと行列Mを下記のように定義すると
\[
\begin{array}{rl}
r &= \begin{bmatrix} x & u \end{bmatrix}^{\intercal}\\
M &= \begin{bmatrix} A & B \\ 0 & 0 \end{bmatrix}
\end{array}
\]
上記方程式は、\( \dot{r} = M r \) であって、この方程式の解は、\( r(t) = e^{Mt} r(0) \) .
これより、\(t_{k+1} = t_{k} + \Delta \text{t} \) と置くと、
\(x_{k+1} = x(t_{k+1}) \) であるので
\[
\begin{bmatrix} x_{k+1} \\ u_{k+1} \end{bmatrix} = e^{M \Delta \text{t}} \begin{bmatrix} x_k \\ u_k \end{bmatrix}
\]
以上を開いて
\(M_d = e^{M \Delta \text{t}} \) を計算してやると
\[
M_d = \begin{bmatrix} A_d & B_d \\ 0 & I \end{bmatrix}
\]
となるので、これをありがたく使う．
このようにステップ応答でZOHが精度無双したのは、導出の前提(制御入力ZOH)とステップ応答が同じ解を求めていることによる．
伝達関数の離散化についてはZ変換表も使えるので、プラントモデルの離散化：Z変換表も是非どうぞ．
FOH: First Order Hold
(レビュー必要)
ZOHと似たようなことをする. まずZOHで求めた解を眺めると
\[
\begin{bmatrix} x_{k+1} \\ u_{k+1} \end{bmatrix} = \begin{bmatrix} A_d & B_d \\ 0 & I \end{bmatrix} \begin{bmatrix} x_k \\ u_k \end{bmatrix}
\]
\(u_{k+1}=u_{k}\) となっていることがわかる．
DDP(Differential Dynamic Programming)などの制御だと低次制御(比例制御)と高次制御(MPC制御)のカスケードにすることがあり、
\(u_{k+1}=u_{k}\) とは言えなくなるほど誤差が大きい場合がある．
そこで、\( \dot{u} = 0 \) をちょっと見直す．
もし、 \(k\) の時点で \(u_{k+1}\) の値がわかっているなら、 \(\dot{u} = \frac{u_{k+1} - u_k}{\Delta \text{t}} \) と置いた方が精度が高まる．これがFOHの考え方．
ところが実際にはわからない事が多いので、scipyの実装では固定的に \(\dot{u} = \mathbb{1}/ \Delta \text{t} \) と置いている．
これは投機的な予測になるので当たれば精度が高まるが、その逆ではズレが大きくなると思う．
また、\(u_k\) が充分に大きいと \(u_{k+1} \approx u_k + 1 \) となり、ZOHとあまり変わらないのではという予感も．
兎にも角にも、Scipyの実装にならえば、
\[
\begin{array}{rl}
\begin{bmatrix} \dot{x} \\ \dot{u} \\ 0_n \end{bmatrix} &= \begin{bmatrix} A & B & 0 \\ 0 & 0 & \mathbb{1}_{n \times n}/ \Delta \text{t} \\ 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} x \\ u \\ \mathbb{1}_n \end{bmatrix} \\
s.t. & u \in \mathcal{R}^{n}
\end{array}
\]
と方程式をたてて、あとはZOHと同じフローで \(A_d, B_d \) を得る．
参考
