pyaugmecon 1.0.8

An AUGMECON based multi-objective optimization solver for Pyomo.

An AUGMECON based multi-objective optimization solver for Pyomo

This is a Python and Pyomo based implementation of the augmented ε-constraint (AUGMECON) method and its variants which can be used to solve multi-objective optimization problems, i.e., to return an approximation of the exact Pareto front.

It currently supports:

• Inner loop early exit (AUGMECON)
• Bypass coefficient (AUGMECON2)
• Flag array (AUGMECON-R)
• Processs-based parallelization

GAMS implementations of the method and its variants were provided by the authors of the papers. To the best of our knowledge, this is the first publicly available Python implementation.

Contents

• Installation
• Usage
• Tests
• Useful resources
• Notes
• Known limitations

Installation

PyAUGMECON can be installed from PyPI using pip install pyaugmecon. Detailed installation instructions for both Anaconda and PyPI installations are available below.

Requirements

• Python
• Pyomo
• Numpy
• Pandas
• Cloudpickle
• Pymoo
• Gurobi (other sovlers supported as well)

The Gurobi Python bindings are significantly faster than using the executable. So even if you have already installed Gurobi, please still install the Python version for an optimal experience.

Anaconda installation (advised)

This installation is advised as the PyPI installation of Gurobi does not include the licensing tools. Only Gurobi needs to be installed as other tools are by default included in Anaconda or will be automatically installed as dependencies of PyAUGMECON.

# Install Anaconda from https://www.anaconda.com

# Install Gurobi
conda config --add channels http://conda.anaconda.org/gurobi
conda install gurobi

# Install PyAUGMECON, and dependencies
pip install pyaugmecon

PyPI installation

For a PyPI installation, only Gurobi needs to be installed, other requirements will automatically be installed as dependencies of PyAUGMECON. As mentioned above, the PyPI version of Gurobi does not include licensing tools, see Gurobi documentation on how to install these.

# Install Gurobi, PYAUGMECON, and dependencies
pip install gurobipy pyaugmecon

Other solvers

The installation instructions above describe an installation with Gurobi, since this is the solver that has been used for testing. However, other solvers should work as well. The only requirement is that the solver is supported by Pyomo. An example of the GLPK solver can be found in the test folder.

Usage

First, an optimization model generator function needs to be created to pass the model to PyAUGMECON. This function should return an unsolved instance of the optimization problem, see Pyomo model. Essentially, the only difference in comparison with creating a single objective Pyomo optimization model is the fact that multiple objectives are defined using Pyomo's ObjectiveList().

The rest of the process is automatically being taken care of by instantiating a PyAugmecon() object and solving it afterward with PyAugmecon.solve().

PyAugmecon parameters

Instantiating a PyAugmecon(model, opts, solver_opts) object requires the following parameters:

Name	Description	Required
model	Function that returns an unsolved instance of the optimization problem, see Pyomo model	Yes
opts	Dictionary of PyAUGMECON related options, see PyAugmecon options	Yes
solver_opts	Dictionary of solver (Gurobi) options, see solver options	No

Example 3kp40

The following snippet shows how to solve the 3kp40 model from the tests folder:

from pyaugmecon import PyAugmecon
from tests.optimization_models import three_kp_model

# Multiprocessing requires this If statement (on Windows)
if __name__ == "__main__":
    model_type = "3kp40"

    # AUGMECON related options
    opts = {
        "name": model_type,
        "grid_points": 540,
        "nadir_points": [1031, 1069],
    }

    pyaugmecon = PyAugmecon(three_kp_model(model_type), opts)  # instantiate  PyAugmecon
    pyaugmecon.solve()  # solve PyAugmecon multi-objective optimization problem
    sols = pyaugmecon.get_pareto_solutions()  # get all pareto solutions
    payoff = pyaugmecon.get_payoff_table()  # get the payoff table
    decision_vars = pyaugmecon.get_decision_variables(sols[0])  # get the decision variables

PyAugmecon methods

Name	Description
get_pareto_solutions()	Get a list of Pareto-optimal solutions (list[tuple])
get_decision_variables(pareto_solution)	Get a dictionary of decision variables for a given Pareto-optimal solution. Where the key represents the decision variable name and the value is a pd.Series with the values.
get_payoff_table()	Get the payoff table from the model.

PyAugmecon attributes

After solving the model with PyAugmecon.solve(), the following object attributes are available:

