如何在 Gekko 中初始化 3D 数组变量？

Question

我正在尝试从三维主时间表模型解决一个步骤，该模型涉及时间段 (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 个地点。

Answer 1

使用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 帮助求解器找到初始整数解以进行更好的修剪。如果次优解决方案是可以的，间隙容差也有助于加快解决方案。下面是完整的脚本。考虑使用remote=False 在本地运行而不是使用公共服务器，尤其是对于大型优化问题。

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)