tried different formulation but the code won't run

Question

I want to implement CLR model in pyomo. It seems the runtime is endless because when I run it nothing will happen and just a star comes in and won't end running. I've tried with different formulation of const5 but it doesn't work. even I've tried solving with only const1 and const2 and it worked. and solving by only const5 it works but it doesn't work when I add all together. I would appreciate if someone can help me. this is the mathematical model:

I have the values of X,Y and in M in excel file and they are correct. this is my code:

import pyomo.environ as pyo
import pandas as pd


model = pyo.ConcreteModel()


model.Iset = pyo.Set(initialize=range(1, 44))
model.Bset = pyo.Set(initialize=[1,2])

X = pd.read_excel("D:\\Project 2\\DebrisFlow.xlsx", sheet_name='X', header=0, index_col=0)
X=X.values
X={(i):X[i] for i in model.Iset}

Y = pd.read_excel("D:\\Project 2\\DebrisFlow.xlsx", sheet_name='Y', header=0, index_col=0)
Y=Y.values
Y={(i):Y[i] for i in model.Iset}

M = pd.read_excel("D:\\Project 2\\DebrisFlow.xlsx", sheet_name='w', header=0, index_col=0)
M=M.values
M={(i):M[i] for i in model.Iset}

model.X=pyo.Param(model.Iset,initialize=X)
model.Y=pyo.Param(model.Iset,initialize=Y)
model.M=pyo.Param(model.Iset,initialize=M)

model.e = pyo.Var(model.Iset, domain=pyo.NonNegativeReals)
model.c = pyo.Var(model.Bset, domain=pyo.Reals)
model.d = pyo.Var(model.Bset ,domain=pyo.Reals)
model.delta = pyo.Var(model.Iset, model.Bset, domain=pyo.Binary)

def obj(model):
    return sum(model.e[i] for i in model.Iset)

model.obj = pyo.Objective(rule=obj,sense=pyo.minimize)


def const1_rule(model,i,b):
    return model.Y[i]-(model.c[b]*model.X[i])-model.d[b]<=model.e[i]+model.M[i]-(model.M[i]*model.delta[i,b])

model.const1 = pyo.Constraint(model.Iset, model.Bset, rule=const1_rule)


def const2_rule(model,i,b):
    return (model.c[b]*model.X[i])+model.d[b]-model.Y[i]<=model.e[i]+model.M[i]-(model.M[i]*model.delta[i,b])

model.const2 = pyo.Constraint(model.Iset, model.Bset, rule=const2_rule)


def const3_rule(model,b):
    return pyo.inequality(0.0537, model.c[b], 0.1245)

model.const3 = pyo.Constraint(model.Bset,rule=const3_rule)

def const4_rule(model,b):
    return pyo.inequality(-26.0627, model.d[b], 10.3121)

model.const4 = pyo.Constraint(model.Bset, rule=const4_rule)

def const5_rule(model,i):
    return sum(model.delta[i, b] for b in model.Bset) == 1

model.const5 = pyo.Constraint(model.Iset, rule=const5_rule)




opt = pyo.SolverFactory('glpk')
result=opt.solve(model,'glpk')

print (str(result.solver))

print(pyo.value(model.obj))

Answer 1

I think your model is fine. It has elastic constraints and I don't see anything wrong with it. I think that for whatever reason glpk is just struggling with it.

When I patched in the "small data" below, your model solves almost instantly with glpk. If I go to |I| ~20 it bogs down, and 44 seems like an eternity. So you are probably just waiting on the solver, which is odd because this is a MILP, but not a large one.

Anyhow, if you have another solver, try it. When I use cbc the solve is about 1 second for |I| ~44

You can toggle between the small data and the larger stuff by commenting out sections below.

Code:

import pyomo.environ as pyo
from random import random
# import pandas as pd

i_size = 45

model = pyo.ConcreteModel()


model.Iset = pyo.Set(initialize=range(1, i_size))
model.Bset = pyo.Set(initialize=[1,2])

# X = pd.read_excel("D:\\Project 2\\DebrisFlow.xlsx", sheet_name='X', header=0, index_col=0)
# X=X.values
# X={(i):X[i] for i in model.Iset}

# Y = pd.read_excel("D:\\Project 2\\DebrisFlow.xlsx", sheet_name='Y', header=0, index_col=0)
# Y=Y.values
# Y={(i):Y[i] for i in model.Iset}

# M = pd.read_excel("D:\\Project 2\\DebrisFlow.xlsx", sheet_name='w', header=0, index_col=0)
# M=M.values
# M={(i):M[i] for i in model.Iset}

# small test...  works / optimal
Y = { 1: 13, 2: 10, 3: 5}
X = { 1: 10, 2: 15, 3: 2}
M = { 1: 99, 2: 99, 3: 99}

# big test:
Y = {k: random()*100 for k in range(1, i_size)}
X = {k: random()*100 for k in range(1, i_size)}
M = {k: 200 for k in range(1, i_size)}


model.X=pyo.Param(model.Iset,initialize=X)
model.Y=pyo.Param(model.Iset,initialize=Y)
model.M=pyo.Param(model.Iset,initialize=M)

model.e = pyo.Var(model.Iset, domain=pyo.NonNegativeReals)
model.c = pyo.Var(model.Bset, domain=pyo.Reals)
model.d = pyo.Var(model.Bset ,domain=pyo.Reals)
model.delta = pyo.Var(model.Iset, model.Bset, domain=pyo.Binary)

def obj(model):
    return sum(model.e[i] for i in model.Iset)

model.obj = pyo.Objective(rule=obj,sense=pyo.minimize)


def const1_rule(model,i,b):
    return model.Y[i]-(model.c[b]*model.X[i])-model.d[b]<=model.e[i]+model.M[i]-(model.M[i]*model.delta[i,b])

model.const1 = pyo.Constraint(model.Iset, model.Bset, rule=const1_rule)


def const2_rule(model,i,b):
    return (model.c[b]*model.X[i])+model.d[b]-model.Y[i]<=model.e[i]+model.M[i]-(model.M[i]*model.delta[i,b])

model.const2 = pyo.Constraint(model.Iset, model.Bset, rule=const2_rule)


def const3_rule(model,b):
    return pyo.inequality(0.0537, model.c[b], 0.1245)

model.const3 = pyo.Constraint(model.Bset,rule=const3_rule)

def const4_rule(model,b):
    return pyo.inequality(-26.0627, model.d[b], 10.3121)

model.const4 = pyo.Constraint(model.Bset, rule=const4_rule)

def const5_rule(model,i):
    return sum(model.delta[i, b] for b in model.Bset) == 1

model.const5 = pyo.Constraint(model.Iset, rule=const5_rule)




opt = pyo.SolverFactory('glpk')
result=opt.solve(model) #,'glpk')

print (str(result.solver))

print(pyo.value(model.obj))

model.delta.display()