% Example 1 in Relative Gain Array for Norm-bounded Uncertain Systems
%              - Kariwala, Skogestad and Forbes (Submitted to I&ECR)

% The gain matrix corresponds to the Wood-Berry distillation column 

% Written by Vinay Kariwala, NTNU, Norway
%            July 4, 2005

% For the alternate uncertainty description, 
% run this program with Delta = [8 0; 1 1] on line 14

G = [12.8 -18.9;6.6 -19.4];

Delta = [2 0;2 0;2 0;2 0;1 1];

Ginv = inv(G);

RGA = G.*inv(G');

alpha = 0.17;

V = alpha*[1 0 1 0;0 1 0 1]';
W = [abs(G(1,:)) zeros(1,2);zeros(1,2) abs(G(2,:))];

e = [1 0];
E = e'*e;

%Permutation matrix to permute Delta_1 to the standard repeated scalar
%uncertainty form 

temp = [1:4 1:4];

[Y,I] = sort(temp);

P = zeros(8,8);

for i = 1:8
    P(I(i),i) = 1;
end

P = [P zeros(8,1);zeros(1,8), 1];

M = [zeros(4,4),-V*E*Ginv*W, V*E*Ginv*e';
     zeros(4,4), -V*Ginv*W, V*Ginv*e';
     e*W, -e*G*E*Ginv*W, e*G*E*Ginv*e'];

M = P'*M*P;

M11 = M(1:end-1,1:end-1); M12 = M(1:1:end-1,end);
M21 = M(end,1:1:end-1); M22 = M(end,end);

% Bisection method 

LB = 0.01; UB = 100000; 
objLB = skewedMu(M11, M12, M21, M22, Delta, LB);
objUB = skewedMu(M11, M12, M21, M22, Delta, UB);
tol = 1e-3;

while (UB-LB) > tol

    Sol = 0.5*(LB+UB);
    objSol = skewedMu(M11, M12, M21, M22, Delta, Sol);
    
    if objSol*objLB < 0
        UB = Sol;
        objUB = objSol;
    else
        LB = Sol;
        objLB = objSol;
    end
end
   
gamma = LB;
