1、精選優(yōu)質文檔-傾情為你奉上電力系統(tǒng)短路故障分析的MATLAB輔助程序設計電力系統(tǒng)短路故障可分為三相對稱短路故障(three-phase balanced faults)和不對稱短路故障(unbalanced faults )。不對稱短路故障又分為單相接地短路故障(single line-to-ground fault)、兩相短路故障 (line-to-line fault)以及兩相接地短路故障(double line-to-ground fault)。根據(jù)故障分析結果可以對繼電保護裝置、自動裝置進行整定計算，我們可以建立算法來形成節(jié)點阻抗矩陣,利用節(jié)點阻抗矩陣來計算短路故障情況下的節(jié)點電壓和線

2、路電流。一、三相對稱短路故障進行三相短路計算需要兩個程序：zbuild / zbuildpi和symfault程序 ，zbuild、zbuildpi用來在MATLAB中形成節(jié)點阻抗矩陣，symfault用來計算三相對稱故障。Zbus=zbuild(zdata) 這里的參數(shù)zdata是一個(e×4)階矩陣，e是拓撲圖的總支路數(shù)目。第一列和第二列為元素兩端的節(jié)點編號，第三列和第四列分別是線路的電阻、電抗的標幺值。連接在0節(jié)點和發(fā)電機節(jié)點之間的發(fā)電機阻抗可能是次暫態(tài)電抗、暫態(tài)電抗或同步電抗，而且這個矩陣中還包含并聯(lián)電抗器和負荷阻抗。Zbus=zbuildpi(linedata,gendat

3、a,yload) 這個函數(shù)與潮流計算程序是相容的，第一個參數(shù)linedata與潮流計算程序中的文件是一致的。第一列和第二列為節(jié)點編號；第三列到第五列分別是線路的電阻、電抗以及1/2線路電納值，這三項都為在統(tǒng)一基準容量下的標幺值；最后一列是變壓器分接頭位置，對線路來說，必須輸入1；線路無輸入順序。發(fā)電機參數(shù)不包含在Linedata參數(shù)中，而是包含在第二個參數(shù)gendata中，gendata是一個g×4階矩陣,g是發(fā)電機總數(shù)。第一列和第二列為0節(jié)點、發(fā)電機節(jié)點編號，第三列和第四列為發(fā)電機的暫態(tài)電阻和暫態(tài)電抗。最后一個參數(shù)yload是可選擇的，這個矩陣有兩列，第一列為節(jié)點編號，第二列為復數(shù)

4、導納值，yload可以由潮流程序lfguss，lfnewton或者decouple自動生成。Zbuild和 zbuildpi兩個函數(shù)可以通過建立算法求出節(jié)點阻抗矩陣。首先添加所有與參考節(jié)點相連的樹支，然后添加其余的樹支，最后添加共軛連支。程序symfault(zdata,Zbus,V)用來計算三相對稱故障，程序要求輸入zdata和Zbus兩個矩陣，第三個參數(shù)V是可選的。如果V不存在，程序將默認故障前所有的節(jié)點電壓標幺值為1.0，如果變量V存在，那么V包括節(jié)點編號和復數(shù)電壓值。電壓向量V也可以由潮流計算程序自動生成。當symfault程序運行時，用戶要輸入故障節(jié)點編號和故障阻抗，運行可得到總的故

