<<

. 16
( 17)



>>

end

[cost] = CostGeneration(cost,A,B,C,PGEN,ngn);
[suma] = CheckConvexi¬cation(nbb,Cc,VM,VOLTAGE);

%% ******************************************************************
%% If the total cost is within tolerance *
%% and the optimality conditions are satis¬ed, then *
%% the process can be stopped. *
%% ******************************************************************

if((in0==0)&(in1==0)&(in2==0)&(in3==0)&(in4==0)&(RelVol==0)& ...
(RelGen==0))
eps=LastCost-cost;
if(eps<0),eps=-eps;end
if(TolEps==1e-8)
if(eps<TolEps),Optimo=0;end
end
end

LastCost=cost;

if((in0==0)&(in1==0)&(in2==0)&(in3==0)&(in4==0)&(RelVol==0)...
&(RelGen==0)&(Optimo==0))
VM
VA=VA*180/pi
PGEN
LambdaP
LambdaQ
fprintf(™\n=== Objective Function Value ===™);
fprintf(™\n f = %%12.8f $/hr ™, cost);
fprintf(™\n™);
fprintf(™\n™);
fprintf(™\n************************************************™);
fprintf(™\n* End of main iteration *™);
fprintf(™\n* *™);
fprintf(™\n************************************************™);
break; %% This instruction breaks the main iterations loop
end
374 APPENDIX C

%% ******************************************************************
%% It checks for the Dz being too small. If true then *
%% Ckv=const, Ckg=const and IntIter=30 *
%% ******************************************************************
if((MaxDz<0.00001)&(in1==0)&(in2==0)&(in3==0)&(in4==0)&...
(RelVol==0)& (RelGen==0))
NIntIter==30;
Ckv=Ckv/1.3;
Ckg=Ckg/1.3;
end


end %% End of the Main PROGRAM



function [GKK,BKK,GKM,BKM] = YBus(tlsend,tlrec,tlresis,tlreac,...
tlsuscep,tlcond,ntl,nbb);
%% Transmission lines contribution


GKK=zeros(ntl,1);
BKK=zeros(ntl,1);
GKM=zeros(ntl,1);
BKM=zeros(ntl,1);
for ii = 1: ntl
denom = tlresis(ii)^2+tlreac(ii)^2;
GKK(ii) = GKK(ii) + tlresis(ii)/denom + 0.5*tlcond(ii);
BKK(ii) = BKK(ii) - tlreac(ii)/denom + 0.5*tlsuscep(ii);
GKM(ii) = GKM(ii) - tlresis(ii)/denom;
BKM(ii) = BKM(ii) + tlreac(ii)/denom;
end
return;%% End of Ybus



function [PGEN,Ckg,LAMBDAP] = InitialDispatch(nbb,ngn,A,B,C,PLOAD,...
PGEN,PMIN,PMAX,Ckg,LambdaP);
sum1=0.0;sum2=0.0;sum3=0.0;sum4=0;sum5=0;
lambda=0.0;
Ckg=0;
for(ii= 1:ngn)
sum1=sum1+B(ii)/C(ii);
sum2=sum2+1/C(ii);
if(C(ii)>=Ckg),Ckg=C(ii);end
end
for(ii= 1:nbb)
sum3=sum3+PLOAD(ii);
end
sum3=1.03*sum3;
375
COMPUTER PROGRAM FOR OPTIMAL POWER FLOW SOLUTIONS

for( ii=1:ngn)
sum4=sum4+PMIN(ii);
sum5=sum5+PMAX(ii);
end
if(sum3>sum5|sum3<sum4)
fprintf(™\n************************************************™);
fprintf(™\n* THERE IS NO SOLUTION *™);
fprintf(™\n* LOAD<PMIN OR LOAD>PMAX *™);
fprintf(™\n* *™);
fprintf(™\n************************************************™);
end


lambda=(2*sum3+sum1)/sum2;
for(ii= 1:nbb)
LambdaP(ii)=lambda;
end
for(ii= 1:ngn)
sum1=1.0/(2*C(ii))*(lambda-B(ii));
PGEN(ii)=sum1;
end
return; %%End of InitialDispatch



function[cost] = CostGeneration(cost,A,B,C,PGEN,ngn)
cost=0.0;
for( ii = 1: ngn)
cost=cost+A(ii)+B(ii)*PGEN(ii)+C(ii)*PGEN(ii)*PGEN(ii);
end
return; %%End of CostGeneration function


function [Pbus,Qbus] = CalculatedPowers(nbb,VM,VA,GKK,BKK,GKM,BKM,
ntl,tlsend,tlrec);
Pbus= zeros(nbb,1);
Qbus = zeros(nbb,1);
V =zeros(2,1);
A =zeros(2,1);
for ii = 1: ntl
send = tlsend(ii); rece = tlrec(ii);
V(1)= VM(send); V(2)= VM(rece);
A(1)= VA(send); A(2)= VA(rece);
angle=A(1)-A(1);
Pbus(send)=Pbus(send)+V(1)*V(1)*(GKK(ii)*cos(angle)+...
BKK(ii)*sin(angle));
Qbus(send)=Qbus(send)+V(1)*V(1)*(GKK(ii)*sin(angle)-BKK(ii)...
*cos(angle));angle=A(1)-A(2);
376 APPENDIX C

Pbus(send)=Pbus(send)+V(1)*V(2)*(GKM(ii)*cos(angle)+...
BKM(ii)*sin(angle));
Qbus(send)=Qbus(send)+V(1)*V(2)*(GKM(ii)*sin(angle)-BKM(ii)...
*cos(angle));
angle=A(2)-A(2);
Pbus(rece)=Pbus(rece)+V(2)*V(2)*(GKK(ii)*cos(angle)+...
BKK(ii)*sin(angle));
Qbus(rece)=Qbus(rece)+V(2)*V(2)*(GKK(ii)*sin(angle)-BKK(ii)...
*cos(angle));
angle=A(2)-A(1);
Pbus(rece)=Pbus(rece)+V(2)*V(1)*(GKM(ii)*cos(angle)+...
BKM(ii)*sin(angle));
Qbus(rece)=Qbus(rece)+V(2)*V(1)*(GKM(ii)*sin(angle)-BKM(ii)...
*cos(angle));
end
return; %%End of CalculatedPowers function

function [QGEN] = ReactivePowerGenerators(QGEN,ngn,Qbus,QLOAD,...
genbus,bustype);
for ii = 1: ngn
bgen=genbus(ii);
if(bustype(bgen)==1 | bustype(bgen)==2)
QGEN(ii)=Qbus(bgen)+QLOAD(bgen);
end
end
return; %%End of ReactivePowerGenerators function

function [Violation,ActivedQ,statusgen,SetGenQ,IndGenQ] = Voltage...
(nbb,ngn,genbus,QGEN,QLOAD,Violation,ActivedQ,vmin,vmax,TolVoltage,...
VM, Qbus,statusgen,bustype,QMAX,QMIN,SetGenQ,IndGenQ,iterOpf);
if (iterOpf >=2 )
for ii = 1: ngn
bgen=genbus(ii);
if(ActivedQ(ii)==1)
bgen=bgen;
else
Violation(bgen)=1;
if(((vmin(bgen)-TolVoltage)<VM(bgen))&(VM(bgen)<(vmax(bgen)+...
TolVoltage)))
Violation(bgen)=0;
end
if((VM(bgen)<0.5)j(VM(bgen)>1.5))
fprintf(™\n************************************************™);
fprintf(™\n* *™);
fprintf(™\n* UNFEASIBLE SOLUTION *™);
fprintf(™\n* *™);
fprintf(™\n************************************************™);
end
377
COMPUTER PROGRAM FOR OPTIMAL POWER FLOW SOLUTIONS

if(Violation(bgen)==0)
Rpower=QLOAD(bgen)+Qbus(bgen);
if(statusgen(ii)==1)
if((Rpower>(QMAX(ii)+0.001)) |(Rpower<(QMIN(ii)...
-0.001)))
QGEN(ii)=-QGEN(ii);
bustype(bgen)=3; %% load bus
if(Rpower>QMAX(ii))
QGEN(ii)=QMAX(ii);
statusgen(ii)=0;
else
QGEN(ii)=QMIN(ii);
statusgen(ii)=2;
end
end
end
end
end
if(bustype(bgen)==1 | bustype(bgen)==2)
SetGenQ=SetGenQ+1;
IndGenQ(SetGenQ)=ii;
end
end
else
for(ii = 1: ngn )
SetGenQ=SetGenQ+1;
IndGenQ(SetGenQ)=ii;
end
end
return; %% End of Voltage function



