Why does the code terminate with a "Solution Not Found" error and "EXIT: Converged to a point of local infeasibility. Problem may be infeasible"?

Question

I cannot seem to figure out why IPOPT cannot find a solution to this. Initially, I thought the problem was totally infeasible but when I reduce the value of col_total to any number below 161000 or comment out the last constraint equation that contains col_total, it solves and EXITs with an Optimal Solution Found and a final objective value function of -161775.256826753. I have solved the same Maximization problem using Artificial Bee Colony and Particle Swamp Optimization techniques, and they solve and return optimal objective value function at least 225000 and 226000 respectively. Could it be that another solver is required? I have also tried APOPT, BPOPT, and IPOPT and have tinkered around with the tolerance values, but no combination none seems to work just yet. The code is posted below. Any guidance will be hugely appreciated.

from gekko import GEKKO  
import numpy as np  

distances = np.array([[[0, 0],[0,0],[0,0],[0,0]],\
                      [[155,0],[0,0],[0,0],[0,0]],\
                      [[310,0],[155,0],[0,0],[0,0]],\
                      [[465,0],[310,0],[155,0],[0,0]],\
                      [[620,0],[465,0],[310,0],[155,0]]])

        
alpha = 0.5 / np.log(30/0.075)  
diam = 31  
free = 7  
rho = 1.2253  
area = np.pi * (diam / 2)**2  
min_v = 5.5  
axi_max = 0.32485226746  
col_total = 176542.96546512868    
rat = 14  
nn = 5
u_hub_lowerbound = 5.777777777777778
c_pow = 0.59230249
p_max = 0.5 * rho * area * c_pow * free**3 

# Initialize Model  
m = GEKKO(remote=True)  
#initialize variables, Set lower and upper bounds  
x = [m.Var(value = 0.03902278, lb = 0, ub = axi_max) \
     for i in range(nn)]  

# i = 0  
b = 1  
c = 0  
v_s = list()  
for i in range(nn-1):       # Loop runs for nn-1 times  
        # print(i)
        # print(i,b,c)
        squared_defs = list()  
        while i < b:  

                d = distances[b][c][0]  
                r = distances[b][c][1]  
                ss = (2 * (alpha * d) / diam)  
                tt = r / ((diam/2) + (alpha * d))
                squared_defs.append((2 * x[i] / (1 + ss**2)) * np.exp(-(tt**2)) ** 2)  

                i+=1  
                c+=1  

        #Equations  
        m.Equation((free * (1 - (sum(squared_defs))**0.5)) - rat <= 0)
        m.Equation((free * (1 - (sum(squared_defs))**0.5)) - u_hub_lowerbound >= 0)  
        v_s.append(free * (1 - (sum(squared_defs))**0.5))  
        squared_defs.clear()  
        b+=1   
        c=0  


# Inserts free as the first item on the v_s list to 
# increase len(v_s) to nn, so that 'v_s' and 'x' 
# are of same length  
v_s.insert(0, free)                 
 
gamma = list()  
for i in range(len(x)):  

        bet = (4*x[i]*((1-x[i])**2) * rho * area) / 2 
        gam = bet * v_s[i]**3
        gamma.append(gam)
        #Equations
        m.Equation(x[i] - axi_max <= 0)  
        m.Equation((((4*x[i]*((1-x[i])**2) * rho * area) / 2) \
                    * v_s[i]**3) - p_max <= 0)  
        m.Equation((((4*x[i]*((1-x[i])**2) * rho * area) / 2) * \
                      v_s[i]**3) > 0)
          

#Equation  
m.Equation(col_total - sum(gamma) <= 0)  

#Objective  
y = sum(gamma)  
m.Maximize(y)  # Maximize  

#Set global options  
m.options.IMODE = 3 #steady state optimization  
#Solve simulation  
m.options.SOLVER = 3
m.solver_options = ['linear_solver ma27','mu_strategy adaptive','max_iter 2500', 'tol 1.0e-5' ]  
m.solve()

@John Hedengren, the solution you provided to the initial question worked well in the bigger model but once I changed the values of axi_max, introduced constants: c_pow and u_hub_lowerbound, and changed the second constraint equation, it started to misbehave — Fonnie, May 06 '21 at 09:54
How does it misbehave? Does it reach maximum iterations or does it report an infeasible solution? Try another solver such as `m.options.SOLVER=1` (APOPT) to verify an infeasibility. — John Hedengren, May 06 '21 at 15:40
Yea, with `APOPT`, this standalone script runs till `maximum iterations of 250` and outputs an objective function of approximately **176523** which is not far off from `col_total` and still returns " @error: Solution Not Found". However, the script connected to my system model (which contains the same thing on the standalone script, except that some parameters are inherited from other scripts) reaches maximum iterations and returns an objective of approximately 224230 which is not far off from the values produced by `PSO` and `ABC` heuristic algorithms, but returns @error: Solution Not Found. — Fonnie, May 06 '21 at 16:31
Iter Objective Convergence 250 -2.29336E+05 5.40672E-11 Maximum iterations --------------------------------------------------- Solver : APOPT (v1.0) Solution time : 0.493200000011711 sec Objective : -229336.244661586 Unsuccessful with error code 0 --------------------------------------------------- Creating file: infeasibilities.txt Use command apm_get(server,app,'infeasibilities.txt') to retrieve file @error: Solution Not Found — Fonnie, May 06 '21 at 22:43
`APOPT` solves and produces a larger `Maximized Objective` as shown above, but still returns a `Solution Not Found` error. Could adjusting any of the tolerance values avert this error? — Fonnie, May 06 '21 at 22:46
Yes, the `minlp_gap_tol` can help. Adding integer constraints always makes the objective function the same or worse. https://apmonitor.com/wiki/index.php/Main/IntegerProgramming — John Hedengren, May 07 '21 at 17:22
I have found a solution, finally. I read your response here-> https://stackoverflow.com/questions/67146271/gekko-parameter-estimation-with-custom-objective-function-error-code-13/67168507#67168507, and increased my lower bound (`lb`) from `0 to 0.01`, to ensure that no variable picks the value 0 after initialization. With that adjustment, `APOPT, BPOPT, and IPOPT` now solve with a Maximized Objective of approximately **238767** which is significantly larger than the results produced by `PSO and ABC` heuristic algorithms, as expected. — Fonnie, May 07 '21 at 18:53
I stumbled on that response after increasing `max_iter to 10000` and using the `linear solver mumps` with `IPOPT`. It exited with the error `EXIT: Invalid number in NLP function or derivative detected` and an error code of `-13`. I searched that error and it took me to your response on the link posted above and another on https://stackoverflow.com/questions/64491175/python-optimization-using-gekko/64520777#64520777. — Fonnie, May 07 '21 at 19:03

