技 術 論 文

フリーソフトウェアScilab/Scicosによる数値計算
Numerical Computations with Scilab/Scicos
矢田部 学＊
Manabu Yatabe

数理モデルを開発する場合、その妥当性を数値計算により確認することが重要である。この際、数値計算に特化したツールを用いると効率的である。モデル検証のために、我々のグループは、多数の数学関数を備えているフリーソフトウェアのScliab/Scicosを使用している。本稿では、Scilab/Scicosによる２つの基本的な数値計算例について紹介する。
In developing mathematical models, it is essential to confirm the validity of the formulations bynumerical computations. Use of specialized software is effective for this purpose. To verify themathematical models, members of our group have employed Scilab/Scicos, a free scientificsoftware package having numerous mathematical functions. In this report, the author shows twobasic numerical computations using Scilab/Scicos.

1．まえがき令や関数が豊富に用意されており、比較的簡単に数値計算プログラムを作成することが可能である。

鎌倉事業部第一技術部では防衛関係のM&S（Model-

また、ScicosはScilabに備えられているGUIを有する

ing and Simulation）業務を行っている。M&S業務におtoolboxで、この機能を用いるとブロック線図を用いていては、基礎原理から新たにモデルを定式化したり、既直截的にモデルを作りシミュレーションを行うことがで存のモデルを考え直す作業が発生する。その場合、モデきる。これはMATLABのSimulink※1に類似している。

ルの妥当性を数値計算により検証することが重要であ

Scilab/ScicosはWindows、LinuxおよびMacOSに対応る。この目的のためには、豊富な数学関数やtoolboxをし、インターネット上からダウンロードして自由に使用備えたMATLAB※1などの数値計算に特化したソフトウすることできる（http://www.scilab.org/）。

ェアの利用が便利である。
しかし、複数の計算機を用いてグループで作業を行う場合、ライセンスの問題が発生する。このような場合に使用方法についてはWeb検索、一般の書籍（例えば⑴⑵⑶⑷）、およびScilabのHELP機能などから必要な情報を入手できる。

は、フリーのソフトウェアの使用が考えられる。そのた

なお、2009年12月リリースのScilab5.2からScicosは、め、我々のグループではモデルの検証ツールとしてフリそれを発展させたXcosに変更されている（http://www.

ーソフトウェアのScilab/Scicosを用いている。

scilab.org/products/xcos）。ただし、2010年７月時点で本稿では、数理モデルの検証に有用なScilab/Scicosによる数値計算について紹介する。

は、古いバージョンのScilab/Scicosもダウンロード可能である。

2．Scilab/Scicosとは3．数値計算例

Scilabは、フランスのINRIA（Institut national deScilabによる数値計算のプログラムとScicosのプロッ

recherche en informatique et en automatique、国立情ク線図によるシミュレーションの例を紹介する。ここで報学自動制御研究所）とENPC（École nationale desは、2007年10月 に リ リ ー ス さ れ たScilab 4.1.2 お よ び

ponts et chaussées、国立土木学校）により開発された

Scicos 4.2 を用いた。

数値計算用のソフトウェアである。これは、線形代数、※１

制御系、信号処理、グラフィック機能などに対応する命49＊鎌倉事業部 第一技術部

MATLAB、Simulinkは、米国The MathWorks社の登録商標である。

3.1

式⑴は、

Scilabによるプログラミング

図１のような弾性定数 のバネ、減衰定数 のダンパー、および質量とおくことにより、１階の連立微分方程式の質点よりなる系を考える。この系と速度に比例する力には、変位に比例する力が作用するので、運動方程式はと表される。ただし、⑵⑴
とおおよびとなり、Runge−Kutta法が適用できる。この解法に従ったScilabのプログラムを図２に示す。ここで、Runge−Kutta法のメイン処理が関数rk4であり、微分方程式いた。また、ドットは時間微分を表す。
この方程式の解を数値計算により求める。Scilabには、の 定 義 が 関 数diffEqsで あ る。 こ の 例 に 示 す よ う に、常微分方程式を解く関数odeが用意されているが、プロScilabはベクトルや行列を直接取り扱うことができるた

グラム例を示すという目的から、ここでは敢えて４次のめに、成分ごとの処理をする必要がなくプログラムが簡Runge−Kutta法⑷⑸によるプログラムを作成する。

潔になる。ただし、Scilabでは、関数rk4で使用するユーザー定義の関数diffEqsを命令getfにより、使用前にロードする必要がある（図２参照）。これはMATLABなどとは異なる。
このプログラムを実行した結果を図３に示す。これより、変位（赤）と速度（青）が時間とともに減衰して行

図１

く様子が分かる。なお、ここでは次のパラメータを設定減衰振動モデルした。

関数rk4（Runge-Kutta法の処理）関数diﬀEqs（微分方程式の定義）

