From f50bbd343594ce5620378322bf9ac43066d87f13 Mon Sep 17 00:00:00 2001
From: facat <dugg@21cn.com>
Date: Sun, 14 Apr 2013 16:30:53 +0800
Subject: [PATCH] =?UTF-8?q?=E5=8A=A0=E5=8F=91=E7=94=B5=E6=9C=BA=E6=B3=A8?=
 =?UTF-8?q?=E5=85=A5=E5=8A=9F=E7=8E=87?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Signed-off-by: facat <dugg@21cn.com>
---
 run.asv | 14 ++++++++++----
 run.m   | 14 ++++++++++----
 2 files changed, 20 insertions(+), 8 deletions(-)

diff --git a/run.asv b/run.asv
index da0a6aa..3e7cdaa 100644
--- a/run.asv
+++ b/run.asv
@@ -65,14 +65,20 @@ 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(SEVolt)*conj(cmpY*SEVolt);
+zeroInjP=real(PQ(zerosInjectionIndex));%% 0注入节点
+zeroInjQ=imag(PQ(zerosInjectionIndex));%% 0注入节点
+%% 发电机注入功率
+% 先找到只有发电机的节点
+PDQDi=union(PDi,QDi);
+onlyPG=setdiff(PGi,PDQDi);
+onlyQG=setdiff(QGi,PDQDi);
+PG-real(PQ())
 %% YALMIP求解
-PQ=diag(SEVolt)*cmpY*SEVolt;
-zeros
-Constraints=[SEVAngel(Balance)==0,];
+Constraints=[SEVAngel(Balance)==0,zeroInjP==0,zeroInjP==0];
 % Constraints=[[zeros(length(c)) A' -eye(length(lbounds))]*x==-c;[A zeros(length(b)) zeros(length(b),length(lbounds))]*x<=b;10>=x>=0];
 options = sdpsettings('verbose',2,'solver','ipopt');
 solvesdp(Constraints,Objective,options)
diff --git a/run.m b/run.m
index 3ee6150..f860291 100644
--- a/run.m
+++ b/run.m
@@ -1,5 +1,5 @@
 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');
 %% 量测量
 % 电压 节点电流 支路电流 节点功率 支路功率
 %%
@@ -65,14 +65,20 @@ 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) );
-%% YALMIP求解
-PQ=diag(SEVolt)*cmpY*SEVolt;
+PQ=diag(SEVolt)*conj(cmpY*SEVolt);
 zeroInjP=real(PQ(zerosInjectionIndex));%% 0注入节点
 zeroInjQ=imag(PQ(zerosInjectionIndex));%% 0注入节点
+%% 发电机注入功率
+% 先找到只有发电机的节点
+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求解
 Constraints=[SEVAngel(Balance)==0,zeroInjP==0,zeroInjP==0];
 % Constraints=[[zeros(length(c)) A' -eye(length(lbounds))]*x==-c;[A zeros(length(b)) zeros(length(b),length(lbounds))]*x<=b;10>=x>=0];
 options = sdpsettings('verbose',2,'solver','ipopt');