score 0 · Accepted Answer · answered Apr 21 '21 at 17:01

Built the equations without .value in the expressions. The x[i].value is only needed at the end to view the solution after the solution is complete or to initialize the value of x[i]. The expression m.Maximize(y) is more readable than m.Obj(-y) although they are equivalent.

from gekko import GEKKO  
import numpy as np  

distances = np.array([[[0, 0],[0,0],[0,0],[0,0]],\
                      [[155,0],[0,0],[0,0],[0,0]],\
                      [[310,0],[155,0],[0,0],[0,0]],\
                      [[465,0],[310,0],[155,0],[0,0]],\
                      [[620,0],[465,0],[310,0],[155,0]]])  
        
alpha = 0.5 / np.log(30/0.075)  
diam = 31  
free = 7  
rho = 1.2253  
area = np.pi * (diam / 2)**2  
min_v = 5.5  
axi_max = 0.069262150781  
col_total = 20000  
p_max = 4000  
rat = 14  
nn = 5  

# Initialize Model  
m = GEKKO(remote=True)  
#initialize variables, Set lower and upper bounds  
x = [m.Var(value = 0.03902278, lb = 0, ub = axi_max) \
     for i in range(nn)]  

i = 0  
b = 1  
c = 0  
v_s = list()  
for turbs in range(nn-1):       # Loop runs for nn-1 times  
    squared_defs = list()  
    while i < b:  
        d = distances[b][c][0]  
        r = distances[b][c][1]  
        ss = (2 * (alpha * d) / diam)  
        tt = r / ((diam/2) + (alpha * d))  
        squared_defs.append((2 * x[i] / (1 + ss**2)) \
                            * m.exp(-(tt**2)) ** 2)  

        i+=1  
        c+=1  

    #Equations  
    m.Equation((free * (1 - (sum(squared_defs))**0.5)) - rat <= 0)
    m.Equation(min_v - (free * (1 - (sum(squared_defs))**0.5)) <= 0 )  
    v_s.append(free * (1 - (sum(squared_defs))**0.5))  
    squared_defs.clear()  
    b+=1  
    a=0  
    c=0  

# Inserts free as the first item on the v_s list to
#   increase len(v_s) to nn, so that 'v_s' and 'x'
#   are of same length  
v_s.insert(0, free)                 
beta = list()  
gamma = list()  
for i in range(len(x)):  
    bet = (4*x[i]*((1-x[i])**2) * rho * area) / 2 
    gam = bet * v_s[i]**3
    #Equations  
    m.Equation((((4*x[i]*((1-x[i])**2) * rho * area) / 2) \
                * v_s[i]**3) - p_max <= 0)  
    m.Equation((((4*x[i]*((1-x[i])**2) * rho * area) / 2) \
                * v_s[i]**3) > 0)
    gamma.append(gam)  

#Equation  
m.Equation(col_total - sum(gamma) <= 0)  

#Objective  
y = sum(gamma)  
m.Maximize(y)  # Maximize  

#Set global options  
m.options.IMODE = 3 #steady state optimization  
#Solve simulation  
m.options.SOLVER = 3  
m.solve()

This gives a successful solution with maximized objective 20,000:

Number of Iterations....: 12

                                   (scaled)                 (unscaled)
Objective...............:  -4.7394814741924645e+00   -1.9999999999929641e+04
Dual infeasibility......:   4.4698510326511536e-07    1.8862194343304290e-03
Constraint violation....:   3.8275766582203308e-11    1.2941979026166479e-07
Complementarity.........:   2.1543608536533588e-09    9.0911246952931704e-06
Overall NLP error.......:   4.6245685940749926e-10    1.8862194343304290e-03


Number of objective function evaluations             = 80
Number of objective gradient evaluations             = 13
Number of equality constraint evaluations            = 80
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 13
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 12
Total CPU secs in IPOPT (w/o function evaluations)   =      0.010
Total CPU secs in NLP function evaluations           =      0.011

EXIT: Optimal Solution Found.
 
 The solution was found.
 
 The final value of the objective function is   -19999.9999999296     
 
 ---------------------------------------------------
 Solver         :  IPOPT (v3.12)
 Solution time  :   3.210000000399305E-002 sec
 Objective      :   -19999.9999999296     
 Successful solution
 ---------------------------------------------------

Awesome!!! It works perfectly at my end too. I'll slot that into the bigger model and see as it goes. Thanks for the timely response. Well appreciated. — Fonnie, Apr 21 '21 at 18:00

Why does the code terminate with a "Solution Not Found" error and "EXIT: Converged to a point of local infeasibility. Problem may be infeasible"?

1 Answers1

Linked