mathematical-optimizationpyomomixed-integer-programmingipopt

Pyomo constraints to enforce integer variables


In the problem described below, I try to find the optimal mix of batch size and number of production runs to meet at minimum certain outcome value. While I am quite sure that for the given data constellation the problem is indeed infeasible, I want to avoid that solver returns infeasible due to "human error" ;)

Running the code below with ipopt, it gives the following feasible result:

Optimal number of batches (x): 0.9999999903937844

Optimal batch size (a): 6.363636302503406

Optimal total production: 6.3636362413729435

However, the values of variables a and x should be integer.

  1. How can I refine my constraints as I learned that Pyomo is not accepting constraints which return boolean?
  2. How can I formulate the constraint that 0.2a and 0.8a are integer? If boolean expressions were allowed, it would be something like float(0.2*a).is_integer().

The code:

model.x = Var(name="Number of batches", domain=NonNegativeIntegers, initialize=10)                    
model.a = Var(name="Batch Size", domain=NonNegativeIntegers, bounds=(5,20))

# Objective function
def total_production(model):
    return model.x * model.a
model.total_production = Objective(rule=total_production, sense=minimize)

# Constraints
# Minimum production of the two output products
def first_material_constraint_rule(model):
    return sum(0.2 * model.a * i for i in range(1, value(model.x)+1)) >= 70
model.first_material_constraint = Constraint(rule=first_material_constraint_rule)

def second_material_constraint_rule(model):
    return sum(0.8 * model.a * i for i in range(1, value(model.x)+1)) >= 90
model.second_material_constraint = Constraint(rule=second_material_constraint_rule)

# At least one production run
def min_production_rule(model):
    return model.x >= 1
model.min_production = Constraint(rule=min_production_rule)

Solution

  • I don't have pyomo so I demonstrate with differential_evolution. Maintaining that 0.2a, 0.8a and a are all integers is easy by defining the decision variable as 0.2*a, and deriving a within the objective and constraint functions. Also, don't run a sum within the constraint; that evaluates to x(x + 1)/2:

    import numpy as np
    from scipy.optimize import differential_evolution, Bounds, NonlinearConstraint
    
    
    def total_production(xa02: np.ndarray) -> float:
        x, a02 = xa02
        a = a02/0.2
        return x*a
    
    
    def first_material_constraint_rule(xa02: np.ndarray) -> float:
        x, a02 = xa02
        return a02 * x*(x + 1)/2
    
    
    def second_material_constraint_rule(xa02: np.ndarray) -> float:
        x, a02 = xa02
        a = a02/0.2
        return 0.8*a * x*(x + 1)/2
    
    
    # Variables: x, non-negative integer, number of batches;
    #         0.2a, non-negative integer, where a is the batch size.
    result = differential_evolution(
        func=total_production,
        x0=(10, 1),
        integrality=(True, True),
        bounds=Bounds(
            # 1 <= x <= large, 5 <= a <= 20
            lb=(1, 5*0.2),
            ub=(1000, 20*0.2),
        ),
        constraints=(
            NonlinearConstraint(
                fun=first_material_constraint_rule, lb=70, ub=np.inf),
            NonlinearConstraint(
                fun=second_material_constraint_rule, lb=90, ub=np.inf),
        ),
    )
    assert result.success
    x, a02 = result.x
    a = a02/0.2
    print('x =', x)
    print('a =', a)
    print('0.2a =', a02)
    print('0.8a =', a*0.8)
    print('total production =', result.fun)
    print('first material =', sum(0.2*a*i for i in range(1, int(round(x))+1)))
    print('second material =', sum(0.8*a*i for i in range(1, int(round(x))+1)))
    
    x = 12.0
    a = 5.0
    0.2a = 1.0
    0.8a = 4.0
    total production = 60.0
    first material = 78.0
    second material = 312.0