在gekko中添加条件变量导致无解

问题描述 投票：0回答：1

我正在使用gekko来优化某个功能。当我使用像 

m.Obj(0)
这样的虚拟目标函数来测试可行性时，求解器能够找到可行的解决方案。但是，当我添加主要目标函数（在下面注释掉）时，求解器无法找到解决方案。下面是一些功能齐全的测试代码：

# ## Imports
from gekko import GEKKO
import numpy as np

## Set fixed variabls

h_matrix = np.array(
        [
        [119.4, 119.4, 119.4, 119.4, 119.4, 119.4, 111.4, 111.4, 111.4, 111.4],
        [119.4, 119.4, 119.4, 119.4, 119.4, 119.4, 111.4, 111.4, 111.4, 111.4],
        [119.4, 119.4, 119.4, 119.4, 119.4, 111.4, 111.4, 111.4, 111.4, 111.4]
        ]
    )

z_matrix = np.array(
        [
        [383.91, 383.91, 383.91, 383.91, 383.91, 383.91, 254.49, 254.49, 254.49, 254.49],
        [383.91, 383.91, 383.91, 383.91, 383.91, 383.91, 254.49, 254.49, 254.49, 254.49],
        [383.91, 383.91, 383.91, 383.91, 383.91, 254.49, 254.49, 254.49, 254.49, 254.49]
        ]
    )

w = np.array([47.93, 66.37])
h = np.array([12.10, 8.6])
t = np.array([104, 48])


## Reshape the matrices to make calulations easier
h_matrix_reshaped = h_matrix.reshape(30,1)
z_matrix_reshaped = z_matrix.reshape(30,1)


# ## Initialize
# Initialize the model

m = GEKKO(remote=False)

## Fixed variables
h_constants = [m.Param(value=h_matrix_reshaped[i][0]) for i in range(30)]
z_constants = [m.Param(value=z_matrix_reshaped[i][0]) for i in range(30)]

w_base = [m.Param(value=w[i]) for i in range(w.shape[0])]
h_base = [m.Param(value=h[i]) for i in range(h.shape[0])]
t_base = [m.Param(value=t[i]) for i in range(t.shape[0])]

h_cm = m.Param(value=220)
ho_cm = m.Param(value=14.4)
w_kg = m.Param(value=21.2)
s_constraint = m.Param(value=40)

# ### Set up x var (main integer variable)
## Initialize x array

x = np.empty((30, t.shape[0]), dtype=object)
for i in range(30):
    for j in range(t.shape[0]):
        x[i][j] = m.Var(value=0, lb=0, integer=True)
        

# ### Set up Constraints

## Total constraint
for j in range(len(t_base)):
    t_contraints = sum(x[i][j] for i in range(30))
    m.Equation(t_contraints == t_base[j])


## Weight contraints
for i in range(30):
    w_constraints = sum(x[i][j]*w_base[j] for j in range(len(w_base))) + w_kg
    m.Equation(w_constraints <= z_constants[i])

## Height constraints
for i in range(30):
    h_constraints = sum(x[i][j]*h_base[j] for j in range(len(h_base))) + ho_cm
    m.Equation(h_constraints <= h_cm)


## Neighbor constraints
for i in range(9):
    # set neighbor constraints horizontally over first row
    neighbor_constraints_1 = m.abs3((sum(x[i][j]*h_base[j] for j in range(len(h_base))) +
                               h_constants[i]) - (sum(x[i+1][j]*h_base[j] for j in range(len(h_base))) +
                               h_constants[i+1]))
    m.Equation(neighbor_constraints_1 <= s_constraint)
    
for i in range(10,19):  
    # set neighbor constraints horizontally over second row
    neighbor_constraints_2 = m.abs3((sum(x[i][j]*h_base[j] for j in range(len(h_base))) +
                               h_constants[i]) - (sum(x[i+1][j]*h_base[j] for j in range(len(h_base))) +
                               h_constants[i+1]))
    m.Equation(neighbor_constraints_2 <= s_constraint)
    
for i in range(20,29):  
    # set neighbor constraints horizontally over second row
    neighbor_constraints_3 = m.abs3((sum(x[i][j]*h_base[j] for j in range(len(h_base))) +
                               h_constants[i]) - (sum(x[i+1][j]*h_base[j] for j in range(len(h_base))) +
                               h_constants[i+1]))
    m.Equation(neighbor_constraints_3 <= s_constraint)
    
for i in range(10):
    # set neighbor constrainst vertically A with B 
    neighbor_constraints_4 = m.abs3((sum(x[i][j]*h_base[j] for j in range(len(h_base))) +
                                h_constants[i]) - (sum(x[i+10][j]*h_base[j] for j in range(len(h_base))) +
                                h_constants[i+10]))
    m.Equation(neighbor_constraints_4 <= s_constraint)

