Проблема со скоростью решения задачи SciPy.Minimaze (SLSQP)

Я пытаюсь написать решатель для производственной задачи, который должен распределять рабочее время по дням и рабочим местам, минимизируя переработку. В водных данных у меня 6000 переменных и 7000 ограничений. Он работает неприлично медленно. Есть ли какие либо способы ускорить код?

import numpy as np
import pandas as pd
from scipy.optimize import minimize
from scipy.optimize import Bounds
from datetime import datetime

BD = pd.read_excel('Загрузка.xlsx', sheet_name='БД')
PZ = pd.read_excel('Загрузка.xlsx', sheet_name='ПЗ')
TIME = pd.read_excel('Загрузка.xlsx', sheet_name='Время работы')

x0 = np.zeros(n)

#Определяем индексы BD для каждой стоки TIME
TIME_ind = []
for n_Price in TIME.index:
    tab_BD_ind = BD[(BD['День'] == TIME.at[n_Price, 'День']) & (BD['Рабочий центр пр-во'] == TIME.at[n_Price, 'Рабочий центр пр-во'])].index.values.tolist()
    TIME_ind.append(tab_BD_ind)

#Целевая функция
def objective(x):
    res = 0
    for n_Price in TIME.index:
        res += sum(x[TIME_ind[n_Price]] * TIME.at[n_Price, 'Цена часа работы']) #Стоимость работы
        res += max(sum(x[TIME_ind[n_Price]]) - TIME.at[n_Price, 'Опт загрузка'],0) * TIME.at[n_Price, 'Цена часа переработки']
    return res

#Определяем индексы для производства и потребления деталей ПЗ
PZ_ind_prvo = []
PZ_ind_potr = []
for n_PZ in PZ.index:
    tab_BD_ind_prvo = BD[(BD['Номенклатура пр-во'] == PZ.at[n_PZ, 'Номенклатура']) & (BD['Характеристика пр-во'] == PZ.at[n_PZ, 'Характеристика']) 
                    & (BD['Рабочий центр пр-во'] == PZ.at[n_PZ, 'Рабочий центр пр-во'])].index.values.tolist()
    tab_BD_ind_potr = BD[(BD['Номенклатура потр-е'] == PZ.at[n_PZ, 'Номенклатура']) & (BD['Характеристика потр-е'] == PZ.at[n_PZ, 'Характеристика']) 
                    & (BD['Рабочий центр потр-е'] == PZ.at[n_PZ, 'Рабочий центр пр-во'])].index.values.tolist()
    PZ_ind_prvo.append(tab_BD_ind_prvo)
    PZ_ind_potr.append(tab_BD_ind_potr)

#Добавляем ограничения по ПЗ
for n_PZ in PZ.index:
    constrants.append({'type' : 'ineq', 'fun' : lambda x, n_PZ = n_PZ: sum(x[PZ_ind_prvo[n_PZ]] * BD.loc[PZ_ind_prvo[n_PZ], 'Объем производства в час'].values.tolist()) 
                       - sum(x[PZ_ind_potr[n_PZ]] * BD.loc[PZ_ind_potr[n_PZ], 'Объем потребления в час'].values.tolist()) 
                      - PZ.loc[n_PZ, 'ПЗ']})

#Определяем индексы для остатков (производство, потребление)
ost_ind_prvo = []
ost_ind_potr = []

tab_nom = BD[['День','Номенклатура пр-во', 'Характеристика пр-во', 'Рабочий центр пр-во']].drop_duplicates().reindex()

for n_nom in tab_nom.index:
    tab_BD_ind_prvo = BD[(BD['День'] <= tab_nom.at[n_nom, 'День']) & (BD['Номенклатура пр-во'] == tab_nom.at[n_nom, 'Номенклатура пр-во']) 
                        & (BD['Характеристика пр-во'] == tab_nom.at[n_nom, 'Характеристика пр-во']) 
                        & (BD['Рабочий центр пр-во'] == tab_nom.at[n_nom, 'Рабочий центр пр-во'])].index.values.tolist()
    tab_BD_ind_potr = BD[(BD['День'] <= tab_nom.at[n_nom, 'День']) & (BD['Номенклатура потр-е'] == tab_nom.at[n_nom, 'Номенклатура пр-во']) 
                        & (BD['Характеристика потр-е'] == tab_nom.at[n_nom, 'Характеристика пр-во']) 
                        & (BD['Рабочий центр потр-е'] == tab_nom.at[n_nom, 'Рабочий центр пр-во'])].index.values.tolist()
    ost_ind_prvo.append(tab_BD_ind_prvo)
    ost_ind_potr.append(tab_BD_ind_potr)

#Добавляем ограничения по остатку
for n_ost in tab_nom.index:
    constrants.append({'type' : 'ineq', 'fun' : lambda x, n_ost = n_ost: sum(x[ost_ind_prvo[n_ost]] * BD.loc[ost_ind_prvo[n_ost], 'Объем производства в час'].values.tolist()) 
                       - sum(x[ost_ind_potr[n_ost]] * BD.loc[ost_ind_potr[n_ost], 'Объем потребления в час'].values.tolist())})

#Определяем индексы для дней и станков
days_ind = []
tab_days = BD[['День', 'Рабочий центр пр-во']].drop_duplicates().reset_index()
for n_tab_day in tab_days.index:
    tab_24_ind = BD[(BD['День'] == tab_days.at[n_tab_day, 'День']) & (BD['Рабочий центр пр-во'] == tab_days.at[n_tab_day, 'Рабочий центр пр-во'])].index.values.tolist()
    days_ind.append(tab_24_ind)

#Добавляем ограничения по 24 часам
for n_days in tab_days.index:
    constrants.append({'type' : 'ineq', 'fun' : lambda x, n_days = n_days: (24 - sum(x[days_ind[n_days]]))})

#CALLBACK
n_iter = 0
def call_back(kx):
    global n_iter
    n_iter += 1
    print('n_iter: {0} time: {1} x: {2} fun: {3}'.format(n_iter , datetime.now(), sum(kx), objective(kx)))

#Решаем
solution = minimize(objective, x0, method='SLSQP', constraints=constrants, options={'maxiter' : 10}, callback = call_back)