制御アルゴリズムを MATLABで書く～C/GMRES法による非線形モデル予測制御～

制御アルゴリズムを MATLAB

で書く

～ C/GMRES 法による非線形モデル予測制御 ～

MathWorks Japan

アプリケーションエンジニアリング部 制御 赤阪 大介

アジェンダ

 はじめに

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

 MATLAB で非線形モデル予測制御 (NMPC) を書く

 まとめ

近年、高度アルゴリズムの適用への動きが飛躍的に進んでいる

製品・システムとしての価値の追求が進む中で 求められる機能・性能はより 高度 かつ 複雑 に

複雑な系に対して、より高性能な制御を実現 アルゴリズム

フィードバック制御 技術

センサ情報処理 技術

データ解析 技術

モデル予測制御 (MPC)

お伝えしたいこと

日本の産業界で MPC 適用へ向けた取り組みを より活発化させてゆきませんか？

制御アルゴリズムの統合開発環境として MATLAB を是非ご活用ください

提起

提案

アジェンダ

 はじめに

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

 MATLAB で非線形モデル予測制御 (NMPC) を書く

 まとめ

MPC の鍵は 「予測」 と 「最適化」

プラント

最適化器

予測モデル 目標

指令

MPC

制御

入力 制御

出力

r(t) u(t) y(t)

測定

プラントの数理モデル を利用して未来の挙動を 予測 し

各時刻で 最適化問題 を解きながら制御入力を決定する方法

目的：

ラップタイムの最小化 制約：

•

•

•

•

MPC は最適化ベースの制御

最適制御問題 と捉えて解く

予期せぬことが 起こる可能性も

MPC は走りながら、計測しながら、予測して考えながら 常に 経路の再計画 や 車両の動き方の決定 を行う

先を予測し

車両運動制御 ラップタイムを

短縮したい

プラント

最適化器

予測モデル 目標

指令

MPC

制御

入力 制御

出力

r(t) u(t) y(t)

測定

予測 & 最適化により制御入力を決定

現時刻 未来 上限

下限

次の時間ステップの 入力として使用

現時刻 未来 上限

下限 目標値

例: 誤差面積が最小に なるように入力を決定

予測区間

MPC がなぜ求められているか？

多入出力系のような 複雑な対象 へ適用しやすい 複雑な系に対して、より高度な制御を実現する手段

制御系の性能を最大限引き出すための基盤技術

1. 汎用性 2. 制約遵守 3. 協調上手 4. 最適性

汎用的な手法ゆえ、幅広い応用先 が考えられる 物理・性能・安全などの 制約 を守りながら制御

最適制御 に実用的な フィードバック制御 の補正効果

MPC の適用先は ?

ナノ秒

マイクロ秒 ミリ秒

秒 分 時 日 週

コンピュータ制御 パワエレ・電力系 トラクション制御 ビル

石油精製

人員配置計画

電車スケジューリング 生産計画

計算機 (HW) とアルゴリズム (SW) の進化に伴い適用範囲は拡大

例) 自動車分野で想定される適用先

 エンジン・パワートレイン制御

 車両運動制御

 エネルギーマネジメント

 車群制御 etc

サンプリング時間

産業分野の広がり

エコ・自動運転支援・V2X

MPC の可能性 ～ 車両のダブルレーンチェンジ

MPC

PID

0 m 20 m 40 m 60 m 80 m 100 m 120 m

目標軌道 操舵角制限 ±10°

目標車速 10m/s 滑りやすい路面 (雪)

6m

プラント

最適化器

予測モデル 目標

指令

MPC

制御

入力 制御

出力

r(t) u(t) y(t)

測定

モデル予測制御は最適制御問題に帰着する

最小化:

制約条件

制御入力のパターン 探索:

目的関数

u

現時刻

𝑢

𝑢

𝑢

予測区間

t

MPC で解くべき問題

① 制御目標 の定式化

② 数値計算アルゴリズム

解くべき最適化問題を一般的に定式化すると …

目的関数

制約条件 プラントの予測モデル

𝑥(𝑡) = 𝑓(𝑥 𝑡 , 𝑢 𝑡 )

