facat
/
threephasese-onlyenquationconstraint-gn
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

这个版本是只有电压幅值和相角为变量，包含0注入约束。 

Signed-off-by: dmy@lab <dmy@lab.lab>
This commit is contained in:
dmy@lab 2015-06-02 17:15:04 +08:00
parent ad81ecb6cc
commit 1714aed115
7 changed files with 171 additions and 99 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
AssignXX.m
View File

@ -1,18 +1,13 @@
function [deltX,deltY]=AssignXX(XX,ContrlCount,Loadi,Balance,Busnum) function [deltX,deltY]=AssignXX(XX,ContrlCount,Loadi,Balance,Busnum)
% deltX=XX(1:14); deltX=XX(1:Busnum*3*2);
% deltY=XX(15:24); t=0;
% deltZ=XX(25:38);
% deltW=XX(39:52);
% deltL=XX(53:66);
% deltU=XX(67:80);
deltX=XX(1:ContrlCount);
t=size(Loadi,1)*2;
deltX(t+(Balance-1)*3+1)=0; deltX(t+(Balance-1)*3+1)=0;
deltX(t+(Balance-1)*3+2)=0; deltX(t+(Balance-1)*3+2)=0;
deltX(t+(Balance-1)*3+3)=0; deltX(t+(Balance-1)*3+3)=0;
t=size(Loadi,1)*2+Busnum*3; t=Busnum*3;
deltX(t+(Balance-1)*3+1)=0; deltX(t+(Balance-1)*3+1)=0;
deltX(t+(Balance-1)*3+2)=0; deltX(t+(Balance-1)*3+2)=0;
deltX(t+(Balance-1)*3+3)=0; deltX(t+(Balance-1)*3+3)=0;
deltY=XX(ContrlCount+1:end); t=Busnum*3*2;
deltY=XX(t+1:end);
end end

6
FormH.m
View File

@ -6,8 +6,10 @@ function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi)
%% %%
%PD(Loadi)=QD(Loadi)./tan(acos(0.98)); %PD(Loadi)=QD(Loadi)./tan(acos(0.98));
AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum*3,Busnum*3); AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum*3,Busnum*3);
dP=PG-sparse(Loadi,1,PD,Busnum*3,1)-diag(Volt)*Y.*cos(AngleIJ)*Volt; % dP=-sparse(Loadi,1,PD,Busnum*3,1)-diag(Volt)*Y.*cos(AngleIJ)*Volt;
dQ=QG-sparse(Loadi,1,QD,Busnum*3,1)-diag(Volt)*Y.*sin(AngleIJ)*Volt; % dQ=-sparse(Loadi,1,QD,Busnum*3,1)-diag(Volt)*Y.*sin(AngleIJ)*Volt;
dP=diag(Volt)*Y.*cos(AngleIJ)*Volt;
dQ=diag(Volt)*Y.*sin(AngleIJ)*Volt;
Mat_H=[dP;dQ;]; Mat_H=[dP;dQ;];
Mat_H=sparse(Mat_H); Mat_H=sparse(Mat_H);

10
Modification.m
View File

