无法让 Gekko 解决中等规模的 MINLP 问题答案

【问题标题】：Cannot get Gekko to work on moderate size MINLP problem无法让 Gekko 解决中等规模的 MINLP 问题
【发布时间】：2021-07-15 03:33:15
【问题描述】：

我有一个问题，希望最大化一个比率目标函数。

对象：sum(Ax)/Bx

对于具有相同维度的矩阵A 和B（随机零和一）。问题是不受约束的。所以基本上我想找到一个由零和一组成的x，使给定的比率最大化。

这是我的问题的最小可重现示例（下面的第 2 项）：

import numpy as np
import pandas as pd
from gekko import GEKKO

# Number of columns and rows for the given matrices
C = 43
R = 20

# Name columns since I will need to work with dataframes
names = ['n'+str(i) for i in range(C)]
names2 = ['d'+str(i) for i in range(C)]

# Toy example of random distribution of zeros and ones
num = np.random.randint(2, size=(R, C))
den = np.random.randint(2, size=(R, C))
df_num = pd.DataFrame(columns=names, data=num)
df_den = pd.DataFrame(columns=names2, data=den)

# Start Gekko model
m = GEKKO(remote=False)
x = m.Array(m.Var,(df_num.shape[1]),lb=0,ub=1,integer=True)
for i in x:
    i.value = np.random.randint(2, size=1)[0]
m.solver_options = ['minlp_as_nlp 1']
ival1 = m.Intermediate(df_num.mul(x).sum(1).sum())
ival2 = m.Intermediate(df_den.mul(x).sum(1).sum())
m.Obj(-np.divide(ival1, ival2))
m.options.SOLVER = 1 # APOPT solver
m.solve(disp=True)

我尝试按照here 和here 以及其他一些地方给出的说明进行操作，因此我使用minlp_as_nlp 1 选项来放松整数限制，并使用中间变量来查看它是否有帮助，但是，我无法得到工作的模型。

基本上会发生以下两个问题：

1.未找到解（x 全为零）

那时我遇到了一个非常小的问题（如上面的问题，但有 42 列或更少，但是，我的实际问题范围从 10k x 5e2 到 7e6 x 5e2）。我希望x 能够找到最终值，但它全为零。这是 APM 打印出来的：

----------------------------------------------------------------
 APMonitor, Version 0.9.2
 APMonitor Optimization Suite
 ----------------------------------------------------------------
 
 
 --------- APM Model Size ------------
 Each time step contains
   Objects      :            0
   Constants    :            0
   Variables    :           42
   Intermediates:            2
   Connections  :            0
   Equations    :            3
   Residuals    :            1
 
 Number of state variables:             42
 Number of total equations: -            0
 Number of slack variables: -            0
 ---------------------------------------
 Degrees of freedom       :             42
 
 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:      0.04 NLPi:   18 Dpth:    0 Lvs:    0 Obj: -2.17E+00 Gap:  0.00E+00
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   4.780000000027940E-002 sec
 Objective      :   -2.16666666666667     
 Successful solution
 ---------------------------------------------------

2。运行时错误

当我将列数增加到 43（如上或更多）时，我收到以下错误（已编辑以适应空间）：

----------------------------------------------------------------
 APMonitor, Version 0.9.2
 APMonitor Optimization Suite
 ----------------------------------------------------------------
 
 
 @error: Max Equation Length
 Error with line number:           48
 i28=(((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((0)*(int_v1))
 +((0)*(int_v2)))+((0)*(int_v3)))+((1)*(int_v4)))+((1)*(int_v5)))+((1)*(int_v6))
 )+((0)*(int_v7)))+((1)*(int_v8)))+((0)*(int_v9)))+((1)*(int_v10)))+((0)*(int_v1
.
.
.
_v31)))+((0)*(int_v32)))+((0)*(int_v33)))+((1)*(int_v34)))+((1)*(int_v35)))+((1
 )*(int_v36)))+((1)*(int_v37)))+((0)*(int_v38)))+((0)*(int_v39)))+((0)*(int_v40)
 ))+((0)*(int_v41)))+((0)*(int_v42)))+((1)*(int_v43))))
 
 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
 
 STOPPING...
 

---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
<ipython-input-50-5885f6e9b4ca> in <module>
      6 m.Obj(-np.divide(ival1, ival2))
      7 m.options.SOLVER = 1
----> 8 m.solve(disp=True)

/anaconda_env/personal/rafaeldasilv/py3/lib/python3.7/site-packages/gekko/gekko.py in solve(self, disp, debug, GUI, **kwargs)
   2128                 print("Error:", errs)
   2129             if (debug >= 1) and record_error:
-> 2130                 raise Exception(apm_error)
   2131 
   2132         else: #solve on APM server

Exception: @error: Max Equation Length
 Error with line number:           48
 i28=(((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((0)*(int_v1))
 +((0)*(int_v2)))+((0)*(int_v3)))+((1)*(int_v4)))+((1)*(int_v5)))+((1)*(int_v6))
 )+((0)*(int_v7)))+((1)*(int_v8)))+((0)*(int_v9)))+((1)*(int_v10)))+((0)*(int_v1
 1)))+((1)*(int_v12)))+((1)*(int_v13)))+((0)*(int_v14)))+((1)*(int_v15)))+((1)*(
.
.
.
 )+((1)*(int_v27)))+((0)*(int_v28)))+((1)*(int_v29)))+((0)*(int_v30)))+((0)*(int
 _v31)))+((0)*(int_v32)))+((0)*(int_v33)))+((1)*(int_v34)))+((1)*(int_v35)))+((1
 )*(int_v36)))+((1)*(int_v37)))+((0)*(int_v38)))+((0)*(int_v39)))+((0)*(int_v40)
 ))+((0)*(int_v41)))+((0)*(int_v42)))+((1)*(int_v43))))
 
 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
 
 STOPPING...

