1 Introduction

Robust optimization is a common way of managing optimization under uncertainty in process systems engineering: applications range from production scheduling to flexible chemical process design (Janak and Floudas 2005; Li and Ierapetritou 2008; Zhang et al. 2015; Ning and You 2017; Shang and You 2018; Grossmann et al. 2016; Pistikopoulos et al. 2021). New developments in robust optimization include distributionally and adjustable robust optimization to reduce solution conservatism (Grossmann et al. 2016), data-driven robust optimization to model uncertainty sets based on available data (Bertsimas et al. 2018), and approximate robust optimization to make non-linear problems more tractable (Houska and Diehl 2013). This large number of alternative approaches and the required domain knowledge can discourage practitioners from making the transition from deterministic optimization to optimization under uncertainty. Furthermore, the lack of a platform that allows the easy implementation and application of new algorithms means that it is difficult for researchers to compare different approaches (Marc and Schöbel 2016).

This paper introduces ROmodel, a Python package extending the popular, Python-based algebraic modeling language Pyomo (Hart et al. 2011, 2017) to facilitate modeling of robust optimization problems and implementation of robust optimization algorithms. ROmodel aims to (i) make robust optimization more accessible to practitioners and facilitate moving from deterministic to robust optimization, (ii) demonstrate how robust optimization problems can be modeled in close analogy to their mathematical formulation, and (iii) provide an open-source platform for researchers to implement and compare new robust reformulations and uncertainty sets. To this end, ROmodel introduces intuitive modeling which closely resembles the optimization problem’s underlying mathematical formulation and will feel familiar to Pyomo users. ROmodel uses the richness of Pyomo’s solver interfaces and methods for model transformations and Python’s data processing capabilities. It supports both automatic reformulation and cutting plane algorithms. Uncertainty sets can be chosen from a library of common geometries, or custom defined using Pyomo constraints. ROmodel also support adjustable robust optimization through linear decision rules. Because ROmodel and Pyomo are open-source, ROmodel can be extended to incorporate additional uncertainty set geometries and reformulations. As an example, we implemented uncertainty sets based on (warped) Gaussian processes for black-box constrained problems.

Robust optimization has been applied to many types of engineering problems. It has, for example, extensively been applied to scheduling and planning problems in process systems engineering (Janak and Floudas 2005; Li and Ierapetritou 2008; Ning and You 2017; Shang and You 2018; Wiebe et al. 2018), but also to robust design problems (Yuan et al. 2018; Gong and You 2018) and robust control (Diehl et al. 2008). Types of uncertainties considered have ranged from uncertain demands (Li et al. 2012) and uncertain prices (Wiebe et al. 2020) to uncertain rates of equipment degradation (Wiebe et al. 2018, 2020) and uncertain product/reagent purities and consumption rates (Zhang et al. 2018; Wiebe et al. 2019). A large number of these applications use ellipsoidal and/or polyhedral uncertainty sets which are available as library sets in ROmodel (Li et al. 2012; Shang et al. 2017; Zhang et al. 2018; Wiebe et al. 2018, 2019). Other geometries could be accomodated using ROmodel’s capabilities for modeling generic convex sets using Pyomo constraints or by extending ROmodel to new, custom uncertainty sets. In many engineering applications, the optimal choice of uncertainty set geometry and solution method (reformulation or cutting plane) is not obvious before solving the problem. With ROmodel, engineers can implement their robust model once and then solve it using different uncertainty set geometries and solvers. While most engineering applications have been LP and MILP problems, robust optimization has also been applied to non-linear engineering problems (Diehl et al. 2008; Li et al. 2012; Gong and You 2018; Wiebe et al. 2019). ROmodel can solve both MILP and MINLP problems depending on the Pyomo solver interface used.

