Scipy optimize minimize always returns initial guess (SLSQP)

Question

Like the title explains, my program always returns the initial guess.

For context, the program is trying to find the best way to allocate some product across multiple stores. Each stores has a forecast of what they are expected to sell in the following days (sales_data). This forecast does not necessarily have to be integers, or above 1 (it rarely is), it is an expectation in the statistical sense. So, if a store has sales_data = [0.33, 0.33, 0.33] the it is expected that after 3 days, they would sell 1 unit of product.

I want to minimize the total time it takes to sell the units i am allocating (i want to sell them the fastest) and my constraints are that I have to allocate the units that I have available, and I cannot allocate a negative number of product to a store. I am ok having non-integer allocations for now. For my initial allocations i am dividing the units I have available equally among all stores.

Below is a shorter version of my code where I am having the problem:

import numpy, random
from scipy.optimize import curve_fit, minimize

unitsAvailable = 50
days = 15

class Store:

    def __init__(self, num):
        self.num = num

        self.sales_data = []


stores = []
for i in range(10):
    # Identifier
    stores.append(Store(random.randint(1000, 9999)))
    # Expected units to be sold that day (It's unlikey they will sell 1 every day)
    stores[i].sales_data = [random.randint(0, 100) / 100 for i in range(days)]
    print(stores[i].sales_data)


def days_to_turn(alloc, store):
    day = 0

    inventory = alloc
    while (inventory > 0 and day < days):
        inventory -= store.sales_data[day]
        day += 1
    return day

def time_objective(allocations):
    time = 0
    for i in range(len(stores)):
        time += days_to_turn(allocations[i], stores[i])
    return time

def constraint1(allocations):

    return unitsAvailable - sum(allocations)

def constraint2(allocations):

    return min(allocations) - 1

cons = [{'type':'eq', 'fun':constraint1}, {'type':'ineq', 'fun':constraint2}]
guess_allocs = []

for i in range(len(stores)):
    guess_allocs.append(unitsAvailable / len(stores))

guess_allocs = numpy.array(guess_allocs)

print('Optimizing...')

time_solution = minimize(time_objective, guess_allocs, method='SLSQP', constraints=cons, options={'disp':True, 'maxiter': 500})

time_allocationsOpt = [max([a, 0]) for a in time_solution.x]

unitsUsedOpt = sum(time_allocationsOpt)
unitsDaysProjected = time_solution.fun

for i in range(len(stores)):
    print("----------------------------------")
    print("Units to send to Store %s: %s" % (stores[i].num, time_allocationsOpt[i]))
    print("Time to turn allocated: %d" % (days_to_turn(time_allocationsOpt[i], stores[i])))

print("----------------------------------")
print("Estimated days to be sold: " + str(unitsDaysProjected))
print("----------------------------------")
print("Total units sent: " + str(unitsUsedOpt))
print("----------------------------------")

The optimization finishes successfully, with only 1 iteration, and no matter how i change the parameters, it always returns the initial guess_allocs.

Any advice?

Answer 1

The objective function does not have a gradient because it returns discrete multiples of days. This is easily visualized:

import numpy as np
import matplotlib.pyplot as plt

y = []
x = np.linspace(-4, 4, 1000)
for i in x:
    a = guess_allocs + [i, -i, 0, 0, 0, 0, 0, 0, 0, 0]    
    y.append(time_objective(a))

plt.plot(x, y)
plt.xlabel('relative allocation')
plt.ylabel('objective')
plt.show()

If you want to optimize such a function you cannot use gradient based optimizers. There are two options: 1) Find a way to make the objective function differentiable. 2) Use a different optimizer. The first is hard. For the second, let's try dual annealing. Unfortunately, it does not allow constraints so we need to modify the objective function.

Constraining N numbers to a constant sum is the same as having N-1 unconstrained numbers and setting the Nth number to constant - sum.

import scipy.optimize as spo

bounds = [(0, unitsAvailable)] * (len(stores) - 1)

def constrained_objective(partial_allocs):
    if np.sum(partial_allocs) > unitsAvailable:
        # can't sell more than is available, so make the objective infeasible
        return np.inf
    # Partial_alloc contains allocations to all but one store.
    # The final store gets allocated the remaining units.
    allocs = np.append(partial_allocs, unitsAvailable - np.sum(partial_allocs))
    return time_objective(allocs)

time_solution = spo.dual_annealing(constrained_objective, bounds, x0=guess_allocs[:-1])
print(time_solution)

This is a stochastic optimization method. You may want to run it multiple times to see if it can do better, or play with the optional parameters...

Finally, I think there is a problem with the objective function:

for i in range(len(stores)):
    time += days_to_turn(allocations[i], stores[i])

This says that the stores do not sell at the same time but only one after another. Does each store wait with selling until the previous store runs out of items? I think not. Instead, they will sell simultaneously and the time it takes for all units to be sold is the time of the store that takes longest. Try this instead:

for i in range(len(stores)):
    time = max(time, days_to_turn(allocations[i], stores[i]))