Model

PyEPO is an implementation that aims to support an end-to-end predict-then-optimize with linear objective function and unknown cost coefficients. The core component is the differentiable optimization solver, which is involved in updating the gradient of the cost coefficients with respect to the optimal solution.

optModel is a module of PyEPO library. It is not a solver but serves as a container of a solver or an algorithm. This design allows for flexibility in the selection of solvers and algorithms by users. Other modules of PyEPO can use optModel for tasks such as training and testing.

PyEPO contains several pre-defined optimization models with GurobiPy and Pyomo. It includes the shortest path problem (GurobiPy & Pyomo), the knapsack problem (GurobiPy & Pyomo), and the traveling salesman problem (GurobiPy).

To build optimization models with PyEPO, users do not need specific costs of objective functions since the cost vector is unknown but will be estimated from data.

For more information and details about the Optimization Model, please see the 01 Optimization Model:

User-defined Models

User can build optimization problem with linear objective function. Our API is also designed to support users to define their own problems based on GurobiPy and Pyomo. Besides the API of GurobiPy & Pyomo, users can also build problems from scratch with whatever solvers and algorithms they want to use.

optModel treats these solvers as black boxes and provides interfaces _getModel, setObj, and solve.

_getModel: Build and return optimization solver and corresponding decision variables.
setObj: Give a cost vector to set the objective function.
solve: Solve optimization problem and return optimal solution and objective value.

User-defined GurobiPy Models

User-defined models with GurobiPy can be easily defined by the inheritance of the abstract class pyepo.model.grb.optGrbModel.

For optGrbModel, users does not need specify the sense modelSense as EPO.MINIMIZE or EPO.MINIMIZE. Both the minimization and maximization models are correctly recognized and run by pyepo.

class pyepo.model.grb.optGrbModel

This is an abstract class for Gurobi-based optimization model

_model

Gurobi model

Type:: GurobiPy model

__init__()

abstract _getModel()

An abstract method to build a model from a optimization solver

Returns:: optimization model and variables
Return type:: tuple

property num_cost: number of cost to be predicted

relax(): A unimplemented method to relax MIP model

setObj(c)

A method to set objective function

Parameters:: c (np.ndarray / list) – cost of objective function

solve()

A method to solve model

Returns:: optimal solution (list) and objective value (float)
Return type:: tuple

For example, users can build models for the following problem:

\[\begin{split}\begin{aligned} \max_{x} & \sum_{i=0}^4 c_i x_i \\ s.t. \quad & 3 x_0 + 4 x_1 + 3 x_2 + 6 x_3 + 4 x_4 \leq 12 \\ & 4 x_0 + 5 x_1 + 2 x_2 + 3 x_3 + 5 x_4 \leq 10 \\ & 5 x_0 + 4 x_1 + 6 x_2 + 2 x_3 + 3 x_4 \leq 15 \\ & \forall x_i \in \{0, 1\} \end{aligned}\end{split}\]

In the general case, users only need to implement _getModel method with GurobiPy.

import random

import gurobipy as gp
from gurobipy import GRB

from pyepo.model.grb import optGrbModel

class myModel(optGrbModel):

    def _getModel(self):
        # create a model
        m = gp.Model()
        # variables
        x = m.addVars(5, name="x", vtype=GRB.BINARY)
        # model sense
        m.modelSense = GRB.MAXIMIZE
        # constraints
        m.addConstr(3 * x[0] + 4 * x[1] + 3 * x[2] + 6 * x[3] + 4 * x[4] <= 12)
        m.addConstr(4 * x[0] + 5 * x[1] + 2 * x[2] + 3 * x[3] + 5 * x[4] <= 10)
        m.addConstr(5 * x[0] + 4 * x[1] + 6 * x[2] + 2 * x[3] + 3 * x[4] <= 15)
        return m, x

myoptmodel = myModel()
cost = [random.random() for _ in range(myoptmodel.num_cost)] # random cost vector
myoptmodel.setObj(cost) # set objective function
myoptmodel.solve() # solve

User-defined Pyomo Models

