この投稿では、四つのプロペラを持つクアッドコプター（ドローン）の 2 D（平面）シミュレーションを構築します。
運動方程式を導出し、それらを状態空間形式に書き換えて Python で数値解析する方法を示します。

---

なぜクアッドコプターをシミュレートするのか？

制御器を設計したり、強化学習エージェントを訓練したりする場合、まずソフトウェア上で動作できるモデルが必要です。
この投稿では、2 D で最初にシミュレーションステップをゼロから構築し、制御へ進む前にダイナミクスを理解します。

---

問題設定と座標系

下図の自由体図（free‑body diagram）から始めます。
水平軸を (y)、垂直軸を (z) とし、右手則で (x) 軸はページ外に向かうようにします。

クアッドコプター本体は質量 (m) を持つ剛体で、中心点 (C) から長さ (\ell) の二本のアームを伸ばしています。
(yz)-平面内で角度 (\phi) に回転し、(\phi>0) は (y) 軸から (z) 軸への逆時計回りの回転です。
重力は下向きに作用し、その大きさは (mg) です。
各アームには推力を生むプロペラがあり、推力は本体に対して垂直であるとします。これらを (F_1) と (F_2) と表記します。

---

運動方程式の導出

自由体図からニュートン–オイラー剛体運動方程式を得ます：
(y) 軸と (z) 軸に沿った二つの直線運動方程式と、角度 (\phi) に関する回転方程式です。
(F_y) と (F_z) はそれぞれの軸方向の合力成分、(\tau) は合せん断モーメント、(I) は重心を通る慣性モーメントです。時間微分は点で示します。

[ \begin{aligned} m,\ddot y &= \sum F_y\ m,\ddot z &= \sum F_z\ I,\ddot\phi &= \sum\tau \end{aligned} ]

推力を (y) と (z) の方向に分解すると

[ \begin{aligned} m,\ddot y &= -(F_1+F_2)\sin\phi\ m,\ddot z &= (F_1+F_2)\cos\phi - mg\ I,\ddot\phi &= (F_1-F_2),\ell \end{aligned} ]

---

状態空間形式への変換

数値解析を行うために、まず状態空間形式に書き直します。
状態 (\mathbf x) をシステムを記述する最小集合の変数と定義し、入力 (\mathbf u) を制御対象の変数とします。

制御可能な量はプロペラ推力 (F_1) と (F_2) です。
直接用いる代わりに次のように入力ベクトルを定義すると、表現がシンプルになります。

[ \mathbf u = \begin{bmatrix} u_1\u_2 \end{bmatrix}

\begin{bmatrix} F_1+F_2\[2pt] F_1-F_2 \end{bmatrix}, ]

これを代入し、微分項を整理すると

[ \begin{aligned} \ddot y &= -,\frac{u_1}{m}\sin\phi\ \ddot z &= ;\frac{u_1}{m}\cos\phi - g\ \ddot\phi &= ;\frac{u_2}{I},\ell \end{aligned} ]

状態空間モデルは一階微分方程式であるため、状態ベクトルを次のように定義します。

[ \mathbf x = \begin{bmatrix} y\ z\ \phi\ \dot y\ \dot z\ \dot\phi \end{bmatrix}. ]

すると完全な一階システムは

[ \boxed{ \dot{\mathbf x} = \begin{bmatrix} \dot y\[2pt] \dot z\[2pt] \dot\phi\[4pt] -\dfrac{u_1}{m}\sin\phi\[6pt] \dfrac{u_1}{m}\cos\phi - g\[6pt] \dfrac{u_2}{I},\ell \end{bmatrix} }. ]

---

Python でのシミュレーション

import numpy as np
from numpy.typing import NDArray

# 物理パラメータ
m = 0.8          # [kg]
g = 9.81         # [m/s²]
l = 0.5          # [m]
I = 1e-3         # [kg·m²]

def dynamics(x: NDArray, u: NDArray) -> NDArray:
    y, z, phi, y_dot, z_dot, phi_dot = x
    return np.array([
        y_dot,
        z_dot,
        phi_dot,
        -u[0]/m * np.sin(phi),
        u[0]/m * np.cos(phi) - g,
        u[1]/I * l
    ])