function [ANGLE,VOLTAGE,LAMBDAP,LAMBDAQ,PGENERATED] = VectorAux(nbb,...
ngn,VM,VA,LambdaP,LambdaQ,PGEN,vmax,vmin,ANGLE,VOLTAGE,LAMBDAP,...
LAMBDAQ,PGENERATED);
for ii = 1: nbb
ANGLE(ii)=VA(ii);
VOLTAGE(ii)=VM(ii);
LAMBDAP(ii)=LambdaP(ii);
LAMBDAQ(ii)=LambdaQ(ii);
if(VM(ii) >=vmax(ii))
VOLTAGE(ii)=vmax(ii);
elseif(VM(ii)<=vmin(ii))
VOLTAGE(ii)=vmin;
end
end
for ii = 1: ngn
378 APPENDIX C

PGENERATED(ii)=PGEN(ii);
end
return; %%End of VectorAux function


function [Hessian,grad] = MatrixW(nbb,bustype,ngn,GKK,BKK,GKM,BKM,...
tlsend,tlrec,LambdaP,LambdaQ,VM,VA,ntl);
Hessian = zeros(4*nbb+ngn,4*nbb+ngn);
grad = zeros(4*nbb+ngn,1);
for jj = 1: ntl
send= tlsend(jj);
rece= tlrec(jj);
%% Load-Load
if(bustype(send)== 3 & bustype(rece)==3)
nb=nbb;
i1=1;
i2=2;


for ii = 1: 2
v1= VM(send); v2= VM(rece);
A1= VA(send); A2= VA(rece);
LP1=LambdaP(send); LP2=LambdaP(rece);
LQ1=LambdaQ(send); LQ2=LambdaQ(rece);
difAng12=A1-A2;
difAng21=A2-A1;
G12=GKM(jj); B12=BKM(jj); G11=GKK(jj); B11=BKK(jj);
G21=G12; B21=B12; G22=G11; B22=B11;
Hkm=(G12*sin(difAng12)-B12*cos(difAng12));
Nkm=(G12*cos(difAng12)+B12*sin(difAng12));
Hmk=(G21*sin(difAng21)-B21*cos(difAng21));
Nmk=(G21*cos(difAng21)+B21*sin(difAng21));


%% Diagonal element
Hessian(send,send)=Hessian(send,send)-LP1*v1*v2*Nkm-LP2*v1...
*v2*Nmk -LQ1*v1*v2*Hkm-LQ2*v1*v2*Hmk;
Hessian(send,1*nb+send)=Hessian(send,1*nb+send)-...
LP1*v2*Hkm+LP2*v2*Hmk+LQ1*v2*Nkm-LQ2*v2*Nmk;
Hessian(send,2*nb+send)=Hessian(send,2*nb+send)-v1*v2*Hkm;
Hessian(send,3*nb+send)=Hessian(send,3*nb+send)+v1*v2*Nkm;


Hessian(nb+send,send)=Hessian(nb+send,send)-...
LP1*v2*Hkm+LP2*v2*Hmk+LQ1*v2*Nkm-LQ2*v2*Nmk;
Hessian(nb+send,nb+send)=Hessian(nb+send,nb+send)...
+LP1*2*G11-LQ1*2*B11;
Hessian(nb+send,2*nb+send)=Hessian(nb+send,2*nb+send)...
+2*v1*G11+v2*Nkm;
379
COMPUTER PROGRAM FOR OPTIMAL POWER FLOW SOLUTIONS

Hessian(nb+send,3*nb+send)=Hessian(nb+send,3*nb+send)-...
2*v1*B11+v2*Hkm;


Hessian(2*nb+send,send)=Hessian(2*nb+send,send)-v1*v2*Hkm;
Hessian(2*nb+send,nb+send)=Hessian(2*nb+send,nb+send)+2*v1*G11...
+v2*Nkm;
Hessian(2*nb+send,2*nb+send)=Hessian(2*nb+send,2*nb+send);
Hessian(2*nb+send,3*nb+send)=Hessian(2*nb+send,3*nb+send);


Hessian(3*nb+send,send)=Hessian(3*nb+send,send)+v1*v2*Nkm;
Hessian(3*nb+send,nb+send)=Hessian(3*nb+send,nb+send)-...
2*v1*B11+v2*Hkm;
Hessian(3*nb+send,2*nb+send)=Hessian(3*nb+send,2*nb+send);
Hessian(3*nb+send,3*nb+send)=Hessian(3*nb+send,3*nb+send);


%% Off diagonal element
Hessian(send,rece)=Hessian(send,rece)+LP1*v1*v2*Nkm+...
LP2*v1*v2*Nmk+LQ1*v1*v2*Hkm+LQ2*v1*v2*Hmk;
Hessian(send,nb+rece)=Hessian(send,nb+rece)-...
LP1*v1*Hkm+LP2*v1*Hmk+LQ1*v1*Nkm-LQ2*v1*Nmk;
Hessian(send,2*nb+rece)=Hessian(send,2*nb+rece)+v1*v2*Hmk;
Hessian(send,3*nb+rece)=Hessian(send,3*nb+rece)-v1*v2*Nmk;


Hessian(nb+send,rece)=Hessian(nb+send,rece)+LP1*v2...
*Hkm-LP2*v2*Hmk-LQ1*v2*Nkm+LQ2*v2*Nmk;
Hessian(nb+send,nb+rece)=Hessian(nb+send,nb+rece)+LP1*Nkm+...
LP2*Nmk+LQ1*Hkm+LQ2*Hmk;
Hessian(nb+send,2*nb+rece)=Hessian(nb+send,2*nb+rece)+v2*Nmk;
Hessian(nb+send,3*nb+rece)=Hessian(nb+send,3*nb+rece)+v2*Hmk;


Hessian(2*nb+send,rece)=Hessian(2*nb+send,rece)+v1*v2*Hkm;
Hessian(2*nb+send,nb+rece)=Hessian(2*nb+send,nb+rece)+v1*Nkm;
Hessian(2*nb+send,2*nb+rece)=Hessian(2*nb+send,2*nb+rece);
Hessian(2*nb+send,3*nb+rece)=Hessian(2*nb+send,3*nb+rece);


Hessian(3*nb+send,rece)=Hessian(3*nb+send,rece)-v1*v2*Nkm;
Hessian(3*nb+send,nb+rece)=Hessian(3*nb+send,nb+rece)+v1*Hkm;
Hessian(3*nb+send,2*nb+rece)=Hessian(3*nb+send,2*nb+rece);
Hessian(3*nb+send,3*nb+rece)=Hessian(3*nb+send,3*nb+rece);

grad(send)=grad(send)-(-LP1*v1*v2*Hkm+LP2*v1*v2*Hmk+LQ1*v1*v2...
*Nkm-LQ2*v1*v2*Nmk);
grad(nb+send)=grad(nb+send)-...
(+LP1*(2*v1*G11+v2*Nkm)+LP2*v2*Nmk+LQ1*(2*v1*B11+v2*Hkm)...
+LQ2*v2*Hmk);
380 APPENDIX C

grad(2*nb+send)=0.0;
grad(3*nb+send)=0.0;


itemp=send; send=rece; rece=itemp;
itemp=i1; i1=i2; i2=itemp;
rtemp=A1; A1=A2; A2=rtemp;
end
end


%% PV-load buses or Slack-load buses
if((bustype(send)== 2 & bustype(rece)==3)|...
(bustype(send)== 3 & bustype(rece)==2)|...
(bustype(send)== 1 & bustype(rece)==3)|...
(bustype(send)== 3 & bustype(rece)==1))
nb=nbb;
i1=1;
i2=2;


for ii = 1: 2
v1= VM(send); v2= VM(rece);
A1= VA(send); A2= VA(rece);
LP1=LambdaP(send); LP2=LambdaP(rece);
LQ1=LambdaQ(send); LQ2=LambdaQ(rece);
difAng12=A1-A2; difAng21=A2-A1;
G12=GKM(jj); B12=BKM(jj); G11=GKK(jj); B11=BKK(jj);
G21=G12; B21=B12; G22=G11; B22=B11;

