多目的ゴール到達の最適化

この例では、多目的ゴール到達法を使って極配置の問題を解く方法を説明します。このアルゴリズムは、関数 fgoalattain によって実装します。

システムの変化を記述する方程式

2 入力、2 出力の不安定なプラントを考えます。システム $x (t)$ の変化を記述する方程式は以下のとおりです。

$\frac{d x}{d t} = A x (t) + B u (t),$

ここで、 $u (t)$ は入力 (制御) 信号です。このシステムの出力は以下のとおりです。

$y (t) = C x (t) .$

行列 $A$ 、 $B$ 、 $C$ は以下のとおりです。

A =  [ -0.5  0  0;  0  -2  10;  0  1  -2 ];
B =  [ 1  0;  -2  2;  0  1 ];
C =  [ 1  0  0;  0  0  1 ];

最適化の目的

次のように、制御信号 $u (t)$ が出力 $y (t)$ に比例すると設定する場合について考えます。

$u (t) = K y (t)$

ここで $K$ は行列です。

これは、システム $x (t)$ の変化が以下のようになることを意味します。

$\frac{d x}{d t} = A x (t) + B K C x (t) = (A + B K C) x (t) .$

最適化の目的は、以下の 2 つの特性をもつように $K$ を設計することです。

1. $(A + B K C)$ の固有値の実数部が [-5, -3, -1] より小さい (これは、制御分野で極配置と呼ばれています)。

2. abs( $K$ ) <= 4 ( $K$ の各要素が -4 ～ 4 の範囲である)。

最適化を解くために、まず多目的ゴールを設定します。

goal = [-5, -3, -1];

ゴールに等しい重みを設定して、ゴールの劣到達または過到達が同じ割合になるようにします。

weight = abs(goal);

出力フィードバックコントローラーを初期化します。

K0 = [ -1 -1; -1 -1];

コントローラーに上限と下限を設定します。

lb = repmat(-4,size(K0))

lb = 2×2

    -4    -4
    -4    -4

ub = repmat(4,size(K0))

反復ごとに出力が与えられるように、最適化表示パラメーターを設定します。

options = optimoptions('fgoalattain','Display','iter');

閉ループシステムの固有値を返すベクトル値関数 eigfun を作成します。この関数には、追加のパラメーター (つまり、行列 $A$ 、 $B$ 、および $C$ ) が必要です。これらのパラメーターを渡す最も都合の良い方法は、次のように無名関数を介して渡すというものです。

eigfun = @(K) sort(eig(A+B*K*C));

最適化ソルバーの呼び出し

最適化を開始するために fgoalattain を呼び出します。

[K,~,attainfactor] = ...
        fgoalattain(eigfun,K0,goal,weight,[],[],[],[],lb,ub,[],options);

                 Attainment        Max     Line search     Directional 
 Iter F-count        factor    constraint   steplength      derivative   Procedure 
    0      6              0       1.88521                                            
    1     13          1.031       0.02998            1           0.745     
    2     20         0.3525       0.06863            1          -0.613     
    3     27        -0.1706        0.1071            1          -0.223    Hessian modified  
    4     34        -0.2236       0.06654            1          -0.234    Hessian modified twice  
    5     41        -0.3568      0.007894            1         -0.0812     
    6     48        -0.3645      0.000145            1          -0.164    Hessian modified  
    7     55        -0.3645             0            1        -0.00515    Hessian modified  
    8     62        -0.3675     0.0001548            1        -0.00812    Hessian modified twice  
    9     69        -0.3889      0.008327            1        -0.00751    Hessian modified  
   10     76        -0.3862             0            1         0.00568     
   11     83        -0.3863     4.469e-13            1          -0.998    Hessian modified twice  

Local minimum possible. Constraints satisfied.

fgoalattain stopped because the size of the current search direction is less than
twice the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.

解における制御パラメーターの値は以下のとおりです。

K = 2×2

   -4.0000   -0.2564
   -4.0000   -4.0000

