Pyomo Error: No value for uninitialized NumericValue object Variable

Question

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!