宇宙ステーションの巨大なトラス構造を組み立てる場面を想像してください。1台のロボットアームが黙々とモジュールを運搬し、取り付け位置に固定し、次のモジュールを取りに行く — この作業を数百回繰り返す。1台では途方もない時間がかかります。では、10台のロボットが同時に、互いに衝突せず、かつ作業の重複もなく組み立てを進められたらどうでしょうか。作業時間は劇的に短縮されるはずです。
しかし、10台のロボットに「さあ、協力して組み立てなさい」と命じるだけではうまくいきません。どのロボットがどのモジュールを運ぶのか。ロボット同士がぶつからないようにどうフォーメーションを組むのか。通信が途切れたらどうするのか。1台が故障したら残りでどうリカバーするのか。これらの問いに体系的に答えるのがマルチロボット協調（Multi-Robot Coordination）の理論です。
マルチロボット協調を理解すると、以下のような応用に直結します。
- 軌道上組立: 大型宇宙構造物（宇宙太陽光発電衛星、大型アンテナ、宇宙ステーション）を複数ロボットで効率的に組み立てる
- デブリ除去フリート: 複数の除去衛星が協調して、優先度の高いデブリを同時並行で除去する
- コンステレーション保守: 数百〜数千機の衛星群の中で、複数のサービサー衛星が分担して燃料補給や修理を行う
- 惑星探査: 複数のローバーやドローンが協調して広域探査を行う（火星のマルチローバーミッション構想）
本記事の内容
- 集中制御 vs 分散制御の基本思想と数学的定式化
- グラフ理論による通信トポロジーの表現
- コンセンサスアルゴリズム（平均合意）の理論と収束証明
- タスク割当問題（ハンガリアン法、オークションベースアルゴリズム）
- フォーメーション制御（リーダーフォロワー法、仮想構造体法）
- 通信制約下での協調戦略
- 軌道上組立・デブリ除去フリートへの応用
- Pythonによるコンセンサス・フォーメーション制御のシミュレーション
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- スペースデブリ除去のロボティクス — 非協力物体の捕獲戦略と宇宙デブリ問題
- PID制御の理論と実装 — フィードバック制御の基礎
- グラフ理論の基礎 — 通信トポロジーの記述に必要なグラフの概念
集中制御と分散制御 — 2つのアーキテクチャ
まず身近な例で考える
複数のロボットを協調させるには、大きく分けて2つの方法があります。日常のアナロジーで考えてみましょう。
1つ目は集中制御です。サッカーの試合で、監督がリアルタイムに全選手に無線で指示を出す場面を想像してください。「A選手は左に走れ、B選手はパスを出せ、C選手はゴール前に待機しろ」— 監督が全体を見て最適な指示を出します。しかし、無線が途切れたら？ 監督が判断を間違えたら？ 全員が止まってしまいます。
2つ目は分散制御です。渡り鳥の群れを思い浮かべてください。リーダーの鳥が全体に指令を出しているわけではありません。各鳥は「隣の鳥との距離を保つ」「前の鳥の方向に合わせる」「群れの中心に向かう」という単純なローカルルールだけに従って飛んでいます。それだけで、群れ全体として美しいV字フォーメーションが自然に形成されるのです。
集中制御の数学的定式化
集中制御では、中央の制御局がすべてのロボットの状態を集約し、最適な指令を計算して各ロボットに配信します。
$N$ 台のロボットが存在するとき、ロボット $i$ の状態を $\bm{x}_i \in \mathbb{R}^n$ とします。全体の状態ベクトルは各ロボットの状態を縦に並べたものです。
$$ \bm{X} = \begin{bmatrix} \bm{x}_1 \\ \bm{x}_2 \\ \vdots \\ \bm{x}_N \end{bmatrix} \in \mathbb{R}^{nN} $$
中央制御局は全体の状態 $\bm{X}$ を入力として、全ロボットへの制御入力を一括で計算します。
$$ \bm{U} = \begin{bmatrix} \bm{u}_1 \\ \bm{u}_2 \\ \vdots \\ \bm{u}_N \end{bmatrix} = f(\bm{X}) $$
ここで $f$ は中央の制御則（例えば最適制御問題の解）です。この方式のメリットは、全体の情報を使って大域的に最適な解を求められる点にあります。
しかし、集中制御には本質的な問題があります。
| 問題 | 具体的な影響 |
|---|---|
| 通信ボトルネック | 全ロボットが中央と通信するため、ロボット数に比例して通信量が増大 |
| 単一障害点 | 中央制御局が故障すると全体が停止 |
| スケーラビリティ | $nN$ 次元の最適化問題を毎ステップ解く必要があり、計算量が $O(N^3)$ 以上に膨れ上がる |
| 通信遅延 | 宇宙空間では光速遅延があり、リアルタイム集中制御が困難 |
分散制御の数学的定式化
分散制御では、各ロボットは自分自身と近傍のロボットの情報のみに基づいて制御入力を決定します。ロボット $i$ の制御則は次のように書けます。
$$ \bm{u}_i = g_i(\bm{x}_i, \{\bm{x}_j : j \in \mathcal{N}_i\}) $$
ここで $\mathcal{N}_i$ はロボット $i$ の近傍集合、つまりロボット $i$ と直接通信できるロボットの集合です。各ロボットは局所的な情報だけで行動を決めるため、中央制御局は不要です。
分散制御の利点を整理すると次のようになります。
| 利点 | 説明 |
|---|---|
| ロバスト性 | 1台が故障しても他のロボットは自律的に動作を継続できる |
| スケーラビリティ | ロボット数を増やしても各ロボットの計算量は近傍の数にのみ依存する |
| 通信効率 | 局所的な通信で済むため、全体の通信量が抑えられる |
| 適応性 | 環境の変化にローカルに対応でき、全体の再計画が不要 |
一方で、分散制御には「各ロボットが局所情報しか持たないため、大域的な最適性が保証されない」という根本的なトレードオフがあります。このトレードオフに対処するための理論的基盤が、次に解説するグラフ理論と合意アルゴリズムです。
ここまでで、集中制御と分散制御の基本的な枠組みがわかりました。分散制御を実現するためには、「誰が誰と通信できるのか」を数学的に記述する道具が必要です。それがグラフ理論です。
通信トポロジーとグラフ理論
なぜグラフが必要なのか
分散制御では、各ロボットは「近傍」のロボットとだけ情報をやり取りします。しかし、この「近傍」の構造は固定的ではありません。ロボットが移動すると通信範囲が変わり、近傍関係が動的に変化します。この通信構造を厳密に記述するために、グラフ理論の枠組みを使います。
通信グラフの定義
$N$ 台のロボットの通信関係をグラフ $\mathcal{G} = (\mathcal{V}, \mathcal{E})$ で表現します。
- 頂点集合 $\mathcal{V} = \{1, 2, \dots, N\}$: 各頂点がロボット1台に対応
- 辺集合 $\mathcal{E} \subseteq \mathcal{V} \times \mathcal{V}$: ロボット $i$ と $j$ が通信可能なら辺 $(i, j) \in \mathcal{E}$
通信が双方向の場合、グラフは無向グラフになります。片方向の場合は有向グラフです。宇宙ロボットの通信は通常、同じ通信機器を使った双方向通信なので、ここでは無向グラフを基本とします。
隣接行列とラプラシアン行列
通信グラフの構造を行列として表現する方法を導入します。これがマルチロボット制御の数学的基盤となります。
隣接行列 $\bm{A} \in \mathbb{R}^{N \times N}$ は、ロボット間の通信関係を表す行列です。
$$ a_{ij} = \begin{cases} 1 & \text{if } (i, j) \in \mathcal{E} \\ 0 & \text{otherwise} \end{cases} $$
無向グラフの場合、$\bm{A}$ は対称行列（$a_{ij} = a_{ji}$）になります。
次数行列 $\bm{D} \in \mathbb{R}^{N \times N}$ は、各ロボットの通信相手の数を対角に並べた行列です。
$$ \bm{D} = \text{diag}(d_1, d_2, \dots, d_N), \quad d_i = \sum_{j=1}^{N} a_{ij} $$
$d_i$ はロボット $i$ の次数（degree）で、通信相手の数を表します。
そして最も重要な行列がグラフラプラシアン $\bm{L}$ です。次数行列と隣接行列の差として定義されます。
$$ \bm{L} = \bm{D} – \bm{A} $$
ラプラシアン行列の成分を具体的に書くと次のようになります。
$$ l_{ij} = \begin{cases} d_i & \text{if } i = j \\ -1 & \text{if } (i, j) \in \mathcal{E} \\ 0 & \text{otherwise} \end{cases} $$
ラプラシアン行列の重要な性質
グラフラプラシアンはマルチロボット制御を支える中心的な数学的対象です。以下の性質を押さえておきましょう。
性質1: 各行の和がゼロ
$$ \sum_{j=1}^{N} l_{ij} = d_i – \sum_{j=1}^{N} a_{ij} = d_i – d_i = 0 $$
これは $\bm{L} \bm{1} = \bm{0}$（$\bm{1}$ は全成分が1のベクトル）と書けます。つまり $\bm{1}$ はラプラシアンの固有値0に対応する固有ベクトルです。
性質2: 半正定値
任意のベクトル $\bm{z} \in \mathbb{R}^N$ に対して、次のような二次形式の変形ができます。
$$ \bm{z}^\top \bm{L} \bm{z} = \bm{z}^\top (\bm{D} – \bm{A}) \bm{z} $$
隣接行列と次数行列の定義を代入して展開すると、
$$ \bm{z}^\top \bm{L} \bm{z} = \sum_{i=1}^{N} d_i z_i^2 – \sum_{i=1}^{N}\sum_{j=1}^{N} a_{ij} z_i z_j $$
$a_{ij} = a_{ji}$ かつ $(i,j) \in \mathcal{E}$ の辺について整理すると、
$$ \bm{z}^\top \bm{L} \bm{z} = \frac{1}{2}\sum_{(i,j) \in \mathcal{E}} (z_i – z_j)^2 \geq 0 $$
この式は極めて直感的です。ラプラシアンの二次形式は、辺で結ばれたロボット間の状態の差の二乗和に等しいのです。差が大きいほど値が大きくなるので、ラプラシアンは「ネットワーク全体の不一致度」を測る尺度として機能します。
性質3: 連結性との関係
ラプラシアンの固有値を小さい順に $0 = \lambda_1 \leq \lambda_2 \leq \dots \leq \lambda_N$ と並べたとき、$\lambda_2 > 0$ であることとグラフが連結であること（すべてのロボット間にパスが存在すること）は同値です。この $\lambda_2$ は代数的連結度（algebraic connectivity） またはフィードラー値と呼ばれ、ネットワークの連結の「強さ」を定量化します。
ラプラシアン行列の性質が、次に解説するコンセンサスアルゴリズムの収束保証に直接結びつきます。
コンセンサスアルゴリズム — 分散的な合意形成
全体の意思を局所通信だけで決める
分散制御において最も基本的な問題は、各ロボットがそれぞれ異なる初期値（位置、推定値、観測値など）を持っているとき、全体の平均値に合意することです。これを平均合意（average consensus） と呼びます。
日常的な例で考えてみましょう。10人が円形に座っていて、それぞれが手元の温度計で気温を測りました。しかし測定にはばらつきがあります。全員の平均を知りたいのですが、全員が互いに話せるわけではなく、隣の人としか会話できません。それでも、「隣の人の値と自分の値を平均して更新する」操作を繰り返すだけで、全員が真の平均に収束していく — これがコンセンサスアルゴリズムの核心です。
連続時間コンセンサスプロトコル
ロボット $i$ のスカラー状態 $x_i(t)$ を考えます。各ロボットが近傍との差に比例して状態を更新するプロトコルを記述します。
$$ \dot{x}_i(t) = \sum_{j \in \mathcal{N}_i} a_{ij} \left( x_j(t) – x_i(t) \right) $$
この式の意味は明快です。自分の値 $x_i$ よりも近傍 $x_j$ の値が大きければ $x_i$ を増やし、小さければ減らす。つまり、隣のロボットの値に近づくように自分の値を調整しているだけです。
全ロボットの状態をベクトル $\bm{x}(t) = [x_1(t), x_2(t), \dots, x_N(t)]^\top$ としてまとめると、このプロトコルは行列形式で簡潔に書けます。
$$ \dot{\bm{x}}(t) = -\bm{L}\bm{x}(t) $$
ここで $\bm{L}$ はグラフラプラシアンです。この方程式は線形常微分方程式であり、解は行列指数関数で表されます。
$$ \bm{x}(t) = e^{-\bm{L}t}\bm{x}(0) $$
収束の証明
コンセンサスプロトコルが全ロボットを平均値に収束させることを証明しましょう。
まず、ラプラシアンの固有値分解を利用します。$\bm{L}$ は対称半正定値なので、固有値 $0 = \lambda_1 < \lambda_2 \leq \dots \leq \lambda_N$（グラフが連結の場合）と正規直交固有ベクトル $\bm{v}_1, \bm{v}_2, \dots, \bm{v}_N$ が存在します。特に $\lambda_1 = 0$ に対応する固有ベクトルは $\bm{v}_1 = \frac{1}{\sqrt{N}}\bm{1}$ です。
初期状態 $\bm{x}(0)$ を固有ベクトルで展開すると、
$$ \bm{x}(0) = \sum_{k=1}^{N} c_k \bm{v}_k, \quad c_k = \bm{v}_k^\top \bm{x}(0) $$
行列指数関数 $e^{-\bm{L}t}$ を各固有値に適用すると、時刻 $t$ の状態は次のようになります。
$$ \bm{x}(t) = \sum_{k=1}^{N} c_k e^{-\lambda_k t} \bm{v}_k $$
$t \to \infty$ で $\lambda_k > 0$（$k \geq 2$）の成分は指数的に減衰し、$\lambda_1 = 0$ の成分のみが残ります。
$$ \lim_{t \to \infty} \bm{x}(t) = c_1 \bm{v}_1 = \left(\bm{v}_1^\top \bm{x}(0)\right) \bm{v}_1 $$
$\bm{v}_1 = \frac{1}{\sqrt{N}}\bm{1}$ を代入して計算すると、
$$ c_1 = \bm{v}_1^\top \bm{x}(0) = \frac{1}{\sqrt{N}} \sum_{i=1}^{N} x_i(0) $$
したがって、収束先は各成分が等しくなります。
$$ \lim_{t \to \infty} x_i(t) = \frac{1}{\sqrt{N}} \cdot \frac{1}{\sqrt{N}} \sum_{j=1}^{N} x_j(0) = \frac{1}{N}\sum_{j=1}^{N} x_j(0) = \bar{x} $$
つまり、全ロボットの状態が初期値の平均 $\bar{x}$ に収束することが示されました。しかも収束速度は第2固有値 $\lambda_2$ で決まります。$\lambda_2$ が大きいほど（ネットワークの連結度が高いほど）、合意に到達するまでの時間が短くなります。
離散時間コンセンサスプロトコル
実際のディジタル実装では、連続時間ではなく離散的に更新を行います。離散時間コンセンサスプロトコルは次のように書けます。
$$ x_i[k+1] = x_i[k] + \epsilon \sum_{j \in \mathcal{N}_i} a_{ij}\left(x_j[k] – x_i[k]\right) $$
ここで $\epsilon > 0$ はステップサイズで、収束のためには $\epsilon < \frac{1}{\max_i d_i}$ を満たす必要があります。ベクトル表記では次のようになります。
$$ \bm{x}[k+1] = (\bm{I} – \epsilon \bm{L})\bm{x}[k] = \bm{P}\bm{x}[k] $$
$\bm{P} = \bm{I} – \epsilon \bm{L}$ はペロン行列と呼ばれ、各行の和が1の確率的行列（doubly stochastic matrix）になります。
Pythonでコンセンサスをシミュレーションする
理論の確認として、6台のロボットがリング状の通信グラフ上で平均合意に到達する過程をシミュレーションしてみましょう。
import numpy as np
import matplotlib.pyplot as plt
# ロボット数とパラメータ
N = 6
dt = 0.01 # 時間刻み
T = 10.0 # シミュレーション時間
steps = int(T / dt)
# リング状通信グラフの隣接行列
A = np.zeros((N, N))
for i in range(N):
A[i, (i + 1) % N] = 1
A[i, (i - 1) % N] = 1
# ラプラシアン行列
D = np.diag(A.sum(axis=1))
L = D - A
# ラプラシアンの固有値
eigvals = np.sort(np.linalg.eigvalsh(L))
print(f"ラプラシアンの固有値: {np.round(eigvals, 4)}")
print(f"代数的連結度 λ₂ = {eigvals[1]:.4f}")
# 初期状態（各ロボットの値をランダムに設定）
np.random.seed(42)
x0 = np.random.uniform(0, 10, N)
x_avg = np.mean(x0)
print(f"初期値: {np.round(x0, 2)}")
print(f"平均値: {x_avg:.4f}")
# コンセンサスのシミュレーション（オイラー法）
x_history = np.zeros((steps + 1, N))
x_history[0] = x0.copy()
for k in range(steps):
x_history[k + 1] = x_history[k] - dt * L @ x_history[k]
# 可視化
time = np.linspace(0, T, steps + 1)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 左: 各ロボットの状態推移
for i in range(N):
axes[0].plot(time, x_history[:, i], label=f"Robot {i+1}")
axes[0].axhline(y=x_avg, color='k', linestyle='--', alpha=0.7, label=f"Average = {x_avg:.2f}")
axes[0].set_xlabel("Time [s]")
axes[0].set_ylabel("State value")
axes[0].set_title("Consensus Protocol on Ring Graph")
axes[0].legend(loc="upper right", fontsize=8)
axes[0].grid(True, alpha=0.3)
# 右: 不一致度の推移
disagreement = np.array([np.sum((x_history[k] - x_avg)**2) for k in range(steps + 1)])
axes[1].semilogy(time, disagreement)
axes[1].set_xlabel("Time [s]")
axes[1].set_ylabel("Disagreement (sum of squared errors)")
axes[1].set_title("Convergence of Disagreement")
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
左のグラフから、6台のロボットがそれぞれ異なる初期値（0〜10の範囲）から出発し、時間とともに全員が平均値（黒の破線）に収束していく様子がわかります。リング状グラフでは各ロボットは隣の2台としか通信していないにもかかわらず、局所的な情報交換だけで大域的な合意が達成されています。
右のグラフは不一致度（各ロボットの状態と平均値の差の二乗和）の時間推移を対数スケールで描いたものです。直線的に減少しており、これは不一致度が指数的に減衰していること — すなわち収束速度が $e^{-\lambda_2 t}$ に従っていることを示しています。リングのラプラシアンの第2固有値 $\lambda_2$ は比較的小さいため、収束は緩やかです。完全グラフや星型グラフなど、連結度の高いトポロジーを使えば収束は速くなりますが、通信コストとのトレードオフがあります。
コンセンサスアルゴリズムは「全員が同じ値に合意する」問題を解きました。しかし、マルチロボット協調では「全員が同じことをする」のではなく、「各ロボットに異なるタスクを適切に振り分ける」ことが必要です。次に、タスク割当の問題を扱います。
タスク割当 — 誰が何をやるか
タスク割当問題とは
軌道上に10個のモジュールがあり、10台のロボットがそれぞれ1つずつモジュールを運搬するとします。ロボット $i$ がモジュール $j$ を運搬するコスト（消費燃料、所要時間など）を $c_{ij}$ とすると、全体のコストを最小化する割当を見つけたいという問題です。
これは数学的には線形割当問題（Linear Assignment Problem, LAP） に帰着されます。$N$ 台のロボットと $N$ 個のタスクの割当を決定変数 $x_{ij} \in \{0, 1\}$ で表現すると、次のように定式化できます。
$$ \min_{\bm{X}} \sum_{i=1}^{N}\sum_{j=1}^{N} c_{ij} x_{ij} $$
制約条件として、各ロボットはちょうど1つのタスクを担当し、各タスクはちょうど1台のロボットに割り当てられます。
$$ \sum_{j=1}^{N} x_{ij} = 1 \quad \forall i \quad (\text{各ロボットは1タスクのみ}) $$
$$ \sum_{i=1}^{N} x_{ij} = 1 \quad \forall j \quad (\text{各タスクは1ロボットのみ}) $$
$$ x_{ij} \in \{0, 1\} $$
この問題はNP困難に見えるかもしれませんが、実は多項式時間で厳密に解けることが知られています。その代表的アルゴリズムがハンガリアン法です。
ハンガリアン法の直感とアルゴリズム
ハンガリアン法は1955年にハロルド・クーンが発表したアルゴリズムで、$O(N^3)$ の計算量で最適割当を求めます。名前の由来は、このアルゴリズムの理論的基盤となったハンガリーの数学者エグレヴァーリとケーニッヒの定理に敬意を表したものです。
ハンガリアン法の直感を、具体的な数値例で説明しましょう。3台のロボットと3つのタスクを考えます。コスト行列が次のように与えられているとします。
$$ \bm{C} = \begin{bmatrix} 4 & 2 & 8 \\ 6 & 4 & 2 \\ 3 & 7 & 5 \end{bmatrix} $$
ハンガリアン法の手順は次の通りです。
ステップ1: 行縮約 — 各行の最小値を引きます。
各行の最小値は $[2, 2, 3]$ です。これを各行から引くと、
$$ \bm{C}’ = \begin{bmatrix} 2 & 0 & 6 \\ 4 & 2 & 0 \\ 0 & 4 & 2 \end{bmatrix} $$
ステップ2: 列縮約 — 各列の最小値を引きます。
各列の最小値は $[0, 0, 0]$ なので、この例では変化しません。
ステップ3: 最小本数の直線でゼロを覆う — 行と列の線を引いて、すべてのゼロを覆うのに必要な最小線の数を求めます。もしその数が $N$ に等しければ、ゼロ位置から最適割当が得られます。等しくなければ追加のステップで行列を変形し、ステップ3に戻ります。
この例では3本の線でゼロを覆えるので、最適割当は次のようになります。ロボット1にタスク2（コスト2）、ロボット2にタスク3（コスト2）、ロボット3にタスク1（コスト3）で、合計コストは $2 + 2 + 3 = 7$ です。
Pythonによるハンガリアン法の実装
SciPyの linear_sum_assignment を使って、軌道上タスク割当をシミュレーションしてみましょう。
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import linear_sum_assignment
# 5台のロボットと5つのタスク（軌道上モジュール）の位置
np.random.seed(123)
N = 5
robot_positions = np.random.uniform(-10, 10, (N, 2))
task_positions = np.random.uniform(-10, 10, (N, 2))
# コスト行列（ユークリッド距離）
cost_matrix = np.zeros((N, N))
for i in range(N):
for j in range(N):
cost_matrix[i, j] = np.linalg.norm(robot_positions[i] - task_positions[j])
print("コスト行列:")
print(np.round(cost_matrix, 2))
# ハンガリアン法で最適割当を求める
row_ind, col_ind = linear_sum_assignment(cost_matrix)
optimal_cost = cost_matrix[row_ind, col_ind].sum()
print(f"\n最適割当: {list(zip(row_ind, col_ind))}")
print(f"最適総コスト: {optimal_cost:.2f}")
# 比較: ランダム割当の平均コスト
n_random = 10000
random_costs = []
for _ in range(n_random):
perm = np.random.permutation(N)
random_costs.append(sum(cost_matrix[i, perm[i]] for i in range(N)))
print(f"ランダム割当の平均コスト: {np.mean(random_costs):.2f}")
# 可視化
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
# ロボットとタスクをプロット
ax.scatter(robot_positions[:, 0], robot_positions[:, 1],
c='cyan', s=150, marker='s', edgecolors='white', zorder=5, label='Robots')
ax.scatter(task_positions[:, 0], task_positions[:, 1],
c='orange', s=150, marker='^', edgecolors='white', zorder=5, label='Tasks')
# 最適割当の接続線を描画
for i, j in zip(row_ind, col_ind):
ax.plot([robot_positions[i, 0], task_positions[j, 0]],
[robot_positions[i, 1], task_positions[j, 1]],
'g-', linewidth=2, alpha=0.7)
ax.annotate(f"R{i+1}", robot_positions[i], fontsize=10,
ha='center', va='bottom', color='white')
ax.annotate(f"T{j+1}", task_positions[j], fontsize=10,
ha='center', va='bottom', color='white')
ax.set_xlabel("x [m]")
ax.set_ylabel("y [m]")
ax.set_title(f"Optimal Task Assignment (Hungarian Method)\nTotal cost = {optimal_cost:.2f}")
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
plt.tight_layout()
plt.show()
グラフを見ると、各ロボット（四角）が最も近いタスク（三角）に割り当てられている場合と、全体の最適化のためにあえて遠いタスクに割り当てられている場合があることがわかります。これが「局所最適」と「大域最適」の違いです。各ロボットが個別に最も近いタスクを選ぶ貪欲法では、割当が衝突したり全体コストが最適にならなかったりしますが、ハンガリアン法は全体を俯瞰して最適な組み合わせを見つけます。また、ランダム割当の平均コストと比較すると、最適割当が大幅にコストを削減していることが数値的に確認できます。
オークションベースアルゴリズム — 分散的なタスク割当
ハンガリアン法は集中的なアルゴリズムであり、全ロボットのコスト行列が1箇所に集約されている必要があります。軌道上で通信制約がある場合には、分散的にタスク割当を行う手法が求められます。
Bertsekas（1988）が提案したオークションアルゴリズムは、経済学の競り（オークション）のメタファーに基づく分散割当手法です。
基本的な仕組みは次の通りです。
- 入札フェーズ: 各ロボットは未割当のタスクの中から最も価値の高い（コストの低い）タスクを選び、そのタスクに「入札」する。入札額は、最良タスクと次善タスクの価値の差に小さな増分 $\epsilon$ を加えたものです。
$$ b_i = (v_{i,j^*} – p_{j^*}) – (v_{i,j’} – p_{j’}) + \epsilon $$
ここで $v_{i,j}$ はロボット $i$ にとってのタスク $j$ の価値、$p_j$ はタスク $j$ の現在価格、$j^*$ は最良タスク、$j’$ は次善タスクです。
-
割当フェーズ: 各タスクは最高入札者に割り当てられ、その入札額だけ価格が上昇する。
-
繰り返し: すべてのロボットにタスクが割り当てられるまで1〜2を繰り返す。
オークションアルゴリズムの利点は、各ロボットが局所的な情報（自分のコストと現在の価格情報）だけで入札を決定できる点にあります。完全な分散実装では、価格情報は近傍間のメッセージパッシングで伝搬させます。計算量は $O(N^2 / \epsilon)$ であり、$\epsilon$ を小さくするほど最適解に近づきますが収束は遅くなります。
ここまでで、各ロボットにタスクを割り振る方法がわかりました。しかし、タスクを実行するためにはロボットが適切な位置関係を保ちながら移動する必要があります。次は、複数ロボットが隊形を維持しながら移動する「フォーメーション制御」を扱います。
フォーメーション制御 — 隊形を維持しながら移動する
なぜフォーメーションが重要か
軌道上で大型構造物を搬送する場面を考えましょう。3台のロボットが巨大なトラス部材を3点で支えながら移動する場合、ロボット間の相対位置が変わってしまうと部材に不均一な力がかかり、破損やロボットの制御不安定につながります。また、デブリ除去フリートが一定の間隔を保ちながらスキャンする場合も、フォーメーションの維持が不可欠です。
フォーメーション制御には主に3つのアプローチがあります。
| 手法 | 特徴 | 適用場面 |
|---|---|---|
| リーダーフォロワー | 1台のリーダーに他が追従 | 経路が明確な搬送タスク |
| 仮想構造体 | 仮想的な剛体を定義しその上に配置 | 厳密な形状維持が必要な組立作業 |
| ビヘイビアベース | 複数の行動ルールを重み付け合成 | 環境適応が必要な探査タスク |
リーダーフォロワー法
リーダーフォロワー法は最も直感的なフォーメーション制御手法です。1台のリーダーロボットが目標軌道に沿って移動し、残りのフォロワーはリーダーとの相対位置関係を維持するように制御されます。
リーダーの位置を $\bm{p}_L(t)$ とし、フォロワー $i$ の目標位置をリーダーからのオフセットで定義します。
$$ \bm{p}_i^{\text{des}}(t) = \bm{p}_L(t) + \bm{d}_i $$
ここで $\bm{d}_i$ はフォロワー $i$ のリーダーに対するオフセットベクトルです。例えば、正三角形フォーメーションなら $\bm{d}_1 = [r, 0]^\top$, $\bm{d}_2 = [-r/2, r\sqrt{3}/2]^\top$ のように設定します。
フォロワー $i$ の制御則は、目標位置と現在位置の偏差を用いた比例微分（PD）制御で設計できます。
$$ \bm{u}_i = K_p (\bm{p}_i^{\text{des}} – \bm{p}_i) + K_d (\dot{\bm{p}}_i^{\text{des}} – \dot{\bm{p}}_i) $$
$K_p, K_d > 0$ はそれぞれ比例ゲインと微分ゲインです。
リーダーフォロワー法はシンプルで実装が容易ですが、本質的な弱点があります。リーダーが故障するとフォーメーション全体が崩壊するという単一障害点の問題です。また、フォロワーはリーダーの動きに「反応して」追従するため、リーダーの急な方向転換に対して遅れが生じ、フォーメーションが一時的に崩れます。
仮想構造体法
仮想構造体法（Virtual Structure）は、より厳密なフォーメーション維持を実現するアプローチです。フォーメーション全体を1つの仮想的な剛体として扱い、その剛体の位置・姿勢・速度を制御します。各ロボットは仮想剛体上の固定点に対応づけられます。
仮想構造体の状態を、重心位置 $\bm{p}_c(t)$ と姿勢角 $\theta_c(t)$ で表します。ロボット $i$ の仮想構造体上の基準位置を $\bm{r}_i^b$（ボディフレーム座標）とすると、目標位置は回転行列を用いて次のように計算されます。
$$ \bm{p}_i^{\text{des}}(t) = \bm{p}_c(t) + \bm{R}(\theta_c(t)) \bm{r}_i^b $$
ここで、2次元の回転行列は次の通りです。
$$ \bm{R}(\theta) = \begin{bmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{bmatrix} $$
仮想構造体の動力学は、目標軌道に追従するように設計します。
$$ \ddot{\bm{p}}_c = K_p^c (\bm{p}_c^{\text{ref}} – \bm{p}_c) + K_d^c (\dot{\bm{p}}_c^{\text{ref}} – \dot{\bm{p}}_c) $$
各ロボットは仮想構造体上の目標位置に追従するPD制御を実行します。仮想構造体法の強みは、フォーメーションの幾何的形状を厳密に維持できることです。軌道上で精密な組立作業を行う場合には、この厳密性が不可欠です。
一方で弱点もあります。仮想構造体の状態を全ロボットが共有する必要があるため、完全な分散実装が難しい点です。通信遅延があると仮想構造体の状態の同期にずれが生じ、各ロボットの目標位置にも不整合が発生します。
コンセンサスベースのフォーメーション制御
リーダーフォロワー法と仮想構造体法のそれぞれの弱点を補う第三のアプローチとして、先ほど学んだコンセンサスアルゴリズムを拡張したフォーメーション制御があります。
各ロボットの目標は「近傍とのオフセットを所望の値に一致させる」ことです。ロボット $i$ の位置を $\bm{p}_i$、フォーメーション上での目標オフセットを $\bm{d}_{ij}$ とすると、制御則は次のようになります。
$$ \dot{\bm{p}}_i = \sum_{j \in \mathcal{N}_i} \left[(\bm{p}_j – \bm{p}_i) – \bm{d}_{ij}\right] $$
変数 $\bm{q}_i = \bm{p}_i – \bm{d}_i$（ただし $\bm{d}_i$ は共通参照点からのオフセット）に変換すると、これは標準的なコンセンサスプロトコルに帰着されます。
$$ \dot{\bm{q}}_i = \sum_{j \in \mathcal{N}_i} (\bm{q}_j – \bm{q}_i) = -\sum_{j=1}^{N} l_{ij} \bm{q}_j $$
したがって、先に証明したコンセンサスの収束定理がそのまま適用でき、$\bm{q}_i \to \bar{\bm{q}}$（全員の平均）に収束します。これをもとの変数に戻すと $\bm{p}_i \to \bar{\bm{q}} + \bm{d}_i$ となり、各ロボットが所望のオフセット位置に収束することが保証されます。この方法では特定のリーダーが不要であり、完全な分散実装が可能です。
Pythonによるフォーメーション制御のシミュレーション
リーダーフォロワー法と仮想構造体法の両方を実装し、円形軌道に沿ったフォーメーション移動をシミュレーションします。
import numpy as np
import matplotlib.pyplot as plt
# --- パラメータ ---
N_followers = 4
dt = 0.01
T = 20.0
steps = int(T / dt)
# フォーメーションオフセット（正五角形の頂点、リーダーが先頭）
angles = np.linspace(0, 2 * np.pi, N_followers + 1, endpoint=False)
formation_radius = 2.0
offsets = np.column_stack([formation_radius * np.cos(angles[1:]),
formation_radius * np.sin(angles[1:])])
# --- リーダーの軌道（円形） ---
omega = 0.3 # 角速度
R_orbit = 8.0 # 軌道半径
def leader_trajectory(t):
px = R_orbit * np.cos(omega * t)
py = R_orbit * np.sin(omega * t)
vx = -R_orbit * omega * np.sin(omega * t)
vy = R_orbit * omega * np.cos(omega * t)
return np.array([px, py]), np.array([vx, vy])
# --- リーダーフォロワー法 ---
Kp, Kd = 2.0, 3.0
# 初期位置（ランダムにずらす）
np.random.seed(7)
positions_lf = np.random.uniform(-3, 3, (N_followers, 2))
velocities_lf = np.zeros((N_followers, 2))
trajectory_lf = np.zeros((steps + 1, N_followers, 2))
leader_traj = np.zeros((steps + 1, 2))
trajectory_lf[0] = positions_lf.copy()
for k in range(steps):
t = k * dt
p_L, v_L = leader_trajectory(t)
leader_traj[k] = p_L
for i in range(N_followers):
p_des = p_L + offsets[i]
v_des = v_L
acc = Kp * (p_des - positions_lf[i]) + Kd * (v_des - velocities_lf[i])
velocities_lf[i] += acc * dt
positions_lf[i] += velocities_lf[i] * dt
trajectory_lf[k + 1] = positions_lf.copy()
leader_traj[steps] = leader_trajectory(T)[0]
# --- 仮想構造体法 ---
positions_vs = np.random.uniform(-3, 3, (N_followers, 2))
velocities_vs = np.zeros((N_followers, 2))
p_c = np.array([0.0, 0.0]) # 仮想構造体の重心
v_c = np.array([0.0, 0.0])
theta_c = 0.0
Kp_c, Kd_c = 1.5, 2.5
Kp_f, Kd_f = 3.0, 4.0
trajectory_vs = np.zeros((steps + 1, N_followers, 2))
vs_center_traj = np.zeros((steps + 1, 2))
trajectory_vs[0] = positions_vs.copy()
vs_center_traj[0] = p_c.copy()
for k in range(steps):
t = k * dt
p_ref, v_ref = leader_trajectory(t)
# 仮想構造体の重心を目標軌道に追従
acc_c = Kp_c * (p_ref - p_c) + Kd_c * (v_ref - v_c)
v_c += acc_c * dt
p_c += v_c * dt
theta_c = np.arctan2(v_c[1], v_c[0])
R_mat = np.array([[np.cos(theta_c), -np.sin(theta_c)],
[np.sin(theta_c), np.cos(theta_c)]])
for i in range(N_followers):
p_des = p_c + R_mat @ offsets[i]
v_des = v_c
acc = Kp_f * (p_des - positions_vs[i]) + Kd_f * (v_des - velocities_vs[i])
velocities_vs[i] += acc * dt
positions_vs[i] += velocities_vs[i] * dt
trajectory_vs[k + 1] = positions_vs.copy()
vs_center_traj[k + 1] = p_c.copy()
# --- 可視化 ---
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# リーダーフォロワー法
ax = axes[0]
ax.plot(leader_traj[:, 0], leader_traj[:, 1], 'w--', alpha=0.4, label='Leader path')
for i in range(N_followers):
ax.plot(trajectory_lf[:, i, 0], trajectory_lf[:, i, 1], alpha=0.5, linewidth=0.8)
# 最終位置
t_final = steps
p_L_final = leader_traj[t_final]
ax.plot(p_L_final[0], p_L_final[1], 'r*', markersize=15, label='Leader (final)')
for i in range(N_followers):
ax.plot(trajectory_lf[t_final, i, 0], trajectory_lf[t_final, i, 1],
'o', markersize=8, label=f'Follower {i+1}')
ax.set_xlabel("x [m]")
ax.set_ylabel("y [m]")
ax.set_title("Leader-Follower Formation Control")
ax.legend(fontsize=7, loc='upper right')
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
# 仮想構造体法
ax = axes[1]
ax.plot(vs_center_traj[:, 0], vs_center_traj[:, 1], 'w--', alpha=0.4, label='VS center path')
for i in range(N_followers):
ax.plot(trajectory_vs[:, i, 0], trajectory_vs[:, i, 1], alpha=0.5, linewidth=0.8)
for i in range(N_followers):
ax.plot(trajectory_vs[t_final, i, 0], trajectory_vs[t_final, i, 1],
'o', markersize=8, label=f'Robot {i+1}')
ax.set_xlabel("x [m]")
ax.set_ylabel("y [m]")
ax.set_title("Virtual Structure Formation Control")
ax.legend(fontsize=7, loc='upper right')
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
左のリーダーフォロワー法の結果を見ると、フォロワーがリーダーの円軌道に追従しつつ、フォーメーション形状を維持している様子がわかります。初期のランダム配置からフォーメーションが形成されるまでの過渡期では軌跡が乱れていますが、PD制御のゲインが十分大きいため、速やかに目標位置に収束しています。
右の仮想構造体法では、フォーメーション全体が仮想的な剛体として回転しながら移動しています。円軌道の曲率に応じてフォーメーションの向きが変わり、各ロボットが構造体上の固定点に追従する様子が見て取れます。リーダーフォロワー法と比較すると、フォーメーションの回転が考慮されているため、曲線軌道上でより自然な隊形配置が実現されています。
両手法とも所望のフォーメーションに収束していますが、実際の軌道上応用では通信遅延やセンサノイズの影響が加わります。次のセクションでは、通信に制約がある場合の協調戦略について考えます。
通信制約下での協調 — 途切れてもなんとかする
軌道上通信の特殊性
軌道上でのマルチロボット通信には、地上とは根本的に異なる制約が存在します。
距離と指向性: ロボット間の距離が数km以上に及ぶ場合、無指向性アンテナでの通信には高い送信出力が必要です。出力を抑えるために指向性アンテナを使うと、通信を維持するための姿勢制御が新たな制約になります。
間欠的な通信: 軌道力学的な制約（地球の影に入る、他の構造物に遮蔽される）により、通信が周期的に途絶するのは日常的なことです。ISSのロボットアームですら、TDRS経由の通信が使えない「ブラックアウト」時間帯があります。
帯域制限: 宇宙空間での通信帯域は貴重な資源であり、各ロボットが大量のデータ（高解像度画像、詳細な状態ベクトルなど）を常時やり取りすることは現実的ではありません。
イベント駆動型通信
通信帯域を節約するために、各ロボットが「必要なときだけ」通信するイベント駆動型（event-triggered）通信が研究されています。ロボット $i$ は自分の状態 $\bm{x}_i(t)$ を直近に送信した値 $\hat{\bm{x}}_i$ と比較し、誤差がしきい値を超えたときのみ送信します。
$$ \|\bm{x}_i(t) – \hat{\bm{x}}_i\| > \sigma_i(t) $$
ここで $\sigma_i(t)$ は動的なしきい値関数で、例えば近傍との不一致度に比例して設定します。
$$ \sigma_i(t) = c_0 + c_1 \sum_{j \in \mathcal{N}_i} \|\hat{\bm{x}}_j – \hat{\bm{x}}_i\| $$
$c_0 > 0$ は最低限のしきい値、$c_1 > 0$ は調整パラメータです。近傍との不一致が大きい（合意が遠い）ときはしきい値が緩くなり、不一致が小さい（合意に近い）ときはしきい値が厳しくなります。この設計により、合意の初期段階では頻繁に通信し、合意に近づくにつれて通信頻度が自動的に下がります。
理論的に重要なのは、このイベント駆動型プロトコルでも「ゼノ挙動」（有限時間に無限回の通信イベントが発生すること）が起きない条件が保証されている点です。$c_0 > 0$ を正に保つことで、通信イベント間に必ず正の最小時間間隔が存在することが示されます。
通信トポロジーの切り替え
ロボットの移動や通信途絶により、通信グラフが時間とともに変化する状況をスイッチングトポロジーと呼びます。時刻 $t$ における通信グラフを $\mathcal{G}(t)$、対応するラプラシアンを $\bm{L}(t)$ と書くと、コンセンサスプロトコルは次のようになります。
$$ \dot{\bm{x}}(t) = -\bm{L}(t)\bm{x}(t) $$
スイッチングトポロジー下でコンセンサスが達成される十分条件として、次の結果が知られています。
定理（Jadbabaie et al., 2003）: 任意の長さ $T > 0$ の時間区間 $[t, t+T]$ において、区間内に出現するグラフの合併（union graph）が連結であるならば、コンセンサスは漸近的に達成される。
つまり、各時刻で通信グラフが連結である必要はなく、十分長い時間で見たときに全体が「つながっている」ことが保証されればよいのです。これは軌道上ロボットにとって極めて好都合な結果です。地球の影で通信が一時的に途切れても、周回ごとに通信が回復するならば、合意は時間をかけて必ず達成されます。
Pythonで通信途絶下のコンセンサスをシミュレーションする
通信トポロジーが周期的に切り替わる場合のコンセンサスを確認しましょう。
import numpy as np
import matplotlib.pyplot as plt
N = 6
dt = 0.01
T = 30.0
steps = int(T / dt)
# 2つの通信グラフ：連結なリングと分断されたグラフ
# グラフ1: リング（連結）
A1 = np.zeros((N, N))
for i in range(N):
A1[i, (i + 1) % N] = 1
A1[i, (i - 1) % N] = 1
L1 = np.diag(A1.sum(axis=1)) - A1
# グラフ2: 2つの独立クリーク（非連結）
A2 = np.zeros((N, N))
for i in range(3):
for j in range(3):
if i != j:
A2[i, j] = 1
for i in range(3, 6):
for j in range(3, 6):
if i != j:
A2[i, j] = 1
L2 = np.diag(A2.sum(axis=1)) - A2
# スイッチング周期
T_switch = 5.0 # 5秒ごとに切り替え
# 初期状態
np.random.seed(42)
x0 = np.random.uniform(0, 10, N)
x_avg = np.mean(x0)
# シミュレーション
x_sw = np.zeros((steps + 1, N))
x_sw[0] = x0.copy()
topology_log = np.zeros(steps + 1)
for k in range(steps):
t = k * dt
# 周期的に切り替え: 前半はリング、後半は分断
cycle_phase = (t % (2 * T_switch)) / (2 * T_switch)
if cycle_phase < 0.5:
L_current = L1
topology_log[k] = 1 # 連結
else:
L_current = L2
topology_log[k] = 0 # 非連結
x_sw[k + 1] = x_sw[k] - dt * L_current @ x_sw[k]
topology_log[steps] = topology_log[steps - 1]
# 比較: 常に分断されたグラフ
x_disc = np.zeros((steps + 1, N))
x_disc[0] = x0.copy()
for k in range(steps):
x_disc[k + 1] = x_disc[k] - dt * L2 @ x_disc[k]
# 可視化
time = np.linspace(0, T, steps + 1)
fig, axes = plt.subplots(2, 1, figsize=(12, 8), gridspec_kw={'height_ratios': [3, 1]})
# 上: 状態推移
ax = axes[0]
colors = plt.cm.tab10(np.linspace(0, 1, N))
for i in range(N):
ax.plot(time, x_sw[:, i], color=colors[i], linewidth=1.5,
label=f"Robot {i+1} (switching)")
ax.plot(time, x_disc[:, i], color=colors[i], linewidth=1.0,
linestyle=':', alpha=0.5)
ax.axhline(y=x_avg, color='white', linestyle='--', alpha=0.5,
label=f"True average = {x_avg:.2f}")
ax.set_ylabel("State value")
ax.set_title("Consensus under Switching Topology")
ax.legend(fontsize=7, loc='upper right', ncol=2)
ax.grid(True, alpha=0.3)
# 下: トポロジー状態
ax2 = axes[1]
ax2.fill_between(time, topology_log, alpha=0.3, color='cyan', step='mid')
ax2.set_xlabel("Time [s]")
ax2.set_ylabel("Connected?")
ax2.set_yticks([0, 1])
ax2.set_yticklabels(["Disconnected", "Connected"])
ax2.set_title("Communication Topology State")
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
上のグラフ（実線）は、スイッチングトポロジー下での各ロボットの状態推移です。下のグラフが通信状態を示しており、周期的に連結（リング）と非連結（2クリーク）が切り替わっています。非連結期間中は、2つのクリーク内でそれぞれ局所的な合意が進みますが、クリーク間では情報が伝わりません。しかし、連結期間に入ると2つのクリーク間の情報がリングを通じて伝搬し、全体の合意が進みます。
点線は常に非連結（2クリーク）の場合の結果です。全体の平均に収束せず、2つの異なる値に分離していることがわかります。スイッチングトポロジーの理論（合併グラフの連結性条件）が正しく機能していることが、この対比から明確に確認できます。
通信制約下でも協調が可能であることを確認しました。ここからは、これらの理論を軌道上の具体的な応用場面に当てはめていきましょう。
軌道上組立への応用 — 大型構造物をロボットで作る
なぜ軌道上で組み立てるのか
直径数十メートルの大型アンテナや宇宙太陽光発電衛星は、地上で完成品としてロケットに搭載することは物理的に不可能です。ロケットのフェアリング（先端カバー）の直径は最大でも5〜6m程度であり、大型構造物は折り畳むか、分割して打ち上げ、軌道上で組み立てる必要があります。
国際宇宙ステーション（ISS）はまさに軌道上組立の代表例で、40回以上のミッションで組み立てられました。しかしISSの組立はEVA（宇宙飛行士の船外活動）とカナダアーム2を組み合わせた手作業が主であり、膨大な時間とコストがかかりました。次世代の大型構造物を効率的に組み立てるには、複数のロボットが自律的に協調する能力が不可欠です。
組立タスクの分解
軌道上組立をマルチロボット協調の枠組みで捉えると、次のようなサブ問題に分解できます。
1. タスクグラフの構築: 組立作業には順序依存性があります。例えば「トラス部材Aを取り付けてからパネルBを展開する」という順序は守らなければなりません。この依存関係は有向非巡回グラフ（DAG） で表現できます。ノードが作業タスク、辺が依存関係です。
2. タスク割当: DAGの中で並列実行可能なタスクを特定し、利用可能なロボットに割り当てます。先ほど解説したハンガリアン法やオークションアルゴリズムが適用できます。
3. フォーメーション制御: 大型部材の搬送では、複数ロボットが部材の複数点を把持してフォーメーションを維持しながら移動します。仮想構造体法が適しています。
4. 衝突回避: 複数ロボットが同時に作業する場合、ロボット同士の衝突を避ける必要があります。人工ポテンシャル場法では、ロボット間に仮想的な斥力を設定します。ロボット $i$ と $j$ の間の斥力ポテンシャルは次のように定義できます。
$$ U_{ij}^{\text{rep}} = \begin{cases} \frac{1}{2} k_{\text{rep}} \left(\frac{1}{\|\bm{p}_i – \bm{p}_j\|} – \frac{1}{d_0}\right)^2 & \text{if } \|\bm{p}_i – \bm{p}_j\| < d_0 \\ 0 & \text{otherwise} \end{cases} $$
ここで $d_0$ は斥力が発生する最小距離、$k_{\text{rep}}$ は斥力の強さです。この斥力から導かれる力は次のようになります。
$$ \bm{F}_{ij}^{\text{rep}} = -\nabla_{\bm{p}_i} U_{ij}^{\text{rep}} $$
引力（目標位置への引力）と斥力（他ロボットとの衝突回避）を合成することで、各ロボットは安全にタスク位置へ移動できます。
組立シミュレーション
簡易的な軌道上組立シナリオとして、4台のロボットが4つのモジュールを供給ステーションから組立位置に搬送する問題をシミュレーションします。
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import linear_sum_assignment
# パラメータ
N_robots = 4
dt = 0.02
T = 15.0
steps = int(T / dt)
# 供給ステーション位置（モジュールの初期位置）
supply_positions = np.array([[-8, -5], [-8, -2], [-8, 2], [-8, 5]], dtype=float)
# 組立位置（正方形構造の四隅）
assembly_positions = np.array([[5, -3], [5, 3], [9, -3], [9, 3]], dtype=float)
# ロボットの初期位置
np.random.seed(55)
robot_positions = np.array([[-3, -4], [-2, 0], [-3, 3], [-1, -1]], dtype=float)
# タスク割当（ハンガリアン法で供給→組立の割当を決定）
cost = np.zeros((N_robots, N_robots))
for i in range(N_robots):
for j in range(N_robots):
# ロボットiが供給位置jに行き、組立位置jに運ぶ総距離
cost[i, j] = (np.linalg.norm(robot_positions[i] - supply_positions[j])
+ np.linalg.norm(supply_positions[j] - assembly_positions[j]))
row_ind, col_ind = linear_sum_assignment(cost)
# 各ロボットのウェイポイント: 供給位置 → 組立位置
waypoints = []
for i in range(N_robots):
j = col_ind[i]
waypoints.append([supply_positions[j], assembly_positions[j]])
# シミュレーション（ポテンシャル場ベースのナビゲーション）
k_att = 1.0 # 引力ゲイン
k_rep = 5.0 # 斥力ゲイン
d_safe = 1.5 # 安全距離
v_max = 2.0 # 最大速度
positions = robot_positions.copy()
phase = np.zeros(N_robots, dtype=int) # 0: 供給へ移動, 1: 組立へ移動, 2: 完了
trajectory = np.zeros((steps + 1, N_robots, 2))
trajectory[0] = positions.copy()
for k in range(steps):
for i in range(N_robots):
if phase[i] >= 2:
continue
# 現在の目標
target = waypoints[i][phase[i]]
# 引力
diff = target - positions[i]
dist = np.linalg.norm(diff)
if dist < 0.3:
phase[i] += 1
continue
f_att = k_att * diff
# 斥力（他ロボットとの衝突回避）
f_rep = np.zeros(2)
for j in range(N_robots):
if i == j:
continue
dij = positions[i] - positions[j]
dist_ij = np.linalg.norm(dij)
if dist_ij < d_safe and dist_ij > 0.01:
f_rep += k_rep * (1.0 / dist_ij - 1.0 / d_safe) / (dist_ij ** 2) * (dij / dist_ij)
# 合力
f_total = f_att + f_rep
speed = np.linalg.norm(f_total)
if speed > v_max:
f_total = f_total / speed * v_max
positions[i] += f_total * dt
trajectory[k + 1] = positions.copy()
# 可視化
fig, ax = plt.subplots(1, 1, figsize=(12, 7))
# 供給ステーションと組立位置
ax.scatter(supply_positions[:, 0], supply_positions[:, 1],
c='yellow', s=200, marker='D', edgecolors='white', zorder=5, label='Supply')
ax.scatter(assembly_positions[:, 0], assembly_positions[:, 1],
c='lime', s=200, marker='s', edgecolors='white', zorder=5, label='Assembly')
# 軌跡と最終位置
colors = ['#00BFFF', '#FF6B6B', '#FFD93D', '#6BCB77']
for i in range(N_robots):
ax.plot(trajectory[:, i, 0], trajectory[:, i, 1],
color=colors[i], alpha=0.5, linewidth=1.0)
ax.plot(trajectory[0, i, 0], trajectory[0, i, 1],
'o', color=colors[i], markersize=10, label=f'Robot {i+1} start')
ax.plot(trajectory[-1, i, 0], trajectory[-1, i, 1],
'*', color=colors[i], markersize=15)
# 割当関係を矢印で表示
for i in range(N_robots):
j = col_ind[i]
ax.annotate("", xy=supply_positions[j], xytext=robot_positions[i],
arrowprops=dict(arrowstyle='->', color=colors[i], alpha=0.3, lw=1.5))
ax.annotate("", xy=assembly_positions[j], xytext=supply_positions[j],
arrowprops=dict(arrowstyle='->', color=colors[i], alpha=0.3, lw=1.5))
ax.set_xlabel("x [m]")
ax.set_ylabel("y [m]")
ax.set_title("Multi-Robot Orbital Assembly Simulation")
ax.legend(fontsize=7, loc='lower right')
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
plt.tight_layout()
plt.show()
シミュレーション結果を見ると、4台のロボットがそれぞれ異なる初期位置から出発し、割り当てられた供給ステーション（黄色のダイヤモンド）でモジュールを回収した後、組立位置（緑色の四角）に搬送している様子がわかります。途中、ロボット同士が接近する場面では斥力ポテンシャルが働いて自動的に回避軌道を取っており、衝突なく全モジュールの搬送が完了しています。
注目すべきは、ハンガリアン法による割当が単純な「最近接」割当とは異なる結果を生むことがある点です。個別に見れば遠回りに見える割当でも、全体のコストが最小になるよう最適化されています。
軌道上組立に必要な技術要素が揃いました。次は、もう一つの重要な応用先であるデブリ除去フリートについて見ていきましょう。
デブリ除去フリート — 複数衛星による軌道清掃
フリート運用の必要性
先述のように、ケスラーシンドロームの進行を食い止めるには年間5〜10個の大型デブリを能動的に除去する必要があります。しかし、現在計画されているClearSpace-1のような単体ミッションでは、1回のミッションで1個のデブリを除去するのに数年の開発と数百億円のコストがかかります。これでは到底間に合いません。
解決策の一つがデブリ除去フリート — 複数の小型サービサー衛星を軌道上に配備し、協調的にデブリを除去するコンセプトです。フリート運用では以下の問題をマルチロボット協調の枠組みで解く必要があります。
デブリ優先度の動的割当
軌道上には数万個の追跡可能デブリが存在しますが、すべてを同時に除去することは不可能です。衝突確率や質量に基づいて優先度を計算し、限られたサービサーを効率的に割り当てる問題は、動的なタスク割当問題として定式化できます。
デブリ $j$ の優先度スコアを次のように定義します。
$$ s_j = w_1 \cdot P_{\text{collision},j} + w_2 \cdot m_j + w_3 \cdot \frac{1}{h_j} $$
ここで $P_{\text{collision},j}$ はデブリ $j$ の衝突確率、$m_j$ は質量（大きいほど衝突時の破片生成が多い）、$h_j$ は軌道高度（低いほど大気抵抗で自然減衰しにくい高度帯では優先度が上がる）、$w_1, w_2, w_3$ は重み係数です。
サービサー $i$ がデブリ $j$ を除去するためのコスト（$\Delta v$ — 必要な速度変化量）を $c_{ij}$ とすると、優先度つきの割当問題は次のように定式化できます。
$$ \max_{\bm{X}} \sum_{i=1}^{M}\sum_{j=1}^{K} s_j \cdot x_{ij} – \lambda \sum_{i=1}^{M}\sum_{j=1}^{K} c_{ij} \cdot x_{ij} $$
第1項は除去するデブリの優先度の合計を最大化し、第2項（$\lambda$ は正の重み）は $\Delta v$ コストを最小化します。これは標準的な割当問題の一般化であり、ハンガリアン法を拡張した手法やオークションアルゴリズムで解くことができます。
協調的な軌道遷移
フリート内のサービサーは、デブリに到達するために軌道遷移を行います。複数のサービサーが同時に軌道を変更する場合、互いの軌道が交差しないように調整する必要があります。これはフォーメーション制御と衝突回避の組み合わせ問題であり、先ほど解説した人工ポテンシャル場法が有効です。
さらに、サービサー間で燃料（$\Delta v$ 予算）の融通も重要です。あるサービサーの $\Delta v$ 残量が少なくなった場合、近傍のサービサーがタスクを引き継ぐ「再割当」が必要になります。これはコンセンサスアルゴリズムの応用として、各サービサーの $\Delta v$ 残量を分散的に共有し、負荷バランスを取ることで実現できます。
Pythonによるデブリ除去フリートの割当シミュレーション
優先度付きのデブリ割当を簡易的にシミュレーションしてみましょう。
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import linear_sum_assignment
np.random.seed(2026)
# サービサーとデブリの設定
M = 3 # サービサー数
K = 8 # デブリ数
# 軌道高度（km）
servicer_alt = np.random.uniform(700, 900, M)
debris_alt = np.random.uniform(600, 1000, K)
# 軌道傾斜角（度）
servicer_inc = np.random.uniform(96, 100, M)
debris_inc = np.random.uniform(95, 101, K)
# デブリの優先度（衝突確率と質量に基づく）
collision_prob = np.random.exponential(0.001, K)
mass = np.random.uniform(100, 5000, K)
priority = 0.5 * collision_prob / collision_prob.max() + 0.5 * mass / mass.max()
# Δv コスト（高度差と傾斜角差の関数として簡易計算）
delta_v_cost = np.zeros((M, K))
for i in range(M):
for j in range(K):
dh = abs(servicer_alt[i] - debris_alt[j])
di = abs(servicer_inc[i] - debris_inc[j])
delta_v_cost[i, j] = 0.05 * dh + 50 * np.radians(di) # km/s 近似
# 利益行列（優先度 - コストの重み付き）
lam = 0.3
profit_matrix = priority[np.newaxis, :] - lam * delta_v_cost
# 各サービサーは1つのデブリのみ担当（M < K なので、K個からM個を選ぶ）
# ダミー行を追加してK x Kにする
padded_profit = np.zeros((K, K))
padded_profit[:M, :] = profit_matrix
row_ind, col_ind = linear_sum_assignment(-padded_profit) # 最大化なので符号反転
# 実際の割当（ダミーを除外）
assignments = [(row_ind[i], col_ind[i])
for i in range(K) if row_ind[i] < M]
print("=== デブリ除去フリート割当結果 ===")
total_priority = 0
total_dv = 0
for srv, deb in assignments:
print(f" Servicer {srv+1} → Debris {deb+1} "
f"(priority={priority[deb]:.3f}, Δv={delta_v_cost[srv, deb]:.2f} km/s)")
total_priority += priority[deb]
total_dv += delta_v_cost[srv, deb]
print(f" 除去優先度合計: {total_priority:.3f}")
print(f" Δv合計: {total_dv:.2f} km/s")
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import linear_sum_assignment
# （前のコードブロックの変数を再利用する前提で再度設定）
np.random.seed(2026)
M, K = 3, 8
servicer_alt = np.random.uniform(700, 900, M)
debris_alt = np.random.uniform(600, 1000, K)
servicer_inc = np.random.uniform(96, 100, M)
debris_inc = np.random.uniform(95, 101, K)
collision_prob = np.random.exponential(0.001, K)
mass = np.random.uniform(100, 5000, K)
priority = 0.5 * collision_prob / collision_prob.max() + 0.5 * mass / mass.max()
delta_v_cost = np.zeros((M, K))
for i in range(M):
for j in range(K):
dh = abs(servicer_alt[i] - debris_alt[j])
di = abs(servicer_inc[i] - debris_inc[j])
delta_v_cost[i, j] = 0.05 * dh + 50 * np.radians(di)
lam = 0.3
profit_matrix = priority[np.newaxis, :] - lam * delta_v_cost
padded_profit = np.zeros((K, K))
padded_profit[:M, :] = profit_matrix
row_ind, col_ind = linear_sum_assignment(-padded_profit)
assignments = [(row_ind[i], col_ind[i]) for i in range(K) if row_ind[i] < M]
# 可視化（高度-傾斜角平面上でのサービサー-デブリ割当）
fig, ax = plt.subplots(1, 1, figsize=(10, 7))
# デブリ（サイズを優先度に比例）
scatter = ax.scatter(debris_inc, debris_alt, s=priority * 500 + 50,
c=priority, cmap='YlOrRd', alpha=0.8,
edgecolors='white', zorder=5)
for j in range(K):
ax.annotate(f"D{j+1}", (debris_inc[j], debris_alt[j]),
fontsize=8, ha='center', va='bottom', color='white')
# サービサー
ax.scatter(servicer_inc, servicer_alt, s=200, c='cyan',
marker='s', edgecolors='white', zorder=6, label='Servicers')
for i in range(M):
ax.annotate(f"S{i+1}", (servicer_inc[i], servicer_alt[i]),
fontsize=9, ha='center', va='bottom', color='cyan')
# 割当線
for srv, deb in assignments:
ax.plot([servicer_inc[srv], debris_inc[deb]],
[servicer_alt[srv], debris_alt[deb]],
'g-', linewidth=2, alpha=0.7)
plt.colorbar(scatter, label='Debris Priority Score')
ax.set_xlabel("Orbital Inclination [deg]")
ax.set_ylabel("Orbital Altitude [km]")
ax.set_title("Debris Removal Fleet Assignment")
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
この可視化では、高度-傾斜角平面上にデブリ（円、サイズと色が優先度に対応）とサービサー（水色の四角）をプロットし、最適割当を緑の線で結んでいます。3台のサービサーが8個のデブリの中から優先度と $\Delta v$ コストのバランスを考慮して最適な3個を選択していることがわかります。
結果を分析すると、必ずしも最も優先度が高い3個のデブリが選ばれているわけではない点が重要です。優先度が高くても $\Delta v$ コストが大きすぎるデブリは回避され、適度な優先度でアクセスしやすいデブリが選ばれます。これはフリート全体のリソース効率を反映しており、実際のミッション計画でもこのようなトレードオフ分析が不可欠です。
実用的な設計上の考慮事項
異機種フリートの協調
実際の軌道上ミッションでは、すべてのロボットが同一仕様であることは稀です。例えば、大型モジュールを搬送する力のあるロボットと、精密な接合作業を行う器用なロボットが混在する異機種フリート（heterogeneous fleet） が現実的です。
異機種フリートのタスク割当では、コスト行列がロボットの能力に依存します。ロボット $i$ がタスク $j$ を実行できない場合は $c_{ij} = \infty$ とすることで、ハンガリアン法の枠組み内で制約を表現できます。
耐故障性とグレースフルデグラデーション
宇宙環境では放射線によるシングルイベントアップセット（SEU）やマイクロメテオライトの衝突によるセンサ故障が避けられません。マルチロボットシステムの設計では、1台の故障がフリート全体を停止させない耐故障性（fault tolerance） が必須です。
分散制御の枠組みでは、故障したロボットを通信グラフから除外するだけで対応できます。ラプラシアン行列から対応する行と列を削除し、残りのロボットでコンセンサスを再実行します。タスク割当も、故障ロボットのタスクを残りのロボットに再割当することで対応できます。
重要なのは、$N$ 台中 $f$ 台が故障しても、$N – f$ 台で元のタスクの一定割合以上を遂行できる「グレースフルデグラデーション」特性を設計段階から保証することです。これは、タスクグラフの冗長性（同じタスクを実行できるロボットが複数存在すること）と通信グラフの $k$-連結性（$k-1$ 台の故障まで連結性を維持できること）の両方で実現されます。
エネルギーと燃料の制約
軌道上ロボットは太陽電池と蓄電池で動作し、軌道遷移にはスラスタの推進剤を消費します。エネルギーと推進剤は有限資源であり、協調計画にこれらの制約を組み込む必要があります。
各ロボットの $\Delta v$ 残量を $\Delta v_i^{\text{rem}}$ として、タスク割当の制約に追加します。
$$ \sum_{j \in \text{tasks}} c_{ij} \cdot x_{ij} \leq \Delta v_i^{\text{rem}} \quad \forall i $$
これはナップサック問題の一種であり、動的計画法や分枝限定法で解くことができます。また、コンセンサスアルゴリズムを用いて、フリート全体の $\Delta v$ 残量を分散的に均等化し、特定のサービサーだけが早期に燃料切れを起こさないように負荷分散を行うこともできます。
発展的トピック — 大規模スウォームへの道
バイオインスパイアードアルゴリズム
生物の群れ行動に着想を得たスウォームインテリジェンス（Swarm Intelligence） は、大規模なマルチロボットシステムに有効なアプローチです。
レイノルズの3ルール: 1987年にクレイグ・レイノルズが提案した3つの単純なルール — 分離（近すぎたら離れる）、整列（近傍の平均速度に合わせる）、結合（近傍の平均位置に向かう）— だけで、鳥の群れに似た複雑な集団行動が創発します。これをロボットに適用すると、個々のロボットの制御則は極めてシンプルでありながら、フリート全体として衝突回避・方向合わせ・群れ維持が自然に実現されます。
蟻コロニー最適化: アリがフェロモンを使って最短経路を発見する行動に着想を得た手法で、巡回セールスマン問題（デブリ除去の順序最適化）に適用できます。
粒子群最適化（PSO）: 鳥の群れの餌探し行動をモデル化した手法で、各粒子（ロボット）が自分の過去最良解と群れの過去最良解を参考にして位置を更新します。フォーメーションの最適配置計算などに利用されています。
自己再構成ロボット
究極のマルチロボットシステムとして、自己再構成ロボット（Self-Reconfigurable Robot） が研究されています。これは、多数の同一形状のモジュールが互いに結合・分離を繰り返して、全体の形状を変化させるシステムです。例えば、搬送時には長い列状に変形して狭い通路を通り、組立時には腕状に変形して部材を操作するといった柔軟な運用が可能になります。
軌道上での自己再構成ロボットは、NASAやMITなどが研究を進めており、将来的には宇宙ステーションの自律的な拡張・修復に活用されると期待されています。
Pythonでレイノルズルールによるスウォームシミュレーション
最後に、レイノルズの3ルールに基づくスウォームシミュレーションを実装し、群ロボットの創発的な集団行動を確認しましょう。
import numpy as np
import matplotlib.pyplot as plt
# パラメータ
N = 30 # ロボット数
dt = 0.05
T = 20.0
steps = int(T / dt)
domain = 50.0 # 空間サイズ
# レイノルズのルール重み
w_sep = 2.0 # 分離
w_ali = 1.0 # 整列
w_coh = 0.8 # 結合
# 各ルールの作用範囲
r_sep = 3.0 # 分離距離
r_neighbor = 10.0 # 近傍半径
v_max = 3.0 # 最大速度
# 初期状態
np.random.seed(42)
positions = np.random.uniform(0, domain, (N, 2))
velocities = np.random.uniform(-1, 1, (N, 2))
# 軌跡記録
traj = np.zeros((steps + 1, N, 2))
traj[0] = positions.copy()
for k in range(steps):
new_velocities = velocities.copy()
for i in range(N):
# 近傍の取得
diffs = positions - positions[i]
dists = np.linalg.norm(diffs, axis=1)
dists[i] = np.inf # 自分自身を除外
neighbors = np.where(dists < r_neighbor)[0]
if len(neighbors) == 0:
continue
# 1. 分離（Separation）: 近すぎる仲間から離れる
sep_force = np.zeros(2)
close = neighbors[dists[neighbors] < r_sep]
if len(close) > 0:
for j in close:
sep_force -= (positions[j] - positions[i]) / (dists[j] + 1e-6)
# 2. 整列（Alignment）: 近傍の平均速度に合わせる
ali_force = velocities[neighbors].mean(axis=0) - velocities[i]
# 3. 結合（Cohesion）: 近傍の重心に向かう
coh_force = positions[neighbors].mean(axis=0) - positions[i]
# 合成
new_velocities[i] += (w_sep * sep_force + w_ali * ali_force + w_coh * coh_force) * dt
# 速度制限
speed = np.linalg.norm(new_velocities[i])
if speed > v_max:
new_velocities[i] = new_velocities[i] / speed * v_max
velocities = new_velocities
positions = positions + velocities * dt
# 周期境界条件
positions = positions % domain
traj[k + 1] = positions.copy()
# 可視化（3つのスナップショット）
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
snapshots = [0, steps // 2, steps]
titles = ["t = 0 s", f"t = {T/2:.0f} s", f"t = {T:.0f} s"]
for ax, t_idx, title in zip(axes, snapshots, titles):
pos = traj[t_idx]
# 速度方向を矢印で表示（最終ステップに近い速度を推定）
if t_idx < steps:
vel = (traj[min(t_idx + 1, steps)] - traj[t_idx]) / dt
else:
vel = (traj[t_idx] - traj[t_idx - 1]) / dt
ax.quiver(pos[:, 0], pos[:, 1], vel[:, 0], vel[:, 1],
color='cyan', alpha=0.7, scale=50)
ax.scatter(pos[:, 0], pos[:, 1], c='orange', s=30, zorder=5, edgecolors='white')
ax.set_xlim(0, domain)
ax.set_ylim(0, domain)
ax.set_aspect('equal')
ax.set_title(title)
ax.grid(True, alpha=0.2)
fig.suptitle("Reynolds Flocking Rules — Swarm Robot Simulation", fontsize=13)
plt.tight_layout()
plt.show()
3つのスナップショットから、スウォームの時間発展が読み取れます。初期状態（$t=0$）ではロボットがランダムに散らばり、速度方向もバラバラです。しかし中間時点（$t=10$s）では、近傍同士が速度方向を揃え始め、小さなクラスターが形成されています。最終状態（$t=20$s）では、ロボットの大部分が同じ方向に移動する1〜数個の群れにまとまっています。
注目すべきは、各ロボットは「群れを作れ」という明示的な指令を受けていないにもかかわらず、3つの単純なローカルルール（分離・整列・結合）だけで群れ行動が創発している点です。これがスウォームインテリジェンスの本質であり、大規模マルチロボットシステムの設計原理として極めて強力です。1台1台に複雑な制御則をプログラムするのではなく、シンプルなローカルルールの設計だけで、スケーラブルな集団行動が実現できるのです。
まとめ
本記事では、マルチロボット協調の理論と軌道上応用について体系的に解説しました。
- 集中制御 vs 分散制御: 集中制御は大域的な最適性を保証できるが、単一障害点やスケーラビリティの問題がある。分散制御はロバストかつスケーラブルだが、大域最適性は保証されない
- グラフラプラシアン: 通信トポロジーをグラフで表現し、ラプラシアン行列の固有値（特に代数的連結度 $\lambda_2$）が合意の速度とネットワーク性能を決定する
- コンセンサスアルゴリズム: 局所通信だけで全ロボットが平均値に合意できることを理論的に証明し、スイッチングトポロジー下でも合併グラフの連結性条件のもとで合意が達成される
- タスク割当: ハンガリアン法（集中的、$O(N^3)$）とオークションアルゴリズム（分散的）でロボットとタスクの最適マッチングを計算できる
- フォーメーション制御: リーダーフォロワー法、仮想構造体法、コンセンサスベース法の3手法を解説し、用途に応じた使い分けを示した
- 軌道上応用: 大型構造物の組立とデブリ除去フリートに対して、タスク割当・フォーメーション制御・衝突回避を統合的に適用する方法を示した
マルチロボット協調は、宇宙開発の次の段階 — 大型構造物の自律組立、デブリ環境の持続的な清掃、メガコンステレーションの保守 — を支える基盤技術です。個々のロボットの性能向上だけでなく、「チームとしての賢さ」を設計する能力が、今後ますます重要になっていきます。
次のステップとして、以下の記事も参考にしてください。