Scipy minimize successfully terminates, but doesn’t fulfill inequality constraints

Question

I'm trying to minimize the function 0.5*(x^2+y^2) subject to a series (N=20) of inequality constraints of the form x a1+y a2+a3 z >= 1. The solution should be around x=0.50251, y=-0.5846, z=0.36787. The routine terminates with the message 'Optimization terminated successfully' but more than half of the constraints are not respected. I tried also different solvers, with the same result.

Scaling the objective function changes the solution, but does not converge to the expected solution.

from scipy.optimize import minimize
import numpy as np


Pct=np.array([[-0.664,  3.179],[ 0.231, -2.044],[-2.493,  3.25 ],[ 0.497, -0.654],[-1.27,   1.248],[-1.185,  1.814],[-1.843,  4.386],[-1.616,  1.401],[ 0.052, -1.232],[-3.145,  0.404],[ 0.672, -1.655],[ 2.202, -1.888],[ 4.084, -1.067],[ 1.006, -1.671],[-2.255,  1.51 ],[-1.264,  1.663],[ 1.897, -2.217],[ 1.843, -1.276],[-1.693,  1.623],[ 2.297, -1.709]])
Sid=np.array([-1,  1, -1,  1, -1, -1, -1, -1,  1, -1,  1,  1,  1,  1, -1, -1,  1,  1, -1,  1])

# func to be minimized 
def OptFunc(x):
  return 0.5*(x[0]**2+x[1]**2)
def JacOptFunc(x):
  return np.array([x[0],x[1],0.0])  

# Constraints
c=[]
for i in range(len(Sid)):
  c+=[{'type': 'ineq', 'fun': lambda x:  Sid[i]*(x[0]*Pct[i,0]+x[1]*Pct[i,1]+x[2])-1 }]
cons=tuple(c)    

# start optimization
res = minimize(OptFunc,(0.3,-0.2,0.1),constraints=cons,method='SLSQP',jac=JacOptFunc)

#expected solution should be around
# [0.5025062702615434, -0.584685257866671, 0.36787016514022236]
print("-->",res.message)
print("solution ",res.x,flush=True)

print("Check Constraints")
cons=list(cons)
for i in range(len(cons)):
  lokfun=c[i]['fun']
  print("Constraint # ",i," value: ",lokfun(res.x))

The expected result is around x=0.50251, y=-0.5846, z=0.36787 but I get the following output:

--> Optimization terminated successfully.
solution  [-1.14580677e-04 -1.16285140e-04  1.00006446e+00]

Check Constraints
Constraint #  0  value:  -1.9997708716077622
Constraint #  1  value:  0.0002756791862408292
Constraint #  2  value:  -1.999972183420499
Constraint #  3  value:  8.356438220613605e-05
Constraint #  4  value:  -2.0000648541023893
Constraint #  5  value:  -1.9999892973558606
Constraint #  6  value:  -1.9997656060620763
Constraint #  7  value:  -2.000086707390163
Constraint #  8  value:  0.00020176559401496874
Constraint #  9  value:  -2.0003778375289833
Constraint #  10  value:  0.00017991418852214558
Constraint #  11  value:  3.1700190727068644e-05
Constraint #  12  value:  -0.0002794107423930159
Constraint #  13  value:  0.00014350480474445426
Constraint #  14  value:  -2.000147249362345
Constraint #  15  value:  -2.0000159082853974
Constraint #  16  value:  0.00010490510804150865
Constraint #  17  value:  1.6681482228886324e-06
Constraint #  18  value:  -2.0000697148012767
Constraint #  19  value:  -1.354516498963676e-11

Answer 1

I know little about scipy.optimize, but I can spot a problem with

for i in range(len(Sid)):
  c+=[{'type': 'ineq', 'fun': lambda x:  Sid[i]*(x[0]*Pct[i,0]+x[1]*Pct[i,1]+x[2])-1 }]

The problem is that Python closures are late-binding, which means that the value of i in each constraint is actually evaluated after the loop completed. In effect, you are actually imposing the same (last) constraint 20 times. See https://docs.python-guide.org/writing/gotchas/#late-binding-closures

A possible fix:

from scipy.optimize import minimize
import numpy as np


Pct=np.array([[-0.664,  3.179],[ 0.231, -2.044],[-2.493,  3.25 ],[ 0.497, -0.654],[-1.27,   1.248],[-1.185,  1.814],[-1.843,  4.386],[-1.616,  1.401],[ 0.052, -1.232],[-3.145,  0.404],[ 0.672, -1.655],[ 2.202, -1.888],[ 4.084, -1.067],[ 1.006, -1.671],[-2.255,  1.51 ],[-1.264,  1.663],[ 1.897, -2.217],[ 1.843, -1.276],[-1.693,  1.623],[ 2.297, -1.709]])
Sid=np.array([-1,  1, -1,  1, -1, -1, -1, -1,  1, -1,  1,  1,  1,  1, -1, -1,  1,  1, -1,  1])

# func to be minimized 
def OptFunc(x):
    return 0.5*(x[0]**2+x[1]**2)
def JacOptFunc(x):
    return np.array([x[0],x[1],0.0])  

# Constraints
def constraint_maker(i=0):  # i MUST be an optional keyword argument, else it will not work
    def constraint(x):
        return Sid[i]*(x[0]*Pct[i,0]+x[1]*Pct[i,1]+x[2])-1 
    return constraint

c=[]
for i in range(len(Sid)):
    c+=[{'type': 'ineq', 'fun': constraint_maker(i)}]
cons=tuple(c)

# start optimization
res = minimize(OptFunc,(0.3,-0.2,0.1),constraints=cons,method='SLSQP',jac=JacOptFunc)

#expected solution should be around
# [0.5025062702615434, -0.584685257866671, 0.36787016514022236]
print("-->",res.message)
print("solution ",res.x)

print("Check Constraints")
cons=list(cons)
for i in range(len(cons)):
    lokfun=c[i]['fun']
    print("Constraint # ",i," value: ",lokfun(res.x))

results in

--> Optimization terminated successfully.
solution  [ 0.52374351 -0.56495542  0.37021863]
Check Constraints
Constraint #  0  value:  0.7735403550593944
Constraint #  1  value:  0.6459722649608017
Constraint #  2  value:  1.7715790719554194
Constraint #  3  value:  8.137268636687622e-11
Constraint #  4  value:  -2.2235047136831554e-10
Constraint #  5  value:  0.27524657110337936
Constraint #  6  value:  2.0729351509689136
Constraint #  7  value:  0.2676534344356165
Constraint #  8  value:  0.09347837249122604
Constraint #  9  value:  0.5051967055706261
Constraint #  10  value:  0.6571754935710583
Constraint #  11  value:  1.5901376792721638
Constraint #  12  value:  2.1119945643862095
Constraint #  13  value:  0.8411451130595076
Constraint #  14  value:  0.6639056792092357
Constraint #  15  value:  0.23131403951409935
Constraint #  16  value:  1.6162662427554526
Constraint #  17  value:  1.0563610395273058
Constraint #  18  value:  0.43340178883510116
Constraint #  19  value:  1.5387662919992176