%% A program to solve the non-linear equations of the optimization problem 
%% by  a variant of Gauss Siedel Method

clear all;
%Constants
z=0.02 + 0.1i; Pl1=4;Ql1=1;Pl2=1.5;Ql2=0.75;Z=abs(z);f=angle(z); 
%Initial Values
P1=1; P2=1;ld1=1550;ld2=1550;V1=1.02;V2=1.02;d=-10*pi/180;
iter=1;
error=1000;
Pm = (V1*V2/Z ); S1 = V1*V1/Z ; S2 =V2*V2/Z;
  while ( (max(error)>0.00001)&(iter < 1000 ) )
      Pm = (V1*V2/Z ); S1 = V1*V1/Z ; S2 =V2*V2/Z; 
      
      ld1old=ld1;
      ld2old=ld2;
      P1old=P1;
      P2old=P2;
      dold=d;
      V1old=V1;
      V2old=V2;
      
      %% The nonlinear equations obtained by setting partial derivatives to zero
      %% are re-written in the form x=f(x)
      
ld1= (50 + 600*P1)+100*exp(100*(P1-3));
ld2= (50 + 500*P2)+100*exp(100*(P2-3)) ;
P1=-(- Pl1- S1*cos(f) + Pm*cos(f+d));
P2=-(- Pl2- S2*cos(f) + Pm*cos(f-d));
d=real(asin(ld2*sin(f-d)/ld1))-f; 
V1=(ld1*(V2*cos(f+d))+ld2*(V2*cos(f-d))-100*exp(100*(V1-1.025))*Z+100*exp(100*(0.975-V1))*Z )/ld1/2/cos(f);
V2=(ld1*(V1*cos(f+d))+ld2*(V1*cos(f-d))-100*exp(100*(V2-1.025))*Z+100*exp(100*(0.975-V2))*Z )/ld2/2/cos(f);

%% Error 
error=abs([ld1old-ld1;ld2old-ld2;P1old-P1;P2old-P2;dold-d;V1old-V1;V2old-V2;]);

%% To avoid divergence we make an augmented update

ld1=0.95*ld1old+0.05*ld1;
ld2=0.95*ld2old+0.05*ld2;
P1=0.95*P1old+0.05*P1;
P2=0.95*P2old+0.05*P2;
d=0.95*dold+0.05*d;
V1=0.95*V1old+0.05*V1;
V2=0.95*V2old+0.05*V2;

iter=iter+1;
end

ld1
ld2
P1
P2
d
V1
V2
PFlow12= S1*cos(f) - Pm*cos(f+d)
PFlow21= S2*cos(f) - Pm*cos(f-d)
Q1=Ql1+S1*sin(f) - Pm*sin(f+d) 
Q2=Ql2+S2*sin(f) - Pm*sin(f-d)
C1plusC2=50*P1 + 300*P1*P1 + 50*P2 + 250*P2*P2

