GEKKO 异常 : @error: Max Equation Length (Number of variables greater than 100k)

Question

我需要对 100k 到 500k 变量进行优化，但它给我最大方程长度达到错误。谁能帮我解决这个问题？时间不是限制，只要跑个3-4小时就可以了。

df1 = df_opt.head(100000).copy()

#initialize model
m= GEKKO()
m.options.SOLVER=1

#initialize variable
x =  np.array([m.Var(lb=0,ub=100,integer=True) for i in range(len(df1))])

#constraints
m.Equation(m.sum(x)<=30000)

#objective
responsiveness = np.array([m.Const(i) for i in df1['responsivness'].values])
affinity_score = np.array([m.Const(i) for i in df1['affinity'].values])
cost = np.array([m.Const(i) for i in df1['cost'].values])

expr = np.array([m.log(i) - k * j \
    for i,j,k in zip((1+responsiveness * affinity_score * x),x,cost)])

m.Obj(-(m.sum(expr)))

#optimization
m.solve(disp=False)

Answer 1

创建问题时，有一个 Minimal Example that is complete 很重要。这是一个创建带有 n 行的随机 DataFrame 的修改。

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

n = 10
df1 = pd.DataFrame({'responsivness':np.random.rand(n),\
                    'affinity':np.random.rand(n),\
                    'cost':np.random.rand(n)})
print(df1.head())

#initialize model
m= GEKKO(remote=False)
m.options.SOLVER=1

#initialize variable
x =  np.array([m.Var(lb=0,ub=100,integer=True) for i in range(len(df1))])

#constraints
m.Equation(m.sum(x)<=30000)

#objective
responsiveness = np.array([m.Const(i) for i in df1['responsivness'].values])
affinity_score = np.array([m.Const(i) for i in df1['affinity'].values])
cost = np.array([m.Const(i) for i in df1['cost'].values])

expr = np.array([m.log(i) - k * j \
    for i,j,k in zip((1+responsiveness * affinity_score * x),x,cost)])

m.Obj(-(m.sum(expr)))

#optimization
m.solve(disp=True)

这成功解决了选择随机数的 n=10。

 --------- APM Model Size ------------
 Each time step contains
   Objects      :            0
   Constants    :           30
   Variables    :           11
   Intermediates:            0
   Connections  :            0
   Equations    :            2
   Residuals    :            2
 
 Number of state variables:             11
 Number of total equations: -            1
 Number of slack variables: -            1
 ---------------------------------------
 Degrees of freedom       :              9
 
 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:      0.00 NLPi:   20 Dpth:    0 Lvs:    3 Obj: -1.35E+00 Gap:       NaN
--Integer Solution:  -1.34E+00 Lowest Leaf:  -1.35E+00 Gap:   4.73E-03
Iter:     2 I:  0 Tm:      0.00 NLPi:    2 Dpth:    1 Lvs:    3 Obj: -1.34E+00 Gap:  4.73E-03
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   1.519999999436550E-002 sec
 Objective      :   -1.34078995171088     
 Successful solution
 ---------------------------------------------------

底层模型 gk_model0.apm 可以通过导航到 m.path 或使用 m.open_folder() 来访问。

Model
Constants
    i0 = 0.14255660947333681
    i1 = 0.9112789578520111
    i2 = 0.10526966142004568
    i3 = 0.6255161023214897
    i4 = 0.2434604974789274
    i5 = 0.812768922376058
    i6 = 0.555163868440599
    i7 = 0.7286240480266872
    i8 = 0.39643651685899695
    i9 = 0.4664238475079081
    i10 = 0.588654005219946
    i11 = 0.7807594551372589
    i12 = 0.623910408858981
    i13 = 0.19421798736230456
    i14 = 0.3061420839190525
    i15 = 0.07764492888189267
    i16 = 0.7276569154297892
    i17 = 0.5630014016669598
    i18 = 0.9633171115575193
    i19 = 0.23310692223695684
    i20 = 0.008089496373502647
    i21 = 0.7533529530133879
    i22 = 0.4218710975774087
    i23 = 0.03329287687223692
    i24 = 0.9136665338169284
    i25 = 0.7528330460265494
    i26 = 0.0810779357870034
    i27 = 0.4183140612726107
    i28 = 0.4381547602657835
    i29 = 0.907339329732971
End Constants
Variables
    int_v1 = 0, <= 100, >= 0
    int_v2 = 0, <= 100, >= 0
    int_v3 = 0, <= 100, >= 0
    int_v4 = 0, <= 100, >= 0
    int_v5 = 0, <= 100, >= 0
    int_v6 = 0, <= 100, >= 0
    int_v7 = 0, <= 100, >= 0
    int_v8 = 0, <= 100, >= 0
    int_v9 = 0, <= 100, >= 0
    int_v10 = 0, <= 100, >= 0
