Zenn
🙆‍♀️

モデル予測制御（MPC）の基本

に公開
Python
数理最適化
制御工学
モデル予測制御
tech

モデル予測制御(MPC)の考え方

はじめに

最近、技術系のメディアやコミュニティで「MCP」というキーワードを目にする機会が増えましたね。MCPという三文字は今や、Model Context Protocolのことを指すというのが普通になってきました。これはLLMアプリケーションが外部サービスとの連携を標準化するためのプロトコルということで、ソフトウェアアーキテクチャの観点からも非常に興味深い動きだと感じています。

一方で、自分自身はMCPをゴリゴリ扱ってLLMサービスを拡充する機会というのは多くありません。仕事的にも興味の範囲的にも、AI分野を制御と絡めて議論するようなケースが非常に多く、MCPをMPCと空目する機会が増えてまいりました。もちろん両者は直接関係ありませんが、この偶然を機に、現代制御理論の重要な一角を占めるMPCのコンセプトについて、改めて自身の理解を整理しつつ、ご紹介・宣伝したいと思い、久々に記事を書いてみます。

本記事では、MPCの基本的な考え方から、その定式化、さらには簡単な実装例まで、制御工学のバックグラウンドがない方にもそのエッセンスが伝わるよう、順を追って解説していくことを目指します。少々長丁場になるかもしれませんが、お付き合いいただけますと幸いです。

モデル予測制御（MPC）って、ざっくり言うと？

モデル予測制御（MPC）の核心を一言で表現するならば、「システムの未来の挙動をモデルに基づいて予測し、その予測結果に基づいて現在の最適な制御入力を決定する」というフィードバック制御戦略、と言えるでしょうか。これは「Receding Horizon Control（後退ホライズン制御）」とも呼ばれる考え方に基づいています。

もう少し具体的に、MPCが内部で行っているプロセスを見てみましょう。それは、非常に短い時間ステップごとに、以下のサイクルを繰り返すことに相当します。

  1. システムの数学モデルを用いて、現在から有限な未来（予測ホライズン）までの状態遷移を予測する。
    • ここでは、考えられる制御入力のシーケンスを適用した場合に、システムが将来どのように振る舞うかをシミュレーションします。
  2. 予測された未来の軌道の中から、事前に定義された評価関数（目的関数）を最も良くする（通常は最小化する）制御入力シーケンスを「最適化計算」によって探索する。
    • 評価関数には、目標値からの逸脱ペナルティや制御入力のエネルギー消費などが含まれます。
    • 重要な点として、アクチュエータの物理的限界や安全上の要求といった制約条件（例：入力の上限・下限、状態量の許容範囲）を満たす範囲内で最適な解を探します。これがMPCの際立った利点の一つです。
  3. 最適化によって得られた制御入力シーケンスのうち、最初のステップの入力値のみを、実際にシステムに適用する。
  4. 時間が1ステップ進んだ後、最新のシステム状態に基づいて、再びステップ1からの予測と最適化のプロセスを繰り返す。

このように、常に最新の状態情報に基づいて未来を見通し、計画を立て直していくアプローチにより、モデルの不確かさや外乱の影響が存在する現実のシステムに対しても、ロバストかつ高性能な制御を実現することが期待されます。

もう少し詳しく：MPCの仕組み（数式と共に）

さて、ここからは少しだけ数式を交えて、MPCが内部でどんな計算をしているのか見ていきましょう。制御理論に慣れていないと少し難しく感じるかもしれませんが、基本的なアイデアは先ほどの説明と同じです。

一般的な定式化

まず、制御対象のシステムの振る舞いを数式で表した「モデル」が必要です。離散時間システムを考えると、現在の状態 x_k と制御入力 u_k から、次の時刻の状態 x_{k+1} が決まるとします。これは以下のような式で表現できます。

x_{k+1} = f(x_k, u_k)

ここで、x_k は時刻 k でのシステムの状態（位置、速度、温度など）、u_k は時刻 k でシステムに加える操作（モーターへの電圧、バルブの開度など）、f はシステムの状態がどう変化するかを表す関数（システムダイナミクス）です。

次に、MPCでは「未来予測」を行います。現在の時刻 k から N ステップ先までの未来（これを予測ホライズンと呼びます）の状態と制御入力を考えます。そして、この未来予測期間における「制御の良さ」を評価するための評価関数 J を定義します。一般的には以下のような形をしています。

