如何在 Gekko 中初始化 3D 数组变量？
How to initialize an 3D array variable in Gekko?
我正在尝试根据三维 主时间表模型 求解一个步骤，其中涉及时段 (5)、课程 (19) 和地点 (8)。
所以我在 Gekko 中用 3D 数组初始化这些变量时遇到了问题。如果没有这个初始化，算法在超过 15 分钟 运行 和 1000 次迭代后不会收敛。
当我尝试初始化时，出现此错误：
”
引发异常（响应）
异常：@error：方程定义
没有等式 (=) 或不等式 (>,<) 的等式
真的
停止...
“
我该如何解决这个问题？遵循我的代码版本：
import numpy as np
from gekko import GEKKO
# Input data
# Schedule of periods and courses
sched = np.array([ [0, 1, 0, 0, 1], [0, 0, 1, 1, 0], [0, 0, 1, 1, 0], \
[0, 0, 0, 0, 1], [1, 0, 0, 0, 1], [0, 0, 0, 1, 1], [0, 1, 1, 0, 0], \
[1, 0, 0, 1, 0], [0, 1, 0, 0, 1], [1, 1, 0, 0, 0], [0, 1, 1, 0, 0], \
[0, 1, 1, 0, 0], [1, 0, 0, 1, 0], [1, 0, 0, 1, 0], [0, 0, 1, 0, 1], \
[1, 0, 1, 0, 0], [0, 1, 0, 1, 0], [0, 0, 1, 1, 0], [0, 1, 0, 0, 1] ], dtype=np.int64)
# Initial allocation of all periods, courses and locations
alloc=np.array([0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,1,1,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,\
0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], dtype=np.int64)
# Number of students enrolled in each course
enrol = np.array([ 60, 60, 60, 40, 40, 110, 120, 50, 60, 55, 50, \
55, 40, 64, 72, 50, 50, 55, 55], dtype=np.float64)
# Capacity of each location (classroom)
capac = np.array([ 60, 60, 120, 60, 80, 60, 60, 65], dtype=np.float64)
# Total costs of using each location
costs = np.array([ 9017.12, 9017.12, 12050.24, 9017.12, 9413.68, 9017.12, \
9017.12, 9188.96 ])
# Estimated cost of each location by period and student
ecost = np.repeat(np.array([[costs[i]*pow(enrol[j]*5,-1) for j in range(19)] for i in range(8)]), 5)
# The model construction
m = GEKKO()
# Constant arrays
x = m.Array(m.Const,(19,5))
y = m.Array(m.Const,(8,19,5))
N = m.Array(m.Const,(19))
C = m.Array(m.Const,(8))
Ec = m.Array(m.Const,(8,19,5))
Ecy = m.Array(m.Const,(8,19,5))
Alt = m.Array(m.Const,(8,19,5))
for k in range(5):
for j in range(19):
N[j] = enrol[j]
x[j,k] = sched[j,k]
for i in range(8):
C[i] = capac[i]
Ec[i,j,k] = ecost[k+j*5+i*19*5]
y[i,j,k] = alloc[k+j*5+i*19*5]
Ecy[i,j,k] = Ec[i,j,k]*y[i,j,k]
if sched[j,k]==1:
Alt[i,j,np.where(sched[j,:]==1)[0][0]]=-sched[j,k]*(1-sum(sched[j,:]))
if sum(sched[j,:])==2:
Alt[i,j,np.where(sched[j,:]==1)[0][1]]=sched[j,k]*(1-sum(sched[j,:]))
else:
Alt[i,j,k]=0
# Initialize the variable z with the initial value y:
# These commented approaches produce the error.
z = m.Array(m.Var,(8,19,5),lb=0,ub=1,integer=True)
#for i in range(8):
# for j in range(19):
# for k in range(5):
# z[i,j,k] = y[i,j,k]
# nor
#z = m.Array(m.Var,(8,19,5),value=y,lb=0,ub=1,integer=True)
# Intermediate equations
Ecz = m.Array(m.Var,(8,19,5),lb=0)
Altz = m.Array(m.Var,(8,19))
for i in range(8):
for j in range(19):
Altz[i,j]=m.Intermediate(m.sum(Alt[i,j,:]*z[i,j,:]))
for k in range(5):
Ecz[i,j,k]=m.Intermediate(Ec[i,j,k]*z[i,j,k])
# Constraints
m.Equation(m.sum(m.sum(m.sum(Ecz)))<=m.sum(m.sum(m.sum(Ecy))))
for j in range(19):
for k in range(5):
m.Equation(m.sum(z[:,j,k])==x[j,k])
for i in range(8):
for k in range(5):
m.Equation(m.sum(z[i,:,k])==m.sum(y[i,:,k]))
for i in range(8):
for j in range(19):
m.Equation(m.sum((C[i]/N[j]-x[j,:])*z[i,j,:])>=0)
# Objective: to minimize the quantity of courses allocated in different locations
# Example: with the solution y, I have 12 courses in different locations in the periods
# print(sum([sum(Alt[i,j,:]*y[i,j,:])**2 for j in range(19) for i in range(8)])/2)
for i in range(8):
for j in range(19):
m.Obj(Altz[i,j]**2/2)
# Options and final results
m.options.SOLVER=1
m.options.IMODE=2
m.solve()
print(z)
print(m.options.OBJFCNVAL)
注意：我的原题有20节课，171门课，18个地点。
使用 z[i,j,k].value = y[i,j,k] 给出 z 的初步猜测。使用 z[i,j,k] = y[i,j,k] 将 z 条目重新定义为浮点数而不是 gekko 变量类型。
另一个问题是变量 Ecz 和 Altz 被定义为变量 m.Var，然后被覆盖为中间体。相反，尝试分配它们并将它们分配为中间体：
Ecz = np.empty((8,19,5),dtype=object)
Altz = np.empty((8,19),dtype=object)
使用flatten()简化3维数组所有元素的求和。
m.Equation(m.sum(Ecz.flatten())<=sum(Ecy.flatten()))
可以将常量数组定义为 numpy 数组，以避免 Gekko 进行额外的符号处理。这加快了模型编译时间，但对最终解决方案没有影响。
x = np.empty((19,5))
y = np.empty((8,19,5))
N = np.empty((19))
C = np.empty((8))
Ec = np.empty((8,19,5))
Ecy = np.empty((8,19,5))
Alt = np.empty((8,19,5))
为了优化，IMODE应该是3。 IMODE=2 用于参数回归。 IMODE=2 也应该适用于此问题，但 3 是正确的选项，因为您没有尝试适应数据。
m.options.IMODE=3
尝试使用 IPOPT 获得初始非整数解，然后使用 APOPT 找到整数解。
m.solver_options = ['minlp_gap_tol 1.0e-2',\
'minlp_maximum_iterations 10000',\
'minlp_max_iter_with_int_sol 500',\
'minlp_branch_method 1']
混合整数非线性规划 (MINLP) 问题的求解可能具有挑战性，因此您可能需要使用一些 solver options 来加快求解速度。尝试 minlp_branch_method 1 帮助求解器找到初始整数解，以更好地进行 p运行ing。如果可以使用次优解决方案，间隙容差还可以帮助加快解决方案的速度。下面是完整的脚本。考虑在本地使用 remote=False 到 运行 而不是使用 public 服务器，尤其是对于大型优化问题。
import numpy as np
from gekko import GEKKO
# Input data
# Schedule of periods and courses
sched = np.array([ [0, 1, 0, 0, 1], [0, 0, 1, 1, 0], [0, 0, 1, 1, 0], \
[0, 0, 0, 0, 1], [1, 0, 0, 0, 1], [0, 0, 0, 1, 1], [0, 1, 1, 0, 0], \
[1, 0, 0, 1, 0], [0, 1, 0, 0, 1], [1, 1, 0, 0, 0], [0, 1, 1, 0, 0], \
[0, 1, 1, 0, 0], [1, 0, 0, 1, 0], [1, 0, 0, 1, 0], [0, 0, 1, 0, 1], \
[1, 0, 1, 0, 0], [0, 1, 0, 1, 0], [0, 0, 1, 1, 0], [0, 1, 0, 0, 1] ], dtype=np.int64)
# Initial allocation of all periods, courses and locations
alloc=np.array([0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,1,1,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,\
0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], dtype=np.int64)
# Number of students enrolled in each course
enrol = np.array([ 60, 60, 60, 40, 40, 110, 120, 50, 60, 55, 50, \
55, 40, 64, 72, 50, 50, 55, 55], dtype=np.float64)
# Capacity of each location (classroom)
capac = np.array([ 60, 60, 120, 60, 80, 60, 60, 65], dtype=np.float64)
# Total costs of using each location
costs = np.array([ 9017.12, 9017.12, 12050.24, 9017.12, 9413.68, 9017.12, \
9017.12, 9188.96 ])
# Estimated cost of each location by period and student
ecost = np.repeat(np.array([[costs[i]*pow(enrol[j]*5,-1) for j in range(19)] for i in range(8)]), 5)
# The model construction
m = GEKKO(remote=True)
# Constant arrays
x = np.empty((19,5))
y = np.empty((8,19,5))
N = np.empty((19))
C = np.empty((8))
Ec = np.empty((8,19,5))
Ecy = np.empty((8,19,5))
Alt = np.empty((8,19,5))
for k in range(5):
for j in range(19):
N[j] = enrol[j]
x[j,k] = sched[j,k]
for i in range(8):
C[i] = capac[i]
Ec[i,j,k] = ecost[k+j*5+i*19*5]
y[i,j,k] = alloc[k+j*5+i*19*5]
Ecy[i,j,k] = Ec[i,j,k]*y[i,j,k]
if sched[j,k]==1:
Alt[i,j,np.where(sched[j,:]==1)[0][0]]=-sched[j,k]*(1-sum(sched[j,:]))
if sum(sched[j,:])==2:
Alt[i,j,np.where(sched[j,:]==1)[0][1]]=sched[j,k]*(1-sum(sched[j,:]))
else:
Alt[i,j,k]=0
# Initialize the variable z with the initial value y:
# These commented approaches produce the error.
z = m.Array(m.Var,(8,19,5),lb=0,ub=1,integer=True)
for i in range(8):
for j in range(19):
for k in range(5):
z[i,j,k].value = y[i,j,k]
# nor
#z = m.Array(m.Var,(8,19,5),value=y,lb=0,ub=1,integer=True)
# Intermediate equations
Ecz = np.empty((8,19,5),dtype=object)
Altz = np.empty((8,19),dtype=object)
for i in range(8):
for j in range(19):
Altz[i,j]=m.Intermediate(m.sum(Alt[i,j,:]*z[i,j,:]))
for k in range(5):
Ecz[i,j,k]=m.Intermediate(Ec[i,j,k]*z[i,j,k])
# Constraints
m.Equation(m.sum(Ecz.flatten())<=sum(Ecy.flatten()))
for j in range(19):
for k in range(5):
m.Equation(m.sum(z[:,j,k])==x[j,k])
for i in range(8):
for k in range(5):
m.Equation(m.sum(z[i,:,k])==m.sum(y[i,:,k]))
for i in range(8):
for j in range(19):
m.Equation(m.sum((C[i]/N[j]-x[j,:])*z[i,j,:])>=0)
# Objective: to minimize the quantity of courses allocated in different locations
# Example: with the solution y, I have 12 courses in different locations in the periods
# print(sum([sum(Alt[i,j,:]*y[i,j,:])**2 for j in range(19) for i in range(8)])/2)
for i in range(8):
for j in range(19):
m.Obj(Altz[i,j]**2/2)
# Options and final results
m.options.IMODE=3
# Initialize with IPOPT
m.options.SOLVER=3
m.solve()
# Integer solution with APOPT
m.options.SOLVER=1
m.solver_options = ['minlp_gap_tol 1.0e-2',\
'minlp_maximum_iterations 10000',\
'minlp_max_iter_with_int_sol 500',\
'minlp_branch_method 1']
m.solve()
print(z)
print(m.options.OBJFCNVAL)
我正在尝试根据三维 主时间表模型 求解一个步骤，其中涉及时段 (5)、课程 (19) 和地点 (8)。
所以我在 Gekko 中用 3D 数组初始化这些变量时遇到了问题。如果没有这个初始化，算法在超过 15 分钟 运行 和 1000 次迭代后不会收敛。
当我尝试初始化时，出现此错误：
” 引发异常（响应） 异常：@error：方程定义 没有等式 (=) 或不等式 (>,<) 的等式 真的 停止... “
我该如何解决这个问题？遵循我的代码版本：
import numpy as np
from gekko import GEKKO
# Input data
# Schedule of periods and courses
sched = np.array([ [0, 1, 0, 0, 1], [0, 0, 1, 1, 0], [0, 0, 1, 1, 0], \
[0, 0, 0, 0, 1], [1, 0, 0, 0, 1], [0, 0, 0, 1, 1], [0, 1, 1, 0, 0], \
[1, 0, 0, 1, 0], [0, 1, 0, 0, 1], [1, 1, 0, 0, 0], [0, 1, 1, 0, 0], \
[0, 1, 1, 0, 0], [1, 0, 0, 1, 0], [1, 0, 0, 1, 0], [0, 0, 1, 0, 1], \
[1, 0, 1, 0, 0], [0, 1, 0, 1, 0], [0, 0, 1, 1, 0], [0, 1, 0, 0, 1] ], dtype=np.int64)
# Initial allocation of all periods, courses and locations
alloc=np.array([0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,1,1,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,\
0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], dtype=np.int64)
# Number of students enrolled in each course
enrol = np.array([ 60, 60, 60, 40, 40, 110, 120, 50, 60, 55, 50, \
55, 40, 64, 72, 50, 50, 55, 55], dtype=np.float64)
# Capacity of each location (classroom)
capac = np.array([ 60, 60, 120, 60, 80, 60, 60, 65], dtype=np.float64)
# Total costs of using each location
costs = np.array([ 9017.12, 9017.12, 12050.24, 9017.12, 9413.68, 9017.12, \
9017.12, 9188.96 ])
# Estimated cost of each location by period and student
ecost = np.repeat(np.array([[costs[i]*pow(enrol[j]*5,-1) for j in range(19)] for i in range(8)]), 5)
# The model construction
m = GEKKO()
# Constant arrays
x = m.Array(m.Const,(19,5))
y = m.Array(m.Const,(8,19,5))
N = m.Array(m.Const,(19))
C = m.Array(m.Const,(8))
Ec = m.Array(m.Const,(8,19,5))
Ecy = m.Array(m.Const,(8,19,5))
Alt = m.Array(m.Const,(8,19,5))
for k in range(5):
for j in range(19):
N[j] = enrol[j]
x[j,k] = sched[j,k]
for i in range(8):
C[i] = capac[i]
Ec[i,j,k] = ecost[k+j*5+i*19*5]
y[i,j,k] = alloc[k+j*5+i*19*5]
Ecy[i,j,k] = Ec[i,j,k]*y[i,j,k]
if sched[j,k]==1:
Alt[i,j,np.where(sched[j,:]==1)[0][0]]=-sched[j,k]*(1-sum(sched[j,:]))
if sum(sched[j,:])==2:
Alt[i,j,np.where(sched[j,:]==1)[0][1]]=sched[j,k]*(1-sum(sched[j,:]))
else:
Alt[i,j,k]=0
# Initialize the variable z with the initial value y:
# These commented approaches produce the error.
z = m.Array(m.Var,(8,19,5),lb=0,ub=1,integer=True)
#for i in range(8):
# for j in range(19):
# for k in range(5):
# z[i,j,k] = y[i,j,k]
# nor
#z = m.Array(m.Var,(8,19,5),value=y,lb=0,ub=1,integer=True)
# Intermediate equations
Ecz = m.Array(m.Var,(8,19,5),lb=0)
Altz = m.Array(m.Var,(8,19))
for i in range(8):
for j in range(19):
Altz[i,j]=m.Intermediate(m.sum(Alt[i,j,:]*z[i,j,:]))
for k in range(5):
Ecz[i,j,k]=m.Intermediate(Ec[i,j,k]*z[i,j,k])
# Constraints
m.Equation(m.sum(m.sum(m.sum(Ecz)))<=m.sum(m.sum(m.sum(Ecy))))
for j in range(19):
for k in range(5):
m.Equation(m.sum(z[:,j,k])==x[j,k])
for i in range(8):
for k in range(5):
m.Equation(m.sum(z[i,:,k])==m.sum(y[i,:,k]))
for i in range(8):
for j in range(19):
m.Equation(m.sum((C[i]/N[j]-x[j,:])*z[i,j,:])>=0)
# Objective: to minimize the quantity of courses allocated in different locations
# Example: with the solution y, I have 12 courses in different locations in the periods
# print(sum([sum(Alt[i,j,:]*y[i,j,:])**2 for j in range(19) for i in range(8)])/2)
for i in range(8):
for j in range(19):
m.Obj(Altz[i,j]**2/2)
# Options and final results
m.options.SOLVER=1
m.options.IMODE=2
m.solve()
print(z)
print(m.options.OBJFCNVAL)
注意：我的原题有20节课，171门课，18个地点。
使用 z[i,j,k].value = y[i,j,k] 给出 z 的初步猜测。使用 z[i,j,k] = y[i,j,k] 将 z 条目重新定义为浮点数而不是 gekko 变量类型。
另一个问题是变量 Ecz 和 Altz 被定义为变量 m.Var，然后被覆盖为中间体。相反，尝试分配它们并将它们分配为中间体：
Ecz = np.empty((8,19,5),dtype=object)
Altz = np.empty((8,19),dtype=object)
使用flatten()简化3维数组所有元素的求和。
m.Equation(m.sum(Ecz.flatten())<=sum(Ecy.flatten()))
可以将常量数组定义为 numpy 数组，以避免 Gekko 进行额外的符号处理。这加快了模型编译时间，但对最终解决方案没有影响。
x = np.empty((19,5))
y = np.empty((8,19,5))
N = np.empty((19))
C = np.empty((8))
Ec = np.empty((8,19,5))
Ecy = np.empty((8,19,5))
Alt = np.empty((8,19,5))
为了优化，IMODE应该是3。 IMODE=2 用于参数回归。 IMODE=2 也应该适用于此问题，但 3 是正确的选项，因为您没有尝试适应数据。
m.options.IMODE=3
尝试使用 IPOPT 获得初始非整数解，然后使用 APOPT 找到整数解。
m.solver_options = ['minlp_gap_tol 1.0e-2',\
'minlp_maximum_iterations 10000',\
'minlp_max_iter_with_int_sol 500',\
'minlp_branch_method 1']
混合整数非线性规划 (MINLP) 问题的求解可能具有挑战性，因此您可能需要使用一些 solver options 来加快求解速度。尝试 minlp_branch_method 1 帮助求解器找到初始整数解，以更好地进行 p运行ing。如果可以使用次优解决方案，间隙容差还可以帮助加快解决方案的速度。下面是完整的脚本。考虑在本地使用 remote=False 到 运行 而不是使用 public 服务器，尤其是对于大型优化问题。
import numpy as np
from gekko import GEKKO
# Input data
# Schedule of periods and courses
sched = np.array([ [0, 1, 0, 0, 1], [0, 0, 1, 1, 0], [0, 0, 1, 1, 0], \
[0, 0, 0, 0, 1], [1, 0, 0, 0, 1], [0, 0, 0, 1, 1], [0, 1, 1, 0, 0], \
[1, 0, 0, 1, 0], [0, 1, 0, 0, 1], [1, 1, 0, 0, 0], [0, 1, 1, 0, 0], \
[0, 1, 1, 0, 0], [1, 0, 0, 1, 0], [1, 0, 0, 1, 0], [0, 0, 1, 0, 1], \
[1, 0, 1, 0, 0], [0, 1, 0, 1, 0], [0, 0, 1, 1, 0], [0, 1, 0, 0, 1] ], dtype=np.int64)
# Initial allocation of all periods, courses and locations
alloc=np.array([0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,1,1,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,\
0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], dtype=np.int64)
# Number of students enrolled in each course
enrol = np.array([ 60, 60, 60, 40, 40, 110, 120, 50, 60, 55, 50, \
55, 40, 64, 72, 50, 50, 55, 55], dtype=np.float64)
# Capacity of each location (classroom)
capac = np.array([ 60, 60, 120, 60, 80, 60, 60, 65], dtype=np.float64)
# Total costs of using each location
costs = np.array([ 9017.12, 9017.12, 12050.24, 9017.12, 9413.68, 9017.12, \
9017.12, 9188.96 ])
# Estimated cost of each location by period and student
ecost = np.repeat(np.array([[costs[i]*pow(enrol[j]*5,-1) for j in range(19)] for i in range(8)]), 5)
# The model construction
m = GEKKO(remote=True)
# Constant arrays
x = np.empty((19,5))
y = np.empty((8,19,5))
N = np.empty((19))
C = np.empty((8))
Ec = np.empty((8,19,5))
Ecy = np.empty((8,19,5))
Alt = np.empty((8,19,5))
for k in range(5):
for j in range(19):
N[j] = enrol[j]
x[j,k] = sched[j,k]
for i in range(8):
C[i] = capac[i]
Ec[i,j,k] = ecost[k+j*5+i*19*5]
y[i,j,k] = alloc[k+j*5+i*19*5]
Ecy[i,j,k] = Ec[i,j,k]*y[i,j,k]
if sched[j,k]==1:
Alt[i,j,np.where(sched[j,:]==1)[0][0]]=-sched[j,k]*(1-sum(sched[j,:]))
if sum(sched[j,:])==2:
Alt[i,j,np.where(sched[j,:]==1)[0][1]]=sched[j,k]*(1-sum(sched[j,:]))
else:
Alt[i,j,k]=0
# Initialize the variable z with the initial value y:
# These commented approaches produce the error.
z = m.Array(m.Var,(8,19,5),lb=0,ub=1,integer=True)
for i in range(8):
for j in range(19):
for k in range(5):
z[i,j,k].value = y[i,j,k]
# nor
#z = m.Array(m.Var,(8,19,5),value=y,lb=0,ub=1,integer=True)
# Intermediate equations
Ecz = np.empty((8,19,5),dtype=object)
Altz = np.empty((8,19),dtype=object)
for i in range(8):
for j in range(19):
Altz[i,j]=m.Intermediate(m.sum(Alt[i,j,:]*z[i,j,:]))
for k in range(5):
Ecz[i,j,k]=m.Intermediate(Ec[i,j,k]*z[i,j,k])
# Constraints
m.Equation(m.sum(Ecz.flatten())<=sum(Ecy.flatten()))
for j in range(19):
for k in range(5):
m.Equation(m.sum(z[:,j,k])==x[j,k])
for i in range(8):
for k in range(5):
m.Equation(m.sum(z[i,:,k])==m.sum(y[i,:,k]))
for i in range(8):
for j in range(19):
m.Equation(m.sum((C[i]/N[j]-x[j,:])*z[i,j,:])>=0)
# Objective: to minimize the quantity of courses allocated in different locations
# Example: with the solution y, I have 12 courses in different locations in the periods
# print(sum([sum(Alt[i,j,:]*y[i,j,:])**2 for j in range(19) for i in range(8)])/2)
for i in range(8):
for j in range(19):
m.Obj(Altz[i,j]**2/2)
# Options and final results
m.options.IMODE=3
# Initialize with IPOPT
m.options.SOLVER=3
m.solve()
# Integer solution with APOPT
m.options.SOLVER=1
m.solver_options = ['minlp_gap_tol 1.0e-2',\
'minlp_maximum_iterations 10000',\
'minlp_max_iter_with_int_sol 500',\
'minlp_branch_method 1']
m.solve()
print(z)
print(m.options.OBJFCNVAL)