Defining dynamic constraints for scipy optimize in Python

Question

I wanted to abstract the following function that calculates minimum value of a objective function and values when we can get this minimal value for arbitrary number of g's.

I started with simple case of two variables, which works fine

import numpy as np
from scipy.optimize import minimize

def optimize(g_0, s, g_max, eff, dist):
    objective = lambda x: s[0] * x[0] + s[1] * x[1]

    cons = [
        {'type': 'ineq', 'fun': lambda x: x[0] - ((dist[0] + dist[1]) / eff - g_0)},                   # g_1 > (dist[0] + dist[1]) / eff - g_0
        {'type': 'ineq', 'fun': lambda x: x[1] - ((dist[0] + dist[1] + dist[2]) / eff - x[0] - g_0)},  # g_2 > (dist[0] + dist[1] + dist[2]) / eff - g_1 - g_0
        {'type': 'ineq', 'fun': lambda x: g_max - (dist[0] / eff - g_0) - x[0]},                       # g_1 < g_max - (dist[0] / eff - g_0)
        {'type': 'ineq', 'fun': lambda x: g_max - (g_0 + x[0] - (dist[0] - dist[1]) / eff) - x[1]},    # g_2 < g_max - (g_0 + g_1 - (dist[0] - dist[1]) / eff)
    ]

    # General constraints for all g
    for i in range(len(s)):
        cons.append({'type': 'ineq', 'fun': lambda x, i=i: x[i]})                   # g_i > 0
        cons.append({'type': 'ineq', 'fun': lambda x, i=i: g_max - x[i]})           # g_i < g_max

    # Bounds for the variables (g_1 and g_2)
    g1_lower_bound = max(0, (dist[0] + dist[1]) / eff - g_0)
    g1_upper_bound = min(g_max, g_max - (dist[0] / eff - g_0))

    # Initial guess for the variables
    x0 = [g1_lower_bound, max(0, ((dist[0] + dist[1] + dist[2]) / eff - g1_lower_bound - g_0) + 1)]

    solution = minimize(objective, x0, method='SLSQP', bounds=[(g1_lower_bound, g1_upper_bound), (0, g_max)], constraints=cons)

    g_1, g_2 = map(round, solution.x)
    return g_1, g_2, round(solution.fun, 2)


g_0 = 80
s = [4.5, 3]
g_max = 135
eff = 5
dist = [400, 500, 600]
optimal_g1, optimal_g2, minimum_value = optimize(g_0, s, g_max, eff, dist)
print(f"Optimal values: g_1 = {optimal_g1}, g_2 = {optimal_g2}")
print(f"Minimum value of the objective function: {minimum_value}")

Then I started to abstract it for any arbitrary numbers of g (i>=1). Here's the generale rule for inequality

So far I came to this code, but when I tried to send the exact same parametrs it gives different result

import numpy as np
from scipy.optimize import minimize

def optimize(g_0, s, g_max, eff, dist):
    objective = lambda x: sum(s[i] * x[i] for i in range(len(x)))

    cons = []
    for i in range(len(s)):
        cons.append({'type': 'ineq', 'fun': lambda x, i=i: x[i] - (sum(dist[:i+2]) / eff - sum([g_0] + x[:i]))})  # g_i > sum of dists from dist[0] to dist[i] / eff - sum of g from g_0 to g_(i-1)

        cons.append({'type': 'ineq', 'fun': lambda x, i=i: g_max - (sum([g_0] + x[:i]) - (sum(dist[:i+1]) / eff)) - x[i]})  # g_i < g_max - (sum of g from g_0 to g_(i-1) - sum of dists from dist[0] to dist[i] / eff)
        
        # General constraints to ensure each g is between 0 and g_max
        cons.append({'type': 'ineq', 'fun': lambda x, i=i: x[i]})                   # g_i > 0
        cons.append({'type': 'ineq', 'fun': lambda x, i=i: g_max - x[i]})           # g_i < g_max

    g1_lower_bound = max(0, (dist[0] + dist[1]) / eff - g_0)
    g1_upper_bound = min(g_max, g_max - (dist[0] / eff - g_0))

    # Initial guess for the variables
    x0 = [g1_lower_bound, max(0, ((dist[0] + dist[1] + dist[2]) / eff - g1_lower_bound - g_0) + 1)]
    solution = minimize(objective, x0, method='SLSQP', bounds=[(0, g_max) for _ in range(len(s))], constraints=cons)

    g_values = list(map(round, solution.x))
    return g_values, round(solution.fun, 2)

