このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。
多目的ゴール到達の最適化
この例では、多目的ゴール到達法を使って極配置の問題を解く方法を説明します。このアルゴリズムは、関数 fgoalattain
によって実装します。
システムの変化を記述する方程式
2 入力、2 出力の不安定なプラントを考えます。システム の変化を記述する方程式は以下のとおりです。
ここで、 は入力 (制御) 信号です。このシステムの出力は以下のとおりです。
行列 、、 は以下のとおりです。
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 ];
最適化の目的
次のように、制御信号 が出力 に比例すると設定する場合について考えます。
ここで は行列です。
これは、システム の変化が以下のようになることを意味します。
最適化の目的は、以下の 2 つの特性をもつように を設計することです。
1. の固有値の実数部が [-5, -3, -1] より小さい (これは、制御分野で極配置と呼ばれています)。
2. abs() <= 4 ( の各要素が -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))
ub = 2×2
4 4
4 4
反復ごとに出力が与えられるように、最適化表示パラメーターを設定します。
options = optimoptions('fgoalattain','Display','iter');
閉ループ システムの固有値を返すベクトル値関数 eigfun を作成します。この関数には、追加のパラメーター (つまり、行列 、、および ) が必要です。これらのパラメーターを渡す最も都合の良い方法は、次のように無名関数を介して渡すというものです。
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
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 の解でシステムの変化を示す
ここでは、計算されたフィードバック行列 を使用して、点 x(0) = [1;1;1] から始めて、システム が時刻 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
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 の解で新しいシステムの変化を示す
ここでは、計算された新しいフィードバック行列 を使用して、点 x(0) = [1;1;1] から始めて、システム が時刻 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)');