%KK4sim.m
clear all

global W V Y YY tt U

%load results/NyestKKopt/opt010
%load results/NyestKKopt/sim01

load opt010
load sim01

format compact
format long g

%Pådrag 
    %u = [F_CH4; F_LUFT; w_8; w_15; w_12; w_21]
    %u = [F_CH4; w_8; w_15]
    u = U(1:6); %Henter verdier fra optimalisering ved nom.verdier for forstyrrelsene
    u(5) = 1000;    %kg/s, "åpen" ventil er aktiv beskrankning
    u(6) = 0.1;       %kg/s, stengt ventil er aktiv beskrankning
    u = u(:);
    u_0 = u;            %Benyttes som bias i regulatorene    
   
    %u(1) = 0.9.*u(1);   %Sprang i pådrag, 10%-økning
       
%Forstyrrelser
    %d = [p_CH4; p_LUFT; T_CH4; T_LUFT]
    d = [20E5; 1E5; 283; 283];  %Pa (d(1:2)) og K (d(3:4)) (nom. verdier hentet fra Skogestad(2000))
            
%Initialbetingelser(verdier funnet vha fmincon)
    %X = [X_BK_0; T_OH_0; T_FD_0; h_DR_0; T_K_0; h_KT_0; T_DE_0; T_EK_0; T_FV_0] 
    X = U(7:40); %Henter verdier fra optimalisering ved nom.forstyrrelser
    X(16) = 1.0;    %Setter h_DR = 1.0 m      
    X(21) = 1.0;    %Setter h_KT = 1.0 m
    X = X(:);       %Setter X på kolonneform
    
%Antall celler, N, i varmevekslerne
N_OH = 3;       %Antall celler i overheteren
N_FD = 3;       %Antall celler i fordamperen
N_K = 3;        %Antall celler i kondensatoren
N_EK = 3;       %Antall celler i ekonomiseren
N_FV = 3;       %Antall celler i forvarmeren

%Konstante tilstandsverdier som benyttes til å beregne hold-up i KK2.m 
T_OH_RG_0 = uX010(12:14);   %K, Temp(overheter,røykgass,celle 1, 2 og 3)                                       
T_OH_D_0 = uX010(15:17);    %K, Temp(overheter,damp,celle 1, 2 og 3) 
T_FD_RG_0 = uX010(18:20);   %K, Temp(fordamper, røykgass, celle 1, 2 og 3)
T_EK_RG_0 = uX010(29:31);   %K, Temp(ekonomiser, røykgass, celle 1, 2 og 3)
T_FV_RG_0 = uX010(35:37);   %K, Temp(forvarmer, røykgass, celle 1, 2 og 3)

%Settpunkt
T_BKs = 1773;   %K
h_KTs = 1.0;    %m
h_DRs = 1.0;    %m
T_DEs = 374;    %K
w_21s = 0.1;    %kg/s

if 1==1
%Integrert reguleringsavvik legges inn som tilstander
e_T_BK_int = T_BKs - X(5);  %Integrert reguleringsavvik 
e_h_DR_int = h_DRs - X(16); %Integrert reguleringsavvik 
e_h_KT_int = h_KTs - X(21); %Integrert reguleringsavvik 
e_T_DE_int = T_DEs - X(22); %Integrert reguleringsavvik 
e_w_21_int = w_21s - u(6);      %Integrert reguleringsavvik


X = [X; e_T_BK_int; e_h_DR_int; e_h_KT_int; e_T_DE_int; e_w_21_int];
end

% ODE15s
    t = 3E4;      %s
    %t = 100;
    [tid,X] = ode15s('KK4',[0 t],X,[],u,d,T_BKs,h_DRs,h_KTs,T_DEs,w_21s,u_0,N_OH,N_FD,N_K,N_EK,N_FV,T_OH_RG_0,T_OH_D_0,T_FD_RG_0,T_EK_RG_0,T_FV_RG_0)


V
%save results/NyestKKopt/sim11 U V
% Figurer
if 1==0
figure(11)
    subplot(321)
    plot(tid, X(:,1)), ylabel('mol O2'), title('Hold-up BK');
    subplot(322)
    plot(tid, X(:,2)), ylabel('mol N2'), title('Hold-up BK');
    subplot(323)
    plot(tid, X(:,3)), ylabel('mol CO2'), title('Hold-up BK');
    subplot(324)
    plot(tid, X(:,4)), ylabel('mol H2O'), title('Hold-up BK');
    subplot(325)
    plot(tid, X(:,5)), xlabel('Tid [s]'), ylabel('Kelvin'), title('Temperatur BK');
