はじめに
近年注目を集めている制御技術の一つとして、実時間最適制御のMPC（モデル予測制御）があります。
このMPCですが、「モデル予測」とあるように制御をしていく中で制御対象の数式モデルを用いた有限な未来の応答を制御周期毎に予測しながら各種制約条件を満たす最適な操作量を計算します。
しかし、MPCの適用において、いかに程良く制御対象を模した数式モデルを得られるかがハードルとなるケースがよくあります。
例として、明示的な原理式が存在しない、制御対象のパラメーターが未知など、数式モデルの構築に当たっては悩みが多くありますが、この数式モデルをニューラルネットワークによって代替してしまおうという考え方が古くから提案されてきました。
今回はMPCの予測モデルとしてニューラルネットワークを用いた設計をMATLAB及びSimulinkの環境を使用して試してみたいと思います。
住宅暖房システムのエネルギーマネジメント
題材として、MathWorksから例題として提供されている住宅暖房システムに対するエネルギーマネジメントシステムの内容を私のなりの解釈も交えて解説します。
今回の住宅の暖房システムでは、室内の温度をヒーターを使って温めます。
ただし、室温は外気との間で壁、窓、屋根を通して熱交換を伴うため、コントローラーは外気温の変動に対し室温をロバストに制御しなければなりません。
そのような中で今回の制御の目的は以下で与えられます。
・室内温度を快適な温度範囲である 20～22 [degC] に収まるようにする
・ヒーター使用に伴う電気料金を最小化する
今回の例題では、電気料金の単価が24時間を通して以下の様にプロファイルされるとします。
特に10時から21時に亘っては、価格が上昇します。
以上のような問題設定の下で、MPCは所定の温度制約を守りつつ、オンピーク時における電気料金を最小化するような制御を行わなければなりません。
Neural State-Spaceを用いた住宅暖房環境の状態空間システムの推定
室内環境の熱システムをAIモデルの一種であるNeural State-Spaceを用いて学習していきます。
MATLABのSystem Identification Toolboxでは、MLP（Multi Layer Perceptron）に基づく状態空間システムの学習をサポートする機能 Neural State-Space がR2022bより導入されました。
これを使ってMPCのための予測モデルをまずは推定していきます。
予測モデルのトレーニングおよび検証のためのデータセットはデモ中で用意されています。
その中身としては、計3996点分の実験データ（6時間分）が収録されており、データセットに含まれる各時系列には1つの状態量（室内温度）と2つの入力（ヒーター入力と外気温度）が含まれている形となります。
ここで、ヒーター入力は操作量（制御入力）であり、外気温度は測定された外乱となります。
データセットではトレーニング用に5分間隔（300秒）で収録された72点分のデータが複数と、検証用に24時間分の1つデータがMATファイルに格納されています。
下図は検証用データ（24時間）のトレンドです。
このデータを使って、システムの挙動を示す状態空間システムをNeural State-Spaceによって学習していきます。
% サンプリング時間
Ts = 300;
% 状態数：1、入力：2、離散時間TsのNeural State-Spaceを構築
obj = idNeuralStateSpace(1, NumInputs=2, Ts=Ts);
% 状態関数fについて、2つの隠れ層かつサイズがそれぞれ32のMLPとして推定する
obj.StateNetwork = createMLPNetwork(obj, "state", LayerSizes=[32 32]);
% トレーニングオプションとしてBPのソルバーにAdamを指定
options = nssTrainingOptions("adam");
% 最大エポック数を指定する
options.MaxEpochs = 300;
% ミニバッチの一つあたりのサイズを指定する
options.MiniBatchSize = 999;
% 学習率を指定する
options.LearnRate = 0.002;
% nlssestでトレーニング開始、トレーニングデータの最後の値を検証用に使う
nss = nlssest(U,X,obj,options,"UseLastExperimentForValidation",true);
Neural State-Spaceを設計するためにはまずオブジェクトを作る必要があります。
そのためのAPIがidNeuralStateSpaceというコマンドです。
本コマンドには、学習させたいシステムの形として、状態数や入出力数、サンプル時間、状態関数及び出力関数に対するネットワーク構造、レイヤーサイズ、ソルバー設定等のプロパティを指定することが出来ます。
なお、Neural State-Spaceはシステムの状態がすべて観測可能な場合についてのみ適用することが出来るため、適用できる対象が限定されます。
また出力の数に関しても状態量の数(n)≦出力の数(l)という制約があります。残りのl-n個の出力はn個の状態量によって特徴づけられる出力値となります。
以上でネットワークの構造が決まったらあとは学習用のコマンドnlssestにトレーニング＆検証を含んだデータセットとともに渡すことで学習が始まります。
なお、学習中は精度を表す損失値と検証データに対する結果を表示してくれるため、学習の進捗具合を確認することができます。
デモでは学習済みのネットワークが既にMATファイルとして提供されているため、ネットワークのトレーニング過程をスキップすることが可能になっています。
学習が終了したNeural State-SpaceオブジェクトはSimulinkのブロック上から参照させてシミュレーションすることもできます。
上図はNeural State-Space Modelブロックを使って学習済みのオブジェクトをMATLABのワークスペースより参照させて実行した結果になりますが、まずまずの予測精度を有していることが分かります。
なお、検証用のデータはExcelに保存し、From Spreadsheetブロックを使ってモデル上に流し込んでいます。
Neural State-Space からMPCのための予測モデルを生成
前節で設計したNeural State-SpaceのオブジェクトからMPCで使用するための予測モデルの関数を生成していきます。
専用のコマンドgenerateMATLABFunctionがあるので、これを使って関数のMファイルを自動生成できます。
コマンドの引数に生成する状態関数のファイル名"stateFcn"を指定します。
なお、今回は出力＝状態としているため、出力関数に相当するネットワークモデルは必要ありません。
generateMATLABFunction(nss,"stateFcn"); % stateFcnという名前で状態関数のファイルを生成
実行すると以下のような状態関数及びそのヤコビアンに関するコードが自動で生成されます。
function x1 = stateFcn(x,u)
%% auto-generated state function of neural state space system
%# codegen
persistent StateNetwork
MATname = 'stateFcnData';
if isempty(StateNetwork)
StateNetwork = coder.load(MATname);
end
out = [x;u];
% hidden layer #1
out = StateNetwork.fc1.Weights*out + StateNetwork.fc1.Bias;
out = deep.internal.coder.tanh(out);
% hidden layer #2
out = StateNetwork.fc2.Weights*out + StateNetwork.fc2.Bias;
out = deep.internal.coder.tanh(out);
% output layer
dx = StateNetwork.output.Weights*out + StateNetwork.output.Bias;
x1 = x + dx;
function [A, B] = stateFcnJacobian(x,u)
%% auto-generated state Jacobian function of neural state space system
%# codegen
persistent StateNetwork
MATname = 'stateFcnData';
if isempty(StateNetwork)
StateNetwork = coder.load(MATname);
end
out = [x;u];
J = eye(length(out));
% hidden layer #1
Jfc = StateNetwork.fc1.Weights;
out = StateNetwork.fc1.Weights*out + StateNetwork.fc1.Bias;
Jac = deep.internal.coder.jacobian.tanh(out);
out = deep.internal.coder.tanh(out);
J = Jac*Jfc*J;
% hidden layer #2
Jfc = StateNetwork.fc2.Weights;
out = StateNetwork.fc2.Weights*out + StateNetwork.fc2.Bias;
Jac = deep.internal.coder.jacobian.tanh(out);
out = deep.internal.coder.tanh(out);
J = Jac*Jfc*J;
% output layer
J = StateNetwork.output.Weights*J;
J(:,1:1) = J(:,1:1) + eye(1);
% generate Jacobian matrices
A = J(:,1:1);
B = J(:,2:3);
内容を見てみると、重みやバイアス、活性化関数等のネットワーク構造によってダイナミクスが表現されていることが分かります。
MPCの設計
前節までで予測モデルの準備ができたので、続いてMPCの設計を進めていきます。
今回モデルがニューラルネットワークの構造となるため、一種の非線形システムとなります。
ゆえに、MPCは非線形なMPCとして解いていくことになります。
今回はModel Predictive Control Toolboxから提供されているMultistage Nonlinear MPCを使って非線形MPCを設計していきます。
ここで補足しておくと、本Toolboxからは非線形MPCとして、Multistage Nonlinear MPCと(Generic) Nonlinear MPCという2つの非線形MPCが登場しています。
同じ非線形MPCのフレームワークを提供している２つの機能ですが、違いとして(Generic) Nonlinear MPCは、予測ホライズンに亘って一つの評価関数で最適化問題を構築するのに対し、Multistage Nonlinear MPCは、各ステップ（ステージ）毎に評価関数を個別に指定できる構成となっています。
その他の違いは、以下のヘルプに詳細があります。
今回はこのMultistage Nonlinear MPCを使って以下のような非線形最適制御問題を解いていきます。
% 今回は予測ホライズンを48ステップ（4時間）に設定します。
p = (4*3600)/Ts;
% 非線形MPCコントローラ msobj を作成します。
% 今回のモデルでは、状態関数の最初の入力(MD:測定外乱)が外気温度で、2番目の入力(MV)がヒーター入力に該当します。
msobj = nlmpcMultistage(p,1,'MV',2,'MD',1);
% MPCの制御周期は300秒（5分）です。
msobj.Ts = Ts;
% MPCによる操作量の決定において、操作量の変化率も加味した定式化を有効にします。
msobj.UseMVRate = true;
% 状態量、操作量、外乱にそれぞれ名前を指定します。（任意）
msobj.States.Name = "T_room";
msobj.MV.Name = "P_heat";
msobj.MD.Name = "T_atm";
% 非線形MPCの予測モデルとして生成したNeural State-Spaceによる状態関数を指定します。
msobj.Model.StateFcn = "stateFcn";
% この予測モデルの形式が離散時間であることを指定します。
msobj.Model.IsContinuousTime = false;
% コスト関数を指定します。
% コスト関数はステージ毎にそれぞれ独立して設定することが可能です。
% なお、今回はヒーター操作に伴う電気料金を主な評価値に据えています。
% 最初のステージ
msobj.Stages(1).CostFcn = "myCostFcnwSlackStageFirst";
msobj.Stages(1).ParameterLength = 1; % 外部パラメーターの数
% 最後のステージ
msobj.Stages(p+1).CostFcn = "myCostFcnwSlackStageLast";
msobj.Stages(p+1).ParameterLength = 1; % 外部パラメーターの数
msobj.Stages(p+1).SlackVariableLength = 2; % スラックの数
% 中間のステージ
for ct = 2:p
msobj.Stages(ct).SlackVariableLength = 2; % スラックの数
msobj.Stages(ct).ParameterLength = 1; % 外部パラメーターの数
msobj.Stages(ct).CostFcn = "myCostFcnwSlack";
end
% ステージ毎に室内温度に関する不等式制約を指定します。
for ct = 2:p+1
msobj.Stages(ct).IneqConFcn = "myIneqConFunctionwSlack";
end
% 操作変数の上下限制約として最小値（0）と最大値（1）を設定します。
msobj.MV.Min = 0;
msobj.MV.Max = 1;
コスト関数はステージの最初、中間、最後の3つで以下の様に使い分けて指定しています。
なお、priceは電気料金の単価情報であり、MPCの外部パラメーターという扱いです。
function J = myCostFcnwSlackStageFirst(stage, x, mv, dmv, price)
% stage : Stage number
% x : State
% mv : Input
% dmv : Input rate
% price : extra parameter about price
J = 1e4*price*mv(2) + 0.05*(dmv^2);
end
function J = myCostFcnwSlack(stage, x, mv, dmv, e, price)
% stage : Stage number
% x : State
% mv : Input
% dmv : Input rate
% e : Slack
% price : extra parameter about price
J = 1e4*price*mv(2) + 1e5*(e(1) + e(2)) + 0.05*(dmv^2);
end
function J = myCostFcnwSlackStageLast(stage, x, mv, dmv, e, price)
% stage : Stage number
% x : State
% mv : Input
% dmv : Input rate
% e : Slack
% price : extra parameter about price
J = 1e5*(e(1) + e(2));
end
不等式制約については、室内温度に関して快適温度の範囲＋αということで$20-e≦x≦22+e$と与えます。（ソフト制約）
function C = myIneqConFunctionwSlack(stage, x, mv, dmv, e, price)
% stage : Stage number
% x : State
% mv : Input
% dmv : Input rate
% e : Slack
% price : extra parameter about price
C = [-x+20-e(1); x-22-e(2)];
end
非線形MPCでは非線形計画問題（NLP）を解くためにSQP（逐次2次計画法）を使用します。
ゆえに毎イタレーションにおいてコスト関数やモデル等のヤコビアン（勾配）の情報が必要となりますが、それらの情報をユーザーが指定しない限り、ツールでは差分近似処理によって逐次的に求めようとします。
そのため、計算におけるオーバーヘッドとなり、最適化計算の低速化を招くことになります。
そこで、Multispage Nonlinear MPCではR2023aより自動微分(AD：Automatic Differentiation)の機能をサポートするようになりました。
自動微分とはシンボリックな数式処理に基づき、関数の勾配量を自動で正確に求めるための手法です。
Multispage Nonlinear MPCでは、評価関数、予測モデル、制約条件に関するヤコビアンをADによって自動生成する機能をサポート関数generateJacobianFunctionによって使用することができるようになっています。
msobj = generateJacobianFunction(msobj,"state"); % 状態関数に関するAD
msobj = generateJacobianFunction(msobj,"cost"); % 評価関数に関するAD
msobj = generateJacobianFunction(msobj,"ineqcon"); % 不等式制約に関するAD
シミュレーション
それでは、Simulinkを使ってシミュレーションを行ってみます。
なお、住宅暖房のプラントモデルとして、前節で設計したNeural State-Spaceのモデルを使用しています。
下図がシミュレーション結果です。
上段から温度、ヒーター操作量、電気料金の単価の24時間の履歴となっています。
温度の結果を見ると、外気温の変動があるにも関わらず、室内温度は快適な温度範囲である20～22[degC]内にほぼ収まっていることが分かります。（途中の逸脱はソフト制約による影響）
また、日中10時頃を境に電気料金の単価が上昇しますが、事前にヒーター操作量を増加させて室内温度を22[degC]付近まで上げておくことがこの場合における最適なパターンであるとMPCが判断していることが分かります。人の直感とも一致した解です。
最後に
ニューラルネットワークは非線形システムをユニバーサルに近似する手段として、古くから提案されています。
また制御系設計においても今回のような予測モデルの用途に適用したり、コントローラーの学習やゲインのチューニングにおいて用いる等の研究がされてきました。
しかしながら、実用に向けた展開がなかなか進んできていないのが事実かと思います。
ニューラルネットワークもさることながら先進手法に大事なことは信頼性およびロバスト性です。
様々なシナリオを想定したシミュレーションによる入念な検証に加えて、最終的には実機による試験が重要であることは間違いありません。
参考