[Python] SHGO optimisation ending upon None. Tried all three sampling methods: sobol, halton,...

I am minimising the objective function which is -100 * fill rate(floors) together with L2 penalty on epsilon and gamma: -100 * fill_rate(floors) + l1 * ||epsilon|| + l2 * ||gamma|| I have constraints which involves evaluation of revenue function which is computed via sampling(Monte Carlo). Hence the revenue function has variability. Constraint on revenue is: 0.95 * budget <= revenue(floors) <= budget Constraint on gamma also comes from the equation: 0 <= participation rate(floor) + gamma(floor) <= 1 Moreover we have bounds on floors, epsilon and gamma. Each floors, epsilon, gamma are 3 dim array. We want to find the optimum floors having minimal objection function value satisfying constraints and bounds. Epsilon and gamma are slack variables. It may be possible there are multiple optimum floors exists. At least one optimal floor is there is also shown in the code. Im using SHGO for optimisation but output I am getting is None. Same code with optimisation 'COBYLA' and 'BOBYQA' is working. Any idea whether the variability of revenue evaluation is creating problem for SHGO optimisation. Any suggestion to make it work for SHGO are welcome. Thanks in advance!

import time
import numpy as np
from scipy.stats import norm, bernoulli
from shgo import shgo

"""
Supply = BQ1 U BQ2 U BQ3
"""
# means of the bqs
mean1 = 5
mean2 = 6
mean3 = 3
# variance of bqs
std1 = 5
std2 = 4
std3 = 1
# probabilities that bqs wont participate in the auction due to some reasons
err1 = 0.6
err2 = 0.7
err3 = 0.2
# current floors of each bqs
# do the auction and collect data to get ao, fr, cpm


def get_bq_samples(pr, loc, scale, size):
"""
get real samples of bqs when they participate in the auction and collect the data. The real distributions of bqs are
unknown to us. Once we have collected data, we will use that to evaluate our objective function
"""
samples1 = norm.rvs(loc=loc, scale=scale, size=size)
samples2 = bernoulli.rvs(p=pr, size=size)
return np.multiply(samples1, samples2)


def simulate_auction_collect_bids(bq_budgets, floors, num_samples):
"""
Simulate the first price auction for 3 advertisers(bqs) and collect the bids data
"""
X1, X2, X3 = [], [], []
R1, R2, R3 = 0, 0, 0
bids_above_floor1 = 0
bids_above_floor2 = 0
bids_above_floor3 = 0
samples1 = get_bq_samples(err1, mean1, std1, num_samples)
samples2 = get_bq_samples(err2, mean2, std2, num_samples)
samples3 = get_bq_samples(err3, mean3, std3, num_samples)

bq1_budget_exhausted = False
bq2_budget_exhausted = False
bq3_budget_exhausted = False
served_ao = 0
supply_rev = 0
for b1, b2, b3 in zip(samples1, samples2, samples3):
if bq1_budget_exhausted and bq2_budget_exhausted and bq3_budget_exhausted:
print("Budget exhausted but inventory is there")
X1.append(0)
X2.append(0)
X3.append(0)
continue
ab1 = b1
ab2 = b2
ab3 = b3
bids_above_floor1 += 1
bids_above_floor2 += 1
bids_above_floor3 += 1
if not bq1_budget_exhausted:
X1.append(b1)
else:
X1.append(0)
if not bq2_budget_exhausted:
X2.append(b2)
else:
X2.append(0)
if not bq3_budget_exhausted:
X3.append(b3)
else:
X3.append(0)
if (b1 < floors[0]) or (bq1_budget_exhausted == True):
ab1 = 0
bids_above_floor1 -= 1
if (b2 < floors[1]) or (bq2_budget_exhausted == True):
ab2 = 0
bids_above_floor2 -= 1
if (b3 < floors[2]) or (bq3_budget_exhausted == True):
ab3 = 0
bids_above_floor3 -= 1
actual_bids = np.array([ab1, ab2, ab3])
Y = np.amax(actual_bids)
Z = np.argmax(actual_bids)
if Y > 0:
served_ao += 1
supply_rev += Y
if Z == 0:
R1 += Y
elif Z == 1:
R2 += Y
elif Z == 2:
R3 += Y
if R1 >= bq_budgets[0]:
bq1_budget_exhausted = True
if R2 >= bq_budgets[1]:
bq2_budget_exhausted = True
if R3 >= bq_budgets[2]:
bq3_budget_exhausted = True

