【发布时间】：2021-02-04 23:26:31
【问题描述】：
我试图在我的 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
【问题讨论】：
-
另一个想法：有没有办法将 GEKKO 运算符转换为整数？如果是这样，那么我可以使用该整数作为查找表列表的索引。
标签： python arrays optimization lookup-tables gekko