[Python] Need hints for tuning scipy .linprog/.milp

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 * prods_vars for i in prods]),"total_cost")
prob += (pulp.lpSum([exps * prods_vars for i in prods]) >= target_value, "lowerBound")
prob += (pulp.lpSum([exps * prods_vars 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.
  

Continue reading...