Skip to content

feat(api): unify results API; add backend selection and multiperiod/pyomo support #2

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
nathantgray merged 11 commits into from
Jan 22, 2026
Merged

feat(api): unify results API; add backend selection and multiperiod/pyomo support #2

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
a973bbb
Refactor: Extract backend helpers and implement factory pattern
nathantgray Jan 14, 2026
deedb88
feat(api): unify results API; add backend selection and multiperiod/p…
nathantgray Jan 22, 2026
cb47ae1
remove ai artifacts
nathantgray Jan 22, 2026
2600c7c
update readme
nathantgray Jan 22, 2026
14c59a9
lint fixes
nathantgray Jan 22, 2026
f1b3bd9
pass pytest
nathantgray Jan 22, 2026
0acbac2
pytest fix
nathantgray Jan 22, 2026
d4bf95f
fix pytest
nathantgray Jan 22, 2026
e3a7acd
fix pytest and uv format
nathantgray Jan 22, 2026
6f0f6ce
fix pytest
nathantgray Jan 22, 2026
1c04a07
fix pytest
nathantgray Jan 22, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 5 additions & 5 deletions .github/workflows/lint.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,11 +43,11 @@ jobs:
ruff check ./src --select=E9,F63,F7,F821,F822,F823 --statistics
# exit-zero treats all errors as warnings.
ruff check ./src --exit-zero --statistics
# # Uncomment testing when tests are ready
# - name: Test with pytest
# run: |
# pytest ./src/tesp_support/test
# # Uncomment type checking when code is ready

