カスタムの制約と目的関数を使用した操作点の計算

スクリプトを開く

通常、最適化ベースの探索を使用して Simulink® モデルの定常状態の操作点を計算するには、既知の固定値または範囲を指定してモデルの状態、入力、または出力を制約します。ただし、一部のシステムやアプリケーションでは、最適化探索のパラメーター定義においてさらに柔軟性が必要となります。

そのようなシステムの場合、カスタムの制約または追加の最適化目的関数、あるいはその両方を指定することができます。ソフトウェアは、定常状態の操作点を計算するときに、標準の状態、入力、および出力仕様に加えてこれらのカスタムの制約と目的関数を適用します。

カスタムの等式制約と不等式制約は、モデルの状態、入力、および出力の代数的組み合わせとして指定できます。これらの制約では、入力、出力、および状態間の既知の関係を指定することで操作点の探索空間を制限できます。たとえば、モデルのある状態が他の 2 つの状態の和となるように指定できます。

また、カスタムのスカラー目的関数を、モデルの状態、入力、および出力の代数的組み合わせとして指定することもできます。目的関数を使用して、アプリケーションの要件に基づき定常状態の操作点を最適化できます。たとえば、可能性のある平衡点がモデルに複数あるとします。目的関数を指定して、入力エネルギーが最小の状態安定点を検出できます。

複雑なモデルの場合、モデルの入力、出力、および状態のサブセットを選択してカスタムのコスト関数および制約関数へ渡す、カスタムのマッピング関数を指定できます。

モデルを平衡化する場合、次の方法でカスタムの最適化関数を指定できます。

コマンドライン: operspec を使用して操作点の仕様を作成し、仕様の CustomConstrFcn、CustomCostFcn、および CustomMappingFcn プロパティを使ってカスタム関数を指定します。
定常状態マネージャーの使用: [仕様] タブで、[平衡化オプション] をクリックします。[平衡化オプション] ダイアログボックスの [カスタム最適化関数] セクションで、関数名を指定します。
モデル線形化器の使用: [線形解析] タブの [操作点] ドロップダウンリストで、[モデルの平衡化] をクリックします。[モデルの平衡化] ダイアログボックスで、[オプション] タブの [カスタム最適化関数] セクションで関数名を指定します。

次の例では、カスタム最適化関数を作成する方法と、コマンドラインでこれらのカスタム関数を使用してモデルを平衡化する方法を説明します。

Simulink モデル

この例では、オリフィスで相互に連結された 3 つのタンクをもつモデルを使用します。

mdl = 'scdTanks';
open_system(mdl)

Tank1 と Tank2 間のフローは望ましい状態です。Tank2 と Tank3 間のフローは望ましくない回避不可能な漏れです。

このシステムの予想される定常状態は次のようになります。

Tank1 と Tank2 の圧力が同じである。
Tank2 と Tank3 の間に負荷を補償するほぼ一定の圧力差 1 がある。

Tank1 と Tank2 の連結が弱いため、Tank1 と Tank2 の圧力が等しくなるようにモデルを平衡化するのは困難です。

カスタマイズなしでのモデルの平衡化

モデルの既定の操作点の仕様を作成します。仕様は 3 つすべてのタンク圧力を、平衡化された操作点において定常状態になければならない自由な状態として設定します。

opspec = operspec(mdl);

モデルを平衡化するオプションセットを作成し、操作点探索レポートのコマンドウィンドウを非表示に設定します。具体的な平衡化オプションはアプリケーションによって異なります。この例では、非線形最小二乗最適化を使用します。

opt = findopOptions('OptimizerType','lsqnonlin');
opt.DisplayReport = 'off';

モデルを平衡化し、平衡化されたタンク圧力を表示します。

[op0,rpt0] = findop(mdl,opspec,opt);
op0.States

ans = 

     <strong>x</strong>  
    <strong>____</strong>

(1.) scdTanks/Inertia
       0
