# Lecture 12, Methods for multiobjective optimization 

If you want to know more about the topic of this lecture, I urge you to read Professor Miettinen's book Nonlinear Multiobjective Optimization



## Classification of methods

Methods for multiobjective optimization are often charactermized by the involvement of the decision maker in the process.

The types of methods are
* **no preference methods**, where the decision maker does not play a role,
* **a priori methods**, where the decision maker gives his/her preference information at first and then the optimization method find the best match to that preference information,
* **a posteriori methods**, where the optimization methods try to characterize all/find a good represenatation of the Pareto optimal solutions and the decision maker chooses the most preferred one of those,
* **interactive methods**, where the optimization method and the decision maker alternate in iteratively search for the most preferred solution.

##  Our example problem for this lecture

We study a hypothetical decision problem of buying a car, when you can choose to have a car with power between (denoted by $p$) 50 and 200 kW and average consumption (denoted by $c$) per 100 km between 3 and 10 l. However, in addition to the average consumption and power, you need to decide the volume of the cylinders (v), which may be between 1000 $cm^3$ and $4000cm^3$. Finally, the price of the car follows now a function 
$$
\left(\sqrt{\frac{p-50}{50}}\\
+\left(\frac{p-50}{50}\right)^2+0.3(10-c)\\ +10^{-5}\left(v-\left(1000+3000\frac{p-50}{150}\right)\right)^2\right)10000\\+5000
$$
in euros. This problem can be formulated as a multiobjective optimization problem
$$
\begin{align}
\min \quad & \{c,-p,P\},\\
\text{s.t. }\quad
&50\leq p\leq 200\\
&3\leq c\leq 10\\
&1000\leq v\leq 4000,\\
\text{where }\quad&P = \left(\sqrt{\frac{p-50}{50}}+\frac{p-50}{50}^2+\frac{10-c}{10} +\right.\\
&\left.\frac{p-50}{50}\frac{10-c}{10}\right)10000+5000
\end{align}
$$


```python
#Let us define a Python function which returns the value of this
import math
def car_problem(c,p,v):
#    import pdb; pdb.set_trace()
    return [#Objective function values
        c,-p,
        (math.sqrt((p-40.)/50.)+((p-50.)/50.)**2+
        0.3*(10.-c)+0.00001*(v-(1000.+3000.*(p-50.)/150.))**2)*10000.
        +5000.] 
```


```python
print("Car with 3 l/100km consumption, 50kW and 1000$cm^3$ engine would cost "
      +str(car_problem(3,50,1000)[2])+"€")
print("Car with 3 l/100km consumption, 100kW and 2000$cm^3$ engine would cost "
      +str(car_problem(3,100,2000)[2])+"€")
print("Car with 3 l/100km consumption, 100kW and 1000$cm^3$ engine would cost "
      +str(car_problem(3,100,1000)[2])+"€")
```

    Car with 3 l/100km consumption, 50kW and 1000$cm^3$ engine would cost 30472.135955€
    Car with 3 l/100km consumption, 100kW and 2000$cm^3$ engine would cost 46954.4511501€
    Car with 3 l/100km consumption, 100kW and 1000$cm^3$ engine would cost 146954.45115€


## Normalization of the objectives

**In many of the methods, the normalization of the objectives is necessary. **

We can normalize the objectives using the nadir and ideal and setting the normalized objective as
$$ \tilde f_i = \frac{f_i-z_i^{ideal}}{z_i^{nadir}-z_i^{ideal}}$$

## Calculating the ideal

**Finding the ideal for problems is usually easy, if you can optimize the objective functions separately.**

For the car problem, ideal can be computed easily using the script:


