Multiobjective Optimization#

This section provides code implementations for the concepts covered in class for multiobjective optimization. Multiobjective optimization provides methods to solve optimization problems with multiple competing objectives. The application of these methods to simple analytical examples is also provided in this section. The multiobjective methods used in class and in this section of the Jupyter book are based on differential evolution (DE) covered earlier in the class. This section has the following two subsections:

  1. Multiobjective optimization using differential evolution

  2. Multiobjective optimization using Kriging models

The code blocks below introduce the multiobjective optimization examples used in this section. The next code block imports the required packages for this section.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as colors
from smt.sampling_methods import LHS
from pymoo.util.nds.non_dominated_sorting import NonDominatedSorting
import pymoo.gradient.toolbox as anp
from pymoo.core.problem import Problem

Branin-Currin optimization problem#

The first example used in this section is an unconstrained multiobjective problem with two design variables and two objective functions. The two functions are the rescaled Branin function and the Currin function. The optimization problem statement is given as

\[\begin{split} \begin{gathered} \min f(\textbf{x}) = \begin{cases} f_1(\textbf{x}) = \frac{1}{51.95}(( \bar{x_2} - \frac{5.1}{4\pi^2} \bar{x_1}^2 + \frac{5}{\pi}\bar{x_1} - 6 )^2 + 10 ( 1-\frac{1}{8\pi} )\cos \bar{x_1} - 44.81) \\ \\ f_2(\textbf{x}) = [1-\exp(\frac{-1}{2x_2})]\frac{2300x_1^3+1900x_1^2+2092x_1+60}{100x_1^3+500x_1^2+4x_1+20} \end{cases} \\ 0.0 \leq x_1, x_2 \leq 1.0 \\ \text{where} \quad \bar{x_1} = 15x_1 - 5, \bar{x_2} = 15x_2 \end{gathered} \end{split}\]

The block of code below defines the two functions.

# Defining the objective functions
def branin(x):

    dim = x.ndim
    if dim == 1:
        x = x.reshape(1,-1)

    x1 = 15*x[:,0] - 5
    x2 = 15*x[:,1]

    b = 5.1 / (4*np.pi**2)
    c = 5 / np.pi
    t = 1 / (8*np.pi)

    y = (1/51.95)*((x2 - b*x1**2 + c*x1 - 6)**2 + 10*(1-t)*np.cos(x1) + 10 - 44.81)

    if dim == 1:
        y = y.reshape(-1)

    return y

def currin(x):

    dim = x.ndim
    if dim == 1:
        x = x.reshape(1,-1)

    x1 = x[:,0]
    x2 = x[:,1]
    
    factor = 1 - np.exp(-1/(2*x2))
    num = 2300*x1**3 + 1900*x1**2 + 2092*x1 + 60
    den = 100*x1**3 + 500*x1**2 + 4*x1 + 20
    y = factor*num/den
    
    if dim == 1:
        y = y.reshape(-1)

    return y

The next few code blocks will locate the Pareto front for the problem and plot the Pareto front in the objective space. This is done by generating a mesh of points within the bounds of the problem and sorting the points to locate the non-dominated points. The code block uses the non-dominated sorting algorithm of pymoo to find the non-dominated points in the objective space. The non-dominated points represent the Pareto front of the problem. The non-dominated sorting algorithm may not always perfectly provide all of the non-dominated points of the problem and it is a good idea to make sure that all the points provided are actually non-dominated solutions. It is also sometimes necessary to adjust the parameter n_stop_if_ranked which indicates approximately how many points should survive in the initial population of points at the end of the sorting algorithm. Raising or lowering this value can improve the representation of the Pareto front.

# Generating a grid of points
num_points = 100

# Defining x and y values
x = np.linspace(1e-6,1,num_points)
y = np.linspace(1e-6,1,num_points)

# Creating a mesh
X, Y = np.meshgrid(x, y)

# Finding the front through non-dominated sorting
z1 = branin(np.hstack((X.reshape(-1,1),Y.reshape(-1,1))))
z2 =  currin(np.hstack((X.reshape(-1,1),Y.reshape(-1,1))))
nds = NonDominatedSorting()
F = np.column_stack((z1,z2))
pareto = nds.do(F, n_stop_if_ranked=50)

# Plotting the non-dominated solutions 
Z_pareto = F[pareto[0]]
fig, ax = plt.subplots(1,1,figsize=(8,6))
ax.scatter(Z_pareto[1:,0], Z_pareto[1:,1], color="red", label="Pareto Front")
ax.set_xlabel("$f_1$", fontsize=14)
ax.set_ylabel("$f_2$", fontsize=14)
ax.grid()
ax.legend(loc="upper right", fontsize = 14)

<matplotlib.legend.Legend at 0x13c52c250>
../_images/e13715a4a8bcb1b3dc96d0bed97dfddf9af854d7c50915ab3869f3846e020f28.png

The above plot shows the non-dominated solutions of the problem in the objective space. The curve passing through these points is the Pareto front of the problem. The Pareto front is convex in shape for this problem. This plot gives good insight into the Pareto front of the problem, however, simply sorting the points from the mesh will not reveal the true front and directly solving the optimization problem will be a better method for a complex problem such as this.

Constrained multiobjective optimization problem#

The second example used in this section is a constrained multiobjective problem with two design variables, two objective functions and two constraints. The optimization problem statement is given as

\[\begin{split} \begin{gathered} \min f(\textbf{x}) = \begin{cases} f_1(\textbf{x}) = 4x_1^2+4x_2^2 \\ f_2(\textbf{x}) = (x_1-5)^2 + (x_2-5)^2 \end{cases} \\ \textrm{subject to} \quad g_1(\textbf{x}) = (x_1-5)^2 + x_2^2 -25 \leq 0\\ g_2(\textbf{x}) = 7.7 - ((x_1-8)^2 + (x_2+3)^2) \leq 0\\ -20 \leq x_1,x_2 \leq 20\\ \end{gathered} \end{split}\]

The block of code below defines the objective and constraint functions for the problem.

# Defining the objective functions
def f1(x):
    
    dim = x.ndim
    if dim == 1:
        x = x.reshape(1,-1)

    y = 4*x[:,0]**2 + 4*x[:,1]**2
    
    return y
    
def f2(x):

    dim = x.ndim
    if dim == 1:
        x = x.reshape(1,-1)

    y = (x[:,0]-5)**2 + (x[:,1]-5)**2

    return y

def g1(x):
    
    dim = x.ndim
    if dim == 1:
        x = x.reshape(1,-1)
    
    g = (x[:,0]-5)**2 + x[:,1]**2 - 25

    return g

def g2(x):

    dim = x.ndim
    if dim == 1:
        x = x.reshape(1,-1)
    
    g = 7.7 - ((x[:,0]-8)**2 + (x[:,1]+3)**2)

    return g

The code block below uses the non-dominated sorting algorithm of pymoo to find the non-dominated points in the objective space. The non-dominated points represent the Pareto front of the problem. The obtained Pareto front is then plotted in the objective space of the problem.

# Generating a grid of points
num_points = 100

# Defining x and y values
x = np.linspace(-20,20,num_points)
y = np.linspace(-20,20,num_points)

# Creating a mesh
X, Y = np.meshgrid(x, y)

# Finding the front through non-dominated sorting
z1 = f1(np.hstack((X.reshape(-1,1),Y.reshape(-1,1))))
z2 = f2(np.hstack((X.reshape(-1,1),Y.reshape(-1,1))))
const1 = g1(np.hstack((X.reshape(-1,1),Y.reshape(-1,1))))
const2 = g2(np.hstack((X.reshape(-1,1),Y.reshape(-1,1))))
z1 = z1[(const1<0) & (const2<0)]
z2 = z2[(const1<0) & (const2<0)]
nds = NonDominatedSorting()
F = np.column_stack((z1,z2))
pareto = nds.do(F, n_stop_if_ranked=50)
Z_pareto = F[pareto[0]]

# Plotting the contours
fig, ax = plt.subplots(1,1,figsize=(8,6))
ax.scatter(Z_pareto[:,0], Z_pareto[:,1], color="red", label="Pareto Front")
ax.set_xlabel("$f_1$", fontsize=14)
ax.set_ylabel("$f_2$", fontsize=14)
ax.grid()
ax.legend(loc="upper right", fontsize = 14)

<matplotlib.legend.Legend at 0x13c9bf970>
../_images/c2f152be5e982cc40aa613a40b00d5101f8ce0f7305850b98a8d02b689da439a.png

The Pareto front of the problem is shown in the plot above. Since the functions involved in the problem are simple, it is likely that the sorted non-dominated points represent a fairly accurate Pareto front for the problem. The Pareto front is continuous and convex for this problem.