import numpy as np
from scipy.optimize import minimize
from scipy.optimize import Bounds

bnds = Bounds([0, 0, 0, 0, 0, 85, 90, 3, 1.2, 145], [2000, 16000, 120, 5000, 2000, 93, 95, 12, 4, 162])  # Задаю нижние и верхние ограничения из таблицы статьи.

eq_cons = {'type': 'eq',
           'fun': lambda x: np.array([1.22 * x[3] - x[0] - x[4], ((98000 * x[2]) / x[3] * x[8] + 1000 * x[2]) - x[5], ((x[1] + x[4]) / x[0]) - x[7]]),
           'jac': lambda x: np.array([[1.22, -1, -1], [98000, 1, 1, 1000, -1], [1, 1, 1, -1]])
           }
ineq_cons = {'type': 'ineq',
            'fun': lambda x: np.array([x[0] * (1.12 + 0.13167 * x[7] - 0.0067 * x[7] ** 2) - 0.99 * x[3],
                                       -(x[0] * (1.12 + 0.13167 * x[7] - 0.0067 * x[7] ** 2) - (100 / 99) * x[3]),
                                       (86.35 + 1.098 * x[7] - 0.038 * x[7] ** 2 + 0.325 * (x[5] - 89)) - 0.99 * x[6],
                                       -(86.35 + 1.098 * x[7] - 0.038 * x[7] ** 2 + 0.325 * (x[5] - 89)) + (100 / 99) *
                                       x[6],
                                       (35.82 - 0.222 * x[9]) - 0.9 * x[8],
                                       -(35.82 - 0.222 * x[9]) - (10 / 9) * x[8],
                                       (-133 + 3 * x[6]) - 0.99 * x[9],
                                        -(133 + 3 * x[6]) + (100 / 99) * x[9]]),
             'jac': lambda x: np.array([[1, 0.13167, -0.0134, -0.99],
                                        [-1, -0.13167, 0.0134, -(100 / 99)],
                                        [1.098, -0.076, 0.325, -0.99],
                                        [-1.098, 0.076, -0.325, + (100 / 99)],
                                        [-0.222, -0.9],
                                        [0.222, -(10 / 9)],
                                        [3, -0.99],
                                        [-3, (100 / 99)]])
              }

x0 = np.array([2000, 0])


def f(x):
    return 0.063 * x[3] * x[6] - 5.04 * x[0] - 0.035 * x[1] - 10 * x[2] - 3.36 * x[4]


res = minimize(f, x0, method='SLSQP', constraints=[eq_cons, ineq_cons], bounds=bnds, options={'ftol': 1e-9, 'disp': True})
print(res.x)