J = \sum_{i=0}^{N-1} L(x_{k+i|k}, u_{k+i|k}) + \Phi(x_{k+N|k})

ここで、x_{k+i|k} は時刻 k で予測した i ステップ先の状態、u_{k+i|k} は同じく時刻 k で計画する i ステップ先の制御入力です。

  • L(x, u) は各時刻での「コスト」を表す関数（ステージコスト）で、目標値からのずれや制御入力の大きさなどを評価します。これを小さくしたい。
  • \Phi(x) は予測ホライズンの終端 N ステップ目での状態に対するコスト（終端コスト）で、最終的にシステムを安定な状態に導くために設定されたりします。これも小さくしたい。

さらに、現実のシステムには制約条件があります。例えば、制御入力 u の大きさには上限と下限があったり（u_{min} \le u_k \le u_{max}）、システムの状態 x も特定の範囲内に収めなければならなかったりします（x_{min} \le x_k \le x_{max}）。これらの制約も考慮に入れる必要があります。

MPCが各時刻 k で行うことは、これらのモデル、評価関数、制約条件のもとで、評価関数 J を最小にするような未来の制御入力のシーケンス \{ u_{k|k}, u_{k+1|k}, \dots, u_{k+N-1|k} \} を見つけ出すことです。数式で書くと、以下のような最適化問題を解くことになります。

\begin{aligned} \min_{u_{k|k}, \dots, u_{k+N-1|k}} \quad & \sum_{i=0}^{N-1} L(x_{k+i|k}, u_{k+i|k}) + \Phi(x_{k+N|k}) \\ \text{s.t.} \quad & x_{k+i+1|k} = f(x_{k+i|k}, u_{k+i|k}), && i=0, \dots, N-1 \\ & x_{k|k} = x_{current} && \\ & u_{min} \le u_{k+i|k} \le u_{max}, && i=0, \dots, N-1 \\ & x_{min} \le x_{k+i|k} \le x_{max}, && i=1, \dots, N \\ & \dots && \text{(その他の制約)} \end{aligned}

この最適化問題を解いて得られた未来の制御入力シーケンスのうち、最初の制御入力 u_{k|k} のみを実際にシステムに適用します。そして、次の時刻 k+1 になったら、新しい状態 x_{k+1} を観測（または推定）し、再び同じ手順で N ステップ先までを予測・最適化する、という流れを繰り返します。

特殊なケース：線形MPC（Linear MPC）

さて、上記の一般的な定式化は、関数 f や L, \Phi が複雑な形をしていると、最適化問題を解くのが非常に難しくなることがあります。計算時間がかかりすぎたり、そもそも解が見つからなかったり…。

ここからが少し注意が必要なのですが、これはあくまでMPCの中でも特殊で幸運なケースです。 もし、扱いたいシステムが線形で、評価関数が二次形式、制約条件も線形で書ける場合、話は少し簡単になります。

  • 線形システムモデル: x_{k+1} = Ax_k + Bu_k
    • A, B は定数行列です。
  • 二次形式の評価関数: J = \sum_{i=0}^{N-1} (x_{k+i|k}^T Q x_{k+i|k} + u_{k+i|k}^T R u_{k+i|k}) + x_{k+N|k}^T P x_{k+N|k}
    • Q, R, P は重みを表す半正定値行列（通常は対角行列などで設定）です。
  • 線形不等式制約: C x_{k+i|k} + D u_{k+i|k} \le e （入力制約 u_{min} \le u_{k+i|k} \le u_{max} などもこの形式で書けます）

これらの条件が満たされる場合、先ほどの最適化問題は**二次計画問題（Quadratic Programming, QP）**と呼ばれる、比較的よく研究されていて解きやすい（効率的なアルゴリズムが存在する）クラスの最適化問題に帰着できます。

QP問題は、一般的に以下のような形式で表されます。

\begin{aligned} \min_{z} \quad & \frac{1}{2} z^T H z + g^T z \\ \text{s.t.} \quad & L \le A_{cons} z \le U \end{aligned}

線形MPCの場合、最適化変数 z は未来の制御入力列 \{ u_{k|k}, \dots, u_{k+N-1|k} \} などに対応し、行列 H, A_{cons} やベクトル g, L, U は、システムのモデル(A, B)、評価関数の重み(Q, R, P)、制約(C, D, e)、そして現在の状態 x_k から計算できます。（詳細な変換は少し煩雑なのでここでは省略しますが、興味のある方は調べてみてください。）

繰り返しになりますが、このようにQP問題として解けるのは線形MPCという特殊な場合です。 世の中の多くのシステムは非線形なので、一般的なMPCではもっと複雑な最適化計算が必要になることが多い、という点は覚えておく必要があると思います。

Pythonで動かしてみる：簡単な線形MPCの例

理論だけだとイメージが湧きにくいかもしれませんので、簡単な線形システムに対してPythonで線形MPCを実装してみましょう。ここでは、状態が2次元の減衰振動系を考えます。これは、質量 m=1、バネ定数 k=0.75、ダンピング係数 c=0.25 の物体に力（入力 u、今回は加速度指令と解釈）を加えたときの、位置 p と速度 v の変化を表すモデル m \ddot{p} + c \dot{p} + k p = u です。

状態 x = [p, v]^T として離散化（時間刻み dt=0.1、前進オイラー法）すると、システムモデルは以下のようになります。
x_{k+1} = A x_k + B u_k
A = \begin{bmatrix} 1 & dt \\ -k \cdot dt & 1 - c \cdot dt \end{bmatrix} = \begin{bmatrix} 1 & 0.1 \\ -0.075 & 0.975 \end{bmatrix}
B = \begin{bmatrix} 0 \\ dt \end{bmatrix} = \begin{bmatrix} 0 \\ 0.1 \end{bmatrix}

目標: 状態 x を初期値 x_0 = [-2.0, 0.0]^T （位置-2、速度0）から目標値 x_{ref} = [5.0, 0.0]^T （位置5、速度0）に素早く整定させたい。
制約: 制御入力 u は -5 から 5 の範囲内 (-5 \le u_k \le 5)。

評価関数は J = \sum_{i=0}^{N-1} \left( (x_{k+i|k}-x_{ref})^T Q (x_{k+i|k}-x_{ref}) + R u_{k+i|k}^2 \right) + (x_{k+N|k}-x_{ref})^T P (x_{k+N|k}-x_{ref}) とします。予測ホライズン N=20 とします。

import numpy as np
import cvxpy as cp
import matplotlib.pyplot as plt

# システムパラメータ (減衰振動系, dt=0.1, m=1, k=0.75, c=0.25)
dt = 0.1
k_spring = 0.75  # バネ定数
c_damper = 0.25  # ダンピング係数
A = np.array([[1, dt],
              [-k_spring * dt, 1 - c_damper * dt]])
B = np.array([[0],
              [dt]]) # 入力は加速度指令 u
n = A.shape[0]  # 状態次元 (p, v)
m = B.shape[1]  # 入力次元 (u)

# MPCパラメータ
N = 20  # 予測ホライズン
Q = np.diag([1.0, 0.1])  # 状態コスト重み (位置誤差を重視)
R = np.diag([0.01])      # 入力コスト重み
P = Q                   # 終端コスト重み (ここではQと同じにする)

u_max = 5.0
u_min = -5.0

# 目標値
x_ref = np.array([5.0, 0.0]) # 位置 5, 速度 0

# シミュレーション設定
T = 100  # シミュレーションステップ数 (10秒間)
x0 = np.array([-2.0, 0.0])  # 初期状態 (位置 -2, 初速 0)

# 状態と入力の履歴を保存する配列
x_hist = np.zeros((n, T + 1))
u_hist = np.zeros((m, T))
x_hist[:, 0] = x0

