ROmodel: Modeling robust optimization problems in Pyomo

This paper introduces ROmodel, an open source Python package extending the modeling capabilities of the algebraic modeling language Pyomo to robust optimization problems. ROmodel helps practitioners transition from deterministic to robust optimization through modeling objects which allow formulating robust models in close analogy to their mathematical formulation. ROmodel contains a library of commonly used uncertainty sets which can be generated using their matrix representations, but it also allows users to define custom uncertainty sets using Pyomo constraints. ROmodel supports adjustable variables via linear decision rules. The resulting models can be solved using ROmodels solvers which implement both the robust reformulation and cutting plane approach. ROmodel is a platform to implement and compare custom uncertainty sets and reformulations. We demonstrate ROmodel's capabilities by applying it to six case studies. We implement custom uncertainty sets based on (warped) Gaussian processes to show how ROmodel can integrate data-driven models with optimization.

model also allows access to a much larger number of deterministic solvers than these other solvers.
In contrast to AIMMS, which has some capabilities for modeling robust optimization problems [1], ROModel is open source and allows possibilities for extension. JumPeR [8], which extends Julia's modeling language JumP to robust optimization problems, is the most similar to ROModel. In contrast to JumPeR, ROModel is more tightly tied to its respective algebraic modeling language, e.g. developing new convex uncertainty sets in ROmodel only requires adding Pyomo constraints while JumPeR would require adding Julia functionality. ROmodel also automatically recognizes uncertainy set geometries and applies the corresponding transformations without the user having to specify which geometry they are using.
A further advantage of ROmodel is that it is based on Python, which is popular in data analytics and machine learning. ROmodel therefore allows data-based techniques to be integrated seemlesly with robust optimization methods. As an example, we implement (warped) Gaussian process-based uncertainty sets in ROmodel [26]. These sets seemlesly integrate Gaussian processes trained in the Python library GPy [10] into robust optimization problems. ROmodel is open source and available on Github [28].
The rest of this paper is structured as follows. Section 2 introduces ROmodel's new modeling objects and shows how they can be used to model robust optimization problems. Section 3 introduces ROmodel's three solvers: a reformulation based solver, a cutting plane solver, and a nominal solver for obtaining nominal solutions of robust problems. Section 4 discusses how ROmodel can be extended and demonstrates our implementation of Gaussian process-based uncertainty sets for black-box constrained problems. Section 5 introduces six case studies and presents results.

Modeling
Consider the following generic deterministic optimization problem min x∈X ,y f (x, y,ξ) (1a) where x ∈ R n and y ∈ R m are decision variables andξ ∈ R k is a vector of (nominal) parameters.
If the parameter vectorξ is not known exactly, we can construct the following robust version of Here x ∈ R n are "here and now" variables, determined before the uncertainty is revealed, while y(ξ) are adjustable variables, determined after the uncertainty is revealed. The uncertain parameter vector ξ is bounded by the uncertainty set U(x), which may depend on x and which contains the nominal valuesξ from Problem 1. Optimization Problem 2, a generic robust optimization problem with uncertainty in the objective and constraints, is the basis for ROmodel. Note that we are limiting ourselves here to one uncertain constraint for simplicity of notation only, ROmodel can handle multiple robust constraints. Also note that there can be an arbitrary number of deterministic constraints which define the set X .
ROmodel introduces three new modeling objects to represent robust optimization problems like Problem 2 within Pyomo: can be determined after some of the uncertainty has been resolved, i.e. y(ξ).
These three new modeling objects are sufficient for modeling quite generic robust optimization problems. They are set up to allow modeling problems in an intuitive way which is closely related to their mathematical formulation (Problem 2). We discuss each modeling object and how it is used to construct robust optimization problems in the subsequent sections.

Uncertain Parameters
Indexed uncertain parameters are constructed in analogy to Pyomo's Var type for variables: The nominal argument specifies a list of nominal valuesξ for the uncertain paramters ξ. The uncset argument specifies the uncertainty set to use for these parameters. The two approaches for constructing the uncertainty set m.U are dicussed in the next two chapters.

Generic uncertainty sets
Generic uncertainty sets are constructed with the UncSet class. This class inherits from Pyomo's Block class. Users can construct generic uncertainty sets by adding Pyomo constraints to an UncSet object in the same way as they would add constraints to a Block object in Pyomo. The following example shows how a polyhedral set can be modeled using the UncSet class: solver . solve ( m )

