　はじめまして、りょーつといいます。高専出身の大学院2年生です。研究の専門は力学や機構学で、Qiitaでは主に制御工学や数学に関する記事を書いています。今年度のキャチロボバトルコンテストで使用予定のモデル予測制御についてまとめていこうと思います。いちおうこちらの記事の続きになっています。

目次

1.はじめに
 2.モータのモデル予測制御
 3.モデル予測を用いた追従制御系
 4.おわりに

1. はじめに

1.1 本記事の目的

　今回の記事では、モデル予測制御を使ってモータの角速度を制御する方法について簡単にまとめていこうと思います。2年ほど前にモデル予測制御の記事を書いたのですが、内容がゴリゴリしすぎていて読みづらかったので、簡易版を作ろうとしています。モデル予測制御はちゃんと組めばとても性能がいい(経験上)ので、ロボコンでも普及するといいなぁと思い、この記事を書いています。なお、中の人の専門は機構学であり、制御工学は趣味でやっている程度なので、この記事の内容は独断と偏見の塊であることをあらかじめ明記しておきます。

　モデル予測制御(Model Predictive Control：MPC)は、名前のとおり、モデルを使って未来を予測し、いい感じに制御する手法です。非線形系に対してもそれっぽい性能を発揮するため、車の制御とかに取り入れられている(?)というのを聞いたことがあります。

　モデル予測制御の文献を漁ると、行列や数式がゴリゴリ出てくる上に大学の講義では扱われないため、手を出しづらいと思います。本記事ではできる限り数式を単純化して解説し、全国のロボコニストがモデル予測制御で遊べるようにすることを目的とします。

1.2. 時間の無い方へ

　とにかく早く実装したいという方のために本節には結果をずらっとまとめました。まず本記事では、モータの入力電圧$u(t)$と角速度$\omega(t)$の入出力関係が1次遅れ系で表現できることを仮定し、モータの運動方程式が以下の式で表されると考えます。

\frac{d\omega}{dt}
+
a \omega
=
b u
\tag{1}

$a$と$b$はモータ固有の定数です。手元にあるモータの$a$と$b$を調べる手法についてはこちらの記事を参照ください。(1)式で表されるモータを、制御周期$\Delta t$で目標角速度$\omega_{ref}$になるように制御するとき、以下式で計算できる入力$u$を与えるとよいです。

u
=
\frac{a}{b} \omega
+
\frac{1}{b \Delta t} (\omega_{ref} - \omega)
\tag{2}

(2)式は(1)式のモデルを用いた予測により、1制御周期分未来の状態を予測した結果計算される最適な入力です。理由や詳細が気になる方は2章以降を読んでみてください。

2. モータのモデル予測制御

　本章では、モデル予測制御を使ってモータを速度制御する際の、入力の決定方法についてまとめます。簡潔にいうと(2)式の導出を行います。

　まずはモデルから未来を予測するために、(1)式の離散化を行います。制御周期(=サンプリング周期)を$\Delta t$とすると、(1)式は以下のように離散化されます。

\omega_{n+1}
=
(1 - a \Delta t) \omega_n
+
b \Delta t \ u_n
\tag{3}

ここで$\omega_n$は時刻$t = n\Delta t$における角速度$\omega(n \Delta t)$を表しています。$\omega_{n+1}$、$u_n$も同じです。
　ここで注意するべきは(3)式の解釈です。(3)式は「ある時刻において、$\omega_n$の角速度で回転しているモータに対して$u_n$を制御入力として与えると、1制御周期後($\Delta t$ [s]後)に角速度が$\omega_{n+1}$になるよ」ということを意味しています。つまり、(3)式は1制御周期分の未来を予測する式であると解釈できます。

　1制御周期分の未来を予測することができたので、あとは目標角速度$\omega_{ref}$と$\omega_{n+1}$の差が最小になるような$u_n$を計算するだけです。評価関数$J$を以下のように定義します。

J
=
{(\omega_{ref} - \omega_{n+1})}^2
\tag{4}

(4)式を$u_n$で微分すると

\frac{dJ}{d u_n}
=
2(\omega_{ref} - \omega_{n+1})
\frac{d \omega_{n+1}}{d u_n}
=
-2(\omega_{ref} - \omega_{n+1}) b \Delta t
\tag{5}

になります。評価関数が極小になるとき、$\frac{dJ}{d u_n} = 0$になるため、最適入力$u_n$は$\omega_{ref} - \omega_{n+1} = 0$の解となります。(3)式に$\omega_{ref} - \omega_{n+1} = 0$を代入して実際に解いてみましょう。

\omega_{ref}
=
(1 - a \Delta t) \omega_n
+
b \Delta t \ u_n

b \Delta t \ u_n
=
a \Delta t \omega_n
+
(\omega_{ref} - \omega)

\therefore
\ u_n
=
\frac{a}{b} \omega_n
+
\frac{1}{b \Delta t} (\omega_{ref} - \omega_n)
\tag{6}

たしかに(2)式の結果が得られました。

　非常にシンプルな式で記述できたのは、未来を1制御周期分のみ予測したためであり、もっと先の未来(複数ステップの未来)を予測して制御をするのであれば、行列などを使ってゴリゴリする必要があります。
　ただし計算が複雑化するだけで解答方針自体は変わらないので、興味がある方はこちらの記事を参考にチャレンジしてみてほしいです！

