facat
/
pes2014
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Compare commits

base: facat:master
Branches Tags
facat:master
facat:20140513
facat:thesis
facat:withVolt
facat:dev
..
compare: facat:thesis
Branches Tags
facat:master
facat:20140513
facat:thesis
facat:withVolt
facat:dev

2 Commits
master ... thesis

Author SHA1 Message Date
dmy 77a46db096 用了赵晓慧的办法，没什么效果。 
Signed-off-by: dmy <dugg@21cn.com>
 2014-05-13 21:47:02 +08:00
dmy bfd2cd8240 改好了目标函数，收敛了。 
Signed-off-by: dmy <dugg@21cn.com>
 2014-05-13 16:53:38 +08:00
70 changed files with 641 additions and 1873 deletions
Show Stats Expand all files Collapse all files

2
.gitignore vendored
View File

@ -1,2 +0,0 @@
*.asv
*.7z

0
@ForThesis1/AddPDQDPGQG.m → @ForThesis/AddPDQDPGQG.m
View File

0
@ForThesis1/ForThesis.m → @ForThesis/ForThesis.m
View File

41
@ForThesis/MaxBranchDeviation.asv Normal file
View File

@ -0,0 +1,41 @@
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,PD,QD,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);
t2=(dispTransloss0(:,3) - dispTransloss(:,3))./dispTransloss0(:,3);
t11=dispLineloss0(:,3)>1e-5;% 太小的值不计算
t22=dispTransloss0(:,3)>1e-5;% 太小的值不计算
t3=abs([t1(t11);t2(t22)]);
output_args=max(t3(t3~=Inf));
end

0
@ForThesis1/MaxBranchDeviation.m → @ForThesis/MaxBranchDeviation.m
View File

7
@ForThesis/MaxDeviation.asv Normal file
View File

@ -0,0 +1,7 @@
function MaxDeviation(this,PG0,QG0,PD0,QD0)
PD0Array=repmat(PD0,1,this.sampleNum);
QD0Array=repmat(QD0,1,this.sampleNum);
PDMaxDev=max(abs(this.PDArray-PD0Array),[],2)
QDMaxDev=max(abs(this.QDArray-QD0Array),[],2)
PG0Array=repmat()
end

0
@ForThesis1/MaxDeviation.m → @ForThesis/MaxDeviation.m
View File

0
@ForThesis1/MeanPD.m → @ForThesis/MeanPD.m
View File

0
@ForThesis1/MeanPG.m → @ForThesis/MeanPG.m
View File

0
@ForThesis1/MeanQD.m → @ForThesis/MeanQD.m
View File

0
@ForThesis1/MeanQG.m → @ForThesis/MeanQG.m
View File

0
@ForThesis1/PercentOfPass.m → @ForThesis/PercentOfPass.m
View File

0
@ForThesis1/SqareDeviation.m → @ForThesis/SqareDeviation.m
View File

0
@ForThesis1/StatBranchDeviation.m → @ForThesis/StatBranchDeviation.m
View File

15
@ForThesis/StatDeviation.asv Normal file
View File

@ -0,0 +1,15 @@
function [ output_args ] = StatDeviation( this,PG0,QG0,PD0,QD0 )%ͳ¼ÆÎó²î
%STATDEVIATION Summary of this function goes here
% Detailed explanation goes here
PD0Array=repmat(PD0,1,this.sampleNum);
QD0Array=repmat(QD0,1,this.sampleNum);
PDMaxDev=max(abs(this.PDArray-PD0Array),[],2);
QDMaxDev=max(abs(this.QDArray-QD0Array),[],2);
PG0Array=repmat(PG0,this.sampleNum,1);
QG0Array=repmat(QG0,this.sampleNum,1);
PGMaxDev=max(abs(PG0Array));
QGMaxDev=max(abs(QG0Array));
end

0
@ForThesis1/StatDeviation.m → @ForThesis/StatDeviation.m
View File

83
AngelBar.m
View File

@ -1,83 +0,0 @@
function createfigure(xvector1, yvector1, xvector2, yvector2, xvector3, yvector3)
%CREATEFIGURE(XVECTOR1,YVECTOR1,XVECTOR2,YVECTOR2,XVECTOR3,YVECTOR3)
% XVECTOR1: bar xvector
% YVECTOR1: bar yvector
% XVECTOR2: bar xvector
% YVECTOR2: bar yvector
% XVECTOR3: bar xvector
% YVECTOR3: bar yvector
% Auto-generated by MATLAB on 12-Feb-2015 13:48:54
% Create figure
figure1 = figure('Name','相角直方图','Color',[1 1 1]);
% Create axes
axes1 = axes('Parent',figure1,'XTick',[-0.0002 -0.0001 0 0.0001 0.0002],...
'Position',[0.0920853080568723 0.470619469026549 0.20649289099526 0.372300884955752],...
'FontSize',13);
% Uncomment the following line to preserve the X-limits of the axes
xlim(axes1,[-0.0003 0.0003]);
% Uncomment the following line to preserve the Y-limits of the axes
ylim(axes1,[0 11]);
box(axes1,'on');
hold(axes1,'all');
% Create bar
bar(xvector1,yvector1,'FaceColor',[0 0 0],'Parent',axes1);
% Create xlabel
xlabel('偏差/rad','FontSize',13);
% Create ylabel
ylabel('数量','FontSize',13);
% Create title
title('情形A','FontSize',13);
% Create axes
axes2 = axes('Parent',figure1,'XTick',[-0.002 -0.001 0 0.001 0.002],...
'Position',[0.350370561164917 0.466814159292035 0.20649289099526 0.372300884955752],...
'FontSize',13);
% Uncomment the following line to preserve the X-limits of the axes
xlim(axes2,[-0.002 0.002]);
% Uncomment the following line to preserve the Y-limits of the axes
ylim(axes2,[0 7]);
box(axes2,'on');
hold(axes2,'all');
% Create bar
bar(xvector2,yvector2,'FaceColor',[0 0 0],'Parent',axes2);
% Create xlabel
xlabel('偏差/rad','FontSize',13);
% Create title
title('情形B','FontSize',13);
% Create ylabel
ylabel('数量','FontSize',13);
% Create axes
axes3 = axes('Parent',figure1,'XTick',[-0.002 -0.001 0 0.001 0.002],...
'Position',[0.606286146026517 0.463982300884956 0.20649289099526 0.372300884955752],...
'FontSize',13);
% Uncomment the following line to preserve the X-limits of the axes
xlim(axes3,[-0.002 0.002]);
% Uncomment the following line to preserve the Y-limits of the axes
ylim(axes3,[0 10]);
box(axes3,'on');
hold(axes3,'all');
% Create bar
bar(xvector3,yvector3,'FaceColor',[0 0 0],'Parent',axes3);
% Create xlabel
xlabel('偏差/rad','FontSize',13);
% Create ylabel
ylabel('数量','FontSize',13);
% Create title
title('情形C','FontSize',13);

73
CalTimeBar.m
View File

@ -1,73 +0,0 @@
function createfigure(ymatrix1)
%CREATEFIGURE(YMATRIX1)
% YMATRIX1: bar matrix data
% Auto-generated by MATLAB on 11-Feb-2015 20:23:32
% Create figure
figure1 = figure('Name','最大不平衡量','Color',[1 1 1]);
% Create axes
axes1 = axes('Parent',figure1,'XTickLabel',{'情形A','情形B','情形C'},...
'XTick',[1 2 3],...
'Position',[0.13 0.11 0.29654028436019 0.391061571125265],...
'FontName','Times New Roman');
hold(axes1,'all');
% Create multiple lines using matrix input to bar
bar1 = bar(ymatrix1,'BarLayout','stacked');
set(bar1(1),...
'FaceColor',[0.831372559070587 0.815686285495758 0.7843137383461],...
'DisplayName','阶段1');
set(bar1(2),...
'FaceColor',[0.941176474094391 0.941176474094391 0.941176474094391],...
'DisplayName','fdsfsd');
% Create xlabel
xlabel('情形','FontSize',12,'FontName','黑体');
% Create ylabel
ylabel('时间/ms','FontSize',12,'FontName','黑体');
% Create textbox
annotation(figure1,'textbox',...
[0.166919431279624 0.371781316348195 0.0888625592417062 0.059447983014862],...
'String',{'29.68ms'},...
'FontName','Times New Roman',...
'LineStyle','none');
% Create textbox
annotation(figure1,'textbox',...
[0.16706161137441 0.262269639065817 0.0888625592417062 0.059447983014862],...
'String',{'70.16ms'},...
'FontName','Times New Roman',...
'LineStyle','none');
% Create textbox
annotation(figure1,'textbox',...
[0.243791469194317 0.353904458598726 0.0888625592417062 0.059447983014862],...
'String',{'31.66ms'},...
'FontName','Times New Roman',...
'LineStyle','none');
% Create textbox
annotation(figure1,'textbox',...
[0.243933649289103 0.244392781316348 0.0888625592417062 0.059447983014862],...
'String',{'68.48ms'},...
'FontName','Times New Roman',...
'LineStyle','none');
% Create textbox
annotation(figure1,'textbox',...
[0.314881516587683 0.343288747346072 0.0888625592417062 0.059447983014862],...
'String',{'30.08ms'},...
'FontName','Times New Roman',...
'LineStyle','none');
% Create textbox
annotation(figure1,'textbox',...
[0.315023696682469 0.233777070063694 0.0888625592417062 0.059447983014862],...
'String',{'65.16ms'},...
'FontName','Times New Roman',...
'LineStyle','none');

40
CreateMaxErrorFigure.m
View File

@ -1,40 +0,0 @@
function createfigure(ymatrix1)
%CREATEFIGURE(YMATRIX1)
% YMATRIX1: bar matrix data
% Auto-generated by MATLAB on 11-Feb-2015 21:08:10
% Create figure
figure1 = figure('Color',[1 1 1]);
% Create axes
axes1 = axes('Parent',figure1,...
'XTick',[0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 35],...
'Position',[0.127153700189753 0.112232142857143 0.430721062618596 0.385535714285714],...
'FontName','Times New Roman');
% Uncomment the following line to preserve the X-limits of the axes
xlim(axes1,[1 35]);
% Uncomment the following line to preserve the Y-limits of the axes
ylim(axes1,[0 0.0021]);
hold(axes1,'all');
% Create multiple lines using matrix input to bar
bar1 = bar(ymatrix1,'BarWidth',1,'Parent',axes1);
set(bar1(1),'FaceColor',[0 0 0],'DisplayName','电压幅值');
set(bar1(2),...
'FaceColor',[0.941176474094391 0.941176474094391 0.941176474094391],...
'DisplayName','电压相角');
% Create xlabel
xlabel('节点号','FontSize',12,'FontName','黑体');
% Create ylabel
ylabel('最大绝对偏差/(p.u，弧度)','FontSize',12,'FontName','黑体');
% Create legend
legend1 = legend(axes1,'show');
set(legend1,...
'Position',[0.367488931056294 0.487351190476188 0.115749525616698 0.09375],...
'FontSize',12,...
'FontName','黑体');

115
DeviationFigure.m
View File

@ -1,115 +0,0 @@
function createfigure(X1, YMatrix1, YMatrix2, YMatrix3, YMatrix4)
%CREATEFIGURE(X1,YMATRIX1,YMATRIX2,YMATRIX3,YMATRIX4)
% X1: vector of x data
% YMATRIX1: matrix of y data
% YMATRIX2: matrix of y data
% YMATRIX3: matrix of y data
% YMATRIX4: matrix of y data
% Auto-generated by MATLAB on 11-Feb-2015 16:31:59
% Create figure
figure1 = figure('Color',[1 1 1]);
% Create axes
axes1 = axes('Parent',figure1,'XTick',[1 4 7 10 13 16 19 22 25 28 31 33],...
'Position',[0.140663507109006 0.572508000217951 0.537061611374407 0.365998852384935],...
'FontName','Times New Roman');
hold(axes1,'all');
% Create multiple lines using matrix input to plot
plot1 = plot(X1,YMatrix1,'Parent',axes1,'LineStyle',':','Color',[0 0 0]);
set(plot1(1),'Marker','diamond','DisplayName','情形A');
set(plot1(2),'Marker','square','DisplayName','情形B');
set(plot1(3),'Marker','o','DisplayName','情形C');
% Create xlabel
xlabel('节点号','FontName','黑体');
% Create ylabel
ylabel('偏差/p.u','FontName','Times New Roman');
% Create axes
axes2 = axes('Parent',figure1,'XTick',[2 5 8 11 14 17 20 23 26 29 32],...
'Position',[0.143033175355449 0.120301205914361 0.537914691943128 0.365998852384935],...
'FontName','Times New Roman');
% Uncomment the following line to preserve the X-limits of the axes
% xlim(axes2,[2 33]);
hold(axes2,'all');
% Create multiple lines using matrix input to plot
plot2 = plot(X1,YMatrix3,'Parent',axes2,'LineStyle',':','Color',[0 0 0]);
set(plot2(1),'Marker','diamond','DisplayName','算例A');
set(plot2(2),'Marker','square','DisplayName','算例B');
set(plot2(3),'Marker','o','DisplayName','算例C');
% Create xlabel
xlabel('节点号','FontName','Times New Roman');
% Create ylabel
ylabel('相对偏差%','FontName','Times New Roman');
% Create title
title('有功');
% Create axes
axes3 = axes('Parent',figure1,...
'XTickLabel',{'1','4','7','10','13','16','19','22','25','28','31','33'},...
'XTick',[1 4 7 10 13 16 19 22 25 28 31 33],...
'Position',[0.440758293838865 0.574607379571982 0.537914691943128 0.365998852384935]);
hold(axes3,'all');
% Create multiple lines using matrix input to plot
plot3 = plot(X1,YMatrix2,'Parent',axes3,'LineStyle',':','Color',[0 0 0]);
set(plot3(1),'Marker','diamond','DisplayName','情形A');
set(plot3(2),'Marker','square','DisplayName','情形B');
set(plot3(3),'Marker','o','DisplayName','情形C');
% Create xlabel
xlabel('节点号');
% Create ylabel
ylabel('偏差/弧度');
% Create axes
axes4 = axes('Parent',figure1,'XTick',[2 5 8 11 14 17 20 23 26 29 32],...
'Position',[-0.490853080568712 -0.0598281717337482 0.537914691943128 0.365998852384935]);
% Uncomment the following line to preserve the X-limits of the axes
% xlim(axes4,[2 33]);
hold(axes4,'all');
% Create multiple lines using matrix input to plot
plot4 = plot(X1,YMatrix4,'Parent',axes4,'LineStyle',':','Color',[0 0 0]);
set(plot4(1),'Marker','diamond','DisplayName','算例A');
set(plot4(2),'Marker','square','DisplayName','算例B');
set(plot4(3),'Marker','o','DisplayName','算例C');
% Create xlabel
xlabel('节点号');
% Create ylabel
ylabel('相对偏差%');
% Create legend
legend1 = legend(axes1,'show');
set(legend1,...
'Position',[0.155301294649576 0.633902089695378 0.112559241706161 0.12206572769953],...
'FontName','黑体');
% Create legend
legend2 = legend(axes2,'show');
set(legend2,...
'Position',[0.533220206544619 0.178635859349444 0.112559241706161 0.12206572769953],...
'FontName','黑体');
% Create legend
legend3 = legend(axes3,'show');
set(legend3,...
'Position',[0.450281817919025 0.807479462339755 0.117298578199052 0.12206572769953],...
'FontName','黑体');
% Create legend
legend4 = legend(axes4,'show');
set(legend4,...
'Position',[-0.098296381133102 -0.0539564029130473 0.117298578199052 0.122065727699531]);