User-defined models with Pyomo can be easily defined by the inheritance of the abstract class pyepo.model.omo.optOmoModel.

class pyepo.model.omo.optOmoModel(solver='glpk')

This is an abstract class for Pyomo-based optimization model

_model

Pyomo model

Type:: PyOmo model

solver

optimization solver in the background

Type:: str

Parameters:: solver (str) – optimization solver in the background

__init__(solver='glpk')

Parameters:: solver (str) – optimization solver in the background

abstract _getModel()

An abstract method to build a model from a optimization solver

Returns:: optimization model and variables
Return type:: tuple

property num_cost: number of cost to be predicted

relax(): A unimplemented method to relax MIP model

setObj(c)

A method to set objective function

Parameters:: c (np.ndarray / list) – cost of objective function

solve()

A method to solve model

Returns:: optimal solution (list) and objective value (float)
Return type:: tuple

Let’s build models for the problem again with Pyomo:

\[\begin{split}\begin{aligned} \max_{x} & \sum_{i=0}^4 c_i x_i \\ s.t. \quad & 3 x_0 + 4 x_1 + 3 x_2 + 6 x_3 + 4 x_4 \leq 12 \\ & 4 x_0 + 5 x_1 + 2 x_2 + 3 x_3 + 5 x_4 \leq 10 \\ & 5 x_0 + 4 x_1 + 6 x_2 + 2 x_3 + 3 x_4 \leq 15 \\ & \forall x_i \in \{0, 1\} \end{aligned}\end{split}\]

In the general case, users only need to implement _getModel method with Pyomo.

Warning

Unlike the optGrbModel, the optOmoModel need to set modelSense in the _getModel.

import random

from pyomo import environ as pe

from pyepo.model.omo import optOmoModel
from pyepo import EPO

class myModel(optOmoModel):

    def _getModel(self):
        # sense
        self.modelSense = EPO.MAXIMIZE
        # create a model
        m = pe.ConcreteModel()
        # variables
        x = pe.Var([0,1,2,3,4], domain=pe.Binary)
        m.x = x
        # constraints
        m.cons = pe.ConstraintList()
        m.cons.add(3 * x[0] + 4 * x[1] + 3 * x[2] + 6 * x[3] + 4 * x[4] <= 12)
        m.cons.add(4 * x[0] + 5 * x[1] + 2 * x[2] + 3 * x[3] + 5 * x[4] <= 10)
        m.cons.add(5 * x[0] + 4 * x[1] + 6 * x[2] + 2 * x[3] + 3 * x[4] <= 15)
        return m, x

myoptmodel = myModel(solver="gurobi")
cost = [random.random() for _ in range(myoptmodel.num_cost)] # random cost vector
myoptmodel.setObj(cost) # set objective function
myoptmodel.solve() # solve

User-defined Models from Scratch

pyepo.model.opt.optModel provides an abstract class for users to create an optimization model with any solvers or algorithms. By overriding _getModel, setObj, solve, and num_cost, user-defined optModel can work for end-to-end training.

class pyepo.model.opt.optModel

This is an abstract class for optimization model

_model

Gurobi model

Type:: GurobiPy model

__init__()

abstract _getModel()

An abstract method to build a model from a optimization solver

Returns:: optimization model and variables
Return type:: tuple

property num_cost: number of cost to be predicted

abstract setObj(c)

An abstract method to set objective function

Parameters:: c (ndarray) – cost of objective function

abstract solve()

An abstract method to solve model

Returns:: optimal solution (list) and objective value (float)
Return type:: tuple

Warning

The optModel need to set modelSense in the _getModel. If not set, the default is to minimize.

For example, we can use networkx to solve the previous shortest path problem using the Dijkstra algorithm. And pyepo.model.opt.optModel allows users to create a model in the following way:

import random

import numpy as np
import networkx as nx

from pyepo.model.opt import optModel

