clc
clear all
% **** System Inputs ****
par = v1_gaslift_parameters();  % Parameter Constants for Model
% **** Automatic ********
zs=0.1:0.05:1;
ns = length(zs);

y1min = zeros(ns,1);
y1max = zeros(ns,1);
y1init = zeros(ns,1);
y1end = zeros(ns,1);

y2min = zeros(ns,1);
y2max = zeros(ns,1);
y2init = zeros(ns,1);
y2end = zeros(ns,1);

y3min = zeros(ns,1);
y3max = zeros(ns,1);
y3init = zeros(ns,1);
y3end = zeros(ns,1);

y4min = zeros(ns,1);
y4max = zeros(ns,1);
y4init = zeros(ns,1);
y4end = zeros(ns,1);

y5min = zeros(ns,1);
y5max = zeros(ns,1);
y5init = zeros(ns,1);
y5end = zeros(ns,1);

y6min = zeros(ns,1);
y6max = zeros(ns,1);
y6init = zeros(ns,1);
y6end = zeros(ns,1);

dx=[-3;0;0;];
%% ***** Inputs to model ********
% Manipulated Variables
z2=0.8;                         % Injection Valve Opening(u2)
% Disturbances
GOR = 0;                     % Mass GOR from reservoir, NB: MASS, increasing gives: lower BH, stable
P_r = 160;                      % Rservoir pressure [bar]
P_gs = 150;                     % Gass source pressure [bar]
P0 = 15;                        % Separator Pressure [bar]
d0 = [GOR;P_r;P_gs;P0];         % Disturbance Vector
%Constans in Model
M=eye(3);
tf = 6;
options=odeset('Mass',M,'BDF','on','AbsTol',1e-12,'RelTol',1e-12,'MaxStep',1);
%%
for s = 1:ns
    % **** System Inputs ****
    z1 = zs(s);
    u0=[z1;z2;d0];
    % ******** Initializing *************
    [x0,flag,y0] = v1_gaslift_initialize(u0,par);
    % **********************
    %% Check if unstable system
    [A B C D]=central_difference(x0,u0,par); % State Space Model
    [eigval]=real(eig(A));
    unstable=0;
    for i=1:length(eigval)
        if eigval(i) >=0
            unstable=1;
        end
    end
    if(unstable) %If unstable system
        disp(['Z=' num2str(z1)])
        disp('UNSTABLE')
        [t,x]=ode15s(@v1_gaslift_model,[0;3600*tf],x0 + dx,options,u0,'derivatives',par);
        % *** Calculating measurements ****
        n = length(t);
        y = zeros(6,n);
        for k = 1:n
            yt = v1_gaslift_model(t(k),x(k,:),u0,'measurements',par);
            y(:,k) = [yt([5 6 4 2 3]); 1000*yt(6)/yt(7)];
        end
        n1 = floor(n/2);
        y1min(s) = min(y(1,n1:n));
        y1max(s) = max(y(1,n1:n));
        
        y2min(s) = min(y(2,n1:n));
        y2max(s) = max(y(2,n1:n));
        
        y3min(s) = min(y(3,n1:n));
        y3max(s) = max(y(3,n1:n));
        
        y4min(s) = min(y(4,n1:n));
        y4max(s) = max(y(4,n1:n));
        
        y5min(s) = min(y(5,n1:n));
        y5max(s) = max(y(5,n1:n));
        
        y6min(s) = min(y(6,n1:n));
        y6max(s) = max(y(6,n1:n));
        
        clear y x t;
    else
        disp(['Z=' num2str(z1)])
        disp('STABLE')
        y1min(s) = y0(5);
        y1max(s) = y0(5);
        
        y2min(s) = y0(6);
        y2max(s) = y0(6);
        
        y3min(s) = y0(4);
        y3max(s) = y0(4);
        
        y4min(s) = y0(2);
        y4max(s) = y0(2);
        
        y5min(s) = y0(3);
        y5max(s) = y0(3);
        
        y6min(s) = 1000*y0(6)/y0(7);
        y6max(s) = 1000*y0(6)/y0(7);
    end
end

save('Matlab_bif_data.mat','zs','y1min','y1max','y2min','y2max','y3min','y3max','y4min','y4max','y5min','y5max','y6min','y6max');
Steady_State_auto;