MATLAB/Simscape Multibodyでの台車型倒立振子シミュレーション

[MATLAB/Simulink Advent Calender 2022](https://qiita.com/advent-calendar/2022/matlab)の4日目になります。
3日目は、@kentaPtさんの[「iPhone LiDARで取得した3次元点群を地図上にプロットしてみよう」](https://kentapt.hatenablog.com/entry/matlabQiita)でした。MATLABで3次元データも扱えるとは、大体の工学分野はMATLABだけでいいのではないか、という気がしてきます。myスマホがiphone SE(第一世代)なので、試せないのが残念です。

# 背景
MATLABツールボックスの一つであるSimscape Multibodyについて、MathWorksが公開している[入門文書](https://jp.mathworks.com/help/sm/getting-started-with-simmechanics.html?s_tid=CRUX_lftnav)は、最終的にリンク機構を例に説明されていました。
しかし、制御工学勢からするとリンク機構はなじみが薄いです。そこで、みんな大好き倒立振子の実装を通して、別視点からの入門記事を提供します。
フライホイール型の倒立振子は、[コチラ](https://underbirdworks.blogspot.com/2021/12/simscape-multibody.html)で実装された方がいらっしゃるため、今回は台車型の倒立振子を対象とします。

![fig1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/287d4ee0-1f9d-c798-08b6-451919d2ffbf.png)

本記事の対象読者は、MathWorksが公開している次の2つのチュートリアルは実践されたものとしています。
1. [単純なリンクのモデル化](https://jp.mathworks.com/help/sm/gs/model-simple-link.html)
1. [単純な振子のモデル化](https://jp.mathworks.com/help/sm/gs/model-pendulum.html)


# 実装過程
## 前知識
Simscape MultibodyはZ軸を基準としていることが多いです。例えば、回転接続であるRevolute Jointは、問答無用でZ軸回りの回転となります。そのため、各オブジェクトのどこが何軸かを意識しながら、モデル作成を進めることが上達の近道です。

## 基本要素への分解
シミュレーション対象を構成要素に分解していきます。オブジェクトは振子と台車になります。また、オブジェクト間はそれぞれ回転軸および直線軸で接続されています。
+ 振子　→台車と回転軸で接続する
+ 台車　→振子と回転軸で接続する、地上と直線軸で接続する
![fig2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/8cafca6c-9e8d-bca7-b2dc-a26d43cab5c3.png)

## 方針
* はじめに、振子オブジェクトを作成します。
* 次に、単純な振子を作成し、重力によって減衰振動することを確認します。
* 平行して、台車を作成します。
* 最後に台車と振子を繋げます。

## 振子を作る
* smnewコマンドで新しいシミュレーション環境を作成します
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/23c6aebe-c303-ca94-acc1-00d49b6d3545.png)

* チュートリアルで実施したように、まずは棒を作ります。ここでは、アルミ材料を想定したパラメータ設定をしました。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/b12b4c31-509f-4073-4ebd-e737161c88d3.png)

* 今回の振子オブジェクトは、片側しか接続しないため、座標(Rigid Transform)も片端にだけ設定します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/5e9e80ad-a02e-981f-c30a-a97b0c8ddb2c.png)

* 位置だけでなく、角度も同じRigid Transformで変換します。Rotatin About設定がFollower Axesの場合は、座標変換前の座標系で3軸を回転します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/0a89b467-1f79-b128-5c69-a8e5008b1f40.png)

* 地面との間を回転軸(Revolute Joint)で接続します。回転軸はZ軸回りでの回転となるので、回転軸と振子の回転とが合うようにRevolute Jointの前後で座標変換②と戻し③をしています。F⇒Bにした場合は、同じブロックでB⇒Fをすれば、結果として元の座標系に戻ります。座標系が見ずらいため、①のブロックで原点からx方向に20mmほどズラしています。
* 回転軸のダンピング係数(Damping Coefficient)は1e-4 N*m/(deg/s)に設定しました
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/c3384e44-a492-ecab-76b9-246784240755.png)
※Revolute Jointの向きを逆にしていました。以降のモデルでは修正しています。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/5000d672-0f8d-dd2a-3722-834e0649de02.png)

* この時点でのシミュレーション例
![CartPendforQiita.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/eae10b56-9426-523c-2fb9-afa17da7c739.gif)

* 振子オブジェクトが完成したので、今後のことを考えてサブシステム化します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/f947d6f6-cd42-0363-3dfe-871e0a1ed33f.png)