Hkm=(G12*sin(difAng12)-B12*cos(difAng12));
Nkm=(G12*cos(difAng12)+B12*sin(difAng12));
Hmk=(G21*sin(difAng21)-B21*cos(difAng21));
Nmk=(G21*cos(difAng21)+B21*sin(difAng21));

if(bustype(send)== 2 | bustype(send)== 1)
LQ1=0;
else
LQ2=0;
end



%% Diagonal elements
Hessian(send,send) = Hessian(send,send)-LP1*v1*v2...
*Nkm-LP2*v1*v2*Nmk-LQ1*v1*v2*Hkm-LQ2*v1*v2*Hmk;
Hessian(send,1*nb+send) = Hessian(send,1*nb+send)-...
LP1*v2*Hkm+LP2*v2*Hmk+LQ1*v2*Nkm-LQ2*v2*Nmk;
Hessian(send,2*nb+send) = Hessian(send,2*nb+send)-v1*v2*Hkm;
381
COMPUTER PROGRAM FOR OPTIMAL POWER FLOW SOLUTIONS

Hessian(send,3*nb+send) = Hessian(send,3*nb+send)+v1*v2*Nkm;
Hessian(nb+send,send) = Hessian(nb+send,send)-...
LP1*v2*Hkm+LP2*v2*Hmk+LQ1*v2*Nkm-LQ2*v2*Nmk;
Hessian(nb+send,nb+send) = Hessian(nb+send,nb+send)+...
LP1*2*G11-LQ1*2*B11;
Hessian(nb+send,2*nb+send) = Hessian(nb+send,2*nb+send)+...
2*v1*G11+v2*Nkm;
Hessian(nb+send,3*nb+send) = Hessian(nb+send,3*nb+send)-...
2*v1*B11+v2*Hkm;


Hessian(2*nb+send,send) = Hessian(2*nb+send,send)-v1*v2*Hkm;
Hessian(2*nb+send,nb+send) = Hessian(2*nb+send,nb+send)+...
2*v1*G11+v2*Nkm;
Hessian(2*nb+send,2*nb+send)=Hessian(2*nb+send,2*nb+send);
Hessian(2*nb+send,3*nb+send)=Hessian(2*nb+send,3*nb+send);


Hessian(3*nb+send,send)=Hessian(3*nb+send,send)+v1*v2*Nkm;
Hessian(3*nb+send,nb+send)=Hessian(3*nb+send,nb+send)-...
2*v1*B11+v2*Hkm;
Hessian(3*nb+send,2*nb+send)=Hessian(3*nb+send,2*nb+send);
Hessian(3*nb+send,3*nb+send)=Hessian(3*nb+send,3*nb+send);


%% Off-Diagonal elements
Hessian(send,rece)=Hessian(send,rece)+LP1*v1*v2*Nkm+LP2*v1*v2* ...
Nmk+LQ1*v1*v2*Hkm+LQ2*v1*v2*Hmk;
Hessian(send,nb+rece)=Hessian(send,nb+rece)-...
LP1*v1*Hkm+LP2*v1*Hmk+LQ1*v1*Nkm-LQ2*v1*Nmk;
Hessian(send,2*nb+rece)=Hessian(send,2*nb+rece)+v1*v2*Hmk;
Hessian(send,3*nb+rece)=Hessian(send,3*nb+rece)-v1*v2*Nmk;


Hessian(nb+send,rece)=Hessian(nb+send,rece)+LP1*v2...
*Hkm-LP2*v2*Hmk-LQ1*v2*Nkm+LQ2*v2*Nmk;
Hessian(nb+send,nb+rece)=Hessian(nb+send,nb+rece)+LP1*Nkm+LP2*...
Nmk+LQ1*Hkm+LQ2*Hmk;
Hessian(nb+send,2*nb+rece)=Hessian(nb+send,2*nb+rece)+v2*Nmk;
Hessian(nb+send,3*nb+rece)=Hessian(nb+send,3*nb+rece)+v2*Hmk;


Hessian(2*nb+send,rece)=Hessian(2*nb+send,rece)+v1*v2*Hkm;
Hessian(2*nb+send,nb+rece)=Hessian(2*nb+send,nb+rece)+v1*Nkm;
Hessian(2*nb+send,2*nb+rece)=Hessian(2*nb+send,2*nb+rece);
Hessian(2*nb+send,3*nb+rece)=Hessian(2*nb+send,3*nb+rece);


Hessian(3*nb+send,rece)=Hessian(3*nb+send,rece)-v1*v2*Nkm;
Hessian(3*nb+send,nb+rece)=Hessian(3*nb+send,nb+rece)+v1*Hkm;
Hessian(3*nb+send,2*nb+rece)=Hessian(3*nb+send,2*nb+rece);
382 APPENDIX C

Hessian(3*nb+send,3*nb+rece)=Hessian(3*nb+send,3*nb+rece);

grad(send)=grad(send)-(-LP1*v1*v2*Hkm+LP2*v1*v2*Hmk+ ...
LQ1*v1*v2*Nkm-LQ2*v1*v2*Nmk);
grad(nb+send)=grad(nb+send)(+LP1*(2*v1*G11+v2*Nkm)+LP2*v2*Nmk+...
LQ1*(-2*v1*B11+v2*Hkm)+LQ2*v2*Hmk);
grad(2*nb+send)=0.0;
grad(3*nb+send)=0.0;


itemp=send; send=rece; rece=itemp;
itemp=i1; i1=i2; i2=itemp;
rtemp=A1; A1=A2; A2=rtemp;
end
end
%% Slack-PV or PV-PV
if((bustype(send)== 1 & bustype(rece)==2)|...
(bustype(send)== 2 & bustype(rece)==1)|...
(bustype(send)== 2 & bustype(rece)==2))
nb=nbb;
i1=1;
i2=2;
for ii = 1: 2
v1= VM(send); v2= VM(rece);
A1= VA(send); A2= VA(rece);
LP1=LambdaP(send); LP2=LambdaP(rece);
LQ1=0; LQ2=0;
difAng12=A1-A2; difAng21=A2-A1;
G12=GKM(jj); B12=BKM(jj); G11=GKK(jj); B11=BKK(jj);
G21=G12; B21=B12; G22=G11; B22=B11;
Hkm=(G12*sin(difAng12)-B12*cos(difAng12));
Nkm=(G12*cos(difAng12)+B12*sin(difAng12));
Hmk=(G21*sin(difAng21)-B21*cos(difAng21));
Nmk=(G21*cos(difAng21)+B21*sin(difAng21));


%% Diagonal elements
Hessian(send,send)=Hessian(send,send)-LP1*v1*v2*Nkm-LP2*v1*v2*...
Nmk-LQ1*v1*v2*Hkm-LQ2*v1*v2*Hmk;
Hessian(send,1*nb+send)=Hessian(send,1*nb+send)-...
LP1*v2*Hkm+LP2*v2*Hmk+LQ1*v2*Nkm-LQ2*v2*Nmk;
Hessian(send,2*nb+send)=Hessian(send,2*nb+send)-v1*v2*Hkm;
Hessian(send,3*nb+send)=Hessian(send,3*nb+send)+v1*v2*Nkm;


Hessian(nb+send,send)=Hessian(nb+send,send)-...
LP1*v2*Hkm+LP2*v2*Hmk+LQ1*v2*Nkm-LQ2*v2*Nmk;
Hessian(nb+send,nb+send)=Hessian(nb+send,nb+send)+LP1*2*G11-...
LQ1*2*B11;
383
COMPUTER PROGRAM FOR OPTIMAL POWER FLOW SOLUTIONS

Hessian(nb+send,2*nb+send)=Hessian(nb+send,2*nb+send)+2*v1*G11...
+v2*Nkm;
Hessian(nb+send,3*nb+send)=Hessian(nb+send,3*nb+send)-...
2*v1*B11+v2*Hkm;


Hessian(2*nb+send,send)=Hessian(2*nb+send,send)-v1*v2*Hkm;
Hessian(2*nb+send,nb+send)=Hessian(2*nb+send,nb+send)+2*v1*G11...
+v2*Nkm;
Hessian(2*nb+send,2*nb+send)=Hessian(2*nb+send,2*nb+send);
Hessian(2*nb+send,3*nb+send)=Hessian(2*nb+send,3*nb+send);