end

if 1==1
tid1 = tt(1);
yy = YY(1,:);
    for i = [2:length(tt)];
        if tt(i-1) < tt(i);
        tid1 = [tid1; tt(i)];
        yy = [yy; YY(i,:)];   
        end
    end
end
if 1==0
figure(1)
    subplot(311)
    plot(tid1, yy(:,5)), ylabel('w(21), [kg/s]'), title('Settpunktsendring for w(21) fra 0.10 til 0.15 m')%, xlabel('Tid [s]');
    subplot(312)
    plot(tid, X(:,22)), ylabel('Temperatur i DE, [K]')
    subplot(313)
    plot(tid1, yy(:,1)), ylabel('F(CH4), [mol/s]'), xlabel('Tid [s]')%, title('F(CH4)');
figure(2)
    subplot(211)
    plot(tid, X(:,5)), ylabel('Temperatur BK, [K]'), title('+ 10%-sprang i F(LUFT)');
    subplot(212)
    plot(tid1, yy(:,2)), ylabel('F(LUFT), [mol/s]'), xlabel('Tid [s]');%, title('F(LUFT)')
figure(3)
    subplot(211)
    plot(tid, X(:,16)), ylabel('Væskenivå i DR, [m]')%, title('Settpunktsendring for h(DR) fra 1.0m til 1.2m');
    subplot(212)
    plot(tid1, yy(:,4)), ylabel('w(15), [kg/s]'), xlabel('Tid [s]');
figure(4)
    subplot(211)
    plot(tid, X(:,21)), ylabel('Væskenivå i KT, [m]')%, title('Settpunktsendring for h(KT) fra 1.0m til 1.2m');
    subplot(212)
    plot(tid1, yy(:,3)), ylabel('w(8), [kg/s]'), xlabel('Tid [s]');
end
if 1==0
figure(5)
    subplot(121)
    plot(tid, X(:,5)), ylabel('Temperatur BK, [K]'), title('Forstyrrelse P(CH4) = 16 bara')%, title('+ 10%-sprang i F(CH4)');
    subplot(122)
    plot(tid1, yy(:,2)), ylabel('F(LUFT), [mol/s]') %, title('')%, xlabel('Tid [s]');
    
figure(6)
    subplot(221)
    plot(tid, X(:,16)), ylabel('Væskenivå i DR, [m]') %title('+ 10%-sprang i w(15)');
    subplot(222)
    plot(tid1, yy(:,4)), ylabel('w(15), [kg/s]');
    subplot(223)
    plot(tid, X(:,21)), ylabel('Væskenivå i KT, [m]'), xlabel('Tid [s]') %title('+ 10%-sprang i w(8)');
    subplot(224)
    plot(tid1, yy(:,3)), ylabel('w(8), [kg/s]'), xlabel('Tid [s]');

figure(7)
    subplot(211)
    plot(tid, X(:,22)), ylabel('Temperatur i DE, [K]')%, title('+ 10%-sprang i w(21)')%, xlabel('Tid [s]');
    subplot(212)
    plot(tid1, yy(:,5)), ylabel('w(21), [kg/s]'), xlabel('Tid [s]');
end
figure(8)
    subplot(521)
    plot(tid1, yy(:,1)), ylabel('F(CH4), [mol/s]'), title('Forstyrrelse T(LUFT) = 288 K')
    subplot(522)
    plot(tid, X(:,22)), ylabel('T(DE), [K]')
    subplot(523)
    plot(tid1, yy(:,2)), ylabel('F(LUFT), [mol/s]')
    subplot(524)
    plot(tid, X(:,5)), ylabel('T(BK), [K]')
    subplot(525)
    plot(tid1, yy(:,3)), ylabel('w(8), [kg/s]')
    subplot(526)
    plot(tid, X(:,21)), ylabel('h(KT), [m]')
    subplot(527)
    plot(tid1, yy(:,4)), ylabel('w(15), [kg/s]')
    subplot(528)
    plot(tid, X(:,16)), ylabel('h(DR), [m]'), xlabel('Tid [s]')
    subplot(529)
    plot(tid1, yy(:,5)), ylabel('w(21), [kg/s]'), xlabel('Tid [s]')
    
figure(9)
    plot(tid1, yy(:,6)), ylabel('J, [NOK/h]'), xlabel('Tid [s]'), title('Forstyrrelse T(LUFT) = 288 K')