力学系の基礎と古典力学と制御工学

1.1 ニュートン力学と微分方程式

運動方程式 $F = ma$

ニュートンの第2法則は、質点の加速度 $a = \frac{d^2x}{dt^2}$ がその質量 $m$ に比例し、外力 $F$ によって決まることを表す。これを微分方程式として記述すれば：

$m \frac{d^2x}{dt^2} = F(x, \dot{x}, t)$

これは**常微分方程式（ODE）**の形をとり、解析的解や数値解法を用いることで未来の運動を予測可能とする。

質点系と剛体運動

• 質点系：複数の質点を含む系は、重心運動と相対運動に分解して記述できる。
• 剛体運動：並進運動と回転運動を統合的に記述し、回転行列や慣性テンソルによってそのダイナミクスを表す。

一般化座標と拘束条件

自由度を削減するため、位置を記述する座標を ${q_i}$ で表現。拘束条件は：

$f(q_1, q_2, ..., q_n, t) = 0$

のように関数関係で与えられ、力学系を抽象的に扱う際の基盤となる。

---

1.2 ラグランジュの運動方程式

仮想仕事とダランベールの原理

拘束付き運動では、ニュートン力学が直接使えない場合がある。ダランベールの原理は、仮想変位に対する力の合計仕事がゼロであることを前提に：

$\sum_i (F_i - m_i \ddot{q}_i) \delta q_i = 0$

を導き、運動方程式の導出に役立つ。

ラグランジアン $L = T - V$

• $T$：運動エネルギー
• $V$：位置エネルギー

これらの差によって構成される：

$L(q, \dot{q}, t) = T(q, \dot{q}) - V(q)$

オイラー・ラグランジュ方程式

作用積分 $\int L dt$ を極値にする条件として：

$\frac{d}{dt} \left( \frac{\partial L}{\partial \dot{q}_i} \right) - \frac{\partial L}{\partial q_i} = 0$

この形式は、座標変換に対して不変であり、汎用的な力学系に適用できる。

---

1.3 ハミルトン力学と正準形式

共役変数とハミルトニアン

運動量：

$p_i = \frac{\partial L}{\partial \dot{q}_i}$

ハミルトニアンは：

$H(q, p) = \sum_i p_i \dot{q}_i - L(q, \dot{q})$

正準方程式

ハミルトン力学系の基本方程式：

$\frac{dq_i}{dt} = \frac{\partial H}{\partial p_i}, \quad \frac{dp_i}{dt} = -\frac{\partial H}{\partial q_i}$

ポアソン括弧と保存則

任意の物理量 $f(q, p)$ の時間発展は：

$\frac{df}{dt} = {f, H} + \frac{\partial f}{\partial t}$

ポアソン括弧 ${f, g}$ がゼロであれば、$f$ と $g$ は可換であり、保存量の候補となる。

Noetherの定理と対称性

ラグランジアンがある連続対称変換に不変である場合、その対称性に対応する保存則（運動量保存、エネルギー保存など）が存在する。

# Program Name: mechanics_lagrange_hamilton.py
# Creation Date: 20250605
# Overview: Simulate classical mechanics using Newtonian, Lagrangian, and Hamiltonian formalisms
# Usage: Run in Python environment with required libraries installed

# -------------------- Install Required Libraries / ライブラリのインストール --------------------
!pip install numpy scipy matplotlib --quiet

# -------------------- Import Libraries / ライブラリのインポート --------------------
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# -------------------- Parameters / パラメータの一元管理 --------------------
m = 1.0          # 質量 [kg] / Mass [kg]
k = 1.0          # ばね定数 [N/m] / Spring constant [N/m]
x0 = 1.0         # 初期位置 [m] / Initial position [m]
v0 = 0.0         # 初期速度 [m/s] / Initial velocity [m/s]
t_span = (0, 20) # 時間範囲 [s] / Time span [s]
t_eval = np.linspace(t_span[0], t_span[1], 1000)  # 解の評価点 / Evaluation points

# -------------------- Newtonian Mechanics / ニュートン力学による運動方程式 --------------------
# dx/dt = v, dv/dt = -k/m * x

def newton_eq(t, y):
    # y[0] = x, y[1] = v
    dxdt = y[1]
    dvdt = -k / m * y[0]
    return [dxdt, dvdt]

sol_newton = solve_ivp(newton_eq, t_span, [x0, v0], t_eval=t_eval)

# -------------------- Lagrangian Mechanics / ラグランジュ形式 --------------------
# L = T - V = 1/2 m v^2 - 1/2 k x^2
# Euler-Lagrange 方程式により Newton と同一式が導出されるため再利用可能

def lagrangian_eq(t, y):
    # Lagrangian: L = 1/2 m v^2 - 1/2 k x^2
    return newton_eq(t, y)

sol_lagrange = solve_ivp(lagrangian_eq, t_span, [x0, v0], t_eval=t_eval)

# -------------------- Hamiltonian Mechanics / ハミルトン形式 --------------------
# p = ∂L/∂v = mv, H = p^2/2m + 1/2 k x^2
# dq/dt = ∂H/∂p = p/m, dp/dt = -∂H/∂q = -kx

def hamilton_eq(t, y):
    # y[0] = q (位置), y[1] = p (運動量)
    dqdt = y[1] / m
    dpdt = -k * y[0]
    return [dqdt, dpdt]

p0 = m * v0  # 初期運動量 / Initial momentum
sol_hamilton = solve_ivp(hamilton_eq, t_span, [x0, p0], t_eval=t_eval)

