シミュレーション
この資料について
本資料ではシステム制御の学習や研究に必要な動的システムのシミュレーション手法について説明します．シミュレーション環境は基本的にGNU Octaveを対象とします．
シミュレーションとは，ある現象をそれを模擬した装置や計算機を利用して擬似的に再現し，予測や評価に利用することといえます．最近ではデジタル計算機を用いた数値シミュレーションが一般的です． 数値シミュレーションを行うためには，現象の数学的モデルが必要になります．特に工学的な現象の場合，時間に関する微分方程式で表現できることが多いでしょう．その微分方程式を解くことができればシミュレーションが行えます．しかし，微分方程式の解析的な解（数式で表現される解）を求めることは一般的には難しいため，数値解（時間に対する数値の列）を求める手法が用いられます．
この資料ではまず，最も基本的なオイラー法やより精度の高いルンゲ・クッタ法による微分方程式の数値解法をOctaveでのスクリプト例を示しながら説明します．次に，システム制御では基本となる線形システムの状態空間表現について述べ，オイラー法やルンゲ・クッタ法でシミュレ－ションするOctaveスプリプトを示します． 最後に，より高度な数値解法が実装されているOctaveのlsode関数の利用方法について説明します．基礎はいいので手っ取り早くOctaveによるシミュレーションをしたいのなら，「Octave関数を利用したシミュレーション」の節だけを読めばいいでしょう．
オイラー法によるシミュレーション
一階常微分方程式
微分方程式とは，ある未知関数とその関数の導関数との関係が表されている方程式です．そのうち特に1変数の導関数で表されているものを常微分方程式といいます．システム制御に関わるシミュレーションでは時間に関する導関数で表される場合がほとんどでしょうから，常微分方程式が対象になります．さらに，導関数の階数（次数）が微分方程式の階数となります．
ということで，これから取り扱う最も簡単な微分方程式は以下のような一階の常微分方程式となります．ここで&math(x(t));が未知関数，&math(x_0);はその初期値です．
- math(\frac{dx(t)}{dt} = f(x(t), t) \, , \quad x(0) = x_0)
もちろん実際の問題はより高階の常微分方程式が必要になりますが，この一階の微分方程式の数値解法を応用することで対応できます．
オイラー法
微分方程式の数値解法の基本は，ある時刻$t_k$から微少な刻み幅&math(\Delta t);（サンプリング時間またはサンプリング周期ともいう）だけ時間がたった後の関数値&math(x(t_k + \Delta t));を逐次求めていくことです．これをいかに精度よく（できるだけ少ない計算量で）求めるかが数値解法の善し悪しと言えます．
オイラー法は以下のように&math(x(t_k + \Delta t));のテイラー展開の最初の2項だけに基づきます．
- math("\begin{array}{l l l} x(t_k + \Delta t) &=& x(t_k) + \frac{dx}{dt}\bigg|_{t=t_k} \, \Delta t + \frac{1}{2} \frac{d^2x}{dt^2}\bigg|_{t=t_k} \, \Delta t ^2 + \cdots \\ &\approx& x(t_k) + \frac{dx}{dt}\bigg|_{t=t_k} \, \Delta t \\ &\approx& x(t_k) + f(x(t_k), t_k) \, \Delta t \end{array}")
ここで，表記上の簡便性から，&math(x(t_k));を&math(x_k);と表記することにし，&math(t_k + \Delta t);を&math(t_{k+1});, &math(x(t_k + \Delta t));を&math(x_{k+1});などと表記することにします．この時オイラー法は以下のように表すことができます．
- math("x_{k+1} = x_k + f(x_k, t_k) \, \Delta t")
オイラー法による近似をグラフで示すと図のようになります．
- ref(Euler.png)
次のサンプリング時刻における&math(x(t));の値，すなわち&math(x_{k+1});を接線で一次近似して求めていると解釈できます．もとの微分方程式から，接線の傾きが&math(f(x_k, t_k));ですから，それに&math(\Delta t);を乗じた物を&math(x_{k});に加えることで&math(x_{k+1});を求めています．図をみると誤差（error）が大きくて，本当に使えるのか心配ですが，&math(\Delta t);をできるだけ小さくすれば誤差を少なくできます．しかし，&math(\Delta t);を小さくすると，ある期間（例えば10秒間とか）のシミュレーションをするのにそれだけ計算回数が増えることになるので注意が必要です．
図では&math(x_k);が真値と一致している場合で書きましたが，実際はそうとも限りません．初期値&math(x_0);からスタートして，オイラー法を繰り返し使った近似計算は次図のようなイメージになります．&math(\Delta t);をできるだけ小さくとって，真値から離れすぎないようにする必要があります．なお，この図では近似のための接線を真値のグラフ&math(x(t));に接するように引いていますが，&math(x_k);が真値からずれている場合は図の通りにはなりません．あくまでイメージです．
- ref(Euler2.png)
Octaveによる実現
次の微分方程式を例題に，Octaveでオイラー法によるシミュレーションを行うスクリプトを考えてみます．
- math(\frac{dx}{dt} = - x + 1, \quad x_0 = 0)
オイラー法を適用すれば
- math(x_{k+1} = x_k + (-x_k + 1) \, \Delta t)
とかけるので，繰り返しこの計算を行うOctaveスクリプトの例は以下のようになります．
- オイラー法による1階微分方程式のシミュレーション
dt = 0.01; % サンプリングタイム
x_k = 0; % xの初期値
x = []; % シミュレーション結果格納用変数
tim =[]; % シミュレーション時刻格納用変数
for t = 0 : dt : 10 % 時刻tは0から10までdt刻み
x_k1 = x_k + ( - x_k + 1) * dt ; % x_{k+1} の計算
x = [x x_k]; % xの末尾にx_kを追加
tim = [tim t]; % timの末尾にtを追加
x_k = x_k1; % 次回の計算のため x_{k} ← x_{k+1}に繰り上げ
endfor
grid on;
axis([0 10 0 1.5]);
plot(tim, x);
xlabel('time (s)');
ylabel('x');
Octaveにおける変数名の都合上，&math(x_k);をx_k，&math(x_{k+1});をx_k1としています． また，シミュレーション結果と時刻を格納する変数をそれぞれx，timとして用意し，&math(x_k);と&math(t);を行の末尾に追加しています．また，forループで繰り返す際に，次回の計算のために今回の&math(x_{k+1});を次回の&math(x_k);として繰り上げておくことが必要です．
高次システムのシミュレーション
高次微分方程式の取り扱い
これまでは1階(次)の微分方程式を扱ってきましたが，現実には2次以上のシステムのシミュレーションが必要になることが多いと思います．その場合は，変数変換をすることで連立多元1階常微分方程式に変形して考えます．
例えば，次のような微分方程式を考えます．これは，長さ&math(l);の軽い棒の先に質量&math(m);の重りがついた単振り子を初期角度&math(\theta_0);の位置から初速度0で放す状況を表しています（下図参照）．ここで振り子は空気その他から速度に比例した抵抗を受けるものとし，その係数を&math(c);としています．
- math(\frac{d^2 \theta(t)}{d t^2} = - \frac{g}{l} \sin \theta(t) - \frac{c}{m} \frac{d\theta(t)}{dt} \quad , \theta(0) = \theta_0 \quad , \frac{d\theta(t)}{dt} \bigg|_{t=0} = 0)
- ref(pendulum.png)
このような2次のシステムの場合，次のような2個の変数を新たに定義します．これらを状態変数と呼びます．
- math(\begin{array}{l l l}x_1(t) &=& \theta(t)\\x_2(t) &=& \frac{d \theta(t)}{dt}\end{array})
状態変数を用いてもとの微分方程式を表すと以下のように1階の微分方程式が2本連立したものになります．
- math(\begin{array}{l l l}\frac{d x_1(t)}{dt} &=& x_2(t) \quad , x_1(0) = \theta_0\\\frac{d x_2(t)}{dt} &=& - \frac{g}{l} \sin x_1(t) - \frac{c}{m} x_2(t) \quad , x_2(0) = 0\end{array})
さらに，2つの状態変数を組にして，
- math(\boldsymbol{x}(t) = \left(\begin{array}{c}x_1(t)\\x_2(t)\end{array}\right) \quad , \ \boldsymbol{x}_0 = \left(\begin{array}{c}\theta_0\\0\end{array}\right))
のようにベクトル&math(\boldsymbol{x}(t));で表現することで，最初の1階常微分方程式と同じ見た目になります．
- math(\frac{d \boldsymbol{x}(t)}{dt} = \boldsymbol{f}(\boldsymbol{x}(t), t) \quad , \ \boldsymbol{x}(0) = \boldsymbol{x}_0)
ここで，&math(\boldsymbol{f}(\boldsymbol{x}(t), t));は以下となります．
- math(\boldsymbol{f}(\boldsymbol{x}(t), t) = \left(\begin{array}{c}x_2(t)\\- \frac{g}{l} \sin x_1(t) - \frac{c}{m} x_2(t) \end{array}\right))
オイラー法も形式的に次のように表すことができます．
- math(\boldsymbol{x}_{k+1} = \boldsymbol{x}_k + \boldsymbol{f}(\boldsymbol{x}_k, t_k) \, \Delta t)
この状態変数表現に基づいたOctaveによるシミュレーション・スクリプトは次のようになります． 基本的には1次のシステムに対するオイラー法のスクリプトと同じですが，x_k やx_k1が2行×1列の列ベクトルになっていることが違いです．それに伴い，結果を格納する変数xが最終的に2行×サンプリング点数の行列になります．よって，シミュレーション結果として&math(\theta);だけをグラフ表示するには1行目だけを取り出す必要があり，そのため，plot関数の引数をx(1, : )としています．
なお，この例ではサンプリングタイム&math(\Delta t);を0.01秒としていますが，これを0.05秒にするとシミュレーション結果が大きく異なることに気がつくと思います．さらに&math(\Delta t);を大きくしていくと，それは誤差が大きくなるというようなレベルではなく，全く異なったシミュレーション結果になってしまいます．このようにサンプリングタイムの選択は重要で，基本的には十分小さくとる必要があります．
- オイラー法による単振り子のシミュレーション
dt = 0.01; % サンプリングタイム
m = 1; % 質量 1 kg
l = 0.5; % 振り子長さ 0.5 m
c = 0.5; % 空気抵抗係数
x_k = [30/180*pi ;
0 ]; % 初期角度30度, 初期速度0
x = []; % シミュレーション結果格納用変数
tim =[]; % シミュレーション時刻格納用変数
for t = 0 : dt : 10 % 時刻tは0から10までdt刻み
x_k1 = x_k + [ x_k(2) ;
-9.8/l * sin(x_k(1)) - c/m * x_k(2) ] * dt; % x_{k+1} の計算
x = [x x_k]; % xの末尾にx_kを追加
tim = [tim t]; % timの末尾にtを追加
x_k = x_k1; % 次回の計算のため x_{k} ← x_{k+1}に繰り上げ
endfor
grid on;
axis([0 10 -1 1]);
plot(tim, x(1, :)); % x_1だけをグラフ表示
xlabel('time (s)');
ylabel('theta (rad)');
ルンゲ・クッタ法によるシミュレーション
2次のルンゲ・クッタ法（ホイン法）
オイラー法は&math(t_k);における導関数の値&math(f(x_k, t_k));を使って&math(x_{k+1});を近似していました．ここでは．より精度を高めるために，複数点における導関数の値を利用することを考えます．
先の図のように&math(t_k);における接線の傾き&math(f(x_k, t_k));だけを考えるのではなく，&math(t_{k+1});における傾き&math(f(x_{k+1}, t_{k+1}));も利用して，それらの平均をその区間での傾きとして利用すれば，精度が高まりそうです．その場合，&math(x_{k+1});は以下のように考えられます．
- math(x_{k+1} = x_k + \frac{1}{2} ( f(x_k, t_k) + f(x_{k+1}, t_{k+1}) ) \Delta t)
しかし，&math(x_{k+1});を求める計算に，&math(x_{k+1});が含まれているのは無理があるので，&math(x_{k+1});をオイラー法で近似した
- math(x_{k+1} = x_k + f(x_k, t_k) \Delta t)
を代入すれば，
- math(x_{k+1} = x_k + \frac{1}{2} ( f(x_k, t_k) + f(x_k + f(x_k, t_k) \Delta t, t_{k+1}) ) \Delta t)
とかけます．さらに，逐次計算しやすいように&math(k_1);, &math(k_2);という変数を定義して整理すれば，
- math(\begin{array}{l l l}k_1 &=& f(x_k, t_k) \Delta t \\k_2 &=& f(x_k + k_1, t_{k+1}) \Delta t \\x_{k+1} &=& x_k + \frac{1}{2} ( k_1 + k_2 )\end{array})
となります．これは2次のルンゲ・クッタ法の一種で，ホイン法または修正オイラー法と呼ばれるものです．2次のルンゲ・クッタ法としては，導関数を取る点として区間の中点を考えた中点法と呼ばれるものもあります．
なお，オイラー法と同様に高階の微分方程式に対しても状態空間表現することでもちろん適用できます．
- ref(Rungehoine.png)
- 2次のルンゲ・クッタ法（ホイン法）の原理
前節でオイラー法を用いてシミュレーションした1階の微分方程式と，2階の微分方程式をそれぞれホイン法でシミュレーションするOctaveスクリプトは以下のようになります．2階の微分方程式の場合は&math(k_1);, &math(k_2);もベクトルになることに注意してください．
また，オイラー法に比べてサンプリングタイムを大きくとっても精度よくシミュレーションできることを確かめてみてください．
- ホイン法による1階微分方程式のシミュレーション
dt = 0.01; % サンプリングタイム
x_k = 0; % xの初期値
x = []; % シミュレーション結果格納用変数
tim =[]; % シミュレーション時刻格納用変数
for t = 0 : dt : 10 % 時刻tは0から10までdt刻み
k1 = ( - x_k + 1) * dt; % k1 の計算
k2 = ( - (x_k + k1) + 1) * dt; % k2 の計算
x_k1 = x_k + ( k1 + k2) / 2 ; % x_{k+1} の計算
x = [x x_k]; % xの末尾にx_kを追加
tim = [tim t]; % timの末尾にtを追加
x_k = x_k1; % 次回の計算のため x_{k} ← x_{k+1}に繰り上げ
endfor
grid on;
axis([0 10 0 1.5]);
plot(tim, x);
xlabel('time (s)');
ylabel('x ');
- ホイン法による単振り子のシミュレーション
dt = 0.01; % サンプリングタイム
m = 1; % 質量 1 kg
l = 0.5; % 振り子長さ 0.5 m
c = 0.5; % 空気抵抗係数
x_k = [30/180*pi ;
0 ]; % 初期角度30度, 初期速度0
x = []; % シミュレーション結果格納用変数
tim =[]; % シミュレーション時刻格納用変数
for t = 0 : dt : 10 % 時刻tは0から10までdt刻み
k1 = [ x_k(2) ;
-9.8/l * sin(x_k(1)) - c/m * x_k(2) ] * dt; % k1 の計算
x_tmp = x_k + k1; % x_k1の近似値を一旦計算
k2 = [ x_tmp(2) ;
-9.8/l * sin(x_tmp(1)) - c/m * x_tmp(2) ] * dt; % k2 の計算
x_k1 = x_k + (k1 + k2) /2 ; % x_{k+1} の計算
x = [x x_k]; % xの末尾にx_kを追加
tim = [tim t]; % timの末尾にtを追加
x_k = x_k1; % 次回の計算のため x_{k} ← x_{k+1}に繰り上げ
endfor
grid on;
axis([0 10 -1 1]);
plot(tim, x(1, :)); % x_1だけをグラフ表示
xlabel('time (s)');
ylabel('theta (rad)');
4次のルンゲ・クッタ法
2次のルンゲ・クッタ法では2点における導関数の値を使いました．同様に3点，4点と点数を増やしていくことで原理的には精度がよくなります．しかし実用上，4点を利用した4次のルンゲ・クッタ法で十分なようです．そのため，単にルンゲ・クッタ法というと4次の場合を指すことが多いです．
4次のルンゲ・クッタ法の計算手順は以下のようになります．
- math(\begin{array}{l l l}\boldsymbol{k}_1 &=& \boldsymbol{f}(\boldsymbol{x}_k, t_k) \Delta t \\\boldsymbol{k}_2 &=& \boldsymbol{f}(\boldsymbol{x}_k + \frac{\boldsymbol{k}_1}{2}, t_k + \frac{\Delta t}{2}) \Delta t \\\boldsymbol{k}_3 &=& \boldsymbol{f}(\boldsymbol{x}_k + \frac{\boldsymbol{k}_2}{2}, t_k + \frac{\Delta t}{2}) \Delta t \\\boldsymbol{k}_4 &=& \boldsymbol{f}(\boldsymbol{x}_k + \boldsymbol{k}_3, t_{k+1}) \Delta t \\\boldsymbol{x}_{k+1} &=& \boldsymbol{x}_k + \frac{1}{6} ( \boldsymbol{k}_1 + 2 \boldsymbol{k}_2 + 2 \boldsymbol{k}_3 + \boldsymbol{k}_4)\end{array})
先の単振り子のシミュレーションをルンゲ・クッタ法で計算するOctaveスクリプトを以下に示します．一つ目は上記の計算手順をストレートにベタ打ちしたもの．二つ目はforループと条件判断を利用して書いたものです．
- ルンゲ・クッタ法による単振り子のシミュレーション（ベタ打ち）
dt = 0.01; % サンプリングタイム
m = 1; % 質量 1 kg
l = 0.5; % 振り子長さ 0.5 m
c = 0.5; % 空気抵抗係数
x_k = [30/180*pi ;
0 ]; % 初期角度30度, 初期速度0
x = []; % シミュレーション結果格納用変数
tim =[]; % シミュレーション時刻格納用変数
for t = 0 : dt : 10 % 時刻tは0から10までdt刻み
k1 = [ x_k(2) ;
-9.8/l * sin(x_k(1)) - c/m * x_k(2) ] * dt; % k1 の計算
x_tmp = x_k + k1/2;
k2 = [ x_tmp(2) ;
-9.8/l * sin(x_tmp(1)) - c/m * x_tmp(2) ] * dt; % k2 の計算
x_tmp = x_k + k2/2;
k3 = [ x_tmp(2) ;
-9.8/l * sin(x_tmp(1)) - c/m * x_tmp(2) ] * dt; % k3 の計算
x_tmp = x_k + k3;
k4 = [ x_tmp(2) ;
-9.8/l * sin(x_tmp(1)) - c/m * x_tmp(2) ] * dt; % k4 の計算
x_k1 = x_k + (k1 + 2 * k2 + 2 * k3 + k4 ) /6 ; % x_{k+1} の計算
x = [x x_k]; % xの末尾にx_kを追加
tim = [tim t]; % timの末尾にtを追加
x_k = x_k1; % 次回の計算のため x_{k} ← x_{k+1}に繰り上げ
endfor
grid on;
axis([0 10 -1 1]);
plot(tim, x(1, :)); % x_1だけをグラフ表示
xlabel('time (s)');
ylabel('theta (rad)');
- ルンゲ・クッタ法による単振り子のシミュレーション（ループ利用）
dt = 0.01; % サンプリングタイム
m = 1; % 質量 1 kg
l = 0.5; % 振り子長さ 0.5 m
c = 0.5; % 空気抵抗係数
x_k = [30/180*pi ;
0 ]; % 初期角度30度, 初期速度0
x = []; % シミュレーション結果格納用変数
tim =[]; % シミュレーション時刻格納用変数
for t = 0 : dt : 10 % 時刻tは0から10までdt刻み
x_k1 = x_k;
x_tmp = x_k;
for j =1 : 1 : 4
k_j = [ x_tmp(2) ;
-9.8/l * sin(x_tmp(1)) - c/m * x_tmp(2) ] * dt; % k_jの計算
if j == 3
x_tmp =x_k + k_j ;
else
x_tmp = x_k + k_j /2;
endif
if (j == 1) || (j == 4)
x_k1 = x_k1 + k_j /6;
else
x_k1 = x_k1 + k_j /3;
endif
endfor
x = [x x_k]; % xの末尾にx_kを追加
tim = [tim t]; % timの末尾にtを追加
x_k = x_k1; % 次回の計算のため x_{k} ← x_{k+1}に繰り上げ
endfor
grid on;
axis([0 10 -1 1]);
plot(tim, x(1, :)); % x_1だけをグラフ表示
xlabel('time (s)');
ylabel('theta (rad)');
線形システムのシミュレーション
システムの状態方程式
前に見たように，2階以上の常微分方程式は連立多元1階常微分方程式に変形することでオイラー法やルンゲ・クッタ法を適用することができます．ここでは，より一般的に動的システムを記述することを考えます．
一般にシステム制御で取り扱う対象は，ある変数を入力，別のある変数を出力とした入出力システムと考えます． これまでの例では明示的に入出力を考えていませんでしたが，入力を&math(\boldsymbol{u}(t));，出力を&math(\boldsymbol{y}(t));としてシステムを表現すると一般的に次のように書くことができます．
- math(\begin{array}{l l l}\frac{d \boldsymbol{x}(t)}{dt} &=& \boldsymbol{f}(\boldsymbol{x}(t), \boldsymbol{u}(t), t) \\\boldsymbol{y}(t) &=& \boldsymbol{g}(\boldsymbol{x}(t), \boldsymbol{u}(t), t)\end{array})
一つ目の式は，先に述べた連立多元1階常微分方程式において入力&math(\boldsymbol{u}(t));を明示的に表記したものです．二つ目は，出力方程式と呼ばれるもので，関数&math(\boldsymbol{g}());を出力関数と呼びます．
ここで，システムのパラメータ（係数）が時間によって変化せず（時不変），関数&math(\boldsymbol{f}());と&math(\boldsymbol{g}());が状態変数と入力変数の一次結合（定数倍して加算）だけで表される場合，そのシステムは線形時不変であるといいます．このとき，システムを表す方程式は以下のように行列を使って簡単に表されます．
- math(\begin{array}{l l l}\frac{d \boldsymbol{x}(t)}{dt} &=& \boldsymbol{A} \boldsymbol{x}(t) + \boldsymbol{B} \boldsymbol{u}(t) \\ \boldsymbol{y}(t) &=& \boldsymbol{C} \boldsymbol{x}(t) + \boldsymbol{D} \boldsymbol{u}(t)\end{array})
簡単な例として，粘性係数&math(D);のダンパーと弾性係数&math(K);のバネとで壁面に接続されている質量&math(M);の物体を考え，入力を物体への力&math(u(t));，出力を物体の変位&math(y(t));としたシステムの状態方程式を求めてみます．
- 慣性・粘性・弾性系
- ref(MDKsystem.png)
物体の運動方程式は，
- math(M \frac{d^2 y(t)}{dt^2} = - D \frac{dy(t)}{dt} - K y(t) + u(t))
となります．これを状態変数&math(x_1(t));, &math(x_2(t));を
- math(\begin{array}{l l l}x_1(t) &=& y(t) \\x_2(t) &=& \frac{d y(t)}{dt}\end{array})
とおくことで以下のような連立微分方程式に書き換えられます．
- math(\begin{array}{l l l}\frac{d x_1(t)}{dt} &=& x_2(t) \\\frac{d x_2(t)}{dt} &=& - \frac{K}{M} x_1(t) - \frac{D}{M} x_2(t) + \frac{1}{M} u(t)\end{array})
これを行列を用いて表現すると以下のようになります．
- math(\frac{d}{dt} \left(\begin{array}{c}x_1(t) \\x_2(t)\end{array}\right) =\left(\begin{array}{cc}0 & 1 \\-\frac{K}{M} & -\frac{D}{M}\end{array}\right)\left(\begin{array}{c}x_1(t) \\x_2(t)\end{array}\right)+ \left(\begin{array}{c}0 \\\frac{1}{M}\end{array}\right)u(t))
状態変数を組み合わせたものを状態ベクトル&math(\boldsymbol{x}(t));とすれば，
- math(\frac{d\boldsymbol{x}(t)}{dt} =\left(\begin{array}{cc}0 & 1 \\-\frac{K}{M} & -\frac{D}{M}\end{array}\right)\boldsymbol{x}(t)+ \left(\begin{array}{c}0 \\\frac{1}{M}\end{array}\right)u(t))
となり，これが状態方程式になります．出力は，状態量&math(x_1(t));そのものですから，
- math(y(t) = \left(\begin{array}{cc}1 & 0\end{array}\right)\boldsymbol{x}(t))
となります．
以上から，状態空間表現の行列&math(\boldsymbol{A},\boldsymbol{B},\boldsymbol{C},\boldsymbol{D});は以下のようになります．
- math(\boldsymbol{A} = \left(\begin{array}{cc}0 & 1 \\-\frac{K}{M} & -\frac{D}{M}\end{array}\right), \,\boldsymbol{B} = \left(\begin{array}{c}0 \\\frac{1}{M}\end{array}\right), \, \boldsymbol{C} = \left(\begin{array}{cc}1 & 0\end{array}\right), \, \boldsymbol{D} = 0)
オイラー法によるシミュレーション
状態空間表現された線形システムをオイラー法でシミュレーションすることを考えます．出力方程式はダイナミクスを持たない（微分方程式ではない）ので，オイラー法を適用するのは状態方程式のほうだけです．
状態方程式に対するオイラー法の演算は以下のように書き表せます．
- math(\begin{array}{l l l}\boldsymbol{x}_{k+1} &=& \boldsymbol{x}_k + ( \boldsymbol{A} \boldsymbol{x}_k + \boldsymbol{B} \boldsymbol{u}_k ) \Delta t \\\boldsymbol{y}_{k} &=& \boldsymbol{C} \boldsymbol{x}_k + \boldsymbol{D} \boldsymbol{u}_k \end{array})
この式に基づいて先の慣性・粘性・弾性系のシミュレーションを行うOctaveスクリプトの例を示します．入力&math(\boldsymbol{u}(t));を正弦波とした強制振動を想定しています．
- オイラー法による慣性・粘性・弾性系のシミュレーション
dt = 0.01; % サンプリングタイム
M = 1; % 質量 1 kg
D = 2; % 粘性係数 2 Ns/m
K = 1; % 弾性係数 1 N/m
x_k = [0 ;
0 ]; % 初期位置0, 初期速度0
x = []; % シミュレーション結果格納用変数（状態量）
u =[]; % シミュレーション結果格納用変数（入力）
y =[]; % シミュレーション結果格納用変数（出力）
tim =[]; % シミュレーション結果格納用変数（時刻）
A = [ 0, 1 ;
-K/M, -D/M];
B = [ 0 ;
1/M];
C = [1, 0];
D = 0;
for t = 0 : dt : 10 % 時刻tは0から10までdt刻み
u_k = 10 * sin(2*pi*1*t); % u_kの計算．振幅10 N，周波数1 Hzの正弦波
x_k1 = x_k + (A * x_k + B * u_k) * dt; % x_{k+1} の計算
y_k = C * x_k + D * u_k; % y_k の計算
x = [x x_k]; % xの末尾にx_kを追加
u = [u u_k]; % uの末尾にu_kを追加
y = [y y_k]; % yの末尾にy_kを追加
tim = [tim t]; % timの末尾にtを追加
x_k = x_k1; % 次回の計算のため x_{k} ← x_{k+1}に繰り上げ
endfor
grid on;
axis([0 10 -1 1]);
plot(tim, y); % yをグラフ表示
xlabel('time (s)');
ylabel('y');
ルンゲ・クッタ法によるシミュレーション
状態空間表現された線形システムを(4次の)ルンゲ・クッタ法でシミュレーションする演算式は以下のようになります．
- math(\begin{array}{l l l}\boldsymbol{k}_1 &=& \left( \boldsymbol{A} \boldsymbol{x}_k + \boldsymbol{B} \boldsymbol{u}_k \right) \Delta t \\\boldsymbol{k}_2 &=& \left( \boldsymbol{A} (\boldsymbol{x}_k + \frac{\boldsymbol{k}_1}{2}) + \boldsymbol{B} \boldsymbol{u}(t_k + \frac{\Delta t}{2}) \right) \Delta t \\\boldsymbol{k}_3 &=& \left( \boldsymbol{A} (\boldsymbol{x}_k + \frac{\boldsymbol{k}_2}{2}) + \boldsymbol{B} \boldsymbol{u}(t_k + \frac{\Delta t}{2}) \right) \Delta t \\\boldsymbol{k}_4 &=& \left( \boldsymbol{A} (\boldsymbol{x}_k + \boldsymbol{k}_3) + \boldsymbol{B} \boldsymbol{u}_{k+1} \right) \Delta t \\\boldsymbol{x}_{k+1} &=& \boldsymbol{x}_k + \frac{1}{6} ( \boldsymbol{k}_1 + 2 \boldsymbol{k}_2 + 2 \boldsymbol{k}_3 + \boldsymbol{k}_4)\\\boldsymbol{y}_{k} &=& \boldsymbol{C} \boldsymbol{x}_k + \boldsymbol{D} \boldsymbol{u}_k \end{array})
ルンゲ・クッタ法の場合，時刻&math(t_k);だけでなく，&math(t_k + \frac{\Delta t}{2});や&math(t_{k+1});での入力&math(\boldsymbol{u}(t));が必要になることに注意します．
ルンゲ・クッタ法で慣性・粘性・弾性系のシミュレーションを行うOctaveスクリプトの例は以下の通りです．
- ルンゲ・クッタ法による慣性・粘性・弾性系のシミュレーション
dt = 0.01; % サンプリングタイム
M = 1; % 質量 1 kg
D = 2; % 粘性係数 2 Ns/m
K = 1; % 弾性係数 1 N/m
x_k = [0 ;
0 ]; % 初期位置0, 初期速度0
x = []; % シミュレーション結果格納用変数（状態量）
u =[]; % シミュレーション結果格納用変数（入力）
y =[]; % シミュレーション結果格納用変数（出力）
tim =[]; % シミュレーション結果格納用変数（時刻）
A = [ 0, 1 ;
-K/M, -D/M];
B = [ 0 ;
1/M];
C = [1, 0];
D = 0;
for t = 0 : dt : 10 % 時刻tは0から10までdt刻み
x_k1 = x_k;
x_tmp = x_k;
u_k = 10 * sin(2*pi*1*t); % u_kの計算．振幅10 N，周波数1 Hzの正弦波
u_k1 = 10 * sin(2*pi*1*(t+dt)); % u_{k+1}の計算．
u_k12 = 10 * sin(2*pi*1*(t+dt/2)); % u_{k+1/2}の計算．
u_tmp = u_k;
for j =1 : 1 : 4
k_j = ( A * x_tmp + B * u_k ) * dt; % k_j の計算
if j == 3
x_tmp = x_k + k_j ;
u_tmp = u_k1;
else
x_tmp = x_k + k_j /2;
u_tmp = u_k12;
endif
if (j == 1) || (j == 4)
x_k1 = x_k1 + k_j /6;
else
x_k1 = x_k1 + k_j /3;
endif
endfor
y_k = C * x_k + D * u_k; % y_k の計算
x = [x x_k]; % xの末尾にx_kを追加
u = [u u_k]; % uの末尾にu_kを追加
y = [y y_k]; % yの末尾にy_kを追加
x_k = x_k1; % 次回の計算のため x_{k} ← x_{k+1}に繰り上げ
endfor
grid on;
axis([0 10 -1 1]);
plot(tim, y); % yをグラフ表示
xlabel('time (s)');
ylabel('y ');
Octave関数を利用したシミュレーション
多段解法
これまで微分方程式の数値解法としてオイラー法とルンゲ・クッタ法について述べましたが，両手法とも&math(x_{k+1});を求めるために状態量としては&math(x_k);だけを利用する方法でした（一段階法）．しかし，&math(x_{k-1});や&math(x_{k-2});といった，より過去の情報も利用する多段階法という手法があります．複数の過去の情報を利用することで計算精度を上げることができるだけでなく，数値解法の安定性の点でも優れることがわかっています．すなわち，オイラー法やルンゲ・クッタ法では，サンプリング周期&math(\Delta t);を大きくとりすぎるとシミュレーション結果が発散してしまうことがありましたが，多段階法では発散しにくいのです．この性質は，時間のかかるシミュレーションをする前に&math(\Delta t);を大きめにとってラフなシミュレーションをすることができるので便利です．
lsode関数の使用方法
一般に多段階法の計算アルゴリズムは複雑なので，1からプログラムを書くのはたいへんです．幸い，Octaveには微分方程式を数値的に解くためにlsodeという関数が用意されていて，この関数は代表的な多段解法である後退微分法（BDF法）やアダムス法を利用しています（オプションで指定できます．デフォルトはBDF法）．
以下のような常微分方程式の数値解をlsode関数で求めるための基本的な手順は次のようになります．
- math(\frac{d \boldsymbol{x}(t)}{dt} = \boldsymbol{f}(\boldsymbol{x}(t), t) \quad , \ \boldsymbol{x}(0) = \boldsymbol{x}_0)
- 1. まず，微分方程式を適当な関数名（例えばfnc）でOctaveのユーザ関数として定義します．もし解きたい微分方程式が2階以上なら，関数の引数xと戻り値xdotはそれぞれベクトルになります．
function xdot = fnc(x, t)
xdot = ・・・
endfunction
- 2. 初期値を例えばx_0をいう変数名でセットします．
- 3. 数値解を求めたい時間範囲を，linspace関数などを利用してベクトル変数（例えば変数名t）として用意します．
- 4. 以下のようにlsode関数を呼び出すと，変数xに数値解が代入されます．このとき，xは列の数が状態変数の個数（次元），行の数が時間変数tの点数になります．
x = lsode (fcn, x_0, t);
先の単振り子の自由応答のシミュレーションをlsode関数を使って求めるスクリプト例を示します．
Octaveでユーザ関数を定義するには，functionとendfunctionとで定義式を囲み，関数名と同じ名前で拡張子が.mのファイル（この例だとfnc.m）として保存するのが基本です．そして，そのファイルをカレントディレクトリまたはpathの通ったディレクトリに置いておけば，コマンドラインやその他のスクリプトファイルから呼び出すことができます．しかし，シミュレーション用のスクリプトと別ファイルでユーザ関数を定義することになり，管理が面倒になる可能性があります．
そこで，1つのファイルにユーザ関数定義とシミュレーションのスクリプトを合わせることを考えます．ここで注意が必要です．Octaveはユーザ関数の定義はそれが呼び出される前にしておかなければなりません．よってファイルを１つにまとめる場合，functionとendfunctionで囲まれた部分はlsode関数より前になければなりません．しかし，ファイルの先頭に書いてしまうと，キーワードfunctionで始まるファイルであるということでOctaveがおせっかいにもファイル名が関数名と合わないとか言ってきてうまく実行できません．これを回避するには，関数の定義を先頭ではなくて，かつ，lsode関数より前の適当なところに書くか，ファイルの先頭に例えば1;といった無意味なコマンドを挿入し，その直後に関数定義をする方法をとります．ここで示したスクリプト例では前者の方法をとっています．
- lsode関数を使った単振り子のシミュレーション1
dt = 0.01; % サンプリングタイム
x_0 = [30/180*pi ;
0 ]; % 初期角度30度, 初期速度0
t = linspace(0,10,10/dt+1); % 0秒と10秒の両端が入るので点数は1つ余計になることに注意
function xdot = fnc(x, t)
m = 1; % 質量 1 kg
l = 0.5; % 振り子長さ 0.5 m
c = 0.5; % 空気抵抗係数
xdot = [ x(2) ;
-9.8/l * sin(x(1)) - c/m * x(2) ] ;
endfunction
x = lsode("fnc", x_0, t);
grid on;
axis([0 10 -1 1]);
plot(t, x(:, 1)); % thetaすなわちx(1)をグラフ表示．xは列方向が時間軸になっているの注意
xlabel('time (s)');
ylabel('theta ');
なお，このスクリプト例では関数定義の中で質量&math(m);などの定数を定義しています．これらの定数は関数の外側で定義できた方が便利な場合が多いです．通常，関数内部で宣言した変数はローカル変数となり，外部とは名前が同じでも内容は異なるということになってしまうので，グローバルな変数として宣言する必要があります．それにはglobal命令を使います．その際に，関数内部と外部の両方で宣言する必要があるので注意が必要です．
- lsode関数を使った単振り子のシミュレーション2
global m l c; % 関数内部でも使うのでグローバル変数として宣言
dt = 0.01; % サンプリングタイム
m = 1; % 質量 1 kg
l = 0.5; % 振り子長さ 0.5 m
c = 0.5; % 空気抵抗係数
x_0 = [30/180*pi ;
0 ]; % 初期角度30度, 初期速度0
t = linspace(0,10,10/dt+1); % 0秒と10秒の両端が入るので点数は1つ余計になることに注意
function xdot = fnc(x, t)
global m l c; % ここでも再びグローバル変数として宣言
xdot = [ x(2) ;
-9.8/l * sin(x(1)) - c/m * x(2) ] ;
endfunction
x = lsode("fnc", x_0, t);
grid on;
axis([0 10 -1 1]);
plot(t, x(:, 1)); % thetaすなわちx(1)をグラフ表示．xは列方向が時間軸になっているの注意
xlabel('time (s)');
ylabel('theta ');
lsode関数による線形システムのシミュレーション
次に，状態空間表現された線形システムのシミュレーションをlsode関数を用いて行うスクリプト例を示します．システムが，次のような状態方程式と出力方程式で表されているとします．
- math(\begin{array}{l l l}\frac{d \boldsymbol{x}(t)}{dt} &=& \boldsymbol{A} \boldsymbol{x}(t) + \boldsymbol{B} \boldsymbol{u}(t) \\\boldsymbol{y}(t) &=& \boldsymbol{C} \boldsymbol{x}(t) + \boldsymbol{D} \boldsymbol{u}(t)\end{array})
lsode関数で計算するのは状態方程式のほうだけです．その結果得られた状態変数ベクトル&math(\boldsymbol{x});を用いて出力方程式に基づいて出力ベクトル&math(\boldsymbol{y});を計算します．
ここでは前述の慣性・粘性・弾性系の強制振動をlsode関数を用いてシミュレーションするスクリプト例を示します．このスクリプトで注意すべきなのは，入力&math(u(t));の取り扱いです．&math(u(t));はlsodeの内部では変数tの代数式として記述しますが，lsodeで得られた状態ベクトルを用いて出力yを計算する際には時間軸ベクトルtに対応したベクトルuとして必要になります．Octaveではsin関数などが引数にベクトルを受けると，要素ごとに関数をかけて同じサイズのベクトルを返しますのでそれを利用しています．&math(u(t));の内容によっては，要素ごとの積算 .* や要素ごとの除算 ./ といった演算子が必要になるので注意です．
- lsode関数を使った慣性・粘性・弾性系のシミュレーション
global A B; % 関数内部でも使うのでグローバル変数として宣言
dt = 0.01; % サンプリングタイム
x_0 = [0;
0 ]; % 初期位置0, 初期速度0
M = 1; % 質量 1 kg
D = 2; % 粘性係数 2 Ns/m
K = 1; % 弾性係数 1 N/m
A = [ 0 , 1 ;
-K/M, -D/M];
B = [ 0 ;
1/M];
C = [1, 0];
D = 0;
t = linspace(0,10,10/dt+1);
function xdot = fnc(x, t)
global A B; % ここでも再びグローバル変数として宣言
u = 10 * sin(2*pi*1*t);
xdot = A*x + B*u;
endfunction
x = lsode("fnc", x_0, t);
u = 10 * sin(2*pi*1*t); % 時間ベクトルtに対応した入力ベクトルuの計算
y = C * x '+ D * u; %出力方程式に基づいてyの計算
grid on;
axis([0 10 -1 1]);
plot(t, y); % yをグラフ表示
xlabel('time (s)');
ylabel('y ');
線形離散時間システムとしてのシミュレーション
これまで常微分方程式の数値解法に基づいて動的システムのシミュレーション方法を考えてきましたが，シミュレーション対象が線形システムで，入力がサンプリング期間中一定値と考えてよい場合，線形離散時間システムとして扱うことで簡単にシミュレーションを行うことができます．
前述の状態空間表現された線形連続時間システムの入力側にサンプラと零次ホールドをつけ，出力をサンプルしたシステムは，次のように離散時間システムとして表すことができます．
- math(\begin{array}{l l l}\boldsymbol{x}_{k+1} &=& \boldsymbol{\Phi} \boldsymbol{x}_k + \boldsymbol{\Gamma} \boldsymbol{u}_k \\\boldsymbol{y}_k &=& \boldsymbol{C} \boldsymbol{x}_k + \boldsymbol{D} \boldsymbol{u}_k\end{array})
ここで，
- math(\begin{array}{l l l}\boldsymbol{\Phi} &=& e^{\boldsymbol{A} \Delta t} \\\boldsymbol{\Gamma} &=& \int_0 ^{\Delta t} e^{\boldsymbol{A} \tau}\boldsymbol{B} d\tau\end{array})
です．
ここで，&math(\boldsymbol{\Phi});と&math(\boldsymbol{\Gamma});が求まれば，離散時間システムの式を直接演算するだけで行えるので簡単です．
連続時間システムから離散時間システムへの変換には，Octaveに用意されているc2d関数が利用できます．この関数はサンプル・ホールド付きのz変換か双一次変換のいずれかで離散化することができます（デフォルトはサンプル・ホールド）．シミュレーションにはサンプル・ホールド付きのz変換を使います．
次のスクリプトは，先の慣性・粘性・弾性系の強制振動を離散時間システムに変換した上でシミュレーションするものです．まず，ss関数で4つの行列による状態空間表現をOctave上のシステム行列表現に変換し，その上でc2d関数で離散化しています．シミュレーションの計算自体は離散時間システムの式を繰り返し演算しているだけです．
- 離散時間システムに変換した慣性・粘性・弾性系のシミュレーション
dt = 0.01; % サンプリングタイム
M = 1; % 質量 1 kg
D = 2; % 粘性係数 2 Ns/m
K = 1; % 弾性係数 1 N/m
x_k = [0 ;
0 ]; % 初期位置0, 初期速度0
x = []; % シミュレーション結果格納用変数（状態量）
u =[]; % シミュレーション結果格納用変数（入力）
y =[]; % シミュレーション結果格納用変数（出力）
tim =[]; % シミュレーション結果格納用変数（時刻）
A = [ 0, 1 ;
-K/M, -D/M];
B = [ 0 ;
1/M];
C = [1, 0];
D = 0;
sys = ss(A, B, C, D); % systemマトリックスに変換
sys_d = c2d(sys, dt); % 離散時間システムに変換（サンプルホールド付きz変換）
[A_d, B_d, C_d, D_d] = sys2ss(sys_d); % 離散時間システムの各行列取り出し
for t = 0 : dt : 10 % 時刻tは0から10までdt刻み
u_k = 10 * sin(2*pi*1*t); % u_kの計算．振幅10 N，周波数1 Hzの正弦波
x_k1 = A_d * x_k + B_d * u_k; % x_{k+1} の計算
y_k = C_d * x_k + D_d * u_k; % y_k の計算
x = [x x_k]; % xの末尾にx_kを追加
u = [u u_k]; % uの末尾にu_kを追加
y = [y y_k]; % yの末尾にy_kを追加
tim = [tim t]; % timの末尾にtを追加
x_k = x_k1; % 次回の計算のため x_{k} ← x_{k+1}に繰り上げ
endfor
grid on;
axis([0 10 -1 1]);
plot(tim, y); % yをグラフ表示
xlabel('time (s)');
ylabel('y');