Scipy.optimize - minimize not respecting constraints

Question

Using the code below to to understand how Scipy optmization/minimization works. The results are not matching what I am expecting.

"""
Minimize: f = 2*x[0]*x[1] + 2*x[0] - x[0]**2 - 2*x[1]**2

Subject to: -2*x[0] + 2*x[1] <= -2
             2*x[0] - 4*x[1] <= 0
               x[0]**3 -x[1] == 0

where: 0 <= x[0] <= inf
       1 <= x[1] <= inf
"""
import numpy as np
from scipy.optimize import minimize

def objective(x):
    return 2.0*x[0]*x[1] + 2.0*x[0] - x[0]**2 - 2.0*x[1]**2

def constraint1(x):
    return +2.0*x[0] - 2.0*x[1] - 2.0

def constraint2(x):
    return -2.0*x[0] + 4.0*x[1] 

def constraint3(x):
    sum_eq = x[0]**3.0 -x[1]
    return sum_eq

# initial guesses
n = 2
x0 = np.zeros(n)
x0[0] =  10.0
x0[1] = 100.0

# show initial objective
print('Initial SSE Objective: ' + str(objective(x0)))

# optimize
#b = (1.0,None)
bnds = ((0.0,1000.0), (1.0,1000.0))
con1 = {'type': 'ineq', 'fun': constraint1} 
con2 = {'type': 'ineq', 'fun': constraint2} 
con3 = {'type': 'eq', 'fun': constraint3}
cons = ([con1, con2, con3])
solution = minimize(objective,
                    x0,
                    method='SLSQP',
                    bounds=bnds,
                    constraints=cons)
x = solution.x

print(solution)

# show final objective
print('Final SSE Objective: ' + str(objective(x)))

# print solution
print('Solution')
print('x1 = ' + str(x[0]))
print('x2 = ' + str(x[1]))
print('\n')
print('x', x)
print('constraint1', constraint1(x))
print('constraint2', constraint2(x))
print('constraint3', constraint3(x))

When I run, this is what Python throws on its output console:

Initial SSE Objective: -18080.0
     fun: 2.0
     jac: array([ 0.00000000e+00, -2.98023224e-08])
 message: 'Optimization terminated successfully.'
    nfev: 122
     nit: 17
    njev: 13
  status: 0
 success: True
       x: array([2., 1.])
Final SSE Objective: 2.0

Solution
x1 = 2.0000000000010196
x2 = 1.0000000000012386


x [2. 1.]
constraint1 -4.3787196091216174e-13
constraint2  2.915001573455811e-12
constraint3 7.000000000010997

Despite the optimizer says the result was successful, the constraint3 is not respected because the result should be zero. What am I missing?

Answer 1

Your problem is incompatible. You can eliminate the 3rd constraint (which makes your problem simpler in the first place - only a scalar optimization), after this it is a bit more clear to see what is the problem. From constraint 3 and the lower bound on the original x1 follows, that x0 is not feasible from 0 to 1, so the lower bound in the 1D problem should be 1. It is easy to see that constraint 2 will be always positive, when x0 is larger than 1, therefore it will never be satisfied.

When I run your original problem for me it stops with positive directional derivative (and for the rewritten problem with 'Inequality constraints incompatible'). Which SciPy are you using? For me it is 1.4.1.

On the picture below you can see the objective and the remaining constraints for the 1D problem (horizontal axis is the original x0 variable)

""" Minimize: f = 2*x[0]*x1 + 2*x[0] - x[0]**2 - 2*x1**2

Subject to: -2*x[0] + 2*x[1] <= -2
             2*x[0] - 4*x[1] <= 0
               x[0]**3 -x[1] == 0

where: 0 <= x[0] <= inf
       1 <= x[1] <= inf
"""
import numpy as np
from scipy.optimize import minimize

def objective(x):
    return 2*x**4 + 2*x - x**2 - 2*x**6

def constraint1(x):
    return x - x**3 - 1

def constraint2(x):
    return 2 * x**3 - x
#
# def constraint3(x):
#     sum_eq = x[0]**3.0 -x[1]
#     return sum_eq

# initial guesses
n = 1
x0 = np.zeros(n)
x0[0] =  2.
# x0[1] = 100.0

# show initial objective
print('Initial SSE Objective: ' + str(objective(x0)))

# optimize
#b = (1.0,None)
bnds = ((1.0,1000.0),)
con1 = {'type': 'ineq', 'fun': constraint1}
con2 = {'type': 'ineq', 'fun': constraint2}
# con3 = {'type': 'eq', 'fun': constraint3}
cons = [
    # con1,
    con2,
    # con3,
        ]
solution = minimize(objective,
                    x0,
                    method='SLSQP',
                    bounds=bnds,
                    constraints=cons)
x = solution.x

print(solution)

# show final objective
print('Final SSE Objective: ' + str(objective(x)))

# print solution
print('Solution')
print('x1 = ' + str(x[0]))
# print('x2 = ' + str(x[1]))
print('\n')
print('x', x)
print('constraint1', constraint1(x))
print('constraint2', constraint2(x))
# print('constraint3', constraint3(x))

x_a = np.linspace(1, 2, 200)
f = objective(x_a)
c1 = constraint1(x_a)
c2 = constraint2(x_a)

import matplotlib.pyplot as plt

plt.figure()

plt.plot(x_a, f, label="f")
plt.plot(x_a, c1, label="c1")
plt.plot(x_a, c2, label="c2")
plt.legend()
plt.show()