Minimizing multi-objective function using DEAP

Question

I'm trying to perform multi-objective optimisation by minimizing a custom function using the DEAP library. Although I'm getting decent results when minimizing for several objectives (targets), for more than 3 or 4 it fails to converge. Typically it will minimize the first objective to 0, whilst leaving the other objectives bouncing around (not minimizing).

I built a meta-model (ridge regression) to describe some simulation data using sci-kit library, so my model is based on coefficients and intercepts (included in my code). New predictions are based on ~150 inputs that are being varied uniformly.

There is a year option that minimizes for 3 targets, and a month option that minimizes for 8 targets.

I've included my code as a gist as it is quite large. Please find it HERE.

Question: Anyone has any idea what the reason might be for the remaining objectives not being minimized? I've tried playing around with the selection, mutation and crossover processes, but no luck yet. Or could it possibly be related to the model itself? I've also tried different weights for the fitness, but for some reason it doesn't seem to make a difference.

Results for the year targets:

Results for the monthly targets:

Answer 1

Just to answer my own question.

It seems that I wasn't returning the right kind of values during evaluation.

Changing to the RMSE of the differences instead of absolute differences between objectives and predictions did the trick:

def EvaluateObjective(individual):
    prediction = calculate(individual, for_sensitivity)
    prediction = [int(i) for i in prediction]

    # diff = []
    # for y in range(len(targets)):
    #     output = math.sqrt((targets[y] - prediction[y]) ** 2)
    #     #output = abs(targets[y] - prediction[y])
    #     diff.append(output)

    rmse = np.sqrt((sum((i - j)**2 for i, j in zip(prediction, targets)) / len(targets)))

    return (rmse,)

Answer 2

You gave me the solution to the exact same issue I've been struggling with. Nice method, a little trick made my program work too!

I'm pretty sure there must be lots of deap users who tried to use more than two weights like weights=(-1.0, -1.0, 1.0) as I did.

I will post the simple example of 3 parameters (2 parameters minimized, 1 parameter maximized.)

• The example is about "How to load as many items as possible with conditions of maximum weight, maximum size"

• Conditions :

  1. Minimize weight sum.
  2. Minimize size sum.
  3. Maximize a sum of values.

from numpy import array
import numpy
import random
from deap import base, creator, tools, algorithms

###  Multi-objective Optimization Problem  ###

IND_INIT_SIZE = 5

MAX_WEIGHT = 2000 # kg
MAX_SIZE = 1500 # m**3


# Create the item dictionary:

r = array([[213, 508,  22],  # 1st arg : weight / 2nd arg : size / 3rd arg : value
       [594, 354,  50],
       [275, 787,  43],
       [652, 218,  46],
       [728, 183,  43],
       [856, 308,  33],
       [727, 482,  45],
       [762, 683,  26],
       [707, 450,  19],
       [909, 309,  45],
       [979, 247,  42],
       [259, 705,  42],
       [260, 543,  14],
       [899, 825,  17],
       [446, 360,  35],
       [491, 818,  47],
       [647, 404,  17],
       [604, 623,  32],
       [900, 840,  45],
       [374, 127,  33]] )


NBR_ITEMS = r.shape[0]

items = {}
# Create random items and store them in the items' dictionary.
for i in range(NBR_ITEMS):
    items[i] = ( r[i][0] , r[i][1] , r[i][2] )


creator.create("Fitness", base.Fitness, weights=(-1.0, 1.0 ))  # Note here <- I used only two weights!  (at first, I tried weights=(-1.0 , -1.0, 1.0)) but it crashes. With deap, you cannot do such a thing.

creator.create("Individual", set, fitness=creator.Fitness)

toolbox = base.Toolbox()

# Attribute generator
toolbox.register("attr_item", random.randrange, NBR_ITEMS)

# Structure initializers
toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_item, n=IND_INIT_SIZE) #
toolbox.register("population", tools.initRepeat, list, toolbox.individual)


def evaluation(individual):
    weight = 0.0
    size =0.0
    value = 0.0

    # Maximize or Minimize Conditions
    for item in individual:
        weight += items[item][0]  # It must be minimized.
        size += items[item][1]  # It must be minimized.
        value += items[item][2]  # It must be maximized.

    # Limit Conditions
    if weight > MAX_WEIGHT or size > MAX_SIZE:
        return 10000, 0

    if value == 0:
        value = 0.0000001

    MinFitess_score = weight + size   # NOTE : Minimize weight, size
    MaxFitenss_score = value  # NOTE : Maximize weight, size

    return MinFitess_score , MaxFitenss_score,



def cxSet(ind1, ind2):
    """Apply a crossover operation on input sets. The first child is the
    intersection of the two sets, the second child is the difference of the
    two sets.
    """
    temp = set(ind1)  # Used in order to keep type
    ind1 &= ind2  # Intersection (inplace)
    ind2 ^= temp  # Symmetric Difference (inplace)
    return ind1, ind2


def mutSet(individual):
    """Mutation that pops or add an element."""
    if random.random() < 0.5:
        if len(individual) > 0:  # We cannot pop from an empty set
            individual.remove(random.choice(sorted(tuple(individual))))
    else:
        individual.add(random.randrange(NBR_ITEMS))
    return individual,  # NOTE comma(,) , if there's no comma, an error occurs.



toolbox.register("mate", cxSet)
toolbox.register("mutate", mutSet)
toolbox.register("select", tools.selNSGA2) # NSGA-2 applies to multi-objective problems such as knapsack problem
toolbox.register("evaluate", evaluation)


def main():
    ngen = 300  # a number of generation  < adjustable value >

    pop = toolbox.population(n= 300)
    hof = tools.ParetoFront() # a ParetoFront may be used to retrieve the best non dominated individuals of the evolution
    stats = tools.Statistics(lambda ind: ind.fitness.values)
    stats.register("avg", numpy.mean, axis=0)
    stats.register("std", numpy.std, axis=0)
    stats.register("min", numpy.min, axis=0)
    stats.register("max", numpy.max, axis=0)

    algorithms.eaSimple(pop, toolbox, 0.7, 0.2, ngen=ngen, stats=stats, halloffame=hof, verbose=True)

    return hof, pop


if __name__ == "__main__":
    hof, pop = main()

    print(hof) # non-dominated individuals' list  # the fittest value is placed on the most right side.

The ideal result:

1. Individual({1, 2, 19, 4}) or
2. Individual({1, 2, 19, 3})

as their total scores are quite similar. You will get one of the results.