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)]
          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)


# 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.

  • $\begingroup$ Can you change this into a minimum working example and provide a sample dataset? $\endgroup$ – Richard Jun 25 '20 at 0:13
  • $\begingroup$ @Richard I would normally always try to post as an MRE, however in this case it is very difficult owing to the size of the logic (>2000 lines) with no quick way to refactor it all in order to produce this - apologies. $\endgroup$ – cmp Jun 25 '20 at 7:06
  • $\begingroup$ Don't worry about the binary variable value of 1000. The solver has declared the problem infeasible, so any variable value is likely meaningless. Some solvers/modeling systems just return NaN for infeasible problems. $\endgroup$ – Mark L. Stone Jun 25 '20 at 19:41
  • $\begingroup$ Read yalmip.github.io/debugginginfeasible . Paragraph 1 is YAMLIP-specific and doesn't apply to you, but the rest does. $\endgroup$ – Mark L. Stone Jun 25 '20 at 20:08
  • 1
    $\begingroup$ Thanks @MarkL.Stone. So I have managed to work out which constraint caused the issue, although now the binary, whilst working, allows the ships to load on the first day without sufficient material at the port - which suggests its an issue with the formulation. $\endgroup$ – cmp Jun 26 '20 at 13:06

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.

  • $\begingroup$ Thank you for your help once again Wesley - I have already formulated a constraint as per your recommendation above - I have added this to my OP. This constraint I believe works fine, as the outputs behave as expected (inputs from rail are reflected in port inventory day to day). However, when I set the input inventory to zero, for a particular grade, the vessel_load_start_date[vessel, date]still outputs 1, instead of waiting until sufficient supply is available at the port in order for the vessel to load? $\endgroup$ – cmp Jul 6 '20 at 16:04
  • $\begingroup$ Is the variable for sales demand vars that corresponds to vessel_load_start_date positive or is it zero? If it is zero, then the inventory may be controlled correctly. This would represent an issue that could be resolved by using equality rather than less than or equal, but only if you load exactly the amount of inventory demanded. Another thing to check is if the vessel_load_start_date for each vessel is 1 for only one date. $\endgroup$ – Wesley Dyk Jul 6 '20 at 16:23
  • $\begingroup$ Yes, the current output of the vessel_sales_demand_vars is correct. I am currently testing with one vessel requiring two products. The vessel_load_start_date output is also correct insofar as it is 1 for the first date in the horizon, however the output should initially be zero, as the there is no product initially available for one of the grades? $\endgroup$ – cmp Jul 6 '20 at 16:31
  • $\begingroup$ Does the instance become infeasible if you remove the insufficient_port_supply term from the constraints? $\endgroup$ – Wesley Dyk Jul 6 '20 at 16:55
  • $\begingroup$ Interestingly it doesnt - I am unsure why because it should become infeasible by removing this constraint. This constraint has a very heavy coefficient to prevent the model from simply using this rather than charging demurrage though. $\endgroup$ – cmp Jul 6 '20 at 17:00

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.

  • $\begingroup$ Thank you for this @Richard. I suspect that the reason for the infeasible solution is because of how the binary variable is defined rather than solution infeasibility - I think i might have made a mistake in its definition. $\endgroup$ – cmp Jul 2 '20 at 11:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.