• cxhhh
    了解作者
  • matlab
    开发工具
  • 4KB
    文件大小
  • rar
    文件格式
  • 0
    收藏次数
  • 10 积分
    下载积分
  • 3
    下载次数
  • 2018-11-27 19:55
    上传日期
利用牛顿拉夫逊法进行10节点10支路的潮流计算
牛顿拉夫逊法潮流计算.rar
  • Untitled5.m
    12.5KB
内容介绍
function [ output_args ] = Untitled5( input_args ) %UNTITLED5 此处显示有关此函数的摘要 % 此处显示详细说明 %本程序的功能是用牛顿——拉夫逊法进行潮流计算 % B1矩阵:1、支路首端号; 2、末端号; 3、支路阻抗; 4、线路对地电纳 (或变压器导纳); % 5、支路的变比; 6、支路首端处于K侧为1,1侧为0; % 7、线路/变压器标识(0/1)变压器参数当支路首端处于K侧标识为1时归算至末端侧,0归算至首端侧 % B2矩阵:1、该节点发电机功率; 2、该节点负荷功率; % 3、PQ节点电压初始值,或平衡节点及PV节点电压的给定值 % 4、节点所接无功补偿并联电容(感)的电纳 % 5、节点分类标号:1为平衡节点(应为1号节点);2为PQ节点;3为PV节点; clear; isb=1; %input('请输入平衡母线节点号:isb='); pr=1e-5; %input('请输入误差精度:pr='); %--------------------------------------------------- n=10;%input('请输入节点数:n='); nl=10;%input('请输入支路数:nl='); B1=[1 2 3.4+12.8i 1.4e-4i 1 0 0; 1 4 5.1+19.2i 2.1e-4i 1 0 0; 2 3 4.25+16i 1.75e-4i 1 0 0; 4 5 4.25+16i 7e-4i 1 0 0; 1 3 5.1+19.2i 2.1e-4i 1 0 0; 6 4 5.95+22.4i 9.8e-4i 1 0 0; 2 7 1.78+53.89i 0 38.5/231 0 1; 3 8 1.49+48.02i 0 11/231 0 1; 4 9 1.49+48.02i 0 11/231 0 1; 5 10 2.46+70.17i 0 38.5/231 0 1]; %input('请输入由支路参数形成的矩阵: B1='); B2=[0 0 225.5 0 1; 0 0 220 0 2; 0 0 220 0 2; 0 0 220 0 2; 0 0 220 0 2; 120 0 231 0 3; 0 61.11+37.87i 35 0 2; 0 47.53+29.46i 10 0 2; 0 54.32+33.66i 10 0 2; 0 40.74+25.25i 35 0 2]; %input('请输入各节点参数形成的矩阵: B2='); %------------------------------------------------------------- %n=4;%input('请输入节点数:n='); %nl=4;%input('请输入支路数:nl='); %B1=[1 2 4+16i 0 1 0 0; %1 3 4+16i 0 1 0 0; %2 3 2+8i 0 1 0 0; %2 4 1.49+48.02i 0 11/110 0 1] %input('请输入由支路参数形成的矩阵: B1='); %B2=[0 0 115 0 1; %0 0 110 0 2; %0 20+4i 110 0 2; %0 10+6i 10 0 2] %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 %从1到n1(总支路数) if B1(i,7)==1 %-----------如果是变压器支路-------- if B1(i,6)==0 %左节点(首端)处于1侧 p=B1(i,1);q=B1(i,2); else %左节点(首端)处于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); %对角元K侧 Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4); %对角元1侧+励磁导纳 else %------------否则为线路支路-------------------- p=B1(i,1);q=B1(i,2); Y(p,q)=Y(p,q)-1./B1(i,3); %非对角元 Y(q,p)=Y(p,q); %非对角元 Y(q,q)=Y(q,q)+1./B1(i,3)+B1(i,4)./2.0; %对角元j侧+线路电纳的一半 Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2.0; %对角元i侧+线路电纳的一半 end end 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)=abs(B2(i,3)); %PV、平衡节点及PQ节点电压模值 end for i=1:n %给定各节点注入功率 S(i)=B2(i,1)-B2(i,2); %i节点注入功率SG-SL B(i,i)=B(i,i)+B2(i,4); %i节点无功补偿量(电纳值) end %==================用牛顿-拉夫逊法迭代求解非线性代数方程(功率方程)======================= P=real(S);Q=imag(S); %分解出各节点注入的有功和无功功率 ICT1=0;IT2=1;N0=2*n;N1=N0+1;a=0; %迭代次数ICT1、a;不满足收敛要求的节点数IT2 while IT2~=0 % N0=2*n 雅可比矩阵的阶数;N1=N0+1扩展列 IT2=0;a=a+1; JZ=['Jacobi矩阵第(',num2str(a),')次消去运算'];JZ1=['Jacobi矩阵第(',num2str(a),')次回代运算'];JZ0=['功率方程第(',num2str(a),')次差值:']; %----------------求取各个节点的功率及功率偏差及PV节点的电压偏差-------------------- for i=1:n %n个节点2n行(每节点两个方程P和Q或U) p=2*i-1;m=p+1;C(i)=0;D(i)=0; for j1=1:n %第i行共n列(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 %求i节点有功和无功功率P',Q'的计算值 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) V2=e(i)^2+f(i)^2; %电压模平方 %===求取功率差及PV节点电压模平方差 ========= if i~=isb %非平衡节点(PQ或PV节点) if B2(i,5)~=3 %非PV节点(只能是PQ节点) J(m,N1)=P(i)-P1; %PQ节点有功功率差J(m,N1)扩展列△P J(p,N1)=Q(i)-Q1; %PQ节点无功功率差J(p,N1)扩展列△Q else %PV节点================== J(m,N1)=P(i)-P1; %PV节点有功功率差J(m,N1)扩展列△P J(p,N1)=V(i)^2-V2; %PV节点电压模平方差J(p,N1)扩展列△U end end %(if i~=isb) 非平衡节点(PQ或PV节点) end %(for i=1:n) n个节点2n行(每节点两个方程P和Q或U) for m=1:N0 JJN1(m)=J(m,N1); end disp(JZ0);disp(JJN1); %-------------判断功率偏差量及PV节点的电压偏差量是否满足要求----------------- for k=3:N0 %除去平衡节点1、2号以外的所有节点 DET=abs(J(k,N1)); if DET>=pr %PQ节点的功率偏差量及PV节点的电压偏差量是否满足要求 IT2=IT2+1; %不满足要求的节点数加1 end end ICT2(a)=IT2; %不满足要求的节点数;a为迭代次数 ICT1=ICT1+1; %迭代次数 if ICT2(a)==0; %当前不满足要求的节点数为零 break %退出迭代运算 end %--------------------以上为求取各个节点的功率及功率偏差及PV节点的电压偏差------------- %================= 求取Jacobi矩阵形成修正方程 =================== for i=2:n %n个节点2n行(每节点两个方程P和Q或U) if i~=isb %非平衡节点(PQ或PV节点) if B2(i,5)~=3 %下面是针对PQ节点来求取Jacobi矩阵的元素 =========== C(i)=0;D(i)=0; for j1=1:n %第i行共n列(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 for j1=2:n %第i行共n列(2n个Jacobi矩阵元素dP/de及dP/df或dQ/de及dQ/df) if j1~=isb&j1~=i %非平衡节点&非对角元 X1=-G(i,j1)*e(i)-B(i,j1)*f(i); % X1=dP/de=-dQ/df=-X4 X2=B(i,j1)*e(i)-G(i,j1)*f(i); % X2=dP/df=dQ/de=X3 X3=X2; % X2=dp/df X3=dQ/de X4=-X1; % X1=dP/de X4=dQ/df p=2*i-1;q=2*j1-1; J(p,q)=X3;m=p+1; % X3=dQ/de J(p,N)=DQ节点无功功率差 J(p,N)=DQ; J(m,q)=X1;q=q+1; % X1=dP/de J(m,N)=DP节点有功功率差 J(m,N)=DP; J(p,q)=X4;J(m,q)=X2; % X4=dQ/df X2=dp/df elseif j1==i&j1~=isb %非平衡节点&对角元 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/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;%扩展列△Q J(p,N)=DQ; m=p+1; J(m,q)=X1;q=q+1;J(p,q)=X4;%扩展列△P J(m,N)=DP; J(m,q)=X2; end end else %if B2(i,5)~=3 否则(即为PV节点) %=============== 下面是针对PV节点来求取Jacobi矩阵的元素 =========== for j1=1:n if j1~=isb&j1~=i %非平衡节点&非对角元 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; % PV节点电压误差J(p,N)=DV; m=p+1; J(m,q)=X1;q=q+1;J(p,q)=X6; % PV节点有功误差J(m,N)=DP; J(m,q)=X2; elseif j1==i&j1~=isb %非平衡节点&对角元 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/df X5=-2*e(i); X6=-2*f(i); p=2*i-1;q=2*j1-1;J(p,q)=X5; % PV节点电压误差J(p,N)=DV; m=p+1; J(m,q)=X1;q=q+1;J(p,q)=X6; % PV节点有功误差J(m,N)=DP; J(m,q)=X2; end end end %(if B2(i,5)~=3 else) end %(if i~=isb) end %(for i=1:n)n个节点2n行(每节点两个方程P和Q或U) JZ0=['形成的第(',num2str(a),')次Jacobi矩阵:']; %disp(JZ0);disp(J); %=============================== 以上为形成完整的Jacobi矩阵 ================================ %====下面用高斯消去法对由Jacobi矩阵形成的修正方程进行求解(按列消去、回代) ========== for k=3:N0 % N0=2*n (从第三行开始,第一、二行是平衡节点) for k1=k+1:N1 % 从k+1列的Jacobi元素到扩展列的△P、△Q 或 △U J(k,k1)=J(k,k1)./J(k,k);% 用K行K列对角元素去除K行K列后的非对角元素进行规格化 end J(k,k)=1; % 对角元规格化K行K列对角元素赋1 %================== 按列消去运算 ================================== for k2=k+1:N0 % 从k+1行到2*n最后一行 for k3=k+1:N1 %从k2+1列到扩展列消去k+1行后各行下三角元素 J(k2,k3)=J(k2,k3)-J(k2,k)*J(k,k3);%消去运算 end %用当前行K3列元素减去当前行K列元素乘以第k行K3列元素 J(k2,k)=0; %当前行第k列元素已消为0 end end %JZ=['Jacobi矩阵第(',num2str(a),')次消去运算'];JZ1=['Jacobi矩阵第(',num2str(a),')次回代运算']
评论
    相关推荐