close all
clc
clear all
%Try to optimize distillation column "A" for different disturbances!

load('x0_F15_p007.mat') %Initial values at F=1.5 and pV=0.07

Jcol
D0  = D;
B0  = B;
xD0 = xD;
xB0 = xB;
x0  = x; 
LT0 = u_opt(1);
VB0 = u_opt(2);
F0  = F;
pV0 = pV;
% Define nominal inputs!

% Nominal inputs from optimization file
% LT0=2.70629;                          % Reflux
% VB0= 3.20629;                          % Boilup
%  D0=0.5;                               % Distillate
%  B0=0.5;                               % Bottoms
% F0=1.0;                               % Feedrate
% x0 = 0.5*ones(1,82)' ; %This is X in cola_lv (this one is given in col_constr)

%Constants
zF=0.5;                              % Feed composition
qF=1.0;                              % Feed liquid fraction
Uinit = [LT0 VB0 D0 B0 F0 zF qF]';


%global NT d
NT = 41 ;

% Inputs
u0 = Uinit(1:2) ; %L, V

%Prices
pf = 1 ;  pB = 2 ; pD = 1 ;

%% Setpoints and min values

xBsp    = 0.95;
xDmin   = 0.95; 

%% Tuning
TuningCol


%% Matrices to extract values in simulation
%xprime is of size 41*2+1=83
%Stage 1 is the reboiler
Extract_J       =zeros(83,1)';
Extract_J(end)  =1;

Extract_x       =ones(83,1)';
Extract_x(end)  =0;

%Now we have 82 states. 
Extract_xB      =zeros(82,1)'; %horizontal
Extract_xB(1)   =1;

Extract_xD      =zeros(82,1)'; %horizontal
Extract_xD(NT)  = 1; 



%% Steps for finding open loop tf
LTstep=  0*LT0*.2;
LTstep_t = 10; %time

VBstep=  0*VB0*.1;
VBstep_t = 10; %time

F0step=  1*F0*.1;
F0step_t = 10; %time

%% Steps for simulation

% F0step=  0.02;
% F0step_t = 10; %time
% 
% F0step2=  0*(0.15-F0step); % To get to 1.65 in RII
% F0step2_t = 30; %time
% 
 pVstep=0*(0.13-pV0);
 pVstept= 50;
%% Simulate
tend=150;
a = sim('OneColA','SimulationMode','normal'); 
%%
xmin=0.95*ones(length(a.time));
Vmax_plot=4*ones(length(a.time));
Lmax_plot=10*ones(length(a.time));

figure
subplot(4,1,1)
plot(a.time,a.xB)
hold on
%plot(a.t,a.xB)
plot(a.time, xmin,'--','Color',[  0.9290    0.6940    0.1250])
%xlabel('Time (s)')
ylabel('x_B, (mol/mol)')
%ylim([0.9 1])
%legend('x_B')
%legend boxoff
ytickformat('%.2f')

subplot(4,1,2)
plot(a.time,a.xD)
hold on
%plot(a.time, a.xDsp,'k--')
plot(a.time, xmin,'--','Color',[  0.9290    0.6940    0.1250])
%xlabel('Time (s)')
ylabel('x_D, (mol/mol)')
%ylim([0.9 1])
%legend('x_D')
%legend boxoff
ytickformat('%.2f')

subplot(4,1,3)
plot(a.time,a.VB)
hold on
plot(a.time, Vmax_plot,'--','Color',[  0.9290    0.6940    0.1250])
%xlabel('Time (s)')
ylabel('Boil-up, V (mol/s)')
%ylim([0.9 1])
%legend('V', 'V^{max}')
%legend boxoff
%ytickformat('%.2f')

subplot(4,1,4)
plot(a.time,a.LT)
hold on
%plot(a.time, Lmax_plot,'--','Color',[  0.9290    0.6940    0.1250])
%xlabel('Time (s)')
ylabel('Reflux, L (mol/s)')
%ylim([0.9 1])
%legend('L', 'L^{max}')
%legend boxoff

%%
figure
subplot(3,1,1)
plot(a.time, a.Fin)
ylabel('F (mol/s)')

subplot(3,1,2)
plot(a.time, a.pV)
ylabel('p_V ($/mol)')

subplot(3,1,3)
plot(a.time,a.Jactual)
xlabel('Time (s)')
ylabel('J ($/s)')



%% Save results
%filename = 'ColSim.mat';
%save(filename)
% x = load('ColSim.mat');
% 
% names = fieldnames(x);
% for iname = 1:length(names)
%   x.(['ColSim', names{iname}]) = x.(names{iname});
%   x = rmfield(x, names{iname});
% end
% %save('newmymat.mat', '-struct', 's');
% save('ColSim.mat', '-struct', 'x');
