Python Pulp linear programming with dynamic constraint

I currently use Solver in excel to find somewhat optimal solution for manufacturing. Here is the current setup:

It regards manufacturing shoes on a rotary machine, that is, production is done in batches of repetitions. For example one batch would be '10x A1' (see A1 in table) which would yield 10x size 36, 20x size 37... 10x size 41.

There are some prefixed setups; A1, A2; R7... as you see in the table above.

Then there is the requested variable (or rather a list of variables) that basically says what the customer requested, quantities per size.

The objective function is to find a set of repetitions such that it matches requested quantities as closely as possible. Hence in the solver (sorry for non-English screenshot) you can see the goal is N21 (that is the sum of absolute differences per size). The variables are N2:N9 - that's the repetitions per setup and the only constraint is that N2:N9 is an integer.

How can I model this behaviour with python? My start:

from collections import namedtuple

from pulp import *


class Setup(namedtuple('IAmReallyLazy', 'name ' + ' '.join(f's{s}' for s in range(36, 47)))):
    # inits with name and sizes 's36', 's37'... 's46'
    repetitions = 0


setups = [
    Setup('A1', 1, 2, 3, 3, 2, 1, 0, 0, 0, 0, 0),
    Setup('A2', 0, 1, 2, 3, 3, 2, 1, 0, 0, 0, 0),
    Setup('R7', 0, 0, 1, 1, 1, 1, 2, 0, 0, 0, 0),
    Setup('D1', 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0),
    # and others
]

setup_names = [s.name for s in setups]

requested = {
    's36': 100,
    's37': 250,
    's38': 300,
    's39': 450,
    's40': 450,
    's41': 250,
    's42': 200,
}


def get_quantity_per_size(size: str) -> int:
    return sum([getattr(setup, size) * setup.repetitions for setup in setups])


def get_abs_diff(size: str) -> int:
    requested_size = requested.get(size, 0)
    return abs(get_quantity_per_size(size) - requested_size)


problem = LpProblem('Optimize Batches', LpMinimize)
# goal is to minimise the sum(get_abs_diff(f's{size}') for size in range(36, 47))
# variables are [setup.repetitions for setup in setups]
# constraints are all([isinstance(setup.repetitions, int) for setup in setups])

In an ideal world if there is more than one optimal solution, the one with most spread abs diff should be selected (ie. the one with the smallest highest difference). That is, if one solution has abs diff of 10 per size and 10 sizes (total 100) and other has 20 + 80 = 100, the first one is more optimal for customer.

Another constraint should be min(setup.repetitions for setup in setups if setup.repetitions > 0) > 9 basically the repetitions constraint should be:

Is integer
Is either 0 or greater than 9 - This is not possible in linear programming from what I've understood so far though.

Solution

A few things here. First, if you use abs() then the problem will be nonlinear. Instead, you should introduce new variables called, say, over_mfg and under_mfg, that represent the number of units produced above of the target and the number below the target, respectively. You can declare these like this:

over_mfg = LpVariable.dicts("over_mfg", sizes, 0, None, LpInteger)
under_mfg = LpVariable.dicts("under_mfg", sizes, 0, None, LpInteger)

I declared a list called sizes, which is used in the definitions above:

min_size = 36
max_size = 46
sizes = ['s' + str(s) for s in range(min_size, max_size+1)]

You also need variables indicating the repetitions of each setup:

repetitions = LpVariable.dicts("repetitions", setup_names, 0, None, LpInteger)

Your objective function is then declared as:

problem += lpSum([over_mfg[size] + under_mfg[size] for size in sizes])

(Note that in pulp you use lpSum rather than sum.) Now, you need constraints that say that over_mfg is the excess and under_mfg is the shortfall:

for size in sizes:
    problem += over_mfg[size] >= lpSum([repetitions[setup.name] * getattr(setup, size) for setup in setups]) - requested[size], "DefineOverMfg" + size
    problem += under_mfg[size] >= requested[size] - lpSum([repetitions[setup.name] * getattr(setup, size) for setup in setups]), "DefineUnderMfg" + size

Note also that I didn't use your get_quantity_per_size() and get_abs_diff() functions. These will also confuse pulp since it won't realize that these are simple linear functions.

Here is my complete code:

from collections import namedtuple

from pulp import *


class Setup(namedtuple('IAmReallyLazy', 'name ' + ' '.join(f's{s}' for s in range(36, 47)))):
    # inits with name and sizes 's36', 's37'... 's46'
    repetitions = 0


