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