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

% The gain matrix corresponds to a distillation column (Ref: Ogunnaike et
% al, AICHE Journal, 1983)

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

G = [0.66 -0.61 -0.0049;1.11 -2.36 -0.012;-33.68 46.2 0.87];
Ginv = inv(G);

RGA = G.*inv(G');

% Only diagonal pairings are permissible

% Define the principal submatrices 
% Gii is G with ith row and ith column removed
G11 = G(2:3,2:3); G22 = [G(1,1) G(1,3);G(3,1) G(3,3)]; G33 = G(1:2,1:2); 

% Check the relative gains of submatrices
RGA11 = G11.*inv(G11');
RGA22 = G22.*inv(G22');
RGA33 = G33.*inv(G33');

% Diagonal pairings are perfect for integrity as corresponding relative
% gains for all principal submatrices are positive

% Now find the smallest allowable perturbation of the form g_ij + alpha |g_ij|

V = abs(G);
W = eye(3);

for i = 1:3
    e = zeros(1,3); e(i) = 1; 
    E = e'*e;
    
    N = [-inv(RGA(i,i))*V*E*Ginv*E*W, -V*E*Ginv*(eye(3)-inv(RGA(i,i))*E*G*E*Ginv)*W
     -inv(RGA(i,i))*V*Ginv*E*W -V*Ginv*(eye(3)-inv(RGA(i,i))*E*G*E*Ginv)*W];
 
    specrad(i) = max(abs(eig(N)));  % specrad = 1, 3.3114, 1
end

alpha = 1/max(specrad); % alpha = 0.302

% Check that Delta makes G22 singular.

verifyG22 = det(G22 - alpha*abs(G22))

% Everything else is non-singular (somewhat crude)

delta = -alpha:0.01:alpha;

for j = 1:length(delta)
    verifyG(j) = det(G + delta(j)*abs(G));
    verifyG11(j) = det(G11 + delta(j)*abs(G11));
    verifyG33(j) = det(G33 + delta(j)*abs(G33));
end

min(verifyG)/det(G), min(verifyG11)/det(G11), min(verifyG33)/det(G33) 

% All positive, implying that perturbed G, G11 and G33 are non-singular 