把一次迭代过程的代码复制过来了，还没检查是不是有复制错的地方。
Signed-off-by: dugg@lab-desk <dugg@lab-desk>
This commit is contained in:
dugg@lab-desk
parent
276db0e00e
commit
99aab5404c
|
|
@ -1,4 +1,5 @@
|
|||
*.asv
|
||||
/feeder13
|
||||
*.mexw32
|
||||
/DistributionNetwork-OnlyFortescue
|
||||
#/DistributionNetwork-OnlyFortescue
|
||||
*.7z
|
||||
|
|
@ -0,0 +1,23 @@
|
|||
function [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX1(XX,ContrlCount,RestraintCount,Busnum)
|
||||
% deltX=XX(1:14);
|
||||
% deltY=XX(15:24);
|
||||
% deltZ=XX(25:38);
|
||||
% deltW=XX(39:52);
|
||||
% deltL=XX(53:66);
|
||||
% deltU=XX(67:80);
|
||||
deltX=XX(1:ContrlCount);
|
||||
k1=ContrlCount+2*Busnum;
|
||||
deltY=XX(ContrlCount+1:k1);
|
||||
k2=k1;
|
||||
k1=k2+RestraintCount;
|
||||
deltZ=XX(k2+1:k1);
|
||||
k2=k1;
|
||||
k1=k2+RestraintCount;
|
||||
deltW=XX(k2+1:k1);
|
||||
k2=k1;
|
||||
k1=k2+RestraintCount;
|
||||
deltL=XX(k2+1:k1);
|
||||
k2=k1;
|
||||
k1=k2+RestraintCount;
|
||||
deltU=XX(k2+1:k1);
|
||||
end
|
||||
|
|
@ -0,0 +1 @@
|
|||
Subproject commit ecaed830b547d025369bd6137707f6a96856556a
|
||||
|
|
@ -0,0 +1,9 @@
|
|||
function AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW)
|
||||
tdeltZinx=find(deltZ<0);
|
||||
tdeltWinx=find(deltW>0);
|
||||
t1=-Init_Z(tdeltZinx)./deltZ(tdeltZinx);
|
||||
t2=-Init_W(tdeltWinx)./deltW(tdeltWinx);
|
||||
t3=[t1,t2];
|
||||
t4=min(t3);
|
||||
AlphaD=0.9995*min([t4 1]);
|
||||
end
|
||||
|
|
@ -0,0 +1,9 @@
|
|||
function AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU)
|
||||
tdeltLinx=find(deltL<0);
|
||||
tdeltUinx=find(deltU<0);
|
||||
t1=-Init_L(tdeltLinx)./deltL(tdeltLinx);
|
||||
t2=-Init_U(tdeltUinx)./deltU(tdeltUinx);
|
||||
t3=[t1,t2];
|
||||
t4=min(t3);
|
||||
AlphaP=0.9995*min([t4 1]);
|
||||
end
|
||||
|
|
@ -0,0 +1,7 @@
|
|||
function Mat_G=FormG(Volt,PD,QD)
|
||||
Mat_G=[
|
||||
sparse(PD);
|
||||
sparse(QD);
|
||||
Volt;
|
||||
];
|
||||
end
|
||||
|
|
@ -0,0 +1,14 @@
|
|||
function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle)
|
||||
%%
|
||||
%QDcos=textread('D:\Project\×îС»¯³±Á÷\×îС³±Á÷ËãÀý\Ïɺ£919PDQDglys.txt');
|
||||
%QD(QD~=0)=PD(QD~=0)./tan(QDcos);
|
||||
%QD(QD_NON_ZERO_IND)=QD_NON_ZERO;
|
||||
%%
|
||||
%PD(Loadi)=QD(Loadi)./tan(acos(0.98));
|
||||
AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum);
|
||||
dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt;
|
||||
dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt;
|
||||
|
||||
Mat_H=[dP;dQ;];
|
||||
|
||||
end
|
||||
|
|
@ -0,0 +1,10 @@
|
|||
function Lw=FormLw(Mat_G,Init_U,Busnum,Loadi)
|
||||
VoltU=sparse(1.1*ones(Busnum,1));
|
||||
PDU=sparse(1.5*ones(length(Loadi),1));
|
||||
QDU=sparse(1.5*ones(length(Loadi),1));
|
||||
t1=[PDU;
|
||||
QDU;
|
||||
VoltU];
|
||||
t2=Mat_G+Init_U-t1;
|
||||
Lw=t2;
|
||||
end
|
||||
|
|
@ -0,0 +1,4 @@
|
|||
function Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W)
|
||||
t1=deltF-deltH*Init_Y-deltG*(Init_Z+Init_W);
|
||||
Lx=t1;
|
||||
end
|
||||
|
|
@ -0,0 +1,11 @@
|
|||
function LxComa=FormLxComa(deltF,deltG,deltH,Init_L,Luu,Lul,Init_Z,Init_Y,Lz,Init_U,Init_W,Lw,Lx)
|
||||
%t1=deltG*(Init_Z'+Init_W');%%
|
||||
t2=Lul+diag(Init_Z)*Lz;
|
||||
t3=inv(diag(Init_L));
|
||||
t4=t3*t2;%
|
||||
t5=Luu-diag(Init_W)*Lw;
|
||||
t6=inv(diag(Init_U));
|
||||
t7=t6*t5;%
|
||||
t8=deltG*(t4+t7);%%
|
||||
LxComa=Lx+t8;
|
||||
end
|
||||
|
|
@ -0,0 +1,12 @@
|
|||
function Lz=FormLz(Mat_G,Init_L,Busnum,Loadi)
|
||||
VoltL=sparse(0.9*ones(Busnum,1));
|
||||
PDL=sparse(0*ones(length(Loadi),1));
|
||||
QDL=sparse(0*ones(length(Loadi),1));
|
||||
t1=[PDL;
|
||||
QDL;
|
||||
VoltL
|
||||
];
|
||||
t2=Mat_G-Init_L-t1;
|
||||
Lz=t2;
|
||||
|
||||
end
|
||||
|
|
@ -0,0 +1,10 @@
|
|||
function YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx)
|
||||
YY=[
|
||||
Lx;
|
||||
-Ly;
|
||||
-Lz;
|
||||
-Lw;
|
||||
-Lul;
|
||||
-Luu;
|
||||
];
|
||||
end
|
||||
|
|
@ -0,0 +1,18 @@
|
|||
function [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel,ContrlCount,Balance,Busnum,PD,QD,Loadi)
|
||||
AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU);
|
||||
AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW);
|
||||
Init_Z=Init_Z+AlphaD*deltZ;
|
||||
Init_L=Init_L+AlphaP*deltL;
|
||||
Init_W=Init_W+AlphaD*deltW;
|
||||
Init_U=Init_U+AlphaP*deltU;
|
||||
Init_Y=Init_Y+AlphaD*deltY;
|
||||
t=deltX(1:size(Loadi,1)*2);
|
||||
PD(Loadi)=PD(Loadi)+AlphaP*t(1:length(Loadi));
|
||||
QD(Loadi)=QD(Loadi)+AlphaP*t(length(Loadi)+1:length(Loadi)*2);
|
||||
t=deltX(size(Loadi,1)*2+1:ContrlCount);
|
||||
t(Busnum+Balance)=0;
|
||||
balVolt=Volt(Balance);
|
||||
Volt=Volt+AlphaP*t(1:Busnum);
|
||||
Volt(Balance)=balVolt;
|
||||
UAngel=UAngel+AlphaP*t(Busnum+1:2*Busnum);
|
||||
end
|
||||
|
|
@ -0,0 +1,46 @@
|
|||
function XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,Lx,Balance,Busnum,Loadi)
|
||||
LxComa=FormLxComa(deltF,deltG,deltH,Init_L,Luu,Lul,Init_Z,Init_Y,Lz,Init_U,Init_W,Lw,Lx);
|
||||
H=-deltdeltF+ddh;
|
||||
t1=diag(Init_Z./Init_L-Init_W./Init_U);
|
||||
t2=-deltG*( t1 )*deltG';
|
||||
aa=[
|
||||
(H+t2),deltH;
|
||||
deltH',zeros(length(Init_Y));
|
||||
];
|
||||
yy=[LxComa;-Ly];
|
||||
%% ƽºâ½Úµãµçѹ²»±ä
|
||||
t=size(Loadi,1)*2;
|
||||
aa(t+Balance,:)=0;
|
||||
aa(:,t+Balance)=0;
|
||||
%aa(t+Balance,t+Balance)=1;
|
||||
aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance),1),ContrlCount+2*Busnum,ContrlCount+2*Busnum);
|
||||
deltG(t+Balance,:)=0;
|
||||
%%
|
||||
t=size(Loadi,1)*2+Busnum*1;
|
||||
aa(t+Balance,:)=0;
|
||||
aa(:,t+Balance)=0;
|
||||
%aa(t+Balance,t+Balance)=1;
|
||||
aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance),1),ContrlCount+2*Busnum,ContrlCount+2*Busnum);
|
||||
deltG(t+Balance,:)=0;
|
||||
%%
|
||||
dxdy=aa\yy;
|
||||
%% KLU
|
||||
%spy(aa)
|
||||
%dxdy = klu(aa,'\',full(yy));
|
||||
%%
|
||||
dX=dxdy(1:ContrlCount);
|
||||
dY=dxdy(ContrlCount+1:ContrlCount+2*Busnum);
|
||||
dL=Lz+deltG'*dX;
|
||||
dU=-Lw-deltG'*dX;
|
||||
dZ=-diag(Init_L)\Lul-diag(Init_L)\diag(Init_Z)*dL;
|
||||
dW=-diag(Init_U)\Luu-diag(Init_U)\diag(Init_W)*dU;
|
||||
XX=[
|
||||
dX;
|
||||
dY;
|
||||
dZ;
|
||||
dW;
|
||||
dL;
|
||||
dU;
|
||||
|
||||
];
|
||||
end
|
||||
|
|
@ -0,0 +1,4 @@
|
|||
function ddg=func_ddg()
|
||||
|
||||
ddg=0;
|
||||
end
|
||||
|
|
@ -0,0 +1,48 @@
|
|||
function ddh=func_ddh(Volt,Init_Y,Busnum,Y,UAngel,r,c,Angle,Loadi,ContrlCount)
|
||||
mat_AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum);
|
||||
mat_INV_AngleIJ=mat_AngleIJ';
|
||||
yP=Init_Y(1:length(Init_Y)/2);%暂时改这里 20111227
|
||||
yQ=Init_Y(length(Init_Y)/2+1:length(Init_Y));%暂时改这里 20111227
|
||||
t1=-sparse(1:Busnum,1:Busnum,Y.*cos(mat_INV_AngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP,Busnum,Busnum);
|
||||
t2=sparse(1:Busnum,1:Busnum,sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP,Busnum,Busnum)*Y.*cos(mat_AngleIJ);
|
||||
t3=(t1+t2)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum);
|
||||
t4=-(sparse(1:Busnum,1:Busnum,Y.*cos(mat_AngleIJ)*Volt,Busnum,Busnum) -sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*Y.*cos(mat_INV_AngleIJ) )*sparse(1:Busnum,1:Busnum,sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP,Busnum,Busnum);
|
||||
ddPdTdT=t3+t4;%ok1
|
||||
t1=(-sparse(1:Busnum,1:Busnum,Y.*sin(mat_AngleIJ)*Volt,Busnum,Busnum)+sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*Y.*sin(mat_INV_AngleIJ) )*sparse(1:Busnum,1:Busnum,yP,Busnum,Busnum);
|
||||
t2= -sparse(1:Busnum,1:Busnum, sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP ,Busnum,Busnum)*Y.*sin(mat_AngleIJ)+sparse(1:Busnum,1:Busnum,Y.*sin(mat_INV_AngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP,Busnum,Busnum);
|
||||
ddPdVdT=t1+t2;%ok1
|
||||
t1=diag( Y.*sin(mat_INV_AngleIJ)*diag(Volt)*yP);
|
||||
t2=diag(yP)*Y.*sin(mat_AngleIJ)*diag(Volt);
|
||||
t3=-diag(yP)*diag(Y.*sin(mat_AngleIJ)*Volt);
|
||||
t4=-Y.*sin(mat_INV_AngleIJ)*diag( diag(Volt)*yP );
|
||||
ddPdTdV=t1+t2+t3+t4;%存疑与我的不一样
|
||||
t1=Y.*cos(mat_INV_AngleIJ)*diag(yP);
|
||||
t2=diag(yP)*Y.*cos(mat_AngleIJ);
|
||||
ddPdVdV=t1+t2;
|
||||
t1=-diag(Y.*sin(mat_AngleIJ)*Volt);
|
||||
t2=diag(Volt)*Y.*sin(mat_INV_AngleIJ);
|
||||
t3=(t1+t2)*diag( diag(Volt)*yQ );
|
||||
t4=-diag( diag(Volt)*yQ )*Y.*sin(mat_AngleIJ);
|
||||
t5=diag(Y.*sin(mat_INV_AngleIJ)*diag(Volt)*yQ);
|
||||
t6=-(t4+t5)*diag(Volt);
|
||||
ddQdTdT=t3+t6;%ok1
|
||||
t1=(diag(Y.*cos(mat_AngleIJ)*Volt)-diag(Volt)*Y.*cos(mat_INV_AngleIJ) )*diag(yQ);
|
||||
t2=+diag( diag(Volt)*yQ )*Y.*cos(mat_AngleIJ)-diag(Y.*cos(mat_INV_AngleIJ)*diag(Volt)*yQ);
|
||||
ddQdVdT=t1+t2;
|
||||
t1=Y.*cos(mat_INV_AngleIJ)*diag(diag(Volt)*yQ);
|
||||
t2=diag(yQ)*diag(Y.*cos(mat_AngleIJ)*Volt);
|
||||
t3=-diag(Y.*cos(mat_INV_AngleIJ)*diag(Volt)*yQ);
|
||||
t4=-diag(yQ)*Y.*cos(mat_AngleIJ)*diag(Volt);
|
||||
ddQdTdV=t1+t2+t3+t4;
|
||||
t1=Y.*sin(mat_INV_AngleIJ)*diag(yQ);
|
||||
t2=diag(yQ)*Y.*sin(mat_AngleIJ);
|
||||
ddQdVdV=t1+t2;
|
||||
t=[ddPdVdV+ddQdVdV,ddPdVdT+ddQdVdT;
|
||||
ddPdTdV+ddQdTdV,ddPdTdT+ddQdTdT;
|
||||
];
|
||||
sizeLoadi=size(Loadi,1)*2;
|
||||
ddh=[
|
||||
sparse(sizeLoadi,ContrlCount);
|
||||
sparse(2*Busnum,sizeLoadi),-t;
|
||||
];
|
||||
end
|
||||
|
|
@ -0,0 +1,8 @@
|
|||
function deltF=func_deltF(SEPD,ContrlCount)
|
||||
t=sparse(ones(length(SEPD),1));
|
||||
deltF=[
|
||||
t;
|
||||
sparse(ContrlCount-length(t),1);
|
||||
];
|
||||
|
||||
end
|
||||
|
|
@ -0,0 +1,22 @@
|
|||
function deltG=func_deltG(busNum,Loadi)
|
||||
dgV_dPD=sparse(length(Loadi),busNum);
|
||||
dgV_dQD=sparse(length(Loadi),busNum);
|
||||
dgV_dx=[sparse(1:busNum,1:busNum,ones(busNum,1),busNum,busNum);
|
||||
sparse(busNum,busNum);
|
||||
];
|
||||
%
|
||||
dgPD_dPD=sparse(eye(length(Loadi)));
|
||||
dgPD_dQD=sparse(length(Loadi),length(Loadi));
|
||||
dgPD_dx=sparse(busNum*2,length(Loadi));
|
||||
%
|
||||
dgQD_dPD=sparse(length(Loadi),length(Loadi));
|
||||
dgQD_dQD=sparse(eye(length(Loadi)));
|
||||
dgQD_dx=sparse(busNum*2,length(Loadi));
|
||||
%%
|
||||
|
||||
deltG=[
|
||||
dgPD_dPD,dgQD_dPD,dgV_dPD;
|
||||
dgPD_dQD,dgQD_dQD,dgV_dQD;
|
||||
dgPD_dx,dgQD_dx,dgV_dx;
|
||||
];
|
||||
end
|
||||
|
|
@ -0,0 +1,7 @@
|
|||
function deltH=func_deltH(Busnum,Volt,Y,UAngel,r,c,Angle,Loadi)
|
||||
dH_dPD=[sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum) sparse(size(Loadi,1),Busnum)];
|
||||
dH_dQD=[sparse(size(Loadi,1),Busnum) sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum)];
|
||||
dH_dx = jacobian(Busnum,Volt,Y,Angle,UAngel,r,c); %形成雅克比矩阵
|
||||
deltH=[dH_dPD;dH_dQD;dH_dx'];
|
||||
|
||||
end
|
||||
|
|
@ -0,0 +1,5 @@
|
|||
function deltdeltF=func_deltdeltF(SEPD,ContrlCount)
|
||||
deltdeltF=[
|
||||
0*zeros(ContrlCount);
|
||||
];
|
||||
end
|
||||
|
|
@ -0,0 +1,20 @@
|
|||
function [Jacob]=jacobian(Busnum,Volt,Y,Angle,UAngel,r,c)
|
||||
%**************************************************************************
|
||||
% 程序功能 : 子函数——形成雅可比矩阵Jacobian
|
||||
% 编 者：
|
||||
% 编制时间：2010.12
|
||||
%**************************************************************************
|
||||
%%参照图书馆6楼的书编写
|
||||
%% 分别求雅克比矩阵的子阵H,L,N,J及有功无功分量P,Q
|
||||
AngleIJ=UAngel(r)-UAngel(c)-Angle;
|
||||
mat_AngleIJ=sparse(r,c,AngleIJ,Busnum,Busnum);
|
||||
mat_IvAngleIJ=mat_AngleIJ';
|
||||
H=sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*Y.*sin(mat_IvAngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)-sparse(1:Busnum,1:Busnum,Y.*sin(mat_AngleIJ)*Volt,Busnum,Busnum)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum);
|
||||
N=-sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*Y.*cos(mat_IvAngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)+sparse(1:Busnum,1:Busnum,Y.*cos(mat_AngleIJ)*Volt,Busnum,Busnum)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum);
|
||||
J=sparse(1:Busnum,1:Busnum,Y.*cos(mat_AngleIJ)*Volt,Busnum,Busnum)+Y.*cos(mat_IvAngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum);
|
||||
L=sparse(1:Busnum,1:Busnum,Y.*sin(mat_AngleIJ)*Volt,Busnum,Busnum)+Y.*sin(mat_IvAngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum);
|
||||
t1=[J,L;
|
||||
H,N;
|
||||
]';
|
||||
Jacob=-t1;
|
||||
end
|
||||
168
run.m
168
run.m
|
|
@ -1,9 +1,9 @@
|
|||
clc
|
||||
clear
|
||||
lineZ=readLineZ('D:\share\feeder123\lineParameter.txt');
|
||||
lineZ=readLineZ('feeder13\lineParameter.txt');
|
||||
[ fsY0, fsY1, fsY2,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP ...
|
||||
phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY ...
|
||||
cap]=dataRead(lineZ,'D:\share\feeder123\data.txt');
|
||||
cap]=dataRead(lineZ,'feeder13\data1.txt');
|
||||
a=exp(1j*2*pi/3);
|
||||
Tp2f=1/3*[1 1 1;
|
||||
1 a a^2;
|
||||
|
|
@ -26,6 +26,7 @@ QG=sparse(busNum,1);
|
|||
QGi=[Balance];
|
||||
PD=Pabc/3;
|
||||
QD=Qabc/3;
|
||||
Loadi=find(PD~=0);
|
||||
maxD=100000;% 最大不平衡量
|
||||
EPS=1e-5;
|
||||
k=0;
|
||||
|
|
@ -35,7 +36,6 @@ fsY00=fsY0;
|
|||
fsY22=fsY2;
|
||||
Vf2=sparse(busNum,1);
|
||||
If2=sparse(busNum,1);
|
||||
% Vf0=sparse(busNum,1);
|
||||
Vf0=sparse(busNum,1);
|
||||
If0=sparse(busNum,1);
|
||||
%准备序矩阵
|
||||
|
|
@ -115,58 +115,110 @@ ub=checkSSatisfied(Balance,phaseABCY,VoltpABC, ...
|
|||
phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ );
|
||||
fprintf('最大不平衡量为%f\n\n',full(max(abs(ub))))
|
||||
fprintf('开始牛顿法迭代\n');
|
||||
%%用牛顿法求解
|
||||
[r,c,GB]=find(phaseABCY);
|
||||
Y=abs(phaseABCY);
|
||||
Yangle=angle(GB);
|
||||
Vp3=sparse(ones(busNum*3,1));%给电压赋初值
|
||||
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(PQi)*3,1);
|
||||
PQi3P(1:3:end)=(PQi-1)*3+1;
|
||||
PQi3P(2:3:end)=(PQi-1)*3+2;
|
||||
PQi3P(3:3:end)=(PQi-1)*3+3;
|
||||
PG=0;
|
||||
QG=0;
|
||||
PD3P=sparse(busNum*3,1);
|
||||
QD3P=sparse(busNum*3,1);
|
||||
PD3P(1:3:end)=phaseASpotLoadP;
|
||||
PD3P(2:3:end)=phaseBSpotLoadP;
|
||||
PD3P(3:3:end)=phaseCSpotLoadP;
|
||||
QD3P(1:3:end)=phaseASpotLoadQ;
|
||||
QD3P(2:3:end)=phaseBSpotLoadQ;
|
||||
QD3P(3:3:end)=phaseCSpotLoadQ;
|
||||
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;
|
||||
Vp3a((Balance-1)*3+1)=0;
|
||||
Vp3a((Balance-1)*3+2)=-120/180*pi;
|
||||
Vp3a((Balance-1)*3+3)=+120/180*pi;
|
||||
k=0;
|
||||
maxD=10000;
|
||||
tic
|
||||
while(k<=kmax && maxD> EPS)
|
||||
k=k+1;
|
||||
[dP, dQ, YdotSinVolt, YdotCosVolt, diag_Volt_YdotSin, diag_Volt_YdotCos]=Unbalance(Balance3P,busNum*3, ...
|
||||
PQi3P,PG,QG,QGi3P,PD3P,QD3P,Vp3m,Vp3a,Y,Yangle,r,c,0,0,0,0);
|
||||
maxD=max(abs([dP;dQ;]));
|
||||
jaco=Jacobi(Balance3P,busNum*3,QGi3P,Vp3m,YdotSinVolt,YdotCosVolt,diag_Volt_YdotSin,diag_Volt_YdotCos);%雅克比矩阵
|
||||
[dV, dVangle]=Solv(busNum*3,jaco,dP,dQ);%解出修正量
|
||||
[Vp3m, Vp3a]=Modify(Vp3m,Vp3a,dV,dVangle,1);
|
||||
fprintf('第 %d 次迭代， 最大不平衡量为 %f\n',k,full(maxD));
|
||||
fprintf('迭代时间%f\n',toc);
|
||||
end
|
||||
NewtonToc=toc;
|
||||
fprintf('牛顿法计算时间 %f\n',NewtonToc);
|
||||
fprintf('加速比为%f\n',NewtonToc/FortiscueToc);
|
||||
VoltpA=Vp3m(1:3:end).*exp(1j*Vp3a(1:3:end));
|
||||
VoltpB=Vp3m(2:3:end).*exp(1j*Vp3a(2:3:end));
|
||||
VoltpC=Vp3m(3:3:end).*exp(1j*Vp3a(3:3:end));
|
||||
|
||||
%% 用牛顿法求解begin
|
||||
% [r,c,GB]=find(phaseABCY);
|
||||
% Y=abs(phaseABCY);
|
||||
% Yangle=angle(GB);
|
||||
% Vp3=sparse(ones(busNum*3,1));%给电压赋初值
|
||||
% 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(PQi)*3,1);
|
||||
% PQi3P(1:3:end)=(PQi-1)*3+1;
|
||||
% PQi3P(2:3:end)=(PQi-1)*3+2;
|
||||
% PQi3P(3:3:end)=(PQi-1)*3+3;
|
||||
% PG=0;
|
||||
% QG=0;
|
||||
% PD3P=sparse(busNum*3,1);
|
||||
% QD3P=sparse(busNum*3,1);
|
||||
% PD3P(1:3:end)=phaseASpotLoadP;
|
||||
% PD3P(2:3:end)=phaseBSpotLoadP;
|
||||
% PD3P(3:3:end)=phaseCSpotLoadP;
|
||||
% QD3P(1:3:end)=phaseASpotLoadQ;
|
||||
% QD3P(2:3:end)=phaseBSpotLoadQ;
|
||||
% QD3P(3:3:end)=phaseCSpotLoadQ;
|
||||
% 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;
|
||||
% Vp3a((Balance-1)*3+1)=0;
|
||||
% Vp3a((Balance-1)*3+2)=-120/180*pi;
|
||||
% Vp3a((Balance-1)*3+3)=+120/180*pi;
|
||||
% k=0;
|
||||
% maxD=10000;
|
||||
% tic
|
||||
% while(k<=kmax && maxD> EPS)
|
||||
% k=k+1;
|
||||
% [dP, dQ, YdotSinVolt, YdotCosVolt, diag_Volt_YdotSin, diag_Volt_YdotCos]=Unbalance(Balance3P,busNum*3, ...
|
||||
% PQi3P,PG,QG,QGi3P,PD3P,QD3P,Vp3m,Vp3a,Y,Yangle,r,c,0,0,0,0);
|
||||
% maxD=max(abs([dP;dQ;]));
|
||||
% jaco=Jacobi(Balance3P,busNum*3,QGi3P,Vp3m,YdotSinVolt,YdotCosVolt,diag_Volt_YdotSin,diag_Volt_YdotCos);%雅克比矩阵
|
||||
% [dV, dVangle]=Solv(busNum*3,jaco,dP,dQ);%解出修正量
|
||||
% [Vp3m, Vp3a]=Modify(Vp3m,Vp3a,dV,dVangle,1);
|
||||
% fprintf('第 %d 次迭代， 最大不平衡量为 %f\n',k,full(maxD));
|
||||
% fprintf('迭代时间%f\n',toc);
|
||||
% end
|
||||
% NewtonToc=toc;
|
||||
% fprintf('牛顿法计算时间 %f\n',NewtonToc);
|
||||
% fprintf('加速比为%f\n',NewtonToc/FortiscueToc);
|
||||
% VoltpA=Vp3m(1:3:end).*exp(1j*Vp3a(1:3:end));
|
||||
% VoltpB=Vp3m(2:3:end).*exp(1j*Vp3a(2:3:end));
|
||||
% VoltpC=Vp3m(3:3:end).*exp(1j*Vp3a(3:3:end));
|
||||
%% 用牛顿法求解end
|
||||
%% 开始进入状态估计
|
||||
%状态量
|
||||
SEVoltpA=sparse(ones(busNum,1));
|
||||
SEVoltpB=sparse(ones(busNum,1)).*exp(1j*-120/180*pi);
|
||||
SEVoltpC=sparse(ones(busNum,1)).*exp(1j*+120/180*pi);
|
||||
SEphaseASpotLoadP=zeros(length(phaseASpotLoadP),1);
|
||||
SEphaseBSpotLoadP=zeros(length(phaseBSpotLoadP),1);
|
||||
SEphaseCSpotLoadP=zeros(length(phaseCSpotLoadP),1);
|
||||
SEphaseASpotLoadQ=zeros(length(phaseASpotLoadQ),1);
|
||||
SEphaseBSpotLoadQ=zeros(length(phaseBSpotLoadQ),1);
|
||||
SEphaseCSpotLoadQ=zeros(length(phaseCSpotLoadQ),1);
|
||||
SEVmf1=sparse(ones(busNum,1));
|
||||
SEVaf1=sparse(zeros(busNum,1));
|
||||
SEPD=sparse(length(Loadi),1);
|
||||
SEQD=sparse(length(Loadi),1);
|
||||
KK=0;
|
||||
plotGap=zeros(1,60);
|
||||
%初始化
|
||||
%状态量为 SEPD SEQD SEVmf1 SEVaf1
|
||||
RestraintCount=length(SEVmf1)+length(Loadi)*2;
|
||||
ContrlCount=length(SEVmf1)*2+length(Loadi)*2;
|
||||
CenterA=0.1;
|
||||
Init_Z=sparse(ones(RestraintCount,1));
|
||||
Init_W=sparse(-1*ones(RestraintCount,1));
|
||||
Init_L=1*sparse(ones(RestraintCount,1));
|
||||
Init_U=1*sparse(ones(RestraintCount,1));
|
||||
Init_Y=sparse(2*busNum,1);%等式约束乘子
|
||||
Gap=(Init_L'*Init_Z-Init_U'*Init_W);
|
||||
Init_u=Gap/2/RestraintCount*CenterA;
|
||||
deltH=func_deltH(busNum,SEVmf1,fsY1amp,SEVaf1,r,c,fsY1ang,Loadi);
|
||||
deltG=func_deltG(busNum,Loadi);
|
||||
L_1Z=diag(Init_Z./Init_L);
|
||||
U_1W=diag(Init_W./Init_U);
|
||||
deltdeltF=func_deltdeltF(SEPD,ContrlCount);
|
||||
ddh=func_ddh(SEVmf1,Init_Y,busNum,fsY1amp,SEVaf1,r,c,fsY1ang,Loadi,ContrlCount);
|
||||
ddg=func_ddg();
|
||||
deltF=func_deltF(SEPD,ContrlCount);
|
||||
Luu=Init_U.*Init_W+Init_u*ones(RestraintCount,1);
|
||||
Lul=Init_L.*Init_Z-Init_u*ones(RestraintCount,1);
|
||||
Mat_G=FormG(SEVaf1,SEPD,SEQD);
|
||||
Mat_H=FormH(busNum,SEVmf1,PG,SEPD,QG,SEQD,fsY1amp,SEVaf1,r,c,fsY1ang);
|
||||
Ly=Mat_H;
|
||||
Lz=FormLz(Mat_G,Init_L,busNum,Loadi);
|
||||
Lw=FormLw(Mat_G,Init_U,busNum,Loadi);
|
||||
Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
|
||||
YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
|
||||
%% 开始解方程
|
||||
fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1));
|
||||
XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,Lx,Balance,busNum,Loadi);
|
||||
[deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,busNum);
|
||||
[Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,SEVmf1,SEVaf1,ContrlCount,Balance,busNum,PD,QD,Loadi);
|
||||
Gap=(Init_L'*Init_Z-Init_U'*Init_W);
|
||||
Loading…
Reference in New Issue