Pyomo 错误：未初始化的 NumericValue 对象变量没有值

Question

我已经构建了一个多目标模型，并且正在使用 pyomo 的增强 epsilon 方法包来解决这个问题，该方法包可以在 Github 中找到。

该模型为我提供了目标值，但是当我尝试打印变量的某些单独值时，它给出了“未初始化的 NumericValue 对象没有值”错误。我做了一些研究，发现这可能表明我的模型不可行，但奇怪的是，我有我想要的客观值。这是我的代码，没有数据，因为我有一长串数据：

#Uncertain Demand Model
from pyomo.environ import *
import numpy as np
from pathlib import Path
from pyaugmecon import PyAugmecon
import pandas as pd

*Data entered here*

def model_sto():
     model = ConcreteModel()
    
    #Sets
    E=RangeSet(6)
    P=RangeSet(8)
    L=RangeSet(2)
    T=RangeSet(13) 
    S=RangeSet(4) 
    N=RangeSet(3)
    A=RangeSet(0,53) 
    model.E = Set(initialize=E) #Raw material
    model.P = Set(initialize=P) #Products
    model.L = Set(initialize=L) #Production line
    model.T = Set(initialize=T) #Planning periods in week
    model.S = Set(initialize=S) #Seasons in months
    model.N = Set(initialize=N) #Scenarios
    model.A = Set(initialize=A) #Age of raw material inventory

    #Production
    model.f = Param(model.N, model.S, model.T, model.P, initialize=f) #Demand forecast 
    model.error = Param(model.N, model.S, model.T, model.P, initialize=error) #Forecast error
    model.b = b #Lost sales penalty
    model.prate= Param(model.P, model.L, initialize = prate) #Production rate of product p on line l
    model.lcap=Param(model.L, initialize = lcap) #capacity of production line l in terms of time in week
    model.rcon=Param(model.E, model.P, initialize =rcon) #raw material consumption rate 
    model.pa =Param(model.P, model.L, initialize = pa)  #production assignment binary variable
    model.oee =Param(model.L, initialize = oee) #overall equipment effificiency, not relevant for this model
    model.i0 = 100 #initial inventory
    model.u0 = 100 #initial inventory

    #Inventory 
    model.hi = Param(model.E, initialize = hi)
    model.hp = Param(model.P, initialize = hp)

    #Environmental Impact - Not relevant for this case
    model.EI = Param(model.E, initialize = EI) 
    model.EP = Param(initialize = EP) 
    model.ES = Param(initialize = ES) 


    #Decision variables 
    model.q = Var(model.P, model.L, model.T, model.S, within=NonNegativeReals) #production quantity 
    model.r = Var(model.P, model.L, model.T, model.S, model.N, within=NonNegativeReals) #recourse production quantity
    model.w = Var(model.E, model.T, model.S, model.N, within=NonNegativeReals) #Raw material orders
    model.i = Var(model.P, model.T, model.S, model.N, within=NonNegativeReals) #Inventory level of final products
    model.u = Var(model.A, model.E, model.T, model.S, model.N, within=NonNegativeReals) #Inventory level of raw material
    model.y = Var(model.P, model.T, model.S, model.N, within=NonNegativeReals) #Lost sales 
    model.sales=Var(model.P, model.T, model.S, model.N, within=NonNegativeReals) #Sales


    ##Set of constraints

    model.ct1productInventory=ConstraintList() #Final product inventory constraint
    for n in model.N: 
        for s in model.S: 
            for t in model.T: 
                for p in model.P:
                    if t == model.T.first():
                        lhs = model.i[p,t,s,n]
                        rhs = model.i0 + sum(model.q[p,l,t,s] + model.r[p,l,t,s,n] for l in model.L) - model.sales[p,t,s,n]
                    if t!= model.T.first():
                        lhs = model.i[p,t,s,n]
                        rhs = model.i[p,t-1,s,n] + sum(model.q[p,l,t,s] + model.r[p,l,t,s,n] for l in model.L) - model.sales[p,t,s,n]
                    model.ct1productInventory.add (lhs == rhs) 

    model.ct2demand = ConstraintList() #Forecast + forecast error fulfilled from lost sales and sales
    for n in model.N: 
        for s in model.S:
            for t in model.T: 
                for p in model.P:
                    lhs = model.f[n,s,t,p] + model.error[n,s,t,p]*0.1
                    rhs = model.y[p,t,s,n] + model.sales[p,t,s,n] 
                    model.ct2demand.add (lhs == rhs) 

    model.ct3productionCap = ConstraintList() #Production capacity constraint
    for n in model.N: 
        for s in model.S: 
            for t in model.T: 
                for l in model.L:
                    lhs = sum((model.q[p,l,t,s] + model.r[p,l,t,s,n])/(model.prate[p,l]*7*24) for p in model.P)
                    rhs = model.lcap[l]*model.pa[p,l]
                    model.ct3productionCap.add(lhs <= rhs) 

    model.ct4rawMaterInventory = ConstraintList()
    for n in model.N: 
        for s in model.S:
            for t in model.T: 
                for e in model.E:
                    if t==model.T.first(): #Initial inventory constraint
                        lhs = sum(model.u[a,e,t,s,n] for a in model.A) 
                        rhs = model.u0 + model.w[e,t,s,n] - sum((model.q[p,l,t,s]+model.r[p,l,t,s,n])*model.rcon[e,p] for p in model.P for l in model.L)
                    if t!=model.T.first():
                        if e == model.E.first(): #Shelf-life constraint for Milk 
                            lhs = sum(model.u[a,e,t,s,n] for a in A if a <=2)
                            rhs = sum(model.u[a,e,t-1,s,n] for a in A if a <=1) + model.w[e,t,s,n] - sum((model.q[p,l,t,s]+model.r[p,l,t,s,n])*model.rcon[e,p] for p in model.P for l in model.L)
                        if e != model.E.first(): #No shelf-life constraint for all other ingredients
                            lhs = sum(model.u[a,e,t,s,n] for a in A)
                            rhs = sum(model.u[a,e,t-1,s,n] for a in A) + model.w[e,t,s,n] - sum((model.q[p,l,t,s]+model.r[p,l,t,s,n])*model.rcon[e,p] for p in model.P for l in model.L)    
                    model.ct4rawMaterInventory.add(lhs == rhs)


    objective1 = sum((sum((model.hp[p]*model.i[p,t,s,n])*0.33333 for p in model.P) + sum((model.hi[e]*model.u[a,e,t,s,n])*0.33333 for a in A for e in model.E) + sum((model.b*model.y[p,t,s,n])*0.33333 for p in model.P)) for t in model.T for s in model.S for n in model.N)
    objective2 = sum(sum((model.w[e,t,s,n]*model.EI[e])*0.33333 for e in model.E) + sum((model.q[p,l,t,s]*model.EP) +(model.r[p,l,t,s,n]*model.EP*0.33333) for p in model.P for l in model.L) + sum((model.i[p,t,s,n]*model.ES)*0.33333 for p in model.P) + sum(model.u[a,e,t,s,n]*model.EI[e] for e in model.E if e == 1 for a in A if a == 2) for t in model.T for s in model.S for n in N)

    model.obj_list=ObjectiveList()
    model.obj_list.add(expr=objective1, sense = minimize)
    model.obj_list.add(expr=objective2, sense = minimize) 

    #By default, deactive all the objectives for PyAugmecon package
    for o in range(len(model.obj_list)):
        model.obj_list[o + 1].deactivate()
    
    return model