```python
#Calculating the ideal
from scipy.optimize import minimize
import ad
def calc_ideal(f):
    ideal = [0]*3 #Because three objectives
    solutions = [] #list for storing the actual solutions, which give the ideal
    bounds = ((3,10),(50,200),(1000,4000)) #Bounds of the problem
    for i in range(3):
        res=minimize(
            #Minimize each objective at the time
            lambda x: f(x[0],x[1],x[2])[i], [3,50,1000], method='SLSQP'
            #Jacobian using automatic differentiation
            ,jac=ad.gh(lambda x: f(x[0],x[1],x[2])[i])[0]
            #bounds given above
            ,bounds = bounds
            ,options = {'disp':True, 'ftol': 1e-20, 'maxiter': 1000})
        solutions.append(f(res.x[0],res.x[1],res.x[2]))
        ideal[i]=res.fun
    return ideal,solutions
```


```python
ideal, solutions= calc_ideal(car_problem)
print ("ideal is "+str(ideal))
```

    Optimization terminated successfully.    (Exit mode 0)
                Current function value: 3.0
                Iterations: 1
                Function evaluations: 12
                Gradient evaluations: 1
    Positive directional derivative for linesearch    (Exit mode 8)
                Current function value: -200.000007903
                Iterations: 14
                Function evaluations: 10
                Gradient evaluations: 10
    Positive directional derivative for linesearch    (Exit mode 8)
                Current function value: 9371.55596176
                Iterations: 10
                Function evaluations: 6
                Gradient evaluations: 6
    ideal is [3.0, -200.00000790345445, 9371.5559617552626]


## Pay-off table method

**Finding the nadir value is however, usually much harder.**

Usually, the nadir value is estimated using the so-called pay-off table method.

The pay-off table method does not guarantee to find the nadir for problems with more than two objectives. (One of your exercises this week will be to show this.) 

The method is, however, a generally accepted way of approximating the nadir vector.

In the pay-off table method:
1. the objective values for attaining the individual minima are added in table
2. the nadir is estimated by each objectives maxima in the table.

### The nadir for the car selection problem
The table now becomes by using the *solutions* that we returned while calculating the ideal


```python
for solution in solutions:
    print solution
```

    [3.0, -50.0, 30472.13595499958]
    [3.0, -200.00000790345445, 1033888.6485874163]
    [10.033526664414772, -50.0, 9371.5559617552626]


Thus, the esimation of the nadir vector is 
$$(10,-50,1033888.543820).$$

This is actually the real Nadir vector for this problem.

### Normalized car problem


```python
#Let us define a Python function which returns the value of this
import math
def car_problem_normalized(c,p,v):
    z_ideal = [3.0, -200.0, 9472.1359549995796]
    z_nadir = [10,-50,1033888.543820]
#    import pdb; pdb.set_trace()
    z = car_problem(c,p,v) 
    return [(zi-zideali)/(znadiri-zideali) for 
            (zi,zideali,znadiri) in zip(z,z_ideal,z_nadir)]
```


```python
print("Normalized value of the car problem at (3,50,1000) is "
      +str(car_problem_normalized(3,50,1000)))
print("Normalized value of the car problem at (3,125,2500) is "
      +str(car_problem_normalized(3,125,2500)))
print("Normalized value of the car problem at (10,100,1000) is "
      +str(car_problem_normalized(10,100,1000)))
```

    Normalized value of the car problem at (3,50,1000) is [0.0, 1.0, 0.020499476422645723]
    Normalized value of the car problem at (3,125,2500) is [0.0, 0.5, 0.05082529765792966]
    Normalized value of the car problem at (10,100,1000) is [1.0, 0.6666666666666666, 0.11370602257129601]


** So, value 1 now indicates the worst value on the Pareto frontier and value 0 indicates the best values**

Let's set the ideal and nadir for later reference:


```python
z_ideal = [3.0, -200.0, 9472.1359549995796]
z_nadir = [10.,-50,1033888.543820]
```

**From now on, we will deal with the normalized problem, although, we write just $f$.** The aim of this is to simplify presentation.

## No preference methods

* Usually only for situations, where the decision maker is not available or does not want to get involved
* These just compute a single Pareto optimal solution, which is in somehow mathematically thought as the best compromise

