Constrained optimisation in scipy enters restricted area

Question

I am trying to solve multivariate optimisation problem using Python with scipy.

Let me define environment I am working in:

Searched parameters:

and the problem itself:

(In my case logL function is complex, so I will substitute it with the trivial one, generating similar issue. Therefore in this example I am not using function parameters fully, but I am including those, for problem consistency).

I am using following convention on storing parameters in single, flat array:

Here is the script, that was supposed to solve my problem.

import numpy as np
from scipy import optimize as opt
from pprint import pprint
from typing import List

_d = 2
_tmax = 500.0
_T = [[1,2,3,4,5], [6,7,8,9]]

def logL(args: List[float], T : List[List[float]], tmax : float):
    # simplified - normaly using T in computation, here only to determine dimension
    d = len(T)
    # trivially forcing args to go 'out-of constrains'
    return -sum([(args[2 * i] + args[2 * i + 1] * tmax)**2 for i in range(d)]) 

def gradientForIthDimension(i, d, t_max):
    g = np.zeros(2 * d + 2 * d**2)
    g[2 * i] = 1.0
    g[2 * i + 1] = t_max + 1.0

    return g

def zerosWithOneOnJth(j, l):
    r = [0.0 for _ in range(l)]
    r[j] = 1.0
    
    return r

new_lin_const = {
    'type': 'ineq',
    'fun' : lambda x: np.array(
        [x[2 * i] + x[2 * i + 1] * (_tmax + 1.0) for i in range(_d)] 
        + [x[j] for j in range(2*_d + 2*_d**2) if j not in [2 * i + 1 for i in range(_d)]]
        ),
    'jac' : lambda x: np.array(
        [gradientForIthDimension(i, _d, _tmax) for i in range(_d)]
         + [zerosWithOneOnJth(j, 2*_d + 2*_d**2) for j in range(2*_d + 2*_d**2) if j not in [2 * i + 1 for i in range(_d)]]
        )
}

and finally optimisation

logArgs = [2 for _ in range(2 * (_d ** 2) + 2 * _d)]

# additional bounds, not mentioned in a problem, but suppose a'priori knowledge
bds = [(0.0, 10.0) for _ in range(2 * (_d ** 2) + 2 * _d)]
for i in range(_d):
    bds[2*i + 1] = (-10.0, 10.0)

res = opt.minimize(lambda x, args: -logL(x, args[0], args[1]),
constraints=new_lin_const, x0 = logArgs, args=([_T, _tmax]), method='SLSQP', options={'disp': True}, bounds=bds)

But when checking for result, i am getting:

pprint(res)
#     fun: 2.2124712864600578e-05
#     jac: array([0.00665204, 3.32973738, 0.00665204, 3.32973738, 0.        ,
#       0.        , 0.        , 0.        , 0.        , 0.        ,
#       0.        , 0.        ])
# message: 'Optimization terminated successfully'
#    nfev: 40
#     nit: 3
#    njev: 3
#  status: 0
# success: True
#       x: array([ 1.66633206, -0.00332601,  1.66633206, -0.00332601,  2.        ,
#        2.        ,  2.        ,  2.        ,  2.        ,  2.        ,
#        2.        ,  2.        ])

particullary:

print(res.x[0] + res.x[1]*(501.0))
# -3.2529534621517087e-13

so result is out of constrained area... I was trying to follow documentation, but for me it does not work. I will be happy to hear any advice on what is wrong.

Answer 1

First of all, please stop posting the same question multiple times. This question is basically the same as your other one here. Next time, just edit your question instead of posting a new one.

That being said, your code is needlessly complicated given that your optimization problem is quite simple. It should be your goal that reading your code is as simple as reading the mathematical optimization problem. A more than welcome side effect is that it's much easier to debug your code then in case it's not working as expected.

