knapsack problem with non-linear constraint

Question

I have a basic knapsack problem where I need to fit the most weight possible in a bin:

import pyomo.environ as pyo


w = {'hammer': 5, 'wrench': 7, 'screwdriver': 4, 'towel': 3}

limit = 14

M = pyo.ConcreteModel()

M.ITEMS = pyo.Set(initialize=w.keys())

M.x = pyo.Var(M.ITEMS, within=pyo.Binary)

M.value = pyo.Objective(expr=sum(M.x[i] for i in M.ITEMS), sense=pyo.maximize)

M.weight = pyo.Constraint(expr=sum(w[i] * M.x[i] for i in M.ITEMS) <= limit)

solver = pyo.SolverFactory("cbc")
solver.solve(M)
M.x.display()

Which works fine:

    Key         : Lower : Value : Upper : Fixed : Stale : Domain
         hammer :     0 :   1.0 :     1 : False : False : Binary
    screwdriver :     0 :   1.0 :     1 : False : False : Binary
          towel :     0 :   1.0 :     1 : False : False : Binary
         wrench :     0 :   0.0 :     1 : False : False : Binary

My problem is that I need to include a value component and I must ensure that the average value of all selected items is above a certain fixed threshold. I implemented it like this:

v = {'hammer': 1, 'wrench': 3, 'screwdriver': 1, 'towel': 2}
w = {'hammer': 5, 'wrench': 7, 'screwdriver': 4, 'towel': 3}

limit = 14

M = pyo.ConcreteModel()

M.ITEMS = pyo.Set(initialize=w.keys())

M.x = pyo.Var(M.ITEMS, within=pyo.Binary)

M.value = pyo.Objective(expr=sum(M.x[i] for i in M.ITEMS), sense=pyo.maximize)

M.weight = pyo.Constraint(expr=sum(w[i] * M.x[i] for i in M.ITEMS) <= limit)

M.above_average = pyo.Constraint(
    expr=sum(M.x[i] * v[i] for i in M.ITEMS) / sum(M.x[i] for i in M.ITEMS) >= 2
)


solver = pyo.SolverFactory("cbc")
solver.solve(M)
M.x.display()

But that caused the following error:

ValueError: Cannot write legal LP file.  Constraint 'above_average' has a body with nonlinear terms.

I understand that the constraint that I added is non-linear, therefore the error message is correct. My question is how I could solve that with pyomo (or equivalent) or workaround this problem.

Answer 1

Hold the phone...

You can keep this linear. Just sum the selection variables and multiply by the min average requirement. No division required.

import pyomo.environ as pyo

v = {'hammer': 1, 'wrench': 3, 'screwdriver': 1, 'towel': 2}
w = {'hammer': 5, 'wrench': 7, 'screwdriver': 4, 'towel': 3}

limit = 14

M = pyo.ConcreteModel()

M.ITEMS = pyo.Set(initialize=w.keys())

M.x = pyo.Var(M.ITEMS, within=pyo.Binary)

M.value = pyo.Objective(expr=sum(M.x[i] for i in M.ITEMS), sense=pyo.maximize)

M.weight = pyo.Constraint(expr=sum(w[i] * M.x[i] for i in M.ITEMS) <= limit)

# M.above_average = pyo.Constraint(
#     expr=sum(M.x[i] * v[i] for i in M.ITEMS) / sum(M.x[i] for i in M.ITEMS) >= 2
# )

min_avg_val = 2

M.above_average = pyo.Constraint(
    expr=sum(M.x[i] * v[i] for i in M.ITEMS) >= sum(M.x[i] * min_avg_val for i in M.ITEMS) )


solver = pyo.SolverFactory("cbc")
solver.solve(M)
M.x.display()

Yields:

x : Size=4, Index=ITEMS
    Key         : Lower : Value : Upper : Fixed : Stale : Domain
         hammer :     0 :   0.0 :     1 : False : False : Binary
    screwdriver :     0 :   1.0 :     1 : False : False : Binary
          towel :     0 :   1.0 :     1 : False : False : Binary
         wrench :     0 :   1.0 :     1 : False : False : Binary

Answer 2

Here is a workaround for your nonlinear problem:

solver_manager = SolverManagerFactory('neos')
#solver = pyo.SolverFactory("BARON")
solution = solver_manager.solve(M, solver = 'couenne')
#solver.solve(M)
for i in M.ITEMS:
    print(i,':',pyo.value(M.x[i]))

You can either use a local nonlinear solver like BARON, or some solvers in NEOS server. I comment out the local solver option in the code, but you can easily use NEOS server to solve it. Please note that you need to add the following line to the beginning of your code:

from pyomo.opt.parallel import SolverManagerFactory

solution of the model is(using cuenne in NEOS):

hammer : 3.502043595574344e-08
wrench : 1.0
screwdriver : 1.0
towel : 1.0

and with BARON:

hammer 0.0
wrench 1.0
screwdriver 1.0
towel 1.0