Simple linear programme for a reservoir system with a target level

res.py

'''
A simple model of a reservoir, which supplies water to a demand center and can
request water from another supply. The objective is to minimise the deficit in
supply to the demand center, then to minimise the difference between the
reservoir level at the end of the timestep and the target level.

variables:
  - volume demanded by reservoir
  - deficit between target volume and final volume
  - volume supplied by reservoir
  - deficit between requested volume and supplied volume
'''

import glpk

lp = glpk.LPX()

RES_MAX_VOL = 100
RES_CURRENT_VOL = 50
RES_TARGET_VOL = 70
LIMIT_RES_DEMAND = 15
COST_SUPPLY_RESERVOIR = -0.5
DEMAND = 10
LIMIT_SUPPLY = 15
COST_SUPPLY = 5
COST_DEMAND_FAIL = -1000

# variable representing the volume demanded by the reservoir
# this is somewhat simplified, as we're only indirectly modelling the supply here
lp.cols.add(1)
col_res_demand = lp.cols[-1]
col_res_demand.name = 'col_res_demand'
col_res_demand.bounds = 0, LIMIT_RES_DEMAND
lp.obj[col_res_demand.index] = COST_SUPPLY_RESERVOIR

# variable representing the deficit between the reservoir volume at the end of
# the timestep and the target volume
lp.cols.add(1)
col_volume_target_deficit = lp.cols[-1]
col_volume_target_deficit.name = 'col_volume_target_deficit'
col_volume_target_deficit.bounds = 0, None
lp.obj[col_volume_target_deficit.index] = -1

# variable representing the supply from the reservoir to the demand center
lp.cols.add(1)
col_supply = lp.cols[-1]
col_supply.name = 'col_supply'
col_supply.bounds = 0, LIMIT_SUPPLY

# constrain the change to the reservoir level
# - cannot decrease more than the current volume (i.e. go below zero)
# - cannot increase more than the difference between the maximum volume and then
#   current volume (i.e. go above 100%)
lp.rows.add(1)
row = lp.rows[-1]
row.matrix = [(col_res_demand.index, 1), (col_supply.index, -1)]
row.bounds = -RES_CURRENT_VOL, RES_MAX_VOL - RES_CURRENT_VOL

# the deficit between the target volume and the volume at the end of the timestep
# is equal to the target volume minus (the current volume plus the change in volume),
# where the change in volume is equal to the reservoir demand minus supply
lp.rows.add(1)
row = lp.rows[-1]
row.matrix = [(col_volume_target_deficit.index, 1), (col_res_demand.index, 1), (col_supply.index, -1)]
row.bounds = RES_TARGET_VOL - RES_CURRENT_VOL, None

# variable representing the deficit between the volume requested by the demand
# and the volume supplied by the reservoir
lp.cols.add(1)
col_demand_deficit = lp.cols[-1]
col_demand_deficit.name = 'col_demand_deficit'
col_demand_deficit.bounds = 0, DEMAND
lp.obj[col_demand_deficit.index] = COST_DEMAND_FAIL

# the demand deficit is equal to the demand requested minus the volume supplied
lp.rows.add(1)
row = lp.rows[-1]
row.matrix = [(col_supply.index, 1), (col_demand_deficit.index, 1)]
row.bounds = DEMAND, DEMAND

# cannot supply the demand center more than it requests
lp.rows.add(1)
row_demand = lp.rows[-1]
row_demand.bounds = 0, DEMAND

# maximize the benefit (i.e. minimize the cost)
lp.obj.maximize = True

lp.simplex()
assert(lp.status == 'opt')
print('obj val: {}'.format(lp.obj.value))
for col in lp.cols:
    print('{}: {}'.format(col.name, col.primal))

print('start volume: {}'.format(RES_CURRENT_VOL))
print('target volume: {}'.format(RES_TARGET_VOL))
print('new volume: {}'.format(RES_CURRENT_VOL + col_res_demand.primal - col_supply.primal))
print('change in volume: {}'.format(col_res_demand.primal - col_supply.primal))