eneros-powerflow

eneros-powerflow provides a power system Load Flow solver, using Newton-Raphson as the primary method, with support for PQ decoupling and DC Power Flow as fast approximations. Based on eneros-linalg’s sparse LU decomposition, the typical solve time for a 1000-node network is < 10ms. Load Flow results serve as input for upper-layer analyses such as constraint validation, OPF, and Twin models.

Responsibilities

eneros-powerflow handles the following core responsibilities:

• YBus Construction: Builds sparse CSR-format nodal admittance matrix based on the network topology provided by eneros-topology
• Newton-Raphson Solving: Uses full Jacobian matrix iterative solving, automatically handling PQ / PV / Slack bus types
• Fast Approximation: Provides two lightweight solvers: Fast Decoupled and DC Power Flow
• Convergence Diagnostics: Power residual < 1e-8 convergence criterion; returns residual curve for analysis on divergence
• Result Caching: Reuses previous results when topology is unchanged, avoiding redundant computation
• Branch Flow: Calculates active power, reactive power, losses, and loading percentage for each branch after solving

Key Types and Interfaces

PowerFlowSolver

use std::time::Duration;
use eneros_core::{BusId, BranchId, Result};
use eneros_topology::NetworkGraph;

pub struct PowerFlowSolver {
    method: PowerFlowMethod,
    tolerance: f64,
    max_iterations: usize,
    enable_cache: bool,
}

#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum PowerFlowMethod {
    NewtonRaphson,
    FastDecoupled,
    DcPowerFlow,
}

impl Default for PowerFlowSolver {
    fn default() -> Self {
        Self {
            method: PowerFlowMethod::NewtonRaphson,
            tolerance: 1e-8,
            max_iterations: 50,
            enable_cache: true,
        }
    }
}

impl PowerFlowSolver {
    pub fn new() -> Self {
        Self::default()
    }

    pub fn method(mut self, m: PowerFlowMethod) -> Self {
        self.method = m;
        self
    }

    pub fn tolerance(mut self, t: f64) -> Self {
        self.tolerance = t;
        self
    }

    pub fn max_iterations(mut self, n: usize) -> Self {
        self.max_iterations = n;
        self
    }

    pub fn enable_cache(mut self, on: bool) -> Self {
        self.enable_cache = on;
        self
    }

    pub fn solve(&self, network: &NetworkGraph) -> Result<PowerFlowResult> {
        let ybus = YBusMatrix::build(network);
        let start = std::time::Instant::now();
        let inner = match self.method {
            PowerFlowMethod::NewtonRaphson => Self::solve_nr(network, &ybus, self.tolerance, self.max_iterations)?,
            PowerFlowMethod::FastDecoupled => Self::solve_fd(network, &ybus, self.tolerance, self.max_iterations)?,
            PowerFlowMethod::DcPowerFlow => Self::solve_dc(network, &ybus)?,
        };
        let duration = start.elapsed();
        let branch_flows = Self::compute_branch_flows(network, &inner.buses);
        Ok(PowerFlowResult {
            converged: inner.converged,
            iterations: inner.iterations,
            residual: inner.residual,
            duration,
            buses: inner.buses,
            branch_flows,
        })
    }
}

YBusMatrix

use eneros_linalg::SparseMatrix;

pub struct YBusMatrix {
    pub data: SparseMatrix<Complex64>,
    pub bus_order: Vec<BusId>,
}

impl YBusMatrix {
    pub fn build(network: &NetworkGraph) -> Self {
        let buses: Vec<BusId> = collect_bus_ids(network);
        let n = buses.len();
        let mut matrix = SparseMatrix::new(n, n);
        let index: HashMap<BusId, usize> = buses
            .iter()
            .enumerate()
            .map(|(i, id)| (id.clone(), i))
            .collect();

        for branch in network.iter_branches() {
            if branch.status != BranchStatus::Closed {
                continue;
            }
            let i = index[&branch.from];
            let j = index[&branch.to];
            let z = Complex64::new(branch.params.r_pu, branch.params.x_pu);
            let y = Complex64::new(0.0, 0.0) - z.inv();
            matrix.add(i, j, -y);
            matrix.add(j, i, -y);
            matrix.add(i, i, y + Complex64::new(0.0, branch.params.b_pu / 2.0));
            matrix.add(j, j, y + Complex64::new(0.0, branch.params.b_pu / 2.0));
        }

        Self { data: matrix, bus_order: buses }
    }
}

PowerFlowResult

use eneros_core::{Voltage, BusId, BranchId};

pub struct PowerFlowResult {
    pub converged: bool,
    pub iterations: usize,
    pub residual: f64,
    pub duration: Duration,
    pub buses: Vec<BusSolution>,
    pub branch_flows: Vec<BranchFlow>,
}

#[derive(Debug, Clone)]
pub struct BusSolution {
    pub id: BusId,
    pub voltage: Voltage,
    pub p_pu: f64,
    pub q_pu: f64,
}

#[derive(Debug, Clone)]
pub struct BranchFlow {
    pub from: BusId,
    pub to: BusId,
    pub p_from_pu: f64,
    pub q_from_pu: f64,
    pub p_to_pu: f64,
    pub q_to_pu: f64,
    pub loading_percent: f64,
    pub loss_pu: f64,
}

impl PowerFlowResult {
    pub fn bus_by_id(&self, id: &BusId) -> Option<&BusSolution> {
        self.buses.iter().find(|b| &b.id == id)
    }