有没有人猜到我在这里可能遗漏了什么，或者即使这是正确的方法？谢谢！

编辑：运行建议的代码更改后的解决方案（R=400，C=76）：

----------------------------------------------------------------
 APMonitor, Version 0.9.2
 APMonitor Optimization Suite
 ----------------------------------------------------------------
 
 
 --------- APM Model Size ------------
 Each time step contains
   Objects      :          802
   Constants    :            0
   Variables    :        61678
   Intermediates:            0
   Connections  :        62402
   Equations    :        60801
   Residuals    :        60801
 
 Number of state variables:          61678
 Number of total equations: -        61602
 Number of slack variables: -            0
 ---------------------------------------
 Degrees of freedom       :             76
 
 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:      1.60 NLPi:    3 Dpth:    0 Lvs:    3 Obj: -1.00E+00 Gap:       NaN
--Integer Solution:   0.00E+00 Lowest Leaf:  -1.00E+00 Gap:   1.00E+00
Iter:     2 I:  0 Tm:      0.44 NLPi:    2 Dpth:    1 Lvs:    2 Obj:  0.00E+00 Gap:  1.00E+00
Iter:     3 I:  0 Tm:      1.37 NLPi:    2 Dpth:    1 Lvs:    4 Obj: -1.00E+00 Gap:  1.00E+00
--Integer Solution:  -1.00E+00 Lowest Leaf:  -1.00E+00 Gap:   2.82E-06
Iter:     4 I:  0 Tm:      0.45 NLPi:    2 Dpth:    2 Lvs:    4 Obj: -1.00E+00 Gap:  2.82E-06
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :    4.02930000000015      sec
 Objective      :  -0.999995215333898     
 Successful solution
 ---------------------------------------------------
 

Runtime: 1568.80s

【问题讨论】：

标签： python nonlinear-optimization gekko

【解决方案1】：

尝试使用 Gekko 求和来避免冗长的符号表达式。

m.Maximize(m.sum(Ax)/m.sum(Bx))

没有找到解，因为分母变为零，创建了一个无穷大的目标。这可以通过在分母上添加一个小常数来避免，例如：

m.Maximize(m.sum(Ax)/(m.sum(Bx)+1e-3))

没有m.sum() 的当前版本会生成一个超过 15,000 个字符的目标函数表达式。

我还切换到 A@x 和 B@x 来简化矩阵乘法表达式，但您原来的 Pandas 表达式也可以。

import numpy as np
import pandas as pd
from gekko import GEKKO

# Number of columns and rows for the given matrices
C = 43
R = 20

# Name columns since I will need to work with dataframes
names = ['n'+str(i) for i in range(C)]
names2 = ['d'+str(i) for i in range(C)]

# Toy example of random distribution of zeros and ones
A = np.random.randint(2, size=(R, C))
B = np.random.randint(2, size=(R, C))

# Start Gekko model
m = GEKKO(remote=False)
x = m.Array(m.Var,C,lb=0,ub=1,integer=True)
for i in x:
    i.value = np.random.randint(2, size=1)[0]
m.solver_options = ['minlp_as_nlp 1']

# Calculate Ax and Bx
Ax = A@x
Bx = B@x

# Maximize ratio
m.Maximize(m.sum(Ax)/(m.sum(Bx)+1e-3))

# Solve
m.options.SOLVER = 1 # APOPT solver
m.solve(disp=True)

解决方案：

 Iter    Objective  Convergence
   30 -2.19956E+00  1.30104E-18
   31 -2.19956E+00  6.93889E-18
   32 -2.19956E+00  6.93889E-18
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :  0.13 sec
 Objective      :  -2.1995600879824035
 Successful solution
 ---------------------------------------------------

编辑：替代矩阵乘法

如果有更大的矩阵，则使用 Gekko 求和进行矩阵乘法。

# Calculate Ax and Bx
#Ax = A@x
#Bx = B@x
Ax = [m.sum([A[i,j]*x[j] for j in range(C)]) for i in range(R)]
Bx = [m.sum([A[i,j]*x[j] for j in range(C)]) for i in range(R)]

这可以通过C=43 和R=20 成功求解，但也适用于更大的矩阵。

【讨论】：

谢谢@johnh 我认为为什么以及何时使用 sum 让我感到困惑，因为我认为点积运算可以做到这一点。我对documentation 的看法是，在大多数情况下，使用 sum 是可选的。我已经运行了您的代码，但仍然收到错误 2。我注意到我的 APM 监视器是 0.9.2，而您的 APOPT 是 1.0（它们是否相关？），我的 gekko 版本是 0.2.8，我尝试升级它，但它似乎保持最新。还有什么我可以尝试的吗？
它可能只是选择的随机数 - 这些对于脚本来说不一定是确定性的。也许尝试使用固定的A 和B 矩阵，或者使用固定的种子值使随机化器具有确定性？
我尝试了大约 100 次以查看随机数是否会成为问题，但这并没有帮助。但即便如此，图书馆不应该找到解决方案吗？我想我无法理解为什么矩阵中的稀疏性会影响错误。
是的，它适用于 43 列和 20 行。我还尝试了 100 列和 20 行，并且没有问题。我将原始问题的大小减少到 11k 行 x 75 列，它目前正在运行，我会及时通知你，谢谢 John。有没有什么地方可以找到 m.sum 的这些功能？
我第一次尝试运行完整的问题大约花了 3 天时间，但我不小心中断了它。我已经运行了几个较小的问题来估计完整问题的运行时间（运行时间基本上是一个几何序列，行数每增加一倍，运行时间就会增加 4 倍）所以对于我原来的问题，答案会来大约 18 天后。