facat
/
stateestimation-self-derivatives-measurementdistribution
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

准备自己写代码求导 

Signed-off-by: facat <facat@facat.cn>
This commit is contained in:
facat 2013-08-14 10:46:45 +08:00
parent 202e7e62af
commit 006adb9fc0
14 changed files with 68 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
@SEOpti/SEOpti.m
View File

@ -25,6 +25,10 @@ classdef SEOpti
transI=NaN; transI=NaN;
transJ=NaN; transJ=NaN;
newwordParameter=NaN; newwordParameter=NaN;
BalanceVolt=NaN;
%真实值
rVolt=NaN;
rVAngel
end end
methods methods

4
@SEOpti/equ.m
View File

@ -11,7 +11,9 @@ PQ=diag(cmpSEV)*conj(this.cmpY*cmpSEV);
out_arg=[real(PQ(this.zerosInjectionIndex));imag(PQ(this.zerosInjectionIndex));]; out_arg=[real(PQ(this.zerosInjectionIndex));imag(PQ(this.zerosInjectionIndex));];
% out_arg=[PQ1(this.zerosInjectionIndex);PQ2(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)]; 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); out_arg=full(out_arg);
this.mnle=zeros(length(out_arg),1); this.mnle=zeros(length(out_arg),1);
this.mnlrhs=this.mnle; this.mnlrhs=this.mnle;

31
@SEOpti/fun.asv
View File

@ -4,28 +4,43 @@ function [ Objective ] = fun(this, x )
Objective=0; Objective=0;
SEVolt=x(1:this.Busnum); SEVolt=x(1:this.Busnum);
SEVAngel=x(this.Busnum+1:2*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); %复数电压 cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压
cmpSEBranchI=(cmpSEV(this.lineI)-cmpSEV(this.lineJ))./(this.lineR+1j*this.lineX);%复数支路电流 cmpSEBranchI=(cmpSEV(this.lineI)-cmpSEV(this.lineJ))./(this.lineR+1j*this.lineX);%复数支路电流
SEBranchI=abs(cmpSEBranchI);% 支路电流幅值 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)); SEBranchP=real((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(cmpSEBranchI));
SEBranchQ=imag((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((SEBranchP-this.mBranchP).^2./(SEBranchPSigma.^2));
Objective=Objective+sum((SEBranchQ-this.mBranchQ).^2/this.sigma^2); % Objective=Objective+sum((SEBranchQ-this.mBranchQ).^2./(SEBranchQSigma.^2));
%% 变压器 %% 变压器
newwordParameter=this.newwordParameter; newwordParameter=this.newwordParameter;
cmpY=this.cmpY; cmpY=this.cmpY;
transP=TransPower( newwordParameter,cmpY,SEVolt,SEVAngel ); 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 ); 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注入节点 %% 0注入节点
PQ=diag(cmpSEV)*conj(this.cmpY*cmpSEV); 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+sum(((this.mPG(this.onlyPG)-real(PQ(this.onlyPG)))/this.sigma).^2);
% 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.mQG(this.onlyQG)-imag(PQ(this.onlyQG)))./this.sigma).^2);
end end

25
@SEOpti/fun.m
View File