Hessian(3*nb+send,send)=Hessian(3*nb+send,send)+v1*v2*Nkm;
Hessian(3*nb+send,nb+send)=Hessian(3*nb+send,nb+send)...
-2*v1*B11+v2*Hkm;
Hessian(3*nb+send,2*nb+send)=Hessian(3*nb+send,2*nb+send);
Hessian(3*nb+send,3*nb+send)=Hessian(3*nb+send,3*nb+send);


%% off-diagonal elements
Hessian(send,rece)=Hessian(send,rece)+LP1*v1*v2*Nkm+...
LP2*v1*v2*Nmk+LQ1*v1*v2*Hkm+LQ2*v1*v2*Hmk;
Hessian(send,nb+rece)=Hessian(send,nb+rece)-...
LP1*v1*Hkm+LP2*v1*Hmk+LQ1*v1*Nkm-LQ2*v1*Nmk;
Hessian(send,2*nb+rece)=Hessian(send,2*nb+rece)+v1*v2*Hmk;
Hessian(send,3*nb+rece)=Hessian(send,3*nb+rece)-v1*v2*Nmk;


Hessian(nb+send,rece)=Hessian(nb+send,rece)+LP1*v2...
*Hkm-LP2*v2*Hmk-LQ1*v2*Nkm+LQ2*v2*Nmk;
Hessian(nb+send,nb+rece)=Hessian(nb+send,nb+rece)+LP1...
*Nkm+LP2*Nmk+LQ1*Hkm+LQ2*Hmk;
Hessian(nb+send,2*nb+rece)=Hessian(nb+send,2*nb+rece)+v2*Nmk;
Hessian(nb+send,3*nb+rece)=Hessian(nb+send,3*nb+rece)+v2*Hmk;


Hessian(2*nb+send,rece)=Hessian(2*nb+send,rece)+v1*v2*Hkm;
Hessian(2*nb+send,nb+rece)=Hessian(2*nb+send,nb+rece)+v1*Nkm;
Hessian(2*nb+send,2*nb+rece)=Hessian(2*nb+send,2*nb+rece);
Hessian(2*nb+send,3*nb+rece)=Hessian(2*nb+send,3*nb+rece);


Hessian(3*nb+send,rece)=Hessian(3*nb+send,rece)-v1*v2*Nkm;
Hessian(3*nb+send,nb+rece)=Hessian(3*nb+send,nb+rece)+v1*Hkm;
Hessian(3*nb+send,2*nb+rece)=Hessian(3*nb+send,2*nb+rece);
Hessian(3*nb+send,3*nb+rece)=Hessian(3*nb+send,3*nb+rece);

grad(send)=grad(send)-(-LP1*v1*v2*Hkm+LP2*v1*v2*Hmk+LQ1*v1*v2...
*Nkm-LQ2*v1*v2*Nmk);
grad(nb+send)=grad(nb+send)-(+LP1*(2*v1*G11+v2*Nkm)+...
384 APPENDIX C

LP2*v2*Nmk+LQ1*(-2*v1*B11+v2*Hkm)+LQ2*v2*Hmk);
grad(2*nb+send)=0.0;
grad(3*nb+send)=0.0;

itemp=send; send=rece; rece=itemp;
itemp=i1; i1=i2; i2=itemp;
rtemp=A1; A1=A2; A2=rtemp;
end
end
end
return; %%End of MatrixW function

function [Hessian,grad] = MatrixWGen(nbb,ntl,ngn,genbus,B,C,LambdaP,...
PGEN,Hessian,grad);
for jj = 1: ngn
Gbus=genbus(jj);
Hessian(4*nbb+jj,4*nbb+jj)=Hessian(4*nbb+jj,4*nbb+jj)+2*C(jj);
%% // - - - - - Out of Diagonal
Hessian(4*nbb+jj,2*nbb+Gbus)=Hessian(4*nbb+jj,2*nbb+Gbus)-1.0;
Hessian(2*nbb+Gbus,4*nbb+jj)=Hessian(2*nbb+Gbus,4*nbb+jj)-1.0;
LP1=LambdaP(Gbus);
PGenI=PGEN(jj);
grad(4*nbb+jj)= grad(4*nbb+jj)-(B(jj)+2*C(jj)*PGenI-LP1);
end
return; %%End of MatrixWGen function

function [Hessian,grad] = Mismatch(nbb,ngn,genbus,PGEN,QGEN,PLOAD,...
QLOAD,Pbus,Qbus,Hessian,grad);
AUXP = zeros(nbb,1);
AUXQ = zeros(nbb,1);
for jj = 1: ngn
AUXP(genbus(jj)) = PGEN(jj);
AUXQ(genbus(jj)) = QGEN(jj);
end
for ii= 1: nbb
grad(2*nbb+ii)=0;
grad(3*nbb+ii)=0;
grad(2*nbb+ii)=grad(2*nbb+ii)-(Pbus(ii)-AUXP(ii)+PLOAD(ii));
grad(3*nbb+ii)=grad(3*nbb+ii)-(Qbus(ii)-AUXQ(ii)+QLOAD(ii));
end
return; %%End of Mismatch function

function [Hessian,grad] = PenaltyFunctionQ(nbb,BigNumber,Hessian,...
grad,genbus,IndGenQ,SetGenQ,LambdaQ,ActivedQ);
for(ii = 1: SetGenQ)
in1=IndGenQ(ii);
send=genbus(in1);
385
COMPUTER PROGRAM FOR OPTIMAL POWER FLOW SOLUTIONS

Hessian(3*nbb+send,3*nbb+send)=Hessian(3*nbb+send,3*nbb+send)+...
2*BigNumber;
ActivedQ(in1)=2;
LQ1=LambdaQ(send);
grad(3*nbb+send)=grad(3*nbb+send)-2*BigNumber*LQ1;
end
return; %%End of PenaltyFunctionQ function

function [Hessian,grad] = Convexi¬cacion(nbb,Hessian,grad,Cc,vmax,...
vmin,VOLTAGE,StatusVoltage,VM);
%% Slack bus must be the ¬rst node
Hessian(1,1)=Hessian(1,1)+10e10;
%% Voltage magnitudes
for ii= 1: nbb
Vol=VM(ii);
VOLD=VOLTAGE(ii);
VolMax=vmax(ii);
VolMin=vmin(ii);
StatV=StatusVoltage(ii);
Hessian(nbb+ii,nbb+ii)=Hessian(nbb+ii,nbb+ii)+Cc;
if(StatV==0)
grad(nbb+ii)=grad(nbb+ii)-(Cc*(Vol-VolMax));
elseif(StatV==1)
grad(nbb+ii)=grad(nbb+ii)-(Cc*(Vol-VOLD));
elseif(StatV==2)
grad(nbb+ii)=grad(nbb+ii)-(Cc*(Vol-VolMin));
end
end
return; %%End of Convexi¬cation function


function [Hessian,grad,ActivedV,ActivedP] = AugmentedLagrangian(nbb,...
Hessian,grad,Ckv,Ckg,vmax,vmin,VOLTAGE,StatusVoltage,VM,ActivedV,...
MiuBus,ngn,ActivedP,PGEN,StatP,MiuGen,PMAX,PMIN);
%% NODES
for(ii= 1: nbb)
if(ActivedV(ii)˜=0)
Vol=VM(ii);
VOLD=VOLTAGE(ii);
VolMax=vmax(ii);
VolMin=vmin(ii);
StatV=StatusVoltage(ii);
MiuB=MiuBus(ii);
send= ii;
Hessian(nbb+send,nbb+send)=Hessian(nbb+send,nbb+send)+Ckv;
if(StatV==0)
grad(nbb+send)=grad(nbb+send)-(MiuB+Ckv*(Vol-VolMax));
end
386 APPENDIX C

if(StatV==2)
grad(nbb+send)=grad(nbb+send)-(MiuB+Ckv*(Vol-VolMin));
end
if(ActivedV(send)==0)
ActivedV(send)=1;
end
end
end


%% GENERATORS
for(ii= 1: ngn)
if(ActivedP(ii)˜=0)
for (ii = 1: ngn)
send=ii;
PGenI=PGEN(ii);
Stat=StatP(ii);
MiuG= MiuGen(ii);
Max_PGen=PMAX(ii);
Min_PGen=PMIN(ii);
Hessian(4*nbb+send,4*nbb+send) = Hessian(4*nbb+send,4*nbb...
+send)+Ckg;
if(Stat==0)
grad(4*nbb+send)=grad(4*nbb+send)-(MiuG+Ckg*...
(PGenI-Max_PGen));
end
if(Stat==2)
grad(4*nbb+send)=grad(4*nbb+send)-(MiuG+Ckg*...
(PGenI-Min_PGen));
end
if(ActivedP(ii)==0)
ActivedP(ii)=1;
end
end
end
end
return; %%End of AugmentedLagrangian function