# --- MPC計算関数 ---
def solve_mpc(x_current, x_ref):
    # 最適化変数を定義 (未来の入力列と状態列)
    U = cp.Variable((m, N))
    X = cp.Variable((n, N + 1))
    
    # 目標値ベクトルをN+1個並べたものを作成 (cvxpyのブロードキャストを助ける)
    # x_ref_seq = np.tile(x_ref.reshape(-1, 1), (1, N + 1))
    # ↑ cvxpy は numpy array との演算でブロードキャスト可能なので不要

    # コスト関数を初期化
    cost = 0
    # 制約リストを初期化
    constraints = [X[:, 0] == x_current] # 初期状態制約
    
    # 予測ホライズンにわたってコストと制約を定義
    for t in range(N):
        # 状態コスト (目標値からの偏差の二次形式)
        cost += cp.quad_form(X[:, t] - x_ref, Q)
        # 入力コスト
        cost += cp.quad_form(U[:, t], R)
        
        # システムダイナミクス制約
        constraints += [X[:, t + 1] == A @ X[:, t] + B @ U[:, t]]
        
        # 入力制約
        constraints += [u_min <= U[:, t], U[:, t] <= u_max]
        
    # 終端コスト (目標値からの偏差の二次形式)
    cost += cp.quad_form(X[:, N] - x_ref, P)
        
    # 最適化問題の定義
    problem = cp.Problem(cp.Minimize(cost), constraints)
    
    # 問題を解く (QPソルバー OSQP を使用)
    problem.solve(solver=cp.OSQP, verbose=False)
    
    # 解が見つかったかチェック
    if problem.status == cp.OPTIMAL or problem.status == cp.OPTIMAL_INACCURATE:
        # 最初の制御入力を返す
        return U[:, 0].value
    else:
        # 解けなかった場合
        print(f"時刻 {k}: 最適化問題が解けませんでした。status: {problem.status}")
        return np.zeros(m) # フォールバックとしてゼロ入力

# --- シミュレーションループ ---
xk = x0
for k in range(T):
    # 現在の状態 xk と目標値 x_ref に基づいてMPC問題を解く
    uk = solve_mpc(xk, x_ref)
    
    # uk が None でないか確認 (ソルバーが失敗した場合など)
    if uk is None:
        print(f"時刻 {k} でソルバーが解を返しませんでした。シミュレーションを停止します。")
        # 実際の履歴に合わせて T を調整
        T = k 
        x_hist = x_hist[:, :T+1]
        u_hist = u_hist[:, :T]
        break
    
    # 履歴に保存
    u_hist[:, k] = uk
    
    # システムを1ステップ進める
    xk_plus1 = A @ xk + B @ uk
    x_hist[:, k + 1] = xk_plus1
    
    # 状態を更新
    xk = xk_plus1

# --- 結果のプロット ---
plt.figure(figsize=(10, 8))

# 状態 (位置と速度) のプロット
plt.subplot(2, 1, 1)
plt.plot(np.arange(T + 1) * dt, x_hist[0, :], '-o', markersize=3, label='Position p (State x1)')
plt.plot(np.arange(T + 1) * dt, x_hist[1, :], '-o', markersize=3, label='Velocity v (State x2)')
plt.axhline(x_ref[0], color='r', linestyle='--', label='Target Position')
plt.axhline(x_ref[1], color='g', linestyle='--', label='Target Velocity')
plt.ylabel('States')
plt.title('Linear MPC Simulation (Damped Oscillator)')
plt.grid(True)
plt.legend()

# 入力のプロット
plt.subplot(2, 1, 2)
plt.step(np.arange(T) * dt, u_hist[0, :], '-o', where='post', markersize=3, label='Input u')
plt.plot(np.arange(T) * dt, np.ones(T) * u_max, 'r--', label='u_max')
plt.plot(np.arange(T) * dt, np.ones(T) * u_min, 'r--', label='u_min')
plt.ylabel('Input')
plt.xlabel('Time [s]')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

このコードを実行すると、初期状態で位置-2、速度0だったシステムが、入力制約を守りながら目標位置5に振動を抑えつつ近づき、最終的に速度が0になるように制御される様子が確認できるはずです。システムの固有の振動特性に対し、MPCが予測に基づいて制御入力を調整し、目標達成を目指す様子が見て取れるかと思います。（定常応答においては入力値がばね係数にしたがった引っ張りに抵抗する形で力を釣り合わせていることと、過渡応答に関してはオーバーシュートすることなく上手に力を与えていることを見てください）

あくまで簡単な例ですが、MPCが「制約を守りながら」「未来を予測して」「最適化計算によって」制御入力を決定している雰囲気が、少しでも伝われば幸いです。その他、パラメータをいじったり、目標値を一定値ではなく時間変化させてみたり、試すと面白いと思います。

少し補足：この例の単純化について