Name	Description
sols	Full unprocessed solutions, only checked for uniqueness
unique_sols	Unique solutions, rounded to round_decimals and checked for uniqueness
unique_pareto_sols	Unique Pareto solutions, only dominant solutions, rounded to round_deicmals and checked for uniqueness
num_sols	Number of solutions
num_unique_sols	Number of unique solutions
num_unique_pareto_sols	Number of unique Pareto solutions
model.models_solved	Number of models solved (excluding solves for payoff table)
model.infeasibilites	Number of models solved that were infeasible
model.to_solve	Total number of models to solve (including payoff table)
model.model	Pyomo model instance
model.payoff	Payoff table
model.e	Gridpoints of p-1 objective functions that are used as constraints
hv_indicator	The hypervolume indicator of the unique Pareto solutions, see Pymoo documentation for more details
runtime	Total runtime in seconds

PyAugmecon solutions details

The solutions are stored in the sols, unique_sols, and unique_pareto_sols attributes. These differ in the following ways:

• sols contains all solutions, including duplicates and non-Pareto solutions
• unique_sols contains all unique solutions, including non-Pareto solutions
• unique_pareto_sols contains all unique Pareto solutions

The solutions are stored dictionary, where each dictionary contains the objective values (keys) and the decision variables (values). The objective values are stored as a tuple, where the first element is the objective value of the first objective function, the second element is the objective value of the second objective function, and so on. The decision variables are stored as a dictionary, where the keys are the names of the decision variables and the values are the values of the decision variables.

This is an example of the unique_pareto_sols attribute:

{
    (3, 5): {'x': pd.Series([1, 2]), 'y': pd.Series([2, 3])},
    (4, 4): {'x': pd.Series([3, 4]), 'y': pd.Series([1, 2])},
    (5, 3): {'x': pd.Series([2, 4]), 'y': pd.Series([1, 3])},
}

PyAugmecon options

The following PyAUGMECON related options can be passed as a dictionary to the solver:

Option	Description	Default
name	Name of the model, used in logging output	Undefined
grid_points	Number of grid points (if X grid points in GAMS, X+1 grid points are needed in order to exactly replicate results)
nadir_points	If the exact nadir points are known — these can be provided, otherwise they will be taken from the payoff table	True
early_exit	Use inner loop early-exit (AUGMECON) — exit the inner loop when the model result becomes infeasible	True
bypass_coefficient	Use bypass coefficient (AUGMECON2) — utilize slack variables to skip redundant iterations	True
flag_array	Use flag array (AUGMECON-R) — store early exit and bypass coefficients for upcoming iterations	True
penalty_weight	The penalty by which other terms in the objective are multiplied, usually between 10e-3 and 10e-6	1e-3
nadir_ratio	For problems with three or more objective functions, the payoff table minima do not guarantee the exact nadir. This factor scales the payoff minima. By providing a value lower than one, lower bounds of the minima can be estimated	1
logging_folder	Folder to store log files (relative to the root directory)	logs
cpu_count	Specify over how many processes the work should be divided. Use 1 to disable parallelization	multiprocessing.cpu_count()
redivide_work	Have processes take work from unfinished processes when their queue is empty	True
shared_flag	Share the flag array between processes. If false, each process will have its own flag array	True
round_decimals	The number of decimals to which the solutions are rounded before checking for uniqueness	9
pickle_file	Filename of the pickled Pyomo model for cloudpickle	model.p
solver_name	Name of the solver provided to Pyomo	gurobi
solver_io	Name of the solver interface provided to Pyomo. As the Gurobi Python bindings are significantly faster than using the executable, they are set by default. If you still prefer to use the executable, set this option to None	python
output_excel	Create an excel with solver output in the logging_folder containing the payoff table, e-points, solutions, unique solutions, and unique Pareto solutions	True
process_logging	Outputs solution process information from sub-processes to the log file. This significantly reduces performances and does not work on Windows. Don't enable this unless necessary	False
process_timeout	Gracefully stops all processes after process_timeout in seconds. As processes are stopped gracefully they will be allowed to finish their assigned work and not stop exactly after process_timeout	None

Solver options

To solve the single-objective optimization problems generated by the AUGMECON method they are passed to Gurobi, an external solver. All Gurobi solver parameters can be passed to Gurobi with this dictionary. Two parameters are already set by default, but can be overridden:

Option	Description	Default
MIPGap	See Gurobi documentation: MIPGap	0.0
NonConvex	See Gurobi documentation NonConvex	2

Pyomo model

A multi-objective Pyomo model is similar to a single-objective model, except that the objectives should be defined using Pyomo's ObjectiveList() and attached to the model by an attribute named obj_list. Multiple example models can be found in the tests folder and a simple model can be found below:

