電力網絡潮流計算

TJU電力系統分析1課程編程大作業
以下程序由自己編寫，如有錯誤還望指出糾正

題目要求

上圖中的公式均為課本中相應的公式




A

程序

牛頓拉夫遜法

程序框圖

一共包含三個文件

實現方法如下

%system_branch_data.m %main.m %Jacobi.m %system_branch_data.m 1 4 0 0.0576 0 2 7 0 0.0625 0 3 9 0 0.0586 0 4 5 0.01 0.085 0.088 4 6 0.017 0.092 0.079 5 7 0.032 0.161 0.153 6 9 0.039 0.17 0.179 7 8 0.0085 0.072 0.0745 8 9 0.0119 0.1008 0.1045 %main.m %極坐標系下的牛頓-拉夫遜潮流算法 clear; clc tic; net_info = importdata("system_branch_data.m"); [row, column] = size(net_info); f_branch = net_info(:,1);%支路首節點 e_branch = net_info(:,2);%支路末節點 r_branch = net_info(:,3);%支路電阻 x_branch = net_info(:,4);%支路電抗 Yn = zeros(row); %導納矩陣 for i = 1:rownode_i = net_info(i, 1); node_j = net_info(i, 2);Yn(node_i, node_i) = Yn(node_i, node_i) + 1 / (net_info(i, 3) + 1j * net_info(i, 4)) + 1j * net_info(i, 5);Yn(node_j, node_j) = Yn(node_j, node_j) + 1 / (net_info(i, 3) + 1j * net_info(i, 4)) + 1j * net_info(i, 5);Yn(node_i, node_j) = -1 / (net_info(i, 3) + 1j * net_info(i, 4));Yn(node_j, node_i) = -1 / (net_info(i, 3) + 1j * net_info(i, 4)); end disp('節點導納矩陣（非零元，一行一行地輸出）') for i = 1:rowfor j = i:rowmat = Yn(i,j);if(mat~=0)matendend end%初值 %按照平衡節點，PV節點，PQ節點順序依次賦予初值 global Num_PV; global Num_PQ; Num_PV=2; Num_PQ=6; U_init = [1.04, 1.025, 1.025, 1, 1, 1, 1, 1, 1]'; Theta_init = [0,0,0,0,0,0,0,0,0]'; P_init = [1.63, 0.85, 0, -1.25, -0.9, 0, -1, 0]'; Q_init = [0, -0.5, -0.3, 0, -0.35, 0]'; P_node = [0,0,0,0,1.25,0.9,0,1.0,0]'; Q_node = [0,0,0,0,0.5,0.3,0,0.35,0]'; global n; global m; n=Num_PQ+Num_PV;%n=8 m=Num_PQ;%m=6 Theta_ij=zeros(n+1,n+1); G=real(Yn); B=imag(Yn); cnt=0; cntmax=10; Test=ones(n+m,1);while cnt<=cntmax && max(abs(Test))>0.00001U_temp1=zeros(n+1,1);U_temp2=zeros(n+1,1);for i=1:n+1for j=1:n+1Theta_ij(i,j)=Theta_init(i)-Theta_init(j);U_temp1(i)=U_temp1(i)+U_init(j)*(G(i,j)*sin(Theta_ij(i,j))-B(i,j)*cos(Theta_ij(i,j)));U_temp2(i)=U_temp2(i)+U_init(j)*(G(i,j)*cos(Theta_ij(i,j))+B(i,j)*sin(Theta_ij(i,j)));endendDelta_P=P_init-U_init(2:n+1).*U_temp2(2:n+1);Delta_Q=Q_init-U_init(2+Num_PV:n+1).*U_temp1(2+Num_PV:n+1);P_now=U_init(1:n+1).*U_temp2(1:n+1);Q_now=U_init(1:n+1).*U_temp1(1:n+1); J = Jacobi(P_now,Q_now,Yn,U_init,Theta_init);J_inv=inv(J);Temp=-J_inv*[Delta_P;Delta_Q];Delta_theta=Temp(1:n);Delta_U=(diag(U_init(1:m))*Temp(n+1:m+n));Theta_init(2:n+1)=Theta_init(2:n+1)+Delta_theta;U_init(2+Num_PV:n+1)=U_init(2+Num_PV:n+1)+Delta_U.*U_init(2+Num_PV:n+1); Test=[Delta_P;Delta_Q];cnt=cnt+1; end if cnt>cntmaxmsgbox('無收斂解');elsedisp("節點電壓幅值")U_initdisp("節點電壓相角")Theta_init*180/pi endTheta1 = Theta_init*180/pi; Theta_final=Theta_init; U_final=U_init.*cos(Theta_final)+1j*U_init.*sin(Theta_final); I_final=Yn*U_final; S_final=U_final.*conj(I_final); P_final=real(S_final); Q_final=imag(S_final); Y=Yn; for i=1:Num_PV+1Y(i,i)=Y(i,i)+1/(0.3*1j); end %第二個輸出：增加發電機導納后的自導納 disp("增加發電機導納后的自導納"); for i=1:rowY(i,i) end Z1=inv(Y); for i=Num_PV+2:n+1Y(i,i)=Y(i,i)+(P_node(i)-1j*Q_node(i))/(U_init(i)^2); end %第三個輸出：增加負荷節點導納后的自導納 disp("增加負荷節點導納后的自導納"); for i=1:rowY(i,i) endZ=inv(Y);%第四個輸出：節點阻抗矩陣第四列 disp("節點阻抗矩陣第四列"); disp(Z(:,4));%節點4精確短路電流 I_f = abs(U_init(4)./Z(4,4)); %第五個輸出：精確短路電流幅值 disp("精確短路電流幅值"); I_f %第六個輸出：精確短路電流相角 disp("精確短路電流相角"); I_f_Theta = (Theta_final(4)-angle(Z(4,4))); I_f_Theta*180/pi%近似計算 I_f1 = 1/abs(Z1(4,4)); %第七個輸出：近似短路電流幅值 disp("近似短路電流幅值"); I_f1 I_f1_Theta = (-angle(Z1(4,4))); %第八個輸出：近似短路電流相角 disp("近似短路電流相角"); I_f1_Theta*180/pi%精確與近似誤差比較 I_f_errAmp = abs((I_f-I_f1))/I_f; %第九個輸出：電流幅值誤差 disp("電流幅值誤差"); I_f_errAmp I_f_errAng = abs((I_f_Theta-I_f1_Theta)/I_f_Theta); %第十個輸出：電流相角誤差 disp("電流相角誤差"); I_f_errAng%短路后各點電壓for i = 1:rowU_f(i)=U_init(i)*(cos(Theta_final(i))+1j*sin(Theta_final(i)))-Z(i,4)*U_init(4)*(cos(Theta_final(4))+1j*sin(Theta_final(4)))/Z(4,4);U_f1(i)=1-Z1(i,4)*(I_f1*cos(I_f1_Theta)+1j*I_f1*sin(I_f1_Theta));end %第十一個輸出：精確計算后的短路后各節點電壓幅值 disp('精確計算后的短路后各節點電壓幅值'); abs(U_f) %第十二個輸出：近似計算后的短路后各節點電壓幅值 disp('近似計算后的短路后各節點電壓幅值'); abs(U_f1) %第十三個輸出：誤差電壓比較for i=1:rowif(abs(U_f(i))>0.0001)U_f_err(i)=abs((abs(U_f(i))-abs(U_f1(i)))/abs(U_f(i)));elseU_f_err(i) = 0.0000;end end disp('誤差電壓比較'); abs(U_f_err)%第十四個輸出：各支路精確電流 disp('各支路精確電流'); I_f_branch = zeros(row); for i=1:rowI_f_branch(f_branch(i),e_branch(i)) = (U_f(f_branch(i))-U_f(e_branch(i)))/((r_branch(i)+1j*x_branch(i)));if(I_f_branch(f_branch(i),e_branch(i))~=0)disp(['I_f_branch_',num2str(f_branch(i)),num2str(e_branch(i))]);abs(I_f_branch(f_branch(i),e_branch(i)))angle(I_f_branch(f_branch(i),e_branch(i)))*180/piend end %Jacobi.m function J = Jacobi(P,Q,Yn,U,Theta) %求解雅可比矩陣 % % Syntax: J = Jacobi(P,Q,Yn,U,Theta) global n global Num_PV G=real(Yn); B=imag(Yn); for i=2:n+1for j=2:n+1if i~=jTheta_ij(i,j)=Theta(i)-Theta(j);H(i-1,j-1)=-U(i)*U(j)*(G(i,j)*sin(Theta_ij(i,j))-B(i,j)*cos(Theta_ij(i,j)));elseH(i-1,j-1)=U(i)*U(j)*B(i,j)+Q(i);endend end for i=2:n+1for j=(2+Num_PV):n+1if i~=jTheta_ij(i,j)=Theta(i)-Theta(j);N(i-1,j-1-Num_PV)=-U(i)*U(j)*(G(i,j)*cos(Theta_ij(i,j))+B(i,j)*sin(Theta_ij(i,j)));elseN(i-1,j-1-Num_PV)=-U(i)*U(j)*G(i,j)-P(i);endend end for i=2+Num_PV:n+1for j=2:n+1if i~=jTheta_ij(i,j)=Theta(i)-Theta(j);K(i-1-Num_PV,j-1)=U(i)*U(j)*(G(i,j)*cos(Theta_ij(i,j))+B(i,j)*sin(Theta_ij(i,j)));elseK(i-1-Num_PV,j-1)=U(i)*U(j)*G(i,j)-P(i);endend end for i=2+Num_PV:n+1for j=2+Num_PV:n+1if i~=jTheta_ij(i,j)=Theta(i)-Theta(j);L(i-1-Num_PV,j-1-Num_PV)=-U(i)*U(j)*(G(i,j)*sin(Theta_ij(i,j))-B(i,j)*cos(Theta_ij(i,j)));elseL(i-1-Num_PV,j-1-Num_PV)=U(i)*U(j)*B(i,j)-Q(i);endend end J=[H,N;K,L];end

快速解耦法

程序框圖

%相應資源后續上傳

總結

以上是生活随笔為你收集整理的编写程序计算算例系统的潮流及三相短路电流(Matlab)的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。