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: 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]).