%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%_________________Starter of the Distillation column_______________________
%<<<<                                                                 >>>>%
%  This file can run both with open loop, to tune the controller based on %
%  which operation point is chosen, but also with closed loop to simulate %
%  how controller handles disturbances in feed rate and feed composition. %
%  To choose between open and closed loop open simulink file              %
%  "colas_nonlin_operation_All" and change the position of the            %
%  "Manual Switch" by dobbelklicking on it. Then run this file            %
%  "Starter_All".
%<<<<                                                                 >>>>%


clc
clear all
clf
N=41;   %Number of stages
i=1:N; n(i)=i; % Matrix with stages in increasing order
Tb1=261.45; %Boiling point of iso-butane
Tb2=272.65; %Boiling point of n-butane
deltaTb= Tb1-Tb2; %Boiling point difference

%Initial delay in boiler
Timedelay=2;


%TC
TC=15; %Location of TC, numbering from the bottum of the column
T_under=TC-1; %Stages under TC
T_control=1;
T_over=N-TC; %Stages over TC 

%Initial values 
tot_time=5;
Steptime=0; Stepsize=0; %Reboiler
Steptime_setpunkt=0; Stepsize_setpunkt=0; %Setpoint in TC
dF_time=0; dF_size=0; dzF_time=0; dzF_size=0; % Disturbances

global OP F0 V0 L0

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   Choose operation point                                %
OP=2;  %<<-------                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if OP==1
    %Operation point 1
    load cola_init_1
    load Parameters
    P=P(1); I=I(1); T_s=T_s1(1);
elseif OP==2
    %Operation point 2
    load cola_init_2
    load Parameters
    P=P(2); I=I(2); T_s=T_s1(2);
elseif OP==3
    %Operation point 3
    load cola_init_3
    load Parameters
    P=P(3); I=I(3); T_s=T_s1(3);
elseif OP==4
    %Operation point 4a
    load cola_init_4a
    load Parameters
    P=P(4); I=I(4); T_s=T_s1(4);
elseif OP==5
    %Operation point 4b
    load cola_init_4b
    load Parameters
    P=P(5); I=I(5); T_s=T_s1(5);
elseif OP==6
    %Operation point 4c
    load cola_init_4c
    load Parameters
    P=P(6); I=I(6); T_s=T_s1(6);
elseif OP==7
    %Operation point 5
    load cola_init_5
    load Parameters
    P=P(7); I=I(7); T_s=T_s1(7);
elseif OP==8
    %Operation point 6
    load cola_init_6
    load Parameters
    P=P(8); I=I(8); T_s=T_s1(8);
elseif OP==9
    %Operation point 7
    load cola_init_7
    load Parameters
    P=P(9); I=I(9); T_s=T_s1(9);
elseif OP==10
    %Operation point 8
    load cola_init_8
    load Parameters
    P=P(10); I=I(10); T_s=T_s1(10);
else 
    disp('Error in determing the operation point')
end

F0=Uinit(5);
V0=Uinit(2);
L0=Uinit(1);
% Runs simulation to fin if the loop is closed or opened
sim('colas_nonlin_operation_All')
Timedelay= Timedelay*3.20629/Setpoint_V(1);
%%