g_0 = 80
s = [4.5, 3]
g_max = 135
eff = 5
dist = [400, 500, 600]
optimal_g_values, minimum_value = optimize(g_0, s, g_max, eff, dist)
print(f"Optimal values: g = {optimal_g_values}")
print(f"Minimum value of the objective function: {minimum_value}")

The first function gives following result, which is also correct one:

Optimal values: g_1 = 100, g_2 = 120 Minimum value of the objective function: 810.0

But the second, after trying to make it for arbitrary len of g it returns:

Optimal values: g = [135, 85] Minimum value of the objective function: 862.5

I checked the lambda functions in second function and they seems to be correct. Also when I try to execute this code:

import numpy as np
from scipy.optimize import minimize

def optimize(g_0, s, g_max, eff, dist):
    objective = lambda x: sum(s[i] * x[i] for i in range(len(x)))

    cons = [
        {'type': 'ineq', 'fun': lambda x: x[0] - ((dist[0] + dist[1]) / eff - g_0)},                   # g_1 > (dist[0] + dist[1]) / eff - g_0
        {'type': 'ineq', 'fun': lambda x: x[1] - ((dist[0] + dist[1] + dist[2]) / eff - x[0] - g_0)},  # g_2 > (dist[0] + dist[1] + dist[2]) / eff - g_1 - g_0
    ]
    
    for i in range(len(s)):
        #cons.append({'type': 'ineq', 'fun': lambda x, i=i: x[i] - (sum(dist[:i+2]) / eff - sum([g_0] + x[:i]))})  # g_i > sum of dists from dist[0] to dist[i] / eff - sum of g from g_0 to g_(i-1)

        cons.append({'type': 'ineq', 'fun': lambda x, i=i: g_max - (sum([g_0] + x[:i]) - (sum(dist[:i+1]) / eff)) - x[i]})  # g_i < g_max - (sum of g from g_0 to g_(i-1) - sum of dists from dist[0] to dist[i] / eff)
        
        # General constraints to ensure each g is between 0 and g_max
        cons.append({'type': 'ineq', 'fun': lambda x, i=i: x[i]})                   # g_i > 0
        cons.append({'type': 'ineq', 'fun': lambda x, i=i: g_max - x[i]})           # g_i < g_max

    g1_lower_bound = max(0, (dist[0] + dist[1]) / eff - g_0)
    g1_upper_bound = min(g_max, g_max - (dist[0] / eff - g_0))

    # Initial guess for the variables
    x0 = [g1_lower_bound, max(0, ((dist[0] + dist[1] + dist[2]) / eff - g1_lower_bound - g_0) + 1)]
    solution = minimize(objective, x0, method='SLSQP', bounds=[(0, g_max) for _ in range(len(s))], constraints=cons)

    g_values = list(map(round, solution.x))
    return g_values, round(solution.fun, 2)

g_0 = 80
s = [4.5, 3]
g_max = 135
eff = 5
dist = [400, 500, 600]
optimal_g_values, minimum_value = optimize(g_0, s, g_max, eff, dist)
print(f"Optimal values: g = {optimal_g_values}")
print(f"Minimum value of the objective function: {minimum_value}")

It gives the correct answer:

Optimal values: g_1 = 100, g_2 = 120 Minimum value of the objective function: 810.0

