问题描述
我试图在我的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,15)
T_data = [1023,273]
A240_data = [870,1017,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,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,n)) # Quenching Cool
QH = m.Array(m.Var,n)) # Quenching Heat
T = m.Array(m.Var,n)) # Temperature
A240 = m.Array(m.Var,n)) # prediction results
s240 = m.Array(m.Var,n)) # prediction results
B240 = m.Array(m.Var,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],B240_data)
m.cspline(T[j][0],s240[j][0],s240_data)
m.cspline(tc[j][0],QC[j][0],tc_data,QC_data)
m.cspline(tc[j][0],QH[j][0],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],QC_data)
m.cspline(tc[j][i],QH[j][i],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],A240_data)
m.cspline(T[j][i],B240[j][i],B240_data)
m.cspline(T[j][i],s240[j][i],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
解决方法
cubic spline可能适用于您的应用程序。这会根据您来自T_data和A240_data的数据创建一个连续函数。如果T成为变量以避免外推误差，则可以考虑使用边界。
import numpy as np
from gekko import GEKKO
m = GEKKO()
m.options.SOLVER= 1
n = 4
p = 2
Tr = 300
# Cubic spline
T_data = [1023,973,923,873,823,723,623,523,423,323,273]
A240_data = [870,918,956,980,1002,1015,1019,1021,1021.7,1022.1,1022.4]
T = m.Array(m.Param,(n,p))
A240 = m.Array(m.Var,p))
for j in range(p):
for i in range(1,n):
T[i,j].value = Tr+30*i
m.cspline(T[i,j],A240[i,T_data,A240_data)
m.solve()
求解摘要
--------- APM Model Size ------------
Each time step contains
Objects : 6
Constants : 0
Variables : 16
Intermediates: 0
Connections : 12
Equations : 0
Residuals : 0
Number of state variables: 8
Number of total equations: - 6
Number of slack variables: - 0
---------------------------------------
Degrees of freedom : 2
----------------------------------------------
Steady State Optimization with APOPT Solver
----------------------------------------------
Iter Objective Convergence
0 1.25109E-09 1.02206E+03
1 2.44846E-23 5.37222E-17
2 2.44846E-23 5.37222E-17
Successful solution
---------------------------------------------------
Solver : APOPT (v1.0)
Solution time : 1.350000000093132E-002 sec
Objective : 0.000000000000000E+000
Successful solution
---------------------------------------------------