%     REATOR CSTR - DETERMINAÇÃO DO CONTROLE ÓTIMO 
%
%     8.  OTIMIZAÇÃO DE PROCESSOS DINÂMICOS
%
%     8.3 Controle Ótimo
% 					
%
%     Giovani Tonel (giovani.tonel@ufrgs.br) - Outubro de 2006 
%
% Matlab

clc
help run_cstr

disp(' ')
disp(' ')
disp('                Você vai necessitar do arquivo reator_simulink_optim.mdl     ');
disp('                e do arquivo reator_s_function_optim.m !                     ');
disp(' ')
disp('                 Pausa: Pressione qualquer tecla p continuar...'); pause;
disp(' ')

clc
clear all

ip=0; % ip=1 necessita do arquivo CSTR_controller.rlt


addpath D:\Giovani\Otimiza\dynamic_otimization

% res=emso2ml('CSTR_controller.rlt');

disp(' ');
disp('        =====================================================');
disp('        xxxxxx                                       xxxxxxxx');
disp('        xxxxxx   Dados exportados do Simulador EMSO  xxxxxxxx');
disp('        xxxxxx                                       xxxxxxxx');
disp('        =====================================================');

disp(' ');
disp('                           Reação A --> B');
disp(' ');
%=========================================================================
% Bruna's implementation

	syms Fe Cv Te u1ss hss Area Cae Cass ko Ea R Tss U At Twss ro Cp Hr u2ss Tmax Tmin

% Equations   
  	dh = (Fe - Cv * (1 - u1ss) * sqrt (hss))/Area;
   
   dCa = (Cae - Cass) * Fe / (Area * hss) - ko * Cass * exp(-Ea/(R*Tss));
   
   dT = (Te - Tss) * Fe /(Area * hss) - U * At*(Tss - (u2ss *(Tmax - Tmin)+Tmin))/...
      (ro * Area * hss * Cp) + (-Hr)*(ko * exp(-Ea/(R*Tss)) * Cass)/(ro * Cp );

      
   % Making substituitions
   
   % Tw = u2 *(Tmax-Tmin)+Tmin
   % dT = subs (dT, Twss, (u2ss *(Tmax - Tmin)+Tmin));
   
% Jacobian matrix   
   
  Aa = jacobian ([dh; dCa; dT], [hss Cass Tss]);
   
  Bb = jacobian ([dh; dCa; dT], [u1ss u2ss]);
   
  C =  (jacobian ([hss Tss], [hss Cass Tss]));
  
  D =  zeros(2,2);
  %=========================================================================
  



%=======================================================================================
%                                 Extraindo os dados do CSTR
%=======================================================================================

%Valores limites para normalizações
Lmax=5;			%m
Lmin=0;			%m
Tmax=700;		%K
Tmin=230;		%K



%estacionario 
% Tss=		res.CSTR.T(end); 	% ou res.CSTR.Outlet.T(end) 
% Cass=	res.CSTR.Ca(end);		% ou res.CSTR.Outlet.Ca(end)
% hss=		res.CSTR.h(end); 	% altura
u1ss=	0;									% res.PIDL.PID_basic.Ports.output(end)
u2ss=	1;									%	res.PIDT.PID_basic.Ports.output(end)
Twss=	u2ss*(Tmax-Tmin)+Tmin;% res.CSTR.Tw(end)	( 700 K)


hss 	= 2.5; 			%m"; 			= hini
Cass = 1.7244; 	%kmol/m^3;	= Caini
Tss 	= 778.56; 	%K";				= Tini

xss=[2.5 1.7244 778.56]';
	
%Variáveís da corrente de alimentação
Cae = 300;  % kmol/m^3
Fe = 3.5; 	% m^3/h
Te = 300;  	% K

%Cae	=		res.CSTR.Inlet.Ca(1); % Cae=300kmol/m^3
%Fe	=		res.CSTR.Inlet.F(1);
%Te	=		res.CSTR.Inlet.T(1);

%Condiçõs inicias necessaqrias para resolver as equções diferencias
hini = 2.5; %m";
Caini = 50; %kmol/m^3;
Tini = 650; %K";
   

%imprimindo
disp(' ');
disp('       ==================================================================');
dados=[Tss Twss Cae Cass Fe hss];
printmat(dados,'CSTR - Estacionário','Valores','Tss Twss Cae Cass Fe hss');
 
disp(' ');
disp('       ==================================================================');
dados2=[Cae Fe Te];
printmat(dados2,'CSTR- Entrada','Valores','Cae Fe Te');

disp(' ');
disp('       ==================================================================');
dados3=[hini Caini Tini];
printmat(dados3,'CSTR- Condições Inicias','Valores','hini Caini Tini');




% ==================================================================================


if ip,
   
res=emso2ml('CSTR_controller.rlt');

figure(1)
hold on
	
t=res.time/3600;
 g1=plot(t,res.CSTR.T,'r-');
 g2=plot(t,res.CSTR.h*200,'b-');
% g1= plot(t(1:250),res.PIDT.PID_basic.Ports.input(1:250),'r-');
% g2= plot(t(1:250),res.PIDT.PID_basic.Ports.output(1:250),'b'); 

legend([g1 g2],'Temperatura','Nivel',0);
xlabel('Tempo[h]');

ylabel('[K] [m*200]');
title('Reator CSTR');
text(1.6e5,600,'Step')
hold off

   
t=res.time/3600;

figure(3); plot(t(1:5),res.PIDL.PID_basic.Ports.output(1:5),'r'); legend('u1'); title('Ação do PIDL'); xlabel('Tempo [horas]');