6
DrawGap.asv Normal file
View File

@ -0,0 +1,6 @@
function DrawGap(plotGap)
x=find(plotGap);
ts=size(x,2);
end

20
DrawLoadProfile.m
View File

@ -1,20 +0,0 @@
function [ ] = DrawLoadProfile( )
figure();
data=dlmread('LoadProfile.txt');
hold on;
y1=data(:,3);
y2=data(:,4);
x=1:96;
k=1;
for i=0:96-1
x1(k)=x(length(x)-i);
Z2(k)=y2(length(x)-i);
k=k+1;
end
h=patch([x x1],[y1' Z2],[230 230 230]./255)
set(h,'EdgeColor','None');
plot(1:96,data(:,2),'k');
plot(1:96,data(:,3),'--k');
plot(1:96,data(:,4),'--k');
end

12
FormAA.asv Normal file
View File

@ -0,0 +1,12 @@
function AA=FormAA(L_1Z,deltG,U_1W,Hcoma,deltH)
tOnes=eye(14);
tZeros=zeros(14);
AA=[
tOnes,L_1Z,tZeros,tZeros,tZeros,zeros(14,10);
tZeros,tOnes,tZeros,tZeros,-deltG',zeros(14,10);
tZeros,tZeros,tOnes,U_1W,tZeros,zeros(14,10);
tZeros,tZeros,tZeros,tOnes,deltG,zeros(14,10);
tZeros,tZeros,tZeros,tZeros,Hcoma,deltH;
zeros(10,14),zeros(10,14),zeros(10,14),zeros(10,14),deltH',zeros(10,10);
];
end

7
FormAlphaP.asv Normal file
View File

@ -0,0 +1,7 @@
function AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU)
ti=deltL(delt)
t1=-Init_L./deltL';
t2=-Init_U./deltU';
t3=[t1,t2];
t4=t3()
end

11
FormG.asv Normal file
View File

@ -0,0 +1,11 @@
function Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,QD,Loadi)
Mat_G=[
PG(PGi);
QG(PVi);
sparse(PD(Loadi));
sparse(QD(Loadi));
Volt';
];
end

17
FormG.m
View File

@ -1,17 +1,10 @@
function Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD,uPD,lPD,uQD,lQD,uVolt,lVolt)
function Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,QD,Loadi,Vbi)
Mat_G=[
sparse(PD(Loadi))-mPD(Loadi)-bigM*PDbi-uPD(Loadi);
sparse(PD(Loadi))-mPD(Loadi)+bigM*PDbi+lPD(Loadi);
sparse(QD(Loadi))-mQD(Loadi)-bigM*QDbi-uQD(Loadi);
sparse(QD(Loadi))-mQD(Loadi)+bigM*QDbi+lQD(Loadi);
Volt'-mVolt'-bigM*Vbi-uVolt;
Volt'-mVolt'+bigM*Vbi+lVolt;
sparse(PD(Loadi));
sparse(QD(Loadi));
Volt'-Vbi;
Volt'+Vbi;
Vbi;
PDbi;
QDbi;
(0.1-Vbi).*Vbi;
(0.1-PDbi).*PDbi;
(0.1-QDbi).*QDbi;
];
end

31
FormH.asv Normal file
View File

@ -0,0 +1,31 @@
function Mat_H=FormH(Busnum,GB,AngleIJMat,Volt,PG,PD,QG,QD,Y)
t1=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
t2=Volt'*Volt;
t3=t1.*t2;
t4=sum(-t3,2);%P
t5=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat);
t6=t2.*t5;
t7=sum(-t6,2);%Q
t8=PG-PD;
t9=QG-QD;
%Mat_H=([(PG-PD)',(QG-QD)'])'+([t4',t7'])';
Mat_H(1:2:2*Busnum)=t8(1:Busnum)+t4(1:Busnum);
Mat_H(2:2:2*Busnum)=t9(1:Busnum)+t7(1:Busnum);
Mat_H=Mat_H';
%%
QDcos=textread('300glys.txt');
t=QD(PD==0 &&);
aa=QD;
QD(QD~=0)=PD(QD~=0)./tan(QDcos);
QD(PD==0)=t;
%%
%%%%一下是学姐给的公式
AngleIJ=AngleIJMat-angle(GB);
%dP=PG-PD-diag(Volt)*Y*cos(AngleIJ)*Volt';
dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt';
dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt';
%Mat_H(1:2:2*Busnum)=dP(1:Busnum);暂时改一下 20111227
%Mat_H(2:2:2*Busnum)=dQ(1:Busnum);暂时改一下 20111227
Mat_H=[dP;dQ;];
end

2
FormH.m
View File

@ -2,5 +2,7 @@ function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi)
AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt';
dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt';
Mat_H=[dP;dQ;];
end

16
FormLw.m
View File

@ -1,7 +1,11 @@
function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,eps)
function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,noDataTransCapacity)
KK=999;
%PU=GenU(:,2);%发电机有功上界
VoltU=(1.1)*ones(1,Busnum);
%VoltU=10*ones(1,Busnum);
PDU=PD0(Loadi);
% PDU=noDataTransCapacity;
PDU(PDU>0)=1.200*PDU(PDU>0);
PDU(PDU<0)=0.800*PDU(PDU<0);
PDU(PDU==0)=0.400;
@ -9,6 +13,7 @@ realPD=PD0(Loadi);
indPD=find(realPD>0);
PDU(indPD(3:12:end))=1.55*realPD(indPD(3:12:end));
PDU(indPD(9:12:end))=1.05*realPD(indPD(9:12:end));
%PDU=10*ones(length(Loadi),1);
QDU=QD0(Loadi);
QDU(QDU>0)=1.200*QDU(QDU>0);
QDU(QDU<0)=0.800*QDU(QDU<0);
@ -17,12 +22,9 @@ realQD=QD0(Loadi);
indQD=find(realQD>0);
QDU(indQD(3:12:end))=1.55*realQD(indQD(3:12:end));
QDU(indQD(9:12:end))=1.05*realQD(indQD(9:12:end));
t1=([0*PDU',1*ones(1,length(PDU)),0*QDU',1*ones(1,length(QDU)),0*VoltU,1*ones(1,length(VoltU)),0.1*ones(1,Busnum+length(Loadi)*2)])';
t1=[t1;
sparse(eps*ones(Busnum,1));
sparse(eps*ones(length(Loadi),1));
sparse(eps*ones(length(Loadi),1));
];
% PF=0.85;
% QDU=1.0*PD(Loadi).*sqrt(1 -PF.^2)./PF;
t1=([PDU',QDU',1.1*VoltU,1.1*VoltU,1*ones(1,Busnum)])';
t2=Mat_G+Init_U'-t1;
Lw=t2;

5
FormLxComa.m
View File

@ -4,9 +4,8 @@ t2=Lul+diag(Init_Z)*Lz;
t3=inv(diag(Init_L));
t4=t3*t2;%
t5=Luu-diag(Init_W)*Lw;
% t6=inv(diag(Init_U));
% t7=t6*t5;%
t7=diag(Init_U)\t5;
t6=inv(diag(Init_U));
t7=t6*t5;%
t8=deltG*(t4+t7);%%
LxComa=Lx+t8;
end

10
FormLz.m
View File

@ -1,4 +1,4 @@
function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,PD0,QD0,Loadi,KK,PF,eps)
function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,PD0,QD0,Loadi,KK,PF)
KK=999;
VoltL=(0.9)*ones(1,Busnum);
@ -21,12 +21,8 @@ indQD=find(realQD>0);
QDL(indQD(3:12:end))=0.95*realQD(indQD(3:12:end));
QDL(indQD(9:12:end))=0.95*realQD(indQD(9:12:end));
% QDL=0*PD(Loadi).*sqrt((1-PF.^2))./PF;
t1=([-1*ones(1,length(PDL)),0*PDL',-1*ones(1,length(QDL)),0*QDL',-1*ones(1,length(VoltL)),0*VoltL,0*ones(1,Busnum+length(Loadi)*2)])';
t1=[t1;
sparse(-eps*ones(Busnum,1));
sparse(-eps*ones(length(Loadi),1));
sparse(-eps*ones(length(Loadi),1));
];
%t1=([PDL',QDL',VoltL,VoltL,0*ones(1,Busnum)])';
t1=([PDL',QDL',-20000*VoltL,-20000*VoltL,0*ones(1,Busnum)])';
t2=Mat_G-Init_L'-t1;
Lz=t2;

5
FormYY.asv Normal file
View File

@ -0,0 +1,5 @@
function FormYY(Init_L,Lul.Lz,Init_U,Luu,Lw,Lz,LxComa)
t=[
-
]
end

18
LoadProfile.m
View File

@ -1,18 +0,0 @@
txt=dlmread('E:\书籍\论文\配电网状态估计\负荷.txt');
txt=txt./10000;
x=1:length(txt);
plot(x,txt,'k');
hold on
plot(x,txt*0.8,'k--');
plot(x,txt*1.2,'k--');
legend('Load profile','Lower/Upper bound')
ylabel('Load /kW')
xlabel('Time /h')
k=1;
for i=0:length(x)-1
x1(k)=x(length(x)-i);
Z2(k)=txt(length(x)-i)*1.2;
k=k+1;
end
fill([x x1],[txt'*0.8 Z2],[240, 240, 240]./255,'edgealpha',0)

9
MaxBoundErrorFigure.m
View File

@ -1,9 +0,0 @@
function MaxBoundErrorFigure()
figure();
data=dlmread('maxBoundError.txt');
% data(1,:)=0;
maxV=data(:,1);
maxA=data(:,2);
% bar([maxV';maxA']');
CreateMaxErrorFigure([maxV';maxA']');
end

9
MaxErrorFigure.m
View File

@ -1,9 +0,0 @@
function MaxErrorFigure()
figure();
data=dlmread('maxError.txt');
data(1,:)=0;
maxV=data(:,1);
maxA=data(:,2);
% bar([maxV';maxA']');
CreateMaxErrorFigure([maxV';maxA']');
end

18
Modification.asv Normal file
View File

@ -0,0 +1,18 @@
function [Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel]=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)
AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU);
fprintf('AlphaP %f\n',AlphaP);
AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW);
fprintf('AlphaD %f\n',AlphaD);
Init_Z=Init_Z+AlphaD*deltZ';
Init_L=Init_L+AlphaP*deltL';
Init_W=Init_W+AlphaD*deltW';
Init_U=Init_U+AlphaP*deltU';
Init_Y=Init_Y+AlphaD*deltY';
PG(PVi)=PG(PGi)+deltX(size(PGi,1));
QG(PVi)=QG(PVi)+deltX(size(PGi,1)+1,1:(size(PVi,1)+size(PGi,1)));
t=deltX(size(PVi,1)+size(PGi,1)+1,ContrlCount)';
t(2*Balance-1)=0;
Volt=Volt+AlphaP*t(2:2:2*Busnum);
UAngel=UAngel+AlphaP*t(1:2:2*Busnum);
end

14
Modification.m
View File

@ -1,16 +1,8 @@
function [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)
function [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD,Vbi]=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)
AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU);
%fprintf('AlphaP %f\n',full(AlphaP));
AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW);
%fprintf('AlphaD %f\n',full(AlphaD));
% if AlphaP<1e-4
% AlphaP=0.1;
% end
%
% if AlphaD<1e-4
% AlphaD=0.1;
% end
Init_Z=Init_Z+AlphaD*deltZ';
Init_L=Init_L+AlphaP*deltL';
@ -32,7 +24,5 @@ balVolt=Volt(Balance);
Volt=Volt+AlphaP*t(1:Busnum);
Volt(Balance)=balVolt;
UAngel=UAngel+AlphaP*t(Busnum+1:2*Busnum);
Vbi=Vbi+AlphaP*t(2*Busnum+1:2*Busnum+Busnum)';
PDbi=PDbi+AlphaP*t(3*Busnum+1:3*Busnum+length(Loadi))';
QDbi=QDbi+AlphaP*t(3*Busnum+length(Loadi)+1:end)';
Vbi=Vbi+AlphaP*t(2*Busnum+1:end)';
end

151
OPF.asv Normal file
View File

@ -0,0 +1,151 @@
tic
clc
clear
%% 存在问题
% 变压器变比的位置没有考虑，由于现在用的变比都是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:\算例\东际911_2751267_2012-09-05\newFIle20-使用22.txt');
% 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;
%%
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;
kmax=60;
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));
Vbi=sparse(1*ones(Busnum,1));
kInit_Z=zeros(length(Init_Z),4);
kInit_L=zeros(length(Init_L),4);
kInit_W=zeros(length(Init_W),4);
kInit_U=zeros(length(Init_U),4);
kInit_Y=zeros(length(Init_Y),4);
kPG=zeros(length(PG),4);
kQG=zeros(length(QG),4);
kVolt=zeros(length(Volt),4);
kUAngel=zeros(length(UAngel),4);
kPD=zeros(length(PD),4);
kQD=zeros(length(QD),4);
kVbi=zeros(length(Vbi),4);
kdeltZ=zeros(length(Init_Z),4);
kdeltL=zeros(length(Init_L),4);
kdeltW=zeros(length(Init_W),4);
kdeltU=zeros(length(Init_U),4);
kdeltX=zeros(ContrlCount,4);
kdeltY=zeros(length(Init_Y),4);
while(abs(Gap)>Precision)
if KK>kmax
break;
end
plotGap(KK+1)=Gap;
for I=1:4
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);
%%
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
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);
Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,QD,Loadi,Vbi);
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);
Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,RealPD,RealQD,Loadi,KK,PF);
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,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,RestraintCount,Lx,Balance,PVi,PGi,Busnum,Loadi);
%%取各分量
[kdeltZ(:,I),kdeltL(:,I),kdeltW(:,I),kdeltU(:,I),kdeltX(:,I),kdeltY(:,I)]=AssignXX(XX,ContrlCount,RestraintCount,Busnum);
if I==1
stepT=0.5;
end
if I==2
stepT=0.5;
end
if I==3
stepT=1;
end
if I<=3
[kInit_Z(:,I),kInit_L(:,I),kInit_W(:,I),kInit_U(:,I),kInit_Y(:,I), ...
kPG(:,I),kQG(:,I),kVolt(:,I),kUAngel(:,I),kPD(:,I),kQD(:,I),kVbi(:,I)]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,stepT*kdeltZ(:,I),stepT*kdeltL(:,I),stepT*kdeltW(:,I),stepT*kdeltU(:,I),stepT*kdeltX(:,I),stepT*kdeltY(:,I),PG,QG,Volt,UAngel,PVi,ContrlCount,Balance,Busnum,PGi,PD,QD,Loadi,Vbi);
Init_Z=kInit_Z(:,I)';
Init_L=kInit_L(:,I)';
Init_W=kInit_W(:,I)';
Init_U=kInit_U(:,I)';
Init_Y=kInit_Y(:,I)';
PG=kPG(:,I);
QG=kQG(:,I);
Volt=kVolt(:,I)';
UAngel=kUAngel(:,I)';
PD=kPD(:,I);
QD=kQD(:,I);
Vbi=kVbi(:,I);
end
if I==4
deltZ=(kdeltZ(:,1)+kdeltZ(:,2)*2+kdeltZ(:,3)*2+kdeltZ(:,4))/6;
deltL=(kdeltL(:,1)+kdeltL(:,2)*2+kdeltL(:,3)*2+kdeltL(:,4))/6;
deltW=(kdeltW(:,1)+kdeltW(:,2)*2+kdeltW(:,3)*2+kdeltW(:,4))/6;
deltU=kdeltU(:,1)+kdeltU(:,2)*2+kdeltU(:,3)*2+kdeltU(:,4));
deltX=kdeltX(:,1)+kdeltX(:,2)*2+kdeltX(:,3)*2+kdeltX(:,4);
deltY=kdeltY(:,1)+kdeltY(:,2)*2+kdeltY(:,3)*2+kdeltY(:,4);
[Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD,Vbi]=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);
end
end
Gap=(Init_L*Init_Z'-Init_U*Init_W');
KK=KK+1;
end
fprintf('%f\n',sum(full(Vbi)));
toc