There are other software packages for solving robust optimization problems, e.g. ROME (Goh and Sim 2011), RSOME (Chen et al. 2020), ROC++ (Vayanos et al. 2020), and PyROS Isenberg et al. (2020). All four approaches are designed as robust optimization solvers and, except for PyROS which is based on Pyomo, do not rely on a general purpose algebraic modeling language. In contrast, ROmodel focuses on modeling robust optimization problems. By building on Pyomo, ROmodel simplifies the transition from deterministic to robust optimization, since both the deterministic and robust model can be implemented in the same environment. ROmodel also allows access to a much larger number of deterministic solvers than these other solvers. Unlike PyROS, which uses Pyomo’s Var component to model uncertain parameters and adjustable variables, ROmodel introduces new, dedicated Pyomo modeling components for robust problems allowing for more intuitive modeling. ROmodel also supports both the cutting plane and reformulation methods as well as modeling custom uncertainty sets with Pyomo constraints. In contrast to AIMMS, which has some capabilities for modeling robust optimization problems (AIMMS 2021), ROModel is open source and allows possibilities for extension. Ju MPeR Dunning (2016), which extends Julia’s modeling language Ju MP to robust optimization problems, is the most similar to ROModel. In contrast to JuMPeR, which uses polyhedral sets by default, ROModel automatically recognizes uncertainty set geometries and applies appropriate reformulations, if available, without the user having to specify which geometry they are using. ROmodel also supports both uncertain constraints and objectives while JuMPeR only supports uncertain constraints.

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 Wiebe et al. (2020). These sets seemlesly integrate Gaussian processes trained in the Python library GPy (2012) into robust optimization problems. ROmodel is open source and available on Wiebe et al. (2020).

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.

2 Modeling

Consider the following generic deterministic optimization problem

$$\begin{aligned}&\min \limits _{\varvec{x}\in \mathcal {X} , \varvec{y}}\quad&f(\varvec{x}, \varvec{y}, \bar{\varvec{\xi }}) \end{aligned}$$
(1a)
$$\begin{aligned}&\text {s.t} \quad&g(\varvec{x}, \varvec{y}, \bar{\varvec{\xi }}) \le 0 \end{aligned}$$
(1b)

where \(\varvec{x}\in \mathbb {R}^n\) and \(\varvec{y}\in R^m\) are decision variables and \(\bar{\varvec{\xi }} \in \mathbb {R}^k\) is a vector of (nominal) parameters. If the parameter vector \(\bar{\varvec{\xi }}\) is not known exactly, we can construct the following robust version of Problem 1:

$$\begin{aligned}&\min \limits _{\varvec{x}\in \mathcal {X} , \varvec{y}(\varvec{\xi })}\max \limits _{\varvec{\xi }\in \mathcal {U}(\varvec{x})}\quad&f(\varvec{x}, \varvec{y}(\varvec{\xi }), \varvec{\xi }) \end{aligned}$$
(2a)
$$\begin{aligned}&\text {s.t} \quad&g(\varvec{x}, \varvec{y}(\varvec{\xi }), \varvec{\xi }) \le 0&\forall \varvec{\xi }\in \mathcal {U}(\varvec{x}) \end{aligned}$$
(2b)

Here \(\varvec{x}\in \mathbb {R}^n\) are“here and now”variables, determined before the uncertainty is revealed, while \(\varvec{y}(\varvec{\xi })\) are adjustable variables, determined after the uncertainty is revealed. The uncertain parameter vector \(\varvec{\xi }\) is bounded by the uncertainty set \(\mathcal {U}(\varvec{x})\), which may depend on \(\varvec{x}\) and which contains the nominal values \(\bar{\varvec{\xi }}\) 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 \(\mathcal {X}\).

ROmodel introduces three new modeling objects to represent robust optimization problems like Problem 2 within Pyomo:

  1. 1

    UncParam: A class similar to Pyomo’s Param and Var class used to model uncertain parameters \(\varvec{\xi }\). One can supply a nominal argument, which defines the vector of nominal values \(\bar{\varvec{\xi }}\) used to replace the uncertain parameters when solving the deterministic (nominal) problem (Eq. 1).

  2. 2

    UncSet: Each UncParam object has an UncSet object associated with it. The UncSet class, based on Pyomo’s Block class, models the uncertainty sets \(\mathcal {U}\). Uncertainty sets \(\mathcal {U}\) can be defined by (i) constructing generic sets by adding Pyomo constraints to the UncSet object, or (ii) through a library of common uncertainty set geometries, using their matrix representation as an input.

  3. 3

    AdjustableVar: A class similar to Pyomo’s Var class used to model adjustable variables which can be determined after some of the uncertainty has been resolved, i.e. \(\varvec{y}(\varvec{\xi })\).

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.

2.1 Uncertain Parameters

Indexed uncertain parameters are constructed in analogy to Pyomo’s Var type for variables:

