如何在mystic中为数组方程正确生成约束条件

1. 当前使用的符号约束解析逻辑存在变量类型适配错误：mystic的symbolic模块默认将传入的p识别为单个标量变量，而非长度为100的数组，你写的sum(p)没有被正确解析为100个数组元素的求和，生成的约束函数cf本身是无效的，因此优化过程完全没有遵循指定的约束。
2. 目标函数中缺少数值保护逻辑：当p的元素为0时，log(p)会返回nan，会干扰优化迭代的正确性。

推荐使用更不易出错的二次惩罚项来处理等式约束，无需处理复杂的符号数组解析逻辑，完整实现代码如下：

from pylab import *
from mystic.solvers import diffev2
from mystic.penalty import quadratic_equality

x = linspace(0, 5, 100)
dx = x[1]-x[0]

# 目标函数加极小epsilon避免log(0)报错
def objective(p):
    return -sum(p * log(p + 1e-12)) * dx

bounds = [(0, 1)] * len(x)

# 定义两个等式约束
def constraint1(p):
    return sum(p) * dx - 1

def constraint2(p):
    return sum(p * x) * dx - 0.2

# 生成二次惩罚项，权重设置足够大保证约束满足
@quadratic_equality(constraint1, k=1e6)
@quadratic_equality(constraint2, k=1e6)
def penalty(p):
    return 0.0

# 调大种群数量和迭代次数保证收敛
res = diffev2(
    objective, 
    x0=[1/len(x)]*len(x), 
    bounds=bounds, 
    penalty=penalty,
    npop=200,
    maxiter=1000,
    gtol=100
)

print("sum(p)*dx =", sum(res)*dx)
print("sum(p*x)*dx =", sum(res*x)*dx)

调整后输出的两个约束值会非常接近1和0.2，误差在1e-3量级以内，如果需要更高的约束精度，可以进一步调大惩罚项的k值，或者调整算法的收敛阈值参数。

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