Need hints for tuning scipy .linprog/.milp

Question

I have a task to optimize the cost of a product, for example, by finding the optimal combination of cans of paint of different volumes and prices for coloring a given square.

The simplest example is:

pdict = {'Can_9ltr':[9, 0.17, 4870],
            'Can_4.5ltr':[4.5, 0.17, 2910],
            'Can_1ltr':[1, 0.17, 632],
            'Can_2.25ltr':[2.25, 0.17, 1790]}

where list[0] is volume of particular can, list[1] is paint consumption per square meter and list[2] is price per can.

In fact, I got a fully working solution using the pulp library:

{'Target_value': 30,
 'Optimal_set': {'Can_1ltr': {'CPU': 632, 'Amount': 1.0, 'Cost': 632.0},
  'Can_2.25ltr': {'CPU': 1790, 'Amount': 0.0, 'Cost': 0.0},
  'Can_4.5ltr': {'CPU': 2910, 'Amount': 1.0, 'Cost': 2910.0},
  'Can_9ltr': {'CPU': 4870, 'Amount': 0.0, 'Cost': 0.0}},
 'Cost_total': 3542.0}

Code below:

import pulp
import math

pdict = {'Can_9ltr':[9, 0.17, 4870],
            'Can_4.5ltr':[4.5, 0.17, 2910],
            'Can_1ltr':[1, 0.17, 632],
            'Can_2.25ltr':[2.25, 0.17, 1790]}
square = 30

def calculator(prod_dict, target_value):
    prods = list(prod_dict.keys())
    costs = {}
    for key in prod_dict:
        costs[key] = prod_dict[key][2]
    exps = {}
    for key in prod_dict:
        exps[key] = prod_dict[key][0]/prod_dict[key][1]
    upperBound = math.ceil(target_value/max([prod_dict[a][0]/prod_dict[a][1] for a in prod_dict]))\
    *max([prod_dict[a][0]/prod_dict[a][1] for a in prod_dict])
    prob = pulp.LpProblem('optimal_count', pulp.LpMinimize)
    prods_vars = pulp.LpVariable.dicts('prods', prods, lowBound=0, cat=pulp.LpInteger)

    prob += (pulp.lpSum([costs[i] * prods_vars[i] for i in prods]),"total_cost")
    prob += (pulp.lpSum([exps[i] * prods_vars[i] for i in prods]) >= target_value, "lowerBound")
    prob += (pulp.lpSum([exps[i] * prods_vars[i] for i in prods]) <= upperBound, "upperBound")

    prob.solve()
    optimal_set = {}
    for v in prob.variables():
        optimal_set[v.name.replace('prods_', '')] = {'CPU': prod_dict[v.name.replace('prods_', '')][2], 
                                                     'Amount': v.varValue, 
                                                     'Cost': prod_dict[v.name.replace('prods_', '')][2]*v.varValue}

    result_dict = {}
    result_dict['Target_value'] = target_value
    result_dict['Optimal_set'] = optimal_set
    result_dict['Cost_total'] = sum([optimal_set[a]['Cost'] for a in optimal_set])

    return result_dict

However, I would like to reproduce this solution using scipy.optimize.linprog or scipy.optimize.milp, but I have encountered difficulty understanding the settings for these methods.

At the moment, all I have managed to get with linprog and milp is a solution that satisfies only one requirement - sufficiency for coloring a given area

linprog code:

from scipy.optimize import linprog
import math
pdict = {'Can_9ltr':[9, 0.17, 4870],
            'Can_4.5ltr':[4.5, 0.17, 2910],
            'Can_1ltr':[1, 0.17, 632],
            'Can_2.25ltr':[2.25, 0.17, 1790]}
square = 30
c = []
for key in pdict:
    c.append(pdict[key][2])
exps = []
for key in pdict:
    exps.append(pdict[key][0]/pdict[key][1])
A_ub = [[(1/c[0])*exps[0], (1/c[1])*exps[1], (1/c[2])*exps[2], (1/c[3])*exps[3]],
        [(1/c[0])*exps[0]*-1, (1/c[1])*exps[1]*-1, (1/c[2])*exps[2]*-1, (1/c[3])*exps[3]*-1]]
b_ub = [math.ceil(square/max(exps))*max(exps), square*-1]
integrality=[1, 1, 1, 1]
res = linprog(c, A_ub, b_ub, integrality=integrality)
print(list(res.x))

milp code:

from scipy.optimize import milp, LinearConstraint
import math

pdict = {'Can_9ltr':[9, 0.17, 3980],
            'Can_4.5ltr':[4.5, 0.17, 2536],
            'Can_1ltr':[1, 0.17, 517],
            'Can_2.25ltr':[2.25, 0.17, 1470]}
square = 30
c = []
for key in pdict:
    c.append(pdict[key][2])
exps = []
for key in pdict:
    exps.append(pdict[key][0]/pdict[key][1])    
Ac = [[(1/c[0])*exps[0], (1/c[1])*exps[1], (1/c[2])*exps[2], (1/c[3])*exps[3]],
        [(1/c[0])*exps[0], (1/c[1])*exps[1], (1/c[2])*exps[2], (1/c[3])*exps[3]]]
b_u = [math.ceil(square/max(exps))*max(exps)]
b_l = [square]
integrality = [1, 1, 1, 1]
constraints = LinearConstraint(Ac, b_l, b_u)
res = milp(c=c, constraints=constraints, integrality=integrality)
print(list(res.x))

Solution same for linprog and milp:

[0.0, 0.0, 3224.0, 0.0]

Two key issues:

I can't figure out how to force integer coefficient - 3224.0 is 5.1 of Can_1ltr unit cost. Option integrality = [1, 1, 1, 1] - I'm not sure if I'm using this correctly, or if it's working as I expected.
I don't understand why both methods return solution with only one unit taken into account instead of minimizing whole function.

When you call calculator, what is the value of the target_value parameter? 30 (square)? — Reinderien
– Reinderien, Commented Dec 19, 2024 at 22:55

Reinderien · Accepted Answer · 2024-12-19 23:58:34Z

0

This is equivalent. Pandas is optional but helps understand the solution.

import numpy as np
import pandas as pd
import math

from scipy.optimize import milp, Bounds, LinearConstraint

df = pd.DataFrame(
    data={
        'volume': (9, 4.5, 1, 2.25),  # litres
        'consume': (0.17, 0.17, 0.17, 0.17),  # ?? per m2
        'price': (4870, 2910, 632, 1790),  # ?? per can
    },
)
df['exp'] = df['volume']/df['consume']
n = len(df)
target_value = 30
max_prod = df['exp'].max()
max_exp = math.ceil(target_value / max_prod)*max_prod

result = milp(
    c=df['price'],
    integrality=np.ones(n, dtype=np.uint8),
    bounds=Bounds(lb=np.zeros(n)),
    constraints=LinearConstraint(A=df['exp'], lb=target_value, ub=max_exp),
)
assert result.success, result.message
df['prod'] = result.x
print(df)

   volume  consume  price        exp  prod
0    9.00     0.17   4870  52.941176   0.0
1    4.50     0.17   2910  26.470588   1.0
2    1.00     0.17    632   5.882353   1.0
3    2.25     0.17   1790  13.235294   0.0

answered Dec 19, 2024 at 23:58

Reinderien

16.7k9 gold badges56 silver badges92 bronze badges

Sign up to request clarification or add additional context in comments.

3 Comments

Артур Дутов Dec 20, 2024 at 13:31

Thank you, it works! Seems like I’ve tricked myself with constraints and missed the importance of bounds option. May I ask you to clarify is it the only and specific solution to build integrality variable with np methods?

Reinderien Dec 20, 2024 at 13:37

These are not numpy methods; they're scipy methods. This is essentially the only method to ask for integral variables in scipy LP support.

Reinderien Dec 20, 2024 at 13:41

If you're asking about the form of the integrality argument, there are other ways to pass scalars that get broadcast, but this form with this dtype is the form that makes scipy do the least amount of work.

Collectives™ on Stack Overflow

Need hints for tuning scipy .linprog/.milp

1 Answer 1

3 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

3 Comments

Your Answer

Sign up or log in

Post as a guest

Related