## 台車を作る
* 次は台車単体を実装します。ひとまず振子をコメントアウトします。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/cfbba872-3574-9d9f-2618-73aabc4736d0.png)

* 直線案内(Prismatic Joint)で台車と地面を接続します。Prismatic JointがZ方向のみに移動可能な直線運動を規定するので、Jointの前後に座標変換を加えて、固定座標系のX軸方向の拘束とします。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/77375add-654c-b47e-63ed-ea7fb7d916fc.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/6befaed1-dcc0-751a-5f88-385a74f6c54e.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/b46ce1fe-cbaa-d49c-f218-c43f98bdea10.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/db19cdd1-468b-4cbf-5801-bdcfe8465008.png)
※Rigid Transformの設定は、どちらも同じ設定で、方向(接続端子)のみが違います。

* External Force and Torqueブロックによって、台車に力を加えます。
ここでは、正弦波状にトルクを加えて、台車が動くことを確認しました。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/4626aa0f-32b3-74e9-5dcc-7f6ff5d3ca0c.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/9cfbb3c7-026a-ad18-7805-500040d91110.png)

![CartPendforQiita2.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/eb33bc41-639a-cc0b-f9ec-1c321b167be0.gif)

## 振子と台車をつなげる
台車と振子ができたので、２つを繋げます(実際には振子のコメントアウトを解除します)。正弦波の振幅を大きくしてあげれば、台車の移動に伴って、振子も揺れていることが確認できます。この時に振子の初期角度は0度に設定しておきます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/d0f5819b-ff31-0e2a-a862-1a0b281e5de5.png)

## 制御装置を付与する
* 制御対象を作れたので、次は状態フィードバック系を実装します。状態量として、台車の位置と速度、振子の角度と角速度をフィードバックします。
* Prismatic Jointの設定でSensingのPositionとVelocityを有効化します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/5054546e-2cbd-7aa6-0928-ddc7736fc742.png)

* 振子のRevolute Jointも同様にSensing設定でPositionとVelocityを設定します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/42b6f8a4-a1c7-cf1a-aaab-bb1d1d833d36.png)

* 振子のブロックから、角度と角速度を外部に出力します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/3e548908-b8bc-f5b1-037c-f04d6db7af2d.png)

* 状態フィードバック制御器をSimulinkで作ります。振子角度は、鉛直下向きが0radなので、鉛直上向きが0radとなるようにΠ[rad]だけオフセットしています(Revolute Jointの前後の座標変換をうまくやれば不要になる)。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/25bfc2e4-45cf-f73d-149d-0fb9ca92857d.png)

## 制御器を調整する
* シミュレーション結果やグラフを見ながらゲインを調整します。
1. 振子角度ゲイン
1. 振子角速度ゲイン
1. 台車位置のゲイン
1. 台車速度のゲイン
* 調整したゲインの一例と実際の制御の様子は以下の通りです
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/e8519697-c576-fe07-26e7-0023a8f838a8.png)

![CartPendforQiita3.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/8bb0693f-99a4-d219-2056-3d307d944438.gif)

※ソルバー設定は、チュートリアルを参考に以下の通りに設定
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203155/7db42a8d-e759-4b80-d2d7-ee8e7ff9a17f.png)


# 感想
Simscape Multibodyは座標系の概念に慣れ親しめれば、状態方程式や数式モデルレスでシミュレーション環境を構築できます。さらにシミュレーション結果の動画をすぐに得られる点も素晴らしいです。

当初は、参考文献[2]に記載の資料を読んで、「正統的技術系記事」を書きたいとは思いました。が、残念ながら、ジャンクフード群に分類される記事になってしまいました。

[MATLAB/Simulink Advent Calender 2022](https://qiita.com/advent-calendar/2022/matlab)の25日に、川田先生も倒立振子ネタを投稿されるようで、つよつよな方より先に投稿できて助かりました。

# 参考文献
[1][Simscape Multibody](https://jp.mathworks.com/help/sm/index.html?s_tid=CRUX_lftnav)
[2][テクノロジーマップ、技術カタログの在り方について](https://www.digital.go.jp/assets/contents/node/basic_page/field_ref_resources/b1ae707f-8dfa-4b40-8003-e8eec74f6a02/32230edd/20221003_meeting_technology_based_regulatory_reform_outline_05.pdf)