Gauss-Seidel法求解电力系统节点电压时赋值错误求助

在电力系统分析项目中，使用Gauss-Seidel法求解节点电压时，执行以下代码行：

V(1,i) = (1/Y(i, i)) * (conj( S(1,i) / V(1,i) ) - Y(i, :) .* V(1,i));

出现错误：Unable to perform assignment because the size of the left side is 1-by-1 and the size of the right side is 1-by-6

完整原始代码如下：

%admittance matrix
Y = [-6.8 + 1i*42.033 , -6.748+1i*36.725 , 0 , -0.0532+1i*5.29 , 0 , 0 ;
    -6.748+1i*36.725, -6.75+1i*42.09   , -0.08+1i*5.35, 0 , 0 , 0 ;
    0 , -0.08+1i*5.35 , -0.0795+1i*5.3497 , 0 , 0 , 0 ;
    -0.0532+1i*5.29 , 0 , 0 , -12.3428+1i*134.8 , -12.29+1i*129.2 , 0 ;
    0 , 0 , 0 , -12.29+1i*129.2, -0.265+1i*12.19763 , -0.26+1i*13.08 ;
    0 , 0 , 0 , 0 , -0.2654+1i*13.08 , -0.265+1i*13.079 ] ;
S2 = 0;
S3 = 150e6*(0.95 + 1i*sin(acos(0.95)));
S4 = 0;
S5 = 0;
S6 = 300e6*(0.95 + 1i*sin(acos(0.95)));
S = [0, S2, S3, S4, S5, S6];
% B) Gauss-Seidel method to solve for nodal voltages
V = [1.05, 1.05, 1.05, 1.05, 1.05, 1.05]; % Initial guess
max_iter = 1000;
tolerance = 1e-6;

for iter = 1:max_iter
    
    for i = 2:6
        dVmax = 0;

        V_old = 1.05;

        V(1,i) = (1/Y(i, i)) * (conj( S(1,i) / V(1,i) ) - (Y(i, :) .* V(1,i)));
        dV = abs(V(i) - V_old );
        if (dV > dVmax)
            dVmax = dV;
            V_old = V(i);
        else
            V_old = V(i);
        end
    end
    if max(abs(V - V_old)) < tolerance
        break;
    end
display (V);
end

错误原因分析 #

1. 核心公式错误：Gauss-Seidel节点电压计算公式中，需要对第i行导纳矩阵与当前电压向量的乘积求和，再减去自身项Y(i,i)*V(i)，但原始代码中Y(i, :) .* V(1,i)是逐元素相乘得到1x6向量，导致右侧结果为向量，无法赋值给左侧的标量V(1,i)。
2. 迭代逻辑错误：V_old被固定初始化为1.05，未保存每次迭代前的实际电压值；dVmax在每次内层循环重置为0，无法跟踪全局最大电压变化；收敛判断条件使用了错误的变量对比逻辑。

修正后的代码 #

% 导纳矩阵
Y = [-6.8 + 1i*42.033 , -6.748+1i*36.725 , 0 , -0.0532+1i*5.29 , 0 , 0 ;
    -6.748+1i*36.725, -6.75+1i*42.09   , -0.08+1i*5.35, 0 , 0 , 0 ;
    0 , -0.08+1i*5.35 , -0.0795+1i*5.3497 , 0 , 0 , 0 ;
    -0.0532+1i*5.29 , 0 , 0 , -12.3428+1i*134.8 , -12.29+1i*129.2 , 0 ;
    0 , 0 , 0 , -12.29+1i*129.2, -0.265+1i*12.19763 , -0.26+1i*13.08 ;
    0 , 0 , 0 , 0 , -0.2654+1i*13.08 , -0.265+1i*13.079 ] ;
S2 = 0;
S3 = 150e6*(0.95 + 1i*sin(acos(0.95)));
S4 = 0;
S5 = 0;
S6 = 300e6*(0.95 + 1i*sin(acos(0.95)));
S = [0, S2, S3, S4, S5, S6];
% Gauss-Seidel法求解节点电压
V = [1.05, 1.05, 1.05, 1.05, 1.05, 1.05]; % 初始猜测值
max_iter = 1000;
tolerance = 1e-6;

for iter = 1:max_iter
    V_old = V; % 保存当前迭代前的电压向量
    dVmax = 0; % 初始化最大电压变化量
    
    for i = 2:6 % 遍历PQ节点（节点1为平衡节点，电压固定）
        % 修正核心公式：求和后减去自身项，确保结果为标量
        sum_YV = sum(Y(i, :) .* V) - Y(i,i)*V(i);
        V(i) = (1/Y(i,i)) * (conj(S(i)/V_old(i)) - sum_YV);
        
        % 更新全局最大电压变化量
        dV = abs(V(i) - V_old(i));
        if dV > dVmax
            dVmax = dV;
        end
    end
    
    % 收敛判断
    if dVmax < tolerance
        fprintf('迭代收敛，迭代次数：%d\n', iter);
        break;
    end
    
    if iter == max_iter
        fprintf('达到最大迭代次数，未收敛\n');
    end
end

disp('节点电压结果：');
disp(V);

关键修改点 #

• 修正核心计算公式：使用sum(Y(i, :) .* V) - Y(i,i)*V(i)计算除自身外的导纳-电压乘积和，确保结果为标量，匹配左侧赋值要求。
• 迭代前保存电压向量：每次外层循环开始时将当前V赋值给V_old，用于计算电压变化量和收敛判断。
• 全局跟踪最大变化：将dVmax初始化移到内层循环外，确保能记录整个迭代周期的最大电压变化。
• 优化收敛判断：用全局最大电压变化量dVmax与阈值对比，符合Gauss-Seidel法的收敛判定逻辑。

内容的提问来源于stack exchange，提问作者A Burial At Ornans