for i in range(10,20):
    # set neighbor constrainst vertically B with C
    neighbor_constraints_5 = m.abs3((sum(x[i][j]*h_base[j] for j in range(len(h_base))) +
                                h_constants[i]) - (sum(x[i+10][j]*h_base[j] for j in range(len(h_base))) +
                                h_constants[i+10]))
    m.Equation(neighbor_constraints_5 <= s_constraint)
    



# ### Mix Score section below ##################

## Create a binary variable/array b that identifies if x[i][j] is non-zero or not
## We will use the count of these b[i][j] values in our objective function

## Constraint to set b[i][k] = 1 if x[i][k] > 0, and 0 otherwise
## Use if3 to set b directly based on x values
b = np.empty((30, len(t_base)), dtype=object)
epsilon = 1e-2  # Small margin allows floating-point considerations

for i in range(30):
    for j in range(len(t_base)):
        b[i][j] = m.if3(x[i][j] - epsilon, 0, 1)

# # Calculation of count(i) for each row
counts = [m.Intermediate(m.sum(b[i])) for i in range(30)]

# # x_sums for sum of each row in x
x_sums = [m.Intermediate(m.sum(x[i])) for i in range(30)]

### Mix Score section above #############################

# ## Run Solver ##

## Set a dummy objective just to identify solutions that are feasible
m.Obj(0)

# Define the main objective function
# mix_score = [counts[i] / m.max2(x_sums[i], 1e-3) for i in range(30)]
# m.Obj(m.sum(mix_score))  


# Set the solver options
m.options.SOLVER = 1  # APOPT solver for non-linear programs

# Increase max iterations because we don't care much about time
m.solver_options = ['minlp_gap_tol 1.0e-4',\
                    'minlp_maximum_iterations 50000',\
                    'minlp_max_iter_with_int_sol 40000']


## Solve
m.solve(disp=True)

在“混合分数”部分中，您将看到我在其中包含定义为 

b[i][j]
的条件变量，这些变量为 1 或 0，具体取决于 
x[i][j]
的值（如果 
x[i][j]
非零则为 1，否则为 0） 。我想将此二进制变量应用于我在注释掉的目标函数中使用的计数函数

# diversity_score = [counts[i] / m.max2(x_sums[i], 1e-3) for i in range(30)]

# m.Obj(m.sum(diversity_score))

当我实现目标函数时我很惊讶

m.Obj(m.sum(diversity_score))

求解器无法找到解决方案，尽管我知道它可以识别至少一个可行的解决方案（通过运行

m.Obj(0)
）。我认为求解器至少应该为使用 
x[i][j]
时找到的可行 
m.Obj(0)
的目标函数返回一个值，但事实并非如此。如有任何帮助，我们将不胜感激。

python optimization gekko
1个回答
0
投票

尝试将 

m.max3()
与二进制开关功能一起使用，而不是使用 MPCC 版本的开关功能 
m.max2()
。有关基于连续梯度的开关函数的更多信息，请参阅工程优化在线课程。

# Define the main objective function
mix_score = [counts[i] / m.max3(x_sums[i], 1e-3) for i in range(30)]
m.Obj(m.sum(mix_score))

对 

m.max3()
进行更改会产生成功的解决方案。

--Integer Solution:   7.34E+00 Lowest Leaf:   7.33E+00 Gap:   2.49E-04
Iter:  3225 I:  0 Tm:      0.01 NLPi:    5 Dpth:   98 Lvs: 1026 Obj:  7.34E+00 Gap:  2.49E-04
Iter:  3226 I: -1 Tm:      0.02 NLPi:    5 Dpth:  100 Lvs: 1025 Obj:  7.33E+00 Gap:  2.49E-04
Iter:  3227 I: -1 Tm:      0.03 NLPi:    5 Dpth:  100 Lvs: 1024 Obj:  7.33E+00 Gap:  2.49E-04
Iter:  3228 I: -1 Tm:      0.01 NLPi:    2 Dpth:  100 Lvs: 1023 Obj:  7.33E+00 Gap:  2.49E-04
Iter:  3229 I: -1 Tm:      0.03 NLPi:    5 Dpth:  103 Lvs: 1022 Obj:  7.33E+00 Gap:  2.49E-04
Iter:  3230 I:  0 Tm:      0.03 NLPi:    7 Dpth:  103 Lvs: 1021 Obj:  7.34E+00 Gap:  2.49E-04
Iter:  3231 I:  0 Tm:      0.02 NLPi:    5 Dpth:  103 Lvs: 1020 Obj:  7.34E+00 Gap:  2.49E-04
Iter:  3232 I: -1 Tm:      0.04 NLPi:    7 Dpth:   87 Lvs: 1019 Obj:  7.34E+00 Gap:  2.49E-04
Iter:  3233 I:  0 Tm:      0.02 NLPi:    5 Dpth:   87 Lvs: 1021 Obj:  7.34E+00 Gap:  2.49E-04
Iter:  3234 I:  0 Tm:      0.02 NLPi:    6 Dpth:   87 Lvs: 1020 Obj:  7.34E+00 Gap:  2.49E-04
Iter:  3235 I: -1 Tm:      0.03 NLPi:    7 Dpth:   94 Lvs: 1019 Obj:  7.34E+00 Gap:  2.49E-04
Iter:  3236 I: -1 Tm:      0.02 NLPi:    4 Dpth:   94 Lvs: 1018 Obj:  7.34E+00 Gap:  2.49E-04
Iter:  3237 I:  0 Tm:      0.02 NLPi:    5 Dpth:   94 Lvs: 1017 Obj:  7.35E+00 Gap:  2.49E-04
Iter:  3238 I: -1 Tm:      0.05 NLPi:    9 Dpth:   88 Lvs: 1016 Obj:  7.34E+00 Gap:  2.49E-04
Iter:  3239 I:  0 Tm:      0.03 NLPi:    5 Dpth:   88 Lvs: 1015 Obj:  7.34E+00 Gap:  2.49E-04
Iter:  3240 I:  0 Tm:      0.02 NLPi:    5 Dpth:   88 Lvs: 1014 Obj:  7.34E+00 Gap:  2.49E-04
 No additional trial points, returning the best integer solution
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :    110.606499999994      sec
 Objective      :    7.33571428571429     
 Successful solution
 ---------------------------------------------------

这是完整的脚本：

# ## Imports
from gekko import GEKKO
import numpy as np

## Set fixed variabls

h_matrix = np.array(
        [
        [119.4, 119.4, 119.4, 119.4, 119.4, 119.4, 111.4, 111.4, 111.4, 111.4],
        [119.4, 119.4, 119.4, 119.4, 119.4, 119.4, 111.4, 111.4, 111.4, 111.4],
        [119.4, 119.4, 119.4, 119.4, 119.4, 111.4, 111.4, 111.4, 111.4, 111.4]
        ]
    )

z_matrix = np.array(
        [
        [383.91, 383.91, 383.91, 383.91, 383.91, 383.91, 254.49, 254.49, 254.49, 254.49],
        [383.91, 383.91, 383.91, 383.91, 383.91, 383.91, 254.49, 254.49, 254.49, 254.49],
        [383.91, 383.91, 383.91, 383.91, 383.91, 254.49, 254.49, 254.49, 254.49, 254.49]
        ]
    )

w = np.array([47.93, 66.37])
h = np.array([12.10, 8.6])
t = np.array([104, 48])


## Reshape the matrices to make calulations easier
h_matrix_reshaped = h_matrix.reshape(30,1)
z_matrix_reshaped = z_matrix.reshape(30,1)


# ## Initialize
# Initialize the model

m = GEKKO(remote=False)

## Fixed variables
h_constants = [m.Param(value=h_matrix_reshaped[i][0]) for i in range(30)]
z_constants = [m.Param(value=z_matrix_reshaped[i][0]) for i in range(30)]

w_base = [m.Param(value=w[i]) for i in range(w.shape[0])]
h_base = [m.Param(value=h[i]) for i in range(h.shape[0])]
t_base = [m.Param(value=t[i]) for i in range(t.shape[0])]

h_cm = m.Param(value=220)
ho_cm = m.Param(value=14.4)
w_kg = m.Param(value=21.2)
s_constraint = m.Param(value=40)

# ### Set up x var (main integer variable)
## Initialize x array

x = np.empty((30, t.shape[0]), dtype=object)
for i in range(30):
    for j in range(t.shape[0]):
        x[i][j] = m.Var(value=0, lb=0, integer=True)
        

# ### Set up Constraints

## Total constraint
for j in range(len(t_base)):
    t_contraints = sum(x[i][j] for i in range(30))
    m.Equation(t_contraints == t_base[j])


## Weight contraints
for i in range(30):
    w_constraints = sum(x[i][j]*w_base[j] for j in range(len(w_base))) + w_kg
    m.Equation(w_constraints <= z_constants[i])

## Height constraints
for i in range(30):
    h_constraints = sum(x[i][j]*h_base[j] for j in range(len(h_base))) + ho_cm
    m.Equation(h_constraints <= h_cm)


## Neighbor constraints
for i in range(9):
    # set neighbor constraints horizontally over first row
    neighbor_constraints_1 = m.abs3((sum(x[i][j]*h_base[j] for j in range(len(h_base))) +
                               h_constants[i]) - (sum(x[i+1][j]*h_base[j] for j in range(len(h_base))) +
                               h_constants[i+1]))
    m.Equation(neighbor_constraints_1 <= s_constraint)
    