def two_objective_model():
    model = ConcreteModel()

    # Define variables
    model.x1 = Var(within=NonNegativeReals)
    model.x2 = Var(within=NonNegativeReals)

    # --------------------------------------
    #   Define the objective functions
    # --------------------------------------

    def objective1(model):
        return model.x1

    def objective2(model):
        return 3 * model.x1 + 4 * model.x2

    # --------------------------------------
    #   Define the regular constraints
    # --------------------------------------

    def constraint1(model):
        return model.x1 <= 20

    def constraint2(model):
        return model.x2 <= 40

    def constraint3(model):
        return 5 * model.x1 + 4 * model.x2 <= 200

    # --------------------------------------
    #   Add components to the model
    # --------------------------------------

    # Add the constraints to the model
    model.con1 = Constraint(rule=constraint1)
    model.con2 = Constraint(rule=constraint2)
    model.con3 = Constraint(rule=constraint3)

    # Add the objective functions to the model using ObjectiveList(). Note
    # that the first index is 1 instead of 0!
    model.obj_list = ObjectiveList()
    model.obj_list.add(expr=objective1(model), sense=maximize)
    model.obj_list.add(expr=objective2(model), sense=maximize)

    # By default deactivate all the objective functions
    for o in range(len(model.obj_list)):
        model.obj_list[o + 1].deactivate()

    return model

Tests

To test the correctness of the proposed approach, the following optimization models have been tested both using the GAMS and the proposed implementations:

• Sample bi-objective problem of the original paper
• Sample three-objective problem of the original implementation/presentation
• Multi-objective multidimensional knapsack (MOMKP) problems: 2kp50, 2kp100, 2kp250, 3kp40, 3kp50, 4kp40 and 4kp50
• Multi-objective economic dispatch (minimize cost, emissions and unmet demand)

These models can be found in the tests folder and run with PyTest to test for the correctness of the payoff table and the Pareto solutions. Example: pytest tests/test_3kp40.py

Despite the extensive testing, the PyAUGMECON implementation is provided without any warranty. Use it at your own risk!

Useful resources

Original AUGMECON

• G. Mavrotas, “Effective implementation of the ε-constraint methodin Multi-Objective Mathematical Programming problems,” Applied Mathematics and Computation, vol. 213, no. 2, pp. 455–465, 2009
• GAMS implementation

Improved AUGMECON2

• Mavrotas and K. Florios, “An improved version of the augmented ε-constraint method (AUGMECON2) for finding the exact pareto set in multi-objective integer programming problems,” Applied Mathematics and Computation, vol. 219, no. 18, pp. 9652–9669, 2013.
• GAMS implementation

Further improved AUGMECON-R

• A. Nikas, A. Fountoulakis, A. Forouli, and H. Doukas, A robust augmented ε-constraint method (AUGMECON-R) for finding exact solutions of multi-objective linear programming problems. Springer Berlin Heidelberg, 2020, no. 0123456789.
• GAMS implementation

Other resources

• The following PhD thesis is also very useful
• This presentation provides a quick overview

Notes

• The choice of the objective function(s) to add as constraints affects the mapping of the Pareto front. Empirically, having one of the objective functions with a large range in the constraint set tends to result in a denser representation of the Pareto front. Nevertheless, despite the choice of which objectives to add in the constraint set, the solutions belong to the same Pareto front (obviously).

• If X grid points in GAMS, X+1 grid points are needed in PyAUGMECON in order to exactly replicate results

• For some models, it is beneficial to runtime to reduce the number of solver (Gurobi) threads. More details: Gurobi documentation

Known limitations

• For relatively small models (most of the knapsack problems), disabling redivide_work and shared_flag should lead to slightly lower runtimes. This is due to the overhead of sharing the flag array between processes that don't have a shared memory space. Redividing the work results in additional models solved with duplicate solutions, leading to longer runtime. By sharing the flag array in a different way and dividing the work more smartly these limitations might be solved in the future.

• To parallelize the solution process, the grid is divided into blocks with a minimum length of the inner loop (number of grid points) over the processes. With fewer than three objectives, there is only one dimension to iterate over, and as such no parallelization possible. This is something that can be improved in the future.

Changelog

See Changelog

Citing

If you use PyAUGMECON for academic work, please cite the following DOI (Zenodo).

Credit

This software was developed at the Electricity Markets & Power System Optimization Laboratory (EMPSOLab), Electrical Energy Systems Group, Department of Electrical Engineering, Eindhoven University of Technology.

Contributors:

• Wouter Bles (current version of the package)
• Nikolaos Paterakis (initial implementation)