Gekko 返回错误的成功解决方案 [英] Gekko returning incorrect successful solution

问题描述

以下代码返回Successful Solution 目标:0..但这不是最优解.最优解是-6.通过阅读其他问题，我认为在目标函数中使用非 Gekko 函数存在问题，但我使用的唯一非 Gekko 函数是 np.matmul.np.matmul 是否适用于 Gekko 数组?注意 sigma_post 是一个 nxn numpy 单位矩阵.

The following code returns Successful Solution Objective: 0.. But it is not the optimal solution. The optimal solution is -6. From reading other issues I believe it's a problem with using non-Gekko functions in the objective function but the only non-Gekko function I use is np.matmul. Does np.matmul work with gekko arrays? Note sigma_post is an nxn numpy identity matrix.

m = GEKKO(remote=False)
m.options.max_iter=1000
#m.options.ipopt_integer_tol=1
#m.solver_options = ['minlp_integer_tol 50']
#m.solver_options = ['minlp_max_iter_with_int_sol 1000',
#                    'minlp_branch_method 1']
                    

N = 2
b = m.Array(m.Var,(N,n), lb=0, ub=1, integer=True)

for i in range(N):
    for j in range(n):
        if j in [qb_index_range[0], rb_index_range[0], wr_index_range[0]]:
            b[i][j].value = 1
        else:
            b[i][j].value = 0

print('b: ', b)

# CONSTRAINT: Each Lineup must be less than budget
z = np.array([None]*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]*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]*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 with two lineups
sigma_1 = np.array([[None]*N for i in range(N)])
sig = np.matmul(np.matmul(b, sigma_post), b.T)

for i in range(N):
    for j in range(N):
        sigma_1[i][j] = m.Intermediate(sig[i][j])
        
        
m.Obj(-(sigma_1[0][0] + sigma_1[1][1]- 2*sigma_1[1][0]))

m.options.SOLVER = 1 

m.solve(debug=0)

为了透明，理想情况下我关心的目标函数如下，但上面详述的更简单的目标函数会导致问题，所以我决定从那里开始故障排除.下面的目标抛出 Warning: no more possible trial points and no integer solution Maximum iterations 对于 mu 的某些值，但 mu 不存在于约束.非常感谢您的建议！

To be transparent, ideally the objective function I care about is below but the simpler objective function detailed above is causing issues so I decided to start troubleshooting there. The below objective throws Warning: no more possible trial points and no integer solution Maximum iterations for some values of mu but mu is not present in the constraints. Thanks so much for any advice!

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])

# OBJECTIVE: Maximize 
mu_1 = np.array([None]*N)
for i in range(N):
    mu_1[i] = m.Intermediate(np.matmul(b[i, :], mu))


inter = m.if2(theta(sigma_1, m)-eps, .5*mu_1[0]+.5*mu_1[1], 
             (mu_1[0]*normal_cdf((mu_1[0]-mu_1[1])/theta(sigma_1, m), m) + \
              mu_1[1]*normal_cdf((mu_1[1]-mu_1[0])/theta(sigma_1, m), m) + \
              theta(sigma_1, m)*normal_pdf((mu_1[0]-mu_1[1])/theta(sigma_1, m), m)))

m.Obj(-inter)

推荐答案

使用 np.matmul 或任何其他允许对象而不是数字值的函数都没有问题.需要对象是因为 b 是 Gekko 类型值的数组，需要使用自动微分计算导数.您还可以使用新的 @ 运算符来简化表达式.您的原始问题陈述不完整，缺少许多定义.我添加了一些示例值，以便脚本可以在没有定义错误的情况下运行.以下是帮助重现错误的指南.

There is no problem to use np.matmul or any other function that allows objects instead of only numeric values. Objects are needed because b is an array of Gekko type values that are needed to compute the derivatives with automatic differentiation. You can also use the new @ operator that simplifies the expressions. Your original problem statement was incomplete with many missing definitions. I added a few sample values so that the script can run without definition errors. Here are guidelines to help reproduce the error.

N = 2
n = 3
qb_index_range = [0,2]
rb_index_range = [0,2]
wr_index_range = [0,2]
info_df = pd.DataFrame({'cost':np.ones(n)})
budget = 100
sigma_post = np.random.rand(n,n)

这是一个使用 np.matmul() 的例子，它也可以是点积 np.dot().

Here is an example of using np.matmul() that can also be the dot product np.dot().

sigma_1 = np.matmul(np.matmul(b,sigma_post), b.T)

这也可以用矩阵乘法运算符来写.

This can also be written with the matrix multiplication operator.

sigma_1 = b@sigma_post@b.T

这是完整的脚本.

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

m = GEKKO(remote=False)
m.options.max_iter=1000
                    
N = 2
n = 3
b = m.Array(m.Var,(N,n), lb=0, ub=1, integer=True)
qb_index_range = [0,2]
rb_index_range = [0,2]
wr_index_range = [0,2]
info_df = pd.DataFrame({'cost':np.ones(n)})
budget = 100
sigma_post = np.eye(n)

for i in range(N):
    for j in range(n):
        if j in [qb_index_range[0], rb_index_range[0], wr_index_range[0]]:
            b[i][j].value = 1
        else:
            b[i][j].value = 0

# CONSTRAINT: Each Lineup must be less than budget
z = [None]*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 = [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]*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]*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 with two lineups
#sigma_1 = np.matmul(np.matmul(b,sigma_post), b.T)
sigma_1 = b@sigma_post@b.T
                
m.Maximize(sigma_1[0][0] + sigma_1[1][1]- 2*sigma_1[1][0])

m.options.SOLVER = 1 

m.solve(debug=0,disp=False)

print(b)

这产生了一个成功的解决方案.无法验证正确的解决方案，因为原始问题陈述不完整.

This produces a successful solution. The correct solution cannot be verified because the original problem statement is not complete.

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

这篇关于Gekko 返回错误的成功解决方案的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！