PQ分解法潮流计算程序(最新整理)

1、% 程序功能:14 节点潮流计算 pq 分解法clc;% node% node衡节点infotypepq:1pv:2ph:3所有发电机节点均被视为 pv 节点，节点 1 为平node=% typevangelpgqgplqlshuntc（节点并联电容）31.060000000021.0450040.0000021.700012.7000021.010000094.200019.0000011.000000047.8000-3.9000011.00000007.60001.6000021.070000011.20007.5000011.000000000021.090000000011.0000

2、00029.500016.60000.190011.00000009.00005.8000011.00000003.50001.8000011.00000006.10001.6000011.000000013.50005.8000011.000000014.90005.00000% 功率为有名值，基准功率为 100mva% 线路对地电纳参数为 b，不是 b/2%branchidatajieee14noderxbk（变比）branch=21.01938.05917.0264*21.000;32.04699.19797.0219*21.000;42.05811.17632.0187*21.000;

3、51.05403.22304.0246*21.000;52.05695.17388.0170*21.000;43.06701.17103.0173*21.000;54.01335.04211.0064*21.000;65.00000.25202.00000.932;%变压器支路74.00000.20912.00000.978;%变压器支路87.00000.17615.00001.000;94.00000.55618.00000.969;%变压器支路97.00000.11001.00001.000;109.03181.08450.00001.000;116.09498.19890.00001.0

4、00;126.12291.15581.00001.000;136.06615.13027.00001.000;149.12711.27038.00001.000;1110.08205.19207.00001.000;1312.22092.19988.00001.000;1413.17093.34802.00001.000;x = size(node); %x=14,8 size 矩阵的行数和列数n = x(1);m = length(find(node(:,1)=1); %length 列数 y=x(:,i)的作用就就是把矩阵 x 的第 i 列取出来赋值到 y 里面，那么得到的 y 就是一个列

5、向量。 find 稀疏矩阵非 0 值的个数。y = size(branch);%y=20,6g = zeros(n);%zeros 零矩阵的生成，n 阶方阵b = zeros(n); for i=1:y(1)st= branch(i,1);%起始节点ed= branch(i,2);%终止节点r= branch(i,3);x= branch(i,4);bb= branch(i,5);k= 1/branch(i,6);%变比理想变压器变换为 形等值电路时，线路导纳要除以 k，对地导纳除以 k 或者 k2，该处已经形成倒数。rx2 = r2+x2;g(st,ed) = -k*r/rx2; g(ed,

6、st) = -k*r/rx2; b(st,ed) = k*x/rx2; b(ed,st) = k*x/rx2;g(st,st) = g(st,st) + r/rx2; b(st,st) = b(st,st) - x/rx2;b(st,st) = b(st,st) + bb/2;g(ed,ed) = g(ed,ed) + k2*r/rx2;%根据变压器等值变换电路，系数应该为（1-k)/k2 or(k-1)/k， 这 里 只 有 k2？ b(ed,ed) = b(ed,ed) - k2*x/rx2; b(ed,ed) = b(ed,ed) + bb/2;end;%形成节点导纳矩阵for i=1:

7、nb(i,i) = b(i,i) + node(i,8);%node（i,8)为节点并联电容end;delta_p= zeros(n,1);delta_q= zeros(n,1); delta_v_ample = zeros(n,1);%电压幅值delta_v_angle = zeros(n,1);%电压相角pp = find(node(:,1)3);% pq/pv 节点数（13）qq = find(node(:,1)=1); % pq 节点数（9）%v_ample = node(:,2);%v_angle = node(:,3)*pi/180;v_ample= node(:,2);%各节点电