3. モデル予測を用いた追従制御系

　本章では、外乱やモデル化誤差による偏差を抑えるために、モデル予測制御へ積分項を追加する方法についてまとめます。少しゴリゴリした内容になるので、2章の内容が物足りないという方に読んでみていただきたいです。

　使用するモータのモデルは(1)式と同じで、(3)式のように離散化します。これと並行して、目標角速度$\omega_{ref}$と角速度$\omega$の間の偏差を積分したパラメータ$e$を以下のように定義します。

e_n
=
\sum_{k = 0}^n
(\omega_{ref} - \omega_k)\Delta t
\tag{7}

(7)式を漸化式の形に展開すると

e_{n+1}
=
e_n
+
(\omega_{ref} - \omega_{n+1})\Delta t
\tag{8}

になります。モデル予測制御に積分項を追加するためには(4)式の評価関数$J$を以下のように修正すればよいと考えられます。

J
=
q_1 {(\omega_{ref} - \omega_{n+1})}^2
+
q_2 {e_{n+1}}^2
\tag{9}

なお、$q_1$、$q_2$は評価関数の各項の重み付け係数です。2章と同様に評価関数$J$の極小値を探しましょう。(9)式を入力$u_n$で微分すると

J
=
- 2 q_1 (\omega_{ref} - \omega_{n+1})\frac{d \omega_{n+1}}{d u_n}
+
2 q_2 \  e_{n+1} \frac{d e_{n+1}}{d u_n}

=
- 2 q_1 (\omega_{ref} - \omega_{n+1})\frac{d \omega_{n+1}}{d u_n}
- 2 q_2 \Delta t \  e_{n+1} \frac{d \omega_{n+1}}{d u_n}

=
- 2 
\big(
q_1 (\omega_{ref} - \omega_{n+1})
+
q_2 \Delta t \  e_{n+1} 
\big)
\frac{d \omega_{n+1}}{d u_n}
\tag{10}

　評価関数が極小になるとき、$\frac{dJ}{d u_n} = 0$になるため、これを解きたいと思います。(10)式に$\frac{dJ}{d u_n} = 0$を代入して整理すると

- 2 
\big(
q_1 (\omega_{ref} - \omega_{n+1})
+
q_2 \Delta t \  e_{n+1} 
\big)
\frac{d \omega_{n+1}}{d u_n}
=
0

q_1 (\omega_{ref} - \omega_{n+1})
+
q_2 \Delta t \  e_{n+1} 
=
0

q_1 (\omega_{ref} - \omega_{n+1})
+
q_2 \Delta t \  \big( e_n + (\omega_{ref} - \omega_{n+1})\Delta t  \big)
=
0

q_2 \Delta t \ e_n
+
\big(
q_1 + q_2 {\Delta t}^2 \big)(\omega_{ref} - \omega_{n+1})
=
0

\therefore
\omega_{n+1}
=
\omega_{ref}
+
\frac{q_2 \Delta t}{q_1 + q_2 {\Delta t}^2}e_n
\tag{11}

最後に(3)式と(11)式を連立して$\omega_{n+1}$を消去し、(6)式と同様に$u_n$について解けば以下のようにして最適入力が得られます。

(1-a \Delta t)\omega_n + b \Delta t \ u_n
=
\omega_{ref}
+
\frac{q_2 \Delta t}{q_1 + q_2 {\Delta t}^2}e_n

\therefore
u_n
=
\frac{a}{b} \omega_n
+
\frac{1}{b \Delta t} (\omega_{ref} - \omega_n)
+
\frac{q_2 \Delta t}{b \Delta t(q_1 + q_2 {\Delta t}^2)}e_n
\tag{12}

ここで(12)式の右辺第3項の係数をゲイン$k_I$と置くことで、モデル予測制御器に積分項が追加されたような入力になります。

\frac{q_2 \Delta t}{b \Delta t(q_1 + q_2 {\Delta t}^2)}
=
\frac{\Delta t}{b \Delta t(\frac{q_1}{q_2} + {\Delta t}^2)}
=
k_I
\tag{13}

u_n
=
\frac{a}{b} \omega_n
+
\frac{1}{b \Delta t} (\omega_{ref} - \omega_n)
+
k_I \ e_n
\tag{14}

　2章と同様、(14)式のように入力$u_n$がシンプルな式で導出されたのは未来を1制御周期分のみ予測したためであり、もっと先の未来(複数ステップの未来)を予測して制御をするのであれば、行列などを使ってゴリゴリする必要があります。
　ただし計算が複雑化するだけで解答方針自体は変わりません。興味がある方はぜひチャレンジしてみてほしいです！

4. おわりに

　本稿では初めての方にもノリを知っていただくため、1ステップのみ予測するモデル予測制御の理論的な部分についてまとめました。実際に(14)式に沿ってモータを回してみたところ、1ステップのみの予測でもかなりうまく制御できたので、Arduinoのプログラムコードや実験結果についても近いうちにまとめたいなぁと考えています。

　おそらく中の人のチームはキャチロボで(14)式の制御を使用予定なので、大会当日はそこら辺にも注目してみていただけると嬉しいです。最後まで読んでいただき、ありがとうございました！