function [Dz,VA,VM,LambdaP,LambdaQ,PGEN] = Actualisation(Dz,nbb,VA,VM,...
LambdaP,LambdaQ,bustype,ngn,PGEN);
for(ii= 1: nbb)
p1 = Dz(ii);
p2 = Dz(nbb+ii);
p3 = Dz(2*nbb+ii);
p4 = Dz(3*nbb+ii);
%% Load bus
387
COMPUTER PROGRAM FOR OPTIMAL POWER FLOW SOLUTIONS

if(bustype(ii)==3)
VA(ii)=VA(ii)+p1;
VM(ii)=VM(ii)+p2;
LambdaP(ii)=LambdaP(ii)+p3;
LambdaQ(ii)=LambdaQ(ii)+p4;
%%%% Slack bus
elseif(bustype(ii)==1)
VA(ii)=0.0; %%%% It must be cero
VM(ii)=VM(ii)+p2;
LambdaP(ii)=LambdaP(ii)+p3;
LambdaQ(ii)=LambdaQ(ii)+0.0; %%%% It must be cero
%%%% Generator bus
elseif(bustype(ii)==2)
VA(ii)=VA(ii)+p1;
VM(ii)=VM(ii)+p2;
LambdaP(ii)=LambdaP(ii)+p3;
LambdaQ(ii)=LambdaQ(ii)+0.0; %%%% It must be cero
end
end

for(ii= 1: ngn)
PGEN(ii)=PGEN(ii)+Dz(4*nbb+ii);
end
return; %%End of Actualisation function



function [Hessian,grad] = Reset_Hessian_grad(Hessian,grad,ngn,nbb);
for ii = 1: (4*nbb+ngn)
for jj = 1: (4*nbb+ngn)
Hessian(ii,jj)=0;
end
grad(ii)=0;
end
return; %%End of Reset_Hessian_grad function



function [Optimo] =ReviewNodes(Optimo,grad,nbb,bustype,ActivedV,TOL);
for (ii = 1: nbb)
p1=grad(ii);
p2=grad(nbb+ii);
p3=grad(2*nbb+ii);
p4=grad(3*nbb+ii);
if(ActivedV(ii)˜=0)
p2=1e-10;
else
p2=1e-10;
end
388 APPENDIX C

if(p1<0),p1=-p1;end
if(p2<0),p2=-p2;end
if(p3<0),p3=-p3;end
if(p4<0),p4=-p4;end


if(bustype(ii)==3) %%%%%% Load bus
if(p1>TOL),Optimo=1;end
if(p2>TOL),Optimo=1;end
if(p3>TOL),Optimo=1;end
if(p4>TOL),Optimo=1;end
end
if(bustype(ii)==1) %%%% Slack bus
if(p1>TOL),Optimo=1;end
if(p2>TOL),Optimo=1;end
if(p3>TOL),Optimo=1;end
if(p4>TOL),Optimo=1;end
end
if(bustype(ii)==2) %%%% Generator bus
if(p1>TOL),Optimo=1;end
if(p2>TOL),Optimo=1;end
if(p3>TOL),Optimo=1;end
if(p4>TOL),Optimo=1;end
end
end
return; %%End of ReviewNodes function



function [Optimo] = ReviewGen(Optimo,grad,ngn,ActivedP,TOL,nbb);
for (ii = 1: ngn)
p1=grad(4*nbb+ii);
if(ActivedP(ii)˜=0), p1=1e-10;end
if(p1<0),p1=-p1;end
if(p1>TOL)Optimo=1;end
end
return; %%End of ReviewGen function



function [MaxDz] = Norma(MaxDz,grad,nbb,ngn);
MaxDz=0;
for(ii = 1: nbb)
value=grad(ii);
if(value<0),value=-value;end
if(value>MaxDz),MaxDz=value;end
end
for(ii = (2*nbb+1): (4*nbb+ngn))
value=grad(ii);
389
COMPUTER PROGRAM FOR OPTIMAL POWER FLOW SOLUTIONS

if(value<0),value=-value;end
if(value>MaxDz),MaxDz=value;end
end
return; %%End of Norma function




function [StatusVoltage,SetVol,NumNode,ValueVoltage]...
= Check_Limits_V(nbb,TolVoltage,vmin,vmax,VM,StatusVoltage,Dz,...
SetVol,NumNode,ValueVoltage,ActivedV);
for (ii = 1: nbb)
if(ActivedV(ii)==0)
Movement=Dz(nbb+ii);
Kindex=0.0;
if(((vmin(ii)-TolVoltage)<VM(ii)) &(VM(ii)<(vmax(ii) + ...
TolVoltage)))
StatusVoltage(ii)=1;
end
if(VM(ii)>(vmax(ii)+TolVoltage))
StatusVoltage(ii)=0;
Kindex=(VM(ii)-vmax(ii))/Movement;
if(Kindex<0),Kindex=-Kindex;end
end
if(VM(ii)<(vmin(ii)-TolVoltage))
StatusVoltage(ii)=2;
Kindex=(vmin(ii)-VM(ii))/Movement;
if(Kindex<0),Kindex=-Kindex;end
end
if(Kindex > 0)
SetVol=SetVol+1;
NumNode(SetVol)= ii;
ValueVoltage(SetVol)=Kindex;
end
end
end
return; %%End of Check_Limits_V function




function [StatP,SetGenP,NumGenerator,ValueGenerator]...
= Check_Limits_P(ngn,TolPower,PMIN,PMAX,PGEN,StatP,SetGenP,...
NumGenerator,ValueGenerator,ActivedP);
for (ii = 1: ngn)
if(ActivedP(ii)==0)
Movement=1;
Kindex=0.0;
390 APPENDIX C

if(((PMIN(ii)-TolPower)<PGEN(ii)) &(PGEN(ii)<(PMAX(ii)...
+TolPower)))
StatP(ii)=1;
end
if(PGEN(ii)>(PMAX(ii)+TolPower))
StatP(ii)=0;
Kindex=(PGEN(ii)-PMAX(ii))/Movement;
if(Kindex<0),Kindex=-Kindex;end
end
if(PGEN(ii)<(PMIN(ii)-TolPower))
StatP(ii)=2;
Kindex=(PMIN(ii)-PGEN(ii))/Movement;
if(Kindex<0),Kindex=-Kindex;end
end
if(Kindex > 0)
SetGenP=SetGenP+1;
NumGenerator(SetGenP)= ii;
ValueGenerator(SetGenP)=Kindex;
end
end
end
return; %%End of Check_Limits_P function




function [Hessian,grad,ActivedV] = AugmentedLagrangianV(nbb,...
Hessian,grad,Ckv,vmax,vmin,StatusVoltage,VM,ActivedV,MiuBus,...
SetVol,NumNode,ValueVoltage,MaxVol,EnforceTol);
for(jj = 1: SetVol)
in1=NumNode(jj);
if(((ValueVoltage(jj)/MaxVol)>=EnforceTol)|(VM(in1)<0.80)...
|(VM(in1)>1.20))
ii=in1;
Vol=VM(ii);
VolMax=vmax(ii);
VolMin=vmin(ii);
StatV=StatusVoltage(ii);
MiuB=MiuBus(ii);
send= ii;
Hessian(nbb+send,nbb+send)=Hessian(nbb+send,nbb+send)+Ckv;
if(StatV==0)
grad(nbb+send)=grad(nbb+send)-(MiuB+Ckv*(Vol-VolMax));
end
if(StatV==2)
grad(nbb+send)=grad(nbb+send)-(MiuB+Ckv*(Vol-VolMin));
end
if(ActivedV(send)==0)
391
COMPUTER PROGRAM FOR OPTIMAL POWER FLOW SOLUTIONS

ActivedV(send)=1;
end
end
end
return; %%End of AugmentedLagrangianV function