dynamics

が定義できたら、オイラー法で一次微分方程式を解きます。

t_start, t_stop, n_steps = 0.0, 15.0, 1_000
t   = np.linspace(t_start, t_stop, n_steps)
dt  = t[1] - t[0]

x      = np.zeros((n_steps, 6))          # 状態軌跡
u      = np.zeros((n_steps, 2))          # 入力軌跡

# 重さよりわずかに大きい定数推力 → ゆっくり上昇
u[:, 0] = m * g + 0.01        # u1
u[:, 1] = 0.0                # u2 (トルクなし)

for i in range(n_steps-1):
    x[i+1] = x[i] + dynamics(x[i], u[i]) * dt

これが、運動方程式を導出し、状態空間形式に変換して Python で数値積分する一連のシミュレーションパイプラインです。

---

シミュレーションの実行

トルクゼロの場合

(u_2=0) のとき、クアッドコプターは回転しません。
垂直加速度は一定であるため、(z(t)) は二次関数的に増加し、(\dot z(t)) は線形に増加します。
(y) と (\phi) は定数のままです。

---

トルクありの場合

今度は非ゼロトルク入力 (u_2) を与え、合計推力 (u_1) も増やします。
この単純モデルには地面接触を考慮していないため、クアッドコプターが (z=0) 以下に到達したらシミュレーションを停止します。

# 新しい入力
u[:, 0] = m * g + 0.5     # より大きな推力
u[:, 1] = 5e-5            # 小さなトルク

eps = 1e-8
for i in range(n_steps-1):
    if x[i, 1] < -eps:   # クアッドコプターが墜落した
        t = t[:i]
        x = x[:i]
        u = u[:i]
        break
    x[i+1] = x[i] + dynamics(x[i], u[i]) * dt

角度はゆっくりと 0 から約 3 rad（ほぼ逆転）へ増加します。
最初は高さが上がりますが、(\phi) が大きくなるにつれて推力の垂直成分が減少し、((u_1/m)\cos\phi < g) になると重力が支配的になり、機体は落下します。

---

付録：プロットコード

import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator

fig, axes = plt.subplots(3, 1, sharex=True, figsize=(9, 7))

# 位置状態
ax0, ax0r = axes[0], axes[0].twinx()
ax0.plot(t, x[:, 0], label='y', color='C0')
ax0.plot(t, x[:, 1], label='z', color='C1')
ax0r.plot(t, x[:, 2], label=r'$\phi$', color='C2')
ax0.set_ylabel('Position [m]')
ax0r.set_ylabel('Angle [rad]')
ax0.set_title('位置状態')
lines0 = ax0.get_lines() + ax0r.get_lines()
ax0.legend(handles=lines0, labels=[l.get_label() for l in lines0])

# 速度状態
ax1, ax1r = axes[1], axes[1].twinx()
ax1.plot(t, x[:, 3], label=r'$\dot y$', color='C0')
ax1.plot(t, x[:, 4], label=r'$\dot z$', color='C1')
ax1r.plot(t, x[:, 5], label=r'$\dot\phi$', color='C2')
ax1.set_ylabel('Linear velocity [m/s]')
ax1r.set_ylabel('Angular velocity [rad/s]')
ax1.set_title('速度状態')
lines1 = ax1.get_lines() + ax1r.get_lines()
ax1.legend(handles=lines1, labels=[l.get_label() for l in lines1])

# 入力
axes[2].plot(t, u[:, 0], label='u₁')
axes[2].plot(t, u[:, 1], label='u₂')
axes[2].set_title('入力')
axes[2].set_ylabel('Input [N]')
axes[2].legend()

for ax in axes:
    ax.xaxis.set_major_locator(MultipleLocator(1))
    ax.set_xlabel('time [s]')
    ax.grid(True)

fig.suptitle('2‑D クアッドコプターシミュレーション')
fig.tight_layout()
plt.show()

注: 本投稿は完全に人間が執筆し、GPT‑5.4（Codex を経由）で校正／軽微編集を行いました。