class myShortestPathModel(optModel):

    def __init__(self, grid):
        """
        Args:
            grid (tuple): size of grid network
        """
        self.grid = grid
        self.arcs = self._getArcs()
        super().__init__()

    def _getArcs(self):
        """
        A method to get list of arcs for grid network

        Returns:
            list: arcs
        """
        arcs = []
        for i in range(self.grid[0]):
            # edges on rows
            for j in range(self.grid[1] - 1):
                v = i * self.grid[1] + j
                arcs.append((v, v + 1))
            # edges in columns
            if i == self.grid[0] - 1:
                continue
            for j in range(self.grid[1]):
                v = i * self.grid[1] + j
                arcs.append((v, v + self.grid[1]))
        return arcs

    def _getModel(self):
        """
        A method to build model

        Returns:
            tuple: optimization model and variables
        """
        # build graph as optimization model
        g = nx.Graph()
        # add arcs as variables
        g.add_edges_from(self.arcs, cost=0)
        return g, g.edges

    def setObj(self, c):
        """
        A method to set objective function

        Args:
            c (ndarray): cost of objective function
        """
        for i, e in enumerate(self.arcs):
            self._model.edges[e]["cost"] = c[i]

    def solve(self):
        """
        A method to solve model

        Returns:
            tuple: optimal solution (list) and objective value (float)
        """
        # dijkstra
        path = nx.shortest_path(self._model, weight="cost", source=0, target=self.grid[0]*self.grid[1]-1)
        # convert path into active edges
        edges = []
        u = 0
        for v in path[1:]:
            edges.append((u,v))
            u = v
        # init sol & obj
        sol = np.zeros(self.num_cost)
        obj = 0
        # convert active edges into solution and obj
        for i, e in enumerate(self.arcs):
            if e in edges:
                sol[i] = 1 # active edge
                obj += self._model.edges[e]["cost"] # cost of active edge
        return sol, obj

# solve model
grid = (5,5)
myoptmodel = myShortestPathModel(grid)
cost = [random.random() for _ in range(myoptmodel.num_cost)] # random cost vector
myoptmodel.setObj(cost) # set objective function
sol, obj = myoptmodel.solve() # solve
# print res
print('Obj: {}'.format(obj))
for i, e in enumerate(myoptmodel.arcs):
    if sol[i] > 1e-3:
        print(e)

Pre-defined Models

Pre-defined models are classic optimization problems, including shortest path, multidimensional knapsack, and traveling salesman.

Shortest Path

It is a (h,w) grid network and the goal is to find the shortest path from northwest to southeast. In our examples, the grid size of network is (5,5).

In PyEPO, the shortest path problem is built as a linear program and formulated as a minimum cost flow problem.

Shortest Path GurobiPy Model

The optModel is built from pyepo.model.grb.shortestPathModel, in which API uses GurobiPy to model the shortest path problem.

class pyepo.model.grb.shortestPathModel(grid)

This class is optimization model for shortest path problem

_model

Gurobi model

Type:: GurobiPy model

grid

Size of grid network

Type:: tuple of int

arcs

List of arcs

Type:: list

Parameters:: grid (tuple of int) – size of grid network

__init__(grid)

Parameters:: grid (tuple of int) – size of grid network

property num_cost: number of cost to be predicted

setObj(c)

A method to set objective function

Parameters:: c (np.ndarray / list) – cost of objective function

solve()

A method to solve model

Returns:: optimal solution (list) and objective value (float)
Return type:: tuple

import pyepo

grid = (5,5) # network grid
optmodel = pyepo.model.grb.shortestPathModel(grid) # build model

Users can use setObj to assign a specific cost vector and use solve to optimize. However, setObj or solve methods do not require manual calls during training.

import random
cost = [random.random() for _ in range(optmodel.num_cost)] # random cost vector
optmodel.setObj(cost) # set objective function
optmodel.solve() # solve

Shortest Path Pyomo Model

The optModel is built from pyepo.model.omo.shortestPathModel, in which API uses Pyomo to model the shortest path problem.

class pyepo.model.omo.shortestPathModel(grid, solver='glpk')

This class is optimization model for shortest path problem

_model

Pyomo model

Type:: PyOmo model

