编写程序计算算例系统的潮流及三相短路电流(Matlab)

电力网络潮流计算

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:row
    node_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:row
    for j = i:row
        mat = Yn(i,j);
        if(mat~=0)
            mat
        end
    end
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.00001
    U_temp1=zeros(n+1,1);
    U_temp2=zeros(n+1,1);
    for i=1:n+1
        for j=1:n+1
            Theta_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)));
        end
    end
    Delta_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>cntmax
    msgbox('无收敛解');
 else
     disp("节点电压幅值")
     U_init
     disp("节点电压相角")
     Theta_init*180/pi
end
 
Theta1 = 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+1
    Y(i,i)=Y(i,i)+1/(0.3*1j);
end
%第二个输出：增加发电机导纳后的自导纳
disp("增加发电机导纳后的自导纳");
for i=1:row
    Y(i,i)
end
Z1=inv(Y);
for i=Num_PV+2:n+1
    Y(i,i)=Y(i,i)+(P_node(i)-1j*Q_node(i))/(U_init(i)^2);
end
%第三个输出：增加负荷节点导纳后的自导纳
disp("增加负荷节点导纳后的自导纳");
for i=1:row
    Y(i,i)
end

Z=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:row
    U_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:row
    if(abs(U_f(i))>0.0001)
    U_f_err(i)=abs((abs(U_f(i))-abs(U_f1(i)))/abs(U_f(i)));
    else
        U_f_err(i) = 0.0000;
    end
end
disp('误差电压比较');
abs(U_f_err)

%第十四个输出：各支路精确电流
disp('各支路精确电流');
I_f_branch = zeros(row);
for i=1:row
    I_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/pi
    end
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+1
    for j=2:n+1
        if i~=j
            Theta_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)));
        else
            H(i-1,j-1)=U(i)*U(j)*B(i,j)+Q(i);
        end
    end
end
for i=2:n+1
    for j=(2+Num_PV):n+1
         if i~=j
            Theta_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)));
         else
            N(i-1,j-1-Num_PV)=-U(i)*U(j)*G(i,j)-P(i);
         end
    end
end
for i=2+Num_PV:n+1
    for j=2:n+1
        if i~=j
            Theta_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)));
        else
            K(i-1-Num_PV,j-1)=U(i)*U(j)*G(i,j)-P(i);
        end
    end
end
for i=2+Num_PV:n+1
    for j=2+Num_PV:n+1
        if i~=j
            Theta_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)));
        else
            L(i-1-Num_PV,j-1-Num_PV)=U(i)*U(j)*B(i,j)-Q(i);
        end
    end
end  
J=[H,N;K,L];

end

快速解耦法

程序框图

%相应资源后续上传