@ -4,26 +4,39 @@ function [ Objective ] = fun(this, x )
Objective=0; Objective=0;
SEVolt=x(1:this.Busnum); SEVolt=x(1:this.Busnum);
SEVAngel=x(this.Busnum+1:2*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); %复数电压 cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压
cmpSEBranchI=(cmpSEV(this.lineI)-cmpSEV(this.lineJ))./(this.lineR+1j*this.lineX);%复数支路电流 cmpSEBranchI=(cmpSEV(this.lineI)-cmpSEV(this.lineJ))./(this.lineR+1j*this.lineX);%复数支路电流
SEBranchI=abs(cmpSEBranchI);% 支路电流幅值 SEBranchI=abs(cmpSEBranchI);% 支路电流幅值
% (SEBranchI-this.mBranchI).^2./(this.sigma^2) % (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=real((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(cmpSEBranchI));
SEBranchP=SEBranchP(indBranchPS);
SEBranchQ=imag((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(cmpSEBranchI)); SEBranchQ=imag((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(cmpSEBranchI));
% (SEBranchP-this.mBranchP).^2/(this.sigma^2)./(this.mBranchP.^2); SEBranchQ=SEBranchQ(indBranchQS);
Objective=Objective+sum((SEBranchP-this.mBranchP).^2/(this.sigma^2)); Objective=Objective+sum((SEBranchP-this.mBranchP).^2./(SEBranchPSigma.^2));
Objective=Objective+sum((SEBranchQ-this.mBranchQ).^2/(this.sigma^2)); Objective=Objective+sum((SEBranchQ-this.mBranchQ).^2./(SEBranchQSigma.^2));
%% 变压器 %% 变压器
newwordParameter=this.newwordParameter; newwordParameter=this.newwordParameter;
cmpY=this.cmpY; cmpY=this.cmpY;
transP=TransPower( newwordParameter,cmpY,SEVolt,SEVAngel ); 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 ); 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注入节点 %% 0注入节点
PQ=diag(cmpSEV)*conj(this.cmpY*cmpSEV); PQ=diag(cmpSEV)*conj(this.cmpY*cmpSEV);
%% 发电机注入功率 %% 发电机注入功率

7
@SEOpti/init.m
View File

@ -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 %INIT Summary of this function goes here
% Detailed explanation goes here % Detailed explanation goes here
this.mVolt=mVolt; this.mVolt=mVolt;
@ -22,7 +22,10 @@ function [ this ] = init(this,Busnum, mVolt,sigma,newwordParameter,zerosInjectio
this.mTransP=mTransP; this.mTransP=mTransP;
this.mTransQ=mTransQ; this.mTransQ=mTransQ;
this.newwordParameter=newwordParameter; this.newwordParameter=newwordParameter;
this.BalanceVolt=BalanceVolt;
%Õæʵֵ
this.rVolt=rVolt;
this.rVAngel=rVAngel;
% this.Y=Y; % this.Y=Y;
% this.Yangle=Yangle; % this.Yangle=Yangle;
% this.r=r; % this.r=r;

BIN
mBranchI.mat
View File

Binary file not shown.

BIN
mBranchP.mat
View File

Binary file not shown.

BIN
mBranchQ.mat
View File

Binary file not shown.

BIN
mPG.mat
View File

Binary file not shown.

BIN
mQG.mat
View File

Binary file not shown.

BIN
mTransP.mat
View File

Binary file not shown.

BIN
mTransQ.mat
View File

Binary file not shown.

BIN
mVolt.mat
View File

Binary file not shown.

25
run.m
View File

@ -10,10 +10,11 @@ addpath('.\Powerflow')
% 电压 相角 % 电压 相角
%% %%
%% 开始生成量测量 %% 开始生成量测量
sigma=0.03;% ±ê×¼²î sigma=0.05;% ±ê×¼²î
%% 电压 %% 电压
%电压幅值 %电压幅值
rVolt=Volt; %幅值 rVolt=Volt; %幅值
BalanceVolt=Volt(Balance);
mVolt=rVolt.*(normrnd(0,sigma,length(Volt),1)+1);%电压量测量 mVolt=rVolt.*(normrnd(0,sigma,length(Volt),1)+1);%电压量测量
rVAngel=Vangle; rVAngel=Vangle;
%% 电流 %% 电流
@ -36,8 +37,10 @@ rBranchI=abs(cmpBranchI);% ֧·
mBranchI=rBranchI.*(normrnd(0,sigma,length(rBranchI),1)+1);%支路电流量测量 mBranchI=rBranchI.*(normrnd(0,sigma,length(rBranchI),1)+1);%支路电流量测量
%% 支路功率 %% 支路功率
rBranchP=real((cmpV(lineI)-cmpV(lineJ)).*conj(cmpBranchI)); rBranchP=real((cmpV(lineI)-cmpV(lineJ)).*conj(cmpBranchI));
rBranchP=rBranchP(abs(rBranchP)>1e-5);
mBranchP=rBranchP.*(normrnd(0,sigma,length(rBranchP),1)+1);%支路功率量测量 mBranchP=rBranchP.*(normrnd(0,sigma,length(rBranchP),1)+1);%支路功率量测量
rBranchQ=imag((cmpV(lineI)-cmpV(lineJ)).*conj(cmpBranchI)); rBranchQ=imag((cmpV(lineI)-cmpV(lineJ)).*conj(cmpBranchI));
rBranchQ=rBranchQ(abs(rBranchQ)>1e-5);
mBranchQ=rBranchQ.*(normrnd(0,sigma,length(rBranchQ),1)+1);%支路功率量测量 mBranchQ=rBranchQ.*(normrnd(0,sigma,length(rBranchQ),1)+1);%支路功率量测量
%% 注入功率 %% 注入功率
rPD=PD; rPD=PD;
@ -82,15 +85,15 @@ seOpti=SEOpti();
% save('mTransP','mTransP'); % save('mTransP','mTransP');
% save('mTransQ','mTransQ'); % save('mTransQ','mTransQ');
%% load %% load
load('mVolt'); % load('mVolt');
load('mPG'); % load('mPG');
load('mQG'); % load('mQG');
load('mBranchI'); % load('mBranchI');
load('mBranchP'); % load('mBranchP');
load('mBranchQ'); % load('mBranchQ');
load('mTransP'); % load('mTransP');
load('mTransQ'); % load('mTransQ');
seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ,mTransP,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 = optiset('solver','ipopt');
opts.maxiter=85500; opts.maxiter=85500;
opts.maxtime=3000; opts.maxtime=3000;
@ -100,7 +103,7 @@ opts.tolrfun=1e-4;
opts.tolafun=1e-4; opts.tolafun=1e-4;
opts.warnings='all'; opts.warnings='all';
opts.display='off'; opts.display='off';
x0=[ones(Busnum,1);-0.5*ones(Busnum,1)]; x0=[ones(Busnum,1);-0.2*ones(Busnum,1)];
% x0=[rVolt;rVAngel]; % x0=[rVolt;rVAngel];
[~,seOpti]=seOpti.equ(x0); [~,seOpti]=seOpti.equ(x0);
nlrhs=seOpti.nlrhs(); nlrhs=seOpti.nlrhs();