Optimal power flow#

This tutorial uses one small network to show four optimal power flow studies: economic dispatch, preventive security-constrained dispatch, corrective security-constrained dispatch, and multi-period dispatch. Each study reports locational marginal prices (LMP) and separates them into energy, congestion, and security components.

The example is built in memory so the notebook runs without external files. In practice, load a case with load_case, which the file-name suffix dispatches. MATPOWER files can include cost data. You can also provide costs in the study definition.

Building the example#

A three-bus network is used. Bus 1 carries a cheap generator and is the reference; bus 2 carries an expensive generator; bus 3 carries the load. Bus 1 and bus 3 are joined by two parallel circuits, so the loss of one circuit still leaves the load connected. The generator cost curves are affine, so the marginal cost is a constant per unit (10 and 40 $/MWh).

The ramp value gives the megawatt change available to each generator for corrective redispatch or movement between periods.

from gpf import (Bus, Branch, Generator, GenCost, Load, GpfCase,
                 solve_dc_opf, solve_dc_scopf, solve_sced)


def build_case(ramp=None):
    """A three-bus test system with two parallel 1-3 circuits (60 MW normal,
    90 MW emergency). ``ramp`` sets each unit's MW recourse/inter-period budget."""
    buses = [Bus(ibus=1, ide=3, baskv=138.0),   # reference, cheap generator
             Bus(ibus=2, ide=2, baskv=138.0),   # expensive generator
             Bus(ibus=3, ide=1, baskv=138.0)]   # load
    gens = [Generator(ibus=1, machid="1", pg=100.0, pt=250.0, pb=0.0,
                      qt=200, qb=-200, vs=1.0, stat=1),
            Generator(ibus=2, machid="1", pg=50.0, pt=150.0, pb=0.0,
                      qt=200, qb=-200, vs=1.0, stat=1)]
    branches = [
        Branch(ibus=1, jbus=2, ckt="1", r=0.005, x=0.1, rate1=200, rate3=200, stat=1),
        Branch(ibus=1, jbus=3, ckt="1", r=0.005, x=0.1, rate1=60, rate3=90, stat=1),
        Branch(ibus=1, jbus=3, ckt="2", r=0.005, x=0.1, rate1=60, rate3=90, stat=1),
        Branch(ibus=2, jbus=3, ckt="1", r=0.005, x=0.1, rate1=200, rate3=200, stat=1),
    ]
    loads = [Load(ibus=3, loadid="1", pl=180.0, ql=0.0, stat=1)]
    kw = dict(ramp=ramp) if ramp is not None else {}
    costs = [GenCost(ibus=1, machid="1", model=2, cost=[10.0, 0.0], **kw),
             GenCost(ibus=2, machid="1", model=2, cost=[40.0, 0.0], **kw)]
    return GpfCase(sbase=100.0, rev=33, basfrq=60.0, buses=buses,
                   generators=gens, branches=branches, loads=loads,
                   fixed_shunts=[], transformers=[], gen_costs=costs)


case = build_case()
len(case.buses), len(case.branches), len(case.generators), len(case.gen_costs)

(3, 4, 2, 2)

Base dispatch and locational prices#

The base DC-OPF picks the least-cost dispatch on the linear network. Serving the whole 180 MW load from the cheap unit would overload the two 1-3 circuits. The dispatch therefore uses 120 MW from the cheap unit and 60 MW from the expensive unit. Each circuit reaches its 60 MW limit, and the total cost is 3600 $/hr.

base = solve_dc_opf(case)
assert abs(base.objective - 3600.0) < 1e-4
round(base.objective, 1), {g: round(mw, 1) for g, mw in base.pg.items()}

(3600.0, {(1, '1'): 120.0, (2, '1'): 60.0})

Both parallel circuits reach their 60 MW limits. The 1-2 tie carries no flow.

{br: round(mw, 1) for br, mw in base.flow.items()}

{(1, 2, '1'): 0.0, (1, 3, '1'): 60.0, (1, 3, '2'): 60.0, (2, 3, '1'): 60.0}