## Notation

For short, let us denote the feasible set $\{x\in\mathbb R^n: g_j(x) \geq 0 \text{ for all }j=1,\ldots,J \text{ and } h_k(x) = 0\text{ for all }k=1,\ldots,K\}$ by $S$.

### Method of Global criterion

Involved minimization of the p-norm of $f(x)-z^{ideal}$, that is we solve the problem
$$
\min_{x\in S} \|f(x) - z^{ideal}\|_p.
$$



**An optimal solution to this problem is Pareto optimal**

### Applying to our problem studied



```python
import numpy as np
def global_criterion_method(f,p):
    #Ideal for any normalized problem is (0,0,0)
    ideal = [0,0,0]
    bounds = ((3,10),(50,200),(1000,4000)) #Bounds of the problem
    #Objective is the norm of objective function values minus
    #the ideal
    obj = lambda x: np.linalg.norm(np.array(f(x[0],x[1],x[2]))
                                   -np.array(ideal),ord=p)
    res=minimize(
            #Minimize p distance from the ideal
            obj, [3,50,1000], method='SLSQP'
            #Jacobian using automatic differentiation
            ,jac=ad.gh(obj)[0]
            #bounds given above
            ,bounds = bounds,options = {'disp':True})
    return res.x
```


```python
%pdb
```

    Automatic pdb calling has been turned ON
    



```python
solution = global_criterion_method(car_problem_normalized,
                                           2)
print "variable values are ",solution
f_global_criterion = car_problem(solution[0],solution[1],solution[2])
print "objective values are ",f_global_criterion
```

    Optimization terminated successfully.    (Exit mode 0)
                Current function value: 0.549669334387
                Iterations: 15
                Function evaluations: 15
                Gradient evaluations: 15
     variable values are  [    3.05211179   138.12821175  1000.00729216]
    objective values are  [3.0521117947829333, -138.12821174914697, 381579.84409750014]


## A posteriori methods

* A posteriori methods generate a representation of the Pareto optimal solutions, or the complete set of Pareto optimal solutions
* Benefits
  * The solutions can be visualized for problems with 2 or 3 objectives so the decision making is possible
  * When succesful, they give an understanding of the Pareto front
* Drawbacks
  * Approximating the Pareto optimal set often time-consuming
  * Decision making from a large representation may be very difficut


### The weighting method

Based on solving optimization problem
$$
\min_{x\in S} \sum_{i=1}^kw_if_i(x)
$$
for different weights $w_i\geq0$, $i=1,\ldots,k$. 

**The idea is to generate weights evenly and then have evenly spread solutions.**

**An optimal solution the weighted problem is Pareto optimal, if all the weights $w_i$ are $>0$.**



### Application to our problem


```python
import numpy as np
def weighting_method(f,w):
    points = []
    bounds = ((3,10),(50,200),(1000,4000)) #Bounds of the problem
    for wi in w:
        res=minimize(
            #weighted sum
            lambda x: sum(np.array(wi)*np.array(f(x[0],x[1],x[2]))), 
            [3,50,1000], method='SLSQP'
            #Jacobian using automatic differentiation
            ,jac=ad.gh(lambda x: sum(np.array(wi)*np.array(f(x[0],x[1],x[2]))))[0]
            #bounds given above
            ,bounds = bounds,options = {'disp':False})
        points.append(res.x)
    return points
```


```python
w = np.random.random((500,3)) #500 random weights
repr = weighting_method(car_problem_normalized,w)
```


```python
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
f_repr_ws = [car_problem(repri[0],repri[1],repri[2]) for repri in repr]
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter([f[0] for f in f_repr_ws],[f[1] for f in f_repr_ws],[f[2] for f in f_repr_ws])
plt.show()
```

### Epsilon-constraint method



Based on solving optimization problem
$$
\begin{align}
\min \quad &f_j(x)\\
\text{s.t. }\quad &x\in S\\
&f_i(x)\leq \epsilon_i \text{ for all }i\neq j
\end{align}
$$
for different bounds $\epsilon_i, $i\neq j$. 

