pythongloballocalpyomoglpk

Pyomo glpk solver doesn't give me the optimum

I hope someone can help me. I am practicing with optimization modelling and I am solving the following LP problem using pyomo glpk:

max z = 4x1 + 3x2

Subject to:

• x1 + x2 <= 40
• 2x1 + x2 <= 60
• x1, x2 >= 0

The code I have is as follows:

``````# Defining the model
model = pyo.ConcreteModel()

# Decision variables
model.x1 = pyo.Var(within = pyo.NonNegativeReals)
x1 = model.x1
model.x2 = pyo.Var(within = pyo.NonPositiveReals)
x2 = model.x2

# Objective function
model.Obj = pyo.Objective(expr = 4*x1+3*x2, sense = pyo.maximize)

# Constraints
model.Const1 = pyo.Constraint(expr = x1+x2<=40)
model.Const2 = pyo.Constraint(expr = 2*x1+x2<=60)

# Run the solver
optm = SolverFactory('glpk')
results = optm.solve(model)

# Show the results
print(results)
print('Objective function = ', model.Obj())
print('x1 = ', x1())
print('x2 = ', x2())
``````

And the results I get are:

``````Problem:
- Name: unknown
Lower bound: 120.0
Upper bound: 120.0
Number of objectives: 1
Number of constraints: 3
Number of variables: 3
Number of nonzeros: 5
Sense: maximize
Solver:
- Status: ok
Termination condition: optimal
Statistics:
Branch and bound:
Number of bounded subproblems: 0
Number of created subproblems: 0
Error rc: 0
Time: 0.012318611145019531
Solution:
- number of solutions: 0
number of solutions displayed: 0

Objective function =  120.0
x1 =  30.0
x2 = 0.0
``````

However, the result should be:

``````Object function = 140.0
x1 = 20.0
x2 = 20.0
``````

Since I only use linear equations, I believe it is both convex and concave, not sure if local optima exist in this case?

Otherwise, can anyone tell me what I am doing wrong?

Thank you so much in advance for your help!

Solution

• You are on the right track. You have an unfortunate typo that is biting you. You declared the domain of `x2` to be non-positive where you clearly intended `pyo.NonNegativeReals`

If you are having odd behavior, always `pprint` and/or `display` your model. Errors tend to stand out pretty quickly. `pprint` shows the construction, `display` is similar, but shows the evaluation of the expressions with the values.

2 other minor nits... I would not rename your variables, just type them out. Also, I believe value(var) is the preferred way to access the values. Here is a working version with a couple edits.

``````import pyomo.environ as pyo

# Defining the model
model = pyo.ConcreteModel()

# Decision variables
model.x1 = pyo.Var(within = pyo.NonNegativeReals)
# x1 = model.x1
model.x2 = pyo.Var(within = pyo.NonNegativeReals)
# x2 = model.x2

# Objective function
model.Obj = pyo.Objective(expr = 4*model.x1+3*model.x2, sense = pyo.maximize)

# Constraints
model.Const1 = pyo.Constraint(expr = model.x1+model.x2<=40)
model.Const2 = pyo.Constraint(expr = 2*model.x1+model.x2<=60)

# Run the solver
optm = pyo.SolverFactory('glpk')
results = optm.solve(model)

model.display()

# Show the results
print(results)
print('Objective function = ', pyo.value(model.Obj))
print('x1 = ', pyo.value(model.x1))
print('x2 = ', pyo.value(model.x2))
``````