gurobipy: Optimal way of writing piecewise linear constraints (SOS2 constraints) using python API

Viewed 171

I formulated simple MILP problem via gurobipy (energy storage optimisation). I wanted to introduce non-linear relationship between maximal and minimal flow allowed to go into storage and inventory (how much of a medium is already stored in the storage). I did it via following ugly-looking for loop (inventory and flow are 1D array variables obtained via gurobipy.Model.addMVar() method, and m is my gurobipy.Model class):

# PWL - injection rate constrain
ir_x = [0.0, 8000.0, 12500.0]
ir_y = [5000.0, 5000.0, 2500.0]
for a,b in zip(inventory.tolist(), injection_rate.tolist()):
    m.addGenConstrPWL(a, b, ir_x, ir_y)
m.addConstr(flow <= injection_rate)

Is there some clever way to do this without actually using for loop?

Whole working example:

import gurobipy as gp
import numpy as np

LEN = 10 # months
MAXIMAL_CAPACITY = 50000.0 # [MWh]
MINIMAL_CAPACITY = 0.0 # [MWh]
INITIAL_STORAGE = 0.0 # [MWh]
END_STORAGE = 0.0 # [MWh]

# create model
m = gp.Model(name="Simple Energy Storage")

# variables
flow = m.addMVar(LEN, lb=float("-inf"), name="flow")
inventory = m.addMVar(LEN, lb=float("-inf"), name="inventory")
injection_rate = m.addMVar(LEN, lb=float("-inf"), name="injection rate")
withdrawal_rate = m.addMVar(LEN, lb=float("-inf"), name="withdrawal rate")

# create inventory condition (cumulative sum of flow)
mat = np.tril(np.ones(shape=(LEN, LEN)))
m.addConstr(mat @ flow == inventory)

# inventory constrains
m.addConstr(inventory <= MAXIMAL_CAPACITY)
m.addConstr(inventory >= MINIMAL_CAPACITY)
m.addConstr(inventory[0] == INITIAL_STORAGE)
m.addConstr(inventory[-1] == END_STORAGE)

# PWL - injection rate constrain
ir_x = [0.0, 8000.0, 12500.0]
ir_y = [5000.0, 5000.0, 2500.0]
for a,b in zip(inventory.tolist(), injection_rate.tolist()):
    m.addGenConstrPWL(a, b, ir_x, ir_y)
m.addConstr(flow <= injection_rate)

# PWL - withdrawal rate constrain
wr_x = [0.0, 6000.0, 12500.0]
wr_y = [3000.0, 6500.0, 6500.0]
for a,b in zip(inventory.tolist(), withdrawal_rate.tolist()):
    m.addGenConstrPWL(a, b, wr_x, wr_y)
m.addConstr(flow >= -withdrawal_rate)

# objective function
m.setObjective(inventory[5], gp.GRB.MAXIMIZE)

m.optimize()
0 Answers
Related