如何为scipy的dual_annealing设置简单线性约束？

在scipy中实现约束的方法 #

scipy的dual_annealing本身不直接支持线性不等式约束，最稳妥的方式是通过变量替换将约束转化为变量的边界限制，从而直接适配原函数的求解逻辑。

针对你的约束xᵢ ≥ xᵢ₋₁ + 0.5（i≥1），我们引入新变量y替换原变量x：

• 令y₀ = x₀
• 对i≥1，令yᵢ = xᵢ - xᵢ₋₁ - 0.5，显然yᵢ ≥ 0

通过累加可将x用y表示：

• x₀ = y₀
• x₁ = y₀ + y₁ + 0.5
• x₂ = y₀ + y₁ + y₂ + 1.0
• ...
• xₙ₋₁ = y₀ + y₁ + ... + yₙ₋₁ + 0.5*(n-1)

随后将原x的边界转化为y的边界：

• y₀ ≥ 0（对应x₀ ≥ 0）
• 所有yᵢ的总和需满足y₀+...+yₙ₋₁ ≤ upper_bound - 0.5*(n-1)（对应xₙ₋₁ ≤ upper_bound），因此每个yᵢ的上界可设为该总和的最大值（若该值≤0，说明约束与上界冲突，问题无解）

最后将原目标函数改写为接受y的函数，先计算对应x再代入原函数。示例代码如下：

import numpy as np
from scipy.optimize import dual_annealing

def original_fun(x):
    # 替换为你的实际目标函数
    return np.sum(x**2)

def transformed_fun(y):
    num_points = len(y)
    x = np.zeros(num_points)
    x[0] = y[0]
    for i in range(1, num_points):
        x[i] = x[i-1] + y[i] + 0.5
    return original_fun(x)

upper_bound = 20
num_points = 30
max_sum_y = upper_bound - 0.5 * (num_points - 1)

if max_sum_y <= 0:
    raise ValueError("约束与上界冲突，无解")

# 设置y的边界
y_bounds = [(0, max_sum_y) for _ in range(num_points)]
res = dual_annealing(transformed_fun, y_bounds, maxiter=1000)

# 转换回原变量x
optimal_y = res.x
optimal_x = np.zeros(num_points)
optimal_x[0] = optimal_y[0]
for i in range(1, num_points):
    optimal_x[i] = optimal_x[i-1] + optimal_y[i] + 0.5

print("最优x值：", optimal_x)
print("最优目标函数值：", res.fun)

替代全局优化库推荐 #

若scipy的方案无法满足需求，以下库支持带约束的全局优化：

1. Pyomo #

Pyomo是建模优化问题的框架，可直接声明约束，支持Couenne、Bonmin等全局求解器：

from pyomo.environ import ConcreteModel, Var, Objective, Constraint, SolverFactory

model = ConcreteModel()
model.n = num_points
model.x = Var(range(model.n), bounds=(0, upper_bound))
model.obj = Objective(expr=original_fun(model.x))
# 添加递增约束
model.increasing_constraint = Constraint(range(1, model.n), rule=lambda m, i: m.x[i] >= m.x[i-1] + 0.5)

# 使用Couenne求解器（需单独安装）
solver = SolverFactory('couenne')
result = solver.solve(model)
optimal_x = [model.x[i]() for i in range(model.n)]

2. Optuna #

Optuna是超参数调优框架，可通过惩罚项处理约束：不满足约束时返回极大值，引导优化器避开无效解：

import optuna

def objective(trial):
    x = [trial.suggest_float(f"x{i}", 0, upper_bound) for i in range(num_points)]
    # 检查约束
    for i in range(1, num_points):
        if x[i] < x[i-1] + 0.5:
            return float('inf')
    return original_fun(x)

study = optuna.create_study(direction="minimize")
study.optimize(objective, n_trials=1000)
optimal_x = [study.best_params[f"x{i}"] for i in range(num_points)]

3. Nevergrad #

Facebook开源的全局优化库，支持直接传入约束函数：

import nevergrad as ng

def constraint(x):
    for i in range(1, len(x)):
        if x[i] < x[i-1] + 0.5:
            return False
    return True

param = ng.p.Array(shape=(num_points,)).set_bounds(lower=0, upper=upper_bound)
optimizer = ng.optimizers.OnePlusOne(parametrization=param, budget=1000)
optimizer.parametrization.constraint = constraint

result = optimizer.minimize(original_fun)
optimal_x = result.value

4. NLopt #

专门的优化库，支持多种全局算法，可直接设置线性不等式约束：

import nlopt

def objective(x, grad):
    if grad.size > 0:
        grad[:] = 2*x  # 替换为原函数的梯度（若可导）
    return original_fun(x)

opt = nlopt.opt(nlopt.GN_DIRECT_L_RAND, num_points)
opt.set_lower_bounds([0]*num_points)
opt.set_upper_bounds([upper_bound]*num_points)
opt.set_min_objective(objective)

# 添加约束：x_i - x_{i-1} >= 0.5
for i in range(1, num_points):
    coeffs = np.zeros(num_points)
    coeffs[i-1] = 1
    coeffs[i] = -1
    opt.add_inequality_constraint(lambda x, grad, coeffs=coeffs: coeffs @ x - 0.5, 1e-8)

opt.set_maxeval(1000)
optimal_x = opt.optimize(np.zeros(num_points))

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