把添加噪音后的量测值保存下来，以便下次使用。

Signed-off-by: facat <dugg@21cn.com>
2013-03-05 10:43:18 +08:00 · 2013-03-05 10:43:18 +08:00 · 520d749263
parent 4944e38cbd
commit 520d749263
22 changed files with 94 additions and 53 deletions
--- a/0PD0.mat
+++ b/0PD0.mat
--- a/0QD0.mat
+++ b/0QD0.mat
--- a/12PD0.mat
+++ b/12PD0.mat
--- a/12QD0.mat
+++ b/12QD0.mat
--- a/16PD0.mat
+++ b/16PD0.mat
--- a/16QD0.mat
+++ b/16QD0.mat
--- a/20PD0.mat
+++ b/20PD0.mat
--- a/20QD0.mat
+++ b/20QD0.mat
--- a/4PD0.mat
+++ b/4PD0.mat
--- a/4QD0.mat
+++ b/4QD0.mat
--- a/8PD0.mat
+++ b/8PD0.mat
--- a/8QD0.mat
+++ b/8QD0.mat
--- a/@ForThesis/MaxBranchDeviation.m
+++ b/@ForThesis/MaxBranchDeviation.m
@ -1,7 +1,34 @@
-function [ output_args ] = MaxBranchDeviation( ~, Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel)
+function [ output_args ] = MaxBranchDeviation( ~, Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel,FileName,PD0,QD0)
 %% 最大支路功率偏差
 % 支路功率包括线路和变压器
-%% 
+%% 先用加了误差的负荷功率计算潮流值
+[Busnum,Balance,PQstandard,Precision,~,~,~,~,~,kmax,~ ,...
+    ~,~,~,~,~,~,~,Pointpoweri,PG,QG,~,~,PVi,PVu,~,~,~,~,~,~,~]= openfile2(FileName);
+PD=PD0;
+QD=QD0;
+%% 形成节点导纳矩阵
+[~,Y,r,c,Angle] = admmatrix(Busnum,Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,...
+    Transforx,Transfork0,Branchi,Branchg,Branchb);
+[P0,Q0,U,Uangle] = Initial(PG,PD,PQstandard,Pointpoweri,QG,QD,Busnum);                         %求功率不平衡量
+%disp('迭代次数i  最大不平衡量');
+%% 循环体计算
+for i = 0:kmax
+    [Jacob,PQ,U,Uangle] = jacobian(Busnum,Balance,PVi,PVu,U,Uangle,Y,Angle,P0,Q0,r,c);  %形成雅克比矩阵
+   % disp('第一次雅克比');
+    m = max(abs(PQ));
+    m=full(m);
+    %fprintf('       %u         %.8f \n',i,m);
+    if m > Precision                                                                                                                                 %判断不平衡量是否满足精度要求
+        [Uangle,U] = solvefun(Busnum,Jacob,PQ,Uangle,U);                                                           %求解修正方程，更新电压变量
+    else
+        disp(['收敛，迭代次数为',num2str(i),'次']);
+        break                                                                                                                                             %若满足精度要求，则计算收敛
+    end
+end
+
+Volt0=U;
+UAngel0=Uangle;
+
 [dispLineloss0 dispTransloss0]=Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0);
 [dispLineloss dispTransloss]=Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt,UAngel);
 t1=(dispLineloss0(:,3) - dispLineloss(:,3))./dispLineloss0(:,3);
--- a/@ForThesis/MaxDeviation.m
+++ b/@ForThesis/MaxDeviation.m
@ -7,5 +7,5 @@ PG0Array=repmat(PG0,this.sampleNum,1);
 QG0Array=repmat(QG0,this.sampleNum,1);
 PGMaxDev=max(abs( (PG0Array-this.PGArray)./PG0Array  ));
 QGMaxDev=max( abs( (QG0Array-this.QGArray)./QG0Array  ) );
-output_arg=max([PDMaxDev;QDMaxDev;PGMaxDev;QGMaxDev]);
+output_arg=max([PDMaxDev;QDMaxDev;PGMaxDev;QGMaxDev;]);
 end
