% Solving network for each forecast time

clear all
clc

ef=0.001;
m=2;
Cacc=1.7;

load GasNetworkProperties

% Applying symmetry of matrices

for i=1:n
    for j=1:n
        B(j,i)=B(i,j);
    end
end

for i=1:n
    for j=1:n
        L(j,i)=L(i,j);
    end
end

for i=1:n
    for j=1:n
        KJ(j,i)=KJ(i,j);
    end
end

for i=1:n
    for j=1:n
        D(j,i)=D(i,j);
    end
end

% Calculating KP from Weymouth equattion

Tb=300;
Pb=101;
G=0.65;
Z=0.9;
T=288.5;
for i=1:n
    for j=1:n
        if B(i,j)==1
            KP(i,j)=71358.28*G*T*Z*((Pb/Tb)^2)*L(i,j)/(D(i,j)^(16/3));
        else
            KP(i,j)=0;
        end
    end
end

% Total resistance

K=KP+KJ;

load EX3Solved

% Solving the gas network

for cons=1:1200
    cons
    Q0=1*X1(cons,:);
% Initial guess for nodal pressure

for i=1:n
    Ps(i)=(Pmax*(n-i+1)/n)^2;
end

% Take the pressures at CGS and TBS stations as boundary conditions.

for i=1:NCGS
    Ps(i)=(PCGS(i))^2;
end
for i=NCGS+1:NCGS+NTBS
    Ps(i)=(PTBS(i-NCGS))^2;
end
Ps=Ps';

% Calculating nodal pressure.

df=10*ef;
J1=0*rand(n,n);
fe1=0*rand(1,n);
cc=0;
while df>ef    
%     Calculating Jacobian matrix    
    for i=1:n
        for j=1:n
            if i~=j && B(i,j)==1 && Ps(i)~=Ps(j)
                J1(i,j)=Cacc/(m*((K(i,j)*((abs(Ps(i)-Ps(j)))^(m-1)))^(1/m)));
            elseif i==j
                J1(i,j)=0;
                for j0=1:n
                    if B(i,j0)==1 && Ps(i)~=Ps(j0)
                        J1(i,j)=J1(i,j)-Cacc/(m*((K(i,j0)*((abs(Ps(i)-Ps(j0)))^(m-1)))^(1/m)));
                    else
                    end
                end
            end
        end
    end
    
    for i=NTBS+NCGS+1:n
        for j=NTBS+NCGS+1:n
            i0=i-NTBS-NCGS;
            j0=j-NTBS-NCGS;
            J(i0,j0)=J1(i,j);
        end
    end
    
%     Calculation of the error in objective function, f.
    
    for i=1:n
        fe1(i)=-Q0(i);
        for j=1:n
            if B(i,j)==1 && Ps(i)~=Ps(j)
                e(i,j)=(Ps(j)-Ps(i))/abs(Ps(i)-Ps(j));
                fe1(i)=fe1(i)+e(i,j)*(abs(Ps(i)-Ps(j))/K(i,j))^(1/m);
            else
            end
        end
    end
    for i=NTBS+NCGS+1:n
        i0=i-NTBS-NCGS;
        fe(i0)=fe1(i);
    end
    cc=cc+1;
    df=norm(fe);
    
%     Calculating nodal pressure
    
    dPs=-inv(J)*fe';
    
    for i=NTBS+NCGS+1:n
        i0=i-NTBS-NCGS;
        Ps1(i0,1)=Ps(i);
    end
    
    Ps1=Ps1+dPs;
        
    for i=NTBS+NCGS+1:n
        Ps(i)=Ps1(i-NTBS-NCGS);
    end
    
%     for i=1:NCGS
%         Ps(i)=(PCGS(i))^2;
%     end
%     for i=NCGS+1:NCGS+NTBS
%         Ps(i)=(PTBS(i-NCGS))^2;
%     end
end
for i=1:n
    P(i)=sqrt(Ps(i));
end
Pt(cons,:)=P;
end

% Plotting outputs

% Excluding CGS constant and high pressure to improve resolution of the graphs

for j=2:n
Pm(:,j-1)=Pt(:,j);
end

a0=1:length(Pm(1,:));
b0=1:length(Pm(:,1));
[a,b]=meshgrid(a0,b0);
mesh(a,b,Pm)

plot3(XN,YN,Pt(254,:))
for i=1:n
    hold on
    plot(XN(i),YN(i),'o','LineWidth',1,'MarkerEdgeColor','k','MarkerFaceColor','r','MarkerSize',4)
end
for i=1:n
    for j=1:i
        if B(i,j)==1
            [x,y]=gnetplot(XN(i),XN(j),YN(i),YN(j));
            hold on
            plot(x,y,'k-','LineWidth',1,'MarkerEdgeColor','k')
        else
        end
    end
end
