pythonpython-3.xnonlinear-optimizationmystic

How to specify multiple constraints in function form in Mystic?


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?


Solution

  • 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.