Mystic非线性不等式约束惩罚函数域外求值问题问询

问题背景 #

我希望使用mystic求解器求解以下带非线性约束的非线性优化问题，代码如下：

import numpy as np
import matplotlib.pyplot as plt
from math import sqrt
from mystic.solvers import diffev2, fmin, fmin_powell
from mystic.monitors import VerboseMonitor
from mystic.penalty import quadratic_inequality, quadratic_equality

def pos_scale(c, q):
    return 1.0 / (1 + c*sqrt(q))

def omega_scaled(w, c, q):
    return min(w, pos_scale(c, q))

def constraints(q1, q2, w1, w2, c1, c2, fx1, fx2):
    #print('{} {}'.format(q1, q2))
    v1 = omega_scaled(w1, c1, q1)*q1*fx1
    v2 = omega_scaled(w2, c2, q2)*q2*fx2
    return v1 + v2

constraints_f = lambda q1, q2: constraints(q1, q2, 0.95, 0.92, 0.06, 0.05, 10000, 1000)
constraints_v = np.vectorize(constraints_f)

def cost(q1, q2, w1, w2, c1, c2, fx1, fx2):
    v1 = (1-omega_scaled(w1, c1, q1))*q1*fx1
    v2 = (1-omega_scaled(w2, c2, q2))*q2*fx2
    return v1 + v2

cost_f = lambda q1, q2: cost(q1, q2, 0.95, 0.92, 0.06, 0.04, 10000, 1000)
cost_v = np.vectorize(cost_f)

objective = lambda x: cost_v(x[0], x[1]).item()

def penalty_value(x, target):
    return target - constraints_f(x[0], x[1])

@quadratic_inequality(penalty_value, kwds={'target': 200000.0})
def penalty(x):
  return 0.0

我采用二次惩罚项对约束建模，约束要求定义域为正象限。求解代码和输出如下：

mon = VerboseMonitor(10)
bounds=[(0, 50), (0, 300)]
result = fmin(objective, x0=[15, 150], bounds=bounds, penalty=penalty, 
              npop=10, gtol=200, disp=False, full_output=True, itermon=mon, maxiter=500)

result

输出：

Generation 0 has ChiSquare: 77606.160271
Generation 10 has ChiSquare: 62080.449073
Generation 20 has ChiSquare: 55726.285526
Generation 30 has ChiSquare: 55505.829370
Generation 40 has ChiSquare: 55478.612377
Generation 50 has ChiSquare: 55475.462051
Generation 60 has ChiSquare: 55474.597220
Generation 70 has ChiSquare: 55474.532390
Generation 80 has ChiSquare: 55474.530891
Generation 90 has ChiSquare: 55474.530773
STOP("CandidateRelativeTolerance with {'xtol': 0.0001, 'ftol': 0.0001}")
(array([21.50326424, 42.0783277 ]), 55474.53077292251, 98, 177, 0)

当我使用合理的初始值时求解器可以找到最优解，但使用其他初始值、或是切换为Powell求解器等其他求解器时，步长搜索阶段会调用有效域外的约束函数。
请问如何才能最好地强制惩罚项中的约束函数仅在我提供给求解器的边界范围内被求值？该边界校验是否应该由求解器本身完成？还是我需要在约束函数中自行处理该逻辑？
如需可视化求解结果可参考如下代码：

fig = plt.figure(figsize=(12,6))
left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
ax = fig.add_axes([left, bottom, width, height])

q1 = np.linspace(0.1, 50, 100)
q2 = np.linspace(1, 300, 100)
X, Y = np.meshgrid(q1, q2)

Z = constraints_v(X, Y)
cp1 = plt.contour(X, Y, Z, 20, colors='black', linestyles='dashed')
cp2 = plt.contour(X, Y, Z, [200000], colors='white', linestyles='solid')
plt.clabel(cp2, inline=True, fontsize=12)

Z = cost_v(X, Y)
cp3 = plt.contourf(X, Y, Z, 25)
plt.colorbar(cp3)

sol = list(result[0])
plt.plot(sol[0], sol[1], 'go--', linewidth=2, markersize=14)

#

解决方案 #

首先明确边界校验的职责：mystic的差分进化类群体求解器原生会限制候选点在边界范围内，但Powell等局部搜索求解器默认不会强制截断步长，出现域外求值是预期行为，你可以通过以下三种方案解决，优先选择改动最小的方案一：

方案一：开启求解器严格边界约束（推荐） #

所有mystic求解器都支持strict_bounds=True参数，开启后求解器会自动将步长搜索产生的所有候选点截断到你指定的边界范围内，不会产生域外点，无需修改其他业务代码：

# 以Powell求解器为例，fmin等其他求解器参数通用
result = fmin_powell(objective, x0=[15, 150], bounds=bounds, penalty=penalty, 
              gtol=200, disp=False, full_output=True, itermon=mon, maxiter=500,
              strict_bounds=True)

方案二：函数内增加边界校验兜底 #

如果需要兼容旧版本mystic，或者有自定义边界逻辑，可以在目标函数和惩罚项函数入口增加边界判断，域外点直接返回极高的惩罚值，避免后续逻辑报错：

# 改造目标函数
def objective(x):
    q1, q2 = x
    # 超出边界直接返回无穷大，求解器会自动抛弃这类候选点
    if q1 < 0 or q1 >50 or q2 <0 or q2>300:
        return 1e18
    return cost_v(q1, q2).item()

# 改造惩罚项计算函数
def penalty_value(x, target):
    q1, q2 = x
    if q1 < 0 or q1 >50 or q2 <0 or q2>300:
        return 1e18
    return target - constraints_f(x[0], x[1])

方案三：使用原生约束工具包装边界 #

你也可以用mystic自带的约束工具生成边界约束，和现有惩罚项配合传入求解器：

from mystic.constraints import impose_bounds
# 生成边界约束
bounds_constraint = impose_bounds(lb=(0,0), ub=(50,300))
# 求解时传入constraints参数
result = fmin(objective, x0=[15, 150], bounds=bounds, penalty=penalty,
              constraints=bounds_constraint,
              npop=10, gtol=200, disp=False, full_output=True, itermon=mon, maxiter=500)

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