Python and GEKKO optimisation of a blackbox function

Question

I am using GEKKO for fitting purposes trying to optimise functions which are explicitly defined - so I have a fully functional form and can create equation objects for optimisation purposes.

But now I have a different problem. I can't create equations because of the complexed functional dependence. But I have a python function that calculates the output using some inputs - optimisation parameters and some other that can be interpreted as fixed or known.

---

The key moments: I have the experimental data and a complexed model that is described in f1(set_of_parameters) - python function. f1 - is nonlinear, nonconvex and it can't be expressed as one simple equation - it has a lot of conditional parameters and a lot of branches the calls of other python functions inside, etc. So actually f1 can't be converted to a gekko model equation.

And I need to find such parameters - set_of_optimal_parameters, which will lead to the minimum of a distance so that f1(set_of_optimal_parameters) will be as close as possible to the experimental data I have, so I will find a set_of_optimal_parameters.

For each parameter of a set, I have initial values and boundaries and even some constraints. So I need to do something like this:

---

m = GEKKO()

#parameters I need to find
param_1 = m.FV(value = val_1, lb = lb_1, ub=rb_1)
param_2 = m.FV(value = val_2, lb = lb_2, ub=rb_2)
...
param_n = m.FV(value = val_n, lb = lb_n, ub=rb_n) #constructing the input for the function f1()
params_dataframe = ....()# some function that collects all the parameters and arranges all of them to a proper form to an input of f1()

#exp data description
x = m.Param(value = xData)
z = m.Param(value = yData)

y = m.Var()

#model description - is it possible to use other function inside equation? because f1 is very complexed with a lot of branches and options.. I don't really want to translate it in equation form..

m.Equation(
y==f1(params_dataframe)
)

#add some constraints
min = m.Param(value=some_val_min)
m.Equation(min <= (param_1+param_2) / (param_1+param_2)**2))

# trying to solve and minimize the sum of squares
m.Minimize(((y-z))**2)

# Options for solver
param_1.STATUS = 1
param_2.STATUS = 1
...
param_n.STATUS = 1

m.options.IMODE = 2
m.options.SOLVER = 1
m.options.MAX_ITER = 1000

m.solve(disp=1)

Is it possible to use GEKKO this way or it's not allowed? and why?

Answer 1

Gekko compiles equations into byte-code and requires all equations in Gekko format so that it can overload equation operators to provide exact first and second derivatives in sparse form. Black-box functions do not provide the necessary first and second derivatives, but they can provide function evaluations for finite differences (derivative approximations) or for surrogate functions.

To answer your question directly, you can't use f1(params) in a Gekko problem. If you need an optimizer to evaluate arbitrary black box functions, an optimizer such as scipy.optimize.minimize() is a good choice.

If you would still like to use Gekko, there are several options to built a surrogate model for f1 that has continuous first and second derivatives. The surrogate model depends on the number of params:

• 1D: use cspline()
• 2D: use bspline()
• 3D+: use Machine learning such as Gaussian Processes, Neural Network, Linear Regression, etc.

Here is an example that create a surrogate model for y=f(x) where f(x)=3*np.sin(x) - (x-3). This equation could be modeled directly in Gekko, but it serves as an example of creating a cspline() object that approximates the function and finds the minimum.

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt

"""
minimize y
s.t.     y = f(x)
using cubic spline with random sampling of data
"""

# function to generate data for cspline
def f(x):
    return 3*np.sin(x) - (x-3)

x_data = np.random.rand(50)*10+10
y_data = f(x_data)

c = GEKKO()
x = c.Var(value=np.random.rand(1)*10+10)
y = c.Var()
c.cspline(x,y,x_data,y_data,True)
c.Obj(y)
c.options.IMODE = 3
c.options.CSV_READ = 0
c.options.SOLVER = 3
c.solve(disp=True)

if c.options.SOLVESTATUS == 1:
    plt.figure()
    plt.scatter(x_data,y_data,5,'b')
    plt.scatter(x.value,y.value,200,'r','x')
else:
    print ('Failed!')
    print(x_data,y_data)
    plt.figure()
    plt.scatter(x_data,y_data,5,'b')
plt.show()