7

I'm trying to optimise an exponential objective function with GEKKO, but I don't know if the selected solver is the best one for these kind of problems.

Is the selected one a valid choice??

import numpy as np

'GEKKO MODELING'
from gekko import GEKKO
m = GEKKO()
m.options.SOLVER=1  # APOPT is an MINLP solver

# Initialize variables
x = []
x1 = m.Var(value=20,lb=20, ub=6555)  #integer=True
x2 = m.Var(value=0,lb=0,ub=10000)  #integer=True
x3 = m.sos1([30, 42, 45, 55])

x = [x1, x2, x3]
# Equations
m.Equation((x1 * x2* x3) * 10 ** (-6)>=50)

def fun(x):
    return 44440 + ((np.pi * x[0] * x[1] * x[2]) * 10 ** (-4))**0.613

x = [400,300,19]

'GEKKO Optimization'
m.Obj(fun(x))

m.solve(disp=False) # Solve
print('Results')
print('x1: ' + str(x1.value))
print('x2: ' + str(x2.value))
print('x3: ' + str(x3.value))

print('Objective: ' + str(m.options.objfcnval))
srgam
  • 366
  • 1
  • 13
  • The objective can be simplified a lot. – Erwin Kalvelagen Aug 20 '19 at 14:08
  • Do you have some advice? – srgam Aug 20 '19 at 14:11
  • (1) Get rid of the constant 44440. (2) Get rid of **0.613. (3) get rid of pi (4) get rid of 10^-4. – Erwin Kalvelagen Aug 20 '19 at 14:33
  • 1
    The title (exponential function) implies that you have 0.613^x in your objective instead of x^0.613. The same responses should apply to both but many times an exponential function is more "nonlinear" meaning that special precautions must sometimes be used (constraints, scaling, etc) to avoid solver divergence and help it find a solution. – John Hedengren Aug 20 '19 at 16:06

1 Answers1

5

A problem with your script is that you are redefining the value of x = [400,300,19] before you call your objective function. The objective function should be called with your original definition x = [x1, x2, x3] so that it can optimize these variables. One other change is that the value of x3 is by default equal to zero. Setting it away from zero x3.value=1.0 allows APOPT and IPOPT solvers to converge because you previously started right on the border of an imaginary number objective with x3<0.

import numpy as np
from gekko import GEKKO
m = GEKKO()
x = []
x1 = m.Var(value=20,lb=20, ub=6555)  #integer=True
x2 = m.Var(value=1,lb=1,ub=10000)  #integer=True
x3 = m.sos1([30, 42, 45, 55])
x3.value = 1.0
x = [x1, x2, x3]
m.Equation((x1 * x2* x3) * 1e-6 >= 50)
def fun(x):
    return 44440 + ((np.pi * x[0] * x[1] * x[2]) * 1e-4)**0.613
m.Obj(fun(x))

# Change to True to initialize with IPOPT
init = False
if init:
    m.options.SOLVER=3  
    m.solve(disp=False) # Solve

m.options.SOLVER=1
m.solve(disp=True) # Solve

print('Results')
print('x1: ' + str(x1.value))
print('x2: ' + str(x2.value))
print('x3: ' + str(x3.value))
print('Objective: ' + str(m.options.objfcnval))

For the solver suggestion, here is a list of publicly available solvers in Gekko. There are additional commercially available solver options in Gekko but I'll stick with just the publicly accessible ones (APOPT, BPOPT, and IPOPT) for this response. Any nonlinear programming solver should be able to handle a nonlinear objective such as x**0.613. Your problem also includes a Special Ordered Set, Type 1 (m.sos1) and so your problem is not just a Nonlinear Programming (NLP) problem but also includes binary variables for the sos1. This means that you'll need to use a Mixed Integer Nonlinear Programming (MINLP) solver. The APOPT solver is the only MINLP solver publicly available in Gekko and it is automatically selected for you when you create an sos1 object. If you would like to try to solve the MINLP problem with an NLP solver (such as IPOPT), then you'll need to specify the solver after you create the m.sos1 object.

m.options.SOLVER = 3

This could lead to an incorrect solution because x3 can only be one of the following: 30, 42, 45, 55. IPOPT finds a minimum solution of x3==47.079550873 so in this case, it did not return an integer solution. If you want to guarantee an integer solution, you'll need to use APOPT.

 Successful solution

 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   4.279999999562278E-002 sec
 Objective      :    44813.4405591393     
 Successful solution
 ---------------------------------------------------

Results
x1: [677.59896405]
x2: [2459.665311]
x3: [30.0]
Objective: 44813.440559

If you need to change some of the tuning parameters of the MINLP APOPT solver then you can use something like the following:

m.solver_options = ['minlp_gap_tol 1.0e-2',\
                    'minlp_maximum_iterations 10000',\
                    'minlp_max_iter_with_int_sol 500']

There is additional information on the APOPT solver options.

John Hedengren
  • 12,068
  • 1
  • 21
  • 25
  • And what about include penalties in my objective function? How is the best way for it with GEKKO? [Penalty in an objective Function. GEKKO](https://scicomp.stackexchange.com/questions/33279/how-to-include-penalty-in-a-objective-function-with-python-gekko/33281#33281) – srgam Aug 20 '19 at 17:43
  • 1
    You can add many m.Obj() terms. Gekko adds them up and minimizes the sum because most solvers require just one objective function. You can add penalties with something like m.Obj(-x[2]) for maximize and m.Obj(x[2]) for minimize. If you want it to encourage a certain value then use something like m.Obj((x[2]-target)**2). – John Hedengren Aug 20 '19 at 20:51
  • But what about if I have this case? [See the new comment](https://scicomp.stackexchange.com/a/33281/32550) My penalties are function of the variables. I calculate a physical property from variable values (as an example, I recalculate temperature), and then I evaluate if my temperature is higher than my fixed maximum. @JohnHedengren – srgam Aug 21 '19 at 10:41