function[Hessian,grad,ActivedP] = AugmentedLagrangianG(ngn,nbb,...
Hessian,grad,Ckg,PMIN,PMAX,StatP,PGEN,ActivedP,MiuGen,SetGenP,...
NumGenerator,ValueGenerator,MaxGen,EnforceTol);
for(jj = 1: SetGenP)
in1=NumGenerator(jj);
if((ValueGenerator(jj)/MaxGen)>=EnforceTol)
ii=in1;
send=ii;
PGenI=PGEN(ii);
Stat=StatP(ii);
MiuG= MiuGen(ii);
Max_PGen=PMAX(ii);
Min_PGen=PMIN(ii);
Hessian(4*nbb+send,4*nbb+send)=Hessian(4*nbb+send,4*nbb+send)...
+Ckg;
if(Stat==0)
grad(4*nbb+send)=grad(4*nbb+send)-(MiuG+Ckg*(PGenI- ...
Max_Pgen));
end
if(Stat==2)
grad(4*nbb+send)=grad(4*nbb+send)-(MiuG+Ckg*(PGenI-...
Min_PGen));
end
if(ActivedP(ii)==0)
ActivedP(ii)=1;
end
end
end
return; %%End of AugmentedLagrangianG function


function [Hessian,grad,ActivedV] = AugmentedLagrangian_IV(nbb,...
Hessian,grad,Ckv,vmax,vmin,StatusVoltage,VM,ActivedV,MiuBus,ii);
Vol=VM(ii);
VolMax=vmax(ii);
VolMin=vmin(ii);
StatV=StatusVoltage(ii);
MiuB=MiuBus(ii);
send= ii;
Hessian(nbb+send,nbb+send)=Hessian(nbb+send,nbb+send)+Ckv;
if(StatV==0)
392 APPENDIX C

grad(nbb+send)=grad(nbb+send)-(MiuB+Ckv*(Vol-VolMax));
end
if(StatV==2)
grad(nbb+send)=grad(nbb+send)-(MiuB+Ckv*(Vol-VolMin));
end
if(ActivedV(send)==0);ActivedV(send)=1;end
return; %%End of AugmentedLagrangian_IV function


function [MiuBus,NumMiu_V,ValMiu_V,RelVol] = IdentifyMiuBus(...
StatusVoltage,MiuBus,Ckv,VM,vmax,vmin,ActivedV,RelVol,NumMiu_V,...
ValMiu_V,ii);
Tolerance=0.0001;
StatV=StatusVoltage(ii);
if(MiuBus(ii)==0),Tolerance=0;end
temp=MiuBus(ii);
if(StatV==0)
MiuBus(ii)=MiuBus(ii)+Ckv*(VM(ii)-vmax(ii)+Tolerance);
end
if(StatV==2)
MiuBusS(ii)=MiuBus(ii)+Ckv*(VM(ii)-vmin(ii)-Tolerance);
end
Kindex=0.0;
Kindex1=0.0;
if(StatV==1),Kindex1=Kindex; end
if(ActivedV(ii)==2)
if(StatV==0)
if(MiuBus(ii)<0),Kindex1=MiuBus(ii); end
if(MiuBus(ii)>=0),Kindex1=0; end
if(Kindex1<0),Kindex1=-Kindex1; end
end
if(StatV==2)
if(MiuBus(ii)>0),Kindex1=MiuBus(ii); end
if(MiuBus(ii)<=0),Kindex1=0; end
if(Kindex1<0),Kindex1=-Kindex1; end
end
end
MiuBus(ii)=temp;
if(Kindex1>0)
RelVol=RelVol+1;
NumMiu_V(RelVol)=ii;
ValMiu_V(RelVol)=Kindex;
end
return; %%End of IdentifyMiuBus function


function [Hessian,grad,ActivedP] = AugmentedLagrangian_IG(ngn,...
Hessian,grad,Ckg,PMIN,PMAX,StatP,PGEN,ActivedP,MiuGen,ii);
393
COMPUTER PROGRAM FOR OPTIMAL POWER FLOW SOLUTIONS

send=ii;
PGenI=PGEN(ii);
Stat=StatP(ii);
MiuG= MiuGen(ii);
Max_PGen=PMAX(ii);
Min_PGen=PMIN(ii);

Hessian(4*nbb+send,4*nbb+send)=Hessian(4*nbb+send,4*nbb+send)+Ckg;
if(Stat==0)
grad(4*nbb+send)=grad(4*nbb+send)-(MiuG+Ckg*(PGenI-Max_PGen));
end
if(Stat==2)
grad(4*nbb+send)=grad(4*nbb+send)-(MiuG+Ckg*(PGenI-Min_PGen));
end

if(ActivedP(ii)==0);ActivedP(ii)=1;end
return; %%End of AugmentedLagrangian_IG function



function [MiuBus,NumMiu_P,ValMiu_P,RelGen] = IdentifyMiuGen...
(StatP,MiuGen,Ckg,PGEN,PMAX,PMIN,ActivedP,RelGen,NumMiu_P,...
ValMiu_P,ii);
Tolerance=0.0000001;
Stat=StatP(ii);
if(MiuGen(ii)==0),Tolerance=0;end
temp=MiuGen(ii);
Kindex=0.0;
Kindex1=0.0;
if(Stat==0),MiuGen(ii)=MiuGen(ii)+Ckg*(PGen+Tolerance-Max_PGen);end
if(Stat==2),MiuGen(ii)=MiuGen(ii)+Ckg*(PGen-Tolerance-Min_PGen);end
if(Stat==1),Kindex1=Kindex;end
if(ActivedP(ii)==2)
if(Stat==0)
if(MiuGen(ii)<0),Kindex1=MiuGen(ii);end
if(MiuGen(ii)>0),Kindex1=0;end
if(Kindex1<0),Kindex1=-Kindex1;end
end
if(Stat==2)
if(MiuGen(ii)>0),Kindex1=MiuGen(ii);end
if(MiuGen(ii)<0),Kindex1=0;end
if(Kindex1<0),Kindex1=-Kindex1;end
end
end
MiuGen(ii)=temp;
if(Kindex1>0)
RelGen=RelGen+1;
NumMiu_P(RelVol)=ii;
394 APPENDIX C

ValMiu_P(RelVol)=Kindex;
end
return; %%End of IdentifyMiuGen function



function [Hessian,grad,ActivedV,StatusVoltage] =...
ReleasingAumentedLagrangianV(VM,vmax,vmin,StatusVoltage,MiuBus,...
Hessian,grad,Ckv,ActivedV,StatusVoltage,send);
Vol=VM(send);
VolMax=vmax(send);
VolMin=vmin(send);
StatV=StatusVoltage(send);
MiuB=MiuBus(send);
Hessian(nbb+send,nbb+send)=Hessian(nbb+send,nbb+send)-Ckv;
if(StatV==0),grad(nbb+send)=grad(nbb+send)+MiuB+Ckv*(Vol-VolMax);end
if(StatV==2),grad(nbb+send)=grad(nbb+send)+MiuB+Ckv*(Vol-VolMin);end
ActivedV(send)=0;
StatusVoltage(send)=1;
return; %%End of ReleasingAumentedLagrangianV function




function [Hessian,grad,ActivedP,StatP]=ReleasingAumentedLagrangianG...
(PGEN,PMAX,PMIN,StatP,MiuGen,Hessian, grad,Ckg,ActivedP,StatP,send);
PGenI=PGEN(send);
Stat=StatP(send);
MiuG=MiuGen(send);
Hessian(4*nbb+send,4*nbb+send)=Hessian(4*nbb+send,4*nbb+send)-Ck;
if(Stat==0)
grad(4*nbb+send)=grad(4*nbb+send)+MiuG+Ckg*(PGenI-PMAX(send));
end
if(Stat==2),
grad(4*nbb+send)=grad(4*nbb+send)+MiuG+Ckg*(PGenI-PMIN(send));
end
ActivedP(send)=0;
StatP(send)=1;
return; %%End of ReleasingAumentedLagrangianG function



function [MiuBus] = MultiplierBus(nbb,Ckv,MiuBus,VM,vmax,vmin);
TolVoltage=0.00000001;
for(ii = 1:nbb)
if((MiuBus(ii)+Ckv*(VM(ii)-vmax(ii)+TolVoltage)) >= 0)
MiuBus(ii)=MiuBus(ii)+Ckv*(VM(ii)-vmax(ii)+TolVoltage);
elseif((MiuBus(ii)+Ckv*(VM(ii)-vmin(ii)-TolVoltage)) <= 0)
MiuBus(ii)=MiuBus(ii)+Ckv*(VM(ii)-vmin(ii)-TolVoltage);
395
COMPUTER PROGRAM FOR OPTIMAL POWER FLOW SOLUTIONS

else
MiuBus(ii)=0;
end
end
return; %%End of MultiplierBus function



function [MiuGen] = MultiplierGen(ngn,Ckg,MiuGen,PGEN,PMAX,PMIN);
TolPower=0.0000001;
for(ii = 1:ngn)
if((MiuGen(ii)+Ckg*(PGEN(ii)-PMAX(ii)+TolPower)) >= 0)
MiuGen(ii)=MiuGen(ii)+Ckg*(PGEN(ii)-PMAX(ii)+TolPower);
elseif((MiuGen(ii)+Ckg*(PGEN(ii)-PMIN(ii)-TolPower)) <= 0)
MiuGen(ii)=MiuGen(ii)+Ckg*(PGEN(ii)-PMIN(ii)-TolPower);
else
MiuGen(ii)=0;
end
end
return; %%End of MultiplierGen function



function [ChangeStat,Kindex1,statusgen,LambdaQ,ActivedQ,bustype] =...
IdentifyConstQ(ChangeStat,genbus,statusgen,LambdaQ,ActivedQ,...
bustype, ii);
send=genbus(ii);
ChangeStat=statusgen(ii);
Kindex1=0;
if(statusgen(ii)==0)
if(LambdaQ(send)<0)
ActivedQ(ii)=2;
Kindex1=0;
bustype(send)=2;
statusgen(ii)=1;
LambdaQ(send)=0;
else
ActivedQ(ii)=1;
Kindex1=1;
end
end
if(statusgen(ii)==2)
if(LambdaQ(send)>0)
ActivedQ(ii)=2;
Kindex1=0;
bustype(send)=2;
statusgen(ii)=1;
LambdaQ(ii)=0;
396 APPENDIX C

else
ActivedQ(ii)=1;
Kindex1=1;
end
end
if(ChangeStat==statusgen(ii))
ChangeStat=0;
else
ChangeStat=1;
end
return; %%End of IdentifyConstQ function



function [GenViolado] = CheckQGenLimits(GenViolado,QLOAD,Qbus,...
QMIN,QMAX,genbus,TolPower,ii);
GenViolado=1;
Qpower=QLOAD(genbus(ii))+Qbus(genbus(ii));
if(((QMIN(genbus(ii))-TolPower)<Qpower)(QMAX(genbus(ii))+TolPower)))
GenViolado=0;
end
return; %%End of CheckQGenLimits function



function [in1,in2] = MatrixWVoltageMiu(in1,in2,TolVoltage,VM,vmax,...
vmin,MiuBus,ActivedV,StatusVoltage,ii);
if(((vmin(ii)-TolVoltage)<VM(ii))&(VM(ii)<(vmax(ii)+TolVoltage)))
in1=0;
else
in1=1;
end
if((VM(ii)<0.5)|(VM(ii)>1.5))
fprintf(™\n************************************************™);
fprintf(™\n* *™);
fprintf(™\n* UNFEASIBLE SOLUTION *™);
fprintf(™\n* *™);
fprintf(™\n************************************************™);
end
delta=0;
if(ActivedV(ii)˜=0)
if(StatusVoltage(ii)==0)
delta=VM(ii)-vmax(ii);
if(delta<0),delta=-delta; end
if(delta>0.0001),in2=1; end
end
if(StatusVoltage(ii)==2)
delta=VM(ii)-vmax(ii);
397
COMPUTER PROGRAM FOR OPTIMAL POWER FLOW SOLUTIONS

if(delta<0),delta=-delta; end
if(delta>0.0001),in2=1; end
end
end
if(ActivedV(ii)==0)
if(MiuBus(ii)>0),in2=1;end
if(MiuBus(ii)<0),in2=1;end
end
return; %%End of MatrixWVoltageMiu function


function [in3,in4] = MatrixWGenMiu(in3,in4,TolPower,PGEN,PMAX,PMIN,...
MiuGen,ActivedP,StatP,ii);
if(((PMIN(ii)-TolPower)<PGEN(ii))(PMAX(ii)+TolPower)))
in3=0;
else
in3=1;
end
delta=0;
if(ActivedP(ii)˜=0)
if(StatP(ii)==0)
delta=PGEN(ii)-PMAX(ii);
if(delta<0),delta=-delta; end
if(delta>0.0001),in4=1; end
end
if(StatP(ii)==2)
delta=PGEN(ii)-PMAX(ii);
if(delta<0),delta=-delta; end
if(delta>0.0001),in4=1; end
end
end
if(ActivedP(ii)==0)
if(MiuGen(ii)>0),in4=1;end
if(MiuGen(ii)<0),in4=1;end
end
return; %%End of MatrixWGenMiu function


function [sum] = CheckConvexi¬cation(nbb,Cc,VM,VOLTAGE);
sum=0;
for( ii = 1: nbb)
sum =sum+Cc*(VM(ii)-VOLTAGE(ii))*(VM(ii)-VOLTAGE(ii));
end
fprintf(™\n=== Convexi¬cation Value ===™);
fprintf(™\n C = %%12.8f ™, sum);
fprintf(™\n™);

return; %%End of CheckConvexi¬cation function
Index

Convergence, 6, 48, 94, 98, 101“2, 111, 129“
a axis, 82
31, 141, 415“6, 150, 154, 166, 170, 178,
d axis, 82, 84
188, 198, 213, 253, 234, 248, 251“2, 262,
q axis, 83“4
265, 276, 279, 280, 284, 287, 288“9, 291,
295, 299, 301, 306, 308
ABCD parameters, 62“3
Converters, 4“5
Active power, 3“4, 18, 28, 34“8, 40, 94, 96“7,
Current source converters, 29
111, 118“9, 129, 131“3, 134, 140“4, 147,
Custom power, 9“11
150, 171“3, 178, 181“2, 187, 199“201, 215,
218, 225, 227, 237, 253, 255“7, 262, 264,
Damper winding, 82, 84, 86
266, 268“70, 272, 275, 278“80, 282,
DC-AC converters, 28“9, 231, 301
286“97, 299, 301“2, 304, 306“7, 315,
DC capacitor, 5, 34“5,
317, 318, 320, 322“3, 328, 335“7, 339
Delta-connected TCR, 14, 17, 249“50
Advanced series capacitor (ASC), 42
Delta-delta, 70, 76, 78“81
Amplitude modulation ratio, 32
Dependent variables, 147, 268“9
Angular stability, 3
Deregulation, 4, 93, 311, 331, 339
Annular conductors, 47
Diode, 30
Distribution, 1“2, 4, 9“11, 43“4, 47, 129, 153,
Bank of capacitors, 16“8, 231“49
231, 267, 290, 307
Bank of transformers, 89, 236
Double circuit transmission line, 54“5
Bulk power load, 87
Dynamic stability, 2“3
Bulk power supply, 28, 56
Bundle conductors, 56, 59
Economic dispatch , 267, 278“80, 299,
305
Carrier frequency, 30 Electrical length, 3, 10“1, 18, 20, 45
Carrier signal, 30, 32 Electricity supply industry, 1, 4, 41, 339
Compensated transmission line, 228, 264“5, 297 Electromagnetic transients, 2“3
Compensation, 3, 29, 33, 70, 159, 178, 232, 249, EMF, 82, 86
265, 299, 336 Equality constraints, 268“72, 275, 289, 292,
Series, 4, 40, 43, 60“1, 63, 89, 178, 229, 319 296, 301
Shunt, 5, 11, 44, 61, 63, 89, 104, 319 Equivalent p-circuit, 35, 44
Complex depth, 49
Conduction angle, 12“3, 20 FACTS
Controllable reactance, 41 Concepts, 2, 9
Controllable susceptance, 5, 12, 149 Controllers, 2“6, 9, 11, 150, 153“4, 227“9,
Control variables, 97, 301, 306, 268“70, 278 232, 249, 265, 307“8, 329

FACTS: Modelling and Simulation in Power Networks.
´ ´
Enrique Acha, Claudio R. Fuerte-Esquivel, Hugo Ambriz-Perez and Cesar Angeles-Camacho
# 2004 John Wiley & Sons, Ltd ISBN: 0-470-85271-2
400 INDEX

