linear model snippet

model.py

# the LP model
import mip 

mdl = mip.Model(sense=mip.MINIMIZE, solver_name=mip.CBC)

# sets
days = list(range(len(demand[1])))
stores = list(demand_parameters.keys())

# variables
e = {} # the shiped units the day t to store i 
u = {} # the the number of unfulfilled orders the day t to store i 
s = {} # inventory of units the day t in store i
b_u = {} # binary variable to codificate unfulfilled var
M = 1e7 # big M

for i in stores:
    for t in days:
        e[(i,t)] = mdl.add_var(name='e_{}_{}'.format(i,t), var_type=mip.CONTINUOUS, lb=0)
        u[(i,t)] = mdl.add_var(name='u_{}_{}'.format(i,t), var_type=mip.CONTINUOUS, lb=0)
        s[(i,t)] = mdl.add_var(name='s_{}_{}'.format(i,t), var_type=mip.CONTINUOUS, lb=0)
        b_u[(i,t)] = mdl.add_var(name='b_u_{}_{}'.format(i,t), var_type=mip.BINARY)
        
# initial inventory 
for i in stores:
    mdl.add_constr(s[(i,0)] == 20, name='initial_stock_s{}'.format(i))
        
# IO ecuation definition
for t in days[1:]:
    for i in stores:
        mdl.add_constr(s[(i,t)] == s[(i,t-1)] + e[(i,t-1)] - float(demand[i][t-1]) + u[(i,t-1)], 
                       name='io_s{}_t{}'.format(i,t))

for t in days:
    for i in stores:        
        # codification unfulfilled orders     
        mdl.add_constr(u[(i,t)] >= - s[(i,t)] - e[(i,t)] + float(demand[i][t]), 
                       name='cod_u_lb_s{}_t{}'.format(i,t))
        
        mdl.add_constr(u[(i,t)] <= - s[(i,t)] - e[(i,t)] + float(demand[i][t]) + M*b_u[(i,t)], 
                       name='cod_u_ub_s{}_t{}'.format(i,t))
        
        mdl.add_constr(u[(i,t)] <= M*(1-b_u[(i,t)]), 
                       name='cod_u_ub_b_s{}_t{}'.format(i,t))                                                                                            
        #print(crt)
# max shipments per day
for t in days:
    mdl.add_constr(mip.xsum(e[(i,t)] for i in stores)<= 120, 
                   name='shipment_cap_t{}'.format(t))

# lineal objective
mdl.objective = mip.xsum(u[(i,t)] for i in stores for t in days)
    
status = mdl.optimize(max_seconds=300)
print(status)