Pyomo: ошибка Нет значения для неинициализированного объекта NumericValue Pyomo

У меня проблема с новым алгоритмом, здесь я показываю свой код, чтобы лучше его понять.

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)

Проблема возникает, когда я пытаюсь добавить последнее ограничение 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)

тогда оптимизатор показывает следующую ошибку:

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

Когда ограничение не добавлено, ошибка исчезает, пожалуйста, есть идеи? Я много чего пробовал и не знаю, почему возникла эта ошибка. Заранее спасибо.