关于matpower中目标函数修改的问题
如何修改matpower的目标函数或者重建一个目标函数重新运行一下？求大神指教！！！ 有没有例程贴上来看一下 你要把它的程序读懂，相应部分都要修改的 回复 2# redplum首先谢谢指教
下面是这个程序，指教一下那里需要修改啊clear;clc;errArr=[];%%%³õʼ»¯£¡£¡£¡initial;% Start clockt1 = 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=; SIGMA=0; MU=SIGMA*ROU/(2*length(sl)); %ÖÐÐIJÎÊýÖÃÁã%
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
% Calh 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
% CalgEND
%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=; 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 LMU_MAX0=g+su-gmax;
%Cal Lsl-------------------------5 Lsl0=diag(MU_MIN)*diag(sl)*ones(length(sl),1)-MU*ones(length(sl),1);
%CAl Lsu-------------------------6 Lsu0=diag(MU_MAX)*diag(su)*ones(length(su),1)-MU*ones(length(su),1);
%Calculate Newton Iteration Îó²îµü´úÁ¿ END!!!
%%
%Calculate Newton Iteration ·ÂÉäÐÞÕýÁ¿
%1st Step: ·ÂÉäÐÞÕýÁ¿delXaf,·ÂÉäÐÞÕýÁ¿delLamaf
%S1:H temp=0;%֮ǰÒÑÓùýÁÙʱ±äÁ¿temp£¬ÔÚ´ËÇåÁã
for i=1:60 temp=temp+Lam(i)*Hh(:,:,i);
end tempgMUg1=0; tempgMUg1=Jg*diag(MU_MIN)*Jg';
for i=1:42 tempgMUg1(i,:)=tempgMUg1(i,:)/sl(i);
end tempgMUg2=0; tempgMUg2=Jg*diag(MU_MAX)*Jg';
for i=1:42 tempgMUg2(i,:)=tempgMUg2(i,:)/su(i);
end H=Hf-temp+tempgMUg1+tempgMUg2;
%S2:·ÂÉäKESAaf tempgMUg1=0; tempgMUg1=diag(MU_MIN)*sl+diag(MU_MIN)*LMU_MIN0;
for i=1:42 tempgMUg1(i)=tempgMUg1(i)/sl(i);
end tempgMUg1=Jg*tempgMUg1; tempgMUg2=0; tempgMUg2=diag(MU_MAX)*su-diag(MU_MAX)*LMU_MAX0;
for i=1:42 tempgMUg2(i)=tempgMUg2(i)/su(i);
end tempgMUg2=Jg*tempgMUg2; KESAaf=LX0+tempgMUg1-tempgMUg2;
%S3:JACOBIAN JACOBIAN=;
%S4£ºRESULT RESULT=-;
%S5:Cal ·ÂÉäÐÞÕýÁ¿delXaf,·ÅÉäÐÞÕýÁ¿delLamaf delXaf=-Jh'\LLam0; delLamaf=Jh\(H*delXaf+KESAaf);% temp=JACOBIAN\RESULT;% for i=1:72% delX(i)=temp(i);% end% for i=1:60% delLam(i)=temp(i+72);% end
%2nd Step: ·ÅÉäËɳڱäÁ¿ÐÞÕýÁ¿delslaf, delsuaf
%S1: delslaf delslaf=Jg'*delXaf+LMU_MIN0;
%S2: delsuaf delsuaf=-Jg'*delXaf-LMU_MAX0;
%3rd Step: ·ÅÉäÀ ¸ñÀÊÈÕ³Ë×ÓÐÞÕýÁ¿delMU_MINaf,delMU_MAXaf
%S1: delMU_MINaf temp=0; temp=-Lsl0-diag(MU_MIN)*delslaf;
for i=1:42 temp(i)=temp(i)/sl(i);
end temp(13)=0; delMU_MINaf=temp;
%S2: delMU_MAXaf temp=0; temp=-Lsu0-diag(MU_MAX)*delsuaf;
for i=1:42 temp(i)=temp(i)/su(i);
end temp(13)=0; delMU_MAXaf=temp;
%Calculate Newton Iteration ·ÂÉäÐÞÕýÁ¿ END!!!!!!!!!!!!!!
%%
%¼ÆËã·ÂÉäÔ Ê¼²½³¤ºÍ¶Ôż²½³¤ STEPpaf£¬ STEPdaf
%S1: STEPpaf
for i=1:42
if (delslaf(i)~=0) temp1=-sl(i)/delslaf(i);
else temp1=Inf;
end
if (delsuaf(i)~=0) temp2=-su(i)/delsuaf(i);
else temp2=Inf;
end
if temp1<temp2 min=temp1;
else min=temp2;
end
if min>1 min=1;
end
end STEPpaf=0.9995*min;
%S2: STEPdaf
for i=1:42 temp1=-MU_MIN(i)/delMU_MINaf(i); temp2=-MU_MAX(i)/delMU_MAXaf(i);
if temp1<temp2 min=temp1;
else min=temp2;
end
if (i==13) min=0;
end
if min>1 min=1;
end
end STEPdaf=0.9995*min;
%¼ÆËã·ÂÉäÔ Ê¼²½³¤ºÍ¶Ôż²½³¤ STEPp£¬ STEPd END!!!!!
%¼ÆËã·ÂÉä¶ÔżÒòÊý¼°³Í·£Òò×Ó% ROUaf=(sl+STEPpaf*delslaf)'*(MU_MIN+STEPdaf*delMU_MINaf)+(su+STEPpaf*delsuaf)'*(MU_MAX+STEPdaf*delMU_MAXaf);
if (ROUaf/ROU)^2<0.2 SIGMA=(ROUaf/ROU)^2;
else SIGMA=0.2;
end MUaf=SIGMA*ROUaf/(2*length(sl));
%¼ÆËãÍê±Ï% Lsl0=diag(MU_MIN)*diag(sl)*ones(length(sl),1)-MUaf*ones(length(sl),1)-diag(delslaf)*delMU_MINaf; Lsu0=diag(MU_MAX)*diag(su)*ones(length(su),1)-MUaf*ones(length(su),1)-diag(delsuaf)*delMU_MAXaf;
%Calculate Newton Iteration ÐÞÕýÁ¿
%1st Step: УÕýÐÞÕýÁ¿delX,ÐÞÕýÁ¿delLam
%S1:H temp=0;%֮ǰÒÑÓùýÁÙʱ±äÁ¿temp£¬ÔÚ´ËÇåÁã
for i=1:60 temp=temp+Lam(i)*Hh(:,:,i);
end tempgMUg1=0; tempgMUg1=Jg*diag(MU_MIN)*Jg';
for i=1:42 tempgMUg1(i,:)=tempgMUg1(i,:)/sl(i);
end tempgMUg2=0; tempgMUg2=Jg*diag(MU_MAX)*Jg';
for i=1:42 tempgMUg2(i,:)=tempgMUg2(i,:)/su(i);
end H=Hf-temp+tempgMUg1+tempgMUg2;
%S2:УÕýKESA tempgMUg1=0; tempgMUg1=Lsl0+diag(MU_MIN)*LMU_MIN0+diag(delslaf)*delMU_MINaf;
for i=1:42 tempgMUg1(i)=tempgMUg1(i)/sl(i);
end tempgMUg1=Jg*tempgMUg1; tempgMUg2=0; tempgMUg2=Lsu0-diag(MU_MAX)*LMU_MAX0+diag(delsuaf)*delMU_MAXaf;
for i=1:42 tempgMUg2(i)=tempgMUg2(i)/su(i);
end tempgMUg2=Jg*tempgMUg2; KESA=LX0+tempgMUg1-tempgMUg2;
%S3:JACOBIAN JACOBIAN=;
%S4£ºRESULT RESULT=-;
%S5:Cal УÕýÐÞÕýÁ¿delX,ÐÞÕýÁ¿delLam delX=-Jh'\LLam0; delLam=Jh\(H*delX+KESA);% temp=JACOBIAN\RESULT;% for i=1:72% delX(i)=temp(i);% end% for i=1:60% delLam(i)=temp(i+72);% end
%2nd Step: УÕýËɳڱäÁ¿ÐÞÕýÁ¿delsl, delsu
%S1: delsl delsl=Jg'*delX+LMU_MIN0;
%S2: delsu delsu=-Jg'*delX-LMU_MAX0;
%3rd Step: УÕýÀ ¸ñÀÊÈÕ³Ë×ÓÐÞÕýÁ¿delMU_MIN,delMU_MAX
%S1: delMU_MIN temp=0; temp=-diag(MU_MIN)*sl+MUt*ones(42,1)-diag(delMU_MINaf)*delslaf-diag(MU_MIN)*delslaf;
for i=1:42 temp(i)=temp(i)/sl(i);
end temp(13)=0; delMU_MIN=temp;
%S2: delMU_MAX temp=0; temp=-diag(MU_MAX)*su+MUt*ones(42,1)-diag(delMU_MAXaf)*delsuaf-diag(MU_MAX)*delsuaf;
for i=1:42 temp(i)=temp(i)/su(i);
end temp(13)=0; delMU_MAX=temp;
%Calculate Newton Iteration ÐÞÕýÁ¿ END!!!!!!!!!!!!!!
%%
%¼ÆËãУÕýÔ Ê¼²½³¤ºÍ¶Ôż²½³¤ STEPp£¬ STEPd
%S1: STEPp
for i=1:42
if (delslaf(i)~=0) temp1=-sl(i)/delsl(i);
else temp1=Inf;
end
if (delsuaf(i)~=0) temp2=-su(i)/delsu(i);
else temp2=Inf;
end
if temp1<temp2 min=temp1;
else min=temp2;
end
if min>1 min=1;
end
end STEPp=0.9995*min;
%S2: STEPd
for i=1:42 temp1=-MU_MIN(i)/delMU_MIN(i); temp2=-MU_MAX(i)/delMU_MAX(i);
if temp1<temp2 min=temp1;
else min=temp2;
end
if (i==13) min=0;
end
if min>1 min=1;
end
end STEPd=0.9995*min;
%¸üÐÂÔ Ê¼±äÁ¿ºÍ¶Ôż±äÁ¿ X=X+STEPp*delX; Lam=Lam+STEPd*delLam; sl=sl+STEPp*delsl; su=su+STEPp*delsu; MU_MIN=MU_MIN+STEPd*delMU_MIN; MU_MAX=MU_MAX+STEPd*delMU_MAX; Pg=X(1:6);Qg=X(7:12);V=X(13:42);Vth=X(43:72);
%¸üÐÂÔ Ê¼±äÁ¿ºÍ¶Ôż±äÁ¿ END!!!
%¼ÆËã¶Ôż¼ä϶% ROU=sl'*MU_MIN+su'*MU_MAX; MUt=SIGMA*ROU/(2*length(sl));
%%
%ÅжÏÊÇ·ñ³¬¹ýµü´ú´ÎÊý ik=ik+1;
if ik>IterNumMAX disp('IterNum ERROR!!!!!');
break;
end
endet = etime(clock, t1);%%%Êä³ö²¿·Ö
if (ik<=IterNumMAX)
%%
%ÇóµÃ×î´ó×îСµÈʽÎó²î F=0;for i=1:6 F=F+gencost(i,5)*(X(i)*baseMVA)^2+gencost(i,6)*X(i)*baseMVA;end
max=-10;min=10;for i=1:ik temp=LLam0(i);
if (temp>max) max=temp;
end
if (temp<min) min=temp;
endendmaxik%%%Êä³ö½á¹û fprintf('\nConverged in %.2f seconds', et); fprintf('\nObjective Function Value = %.2f $/hr', F); fprintf('\n================================================================================'); fprintf('\n| Bus Data |'); fprintf('\n================================================================================'); fprintf('\n Bus Voltage Generation Load '); fprintf('Lambda($/MVA-hr)'); fprintf('\n# Mag(pu) Ang(deg) P (MW) Q (MVAr) P (MW) Q (MVAr)'); fprintf(' P Q '); fprintf('\n----- ------- ----------------------------------------'); fprintf('--------------');
for i = 1:30 fprintf('\n%5d%7.3f%9.3f', i,X(i+12),X(i+42)*180/pi);
if (i<=6) fprintf('%10.2f%10.2f', X(i)*baseMVA, X(i+6)*baseMVA);
else fprintf(' - -');
end fprintf('%10.2f %9.2f', bus(i, 3) , bus(i, 4) ); fprintf('%9.3f %9.3f', Lam(i),Lam(i+30));
if (i==30) fprintf('\n --------------------------------'); fprintf('\n Total: %9.2f %9.2f',(X(1)+X(2)+X(3)+X(4)+X(5)+X(6))*baseMVA,(X(7)+X(8)+X(9)+X(10)+X(11)+X(12))*baseMVA) fprintf('\n');
end
end
%%%%Êä³öÇúÏßÄâºÏmainiterNum=';%axis();hold on;f = fit(iterNum,errArr,'spline');f=feval(f,iterNum);plot(iterNum,errArr,'o',iterNum,f,'-');end 回复 3# qcsun
谢谢，忙看一下哪里需要修改嘛，程序在二楼 没看懂这个是做什么的。 matpower里有自带的m文件，直接调用即可
页:
[1]
