はじめに
python で作った PID 制御の関数がどの程度の精度で動作するのか、理論的な計算値と比較してみました。
前提条件
PID 制御器の伝達関数 $C(s)$ は以下とする。
$$C(s) = K_p + K_i/s + K_ds$$
制御対象の伝達関数 $G(s)$ は以下とする。
$$G(s) = 1/s$$
これは 測定値 $y(t)$ の変化率が 操作量 $u(t)$ に比例するというモデルに相当する。
$$\dot{y}(t) = u(t)$$
以上から、システム全体としての伝達関数 $K(s)$ は以下となる。
$$K(s) = \frac{C(s)G(s)}{C(s)G(s) + 1}$$
これに ステップ関数 $R(s) = 1/s$ を入力した時の応答 $Y(s)$ は以下となる。
$$Y(s) = K(s)R(s) = \frac{C(s)G(s)}{C(s)G(s) + 1}\cdot\frac{1}{s} = \frac{K_ds^2 + K_ps +Ki}{(K_d + 1)s^2 + K_ps +K_i}\cdot\frac{1}{s}$$
計算方法
以下に示す 3 通りの方法で計算して、結果を比較する。
方法 1: 数学的に解く方法
$Y(s)$ を逆ラプラス変換して $y(t)$ を求められるが、上記の条件であれば数式の形で比較的簡単に解くことができる。
そこで、まずこの方法で $y(t)$ を求めた結果を $Y_1$ としてプロットする。
方法 2: Python-Control を使う方法
python のパッケージで Python-Control という便利な機能があるので、これを使って計算してみることにする。
これは上記で定義した 伝達関数 $C(s)$, $G(s)$ を引数として関数を呼び出すと $y(t)$ が (数値計算結果として) 求められるというもので、
この結果を $Y_2$ とする。
方法 3: 近似的な数値計算
$\Delta{t}$ を十分微小な時間変化と仮定して、微分・積分を近似値として求める。
以下では、制御対象の測定量 $y(t)$ を離散的な数列 $y_i$ として表現する。
同様に、偏差量 $e(t)$ および 操作量 $u(t)$ を各々 $e_i$, $u_i$ と表現する。
また、目標値 $r(t)$ は $r(t) = 1$ とする。
経過時間 $t$ は $t = i\Delta{t}$ であるため、
$$y_i = y(i\Delta{t})$$
と相互変換できる。
初期条件として $i = 0$ のとき
$$y_0 = 0$$
$$e_0 = 0$$
とする。
$i \geqq 1$ のときは以下の手順で $y_i$ を計算する。
まず、フィードバックされた測定量 $y_{i-1}$ から偏差量 $e_i$ を求める。
$$e_i = 1 - y_{i-1}$$
次に、操作量 $u_i$ を以下の近似式によって求める。
$$u_i = K_p + K_i\bigg(\sum_{k=1}^{i} e_k\bigg)\Delta{t} + K_d\frac{e_i - e_{i-1}}{\Delta{t}}$$
すると、 $\dot{y}(t) = u(t)$ という条件から $y_i$ は以下の式で近似できる。
$$y_i = y_{i-1} + u_i\Delta{t}$$
以上の手順を繰り返すことで、逐次的に $y_i$ を求めることができる。
この方法によって求めた結果を $Y_3$ とする。
実行結果
実際に行ってみると、思ったより合っていることが確認できました。
$Y_1$, $Y_2$ は厳密な解なので、ある意味、合ってないと逆におかしいです。
$Y_3$ はあくまで近似計算なので多少誤差があることを予想しておりましたが、以外と正確なことがわかります。
実行方法について
pidctest.py を Jupyter Lab のセルに読み込んで実行します。
また、Python-Control を使うため、以下の通りパッケージをインストールの必要があります。
conda install -c conda-forge control slycot
または
pip install control slycot
以下にJupyter Labで実行するpid-control-test.ipynbのリンクを記載します。