Gekko 错误地找不到整数解

如何解决Gekko 错误地找不到整数解

以下 MINLP 问题返回 Warning: no more possible trial points and no integer solution。我认为这是错误的，因为至少 b = [[1,1,1],[1,0]] 是一个可行的点。

当切换到 m.options.soLVER=3 时，我确实得到了一个合理的非整数解。当我使用更简单的目标函数（即 m.Obj(-sig[0][0])）测试完全相同的脚本时，使用 m.options.soLVER=1，Gekko 找到了整数解。

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

info_df = pd.DataFrame({'pos': ['qb','qb','rb','wr','wr'],'cost': [30,40,15,20,30]})
budget = 80
mu_post = np.array([15,16,8,9,10,14,15])
n = mu_post.shape[0]
sigma_post = np.identity(n)
for i in range(n):
    sigma_post[i][i] = i+1

def get_football_position_range(pos,df):
    return (df[df['pos'] == pos].index[0],df[df['pos'] == pos].index[-1])

qb_index_range = get_football_position_range('qb',info_df)
rb_index_range = get_football_position_range('rb',info_df)
wr_index_range = get_football_position_range('wr',info_df)

# Number of lineups
N = 2

pi = 3.14159
eps = 1.0E-6

def normal_cdf(x,m):
    return 1/(1+m.exp(-1.65451*x))
    
def normal_pdf(x,m):
    return (1/((2*pi)**(.5)))*m.exp((x**2)/2)
    
def theta(s,m):
    return m.sqrt(s[0][0]+s[1][1] - 2*s[0][1])

#################################################
#Integer Optimization Program
m = GEKKO(remote=False)

b = m.Array(m.Var,(N,n),lb=0,ub=1,integer=True,value = 1e-3)

# CONSTRAINT: Each Lineup must be less than budget
z = np.array([None for i in range(N)])
for i in range(N):
    z[i] = m.Intermediate(sum(b[i,:]*list(info_df['cost'])))
    
m.Equations([z[i] <= budget for i in range(N)])


# CONSTRAINT: Each Lineup has one QB
z_1 = np.array([None]*N)
for i in range(N):
    z_1[i] = m.Intermediate(sum(b[i,qb_index_range[0]: qb_index_range[1]+1]))

m.Equations([z_1[i] == 1 for i in range(N)])


# CONSTRAINT: Each Lineup has one RB
z_2 = np.array([None for i in range(N)])
for i in range(N):
    z_2[i] = m.Intermediate(sum(b[i,rb_index_range[0]: rb_index_range[1]+1]))

m.Equations([z_2[i] == 1 for i in range(N)])


# CONSTRAINT: Each Lineup has one WR
z_3 = np.array([None for i in range(N)])
for i in range(N):
    z_3[i] = m.Intermediate(sum(b[i,wr_index_range[0]: wr_index_range[1]+1]))

m.Equations([z_3[i] == 1 for i in range(N)])



# OBJECTIVE: Maximize 
mu = b@mu_post
sig = b@sigma_post@b.T