547
OPF.m
View File

@ -1,406 +1,151 @@
tic
clc
clear
close all
arrayA=zeros(21,10);
sumCaseA_SE=0;
sumCaseB_SE=0;
sumCaseC_SE=0;
badDataResult=zeros(10,33);
badDataLocation=zeros(34,10);
nodeMaxDVolt=zeros(33,1);
nodeMaxDVAngle=zeros(33,1);
for badDataNode=1:1
loopN=1;
maxDVolt=0;
maxDVAngle=0;
while 1
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:\算例\feeder33\feeder33.txt');
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))';
%加载保存的变量
PD0=load('PD0');
PD0=PD0.PD0;
QD0=load('QD0');
QD0=QD0.QD0;
mVolt=load('mVolt');
mVolt=mVolt.mVolt;
% mVolt(3)=rVolt(3)*(1-sigma*6);
%% 画Case A的图
% figure('Color',[1 1 1]);
[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapA,isConverge]=subOPF([],PD0,QD0,mVolt,sigma);%全部有
% badDataResult(I,badDataNode)=sum(Vbi);
% badDataLocation(1:33,I)=Vbi;
% badDataLocation(34,I)=sum(abs((rVolt-Volt)./rVolt./length(rVolt)))+sum(abs( (UAngel(2:33)-rUAngel(2:33))./rUAngel(2:33)./length(rUAngel(2:33))));
break;
if isConverge==0
continue;
end
% 加最大偏差
maxDVolt_=max(abs((rVolt-Volt)));
if maxDVolt_>maxDVolt
maxDVolt=maxDVolt_;
end
maxDVAngle_=max(abs((UAngel(2:33)-rUAngel(2:33))));
if maxDVAngle_>maxDVAngle
maxDVAngle=maxDVAngle_;
end
nodeMaxDVolt_t=abs((rVolt-Volt))';
nodeMaxDVAngle_t=abs((UAngel-rUAngel))';
nodeMaxDVolt_t([1:18,23:33])=0;
nodeMaxDVAngle_t([1:18,23:33])=0;
nodeMaxDVolt(nodeMaxDVolt<nodeMaxDVolt_t)=nodeMaxDVolt_t(nodeMaxDVolt<nodeMaxDVolt_t);
nodeMaxDVAngle(nodeMaxDVAngle<nodeMaxDVAngle_t)=nodeMaxDVAngle_t(nodeMaxDVAngle<nodeMaxDVAngle_t);
if loopN>=500
% nodeMaxDVolt(badDataNode)=maxDVolt;
% nodeMaxDVAngle(badDataNode)=maxDVAngle;
break;
end
loopN=loopN+1;
%% 存在问题
% 变压器变比的位置没有考虑，由于现在用的变比都是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:\算例\东际911_2751267_2012-09-05\newFIle20-使用22.txt');
% 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;
%%
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;
kmax=60;
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));
Vbi=sparse(1*ones(Busnum,1));
kInit_Z=zeros(length(Init_Z),4);
kInit_L=zeros(length(Init_L),4);
kInit_W=zeros(length(Init_W),4);
kInit_U=zeros(length(Init_U),4);
kInit_Y=zeros(length(Init_Y),4);
kPG=zeros(length(PG),4);
kQG=zeros(length(QG),4);
kVolt=zeros(length(Volt),4);
kUAngel=zeros(length(UAngel),4);
kPD=zeros(length(PD),4);
kQD=zeros(length(QD),4);
kVbi=zeros(length(Vbi),4);
kdeltZ=zeros(length(Init_Z),4);
kdeltL=zeros(length(Init_L),4);
kdeltW=zeros(length(Init_W),4);
kdeltU=zeros(length(Init_U),4);
kdeltX=zeros(ContrlCount,4);
kdeltY=zeros(length(Init_Y),4);
while(abs(Gap)>Precision)
if KK>kmax
break;
end
plotGap(KK+1)=Gap;
for I=1:4
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);
%%
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
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);
Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,QD,Loadi,Vbi);
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);
Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,RealPD,RealQD,Loadi,KK,PF);
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,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,RestraintCount,Lx,Balance,PVi,PGi,Busnum,Loadi);
%%取各分量
[kdeltZ(:,I),kdeltL(:,I),kdeltW(:,I),kdeltU(:,I),kdeltX(:,I),kdeltY(:,I)]=AssignXX(XX,ContrlCount,RestraintCount,Busnum);
if I==1
stepT=0.5;
end
if I==2
stepT=0.5;
end
if I==3
stepT=1;
end
if I<=3
[kInit_Z(:,I),kInit_L(:,I),kInit_W(:,I),kInit_U(:,I),kInit_Y(:,I), ...
kPG(:,I),kQG(:,I),kVolt(:,I),kUAngel(:,I),kPD(:,I),kQD(:,I),kVbi(:,I)]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,stepT*kdeltZ(:,I),stepT*kdeltL(:,I),stepT*kdeltW(:,I),stepT*kdeltU(:,I),stepT*kdeltX(:,I),stepT*kdeltY(:,I),PG,QG,Volt,UAngel,PVi,ContrlCount,Balance,Busnum,PGi,PD,QD,Loadi,Vbi);
Init_Z=kInit_Z(:,I)';
Init_L=kInit_L(:,I)';
Init_W=kInit_W(:,I)';
Init_U=kInit_U(:,I)';
Init_Y=kInit_Y(:,I)';
PG=kPG(:,I);
QG=kQG(:,I);
Volt=kVolt(:,I)';
UAngel=kUAngel(:,I)';
PD=kPD(:,I);
QD=kQD(:,I);
Vbi=kVbi(:,I);
end
if I==4
deltZ=(kdeltZ(:,1)+kdeltZ(:,2)*2+kdeltZ(:,3)*2+kdeltZ(:,4))/6;
deltL=(kdeltL(:,1)+kdeltL(:,2)*2+kdeltL(:,3)*2+kdeltL(:,4))/6;
deltW=(kdeltW(:,1)+kdeltW(:,2)*2+kdeltW(:,3)*2+kdeltW(:,4))/6;
deltU=(kdeltU(:,1)+kdeltU(:,2)*2+kdeltU(:,3)*2+kdeltU(:,4))/6;
deltX=(kdeltX(:,1)+kdeltX(:,2)*2+kdeltX(:,3)*2+kdeltX(:,4))/6;
deltY=(kdeltY(:,1)+kdeltY(:,2)*2+kdeltY(:,3)*2+kdeltY(:,4))/6;
[Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD,Vbi]=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);
end
end
Gap=(Init_L*Init_Z'-Init_U*Init_W');
KK=KK+1;
end
% end
subplot(4,1,1,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33]);
CaseAREV=(Volt-rVolt);%Relative Error of Voltage in Case A
% CaseAREV=CaseAREV(2:end)*100;
%真实值的
plot(1:length(CaseAREV),(CaseAREV),'k.:','Marker','diamond');
%测量值的
% plot(1:length(CaseAREV),abs((mVolt-rVolt)*100),'c.:','Marker','diamond');
box off;
set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33])
xlabel('节点号');
ylabel('误差%');
subplot(4,1,2);
% CaseAREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case A
CaseAREA=(UAngel-rUAngel);%Relative Error of Angle in Case A
CaseAREA(1)=0;
%真实值的
plot(1:length(CaseAREA),(CaseAREA),'k:','Marker','diamond');
box off;
set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33])
xlabel('节点号');
ylabel('误差%');
subplot(4,1,3);
% CaseAREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case A
CaseAREP=(PD-RealPD)./RealPD*100;%Relative Error of PD in Case A
CaseAREP(1)=0;
%真实值的
plot(1:length(CaseAREP),(CaseAREP),'k:','Marker','diamond');
%测量值的
% plot(1:length(CaseAREV),abs((PD0-RealPD)./(RealPD+0.00001)*100),'c.:','Marker','diamond');
box off;
set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33])
xlabel('节点号');
ylabel('误差%');
subplot(4,1,4);
% CaseAREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case A
CaseAREQ=(QD-RealQD)./RealQD*100;%Relative Error of QD in Case A
CaseAREQ(1)=0;
%真实值的
plot(1:length(CaseAREQ),(CaseAREQ),'k:','Marker','diamond');
%测量值的
% plot(1:length(CaseAREV),abs((QD0-RealQD)./(RealQD+0.00001)*100),'c.:','Marker','diamond');
box off;
set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33])
xlabel('节点号');
ylabel('误差%');
%计算Case A的误差
CaseAE=sqrt((sum(CaseAREV.^2)+sum(CaseAREA.^2)+sum(CaseAREP.^2)+sum(CaseAREQ.^2))/132);
objA=full(sum(Vbi)+sum(PDbi)+sum(QDbi));
notZeros=find(PD0~=0);
CaseA_SE=sum(((RealPD(notZeros)-PD(notZeros))./(RealPD(notZeros)*1)).^2)+sum(((RealQD(notZeros)-QD(notZeros))./(RealQD(notZeros)*1)).^2)+sum(((Volt-rVolt)./(rVolt*1)).^2);
CaseA_SE=(CaseA_SE/(length(notZeros)*2+length(Volt)))^.5;
sumCaseA_SE=sumCaseA_SE+CaseA_SE;
% arrayA(1:19,I)=Vbi;
% arrayA(21,I)=CaseAE*1000;
%% Case B
% [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapB]=subOPF(setdiff(1:Busnum,[18,21,22,29]),PD0,QD0,mVolt,sigma);%
[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapB]=subOPF(setdiff(1:Busnum,[2,3,5,20,24,27,28,10,11,12,13]),PD0,QD0,mVolt,sigma);%
subplot(4,1,1);
hold on;
% CaseBREV=(Volt-rVolt)*100;%Relative Error of Voltage in Case B
CaseBREV=(Volt-rVolt);%Relative Error of Voltage in Case B
plot(1:length(CaseBREV),(CaseBREV),'b.:','Marker','square');
subplot(4,1,2);
hold on;
% CaseBREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case B
CaseBREA=(UAngel-rUAngel);%Relative Error of Angle in Case B
plot(1:length(CaseBREA),(CaseBREA),'b:','Marker','square');
subplot(4,1,3);
hold on;
% CaseBREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case B
CaseBREP=(PD-RealPD)./RealPD*100;%Relative Error of PD in Case B
RealPD(1)=0;
plot(1:length(CaseBREP),(CaseBREP),'b:','Marker','square');
subplot(4,1,4);
hold on;
% CaseBREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case B
CaseBREQ=(QD-RealQD)./RealQD*100;%Relative Error of QD in Case B
CaseBREQ(1)=0;
plot(1:length(CaseBREQ),(CaseBREQ),'b:','Marker','square');
CaseBE=sqrt((sum(CaseBREV.^2)+sum(CaseBREA.^2)+sum(CaseBREP.^2)+sum(CaseBREQ.^2))/132);
objB=full(sum(Vbi)+sum(PDbi)+sum(QDbi));
noMeasurei=[2,3,5,20,24,27,28,10,11,12,13];
Measurei=setdiff(2:33,[2,3,5,20,24,27,28,10,11,12,13]);
CaseB_SE=sum(((RealPD(Measurei)-PD(Measurei))./(RealPD(Measurei)*1)).^2)+sum(((RealQD(Measurei)-QD(Measurei))./(RealQD(Measurei)*1)).^2)+sum(((Volt(Measurei)-rVolt(Measurei))./(rVolt(Measurei)*1)).^2);
CaseB_SE=CaseB_SE+sum(((RealPD(noMeasurei)-PD(noMeasurei))./(RealPD(noMeasurei)*1)).^2)+sum(((RealQD(noMeasurei)-QD(noMeasurei))./(RealQD(noMeasurei)*1)).^2)+sum(((Volt(noMeasurei)-rVolt(noMeasurei))./(rVolt(noMeasurei)*1)).^2);
CaseB_SE=(CaseB_SE/(length(notZeros)+length(noMeasurei) +length(Volt)))^.5;
sumCaseB_SE=sumCaseB_SE+CaseB_SE;
%% Case C
[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapC]=subOPF([1:33],PD0,QD0,mVolt,sigma);%
subplot(4,1,1);
hold on;
% CaseCREV=(Volt-rVolt)*100;%Relative Error of Voltage in Case C
CaseCREV=(Volt-rVolt);
plot(1:length(CaseCREV),(CaseCREV),'r.:','Marker','o');
subplot(4,1,2);
hold on;
% CaseCREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case C
CaseCREA=(UAngel-rUAngel);
plot(1:length(CaseCREA),(CaseCREA),'r:','Marker','o');
subplot(4,1,3);
hold on;
% CaseCREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case C
CaseCREP=(PD-RealPD)./RealPD*100;%Relative Error of PD in Case C
CaseCREP(1);
plot(1:length(CaseCREP),(CaseCREP),'r:','Marker','o');
subplot(4,1,4);
hold on;
% CaseCREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case C
CaseCREQ=(QD-RealQD)./RealQD*100;%Relative Error of QD in Case C
CaseCREQ(1)=0;
plot(1:length(CaseCREQ),(CaseCREQ),'r:','Marker','o');
% 画legend
subplot(4,1,1);
% title('Voltage');
ld=legend('算例A','算例B','算例C');
set(ld,'Position',[0.847865087908145 0.786094477711244 0.0543595263724435 0.0605455755156354]);
subplot(4,1,2);
% title('Voltage Angle');
ld=legend('算例A','算例B','算例C');
set(ld,'Position',[0.847865087908145 0.586094477711244 0.0543595263724435 0.0605455755156354]);
subplot(4,1,3);
% title('Active load power');
ld=legend('算例A','算例B','算例C');
set(ld,'Position',[0.847865087908145 0.386094477711244 0.0543595263724435 0.0605455755156354]);
subplot(4,1,4);
% title('Reactive load power');
ld=legend('算例A','算例B','算例C');
set(ld,'Position',[0.847865087908145 0.186094477711244 0.0543595263724435 0.0605455755156354]);
CaseCE=sqrt((sum(CaseCREV.^2)+sum(CaseCREA.^2)+sum(CaseCREP.^2)+sum(CaseCREQ.^2))/132);
objC=full(sum(Vbi)+sum(PDbi)+sum(QDbi));
fprintf('%f\n',sum(full(Vbi)));
toc
CaseC_SE=sum(((RealPD(notZeros)-PD(notZeros))./(RealPD(notZeros)*1)).^2)+sum(((RealQD(notZeros)-QD(notZeros))./(RealQD(notZeros)*1)).^2)+sum(((Volt-rVolt)./(rVolt*1)).^2);
CaseC_SE=(CaseC_SE/(length(notZeros)*2+length(Volt)))^.5;
sumCaseC_SE=sumCaseC_SE+CaseC_SE;
% fprintf('目标函数值 %.2f\n',full(obj));
fprintf('Case A Case B Case C 的误差\n')
fprintf('%f %f %f \n',CaseAE,CaseBE,CaseCE);
fprintf('三个Case目标值\n')
fprintf('%f\t%f\t%f \n',objA,objB,objC)
%% 画测量值
% subplot(4,1,1);
% plot(1:Busnum,mVolt-rVolt,'k.:','Marker','pentagram')
% subplot(4,1,3);
% plot(1:Busnum,(PD0-RealPD)./(RealPD+0.00001),'k:','Marker','pentagram');
% subplot(4,1,4);
% plot(1:Busnum,(QD0-RealQD)./(RealQD+0.00001),'k:','Marker','pentagram');
%% 直方图
% 电压
figure('Name','电压直方图')
split_number=20;
%Case A
subplot(1,3,1)
y=CaseAREV;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
yy=hist(y,x); %计算各个区间的个数
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
bar(x,yy) %画出概率密度分布图
voltBarCaseAX=x;
voltBarCaseAY=yy;
xlabel('Error');
ylabel('Number of buses');
title('算例A');
% ylim([0 4])
%Case B
subplot(1,3,2)
y=CaseBREV;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
yy=hist(y,x); %计算各个区间的个数
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
bar(x,yy) %画出概率密度分布图
voltBarCaseBX=x;
voltBarCaseBY=yy;
xlabel('Error');
ylabel('Number of buses');
title('算例B');
% ylim([0 4])
%Case C
subplot(1,3,3)
y=CaseCREV;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
yy=hist(y,x); %计算各个区间的个数
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
bar(x,yy) %画出概率密度分布图
voltBarCaseCX=x;
voltBarCaseCY=yy;
xlabel('Error');
ylabel('Number of buses');
title('算例C');
% ylim([0 4])
figure('Name','相角直方图')
% 相角
split_number=20;
%Case A
subplot(2,2,1)
y=CaseAREA;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
yy=hist(y,x); %计算各个区间的个数
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
bar(x,yy) %画出概率密度分布图
angelBarCaseAX=x;
angelBarCaseAY=yy;
% ylim([0 4])
%Case B
subplot(2,2,2)
y=CaseBREA;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
yy=hist(y,x); %计算各个区间的个数
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
bar(x,yy) %画出概率密度分布图
angelBarCaseBX=x;
angelBarCaseBY=yy;
% ylim([0 4])
%Case C
subplot(2,2,3)
y=CaseCREA;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
yy=hist(y,x); %计算各个区间的个数
angelBarCaseCX=x;
angelBarCaseCY=yy;
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
bar(x,yy) %画出概率密度分布图
%PD
figure('Name','有功直方图')
split_number=20;
%Case A
subplot(2,2,1)
y=CaseAREP;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
yy=hist(y,x); %计算各个区间的个数
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
bar(x,yy) %画出概率密度分布图
PDBarCaseAX=x;
PDBarCaseAY=yy;
% ylim([0 4])
%Case B
subplot(2,2,2)
y=CaseBREP;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
yy=hist(y,x); %计算各个区间的个数
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
bar(x,yy) %画出概率密度分布图
PDBarCaseBX=x;
PDBarCaseBY=yy;
% ylim([0 4])
%Case C
subplot(2,2,3)
y=CaseCREP;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
yy=hist(y,x); %计算各个区间的个数
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
bar(x,yy) %画出概率密度分布图
PDBarCaseCX=x;
PDBarCaseCY=yy;
%QD
figure('Name','无功直方图')
split_number=20;
%Case A
subplot(2,2,1)
y=CaseAREQ;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
yy=hist(y,x); %计算各个区间的个数
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
bar(x,yy) %画出概率密度分布图
QDBarCaseAX=x;
QDBarCaseAY=yy;
% ylim([0 4])
%Case B
subplot(2,2,2)
y=CaseBREQ;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
yy=hist(y,x); %计算各个区间的个数
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
bar(x,yy) %画出概率密度分布图
QDBarCaseBX=x;
QDBarCaseBY=yy;
% ylim([0 4])
%Case C
subplot(2,2,3)
y=CaseCREQ;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
yy=hist(y,x); %计算各个区间的个数
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
bar(x,yy) %画出概率密度分布图
QDBarCaseCX=x;
QDBarCaseCY=yy;
%画收敛曲线
fz=find(abs(plotGapA)==0);
% fz=fz(1);
figure('Name','互补曲线')
plot(1:fz-1,plotGapA(1:fz-1));
figure('Name','最大不平衡量');
%% 最大不平衡量
% maxDismatchPQ = [0.3123e-10 0.1497e-10 0.7351e-10; 0.6854e-10 0.1973e-10 0.5824e-10];
% bar(maxDismatchPQ);
%% 计算时间
calTime=[70.16 29.68; 68.48 31.661; 65.156 30.08;];
bar(calTime,'stacked');
% figure();
% DeviationFigure(1:33,[CaseAREV;CaseBREV;CaseCREV],[CaseAREA;CaseBREA;CaseCREA],[CaseAREP,CaseBREP,CaseCREP],[CaseAREQ,CaseBREQ,CaseCREQ]);%,[CaseAREA;CaseBREA;CaseCREA],[CaseAREV;CaseBREV;CaseCREV],[CaseAREV;CaseBREV;CaseCREV]);
% VoltBar(voltBarCaseAX,voltBarCaseAY,voltBarCaseBX,voltBarCaseBY,voltBarCaseCX,voltBarCaseCY);
% AngelBar(angelBarCaseAX,angelBarCaseBY,angelBarCaseCX,angelBarCaseAY,angelBarCaseBX,angelBarCaseCY);
% PDBar(PDBarCaseAX,PDBarCaseAY,PDBarCaseBX,PDBarCaseBY,PDBarCaseCX,PDBarCaseCY);
% QDBar(QDBarCaseAX,QDBarCaseAY,QDBarCaseBX,QDBarCaseBY,QDBarCaseCX,QDBarCaseCY);
% MaxErrorFigure()
% MaxBoundErrorFigure();
% DrawLoadProfile();

