facat
/
threephasese-huber
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
threephasese-huber/OPF.m

353 lines
12 KiB
Mathematica
Raw Permalink Normal View History

第一次提交 Signed-off-by: dmy@lab <dmy@lab.lab>
2015-06-02 17:24:13 +08:00
function [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,AME_mVolt,AME_mPD,AME_mQD,isConverged,totalTime]=OPF()
tic
clc
clear
完成了Huber方法。 Signed-off-by: dmy@lab <dmy@lab.lab>
2015-06-03 10:35:15 +08:00
lineZ=readLineZ('./../DistributionNetwork-Power2Current/modified-feeder69/lineParameter.txt');
第一次提交 Signed-off-by: dmy@lab <dmy@lab.lab>
2015-06-02 17:24:13 +08:00
[ fsY0, fsY1, fsY2,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP ...
phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY ...
完成了Huber方法。 Signed-off-by: dmy@lab <dmy@lab.lab>
2015-06-03 10:35:15 +08:00
cap]=dataRead(lineZ,'./../DistributionNetwork-Power2Current/modified-feeder69/data.txt');
第一次提交 Signed-off-by: dmy@lab <dmy@lab.lab>
2015-06-02 17:24:13 +08:00
% phaseASpotLoadP(phaseASpotLoadP==0)=0.002;
% phaseBSpotLoadP(phaseBSpotLoadP==0)=0.002;
% phaseCSpotLoadP(phaseCSpotLoadP==0)=0.002;
% phaseASpotLoadQ(phaseASpotLoadQ==0)=0.002;
% phaseBSpotLoadQ(phaseBSpotLoadQ==0)=0.002;
% phaseCSpotLoadQ(phaseCSpotLoadQ==0)=0.002;
%<EFBFBD><EFBFBD><EFBFBD>ɲ<EFBFBD>ƽ<EFBFBD><EFBFBD>
% for I=1:length(phaseASpotLoadP)
% roll=mod(round(rand()*10),3);
% if roll==0
% phaseBSpotLoadP(I)=phaseBSpotLoadP(I)+phaseASpotLoadP(I)*.45;
% phaseCSpotLoadP(I)=phaseCSpotLoadP(I)+phaseASpotLoadP(I)*.45;
% phaseASpotLoadP(I)=phaseASpotLoadP(I)*.10;
%
% phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)+phaseASpotLoadQ(I)*.45;
% phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)+phaseASpotLoadQ(I)*.45;
% phaseASpotLoadQ(I)=phaseASpotLoadQ(I)*.10;
% end
% if roll==1
% phaseASpotLoadP(I)=phaseASpotLoadP(I)+phaseBSpotLoadP(I)*.45;
% phaseCSpotLoadP(I)=phaseCSpotLoadP(I)+phaseBSpotLoadP(I)*.45;
% phaseBSpotLoadP(I)=phaseBSpotLoadP(I)*.10;
%
% phaseASpotLoadQ(I)=phaseASpotLoadQ(I)+phaseBSpotLoadQ(I)*.45;
% phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)+phaseBSpotLoadQ(I)*.45;
% phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)*.10;
% end
% if roll==2
% phaseASpotLoadP(I)=phaseASpotLoadP(I)+phaseCSpotLoadP(I)*.45;
% phaseBSpotLoadP(I)=phaseBSpotLoadP(I)+phaseCSpotLoadP(I)*.45;
% phaseCSpotLoadP(I)=phaseCSpotLoadP(I)*.10;
%
% phaseASpotLoadQ(I)=phaseASpotLoadQ(I)+phaseCSpotLoadQ(I)*.45;
% phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)+phaseCSpotLoadQ(I)*.45;
% phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)*.10;
% end
% end
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>begin
a=exp(1j*2*pi/3);
Tp2f=1/3*[1 1 1;
1 a a^2;
1 a^2 a];
Tp2f=sparse(Tp2f);
Tf2p=inv(Tp2f);
fsY1amp=abs(fsY1);
[r,c,fsY1ang]=find(fsY1);
fsY1ang=angle(fsY1ang);
Pabc=phaseASpotLoadP+phaseBSpotLoadP+phaseCSpotLoadP;
Qabc=phaseASpotLoadQ+phaseBSpotLoadQ+phaseCSpotLoadQ;
busNum=length(phaseASpotLoadP);
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѹ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ
Vmf1=sparse(ones(busNum,1));
Vaf1=sparse(zeros(busNum,1));
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
PQi=nodeNum;
PG=sparse(busNum,1);
QG=sparse(busNum,1);
QGi=[Balance];
PD=Pabc/3;
QD=Qabc/3;
Loadi=find(PD~=0);
maxD=100000;% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
EPS=1e-6;
k=0;
kmax=20;
fsY11=fsY1;
fsY00=fsY0;
fsY22=fsY2;
Vf2=sparse(busNum,1);
If2=sparse(busNum,1);
Vf0=sparse(busNum,1);
If0=sparse(busNum,1);
%׼<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%ƽ<EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>0<EFBFBD><EFBFBD>1
fsY2(Balance,:)=0;
fsY2(:,Balance)=0;
fsY2=fsY2+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
%ƽ<EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>0<EFBFBD><EFBFBD>1
fsY0(Balance,:)=0;
fsY0(:,Balance)=0;
fsY0=fsY0+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
%%LU<EFBFBD>ֽ<EFBFBD>
[fsY0L,fsY0U,fsY0P,fsY0Q,fsY0R]=lu(fsY0);
[fsY2L,fsY2U,fsY2P,fsY2Q,fsY2R]=lu(fsY2);
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
tic
VoltpA=sparse(ones(busNum,1));
VoltpB=sparse(ones(busNum,1)).*exp(1j*-120/180*pi);
VoltpC=sparse(ones(busNum,1)).*exp(1j*+120/180*pi);
while(k<=kmax+10 && maxD> EPS)
k=k+1;
%<EFBFBD>Ѳ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ݿ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
SA=VoltpA.*conj(VoltpA.*sparse(cap.capNode,1,1j*cap.capB(:,1),busNum,1));
SB=VoltpB.*conj(VoltpB.*sparse(cap.capNode,1,1j*cap.capB(:,2),busNum,1));
SC=VoltpC.*conj(VoltpC.*sparse(cap.capNode,1,1j*cap.capB(:,3),busNum,1));
%<EFBFBD>Ȳ<EFBFBD>Ҫ<EFBFBD><EFBFBD><EFBFBD>ݣ<EFBFBD><EFBFBD>Ѿ<EFBFBD><EFBFBD>ӵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĵ<EFBFBD><EFBFBD>ɾ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ˡ<EFBFBD>20150405 By <EFBFBD><EFBFBD><EFBFBD><EFBFBD>Զ
SA=0;
SB=0;
SC=0;
iterPD=PD+real(SA+SB+SC)/3;
iterQD=QD+imag(SA+SB+SC)/3;
iterPhaseASpotLoadP=phaseASpotLoadP+real(SA);
iterPhaseBSpotLoadP=phaseBSpotLoadP+real(SB);
iterPhaseCSpotLoadP=phaseCSpotLoadP+real(SC);
iterPhaseASpotLoadQ=phaseASpotLoadQ+imag(SA);
iterPhaseBSpotLoadQ=phaseBSpotLoadQ+imag(SB);
iterPhaseCSpotLoadQ=phaseCSpotLoadQ+imag(SC);
[dP, dQ, YdotSinVolt, YdotCosVolt, diag_Volt_YdotSin, diag_Volt_YdotCos]=Unbalance(Balance,busNum, ...
PQi,PG,QG,QGi,iterPD,iterQD,Vmf1,Vaf1,fsY1amp,fsY1ang,r,c,Vf2,If2,Vf0,If0);%<EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
maxD=max(abs([dP;dQ;]));
jaco=Jacobi(Balance,busNum,QGi,Vmf1,YdotSinVolt,YdotCosVolt,diag_Volt_YdotSin,diag_Volt_YdotCos);%<EFBFBD>ſ˱Ⱦ<EFBFBD><EFBFBD><EFBFBD>
[dV, dVangle]=Solv(busNum,jaco,dP,dQ);%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
[Vmf1, Vaf1]=Modify(Vmf1,Vaf1,dV,dVangle,1);
fprintf('<EFBFBD><EFBFBD> %d <EFBFBD>ε<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϊ %f\n',k,full(maxD));
%ת<EFBFBD><EFBFBD>Ϊ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѹ
VoltpABC=Tp2f\conj([ Vf0'; (Vmf1.*exp(1j*Vaf1))'; Vf2']);%<EFBFBD><EFBFBD>Tp2f\ <EFBFBD><EFBFBD><EFBFBD><EFBFBD>Tf2p*
VoltpA=conj(VoltpABC(1,:)');
CurpA=-conj((iterPhaseASpotLoadP+1j*iterPhaseASpotLoadQ)./VoltpA);
VoltpB=conj(VoltpABC(2,:)');
CurpB=-conj((iterPhaseBSpotLoadP+1j*iterPhaseBSpotLoadQ)./VoltpB);
VoltpC=conj(VoltpABC(3,:)');
CurpC=-conj((iterPhaseCSpotLoadP+1j*iterPhaseCSpotLoadQ)./VoltpC);
f012=Tp2f*conj([CurpA';CurpB';CurpC']);
If0=conj(f012(1,:)');
If1=conj(f012(2,:)');
If2=conj(f012(3,:)');
If0(Balance)=0;
If2(Balance)=0;
%Vf0=fsY0\If0;
Vf0=fsY0Q*(fsY0U\(fsY0L\(fsY0P*(fsY0R\If0))));
%Vf2=fsY2\If2;
Vf2=fsY2Q*(fsY2U\(fsY2L\(fsY2P*(fsY2R\If2))));
fprintf('<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD>%f\n',toc);
%
end
FortiscueToc=toc;
fprintf('Fortiscue<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD> %f\n',FortiscueToc);
Vf1=Vmf1.*exp(1j*Vaf1);
(Vf0.*conj(fsY00*Vf0)+Vf1.*conj(fsY11*Vf1)+Vf2.*conj(fsY22*Vf2))*3;%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ݵĹ<EFBFBD><EFBFBD><EFBFBD>
conj(Tf2p*[If0(2);If1(2);If2(2)]).*(Tf2p*[Vf0(2);Vf1(2);Vf2(2)]);
IpABC=Tf2p*conj([If0';If1';If2']);
%ת<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѹ
VoltpABC=Tf2p*conj([ Vf0'; Vf1'; Vf2']);
disp([' A B C'])
full(abs(VoltpABC'))
fprintf('<EFBFBD>ڵ<EFBFBD><EFBFBD>Ŷ<EFBFBD>Ӧ\n');
disp([setIJ,nodeNum ])
%%<EFBFBD><EFBFBD><EFBFBD>鷴<EFBFBD>ƻ<EFBFBD>ȥ<EFBFBD>Ĺ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ƿ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
ub=checkSSatisfied(Balance,phaseABCY,VoltpABC, ...
phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP, ...
phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ );
PGQG=CalPGQG(Balance,phaseABCY,VoltpABC,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP,phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ );
fprintf('<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϊ%f\n\n',full(max(abs(ub))))
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>end
busNum=length(nodeNum);
Busnum=busNum;
% PQi=setxor(nodeNum,Balance);
QGi=[Balance];
% fprintf('<EFBFBD><EFBFBD>ʼţ<EFBFBD>ٷ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>\n');
[r,c,GB]=find(phaseABCY);
Y=abs(phaseABCY);
Yangle=angle(GB);
Vp3=sparse(ones(busNum*3,1));%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѹ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ
Vp3(2:3:end)=Vp3(2:3:end)*exp(1j*-120/180*pi);
Vp3(3:3:end)=Vp3(3:3:end)*exp(1j*+120/180*pi);
PQi3P=zeros(length(Loadi)*3,1);
PQi3P(1:3:end)=(Loadi-1)*3+1;
PQi3P(2:3:end)=(Loadi-1)*3+2;
PQi3P(3:3:end)=(Loadi-1)*3+3;
Loadi=PQi3P;
%<EFBFBD>ѷ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ҳ<EFBFBD><EFBFBD>Ϊ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
Loadi=[Loadi;(Balance-1)*3+1;(Balance-1)*3+2;(Balance-1)*3+3];
PD3P=sparse(Busnum*3,1);
QD3P=sparse(Busnum*3,1);
phaseASpotLoadP(Balance)=-real(PGQG((Balance-1)*3+1));
phaseBSpotLoadP(Balance)=-real(PGQG((Balance-1)*3+2));
phaseCSpotLoadP(Balance)=-real(PGQG((Balance-1)*3+3));
phaseASpotLoadQ(Balance)=-imag(PGQG((Balance-1)*3+1));
phaseBSpotLoadQ(Balance)=-imag(PGQG((Balance-1)*3+2));
phaseCSpotLoadQ(Balance)=-imag(PGQG((Balance-1)*3+3));
% PD3P(1:3:end)=phaseASpotLoadP*0.9;
% PD3P(2:3:end)=phaseBSpotLoadP*0.9;
% PD3P(3:3:end)=phaseCSpotLoadP*0.9;
% QD3P(1:3:end)=phaseASpotLoadQ*0.9;
% QD3P(2:3:end)=phaseBSpotLoadQ*0.9;
% QD3P(3:3:end)=phaseCSpotLoadQ*0.9;
PD3P=PD3P(Loadi);
QD3P=QD3P(Loadi);
QGi3P=zeros(length(QGi)*3,1);
QGi3P(1:3:end)=(QGi-1)*3+1;
QGi3P(2:3:end)=(QGi-1)*3+2;
QGi3P(3:3:end)=(QGi-1)*3+3;
Vp3m=abs(Vp3);
Vp3a=angle(Vp3);
Balance3P=zeros(length(Balance)*3,1);
Balance3P(1:3:end)=(Balance-1)*3+1;
Balance3P(2:3:end)=(Balance-1)*3+2;
Balance3P(3:3:end)=(Balance-1)*3+3;
% PGA=sum(PD3P(1:3:end))*1.03;
% PGB=sum(PD3P(2:3:end))*1.03;
% PGC=sum(PD3P(3:3:end))*1.03;
% QGA=sum(QD3P(1:3:end))*1.03;
% QGB=sum(QD3P(2:3:end))*1.03;
% QGC=sum(QD3P(3:3:end))*1.03;
% PG3P=sparse( (Balance-1)*3+1:(Balance-1)*3+3,1,[PGA,PGB,PGC],Busnum*3,1);
% QG3P=sparse( (Balance-1)*3+1:(Balance-1)*3+3,1,[QGA,QGB,QGC],Busnum*3,1);
PG3P=real(PGQG);
QG3P=imag(PGQG);
Vp3a((Balance-1)*3+1)=0;
Vp3a((Balance-1)*3+2)=-120/180*pi;
Vp3a((Balance-1)*3+3)=+120/180*pi;
%׼<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͷ<EFBFBD><EFBFBD><EFBFBD>
%<EFBFBD><EFBFBD>ʵֵ
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѹ<EFBFBD><EFBFBD>ֵ
rVoltABCV=zeros(busNum*3,1);
rVoltABCV(1:3:end)=abs(VoltpABC(1,:));
rVoltABCV(2:3:end)=abs(VoltpABC(2,:));
rVoltABCV(3:3:end)=abs(VoltpABC(3,:));
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѹ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
rVoltABCA=zeros(busNum*3,1);
rVoltABCA(1:3:end)=angle(VoltpABC(1,:));
rVoltABCA(2:3:end)=angle(VoltpABC(2,:));
rVoltABCA(3:3:end)=angle(VoltpABC(3,:));
%<EFBFBD><EFBFBD><EFBFBD>ฺ<EFBFBD><EFBFBD>
rPD3P=zeros(busNum*3,1);
rPD3P(1:3:end)=phaseASpotLoadP;
rPD3P(2:3:end)=phaseBSpotLoadP;
rPD3P(3:3:end)=phaseCSpotLoadP;
rQD3P=zeros(busNum*3,1);
rQD3P(1:3:end)=phaseASpotLoadQ;
rQD3P(2:3:end)=phaseBSpotLoadQ;
rQD3P(3:3:end)=phaseCSpotLoadQ;
rPD3P=rPD3P(Loadi);
rQD3P=rQD3P(Loadi);
%% û<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% noLoadi=[1,5,6,10,11];
noLoadi=[1,8,10,12];
% noLoadi=Loadi;
% noLoadi=[1,11];
noLoadi=[5];
% noLoadi=[5,11,19,27,37,43,54,61];
noPQi3P=zeros(length(noLoadi)*3,1);
noPQi3P(1:3:end)=(noLoadi-1)*3+1;
noPQi3P(2:3:end)=(noLoadi-1)*3+2;
noPQi3P(3:3:end)=(noLoadi-1)*3+3;
% noPQi3P=Loadi;
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
sigma=0.03;
VoltSigma=(1+normrnd(0,sigma/3,length(rVoltABCV),1));
mVoltABCV=rVoltABCV.*VoltSigma;
% mVoltABCV(noPQi3P)=rVoltABCV(noPQi3P).*(1+normrnd(0,0.10,length(noPQi3P),1));
PD3PSigma=(1+normrnd(0,sigma,length(rPD3P),1));
mPD3P=rPD3P.*PD3PSigma;
QD3PSigma=(1+normrnd(0,sigma,length(rQD3P),1));
mQD3P=rQD3P.*QD3PSigma;
mPD3P(ismember(Loadi,noPQi3P))=rPD3P(ismember(Loadi,noPQi3P)).*(1+unifrnd(-0.15,0.15,length(noPQi3P),1));
mQD3P(ismember(Loadi,noPQi3P))=rQD3P(ismember(Loadi,noPQi3P)).*(1+unifrnd(-0.15,0.15,length(noPQi3P),1));
%<EFBFBD><EFBFBD><EFBFBD>ⷽ<EFBFBD><EFBFBD>
wVolt=1./(abs(mVoltABCV*sigma/3).^2);
wPD=1./(abs(mPD3P*sigma).^2);
wQD=1./(abs(mQD3P*sigma).^2);
wVolt(setdiff(1:length(wVolt),Loadi))=0;%ֻ<EFBFBD>и<EFBFBD><EFBFBD>ɴ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>е<EFBFBD>ѹ<EFBFBD><EFBFBD><EFBFBD>⡣
完成了Huber方法。 Signed-off-by: dmy@lab <dmy@lab.lab>
2015-06-03 10:35:15 +08:00
% wVolt( noPQi3P)=0;
第一次提交 Signed-off-by: dmy@lab <dmy@lab.lab>
2015-06-02 17:24:13 +08:00
wPD(ismember( Loadi,noPQi3P))=1./(abs(mPD3P(ismember( Loadi,noPQi3P))*0.15).^2);
wQD(ismember( Loadi,noPQi3P))=1./(abs(mQD3P(ismember( Loadi,noPQi3P))*0.15).^2);
%%
% RestraintCount=size(Loadi,1)*2+length(rVoltABCV); %Լ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% RestraintCount=size(Loadi,1)*2; %Լ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
Init_Y=sparse(length(setdiff(1:Busnum*3,Loadi))*2,1);%<EFBFBD><EFBFBD>ѧ<EFBFBD><EFBFBD>һ<EFBFBD><EFBFBD>
KK=0;
ContrlCount=size(Loadi,1)*2+Busnum*6;
kmax=200;
%%
Precision=1e-5;
CenterA=0.1;
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
maxD=100;
tic
Volt=1*Vp3m;
UAngel=Vp3a;
while(maxD>Precision)
if KK>kmax
break;
end
%% <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>OPF
%% <EFBFBD>γɵ<EFBFBD>ʽԼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ſ˱<EFBFBD>
dH_dx = jacobian_M(Busnum,Volt,Y,Yangle,UAngel,r,c); %<EFBFBD>γ<EFBFBD><EFBFBD>ſ˱Ⱦ<EFBFBD><EFBFBD><EFBFBD>
%% <EFBFBD>γɲ<EFBFBD><EFBFBD><EFBFBD>ʽԼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ſ˱<EFBFBD>
%% <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ddg
%% <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>deltF
deltF=func_deltF(wVolt,wPD,wQD,mPD3P,PD3P,QD3P,mQD3P,Volt,mVoltABCV,Busnum,Loadi);
%% <EFBFBD>γɷ<EFBFBD><EFBFBD>̾<EFBFBD><EFBFBD><EFBFBD>
Mat_H=FormH(Busnum,Volt,PG3P,PD3P,QG3P,QD3P,Y,UAngel,r,c,Yangle,Loadi);
%% <EFBFBD><EFBFBD>ʼ<EFBFBD>ⷽ<EFBFBD><EFBFBD>
% fprintf('<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> %d Gap %f\n',KK+1,plotGap(KK+1));
XX=SolveIt(ContrlCount,Balance,Busnum,Loadi,deltF,mPD3P,PD3P,mQD3P,QD3P,mVoltABCV,Volt,Mat_H,wVolt,wPD,wQD,dH_dx);
%%ȡ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
[deltX,deltY]=AssignXX(XX,ContrlCount,Loadi,Balance,Busnum);
[Init_Y,PG,QG,Volt,UAngel,PD3P,QD3P]=Modification(Init_Y,deltX,deltY,PG,QG,Volt,UAngel,ContrlCount,Balance,Busnum,PD3P,QD3P,Loadi);
maxD=max(abs([deltX]));
fprintf('%f\n',full(maxD));
KK=KK+1;
end
totalTime=toc
(rVoltABCV-Volt)./rVoltABCV*100;
PD3P=-Mat_H(Loadi);
QD3P=-Mat_H(Loadi+Busnum*3);
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͳ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%Ŀ<EFBFBD>꺯<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ
JMeasurement=sum(((mVoltABCV-Volt)./mVoltABCV./sigma).^2)+sum(((mPD3P-PD3P)./mPD3P./sigma).^2)+sum(((mQD3P-QD3P)./mQD3P./sigma).^2);
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͳ<EFBFBD><EFBFBD>
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
Busnum=busNum;
mCount=Busnum*3+length(Loadi)*3*2;
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
AME_Volt=sum(sum(abs( abs(rVoltABCV)-abs(Volt))))/3/busNum;
AME_VAngle=sum(sum(abs( rVoltABCA-UAngel)))/3/busNum;
AME_PD=sum(sum(abs(rPD3P-PD3P)))/length(Loadi);
AME_QD=sum(sum(abs(rQD3P-QD3P)))/length(Loadi);
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ<EFBFBD><EFBFBD>
AME_mVolt=sum(sum(abs( mVoltABCV-rVoltABCV)))/3/busNum;
AME_mPD=sum(sum(abs(rPD3P-mPD3P)))/length(Loadi);
AME_mQD=sum(sum(abs(rQD3P-mQD3P)))/length(Loadi);
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϣ
isConverged=1;
if KK>=kmax
isConverged=0;
end
if abs(maxD)>Precision
isConverged=0;
end
fprintf('<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>%d\n',KK);
toc
end