if any(RUN)>=1
%------------------- Runing with closed loop ------------------------------
    Steptime_setpunkt=10; Stepsize_setpunkt=0.0;
    dF_time=10; dF_size=0.1; dzF_time=150; dzF_size=-0.05;
    tot_time=300;
    sim('colas_nonlin_operation_All')
    Y=T_on15th; % To be analysed
    y1_1=y1;
    y2_1=y2;
    V1=V;
    t1=t;
    
    %Plot of Temperatureprofiles
    figure (1)
    subplot(1,2,1)
    T_grader=T-273.25;
    figure(1)
    plot(n,T_grader(1,:),n,T_grader(end,:),TC,T_grader(end,TC), 'o')
    legend('Temperatureprofile before disturbance',...
        'Temperatureprofile after disturbance',...
        'Location of temperature controller')
    title('Temperatureprofiles','BackgroundColor',[.99 .99 .99])
    xmin=n(1);
    xmax=n(end);
    ymin=min(T_grader(1,end),T_grader(end,end));
    ymax=max(T_grader(1,1),T_grader(end,1));    
    axis([xmin xmax ymin ymax])
    grid on
    xlabel('Stage number')
    ylabel('Temperature [oC]')

    %Plot of steprespons
    subplot(1,2,2)
    plot(t, Setpoint_T, t, Y)
    legend('Step in setpoint','Step respons')
    title('Behaviour of the temperature contoller with change in setpoint',...
        'BackgroundColor',[.99 .99 .99])
    grid on
    xlabel('Time [min]')
    ylabel('Temperature controler [oC]')


    %Simulation with tuning done in Operation point 1
    clear P I
    load Parameters 
    P=P(1); I=I(1);
    sim('colas_nonlin_operation_All')
    Y2=T_on15th;
    y1_2=y1;
    y2_2=y2;
    V2=V;
    t2=t;

    %Simulation with Global tuning (tuning from operation point 4a)
    clear P I
    load Parameters 
    P=P(4); I=I(4);
    sim('colas_nonlin_operation_All')
    Y3=T_on15th;
    y1_3=y1;
    y2_3=y2;
    V3=V;
    t3=t;


