基于电力系统潮流计算的IEEE14课程设计matlab程序（含牛顿拉夫逊法和PQ分解法进行极坐标分析）
资源内容介绍
基于IEEE14设计的matlab程序代码 clcclear%% 读取数据mpc = IEEE14;%% 初始化%初始节点电压[baseMVA,bus,gen,branch]=deal(mpc.baseMVA,mpc.bus,mpc.gen,mpc.branch);%复制各基准值n=size(bus,1); %总节点数b=0;c=0;m=0;y=zeros(n,n);U=bus(:,8)';%电压初始值由此确定%电压相角cita=bus(:,9)';cita=(deg2rad(cita)); %角度转换成弧度P_load=bus(:,3);Q_load=bus(:,4);S_load=zeros(n,1);S=zeros(n,1);S_gen=zeros(n,1);P_gen=zeros(n,1);Q_gen=zeros(n,1);%生成发电机节点初始量for k=1:length(gen(:,1))%计算发电机节点数 P_gen(gen(k,1))=gen(k,2);%将P，Q数值赋予各矩阵 Q_gen(gen(k,1))=gen(k,3);endP=zeros(1,n); %P,Q为原始数据，Pi,Qi为计算结果Q=zeros(1,n);for k=1:n P(k)=P_gen(k)-P_load(k);%生成支路功率 Q(k)=Q_gen(k)-Q_load(k);end[P,Q]=deal(P/baseMVA,Q/baseMVA);%多变量赋值%% 求导纳矩阵Ybus=zeros(n,n);Y=creat_Y(mpc);G=real(Y);%G为实部B=imag(Y);%y=zeros(n,n);for i=1:n for j=1:n if i~=j y(i,j)=-Y(i,j); else y(i,j)=sum(Y(i,:)); end endendfor i=1:n if bus(i,2)==2%PV节点数 b=b+1; PV(1,b)=i; endendfor i=1:n if bus(i,2) == 1 m=m+1; %PQ节点数 PQ(1,m)=i; endendfor i=1:n if bus(i,2)==3%平衡节点数 c=c+1; PH=i; end enda1=sortrows(bus,2);%重新进行排序bus=a1;%更新节点顺序U1=bus(:,8)';%求根据PQ节点排序在前的节点矩阵电压幅值dS=[];%定义功率不平衡量show_index=input('是否在命令行展示计算过程？1-是，0-否\n');index=input('请选择潮流计算方法，1-极坐标下牛顿法，2-P-Q分解法\n');%% 迭代求解潮流fprintf('节点数：%d\n',n);disp(n);fprintf('PQ节点数：%d\n',m);disp(PQ);fprintf('PV节点数：%d\n',b);disp(PV);fprintf('PH节点数：%d\n',c);disp(PH);disp(bus);it=1;%初始的迭代次数x=1;indexdiea=input('请输入最大迭代次数\n');it_max=indexdiea;indexepr=input('请输入迭代收敛精度\n');epr=indexepr;%迭代收敛精度%% 计算功率不平衡量[dP,dQ,Pi,Qi]=Unbalanced(n,m,P,Q,U,G,B,cita,bus);if show_index==1 disp('初始条件：');disp('各节点有功：');disp(P); disp('各节点无功：');disp(Q); disp('各节点电压幅值：');disp(U); cita=(deg2rad(cita)); %角度转换成弧度 disp('各节点电压相角(度)：');disp(rad2deg(cita)); %显示依然使用角度 disp('节点导纳矩阵：'); disp(y); disp('有功不平衡量：'); disp(dP); disp('无功不平衡量：'); disp(dQ);endJ=zeros(n+m-1,n+m-1);% index=1;while it<it_max fprintf('第%d次迭代\n',it); %求雅可比矩阵 if index==1 J=Jacobi( n,m,U,U1,cita,B,G,Pi,Qi,bus ); if show_index==1 disp('雅可比矩阵（牛顿拉夫逊法）'); disp(J); end elseif index==2 else disp('输入有误，PQ分解法不需雅可比矩阵') end %求节点电压修正量 if index==1 [dU,dcita]=Correct( n,m,U1,dP,dQ,J ); elseif index==2 [dU,dcita]=PQ_LJ( n,m,dP,dQ,U,U1,B,bus ); else disp('输入有误') end if show_index==1 disp('电压、相角修正量：'); disp(dU); disp(rad2deg(dcita)); end cita=cita+dcita;%更新节点电压相角 for i=1:m U(bus(i,1))=U(bus(i,1))+dU(i);%更新节点电压幅值 end %计算功率不平衡量 [dP,dQ,Pi,Qi]=Unbalanced( n,m,P,Q,U,G,B,cita,bus ); dS(x,1)=max(abs(dP)); dS(x,2)=max(abs(dQ)); if (max(abs(dP))<epr && max(abs(dQ))<epr ) disp('潮流计算收敛'); break end it=it+1; x=x+1;endif(it==it_max&&(max(abs(dP))>epr && max(abs(dQ))>epr)) disp('潮流计算不收敛');elseSij = line_power( n,y,U,cita );%支路功率的计算Pij=real(Sij);Qij=imag(Sij);S=Pi'+1i*Qi';%各节点有功和无功功率Bd=branch(:,5);f=0;h=size(branch,1);Sd=zeros(n,n);for i=1:n%各个节点对地电容功率计算 for j=1:n for p=1:h if (i==branch(p,1)&&j==branch(p,2)) f=f+1; Sd(i,j)=1i*Bd(f,1)*(U(1,i)^2)/2; Sd(j,i)=1i*Bd(f,1)*(U(1,j)^2)/2; end end endendS_load=deal(P_load/baseMVA)+deal(Q_load/baseMVA)*1i;%节点负荷功率的提取k=size(gen,1);for i=1:n%发电机输出功率的提取 for j=1:k if i==gen(j,1) S_gen(i,1)=S(i,1)+S_load(i,1); end endendSij=sparse(Sij);Sij1=Sij-Sd;S1=Sij+conj(Sij');%先对折相加，然后去对角线上的元素S1=triu(S1);%算上对地电容求取的支路功率fprintf('迭代总次数：%d\n', it);disp('节点导纳矩阵：');disp(y);%节点导纳矩阵disp('节点电压幅值：');disp(sparse(U'));%节点电压的幅值disp('节点电压相角：');disp(rad2deg(cita)');%节点电压的相角deg_cita = rad2deg(cita);deg_cita_range = max(deg_cita) - min(deg_cita);%计算极差disp('相角的极差为：');disp((deg_cita_range'));disp('节点注入有功计算结果：');disp(Pi');disp('节点注入无功计算结果：');disp(Qi');Pi_max = max(Pi);Pi_min = min(Pi);disp('有功功率最大值：');disp(Pi_max);disp('有功功率最小值：');disp(Pi_min);Qi_max = max(Qi);Qi_min = min(Qi);disp('无功功率最大值：');disp(Qi_max);disp('无功功率最小值：');disp(Qi_min);disp('各节点负荷功率：');disp(sparse(S_load));disp('各节点发电机功率');disp(sparse(S_gen));disp('支路功率计算结果：');disp(sparse(Sij1))%生成支路功率disp('网络损耗');disp(sparse(S1));disp('对地电容功率');disp(sparse(Sd));subplot(2,3,1);%画节点电压图plot([1:length(U)],U,'.b-','LineWidth',0.8);xlabel('节点号')ylabel('节点电压');subplot(2,3,2);%画节点相角图plot([1:length(deg_cita)],deg_cita,'*m-','LineWidth',0.8);xlabel('节点号')ylabel('节点电压相角');subplot(2,3,3);%画节点有功功率图plot([1:length(Pi)],Pi,'+k-','LineWidth',0.8);xlabel('节点号')ylabel('节点有功功率');subplot(2,3,4);%画节点无功功率的图plot([1:length(Qi)],Qi,'.g-','LineWidth',0.8);xlabel('节点号')ylabel('节点无功功率');subplot(2,3,5);plot([1:length(dS(:,1))],dS(:,1),'.c-','LineWidth',0.001);xlabel('迭代次数')ylabel('有功功率不平衡量');subplot(2,3,6);plot([1:length(dS(:,2))],dS(:,2),'.r-','LineWidth',0.001);xlabel('迭代次数')ylabel('无功功率不平衡量');end