function PipeNetworkSolver
clc; clear all;

Parameters.Demand=[0 300 150 10]; %[m3/h]
Parameters.NodalHeight=[0 0 10 30]; %[m]
Parameters.FrictionCoefficients=[0 0.018 0.020 0.025 0]; %[-]
Parameters.PipeLength=[0 2000 1000 3400 0]; %[m]
Parameters.PipeDiameter=[0 400 300 200 0]; %[mm]
Parameters.AmbientPressure=1; %[bar]
Parameters.Density=1000; %[kg/m3]
Parameters.Gravity=9.81; %[m/s2]
Parameters.ReservoirBottomHeight=60; %[m]
Parameters.ReservoirWaterLevel=2; %[m]

InitialGuess=[3 3 3 3 250 250 250 250 250]'; %[bar] or [m3/h]



Results=SolveThePipeNetworkSystem(InitialGuess,Parameters);
fprintf('The data:\n');
    for k=1:length(Results)
        disp(Results(k));
    end



function f=SystemOfEquations(x,Parameters)
f=zeros(9,1);

f(1)=x(5)-x(6)-x(7)-Parameters.Demand(2);
f(2)=x(7)-x(8)-Parameters.Demand(3);
f(3)=x(8)+x(6)-x(9)-Parameters.Demand(4);

f(4)=x(1)-Parameters.AmbientPressure;

f(5)=(x(2)-x(1)) - Parameters.Density*Parameters.Gravity*PumpCharacteristicCurve(x(5))/100000;

f(6)=(x(4)-x(2)) + Parameters.Density*Parameters.Gravity*(Parameters.NodalHeight(4)-Parameters.NodalHeight(2))/100000 + ...
     8*Parameters.FrictionCoefficients(2)*Parameters.PipeLength(2)*Parameters.Density*x(6)*abs(x(6))/Parameters.PipeDiameter(2)^5/pi^2 * 771.6;
f(7)=(x(3)-x(2)) + Parameters.Density*Parameters.Gravity*(Parameters.NodalHeight(3)-Parameters.NodalHeight(2))/100000 + ...
     8*Parameters.FrictionCoefficients(3)*Parameters.PipeLength(3)*Parameters.Density*x(7)*abs(x(7))/Parameters.PipeDiameter(3)^5/pi^2 * 771.6;
f(8)=(x(4)-x(3)) + Parameters.Density*Parameters.Gravity*(Parameters.NodalHeight(4)-Parameters.NodalHeight(3))/100000 + ...
     8*Parameters.FrictionCoefficients(4)*Parameters.PipeLength(4)*Parameters.Density*x(8)*abs(x(8))/Parameters.PipeDiameter(4)^5/pi^2 * 771.6;

f(9)=(Parameters.AmbientPressure-x(4)) + Parameters.Density*Parameters.Gravity*(Parameters.ReservoirBottomHeight+Parameters.ReservoirWaterLevel-Parameters.NodalHeight(4))/100000;

function H = PumpCharacteristicCurve(Q)
% H[m], Q[m3/h]

H=250-0.0002*Q^2;



function x=SolveThePipeNetworkSystem(InitialGuess,Parameters)
ErrorTolerance=1e-6;
MaximumIteration=100;

[Error, Index]=max(abs(SystemOfEquations(InitialGuess,Parameters)));

NumberOfSteps=1; xOld=InitialGuess;
while Error>ErrorTolerance
    Jacobian=NumericalJacobian(xOld,Parameters);
    xNew=xOld-inv(Jacobian)*SystemOfEquations(xOld,Parameters);
    
    [Error, Index]=max(abs(SystemOfEquations(xNew,Parameters)));
        fprintf('%2d: Error=%5.3e, %2d\n',NumberOfSteps,Error,Index);
    
    xOld=xNew; x=xNew;
    NumberOfSteps=NumberOfSteps+1;

    if NumberOfSteps==MaximumIteration+1
        fprintf('Too many number of iteration, please try an other initial guess!!!\n');
        Error=ErrorTolerance/2;
    end
    if Error<ErrorTolerance
        fprintf('Iteration converged!!!\n');
    end
end

function dfdx=NumericalJacobian(x,Parameters)
dfdx=zeros(length(x));
f=SystemOfEquations(x,Parameters);
    for i=1:length(x)
       x1=x;
       dx=x1(i)*0.001;
       if abs(dx)<1e-6
           dx=1e-6;
       end
       x1(i)=x1(i)+dx; 
       f1=SystemOfEquations(x1,Parameters);
       dfdx(:,i)=(f1-f)/dx;
    end