基于电力系统潮流计算的IEEE14课程设计matlab程序（含牛顿拉夫逊法和PQ分解法进行极坐标分析）

clc clear %% 读取数据 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); end P=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 end end for i=1:n if bus(i,2)==2%PV节点数 b=b+1; PV(1,b)=i; end end for i=1:n if bus(i,2) == 1 m=m+1; %PQ节点数 PQ(1,m)=i; end end for i=1:n if bus(i,2)==3%平衡节点数 c=c+1; PH=i; end end a1=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); end J=zeros(n+m-1,n+m-1); % index=1; while itepr && max(abs(dQ))>epr)) disp('潮流计算不收敛'); else Sij = 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 end end S_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 end end Sij=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