Skip to content

OPF formulation work #84

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
bmeyers merged 3 commits into from
Dec 16, 2025
Merged

OPF formulation work #84

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
0841900
update to dc-opf notebook
bmeyers Aug 26, 2025
24a3628
opf work
bmeyers Sep 2, 2025
4e428fd
Merge branch 'main' of github.com:slacgismo/regrow into opf
bmeyers Dec 16, 2025
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
552 changes: 552 additions & 0 deletions dc-opf/SD_relax_test.py
View file
Open in desktop

Large diffs are not rendered by default.

585 changes: 585 additions & 0 deletions dc-opf/SD_relax_test_real.py
View file
Open in desktop

Large diffs are not rendered by default.

559 changes: 559 additions & 0 deletions dc-opf/SD_relax_test_v2.py
View file
Open in desktop

Large diffs are not rendered by default.

635 changes: 635 additions & 0 deletions dc-opf/SOC_relax_test.py
View file
Open in desktop

Large diffs are not rendered by default.

939 changes: 939 additions & 0 deletions dc-opf/case240.py
View file
Open in desktop

Large diffs are not rendered by default.

287 changes: 178 additions & 109 deletions dc-opf/network_flow_example_v4.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,9 +1,129 @@
import marimo

__generated_with = "0.14.12"
__generated_with = "0.14.17"
app = marimo.App(width="medium")


@app.cell
def _(cp, np):
def make_problem(incidence_mat, generator_dict, load_dict, flow_dict):
def indices_not_in_list(array_length, given_indices):
all_indices = set(range(array_length))
given_indices_set = set(given_indices)
not_in_list_indices = all_indices - given_indices_set
return list(not_in_list_indices)
_B = incidence_mat
num_buses = _B.shape[0]
num_edges = _B.shape[1]
gen_ixs = np.array(list(generator_dict.keys())) - 1
nogen_ixs = indices_not_in_list(num_buses, gen_ixs)
p_flows = cp.Variable((num_edges), name='line flows')
p_gen = cp.Variable((num_buses), nonneg=True, name='p_gen')
p_load = cp.Variable((num_buses), nonneg=True, name='p_load')
l_min = cp.Parameter(
shape=num_buses,
value=[_v['l_min'] for _k, _v in load_dict.items()],
name='l_min'
)
l_upper = cp.Parameter(
shape=num_buses,
value=[_v['l_upper'] for _k, _v in load_dict.items()],
name='l_upper'
)
# l_cost = cp.Parameter(
# shape=num_buses,
# value=[_v['cost'] for _k, _v in load_dict.items()],
# name='l_cost'
# )
l_cost = [_v['cost'] for _k, _v in load_dict.items()]
g_min = cp.Parameter(
shape=len(gen_ixs),
value=np.array([_v['p_min'] for _k, _v in generator_dict.items()]),
name='g_min'
)
g_max = cp.Parameter(
shape=len(gen_ixs),
value=np.array([_v['p_max'] for _k, _v in generator_dict.items()]),
name='g_max'
)
c0 = cp.Parameter(
shape=len(gen_ixs),
value=np.array([_v['c0'] for _k, _v in generator_dict.items()]),
name='c0'
)
c1 = cp.Parameter(
shape=len(gen_ixs),
value=np.array([_v['c1'] for _k, _v in generator_dict.items()]),
name='c1'
)
# c2 = cp.Parameter(
# shape=len(gen_ixs),
# value=np.array([_v['c2'] for _k, _v in generator_dict.items()]),
# name='c2'
# )
c2 = np.array([_v['c2'] for _k, _v in generator_dict.items()])
r = cp.Parameter(
shape=num_edges,
value=np.array([_v['f_resistance'] for _k, _v in flow_dict.items()]),
name='line resistance'
)
f_max = cp.Parameter(
shape=num_edges,
value=np.array([_v['f_max'] for _k, _v in flow_dict.items()]),
name='line capacities'
)
# generator costs
cost = cp.sum(c0 + cp.multiply(c1, p_gen[gen_ixs]) + cp.multiply(c2, cp.square(p_gen[gen_ixs])))
# load curtailment costs
cost += cp.sum(
cp.multiply(
cp.neg(p_load - l_upper),
l_cost
)
)
# line penalties
cost += cp.sum_squares(cp.multiply(r, p_flows))
constraints = [
_B @ p_flows + p_gen - p_load == 0,
p_load >= l_min,
cp.abs(p_flows) <= f_max,
p_gen[nogen_ixs] == 0,
p_gen[gen_ixs] <= g_max,
p_gen[gen_ixs] >= g_min
]
problem = cp.Problem(cp.Minimize(cost), constraints)
return problem
return (make_problem,)