**The idea is to generate $\epsilon$ evenly within the bounds of the ideal and nadir vectors and then have evenly spread solutions.**

**A solution is Pareto optimal $x^*$, if it is the solution to the epsilon constraint problem for all $j=1,\ldots,k$ and $\epsilon = f(x^*)$.**



### Application to our problem


```python
import numpy as np
from scipy.optimize import minimize
import ad
def e_constraint_method(f,eps,z_ideal,z_nadir):
    points = []
    for epsi in eps:
        bounds = ((3,epsi[0]*(z_nadir[0]-z_ideal[0])+z_ideal[0]),
                  (-1.*(epsi[1]*(z_nadir[1]-z_ideal[1])+z_ideal[1]),
                   200),(1000,4000)) #Added bounds for two first objectives
        res=minimize(
            #Third objective
            lambda x: f(x[0],x[1],x[2])[2], 
            [3,200,1000], method='SLSQP'
            #Jacobian using automatic differentiation
            ,jac=ad.gh(lambda x: f(x[0],x[1],x[2])[2])[0]
            #bounds given above
            ,bounds = bounds,options = {'disp':False})
        if res.success:
            points.append(res.x)
    return points
```


```python
eps = np.random.random((100,2))
repr = e_constraint_method(car_problem_normalized,eps,z_ideal,z_nadir)
```


```python
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
f_repr_eps = [car_problem(repri[0],repri[1],repri[2]) for repri in repr]
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter([f[0] for f in f_repr_eps],[f[1] for f in f_repr_eps],[f[2] for f in f_repr_eps])
plt.show()
```

## Comparison of the weighted sum method and the epsilon constraint method


```python
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
f_repr_eps = [car_problem(repri[0],repri[1],repri[2]) for repri in repr]
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
#Mark solutions to epsilon constsraint problem using crosses
ax.scatter([f[0] for f in f_repr_eps],[f[1] for f in f_repr_eps],
           [f[2] for f in f_repr_eps],marker='x')
#Mark solutions to weighted sum problem using dots
ax.scatter([f[0] for f in f_repr_ws],[f[1] for f in f_repr_ws],
           [f[2] for f in f_repr_ws],marker='o')
plt.show()
```

**The weighting method can find all the Pareto optimal solutions only, when the objective functions are convex and the feasible set $S$ is convex.** 

**The weighting method produces very unevenly spread Pareto optimal solutions, even when the problem is convex. ** 

**The epsilon constraint method, however, adds constraints to the problem, which may make it much harder to solve**

## A priori methods

* A priori methods ask for preferences from the decision maker, and then find the Pareto optimal solution that best matches these preferences
* Benefits
  * If the decision maker knows what he/she wants and understands the preference information asked for, then application is fast
* Drawbacks
  * The decision maker may not know what he/she wants, because he does not know the Pareto optimal solutions
  * The decision maker may not understand how the preferences he/she gives affect the solutions found

## Achievement scalarizing problem

There are multiple versions of the achievement scalarizing problem, but all of them are based on a refence point.

A reference point
$$z^{ref} = (z^{ref}_1,\ldots,z^{ref}_k)$$
contains preferable values (so-called aspiration levels) for the objectives.

Then the achievement scalarizing problem maps this point and a feasible solution to the multiobjective problem to a scalar (i.e., scalarizes it). One of the most commonly used is
$$
\min_{x\in S}\left( \max_{i=1}^k(f_i(x)-z_i^{ref}) +\rho\sum_{i=1}^nf_i(x)\right)
$$
where $\rho>0$ is a small value. The second part is called an augmentation term.

** The solution to the problem is guaranteed to be Pareto optimal **

** Any (properly) Pareto optimal solution can be found with some reference point**



### Application to our car problem