ちなみに、今回のシミュレーションは非常にシンプルな状況設定であることにも注意が必要です。特に重要な点として、

  1. モデル化誤差がない: MPCが予測に使う内部モデル（x_{k+1} = 0.9 x_k + 0.5 u_k）と、シミュレーション対象のシステム（xk_plus1 = A @ xk + B @ uk の部分）が完全に一致しています。現実には、モデルは必ずしも完璧ではなく、モデル化誤差が存在します。MPCの性能は内部モデルの精度に大きく依存するため、対象システムを正確にモデル化するシステム同定の技術が重要になってきます。一方で正確性を重んじると非線形モデルの作りこみが作業という泥沼にハマってしまいます。後に話しますが、非線形予測モデルを利用すると上記のコード例ほど単純には問題を解けません。モデル化の複雑さと計算量にはトレードオフが存在することを意識しなければなりません。また、場合によっては確率的な挙動を考慮した統計モデルを用いることもあります。

  2. 状態が直接観測できる: シミュレーションループ内で、次の時刻の状態 xk_plus1 が計算後すぐに xk として次のMPC計算に利用できています。しかし、実システムでは、システムの状態（例えば内部温度や化学反応の進行度など）を直接、正確に、遅れなく測定できるとは限りません。測定できない状態変数は、センサー情報などから**状態推定器（オブザーバ）**を用いて推定する必要が生じます。カルマンフィルタなどがその代表例ですね。

これらの要素（モデル化誤差、状態推定）を考慮に入れると、MPCシステムの設計はさらに複雑で奥深いものになります。今回は、あくまでMPCの基本的な計算の流れを理解するための、理想的なケーススタディと捉えていただければと思います。

もっと一般的なMPCへ：非線形MPCとその解法

さて、線形MPCはQPという比較的扱いやすい問題に帰着できましたが、現実の多くのシステムは非線形性を持っています。システムモデル f(x, u) が非線形であったり、評価関数 L, \Phi が二次形式でなかったりする場合、MPCは各時刻で一般的な非線形計画法（Nonlinear Programming, NLP）問題を解く必要が出てきます。これが非線形MPC（NMPC）**です。

NLP問題は、目的関数や制約条件が非線形であるため、一般的に解くのが難しくなります。大域的最適解を見つける保証はなく、局所解に陥る可能性もありますし、計算コストも高くなりがちです。NMPCを解くためのアプローチやアルゴリズムは多岐にわたりますが、ここでは代表的な考え方をいくつかご紹介したいと思います。

1. NMPC問題の定式化と解法アプローチ

NMPCで解くべき最適化問題は、基本的には以下の形式をしています。

\begin{aligned} \min_{U} \quad & J(X, U) \\ \text{s.t.} \quad & x_{k+1} = f(x_k, u_k), && k=0, \dots, N-1 \\ & g(x_k, u_k) = 0, && k=0, \dots, N-1 \\ & h(x_k, u_k) \le 0, && k=0, \dots, N-1 \end{aligned}

ここで U = \{u_0, ..., u_{N-1}\} は制御入力シーケンス、X = \{x_0, ..., x_N\} はそれによって決まる状態シーケンスです。この問題を解くためのアプローチとして、大きく以下のような方向性があります。

  • 制約付き最適化問題として直接解く:
    最適化アルゴリズムの中で、制約条件を直接扱いながら解を探索する方法です。後述する逐次二次計画法（SQP）や内点法（IPM）などがこのカテゴリに含まれます。各反復ステップで問題を線形化したり、制約を満たすように探索方向を調整したりします。

  • 制約なし最適化問題への変換:
    制約条件を目的関数にペナルティ項として組み込むことで、見かけ上、制約なしの問題に変換して解くアプローチです。単純なペナルティ法では数値的な問題が生じやすいため、 拡張ラグランジュ法（Augmented Lagrangian Method） などがよく用いられます。これはラグランジュ関数に二次のペナルティ項を追加するもので、より安定して元の問題の解に収束させることが期待できます。

  • KKT条件に基づく解法:
    最適性の必要条件であるカルーシュ・キューン・タッカー（KKT）条件を直接解くアプローチです。KKT条件は、ラグランジュ関数に関する勾配条件、元の問題の制約条件、そして相補性条件からなる連立方程式（と不等式）になります。この（しばしば大規模な）非線形連立方程式をニュートン法などで解くことになります。この際、ニュートン法の各ステップで現れる線形連立方程式を解くために、後述するGMRESなどの反復解法が利用されることがあります（特にC/GMRESとして知られる手法）。

  • 関数近似アプローチ（コロケーション法など）:
    これは少し異なる視点で、連続時間系のNMPCなどでよく使われます。状態変数や制御入力の時間軌道を、適当な基底関数（例えば多項式）の線形結合で近似します。そして、システムの微分方程式や制約条件を、いくつかの代表的な時刻（コロケーション点）でのみ満たすように要求します。これにより、無限次元の最適制御問題が、有限次元のNLP問題に変換され、それをSQPやIPMなどで解く、という流れです。直接多重シューティング法なども関連するアプローチですね。

