【優(yōu)化充電】基于粒子群算法求解IEEE3環(huán)境下4個充電站電動汽車群有序充電優(yōu)化策略問題
??作者簡介：熱愛科研的Matlab仿真開發(fā)者，修心和技術同步精進，
代碼獲取、論文復現及科研仿真合作可私信。
??個人主頁：Matlab科研工作室
??個人信條：格物致知。
更多Matlab完整代碼及仿真定制內容點擊??
?? ? ???? ? ??? ? ??? ? ? ?
?? ? ? ? ? ? ??? ? ? ? ? ? ???? ? ???? ? ? ?
?? 內容介紹
粒子群算法（PSO）是一種用于解決優(yōu)化問題的啟發(fā)式算法，它模擬了鳥群或魚群的行為，通過不斷調整粒子的位置來尋找最優(yōu)解。在電動汽車充電優(yōu)化領域，粒子群算法被廣泛應用于解決充電站位置優(yōu)化、充電樁調度等問題。本文將介紹粒子群算法在IEEE3環(huán)境下4個充電站電動汽車群有序充電優(yōu)化策略問題的研究。
首先，我們需要了解IEEE3環(huán)境下的充電站電動汽車群有序充電優(yōu)化策略問題。在這個問題中，我們需要考慮4個充電站和一組電動汽車，如何安排充電順序，以最大程度地減少充電時間和充電成本。這是一個典型的優(yōu)化問題，需要找到最佳的充電順序和充電站位置，以滿足汽車的充電需求，同時最大限度地減少充電時間和成本。
接下來，我們將介紹粒子群算法在解決這個問題上的應用。粒子群算法通過不斷調整粒子的位置和速度來搜索最優(yōu)解。在這個問題中，我們可以將每個粒子看作是一種充電順序的解，通過不斷更新粒子的位置和速度，最終找到最優(yōu)的充電順序和充電站位置。粒子群算法具有全局尋優(yōu)能力和較好的收斂性，能夠有效地解決這個優(yōu)化問題。
在實際應用中，我們可以將粒子群算法與充電站電動汽車群有序充電優(yōu)化策略問題相結合，通過編寫相應的優(yōu)化算法，實現對充電順序和充電站位置的優(yōu)化。通過實驗和仿真，我們可以驗證粒子群算法在這個問題上的有效性和性能優(yōu)勢，為實際應用提供可靠的優(yōu)化方案。
綜上所述，粒子群算法在IEEE3環(huán)境下4個充電站電動汽車群有序充電優(yōu)化策略問題的研究具有重要意義。通過將粒子群算法與充電優(yōu)化問題相結合，可以有效地解決充電順序和充電站位置的優(yōu)化問題，為電動汽車充電領域的發(fā)展提供有力支持。希望本文的介紹能夠對相關領域的研究和實踐工作有所幫助，為粒子群算法在電動汽車充電優(yōu)化領域的應用提供新的思路和方法。
?? 部分代碼
%相關原始數據格式說明如下：%n——節(jié)點個數；n1——支路條數；isb——平衡節(jié)點號；H——PQ節(jié)點個數（為后面形成PVU存儲PV節(jié)點初始電壓用）；pr——誤差精度。%B1——支路參數矩陣，其中第一列和第二列是起始節(jié)點編號和終點節(jié)點編號,第三列、第四列、第五列、第六列分別為：支路電阻、電抗、變壓器變比、電納。（不考慮電導）%B2——節(jié)點參數矩陣，其中第一列和第二列為節(jié)點編號和節(jié)點類型；第三列到第六列分別為：注入有功、注入無功、電壓幅值、電壓相位。%節(jié)點類型分類如下：“0”為平衡節(jié)點,“1”為PQ，“2”為PV節(jié)點；“3”為PQ（V）節(jié)點，“4”為PI節(jié)點。function ?[V_flu,Ploss]=IEEE33(car_1,car_2,car_3,P_flex)n=33 ; ? ? ? n1=32;isb=1;H=32; ? ? ? ? ? ? ? ?%%%%%%%%%%%%%18節(jié)點加DG ? ?PQV處理pr=0.01; v_amp=0;B1=[1 2 0.00922 0.0047i 1 0; ? ?2 3 0.00493 0.02511i 1 0; ? ?3 4 0.0366 0.01864i 1 0; ? ?4 5 0.03811 0.01941i 1 0; ? ?5 6 0.0819 0.0707i 1 0; ? ?6 7 0.01872 0.06188i 1 0; ? ?7 8 0.07114 0.02351i 1 0; ? ?8 9 0.103 0.074i 1 0; ? ?9 10 0.1044 0.074i 1 0; ? ?10 11 0.01966 0.0065i 1 0; ? ?11 12 0.03744 0.01238i 1 0; ? ?12 13 0.1468 0.1155i 1 0; ? ?13 14 0.05416 0.07129i 1 0; ? ?14 15 0.05910 0.0526i 1 0; ? ?15 16 0.07463 0.05450i 1 0; ? ?16 17 0.1289 0.1721i 1 0; ? ?17 18 0.0732 0.0574i 1 0; ? ?2 19 0.0164 0.01565i 1 0; ? ?19 20 0.15042 0.13554i 1 0; ? ?20 21 0.04095 0.04784i 1 0; ? ?21 22 0.07089 0.09373i 1 0; ? ?3 23 0.04512 0.03083i 1 0; ? ?23 24 0.08980 0.07091i 1 0; ? ?24 25 0.08960 0.07011i 1 0; ? ?6 26 0.0203 0.01034i 1 0; ? ?26 27 0.02842 0.01447i 1 0; ? ?27 28 0.1059 0.09337i 1 0; ? ?28 29 0.08042 0.07006i 1 0; ? ?29 30 0.05075 0.02585i 1 0; ? ?30 31 0.09744 0.0963i 1 0; ? ?31 32 0.03105 0.03619i 1 0; ? ?32 33 0.03410 0.05302i 1 0];B2=[1 0 0 0 1.05 0; ? ?2 1 -0.01 -0.006 1 0; ? ?3 1 -0.009 -0.004 1 0; ? ?4 1 -0.012 -0.008 1 0; ? ?5 1 -0.006 -0.003 1 0; ? ?6 1 -0.006 -0.002 1 0; ? ?7 1 -0.02 -0.01 1 0; ? ?8 1 -0.02 -0.01 1 0; ? ?9 1 -0.006 -0.002 1 0; ? ?10 1 -0.006 -0.0035 1 0; ? ?11 1 -0.0045 -0.003 1 0; ? ?12 1 -0.006 -0.0035 1 0; ? ?13 1 -0.006 -0.0035 1 0; ? ?14 1 -0.012 -0.008 1 0; ? ?15 1 -0.006 -0.001 1 0; ? ?16 1 -0.006 -0.002 1 0; ? ?17 1 -0.006 -0.002 1 0; ? ?18 1 -0.009 -0.004 1 0; ? ?19 1 -0.009 -0.004 1 0; ? ?20 1 -0.009 -0.004 1 0; ? ?21 1 -0.009 -0.004 1 0; ? ?22 1 -0.009 -0.004 1 0; ? ?23 1 -0.009 -0.005 1 0; ? ?24 1 -0.042 -0.02 1 0; ? ?25 1 -0.042 -0.02 1 0; ? ?26 1 -0.006 -0.0025 1 0; ? ?27 1 -0.006 -0.0025 1 0; ? ?28 1 -0.006 -0.002 1 0; ? ?29 1 -0.012 -0.007 1 0; ? ?30 1 -0.02 -0.06 1 0; ? ?31 1 -0.015 -0.007 1 0; ? ?32 1 -0.021 -0.01 1 0; ? ?33 1 -0.006 -0.004 1 0];for i=1:33 ? ?B2(i,3)= ?B2(i,3)*P_flex; %負荷時變系數 ? ?B2(i,4)= ?B2(i,4)*P_flex;endY=zeros(n); ? ? ? ? ? ? ? %zeros就是生成一個全0的矩陣Times=1; ? ? ? ? ? ? ? ? ?%置迭代次數為初始值global positionB2(position(1),3)=B2(position(1),3)-car_1/10000; ?%電動汽車接入B2(position(1),4)=B2(position(1),4)-0.484*car_1/10000; ?%功率因素為0.9B2(position(2),3)=B2(position(2),3)-car_2/10000; ?%電動汽車接入B2(position(2),4)=B2(position(2),4)-0.484*car_2/10000; ?%功率因素為0.9B2(position(3),3)=B2(position(3),3)-car_3/10000; ?%電動汽車接入B2(position(3),4)=B2(position(3),4)-0.484*car_3/10000; ?%功率因素為0.9for i=1:n1 ? ? ? ?p=B1(i,1); ? ? ? ?q=B1(i,2); ? ? ? ?Y(p,q)=Y(p,q)-1/((B1(i,3)+B1(i,4))*B1(i,5)); ? ? ? ?Y(q,p)=Y(p,q); ? ? ? ?Y(p,p)=Y(p,p)+1/(B1(i,3)+B1(i,4))+0.5*B1(i,6);% ? ? ? ?Y(q,q)=Y(q,q)+1/((B1(i,3)+B1(i,4))*B1(i,5)^2)+0.5*B1(i,6);%高壓側阻抗乘以變比平方 ?輸入時注意低壓側在前end%disp('節(jié)點導納矩陣：') ;Y;G=real(Y);B=imag(Y);OrgS=zeros(2*n-2,1);DetaS=zeros(2*n-2,1); ? %將OrgS、DetaS初始化%創(chuàng)建OrgS，用于存儲初始功率參數Q=0;PQV=0;x=1.655; ? ? ? ? %%%x1=6.7,x2=9.85,x=x1+x2；其中x1為定子漏抗，x2為轉子漏抗xp=18.8; ? ? ? ? %%%xc=,xm=,xp=xc*xm/(xc-xm)；其中xc為機端并聯電容器電抗，xm為激磁電抗h=0;for i=1:n ? ?%對PQ(V)節(jié)點的處理 ? ?h=h+1; ? ?if i~=isb&&B2(i,2)==3 ? ? ? ? ? Q(i)=-(B2(i,5))^2/xp+(-(B2(i,5))^2+sqrt((B2(i,5))^4-4*(B2(i,3))^2*x^2))/2*x; ? ? ? ? B2(i,4)=Q(i); ? ? ? ? B2(i,2)=1; ? ? ? ? PQV=h; ? ?end ?end ? ? ? %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Ig=0.01;Q=0;PI=0;h=0;for i=1:n ? ?%對PI節(jié)點的處理 ? ?h=h+1; ? ?if i~=isb&&B2(i,2)==4 ? ? ? ? ? Q(i)=sqrt(Ig^2*((B2(i,5))^2)-B(i,3)^2); ? ? ? ?%e=B2(i,5)，f=0，e^2+f^2=B2(i,5))^2，其中e和f為光伏發(fā)電系統接入節(jié)點電壓的實部和虛部 ? ? ? ? B2(i,4)=Q(i); ? ? ? ? B2(i,2)=1; ? ? ? ? PI=h; ? ?end ?end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%h=0;j=0;for i=1:n ? ? ? ? ? ?%對PQ節(jié)點的處理 ? ?if i~=isb&&B2(i,2)==1 ? ? ? ? ?h=h+1; ? ? ? ?for j=1:n ? ? ? ? ? ?OrgS(2*h-1,1)=OrgS(2*h-1,1)+B2(i,5)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))+B2(i,6)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5)); ? %Pi 書P57頁11-45 ? ? ? ? ? ?OrgS(2*h,1)=OrgS(2*h,1)+B2(i,6)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))-B2(i,5)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5)); ? ? ? %Qi ?同上 ? ? ? ?end ? ?endendfor i=1:n ? ? ? ? ? %對PV節(jié)點的處理 ? ?if i~=isb&&B2(i,2)==2 ? ? ? ?h=h+1; ? ? ? ?for j=1:n ? ? ? ? ? ?OrgS(2*h-1,1)=OrgS(2*h-1,1)+ B2(i,5)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))+B2(i,6)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5)); %同上 ? ? ? ? ? ?OrgS(2*h,1)=OrgS(2*h,1)+ B2(i,6)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))-B2(i,5)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5)); ? ?%同上 ? ? ? ?end ? ?endendOrgS;%創(chuàng)建PVU 用于存儲PV節(jié)點的初始電壓PVU=zeros(n-H-1,2);t=0;for i=1:n ? ?if B2(i,2)==2 ? ? ? ?t=t+1; ? ? ? ?PVU(t,1)=B2(i,5); ? ? ? ?PVU(t,2)=B2(i,6); ? ?endend%disp('PV節(jié)點初始值:電壓、相位ei、fi：') ? ; PVU;%創(chuàng)建DetaS，用于存儲有功功率、無功功率和電壓幅值的不平衡量 書p58 11-46和11-47h=0;for i=1:n ? ? ? ? ? %對PQ節(jié)點的處理 ? ?if i~=isb&&B2(i,2)==1 ? ? ? ?h=h+1; ? ? ? ?DetaS(2*h-1,1)=B2(i,3)-OrgS(2*h-1,1); ? ? ? ?DetaS(2*h,1)=B2(i,4)-OrgS(2*h,1); ? ?endendt=0;for i=1:n ? ? ? ? ? %對PV節(jié)點的處理 ? ?if i~=isb&&B2(i,2)==2 ? ? ? ?h=h+1; ? ? ? ?t=t+1; ? ? ? ?DetaS(2*h-1,1)=B2(i,3)-OrgS(2*h-1,1); ? ? ? ?DetaS(2*h,1)=PVU(t,1)^2+PVU(t,2)^2-B2(i,5)^2-B2(i,6)^2; ? ?endend%disp('P、Q、V不平衡量：') ?;DetaS;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%創(chuàng)建I，用于存儲節(jié)點電流參數（計算雅克比矩陣用 具體見紙板公式推導）I=zeros(n-1,1);h=0;for i=1:n ? if i~=isb ? ? ? ?h=h+1; ? ? ? ?I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,5)+B2(i,6)*sqrt(-1)); ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?end ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? endI;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%創(chuàng)建Jacbi(雅可比矩陣) ? ? 各元素是書中各元素求法乘個-1 所以△W=J△VJacbi=zeros(2*n-2);h=0;k=0;for i=1:n ? ? ? %對PQ節(jié)點的處理 ? ?if B2(i,2)==1 ? ? ? ?h=h+1; ? ? ? ?for j=1:n ? ? ? ? ? ?if j~=isb ? ? ? ? ? ? ? ?k=k+1; ? ? ? ? ? ? ? ?if i==j ? ? %對角元素的處理 ? ? ? ? ? ? ? ? ? ?Jacbi(2*h-1,2*k-1)=-B(i,j)*B2(i,5)+G(i,j)*B2(i,6)+imag(I(h,1)); ? ? ? ? ? ? ? ? ? ? ? ? ?Jacbi(2*h-1,2*k)=G(i,j)*B2(i,5)+B(i,j)*B2(i,6)+real(I(h,1)); ? ? ? ? ? ? ? ? ? ?Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1)); ? ? ? ? ? ? ? ? ? ?Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1)); ? ? ? ? ? ? ? ?else ? ? ? ?%非對角元素的處理 ? ? ? ? ? ? ? ? ? ?Jacbi(2*h-1,2*k-1)=-B(i,j)*B2(i,5)+G(i,j)*B2(i,6); ? ? ? ? ? ? ? ? ? ?Jacbi(2*h-1,2*k)=G(i,j)*B2(i,5)+B(i,j)*B2(i,6); ? ? ? ? ? ? ? ? ? ?Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k); ? ? ? ? ? ? ? ? ? ?Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1); ? ? ? ? ? ? ? ?end ? ? ? ? ? ? ? ?if k==(n-1) %將用于內循環(huán)的指針置于初始值，以確保雅可比矩陣換行 ? ? ? ? ? ? ? ? ? ?k=0; ? ? ? ? ? ? ? ?end ? ? ? ? ? ?end ? ? ? ?end ? ?endendk=0;for i=1:n ? ? ? %對PV節(jié)點的處理 ? ?if B2(i,2)==2 ? ? ? ?h=h+1; ? ? ? ?for j=1:n ? ? ? ? ? ?if j~=isb ? ? ? ? ? ? ? ?k=k+1; ? ? ? ? ? ? ? ?if i==j ? ? %對角元素的處理 ? ? ? ? ? ? ? ? ? ?Jacbi(2*h-1,2*k-1)= -B(i,j)*B2(i,5)+G(i,j)*B2(i,6)+imag(I(h,1)); ? ? ? ? ? ? ? ? ? ?Jacbi(2*h-1,2*k)= G(i,j)*B2(i,5)+B(i,j)*B2(i,6)+real(I(h,1)); ? ? ? ? ? ? ? ? ? ?Jacbi(2*h,2*k-1)=2*B2(i,6); ? ? ? ? ? ? ? ? ? ?Jacbi(2*h,2*k)=2*B2(i,5); ? ? ? ? ? ? ? ?else ? ? ? ?%非對角元素的處理 ? ? ? ? ? ? ? ? ? ?Jacbi(2*h-1,2*k-1)= -B(i,j)*B2(i,5)+G(i,j)*B2(i,6); ? ? ? ? ? ? ? ? ? ?Jacbi(2*h-1,2*k)= G(i,j)*B2(i,5)+B(i,j)*B2(i,6); ? ? ? ? ? ? ? ? ? ?Jacbi(2*h,2*k-1)=0; ? ? ? ? ? ? ? ? ? ?Jacbi(2*h,2*k)=0; ? ? ? ? ? ? ? ?end ? ? ? ? ? ? ? ?if k==(n-1) ? ? %將用于內循環(huán)的指針置于初始值，以確保雅可比矩陣換行 ? ? ? ? ? ? ? ? ? ?k=0; ? ? ? ? ? ? ? ?end ? ? ? ? ? ?end ? ? ? ?end ? ?endend%disp('雅克比矩陣為：') ?;Jacbi; %列參數跟書上的相反 ? 這里面第一列是對電壓相位的偏導 第二列是對幅值的偏導%求解修正方程，獲取節(jié)點電壓的不平衡量DetaU=zeros(2*n-2,1);DetaU=Jacbi\DetaS; ? ? ? %disp('求解修正方程得：') ?;DetaU;%修正節(jié)點電壓j=0;for i=1:n ? ? ? %對PQ節(jié)點處理 ? ?if B2(i,2)==1 ? ? ? ?j=j+1; ? ? ? ?B2(i,5)=B2(i,5)+DetaU(2*j,1); ? ? ? ?B2(i,6)=B2(i,6)+DetaU(2*j-1,1); ? ?endendfor i=1:n ? ? ? %對PV節(jié)點的處理 ? ?if B2(i,2)==2 ? ? ? ?j=j+1; ? ? ? B2(i,5)=B2(i,5)+DetaU(2*j,1); ? ? ? B2(i,6)=B2(i,6)+DetaU(2*j-1,1); ? ?endendfor i=1:n ? ? ? ? ? ?%對PQ(V)節(jié)點的處理 ? ?if i==PQV ? ? ? B2(i,2)=3; ? ?endendfor i=1:n ? ? ? ? ? ?%對PI節(jié)點的處理 ? ?if i==PI ? ? ? B2(i,2)=4; ? ?endend%disp('修正后的B2陣：') ?;B2;%開始循環(huán)**********************************************************************while max(abs(DetaU))>pr%先取絕對值在找出最大值來與誤差精度比較for i=1:n ? ? ? ? ? ?%對PQ(V)節(jié)點的處理 ? ?if i~=isb&&B2(i,2)==3 ? ? ? ? ? Q(i)=-(B2(i,5))^2/xp+(-(B2(i,5))^2+sqrt((B2(i,5))^4-4*(B2(i,3))^2*x^2))/2*x; ? ? ? ? B2(i,4)=Q(i); ? ?end ? ?if B2(i,2)==3 ? ? ? B2(i,2)=1; ? ?endendfor i=1:n ? ?%對PI節(jié)點的處理 ? ?if i~=isb&&B2(i,2)==4 ? ? ? ? ? Q(i)=sqrt(Ig^2*((B2(i,5))^2)-B(i,3)^2); ? ? ? ? ? ? ? ?B2(i,4)=Q(i); ? ?end ? ?if B2(i,2)==4 ? ? ? B2(i,2)=1; ? ?endendOrgS=zeros(2*n-2,1); ? ? ? ?h=0;j=0;for i=1:n ? ? ? ? ? ?%對PQ節(jié)點的處理 ? ?if i~=isb&&B2(i,2)==1 ? ? ? ? ?h=h+1; ? ? ? ?for j=1:n ? ? ? ? ? ?OrgS(2*h-1,1)=OrgS(2*h-1,1)+B2(i,5)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))+B2(i,6)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5)); ? ? ? ? %Pi ? ? ? ? ? ?OrgS(2*h,1)=OrgS(2*h,1)+B2(i,6)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))-B2(i,5)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5)); ? ? ? ? ? ?%Qi ? ? ? ?end ? ?endendfor i=1:n ? ? ? ? ? %對PV節(jié)點的處理 ? ?if i~=isb&&B2(i,2)==2 ? ? ? ?h=h+1; ? ? ? ?for j=1:n ? ? ? ? ? ?OrgS(2*h-1,1)=OrgS(2*h-1,1)+ B2(i,5)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))+B2(i,6)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5)); ? ? ? ? ? ?OrgS(2*h,1)=OrgS(2*h,1)+ B2(i,6)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))-B2(i,5)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5)); ? ? ? ?end ? ?endend%disp('修正后的迭代計算PQ、PV節(jié)點參數：') ; OrgS;%創(chuàng)建DetaSh=0;for i=1:n ? ? ? ? ? %對PQ節(jié)點的處理 ? ?if i~=isb&&B2(i,2)==1 ? ? ? ?h=h+1; ? ? ? ?DetaS(2*h-1,1)=B2(i,3)-OrgS(2*h-1,1); ? ? ? ?DetaS(2*h,1)=B2(i,4)-OrgS(2*h,1); ? ?endendt=0;for i=1:n ? ? ? ? ? %對PV節(jié)點的處理 ? ?if i~=isb&&B2(i,2)==2 ? ? ? ?h=h+1; ? ? ? ?t=t+1; ? ? ? ?DetaS(2*h-1,1)=B2(i,3)-OrgS(2*h-1,1); ? ? ? ?DetaS(2*h,1)=PVU(t,1)^2+PVU(t,2)^2-B2(i,5)^2-B2(i,6)^2; ? ?endend%disp('修正后的迭代計算PQ、PV節(jié)點不平衡量：') ; DetaS;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%創(chuàng)建II=zeros(n-1,1);h=0;for i=1:n ? ?if i~=isb ? ? ? ?h=h+1; ? ? ? ?I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,5)+B2(i,6)*sqrt(-1)); ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?end ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? endI;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%創(chuàng)建JacbiJacbi=zeros(2*n-2);h=0;k=0;for i=1:n ? ? ? %對PQ節(jié)點的處理 ? ?if B2(i,2)==1 ? ? ? ?h=h+1; ? ? ? ?for j=1:n ? ? ? ? ? ?if j~=isb ? ? ? ? ? ? ? ?k=k+1; ? ? ? ? ? ? ? ?if i==j ? ? %對角元素的處理 ? ? ? ? ? ? ? ? ? ?Jacbi(2*h-1,2*k-1)=-B(i,j)*B2(i,5)+G(i,j)*B2(i,6)+imag(I(h,1)); ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?Jacbi(2*h-1,2*k)=G(i,j)*B2(i,5)+B(i,j)*B2(i,6)+real(I(h,1)); ? ? ? ? ? ? ? ? ? ?Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1)); ? ? ? ? ? ? ? ? ? ?Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1)); ? ? ? ? ? ? ? ?else ? ? ? ?%非對角元素的處理 ? ? ? ? ? ? ? ? ? ?Jacbi(2*h-1,2*k-1)=-B(i,j)*B2(i,5)+G(i,j)*B2(i,6); ? ? ? ? ? ? ? ? ? ?Jacbi(2*h-1,2*k)=G(i,j)*B2(i,5)+B(i,j)*B2(i,6); ? ? ? ? ? ? ? ? ? ?Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k); ? ? ? ? ? ? ? ? ? ?Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1); ? ? ? ? ? ? ? ?end ? ? ? ? ? ? ? ?if k==(n-1) %將用于內循環(huán)的指針置于初始值，以確保雅可比矩陣換行 ? ? ? ? ? ? ? ? ? ?k=0; ? ? ? ? ? ? ? ?end ? ? ? ? ? ?end ? ? ? ?end ? ?endendk=0;for i=1:n ? ? ? %對PV節(jié)點的處理 ? ?if B2(i,2)==2 ? ? ? ?h=h+1; ? ? ? ?for j=1:n ? ? ? ? ? ?if j~=isb ? ? ? ? ? ? ? ?k=k+1; ? ? ? ? ? ? ? ?if i==j ? ? %對角元素的處理 ? ? ? ? ? ? ? ? ? ?Jacbi(2*h-1,2*k-1)= -B(i,j)*B2(i,5)+G(i,j)*B2(i,6)+imag(I(h,1)); ? ? ? ? ? ? ? ? ? ?Jacbi(2*h-1,2*k)= G(i,j)*B2(i,5)+B(i,j)*B2(i,6)+real(I(h,1)); ? ? ? ? ? ? ? ? ? ?Jacbi(2*h,2*k-1)=2*B2(i,6); ? ? ? ? ? ? ? ? ? ?Jacbi(2*h,2*k)=2*B2(i,5); ? ? ? ? ? ? ? ?else ? ? ? ?%非對角元素的處理 ? ? ? ? ? ? ? ? ? ?Jacbi(2*h-1,2*k-1)= -B(i,j)*B2(i,5)+G(i,j)*B2(i,6); ? ? ? ? ? ? ? ? ? ?Jacbi(2*h-1,2*k)= G(i,j)*B2(i,5)+B(i,j)*B2(i,6); ? ? ? ? ? ? ? ? ? ?Jacbi(2*h,2*k-1)=0; ? ? ? ? ? ? ? ? ? ?Jacbi(2*h,2*k)=0; ? ? ? ? ? ? ? ?end ? ? ? ? ? ? ? ?if k==(n-1) ? ? %將用于內循環(huán)的指針置于初始值，以確保雅可比矩陣換行 ? ? ? ? ? ? ? ? ? ?k=0; ? ? ? ? ? ? ? ?end ? ? ? ? ? ?end ? ? ? ?end ? ?endend%disp('修正后的雅克比矩陣：') ?;Jacbi;DetaU=zeros(2*n-2,1);DetaU=Jacbi\DetaS;DetaU;%修正節(jié)點電壓j=0;for i=1:n ? ? ? %對PQ節(jié)點處理 ? ?if B2(i,2)==1 ? ? ? ?j=j+1; ? ? ? ?B2(i,5)=B2(i,5)+DetaU(2*j,1); ? ? ? ?B2(i,6)=B2(i,6)+DetaU(2*j-1,1); ? ?endendfor i=1:n ? ? ? %對PV節(jié)點的處理 ? ?if B2(i,2)==2 ? ? ? ?j=j+1; ? ? ? B2(i,5)=B2(i,5)+DetaU(2*j,1); ? ? ? B2(i,6)=B2(i,6)+DetaU(2*j-1,1); ? ?endendfor i=1:n ? ? ? ? ? ?%對PQ(V)節(jié)點的處理 ? ?if i==PQV ? ? ? B2(i,2)=3; ? ?endendfor i=1:n ? ? ? ? ? ?%對PI節(jié)點的處理 ? ?if i==PI ? ? ? B2(i,2)=4; ? ?endendTimes=Times+1; ? ? ?%迭代次數加1enddisp('初始潮流計算結果如下：');disp('迭代次數：');Timesdisp('節(jié)點電壓幅值如下（按節(jié)點由小到大的順序）：');B2(:,5)%創(chuàng)建Sb，用于存儲平衡節(jié)點功率Sb=zeros(1);for i=1:n ? ?if i==isb ? ? ? ?for j=1:n ? ? Sb=Sb+(B2(i,5)+sqrt(-1)*B2(i,6))*conj(Y(i,j))*conj(B2(j,5)+sqrt(-1)*B2(j,6)); ? ? ? ?end ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? endenddisp('初始平衡節(jié)點功率：')Sb%求解各條支路的功率Sij=zeros(n);for i=1:n ? ? for j=1:n ? ? Sij(i,j)=Sij(i,j)+(B2(i,5)+sqrt(-1)*B2(i,6))^2*conj(B1(1,4))+(B2(i,5)+sqrt(-1)*B2(i,6))*conj(Y(i,j))*(conj((B2(i,5)+sqrt(-1)*B2(i,6)))-conj((B2(j,5)+sqrt(-1)*B2(j,6)))); ? ? ?endend%disp('線路功率分布：')Sij;Sij_P=real(Sij);Sij_Q=imag(Sij);for i=1:32 ? ? P_flow(i)= ?abs( Sij_P(B1(i,1),B1(i,2)));endfor i=1:33 ? ?for j=1:33 ? ? ? ?if Sij_P(i,j)<0.0001 ? ? ? ? ? ? Sij_P(i,j)=0; ? ? ? ?end ? ? ? ? if Sij_Q(i,j)<0.0001 ? ? ? ? ? ? Sij_Q(i,j)=0; ? ? ? ? end ? ?endend%%%求解系統網損Ploss=0;for i=1:n ? ? ? ? ? ? ? ? ? ?for j=1:n ? ? ? ? ? ?Ploss=Ploss+B2(i,5)*B2(j,5)*(G(i,j)*cos(B2(i,6)-B2(j,6))+B(i,j)*sin(B2(i,6)-B2(j,6))); ? ? ? ?endendPloss=Ploss*10000; V_amp=B2(:,5); v_d=0;for i=1:n ?%電壓罰函數 ? ?if V_amp<0.9 ? ? ? ?v_d=v_d+abs(V_amp-0.9); ? ?end ? ?if V_amp>1.1 ? ? ? ?v_d=v_d+abs(V_amp-1.1); ? ?end ? ?V_33(i)=abs(1-V_amp(i)); %求與額定電壓的差值endPloss=Ploss+5*v_d; %目標函數加懲罰項V_flu=sum(V_33);end
?? 運行結果