function [ ] = rk4( )function [xdot] = diﬀEqs(t,x);////// 減衰振動の微分方程式

４次のRunge-Kutta法で微分方程式dx/dt=f(x,t)を初期条件(x0,t0)で解く

getf ('diﬀEqs.sci');

微分方程式を定義した関数のロード

A = [ 0 1; -omg2−2*gmm];
xdot = A*x;
endfunction

global omg2 gmm;
omg2 = 2.0; // ばね定数/質量gmm = 0.2; // 減衰定数/質量

微分方程式

x = [0; sqrt(omg2)]; // 初期条件
tmax = 30; // 評価終了時刻
dt = 0.1; // 積分刻み
i = 1;
T(i,1) = 0; // 初期時刻 ⇒ T(1,1)
X(i,1) = x(1); X(i,2) = x(2); // 初期値 x ⇒ X(1,:)
for t = dt: dt: tmax
Runge-Kutta法の処理

i = i+1;
T(i,1) = t; X(i,1) = x(1); X(i,2) = x(2); // プロット用データ
end

Runge-Kutta法⇔[k1] = dt*diﬀEqs(t,x);[k2] = dt*diﬀEqs(t+dt/2,x+k1/2);[k3] = dt*diﬀEqs(t+dt/2,x+k2/2);[k4] = dt*diﬀEqs(t+dt,x+k3);x = x+(k1+2*k2+2*k3+k4)/6;// 結果のプロット
clf( );
plot2d(T,X);
xset("font size",3);
xtitle("Dampling Oscillation","Time","Displacement & Velocity");
endfunction

図２

Runge−Kutta法のプログラム例

MSS技報・Vol.2150

技 術 論 文

Damped Oscillation

システム

観測ノイズ

Displacement（m）/ Velocity
（m/s）

1.5

測定

Displacement

1.0

変位 を観測0.50.0Kalmanフィルター

システムノイズ−0.5

変位 と速度を推定−1.0

Velocity−1.5051015202530

図４

Kalmanフィルター推定のイメージ

Time
（s）

図３

プログラムの実行結果

として⑷と書ける。ただし、

はね定数/質量：減衰定数/質量：初期変位：である。いまの場合、とはスカラー量であるが、初期速度：一般性を持たせるためベクトル的にボールド体で表記した。
式⑶と式⑷で表されるシステムについてのKalmanフ
ィルター処理は以下のとおりである（例えば⑹⑺⑻⑼⑽3.2

Scicosによるシミュレーションなどを参照）。

Scicosのブロック線図を用いたシミュレーションの例
として、図１に述べた減衰振動のシステムの変位を測定して、その時点の変位と速度を推定するKalmanフィル
ターを考える（図４）。ただし、システムはモデル化で⑸考慮されていない擾乱を考慮するためにシステムノイズが加わるものとする。これは測定時の観測ノイズとは異なるものである。
図４の系にKalmanフィルターのアルゴリズムを適用するための定式化について述べる。式⑵より

ここで、式⑶と式⑷で述べた記号以外の新たなものは（Kalmanゲイン）、（推定誤差の共分散行列）、（システムノイズの共分散行列）、および（観測ノイズ の 共 分 散 行 列） で あ る。 上 付 き 添 え 字 は 転 置 行である。この解は、時刻、時間刻みをにおける状態ベクトルをとして、離散的に考え列、−１は逆行列を意味する。
また、推定値であることを強調するために状態ベクト
ルにハットをつけ、状態ベクトルと共分散行列についてるとは観測更新前の値を（−）
、観測更新後の値を（+）をつけて区別した。と書ける。ここで、以上に基づいてKalmanフィルターのシミュレーションを行う。このときのScicosのブロック線図を図５に示を遷移行列とする。ただし、 は２×２単位行列であす。このブロック線図の意味は次のとおりである。左かる。このとき、状態方程式は、擾乱（システムノイズ）

ら 右 に 眺 め る と、 先 ず、（I）シ ス テ ム の 信 号 が 生 成（System）され、（II）それを測定（Measurement）して、を考慮して⑶と表される。
他方、測定については変位しか観測できないとすると、 時 刻51に お け る 観 測 方 程 式 は、 観 測 ノ イ ズ を（III）
Kalmanフィルター処理（Kalman Filter）を行う。
ここで、図５の各処理ブロック（System、Measurement、Kalman Filter）はユーザー定義のブロックであり、図６に示すブロック線図をSuper Block機能を用いKalman Filter Simulation

randomgenerator

NoiseSystem

Measurement

図５

write tooutput ﬁleKalman Filter

ScicosによるKalmanフィルターシミュレーション

1

SystemMeasurement

randomgeneratorFH

MATMUL

1

MATMUL

1

1

1/z

1

Kalman Filter
1

1