(2.) scdTanks/Tank1
       9
(3.) scdTanks/Tank2
     9.5
(4.) scdTanks/Tank3
    10.5

Tank1 と Tank2 の平衡化された圧力は一致しません。したがって、既定の操作点の仕様では、期待される定常状態の要件を満たす操作点を見つけることができません。制約の許容誤差 opt.OptimizationOptions.TolCon を小さくした場合、Tank2 と Tank3 の間に漏れがあるため、実行可能な定常状態解が得られません。

カスタムの制約の追加

カスタムの制約を指定するには、現在の作業フォルダーか MATLAB パス上に、次の入力引数をもつ関数を定義します。

x - 操作点の仕様の状態。ベクトルとして指定します。
u - 操作点の仕様の入力。ベクトルとして指定します。
y - 操作点の仕様の出力。ベクトルとして指定します。

出力引数は次のとおりです。

c_ineq - 平衡化中に c_ineq <= 0 を満たさなければならない不等式制約。ベクトルとして返されます。
c_eq - 平衡化中に c_eq = 0 を満たさなければならない等式制約。ベクトルとして返されます。

c_ineq と c_eq の各要素は単一の制約を指定します。アプリケーションの具体的な制約を、状態、入力、および出力の代数的な組み合わせとして定義します。カスタムの等式制約や不等式制約がない場合は、対応する出力引数を [] として返します。

この例では、予期される定常状態の条件を満たすため、次のカスタムの制約関数を定義します。

function [c_ineq,c_eq] =  myConstraints(x,u,y)
    c_ineq = [];
    c_eq = [x(2)-x(3);     % Tank1 pressure - Tank2 pressure
            x(3)-x(4)+1];  % Tank2 pressure - Tank3 pressure + 1
end

c_eq の最初のエントリは、Tank1 と Tank2 の圧力が同じ値になるよう制約します。2 番目の等式制約は、Tank2 と Tank3 間の圧力低下を定義します。

操作点の仕様にカスタムの制約関数を追加します。

opspec.CustomConstrFcn = @myConstraints;

カスタムの制約を含む修正された操作点の仕様を使用してモデルを平衡化し、平衡化された状態値を表示します。

[op1,rpt1] = findop(mdl,opspec,opt);
op1.States

ans = 

      <strong>x</strong>   
    <strong>______</strong>

(1.) scdTanks/Inertia
         0
(2.) scdTanks/Tank1
    9.3333
(3.) scdTanks/Tank2
    9.3333
(4.) scdTanks/Tank3
    10.333

カスタムの制約関数を使ってモデルを平衡化すると、予期したとおりに最初と 2 番目のタンクの圧力が等しい操作点が生成されます。また、想定のとおり、3 番目と 2 番目のタンク間には 1 の圧力差があります。

指定された制約の最終値を調べるには、操作点探索レポートの CustomEqualityConstr プロパティと CustomInequalityConstr プロパティをチェックできます。

rpt1.CustomEqualityConstr

ans =

   1.0e-06 *

   -0.0001
   -0.1540

ゼロに近い値は等式制約が満たされていることを示します。

カスタムの目的関数の追加

カスタムの目的関数を指定するには、カスタムの制約関数と同じ入力引数 (x、u、および y) と、出力引数 F をもつ関数を定義します。F は、平衡化中に最小化される目的関数値で、スカラーとして返されます。

アプリケーション用の目的関数を、状態、入力、および出力の代数的な組み合わせとして定義します。

この例では、Tank3 の圧力を [16,20] の範囲内に保持すると仮定します。ただし、この条件が常に可能であるとは限りません。したがって、厳密な制約を課すのではなく、圧力が [16,20] の範囲に収まらない場合にペナルティを課す目的関数を追加します。そのためには、次のカスタムの目的関数を定義します。

function F =  myObjective(x,u,y)
    F = max(x(4)-20, 0) + max(16-x(4), 0);
