diff --git a/@SEOpti/SEOpti.m b/@SEOpti/SEOpti.m index d28e6d5..1963266 100644 --- a/@SEOpti/SEOpti.m +++ b/@SEOpti/SEOpti.m @@ -25,6 +25,10 @@ classdef SEOpti transI=NaN; transJ=NaN; newwordParameter=NaN; + BalanceVolt=NaN; + %真实值 + rVolt=NaN; + rVAngel end methods diff --git a/@SEOpti/equ.m b/@SEOpti/equ.m index 3e286bb..3b5d007 100644 --- a/@SEOpti/equ.m +++ b/@SEOpti/equ.m @@ -11,7 +11,9 @@ PQ=diag(cmpSEV)*conj(this.cmpY*cmpSEV); out_arg=[real(PQ(this.zerosInjectionIndex));imag(PQ(this.zerosInjectionIndex));]; % out_arg=[PQ1(this.zerosInjectionIndex);PQ2(this.zerosInjectionIndex);]; % out_arg=[out_arg;x(this.Balance+this.Busnum)]; -out_arg=[out_arg;x(this.Balance+this.Busnum)]; +out_arg=[out_arg;x(this.Balance+this.Busnum)];%相角等于0 +BalanceVolt=this.BalanceVolt; +out_arg=[out_arg;SEVolt(this.Balance)-BalanceVolt];%平衡节点电压 out_arg=full(out_arg); this.mnle=zeros(length(out_arg),1); this.mnlrhs=this.mnle; diff --git a/@SEOpti/fun.asv b/@SEOpti/fun.asv index e3682ed..c41e226 100644 --- a/@SEOpti/fun.asv +++ b/@SEOpti/fun.asv @@ -4,28 +4,43 @@ function [ Objective ] = fun(this, x ) Objective=0; SEVolt=x(1:this.Busnum); SEVAngel=x(this.Busnum+1:2*this.Busnum); -Objective=sum(((SEVolt-this.mVolt)./this.sigma).^2);%电压 +rVolt=this.rVolt; +sigma=this.sigma; +VoltSigma=rVolt*sigma; +Objective=sum(((SEVolt-this.mVolt)./VoltSigma).^2);%电压 %% 支路电流 cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压 cmpSEBranchI=(cmpSEV(this.lineI)-cmpSEV(this.lineJ))./(this.lineR+1j*this.lineX);%复数支路电流 SEBranchI=abs(cmpSEBranchI);% 支路电流幅值 -Objective=Objective+sum((SEBranchI-this.mBranchI.^2).^2./this.sigma^2);%%电流 +% (SEBranchI-this.mBranchI).^2./(this.sigma^2) +% Objective=Objective+sum((SEBranchI-this.mBranchI).^2./(this.sigma^2));%%电流 %% 线路功率 +rVAngel=this.rVAngel; +cmpV=rVolt.*exp(1j*rVAngel); +rBranchPQ=(cmpV(this.lineI)-cmpV(this.lineJ)).*conj( (cmpV(this.lineI)-cmpV(this.lineJ))./(this.lineR+1j*this.lineX) ); +SEBranchPSigma=real(rBranchPQ)*sigma; +indBranchPS=find(SEBranchPSigma>1e-5); +SEBranchPSigma=SEBranchPSigma(indBranchPS); +SEBranchQSigma=imag(rBranchPQ)*sigma; +indBranchPS=find(SEBranchPSigma>1e-5); +SEBranchQSigma=SEBranchQSigma(abs(SEBranchQSigma)>1e-5); SEBranchP=real((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(cmpSEBranchI)); SEBranchQ=imag((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(cmpSEBranchI)); -Objective=Objective+sum((SEBranchP-this.mBranchP).^2/this.sigma^2); -Objective=Objective+sum((SEBranchQ-this.mBranchQ).^2/this.sigma^2); +Objective=Objective+sum((SEBranchP-this.mBranchP).^2./(SEBranchPSigma.^2)); +% Objective=Objective+sum((SEBranchQ-this.mBranchQ).^2./(SEBranchQSigma.^2)); + + %% 变压器 newwordParameter=this.newwordParameter; cmpY=this.cmpY; transP=TransPower( newwordParameter,cmpY,SEVolt,SEVAngel ); -Objective=Objective+sum((transP-this.mTransP).^2/this.sigma^2); +Objective=Objective+sum((transP-this.mTransP).^2/(this.sigma^2)); transQ=TransReactivePower( newwordParameter,cmpY,SEVolt,SEVAngel ); -Objective=Objective+sum((transQ-this.mTransQ).^2/this.sigma^2); +Objective=Objective+sum((transQ-this.mTransQ).^2/(this.sigma^2)); %% 0注入节点 PQ=diag(cmpSEV)*conj(this.cmpY*cmpSEV); %% 发电机注入功率 -Objective=Objective+(this.mPG(this.onlyPG)-real(PQ(this.onlyPG)))/this.sigmas.onlyPG)-real(PQ(this.onlyPG))); -% Objective=Objective+(this.mQG(this.onlyQG)-imag(PQ(this.onlyQG)))'*(1./this.sigma^2*eye(length(this.mQG(this.onlyQG))))*(this.mQG(this.onlyQG)-imag(PQ(this.onlyQG))); +% Objective=Objective+sum(((this.mPG(this.onlyPG)-real(PQ(this.onlyPG)))/this.sigma).^2); +% Objective=Objective+sum(((this.mQG(this.onlyQG)-imag(PQ(this.onlyQG)))./this.sigma).^2); end diff --git a/@SEOpti/fun.m b/@SEOpti/fun.m index ebe5058..9422f88 100644 --- a/@SEOpti/fun.m +++ b/@SEOpti/fun.m @@ -4,26 +4,39 @@ function [ Objective ] = fun(this, x ) Objective=0; SEVolt=x(1:this.Busnum); SEVAngel=x(this.Busnum+1:2*this.Busnum); -Objective=sum(((SEVolt-this.mVolt)./this.sigma).^2);%电压 +rVolt=this.rVolt; +sigma=this.sigma; +VoltSigma=rVolt*sigma; +Objective=sum(((SEVolt-this.mVolt)./VoltSigma).^2);%电压 %% 支路电流 cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压 cmpSEBranchI=(cmpSEV(this.lineI)-cmpSEV(this.lineJ))./(this.lineR+1j*this.lineX);%复数支路电流 SEBranchI=abs(cmpSEBranchI);% 支路电流幅值 % (SEBranchI-this.mBranchI).^2./(this.sigma^2) -Objective=Objective+sum((SEBranchI-this.mBranchI).^2./(this.sigma^2));%%电流 +% Objective=Objective+sum((SEBranchI-this.mBranchI).^2./(this.sigma^2));%%电流 %% 线路功率 +rVAngel=this.rVAngel; +cmpV=rVolt.*exp(1j*rVAngel); +rBranchPQ=(cmpV(this.lineI)-cmpV(this.lineJ)).*conj( (cmpV(this.lineI)-cmpV(this.lineJ))./(this.lineR+1j*this.lineX) ); +SEBranchPSigma=real(rBranchPQ)*sigma; +indBranchPS=find(real(rBranchPQ)>1e-5);%不考虑太小的数据,以免权重太大。例如除以0 +SEBranchPSigma=SEBranchPSigma(indBranchPS); +SEBranchQSigma=imag(rBranchPQ)*sigma; +indBranchQS=find(imag(rBranchPQ)>1e-5);%不考虑太小的数据,以免权重太大。例如除以0 +SEBranchQSigma=SEBranchQSigma(indBranchQS); SEBranchP=real((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(cmpSEBranchI)); +SEBranchP=SEBranchP(indBranchPS); SEBranchQ=imag((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(cmpSEBranchI)); -% (SEBranchP-this.mBranchP).^2/(this.sigma^2)./(this.mBranchP.^2); -Objective=Objective+sum((SEBranchP-this.mBranchP).^2/(this.sigma^2)); -Objective=Objective+sum((SEBranchQ-this.mBranchQ).^2/(this.sigma^2)); +SEBranchQ=SEBranchQ(indBranchQS); +Objective=Objective+sum((SEBranchP-this.mBranchP).^2./(SEBranchPSigma.^2)); +Objective=Objective+sum((SEBranchQ-this.mBranchQ).^2./(SEBranchQSigma.^2)); %% 变压器 newwordParameter=this.newwordParameter; cmpY=this.cmpY; transP=TransPower( newwordParameter,cmpY,SEVolt,SEVAngel ); Objective=Objective+sum((transP-this.mTransP).^2/(this.sigma^2)); transQ=TransReactivePower( newwordParameter,cmpY,SEVolt,SEVAngel ); -Objective=Objective+sum((transQ-this.mTransQ).^2/(this.sigma^2)); +% Objective=Objective+sum((transQ-this.mTransQ).^2/(this.sigma^2)); %% 0注入节点 PQ=diag(cmpSEV)*conj(this.cmpY*cmpSEV); %% 发电机注入功率 diff --git a/@SEOpti/init.m b/@SEOpti/init.m index 6b7f752..87a38d1 100644 --- a/@SEOpti/init.m +++ b/@SEOpti/init.m @@ -1,4 +1,4 @@ -function [ this ] = init(this,Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ,mTransP,mTransQ) +function [ this ] = init(this,Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ,mTransP,mTransQ,BalanceVolt,rVolt,rVAngel) %INIT Summary of this function goes here % Detailed explanation goes here this.mVolt=mVolt; @@ -22,7 +22,10 @@ function [ this ] = init(this,Busnum, mVolt,sigma,newwordParameter,zerosInjectio this.mTransP=mTransP; this.mTransQ=mTransQ; this.newwordParameter=newwordParameter; - + this.BalanceVolt=BalanceVolt; + %真实值 + this.rVolt=rVolt; + this.rVAngel=rVAngel; % this.Y=Y; % this.Yangle=Yangle; % this.r=r; diff --git a/mBranchI.mat b/mBranchI.mat index 502780e..561bacf 100644 Binary files a/mBranchI.mat and b/mBranchI.mat differ diff --git a/mBranchP.mat b/mBranchP.mat index 1ba197c..1cd943b 100644 Binary files a/mBranchP.mat and b/mBranchP.mat differ diff --git a/mBranchQ.mat b/mBranchQ.mat index d3655ba..53d6409 100644 Binary files a/mBranchQ.mat and b/mBranchQ.mat differ diff --git a/mPG.mat b/mPG.mat index c08e25a..020a2af 100644 Binary files a/mPG.mat and b/mPG.mat differ diff --git a/mQG.mat b/mQG.mat index 6df8a65..f1f11f6 100644 Binary files a/mQG.mat and b/mQG.mat differ diff --git a/mTransP.mat b/mTransP.mat index 57d0292..3591b68 100644 Binary files a/mTransP.mat and b/mTransP.mat differ diff --git a/mTransQ.mat b/mTransQ.mat index f7a3b7a..3e32805 100644 Binary files a/mTransQ.mat and b/mTransQ.mat differ diff --git a/mVolt.mat b/mVolt.mat index 492b786..63417b7 100644 Binary files a/mVolt.mat and b/mVolt.mat differ diff --git a/run.m b/run.m index 830d046..cac1b6b 100644 --- a/run.m +++ b/run.m @@ -10,10 +10,11 @@ addpath('.\Powerflow') % 电压 相角 %% %% 开始生成量测量 -sigma=0.03;% 标准差 +sigma=0.05;% 标准差 %% 电压 %电压幅值 rVolt=Volt; %幅值 +BalanceVolt=Volt(Balance); mVolt=rVolt.*(normrnd(0,sigma,length(Volt),1)+1);%电压量测量 rVAngel=Vangle; %% 电流 @@ -36,8 +37,10 @@ rBranchI=abs(cmpBranchI);% 支路 mBranchI=rBranchI.*(normrnd(0,sigma,length(rBranchI),1)+1);%支路电流量测量 %% 支路功率 rBranchP=real((cmpV(lineI)-cmpV(lineJ)).*conj(cmpBranchI)); +rBranchP=rBranchP(abs(rBranchP)>1e-5); mBranchP=rBranchP.*(normrnd(0,sigma,length(rBranchP),1)+1);%支路功率量测量 rBranchQ=imag((cmpV(lineI)-cmpV(lineJ)).*conj(cmpBranchI)); +rBranchQ=rBranchQ(abs(rBranchQ)>1e-5); mBranchQ=rBranchQ.*(normrnd(0,sigma,length(rBranchQ),1)+1);%支路功率量测量 %% 注入功率 rPD=PD; @@ -82,15 +85,15 @@ seOpti=SEOpti(); % save('mTransP','mTransP'); % save('mTransQ','mTransQ'); %% load -load('mVolt'); -load('mPG'); -load('mQG'); -load('mBranchI'); -load('mBranchP'); -load('mBranchQ'); -load('mTransP'); -load('mTransQ'); -seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ,mTransP,mTransQ); +% load('mVolt'); +% load('mPG'); +% load('mQG'); +% load('mBranchI'); +% load('mBranchP'); +% load('mBranchQ'); +% load('mTransP'); +% load('mTransQ'); +seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ,mTransP,mTransQ,BalanceVolt,rVolt,rVAngel); opts = optiset('solver','ipopt'); opts.maxiter=85500; opts.maxtime=3000; @@ -100,7 +103,7 @@ opts.tolrfun=1e-4; opts.tolafun=1e-4; opts.warnings='all'; opts.display='off'; -x0=[ones(Busnum,1);-0.5*ones(Busnum,1)]; +x0=[ones(Busnum,1);-0.2*ones(Busnum,1)]; % x0=[rVolt;rVAngel]; [~,seOpti]=seOpti.equ(x0); nlrhs=seOpti.nlrhs();