solver

optimization solver in the background

Type:: str

grid

size of grid network

Type:: tuple of int

arcs

list of arcs

Type:: list

Parameters:

grid (tuple of int) – size of grid network
solver (str) – optimization solver in the background

__init__(grid, solver='glpk')

Parameters:

grid (tuple of int) – size of grid network
solver (str) – optimization solver in the background

property num_cost: number of cost to be predicted

setObj(c)

A method to set objective function

Parameters:: c (np.ndarray / list) – cost of objective function

solve()

A method to solve model

Returns:: optimal solution (list) and objective value (float)
Return type:: tuple

Pyomo supports a wide variety of solvers in the background (e.g. BARON, CBC, CPLEX, and Gurobi). Thus, pyepo.model.omo.shortestPathModel allows users to call different solvers with class parameter solver.

import pyepo

grid = (5,5) # network grid
optmodel = pyepo.model.omo.shortestPathModel(grid, solver="glpk") # build model with glpk
optmodel = pyepo.model.omo.shortestPathModel(grid, solver="gurobi") # build model with gurobi

You can get the current list of supported solvers using the pyomo command:

pyomo help --solvers

Same as pyepo.model.grb.shortestPathModel, methods setObj and solve can specify objective function and solve the problem.

import random
cost = [random.random() for _ in range(optmodel.num_cost)] # random cost vector
optmodel.setObj(cost) # set objective function
optmodel.solve() # solve

Knapsack

Multi-dimensional knapsack problem is a maximization problem with multiple resource constraints: Given a set of items, the aim is to find a collection that the total weights in is less than or equal to resource capacities and the total value is as large as possible. Let’s define a 3D knapsack problem as follow:

\[\begin{split}\begin{aligned} \max_{x} & \sum_{i=0}^4 c_i x_i \\ s.t. \quad & 3 x_0 + 4 x_1 + 3 x_2 + 6 x_3 + 4 x_4 \leq 12 \\ & 4 x_0 + 5 x_1 + 2 x_2 + 3 x_3 + 5 x_4 \leq 10 \\ & 5 x_0 + 4 x_1 + 6 x_2 + 2 x_3 + 3 x_4 \leq 15 \\ & \forall x_i \in \{0, 1\} \end{aligned}\end{split}\]

Constraints coefficients weights and constraints rhs capacities are required to create models.

Note

The dimension of the knapsack and the number of items are implicitly defined by the shape of weights and capacities.

Knapsack GurobiPy Model

The optModel is built from pyepo.model.grb.knapsackModel, in which API uses GurobiPy to model the knapsack problem.

class pyepo.model.grb.knapsackModel(weights, capacity)

This class is optimization model for knapsack problem

_model

Gurobi model

Type:: GurobiPy model

weights

Weights of items

Type:: np.ndarray / list

capacity

Total capacity

Type:: np.ndarray / listy

items

List of item index

Type:: list

Parameters:

weights (np.ndarray / list) – weights of items
capacity (np.ndarray / list) – total capacity

__init__(weights, capacity)

Parameters:

weights (np.ndarray / list) – weights of items
capacity (np.ndarray / list) – total capacity

property num_cost: number of cost to be predicted

relax(): A method to get linear relaxation model

setObj(c)

A method to set objective function

Parameters:: c (np.ndarray / list) – cost of objective function

solve()

A method to solve model

Returns:: optimal solution (list) and objective value (float)
Return type:: tuple

import pyepo

weights = [[3, 4, 3, 6, 4],
           [4, 5, 2, 3, 5],
           [5, 4, 6, 2, 3]] # constraints coefficients
capacities = [12, 10, 15] # constraints rhs
optmodel = pyepo.model.grb.knapsackModel(weights, capacities) # build model

Similarly, users can use setObj to assign a specific cost vector and use solve to optimize. However, setObj or solve methods do not require manual calls during training.

import random
cost = [random.random() for _ in range(optmodel.num_cost)] # random cost vector
optmodel.setObj(cost) # set objective function
optmodel.solve() # solve

