From 03865e7b640b17c185b305dbb8bfe94c8b21aeb8 Mon Sep 17 00:00:00 2001 From: facat Date: Mon, 22 Apr 2013 21:56:16 +0800 Subject: [PATCH] =?UTF-8?q?=E5=8A=A0=E5=85=A5=E5=8F=98=E5=8E=8B=E5=99=A8?= =?UTF-8?q?=E6=94=AF=E8=B7=AF?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- @SEOpti/SEOpti.m | 5 ++++ @SEOpti/equ.m | 2 +- @SEOpti/fun.asv | 28 +++++++++++++++++++ @SEOpti/fun.m | 22 ++++++++++----- @SEOpti/init.m | 6 +++- TransPower.m | 12 ++++++++ TransReactivePower.m | 12 ++++++++ run.asv | 65 ++++++++++++++++---------------------------- run.m | 59 ++++++++++++---------------------------- 9 files changed, 119 insertions(+), 92 deletions(-) create mode 100644 @SEOpti/fun.asv create mode 100644 TransPower.m create mode 100644 TransReactivePower.m diff --git a/@SEOpti/SEOpti.m b/@SEOpti/SEOpti.m index e757fa3..d28e6d5 100644 --- a/@SEOpti/SEOpti.m +++ b/@SEOpti/SEOpti.m @@ -20,6 +20,11 @@ classdef SEOpti Balance=NaN; mnle=NaN; mnlrhs=NaN; + mTransP=NaN; + mTransQ=NaN; + transI=NaN; + transJ=NaN; + newwordParameter=NaN; end methods diff --git a/@SEOpti/equ.m b/@SEOpti/equ.m index 8e2bfa0..d1e95d5 100644 --- a/@SEOpti/equ.m +++ b/@SEOpti/equ.m @@ -10,7 +10,7 @@ PQ=diag(cmpSEV)*conj(this.cmpY*cmpSEV); % PQ2=diag(SEVolt)*(this.Y.*sin(YAngle))*SEVolt; out_arg=[real(PQ(this.zerosInjectionIndex));imag(PQ(this.zerosInjectionIndex));]; % out_arg=[PQ1(this.zerosInjectionIndex);PQ2(this.zerosInjectionIndex);]; -out_arg=[out_arg;x(this.Busnum+this.Balance)]; +out_arg=[out_arg;x(this.Balance+this.Busnum)]; out_arg=full(out_arg); this.mnle=zeros(2*length(this.zerosInjectionIndex)+1,1); this.mnlrhs=this.mnle; diff --git a/@SEOpti/fun.asv b/@SEOpti/fun.asv new file mode 100644 index 0000000..f23e6b8 --- /dev/null +++ b/@SEOpti/fun.asv @@ -0,0 +1,28 @@ +function [ Objective ] = fun(this, x ) +%FUN Summary of this function goes here +% Detailed explanation goes here +Objective=0; +SEVolt=x(1:this.Busnum); +SEVAngel=x(this.Busnum+1:2*this.Busnum); +Objective=(SEVolt-this.mVolt)'*(1./this.sigma^2*eye(length(this.mVolt)))*(SEVolt-this.mVolt);%电压 +% Objective=0; +% %% 支路电流 +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);%%电流 +%% 线路功率 +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); +%% 变压器 +newwordParameter=th +transP=TransPower( newwordParameter,cmpY,Volt,VAngel ); +%% 0注入节点 +PQ=diag(cmpSEV)*conj(this.cmpY*cmpSEV); +%% 发电机注入功率 +% Objective=Objective+(this.mPG(this.onlyPG)-real(PQ(this.onlyPG)))'*(1./this.sigma^2*eye(length(this.mPG(this.onlyPG))))*(this.mPG(this.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))); +end + diff --git a/@SEOpti/fun.m b/@SEOpti/fun.m index 6e4e726..47fee30 100644 --- a/@SEOpti/fun.m +++ b/@SEOpti/fun.m @@ -1,6 +1,7 @@ function [ Objective ] = fun(this, x ) %FUN Summary of this function goes here % Detailed explanation goes here +Objective=0; SEVolt=x(1:this.Busnum); SEVAngel=x(this.Busnum+1:2*this.Busnum); Objective=(SEVolt-this.mVolt)'*(1./this.sigma^2*eye(length(this.mVolt)))*(SEVolt-this.mVolt);%电压 @@ -8,17 +9,24 @@ Objective=(SEVolt-this.mVolt)'*(1./this.sigma^2*eye(length(this.mVolt)))*(SEVolt % %% 支路电流 cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压 cmpSEBranchI=(cmpSEV(this.lineI)-cmpSEV(this.lineJ))./(this.lineR+1j*this.lineX);%复数支路电流 -SEBranchI=real(cmpSEBranchI).^2+imag(cmpSEBranchI).^2;% 支路电流幅值 -Objective=Objective+(SEBranchI-this.mBranchI.^2)'*(1./this.sigma^2*eye(length(this.mBranchI)))*(SEBranchI-this.mBranchI.^2);%%电流 -%% 支路功率 +SEBranchI=abs(cmpSEBranchI);% 支路电流幅值 +Objective=Objective+sum((SEBranchI-this.mBranchI.^2).^2./this.sigma^2);%%电流 +%% 线路功率 SEBranchP=real((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(cmpSEBranchI)); SEBranchQ=imag((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(cmpSEBranchI)); -Objective=Objective+(SEBranchP-this.mBranchP)'*(1./this.sigma^2*eye(length(this.mBranchP)))*(SEBranchP-this.mBranchP); -Objective=Objective+(SEBranchQ-this.mBranchQ)'*(1./this.sigma^2*eye(length(this.mBranchQ)))*(SEBranchQ-this.mBranchQ); +Objective=Objective+sum((SEBranchP-this.mBranchP).^2/this.sigma^2); +Objective=Objective+sum((SEBranchQ-this.mBranchQ).^2/this.sigma^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); %% 0注入节点 PQ=diag(cmpSEV)*conj(this.cmpY*cmpSEV); %% 发电机注入功率 -Objective=Objective+(this.mPG(this.onlyPG)-real(PQ(this.onlyPG)))'*(1./this.sigma^2*eye(length(this.mPG(this.onlyPG))))*(this.mPG(this.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+(this.mPG(this.onlyPG)-real(PQ(this.onlyPG)))'*(1./this.sigma^2*eye(length(this.mPG(this.onlyPG))))*(this.mPG(this.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))); end diff --git a/@SEOpti/init.m b/@SEOpti/init.m index 7f5b04e..6b7f752 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) +function [ this ] = init(this,Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ,mTransP,mTransQ) %INIT Summary of this function goes here % Detailed explanation goes here this.mVolt=mVolt; @@ -19,6 +19,10 @@ function [ this ] = init(this,Busnum, mVolt,sigma,newwordParameter,zerosInjectio this.mBranchI=mBranchI; this.mBranchP=mBranchP; this.mBranchQ=mBranchQ; + this.mTransP=mTransP; + this.mTransQ=mTransQ; + this.newwordParameter=newwordParameter; + % this.Y=Y; % this.Yangle=Yangle; % this.r=r; diff --git a/TransPower.m b/TransPower.m new file mode 100644 index 0000000..41fd53c --- /dev/null +++ b/TransPower.m @@ -0,0 +1,12 @@ +function [ output_args ] = TransPower( newwordParameter,cmpY,Volt,VAngel ) +%TRANSPOWER Summary of this function goes here +% Detailed explanation goes here +transI=newwordParameter.trans.transI; +transJ=newwordParameter.trans.transJ; +cmpSEV=Volt.*exp(1j*VAngel); %复数电压 +% cmpY=sparse(r,c,Y.*(1j*YAngel),length(Volt),length(Volt)); +cmpSETransI=(cmpSEV(transI)-cmpSEV(transJ)).*cmpY(sub2ind(size(cmpY),transI,transJ)); +SETransP=real((cmpSEV(transI)-cmpSEV(transJ)).*conj(cmpSETransI)); +output_args=SETransP; +end + diff --git a/TransReactivePower.m b/TransReactivePower.m new file mode 100644 index 0000000..7ca1d2f --- /dev/null +++ b/TransReactivePower.m @@ -0,0 +1,12 @@ +function [ output_args ] = TransReactivePower( newwordParameter,cmpY,Volt,VAngel ) +%TRANSREAVTIVEPOWER Summary of this function goes here +% Detailed explanation goes here +transI=newwordParameter.trans.transI; +transJ=newwordParameter.trans.transJ; +cmpSEV=Volt.*exp(1j*VAngel); %复数电压 +% cmpY=sparse(r,c,Y.*(1j*YAngel),length(Volt),length(Volt)); +cmpSETransI=(cmpSEV(transI)-cmpSEV(transJ)).*cmpY(sub2ind(size(cmpY),transI,transJ)); +SETransQ=imag((cmpSEV(transI)-cmpSEV(transJ)).*conj(cmpSETransI)); +output_args=SETransQ; +end + diff --git a/run.asv b/run.asv index 0c69dff..982b7a1 100644 --- a/run.asv +++ b/run.asv @@ -1,8 +1,8 @@ clear clc -yalmip('clear') +% yalmip('clear') addpath('.\Powerflow') -[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee30.dat', '0'); +[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee14.dat', '0'); %% 量测量 % 电压 节点电流 支路电流 节点功率 支路功率 %% @@ -10,7 +10,7 @@ addpath('.\Powerflow') % 电压 相角 %% %% 开始生成量测量 -sigma=0.05;% 标准差 +sigma=1;% 标准差 %% 电压 %电压幅值 rVolt=Volt; %幅值 @@ -30,6 +30,7 @@ lineJ=newwordParameter.line.lineJ; lineR=newwordParameter.line.lineR; lineX=newwordParameter.line.lineX; lineB2=newwordParameter.line; +transI= cmpBranchI=(cmpV(lineI)-cmpV(lineJ))./(lineR+1j*lineX);%复数支路电流 rBranchI=abs(cmpBranchI);% 支路电流幅值 mBranchI=rBranchI.*(normrnd(0,sigma,length(rBranchI),1)+1);%支路电流量测量 @@ -45,7 +46,7 @@ rQD=QD; QDi=find(QD~=0); rPG=PG; PGi=find(PG~=0); -rQG=QG(QG~=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); @@ -53,61 +54,40 @@ mPG=rPG.*(normrnd(0,sigma,length(rPG),1)+1); mQG=rQG.*(normrnd(0,sigma,length(rQG),1)+1); %% 冗余度计算 stateVarCount=2*length(Volt); -measurements=length(mVolt)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG); +measurements=length(mVolt)+length(mBranchI)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG); fprintf('冗余度 %f\n',measurements/stateVarCount); -%% 进入状态估计计算 -% SEVolt=sdpvar(length(Volt),1,'full'); -% SEVAngel=sdpvar(length(Vangle),1,'full'); -% Objective=(SEVolt-mVolt)'*(1./sigma^2*eye(length(mVolt)))*(SEVolt-mVolt);%电压 -%% 支路电流 -% cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压 -% cmpSEBranchI=(cmpSEV(lineI)-cmpSEV(lineJ))./(lineR+1j*lineX);%复数支路电流 -% SEBranchI=abs(cmpSEBranchI);% 支路电流幅值 -% Objective=Objective+(SEBranchI-mBranchI)'*(1./sigma^2*eye(length(mBranchI)))*(SEBranchI-mBranchI);%%电流 -%% 支路功率 -% SEBranchP=real((cmpSEV(lineI)-cmpSEV(lineJ)).*conj(SEBranchI)); -% SEBranchQ=imag((cmpSEV(lineI)-cmpSEV(lineJ)).*conj(SEBranchI)); -% Objective=Objective+(SEBranchP-mBranchP)'*(1./sigma^2*eye(length(mBranchP)))*(SEBranchP-mBranchP); -% Objective=Objective+(SEBranchQ-mBranchQ)'*(1./sigma^2*eye(length(mBranchQ)))*(SEBranchQ-mBranchQ); %% 0注入节点 zerosInjectionIndex=1:length(Volt); zerosInjectionIndex=zerosInjectionIndex( ~(PD~=0|QD~=0|PG~=0|QG~=0) ); -% PQ=diag(cmpSEV)*cmpY*cmpSEV; -% zeroInjP=real(PQ(zerosInjectionIndex));%% 0注入节点 -% zeroInjQ=imag(PQ(zerosInjectionIndex));%% 0注入节点 +%% 变压器功率 +rTransP=TransPower( newwordParameter,cmpY,rVolt,rVAngel ); +rTransQ=TransReactivePower( newwordParameter,cmpY,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); -% Objective=Objective+(PG(onlyPG)-real(PQ(onlyPG)))'*(1./sigma^2*eye(length(PG(onlyPG))))*(PG(onlyPG)-real(PQ(onlyPG))); -% Objective=Objective+(QG(onlyQG)-imag(PQ(onlyQG)))'*(1./sigma^2*eye(length(QG(onlyQG))))*(QG(onlyQG)-imag(PQ(onlyQG))); -%% YALMIP求解 -% assign(SEVolt,1); -% assign(SEVAngel,1); -% Constraints=[SEVAngel(Balance)==0,zeroInjP==0,zeroInjP==0]; -% Constraints=[SEVAngel(Balance)==0,PQ(zerosInjectionIndex)==0]; -% options = sdpsettings('verbose',2,'solver','filtersd','usex0',1); -% solvesdp([],Objective,options) -% double(Objective) -% double(SEVolt) -% double(SEVAngel) %% Opti ToolBox Busnum=length(Volt); seOpti=SEOpti(); -seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ,Y,Yangle,r,c ); +seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ,mTransP,mTransQ); opts = optiset('solver','ipopt'); -opts.maxiter=2500; -opts.tolrfun=1e-3; -opts.tolafun=1e-3; -opts.warnings='warnings'; -opts.display='iter'; +opts.maxiter=85500; +opts.maxtime=3000; +opts.maxfeval=85000; +opts.maxnodes=85000; +opts.tolrfun=1e-4; +opts.tolafun=1e-4; +opts.warnings='all'; +opts.display='off'; x0=[ones(Busnum,1);zeros(Busnum,1)]; [~,seOpti]=seOpti.equ(x0); nlrhs=seOpti.nlrhs(); nle=seOpti.nle(); -Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'nlmix',@seOpti.equ,nlrhs,nle,'options',opts) -% Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'options',opts) +% Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'nlmix',@seOpti.equ,nlrhs,nle,'options',opts) +Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'options',opts) [x,fval,exitflag,info] = solve(Opt,x0); info %% 输出结果 @@ -117,3 +97,4 @@ fprintf('目 fprintf('相对误差\n'); (abs(rVolt-double(SEVolt)))./(rVolt) MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel) +plotAndComparison( rVolt,rVAngel,SEVolt,SEVAngel ) diff --git a/run.m b/run.m index d7cfedb..791e6a0 100644 --- a/run.m +++ b/run.m @@ -1,8 +1,8 @@ clear clc -yalmip('clear') +% yalmip('clear') addpath('.\Powerflow') -[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee30.dat', '0'); +[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee14.dat', '0'); %% 量测量 % 电压 节点电流 支路电流 节点功率 支路功率 %% @@ -10,7 +10,7 @@ addpath('.\Powerflow') % 电压 相角 %% %% 开始生成量测量 -sigma=0.05;% 标准差 +sigma=1;% 标准差 %% 电压 %电压幅值 rVolt=Volt; %幅值 @@ -30,6 +30,7 @@ lineJ=newwordParameter.line.lineJ; lineR=newwordParameter.line.lineR; lineX=newwordParameter.line.lineX; lineB2=newwordParameter.line; + cmpBranchI=(cmpV(lineI)-cmpV(lineJ))./(lineR+1j*lineX);%复数支路电流 rBranchI=abs(cmpBranchI);% 支路电流幅值 mBranchI=rBranchI.*(normrnd(0,sigma,length(rBranchI),1)+1);%支路电流量测量 @@ -55,62 +56,38 @@ mQG=rQG.*(normrnd(0,sigma,length(rQG),1)+1); stateVarCount=2*length(Volt); measurements=length(mVolt)+length(mBranchI)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG); fprintf('冗余度 %f\n',measurements/stateVarCount); -%% 进入状态估计计算 -% SEVolt=sdpvar(length(Volt),1,'full'); -% SEVAngel=sdpvar(length(Vangle),1,'full'); -% Objective=(SEVolt-mVolt)'*(1./sigma^2*eye(length(mVolt)))*(SEVolt-mVolt);%电压 -%% 支路电流 -% cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压 -% cmpSEBranchI=(cmpSEV(lineI)-cmpSEV(lineJ))./(lineR+1j*lineX);%复数支路电流 -% SEBranchI=abs(cmpSEBranchI);% 支路电流幅值 -% Objective=Objective+(SEBranchI-mBranchI)'*(1./sigma^2*eye(length(mBranchI)))*(SEBranchI-mBranchI);%%电流 -%% 支路功率 -% SEBranchP=real((cmpSEV(lineI)-cmpSEV(lineJ)).*conj(SEBranchI)); -% SEBranchQ=imag((cmpSEV(lineI)-cmpSEV(lineJ)).*conj(SEBranchI)); -% Objective=Objective+(SEBranchP-mBranchP)'*(1./sigma^2*eye(length(mBranchP)))*(SEBranchP-mBranchP); -% Objective=Objective+(SEBranchQ-mBranchQ)'*(1./sigma^2*eye(length(mBranchQ)))*(SEBranchQ-mBranchQ); %% 0注入节点 zerosInjectionIndex=1:length(Volt); zerosInjectionIndex=zerosInjectionIndex( ~(PD~=0|QD~=0|PG~=0|QG~=0) ); -% PQ=diag(cmpSEV)*cmpY*cmpSEV; -% zeroInjP=real(PQ(zerosInjectionIndex));%% 0注入节点 -% zeroInjQ=imag(PQ(zerosInjectionIndex));%% 0注入节点 +%% 变压器功率 +rTransP=TransPower( newwordParameter,cmpY,rVolt,rVAngel ); +rTransQ=TransReactivePower( newwordParameter,cmpY,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); -% Objective=Objective+(PG(onlyPG)-real(PQ(onlyPG)))'*(1./sigma^2*eye(length(PG(onlyPG))))*(PG(onlyPG)-real(PQ(onlyPG))); -% Objective=Objective+(QG(onlyQG)-imag(PQ(onlyQG)))'*(1./sigma^2*eye(length(QG(onlyQG))))*(QG(onlyQG)-imag(PQ(onlyQG))); -%% YALMIP求解 -% assign(SEVolt,1); -% assign(SEVAngel,1); -% Constraints=[SEVAngel(Balance)==0,zeroInjP==0,zeroInjP==0]; -% Constraints=[SEVAngel(Balance)==0,PQ(zerosInjectionIndex)==0]; -% options = sdpsettings('verbose',2,'solver','filtersd','usex0',1); -% solvesdp([],Objective,options) -% double(Objective) -% double(SEVolt) -% double(SEVAngel) %% Opti ToolBox Busnum=length(Volt); seOpti=SEOpti(); -seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ); -opts = optiset('solver','filtersd'); -opts.maxiter=15500; +seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ,mTransP,mTransQ); +opts = optiset('solver','ipopt'); +opts.maxiter=85500; opts.maxtime=3000; -opts.maxfeval=15000; -opts.maxnodes=15000; -opts.tolrfun=1e-5; -opts.tolafun=1e-5; +opts.maxfeval=85000; +opts.maxnodes=85000; +opts.tolrfun=1e-4; +opts.tolafun=1e-4; opts.warnings='all'; opts.display='off'; x0=[ones(Busnum,1);zeros(Busnum,1)]; [~,seOpti]=seOpti.equ(x0); nlrhs=seOpti.nlrhs(); nle=seOpti.nle(); -Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'nlmix',@seOpti.equ,nlrhs,nle,'options',opts) -% Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'options',opts) +% Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'nlmix',@seOpti.equ,nlrhs,nle,'options',opts) +Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'options',opts) [x,fval,exitflag,info] = solve(Opt,x0); info %% 输出结果