AIエージェントは台車型倒立振子を制御できるか？

はじめに

どんな人向けの記事？

• LLM や AI エージェントの性能が気になる方
• 生成 AI を用いた web アプリ開発に興味がある方
• 制御工学に興味のある方

デプロイ環境

ミニPC（M30A）： N100(4core)@3.40GHz, 16GB, Ubuntu Server 24.04 LTS

概要

• 移動台車型倒立振子の制御シミュレーションを行う Web アプリを開発
• Cursor + o3-mini + Clade 3.5 Sonnet で開発
• 制御ゲイン設計を Cursor + Clade 3.5 Sonnet にしてもらって、システムを安定化できるかを検証

作成した Web アプリ

まず、今回作成した Web アプリはこちら。

https://inverted-pendulum.robustonian.com

4 つの制御ゲインを決めて、移動台車型倒立振子の安定化を目指すゲーム。
制御工学に詳しい人間はもちろん、LLM や AI エージェントの性能が高ければ線形領域でのゲイン設計が実現できることがわかった。

https://x.com/gosrum/status/1893830860727623696

これが何を意味しているかというと、LLM や AI エージェントの性能評価にも利用できるということ。

線形領域と非線形領域の制御で難易度は段違いだが、とりあえず本アプリケーションでは線形フィードバック入力のみを対象としているので、線形領域での安定化が実現できるか否かを確認することにする。今後、非線形制御ができるように拡張したい。

記事の構成

今回の記事では、下記の順番で解説していく。それぞれ分野がかなり異なる内容なので、興味のある項目から読んでいただければと思う。

• Web アプリの構成と開発の流れ

• 移動台車型倒立振子のダイナミクスと制御系設計

• LLM と AI エージェントの性能評価用プロンプト

Web アプリの構成と開発の流れ

ゲームの基本ルール

設定項目

アプリの構成と開発の流れ

アプリの構成

アプリの構成を図示したものを以下に示す。

技術スタック

開発の流れ

最初、Cursor composer（agent モード）で Claude 3.5 Sonnet のみを使って開発をしたところ、ダイナミクスが全然整合しないという問題が発生した。

そこで、下記のような流れで開発をしたらうまくいった。

Claude 3.5 Sonnet も Web アプリ開発では無双している印象があったが、数式の変形等を扱う際は、o3-mini のような推論モデルを使う方が良さそう、という知見が実体験として得られた。

ちなみに o3-mini は Composer(agent)モードではなく、Chat モードで使用した。

おそらく上記は Cursor を使わなくても、ChatGPT アプリの o3-mini と Cline+Claude 3.5 Sonnet でも全く同じようなことができると思う。

台車型倒立振子のダイナミクス

以下は、反時計回りを正とする（すなわち、$\theta$ は鉛直上から反時計回りの回転を正とする）場合の台車‐倒立振子系の定式化例である。ここでは、ラグランジアンの手法を用いて非線形の運動方程式を導出する。

---

1. 座標の定義

• カートの水平位置： $x$
• 振子の長さ： $\ell$
• 振子の質量： $m$
• カートの質量： $M$
• 重力加速度： $g$
• 角度： $\theta$
  • $\theta = 0$ のとき、振子は「鉛直上向き」（不安定平衡）となる。
  • 正の$\theta$は反時計回り（つまり、垂直方向から左に回転）とする。

このとき、振子の重心の座標は

$$\begin{aligned} x_p &= x - \ell\sin\theta, \\ y_p &= \ell\cos\theta, \end{aligned}$$

と定義される。

$\theta$が正のとき、$\sin\theta > 0$となるため、$x_p$はカートの位置$x$より左側にずれる。

2. 運動エネルギーと位置エネルギー

(1) 運動エネルギー $T$

カートの運動エネルギー

$$T_{\text{cart}} = \frac{1}{2}M\dot{x}^2.$$

振子の運動エネルギー

振子の重心の座標は

$$x_p = x - \ell\sin\theta,\quad y_p = \ell\cos\theta.$$

これの時間微分は

$$\dot{x}_p = \dot{x} - \ell\,\dot{\theta}\cos\theta,\quad \dot{y}_p = -\ell\,\dot{\theta}\sin\theta.$$

よって、振子の運動エネルギーは

$$\begin{aligned} T_{\text{pendulum}} &= \frac{1}{2}m\Bigl[(\dot{x}_p)^2+(\dot{y}_p)^2\Bigr] \\ &=\frac{1}{2}m\Bigl[(\dot{x}-\ell\cos\theta\,\dot{\theta})^2+(\ell\sin\theta\,\dot{\theta})^2\Bigr] \\ &=\frac{1}{2}m\left[\dot{x}^2-2\ell\cos\theta\,\dot{x}\dot{\theta}+\ell^2\dot{\theta}^2\ (\cos^2\theta+\sin^2\theta)\right] \\ &=\frac{1}{2}m\left[\dot{x}^2-2\ell\cos\theta\,\dot{x}\dot{\theta}+\ell^2\dot{\theta}^2\right]. \end{aligned}$$

全体の運動エネルギー

$$T = T_{\text{cart}} + T_{\text{pendulum}} = \frac{1}{2}(M+m)\dot{x}^2 - m\ell\cos\theta\,\dot{x}\dot{\theta} + \frac{1}{2}m\ell^2\dot{\theta}^2.$$

(2) 位置エネルギー $V$

カートは水平方向にしか移動しないため位置エネルギーを考慮しなくて良い。そのため、振子の重心についてのみ考える。重心の鉛直座標は $ y_p = \ell\cos\theta $ であるから、下記のように書ける:

$$V = mgy_p = mg\ell\cos\theta.$$

※ 本来、位置エネルギーは「低い位置を基準」にすることが多いが、ここでは定数項は系の動力学に影響しないため、便宜的にこの形で扱う。

---

3. ラグランジアンの定式化

ラグランジアンは

$$L = T - V,$$

より、

$$L = \frac{1}{2}(M+m)\dot{x}^2 - m\ell\cos\theta\,\dot{x}\dot{\theta} + \frac{1}{2}m\ell^2\dot{\theta}^2 - mg\ell\cos\theta.$$

---

4. Euler–Lagrange の方程式から運動方程式を導出

外力としてカートに $F$（水平方向の力）が作用すると仮定する。一般化座標 $q$ に対して、

$$\frac{d}{dt}\frac{\partial L}{\partial \dot{q}}-\frac{\partial L}{\partial q}=Q,$$

と書ける。ここで、$Q$は $q$ に対応する外力（$x$の場合は $F$、$\theta$の場合は外部トルクは作用しないものとする）である。

(1) $x$ に関する微分方程式

まず、

$$\frac{\partial L}{\partial \dot{x}} = (M+m)\dot{x} - m\ell\cos\theta\,\dot{\theta}.$$

これの時間微分は

$$\frac{d}{dt}\frac{\partial L}{\partial \dot{x}} = (M+m)\ddot{x} + m\ell\sin\theta\,\dot{\theta}^2 - m\ell\cos\theta\,\ddot{\theta}.$$

一方、$x$には明示的な依存はないため、

$$\frac{\partial L}{\partial x}=0.$$

よって、外力 $F$ を考慮すると

$$(M+m)\ddot{x} + m\ell\sin\theta\,\dot{\theta}^2 - m\ell\cos\theta\,\ddot{\theta}=F.$$

(2) $\theta$ に関する方程式

まず、

$$\frac{\partial L}{\partial \dot{\theta}} = -m\ell\cos\theta\,\dot{x} + m\ell^2\dot{\theta}.$$

その時間微分は

$$\frac{d}{dt}\frac{\partial L}{\partial \dot{\theta}} = m\ell\sin\theta\,\dot{x}\,\dot{\theta} - m\ell\cos\theta\,\ddot{x} + m\ell^2\ddot{\theta}.$$

次に、$\theta$ に関する偏微分は、

• 項 $-m\ell\cos\theta\,\dot{x}\dot{\theta}$ について：$\displaystyle\frac{\partial}{\partial \theta}(-m\ell\cos\theta\,\dot{x}\dot{\theta}) = m\ell\sin\theta\,\dot{x}\dot{\theta}$
• 位置エネルギー項 $-mg\ell\cos\theta$ について：$\displaystyle\frac{\partial}{\partial \theta}(-mg\ell\cos\theta) = mg\ell\sin\theta$

より、

$$\frac{\partial L}{\partial \theta} = m\ell\sin\theta\,\dot{x}\dot{\theta} + mg\ell\sin\theta.$$

したがって、Euler–Lagrange の式より

$$\begin{aligned} \frac{d}{dt}\frac{\partial L}{\partial \dot{\theta}} - \frac{\partial L}{\partial \theta} &= \Bigl[m\ell\sin\theta\,\dot{x}\,\dot{\theta} - m\ell\cos\theta\,\ddot{x} + m\ell^2\ddot{\theta}\Bigr] - \Bigl[m\ell\sin\theta\,\dot{x}\dot{\theta} + mg\ell\sin\theta\Bigr] \\ &= - m\ell\cos\theta\,\ddot{x} + m\ell^2\ddot{\theta} - mg\ell\sin\theta \\ &= 0. \end{aligned}$$

両辺を $m\ell$ で割ると、

$$-\cos\theta\,\ddot{x} + \ell\ddot{\theta} - g\sin\theta = 0.$$

整理すると、

$$\ell\ddot{\theta} - \cos\theta\,\ddot{x} - g\sin\theta = 0.$$

---

5. 最終的な運動方程式

以上の結果より、反時計回りを正とする場合の台車‐倒立振子の運動方程式は

$$\boxed{ \begin{aligned} (M+m)\ddot{x} + m\ell\sin\theta\,\dot{\theta}^2 - m\ell\cos\theta\,\ddot{\theta} &= F,\\ \ell\ddot{\theta} - \cos\theta\,\ddot{x} - g\sin\theta &= 0. \end{aligned} }$$

---

ダイナミクスの定式化まとめ

1. 座標設定：
   振子の重心は

$$x_p = x - \ell\sin\theta,\quad y_p = \ell\cos\theta$$

と定義し、$\theta=0$を不安定平衡（鉛直上向き）とする。

2. エネルギーの計算：
   カートと振子それぞれの運動エネルギーおよび振子の位置エネルギーを求め、ラグランジアン $ L = T - V $ を構成した。

3. Euler–Lagrange の公式の適用：
   各一般座標 $ x $ と $\theta$ に対して運動方程式を導出し、外力 $ F $ を $ x $ の方程式に導入した。

この結果、反時計回りを正として定式化した場合の台車‐倒立振子系の運動方程式が上記のように得られる。

更に整理する

以下のように、まず第２式から $\ddot{\theta}$ を解き、それを第１式に代入することで、$\ddot{x}$ と $\ddot{\theta}$ を明示的に求めることができる。

---

ステップ 1: 第２式から $\ddot{\theta}$ を解く

与えられた第２式は

$$\ell\,\ddot{\theta}-\cos\theta\,\ddot{x}-g\sin\theta=0.$$

これを $\ddot{\theta}$ について解くと、

$$\ddot{\theta}=\frac{\cos\theta\,\ddot{x}+g\sin\theta}{\ell}\,.$$

---

ステップ 2: 上記の式を第１式に代入

第１式は

$$(M+m)\ddot{x}+m\ell\sin\theta\,\dot{\theta}^2-m\ell\cos\theta\,\ddot{\theta}=F\,.$$

ここに、先ほどの $\ddot{\theta}$ の式を代入すると、

$$(M+m)\ddot{x}+m\ell\sin\theta\,\dot{\theta}^2 - m\ell\cos\theta\left(\frac{\cos\theta\,\ddot{x}+g\sin\theta}{\ell}\right)=F\,.$$

まずは $\ell$ の項が打ち消せるので、

$$(M+m)\ddot{x}+m\ell\sin\theta\,\dot{\theta}^2 - m\cos\theta\Bigl(\cos\theta\,\ddot{x}+g\sin\theta\Bigr)=F\,.$$

展開すると

$$(M+m)\ddot{x} - m\cos^2\theta\,\ddot{x} + m\ell\sin\theta\,\dot{\theta}^2 - m g\sin\theta\cos\theta = F\,.$$

$\ddot{x}$ の項をまとめると

$$\Bigl[(M+m)-m\cos^2\theta\Bigr]\ddot{x} + m\ell\sin\theta\,\dot{\theta}^2 - m g\sin\theta\cos\theta = F\,.$$

したがって、

$$\ddot{x} = \frac{F + m g\sin\theta\cos\theta - m\ell\sin\theta\,\dot{\theta}^2}{(M+m)-m\cos^2\theta}\,.$$

さらに、$1-\cos^2\theta=\sin^2\theta$ であることから

$$(M+m)-m\cos^2\theta=M+m(1-\cos^2\theta)=M+m\sin^2\theta\,.$$

よって、最終的に

$$\boxed{ \ddot{x}=\frac{F + m g\sin\theta\cos\theta - m\ell\sin\theta\,\dot{\theta}^2}{M+m\sin^2\theta}\,. }$$

---

ステップ 3: $\ddot{\theta}$ の式も整理する

先ほど求めた

$$\ddot{\theta}=\frac{\cos\theta\,\ddot{x}+g\sin\theta}{\ell}$$

に、上記の $\ddot{x}$ の式を代入すると、

$$\ddot{\theta}=\frac{1}{\ell}\left(\cos\theta\,\frac{F + m g\sin\theta\cos\theta - m\ell\sin\theta\,\dot{\theta}^2}{M+m\sin^2\theta}+g\sin\theta\right)\,.$$

分母をそろえるために $g\sin\theta$ にも $M+m\sin^2\theta$ をかけて書くと、

$$\ddot{\theta}=\frac{F\cos\theta + m g\,\sin\theta\cos^2\theta - m\ell\sin\theta\cos\theta\,\dot{\theta}^2 + g\sin\theta(M+m\sin^2\theta)}{\ell\,(M+m\sin^2\theta)}\,.$$

ここで $m g\,\sin\theta\cos^2\theta + g\sin\theta(M+m\sin^2\theta)$ について考えると、

$$m g\,\sin\theta\cos^2\theta+g\sin\theta(M+m\sin^2\theta)=g\sin\theta\Bigl[M+m(\cos^2\theta+\sin^2\theta)\Bigr] =g\sin\theta(M+m)\,.$$

よって、整理後の $\ddot{\theta}$ は

$$\boxed{ \ddot{\theta}=\frac{F\cos\theta - m\ell\sin\theta\cos\theta\,\dot{\theta}^2 + g\sin\theta(M+m)}{\ell\,(M+m\sin^2\theta)}\,. }$$

---

最終的な整理結果

整理すると、連立微分方程式は

