I have built a multiobjective model and I am solving this with an augmented epsilon method package for pyomo, which is available in Github.
The model gives me the objective values, but when I try to print out some individual values of variables, it gives me "No value for uninitialized NumericValue object" Error. I have done some research and saw that this could be a sign that my model is infeasible, but weirdly enough, I have the objective values that I want. here is my code without the data as I have a long list of data:
#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]))
And after running the model, I get :
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]
So the results are given and they are reasonable as they give me an expected result but when I try to print individual variables, I get an error. I tried printing other variables and got the same error. Any guidance or help would be highly appreciated!
Thanks a lot in advance!
This is just a hunch, but I think you have something odd going on within your namespace. model is not defined outside of your function and you aren't "catching" the instance you create when calling the function in this segment:
#Solve using PyAugmecon Package
py_augmecon=PyAugmecon(model_sto(), options)
py_augmecon.solve()
So, I'm not sure how model is a known variable when you inspect the variables later.
Try this (or similar):
my_model = model_sto() # <--note I'm catching the returned model
#Solve using PyAugmecon Package
py_augmecon=PyAugmecon(my_model, options)
py_augmecon.solve()
and then in your latter code, iterate over the values in my_model instead of model
for n in my_model.N:
print(value(my_model.w[e,t,s,n]))
Ahhh, I see now in your error that you have an ipykernel going, likely from a notebook. That explains part of my confusion above. You have an old instance of model floating around that is not connected to what you solved. This can happen in a notebook where you are repeatedly running different segments of code. Reset the kernel frequently if you are re-running things to avoid these types of name collisions, I think if you do, the code will barf and not recognize the model at all. Or you could just peel this out of the notebook and run it independently.