In mathematics, the relaxation of a (Mixed) Integer Linear Programming is the problem that arises by removing the integrality constraint of each variable. As an ILP, optGrbModel allows users to relax ILP with relax method to obtain a relaxation optModel from the original.

optmodel_rel = optmodel.relax() # relax

Knapsack Pyomo Model

The optModel is built from pyepo.model.omo.knapsackModel, in which API uses Pyomo to model the knapsack problem.

class pyepo.model.omo.knapsackModel(weights, capacity, solver='glpk')

This class is optimization model for knapsack problem

_model

Pyomo model

Type:: PyOmo model

solver

optimization solver in the background

Type:: str

weights

weights of items

Type:: np.ndarray

capacity

total capacity

Type:: np.ndarray

items

list of item index

Type:: list

Parameters:

weights (np.ndarray / list) – weights of items
capacity (np.ndarray / list) – total capacity
solver (str) – optimization solver in the background

__init__(weights, capacity, solver='glpk')

Parameters:

weights (np.ndarray / list) – weights of items
capacity (np.ndarray / list) – total capacity
solver (str) – optimization solver in the background

property num_cost: number of cost to be predicted

relax(): A method to get linear relaxation model

setObj(c)

A method to set objective function

Parameters:: c (np.ndarray / list) – cost of objective function

solve()

A method to solve model

Returns:: optimal solution (list) and objective value (float)
Return type:: tuple

import pyepo

weights = [[3, 4, 3, 6, 4],
           [4, 5, 2, 3, 5],
           [5, 4, 6, 2, 3]] # constraints coefficients
capacities = [12, 10, 15] # constraints rhs
# build model with GLPK
optmodel = pyepo.model.omo.knapsackModel(weights, capacities, solver="glpk")
# build model with Gurobi
optmodel = pyepo.model.omo.knapsackModel(weights, capacities, solver="gurobi")

Same as pyepo.model.grb.knapsackModel, users can use setObj, solve, and relax methods.

You can get the current list of supported solvers using the pyomo command:

pyomo help --solvers

Traveling Salesman

The traveling salesman problem (TSP) is the shortest route that visits each city exactly once and returns to the origin city. We consider the symmetric TSP, in which the distance between two cities is the same in each opposite direction. In our examples, the number of nodes is 20.

The TSP can be formulated as an Integer Linear Programming with several formulations. We implemented Dantzig–Fulkerson–Johnson (DFJ) formulation, Gavish–Graves (GG) formulation, and Miller–Tucker–Zemlin (MTZ) formulation.

Note

The implementation of TSP is only based on GurobiPy. Pyomo is not supported.

DFJ formulation

class pyepo.model.grb.tspDFJModel(num_nodes)

This class is optimization model for traveling salesman problem based on Danzig–Fulkerson–Johnson (DFJ) formulation and constraint generation.

_model

Gurobi model

Type:: GurobiPy model

num_nodes

Number of nodes

Type:: int

edges

List of edge index

Type:: list

Parameters:: num_nodes (int) – number of nodes

__init__(num_nodes)

Parameters:: num_nodes (int) – number of nodes

property num_cost: number of cost to be predicted

setObj(c)

A method to set objective function

Parameters:: c (list) – cost vector

solve(): A method to solve model

The number of subtour elimination constraints for DFJ formulation is exponential. Thus, we solved it with column generation. Because of that, the linear relaxation of DFJ is not supported in our implementation.

Same as previous model, the code for traveling salesman problem with DFJ formulation is as follows:

import pyepo
import random

num_nodes = 20 # number of nodes
optmodel = pyepo.model.grb.tspDFJModel(num_nodes) # build model

cost = [random.random() for _ in range(optmodel.num_cost)] # random cost vector
optmodel.setObj(cost) # set objective function
optmodel.solve() # solve

GG formulation

class pyepo.model.grb.tspGGModel(num_nodes)

This class is optimization model for traveling salesman problem based on Gavish–Graves (GG) formulation.

_model

Gurobi model

Type:: GurobiPy model

num_nodes

Number of nodes

Type:: int