%--------------------- PLOTS ----------------------------------
%------------- Temperature controller -------------------------
    if OP==1
        figure (2)
        plot(t3, Setpoint_T, '--','color', [0.5 0.5 0.5])
        hold on
        plot(t2, Y2, 'color', [0.9 0 0])
        hold on
        plot(t1, Y, '--', 'color', [0 0 0.9])
        hold on
        plot(t3, Y3, 'color', [0 0.6 0])
        grid on
        axis([0 tot_time 267 272.5])
        xlabel('Time [min]')
        ylabel('Temperature controller [^oC]')
        legend('Setpoint','Controller tuned at operation point 1',...
            'Controller retuned at current operation point^{ }_{ }_{ }',...
            'Controller with recommended global tuning',...
            'Location', 'Best')

    %------------ Impurity in top product ---------------------------
        figure (3)
        plot(t2, y1_2, 'color', [0.9 0 0])
        hold on
        plot(t1, y1_1, '--', 'color', [0 0 0.9])
        hold on
        plot(t3, y1_3, 'color', [0 0.6 0])
        grid on
        axis([0 tot_time 0 0.035])
        xlabel('Time [min]')
        ylabel('Impurity in top product')
        legend('Controller tuned at operation point 1',...
            'Controller retuned at current operation point^{ }_{ }_{ }',...
            'Controller with recommended global tuning',...
            'Location', 'Best')

    %---------- Impurity in bottom product ---------------------------
        figure (4)
        plot(t2, y2_2, 'color', [0.9 0 0])
        hold on
        plot(t1, y2_1, '--', 'color', [0 0 0.9])
        hold on
        plot(t3, y2_3, 'color', [0 0.6 0])
        grid on
        axis([0 tot_time 0 0.035])
        xlabel('Time [min]')
        ylabel('Impurity in bottom product')
        legend('Controller tuned at operation point 1',...
            'Controller retuned at current operation point^{ }_{ }_{ }',...
            'Controller with recommended global tuning',...
            'Location', 'Best')

      %-------------- Reboiler duty ----------------------------------  
        figure (5)
        plot(t2, V2, 'color', [0.9 0 0])
        hold on
        plot(t1, V1, '--', 'color', [0 0 0.9])
        hold on
        plot(t3, V3, 'color', [0 0.6 0])
        grid on
        axis([0 tot_time 1.5 3.5])
        xlabel('Time [min]')
        ylabel('Reboiler duty')
        legend('Controller tuned at operation point 1',...
            'Controller retuned at current operation point^{ }_{ }_{ }',...
            'Controller with recommended global tuning',...
            'Location', 'Best')  
        
    elseif OP==4
        figure (2)
        plot(t3, Setpoint_T, '--','color', [0.5 0.5 0.5])
        hold on
        plot(t2, Y2, 'color', [0.9 0 0])
        hold on
        plot(t1, Y, 'color', [0 0 0.9])
        hold on
        plot(t3, Y3, '--', 'color', [0 0.6 0])
        grid on
        axis([0 tot_time 267 272.5])
        xlabel('Time [min]')
        ylabel('Temperature controller [^oC]')
        legend('Setpoint','Controller tuned at operation point 1',...
            'Controller retuned at current operation point^{ }_{ }_{ }',...
            'Controller with recommended global tuning',...
            'Location', 'Best')

    %------------ Impurity in top product ---------------------------
        figure (3)
        plot(t2, y1_2, 'color', [0.9 0 0])
        hold on
        plot(t1, y1_1, 'color', [0 0 0.9])
        hold on
        plot(t3, y1_3, '--', 'color', [0 0.6 0])
        grid on
        axis([0 tot_time 0 0.035])
        xlabel('Time [min]')
        ylabel('Impurity in top product')
        legend('Controller tuned at operation point 1',...
            'Controller retuned at current operation point^{ }_{ }_{ }',...
            'Controller with recommended global tuning',...
            'Location', 'Best')

    %---------- Impurity in bottom product ---------------------------
        figure (4)
        plot(t2, y2_2, 'color', [0.9 0 0])
        hold on
        plot(t1, y2_1, 'color', [0 0 0.9])
        hold on
        plot(t3, y2_3, '--', 'color', [0 0.6 0])
        grid on
        axis([0 tot_time 0 0.035])
        xlabel('Time [min]')
        ylabel('Impurity in bottom product')
        legend('Controller tuned at operation point 1',...
            'Controller retuned at current operation point^{ }_{ }_{ }',...
            'Controller with recommended global tuning',...
            'Location', 'Best')

      %-------------- Reboiler duty ----------------------------------  
        figure (5)
        plot(t2, V2, 'color', [0.9 0 0])
        hold on
        plot(t1, V1, 'color', [0 0 0.9])
        hold on
        plot(t3, V3, '--', 'color', [0 0.6 0])
        grid on
        axis([0 tot_time 1.5 3.5])
        xlabel('Time [min]')
        ylabel('Reboiler duty')
        legend('Controller tuned at operation point 1',...
            'Controller retuned at current operation point^{ }_{ }_{ }',...
            'Controller with recommended global tuning',...
            'Location', 'Best')
                
    else   
        figure (2)
        plot(t3, Setpoint_T, '--','color', [0.5 0.5 0.5])
        hold on
        plot(t2, Y2, 'color', [0.9 0 0])
        hold on
        plot(t1, Y, 'color', [0 0 0.9])
        hold on
        plot(t3, Y3, 'color', [0 0.6 0])
        grid on
        axis([0 tot_time 267 272.5])
        xlabel('Time [min]')
        ylabel('Temperature controller [^oC]')
        legend('Setpoint','Controller tuned at operation point 1',...
            'Controller retuned at current operation point^{ }_{ }_{ }',...
            'Controller with recommended global tuning',...
            'Location', 'Best')

    %------------ Impurity in top product ---------------------------
        figure (3)
        plot(t2, y1_2, 'color', [0.9 0 0])
        hold on
        plot(t1, y1_1, 'color', [0 0 0.9])
        hold on
        plot(t3, y1_3, 'color', [0 0.6 0])
        grid on
        axis([0 tot_time 0 0.035])
        xlabel('Time [min]')
        ylabel('Impurity in top product')
        legend('Controller tuned at operation point 1',...
            'Controller retuned at current operation point^{ }_{ }_{ }',...
            'Controller with recommended global tuning',...
            'Location', 'Best')

    %---------- Impurity in bottom product ---------------------------
        figure (4)
        plot(t2, y2_2, 'color', [0.9 0 0])
        hold on
        plot(t1, y2_1, 'color', [0 0 0.9])
        hold on
        plot(t3, y2_3, 'color', [0 0.6 0])
        grid on
        axis([0 tot_time 0 0.035])
        xlabel('Time [min]')
        ylabel('Impurity in bottom product')
        legend('Controller tuned at operation point 1',...
            'Controller retuned at current operation point^{ }_{ }_{ }',...
            'Controller with recommended global tuning',...
            'Location', 'Best')

      %-------------- Reboiler duty ----------------------------------  
        figure (5)
        plot(t2, V2, 'color', [0.9 0 0])
        hold on
        plot(t1, V1, 'color', [0 0 0.9])
        hold on
        plot(t3, V3, 'color', [0 0.6 0])
        grid on
        axis([0 tot_time 1.5 3.5])
        xlabel('Time [min]')
        ylabel('Reboiler duty')
        legend('Controller tuned at operation point 1',...
            'Controller retuned at current operation point^{ }_{ }_{ }',...
            'Controller with recommended global tuning',...
            'Location', 'Best')
    end
