Using spline or piece wise function as input (intermediate) in Gekko

Simple paraphrase of my question is what is the best way to use a table or curve as an input for an optimization problem in Gekko (and/or AP Monitor)?

For example, I have simplified atmospheric density data by altitude. I can fit the data to a a cubic spline in hermite form (slopes at control points are given explicitly, ref wikipedia) with about 5 control points and I have a python function to interpolate, but as far as I know I can’t use a python function in a Gekko model.

My understanding and assumptions on possibilities:

Typically the first choice would be to fit the data to an equation, however in this case it doesn’t fit. I could fit it to a 2 or 3 segment piece-wise function with sufficient accuracy, but even that isn’t always feasible. That would also involve inequality conditions and possibly discontinuities in the derivatives, which are not preferable in an optimization problem to say the least.

I can technically re-write the interpolation routine I have as model equations but it would get really messy really quickly.

It seems like the built in cspline object might work, however the examples only shows fitting the spline to the data. It is my understanding that you can have a model within a model so it seems possible to have an inner nested model that would solve for interpolated points. That seems like overkill and seems likely to have performance issues, especially if a sequential solver is required.

Similarly, AP Monitor has the pwl (piece wise liner) object and a look up object. The documentation shows the implementation of the pwl object in Gekko (basically does linear interpolation with slack variables), but again this seems like it would get messy when more than one of this type of input would be required.

So to restate my question(s) more specifically: Is there a clear best practice on how to do this? If not, are there significant deciding factors (i.e. pros/cons) of each that I’ve missed?

Bonus question – same thing but with 2D data like coefficient of drag dependent on mach and angle of attack? (using bspline object instead of cspline)

My initial assumption was to have pressure, temperature, and density data as csplines and drag and lift coefficients as bsplines. I’m still learning Gekko/AP Monitor and python so I figured it was better to ask now before I spend days trying to troubleshoot something that’s not practical.

from gekko import GEKKO
import numpy as np
#knots and coeffs
m = GEKKO(remote=False)
tx = [ -1, -1, -1, -1,  1,  1,  1,  1]
ty = [ -1, -1, -1, -1,  1,  1,  1,  1]
c = [1.0, 0.33333333, -0.33333333, -1.0, 0.33333333, 
     0.11111111, -0.11111111, -0.33333333, -0.33333333, 
    -0.11111111, 0.11111111, 0.33333333, -1.0, -0.33333333,
     0.33333333, 1.0]
x = m.Var(0.5,-1,1)
y = m.Var(0.5,-1,1)
z = m.Var(2)
m.bspline(x,y,z,tx,ty,c,data=False)
m.Obj(z)
m.solve()

from gekko import GEKKO
import numpy as np
#raw data
m = GEKKO(remote=False)
xgrid = np.linspace(-1, 1, 20)
ygrid = xgrid
z_data = x*y
x = m.Var(0.5,-1,1)
y = m.Var(0.5,-1,1)
z = m.Var(2)
m.bspline(x,y,z,xgrid,ygrid,z_data)
m.Obj(z)
m.solve()