I am facing an issue when scaling up my optimisation model. Here is a minimal example, which results in the stated error, when given enough entities to process. You can adjust the number of entities to be processed at the top of the code snippet.
Ideally, this script should run for tens of thousands of entities.
I appreciate the help! Thank you a lot :-)
import numpy as np
import pandas as pd
import datetime
from tqdm import tqdm
from scipy.stats import gaussian_kde
import math
import os
import matplotlib.pyplot as plt
from scipy.stats import norm, poisson
from gekko import GEKKO
NUM_ENTITIES = 50 # 10 crashes already... just a few work out fine
np.random.seed(321)
di_kde = dict()
for entity in range(NUM_ENTITIES):
arr = np.random.normal(14.75, 6, 500)
kernel = gaussian_kde(arr, bw_method='scott')
di_kde[entity] = (kernel, int(np.amax(arr)), int(np.amin(arr)))
def normcdf_gekko(yi, mean, sd):
y = 0.5*(1 + m.erf((yi-mean)/(sd*m.sqrt(2.0))))
z = m.if3(1-y, 1, y)
return z
def cdf_gauss_nogecko(kdes, maxi):
bw = kdes.neff**(-1./(1+4))
res = np.divide(np.cumsum([sum([norm.pdf(i, val, bw) for val in kdes.dataset.flatten()]) for i in np.linspace(0, maxi, 100, endpoint=True)]), max(np.cumsum([sum([norm.pdf(i, val, bw) for val in kdes.dataset.flatten()]) for i in np.linspace(0, maxi, 100, endpoint=True)])))
return res
def cdf_gauss_singular_gecko(kdes, yi):
bw = kdes.neff**(-1./(1+4))
res = m.sum([normcdf_gekko(yi, val, bw) for val in kdes.dataset.flatten()]) / kdes.dataset.flatten().size
return res
def cdf_gauss_singular_nogecko(kdes, yi):
bw = kdes.neff**(-1./(1+4))
res = np.sum([norm.cdf(yi, val, bw) for val in kdes.dataset.flatten()]) / kdes.dataset.flatten().size
return res
def pdf_gauss_gekko(kdes, yi):
bw = kdes.neff**(-1./(1+4))
res_gecko = m.sum([normpdf_gekko(yi, val, bw) for val in kdes.dataset.flatten()])
return res_gecko
def normpdf_gekko(x, mean, sd):
var = float(sd)**2
denom = m.sqrt((2*math.pi*var))
num = m.exp(-(x-float(mean))**2/(2*var))
return num/denom
cost = 0.1
revenue = 1
print(f'{datetime.datetime.now().strftime("%H:%M:%S ")}Running optimisation under constraint...')
if not os.path.isdir(os.path.abspath(r'.\Logging')):
os.mkdir(os.path.abspath(r'.\Logging'))
m = GEKKO(remote=True) # Initialize gekko
#m.time = np.linspace(0,20,100)
m._path = os.path.abspath(r'.\Logging')
m.options.SOLVER = 3 # APOPT (1) is an MINLP solver
m.options.IMODE = 3 # 3 or 6
m.options.MAX_MEMORY = 10
m.solver_options = ['minlp_maximum_iterations 500', \
# minlp iterations with integer solution
'minlp_max_iter_with_int_sol 10', \
# treat minlp as nlp
'minlp_as_nlp 0', \
# nlp sub-problem max iterations
'nlp_maximum_iterations 50', \
# 1 = depth first, 2 = breadth first
'minlp_branch_method 2', \
# maximum deviation from whole number
'minlp_integer_tol 0.05', \
# covergence tolerance
'minlp_gap_tol 0.01']
print(f'{datetime.datetime.now().strftime("%H:%M:%S ")}Adding variables...')
qi = [m.FV(value=di_kde[k][1] , lb=di_kde[k][2], ub=di_kde[k][1], integer=True, name=f'qi_{k}') for idx, k in enumerate(di_kde.keys())]
print(f'{datetime.datetime.now().strftime("%H:%M:%S ")}Adding constraint equations...')
inter = m.Intermediate(m.sum([cdf_gauss_singular_gecko(di_kde[k][0], qi[idx4]) for idx4, k in enumerate(di_kde.keys())]) / len(di_kde), name='InterConstr')
m.Equation(inter >= 0.9)
print(f'{datetime.datetime.now().strftime("%H:%M:%S ")}Adding maximisation function...')
for ii, kk in enumerate(di_kde.keys()):
m.Maximize(m.sum([(-cost * (m.max3(0, qi[ii]-j) + m.max3(0, j-qi[ii]) + qi[ii]) + revenue * m.min3(qi[ii], j)) * pdf_gauss_gekko(di_kde[kk][0], j) for j in range(di_kde[kk][2], di_kde[kk][1])]))
print(f'{datetime.datetime.now().strftime("%H:%M:%S ")}Start solving...')
m.open_folder()
m.options.DIAGLEVEL=10
output = m.solve(disp=True)
print(f'{datetime.datetime.now().strftime("%H:%M:%S ")}End solving...')
print('Results')
print('qi: ' + str(qi))
print('constraint: ' + str(inter.value))
print('Objective: ' + str(-m.options.objfcnval))
