数式処理ソフトウェアで立式した運動方程式をSimulinkで実行する

この記事は制御工学アドベントカレンダー15日目の記事です．

制御工学で避けては通れないのが運動方程式の立式です．
プラントモデリングではもちろんのこと，現代制御に代表されるモデルベース制御では制御ロジックを構築するために運動方程式（状態方程式）が必要になります．

運動方程式の立式というと，ほぼラグランジュ法（オイラー＝ラグランジュ方程式）一択のような気がしています．

いまいちこのあたりの用語を掴みきれていないのですが，汎函数と呼ばれる“函数を変数とした函数”を変分（汎函数での微分）することで運動方程式を求めるといった流れです．

たとえば，$f(t) = at$は函数と呼ばれ，$\frac{df(t)}{dt}$という操作が微分と呼ばれることはご存知のとおりです．
それと似たような関係で，$I[f(t)] = Af(t)$を汎函数と呼び，$\frac{dI[f(t)]}{df(t)}$を変分と呼ぶようです．

物体の運動の場合，最も必要エネルギの少ない方向に物体が運動するという原理があるらしく（最小作用の原理？），以下の形式で与えられます．

$$I[ y ] = \int^{b}_{a}f( x, y, y' ) dx$$

特に運動の場合は，右辺$f(x,y,y')$をラグランジアン$L$というエネルギ量で与えるようです．
合っている自信はないですが，以下の式のようなイメージです．（時間微分なのでドットで書いています）

$$I[ x ] = \int^{b}_{a}L( t, x, \dot{x} ) dt$$

これはある時刻$a$から$b$までのエネルギの合計となります．
上の式が最小になるようにすれば良いということで，$I[x]$の変分（微分）がゼロとなるような$x$（停留値）を求めれば良いことになります．

この解（停留値）は決まっていて，これをオイラー＝ラグランジュ方程式と呼びます．
$$\frac{d}{dt}\left( \frac{\partial L}{\partial \dot{x}} \right) - \left( \frac{\partial L}{\partial x} \right) = 0$$

解の形式が決まっているのでシーケンシャルに解けて便利なのですが，非常に計算量が多いです．

前置きが長くなりましたがこれを数式処理ソフトウェアで楽にしようというのが今回のコンセプトです．

---

この記事では数式処理ソフトウェアとして以下を使用しています．

• wxMaxima（Windows版Maxima，無償）
• Symbolic Math Toolbox（MATLAB，有償）

恐らく運動方程式を求めたあとはMATLAB/SimulinkとかPyControlとかでシミュレーションすることになると思うのですが，結論としてはwxMaxima優位といった印象です．

---

今回，2重振り子を例として，Simulinkでのシミュレーションまでやってみようかと思います．

ラグランジュ法では物体の運動エネルギ，位置エネルギの関係を与えてあげることが必要になります．
逆に言うとエネルギまで出せてしまえばあとはソフトウェア任せにできます．

事前にやってみたところ，wxMaximaのインストールから結果確認まで45分程度で終わりました．

最適化函数を使って船を横移動させてみる

この記事は制御工学アドベントカレンダー9日目の記事です．
さいきん，仕事の方で船に関わることが多いので船の動きも面白いよ，と布教してみる．
船，とひとくちに言ってもタンカーなどの超大型船から客船，個人所有のクルーザー，漁船までいろいろ種類があります．
また，推進方法によっても船内にエンジンを積んだ船内機艇，エンジンごと取り外し可能な船外機艇，ジェット推進の艇などの種類があります．

ここではクルーザーなどに多く使われる下図のような船外機艇(Outboard motor)に焦点を当ててみます．

Mercury社のウェブページ
図の船は船外機が2つ取り付けられています（2機掛けと呼びます）が，最大で5つ取り付ける方もいるようです．
250馬力を5機掛け．パワーの面でも価格の面でも恐ろしいですね．

---

その船外機艇のなかでも面白いなあと思ったものがMercury社のSkyhookというJoystick操船システムです．
1分34秒あたりを見るとわかりやすいですが，着岸状態から船がほぼ真横に移動しています．

Mercury Joystick Piloting for Outboards
通常，船外機は何機掛けであろうと同じ方向を向かせて舵を切るのですが，この操船システムではそれぞれの船外機を独立させて反対方向に動かすことで横移動を実現しています．

このSkyhookを再現するため，最適化函数を使って遊んでみようというのがこの記事の趣旨です．
環境はMATLAB/Simulink R2015aです．

C/GMRESによる非線形モデル予測制御1

今回は珍しく2部構成で，本記事ではシステムの立式から評価函数の設定まで，次の記事では具体的な解き方について述べます．

第1部となる本記事ではC/GMRESの核となる部分には言及しないため，そちらについては2部以降を見ていただければと思います．

---

非線形モデル予測制御(Nonlinear Model Predictive Control; NMPC)とは，その名の通り”状態方程式が非線形なシステム”に対するモデル予測制御(Model Predictive Control; MPC)です．

最近流行りの制御則で，車の自動運転やロケットの再着陸システムにも使用されているそうです．

MyEnigmaさん
モデル予測制御(Model Predictive Control:MPC)の応用例

そもそも制約有りのMPCは非線形な制御系ですが，”非線形MPC”といった場合には“制御対象が非線形”な系を指すように思います．

線形のモデル予測制御については別途，記事に起こすつもりです．

通常，線形MPCでは，システムの表現に以下のような線形な状態方程式を使うのに対して，
$$\dot{x}(t) = Ax(t)+Bu(t)$$
非線形MPCでは状態方程式の形式を問わず使用でき，以下のような形式で表現されます．
$$\dot{x}(t)=f(x(t),u(t)) = \left[\begin{array}{c} x_{2}(t) \\ ax_{1}(t)+bx_{2}(t)u_{1}(t) \end{array}\right]$$

これはセミアクティブダンパの状態方程式ですが，式中に$x_{2}(t)u_{1}(t)$という“状態変数と操作量の積”が含まれています．少し考えるとわかりますが，線形の状態方程式に書き直すことはできません．

実際，このシステムを制御する際に単純に線形MPCを適用することはできません．
このような，非線形なシステムの制御を目的としたMPCが非線形MPCです．

非線形MPCの中にも種々ありますが，ここではC/GMRESによる非線形MPCについて取り上げます．

ただ，このシステムは減衰が操作量なのでダンパとして動作させること自体はPID等でも可能です．というよりずっとONにしておけばよいのです．しかし，できるだけ少ない操作量で最大の減衰を得よう，などと考え始めるとちゃんと制御する必要があります．

カルマンフィルタの基礎式を代数とベイズ定理から見る

カルマンフィルタの歴史ももう長いことと思いますが，各所で各人なりにまとめられているあたり，なかなか理解の難しいものなんだなあと感じています．
例に漏れず自分もカルマンフィルタの原理がいまいちピンとこなかったので自分なりにまとめてみました．

カルマンフィルタ(Kalman Filter)はフィルタと名付けられてはいますが，使われ方としてはフィルタよりも現代制御のオブザーバ（Observer; 状態推定器）に近いような気がします．
観測した情報から，状態ベクトルが真値に最も近くなるように推定するものです．

概念は以下のところでわかりやすく説明されています．

モノを作りたくなるブログ

シンプルなモデルとイラストでカルマンフィルタを直感的に理解してみる

ただ，数式として少し省略されているところもあったので，改めてまとめるに至りました．
基本的に，代数とベイズの定理のみから丁寧に導出していこうという感じです．

カルマンゲインの式の導出，予測値の算出，分散の更新は天下り的な部分も多いので，今回はこれらの導出に焦点を当ててみました．
統計的にみて，予測値と実プラントの誤差が最小となるように演繹的に求めていくとカルマンフィルタの式に行き着くあたり，統計学の応用先の代表例になっているだけのことはあるなあ，という感想です．

ここでは，状態空間と観測空間が同一の場合のみ扱います．

4ストロークエンジンの物理モデル

4ストロークのガソリンエンジンの物理モデリングをしてみました．
シミュレーションはMATLAB/Simulink，可視化ツールはProcessingです．

参考元はこのあたり．
双方ともに高知工科大学の磯村先生の文献です．

4 サイクルガソリンエンジンの動特性について

並列計算機によるエンジンのシミュレータ開発

---

以前にも倒立振子とかパーティクルフィルタでのグリッドマップマッチングとか2リンク平面マニピュレータとかのモデリングをやってきたのですが，そのどれとも色が違うなあ，ということで．

例によって，ラグランジュ法で運動方程式を立てます．
このあたりは参考文献の上の方が詳しいかな．

シンプルな外乱オブザーバ

外乱オブザーバ(Disturbance observer; DOB)はフィードバックシステムに外乱が加わるような状況で，比較的簡単に外乱を除去し安定化できる強力なツールですが，ウェブにはあまり記述が無いようだったので書いてみました．

• 外乱オブザーバとは？
• なぜ外乱オブザーバが必要なのか
• 基本原理
• 実際の外乱オブザーバ（ノミナルモデル，ローパスフィルタ）

---

そもそも外乱オブザーバとは？

“外乱”でググってみると，以下のような説明が得られます．

ある通信系に、所定の信号系以外から加わり妨害となる信号。 --- 大辞林

大体内容としては間違いないですが，今回扱うのは”通信系”に限らず，制御システムに対してはたらく，妨害となる信号くらいに留めておきます．

外乱オブザーバとは，このような外乱を除去するために，制御系のコントローラに組み込む補償器のひとつです．

---

なぜ外乱オブザーバが必要なのか

次のような1次遅れ系のフィードバックループを考えてみます．

ここで，$C(s)$はコントローラ，$G(s)$は制御対象です．また，目標値を$r$，出力を$y$で表します．
外乱というと，図の1, 2, 3番の箇所に加わることが想定されますが，本稿では1について記述します．

また，以下の理由で2, 3番を省いています．

2番の部分に入力される外乱は，フィードバックループを構成する主な目的のひとつのため，比較的簡単に除去（補償）できる．
3番については，システム自体は目標値通り動作している（外乱によって我々使用者が目標値どおりに動作しているか確認できないだけの状態）ため．

1番は，コントローラは制御対象を最適に動作させるための操作量を正常に生成できているが，それが外乱によって変えられ制御対象に入力されてしまっている状態です．

つまり，PIDなどのコントローラを適切に設計しても，設計通りの操作量が制御対象に入力されないということです．

今回はこの外乱を除去する外乱オブザーバについて述べていきます．

ここでは，"最適な操作量"は設計時に想定している操作量を指します．
最適レギュレータなどの最適とは根本的に異なる内容ですのであしからず．

ニューラルネットで雑な手書き文字認識

以前にもニューラルネットについて言及しましたが，かなり原初的な内容ばかりだったので今回は簡単な応用について．

C言語でニューラルネット（関数近似）

C言語でニューラルネット（2次元平面での分離）

といっても，修士の授業で作った内容そのままなので，いざ実装してみると動かないかも．

ニューラルネットの応用先として，真っ先に上がるのがパターン認識の分野でしょう．

とくに文字認識の分野ではすでにかなりの数，実用化されていると聞きます．

ここでは，最も単純な文字の認識について述べていきたいと思います．

具体的には以下の2点．

• ニューラルネット（誤差逆伝播法）の基本的な部分
• 特徴量（ニューラルネットへの入力）