Nonlinear Programming with Inequality Constraints
==============
Imagine that you work out in a factory that produces some kind of composite products.
Each product has its own composition of parts.
Some parts are bought, and some are made at your factory.
Those products that are made in your factory have their own production route.
You know the number of your machines in production and the fund of their working hours.
You have a sales plan for the year It contains an order for any product that you produce, terms, quantity of required products and other information.
Having previously received the composition of each required product needed for sale, and having processed the information from the sales plan, you get two tables.
Data on how much time and in what batches each product is processed and on which machine is in the df_t_o_p table.
This data is arranged in an mxn matrix(variable name in code M_cons), where m = 87 are rows (machines) and n = 864 are columns (products).
The data table stores has the following data:
-Definition - the name of the children;
-Quantity - the required quantity of the part;
-Dim - units of measure for the part (it doesn't matter at the moment);
-p_i, a_i, b_i - coefficients calculated during data preprocessing;
-Cost - the cost of each part associated with its dimension (That is, there is a kind of automatic normalization here, but we don’t think about it yet.)
Model
Problem statement: It is necessary to find the optimal number of parts that we will produce and store in advance in order to ensure the possibility of fulfilling the entire sales plan.
The goal function is actually the minimization of the storage cost function.
The goal function and constraints are described in the code.
from functools import reduce
import cvxpy
from scipy.optimize import Bounds, minimize
from cvxpy import multiply, minimum, maximum, power
import cvxpy as cp
import numpy as np
import pandas as pd
Path_data = 'C:\\your_path\\data.xlsx'
Path_df_t_o_p = 'C:\\your_path\\df_t_o_p.xlsx'
data_xlsx = pd.ExcelFile(Path_data)
df_t_o_p_xlsx = pd.ExcelFile(Path_df_t_o_p)
_data = pd.DataFrame(data_xlsx.parse(data_xlsx.sheet_names[0]))
_df_t_o_p = pd.DataFrame(df_t_o_p_xlsx.parse(df_t_o_p_xlsx.sheet_names[0])).set_index('Index')
data = {
'Definition': np.array(_data['Definition']),
'q_i' : np.array(_data['Quantity']),
'Dim' : np.array(_data['Dim']),
'p_i' : np.array(_data['p_i']),
'a_i' : np.array(_data['a_i']),
'b_i' : np.array(_data['b_i']),
'c_i' : np.array(_data['Cost'])
}
# M_cons * X <= T_l
# Matrix of ineq_constraints
M_cons = np.array(_df_t_o_p)
# Time fund for each machines mxn = 87x864
# m = 87 - machines
T_l = np.ones(len(_df_t_o_p.index)) * 60 * 360
# n = 864 - details
n = len(data['Definition'])
# 0 <= x_i <= q_i
l_bound = np.zeros(n)
r_bound = data['q_i'] + 0.00000001
bounds = Bounds(l_bound, r_bound)
ineq_constraints = {
'type': 'ineq',
'fun' : lambda x: T_l - M_cons.dot(x),
'jac' : lambda x: -1 * M_cons
}
x0 = data['q_i'] / 2
def f_scipy(x):
'''
My model is:
OBJECTIVE FUNCTION:
f(x) = \sum_i term1 + term 2 + term3: f(x) ---> min
where for all i :
term1 = c_i * a_i * x_i
term2 = c_i * b_i * (x_i - q_i) * (x_i - q_i > 0) = c_i * b_i * max{x_i - q_i, 0}
term3 = c_i^k * (q_i - x_i)^2 * (q_i - x_i > 0) = c_i^k * (q_i - x_i) * max{q_i - x_i, 0}
, k - must be chosen such that if b_i = 0, then x_i Should be equal to q_i in the model solution.
---------------------------------------------------------------------------------------------------------
CONSTRAINTS:
X = (x_0, x_1, .... , x_n)
for all i 0 <= x_i <= q_i
actually the upper bound is Infinity.
According to the logic of our objective function,
the solution should not give x_i > q_i.
Explanations a little later.
M_cons * X <= T_l
:param x: np.array len(x) = n
:return: float
'''
term1 = data['c_i'] * data['a_i'] * x
term2 = data['c_i'] * data['b_i'] * (x - data['q_i']) * (x - data['q_i'] > 0)
term3 = data['c_i']**10 * (data['q_i'] - x)**2 * (data['q_i'] - x > 0)
return np.sum(term1 + term2 + term3)
def model_scipy():
res = minimize(f_scipy, x0, bounds=bounds)
tmp = {'Definition': [], 'x_i': []}
for name, x in zip(data['Definition'], res.x):
tmp['Definition'].append(name)
tmp['x_i'].append(x)
tmp.update({'Quantity': data['q_i'],
'Dim': data['Dim'],
'b_i': data['b_i'],
'c_i': data['c_i']})
df = pd.DataFrame.from_dict(tmp)
if reduce(lambda x, y: x*y, res.x == x0):
print('res.x = x0')
else:
print(f"ineq in x0: \n{T_l - M_cons.dot(x0)}\n\n"
f"ineq in res.x: \n{T_l - M_cons.dot(res.x)}\n\n"
f"f(x0) = {f_scipy(x0)}\t"
f"f(res.x) = {f_scipy(res.x)}")
return df
def f_cvxpy(x):
'''
My model is:
OBJECTIVE FUNCTION:
f(x) = \sum_i term1 + term 2 + term3: f(x) ---> min
where for all i :
term1 = c_i * a_i * x_i
term2 = c_i * b_i * (x_i - q_i) * (x_i - q_i > 0) = c_i * b_i * max{x_i - q_i, 0}
term3 = c_i^k * (q_i - x_i)^2 * (q_i - x_i > 0) = c_i^k * (q_i - x_i) * max{q_i - x_i, 0}
, k - must be chosen such that if b_i = 0, then x_i Should be equal to q_i in the model solution.
---------------------------------------------------------------------------------------------------------
CONSTRAINTS:
X = (x_0, x_1, .... , x_n)
for all i 0 <= x_i <= q_i
actually the upper bound is Infinity.
According to the logic of our objective function,
the solution should not give x_i > q_i.
Explanations a little later.
M_cons * X <= T_l
:param x: np.array len(x) = n
:return: float
'''
result = 0
for i in range(len(data['Definition'])):
a_i = data['a_i'][i]
b_i = data['b_i'][i]
q_i = data['q_i'][i]
c_i = data['c_i'][i]
x_i = x[i]
term1 = multiply(c_i, multiply(a_i, x_i))
term2 = multiply(c_i, multiply(b_i, maximum(x_i - q_i, 0)))
term3 = multiply(power(c_i, 10), power(maximum(q_i - x_i, 0), 2))
result += term1 + term2 + term3
return result
def model_cvxpy():
x = cp.Variable(n)
print(cp.installed_solvers())
object = cp.Minimize(f_cvxpy(x))
prob = cp.Problem(object, [l_bound <= x, x <= r_bound, M_cons @ x <= T_l])
prob.solve(verbose=True, solver=cp.CVXOPT)
print('\nThe optimal values is', prob.value)
print('A solution x is')
print(x.value)
print('A dual solution corresponding to the inequality constraints is')
print(prob.constants[0].dual_value)
I try to use two libraries to solve the optimization problem - scipy and cvxpy.
The point is that my objective function is expected to be convex. Therefore, in the Scipy library, you can use the SLSQP method without fear that we will fall into a local extremum.
For this reason, I tried the second library, since the problem is convex.
My result after calling this code
import pandas as pd
from Model import model_scipy, model_cvxpy
def print_pd_df(df):
with pd.option_context('display.max_rows', None, 'display.max_columns', None, 'display.width', None):
print(df)
if __name__ == '__main__':
print_pd_df(model_scipy())
model_cvxpy()
It turns out like this:
f(x0) = 8.050283944386748e+49 f(res.x) = 1.0724365904778775e+41
Definition x_i Quantity Dim b_i c_i
0 100740 1.476151e-02 2.952192e-02 М2 0.000000e+00 6250.000600
1 101096 2.412800e-01 4.825600e-01 М 0.000000e+00 1426.512333
2 101830 8.800000e-02 1.760000e-01 М 0.000000e+00 578.696400
3 101831 1.240000e-01 2.480000e-01 М 0.000000e+00 1031.928820
4 102336 1.160004e-01 2.320000e-01 М 0.000000e+00 4954.878020
5 102488 3.910975e+04 7.821950e+04 ШТ 2.190669e+03 8.000000
6 102505 3.910975e+04 7.821950e+04 М2 2.190669e+03 34.725000
7 102583 1.955487e+03 3.910975e+03 КГ 2.190669e+03 130.000000
8 102649 9.827521e-02 9.827520e-02 М2 0.000000e+00 56493.969550
9 104034 7.500000e-02 1.500000e-01 М 0.000000e+00 74.950460
10 104518 5.000000e-01 1.000000e+00 ШТ 0.000000e+00 18.091700
11 104595 1.665000e+00 3.330000e+00 М 0.000000e+00 47.517200
12 104596 1.005000e+00 2.010000e+00 М 0.000000e+00 76.753333
13 104597 4.350000e-01 8.700000e-01 М 0.000000e+00 75.168600
14 104598 1.530000e+00 3.060000e+00 М 0.000000e+00 62.760240
15 104650 3.635030e-02 7.270000e-02 М 0.000000e+00 5378.125000
16 104651 1.083059e-02 2.160000e-02 М 0.000000e+00 9637.638850
17 110033 6.000000e-02 1.200000e-01 М 0.000000e+00 12.720611
18 110034 6.000000e-02 1.200000e-01 М 0.000000e+00 11.985950
19 110037 6.000000e-02 1.200000e-01 М 0.000000e+00 11.091833
20 110038 6.000000e-02 1.200000e-01 М 0.000000e+00 10.782600
21 110039 6.000000e-02 1.200000e-01 М 0.000000e+00 8.230700
22 111622 9.500000e-02 1.900000e-01 М 0.000000e+00 34.055200
23 111625 3.150000e-01 6.300000e-01 М 0.000000e+00 80.854000
24 111626 1.047500e+00 2.095000e+00 М 0.000000e+00 70.287500
25 111627 1.000000e+00 2.000000e+00 М 0.000000e+00 43.604667
26 111628 5.875517e+03 1.175103e+04 М 2.190669e+03 35.818400
27 111630 5.000000e-02 1.000000e-01 М 0.000000e+00 303.040000
28 112026 2.542671e+05 5.085342e+05 М 2.190669e+03 0.000000
29 112276 2.738235e+03 5.476470e+03 М 2.190669e+03 569.310000
30 113066 1.700000e+00 3.400000e+00 М 0.000000e+00 227.298333
31 114201 1.467500e+00 2.935000e+00 М 0.000000e+00 59.539600
32 114254 1.799814e+04 3.599627e+04 М 2.190669e+03 39.824875
33 114400 1.661000e-01 3.322000e-01 М 0.000000e+00 1025.462750
34 117365 2.750000e-02 5.500000e-02 М 0.000000e+00 53.830000
35 118160 2.000000e+00 4.000000e+00 ШТ 0.000000e+00 2.364400
36 118169 2.000000e+00 4.000000e+00 ШТ 0.000000e+00 1.933300
37 118266 3.500000e-01 7.000000e-01 М 0.000000e+00 64.716700
38 118268 2.000000e-01 4.000000e-01 М 0.000000e+00 63.644083
39 118613 1.000000e-01 2.000000e-01 М 0.000000e+00 38.017380
40 118743 1.000000e-01 2.000000e-01 М 0.000000e+00 39.895000
41 120619 6.000000e+00 1.200000e+01 М 0.000000e+00 1.500000
42 121514 2.000000e+00 4.000000e+00 ШТ 0.000000e+00 4.150000
43 122299 5.000000e-01 1.000000e+00 ШТ 0.000000e+00 3.872900
44 122710 5.000000e-02 1.000000e-01 М 0.000000e+00 27.270000
45 124114 5.000000e-01 1.000000e+00 ШТ 0.000000e+00 329.237300
46 124200 2.000000e+00 4.000000e+00 ШТ 0.000000e+00 2.975000
47 124201 8.000000e+00 1.600000e+01 ШТ 0.000000e+00 3.518233
48 124361 1.955890e+04 3.911780e+04 М 2.190669e+03 1.075000
49 125242 1.243001e-02 1.243000e-02 М2 0.000000e+00 45359.576350
50 125243 3.383521e-02 3.383520e-02 М2 0.000000e+00 63989.285700
51 125244 3.294901e-02 3.294900e-02 М2 0.000000e+00 39410.770200
52 125852 3.129560e+05 6.259120e+05 ШТ 2.190669e+03 1.380025
53 125853 8.000000e+00 1.600000e+01 ШТ 0.000000e+00 1.339000
54 125856 1.955942e+05 3.911885e+05 ШТ 2.190669e+03 1.347200
55 125857 7.350000e+01 1.470000e+02 ШТ 0.000000e+00 1.081967
56 125858 1.000000e+00 2.000000e+00 ШТ 0.000000e+00 1.341200
57 125859 1.500000e+00 3.000000e+00 ШТ 0.000000e+00 1.367400
58 125860 1.000000e+01 2.000000e+01 ШТ 0.000000e+00 1.591250
59 125863 3.000000e+00 6.000000e+00 ШТ 0.000000e+00 2.283300
60 125864 4.700000e+01 9.400000e+01 ШТ 0.000000e+00 1.948950
61 125867 2.000000e+00 4.000000e+00 ШТ 0.000000e+00 4.916700
62 125870 1.400000e+01 2.800000e+01 ШТ 0.000000e+00 11.322233
63 125903 5.000000e-01 1.000000e+00 ШТ 0.000000e+00 16.070800
64 126780 1.000000e+01 2.000000e+01 ШТ 0.000000e+00 0.250000
65 127259 1.000000e+00 2.000000e+00 ШТ 0.000000e+00 0.741700
66 127601 1.000000e+00 2.000000e+00 ШТ 0.000000e+00 1.754167
67 129605 1.850691e-02 3.700000e-02 М 0.000000e+00 7870.538000
68 130383 8.000000e+00 1.600000e+01 ШТ 0.000000e+00 3.619775
69 130560 1.010400e-01 2.020800e-01 М 0.000000e+00 1126.219600
70 130978 1.026615e+00 2.000000e+00 ШТ 0.000000e+00 12058.999400
71 130998 3.911725e+04 7.823450e+04 ШТ 2.190669e+03 19.000000
72 131072 1.550000e+01 3.100000e+01 ШТ 0.000000e+00 0.193340
73 131073 1.173352e+05 2.346705e+05 ШТ 2.190669e+03 0.391700
74 131075 4.800000e+01 9.600000e+01 ШТ 0.000000e+00 1.241633
75 131108 5.000000e-01 1.000000e+00 ШТ 0.000000e+00 2.950000
76 131194 1.260000e-02 2.520000e-02 М 0.000000e+00 450.270200
0.000000
393 159925 5.000000e-01 1.000000e+00 ШТ 0.000000e+00 380.000000
394 159991 2.000000e+00 4.000000e+00 ШТ 0.000000e+00 0.000000
395 160145 5.000000e+00 1.000000e+01 ШТ 0.000000e+00 0.0
-------------------------------------------------------------------------------
Numerical solver
-------------------------------------------------------------------------------
(CVXPY) Apr 05 01:01:27 PM: Invoking solver CVXOPT to obtain a solution.
pcost dcost gap pres dres k/t
0: 3.1191e+11 -1.7680e+06 7e+57 1e+01 1e-16 1e+00
Traceback (most recent call last):
File "C:\Users\chamo\PycharmProjects\OptomizeModel\main.py", line 12, in <module>
model_cvxpy()
File "C:\Users\chamo\PycharmProjects\OptomizeModel\Model.py", line 140, in model_cvxpy
prob.solve(verbose=True, solver=cp.CVXOPT)
File "C:\Users\chamo\PycharmProjects\OptomizeModel\venv\lib\site-packages\cvxpy\problems\problem.py", line 481, in solve
return solve_func(self, *args, **kwargs)
File "C:\Users\chamo\PycharmProjects\OptomizeModel\venv\lib\site-packages\cvxpy\problems\problem.py", line 1012, in _solve
solution = solving_chain.solve_via_data(
File "C:\Users\chamo\PycharmProjects\OptomizeModel\venv\lib\site-packages\cvxpy\reductions\solvers\solving_chain.py", line 361, in solve_via_data
return self.solver.solve_via_data(data, warm_start, verbose,
File "C:\Users\chamo\PycharmProjects\OptomizeModel\venv\lib\site-packages\cvxpy\reductions\solvers\conic_solvers\cvxopt_conif.py", line 202, in solve_via_data
results_dict = cvxopt.solvers.conelp(c, G, h, dims, A, b,
File "C:\Users\chamo\PycharmProjects\OptomizeModel\venv\lib\site-packages\cvxopt\coneprog.py", line 1033, in conelp
W = misc.compute_scaling(s, z, lmbda, dims, mnl = 0)
File "C:\Users\chamo\PycharmProjects\OptomizeModel\venv\lib\site-packages\cvxopt\misc.py", line 322, in compute_scaling
W['beta'][k] = math.sqrt( aa / bb )
ZeroDivisionError: float division by zero
Process finished with exit code 1
The fact is that if I put a parameter in the function, I will call the minimize method with parameters
res = minimize(f_scipy, x0, method='SLSQP', constraints=[ineq_constraints], bounds=bounds)
That as a result for some reason res.x = x0, I don't understand why. And he's ignoring my conditions for the time fund. Here I need help. Plus, my constraint conditions will then also include indicator functions (there will be added time for equipment changeover if we are going to produce a specific part on the machine, now I just want to run the model so that it gives the expected result).
Among other things, if I run my function as at the very beginning, just by conditions on x_i, then my answer seems to coincide with the original one, or rather, it really matches, except for those places where the penalty (the 3rd term from the target is called that) overlaps the first two terms and everything must be done there.
Moreover, for parts with b_i = 0, x_i should result in q_i after optimization. This is achieved by the desired degree of c_i in the 3rd term in the target for the desired part. But I don't know how to select this degree programmatically.
The solution by the cvxpy library breaks somewhere inside, there is a division by zero, and other methods from this library either do not work, or give infinity in the answer.
I need help, maybe I'm setting the model incorrectly, or is it easier to write the algorithm myself, rather than using ready-made libraries?
How to make scipy.optimize respect inequality constraints?
Please help.
P.S. I am attaching the files data.xlsx
df_t_o_p.xlsx.