拘束条件

𝐶 𝑥(𝑡), 𝑢(𝑡) ≥ 0

𝐽 = 𝜑 𝑥 𝑡

+ 𝑇 +

𝑡₀

𝑡₀+𝑇

𝐿 𝑥 𝜏 , 𝑢 𝜏 𝑑𝜏

𝑥 𝑡

= 𝑥

系のエネルギー消費やコスト

初期条件

目的関数を最小化する制御入力のパターン

範囲は現時刻

𝑡

𝑡

+ 𝑇

探索

終端コスト エネルギーコスト

𝑢

𝑡

𝑡

+ 𝑇

常微分方程式 現時刻での状態 物理・性能・安全

例題: DC モーターの角度制御

目的関数

制約条件

^𝑥

𝑥

𝑥

=

0 1 0

0 −𝑏/𝐽 𝐾

/𝐽 0 −𝐾

/𝐿 −𝑅/𝐿

𝑥

𝑥

𝑥

+

0 0

1/𝐿 𝑢 𝑥

𝑥

𝑥

−5 ≤ 𝑢 ≤ 5

−𝜋/2 ≤ 𝑥

≤ 𝜋/2 𝐽 =

𝑡₀

𝑡₀+𝑇

{𝑞 𝑥

(𝜏) − 𝜋/4

+ 𝑟 𝑢(𝜏)

}𝑑𝜏 𝑢

(範囲制約 ±5V)

角度追従偏差 制御入力の大きさ

プラントの予測モデル

運動・回路方程式

制御性能の評価関数

追従性と入力大きさを評価

拘束条件

角度・電圧の範囲制約

プラント

