diff --git a/.gitignore b/.gitignore
index 180218d4e..0344fed1a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -15,5 +15,6 @@ MANIFEST
 *.pyc
 *.pyo
 *.orig
+*~
 .DS_Store
 *.so
diff --git a/deap/_algorithms/__init__.py b/deap/_algorithms/__init__.py
new file mode 100644
index 000000000..3cabb7672
--- /dev/null
+++ b/deap/_algorithms/__init__.py
@@ -0,0 +1,16 @@
+#    This file is part of DEAP.
+#
+#    DEAP is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU Lesser General Public License as
+#    published by the Free Software Foundation, either version 3 of
+#    the License, or (at your option) any later version.
+#
+#    DEAP is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#    GNU Lesser General Public License for more details.
+#
+#    You should have received a copy of the GNU Lesser General Public
+#    License along with DEAP. If not, see <http://www.gnu.org/licenses/>.
+
+from .moead import MOEAD
diff --git a/deap/_algorithms/moead.py b/deap/_algorithms/moead.py
new file mode 100644
index 000000000..45d92a53a
--- /dev/null
+++ b/deap/_algorithms/moead.py
@@ -0,0 +1,353 @@
+#    This file is part of DEAP.
+#
+#    DEAP is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU Lesser General Public License as
+#    published by the Free Software Foundation, either version 3 of
+#    the License, or (at your option) any later version.
+#
+#    DEAP is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#    GNU Lesser General Public License for more details.
+#
+#    You should have received a copy of the GNU Lesser General Public
+#    License along with DEAP. If not, see <http://www.gnu.org/licenses/>.
+
+import random
+import operator
+
+import numpy as np
+
+from deap import tools
+
+def sort_index(array):
+    """Sort array index in assending order.
+
+    :params array: An array to be sorted.
+    :returns: An index array sorted in assending order for given array.
+    """
+    index = [i for i in range(len(array))]
+    for i in range(1, len(array)):
+        for j in range(i):
+            if array[index[i]] < array[index[j]]:
+                temp = index[i]
+                for k in range(i-1, j-1, -1):
+                    index[k+1] = index[k]
+                index[j] = temp
+                break
+    return index
+
+class MOEAD:
+    """Implementation of MOEA/D.
+    """
+    def __init__(self, population, toolbox, weights, neighbour_size,
+                 scalar_method, cxpb, mutpb, ngen, stats, verbose):
+        self.population = population
+        self.toolbox = toolbox
+        self.weights = np.array(weights)
+        self.neighbour_size = neighbour_size
+        self.scalar_method = scalar_method
+        self.cxpb = cxpb
+        self.mutpb = mutpb
+        self.ngen = ngen
+        self.stats = stats
+        self.verbose = verbose
+
+        if hasattr(toolbox, 'improve'):
+            self.improve = self.toolbox.improve
+        else:
+            self.improve = self.noop_improve
+
+        if self.scalar_method == 'weightedSum':
+            self.scalarize_objective = self.weighted_sum
+        elif self.scalar_method == 'tchebycheff':
+            self.scalarize_objective = self.tchebycheff
+        else:
+            raise ValueError('unknown scalar method')
+
+    def run(self):
+        """Execute MOEA/D.
+
+        :returns: The external population
+        :returns: A class:`~deap.tools.Logbook` with the statistics of the
+                  evolution
+        """
+        self.logbook = tools.Logbook()
+        logbook_header  = ['gen', 'nevals']
+        logbook_header += (self.stats.fields if self.stats else [])
+        self.logbook.header = logbook_header
+
+        # Step 1)
+        self.initialize()
+        npop = len(self.population)
+
+        record = self.stats.compile(self.population) if self.stats else {}
+        self.logbook.record(gen=0, nevals=npop, **record)
+        if self.verbose:
+            print(self.logbook.stream)
+
+        # Step 2)
+        for i in range(self.ngen):
+            for j in range(npop):
+                self.update(j)
+
+            record = self.stats.compile(self.population) if self.stats else {}
+            self.logbook.record(gen=i+1, nevals=npop, **record)
+            if self.verbose:
+                print(self.logbook.stream)
+
+        return self.EP, self.logbook
+
+    def initialize(self):
+        """Step 1) Initialization.
+        """
+
+        # Step 1.1)
+        self.EP = []
+
+        # Step 1.2)
+        self.initialize_neighbour()
+
+        # Step 1.3)
+        # Generate an initial population at this step in the paper.
+        # But this implementation assumes population as argument.
+
+        # Step 1.4)
+        self.ideal_point = []
+        for w in self.population[0].fitness.weights:
+            if w > 0:
+                # Initialize very small value for maximize
+                self.ideal_point.append(-1.0e+30)
+            else:
+                # Initialize very large value for minimize
+                self.ideal_point.append(1.0e+30)
+
+        for ind in self.population:
+            ind.fitness.values = self.toolbox.evaluate(ind)
+            self.update_reference(ind)
+
+    def initialize_neighbour(self):
+        """Step 1.2) Compute Euclidean distances between any two weight vectors
+        and then work out the T closest weight vectors to each weight vector.
+        """
+        self.neighbour_table = []
+        npop = len(self.population)
+        distance_matrix = np.zeros((npop, npop))
+        for i in range(npop):
+            for j in range(i+1, npop):
+                distance = np.linalg.norm(self.weights[i] - self.weights[j])
+                distance_matrix[j, i] = distance_matrix[i, j] = distance
+
+        for i in range(npop):
+            index = sort_index(distance_matrix[i])
+            self.neighbour_table.append(index[:self.neighbour_size])
+
+    def update(self, i):
+        """Step 2) Update.
+
+        :params i: An index of subproblem.
+        """
+        offspring = self.generate_offspring(i)
+        offspring = self.improve(i, offspring)
+        offspring.fitness.values = self.toolbox.evaluate(offspring)
+        self.update_reference(offspring)
+        self.update_neighbours(i, offspring)
+        self.update_EP(offspring)
+
+    def generate_offspring(self, i):
+        """Step 2.1) Reproduction.
+
+        :params i: An index of subproblem.
+        """
+        neighbour_table = self.neighbour_table[i][:]
+        # Select neighbour but not i
+        neighbour_table.remove(i)
+        k = random.choice(neighbour_table)
+        # Select neighbour but not i, k
+        neighbour_table.remove(k)
+        l = random.choice(neighbour_table)
+
+        c1 = self.population[k]
+        c2 = self.population[l]
+
+        choice = random.random()
+        if choice < self.cxpb:
+            c1 = self.toolbox.clone(c1)
+            c2 = self.toolbox.clone(c2)
+            c1, c2 = self.toolbox.mate(c1, c2)
+            del c1.fitness.values
+        choice = random.random()
+        if choice < self.mutpb:
+            c1 = self.toolbox.clone(c1)
+            c1, = self.toolbox.mutate(c1)
+            del c1.fitness.values
+        # Ignore one individual. Bad idea?
+        offspring = c1
+
+        return offspring
+
+    def noop_improve(self, i, offspring):
+        """Step 2.2) Improvement.
+        No improvement default. Register :meth:`toolbox.improve` for problem
+        specific improvement.
+
+        :params i: An index of subproblem.
+        :params offspring: individual to improve.
+        :returns: improved individual.
+        """
+        return offspring
+
+    def update_reference(self, ind):
+        """Step 2.3) Uodate of z.
+        Also use by Step 1.4) initialize z.
+
+        :params ind: An individual that maybe update ideal point.
+        """
+        for i in range(len(ind.fitness.values)):
+            if ind.fitness.weights[i] > 0:
+                if ind.fitness.values[i] > self.ideal_point[i]:
+                    self.ideal_point[i] = ind.fitness.values[i]
+            else:
+                if ind.fitness.values[i] < self.ideal_point[i]:
+                    self.ideal_point[i] = ind.fitness.values[i]
+
+    def update_neighbours(self, i, offspring):
+        """Step 2.4) Update of Neighbouring solutions.
+
+        :params i: An index of subproblem.
+        :params offspring: An individual that maybe update neighbour.
+        """
+        for j in range(self.neighbour_size):
+            weight_index = self.neighbour_table[i][j]
+            sol = self.population[weight_index]
+            d = self.calculate_scalar(weight_index, offspring)
+            e = self.calculate_scalar(weight_index, sol)
+            if d < e:
+                self.population[weight_index] = self.toolbox.clone(offspring)
+
+    def calculate_scalar(self, weight_index, ind):
+        """Calculate scalar value for individual.
+
+        :params weight_index: An index of weight vector.
+        :params ind: An indivudlal to calculate scalar value.
+        :returns: Scalalized fitness value.
+        """
+        weight = self.weights[weight_index]
+        return self.scalarize_objective(weight, ind)
+
+    def weighted_sum(self, weight, ind):
+        """Calculate scalar value using Weighted Sum Approch.
+
+        :params weight: A weight vector.
+        :params ind: An indivudlal to calculate scalar value.
+        :returns: Scalalized fitness value.
+        """
+        return np.sum(weight * ind.fitness.values)
+
+    def tchebycheff(self, weight, ind):
+        """Calculate scalar value using Tchebycheff Approch.
+
+        :params weight: A weight vector.
+        :params ind: An indivudlal to calculate scalar value.
+        :returns: Scalalized fitness value.
+        """
+        max_fun = -1.0e+30
+        for f, z, w in zip(ind.fitness.values, self.ideal_point, weight):
+            diff = abs(f - z)
+            if w == 0:
+                feval = 0.00001 * diff
+            else:
+                feval = w * diff
+            if feval > max_fun:
+                max_fun = feval
+        return max_fun
+
+    def update_EP(self, offspring):
+        """Step 2.5) Update of EP.
+
+        :params offspring: An individual that maybe add to EP.
+        """
+
+        """
+        # Setup condition list for each objective
+        cond = [operator.gt if w > 0 else operator.lt
+                for w in offspring.fitness.weights]
+
+        nobj = len(cond)
+        def dominated(a, b):
+            # Dominated means all conditions are True
+            return sum([
+                c(a, b)
+                for a, b, c in zip(a.fitness.values, b.fitness.values, cond)
+            ]) == nobj
+
+        # Remove from EP if all vectors dominated by offspring
+        self.EP = [ind for ind in self.EP if not dominated(offspring, ind)]
+
+        # Add offspring to EP if no vectors in EP dominate offspring
+        for ind in self.EP:
+            if offspring == ind:
+                return
+            if dominated(ind, offspring):
+                return
+        self.EP.append(offspring)
+        """
+
+        # Above code is clear but very slow.
+        # So, speed up using numpy.
+
+        # np.array for empty list makes bad shape.
+        # So, first check and return.
+        if not self.EP:
+            self.EP.append(offspring)
+            return
+
+        # Setup condition list for each objective(like above but numpy)
+        cond = [np.greater if w > 0 else np.less
+                for w in offspring.fitness.weights]
+
+        nobj = len(offspring.fitness.weights)
+
+        # Make fitness array for EP and offspring
+        epf = np.array([ind.fitness.values for ind in self.EP])
+        osf = np.array(offspring.fitness.values)
+
+        # Evaluate condition for each objective
+        check = np.empty(epf.shape, dtype=np.bool)
+        for i, c in enumerate(cond):
+            check[:, i] = c(osf[i], epf[:, i])
+        # Dominated means all conditions are True
+        dominated = (check.sum(axis=1) == nobj)
+
+        # Remove from EP if all vectors dominated by offspring
+        self.EP = [ind for ind, dom in zip(self.EP, dominated) if not dom]
+
+        # If all EP are removed, simple add offspring.
+        if not self.EP:
+            self.EP.append(offspring)
+            return
+
+        # Don't add same gene.
+        # This is useful for binary code, but slow down for real gene(exact
+        # same gene is very rare).
+        # So, enable only not float gene.
+        if not type(offspring[0]) == float:
+            epg = np.array(self.EP)
+            osg = np.array(offspring)
+            # Count same gene and check that count is equal to length of gene.
+            same = ((epg == osg).sum(axis=1) == len(offspring))
+            # Same gene is already in EP. Don't add.
+            if np.any(same):
+                return
+
+        # Remove dominated idnividual from EP fitness array
+        epf = epf[~dominated]
+
+        # Like above dominate check but reverse operand
+        check = np.empty(epf.shape, dtype=np.bool)
+        for i, c in enumerate(cond):
+            check[:, i] = c(epf[:, i], osf[i])
+        dominated = check.sum(axis=1) == nobj
+        # Add offspring to EP if no vectors in EP dominate offspring
+        if dominated.sum() == 0:
+            self.EP.append(offspring)
diff --git a/deap/algorithms.py b/deap/algorithms.py
index b08897209..c97ca9266 100644
--- a/deap/algorithms.py
+++ b/deap/algorithms.py
@@ -28,6 +28,7 @@
 import random
 
 import tools
+import _algorithms
 
 
 def varAnd(population, toolbox, cxpb, mutpb):
@@ -501,3 +502,45 @@ def eaGenerateUpdate(toolbox, ngen, halloffame=None, stats=None,
             print logbook.stream
 
     return population, logbook
+
+
+def moead(population, toolbox, weights, neighbour_size, scalar_method, cxpb,
+          mutpb, ngen, stats=None, verbose=__debug__):
+    """This is algorithm implements the MOEA/D propsed in [Zhang2007]_.
+
+    :param population: A list of individuals.
+    :param toolbox: A :class:`~deap.base.Toolbox` that contains the evolution
+                    operators.
+    :param weights: A list of weight vectors.
+    :param neighbour_size: The number of neighbour.
+    :param scalar_method: The method to scalarize multi objective.
+                          Specify 'weightedSum' or 'tchebycheff'.
+    :param cxpb: The probability of mating two individuals.
+    :param mutpb: The probability of mutating an individual.
+    :param ngen: The number of generation.
+    :param stats: A :class:`~deap.tools.Statistics` object that is updated
+                  inplace, optional.
+    :param verbose: Whether or not to log the statistics.
+    :returns: The external population
+    :returns: A class:`~deap.tools.Logbook` with the statistics of the
+              evolution
+
+    This function expects the :meth:`toolbox.mate`, :meth:`toolbox.mutate`,
+    and :meth:`toolbox.evaluate` aliases to be registered in the toolbox.
+
+    You can register :meth:`toolbox.improve` for problem specific improvement.
+    ::
+
+        def custom_improve(i, offspring):
+            # Some problem specific improvement for subproblem i
+            return offspring
+
+        toolbox.register('improve', custom_improve)
+
+    .. [Zhang2007] Q. Zhang and H. Li: MOEA/D: A multi-objective evolutionary
+       algorithm based on decomposition, IEEE Trans. on Evolutionary
+       Computation, Vol. 11, No. 6, pp. 712-731, 2007
+    """
+    moead = _algorithms.MOEAD(population, toolbox, weights, neighbour_size,
+                              scalar_method, cxpb, mutpb, ngen, stats, verbose)
+    return moead.run()
diff --git a/deap/tests/test_algorithms.py b/deap/tests/test_algorithms.py
index 62e17a757..7b71ac4d9 100644
--- a/deap/tests/test_algorithms.py
+++ b/deap/tests/test_algorithms.py
@@ -231,3 +231,38 @@ def test_nsga3():
 
     for ind in pop:
         assert not (any(numpy.asarray(ind) < BOUND_LOW) or any(numpy.asarray(ind) > BOUND_UP))
+
+
+@with_setup(setup_func_multi_obj, teardown_func)
+def test_moead():
+    NDIM = 5
+    BOUND_LOW, BOUND_UP = 0.0, 1.0
+    MU = 16
+    NGEN = 100
+
+    toolbox = base.Toolbox()
+    toolbox.register("attr_float", random.uniform, BOUND_LOW, BOUND_UP)
+    toolbox.register("individual", tools.initRepeat, creator.__dict__[INDCLSNAME], toolbox.attr_float, NDIM)
+    toolbox.register("population", tools.initRepeat, list, toolbox.individual)
+
+    toolbox.register("evaluate", benchmarks.zdt1)
+    toolbox.register("mate", tools.cxSimulatedBinaryBounded, low=BOUND_LOW, up=BOUND_UP, eta=20.0)
+    toolbox.register("mutate", tools.mutPolynomialBounded, low=BOUND_LOW, up=BOUND_UP, eta=20.0, indpb=1.0/NDIM)
+
+    pop = toolbox.population(n=MU)
+    weights = []
+    for i in range(MU):
+        weights.append([i / float(MU - 1), (MU - i - 1) / float(MU - 1)])
+
+    EP, logbook = algorithms.moead(
+        pop, toolbox, weights, 4, 'tchebycheff',
+        0.9, 1.0, NGEN, None, False
+    )
+
+    hv = hypervolume(EP, [11.0, 11.0])
+    # hv = 120.777 # Optimal value
+
+    assert hv > HV_THRESHOLD, "Hypervolume is lower than expected %f < %f" % (hv, HV_THRESHOLD)
+
+    for ind in EP:
+        assert not (any(numpy.asarray(ind) < BOUND_LOW) or any(numpy.asarray(ind) > BOUND_UP))
diff --git a/examples/ga/moead.py b/examples/ga/moead.py
new file mode 100644
index 000000000..91719bc48
--- /dev/null
+++ b/examples/ga/moead.py
@@ -0,0 +1,110 @@
+#    This file is part of DEAP.
+#
+#    DEAP is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU Lesser General Public License as
+#    published by the Free Software Foundation, either version 3 of
+#    the License, or (at your option) any later version.
+#
+#    DEAP is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#    GNU Lesser General Public License for more details.
+#
+#    You should have received a copy of the GNU Lesser General Public
+#    License along with DEAP. If not, see <http://www.gnu.org/licenses/>.
+
+import array
+import random
+import json
+
+import numpy
+
+from math import sqrt
+
+from deap import algorithms
+from deap import base
+from deap import benchmarks
+from deap.benchmarks.tools import diversity, convergence, hypervolume
+from deap import creator
+from deap import tools
+
+creator.create("FitnessMin", base.Fitness, weights=(-1.0, -1.0))
+creator.create("Individual", array.array, typecode='d', fitness=creator.FitnessMin)
+
+toolbox = base.Toolbox()
+
+# Problem definition
+# Functions zdt1, zdt2, zdt3, zdt6 have bounds [0, 1]
+BOUND_LOW, BOUND_UP = 0.0, 1.0
+
+# Functions zdt4 has bounds x1 = [0, 1], xn = [-5, 5], with n = 2, ..., 10
+# BOUND_LOW, BOUND_UP = [0.0] + [-5.0]*9, [1.0] + [5.0]*9
+
+# Functions zdt1, zdt2, zdt3 have 30 dimensions, zdt4 and zdt6 have 10
+NDIM = 30
+
+def uniform(low, up, size=None):
+    try:
+        return [random.uniform(a, b) for a, b in zip(low, up)]
+    except TypeError:
+        return [random.uniform(a, b) for a, b in zip([low] * size, [up] * size)]
+
+toolbox.register("attr_float", uniform, BOUND_LOW, BOUND_UP, NDIM)
+toolbox.register("individual", tools.initIterate, creator.Individual, toolbox.attr_float)
+toolbox.register("population", tools.initRepeat, list, toolbox.individual)
+
+toolbox.register("evaluate", benchmarks.zdt1)
+toolbox.register("mate", tools.cxSimulatedBinaryBounded, low=BOUND_LOW, up=BOUND_UP, eta=20.0)
+toolbox.register("mutate", tools.mutPolynomialBounded, low=BOUND_LOW, up=BOUND_UP, eta=20.0, indpb=1.0/NDIM)
+
+def main(seed=None):
+    random.seed(seed)
+
+    NGEN = 250
+    MU = 100
+    CXPB = 0.9
+    MUTPB = 1.0
+    NEIGHBOUR_SIZE = 10
+
+    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)
+    
+    pop = toolbox.population(n=MU)
+    weights = []
+    for i in range(MU):
+        weights.append([i / float(MU - 1), (MU - i - 1) / float(MU - 1)])
+
+    EP, logbook = algorithms.moead(
+        pop, toolbox, weights, NEIGHBOUR_SIZE, 'tchebycheff',
+        CXPB, MUTPB, NGEN, stats, True
+    )
+
+    print("Final population hypervolume is %f" % hypervolume(EP, [11.0, 11.0]))
+
+    return EP, logbook
+        
+if __name__ == "__main__":
+    with open("pareto_front/zdt1_front.json") as optimal_front_data:
+        optimal_front = json.load(optimal_front_data)
+    # Use 500 of the 1000 points in the json file
+    optimal_front = sorted(optimal_front[i] for i in range(0, len(optimal_front), 2))
+    
+    pop, stats = main()
+    pop.sort(key=lambda x: x.fitness.values)
+    
+    print(stats)
+    print("Convergence: ", convergence(pop, optimal_front))
+    print("Diversity: ", diversity(pop, optimal_front[0], optimal_front[-1]))
+    
+    import matplotlib.pyplot as plt
+    import numpy
+    
+    front = numpy.array([ind.fitness.values for ind in pop])
+    optimal_front = numpy.array(optimal_front)
+    plt.scatter(optimal_front[:,0], optimal_front[:,1], c="r")
+    plt.scatter(front[:,0], front[:,1], c="b")
+    plt.axis("tight")
+    plt.show()
diff --git a/examples/ga/moead_onemax.py b/examples/ga/moead_onemax.py
new file mode 100644
index 000000000..e78f2d4ba
--- /dev/null
+++ b/examples/ga/moead_onemax.py
@@ -0,0 +1,86 @@
+#    This file is part of DEAP.
+#
+#    DEAP is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU Lesser General Public License as
+#    published by the Free Software Foundation, either version 3 of
+#    the License, or (at your option) any later version.
+#
+#    DEAP is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#    GNU Lesser General Public License for more details.
+#
+#    You should have received a copy of the GNU Lesser General Public
+#    License along with DEAP. If not, see <http://www.gnu.org/licenses/>.
+
+import array
+import random
+
+import numpy
+import matplotlib.pyplot as plt
+
+from deap import base
+from deap import creator
+from deap import tools
+from deap import algorithms
+
+creator.create("FitnessMax", base.Fitness, weights=(1.0, 1.0))
+creator.create("Individual", array.array, typecode='b', fitness=creator.FitnessMax)
+
+toolbox = base.Toolbox()
+
+GENE_SIZE = 100
+ZERO_COUNT = 10
+
+toolbox.register("attr_bool", random.randint, 0, 1)
+toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_bool, GENE_SIZE)
+toolbox.register("population", tools.initRepeat, list, toolbox.individual)
+
+def calc(individual, check):
+    diff = 0
+    for i in range(len(individual)):
+        if individual[i] != check[i]:
+            diff += 1
+    return diff
+
+CHECK1 = [0] * 0 + [1] * GENE_SIZE
+CHECK2 = [0] * ZERO_COUNT + [1] * (GENE_SIZE - ZERO_COUNT)
+
+def evalMoOneMax(individual):
+    l = len(individual)
+    return l - calc(individual, CHECK1), l - calc(individual, CHECK2)
+
+toolbox.register("evaluate", evalMoOneMax)
+toolbox.register("mate", tools.cxTwoPoint)
+toolbox.register("mutate", tools.mutFlipBit, indpb=0.01)
+
+MU = 100
+NEIGHBOUR_SIZE = 10
+NGEN = 100
+
+CXPB = 0.5
+MUTPB = 0.2
+
+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)
+
+random.seed(64)
+
+pop = toolbox.population(n=MU)
+weights = []
+for i in range(MU):
+    weights.append([i / float(MU - 1), (MU - i - 1) / float(MU - 1)])
+
+EP, logbook = algorithms.moead(pop, toolbox, weights,
+                               neighbour_size=NEIGHBOUR_SIZE,
+                               scalar_method='tchebycheff',
+                               cxpb=CXPB, mutpb=MUTPB, ngen=NGEN,
+                               stats=stats, verbose=True)
+
+front = numpy.array([ind.fitness.values for ind in EP])
+plt.scatter(front[:,0], front[:,1], c="b")
+plt.axis("tight")
+plt.show()