MATLAB,Python,Scilab,Julia比較 その58【状態空間モデル⑯】

MATLAB,Python,Scilab,Julia比較 その58【状態空間モデル⑯】 数値計算
MATLAB,Python,Scilab,Julia比較 その58【状態空間モデル⑯】

バックナンバーはこちら。
https://www.simulationroom999.com/blog/compare-matlabpythonscilabjulia-backnumber/

はじめに

前回までで、状態空間モデルの微分解決を実施。
精度としてはシンプルにオイラー法で。(テイラー1次)

これをプログラムに落とし込むとどうなるか。

登場人物

博識フクロウのフクさん

指差しフクロウ

イラストACにて公開の「kino_k」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=iKciwKA9&area=1

エンジニア歴8年の太郎くん

技術者太郎

イラストACにて公開の「しのみ」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=uCKphAW2&area=1

【再掲】微分解決済みの状態空間モデル

フクさん
フクさん

まずは微分解決済みの状態空間モデルを再掲しよう。

状態方程式
\(\boldsymbol{x}(t+\Delta t)=\boldsymbol{x}(t)+\{A\boldsymbol{x}(t)+B\boldsymbol{u}(t)\}\Delta t \)

出力方程式
\(\boldsymbol{y}(t+\Delta T)=C\boldsymbol{x}(t+\Delta t)+D\boldsymbol{u}(t)\)

太郎くん
太郎くん

これを、プログラム化か・・・。

太郎くん
太郎くん

まずはMATLABからってことになるのかな?

フクさん
フクさん

そうだね。
そこらへんの流れは従来と一緒にしておこう。

MATLABコード

フクさん
フクさん

運動方程式を微分解決済みの状態空間モデル用いてMATLABで書くとこうなる。

statespacemodel.m

function [x,y] = statespacemodel(A, B, C, D, u, dt, x)
    % 様態方程式
    x = x + (A*x + B*u) * dt;
    
    % 出力方程式
    y = C*x + D*u;
end
m=1;

A=[0,0 ; 1,0];
B=[1/m ; 0];
C=[1,0;0,1];
D=[0;0];
dt=0.001;

t=linspace(0,10,10000); % 時間(横)軸
u=zeros(1,10000);	% 入力信号生成
u(1,5000:10000)=1;	% 5秒後に0から1へ
y=zeros(2, length(t));
x=zeros(2,1);

for i = 1:length(t)
    [x,y(:,i)] = statespacemodel(A,B,C,D,u(i),dt,x);
end

hold on;
plot(t,y,'linewidth',3);
plot(t,u,'--b','linewidth',3);
ylim([-1,14]);
grid();
hold off;

コードを見た感想

太郎くん
太郎くん

statespacemodel関数が状態空間モデルを演算しているところか。
たしかに元にした数式と同じだね。

フクさん
フクさん

そうそう。
statespacemodel関数は変更せず、
入力する、各行列によって振る舞いが大きく変わる感じだ。

太郎くん
太郎くん

なるほど。
これは確かに便利な使い方な気はするね。

フクさん
フクさん

扱う変数はベクトルと行列なんで、多変量になってもロジックは変わらない
まぁこういう書き方のルールがあると思えば良いかな。

太郎くん
太郎くん

そのルールに即していれば、楽に処理ができるってことか。

シミュレーション結果

フクさん
フクさん

そしてシミュレーション結果。

微分解決済み状態空間モデルで運動方程式(MATLAB版)、力F、速度v、距離s
太郎くん
太郎くん

結果も想定通りだね。

まとめ

フクさん
フクさん

まとめだよ。

  • MATLABでベクトル、行列演算による状態空間モデルの演算実施。
  • 導出した数式のまんまでコードが組める。
    • このルールに即していれば、さまざまな振る舞いを規定できる。
  • シミュレーション結果も想定通り。

バックナンバーはこちら。

コメント

タイトルとURLをコピーしました