Использование таблиц поиска в Gekko Python Intermediate Equation

Я пытаюсь добавить таблицы поиска в свой код GEKKO как функцию переменной m.Intermediate, но не могу заставить ее работать. Когда я запускаю приведенный ниже код, я получаю сообщение об ошибке:

A240[i][j] = m.Intermediate(A240lookup[T[i][j]])

#IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices 

Я попытался установить T как m.Var и как m.Param (нежелательно), но это не сработало. Решение действительно работает, когда я использую:

A240[i][j] = m.Intermediate(p1(T[i][j]))

#instead of:

A240[i][j] = m.Intermediate(A240lookup[T[i][j]])

Это не весь код, потому что я полагал, что ненужные вещи будут мешать ответу. Я заметил, что APMonitor имеет объект поиска (http://apmonitor.com/wiki/index.php/Main/Objects), но я не знал, сработает ли это и даже как реализовать в моей модели GEKKO.

Есть у кого-нибудь идея заставить эту проблему работать? Ниже мой код:

import numpy as np
from gekko import GEKKO
m = GEKKO()
m.options.SOLVER= 1

n =      4             
p =      2
Tr =     300
temp =11
t = np.linspace(300,1073,15)
Temp = np.array([1023, 973,923,873,823,723, 623, 523, 423, 323, 273])

ATemp240 = np.array([870, 918,956,980,1002,1015,1019,1021,1021.7, 1022.1,1022.4])
p1 = np.poly1d(np.polyfit(Temp,ATemp240,7))

A240lookup = m.Array(m.Param,(1000))

for i in range(1000):
    A240lookup[i].value = p1(i)

F = m.Array(m.Var,(n,p), integer = True) #State in or out
for i in range (n):
    for j in range(p):
        F[i][j].upper = 1
        F[i][j].lower = 0
        
T     = [[[]for j in range(p)]for i in range(n)]
A240  = [[[]for j in range(p)]for i in range(n)]
    
for j in range(p):
      for i in range(1,n):
          T [i][j]   = m.Intermediate(Tr+30*i)
          A240[i][j] = m.Intermediate(A240lookup[T[i][j]])    #I want this
          # A240[i][j] = m.Intermediate(p1(T[i][j]))          #But this is all that works

m.solve(debug=False)

ЗДЕСЬ ПОЛНЫЙ КОД

from mpl_toolkits.axes_grid1 import make_axes_locatable
import numpy as np
import math
import matplotlib.pyplot as plt
from gekko import GEKKO

m = GEKKO()
m.options.SOLVER= 1
m.options.max_iter=3000
# m.options.IMODE=2
# m.options.REDUCE = 3
# m.options.DIAGLEVEL = 1
m.solver_options = ['minlp_print_level 10',\
                    'minlp_maximum_iterations 3000',\
                    'minlp_as_nlp 0',\
                    'minlp_branch_method 2']

# m.options.csv_write=2
n =      120  #time steps (rows)
p =      2    #parts (columns)
ts =     120  #seconds per time step
B =      5    #Block (number of time steps per block)
Tf =     810  #furnace temperature
Tr =     300  #room temperature
H =      .005 #part height (m)
Hsq =    H*H
H2 =     2*H
z =      0
n1 =     -0.6
Pi =     3.141592653539
pi2sqr = (Pi/2)**2
Piz =    Pi*z
aC = 0.000000005
aH = 0.00000005

t = np.linspace(300,1073,15)

T_data = [1023, 973,923,873,823,723, 623, 523, 423, 323, 273]
A240_data = [870, 918,956,980,1002,1017,1019,1021,1021.7, 1022.1,1022.4]
s240_data = [10.4, 11.7, 17, 32, 50, 93, 138, 175, 205, 221, 222]
B240_data = [8.42705, 9.26699, 10.26703, 11.94154, 14.30153, 22.08034, 39.74829, 81.99357, 221.00829, 650, 950]

tc_data = 120*np.arange(0,n)
def qC(t):
    return(sum((4*((-1)**x))/((2*x+1)*Pi)*math.exp((-1*((2*x+1)*Pi/2)**2)*t*aC/(H**2))*math.cos((2*x+1)/2*Pi*z/H)for x in range(n)))
# QC_data =qH(tc_data)
QC_data = np.zeros(n)
for i in range(n):
    QC_data[i] = qC(tc_data[i])
    
def qH(t):
    return(sum((4*((-1)**x))/((2*x+1)*Pi)*math.exp((-1*((2*x+1)*Pi/2)**2)*t*aH/(H**2))*math.cos((2*x+1)/2*Pi*z/H)for x in range(n)))
# QC_data =qH(tc_data)
QH_data = np.zeros(n)
for i in range(n):
    QH_data[i] = qH(tc_data[i])
 
F = m.Array(m.Param,(n,p)) #State in or out
for i in range(n):
    for j in range(p):
        F[i][j].value = 1-((i//B+j)%2)
 
tc   = m.Array(m.Var,(p,n)) # time in current state
QC   = m.Array(m.Var,(p,n)) # Quenching Cool
QH   = m.Array(m.Var,(p,n)) # Quenching Heat
T    = m.Array(m.Var,(p,n)) # Temperature
A240 = m.Array(m.Var,(p,n)) # prediction results
s240 = m.Array(m.Var,(p,n)) # prediction results
B240 = m.Array(m.Var,(p,n)) # prediction results
 
Q     = [[[]for i in range(n)]for j in range(p)] # Quenching
S     = [[[]for i in range(n)]for j in range(p)] # Residual Stress
tr    = [[[]for i in range(n)]for j in range(p)] # Relative Time
Tc    = [[[]for i in range(n)]for j in range(p)] 
var   = [[[]for i in range(n)]for j in range(p)]

for j in range(p):
    m.Equation(tc [j][0] == 0)
    m.Equation(T[j][0] == Tr)
    m.cspline(T[j][0],A240[j][0],T_data,A240_data)
    m.cspline(T[j][0],B240[j][0],T_data,B240_data)
    m.cspline(T[j][0],s240[j][0],T_data,s240_data)
    m.cspline(tc[j][0],QC[j][0],tc_data,QC_data)
    m.cspline(tc[j][0],QH[j][0],tc_data,QH_data) 
    tr   [j][0] = m.Intermediate(0)
    S    [j][0] = m.Intermediate(240)
    Tc   [j][0] = (Tr)

    
    for i in range(1,n):
        m.Equation(tc [j][i] == (tc[j][i-1])*(1-m.abs(F[i][j]-F[i-1][j]))+ts)
        m.cspline(tc[j][i],QC[j][i],tc_data,QC_data)
        m.cspline(tc[j][i],QH[j][i],tc_data,QH_data)
        Q[j][i] = m.Intermediate(QC[j][i]*(1-F[i][j])+(QH[j][i]*F[i][j]))
        Tc    [j][i] = m.Intermediate(Tc[j][i-1]-(m.abs(F[i][j]-F[i-1][j])*(Tc[j][i-1]-T[j][i-1])))
        m.Equation(T [j][i] == Q[j][i]*(Tc[j][i]-(Tr*(1-F[i][j])+(Tf*F[i][j])))+(Tr*(1-F[i][j])+(Tf*F[i][j])))        
        S     [j][i] = m.Intermediate(m.min2(S[j][i-1],(A240[j][i-1]*((tr[j][i-1]+ts+B240[j][i-1])**(n1))+s240[j][i-1])))
        m.cspline(T[j][i],A240[j][i],T_data,A240_data)
        m.cspline(T[j][i],B240[j][i],T_data,B240_data)
        m.cspline(T[j][i],s240[j][i],T_data,s240_data)
        var[j][i] = m.Intermediate(m.min2(S[j][i]-.1,s240[j][i]))
        tr   [j][i] = m.Intermediate(((S[j][i]-var[j][i])/A240[j][i])**(1/n1)-B240[j][i])

m.Minimize(sum(S[j][n-1] for j in range(p)))

m.solve(debug=False)    # solve