python - 如何修复 Python Gekko 最大方程长度错误

Question

我的代码适用于小变量。但是当我在 128*128 变量数组中执行操作时，出现以下错误:

APM model error: string > 15000 characters Consider breaking up the line into multiple equations The may also be due to only using newline character CR instead of CR LF (for Windows) or LF (for MacOS/Linux) To fix this problem, save APM file with appropriate newline characters

我不知道如何修复这个错误。

假设 gekko 无法像那样运行大数组变量。

我把代码放在下面。

from gekko import GEKKO
import numpy as np

m = GEKKO()
m.options.SOLVER=1
# optional solver settings with APOPT
m.solver_options = ['minlp_maximum_iterations max', \
                    # minlp iterations with integer solution
                    'minlp_max_iter_with_int_sol max', \
                    # treat minlp as nlp
                    'minlp_as_nlp min', \
                    # nlp sub-problem max iterations
                    'nlp_maximum_iterations 5000', \
                    # 1 = depth first, 2 = breadth first
                    'minlp_branch_method 1', \
                    # maximum deviation from whole number
                    'minlp_integer_tol 0.0005', \
                    # covergence tolerance
                    'minlp_gap_tol 0.0001' ,\
                    'max_cpu_time 10min']

n = 128 
p = 8 
q = 16 


x = m.Array(m.Var,(n,n),integer=True)

for i in range(n):
    for j in range(n):
        x[i,j].value = 1
        x[i,j].lower = 0
        x[i,j].upper = 1


s = 60

a = m.Array(m.Param,(x.shape[0],))
a_ = [172,282,10, 264, 287, 442, 393, 428, 484, 444, 344, 250, 293,  34, 473, 110, 338,  62,
 250, 205,  81, 336, 249, 199, 328, 447, 408,  82, 357, 334, 181, 133, 378,  79, 292, 103,
 494, 382,  10, 477, 237, 267, 337, 313, 395, 110, 114, 381,  52, 232, 457,  69, 167, 344,
 363, 284, 136 ,240, 258, 449, 119, 317, 370, 404, 197, 102, 428, 238, 321, 103, 282,  37,
  41,  86, 496  ,31, 148, 245,  78, 219,  37, 115, 378, 129,  37, 448, 415, 368, 238, 406,
 408, 100, 112 ,289 ,205, 329, 338, 396, 494, 145, 355,  45,   5, 220, 411, 167,  85, 369,
 426, 442, 406 ,217,  57, 176,  12, 368, 444, 227, 337,  63, 267, 216, 227, 182, 408, 116,
 385, 140]
for i in range(len(a)):
    a[i] = a_[i]


O1 = np.array(range(n))
O_= np.transpose(np.nonzero(O1+1))

O = np.zeros((x.shape[0],x.shape[1]),dtype=int)
for i in range(n):
    for j in range(n):
        O[i,O_[i]] = 1


d = m.Array(m.Param,(p,q))
for i in range(p):
    for j in range(q):
        d[i,j] = (round(p/2)-int(i/2))*4-(j%4)

#condition
m.Equation(m.sum(abs(x[:,:]-O[:,:]))<=2*s)

for i in range(n):
    m.Equation(m.sum(x[:,i])==1)

for i in range(n):
    m.Equation(m.sum(x[i,:])==1)

#Objective
m.Obj(np.sum(a.T*np.sum(x*d.reshape(1,n),axis=1)))

#Set global options
m.options.IMODE = 3

#Solve simulation
m.solve(disp=False,debug=True)

Answer 1

如果你让它求解足够长的时间，下面的版本应该可以工作，尽管解决方案不太可能，因为你有 128*128 = 16384 个二进制变量。有 2^16384 个解，如果整数变量的数量大于几百，混合整数非线性规划求解器(例如 APOPT)不太可能快速找到解。它可以解决具有数千个变量的连续变量，但混合整数更加困难，因为它必须使用分支定界方法来解决它们。如果您可以使问题线性化(用松弛变量代替 abs 函数)，您可能想尝试使用 CPLEX 或 Gurobi 等 MILP 求解器。

以下是针对您的模型的一些建议:

• 请勿将 numpy 用于模型表达式。 np.abs 和 np.sum 等函数无法为求解器提供所需的信息。您应该使用 m.abs2(或 m.abs3)和 m.sum。
• 数组 a 和 d 可以是常量的 numpy 数组，而不是 Gekko 参数。
• 您可以将至少一个表达式分解为 Intermediates帮助解决问题。
• 目标函数是 a.T x[:,:] d 的标量结果。您可以使用 a[i] x[i,j] d[j] 中的多个 m.Obj 函数，而不是使用矩阵乘法产生的一个大型表达式。
• 我建议您从一个较小的问题开始，然后随着您对模型和解决方案的信心增强，逐步扩展到更大的问题。

from gekko import GEKKO
import numpy as np

m = GEKKO()
m.options.SOLVER=1
# optional solver settings with APOPT
m.solver_options = ['minlp_maximum_iterations 1000', \
                    # minlp iterations with integer solution
                    'minlp_max_iter_with_int_sol 100', \
                    # treat minlp as nlp
                    'minlp_as_nlp 1', \
                    # nlp sub-problem max iterations
                    'nlp_maximum_iterations 5000', \
                    # 1 = depth first, 2 = breadth first
                    'minlp_branch_method 1', \
                    # maximum deviation from whole number
                    'minlp_integer_tol 0.0005', \
                    # covergence tolerance
                    'minlp_gap_tol 0.0001']

n = 128 
p = 8 
q = 16 

x = m.Array(m.Var,(n,n),integer=True)

for i in range(n):
    for j in range(n):
        x[i,j].value = 1
        x[i,j].lower = 0
        x[i,j].upper = 1


s = 60

# don't need a gekko array here, constants are okay
#a = m.Array(m.Param,(x.shape[0],))
a = np.array([172,282,10, 264, 287, 442, 393, 428, 484, 444, 344, 250, 293,  34, 473, 110, 338,  62,
 250, 205,  81, 336, 249, 199, 328, 447, 408,  82, 357, 334, 181, 133, 378,  79, 292, 103,
 494, 382,  10, 477, 237, 267, 337, 313, 395, 110, 114, 381,  52, 232, 457,  69, 167, 344,
 363, 284, 136 ,240, 258, 449, 119, 317, 370, 404, 197, 102, 428, 238, 321, 103, 282,  37,
  41,  86, 496  ,31, 148, 245,  78, 219,  37, 115, 378, 129,  37, 448, 415, 368, 238, 406,
 408, 100, 112 ,289 ,205, 329, 338, 396, 494, 145, 355,  45,   5, 220, 411, 167,  85, 369,
 426, 442, 406 ,217,  57, 176,  12, 368, 444, 227, 337,  63, 267, 216, 227, 182, 408, 116,
 385, 140])

# use .value to assign a value instead of over-writing the gekko variable
#   with a floating point number
#for i in range(len(a)):
#    a[i].value = a_[i]

O1 = np.array(range(n))
O_= np.transpose(np.nonzero(O1+1))

O = np.zeros((x.shape[0],x.shape[1]),dtype=int)
for i in range(n):
    for j in range(n):
        O[i,O_[i]] = 1

d = np.empty((p,q)) # doesn't need to be a gekko array
for i in range(p):
    for j in range(q):
        d[i,j] = (round(p/2)-int(i/2))*4-(j%4)

#condition - use m.abs2 or m.abs3 for versions of ABS that
# have continuous first and second derivatives
xabs = [[m.Intermediate(m.abs3(x[i,j]-O[i,j])) for j in range(n)] for i in range(n)]
# matrix summation
msum = m.sum([m.sum([xabs[i][j] for i in range(n)]) for j in range(n)])
m.Equation(msum<=2*s)

for i in range(n):
    m.Equation(m.sum(x[:,i])==1)

for i in range(n):
    m.Equation(m.sum(x[i,:])==1)

#Objective - can add multiple objectives
e = d.reshape(n)
for i in range(n):
    for j in range(n):
        m.Obj(a[i]*x[i,j]*e[j])
#m.Obj(m.sum(a.T*m.sum(x*d.reshape(1,n))))

#Set global options
m.options.IMODE = 3

#Solve simulation
m.solve(disp=True,debug=True)