piecewise.py

"""
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()