Search code examples
pythonoptimizationminimizationscipy-optimize

Choosing variables for scipy.optimize from a pre-defined set


I am trying to minimize a function with scipy.optimize with three input variables, two of which are bounded and one has to be chosen from a set of values. To ensure that the third variable is chosen from a predefined set of values, I introduced the following constraint:

from scipy.optimize import rosin, shgo
import numpy as np

# Set from which the third variable to be optimized can hold 
Z = np.array([-1, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1]) 

def Reson_Test(x):   # arbitrary objective function
     print (x)
     return rosen(x)**2 - np.sin(x[0]) 

def Cond_1(x):
     if x[2] in Z:
          return 1
     else:
          return -1

bounds = [(-512,512),]*3 
conds = ({'type': 'ineq' , 'fun' : Cond_1})


result = shgo(Rosen_Test, bounds, constraints=conds)
print (result)

However, when looking at the print results from Rosen_Test, it is evident that the condition is not being enforced - perhaps condition is not defined correctly?

I was wondering if anyone has any ideas to ensure that the third variable can be chosen from a set.

Note: The use of the shgo method was chosen such that constraints can be introduced and can be changed. Also, I am open to use other optimization packages if this condition is met


Solution

  • The inequality constraints do not work like that.

    As mentioned in the docs they are defined as

    g(x) <= 0
    

    and you need to write g(x) work like that. In your cases that is not the case. You are only returning a single scalar for one dimension. You need to return a vector with three dimensions, of shape (3,).

    In your case you could try to use the equality constraints instead, as this could allow a slightly better hack. But I am still not sure if it will work as these optimizers don't work like that. And the whole thing will probably leave the optimizer with a rather bumpy and discontinuous objective function. You can read on Mixed-Integer Nonlinear Programming (MINLP), maybe start here.

    There are is one more reasons why your approach won't work as expected As optimizers work with floating point numbers it will likely never find a number in your array when optimizing and guessing new solutions.

    This illustrates the issue:

    import numpy as np
    
    Z = np.array([-1, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1])
    
    print(0.7999999 in Z) # False, this is what the optimizer will find
    print(0.8 in Z) # True, this is what you want
    

    Maybe you should try to define your problem in a way that allows to use an inequality constraint on the whole range of Z.

    But let's see how it could work.

    An equality constraint is defined as

    h(x) == 0
    

    So you could use

    def Cond_1(x):
        if x[2] in Z:
            return numpy.zeros_like(x)
        else:
            return numpy.ones_like(x) * 1.0 # maybe multiply with some scalar?
    

    The idea is to return an array [0.0, 0.0, 0.0] that satisfies the equality constraint if the number is found. Else return [1.0, 1.0, 1.0] to show that it is not satisfied.

    Caveats:

    1.) You might have to tune this to return an array like [0.0, 0.0, 1.0] to show the optimizer which dimension you are unhappy about so the optimizer can make better guesses by only adjusting a single dimension.

    2.) You might have to return a larger value than 1.0 to state an non-satisfied equality constraint. This depends on the implementation. The optimizer could think that 1.0 is fine as it is close to 0.0. So maybe you have to try something [0.0, 0.0, 999.0].

    This solves the problem with the dimension. But still will not find any numbers do to the floating point number thing mentioned above.

    But we can try to hack this like

    import numpy as np
    
    Z = np.array([-1, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1])
    
    def Cond_1(x):
    
        # how close you want to get to the numbers in your array
        tolerance = 0.001 
        delta = np.abs(x[2] - Z)
        print(delta)
        print(np.min(delta) < tolerance)
    
        if np.min(delta) < tolerance:
    
            return np.zeros_like(x)
    
        else:
    
            # maybe you have to multiply this with some scalar
            # I have no clue how it is implemented
            # we need a value stating to the optimizer "NOT THIS ONE!!!"
            return np.ones_like(x) * 1.0
    
    
    sol = np.array([0.5123, 0.234, 0.2])
    print(Cond_1(sol)) # True
    
    sol = np.array([0.5123, 0.234, 0.202])
    print(Cond_1(sol)) # False
    

    Here are some recommendations on optimization. To make sure it works in a reliable way try to start the optimization at different initial values. Global optimization algorithms might not have initial values if used with boundaries. The optimizer somehow discretizes the space.

    What you could do to check the reliability of your optimization and get better overall results:

    • Optimize on the complete region [-512, 512] (for all three dimensions)

    • Try 1/2 of that: [-512, 0] and [0, 512] (8 sub-optimizations, 2 for each dimension)

    • Try 1/3 of that: [-512, -171], [-171, 170], [170, 512] (27 sub-optimizations, 3 for each dimension)

    • Now compare the converged results to see if the complete global optimization found the same result

    • If the global optimizer did not find the "real" minima but the sub-optimization:

      • your objective function is too difficult on the whole domain
      • try a different global optimizer
      • tune the parameters (maybe the 999 for the equality constraint)
      • I often use the sub-optimization as part of the normal process, not only for testing. Especially for blackbox problems.

    Please also see these answers:

    Scipy.optimize Inequality Constraint - Which side of the inequality is considered?

    Scipy minimize constrained function