GaussSeidel.C File Reference

Gauss-Seidel power flow solver implementation. More...

#include <cmath>
#include <complex>
#include <iostream>
#include <limits>
#include "GaussSeidel.H"
#include "Logger.H"
#include "Progress.H"

Include dependency graph for GaussSeidel.C:

Go to the source code of this file.

Functions
bool	GaussSeidel (const Eigen::MatrixXcd &Y, Eigen::VectorXd &Vmag, Eigen::VectorXd &delta, const Eigen::VectorXi &type_bus, const Eigen::VectorXd &P, const Eigen::VectorXd &Q, int N, int maxIter, double tolerance, double omega, std::vector< std::pair< int, double > > *iterHistory)
	Solves the power flow equations using the Gauss-Seidel iterative method.

Detailed Description

Gauss-Seidel power flow solver implementation.

Definition in file GaussSeidel.C.

Function Documentation

GaussSeidel()

bool GaussSeidel	(	const Eigen::MatrixXcd &	Y,
		Eigen::VectorXd &	V,
		Eigen::VectorXd &	delta,
		const Eigen::VectorXi &	type_bus,
		const Eigen::VectorXd &	P,
		const Eigen::VectorXd &	Q,
		int	N,
		int	maxIter = 1024,
		double	tolerance = 1E-8,
		double	omega = 1.0,
		std::vector< std::pair< int, double > > *	iterHistory = nullptr
	)

Solves the power flow equations using the Gauss-Seidel iterative method.

Pure solver: no Q-limit checking (handled by outer loop via checkQlimits). PV buses maintain their scheduled voltage magnitude; PQ buses are fully updated.

Parameters

Y	The bus admittance matrix ($$ Y_{bus} $$).
V	(in/out) Voltage magnitudes [p.u.]. PV buses maintain scheduled value.
delta	(in/out) Voltage angles [rad].
type_bus	Bus type vector (1=Slack, 2=PV, 3=PQ).
P	Scheduled net active power injections [p.u.] ($$ P_g - P_l $$).
Q	Scheduled net reactive power injections [p.u.] ($$ Q_g - Q_l $$). Used only for PQ buses.
N	Total number of buses.
maxIter	Maximum number of iterations (default: 1024).
tolerance	Convergence tolerance for bus voltage updates (default: $$ 1 \times 10^{-8} $$).
omega	Relaxation parameter for the Successive Over-Relaxation (SOR) method (default: 1.0; SOR not applied).
iterHistory	Optional pointer to store iteration number and error at each step.

Returns

true if the algorithm converged within the specified number of iterations, false otherwise.

Definition at line 35 of file GaussSeidel.C.

  {
    // Store scheduled voltage magnitudes for PV buses
    Eigen::VectorXd Vmag_sched = Vmag;
 
    // Initialize complex voltage vector
    Eigen::VectorXcd V(N);
    for (int i = 0; i < N; ++i)
        V(i) = std::polar(Vmag(i), delta(i));
 
    int iteration = 0;
    double error = std::numeric_limits<double>::infinity();
 
    if (omega <= 0.0 || omega >= 2.0) {
        LOG_WARN("Invalid input: Relaxation coefficient must be between 0 and 2.");
        LOG_DEBUG("Setting Relaxation coefficient to 1.");
        omega = 1.0;
    }
 
    if (omega < 1.0) {
        LOG_CRITICAL("Under-relaxation enabled (omega < 1), this will slow down convergence.");
    }
    else if (omega == 1.0) {
        LOG_DEBUG("Standard Gauss-Seidel enabled (omega = 1).");
    }
    else if (omega > 1.0) {
        LOG_DEBUG("Over-relaxation enabled (omega > 1), this will accelerate convergence.");
    }
 
    LOG_DEBUG("Relaxation Coefficient :: {}", omega);
 
    while (error >= tolerance && iteration < maxIter) {
        Eigen::VectorXcd dV = Eigen::VectorXcd::Zero(N);
 
        for (int n = 0; n < N; ++n) {
            if (type_bus(n) == 1) continue;  // Skip slack bus
 
            std::complex<double> In = Y.row(n) * V;
 
            if (type_bus(n) == 2) {  // PV Bus
                // Compute Q from current solution (no clamping)
                double Qn = -std::imag(std::conj(V(n)) * In);
 
                // GS update maintaining scheduled voltage magnitude
                std::complex<double> I_excl = In - Y(n, n) * V(n);
                std::complex<double> V_updated = ((P(n) - std::complex<double>(0, 1) * Qn) / std::conj(V(n)) - I_excl) / Y(n, n);
                std::complex<double> V_corrected = Vmag_sched(n) * V_updated / std::abs(V_updated);
                dV(n) = V_corrected - V(n);
                V(n) = V_corrected;
 
            } else if (type_bus(n) == 3) {  // PQ Bus
                std::complex<double> I_excl = In - Y(n, n) * V(n);
                std::complex<double> V_updated = ((P(n) - std::complex<double>(0, 1) * Q(n)) / std::conj(V(n)) - I_excl) / Y(n, n);
                std::complex<double> V_relaxed = V(n) + omega * (V_updated - V(n));
                dV(n) = V_relaxed - V(n);
                V(n) = V_relaxed;
            }
        }
 
        error = dV.norm();
        iteration++;
        if (iterHistory) iterHistory->emplace_back(iteration, error);
        printIterationProgress("Gauss-Seidel", iteration, maxIter, error, tolerance);
    }
 
    if (iteration >= maxIter) {
        printConvergenceStatus("Gauss-Seidel", false, iteration, maxIter, error, tolerance);
        LOG_WARN("Gauss-Seidel did not converge within max iterations ({}).", maxIter);
        LOG_DEBUG("Final error norm was {:.6e}, tolerance is {:.6e}.", error, tolerance);
        return false;
    }
 
    // Extract converged V magnitudes and angles
    for (int i = 0; i < N; ++i) {
        Vmag(i) = std::abs(V(i));
        delta(i) = std::arg(V(i));
    }
 
    printConvergenceStatus("Gauss-Seidel", true, iteration, maxIter, error, tolerance);
    LOG_DEBUG("Gauss-Seidel converged in {} iterations with error norm {:.6e}.", iteration, error);
 
    return true;
}

References LOG_CRITICAL, LOG_DEBUG, LOG_WARN, printConvergenceStatus(), and printIterationProgress().

Here is the call graph for this function: