Output of binary variable greater than one

Question

I am working on a shipping demurrage problem that uses a binary variable to denote the date a specific vessel can be loaded (I have been kindly helped by Wesley on OR before with this).

I am confident that the rest of the model works fine, however when using the binary to determine the load dates, the model is unable to find a solution i.e. is 'infeasible' and the output of the variable is either 1.0 or 10000.0 which is clearly incorrect (and non-binary).

Question: Do I need to linearize the binary variable, vessel_load_start_date[vessel, date] in some way?

Code & Descriptions below:

• port_inventory_vars = Variable: the inventory on the specified date of each product grade.
• vessel_grade_demand_tonnes = Constant: the required amount, in tonnes of each grade required by each vessel.
• vessel_sales_demand_vars[(vessel, grade, date)] = Variable. The date a vessels demand requirements are fully satisfied.
• vessel_load_start_date[vessel, date] = Binary: The date indicating when a vessel can be loaded. NOTE a vessel can only load if the total amount it requires is available in the port inventory, port_inventory_vars[date, grade].
• demurrage_rates[vessel, date] = Constant: The daily demurrage rates per vessel.
• demurrage_charge_vars[vessel, date] = Decision variable: The demurrage charged.

    # PORT STOCKPILE: Port Stockpile Inventory
    for date in dates:
      current_date = PLAN_START_DATE
      date_t_minus_one = datetime.datetime.strptime(date, '%Y-%m-%d') \
          - datetime.timedelta(days=1)
      date_t_minus_one = date_t_minus_one.strftime('%F')
      for grade in grades:
        # Filter plants
        _plants_combo = [
            plant for plant in plants 
            if (plant, grade) in plant_combinations]
        # Get vessel demands for requisite date
        _vessel_demands_combination = [
            (vessel, date) for vessel in vessels for date in dates
            if (vessel, date) in vessel_load_start_date
        ]
        if date == current_date:
          # Current Inv == current inventory + train in - sales demand
          model += port_stockpile_current[grade] \
              + pulp.lpSum(
                    train_consignment_variables[(date, plant, grade)] 
                    for plant in _plants_combo) \
              - pulp.lpSum(
                    vessel_sales_demand_vars[(vessel, grade, date)]
                    for vessel, date in _vessel_demands_combination) \
              + insufficient_port_supply[(date, grade)] \
              == port_inventory_vars[(date, grade)]
        else:
          model += port_inventory_vars[(f'{date_t_minus_one}', grade)] \
              + pulp.lpSum(
                    train_consignment_variables[(date, plant, grade)] 
                    for plant in _plants_combo) \
              - pulp.lpSum(
                    vessel_sales_demand_vars[(vessel, grade, date)]
                    for vessel, date in _vessel_demands_combination) \
              + insufficient_port_supply[(date, grade)] \
              == port_inventory_vars[(date, grade)]

    # Port stockpile total inventory tonnage must be <= 2.1M tonnes
    for date, grade in port_inventory_vars:
      model += pulp.lpSum(port_inventory_vars[(date, grade)]) <= 2100000
    
    # Control vessel loading
    for grade in grades:
      for vessel, date in vessel_load_start_date:
        model += vessel_sales_demand_vars[(vessel, grade, date)] - vessel_grade_demand_tonnes[vessel, grade] * vessel_load_start_date[vessel, date] <= 0
        model += vessel_sales_demand_vars[(vessel, grade, date)] <= vessel_load_start_date[vessel, date] * vessel_grade_demand_tonnes[vessel, grade] 
    
    # Vessel sales requirements must be satisfied by sales vars
    for vessel, grade in vessel_grade_requirements:
      for vessel, date in vessel_load_start_date:
        _dates = [
            tup[1] for tup in vessel_load_start_date
            if tup[0] == vessel 
        ] 
      model += pulp.lpSum(vessel_sales_demand_vars[vessel, grade, date] for date in _dates) == vessel_grade_demand_tonnes[vessel, grade] 
     
    # Demurrage charges per vessel
    for vessel, date in vessel_load_start_date:
      model += vessel_load_start_date[vessel, date] * demurrage_rates[vessel, date] == demurrage_charge_vars[vessel, date]

Current Model Outputs

    # Vessel load start date vars
>>> for vessel, date in vessel_load_start_date:
      print(vessel, date, ':', vessel_load_start_date[vessel, date].varValue

CEYLON BREEZE 2020-05-28 : 1.0
CEYLON BREEZE 2020-05-29 : 0.0
CEYLON BREEZE 2020-05-30 : 0.0

    # Demurrage Vars
>>> for vessel in demurrage_charge_vars:
      print(vessel, ':', demurrage_charge_vars[vessel].varValue)

CEYLON BREEZE : 0.0

# vessel sales demand vars
>>> for vessel, grade, date in vessel_sales_demand_vars:
      print(vessel, grade, date,':', vessel_sales_demand_vars[vessel, grade, date].varValue)

CEYLON BREEZE ZBL 2020-05-28 : 10000.0
CEYLON BREEZE ZBL 2020-05-29 : 0.0
CEYLON BREEZE ZBL 2020-05-30 : 0.0
CEYLON BREEZE MFA 2020-05-28 : 0.0
CEYLON BREEZE MFA 2020-05-29 : 0.0
CEYLON BREEZE MFA 2020-05-30 : 0.0
CEYLON BREEZE PRE 2020-05-28 : 0.0
CEYLON BREEZE PRE 2020-05-29 : 0.0
CEYLON BREEZE PRE 2020-05-30 : 0.0
CEYLON BREEZE AAE 2020-05-28 : 0.0
CEYLON BREEZE AAE 2020-05-29 : 0.0
CEYLON BREEZE AAE 2020-05-30 : 0.0
CEYLON BREEZE ACC 2020-05-28 : 10000.0
CEYLON BREEZE ACC 2020-05-29 : 0.0
CEYLON BREEZE ACC 2020-05-30 : 0.0

>>> for (date, grade) in port_inventory_vars:
      print(date, grade, ':', port_inventory_vars[(date, grade)].varValue)
2020-05-28 ZBL : 215200.0
2020-05-28 MFA : 216800.0
2020-05-28 PRE : 222000.0
2020-05-28 AAE : 200000.0
2020-05-28 ACC : 10000.0
2020-05-29 ZBL : 205200.0
2020-05-29 MFA : 216800.0
2020-05-29 PRE : 306000.0
2020-05-29 AAE : 200000.0
2020-05-29 ACC : 0.0
2020-05-30 ZBL : 195200.0
2020-05-30 MFA : 216800.0
2020-05-30 PRE : 306000.0
2020-05-30 AAE : 200000.0
2020-05-30 ACC : 32000.0

Any help gratefully received as I have been scratching my head on why the solution is infeasible. I am using dummy data and over a very short timeframe to help troubleshoot this issue.

Answer 1

The latest issue has to do with not having an inventory control constraint. You need to have a constraint like:

port_inventory_vars[(date, grade)] == port_inventory_vars[(date-1, grade)] - pulp.LpSum(vessel_sales_demand_vars[(vessel, grade, date)] for vessel in vessels) + ...

The ellipsis indicates where you can put additional terms that describe how inventory increases at the port. As long as the variables used in this constraint are non-negative, this should restrict the loading of goods until there is sufficient inventory at the port.

Answer 2

When dealing with infeasibility, I like to do two things:

a) Create the Irreducible Infeasible Subset (IIS). I don't think PuLP directly allows you to create that, however you could export your model instance and then use a (commercial) solver (e.g. Gurobi) to do so (see here for docs). This will allow you to narrow down where the infeasibility lies.

b) Add slack variables to your constraints: slack variables are "dummy" variables that you add to a constraint to ensure that it can be feasible. Consider this inequality constraint:

$$x + y \leq 5$$

with $x,y\geq 3$. This would of course be infeasible. Now however you can add a slack variable $s$ such that:

$$x + y \leq 5 + s$$

This will make the problem feasible for $s\geq 1$. Note that slack variables typically have very high coefficients in the objective function to force them to be as small as possible. This is actually also implemented in Gurobi (see here), and you may find this article as well as this page about it also helpful.