125 lines
4.9 KiB
Matlab
125 lines
4.9 KiB
Matlab
function [mut,mub,dNCPT,dNCPT_dt,d2NCPT_dt2,dT_Zt,transGmut,transBmut,dNCPcap,dNCP_dcap,d2NCP_dcap2,capBi_dBci,H1,N1,diagE,diagF,hx,Gap,max_hx,fx,...
|
||
Pg,Qg,xz1,dmaxOut1,xz2,dmaxOut2,xz3,dmaxOut3,xz4,dmaxOut4,PGR,Tn,kt0,capKn,q0,Tn1,capKn1,T,G]=midResult(nodeNum,pgNode,opfGoal,z,l,u,w,Pg,Qg,...
|
||
e,f,G,B,a,b,c,PG,Pd,Qd,k,xz1,dmaxOut1,xz2,dmaxOut2,xz3,dmaxOut3,xz4,dmaxOut4,times,mut,mub,oneTrans,oneCap,para3,accuracy,T,capK,transKmin,transKmax,...
|
||
capKmin,capBi,mu,transNum,transG,transB,Tn,Tn1,capKn,capKn1,kt0,q0,mismatch)
|
||
%程序功能：计算每次迭代后的目标函数值和等式约束的最大不平衡量，并根据每步迭代结果修正光滑参数
|
||
%编写时间：2010年11月
|
||
|
||
Gap=l'*z-u'*w; %计算互补间隙
|
||
PGR=PG;
|
||
if(opfGoal==2)
|
||
sumP=sum(Pg);
|
||
sumPd=sum(Pd);
|
||
fx=sumP-sumPd;
|
||
else if(opfGoal==1)
|
||
fx=sum(a(pgNode).*PGR(pgNode).*PGR(pgNode)+b(pgNode).*PGR(pgNode)+c(pgNode));
|
||
else if(opfGoal==3)
|
||
EE=e*ones(1,nodeNum);
|
||
FF=f*ones(1,nodeNum);
|
||
fx=-0.5*sum(sum(G.*((EE-EE').^2+(FF-FF').^2),2));
|
||
else if(opfGoal==4)
|
||
EE=e*ones(1,nodeNum);
|
||
FF=f*ones(1,nodeNum);
|
||
fx=0.5*sum(sum(B.*((EE-EE').^2+(FF-FF').^2),2));
|
||
end
|
||
end
|
||
end
|
||
end
|
||
H1=G*e-B*f;
|
||
N1=G*f+B*e;
|
||
diagE=sparse(1:nodeNum,1:nodeNum,e,nodeNum,nodeNum);
|
||
diagF=sparse(1:nodeNum,1:nodeNum,f,nodeNum,nodeNum);
|
||
P=Pg-Pd;
|
||
Q=Qg-Qd;
|
||
dP=P-diagE*H1-diagF*N1;
|
||
dQ=Q-diagF*H1+diagE*N1;
|
||
firstOK=0;
|
||
hx=[dP;dQ];
|
||
tmpmax_hx=max(abs(hx));
|
||
if Gap<accuracy&&tmpmax_hx<mismatch&& times==1 %判断互补间隙和最大不平衡量是否达到精度
|
||
firstOK=1;
|
||
times=2;
|
||
maxTap=5; %变压器最大调节档位
|
||
K=1./T;
|
||
Kmax=1./transKmin;
|
||
Kmin=1./transKmax;
|
||
kt0=(Kmax-Kmin)/(maxTap-1); %变压器单位调节幅度
|
||
Kn=fix((K-Kmin)./kt0).*kt0+Kmin; %获取变压器变比连续最优解左右两侧的档位所对应的变比
|
||
Kn1=Kn+kt0;
|
||
Tn=1./Kn1;
|
||
Tn1=1./Kn;
|
||
|
||
q0=1; %电容器档位的单位调节幅度
|
||
capQmin=capKmin.*capBi;
|
||
capKn=fix((capK-capQmin)./capBi).*capBi+capQmin; %获取电容电抗器投入组数连续最优解左右两侧的离散值
|
||
capKn1=capKn+capBi;
|
||
end
|
||
if times==2
|
||
%% 修正陡度参数μt
|
||
tmpmu=para3;
|
||
tmpmut=mut;
|
||
tmpmub=mub;
|
||
Ezu=exp(-(T-Tn)./tmpmut);
|
||
Ezu_1=exp(-(Tn1-T)./tmpmut);
|
||
dNCPT=-tmpmut.*log(Ezu+Ezu_1); %变压器变比对应的凝聚函数常数项
|
||
Ezbu=exp((capKn-capK)./tmpmub);
|
||
Ezbu_1=exp((capK-capKn1)./tmpmub);
|
||
dNCPcap=-tmpmub.*log(Ezbu+Ezbu_1); %电容电抗器投入组数对应的凝聚函数常数项
|
||
|
||
dNCP=[dNCPT;dNCPcap];%只有一个光滑参数
|
||
oneTC=[oneTrans;oneCap];
|
||
mutb=[mut;mub];
|
||
mutb=min(max(mu,max(abs(dNCP)))*oneTC,tmpmu*mutb.*oneTC); %修正光滑参数
|
||
mut=mutb(1:transNum);
|
||
mub=mutb(transNum+1:end);
|
||
%% 计算凝聚函数相关参数
|
||
Ezu=exp(-(T-Tn)./mut);
|
||
Ezu_1=exp(-(Tn1-T)./mut);
|
||
dNCPT=-mut.*log(Ezu+Ezu_1); %变压器变比对应的凝聚函数常数项
|
||
|
||
dT_Zt=1;
|
||
dNCPT_dt=(Ezu-Ezu_1)./(Ezu+Ezu_1); %凝聚函数对tij求一次偏导
|
||
d2NCPT_dt2=-(oneTrans-dNCPT_dt.^2)./mut; %凝聚函数对tij求二次偏导,原来的公式推导错误d2NCPT_dt2=-(oneTrans+dNCPT_dt.^2)./mut;3月7日修改
|
||
transGmut=transG;
|
||
transBmut=transB;
|
||
|
||
Ezbu=exp((capKn-capK)./mub);
|
||
Ezbu_1=exp((capK-capKn1)./mub);
|
||
dNCPcap=-mub.*log(Ezbu+Ezbu_1); %电容电抗器投入组数对应的凝聚函数常数项
|
||
|
||
dcapK_Zb=1;
|
||
dNCP_dcap=(Ezbu-Ezbu_1)./(Ezbu+Ezbu_1); %凝聚函数对capK求一次偏导
|
||
d2NCP_dcap2=-(oneCap-dNCP_dcap.^2)./mub; %凝聚函数对capK求二次偏导
|
||
capBi_dBci=capBi.*dcapK_Zb;
|
||
end
|
||
if times==1
|
||
dNCPT=[];dNCPT_dt=[];d2NCPT_dt2=[];dT_Zt=[];transGmut=[];transBmut=[];dNCPcap=[];dNCP_dcap=[];d2NCP_dcap2=[];capBi_dBci=[];
|
||
end
|
||
hx=[dP;dQ;dNCPT;dNCPcap];
|
||
fx=full(fx);
|
||
max_hx=full(max(abs(hx)));
|
||
if firstOK==1
|
||
max_hx=full(tmpmax_hx);
|
||
end
|
||
fprintf('%4.0f %4.4e %4.4e %7.6f\n',k,Gap,max_hx,fx); %输出每次迭代结果
|
||
%% 保留每次迭代的互补间隙和最大不平衡量
|
||
if times==1
|
||
xz1=[xz1 k];
|
||
dmaxOut1=[dmaxOut1 Gap];
|
||
xz3=[xz3 k];
|
||
dmaxOut3=[dmaxOut3 max_hx];
|
||
end
|
||
if times==2
|
||
if firstOK==1
|
||
xz1=[xz1 k];
|
||
dmaxOut1=[dmaxOut1 Gap];
|
||
xz3=[xz3 k];
|
||
dmaxOut3=[dmaxOut3 max_hx];
|
||
firstOK=0;
|
||
else
|
||
xz2=[xz2 k];
|
||
dmaxOut2=[dmaxOut2 Gap];
|
||
xz4=[xz4 k];
|
||
dmaxOut4=[dmaxOut4 max_hx];
|
||
end
|
||
end |