--- a/@ForThesis/StatBranchDeviation.m
+++ b/@ForThesis/StatBranchDeviation.m
@ -14,6 +14,6 @@ t22=dispTransloss0(:,3)>1e-5;% ̫С
 t3=[t1(t11);t2(t22)];
 t4=t3.^2;
 t5=sum(t4).^.5;
-output_args=t5
+output_args=t5;
 end

--- a/FormLw.m
+++ b/FormLw.m
@ -6,13 +6,13 @@ VoltU=(1.1)*ones(1,Busnum);
 %VoltU=10*ones(1,Busnum);
 PDU=PD0(Loadi);
 % PDU=noDataTransCapacity;
-PDU(PDU>0)=1.300*PDU(PDU>0);
-PDU(PDU<0)=0.700*PDU(PDU<0);
+PDU(PDU>0)=1.200*PDU(PDU>0);
+PDU(PDU<0)=0.800*PDU(PDU<0);
 PDU(PDU==0)=0.400;
 %PDU=10*ones(length(Loadi),1);
 QDU=QD0(Loadi);
-QDU(QDU>0)=1.300*QDU(QDU>0);
-QDU(QDU<0)=0.700*QDU(QDU<0);
+QDU(QDU>0)=1.200*QDU(QDU>0);
+QDU(QDU<0)=0.800*QDU(QDU<0);
 QDU(QDU==0)=0.400;
 % PF=0.85;
 % QDU=1.0*PD(Loadi).*sqrt(1 -PF.^2)./PF;
--- a/FormLz.m
+++ b/FormLz.m
@ -4,13 +4,13 @@ KK=999;
 VoltL=(0.9)*ones(1,Busnum);
 %VoltL=-10*ones(1,Busnum);
 PDL=PD0(Loadi);
-PDL(PDL>0)=0.700*PDL(PDL>0);
-PDL(PDL<0)=1.300*PDL(PDL<0);
+PDL(PDL>0)=0.800*PDL(PDL>0);
+PDL(PDL<0)=1.200*PDL(PDL<0);
 PDL(PDL==0)=-0.400;
 %PDL=-10*ones(length(Loadi),1);
 QDL=QD0(Loadi);
-QDL(QDL>0)=0.700*QDL(QDL>0);
-QDL(QDL<0)=1.300*QDL(QDL<0);
+QDL(QDL>0)=0.800*QDL(QDL>0);
+QDL(QDL<0)=1.200*QDL(QDL<0);
 QDL(QDL==0)=-0.400;
 % QDL=0*PD(Loadi).*sqrt((1-PF.^2))./PF;
 t1=([PDL',QDL',VoltL])';
--- a/Lineloss.asv
+++ b/Lineloss.asv
@ -1,10 +1,10 @@
 %% 计算线损
-function Lineloss(Linei,Linej,Liner,Linex,Lineb2,Transi,Transj,Transr,Transx,k0,Volt,Angle)
+function [dispLineloss dispTransloss]=Lineloss(Linei,Linej,Liner,Linex,Lineb2,Transi,Transj,Transr,Transx,Branchi,Branchg,Branchb,k0,Volt,Angle)
 %format long
-fprintf('功率为有名值\n');
-fprintf('节点号\t节点号\t有功')
-cmpVolt=Volt.*cos(Angle)+1i*Volt.*sin(Angle);
-cmpVolt=cmpVolt';
+% fprintf('功率为有名值\n');
+% fprintf('节点号\t节点号\t有功损耗 MW\t无功损耗 MVar')
+cmpVolt=Volt'.*cos(Angle')+1i*Volt'.*sin(Angle');
+% cmpVolt=cmpVolt';
 y0=1i*Lineb2;
 %yj0=1i*standardinput(:,7);
 yij=1./(Liner+1i*Linex);