supply_fr = served_ao / num_samples
bq1_fr = bids_above_floor1 / num_samples
bq2_fr = bids_above_floor2 / num_samples
bq3_fr = bids_above_floor3 / num_samples
return supply_rev, supply_fr, [R1, R2, R3], [bq1_fr, bq2_fr, bq3_fr], np.array([X1, X2, X3])


def get_participation_rate(max_bids, floors):
"""
Given bids data, get the participation rate for all bqs
"""
N = len(max_bids[0])
pr = []
for bq_bids, f in zip(max_bids, floors):
rate = np.count_nonzero(np.where(bq_bids < f, 0, bq_bids))
rate = rate / N
pr.append(rate)
return pr


def get_delta_samples(max_bids, floors, gamma, num_samples):
"""
Given bids data, get the bernoulli samples with the participation rate determined by the bids observed
"""
pr_arr = get_participation_rate(max_bids, floors)
temp = []
for pr, g in zip(pr_arr, gamma):
if 0 <= pr+g <= 1:
temp.append(bernoulli.rvs(p=pr+g, size=num_samples))
else:
print("Optimisation not respecting the gamma constraint!")
t = max(min(pr+g, 1),0)
temp.append(bernoulli.rvs(p=0, size=num_samples))
return np.array(temp)


def get_max_bids_samples(max_bids, floors, epsilon, num_samples):
"""
Get max-bids samples using the observed bids data
"""
temp = []
for bq_bids, f, e in zip(max_bids, floors, epsilon):
valid_bids = bq_bids[bq_bids >= f]
valid_bids = np.add(valid_bids, e)
temp.append(np.random.choice(valid_bids, size=num_samples))
return np.array(temp)


def get_revenue(bids, X, ao, num_samples, truncate=True):
"""
Get the revenue as a function of floor using the observed bids data. We need this into constraint equation
"""
floors = X[:3]
epsilon = X[3:6]
gamma = X[6:9]
max_bid_samples = get_max_bids_samples(bids, floors, epsilon, num_samples)
delta_samples = get_delta_samples(bids, floors, gamma, num_samples)
X = np.multiply(max_bid_samples, delta_samples)
Y = np.amax(X, axis=0)
Z = np.argmax(X, axis=0)
R1 = np.sum(Y[np.argwhere(Z == 0)])
R2 = np.sum(Y[np.argwhere(Z == 1)])
R3 = np.sum(Y[np.argwhere(Z == 2)])
if truncate:
return [round(np.mean(Y), 3) * ao, round(R1 / num_samples, 3) * ao, round(R2 / num_samples, 3) * ao,
round(R3 / num_samples, 3) * ao]
else:
return [np.mean(Y) * ao, R1 * ao / num_samples, R2 * ao / num_samples, R3 * ao / num_samples]


def get_fr(bids, X):
"""
Get fill rate as a function of floor using the observed bids data. fr is part of optimisation objective function
"""
floors = X[:3]
gamma = X[6:9]
prs = get_participation_rate(bids, floors)
fr = 1
for pr, g in zip(prs, gamma):
if 0 <= pr + g <= 1:
fr = fr * (1 - (pr+g))
else:
t = max(min(pr+g, 1),0)
t=0
fr = fr*(1-t)
return 1-fr


def l2_norm(arr):
out = 0
for a in arr:
out += a * a
return out/len(arr)


def find_optimum_floors_shgo(bids, budget, num_samples):
"""
We want to minimise out objective function which is -100 * fr + l1 * ||epsilon|| + l2 * ||gamma||, subject to revenue
constraint: 0.95 * budget <= revenue(floor) <= budget
and 0 <= participation rate + gamma <= 1
together with bounds on floors, epsilon and gamma
"""
# X: input vector for optimisation [floor1, floor2, floor3, epsilon1, epsilon2, epsilon3, gamma1, gamma2, gamma3]
# bounds contains the bounds for X vector
AO = len(bids[0])
bounds = [(2.5, 3.5), (3.5, 4.5), (1.0, 2.0), (-0.15, 0.15), (-0.15, 0.15), (-0.15, 0.15), (-0.1, 0.1), (-0.1, 0.1), (-0.1, 0.1)]

# l1 l2 are hyperparameters for the objective function
l1 = 50
l2 = 500

def objective_function(X):
# objective function is the fill rate minus penalties on epsilon and gamma
ep = X[3:6]
gm = X[6:9]
fr = get_fr(bids, X)
return -100*fr + l1 * l2_norm(ep) + l2 * l2_norm(gm)

# budget constraint: we want the revenue to remain between 0.95 of budget and budget
def get_lower_budget_constraint(X, i):
rev_arr = get_revenue(bids, X, AO, num_samples, False)
return rev_arr[1:] - 0.95 * budget

def get_upper_budget_constraint(X, i):
rev_arr = get_revenue(bids, X, AO, num_samples, False)
return budget - rev_arr[1:]

# participation rate constraint: we want participation rate to remain between 0.1 and 0.8
def get_lower_delta_constraint(X, i):
flr = X[:3]
gm = X[6:9]
pr = np.add(get_participation_rate(bids, flr),gm)
return pr

def get_upper_delta_constraint(X, i):
flr = X[:3]
gm = X[6:9]
pr = np.add(get_participation_rate(bids, flr), gm)
return 1-pr

constraints = []
lower_delta_constraint = [{'type': 'ineq', 'fun': lambda x, i=i: get_lower_delta_constraint(x, i)} for i in
range(3)]
constraints.extend(lower_delta_constraint)
upper_delta_constraint = [{'type': 'ineq', 'fun': lambda x, i=i: get_upper_delta_constraint(x, i)} for i in
range(3)]
constraints.extend(upper_delta_constraint)
lower_budget_constraint = [{'type': 'ineq', 'fun': lambda x, i=i: get_lower_budget_constraint(x,i)} for i in range(3)]
constraints.extend(lower_budget_constraint)
upper_budget_constraint = [{'type': 'ineq', 'fun': lambda x, i=i: get_upper_budget_constraint(x, i)} for i in range(3)]
constraints.extend(upper_budget_constraint)

res = shgo(func = objective_function, bounds = bounds, sampling_method='sobol', constraints=constraints)

print(res.x)
print(res.fun)

if res.x is not None:
print(get_revenue(bids, res.x, AO, num_samples, False))
print(get_fr(bids, res.x))
print(res.x[:3])


if __name__ == '__main__':
ao = 10000
budgets = np.array([25000, 35000, 2000])
initial_floors = np.array([3, 4, 1.5])


# do the simulate and collect the data
s_rev, s_fr, bq_rev, bq_fr, all_bids = simulate_auction_collect_bids(budgets, initial_floors, ao)
print(f"Supply revenue: {s_rev}")
print(f"Supply fill rate: {s_fr}")
print(f"Each BQ revenue: {bq_rev}")
print(f"Each BQ fill rate: {bq_fr}")

# one solution satisfying all constraints
initial_X = np.append(initial_floors, [0, -0.1, -0.1, 0.0, 0.03, 0.018])
print(f"One solution with X: {initial_X}")
revenue_init_flr = get_revenue(all_bids, initial_X, ao, 100000, truncate=False)
print(np.divide(np.array(revenue_init_flr[1:]), budgets) * 100)
fr_init_flr = get_participation_rate(all_bids, initial_floors)
print(fr_init_flr)

#finding optimum via shgo
print("Starting shgo optimisation, to find optimum floors")
find_optimum_floors_shgo(all_bids, budgets, 100000)

Continue reading...