edges

List of edge index

Type:: list

Parameters:: num_nodes (int) – number of nodes

__init__(num_nodes)

Parameters:: num_nodes (int) – number of nodes

property num_cost: number of cost to be predicted

relax(): A method to get linear relaxation model

setObj(c)

A method to set objective function

Parameters:: c (list) – cost vector

solve(): A method to solve model

Same as previous model, the code for traveling salesman problem with GG formulation is as follows:

import pyepo
import random

num_nodes = 20 # number of nodes
optmodel = pyepo.model.grb.tspGGModel(num_nodes) # build model

cost = [random.random() for _ in range(optmodel.num_cost)] # random cost vector
optmodel.setObj(cost) # set objective function
optmodel.solve() # solve

optmodel.relax() # relax

MTZ formulation

class pyepo.model.grb.tspMTZModel(num_nodes)

This class is optimization model for traveling salesman problem based on Miller-Tucker-Zemlin (MTZ) formulation.

_model

Gurobi model

Type:: GurobiPy model

num_nodes

Number of nodes

Type:: int

edges

List of edge index

Type:: list

Parameters:: num_nodes (int) – number of nodes

__init__(num_nodes)

Parameters:: num_nodes (int) – number of nodes

property num_cost: number of cost to be predicted

relax(): A method to get linear relaxation model

setObj(c)

A method to set objective function

Parameters:: c (list) – cost vector

solve(): A method to solve model

Same as previous model, the code for traveling salesman problem with MTZ formulation is as follows:

import pyepo
import random

num_nodes = 20 # number of nodes
optmodel = pyepo.model.grb.tspMTZModel(num_nodes) # build model

cost = [random.random() for _ in range(optmodel.num_cost)] # random cost vector
optmodel.setObj(cost) # set objective function
optmodel.solve() # solve

optmodel.relax() # relax

Portfolio

Portfolio optimization is a crucial concept in financial management and investment theory, primarily concerned with the process of selecting the best portfolio (asset distribution) out of the set of all portfolios being considered according to some objective. Our model seeks highest expected return for a given level of risk.

\[\begin{split}\begin{aligned} \max_{x} & \sum_{i} r_i x_i \\ s.t. \quad & \sum_{i} x_i = 1 \\ & \mathbf{x}^{\intercal} \mathbf{\Sigma} \mathbf{x} \leq \gamma \bar{\Sigma}\\ & \forall x_i \geq 0 \end{aligned}\end{split}\]

Portfolio GurobiPy Model

The optModel is built from pyepo.model.grb.portfolioModel, in which API uses GurobiPy to model the portfolio optimization.

class pyepo.model.grb.portfolioModel(num_assets, covariance, gamma=2.25)

This class is an optimization model for portfolio problem

_model

Gurobi model

Type:: GurobiPy model

num_assets

number of assets

Type:: int

covariance

covariance matrix of the returns

Type:: numpy.ndarray

risk_level

risk level

Type:: float

Parameters:

num_assets (int) – number of assets
covariance (numpy.ndarray) – covariance matrix of the returns
gamma (float) – risk level parameter

__init__(num_assets, covariance, gamma=2.25)

Parameters:

num_assets (int) – number of assets
covariance (numpy.ndarray) – covariance matrix of the returns
gamma (float) – risk level parameter

property num_cost: number of cost to be predicted

setObj(c)

A method to set objective function

Parameters:: c (np.ndarray / list) – cost of objective function

solve()

A method to solve model

Returns:: optimal solution (list) and objective value (float)
Return type:: tuple

import pyepo
m = 50 # number of assets
cov = np.cov(np.random.randn(10, m), rowvar=False) # covariance matrix
optmodel = pyepo.model.grb.portfolioModel(m, cov) # build model

Similarly, users can use setObj to assign a specific cost vector and use solve to optimize. However, setObj or solve methods do not require manual calls during training.

import random
revenue = [random.random() for _ in range(optmodel.num_cost)] # random cost vector
optmodel.setObj(revenue) # set objective function
optmodel.solve() # solve