From 382763f559065e83107018a6a58f382336bdf31b Mon Sep 17 00:00:00 2001 From: "dugg@lab-desk" Date: Mon, 21 Jul 2014 21:49:50 +0800 Subject: [PATCH] =?UTF-8?q?1.=E4=B8=BA=E4=BA=86=E5=A4=9A=E6=AC=A1=E8=AE=A1?= =?UTF-8?q?=E7=AE=97=EF=BC=8C=E6=8A=8A=E5=8E=9F=E6=9D=A5=E7=9A=84=E5=86=85?= =?UTF-8?q?=E7=82=B9=E6=B3=95=E5=BE=AA=E7=8E=AF=E5=8D=95=E7=8B=AC=E6=94=BE?= =?UTF-8?q?=E5=88=B0=E4=B8=80=E4=B8=AA=E6=96=87=E4=BB=B6=E4=B8=AD=E3=80=82?= =?UTF-8?q?=202.=E6=B2=A1=E6=9C=89=E9=87=8F=E6=B5=8B=E9=87=8F=E7=9A=84?= =?UTF-8?q?=E6=97=B6=E5=80=99=EF=BC=8C=E8=BF=98=E6=98=AF=E8=A6=81=E7=BA=A6?= =?UTF-8?q?=E6=9D=9F=E4=B8=80=E4=B8=8B=E7=9A=84=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: dugg@lab-desk --- OPF.m | 223 ++++++++++++++--------------------------------------- OPF_Init.m | 4 +- subOPF.m | 190 +++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 252 insertions(+), 165 deletions(-) create mode 100644 subOPF.m diff --git a/OPF.m b/OPF.m index d5419ff..1f530cc 100644 --- a/OPF.m +++ b/OPF.m @@ -1,174 +1,71 @@ -tic clc clear -%% 存在问题 -% 变压器变比的位置没有考虑，由于现在用的变比都是1，所以没有影响。 20130123 -%% -thesis=ForThesis(1,62); +close [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:\算例\17\17.csv'); -% pf('E:\算例\柳金Ⅰ926_21671693_2012-09-06\newFIle16.txt'); -%pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt'); -%pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt'); -%pf('c:/file31.txt'); -%% 计算功率因数 -Loadi=QD~=0 | PD~=0; -PF=sqrt(PD(Loadi).^2./(QD(Loadi).^2+PD(Loadi).^2)); -%% -Volt; -UAngel*180/3.1415926; -%% 通过潮流计算PG -AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); -PGBal=PD+diag(Volt)*Y.*cos(AngleIJ)*Volt'; -QGBal=QD+diag(Volt)*Y.*sin(AngleIJ)*Volt'; -%% 初值-即测量值 -PG0=PG; -QG0=QG; -PD0=PD; -QD0=QD; -Volt0=Volt; -UAngel0=UAngel; -rVolt=Volt0; -%% -PG0(Balance)=PGBal(Balance); -PG(Balance)=PGBal(Balance); -QG0(Balance)=QGBal(Balance); -QG0(PVi)=QGBal(PVi); -QG(PVi)=QGBal(PVi); -%% 真实值 -RealPG=PG0; -RealQG=QG0; -RealPD=PD0; -RealQD=QD0; -%% -[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,RealPD,RealQD,QD,PD); -Gap=(Init_L*Init_Z'-Init_U*Init_W'); -KK=0; -plotGap=zeros(1,60); -ContrlCount=size(Loadi,1)*2+Busnum*2+Busnum+length(Loadi)*2; -kmax=100; -Precision=Precision/1; -%% 加误差 sigma=0.01; +RealPD=PD; +RealQD=QD; +rVolt=Volt; +Loadi=find(PD~=0); +PD0=sparse(Busnum,1); +QD0=sparse(Busnum,1); PD0(Loadi)=RealPD(Loadi).*(1+normrnd(0,sigma,length(Loadi),1)); QD0(Loadi)=RealQD(Loadi).*(1+normrnd(0,sigma,length(Loadi),1)); mVolt=rVolt.*(1+normrnd(0,sigma,length(rVolt),1))'; -mPD=PD0; -mQD=QD0; -%加估计上下界 -lPD=abs(RealPD*3*sigma); -uPD=abs(RealPD*3*sigma); -lQD=abs(RealQD*3*sigma); -uQD=abs(RealQD*3*sigma); -lVolt=abs(rVolt'*3*sigma); -uVolt=abs(rVolt'*3*sigma); -%错误数据 -%mVolt(2)=5; -bigM=10; -Vbi=sparse(0.5*ones(Busnum,1)); -PDbi=sparse(0.5*ones(length(Loadi),1)); -QDbi=sparse(0.5*ones(length(Loadi),1)); -eps=10; -% 第一遍，算连续的值 -fprintf('第1次迭代，算连续量。\n'); -while(abs(Gap)>Precision*1) - if KK>kmax - break; - end - plotGap(KK+1)=Gap; - Init_u=Gap/2/RestraintCount*CenterA; - AngleIJMat=0; - %% 开始计算OPF - %% 形成等式约束的雅克比 - deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi); - %% 形成不等式约束的雅克比 - deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD,Vbi,PDbi,QDbi); - %% - L_1Z=diag(Init_Z./Init_L); - U_1W=diag(Init_W./Init_U); - %% 形成海森阵 - deltdeltF=func_deltdeltF(PVi,wPG,wQG,wPD,wQD,ContrlCount); - %% 形成ddHy - ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount); - %% 开始构建ddg - ddgzw=func_ddg(Busnum,ContrlCount,Loadi,Init_Z,Init_W); - %% 开始构建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); - Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD,uPD,lPD,uQD,lQD,uVolt,lVolt); - Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi); - Ly=Mat_H; - Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,RealPD,RealQD,Loadi,KK,PF,eps); - Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,eps); - Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); - YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx); - %% 开始解方程 - fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1)); - XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddgzw,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,RestraintCount,Lx,Balance,PVi,PGi,Busnum,Loadi); - %%取各分量 - [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,Busnum); - [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD,Vbi,PDbi,QDbi]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel,PVi,ContrlCount,Balance,Busnum,PGi,PD,QD,Loadi,Vbi,PDbi,QDbi); - Gap=(Init_L*Init_Z'-Init_U*Init_W'); - KK=KK+1; -end -% 第二遍，算离散 -Gap=1000; -KK=0; -fprintf('\n'); -fprintf('第2次迭代，算离散量。\n'); -while(abs(Gap)>Precision*1) - if KK>kmax - break; - end - plotGap(KK+1)=Gap; - Init_u=Gap/2/RestraintCount*CenterA; - AngleIJMat=0; - %% 开始计算OPF - %% 形成等式约束的雅克比 - deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi); - %% 形成不等式约束的雅克比 - deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD,Vbi,PDbi,QDbi); - %% - L_1Z=diag(Init_Z./Init_L); - U_1W=diag(Init_W./Init_U); - %% 形成海森阵 - deltdeltF=func_deltdeltF(PVi,wPG,wQG,wPD,wQD,ContrlCount); - %% 形成ddHy - ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount); - %% 开始构建ddg - ddgzw=func_ddg(Busnum,ContrlCount,Loadi,Init_Z,Init_W); - %% 开始构建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); - Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD,uPD,lPD,uQD,lQD,uVolt,lVolt); - Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi); - Ly=Mat_H; - if KK>8 - eps=eps*0.3; - if abs(eps)<1e-6 - eps=1e-5; - end - eps; - end - Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,RealPD,RealQD,Loadi,KK,PF,eps); - Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,eps); - Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); - YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx); - %% 开始解方程 - fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1)); - XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddgzw,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,RestraintCount,Lx,Balance,PVi,PGi,Busnum,Loadi); - %%取各分量 - [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,Busnum); - [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD,Vbi,PDbi,QDbi]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel,PVi,ContrlCount,Balance,Busnum,PGi,PD,QD,Loadi,Vbi,PDbi,QDbi); - Gap=(Init_L*Init_Z'-Init_U*Init_W'); - KK=KK+1; -end -toc -eps +%% 画Case A的图 +figV=figure(); +[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi]=subOPF([],PD0,QD0,mVolt,sigma);%全部有 +subplot(4,1,1); +plot(1:Busnum,Volt-rVolt,'b.:','Marker','diamond'); +subplot(4,1,2); +plot(1:Busnum,UAngel-rUAngel,'b:','Marker','diamond'); +subplot(4,1,3); +plot(1:Busnum,(PD-RealPD)./(RealPD+0.00001),'b:','Marker','diamond'); +subplot(4,1,4); +plot(1:Busnum,(QD-RealQD)./(RealQD+0.00001),'b:','Marker','diamond'); +%% Case B +[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi]=subOPF(cat(2,[1,2,3,4,5,11,12,13,16],[18,19]),PD0,QD0,mVolt,sigma);% +% [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi]=subOPF(cat(2,[1,2,3,4,6,7],[]),PD0,QD0,mVolt,sigma);% +subplot(4,1,1); +hold on; +plot(1:Busnum,Volt-rVolt,'g.:','Marker','square'); +subplot(4,1,2); +hold on; +plot(1:Busnum,UAngel-rUAngel,'g:','Marker','square'); +subplot(4,1,3); +hold on; +plot(1:Busnum,(PD-RealPD)./(RealPD+0.00001),'g:','Marker','square'); +subplot(4,1,4); +hold on; +plot(1:Busnum,(QD-RealQD)./(RealQD+0.00001),'g:','Marker','square'); +%% Case C +[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi]=subOPF([1:17],PD0,QD0,mVolt,sigma);% +subplot(4,1,1); +hold on; +plot(1:Busnum,Volt-rVolt,'r.:','Marker','o'); +subplot(4,1,2); +hold on; +plot(1:Busnum,UAngel-rUAngel,'r:','Marker','o'); +subplot(4,1,3); +hold on; +plot(1:Busnum,(PD-RealPD)./(RealPD+0.00001),'r:','Marker','o'); +subplot(4,1,4); +hold on; +plot(1:Busnum,(QD-RealQD)./(RealQD+0.00001),'r:','Marker','o'); +% 画legend +subplot(4,1,1); +title('Voltage'); +legend('Case A','Case B','Case C'); +subplot(4,1,2); +title('Voltage Angle'); +legend('Case A','Case B','Case C'); +subplot(4,1,3); +title('Active load power'); +legend('Case A','Case B','Case C'); +subplot(4,1,4); +title('Reactive load power'); +legend('Case A','Case B','Case C'); obj=sum(Vbi)+sum(PDbi)+sum(QDbi); fprintf('目标函数值 %.2f\n',full(obj)); diff --git a/OPF_Init.m b/OPF_Init.m index 7f9ddba..816207b 100644 --- a/OPF_Init.m +++ b/OPF_Init.m @@ -29,8 +29,8 @@ wQD(1:2:end)=0; %wD(randPDind)=0;%一些负荷不约束 %wD(7)=0; % wD(11)=0; -PD=0.8*PD0; +PD=0.0*PD0; %powerFacter=0.98; %QD=PD*sqrt((1-powerFacter^2)/powerFacter^2); -QD=.8*QD0; +QD=.0*QD0; end \ No newline at end of file diff --git a/subOPF.m b/subOPF.m new file mode 100644 index 0000000..e13702c --- /dev/null +++ b/subOPF.m @@ -0,0 +1,190 @@ +function [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi]=subOPF(noMeasurei,PD0,QD0,mVolt,sigma) +tic +%% 存在问题 +% 变压器变比的位置没有考虑，由于现在用的变比都是1，所以没有影响。 20130123 +%% +thesis=ForThesis(1,62); +[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:\算例\17\17.csv'); +% pf('E:\算例\柳金Ⅰ926_21671693_2012-09-06\newFIle16.txt'); +%pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt'); +%pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt'); +%pf('c:/file31.txt'); +%% 计算功率因数 +Loadi=QD~=0 | PD~=0; +PF=sqrt(PD(Loadi).^2./(QD(Loadi).^2+PD(Loadi).^2)); +%% +% Volt; +% UAngel*180/3.1415926; +%% 通过潮流计算PG +AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); +PGBal=PD+diag(Volt)*Y.*cos(AngleIJ)*Volt'; +QGBal=QD+diag(Volt)*Y.*sin(AngleIJ)*Volt'; +%% 真实值 +% PG0=PG; +% QG0=QG; +% PD0=PD; +% QD0=QD; +% Volt0=Volt; +% UAngel0=UAngel; +RealPD=PD; +RealQD=QD; +rUAngel=UAngel; +rVolt=Volt; +%% +PG0(Balance)=PGBal(Balance); +PG(Balance)=PGBal(Balance); +QG0(Balance)=QGBal(Balance); +QG0(PVi)=QGBal(PVi); +QG(PVi)=QGBal(PVi); + + +%% +[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,PD,PD00,QD,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,RealPD,RealQD,QD,PD); +Gap=(Init_L*Init_Z'-Init_U*Init_W'); +KK=0; +plotGap=zeros(1,60); +ContrlCount=size(Loadi,1)*2+Busnum*2+Busnum+length(Loadi)*2; +kmax=100; +Precision=Precision/1; +%% 加误差 +% sigma=0.03; +% PD0(Loadi)=RealPD(Loadi).*(1+normrnd(0,sigma,length(Loadi),1)); +% QD0(Loadi)=RealQD(Loadi).*(1+normrnd(0,sigma,length(Loadi),1)); +% mVolt=rVolt.*(1+normrnd(0,sigma,length(rVolt),1))'; +%找DG +DGi=find(PD0<0); +% PD0(DGi)=PD0(DGi).*(1+normrnd(0,sigma*0.1,length(DGi),1)); +% QD0(DGi)=QD0(DGi).*(1+normrnd(0,sigma*0.1,length(DGi),1)); +% mVolt(DGi)=mVolt(DGi).*(1+normrnd(0,sigma*0.1,length(DGi),1))'; +% +mPD=PD0; +mQD=QD0; +%加估计上下界 +lPD=abs(RealPD*3*sigma); +uPD=abs(RealPD*3*sigma); +lQD=abs(RealQD*3*sigma); +uQD=abs(RealQD*3*sigma); +lVolt=abs(mVolt'*3*sigma); +uVolt=abs(mVolt'*3*sigma); +%DG +lPD(DGi)=abs(RealPD(DGi)*3*sigma*0.1); +uPD(DGi)=abs(RealPD(DGi)*3*sigma*0.1); +lQD(DGi)=abs(RealQD(DGi)*3*sigma*0.1); +uQD(DGi)=abs(RealQD(DGi)*3*sigma*0.1); +lVolt(DGi)=abs(mVolt(DGi)'*3*sigma*0.1); +uVolt(DGi)=abs(mVolt(DGi)'*3*sigma*0.1); +%% 不测量的数据 +lPD(noMeasurei)=0.1*mPD(noMeasurei); +uPD(noMeasurei)=0.1*mPD(noMeasurei); +lQD(noMeasurei)=0.1*mQD(noMeasurei); +uQD(noMeasurei)=0.1*mPD(noMeasurei); +lVolt(noMeasurei)=0.1*mVolt(noMeasurei); +uVolt(noMeasurei)=0.1*mVolt(noMeasurei); +%错误数据 +%mVolt(2)=5; +bigM=10; +Vbi=sparse(0.5*ones(Busnum,1)); +PDbi=sparse(0.5*ones(length(Loadi),1)); +QDbi=sparse(0.5*ones(length(Loadi),1)); +eps=10; +% 第一遍，算连续的值 +fprintf('第1次迭代，算连续量。\n'); +while(abs(Gap)>Precision*1) + if KK>kmax + break; + end + plotGap(KK+1)=Gap; + Init_u=Gap/2/RestraintCount*CenterA; + AngleIJMat=0; + %% 开始计算OPF + %% 形成等式约束的雅克比 + deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi); + %% 形成不等式约束的雅克比 + deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD,Vbi,PDbi,QDbi); + %% + L_1Z=diag(Init_Z./Init_L); + U_1W=diag(Init_W./Init_U); + %% 形成海森阵 + deltdeltF=func_deltdeltF(PVi,wPG,wQG,wPD,wQD,ContrlCount); + %% 形成ddHy + ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount); + %% 开始构建ddg + ddgzw=func_ddg(Busnum,ContrlCount,Loadi,Init_Z,Init_W); + %% 开始构建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); + Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD,uPD,lPD,uQD,lQD,uVolt,lVolt); + Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi); + Ly=Mat_H; + Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,RealPD,RealQD,Loadi,KK,PF,eps); + Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,eps); + Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); + YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx); + %% 开始解方程 + fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1)); + XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddgzw,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,RestraintCount,Lx,Balance,PVi,PGi,Busnum,Loadi); + %%取各分量 + [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,Busnum); + [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD,Vbi,PDbi,QDbi]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel,PVi,ContrlCount,Balance,Busnum,PGi,PD,QD,Loadi,Vbi,PDbi,QDbi); + Gap=(Init_L*Init_Z'-Init_U*Init_W'); + KK=KK+1; +end +% 第二遍，算离散 +Gap=1000; +KK=0; +fprintf('\n'); +fprintf('第2次迭代，算离散量。\n'); +while(abs(Gap)>Precision*1) + if KK>kmax + break; + end + plotGap(KK+1)=Gap; + Init_u=Gap/2/RestraintCount*CenterA; + AngleIJMat=0; + %% 开始计算OPF + %% 形成等式约束的雅克比 + deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi); + %% 形成不等式约束的雅克比 + deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD,Vbi,PDbi,QDbi); + %% + L_1Z=diag(Init_Z./Init_L); + U_1W=diag(Init_W./Init_U); + %% 形成海森阵 + deltdeltF=func_deltdeltF(PVi,wPG,wQG,wPD,wQD,ContrlCount); + %% 形成ddHy + ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount); + %% 开始构建ddg + ddgzw=func_ddg(Busnum,ContrlCount,Loadi,Init_Z,Init_W); + %% 开始构建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); + Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD,uPD,lPD,uQD,lQD,uVolt,lVolt); + Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi); + Ly=Mat_H; + if KK>8 + eps=eps*0.3; + if abs(eps)<1e-6 + eps=1e-5; + end + eps; + end + Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,RealPD,RealQD,Loadi,KK,PF,eps); + Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,eps); + Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); + YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx); + %% 开始解方程 + fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1)); + XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddgzw,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,RestraintCount,Lx,Balance,PVi,PGi,Busnum,Loadi); + %%取各分量 + [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,Busnum); + [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD,Vbi,PDbi,QDbi]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel,PVi,ContrlCount,Balance,Busnum,PGi,PD,QD,Loadi,Vbi,PDbi,QDbi); + Gap=(Init_L*Init_Z'-Init_U*Init_W'); + KK=KK+1; +end +toc +end \ No newline at end of file