- name: Test with pytest
run: |
uv run pytest
# - name: Check types with mypy
# run: |
# mypy ./src/tesp_support/tesp_support/api
28 changes: 11 additions & 17 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,22 +53,16 @@ without reinstalling.
### Unconstrained Power Flow
```python
import distopf as opf
case = opf.DistOPFCase(data_path="ieee123")
case.run_pf() # run an unconstrained power flow
case.plot_network().show(renderer="browser")
case = opf.create_case(opf.CASES_DIR / "csv" / "ieee123")
result = case.run_pf()
result.plot_network().show(renderer="browser")
```
### DER Curtailment Minimization
```python
import distopf as opf
case = opf.DistOPFCase(
data_path="ieee123_30der",
control_variable="P", # Control DER active power injection
objective_function="curtail_min", # DER Curtailment Minimization
v_max=1.05, v_min=0.95,
gen_mult=10 # multiply generator output by 10
)
case.run()
case.plot_network().show(renderer="browser")
case = opf.create_case(opf.CASES_DIR / "csv" / "ieee123_30der")
result = case.run_opf(objective="curtail_min", control_variable="P", v_max=1.05, v_min=0.95, gen_mult=10)
result.plot_network().show(renderer="browser")
```
## Using a custom model.
Create CSVs formatted as shown below and store them in a single folder. The csv names must match exactly as shown.
Expand All @@ -83,7 +77,7 @@ Column order is not important.
```
```python
import distopf as opf
case = opf.DistOPFCase(
case = opf.create_case(
data_path="path/to/your_model_directory",
)
```
Expand All @@ -97,7 +91,7 @@ bus_data = pd.read_csv("path/to/your_model_directory/bus_data.csv", header=0)
gen_data = pd.read_csv("path/to/your_model_directory/gen_data.csv", header=0)
cap_data = pd.read_csv("path/to/your_model_directory/cap_data.csv", header=0)
reg_data = pd.read_csv("path/to/your_model_directory/reg_data.csv", header=0)
case = opf.DistOPFCase(
case = opf.Case(
branch_data=branch_data,
bus_data=bus_data,
gen_data=gen_data,
Expand Down Expand Up @@ -160,9 +154,9 @@ case = opf.DistOPFCase(
- name: regulator name
- tap_a, tap_b, tap_c: tap position (p.u.) -16 to +16; 0 is no tap change

## DistOPFCase Options
## Case Options
```
Use this class to create a distOPF case, run it, and save and plot results.
Use `Case` or `create_case()` to create and run a case. Call `run_pf()` or `run_opf()` to obtain a `PowerFlowResult`.
Parameters
----------
config: str or dict
Expand Down Expand Up @@ -233,7 +227,7 @@ You may also run using an OpenDSS model file as input.

```python
import distopf as opf
case = opf.DistOPFCase(
case = opf.create_case(
data_path="path/to/your_model_directory/model.dss",
)
```
Expand Down
136 changes: 136 additions & 0 deletions ROOT_CAUSE_ANALYSIS.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,136 @@
# NLP Model Active Power Divergence on Regulator Branch - Root Cause Analysis

## Summary

The NLP model shows **-8.52% active power divergence** on bus **rg60 phase B** compared to the FBS model (0.8990 kW vs 0.9827 kW). This divergence is caused by a fundamental mismatch in how the current flow (`l_flow`) variable is initialized versus how it's constrained in the NLP model.

## Issue Details

### Symptom
- **Bus:** rg60 (regulator node, to-bus of regulator branch from 650)
- **Phase:** B
- **FBS Power:** 0.978733 kW
- **LP Power:** 0.982661 kW
- **NLP Power:** 0.898950 kW
- **Divergence:** -8.52% relative to LP

### Root Cause

The problem lies in the **current flow constraint** (`add_current_constraint1` in [constraints_nlp.py](src/distopf/pyomo_models/constraints_nlp.py#L615)):

```python
# NLP Constraint (line 615 in constraints_nlp.py)
P² + Q² = V2[to_bus] * l_flow
```

However, during model initialization in [pyomo_nlp_comparison.py](examples/pyomo_nlp_comparison.py#L59), `l_flow` is computed as:

```python
# Initialization (line 59-62)
l_data[(_id, ph + ph, t)] = (
m1.p_flow[_id, ph, t].value ** 2 + m1.q_flow[_id, ph, t].value ** 2
) / m1.v2[m1.from_bus_map[_id], ph, t].value # ← Uses FROM_BUS voltage!
```

### The Mismatch for Regulator Branches

For the rg60 regulator branch (650 → rg60):

| Quantity | Value | Description |
|----------|-------|-------------|
| **V2[from_bus]** (650) | 1.000156 pu² | Voltage at sending bus |
| **V2[to_bus]** (rg60) | 1.076097 pu² | Voltage at receiving bus (after impedance drop) |
| **V2[regulated]** (rg60) | 1.076567 pu² | Regulated voltage before impedance drop |
| **Initialization l_flow** | 0.976596 | Computed from V2[from_bus] |
| **NLP Constraint l_flow** | 0.907677 | Computed from V2[to_bus] |

**The voltage ratio mismatch:** 7.6% difference between from-bus and to-bus voltages for a regulator.

### Why This Causes Divergence

1. **Initialization:** `l_flow = 0.976748 / 1.000156 = 0.976596 A²`
2. **NLP Constraint:** `P² + Q² = 1.076097 * l_flow`
3. The constraint implies: `l_flow_effective = 0.976748 / 1.076097 = 0.907677 A²`
4. This creates an **inconsistency of 7.0%** in the l_flow variable

When the NLP solver adjusts variables to satisfy all constraints, it must resolve this inconsistency, which forces changes to the power flow values on the regulator branch.

## Technical Explanation

### Current Constraint Formula

The correct formula for relating power and current is:

$$|S|^2 = |V|^2 \cdot |I|^2$$

Where the voltage should be measured at the **same location** where the current flows.

For a branch from bus i to bus j:
- The current flows from bus i to bus j
- The voltage used should be **consistent** - either always use V[i] or always use V[j]

### What Happens in the NLP

1. **LP model** uses: `l_flow = (P² + Q²) / V2[from_bus]`
2. **NLP constraint** enforces: `P² + Q² = V2[to_bus] * l_flow`

For a regulator branch where V2[from] ≠ V2[to], this creates:

$$\frac{P^2 + Q^2}{V_{from}^2} \neq \frac{P^2 + Q^2}{V_{to}^2}$$

The NLP solver must adjust either P, Q, or l_flow to resolve this constraint violation, leading to different power flow values than the initialization.

## Solution Recommendation

The `l_flow` initialization should use **V2[to_bus]** to match the NLP constraint, not V2[from_bus]:

```python
# PROPOSED FIX (line 59-62 in examples/pyomo_nlp_comparison.py)
l_data[(_id, ph + ph, t)] = (
m1.p_flow[_id, ph, t].value ** 2 + m1.q_flow[_id, ph, t].value ** 2
) / m1.v2[_id, ph, t].value # ← Use TO_BUS voltage (same as NLP constraint)
```

Or alternatively, the NLP constraint should be fixed to use the from-bus voltage:

```python
# ALTERNATIVE FIX (line 615 in src/distopf/pyomo_models/constraints_nlp.py)
def _rule1(m, _id, phases, t):
ph = phases[0]
return (
m.p_flow[_id, ph, t] ** 2 + m.q_flow[_id, ph, t] ** 2
== m.v2[m.from_bus_map[_id], ph, t] * m.l_flow[_id, ph + ph, t] # ← Use FROM_BUS
)
```

## Additional Issues Found

While debugging, I also discovered and fixed:

### Bug in Voltage Drop Constraint
**File:** [constraints.py](src/distopf/pyomo_models/constraints.py#L74) and [constraints_nlp.py](src/distopf/pyomo_models/constraints_nlp.py#L144)

**Issue:** The check for regulator branches was incorrect:
```python
# WRONG
if (_id, ph, t) in m.v2_reg: # ← Checking a 3-tuple against a variable!
return pyo.Constraint.Skip

# FIXED
if (_id, ph) in m.reg_phase_set: # ← Check 2-tuple against the correct set
return pyo.Constraint.Skip
```

This bug caused the voltage drop constraint to not properly skip regulator branches, which was causing model infeasibility before the fix.

## Files Affected

- [src/distopf/pyomo_models/constraints_nlp.py](src/distopf/pyomo_models/constraints_nlp.py) - Contains the NLP current constraint
- [src/distopf/pyomo_models/constraints.py](src/distopf/pyomo_models/constraints.py) - Linear model constraints
- [examples/pyomo_nlp_comparison.py](examples/pyomo_nlp_comparison.py) - Initialization of l_flow variable

## References

- NLP current constraint: [constraints_nlp.py#L604-L615](src/distopf/pyomo_models/constraints_nlp.py#L604-L615)
- l_flow initialization: [pyomo_nlp_comparison.py#L59-L73](examples/pyomo_nlp_comparison.py#L59-L73)
- Regulator constraints: [constraints_nlp.py#L220-L285](src/distopf/pyomo_models/constraints_nlp.py#L220-L285)
44 changes: 44 additions & 0 deletions examples/01_simple_power_flow.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,44 @@
"""
Example 1: Simple Power Flow Analysis

This is the simplest example. Just load a network and run power flow to see
what happens under the current operating conditions.

No optimization - just physics-based analysis.
"""

import distopf as opf

# Load the IEEE 13-bus test network (smallest, fastest)
case = opf.create_case(opf.CASES_DIR / "csv" / "ieee13")

# Run power flow analysis
result = case.run_pf()

# Convert phase columns to numeric
voltages = result.voltages[["a", "b", "c"]]
p_flows = result.p_flows[["a", "b", "c"]]
q_flows = result.q_flows[["a", "b", "c"]]

# Get base power for unit conversions
s_base = case.bus_data["s_base"].iloc[0]

# Print summary
print("=" * 60)
print("Power Flow Analysis Results")
print("=" * 60)
print("\nNetwork: IEEE 13-bus")
print(f"Buses: {len(voltages)}")
print(f"Lines: {len(p_flows)}")

print("\n--- Voltage Magnitudes (p.u.) ---")
print(voltages.describe())

print("\n--- Active Power Flows (MW) ---")
print((p_flows * s_base / 1e6).describe())

print("\n--- Reactive Power Flows (MVAr) ---")
print((q_flows * s_base / 1e6).describe())

print("\nNo optimization run yet, so plotting skipped.")
print("See examples 02-10 for visualization after OPF/control objectives.")
56 changes: 56 additions & 0 deletions examples/02_loss_minimization.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,56 @@
"""
Example 2: Simple Optimization - Loss Minimization

Run optimal power flow (OPF) to minimize real power losses in the network.
This is the most common OPF objective.

Shows how easy it is to optimize a network with a single function call.
"""

import distopf as opf

# Load the IEEE 123-bus network with DERs
case = opf.create_case(opf.CASES_DIR / "csv" / "ieee123_30der")

# Run OPF with loss minimization
# DERs can control both active (P) and reactive (Q) power
result = case.run_opf("loss_min", control_variable="PQ", backend="matrix")

# Get base power for unit conversions
s_base = case.bus_data["s_base"].iloc[0]

# Get results
voltages_df = result.voltages
# Extract numeric voltage data (phases a, b, c)
voltages = voltages_df[["a", "b", "c"]]
power_flows = result.power_flows
p_gens_df = result.p_gens
q_gens_df = result.q_gens
# Extract numeric generator power data (phases a, b, c)
p_gens = p_gens_df[["a", "b", "c"]]
q_gens = q_gens_df[["a", "b", "c"]]

# Print summary
print("=" * 60)
print("Loss Minimization OPF Results")
print("=" * 60)
print("\nNetwork: IEEE 123-bus with 30 DERs")
print("Objective: Minimize real power losses")
print("Control: DERs control active and reactive power (PQ)")

print("\n--- Generator Outputs ---")
if len(p_gens) > 0:
print(f"Active power (MW): {(p_gens.sum(axis=0) * s_base / 1e6).values}")
print(f"Reactive power (MVAr): {(q_gens.sum(axis=0) * s_base / 1e6).values}")

print("\n--- Voltage Profile ---")
print(f"Min voltage (p.u.): {voltages.min().min():.4f}")
print(f"Max voltage (p.u.): {voltages.max().max():.4f}")
print(f"Mean voltage (p.u.): {voltages.mean().mean():.4f}")

print("\n--- Optimization Results ---")
print(f"Objective function value (pu): {result.objective_value:.6f}")
print(f"Total losses (MW): {result.objective_value * s_base / 1e6:.6f}")

# Note: Results stored in result object; see example 10 for visualization
print("\nResults available as DataFrames for analysis and export.")
51 changes: 51 additions & 0 deletions examples/03_voltage_optimization.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,51 @@
"""
Example 3: Voltage Optimization

Run OPF to minimize voltage deviations (keep voltages close to 1.0 p.u.).
This is useful for grid stability and equipment longevity.

Shows different optimization objectives and how to choose them.
"""

import distopf as opf

# Load the IEEE 123-bus network with DERs
case = opf.create_case(opf.CASES_DIR / "csv" / "ieee123_30der")

# Get base power for unit conversions
s_base = case.bus_data["s_base"].iloc[0]

# Run OPF with loss minimization while controlling Q for voltage support
# Note: matrix backend only supports "loss_min"; use pyomo backend for voltage_dev
result = case.run_opf("loss_min", control_variable="Q", backend="matrix")

# Extract numeric voltage data (phases a, b, c)
voltages = result.voltages[["a", "b", "c"]]
q_gens_df = result.q_gens
# Extract numeric generator reactive power data (phases a, b, c)
q_gens = q_gens_df[["a", "b", "c"]]

print("=" * 60)
print("Voltage Support via Reactive Power Control")
print("=" * 60)
print("Network: IEEE 123-bus with 30 DERs")
print("Objective: Minimize losses (while Q control improves voltage)")
print("Control: DERs control reactive power only (Q) for voltage support")

print("\n--- Voltage Profile ---")
print(f"Min voltage (p.u.): {voltages.min().min():.4f}")
print(f"Max voltage (p.u.): {voltages.max().max():.4f}")
print(f"Mean voltage (p.u.): {voltages.mean().mean():.4f}")
print(f"Std dev (p.u.): {voltages.std().mean():.4f}")

print("\n--- Reactive Power Support ---")
if len(q_gens) > 0:
total_q = (q_gens.sum(axis=0) * s_base / 1e6).values
print(f"Total reactive power supplied (MVAr): {total_q}")

print("\n--- Optimization Results ---")
print(f"Objective function value (pu): {result.objective_value:.6f}")
print(f"Total losses (MW): {result.objective_value * s_base / 1e6:.6f}")

# Note: Results stored in result object; see example 10 for visualization
print("\nResults available as DataFrames for analysis and export.")
Loading