figure a

The nominal argument specifies a list of nominal values \(\bar{\varvec{\xi }}\) for the uncertain paramters \(\varvec{\xi }\). The nominal values are currently only used in conjuction with the nominal solver outlined in Sect. 3.3. ROmodel does not check if the nominal values are contained in the uncertainty set. 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.

2.2 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:

figure b

The uncset and nominal arguments define the uncertainty set and a vector of nominal values for the uncertain parameter m.w. ROmodel’s strategy of modeling uncertainty sets using Pyomo’s Constraint modeling object is analogous to typical robust optimization formulations. Users can also define multiple uncertainty sets and replace one by another:

figure c

2.3 Library uncertainty sets

For commonly used, standard uncertainty sets, the generic approach (Sect. 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 \(\mathcal {U} = \left\{ w \; | \; Pw \le b\right\}\) by passing the matrix P and the right hand side b to the PolyhedralSet class:

figure d

ROmodel creates ellipsoidal sets of the form \((w-\mu )\Sigma ^{-1}(w-\mu ) \le 1\) using the EllipsoidalSet class, the covariance matrix \(\Sigma\) and the mean vector \(\mu\):

figure e

Note that ROmodel does not currently have capabilities for taking intersections of (library) uncertainty sets. Library sets specifying intersections of common uncertainty set geometries could be added to ROmodel. In Sect. 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.

2.4 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:

figure f

If the coefficients c are uncertain, we can model the robust constraint \(c^T x \le 0, \; \forall c \in \mathcal {U}\) as:

figure g

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 uncertain parameters.

2.5 Adjustable variables

ROmodel also has capabilities for modeling adjustable variables. Adjustable variables \(\varvec{y}\left( \varvec{\xi }\right)\) 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:

figure h

The uncertain parameters can also be set individually for each element of the adjustable variables index using the set_uncparams function:

figure i

This means that multi-stage problems, in which each adjustable variable depends on the uncertain parameters from previous stages, can be modeled in ROmodel. 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.

3 Solvers

ROmodel’s UncSet component can be used to model any uncertainty set geometry which can be expressed using Pyomo constraints. However, not every geometry that can be modeled can necessarily also be solved. ROmodel has three solvers: A robust reformulation based solver, a cutting plane based solver, and a nominal solver. This section discusses how each of these solvers works and which types of problems and uncertainty set geometries they can solve.

3.1 Reformulation

The reformulation-based solver, illustrated in Fig. 1, implements standard duality based techniques for reformulating robust optimization problems into deterministic counterparts (Bertsimas and Sim 2004).

Fig. 1
figure 1

Schematic of reformulation solver. For each uncertain constraint, the reformulation solver goes through all known transformations and identifies the geometry of the constraint/uncertainty set. It then applies the corresponding model transformation. If all uncertain constraints can be reformulated, the resulting deterministic counterpart is solved using one of Pyomo’s solver interfaces. If one or more constraints cannot be reformulated, the problem cannot be solved.

First, it detects every constraint containing uncertain parameters. Second, it checks the structure of each uncertain 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 (Ben-Tal and Nemirovski 1999; Bertsimas and Sim 2004). It also implements reformulations for black-box constrained problems (Wiebe et al. 2020). These are discussed in more detail in Sect. 4, which also discusses how ROmodel can be extended to include further reformulations.

3.2 Cutting planes

The cutting plane solver, outlined in Fig. 2, implements an iterative strategy for solving robust optimization problems (Mutapcic and Boyd 2009). It replaces each uncertain constraint and objective by a CutGenerator object which initially just contains the nominal constraint. The solver then iteratively solves the master problem and generates cuts to cut off solutions which are not robustly feasible.

Fig. 2
figure 2

Schematic of cutting plane solver. The cutting plane solver iteratively solves the master problem and adds cutting planes until a robustly feasible solution is found.

A solution \(x^*\) is considered to be robustly feasible when for each uncertain constraint \(g(\varvec{x}^*, \varvec{y}(\varvec{\xi }), \varvec{\xi }) \le 0\), the objective value of the separation problem is smaller than some tolerance \(\epsilon\):

$$\begin{aligned} \max \limits _{\varvec{\xi }\in \mathcal {U}} \;&g(\varvec{x}^*, \varvec{y}(\varvec{\xi }), \varvec{\xi }) \le \epsilon \end{aligned}$$
(3)

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:

figure j

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.

3.3 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:

figure k

The nominal solver allows users to combine their implementations of the nominal and robust problem. An implementation of the robust model can be used to obtain the solution of the nominal problem.

4 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. 2020) propose a robust optimization-based approximation of a class of chance-constraints containing uncertain black-box functions \(g(\cdot )\):

