Multi-objective optimization example Pyomo

Question

Any example for multi-objective optimization in Pyomo?

I am trying to minimize 4 Objectives (Non Linear) and I would like to use pyomo and ipopt. Have also access to Gurobi.

I want to see even very simple example where we try to optimize for two or more objective (one minimization and one maximization) for a list of decision variables (not just one dimension but maybe a vector).

Pyomo book that I have (https://link.springer.com/content/pdf/10.1007%2F978-3-319-58821-6.pdf) does not provide a signle clue.

Answer 1

With Pyomo you have to implement it yourself. I am doing it right now. The best method is the augmented epsilon-constraint method. It will always be efficient and always find the global pareto-optimum. Best example is here: Effective implementation of the epsilon-constraint method in Multi-Objective Mathematical Programming problems, Mavrotas, G, 2009.

Edit: Here I programmed the example from the Paper above in pyomo: It will first maximize for f1 then for f2. Then It'll apply the normal epsilon-constraint and plot the inefficient Pareto-front and then It'll apply the augmented epsilon-constraint, which finally is the method to go with!

from pyomo.environ import *
import matplotlib.pyplot as plt


# max f1 = X1 <br>
# max f2 = 3 X1 + 4 X2 <br>
# st  X1 <= 20 <br>
#     X2 <= 40 <br>
#     5 X1 + 4 X2 <= 200 <br>

model = ConcreteModel()

model.X1 = Var(within=NonNegativeReals)
model.X2 = Var(within=NonNegativeReals)

model.C1 = Constraint(expr = model.X1 <= 20)
model.C2 = Constraint(expr = model.X2 <= 40)
model.C3 = Constraint(expr = 5 * model.X1 + 4 * model.X2 <= 200)

model.f1 = Var()
model.f2 = Var()
model.C_f1 = Constraint(expr= model.f1 == model.X1)
model.C_f2 = Constraint(expr= model.f2 == 3 * model.X1 + 4 * model.X2)
model.O_f1 = Objective(expr= model.f1  , sense=maximize)
model.O_f2 = Objective(expr= model.f2  , sense=maximize)

model.O_f2.deactivate()

solver = SolverFactory('cplex')
solver.solve(model);

print( '( X1 , X2 ) = ( ' + str(value(model.X1)) + ' , ' + str(value(model.X2)) + ' )')
print( 'f1 = ' + str(value(model.f1)) )
print( 'f2 = ' + str(value(model.f2)) )
f2_min = value(model.f2)


# ## max f2

model.O_f2.activate()
model.O_f1.deactivate()

solver = SolverFactory('cplex')
solver.solve(model);

print( '( X1 , X2 ) = ( ' + str(value(model.X1)) + ' , ' + str(value(model.X2)) + ' )')
print( 'f1 = ' + str(value(model.f1)) )
print( 'f2 = ' + str(value(model.f2)) )
f2_max = value(model.f2)


# ## apply normal $\epsilon$-Constraint

model.O_f1.activate()
model.O_f2.deactivate()

model.e = Param(initialize=0, mutable=True)

model.C_epsilon = Constraint(expr = model.f2 == model.e)

solver.solve(model);

print('Each iteration will keep f2 lower than some values between f2_min and f2_max, so ['       + str(f2_min) + ', ' + str(f2_max) + ']')

n = 4
step = int((f2_max - f2_min) / n)
steps = list(range(int(f2_min),int(f2_max),step)) + [f2_max]

x1_l = []
x2_l = []
for i in steps:
    model.e = i
    solver.solve(model);
    x1_l.append(value(model.X1))
    x2_l.append(value(model.X2))
plt.plot(x1_l,x2_l,'o-.');
plt.title('inefficient Pareto-front');
plt.grid(True);


# ## apply augmented $\epsilon$-Constraint

# max   f2 + delta*epsilon <br>
#  s.t. f2 - s = e

model.del_component(model.O_f1)
model.del_component(model.O_f2)
model.del_component(model.C_epsilon)

model.delta = Param(initialize=0.00001)

model.s = Var(within=NonNegativeReals)

model.O_f1 = Objective(expr = model.f1 + model.delta * model.s, sense=maximize)

model.C_e = Constraint(expr = model.f2 - model.s == model.e)

x1_l = []
x2_l = []
for i in range(160,190,6):
    model.e = i
    solver.solve(model);
    x1_l.append(value(model.X1))
    x2_l.append(value(model.X2))
plt.plot(x1_l,x2_l,'o-.');
plt.title('efficient Pareto-front');
plt.grid(True);

Answer 2

To my knowledge, while Pyomo supports the expression of models with multiple objectives, it does not yet have automatic model transformations to generate common multi-objective optimization formulations for you.

That said, you can still create these formulations yourself. Take a look at epsilon-constraint, 1-norm, and infinity norm for some ideas.

Answer 3

Disclaimer: I am the main developer of pymoo, a multi-objective optimization framework in Python.

You might want to consider other frameworks in Python that have a focus on multi-objective optimization. For instance, in pymoo the definition of the rather simple test problem mentioned above is more or less straightforward. You can find an implementation of it below. The results in the design and objectives space look as follows:

pymoo is well documented and provides a getting started guide that demonstrates defining your own optimization problem, obtaining a set of near-optimal solutions and analyzing it: https://pymoo.org/getting_started.html

The focus of the framework is anything related to multi-objective optimization including visualization and decision making.

import matplotlib.pyplot as plt
import numpy as np

from pymoo.algorithms.nsga2 import NSGA2
from pymoo.model.problem import Problem
from pymoo.optimize import minimize
from pymoo.visualization.scatter import Scatter


class MyProblem(Problem):
    def __init__(self):
        """
        max f1 = X1 <br>
        max f2 = 3 X1 + 4 X2 <br>
        st  X1 <= 20 <br>
            X2 <= 40 <br>
            5 X1 + 4 X2 <= 200 <br>
        """

        super().__init__(n_var=2,
                         n_obj=2,
                         n_constr=1,
                         xl=np.array([0, 0]),
                         xu=np.array([20, 40]))

    def _evaluate(self, x, out, *args, **kwargs):
        # define both objectives
        f1 = x[:, 0]
        f2 = 3 * x[:, 0] + 4 * x[:, 1]

        # we have to negate the objectives because by default we assume minimization
        f1, f2 = -f1, -f2

        # define the constraint as a less or equal to zero constraint
        g1 = 5 * x[:, 0] + 4 * x[:, 1] - 200

        out["F"] = np.column_stack([f1, f2])
        out["G"] = g1


problem = MyProblem()

algorithm = NSGA2()

res = minimize(problem,
               algorithm,
               ('n_gen', 200),
               seed=1,
               verbose=True)

print(res.X)
print(res.F)

fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(12, 6))
Scatter(fig=fig, ax=ax1, title="Design Space").add(res.X, color="blue").do()
Scatter(fig=fig, ax=ax2, title="Objective Space").add(res.F, color="red").do()
plt.show()