% Simulation of the new 4d model
clc
clear all
% **** System parameters ****
par = MiniLoop_parameters();
% **** System Inputs ****
z0 = 0.16; % NB! This is choke valve opening [0 1]
% System goes unstable at z= 0.15
QG_in = 4.5; %litr/min
QL_lin = 4; %litr/min
wG_in0 = par.Rog_n*QG_in/(60*1000); %kg/s
wL_in0 = par.rho_L*QL_lin/(60*1000); %kg/s
global u0;
u0 = [z0;wG_in0;wL_in0];
% ******** Finding initial states *************
disp('Initializing...')
tic
[x0,y0,par] =  MiniLoop_initialize(u0,par);
toc
disp('Initialization completed.');
disp('#########################');
%% ********* Time and logging variables ******
Minutes = 15;
tfinal=Minutes*60; % sec
t0=0;
Ts = 1;
t = t0:Ts:tfinal;
nf = length(t)-1;


u1 = u0(1)*ones(1,nf);
u2 = u0(2)*ones(1,nf);
u3 = u0(3)*ones(1,nf);
u = [u1;u2;u3];

% This a deviation form the steady state to make the model unstable
% in unstable operating point
dx = [-1e-4;0;0;0];
x = zeros(4,nf);

x(:,1) = x0+dx;
y = zeros(4,nf);
y(:,1) = MiniLoop_model(t(1),x(:,1),u(:,1),'measurements',par);

%% **********************
disp('Simulation of nonlinear model');
tic
M=eye(4);
options=odeset('Mass',M,'BDF','on','AbsTol',1e-8,'RelTol',1e-8,'MaxStep',0.1);
for k=1:nf-1
    [tt,xx]=ode15s(@MiniLoop_model,[t(k);t(k+1)],x(:,k),options,u(:,k),'derivatives',par);
    x(:,k+1) = xx(end,:)';
    if mod(t(k),tfinal/10)==0 % Show message on every 10% of simulation
        disp([num2str(floor(t(k)/tfinal*100)) '% of simulation completed'])
    end
    % *** Calculating measurements ****
    yt = MiniLoop_model(t(k+1),x(:,k+1),u(:,k+1),'measurements',par);
    y(:,k+1) = yt(1:4);
end
toc
disp('Simulation completed.');
disp('#########################');

%% ********************* Figure 1 **************************
figure(5)
clf
set(gcf,'Color',[1,1,1])
set(gca,'FontSize',10)
subplot(4,1,1)
set(gca,'FontSize',10,'OuterPosition',[0  3/4  1  1/4])
plot(t(1:nf),(x(1,1:nf)),'k');
title('a) mass of air in buffer-tank and pipe (\itx\rm_1)','FontSize',14,'FontName','Times');
ylabel('\itx\rm_1 [kg]','FontSize',12,'FontName','Times')
xlabel('time [sec]','FontSize',12,'FontName','Times')
xlim([0 tfinal])
%ylim([68 80])

subplot(4,1,2)
set(gca,'FontSize',10,'OuterPosition',[0  2/4  1  1/4])
plot(t(1:nf),(x(2,1:nf)),'k');
title('b) mass of water in pipe  (\itx\rm_2)','FontSize',14,'FontName','Times');
ylabel('\itx\rm_2 [kg]','FontSize',12,'FontName','Times')
xlabel('time [sec]','FontSize',12,'FontName','Times')
xlim([0 tfinal])
%ylim([50 62])


subplot(4,1,3)
set(gca,'FontSize',10,'OuterPosition',[0  1/4  1  1/4])
plot(t(1:nf),x(3,1:nf),'k');
title('c) mass of air in riser (\itx\rm_3)','FontSize',14,'FontName','Times');
ylabel('\itx\rm_3 [kg]','FontSize',12,'FontName','Times')
xlabel('time [sec]','FontSize',12,'FontName','Times')
xlim([0 tfinal])
%ylim([5 35])

subplot(4,1,4)
set(gca,'FontSize',10,'OuterPosition',[0  0  1  1/4])
plot(t(1:nf),x(4,1:nf),'k');
title('c) mass of water in riser (\itx\rm_4)','FontSize',14,'FontName','Times');
ylabel('\itx\rm_4 [kg]','FontSize',12,'FontName','Times')
xlabel('time [sec]','FontSize',12,'FontName','Times')
xlim([0 tfinal])
%ylim([5 35])

%% ********************* Figure 2 **************************
figure(6)
clf
set(gcf,'Color',[1,1,1])
set(gca,'FontSize',10)
subplot(3,1,1)
set(gca,'FontSize',10,'OuterPosition',[0  2/3  1  1/3])
plot(t(1:nf),y(1,1:nf),'k');
title('a) Pressure at inlet of pipeline (\itP\rm\bf_1)','FontSize',14,'FontName','Times');
ylabel('\itP\rm_1 [kpa gauge]','FontSize',12,'FontName','Times')
xlabel('time [sec]','FontSize',12,'FontName','Times')
xlim([0 tfinal])
%ylim([68 80])

subplot(3,1,2)
set(gca,'FontSize',10,'OuterPosition',[0  1/3  1  1/3])
plot(t(1:nf),y(2,1:nf),'k');
title('b) Pressure at top of riser (\itP\rm\bf_2)','FontSize',14,'FontName','Times');
ylabel('\itP\rm_2 [kpa gauge]','FontSize',12,'FontName','Times')
xlabel('time [sec]','FontSize',12,'FontName','Times')
xlim([0 tfinal])
%ylim([50 62])


subplot(3,1,3)
set(gca,'FontSize',10,'OuterPosition',[0  0  1  1/3])
plot(t(1:nf),y(3,1:nf),'k');
title('c) Volumetric flow form choke valve (\itQ_{out}\rm\bf)','FontSize',14,'FontName','Times');
ylabel('\itQ_{out}\rm [litr/min]','FontSize',12,'FontName','Times')
xlabel('time [sec]','FontSize',12,'FontName','Times')
xlim([0 tfinal])
%ylim([5 35])
