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