@app.cell
def _(mo):
mo.md(
"""
# DC Optimal Power Flow

This notebook demonstrates the optimal flow problem, over 6 test networks provided by [pypower](https://github.com/rwl/PYPOWER).

- network with $n$ nodes, $m$ edges, given by $n \\times m$ incidence
matrix $A$
- edge power vector $p \\in \\mathbf{R}^m$,
with capacity limits $|p_j| \\leq C_j$
- node generator power $g_i$, load $l_i$,
each with lower and upper limits
- flow conservation $Af + g = l$
- generator cost function is
$\\phi_i(g_i) = a_i g_i + b_i g_i^2$, total $G= \\sum_i \\phi_i (g_i)$
- load shortfall cost is $\\psi_i(l_i) = c_i(l^\\text{tar}_i-l_i)_+$,
total is $S = \\sum_i \\psi_i(l_i)$
- total line loss is $L = \\sum_j r_j p_j^2$
- objective is $G + \\lambda L + \\mu S$

The model automatically populates parameters for the generators, loads, and lines, with the sliders set to the default values for the test case. Users may try adjusting these values with the sliders to see how it changes the optimal flow.
"""
)
return


@app.cell(hide_code=True)
def _(mo, os, re):
# model loading widget
Expand Down Expand Up @@ -69,7 +189,7 @@ def _(lines, np, num_buses, sns):
for _ix, _l in enumerate(lines):
B[_l[0]-1, _ix] = -1
B[_l[1]-1, _ix] = 1
fig_heatmap = sns.heatmap(B, cmap='seismic', center=0)
fig_heatmap = sns.heatmap(B, cmap='bwr', center=0)
return B, fig_heatmap


Expand Down Expand Up @@ -116,7 +236,7 @@ def _(B, flow_dict, generator_dict, load_dict, make_problem):


@app.cell(hide_code=True)
def _(flow_limits, gen_limits, load_limits, mo, problem):
def _(flow_limits, gen_limits, load_limits, mo, np, problem):
# solve with updated parameters
am_solving = True
problem.param_dict['g_min'].value = [_v[0] for _v in gen_limits.value]
Expand All @@ -125,8 +245,9 @@ def _(flow_limits, gen_limits, load_limits, mo, problem):
problem.param_dict['l_upper'].value = [_v[1] for _v in load_limits.value]
problem.param_dict['line capacities'].value = flow_limits.value
obj_val = problem.solve(verbose=False, solver='CLARABEL')
constrained_lines = np.where(~np.isclose(problem.constraints[2].dual_value, 0, atol=1e-2))[0]
mo.md(f'objective value: {obj_val:.2f}')
return (am_solving,)
return am_solving, constrained_lines


@app.cell
Expand All @@ -147,7 +268,8 @@ def _(generator_dict, mo, np):
stop=_max,
label=f'gen {_ix}',
value=[_gen['p_min'], _gen['p_max']],
debounce=True)
debounce=True,
show_value=True)
for _ix, _gen in generator_dict.items()], label='generator limits')
return (gen_limits,)

Expand All @@ -160,7 +282,8 @@ def _(load_dict, mo, np):
stop=_max,
label=f'load {_ix}',
value=[_load['l_min'], _load['l_upper']],
debounce=True)
debounce=True,
show_value=True,)
for _ix, _load in load_dict.items()], label='load limits')
return (load_limits,)

Expand All @@ -173,7 +296,8 @@ def _(flow_dict, mo, np):
stop=_max,
label=f'line {_ix+1}',
value=_line['f_max'],
debounce=True)
debounce=True,
show_value=True,)
for _ix, _line in flow_dict.items()], label='flow limits')
return (flow_limits,)

