LQR(線形二次レギュレータ)の理論と実装

制御系の設計において、「どのような制御が最も良いか」という問いには複数のアプローチがあります。極配置法では「極をどこに置くか」を設計者が直接指定しました。しかし、多次元のシステムでは「良い極配置」を直感的に選ぶのが難しくなります。4つの状態変数を持つ倒立振子の制御で、8つの極をどこに配置するのが最適かを直感で決めるのは困難です。

もっと自然な問いかけは「状態の偏差を小さくしつつ、制御入力も抑えたい。この2つのバランスが最も良い制御は何か？」でしょう。これは日常的な場面でも馴染みのある考え方です。たとえば車を運転するとき、目的地にできるだけ早く着きたい（状態の偏差を小さく）けれど、急ハンドルや急ブレーキは避けたい（制御入力を抑えたい）。この2つの相反する要求の「ちょうどいいバランス」を数学的に厳密に求める方法があります。

LQR（Linear Quadratic Regulator、線形二次レギュレータ）は、この問いに対する厳密な最適解を与えます。「状態の二次形式 + 入力の二次形式」のコスト関数を最小化する状態フィードバックゲインを、リカッチ方程式を解くことで一意に求めます。

LQR制御は1960年代にR.E.カルマンらによって確立され、現代制御理論の中心的な結果の一つです。その美しさは、複雑な最適化問題が1つの代数方程式（リカッチ方程式）を解くだけに帰着されるという点にあります。しかも、得られた制御則は単純な状態フィードバック $\bm{u} = -\bm{K}\bm{x}$ の形をしており、実装も容易です。

LQR制御を理解すると、以下のような応用が開けます。

最適な状態フィードバック: 性能と制御エネルギーのトレードオフを数学的に最適化。設計者は「何を重視するか」を重み行列で指定するだけでよい
ロバスト性の保証: LQR制御はSISO系においてゲイン余裕 $\geq 6$ dB、位相余裕 $\geq 60^\circ$ を自動的に保証する
航空宇宙: 航空機の飛行制御（F-16の飛行制御系設計など）、宇宙機の姿勢制御
ロボティクス: 倒立振子の安定化、二足歩行ロボットの制御、ドローンの飛行制御
LQG制御への発展: カルマンフィルタと組み合わせてノイズ環境下の最適制御へ

本記事の内容

LQRのコスト関数と最適制御問題の定式化
ハミルトン・ヤコビ・ベルマン方程式を経由したリカッチ方程式の導出
重み行列 $\bm{Q}$, $\bm{R}$ の設計指針とBrysonの法則
ロバスト性の保証とその証明のスケッチ
具体的な数値計算例
Pythonでの倒立振子の実装と検証

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

状態フィードバック制御と極配置法 — 状態フィードバックの基礎
状態空間モデルの定義と構成 — 状態空間モデル

LQRのコスト関数

問題設定

LQR問題を定式化する前に、「最適制御」とは何かを直感的に理解しましょう。

ある山の頂上にボールを置くと、ボールは重力に従って転がり落ちます。もしボールが最も早く麓に到達する経路を選べるとしたら、それが「最適制御」です。ただし、経路を制御するにはエネルギー（制御入力）が必要で、使えるエネルギーには限りがあります。LQR問題は、この「到達速度（性能）」と「エネルギー消費（コスト）」のバランスを最適化する問題です。

連続時間線形システム $\dot{\bm{x}} = \bm{A}\bm{x} + \bm{B}\bm{u}$ に対し、次のコスト関数（評価関数）を最小化する状態フィードバック $\bm{u} = -\bm{K}\bm{x}$ を求めます。

$$ \begin{equation} J = \int_0^{\infty} (\bm{x}^T\bm{Q}\bm{x} + \bm{u}^T\bm{R}\bm{u}) \, dt \end{equation} $$

$\bm{Q} \in \mathbb{R}^{n \times n}$: 状態への重み行列（半正定値 $\bm{Q} \geq 0$）
$\bm{R} \in \mathbb{R}^{m \times m}$: 入力への重み行列（正定値 $\bm{R} > 0$）

コスト関数の各項の物理的意味を深堀りしましょう。