else
%------------------- Runing with open loop --------------------------------
%-------------------  For tuning of TC  -----------------------------------
    Steptime=15; Stepsize=0.02;
    tot_time=500;
    sim('colas_nonlin_operation_All')
    Y=T_on15th; % To be analysed
    
    
    %Plot of Temperatureprofiles
    figure (1)
    subplot(2,2,[1 3])
    T_grader=T-273.25;
    figure(1)
    plot(n,T_grader(1,:),n,T_grader(end,:),TC,T_grader(end,TC), 'o')
    legend('Temperatureprofile before disturbance',...
        'Temperatureprofile after disturbance',...
        'Location of temperature controller')
    title('Temperatureprofiles','BackgroundColor',[.99 .99 .99])
    xmin=n(1);
    xmax=n(end);
    ymin=min(T_grader(1,end),T_grader(end,end));
    ymax=max(T_grader(1,1),T_grader(end,1));    
    axis([xmin xmax ymin ymax])
    grid on
    xlabel('Stage number')
    ylabel('Temperature [oC]')

        
    %SIMC tuning
    i=1:length(t); t_simc(i)=i;
    t_start=round(interp1(t,t_simc,Steptime+Timedelay,'linear'));     %or 'nearest'
    k63=Y(t_start)+(.63*(Y(end)-Y(t_start)));
   
    delay=Timedelay;
    t_action= t(t_start+delay:end);
    Y_action=Y(t_start+delay:end);

    %Finding tau_c
    tau_c=1*delay;
    
    
    %Finding timekonstant
    [diff63 I]=min(abs(Y_action-k63));
    tau1_plot=t_action(I);
    %tau12=t(I)
    tau1=tau1_plot-delay-Steptime;


    subplot (2,2,2)
    plot(t, Y, tau1_plot, k63, 'o')
    hold on
    [AX,H1,H2]= plotyy(tau1_plot, k63, t, Setpoint_V);
    set(get(AX(1),'Ylabel'),'String','Stepresponse in ouput');
    set(get(AX(2),'Ylabel'),'String','Step in input');
    grid on
    legend('Step respons', '\tau_1')
    title('Step respons in V-T controller','BackgroundColor',[.99 .99 .99])

    %Finding K_c
    Range_Y=1;
    Range_u=1; 
    k=((Y(end)-Y(1))/Range_Y)/(Stepsize/Range_u);
    K_c=tau1/(k*(tau_c+delay));

    %Finding tau_I
    tau_I=min(tau1, 4*(tau_c+delay));

    subplot (2,2,4)
    text(0.05,0.7, ['SIMC tuning results', char(10)...
        char(10)...
        'delay=',num2str(delay), char(10)...
        '\tau_c=',num2str(tau_c), char(10)...
        '\tau_1=',num2str(tau1), char(10)...
        'k=',num2str(k), char(10)...
        'K_c=' ,num2str(K_c), char(10)...
        '\tau_I=',num2str(tau_I), char(10)...
        'K_c/\tau_I=',num2str(K_c/tau_I)],...
        'BackgroundColor', 'white', 'EdgeColor',[ 0 0 0])
    set(gca,'Visible','off')


    %Assume pure integrator
    k_merket= (Y_action(150)-Y_action(1))/(Stepsize*(t_action(150)-t_action(1)));
    K_c2= 1/(k_merket*(tau_c+delay));
    tau_I2= 4*(tau_c+delay);
    subplot (2,2,4)
    text(0.7,0.7, ['SIMC tuning results 2', char(10)...
        char(10)...
        'delay=',num2str(delay), char(10)...
        '\tau_c=',num2str(tau_c), char(10)...
        'k`=',num2str(k_merket), char(10)...
        'K_c=' ,num2str(K_c2), char(10)...
        '\tau_I=',num2str(tau_I2), char(10)...
        'K_c/\tau_I=',num2str(K_c2/tau_I2)],...
        'BackgroundColor', 'white', 'EdgeColor',[ 0 0 0])
    set(gca,'Visible','off')
    
end