Scifunc
Mux
Measurement Update

2
Mux

1/z
1/z
Scifunc
Time Update

1/z
1/z
1/z
1/z

図６

Super Blockの中で定義したブロック線図

MSS技報・Vol.2152

技 術 論 文てグループ化したものである。図６において、Systemは減衰振動の力学系にシステムノイズが加わった式⑶の処理である。Measurementは式⑷で表される測定装置を模擬する。Kalman Filterは式⑸のKalmanフィルター処理のブロックである。
更に、Kalman Filterでは観測更新（MeasurementUpdate）と時間更新（Time Update）に対してユーザーが任意に処理を定義できる関数Scifuncを用いている。この処理はScilabのスクリプト言語で記述される（図７参照）。これらは式⑸に示した処理である。
図５に示したブロック線図によるシミュレーション結果を図８に示す。図の上段は変位の観測値であり、下段はこの観測値に基づいてKalmanフィルターで推定した

システム
ばね定数/質量減衰定数/質量初期変位
初期速度
システムノイズ
測定
観測間隔
観測ノイズ
Kalmanフィルター
状態ベクトルの初期値
誤差共分散行列の初期値
システムノイズの共分散行列
観測ノイズの共分散行列

変位（赤）と速度（青）である。ここでは、比較のため

4．むすびに減衰振動の理論値（システムノイズ無し、黒破線）を重ね描きしている。変位と速度の真の初期値（0，1.41）

本稿では、M&S業務で数理モデルの検証に利用してから大きく離れたアプオリな初期値（5，5）から推定をいるフリーソフトウェアのScilab/Scicosの数値計算例に開始したにもかかわらず、Kalmanフィルターが変位とついて紹介した。

速度を推定していることが分かる。
なお、このシミュレーションで設定したパラメータはMeasurement

以下のとおりである。

1.5

time = u1; Y = u2;
Xm(1) = u3; Xm(2) = u4;
Pm(1,1) = u5; Pm(1,2) = u6; Pm(2,1) = u7; Pm(2,2) = u8;
if modulo(time,dt̲obs) <= 1e-5 then
K = Pm*H'*inv(H*Pm*H'+R);
Xp = Xm+K*(Y-H*Xm);
Pp = (eye(2,2)-K*H)*Pm;
else
Xp = Xm;
Pp = Pm;end
y1 = Xp(1); y2 = Xp(2);y3 = Pp(1,1); y4 = Pp(1,2); y5 = Pp(2,1); y6 = Pp(2,2);

Scifunc (Time Update)Xp(1) = u1; Xp(2) = u2;
Pp(1,1) = u3; Pp(1,2) = u4; Pp(2,1) = u5; Pp(2,2) = u6;
Xm = F*Xp;
Pm = F*Pp*F'+Q;
y1 = Xm(1); y2 = Xm(2);y3 = Pm(1,1); y4 = Pm(1,2); y5 = Pm(2,1); y6 = Pm(2,2);

0.50.0−0.5−1.0−1.50510152025302530Time（s）
Kalman Filter Estimation
8

Displacement（m）/Velocity（m/s）

Scifunc (Measurement Update)

Displacement（m）

1.06

Displacement420−2−4

Velocity−605101520Time（s）

図７53

Scifuncの記述例

図８

Kalmanフィルターシミュレーションの実行結果

これはフリーソフトウェアなので業務から離れて自己⑹啓発のためにも使える。例えば、物理や工学などの理論

1974を修得する際、身近な数値計算ツールとして利用すれば⑺理解が深まる。興味を持たれた方は、Scilab/Scicosを試⑻みられることをお勧めする。
参考文献⑴

Gelb, A.：Applied Optimal Estimation, MIT Press,
片山徹：応用カルマンフィルタ、朝倉書店、1983

Brown, R. and P. Hwang：Introduction to Random
Signals and Applied Kalman Filtering with Matlab
Exercises and Solutions, John Wiley & Sons, 1996⑼

Grewal, M. and A. Andrews：Kalman Filtering:

櫻井鉄也：MATLAB/Scilabで理解する数値計算、

Theory and Practice Using MATLAB, John Wiley

東京大学出版会、2003& Sons, 2001

⑵ Campbell, S., J. Chancelier and R. Nikoukhah：⑽

西山清：最適フィルタリング、培風館、2001

Modeling and Simulation in Scilab/Scicos,
Springer, 2006⑶

橋本洋志、石井千春：Scilab/Scicosで学ぶシミュレーションの基礎−自然・社会現象から、経済・金融、システム制御まで，オーム社、2008⑷

川田昌克：Scilabで学ぶわかりやすい数値計算法、
森北出版、2008⑸ Woan, G.：The Cambridge Handbook of Physics
Formulas, Cambridge Univ. Press, 2000

MSS技報・Vol.21

54