第1項 $\bm{x}^T\bm{Q}\bm{x}$ は「状態の偏差のペナルティ」です。対角的な $\bm{Q} = \text{diag}(q_1, q_2, \dots, q_n)$ の場合、$\bm{x}^T\bm{Q}\bm{x} = q_1 x_1^2 + q_2 x_2^2 + \cdots + q_n x_n^2$ となり、各状態変数の偏差の2乗に $q_i$ の重みをかけた和になります。$q_i$ が大きい状態変数ほど、偏差が厳しく罰せられます。

第2項 $\bm{u}^T\bm{R}\bm{u}$ は「制御入力のペナルティ」です。$\bm{R}$ が大きいほど大きな入力が罰せられるため、控えめな制御を好む設定になります。物理的には、アクチュエータの電力消費、モーターの発熱、燃料消費などに対応します。

なぜ二次形式（2乗の形）を使うのでしょうか。これには3つの理由があります。第一に、偏差の符号によらずペナルティが正になるため、「ずれの大きさ」を適切に評価できます。第二に、二次形式は微分が線形になるため、最適解が線形フィードバックの形になります。第三に、二次形式のコスト関数を持つ最適制御問題は解析的に解けるという数学的な利点があります。

$\bm{Q}$ を大きくすると偏差の抑制を重視し（応答が速いが入力が大きい）、$\bm{R}$ を大きくすると入力の節約を重視します（応答は遅いが入力が小さい）。重要なのは $\bm{Q}$ と $\bm{R}$ の絶対的な大きさではなく相対的な比であるという点です。$\bm{Q}$ と $\bm{R}$ を同じ定数倍しても、最適ゲイン $\bm{K}$ は変わりません。

重み行列の設計指針

重み行列の選び方は実用上の最大の課題です。以下に主要な設計法を紹介します。

Brysonの法則: 最も広く使われる経験的な方法です。各状態変数と制御入力の許容最大値 $x_{i,\max}$、$u_{j,\max}$ を用いて、

$$ Q_{ii} = \frac{1}{x_{i,\max}^2}, \quad R_{jj} = \frac{1}{u_{j,\max}^2} $$

と設定します。たとえば、倒立振子で角度の許容偏差が $0.1$ rad、カートの位置の許容偏差が $0.5$ m であれば、$Q_{33} = 100$、$Q_{11} = 4$ とします。

物理的考慮に基づく方法: コスト関数の各項に物理的な意味を持たせます。たとえば、$\bm{x}^T\bm{Q}\bm{x}$ をエネルギーに比例する量として設定し、$\bm{u}^T\bm{R}\bm{u}$ を制御コスト（電力消費 $\propto I^2$）として設定します。

出力重みの利用: $\bm{Q} = \bm{C}^T\bm{C}$ と設定すると、出力 $\bm{y} = \bm{C}\bm{x}$ の2乗和を最小化することになります。観測できる量の偏差を重視する自然な選択です。

問題設定が整理できたところで、次にこの最適化問題の解を導出しましょう。

リカッチ方程式の導出

最適制御の解を導くために、ハミルトン・ヤコビ・ベルマン（HJB）方程式のアプローチを使います。

まず、コスト関数の最小値を $V(\bm{x}) = \bm{x}^T\bm{P}\bm{x}$ という二次形式で表すことを仮定します。これは最適値関数と呼ばれ、現在の状態 $\bm{x}$ から出発して最適制御を適用したときの最小コストを表します。

HJB方程式は次の形をしています。

$$ 0 = \min_{\bm{u}} \left[\bm{x}^T\bm{Q}\bm{x} + \bm{u}^T\bm{R}\bm{u} + \left(\frac{\partial V}{\partial \bm{x}}\right)^T (\bm{A}\bm{x} + \bm{B}\bm{u})\right] $$

$V = \bm{x}^T\bm{P}\bm{x}$ を代入すると $\frac{\partial V}{\partial \bm{x}} = 2\bm{P}\bm{x}$ であるから、

$$ 0 = \min_{\bm{u}} \left[\bm{x}^T\bm{Q}\bm{x} + \bm{u}^T\bm{R}\bm{u} + 2\bm{x}^T\bm{P}(\bm{A}\bm{x} + \bm{B}\bm{u})\right] $$

角括弧内を $\bm{u}$ で偏微分してゼロとおくと、最適入力が求まります。

$$ 2\bm{R}\bm{u} + 2\bm{B}^T\bm{P}\bm{x} = \bm{0} $$

$$ \bm{u}^* = -\bm{R}^{-1}\bm{B}^T\bm{P}\bm{x} $$

したがって、最適ゲインは

$$ \begin{equation} \bm{K} = \bm{R}^{-1}\bm{B}^T\bm{P} \end{equation} $$

この最適入力をHJB方程式に代入して整理すると、$\bm{P}$ が満たすべき方程式が得られます。

$$ \bm{x}^T\bm{Q}\bm{x} + \bm{x}^T\bm{P}\bm{B}\bm{R}^{-1}\bm{R}\bm{R}^{-1}\bm{B}^T\bm{P}\bm{x} – 2\bm{x}^T\bm{P}\bm{B}\bm{R}^{-1}\bm{B}^T\bm{P}\bm{x} + 2\bm{x}^T\bm{P}\bm{A}\bm{x} = 0 $$

整理して、$\bm{R}^{-1}\bm{R}\bm{R}^{-1} = \bm{R}^{-1}$ に注意すると、

$$ \bm{x}^T(\bm{Q} – \bm{P}\bm{B}\bm{R}^{-1}\bm{B}^T\bm{P} + \bm{P}\bm{A} + \bm{A}^T\bm{P})\bm{x} = 0 $$

これが任意の $\bm{x}$ で成り立つためには、括弧内がゼロでなければなりません。こうして連続代数リカッチ方程式（CARE: Continuous Algebraic Riccati Equation）が得られます。

$$ \begin{equation} \bm{A}^T\bm{P} + \bm{P}\bm{A} – \bm{P}\bm{B}\bm{R}^{-1}\bm{B}^T\bm{P} + \bm{Q} = \bm{0} \end{equation} $$

この方程式の各項の意味を確認しましょう。

$\bm{A}^T\bm{P} + \bm{P}\bm{A}$: 自由応答（制御なし）によるコストの伝播
$-\bm{P}\bm{B}\bm{R}^{-1}\bm{B}^T\bm{P}$: 最適制御によるコストの削減（常に負）
$\bm{Q}$: 状態偏差によるコストの追加

解の存在条件として、$(\bm{A}, \bm{B})$ が可安定（stabilizable）かつ $(\bm{A}, \bm{Q}^{1/2})$ が可検出（detectable）であれば、CARE は一意の正半定値解を持ちます。より強い条件として $(\bm{A}, \bm{B})$ が可制御かつ $(\bm{A}, \bm{Q}^{1/2})$ が可観測であれば、解は正定値になります。

最適コストの計算

リカッチ方程式の解 $\bm{P}$ を使うと、最適コストは初期状態の二次形式で表されます。

$$ J^* = \bm{x}(0)^T\bm{P}\bm{x}(0) $$

つまり、$\bm{P}$ は「各状態方向の偏差のコスト感度」を表す行列です。$\bm{P}$ の固有値が大きい方向の状態偏差は、コストに大きく寄与します。

ロバスト性の保証

LQR制御は安定余裕の面で優れた性質を持ちます。これはLQR制御の「ボーナス」とも言える重要な結果です。

SISO系のLQR制御は以下を自動的に保証します。

ゲイン余裕: $-6$ dB 以上 $\sim$ $+\infty$（$1/2 \leq |K_{\text{gain}}| < \infty$）
位相余裕: $\pm 60^\circ$ 以上

この結果のスケッチを示しましょう。LQRのループ伝達関数は $L(s) = \bm{K}(s\bm{I} – \bm{A})^{-1}\bm{B}$ ですが、カルマン等式（Kalman equality）として知られる関係式、

$$ (1 + L(j\omega))^*(1 + L(j\omega)) = 1 + \bm{B}^T(-j\omega\bm{I} – \bm{A}^T)^{-1}\bm{Q}(j\omega\bm{I} – \bm{A})^{-1}\bm{B} / R $$

の右辺が常に1以上であることから、$|1 + L(j\omega)| \geq 1$ が成り立ちます。これはナイキスト線図上で $L(j\omega)$ が $(-1, 0)$ を中心とする単位円の外側にあることを意味し、上記の安定余裕が導かれます。

極配置法にはこのような保証がないため、ロバスト性の面でLQRは大きなアドバンテージを持っています。ただし、この保証はSISO系に限られ、MIMO系では一般に成り立たないことに注意が必要です。

理論の基盤が整ったところで、次にPythonで倒立振子の制御を実装して検証しましょう。

Pythonでの実装と検証

import numpy as np
import matplotlib.pyplot as plt
import control as ctrl
from scipy.linalg import solve_continuous_are

# 倒立振子のモデル
# 状態: [x, dx, theta, dtheta]
M = 1.0   # カート質量
m = 0.1   # 振子質量
l = 0.5   # 振子長さ
g = 9.81
b = 0.1   # 摩擦

# 線形化されたモデル
A = np.array([
    [0, 1, 0, 0],
    [0, -b/M, -m*g/M, 0],
    [0, 0, 0, 1],
    [0, b/(M*l), (M+m)*g/(M*l), 0]
])
B = np.array([[0], [1/M], [0], [-1/(M*l)]])
C = np.array([[1, 0, 0, 0], [0, 0, 1, 0]])
D = np.zeros((2, 1))

print(f"開ループ固有値: {np.linalg.eigvals(A)}")

# LQR設計（重みを変えて比較）
configs = [
    ("Light control", np.diag([1, 0, 10, 0]), np.array([[1.0]])),
    ("Balanced", np.diag([10, 1, 100, 1]), np.array([[0.1]])),
    ("Aggressive", np.diag([100, 10, 1000, 10]), np.array([[0.01]])),
]

fig, axes = plt.subplots(2, 2, figsize=(14, 9))

t = np.linspace(0, 5, 1000)
x0 = np.array([0.5, 0, 0.1, 0])  # 初期: 少しずれた位置と角度

for label, Q, R in configs:
    # リカッチ方程式を解く
    P = solve_continuous_are(A, B, Q, R)
    K = np.linalg.inv(R) @ B.T @ P

    A_cl = A - B @ K
    sys_cl = ctrl.StateSpace(A_cl, B, C, D)

    eigs = np.linalg.eigvals(A_cl)
    print(f"\n{label}: K = {K.flatten()}")
    print(f"  閉ループ固有値: {eigs}")

    # 初期値応答をシミュレーション
    from scipy.integrate import solve_ivp
    def cl_dynamics(t, x):
        u = (-K @ x.reshape(-1, 1)).flatten()
        return (A @ x + B.flatten() * u[0])

    sol = solve_ivp(cl_dynamics, [0, 5], x0, t_eval=t, max_step=0.01)

    # カート位置
    axes[0, 0].plot(sol.t, sol.y[0], linewidth=2, label=label)
    # 振子角度
    axes[0, 1].plot(sol.t, np.degrees(sol.y[2]), linewidth=2, label=label)
    # 制御入力
    u_hist = np.array([(-K @ sol.y[:, i]).flatten()[0] for i in range(len(sol.t))])
    axes[1, 0].plot(sol.t, u_hist, linewidth=2, label=label)

# カート位置
axes[0, 0].axhline(0, color="black", linestyle="--", alpha=0.4)
axes[0, 0].set_xlabel("Time [s]", fontsize=11)
axes[0, 0].set_ylabel("Cart position [m]", fontsize=11)
axes[0, 0].set_title("Cart Position", fontsize=12)
axes[0, 0].legend(fontsize=9)
axes[0, 0].grid(True, alpha=0.3)

# 振子角度
axes[0, 1].axhline(0, color="black", linestyle="--", alpha=0.4)
axes[0, 1].set_xlabel("Time [s]", fontsize=11)
axes[0, 1].set_ylabel("Pendulum angle [deg]", fontsize=11)
axes[0, 1].set_title("Pendulum Angle", fontsize=12)
axes[0, 1].legend(fontsize=9)
axes[0, 1].grid(True, alpha=0.3)

# 制御入力
axes[1, 0].axhline(0, color="black", linestyle="--", alpha=0.4)
axes[1, 0].set_xlabel("Time [s]", fontsize=11)
axes[1, 0].set_ylabel("Control force [N]", fontsize=11)
axes[1, 0].set_title("Control Input", fontsize=12)
axes[1, 0].legend(fontsize=9)
axes[1, 0].grid(True, alpha=0.3)

# 極配置の比較
ax = axes[1, 1]
for label, Q, R in configs:
    P = solve_continuous_are(A, B, Q, R)
    K = np.linalg.inv(R) @ B.T @ P
    A_cl = A - B @ K
    poles = np.linalg.eigvals(A_cl)
    ax.scatter(poles.real, poles.imag, s=80, label=label, zorder=5)

ax.axvline(0, color="black", linewidth=0.8, alpha=0.5)
ax.axhline(0, color="black", linewidth=0.8, alpha=0.5)
ax.axvspan(-30, 0, alpha=0.03, color="green")
ax.set_xlabel("Real", fontsize=11)
ax.set_ylabel("Imaginary", fontsize=11)
ax.set_title("Closed-loop Poles (LQR)", fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("lqr_control.png", dpi=150, bbox_inches="tight")
plt.show()

この結果から、LQR制御の特性と重み行列の影響が明確に理解できます。

カート位置（左上）: Aggressive設定（$\bm{Q}$ が大きく $\bm{R}$ が小さい）ではカートが約1秒で原点に戻り、その後ほぼ静止しています。一方、Light control設定では収束に約4秒かかっています。Balanced設定はその中間で、約2秒で収束しています。初期位置 $x_0 = 0.5$ m からの応答を見ると、$\bm{Q}/\bm{R}$ の比が10倍変わるごとに応答速度が概ね2倍変わることがわかります
振子角度（右上）: 初期角度 $\theta_0 = 0.1$ rad（約5.7度）からの安定化が3つの設定で比較できます。Aggressive設定では角度が0.3秒以内にゼロに戻りますが、Light control設定では約3秒かかっています。$\bm{Q}$ の $\theta$ 成分を大きくすることで角度の偏差に高いペナルティを課し、振子の安定化を最優先にできます。倒立振子は不安定システムであるため、角度の制御は特に重要です
制御入力（左下）: 応答が速いほど大きな制御力が必要です。Aggressive設定では初期に約40 Nの力が必要ですが、Light control設定では約5 N程度です。$\bm{R}$ を大きくすることで入力の大きさを制限できますが、応答は遅くなります。これがLQRの本質的なトレードオフです。実システムではアクチュエータ（モーター）の出力限界があるため、$\bm{R}$ の下限は物理的制約から決まります
極配置（右下）: 各設定での閉ループ極の位置が比較できます。すべての極が左半面（安定領域）に配置されており、LQRが安定性を保証していることがわかります。Aggressive設定ほど極が左半面の遠くに配置され、速い応答に対応しています。重要なのは、これらの極配置を設計者が直接指定したのではなく、LQRが重み行列から自動的に最適な位置を決定しているという点です。極配置法では4つの極の位置を試行錯誤で決める必要がありましたが、LQRでは「何を重視するか」を $\bm{Q}$ と $\bm{R}$ で指定するだけで済みます
全体的な考察: 倒立振子は開ループで不安定なシステム（正の固有値を持つ）ですが、LQR制御により安定化されています。3つの設定を比較することで、重み行列が制御性能に与える影響を系統的に理解できます。実際の設計では、まずBrysonの法則で初期値を設定し、シミュレーションを見ながら微調整するという手順が一般的です

次に、LQRの理論をより深く理解するため、手計算で追えるシンプルな数値計算例を見ていきましょう。

数値計算例

例1：1入力1状態系のLQR

最も簡単なLQR問題として、スカラー系 $\dot{x} = ax + bu$ を考えます。ここで $a = 1$（不安定な系）、$b = 1$ とします。

コスト関数は $J = \int_0^\infty (qx^2 + ru^2) \, dt$ で、$q = 10$、$r = 1$ と設定します。

スカラーの場合、リカッチ方程式は次の2次方程式になります。

$$ 2aP – \frac{b^2}{r}P^2 + q = 0 $$

$$ 2 \times 1 \times P – \frac{1}{1}P^2 + 10 = 0 $$

$$ P^2 – 2P – 10 = 0 $$

解の公式より $P = \frac{2 + \sqrt{4 + 40}}{2} = 1 + \sqrt{11} \approx 4.317$（正の解を採用）。

最適ゲインは $K = \frac{b}{r}P = 4.317$、最適制御則は $u = -4.317x$ です。

閉ループの極は $a – bK = 1 – 4.317 = -3.317$ で、安定かつ速い応答です。

初期状態 $x(0) = 1$ のとき、最適コストは $J^* = Px(0)^2 = 4.317$ です。

例2：$Q/R$ 比の影響の定量的分析

上のスカラー系で $Q/R$ の比を変えたとき、最適ゲインと閉ループ極がどう変化するかを系統的に調べましょう。

$q/r$	$P$	$K$	閉ループ極	$J^*$ ($x_0 = 1$)
0.1	1.105	1.105	$-0.105$	1.105
1	1.414	1.414	$-0.414$	1.414
10	4.317	4.317	$-3.317$	4.317
100	11.05	11.05	$-10.05$	11.05
1000	32.62	32.62	$-31.62$	32.62

$q/r$ が10倍になるごとに、ゲインは約 $\sqrt{10} \approx 3.16$ 倍に増加し、閉ループ極は約 $\sqrt{10}$ 倍左に移動します。これは $q/r \gg 1$ の極限で $K \approx \sqrt{q/r}$、閉ループ極 $\approx -\sqrt{q/r}$ となることから理解できます。

最適コスト $J^*$ も $q/r$ の増加とともに増加しますが、$q/r$ の比ほどは増えません。これは、制御入力を強くすることで状態偏差が大幅に減少し、$qx^2$ の項が小さくなるためです。

例3：2状態系のリカッチ方程式

次に、2次元システムの例を考えます。バネ-マス-ダンパ系を状態空間モデルで表します。

$$ \bm{A} = \begin{bmatrix} 0 & 1 \\ -\omega_n^2 & -2\zeta\omega_n \end{bmatrix}, \quad \bm{B} = \begin{bmatrix} 0 \\ 1 \end{bmatrix} $$

$\omega_n = 2$ rad/s（固有角振動数）、$\zeta = 0.1$（減衰比、かなり振動的）として、$\bm{Q} = \text{diag}(10, 1)$、$\bm{R} = [0.1]$ とします。

リカッチ方程式を解くと、$\bm{P} = \begin{bmatrix} 3.287 & 0.497 \\ 0.497 & 0.274 \end{bmatrix}$ が得られます。

最適ゲインは $\bm{K} = \bm{R}^{-1}\bm{B}^T\bm{P} = [4.97, \; 2.74]$ です。

閉ループ行列 $\bm{A} – \bm{B}\bm{K} = \begin{bmatrix} 0 & 1 \\ -8.97 & -3.14 \end{bmatrix}$ の固有値は $s = -1.57 \pm 2.68j$ です。

元のシステムの固有値は $s = -0.2 \pm 1.99j$（振動的で減衰が遅い）でしたが、LQR制御により減衰が大幅に改善されています。減衰比は $\zeta_{\text{cl}} = 1.57 / \sqrt{1.57^2 + 2.68^2} = 0.505$ となり、元の $\zeta = 0.1$ から約5倍に改善されました。

例4：LQRと極配置法の比較

同じ2状態系で極配置法とLQRの結果を比較してみましょう。

極配置法で閉ループ極を $s = -2 \pm 2j$ に設定すると、ゲインは $\bm{K}_{\text{pp}} = [4.0, \; 3.6]$ が得られます（アッカーマンの公式による）。

LQRでは $\bm{K}_{\text{LQR}} = [4.97, \; 2.74]$ でした。

両者のゲインは近い値ですが、違いがあります。LQRの方が位置ゲインは大きく、速度ゲインは小さくなっています。これは、LQRが $\bm{Q}$ の設定（位置の重みが10、速度の重みが1）に基づいて、位置の制御を優先する設計を自動的に行っているためです。

ロバスト性の面では、LQRが自動的に安定余裕を保証するのに対し、極配置法ではゲイン余裕や位相余裕を別途確認する必要があります。

LQRの拡張と限界

離散時間LQR

連続時間LQRの離散時間版も広く使われます。システム $\bm{x}_{k+1} = \bm{A}_d\bm{x}_k + \bm{B}_d\bm{u}_k$ に対し、コスト関数は

$$ J = \sum_{k=0}^{\infty} (\bm{x}_k^T\bm{Q}\bm{x}_k + \bm{u}_k^T\bm{R}\bm{u}_k) $$

最適ゲインは $\bm{K} = (\bm{R} + \bm{B}_d^T\bm{P}\bm{B}_d)^{-1}\bm{B}_d^T\bm{P}\bm{A}_d$ で、$\bm{P}$ は離散代数リカッチ方程式（DARE）の解です。連続時間版との違いは、$\bm{R}^{-1}$ の代わりに $(\bm{R} + \bm{B}_d^T\bm{P}\bm{B}_d)^{-1}$ が現れる点です。

有限時間LQR

上で扱った無限時間LQR（定常問題）に対し、有限時間区間 $[0, T]$ での最適制御を考えることもできます。この場合、リカッチ方程式は微分方程式（リカッチ微分方程式）になり、$\bm{P}(t)$ が時間とともに変化します。終端条件 $\bm{P}(T) = \bm{S}$（終端重み行列）から逆方向に積分して解きます。$T \to \infty$ で定常解に収束するため、無限時間LQRは有限時間LQRの極限と見なせます。

LQRの限界

LQR制御は強力な手法ですが、以下の限界があります。

全状態フィードバックの仮定: すべての状態変数を直接測定できることを前提としています。実際にはカルマンフィルタなどの状態推定器が必要で、LQG制御として統合されます
線形性の仮定: 非線形システムに対しては、動作点周りで線形化して適用するか、非線形拡張（SDRE法など）を用います
MIMO系でのロバスト性: SISO系では安定余裕が保証されますが、MIMO系ではこの保証は一般に成り立ちません。入力チャネルを1つずつ切断した場合のロバスト性は保証されますが、複数チャネルの同時変動に対しては脆弱になり得ます

これらの限界を踏まえた上で、LQRは最適制御理論の出発点であり、LQG制御、$H_\infty$ 制御、モデル予測制御（MPC）など、より高度な制御手法の基盤となっています。

まとめ

本記事では、LQR制御の理論を基礎から応用まで体系的に解説しました。要点を振り返ります。

LQRは二次コスト関数 $J = \int_0^\infty(\bm{x}^T\bm{Q}\bm{x} + \bm{u}^T\bm{R}\bm{u})dt$ を最小化する最適状態フィードバックゲインを求める手法です。「性能」と「制御コスト」のトレードオフを数学的に厳密に最適化します
最適ゲインは $\bm{K} = \bm{R}^{-1}\bm{B}^T\bm{P}$ で求められ、$\bm{P}$ は連続代数リカッチ方程式（CARE）の一意の正定値解です。HJB方程式から導出されるこの結果は、複雑な最適化問題を1つの代数方程式に帰着させます
重み行列 $\bm{Q}$ と $\bm{R}$ の相対比が応答速度と制御入力のトレードオフを決定します。Brysonの法則は実用的な初期設計法として広く使われています
LQR制御はSISO系でゲイン余裕 $\geq 6$ dB、位相余裕 $\geq 60°$ を自動的に保証するロバスト性を持ち、極配置法にはないこの性質がLQRの大きな利点です
最適コストは $J^* = \bm{x}(0)^T\bm{P}\bm{x}(0)$ で初期状態の二次形式として表され、$\bm{P}$ が各方向のコスト感度を表します
LQRの限界（全状態フィードバックの仮定、線形性、MIMO系のロバスト性）を理解することで、LQG制御や $H_\infty$ 制御などへの発展の動機が自然に理解できます

次のステップとして、以下の記事も参考にしてください。

LQG制御 — LQRとカルマンフィルタの統合
カルマンフィルタと最適推定 — 最適オブザーバ
状態フィードバック制御と極配置法 — LQRとの対比