setups = [
    Setup('A1', 1, 2, 3, 3, 2, 1, 0, 0, 0, 0, 0),
    Setup('A2', 0, 1, 2, 3, 3, 2, 1, 0, 0, 0, 0),
    Setup('R7', 0, 0, 1, 1, 1, 1, 2, 0, 0, 0, 0),
    Setup('D1', 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0),
    # and others
]

setup_names = [s.name for s in setups]

min_size = 36
max_size = 46
sizes = ['s' + str(s) for s in range(min_size, max_size+1)]

requested = {
    's36': 100,
    's37': 250,
    's38': 300,
    's39': 450,
    's40': 450,
    's41': 250,
    's42': 200,
    's43': 0,    # I added these for completeness
    's44': 0,
    's45': 0,
    's46': 0
}
    
problem = LpProblem('Optimize Batches', LpMinimize)
# goal is to minimise the sum(get_abs_diff(f's{size}') for size in range(36, 47))
# variables are [setup.repetitions for setup in setups]
# constraints are all([isinstance(setup.repetitions, int) for setup in setups])

repetitions = LpVariable.dicts("repetitions", setup_names, 0, None, LpInteger)
over_mfg = LpVariable.dicts("over_mfg", sizes, 0, None, LpInteger)
under_mfg = LpVariable.dicts("under_mfg", sizes, 0, None, LpInteger)

problem += lpSum([over_mfg[size] + under_mfg[size] for size in sizes])

for size in sizes:
    problem += over_mfg[size] >= lpSum([repetitions[setup.name] * getattr(setup, size) for setup in setups]) - requested[size], "DefineOverMfg" + size
    problem += under_mfg[size] >= requested[size] - lpSum([repetitions[setup.name] * getattr(setup, size) for setup in setups]), "DefineUnderMfg" + size

# Solve problem
problem.solve()

# Print status
print("Status:", LpStatus[problem.status])

# Print optimal values of decision variables
for v in problem.variables():
    if v.varValue is not None and v.varValue > 0:
        print(v.name, "=", v.varValue)

Here is the output:

Status: Optimal
over_mfg_s38 = 62.0
over_mfg_s41 = 62.0
repetitions_A1 = 25.0
repetitions_A2 = 88.0
repetitions_D1 = 110.0
repetitions_R7 = 1.0
under_mfg_s36 = 75.0
under_mfg_s37 = 2.0
under_mfg_s40 = 25.0

So, we manufacture 25 repetitions of A1, 88 of A2, and 110 of D1, and 1 of R7. This gives 25 units of s36 (hence 75 units under the target of 100); 248 units of s37 (2 under target); 362 units of s38 (62 units over the target of 300); and so on.

Now, for your constraint that says you either have to produce 0 of a setup or >9, you can introduce new binary variables indicating whether each setup is produced:

is_produced = LpVariable.dicts("is_produced", setup_names, 0, 1, LpInteger)

Then add these constraints:

M = 1000
min_reps = 9
for s in setup_names:
    problem += M * is_produced[s] >= repetitions[s] # if it's produced at all, must set is_produced = 1
    problem += min_reps * (1 - is_produced[s]) + repetitions[s] >= min_reps

M is a big number; it should be bigger than the largest possible number of repetitions, but no bigger. And I defined min_reps to avoid "hard-coding" the 9 in the constraints. So, these constraints say that (1) if repetitions[s] > 0, then is_produced[s] must equal 1, and (2) either is_produced[s] = 1 or repetitions[s] > 9.

The output:

Status: Optimal
is_produced_A1 = 1.0
is_produced_A2 = 1.0
is_produced_D1 = 1.0
over_mfg_s38 = 63.0
over_mfg_s39 = 1.0
over_mfg_s41 = 63.0
repetitions_A1 = 25.0
repetitions_A2 = 88.0
repetitions_D1 = 112.0
under_mfg_s36 = 75.0
under_mfg_s40 = 24.0

Note that now we don't have any setups with >0 but <9 repetitions.

In an ideal world if there is more than one optimal solution, the one with most spread abs diff should be selected (ie. the one with the smallest highest difference).

This one is trickier, and (at least for now), I'll leave this to an ideal world, or another person's answer.

By the way: There is a Stack Exchange site for Operations Research, where we are all over questions like this one.