But if we apply for the other constrain:

import numpy as np
from scipy.optimize import minimize

def optimize(g_0, s, g_max, eff, dist):
    objective = lambda x: sum(s[i] * x[i] for i in range(len(x)))

    cons = [
        {'type': 'ineq', 'fun': lambda x: g_max - (dist[0] / eff - g_0) - x[0]},                       # g_1 < g_max - (dist[0] / eff - g_0)
        {'type': 'ineq', 'fun': lambda x: g_max - (g_0 + x[0] - (dist[0] - dist[1]) / eff) - x[1]},    # g_2 < g_max - (g_0 + g_1 - (dist[0] - dist[1]) / eff)
    ]
    
    for i in range(len(s)):
        cons.append({'type': 'ineq', 'fun': lambda x, i=i: x[i] - (sum(dist[:i+2]) / eff - sum([g_0] + x[:i]))})  # g_i > sum of dists from dist[0] to dist[i] / eff - sum of g from g_0 to g_(i-1)

        #cons.append({'type': 'ineq', 'fun': lambda x, i=i: g_max - (sum([g_0] + x[:i]) - (sum(dist[:i+1]) / eff)) - x[i]})  # g_i < g_max - (sum of g from g_0 to g_(i-1) - sum of dists from dist[0] to dist[i] / eff)
        
        # General constraints to ensure each g is between 0 and g_max
        cons.append({'type': 'ineq', 'fun': lambda x, i=i: x[i]})                   # g_i > 0
        cons.append({'type': 'ineq', 'fun': lambda x, i=i: g_max - x[i]})           # g_i < g_max

    g1_lower_bound = max(0, (dist[0] + dist[1]) / eff - g_0)
    g1_upper_bound = min(g_max, g_max - (dist[0] / eff - g_0))

    # Initial guess for the variables
    x0 = [g1_lower_bound, max(0, ((dist[0] + dist[1] + dist[2]) / eff - g1_lower_bound - g_0) + 1)]
    solution = minimize(objective, x0, method='SLSQP', bounds=[(0, g_max) for _ in range(len(s))], constraints=cons)

    g_values = list(map(round, solution.x))
    return g_values, round(solution.fun, 2)

g_0 = 80
s = [4.5, 3]
g_max = 135
eff = 5
dist = [400, 500, 600]
optimal_g_values, minimum_value = optimize(g_0, s, g_max, eff, dist)
print(f"Optimal values: g = {optimal_g_values}")
print(f"Minimum value of the objective function: {minimum_value}")

It gives the correct values, but wrong minimal value of objective function:

Optimal values: g = [100, 120] Minimum value of the objective function: 810.65

At this point my best guess is that either I have written wrong lambda function in second function or that something unexpected happening in loop

Answer 1

I see two mistakes. One is a mistake in the way you originally specified the two-variable constraint. The second is a mistake in the way the N-variable constraint is specified.

I also see some opportunities for general improvements.

Let's start with the mistake in the original specification:

        {'type': 'ineq', 'fun': lambda x: g_max - (g_0 + x[0] - (dist[0] - dist[1]) / eff) - x[1]},    # g_2 < g_max - (g_0 + g_1 - (dist[0] - dist[1]) / eff)

Based on the inequality you are trying to implement, the expression (dist[0] - dist[1]) ought to be (dist[0] + dist[1]).

Second, I see a mistake relating to the way that NumPy implements +, in the following two lines of code:

        cons.append({'type': 'ineq', 'fun': lambda x, i=i: x[i] - (sum(dist[:i+2]) / eff - sum([g_0] + x[:i]))})  # g_i > sum of dists from dist[0] to dist[i] / eff - sum of g from g_0 to g_(i-1)

        cons.append({'type': 'ineq', 'fun': lambda x, i=i: g_max - (sum([g_0] + x[:i]) - (sum(dist[:i+1]) / eff)) - x[i]})  # g_i < g_max - (sum of g from g_0 to g_(i-1) - sum of dists from dist[0] to dist[i] / eff)

To explain why, let me start with some simple examples.

If you use +, for Python lists this means "concatenate," or put the second list at the end of the first list.

For example:

>>> [1] + [2, 3]
[1, 2, 3]
>>> [1] + []
[1]
>>> [1, 2] + [3, 4]
[1, 2, 3, 4]

However, in NumPy, + means add. If either of the operands to + is a NumPy array, then the two arrays are added together.

For example:

>>> np.array([1]) + [2, 3]
array([3, 4])
>>> np.array([1]) + []
array([], dtype=float64)
>>> np.array([1, 2]) + [3, 4]
array([4, 6])

These are very different results. The effect of that is that the expression sum([g_0] + x[:i]) does different things depending on whether x is a list or array. When SciPy is optimizing your function, x will always be a NumPy array.

For that reason, I suggest that you replace sum([g_0] + x[:i]) with (g_0 + sum(x[:i])), which provides the same results for both lists and arrays.

Also, these constraints are redundant:

        cons.append({'type': 'ineq', 'fun': lambda x, i=i: x[i]})                   # g_i > 0
        cons.append({'type': 'ineq', 'fun': lambda x, i=i: g_max - x[i]})           # g_i < g_max

These do the same thing as your bounds, and bounds are more efficient than constraints. I would remove these.

Also, given that all of your constraints are a linear function of x and some constants, you might find that scipy.optimize.LinearConstraint is more appropriate. Your constraints (except the redundant ones) are equivalent to the following LinearConstraint:

    A = np.zeros((len(s), len(s)))
    cons_lb = np.zeros(len(s))
    cons_ub = np.zeros(len(s))
    for i in range(len(s)):
        A[i, :i + 1] = 1
        cons_lb[i] = sum(dist[:i+2]) / eff - g_0
        cons_ub[i] = g_max - g_0 + (sum(dist[:i+1]) / eff)
    cons = LinearConstraint(A, cons_lb, cons_ub)

This gives a number of benefits. (SciPy does not need numeric differentiation with LinearConstraint, and evaluating a matrix multiply is much faster than evaluating a number of Python functions.)

Speaking of differentiation, you can also speed this up by providing a jacobian, and using NumPy to calculate your objective function. For 2 variables this does not matter much, but for N variables, numeric differentiation takes N calls to your objective function, so it's best to avoid it for large N.

Here is the final code after improving it:

import numpy as np
from scipy.optimize import minimize, LinearConstraint

def optimize(g_0, s, g_max, eff, dist):
    s = np.array(s)
    objective = lambda x: np.sum(s * x)
    jac = lambda x: s

    A = np.zeros((len(s), len(s)))
    cons_lb = np.zeros(len(s))
    cons_ub = np.zeros(len(s))
    for i in range(len(s)):
        A[i, :i + 1] = 1
        cons_lb[i] = sum(dist[:i+2]) / eff - g_0
        cons_ub[i] = g_max - g_0 + (sum(dist[:i+1]) / eff)
    cons = LinearConstraint(A, cons_lb, cons_ub)

    g1_lower_bound = max(0, (dist[0] + dist[1]) / eff - g_0)
    g1_upper_bound = min(g_max, g_max - (dist[0] / eff - g_0))

    # Initial guess for the variables
    x0 = [g1_lower_bound, max(0, ((dist[0] + dist[1] + dist[2]) / eff - g1_lower_bound - g_0) + 1)]
    solution = minimize(
        objective,
        x0,
        jac=jac,
        method='SLSQP',
        bounds=[(0, g_max) for _ in range(len(s))],
        constraints=cons
    )

    g_values = list(map(round, solution.x))
    return g_values, round(solution.fun, 2)

This is about 3x faster, and fixes the bug with the expression sum([g_0] + x[:i]).