Pyomo constraints to enforce integer variables

Question

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)

Answer 1

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