I have created a model that I believe addresses your objective of minimizing demurrage charges. It utilizes indicator variables for each day and vessel showing when each vessel is finally loaded. The factors for the variables in the objective function are to be calculated based on the demurrage formula you provided.
The first thing to do with this model is to ensure that you have the port inventory working correctly, then add the demurrage calculation for the objective. The formulation I have provided does not deal with the two day loading period. You'll want to fill in details on how inventory increases at the port.
I hope this is helpful.
vessels = ['a', 'b']
sales_days = list(range(20))
sales = {1: {"a": 10}, 3: {"b": 10}}
vessel_qty = {"a": 10, "b": 10}
per_diem_charge = 100
laycan_days = 10
initial_port_inventory = 2
daily_inventory_increase = 1
model = pulp.LpProblem("PortDemurrage")
day_port_inventory = pulp.LpVariable.dict("port_inventory", sales_days,
lowBound=0, cat="Integer")
# set initial inventory
model += day_port_inventory[0] == initial_port_inventory
vessel_start_dates = {vessel: day for day in sales for vessel in sales[day]}
# denote the vessel sales by day the vessel is finally loaded
vessel_sales = pulp.LpVariable.dicts("vessel_sales", [f"{vessel}_{day}" for
vessel in vessels for day in sales_days[vessel_start_dates[vessel]:]],
lowBound=0, cat="Integer")
for day in sales_days[1:]:
day_sales = 0
for vessel in vessels:
if vessel_start_dates[vessel] <= day:
day_sales += vessel_sales[f"{vessel}_{day}"]
model += day_port_inventory[day] == day_port_inventory[day - 1] +
daily_inventory_increase - day_sales
for vessel in vessels:
vessel_start_date = vessel_start_dates[vessel]
possible_days = sales_days[vessel_start_date:]
model += sum(vessel_sales[f"{vessel}_{day}"] for day in possible_days) ==
vessel_qty[vessel]
obj = 0
final_day_vars = []
for vessel in vessels:
demurrage_day_start = vessel_start_dates[vessel]+laycan_days+1
vessel_start_date = vessel_start_dates[vessel]
for day in sales_days[vessel_start_date:]:
vessel_day_var = pulp.LpVariable(f"{vessel}_sales_day_{day}",
lowBound=0, upBound=1, cat="Integer")
final_day_vars.append(vessel_day_var)
if day >= demurrage_day_start:
obj += vessel_day_var * (per_diem_charge * (day -
vessel_start_dates[vessel] - laycan_days))
model += vessel_sales[f"{vessel}_{day}"] <= vessel_day_var *
vessel_qty[vessel]
model += sum(final_day_vars) == len(vessels)
model += obj, "minimize demurrage charge"
Solving this model produces a solution where vessel a is finally loaded on day 8 and vessel b is finally loaded on day 18. The objective value of this solution is 500, since vessel b needs to wait 15 days after arrival.