The energy component equals the reference-bus price of 10 $/MWh. Congestion adds 30 $/MWh at bus 2 and 60 $/MWh at bus 3. The load at bus 3 therefore pays 70 $/MWh. Each row below shows energy, congestion, and total price.

{b: (base.lmp_energy, round(base.lmp_congestion[b], 1), round(base.lmp[b], 1))
 for b in (1, 2, 3)}

{1: ({1: 10.0, 2: 10.0, 3: 10.0}, 0.0, 10.0),
 2: ({1: 10.0, 2: 10.0, 3: 10.0}, 30.0, 40.0),
 3: ({1: 10.0, 2: 10.0, 3: 10.0}, 60.0, 70.0)}

Preventive security-constrained dispatch#

The security-constrained DC-OPF keeps branch flows within emergency ratings after each single outage. If one 1-3 circuit opens, the other circuit must carry more flow. The secure dispatch therefore uses 90 MW from each generator. This raises the cost to 4500 $/hr.

prev = solve_dc_scopf(case, contingency_rating="emergency")
assert prev.converged
assert prev.contingency_rating == "emergency"
assert abs(prev.objective - 4500.0) < 1e-4
round(prev.objective, 1), prev.converged, {g: round(mw, 1) for g, mw in prev.pg.items()}

(4500.0, True, {(1, '1'): 90.0, (2, '1'): 90.0})

The loss of the second 1-3 circuit is binding. It loads the surviving circuit to 100 percent of its emergency rating and has a shadow price of 90 $/MWh.

[(d['contingency'], d['monitored'], round(d['shadow_price'], 1))
 for d in prev.binding]

[('1-3-2', (1, 3, '1'), 90.0)]

The load price remains 70 $/MWh, but the N-1 security limit now causes the price difference. Each row shows energy, congestion, security, and total price.

{b: (prev.lmp_energy, round(prev.lmp_congestion[b], 1),
     round(prev.lmp_security[b], 1), round(prev.lmp[b], 1))
 for b in (1, 2, 3)}

{1: ({1: 10.0, 2: 10.0, 3: 10.0}, 0.0, 0.0, 10.0),
 2: ({1: 10.0, 2: 10.0, 3: 10.0}, 0.0, 30.0, 40.0),
 3: ({1: 10.0, 2: 10.0, 3: 10.0}, 0.0, 60.0, 70.0)}

The full result renders as a report.

prev

Security-constrained dispatch (preventive DC-SCOPF)

Objective: $4,500.00/hr
Security: secure
Energy price: $10.00/MWh
Contingencies: 4 evaluated, 0 omitted by limit, 0 failed

Congested-bus prices ($/MWh)

bus	total	energy	congestion	security	loss
3	70.00	10.00	0.00	60.00	0.00
2	40.00	10.00	0.00	30.00	0.00

Binding contingencies

contingency	monitored	loading %	shadow price
1-3-2	1-3-1	100.0	90.00

Corrective recourse#

Corrective SCOPF permits generator redispatch after an outage. Each generator can move by 50 MW. This flexibility lets the base case use the economic 120/60 MW dispatch at a cost of 3600 $/hr. If an outage occurs, the post-contingency dispatch remains within each generator's ramp limit.

corr = solve_dc_scopf(build_case(ramp=50.0), corrective=True,
                      contingency_rating="emergency")
assert corr.converged
assert abs(corr.objective - base.objective) < 1e-4
round(corr.objective, 1), corr.converged, {g: round(mw, 1) for g, mw in corr.pg.items()}

(3600.0, True, {(1, '1'): 120.0, (2, '1'): 60.0})

For the loss of either 1-3 circuit, the cheap generator reduces output by 30 MW and the expensive generator increases output by 30 MW. Both changes are within the 50 MW limits.

{lbl: {g: round(mw, 1) for g, mw in acts.items()}
 for lbl, acts in corr.remedial.items()}

{'1-3-2': {(1, '1'): -30.0, (2, '1'): 30.0},
 '1-3-1': {(1, '1'): -30.0, (2, '1'): 30.0}}

The result reports the 50 MW ramp limit used for each generator.

corr.stats['ramp_source']

{'explicit': 2, 'ramp_30': 0, 'ramp_10': 0, 'ramp_agc': 0, 'headroom': 0}

Multi-period dispatch with a binding ramp#

Security-constrained economic dispatch (SCED) schedules generation across several periods and applies each generator's ramp limit between periods. This example uses half load in the first period and full load in the second. The expensive generator cannot move from its first-period economic output to the required 60 MW within its 50 MW ramp limit. SCED therefore starts it above its first-period economic output. The total cost is 4800 $/hr.

sced = solve_sced(build_case(ramp=50.0), demand=[0.5, 1.0])
assert abs(sced.objective - 4800.0) < 1e-4
round(sced.objective, 1), [{g: round(mw, 1) for g, mw in per.items()} for per in sced.pg]

(4800.0, [{(1, '1'): 80.0, (2, '1'): 10.0}, {(1, '1'): 120.0, (2, '1'): 60.0}])

The first-period price is 10 $/MWh at each bus. In the second period, the binding ramp limit raises the load price to 130 $/MWh. A single-period dispatch at the same load would set a 70 $/MWh price. The difference reflects the limited ability to increase output between periods. Each period lists the prices at buses 1, 2, and 3.

[tuple(float(round(per[b], 1)) for b in (1, 2, 3)) for per in sced.lmp]

[(10.0, 10.0, 10.0), (10.0, 70.0, 130.0)]

The expensive generator increases output by the full 50 MW ramp limit between the two periods.

[(b['gen'], b['from_period'], b['to_period'], round(b['delta_mw'], 1))
 for b in sced.ramp_binding]

[((2, '1'), 0, 1, 50.0)]

Scenario sweeps and price distributions#

A scenario sweep shows how prices and congestion change across operating conditions. It reports average prices, price percentiles, branch binding frequencies, and path congestion spreads.

Use load_scenario, outage_scenario, and monte_carlo_load_scenarios to define changes to the base case. The 24 scenarios below vary load and include cases with the 1-2 tie out of service.

from gpf import (load_scenario, outage_scenario, monte_carlo_load_scenarios,
                 sweep, FtrPortfolio, FtrPosition)

scenarios = ([load_scenario(m, label=f"load{int(m * 100)}") for m in (0.8, 1.0, 1.2)]
             + [outage_scenario((1, 2, "1"), label="lose_tie_1-2")]
             + monte_carlo_load_scenarios(20, seed=7, sigma=0.1))
result = sweep(case, scenarios, paths=[(1, 3)])
assert result.n_scenarios == 24
assert not result.failed
result.n_scenarios, result.n_networks, result.failed

(24, 2, [])

In most scenarios, the two 1-3 circuits reach their limits. The median and 95th-percentile prices at bus 3 are therefore 70 $/MWh. The 5th-percentile price is closer to the uncongested price. The circuits bind in approximately nine out of ten scenarios, and the median 1-3 congestion spread is 60 $/MWh. The bus 3 price has a standard deviation of approximately 17 $/MWh.

summary = result.distribution.summary()
lmp3 = {s: round(summary["bus_lmp"][3][s], 1) for s in ("mean", "std", "p5", "p50", "p95")}
binding = {f"{k[0]}-{k[1]}({k[2]})": round(summary["branch"][k]["binding_frequency"], 2)
           for k in summary["branch"] if (k[0], k[1]) == (1, 3)}
spread = {s: round(summary["paths"][(1, 3)][s], 1) for s in ("mean", "p50", "p95")}
lmp3, binding, spread

({'mean': 63.8, 'std': 17.3, 'p5': 10.0, 'p50': 70.0, 'p95': 70.0},
 {'1-3(1)': 0.92, '1-3(2)': 0.92},
 {'mean': 53.8, 'p50': 60.0, 'p95': 60.0})

Next steps#

The API reference describes solve_dc_opf, solve_dc_scopf, solve_scopf, and solve_sced. It also describes the DcNetwork and DcSensitivity classes used for PTDF, LODF, and OTDF calculations.

Transmission rights explains how to hedge congestion prices. Sensitivities explains the factors used in security analysis. Transfer capability explains how to measure the available transfer margin.