figure(4); plot(t(1:5),res.PIDT.PID_basic.Ports.output(1:5),'b'); legend('u2'); title('Ação do PIDT'); xlabel('Tempo [horas]');

end

% ======================================================================================




%=======================================================================================
%                                       Outros DADOS
%=======================================================================================

%Parâmetros do CSTR
R 		= 8.3144; 	%kJ/kmol/K
U 		= 300; 			%kJ/h/m^2/K
ro		= 55.56; 		%kmol/m^3 densidade
Cp 	= 70; 			%kJ/kmol/K
Hr  	= -7000; 		%kJ/kmol
Ea  	= 6e4; 			%kJ/kmol
ko  	= 89*3600; 	%1/h
Area = 8; 			%m^2
At  	= 25;			%m^2
Vrss	= 20;			%m^3
Cv		= 2.2136; 	%m^2.5/h


disp(' ');
disp('       ====================================================================');
dados3=[R U ro Cp Hr Ea ko Area At];
printmat(dados3,'Parâmetros do CSTR','Valores','R U ro Cp Hr Ea ko Area At');

disp(' ');
disp('       *********  Assumido Area de troca térmica (At) constante !!!!  **** ')
disp(' ');
disp('       ====================================================================');
dados4=[Lmax Lmin Tmax Tmin];
printmat(dados4,'Valores limites p/ normalizações','Valores','Lmax Lmin Tmax Tmin');


disp(' ');
disp('       ====================================================================');

disp(' ')   
disp('       ***************    Ver Arquivos do Maple!!!!  ********************* ')
disp(' ')
disp('       ---------------------------Matriz A----------------------------------');


A = eval (Aa)

disp(' ')
disp('       ---------------------------Matriz B----------------------------------');



B = eval (Bb)

disp(' ')

disp('       ---------------------------Matriz C----------------------------------');

C=eval(C)

disp(' ')

disp('       ---------------------------Matriz D----------------------------------');

D

disp(' ')

disp('       ---------------------------Matriz Q----------------------------------');

Q=[1 0 0; 0 1 0; 0 0 1]

disp(' ')
disp('       ---------------------------Matriz R----------------------------------');

R_=[1 0 ; 0 1 ]
disp(' ')

% ================================================================================



disp('       ---------------------------------------------------------------------');
disp(' ')
disp('                 Montando o Problema de Otimização!');
disp(' ')
disp('                 Pausa: Pressione qualquer tecla p continuar...'); pause;
disp('       ---------------------------------------------------------------------');
disp(' ')

disp('                 Resolvendo via LQR - Linear-quadratic regulator design for continuous-time systems.');
disp(' ')

disp('       ---------------------------Help LQR----------------------------------');
help lqr
disp(' ')
disp('       ------------------------------------------------------------------------------------------');
disp(' ')
disp('       Para o nosso problema R=W, N=0, S=P e') 
disp('       A-B*K = A_k --> A matriz A_k caracteriza se o sistema será instável ou não!')
disp('                          -1')
disp('       P*A + A''P - (P*B)R  (B''P) + Q = 0 --> dP/dt=0')
disp('       A e B são invariantes no tempo(=ctes)')
disp(' ')
disp('       --------------------------------------------------------------------------------------------');

disp('                 Pausa: Pressione qualquer tecla p continuar...'); pause;

disp(' ')
disp('       ---------------------------Resolvendo LQR----------------------------------');

[K,P,E] = lqr(A,B,Q,R_)

disp('       ---------------------------Matriz Ak----------------------------------');
disp(' ')

Ak= A-B*K


disp(' ')

if E(:,1)~=0,
   disp('       Todos os auto-valores (E) da matriz Ak são negativos, portanto o sistema é instável!!!!')
end
disp(' ')
disp('                 Pausa: Pressione qualquer tecla p continuar...'); pause;

disp(' ')

disp('       Sabe-se que o Controle Ótimo é dado por: ')
disp(' ')
disp('       u=-Kx e dx/dt= Ax-Bu  --> dx/dt= (A-B*K)x = Ak*x')
disp('       Espaço Estados:  dx/dt= Ax(t)-Bu(t)  e y(t)= Cx(t)+ Du(t)')
disp(' ')



disp(' ')
disp('                 Pausa: Pressione qualquer tecla p continuar...'); 
pause;
disp(' ')
disp('       Matrix de Espaço Estados Continua: SYS = SS(A,B,C,D) ')

disp('       ---------------------------Matriz SYS----------------------------------');


sys_ = ss(A,B,C,D)
disp(' ')
disp('       ---------------------------Matriz de Transferência G-------------------');

G=tf(sys_)

if ip,
   figure(5); Step(G); title('Resp. Degrau U(1)=h U(2)=T'); 
end


disp(' ')
disp('       ---------------------------Matriz de Ganhos---------------------------');


%Ganho= det(dcgain(sys_))
K_=dcgain(sys_)

disp(' ')
disp('       ---------------------------Matriz RGA---------------------------------');

K_=dcgain(sys_);
RGA=K_.*(inv(K_).')


disp(' ')
disp('       ----------------Simulando o Controle Ótimo (uot= -K*x)----------------');

disp(' ')

[t, x, uot]=sim('reator_simulink_optim', 1);
disp(' ')

% x=xot.signals.values; % sem desvio

uot= (-K*(x'));

figure(6); plot(t,uot(1,:)','r'); legend('u1*'); title('Controle Ótimo'); xlabel('Tempo [horas]');
figure(7); plot(t,uot(2,:)','b'); legend('u2*'); title('Controle Ótimo'); xlabel('Tempo [horas]');

   
   

disp('       -------------------------THE END !-------------------');