2

How do you input multiple constraints in using the Mystic solver? For instance, I have these two functions which describes my two constraints:

def constraint1(x):
    return x[0]*x[1]*x[2]*x[3]-25.0

def constraint2(x):
    sum_eq = 40.0
    for i in range(4):
        sum_eq = sum_eq - x[i]**2
    return sum_eq

When using diffev(...), what is the right way to input these constraints?

Mike McKerns
  • 33,715
  • 8
  • 119
  • 139
warrented
  • 23
  • 2

1 Answers1

2

I'm the mystic author. I'm not sure what you want... You have written the constraints in the form of soft constraints (i.e. penalties), so I'll go with that.

To couple penalties, you can use a coupler like mystic.coupler.and_, or you can use a penalty decorator, like this:

>>> import mystic as my
>>>
>>> def penalty1(x):
...   return x[0]*x[1]*x[2]*x[3] - 25.0
... 
>>> def penalty2(x):
...   return 40.0 - sum(i*i for i in x)
... 
>>> x
[4, 5, 1, 9]
>>> 
>>> @my.penalty.linear_equality(penalty1)
... @my.penalty.linear_equality(penalty2)
... def penalty(x):
...   return 0.0
... 
>>> penalty(x)
23800.0
>>> 

Then, if all you want to do is solve where the penalty is smallest, you can use diffev on the penalty directly:

>>> my.solvers.diffev(penalty, [10,10,10,10], npop=100)
Optimization terminated successfully.
         Current function value: 0.001920
         Iterations: 167
         Function evaluations: 16800
array([1.92486578, 0.86032337, 2.89649062, 5.21201382])
>>> 

If you want to apply these as a penalty, while solving some other function, then you'd do it like this:

>>> from mystic.models import rosen  # the Rosenbrock test function 
>>> my.solvers.diffev(rosen, [0,0,0,0], penalty=penalty, npop=100)
Optimization terminated successfully.
         Current function value: 2.263629
         Iterations: 182
         Function evaluations: 18300
array([1.24923882, 1.53972652, 2.35201514, 5.5259994 ])

If you wanted hard constraints, that's a different matter:

>>> def constraint1(x, val=29):
...     res = my.constraints.normalize([i*i for i in x], val)
...     return [i**.5 for i in res]
... 
>>> def constraint2(x, val=24):
...     return my.constraints.impose_product(val, x)
... 
>>> constraint3 = my.constraints.integers()(lambda x:x)
>>> 
>>> constraint = my.constraints.and_(constraint1, constraint2, constraint3)
>>> 
>>> x = [1,2,3]
>>> z = constraint(x)
>>> 
>>> assert constraint1(z) == z
>>> assert constraint2(z) == z
>>> assert constraint3(z) == z
>>> z
[2.0, 3.0, 4.0]

They are applied to the solvers with the keyword constraints.

Mike McKerns
  • 33,715
  • 8
  • 119
  • 139