目標角度 45° (

𝜋/4

MPC の最適化問題を解くアプローチは様々

最適制御問題 動的計画法

間接法 直接法

Hamilton-Jacobi-

Bellman 方程式 を解く 変分法から導出される

多点境界値問題 を解く 非線形計画問題 に変換して解く

オフラインでの計算方法

オンラインでの計算方法

間接法の手法

 勾配法 例: ニュートン法、最急降下法、共役勾配法 …

 連続変形法 例: C/GMRES 法 …

なぜ C/GMRES 法か？

MPC の課題である 非線形プラントの扱い や 計算量の低減 を克服

Continuation/GMRES 法は 非線形 MPC の 高速解法

最適制御の時間依存を利用し 反復計算を回避

1. 連続変形法 2. GMRES 法 3. 差分近似

速さのポイント

各サンプリング時刻で連立 1 次方程式を一回解くだけ

連立 1 次方程式の 効率的な解法 を利用

ヤコビ行列とベクトルの積を近似し 演算量を低減

変分法から導出される 「2点境界値問題」 を解く…課題は？

𝑥(𝑡) = 𝑓(𝑥(𝑡), 𝑢(𝑡)) 𝑥

𝑡 = 𝑥 𝑡

𝐻

𝑥(𝑡), 𝜆 𝑡 , 𝑢 𝑡 , 𝜇 𝑡 = 0 𝜆 𝑡 = −𝐻

𝑥(𝑡), 𝜆 𝑡 , 𝑢 𝑡 , 𝜇 𝑡

𝐶 𝑥(𝑡), 𝑢(𝑡) = 0

𝜆 𝑡 + 𝑇 = 𝜑

(𝑥(𝑡 + 𝑇))

𝐻 𝑥, 𝜆, 𝑢, 𝜇 = 𝐿(𝑥, 𝑢) + 𝜆

𝑓(𝑥, 𝑢) + 𝜇

𝐶(𝑥, 𝑢)

ハミルトニアン：

初期条件

終端条件

𝑢 𝑡

𝜇 𝑡

課題 2 必要な数式の導出は?

解決策 数式処理技術を利用した 数式の導出が有効

課題 1 境界値問題をどう解くか?

解決策 予測区間にわたり離散化し 非線形連立方程式に帰着

予測区間で 「離散化」 し差分方程式へ

プラントの予測モデル

ハミルトニアンの偏微分 初期条件 (測定された状態)

随伴方程式 終端条件 (終端コストの偏微分)

制約条件

𝑥

: 𝑖

𝜆

:

𝑢

:

𝜇

:

𝑥

𝑡 = 𝑥

𝑡 + 𝑓 𝑥

𝑡 , 𝑢

𝑡 Δ𝜏 𝑥

𝑡 = 𝑥 𝑡

𝐻

𝑥

𝑡 , 𝜆

𝑡 , 𝑢

𝑡 , 𝜇

𝑡 = 0

𝜆

𝑡 = 𝜆

𝑡 + 𝐻

𝑥

𝑡 , 𝜆

𝑡 , 𝑢

𝑡 , 𝜇

𝑡 Δ𝑡

𝐶 𝑥

𝑡 , 𝑢

𝑡 = 0 𝜆

𝑡 = 𝜑

(𝑥

(𝑡))

2 点境界値問題

𝑢 𝑡

𝜇 𝑡

𝐻 𝑥, 𝜆, 𝑢, 𝜇 = 𝐿(𝑥, 𝑢) + 𝜆

𝑓(𝑥, 𝑢) + 𝜇

𝐶(𝑥, 𝑢)

評価コスト 予測モデル 制約条件

MPC は 「非線形方程式」 を解く問題に帰着する

制御入力

𝒖

𝝁

𝑈 𝑡 = 𝑢

𝑡 , 𝜇

𝑡 , ⋯ , 𝑢

𝑡 , 𝜇

𝑡

𝐹 𝑈 𝑡 , 𝑥 𝑡 , 𝑡 ≔

𝐻

𝑥

𝑡 , 𝜆

𝑡 , 𝑢

𝑡 , 𝜇

𝑡 𝐶 𝑥

𝑡 , 𝑢

𝑡

⋮

𝐻

𝑥

𝑡 , 𝜆

𝑡 , 𝑢

𝑡 , 𝜇

𝑡 𝐶 𝑥

𝑡 , 𝑢

𝑡

= 0

各時刻 𝒕 で測定された状態 𝒙(𝒕) に対して

方程式 𝑭(𝑼, 𝒙, 𝒕) = 𝟎 を解き 𝑼(𝒕) を計算

時間パラメータ t に沿って時刻と共に解 U(t) を追跡できる

C/GMRES 解法: 連続変形法の実時間利用による反復計算の回避

𝐹 𝑈 𝑡 , 𝑥 𝑡 , 𝑡 = 0

基本アルゴリズム

課題 ニュートン法のような反復解法は非効率

𝐹 𝑈, 𝑥, 𝑡 = −𝑎𝐹(𝑈, 𝑥, 𝑡) ^{(𝑎 > 0)}

解決策

𝐹 = 0

𝑈(𝑡)

𝐹 = 𝐹 _𝑈 𝑈 + 𝐹 _𝑥 𝑥 + 𝐹 _𝑡 = −𝑎𝐹

_𝑈

𝐹 𝑈 0 , 𝑥 0 , 0 = 0

𝑈(0)

1.

2.

𝐹

𝑡

𝐹 = 0

𝐹(𝑡) = 𝑒

C/GMRES 解法: GMRES法 と 前進差分近似 による数値計算の効率化

前進差分近似 し、F の評価回数を減らす

解決策

1.

2.

課題 ヤコビ行列

𝐹

𝐹

𝐹

𝐹

𝐹 𝑈 + ℎ𝑈, 𝑥 + ℎ 𝑥, 𝑡 + ℎ − 𝐹 𝑈, 𝑥, ℎ

𝐹 𝑈, 𝑥, 𝑡 ≈ ℎ = −𝑎𝐹(𝑈, 𝑥, 𝑡)

𝐹 _𝑈 𝑈 + 𝐹 _𝑥 𝑥 + 𝐹 _𝑡 = −𝑎𝐹

𝐴 𝑈 = 𝑏

𝑈

要点まとめ

1. MPC は 実時間最適制御

–

–

2. C/GMRES 法 は非線形 MPC の高速解法

–

–

では、アルゴリズム構築や制御系の検証をどう行いますか？

皆様のご業務で MPC を検討する価値がありそうでしょうか？

アジェンダ

 はじめに

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

 MATLAB で非線形モデル予測制御 (NMPC) を書く

 まとめ

なぜ MATLAB で書くのか？

数値シミュレーション技術 数式処理技術

プログラミング環境

行列計算・数学関数

自動コード生成技術 データ分析・評価技術

制御アルゴ開発

アルゴリズム開発

C/GMRES アルゴリズム開発ワークフロー

システム設計・検証

M 関数生成

数式処理 MATLAB 関数作成

MEX 生成 C 生成

C/GMRES 法 を M 言語でプログラム

シミュレーション検証

閉ループ系シミュレーション による制御アルゴの評価 制御アルゴに必要

な数式の導出

高速化

MEX C 関数

（dll）

実装

制御 プラント

シミュレーション検証

 性能指標

 相関/感度

 ロバスト性

 トレードオフ etc 機能・性能評価

MATLAB Function, S-Function ブロック

例題: 車両のダブルレーンチェンジ

𝑌

𝑋

𝜑 𝜔

𝛿 𝛼

𝛼

𝑉

𝑉

𝐹

(𝛼

)

𝐹

(𝛼

)

𝐹

𝑙

𝑙

4 輪車両の等価２輪モデル

𝐹_𝑟𝑦(𝛼_𝑟) 𝛼_𝑟

𝐹_𝑓𝑦(𝛼_𝑓) 𝛼_𝑓

プラントモデル

制約 (操舵角)

𝑋 = 𝑉

cos 𝜑 − 𝑉

sin(𝜑) 𝑌 = 𝑉

sin 𝜑 + 𝑉

cos(𝜑)

𝜑 = 𝜔

𝑥

𝑉 = (2𝐹

− 2𝐹

𝑠𝑖𝑛𝛿 + 𝑚𝑉

𝜔)/𝑚

𝑦

𝑉 = (2𝐹

+ 2𝐹

𝑐𝑜𝑠𝛿 − 𝑚𝑉

𝜔)/𝑚

𝜔 = (2𝐹

𝑙

𝑐𝑜𝑠𝛿 − 2𝐹

𝑙

)/𝐼

制動力

入力2:

操舵角

𝛿

≤ 𝛿 ≤ 𝛿

横力

スリップ角

問題設定

0 m 20 m 40 m 60 m 80 m 100 m 120 m

目標: 進行方向に一定速度で目標軌道 (位置・姿勢) に追従

𝑌 0 𝑋

•

•

目標軌道

滑りやすい路面 (雪)

6m

例題の流れ

1. 数式処理技術の活用

2. C/GMRES 法アルゴリズムの記述

3. 自動コード生成技術の活用

4. ブロック線図表現のシミュレーション

1. 数式処理技術の活用

数式が非線形になる.. 試行を重ねたい.. 自動化したい場合は特に有効

必要な数式の解析的な操作・導出には数式処理が有効

𝐹(𝑈, 𝑥, 𝑡) = 0

𝐻 𝑥, 𝜆, 𝑢, 𝜇 = 𝐿 𝑥, 𝑢 + 𝜆

𝑓 𝑥, 𝑢 + 𝜇

𝐶(𝑥, 𝑢)

𝜕𝐻

𝜕𝑢

ハミルトニアン

𝐻

𝜑

𝐻

, 𝐻

, 𝜑

𝐿 𝑥, 𝑢 = 𝑥 − 𝑥_𝑓 ^𝑇𝑄(𝑥 − 𝑥_𝑓) + 𝑢^𝑇𝑅𝑢 − ∑𝑘_𝑖 log(𝑢_𝑖 − 𝑢_𝑚𝑎𝑥)(𝑢_𝑚𝑖𝑛 − 𝑢_𝑖)

𝑓 𝑥, 𝑢 =

𝑥₄cos 𝑥₃ − 𝑥₅sin(𝑥₃) 𝑥₄sin 𝑥₃ + 𝑥₅cos(𝑥₃)

𝑥₆

(2𝑢₁ − 2𝐹_𝑓𝑦(𝑥)𝑠𝑖𝑛𝑢₂ + 𝑚𝑥₅𝑥₆)/𝑚 (2𝐹_𝑟𝑦(𝑥) + 2𝐹_𝑓𝑦(𝑥)𝑐𝑜𝑠𝑢₂ − 𝑚𝑥₄𝑥₆)/𝑚

(2𝐹 (𝑥)𝑙 𝑐𝑜𝑠𝑢 − 2𝐹 (𝑥)𝑙 )/𝐼 例)

