MPC tool box内の最適化器を​変更することは可能で​しょうか？

こんにちは。

MPC tool box内で二次計画法(QP solver)を用いて操作量を決定していると思います。

そのQP solverを別の最適化手法に置き換えることは可能でしょうか？

また、現在MPC tool box内のoptimizer(MATLAB function)において

外部関数を用いている箇所がありますが、その外部関数内部はどのように見ることができますか？

質問をまとめると

①MPCの最適化器としてQP以外の手法に置き換えることは可能か。

②MPC内部の外部関数へのアクセス。

以上です。宜しくお願い致します。

Answer 1

[Toshinobu Shintai]

①MPCの最適化器としてQP以外の手法に置き換えることは可能か。

可能です。MPCのオブジェクトを作成した後、以下のようなコマンドを実行すると、自作のQPソルバーを使って計算を行うようになります。

% Simulinkでシミュレーション実行する時の計算を置き換える場合
mpcobj.Optimizer.CustomSolver = true;
% コード生成時に置き換える場合（コード生成で外部に機能を渡す時）
mpcobj.Optimizer.CustomSolverCodeGen = true;

自作QPソルバーは以下のような様式で作成しなければなりません。

function [x, status] = mpcCustomSolver(H, g, A, b, x0)

end

上記関数を先頭に記述したファイルを作成します。そのファイル名は必ず「mpcCustomSolver.m」「mpcCustomSolverCodeGen.m」です。

上記二つのファイルの中身は同一で構いません。それぞれシミュレーション実行時、コード生成時に参照されるファイルとなります。

x, H, g, A, b, x0が何かはドキュメントを参照してください。

https://www.mathworks.com/help/mpc/controller-design.html

statusはデバッグ用の外部出力変数なのですが、大変不便なことにスカラーしか指定できません。とりあえず反復計算の回数を出力させることをお勧めします。

②MPC内部の外部関数へのアクセス。

MATLAB Function ブロック内のスクリプトを開き、例えば

[u, cost, useq, status, iAout] = mpcblock_optimizer(...

となっている所で、「mpcblock_optimizer」を右クリックして開くをクリックすると、その関数に飛ぶことができます。

しかし、一部mex化されているものやpファイルで隠蔽化されているものもあります。それらの関数については公開できないものになっています。申し訳ありません。

アルゴリズムについては、上記のリンク先に詳細が書かれていますので、そちらを参照していただけますでしょうか。

Answer 2

[Youhei Ogawa]

ご回答ありがとうございます。

どのように作成すれば良いか分かりました。

また、自作QPソルバーの部分に関して追加で質問ですが、

[H,f,A,b,x0]から、この部分にL1正則化(LASSO)を作成し、CVXを用いてxを算出することは可能でしょうか？

MPC tool boxの最適化器を二次計画法からスパースモデリングに置き換えたいのですが...

別の方法がありましたら、そちらを教えていただくことはできますでしょうか。

宜しくお願い致します。

[Toshinobu Shintai]

私はやったことがないのでわかりません。

一つ注意事項としては、MATLABの関数の中にはSimulink Coderでコード生成できるものとできないものがあります。cvxがコード生成できない関数である場合は組み込めないということになります。

[Youhei Ogawa]

分かりました。ありがとうございます。

[Kohei Iwamoto]

Simulationについてはcvxを使用することが出来る模様です。

MPC Toolboxで2次計画問題を取り扱う場合には、評価関数及び制約式の定式化は以下のURLに記載の通りとなっていますのでご注意下さい。

QP Solver：

https://jp.mathworks.com/help/mpc/ug/qp-solver.html

ご希望のL1正則化用に定式化を変更して上手く機能するかどうかは、実際のところ試してみないと分かりませんので試行錯誤してください。

試したコードを以下に貼り付けておきますので参考にしてください。

function [X, status] = mpcCustomSolver(H,f,A,b,x0)

persistent n;   %Keep variable 'n' in memory

if isempty(n)
    % Run just first time step
    n = evalin('base','Obj.ControlHorizon');    %Extract variable from base workspace
end

cvx_begin quiet
    variable X(n+1)
    minimize (0.5*X'*H*X+f'*X)
    subject to
        A*X >= b;
cvx_end

% Get cvx status
str = cvx_status;

if strcmp('Solved',str)
    status = 1;     % Feasible solution
else
    status = -1;    % Infeasible solution
end

end

[Youhei Ogawa]

親切な対応ありがとうございます。ぜひ参考にさせていただきたいと思います。