506 lines
24 KiB
Matlab
506 lines
24 KiB
Matlab
clear
|
|
clc
|
|
% yalmip('clear')
|
|
addpath('.\Powerflow')
|
|
[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('s1047.dat', '0');
|
|
%% 节点总数
|
|
busNum=length(Volt);
|
|
%% 变量总数
|
|
varNum=busNum*2;
|
|
%% 开始生成量测量
|
|
sigma=0.03;% 标准差
|
|
%% 电压
|
|
%电压幅值
|
|
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));
|
|
cmpBranchI=BranchI( cmpV,lineI,lineJ,lineR,lineX,lineB2 );%复数支路电流
|
|
rBranchI=abs(cmpBranchI);% 支路电流幅值
|
|
mBranchI=rBranchI.*(normrnd(0,sigma,length(rBranchI),1)+1);%支路电流量测量
|
|
%% 支路功率
|
|
rBranchP= BranchP( cmpV,cmpBranchI,lineI,lineB2 );
|
|
mBranchP=rBranchP.*(normrnd(0,sigma,length(rBranchP),1)+1);%支路功率量测量
|
|
%load('mBranchP');
|
|
%mBranchP=[0.284828008908987,-0.0282291720369099,-0.296295663000711]';
|
|
rBranchQ=BranchQ( cmpV,cmpBranchI,lineI,lineB2 );
|
|
mBranchQ=rBranchQ.*(normrnd(0,sigma,length(rBranchQ),1)+1);%支路功率量测量
|
|
%% 变压器功率
|
|
transI=newwordParameter.trans.transI;
|
|
transJ=newwordParameter.trans.transJ;
|
|
% newwordParameter.trans.transK=ones(length(newwordParameter.trans.transK),1);
|
|
transK=newwordParameter.trans.transK;
|
|
% transK=ones(length(transK),1);
|
|
transR=newwordParameter.trans.transR;
|
|
transX=newwordParameter.trans.transX;
|
|
transG=real(1./(transR+1j*transX));
|
|
transB=imag(1./(transR+1j*transX));
|
|
rTransP=TransPower( newwordParameter,rVolt,rVAngel );
|
|
rTransQ=TransReactivePower( newwordParameter,rVolt,rVAngel );
|
|
mTransP=rTransP.*(normrnd(0,sigma,length(rTransP),1)+1);
|
|
mTransQ=rTransQ.*(normrnd(0,sigma,length(rTransQ),1)+1);
|
|
%% 注入功率
|
|
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) );
|
|
% zerosInjectionIndex=zeros(0,0);
|
|
%% 发电机注入功率
|
|
% 先找到只有发电机的节点
|
|
PDQDi=union(PDi,QDi);
|
|
onlyPG=setdiff(PGi,PDQDi);
|
|
onlyQG=setdiff(QGi,PDQDi);
|
|
%% 计算方差
|
|
measureSigma=abs(([rVolt;rBranchP;rBranchQ;rTransP;rTransQ].*sigma));
|
|
measureSigma(measureSigma<1e-6)=mean(measureSigma(measureSigma>1e-6));
|
|
W=sparse(diag(1./measureSigma.^2)) ;
|
|
% W=sparse(1:length(W),1:length(W),1,length(W),length(W));
|
|
%% 冗余度计算
|
|
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);
|
|
%% 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');
|
|
%% 自己写的微分代码
|
|
% 初始化一些值
|
|
%SEVolt=sparse(ones(length(mVolt),1));
|
|
SEVolt=rVolt;
|
|
% SEVolt=rVolt;
|
|
% SEVolt(1:2:end)=1.2;
|
|
% load('SEVolt');
|
|
% load('SEVAngle');
|
|
SEVolt(Balance)=rVolt(Balance);
|
|
% SEVAngle=-0.1*sparse(ones(length(mVolt),1));
|
|
% SEVAngle(Balance)=0;
|
|
SEVAngle=rVAngel;
|
|
% SEVAngle=rVAngel;
|
|
% SEVolt2= load('SEVolt2');
|
|
% SEVolt2=SEVolt2.SEVolt;
|
|
% SEVAngle2= load('SEVAngle2');
|
|
% SEVAngle2=SEVAngle2.SEVAngle;
|
|
maxD=1000;
|
|
Iteration=0;
|
|
optimalCondition=100;
|
|
eps=1e-4;
|
|
mu=0;
|
|
v=2;
|
|
ojbFunDecrease=1000;% 目标函数下降
|
|
% 以下都是Jacobi矩阵
|
|
% while max(abs(g))>1e-5;
|
|
% while maxD>1e-5
|
|
while max(abs(optimalCondition))>eps
|
|
% 电压
|
|
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(SEVAngle(lineI)-SEVAngle(lineJ)) +lineB.*sin(SEVAngle(lineI)-SEVAngle(lineJ))...
|
|
)...
|
|
+2*(lineG).*SEVolt(lineI) ...
|
|
,length(lineI),length(mVolt));%线路的
|
|
dLPij_dVj=sparse(1:length(lineI),lineJ, ...
|
|
-SEVolt(lineI).*( ...
|
|
lineG.*cos(SEVAngle(lineI)-SEVAngle(lineJ)) +lineB.*sin(SEVAngle(lineI)-SEVAngle(lineJ))...
|
|
) ...
|
|
,length(lineI),length(mVolt));%线路的
|
|
dLPij_dThetai=sparse(1:length(lineI),lineI, ...
|
|
SEVolt(lineI).*SEVolt(lineJ).*( ...
|
|
lineG.*sin(SEVAngle(lineI)-SEVAngle(lineJ)) -lineB.*cos(SEVAngle(lineI)-SEVAngle(lineJ))...
|
|
)...
|
|
,length(lineI),length(mVolt));%线路的
|
|
dLPij_dThetaj=sparse(1:length(lineI),lineJ, ...
|
|
-SEVolt(lineI).*SEVolt(lineJ).*( ...
|
|
lineG.*sin(SEVAngle(lineI)-SEVAngle(lineJ)) -lineB.*cos(SEVAngle(lineI)-SEVAngle(lineJ))...
|
|
)...
|
|
,length(lineI),length(mVolt));%线路的
|
|
|
|
dLPij=[dLPij_dVi+dLPij_dVj,dLPij_dThetai+dLPij_dThetaj];
|
|
dLQij_dVi=sparse(1:length(lineI),lineI, ...
|
|
-SEVolt(lineJ).*( ...
|
|
lineG.*sin(SEVAngle(lineI)-SEVAngle(lineJ)) -lineB.*cos(SEVAngle(lineI)-SEVAngle(lineJ))...
|
|
)...
|
|
-2*(lineB+lineB2).*SEVolt(lineI) ...
|
|
,length(lineI),length(mVolt));%线路的
|
|
dLQij_dVj=sparse(1:length(lineI),lineJ, ...
|
|
-SEVolt(lineI).*( ...
|
|
lineG.*sin(SEVAngle(lineI)-SEVAngle(lineJ)) -lineB.*cos(SEVAngle(lineI)-SEVAngle(lineJ))...
|
|
) ...
|
|
,length(lineI),length(mVolt));%线路的
|
|
dLQij_dThetai=sparse(1:length(lineI),lineI, ...
|
|
-SEVolt(lineI).*SEVolt(lineJ).*( ...
|
|
lineG.*cos(SEVAngle(lineI)-SEVAngle(lineJ)) +lineB.*sin(SEVAngle(lineI)-SEVAngle(lineJ))...
|
|
) ...
|
|
,length(lineI),length(mVolt));%线路的
|
|
dLQij_dThetaj=sparse(1:length(lineI),lineJ, ...
|
|
SEVolt(lineI).*SEVolt(lineJ).*( ...
|
|
lineG.*cos(SEVAngle(lineI)-SEVAngle(lineJ)) +lineB.*sin(SEVAngle(lineI)-SEVAngle(lineJ))...
|
|
) ...
|
|
,length(lineI),length(mVolt));%线路的
|
|
dLQij=[dLQij_dVi+dLQij_dVj,dLQij_dThetai+dLQij_dThetaj];
|
|
% 变压器支路
|
|
dTPij_dVi=sparse(1:length(transI),transI, ...
|
|
-SEVolt(transJ)./transK.*( ...
|
|
transG.*cos(SEVAngle(transI)-SEVAngle(transJ)) +transB.*sin(SEVAngle(transI)-SEVAngle(transJ))...
|
|
)...
|
|
+2*transG.*SEVolt(transI)./(transK.^2) ...
|
|
,length(transI),length(mVolt));%变压器
|
|
dTPij_dVj=sparse(1:length(transI),transJ, ...
|
|
-SEVolt(transI)./transK.*( ...
|
|
transG.*cos(SEVAngle(transI)-SEVAngle(transJ)) +transB.*sin(SEVAngle(transI)-SEVAngle(transJ))...
|
|
) ...
|
|
,length(transI),length(mVolt));%变压器
|
|
dTPij_dThetai=sparse(1:length(transI),transI, ...
|
|
SEVolt(transI)./transK.*SEVolt(transJ).*( ...
|
|
transG.*sin(SEVAngle(transI)-SEVAngle(transJ)) -transB.*cos(SEVAngle(transI)-SEVAngle(transJ))...
|
|
)...
|
|
,length(transI),length(mVolt));%变压器
|
|
dTPij_dThetaj=sparse(1:length(transI),transJ, ...
|
|
-SEVolt(transI)./transK.*SEVolt(transJ).*( ...
|
|
transG.*sin(SEVAngle(transI)-SEVAngle(transJ)) -transB.*cos(SEVAngle(transI)-SEVAngle(transJ))...
|
|
)...
|
|
,length(transI),length(mVolt));%变压器
|
|
dTQij_dVi=sparse(1:length(transI),transI, ...
|
|
-SEVolt(transJ)./transK.*( ...
|
|
transG.*sin(SEVAngle(transI)-SEVAngle(transJ)) -transB.*cos(SEVAngle(transI)-SEVAngle(transJ))...
|
|
)...
|
|
-2*transB.*SEVolt(transI)./(transK.^2) ...
|
|
,length(transI),length(mVolt));%变压器
|
|
dTQij_dVj=sparse(1:length(transI),transJ, ...
|
|
-SEVolt(transI)./transK.*( ...
|
|
transG.*sin(SEVAngle(transI)-SEVAngle(transJ)) -transB.*cos(SEVAngle(transI)-SEVAngle(transJ))...
|
|
) ...
|
|
,length(transI),length(mVolt));%变压器
|
|
dTQij_dThetai=sparse(1:length(transI),transI, ...
|
|
-SEVolt(transI)./transK.*SEVolt(transJ).*( ...
|
|
transG.*cos(SEVAngle(transI)-SEVAngle(transJ)) +transB.*sin(SEVAngle(transI)-SEVAngle(transJ))...
|
|
) ...
|
|
,length(transI),length(mVolt));%变压器
|
|
dTQij_dThetaj=sparse(1:length(transI),transJ, ...
|
|
SEVolt(transI)./transK.*SEVolt(transJ).*( ...
|
|
transG.*cos(SEVAngle(transI)-SEVAngle(transJ)) +transB.*sin(SEVAngle(transI)-SEVAngle(transJ))...
|
|
) ...
|
|
,length(transI),length(mVolt));%变压器
|
|
%% 考虑等式约束
|
|
% 等式约束的Jacobi
|
|
% r=newwordParameter.r;
|
|
% c=newwordParameter.c;
|
|
% Yangle=newwordParameter.Yangle;
|
|
% VAngleIJ=sparse(r,c,SEVAngle(r)-SEVAngle(c) -Yangle,length(mVolt),length(mVolt)) ;
|
|
% YdotSin=Y.* ( spfun(@sin,VAngleIJ) );
|
|
% YdotCos=Y.* ( spfun (@cos, VAngleIJ ) );
|
|
% diag_Volt_YdotCos=diag(SEVolt)*YdotCos;
|
|
% diag_Volt_YdotSin=diag(SEVolt)*YdotSin;
|
|
% YdotCosVolt=YdotCos*Volt;
|
|
% YdotSinVolt=YdotSin*Volt;
|
|
% diag_Volt_YdotCosVolt=diag_Volt_YdotCos*Volt;
|
|
% diag_Volt_YdotSinVolt=diag_Volt_YdotSin*Volt;
|
|
% diag_YdotSinVolt_=diag(YdotSinVolt);
|
|
% diag_YdotCosVolt_=diag(YdotCosVolt);
|
|
% dPdTyta=diag_Volt_YdotSin*diag(SEVolt)-diag_YdotSinVolt_*diag(SEVolt); % 简化第三次
|
|
% dQdTyta=-diag_Volt_YdotCos*diag(SEVolt)+diag_YdotCosVolt_*diag(SEVolt);%dQ/dThyta
|
|
% dPdV=diag_YdotCosVolt_+diag_Volt_YdotCos;%dP/dV
|
|
% dQdV=diag_YdotSinVolt_+diag_Volt_YdotSin;%dQ/dV
|
|
% % C 是等式约束 c 的Jacobi
|
|
% C=[dPdV dPdTyta;
|
|
% dQdV dQdTyta];
|
|
% C=C(zerosInjectionIndex,:);
|
|
% % 形成等式约束 c
|
|
% nodeP=diag_Volt_YdotCosVolt;
|
|
% nodeQ=diag_Volt_YdotSinVolt;
|
|
% nodePQ=[nodeP;nodeQ];
|
|
% c=nodePQ(zerosInjectionIndex);
|
|
%% 进入迭代
|
|
SEBranchI=BranchI( SEVolt.*exp(1j*SEVAngle),lineI,lineJ,lineR,lineX,lineB2 );%复数支路电流
|
|
SEBranchP=BranchP( SEVolt.*exp(1j*SEVAngle),SEBranchI,lineI,lineB2 );
|
|
SEBranchQ=BranchQ( SEVolt.*exp(1j*SEVAngle),SEBranchI,lineI,lineB2 );
|
|
SETransP=TransPower( newwordParameter,SEVolt,SEVAngle );
|
|
SETransQ=TransReactivePower( newwordParameter,SEVolt,SEVAngle );
|
|
h=[SEVolt;SEBranchP;SEBranchQ;SETransP;SETransQ];
|
|
|
|
% load('mBranchP');
|
|
% mBranchQ=mBranchP;
|
|
% mTransP=mBranchP(1:length(mTransP));
|
|
% mTransQ=mBranchP(1:length(mTransQ));
|
|
% mVolt=mBranchP(1:length(mVolt));
|
|
%mBranchP=[0.284828008908987,-0.0282291720369099,-0.296295663000711]';
|
|
% load('mTransP');
|
|
% mTransP=mTransP';
|
|
% mTransQ=mTransP;
|
|
z=[mVolt;mBranchP;mBranchQ;mTransP;mTransQ];
|
|
Whz=W*(h-z);%形成海森矩阵时使用
|
|
|
|
%% 一阶导数
|
|
H=[dV_dV,dV_dTyta;
|
|
dLPij_dVi+dLPij_dVj,dLPij_dThetai+dLPij_dThetaj ;
|
|
dLQij_dVi+dLQij_dVj,dLQij_dThetai+dLQij_dThetaj ;
|
|
dTPij_dVi+dTPij_dVj,dTPij_dThetai+dTPij_dThetaj;
|
|
dTQij_dVi+dTQij_dVj,dTQij_dThetai+dTQij_dThetaj];%jacobi
|
|
%% 二阶导数
|
|
% 线路二阶导数
|
|
% 对电压求导
|
|
d2LPij_dVi2=sparse(lineI,lineI, ...
|
|
2*lineG.*Whz(length(mVolt)+1:length(mVolt)+length(lineI)),...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
d2LPij_dVidVj=sparse(lineI,lineJ, ...
|
|
Whz(length(mVolt)+1:length(mVolt)+length(lineI))...
|
|
.* ( -(cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG + sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB)), ...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
d2LPij_dVj2=0;%线路有功海森矩阵
|
|
% 对相角求导
|
|
d2LPij_dTi2=sparse(lineI+busNum,lineI+busNum, ...
|
|
Whz(length(mVolt)+1:length(mVolt)+length(lineI))...
|
|
.*SEVolt(lineI).*SEVolt(lineJ).*( cos( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineG + sin( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineB), ...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
d2LPij_dTidTj=sparse(lineI+busNum,lineJ+busNum, ...
|
|
Whz(length(mVolt)+1:length(mVolt)+length(lineI))...
|
|
.*(-SEVolt(lineI).*SEVolt(lineJ).*( cos( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineG + sin( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineB) ), ...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
d2LPij_dTj2=sparse(lineJ+busNum,lineJ+busNum, ...
|
|
Whz(length(mVolt)+1:length(mVolt)+length(lineI))...
|
|
.*SEVolt(lineI).*SEVolt(lineJ).*( cos( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineG + sin( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineB), ...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
%%
|
|
% 对电压相角求导
|
|
d2LPij_dVidTi=sparse(lineI,lineI+busNum, ...
|
|
Whz(length(mVolt)+1:length(mVolt)+length(lineI))...
|
|
.* (SEVolt(lineJ).* (sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG - cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
d2LPij_dVidTj=sparse(lineI,lineJ+busNum, ...
|
|
Whz(length(mVolt)+1:length(mVolt)+length(lineI))...
|
|
.* (-SEVolt(lineJ).* (sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG - cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
d2LPij_dVjdTi=sparse(lineJ,lineI+busNum, ...
|
|
Whz(length(mVolt)+1:length(mVolt)+length(lineI))...
|
|
.* (SEVolt(lineI).* (sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG - cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
d2LPij_dVjdTj=sparse(lineJ,lineJ+busNum, ...
|
|
Whz(length(mVolt)+1:length(mVolt)+length(lineI))...
|
|
.* (-SEVolt(lineI).* (sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG - cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
%%
|
|
% 对电压求导
|
|
d2LQij_dVi2=sparse(lineI,lineI, ...
|
|
-2*(lineB+lineB2).*Whz(length(mVolt)+length(lineI)+1:length(mVolt)+length(lineI)*2),...
|
|
varNum,varNum);%线路无功海森矩阵
|
|
d2LQij_dVidVj=sparse(lineI,lineJ, ...
|
|
Whz(length(mVolt)+length(lineI)+1:length(mVolt)+length(lineI)*2) ...
|
|
.*( -(sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG - cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB)), ...
|
|
varNum,varNum);%线路无功海森矩阵
|
|
d2LQij_dVj2=0;%线路无功海森矩阵
|
|
% 对相角求导
|
|
d2LQij_dTi2=sparse(lineI+busNum,lineI+busNum, ...
|
|
Whz(length(mVolt)+length(lineI)+1:length(mVolt)+length(lineI)*2)...
|
|
.*SEVolt(lineI).*SEVolt(lineJ).*( sin( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineG - cos( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineB), ...
|
|
varNum,varNum);%线路无功海森矩阵
|
|
d2LQij_dTidTj=sparse(lineI+busNum,lineJ+busNum, ...
|
|
Whz(length(mVolt)+length(lineI)+1:length(mVolt)+length(lineI)*2)...
|
|
.*(-SEVolt(lineI).*SEVolt(lineJ).*( sin( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineG - cos( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineB)), ...
|
|
varNum,varNum);%线路无功海森矩阵
|
|
d2LQij_dTj2=sparse(lineJ+busNum,lineJ+busNum, ...
|
|
Whz(length(mVolt)+length(lineI)+1:length(mVolt)+length(lineI)*2)...
|
|
.*SEVolt(lineI).*SEVolt(lineJ).*( sin( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineG - cos( SEVAngle(lineI) -SEVAngle(lineJ) ).*lineB), ...
|
|
varNum,varNum);%线路无功海森矩阵
|
|
% 对电压相角求导
|
|
d2LQij_dVidTi=sparse(lineI,lineI+busNum, ...
|
|
Whz(length(mVolt)+length(lineI)+1:length(mVolt)+length(lineI)*2)...
|
|
.* (-SEVolt(lineJ).* (cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG + sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
d2LQij_dVidTj=sparse(lineI,lineJ+busNum, ...
|
|
Whz(length(mVolt)+length(lineI)+1:length(mVolt)+length(lineI)*2)...
|
|
.* (SEVolt(lineJ).* (cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG + sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
d2LQij_dVjdTi=sparse(lineJ,lineI+busNum, ...
|
|
Whz(length(mVolt)+length(lineI)+1:length(mVolt)+length(lineI)*2)...
|
|
.* (-SEVolt(lineI).* (cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG + sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
d2LQij_dVjdTj=sparse(lineJ,lineJ+busNum, ...
|
|
Whz(length(mVolt)+length(lineI)+1:length(mVolt)+length(lineI)*2)...
|
|
.* (SEVolt(lineI).* (cos(SEVAngle(lineI) -SEVAngle(lineJ)).*lineG + sin(SEVAngle(lineI) -SEVAngle(lineJ)).*lineB) ), ...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
%% 变压器二阶导数
|
|
% 对电压求导
|
|
d2TPij_dVi2=sparse(transI,transI, ...
|
|
2*transG./transK./transK.*Whz(length(mVolt)+2*length(lineI)+1:length(mVolt)+2*length(lineI)+length(transI)),...
|
|
varNum,varNum);%变压器有功海森矩阵
|
|
d2TPij_dVidVj=sparse(transI,transJ, ...
|
|
Whz(length(mVolt)+2*length(lineI)+1:length(mVolt)+2*length(lineI)+length(transI))...
|
|
.* ( -(cos(SEVAngle(transI) -SEVAngle(transJ)).*transG + sin(SEVAngle(transI) -SEVAngle(transJ)).*transB)./transK), ...
|
|
varNum,varNum);%变压器有功海森矩阵
|
|
d2TPij_dVj2=0;%变压器有功海森矩阵
|
|
% 对相角求导
|
|
d2TPij_dTi2=sparse(transI+busNum,transI+busNum, ...
|
|
Whz(length(mVolt)+2*length(lineI)+1:length(mVolt)+2*length(lineI)+length(transI))...
|
|
.*SEVolt(transI).*SEVolt(transJ)./transK.*( cos( SEVAngle(transI) -SEVAngle(transJ) ).*transG + sin( SEVAngle(transI) -SEVAngle(transJ) ).*transB), ...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
d2TPij_dTidTj=sparse(transI+busNum,transJ+busNum, ...
|
|
Whz(length(mVolt)+2*length(lineI)+1:length(mVolt)+2*length(lineI)+length(transI))...
|
|
.*(-SEVolt(transI).*SEVolt(transJ)./transK.*( cos( SEVAngle(transI) -SEVAngle(transJ) ).*transG + sin( SEVAngle(transI) -SEVAngle(transJ) ).*transB) ), ...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
d2TPij_dTj2=sparse(transJ+busNum,transJ+busNum, ...
|
|
Whz(length(mVolt)+2*length(lineI)+1:length(mVolt)+2*length(lineI)+length(transI))...
|
|
.*SEVolt(transI).*SEVolt(transJ)./transK.*( cos( SEVAngle(transI) -SEVAngle(transJ) ).*transG + sin( SEVAngle(transI) -SEVAngle(transJ) ).*transB), ...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
% 对电压相角求导
|
|
d2TPij_dVidTi=sparse(transI,transI+busNum, ...
|
|
Whz(length(mVolt)+length(lineI)*2+1:length(mVolt)+length(lineI)*2+length(transI)) ...
|
|
.* (SEVolt(transJ)./transK.* (sin(SEVAngle(transI) -SEVAngle(transJ)).*transG - cos(SEVAngle(transI) -SEVAngle(transJ)).*transB) ), ...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
d2TPij_dVidTj=sparse(transI,transJ+busNum, ...
|
|
Whz(length(mVolt)+length(lineI)*2+1:length(mVolt)+length(lineI)*2+length(transI)) ...
|
|
.* (-SEVolt(transJ)./transK.* (sin(SEVAngle(transI) -SEVAngle(transJ)).*transG - cos(SEVAngle(transI) -SEVAngle(transJ)).*transB) ), ...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
d2TPij_dVjdTi=sparse(transJ,transI+busNum, ...
|
|
Whz(length(mVolt)+length(lineI)*2+1:length(mVolt)+length(lineI)*2+length(transI)) ...
|
|
.* (SEVolt(transI)./transK.* (sin(SEVAngle(transI) -SEVAngle(transJ)).*transG - cos(SEVAngle(transI) -SEVAngle(transJ)).*transB) ), ...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
d2TPij_dVjdTj=sparse(transJ,transJ+busNum, ...
|
|
Whz(length(mVolt)+length(lineI)*2+1:length(mVolt)+length(lineI)*2+length(transI)) ...
|
|
.* (-SEVolt(transI)./transK.* (sin(SEVAngle(transI) -SEVAngle(transJ)).*transG - cos(SEVAngle(transI) -SEVAngle(transJ)).*transB) ), ...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
% 对电压求导
|
|
d2TQij_dVi2=sparse(transI,transI, ...
|
|
-2*transB./transK./transK.*Whz(length(mVolt)+length(lineI)*2+length(transI)+1:length(mVolt)+length(lineI)*2+length(transI)*2),...
|
|
varNum,varNum);%变压器无功海森矩阵
|
|
d2TQij_dVidVj=sparse(transI,transJ, ...
|
|
Whz(length(mVolt)+length(lineI)*2+length(transI)+1:length(mVolt)+length(lineI)*2+length(transI)*2) ...
|
|
.*( -(sin(SEVAngle(transI) -SEVAngle(transJ)).*transG - cos(SEVAngle(transI) -SEVAngle(transJ)).*transB)./transK), ...
|
|
varNum,varNum);%变压器无功海森矩阵
|
|
d2TQij_dVj2=0;%变压器无功海森矩阵
|
|
% 对相角求导
|
|
d2TQij_dTi2=sparse(transI+busNum,transI+busNum, ...
|
|
Whz(length(mVolt)+length(lineI)*2+length(transI)+1:length(mVolt)+length(lineI)*2+length(transI)*2)...
|
|
.*SEVolt(transI).*SEVolt(transJ)./transK.*( sin( SEVAngle(transI) -SEVAngle(transJ) ).*transG - cos( SEVAngle(transI) -SEVAngle(transJ) ).*transB), ...
|
|
varNum,varNum);%变压器无功海森矩阵
|
|
d2TQij_dTidTj=sparse(transI+busNum,transJ+busNum, ...
|
|
Whz(length(mVolt)+length(lineI)*2+length(transI)+1:length(mVolt)+length(lineI)*2+length(transI)*2)...
|
|
.*(-SEVolt(transI).*SEVolt(transJ)./transK.*( sin( SEVAngle(transI) -SEVAngle(transJ) ).*transG - cos( SEVAngle(transI) -SEVAngle(transJ) ).*transB)), ...
|
|
varNum,varNum);%变压器无功海森矩阵
|
|
d2TQij_dTj2=sparse(transJ+busNum,transJ+busNum, ...
|
|
Whz(length(mVolt)+length(lineI)*2+length(transI)+1:length(mVolt)+length(lineI)*2+length(transI)*2)...
|
|
.*SEVolt(transI).*SEVolt(transJ)./transK.*( sin( SEVAngle(transI) -SEVAngle(transJ) ).*transG - cos( SEVAngle(transI) -SEVAngle(transJ) ).*transB), ...
|
|
varNum,varNum);%变压器无功海森矩阵
|
|
% 对电压相角求导
|
|
d2TQij_dVidTi=sparse(transI,transI+busNum, ...
|
|
Whz(length(mVolt)+length(lineI)*2+length(transI)+1:length(mVolt)+length(lineI)*2+length(transI)*2)...
|
|
.* (-SEVolt(transJ)./transK.* (cos(SEVAngle(transI) -SEVAngle(transJ)).*transG + sin(SEVAngle(transI) -SEVAngle(transJ)).*transB) ), ...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
d2TQij_dVidTj=sparse(transI,transJ+busNum, ...
|
|
Whz(length(mVolt)+length(lineI)*2+length(transI)+1:length(mVolt)+length(lineI)*2+length(transI)*2)...
|
|
.* (SEVolt(transJ)./transK.* (cos(SEVAngle(transI) -SEVAngle(transJ)).*transG + sin(SEVAngle(transI) -SEVAngle(transJ)).*transB) ), ...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
d2TQij_dVjdTi=sparse(transJ,transI+busNum, ...
|
|
Whz(length(mVolt)+length(lineI)*2+length(transI)+1:length(mVolt)+length(lineI)*2+length(transI)*2)...
|
|
.* (-SEVolt(transI)./transK.* (cos(SEVAngle(transI) -SEVAngle(transJ)).*transG + sin(SEVAngle(transI) -SEVAngle(transJ)).*transB) ), ...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
d2TQij_dVjdTj=sparse(transJ,transJ+busNum, ...
|
|
Whz(length(mVolt)+length(lineI)*2+length(transI)+1:length(mVolt)+length(lineI)*2+length(transI)*2)...
|
|
.* (SEVolt(transI)./transK.* (cos(SEVAngle(transI) -SEVAngle(transJ)).*transG + sin(SEVAngle(transI) -SEVAngle(transJ)).*transB) ), ...
|
|
varNum,varNum);%线路有功海森矩阵
|
|
%% 形成大的海森矩阵
|
|
%海森矩阵对角，即
|
|
% [ d2h_dV2 0
|
|
% 0 d2h_dT2]
|
|
hessian14=d2LPij_dVi2+d2LPij_dVidVj+d2LPij_dVidVj'+d2LPij_dVj2 ...
|
|
+d2LPij_dTi2+d2LPij_dTidTj+d2LPij_dTidTj'+d2LPij_dTj2 ...
|
|
+d2LQij_dVi2+d2LQij_dVidVj+d2LQij_dVidVj'+d2LQij_dVj2 ...
|
|
+d2LQij_dTi2+d2LQij_dTidTj+d2LQij_dTidTj'+d2LQij_dTj2 ...;
|
|
+d2TPij_dVi2+d2TPij_dVidVj+d2TPij_dVidVj'+d2TPij_dVj2 ...
|
|
+d2TPij_dTi2+d2TPij_dTidTj+d2TPij_dTidTj'+d2TPij_dTj2 ...
|
|
+d2TQij_dVi2+d2TQij_dVidVj+d2TQij_dVidVj'+d2TQij_dVj2 ...
|
|
+d2TQij_dTi2+d2TQij_dTidTj+d2TQij_dTidTj'+d2TQij_dTj2;
|
|
hessian2=d2LPij_dVidTi+d2LPij_dVidTj+d2LPij_dVjdTi+d2LPij_dVjdTj ...
|
|
+d2LQij_dVidTi+d2LQij_dVidTj+d2LQij_dVjdTi+d2LQij_dVjdTj ...
|
|
+d2TPij_dVidTi+d2TPij_dVidTj+d2TPij_dVjdTi+d2TPij_dVjdTj ...
|
|
+d2TQij_dVidTi+d2TQij_dVidTj+d2TQij_dVjdTi+d2TQij_dVjdTj;
|
|
hessian3=hessian2';
|
|
if Iteration<=-1
|
|
G=H'*W*H;
|
|
else
|
|
G=H'*W*H+hessian14+hessian2+hessian3;
|
|
end
|
|
g=H'*Whz;
|
|
|
|
% hVolt=H(1:busNum,:)'*H(1:busNum,:);
|
|
a=G;
|
|
b=-g;
|
|
% 平衡节点相角恒定
|
|
a(length(mVolt)+Balance,:)=0;
|
|
a(:,length(mVolt)+Balance)=0;
|
|
a=a+sparse(length(mVolt)+Balance,length(mVolt)+Balance,1,size(a,1),size(a,1));
|
|
b(length(mVolt)+Balance)=0;
|
|
% 平衡节点电压恒定;
|
|
a(Balance,:)=0;
|
|
a(:,Balance)=0;
|
|
a=a+sparse(Balance,Balance,1,size(a,1),size(a,1));
|
|
b(Balance)=0;
|
|
dX=a\b;
|
|
dXStep=1;
|
|
% dXStep=Armijo(z,newwordParameter,W,SEVolt,SEVAngle,dX,g );
|
|
maxD=max(abs(dX(1:length(mVolt))))
|
|
fprintf('max abs g:%f\n',full(max(abs(g))));
|
|
SEVolt=SEVolt+dX(1:length(mVolt))*dXStep;
|
|
SEVAngle=SEVAngle+dX(length(mVolt)+1:end)*dXStep;
|
|
Iteration=Iteration+1;
|
|
optimalCondition=-g;
|
|
optimalCondition(Balance)=0;
|
|
optimalCondition(Balance+length(mVolt))=0;
|
|
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,SEVAngle)
|
|
plotAndComparison( rVolt,rVAngel,SEVolt,SEVAngle )
|
|
% 检查最优性条件
|
|
|
|
if any(abs(optimalCondition)>eps)
|
|
fprintf('最优性条件不满足\n')
|
|
else
|
|
fprintf('最优性条件满足\n')
|
|
end |