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

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

はじめに
登場人物
【再掲】微分解決済みの状態空間モデル
MATLABコード
コードを見た感想
シミュレーション結果
まとめ

はじめに

前回までで、状態空間モデルの微分解決を実施。
精度としてはシンプルにオイラー法で。(テイラー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;