基于牛顿拉夫逊法潮流计算的matlab实验报告定稿

基于牛顿拉夫逊法潮流计算的matlab实验报告一、实验目的应用MATLAB语言编写具有一定通用性的牛顿-拉夫逊法潮流计算程序。要求：（1）潮流计算方法为牛顿-拉夫逊法。（2）编程语言为MATLAB。（3）程序具有较强通用性。二、程序流程图所用公式其中开始形成节点导纳矩阵开始形成节点导纳矩阵输入原始数据设节点电压，i=1,2…,n,is置迭代次数置节点号i=1按式（3-3），（3-4）计算雅克比矩阵元素按式（3-2）计算节点的，，节点的，求解修正方程式，得，雅克比矩阵是否已全部形成？计算平衡节点及PV节点功率求，迭代次数k=k+1i=i+1？潮流计算完成计算各节点电压的新值：求解问题及其结果IEEE-美国新英格兰10机39节点测试系统系统单线图二、系统参数1）系统容量基准值为100MVA。2)负荷数据见表D-1表D-1负荷数据节点号有功（MW）无功（Mvar）节点号有功（MW）无功（Mvar）BUS-3BUS-4BUS-7BUS-8BUS-12BUS-15BUS-16BUS-18BUS-20322.0500.0233.8522.08.5320.0329.0158.0680.02.4184.084.0176.088.0153.032.330.0103.0BUS-21BUS-23BUS-24BUS-25BUS-26BUS-27BUS-28BUS-29BUS-39274.0247.5308.0224.0139.0281.0206.0283.51104.0115.084.6-75.527.626.93）发电机数据见表D-2表D-2发电机数据发电机节点号有功（MW）无功（Mvar）电压（p.u.）PVPQPQPVPVPQPVPVPVBUS-30BUS-32BUS-33BUS-34BUS-35BUS-36BUS-37BUS-38BUS-39250.00650.00632.00508.00650.00560.00540.00830.001000.00175.90103.3596.881.047501.012301.049301.027801.026501.03000平衡节点BUS-310.0(电压幅角)1.0（幅值）4）线路参数见表D-3表D-3线路参数序号节点I节点JR(p.u.)X(p.u.)B/2(p.u.)LN1LN2LN3LN4LN5LN6LN7LN8LN9LN10LN11LN12LN13LN14LN15LN16LN17LN18LN19LN20LN21LN22LN23LN24LN25LN26LN27LN28LN29LN30LN31LN32LN33LN34BUS-2BUS-39BUS-3BUS-25BUS-4BUS-18BUS-5BUS-14BUS-6BUS-8BUS-7BUS-11BUS-8BUS-9BUS-39BUS-11BUS-13BUS-14BUS-15BUS-16BUS-17BUS-19BUS-21BUS-24BUS-18BUS-27BUS-22BUS-23BUS-24BUS-26BUS-27BUS-28BUS-29BUS-29BUS-1BUS-1BUS-2BUS-2BUS-3BUS-3BUS-4BUS-4BUS-5BUS-5BUS-6BUS-6BUS-7BUS-8BUS-9BUS-10BUS-10BUS-13BUS-14BUS-15BUS-16BUS-16BUS-16BUS-16BUS-17BUS-17BUS-21BUS-22BUS-23BUS-25BUS-26BUS-26BUS-26BUS-280.003500.001000.001300.007000.001300.001100.000800.000800.000200.000800.000600.000700.000400.002300.001000.000400.000400.000900.001800.000900.000700.001600.000800.000300.000700.001300.000800.000600.002200.003200.001400.004300.005700.001400.041100.025000.015100.008600.021300.013300.012800.012900.002600.011200.009200.008200.004600.036300.025000.004300.004300.010100.021700.009400.008900.019500.013500.005900.008200.017300.014000.009600.035000.032300.014700.047400.062500.015100.349350.37500.128600.073000.110700.106900.067100.069100.021700.073800.056500.069450.039000.190200.600000.036450.036450.086150.183000.085500.067100.152000.127400.034000.065950.160800.128250.092300.180500.256500.119800.390100.514500.12450LN35:BUS-4接有并联电容器，B4＝1.0000LN36:BUS-5接有并联电容器，B4＝2.00005)变压器参数见表D-4表D-4变压器参数序号节点I节点JR(p.u.)X(p.u.)变比(%)TR37TR38TR39TR40TR41TR42TR43TR44TR45TR46TR47TR48BUS-11BUS-13BUS-30BUS-31BUS-32BUS-34BUS-33BUS-35BUS-36BUS-37BUS-38BUS-20BUS-12BUS-12BUS-2BUS-6BUS-10BUS-20BUS-19BUS-22BUS-23BUS-25BUS-29BUS-190.001600.001600.000000.000000.000000.000900.000700.000000.000500.000600.000800.000700.043500.043500.018100.025000.020000.018000.014200.014300.027200.023200.015600.01380100.6100.6102.5107.0107.0100.9107.0102.5100.0102.5102.5106.0%IEEE-美国新英格兰10机39节点测试系统%123456%busvoltanglepqtypebus=[11.00000.000.000.00121.00000.000.000.00131.00000.00-3.22-0.024141.00000.00-5.00-1.84151.00000.000.000.00161.00000.000.000.00171.00000.00-2.338-0.84181.00000.00-5.22-1.76191.00000.000.000.001101.00000.000.000.001111.00000.000.000.001121.00000.00-0.085-0.881131.00000.000.000.001141.00000.000.000.001151.00000.00-3.20-1.531161.00000.00-3.29-0.3231171.00000.000.000.001181.00000.00-1.58-0.301191.00000.000.000.001201.00000.00-6.80-1.031211.00000.00-2.74-1.151221.00000.000.000.001231.00000.00-2.475-1.151241.00000.00-3.08-0.9221251.00000.00-2.24-0.4721261.00000.00-1.39-0.171271.00000.00-2.81-0.7551281.00000.00-2.06-0.2761291.00000.00-2.835-0.2691301.04750.002.500.002311.00000.000.000.003321.00000.006.501.7591331.00000.006.321.03351341.01230.005.080.002351.04930.006.500.002361.00000.005.600.96881371.02780.005.400.002381.02650.008.300.002391.03000.00-1.040.002];%1234567%line:frombustobusR,X,G,B/2Kline=[2 1 0.00350 0.04110 00.349350;39 1 0.00100 0.02500 00.375000;3 2 0.00130 0.01510 00.128600;25 2 0.00700 0.00860 00.073000;4 3 0.00130 0.02130 00.110700;18 3 0.00110 0.01330 00.106900;5 4 0.00080 0.01280 00.067100;14 4 0.00080 0.01290 00.069100;6 5 0.00020 0.00260 00.021700;8 5 0.00080 0.01120 00.073800;7 6 0.00060 0.00920 00.056500;11 6 0.00070 0.00820 00.069450;8 7 0.00040 0.0046000.039000;9 8 0.00230 0.03630 00.190200;39 9 0.00100 0.02500 00.600000;11 10 0.00040 0.00430 00.036450;13 10 0.00040 0.00430 00.036450;14 13 0.00090 0.01010 00.086150;15 14 0.00180 0.02170 00.183000;16 15 0.00090 0.00940 00.085500;17 16 0.00070 0.00890 00.067100;19 16 0.00160 0.01950 00.152000;21 16 0.00080 0.01350 00.127400;24 16 0.00030 0.00590 00.034000;18 17 0.00070 0.00820 00.065950;27 17 0.00130 0.01730 00.160800;22 21 0.00080 0.01400 00.128250;23 22 0.00060 0.00960 00.092300;24 23 0.00220 0.03500 00.180500;26 25 0.00320 0.03230 00.256500;27 26 0.00140 0.01470 00.119800;28 26 0.00430 0.04740 00.390100;29 26 0.00570 0.06250 00.514500;29 28 0.00140 0.01510 00.124500;400001.00000;500002.00000;11 12 0.00160 0.04350 00100.60000/100;13 12 0.00160 0.04350 00100.60000/100;30 2 0.00000 0.01810 00102.50000/100;31 6 0.00000 0.02500 00107.00000/100;32 10 0.00000 0.02000 00107.00000/100;34 20 0.00090 0.01800 00100.90000/100;33 19 0.00070 0.01420 00107.00000/100;35 22 0.00000 0.01430 00102.50000/100;36 23 0.00050 0.02720 00100.00000/100;37 25 0.00060 0.02320 00102.50000/100;38 29 0.00080 0.01560 00102.50000/100;20 19 0.00070 0.01380 00106.00000/100];计算结果牛顿－拉夫逊法潮流计算结果节点计算结果：n节点节点电压节点相角（角度）节点注入功率11.049185-8.8749910.000000+j0.00000021.053167-6.3671800.000000+j0.00000031.041493-9.207297-3.220000+j-0.02400041.036574-10.042585-5.000000+j-1.84000051.044652-8.9592370.000000+j0.00000061.043883-8.2931040.000000+j0.00000071.032645-10.342431-2.338000+j-0.84000081.031177-10.811816-5.220000+j-1.76000091.042715-10.5956480.000000+j0.000000101.046426-6.0104760.000000+j0.000000111.044322-6.7924620.000000+j0.000000121.030736-6.795388-0.085000+j-0.880000131.042351-6.6754910.000000+j0.000000141.036310-8.2323370.000000+j0.000000151.018517-8.519794-3.200000+j-1.530000161.025492-7.051856-3.290000+j-0.323000171.032750-8.0771180.000000+j0.000000181.034779-8.936485-1.580000+j-0.300000191.044862-2.3821690.000000+j0.000000200.988148-3.811032-6.800000+j-1.030000211.024926-4.596980-2.740000+j-1.150000221.042650-0.0705120.000000+j0.000000231.032952-0.245457-2.475000+j-1.150000241.021125-6.906503-3.080000+j-0.922000251.060163-4.952002-2.240000+j-0.472000261.052697-6.205207-1.390000+j-0.170000271.037683-8.217337-2.810000+j-0.755000281.050444-2.695196-2.060000+j-0.276000291.0501630.063077-2.835000+j-0.269000301.0043921.5947816.500000+j1.759000310.9916322.8925726.320000+j1.033500321.0505397.7977865.600000+j0.968800331.047500-3.9575982.500000+j1.211174341.0123001.3857745.080000+j1.826359351.0493004.9253246.500000+j2.637566361.0278001.8194765.400000+j-0.108224371.0265007.1255798.300000+j0.214225381.030000-10.390696-1.040000+j-2.291639391.0000000.0000005.628660+j1.384403线路计算结果：n节点I节点J线路功率S(I,J)线路功率S(J,I)线路损耗dS(I,J)211.178698+j-0.360055-1.174311+j-0.3604810.004386+j-0.7205363916.405845+j-2.096152-6.361848+j2.4082870.043997+j0.31213532-3.633961+j-0.5426133.649983+j0.4465770.016021+j-0.0960362522.370242+j-1.109311-2.328681+j0.9973560.041562+j-0.11195543-0.750370+j-0.3071720.751094+j0.0800140.000724+j-0.2271591830.337560+j-0.663855-0.337133+j0.4385990.000427+j-0.225256541.635254+j0.499000-1.633054+j-0.6091190.002200+j-0.1101191442.621711+j-0.216428-2.616576+j0.1507770.005135+j-0.065651654.826035+j-0.675350-4.821682+j0.6846070.004353+j0.00925785-3.178130+j-1.0418363.186428+j0.9989890.008297+j-0.04284776-4.249274+j-0.9695594.259899+j1.0106570.010625+j0.0410981163.465003+j-0.270003-3.457273+j0.2091360.007730+j-0.06086687-1.909893+j-0.1967321.911274+j0.1295590.001381+j-0.067173980.132235+j0.116464-0.131977+j-0.5214320.000258+j-0.4049683997.617154+j-1.902126-7.557438+j2.1426870.059717+j0.2405611110-3.483660+j-0.2030643.488121+j0.1713520.004461+j-0.0317121310-3.008372+j-0.7304893.011879+j0.6886800.003508+j-0.0418091413-2.934129+j-0.4114632.941429+j0.3072640.007300+j-0.1041991514-0.311115+j-0.9985560.312417+j0.6278910.001303+j-0.37066516152.896296+j0.430232-2.888885+j-0.5314440.007411+j-0.1012121716-2.048841+j0.9507402.052282+j-1.0491220.003441+j-0.09838319164.542969+j0.681545-4.511670+j-0.6258730.031300+j0.05567221163.324778+j-0.302389-3.316338+j0.1770060.008440+j-0.12538324160.410793+j-0.811601-0.410571+j0.7447570.000222+j-0.0668441817-1.917560+j0.3638551.920087+j-0.4752080.002527+j-0.1113532717-0.128621+j0.1326480.128754+j-0.4755310.000133+j-0.34288422216.093176+j1.070437-6.064778+j-0.8476110.028398+j0.2228262322-0.406149+j-1.1160630.406824+j0.9280400.000675+j-0.1880242423-3.490793+j-0.1103993.516516+j0.1388370.025723+j0.0284382625-0.771398+j-0.4428810.773189+j-0.1115800.001791+j-0.5544602726-2.681379+j-0.8876482.691475+j0.7319000.010096+j-0416063+j-0.565082-1.408178+j-0.2107470.007885+j-0.77583029261.921038+j-0.679443-1.901899+j-0.2482720.019138+j-0.92771529283.491624+j-0.395924-3.476063+j0.2890820.015561+j-0.106842400.000000+j-1.0744850.000000+j0.0000000.000000+j-1.074485500.000000+j-2.1825960.000000+j0.0000000.000000+j-2018656+j0.473066-0.018327+j-0.4641260.000329+j0.00894013120.066943+j0.423225-0.066673+j-0.4158740.000270+j0.0073513027.897633+j-0.731582-7.897633+j1.8602770.000000+j1.1286953167.506817+j1.371343-7.506817+j0.1091530.000000+j1.480496321012.260592+j5.296517-12.260592+j-2.0640070.000000+j3.23250934205.080000+j1.826359-5.054406+j-1.3144730.025594+j0.5118863319-1.716763+j5.3489101.736896+j-4.9405040.020133+j0.40840535226.500000+j2.637566-6.500000+j-1.9984770.000000+j0.63908936231.402814+j-0.195113-1.401865+j0.2467630.000949+j0.05165037259.586236+j0.419689-9.533808+j1.6075170.052428+j2.0272063829-12.165903+j2.10659312.280860+j0.1350620.114957+j2.2416552019-1.745594+j0.2844731.747837+j-0.2402650.002242+j0.044208结果分析：此程序的运行结果和试验程序给出的结果是一致的。说明程序无误，但在精确度上有微小差异，这主要是和导纳矩阵的精确度以及显示精度有关。心得：本程序分模块进行，先是排序，再是求导纳阵，然后求雅阁比，再进行迭代运算，程序本身很简洁明了，运行的时候只需要在matlab里输入main就行了，然后打开BUS和line所在的.m文件，结果就会自动存在result文件中了，通过编写牛顿拉夫逊法matlab潮流计算程序复习了潮流计算的知识，也实现了计算机算法附录：实验源程序：Main函数：clear[dfile,pathname]=uigetfile('*.m','SelectDataFile');ifpathname==0error('youmustselectavaliddatafile')elselfile=length(dfile);%stripoff.meval(dfile(1:lfile-2));end[nb,mb]=size(bus);[nl,ml]=size(line);%计算bus和line矩阵的行数和列数[bus,line,nPQ,nPV,nodenum]=Num(bus,line);%对节点重新排序的子程序Y=y(bus,line)%计算节点导纳矩阵的子程序myf=fopen('Result.m','w');fprintf(myf,'计算结果');fclose(myf);%在当前目录下生成“Result.m”文件，写入节点导纳矩阵formatlongEPS=1.0e-10;%设定误差精度fort=1:100%开始迭代计算，设定最大迭代次数为100，以便不收敛情况下及时跳出[dP,dQ]=dPQ(Y,bus,nPQ,nPV);%计算功率偏差dP和dQ的子程序J=Jac(bus,Y,nPQ);%计算雅克比矩阵的子程序UD=zeros(nPQ,nPQ);fori=1:nPQUD(i,i)=bus(i,2);%生成电压对角矩阵endenddAngU=J\[dP;dQ];dAng=dAngU(1:nb-1,1);%计算相角修正量dU=UD*(dAngU(nb:nb+nPQ-1,1));%计算电压修正量bus(1:nPQ,2)=bus(1:nPQ,2)-dU;%修正电压bus(1:nb-1,3)=bus(1:nb-1,3)-dAng;%修正相角if(max(abs(dU))<EPS)&(max(abs(dAng))<EPS)breakend%判断是否满足精度误差，如满足则跳出，否则返回继续迭代endbus=PQ(bus,Y,nPQ,nPV);%计算每个节点的有功和无功注入的子程序[bus,line]=ReNum(bus,line,nodenum);%对节点恢复编号的子程序YtYm=YtYm(line);%计算线路的等效Yt和Ym的子程序，以计算线路潮流bus_res=bus_res(bus);%计算节点数据结果的子程序S_res=S_res(bus,line,YtYm);%计算线路潮流的子程序myf=fopen('Result.m','a');fprintf(myf,'牛顿－拉夫逊法潮流计算结果节点计算结果：n节点节点电压节点相角（角度）节点注入功率\n');fori=1:nbfprintf(myf,'%2.0f',bus_res(i,1));fprintf(myf,'%10.6f',bus_res(i,2));fprintf(myf,'%10.6f',bus_res(i,3));fprintf(myf,'%10.6f+j%10.6f\n',real(bus_res(i,4)),imag(bus_res(i,4)));endfprintf(myf,'n线路计算结果：n节点I节点J线路功率S(I,J)线路功率S(J,I)线路损耗dS(I,J)\n');fori=1:nlfprintf(myf,'%2.0f',S_res(i,1));fprintf(myf,'%2.0f',S_res(i,2));fprintf(myf,'%10.6f+j%10.6f',real(S_res(i,3)),imag(S_res(i,3)));fprintf(myf,'%10.6f+j%10.6f',real(S_res(i,4)),imag(S_res(i,4)));fprintf(myf,'%10.6f+j%10.6f\n',real(S_res(i,5)),imag(S_res(i,5)));endfclose(myf);%迭代结束后继续在“Result.m”写入节点计算结果和线路计算结果程序结束"Num.m"作用为对节点重排序，并修改相应的线路数据function[bus,line,nPQ,nPV,nodenum]=Num(bus,line)[nb,mb]=size(bus);[nl,ml]=size(line);nSW=0;%numberofswingbuscounternPV=0;%numberofPVbuscounternPQ=0;%numberofPQbuscounterfori=1:nb,%nb为总节点数type=bus(i,6);iftype==3,nSW=nSW+1;%incrementswingbuscounterSW(nSW,:)=bus(i,:);elseiftype==2,nPV=nPV+1;%incrementPVbuscounterPV(nPV,:)=bus(i,:);elsenPQ=nPQ+1;%incrementPQbuscounterPQ(nPQ,:)=bus(i,:);endend;bus=[PQ;PV;SW];newbus=[1:nb]';nodenum=[newbusbus(:,1)];bus(:,1)=newbus;fori=1:nlforj=1:2fork=1:nbifline(i,j)==nodenum(k,2)line(i,j)=nodenum(k,1);breakendendendend"y.m"作用为计算节点导纳矩阵functionY=y(bus,line)[nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);fork=1:nlI=line(k,1);%读入线路参数J=line(k,2);Zt=line(k,3)+j*line(k,4);Yt=1/Zt;Ym=line(k,5)+j*line(k,6);K=line(k,7);if(K==0)&(J~=0)%普通线路:K=0;Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+Yt+Ym;Y(I,J)=Y(I,J)-Yt;Y(J,I)=Y(I,J);endif(K==0)&(J==0)%对地支路:K=0,J=0,R=X=0;Y(I,I)=Y(I,I)+Ym;endifK>0%变压器线路:Zt和Ym为折算到i侧的值,K在j侧Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+Yt/K/K;Y(I,J)=Y(I,J)-Yt/K;Y(J,I)=Y(I,J);endifK<0%变压器线路:Zt和Ym为折算到K侧的值,K在i侧Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+K*K*Yt;Y(I,J)=Y(I,J)+K*Yt;Y(J,I)=Y(I,J);endend"dPQ.m"作用为计算功率偏差function[dP,dQ]=dPQ(Y,bus,nPQ,nPV)%nPQ、nPV为相应节点个数n=nPQ+nPV+1;%总节点个数dP=bus(1:n-1,4);dQ=bus(1:nPQ,5);%对dP和dQ赋初值PV节点不需计算dQ平衡节点不参与计算fori=1:n-1forj=1:ndP(i,1)=dP(i,1)-bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));ifi<nPQ+1dQ(i,1)=dQ(i,1)-bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));endendend%利用循环计算求取dP和dQ"Jac.m"作用为计算雅克比矩阵functionJ=Jac(bus,Y,nPQ)[nb,mb]=size(bus);H=zeros(nb-1,nb-1);N=zeros(nb-1,nPQ);K=zeros(nPQ,nb-1);L=zeros(nPQ,nPQ);%将雅克比矩阵分块，即：J=[HN;KL]，并初始化Qi=zeros(nb-1,1);Pi=zeros(nb-1,1);fori=1:nb-1forj=1:nbPi(i,1)=Pi(i,1)+bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));Qi(i,1)=Qi(i,1)+bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));endend%初始化并计算Qi和Pifori=1:nb-1forj=1:nb-1ifi~=jH(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));elseH(i,j)=Qi(i,1)+imag(Y(i,j))*((bus(i,2))^2);end%分别计算H矩阵的对角及非对角元素ifj<nPQ+1ifi~=jN(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));elseN(i,j)=-Pi(i,1)-real(Y(i,j))*((bus(i,2))^2);endend%分别计算N矩阵的对角及非对角元素ifi<nPQ+1ifi~=jK(i,j)=bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));elseK(i,j)=-Pi(i,1)+real(Y(i,j))*((bus(i,2))^2);end%分别计算K矩阵的对角及非对角元素ifj<nPQ+1ifi~=jL(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));elseL(i,j)=-Qi(i,1)+imag(Y(i,j))*((bus(i,2))^2);endend%分别计算L矩阵的对角及非对角元素endendendJ=[HN;KL];%生成雅克比矩阵"PQ.m"作用为计算每个节点的功率注入functionbus=PQ(bus,Y,nPQ,nPV)n=nPQ+nPV+1;%n为总节点数fori=nPQ+1:n-1forj=1:nbus(i,5)=bus(i,5)+bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));endend%利用公式计算PV节点的无功注入forj=1:nbus(n,4)=bus(n,4)+bus(n,2)*bus(j,2)*(real(Y(n,j))*cos(bus(n,3)-bus(j,3))+imag(Y(n,j))*sin(bus(n,3)-bus(j,3)));bus(n,5)=bus(n,5)+bus(n,2)*bus(j,2)*(real(Y(n,j))*sin(bus(n,3)-bus(j,3))-imag(Y(n,j))*cos(bus(n,3)-bus(j,3)));end%计算平衡节点的无功及有功注入"ReNum.m"作用为对节点和线路数据恢复编号function[bus,line]=ReNum(bus,line,nodenum)[nb,mb]=size(bus);[nl,ml]=size(line);bus_temp=zeros(nb,mb);%bus_temp矩阵用于临时存放bus矩阵的数据k=1;fori=1:nbforj=1:nbifbus(j,1)==kbus_temp(k,:)=bus(j,:);k=k+1;endendend%利用bus矩阵的首列编号重新对bus矩阵排序并存入bus_temp矩阵中bus=bus_temp;%重新赋值回bus，完成bus矩阵的重新编号fori=1:nlforj=1:2fork=1:nbifline(i,j)==nodenum(k,1)line(i,j)=nodenum(k,2);breakendendendend%恢复line的编号"YtYm.m"作用为计算线路的等效Yt和Ym，以计算线路潮流functionYtYm=YtYm(line)[nl,ml]=size(line);YtYm=zeros(nl,5);%对YtYm矩阵赋初值0YtYm(:,1:2)=line(:,1:2);%矩阵前两列为线路两段节点编号，后三列分别为线路等效Yt，i侧的等效Ym，j侧的等效Ymj=sqrt(-1);fork=1:nlI=line(k,1);J=line(k,2);Zt=line(k,3)+j*line(k,4);ifZt~=0Yt=1/Zt;endYm=line(k,5)+j*line(k,6);K=line(k,7);if(K==0)&(J~=0)%普通线路:K=0YtYm(k,3)=Yt;YtYm(k,4)=Ym;YtYm(k,5)=Ym;endif(K==0)&(J==0)%对地支路:K=0,J=0,R=X=0YtYm(k,4)=Ym;endifK>0%变压器线路:Zt和Ym为折算到i侧的值,K在j侧YtYm(k,3)=Yt/K;YtYm(k,4)=Ym+Yt*(K-1)/K;YtYm(k,5)=Yt*(1-K)/K/K;endifK<0%变压器线路:Zt和Ym为折算到K侧的值,K在i侧YtYm(k,3)=-Yt*K;YtYm(k,4)=Ym+Yt*(1+K);YtYm(k,5)=Yt*(K^2+K);endend"bus_res.m"计算并返回节点数据结果functionbus_res=bus_res(bus)[nb,mb]=size(bus);bus_res=zeros(nb,4);%bus_res矩阵储存着节点计算结果bus_res(:,1:2)=bus(:,1:2);bus_res(:,3)=bus(:,3)*180/pi;%相角采用角度制bus_res(:,4)=bus(:,4)+(sqrt(-1))*bus(:,5);%注入功率end"S_res.m"计算并返回线路潮流functionS_res=S_res(bus,line,YtYm)[nl,ml]=size(line);S_res=zeros(nl,5);%S_res矩阵储存着线路潮流计算结果S_res(:,1:2)=line(:,1:2);%前两列为节点编号fork=1:nlI=S_res(k,1);J=S_res(k,2);if(J~=0)&(I~=0)S_res(k,3)=bus(I,2)^2*(conj(YtYm(k,3))+conj(YtYm(k,4)))-bus(I,2)*bus(J,2)*(cos(bus(I,3))+j*sin(bus(I,3)))*(conj(cos(bus(J,3))+j*sin(bus(J,3))))*conj(YtYm(k,3));S_res(k,4)=bus(J,2)^2*(conj(YtYm(k,3))+conj(YtYm(k,5)))-bus(I,2)*bus(J,2)*(conj(cos(bus(I,3))+j*sin(bus(I,3))))*(cos(bus(J,3))+j*sin(bus(J,3)))*conj(YtYm(k,3));S_res(k,5)=S_res(k,3)+S_res(k,4);%利用公式计算非接地支路的潮流elseif(J==0)S_res(k,3)=bus(I,2)^2*conj(YtYm(k,4));S_res(k,5)=S_res(k,3)+S_res(k,4);elseS_res(k,4)=bus(J,2)^2*conj(YtYm(k,5));S_res(k,5)=S_res(k,3)+S_res(k,4);%利用公式计算接地支路的潮流endendendend作业y=pifun(10000)y=3.1415>>clearall>>tic;>>y=pifun(100000);>>toc;Elapsedtimeis0.014727seconds.functiony=pifun(n)y=0;fori=1:ny=y+(-1)^(i-1)*1/(2*i-1);endy=y*4;end>>a=0:pi/12:2*pi;>>x=5*cos(a);>>y=3*sin(a);>>plot(x,y)clearformatlongmax1=100;eps1=1.0e-10;eps2=1.0e-10;>>x0=[1.01.01.0]';x=x0;fori=1:max1,A=Fd(x);b=F(x);dx=A\b;%矩阵左除，即b除以Ax=x-dx;sprintf('第%d次迭代结果：\n',i)[xdxF(x)]%在屏幕上输出每次的x(i),dx(i),F(x(i))if(max(abs(dx))<eps1)&(max(abs(F(x)))<eps2)breakendendans=第1次迭代结果：ans=0.9271660275216500.0728339724783501.3111821259085110.4608255403622090.539174459637791-24.062491382655779-0.5300540795902311.5300540795902310.147573280945995ans=第2次迭代结果：ans=0.5018641043912150.4253019231304350.0104813987548680.1873147780763010.273510762285908-5.878576435100477-0.528121534118709-0.0019325454715220.013528291912525ans=第3次迭代结果：ans=0.5005844096414050.0012796947498100.0022669049922680.0610447317778920.126270046298410-1.291470136051974-0.525085549559423-0.0030359845592860.002837480507536ans=第4次迭代结果：ans=0.5001071124267420.0004772972146630.0003396889614650.0115643149629220.049480416814970-0.198312653023639-0.523882243946242-0.0012033056131810.000407361167536ans=第5次迭代结果：ans=0.5000054898064050.0001016226203370.0000165185979270.0005989530207930.010965361942129-0.009739343772851-0.523613482331101-0.0002687616151410.000019534404784ans=第6次迭代结果：ans=0.5000000162962940.0000054735101110.0000000488893170.0000017817132920.000597171307501-0.000028885