ページ2分割版

8 その他・断片的なこと
制御理論を学ぶ学生向けの Matlab・Scilab 入門
はじめに
1
1.1
神戸大学工学部
森 耕平
この文書のこと
この文章は，2005 年 4 月から 2005 年 5 月にかけて森耕平が書いた Matlab 入門
用の解説に加筆修正を施したものです． Matlab と Scilab の基本演算から常微分
方程式の扱いまでをさらっと述べます．
- 2007/10/16
1.2
著作権・配布条件関連情報
Matlab とは何か
Matlab(マトラブ : Matrix laboratory) 1 とは,
個人的な用途のために勝手に利用してもかまいません．修正を行っていない
ならば再配布してもかまいません (PS や PDF のみの配布も可)．加筆や修正
を行ったものであっても，著作権情報が削除されておらず，配布条件が修正
されておらず，LATEX のソースコードと共に配布され，加筆修正した旨が加
筆修正後の文書に明記されているならば再配布してかまいません．
• 機能が非常に豊富な電卓
• 数値計算に特化したプログラミング環境
• 工学・産業上の実際の設計業務や解析業務を支援するソフトウェア
として便利な, 非常に有名なソフトウェアである．情報知能工学科の演習室のもの
は Mac OS X 用であるが，もちろん，Windows 版や Linux 版も存在する．Matlab
の特徴として以下のものが挙げられる．
2005 年 11 月 8 日 (火)
森 耕平
[email protected]
Contents
1 はじめに
1.1 この文書のこと . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Matlab とは何か . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 本資料の目的・表記規則など . . . . . . . . . . . . . . . . . . . . .
2
2
2
3
2 英語
3
3 基本操作 (主に行列演算)
3.1 基本表記，および，数学・UNIX・C の知識から直観的にわかりそ
うな基本演算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 線形代数で習ったであろう単項演算・引数が一つの関数 . . . . . .
3.3 行列やベクトルに対する特殊演算 . . . . . . . . . . . . . . . . . .
3.4 オンラインマニュアル . . . . . . . . . . . . . . . . . . . . . . . .
3
4 グラフ表示
8
5 連立一次方程式を解く
9
3
6
7
8
6 常微分方程式 (線形システム) を解く
6.1 一般の微分方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2 線形常微分方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
9
10
7 少しだけ古典制御
11
1
12
• プログラミング言語や電卓としてのデータ型は行列が基本となっている．最
初から行列計算の洗練されたルーチンが組み込まれており, 行列計算を多用
するプログラムを容易に作成できる．固有値の計算などの基本的な演算は
下手な C プログラムよりも高速に行われる．
• プログラミング言語としては非常に高級な2 作りになっている．そのため，
C 言語のように何でもできるというわけではない反面，数値計算関係の複
雑な機能を素早く実装することが容易である．
• 制御シミュレーションなどのための Toolbox(関数や派生ソフトウェアを集
めたもの) が豊富に準備されている．そのため，定型的な計算の多くは自ら
プログラムを書かなくても行える．
• 非常に高価であるが，それでも企業で解析・設計業務に使用されている．
1 日本では
(株) サイバネットシステムが Matlab を販売している．詳細は以下の URL を参照の
こと．マニュアルなどのサポート情報も充実している．
http://www.cybernet.co.jp/
http://www.cybernet.co.jp/matlab/
また，Matlab とよく似た機能を持つソフトウェアとして Octave， Scilab，MaTX などが知られて
いる．これらのうち Octave 以外は GNU の定義するフリーソフトではないが，個人利用の範囲では
無料で使用でき，ソースコードの閲覧も可能である．また，Octave と Scilab の操作方法や文法は
Matlab によく似ている．
http://www.octave.org/
http://scilabsoft.inria.fr/
http://www.matx.org/
2 ハードウェアや原理から遠い．言語に組み込まれた機能が豊富．
2
1.3
本資料の目的・表記規則など
本資料では線形システム理論 (システム解析学及び演習) の初歩の段階と関連した
計算を行うために最低限知っておくべき Matlab/Scilab の操作方法を実例を中心
に説明する．Matlab や Scilab に触れる最初の数時間での到達目標は以下に述べ
る程度である．
1. Matlab/Scilab を行列計算用の電卓として楽に使えるようになる．
• 終了は「quit」．
3 行 3 列の行列 A の定義．改行が「行列としての行の変更」を表している．
>>A=[0 1 2
2 2 4
-6 2 1]
A =
2. Matlab/Scilab を用いて線形時不変システム (システム解析学及び演習の講
義で扱っているやつ) の定型的なシミュレーションができるようになる．
資料中の B=[1;0;-2] などのタイプライタ風の文字は実際の画面表示であり，
そのうち青色のものはキーボードからの入力である．赤色は実行例に追加された
補足説明である．基本的に Matlab の文法でコマンドを記述し，Scilab の文法に
よる表現は両者が違う場合のみ緑色でコマンドや補足説明を記述する．この資料
中ではコマンド入力の結果 (応答) として表示されるグラフや数値の一部が省略さ
れていることに注意して読むこと．
2
英語
Matlab/Scilab の関数やコマンドの多くは，その機能を表す英単語から連想され
るものである．それらの英単語の一部をメモしておく．また，これらの最小限の
単語を知っていれば英語のマニュアルを読むのに (それほど) 苦労しない．
(matrix:行列:「マトリックス」よりも「メトリックス」に近い) (vector:ベクト
ル) (scalar:スカラー:英語の発音は「スケーラー」に近い) (function:関数) (row:行)
(column:列) (eigenvalue:固有値) (determinant:行列式) (matrix equation:行列の
等式:行列式ではない) (origin:原点 ≡ 0 ≡ (0, 0, 0)) (nonsingular≡invertible:逆行
列を持つ) (stable:安定な) (polynomial:多項式:変数の非負の整数乗の重み付き和)
(ordinary differential equation:常微分方程式) (numerator:分子) (denominator:分
母) (rational function:有理関数:分子も分母も多項式の関数) (true:真) (false:偽)
3
3.1
基本操作 (主に行列演算)
基本表記，および，数学・UNIX・C の知識から直観的にわ
かりそうな基本演算
「>>」「-->」はプロンプトであり，Matlab がコマンド入力を受け付
けていることを表している．Matlab は C 言語と同じく大文字と小文
字を別の文字として認識する．基本演算のかなりの部分は数学とよく
似た表記が採用されているため，数式を記述するのに近い感覚で操作
ができる (ただし「=」は等号ではなく代入である)．複素数も利用で
きる．
• 行列は行と列の大きさの整合性がとれていないと足し算，引き
算，掛け算ができないことに注意．
3
0
2
-6
1
2
2
2
4
1
変数名を入力すると値が表示される
>>A
A =
0
2
-6
1
2
2
2
4
1
関数を用いて表示
>>disp(A)
0
2
-6
1
2
2
2
4
1
スカラ変数の定義
>>a=2
最後に「;」を入れると結果が表示されない (計算は行われている)．
>>b=-10;
行の変更は改行でも「;」でもよい．
>>x= [ -2 ; -3 ]
等号 (==) などの論理演算．結果の「1」「T」は真，
「0」「F」は偽．
>>x == -2
>>x == [ 0; -3]
>>a<b
>>a<=b
>>a>b
論理和．2005 年 5 月 2 日現在，演習室では「|」を直接入力できない．
「help or」
と入力すると表示される説明文中に「|」が表示されるので，それをマウスでコ
ピーすれば入力できる．
>>(a<b)|(a>b)
論理積
>>(a<b)&(a>b)
円周率 (Matlab のみ)
>>pi
4
「i」「%i」も「j」(j は Matlab のみ) も虚数単位
>>i
>>i*i
>>j*j
>>i*j
変数の上書きは自由．
>>A=[0 1 3
-1 0 1]
>>A= [ 3 3 3 ; 2 2 2 ]
定義されている変数の一覧．
>>who
変数の詳細．( Matlab のみ．Scilab では別機能．)
>>whos
ans は直前の演算結果
>>ans
虚数部を表すときに掛け算の記号「*」はなくてもよい．(Matlab のみ)
>>A=[ i 2-5i 3 ; 4 5i -2-6*j ; 1 2 3 ]
共役転置 (実数行列なら単なる転置)
>>A’
%から改行まではコメント扱いされる．剰余ではない．(//)
>>a*A % #include<stdio.h>
冪乗．(Scilab だと a^(-b) のように括弧が必要)
>>a^-b
>>A^-b
割り算
>>a/2.3
0 で割っちゃダメ
>>a/0
0 から 1 までの疑似一様乱数．(rand(1,1))
>>rand
>>rand
>>rand
指数表記．以下の例は「a=1.25*10^-2」と同じ．
>>a=1.25e-2
ファイル一覧，ディレクトリ変更，現在のディレクトリ表示，テキストファイルの
表示．当然，work というディレクトリがなければ失敗する．(pwd 以外は Matlab
のみ)
>>ls
>>cd work
>>pwd
UNIX コマンドの実行．文字列は「’」で囲む．(Scilab だと，たぶん，Scilab を
起動したターミナルに結果表示)
>>unix(’ls’)
5
3.2
線形代数で習ったであろう単項演算・引数が一つの関数
>>A = [ 1 0 0 0 ; 5 6 7 8 ; 9 10 11 12 ; 0 0 0 16 ];
行列指数関数．
expm(A) = eA =
∞
X
1 k
A
k!
k=0
ea をマクローリン展開すると以下のようになる．
ea =
∞
X
1 k
a
k!
k=0
行列指数関数 eA はこれを一般化したものである．
>>expm(A)
行列式．行列式の値が 0 ではないことと逆行列が存在することは同じ意味．
>>det(A)
要素ごとの平方根
>>sqrt(A)
行列平方根．A の行列平方根とは，X*X=A を満たす X のこと．
>>sqrtm(A)
>>sqrtm(A)*sqrtm(A)-A
逆行列
>>inv(A)
[ 2 5 9] を対角要素とする対角行列
>>diag([ 2 5 9])
行列の対角要素を並べたベクトル
>>diag(A)
固有値．A の固有値とは，0 ではない何らかのベクトル v に対して A*v=a*v を満
たすスカラー値 a のこと．(spec)
>>eig(A)
トレース．対角要素の和でもあり，固有値の和でもある．
>>trace(A)
ランク．
「行列を構成する行 (列) ベクトルのうち一次独立なものの本数」
>>rank(A)
要素が全て 0 の行列
>>zeros(2,3)
要素が全て 1 の行列
>>ones(1,2)
単位行列 (eye(3,3))
>>eye(3)
代入，足し算，引き算などの二つの行列に対する演算と関数．
>>s = 100
>>t = [ 100 ]
>>A = [ 1 2 3 4 ; 5 6 7 8 ; 9 10 11 12 ; 13 14 15 16 ]
6
>>F = [ 2 1 3 5 ; 0 2 0 -1 ; 1 0 2 0 ; 7 0 0 2 ]
右から F の逆行列を掛ける
>>A*inv(F)
>>A/F
左から F の逆行列を掛ける
>>inv(F)*A
>>F\A
\ は「help slash」と入力すると表示される．
固有値と固有ベクトル．D は A の固有値を対角要素とする対角行列．V は第 i 列が固
有値 D(i,i) に対応する固有ベクトルとなる行列．したがって，D(i,i)*V(:,i)=A*V(:,i)
が成立する．
「:」や「(i,i)」の意味については次項を参照．
>>[V,D]=eig(A)
>>D(2,2)*V(:,2)
>>A*V(:,2)
3.3
行列やベクトルに対する特殊演算
行列の要素を参照するための様々な演算が定義されている．C 言語の
「:」は行や列を範囲で
配列とは異なり，要素の添字は 1 から始まる．
指定するための記号である．
1 行 1 列から 1 行 2 列に値を代入
>>A
>>A(1,1:2) = [ 999 333]
>>u=[ 1 2 -6 ];
要素の和と積
>>sum(u)
>>prod(u)
要素の平均，最大，最小
>>mean(u)
>>max(u)
>>min(u)
行数と列数
>>size(A)
行の和と列の和
>>sum(A,1)
>>sum(A,2)
要素ごとの積，商，冪乗
>>c=[-1 1 2]; b=5;
>>c .* b
>>c ./ b
>>c .^b
要素へのアクセス
>>A(1,2)
>>A(2,3)
>>A(3,4)
2 列目
>>A(:,2)
3 列目
>>A(:,3)
2 列目から 3 列目まで
>>A(:,2:3)
全体
>>A(:,:)
t の定義. 0 から 10 までの区間を幅 0.5 の区間に分割したものがベクトル t．
>>t=[0:0.5:10]
>>u=[ 1 2 -6 ];
>>u(2:3)
>>u([1 3])=[ 2 6]
>>t(2:5)
>>A = [ 1 2 3 4 ; 5 6 7 8 ; 9 10 11 12 ]
>>B = [ -12 -11 -10 -9; -8 -7 -6 -5; -4 -3 -2 -1]
行列を用いて行列を定義
>>C = [ A(:,1) B(:,2:3) ]
>>D = [ A C ]
要素ごとの sin
>>sin(D)
3.4
オンラインマニュアル
GUI から操作するのが最もとっつきやすい．マウスで操作してもよいし，コマン
ドラインから呼び出してもよい．(Scilab だと help で GUI ヘルプ)
>>helpwin
>>helpdesk
コマンドラインから操作する場合は，例えば，以下の通り．
>>help load
次のようにすると関連コマンドを探してくれる (少し時間がかかるかもしれな
い)．(Matlab のみ)
>>lookfor expm
4
グラフ表示
グラフを表示する.
ここでの引数はベクトルである．表示されるグラフの横軸は「ベクトルの何番目
の要素か」であり，縦軸は「要素の値」である．
>>plot( [1 2 3 4 5 4 3 2 1 ] )
>>plot( sin(-1:0.2:2*pi) )
次のように入力した場合，第一引数のベクトルは横軸の座標，第二引数のベクト
ルは縦軸の座標となる．そのため，(10,1), (20,2), (30,3),...,(10,9) を結
ぶラインが描かれる．
>>plot( [10 20 30 40 50 40 30 20 10 ] , [1 2 3 4 5 6 7 8 9 ] )
7
8
重ね書きのためにグラフを保持．(Matlab のみ)
>>hold on
>>plot( [10 20 30 40 50 40 30 20 10 ] , [5 3 5 3 5 3 5 3 5 ] )
解除．
>>hold off
>>plot( cos(-1:0.2:2*pi),2*sin(-1:0.2:2*pi) )
5
>>t = 0:0.1:10;
>>y0 = [2;0];
>>t0 = 0;
>>getf(’myfunc.m’);
>>y=ode(y0,t0,t,myfunc);
>>plot2d([t t],[y(1,:) y(2,:)])
例えば，以下のような関数にすれば結果を確認しやすい．変化率が cos なので
微分方程式の解 (積分したもの) は sin になる．
連立一次方程式を解く
連立一次方程式
Ax = b
(1)
は，A が正則である場合には，
>>A= [ 1 2 ; 3 4], b= [1 ; 5]
>>x=inv(A)*b
で解ける．正則でない場合でも，連立一次方程式が解を持つ場合，つまり，
>>rank(A) - rank( [ A b ])
の結果が 0 の場合には以下のようにすれば解ける．
>>linsolve(A,b)
上記は解を一つ求めるだけ．次は全ての解を求める．
>>[x0,kerA]=linsolve(A,b)
6
6.1
常微分方程式 (線形システム) を解く
一般の微分方程式
function ydot = mycos(t,y)
ydot = cos(t)
6.2
線形常微分方程式
微分方程式の特別な場合であり，次のように書ける．
d
x = Ax
dt
(2)
外部入力信号を考えると次のようになる．
d
x = Ax + Bu
dt
(3)
y = Cx + Du
(4)
Scilab の場合，まず，上記の関数のコメント記号を Scilab 式の「//」に変更
する．また，getf() により関数を読み込む必要がある．
>>
このように表される微分方程式の扱いはさらに簡潔になる．とくに断っていなけ
れば Matlab のみです．
係数行列の定義．最後は係数行列から線形システムへの変換．
>>A=[0 1 0 ; 0 0 1 ; -1 -2 -3];
>>B=[0 0 ; -1 0; 0 2];
>>C=[0 1 1; -1 0 0; 1 1 1];
>>D=zeros(3,2);
システムの定義．Scilab では syslin(’c’,A,B,C,D)．
>>sys = ss(A,B,C,D)
u(t) がインパルス，ステップの場合．
>>impulse(sys)
>>step(sys)
x(t) の初期値を指定した場合．結果を変数に保存することも可能．
>>[x,t]=initial(sys,[2;3;5]);
>>plot3(x(:,1),x(:,2),x(:,3));
諸条件を指定したシミュレーション．
>>t=0:0.01:10;
>>u=[sin(t);cos(t)];
>>x0= [ 1 ;-5; 0];
Scilab だと，y=csim(u,t,sys,x0); plot3d( y(1,:),y(2,:),y(3,:) ); など．
>>lsim(sys,u,t,x0);
9
10
まず，以下の内容のファイルを myfunc.m という名前で作る．
% 微分方程式の解法の説明用の関数
% このコメント部分が help myfunc で表示される．
% 日本語だとたぶん文字化け
function ydot = myfunc(t,y)
y1=y(2);
y2=2*(1-y(1)^2)*y(2)-y(1);
ydot = [ y1; y2];
>>ts=0:0.1:10;
>>y0 = [2;0];
>>[t,y] = ode45(’myfunc’,ts,y0);
>>plot(t,y(:,1),t,y(:,2))
7
少しだけ古典制御
ここは Scilab では試していません．
1 入力 1 出力のシステムの入力を u(s)，システムの出力を y(s) とすると，伝
達関数 G(s) は以下のように定義される．
G(s) =
y(s)
u(s)
システムが線形でプロパーである場合，伝達関数は次の形で表現できる．
G(s) =
al sl + al−1 sl−1 + · · · + a1 s + a0
bm sm + bm−1 sm−1 + · · · + b1 s + b0
>>impulse(ss3)
>>step(ss3)
>>tf(ss3)
>>tf3
入力と初期状態を指定してシミュレーション．
>>t=0:0.02:30;
>>u=sin(t);
>>x0=[ 10 -3 0];
>>lsim(tf3,u,t,x0)
8
その他・断片的なこと
• 附属のエディタで AfileText を編集する方法．edit myfunc で myfunc と
いう関数を編集．
ここで，l ≤ m である．
分子分母の多項式の係数を指定して伝達関数を定義する．第一引数が分子，第二
引数が分母．係数は次数の高い項のものから並べる．
>>tf1=tf([1],[1 0])
>>tf2=tf([2 3 0],[4 3 2 1 ])
指定された零点，極，ゲインを持つ伝達関数の分子と分母の多項式の係数を計算
する．第一引数が零点，第二引数が極，第三引数がゲイン．
>>[num den] = zp2tf( [3 1]’ , [-3 -2 -5 -1]’ , 1.11)
>>tf3=tf( num, den )
極と零点．
>>pole(tf3)
>>zero(tf3)
ボード線図とナイキスト線図．
>>bode(tf3)
>>nyquist(tf3)
インパルス応答とステップ応答．
>>impulse(tf3)
>>step(tf3)
>>edit AfileText
• コマンドライン操作のショートカットは UNIX(bash や tcsh) に近い
状態空間表現が
d
x = Ax + bu , y = cx
dt
であるシステムの伝達関数は
G(s) = c(sI − A)−1 b
で与えられる．特定の状態空間表現に対する伝達関数は一意に定まるが，特定の
伝達関数に対する状態空間表現は一意ではないことに注意．
状態空間表現への変換．
>>ss3=ss(tf3)
>>pole(ss3)
>>zero(ss3)
11
12