MATLABの使い方10

リアプノフ方程式を利用した漸近安定の判断

PA + A'P = -Q のPを求める

A = 線形システムの行列A

Q = eye(Aのサイズと同じ)

P = lyap(A',Q)
eig(P)

使用した関数

• lyap(M,Q) : PM' + MP = -Q という設定で、正定対称行列が返ってくる : M = A' という ように この関数のMは線形システムのAの転置となっているので注意

MATLABの使い方9

可観測性

状態空間表現からA,Cを使用する
Vo = obsv(A,C)
rank(Vc) %このランクがVoの列と同じなら可観測
det(Vc) %この値が0でなければ可観測(1入力の場合)

使用した関数

• obsv(A,C) : k可観測行列を求める A,Cは状態空間表現から

オブザーバーゲインの求め方

双対システムの特徴を利用し A+B*K用のアルゴリズムをA=A',B=C'に置き換えることで利用する

状態空間表現からA,Cを使用する
L = -acker(A',C',q)'

上記の作業でオブザーバーゲインを求めることができる
注意
A+B*K 置き換えで A'+C'L' なのでこれで求められるものはL'なので最後に転置することでLがわかる

MATLABの使い方8

可制御性

可制御を求める方法

状態空間表現からA、Bを使用する
A = 任意の行列
B = 任意の行列
Vc = ctrb(A,B)
rank(Vc) %この結果がVcの行サイズと同じなら可制御
[n,m] = size(Vc) % n:行サイズ, m:列サイズ
rank(Vc) == n %結果が1なら可制御,0なら不可制御
det(Vc) %この値が0でなければ可制御(1入力の場合)

使用した関数

• ctrb(A,B) : 可制御行列を求める A,Bは状態空間表現のA行列、B行列
• rank(行列) : 行列のランクを求められる => 行フルランクなら可制御(行列の行数とランクが同じ場合)

極配置法(1)[直接的な方法(1入力システム)]

%行列の定義-----------------------
A = [-3 1;2 -2];
B = [2;0];
%変数(シンボリック)の定義---------
syms lambda k1 k2 p1 p2
----------------------------------
%ゲイン行列Kと併合システムのAcl行列の定義----
K = [k1, k2];
Acl = A + B*K;
%--------------------------------------------
%特性多項式(|λI-Acl|を求める[k1,k2]を利用したパターン)
eq1 = det(lambda*eye(2) - Acl); %行列式
eq1 = collect(eq1,lambda)       %lambdaのべき乗で係数をまとめる
coe1 = coeffs(eq1,lambda)       %lambdaの係数をすべて返す(次元の低いものから返される)
%--------------------------------------------
%多項式((λ-p1)(λ-p2))のパターン
eq2 = (lambda - p1)*(lambda - p2);
eq2 = collect(eq2,lambda)
coe2 = coeffs(eq2,lambda)
%--------------------------------------------
%方程式を指定したシンボリックに対して解く------------------
[k1,k2] = solve([coe1(1)-coe2(1),coe1(2)-coe2(2)],[k1,k2])
%----------------------------------------------------------
%k1,k2の関数内のp1とp2に指定の極を与えてそれに当てはまるk1,k2を求める---
K = subs([k1,k2],{p1,p2},{-8+4j,-8-4j})
%-----------------------------------------------------------------------
%Aclの固有値を求める----
eig(A+B*K)
%-----------------------

使用した関数

• collect(方程式,シンボリック) : 指定したシンボリックに関してべき乗で方程式内の係数をまとめる(今回の場合lambdaを指定しているのでlambdaのべき乗の式になる) >lambda^2 +(係数)*lambda のような形になる(デフォルトのシンボリックはx)
• coeffs(方程式,シンボリック): 指定した方程式のシンボリックの係数を次元の低いものから配列にして返す(通常みる多項式の後ろからということ)
• solve(数式,シンボリック) : 数式について指定したシンボリックで解を求める(数式とシンボリック)を複数やる場合は対応するもの同士を同じ配列の場所に入れることで複数の解を求めることができる
• eig(行列) : 行列の固有値を求める

極配置(2)[可制御標準形に基づく方法(1入力システム)]

%ステップ1---------------------------
coe = poly(A)               %特性多項式det(λI-A)の係数(高次から低次の順で返される)
a1 = coe(2); a0 = coe(3);
%------------------------------------
%ステップ2---------------------------
Mc = [a1 1;1 0]            %正則行列Mcの計算
Vc = ctrb(A,B)             %可制御行列の計算
%ステップ3---------------------------
Tc = inv(Vc*Mc)            %可制御標準形への変換行列
%ステップ4---------------------------
p1 = -8+4j; p2=-8-4j;
Delta = conv([1 -p1],[1 -p2]) %多項式の計算(因数分解されたものの展開)
d1 = Delta(2); d0 =Delta(3);
%------------------------------------
%ステップ5---------------------------
Kc = [a0-d0 a1-d1]
K  = Kc*Tc
%A+B*Kの固有値を
eig(A + B*K)
%-----------------------------------

使用した関数

• poly(行列) : 引数に行列を与えた場合、特性多項式(|λI-A|)の係数を返す=>λ2+係数*λのような形(λの先頭の係数から返されるため今回の場合先頭のλ2の1から順に返される)　順番は高次から低次となっている
• inv(行列) : 行列の逆行列を求める
• conv(配列, 配列) : 畳み込みおよび多項式の乗算をする(多項式の分配法則的なこと) 配列には計算をする多項式の係数が入った配列を与える

極配置(3)[アッカーマンの極配置アルゴリズム(1入力システム)]

A = [-3 1;2 -2];
B = [2;0];

p = [-8+4j;-8-4j];
K = - acker(A,B,p)

eig(A+B*K)

使用した関数

• acker(A,B,p) : 1入力システムのみ対応 A,B行列と指定したい極pを配列で渡す
  注意
  ackerから返ってくるゲインはu=-Kxを想定しているのでKに代入する際に負の値をかける必要がある

極配置(4)[多入力システム]

A = [0 1 0 0; -1 -1 1 1; 0 0 0 1;2 2 -2 -2];
B = [0 0 ;1/2 0; 0 0; 0 1];

p(1) = -2; p(2) = -2;
p(3) = -4; p(4) = -4;

K = -place(A,B,p)
eig(A+B*K)

使用した関数

• place(A,B,p) : 多入システムに対応　基本的にackerと同じ
  注意
  ackerと同様でKに代入する際にマイナスをかける(符号を反対にする)必要がある

基本的にackerとplaceで極配置できるという認識でOK

MATLABの使い方7

ヘビサイドの公式を利用した係数行列の求め方

行列などで複数の係数行列を求めるときに便利

%(sI-A)^-1 の係数行列を求めるのを目標とする

A = [0 1; -10 -11]; %ここは任意でOK

syms s %ラプラス演算子sを定義

p = eig(A); %Aの固有値を求める
p1 = p(2)
p2 = p(1)
I  = eye(2) %2×2の単位行列を作成
Q = inv(s*I - A) % (s*I -A)の逆行列を求める
%ヘビサイド
eq1 = simplify((s-p1)*Q); %(s-p1)Q(s)を簡略化=>分母分子の共通因子を約分
eq2 = simplify((s-p2)*Q); %上記と同じ

K1 = subs(eq1, s, p1) %eq1のsにp1を代入する
K2 = subs(eq2, s, p2)

注意 上記の作業でK1などの値が分数で表示された場合
double(K1)とすることで小数表示にできる

使用した関数

• eig(行列) : 行列の固有値を求める
• eye(数値) : 数値×数値の単位行列を作成する
• inv(行列) : 逆行列を求める
• simplify(数式) : 約分をする
• subs(シンボリック式, 代入させる変数, 代入する値) : シンボリック式の変数に値をいれて計算する

その他の行列の操作

• det(行列) : 行列式を求める

その他の代数計算系

• factor(数式,シンボリック) : 因数分解をする => シンボリックのところには対応する変数を入れる(デフォルトはx)

遷移行列を求める

expm(行列) : 行列の部分には基本的にA*tだがtにシンボリックを使用することもできる
注意
シンボリックを使用してsubsで計算した場合複素数表記になってしまう問題がある

sin,cosで表記する方法

A = 行列
syms s
exp_At = ilaplace(inv(s*eye(2)-A)) %ラプラス変換で求める方法を利用する

注意
この方法を利用してもsubsの計算がされない問題があるので注意する

時間応答

パターンとして零入力応答,零状態応答,任意の時間応答を求める方法

• initial(状態空間表現,初期値,時間配列) : 零入力応答のグラフをplotする 状態空間表現はssで作成したも、初期値は縦行列(状態変数),tは配列を0:0.01:5のような感じで作成したもの
• y=initial(上記と同じ) : yに時間応答の値を代入する
• step(状態空間表現,時間配列) : 単位ステップ応答のグラフをplotする => initialと同じでyに代入も可能
• impulse(状態空間表現,時間配列) : インパルス応答のグラフをするをplot => yに代入可能
• lsim(状態空間表現,入力値,時間配列,初期値) : 入力値は時間配列と同じサイズの配列にする必要があるため別で入力値の配列を作成する
  また、lsimでのplotだと入力値と応答の両方がグラフに表示される => yに代入可能でyには応答のデータのみ入る
  補足
  上記の操作はすべて多入力多出力に対応しており、出力yに配列としてデータが返される

漸近安定,BIBO安定

• 漸近安定を求める => eigを利用して固有値を求め,実部がすべて負なら漸近安定
• BIBO安定 => zpkを利用して極零点をだし、極がすべて負ならBIBO安定

MATLABの使い方6

現代制御編

状態空間表現の定義

• ss(A,B,C,D) : ss関数に各行列を与えることで状態空間表現を作成できる

状態空間表現から係数行列を取り出す

• [A,B,C,D] = ssdata(状態空間表現) : ssdata関数を使用することで係数行列を取り出すことができる

状態空間表現から伝達関数表現への変換

• tf(状態空間表現)
• zpk(状態空間表現) : こちらは零極ゲインの形になっているものに変換する 上記の操作後は古典制御でのtfやzpkの使いかたでOK

伝達関数表現から状態空間表現への変換

• ss(伝達関数表現) : この変換で得られるものは可制御標準形ではない
• 可制御標準形の求め方

[numP,denP] = tfdata(伝達関数表現)
[A, B, C, D] = tf2ss(numP,denP)
ss_P = ss(A,B,C,D) %tf2ssで作成した行列から状態空間表現を作成する

%ここまでだと可制御標準形と状態空間表現の要素の順番が逆になっている

T = [0 0 1; 0 1 0; 1 0 0]
ss_Pb = ss2ss(ss_P,T) %状態空間表現の変換
%ここで伝達関数表現から可制御標準形を求めることができた

使用した関数

1. tfdata() : 伝達関数から分子要素分母要素を取り出す
2. tf2ss() : 伝達関数表現(の分子分母要素)から状態空間表現の各行列を求める
3. ss2ss() : 状態空間表現を変換する 第二引数には変換行列を与える

最小実現

• ss(状態空間表現,'min') : 通常の状態空間表現から最小実現の状態空間表現を求める
  使用方法

ss_P = ss(A,B,C,D) %通常と同じように行列から状態空間表現を作成する
ss_P = ss(ss_P,'min') %最小実現の状態空間表現を得る

MATLABの使い方5

時間付き構造体の作成方法

時間配列と、値のデータ配列から時間付き構造体を作成する

• timeseries(値データ配列,時間配列) : これで時間付き構造体になったものが返される

上記で作成したデータを.mat形式で保存する場合

• save('保存するファイル名.mat','保存する時間付き構造体') : 左のような関数を実行することで保存することができる

ワークスペース変数の保存

• save('保存するファイル名.mat','変数1','変数2',...,'変数n')

MATLABの使い方4

ボード線図

bode(伝達関数) : でボード線図を作成できる

使用例

• bode(伝達関数,周波数範囲) : 周波数範囲は w = logspace(-2,2,100) のように 指定する. 今回の場合 -2 ~ 2 なので 10^-2 ~ 10^2 の範囲となる
• [Gg,Gp] = bode(伝達関数,周波数範囲) : Ggの方にゲイン特性,Gpの方に位相特性が入る その後の表示方法は

figure(1); semilogx(周波数範囲,20*log*10(Gg(:,:))); %ここでゲイン特性(dB表示)をグラフ1に描画
figure(2); semilogx(周波数範囲,Gp(:,:)); %ここで位相特性をグラフ2に描画

• [Gg,Gp] = bode(伝達関数の分子要素,伝達関数の分母要素,周波数範囲):　以上のコマンドでもGg,Gpを求めることが可能

周波数応答

lsimでsin波を指定すればOK
定常値はbodeを使用して公式の式にその値を使用することで求められる

使用方法

sysP = tf(numP,denP); %伝達関数の作成
w = 2; %入力周波数の設定
t = 0:0.01:20; %時間の設定
[Gg,Gp] = bode(sysP,w);
Gp_rad = Gp*pi/180 %位相差の単位を[deg]から[rad]へ

%入力に対する出力y(t)の定常値yapp(t)の作成----

yapp = Gg * sin(w*t + Gp_rad);

%---------------------------------------------
plot(t,yapp) %定常値のグラフを表示

u = sin(w*t);
lsim(sys,u,t) %lsimを利用した任意の入力に対する時間応答のグラフを出力

ポイント
同じグラフに表示させたい場合
hold on を使用することでこの後に出てくるグラフ出力は同じグラフに描画される
hold off を使用することでhold onの機能を切ることもできる

共振ピーク,共振周波数

[Mp,i] = max(ゲイン特性) : ゲイン特性にはbodeで求めたものを使用する Mpには共振ピークの値,iにはbodeで使用した周波数配列のインデントの番号が入る

使用方法

sysP = tf(numP,denP)
w = logspace(-1,1,10000) %周波数wの指定
[Gg,Gp] = bode(sysP,w)
[Mp, i] = max(Gg) %Mp:共振ピーク i:周波数配列のインデントが入る
wp＝w(i); %wp : 共振周波数

ベクトル軌跡,ナイキスト軌跡

nyquist(伝達関数) : これでナイキスト軌跡のグラフが表示される

MATLABの使い方3

ブロック線図の結合

結合させる伝達関数を
sysP1 = tf([1],[1 3 2])
sysP2 = tf([1],[10, 1])
とする

• 直列結合 : sysP1 * sysP2
• 並列結合 : sysP1 + sysP2 または、 sysP1 - sysP2
• フィードバック結合 : feedback(sysP1, sysP2) =>　分子がsysP1 分母が1+sysP1*sysP2となる
  イメージはフィードバックのない上の部分にsysP1があり、フィードバックの部分にsysP2がある感じ

フィードバック制御系の目標値応答,入力外乱応答

コントローラーを
sysC = 2
制御系の伝達関数を
sysP = tf([5],[1 2 2])
とする

• Gyr(s) (目標値から出力) : feedback(sysP*sysC, 1)
• Gyd(s) (入力外乱から出力) : feedback(sysP, sysC)
• Ger(s) (目標値から偏差) : 1 - Gyr(s) または feedback(1,sysP*sysC)
• Ged(s) (入力外乱から偏差) : - Gyd

根軌跡

• rlocus(伝達関数) : 根軌跡を描くことができる
  グラフィカルな調整がしたい場合は、制御システムアナライザーを使用する

MATLABの使い方2

ラプラス変換・逆ラプラス変換

• laplace(f,var,transVar) : ラプラス変換をする
  f -> シンボリック式、シンボリック関数、あるいはシンボリック式な関数のベクトルor行列 var -> 独立変数を表すシンボリック変数(default:t) 時間変数のこと transVar -> 変換変数を表すシンボリック変数またはシンボリック式(default:s)ラプラス演算子

• ilaplace(F, var, transVar) : 逆ラプラス変換をする
  F -> 上記のラプラス変換のfと同じ var -> ラプラス演算子を表すシンボリックのこと(default:s) transVar -> 時間変数(default:t)

ラプラス変換・逆ラプラス変換の注意

伝達関数のために行列の要素ごとの配列を作成し、tfで伝達関数を作成したが、その伝達関数をそのままlaplaceに代入することはできないので注意が必要

対処法
分子要素をnumY , 分母要素をdenY　とする

ex1)逆ラプラス変換をtfに対してらくにする
  new_numY = 0;
  syms s;
  for i = 1:length(numY)
    new_numY = new_numY + numY(i)*s^(length(numY)-1);
  end
  %上記の作業でnew_numYにシンボリック関数として伝達関数の分子要素が保存された
  %同じ作業をdenYにも行いnew_denYを作成
  ilaplace(new_numY/new_denY)
  
