Octaveでシミュレーションをやってみたい
そもそものモチベーションは植物生理学をもとにしたシミュレーションをやってみたいからだ
作業環境のOctaveのバージョン
Octave
>> version
ans = 6.3.0
Octaveを使ってシミュレーションした記事はこちらを参考にした
なお、本記事に出てくるOctaveのコードは、上記記事から引用したものである。
シミュレーションは、入力に対して出力がどうなるか、を計算して結果を得るのが目的なので、まずは伝達関数によるシンプルなシミュレーションを行うことにした。
そのパッケージであるControlを手に入れる。
OctaveにControlパッケージをインストールする
Octave
>> pkg install -forge control
伝達関数は、制御工学において、入出力の関係を表す関数で、その入出力システムの挙動や安定性を評価するもの。
$P(s)$のように表現し、$s$は複素数。どのような入力に対しても、伝達関数を用いれば出力が求まるためには、変数は時間$t$ではなく、複素数によって表現されるs領域$s$を用いる。
つまり、(出力)=(伝達関数)×(入力)の形式である。
出力を制御するためには、伝達関数が肝であることが分かる。
基本関数の特性については、東北大学の張山先生による以下のテキストが参考になる：
基本伝達関数の特性
「一次遅れ＋むだ時間」系のモデルの記述式は以下：
P(s)=\frac{K}{1+Ts} e^{-Ls}
それぞれの文字の意味は次の通り：
| 文字 | 意味 |
|---|---|
| $P$ | 伝達関数 |
| $s$ | s領域 |
| $K$ | ゲイン |
| $T$ | 時定数（小さいほど入力に対して出力が早くなる） |
| $L$ | むだ時間 |
まずはモデリングをする上で、定数を以下の通り設定することにする
| 文字 | 値 |
|---|---|
| $K$ | 0.5 |
| $T$ | 40.0 |
| $L$ | 8.0 |
Octave
% 推定した制御対象のモデル
K = 2.0;
T = 40.0;
L = 8.0;
P = tf(K,[T,1])*tf([-L/2,1],[L/2,1]); % Octaveの場合（一次で近似）
制御器の設計をする
不完全微分付きPIDコントローラを以下の式で表現する
$$C(s)=k_p\left( 1+\frac{1}{{T_i}s}+\frac{T_ds}{1+0.125 {T_d}s} \right)$$
| 文字 | 意味 |
|---|---|
| $C$ | PID制御伝達関数 |
| $s$ | s領域 |
| $k_p$ | 比例ゲイン（$p:$Proportional） |
| $T_i$ | 積分時間（$i:$Integral） |
| $T_d$ | 微分時間（$d:$Differential） |
PIDコントローラの設計式の詳細は日本工業新聞の以下の記事を参考にした：
Octaveコードは以下
Octave
kp = 0.6*T/K/L;
Ti = 1.2*T;
Td = 0.5*L;
C = kp*(1+tf(1,[1,0])/Ti+Td*tf([1,0],[0.125*Td,1]));
ここで、シミュレーションの結果をグラフに表示すると
Octave
>> pkg load control
% 推定した制御対象のモデル
K = 2.0;
T = 40.0;
L = 8.0;
P = tf(K,[T,1])*tf([-L/2,1],[L/2,1]); % Octaveの場合（一次で近似）
% PIDコントローラの設計
kp = 0.6*T/K/L;
Ti = 1.2*T;
Td = 0.5*L;
C = kp * ( 1+tf(1,[1,0])/Ti + Td * tf([1,0],[0.125*Td,1]))
margin(P*C);
step(P*C/(1+P*C));
step(1/(1+P*C));
グラフが表示された。