8
OPF_Init.m
View File

@ -1,9 +1,7 @@
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(Loadi,1)*4+Busnum*2+Busnum+length(Loadi)*2; %约束条件数,放开所有QD
%互补约束数量
RestraintCount=RestraintCount+Busnum+length(Loadi)*2;
RestraintCount=size(Loadi,1)*2+Busnum*2+Busnum; %约束条件数,放开所有QD
t_Bal_volt=Volt(Balance);
Volt=sparse(1*ones(1,Busnum));
Volt(Balance)=t_Bal_volt;
@ -29,8 +27,8 @@ wQD(1:2:end)=0;
%wD(randPDind)=0;%一些负荷不约束
%wD(7)=0;
% wD(11)=0;
PD=0.0*PD0;
PD=0.8*PD0;
%powerFacter=0.98;
%QD=PD*sqrt((1-powerFacter^2)/powerFacter^2);
QD=.0*QD0;
QD=.8*QD0;
end

BIN
PD0.mat
View File

Binary file not shown.

83
PDBar.m
View File

@ -1,83 +0,0 @@
function createfigure(xvector1, yvector1, xvector2, yvector2, xvector3, yvector3)
%CREATEFIGURE(XVECTOR1,YVECTOR1,XVECTOR2,YVECTOR2,XVECTOR3,YVECTOR3)
% XVECTOR1: bar xvector
% YVECTOR1: bar yvector
% XVECTOR2: bar xvector
% YVECTOR2: bar yvector
% XVECTOR3: bar xvector
% YVECTOR3: bar yvector
% Auto-generated by MATLAB on 11-Feb-2015 17:32:07
% Create figure
figure1 = figure('Name','相角直方图','Color',[1 1 1]);
% Create axes
axes1 = axes('Parent',figure1,...
'Position',[0.0920853080568723 0.470619469026549 0.20649289099526 0.372300884955752],...
'FontSize',13);
% Uncomment the following line to preserve the X-limits of the axes
% xlim(axes1,[-4 4]);
% Uncomment the following line to preserve the Y-limits of the axes
% ylim(axes1,[0 6]);
box(axes1,'on');
hold(axes1,'all');
% Create bar
bar(xvector1,yvector1,'FaceColor',[0 0 0],'Parent',axes1);
% Create xlabel
xlabel('相对误差%','FontSize',13);
% Create ylabel
ylabel('数量','FontSize',13);
% Create title
title('情形A','FontSize',13);
% Create axes
axes2 = axes('Parent',figure1,'XTick',[-12 -6 0 6 12],...
'Position',[0.350370561164917 0.466814159292035 0.20649289099526 0.372300884955752],...
'FontSize',13);
% Uncomment the following line to preserve the X-limits of the axes
% xlim(axes2,[-12 12]);
% Uncomment the following line to preserve the Y-limits of the axes
% ylim(axes2,[0 8]);
box(axes2,'on');
hold(axes2,'all');
% Create bar
bar(xvector2,yvector2,'FaceColor',[0 0 0],'Parent',axes2);
% Create xlabel
xlabel('相对误差%','FontSize',13);
% Create title
title('情形B','FontSize',13);
% Create ylabel
ylabel('数量','FontSize',13);
% Create axes
axes3 = axes('Parent',figure1,'XTick',[-12 -6 0 6 12],...
'Position',[0.606286146026517 0.463982300884956 0.20649289099526 0.372300884955752],...
'FontSize',13);
% Uncomment the following line to preserve the X-limits of the axes
% xlim(axes3,[-12 12]);
% Uncomment the following line to preserve the Y-limits of the axes
% ylim(axes3,[0 5]);
box(axes3,'on');
hold(axes3,'all');
% Create bar
bar(xvector3,yvector3,'FaceColor',[0 0 0],'Parent',axes3);
% Create xlabel
xlabel('相对误差%','FontSize',13);
% Create ylabel
ylabel('数量','FontSize',13);
% Create title
title('情形C','FontSize',13);

BIN
QD0.mat
View File

Binary file not shown.

77
QDBar.m
View File

@ -1,77 +0,0 @@
function createfigure(xvector1, yvector1, xvector2, yvector2, xvector3, yvector3)
%CREATEFIGURE(XVECTOR1,YVECTOR1,XVECTOR2,YVECTOR2,XVECTOR3,YVECTOR3)
% XVECTOR1: bar xvector
% YVECTOR1: bar yvector
% XVECTOR2: bar xvector
% YVECTOR2: bar yvector
% XVECTOR3: bar xvector
% YVECTOR3: bar yvector
% Auto-generated by MATLAB on 11-Feb-2015 17:36:00
% Create figure
figure1 = figure('Name','相角直方图','Color',[1 1 1]);
% Create axes
axes1 = axes('Parent',figure1,...
'Position',[0.0920853080568723 0.470619469026549 0.20649289099526 0.372300884955752],...
'FontSize',13);
% Uncomment the following line to preserve the X-limits of the axes
xlim(axes1,[-3 3]);
box(axes1,'on');
hold(axes1,'all');
% Create bar
bar(xvector1,yvector1,'FaceColor',[0 0 0],'Parent',axes1);
% Create xlabel
xlabel('相对误差%','FontSize',13);
% Create ylabel
ylabel('数量','FontSize',13);
% Create title
title('情形A','FontSize',13);
% Create axes
axes2 = axes('Parent',figure1,'XTick',[-12 -6 0 6 12],...
'Position',[0.350370561164917 0.466814159292035 0.20649289099526 0.372300884955752],...
'FontSize',13);
% Uncomment the following line to preserve the X-limits of the axes
xlim(axes2,[-8 8]);
box(axes2,'on');
hold(axes2,'all');
% Create bar
bar(xvector2,yvector2,'FaceColor',[0 0 0],'Parent',axes2);
% Create xlabel
xlabel('相对误差%','FontSize',13);
% Create title
title('情形B','FontSize',13);
% Create ylabel
ylabel('数量','FontSize',13);
% Create axes
axes3 = axes('Parent',figure1,'XTick',[-12 -6 0 6 12],...
'Position',[0.606286146026517 0.463982300884956 0.20649289099526 0.372300884955752],...
'FontSize',13);
% Uncomment the following line to preserve the X-limits of the axes
xlim(axes3,[-9 9]);
box(axes3,'on');
hold(axes3,'all');
% Create bar
bar(xvector3,yvector3,'FaceColor',[0 0 0],'Parent',axes3);
% Create xlabel
xlabel('相对误差%','FontSize',13);
% Create ylabel
ylabel('数量','FontSize',13);
% Create title
title('情形C','FontSize',13);

4
SolveIt.m
View File

@ -1,6 +1,6 @@
function 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)
function XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,RestraintCount,Lx,Balance,PVi,PGi,Busnum,Loadi)
LxComa=FormLxComa(deltF,deltG,deltH,Init_L,Luu,Lul,Init_Z,Init_Y,Lz,Init_U,Init_W,Lw,Lx);
H=-(deltdeltF-ddh-ddgzw);
H=-deltdeltF+ddh;
%t1=diag(Init_L.\Init_Z-Init_U.\Init_W);
t1=diag(Init_Z./Init_L-Init_W./Init_U);
t2=-deltG*( t1 )*deltG';

79
VoltBar.m
View File

@ -1,79 +0,0 @@
function createfigure(xvector1, yvector1, xvector2, yvector2, xvector3, yvector3)
%CREATEFIGURE(XVECTOR1,YVECTOR1,XVECTOR2,YVECTOR2,XVECTOR3,YVECTOR3)
% XVECTOR1: bar xvector
% YVECTOR1: bar yvector
% XVECTOR2: bar xvector
% YVECTOR2: bar yvector
% XVECTOR3: bar xvector
% YVECTOR3: bar yvector
% Auto-generated by MATLAB on 11-Feb-2015 17:21:34
% Create figure
figure1 = figure('Name','电压直方图','Color',[1 1 1]);
% Create axes
axes1 = axes('Parent',figure1,'XTick',[-0.0002 -0.0001 0 0.0001 0.0002],...
'Position',[0.0920853080568723 0.470619469026549 0.20649289099526 0.372300884955752],...
'FontSize',13);
% Uncomment the following line to preserve the Y-limits of the axes
ylim(axes1,[0 7]);
box(axes1,'on');
hold(axes1,'all');
% Create bar
bar(xvector1,yvector1,'FaceColor',[0 0 0],'Parent',axes1);
% Create xlabel
xlabel('误差/p.u','FontSize',13);
% Create ylabel
ylabel('数量','FontSize',13);
% Create title
title('情形A','FontSize',13);
% Create axes
axes2 = axes('Parent',figure1,...
'Position',[0.350370561164917 0.466814159292035 0.20649289099526 0.372300884955752],...
'FontSize',13);
% Uncomment the following line to preserve the X-limits of the axes
xlim(axes2,[-0.004 0.004]);
box(axes2,'on');
hold(axes2,'all');
% Create bar
bar(xvector2,yvector2,'FaceColor',[0 0 0],'Parent',axes2);
% Create xlabel
xlabel('误差/p.u','FontSize',13);
% Create title
title('情形B','FontSize',13);
% Create ylabel
ylabel('数量','FontSize',13);
% Create axes
axes3 = axes('Parent',figure1,...
'Position',[0.606286146026517 0.463982300884956 0.20649289099526 0.372300884955752],...
'FontSize',13);
% Uncomment the following line to preserve the X-limits of the axes
xlim(axes3,[-0.004 0.004]);
% Uncomment the following line to preserve the Y-limits of the axes
ylim(axes3,[0 7]);
box(axes3,'on');
hold(axes3,'all');
% Create bar
bar(xvector3,yvector3,'FaceColor',[0 0 0],'Parent',axes3);
% Create xlabel
xlabel('误差/p.u','FontSize',13);
% Create ylabel
ylabel('数量','FontSize',13);
% Create title
title('情形C','FontSize',13);

