<jupyter_start><jupyter_text>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![Nonlinear Multiobjective Optimization](images/Miettinen2.gif)  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 lectureWe 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}$$<jupyter_code>#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.] 
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])+"€")<jupyter_output>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€<jupyter_text>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:<jupyter_code>#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
ideal, solutions= calc_ideal(car_problem)
print ("ideal is "+str(ideal))<jupyter_output>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]<jupyter_text>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 table2. the nadir is estimated by each objectives maxima in the table.  The nadir for the car selection problemThe table now becomes by using the *solutions* that we returned while calculating the ideal<jupyter_code>for solution in solutions:
    print solution<jupyter_output>[3.0, -50.0, 30472.13595499958]
[3.0, -200.00000790345445, 1033888.6485874163]
[10.033526664414772, -50.0, 9371.5559617552626]<jupyter_text>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<jupyter_code>#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)]
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)))<jupyter_output>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]<jupyter_text>** 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:<jupyter_code>z_ideal = [3.0, -200.0, 9472.1359549995796]
z_nadir = [10.,-50,1033888.543820]<jupyter_output><empty_output><jupyter_text>**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 NotationFor 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 criterionInvolved 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.$$![alt text](images/mgc.svg "Method of global criterion")**An optimal solution to this problem is Pareto optimal**  Applying to our problem studied<jupyter_code>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
%pdb
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<jupyter_output>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]<jupyter_text>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 methodBased 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$.**![alt text](images/ws.svg "Weighting method")  Application to our problem<jupyter_code>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
w = np.random.random((500,3)) #500 random weights
repr = weighting_method(car_problem_normalized,w)
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()<jupyter_output><empty_output><jupyter_text>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^*)$.**![alt text](images/eps.svg "Epsilon constraint method")  Application to our problem<jupyter_code>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
eps = np.random.random((100,2))
repr = e_constraint_method(car_problem_normalized,eps,z_ideal,z_nadir)
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()<jupyter_output><empty_output><jupyter_text>Comparison of the weighted sum method and the epsilon constraint method<jupyter_code>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()<jupyter_output><empty_output><jupyter_text>**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 problemThere 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**![alt text](images/ach.svg "Achievement scalarizing method")  Application to our car problem<jupyter_code>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
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])
car_problem(10,149.70923875,3000)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11373377957841584, 'fraction_numerical': 0.06458810404556398, 'mean_word_length': 4.067539422790836, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 38, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '639104', 'license': 'CC-BY-3.0', 'n_tokens_mistral': 6145, 'n_tokens_neox': 5271, 'n_words': 1814, 'path': 'Lecture 12, Methods for multiobjective optimization.ipynb', 'repo_name': 'maeehart/TIES483', 'size': '17031'}