I have a problem with a new algorithm, here I show my code to understand it better.
import pyomo.environ
from pyomo.core import *
from pyomo.opt import SolverFactory
model = ConcreteModel()
model.JOBS = RangeSet(1,3)
model.STATIONS = RangeSet(1,5)
model.PERIODS = RangeSet(1,40)
model.hmaxt = Param(model.PERIODS, initialize={1:20, 2:20, 3:20, 4:20, 5:20, 6:20, 7:20, 8:20, 9:20, 10:20, 11:20, 12:20, 13:20, 14:20, 15:20, 16:20, 17:20, 18:20, 19:20, 20:20, 21:20, 22:20, 23:20, 24:20, 25:20, 26:20, 27:20, 28:20, 29:20, 30:20, 31:20, 32:20, 33:20, 34:20, 35:20, 36:20, 37:20, 38:20, 39:20, 40:20})
model.hmaxst = Param(model.STATIONS, initialize={1:3, 2:4, 3:4, 4:3, 5:3})
model.Y = Param(initialize=8)
model.L = Param(model.JOBS, model.STATIONS, initialize={(1,1):48, (1,2):40, (1,3):72, (1,4):56, (1,5):40, (2,1):40, (2,2):32, (2,3):88, (2,4):80, (2,5):72, (3,1):88, (3,2):64, (3,3):40, (3,4):64, (3,5):32})
model.x = Var(model.JOBS, model.STATIONS, model.PERIODS, within=Boolean)
model.y = Var(model.JOBS, model.STATIONS, model.PERIODS, within=Boolean)
model.h = Var(model.JOBS, model.STATIONS, model.PERIODS, within=PositiveIntegers)
def obj_rule(model):
return sum(t*model.y[1,5,t] for t in model.PERIODS)
model.obj = Objective(rule=obj_rule)
def rest1_rule(model, i, j):
return sum(model.x[i,j,t] for t in model.PERIODS) == 1
model.rest1 = Constraint(model.JOBS, model.STATIONS, rule=rest1_rule)
def rest2_rule(model, i, j):
return sum(model.y[i,j,t] for t in model.PERIODS) == 1
model.rest2 = Constraint(model.JOBS, model.STATIONS, rule=rest2_rule)
def rest3_rule(model, i, j):
return sum(t*model.x[i,j,t] for t in model.PERIODS) <= sum(t*model.y[i,j,t] for t in model.PERIODS)
model.rest3 = Constraint(model.JOBS, model.STATIONS, rule=rest3_rule)
@model.Constraint(model.JOBS, model.STATIONS)
def rest4_rule(model, i, j):
if j == 2 or j == 5:
return Constraint.NoConstraint
else:
return sum(t*model.y[i,j,t] for t in model.PERIODS) + 1 <= sum(t*model.x[i,j+1,t] for t in model.PERIODS)
def rest5_rule(model, i):
return sum(t*model.y[i,2,t] for t in model.PERIODS) + 1 <= sum(t*model.x[i,5,t] for t in model.PERIODS)
model.rest5 = Constraint(model.JOBS, rule=rest5_rule)
@model.Constraint(model.JOBS, model.STATIONS)
def rest6_rule(model, i, j):
if i == 1:
return Constraint.NoConstraint
else:
return sum(t*model.y[i-1,j,t] for t in model.PERIODS) + 1 <= sum(t*model.x[i,j,t] for t in model.PERIODS)
@model.Constraint(model.JOBS, model.STATIONS)
def rest14_rule(model, i, j):
if j == 2 or j == 5 or i == 3:
return Constraint.NoConstraint
else:
return sum(t*model.x[i+1,j,t] for t in model.PERIODS) >= sum(t*model.x[i,j+1,t] for t in model.PERIODS)
@model.Constraint(model.JOBS)
def rest15_rule(model, i):
if i == 3:
return Constraint.NoConstraint
else:
return sum(t*model.x[i+1,2,t] for t in model.PERIODS) >= sum(t*model.x[i,5,t] for t in model.PERIODS)
def rest12_rule(model, t):
return sum(sum(model.h[i,j,t] for i in model.JOBS) for j in model.STATIONS) <= model.hmaxt[t]
model.rest12 = Constraint(model.PERIODS, rule=rest12_rule)
def rest13_rule(model, i, j):
return sum((model.Y)*model.h[i,j,t] for t in model.PERIODS) >= model.L[i,j]
model.rest13 = Constraint(model.JOBS, model.STATIONS, rule=rest13_rule)
def rest9_rule(model, i, j, t):
return model.h[i,j,t] <= model.hmaxst[j] * (sum(model.x[i,j,k] for k in sequence(t)) - sum(model.y[i,j,k] for k in sequence(t - 1)))
model.rest9 = Constraint(model.JOBS, model.STATIONS, model.PERIODS, rule=rest9_rule)
opt = SolverFactory("glpk")
model.x[1,1,1] = 1
model.x[1,3,1] = 1
model.x[1,1,1].fixed = True
model.x[1,3,1].fixed = True
model.preprocess()
print("\nOptimal Solution found\n" + '-'*80)
results = opt.solve(model, tee=True)
model.solutions.load_from(results)
def pyomo_postprocess(options=None, instance=None, results=None):
model.x.display(), model.obj.display()
pyomo_postprocess(None, None, results)
print("\nResumen de la solución encontrada\n"+ '-'*80)
for t in model.PERIODS:
for j in model.STATIONS:
for i in model.JOBS:
if model.x[i,j,t].value == 1:
print("The job",i,"of the station",j,"starts on period:",t)
The problem comes when I try to add the last constraint rest9_rule:
def rest9_rule(model, i, j, t):
return model.h[i,j,t] <= model.hmaxst[j] * (sum(model.x[i,j,k] for k in sequence(t)) - sum(model.y[i,j,k] for k in sequence(t - 1)))
model.rest9 = Constraint(model.JOBS, model.STATIONS, model.PERIODS, rule=rest9_rule)
then, the optimizer shows the following Error:
GLPSOL: GLPK LP/MIP Solver, v4.65
Parameter(s) specified in the command line:
--write C:\Users\UX430U\tmpx9_ksgcn.glpk.raw --wglp C:\Users\UX430U\tmpxp64l_x_.glpk.glp
--cpxlp C:\Users\UX430U\tmp1yt7867r.pyomo.lp
Reading problem data from 'C:\Users\UX430U\tmp1yt7867r.pyomo.lp'...
C:\Users\UX430U\tmp1yt7867r.pyomo.lp:35244: warning: lower bound of variable 'x1' redefined
C:\Users\UX430U\tmp1yt7867r.pyomo.lp:35244: warning: upper bound of variable 'x1' redefined
731 rows, 1801 columns, 30601 non-zeros
1800 integer variables, 1200 of which are binary
36444 lines were read
Writing problem data to 'C:\Users\UX430U\tmpxp64l_x_.glpk.glp'...
34508 lines were written
GLPK Integer Optimizer, v4.65
731 rows, 1801 columns, 30601 non-zeros
1800 integer variables, 1200 of which are binary
Preprocessing...
PROBLEM HAS NO PRIMAL FEASIBLE SOLUTION
Time used: 0.0 secs
Memory used: 3.9 Mb (4069204 bytes)
Writing MIP solution to 'C:\Users\UX430U\tmpx9_ksgcn.glpk.raw'...
2541 lines were written
obj : Size=1, Index=None, Active=True
ERROR: evaluating expression: No value for uninitialized NumericValue object
y[1,5,1]
(expression: y[1,5,1] + 2*y[1,5,2] + 3*y[1,5,3] + 4*y[1,5,4] +
5*y[1,5,5] + 6*y[1,5,6] + 7*y[1,5,7] + 8*y[1,5,8] + 9*y[1,5,9] +
10*y[1,5,10] + 11*y[1,5,11] + 12*y[1,5,12] + 13*y[1,5,13] +
14*y[1,5,14] + 15*y[1,5,15] + 16*y[1,5,16] + 17*y[1,5,17] +
18*y[1,5,18] + 19*y[1,5,19] + 20*y[1,5,20] + 21*y[1,5,21] +
22*y[1,5,22] + 23*y[1,5,23] + 24*y[1,5,24] + 25*y[1,5,25] +
26*y[1,5,26] + 27*y[1,5,27] + 28*y[1,5,28] + 29*y[1,5,29] +
30*y[1,5,30] + 31*y[1,5,31] + 32*y[1,5,32] + 33*y[1,5,33] +
34*y[1,5,34] + 35*y[1,5,35] + 36*y[1,5,36] + 37*y[1,5,37] +
38*y[1,5,38] + 39*y[1,5,39] + 40*y[1,5,40])
ERROR: evaluating object as numeric value: obj
(object: <class 'pyomo.core.base.objective.SimpleObjective'>)
No value for uninitialized NumericValue object y[1,5,1]
Key : Active : Value
None : None : None
When the constraint is not added the error disappear, please any idea? I have tried many things and I dont know why of this error. Thanks in advance.
It looks like your problem is infeasible: PROBLEM HAS NO PRIMAL FEASIBLE SOLUTION.
You may want to model.pprint() or model.constraint.pprint() to verify that the constraints look the way that you expect.
As for the specific error message, because glpk returned with infeasible, no variable value was assigned to y[1,5,1], so evaluation of the objective value function during model.obj.display() caused the error that you see.