# -------------------- Plot Results / 結果の可視化 --------------------
plt.plot(sol_newton.t, sol_newton.y[0], label="Newtonian")
plt.plot(sol_lagrange.t, sol_lagrange.y[0], '--', label="Lagrangian")
plt.plot(sol_hamilton.t, sol_hamilton.y[0], ':', label="Hamiltonian")
plt.title("Comparison of Newtonian, Lagrangian, and Hamiltonian Mechanics")
plt.xlabel("Time [s]")
plt.ylabel("Position [m]")
plt.legend()
plt.grid(True)
plt.show()

第2章：力学系と制御理論の接続

2.1 状態空間と制御可能性

状態方程式化

ハミルトン系の変数 $(q, p)$ をまとめて $x = [q; p]$ とすれば、状態方程式：

$\dot{x} = f(x, u)$

可制御性と可観測性

線形系：$\dot{x} = Ax + Bu$ において、可制御性行列：

$\mathcal{C} = [B, AB, A^2B, ..., A^{n-1}B]$

の階数が $n$ ならば可制御。出力 $y = Cx$ に基づく可観測性も同様に定義される。

線形近似

非線形力学系を平衡点まわりで線形化することで、制御設計を行いやすくする。

---

2.2 ラグランジュ系の制御

外力の導入と拡張

制御入力 $u_i$ を右辺に導入：

$\frac{d}{dt} \left( \frac{\partial L}{\partial \dot{q}_i} \right) - \frac{\partial L}{\partial q_i} = Q_i^{\text{ext}}(u)$

トルク制御と運動計画

関節駆動型のメカニズム（ロボットアーム等）において、軌道計画 $q_d(t)$ に基づき制御を設計。

多入力系への拡張

複数の一般化座標を持つ系に対して、それぞれの運動方程式を構成し、制御器を多入力系として設計する。

---

2.3 ハミルトン形式と最適制御

ハミルトニアンとPontryaginの最適性原理

コスト関数：

$J = \int_0^T L(x,u) dt$

最小化のためのハミルトニアン：

$\mathcal{H}(x, u, \lambda) = L(x, u) + \lambda^\top f(x, u)$

最適性条件

• 状態：$\dot{x} = \frac{\partial \mathcal{H}}{\partial \lambda}$
• 共役状態：$\dot{\lambda} = -\frac{\partial \mathcal{H}}{\partial x}$
• 最適性：$\frac{\partial \mathcal{H}}{\partial u} = 0$

制御則の比較

• ラグランジュ型：作用汎関数の極値問題として記述
• ハミルトン型：状態・共役状態の動的最適化問題として記述

軌道最適化などの応用においては、ハミルトン形式の方が自然に拡張される。

# Program Name: control_theory_mechanics.py
# Creation Date: 20250605
# Overview: Connect classical mechanics with modern control theory via state-space modeling
# Usage: Run in Python environment with required libraries installed

# -------------------- Install Required Libraries / ライブラリのインストール --------------------
!pip install numpy scipy matplotlib control --quiet

# -------------------- Import Libraries / ライブラリのインポート --------------------
import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg as la
from control import ss, ctrb, obsv, lqr, forced_response

# -------------------- Parameters / パラメータの一元管理 --------------------
m = 1.0      # 質量 / Mass [kg]
k = 1.0      # ばね定数 / Spring constant [N/m]
x0 = 1.0     # 初期位置 / Initial position [m]
v0 = 0.0     # 初期速度 / Initial velocity [m/s]
u_func = lambda t: 0.0  # 入力関数 / Input function

# -------------------- State-Space Representation / 状態空間表現 --------------------
# 二次元状態 x = [位置; 速度] / Define state vector x = [position; velocity]
A = np.array([[0, 1],
              [-k/m, 0]])
B = np.array([[0],
              [1/m]])
C = np.array([[1, 0]])  # 出力は位置のみ / Output is position
D = np.array([[0]])

sys = ss(A, B, C, D)

# -------------------- Controllability and Observability / 可制御性と可観測性 --------------------
Ctrb = ctrb(A, B)  # 可制御性行列 / Controllability matrix
Obsv = obsv(A, C)  # 可観測性行列 / Observability matrix
rank_ctrb = np.linalg.matrix_rank(Ctrb)  # 可制御性階数 / Rank of controllability matrix
rank_obsv = np.linalg.matrix_rank(Obsv)  # 可観測性階数 / Rank of observability matrix

# -------------------- LQR Optimal Control / LQR 最適制御 --------------------
Q = np.eye(2)       # 状態コスト / State cost matrix
R = np.array([[1]]) # 入力コスト / Input cost matrix
K, _, _ = lqr(sys, Q, R)  # 最適フィードバックゲイン / Optimal feedback gain

# 状態フィードバック制御系（閉ループ）/ Closed-loop system
A_cl = A - B @ K
sys_cl = ss(A_cl, B, C, D)

# -------------------- Simulate Step Response / ステップ応答による確認 --------------------
T = np.linspace(0, 20, 500)
U = np.zeros_like(T)  # 入力ゼロ / Zero input
X0 = [x0, v0]         # 初期状態 / Initial state

T, yout = forced_response(sys_cl, T, U, X0)  # forced_responseの出力を2変数で受け取る

# -------------------- Plot Results / 結果の可視化 --------------------
plt.plot(T, yout)
plt.title("Closed-loop Response with LQR Control")
plt.xlabel("Time [s]")
plt.ylabel("Position [m]")
plt.grid(True)
plt.show()

# -------------------- Display Rank Results / 可制御・可観測の結果表示 --------------------
print(f"Controllability Rank: {rank_ctrb} / 可制御性階数")
print(f"Observability Rank: {rank_obsv} / 可観測性階数")
print(f"LQR Gain K: {K} / LQR制御ゲイン")