Library uncertainty sets
For commonly used, standard uncertainty sets, the generic approach (Section 2.2) is unnecessarily complicated. Therefore, ROmodel implements custom classes which can define an uncertainty set using its matrix representation. ROmodel implements polyhedral and ellipsoidal sets. The user can define polyhedral sets of the form U = {w | P w ≤ b} by passing the matrix P and the right hand side b to the PolyhedralSet class: In Section 4 we discuss how additional library sets can be added to ROmodel using data-driven uncertainty sets based on (warped) Gaussian processes as an example.

Constructing uncertain constraints
After defining uncertain parameters and an uncertainty set, the user can construct uncertain constraints implicitly by using the uncertain parameters in a Pyomo constraint. Consider the following deterministic Pyomo constraint: If the coefficients c are uncertain, we can model the robust constraint c T x ≤ 0, ∀c ∈ U as: The uncertain parameter m.c can be used in Pyomo constraints in the same way as Pyomo Var or Param objects. ROmodel's solvers automatically recognize constraints and objectives containing containing uncertain parameters.

Adjustable variables
ROmodel also has capabilities for modeling adjustable variables. Adjustable variables y (ξ) are variables which can be determined after some (or all) of the uncertainty has been revealed. Defining an adjustable variable is analogous to defining a regular variable in Pyomo, with an additional uncparam argument specifying a list of uncertain parameters which the adjustable variable depends on: The uncertain parameters can also be set individually for each element of the adjustable variables index using the set_uncparams function: ROmodel only implements linear decision rules for solving adjustable robust optimization problems.
If a model contains adjustable variables in a constraint or objective, ROmodel automatically replaces it by a linear decision rule based on the specified uncertain parameters.

Solvers
ROmodel has three solvers: A robust reformulation based solver, a cutting plane based solver, and a nominal solver.

Reformulation
The reformulation-based solver, illustrated in Fig. 1  If one or more constraints cannot be reformulated, the problem cannot be solved.
constraint and the corresponding uncertainty set to determine if a known reformulation is applicable.
Finally, it applies a model transformation, generating the deterministic counterpart of each robust constraint. The deterministic counterpart is then solved using an appropriate solver available in Pyomo. The structure of the optimization problem and the uncertainty set geometry determine which solvers are applicable. If no applicable transformation can be identified for one or more constraints, the problem cannot be solved and ROmodel will raise an error.
ROmodel implements standard reformulations for ellipsoidal and polyhedral uncertainty sets and linear uncertain constraints [3,5]. It also implements reformulations for black-box constrained problems [26]. These are discussed in more detail in Section 4, which also dicussed how ROmodel can be extended to include further reformulations.

Cutting planes
The cutting plane solver, outlined in Fig. 2, implements an iterative strategy for solving robust optimization problems [19]. x * is considered to be robustly feasible when for each uncertain constraint g(x * , y(ξ), ξ) ≤ 0, the objective value of the separation problem is smaller than some tolerance : ROmodel's cutting plane solver can generally be applied to any convex uncertainty set. The Pyomo solvers for solving the master and separation problems can be set individually using: The solvers need to be appropriate for the corresponding problem, i.e., the above choice would be sensible if the master problem is a mixed-integer linear problem and the uncertain constraints and uncertainty set are continuous and convex. If the solvers are not appropriate, Pyomo will raise an error.

Nominal
ROmodel also includes a nominal solver. This solver replaces all occurrences of the uncertain parameters by their nominal values and solves the resulting deterministic problem:

Extending ROmodel for black-box constrained problems
ROmodel can be extended to incorporate additional reformulations and uncertainty set geometries.
This section outlines how ROmodel can be extended using (warped) Gaussian process-based uncertainty sets for black-box constrained problems as an example. This example showcases the ease with which ROmodel can integrate Python's machine learning and data analytics capabilities with Pyomo's mathematical optimization modeling.
Wiebe et al. [26] propose a robust optimization-based approximation of a class of chanceconstraints containing uncertain black-box functions g(·): The approach models the black-box function and associated uncertainty using a (warped) Gaussian process as a stochastic model. The standard Gaussian process is well known and commonly used as a surrogate model [6]. Warped Gaussian processes are a more flexible variant of standard Gaussian process in which observations are mapped into a latent space using a non-linear, often neural netstyle warping function [22]. If a standard Gaussian process models g(·), the chance constraint Eq. 4 can be reformulated exactly. For the warped Gaussian process, Wiebe et al. [26] propose an approximation based on Wolfe duality.
In order to make these approaches available in ROmodel, we need to (i) implement two library uncertainty sets, GPSet and WarpedGPSet, which collect the relevant data, and (ii) implement two corresponding model transformations which perform the reformulations for standard and warped Gaussian process-based sets. The implementation is based on the Python module ROGP [25], which is includes Gaussian process models trained in the Python library GPY [10] in Pyomo models.

