Implementing binary constraint in PuLP Objective Function

Question

I'm trying to determine the maximum revenue that can be earned from a battery connected to the grid using linear programming. The battery can earn revenues in two markets, the energy market and the frequency market. My model is throwing an error when I include a binary constraint in the objective function (TypeError: Non-constant expressions cannot be multiplied).

My Objective function is:

• N is the time horizon of the optimisation
• is the energy price at time t
• are the allocated discharge and charge power at time t
• is the frequency price at time t
• is the allocated frequency power at time t

The battery should only be active in one market (energy or frequency) at each time period t. So needs a constraint that looks something like this:

where is a binary variable that activates activity x.

Ok, so that's what I trying to achieve. I'm struggling to create such a constraint in pulp that essentially switches off participation in one of the markets if the value in the other is higher (all other constraints being met). In my battery class, I've created decision variables for each of the power activities and also for their on/off status.

self.charge = \
pulp.LpVariable.dicts(
    "charging_power",
    ('c_t_' + str(i) for i in range(0,time_horizon)),
    lowBound=0, upBound=max_charge_power_capacity,
    cat='Continuous')

self.discharge = \
pulp.LpVariable.dicts(
    "discharging_power",
    ('d_t_' + str(i) for i in range(0,time_horizon)),
    lowBound=0, upBound=max_discharge_power_capacity,
    cat='Continuous')

self.freq = \
pulp.LpVariable.dicts(
    "freq_power",
    ('f_t_' + str(i) for i in range(0,time_horizon)),
    lowBound=0, upBound=max_freq_power_capacity,
    cat='Continuous')

self.charge_status = \
pulp.LpVariable.dicts(
    "charge_status",
    ('c_status_t_' + str(i) for i in range(0,time_horizon)),
    cat='Binary')

self.discharge_status = \
pulp.LpVariable.dicts(
    "discharge_status",
    ('d_status_t_' + str(i) for i in range(0,time_horizon)),
    cat='Binary')

self.freq_status = \
pulp.LpVariable.dicts(
    "freq_status",
    ('ds3_status_t_' + str(i) for i in range(0,time_horizon)),
    cat='Binary')

In my objective function, I included these binary variables.

    
self.model = pulp.LpProblem("Max Profit", pulp.LpMaximize)

self.model +=\
pulp.lpSum(
    [self.charge['c_t_' + str(i)]*-1*prices[i] *
     self.charge_status['c_status_t_' + str(i)] for i in range(0,self.time_horizon)]
    + [self.discharge['d_t_' + str(i)]*prices[i] * 
     self.discharge_status['d_status_t_' + str(i)] for i in range(0,self.time_horizon)]
    + [self.freq['f_t_' + str(i)]*freq_prices[i] *
     self.freq_status['freq_status_t_' + str(i)] for i in range(0,self.time_horizon)]
)

The constraint for these binary variables, I set up as follows:

 for hour_of_sim in range(1,self.time_horizon+1):     
    self.model += \
        pulp.lpSum([self.charge_status['c_status_t_' + str(i)] for i in range(0,self.time_horizon)] +\
                   [self.discharge_status['d_status_t_' + str(i)] for i in range(0,self.time_horizon)] +\
                   [self.freq_status['freq_status_t_' + str(i)] for i in range(0,self.time_horizon)]
                  ) <= 1

When I try to solve, I get a

TypeError: Non-constant expressions cannot be multiplied

on the objective function. Doesn't like my binary variables, runs if they are removed. There must be an alternative way of setting this up which is escaping me?

Answer 1

The comment is correct... you are violating "linearity" by multiplying 2 variables together. Fortunately, this is easy to linearize. You have a binary variable controlling the mode, so the key element you are looking for (google it) is a Big-M constraint, where you use the binary variable multiplied by a max value (or just something sufficiently large) to limit the other variable either to the max, or clamp it to zero.

An example below. I also re-arranged things a bit. You might find this style more readable. Two main things on style:

1. You are constantly re-creating the indices you are using which is really painful to read and error prone. Just make them and re-use them...and you don't need to get complicated with the index set values

2. You can easily double-index this model, which I think is more clear than making multiple sets of variables. You essentially have 2 sets you are working with: Time periods, and Op modes. Just make those sets and double index.

Example

# battery modes

import pulp

# some data
num_periods = 3
rate_limits = { 'energy'    : 10,
                'freq'      : 20}
price = 2  # this could be a table or double-indexed table of [t, m] or ....

# SETS
M = rate_limits.keys()   # modes.  probably others...  discharge?
T = range(num_periods)   # the time periods

TM = {(t, m) for t in T for m in M}

model = pulp.LpProblem('Batts', pulp.LpMaximize)

# VARS
model.batt = pulp.LpVariable.dicts('batt_state', indexs=TM, lowBound=0, cat='Continuous')
model.op_mode = pulp.LpVariable.dicts('op_mode', indexs=TM, cat='Binary')

# Constraints

# only one op mode in each time period...
for t in T:
    model += sum(model.op_mode[t, m] for m in M) <= 1

# Big-M constraint. limit rates for each rate, in each period.
# this does 2 things:  it is equivalent to the upper bound parameter in the var declaration
#                      It is a Big-M type of constraint which uses the binary var as a control <-- key point
for t, m in TM:
    model += model.batt[t, m] <= rate_limits[m] * model.op_mode[t, m]

# OBJ
model += sum(model.batt[t, m] * price for t, m in TM)

print(model)

#  solve...

Yields:

Batts:
MAXIMIZE
2*batt_state_(0,_'energy') + 2*batt_state_(0,_'freq') + 2*batt_state_(1,_'energy') + 2*batt_state_(1,_'freq') + 2*batt_state_(2,_'energy') + 2*batt_state_(2,_'freq') + 0
SUBJECT TO
_C1: op_mode_(0,_'energy') + op_mode_(0,_'freq') <= 1

_C2: op_mode_(1,_'energy') + op_mode_(1,_'freq') <= 1

_C3: op_mode_(2,_'energy') + op_mode_(2,_'freq') <= 1

_C4: batt_state_(2,_'freq') - 20 op_mode_(2,_'freq') <= 0

_C5: batt_state_(2,_'energy') - 10 op_mode_(2,_'energy') <= 0

_C6: batt_state_(1,_'freq') - 20 op_mode_(1,_'freq') <= 0

_C7: batt_state_(1,_'energy') - 10 op_mode_(1,_'energy') <= 0

_C8: batt_state_(0,_'freq') - 20 op_mode_(0,_'freq') <= 0

_C9: batt_state_(0,_'energy') - 10 op_mode_(0,_'energy') <= 0

VARIABLES
batt_state_(0,_'energy') Continuous
batt_state_(0,_'freq') Continuous
batt_state_(1,_'energy') Continuous
batt_state_(1,_'freq') Continuous
batt_state_(2,_'energy') Continuous
batt_state_(2,_'freq') Continuous
0 <= op_mode_(0,_'energy') <= 1 Integer
0 <= op_mode_(0,_'freq') <= 1 Integer
0 <= op_mode_(1,_'energy') <= 1 Integer
0 <= op_mode_(1,_'freq') <= 1 Integer
0 <= op_mode_(2,_'energy') <= 1 Integer
0 <= op_mode_(2,_'freq') <= 1 Integer