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:
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!
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))