Implementing new library sets
Implementing a new library set mainly requires a new Python class collecting the necessary data.
For the standard and warped Gaussian process set, this data consist of three arguments: Constraints which use this type of uncertainty set need to be linear in the uncertain parameter.
Note that the indices of m.z and m.w need to be identical in the formulation above. If the black-box function depends on more than one variable, the Gaussian process-based sets can alternatively be specified using a dictionary: Note that ROmodel's cutting plane solver is not applicable to the Gaussian process-based sets because the sets are decision dependent. Attempting to solve a problem with one of theses sets therefore results in an error. When implementing new library sets which can be solved using cutting planes, an additional Python function generate_cons_from_lib, which generates Pyomo constraints for the uncertainty set based on the data collected by the library set, is required. For an example, see romodel/uncset/ellipsoidal.py on the ROmodel Github [28].

Implementing new reformulations
For ROmodel to be able to solve models containing the two new Gaussian process-based sets, we need to implement the corresponding reformulations. Adding new reformulations to ROmodel generally requires two Python functions: (i) a function _check_applicability which detects whether a constraint and uncertainty set have the required structure, and (ii) a function _reformulate which generates the robust counterpart. The former function is only required if the reformulation is supposed to work with generically constructed uncertainty sets as described in Section 2.2. For library sets like the Gaussian process-based sets, only the latter function is required. This function takes data describing the constraints and uncertainty set as an input and returns a Pyomo block containing the deterministic counterpart. For a full example see the implemented reformulations in romodel/reformulate/ on the ROmodel Github [28].

Results
We use ROmodel to model and solve six case studies: where the decision which facilities to build has to be made under demand uncertainty, while the decision from which facility to supply individual customers can be made once the uncertainty is resolved, 5. A production planning in which the price at which products can be sold depends on the amount produced through an uncertain black-box function modelled by a (warped) Gaussian process [26], 6. And a drill scheduling problem in which the equipment used to drill a well degrades at a rate which depends other drill parameters through a black-box function [26].
All examples except for the drill scheduling example are included with ROmodel and can be used as follows: The implementation of the drill scheduling example is separately available on Github [24]. We solve the portfolio, knapsack, pooling, and facility location problems with both the reformulation and cutting plane solver for ellipsoidal and polyhedral uncertainty sets and using both the library approach to generating uncertainty sets as well as the generic, Pyomo constraint-based approach.
We solve the production planning and drill scheduling problems using the reformulation solver with uncertainty sets based on both standard and warped Gaussian processes. We solve 30 instances with different uncertainty set sizes for each case study.  Table 1: Median time in milliseconds taken to solve the six example problems with different uncertainty set geometries. Times are shown for the reformulation and cutting plane solvers and both combined (Overall). The cutting plane solver is not applicable (NA) for Gaussian process-based sets and the reformulation solver cannot solve the facility problem with ellipsoidal sets to optimality (-). Table 1 shows the median time in milliseconds taken to solve each problem for a given uncertainty set geometry and solver. The reformulation solver generally outperforms the cutting plane solver with median times of 74ms and 330ms respectively. An exception is the the non-linear, nonconvex pooling problem with an ellipsoidal set. For this instance, the cutting plane solver achieves significantly better results, which is in line with previous work on robust pooling problems [27].
Similarly, for the facility location problem with an ellipsoidal set, the reformulation approach does not solve the problem to optimality within a 10 minute time frame, while the cutting plane solver  By construction, the ellipsoidal sets tested always fully contain the polyhedral sets for a given α.
Correspondingly, they are always more conservative, i.e., for a given α the objective value achieved using the ellipsoidal set is larger for minimization and smaller for maximization problems than the value achieved using the polyhedral set. For the Gaussian process-based sets, the standard approach is always less conservative than the warped approach. However, the limited ability of the standard Gaussian process to model non-Gaussian noise may mean that the actual probability of constraint violation is larger than the intended confidence would suggest. For a more detailed comparison of these two approaches see [26].

Conclusion
ROmodel formulates robust versions of common optimization problems. The modeling environment it provides makes (adjustable) robust optimization methods more readily available to practitioners and makes trying different solution approaches and uncertainty sets very easy. ROmodel is open source and available free of charge and could play a vital role as a platform for prototyping novel robust optimization algorithms and comparing them to existing approaches.