% % 牛顿-拉夫逊法潮流计算(标幺值) - 专家修正版
clear; clc; close all;
% 1. 基础参数设置
n = 10; % 节点数
tol = 1e-6; % 收敛判据
max_iter = 50; % 最大迭代次数
S_base = 100; % 基准容量(MVA)
% 2. 节点数据(精细调整)
bus_data = [
1, 3, 6.0, 1.0, 1.01641, 0.0; % G1: PQ节点
2, 2, 6.0, 0.0, 1.00000, 0.0; % G2: PV节点
3, 3, 6.0, 1.0, 1.01641, 0.0; % G3: PQ节点
4, 3, 0.0, 0.0, 1.05872, 0.0; % B1-500: PQ节点
5, 3, -3.0, -2.0, 1.04097, 0.0; % B2-220: 负荷PQ节点
6, 3, 0.0, 0.0, 0.95738, 0.0; % B3-500: PQ节点
7, 3, 0.0, 0.0, 1.02020, 0.0; % B4-500: PQ节点
8, 3, 0.0, 0.0, 1.00654, 0.0; % B5-500: PQ节点
9, 3, -3.0, -2.0, 1.05709, 0.0; % B6-220: 负荷PQ节点
10,1, 0.0, 0.0, 1.00000, 0.0; % S1: Slack节点
];
% 3. 导纳矩阵计算(精确参数)
branch_params = [
4, 7, 0.00377, 0.0406; % AC1
4, 7, 0.00377, 0.0406; % AC2
8, 6, 0.00377, 0.0406; % AC3
8, 6, 0.00377, 0.0406; % AC4
7, 8, 0.00377, 0.0406; % AC5
7, 8, 0.00377, 0.0406; % AC6
1, 4, 0.0, 0.01; % T2W9(变压器)
2, 5, 0.0, 0.00667; % T2W10(变压器)
4, 5, 0.0, 0.02; % T2W11(变压器)
10,6, 0.0, 0.00467; % T2W12(变压器)
3, 4, 0.0, 0.01; % T2W13(变压器)
7, 9, 0.0, 0.00667; % T2W14(变压器)
];
% 专业导纳矩阵计算
Y = zeros(n, n, 'like', 1j); % 初始化复数矩阵
for k = 1:size(branch_params,1)
i = branch_params(k,1);
j = branch_params(k,2);
R = branch_params(k,3);
X = branch_params(k,4);
% 变压器支路特殊处理
if abs(R) < 1e-10
% 变压器导纳模型
Y_tr = 1/(1j*X);
Y(i,i) = Y(i,i) + Y_tr;
Y(j,j) = Y(j,j) + Y_tr;
Y(i,j) = Y(i,j) - Y_tr;
Y(j,i) = Y(j,i) - Y_tr;
else
% 线路导纳模型
Z = R + 1j*X;
Y_line = 1/Z;
Y(i,i) = Y(i,i) + Y_line;
Y(j,j) = Y(j,j) + Y_line;
Y(i,j) = Y(i,j) - Y_line;
Y(j,i) = Y(j,i) - Y_line;
end
end
% 4. 使用真实初始电压值
U = complex(bus_data(:,5) .* cos(bus_data(:,6)), ...
bus_data(:,5) .* sin(bus_data(:,6)));
% 5. 牛顿-拉夫逊迭代(专业实现)
fprintf('潮流计算迭代过程:\n');
fprintf('迭代次数\t最大功率误差\n');
converged = false;
for iter = 1:max_iter
% 计算节点注入功率
I_inj = Y * U;
S_calc = U .* conj(I_inj);
P_calc = real(S_calc);
Q_calc = imag(S_calc); % 保存Q_calc用于后续计算
% 计算功率误差(专业处理节点类型)
P_error = bus_data(:,3) - P_calc;
Q_error = bus_data(:,4) - Q_calc;
% 节点类型约束
slack_nodes = (bus_data(:,2) == 1);
PV_nodes = (bus_data(:,2) == 2);
P_error(slack_nodes) = 0;
Q_error(slack_nodes | PV_nodes) = 0;
% 检查收敛
max_error = max([abs(P_error); abs(Q_error)]);
fprintf('%d\t\t%.6e\n', iter, max_error);
if max_error <= tol
converged = true;
break;
end
% 构建雅可比矩阵(修正Q_calc问题)
[J, p_idx, q_idx] = build_jacobi_matrix(Y, U, bus_data, Q_calc); % 传递Q_calc
% 构建误差向量
delta_S = [P_error(p_idx); Q_error(q_idx)];
% 求解修正方程(增强稳定性)
if rcond(J) < 1e-10
damping = 1e-4;
J_reg = J' * J + damping * eye(size(J));
delta_X = J_reg \ (J' * (-delta_S));
else
delta_X = J \ (-delta_S);
end
% 更新节点电压
U = update_node_voltage(U, delta_X, bus_data, p_idx, q_idx);
end
% 6. 输出结果(专业格式)
if converged
fprintf('\n潮流计算收敛于%d次迭代!\n', iter);
% 节点电压结果
fprintf('\n节点电压结果(标幺值):\n');
fprintf('节点\t名称\t\t电压幅值\t电压相角(°)\n');
node_names = {'G1','G2','G3','B1-500','B2-220','B3-500','B4-500','B5-500','B6-220','S1'};
for i = 1:n
U_mag = abs(U(i));
U_angle = angle(U(i)) * 180/pi;
fprintf('%d\t%s\t\t%.5f\t\t%.2f\n', i, node_names{i}, U_mag, U_angle);
end
% 支路功率结果
fprintf('\n支路功率结果(标幺值,始端→末端):\n');
branch_names = {'AC1','AC2','AC3','AC4','AC5','AC6','T2W9','T2W10','T2W11','T2W12','T2W13','T2W14'};
for k = 1:size(branch_params,1)
i = branch_params(k,1);
j = branch_params(k,2);
% 精确计算支路功率
if abs(branch_params(k,3)) < 1e-10 % 变压器
V_i = U(i);
V_j = U(j);
I_ij = (V_i - V_j) / (1j * branch_params(k,4));
S_ij = V_i * conj(I_ij);
else % 线路
V_i = U(i);
V_j = U(j);
I_ij = (V_i - V_j) / complex(branch_params(k,3), branch_params(k,4));
S_ij = V_i * conj(I_ij);
end
fprintf('%s\t%s→%s\tP=%.3f\tQ=%.3f\n', ...
branch_names{k}, node_names{i}, node_names{j}, real(S_ij), imag(S_ij));
end
% 系统功率平衡检查
fprintf('\n系统功率平衡检查:\n');
total_P_gen = sum(bus_data(bus_data(:,3) > 0, 3));
total_P_load = -sum(bus_data(bus_data(:,3) < 0, 3));
total_P_loss = total_P_gen - total_P_load;
fprintf('总发电功率: %.3f pu (%.1f MW)\n', total_P_gen, total_P_gen * S_base);
fprintf('总负荷功率: %.3f pu (%.1f MW)\n', total_P_load, total_P_load * S_base);
fprintf('总网损功率: %.3f pu (%.1f MW)\n', total_P_loss, total_P_loss * S_base);
else
fprintf('\n警告: 潮流计算在%d次迭代后未收敛!\n', max_iter);
end
% ======== 专业辅助函数 ========
% 修正的雅可比矩阵函数(添加Q_calc参数)
function [J, p_idx, q_idx] = build_jacobi_matrix(Y, U, bus_data, Q_calc)
n_nodes = length(U);
G = real(Y);
B = imag(Y);
V = abs(U);
theta = angle(U);
% 确定参与方程的节点
p_idx = find(bus_data(:,2) ~= 1); % P方程节点(非松弛)
q_idx = find(bus_data(:,2) == 3); % Q方程节点(PQ节点)
np = length(p_idx);
nq = length(q_idx);
J = zeros(np + nq); % 正确预分配
% 1. H子矩阵: ∂P/∂θ
for i = 1:np
m = p_idx(i);
for j = 1:np
n = p_idx(j);
if m == n
% 对角元素 - ∑V_mV_k(G_mk sinθ_mk - B_mk cosθ_mk)
H_ii = -V(m)^2 * B(m,m);
for k = 1:n_nodes
if k ~= m
angle_diff = theta(m) - theta(k);
H_ii = H_ii - V(m)*V(k)*(G(m,k)*sin(angle_diff) - B(m,k)*cos(angle_diff));
end
end
J(i,j) = H_ii;
else
% 非对角元素 V_mV_n(G_mn sinθ_mn - B_mn cosθ_mn)
angle_diff = theta(m) - theta(n);
J(i,j) = V(m)*V(n)*(G(m,n)*sin(angle_diff) - B(m,n)*cos(angle_diff));
end
end
end
% 2. N子矩阵: ∂P/∂V
for i = 1:np
m = p_idx(i);
for j = 1:nq
n = q_idx(j);
if m == n
% 对角元素 ∑V_k(G_mk cosθ_mk + B_mk sinθ_mk)
N_ii = V(m)*G(m,m);
for k = 1:n_nodes
if k ~= m
angle_diff = theta(m) - theta(k);
N_ii = N_ii + V(k)*(G(m,k)*cos(angle_diff) + B(m,k)*sin(angle_diff));
end
end
J(i, np+j) = N_ii;
else
% 非对角元素 V_m(G_mn cosθ_mn + B_mn sinθ_mn)
angle_diff = theta(m) - theta(n);
J(i, np+j) = V(m)*(G(m,n)*cos(angle_diff) + B(m,n)*sin(angle_diff));
end
end
end
% 3. M子矩阵: ∂Q/∂θ(关键修正)
for i = 1:nq
m = q_idx(i);
for j = 1:np
n = p_idx(j);
if m == n
% 对角元素(标准公式)V_m^2 G_mm - Q_calc(m)
% 使用传入的Q_calc值
J(np+i, j) = V(m)^2 * G(m,m) - Q_calc(m);
else
% 非对角元素 -V_mV_n(G_mn cosθ_mn + B_mn sinθ_mn)
angle_diff = theta(m) - theta(n);
J(np+i, j) = -V(m)*V(n)*(G(m,n)*cos(angle_diff) + B(m,n)*sin(angle_diff));
end
end
end
% 4. L子矩阵: ∂Q/∂V
for i = 1:nq
m = q_idx(i);
for j = 1:nq
n = q_idx(j);
if m == n
% 对角元素 ∑V_k(G_mk sinθ_mk - B_mk cosθ_mk) - V_m B_mm
L_ii = -V(m)*B(m,m);
for k = 1:n_nodes
if k ~= m
angle_diff = theta(m) - theta(k);
L_ii = L_ii + V(k)*(G(m,k)*sin(angle_diff) - B(m,k)*cos(angle_diff));
end
end
J(np+i, np+j) = L_ii;
else
% 非对角元素 V_m(G_mn sinθ_mn - B_mn cosθ_mn)
angle_diff = theta(m) - theta(n);
J(np+i, np+j) = V(m)*(G(m,n)*sin(angle_diff) - B(m,n)*cos(angle_diff));
end
end
end
end % 添加结束语句
% 电压更新函数
function U_new = update_node_voltage(U_old, delta_X, bus_data, p_idx, q_idx)
np = length(p_idx);
nq = length(q_idx);
% 复制当前电压状态
V = abs(U_old);
theta = angle(U_old);
% 1. 更新相角(δ)
for i = 1:np
node_idx = p_idx(i);
theta(node_idx) = theta(node_idx) + delta_X(i);
end
% 2. 更新电压幅值(V)
for i = 1:nq
node_idx = q_idx(i);
delta_V = delta_X(np + i);
V(node_idx) = V(node_idx) + delta_V;
% 电压幅值约束(0.95-1.05 pu范围)
if V(node_idx) < 0.95
V(node_idx) = 0.95;
elseif V(node_idx) > 1.05
V(node_idx) = 1.05;
end
end
% 3. 重构复数电压
U_new = V .* exp(1j * theta);
% 4. 保持PV节点电压幅值不变
PV_nodes = bus_data(:,2) == 2;
U_new(PV_nodes) = abs(U_new(PV_nodes)) .* (bus_data(PV_nodes,5) ./ abs(bus_data(PV_nodes,5))) .* exp(1j * angle(U_new(PV_nodes)));
% 5. 保持Slack节点恒定
slack_nodes = bus_data(:,2) == 1;
U_new(slack_nodes) = bus_data(slack_nodes,5) .* exp(1j * bus_data(slack_nodes,6));
end % 添加结束语句 潮流计算迭代过程:
迭代次数 最大功率误差
1 1.189029e+01
2 2.531832e+01
3 2.677794e+01
4 3.255412e+01
5 7.154282e+01
6 1.419108e+02
7 2.287291e+02
8 3.561172e+02
9 3.217857e+02
10 3.941556e+02
11 4.699058e+02
12 3.120394e+02
13 4.073019e+02
14 5.331827e+02
15 5.947108e+02
16 6.026744e+02
17 5.795085e+02
18 5.285860e+02
19 3.280251e+02
20 4.931953e+02
21 5.653989e+02
22 5.311069e+02
23 5.247585e+02
24 5.123802e+02
25 5.058910e+02
26 4.839999e+02
27 5.928759e+02
28 6.378344e+02
29 6.367423e+02
30 6.388663e+02
31 6.207253e+02
32 5.416889e+02
33 5.519267e+02
34 3.240033e+02
35 3.729455e+02
36 3.632753e+02
37 3.677791e+02
38 3.230372e+02
39 4.078947e+02
40 4.273080e+02
41 3.245169e+02
42 3.908748e+02
43 4.134421e+02
44 2.704812e+02
45 3.173655e+02
46 3.557677e+02
47 4.275630e+02
48 4.231329e+02
49 4.345852e+02
50 4.071010e+02
警告: 潮流计算在50次迭代后未收敛!
>>
最新发布