GaussSeidel.H File Reference

Declaration of the Gauss-Seidel load flow solver for power system analysis. More...

#include <Eigen/Dense>
#include <utility>
#include <vector>

Include dependency graph for GaussSeidel.H:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions
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.

Detailed Description

Declaration of the Gauss-Seidel load flow solver for power system analysis.

The Gauss-Seidel method is an iterative algorithm to solve the power flow equations for bus voltages in a power system. It uses the bus admittance matrix $$ Y_{bus} $$ and iteratively updates the bus voltages $$ V $$ to satisfy the power balance equations:

$$ S_i = V_i \sum_{j=1}^N Y_{ij}^* V_j^* $$

where:

• $$ S_i $$ is the specified complex power at bus $$ i $$
• $$ V_i $$ is the voltage at bus $$ i $$
• $$ Y_{ij} $$ is the element of the admittance matrix

The Gauss-Seidel voltage update for the k-th iteration is:

$$ V_i^{(k+1)} = \frac{1}{Y_{ii}} \left( \frac{S_i^*}{V_i^{*(k)}} - \sum_{j=1, j \neq i}^N Y_{ij} V_j^{(k+1 \text{ or } k)} \right) $$

For voltage-controlled buses (PV buses), the reactive power injection $$ Q_k $$ is unknown initially but can be computed at iteration i by:

$$ Q_k^{(i)} = V_k^{(i)} \sum_{n=1}^N |Y_{kn}| V_n^{(i)} \sin \left( \delta_k^{(i)} - \delta_n^{(i)} - \theta_{kn} \right) $$

The total reactive power generation at bus k is:

$$ Q_{Gk} = Q_k + Q_{Lk} $$

where $$( Q_{Lk} )$$ is the reactive load demand at the bus.

If the calculated $$( Q_{Gk} )$$ remains within its specified limits (e.g., $$ Q_{Gk}^{min} \leq Q_{Gk} \leq Q_{Gk}^{max} $$), the voltage magnitude $$ V_k $$ is fixed at its specified value, and only the voltage angle $$ \delta_k $$ is updated using the Gauss-Seidel formula. This means:

• The updated voltage magnitude $$ V_k^{(i+1)} $$ is set to the specified magnitude,
• Only the angle $$ \delta_k^{(i+1)} $$ changes according to the iterative update.

If $$( Q_{Gk} )$$ exceeds its limits during any iteration, the bus type is switched from voltage-controlled (PV) to load bus (PQ), with $$ Q_{Gk} $$ fixed at the violated limit (either $$ Q_{Gk}^{max} $$ or $$ Q_{Gk}^{min} $$). Under this condition, the voltage magnitude $$ V_k $$ is no longer fixed and is recalculated by the load flow program.

To accelerate convergence, a relaxation (or acceleration) coefficient $$ \omega $$ (typically $$(0 < \omega \leq 1) $$) can be applied to the voltage update:

$$ V_i^{(k+1)} = \omega \times V_i^{(k+1)} + (1 - \omega) \times V_i^{(k)} $$

The iteration proceeds until the maximum change in bus voltages between iterations is below a specified tolerance, or the maximum number of iterations is reached.

Definition in file GaussSeidel.H.

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: