I am trying to add a lookup tables in my GEKKO code as a function of an m.Intermediate variable, but I can't get it to work. When I run the below code, I get the error:
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
I tried setting T as an m.Var and as an m.Param (not preferrable) as well, but it did not work. The solution does work when I use:
A240[i][j] = m.Intermediate(p1(T[i][j]))
#instead of:
A240[i][j] = m.Intermediate(A240lookup[T[i][j]])
This is not the entire code because I figured the unneccesary stuff would get in the way of an answer. I did notice that APMonitor has a lookup object (http://apmonitor.com/wiki/index.php/Main/Objects), but I did not know if that would work or even how to implement in my GEKKO model.
Does anyone have an idea to get this problem to work? Below is my code:
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)
HERE IS THE FULL CODE
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
A cubic spline could work for your application. This creates a continuous function from your data from
T_dataandA240_data. You may consider bounds onTif it becomes a variable to avoid extrapolation error.Solver summary