model_type = 'model_sto'

    # AUGMECON related options
options = {
        'name': model_type,
        'grid_points': 40,
    
        }

solver_options = {}

#Solve using PyAugmecon Package
py_augmecon=PyAugmecon(model_sto(), options)
py_augmecon.solve()
# Options passed to Gurobi

print(py_augmecon.model.payoff) # this prints the payoff table
print(py_augmecon.unique_pareto_sols) # this prints the unique Pareto optimal solutions    

for n in model.N:
    print(value(model.w[e,t,s,n])) 

运行模型后，我得到：

Solved 40 models for 40 unique Pareto solutions in 295.36 seconds               
[[1.15809852e+08 5.72221560e+10]
 [2.54428683e+08 5.31380259e+05]]
[[1.68191673e+08 3.22793965e+10]
 [1.42790323e+08 4.40171657e+10]
 [2.27040732e+08 7.33663710e+09]
 [1.64729964e+08 3.37466177e+10]
 [1.88961929e+08 2.34760697e+10]
 [1.27744246e+08 5.13532714e+10]
 [1.21725815e+08 5.42877137e+10]
 [1.75115092e+08 2.93449542e+10]
 [2.09732185e+08 1.46727428e+10]
 [2.40887569e+08 1.46775252e+09]
 [2.02808766e+08 1.76071851e+10]
 [1.85500220e+08 2.49432908e+10]
 [1.36771892e+08 4.69516080e+10]
 [1.57845901e+08 3.66810600e+10]
 [1.15809852e+08 5.72221560e+10]
 [1.61268255e+08 3.52138388e+10]
 [2.54428682e+08 5.31380259e+05]
 [1.78576801e+08 2.78777331e+10]
 [1.54827184e+08 3.81482811e+10]
 [1.18734613e+08 5.57549348e+10]
 [2.33964150e+08 4.40219481e+09]
 [1.45799538e+08 4.25499445e+10]
 [2.16655604e+08 1.17383005e+10]
 [1.92423639e+08 2.20088485e+10]
 [2.23579022e+08 8.80385824e+09]
 [1.39781107e+08 4.54843868e+10]
 [1.51817969e+08 3.96155023e+10]
 [1.33762677e+08 4.84188291e+10]
 [1.71653383e+08 3.08121754e+10]
 [2.37425860e+08 2.93497367e+09]
 [1.48808754e+08 4.10827234e+10]
 [1.24735030e+08 5.28204925e+10]
 [2.06270476e+08 1.61399640e+10]
 [2.30502441e+08 5.86941595e+09]
 [1.82038511e+08 2.64105120e+10]
 [2.20117313e+08 1.02710794e+10]
 [1.95885348e+08 2.05416274e+10]
 [1.99347057e+08 1.90744062e+10]
 [1.30753461e+08 4.98860503e+10]
 [2.13193894e+08 1.32055217e+10]]
ERROR: evaluating object as numeric value: w[6,13,4,1]
        (object: <class 'pyomo.core.base.var._GeneralVarData'>)
    No value for uninitialized NumericValue object w[6,13,4,1]
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_20860/2659247016.py in <module>
   2737 
   2738 for n in model.N:
-> 2739     print(value(model.w[e,t,s,n]))

pyomo\core\expr\numvalue.pyx in pyomo.core.expr.numvalue.value()

pyomo\core\expr\numvalue.pyx in pyomo.core.expr.numvalue.value()

ValueError: No value for uninitialized NumericValue object w[6,13,4,1]

所以给出了结果，它们是合理的，因为它们给了我预期的结果，但是当我尝试打印各个变量时，我收到错误。我尝试打印其他变量并得到相同的错误。任何指导或帮助将不胜感激！

预先非常感谢！

Answer 1

这只是一种预感，但我认为您的命名空间内发生了一些奇怪的事情。 model 未在您的函数外部定义，并且您在调用此段中的函数时没有“捕获”您创建的实例：

#Solve using PyAugmecon Package
py_augmecon=PyAugmecon(model_sto(), options)
py_augmecon.solve()

因此，我不确定 model 如何当您稍后检查变量时，> 是一个已知变量。

尝试这个（或类似的）：

my_model = model_sto()   # <--note I'm catching the returned model
#Solve using PyAugmecon Package
py_augmecon=PyAugmecon(my_model, options)
py_augmecon.solve()

然后在后面的代码中，迭代 my_model 中的值而不是 model

for n in my_model.N:
    print(value(my_model.w[e,t,s,n])) 

增强...

啊哈，我现在在你的错误中看到你有一个ipykernel，可能来自笔记本。这解释了我上面的部分困惑。您有一个旧的 model 实例，它与您解决的问题无关。这种情况可能发生在您重复运行不同代码段的笔记本中。如果您重新运行程序以避免这些类型的名称冲突，请经常重置内核，我认为如果这样做，代码将会呕吐并且根本无法识别模型。或者您可以将其从笔记本中剥离出来并独立运行。