    pub fn total_loss_pu(&self) -> f64 {
        self.branch_flows.iter().map(|f| f.loss_pu).sum()
    }

    pub fn max_loading(&self) -> Option<&BranchFlow> {
        self.branch_flows.iter().max_by(|a, b| {
            a.loading_percent.partial_cmp(&b.loading_percent).unwrap()
        })
    }
}

Core API

Usage Examples

Basic Load Flow Calculation

use eneros_powerflow::{PowerFlowSolver, PowerFlowMethod};
use eneros_topology::{NetworkGraph, Bus, Branch, BranchParams, BusType};
use eneros_core::Voltage;

let mut net = NetworkGraph::new();
net.add_bus(Bus::new("1", BusType::Slack, Voltage::new(1.05, 0.0)));
net.add_bus(Bus::new("2", BusType::PQ,   Voltage::new(1.00, 0.0)));
net.add_bus(Bus::new("3", BusType::PV,   Voltage::new(1.02, 0.0)));
net.add_branch(Branch::new("1", "2", BranchParams::line(0.1, 0.2, 0.02)));
net.add_branch(Branch::new("2", "3", BranchParams::line(0.05, 0.1, 0.01)));
net.add_branch(Branch::new("1", "3", BranchParams::line(0.08, 0.15, 0.015)));

let solver = PowerFlowSolver::new()
    .method(PowerFlowMethod::NewtonRaphson)
    .tolerance(1e-8)
    .max_iterations(50);

let result = solver.solve(&net)?;

println!(
    "Converged: {} Iterations: {} Residual: {:.2e} Duration: {:?}",
    result.converged, result.iterations, result.residual, result.duration
);

for bus in &result.buses {
    println!(
        "Bus {}: V={:.4} θ={:.4}rad P={:.4} Q={:.4}",
        bus.id, bus.voltage.magnitude.0, bus.voltage.angle.as_radians(),
        bus.p_pu, bus.q_pu
    );
}

Branch Flow and Loading

let result = solver.solve(&net)?;

for flow in &result.branch_flows {
    println!(
        "Branch {}-{}: P={:.4} Q={:.4} loading={:.1}% loss={:.4}",
        flow.from, flow.to,
        flow.p_from_pu, flow.q_from_pu,
        flow.loading_percent, flow.loss_pu
    );
}

println!("Total network loss: {:.6} pu", result.total_loss_pu());

if let Some(max) = result.max_loading() {
    println!("Maximum loading branch: {}-{} ({:.1}%)", max.from, max.to, max.loading_percent);
}

DC Power Flow Fast Estimation

use eneros_powerflow::PowerFlowMethod;

let dc_solver = PowerFlowSolver::new()
    .method(PowerFlowMethod::DcPowerFlow);

let dc_result = dc_solver.solve(&net)?;
println!("DC Power Flow duration: {:?}", dc_result.duration);

Fast Decoupled Solving

let fd_solver = PowerFlowSolver::new()
    .method(PowerFlowMethod::FastDecoupled)
    .tolerance(1e-6)
    .max_iterations(30);

let fd_result = fd_solver.solve(&net)?;
println!(
    "Decoupled Load Flow: converged={} iterations={} duration={:?}",
    fd_result.converged, fd_result.iterations, fd_result.duration
);

Result Cache Invalidation

let solver = PowerFlowSolver::new().enable_cache(true);

let r1 = solver.solve(&net)?;
let r2 = solver.solve(&net)?;
println!("Second call reuses cache: {:?}", r2.duration);

// Cache automatically invalidates after topology changes
// net.add_branch(...);
let r3 = solver.solve(&net)?;
println!("Recalculated after topology change: {:?}", r3.duration);

Convergence and Diagnostics

The Newton-Raphson method solves the linear system J · Δx = -ΔS at each iteration, where:

• J: Jacobian matrix (order 2N-2, sparse)
• Δx: Voltage correction [Δθ, ΔV]
• ΔS: Power residual [ΔP, ΔQ]

Convergence condition: max(|ΔP|, |ΔQ|) < tolerance

Divergence diagnostics return the residual at each iteration for troubleshooting:

match solver.solve(&net) {
    Ok(result) if !result.converged => {
        println!(
            "Not converged: iterated {} times, final residual {:.2e}",
            result.iterations, result.residual
        );
    }
    Ok(result) => println!("Converged at iteration {}", result.iterations),
    Err(e) => println!("Solve error: {}", e),
}

Performance Metrics

Test environment: 4 cores / 8GB / Ubuntu 22.04 / Rust 1.78 release build.

YBus Construction Performance

Configuration Parameter Table

Dependencies

Own Dependencies

[dependencies]
eneros-core = { path = "../core", version = "0.47.0" }
eneros-linalg = { path = "../linalg", version = "0.47.0" }
eneros-topology = { path = "../topology", version = "0.47.0" }
num-complex = "0.4"
tracing = "0.1"

Depended On By

• eneros-analysis — OPF and state estimation use Load Flow results as initial values
• eneros-constraint — Validates whether Load Flow results exceed limits
• eneros-twin — Twin triggers Load Flow in What-If analysis
• eneros-simulator — Scenario simulator calls Load Flow
• eneros-network — Topology-Power Flow unified pipeline

Version and Compatibility

Related Documentation

• Load Flow in Practice — Tutorial
• eneros-topology — Topology modeling
• eneros-constraint — Load Flow result validation
• eneros-linalg — Sparse matrices
• Crate Index — Complete list