[Python] Why is the quadratic term recognized as non-linear, non-quadratic?

I am creating a Pyomo model which has a quadratic term in the objective:

import pandas as pd
import pyomo.environ as pyo
gen_df = pd.DataFrame({'generator': {0: 1, 1: 2, 2: 3},
'a1': {0: 5, 1: 6, 2: 7},
'a2': {0: 0.1, 1: 0.2, 2: 0.15},
'pmin': {0: 200, 1: 180, 2: 10},
'pmax': {0: 600, 1: 300, 2: 200}})
def build_model():
# Create a model
model = pyo.ConcreteModel()

# Decision variables
model.I = pyo.Set(initialize=[1,2,3])
model.p = pyo.Var(model.I, domain=pyo.NonNegativeReals) # Continuous variable: how much of i to produce

# Objective Function: Min costs
def gen_cost(model):
cost = 0
for i in model.I:
a1i = generator_df.loc[generator_df['generator']==i,['a1']].values[0]
a2i = generator_df.loc[generator_df['generator']==i,['a2']].values[0]
cost += a1i * model.p
cost += a2i * model.p* model.p
return cost

def objective(model):
return (
gen_cost(model)
)
model.obj = pybjective(rule=objective, sense=pyo.minimize)

# Constraints
def min_constraint(model, i):
pmin = generator_df.loc[generator_df['generator'] == i, 'pmin'].values[0]
return model.p >= pmin
def max_constraint(model, i):
pmax = generator_df.loc[generator_df['generator'] == i, 'pmax'].values[0]
return model.p <= pmax
model.min_constraint = pyo.Constraint(model.I,rule = min_constraint)
model.max_constraint = pyo.Constraint(model.I,rule = max_constraint)

return model
# Solve
solvername='glpk'
model = build_model()
solver=pyo.SolverFactory(solvername)
print("Poly Degree: " + str(model.obj.expr.polynomial_degree()))
result = solver.solve(model, tee = True)


But get this runtime error:

RuntimeError: Cannot write legal LP file. Objective 'obj' has nonlinear terms that are not quadratic.


While I get 2 as a result of model.obj.expr.polynomial_degree())
​

Is this an issue with the function or with the solver?

Python version 3.11.13 Pandas version 2.3.0 Pyomo version 6.5.0 glpk built with conda

I tried removing the cost += a2i * model.p* model.p term and the solver worked. I then tried removing the cost += a1i * model.p term and the solver gave the same error

Continue reading...