s=tf('s');

%Sigurd exapmle 3.11
%  G=0.01*exp(-5*s)/((s+1.72e-4)*(4.32*s+1))*[-34.54*(s+0.0572) 1.913; -30.22*s -9.188*(s+6.95e-4)]

% Identification based on MPC Case 90%
%  G= [1.65/(2.62*s+1) 0.00968/(0.0483*s+1); 
%     -1.78*(-1.31*s+1)/(54.9*s^2+25.3*s+1)...
%     0.14*(20.7*s+1)/(902*s^2+46.2*s+1)];

%Identification based on Case MPC 95.26%
G= [1.74*(3.94*s+1)/(19.7*s^2+11.4*s+1) 0.038*(1-20*s)/((1+757*s)*(1+32*s));%-0.0265*exp(-2*s)/(102*s^2+20.2*s+1);%;%0.038*exp(-33.33*s)/(1+90.75*s); % 
  -1.02*(9.32*s+1)*exp(-2*s)/(218*s^2+17.3*s+1)...
  0.0754*(6.03*s+1)/(205*s^2+18.8*s+1)];

% Now compute frequency-depende RGA (or rather; RGA-number)
omega =logspace(-3,2,100);   % frequency range
for i=1:length(omega)
    Gf = freqresp(G,omega(i));
    RGAw(:,:,i) = Gf.*inv(Gf).';
    RGAno(i) = sum(sum(abs(RGAw(:,:,i)-eye(2))));
    RGAno_offdiag (i) = sum(sum(abs(RGAw(:,:,i)-[0 1;1 0])));
   %PRGA=diag(diag(Gf))*inv(Gf);
   %PRGA(:,:,i)=PRGA;
end
%RGA=frd(RGAw,omega)
figure(1)
semilogx(omega, RGAno);
hold on
semilogx(omega, RGAno_offdiag,'r'); xlabel('Frequency [rad/min]'),ylabel('||RGA - I||_{sum}');gtext({'Diagonal pairing alt.1';'Off-diagonal pairing alt.2'},'FontSize',15)%legend('Diagonal pairing','Off-diagonal pairing')