High voltage, 2, 4, 5, 9, 11, 14, 16, 29, 30,
Modelling, 6, 227, 232, 265, 308
38“9, 43, 54, 111, 154, 156, 166, 216, 219,
Series controllers, 143, 173
225“6
Simulation, 6
High-voltage transmission, 44, 54, 62, 70, 76,
Technology, 2, 4, 6
89, 111“2
Fast decoupled, 111“3, 115, 117, 150
HVDC, 5, 9“11, 29, 38, 40“1, 43, 154, 216“9,
Feasibility region, 140, 143
225“7, 315
Field winding, 82“4
Filters, 11, 14
IGBT, 28“30
Firing angle, 12“4, 16, 20, 22“3, 26, 155,
Induction motor, 87
158“9, 162“3, 166“7, 169“71, 180, 182,
Inductive effects, 45“6
187“8, 190“1, 227, 229, 249, 251“3,
Inequality constraints, 268, 270“2, 275, 277,
255“7, 291“7, 299
279“80, 284, 287“8, 292, 296, 298, 301,
First derivative, 293, 299, 304
304
Fluctuating loads, 9
Interconnection, 1, 9, 267
Frame of reference, 15, 18, 19, 28, 41, 44, 61,
Internal impedances, 47“8
67, 79, 80, 83, 89, 90, 102, 154“5, 229, 232,
IPC, 4, 5
236“7, 249, 283, 307
Iron core, 70“1
Frequency dependency, 45“8, 50, 52
Iteration, 19, 94, 98“102, 111“2, 117, 120, 122,
Frequency modulation ratio, 32
129“31, 133, 135, 141“2, 144, 146,
Frequency switching, 28, 30, 33
149“50, 154, 158, 162, 167, 170, 178, 182,
Fundamental frequency, 11“5, 19“20, 24“6,
188“9, 198, 213, 251“6, 218, 234, 246“7,
30“4, 40“1, 44“5, 50, 60“1, 70, 89, 180,
251“2, 255, 257, 262, 271, 276“81, 285“6,
189, 200, 216, 253“4, 256
290“1, 295“6, 299, 300, 306“8
Gauss-Seidel, 98
Jacobian, 98“100, 102“4, 111“2, 117, 120,
Generation, 3“4, 35, 43, 94, 96, 102, 118, 147,
122, 132“4, 144“5, 148, 150, 154“5, 157,
153, 199, 215, 232, 261, 268“70, 272, 275,
167, 173, 182, 191“2, 202, 218, 234“5,
279“82, 286“7, 290, 292“6, 299, 301,
251“2, 260, 262, 271“2
306“7, 312“3, 322, 324, 331, 336, 340
Generator PQ bus, 97, 102
Kayenta substation, 20, 24, 26, 189
Generator PV bus, 97, 101“3, 121
Kron™s reduction, 52, 79, 80
Geometric imbalances, 44, 46, 50, 60, 70, 89,
Kuhn and Tucker, 276, 283
231“2, 236, 252, 264
Geometric mean radius, 47, 52
Lagrange multipliers, 270“4, 280“2, 284“5,
Global iteration, 276“7
287, 289, 292“5, 297“9, 301“2, 305“7
Global optimal, 268, 272
Lagrangian, 270“5, 279, 283“4, 286“8, 291“3,
Gradient, 272“3, 276, 279, 283“4, 288, 293,
296“8, 301“2
304
Linear transformations, 15
Ground return, 45“6, 48“50, 60, 70, 89
Line impedance control, 1, 5“6, 52
Ground wires, 44, 50“5, 59, 79, 89
Line transpose, 62
GTO, 5, 28“9, 41
Load angle, 90, 236
Load ¬‚ow analysis, 6, 82, 157“8, 166
Harmonic, 12, 14, 19, 20, 24, 30, 32“3, 41, 70,
Load PQ bus, 97, 268, 270, 277
87, 232, 239
Load tap changer (LTC), 2, 4, 72, 119“20, 122,
Harmonic cancellation, 14, 41
127, 130“1, 145“9, 154, 158, 283“7
Harmonic domain, 19
Loads, 1, 3, 5, 9, 11, 43“4, 56, 70, 86“90, 94,
Harmonic voltages, 9, 31“2
97, 102, 156, 231, 236, 246“7, 265, 267“9,
Hessian, 271“2
274, 297, 301, 313, 315“7, 319“20, 322“4,
High current, 2, 4, 9
330“1, 335“6
High power, 9, 29, 30, 232
Long-distance, 1, 43, 60“1, 67, 231
High-speed control, 2, 6, 10
401
INDEX

Penalty function, 270, 275, 277, 280, 284, 288,
Long-line, 46, 60“1, 66, 89, 236, 319
293, 299, 304
Low-voltage distribution, 9, 11, 44, 231
Penalty weighting factors, 276“80
Lumped parameters, 45, 61, 66
Per-unit system, 55, 72
Phase angle regulator, 2, 4, 35
Magnetic circuit, 70
Phase coordinates, 15“6, 18, 35, 38, 76, 231“2,
Magnetising branch, 73, 75, 119, 121, 132,
236, 253, 261“2, 264“5
134
Phase domain, 11, 15, 28, 41, 257
Magnetising current, 70
Matlab1, 56, 59, 61, 63, 66, 69, 89, 104, Phase-shifting transformers, 3, 4, 35, 73, 75,
81, 89, 112, 119, 132“5, 140“1, 143, 145,
112“3, 115, 122, 135, 150, 159, 162“3,
150, 154, 200, 268, 283, 286“7,
173, 182, 192, 203, 218“9, 236“7,
289“91
265, 281
Power converters, 4, 29, 231, 301
Multilimb transformers, 76
Power electronics, 1, 2, 9, 10, 157, 232, 335
Multiplier method, 275“6, 280, 282, 285,
Power ¬‚ow, 1, 3, 6, 11, 33“5, 40“1, 44, 56,
308
70, 75, 90, 93“5, 97“9, 101, 105, 111“3,
Mutual potential coef¬cients, 47
115, 117, 119“22, 129, 132“4, 135,
140“3, 149“50, 153“4, 157, 159, 162“3,
Net active power, 96“7, 129
167, 171, 173, 181“2, 187, 192, 198“9,
Net reactive power, 96“7, 118, 129
203, 215, 217“9, 225, 229, 231“4, 236,
Network, 1“3, 5“6, 9“11, 13“4, 16, 20“1, 35,
244, 249, 250, 255, 257, 261“2, 264“5,
37, 43“4, 55“6, 70, 93“4, 97, 101“2, 104,
267“70, 272“3, 277“82, 285“90, 294,
112, 115, 118, 127, 129“30, 140, 145, 150,
296“7, 301, 305“7, 311“3, 315“7,
153“5, 157, 169“70, 179, 187, 198“9, 201,
319, 321“3, 325, 329, 331, 335“7,
213, 215, 217“8, 225, 227“9, 244“7,
340
252“3, 257, 262, 264“5, 267“9, 274,
Analysis, 6, 19, 26, 72, 87, 93, 98, 111, 145,
282“3, 285, 290“2, 294“5, 297, 299,
153, 155, 156, 191, 227“8, 236, 257, 283
306“8, 311“12, 315“16, 322, 234“7,
Control, 2, 3, 18, 34, 37, 41, 140, 143“4, 200,
329, 331, 334“7, 339“40
213, 225, 257, 264, 297, 299
Newton“Raphson, 1, 6, 98“9, 101“2, 104,
Loops, 1
111“2, 115, 117, 122, 129, 135, 145, 150,
Model, 119, 132, 154“5, 171, 178, 180, 189,
153“5, 157, 159, 162“3, 166, 171“3, 180,
191, 201, 228, 249, 253, 255, 261
182, 192, 203, 216, 218“9, 227“8, 232,
Overvoltage, 4, 274
234, 236“7, 253, 264“5, 281, 335
Power interruptions, 9
Newton™s method, 158, 268, 270“2, 275, 278,
Powers Mismatch, 94, 96, 99, 101“3, 111, 129,
280, 282“7, 291“3, 296“7, 299, 301“2,
133, 141, 146, 150, 159, 163, 173, 178, 182,
306“8
188“9, 202, 213, 215“6, 218, 234, 246“7,
Nodal admittance, 15“6, 47, 75“6, 78, 112, 135,

<<

. 16
( 17)



>>