% Task  : Transformed MVs. SISO flash tank casestudy
% compare
% Date  : 24. April 2020
% Author: Simen Bjorvand
% TODO  : 
% Var   : MVs=[]; CVs=[]; States=[]  DV=[].

clc;clear;%close all
% profile on
%% Plotting options
set(0,'DefaultTextFontName','Times',...
'DefaultTextFontSize',16,...
'DefaultAxesFontName','Times',...
'DefaultAxesFontSize',16,...
'DefaultLineLineWidth',4,...
'DefaultLineMarkerSize',7.75,...
'DefaultStairLineWidth',2)
%% Colors
blue    = [0 0.4470 0.7410];
red     = [0.83 0 0.33];
green   = [0 0.5 0];
indigo  = [0.2930 0 0.5078];
lilac   = [0.7383    0.5898    0.7930];
hotpink = [0.8945    0.4102    0.5898];
cerulean = [0.1094    0.6641    0.8359];

%% Parameters
%initial values for states
par.ME0 = 1e4;
par.MR0 = 1e4; 
par.xAE0 = 3/14;
par.xWE0 = 1/14; 

%initial values for input and disturbance variables
par.F0 = 100; 

par.S0 = 100;
par.xAF0 = 0.3;

%other parameters
par.K = 1/3;
par.Kc1 = -1;
par.Kc2 = -1;
par.E0 = 140;
par.R0 = 60;

%tau, A and transformed input bias
v0 = 0;        %update v0 in some help function later


%% calculate A = df/dx (x*,u*,d*)



%% Set up simulation
par.p   = [par.K     %1
           par.Kc1   %2                                
           par.Kc2   %3
           par.E0    %4
           par.R0    %5
           par.ME0   %6
           par.MR0   %7
           0         %8 parameter indicating if an alternative CV should be used
           0];       %9 parameter indicating the structure % Parameter list
       
x0      = [par.ME0   %1
           par.MR0   %2
           par.xAE0  %3
           par.xWE0];%4                                    % initial states
       
u0      = par.F0;                                       % physical MVs bias
  
d0      = [par.S0;par.xAF0];                                     % DVs bias

y0      = par.xAE0;


%% variables which will be varied to make up the 30 cases
tsim    = 400;                                          % [s], simulation time


tstep_d = tsim*0.1;                                     % [s], step time for disturbances



 
d1 = 0.9*par.S0;                                           % change in DV 1
d2 = 1.0*par.xAF0;                                         % change in DV 2
%name = '/v0xAF03';
name = '/v0S-10';

ds = [d1;d2];


%input gain error
Ku = [1.0;1.0];


filepath = '/Users/simenbjorvand/Documents/OneDrive - NTNU/master/ExtractionProcess/figures/openLoop';
%% Simulation. Get PID tuning, and run simulations
fig=figure('visible','off');
fig.Position = [100,100,450,450];

IAE = [];
%% v1
par = tunings(par,1,true);
out = sim('ExtractionOpenLoop');

subplot(211)
plot(out.time,out.x(:,3),'DisplayName','v_{01}','color',blue)
ylim([0.2135,0.215])
ylabel('y=x_{AE} [-]')
hold on
legend('Box','off','Orientation','horizontal')
subplot(212)
plot(out.time,out.u,'color',blue)
ylabel('u=F [mol/s]')
xlabel('Time,[s]')
ylim([85,120])
hold on
IAE = [IAE;out.IAE(end)];


%% v2
par = tunings(par,2,true);
out = sim('ExtractionOpenLoop');
IAE = [IAE;out.IAE(end)];

subplot(211)
plot(out.time,out.x(:,3),'DisplayName','v_{02}','color',green)
subplot(212)
plot(out.time,out.u,'color',green)
%% v3
par = tunings(par,3,true);
out = sim('ExtractionOpenLoop');
IAE = [IAE;out.IAE(end)];

subplot(211)
plot(out.time,out.x(:,3),':','DisplayName','v_{03}','color',indigo)
plot(out.time,out.ys,'--','DisplayName','y_s','color',red)
subplot(212)
plot(out.time,out.u,':','color',indigo)

 
fullname = [filepath,name,'.eps'];
saveas(fig,fullname,'epsc');

% fig=figure('visible','on');
% %% v4
% par = tunings(par,4,true);
% out = sim('ExtractionOpenLoop');
% 
% subplot(211)
% plot(out.time,out.x(:,3),'DisplayName','v_{FL1}')
% legend
% hold on
% subplot(212)
% plot(out.time,out.u)
% hold on
% %% v5
% par = tunings(par,5,true);
% out = sim('ExtractionOpenLoop');
% 
% subplot(211)
% plot(out.time,out.x(:,3),'DisplayName','v_{FL2}')
% subplot(212)
% plot(out.time,out.u)
% 
% %% v6
% par = tunings(par,6,true);
% out = sim('ExtractionOpenLoop');
% 
% subplot(211)
% plot(out.time,out.x(:,3),'DisplayName','v_{FL3}')
% subplot(212)
% plot(out.time,out.u)


% fig=figure('visible','on');
% %% v7
% par = tunings(par,7,true);
% out = sim('ExtractionOpenLoop');
% 
% subplot(211)
% plot(out.time,out.x(:,3),'DisplayName','v7')
% legend
% hold on
% subplot(212)
% plot(out.time,out.u)
% hold on
% %% v8
% par = tunings(par,8,true);
% out = sim('ExtractionOpenLoop');
% 
% subplot(211)
% plot(out.time,out.x(:,3),'DisplayName','v8')
% subplot(212)
% plot(out.time,out.u)
% %% v9
% par = tunings(par,9,true);
% out = sim('ExtractionOpenLoop');
% 
% subplot(211)
% plot(out.time,out.x(:,3),'DisplayName','v9')
% subplot(212)
% plot(out.time,out.u)
% %% v10
% par = tunings(par,10,true);
% out = sim('ExtractionOpenLoop');
% 
% subplot(211)
% plot(out.time,out.x(:,3),'DisplayName','v10')
% subplot(212)
% plot(out.time,out.u)