どのアプローチが良いかは、問題の構造やサイズ、要求される計算速度などによって変わってきます。意外と奥が深い世界です。

2. 代表的な最適化アルゴリズム

上記のようなアプローチの中で、具体的な最適化計算を実行するために、以下のようなアルゴリズムがよく利用されます。

  • 逐次二次計画法（Sequential Quadratic Programming, SQP）:
    NLP問題を、各反復ステップでQP（二次計画問題）の部分問題に近似して解く方法です。ラグランジュ関数の二次近似と制約の線形近似に基づいてQP部分問題を構築し、その解を用いて探索点を更新していきます。制約が有効（アクティブ）になる状況をうまく扱える傾向があると言われています。

  • 内点法（Interior Point Method, IPM）:
    不等式制約 h(x, u) \le 0 を、バリア関数（例えば -\mu \log(-h(x,u)) ）を用いて目的関数に組み込み、制約の境界内部（Interior）を探索しながら最適解に近づいていく方法です。バリアパラメータ \mu を徐々にゼロに近づけていきます。大規模な問題に対して効率が良い場合が多いとされ、近年のNLPソルバーで広く採用されています。

  • アクティブセット法（Active Set Method）:
    最適解において有効（等号成立）になっていると予測される不等式制約の集合（アクティブセット）を管理しながら解を探索する手法です。各反復で、アクティブセットに関する等式制約付きの（より簡単な）問題を解き、ラグランジュ乗数などを評価してアクティブセットを更新します。SQPのQP部分問題を解く際などに用いられることがあります。

  • C/GMRES（Constrained / Generalized Minimum Residual Method）:
    先述のKKT条件に基づく解法などで現れる、大規模な線形連立方程式を解くための反復解法の一つです。特に、最適制御問題から導かれるKKTシステムを効率的に解くために、GMRES法を制約付き問題に適用する形で使われることがあります。

これらのアルゴリズムも、それぞれに得意な問題の性質や、収束特性、計算コストなどが異なります。NMPCを実装する際には、問題の特性を見極めて適切なアルゴリズム（を実装したソルバー）を選択することが重要になってきますね。

もっと一般的なMPCへ：最適化ソルバーの話

さて、前のセクションではNMPC問題を解くための様々なアプローチ（拡張ラグランジュ法やKKT条件に基づく方法など）や、具体的な最適化アルゴリズム（SQP、IPMなど）を見てきました。これらのアルゴリズムは非常に強力ですが、実際にこれらをゼロから自分で実装するのは、多くの場合、現実的ではありません。

幸いなことに、これらの複雑な最適化アルゴリズムを実装したソフトウェアパッケージ、すなわち最適化ソルバーが多数開発されています。NMPCのように、一般的に非線形計画法（NLP）の問題を解く必要がある場合、我々はそのためのNLPソルバーを利用することになります。

ただし、一口にNLPソルバーと言っても、実装されているアルゴリズム（SQP系かIPM系か、など）や、得意とする問題の規模、計算速度、ライセンス形態（オープンソースか商用か）などが異なります。そのため、目的や状況に応じて適切なソルバーを選択することが重要になってきます。

代表的なNLPソルバーとしては、以下のようなものがあります。

  • IPOPT: オープンソースで広く使われている内点法のソルバー。比較的安定していると言われています。
  • SNOPT: 大規模な問題を効率的に解くことができる逐次二次計画法（SQP）に基づいたソルバー。商用ライセンスが必要です。
  • WORHP: ドイツ航空宇宙センター(DLR)で開発されたSQP系のソルバー。
  • NLopt: 様々なアルゴリズム（勾配ベース、勾配フリーなど）を統一的なインターフェースで提供するライブラリ。

また、NMPC問題を効率的に定式化し、これらのソルバーを呼び出すためのフレームワークも存在します。例えば、

  • CasADi: Pythonなどから利用できる強力な自動微分・最適化モデリングツール。様々なNLPソルバーと連携できます。
  • ACADOS: 高速なNMPC実装に特化したツールキット。リアルタイム制御用途などで使われます。

どのソルバーやツールを使うかは、対象とする問題の特性（非線形性の強さ、問題規模、要求される計算速度など）によって慎重に選択する必要があります。このあたりも、MPCを実用する上では意外と重要なノウハウだったりしますね。このソルバー周りの話も、また別の機会に詳しくまとめてみても面白いかもしれません。