5、障電流，節(jié)點電壓幅值以及故障情況下的線路電流。在三相短路計算中，zbuild 和symfault程序，zbuildpi和symfault程序都可以進行計算，下面是三相短路計算使用的程序代碼：（1）Zbuild.m程序代碼：function Zbus = zbuild(linedata)nl = linedata(:,1); nr = linedata(:,2); R = linedata(:,3);X = linedata(:,4);nbr=length(linedata(:,1); nbus = max(max(nl), max(nr);for k=1:nbr if R(k) = inf |

6、 X(k) =inf R(k) = ; X(k) = ;%無窮 else, endend ZB = R + j*X;Zbus = zeros(nbus, nbus);tree=0; %從參考總線0上添加一個分支 for I = 1:nbr ntree(I) = 1; if nl(I) = 0 | nr(I) = 0 if nl(I) = 0 n = nr(I); elseif nr(I) = 0 n = nl(I); end if abs(Zbus(n, n) = 0 Zbus(n,n) = ZB(I);tree=tree+1; %new else Zbus(n,n) = Zbus(n,n)*

7、ZB(I)/(Zbus(n,n) + ZB(I); end ntree(I) = 2; else,end end% 添加一個新總線分支到現(xiàn)有總線上while tree < nbus for n = 1:nbus nadd = 1; if abs(Zbus(n,n) = 0 for I = 1:nbr if nadd = 1; if nl(I) = n | nr(I) = n if nl(I) = n k = nr(I); elseif nr(I) = n k = nl(I); end if abs(Zbus(k,k) = 0 for m = 1:nbus if m = n Zbus(m,

8、n) = Zbus(m,k); Zbus(n,m) = Zbus(m,k); else, end end Zbus(n,n) = Zbus(k,k) + ZB(I); tree=tree+1; nadd = 2; ntree(I) = 2; else, end else, end else, end end else, endend end %增加兩個原有總線間的支路阻抗 for n = 1:nbus for I = 1:nbr if ntree(I) = 1 if nl(I) = n | nr(I) = n if nl(I) = n k = nr(I); elseif nr(I) = n k

9、 = nl(I); end DM = Zbus(n,n) + Zbus(k,k) + ZB(I) - 2*Zbus(n,k); for jj = 1:nbus AP = Zbus(jj,n) - Zbus(jj,k); for kk = 1:nbus AT = Zbus(n,kk) - Zbus(k, kk); DELZ(jj,kk) = AP*AT/DM; end end Zbus = Zbus - DELZ; ntree(I) = 2; else,end else,end end end（2）Zbuildpi.m程序代碼：% This program forms the complex b

10、us impedance matrix by the method% of building algorithm. Bus zero is taken as reference.% This program is compatible with power flow data. function Zbus, linedata = zbuildpi(linedata, gendata, yload)ng=length(gendata(:,1);nlg=gendata(:,1);nrg=gendata(:,2); zg= gendata(:,2) + j*gendata(:,3);nl = lin

11、edata(:,1); nr = linedata(:,2); R = linedata(:,3);X = linedata(:,4);nbr=length(linedata(:,1); nbus = max(max(nl), max(nr);nc = length(linedata(1,:);for k=1:nbr if R(k) = inf | X(k) = inf R(k) = ; X(k) = ; else, endend if nc > 4BC = linedata(:,5); for n = 1:nbus yc(n) = 0; nlc(n) = 0; nrc(n) = n;

12、for k = 1:nbr if nl(k) = n | nr(k) = n yc(n) = yc(n) + j*BC(k); else, end end endelseif nc=4 yc= zeros(1, nbr); endnlc=nlc' nrc=nrc' yc=yc.'ZB = R + j*X;if exist('yload') = 1 yload = yload.' yc =yc + yload;else, endm = 0;for n = 1:nbus if abs(yc(n) =0 m=m+1; nlcc(m) = nlc(n);

13、 nrcc(m) = nrc(n); zc(m) = 1/yc(n); else, endendnlcc=nlcc' nrcc=nrcc' zc=zc.'nl=nlg; nlcc; nl; nr = nrg; nrcc; nr; ZB = zg; zc; ZB;linedata=nl nr real(ZB) imag(ZB);nbr= length(nl);Zbus = zeros(nbus, nbus);tree=0; %從參考總線0上添加一個分支 for I = 1:nbr ntree(I) = 1; if nl(I) = 0 | nr(I) = 0 if nl(I

14、) = 0 n = nr(I); elseif nr(I) = 0 n = nl(I); end if abs(Zbus(n, n) = 0 Zbus(n,n) = ZB(I); tree=tree+1; else Zbus(n,n) = Zbus(n,n)*ZB(I)/(Zbus(n,n) + ZB(I); end ntree(I) = 2; else,end end % 添加一個新總線分支到現(xiàn)有總線上while tree < nbus for n = 1:nbus nadd = 1; if abs(Zbus(n,n) = 0 for I = 1:nbr if nadd = 1 if

15、nl(I) = n | nr(I) = n if nl(I) = n k = nr(I); elseif nr(I) = n k = nl(I); end if abs(Zbus(k,k) = 0 for m = 1:nbus if m = n Zbus(m,n) = Zbus(m,k); Zbus(n,m) = Zbus(m,k); else, end end Zbus(n,n) = Zbus(k,k) + ZB(I); tree=tree+1; nadd = 2; ntree(I) = 2; else, end else, end else, end end else, endendend

16、 %增加兩個原有總線間的支路阻抗 for n = 1:nbus for I = 1:nbr if ntree(I) = 1 if nl(I) = n | nr(I) = n if nl(I) = n k = nr(I); elseif nr(I) = n k = nl(I); end DM = Zbus(n,n) + Zbus(k,k) + ZB(I) - 2*Zbus(n,k); for jj = 1:nbus AP = Zbus(jj,n) - Zbus(jj,k); for kk = 1:nbus AT = Zbus(n,kk) - Zbus(k, kk); DELZ(jj,kk) =

17、AP*AT/DM; end end Zbus = Zbus - DELZ; ntree(I) = 2; else,end else,end end end（3）symfault.m程序代碼：% 此程序用來計算電網(wǎng)的三相對稱故障% 計算前需要用戶生成節(jié)點阻抗矩陣% 節(jié)點阻抗矩陣可由函數(shù)Zbus = zbuild(zdata)生成%此程序需要用戶按提示輸入短路節(jié)點編號和過度電阻Zf% 向量V是可選參數(shù)，包含節(jié)點編號和復數(shù)電壓%V可由潮流程序lfgauss, lfnewton,decouple自動生成. % 如果V不存在，程序將默認故障前所有的節(jié)點電壓標幺值為1.0% 程序可得到總的故障電流, 節(jié)點

18、電壓幅值及輸電線路的電流function symfaul(zdata, Zbus, V) nl = zdata(:,1); nr = zdata(:,2); R = zdata(:,3);X = zdata(:,4);nc = length(zdata(1,:); if nc > 4 BC = zdata(:,5); elseif nc =4, BC = zeros(length(zdata(:,1), 1); endZB = R + j*X;nbr=length(zdata(:,1); nbus = max(max(nl), max(nr);if exist('V')

19、= 1 if length(V) = nbus V0 = V; else, endelse, V0 = ones(nbus, 1) + j*zeros(nbus, 1);endfprintf('Three-phase balanced fault analysis n')ff = 999;while ff > 0nf = input('Enter Faulted Bus No. -> '); rtn=isempty(nf); if rtn=1; nf=-1; end while nf <= 0 | nf > nbus fprintf(&#

20、39;Faulted bus No. must be between 1 & %g n', nbus) nf = input('Enter Faulted Bus No. -> '); rtn=isempty(nf); if rtn=1; nf=-1; end endrtz=1; while rtz=1 fprintf('nEnter Fault Impedance Zf = R + j*X in ') Zf = input('complex form (for bolted fault enter 0). Zf = ');

21、 rtz=isempty(Zf); endfprintf(' n')fprintf('Balanced three-phase fault at bus No. %gn', nf) If = V0(nf)/(Zf + Zbus(nf, nf);Ifm = abs(If); Ifmang=angle(If)*180/pi;fprintf('Total fault current = %8.4f per unit nn', Ifm)%fprintf(' p.u. nn', Ifm)fprintf('Bus Voltages d

22、uring fault in per unit nn')fprintf(' Bus Voltage Anglen')fprintf(' No. Magnitude degreesn') for n = 1:nbus if n=nf Vf(nf) = V0(nf)*Zf/(Zf + Zbus(nf,nf); Vfm = abs(Vf(nf); angv=angle(Vf(nf)*180/pi; else, Vf(n) = V0(n) - V0(n)*Zbus(n,nf)/(Zf + Zbus(nf,nf); Vfm = abs(Vf(n); angv=an

23、gle(Vf(n)*180/pi; end fprintf(' %4g', n), fprintf('%13.4f', Vfm),fprintf('%13.4fn', angv) end fprintf(' n') fprintf('Line currents for fault at bus No. %gnn', nf)fprintf(' From To Current Anglen')fprintf(' Bus Bus Magnitude degreesn') for n= 1:

24、nbus for I = 1:nbr if nl(I) = n | nr(I) = n if nl(I) =n k = nr(I); elseif nr(I) = n k = nl(I); end if k=0 Ink = (V0(n) - Vf(n)/ZB(I); Inkm = abs(Ink); th=angle(Ink); if real(Ink) > 0 fprintf(' G '), fprintf('%7g',n), fprintf('%12.4f', Inkm) fprintf('%12.4fn', th*18

25、0/pi) elseif real(Ink) =0 & imag(Ink) < 0 fprintf(' G '), fprintf('%7g',n), fprintf('%12.4f', Inkm) fprintf('%12.4fn', th*180/pi) else, end Ign=Ink; elseif k = 0 Ink = (Vf(n) - Vf(k)/ZB(I)+BC(I)*Vf(n); Inkm = abs(Ink); th=angle(Ink); if real(Ink) > 0 fprintf

26、('%7g', n), fprintf('%10g', k), fprintf('%12.4f', Inkm), fprintf('%12.4fn', th*180/pi) elseif real(Ink) =0 & imag(Ink) < 0 fprintf('%7g', n), fprintf('%10g', k), fprintf('%12.4f', Inkm), fprintf('%12.4fn', th*180/pi) else, end el

27、se, end else, end end if n=nf fprintf('%7g',n), fprintf(' F'), fprintf('%12.4f', Ifm) fprintf('%12.4fn', Ifmang) else, endend resp=0; while strcmp(resp, 'n')=1 & strcmp(resp, 'N')=1 & strcmp(resp, 'y')=1 & strcmp(resp, 'Y')=

28、1 resp = input('Another fault location? Enter ''y'' or ''n'' within single quote -> '); if strcmp(resp, 'n')=1 & strcmp(resp, 'N')=1 & strcmp(resp, 'y')=1 & strcmp(resp, 'Y')=1 fprintf('n Incorrect reply, try

29、again nn'), end end if resp = 'y' | resp = 'Y' nf = 999; else ff = 0; endend 二、不對稱短路故障不對稱故障包括單相接地短路(single line-to-ground fault)，兩相短路(line-to-line fault )，和兩相接地短路(double line-to-ground fault)。這里介紹三個程序來分析不對稱短路故障：單相接地lgfault(zdata0,zbus0,zdata1,zbus1,zdata2,zbus2,V) 兩相短路llfault(zda

30、ta1,zbus1,zdata2,zbus2,V) 兩相接地短路dlgfault(zdata0,zbus0,zdata1,zbus1,zdata2,zbus2,V)對于lgfault和dlgfault需要輸入正序、零序、負序節(jié)點阻抗矩陣Zbus1、Zbus0和Zbus2，而llfault只需輸入正序和負序節(jié)點阻抗矩陣Zbus1、Zbus2，最后參數(shù)V是可以選擇的。如果用戶未輸入V，程序則默認所有故障前的電壓為1.0pu。如果用戶輸入?yún)?shù)V，則需要輸入節(jié)點編號以及節(jié)點電壓的復數(shù)值。變量V可由潮流計算程序獲得。零序網(wǎng)絡和正序網(wǎng)絡的節(jié)點阻抗矩陣分別通過函數(shù)Zbus0=zbuild（zdata0）和Z

31、bus1=zbuild（zdata1）獲得。參數(shù)zdata1包含正序網(wǎng)絡阻抗，參數(shù)zdata0包含零序網(wǎng)絡阻抗。參數(shù)zdata1、zdata2、zdata0均為e×4階矩陣，e為元件的個數(shù)。第一列和第二列是節(jié)點編號，第三列和第四列分別是線路電阻和電抗（均為標幺值）。節(jié)點0為參考節(jié)點，發(fā)電機節(jié)點到節(jié)點0的阻抗為發(fā)電機阻抗，可以是發(fā)電機的次暫態(tài)電抗、暫態(tài)電抗或者同步電抗，而且矩陣還包括并聯(lián)電容和負荷阻抗等。負序網(wǎng)絡和正序網(wǎng)絡有同樣的拓撲結構。線路和變壓器的負序阻抗等于正序阻抗，而發(fā)電機的負序阻抗不等于正序阻抗，但在大型電力系統(tǒng)故障分析中通常按相等來計算。零序網(wǎng)絡的拓撲結構不同于正序網(wǎng)絡，

32、零序網(wǎng)絡的構建取決于變壓器的繞組連接方式，除了兩側中性點都接地的Y-Y連接的變壓器，對于其他連接方式的變壓器，一次側和二次側的零序網(wǎng)絡都是隔離開的。在這樣的連接方式的零序網(wǎng)絡中,，用inf填寫對應的電阻和電抗列。對于Y-（中性點接地）連接的變壓器，Y側的阻抗到節(jié)點0的數(shù)據(jù)需要填寫。當中性點接地阻抗為Xn時，需要填寫為3Xn。 在不對稱故障計算時，還需要一個函數(shù)用來求節(jié)點阻抗矩陣，Zbus=zbuildpi（linedata，gendata，yload），函數(shù)與潮流計算程序是相容的。第一個參數(shù)linedata與潮流程序所需要的數(shù)據(jù)相同。第一列和第二列是節(jié)點編號；第三列到第五列分別是線路的電阻、電

33、抗、1/2線路電納（均為標幺值，以指定的MVA為基準功率）；最后一列為變壓器分接頭位置，對于線路，在這一列中必須填寫1。在潮流計算程序中，發(fā)電機電抗不包含在參數(shù)linedata中，而是包含在參數(shù)gendata中。 參數(shù)gendata是一個(g×4)階矩陣，其四列分別為參考節(jié)點0、發(fā)電機節(jié)點編號、發(fā)電機電阻和電抗。最后一個參數(shù)yload是可選的，這是一個兩列的矩陣，包括節(jié)點編號和負荷復數(shù)導納值。這些數(shù)據(jù)可以從潮流計算程序lfgauss，lfnewton，decouple執(zhí)行過程中自動生成。程序運行時會提示用戶輸入故障節(jié)點編號和故障阻抗Zf，運行結果可得到故障電流，節(jié)點電壓和線路電流。不

34、對稱短路故障計算使用的程序如下：（1）lgfault.m程序代碼：% The program lgfault is designed for the single line-to-ground% fault analysis of a power system network. The program requires% the positive-, negative- and zero-sequence bus impedance matrices,% Zbus1 Zbus2,and Zbus0.The bus impedances matrices may be defined% by t

35、he user, obtained by the inversion of Ybus or it may be% determined either from the function Zbus = zbuild(zdata)% or the function Zbus = zbuildpi(linedata, gendata, yload).% The program prompts the user to enter the faulted bus number% and the fault impedance Zf. The prefault bus voltages are% defi

36、ned by the reserved Vector V. The array V may be defined or% it is returned from the power flow programs lfgauss, lfnewton,% decouple or perturb. If V does not exist the prefault bus voltages% are automatically set to 1.0 per unit. The program obtains the% total fault current, bus voltages and line

37、currents during the fault. function lgfault(zdata0, Zbus0, zdata1, Zbus1, zdata2, Zbus2, V)if exist('zdata2') = 1zdata2=zdata1;else, endif exist('Zbus2') = 1Zbus2=Zbus1;else, endnl = zdata1(:,1); nr = zdata1(:,2);nl0 = zdata0(:,1); nr0 = zdata0(:,2);nbr=length(zdata1(:,1); nbus = max

38、(max(nl), max(nr);nbr0=length(zdata0(:,1);R0 = zdata0(:,3); X0 = zdata0(:,4);R1 = zdata1(:,3); X1 = zdata1(:,4);R2 = zdata1(:,3); X2 = zdata1(:,4); for k=1:nbr0 if R0(k)=inf | X0(k) =inf R0(k) = ; X0(k) = ; else, endendZB1 = R1 + j*X1; ZB0 = R0 + j*X0;ZB2 = R2 + j*X2; if exist('V') = 1 if le

39、ngth(V) = nbus V0 = V; else, endelse, V0 = ones(nbus, 1) + j*zeros(nbus, 1);endfprintf('nLine-to-ground fault analysis n')ff = 999;while ff > 0nf = input('Enter Faulted Bus No. -> ');rtn=isempty(nf); if rtn=1; nf=-1; end while nf <= 0 | nf > nbus fprintf('Faulted bus

40、No. must be between 1 & %g n', nbus) nf = input('Enter Faulted Bus No. -> '); rtn=isempty(nf); if rtn=1; nf=-1; end endrtz=1; while rtz=1 fprintf('nEnter Fault Impedance Zf = R + j*X in ') Zf = input('complex form (for bolted fault enter 0). Zf = '); rtz=isempty(Zf

41、); endfprintf(' n')fprintf('Single line to-ground fault at bus No. %gn', nf)a =cos(2*pi/3)+j*sin(2*pi/3);sctm = 1 1 1; 1 a2 a; 1 a a2;Ia0 = V0(nf)/(Zbus1(nf,nf)+Zbus2(nf, nf)+ Zbus0(nf, nf)+3*Zf); Ia1=Ia0; Ia2=Ia0;I012=Ia0; Ia1; Ia2;Ifabc = sctm*I012;Ifabcm = abs(Ifabc);fprintf('

42、Total fault current = %9.4f per unitnn', Ifabcm(1)fprintf('Bus Voltages during the fault in per unit nn')fprintf(' Bus -Voltage Magnitude- n')fprintf(' No. Phase a Phase b Phase c n') for n = 1:nbusVf0(n)= 0 - Zbus0(n, nf)*Ia0;Vf1(n)= V0(n) - Zbus1(n, nf)*Ia1;Vf2(n)= 0 -

43、Zbus2(n, nf)*Ia2;Vabc = sctm*Vf0(n); Vf1(n); Vf2(n);Va(n)=Vabc(1); Vb(n)=Vabc(2); Vc(n)=Vabc(3);fprintf(' %5g',n)fprintf(' %11.4f', abs(Va(n),fprintf(' %11.4f', abs(Vb(n)fprintf(' %11.4fn', abs(Vc(n)endfprintf(' n')fprintf('Line currents for fault at bus N

44、o. %gnn', nf)fprintf(' From To -Line Current Magnitude- n')fprintf(' Bus Bus Phase a Phase b Phase c n')for n= 1:nbus for I = 1:nbr if nl(I) = n | nr(I) = n if nl(I) =n k = nr(I); elseif nr(I) = n k = nl(I); end if k = 0 Ink1(n, k) = (Vf1(n) - Vf1(k)/ZB1(I); Ink2(n, k) = (Vf2(n)

45、- Vf2(k)/ZB2(I); else, end else, end end for I = 1:nbr0 if nl0(I) = n | nr0(I) = n if nl0(I) =n k = nr0(I); elseif nr0(I) = n k = nl0(I); end if k = 0 Ink0(n, k) = (Vf0(n) - Vf0(k)/ZB0(I); else, end else, end end for I = 1:nbr if nl(I) = n | nr(I) = n if nl(I) =n k = nr(I); elseif nr(I) = n k = nl(I

46、); end if k = 0 Inkabc = sctm*Ink0(n, k); Ink1(n, k); Ink2(n, k); Inkabcm = abs(Inkabc); th=angle(Inkabc); if real(Inkabc(1) > 0 fprintf('%7g', n), fprintf('%10g', k), fprintf(' %11.4f', abs(Inkabc(1),fprintf(' %11.4f', abs(Inkabc(2) fprintf(' %11.4fn', abs

47、(Inkabc(3) elseif real(Inkabc(1) =0 & imag(Inkabc(1) < 0 fprintf('%7g', n), fprintf('%10g', k), fprintf(' %11.4f', abs(Inkabc(1),fprintf(' %11.4f', abs(Inkabc(2) fprintf(' %11.4fn', abs(Inkabc(3) else, end else, end else, end end if n=nf fprintf('%7

48、g',n), fprintf(' F'), fprintf(' %11.4f', Ifabcm(1),fprintf(' %11.4f', Ifabcm(2) fprintf(' %11.4fn', Ifabcm(3) else, endendresp=0; while strcmp(resp, 'n')=1 & strcmp(resp, 'N')=1 & strcmp(resp, 'y')=1 & strcmp(resp, 'Y')=

49、1 resp = input('Another fault location? Enter ''y'' or ''n'' within single quote -> '); if strcmp(resp, 'n')=1 & strcmp(resp, 'N')=1 & strcmp(resp, 'y')=1 & strcmp(resp, 'Y')=1 fprintf('n Incorrect reply, try

50、again nn'), end end if resp = 'y' | resp = 'Y' nf = 999; else ff = 0; endend （2）llfault.m程序代碼：% The program llfault is designed for the line-to-line% fault analysis of a power system network. The program requires% the positive- and negative-sequence bus impedance matrices,% Zbus1, and Zbus2.The bus impedances matrices may be defined% by the user, obtained by the inversion of Ybus or it may be% determined either from the function Zbus = zbuild(zdata)% or the function Zbus = zbuildpi(linedata, gendata,