【優化求解】基於matlab內點法求解實時電價最優問題【含Matlab原始碼 1161期】
阿新 • • 發佈:2021-07-29
一、簡介
法。而且無論是面對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
完整程式碼或代寫加QQ1564658423