Calculating locational marginal pricing in Python

When a computed nodal price disagrees with the ISO’s published LMP by a few cents, the culprit is almost always a mishandled dual variable — a congestion shadow price with the wrong sign, or a loss component quantized differently than the market operator’s clearing engine — and this page shows how to solve the DC optimal power flow, extract those duals, and reassemble a settlement-exact locational marginal price in Python. It is a step in Pricing Logic Implementation, the layer that turns market clearing results into auditable financial obligations before any volume is monetized.

The diagram below traces the solve-to-LMP flow this page implements: a DCOPF solve yields dual variables, which are extracted and combined into the three additive LMP components that sum to the nodal price.

Mathematical Decomposition #

The LMP formulation follows a standardized additive structure derived from the dual variables of a security-constrained economic dispatch (SCED) or DC optimal power flow (DCOPF) model. At any market interval, the price at node \(i\) decomposes as:

$$LMP_i = \lambda_{energy} + \sum_{k} \mu_k \cdot \text{PTDF}{k,i} + \lambda{energy} \cdot \text{MLF}_i$$

where \(\lambda_{energy}\) is the shadow price of the global power-balance constraint (the system-wide energy component common to every node), \(\mu_k\) is the dual variable of binding transmission constraint \(k\) (the congestion component), \(\text{PTDF}_{k,i}\) is the power transfer distribution factor mapping node-\(i\) injections to line-\(k\) flows, and \(\text{MLF}_i\) is the marginal loss factor capturing incremental transmission losses relative to the reference bus. This is the node-level realization of the same three-part decomposition the broader pricing engine defines as \(LMP_n = \lambda + \mu_n + \nu_n\).

The optimization minimizes total generation cost subject to physical and operational limits. The constraint matrix must enforce strict node-index alignment between generator bids, load forecasts, and transmission limits; misaligned topology snapshots, stale shift-factor files, or unvalidated telemetry gaps are the most frequent sources of settlement discrepancies and downstream billing disputes. The interval grid these prices resolve against is produced by the Settlement Cycle Mapping engine, and the field-level shape of each feed is enforced by the ISO/RTO Data Format Standards.

Prerequisites #

• Python packages: cvxpy (convex modelling layer), numpy, scipy for scipy.sparse matrix operations, and a conic solver — clarabel is used below; ecos or a licensed mosek are drop-in alternatives via the solver= argument.
• Data dependencies: a generator bid vector ($/MWh per node), a sparse PTDF / shift-factor matrix aligned to the same node ordering, a per-line thermal limit vector, a system load-balance scalar for the interval, and a marginal-loss-factor vector referenced to the market’s slack bus. Every array must share one canonical node index; the PTDF file and the bid vector drifting apart by even one node silently mis-prices the whole interval.
• Permissions / access: read-only credentials for the market operator’s data portal (PJM Data Miner, MISO Market Reports, CAISO OASIS, or ERCOT MIS) to pull the published nodal LMPs used in the reconciliation step. No write scope is required — this engine computes and verifies, it does not submit.

Implementation #

A production LMP solver must prioritize memory efficiency, deterministic execution, and explicit error boundaries. Leveraging cvxpy for the convex program and scipy.sparse for the PTDF keeps it scalable across large transmission networks. The implementation below performs deterministic dual extraction, audit-safe logging, and settlement-compliant rounding, deferring the Decimal cast to the final quantization so the solver core stays in fast floating point while every persisted value is exact.

import cvxpy as cp
import numpy as np
import scipy.sparse as sp
import logging
from decimal import Decimal, ROUND_HALF_UP
from typing import Optional, Dict, Any

logger = logging.getLogger("lmp_engine")
logger.setLevel(logging.INFO)