ex2)逆ラプラス変換の結果から値を代入し、グラフを作成したい
逆ラプラス変換の結果をsysPとする
  
  T = 0:0.01:10; %時間の値
  data_sysP = subs(sysP, t, T);
  %上記の作業でsysP内の時間変数tを実数Tで置き換える
  plot(T,data_sysP)
  %これでグラフ表示完了

ex3)ex2とは違う方法で逆ラプラス変換の結果からのグラフ作成
  func_sysP = symfun(sysP, [t])
  %上記の作業で逆ラプラス変換の結果がtを引数にする関数にできる
  data_sysP = func_sysP(T);
  plot(T,data_sysP)
  %これでグラフ表示完了

部分分数分解

residue(分子係数,分母係数) : 部分分数分解の結果を分子要素と分母要素で返してくれる

ex)
numY = [4 5]; denY = [1 3 2];
[k p] = residue(numY, denY)
k -> 分子の要素k
p -> 分母要素p (なので実際の式は(s-p)となる)

重解がある場合

重解の累乗が小さい方(1乗->2乗->...)から返される
分母要素のpは同じ値が返ってくる

ex) k_12/(s-p_1) + k_11/(s-p_1)
の場合に上記の例のように [k p] = residue()した場合
k = [k_11, k_12]
p = [p_1, p_1]
といった形で値が返ってくる

インパルス応答

impulse() : インパルス応答を得ることができる　引数には基本的にtfやzpkで作成した伝達関数

使用方法

• impluse(sys) : 伝達関数のみを引数に与える インパルス応答のグラフが作成される
• impluse(sys,Tfinal) : Tfinalに指定した時間までのインパルス応答のグラフが作成される(0~Tfinal)の範囲
• impluse(sys,t): tをt = 0:0.01:10などで作成し、その範囲のインパルス応答を得る
• y = impluse(sys,t) : インパルス応答の計算結果をyに出力しグラフは表示されない
• impluse(sys1,sys2,...sysN,t): 指定した時間の範囲の複数の伝達関数のインパルス応答をグラフに描画する

ステップ応答

step() : ステップ応答を得ることができる. 使用方法は上記のimpulse()と同じ

任意の入力に対する応答

lsim() : 任意の入力に対しての応答を得ることができる

使用方法

• lsim(sys,u,t) : 任意の入力u にたいしての時間応答をプロットする(tは上記のimpluseと同じ)
• lsim(sys,u,t,x0) : x0は初期条件を指定できる
• lsim(sys,u,t,method): sysが連続システムである場合に、入力値がサンプル間でどのように内挿されるかを明示的に指定できる
  methodの種類
  • 'zoh' : ゼロ次ホールドを使用
  • 'foh' : 線形内挿(1次ホールド)を使用
• lsim(sys1,sys2,...,sysN,u,t) : 複数の伝達関数でもOK
• lsim(sys1,PlotStyle1,...,sysN,PlotStyleN,u,t) : プロットにおける各システム応答のスタイルを指定できる
  ex) lsim(sys1, 'y:', sys2,'g--',u,t,x0)

MATLABの使い方1

伝達関数表現

• tf : 伝達関数を作成する

ex1)
 numP = [4,8];
 denP = [1,2,-15,0];
 sysP = tf(numP,denP)
 >>>
 sysP = 
        4 s + 8
  ------------------
  s^3 + 2 s^2 - 15 s

ex2)
 s = tf('s');
 sysP2 = (4*s + 8)/(s^3 + 2*s^2 - 15*s) 
 >>>
 sysP2 =
 
       4 s + 8
  ------------------
  s^3 + 2 s^2 - 15 s

逆に伝達関数が与えらており、その分子要素、分母要素が知りたいとき

%伝達関数をsysPとする
[numP, denP] = tfdata(sysP, 'v') %分子、分母の順で返される

• zpk : 極、零点表現した伝達関数を作成する

ex1)
 z = [-2]; %零点
 p = [-5 0 3]; %極
 K = 4; % ゲイン
 sysP3 = zpk(z,p,K)
 >>>
 sysP3 =
 
     4 (s+2)
  -------------
  s (s+5) (s-3)

上記の二つを利用してtfで作成した伝達関数をzpkの中に入れることで極、零点表現した伝達関数に変換できる