有制約多目的最適化問題で妥協解を求める

t12の計算時にペンローズ擬似逆行列を用いるのがポイント。
パレートフロントはSTOMで希求水準をたくさん設定して求めるのが正攻法なのだろうが、横着してMatlabのboundary()コマンドを利用した。そのため実際はパレート最適値ではない点も拾ってグラフに描画されている。
値を代入するタイミング（つまりsubsコマンドを呼び出すタイミング）によって値が異なるのが気持ち悪い。
見た目は動いてるように見えるが、最初のトレードオフ比Tの計算で、非対角項（つまりトレードオフ比）が負になっている（教科書では正の値）。

（教科書の値）
T =
   -1.0      0.711
   1.405   -1.0
T2 =
   -1.0  1.0
    1.0  -1.0
（計算値）
T =
   -1.0000   -0.6977
   -1.4334   -1.0000
T2 =
   -1.0000    1.0375
    0.9639   -1.0000

clear; close all hidden;

f1 = @(x) x(:,1);
f2 = @(x) x(:,2);
g1 = @(x) 1-x(:,1).^2-x(:,2).^2+0.1*cos(16*atan2(x(:,1),x(:,2)));
g2 = @(x) (x(:,1)-0.5).^2+(x(:,2)-0.5).^2-0.5;

f1i = 0;
f2i = 0;
f1a = 0.4;
f2a = 0.25;

% dx = 0.01;
dx = 0.005;
x1 = 0:dx:pi;
x2 = 0:dx:pi;
length(x1)
[X1,X2] = meshgrid(x1,x2);
X1vec = reshape(X1.',1,[]);
X2vec = reshape(X2.',1,[]);
Xvec = [X1vec;X2vec]';
Xpossible = Xvec(find(g1(Xvec)<=0 & g2(Xvec)<=0),:);

figure(1);
plot(Xpossible(:,1),Xpossible(:,2),'o');

% 希求水準f1a=0.4,f2a=0.25に対するパレート最適値（赤丸の点）
w1 = 1/(f1a-f1i);
w2 = 1/(f2a-f2i);
C = max(w1*(f1(Xpossible)-f1a),w2*(f2(Xpossible)-f2a));
xopt = Xpossible(find(C==min(C)),:);
x1opt = xopt(:,1);
x2opt = xopt(:,2);

f1opt = f1(xopt)
f2opt = f2(xopt)
hold on;
plot(f1opt,f2opt,'ro');
plot(f1a,f2a,'xk');

%トレードオフ比t12の計算
syms f1_symbol(x1,x2) f2_symbol(x1,x2) g1_symbol(x1,x2) g2_symbol(x1,x2);
f1_symbol(x1,x2) = @(x1,x2) x1;
f2_symbol(x1,x2) = @(x1,x2) x2;
g1_symbol(x1,x2) = @(x1,x2) 1-x1^2-x2^2+0.1*cos(16*atan2(x1,x2));
g2_symbol(x1,x2) = @(x1,x2) (x1-0.5)^2+(x2-0.5)^2-0.5;
grad_f1 = gradient(f1_symbol);
grad_f2 = gradient(f2_symbol);
grad_g1 = gradient(g1_symbol);
grad_g2 = gradient(g2_symbol);

N = [grad_g1,grad_g2];
N = simplify(N);
P = eye(2) - N*inv(transpose(N)*N)*transpose(N);%Pの計算はペンローズ擬似逆行列にする
grad_F1 = P*grad_f1;
grad_F2 = P*grad_f2;
grad_F1v = eval(subs(grad_F1,{x1 x2},{x1opt x2opt}));%計算負荷が高すぎるのでここで値を代入する
grad_F2v = eval(subs(grad_F2,{x1 x2},{x1opt x2opt}));%計算負荷が高すぎるのでここで値を代入する
t12 = -transpose(grad_F2v)*pinv(grad_F1v*...
          transpose(grad_F1v))*grad_F1v;%t12の計算もペンローズ逆行列で計算する

t21 = 1/t12;
T = [-1 t12;
    t21 -1]

% パレートフロントを求める
% Matlabのboundary()コマンドを利用して外周の点を広い、f1<1, f2<1となるようにフィルタをかける。
k = boundary(Xpossible(:,1),Xpossible(:,2));
Xk = Xpossible(k,:);
Xb1 = Xk(:,1);
Xb2 = Xk(:,2);
Xc = Xk(find(Xb1<1 & Xb2<1),:);
Xc1 = Xc(:,1);
Xc2 = Xc(:,2);
plot(Xc1,Xc2,'xk');
big;

% トレードオフ比＝１となるようなx1,x2の値（希求水準）を求める

Ds = [];
hw = waitbar(0,'Running...');
for idx=1:length(Xc)
    Xp = Xc(idx,:);
    Xp1 = Xp(1);
    Xp2 = Xp(2);
    try
    grad_F1v = eval(subs(grad_F1,{x1 x2},{Xp1 Xp2}));%計算負荷が高すぎるのでここで値を代入する
    grad_F2v = eval(subs(grad_F2,{x1 x2},{Xp1 Xp2}));%計算負荷が高すぎるのでここで値を代入する
    t12 = -transpose(grad_F2v)*pinv(grad_F1v*...
         transpose(grad_F1v))*grad_F1v;%t12の計算もペンローズ逆行列で計算する
    D = (t12-1)^2;
    Ds(end+1) = D;
%     disp(sprintf('%d %.3g %.3g %.3g %.3g',idx, Xp1, Xp2, t12, D));
    waitbar(idx/length(Xc),hw);
    end
end


figure(2);
plot(Ds,'-');
ylim([-0.1,1]);

Xopt2 = Xc(find(Ds==min(Ds)),:)
f1opt2 = f1(Xopt2)
f2opt2 = f2(Xopt2)

grad_F1v = eval(subs(grad_F1,{x1 x2},{Xopt2(1) Xopt2(2)}));%計算負荷が高すぎるのでここで値を代入する
grad_F2v = eval(subs(grad_F2,{x1 x2},{Xopt2(1) Xopt2(2)}));%計算負荷が高すぎるのでここで値を代入する
t12 = -transpose(grad_F2v)*pinv(grad_F1v*...
          transpose(grad_F1v))*grad_F1v;%t12の計算もペンローズ逆行列で計算する

t21 = 1/t12;
T2 = [-1 t12;
    t21 -1]

figure(1);
plot(f1opt2,f2opt2,'rx');
axis equal;

big;