pythonoptimizationnonlinear-optimizationmystic

constrain f(x) with g(x) when minimizing f(x) in mystic?


Hi I am new to mystic so I apologize if I ask similar questions that have been answered before.

Say I have an input x (a list of length n), and a function f(x) that will map the input to a m-dimensional output space. I also have a m-dimensional constraint list that also depends on the input x (say it is g(x)), and I want to make sure the output is less than the constraint for each element of the m-dimensional list. How should I specify this in mystic?

simplified example:

x = [1,2,3]
output = [6, 12] # with some f(x)
constraints = [5, 10] # with some g(x)

Solution

  • I'm the mystic author. The easiest way I'd constrain one function's output with another is by using a penalty (soft constraint).

    I'll show a simple case, similar to what I think you are looking for:

    """
    2-D input x
    1-D output y = f(x)
    where y > g(x)
    
    with f(x) = x0^(sin(x0)) + x1^4 + 6x1^3 - 5x1^2 - 40x1 + 35
    and g(x) = 7x1 - x0 + 5
    in x = [0,10]
    """
    

    Let's find the minimum of f, larger than g.

    >>> import numpy as np
    >>> import mystic as my
    >>> 
    >>> def f(x):
    ...     x0,x1 = x
    ...     return x0**np.sin(x0) + x1**4 + 6*x1**3 - 5*x1**2 - 40*x1 + 35
    ... 
    >>> def g(x):
    ...     x0,x1 = x
    ...     return 7*x1 - x0 + 5
    ... 
    >>> def penalty(x):
    ...     return g(x) - f(x)
    ... 
    >>> @my.penalty.quadratic_inequality(penalty, k=1e12)
    ... def p(x):
    ...     return 0.0
    ... 
    >>> mon = my.monitors.Monitor()
    >>> my.solvers.fmin(f, [5,5], bounds=[(0,10)]*2, penalty=p, itermon=mon, disp=1)
    Optimization terminated successfully.
             Current function value: 11.898373
             Iterations: 89
             Function evaluations: 175
    STOP("CandidateRelativeTolerance with {'xtol': 0.0001, 'ftol': 0.0001}")
    array([7.81765653, 2.10228969])
    >>> 
    >>> my.log_reader(mon)
    

    log_reader(mon) Make a plot to visualize the constrained solver trajectory.

    >>> fig = my.model_plotter(f, mon, depth=True, scale=1.0, bounds="0:10, 0:10", out=True)
    >>> 
    >>> from matplotlib import cm
    >>> import matplotlib.pyplot as plt
    >>> 
    >>> x,y = my.scripts._parse_axes("0:10, 0:10", grid=True)
    >>> x, y = np.meshgrid(x, y)
    >>> z = 0*x
    >>> s,t = x.shape
    >>> for i in range(s):
    ...     for j in range(t):
    ...         xx,yy = x[i,j], y[i,j]
    ...         z[i,j] = g([xx,yy])
    ... 
    >>> z = np.log(4*z*1.0+1)+2 # scale=1.0
    >>> ax = fig.axes[0]
    >>> ax.contour(x, y, z, 50, cmap=cm.cool)
    <matplotlib.contour.QuadContourSet object at 0x12bc0d0d0>
    >>> plt.show()
    

    model_plotter(f, mon, ...)

    From the above plot, you can sort of see that the solver started minimizing f (color=jet), then hit g (color=cool), and traced along the intersection until it hit the minimum.