8、压的初始值，列向量v_ample(qq) = 1.0;v_angle= zeros(n,1);%n 为 14 的全 0 列向量pis = node(:,4)/100; %等值电源有功qis = node(:,5)/100; %等值电源无功pld = node(:,6)/100; %等值负荷有功qld = node(:,7)/100; %等值负荷无功pis(pp) = pis(pp) - pld(pp);%pq/pv 节点注入有功功率qis(qq) = qis(qq) - qld(qq);%pq 节点注入无功功率kp = 0;kq = 0;ep = 1e-7;%误差k= 0;while(kp*k

9、q) %“”是取非符号，while 检验条件，不为 0 是真，执行循环体，0 为假，跳出 while 循环。该语句表达 kp 和 kq 中有一个等于 0 就执行，两个都不等于 0 就结束循环k =k + 1;%-delta_v_angle_ij = v_angle*ones(1,n) - ( v_angle*ones(1,n) );% 计算相角差矩阵，ones 是生成全 1 阵，n 为 14 的行向量s= g.*cos(delta_v_angle_ij) + b.*sin(delta_v_angle_ij); % sij = gij *cos(delta_ij) + bij *sin(delt

10、a_ij)p= v_ample.*(s*v_ample);% p为计算出来的母线功率，“.*用法举例：a=1 2 3,b=4 5 6，a*b=错误结果，因为不符合矩阵相乘的条件，a.*b=1*4 2*5 3*6=4 10 18delta_p(pp)= pis(pp) - p(pp);% 有功不平衡量%pis(pp), p(pp),pis(pp) - p(pp),(pis(pp) - p(pp)./pis(pp)*100if max(abs(delta_p)ep% 有功不平衡量小于允许范围就做标记，abs 数值的绝对值与复数的幅值kp = 1;end;if kp=0% 有功不平衡量大于允许范围就

11、作修正delta_v_angle(pp)=-(b(pp,pp)(delta_p(pp)./v_ample(pp)./v_ample(pp); % 解修正方程，例如 a 的逆矩阵是 a1，则 b/a = b*a1，ab = a1*b， 如果是数字计算那就是等效的v_angle(pp)= v_angle(pp) + delta_v_angle(pp);% 修正相角kq = 0;% 这是必须的，认为无功不平衡量还是不满足end;%-delta_v_angle_ij = v_angle*ones(1,n) - ( v_angle*ones(1,n) );% 计算相角差矩阵t= g.*sin(delta

12、_v_angle_ij) - b.*cos(delta_v_angle_ij); % tij = gij *sin(delta_ij) - bij *cos(delta_ij)q= v_ample.*(t*v_ample);% q为计算出来的母线功率delta_q(qq)= qis(qq) - q(qq);% 无功不平衡量%qis(qq), q(qq),qis(qq) - q(qq),(qis(qq) - q(qq)./qis(qq)*100if max(abs(delta_q)ep% 无功不平衡量小于允许范围就做标记kq = 1;end;if kq=0% 无功不平衡量大于允许范围就作修正de

13、lta_v_ample(qq) = -b(qq,qq)(delta_q(qq)./v_ample(qq);% 解修正方程正幅值v_ample(qq)= v_ample(qq) + delta_v_ample(qq);% 修kp = 0;% 这是必须的，认为有功不平衡量还是不满足end;%-delta_p(pp),delta_q(qq);v_ample(qq),v_angle(pp)*180/pi;%每次迭代结果end; v_ample;node(:,2);v_ample-node(:,2);v_angle*180/pi;node(:,3);v_angle*180/pi-node(:,3);p

14、= p+pld;q = q+qld;v = v_ample,v_angle*180/pi,p,q % 母线全部信息，幅值、相角、发电有功、发电无功sij_flow = zeros(y(1),6); % y(1)为支路总数,sij_flow 为所有支路上流动的功率sij_loss = zeros(y(1),4); % y(1)为支路总数,sij_loss 为所有支路上损耗的功率for i = 1:y(1)ii = branch(i,1); jj = branch(i,2); bb = branch(i,5); kk = branch(i,6);vi = v_ample(ii);vj = v_am

15、ple(jj);ag = delta_v_angle_ij(ii,jj);% ag 为 i,j 节点相角差sij_flow(i,1)=ii;sij_flow(i,2)=jj;sij_flow(i,3)=-g(ii,jj)*(vj2-vi*vj*cos(ag) + vi*vj*b(ii,jj)*sin(-ag);sij_flow(i,4)=b(ii,jj)*(vj2-vi*vj*cos(ag) + vi*vj*g(ii,jj)*sin(-ag)-vj2*bb/2;sij_flow(i,5)=g(ii,jj)*(vi2-vi*vj*cos(ag) - vi*vj*b(ii,jj)*sin(ag);

16、sij_flow(i,6)=-b(ii,jj)*(vi2-vi*vj*cos(ag)-vi*vj*g(ii,jj)*sin(ag)+vi2*bb/2;sij_loss(i,1) =ii;sij_loss(i,2) =jj;sij_loss(i,3) =sij_flow(i,3) - sij_flow(i,5); sij_loss(i,4) =sij_flow(i,4) - sij_flow(i,6);end;s=sij_flow,sij_loss(:,3:4)%全部线路的功率ploss_sum_line = sum(s(:,7); % 合计每条线路上的有功损耗求网损 （s 阵说明：第一、二列为

17、节点序号，第三四列和第五六列分别为支路的首端和末端功率，第七列为该支路有功损耗，第八列为该支路无功损耗）ploss_sum_geld = sum(p)-sum(pld) % 总的母线注入功率和总的负荷功率之差求网损fid=fopen(潮流计算结果.txt,wt);fprintf(fid,节点结果:);fprintf(fid,n);fprintf(fid,电压幅值相角发电有功发电无功); fprintf(fid,n);fprintf(fid,%ft%ft%ft%fn,v);fprintf(fid,线路功率结果:);fprintf(fid,n);fprintf(fid,节点 i节点 jpijqij

18、pji qjiline_p_lossline_q_loss);fprintf(fid,n);fprintf(fid,%dt%dt%ft%ft%ft%ft%ft%fn,s); fprintf(fid,网损 ploss_sum_geld =%f,ploss_sum_geld);fclose(fid)“”“”at the end, xiao bian gives you a passage. minand once said, people who learn to learn are very happy people. in every wonderful life, learning is an eternal theme. as a professional clerical and teaching position, i understand the importance of continuous learning, life is diligent, nothing can be gained, only continuous learning can achieve better self. only by constantly learning and mastering the la