From ea40031257eb8b3f3b23cbc126c434712e8e0132 Mon Sep 17 00:00:00 2001 From: facat Date: Sat, 17 Aug 2013 21:40:53 +0800 Subject: [PATCH] =?UTF-8?q?=E5=8A=A0=E5=85=A5Armijo=E6=90=9C=E7=B4=A2,?= =?UTF-8?q?=E6=95=88=E6=9E=9C=E4=B8=8D=E5=A5=BD=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: facat --- Armijo.m | 17 +++++++++ objfun.m | 18 ++++++++++ run.m | 107 ++++++++++++++++++++++++++++--------------------------- 3 files changed, 89 insertions(+), 53 deletions(-) create mode 100644 Armijo.m create mode 100644 objfun.m diff --git a/Armijo.m b/Armijo.m new file mode 100644 index 0000000..3ebce72 --- /dev/null +++ b/Armijo.m @@ -0,0 +1,17 @@ +function [ output_args ] = Armijo(z,newwordParameter,W,SEVolt,SEVAngle,d,g ) +%Armijo搜索 + +a=0.9; +b=0.1; +output_args=a; +for m=1:100 + objfun( SEVolt,SEVAngle,W,z,newwordParameter ) + if objfun( SEVolt+a^m*d(1:length(SEVolt)),SEVAngle+a^m*d(length(SEVolt)+1:end),W,z,newwordParameter )<=objfun( SEVolt,SEVAngle,W,z,newwordParameter )+b*a^m*g'*d + output_args= a^m; + return; + end +end + + +end + diff --git a/objfun.m b/objfun.m new file mode 100644 index 0000000..b5d1b18 --- /dev/null +++ b/objfun.m @@ -0,0 +1,18 @@ +function [ output_args ] = objfun( SEVolt,SEVAngle,W,z,newwordParameter ) +%目标函数 +%% +lineI=newwordParameter.line.lineI; +lineJ=newwordParameter.line.lineJ; +lineR=newwordParameter.line.lineR; +lineX=newwordParameter.line.lineX; +lineB2=newwordParameter.line.lineB2; +SEBranchI=BranchI( SEVolt.*exp(1j*SEVAngle),lineI,lineJ,lineR,lineX );%复数支路电流 +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]; +output_args=(z-h)'*W*(z-h); +end + diff --git a/run.m b/run.m index ab9b9fd..21a8167 100644 --- a/run.m +++ b/run.m @@ -2,15 +2,9 @@ clear clc % yalmip('clear') addpath('.\Powerflow') -[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee4.dat', '0'); -%% 量测量 -% 电压 节点电流 支路电流 节点功率 支路功率 -%% -%% 状态量 -% 电压 相角 -%% +[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('S1047.dat', '0'); %% 开始生成量测量 -sigma=0.05;% 标准差 +sigma=0.09;% 标准差 %% 电压 %电压幅值 rVolt=Volt; %幅值 @@ -41,6 +35,18 @@ rBranchP= BranchP( cmpV,cmpBranchI,lineI,lineB2 ); mBranchP=rBranchP.*(normrnd(0,sigma,length(rBranchP),1)+1);%支路功率量测量 rBranchQ=BranchQ( cmpV,cmpBranchI,lineI,lineB2 ); mBranchQ=rBranchQ.*(normrnd(0,sigma,length(rBranchQ),1)+1);%支路功率量测量 +%% 变压器功率 +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)); +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); @@ -57,29 +63,18 @@ 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)); -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); %% 发电机注入功率 % 先找到只有发电机的节点 PDQDi=union(PDi,QDi); onlyPG=setdiff(PGi,PDQDi); onlyQG=setdiff(QGi,PDQDi); +%% 计算方差 +W=(diag([rVolt;rBranchP;rBranchQ;rTransP;rTransQ])*sigma).^2; +W=sparse(1:length(W),1:length(W),400,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); -%% Opti ToolBox -Busnum=length(Volt); %% save % save('mVolt','mVolt'); % save('mPG','mPG'); @@ -101,118 +96,116 @@ Busnum=length(Volt); %% 自己写的微分代码 % 以下都是Jacobi矩阵 SEVolt=sparse(ones(length(mVolt),1)); -SEVAngel=sparse(zeros(length(mVolt),1)); +SEVAngle=sparse(zeros(length(mVolt),1)); maxD=1000; Iteration=0; -while maxD>1e-3; +g=100; +while max(abs(g))>1e-5; +% while g>1e-5 % 电压 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))... + 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(SEVAngel(lineI)-SEVAngel(lineJ)) +lineB.*sin(SEVAngel(lineI)-SEVAngel(lineJ))... + 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(SEVAngel(lineI)-SEVAngel(lineJ)) -lineB.*cos(SEVAngel(lineI)-SEVAngel(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(SEVAngel(lineI)-SEVAngel(lineJ)) -lineB.*cos(SEVAngel(lineI)-SEVAngel(lineJ))... + lineG.*sin(SEVAngle(lineI)-SEVAngle(lineJ)) -lineB.*cos(SEVAngle(lineI)-SEVAngle(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))... + 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(SEVAngel(lineI)-SEVAngel(lineJ)) -lineB.*cos(SEVAngel(lineI)-SEVAngel(lineJ))... + 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(SEVAngel(lineI)-SEVAngel(lineJ)) +lineB.*sin(SEVAngel(lineI)-SEVAngel(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(SEVAngel(lineI)-SEVAngel(lineJ)) +lineB.*sin(SEVAngel(lineI)-SEVAngel(lineJ))... + lineG.*cos(SEVAngle(lineI)-SEVAngle(lineJ)) +lineB.*sin(SEVAngle(lineI)-SEVAngle(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))... + 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(SEVAngel(transI)-SEVAngel(transJ)) +transB.*sin(SEVAngel(transI)-SEVAngel(transJ))... + 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(SEVAngel(transI)-SEVAngel(transJ)) -transB.*cos(SEVAngel(transI)-SEVAngel(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(SEVAngel(transI)-SEVAngel(transJ)) -transB.*cos(SEVAngel(transI)-SEVAngel(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(SEVAngel(transI)-SEVAngel(transJ)) -transB.*cos(SEVAngel(transI)-SEVAngel(transJ))... + 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(SEVAngel(transI)-SEVAngel(transJ)) -transB.*cos(SEVAngel(transI)-SEVAngel(transJ))... + 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(SEVAngel(transI)-SEVAngel(transJ)) +transB.*sin(SEVAngel(transI)-SEVAngel(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(SEVAngel(transI)-SEVAngel(transJ)) +transB.*sin(SEVAngel(transI)-SEVAngel(transJ))... + transG.*cos(SEVAngle(transI)-SEVAngle(transJ)) +transB.*sin(SEVAngle(transI)-SEVAngle(transJ))... ) ... ,length(transI),length(mVolt));%变压器 %% 进入迭代 - % 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 ; dTPij_dVi+dTPij_dVj,dTPij_dThetai+dTPij_dThetaj; dTQij_dVi+dTQij_dVj,dTQij_dThetai+dTQij_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 ); - SETransP=TransPower( newwordParameter,SEVolt,SEVAngel ); - SETransQ=TransReactivePower( newwordParameter,SEVolt,SEVAngel ); + SEBranchI=BranchI( SEVolt.*exp(1j*SEVAngle),lineI,lineJ,lineR,lineX );%复数支路电流 + 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]; z=[mVolt;mBranchP;mBranchQ;mTransP;mTransQ]; - W=sparse(1:length(h),1:length(h),1/sigma.^2,length(h),length(h)); G=H'*W*H; g=-H'*W*(z-h); % 平衡节点相角恒定; @@ -221,11 +214,13 @@ while maxD>1e-3; G=G+sparse(length(mVolt)+Balance,length(mVolt)+Balance,1,length(mVolt)*2,length(mVolt)*2); g(length(mVolt)+Balance)=0; dX=G\-g; + dXStep=1; +% dXStep=Armijo(z,newwordParameter,W,SEVolt,SEVAngle,dX,g ); maxD=max(abs(dX)) % 更新变量 - SEVolt=SEVolt+dX(1:length(mVolt)); + SEVolt=SEVolt+dX(1:length(mVolt))*dXStep; Iteration=Iteration+1; - SEVAngel=SEVAngel+dX(length(mVolt)+1:end); + SEVAngle=SEVAngle+dX(length(mVolt)+1:end)*dXStep; end %% 输出结果 fprintf('迭代%d次\n',Iteration); @@ -233,5 +228,11 @@ 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 ) \ No newline at end of file +MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngle) +plotAndComparison( rVolt,rVAngel,SEVolt,SEVAngle ) +% 检查最优性条件 +if any(abs(g)>1e-5) + fprintf('最优性条件不满足\n') +else + fprintf('最优性条件满足\n') +end \ No newline at end of file