Gekko in pandas workflow

Question

I have an optimisation problem I'd like to solve with Gekko. Currently I'm using genetic algos.

I have a function to optimise that starts from an initial set and after a bunch of transformation returns a real number. This function is not linear, does not have analytical form, no idea about first derivatives and contains real and integer parameters. I'm trying to use APOPT solver of Gekko (with Gekko.options.SOLVER=1). This is my first time with Gekko.

My problem is the following: one of integer parameter is the size of a moving windows used in pandas.

When I create a Gekko variable, it encodes the variable in a string for internal computation and when this variable in passed to pandas the program fails because of this.

I tried to pass to pandas size.value.value to cast it as an int as a GKVariable variable, but it doesn't work.

Any suggestion for this problem?

Thanks in advance.

Here a minimal example:

from gekko import GEKKO
import numpy as np
import pandas as pd

np.random.seed(seed=1)

def fake_function(windows: int):
    df = pd.DataFrame({'B': np.random.randint(10, size=(20,))})
    return  df.rolling(windows).sum().values[0]

m = GEKKO() # Initialize gekko
m.options.SOLVER=1  # APOPT is an MINLP solver
gekko_param = m.Var(lb=0, ub=10, integer=True)
m.Obj(fake_function(gekko_param))
m.solve(disp=True) # Solve
print(f'Objective: {m.options.objfcnval}')

Answer 1

The problem you appear to be having is with the windowing and constructing an Objective function that "sums up some arbitrarily complicated, non-linear stuff related to a data frame and variables under optimization over several windows." I hope I got that right. :)

So, the causal factor here is twofold.

1. You are trying to use a gekko variable as a basis for iteration. It is an unknown value that is determined by the solver after your model is built, not during construction.

2. You are attempting to pass in an iterative or callback approach to defining your objective function. This must be computed in advance and be an expression in terms of numbers and model variables.

The thing that you need to do is pre-compute an expression of the objective value in terms of whatever calculations you need to do and the model variables. (See my example below, derived from yours). This will likely mean pre-computing all of the windows and generating whatever terms are needed, including model variables (as I show) and turning the result of that over to your solver, in the form of an expression.

Remember, the solver does not have access to your functions. It needs an expression in terms of finite numbers and model variables to chew on. Give it what it needs... ;)

Caveat: I'm not a regular GEKKO user, and I'm not sure what the warning is below, but this appears to be "breathing."

Code:

from gekko import GEKKO
import numpy as np
import pandas as pd

df = pd.DataFrame({'B': [2, -5, 2, -4, 8, 11]})

m = GEKKO() # Initialize gekko
m.options.SOLVER=1  # APOPT is an MINLP solver
my_scalar = m.Var(lb=-5, ub=5, integer=True)

def fake_function(window_size:int, scalar):
    return df['B'][:window_size].sum() * scalar

# build the OBJ function incrementally from "windows"
expr = 0
for i in range(1, len(df)):
    t = fake_function(i, my_scalar)
    print(f'adding: {t}')
    expr += t
print(f'The full OBJ: {expr}')

m.Obj(expr)
m.solve(disp=True) # Solve
print(f'Objective: {m.options.objfcnval}')
print(f'The value of the scalar is: {my_scalar.value}')

Yields:

adding: ((2)*(int_v1))
adding: ((-3)*(int_v1))
adding: ((-1)*(int_v1))
adding: ((-5)*(int_v1))
adding: ((3)*(int_v1))
The full OBJ: (((((0+((2)*(int_v1)))+((-3)*(int_v1)))+((-1)*(int_v1)))+((-5)*(int_v1)))+((3)*(int_v1)))
apm 73.222.107.33_gk_model0 <br><pre> ----------------------------------------------------------------
 APMonitor, Version 1.0.1
 APMonitor Optimization Suite
 ----------------------------------------------------------------
 
 
 Warning: there is insufficient data in CSV file 73.222.107.33_gk_model0.csv
 
 --------- APM Model Size ------------
 Each time step contains
   Objects      :            0
   Constants    :            0
   Variables    :            1
   Intermediates:            0
   Connections  :            0
   Equations    :            1
   Residuals    :            1
 
 Number of state variables:              1
 Number of total equations: -            0
 Number of slack variables: -            0
 ---------------------------------------
 Degrees of freedom       :              1
 
 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:      0.00 NLPi:    3 Dpth:    0 Lvs:    0 Obj: -2.00E+01 Gap:  0.00E+00
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   1.430000000982545E-002 sec
 Objective      :   -20.0000000000000     
 Successful solution
 ---------------------------------------------------
 
Objective: -20.0
The value of the scalar is: [5.0]

Answer 2

If the function derivatives aren't continuous or the function can't be evaluated at non-integer values then an approximation can be generated to provide them. The cspline (1D) or bspline (2D) functions are available to create a continuous function for optimization with gradient-based solvers. Here is the optimal with a slightly modified objective to maximize versus minimize based on the window size.

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

np.random.seed(seed=1)
df = pd.DataFrame({'B': np.random.randint(10, size=(20,))})

def fake_function(windows: int):
    return df.rolling(windows).sum().dropna().sum().values[0]

# generate data from function
wnd = np.linspace(1,15,15).astype(int) 
obj = []
for wi in wnd:
    print(fake_function(wi))
    obj.append(fake_function(wi))

m = GEKKO() # Initialize gekko
w = m.Var(lb=1, ub=15, integer=True)
ff = m.Var()
m.cspline(w,ff,wnd,obj,True) # create a continuous function
# maximize fake function (sum) based on window
m.Maximize(ff)
m.options.SOLVER=1  # APOPT is an MINLP solver
m.solve(disp=True) # Solve
print(f'Objective: {m.options.objfcnval}')

# plot results
plt.figure(figsize=(6,3))
plt.plot(wnd,obj,'k:',label='Fake Function')
plt.plot(w.value[0],ff.value[0],'ro',label='Optimal Pt')
plt.xlabel('Window'); plt.ylabel('Objective'); plt.legend()
plt.tight_layout(); plt.savefig('objective.png',dpi=300)