问题描述
对象：sum(Ax)/Bx
对于具有相同维度的矩阵 A 和 B（随机零和一）。问题是不受约束的。所以基本上我想找到一个由 0 和 1 组成的 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,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
解决方法
尝试使用 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,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 求解成功，但也适用于更大的矩阵。