NewtonRaphson.H File Reference

Declaration of the Newton-Raphson load flow solver for power system analysis. More...

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

Include dependency graph for NewtonRaphson.H:

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

Go to the source code of this file.

Functions
bool	NewtonRaphson (const Eigen::MatrixXd &G, const Eigen::MatrixXd &B, const Eigen::VectorXd &Ps, const Eigen::VectorXd &Qs, Eigen::VectorXd &V, Eigen::VectorXd &delta, int n_bus, int n_pq, const std::vector< int > &pq_bus_id, int maxIter=1024, double tolerance=1E-8, std::vector< std::pair< int, double > > *iterHistory=nullptr)
	Solves the power flow equations using the Newton-Raphson iterative method.

Detailed Description

Declaration of the Newton-Raphson load flow solver for power system analysis.

The Newton-Raphson method is a robust and widely-used iterative algorithm for solving the nonlinear power flow equations in electrical power systems.

It uses the bus admittance matrix $$Y_{bus}$$ and iteratively updates the bus voltages $$V$$ by solving the set of nonlinear equations:

$$\Delta x^{(k)} = -[J(x^{(k)})]^{-1} F(x^{(k)})$$

where:

• $$x$$ is the state vector (bus voltage magnitudes and angles),
• $$F(x)$$ is the vector of power mismatch equations,
• $$J(x)$$ is the Jacobian matrix of partial derivatives with respect to $$x$$,
• $$\Delta x$$ is the update to the state vector at iteration $$k$$.

The power injection equations at each bus are:

Active power (P):

$$P_i = \sum_{j=1}^{n} |V_i||V_j|(G_{ij}\cos(\theta_{ij}) + B_{ij}\sin(\theta_{ij}))$$

Reactive power (Q):

$$Q_i = \sum_{j=1}^{n} |V_i||V_j|(G_{ij}\sin(\theta_{ij}) - B_{ij}\cos(\theta_{ij}))$$

where:

• $$G_{ij} + jB_{ij}$$ is the element of the admittance matrix $$Y_{bus}$$,
• $$\theta_{ij} = \delta_i - \delta_j$$ is the voltage angle difference between buses $$i$$ and $$j$$,
• $$V_i, V_j$$ are the voltage magnitudes at buses $$i$$ and $$j$$.

The Jacobian matrix has the general block structure:

$$J = \begin{bmatrix} \frac{\partial P}{\partial \delta} & \frac{\partial P}{\partial |V|} \\ \frac{\partial Q}{\partial \delta} & \frac{\partial Q}{\partial |V|} \end{bmatrix}$$

with submatrices defined as follows:

$$J_1(kn) = \frac{\partial P_k}{\partial \delta_n} = \begin{cases} -V_k \sum_{n=1}^N V_n |Y_{kn}| \sin(\delta_k - \delta_n - \theta_{kn}) & \text{if } k = n \ V_k V_n |Y_{kn}| \sin(\delta_k - \delta_n - \theta_{kn}) & \text{if } k \ne n \end{cases}$$

$$J_2(kn) = \frac{\partial P_k}{\partial |V_n|} = \begin{cases} V_k |Y_{kk}| \cos(\theta_{kk}) + \sum_{n \ne k} V_n |Y_{kn}| \cos(\delta_k - \delta_n - \theta_{kn}) & \text{if } k = n \ V_k |Y_{kn}| \cos(\delta_k - \delta_n - \theta_{kn}) & \text{if } k \ne n \end{cases}$$

$$J_3(kn) = \frac{\partial Q_k}{\partial \delta_n} = \begin{cases} V_k \sum_{n=1}^N V_n |Y_{kn}| \cos(\delta_k - \delta_n - \theta_{kn}) & \text{if } k = n \ -V_k V_n |Y_{kn}| \cos(\delta_k - \delta_n - \theta_{kn}) & \text{if } k \ne n \end{cases}$$

$$J_4(kn) = \frac{\partial Q_k}{\partial |V_n|} = \begin{cases} -V_k |Y_{kk}| \sin(\theta_{kk}) + \sum_{n \ne k} V_n |Y_{kn}| \sin(\delta_k - \delta_n - \theta_{kn}) & \text{if } k = n \ V_k |Y_{kn}| \sin(\delta_k - \delta_n - \theta_{kn}) & \text{if } k \ne n \end{cases}$$

The process continues until all mismatches are within a specified tolerance or the maximum number of iterations is reached.

Definition in file NewtonRaphson.H.

Function Documentation

NewtonRaphson()

bool NewtonRaphson	(	const Eigen::MatrixXd &	G,
		const Eigen::MatrixXd &	B,
		const Eigen::VectorXd &	Ps,
		const Eigen::VectorXd &	Qs,
		Eigen::VectorXd &	V,
		Eigen::VectorXd &	delta,
		int	n_bus,
		int	n_pq,
		const std::vector< int > &	pq_bus_id,
		int	maxIter = 1024,
		double	tolerance = 1E-8,
		std::vector< std::pair< int, double > > *	iterHistory = nullptr
	)

Solves the power flow equations using the Newton-Raphson iterative method.

This is a pure solver function that takes G, B matrices, scheduled power, initial voltage/angle guesses, and bus classification. It modifies V and delta in-place.

Parameters

G	Conductance matrix (real part of $$ Y_{bus} $$).
B	Susceptance matrix (imaginary part of $$ Y_{bus} $$).
Ps	Scheduled active power injections (Pg - Pl) [p.u.].
Qs	Scheduled reactive power injections (Qg - Ql) [p.u.].
V	(in/out) Voltage magnitudes [p.u.].
delta	(in/out) Voltage angles [rad].
n_bus	Total number of buses.
n_pq	Number of PQ buses.
pq_bus_id	0-based indices of PQ buses.
maxIter	Maximum number of iterations (default: 1024).
tolerance	Convergence tolerance for power mismatches (default: $$ 1 \times 10^{-8} $$).
iterHistory	Optional pointer to store iteration number and error at each step.

Returns

true if converged, false otherwise.

Definition at line 39 of file NewtonRaphson.C.

  {
    // Compute initial mismatch
    Eigen::VectorXd P(n_bus), Q(n_bus);
    Eigen::VectorXd mismatch = powerMismatch(Ps, Qs, G, B, V, delta, n_bus, pq_bus_id, P, Q);
 
    double error = mismatch.cwiseAbs().maxCoeff();
    int iter = 0;
 
    if (iterHistory) {
        iterHistory->clear();
        iterHistory->emplace_back(0, error);
    }
 
    while (error >= tolerance) {
        if (iter >= maxIter) {
            printConvergenceStatus("Newton-Raphson", false, iter, maxIter, error, tolerance);
            LOG_WARN("Newton-Raphson did not converge within {} iterations.", maxIter);
            LOG_DEBUG("Final max mismatch was {:.6e}, tolerance is {:.6e}.", error, tolerance);
            return false;
        }
        iter++;
 
        // Build Jacobian
        Eigen::MatrixXd J = computeJacobian(V, delta, n_bus, n_pq, pq_bus_id, G, B, P, Q);
 
        // Solve J * correction = mismatch
        Eigen::VectorXd correction = J.colPivHouseholderQr().solve(mismatch);
 
        // Update delta for non-slack buses (indices 1..N-1)
        for (int i = 1; i < n_bus; ++i) {
            delta(i) += correction(i - 1);
        }
 
        // Update V for PQ buses only
        for (int k = 0; k < n_pq; ++k) {
            V(pq_bus_id[k]) += correction(n_bus - 1 + k);
        }
 
        // Recompute mismatch
        mismatch = powerMismatch(Ps, Qs, G, B, V, delta, n_bus, pq_bus_id, P, Q);
 
        error = mismatch.cwiseAbs().maxCoeff();
        if (iterHistory) iterHistory->emplace_back(iter, error);
        printIterationProgress("Newton-Raphson", iter, maxIter, error, tolerance);
        LOG_DEBUG("NR iteration {}: max mismatch = {:.16e}", iter, error);
    }
 
    printConvergenceStatus("Newton-Raphson", true, iter, maxIter, error, tolerance);
    LOG_DEBUG("Newton-Raphson converged in {} iterations with max mismatch {:.6e}", iter, error);
    return true;
}

References computeJacobian(), LOG_DEBUG, LOG_WARN, powerMismatch(), printConvergenceStatus(), and printIterationProgress().

Here is the call graph for this function: