如何在 mystic 中正确生成数组方程的约束？

Question

我有一个优化问题，我想在约束 sum( p)dx == 1 和 sum(p.x)dx == 0.2，其中 p 和 x 是数组，dx 是标量。我尝试以这种方式使用 mystic 来实现：

from pylab import *
from mystic.symbolic import generate_constraint, generate_solvers, simplify
from mystic.solvers import diffev2

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

def objective(p):
    return -sum(p*log(p))*dx

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

equations = '''
sum(p)*dx == 1
sum(p*x)*dx == 0.2
'''

eqns = simplify(equations, variables=["p"], locals={"dx":dx, "x":x})

cf = generate_constraint(generate_solvers(eqns))

res = diffev2(objective, x0=[1/len(x)]*len(x), bounds=bounds, constraints=cf)

print(sum(res)*dx)
print(sum(res*x)*dx)

但显然效果不佳，因为returns:

Optimization terminated successfully.
         Current function value: 0.145182
         Iterations: 264
         Function evaluations: 26500
0.030782323656083955
0.07399192217757057

因此违反了约束。我应该如何正确解决我的问题？

Answer 1

我是 mystic 作者。主要问题是您假设 mystic.symbolic 可以理解矢量符号，但事实并非如此。

>>> from mystic.symbolic import generate_constraint, generate_solvers, simplify
>>> from mystic.solvers import diffev2
>>> from numpy import *
>>> x = linspace(0, 5, 100)
>>> dx = x[1]-x[0]
>>> def objective(p):
...     return -sum(p*log(p))*dx
... 
>>> bounds = [(0,1)]*len(x)
>>> equations = '''
... sum(p)*dx == 1
... sum(p*x)*dx == 0.2
... '''
>>> eqns = simplify(equations, variables=["p"], locals={"dx":dx, "x":x})
>>> print(eqns)
p == 0.0158400000000000
p == 19.8000000000000
>>>

这显然不是你想要的。所以，如果你想使用符号约束，那么你必须逐个元素地写出来。

>>> equations = 'dx * (' + ' + '.join('p%s' % i for i,j in enumerate(x)) + ') ==
 1\n'
>>> equations += 'dx * (' + ' + '.join('p%s * x%s' % (i,i) for i,j in enumerate(x)) + ') == .2\n'
>>> locals = {'x%s' % i:j for i,j in enumerate(x)}
>>> locals['dx'] = dx
>>> eqns = simplify(equations, variables=['p%s' % i for i,j in enumerate(x)], locals=locals)
>>> eqns.split('\n',1)[0]
'p10 == -0.1*p1 - 1.1*p11 - 1.2*p12 - 1.3*p13 - 1.4*p14 - 1.5*p15 - 1.6*p16 - 1.7*p17 - 1.8*p18 - 1.9*p19 - 0.2*p2 - 2.0*p20 - 2.1*p21 - 2.2*p22 - 2.3*p23 - 2.4*p24 - 2.5*p25 - 2.6*p26 - 2.7*p27 - 2.8*p28 - 2.9*p29 - 0.3*p3 - 3.0*p30 - 3.1*p31 - 3.2*p32 - 3.3*p33 - 3.4*p34 - 3.5*p35 - 3.6*p36 - 3.7*p37 - 3.8*p38 - 3.9*p39 - 0.4*p4 - 4.0*p40 - 4.1*p41 - 4.2*p42 - 4.3*p43 - 4.4*p44 - 4.5*p45 - 4.6*p46 - 4.7*p47 - 4.8*p48 - 4.9*p49 - 0.5*p5 - 5.0*p50 - 5.1*p51 - 5.2*p52 - 5.3*p53 - 5.4*p54 - 5.5*p55 - 5.6*p56 - 5.7*p57 - 5.8*p58 - 5.9*p59 - 0.6*p6 - 6.0*p60 - 6.1*p61 - 6.2*p62 - 6.3*p63 - 6.4*p64 - 6.5*p65 - 6.6*p66 - 6.7*p67 - 6.8*p68 - 6.9*p69 - 0.7*p7 - 7.0*p70 - 7.1*p71 - 7.2*p72 - 7.3*p73 - 7.4*p74 - 7.5*p75 - 7.6*p76 - 7.7*p77 - 7.8*p78 - 7.9*p79 - 0.8*p8 - 8.0*p80 - 8.1*p81 - 8.2*p82 - 8.3*p83 - 8.4*p84 - 8.5*p85 - 8.6*p86 - 8.7*p87 - 8.8*p88 - 8.9*p89 - 0.9*p9 - 9.0*p90 - 9.1*p91 - 9.2*p92 - 9.3*p93 - 9.4*p94 - 9.5*p95 - 9.6*p96 - 9.7*p97 - 9.8*p98 - 9.9*p99 + 7.8408'

这看起来有点像您想要的了。

>>> from mystic.constraints import and_
>>> from mystic.monitors import VerboseMonitor
>>> cf = generate_constraint(generate_solvers(eqns, variables=['p%s' % i for i,j in enumerate(x)]), join=and_)
>>> res = diffev2(objective, x0=[1/len(x)]*len(x), bounds=bounds, constraints=cf, itermon=mon, npop=400)
Generation 0 has ChiSquare: inf
Generation 1 has ChiSquare: inf
Generation 2 has ChiSquare: inf
Generation 3 has ChiSquare: inf
Generation 4 has ChiSquare: inf
Generation 5 has ChiSquare: inf
Generation 6 has ChiSquare: inf
Generation 7 has ChiSquare: inf
Generation 8 has ChiSquare: inf
Generation 9 has ChiSquare: inf
Generation 10 has ChiSquare: inf
Generation 11 has ChiSquare: inf
Generation 12 has ChiSquare: inf
Generation 13 has ChiSquare: inf
Generation 14 has ChiSquare: inf
Generation 15 has ChiSquare: inf
Generation 16 has ChiSquare: inf
Generation 17 has ChiSquare: inf
Generation 18 has ChiSquare: inf
Generation 19 has ChiSquare: inf
Generation 20 has ChiSquare: inf
Generation 21 has ChiSquare: inf
Generation 22 has ChiSquare: inf
Generation 23 has ChiSquare: inf
Generation 24 has ChiSquare: inf
Generation 25 has ChiSquare: inf
Generation 26 has ChiSquare: inf
Generation 27 has ChiSquare: inf
Generation 28 has ChiSquare: inf
Generation 29 has ChiSquare: inf
Generation 30 has ChiSquare: inf
STOP("VTRChangeOverGeneration with {'ftol': 0.005, 'gtol': 1e-06, 'generations': 30, 'target': 0.0}")
Optimization terminated successfully.
         Current function value: inf
         Iterations: 30
         Function evaluations: 12400

您可能想要比这更多的人口，而且我在结果中看到很多 inf，您可能想回到约束并添加严格的界限。后者可以通过将 tight=True 传递给求解器来完成，但它比首次定义方程时明确添加边界要慢。像这样：

>>> from mystic.symbolic import symbolic_bounds
>>> equations_ = symbolic_bounds(*zip(*bounds))

然后您需要合并和简化 equations 和 equations_ 作为约束函数的基础。

作为替代方案，您可以使用惩罚...（一旦您拥有简化的 eons，构建如下：

>>> from mystic.symbolic import generate_penalty, generate_conditions
>>> from mystic.coupler import and_ as _and
>>> pf = generate_penalty(generate_conditions(eqns, variables=['p%s' % i for i,j
 in enumerate(x)]), join=_and)
>>> mon = VerboseMonitor(1000)
>>> res = fmin(objective, x0=[1/len(x)]*len(x), bounds=bounds, penalty=pf, itermon=mon, maxiter=10000)
Generation 0 has ChiSquare: 36179.905048
Generation 1000 has ChiSquare: 34612.689290
Generation 2000 has ChiSquare: 32100.152737
Generation 3000 has ChiSquare: 30305.589205
Generation 4000 has ChiSquare: 29627.956352
Generation 5000 has ChiSquare: 29230.517427
Generation 6000 has ChiSquare: 28933.796220
Generation 7000 has ChiSquare: 28804.810192
Generation 8000 has ChiSquare: 28664.350695
Generation 9000 has ChiSquare: 28572.124583
Generation 10000 has ChiSquare: 28547.949725
Warning: Maximum number of iterations has been exceeded
STOP("EvaluationLimits with {'evaluations': 20000, 'generations': 10000}")

惩罚是软约束，看起来上面还有一些违反约束的地方，但是正在朝着好的解决方案发展。

无论如何，上面的内容应该会给你足够的选择来克服你被困的地方，并找到解决方案。

如何在 mystic 中正确生成数组方程的约束？

How to properly generate constraints for array equations in mystic?

python

optimization

mystic