@ -1,12 +1,9 @@
function [Init_Y,PG,QG,Volt,UAngel,PD,QD]=Modification(Init_Y,deltX,deltY,PG,QG,Volt,UAngel,ContrlCount,Balance,Busnum,PD,QD,Loadi) function [Init_Y,PG,QG,Volt,UAngel,PD,QD]=Modification(Init_Y,deltX,deltY,PG,QG,Volt,UAngel,ContrlCount,Balance,Busnum,PD,QD,Loadi)
Init_Y=Init_Y+deltY'; Init_Y=Init_Y+deltY;
t=deltX(1:size(Loadi,1)*2); t=deltX(1:size(Loadi,1)*2);
PD=PD+t(1:length(Loadi)); PD=PD+t(1:length(Loadi));
QD=QD+t(length(Loadi)+1:length(Loadi)*2); QD=QD+t(length(Loadi)+1:length(Loadi)*2);
t=deltX(size(Loadi,1)*2+1:ContrlCount)'; t=deltX(1:Busnum*3*2)';
t(Busnum*3+(Balance-1)*3+1)=0;
t(Busnum*3+(Balance-1)*3+2)=0;
t(Busnum*3+(Balance-1)*3+3)=0;
%Volt=Volt+AlphaP*t(2:2:2*Busnum);暂时改一下20111227 %Volt=Volt+AlphaP*t(2:2:2*Busnum);暂时改一下20111227
%UAngel=UAngel+AlphaP*t(1:2:2*Busnum);暂时改一下20111227 %UAngel=UAngel+AlphaP*t(1:2:2*Busnum);暂时改一下20111227
% balVolt=Volt( (Balance-1)*3+1 ); % balVolt=Volt( (Balance-1)*3+1 );
@ -14,5 +11,8 @@ Volt=Volt+t(1:Busnum*3)';
Volt( (Balance-1)*3+1 )=1; Volt( (Balance-1)*3+1 )=1;
Volt( (Balance-1)*3+2 )=1; Volt( (Balance-1)*3+2 )=1;
Volt( (Balance-1)*3+3 )=1; Volt( (Balance-1)*3+3 )=1;
t(Busnum*3+(Balance-1)*3+1)=0;
t(Busnum*3+(Balance-1)*3+2)=0;
t(Busnum*3+(Balance-1)*3+3)=0;
UAngel=UAngel+t(Busnum*3+1:2*Busnum*3)'; UAngel=UAngel+t(Busnum*3+1:2*Busnum*3)';
end end

95
OPF.m
View File

