本文详细介绍了牛顿法在电力系统潮流计算中的应用，包括算法流程图、计算过程和代码实现。通过输入电网及注入节点数据，计算导纳阵，判断节点类型，计算雅可比矩阵，输出最终结果，展示了完整的代码示例和运行结果。
目录
潮流计算
计算算例
算法流程图
如下：
计算过程 及代码
1、输入电网及注入节点数据
2、计算导纳阵
%% 输入节点之间阻抗和导纳
Y0=zeros(5,15);
Y0=[0,0,0,complex(0,0.25),1/complex(0.04,0.25),complex(0,0.25),0,1/complex(0.1,0.35),0,0,0,0,0,0,0;
complex(0,0.25),1/complex(0.04,0.25),complex(0,0.25),0,0,0,complex(0,0.25),1/complex(0.08,0.3),complex(0,0.25),(1-1.05)/complex(0,1.05*1.05*0.015),1/complex(0,1.05*0.015),(1.05-1)/complex(0,1.05*0.015),0,0,0;
0,1/complex(0.1,0.35),0,complex(0,0.25),1/complex(0.08,0.3),complex(0,0.25),0,0,0,0,0,0,(1-1.05)/complex(0,1.05*1.05*0.03),1/complex(0,1.05*0.03),(1.05-1)/complex(0,1.05*0.03);
0,0,0,(1.05-1)/complex(0,1.05*0.015),1/complex(0,1.05*0.015),(1-1.05)/complex(0,1.05*1.05*0.015),0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,(1.05-1)/complex(0,1.05*0.03),1/complex(0,1.05*0.03),(1-1.05)/complex(0,1.05*1.05*0.03),0,0,0,0,0,0];
%虚数可以直接用2+i*4表示
%% 计算导纳阵Y
N=5;
Y=zeros(N,N);
for i=1:N
for j=1:N
if i==j
Y(i,j)=sum(Y0(i,:))-Y0(i,3)-Y0(i,6)-Y0(i,9)-Y0(i,12)-Y0(i,15);
else
Y(i,j)=-Y0(i,j*3-1);
end
end
end3、判断节点类型
%% 判断节点类型
for i=1:N
if data(i,2)==0 && data(i,3)==0 %如果PG、QG为0，那就是PQ节点
fprintf('%d节点是PQ节点。\n',i);
[p_q_derta(2*i-1),p_q_derta(2*i)]=PQ(p_q_matrix(2*i-1),p_q_matrix(2*i),Y,e_f_matrix,i);
elseif isnan(data(i,2)) && isnan(data(i,3)) %如果PG、QG为NAN，PL、QL为0，那就是V0节点
fprintf('%d节点是V0节点。\n',i);
else
fprintf('%d节点是PV节点。\n',i);
[p_q_derta(2*i-1),p_q_derta(2*i)]=PV(p_v_matrix(2*i-1),p_v_matrix(2*i),Y,e_f_matrix,i);
end
end4、计算雅可比矩阵并输出
%% 计算雅可比矩阵
for j=1:4
for i=1:2*N
syms (['e_f',num2str(i)]);
end
for i=1:N
if data(i,2)==0 && data(i,3)==0
p_derta_dd=p_q_matrix(2*i-1)-...
['e_f',num2str(2*i-1)]*((real(Y(i,1))*e_f1-imag(Y(i,1))*e_f2)+...
(real(Y(i,2))*e_f3-imag(Y(i,2))*e_f4)+...
(real(Y(i,3))*e_f5-imag(Y(i,3))*e_f6)+...
(real(Y(i,4))*e_f7-imag(Y(i,4))*e_f8)+...
(real(Y(i,5))*e_f9-imag(Y(i,5))*e_f10))-...
['e_f',num2str(2*i)]*((real(Y(i,1))*e_f2+imag(Y(i,1))*e_f1)+...
(real(Y(i,2))*e_f4+imag(Y(i,2))*e_f3)+...
(real(Y(i,3))*e_f6+imag(Y(i,3))*e_f5)+...
(real(Y(i,4))*e_f8+imag(Y(i,4))*e_f7)+...
(real(Y(i,5))*e_f10+imag(Y(i,5))*e_f9));
J1=jacobian(p_derta_dd,[e_f1,e_f2,e_f3,e_f4,e_f5,e_f6,e_f7,e_f8,e_f9,e_f10]);
q_derta_dd=p_q_matrix(2*i)-...
['e_f',num2str(2*i)]*((real(Y(i,1))*e_f1-imag(Y(i,1))*e_f2)+...
(real(Y(i,2))*e_f3-imag(Y(i,2))*e_f4)+...
(real(Y(i,3))*e_f5-imag(Y(i,3))*e_f6)+...
(real(Y(i,4))*e_f7-imag(Y(i,4))*e_f8)+...
(real(Y(i,5))*e_f9-imag(Y(i,5))*e_f10))+...
['e_f',num2str(2*i-1)]*((real(Y(i,1))*e_f2+imag(Y(i,1))*e_f1)+...
(real(Y(i,2))*e_f4+imag(Y(i,2))*e_f3)+...
(real(Y(i,3))*e_f6+imag(Y(i,3))*e_f5)+...
(real(Y(i,4))*e_f8+imag(Y(i,4))*e_f7)+...
(real(Y(i,5))*e_f10+imag(Y(i,5))*e_f9));
J2=jacobian(q_derta_dd,[e_f1,e_f2,e_f3,e_f4,e_f5,e_f6,e_f7,e_f8,e_f9,e_f10]);
J0=[J1;J2];
if i==1
J0_0=J0;
else
J0_0=[J0_0;J0];
end
elseif isnan(data(i,2)) && isnan(data(i,3))
J0=zeros(2,10);
if i==1
J0_0=J0;
else
J0_0=[J0_0;J0];
end
else
p_derta_dd=p_q_matrix(2*i-1)-...
['e_f',num2str(2*i-1)]*((real(Y(i,1))*e_f1-imag(Y(i,1))*e_f2)+...
(real(Y(i,2))*e_f3-imag(Y(i,2))*e_f4)+...
(real(Y(i,3))*e_f5-imag(Y(i,3))*e_f6)+...
(real(Y(i,4))*e_f7-imag(Y(i,4))*e_f8)+...
(real(Y(i,5))*e_f9-imag(Y(i,5))*e_f10))-...
['e_f',num2str(2*i)]*((real(Y(i,1))*e_f2+imag(Y(i,1))*e_f1)+...
(real(Y(i,2))*e_f4+imag(Y(i,2))*e_f3)+...
(real(Y(i,3))*e_f6+imag(Y(i,3))*e_f5)+...
(real(Y(i,4))*e_f8+imag(Y(i,4))*e_f7)+...
(real(Y(i,5))*e_f10+imag(Y(i,5))*e_f9));
J1=jacobian(p_derta_dd,[e_f1,e_f2,e_f3,e_f4,e_f5,e_f6,e_f7,e_f8,e_f9,e_f10]);
v_derta_dd=p_v_matrix(2*i)^2-(eval(['e_f',num2str(2*i-1)])^2+eval(['e_f',num2str(2*i)])^2);
J2=jacobian(v_derta_dd,[e_f1,e_f2,e_f3,e_f4,e_f5,e_f6,e_f7,e_f8,e_f9,e_f10]);
J0=[J1;J2];
if i==1
J0_0=J0;
else
J0_0=[J0_0;J0];
end
end
end
for i=1:2*N
eval(['e_f',num2str(i),'=','e_f_matrix(i)',';']);
end
disp('雅可比矩阵为：')
J0_data=vpa(subs(J0_0),5);
disp(J0_data)
%% derta e f ——> e-f新 derta p q ——> p q新
e_f_derta=vpa(linsolve(J0_data,reshape(p_q_derta,10,1)),5);
e_f_matrix=vpa(reshape(e_f_matrix,10,1)-e_f_derta,5);
for w=1:N
if data(w,2)==0 && data(w,3)==0
[p_q_derta(2*w-1),p_q_derta(2*w)]=PQ(p_q_matrix(2*w-1),p_q_matrix(2*w),Y,e_f_matrix,w);
elseif isnan(data(w,2)) && isnan(data(w,3))
else
[p_q_derta(2*w-1),p_q_derta(2*w)]=PV(p_v_matrix(2*w-1),p_v_matrix(2*w),Y,e_f_matrix,w);
end
end
p_q_derta=vpa(p_q_matrix-p_q_derta,5);
p_q_derta=vpa(p_q_matrix-p_q_derta,5);
%% 输出新e-f阵和新p-q阵
disp('新的derta-ef阵')
disp(e_f_derta)
disp('新的ef阵')
disp(e_f_matrix)
disp('新的derta-pd阵')
disp(p_q_derta)
e_f_end=horzcat(e_f_end,e_f_matrix);
p_q_derta_end=horzcat(p_q_derta_end,p_q_derta);
end5、输出最终所有所需矩阵
%% 输出最后的所有矩阵变化情况
disp('最终节点功率误差的变化情况')
p_q_derta_end= reshape(p_q_derta_end,10,5)';
B = round(p_q_derta_end.* 1e5) ./ 1e5; % 将矩阵中的每个元素保留 5 位小数
C = reshape(B, [5 10]); % 将矩阵转换为一个 5 行 10 列的新矩阵
formatSpec = '%10.5f'; % 设置输出格式
for i = 1:5
fprintf([repmat(formatSpec, 1, 10) '\n'], C(i, :)); % 输出第 i 行元素，每行输出 10 个元素，用空格分隔，再换行
end
disp('最终节点电压的变化情况')
e_f_end= reshape(e_f_end,10,4)';
B0 = round(e_f_end.* 1e5) ./ 1e5; % 将矩阵中的每个元素保留 5 位小数
C0 = reshape(B0, [4 10]); % 将矩阵转换为一个 4 行 10 列的新矩阵
formatSpec = '%10.5f'; % 设置输出格式
for i = 1:4
fprintf([repmat(formatSpec, 1, 10) '\n'], C0(i, :)); % 输出第 i 行元素，每行输出 10 个元素，用空格分隔，再换行
end
for m=1:4
z = e_f_matrix(2*m-1) - e_f_matrix(2*m)*1i;
U0_zeita(2*m-1) = abs(z); % 计算幅值
U0_zeita(2*m) = -rad2deg(angle(z)); % 计算幅角（以度为单位）
end
disp("U_zeita矩阵")
disp(U0_zeita)6、用到的两个函数PQ、PV节点
%% 函数 PQ节点
function [P,Q]=PQ(P0,Q0,Y,e_f,i)
P=P0-...
e_f(2*i-1)*((real(Y(i,1))*e_f(1)-imag(Y(i,1))*e_f(2))+...
(real(Y(i,2))*e_f(3)-imag(Y(i,2))*e_f(4))+...
(real(Y(i,3))*e_f(5)-imag(Y(i,3))*e_f(6))+...
(real(Y(i,4))*e_f(7)-imag(Y(i,4))*e_f(8))+...
(real(Y(i,5))*e_f(9)-imag(Y(i,5))*e_f(10)))-...
e_f(2*i)*((real(Y(i,1))*e_f(2)+imag(Y(i,1))*e_f(1))+...
(real(Y(i,2))*e_f(4)+imag(Y(i,2))*e_f(3))+...
(real(Y(i,3))*e_f(6)+imag(Y(i,3))*e_f(5))+...
(real(Y(i,4))*e_f(8)+imag(Y(i,4))*e_f(7))+...
(real(Y(i,5))*e_f(10)+imag(Y(i,5))*e_f(9)));
Q=Q0-...
e_f(2*i)*((real(Y(i,1))*e_f(1)-imag(Y(i,1))*e_f(2))+...
(real(Y(i,2))*e_f(3)-imag(Y(i,2))*e_f(4))+...
(real(Y(i,3))*e_f(5)-imag(Y(i,3))*e_f(6))+...
(real(Y(i,4))*e_f(7)-imag(Y(i,4))*e_f(8))+...
(real(Y(i,5))*e_f(9)-imag(Y(i,5))*e_f(10)))+...
e_f(2*i-1)*((real(Y(i,1))*e_f(2)+imag(Y(i,1))*e_f(1))+...
(real(Y(i,2))*e_f(4)+imag(Y(i,2))*e_f(3))+...
(real(Y(i,3))*e_f(6)+imag(Y(i,3))*e_f(5))+...
(real(Y(i,4))*e_f(8)+imag(Y(i,4))*e_f(7))+...
(real(Y(i,5))*e_f(10)+imag(Y(i,5))*e_f(9)));
end
%% PV节点
function [P,V]=PV(P0,V0,Y,e_f,i)
P=P0-...
e_f(2*i-1)*((real(Y(i,1))*e_f(1)-imag(Y(i,1))*e_f(2))+...
(real(Y(i,2))*e_f(3)-imag(Y(i,2))*e_f(4))+...
(real(Y(i,3))*e_f(5)-imag(Y(i,3))*e_f(6))+...
(real(Y(i,4))*e_f(7)-imag(Y(i,4))*e_f(8))+...
(real(Y(i,5))*e_f(9)-imag(Y(i,5))*e_f(10)))-...
e_f(2*i)*((real(Y(i,1))*e_f(2)+imag(Y(i,1))*e_f(1))+...
(real(Y(i,2))*e_f(4)+imag(Y(i,2))*e_f(3))+...
(real(Y(i,3))*e_f(6)+imag(Y(i,3))*e_f(5))+...
(real(Y(i,4))*e_f(8)+imag(Y(i,4))*e_f(7))+...
(real(Y(i,5))*e_f(10)+imag(Y(i,5))*e_f(9)));
V=V0^2-(e_f(2*i-1)^2+e_f(2*i)^2);
end
最终总代码及演示结果
代码展示
data = xlsread('C:\Users\17860\Desktop\节点数据表.xlsx');
%% 输入节点之间阻抗和导纳
Y0=zeros(5,15);
Y0=[0,0,0,complex(0,0.25),1/complex(0.04,0.25),complex(0,0.25),0,1/complex(0.1,0.35),0,0,0,0,0,0,0;
complex(0,0.25),1/complex(0.04,0.25),complex(0,0.25),0,0,0,complex(0,0.25),1/complex(0.08,0.3),complex(0,0.25),(1-1.05)/complex(0,1.05*1.05*0.015),1/complex(0,1.05*0.015),(1.05-1)/complex(0,1.05*0.015),0,0,0;
0,1/complex(0.1,0.35),0,complex(0,0.25),1/complex(0.08,0.3),complex(0,0.25),0,0,0,0,0,0,(1-1.05)/complex(0,1.05*1.05*0.03),1/complex(0,1.05*0.03),(1.05-1)/complex(0,1.05*0.03);
0,0,0,(1.05-1)/complex(0,1.05*0.015),1/complex(0,1.05*0.015),(1-1.05)/complex(0,1.05*1.05*0.015),0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,(1.05-1)/complex(0,1.05*0.03),1/complex(0,1.05*0.03),(1-1.05)/complex(0,1.05*1.05*0.03),0,0,0,0,0,0];
%虚数可以直接用2+i*4表示
%% 计算导纳阵Y
N=5;
Y=zeros(N,N);
for i=1:N
for j=1:N
if i==j
Y(i,j)=sum(Y0(i,:))-Y0(i,3)-Y0(i,6)-Y0(i,9)-Y0(i,12)-Y0(i,15);
else
Y(i,j)=-Y0(i,j*3-1);
end
end
end
%% 设定e f p q初值
e_f_matrix=[1,0,1,0,1,0,1.05,0,1.05,0];
p_q_matrix=zeros(2,N);
p_v_matrix=zeros(2,N);
p_q_matrix(1,:)=-data(:,4).';
p_q_matrix(2,:)=-data(:,5).';
p_v_matrix(1,:)=data(:,2).';
p_v_matrix(2,:)=data(:,6).';
%初始p_q_matrix
%disp('p_q_matrix')
%% 定义derta [e f p q]
p_q_derta=zeros(2,N);
e_f_derta=zeros(1,2*N);
%% 定义最终输出变化情况的矩阵
p_q_derta_end=p_q_derta;
e_f_end=[];
U0_zeita=[];
%% 判断节点类型
for i=1:N
if data(i,2)==0 && data(i,3)==0 %如果PG、QG为0，那就是PQ节点
fprintf('%d节点是PQ节点。\n',i);
[p_q_derta(2*i-1),p_q_derta(2*i)]=PQ(p_q_matrix(2*i-1),p_q_matrix(2*i),Y,e_f_matrix,i);
elseif isnan(data(i,2)) && isnan(data(i,3)) %如果PG、QG为NAN，PL、QL为0，那就是V0节点
fprintf('%d节点是V0节点。\n',i);
else
fprintf('%d节点是PV节点。\n',i);
[p_q_derta(2*i-1),p_q_derta(2*i)]=PV(p_v_matrix(2*i-1),p_v_matrix(2*i),Y,e_f_matrix,i);
end
end
%% 打印输出第一次的节点功率dertaP、dertaQ
disp(p_q_derta);
p_q_derta_end=p_q_derta;
%% 计算雅可比矩阵
for j=1:4
for i=1:2*N
syms (['e_f',num2str(i)]);
end
for i=1:N
if data(i,2)==0 && data(i,3)==0
p_derta_dd=p_q_matrix(2*i-1)-...
['e_f',num2str(2*i-1)]*((real(Y(i,1))*e_f1-imag(Y(i,1))*e_f2)+...
(real(Y(i,2))*e_f3-imag(Y(i,2))*e_f4)+...
(real(Y(i,3))*e_f5-imag(Y(i,3))*e_f6)+...
(real(Y(i,4))*e_f7-imag(Y(i,4))*e_f8)+...
(real(Y(i,5))*e_f9-imag(Y(i,5))*e_f10))-...
['e_f',num2str(2*i)]*((real(Y(i,1))*e_f2+imag(Y(i,1))*e_f1)+...
(real(Y(i,2))*e_f4+imag(Y(i,2))*e_f3)+...
(real(Y(i,3))*e_f6+imag(Y(i,3))*e_f5)+...
(real(Y(i,4))*e_f8+imag(Y(i,4))*e_f7)+...
(real(Y(i,5))*e_f10+imag(Y(i,5))*e_f9));
J1=jacobian(p_derta_dd,[e_f1,e_f2,e_f3,e_f4,e_f5,e_f6,e_f7,e_f8,e_f9,e_f10]);
q_derta_dd=p_q_matrix(2*i)-...
['e_f',num2str(2*i)]*((real(Y(i,1))*e_f1-imag(Y(i,1))*e_f2)+...
(real(Y(i,2))*e_f3-imag(Y(i,2))*e_f4)+...
(real(Y(i,3))*e_f5-imag(Y(i,3))*e_f6)+...
(real(Y(i,4))*e_f7-imag(Y(i,4))*e_f8)+...
(real(Y(i,5))*e_f9-imag(Y(i,5))*e_f10))+...
['e_f',num2str(2*i-1)]*((real(Y(i,1))*e_f2+imag(Y(i,1))*e_f1)+...
(real(Y(i,2))*e_f4+imag(Y(i,2))*e_f3)+...
(real(Y(i,3))*e_f6+imag(Y(i,3))*e_f5)+...
(real(Y(i,4))*e_f8+imag(Y(i,4))*e_f7)+...
(real(Y(i,5))*e_f10+imag(Y(i,5))*e_f9));
J2=jacobian(q_derta_dd,[e_f1,e_f2,e_f3,e_f4,e_f5,e_f6,e_f7,e_f8,e_f9,e_f10]);
J0=[J1;J2];
if i==1
J0_0=J0;
else
J0_0=[J0_0;J0];
end
elseif isnan(data(i,2)) && isnan(data(i,3))
J0=zeros(2,10);
if i==1
J0_0=J0;
else
J0_0=[J0_0;J0];
end
else
p_derta_dd=p_q_matrix(2*i-1)-...
['e_f',num2str(2*i-1)]*((real(Y(i,1))*e_f1-imag(Y(i,1))*e_f2)+...
(real(Y(i,2))*e_f3-imag(Y(i,2))*e_f4)+...
(real(Y(i,3))*e_f5-imag(Y(i,3))*e_f6)+...
(real(Y(i,4))*e_f7-imag(Y(i,4))*e_f8)+...
(real(Y(i,5))*e_f9-imag(Y(i,5))*e_f10))-...
['e_f',num2str(2*i)]*((real(Y(i,1))*e_f2+imag(Y(i,1))*e_f1)+...
(real(Y(i,2))*e_f4+imag(Y(i,2))*e_f3)+...
(real(Y(i,3))*e_f6+imag(Y(i,3))*e_f5)+...
(real(Y(i,4))*e_f8+imag(Y(i,4))*e_f7)+...
(real(Y(i,5))*e_f10+imag(Y(i,5))*e_f9));
J1=jacobian(p_derta_dd,[e_f1,e_f2,e_f3,e_f4,e_f5,e_f6,e_f7,e_f8,e_f9,e_f10]);
v_derta_dd=p_v_matrix(2*i)^2-(eval(['e_f',num2str(2*i-1)])^2+eval(['e_f',num2str(2*i)])^2);
J2=jacobian(v_derta_dd,[e_f1,e_f2,e_f3,e_f4,e_f5,e_f6,e_f7,e_f8,e_f9,e_f10]);
J0=[J1;J2];
if i==1
J0_0=J0;
else
J0_0=[J0_0;J0];
end
end
end
for i=1:2*N
eval(['e_f',num2str(i),'=','e_f_matrix(i)',';']);
end
disp('雅可比矩阵为：')
J0_data=vpa(subs(J0_0),5);
disp(J0_data)
%% derta e f ——> e-f新 derta p q ——> p q新
e_f_derta=vpa(linsolve(J0_data,reshape(p_q_derta,10,1)),5);
e_f_matrix=vpa(reshape(e_f_matrix,10,1)-e_f_derta,5);
for w=1:N
if data(w,2)==0 && data(w,3)==0
[p_q_derta(2*w-1),p_q_derta(2*w)]=PQ(p_q_matrix(2*w-1),p_q_matrix(2*w),Y,e_f_matrix,w);
elseif isnan(data(w,2)) && isnan(data(w,3))
else
[p_q_derta(2*w-1),p_q_derta(2*w)]=PV(p_v_matrix(2*w-1),p_v_matrix(2*w),Y,e_f_matrix,w);
end
end
p_q_derta=vpa(p_q_matrix-p_q_derta,5);
p_q_derta=vpa(p_q_matrix-p_q_derta,5);
%% 输出新e-f阵和新p-q阵
disp('新的derta-ef阵')
disp(e_f_derta)
disp('新的ef阵')
disp(e_f_matrix)
disp('新的derta-pd阵')
disp(p_q_derta)
e_f_end=horzcat(e_f_end,e_f_matrix);
p_q_derta_end=horzcat(p_q_derta_end,p_q_derta);
end
%% 输出最后的所有矩阵变化情况
disp('最终节点功率误差的变化情况')
p_q_derta_end= reshape(p_q_derta_end,10,5)';
B = round(p_q_derta_end.* 1e5) ./ 1e5; % 将矩阵中的每个元素保留 5 位小数
C = reshape(B, [5 10]); % 将矩阵转换为一个 5 行 10 列的新矩阵
formatSpec = '%10.5f'; % 设置输出格式
for i = 1:5
fprintf([repmat(formatSpec, 1, 10) '\n'], C(i, :)); % 输出第 i 行元素，每行输出 10 个元素，用空格分隔，再换行
end
disp('最终节点电压的变化情况')
e_f_end= reshape(e_f_end,10,4)';
B0 = round(e_f_end.* 1e5) ./ 1e5; % 将矩阵中的每个元素保留 5 位小数
C0 = reshape(B0, [4 10]); % 将矩阵转换为一个 4 行 10 列的新矩阵
formatSpec = '%10.5f'; % 设置输出格式
for i = 1:4
fprintf([repmat(formatSpec, 1, 10) '\n'], C0(i, :)); % 输出第 i 行元素，每行输出 10 个元素，用空格分隔，再换行
end
for m=1:4
z = e_f_matrix(2*m-1) - e_f_matrix(2*m)*1i;
U0_zeita(2*m-1) = abs(z); % 计算幅值
U0_zeita(2*m) = -rad2deg(angle(z)); % 计算幅角（以度为单位）
end
disp("U_zeita矩阵")
disp(U0_zeita)
%% 函数 PQ节点
function [P,Q]=PQ(P0,Q0,Y,e_f,i)
P=P0-...
e_f(2*i-1)*((real(Y(i,1))*e_f(1)-imag(Y(i,1))*e_f(2))+...
(real(Y(i,2))*e_f(3)-imag(Y(i,2))*e_f(4))+...
(real(Y(i,3))*e_f(5)-imag(Y(i,3))*e_f(6))+...
(real(Y(i,4))*e_f(7)-imag(Y(i,4))*e_f(8))+...
(real(Y(i,5))*e_f(9)-imag(Y(i,5))*e_f(10)))-...
e_f(2*i)*((real(Y(i,1))*e_f(2)+imag(Y(i,1))*e_f(1))+...
(real(Y(i,2))*e_f(4)+imag(Y(i,2))*e_f(3))+...
(real(Y(i,3))*e_f(6)+imag(Y(i,3))*e_f(5))+...
(real(Y(i,4))*e_f(8)+imag(Y(i,4))*e_f(7))+...
(real(Y(i,5))*e_f(10)+imag(Y(i,5))*e_f(9)));
Q=Q0-...
e_f(2*i)*((real(Y(i,1))*e_f(1)-imag(Y(i,1))*e_f(2))+...
(real(Y(i,2))*e_f(3)-imag(Y(i,2))*e_f(4))+...
(real(Y(i,3))*e_f(5)-imag(Y(i,3))*e_f(6))+...
(real(Y(i,4))*e_f(7)-imag(Y(i,4))*e_f(8))+...
(real(Y(i,5))*e_f(9)-imag(Y(i,5))*e_f(10)))+...
e_f(2*i-1)*((real(Y(i,1))*e_f(2)+imag(Y(i,1))*e_f(1))+...
(real(Y(i,2))*e_f(4)+imag(Y(i,2))*e_f(3))+...
(real(Y(i,3))*e_f(6)+imag(Y(i,3))*e_f(5))+...
(real(Y(i,4))*e_f(8)+imag(Y(i,4))*e_f(7))+...
(real(Y(i,5))*e_f(10)+imag(Y(i,5))*e_f(9)));
end
%% PV节点
function [P,V]=PV(P0,V0,Y,e_f,i)
P=P0-...
e_f(2*i-1)*((real(Y(i,1))*e_f(1)-imag(Y(i,1))*e_f(2))+...
(real(Y(i,2))*e_f(3)-imag(Y(i,2))*e_f(4))+...
(real(Y(i,3))*e_f(5)-imag(Y(i,3))*e_f(6))+...
(real(Y(i,4))*e_f(7)-imag(Y(i,4))*e_f(8))+...
(real(Y(i,5))*e_f(9)-imag(Y(i,5))*e_f(10)))-...
e_f(2*i)*((real(Y(i,1))*e_f(2)+imag(Y(i,1))*e_f(1))+...
(real(Y(i,2))*e_f(4)+imag(Y(i,2))*e_f(3))+...
(real(Y(i,3))*e_f(6)+imag(Y(i,3))*e_f(5))+...
(real(Y(i,4))*e_f(8)+imag(Y(i,4))*e_f(7))+...
(real(Y(i,5))*e_f(10)+imag(Y(i,5))*e_f(9)));
V=V0^2-(e_f(2*i-1)^2+e_f(2*i)^2);
end
结果展示
>> chaoliujisuan
1节点是PQ节点。
2节点是PQ节点。
3节点是PQ节点。
4节点是PV节点。
5节点是V0节点。
-1.6000 -2.0000 -3.7000 5.0000 0
-0.5500 5.6980 2.0490 0 0
雅可比矩阵为：
[-1.3787, -6.5417, 0.62402, 3.9002, 0.75472, 2.6415, 0, 0, 0, 0]
[-6.0417, 1.3787, 3.9002, -0.62402, 2.6415, -0.75472, 0, 0, 0, 0]
[0.62402, 3.9002, -1.4539, -73.679, 0.82988, 3.112, 0, 63.492, 0, 0]
[ 3.9002, -0.62402, -60.283, 1.4539, 3.112, -0.82988, 63.492, 0, 0, 0]
[0.75472, 2.6415, 0.82988, 3.112, -1.5846, -39.087, 0, 0, 0, 31.746]
[ 2.6415, -0.75472, 3.112, -0.82988, -32.389, 1.5846, 0, 0, 31.746, 0]
[ 0, 0, 0, 66.667, 0, 0, 0, -63.492, 0, 0]
[ 0, 0, 0, 0, 0, 0, -2.1, 0, 0, 0]
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
警告: Solution is not unique because the system is rank-deficient.
> 位置：symengine
位置: sym/privBinaryOp (第 1214 行)
位置: sym/linsolve (第 63 行)
位置: chaoliujisuan (第 138 行)
新的derta-ef阵
0.033569
0.033482
-0.10538
-0.3607
-0.058813
0.069
0
-0.45749
0
0
新的ef阵
0.96643
-0.033482
1.1054
0.3607
1.0588
-0.069
1.05
0.45749
1.05
0
新的derta-pd阵
[-0.034734, 2.7753, 0.049045, -3.061, 0]
[-0.072036, 0.91802, -0.37145, -0.2093, 0]
雅可比矩阵为：
[0.048513, -6.8427, 0.73366, 3.7483, 0.81782, 2.5276, 0, 0, 0, 0]
[ -5.2259, 3.1347, 3.7483, -0.73366, 2.5276, -0.81782, 0, 0, 0, 0]
[-0.71702, 4.5363, 26.969, -74.858, -0.2052, 3.7393, -22.902, 70.183, 0, 0]
[ 4.5363, 0.71702, -74.27, -18.137, 3.7393, 0.2052, 70.183, 22.902, 0, 0]
[ 0.98137, 2.7448, 1.0934, 3.2378, -0.67478, -38.833, 0, 0, 2.1905, 33.613]
[ 2.7448, -0.98137, 3.2378, -1.0934, -36.627, 7.6126, 0, 0, 33.613, -2.1905]
[ 0, 0, -29.047, 66.667, 0, 0, 22.902, -70.183, 0, 0]
[ 0, 0, 0, 0, 0, 0, -2.1, -0.91498, 0, 0]
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
警告: Solution is not unique because the system is rank-deficient.
> 位置：symengine
位置: sym/privBinaryOp (第 1214 行)
位置: sym/linsolve (第 63 行)
位置: chaoliujisuan (第 138 行)
新的derta-ef阵
0.095285
0.036405
0.074975
0.030736
0.023675
0.0079805
0.071321
0.065054
0
0
新的ef阵
0.87115
-0.069887
1.0304
0.32997
1.0351
-0.07698
0.97868
0.39244
1.05
0
新的derta-pd阵
[-0.006758, 0.1666, 0.0032799, -0.1705, 0]
[-0.026557, 0.065414, -0.0094778, -0.0093187, 0]
雅可比矩阵为：
[ 0.10565, -6.4126, 0.81619, 3.354, 0.84207, 2.2484, 0, 0, 0, 0]
[ -4.3566, 3.3872, 3.354, -0.81619, 2.2484, -0.84207, 0, 0, 0, 0]
[-0.64393, 4.2247, 22.811, -69.824, -0.17176, 3.4805, -20.95, 65.423, 0, 0]
[ 4.2247, 0.64393, -69.17, -18.396, 3.4805, 0.17176, 65.423, 20.95, 0, 0]
[ 0.98458, 2.6762, 1.0986, 3.1575, -0.9257, -38.376, 0, 0, 2.4438, 32.862]
[ 2.6762, -0.98458, 3.1575, -1.0986, -35.367, 7.8571, 0, 0, 32.862, -2.4438]
[ 0, 0, -24.917, 62.138, 0, 0, 20.95, -65.423, 0, 0]
[ 0, 0, 0, 0, 0, 0, -1.9574, -0.78487, 0, 0]
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
警告: Solution is not unique because the system is rank-deficient.
> 位置：symengine
位置: sym/privBinaryOp (第 1214 行)
位置: sym/linsolve (第 63 行)
位置: chaoliujisuan (第 138 行)
新的derta-ef阵
0.011776
0.0018948
0.0043716
-0.00049323
0.0015996
0.00039477
0.0040511
0.00177
0
0
新的ef阵
0.85937
-0.071781
1.026
0.33046
1.0335
-0.077375
0.97463
0.39067
1.05
0
新的derta-pd阵
[-0.00020051, 0.00068996, 4.0216e-6, -0.00061814, 0]
[-0.00063863, -0.000019927, -0.000020817, -0.000019544, 0]
雅可比矩阵为：
[ 0.13507, -6.386, 0.81623, 3.3069, 0.83819, 2.2159, 0, 0, 0, 0]
[ -4.2298, 3.408, 3.3069, -0.81623, 2.2159, -0.83819, 0, 0, 0, 0]
[-0.64858, 4.2079, 22.694, -69.519, -0.17693, 3.4673, -20.982, 65.145, 0, 0]
[ 4.2079, 0.64858, -68.891, -18.592, 3.4673, 0.17693, 65.145, 20.982, 0, 0]
[ 0.98442, 2.6717, 1.0985, 3.1522, -0.93662, -38.331, 0, 0, 2.4564, 32.811]
[ 2.6717, -0.98442, 3.1522, -1.0985, -35.297, 7.8693, 0, 0, 32.811, -2.4564]
[ 0, 0, -24.804, 61.881, 0, 0, 20.982, -65.145, 0, 0]
[ 0, 0, 0, 0, 0, 0, -1.9493, -0.78133, 0, 0]
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
警告: Solution is not unique because the system is rank-deficient.
> 位置：symengine
位置: sym/privBinaryOp (第 1214 行)
位置: sym/linsolve (第 63 行)
位置: chaoliujisuan (第 138 行)
新的derta-ef阵
0.00021554
0.000039199
0.000027807
-0.000010828
0.00002032
7.5728e-6
0.000012922
-7.2225e-6
0
0
新的ef阵
0.85915
-0.071821
1.026
0.33047
1.0335
-0.077383
0.97462
0.39067
1.05
0
新的derta-pd阵
[-7.0314e-8, 1.3407e-8, -3.623e-10, 3.8685e-9, 0]
[-2.6639e-7, -1.1136e-8, -1.9611e-9, -2.1913e-10, 0]
最终节点功率误差的变化情况
-1.60000 -0.55000 -2.00000 5.69803 -3.70000 2.04902 5.00000 0.00000 0.00000 0.00000
-0.03473 -0.07204 2.77528 0.91802 0.04904 -0.37145 -3.06103 -0.20930 0.00000 0.00000
-0.00676 -0.02656 0.16660 0.06541 0.00328 -0.00948 -0.17050 -0.00932 0.00000 0.00000
-0.00020 -0.00064 0.00069 -0.00002 0.00000 -0.00002 -0.00062 -0.00002 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
最终节点电压的变化情况
0.96643 -0.03348 1.10538 0.36070 1.05881 -0.06900 1.05000 0.45749 1.05000 0.00000
0.87115 -0.06989 1.03041 0.32997 1.03514 -0.07698 0.97868 0.39244 1.05000 0.00000
0.85937 -0.07178 1.02604 0.33046 1.03354 -0.07738 0.97463 0.39067 1.05000 0.00000
0.85915 -0.07182 1.02601 0.33047 1.03352 -0.07738 0.97462 0.39067 1.05000 0.00000
U_zeita矩阵
0.8622 -4.7785 1.0779 17.8535 1.0364 -4.2819 1.0500 21.8433
>>
>> 
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
