MATLABによるマルチエージェント合意制御 差分法・微分法による(下準備)

こんにちは。初投稿です。
いつも興味深く皆さんの記事を拝見しているのですが、今回個人的にやっていたことが「ここに書いてもいいかもしれない」と思ったので、公開してみました。
Markdownにそもそも慣れてないんですが、どうぞ寛大なお心でご覧いただければと思います。

……と思って書いていたら重くて長くて仕方ないので、事前準備編/コード編で分けて公開しようと思います。ここは下準備編です。
コード編はこちら→[差分法・微分法のコード](https://qiita.com/Lgarwn/items/5a723bad795295e15269 "プログラムのコード")

---
# 0. まとめ
1. MATLABでマルチエージェントシステムの合意制御を前進差分法による近似的解法・微分方程式による解法の両方で行ってみた。どっちでも出来た。よかった。
2. 実行時間からすると前進差分法の方がよさそう。結果もサンプリング周期を小さくとればさほど変わんないし。ただ厳密に求めるなら微分方程式の方でもよい。正直「微分方程式の扱いになれる」くらいしかないかもしれん。

# 1. 背景とねらい
## 1.1. 背景
近年、マルチエージェントシステム(以下MAS)という分野で多く研究が行われています。理論的検討をしたり強化学習を入れたり、また多種多様な応用分野で適用されていたり、という感じですね。ドローンの制御なんかにも使われていそうです。
このMASですが、基本的な問題/制御課題の1つに「合意制御」というものがあります。詳細は文献をご覧いただきたいんですが、簡単に言うと「全てのエージェントが何らかの形で定義された目標に対し、自身の状態をそれと同一にする」ような感じです。
漠然としちゃったので具体例ですが、2次元平面上に10台のエージェントがいることを考えると、一番シンプルな合意制御は「各エージェントの重心位置に全エージェントが重なることを目標」として、制御を行います。

## 1.2. ねらい
この合意制御なんですが、(僕が思うに)一般的には微分方程式の形でエージェントの制御則を示すことがほとんどです。例えばエージェントの位置が速度入力によって制御される場合、位置を$x$、入力を$u$とすれば
```math
\dot{x}=u
```
という感じですね。ここで、また私の事情で申し訳ないんですが、このような微分方程式は全て前進差分で近似して、それを数値的に解くことで解を求めていました。今の例で言えば移動した軌跡や入力などの値は、すべて近似によって求めていました。
ですが、MATLABです。なんと、Symbolic Math Toolboxを使えば微分・積分は思いのまま。なんと素晴らしいソフトウェアなんだろうと思ったわけです。
……宣伝というわけじゃありませんよ？(個人的にはかなり長くお世話になっているので感謝の念しかない)

閑話休題ですが、そんな **「いつも行っている差分法による解法」** と、 **「微分方程式としての解法」** を比較してみたいな、と思って書いてみました。これがねらいです。


## 1.3. 諸注意
:::note warn
さっさと本題に入れという話ですが、最初に注意事項を申し上げます。
+ 井の中の蛙大海を知らずだと思います。もっと効率のいいコードなどをご存じの場合はコメントにてご指摘いただければ幸いです(なお筆者は心が弱いので優しくご教示いただけますと大変ありがたいです)。
+ 上に関連して、自己流に用語を使っていたりするかもしれません。これも同様にご指摘ください。
+ Markdownも初めてならQiitaも書くのは初めてだったりします。どうぞ温かい目で見ていただければこれも幸いです。
:::
それでは本題に入りましょうか。

---

# 2. 準備
本格的にプログラムを記述する前に、問題設定をここで定義しておきます。今回はあくまで差分法・微分方程式の比較検討のためなので、前述の合意制御を行うことを考えます。

## 2.1. シミュレーション空間・時間
まず、本稿ではシミュレーション空間は2次元平面であるとします。直交座標系ですね。壁は存在しませんが、フィールドの上下左右端に到達したら反対側に出てくるトーラス形状であるとします。
またシミュレーションは全部で$t_{end}$秒だけ行われるものとし、差分法におけるサンプリング周期は$t_s$であるとします。
この前提の下で話を進めていきます。

## 2.2. エージェントの基礎設計
今回は$n$台のエージェントが2次元平面上にランダムな初期位置で存在するものとします。このすべてのエージェントの集合を$G$とし、$|G|=n$とします。
またエージェント$i(i\in G)$の$x,y$座標をそれぞれ$x_i,y_i$とし、与える入力を$u_{xi},u_{yi}$としましょう。さらに$x_i,y_i,u_{xi},u_{yi}$を全エージェント分縦にまとめた列ベクトル$X=[x_1,x_2,...,x_n]^T,Y=[y_1,y_2,...,y_n]^T,U_X=[u_{x1},u_{x2},...,u_{xn}]^T,U_Y=[u_{y1},u_{y2},...,u_{yn}]^T$を定義します。

## 2.3. MASの基礎
ここで、$n$台のエージェントは簡単のため、全てのエージェントが自分以外のエージェントを認識できるものとします。つまり、あるエージェント$i(i\in G)$はエージェント$j (j\neq i,j\in G)$を認識できる、というわけですね。
このとき、MASでは「近傍」という概念をよく導入します。これは「自分以外のエージェントのうち、自分が認識できるもの」の集合であり、エージェント$i$の近傍は$N_i$で定義します。また近傍に属するエージェントの数を$|N_i|$で表します(定義より$|N_i|\leq|G|$ですね、だからどうだっていう話ですが)。
この近傍という概念を用いると、いくつかMASに特有の(他の分野でも用いられてそうですが)行列をいくつか定義できます。詳細は各文献をご参照いただくとして、それらを簡単に述べると下記のようになります。

##### 次数行列
あるエージェントについて、「他のエージェントから認識されている数」を指します。
$n\times n$の行列であり、対角要素に「認識されている数」が入り、非対角要素は0です。例えば今回は、すべてのエージェントは自分以外のエージェントを認識できるので、対角要素に$n-1$、非対角要素に0が入った行列になりますね。
文字としては$D$を使って、今回の例で言うと(1-1)式のようになります。
~~~math
D=
\begin{bmatrix}
n-1      & \cdots & 0 \\
\vdots & \ddots & \vdots \\
0      & \cdots & n-1
\end{bmatrix}
\tag{1-1}
~~~

##### 接続行列 
あるエージェントについて、「他のエージェントを認識しているかどうか」を指します。隣接行列とも言いますね。
$n\times n$の行列であり、要素は以下のようになります。
~~~math
a_{ij} = \left\{
\begin{array}{ll}
0 & (i=j) \\
0 & (i\neq j \quad and \quad j\notin N_i) \\
1 & (i\neq j \quad and \quad j\in N_i)
\end{array}
\right.
~~~
具体的に言うと、対角要素は0、$(i,j)$要素$(i\neq j)$には「エージェント$i$が$j$を認識していれば1、そうでなければ0」が入ります。
例えば今回は、すべてのエージェントは自分以外のエージェントを認識できるので、対角要素に0、非対角要素にすべて1が入った行列になります。
文字としては$A$を使って、今回の例で言うと(1-2)式のようになります。
~~~math
A=
\begin{bmatrix}
0 & \cdots & 1 \\
\vdots & \ddots & \vdots \\
1 & \cdots & 0
\end{bmatrix}
\tag{1-2}
~~~

##### グラフラプラシアン
上記の接続行列・次数行列から定義される行列です。具体的には、グラフラプラシアンを$L$として
$$
L=D-A
$$
で定義されます。本稿では、上記(1-1)、(1-2)式から(1-3)式のように求められます。
~~~math
\begin{align}
\qquad L&=D-A \\

\Leftrightarrow \quad &=
\begin{bmatrix}
n-1      & \cdots & 0 \\
\vdots & \ddots & \vdots \\
0      & \cdots & n-1
\end{bmatrix}
-
\begin{bmatrix}
0 & \cdots & 1 \\
\vdots & \ddots & \vdots \\
1 & \cdots & 0
\end{bmatrix} \\

\Leftrightarrow \quad &=
\begin{bmatrix}
n-1 & \cdots & -1 & \cdots & -1  \\
-1 & n-1 & -1 & \cdots & -1 \\ \\
\vdots & \vdots & \vdots & \vdots & \vdots \\
-1 & \cdots & -1 &\cdots & n-1
\end{bmatrix} \tag{1-3} \\
\end{align}

~~~
※これが実は合意制御では大切な役割をします(今回は触れませんが)。

## 2.4. エージェントの入力
今、全エージェントは以下に示す入力によって移動を行うものとします。
~~~ math
\begin{align}
\dot{X}=U_X \tag{1-1} \\
\dot{Y}=U_Y \tag{1-2}
\end{align}
~~~
この入力ですが、合意制御の場合は(1-3)、(1-4)式とします。
~~~ math
\begin{align}
U_X = -LX \tag{1-3} \\
U_Y = -LY \tag{1-4}
\end{align}
~~~
ここで$L$は$n\times n$の行列、$X,Y$は$n\times 1$の列ベクトルであることに注意してください。
つまり、(1-1)・(1-3)式および(1-2)・(1-4)式より、次の(1-5)・(1-6)式を得ます。
~~~ math
\begin{align}
\dot{X} = -LX \tag{1-5} \\
\dot{Y} = -LY \tag{1-6}
\end{align}
~~~

ここで、$\dot{X}$および$\dot{Y}$の$i$番目の成分に着目して前進差分を行います。なお、以降は簡単のため$x$方向成分のみを考え、着目する成分を$x_i$とします(書くのが面倒になったともいう)。
前進差分を行うと、
~~~ math
\dot{x_i}[k] ~= \frac{x_i[k+1]-x_i[k]}{T_s} \tag{1-7}
~~~
を得ます。なお$k$はステップ数であり、$k$ステップ目からサンプリング周期$T_s$だけ経つと次の$k+1$ステップ目に進むとしました。
このとき、$\dot{x_i}[k]=u_{xi}[k]$ですから(1-7)式は以下のようになりますね。
~~~ math
u_{xi}[k] = \frac{x_i[k+1]-x_i[k]}{T_s} \tag{1-8}
~~~
ここで$u_{xi}[k]$は$k$ステップ目の入力です。なお、サンプリング周期が十分に小さいとして近似しました。

では(1-8)式を変形して、位置$x_i$についての差分式にしましょう。右辺の分母分子に$T_s$をかけて整理すれば、(1-9)式を得ます。
~~~ math
x_i[k+1]-x_i[k] = u_{xi}[k] \times T_s\tag{1-9}
~~~
ここまでを$y$軸方向成分についても行って、(1-10)式を得ます。
~~~ math
y_i[k+1]-y_i[k] = u_{yi}[k] \times T_s\tag{1-10}
~~~

ここでは(数学的な式変形の厳密性に自信がなかったこともあって)$i$番目の成分=エージェントに着目しましたが、これを全台数分で行います。プログラムではベクトル・行列処理です。

前フリが長くなってしまいましたがここまでで下準備は完了です。次ではコードを示します。

---

……というところですみません、冒頭に書いた通り以降は次の記事に示します。
URLはこちら→[差分法・微分法のコード](https://qiita.com/Lgarwn/items/5a723bad795295e15269 "プログラムのコード")


# 参考文献
桜間 一徳, マルチエージェントシステムの制御 : III 合意制御(1), システム／制御／情報, 2013, 57 巻, 9 号, p. 386-396