Skip to content

Instantly share code, notes, and snippets.

@Ogaday
Created March 15, 2024 17:12
Show Gist options
  • Save Ogaday/0eb613124aab82b720418145d7bbb364 to your computer and use it in GitHub Desktop.
Save Ogaday/0eb613124aab82b720418145d7bbb364 to your computer and use it in GitHub Desktop.
"""
Using piecewise transformations to overcome nonlinearity in storage asset optimisation.
"""
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pyomo.environ as pyo
from IPython.display import display
# Piecewise transformation to scale trades based on efficiency
rng = np.random.default_rng(seed=42)
model = pyo.ConcreteModel()
N = 48
# Index
model.T = pyo.Set(initialize=list(range(N)), ordered=True)
prices = rng.uniform(size=(N,))
# Params
model.energy_lb = pyo.Param(initialize=2)
model.energy_ub = pyo.Param(initialize=10)
model.initial_energy = pyo.Param(initialize=5)
model.buy_prices = pyo.Param(model.T, initialize=prices + 0.005)
model.sell_prices = pyo.Param(model.T, initialize=prices - 0.005)
# Variables
## Technically bids and offers are not necessary because we calculate revenue directly from scaled trades
model.buys = pyo.Var(model.T, initialize=[0] * N, bounds=(0, 3))
model.sells = pyo.Var(model.T, initialize=[0] * N, bounds=(0, 3)) # Min/max trades are actually scaled by efficiency
model.net_trades = pyo.Var(model.T, initialize=[0] * N, bounds=(-3, 3))
model.scaled_trades = pyo.Var(model.T, initialize=[0] * N, bounds=(-3, 3))
model.revenue = pyo.Var(model.T)
## Add existing trades by creating a bounds rule for max import and export
## OR
## Use intermediary param/constraint combo like net trades to sum existing trades and new trades.
# Piecewise scaling
bins = [-3, 3]
def scale_trades(model, t, x):
if x < 0:
return x * 0.95
else:
return x
def get_revenue(model, t, x):
if x < 0:
return x * model.buy_prices[t]
else:
return x * model.sell_prices[t]
model.scaler = pyo.Piecewise(
model.T,
model.scaled_trades, # codomain
model.net_trades, # domain
f_rule=scale_trades,
pw_pts=bins,
pw_constr_type="EQ", # Force function output to lie
# pw_repn='DCC',
)
model.revenue_scaler = pyo.Piecewise(
model.T,
model.revenue, # codomain
model.scaled_trades, # domain
f_rule=get_revenue,
pw_pts=bins,
pw_constr_type="EQ",
# pw_repn='DCC',
)
# Constraints
def net_trades(model, t):
return model.net_trades[t] == model.sells[t] - model.buys[t]
def energy_balance(model):
return sum(model.scaled_trades[t] for t in model.T) == 0
def energy_lower(model, t):
return model.initial_energy - sum(model.scaled_trades[inner_t] for inner_t in model.T if inner_t <= t) >= model.energy_lb
def energy_upper(model, t):
return model.initial_energy - sum(model.scaled_trades[inner_t] for inner_t in model.T if inner_t <= t) <= model.energy_ub
model.net_trades_constraint = pyo.Constraint(model.T, rule=net_trades)
model.energy_balance = pyo.Constraint(rule=energy_balance)
model.energy_lower = pyo.Constraint(model.T, rule=energy_lower)
model.energy_upper = pyo.Constraint(model.T, rule=energy_upper)
model.objective = pyo.Objective(
expr=lambda model: pyo.summation(model.revenue),
sense=pyo.maximize
)
solver = pyo.SolverFactory("glpk")
solver.solve(model).write()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment