在Gekko中添加条件变量后求解器无可行解问题求助

我用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)

在Mix Score部分，我定义了二进制变量b[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))

明明存在可行解，但启用这个目标函数后求解器找不到解，求解器至少应该返回那个可行解对应的目标值，求解决方法。

问题分析与解决方案 #

1. 替换非光滑的if3为线性化逻辑约束 #

m.if3会生成非光滑的分段函数，MINLP求解器APOPT处理这类结构时容易出现收敛困难。改用大M法定义二进制变量b[i][j]，把逻辑约束转化为线性约束，更适合求解器：

# 替换原有的b变量定义代码
b = np.empty((30, len(t_base)), dtype=object)
epsilon = 1e-2
M = 1e6  # 足够大的常数，大于x[i][j]可能的最大值

for i in range(30):
    for j in range(len(t_base)):
        b[i][j] = m.Var(value=0, lb=0, ub=1, integer=True)
        # 约束1：如果x[i][j] > 0，b必须为1
        m.Equation(x[i][j] <= M * b[i][j])
        # 约束2：如果b=1，x[i][j]至少大于epsilon（避免浮点误差）
        m.Equation(x[i][j] >= epsilon * b[i][j])

2. 调整目标函数的数值稳定性 #

原目标中的m.max2(x_sums[i], 1e-3)里的1e-3过小，可能导致数值奇异。换成更大的小常数，比如1e-1：

mix_score = [counts[i] / m.max2(x_sums[i], 1e-1) for i in range(30)]
m.Obj(m.sum(mix_score))

如果你的目标是最大化多样性，也可以考虑调整优化方向（比如最大化最小的counts[i]/x_sums[i]），避免求解器过度关注某些项。

3. 优化求解器参数 #

原参数可能不足以应对带目标的复杂问题，调整如下：

m.solver_options = ['minlp_gap_tol 1.0e-3',  # 放宽最优性间隙，优先找可行解
                    'minlp_maximum_iterations 60000',
                    'minlp_max_iter_with_int_sol 50000',
                    'minlp_integer_tol 1e-2']  # 放宽整数变量的容忍度

4. 用可行解作为初始点 #

先运行m.Obj(0)得到可行解，再将该解作为初始值带入带目标的模型，引导求解器从可行区域开始搜索：

# 第一步：求解可行解
m.Obj(0)
m.solve(disp=False)
# 保存x的初始值
x_initial = [[x[i][j].value[0] for j in range(len(t_base))] for i in range(30)]

# 第二步：启用目标函数
m.Obj(0)  # 先清空原目标
mix_score = [counts[i] / m.max2(x_sums[i], 1e-1) for i in range(30)]
m.Obj(m.sum(mix_score))

# 重置x的初始值为可行解
for i in range(30):
    for j in range(len(t_base)):
        x[i][j].value = x_initial[i][j]

# 再次求解
m.solve(disp=True)

内容的提问来源于stack exchange，提问作者jim