2 次遅れ系のパラメータをダイヤルで変化させながら Simulink でステップ応答を確認！

# <font color="#33CC33">1.　はじめに</font>

　「制御工学」で学習する 2 次遅れ系
```math
\begin{align}
    y(s) = P(s)u(s),\ 
    P(s) = \dfrac{K{\omega}_{\rm n}^{2}}
            {s^2 + 2\zeta{\omega}_{\rm n}s + {\omega}_{\rm n}^{2}}
    \tag{1}
\end{align}
```
を考えます．ここで，

* $\zeta > 0$：減衰係数 $\cdots\cdots$ 安定度のパラメータ
* ${\omega}_{\rm n} > 0$：固有角周波数 $\cdots\cdots$ 速応性のパラメータ
* $K \neq 0$：ゲイン $\cdots\cdots$ 定常特性のパラメータ

です．また，信号 $u(t)$, $y(t)$ のラプラス変換を $u(s) = {\cal L}\bigl[u(t)\bigr]$, $y(s) = {\cal L}\bigl[y(t)\bigr]$ のように小文字で記述しています．

　本記事では，$\zeta$, ${\omega}_{\rm n}$, $K$ を変化させたときに，$(1)$ 式の単位ステップ応答 $y(t)$ がどのように変化するのかを，Simulink で確認する例を説明します．

　なお，<font color="#E84746"><b>Simulink モデルは GitHub で公開</b></font>しています（R2023a, R2022a, R2020a 用）．

 https://github.com/Prof-Carter/step_2nd_system

 

 # <font color="#33CC33">2.　2 次遅れ系の Simulink での表現</font>

　単に，Simulink を利用して 2 次遅れ系のステップ応答を描画するのであれば，Simulink ブロック "$\tt Transfer\ Fcn$" を利用して

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2552191/1b16d1d6-fb9a-f6be-366d-76d00d2b89f3.png)

のようにすれば良いです．しかし，今回は，<font color="#E84746"><b>シミュレーションを行っている最中に $\zeta$, ${\omega}_{\rm n}$, $K$ の値を変化させて，その振る舞いの変化を確認する</b></font>ことを目的としているために，

* 2 次遅れ系を微分方程式で書き換え，Simulink ブロック "$\tt Fcn$" で記述する
* Simulink ブロック "$\tt Fcn$" で微分方程式を記述する際にパラメータ $\zeta$, ${\omega}_{\rm n}$, $K$ を信号 $u(t)$, $y(t)$, $\dot{y}(t)$ と同様に変数として扱う

というアプローチを考えます．まぁ，大したことを行うわけではありません．

 　2 次遅れ系 $(1)$ 式を書き換えると，
```math
\begin{align*}
    &
    (s^2 + 2\zeta{\omega}_{\rm n}s + {\omega}_{\rm n}^{2})y(s) 
        = K{\omega}_{\rm n}^{2}u(s)
    \\
    &\quad \Longrightarrow \quad
    s^2y(s) 
        + 2\zeta{\omega}_{\rm n}sy(s) 
        + {\omega}_{\rm n}^{2}y(s)
        = K{\omega}_{\rm n}^{2}u(s)
    \tag{2}
\end{align*}
```
なので，$(1)$ 式は 2 階の線形微分方程式
```math
\begin{align*}
    &
    \ddot{y}(t) 
        + 2\zeta{\omega}_{\rm n}\dot{y}(t) 
        + {\omega}_{\rm n}^{2}y(t)
        = K{\omega}_{\rm n}^{2}u(t)
    \\
    &\quad \Longrightarrow \quad
    \ddot{y}(t) 
        = K{\omega}_{\rm n}^{2}u(t)
        - 2\zeta{\omega}_{\rm n}\dot{y}(t) 
        - {\omega}_{\rm n}^{2}y(t)
    \tag{3}
\end{align*}
```
のように記述できます．そこで，Simulink ブロック "$\tt Fcn$" を利用して $(3)$ 式を表現します．"$\tt Fcn$" の使い方については

https://qiita.com/Carter/items/3f11c2665af9196a8295

で詳しく説明しています（<font color="#E84746">R2020b から R2022b では "$\tt Fcn$" を Simulink ライブラリブラウザからは検索できませんので，上記の記事を参照して simulink3 から "$\tt Fcn$" を移動してください</font>）．

　さて，"$\tt Subsystem$" の中に以下のように Simulink ブロックを配置し，パラメータや名前を設定した上で結線します．

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2552191/735a575a-7eac-e6e2-70a5-11b08b315056.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2552191/aba70844-93e0-8d35-c2b4-58c8d8212a72.png)


"$\tt Fcn$" の入力は "$\tt Mux$" でベクトル化された信号であり，それぞれ

* $\tt u(1)$：$y(t)$
* $\tt u(2)$：$\dot{y}(t)$
* $\tt u(3)$：$u(t)$
* $\tt u(4)$：$\zeta$
* $\tt u(5)$：${\omega}_{\rm n}$
* $\tt u(6)$：$K$

