clear all
clc

%Seperator conditions
Ps=20;  

%  a=0.0001;
%  b=0.01;
%  c=0.01;
%  C= [a b c]'

%Well Pressures
Pw1=40;    
Pw2=40;
Pw3=40;
Pw4=40;

%Load welldata
load welldata1

%well1=[well1(:,1) well1(:,2) 0.93*well1(:,3) well1(:,4)];         % Well1 - implement -7% disturbance in oil flow
%well2=[well2(:,1) well2(:,2) 0.93*well2(:,3) well2(:,4)];         % Well2 - implement -7% disturbance in oil flow.
%well3=[well3(:,1) well3(:,2) 0.93*well3(:,3) well3(:,4)];         % Well3 - implement -7% disturbance in oil flow.
%well4=[well4(:,1) well4(:,2) 0.93*well4(:,3) well4(:,4)];         % Well4 - implement -7% disturbance in oil flow.

Qw1=wpcn1(Pw1,well1);          % Well flow [g o wt]' (from well preformance curve)
Qw2=wpcn2(Pw2,well2);          
Qw3=wpcn3(Pw3,well3);
Qw4=wpcn4(Pw4,well4);



Q=Qw1 + Qw2 + Qw3 + Qw4;    % Flow in pipeline [g o wt]'

Qo=Q(2);

Qtot=sum(Q) ;               % Total flow in pipeline

dP=Pdrop2(Q);               % Pressure loss model

Pm= Ps + dP;                % Manifold pressure (dP from model)



%============= OPTIMIZATION ==============================================

%First guess

% %  Pw1 Pw2 Pw3 Pw4 Qw1(g) Qw1(o) Qw1(w) Qw2(g) Qw2(o) Qw2(w) Qw3(g) Qw3(o) Qw3(w) Qw4(g) Qw4(o) Qw4(w)  Q(g)  Q(o)  Q(w)    dP    Pm  
x0=[ 100  50  70 50 1920    82     410    1550    440    660    750 340    510   1550    440    660    5770  1302  2240   -3.52  23.52 ];

%Optimal values

%x0=[ 109 82 63.1 82 59.7 4.8 23.9 1148.4 330.1 495.1 493.4 390.6 585.9 1148.4 330.1 4955.1 3300 1055.5 1600 34.2 54.2]; % optimal values 


A=[];
b=[];
Aeq=[];
beq=[];
lb=[20 20 20 20]';
ub=[109 100 100 100]';
opt=optimset('display','iter-detailed')

[P,funval,exitflag,output,lambda] = fmincon(@objF2,x0,A,b,Aeq,beq,lb,ub,@nonlncon2,opt)

optpress=P(1:4)                 % Displays optimal pressure
%nomval=P'