diff --git a/20PD0.mat b/20PD0.mat new file mode 100644 index 0000000..fd7c1df Binary files /dev/null and b/20PD0.mat differ diff --git a/20QD0.mat b/20QD0.mat new file mode 100644 index 0000000..d8b2e1a Binary files /dev/null and b/20QD0.mat differ diff --git a/@Opti/Opti.asv b/@Opti/Opti.asv new file mode 100644 index 0000000..2aea53b --- /dev/null +++ b/@Opti/Opti.asv @@ -0,0 +1,34 @@ +classdef Opti + %OPTI Summary of this class goes here + % Detailed explanation goes here + + properties + PD0=NaN; + QD0=NaN; + Volt0=NaN; + wPD=NaN; + wQD=NaN; + wVolt=NaN; + PDi=NaN; + QDi=NaN; + Y=NaN; + Angel=NaN; + r=NaN; + c=NaN; + PG=NaN; + QG=NaN; + cl=NaN;%opti toolbox大于小于等于 + cu=NaN;%等式、不等式约束右边的值 + Balance=NaN; + rPD=NaN; + rQD=NaN; + m + end + + methods + function this=Opti() + end + end + +end + diff --git a/@Opti/Opti.m b/@Opti/Opti.m index 0147958..9e36a51 100644 --- a/@Opti/Opti.m +++ b/@Opti/Opti.m @@ -9,6 +9,7 @@ classdef Opti wPD=NaN; wQD=NaN; wVolt=NaN; + wLoadCurrent=NaN; PDi=NaN; QDi=NaN; Y=NaN; @@ -22,6 +23,7 @@ classdef Opti Balance=NaN; rPD=NaN; rQD=NaN; + mLoadCurrent=NaN; end methods diff --git a/@Opti/equ.m b/@Opti/equ.m index 0101ed6..0ad51e0 100644 --- a/@Opti/equ.m +++ b/@Opti/equ.m @@ -23,7 +23,16 @@ busNum=length(Volt0); VMatrix=sparse(r,c,VAngel(r)-VAngel(c)-Angel,busNum,busNum); dP=PG-PD_-diag(Volt)*(Y.*cos(VMatrix))*Volt; dQ=QG-QD_-diag(Volt)*(Y.*spfun(@sin,VMatrix))*Volt; -output_args=[dP;dQ;VAngel(Balance)]; +output_args=[dP;dQ]; +%% 使用零注入begin +% zeroInj=setdiff(1:busNum,union(Balance,union(PDi,QDi))); +% output_args=[output_args(zeroInj);output_args(length(PDi)+zeroInj)]; +%% 使用零注入end + +%电压恒定 +output_args=[output_args;Volt(Balance)-1]; +%相角恒定 +output_args=[output_args;VAngel(Balance)]; this.cu=zeros(length(output_args),1); this.cl=this.cu; %% 开始增加不等式约束-电压\PD\QD @@ -45,5 +54,6 @@ this.cl=[this.cl;0.8*rQD]; output_args=full(output_args); this.cu=full(this.cu); this.cl=full(this.cl); + end diff --git a/@Opti/init.m b/@Opti/init.m index 3ebf24c..01d996f 100644 --- a/@Opti/init.m +++ b/@Opti/init.m @@ -1,4 +1,4 @@ -function [ this ] = init( this,Volt0,PDi,QDi,wPD,wQD,PD0,QD0,Y,Angel,r,c,PG,QG,Balance,rPD,rQD ) +function [ this ] = init( this,Volt0,PDi,QDi,wPD,wQD,wVolt,wLoadCurrent,PD0,QD0,rPD,rQD,Y,Angel,r,c,PG,QG,Balance,mLoadCurrent) %INIT Summary of this function goes here % Detailed explanation goes here this.PDi=PDi; @@ -17,5 +17,8 @@ this.Balance=Balance; this.Volt0=Volt0; this.rPD=rPD; this.rQD=rQD; +this.wVolt=wVolt; +this.wLoadCurrent=wLoadCurrent; +this.mLoadCurrent=mLoadCurrent; end diff --git a/@Opti/obj.asv b/@Opti/obj.asv index a1efcc8..5a5d070 100644 --- a/@Opti/obj.asv +++ b/@Opti/obj.asv @@ -1,15 +1,22 @@ -function [ output_args ] = obj( this ) +function [ output_args ] = obj( this,x ) %OBJ Summary of this function goes here % Detailed explanation goes here +PDi=this.PDi; +QDi=this.QDi; +PD=x(1:length(PDi)); +QD=x(length(PDi)+1:length(PDi)+length(QDi)); wPD=this.wPD; wQD=this.wQD; -PD=this.PD; -QD=this.QD; -wVolt=t -t4=wPD.*(PD(PD0~=0)-PD0(PD0~=0)).^2; -t5=wQD.*(QD(QD0~=0)-QD0(QD0~=0)).^2; -t6=wVolt.*(Volt-Volt0').^2; -out_arg=sum(t4)+sum(t5)+sum(t6); - +% wVolt=this.wVolt; +PD0=this.PD0; +QD0=this.QD0; +t4=(wPD(PDi).*(PD-PD0(PDi))).^2; +t5=(wQD(QDi).*(QD-QD0(QDi))).^2; +Volt0=this.Volt0; +SEVolt=x(length(PDi)+length(QDi)+1:length(PDi)+length(QDi)+length(Volt0)); +wVolt +t6=(wVolt.*(SEVolt-Volt0)).^2; +output_args=sum(t4)+sum(t5)+sum(t6); +output_args=full(output_args); end diff --git a/@Opti/obj.m b/@Opti/obj.m index 67d97e0..8597e41 100644 --- a/@Opti/obj.m +++ b/@Opti/obj.m @@ -12,7 +12,19 @@ PD0=this.PD0; QD0=this.QD0; t4=(wPD(PDi).*(PD-PD0(PDi))).^2; t5=(wQD(QDi).*(QD-QD0(QDi))).^2; -output_args=sum(t4)+sum(t5); +Volt0=this.Volt0; +SEVolt=x(length(PDi)+length(QDi)+1:length(PDi)+length(QDi)+length(Volt0)); +wVolt=this.wVolt; +t6=(wVolt.*(SEVolt-Volt0)).^2; +Balance=this.Balance; +t6(Balance)=0;%不计算平衡节点的 +% 电流 +SEVAngel=x(length(PDi)+length(QDi)+length(SEVolt)+1:end); +loadCurrent=LoadCurrent( SEVolt,SEVAngel,PD0(PDi),QD0(QDi),PDi,QDi ); +wLoadCurrent=this.wLoadCurrent; +mLoadCurrent=this.mLoadCurrent; +t7=(wLoadCurrent.*(loadCurrent-mLoadCurrent)).^2; +output_args=sum(t4)+sum(t5)+sum(t7); output_args=full(output_args); end diff --git a/LoadCurrent.asv b/LoadCurrent.asv new file mode 100644 index 0000000..924bfc1 --- /dev/null +++ b/LoadCurrent.asv @@ -0,0 +1,6 @@ +function [ output_args ] = LoadCurrent( Volt,PD,QD,PDi,QDi ) +%% 计算负荷电流 + + +end + diff --git a/LoadCurrent.m b/LoadCurrent.m new file mode 100644 index 0000000..7f74a65 --- /dev/null +++ b/LoadCurrent.m @@ -0,0 +1,8 @@ +function [ output_args ] = LoadCurrent( Volt,VAngel,PD,QD,PDi,QDi ) +%% 计算负荷电流 +cmpVolt=Volt.*exp(1j*VAngel); +t1=abs(conj((PD+1j*QD)./cmpVolt(PDi))); +output_args=t1; + +end + diff --git a/MaxDeviation.m b/MaxDeviation.m index 2304770..3c7a897 100644 --- a/MaxDeviation.m +++ b/MaxDeviation.m @@ -1,6 +1,8 @@ function [output_arg]=MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel,rPD,rQD,PD,QD) -t1=[rVolt;rVAngel;rPD;rQD]; -t2=[SEVolt;SEVAngel;PD;QD]; +% t1=[rVolt;rVAngel;rPD;rQD]; +% t2=[SEVolt;SEVAngel;PD;QD]; +t1=[rPD;rQD]; +t2=[PD;QD]; t3=(t1(t1~=0)-double(t2(t1~=0)))./t1(t1~=0); t4=abs(t3); diff --git a/OPF_Init.m b/OPF_Init.m index 805f491..4e66492 100644 --- a/OPF_Init.m +++ b/OPF_Init.m @@ -1,4 +1,4 @@ -function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,wVolt,PD,PD0,QD,randPDind,Loadi,notLoadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD) +function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,wVolt,wLoadCurrent,PD,PD0,QD,randPDind,Loadi,notLoadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD) Loadi=find(QD~=0 | PD~=0); PDi=find(PD~=0); QDi=find(QD~=0); @@ -34,8 +34,11 @@ wQG=1*ones(size(PVi,1),1); %randPDind=randInt(1:10); randPDind=0; wPD=1/0.05*ones(length(PD),1); +% wPD(1:2:end)=0; wQD=1/0.05*ones(length(QD),1); +% wQD(1:2:end)=0; wVolt=1/0.05*ones(Busnum,1); +wLoadCurrent=1/0.05*ones(length(PDi),1); %wD(randPDind)=0;%一些负荷不约束 %wD(7)=0; % wD(11)=0; diff --git a/PD0.mat b/PD0.mat index 9d99aa7..39a84c6 100644 Binary files a/PD0.mat and b/PD0.mat differ diff --git a/QD0.mat b/QD0.mat index f3e9de3..ae05a93 100644 Binary files a/QD0.mat and b/QD0.mat differ diff --git a/Run_YALMIP.asv b/Run_YALMIP.asv index 583fbe1..653d04b 100644 --- a/Run_YALMIP.asv +++ b/Run_YALMIP.asv @@ -1,11 +1,12 @@ clc clear -yalmip('clear') +clear +% yalmip('clear') tic [kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB, ... Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,PG,QG,PD,QD,CenterA,PGi,PVQU,PVQL, ... Liner,Linex,Lineb,Transforr,Transforx,Transfork0]= ... - pf('E:/算例/东际911_2751267_2012-09-05/newFIle20.txt'); + pf('E:/算例/新民Ⅰ906_2729823_2012-09-06/newFIle20.txt'); %% 潮流等式 AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); @@ -18,7 +19,6 @@ PD0=PD; QD0=QD; PDReal=PD;%真值 QDReal=QD;%真值 -%PD0(12)=PD0(12)+0.001; PG0(Balance)=PGBal(Balance); QG0(Balance)=QGBal(Balance); QG0(PVi)=QGBal(PVi); @@ -36,19 +36,19 @@ xVolt=Volt; xUAngel=UAngel; % VMatrix=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); % dP=PG-PD-diag(xVolt)*(Y.*cos(VMatrix))*xVolt'; -[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,wVolt,PD,PD0,QD,randPDind,Loadi,notLoadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD); -%% 定义变量 +rPD=PD; +rQD=QD; +rVolt=Volt'; +rVAngel=UAngel'; BalVolt=Volt(Balance); -% Volt=sdpvar(Busnum,1); -% UAngel=sdpvar(Busnum,1); -% PG=sdpvar(Busnum,1); -% QG=sdpvar(Busnum,1); -% PD=sdpvar(Busnum,1); -% QD=sdpvar(Busnum,1); -AngleIJ=sdpvar(Busnum,Busnum,'full'); +[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,wVolt,wLoadCurrent,PD,PD0,QD,randPDind,Loadi,notLoadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD); + %% 加噪音 -load('PD0'); -load('QD0'); +load('20PD0'); +load('20QD0'); +mPD=PD0; +mQD=QD0; +mVolt=(1+normrnd(0,0.05,length(rVolt),1)).*rVolt; % PD0=(1+normrnd(0,0.05,length(PD0),1)).*PD0; % save('PD0','PD0'); % QD0=(1+normrnd(0,0.05,length(QD0),1)).*QD0; @@ -61,32 +61,46 @@ PDi=find(PD~=0); QDi=find(QD~=0); % PD0=PD0(PDi); % QD0=QD0(QDi); +%% 电流真实值 +rLoadCurrent=LoadCurrent( rVolt,rVAngel,rPD(PDi),rQD(QDi),PDi,QDi ); +%% 电流测量值 +mLoadCurrent=(1+normrnd(0,0.05,length(rLoadCurrent),1)).*rLoadCurrent; seOpti=Opti(); -seOpti=seOpti.init(Volt0,PDi,QDi,wPD,wQD,PD0,QD0,Y,Angle,r,c,PG,QG,Balance); -opts = optiset('solver','NLOPT'); +seOpti=seOpti.init(mVolt,PDi,QDi,wPD,wQD,wVolt,wLoadCurrent,mPD,mQD,rPD(PDi),rQD(QDi),Y,Angle,r,c,PG,QG,Balance,mLoadCurrent); +opts = optiset('solver','ipopt'); opts.maxiter=85500; -opts.maxtime=3000; +opts.maxtime=30000; opts.maxfeval=85000; opts.maxnodes=85000; opts.tolrfun=1e-4; opts.tolafun=1e-4; opts.warnings='all'; opts.display='iter'; -x0=[zeros(length(PDi)+length(QDi),1); ... +x0=[0.8*rPD(PDi);0.8*rQD(QDi); ... ones(length(Volt),1); ... zeros(length(Volt),1)]; % x0=[PD(PDi);QD(QDi);xVolt';xUAngel']; [~,seOpti]=seOpti.equ(x0); -nlrhs=seOpti.Geteb(); -nle=seOpti.Getgle(); -Opt = opti('fun',@seOpti.obj,'ndec',length(Volt)*2+length(PDi)+length(QDi),'nlmix',@seOpti.equ,nlrhs,nle,'options',opts) -% Opt = opti('fun',@seOpti.obj,'ndec',length(PDi)+length(QDi),'options',opts) +cl=seOpti.Getcl(); +cu=seOpti.Getcu(); +Opt = opti('fun',@seOpti.obj,'ndec',length(Volt)*2+length(PDi)+length(QDi),'nl',@seOpti.equ,cl,cu,'options',opts) +% Opt = opti('fun',@seOpti.obj,'ndec',length(Volt)*2+length(PDi)+length(QDi),'options',opts) [x,fval,exitflag,info] = solve(Opt,x0); info fval=seOpti.obj(x); -fprintf('目标函数: %f\n',fval); +fprintf('目标函数: %.20f\n',fval); toc -rVolt=xVolt0'; +rVolt=Volt0'; rVAngel=xUAngel'; -SEVolt=x(length(PDi)+length(QDi)+1:) -MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel) +SEVolt=x(length(PDi)+length(QDi)+1:length(Volt)+length(PDi)+length(QDi)); +SEVAngel=x(length(Volt)+length(PDi)+length(QDi)+1:end); +fprintf('最大偏差\n') +PD=x(1:length(PDi)); +QD=x(length(PDi)+1:length(PDi)+length(QDi)); +MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel,rPD(PDi),rQD(QDi),PD,QD) +fprintf('统计偏差\n') +StatDeviation(rVolt,SEVolt,rVAngel,SEVAngel,rPD(PDi),rQD(QDi),PD,QD) +%% 约束检查 +% seOpti.equ(x); +% sum([SEVolt;PD;QD]>cu(length(SEVolt)*2+2:end)); +% sum([SEVolt;PD;QD]cu(length(SEVolt)*2+2:end)); +% sum([SEVolt;PD;QD]