I have a some materials made from chemical elements having a price for each kilo gram:
A: { cu: 0.02, mg: 0.3, nc: 0.005, price: 1000 }
This means 1 kilo gram of material A is made from 0.02g "cu", 0.3g "mg", 0.005g "nc", and has a price of 1000$. I have some other similar materials:
A: { cu: 0.02, mg: 0.3, nc: 0.005, price: 1000 }
B: { cu: 0.033, mg: 0.2, nc: 0.00, price: 1500 }
C: { cu: 0.05, mg: 0.35, nc: 0.00, price: 1700 }
I want to build a product, lets call it Pr_1, having 10 kilo grams weight, made from materials A, B, and/or C, so the Pr_1 has minimum price. The Pr_1 has some constraints too:
Pr_1: {minCu: 0.00, maxCu: 1.2, minMg: 1.2, maxMg: 1.9, minNc: 0, maxNc: 0}
This means 10 kilo gram of Pr_1 should have at least 0.00g and at most 1.2g "cu", also it should have at least 1.2g and at most 1.9g "mg", and it shouldn't have any "nc". As I mentioned, the Pr_1 has to have minimum price among all possible combinations of materials A, B, and C.
So, Let's denote the weights of Pr_1 made from materials A, B, and C as x_A, x_B, and x_C, respectively. The objective is to minimize the total cost:
Minimize:
1000*x_A + 1500*x_B + 1700*x_C
and we have 3 inequality equations:
0.00 <= 0.02*x_A + 0.033*x_B + 0.05*x_C <= 1.2 // cu constraint
1.2 <= 0.03*x_A + 0.2*x_B + 0.35*x_C <= 1.9 // mg constraint
0.005*x_A + 0.00*x_B + 0.00*x_C = 0 // nc constraint
I can not figure out how can I solve this using Scipy and Python.
from scipy.optimize import linprog
# Coefficients of the objective function (costs)
c = [1000, 1500, 1700]
# Coefficients of the inequality constraints matrix
A = [
[-0.02, -0.033, -0.05], # Cu constraint
[-0.3, -0.2, -0.35], # Mg constraint
[-0.005, 0, 0] # Nc constraint
]
# Right-hand side of the inequality constraints
b = ???? # I don't know what this should be
# Bounds for each variable (x_A, x_B, x_C)
x_bounds = (0, None)
# Solve the linear programming problem
result = linprog(c, A_ub=A, b_ub=b, bounds=[x_bounds, x_bounds, x_bounds], method='highs')
# Print the results
print("Optimal weights for Pr_1 from A, B, and C:", result.x)
print("Minimum cost of Pr_1:", result.fun)
As I said in the comments, with your current constraints the problem is infeasible.
With the same constraints you could produce 9 kg. Also, I believe that you should be more careful and explicit with your units, so I demonstrate dimensional management using Pandas. In the end this produces a much more legible output.
Also, at a guess, your Nc implies nickel which is actually Ni.
import numpy as np
import pandas as pd
from scipy.optimize import milp, Bounds, LinearConstraint
materials = pd.DataFrame(
data={'dollars_per_kg': (1_000, 1_500, 1_700)},
index=pd.Index(name='material_name', data=('A', 'B', 'C')),
)
materials['dollars_per_g'] = 1e-3*materials['dollars_per_kg']
components = pd.DataFrame(
data={
'min_g': (0.0, 1.2, 0),
'max_g': (1.2, 1.9, 0),
},
index=pd.Index(name='component_name', data=('Cu', 'Mg', 'Ni')),
)
material_components = pd.DataFrame(
data={
'g_per_kg': (
0.0200, 0.30, 0.005,
0.0033, 0.20, 0.000,
0.0500, 0.35, 0.000,
),
},
index=pd.MultiIndex.from_product(
iterables=(materials.index, components.index),
),
)
material_components['concentration'] = 1e-3 * material_components['g_per_kg']
mass_kg = 9 # 10 not feasible
mass_g = 1e3*mass_kg
'''
Minimize cost = (A, B, C mass g).(A, B, C dollars per g)
s.t.
A, B, C mass g >= 0
sum(A, B, C mass g) = mass_g
component min <= sum by component(A, B, C mass g) < component max
'''
mass_sum_constraint = LinearConstraint(
A=np.ones(len(materials)),
lb=mass_g, ub=mass_g,
)
component_constraint = LinearConstraint(
A=material_components['concentration'].unstack('material_name'),
lb=components['min_g'],
ub=components['max_g'],
)
result = milp(
c=materials['dollars_per_g'],
bounds=Bounds(lb=0),
constraints=(
mass_sum_constraint,
component_constraint,
),
)
if not result.success:
raise ValueError(result.message)
print(result.message, end='\n\n')
materials['mass_g'] = result.x
materials['dollars'] = result.x * materials['dollars_per_g']
material_components['mass_g'] = material_components['concentration'] * materials['mass_g']
components['mass_g'] = material_components['mass_g'].groupby('component_name').sum()
print('Materials:')
print(materials, end='\n\n')
print('Components:')
print(components, end='\n\n')
Optimization terminated successfully. (HiGHS Status 7: Optimal)
Materials:
dollars_per_kg dollars_per_g mass_g dollars
material_name
A 1000 1.0 -0.0 -0.0
B 1500 1.5 9000.0 13500.0
C 1700 1.7 0.0 0.0
Components:
min_g max_g mass_g
component_name
Cu 0.0 1.2 0.0297
Mg 1.2 1.9 1.8000
Ni 0.0 0.0 0.0000
To find the maximum total mass, change your objective:
import numpy as np
import pandas as pd
from scipy.optimize import milp, Bounds, LinearConstraint
materials = pd.DataFrame(
data={'dollars_per_kg': (1_000, 1_500, 1_700)},
index=pd.Index(name='material_name', data=('A', 'B', 'C')),
)
materials['dollars_per_g'] = 1e-3*materials['dollars_per_kg']
components = pd.DataFrame(
data={
'min_g': (0.0, 1.2, 0),
'max_g': (1.2, 1.9, 0),
},
index=pd.Index(name='component_name', data=('Cu', 'Mg', 'Ni')),
)
material_components = pd.DataFrame(
data={
'g_per_kg': (
0.0200, 0.30, 0.005,
0.0033, 0.20, 0.000,
0.0500, 0.35, 0.000,
),
},
index=pd.MultiIndex.from_product(
iterables=(materials.index, components.index),
),
)
material_components['concentration'] = 1e-3 * material_components['g_per_kg']
mass_kg = 9 # 10 not feasible
mass_g = 1e3*mass_kg
'''
Minimize cost = (A, B, C mass g).(A, B, C dollars per g)
Or: maximize total mass
s.t.
A, B, C mass g >= 0
sum(A, B, C mass g) = mass_g
component min <= sum by component(A, B, C mass g) < component max
'''
mass_sum_constraint = LinearConstraint(
A=np.ones(len(materials)),
lb=mass_g, ub=mass_g,
)
component_constraint = LinearConstraint(
A=material_components['concentration'].unstack('material_name'),
lb=components['min_g'],
ub=components['max_g'],
)
result = milp(
# c=materials['dollars_per_g'], # to minimize price
c=-np.ones_like(materials['dollars_per_g']), # to maximize total mass
bounds=Bounds(lb=0),
constraints=(
# mass_sum_constraint, # to ask for a fixed total mass
component_constraint,
),
)
if not result.success:
raise ValueError(result.message)
print(result.message, end='\n\n')
materials['mass_g'] = result.x
materials['dollars'] = result.x * materials['dollars_per_g']
material_components['mass_g'] = material_components['concentration'] * materials['mass_g']
components['mass_g'] = material_components['mass_g'].groupby('component_name').sum()
print('Materials:')
print(materials, end='\n\n')
print('Components:')
print(components, end='\n\n')
print('Material components:')
print(material_components)
Optimization terminated successfully. (HiGHS Status 7: Optimal)
Materials:
dollars_per_kg dollars_per_g mass_g dollars
material_name
A 1000 1.0 0.0 0.0
B 1500 1.5 9500.0 14250.0
C 1700 1.7 0.0 0.0
Components:
min_g max_g mass_g
component_name
Cu 0.0 1.2 0.03135
Mg 1.2 1.9 1.90000
Ni 0.0 0.0 0.00000
Material components:
g_per_kg concentration mass_g
material_name component_name
A Cu 0.0200 0.000020 0.00000
Mg 0.3000 0.000300 0.00000
Ni 0.0050 0.000005 0.00000
B Cu 0.0033 0.000003 0.03135
Mg 0.2000 0.000200 1.90000
Ni 0.0000 0.000000 0.00000
C Cu 0.0500 0.000050 0.00000
Mg 0.3500 0.000350 0.00000
Ni 0.0000 0.000000 0.00000