Expand Down Expand Up @@ -218,7 +342,7 @@ def _(


@app.cell(hide_code=True)
def _(model_data, model_ui):
def _(model_data, model_ui, np):
# construct dictionaries for problem formulation
def merge_nested_dicts(dict1, dict2):
result = {}
Expand All @@ -233,112 +357,38 @@ def merge_nested_dicts(dict1, dict2):
result[key] = merged_inner_dict

return result
generator_dict1 = {int(_g[0]): {'p_min': _g[9], 'p_max': _g[8]} for _g in model_data['gen']}
# minimum generation value occaisonally negative
generator_dict1 = {int(_g[0]): {'p_min': np.clip(_g[9], 0, np.inf), 'p_max': _g[8]*2} for _g in model_data['gen']}
generator_dict2 = {int(model_data['gen'][_ix, 0]): {'c0': _g[-1], 'c1': _g[-2], 'c2': _g[-3]} for _ix, _g in enumerate(model_data['gencost'])}
generator_dict = merge_nested_dicts(generator_dict1, generator_dict2)
load_dict = {int(_l[0]): {'l_min': _l[2]*0.5, 'l_upper': _l[2], 'cost': 250} for _l in model_data['bus']}
# power demand is occaisonally negative
# load_dict = {int(_l[0]): {'l_min': 0, 'l_upper': np.abs(_l[2]), 'cost': 250} for _l in model_data['bus']}
load_dict = {int(_l[0]): {'l_min': np.abs(_l[2]*0.5), 'l_upper': np.abs(_l[2]), 'cost': 1e4} for _l in model_data['bus']}
_fmax = 10 * np.max([np.max(_b[5:7]) for _b in model_data['branch']])
if model_ui.value != 'case14.py':
flow_dict = {int(_ix): {'f_max': _b[5], 'f_resistance': _b[2]} for _ix, _b in enumerate(model_data['branch'])}
flow_dict = {}
for _ix, _b in enumerate(model_data['branch']):
_fx = np.max(_b[5:7])
if _fx > 0:
flow_dict[_ix] = {'f_max': _fx, 'f_resistance': _b[2]}
else:
flow_dict[_ix] = {'f_max': _fmax, 'f_resistance': _b[2]}
else:
flow_dict = {int(_ix): {'f_max': 175, 'f_resistance': _b[2]} for _ix, _b in enumerate(model_data['branch'])}
line_resistance = model_data['branch'][:, 2]
return flow_dict, generator_dict, load_dict


@app.cell(hide_code=True)
def _(cp, np):
def make_problem(incidence_mat, generator_dict, load_dict, flow_dict):
def indices_not_in_list(array_length, given_indices):
all_indices = set(range(array_length))
given_indices_set = set(given_indices)
not_in_list_indices = all_indices - given_indices_set
return list(not_in_list_indices)
_B = incidence_mat
num_buses = _B.shape[0]
num_edges = _B.shape[1]
gen_ixs = np.array(list(generator_dict.keys())) - 1
nogen_ixs = indices_not_in_list(num_buses, gen_ixs)
p_flows = cp.Variable((num_edges), name='line flows')
p_gen = cp.Variable((num_buses), nonneg=True, name='p_gen')
p_load = cp.Variable((num_buses), nonneg=True, name='p_load')
l_min = cp.Parameter(
shape=num_buses,
value=[_v['l_min'] for _k, _v in load_dict.items()],
name='l_min'
)
l_upper = cp.Parameter(
shape=num_buses,
value=[_v['l_upper'] for _k, _v in load_dict.items()],
name='l_upper'
)
# l_cost = cp.Parameter(
# shape=num_buses,
# value=[_v['cost'] for _k, _v in load_dict.items()],
# name='l_cost'
# )
l_cost = [_v['cost'] for _k, _v in load_dict.items()]
g_min = cp.Parameter(
shape=len(gen_ixs),
value=np.array([_v['p_min'] for _k, _v in generator_dict.items()]),
name='g_min'
)
g_max = cp.Parameter(
shape=len(gen_ixs),
value=np.array([_v['p_max'] for _k, _v in generator_dict.items()]),
name='g_max'
)
c0 = cp.Parameter(
shape=len(gen_ixs),
value=np.array([_v['c0'] for _k, _v in generator_dict.items()]),
name='c0'
)
c1 = cp.Parameter(
shape=len(gen_ixs),
value=np.array([_v['c1'] for _k, _v in generator_dict.items()]),
name='c1'
)
# c2 = cp.Parameter(
# shape=len(gen_ixs),
# value=np.array([_v['c2'] for _k, _v in generator_dict.items()]),
# name='c2'
# )
c2 = np.array([_v['c2'] for _k, _v in generator_dict.items()])
r = cp.Parameter(
shape=num_edges,
value=np.array([_v['f_resistance'] for _k, _v in flow_dict.items()]),
name='line resistance'
)
f_max = cp.Parameter(
shape=num_edges,
value=np.array([_v['f_max'] for _k, _v in flow_dict.items()]),
name='line capacities'
)
# generator costs
cost = cp.sum(c0 + c1 @ p_gen[gen_ixs] + c2 @ cp.square(p_gen[gen_ixs]))
# load curtailment costs
cost += cp.sum(
cp.multiply(
cp.neg(p_load - l_upper),
l_cost
)
)
# line penalties
cost += cp.sum_squares(cp.multiply(r, p_flows))
constraints = [
_B @ p_flows + p_gen - p_load == 0,
p_load >= l_min,
cp.abs(p_flows) <= f_max,
p_gen[nogen_ixs] == 0,
p_gen[gen_ixs] <= g_max,
p_gen[gen_ixs] >= g_min
]
problem = cp.Problem(cp.Minimize(cost), constraints)
return problem
return (make_problem,)


@app.cell(hide_code=True)
def _(am_solving, both, gen_only, lines, load_only, problem):
def _(
am_solving,
both,
constrained_lines,
gen_only,
lines,
load_only,
problem,
):
# power flow graph view 1
am_solving
solution2 = ["flowchart LR"]
Expand All @@ -365,11 +415,28 @@ def _(am_solving, both, gen_only, lines, load_only, problem):
if len(both) > 0:
_str = ",".join([str(int(_i)) for _i in both])
solution2 += "\n class "+_str+" genload;"
if len(constrained_lines) > 0:
_str = ",".join([str(int(_i)) for _i in constrained_lines])
solution2 += "\n linkStyle "+_str+" color:red;"
return (solution2,)


@app.cell
def _(solution):
solution
return


@app.cell(hide_code=True)
def _(am_solving, both, gen_only, lines, load_only, problem):
def _(
am_solving,
both,
constrained_lines,
gen_only,
lines,
load_only,
problem,
):
# power flow graph view 2
am_solving
solution = ["flowchart LR"]
Expand All @@ -383,9 +450,9 @@ def _(am_solving, both, gen_only, lines, load_only, problem):
# else:
# solution.append(f" {_t:.0f}(({_t})) == {-_lf:.2f} ==> {_f:.0f}(({_f}))")
solution = "\n".join(solution)
solution += '''\n
classDef gen stroke-dasharray: 5 5;\n
classDef load fill:#f9f;\n
solution += '''
classDef gen stroke-dasharray: 5 5;
classDef load fill:#f9f;
classDef genload stroke-dasharray: 5 5,fill:#f9f;'''
if len(load_only) > 0:
_str = ",".join([str(int(_i)) for _i in load_only])
Expand All @@ -396,6 +463,9 @@ def _(am_solving, both, gen_only, lines, load_only, problem):
if len(both) > 0:
_str = ",".join([str(int(_i)) for _i in both])
solution += "\n class "+_str+" genload;"
if len(constrained_lines) > 0:
_str = ",".join([str(int(_i)) for _i in constrained_lines])
solution += "\n linkStyle "+_str+" color:red;"
return (solution,)


Expand All @@ -407,7 +477,6 @@ def _():
import cvxpy as cp
import matplotlib.pyplot as plt
import seaborn as sns
from model import Model
return cp, mo, np, os, plt, re, sns


Expand Down
Loading
Loading