Search code examples
pythonoptimizationpyomo

Pyomo: How to include a penalty in the objective function


I'm trying to minimize the cost of manufacturing a product with two machines. The cost of machine A is $30/product and cost of machine B is $40/product.

There are two constraints:

  • we must cover a demand of 50 products per month (x+y >= 50)
  • the cheap machine (A) can only manufacture 40 products per month (x<=40)

So I created the following Pyomo code:

from pyomo.environ import *
model = ConcreteModel()
model.x = Var(domain=NonNegativeReals)
model.y = Var(domain=NonNegativeReals)

def production_cost(m):
    return 30*m.x + 40*m.y

# Objective
model.mycost = Objective(expr = production_cost, sense=minimize)

# Constraints
model.demand = Constraint(expr = model.x + model.y >= 50)
model.maxA = Constraint(expr = model.x <= 40)

# Let's solve it
results = SolverFactory('glpk').solve(model)

# Display the solution
print('Cost=', model.mycost())
print('x=', model.x())
print('y=', model.y())

It works ok, with the obvious solution x=40;y=10 (Cost = 1600)

However, if we start to use the machine B, there will be a fixed penalty of $300 over the cost.

I tried with

def production_cost(m):
  if (m.y > 0):
    return 30*m.x + 40*m.y + 300
  else:
    return 30*m.x + 40*m.y

But I get the following error message

Rule failed when generating expression for Objective mycost with index
    None: PyomoException: Cannot convert non-constant Pyomo expression (0  <
    y) to bool. This error is usually caused by using a Var, unit, or mutable
    Param in a Boolean context such as an "if" statement, or when checking
    container membership or equality. For example,
        >>> m.x = Var() >>> if m.x >= 1: ...     pass
    and
        >>> m.y = Var() >>> if m.y in [m.x, m.y]: ...     pass
    would both cause this exception.

I do not how to implement the condition to include the penalty into the objective function through the Pyomo code.


Solution

  • Since m.y is a Var, you cannot use the if statement with it. You can always use a binary variable using the Big M approach as Airsquid said it. This approach is usually not recommended, since it turns the problem from LP into a MILP, but it is effective.

    You just need to create a new Binary Var:

    model.bin_y = Var(domain=Binary)
    

    Then constraint model.y to be zero if model.bin_y is zero, or else, be any value between its bounds. I use a bound of 100 here, but you can even use the demand:

    model.bin_y_cons = Constraint(expr= model.y <= model.bin_y*100)   
    

    then, in your objective just apply the new fixed value of 300:

    def production_cost(m):
        return 30*m.x + 40*m.y + 300*model.bin_y 
    
    model.mycost = Objective(rule=production_cost, sense=minimize)