？

𝜕𝐻

𝜕𝑥

Symbolic Math Toolbox ～ MATLAB 数式処理ツール

 微積分 (勾配、ヤコビ行列、ヘシアン)

 数式の代入・単純化

 代数計算、微分方程式の求解

Symbolic Math Toolbox MuPAD Notebook

(リッチテキスト形式の専用エディタ)

 変換、特殊関数、線形代数

 可変精度計算

 コード生成 (MATLAB, Simscape言語)

𝑥 𝑡

𝛥𝑥

𝛥𝑡

:

𝑡

2. C/GMRES 法アルゴリズムの記述

数式処理を通じて 作成した関数を使用 初期化関数

初期値 U(0) 計算 C/GMRES 関数 (作成例)

𝑈

fdgmres (𝐹, 𝑈 𝑡 , 𝑥 𝑡 ,

Δ𝑡

, 𝑡, 𝑝𝑎𝑟𝑎𝑚) 𝑈 𝑡 + 𝛥𝑡 = 𝑈 𝑡 + 𝑈(𝑡)𝛥𝑡

^{𝑈 𝑡}

𝑢 𝑡 + 𝛥𝑡 = 𝑃

(𝑈 𝑡 + 𝛥𝑡 )

F 計算用関数 lib

𝐹(𝑈(𝑡), 𝑥, 𝑡)