```python
import numpy as np
from scipy.optimize import minimize
import ad
def asf(f,ref,z_ideal,z_nadir,rho):
    bounds = ((3,10),(50,200),(1000,4000)) #Bounds of the problem
    #Normalizing the reference point
    ref_norm = [(refi-z_ideali)/(z_nadiri-z_ideali) 
                for (refi,z_ideali,z_nadiri) in zip(ref,z_ideal,z_nadir)]
    def obj(x):
        return np.max(np.array(f(x[0],x[1],x[2]))-ref_norm)\
           +rho*np.sum(f(x[0],x[1],x[2]))
    res=minimize(
        #Objective function defined above
        obj, 
        [3,200,1000], method='SLSQP'
        #Jacobian using automatic differentiation
        ,jac=ad.gh(obj)[0]
        #bounds given above
        ,bounds = bounds,options = {'disp':True, 'ftol': 1e-20,
                                    'maxiter': 1000})
    return res
```


```python
rho = 0.000001
#The reference point for the problem
ref =  [5,-200,300000] #To be added at the class
res = asf(car_problem_normalized,ref,z_ideal,z_nadir,rho)
print "Solution is ",res.x
print "Objective function values are ",car_problem(res.x[0],res.x[1],res.x[2])
```

    Positive directional derivative for linesearch    (Exit mode 8)
                Current function value: 1.36797651381e-07
                Iterations: 100
                Function evaluations: 271
                Gradient evaluations: 96
    Solution is  [  3.00000000e+00   2.00000000e+02   3.60349841e+03]
    Objective function values are  [3.000000000000008, -200.00000000000009, 149609.89516963021]



```python
car_problem(10,149.70923875,3000)
```

## Interactive methods

* Interactive methods iteratively search for the preferred solution with decision maker and optimization alternating
* Benefits
  * Decision maker gets to learn about
    * the available solutions, and
    * how preferences affect the solutions found
  * Computation is less intensive, because no need to generate a large representation of Pareto optimal solutions
* Drawbacks
  * Needs active involvement from the decision maker
  * If the proble is computationally expensive, then the decision maker may need to wait a long time between solutions


## Interactive methods (cont)

**Interactive methods are one of the main research areas here at the Industrial optimization research group**

We will study interactive methods using the IND-NIMBUS software developed at the research group.

The IND-NIMBUS software contains many different interactive methods.


### The synchronous NIMBUS Method

Introduced by Profesors Kaisa Miettinen and Marko Mäkelä in 2006.

Is based on classification of objectives into 
* those that should be improved as much as possible,
* those that should be improved until a given limit,
* those that are acceptable at the moment,
* those that can be allowed to worsen until a given limit,
* those that are allowed to move freely at the moment.

Synchronous NIMBUS is based on representing this information as four different single-objective optimization problems that are then solved and solutions are shown to the decision maker.



The Synchronous NIMBUS is implemented as a part of IND-NIMBUS (http://ind-nimbus.it.jyu.fi/) software framework. 

There also exists a WWW-NIMBUS (https://wwwnimbus.it.jyu.fi/) web implementation (has not been maintained in approx 10 years, but works).

4c30b0d814a8de0d443d454db4aa124aef3eb5d2
ipynb
Jupyter Notebook
Lecture 12, Methods for multiobjective optimization.ipynb
maeehart/TIES483
cce5c779aeb0ade5f959a2ed5cca982be5cf2316
CC-BY-3.0
2019-04-26T12:46:14.000Z
2021-11-23T03:38:59.000Z
Lecture 12, Methods for multiobjective optimization.ipynb
maeehart/TIES483
cce5c779aeb0ade5f959a2ed5cca982be5cf2316
CC-BY-3.0
Lecture 12, Methods for multiobjective optimization.ipynb
maeehart/TIES483
cce5c779aeb0ade5f959a2ed5cca982be5cf2316
CC-BY-3.0
2016-01-08T16:28:11.000Z
2021-04-10T05:18:10.000Z
Qwen/Qwen-72B
1. YES
2. YES