From 657fbe570b4bdfdf1bf6188d475bff69a520fe63 Mon Sep 17 00:00:00 2001
From: facat <dugg@21cn.com>
Date: Wed, 17 Apr 2013 15:43:43 +0800
Subject: [PATCH] =?UTF-8?q?=E5=8A=A0=E6=9C=80=E5=A4=A7=E5=81=8F=E5=B7=AE?=
 =?UTF-8?q?=E7=BB=9F=E8=AE=A1?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Signed-off-by: facat <dugg@21cn.com>
---
 @Opti/equ.m    |  2 +-
 MaxDeviation.m |  8 +++++++
 Run_YALMIP.asv | 65 +++++++++++++-------------------------------------
 Run_YALMIP.m   |  7 +++++-
 4 files changed, 32 insertions(+), 50 deletions(-)
 create mode 100644 MaxDeviation.m

diff --git a/@Opti/equ.m b/@Opti/equ.m
index ad1ae8c..2e82476 100644
--- a/@Opti/equ.m
+++ b/@Opti/equ.m
@@ -24,7 +24,7 @@ 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=sparse(output_args);
+output_args=full(output_args);
 this.gle=zeros(length(output_args),1);
 this.eb=this.gle;
 end
diff --git a/MaxDeviation.m b/MaxDeviation.m
new file mode 100644
index 0000000..ad6f95a
--- /dev/null
+++ b/MaxDeviation.m
@@ -0,0 +1,8 @@
+function [output_arg]=MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel)
+t1=[rVolt;rVAngel];
+t2=[SEVolt;SEVAngel];
+t3=(t1(t1~=0)-double(t2(t1~=0)))./t1(t1~=0);
+t4=abs(t3);
+output_arg=max(t4);
+end
+
diff --git a/Run_YALMIP.asv b/Run_YALMIP.asv
index 2c903fd..583fbe1 100644
--- a/Run_YALMIP.asv
+++ b/Run_YALMIP.asv
@@ -47,50 +47,14 @@ BalVolt=Volt(Balance);
 % QD=sdpvar(Busnum,1);
 AngleIJ=sdpvar(Busnum,Busnum,'full');
 %% 加噪音
-PD0=(1+normrnd(0,0.05,length(PD0),1)).*PD0;
-save('PD0','PD0')
-QD0=(1+normrnd(0,0.05,length(QD0),1)).*QD0;
+load('PD0');
+load('QD0');
+% PD0=(1+normrnd(0,0.05,length(PD0),1)).*PD0;
+% save('PD0','PD0');
+% QD0=(1+normrnd(0,0.05,length(QD0),1)).*QD0;
+% save('QD0','QD0');
 %% 目标函数
-% Objective=ObjectiveFun(PG,PG0,PGi,QG,QG0,PVi,PD,PD0,QD,QD0,Volt,Volt0,wPG,wQG,wPD,wQD,wVolt,Loadi);
-%AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum);
-%% 赋初值，可以加快求解速度。
-% assign(Volt(:),1);
-% assign(UAngel(:),0);
-% assign(PD(:),PD0(:));
-% assign(QD(:),QD0(:));
-% %% YALMIP部分
-% dP=PG0-PD-diag(Volt)*Y.*cos( sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum) )*Volt;
-% dQ=QG0-QD-diag(Volt)*Y.*sin( sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum) )*Volt;
-% Loadi=PD0~=0 | QD0~=0 |PG0~=0|QG0~=0;
-% Constraints = [%AngleIJ-sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum)==0, ...
-%     dP(setdiff(1:Busnum,Loadi))==0, ...
-%     dQ(setdiff(1:Busnum,Loadi))==0, ...
-% %     dP==0, ...
-% %     dQ==0, ...
-%      PD(PD0==0)==0, ...
-%      QD(QD0==0)==0, ...
-%     0.9*ones(Busnum,1)<=Volt<=1.1*ones(Busnum,1), ...
-%     Volt(Balance)==BalVolt, ...
-%     UAngel(Balance)==0, ...
-%     0.8*PD0<=PD<=1.2*PD0;
-%     0.8*QD0<=QD<=1.2*QD0;
-%     ];
-% options = sdpsettings('verbose',2,'showprogress',1,'debug',0,'solver','ipopt','usex0','1');
-% sol = solvesdp(Constraints,Objective,options);
-% if sol.problem == 0
-%     fprintf('Volt\n');
-%     dvolt=double(Volt)
-%     fprintf('VoltAngle\n');
-%     dVangle=double(UAngel)
-%     fprintf('ojb\n');
-%     optimalObj=double(Objective)
-%     sol
-% else
-%     display('Hmm, something went wrong!');
-%     sol.info
-%     sol.solveroutput
-%     yalmiperror(sol.problem)
-% end
+
 %% Opti Toolbox
 Busnum=length(Volt);
 PDi=find(PD~=0);
@@ -99,7 +63,7 @@ QDi=find(QD~=0);
 % QD0=QD0(QDi);
 seOpti=Opti();
 seOpti=seOpti.init(Volt0,PDi,QDi,wPD,wQD,PD0,QD0,Y,Angle,r,c,PG,QG,Balance);
-opts = optiset('solver','filtersd');
+opts = optiset('solver','NLOPT');
 opts.maxiter=85500;
 opts.maxtime=3000;
 opts.maxfeval=85000;
@@ -108,10 +72,10 @@ opts.tolrfun=1e-4;
 opts.tolafun=1e-4;
 opts.warnings='all';
 opts.display='iter';
-% x0=[zeros(length(PDi)+length(QDi),1); ...
-%     ones(length(Volt),1); ...
-%     zeros(length(Volt),1)];
-x0=[PD(PDi);QD(QDi);xVolt';xUAngel'];
+x0=[zeros(length(PDi)+length(QDi),1); ...
+    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();
@@ -119,5 +83,10 @@ Opt = opti('fun',@seOpti.obj,'ndec',length(Volt)*2+length(PDi)+length(QDi),'nlmi
 % Opt = opti('fun',@seOpti.obj,'ndec',length(PDi)+length(QDi),'options',opts)
 [x,fval,exitflag,info] = solve(Opt,x0);
 info
+fval=seOpti.obj(x);
 fprintf('目标函数: %f\n',fval);
 toc
+rVolt=xVolt0';
+rVAngel=xUAngel';
+SEVolt=x(length(PDi)+length(QDi)+1:)
+MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel)
diff --git a/Run_YALMIP.m b/Run_YALMIP.m
index bc392eb..3533c5c 100644
--- a/Run_YALMIP.m
+++ b/Run_YALMIP.m
@@ -85,5 +85,10 @@ Opt = opti('fun',@seOpti.obj,'ndec',length(Volt)*2+length(PDi)+length(QDi),'nlmi
 info
 fval=seOpti.obj(x);
 fprintf('目标函数: %f\n',fval);
-
 toc
+rVolt=Volt0';
+rVAngel=xUAngel';
+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')
+MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel)