1

I'm trying to maximise an objective subject to an inequality constraint using Mystic but am struggling to see how to apply the penalty constraints. The problem is non-convex and involves maximising the objective where there is only one variable that will be changing (x). I'm trying Mystic because I've heard that it is good for large scale optimisation, and x is a 1D array contain millions of items (size N).

There are three 1-D arrays of numbers a, b, and c all with N number of values, (the values in a and b are between 0-1). Each item in x will be greater >= 0

def objective(x, a, b, c):   # maximise
    d = 1 / (1 + np.exp(-a * (x)))
    return np.sum(d * (b * c - x))

def constraint(x, f=1000):   # must be >=0
    return f - x.sum()

bounds = [(0, None)]

I've seen examples where the generate_penalty and generate_constraint functions are used and thought I could implement the constraint using the following, but it was unsuccessful:

equations = """
1000 - np.sum(x) >= 0
"""

Generally any advice on how to apply the penalty constraints or advice with using Mystic generally would be appreciated. There are lots of examples on Github but it's difficult to see which of them would be appropriate to draw from. I've implemented a solution with Scipy minimise using SLSQP but it is too slow at the scale required.

Capeboom
  • 425
  • 1
  • 5
  • 13

1 Answers1

1

I think the problem you are asking looks something like the following... although there's not really too many inequality constraints in it -- just one. I'm not using millions of items... as that would take a lot of time, and probably a good bit of parameter tuning... but I used N=100 below.

import numpy as np
import mystic as my
N = 100 #1000 # N=len(x)
M = 1e10 # max of c_i
K = 1000 # max of sum(x)
Q = 4 # 40 # npop = N*Q
G = 200 # gtol

# arrays of fixed values
a = np.random.rand(N)
b = np.random.rand(N)
c = np.random.rand(N) * M

# build objective
def cost_factory(a, b, c, max=False):
    i = -1 if max else 1
    def cost(x):
        d = 1. / (1 + np.exp(-a * x))
        return i * np.sum(d * (b * c - x))
    return cost

objective = cost_factory(a, b, c, max=True)
bounds = [(0., K)] * N

def penalty_norm(x): # < 0
    return np.sum(x) - K

# build penalty: sum(x) <= K
@my.penalty.linear_inequality(penalty_norm, k=1e12)
def penalty(x):
    return 0.0

# uncomment if want hard constraint of sum(x) == K
#@my.constraints.normalized(mass=1000)
def constraints(x):
    return x

And then in to run the script...

if __name__ == '__main__':
    mon = my.monitors.VerboseMonitor(10)
    #from pathos.pools import ThreadPool as Pool
    #from pathos.pools import ProcessPool as Pool
    #p = Pool()
    #Powell = my.solvers.PowellDirectionalSolver

    # use class-based solver interface
    """
    solver = my.solvers.DifferentialEvolutionSolver2(len(bounds), N*Q)
    solver.SetGenerationMonitor(mon)
    solver.SetPenalty(penalty)
    solver.SetConstraints(constraints)
    solver.SetStrictRanges(*my.tools.unpair(bounds))
    solver.SetRandomInitialPoints(*my.tools.unpair(bounds))
    solver.SetTermination(my.termination.ChangeOverGeneration(1e-8,G))
    solver.Solve(objective, CrossProbability=.9, ScalingFactor=.8)
    result = [solver.bestSolution]
    print('cost: %s' % solver.bestEnergy)
    """

    # use one-line interface
    result = my.solvers.diffev2(objective, x0=bounds, bounds=bounds, penalty=penalty, constraints=constraints, npop=N*Q, ftol=1e-8, gtol=G, disp=True, full_output=True, cross=.9, scale=.8, itermon=mon)#, map=p.map)

    # use ensemble of fast local solvers
    #result = my.solvers.lattice(objective, len(bounds), N*Q, bounds=bounds, penalty=penalty, constraints=constraints, ftol=1e-8, gtol=G, disp=True, full_output=True, itermon=mon)#, map=p.map)#, solver=Powell)

    #p.close(); p.join(); p.clear()
    print(np.sum(result[0]))

I've also commented out some use of parallel computing, but it's easy to uncomment.

I think you might have to work pretty hard tuning the solver to get it to find the global maximum for this particular problem. It also needs to have enough parallel elements... due to the large size of N.

However, if you wanted to use symbolic constraints as an input, you'd do it like this:

eqn = ' + '.join("x{i}".format(i=i) for i in range(N)) + ' <= {K}'.format(K=K)
constraint = my.symbolic.generate_constraint(my.symbolic.generate_solvers(my.symbolic.simplify(eqn)))

or, for the soft constraint (i.e. penalty):

penalty = my.symbolic.generate_penalty(my.symbolic.generate_conditions(eqn))
Mike McKerns
  • 33,715
  • 8
  • 119
  • 139