$$\begin{aligned} P\left( \sum \limits _i g(\varvec{y}_i)x_i \le b\right) \ge 1 - \alpha \end{aligned}$$
(4)

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 (Bhosekar and Ierapetritou 2018). 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 net-style warping function (Snelson et al. 2003). If a standard Gaussian process models \(g(\cdot )\), the chance constraint Eq. 4 can be reformulated exactly. For the warped Gaussian process, (Wiebe et al. 2020) propose an approximation based on Wolfe duality. Choosing between the standard and warped approach for a specific application is not trivial. We therefore implemented both approaches in ROmodel, so that practitioners can easily apply both to their application and compare.

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 (Wiebe 2020), which includes Gaussian process models trained in the Python library (GPy 2012) in Pyomo models.

4.1 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:

figure l

The first argument gp_standard/warped is a (warped) Gaussian process object trained in GPy. The second is an indexed Pyomo variable on which the GP depends, i.e. \(\varvec{y}\) in Eq. 4. The third parameter specifies the confidence level \(1-\alpha\) with which the true parameter is contained in the uncertainty set. I.e., in this case the confidence that the true parameter vector is an element of the uncertainty set is at least 95%.

The new sets can be used in the same way as other library sets, e.g.:

figure m

Note that the uncertain parameter m.w corresponds to \(g(\varvec{y}_i)\) in Formulation 4. Constraints which use this type of uncertainty set need to be linear in the uncertain parameter. Note that the indices of m.y 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:

figure n

The dictionary indicates that the uncertain parameter m.w[0] depends on the variables m.y[0, ’a’] and m.y[0, ’b’] through the black-box function \(g(\cdot )\), modeled in GPy by the Gaussian process gp.

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. As an example, the following code shows the implementation of the polyhedral set in ROmodel:

figure o

4.2 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:

  1. 1

    A function _check_applicability which detects whether a constraint and uncertainty set have the required structure, and

  2. 2

    A function _reformulate which generates the robust counterpart.

The function _check_applicability iterates through the constraints of an uncertainty set to identify its geometry. E.g., for polyhedral sets, it checks that all constraints are linear. This function is only required if the reformulation is supposed to work with generically constructed uncertainty sets as described in Sect. 2.2. For library sets like the Gaussian process-based sets, only the _reformulate 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. E.g., the transformation for warped Gaussian process-based sets returns a block containing the Wolfe dual formulation described in Wiebe et al. (2020) while the transformation for polyhedral sets returns a block containing the linear dual. For a full example see the implemented reformulations in romodel/reformulate/ on the (Wiebe et al. 2020).

5 Results

