% Task  : Transformed MVs. extraction process casestudy
% compare
% Date  : 30. 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];
orange  = [255,196,31]/255; 
bluegreen = [37,235,175]/255;
burgundy = [84,20,20]/255;



%% 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
par.tau = par.ME0/par.E0;
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 the structure 
           0];       %9 parameter indicating CV transform  % 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      = stateToOutput(x0,u0,d0,par.p);

%% variables which will be varied to make up the 30 cases
tsim_list    = [200;2000;20000];                     % [s], simulation time

tstep_y_list = tsim_list*0.02;         % [s], step time for setpoint change

tstep_d = 10;                             % [s], step time for disturbances



ys = y0;
yslist = [0.5;0.1];
 
d1 = 1.0*par.S0;                                           % change in DV 1
d2 = 1.0*par.xAF0;                                         % change in DV 2

ds = [d1;d2];


tauC = [10;100;1000];
vlist = [1,9,4]+0;
colorList = {blue,green,hotpink};
namelist = {'tauC10vi1','tauC100vi1','tauC1000vi1'};

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

filepath = '/Users/simenbjorvand/Documents/OneDrive - NTNU/master/ExtractionProcess/figures/closedLoop';
IAE_list = zeros(3*3,2);
% out = cell(4,1);

h = waitbar(0,'Running simulation');
%% Simulation. Get PID tuning, and run simulation
count = 1;
for i = 1:3
    for k = 1:2
        for j = 1:3
            [Kc,KI,par] = tunings(par,tauC(i),vlist(j),false);
            tsim = tsim_list(i);
            tstep_y = tstep_y_list(i);
            ys = yslist(k);
            out = sim('Extraction');
            IAE = out.IAE(end);
            IAE_list((i-1)*3+j,k) = IAE;
            
            %plot some shite down here
            if j == 1
                fig=figure('visible','off');
                fig.Position = [100,100,450,450];
            end
            
            subplot(211)
            plot(out.time,out.x(:,3),'color',colorList{j})
            hold on
            ylabel('y = x_{AE} [-]')
            if k == 1
                ylim([-0.1,0.6])
            else
                ylim([0.05,0.25])
            end
            subplot(212)
            plot(out.time,out.u,'color',colorList{j})
            if k == 1
                ylim([-1000,4000])
            else
                ylim([-500,120])
            end
            ylabel('u = F [mol/s]')
            xlabel('time [s]')
            hold on 
            if j == 3
                step = num2str(ys*10);
                name = ['/',step,namelist{i}];
                subplot(211)
                plot(out.time,out.ys,'--','color',red)
                ldg=legend('v_{01}','v_{L1}','v_{FL1}','y_s');
                ldg.Box = 'off';   
                if k == 1
                    ldg.Location = 'east';
                end
                
                fullname = [filepath,name,'.eps'];
                saveas(fig,fullname,'epsc');
            end
            waitbar(count/(3*3*2));
            count = count+1;
        end
    end
end
close(h)