end

操作点の仕様オブジェクトにカスタムの目的関数を追加します。

opspec.CustomObjFcn = @myObjective;

カスタムの制約と目的関数の両方を使用して操作点を平衡化し、平衡化された状態値を表示します。

[op2,rpt2] = findop(mdl,opspec,opt);
op2.States

ans = 

    <strong>x</strong> 
    <strong>__</strong>

(1.) scdTanks/Inertia
     0
(2.) scdTanks/Tank1
    15
(3.) scdTanks/Tank2
    15
(4.) scdTanks/Tank3
    16

平衡化された操作点では、Tank3 の圧力がカスタムの目的関数で指定された [16,20] の範囲内にあります。

スカラー目的関数の最終値を表示するには、操作点探索レポートの CustomObj プロパティをチェックします。

rpt2.CustomObj

ans =

     0

カスタムのマッピングの追加

複雑なモデルの場合、カスタムのマッピングを定義してモデルの状態、入力、および出力の部分集合を選択し、カスタムの制約および目的関数に渡すことができます。これを行うと、不要な状態、入力、および出力が排除され、制約と目的関数が簡素化されます。

カスタムのマッピングを指定するには、操作点の仕様 opspec を入力引数とし、以下を出力引数とする関数を定義します。

indx - マッピングされた状態のインデックス
indu - マッピングされた入力のインデックス
indy - マッピングされた出力のインデックス

ブロックパスと状態名に基づいて状態、入力、および出力のインデックスを取得するには、getStateIndex、getInputIndex、およびgetOutputIndexを使用します。これらのコマンドの使用は、モデルの状態の追加など、将来のモデル変更に対してロバストです。あるいは、インデックスを手動で指定することもできます。indx、indu、および indy の形式の詳細については、getStateIndex、getInputIndex、および getOutputIndex を参照してください。

カスタムの制約と目的関数によって使用される状態、入力、出力がない場合は、対応する出力引数を [] として返します。

この例では、3 つのタンクの圧力状態のみを含むマッピングを作成します。そのためには、次のカスタムのマッピング関数を定義します。

function [indx,indu,indy] = myMapping(opspec)
    indx = [getStateIndex(opspec,'scdTanks/Tank1');
            getStateIndex(opspec,'scdTanks/Tank2');
            getStateIndex(opspec,'scdTanks/Tank3')];
    indu = [];
    indy = [];
end

操作点の仕様にカスタムのマッピングを追加します。

opspec.CustomMappingFcn = @myMapping;

カスタムのマッピング関数を使用する場合、カスタムの制約および目的関数内の状態、入力、および出力のインデックスは、マッピング関数に指定された順序に対して相対的でなければなりません。カスタムの制約と目的関数を新しいマッピングで更新します。

function [c_ineq,c_eq] =  myConstraintsMap(x,u,y)
    c_ineq = [];
    c_eq = [x(1)-x(2);     % Tank1 pressure - Tank2 pressure
            x(2)-x(3)+1];  % Tank2 pressure - Tank3 pressure + 1
end

function F =  myObjectiveMap(x,u,y)
    F = max(x(3)-20, 0) + max(16-x(3), 0);
end

ここで x、u、および y は、それぞれマッピングされた状態、入力、出力のベクトルです。これらのベクトルは、indx、indu、および indy で指定されるマッピングされた値をそれぞれ含みます。

更新されたカスタム関数を操作点の仕様に追加します。

opspec.CustomConstrFcn = @myConstraintsMap;
opspec.CustomObjFcn = @myObjectiveMap;

カスタムマッピングを使用してモデルを平衡化し、平衡化された状態を表示します。これは op2 での前の結果と一致します。

[op3,rpt3] = findop(mdl,opspec,opt);
op3.States

ans = 

    <strong>x</strong> 
    <strong>__</strong>

(1.) scdTanks/Inertia
     0
(2.) scdTanks/Tank1
    15
(3.) scdTanks/Tank2
    15
(4.) scdTanks/Tank3
    16

カスタム関数への解析勾配の追加

計算の速度や信頼性を高めるため、カスタムの制約および目的関数に解析勾配を追加できます。勾配を追加することで最適化時の関数呼び出しの回数を減らし、最適化の結果の精度を改善できることがあります。勾配を指定する場合、これはカスタムの制約と目的関数の両方に指定しなければなりません (カスタムの平衡化での勾配は Simscape™ モデルではサポートされていません)。

特定の制約または目的関数の勾配を定義するには、その関数の、特定の状態、入力、または出力についての導関数を使用します。たとえば、次の目的関数があるとします。

F = (u(1)+3)^2 + y(1)^2

この場合、u(1) についての F の勾配は次のとおりです。

G = 2*(u(1)+3)

カスタムの制約関数に勾配を追加するには、次の追加の出力引数を指定します。

G_ineq - 不等式制約の勾配配列
G_eq - 等式制約の勾配配列

G_ineq と G_eq の各列に 1 つの制約の勾配が含まれ、列の順序は対応する制約ベクトルの行の順序に一致します。G_ineq と G_eq の両方で、行数は x、u、および y の状態、入力、および出力の合計数に等しくなります。各列には x の状態についての勾配、続いて u の入力と y の出力についての勾配が含まれます。

この例では、カスタムマッピングを使用する制約関数に勾配を追加します。勾配を使用するときにカスタムマッピングを指定する必要はありません。ただし、状態、入力、および出力のマッピングされた部分集合を使用すると、勾配の定義がより簡単です。

function [c_ineq,c_eq,G_ineq,G_eq] =  myConstraintsGrad(x,u,y)
    c_ineq = [];
    c_eq = [x(1)-x(2);     % Tank1 pressure - Tank2 pressure
            x(2)-x(3)+1];  % Tank2 pressure - Tank3 pressure + 1
    
    G_ineq = [];
    G_eq = [1  0;
            -1 1;
            0 -1]; 
end

この関数では G_eq の i 行目に状態 x(i) についての勾配が含まれています。

同様に、カスタムの目的関数に勾配を追加するには、F の勾配を含む追加の出力引数 G を指定します。G は、G_ineq および G_eq の列と同じ形式の列ベクトルとして返されます。

function [F,G] = myObjectiveGrad(x,u,y)
    F = max(x(3)-20, 0) + max(16-x(3), 0);
    
    if x(3) >= 20
        G = [0 0 1]';
    elseif x(3) <= 16
        G = [0 0 -1]';
    else
        G = [0 0 0]';
    end
end

この例の目的関数は区分的に微分可能なので、G の値は Tank3 の圧力の値によって異なります。

更新されたカスタム関数を操作点の仕様に追加します。

opspec.CustomConstrFcn = @myConstraintsGrad;
opspec.CustomObjFcn = @myObjectiveGrad;

最適化アルゴリズムで勾配を有効にするには、Jacobian 最適化オプションを有効にします。

opt.OptimizationOptions.Jacobian = 'on';

定常状態マネージャーまたはモデル線形化器を使用してモデルを平衡化するときに解析ヤコビアンを使用するには、[解析ヤコビアンを有効化] 平衡化オプションを選択します。

勾配をもつカスタム関数を使用してモデルを平衡化し、平衡化された状態を表示します。

[op4,rpt4] = findop(mdl,opspec,opt);
op4.States

ans = 

    <strong>x</strong> 
    <strong>__</strong>

(1.) scdTanks/Inertia
     0
(2.) scdTanks/Tank1
    15
(3.) scdTanks/Tank2
    15
(4.) scdTanks/Tank3
    16

最適化の結果は、勾配をもたない解の結果と同じです。

勾配によって最適化の効率が改善するかどうかをチェックするには、操作点探索レポートを確認します。たとえば、解に対する関数評価の数を比較します。