matlab潮流计算

附录一应用牛顿拉夫森法进行潮流计算的Matlab程序代码%牛拉法计算潮流程序%- *% B1行列： 1、分支头编号2、末端编号3 .旁路阻抗4 .旁路对电存% 5、分支路变化率6、旁路的前端k侧为1、1侧为0% B2矩阵： 1、该节点发电功率2、该节点的负载功率3、节点电压初始值% 4、PV节点电压v给定值5、节点连接的无效补偿设备的容量% 6，节点分类标签： 1是平衡节点(必须是第一个节点) 2或PQ节点3是PV节点%- -clear all；format long；n=input (节点数：nodes=)；nl=input (分支数：lines=)；isb=input (平衡母线节点编号：balance=)；pr=input (误差精度：precision=)；B1=input (由各分支参数构成的矩阵：B1=)；B2=input (由各节点参数形成的矩阵：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%的左节点位于1侧p=B1(i，1) q=B1(i，2 )else %左节点位于k侧p=B1(i，2) q=B1(i，1 )结束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 )； %对角元k侧Y(p，p)=Y(p，p) 1./B1(i，3) B1(i，4 )； %对角元1侧结束求%导纳矩阵disp (导入矩阵Y=)；disp(Y )%- -G=real(Y) B=imag(Y) %分解导纳的实部和虚部for i=1:n%给出各节点初始电压的实部和虚部e(i)=real(B2(i，3 ) )；f(i)=imag(B2(i，3 )；V(i)=B2(i，4) %PV节点电压指定模式值结束for i=1:n%给出各节点的注入功率S(i)=B2(i，1)-B2(i，2) %i节点注入功率SG-SLB(i，i)=B(i，i) B2(i，5) %i节点无效补偿量结束%- -P=real(S) Q=imag(S) %分解了各节点注入的有效功率和无效功率ICT1=0； IT2=1； N0=2*n； N=N0 1； a=0； %重复次数ICT1、a； 不满足收敛请求的节点数IT2while IT2=0 % N0=2*n雅可比矩阵的次数N=N0 1扩展列IT2=0； a=a 1；for i=1:nif i=isb%非平衡节点C(i)=0； D(i)=0；for j1=1:nC(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 )结束P1=C(i)*e(i) f(i)*D(i) %节点功率p计算ei(gij * EJ-bij * FJ )-fi(gij * FJ-bij * EJ )Q1=C(i)*f(i)-e(i)*D(i) %节点功率q计算fi(gij * EJ-bij * FJ )-ei(gij * FJ-bij * EJ )求出%I节点的有效功率和无效功率p，q的计算值V2=e(i)2 f(i)2； %电压模式平方%以下对非PV节点求功率差和雅可比矩阵元素- -if B2(i，6)=3%非PV节点DP=P(i)-P1； %节点的有效功率差DQ=Q(i)-Q1； %节点无功功率差%以上为平衡节点以外节点的功率计算- -求%雅各比行列-for j1=1:nif j1=isbj1=i%非平衡节点非对角元素X1=-G(i，j1)*e(i)-B(i，j1)*f(i) % dP/de=-dQ/dfX2=B(i，j1)*e(i)-G(i，j1)*f(i) % dP/df=dQ/deX3=X2； % X2=dp/df X3=dQ/deX4=-X1； % X1=dP/de X4=dQ/dfp=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节点无功功率差J(m，q)=X1； J(m，N)=DP； q=q 1； % X1=dP/de J(m，N)=DP节点的有效功率差J(p，q)=X4； J(m，q)=X2； % X4=dQ/df X2=dp/dfelseif j1=ij1=isb%非平衡节点对角元素X1=-C(i)-G(i，i)*e(i)-B(i，i)*f(i) % dP/deX2=-D(i) B(i，i)*e(i)-G(i，i)*f(i) % dP/dfX3=D(i) B(i，i)*e(i)-G(i，i)*f(i) % dQ/deX4=-C(i) G(i，i)*e(i) B(i，i)*f(i) % dQ/dfp=2*i-1； q=2*j1-1； J(p，q)=X3； J(p，N)=DQ； %扩展列qm=p 1；J(m，q)=X1； q=q 1； J(p，q)=X4； J(m，N)=DP； %扩展列pJ(m，q)=X2；结束结束else%以下是对PV节点求雅可比矩阵的元素-DP=P(i)-P1； % PV节点有效错误DV=V(i)2-V2； % PV节点电压误差for j1=1:nif j1=isbj1=i%非平衡节点非对角元素X1=-G(i，j1)*e(i)-B(i，j1)*f(i) % dP/deX2=B(i，j1)*e(i)-G(i，j1)*f(i) % dP/dfX5=0； X6=0；p=2*i-1； q=2*j1-1； J(p，q)=X5； J(p，N)=DV； % PV节点电压误差m=p 1；J(m，q)=X1； J(m，N)=DP； q=q 1； J(p，q)=X6； % PV节点有效错误J(m，q)=X2；elseif j1=ij1=isb%非平衡节点对角元素X1=-C(i)-G(i，i)*e(i)-B(i，i)*f(i) % dP/deX2=-D(i) B(i，i)*e(i)-G(i，i)*f(i) % dP/dfX5=-2*e(i )X6=-2*f(i )p=2*i-1； q=2*j1-1； J(p，q)=X5； J(p，N)=DV； % PV节点电压误差m=p 1；J(m，q)=X1； J(m，N)=DP； q=q 1； J(p，q)=X6； % PV节点有效错误J(m，q)=X2；结束结束结束结束结束%以上是雅可比矩阵各要素与扩展列的功率差或电压差-for k=3:N0 % N0=2*n从第n行开始，第1行、第2行为平衡节点)k1=k 1； N1=N； % N=N0 1即N=2*n 1扩展列p、q或u从for k2=k1:N1%列雅可比元素到扩展列p、q或uJ(k，k2)=J(k，k2)./J(k，k) %k行k列对角元素规范化除k行k列以外的非对角元素结束J(k，k)=1； %对角要素正规化k行k列对角要素代入1%迭代运算-if k=3 %不是第3行k 3k4=k-1；for k3=3:k4%从第3行到当前行之前的行k4行删除for k2=k1:N1% k3行后各行上三角元素J(k3，k2)=J(k3，k2)-J(k3，k)*J(k，k2) %消去运算(将当前行的k列要素设为0 )end %将当前行K2列元素中的当前行k列元素乘以k行K2列元素J(k3，k)=0； %当前行中第k列的元素变为0结束如果if k=N0 %来到最后一行break；结束%上一代运算- -for k3=k 1:从n