264 lines
11 KiB
Matlab
264 lines
11 KiB
Matlab
clear
|
||
clc
|
||
% yalmip('clear')
|
||
addpath('.\Powerflow')
|
||
[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee4.dat', '0');
|
||
%% 量测量
|
||
% 电压 节点电流 支路电流 节点功率 支路功率
|
||
%%
|
||
%% 状态量
|
||
% 电压 相角
|
||
%%
|
||
%% 开始生成量测量
|
||
sigma=0.05;% 标准差
|
||
%% 电压
|
||
%电压幅值
|
||
rVolt=Volt; %幅值
|
||
BalanceVolt=Volt(Balance);
|
||
mVolt=rVolt.*(normrnd(0,sigma,length(Volt),1)+1);%电压量测量
|
||
rVAngel=Vangle;
|
||
%% 电流
|
||
%注入电流
|
||
cmpY=Y.*exp(1j*sparse(r,c,Yangle,length(Y),length(Y)));%复数导纳矩阵
|
||
cmpV=Volt.*exp(1j*Vangle); %复数电压
|
||
cmpI=cmpY*cmpV;% 注入电流
|
||
rI=abs(cmpI); %注入电流量测量要的是电流幅值
|
||
mI=rI.*(normrnd(0,sigma,length(rI),1)+1);%电流量测量
|
||
%% 支路电流
|
||
% 支路电流
|
||
lineI=newwordParameter.line.lineI;
|
||
lineJ=newwordParameter.line.lineJ;
|
||
lineR=newwordParameter.line.lineR;
|
||
lineX=newwordParameter.line.lineX;
|
||
lineB2=newwordParameter.line.lineB2;
|
||
lineG=real(1./(lineR+1j*lineX));
|
||
lineB=imag(1./(lineR+1j*lineX));
|
||
% 处理线路电阻或电抗为0的情况，即消除NaN
|
||
% zerosLRXInd=union(find(abs(lineR)<1e-5),find(abs(lineX)<1e-5));
|
||
cmpBranchI=BranchI( cmpV,lineI,lineJ,lineR,lineX );%复数支路电流
|
||
rBranchI=abs(cmpBranchI);% 支路电流幅值
|
||
mBranchI=rBranchI.*(normrnd(0,sigma,length(rBranchI),1)+1);%支路电流量测量
|
||
%% 支路功率
|
||
rBranchP= BranchP( cmpV,cmpBranchI,lineI,lineB2 );
|
||
% rV=abs(cmpV);
|
||
% rA=angle(cmpV);
|
||
% a=rV(lineI).^2.*(lineG)-rV(lineI).*rV(lineJ).*(lineG.*cos(rA(lineI)-rA(lineJ) )+lineB.*sin(rA(lineI)-rA(lineJ) ) );
|
||
% b=-rV(lineI).^2.*(lineB+lineB2)-rV(lineI).*rV(lineJ).*(lineG.*sin(rA(lineI)-rA(lineJ) )-lineB.*cos(rA(lineI)-rA(lineJ) ) );
|
||
% rBranchP=rBranchP(abs(rBranchP)>1e-5);
|
||
mBranchP=rBranchP.*(normrnd(0,sigma,length(rBranchP),1)+1);%支路功率量测量
|
||
rBranchQ=BranchQ( cmpV,cmpBranchI,lineI,lineB2 );
|
||
% rBranchQ=rBranchQ(abs(rBranchQ)>1e-5);
|
||
mBranchQ=rBranchQ.*(normrnd(0,sigma,length(rBranchQ),1)+1);%支路功率量测量
|
||
% 处理线路电阻或电抗为0的情况，即消除NaN
|
||
% mBranchP=mBranchP(setdiff( 1:length(mBranchP),zerosLRXInd));
|
||
% mBranchQ=mBranchQ(setdiff( 1:length(mBranchQ),zerosLRXInd));
|
||
%% 注入功率
|
||
rPD=PD;
|
||
PDi=find(PD~=0);
|
||
rQD=QD;
|
||
QDi=find(QD~=0);
|
||
rPG=PG;
|
||
PGi=find(PG~=0);
|
||
rQG=QG;
|
||
QGi=find(QG~=0);
|
||
mPD=rPD.*(normrnd(0,sigma,length(rPD),1)+1);
|
||
mQD=rQD.*(normrnd(0,sigma,length(rQD),1)+1);
|
||
mPG=rPG.*(normrnd(0,sigma,length(rPG),1)+1);
|
||
mQG=rQG.*(normrnd(0,sigma,length(rQG),1)+1);
|
||
%% 0注入节点
|
||
zerosInjectionIndex=1:length(Volt);
|
||
zerosInjectionIndex=zerosInjectionIndex( ~(PD~=0|QD~=0|PG~=0|QG~=0) );
|
||
%% 变压器功率
|
||
transI=newwordParameter.trans.transI;
|
||
transJ=newwordParameter.trans.transJ;
|
||
transK=newwordParameter.trans.transK;
|
||
transR=newwordParameter.trans.transR;
|
||
transX=newwordParameter.trans.transX;
|
||
transG=real(1./(transR+1j*transX));
|
||
transB=imag(1./(transR+1j*transX));
|
||
% 处理线路电阻或电抗为0的情况，即消除NaN
|
||
% zerosTRXInd=union(find(abs(transR)<1e-5),find(abs(transX)<1e-5));
|
||
rTransP=TransPower( newwordParameter,cmpY,rVolt,rVAngel );
|
||
rTransQ=TransReactivePower( newwordParameter,cmpY,rVolt,rVAngel );
|
||
% rTransP=rTransP( setdiff( 1:length(rTransP),zerosTRXInd) );
|
||
% rTransQ=rTransQ( setdiff( 1:length(rTransQ),zerosTRXInd) );
|
||
mTransP=rTransP.*(normrnd(0,sigma,length(rTransP),1)+1);
|
||
mTransQ=rTransQ.*(normrnd(0,sigma,length(rTransQ),1)+1);
|
||
%% 发电机注入功率
|
||
% 先找到只有发电机的节点
|
||
PDQDi=union(PDi,QDi);
|
||
onlyPG=setdiff(PGi,PDQDi);
|
||
onlyQG=setdiff(QGi,PDQDi);
|
||
%% 冗余度计算
|
||
stateVarCount=2*length(Volt);
|
||
measurements=length(mVolt)+length(mBranchI)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG)+length(mTransP)+length(mTransQ);
|
||
fprintf('冗余度 %f\n',measurements/stateVarCount);
|
||
%% Opti ToolBox
|
||
Busnum=length(Volt);
|
||
%% save
|
||
% save('mVolt','mVolt');
|
||
% save('mPG','mPG');
|
||
% save('mQG','mQG');
|
||
% save('mBranchI','mBranchI');
|
||
% save('mBranchP','mBranchP');
|
||
% save('mBranchQ','mBranchQ');
|
||
% save('mTransP','mTransP');
|
||
% save('mTransQ','mTransQ');
|
||
%% load
|
||
% load('mVolt');
|
||
% load('mPG');
|
||
% load('mQG');
|
||
% load('mBranchI');
|
||
% load('mBranchP');
|
||
% load('mBranchQ');
|
||
% load('mTransP');
|
||
% load('mTransQ');
|
||
%% 自己写的微分代码
|
||
% 以下都是Jacobi矩阵
|
||
SEVolt=sparse(ones(length(mVolt),1));
|
||
SEVAngel=sparse(zeros(length(mVolt),1));
|
||
maxD=1000;
|
||
Iteration=0;
|
||
while maxD>1e-3;
|
||
% 电压
|
||
dV_dV=sparse(1:length(mVolt),1:length(mVolt),1,length(mVolt),length(mVolt));%电压量测量的微分
|
||
dV_dTyta=sparse(length(mVolt),length(mVolt));
|
||
% 线路支路
|
||
dLPij_dVi=sparse(1:length(lineI),lineI, ...
|
||
-SEVolt(lineJ).*( ...
|
||
lineG.*cos(SEVAngel(lineI)-SEVAngel(lineJ)) +lineB.*sin(SEVAngel(lineI)-SEVAngel(lineJ))...
|
||
)...
|
||
+2*(lineG).*SEVolt(lineI) ...
|
||
,length(lineI),length(mVolt));%线路的
|
||
dLPij_dVj=sparse(1:length(lineI),lineJ, ...
|
||
-SEVolt(lineI).*( ...
|
||
lineG.*cos(SEVAngel(lineI)-SEVAngel(lineJ)) +lineB.*sin(SEVAngel(lineI)-SEVAngel(lineJ))...
|
||
) ...
|
||
,length(lineI),length(mVolt));%线路的
|
||
dLPij_dThetai=sparse(1:length(lineI),lineI, ...
|
||
SEVolt(lineI).*SEVolt(lineJ).*( ...
|
||
lineG.*sin(SEVAngel(lineI)-SEVAngel(lineJ)) -lineB.*cos(SEVAngel(lineI)-SEVAngel(lineJ))...
|
||
)...
|
||
,length(lineI),length(mVolt));%线路的
|
||
dLPij_dThetaj=sparse(1:length(lineI),lineJ, ...
|
||
-SEVolt(lineI).*SEVolt(lineJ).*( ...
|
||
lineG.*sin(SEVAngel(lineI)-SEVAngel(lineJ)) -lineB.*cos(SEVAngel(lineI)-SEVAngel(lineJ))...
|
||
)...
|
||
,length(lineI),length(mVolt));%线路的
|
||
|
||
dLQij_dVi=sparse(1:length(lineI),lineI, ...
|
||
-SEVolt(lineJ).*( ...
|
||
lineG.*sin(SEVAngel(lineI)-SEVAngel(lineJ)) -lineB.*cos(SEVAngel(lineI)-SEVAngel(lineJ))...
|
||
)...
|
||
-2*(lineB+lineB2).*SEVolt(lineI) ...
|
||
,length(lineI),length(mVolt));%线路的
|
||
dLQij_dVj=sparse(1:length(lineI),lineJ, ...
|
||
-SEVolt(lineI).*( ...
|
||
lineG.*sin(SEVAngel(lineI)-SEVAngel(lineJ)) -lineB.*cos(SEVAngel(lineI)-SEVAngel(lineJ))...
|
||
) ...
|
||
,length(lineI),length(mVolt));%线路的
|
||
dLQij_dThetai=sparse(1:length(lineI),lineI, ...
|
||
-SEVolt(lineI).*SEVolt(lineJ).*( ...
|
||
lineG.*cos(SEVAngel(lineI)-SEVAngel(lineJ)) +lineB.*sin(SEVAngel(lineI)-SEVAngel(lineJ))...
|
||
) ...
|
||
,length(lineI),length(mVolt));%线路的
|
||
dLQij_dThetaj=sparse(1:length(lineI),lineJ, ...
|
||
SEVolt(lineI).*SEVolt(lineJ).*( ...
|
||
lineG.*cos(SEVAngel(lineI)-SEVAngel(lineJ)) +lineB.*sin(SEVAngel(lineI)-SEVAngel(lineJ))...
|
||
) ...
|
||
,length(lineI),length(mVolt));%线路的
|
||
% 变压器支路
|
||
dTPij_dVi=sparse(1:length(transI),transI, ...
|
||
-SEVolt(transJ)./transK.*( ...
|
||
transG.*cos(SEVAngel(transI)-SEVAngel(transJ)) +transB.*sin(SEVAngel(transI)-SEVAngel(transJ))...
|
||
)...
|
||
+2*transG.*SEVolt(transI)./(transK.^2) ...
|
||
,length(transI),length(mVolt));%变压器
|
||
dTPij_dVj=sparse(1:length(transI),transJ, ...
|
||
-SEVolt(transI)./transK.*( ...
|
||
transG.*cos(SEVAngel(transI)-SEVAngel(transJ)) +transB.*sin(SEVAngel(transI)-SEVAngel(transJ))...
|
||
) ...
|
||
,length(transI),length(mVolt));%变压器
|
||
dTPij_dThetai=sparse(1:length(transI),transI, ...
|
||
SEVolt(transI)./transK.*SEVolt(transJ).*( ...
|
||
transG.*sin(SEVAngel(transI)-SEVAngel(transJ)) -transB.*cos(SEVAngel(transI)-SEVAngel(transJ))...
|
||
)...
|
||
,length(transI),length(mVolt));%变压器
|
||
dTPij_dThetaj=sparse(1:length(transI),transJ, ...
|
||
-SEVolt(transI)./transK.*SEVolt(transJ).*( ...
|
||
transG.*sin(SEVAngel(transI)-SEVAngel(transJ)) -transB.*cos(SEVAngel(transI)-SEVAngel(transJ))...
|
||
)...
|
||
,length(transI),length(mVolt));%变压器
|
||
dTQij_dVi=sparse(1:length(transI),transI, ...
|
||
-SEVolt(transJ)./transK.*( ...
|
||
transG.*sin(SEVAngel(transI)-SEVAngel(transJ)) -transB.*cos(SEVAngel(transI)-SEVAngel(transJ))...
|
||
)...
|
||
-2*transB.*SEVolt(transI)./(transK.^2) ...
|
||
,length(transI),length(mVolt));%变压器
|
||
dTQij_dVj=sparse(1:length(transI),transJ, ...
|
||
-SEVolt(transI)./transK.*( ...
|
||
transG.*sin(SEVAngel(transI)-SEVAngel(transJ)) -transB.*cos(SEVAngel(transI)-SEVAngel(transJ))...
|
||
) ...
|
||
,length(transI),length(mVolt));%变压器
|
||
dTQij_dThetai=sparse(1:length(transI),transI, ...
|
||
-SEVolt(transI)./transK.*SEVolt(transJ).*( ...
|
||
transG.*cos(SEVAngel(transI)-SEVAngel(transJ)) +transB.*sin(SEVAngel(transI)-SEVAngel(transJ))...
|
||
) ...
|
||
,length(transI),length(mVolt));%变压器
|
||
dTQij_dThetaj=sparse(1:length(transI),transJ, ...
|
||
SEVolt(transI)./transK.*SEVolt(transJ).*( ...
|
||
transG.*cos(SEVAngel(transI)-SEVAngel(transJ)) +transB.*sin(SEVAngel(transI)-SEVAngel(transJ))...
|
||
) ...
|
||
,length(transI),length(mVolt));%变压器
|
||
% 处理线路电阻或电抗为0的情况，即消除NaN
|
||
% dLPij_dVi=dLPij_dVi(setdiff( 1:size(dLPij_dVi,1),zerosLRXInd),:);
|
||
% dLPij_dVj=dLPij_dVj(setdiff( 1:size(dLPij_dVj,1),zerosLRXInd),:);
|
||
% dLPij_dThetai=dLPij_dThetai(setdiff( 1:size(dLPij_dThetai,1),zerosLRXInd),:);
|
||
% dLPij_dThetaj=dLPij_dThetaj(setdiff( 1:size(dLPij_dThetaj,1),zerosLRXInd),:);
|
||
% dLQij_dVi=dLQij_dVi(setdiff( 1:size(dLQij_dVi,1),zerosLRXInd),:);
|
||
% dLQij_dVj=dLQij_dVj(setdiff( 1:size(dLQij_dVj,1),zerosLRXInd),:);
|
||
% dLQij_dThetai=dLQij_dThetai(setdiff( 1:size(dLQij_dThetai,1),zerosLRXInd),:);
|
||
% dLQij_dThetaj=dLQij_dThetaj(setdiff( 1:size(dLQij_dThetaj,1),zerosLRXInd),:);
|
||
% 对量测值做同样处理
|
||
|
||
%% 进入迭代
|
||
% H=[dV_dV ,dV_dTyta;
|
||
% dLPij_dVi+dLPij_dVj, dLPij_dThetai+dLPij_dThetaj;
|
||
% dLQij_dVi+dLQij_dVj, dLQij_dThetai+dLQij_dThetaj];%jacobi
|
||
|
||
H=[dV_dV,dV_dTyta;
|
||
dLPij_dVi+dLPij_dVj,dLPij_dThetai+dLPij_dThetaj ;
|
||
dLQij_dVi+dLQij_dVj,dLQij_dThetai+dLQij_dThetaj ];%jacobi
|
||
|
||
SEBranchI=BranchI( SEVolt.*exp(1j*SEVAngel),lineI,lineJ,lineR,lineX );%复数支路电流
|
||
SEBranchP=BranchP( SEVolt.*exp(1j*SEVAngel),SEBranchI,lineI,lineB2 );
|
||
SEBranchQ=BranchQ( SEVolt.*exp(1j*SEVAngel),SEBranchI,lineI,lineB2 );
|
||
% SEBranchP=SEBranchP(setdiff( 1:length(SEBranchP),zerosLRXInd));
|
||
% SEBranchQ=SEBranchQ(setdiff( 1:length(SEBranchQ),zerosLRXInd));
|
||
% h=[SEVolt;SEBranchP;SEBranchQ];
|
||
% z=[mVolt;mBranchP;mBranchQ];
|
||
h=[SEVolt;SEBranchP;SEBranchQ];
|
||
z=[mVolt;mBranchP;mBranchP];
|
||
W=sparse(1:length(h),1:length(h),1/sigma.^2,length(h),length(h));
|
||
G=H'*W*H;
|
||
g=-H'*W*(z-h);
|
||
% 平衡节点相角恒定;
|
||
% G(length(mVolt)+Balance,:)=0;
|
||
% G(:,length(mVolt)+Balance)=0;
|
||
% G=G+sparse(length(mVolt)+Balance,length(mVolt)+Balance,1,length(mVolt)*2,length(mVolt)*2);
|
||
% g(length(mVolt)+Balance)=0;
|
||
dX=G\-g;
|
||
maxD=max(abs(dX))
|
||
% 更新变量
|
||
% SEVolt=SEVolt+dX(1:length(mVolt));
|
||
Iteration=Iteration+1;
|
||
SEVAngel=SEVAngel+dX(1:end);
|
||
end
|
||
%% 输出结果
|
||
fprintf('迭代%d次\n',Iteration);
|
||
fval=full((z-h)'*W*(z-h));
|
||
fprintf('目标函数为:%f\n',fval);
|
||
fprintf('相对误差\n');
|
||
(abs(rVolt-double(SEVolt)))./(rVolt)
|
||
MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel)
|
||
plotAndComparison( rVolt,rVAngel,SEVolt,SEVAngel ) |