@ -13,16 +13,26 @@ Sij=cmpVolt(Linei).*conj(  cmpVolt(Linei) .* y0 + ( cmpVolt(Linei)- cmpVolt(Line
 Sji=cmpVolt(Linej).*conj(  cmpVolt(Linej) .*y0 + ( cmpVolt(Linej)- cmpVolt(Linei) ).*yij );
 %Sij1==Sij2
 deltLineS=Sij+Sji;
+
+%% 另一种计算方式begin
+ss=1*(Volt(Linei)'.^2.*abs(yij).*cos( angle(yij) ) -Volt(Linei)'.*Volt(Linej)'.*cos( Angle(Linei)' - Angle(Linej)' - angle(yij)).*abs(yij));
+ss=(Volt(Linei)'.^2+Volt(Linej)'.^2).*abs(yij).*cos(angle(yij))-2*Volt(Linei)'.*Volt(Linej)'.*cos( Angle(Linei)' - Angle(Linej)').*cos( - angle(yij)).*abs(yij);
+ss=Volt(Linei)'.*Volt(Linej)'.*abs(yij).*cos()
+%% 另一种计算方式end
 dispLineloss=[Linei Linej real(deltLineS)*100 imag(deltLineS)*100];
 %full(dispLineloss)
-dispLineloss=sortrows(dispLineloss,-3);
-full(dispLineloss)
+% dispLineloss=sortrows(dispLineloss,-3);
+full(dispLineloss);
 %% 以下是变压器损耗
 yij=1./(Transr+1i*Transx);
 Sij=cmpVolt(Transi)./k0.*conj(  ( cmpVolt(Transi)./k0- cmpVolt(Transj) ).*yij );
 Sji=cmpVolt(Transj).*conj( ( cmpVolt(Transj)- cmpVolt(Transi)./k0 ).*yij );
 deltTransS=Sij+Sji;
+%% 接地支路损耗
+% 没有考虑变压器变比
+deltTransS =deltTransS+sum(cmpVolt(Branchi).*conj((cmpVolt(Branchi).*(Branchg+1j*Branchb))));
+%% 
 dispTransloss=[Transi Transj real(deltTransS)*100 imag(deltTransS)*100];
-dispTransloss=sortrows(dispTransloss,-3);
-full(dispTransloss)
+% dispTransloss=sortrows(dispTransloss,-3);
+full(dispTransloss);
 end
--- a/Lineloss.m
+++ b/Lineloss.m
@ -3,8 +3,8 @@ function [dispLineloss dispTransloss]=Lineloss(Linei,Linej,Liner,Linex,Lineb2,Tr
 %format long
 % fprintf('功率为有名值\n');
 % fprintf('节点号\t节点号\t有功损耗 MW\t无功损耗 MVar')
-cmpVolt=Volt.*cos(Angle)+1i*Volt.*sin(Angle);
-cmpVolt=cmpVolt';
+cmpVolt=Volt'.*cos(Angle')+1i*Volt'.*sin(Angle');
+% cmpVolt=cmpVolt';
 y0=1i*Lineb2;
 %yj0=1i*standardinput(:,7);
 yij=1./(Liner+1i*Linex);
@ -13,6 +13,12 @@ Sij=cmpVolt(Linei).*conj(  cmpVolt(Linei) .* y0 + ( cmpVolt(Linei)- cmpVolt(Line
 Sji=cmpVolt(Linej).*conj(  cmpVolt(Linej) .*y0 + ( cmpVolt(Linej)- cmpVolt(Linei) ).*yij );
 %Sij1==Sij2
 deltLineS=Sij+Sji;
+
+%% 另一种计算方式begin
+% ss=1*(Volt(Linei)'.^2.*abs(yij).*cos( angle(yij) ) -Volt(Linei)'.*Volt(Linej)'.*cos( Angle(Linei)' - Angle(Linej)' - angle(yij)).*abs(yij));
+% ss=(Volt(Linei)'.^2+Volt(Linej)'.^2).*abs(yij).*cos(angle(yij))-2*Volt(Linei)'.*Volt(Linej)'.*cos( Angle(Linei)' - Angle(Linej)').*cos( - angle(yij)).*abs(yij);
+% ss=abs(yij).*cos(angle(yij)).*(  Volt(Linei)'.^2+Volt(Linej)'.^2-   2*Volt(Linei)'.*Volt(Linej)' .*cos(Angle(Linei)' - Angle(Linej)')  );
+%% 另一种计算方式end
 dispLineloss=[Linei Linej real(deltLineS)*100 imag(deltLineS)*100];
 %full(dispLineloss)
 dispLineloss=sortrows(dispLineloss,-3);
@ -24,7 +30,8 @@ Sji=cmpVolt(Transj).*conj( ( cmpVolt(Transj)- cmpVolt(Transi)./k0 ).*yij );
 deltTransS=Sij+Sji;
 %% 接地支路损耗
 % 没有考虑变压器变比
-deltTransS =deltTransS+sum(cmpVolt(Branchi).*conj((cmpVolt(Branchi).*(Branchg+1j*Branchb))));
+deltTransS =deltTransS+cmpVolt(Branchi).*conj((cmpVolt(Branchi).*(Branchg+1j*Branchb)));
+ss=Volt(Branchi)'.^2.*(Branchg)
 %% 
 dispTransloss=[Transi Transj real(deltTransS)*100 imag(deltTransS)*100];
 dispTransloss=sortrows(dispTransloss,-3);
--- a/OPF.m
+++ b/OPF.m
@ -4,15 +4,13 @@ clear
 %% 存在问题
 % 变压器变比的位置没有考虑，由于现在用的变比都是1，所以没有影响。 20130123
 %% 
-thesis=ForThesis(4,8);
-for II=1:4
+thesis=ForThesis(1,45);
+for II=1:1
    [kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,PG,QG,PD,QD,CenterA,PGi,PVQU,PVQL,Liner,Linex,Lineb,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0]= ...
-        pf('E:\算例\东际911_2751267_2012-09-05\newFIle20.txt');
+        pf('E:\ËãÀý\Áø½ð¢ñ926_21671693_2012-09-06\newFIle20.txt');
    %pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt');
    %pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt');
-    %pf('c:/file31.txt');
-    
-    
+    %pf('c:/file31.txt'); 
    %% 计算功率因数
    Loadi=QD~=0 | PD~=0;
    PF=sqrt(PD(Loadi).^2./(QD(Loadi).^2+PD(Loadi).^2))
@ -53,17 +51,16 @@ for II=1:4
    QD_NON_ZERO=QD(PD==0 & QD~=0);
    QD_NON_ZERO_IND=find(PD==0 & QD~=0);
    %%
-    Precision=Precision/10;
-    
+    Precision=Precision/1;
    %% 加误差
-    PD0(Loadi)=PD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
-    QD0(Loadi)=QD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
+    PD0(Loadi)=PD0(Loadi).*(1+normrnd(0,0.5,length(Loadi),1));
+    QD0(Loadi)=QD0(Loadi).*(1+normrnd(0,0.5,length(Loadi),1));
+    load('0PD0.mat');
+    load('0QD0.mat');
 %     PG0(PGi)=PG0(PGi).*(1+normrnd(0,0.01,length(PGi),1));
 %     QG0(PVi)=QG0(PVi).*(1+normrnd(0,0.01,length(PVi),1));
    %% 读变压器容量
-    %[noDataTransNum noDataTransCapacity noDataTransPowerFactor]=ReadNoDataTrans('C:/b/东际911_2751267_2012-09-05/iPso_东际911_2751267_2012-09-05_变压器无负载.txt');
    noDataTransCapacity=0;
-    
    while(abs(Gap)>Precision)
        if KK>kmax
            break;
@ -87,7 +84,6 @@ for II=1:4
        ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi,PD,QD);
        %%  开始构建deltF
        deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,PG0,QG0,PD0,PD,QD,QD0,Busnum,Loadi);
-        
        %% 形成方程矩阵
        Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
        Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1);
@ -136,12 +132,12 @@ PD(Loadi)=thesis.MeanPD();
 QD(Loadi)=thesis.MeanQD();
 PG(Balance)=thesis.MeanPG();
 QG(PVi)=thesis.MeanQG();
-thesis.MaxDeviation(RealPG(Balance),RealQG(PVi),RealPD(Loadi),RealQD(Loadi))
-thesis.StatDeviation(RealPG(Balance),RealQG(PVi),RealPD(Loadi),RealQD(Loadi))
+thesis.MaxDeviation(PG0(Balance),QG0(PVi),PD0(Loadi),QD0(Loadi))
+thesis.StatDeviation(PG0(Balance),QG0(PVi),PD0(Loadi),QD0(Loadi))
 % thesis.SqareDeviation(RealPG(Balance(1)),RealQG(PVi(1)),RealPD(Loadi),RealQD(Loadi))
 thesis.SqareDeviation(PG0(Balance(1)),QG0(PVi(1)),PD0(Loadi),QD0(Loadi))
 thesis.PercentOfPass(PG0(Balance(1)),QG0(PVi(1)),PD0(Loadi),QD0(Loadi))
-thesis.MaxBranchDeviation(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel)
+thesis.MaxBranchDeviation(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel,'E:\ËãÀý\Áø½ð¢ñ926_21671693_2012-09-06\newFIle0.txt',PD0,QD0)
 thesis.StatBranchDeviation(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel)
 toc

--- a/OPF_Init.asv
+++ b/OPF_Init.asv
@ -1,7 +1,7 @@
-function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wD,PD,PD0,QD,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD)
+function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,PD,PD0,QD,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD)
 Loadi=find(QD~=0 | PD~=0);
 %Loadi=[1:Busnum]';
-RestraintCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)*1+Busnum*1+1; %约束条件数,只放开一个QD
+RestraintCount=size(Loadi,1)*2+Busnum*1; %约束条件数,放开所有QD
 t_Bal_volt=Volt(Balance);
 Volt=sparse(1*ones(1,Busnum));
 Volt(Balance)=t_Bal_volt;
@ -15,19 +15,20 @@ tPU=sparse(GenU(:,2));%
 tQU=sparse(PVQU(:,1));% 无功上限
 tPL=sparse(GenL(:,2));% 发电机有功下限
 tQL=sparse(PVQL(:,1));% 无功下限
-PG(PGi)=(tPU+tPL)/2;
-QG(PVi)=(tQU+tQL)/2;
-wPG=1*ones(size(PGi,1),1);
-wQG=1*zeros(size(PVi,1),1);
+wPG=0;
+wQG=0;
 %randInt=randperm(size(Loadi,1));
 %randPDind=randInt(1:10);
 randPDind=0;
-wD=1*zeros(size(Loadi,1),1);
+wPD=1/.5^2*zeros(size(Loadi,1),1);
+wPD(1:2:end)=0;
+wQD=1/.5^2*ones(size(Loadi,1),1);
+wQD(1:2:end)=0;
 %wD(randPDind)=0;%一些负荷不约束
 %wD(7)=0;
 % wD(11)=0;
-PD=1*PD0;
+PD=0.5*PD0;
 %powerFacter=0.98;
 %QD=PD*sqrt((1-powerFacter^2)/powerFacter^2);
-QD=QD0;
+QD=.5*QD0;
 end
--- a/OPF_Init.m
+++ b/OPF_Init.m
@ -20,10 +20,10 @@ wQG=0;
 %randInt=randperm(size(Loadi,1));
 %randPDind=randInt(1:10);
 randPDind=0;
-wPD=0/.05^2*ones(size(Loadi,1),1);
-% wPD([2 7])=0;
-wQD=0/.05^2*ones(size(Loadi,1),1);
-% wQD([2 7])=0;
+wPD=1/.5^2*ones(size(Loadi,1),1);
+% wPD(1:2:end)=0;
+wQD=1/.5^2*ones(size(Loadi,1),1);
+% wQD(1:2:end)=0;
 %wD(randPDind)=0;%一些负荷不约束
 %wD(7)=0;
 % wD(11)=0;