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

user_9172718 0
Viewed 168

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()
gekko
1 Answers

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
 ---------------------------------------------------
user_23669410
Related