用matlab電力系統(tǒng)潮流計算

1、題目：潮流計算與matlab教學單位 電氣信息學院 姓 名 學 號 年 級 專 業(yè) 電氣工程及其自動化 指導教師 職 稱 副教授 摘 要電力系統(tǒng)穩(wěn)態(tài)分析包括潮流計算和靜態(tài)安全分析。本文主要運用的事潮流計算，潮流計算是電力網絡設計與運行中最基本的運算，對電力網絡的各種設計方案及各種運行方式進行潮流計算，可以得到各種電網各節(jié)點的電壓，并求得網絡的潮流及網絡中的各元件的電力損耗，進而求得電能損耗。本位就是運用潮流計算具體分析，并有matlab仿真。關鍵詞： 電力系統(tǒng) 潮流計算 matlababstractelectric power system steady flow calculation an

2、d analysis of the static safety analysis. this paper, by means of the calculation, flow calculation is the trend of the power network design and operation of the most basic operations of electric power network, various design scheme and the operation ways to tide computation, can get all kinds of ea

3、ch node of the power grid voltage and seek the trend of the network and the network of the components of the power loss, and getting electric power. the standard is to use the power flow calculation and analysis, the specific have matlab simulation.key words: power system; flow calculation; matlab s

4、imulation32 目 錄1 任務提出與方案論證22 總體設計32.1潮流計算等值電路32.2建立電力系統(tǒng)模型32.3模型的調試與運行33 詳細設計43.1 計算前提43.2手工計算74設計圖及源程序114.1matlab仿真114.2潮流計算源程序115 總結19參考文獻20 1 任務提出與方案論證潮流計算是在給定電力系統(tǒng)網絡結構、參數和決定系統(tǒng)運行狀態(tài)的邊界條件的情況下確定系統(tǒng)穩(wěn)態(tài)運行狀態(tài)的一種基本方法,是電力系統(tǒng)規(guī)劃和運營中不可缺少的一個重要組成部分。可以說,它是電力系統(tǒng)分析中最基本、 最重要的計算,是系統(tǒng)安全、 經濟分析和實時控制與調度的基礎。常規(guī)潮流計算的任務是根據給定的運行條件

5、和網路結構確定整個系統(tǒng)的運行狀態(tài)，如各母線上的電壓（幅值及相角）、網絡中的功率分布以及功率損耗等。潮流計算的結果是電力系統(tǒng)穩(wěn)定計算和故障分析的基礎。在電力系統(tǒng)運行方式和規(guī)劃方案的研究中，都需要進行潮流計算以比較運行方式或規(guī)劃供電方案的可行性、可靠性和經濟性。同時，為了實時監(jiān)控電力系統(tǒng)的運行狀態(tài)，也需要進行大量而快速的潮流計算。因此，潮流計算是電力系統(tǒng)中應用最廣泛、最基本和最重要的一種電氣運算。在系統(tǒng)規(guī)劃設計和安排系統(tǒng)的運行方式時，采用離線潮流計算；在電力系統(tǒng)運行狀態(tài)的實時監(jiān)控中，則采用在線潮流計算。是電力系統(tǒng)研究人員長期研究的一個課題。它既是對電力系統(tǒng)規(guī)劃設計和運行方式的合理性、 可靠性及經濟

6、性進行定量分析的依據 ， 又是電力系統(tǒng)靜態(tài)和暫態(tài)穩(wěn)定計算的基礎。潮流計算經歷了一個由手工到應用數字電子計算機的發(fā)展過程，現在的潮流算法都以計算機的應用為前提用計算機進行潮流計算主要步驟在于編制計算機程序，這是一項非常復雜的工作。對系統(tǒng)進行潮流分析，本文利用 matlab中的simpowersystems工具箱設計電力系統(tǒng)，在simulink 環(huán)境下，不僅可以仿真系統(tǒng)的動態(tài)過程，還可以對系統(tǒng)進行穩(wěn)態(tài)潮流分析。2 總體設計simpowersystems使用simulink環(huán)境，可以將該系統(tǒng)中的發(fā)電機、變壓器，線路等模型聯(lián)結起來，形成電力系統(tǒng)仿真模擬圖。在加人測量模塊，并對各元件的參數進行設置后，用