End Variables
Equations
    (((((((((int_v1+int_v2)+int_v3)+int_v4)+int_v5)+int_v6)+int_v7)+int_v8)+int_v9)+int_v10)<=30000
    minimize (-((((((((((log((1+((((i0)*(i10)))*(int_v1))))-((i20)*(int_v1)))+(log((1+((((i1)*(i11)))*(int_v2))))-((i21)*(int_v2))))+(log((1+((((i2)*(i12)))*(int_v3))))-((i22)*(int_v3))))+(log((1+((((i3)*(i13)))*(int_v4))))-((i23)*(int_v4))))+(log((1+((((i4)*(i14)))*(int_v5))))-((i24)*(int_v5))))+(log((1+((((i5)*(i15)))*(int_v6))))-((i25)*(int_v6))))+(log((1+((((i6)*(i16)))*(int_v7))))-((i26)*(int_v7))))+(log((1+((((i7)*(i17)))*(int_v8))))-((i27)*(int_v8))))+(log((1+((((i8)*(i18)))*(int_v9))))-((i28)*(int_v9))))+(log((1+((((i9)*(i19)))*(int_v10))))-((i29)*(int_v10)))))
End Equations

End Model

您可以通过将模型修改为以下方式来避免较大的符号表达式字符串:

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

n = 5000
df1 = pd.DataFrame({'responsiveness':np.random.rand(n),\
                    'affinity':np.random.rand(n),\
                    'cost':np.random.rand(n)})
print(df1.head())

#initialize model
m= GEKKO(remote=False)
m.options.SOLVER=1

#initialize variable
x =  np.array([m.Var(lb=0,ub=100,integer=True) for i in range(len(df1))])

#constraints
m.Equation(m.sum(list(x))<=30000)

#objective
responsiveness = df1['responsiveness'].values
affinity_score = df1['affinity'].values
cost           = df1['cost'].values
[m.Maximize(m.log(i) - k * j) \
    for i,j,k in zip((1+responsiveness * affinity_score * x),x,cost)]

#optimization
m.solve(disp=True)

m.open_folder()

这给出了以下的基础模型，它的符号表达式大小不会随着变量数量的增加而增加。

Model
Variables
    int_v1 = 0, <= 100, >= 0
    int_v2 = 0, <= 100, >= 0
    int_v3 = 0, <= 100, >= 0
    int_v4 = 0, <= 100, >= 0
    int_v5 = 0, <= 100, >= 0
    int_v6 = 0, <= 100, >= 0
    int_v7 = 0, <= 100, >= 0
    int_v8 = 0, <= 100, >= 0
    int_v9 = 0, <= 100, >= 0
    int_v10 = 0, <= 100, >= 0
    v11 = 0
End Variables
Equations
    v11<=30000
    maximize (log((1+((0.16283879947305288)*(int_v1))))-((0.365323493448101)*(int_v1)))
    maximize (log((1+((0.3509872155181691)*(int_v2))))-((0.12162206443479917)*(int_v2)))
    maximize (log((1+((0.20134572143617518)*(int_v3))))-((0.47137701674279087)*(int_v3)))
    maximize (log((1+((0.287818142242232)*(int_v4))))-((0.12042554857067544)*(int_v4)))
    maximize (log((1+((0.48997709502894166)*(int_v5))))-((0.21084485862098745)*(int_v5)))
    maximize (log((1+((0.6178277437136291)*(int_v6))))-((0.42602122419609056)*(int_v6)))
    maximize (log((1+((0.13033555293152563)*(int_v7))))-((0.8796057438355324)*(int_v7)))
    maximize (log((1+((0.5002025885707916)*(int_v8))))-((0.9703263879586648)*(int_v8)))
    maximize (log((1+((0.7095523321888202)*(int_v9))))-((0.8498606490337451)*(int_v9)))
    maximize (log((1+((0.6174815809937886)*(int_v10))))-((0.9390903075640681)*(int_v10)))
End Equations
Connections
    int_v1 = sum_1.x[1]
    int_v2 = sum_1.x[2]
    int_v3 = sum_1.x[3]
    int_v4 = sum_1.x[4]
    int_v5 = sum_1.x[5]
    int_v6 = sum_1.x[6]
    int_v7 = sum_1.x[7]
    int_v8 = sum_1.x[8]
    int_v9 = sum_1.x[9]
    int_v10 = sum_1.x[10]
    v11 = sum_1.y
End Connections
Objects
    sum_1 = sum(10)
End Objects

End Model

我修复了 Gekko 中的一个错误，因此您应该能够在 Gekko 的下一个版本中使用 m.Equation(m.sum(x)<=30000) 而不是将 x 转换为列表。此修改现在适用于以前失败的较大模型。我用 n=5000 测试了它。

 Number of state variables:           5002
 Number of total equations: -            2
 Number of slack variables: -            1
 ---------------------------------------
 Degrees of freedom       :           4999
 
 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:    313.38 NLPi:   14 Dpth:    0 Lvs:    3 Obj: -6.05E+02 Gap:       NaN
--Integer Solution:  -6.01E+02 Lowest Leaf:  -6.05E+02 Gap:   6.60E-03
Iter:     2 I:  0 Tm:      0.06 NLPi:    2 Dpth:    1 Lvs:    3 Obj: -6.01E+02 Gap:  6.60E-03
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :    313.461699999985      sec
 Objective      :   -600.648283994940     
 Successful solution
 ---------------------------------------------------

求解时间增加到 313.46 秒。编译模型的处理时间也更长。您可能希望从较小的模型开始并检查它会增加多少计算时间。我还建议您使用 remote=False 在本地而不是在远程服务器上解决。

如果变量更多，整数优化问题所需的时间可能呈指数级增长，因此您需要确保您不会开始一个需要 30 年才能完成的问题。检查这一点的一个好方法是连续解决更大的问题以了解放大情况。