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

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

はじめに
登場人物
【再掲】DCモータ状態空間モデル
Juliaコード
シミュレーション結果
まとめ

はじめに

前回は、ScilabによるDCモータ状態空間モデルのシミュレーションを実施。

今回は、これのJulia版

登場人物

博識フクロウのフクさん

イラスト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

【再掲】DCモータ状態空間モデル

太郎くん

最後はJuliaでDCモータ状態空間モデルシミュレーションだね。
まぁ雰囲気はMATLABにかなり似てはいるよねー。

フクさん

似てはいるが故の注意点もあるけどね。
まぁ恒例のDCモータ状態空間モデルの再掲だ。

状態方程式

\(
\begin{bmatrix}
\dot{\theta}(t) \\
\dot{\omega}(t) \\
\dot{I}(t)
\end{bmatrix}=
\begin{bmatrix}
0 && 1 && 0 \\
0 && 0 && K/J \\
0 && -K/L && -R/L
\end{bmatrix}
\begin{bmatrix}
\theta(t) \\
\omega(t) \\
I(t)
\end{bmatrix}+
\begin{bmatrix}
0 \\
0 \\
1/L
\end{bmatrix}
E(t)
\)

出力方程式

\(
\boldsymbol{y}=
\begin{bmatrix}
1 && 0 && 0 \\
0 && 1 && 0 \\
0 && 0 && 1 \\
\end{bmatrix}
\begin{bmatrix}
\theta(t) \\
\omega(t) \\
I(t)
\end{bmatrix}+
\begin{bmatrix}
0 \\
0 \\
0
\end{bmatrix}
E(t)
\)

Juliaコード

フクさん

そしてJuliaによるシミュレーションコード

using PyPlot

function statespacemodel(A,B,C,D,u,dt,x)
    # 状態方程式
    x = x + (A*x + B*u) * dt;
    
    # 出力方程式
    y = C*x + D*u;
   
    return x,y;
end

function statespacemodel_motor()
    K=0.016;
    J=0.000000919;
    R=1.34;
    L=0.00012;

    A=[0 1 0 ; 0 0 K/J ; 0 -K/L -R/L];
    B=[0 ; 0; 1/L];
    C=[1.0 0 0; 0 1.0 0; 0 0 1.0];
    D=[0.0 ; 0.0 ; 0.0];
    
    dt = 0.0001;
    t = range(0, 10-dt, step=dt); # 時間(横)軸
    u = zeros(1,100000);          # 入力信号生成
    u[1,50000:100000].=1.0;       # 5秒後に0から1へ
    y = zeros(3,length(t));
    x = zeros(3,1);
    
    for i = 1:length(t)
        x,y[:,i] = statespacemodel(A,B,C,D,reshape(u[:,i],1,1),dt,x);
    end
    subplot(3,1,1);
    plot(t,u[1,:]);
    ylim([-0.05,1.1]);
    xlim([4.95,5.25]);
    subplot(3,1,2);
    plot(t,y[1:2,:]');
    ylim([-0.05,80]);
    xlim([4.95,5.25]);
    subplot(3,1,3);
    plot(t,y[3,:]);
    xlim([4.95,5.25]);
    
    return y;
end

statespacemodel_motor();

太郎くん

予想通りではあるけど、かなりMATLABっぽい。
差は、linspaceがrangeになってるくらいだね。