7、measurement和sink中的儀器可以觀察各元件的電壓、電流、功率的大小。2.1潮流計算等值電路 2.2建立電力系統(tǒng)模型在simulink中按照電力系統(tǒng)原型選擇元件進行建模。所建立的模型和建立的方法在詳細設計中詳述。在電力系統(tǒng)模型的建立工程中主要涉及到的是：元器件的選擇及其參數的設置；發(fā)電機選型；變壓器選擇；線路的選擇；負荷模型的選擇；母線選擇。 2.3模型的調試與運行建立系統(tǒng)模型，并設置好參數以后，就可以在simulink環(huán)境下進行仿真運行。運行的具體結果和分析也在詳細設計中詳述。3 詳細設計3.1 計算前提首先是發(fā)電機的參數計算，先對5個發(fā)電廠簡化為5臺發(fā)電機來計算。發(fā)電機g1：發(fā)電

8、機g2：發(fā)電機g3：發(fā)電機g4：發(fā)電機g5：其次是變電站的參數計算，我們還是對7個變電站簡化為7臺變壓器來計算。變壓器t1：變壓器t2：（雙并聯(lián)）變壓器t3：（四并聯(lián)）變壓器t4：（雙并聯(lián)）變壓器t5：變壓器t6：（兩個三繞組變壓器并聯(lián)）變壓器t7：（雙并聯(lián)）再次是傳輸線參數計算，5條傳輸線的具體計算如下。根據教材查得線路l1：線路l2：線路l3：（雙回路）線路l4：線路l5：（雙回路）3.2手工計算flr1：flr2：flr3：flr4：flr5:計算每一個flr的功率分布和電壓分布計算如下：flr1：flr2：功率分布： 電壓分布：flr3：功率分布： 電壓分布：flr4：功率分布： 電壓

9、分布：flr5:這里我們先將f點和發(fā)電機g5當做電源，經過和構成兩端供電網絡以g點作為運算負荷進行計算。電壓分布：4設計圖及源程序4.1matlab仿真相關的原始數據輸入格式如下：1、b1是支路參數矩陣，第一列和第二列是節(jié)點編號。節(jié)點編號由小到大編寫。2、對于含有變壓器的支路，第一列為低壓側節(jié)點編號，第二列為高壓側節(jié)點編號，將變壓器的串聯(lián)阻抗置于低壓側處理，第三列為支路的串列阻抗參數，第四列為支路的對地導納參數，第五烈為含變壓器支路的變壓器的變比，第六列為變壓器是否是否含有變壓器的參數，其中“1”為含有變壓器，“0”為不含有變壓器。3、b2為節(jié)點參數矩陣，其中第一列為節(jié)點注入發(fā)電功率參數；第二

10、列為節(jié)點負荷功率參數；第三列為節(jié)點電壓參數；第六列為節(jié)點類型參數，其中“1”為平衡節(jié)點，“2”為pq節(jié)點，“3”為pv節(jié)點參數。4、x為節(jié)點號和對地參數矩陣。其中第一列為節(jié)點編號，第二列為節(jié)點對地參數。4.2潮流計算源程序%本程序的功能是用牛頓拉夫遜法進行11節(jié)點潮流計算clear;n=11;%input(請輸入節(jié)點數：n=);nl=11;%input(請輸入支路數：nl=);isb=1;%input(請輸入平衡母線節(jié)點號：isb=);pr=0.00001;%input(請輸入誤差精度：pr=);b1=1 2 0.03512+0.08306i 0.13455i 1 0; 2 3 0.0068+

11、0.18375i 0 1.02381 1; 1 4 0.05620+0.13289i 0.05382i 1 0; 4 5 0.00811+0.24549i 0 1.02381 1; 1 6 0.05620+0.13289i 0.05382i 1 0; 4 6 0.04215+0.09967i 0.04037i 1 0; 6 7 0.0068+0.18375i 0 1.02381 1; 6 8 0.02810+0.06645i 0.10764i 1 0; 1 4 0.05620+0.13289i 0.05382i 1 0; 8 10 0.00811+0.24549i 0 1 1; 8 9 0.0

12、3512+0.08306i 0.13455i 1 0b2=0 0 1.1 1.1 0 1; 0 0 1 0 0 2; 0 0.343+0.21256i 1 0 0 2; 0 0 1 0 0 2; 0 0.204+0.12638i 1 0 0 2; 0 0 1 0 0 2; 0 0.306+0.18962i 1 0 0 2; 0 0 1 0 0 2; 0.5 0 1.1 1.1 0 3; 0 0.343+0.21256i 1 0 0 2; 0 0 1 0 0 2;% b1矩陣：1、支路首端號；2、末端號；3、支路阻抗；4、支路對地電納% 5、支路的變比；6、支路首端處于k側為1，1側為0% b2

