一、简介
法。而且无论是面对LP还是QP，内点法都显示出了相当的极好的性能，例如多项式的算法复杂度。
二、源代码
clear;
clc;
errArr=[];
%%
%初始化！！！
initial;
% Start clock
t1 = clock;
%%
ROU=sl'*MU_MIN+su'*MU_MAX;
MUt=SIGMA*ROU/(2*length(sl));%初始对偶因子与惩罚因子计算%
ik=0;%计迭代次数！！！
%迭代循环过程！！
while(abs(ROU)>=err)
%%
%Calcute h,g matrix
ROU=sl'*MU_MIN+su'*MU_MAX;
errArr=[errArr;ROU;];
SIGMA=0;
MU=SIGMA*ROU/(2*length(sl)); %中心参数置零%
for i=1:30
temp=0;
for j=1:30
temp=temp-V(j)*aY(i,j)*cos(Vth(i)-Vth(j)-Yth(i,j));
end
if (i>6)
tPg=0;
else
tPg=Pg(i);
end
h(i)=tPg-Pd(i)+V(i)*temp;
end
for i=1:30
temp=0;
for j=1:30
temp=temp-V(j)*aY(i,j)*sin(Vth(i)-Vth(j)-Yth(i,j));
end
if (i>6)
tQg=0;
else
tQg=Qg(i);
end
h(i+30)=tQg-Qd(i)+V(i)*temp;
end % Cal h END
for i=1:6
g(i)=Pg(i);
g(i+6)=Qg(i);
end
for i=1:30
g(i+12)=V(i);
end % Cal g END
%Calcute h,g matrix END
%%
%Calculate Jacobian&Hessian matix
%First Step: Jf,Hf
for i=1:6
Jf(i)=2*gencost(i,5)*Pg(i)+gencost(i,6);
Hf(i,i)=2*gencost(i,6);
end
%Second Step: Jh, h为等式约束
for i=1:6 %前6行对Pg求导，由此已求出
Jh(i,i)=1;
end
for i=7:12 %7-12行对Qg求导，由此已求出
Jh(i,i+24)=1;
end
for i=1:30 %形成13-42行的1-60列
for j=1:30
tempVp=0;
tempVq=0;
if (j==i)
for k=1:30
tempVp=tempVp-V(k)*aY(j,k)*cos(Vth(j)-Vth(k)-Yth(j,k));
tempVq=tempVq-V(k)*aY(j,k)*sin(Vth(j)-Vth(k)-Yth(j,k));
end
Jh(12+j,i)=tempVp-aY(j,j)*V(j)*cos(Yth(j,j));
Jh(12+j,30+i)=tempVq+aY(j,j)*V(j)*sin(Yth(j,j));
else
Jh(12+j,i)=-aY(i,j)*V(i)*cos(Vth(i)-Vth(j)-Yth(i,j));
Jh(12+j,30+i)=-aY(i,j)*V(i)*sin(Vth(i)-Vth(j)-Yth(i,j));
end
end
end
for i=1:30 %形成43-72行的1-60列
for j=1:30
tempVp=0;
tempVq=0;
if (j==i)
for k=1:30
tempVp=tempVp+aY(j,k)*V(k)*sin(Vth(j)-Vth(k)-Yth(j,k));
tempVq=tempVq-aY(j,k)*V(k)*cos(Vth(j)-Vth(k)-Yth(j,k));
end
tempVp=tempVp-V(j)*aY(j,j)*sin(-Yth(j,j));
tempVq=tempVq+V(j)*aY(j,j)*cos(-Yth(j,j));
Jh(42+j,i)=V(i)*tempVp;
Jh(42+j,30+i)=V(i)*tempVq;
else
Jh(42+j,i)=-aY(i,j)*V(i)*V(j)*sin(Vth(i)-Vth(j)-Yth(i,j));
Jh(42+j,30+i)=aY(i,j)*V(i)*V(j)*cos(Vth(i)-Vth(j)-Yth(i,j));
end
end
end
%Third Step: Hh
%有功部分
for i=1:30
for j=1:30
for k=j:30
if (j==k&&i~=j)
Hh(j+12,k+12,i)=0; %VV
Hh(j+42,k+42,i)=V(i)*aY(i,j)*V(j)*cos(Vth(i)-Vth(j)-Yth(i,j)); %%thth
elseif (j==k&&i==j)
Hh(j+12,k+12,i)=-2*aY(j,j)*cos(Yth(i,i)); %VV
temp=0; %thth
for l=1:30
temp=temp+aY(j,l)*V(l)*cos(Vth(j)-Vth(l)-Yth(j,l));
end
temp=temp-aY(i,i)*V(i)*cos(-Yth(i,i));
Hh(j+42,k+42,i)=V(i)*temp;
elseif (k==i)
Hh(j+12,k+12,i)=-aY(i,j)*cos(Vth(i)-Vth(j)-Yth(i,j)); %VV
Hh(k+12,j+12,i)=Hh(j+12,k+12,i);
Hh(j+42,k+42,i)=-V(i)*aY(i,j)*V(j)*cos(Vth(i)-Vth(j)-Yth(i,j)); %thth
Hh(k+42,j+42,i)=Hh(j+42,k+42,i);
elseif (j==i)
Hh(j+12,k+12,i)=-aY(i,k)*cos(Vth(i)-Vth(k)-Yth(i,k)); %VV
Hh(k+12,j+12,i)=Hh(j+12,k+12,i);
Hh(j+42,k+42,i)=-V(i)*aY(i,k)*V(k)*cos(Vth(i)-Vth(k)-Yth(i,k)); %thth
Hh(k+42,j+42,i)=Hh(j+42,k+42,i);
end
end
end
end %至此已形成（13-42，13-42）和（42-72，43-72）
for i=1:30
for j=1:30
for k=1:30
if (j==k&&i~=j)
Hh(j+42,k+12,i)=-V(i)*aY(i,j)*sin(Vth(i)-Vth(j)-Yth(i,j)); %thV
elseif (j==k&&i==j)
temp=0; %thV
for l=1:30
temp=temp+aY(j,l)*V(l)*sin(Vth(j)-Vth(l)-Yth(j,l));
end
Hh(j+42,k+12,i)=temp-V(i)*aY(i,i)*sin(-Yth(i,i));
elseif (j==i)
Hh(j+42,k+12,i)=V(i)*aY(i,k)*sin(Vth(i)-Vth(k)-Yth(i,k)); %thV
elseif (k==i)
Hh(j+42,k+12,i)=-V(j)*aY(i,j)*sin(Vth(i)-Vth(j)-Yth(i,j)); %thV
end
end
end
Hh(13:42,43:72,i)=Hh(43:72,13:42,i)';
end %至此已形成（42-72，13-42）和（13-42，43-72）
%无功部分
for i=1:30
for j=1:30
for k=j:30
if (j==k&&i~=j)
Hh(j+12,k+12,i+30)=0; %VV
Hh(j+42,k+42,i+30)=V(i)*aY(i,j)*V(j)*sin(Vth(i)-Vth(j)-Yth(i,j)); %%thth
elseif (j==k&&i==j)
Hh(j+12,k+12,i+30)=2*aY(j,j)*sin(Yth(i,i)); %VV
temp=0; %thth
for l=1:30
temp=temp+aY(j,l)*V(l)*sin(Vth(j)-Vth(l)-Yth(j,l));
end
temp=temp-aY(i,i)*V(i)*sin(-Yth(i,i));
Hh(j+42,k+42,i+30)=V(i)*temp;
elseif (k==i)
Hh(j+12,k+12,i+30)=-aY(i,j)*sin(Vth(i)-Vth(j)-Yth(i,j)); %VV
Hh(k+12,j+12,i+30)=Hh(j+12,k+12,i+30);
Hh(j+42,k+42,i)=-V(i)*aY(i,j)*V(j)*sin(Vth(i)-Vth(j)-Yth(i,j)); %thth
Hh(k+42,j+42,i+30)=Hh(j+42,k+42,i+30);
elseif (j==i)
Hh(j+12,k+12,i+30)=-aY(i,k)*sin(Vth(i)-Vth(k)-Yth(i,k)); %VV
Hh(k+12,j+12,i+30)=Hh(j+12,k+12,i+30);
Hh(j+42,k+42,i+30)=-V(i)*aY(i,k)*V(k)*sin(Vth(i)-Vth(k)-Yth(i,k)); %thth
Hh(k+42,j+42,i+30)=Hh(j+42,k+42,i+30);
end
end
end
end %至此已形成（13-42，13-42）和（42-72，43-72）
for i=1:30
for j=1:30
for k=1:30
if (j==k&&i~=j)
Hh(j+42,k+12,i+30)=V(i)*aY(i,j)*cos(Vth(i)-Vth(j)-Yth(i,j)); %thV
elseif (j==k&&i==j)
temp=0; %thV
for l=1:30
temp=temp-aY(j,l)*V(l)*cos(Vth(j)-Vth(l)-Yth(j,l));
end
Hh(j+42,k+12,i+30)=temp+V(i)*aY(i,i)*cos(-Yth(i,i));
elseif (j==i)
Hh(j+42,k+12,i+30)=-V(i)*aY(i,k)*cos(Vth(i)-Vth(k)-Yth(i,k)); %thV
elseif (k==i)
Hh(j+42,k+12,i+30)=V(j)*aY(i,j)*cos(Vth(i)-Vth(j)-Yth(i,j)); %thV
end
end
end
Hh(13:42,43:72,i+30)=Hh(43:72,13:42,i+30)';
end %至此已形成（42-72，13-42）和（13-42，43-72）
%Hh形成完毕
%Fourth Step: Jg, Hg
Jg=eye(42,42);
Jg=[Jg;zeros(30,42)];
Hg=zeros(72);
%Calculation Jacobian&Hessian matrix END
%%
%Calculate Newton Iteration 误差迭代量
%Cal LX0-------------------------1
LX0=Jf-Jh*Lam+Jg*(-MU_MIN+MU_MAX);
%Cal LLam-------------------------2
LLam0=h;
pferr=max(LLam0);
%Cal LMU_MIN-------------------------3
LMU_MIN0=g-sl-gmin;
%Cal LMU_MAX-------------------------4
三、运行结果
四、备注
版本：2014a