@ -3,10 +3,10 @@ function [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,AME_mVolt,AME_mPD,AME_m
tic tic
clc clc
clear clear
lineZ=readLineZ('./../DistributionNetwork-Power2Current/modified-feeder69\lineParameter.txt'); lineZ=readLineZ('feeder13\lineParameter.txt');
[ fsY0, fsY1, fsY2,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP ... [ fsY0, fsY1, fsY2,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP ...
phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY ... phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY ...
cap]=dataRead(lineZ,'./../DistributionNetwork-Power2Current/modified-feeder69\data-某相没负荷.txt'); cap]=dataRead(lineZ,'feeder13\data1.txt');
% phaseASpotLoadP(phaseASpotLoadP==0)=0.002; % phaseASpotLoadP(phaseASpotLoadP==0)=0.002;
% phaseBSpotLoadP(phaseBSpotLoadP==0)=0.002; % phaseBSpotLoadP(phaseBSpotLoadP==0)=0.002;
% phaseCSpotLoadP(phaseCSpotLoadP==0)=0.002; % phaseCSpotLoadP(phaseCSpotLoadP==0)=0.002;
@ -15,36 +15,36 @@ lineZ=readLineZ('./../DistributionNetwork-Power2Current/modified-feeder69\linePa
% phaseCSpotLoadQ(phaseCSpotLoadQ==0)=0.002; % phaseCSpotLoadQ(phaseCSpotLoadQ==0)=0.002;
%负荷不平衡 %负荷不平衡
for I=1:length(phaseASpotLoadP) % for I=1:length(phaseASpotLoadP)
roll=mod(round(rand()*10),3); % roll=mod(round(rand()*10),3);
if roll==0 % if roll==0
phaseBSpotLoadP(I)=phaseBSpotLoadP(I)+phaseASpotLoadP(I)*.45; % phaseBSpotLoadP(I)=phaseBSpotLoadP(I)+phaseASpotLoadP(I)*.45;
phaseCSpotLoadP(I)=phaseCSpotLoadP(I)+phaseASpotLoadP(I)*.45; % phaseCSpotLoadP(I)=phaseCSpotLoadP(I)+phaseASpotLoadP(I)*.45;
phaseASpotLoadP(I)=phaseASpotLoadP(I)*.10; % phaseASpotLoadP(I)=phaseASpotLoadP(I)*.10;
%
phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)+phaseASpotLoadQ(I)*.45; % phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)+phaseASpotLoadQ(I)*.45;
phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)+phaseASpotLoadQ(I)*.45; % phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)+phaseASpotLoadQ(I)*.45;
phaseASpotLoadQ(I)=phaseASpotLoadQ(I)*.10; % phaseASpotLoadQ(I)=phaseASpotLoadQ(I)*.10;
end % end
if roll==1 % if roll==1
phaseASpotLoadP(I)=phaseASpotLoadP(I)+phaseBSpotLoadP(I)*.45; % phaseASpotLoadP(I)=phaseASpotLoadP(I)+phaseBSpotLoadP(I)*.45;
phaseCSpotLoadP(I)=phaseCSpotLoadP(I)+phaseBSpotLoadP(I)*.45; % phaseCSpotLoadP(I)=phaseCSpotLoadP(I)+phaseBSpotLoadP(I)*.45;
phaseBSpotLoadP(I)=phaseBSpotLoadP(I)*.10; % phaseBSpotLoadP(I)=phaseBSpotLoadP(I)*.10;
%
phaseASpotLoadQ(I)=phaseASpotLoadQ(I)+phaseBSpotLoadQ(I)*.45; % phaseASpotLoadQ(I)=phaseASpotLoadQ(I)+phaseBSpotLoadQ(I)*.45;
phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)+phaseBSpotLoadQ(I)*.45; % phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)+phaseBSpotLoadQ(I)*.45;
phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)*.10; % phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)*.10;
end % end
if roll==2 % if roll==2
phaseASpotLoadP(I)=phaseASpotLoadP(I)+phaseCSpotLoadP(I)*.45; % phaseASpotLoadP(I)=phaseASpotLoadP(I)+phaseCSpotLoadP(I)*.45;
phaseBSpotLoadP(I)=phaseBSpotLoadP(I)+phaseCSpotLoadP(I)*.45; % phaseBSpotLoadP(I)=phaseBSpotLoadP(I)+phaseCSpotLoadP(I)*.45;
phaseCSpotLoadP(I)=phaseCSpotLoadP(I)*.10; % phaseCSpotLoadP(I)=phaseCSpotLoadP(I)*.10;
%
phaseASpotLoadQ(I)=phaseASpotLoadQ(I)+phaseCSpotLoadQ(I)*.45; % phaseASpotLoadQ(I)=phaseASpotLoadQ(I)+phaseCSpotLoadQ(I)*.45;
phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)+phaseCSpotLoadQ(I)*.45; % phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)+phaseCSpotLoadQ(I)*.45;
phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)*.10; % phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)*.10;
end % end
end % end
%% 潮流计算begin %% 潮流计算begin
@ -143,7 +143,6 @@ while(k<=kmax+10 && maxD> EPS)
%Vf2=fsY2\If2; %Vf2=fsY2\If2;
Vf2=fsY2Q*(fsY2U\(fsY2L\(fsY2P*(fsY2R\If2)))); Vf2=fsY2Q*(fsY2U\(fsY2L\(fsY2P*(fsY2R\If2))));
fprintf('迭代时间%f\n',toc); fprintf('迭代时间%f\n',toc);
% %
end end
FortiscueToc=toc; FortiscueToc=toc;
@ -181,14 +180,22 @@ PQi3P(1:3:end)=(Loadi-1)*3+1;
PQi3P(2:3:end)=(Loadi-1)*3+2; PQi3P(2:3:end)=(Loadi-1)*3+2;
PQi3P(3:3:end)=(Loadi-1)*3+3; PQi3P(3:3:end)=(Loadi-1)*3+3;
Loadi=PQi3P; Loadi=PQi3P;
%把发电机也作为负荷
Loadi=[Loadi;(Balance-1)*3+1;(Balance-1)*3+2;(Balance-1)*3+3];
PD3P=sparse(Busnum*3,1); PD3P=sparse(Busnum*3,1);
QD3P=sparse(Busnum*3,1); QD3P=sparse(Busnum*3,1);
PD3P(1:3:end)=phaseASpotLoadP*0.9; phaseASpotLoadP(Balance)=-real(PGQG((Balance-1)*3+1));
PD3P(2:3:end)=phaseBSpotLoadP*0.9; phaseBSpotLoadP(Balance)=-real(PGQG((Balance-1)*3+2));
PD3P(3:3:end)=phaseCSpotLoadP*0.9; phaseCSpotLoadP(Balance)=-real(PGQG((Balance-1)*3+3));
QD3P(1:3:end)=phaseASpotLoadQ*0.9; phaseASpotLoadQ(Balance)=-imag(PGQG((Balance-1)*3+1));
QD3P(2:3:end)=phaseBSpotLoadQ*0.9; phaseBSpotLoadQ(Balance)=-imag(PGQG((Balance-1)*3+2));
QD3P(3:3:end)=phaseCSpotLoadQ*0.9; 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); PD3P=PD3P(Loadi);
QD3P=QD3P(Loadi); QD3P=QD3P(Loadi);
QGi3P=zeros(length(QGi)*3,1); QGi3P=zeros(length(QGi)*3,1);
@ -274,7 +281,7 @@ wQD(ismember( Loadi,noPQi3P))=1./(abs(mQD3P(ismember( Loadi,noPQi3P))*0.15).^2);
%% %%
% RestraintCount=size(Loadi,1)*2+length(rVoltABCV); %约束条件数 % RestraintCount=size(Loadi,1)*2+length(rVoltABCV); %约束条件数
% RestraintCount=size(Loadi,1)*2; %约束条件数 % RestraintCount=size(Loadi,1)*2; %约束条件数
Init_Y=sparse(1,2*Busnum*3);%与学姐一致 Init_Y=sparse(length(setdiff(1:Busnum*3,Loadi))*2,1);%与学姐一致
KK=0; KK=0;
ContrlCount=size(Loadi,1)*2+Busnum*6; ContrlCount=size(Loadi,1)*2+Busnum*6;
kmax=200; kmax=200;
@ -294,7 +301,7 @@ while(maxD>Precision)
end end
%% 开始计算OPF %% 开始计算OPF
%% 形成等式约束的雅克比 %% 形成等式约束的雅克比
deltH=func_deltH(Busnum,Volt,Y,UAngel,r,c,Yangle,Loadi); dH_dx = jacobian_M(Busnum,Volt,Y,Yangle,UAngel,r,c); %形成雅克比矩阵
%% 形成不等式约束的雅克比 %% 形成不等式约束的雅克比
%% 开始构建ddg %% 开始构建ddg
%% 开始构建deltF %% 开始构建deltF
@ -302,8 +309,8 @@ while(maxD>Precision)
%% 形成方程矩阵 %% 形成方程矩阵
Mat_H=FormH(Busnum,Volt,PG3P,PD3P,QG3P,QD3P,Y,UAngel,r,c,Yangle,Loadi); Mat_H=FormH(Busnum,Volt,PG3P,PD3P,QG3P,QD3P,Y,UAngel,r,c,Yangle,Loadi);
%% 开始解方程 %% 开始解方程
% fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1)); % fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1));
XX=SolveIt(deltH,ContrlCount,Balance,Busnum,Loadi,deltF,mPD3P,PD3P,mQD3P,QD3P,mVoltABCV,Volt,Mat_H,wVolt,wPD,wQD); XX=SolveIt(ContrlCount,Balance,Busnum,Loadi,deltF,mPD3P,PD3P,mQD3P,QD3P,mVoltABCV,Volt,Mat_H,wVolt,wPD,wQD,dH_dx);
%%取各分量 %%取各分量
[deltX,deltY]=AssignXX(XX,ContrlCount,Loadi,Balance,Busnum); [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); [Init_Y,PG,QG,Volt,UAngel,PD3P,QD3P]=Modification(Init_Y,deltX,deltY,PG,QG,Volt,UAngel,ContrlCount,Balance,Busnum,PD3P,QD3P,Loadi);
@ -313,6 +320,8 @@ while(maxD>Precision)
end end
totalTime=toc totalTime=toc
(rVoltABCV-Volt)./rVoltABCV*100; (rVoltABCV-Volt)./rVoltABCV*100;
PD3P=-Mat_H(Loadi);
QD3P=-Mat_H(Loadi+Busnum*3);
%% 计算统计量 %% 计算统计量
%目标函数均值 %目标函数均值
JMeasurement=sum(((mVoltABCV-Volt)./mVoltABCV./sigma).^2)+sum(((mPD3P-PD3P)./mPD3P./sigma).^2)+sum(((mQD3P-QD3P)./mQD3P./sigma).^2); JMeasurement=sum(((mVoltABCV-Volt)./mVoltABCV./sigma).^2)+sum(((mPD3P-PD3P)./mPD3P./sigma).^2)+sum(((mQD3P-QD3P)./mQD3P./sigma).^2);