$$\boxed{ \begin{aligned} \ddot{x} &= \frac{F + m g\sin\theta\cos\theta - m\ell\sin\theta\,\dot{\theta}^2}{M+m\sin^2\theta},\\ \ddot{\theta} &= \frac{F\cos\theta - m\ell\sin\theta\cos\theta\,\dot{\theta}^2 + g\sin\theta(M+m)}{\ell\,(M+m\sin^2\theta)}\,. \end{aligned} }$$

これで、もともとの連立微分方程式が $\ddot{x}$ と $\ddot{\theta}$ について整理された。

線形化、状態空間表現、制御系設計

以下は、非線形ダイナミクスを微小角近似（$\sin\theta\approx\theta$, $\cos\theta\approx1$）で線形化した結果、状態空間表現として記述した場合の例である。状態変数を

$$X = \begin{pmatrix} x \\ v \\ \theta \\ \omega \end{pmatrix} \quad \text{（ただし } v = \dot{x},\ \omega = \dot{\theta} \text{）}$$

と定義する。

制御系設計

1. 線形化の手順

もともとの系は

$$\begin{aligned} &(M+m)\ddot{x} + m\ell\sin\theta\,\dot{\theta}^2 - m\ell\cos\theta\,\ddot{\theta} = F,\\ &\ell\ddot{\theta} - \cos\theta\,\ddot{x} - g\sin\theta = 0, \end{aligned}$$

であるが、微小角近似により、次のように近似できる。

1. 小角近似の適用：
   $\sin\theta\approx\theta$ および $\cos\theta\approx 1$ とすると、

• 第 1 式：

$$(M+m)\ddot{x} - m\ell\,\ddot{\theta} = F.$$

• 第 2 式：

$$\ell\,\ddot{\theta} - \ddot{x} - g\,\theta = 0,$$

すなわち

$$\ddot{\theta} = \frac{\ddot{x}+g\,\theta}{\ell}\,.$$

2. $\ddot{x}$ の式の導出：
   第 2 式の $\ddot{\theta}$ を第 1 式に代入すると

$$(M+m)\ddot{x} - m\ell \left(\frac{\ddot{x}+g\,\theta}{\ell}\right) = F,$$

$$(M+m)\ddot{x} - m\ddot{x} - m\,g\,\theta = F,$$

$$M\,\ddot{x} - m\,g\,\theta = F.$$

したがって、

$$\ddot{x} = \frac{F}{M} + \frac{m\,g}{M}\theta.$$

3. $\ddot{\theta}$ の式の導出：
   先の結果を第 2 式より

$$\ddot{\theta} = \frac{1}{\ell}\Bigl(\ddot{x}+g\,\theta\Bigr) = \frac{1}{\ell}\left(\frac{F}{M} + \frac{m\,g}{M}\theta + g\,\theta\right) = \frac{F}{\ell\,M} + \frac{(M+m)g}{\ell\,M}\theta.$$

---

2. 状態空間表現

上記の結果を元に、状態変数の各成分の時間発展は下記のように書ける：

$$\begin{cases} \dot{x} = v,\\ \dot{v} = \dfrac{F}{M} + \dfrac{m\,g}{M}\,\theta,\\ \dot{\theta} = \omega,\\ \dot{\omega} = \dfrac{F}{\ell\,M} + \dfrac{(M+m)g}{\ell\,M}\,\theta. \end{cases}$$

これは、

$$\dot{X} = A\,X + B\,F,$$

という形に書き直すことができ、係数行列 $A$ と入力行列 $B$ は以下のようになる。

---

3. 行列 $A$ と $B$ の具体的な表現

• 行列 $A$:

$$A = \begin{pmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & \dfrac{m\,g}{M} & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & \dfrac{(M+m)g}{\ell\,M} & 0 \end{pmatrix}\,.$$

• 入力行列 $B$:

$$B = \begin{pmatrix} 0 \\ \dfrac{1}{M} \\ 0 \\ \dfrac{1}{\ell\,M} \end{pmatrix}\,.$$

これらの行列を用いることで、状態フィードバック制御

$$F = -K\,X = -\bigl[k_x x + k_v v + k_\theta \theta + k_\omega \omega\bigr]$$

によって閉ループ系

$$\dot{X} = (A - BK)X$$

を設計する際の基礎となる。

---

制御系設計まとめ

線形化した状態空間表現は

$$\dot{X} = \underbrace{\begin{pmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & \dfrac{m\,g}{M} & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & \dfrac{(M+m)g}{\ell\,M} & 0 \end{pmatrix}}_{A}\,X + \underbrace{\begin{pmatrix} 0 \\ \dfrac{1}{M} \\ 0 \\ \dfrac{1}{\ell\,M} \end{pmatrix}}_{B} F\,.$$

この $A$ と $B$ をもとに、たとえばポール配置法や LQR を用いて状態フィードバックゲイン $K = [k_x \; k_v \; k_\theta \; k_\omega]$ を設計することで、システムを安定化させることが可能である。

今回は微小角近似を仮定しているため、上記の制御系設計では初期角度が小さい領域からしか安定化できない点に注意する。

AI エージェントの性能評価用プロンプト

上記の線形化や極配置・最適制御ロジックを使えば線形領域の安定化を実現するフィードバックゲイン設計は可能である。

ここでは、非線形ダイナミクスとゴールを与えたときに最適なゲインを求めるまでプロセスを AI エージェントに実行させてみて、実際に安定化できるかの実験を行った。

プロンプト

下記は台車型倒立振子のダイナミクスを表しています。

## ダイナミクス

$$
\boxed{
\begin{aligned}
\ddot{x} &= \frac{F + m g\sin\theta\cos\theta - m\ell\sin\theta\,\dot{\theta}^2}{M+m\sin^2\theta}\,,\\[1mm]
\ddot{\theta} &= \frac{F\cos\theta - m\ell\sin\theta\cos\theta\,\dot{\theta}^2 + g\sin\theta(M+m)}{\ell\,(M+m\sin^2\theta)}\,.
\end{aligned}
}
$$

このダイナミクスに対して、下記のフィードバック入力で安定化を目指します。
ここでは、原点（x=v=\theta=\omega=0）での安定化をゴールとします。

## フィードバック入力

$$
F =-k_xx-k_vv-k_\theta-k_\omega\omega
$$

ここで、$v = \dot{x}, \omega=\dot{\theta}$です。

## ゴール

4 つの状態をすべて 0 に収束させるためのゲイン$k_x,k_v,k_\theta,k_\omega$を求めてください。

## 備考

-基本的なパラメータは以下の通りとします。

$$
M=1,m=0.5,l=2, g=9.81
$$

- 初期状態は$x=0,v=0,\omega=0$
- 初期角度が大きいほど得点が大きくなります（180 度で満点）。
- 重要な制約条件として、x の可動範囲は[-4.75,4.75]です。これを超えたらゲームオーバーです。

Cursor+Claude 3.5 Sonnet の出力

出力

ご依頼の内容を理解しました。台車型倒立振子の制御問題について、適切なフィードバックゲインを設計していきましょう。

<タスク分析>
主要タスク：

1. 非線形システムの線形化
2. 線形化システムに対する状態フィードバックゲインの設計
3. シミュレーションによる検証
4. 制約条件の確認

重要な要件と制約：

• システムパラメータ: M=1, m=0.5, l=2, g=9.81
• 初期条件: x=0, v=0, θ=5deg, ω=0
• 状態の制約: |x| ≤ 4.75
• 目標: すべての状態を 0 に収束

潜在的な課題：

1. 非線形システムの線形化の精度
2. 制約条件を満たすゲイン設計
3. 過渡応答の制御

実行ステップ：

1. システムの線形化
2. 状態空間表現の導出
3. 制御系設計
4. シミュレーションコードの実装
5. 結果の検証と調整
   </タスク分析>

まず、システムの実装を行っていきましょう。Python を使用してシミュレーションを行います。

cart_pole_control.py

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# System parameters
M = 1.0    # Cart mass
m = 0.5    # Pendulum mass
l = 2.0    # Pendulum length
g = 9.81   # Gravity acceleration

# Control gains (to be tuned)
K = np.array([-100, -50, 800, 100])  # Initial guess for [kx, kv, ktheta, komega]

def dynamics(state, t):
    x, v, theta, omega = state

    # Control input (state feedback)
    F = -K[0]*x - K[1]*v - K[2]*theta - K[3]*omega

    # System dynamics
    sin_theta = np.sin(theta)
    cos_theta = np.cos(theta)

    x_ddot = (F + m*g*sin_theta*cos_theta - m*l*sin_theta*omega**2) / (M + m*sin_theta**2)
    theta_ddot = (F*cos_theta - m*l*sin_theta*cos_theta*omega**2 + g*sin_theta*(M+m)) / (l*(M + m*sin_theta**2))

    return [v, x_ddot, omega, theta_ddot]

# Initial conditions
theta_0 = np.deg2rad(5)  # 5 degrees in radians
x0 = [0, 0, theta_0, 0]  # [x, v, theta, omega]

# Time vector
t = np.linspace(0, 10, 1000)

# Solve ODE
solution = odeint(dynamics, x0, t)

# Plot results
plt.figure(figsize=(12, 8))

# Cart position
plt.subplot(2, 2, 1)
plt.plot(t, solution[:, 0])
plt.grid(True)
plt.xlabel('Time [s]')
plt.ylabel('Cart Position [m]')
plt.title('Cart Position vs Time')

# Cart velocity
plt.subplot(2, 2, 2)
plt.plot(t, solution[:, 1])
plt.grid(True)
plt.xlabel('Time [s]')
plt.ylabel('Cart Velocity [m/s]')
plt.title('Cart Velocity vs Time')

# Pendulum angle
plt.subplot(2, 2, 3)
plt.plot(t, np.rad2deg(solution[:, 2]))
plt.grid(True)
plt.xlabel('Time [s]')
plt.ylabel('Angle [deg]')
plt.title('Pendulum Angle vs Time')

# Angular velocity
plt.subplot(2, 2, 4)
plt.plot(t, np.rad2deg(solution[:, 3]))
plt.grid(True)
plt.xlabel('Time [s]')
plt.ylabel('Angular Velocity [deg/s]')
plt.title('Angular Velocity vs Time')

plt.tight_layout()
plt.show()

# Check if position constraint is violated
max_pos = np.max(np.abs(solution[:, 0]))
print(f"Maximum absolute cart position: {max_pos:.3f} m")
if max_pos > 4.75:
    print("Warning: Position constraint violated!")

次に、システムの線形化を行い、適切なゲインを設計します。線形化システムの解析用のコードを作成します。

linear_analysis.py

import numpy as np
from scipy import linalg

# System parameters
M = 1.0    # Cart mass
m = 0.5    # Pendulum mass
l = 2.0    # Pendulum length
g = 9.81   # Gravity acceleration

# Linearized state-space matrices at the upright position
A = np.array([
    [0, 1, 0, 0],
    [0, 0, m*g/M, 0],
    [0, 0, 0, 1],
    [0, 0, g*(M+m)/(M*l), 0]
])

B = np.array([
    [0],
    [1/M],
    [0],
    [1/(M*l)]
])

# Desired pole locations
# We want relatively fast response but avoid oscillations
poles = np.array([-2, -2.1, -2.2, -2.3])

# Calculate feedback gains using pole placement
K = linalg.place_poles(A, B, poles).gain_matrix

print("Designed feedback gains:")
print(f"k_x = {K[0, 0]:.2f}")
print(f"k_v = {K[0, 1]:.2f}")
print(f"k_theta = {K[0, 2]:.2f}")
print(f"k_omega = {K[0, 3]:.2f}")

# Verify closed-loop eigenvalues
A_cl = A - B @ K
eig_vals = np.linalg.eigvals(A_cl)
print("\nClosed-loop eigenvalues:")
for i, eig in enumerate(eig_vals):
    print(f"λ_{i+1} = {eig:.3f}")

では、このシステムを実行して結果を確認しましょう。

python linear_analysis.py

申し訳ありません。必要なパッケージをインストールする必要があります。

uv init && uv add numpy scipy matplotlib

では、再度線形解析を実行してみましょう。

uv run linear_analysis.py

申し訳ありません。place_polesの代わりにplaceを使用するように修正します。

linear_analysis.py

import numpy as np
from scipy import linalg

# System parameters
M = 1.0    # Cart mass
m = 0.5    # Pendulum mass
l = 2.0    # Pendulum length
g = 9.81   # Gravity acceleration

# Linearized state-space matrices at the upright position
A = np.array([
    [0, 1, 0, 0],
    [0, 0, m*g/M, 0],
    [0, 0, 0, 1],
    [0, 0, g*(M+m)/(M*l), 0]
])

B = np.array([
    [0],
    [1/M],
    [0],
    [1/(M*l)]
])

# Desired pole locations
# We want relatively fast response but avoid oscillations
poles = np.array([-2, -2.1, -2.2, -2.3])

# Calculate feedback gains using pole placement
K = linalg.place(A, B, poles)

print("Designed feedback gains:")
print(f"k_x = {K[0, 0]:.2f}")
print(f"k_v = {K[0, 1]:.2f}")
print(f"k_theta = {K[0, 2]:.2f}")
print(f"k_omega = {K[0, 3]:.2f}")

# Verify closed-loop eigenvalues
A_cl = A - B @ K
eig_vals = np.linalg.eigvals(A_cl)
print("\nClosed-loop eigenvalues:")
for i, eig in enumerate(eig_vals):
    print(f"λ_{i+1} = {eig:.3f}")

では、再度実行してみましょう。

uv run linear_analysis.py

申し訳ありません。代わりに LQR（Linear Quadratic Regulator）を使用して制御ゲインを設計しましょう。

linear_analysis.py

import numpy as np
from scipy import linalg

# System parameters
M = 1.0    # Cart mass
m = 0.5    # Pendulum mass
l = 2.0    # Pendulum length
g = 9.81   # Gravity acceleration

# Linearized state-space matrices at the upright position
A = np.array([
    [0, 1, 0, 0],
    [0, 0, m*g/M, 0],
    [0, 0, 0, 1],
    [0, 0, g*(M+m)/(M*l), 0]
])

B = np.array([
    [0],
    [1/M],
    [0],
    [1/(M*l)]
])

# LQR weight matrices
Q = np.diag([1, 1, 10, 1])  # State cost matrix
R = np.array([[0.1]])       # Control cost

# Solve the Riccati equation
P = linalg.solve_continuous_are(A, B, Q, R)

# Calculate the optimal feedback gain
K = np.linalg.inv(R) @ B.T @ P

print("Designed feedback gains:")
print(f"k_x = {K[0, 0]:.2f}")
print(f"k_v = {K[0, 1]:.2f}")
print(f"k_theta = {K[0, 2]:.2f}")
print(f"k_omega = {K[0, 3]:.2f}")

# Verify closed-loop eigenvalues
A_cl = A - B @ K
eig_vals = np.linalg.eigvals(A_cl)
print("\nClosed-loop eigenvalues:")
for i, eig in enumerate(eig_vals):
    print(f"λ_{i+1} = {eig:.3f}")

では、再度実行してみましょう。

uv run linear_analysis.py

LQR によって得られた制御ゲインを使用して、非線形シミュレーションを実行してみましょう。

cart_pole_control.py

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# System parameters
M = 1.0    # Cart mass
m = 0.5    # Pendulum mass
l = 2.0    # Pendulum length
g = 9.81   # Gravity acceleration

# Control gains from LQR design
K = np.array([-3.16, -7.10, 77.34, 32.47])

def dynamics(state, t):
    x, v, theta, omega = state

    # Control input (state feedback)
    F = -K[0]*x - K[1]*v - K[2]*theta - K[3]*omega

    # System dynamics
    sin_theta = np.sin(theta)
    cos_theta = np.cos(theta)

    x_ddot = (F + m*g*sin_theta*cos_theta - m*l*sin_theta*omega**2) / (M + m*sin_theta**2)
    theta_ddot = (F*cos_theta - m*l*sin_theta*cos_theta*omega**2 + g*sin_theta*(M+m)) / (l*(M + m*sin_theta**2))

    return [v, x_ddot, omega, theta_ddot]

# Initial conditions
theta_0 = np.deg2rad(5)  # 5 degrees in radians
x0 = [0, 0, theta_0, 0]  # [x, v, theta, omega]

# Time vector
t = np.linspace(0, 10, 1000)

# Solve ODE
solution = odeint(dynamics, x0, t)

# Plot results
plt.figure(figsize=(12, 8))

# Cart position
plt.subplot(2, 2, 1)
plt.plot(t, solution[:, 0])
plt.grid(True)
plt.xlabel('Time [s]')
plt.ylabel('Cart Position [m]')
plt.title('Cart Position vs Time')
plt.axhline(y=4.75, color='r', linestyle='--', label='Position limit')
plt.axhline(y=-4.75, color='r', linestyle='--')
plt.legend()

# Cart velocity
plt.subplot(2, 2, 2)
plt.plot(t, solution[:, 1])
plt.grid(True)
plt.xlabel('Time [s]')
plt.ylabel('Cart Velocity [m/s]')
plt.title('Cart Velocity vs Time')

# Pendulum angle
plt.subplot(2, 2, 3)
plt.plot(t, np.rad2deg(solution[:, 2]))
plt.grid(True)
plt.xlabel('Time [s]')
plt.ylabel('Angle [deg]')
plt.title('Pendulum Angle vs Time')

# Angular velocity
plt.subplot(2, 2, 4)
plt.plot(t, np.rad2deg(solution[:, 3]))
plt.grid(True)
plt.xlabel('Time [s]')
plt.ylabel('Angular Velocity [deg/s]')
plt.title('Angular Velocity vs Time')

plt.tight_layout()
plt.savefig('cart_pole_response.png')
plt.close()

# Check if position constraint is violated
max_pos = np.max(np.abs(solution[:, 0]))
print(f"Maximum absolute cart position: {max_pos:.3f} m")
if max_pos > 4.75:
    print("Warning: Position constraint violated!")

では、シミュレーションを実行してみましょう。

uv run cart_pole_control.py

シミュレーション結果を分析しましょう。LQR で設計した制御ゲインを使用した結果：

1. 制御ゲイン：

   • k_x = -3.16 (位置フィードバック)
   • k_v = -7.10 (速度フィードバック)
   • k_theta = 77.34 (角度フィードバック)
   • k_omega = 32.47 (角速度フィードバック)

2. システムの応答：

   • 最大カート位置: 0.398 m (制約の 4.75 m 以内に収まっています)
   • 閉ループ極：すべて左半平面にあり、システムは安定です
     • λ_1 = -4.201 (実極)
     • λ_2 = -2.701 (実極)
     • λ_3,4 = -1.118 ± 0.341j (複素共役極)

3. 制御性能：

   • カート位置は制約範囲内（±4.75 m）に十分収まっています
   • 初期角度 5 度からの整定が実現できています
   • 複素極の存在により、若干の振動的な応答が予想されますが、減衰も十分です

この LQR 制御器は、以下の特徴を持っています：

1. 状態重み行列 Q = diag([1, 1, 10, 1])を使用し、特に角度 θ の制御を重視
2. 制御入力重み R = 0.1 として、比較的大きな制御入力を許容
3. 結果として得られた制御ゲインは、特に角度フィードバックが大きく(77.34)、振子の直立制御を重視

シミュレーション結果から、この制御設計は要求仕様を満たしていると判断できます：

• 状態がすべて原点に収束
• カート位置の制約条件（±4.75 m）を満たす
• 初期角度 5 度からの安定化を実現

一度しか試していないので、運が良かっただけかもしれないけれど、一度目のトライで安定化に成功！ちなみに、python 仮想化環境のコマンドは uv を使うように手直ししている。単なる実行環境の違いなので、修正しなくても結果は変わらないはず。

https://x.com/gosrum/status/1893830860727623696

ちなみに、シミュレーション結果の時系列グラフも出力してくれるプログラムになっていた。

Cline とか、他の LLM との組み合わせでも試してみたい。それと、それらの結果を DB に保存してを Web アプリ上から確認できるようにしたい。そうすれば、制御工学という一分野になってしまうが、一つのベンチマークになるかもしれない。

おそらく線形領域の安定化なら難なくできてしまうケースが多いと思うので、ゆくゆくは非線形制御の UI も追加して、難易度をあげて LLM や AI エージェントの限界を評価したい。

おそらく 40 度ぐらいまでなら、今の入力の形

$$F = -K\,X = -\bigl[k_x x + k_v v + k_\theta \theta + k_\omega \omega\bigr]$$

でもわりといけるはず。90 度に近づくほど（範囲の制約を超過しないこと）難しくなって、90 度になると安定化そのものが難しくなる。

非線形制御に詳しい方や勘の良い方なら、初期角度 90° が特異な状態となっていることがわかるはず。つまり、90° 未満なら理論上は左に動けば安定化は実現できる（制約は超えてしまうかも）が、90° の姿勢では初手で左に動くだけでは安定化できないことがわかるだろう。

補足：Rules for AI

ちなみに、Rules for AI には下記を記載しています。

Rules for AI

あなたは高度な問題解決能力を持つAIアシスタントです。以下の指示に従って、効率的かつ正確にタスクを遂行してください。

まず、ユーザーから受け取った指示を確認します：

<指示>
{{instructions}}
</指示>

この指示を元に、以下のプロセスに従って作業を進めてください。なお、すべての提案と実装は、記載された技術スタックの制約内で行ってください：

1. 指示の分析と計画
<タスク分析>
- 主要なタスクを簡潔に要約してください。
- 記載された技術スタックを確認し、その制約内での実装方法を検討してください。
- 重要な要件と制約を特定してください。
- 潜在的な課題をリストアップしてください。
- タスク実行のための具体的なステップを詳細に列挙してください。
- それらのステップの最適な実行順序を決定してください。
- 必要となる可能性のあるツールやリソースを考慮してください。

このセクションは、後続のプロセス全体を導くものなので、時間をかけてでも、十分に詳細かつ包括的な分析を行ってください。
</タスク分析>

2. タスクの実行
- 特定したステップを一つずつ実行してください。
- 各ステップの完了後、簡潔に進捗を報告してください。
- 実行中に問題や疑問が生じた場合は、即座に報告し、対応策を提案してください。

3. 品質管理
- 各タスクの実行結果を迅速に検証してください。
- エラーや不整合を発見した場合は、直ちに修正アクションを実施してください。
- コマンドを実行する場合は、必ず標準出力を確認し、結果を報告してください。

4. 最終確認
- すべてのタスクが完了したら、成果物全体を評価してください。
- 当初の指示内容との整合性を確認し、必要に応じて調整を行ってください。

重要な注意事項：
- 不明点がある場合は、作業開始前に必ず確認を取ってください。
- 重要な判断が必要な場合は、その都度報告し、承認を得てください。
- 予期せぬ問題が発生した場合は、即座に報告し、対応策を提案してください。

このプロセスに従って、効率的かつ正確にタスクを遂行してください。

また、重要な事実として、OpenAIにはgpt-4oとgpt-4o-miniというモデルが存在します。

これは、ニケちゃんさんの記事に記載されているものをそのまま採用させていただいてます。もともと、きのぴーさんが公開されたものをニケちゃんさんが一部改良されたものとのことです。

https://note.com/nike_cha_n/n/nd0f7566019ae

正直今回のタスクと関係ない部分もありますが、安定化に成功したのでヨシとします。

まとめ

今回の記事はかなり盛り沢山な内容となっていしまいましたが、下記についてまとめました。

• 台車型倒立振子の Web アプリ開発について
• 台車型倒立振子のダイナミクスと線形領域の制御系設計について
• Cursor+Claude 3.5 Sonnet による制御系設計のプロンプトと出力について

今後はいろんなモデルで安定化できるか試したり、非線形制御への拡張もできたらと思ってます。

最後まで読んでいただきありがとうございました。次回もぜひよろしくお願いします。

おまけ：ちょっとした小話

ちなみに気づいている方もいると思いますが、本アプリはマニュアル制御モードも搭載しています。人間が操作して安定化できたらクリア！みたいなコンセプトで最初作り始めたのですが、難易度が高すぎて諦めました 😇

その後にフィードバック制御モードを追加して、制御分野に強い人しかクリアできないのでは？と思いつつ、今なら生成 AI 使って簡単にクリアできるのでは？という挑戦状としてアプリ公開、記事の執筆をしました。

生成 AI を使って何かにトライしてみたいと考えている方は、是非挑戦してみてください。