Problems minimizing a function with 2 constraints - Python

Question

I'm writing a program to minimize a function of several parameters subjected to constraints and bounds. Just in case you want to run the program, the function is given by:

def Fnret(mins):
    Bj, Lj, a, b = mins.reshape(4,N)
    S1 = 0; S2 = 0
    Binf = np.zeros(N); Linf = np.zeros(N);   
    for i in range(N):
        sbi=(Bi/2); sli=(Li/2)
        for j in range(i+1):
            sbi -= Bj[j]
            sli -= Lj[j]
        Binf[i]=sbi
        Linf[i]=sli
    
    for i in range(N):
        S1 += (C*(1-np.sin(a[i]))+T*np.sin(a[i])) * ((2*Bj[i]*Binf[i]+Bj[i]**2)/(np.tan(b[i])*np.cos(a[i]))) +\
            (C*(1-np.sin(b[i]))+T*np.sin(b[i])) * ((2*Bj[i]*Linf[i]+Lj[i]*Bj[i])/(np.sin(b[i])))
    S2 += (gamma*Bj[0]/(6*np.tan(b[0])))*((Bi/2)*(Li/2) + 4*(Binf[0]+Bj[0])*(Linf[0]+Lj[0]) + Binf[0]*Linf[0])
    j=1
    while j<(N):
        S2 += (gamma*Bj[j]/(6*np.tan(b[j])))*(Binf[j-1]*Linf[j-1] + 4*(Binf[j]+Bj[j])*(Linf[j]+Lj[j]) + Binf[j]*Linf[j])
        j += 1
    F = 2*(S1+S2)
    return F

where Bj,Lj,a, and b are the minimization results given by N-sized arrays with N being an input of the program, I double-checked the function and it is working correctly. My constraints are given by:

def Brhs(mins):  # Constraint over Bj
    return np.sum(mins.reshape(4,N)[0]) - Bi

def Lrhs(mins): # Constraint over Lj
    return np.sum(mins.reshape(4,N)[1]) - Li

cons = [{'type': 'eq', 'fun': lambda Bj: 1.0*Brhs(Bj)},
         {'type': 'eq', 'fun': lambda Lj: 1.0*Lrhs(Lj)}]

In such a way that the sum of all components of Bj must be equal to Bi (same thing with Lj). The bounds of the variables are given by:

bounds = [(0,None)]*2*N + [(0,np.pi/2)]*2*N

For the problem to be reproducible, it's important to use the following inputs:

# Inputs:
gamma = 17.     
C = 85.        
T = C         
Li = 1.        
Bi = 0.5           
N = 3           

For the minimization, I'm using the cyipopt library (that is just similar to the scipy.optimize). Then, the minimization is given by:

from cyipopt import minimize_ipopt
x0 = np.ones(4*N) # Initial guess 
res = minimize_ipopt(Fnret, x0, constraints=cons, bounds=bounds)

The problem is that the result is not obeying the conditions I imposed on the constraints (i.e. the sum of Bj or Lj values is different than Bi or Li on the results). But, for instance, if I only write one of the two constraints (over Lj or Bj) it works fine for that variable. Maybe I'm missing something when using 2 constraints and I can't find the error, it seems that it doesn't work with both constraints together. Any help would be truly appreciated. Thank you in advance!

P.S.: In addition, I would like that the function result F to be positive as well. How can I impose this condition? Thanks!

Answer 1

So, based on @joni suggestions, I could find a stationary point respecting the constraints by adopting the trust-constr method of scipy.optimize.minimize library. My code is running as follows:

import numpy as np
from scipy.optimize import minimize

# Inputs:
gamma = 17      
C = 85.        
T = C         
Li = 2.         
Bi = 1.           
N = 3          # for instance

# Constraints:
def Brhs(mins):  
    return np.sum(mins.reshape(4,N)[0]) - Bi/2
def Lrhs(mins):
    return np.sum(mins.reshape(4,N)[1]) - Li/2

# Function to minimize:
def Fnret(mins):
    Bj, Lj, a, b = mins.reshape(4,N)
    S1 = 0; S2 = 0
    Binf = np.zeros(N); Linf = np.zeros(N);   
    for i in range(N):
        sbi=(Bi/2); sli=(Li/2)
        for j in range(i+1):
            sbi -= Bj[j]
            sli -= Lj[j]
        Binf[i]=sbi
        Linf[i]=sli
    
    for i in range(N):
        S1 += (C*(1-np.sin(a[i]))+T*np.sin(a[i])) * ((2*Bj[i]*Binf[i]+Bj[i]**2)/(np.tan(b[i])*np.cos(a[i]))) +\
            (C*(1-np.sin(b[i]))+T*np.sin(b[i])) * ((2*Bj[i]*Linf[i]+Lj[i]*Bj[i])/(np.sin(b[i])))
    S2 += (gamma*Bj[0]/(6*np.tan(b[0])))*((Bi/2)*(Li/2) + 4*(Binf[0]+Bj[0])*(Linf[0]+Lj[0]) + Binf[0]*Linf[0])
    j=1
    while j<(N):
        S2 += (gamma*Bj[j]/(6*np.tan(b[j])))*(Binf[j-1]*Linf[j-1] + 4*(Binf[j]+Bj[j])*(Linf[j]+Lj[j]) + Binf[j]*Linf[j])
        j += 1
    F = 2*(S1+S2)
    return F

eps = 1e-3
bounds = [(0,None)]*2*N + [(0+eps,np.pi/2-eps)]*2*N       # Bounds

cons = ({'type': 'ineq', 'fun': Fnret},
        {'type': 'eq', 'fun':  Lrhs},
        {'type': 'eq', 'fun':  Brhs})

x0 = np.ones(4*N) # Initial guess
res = minimize(Fnret, x0, method='trust-constr', bounds = bounds, constraints=cons, tol=1e-6)

F = res.fun
Bj = (res.x).reshape(4,N)[0]
Lj = (res.x).reshape(4,N)[1]
ai = (res.x).reshape(4,N)[2]
bi = (res.x).reshape(4,N)[3]

Which is essentially the same just changing the minimization technique. From np.sum(Bj) and np.sum(Lj) is easy to see that the results are in agreement with the constraints imposed, which were not working previously.