for i in range(10,19):  
    # set neighbor constraints horizontally over second row
    neighbor_constraints_2 = m.abs3((sum(x[i][j]*h_base[j] for j in range(len(h_base))) +
                               h_constants[i]) - (sum(x[i+1][j]*h_base[j] for j in range(len(h_base))) +
                               h_constants[i+1]))
    m.Equation(neighbor_constraints_2 <= s_constraint)
    
for i in range(20,29):  
    # set neighbor constraints horizontally over second row
    neighbor_constraints_3 = m.abs3((sum(x[i][j]*h_base[j] for j in range(len(h_base))) +
                               h_constants[i]) - (sum(x[i+1][j]*h_base[j] for j in range(len(h_base))) +
                               h_constants[i+1]))
    m.Equation(neighbor_constraints_3 <= s_constraint)
    
for i in range(10):
    # set neighbor constrainst vertically A with B 
    neighbor_constraints_4 = m.abs3((sum(x[i][j]*h_base[j] for j in range(len(h_base))) +
                                h_constants[i]) - (sum(x[i+10][j]*h_base[j] for j in range(len(h_base))) +
                                h_constants[i+10]))
    m.Equation(neighbor_constraints_4 <= s_constraint)

for i in range(10,20):
    # set neighbor constrainst vertically B with C
    neighbor_constraints_5 = m.abs3((sum(x[i][j]*h_base[j] for j in range(len(h_base))) +
                                h_constants[i]) - (sum(x[i+10][j]*h_base[j] for j in range(len(h_base))) +
                                h_constants[i+10]))
    m.Equation(neighbor_constraints_5 <= s_constraint)
    



# ### Mix Score section below ##################

## Create a binary variable/array b that identifies if x[i][j] is non-zero or not
## We will use the count of these b[i][j] values in our objective function

## Constraint to set b[i][k] = 1 if x[i][k] > 0, and 0 otherwise
## Use if3 to set b directly based on x values
b = np.empty((30, len(t_base)), dtype=object)
epsilon = 1e-2  # Small margin allows floating-point considerations

for i in range(30):
    for j in range(len(t_base)):
        b[i][j] = m.if3(x[i][j] - epsilon, 0, 1)

# # Calculation of count(i) for each row
counts = [m.Intermediate(m.sum(b[i])) for i in range(30)]

# # x_sums for sum of each row in x
x_sums = [m.Intermediate(m.sum(x[i])) for i in range(30)]

### Mix Score section above #############################

# ## Run Solver ##

## Set a dummy objective just to identify solutions that are feasible
#m.Obj(0)

# Define the main objective function
mix_score = [counts[i] / m.max3(x_sums[i], 1e-3) for i in range(30)]
m.Obj(m.sum(mix_score))

# Set the solver options
m.options.SOLVER = 1  # APOPT solver for non-linear programs

# Increase max iterations because we don't care much about time
m.solver_options = ['minlp_gap_tol 1.0e-4',\
                    'minlp_maximum_iterations 50000',\
                    'minlp_max_iter_with_int_sol 40000']

## Solve
m.solve(disp=True)

有 APOPT 解算器调整选项，例如 

minlp_branch_method 1
和 
minlp_gap_tol 1.0e-2
，如果您同意 MacBook 上 9.56（18.5 秒）与 7.33（110.6 秒）的目标，这些选项有可能加快求解速度Pro M1 处理器。

Iter:   615 I:  0 Tm:      0.01 NLPi:    4 Dpth:  168 Lvs:  147 Obj:  9.53E+00 Gap:       NaN
Iter:   616 I: -1 Tm:      0.01 NLPi:    3 Dpth:  169 Lvs:  146 Obj:  9.53E+00 Gap:       NaN
Iter:   617 I: -1 Tm:      0.01 NLPi:    2 Dpth:  169 Lvs:  145 Obj:  9.53E+00 Gap:       NaN
Iter:   618 I:  0 Tm:      0.02 NLPi:    4 Dpth:  169 Lvs:  147 Obj:  9.55E+00 Gap:       NaN
Iter:   619 I: -1 Tm:      0.01 NLPi:    3 Dpth:  170 Lvs:  146 Obj:  9.55E+00 Gap:       NaN
--Integer Solution:   9.56E+00 Lowest Leaf:   8.05E+00 Gap:   1.72E-01
Iter:   620 I:  0 Tm:      0.01 NLPi:    4 Dpth:  170 Lvs:  145 Obj:  9.56E+00 Gap:  1.72E-01
 Warning: best integer solution returned after maximum MINLP iterations
 Adjust minlp_max_iter_with_int_sol          100  in apopt.opt to change limit
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :    18.4835999999777      sec
 Objective      :    9.56190468809906     
 Successful solution
 ---------------------------------------------------
最新问题
© www.soinside.com 2019 - 2024. All rights reserved.