を意味しています．そして，これらを用いて $(3)$ 式の右辺
```math
\begin{align*}
    K{\omega}_{\rm n}^{2}u(t)
        - 2\zeta{\omega}_{\rm n}\dot{y}(t) 
        - {\omega}_{\rm n}^{2}y(t)
\end{align*}
```
を "$\tt Fcn$" では
```math
\begin{align*}
    \tt{
    u(6)*u(5)^2*u(3) - 2*u(4)*u(5)*u(2) - u(5)^2*u(1)
    }
\end{align*}
```
のように表現しています．したがって，"$\tt Fcn$" の出力は $\ddot{y}(t)$ となりますので，"$\tt Integrator$" を利用してこれを 1 回積分すると $\dot{y}(t)$，2 回積分すると $y(t)$ が得られます．これらをフィードバックして "$\tt Mux$" に入力しています．

　以上で，2 次遅れ系を Simulink で表現することができました．


# <font color="#33CC33">3.　Simulink モデル</font>

 ## <font color="#33CC33">3.1　完成した Simulink モデル！！</font>

　配布する Simulink モデル "<b>step_2nd_system.slx</b>" を以下に示します．

 ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2552191/0df478f6-c908-2e87-eefe-c6873595e6e1.png)

## <font color="#33CC33">3.2　Simulink ブロック の設定</font>
 
### <font color="#33CC33">3.2.1　Scope の設定</font>

　シミュレーションを行っている間，リアルタイムで信号を表示するために "$\tt Scope$" を利用します．"$\tt Scope$" をダブルクリックするとウィンドウが起動しますので，マウスを右クリックすることで，コンフィギュレーションプロパティを

 ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2552191/40405713-aed2-4b4b-75c3-a570fae3609b.png)

のように設定します．これにより，入力信号 $u(t)$ と出力信号 $y(t)$ を同じウィンドウに 10 秒ごとにリフレッシュしながら表示できるようになります．また，縦軸の範囲を $-1$ から $2.5$ に設定しています．

　同様に，スタイルを

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2552191/b9b58a69-888d-2af7-42a6-0892603557db.png)

とし，入力信号 $u(t)$ と出力信号 $y(t)$ の描画の線の太さや色を設定しています．

 
### <font color="#33CC33">3.2.2　Pulse Generator の設定</font>

　入力信号を，10 秒ごとにステップ状に 0 と 1 の 2 値に変化させます．この 2 値の比は 50% です．また，2.5 秒ずらすことで，ステップ状に変化する "$\tt Scope$" での時間を 10 秒間のうちで 2.5 秒と 7.5 秒にします．
 
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2552191/938369d2-0727-6795-790f-b6a2a0676930.png)

 
### <font color="#33CC33">3.2.3　Knob の設定</font>

　"$\tt Knob$" はダイヤルでパラメータの値を調整するための Simulink ブロックです．今回は，シミュレーションを行っている最中に，"$\tt Knob$" のダイヤルををマウスで回転させることで $\zeta$, ${\omega}_{\rm n}$, $K$ の値を変化できるようにします．
 
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2552191/a06a9ec8-ad3a-ac67-e4fd-83fd07de67a2.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2552191/3419deb9-b20e-53dd-ee02-4414d64248b3.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2552191/dc4ce0d7-7a40-a1cc-769c-49273a0cf9c3.png)


 
## <font color="#33CC33">3.3　モデルコンフィギュレーションパラメータの設定</font>

　サンプリング周期を 0.001 [s] として  $(3)$ 式で示した微分方程式を 4 次のルンゲ=クッタ法で数値的に解くため，以下のように設定します．
 
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2552191/b0a38c5f-0167-a4d3-b860-ae98d782cf99.png)


## <font color="#33CC33">3.4　シミュレーションページングの設定</font>

　Simulink でシミュレーションをすると高速でグラフが描画がされてしまい，変化がわかりにくいです．そこで，<font color="#E84746"><b>シミュレーションページング</b></font>を設定します．

　作成した Simulink モデルの「シミュレーション / 実行」で「シミュレーションページング」を選択すると，下図のウィンドウが開きますので，
 
* 「ページングを有効にしてシミュレーション速度を遅くする」をチェック

としてください．

 ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2552191/55a61ced-59ca-1f0f-3698-d8643f47859a.png)

このとき，標準では，

* 「実経過時間（秒）ごとのシミュレーション時間」が「1」

となっているので，1 秒間のシミュレーションが 1 秒で実行されます（1 倍速）．
今回は，

* 「実経過時間（秒）ごとのシミュレーション時間」を「5」

に変更することで，5 秒間のシミュレーションが 1 秒で実行されるようにします（5 倍速）．

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2552191/9694e9f3-0dd2-b8f0-ef22-1abcf290cc00.png)



# <font color="#33CC33">4.　実行例</font>

　実行した様子の動画を YouTube にアップしていますので，確認してみてください．先に述べたように，実行速度は 5 倍速としています．マウスでダイヤルを回すと，2 次遅れ系の 3 つのパラメータ $\zeta$, ${\omega}_{\rm n}$, $K$ の役割が感じ取れるのではないでしょうか．

 https://youtu.be/JFmkbNfRfhQ


# <font color="#33CC33">5.　おわりに</font>

　ということで，もっと良い方法があるかもしれませんが，ダイヤルでパラメータ ${\omega}_{\rm n}$, $\zeta$, $K$ の値を変化させながら，シミュレーションで結果を確認できる Simulink モデルを作成してみました．

　ちなみに，ここで紹介したした Simulink モデルは，PID 制御の目標値フィルタ（$K = 1$ とした 2 次遅れ要素）のパラメータを実験しながら調整するために作成したものです．

https://twitter.com/carter_control/status/1731268159121059980