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

Not a complete answer, just some hints in arbitrary order:

• Your initial point x0 is not feasible because it contradicts both of your constraints. This can easily be observed by evaluating your constraints at x0. Under the hood, Ipopt typically detects this and tries to find a feasible initial point. However, it's highly recommended to provide a feasible initial point whenever possible.
• Your variable bounds are not well-defined. Evaluating your objective at your bounds yields multiple divisions by zero. For example, the denominator np.tan(b[i]) is zero if and only if b[i] = 0, so 0 isn't a valid value for all of your b[i]s. Proceeding similarly for the other terms, you should obtain 0 < b[i] < pi/2 and 0 ≤ a[i] < pi/2. Here, you can model the strict inequalities by 0 + eps ≤ b[i] ≤ pi/2 - eps and 0 ≤ a[i] ≤ pi/2 - eps, where eps is a sufficiently small positive number.
• If you really want to impose that the objective function is always positive, you can simply add the inequality constraint Fnret(x) >= 0, i.e. {'type': 'ineq', 'fun': Fnret}.

In Code:

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

# (feasible) initial guess
x0 = eps*np.ones(4*N)
x0[[0, N]] = [Bi-(N-1)*eps, Li-(N-1)*eps]

# constraints
cons = [{'type': 'eq', 'fun': Brhs},
        {'type': 'eq', 'fun': Lrhs},
        {'type': 'ineq', 'fun': Fnret}]

res = minimize_ipopt(Fnret, x0, constraints=cons, bounds=bounds, options={'disp': 5})

Last but not least, this still doesn't converge to a stationary point, so chances are that there's indeed no local minimum. From here, you can try experimenting with other (feasible!) initial points and double-check the math of your problem. It's also worth providing the exact gradient and constraint Jacobians.

Answer 2

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.