inter = m.if3(theta(sig,m)-eps,.5*mu[0]+.5*mu[1],(mu[0]*normal_cdf((mu[0]-mu[1])/theta(sig,m),m) + \
              mu[1]*normal_cdf((mu[1]-mu[0])/theta(sig,m) + \
              theta(sig,m)*normal_pdf((mu[0]-mu[1])/theta(sig,m)))

m.Obj(-inter)
m.options.soLVER = 1
m.solve(debug=0,disp=True)

请注意，这是对 Gekko returning incorrect successful solution 的后续问题，未完全指定。

解决方法

问题出在求解器以初始猜测 1e-3 开始时。以下是获得成功解决方案的两种方法。

使用您建议的初始猜测：

bv = np.array([[1,1,1],[1,0]])
for i in range(N):
    for j in range(n):
        b[i,j].value = bv[i,j]

然而，这通常不切实际，尤其是对于大问题。

先用 b 初始化 value=1 并用 IPOPT 求解，然后用 APOPT 求解。

m.Maximize(inter)

m.options.SOLVER = 3
m.solve(debug=0,disp=True)

m.options.SOLVER = 1
m.solve(debug=0,disp=True)
print(b)

这给出了一个使用 IPOPT（非整数解）的成功解法，它改进了 IPOPT（整数解法）求解器寻找解法。

 Number of state variables:    20
 Number of total equations: -  11
 Number of slack variables: -  4
 ---------------------------------------
 Degrees of freedom       :    5
 
 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:      0.06 NLPi:   21 Dpth:    0 Lvs:    3 Obj: -3.90E+01 Gap:       NaN
Iter:     2 I: -9 Tm:      2.44 NLPi:  251 Dpth:    1 Lvs:    2 Obj: -3.90E+01 Gap:       NaN
Iter:     3 I:  0 Tm:      0.01 NLPi:    7 Dpth:    1 Lvs:    3 Obj: -3.90E+01 Gap:       NaN
Iter:     4 I:  0 Tm:      0.02 NLPi:   11 Dpth:    1 Lvs:    4 Obj: -3.90E+01 Gap:       NaN
Iter:     5 I:-11 Tm:      0.00 NLPi:    1 Dpth:    2 Lvs:    3 Obj: -3.90E+01 Gap:       NaN
--Integer Solution:  -3.90E+01 Lowest Leaf:  -3.90E+01 Gap:   0.00E+00
Iter:     6 I:  0 Tm:      0.00 NLPi:    1 Dpth:    2 Lvs:    3 Obj: -3.90E+01 Gap:  0.00E+00
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :  2.6365 sec
 Objective      :  -39.000000008290634
 Successful solution
 ---------------------------------------------------
 

[[[0.0] [1.0] [0.0] [1.0] [0.0] [1.0] [0.0]]
 [[0.0] [1.0] [0.0] [1.0] [0.0] [1.0] [0.0]]]

完整代码如下：

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

info_df = pd.DataFrame({'pos': ['qb','qb','rb','wr','wr'],'cost': [30,40,15,20,30]})
budget = 80
mu_post = np.array([15,16,8,9,10,14,15])
n = mu_post.shape[0]
sigma_post = np.identity(n)
for i in range(n):
    sigma_post[i][i] = i+1

def get_football_position_range(pos,df):
    return (df[df['pos'] == pos].index[0],df[df['pos'] == pos].index[-1])

qb_index_range = get_football_position_range('qb',info_df)
rb_index_range = get_football_position_range('rb',info_df)
wr_index_range = get_football_position_range('wr',info_df)

# Number of lineups
N = 2

pi = 3.14159
eps = 1.0E-6

def normal_cdf(x,m):
    return 1/(1+m.exp(-1.65451*x))
    
def normal_pdf(x,m):
    return (1/((2*pi)**(.5)))*m.exp((x**2)/2)
    
def theta(s,m):
    return m.sqrt(s[0][0]+s[1][1] - 2*s[0][1])

#################################################
#Integer Optimization Program
m = GEKKO(remote=False)

b = m.Array(m.Var,(N,n),lb=0,ub=1,integer=True,value=1)
#bv = np.array([[1,0]])
#for i in range(N):
#    for j in range(n):
#        b[i,j]

# CONSTRAINT: Each Lineup must be less than budget
z = np.array([None for i in range(N)])
for i in range(N):
    z[i] = m.Intermediate(sum(b[i,:]*list(info_df['cost'])))
    
m.Equations([z[i] <= budget for i in range(N)])


# CONSTRAINT: Each Lineup has one QB
z_1 = np.array([None]*N)
for i in range(N):
    z_1[i] = m.Intermediate(sum(b[i,qb_index_range[0]: qb_index_range[1]+1]))

m.Equations([z_1[i] == 1 for i in range(N)])


# CONSTRAINT: Each Lineup has one RB
z_2 = np.array([None for i in range(N)])
for i in range(N):
    z_2[i] = m.Intermediate(sum(b[i,rb_index_range[0]: rb_index_range[1]+1]))

m.Equations([z_2[i] == 1 for i in range(N)])


# CONSTRAINT: Each Lineup has one WR
z_3 = np.array([None for i in range(N)])
for i in range(N):
    z_3[i] = m.Intermediate(sum(b[i,wr_index_range[0]: wr_index_range[1]+1]))

m.Equations([z_3[i] == 1 for i in range(N)])

# OBJECTIVE: Maximize 
mu = b@mu_post
sig = b@sigma_post@b.T

inter = m.if3(theta(sig,m)-eps,.5*mu[0]+.5*mu[1],(mu[0]*normal_cdf((mu[0]-mu[1])/theta(sig,m),m) + \
              mu[1]*normal_cdf((mu[1]-mu[0])/theta(sig,m) + \
              theta(sig,m)*normal_pdf((mu[0]-mu[1])/theta(sig,m)))

m.Maximize(inter)
m.options.SOLVER = 3
m.solve(debug=0,disp=True)
m.options.SOLVER = 1
m.solve(debug=0,disp=True)
print(b)