I set the problem up below for solving using scipy.optimize.linprog. I get the required result for the overlapping case, but only an approximate solution for the non-overlapping case.

• Test for possible scenarios, for ex. a machine becomes unavailable. It's flow has to be redistributed across all other remaining machines, and has to meet all of the above goals too.

In the walk-through below I give all machines the same capacity. You could easily zero one out to mimic a scenario where a machine was unavailable. Similarly, you could set it to a fraction to mimic partial capacity, or make it a function of other factors.

For a more statistical exploration of scenarios, you could make machine capacity probabilistic and report the statistics from multiple trials.

---

Solution walk-through

Supplied data

I use the overlapping case from your examples.

The supplied data comprises:

• Products and machines

products = ['Pa', 'Pb', 'Pc', 'Pd'] machines = ['M1', 'M2', 'M3', 'M4', 'M5', 'M6'] num_products = len(products) num_machines = len(machines) 

• The quantity of each product - one third of the full load in this case.

# Pa Pb Pc Pd product_totals = np.array([6720, 2880, 3840, 5760]) // 3 

• A matrix of machine-product compatibilities

machine_product_compatibility = np.array([ #Pa Pb Pc Pd [1, 0, 1, 0], # M1: Pa, Pc [0, 1, 0, 1], # M2: Pb, Pd [0, 1, 1, 0], # M3: Pb, Pc [1, 0, 0, 1], # M4: Pa, Pd [1, 1, 0, 0], # M5: Pa, Pb [0, 0, 1, 1], # M6: Pc, Pd ]) 

• The cycle time of each machine

# M1 M2 M3 M4 M5 M6 cycle_times = np.array([36, 36, 36, 36, 18, 18]) 

• The capacity of each machine (how much time/seconds it is available for)

#Seconds available per machine T_workday = 24 * 60**2 machine_capacities = np.full(num_machines, T_workday) 

• The utilisation (± tolerance) we want each machine to be at. It is uniform here, and defined with reference to a known 'full load' value.

# For enforcing uniform machine utilisation full_utilisation_reference = np.array([6720, 2880, 3840, 5760]).sum() #all machines were 100% here uniform_ideal = product_totals.sum() / full_utilisation_reference uniformity_tolerance = 0/100 #Accept 0% deviation either side of the ideal uniformity_lb = max(0., uniform_ideal - uniformity_tolerance) uniformity_ub = min(1., uniform_ideal + uniformity_tolerance) 

Framing the problem

Each product will be allotted to one or more machines. Since each product can, in theory, be spread over six machines, we can say that the task is to estimate six weights per product. So we want to estimate num_machines x num_products = 24 weights in total.

# We want to estimate a weight per product, per machine num_weights = num_machines * num_products 

These product weights (percentages) say how much of a product is given to each machine.

Defining constraints

The weights need to encode these constraints:

1. The weights are like percentages, so they must sum to 100% for each product

I set this up using an equality constraint (A_eq, b_eq).

# === Equality constraints === # For each product p, product_weights must sum to 100% # A_eq @ product_weights = b_eq A_eq = np.zeros((num_products, num_weights)) b_eq = np.ones(num_products) for p in range(num_products): #Compatibility matrix, limited to product p compatibilities_p = np.zeros_like(machine_product_compatibility) compatibilities_p[:, p] = machine_product_compatibility[:, p] #Sum of weights for machine p (A_eq @ x) = 1 (b_eq) A_eq[p, :] = compatibilities_p.ravel() 

The constraint works as follows. We use matrix notation to define a constraint per product, and put that into a single matrix. Each row is for a different product.

$$A_{eq} \cdot weights=b_{eq}$$

In this case, $A_{eq}$ is shaped (num_products, num_weights), and $b_{ub}$ is shaped (num_products,) since it must match the number of rows.

Concretely, for a given product $p$, the corresponding row from $A_{eq} \cdot weights$ will select the weights relevant to that product and sum them together: $\Sigma weights_{p}$. The matching element from $b_{eq}$ specifies that the result should equal 100% (hence "equality constraint"): $$\Sigma weights_{p}=100\%$$

2. The usage of each machine should be uniform

I define this using an upper-bound constraint.

# === Uniformity constraints === # Each machine should be utilised at a similar level: # unform_ideal - tolerance < utilisation % < unform_ideal + tolerance A_uniform = np.zeros((num_machines * 2, num_weights)) b_uniform = np.zeros((num_machines * 2,)) for m in range(num_machines): #Compatibility matrix, limited to machine m compatibilities_m = np.zeros_like(machine_product_compatibility) compatibilities_m[m, :] = machine_product_compatibility[m, :] seconds_per_x = compatibilities_m * product_totals[None, :] * cycle_times[m] #Upper bound # utilisation < (ideal + tolerance) * capacity A_uniform[m, :] = seconds_per_x.ravel() b_uniform[m] = uniformity_ub * machine_capacities[m] #Lower bound # utilisation > (ideal - tol)*capacity -> -sum < -(ideal - tol)*capacity A_uniform[2 * m + 1, :] = -seconds_per_x.ravel() b_uniform[2 * m + 1] = -uniformity_lb * machine_capacities[m] A_ub, b_ub = A_uniform, b_uniform 

3. For incompatible product-machine pairs, the weight must be zero

4. Each weight must be bounded between $(0, 1)$

I set these up using boundary constraints:

# === Bounds constraints === # Product weights bounded between (0, 1) for compatible product-machine pairs, else 0 bounds = [ (0., 1.) if machine_product_compatibility[m, p] else (0., 0.) for m in range(num_machines) for p in range(num_products) ] 

Solving for product_weights

We want to estimate product_weights in a way that minimises time utilisation across machines. The utilisation of each machine is given by summing unweighted_utilisation x product_weights per machine.

# Coefficients # Each row (machine) of unweighted_utilisation @ product_weights gives # utilisation distribution over products. # Will find product_weights st utilisation minimised (ie time minimised -> maximal throughput) unweighted_utilisation = ( machine_product_compatibility * cycle_times[:, None] * product_totals[None, :] ) #Estimate product_weights (returned in res.x) res = linprog(unweighted_utilisation.ravel(), A_ub, b_ub, A_eq, b_eq, bounds) display(res) 

Results for this example

The percentage weights come back as:

 #Pa Pb Pc Pd array([[ 36, 0, 0, 0], # M1 [ 0, 0, 0, 42], # M2 [ 0, 0, 62, 0], # M3 [ 36, 0, 0, 0], # M4 [ 29, 100, 0, 0], # M5 [ 0, 0, 38, 58]] # M6 ) 

This suggests that you can reserve machines 1 to 4 for a single product each, and use machines 5 and 6 for two products each.

The total available capacity machine_capacities.sum() (which is seconds in a day $\times$ number of machines), of which a third is utilised, since we have 1/3 of the full load in this example.

available time: 518,400 | utilised time: 172,800 

Fractional utilisation per machine shows each hitting the uniformity requirement:

[0.33 0.33 0.33 0.33 0.33 0.33] 

Non-overlapping case

I also tried for the second non-overlapping case (at 1/10th total). It succeeds if I allow a uniformity tolerance of 4%, by changing uniformity_tolerance to 4/100. The solution I get results in utilisation fractions of [0.14 0.14 0.12 0.06 0.11 0.06], which is close to yours though not identical. This means that you'd need to increase the tolerance until it works, which might not be practical. I think a squared-error minimisation approach might work better, since it will treat constraints as soft conditions rather than hard failure points.

---

Reproducible example

import numpy as np from scipy.optimize import linprog products = ['Pa', 'Pb', 'Pc', 'Pd'] machines = ['M1', 'M2', 'M3', 'M4', 'M5', 'M6'] num_products = len(products) num_machines = len(machines) # Pa Pb Pc Pd product_totals = np.array([6720, 2880, 3840, 5760]) // 3 machine_product_compatibility = np.array([ #Pa Pb Pc Pd [1, 0, 1, 0], # M1: Pa, Pc [0, 1, 0, 1], # M2: Pb, Pd [0, 1, 1, 0], # M3: Pb, Pc [1, 0, 0, 1], # M4: Pa, Pd [1, 1, 0, 0], # M5: Pa, Pb [0, 0, 1, 1], # M6: Pc, Pd ]) # M1 M2 M3 M4 M5 M6 cycle_times = np.array([36, 36, 36, 36, 18, 18]) #Seconds available per machine T_workday = 24 * 60**2 machine_capacities = np.full(num_machines, T_workday) # For enforcing uniform machine utilisation full_utilisation_reference = np.array([6720, 2880, 3840, 5760]).sum() #all machines were 100% here uniform_ideal = product_totals.sum() / full_utilisation_reference uniformity_tolerance = 0/100 #Accept 0% deviation either side of the ideal uniformity_lb = max(0., uniform_ideal - uniformity_tolerance) uniformity_ub = min(1., uniform_ideal + uniformity_tolerance) # We want to estimate a weight per product, per machine num_weights = num_machines * num_products # === Equality constraints === # For each product p, product_weights must sum to 100% # A_eq @ product_weights = b_eq A_eq = np.zeros((num_products, num_weights)) b_eq = np.ones(num_products) for p in range(num_products): #Compatibility matrix, limited to product p compatibilities_p = np.zeros_like(machine_product_compatibility) compatibilities_p[:, p] = machine_product_compatibility[:, p] #Sum of weights for machine p (A_eq @ x) = 1 (b_eq) A_eq[p, :] = compatibilities_p.ravel() # === Uniformity constraints === # Each machine should be utilised at a similar level: # unform_ideal - tolerance < utilisation % < unform_ideal + tolerance A_uniform = np.zeros((num_machines * 2, num_weights)) b_uniform = np.zeros((num_machines * 2,)) for m in range(num_machines): #Compatibility matrix, limited to machine m compatibilities_m = np.zeros_like(machine_product_compatibility) compatibilities_m[m, :] = machine_product_compatibility[m, :] seconds_per_x = compatibilities_m * product_totals[None, :] * cycle_times[m] #Upper bound # utilisation < (ideal + tolerance) * capacity A_uniform[m, :] = seconds_per_x.ravel() b_uniform[m] = uniformity_ub * machine_capacities[m] #Lower bound # utilisation > (ideal - tol)*capacity -> -sum < -(ideal - tol)*capacity A_uniform[2 * m + 1, :] = -seconds_per_x.ravel() b_uniform[2 * m + 1] = -uniformity_lb * machine_capacities[m] A_ub, b_ub = A_uniform, b_uniform # === Bounds constraints === # Product weights bounded between [0, 1] for compatible product-machine pairs, else 0 bounds = [ (0., 1.) if machine_product_compatibility[m, p] else (0., 0.) for m in range(num_machines) for p in range(num_products) ] # Coefficients # Each row (machine) of unweighted_utilisation @ product_weights gives # utilisation distribution over products. # Will find product_weights st utilisation minimised (ie time minimised -> maximal throughput) unweighted_utilisation = ( machine_product_compatibility * cycle_times[:, None] * product_totals[None, :] ) #Estimate product_weights (returned in res.x) res = linprog(unweighted_utilisation.ravel(), A_ub, b_ub, A_eq, b_eq, bounds) display(res) # === View results === #Product weights for each machine product_weights = res.x.reshape(num_machines, num_products) print( 'Percentage product weights for machines (rows) and products (columns):\n', product_weights.__mul__(100).round().astype(int), sep='' ) #Total utilisation print( '\navailable time: {:,}'.format(machine_capacities.sum()), ' | utilised time: {:,}'.format(res.fun) ) #Machine time per product print( '\nUtilisation time per machine (rows) per product (columns):\n', product_weights * unweighted_utilisation, '\n\nFractional utilisation per machine:\n', (product_weights * unweighted_utilisation).sum(axis=1) / machine_capacities, sep='' ) 

This approach frames the constraints as 'soft constraints' whose MSE is minimised using gradient descent.

Concepts

Key constraints

Weights sum to 1 for compatible machines, and individual weights are bounded between (0, 1). These two constraints are strictly enforced using a softmax transform (see the answer by @DevilXD for the original suggestion and commentary).

A soft constraint is mach_overcapacity_loss which tries to enforce machine utilisation $\leq$ machine capacity.

Machine grouping

Machine groups are determined using graph analysis; machines are linked if they share a product, which forms isolated groups of machines when there is no product overlap. The figure below visualises this for the first non-overlapping case.

Uniformity preference

There is also a soft constraint that encourages uniform machine utilisation for each group of machines, and acts more like a 'preference' - you can tune its impact via uniformity_strength. The algorithm discovers the right average for each group, so you don't need to specify any target value it.

Results

Overlapping case (full product volume)

The constraints and uniformity preference are easily respected, resulting in a low loss:

First non-overlapping case (1/10th volume)

The loss is higher here since the uniformity preference is harder to meet, but performance is nonetheless good.

The discrepancy between product totals and products processed is less than 1.

Percent utilisation per machine:

[14.667, 14.667, 14.667, 7.72, 8.16, 5.98] 

Group uniformity looks good. There is a wider spread in the second group:

Second non-overlapping case (1/10th volume)

Convergence quality/characteristics are similar to the preceding case.

Percent utilisation per machine: [14.034, 14.034, 11.932, 8.054, 8.099, 7.873]

There is some variation in terms of utilisation but it's relatively small:

In all cases the algorithm perfectly satisfies the machine capacity constraint (no machines > 100% capacity), which is why that loss component is zero (and doesn't show up on the log-scaled plots).

---

Reproducible example

Imports and supplied data

import numpy as np import torch from torch import nn from collections import defaultdict from matplotlib import pyplot as plt import networkx as nx np.set_printoptions(suppress=True) products = ['Pa', 'Pb', 'Pc', 'Pd'] machines = ['M1', 'M2', 'M3', 'M4', 'M5', 'M6'] num_products = len(products) num_machines = len(machines) # Pa Pb Pc Pd product_totals = torch.tensor([6720, 2880, 3840, 5760]).float() / 10 #Overlapping example machine_product_compatibility = torch.tensor([ #Pa Pb Pc Pd [1, 0, 1, 0], # M1: Pa, Pc [0, 1, 0, 1], # M2: Pb, Pd [0, 1, 1, 0], # M3: Pb, Pc [1, 0, 0, 1], # M4: Pa, Pd [1, 1, 0, 0], # M5: Pa, Pb [0, 0, 1, 1], # M6: Pc, Pd ]).float() #Non-overlapping, ex1 machine_product_compatibility = torch.tensor([ #Pa Pb Pc Pd [1, 0, 0, 0], # M1 [1, 0, 1, 0], # M2 [0, 0, 1, 0], # M3 [0, 0, 0, 1], # M4 [0, 1, 0, 1], # M5 [0, 1, 0, 0], # M6 ]).float() #Non-overlapping, ex2 machine_product_compatibility = torch.tensor([ #Pa Pb Pc Pd [1, 0, 0, 0], # M1 [1, 1, 0, 0], # M2 [0, 1, 0, 0], # M3 [0, 0, 0, 1], # M4 [0, 0, 1, 1], # M5 [0, 0, 1, 0], # M6 ]).float() # M1 M2 M3 M4 M5 M6 cycle_times = torch.tensor([36, 36, 36, 36, 18, 18]).float() #Seconds available per machine T_workday = 24 * 60**2 machine_capacities = torch.full((num_machines,), T_workday).float() seconds_products_matrix = cycle_times.outer(product_totals).float() 

Machine groups

def get_machine_groups(machine_product_compatibility): num_machines = machine_product_compatibility.shape[0] #Build a graph where machines/nodes are connected (edges) if they share a product G = nx.Graph() G.add_nodes_from(range(num_machines)) for i in range(num_machines): for j in range(i + 1, num_machines): if machine_product_compatibility[i].dot(machine_product_compatibility[j]): G.add_edge(i, j) # Extract connected components (groups of machines) machine_groups = [np.sort(list(group)) for group in nx.connected_components(G)] # Visualise graph pos = nx.spring_layout(G, seed=0) f, ax = plt.subplots(figsize=(5, 3), layout='tight') nx.draw_networkx_nodes(G, pos, node_color='lightblue', node_size=400, ax=ax) nx.draw_networkx_edges(G, pos, width=1.5, edge_color='gray', ax=ax) nx.draw_networkx_labels(G, pos, font_size=10, font_weight='bold', ax=ax) ax.set_title('Machine relationships via shared products', fontsize=12) ax.spines[:].set_visible(False) plt.show() return machine_groups machine_groups = get_machine_groups(machine_product_compatibility) print(f'Num. of machine groups is {len(machine_groups)}:', machine_groups) 

Gradient descent

# === Equality constraint === # For each product p, product_weights must sum to 100% # === Machine capacity constraints === # Each machine's utilisation can't exceed its capacity # === Uniformity constraint/preference === # Each machine group should be utilised at a similar level np.random.seed(0) torch.manual_seed(0) logits = nn.Parameter(torch.randn(num_machines, num_products)) optimiser = torch.optim.Adam([logits]) # How strongly to push for uniform utilsation across machines uniformity_strength = 1 def logits_to_weights(logits, mach_prod_compatibility): return ( logits .masked_fill(mach_prod_compatibility.eq(0), -torch.inf) .softmax(dim=0) ) losses = defaultdict(list) for _ in range(num_iter := 10_000): #Each weight is bound (0, 1), and only available for compatible pairs weights = logits_to_weights(logits, machine_product_compatibility) # Current utilisation per machine utilisation_per_mach = ( #Per machine, time spent on each product seconds_products_matrix.mul(weights) #Each machine's total time (add up all product times) .sum(dim=1) #Express as fraction of each machine's capacity .div(machine_capacities) ) #Machine utilisation must not exceed machine's # capacity (ie. normalised utilisation < 1). # Being under-capacity is ok - exlude using relu mach_overcapacity_loss = utilisation_per_mach.sub(1.).relu().square().mean() #Encourage uniformity, per machine group group_uniformity_losses = torch.stack([ #Deviation of machines from group average (utilisation_per_mach[mach_ixs] - utilisation_per_mach[mach_ixs].mean()) #Scalar MSE .square().mean() #Repeat for each group of machine for mach_ixs in machine_groups ]) group_uniformity_loss = group_uniformity_losses.mean().mul(uniformity_strength) #Total loss loss = mach_overcapacity_loss + group_uniformity_loss #Optimisation step optimiser.zero_grad() loss.backward() optimiser.step() # === Record metrics === losses['mach_overcapacity'].append(mach_overcapacity_loss.item()) losses['group_uniformity'].append(group_uniformity_loss.item()) losses['loss'].append(loss.item()) #Stop at relatively low loss if loss < 1e-7: print('Stopping at low loss', loss.item()) break # === Plot === f, ax = plt.subplots(figsize=(5, 2.8)) for name, loss_ in losses.items(): if name == 'loss': continue ax.plot(loss_, label=name) ax.plot(losses['loss'], lw=3, ls=(0, (2, 2)), label='loss') ax.legend(framealpha=0) ax.spines[['top', 'right']].set_visible(False) ax.semilogy(True) ax.set_ylim(None, 1e-2) ax.grid(axis='y', which='major', linestyle=':') ax.set( title='Loss convergence', xlabel='iteration', ylabel='total loss' ) plt.show() 

Plot and report results

# === View results === # Product weights (compatible) for each machine weights = logits_to_weights(logits, machine_product_compatibility).detach() print( '\nPercentage product weights for compatible machines (rows) and products (columns):\n', weights.mul(100).numpy().round(2), sep='' ) #Total utilisation print( '\navailable time: {:,}'.format(machine_capacities.sum()), ' | utilised time: {:,}'.format(seconds_products_matrix.ravel() @ weights.ravel()) ) #Machine time per product print( '\nUtilisation time per machine (rows) per product (columns):\n', weights * seconds_products_matrix, '\n\nPercent utilisation per machine:\n', weights.mul(seconds_products_matrix).sum(dim=1) / machine_capacities * 100, sep='' ) f, ax = plt.subplots(figsize=(5, 3), layout='tight') for mach_ixs in machine_groups: ax.plot(mach_ixs, utilisation_per_mach[mach_ixs].detach(), linestyle='none', marker='d') ax.plot(mach_ixs, [utilisation_per_mach[mach_ixs].detach().mean(),]*len(mach_ixs), color='black') ax.set( xlabel='machine index', ylabel='utilisation', title='Machine utilisation and group averages' ) ax.spines[['top', 'right']].set_visible(False) 

Reported to stdout:

Percentage product weights for compatible machines (rows) and products (columns): [[50.12 0. 0. 0. ] [49.88 0.57 0. 0. ] [ 0. 99.43 0. 0. ] [ 0. 0. 0. 33.56] [ 0. 0. 1.58 66.44] [ 0. 0. 98.42 0. ]] available time: 518,400.0 | utilised time: 55,319.453125 Utilisation time per machine (rows) per product (columns): tensor([[12125.4277, 0.0000, 0.0000, 0.0000], [12066.5742, 58.7491, 0.0000, 0.0000], [ 0.0000, 10309.2510, 0.0000, 0.0000], [ 0.0000, 0.0000, 0.0000, 6958.8979], [ 0.0000, 0.0000, 109.3414, 6888.5508], [ 0.0000, 0.0000, 6802.6587, 0.0000]]) Percent utilisation per machine: tensor([14.0341, 14.0339, 11.9320, 8.0543, 8.0994, 7.8734])