 % 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',2,...
'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
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 = [3,3,2,5,8,10];
colorList = {blue,green,cerulean,hotpink,lilac,indigo};

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


IAE_list = zeros(5*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:6
            [Kc,KI,par] = tunings(par,tauC(i),vlist(j),false);
            if (j == 2 || j == 3)
                [Kc,KI,par] = tunings(par,tauC(i),vlist(j),false,true);
            end
            tsim = tsim_list(i);
            tstep_y = tstep_y_list(i);
            ys = yslist(k);
            out = sim('Extraction');
            IAE = out.IAE(end);
            IAE_list((i-1)*6+j,k) = IAE;
            
            %plot some shite down here
            if j == 1
                fig=figure('visible','off');
            end
            
            subplot(211)
            plot(out.time,out.x(:,3),'color',colorList{j})
            hold on
            ylabel('y = x_{AE} [-]')
            subplot(212)
            plot(out.time,out.u,'color',colorList{j})
            ylabel('u = F [mol/s]')
            xlabel('time [s]')
            hold on 
            if j == 6
                subplot(211)
                plot(out.time,out.ys,'--','color',red)
                legend('v_{3}=v_{A}, I','v_{3}=v_{A}, PI','v_{2}=v_{A}, PI','v_{5}=v_{C}','v_{8}=v_{B}','v_{10}=v_{B}')
                
                name = ['/step',num2str(10*ys),'tauC',num2str(tauC(i))];
                filepath = '/Users/simenbjorvand/Documents/OneDrive - NTNU/master/ExtractionProcess/figures';  %tauc>theta,v5=p_control
                fullname = [filepath,name,'.eps'];
                saveas(fig,fullname,'epsc');
            end
            waitbar(count/(6*3*2));
            count = count+1;
        end
    end
end
close(h)