プラント非依存関数 (汎用)

U(0) Call

制御開始前に

勾配法などで計算

サンプリング時間

Δ𝑡

パラメータ設定用 ファイル

入力

出力 MATLAB で簡潔なプログラミング

3. 自動コード生成技術の活用

M M M

C

MEX MATLAB から 呼び出せるDLL

MPC プラント

Simulink

モデルサブシステム MATLAB

関数ファイル

制御アルゴリズム

C ソースコード

MATLAB Coder

Simulink Coder/

Embedded Coder

シミュレーション高速化

MATLAB Coder/

Embedded Coder ラピッド/組込み実装

1. プロトタイプ試験目的 2. 組込み実装目的

高速化・実装への近道

4. ブロック線図表現のシミュレーション

制御アルゴリズムを Simulink へ

非線形 MPC モデル

MATLAB Function ブロック

M M M

制御パラメータ 状態

時間 目標値

U

U

プラントモデル

スリップ角計算 微分方程式左辺の計算

(数式処理の結果を利用)

状態 制御入力

積分

要点まとめ

 MATLAB は制御アルゴリズムの統合開発環境

1. 数式処理技術の活用

2. C/GMRES 法アルゴリズムの記述

3. 自動コード生成技術の活用

4. ブロック線図表現のシミュレーション

制御アルゴ開発に MATLAB 環境と道具をぜひ有効活用ください 開発者の思考や試行錯誤を助け、開発をより一層快適にします

アジェンダ

 はじめに

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

 MATLAB で非線形モデル予測制御 (NMPC) を書く

 まとめ

まとめ

皆様の MPC 適用への活動を促進する一助となれば幸いです

日本の産業界で MPC 適用へ向けた取り組みを より活発化させてゆきませんか？

制御アルゴリズムの統合開発環境として MATLAB を是非ご活用ください

提起

提案

MathWorks 製品、ワークフローを評価・導入してゆきませんか？

ご清聴いただき誠にありがとうございました

赤阪 大介 MathWorks Japan

[email protected]

© 2015 The MathWorks, Inc. MATLAB and Simulink are registered trademarks of The MathWorks, Inc. See www.mathworks.com/trademarks for a list of additional trademarks. Other product or brand