% file to run cascade implementation for transformed inputs mixing example

SetUpSimulations % initialize values
SetUpPlots

tv1 =100; tv2=100;
q1f = q10; q2f = q20; % bias for input u
tauTs = 1;            % tracking time constant. not used

taucs = 0.01; % closed loop time constant for slave controller

%% Master Flow controller tuning
taucq  = 5*taucs;                          % [s]. closed loop time constant
KIqM   = 1/1.1/(taucq);         % integral gain tune on initial gain = 1.1
vqbias = q0;                                        %[kg/s] controller bias


pairing = 2;
for i = [0,1]
     transf = i; % 0-static1; 1-dynamic
     
     
     
     switch transf
          case 0 % static
               kT     = q0/V;
               tauT   = V/q0;
               vTbias = T0;
               n = 1; % how much faster than the residence time
               switch pairing
                    case 1                                        % u1-v1; u2-v2
                         % slave controller tunings
                         ku1 = 1;                           % gain from u1 to v1
                         ku2 = (T20-T0)/(q10+q20);           % gain from u2 to v2
                         Kparu1 = [1 0]; Kparu2 = [0 1]; %select pairing in Simulink
                    case 2                                        % u1-v2; u2-v1
                         ku1 = (T10-T0)/(q10+q20);          % gain from u1 to v2
                         ku2 = 1;                           % gain from u2 to v1
                         Kparu1 = [0 1]; Kparu2 = [1 0]; % select pairing in Simulink
                    otherwise
                         disp('incorrect pairing')
               end
          case 1 % dynamic
               kT     = -B/A;  % gain of transformed system for T
               tauT   = -1/A;     % time constant of the transformed system
               vTbias = -A/B*T0;                   % controller bias for TC
               n = 2; % how much faster than the residence time
               switch pairing
                    case 1                              % u1-v1; u2-v2
                         ku1 = 1;                       % gain from u1 to v1
                         ku2 = 1/V*1/B*(T20-T0);        % initial gain from u2 to v2
                         Kparu1 = [1 0]; Kparu2 = [0 1];%select pairing in Simulink
                    case 2                              % u1-v2; u2-v1
                         ku1 = 1/V*1/B*(T10-T0);        % initial gain from u1 to v2
                         ku2 = 1;                       % gain from u2 to v1
                         Kparu1 = [0 1]; Kparu2 = [1 0];% select pairing in Simulink
                    otherwise
                         disp('incorrect pairing')
               end
          otherwise
               disp('Incorect transformation')
     end
     
     % Slave controller tunings
     Ku1 = 1/ku1*1/(taucs);     % integral slave controller to find u1
     Ku2 = 1/ku2*1/(taucs);     % integral slave controller to find u2
     
     % Master temperature controller tuning. independent of pairing
     taucT  = V/q0/n;         % [s]. closed loop time constant for  TC
     KcTM    = 1/kT*tauT/(taucT);           % proportional gain for TC
     tauITM  = min(-1/A, 4*(taucT));            % integral time for TC
     
     
     p   = [V;A;B; transf];                     % vector of parameters
     
     if OL == 1 % simulation without outter controller
          vqbias_step =  qstep;
          vTbias_step =  Tstep;
     else % simulation with outer controller
          vqbias_step = vqbias;
          vTbias_step = vTbias;
     end
     
     sim('simTransformedMVs_cascade')
     
     %% Plot results
     figure(1)
     plot(time,d(:,1),'Color','blue')
     
     
     figure(2)
     plot(time,d(:,2),'Color','red')
     
     
     figure(3)
     plot(time, q,'Color',clr(i+1,:),'LineStyle',lin{i+1})
     hold on
     
     
     figure(4)
     plot(time, T,'Color',clr(i+1,:),'LineStyle',lin{i+1})
     hold on
     
     
     figure(5)
     plot(time,v(:,1),'Color',cherry)
     
     
     
     figure(6)
     plot(time,v(:,2),'Color',cherry)
     
     
     
     figure(7)
     plot(time,uplant(:,1),'Color',clr(i+1,:),'LineStyle',lin{i+1})
     hold on
     ylim([1.9 10.1])
     
     
     figure(8)
     plot(time,uplant(:,2),'Color',clr(i+1,:),'LineStyle',lin{i+1})
     hold on
     
     
     
end

figure(1);
ylim([19 31])
legend('d_1=T_1', 'Location','best')
ylabel('Temperature T_1 [^oC]')

figure(2)
ylim([49 61])
legend('d_2=T_2', 'Location','best')
ylabel('Temperature T_2 [^oC]')

figure(3)
ylim([9 15])
plot(time,qsp,'--','Color',cherry)
legend('y_1 = q (static)','y_1 = q (dynamic)','y_1^s = q^s', 'Location','SouthEast')
ylabel('Flowratw q [kg s^{-1}]')

figure(4)
plot(time, Tsp,':','Color',cherry)
legend('y_2=T (static)','y_2=T (dynamic)' ,'y_2^s = T^s', 'Location','best')
ylim([34 46])
ylabel('Temperature T [^oC]')

figure(5)
ylim([9 15])
legend('v_{01}^s','Location','northeast')
ylabel('Tranformed input')

figure(6)
ylim([34 55])
legend('v_{2}^s','Location','northeast')
ylabel('Tranformed input ')

figure(7)
ylim([1 12])
legend('u_1 = q_1 (static)','u_1 = q_1 (dynamic)','Location','northeast')
ylabel('Flowrate q_1 [kg s^{-1}]')

figure(8)
ylim([1 12])
legend('u_2 = q_2 (static)','u_2 = q_2 (dynamic)','Location','southeast')
xlabel('Time')
ylabel('Flowrate q_2 [kg s^{-1}]')




var = {'d1','d2','y1','y2','v1','v2','u1','u2'};
for i = 1:8
     figure(i)
     filename = strcat('mixTc_u1-v2_',var{i},'.eps');
     legend boxoff
     set(gcf,'Position',pos)
     print('-depsc2','-r600',filename)
end