For this purpose, it's highly recommended that you make yourself familiar with numpy and its vectorized operations (as already mentioned in the comments of your previous question). For instance, you don't need loops to implement your objective, the constraint function or the jacobian. Packing all the optimization variables into one large vector x is the right approach. However, you can simply unpack x into its components lambda, gamma, alpha and beta again. This makes it easier for you to write your functions and easier to read, too.

Well, instead of cutting my way through your code, you can find a simplified and working implementation below. By evaluating the functions and comparing the outputs to the evaluated functions in your code snippet, you should get an idea of what's going wrong on your side.

Edit: It seems like most of the algorithms under the hood of scipy.minimize fail to converge to a local minimizer while preserving strict feasibility of the constraints. If you're open to using another package, I'd recommend using the state-of-the-art NLP solver Ipopt. You can use it by means of the cyipopt package and thanks to its minimize_ipopt method, you can use it similar to scipy.optimize.minimize:

import numpy as np
#from scipy.optimize import minimize
from cyipopt import minimize_ipopt as minimize

d = 2
tmax = 500.0
N = 2*d + 2*d**2

def logL(x, d, tmax):
    lambda_, gamma, alpha, beta = np.split(x, np.cumsum([d, d, d**2]))
    return np.sum((lambda_ + tmax*gamma)**2)

def con_fun(x, d, tmax):
    # split the packed variable x = (lambda_, gamma, alpha, beta)
    lambda_, gamma, alpha, beta = np.split(x, np.cumsum([d, d, d**2]))
    return lambda_ + (tmax + 1.0) * gamma 

def con_jac(x, d, tmax):
    jac = np.block([np.eye(d), (tmax + 1.0)*np.eye(d), np.zeros((d, 2*d**2))])
    return jac

constr = {
    'type': 'ineq', 
    'fun': lambda x: con_fun(x, d, tmax), 
    'jac': lambda x: con_jac(x, d, tmax)
}

bounds = [(0, 10.0)]*N + [(-10.0, 10.0)]*N + [(0.0, 10.0)]*2*d**2
x0 = np.full(N, 2.0)

res = minimize(lambda x: logL(x, d, tmax), x0=x0, constraints=constr, 
    method='SLSQP', options={'disp': True}, bounds=bounds)

print(res)

yields

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

     fun: 0.00014085582293562834
    info: {'x': array([ 2.0037865 ,  2.0037865 , -0.00399079, -0.00399079,  2.00700641,
        2.00700641,  2.00700641,  2.00700641,  2.00700641,  2.00700641,
        2.00700641,  2.00700641]), 'g': array([0.00440135, 0.00440135]), 'obj_val': 0.00014085582293562834, 'mult_g': array([-0.01675576, -0.01675576]), 'mult_x_L': array([5.00053270e-08, 5.00053270e-08, 1.00240003e-08, 1.00240003e-08,
       4.99251018e-08, 4.99251018e-08, 4.99251018e-08, 4.99251018e-08,
       4.99251018e-08, 4.99251018e-08, 4.99251018e-08, 4.99251018e-08]), 'mult_x_U': array([1.25309309e-08, 1.25309309e-08, 1.00160027e-08, 1.00160027e-08,
       1.25359789e-08, 1.25359789e-08, 1.25359789e-08, 1.25359789e-08,
       1.25359789e-08, 1.25359789e-08, 1.25359789e-08, 1.25359789e-08]), 'status': 0, 'status_msg': b'Algorithm terminated successfully at a locally optimal point, satisfying the convergence tolerances (can be specified by options).'}
 message: b'Algorithm terminated successfully at a locally optimal point, satisfying the convergence tolerances (can be specified by options).'
    nfev: 15
     nit: 14
    njev: 16
  status: 0
 success: True
       x: array([ 2.0037865 ,  2.0037865 , -0.00399079, -0.00399079,  2.00700641,
        2.00700641,  2.00700641,  2.00700641,  2.00700641,  2.00700641,
        2.00700641,  2.00700641])

and evaluating the constraint function at the found solution yields

In [17]: print(constr['fun'](res.x))
[0.00440135 0.00440135]

Consequently, the constraints are fulfilled.