clear all; close all
% matrizes de ganhos estáticos de uma coluna de destilação binária
% variáveis de saída: composição de topo e fundo do mais leve

LV=[0.878 -0.864;1.082 -1.096]; % entrada: vazão de refluxo de líquido e de vapor
DV=[-0.878 0.014;-1.082 -0.014]; % entrada: vazão de destilado e refluxo de vapor

% funções de transferência
G_LV=tf({LV(1,1), LV(1,2); LV(2,1), LV(2,2)},{[75 1], [75 1]; [75 1], [75 1]});
G_DV=tf({DV(1,1), DV(1,2); DV(2,1), DV(2,2)},{[75 1], [75 1]; [75 1], [75 1]});

% espaço de estados
[A_LV,B_LV,C_LV,D_LV]=ssdata(G_LV);
[A_DV,B_DV,C_DV,D_DV]=ssdata(G_DV);

% ganhos estáticos
G0_LV=-C_LV*(inv(A_LV)*B_LV)+D_LV;
G0_DV=-C_DV*(inv(A_DV)*B_DV)+D_DV;

% estabilidade
eig_LV=eig(A_LV)
eig_DV=eig(A_DV)

% SVD
[U_LV,S_LV,V_LV]=svd(LV)
[U_DV,S_DV,V_DV]=svd(DV)

% condicionamento
cond_LV = cond(LV)
cond_DV = cond(DV)

% ou
S_LV(1,1)/S_LV(2,2)
S_DV(1,1)/S_DV(2,2)

% resposta ao degrau unitário
figure; step(G_LV); title('Resposta ao degrau unitário G_{LV}');
figure; step(G_DV); title('Resposta ao degrau unitário G_{DV}');

% re-escalonamento (condicionamento mínimo)
[gama_LV,L_LV,R_LV,s_LV]=mincond_mu(LV)
[gama_DV,L_DV,R_DV,s_DV]=mincond_mu(DV)
s_G_LV = L_LV * G_LV * R_LV;
s_G_DV = L_DV * G_DV * R_DV;

% SVD do sistema escalonado
[s_U_LV,s_S_LV,s_V_LV]=svd(s_LV)
[s_U_DV,s_S_DV,s_V_DV]=svd(s_DV)

% condicionamento do sistema escalonado
s_cond_LV = cond(s_LV)
s_cond_DV = cond(s_DV)

% resposta ao degrau unitário do sistema escalonado
figure; step(s_G_LV); title('Resposta ao degrau unitário sG_{LV}');
figure; step(s_G_DV); title('Resposta ao degrau unitário sG_{DV}');

% tempo final
tf = 5/max([abs(eig_LV); abs(eig_DV)]);
T = [0:tf]';
n = length(T);

% simulações dinâmicas
rep = ones(n,1);
u1_LV = rep*V_LV(:,1)'; % degrau na direção de maior ganho LV
u2_LV = rep*V_LV(:,2)'; % degrau na direção de menor ganho LV
figure; subplot(121); lsim(G_LV,u1_LV,T);
title('Resposta ao degrau nas direções de maior e menor ganho G_{LV}');
subplot(122); lsim(G_LV,u2_LV,T);
figure; subplot(121); lsim(s_G_LV,u1_LV,T);
title('Resposta ao degrau nas direções de maior e menor ganho sG_{LV}');
subplot(122); lsim(s_G_LV,u2_LV,T);

u1_DV = rep*V_DV(:,1)'; % degrau na direção de maior ganho DV
u2_DV = rep*V_DV(:,2)'; % degrau na direção de menor ganho DV
figure; subplot(121); lsim(G_DV,u1_DV,T);
title('Resposta ao degrau nas direções de maior e menor ganho G_{DV}');
subplot(122); lsim(G_DV,u2_DV,T);
figure; subplot(121); lsim(s_G_DV,u1_DV,T);
title('Resposta ao degrau nas direções de maior e menor ganho sG_{DV}');
subplot(122); lsim(s_G_DV,u2_DV,T);