67
SolveIt.m
View File

@ -1,59 +1,66 @@
function XX=SolveIt(deltH,ContrlCount,Balance,Busnum,Loadi,deltF,PD0,PD,QD0,QD,mVolt,Volt,Mat_H,wVolt,wPD,wQD) function XX=SolveIt(ContrlCount,Balance,Busnum,Loadi,deltF,PD0,PD,QD0,QD,mVolt,Volt,Mat_H,wVolt,wPD,wQD,dH_dx)
% LxComa=FormLxComa(Lx);%Lx=deltF-deltH*Init_Y'; % dH_dx
% Lx=LxComa; dP_x=-dH_dx(Loadi,:);
% H=-deltdeltF+ddh; dQ_x=-dH_dx(Loadi+Busnum*3,:);
% aa=[ dV_x=[sparse(1:Busnum*3,1:Busnum*3,1,Busnum*3,Busnum*3),sparse(Busnum*3,Busnum*3)];
% H,deltH; H=[dP_x;dQ_x;dV_x];
% deltH',zeros(size(Init_Y,2)); c=dH_dx(setdiff(1:Busnum*3,Loadi),:);
% ]; c=[c;dH_dx(setdiff(1:Busnum*3,Loadi)+Busnum*3,:)];
aa=[ aa=[H'*diag([wPD;wQD;wVolt])*H c'
deltF'*diag([wPD ;wQD ;wVolt])*deltF deltH' c zeros(size(c,1),size(c',2))];
deltH zeros(size(deltH,1),size(deltH',2)) % aa=H'*eye(length(H))*H;
]; PD=-Mat_H(Loadi,:);
% yy=[Lx;-Ly]; QD=-Mat_H(Loadi+Busnum*3,:);
% t3=2*wPD.*(PD-PD0);
% t4=2*wQD.*(QD-QD0);
% t5=2*wVolt.*(Volt-mVolt);
deltZ=[(PD0-PD); deltZ=[(PD0-PD);
(QD0-QD); (QD0-QD);
(mVolt-Volt); (mVolt-Volt);
]; ];
yy=[ yy=[
deltF'*diag([wPD; wQD; wVolt])*deltZ; H'*diag([wPD;wQD;wVolt])*deltZ;
-Mat_H ; -Mat_H(setdiff(1:Busnum*3,Loadi));
-Mat_H(setdiff(1:Busnum*3,Loadi)+Busnum*3);
]; ];
% yy=[
% H'*eye(length(H))*deltZ;
% ];
%% 平衡节点电压不变 %% 平衡节点电压不变
t=size(Loadi,1)*2; t=0;
aa(t+(Balance-1)*3+1,:)=0; aa(t+(Balance-1)*3+1,:)=0;
aa(t+(Balance-1)*3+2,:)=0; aa(t+(Balance-1)*3+2,:)=0;
aa(t+(Balance-1)*3+3,:)=0; aa(t+(Balance-1)*3+3,:)=0;
aa(:,t+(Balance-1)*3+1)=0; aa(:,t+(Balance-1)*3+1)=0;
aa(:,t+(Balance-1)*3+2)=0; aa(:,t+(Balance-1)*3+2)=0;
aa(:,t+(Balance-1)*3+3)=0; aa(:,t+(Balance-1)*3+3)=0;
%aa(t+Balance,t+Balance)=1; aa=aa+sparse(t+(Balance-1)*3+1,t+(Balance-1)*3+1,ones(length(Balance),1),size(c,1)+2*Busnum*3,size(c,1)+2*Busnum*3);
aa=aa+sparse(t+(Balance-1)*3+1,t+(Balance-1)*3+1,ones(length(Balance),1),ContrlCount+2*Busnum*3,ContrlCount+2*Busnum*3); aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),size(c,1)+2*Busnum*3,size(c,1)+2*Busnum*3);
aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),ContrlCount+2*Busnum*3,ContrlCount+2*Busnum*3); aa=aa+sparse(t+(Balance-1)*3+3,t+(Balance-1)*3+3,ones(length(Balance),1),size(c,1)+2*Busnum*3,size(c,1)+2*Busnum*3);
aa=aa+sparse(t+(Balance-1)*3+3,t+(Balance-1)*3+3,ones(length(Balance),1),ContrlCount+2*Busnum*3,ContrlCount+2*Busnum*3);
% aa=aa+sparse(t+(Balance-1)*3+1,t+(Balance-1)*3+1,ones(length(Balance),1),2*Busnum*3,2*Busnum*3);
% aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),2*Busnum*3,2*Busnum*3);
% aa=aa+sparse(t+(Balance-1)*3+3,t+(Balance-1)*3+3,ones(length(Balance),1),2*Busnum*3,2*Busnum*3);
%% %%
t=size(Loadi,1)*2+Busnum*3; t=Busnum*3;
aa(t+(Balance-1)*3+1,:)=0; aa(t+(Balance-1)*3+1,:)=0;
aa(t+(Balance-1)*3+2,:)=0; aa(t+(Balance-1)*3+2,:)=0;
aa(t+(Balance-1)*3+3,:)=0; aa(t+(Balance-1)*3+3,:)=0;
aa(:,t+(Balance-1)*3+1)=0; aa(:,t+(Balance-1)*3+1)=0;
aa(:,t+(Balance-1)*3+2)=0; aa(:,t+(Balance-1)*3+2)=0;
aa(:,t+(Balance-1)*3+3)=0; aa(:,t+(Balance-1)*3+3)=0;
%aa(t+Balance,t+Balance)=1; aa=aa+sparse(t+(Balance-1)*3+1,t+(Balance-1)*3+1,ones(length(Balance),1),size(c,1)+2*Busnum*3,size(c,1)+2*Busnum*3);
aa=aa+sparse(t+(Balance-1)*3+1,t+(Balance-1)*3+1,ones(length(Balance),1),ContrlCount+2*Busnum*3,ContrlCount+2*Busnum*3); aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),size(c,1)+2*Busnum*3,size(c,1)+2*Busnum*3);
aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),ContrlCount+2*Busnum*3,ContrlCount+2*Busnum*3); aa=aa+sparse(t+(Balance-1)*3+3,t+(Balance-1)*3+3,ones(length(Balance),1),size(c,1)+2*Busnum*3,size(c,1)+2*Busnum*3);
aa=aa+sparse(t+(Balance-1)*3+3,t+(Balance-1)*3+3,ones(length(Balance),1),ContrlCount+2*Busnum*3,ContrlCount+2*Busnum*3); % aa=aa+sparse(t+(Balance-1)*3+1,t+(Balance-1)*3+1,ones(length(Balance),1),2*Busnum*3,2*Busnum*3);
% aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),2*Busnum*3,2*Busnum*3);
% aa=aa+sparse(t+(Balance-1)*3+3,t+(Balance-1)*3+3,ones(length(Balance),1),2*Busnum*3,2*Busnum*3);
%% %%
dxdy=aa\yy; dxdy=aa\yy;
%% KLU %% KLU
%spy(aa) %spy(aa)
%dxdy = klu(aa,'\',full(yy)); %dxdy = klu(aa,'\',full(yy));
%% %%
dX=dxdy(1:ContrlCount); dX=dxdy(1:Busnum*3*2);
dY=dxdy(ContrlCount+1:end); dY=dxdy(Busnum*3*2+1:end);
XX=[ XX=[
dX; dX;
dY; dY;

60
feeder13/data3.txt Normal file
View File

@ -0,0 +1,60 @@
650 4.16
0
600 632 645 50
602 632 633 50
601 633 634 60
600 645 646 30
601 650 632 200
602 684 652 80
601 632 671 200
600 671 684 30
600 671 692 30
602 684 611 30
601 692 675 50
0
692 0 0 0
0
684 120 90 120 90 160 110
633 120 90 120 90 160 110
632 120 90 120 90 160 110
634 120 90 120 90 160 110
645 40 30 50 40 60 50
646 100 50 70 40 60 40
652 20 10 60 30 40 40
671 385 220 385 220 385 220
675 68 60 290 212 285 190
692 260 200 40 40 70 40
611 80 20 60 40 30 20
0
600 632 645 500
602 632 633 500
601 633 634 600
600 645 646 300
601 650 632 2000
602 684 652 800
601 632 671 2000
600 671 684 300
602 671 680 1000
600 671 692 300
602 684 611 300
601 692 675 500
0
632 0 0 0
0
634 160 110 120 90 120 90
645 123 102 170 125 133 110
646 172 150 230 132 0 211 160
652 128 86 30 10 80 70
671 385 220 385 220 385 220
675 485 190 68 60 290 212
692 142 140 153 121 170 151
611 156 142 185 60 170 80
0
634 160 110 120 90 120 90
645 123 102 170 125 133 110
646 17 15 23 13 21 16
652 128 86 30 10 80 70
675 85 10 68 60 20 12
611 52 42 85 60 70 40
680 12 4 10 6 7 4

17
jacobian_M.m
View File

@ -8,13 +8,12 @@ function [Jacob]=jacobian_M3(Busnum,Volt,Y,Angle,UAngel,r,c)
%% 分别求雅克比矩阵的子阵H,L,N,J及有功无功分量P,Q %% 分别求雅克比矩阵的子阵H,L,N,J及有功无功分量P,Q
AngleIJ=UAngel(r)-UAngel(c)-Angle; AngleIJ=UAngel(r)-UAngel(c)-Angle;
mat_AngleIJ=sparse(r,c,AngleIJ,Busnum*3,Busnum*3); mat_AngleIJ=sparse(r,c,AngleIJ,Busnum*3,Busnum*3);
mat_IvAngleIJ=mat_AngleIJ'; H=sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)*Y.*sin(mat_AngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)-sparse(1:Busnum*3,1:Busnum*3,Y.*sin(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3);
H=sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)*Y.*sin(mat_IvAngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)-sparse(1:Busnum*3,1:Busnum*3,Y.*sin(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3); N=-sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)*Y.*cos(mat_AngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)+sparse(1:Busnum*3,1:Busnum*3,Y.*cos(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3);
N=-sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)*Y.*cos(mat_IvAngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)+sparse(1:Busnum*3,1:Busnum*3,Y.*cos(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3); J=sparse(1:Busnum*3,1:Busnum*3,Y.*cos(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)+Y.*cos(mat_AngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3);
J=sparse(1:Busnum*3,1:Busnum*3,Y.*cos(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)+Y.*cos(mat_IvAngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3); L=sparse(1:Busnum*3,1:Busnum*3,Y.*sin(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)+Y.*sin(mat_AngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3);
L=sparse(1:Busnum*3,1:Busnum*3,Y.*sin(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)+Y.*sin(mat_IvAngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3); t1=[J,H;
t1=[J,L; L,N;
H,N; ];
]'; Jacob=t1;
Jacob=-t1;
end end