13、矩陣：1、該節(jié)點發(fā)電機功率；2、該節(jié)點負荷功率；3、節(jié)點電壓初始值% 4、pv節(jié)點電壓v的給定值；5、節(jié)點所接的無功補償設備的容量% 6、節(jié)點分類標號：1為平衡節(jié)點（應為1號節(jié)點）；2為pq節(jié)點；% 3為pv節(jié)點；%input(請輸入各節(jié)點參數形成的矩陣： b2=);y=zeros(n);e=zeros(1,n);f=zeros(1,n);v=zeros(1,n);sida=zeros(1,n);s1=zeros(nl);% % %-for i=1:nl%支路數 if b1(i,6)=0%左節(jié)點處于1側 p=b1(i,1);q=b1(i,2); else %左節(jié)點處于k側 p=b1(i,2);

14、q=b1(i,1); end y(p,q)=y(p,q)-1./(b1(i,3)*b1(i,5);%非對角元 y(q,p)=y(p,q); %非對角元 y(q,q)=y(q,q)+1./(b1(i,3)*b1(i,5)2)+b1(i,4)./2;%對角元k側 y(p,p)=y(p,p)+1./b1(i,3)+b1(i,4)./2;%對角元1側 end%求導納矩陣disp(導納矩陣 y=);disp(y)%-g=real(y);b=imag(y);%分解出導納陣的實部和虛部 for i=1:n%給定各節(jié)點初始電壓的實部和虛部 e(i)=real(b2(i,3); f(i)=imag(b2(i,3

15、); v(i)=b2(i,4);%pv節(jié)點電壓給定模值 endfor i=1:n%給定各節(jié)點注入功率 s(i)=b2(i,1)-b2(i,2); %i節(jié)點注入功率sg-sl b(i,i)=b(i,i)+b2(i,5);%i節(jié)點無功補償量 end%=p=real(s);q=imag(s); %分解出各節(jié)點注入的有功和無功功率ict1=0;it2=1;n0=2*n;n=n0+1;a=0; %迭代次數ict1、a；不滿足收斂要求的節(jié)點數it2while it2=0 % n0=2*n 雅可比矩陣的階數；n=n0+1擴展列 it2=0;a=a+1; for i=1:n if i=isb%非平衡節(jié)點 c(

16、i)=0;d(i)=0; for j1=1:n c(i)=c(i)+g(i,j1)*e(j1)-b(i,j1)*f(j1);%(gij*ej-bij*fj) d(i)=d(i)+g(i,j1)*f(j1)+b(i,j1)*e(j1);%(gij*fj+bij*ej) end p1=c(i)*e(i)+f(i)*d(i);%節(jié)點功率p計算ei(gij*ej-bij*fj)+fi(gij*fj+bij*ej) q1=c(i)*f(i)-e(i)*d(i);%節(jié)點功率q計算fi(gij*ej-bij*fj)-ei(gij*fj+bij*ej)%求i節(jié)點有功和無功功率p,q的計算值 v2=e(i)2+

17、f(i)2;%電壓模平方 %= 以下針對非pv節(jié)點來求取功率差及jacobi矩陣元素 = if b2(i,6)=3%非pv節(jié)點 dp=p(i)-p1;%節(jié)點有功功率差 dq=q(i)-q1; %節(jié)點無功功率差 %= 以上為除平衡節(jié)點外其它節(jié)點的功率計算 =%= 求取jacobi矩陣 = for j1=1:n if j1=isb&j1=i%非平衡節(jié)點&非對角元 x1=-g(i,j1)*e(i)-b(i,j1)*f(i);% dp/de=-dq/df x2=b(i,j1)*e(i)-g(i,j1)*f(i);% dp/df=dq/de x3=x2; % x2=dp/df x3=dq/de x4=-

18、x1; % x1=dp/de x4=dq/df p=2*i-1;q=2*j1-1; j(p,q)=x3;j(p,n)=dq;m=p+1; % x3=dq/de j(p,n)=dq節(jié)點無功功率差 j(m,q)=x1;j(m,n)=dp;q=q+1; % x1=dp/de j(m,n)=dp節(jié)點有功功率差 j(p,q)=x4;j(m,q)=x2; % x4=dq/df x2=dp/df elseif j1=i&j1=isb%非平衡節(jié)點&對角元 x1=-c(i)-g(i,i)*e(i)-b(i,i)*f(i);% dp/de x2=-d(i)+b(i,i)*e(i)-g(i,i)*f(i);% dp

19、/df x3=d(i)+b(i,i)*e(i)-g(i,i)*f(i); % dq/de x4=-c(i)+g(i,i)*e(i)+b(i,i)*f(i);% dq/df p=2*i-1;q=2*j1-1;j(p,q)=x3;j(p,n)=dq;%擴展列q m=p+1; j(m,q)=x1;q=q+1;j(p,q)=x4;j(m,n)=dp;%擴展列p j(m,q)=x2; end end else %= 下面是針對pv節(jié)點來求取jacobi矩陣的元素 = dp=p(i)-p1;% pv節(jié)點有功誤差 dv=v(i)2-v2;% pv節(jié)點電壓誤差 for j1=1:n if j1=isb&j1=

20、i%非平衡節(jié)點&非對角元 x1=-g(i,j1)*e(i)-b(i,j1)*f(i); % dp/de x2=b(i,j1)*e(i)-g(i,j1)*f(i); % dp/df x5=0;x6=0; p=2*i-1;q=2*j1-1;j(p,q)=x5;j(p,n)=dv; % pv節(jié)點電壓誤差 m=p+1; j(m,q)=x1;j(m,n)=dp;q=q+1;j(p,q)=x6; % pv節(jié)點有功誤差 j(m,q)=x2; elseif j1=i&j1=isb %非平衡節(jié)點&對角元 x1=-c(i)-g(i,i)*e(i)-b(i,i)*f(i);% dp/de x2=-d(i)+b(i,

21、i)*e(i)-g(i,i)*f(i);% dp/df x5=-2*e(i); x6=-2*f(i); p=2*i-1;q=2*j1-1;j(p,q)=x5;j(p,n)=dv; % pv節(jié)點電壓誤差 m=p+1; j(m,q)=x1;j(m,n)=dp;q=q+1;j(p,q)=x6; % pv節(jié)點有功誤差 j(m,q)=x2; end end end end end%= 以上為求雅可比矩陣的各個元素及擴展列的功率差或電壓差 = for k=3:n0 % n0=2*n （從第三行開始，第一、二行是平衡節(jié)點） k1=k+1;n1=n; % n=n0+1 即 n=2*n+1擴展列p、q 或 u

22、for k2=k1:n1% 從k+1列的jacobi元素到擴展列的p、q 或 u j(k,k2)=j(k,k2)./(j(k,k)+eps);% 用k行k列對角元素去除k行k列后的非對角元素進行規(guī)格化 end j(k,k)=1; % 對角元規(guī)格化k行k列對角元素賦1 %= 回代運算 = if k=3 % 不是第三行 k 3 k4=k-1; for k3=3:k4% 用k3行從第三行開始到當前行的前一行k4行消去 for k2=k1:n1% k3行后各行上三角元素 j(k3,k2)=j(k3,k2)-j(k3,k)*j(k,k2);%消去運算（當前行k列元素消為0） end %用當前行k2列元素

23、減去當前行k列元素乘以第k行k2列元素 j(k3,k)=0; %當前行第k列元素已消為0 end if k=n0 %若已到最后一行 break; end %= 前代運算 = for k3=k1:n0 % 從k+1行到2*n最后一行 for k2=k1:n1 % 從k+1列到擴展列消去k+1行后各行下三角元素 j(k3,k2)=j(k3,k2)-j(k3,k)*j(k,k2);%消去運算 end %用當前行k2列元素減去當前行k列元素乘以第k行k2列元素 j(k3,k)=0; %當前行第k列元素已消為0 end else %是第三行k=3 %= 第三行k=3的前代運算 = for k3=k1:n

24、0 %從第四行到2n行（最后一行） for k2=k1:n1 %從第四列到2n+1列（即擴展列） j(k3,k2)=j(k3,k2)-j(k3,k)*j(k,k2);%消去運算（當前行3列元素消為0） end %用當前行k2列元素減去當前行3列元素乘以第三行k2列元素 j(k3,k)=0; %當前行第3列元素已消為0 end end end%=上面是用線性變換方式高斯消去法將jacobi矩陣化成單位矩陣= for k=3:2:n0-1 l=(k+1)./2; e(l)=e(l)-j(k,n); %修改節(jié)點電壓實部 k1=k+1; f(l)=f(l)-j(k1,n); %修改節(jié)點電壓虛部 end

25、 %-修改節(jié)點電壓- for k=3:n0 det=abs(j(k,n); if det=pr %電壓偏差量是否滿足要求 it2=it2+1; %不滿足要求的節(jié)點數加1 end end ict2(a)=it2; %不滿足要求的節(jié)點數 ict1=ict1+1; %迭代次數end%用高斯消去法解w=-j*vdisp(迭代次數：);disp(ict1);disp(沒有達到精度要求的個數：);disp(ict2);for k=1:n v(k)=sqrt(e(k)2+f(k)2); %計算各節(jié)點電壓的模值 sida(k)=atan(f(k)./e(k)*180./pi; %計算各節(jié)點電壓的角度 e(k)

26、=e(k)+f(k)*j; %將各節(jié)點電壓用復數表示end%= 計算各輸出量 =disp(各節(jié)點的實際電壓標幺值e為(節(jié)點號從小到大排列)：);disp(e); %顯示各節(jié)點的實際電壓標幺值e用復數表示disp(-);disp(各節(jié)點的電壓大小v為(節(jié)點號從小到大排列)：);disp(v); %顯示各節(jié)點的電壓大小v的模值disp(-);disp(各節(jié)點的電壓相角sida為(節(jié)點號從小到大排列)：);disp(sida); %顯示各節(jié)點的電壓相角for p=1:n c(p)=0; for q=1:n c(p)=c(p)+conj(y(p,q)*conj(e(q); %計算各節(jié)點的注入電流的共軛

27、值 end s(p)=e(p)*c(p); %計算各節(jié)點的功率 s = 電壓 x 注入電流的共軛值enddisp(各節(jié)點的功率s為(節(jié)點號從小到大排列)：);disp(s); %顯示各節(jié)點的注入功率disp(-);disp(各條支路的首端功率si為(順序同您輸入b1時一致)：);for i=1:nl p=b1(i,1);q=b1(i,2); if b1(i,6)=0 si(p,q)=e(p)*(conj(e(p)*conj(b1(i,4)./2)+(conj(e(p)*b1(i,5). -conj(e(q)*conj(1./(b1(i,3)*b1(i,5); siz(i)=si(p,q); e

28、lse si(p,q)=e(p)*(conj(e(p)*conj(b1(i,4)./2)+(conj(e(p)./b1(i,5). -conj(e(q)*conj(1./(b1(i,3)*b1(i,5); siz(i)=si(p,q); end disp(si(p,q); ssi(p,q)=si(p,q); zf=s(,num2str(p),num2str(q),)=,num2str(ssi(p,q); disp(zf); disp(-);enddisp(各條支路的末端功率sj為(順序同您輸入b1時一致)：);for i=1:nl p=b1(i,1);q=b1(i,2); if b1(i,6)

29、=0 sj(q,p)=e(q)*(conj(e(q)*conj(b1(i,4)./2)+(conj(e(q)./b1(i,5). -conj(e(p)*conj(1./(b1(i,3)*b1(i,5); sjy(i)=sj(q,p); else sj(q,p)=e(q)*(conj(e(q)*conj(b1(i,4)./2)+(conj(e(q)*b1(i,5). -conj(e(p)*conj(1./(b1(i,3)*b1(i,5); sjy(i)=sj(q,p); end disp(sj(q,p); ssj(q,p)=sj(q,p); zf=s(,num2str(q),num2str(p)

30、,)=,num2str(ssj(q,p); disp(zf); disp(-);enddisp(各條支路的功率損耗ds為(順序同您輸入b1時一致)：);for i=1:nl p=b1(i,1);q=b1(i,2); ds(i)=si(p,q)+sj(q,p); disp(ds(i); dds(i)=ds(i); zf=ds(,num2str(p),num2str(q),)=,num2str(dds(i); disp(zf); disp(-);end%本程序的功能是用牛頓拉夫遜法進行10節(jié)點潮流計算%本程序的功能是用牛頓拉夫遜法進行潮流計算clear;n=10;%input(請輸入節(jié)點數：n=)

31、;nl=10;%input(請輸入支路數：nl=);isb=1;%input(請輸入平衡母線節(jié)點號：isb=);pr=0.00001;%input(請輸入誤差精度：pr=);b1=1 2 0.03512+0.08306i 0.13455i 1 0; 2 3 0.0068+0.18375i 0 1.02381 1; 1 4 0.05620+0.13289i 0.05382i 1 0; 4 5 0.00811+0.24549i 0 1.02381 1; 1 6 0.05620+0.13289i 0.05382i 1 0; 4 6 0.04215+0.09967i 0.04037i 1 0; 6 7

32、 0.0068+0.18375i 0 1.02381 1; 6 8 0.02810+0.06645i 0.10764i 1 0; 8 10 0.00811+0.24549i 0 1 1; 8 9 0.03512+0.08306i 0.13455i 1 0b2=0 0 1.1 1.1 0 1; 0 0 1 0 0 2; 0 0.343+0.21256i 1 0 0 2; 0 0 1 0 0 2; 0 0.204+0.12638i 1 0 0 2; 0 0 1 0 0 2; 0 0.306+0.18962i 1 0 0 2; 0 0 1 0 0 2; 0.5 0 1.1 1.1 0 3; 0 0.

33、343+0.21256i 1 0 0 2;% b1矩陣：1、支路首端號；2、末端號；3、支路阻抗；4、支路對地電納% 5、支路的變比；6、支路首端處于k側為1，1側為0% b2矩陣：1、該節(jié)點發(fā)電機功率；2、該節(jié)點負荷功率；3、節(jié)點電壓初始值% 4、pv節(jié)點電壓v的給定值；5、節(jié)點所接的無功補償設備的容量% 6、節(jié)點分類標號：1為平衡節(jié)點（應為1號節(jié)點）；2為pq節(jié)點；% 3為pv節(jié)點；%input(請輸入各節(jié)點參數形成的矩陣： b2=);y=zeros(n);e=zeros(1,n);f=zeros(1,n);v=zeros(1,n);sida=zeros(1,n);s1=zeros(nl)

34、;% % %-for i=1:nl%支路數 if b1(i,6)=0%左節(jié)點處于1側 p=b1(i,1);q=b1(i,2); else %左節(jié)點處于k側 p=b1(i,2);q=b1(i,1); end y(p,q)=y(p,q)-1./(b1(i,3)*b1(i,5);%非對角元 y(q,p)=y(p,q); %非對角元 y(q,q)=y(q,q)+1./(b1(i,3)*b1(i,5)2)+b1(i,4)./2;%對角元k側 y(p,p)=y(p,p)+1./b1(i,3)+b1(i,4)./2;%對角元1側 end%求導納矩陣disp(導納矩陣 y=);disp(y)%-g=real(

35、y);b=imag(y);%分解出導納陣的實部和虛部 for i=1:n%給定各節(jié)點初始電壓的實部和虛部 e(i)=real(b2(i,3); f(i)=imag(b2(i,3); v(i)=b2(i,4);%pv節(jié)點電壓給定模值 endfor i=1:n%給定各節(jié)點注入功率 s(i)=b2(i,1)-b2(i,2); %i節(jié)點注入功率sg-sl b(i,i)=b(i,i)+b2(i,5);%i節(jié)點無功補償量 end%=p=real(s);q=imag(s); %分解出各節(jié)點注入的有功和無功功率ict1=0;it2=1;n0=2*n;n=n0+1;a=0; %迭代次數ict1、a；不滿足收斂要

36、求的節(jié)點數it2while it2=0 % n0=2*n 雅可比矩陣的階數；n=n0+1擴展列 it2=0;a=a+1; for i=1:n if i=isb%非平衡節(jié)點 c(i)=0;d(i)=0; for j1=1:n c(i)=c(i)+g(i,j1)*e(j1)-b(i,j1)*f(j1);%(gij*ej-bij*fj) d(i)=d(i)+g(i,j1)*f(j1)+b(i,j1)*e(j1);%(gij*fj+bij*ej) end p1=c(i)*e(i)+f(i)*d(i);%節(jié)點功率p計算ei(gij*ej-bij*fj)+fi(gij*fj+bij*ej) q1=c(i)*f(i)-e(i)*d(i);%節(jié)點功率q計算fi(gij*ej-bij*fj)-ei(gij*fj+bij*ej)%求i節(jié)點有功和無功功率p,q的計算值 v2=e(i)2+f(i)2;%電壓模平方 %= 以下針對非pv節(jié)點來求取功率差及jacobi矩陣元素 = if b2(i,6)=3%非pv節(jié)點 dp=p(i)-p1;%節(jié)點有功功率差 dq=q(i)-q1; %節(jié)點無功功率差 %= 以上為除平衡節(jié)點外其