31
admmatrix.asv Normal file
View File

@ -0,0 +1,31 @@
function [G,B,GB,Y,r,c,Angle] = admmatrix(Busnum,Linei,Linej,Liner,Linex,Lineb,Transfori...
,Transforj,Transforr,Transforx,Transfork0,Branchi,Branchb)
%**************************************************************************
% 程序功能 : 子函数——形成节点导纳矩阵Y
% 编 者：
% 编制时间：2010.12
%**************************************************************************
%% 支路导纳计算
G = -sparse(Linei,Linej,Liner./(Liner.^2+Linex.^2),Busnum,Busnum) - sparse(Linej,Linei,Liner./(Liner.^2+Linex.^2),Busnum,Busnum);
G = G - sparse(1:Busnum,1:Busnum,sum(G,2)'); % 计算各线路支路电导
B = sparse(Linei,Linej,Linex./(Liner.^2+Linex.^2),Busnum,Busnum)+sparse(Linej,Linei,Linex./(Liner.^2+Linex.^2),Busnum,Busnum);
B = B - sparse(1:Busnum,1:Busnum,sum(B,2)')+sparse(Linei,Linei,Lineb,Busnum,Busnum)+sparse(Linej,Linej,Lineb,Busnum,Busnum);
%% 变压器支路计算
if Transfori>0
mr = Transforr./(Transforr.^2+Transforx.^2); % 计算变压器支路电导
mx = -Transforx./(Transforr.^2+Transforx.^2); % 计算变压器支路电纳
G = G-sparse(Transfori,Transforj,mr./Transfork0,Busnum,Busnum)-sparse(Transforj,Transfori,mr./Transfork0,Busnum,Busnum)...
+sparse(Transfori,Transfori,mr./Transfork0./Transfork0,Busnum,Busnum)+sparse(Transforj,Transforj,mr,Busnum,Busnum);
B = B-sparse(Transfori,Transforj,mx./Transfork0,Busnum,Busnum)-sparse(Transforj,Transfori,mx./Transfork0,Busnum,Busnum)...
+sparse(Transfori,Transfori,mx./Transfork0./Transfork0,Busnum,Busnum)+sparse(Transforj,Transforj,mx,Busnum,Busnum);
end
%% 接地支路计算
if Branchi>0 % 判断有无接地支路
B = B+sparse(Branchi,Branchi,Branchb,Busnum,Busnum);
end
%% 化作极坐标形式
GB = G+B.*1i; %将电导，电纳合并，写成复数形式
Y = abs(GB); %求节点导纳幅值
[r,c] = find(Y);
Angle = angle(GB(GB~=0)); %求节点导纳角度
%Angle=angle(GB);

49
createLoadProfileFigure.m
View File

@ -1,49 +0,0 @@
function createfigure(VertexNormals1, YData1, XData1, Vertices1, Faces1, X1, YMatrix1)
%CREATEFIGURE(VERTEXNORMALS1,YDATA1,XDATA1,VERTICES1,FACES1,X1,YMATRIX1)
% VERTEXNORMALS1: patch vertexnormals
% YDATA1: patch ydata
% XDATA1: patch xdata
% VERTICES1: patch vertices
% FACES1: patch faces
% X1: vector of x data
% YMATRIX1: matrix of y data
% Auto-generated by MATLAB on 11-Feb-2015 22:03:00
% Create figure
figure1 = figure('Color',[1 1 1]);
% Create axes
axes1 = axes('Parent',figure1,...
'XTickLabel',{'6:00','12:00','16:00','24:00','0:00'},...
'XTick',[24 48 72 96 1],...
'Position',[0.13 0.11 0.373795066413662 0.585744680851064],...
'FontName','Times New Roman');
% Uncomment the following line to preserve the X-limits of the axes
xlim(axes1,[1 96]);
hold(axes1,'all');
% Create patch
patch('Parent',axes1,'VertexNormals',VertexNormals1,'YData',YData1,...
'XData',XData1,...
'Vertices',Vertices1,...
'Faces',Faces1,...
'FaceColor',[0.901960784313726 0.901960784313726 0.901960784313726],...
'EdgeColor','none');
% Create multiple lines using matrix input to plot
plot1 = plot(X1,YMatrix1,'Parent',axes1,'LineStyle','--','Color',[0 0 0]);
set(plot1(2),'DisplayName','边界');
set(plot1(3),'DisplayName','负荷曲线','LineStyle','-');
% Create xlabel
xlabel('时刻/时','FontSize',12,'FontName','黑体');
% Create ylabel
ylabel('负荷/kW','FontSize',12,'FontName','黑体');
% Create legend
legend1 = legend(axes1,'show');
set(legend1,...
'Position',[0.155913978494625 0.517021276595745 0.108159392789374 0.168794326241135]);

BIN
figure1.fig
View File

Binary file not shown.

BIN
figure2.fig
View File

Binary file not shown.

19
func_ddg.m
View File

@ -1,17 +1,4 @@
function ddgzw=func_ddg(Busnum,ContrlCount,Loadi,Init_Z,Init_W)
CCVbiZW=Init_Z(length(Loadi)*4+Busnum*2+Busnum+length(Loadi)*2+1:length(Loadi)*4+Busnum*2+Busnum+length(Loadi)*2+Busnum)+Init_W(length(Loadi)*4+Busnum*2+Busnum+length(Loadi)*2+1:length(Loadi)*4+Busnum*2+Busnum+length(Loadi)*2+Busnum);
CCPDbiZW=Init_Z(length(Loadi)*4+Busnum*2+Busnum+length(Loadi)*2+Busnum+1:length(Loadi)*4+Busnum*2+Busnum+length(Loadi)*2+Busnum+length(Loadi))+Init_W(length(Loadi)*4+Busnum*2+Busnum+length(Loadi)*2+Busnum+1:length(Loadi)*4+Busnum*2+Busnum+length(Loadi)*2+Busnum+length(Loadi));
CCQDbiZW=Init_Z(length(Loadi)*4+Busnum*2+Busnum+length(Loadi)*2+Busnum+length(Loadi)+1:end)+Init_W(length(Loadi)*4+Busnum*2+Busnum+length(Loadi)*2+Busnum+length(Loadi)+1:end);
% ddgzw=[
% sparse(length(Loadi)*2+Busnum*2,ContrlCount);
% sparse(Busnum,length(Loadi)*2+Busnum*2),sparse(-2*eye(Busnum))*diag(CCVbiZW),sparse(Busnum,length(Loadi)*2);
% sparse(length(Loadi),length(Loadi)*2+Busnum*2+Busnum),sparse(-2*eye(length(Loadi)))*diag(CCPDbiZW),sparse(length(Loadi),length(Loadi));
% sparse(length(Loadi),length(Loadi)*2+Busnum*2+Busnum+length(Loadi)),sparse(-2*eye(length(Loadi)))*diag(CCQDbiZW);
% ];
ddgzw=[
sparse(length(Loadi)*2+Busnum*2,ContrlCount);
sparse(Busnum,length(Loadi)*2+Busnum*2),sparse(-2*eye(Busnum))*diag(CCVbiZW),sparse(Busnum,length(Loadi)*2);
sparse(length(Loadi),length(Loadi)*2+Busnum*2+Busnum),sparse(-2*eye(length(Loadi)))*diag(CCPDbiZW),sparse(length(Loadi),length(Loadi));
sparse(length(Loadi),length(Loadi)*2+Busnum*2+Busnum+length(Loadi)),sparse(-2*eye(length(Loadi)))*diag(CCQDbiZW);
];
function ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi,PD,QD)
ddg=0;
end

4
func_ddh.m
View File

@ -51,7 +51,7 @@ t=[ddPdVdV+ddQdVdV,ddPdVdT+ddQdVdT;
sizeLoadi=size(Loadi,1)*2;
ddh=[
sparse(sizeLoadi,ContrlCount);
sparse(2*Busnum,sizeLoadi),-t,sparse(2*Busnum,Busnum+length(Loadi)*2);
sparse(Busnum+length(Loadi)*2,ContrlCount);
sparse(2*Busnum,sizeLoadi),-t,sparse(2*Busnum,Busnum);
sparse(Busnum,ContrlCount);
];
end

4
func_deltF.m
View File

@ -1,10 +1,10 @@
function deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,PG0,QG0,PD0,PD,QD,QD0,Busnum,Loadi)
t3=2*wPD.*(PD(Loadi)-PD0(Loadi));
t4=2*wQD.*(QD(Loadi)-QD0(Loadi));
deltF=[sparse(length(Loadi),1);
sparse(length(Loadi),1);
sparse(2*Busnum,1);
sparse(ones(Busnum,1));
sparse(ones(length(Loadi),1));
sparse(ones(length(Loadi),1))
];
end

34
func_deltG.asv Normal file
View File

@ -0,0 +1,34 @@
function deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD)
sizeLoadi=size(Loadi,1);
dg3_dPD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
dg4_dPD=sparse(size(Loadi,1),length(Loadi));
dg5_dPD=sparse(size(Loadi,1),Busnum);
dg6_dPD=dg5_dPD;
dg7_dPD=sparse(sizeLoadi,Busnum);
%%
dg3_dQD=sparse(length(Loadi),length(Loadi));
dg4_dQD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
dg5_dQD=sparse(size(Loadi,1),Busnum);
dg6_dQD=dg5_dQD;
dg7_dQD=sparse(sizeLoadi,Busnum);
%%
dg3_dx=sparse(2*Busnum,sizeLoadi);
dg4_dx=sparse(2*Busnum,length(Loadi));
dg5_dx=[sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum);
sparse(Busnum,Busnum);
];
dg6_dx=dg5_dx;
dg7_dx=sparse(2*Busnum,Busnum);
%%
dg3_dvbi=sparse(Busnum,sizeLoadi);
dg4_dvbi=sparse(Busnum,length(Loadi));
dg5_dvbi=sparse(Busnum,Busnum);
dg6_dvbi=sparse(Busnum,Busnum);
dg6_dvbi=sparse(Busnum,Busnum);
%%
deltG=[dg3_dPD,dg4_dPD,dg5_dPD,dg6_dPD;
dg3_dQD,dg4_dQD,dg5_dQD,dg6_dQD;
dg3_dx,dg4_dx,dg5_dx,dg6_dx;
dg3_dvbi,dg4_dvbi,dg5_dvbi,dg6_dvbi;
];
end

89
func_deltG.m
View File

@ -1,107 +1,34 @@
function deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD,Vbi,PDbi,QDbi)
function deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD)
sizeLoadi=size(Loadi,1);
%dg3 PD-M*b-t-PD0
dg3_dPD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
%dg32 PD+M*B+t+PD0
dg32_dPD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
dg4_dPD=sparse(size(Loadi,1),length(Loadi));
dg42_dPD=sparse(size(Loadi,1),length(Loadi));
dg5_dPD=sparse(size(Loadi,1),Busnum);
dg6_dPD=sparse(size(Loadi,1),Busnum);
dg6_dPD=dg5_dPD;
dg7_dPD=sparse(sizeLoadi,Busnum);
dgPDbi_dPD=sparse(sizeLoadi,sizeLoadi);
dgQDbi_dPD=sparse(sizeLoadi,sizeLoadi);
%互补约束
dgCCVbi_dPD=sparse(length(Loadi),Busnum);
dgCCPDbi_dPD=sparse(length(Loadi),length(Loadi));
dgCCQDbi_dPD=sparse(length(Loadi),length(Loadi));
%%
dg3_dQD=sparse(length(Loadi),length(Loadi));
dg32_dQD=sparse(length(Loadi),length(Loadi));
dg4_dQD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
dg42_dQD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
dg5_dQD=sparse(size(Loadi,1),Busnum);
dg6_dQD=sparse(size(Loadi,1),Busnum);
dg6_dQD=dg5_dQD;
dg7_dQD=sparse(sizeLoadi,Busnum);
dgPDbi_dQD=sparse(sizeLoadi,sizeLoadi);
dgQDbi_dQD=sparse(sizeLoadi,sizeLoadi);
%互补约束
dgCCVbi_dQD=sparse(length(Loadi),Busnum);
dgCCPDbi_dQD=sparse(length(Loadi),length(Loadi));
dgCCQDbi_dQD=sparse(length(Loadi),length(Loadi));
%%
dg3_dx=sparse(2*Busnum,sizeLoadi);
dg32_dx=sparse(2*Busnum,sizeLoadi);
dg4_dx=sparse(2*Busnum,length(Loadi));
dg42_dx=sparse(2*Busnum,length(Loadi));
dg5_dx=[sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum);
sparse(Busnum,Busnum);
];
dg6_dx=[sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum);
sparse(Busnum,Busnum);
];
dg6_dx=dg5_dx;
dg7_dx=sparse(2*Busnum,Busnum);
dgPDbi_dx=sparse(2*Busnum,sizeLoadi);
dgQDbi_dx=sparse(2*Busnum,sizeLoadi);
%互补约束
dgCCVbi_dx=sparse(Busnum*2,Busnum);
dgCCPDbi_dx=sparse(Busnum*2,length(Loadi));
dgCCQDbi_dx=sparse(Busnum*2,length(Loadi));
%%
dg3_dvbi=sparse(Busnum,sizeLoadi);
dg32_dvbi=sparse(Busnum,sizeLoadi);
dg4_dvbi=sparse(Busnum,length(Loadi));
dg42_dvbi=sparse(Busnum,length(Loadi));
dg5_dvbi=sparse(-eye(Busnum,Busnum));
dg6_dvbi=sparse(eye(Busnum,Busnum));
dg7_dvbi=sparse(eye(Busnum,Busnum));
dgPDbi_dvbi=sparse(Busnum,sizeLoadi);
dgQDbi_dvbi=sparse(Busnum,sizeLoadi);
%互补约束
dgCCVbi_dvbi=sparse(diag(0.1-2*Vbi));%互补约束为(1-Vbi).*Vbi
dgCCPDbi_dvbi=sparse(Busnum,length(Loadi));
dgCCQDbi_dvbi=sparse(Busnum,length(Loadi));
%%
dg3_dPDbi=sparse(-eye(sizeLoadi));
dg32_dPDbi=sparse(eye(sizeLoadi));
dg4_dPDbi=sparse(length(Loadi),length(Loadi));
dg42_dPDbi=sparse(length(Loadi),length(Loadi));
dg5_dPDbi=sparse(length(Loadi),Busnum);
dg6_dPDbi=sparse(length(Loadi),Busnum);
dg7_dPDbi=sparse(length(Loadi),Busnum);
dgPDbi_dPDbi=sparse(eye(length(Loadi)));
dgQDbi_dPDbi=sparse(zeros(length(Loadi)));
%互补约束
dgCCVbi_dPDbi=sparse(length(Loadi),Busnum);
dgCCPDbi_dPDbi=sparse(diag(0.1-2*PDbi));
dgCCQDbi_dPDbi=sparse(length(Loadi),length(Loadi));
%%
dg3_dQDbi=sparse(length(Loadi),length(Loadi));
dg32_dQDbi=sparse(length(Loadi),length(Loadi));
dg4_dQDbi=sparse(-eye(sizeLoadi));
dg42_dQDbi=sparse(eye(sizeLoadi));
dg5_dQDbi=sparse(length(Loadi),Busnum);
dg6_dQDbi=sparse(length(Loadi),Busnum);
dg7_dQDbi=sparse(length(Loadi),Busnum);
dgPDbi_dQDbi=sparse(zeros(length(Loadi)));
dgQDbi_dQDbi=sparse(eye(length(Loadi)));
%互补约束
dgCCVbi_dQDbi=sparse(length(Loadi),Busnum);
dgCCPDbi_dQDbi=sparse(length(Loadi),length(Loadi));
dgCCQDbi_dQDbi=sparse(diag(0.1-2*QDbi));
%%
% deltG=[dg3_dPD,dg32_dPD,dg4_dPD,dg42_dPD,dg5_dPD,dg6_dPD,dg7_dPD,dgPDbi_dPD,dgQDbi_dPD,dgCCVbi_dPD,dgCCPDbi_dPD,dgCCQDbi_dPD;
% dg3_dQD,dg32_dQD,dg4_dQD,dg42_dQD,dg5_dQD,dg6_dQD,dg7_dQD,dgPDbi_dQD,dgQDbi_dQD,dgCCVbi_dQD,dgCCPDbi_dQD,dgCCQDbi_dQD;
% dg3_dx,dg32_dx,dg4_dx,dg42_dx,dg5_dx,dg6_dx,dg7_dx,dgPDbi_dx,dgQDbi_dx,dgCCVbi_dx,dgCCPDbi_dx,dgCCQDbi_dx;
% dg3_dvbi,dg32_dvbi,dg4_dvbi,dg42_dvbi,dg5_dvbi,dg6_dvbi,dg7_dvbi,dgPDbi_dvbi,dgQDbi_dvbi,dgCCVbi_dvbi,dgCCPDbi_dvbi,dgCCQDbi_dvbi;
% dg3_dPDbi,dg32_dPDbi,dg4_dPDbi,dg42_dPDbi,dg5_dPDbi,dg6_dPDbi,dg7_dPDbi,dgPDbi_dPDbi,dgQDbi_dPDbi,dgCCVbi_dPDbi,dgCCPDbi_dPDbi,dgCCQDbi_dPDbi;
% dg3_dQDbi,dg32_dQDbi,dg4_dQDbi,dg42_dQDbi,dg5_dQDbi,dg6_dQDbi,dg7_dQDbi,dgPDbi_dQDbi,dgQDbi_dQDbi,dgCCVbi_dQDbi,dgCCPDbi_dQDbi,dgCCQDbi_dQDbi;
% ];
deltG=[dg3_dPD,dg32_dPD,dg4_dPD,dg42_dPD,dg5_dPD,dg6_dPD,dg7_dPD,dgPDbi_dPD,dgQDbi_dPD,dgCCVbi_dPD,dgCCPDbi_dPD,dgCCQDbi_dPD;
dg3_dQD,dg32_dQD,dg4_dQD,dg42_dQD,dg5_dQD,dg6_dQD,dg7_dQD,dgPDbi_dQD,dgQDbi_dQD,dgCCVbi_dQD,dgCCPDbi_dQD,dgCCQDbi_dQD;
dg3_dx,dg32_dx,dg4_dx,dg42_dx,dg5_dx,dg6_dx,dg7_dx,dgPDbi_dx,dgQDbi_dx,dgCCVbi_dx,dgCCPDbi_dx,dgCCQDbi_dx;
dg3_dvbi,dg32_dvbi,dg4_dvbi,dg42_dvbi,dg5_dvbi,dg6_dvbi,dg7_dvbi,dgPDbi_dvbi,dgQDbi_dvbi,dgCCVbi_dvbi,dgCCPDbi_dvbi,dgCCQDbi_dvbi;
dg3_dPDbi,dg32_dPDbi,dg4_dPDbi,dg42_dPDbi,dg5_dPDbi,dg6_dPDbi,dg7_dPDbi,dgPDbi_dPDbi,dgQDbi_dPDbi,dgCCVbi_dPDbi,dgCCPDbi_dPDbi,dgCCQDbi_dPDbi;
dg3_dQDbi,dg32_dQDbi,dg4_dQDbi,dg42_dQDbi,dg5_dQDbi,dg6_dQDbi,dg7_dQDbi,dgPDbi_dQDbi,dgQDbi_dQDbi,dgCCVbi_dQDbi,dgCCPDbi_dQDbi,dgCCQDbi_dQDbi;
deltG=[dg3_dPD,dg4_dPD,dg5_dPD,dg6_dPD,dg7_dPD;
dg3_dQD,dg4_dQD,dg5_dQD,dg6_dQD,dg7_dQD;
dg3_dx,dg4_dx,dg5_dx,dg6_dx,dg7_dx;
dg3_dvbi,dg4_dvbi,dg5_dvbi,dg6_dvbi,dg7_dvbi;
];
end

4
func_deltH.m
View File

@ -3,8 +3,6 @@ dH_dPD=[sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum
dH_dQD=[sparse(size(Loadi,1),Busnum) sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum)];
dH_dx = jacobian_M(Busnum,Volt,Y,Angle,UAngel,r,c); %形成雅克比矩阵
dH_dvbi=sparse(Busnum,2*Busnum);
dH_dPDbi=sparse(length(Loadi),2*Busnum);
dH_dQDbi=sparse(length(Loadi),2*Busnum);
deltH=[dH_dPD;dH_dQD;dH_dx';dH_dvbi;dH_dPDbi;dH_dQDbi];
deltH=[dH_dPD;dH_dQD;dH_dx';dH_dvbi];
end

16
graph/badDataGraph.m
View File

@ -1,16 +0,0 @@
x=1:33;
y=1:33;
x=[x,15,33];
y=[y,17,23];
scatter(x,y);
xlabel('存在坏数据的节点');
ylabel('辨识出坏数据的点');
hold on;
for I=1:33
% plot([1,1,1,1]-1+I,[[1:3]-2+I,I]);
plot([I,I],[0,I]);
end
% plot([1,1,1,1]-1+15,[1:3]-2+17);
% plot([1,1,1,1]-1+33,[1,2,12]-2+23);
plot([15,15],[0,17]);
plot([33,33],[0,23]);

BIN
histA.fig
View File

Binary file not shown.

BIN
histP.fig
View File

Binary file not shown.

BIN
histQ.fig
View File

Binary file not shown.

BIN
histV.fig
View File

Binary file not shown.

258
hs_err_pid4476.log
View File

@ -1,258 +0,0 @@
#
# An unexpected error has been detected by Java Runtime Environment:
#
# EXCEPTION_ACCESS_VIOLATION (0xc0000005) at pc=0x77cb80dc, pid=4476, tid=6716
#
# Java VM: Java HotSpot(TM) Client VM (11.2-b01 mixed mode windows-x86)
# Problematic frame:
# C [ntdll.dll+0x580dc]
#
# If you would like to submit a bug report, please visit:
# http://java.sun.com/webapps/bugreport/crash.jsp
# The crash happened outside the Java Virtual Machine in native code.
# See problematic frame for where to report the bug.
#
--------------- T H R E A D ---------------
Current thread (0x09b9c800): JavaThread "main" [_thread_in_native, id=6716, stack(0x00430000,0x00c30000)]
siginfo: ExceptionCode=0xc0000005, reading address 0xffffffff
Registers:
EAX=0x00000000, EBX=0x00000000, ECX=0x00000008, EDX=0xfffffff8
ESP=0x00c29b0c, EBP=0x00c29b14, ESI=0x09b90000, EDI=0x79cd34e4
EIP=0x77cb80dc, EFLAGS=0x00010246
Top of Stack: (sp=0x00c29b0c)
0x00c29b0c: 00c29c14 79d135a4 00c29b60 7bfc20d6
0x00c29b1c: 09b90000 00000000 00000000 79cd34e4
0x00c29b2c: 79d135a4 00c29c14 00000000 056b4b88
0x00c29b3c: 000007fd 00c29bc0 00c29bc4 00c29b28
0x00c29b4c: 00c2954c 00c29b8c 7bfc240d 7bffa368
0x00c29b5c: ffffffff 00c29b9c 7bfcc0b4 00000000
0x00c29b6c: 79cd34e4 00000003 00c29c14 002b0023
0x00c29b7c: 0053002b 002b002b 00c29b6c 00c2954c
Instructions: (pc=0x77cb80dc)
0x77cb80cc: 02 00 8b 45 10 a8 07 0f 85 41 f7 02 00 8d 50 f8
0x77cb80dc: 80 7a 07 05 0f 84 1d f7 02 00 f6 42 07 3f 0f 84
Stack: [0x00430000,0x00c30000], sp=0x00c29b0c, free space=8166k
Native frames: (J=compiled Java code, j=interpreted, Vv=VM code, C=native code)
C [ntdll.dll+0x580dc]
C [MSVCR71.dll+0x20d6]
C [MSVCR71.dll+0xc0b4]
V [jvm.dll+0x1c8b74]
Java frames: (J=compiled Java code, j=interpreted, Vv=VM code)
j com.mathworks.jmi.NativeMatlab.SendMatlabMessage(Ljava/lang/Object;)Ljava/lang/Object;+0
j com.mathworks.jmi.NativeMatlab.sendMatlabMessage(Ljava/lang/Object;)Ljava/lang/Object;+22
j com.mathworks.jmi.MatlabLooper.sendMatlabMessage(Lcom/mathworks/services/message/MWMessage;)Ljava/lang/Object;+20
j com.mathworks.jmi.Matlab.mtFevalConsoleOutput(Ljava/lang/String;[Ljava/lang/Object;I)Ljava/lang/Object;+58
j com.mathworks.mde.desk.MLDesktop$9.run()V+14
j com.mathworks.jmi.NativeMatlab.dispatchMTRequests(Z)V+50
v ~StubRoutines::call_stub
--------------- P R O C E S S ---------------
Java Threads: ( => current thread )
0x22e17400 JavaThread "Inactive RequestProcessor thread [Was:TimedSoftReference/org.openide.util.TimedSoftReference]" daemon [_thread_blocked, id=1756, stack(0x05f20000,0x05fa0000)]
0x25229400 JavaThread "Thread-1175" [_thread_in_native, id=6188, stack(0x2fd00000,0x30500000)]
0x250c0000 JavaThread "pool-2-thread-1" [_thread_blocked, id=8112, stack(0x05bc0000,0x05c40000)]
0x231b5800 JavaThread "Timer-120" [_thread_blocked, id=8128, stack(0x04710000,0x04790000)]
0x22ff0800 JavaThread "Active Reference Queue Daemon" daemon [_thread_blocked, id=7892, stack(0x1eee0000,0x1ef60000)]
0x22ff3c00 JavaThread "Timer-4" daemon [_thread_blocked, id=5616, stack(0x1ebc0000,0x1ec40000)]
0x22ff3800 JavaThread "Timer-3" daemon [_thread_blocked, id=4204, stack(0x1eaf0000,0x1eb70000)]
0x22ff2c00 JavaThread "Prefs Updater" daemon [_thread_blocked, id=7656, stack(0x24300000,0x24380000)]
0x22ff3000 JavaThread "Timer-1" [_thread_blocked, id=2792, stack(0x23fe0000,0x24060000)]
0x22c6c400 JavaThread "TimerQueue" daemon [_thread_blocked, id=6316, stack(0x23df0000,0x23e70000)]
0x22c66c00 JavaThread "AWT-EventQueue-0" [_thread_blocked, id=2764, stack(0x23d20000,0x23da0000)]
0x22b78000 JavaThread "AWT-Windows" daemon [_thread_in_native, id=7904, stack(0x23300000,0x23380000)]
0x22bffc00 JavaThread "AWT-Shutdown" [_thread_blocked, id=7840, stack(0x23280000,0x23300000)]
0x22bff800 JavaThread "Java2D Disposer" daemon [_thread_blocked, id=8164, stack(0x23200000,0x23280000)]
0x1c7c9800 JavaThread "Timer-0" [_thread_blocked, id=7820, stack(0x1c990000,0x1ca10000)]
0x1c793800 JavaThread "JMI Unnamed Thread" [_thread_in_native, id=7988, stack(0x061c0000,0x069c0000)]
0x09e2c800 JavaThread "Low Memory Detector" daemon [_thread_blocked, id=7900, stack(0x1c680000,0x1c700000)]
0x09e2b800 JavaThread "CompilerThread0" daemon [_thread_blocked, id=8152, stack(0x1f950000,0x20150000)]
0x09e0d400 JavaThread "Attach Listener" daemon [_thread_blocked, id=7876, stack(0x1c5c0000,0x1c640000)]
0x09df8400 JavaThread "Finalizer" daemon [_thread_blocked, id=8136, stack(0x1e3f0000,0x1e470000)]
0x09df6c00 JavaThread "Reference Handler" daemon [_thread_blocked, id=7848, stack(0x0ff80000,0x10000000)]
=>0x09b9c800 JavaThread "main" [_thread_in_native, id=6716, stack(0x00430000,0x00c30000)]
Other Threads:
0x09df1000 VMThread [stack: 0x1db90000,0x1e390000] [id=3524]
0x09e32400 WatcherThread [stack: 0x21e00000,0x22600000] [id=5768]
VM state:not at safepoint (normal execution)
VM Mutex/Monitor currently owned by a thread: None
Heap
def new generation total 29504K, used 23735K [0x10390000, 0x12390000, 0x12390000)
eden space 26240K, 83% used [0x10390000, 0x118fcfb8, 0x11d30000)
from space 3264K, 55% used [0x12060000, 0x12220f70, 0x12390000)
to space 3264K, 0% used [0x11d30000, 0x11d30000, 0x12060000)
tenured generation total 98304K, used 69305K [0x12390000, 0x18390000, 0x18390000)
the space 98304K, 70% used [0x12390000, 0x1673e6f8, 0x1673e800, 0x18390000)
compacting perm gen total 39680K, used 39594K [0x18390000, 0x1aa50000, 0x1c390000)
the space 39680K, 99% used [0x18390000, 0x1aa3a830, 0x1aa3aa00, 0x1aa50000)
No shared spaces configured.
Dynamic libraries:
0x00400000 - 0x00425000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\MATLAB.exe
0x77c60000 - 0x77db7000 C:\Windows\SYSTEM32\ntdll.dll
0x75bf0000 - 0x75d20000 C:\Windows\SYSTEM32\KERNEL32.DLL
0x76360000 - 0x76406000 C:\Windows\SYSTEM32\KERNELBASE.dll
0x74720000 - 0x747c7000 C:\Windows\system32\apphelp.dll
0x7b940000 - 0x7b9ec000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libut.dll
0x78770000 - 0x787ee000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmwfl.dll
0x7ba50000 - 0x7bc63000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmwservices.dll
0x7a0e0000 - 0x7a144000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\mcr.dll
0x00140000 - 0x001c7000 C:\Windows\WinSxS\x86_microsoft.vc80.crt_1fc8b3b9a1e18e3b_8.0.50727.6910_none_d089c358442de345\MSVCP80.dll
0x001d0000 - 0x0026b000 C:\Windows\WinSxS\x86_microsoft.vc80.crt_1fc8b3b9a1e18e3b_8.0.50727.6910_none_d089c358442de345\MSVCR80.dll
0x76900000 - 0x76912000 C:\Windows\SYSTEM32\imagehlp.dll
0x75d20000 - 0x75d26000 C:\Windows\SYSTEM32\PSAPI.DLL
0x7bfa0000 - 0x7bfbf000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\LIBEXPAT.dll
0x7bd50000 - 0x7be38000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\icuuc40.dll
0x7b4c0000 - 0x7b4cc000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\icuio40.dll
0x76490000 - 0x765a6000 C:\Windows\SYSTEM32\USER32.dll
0x76940000 - 0x769ee000 C:\Windows\SYSTEM32\ADVAPI32.dll
0x7d160000 - 0x7d16e000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\boost_date_time-vc80-mt-1_36.dll
0x7d330000 - 0x7d346000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\boost_filesystem-vc80-mt-1_36.dll
0x7d050000 - 0x7d060000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\boost_signals-vc80-mt-1_36.dll
0x7c500000 - 0x7c507000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\boost_system-vc80-mt-1_36.dll
0x7b4a0000 - 0x7b4ac000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\boost_thread-vc80-mt-1_36.dll
0x7b9f0000 - 0x7ba4d000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmx.dll
0x7bc70000 - 0x7bd47000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmwmathutil.dll
0x79e30000 - 0x79e6a000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\mpath.dll
0x7c0a0000 - 0x7c166000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\mlutil.dll
0x78720000 - 0x7874e000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\tbb.dll
0x73d00000 - 0x73d87000 C:\Windows\WinSxS\x86_microsoft.windows.common-controls_6595b64144ccf1df_5.82.9200.16658_none_bf1359a245f1cd12\COMCTL32.dll
0x76870000 - 0x768f9000 C:\Windows\SYSTEM32\comdlg32.dll
0x75030000 - 0x75042000 C:\Windows\SYSTEM32\NETAPI32.dll
0x753a0000 - 0x753f0000 C:\Windows\SYSTEM32\WS2_32.dll
0x76ad0000 - 0x77b96000 C:\Windows\SYSTEM32\SHELL32.dll
0x7e890000 - 0x7e8a9000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmwi18n.dll
0x7b4d0000 - 0x7b539000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\iqm.dll
0x7e850000 - 0x7e871000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmwbridge.dll
0x7b550000 - 0x7b561000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmex.dll
0x79fb0000 - 0x7a017000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\m_dispatcher.dll
0x787f0000 - 0x7894e000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\mcos.dll
0x7be40000 - 0x7bf3f000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmwgui.dll
0x78cf0000 - 0x7914f000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\hg.dll
0x785c0000 - 0x78628000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\hgutils.dll
0x7a6c0000 - 0x7ab7b000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\m_interpreter.dll
0x79d40000 - 0x79d96000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\jmi.dll
0x7a150000 - 0x7a1cc000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\udd_mi.dll
0x7b210000 - 0x7b492000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\uiw.dll
0x78c20000 - 0x78c37000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\mwoles05.DLL
0x79e70000 - 0x79eb9000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\comcli.dll
0x7b600000 - 0x7b60b000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\mlautoregister.dll
0x760b0000 - 0x761ad000 C:\Windows\SYSTEM32\GDI32.dll
0x75e80000 - 0x75f99000 C:\Windows\SYSTEM32\ole32.dll
0x767b0000 - 0x76861000 C:\Windows\SYSTEM32\msvcrt.dll
0x7b4b0000 - 0x7b4b4000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\icudt40.dll
0x7b670000 - 0x7b768000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\icuin40.dll
0x75400000 - 0x75434000 C:\Windows\SYSTEM32\sechost.dll
0x75d30000 - 0x75ddc000 C:\Windows\SYSTEM32\RPCRT4.dll
0x78c90000 - 0x78ca3000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\zlib1.dll
0x7cfb0000 - 0x7d00e000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\xmlcore.dll
0x75800000 - 0x75840000 C:\Windows\SYSTEM32\SHLWAPI.dll
0x75020000 - 0x7502b000 C:\Windows\SYSTEM32\netutils.dll
0x75000000 - 0x7501c000 C:\Windows\SYSTEM32\srvcli.dll
0x74ff0000 - 0x75000000 C:\Windows\SYSTEM32\wkscli.dll
0x75de0000 - 0x75de8000 C:\Windows\SYSTEM32\NSI.dll
0x756c0000 - 0x757f6000 C:\Windows\SYSTEM32\combase.dll
0x7b180000 - 0x7b209000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\udd.dll
0x79720000 - 0x798e2000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\xerces-c_2_7.dll
0x7d280000 - 0x7d31d000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\boost_regex-vc80-mt-1_36.dll
0x78b60000 - 0x78b92000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmat.dll
0x79a40000 - 0x79a70000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmwhardcopy.dll
0x79da0000 - 0x79de3000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libuij.dll
0x78630000 - 0x7871f000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\hgdatatypes.dll
0x78ba0000 - 0x78bd4000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmwlapack.dll
0x00cb0000 - 0x00dbf000 C:\Windows\WinSxS\x86_microsoft.vc80.mfc_1fc8b3b9a1e18e3b_8.0.50727.6195_none_cbf5e994470a1a8f\MFC80.DLL
0x73470000 - 0x734d0000 C:\Windows\SYSTEM32\WINSPOOL.DRV
0x76720000 - 0x767ab000 C:\Windows\SYSTEM32\OLEAUT32.dll
0x79ec0000 - 0x79f1b000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\profiler.dll
0x7b5d0000 - 0x7b5f2000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmwmathrng.dll
0x78c00000 - 0x78c12000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\m_pcodeio.dll
0x79df0000 - 0x79e27000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\m_ir.dll
0x7a1d0000 - 0x7a6bd000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\m_parser.dll
0x78be0000 - 0x78bfa000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\m_pcodegen.dll
0x7e810000 - 0x7e844000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmwspmatrix.dll
0x7b660000 - 0x7b669000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\uinone.dll
0x00f20000 - 0x00f3b000 C:\Windows\WinSxS\x86_microsoft.vc80.atl_1fc8b3b9a1e18e3b_8.0.50727.6195_none_d1cb102c435421de\ATL80.DLL
0x75250000 - 0x7526c000 C:\Windows\SYSTEM32\SspiCli.dll
0x7b770000 - 0x7b939000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libhdf5.dll
0x7b580000 - 0x7b58e000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmwbinder.dll
0x7b540000 - 0x7b54e000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\ir_xfmr.dll
0x7b610000 - 0x7b61a000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\mtok.dll
0x75240000 - 0x75249000 C:\Windows\SYSTEM32\CRYPTBASE.dll
0x75140000 - 0x751b3000 C:\Windows\SYSTEM32\SHCORE.DLL
0x751e0000 - 0x75231000 C:\Windows\SYSTEM32\bcryptPrimitives.dll
0x76920000 - 0x76940000 C:\Windows\system32\IMM32.DLL
0x769f0000 - 0x76acc000 C:\Windows\SYSTEM32\MSCTF.dll
0x03a20000 - 0x03a2a000 C:\Windows\WinSxS\x86_microsoft.vc80.mfcloc_1fc8b3b9a1e18e3b_8.0.50727.6195_none_03ce2c72205943d3\MFC80CHS.DLL
0x76470000 - 0x76481000 C:\Windows\SYSTEM32\profapi.dll
0x73ac0000 - 0x73b48000 C:\Windows\system32\uxtheme.dll
0x70cb0000 - 0x70e70000 C:\Program Files (x86)\360\360safe\safemon\safemon.dll
0x73dc0000 - 0x73fb7000 C:\Windows\WinSxS\x86_microsoft.windows.common-controls_6595b64144ccf1df_6.0.9200.16384_none_893961408605e985\COMCTL32.dll
0x751d0000 - 0x751d8000 C:\Windows\SYSTEM32\VERSION.dll
0x72b50000 - 0x72b90000 C:\Program Files (x86)\TeamViewer\Version9\tv_w32.dll
0x761b0000 - 0x7635f000 C:\Windows\SYSTEM32\SETUPAPI.dll
0x75460000 - 0x754a6000 C:\Windows\SYSTEM32\CFGMGR32.dll
0x75440000 - 0x7545e000 C:\Windows\SYSTEM32\DEVOBJ.dll
0x75fa0000 - 0x76014000 C:\Windows\SYSTEM32\clbcatq.dll
0x747e0000 - 0x748f9000 C:\Windows\system32\propsys.dll
0x74f40000 - 0x74f62000 C:\Windows\SYSTEM32\iphlpapi.dll
0x74f30000 - 0x74f38000 C:\Windows\SYSTEM32\WINNSI.DLL
0x74cd0000 - 0x74d45000 C:\Windows\SYSTEM32\DNSAPI.dll
0x74b90000 - 0x74ba0000 C:\Windows\SYSTEM32\dhcpcsvc6.DLL
0x74df0000 - 0x74e02000 C:\Windows\SYSTEM32\dhcpcsvc.DLL
0x042e0000 - 0x042e3000 C:\Windows\SYSTEM32\icmp.Dll
0x79ae0000 - 0x79d36000 C:\Program Files (x86)\MATLAB\R2010a\sys\java\jre\win32\jre\bin\client\jvm.dll
0x74e40000 - 0x74e61000 C:\Windows\SYSTEM32\WINMM.dll
0x7bfc0000 - 0x7c016000 C:\Program Files (x86)\MATLAB\R2010a\bin\win32\MSVCR71.dll
0x74e10000 - 0x74e3a000 C:\Windows\SYSTEM32\WINMMBASE.dll
0x749f0000 - 0x74b19000 C:\Windows\SYSTEM32\dbghelp.dll
0x79f70000 - 0x79f78000 C:\Program Files (x86)\MATLAB\R2010a\sys\java\jre\win32\jre\bin\hpi.dll
0x72930000 - 0x729a1000 C:\Windows\SYSTEM32\SogouTSF.ime
0x735b0000 - 0x735b6000 C:\Windows\SYSTEM32\MSIMG32.dll
0x73780000 - 0x737a6000 C:\Windows\SYSTEM32\ntmarta.dll
0x79f90000 - 0x79f9c000 C:\Program Files (x86)\MATLAB\R2010a\sys\java\jre\win32\jre\bin\verify.dll
0x78cb0000 - 0x78ccf000 C:\Program Files (x86)\MATLAB\R2010a\sys\java\jre\win32\jre\bin\java.dll
0x79fa0000 - 0x79faf000 C:\Program Files (x86)\MATLAB\R2010a\sys\java\jre\win32\jre\bin\zip.dll
0x10000000 - 0x1038d000 C:\Windows\system32\SogouPy.ime
0x73340000 - 0x73390000 C:\Windows\SYSTEM32\OLEACC.dll
0x72900000 - 0x72923000 C:\Program Files (x86)\SogouInput\Components\SgAppender\1.0.0.207\SgAppender_Dll.dll
0x0fc70000 - 0x0fd2e000 C:\Program Files (x86)\SogouInput\7.1.0.2005\Resource.dll
VM Arguments:
jvm_args: -Xss512k -XX:PermSize=32m -Xms64m -XX:NewRatio=3 -XX:MaxPermSize=64m -Xmx128m -XX:MaxDirectMemorySize=1200000000 -Dsun.java2d.noddraw=true -Dsun.awt.nopixfmt=true -Xshare:off -Xrs -Djava.library.path=C:\Program Files (x86)\MATLAB\R2010a\bin\win32 vfprintf abort
java_command: <unknown>
Launcher Type: generic
Environment Variables:
PATH=C:\Program Files (x86)\MATLAB\R2010a\sys\java\jre\win32\jre\bin;C:\Program Files (x86)\MATLAB\R2010a\sys\webrenderer\windows\corecomponents;C:\Program Files (x86)\MATLAB\R2010a\sys\webrenderer\windows;C:\Program Files\Java\jre7\bin;C:\Program Files (x86)\AMD APP\bin\x86_64;C:\Program Files (x86)\AMD APP\bin\x86;C:\Program Files (x86)\Git\bin;c:/go/bin;C:\Python34\;C:\Python34\Scripts;C:\Program Files (x86)\Common Files\Intel\Shared Libraries\redist\ia32\mpirt;C:\Program Files (x86)\Common Files\Intel\Shared Libraries\redist\ia32\compiler;C:\Windows\system32;C:\Windows;C:\Windows\System32\Wbem;C:\Windows\System32\WindowsPowerShell\v1.0\;C:\Program Files (x86)\MATLAB\R2010a\runtime\win32;C:\Program Files (x86)\MATLAB\R2010a\bin;C:\Program Files (x86)\Bitvise SSH Client;C:\CTEX\UserData\miktex\bin;C:\CTEX\MiKTeX\miktex\bin;C:\CTEX\CTeX\ctex\bin;C:\CTEX\CTeX\cct\bin;C:\CTEX\CTeX\ty\bin;C:\CTEX\Ghostscript\gs9.05\bin;C:\CTEX\GSview\gsview;C:\CTEX\WinEdt;C:\Program Files\TortoiseSVN\bin;C:\Program Files (x86)\ATI Technologies\ATI.ACE\Core-Static;
USERNAME=dmy
OS=Windows_NT
PROCESSOR_IDENTIFIER=Intel64 Family 6 Model 42 Stepping 7, GenuineIntel
--------------- S Y S T E M ---------------
OS: Windows NT 6.2 Build 9200
CPU:total 4 (8 cores per cpu, 2 threads per core) family 6 model 10 stepping 7, cmov, cx8, fxsr, mmx, sse, sse2, sse3, ssse3, ht
Memory: 4k page, physical 4194303k(4194303k free), swap 4194303k(4194303k free)
vm_info: Java HotSpot(TM) Client VM (11.2-b01) for windows-x86 JRE (1.6.0_12-b04), built on Jan 17 2009 09:57:14 by "java_re" with MS VC++ 7.1
time: Wed Sep 03 22:06:16 2014
elapsed time: 6267 seconds

96
loadprofile.txt
View File

@ -1,96 +0,0 @@
57697.32 576.9732 663.51918 490.42722
55270.921 552.70921 635.6155915 469.8028285
53275.549 532.75549 612.6688135 452.8421665
51147.685 511.47685 588.1983775 434.7553225
49132.866 491.32866 565.027959 417.629361
47594.538 475.94538 547.337187 404.553573
45822.028 458.22028 526.953322 389.487238
44106.141 441.06141 507.2206215 374.9021985
42717.381 427.17381 491.2498815 363.0977385
42185.436 421.85436 485.132514 358.576206
40721.922 407.21922 468.302103 346.136337
40139.269 401.39269 461.6015935 341.1837865
38901.576 389.01576 447.368124 330.663396
37535.238 375.35238 431.655237 319.049523
37234.754 372.34754 428.199671 316.495409
36277.407 362.77407 417.1901805 308.3579595
35705.188 357.05188 410.609662 303.494098
35364.886 353.64886 406.696189 300.601531
34644.62 346.4462 398.41313 294.47927
34811.66 348.1166 400.33409 295.89911
34473.493 344.73493 396.4451695 293.0246905
35018.258 350.18258 402.709967 297.655193
35782.508 357.82508 411.498842 304.151318
35520.4 355.204 408.4846 301.9234
37126.94 371.2694 426.95981 315.57899
37852.994 378.52994 435.309431 321.750449
39096.988 390.96988 449.615362 332.324398
40615.023 406.15023 467.0727645 345.2276955
42599.93 425.9993 489.899195 362.099405
44509.964 445.09964 511.864586 378.334694
49173.813 491.73813 565.4988495 417.9774105
53864.107 538.64107 619.4372305 457.8449095
60260.101 602.60101 692.9911615 512.2108585
65195.266 651.95266 749.745559 554.159761
71841.983 718.41983 826.1828045 610.6568555
77802.345 778.02345 894.7269675 661.3199325
88879.392 888.79392 1022.113008 755.474832
100781.676 1007.81676 1158.989274 856.644246
111517.636 1115.17636 1282.452814 947.899906
119678.504 1196.78504 1376.302796 1017.267284
126802.705 1268.02705 1458.231108 1077.822993
131988.024 1319.88024 1517.862276 1121.898204
134124.319 1341.24319 1542.429669 1140.056712
137238.631 1372.38631 1578.244257 1166.528364
138993.205 1389.93205 1598.421858 1181.442243
140492.797 1404.92797 1615.667166 1194.188775
141981.796 1419.81796 1632.790654 1206.845266
142679.116 1426.79116 1640.809834 1212.772486
143682.364 1436.82364 1652.347186 1221.300094
144526.194 1445.26194 1662.051231 1228.472649
146089.742 1460.89742 1680.032033 1241.762807
146519.133 1465.19133 1684.97003 1245.412631
144834.953 1448.34953 1665.60196 1231.097101
145784.897 1457.84897 1676.526316 1239.171625
144503.807 1445.03807 1661.793781 1228.28236
145438.77 1454.3877 1672.545855 1236.229545
145203.708 1452.03708 1669.842642 1234.231518
143780.956 1437.80956 1653.480994 1222.138126
142645.898 1426.45898 1640.427827 1212.490133
143466.557 1434.66557 1649.865406 1219.465735
141922.25 1419.2225 1632.105875 1206.339125
145841.261 1458.41261 1677.174502 1239.650719
144923.155 1449.23155 1666.616283 1231.846818
145430.575 1454.30575 1672.451613 1236.159888
145334.611 1453.34611 1671.348027 1235.344194
146149.143 1461.49143 1680.715145 1242.267716
146482.33 1464.8233 1684.546795 1245.099805
147876.05 1478.7605 1700.574575 1256.946425
149817.698 1498.17698 1722.903527 1273.450433
150415.094 1504.15094 1729.773581 1278.528299
154442.119 1544.42119 1776.084369 1312.758012
158010.804 1580.10804 1817.124246 1343.091834
160774.574 1607.74574 1848.907601 1366.583879
162786.258 1627.86258 1872.041967 1383.683193
161608.648 1616.08648 1858.499452 1373.673508
159454.941 1594.54941 1833.731822 1355.366999
158027.837 1580.27837 1817.320126 1343.236615
157666.495 1576.66495 1813.164693 1340.165208
156460.303 1564.60303 1799.293485 1329.912576
155215.52 1552.1552 1784.97848 1319.33192
152007.853 1520.07853 1748.09031 1292.066751
149345.228 1493.45228 1717.470122 1269.434438
147650.189 1476.50189 1697.977174 1255.026607
145332.662 1453.32662 1671.325613 1235.327627
142143.342 1421.43342 1634.648433 1208.218407
136555.706 1365.55706 1570.390619 1160.723501
134474.784 1344.74784 1546.460016 1143.035664
131788.941 1317.88941 1515.572822 1120.205999
125144.813 1251.44813 1439.16535 1063.730911
116612.327 1166.12327 1341.041761 991.2047795
105557.928 1055.57928 1213.916172 897.242388
90041.36 900.4136 1035.47564 765.35156
79091.837 790.91837 909.5561255 672.2806145
68164.237 681.64237 783.8887255 579.3960145
61305.02 613.0502 705.00773 521.09267
55158.221 551.58221 634.3195415 468.8448785

BIN
mVolt.mat
View File

Binary file not shown.

4
maxBoundError.txt
View File

@ -1,4 +0,0 @@
0.010149012 0.000145532
0.010613753 0.001269963
0.010686824 0.001667111
0.01098127 0.002552894

33
maxError.txt
View File

@ -1,33 +0,0 @@
0.001940128 0.001207832
0.002036094 0.00080888
0.002099042 0.000927452
0.002424885 0.000935623
0.001893664 0.001181456
0.002018829 0.001257517
0.00210676 0.000743724
0.002013494 0.000936731
0.002013003 0.001095999
0.001842267 0.000931757
0.002001811 0.000859772
0.001866919 0.000876312
0.001834921 0.00075514
0.001937562 0.000923904
0.00189722 0.000822962
0.001933966 0.000901072
0.001920252 0.000994934
0.001747249 0.000857745
0.002008556 0.000823106
0.002093043 0.000753621
0.001971762 0.000883769
0.002987941 0.00129286
0.002436167 0.000997165
0.002037594 0.001015502
0.002066403 0.000954605
0.00276576 0.001295141
0.001881581 0.001290146
0.00193955 0.001351513
0.002069665 0.001408808
0.001936768 0.001467117
0.001780607 0.001485602
0.001861682 0.00143783
0.001955221 0.001473342

18
modifyadmmatrix.asv Normal file
View File

@ -0,0 +1,18 @@
function [new_G,new_B,GB,Y,r,c,Angle] = modifyadmmatrix(ii,jj,G,GG,B,BB)
%**************************************************************************
% 程序功能 : 子函数——形成节点导纳矩阵Y
% 编 者：
% 编制时间：2010.12
%**************************************************************************
%% 支路导纳计算
new_G=G;
new_G(ii,jj)=new_G(ii,jj)-G(ii,jj);
new_G(jj,ii)=new_G(jj,ii)-G(jj,ii);
new_G(ii,ii)=new_G(ii,ii)+G(ii,jj);
new_G(jj,jj)=new_G(jj,jj)+G(ii,jj);
%% 化作极坐标形式
GB = G+B.*1i; %将电导，电纳合并，写成复数形式
Y = abs(GB); %求节点导纳幅值
[r,c] = find(Y);
Angle = angle(GB(GB~=0)); %求节点导纳角度

64
openfile.asv Normal file
View File

@ -0,0 +1,64 @@
function [Busnum,Balance,PQstandard,Precision,Linei,Linej,Liner,Linex,Lineb,kmax,Transfori ,...
Transforj,Transforr,Transforx,Transfork0,Branchi,Branchb,Pointpoweri,PG,QG,PD,QD,PVi,PVu,Gen,GenU,GenL,GenC,CenterA,PGi,PVQU,PVQL] = openfile(FileName)
%**************************************************************************
% 程序简介 : 子函数——读取潮流计算所需数据
% 编 者：
% 编制时间 ：2010.12
%**************************************************************************
data = dlmread(FileName); % 一次读入全部数据
Busnum= data(1,1); % 节点数
PQstandard = data(1,3); % 基准容量
kmax = data(1,4); %最大迭代次数
Precision = data(2,1); % 精度
Balance = data(3,2); % 生成1到节点号的列向量
CenterA=data(1,5); %中心参数
LineNum=data(1,2); %支路数
%% 各参数矩阵分块
zeroRow = find(data(:,1)==0); %查找第一列元素为零的行号
line = data(zeroRow(1)+1:zeroRow(2)-1,:); % 形成线路参数矩阵
ground = data(zeroRow(2)+1:zeroRow(3)-1,:); % 形成对地支路参数矩阵
tran = data(zeroRow(3)+1:zeroRow(4)-1,:); % 形成变压器参数矩阵
buspq = data(zeroRow(4)+1:zeroRow(5)-1,:); % 形成节点功率参数矩阵
PV = data(zeroRow(5)+1:zeroRow(6)-1,:); % 形成pv节点功率参数矩阵
Gen=data(zeroRow(6)+1:zeroRow(7)-1,:);
%% 线路参数矩阵分块
Linei = line(:,2); % 节点i
Linej= line(:,3); % 节点j
Liner = line(:,4); % 线路电阻
Linex = line(:,5); % 线路电抗
Lineb = line(:,6); % b/2
%% 对地支路参数矩阵
Branchi = ground(:,1); % 对地支路节点号
Branchb = ground(:,2); % 对地支路的导纳
%% 变压器参数矩阵
Transfori = tran(:,2); % 节点i
Transforj= tran(:,3); % 节点j
Transforr = tran(:,4); % 变压器电阻
Transforx= tran(:,5); % 变压器电抗
Transfork0 = tran(:,6); % 变压器变比
%% 节点功率参数矩阵
Pointpoweri = buspq(:,1);
PG=buspq(:,2); % 发电机有功
QG=buspq(:,3); % 发电机无功
PD=buspq(:,4); % 负荷有功
QD=buspq(:,5); % 负荷无功
%% pv节点功率参数矩阵
PVi = PV(:,1); % PV节点的节点号
PVu = PV(:,2); % PV节点电压
PVQL=PV(:,3);%PV节点无功下限
PVQU=PV(:,4); %PV节点无功上限
%% 发电机参数
%GenU=Gen(:,[1 5 6]);
%GenL=Gen(:,[1 7 8]);
GenC=Gen(:,[1 2:4]);
t=GenC(:,2);
GenC(:,2)=GenC(:,4);
GenC(:,4)=t;
t=Gen(:,[1 5]);
%GenL=[t,PVQL(PVi)];
GenL=t;%有功下界
t=Gen(:,[1 6]);
%GenU=[t,PVQU(PVi)];
GenU=t;%有功上届
PGi=Gen(:,1);%发电机节点号
end

6
openfile2.m
View File

@ -56,12 +56,6 @@ QD=QD/Base;
PD=sparse(PD);
QD=sparse(QD);
%QD=PD*sqrt(1-.85^2)/.85;
%%
%增加大负荷
% PD(22)=PD(22)*65;
%%
PG=sparse(PG);
QG=sparse(QG);
%% pv节点功率参数矩阵

250
subOPF.m
View File

@ -1,250 +0,0 @@
function [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGap,isConverge]=subOPF(noMeasurei,PD0,QD0,mVolt,sigma)
tic
%% 存在问题
% 变压器变比的位置没有考虑，由于现在用的变比都是1，所以没有影响。 20130123
%%
% thesis=ForThesis(1,62);
isConverge=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:\算例\feeder33\feeder33.txt');
% 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);
% Volt=0.99*Volt;
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=170;
Precision=Precision/1;
%% 加误差
%找DG
DGi=find(PD0<0);
%
mPD=PD0;
mQD=QD0;
%负荷曲线有误
% mVolt(3)=0.9*mVolt(3);%15%
% mQD(13)=0.3*mQD(13);
%% 加估计上下界,用真实值
% 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=abs(mPD*3*sigma);
uPD=abs(mPD*3*sigma);
lQD=abs(mQD*3*sigma);
uQD=abs(mQD*3*sigma);
lVolt=abs(mVolt'*3*sigma);
uVolt=abs(mVolt'*3*sigma);
%DG
lPD(DGi)=abs(mPD(DGi)*3*sigma*0.3);
uPD(DGi)=abs(mPD(DGi)*3*sigma*0.3);
lQD(DGi)=abs(mQD(DGi)*3*sigma*0.3);
uQD(DGi)=abs(mQD(DGi)*3*sigma*0.3);
lVolt(DGi)=abs(mVolt(DGi)'*3*sigma*0.3);
uVolt(DGi)=abs(mVolt(DGi)'*3*sigma*0.3);
%% 不测量的数据
% mPD(22)=mPD(22)/65;
lPD(noMeasurei)=0.15*mPD(noMeasurei);%15%
uPD(noMeasurei)=0.15*mPD(noMeasurei);
lQD(noMeasurei)=0.15*mQD(noMeasurei);
uQD(noMeasurei)=0.15*mQD(noMeasurei);
% lPD(22)=0.17*mPD(22);%15%
% uPD(22)=0.17*mPD(22);
% lQD(22)=0.17*mQD(22);
% uQD(22)=0.17*mQD(22);
%电压伪量测错误
mVolt(noMeasurei)=1;
lVolt(noMeasurei)=0.12*mVolt(noMeasurei);%电压在0.93~1.07
uVolt(noMeasurei)=0.02*mVolt(noMeasurei);
%错误数据
%mVolt(2)=5;
% bigM=0.000003;
Vbi=sparse(0.5*ones(Busnum,1));
% Vbi(2)=1;
PDbi=sparse(0.5*ones(length(Loadi),1));
QDbi=sparse(0.5*ones(length(Loadi),1));
eps=1;
% Gap=0;
% 第一遍，算连续的值
fprintf('第1次迭代，算连续量。\n');
tic;
while(abs(Gap)>Precision*10)
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);
bigM=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
fprintf('第一遍计算时间:%f\n',toc)
% 第二遍，算离散
% [~,~,Init_Z,Init_W,Init_L,Init_U,~,~,~,RestraintCount,~,~,~,~,~,~,~,~,~]=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');
Gap=1;
% KK=0;
% eps=0.00001;
fprintf('\n');
% Init_Z=sparse(ones(1,RestraintCount));
% Init_W=sparse(-1*ones(1,RestraintCount));
% Init_L=1*sparse(ones(1,RestraintCount));
% Init_U=1*sparse(ones(1,RestraintCount));
% Init_Y=sparse(1,2*Busnum);%与学姐一致
tic;
fprintf('第2次迭代，算离散量。\n');
while(abs(Gap)>Precision*10)
if KK>kmax
break;
end
% Vbi(23)=0.1;
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);
bigM=1;
eps=Gap*0.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);
%% 开始解方程
% Vbi(Vbi>0.003)=.1;
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
fprintf('第二遍计算时间:%f\n',toc)
%% 计算最大不平衡量
AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
dP=PD+diag(Volt)*Y.*cos(AngleIJ)*Volt';
dP(Balance)=0;
dQ=QD+diag(Volt)*Y.*sin(AngleIJ)*Volt';%不为0 ，因为没有考虑电容的QG
dQ(Balance)=0;
mdP=max(dP)
mdQ=max(dQ)
%判断是否收敛
if min(abs(Vbi-.1))>5e-005
isConverge=0;
end
if KK>=kmax
isConverge=0;
end
toc
end