% ******** Reator de van de Vusse - caso isotérmico **********************
% ******** Autor: Evandro Alvaristo **************************************
% ******** 19dez2013 *****************************************************
% ******** Rotina principal **********************************************
clear all
close all
clc
% -------- Definição das constantes do modelo (Trierweiler - 1997) -------
% -------- Parâmetros das reações ----------------------------------------
global k10 k20 k30 E1 E2 E3
k10 = 1.287e12;                %[1/h]
k20 = 1.287e12;                %[1/h]
k30 = 9.043e9;                 %[L/(h*mol)]
E1 = -9758.3;                  %[Kelvin]
E2 = -9758.3;                  %[Kelvin]
E3 = -8560;                    %[Kelvin]
% -------- Parâmetros do processo ----------------------------------------
global Cain T
Cain = 5.1;                    %[mol/l]
T = 135;                       %[Celsius]
% -------- Estado de referência ------------------------------------------
global FV_ref Ca_ref Cb_ref
%Lado direito da curva 
% FV_ref = 90;                   %[1/h]
% Ca_ref = 2.8147;               %[mol/h]
% Cb_ref = 1.0461;               %[mol/h]
%Novo
FV_ref = 115;                   %[1/h]
Ca_ref = 3.0870;               %[mol/h]
Cb_ref = 0.9768;               %[mol/h]

% -------- Variáveis do controlador (Parâmetros de sintonia) -------------
global Hp Hc Sm SP Ta 
% Hp = 30; Hc = 10;               %Hp = horizonte de predição, Hc = horizonte de controle
% Hp = 25; Hc = 6;
% Hp = 12; Hc = 3;                %Heuristica 1
Hp = 25; Hc = 5;                 %Estratégia de sintonia Seborg (1989)
Ta = 0.004;                    %[h] - tempo de amostragem
Sm = 0.0008;                    %Supressão de movimento
% ------------------------------------------------------------------------

SP = 0.90;                     %[mol/h] - setpoint atingível
%SP = 1.15;                     %[mol/h] - setpoint inatingível 

DeltaFV(1,:) = [5*ones(1,Hp)];
DeltaFV_min = [-30*ones(1,Hc)];
DeltaFV_max = [65*ones(1,Hc)];

% -------- Condições iniciais --------------------------------------------
global Ca_mod Cb_mod FV d Cai Cbi xotimo ganho
Ca0 = Ca_ref; Cb0 = Cb_ref;    %[mol/l] - valores estacionários para FV=16.3 1/h 
Ca_mod(1) = Ca0;               %[mol/l] - Conc. A do modelo inicial
Cb_mod(1) = Cb0;               %[mol/l] - Conc. B do modelo inicial
t_x(1) = 0;                    %[h]
Ca_planta(1) = Ca0;            %[mol/l] - Valor medido da planta (valor inicial) 
Cb_planta(1) = Cb0;            %[mol/l] - Valor medido da planta (valor inicial)

modeloVDV(:,1) = [Ca0 - Ca_ref, Cb0 - Cb_ref];
nsim = 126;          %55           %número de simulações
xotimo(1) = 0;
xotimo_matriz = zeros(nsim+1,Hc);
xotimo_matriz(1,:) = xotimo(1);

FV = FV_ref;
IT = tril(ones(Hc));

for j = 1:nsim

Cb_medido = Cb_planta(j);      %[mol/l] - Valor medido da planta
Cb_mod_atual = Cb_mod(1);      %[mol/l] - Valor calculado pelo modelo
d = Cb_medido-Cb_mod_atual;    %[mol/l] - distúrbio 
disturbio(j,:) = d;

Cb_medido1(j,:) = Cb_medido;
Cb_mod_atual1(j,:) = Cb_mod_atual;

Cai = Ca_mod(1);               %[mol/l] - Valor de Ca medido da planta
Cbi = Cb_mod(1);               %[mol/l] - Valor de Cb medido da planta

% -------- Otimização para a determinação das ações de controle ----------   
OPTIONS = optimset('algorithm','sqp','display','iter-detailed','MaxFunEvals',1e10);
[xotimo,fval,flag,out] = fmincon(@FobjMPCadaptativo,DeltaFV(j,1:Hc),[],[],[],[],DeltaFV_min,DeltaFV_max,[],OPTIONS); %MPC adaptativo

DeltaFV(j+1,:) = [xotimo(2:Hc) xotimo(Hc)*ones(1,Hp-Hc+1)];       %deslocamento do vetor para a próxima iteração

xotimo_matriz(j,:) = FV + (IT*xotimo')';

FV = FV + xotimo(1);                                          %somente a primeira ação de controle é implementada na planta

% -------- Simulação da planta(EDO) com a primeira ação de controle ------   
tspan = [(j-1)*Ta j*Ta];                
OPTS = odeset('RelTol',1e-6,'AbsTol',1e-6);   
[tresp,respp] = ode15s('reatorVDV_planta',tspan,[Ca0 Cb0],OPTS);
Ca0 = respp(end,1); Cb0 = respp(end,2);
Ca_planta(j+1) = respp(end,1);
Cb_planta(j+1) = respp(end,2);

t_x(j+1) = Ta*j;
% % -------- Cálculo do ganho estático -------------------------------------
% g_num = Cb_planta(j+1)-Cb_planta(j);
% g_den = xotimo_matriz(j+1,1)-xotimo_matriz(j,1);
% g = g_num/g_den;
% ganho(j,:) = g;
end

Cb_set = SP*ones(1,length(Cb_planta));

% -------- Critérios de análise comparativa ------------------------------
% -------- Integral of Squared Erros (ISE)
ISE = ((Cb_set - Cb_planta)*(Cb_set - Cb_planta)')*Ta
% -------- Esforço de manipulação (SU)
u1 = xotimo_matriz(1:nsim,1)';
u2 = [FV_ref u1];
for w = 1:nsim
    dU = (u2(w+1) - u2(w))^2;
    dU_v(w,:) = dU;
end
SU=sum(dU_v)*Ta


%Armazenando dados da simulação no Excel
%DADOS=[TEMPO RESPOSTA_Cb SETPOINT ACAO_CONTROLE DISTURBIO]
tempo = t_x(1:nsim)';
respostaCb = Cb_planta(1:nsim)';
setpointCb = Cb_set(1:nsim)';
AcaoControle = xotimo_matriz(1:nsim)';

dadosMPCadaptativo = zeros(nsim,5);
dadosMPCadaptativo = [tempo respostaCb setpointCb AcaoControle disturbio];
[status,msginfo] = xlswrite('MPCadaptativoLDSA.xls',dadosMPCadaptativo,'planMPCadaptativoLDSA','A1');

figure(1)
hold on
box on
plot(t_x,Ca_planta,'k:','LineWidth',1.5);
plot(t_x,Cb_planta,'k-.','LineWidth',1.5);
title('Conc. de A e B resultantes das ações de controle')
 hleg1 = legend('Conc. A','Conc. B');
    set(hleg1,'Location','SouthEast')
       xlabel('\bfTempo(ao longo do tempo de simulação) - h')
       ylabel('\bfConc. A e B - mol/l')

figure(2)
subplot(2,1,1)
hold on
box on
plot(t_x,Cb_planta,'k-','LineWidth',1.5);
plot(t_x,Cb_set,'r--','LineWidth',1.5);
hold off
    hleg1 = legend('Var. Controlada (y1)','SP');
    set(hleg1,'Location','SouthEast')
       xlabel('\bfTempo(ao longo do tempo de simulação) - h')
       ylabel('\bfSaída (y1) e setpoint (SP)')
subplot(2,1,2)
hold on
box on
stairs(t_x(1:nsim),xotimo_matriz(1:nsim),'k-','LineWidth',1.5);
hold off
    hleg1 = legend('Ações calculadas');
    set(hleg1,'Location','SouthEast')
       xlabel('\bfTempo(ao longo do tempo de predição) - h')
       ylabel('\bfAções de controle')
       
figure(3)
hold on
box on
plot(t_x(1:nsim),disturbio,'k-','LineWidth',1.5);
title('Diferença entre Cb-planta e Cb-modelo')
 hleg1 = legend('Distúrbio');
    set(hleg1,'Location','SouthEast')
       xlabel('\bfTempo(ao longo do tempo de simulação) - h')
       ylabel('\bfDistúrbio')
       
figure(4)
hold on
box on
plot(t_x(1:nsim),Cb_medido1,'k-','LineWidth',1.5);
plot(t_x(1:nsim),Cb_mod_atual1,'k:','LineWidth',1.5);
       hleg1 = legend('Cb(planta) - MPC Adaptativo','Cb(modelo) - MPC Adaptativo');
        set(hleg1,'Location','SouthEast')
        xlabel('Tempo de simulação [h]')
        ylabel('Cb [mol/l]')
       
       
% figure(4)
% hold on
% box on
% bar(t_x(2:nsim+1),ganho);
% title('Ganho entre Cb-planta e F/V')
%  hleg1 = legend('Ganho');
%     set(hleg1,'Location','SouthEast')
%        xlabel('\bfTempo(ao longo do tempo de simulação) - h')
%        ylabel('\bfGanho')