閉ループシステムの固有値 eigfun(K) は、次のとおりです(出力引数 fval にも格納されています)。

eigfun(K)

ans = 3×1

   -6.9313
   -4.1588
   -1.4099

到達因子は、ゴール到達のレベルを示します。負の到達因子は過到達を示し、正の到達因子は劣到達を示します。この実行で得られた値 attainfactor は、目的がほぼ 40% 過到達されたことを意味します。

attainfactor

attainfactor = -0.3863

ODE の解でシステムの変化を示す

ここでは、計算されたフィードバック行列 $K$ を使用して、点 x(0) = [1;1;1] から始めて、システム $x (t)$ が時刻 0 から時刻 4 までどのように変化するかを示しています。

まず、微分方程式を解きます。

[Times, xvals] = ode45(@(u,x)((A + B*K*C)*x),[0,4],[1;1;1]);

次に、結果をプロットします。

plot(Times,xvals)
legend('x_1(t)','x_2(t)','x_3(t)','Location','best')
xlabel('t');
ylabel('x(t)');

ゴールを正確に達成するように設定する

固有値をできるだけゴール値 [-5, -3, -1] に近くすることが必要と仮定します。options.EqualityGoalCount を、できる限りゴールに近い (つまり、過到達にしない) 目的の数に設定します。

3 つの目的がすべて、ゴールになるべく近いことが望ましいです。

options.EqualityGoalCount = 3;

最適化ソルバーの呼び出し

最適化ソルバーを呼び出す準備が整いました。

[K,fval,attainfactor,exitflag,output,lambda] = ...
    fgoalattain(eigfun,K0,goal,weight,[],[],[],[],lb,ub,[],options);

                 Attainment        Max     Line search     Directional 
 Iter F-count        factor    constraint   steplength      derivative   Procedure 
    0      6              0       1.88521                                            
    1     13          1.031       0.02998            1           0.745     
    2     20         0.3525       0.06863            1          -0.613     
    3     27         0.1528     -0.009105            1           -0.22    Hessian modified  
    4     34        0.02684       0.03722            1          -0.166    Hessian modified  
    5     41      1.041e-17      0.005702            1          -0.116    Hessian modified  
    6     48      1.459e-18     9.699e-06            1       -7.73e-16    Hessian modified  
    7     55     -5.082e-23     4.969e-11            1       -7.65e-14    Hessian modified  

Local minimum possible. Constraints satisfied.

fgoalattain stopped because the size of the current search direction is less than
twice the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.

この解における制御パラメーターの値は以下のとおりです。

K = 2×2

   -1.5954    1.2040
   -0.4201   -2.9046

今回の閉ループシステムの固有値 (出力 fval でも保持されます) は、以下のようになります。

eigfun(K)

ans = 3×1

   -5.0000
   -3.0000
   -1.0000

到達因子は、ゴール到達のレベルを示します。負の到達因子は過到達を示し、正の到達因子は劣到達を示します。得られた到達因子が低いことは、固有値がほぼ正確にゴールを満たしていることを示しています。

attainfactor

attainfactor = -5.0824e-23

ODE の解で新しいシステムの変化を示す

ここでは、計算された新しいフィードバック行列 $K$ を使用して、点 x(0) = [1;1;1] から始めて、システム $x (t)$ が時刻 0 から時刻 4 までどのように変化するかを示しています。

まず、微分方程式を解きます。

[Times, xvals] = ode45(@(u,x)((A + B*K*C)*x),[0,4],[1;1;1]);

次に、結果をプロットします。

plot(Times,xvals)
legend('x_1(t)','x_2(t)','x_3(t)','Location','best')
xlabel('t');
ylabel('x(t)');

参考

fgoalattain

多目的ゴール到達の最適化

システムの変化を記述する方程式

最適化の目的

最適化ソルバーの呼び出し

ODE の解でシステムの変化を示す

ゴールを正確に達成するように設定する

最適化ソルバーの呼び出し

ODE の解で新しいシステムの変化を示す

参考

関連するトピック