We use ROmodel to model and solve six case studies:

  1. (1)

    A portfolio optimization problem with uncertain returns (Bertsimas and Sim 2004),

  2. (2)

    A knapsack problem with uncertain item weights,

  3. (3)

    A pooling problem instance (Adhya et al. 1999) with uncertain product demands,

  4. (4)

    A capacitated facility location problem as an example for adjustable robust optimization, 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. (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 (Wiebe et al. 2020),

  6. (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 (Wiebe et al. 2020).

For the reformulation approach, we use Gurobi 9.0.3 to solve the reformulation for case studies 1–4 and Ipopt 3.12.12 for case studies 5 and 6. For the cutting plane approach, we use Gurobi 9.0.3 as the main and subsolver for all problems. All experiments were run on an i7-6700 CPU with 8x3.40GHz and 16GB RAM. All examples except for the drill scheduling example are included with ROmodel and can be used as follows:

figure p

The implementation of the drill scheduling example is separately available on (Wiebe 2020). Table 1 shows a summary of the type of problem and number of variables, constraints, and uncertain parameters for each case study.

Table 1 Type of problem, number of variables, constraints, and uncertain parameters for each case study.

Below we discuss two case studies, the Knapsack and Facility location problems, in more detail. Detailed code for the Portfolio, Pooling, and Production planning case studies can be cound on the ROmodel Github at romodel/examples.

5.1 Knapsack problem

Our first case study is a classic Knapsack problem:

$$\begin{aligned} \max \limits _{\varvec{x}} \quad&\sum \limits _i v_i x_i \end{aligned}$$
(5a)
$$\begin{aligned} \text {s.t.} \quad&w_i x_i \le b&\forall \varvec{w} \in \mathcal {U} \end{aligned}$$
(5b)
$$\begin{aligned}&x_i \in \{0, 1\} \end{aligned}$$
(5c)

The goal is to select from a list of items i with value \(v_i\) and weight \(w_i\) such that the total value selected is maximized while the total weight is below some limit b. The binary variable \(x_i\) is 1 if the item is selected and 0 otherwise. We consider the vector of weights \(\varvec{w}\) to be uncertain but bounded by the uncertainty set \(\mathcal {U}\). The following code shows how the problem can be modeled and solved in ROmodel using an ellipsoidal uncertainty set:

figure q

Note that unlike other robut optimization solvers, ROmodel’s modeling components behave exactly like Pyomo’s modeling components. Practitioners who are used to implementing their deterministic models in Pyomo will find the syntax of ROmodel very familiar. Using ROmodel, they can apply robust optimization to their existing Pyomo models by making minimal changes to the model.

The results can be accessed as in any other Pyomo model. E.g. the variable values can be printed using m.x.pprint(). When using the reformulation solver, the reformulation for each uncertain constraint can be accessed using the postfix _counterpart. E.g. for the problem above m.weight_counterpart.pprint() gives the output:

figure r

In this case, the counterpart consists of one additional variable, padding, and two constraints det and rob.

For the cutting plane solvers, the cuts for each uncertain constraint can be accessed using the postfix _generator, e.g. m.weight_generator.pprint() outputs a list of two cuts:

figure s

5.2 Facility location problem

The facility location case study is an adjustable robust optimization problem with the following mathematical formulation:

$$\begin{aligned} \max \limits _{\varvec{x}}\min \limits _{\varvec{d} \in \mathcal {U}} \quad&\sum \limits _i^N \sum \limits _j^M c^t_{ij}y_{ij}(d_j) + \sum \limits _i^N c^f_i x_i \end{aligned}$$
(6a)
$$\begin{aligned} \text {s.t.} \quad&\sum \limits _i^N y_{ij}(d_j) = d_j&\forall j \in [M], \varvec{d} \in \mathcal {U}\end{aligned}$$
(6b)
$$\begin{aligned}&\sum \limits _j^M y_{ij}(d_j) \le m_i x_i&\forall i \in [N], \varvec{d} \in \mathcal {U}\end{aligned}$$
(6c)
$$\begin{aligned}&\varvec{x}_i \in \{0, 1\} \end{aligned}$$
(6d)

Variables \(x_i\) are here and now binary variables which are 1 if facility i is built and 0 otherwise. Variables \(y_{ij}(d_j)\) are adjustable variables which depend on the uncertain customer demand \(d_j\) for each customer j and determine how much of the customer’s demand is satisfied from facility i. The maximum capacity for each facility is given by \(m_i\). The uncertain parameters \(d_i\) are assumed to be bounded by the uncertainty set \(\mathcal {U}\). The code below shows how ROmodel can be used to implement this problem using the generic uncertainty set approach to define a box uncertainty set:

figure tfigure t

The coefficients of the linear decision rules can be accessed through the name of the adjustable variable, the name of the uncertain parameter and the postfix _coef. E.g. for the problem above m.y_demand_coef.pprint() will print:

figure u

5.3 Solution times

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 2 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 median time taken to reformulate the model is shown in parantheses. 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 2 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, non-convex 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 (Wiebe et al. 2019). 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 does. For the production planning and drills scheduling examples only the reformulation solver can be applied. The Wolfe duality-based reformulation for warped Gaussian processes generally takes longer to solve than the chance-constraint reformulation for standard Gaussian processes. Note that most of this time is the time taken by the subsolvers. The transformations which ROmodel performs are generally very quick: the median transformation time across all instances is 3.2 milliseconds, while the maximum transformation time is 1.7 seconds. The median transformation time for each problem and geometry/solver is shown in parantheses in Table 2.

6 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.