実時間計算の難しさと代替アプローチ

ここまで見てきたように、非線形MPC（NMPC）は非常に強力な制御手法ですが、各制御ステップで非線形計画法（NLP）問題を解く必要があるという点は、実用上の大きな課題となり得ます。特に、ロボット制御や自動運転のように、ミリ秒単位での高速な制御周期が要求されるアプリケーションでは、NLPの最適化計算が時間内に完了しない、あるいは安定して収束しないという問題に直面することがあります。

このリアルタイム計算の難しさに対処するために、あるいは異なる視点から同等の制御性能を目指すために、様々な代替アプローチが研究・開発されています。いくつか例を挙げると、

  • 模倣学習（Imitation Learning）ベースのアプローチ:
    あらかじめオフラインでNMPC（あるいは他の専門家による制御）を実行し、その時の状態と最適な制御入力のペアを大量に収集します。そして、そのデータセットを使って、状態から直接最適な制御入力を出力するような関数（例えばニューラルネットワーク）を学習させます。オンラインでは、この学習済み関数を使って制御入力を高速に計算するため、最適化計算は不要になります。ただし、学習データの範囲外の状況に対する性能（汎化性能）が課題となることがあります。

  • 逐次線形化（Successive Linearization）アプローチ:
    非線形のシステムモデルを、現在の状態や予測軌道の周りで逐次的に線形近似し、線形MPC（つまりQP問題）を繰り返し解くというアプローチです。非線形性を完全に捉えることはできませんが、QPソルバーはNLPソルバーよりも一般的に高速かつ安定しているため、計算時間を大幅に削減できる可能性があります。近似精度と計算速度のトレードオフを考慮する必要がありますね。

  • 強化学習（Reinforcement Learning, RL）との連携・融合:
    実は、MPCが解いている最適化問題は、強化学習におけるベルマン方程式を満たすような最適方策（policy）を見つける問題と、本質的に深く関連しています。RLでは、試行錯誤を通じて、状態から行動（制御入力）へのマッピングである方策関数を直接学習します。特にモデルベースRLと呼ばれるアプローチでは、学習したシステムモデルを使って計画を行う点でMPCとの類似性が高まります。近年では、MPCの計画能力とRLの学習能力を組み合わせるような研究も盛んに行われています。

これらのアプローチは、計算コスト、モデルへの依存度、要求されるデータ量、性能保証のしやすさなどの点で、それぞれ一長一短があります。対象とする問題やシステムに応じて、NMPCとこれらの代替手法を比較検討したり、組み合わせたりすることが、現実的な解を得る上で重要になってくるのかもしれません。現在過熱しているフィジカルAI（事実上ロボティクスのAI機能強化？）に関しても、再度MPC・実時間最適化・強化学習などの話題が何らかの形で関与していくものと思われます。このあたりも、また別の機会に深掘りしてみたいテーマです。

まとめ

今回は、最近よく聞く「MCP」という言葉から連想した「モデル予測制御（MPC）」について、基本的な考え方から、少し数式を交えた仕組み、簡単なPythonでの実装例、そして非線形MPCとその解法、ソルバーの話、さらには実時間計算の難しさと代替アプローチまで、駆け足でまとめてみました。

MPCのポイントを改めて整理すると、

  • システムのモデルを使って未来を予測する。
  • 評価関数と制約条件の下で、未来の制御計画を最適化する。
  • 最適化して得られた計画のうち、最初のステップだけを実行し、これを繰り返す。

といった感じでしょうか。特に、制約条件を明示的に扱える点は、実用上非常に大きなメリットだと思います。

一方で、良い制御性能を得るためには精度の良いモデルが必要だったり、各ステップでの最適化計算のコスト（計算時間）が課題になったりすることもあります。また、評価関数の重みや予測ホライズンといった調整パラメータをどう設定するか、というのも、なかなか奥が深い問題です。

とはいえ、MPCは非常に強力で柔軟な制御手法であり、化学プラントの制御から自動車の自動運転、ロボット制御、電力システムの運用まで、本当に幅広い分野で活用されています。

今回はこのあたりにしておきます。

興味がある方が続けて読む本

Zennは自分のサーバーじゃないからアフィリエイトリンクはできないのね。。。（まあよいのだけど）

Discussion