使用mystic求解简单约束优化问题出现次优解且违反约束

我正在处理一个约束优化问题，求解结果和最优解偏差很大。为方便分析，我简化了目标函数，经解析可知：inv2的系数（1.4）高于inv1（1.1），且日期系数为正，所以应该在月末日期优先最大化inv2的资源分配。但求解器给出的结果完全不符合这个逻辑，还出现了总资源未耗尽（违反100000总额约束）的情况。我怀疑代码里有关键小错误，想先排查这些问题，再去研究差分进化算法本身。

问题代码 #

import pandas
from datetime import timedelta
from datetime import datetime
from dateutil.tz import gettz
import mystic
from mystic.symbolic import simplify
from mystic.symbolic import generate_constraint, generate_solvers
from mystic.monitors import VerboseMonitor
from mystic.symbolic import symbolic_bounds


horizon = 30
nowianatimezone = datetime.now(tz=gettz(name='Europe/Stockholm'))
datesbetweendates = pandas.date_range(start=nowianatimezone,
                                      end=nowianatimezone + timedelta(days=horizon),
                                      freq='d')
daysforcalculation = datesbetweendates.day.tolist()

def objectiveforoneday(inv1, inv2, day):
    return 1.1 * inv1 + 1.4 * inv2 + 10 * day
            
def objective(x):
    tbret = 0
    for idx, day in enumerate(daysforcalculation):
        inv1 = x[idx * 2]
        inv2 = x[(idx * 2) + 1]
        tbret += objectiveforoneday(inv1=inv1,
                                    inv2=inv2,
                                    day=day)
    return tbret

constraints = '''
x0 >= 250.0
x1 >= 100.0
x2 >= 250.0
x3 >= 100.0
x4 >= 250.0
x5 >= 100.0
x6 >= 250.0
x7 >= 100.0
x8 >= 250.0
x9 >= 100.0
x10 >= 250.0
x11 >= 100.0
x12 >= 250.0
x13 >= 100.0
x14 >= 250.0
x15 >= 100.0
x16 >= 250.0
x17 >= 100.0
x18 >= 250.0
x19 >= 100.0
x20 >= 250.0
x21 >= 100.0
x22 >= 250.0
x23 >= 100.0
x24 >= 250.0
x25 >= 100.0
x26 >= 250.0
x27 >= 100.0
x28 >= 250.0
x29 >= 100.0
x30 >= 250.0
x31 >= 100.0
x32 >= 250.0
x33 >= 100.0
x34 >= 250.0
x35 >= 100.0
x36 >= 250.0
x37 >= 100.0
x38 >= 250.0
x39 >= 100.0
x40 >= 250.0
x41 >= 100.0
x42 >= 250.0
x43 >= 100.0
x44 >= 250.0
x45 >= 100.0
x46 >= 250.0
x47 >= 100.0
x48 >= 250.0
x49 >= 100.0
x50 >= 250.0
x51 >= 100.0
x52 >= 250.0
x53 >= 100.0
x54 >= 250.0
x55 >= 100.0
x56 >= 250.0
x57 >= 100.0
x58 >= 250.0
x59 >= 100.0
x60 >= 250.0
x61 >= 100.0
x0 <= 16666.666666666668
x1 <= 16666.666666666668
x2 <= 16666.666666666668
x3 <= 16666.666666666668
x4 <= 16666.666666666668
x5 <= 16666.666666666668
x6 <= 16666.666666666668
x7 <= 16666.666666666668
x8 <= 16666.666666666668
x9 <= 16666.666666666668
x10 <= 16666.666666666668
x11 <= 16666.666666666668
x12 <= 16666.666666666668
x13 <= 16666.666666666668
x14 <= 16666.666666666668
x15 <= 16666.666666666668
x16 <= 16666.666666666668
x17 <= 16666.666666666668
x18 <= 16666.666666666668
x19 <= 16666.666666666668
x20 <= 16666.666666666668
x21 <= 16666.666666666668
x22 <= 16666.666666666668
x23 <= 16666.666666666668
x24 <= 16666.666666666668
x25 <= 16666.666666666668
x26 <= 16666.666666666668
x27 <= 16666.666666666668
x28 <= 16666.666666666668
x29 <= 16666.666666666668
x30 <= 16666.666666666668
x31 <= 16666.666666666668
x32 <= 16666.666666666668
x33 <= 16666.666666666668
x34 <= 16666.666666666668
x35 <= 16666.666666666668
x36 <= 16666.666666666668
x37 <= 16666.666666666668
x38 <= 16666.666666666668
x39 <= 16666.666666666668
x40 <= 16666.666666666668
x41 <= 16666.666666666668
x42 <= 16666.666666666668
x43 <= 16666.666666666668
x44 <= 16666.666666666668
x45 <= 16666.666666666668
x46 <= 16666.666666666668
x47 <= 16666.666666666668
x48 <= 16666.666666666668
x49 <= 16666.666666666668
x50 <= 16666.666666666668
x51 <= 16666.666666666668
x52 <= 16666.666666666668
x53 <= 16666.666666666668
x54 <= 16666.666666666668
x55 <= 16666.666666666668
x56 <= 16666.666666666668
x57 <= 16666.666666666668
x58 <= 16666.666666666668
x59 <= 16666.666666666668
x60 <= 16666.666666666668
x61 <= 16666.666666666668
x0 + x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11 + x12 + x13 + x14 + x15 + x16 + x17 + x18 + x19 + x20 + x21 + x22 + x23 + x24 + x25 + x26 + x27 + x28 + x29 + x30 + x31 + x32 + x33 + x34 + x35 + x36 + x37 + x38 + x39 + x40 + x41 + x42 + x43 + x44 + x45 + x46 + x47 + x48 + x49 + x50 + x51 + x52 + x53 + x54 + x55 + x56 + x57 + x58 + x59 + x60 + x61 == 100000'''

constraintsimplified = simplify(constraints,
                                all=True)

generatedconstraints = generate_constraint(generate_solvers(constraintsimplified))

initvals = [0.0 for _ in range(len(daysforcalculation) * 2)]
stepmon=VerboseMonitor(10)
import numpy as np
cost = lambda x: -objective(x)
optimum = mystic.solvers.diffev2(cost,
                                 x0=initvals,
                                 constraints=generatedconstraints,
                                 # bounds=bounds,
                                 #full_output=True,
                                 #stepmon=stepmon,
                                 #disp=True,
                                 npop=100,
                                 maxiter=10000)

optimalvalueslst = optimum[0]

# 求解结果
array([16666.66666667,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ])

排查出的核心问题 #

• 约束执行强度不足：
  mystic的generate_constraint默认生成软约束，不会严格强制执行等式约束（总额=100000），导致求解器可以轻易违反该约束。需要改用硬约束模式，或者在生成约束时指定handler="hard"，确保约束被严格满足。

• 单个变量上限设置错误：
  约束中给每个inv1/inv2都设置了16666.666的上限，这个数值完全不符合逻辑——总资源才100000，却允许单个变量占用近1/6的总资源，且62个变量的上限总和远超100000，这会干扰求解器的资源分配逻辑，应该删除这些无意义的单个变量上限。

• 初始值设置不合理：
  初始值全设为0，与约束要求的inv1>=250、inv2>=100冲突，差分进化算法需要从满足约束的初始种群开始迭代，否则可能陷入局部最优，无法找到正确的分配方案。

• 求解器参数未配置调试信息：
  注释掉了full_output=True和disp=True，无法查看求解过程中的约束满足情况、迭代收敛状态，难以定位问题。打开这些参数可以帮助判断是约束未生效还是算法未收敛。

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