def solve_lmp(
    bid_vector: np.ndarray,          # $/MWh generation offer per node
    ptdf_sparse: sp.spmatrix,        # shape (n_lines, n_nodes), node order == bid order
    line_limits: np.ndarray,         # thermal limit per line (MW)
    load_balance: float,             # system load for the interval (normalized MW)
    loss_factors: np.ndarray,        # marginal loss factor per node, slack-referenced
    max_iterations: int = 10000,
) -> Optional[Dict[str, Any]]:
    """
    Compute LMP components via DCOPF with explicit dual extraction.
    Returns the energy, congestion, loss, and total LMP vectors plus lineage.
    """
    n_nodes = len(bid_vector)
    n_lines = ptdf_sparse.shape[0]
    if ptdf_sparse.shape[1] != n_nodes or len(loss_factors) != n_nodes:
        raise ValueError("PTDF columns, bid vector, and loss factors must share the node index")

    # Decision variable: generation dispatch at each node
    generation = cp.Variable(n_nodes, name="generation")

    # Objective: minimize total generation cost
    objective = cp.Minimize(bid_vector @ generation)

    # Constraints. The two transmission rows bound flow in each direction; their
    # duals are netted below so congestion reflects both forward and reverse limits.
    constraints = [
        generation >= 0,
        generation <= 1.0,                       # normalized capacity (scale to MW in prod)
        cp.sum(generation) == load_balance,      # constraints[2]: system power balance
        ptdf_sparse @ generation <= line_limits, # constraints[3]: forward flow limit
        ptdf_sparse @ generation >= -line_limits,# constraints[4]: reverse flow limit
    ]

    prob = cp.Problem(objective, constraints)

    try:
        prob.solve(solver=cp.CLARABEL, verbose=False, max_iters=max_iterations)
    except cp.SolverError as exc:
        logger.error("SCED solver failed for interval: %s", exc)
        return None

    if prob.status not in (cp.OPTIMAL, cp.OPTIMAL_INACCURATE):
        logger.warning("Suboptimal solver status: %s", prob.status)
        return None

    # --- Dual extraction (shadow prices) ---
    # Energy price: the single dual on the system power-balance equality.
    lambda_energy = float(constraints[2].dual_value)

    # Net the non-negative duals of the forward (<=) and reverse (>=) flow limits so
    # a line binding in either direction contributes the correctly signed congestion
    # price; both are ~zero when the line is uncongested.
    mu_upper = constraints[3].dual_value
    mu_lower = constraints[4].dual_value
    mu_congestion = mu_upper - mu_lower

    # --- Recombine into the three additive LMP components ---
    energy_component = np.full(n_nodes, lambda_energy)
    congestion_component = ptdf_sparse.T @ mu_congestion   # map line duals back to nodes
    loss_component = energy_component * loss_factors
    total_lmp = energy_component + congestion_component + loss_component

    # Settlement-compliant rounding: cast to Decimal only at the boundary.
    def round_settlement(arr: np.ndarray) -> np.ndarray:
        return np.array([
            float(Decimal(str(val)).quantize(Decimal("0.01"), rounding=ROUND_HALF_UP))
            for val in arr
        ])

    return {
        "energy": round_settlement(energy_component),
        "congestion": round_settlement(congestion_component),
        "loss": round_settlement(loss_component),
        "total_lmp": round_settlement(total_lmp),
        "solver_status": prob.status,
        "dual_energy": lambda_energy,
        "binding_lines": np.where(np.abs(mu_congestion) > 1e-6)[0].tolist(),
    }

Dual extraction, sign handling, and rounding #

Extracting duals correctly is where most implementations break. constraints[2].dual_value is the system energy price. Congestion is the net of the non-negative duals on the forward (constraints[3]) and reverse (constraints[4]) flow limits, so a line binding in either direction yields the correctly signed shadow price; taking only the upper dual drops every counter-flow constraint and understates congestion at import-constrained nodes. The transpose ptdf_sparse.T @ mu_congestion then maps line-level congestion costs back to nodal prices without ever materializing a dense matrix.

Settlement reconciliation demands strict precision. Floating-point artifacts introduce basis risk during clearing, so — as across the whole Settlement Calculation & Validation Engines framework — every persisted price uses Python’s decimal module with ROUND_HALF_UP to enforce ISO-compliant rounding rather than binary-float division.

Verification steps #

Confirm the output is correct before any price crosses into the ledger:

1. Shape and keys. solve_lmp returns a dict whose energy, congestion, loss, and total_lmp arrays each have length n_nodes. A length mismatch means the PTDF and bid vector fell out of node alignment.
2. Additive identity. At every node, total_lmp must equal energy + congestion + loss to the quantized cent — the single most important check in the pricing layer.
3. Reconciliation diff against the published LMP. Join the computed total against the ISO’s published nodal price on (node_id, settlement_interval) and assert the residual is within tolerance.

def verify_lmp(result: dict, published_lmp: np.ndarray,
               tol_cents: float = 0.01) -> None:
    """Assert additive identity and reconcile against the ISO's published LMP."""
    reconstructed = result["energy"] + result["congestion"] + result["loss"]
    additive_residual = np.abs(reconstructed - result["total_lmp"]).max()
    assert additive_residual <= tol_cents, f"components do not sum: {additive_residual}"

    recon_diff = np.abs(result["total_lmp"] - published_lmp)
    breaches = np.where(recon_diff > tol_cents)[0].tolist()
    if breaches:
        raise ValueError(f"nodes off published LMP beyond tolerance: {breaches}")
    logger.info("LMP reconciliation clean across %d nodes", len(published_lmp))

Automated reconciliation must flag deviations beyond predefined bands (typically ±$0.01/MWh for energy, ±$0.05/MWh for congestion) and route to an exception workflow for manual review or topology re-validation. For an immutable audit trail, hash the solver inputs and the four output vectors together (e.g. a SHA-256 over the concatenated bytes) and store the digest with the market run identifier, so a replay of a historical interval reproduces both the price and its provenance.

Compliance note #

This implementation must be validated against the market operator’s tariff-defined pricing methodology: PJM Manual 11 / Manual 28, MISO BPM-002, CAISO’s Business Practice Manual for Market Operations, or ERCOT Nodal Protocols Sections 4 and 6, each of which fixes the three-part decomposition, the loss-factor reference, and the rounding convention the code above must mirror. Under FERC’s Open Access rules and ISO/RTO market-design standards, every priced charge must be traceable to a metered interval, a market run identifier, and a published price version — so log the solver inputs, dual outputs, binding constraints, and rounding operations for every interval. Deterministic fallback pricing (reference-bus pricing or archived-interval interpolation) must replace a failed SCED solve rather than a silent default, and PTDF/MLF snapshots must be versioned to the EMS/SCADA state-estimator run they came from, never interpolated across a topology change.

Frequently Asked Questions #

Why does my computed LMP not match the ISO’s published price? #

The three usual causes are a sign error on a congestion dual, a loss component quantized with a different rounding mode than the market operator uses, and a stale PTDF or MLF file that no longer matches the topology of the published interval. Net the forward and reverse flow-limit duals so counter-flow constraints keep their sign, quantize every component with the tariff’s rounding rule, and confirm the PTDF snapshot is versioned to the same state-estimator run as the published LMP before treating a residual as a code bug.

How do I get the congestion component out of a DCOPF dual? #

The congestion component is the line-level dual mapped back to nodes by the PTDF transpose. Extract the non-negative duals on both the forward (<=) and reverse (>=) transmission-limit constraints, subtract the reverse dual from the forward dual to recover a signed per-line shadow price, and compute ptdf_sparse.T @ mu_congestion. Both duals are effectively zero when a line is uncongested, so only binding constraints contribute to nodal congestion.

Which solver should I use for a production LMP DCOPF? #

For a linear DCOPF, an interior-point conic solver such as Clarabel gives deterministic, well-conditioned duals and is what the code above targets; ECOS is a lightweight open alternative, and MOSEK is common where a commercial license and warm-starting are already in place. Whichever you pick, pin the solver version, check prob.status for OPTIMAL before reading any dual, and treat OPTIMAL_INACCURATE as a signal to tighten tolerances rather than to trust the shadow prices blindly.

Related #

• Mapping transmission loss factors to settlement nodes
• Automating imbalance allocation for gas trades
• Pricing Logic Implementation
• Threshold Tuning & Alerts