はじめに
前回の記事では、Docker環境の構築とMonod式によるバイオリアクターODEモデルを実装し、一定フィードではグルコースが枯渇することを可視化しました。
本記事では、 do-mpc を使ってMPCコントローラーを実装します。グルコースを目標値に維持する制御から始め、DOも同時に制御する多変数MPCへと段階的に発展させます。チューニングの試行錯誤も含めて記録しているので、同じように詰まった方の参考になれば幸いです。
本記事でやること
① do-mpcの4つのコンポーネントを理解する
② グルコース単変数制御を実装する
③ フィード上限が制御性能に与える影響を確認する
④ DO制御を追加して多変数MPCにする
⑤ 制御なし vs MPC の比較グラフを作成する
MPCの仕組みをおさらい
予測ホライズンとは
MPCの核心は「未来を予測して今の操作を決める」ことです。
現在（t=0）
↓
MPCが「今から10時間先まで」を予測
↓
「グルコースが目標値から外れないように
フィード流量をどう動かすか」を最適化
↓
最初の0.5時間分の操作量だけを実行
↓
0.5時間後に再計算（繰り返し）
予測する時間幅を 予測ホライズン（n_horizon） と呼びます。今回は n_horizon=20、t_step=0.5時間 なので、常に 10時間先までを見ながら制御します。
do-mpcの4つのコンポーネント
do-mpcは以下の4つのオブジェクトで構成されています。
┌─────────────────────────────────────┐
│ Model（モデル） │
│ ODEの定義・状態変数・操作変数の登録 │
└──────────────────┬──────────────────┘
│
┌─────────┴──────────┐
▼ ▼
┌─────────────────┐ ┌─────────────────────┐
│ MPC │ │ Simulator │
│ コントローラー │ │ 「操作したら │
│ 「最適な操作量 │ │ 状態はどう変わるか」│
│ を計算する」 │ │ を計算する │
└────────┬────────┘ └──────────┬──────────┘
│ 操作量 u │ 次の状態 y
└──────────┬───────────┘
▼
┌─────────────────────┐
│ Estimator │
│ 「現在の状態をMPCに │
│ 伝える」 │
└─────────────────────┘
| コンポーネント | 役割 | 実機での対応 |
|---|---|---|
| Model | ODEの定義 | プロセスの数式モデル |
| MPC | 最適操作量の計算 | コントローラー本体 |
| Simulator | 操作結果の計算 | 培養槽そのもの（PoCでは仮想） |
| Estimator | 現在状態の取得 | センサー計測値 |
ステップ1：グルコース単変数制御
do-mpcモデルの定義
まず 02_mpc_draft.ipynb を作成し、以下を実装します。
import do_mpc
import casadi as ca
import numpy as np
import sys
sys.path.append('/workspace/models')
# ─────────────────────────────
# do-mpcモデル定義
# ─────────────────────────────
model = do_mpc.model.Model('continuous')
# 状態変数の登録
X = model.set_variable('_x', 'X') # 細胞密度 [g/L]
S = model.set_variable('_x', 'S') # グルコース [g/L]
DO = model.set_variable('_x', 'DO') # 溶存酸素 [%]
# 操作変数の登録
F = model.set_variable('_u', 'F') # フィード流量 [L/h]
# パラメータ
mu_max = 0.4; Ks = 0.1; Y_xs = 0.5
KLa = 100.0; DO_sat = 100.0; OUR = 0.05
V = 2.0; S_in = 50.0
# Monod式
mu = mu_max * S / (Ks + S)
D = F / V
# ODEの登録
model.set_rhs('X', mu * X - D * X)
model.set_rhs('S', -mu / Y_xs * X + D * (S_in - S))
model.set_rhs('DO', KLa * (DO_sat - DO) - OUR * X)
model.setup()
print("✅ モデル定義完了")
前回の bioreactor_ode() と同じ数式ですが、model.set_variable() で登録することでdo-mpcが最適化計算に使えるようになります。
前回の ode_modelはscipyで「過去から未来を計算するシミュレーター」として使い、今回のmpc_modelはdo-mpcで「未来を予測して最適な操作を決めるコントローラ」として使うもので、同じODE・数式でも目的が異なります。
MPCコントローラーの設定
# ─────────────────────────────
# MPCコントローラー設定
# ─────────────────────────────
mpc = do_mpc.controller.MPC(model)
mpc.settings.n_horizon = 20 # 予測ホライズン（20ステップ先を予測）
mpc.settings.t_step = 0.5 # 1ステップ = 0.5時間
mpc.settings.n_robust = 0
mpc.settings.store_full_solution = True
# 目的関数：グルコースを1.0 g/Lに近づける
S_target = 1.0
lterm = (S - S_target)**2 # ステージコスト（毎ステップの誤差）
mterm = (S - S_target)**2 # 終端コスト（最後の誤差）
mpc.set_objective(lterm=lterm, mterm=mterm)
# 操作量変化のペナルティ（急激な変化を抑制）
mpc.set_rterm(F=0.1)
# 制約条件
mpc.bounds['lower', '_u', 'F'] = 0.0 # フィード流量の下限
mpc.bounds['upper', '_u', 'F'] = 0.5 # フィード流量の上限（L/h）
mpc.bounds['lower', '_x', 'X'] = 0.0
mpc.bounds['lower', '_x', 'S'] = 0.0
mpc.bounds['lower', '_x', 'DO'] = 0.0
mpc.setup()
print("✅ MPCコントローラー設定完了")
目的関数の意味：
(S - S_target)**2 はグルコース濃度と目標値のズレの二乗です。MPCはこの値を最小にしようとします。ズレが大きいほど値が大きくなるため、「できるだけ目標値に近づける」という意図を数式で表現しています。
ペナルティ set_rterm の意味：
フィード流量を急激に変えすぎないためのペナルティです。実機のポンプへの負担軽減と、GMP観点での急激な操作変化の抑制対応を想定しています。
制約条件 mpc.bounds の意味：
現実的な制約（例えばフィードポンプの能力など）を数式に組み込み、MPCはこの範囲内でしか操作できません。細胞密度などは、物理的にマイナスにはなりません。
シミュレーター・エスティメーターの設定
# シミュレーター
simulator = do_mpc.simulator.Simulator(model)
simulator.set_param(t_step=0.5)
simulator.setup()
# エスティメーター（StateFeedback：状態変数をそのまま使う最もシンプルな設定）
estimator = do_mpc.estimator.StateFeedback(model)
シミュレーターについて：
MPCが決めた操作量を受け取って培養槽の中の変化を計算する役割をもちます。実機の役割をもち、PoCでは実機の代わりにシミュレーターを利用します。
エスティメーターについて：
StateFeedback は全状態変数をそのままMPCに渡します。実機では全変数をリアルタイム計測できないため、カルマンフィルター等を使うのが本番対応ですが、PoCではシンプルなStateFeedbackで十分です。実機でいうと、センサーに相当します。
3つの役割の関係性
┌─────────────────────────────────┐
│ MPCコントローラー │
│ 「次の操作量Fを計算する」 │
└──────────┬──────────────────────┘
│ F（操作量）
▼
┌─────────────────────────────────┐
│ シミュレーター │
│ 「操作したら状態がこう変わった」 │
└──────────┬──────────────────────┘
│ X, S, DO（新しい状態）
▼
┌─────────────────────────────────┐
│ エスティメーター │
│ 「現在の状態をMPCに伝える」 │
└──────────┬──────────────────────┘
│ 現在の状態
└──→ MPCコントローラーへ
シミュレーション実行
# 初期条件
x0 = np.array([[0.1], [5.0], [100.0]]) # [X, S, DO]
mpc.x0 = x0; simulator.x0 = x0; estimator.x0 = x0
mpc.set_initial_guess()
# メインループ（144ステップ = 72時間）
n_steps = 144
X_log = [0.1]; S_log = [5.0]; DO_log = [100.0]
F_log = [0.0]; t_log = [0.0]
for i in range(n_steps):
u0 = mpc.make_step(x0) # MPCが操作量を計算
y_next = simulator.make_step(u0) # シミュレーターが状態を更新
x0 = estimator.make_step(y_next) # エスティメーターが状態を取得
X_log.append(float(x0.flatten()[0]))
S_log.append(float(x0.flatten()[1]))
DO_log.append(float(x0.flatten()[2]))
F_log.append(float(u0.flatten()[0]))
t_log.append((i + 1) * 0.5)
print(f"最終細胞密度 : {X_log[-1]:.2f} g/L")
print(f"最終グルコース: {S_log[-1]:.2f} g/L")
1ループで何が起きているか：
① mpc.make_step(x0)
→ 「今の状態x0を見て10時間先まで予測し
グルコースを1.0 g/Lに近づける最適なFを計算」
② simulator.make_step(u0)
→ 「そのFで0.5時間後の培養槽はどうなるか計算」
③ estimator.make_step(y_next)
→ 「新しい状態をMPCに渡せる形に整える」
→ これを144回繰り返す（=72時間）
結果：F_max=0.5では制御不足
最終細胞密度 : 24.92 g/L
最終グルコース: 0.17 g/L ← 目標1.0 g/Lに届いていない
グラフを見ると、フィード流量Fが 上限の0.5 L/hに張り付いている ことがわかります。つまり「もっとフィードしたいが、上限に阻まれている」状態です。
12〜13時間付近でフィードが一時的に下がっているのは、MPCの「先読み」が働いている証拠です。グルコースが急減する前に「このままフィードを増やし続けると過剰になる」と予測し、一時的に絞っています。PIDのような反応型制御では実現できない動作です。
ステップ2：フィード上限を変えて再チューニング
F_max=0.5が不足だった理由
細胞が増えるほどグルコース消費量が増える
↓
グルコースを1.0 g/Lに維持するには
より多くのフィードが必要
↓
F_max=0.5 L/h では物理的に追いつかない
フィード上限を 0.5 → 1.0 L/h に変更します。
mpc.bounds['upper', '_u', 'F'] = 1.0 # 0.5 → 1.0に変更
結果：F_max=1.0でグルコース制御成功
最終細胞密度 : 24.50 g/L
最終グルコース: 1.00 g/L ✅ 目標値に制御できた
グラフのフィード流量を見ると、上限1.0 L/hに余裕があるため、きれいなS字カーブを描きながら安定して制御できていることがわかります。
ここで得られた実機への示唆
フィードポンプの最大流量が制御性能のボトルネックになる可能性がある。
実機設計時にポンプ能力を確認する必要がある。
ステップ3：DO制御を追加（多変数MPC）
なぜDOを追加するか
前のステップではDOが常に100%のままでした。実際の培養では細胞が酸素を消費するためDOは低下します。DOが低すぎると細胞の代謝に影響するため、目標値（今回は40%）付近に維持したい変数です。
モデルの変更点
DO制御のために2つの変更を加えます。
# ① OURを現実的な値に変更（細胞の酸素消費を大きくする）
OUR = 5.0 # 0.05 → 5.0に変更
# ② 通気量Qを操作変数として追加
Q = model.set_variable('_u', 'Q') # 通気量（DO制御用）
# KLaを通気量Qで動的に変化させる
KLa_base = 10.0
KLa = KLa_base * Q # 通気量が増えるほどKLaが上がる
変更の意味：
| 変更点 | 変更前 | 変更後 | 意味 |
|---|---|---|---|
| OUR | 0.05 | 5.0 | 細胞が酸素をしっかり消費するようになる |
| KLa | 固定値 | KLa_base × Q | 通気量でKLaが変わる = DOを制御できる |
| 操作変数 | F のみ | F + Q | グルコースとDOを独立して制御できる |
目的関数の変更：2変数を同時に制御
S_target = 1.0 # グルコース目標値 [g/L]
DO_target = 40.0 # DO目標値 [%]
# グルコースとDOの両方を目標値に近づける
lterm = (S - S_target)**2 + 0.1 * (DO - DO_target)**2
mterm = (S - S_target)**2 + 0.1 * (DO - DO_target)**2
重み 0.1 の意味：
グルコースとDOのスケールが異なるため、そのまま足すと一方が支配的になります。0.1 という係数でDOの重みを調整しています。
重みを大きくする → DOをより優先して制御する
重みを小さくする → グルコースをより優先する
これがMPCチューニングの核心です。
制約条件の追加
# フィード流量
mpc.bounds['lower', '_u', 'F'] = 0.0
mpc.bounds['upper', '_u', 'F'] = 1.0
# 通気量
mpc.bounds['lower', '_u', 'Q'] = 0.1 # 下限0.1（最低限の通気を維持）
mpc.bounds['upper', '_u', 'Q'] = 2.0 # 上限2.0
# DOの上限（物理的に100%を超えない）
mpc.bounds['upper', '_x', 'DO'] = 100.0
通気量の下限を0ではなく0.1にしているのは、現実の培養槽では最低限の通気は常に維持するためです。
結果：DOが動くようになった
最終細胞密度 : 24.14 g/L
最終グルコース: 1.72 g/L （目標付近）
最終DO : 54.37 % （目標40%に近づいた）
OURを大きくしたことで細胞が酸素を消費し、DOが100%→54%まで低下しました。目標の40%には届いていませんが、これはパラメータの調整（後述）で改善できます。
ステップ4：3パターン比較グラフ
比較グラフの作成
3つのシミュレーション結果を1つのグラフに重ねて比較します。
import matplotlib.pyplot as plt
fig, axes = plt.subplots(4, 1, figsize=(14, 12), sharex=True)
# 細胞密度
axes[0].plot(t_log, X_log, color='gray', linewidth=2,
linestyle='--', label='No Control')
axes[0].plot(t_log2, X_log2, color='steelblue', linewidth=2,
label='MPC F_max=0.5')
axes[0].plot(t_log5, X_log5, color='darkblue', linewidth=2,
label='MPC F+DO')
axes[0].set_ylabel('Cell Density X [g/L]')
axes[0].legend(); axes[0].grid(True)
# グルコース
axes[1].plot(t_log, S_log, color='gray', linewidth=2,
linestyle='--', label='No Control')
axes[1].plot(t_log2, S_log2, color='darkorange', linewidth=2,
label='MPC F_max=0.5')
axes[1].plot(t_log5, S_log5, color='red', linewidth=2,
label='MPC F+DO')
axes[1].axhline(y=1.0, linestyle=':', color='black', label='Target 1.0 g/L')
axes[1].set_ylabel('Glucose S [g/L]')
axes[1].legend(); axes[1].grid(True)
# 溶存酸素
axes[2].plot(t_log, DO_log, color='gray', linewidth=2,
linestyle='--', label='No Control')
axes[2].plot(t_log2, DO_log2, color='lightgreen', linewidth=2,
label='MPC F_max=0.5')
axes[2].plot(t_log5, DO_log5, color='green', linewidth=2,
label='MPC F+DO')
axes[2].axhline(y=40.0, linestyle=':', color='black', label='Target 40%')
axes[2].set_ylabel('Dissolved Oxygen DO [%]')
axes[2].legend(); axes[2].grid(True)
# フィード流量
axes[3].step(t_log, F_log, color='gray', linewidth=2,
linestyle='--', label='No Control')
axes[3].step(t_log2, F_log2, color='steelblue', linewidth=2,
label='MPC F_max=0.5')
axes[3].step(t_log5, F_log5, color='darkblue', linewidth=2,
label='MPC F+DO')
axes[3].set_ylabel('Feed Flow F [L/h]')
axes[3].set_xlabel('Time [h]')
axes[3].legend(); axes[3].grid(True)
plt.suptitle('Bioreactor Simulation - Comparison', fontsize=14)
plt.tight_layout()
plt.savefig('/workspace/results/simulation_comparison.png', dpi=150)
plt.show()
比較結果
| 条件 | 最終細胞密度 | グルコース（72h） | DO（72h） |
|---|---|---|---|
| 制御なし（F=0.05固定） | 21.29 g/L | 0.01 g/L（枯渇）❌ | 99.99% |
| MPC F_max=0.5 L/h | 24.92 g/L | 0.17 g/L（不足）❌ | 99.99% |
| MPC F+DO制御 | 24.14 g/L | 1.72 g/L（目標付近）✅ | 54.37% ✅ |
実機への示唆
このシミュレーションから以下の知見が得られました。
①フィードポンプの最大流量が制御性能を決める
F_max=0.5 L/h → グルコース枯渇
F_max=1.0 L/h → グルコース制御成功
→ 実機のポンプ最大流量を事前に把握することが
MPC設計の前提条件になる
②一定フィードでは必ずグルコースが枯渇する
制御なし → 12時間でグルコース枯渇
MPC → 72時間を通じてグルコースを維持
→ 細胞密度が上がるほど消費量が増えるため
フィードも動的に増やす必要がある
③MPCの「先読み」はPIDにはできない動作
12〜13時間付近でフィードが一時的に減少
→「過剰になる前に絞る」という先読み制御
→ PIDでは反応してからしか調整できない
④グルコースとDOはトレードオフの関係
グルコース制御を優先 → DOが目標に届かない
DO制御の重みを増やす → グルコースがやや不安定
→ 実機データでパラメータ（OUR・KLa）を同定し
適切な重みを設定することが次フェーズの課題
チューニング過程のまとめ
今回の試行錯誤を整理します。
| 試行 | 設定 | 結果 | 判定 |
|---|---|---|---|
| 1 | グルコースのみ制御・F_max=0.5 | グルコース枯渇（0.17 g/L） | ❌ |
| 2 | グルコースのみ制御・F_max=1.0 | グルコース目標値1.0達成 | ✅ |
| 3 | DO追加・OUR=2.0・重み0.01 | DO=81.93%・Q下限張り付き | △ |
| 4 | DO追加・OUR=2.0・重み0.1 | DO=81.74%・改善なし | △ |
| 5 | DO追加・OUR=5.0・重み0.1 | DO=54.37%まで低下 | ✅ |
チューニングで学んだこと：
-
n_horizon（予測ホライズン）を長くするほど計算時間が増える -
set_rtermのペナルティが大きすぎると操作量が動かなくなる - 多変数制御ではスケールの違いに注意して重みを設定する
- 目標値に届かない場合は制約条件が原因のことが多い
残った課題
今回のPoCで見えた技術的課題を整理します。
技術的課題
| 課題 | 原因 | 対応方針 |
|---|---|---|
| DOが目標40%に届かない | OUR・KLaの設定が実機と乖離 | 実機データでパラメータ同定 |
| 通気量Qが制御されていない | モデル構造の問題 | 次フェーズで見直し |
| パラメータが文献値 | 実機とのズレがある | scipy.optimizeでフィッティング |
GMP観点での課題
製薬・バイオ業界でMPCを本番実装する際には、以下も検討が必要です。
- MPCモデルのCSV対象範囲：モデル自体がGAMP5のコンピュータ化システムバリデーション対象になるか
- パラメータ変更の変更管理：チューニングのたびに変更管理が必要か
- 制御ログの電子記録：21 CFR Part 11への対応
まとめ
本記事では以下を実装・確認しました。
- ✅ do-mpcの4コンポーネント（Model / MPC / Simulator / Estimator）を理解
- ✅ グルコース単変数MPCを実装（フィード上限の影響を確認）
- ✅ DO制御を追加して多変数MPCに拡張
- ✅ 3パターン比較グラフで制御効果を可視化
- ✅ 実機への示唆を言語化
次回予告
次の記事 では、シミュレーション結果を OPC UA Pub/Sub over MQTT でデータ送受信し、Medallion 3層アーキテクチャ（Bronze/Silver/Gold） のデータ基盤に蓄積、Grafana でダッシュボード化します。製薬・バイオ業界のGMP観点（ALCOA+・トレーサビリティ）も含めて解説します。
参考
- do-mpc公式ドキュメント
- CasADi公式ドキュメント
- Rawlings, J.B., Mayne, D.Q., Diehl, M. (2017). Model Predictive Control: Theory, Computation, and Design
本記事はPoCとして構築したものであり、実機への適用にはGMP対応・バリデーションが別途必要となります。