# PyROS Solver

PyROS (Pyomo Robust Optimization Solver) is a metasolver capability within Pyomo for solving non-convex, two-stage optimization models using adjustable robust optimization.

It was developed by Natalie M. Isenberg and Chrysanthos E. Gounaris of Carnegie Mellon University, in collaboration with John D. Siirola of Sandia National Labs. The developers gratefully acknowledge support from the U.S. Department of Energy’s Institute for the Design of Advanced Energy Systems (IDAES).

## Methodology Overview

Below is an overview of the type of optimization models PyROS can accomodate.

• PyROS is suitable for optimization models of continuous variables that may feature non-linearities (including non-convexities) in both the variables and uncertain parameters.
• PyROS can handle equality constraints defining state variables, including implicit state variables that cannot be eliminated via reformulation.
• PyROS allows for two-stage optimization problems that may feature both first-stage and second-stage degrees of freedom.

The general form of a deterministic optimization problem that can be passed into PyROS is shown below:

\begin{split}\begin{align*} \displaystyle \min_{\substack{x \in \mathcal{X}, \\ z \in \mathbb{R}^n, y\in\mathbb{R}^a}} & ~~ f_1\left(x\right) + f_2\left(x,z,y; q^0\right) & \\ \displaystyle \text{s.t.} \quad \: & ~~ g_i\left(x, z, y; q^0\right) \leq 0 & \forall i \in \mathcal{I} \\ & ~~ h_j\left(x,z,y; q^0\right) = 0 & \forall j \in \mathcal{J} \\ \end{align*}\end{split}

where:

• $$x \in \mathcal{X}$$ are the “design” variables (i.e., first-stage degrees of freedom), where $$\mathcal{X} \subseteq \mathbb{R}^m$$ is the feasible space defined by the model constraints that only reference these variables
• $$z \in \mathbb{R}^n$$ are the “control” variables (i.e., second-stage degrees of freedom)
• $$y \in \mathbb{R}^a$$ are the “state” variables
• $$q \in \mathbb{R}^w$$ is the vector of parameters that we shall later consider to be uncertain, and $$q^0$$ is the vector of nominal values associated with those.
• $$f_1\left(x\right)$$ are the terms of the objective function that depend only on design variables
• $$f_2\left(x, z, y; q\right)$$ are the terms of the objective function that depend on control and/or state variables
• $$g_i\left(x, z, y; q\right)$$ is the $$i^\text{th}$$ inequality constraint in set $$\mathcal{I}$$ (see Note)
• $$h_j\left(x, z, y; q\right)$$ is the $$j^\text{th}$$ equality constraint in set $$\mathcal{J}$$ (see Note)

Note

• Applicable bounds on variables $$z$$ and/or $$y$$ are assumed to have been incorporated in the set of inequality constraints $$\mathcal{I}$$.
• A key requirement of PyROS is that each value of $$\left(x, z, q \right)$$ maps to a unique value of $$y$$, a property that is assumed to be properly enforced by the system of equality constraints $$\mathcal{J}$$. If such unique mapping does not hold, then the selection of ‘state’ (i.e., not degree of freedom) variables $$y$$ is incorrect, and one or more of the $$y$$ variables should be appropriately redesignated to be part of either $$x$$ or $$z$$.

In order to cast the robust optimization counterpart formulation of the above model, we shall now assume that the uncertain parameters may attain any realization from within an uncertainty set $$\mathcal{Q} \subseteq \mathbb{R}^w$$, such that $$q^0 \in \mathcal{Q}$$. The set $$\mathcal{Q}$$ is assumed to be closed and bounded, while it can be either continuous or discrete.

Based on the above notation, the form of the robust counterpart addressed in PyROS is shown below:

\begin{split}\begin{align*} \displaystyle \min_{x \in \mathcal{X}} & \displaystyle \max_{q \in \mathcal{Q}} & \displaystyle \min_{z \in \mathbb{R}^n, y \in \mathbb{R}^a} \ \ & \displaystyle ~~ f_1\left(x\right) + f_2\left(x, z, y, q\right) & & \\ & & \text{s.t.} \quad \:& \displaystyle ~~ g_i\left(x, z, y, q\right) \leq 0 & & \forall i \in \mathcal{I}\\ & & & \displaystyle ~~ h_j\left(x, z, y, q\right) = 0 & & \forall j \in \mathcal{J} \end{align*}\end{split}

In order to solve problems of the above type, PyROS implements the Generalized Robust Cutting-Set algorithm developed in [GRCSPaper].

When using PyROS, please consider citing the above paper.

## PyROS Required Inputs

The required inputs to the PyROS solver are the following:

• The determinisitic optimization model
• List of first-stage (“design”) variables
• List of second-stage (“control”) variables
• List of parameters to be considered uncertain
• The uncertainty set
• Subordinate local and global NLP optimization solvers

Below is a list of arguments that PyROS expects the user to provide when calling the solve command. Note how all but the model argument must be specified as kwargs.

modelConcreteModel
A ConcreteModel object representing the deterministic model.
first_stage_variableslist(Var)
A list of Pyomo Var objects representing the first-stage degrees of freedom (design variables) in model.
second_stage_variableslist(Var)
A list of Pyomo Var objects representing second-stage degrees of freedom (control variables) in model.
uncertain_paramslist(Param)
A list of Pyomo Param objects in deterministic_model to be considered uncertain. These specified Param objects must have the property mutable=True.
uncertainty_setUncertaintySet
A PyROS UncertaintySet object representing uncertainty in the space of those parameters listed in the uncertain_params object.
local_solverSolver
A Pyomo Solver instance for a local NLP optimization solver.
global_solverSolver
A Pyomo Solver instance for a global NLP optimization solver.

Note

Any variables in the model not specified to be first- or second-stage variables are automatically considered to be state variables.

## PyROS Solver Interface

class pyomo.contrib.pyros.PyROS[source]

PyROS (Pyomo Robust Optimization Solver) implementing a generalized robust cutting-set algorithm (GRCS) to solve two-stage NLP optimization models under uncertainty.

solve(model, first_stage_variables, second_stage_variables, uncertain_params, uncertainty_set, local_solver, global_solver, **kwds)[source]

Solve the model.

Parameters: Keyword Arguments: model (ConcreteModel) – A ConcreteModel object representing the deterministic model, cast as a minimization problem. first_stage_variables (List[Var]) – The list of Var objects referenced in model representing the design variables. second_stage_variables (List[Var]) – The list of Var objects referenced in model representing the control variables. uncertain_params (List[Param]) – The list of Param objects referenced in model representing the uncertain parameters. MUST be mutable. Assumes entries are provided in consistent order with the entries of ‘nominal_uncertain_param_vals’ input. uncertainty_set (UncertaintySet) – UncertaintySet object representing the uncertainty space that the final solutions will be robust against. local_solver (Solver) – Solver object to utilize as the primary local NLP solver. global_solver (Solver) – Solver object to utilize as the primary global NLP solver. time_limit – Optional. Default = None. Total allotted time for the execution of the PyROS solver in seconds (includes time spent in sub-solvers). ‘None’ is no time limit. keepfiles – Optional. Default = False. Whether or not to write files of sub-problems for use in debugging. Must be paired with a writable directory supplied via subproblem_file_directory. tee – Optional. Default = False. Sets the tee for all sub-solvers utilized. load_solution – Optional. Default = True. Whether or not to load the final solution of PyROS into the model object. objective_focus – Optional. Default = ObjectiveType.nominal. Choice of objective function to optimize in the master problems. Choices are: ObjectiveType.worst_case, ObjectiveType.nominal. See Note for details. nominal_uncertain_param_vals – Optional. Default = deterministic model Param values. List of nominal values for all uncertain parameters. Assumes entries are provided in consistent order with the entries of uncertain_params input. decision_rule_order – Optional. Default = 0. Order of decision rule functions for handling second-stage variable recourse. Choices are: ‘0’ for constant recourse (a.k.a. static approximation), ‘1’ for affine recourse (a.k.a. affine decision rules), ‘2’ for quadratic recourse. solve_master_globally – Optional. Default = False. ‘True’ for the master problems to be solved with the user-supplied global solver(s); or ‘False’ for the master problems to be solved with the user-supplied local solver(s). max_iter – Optional. Default = -1. Iteration limit for the GRCS algorithm. ‘-1’ is no iteration limit. robust_feasibility_tolerance – Optional. Default = 1e-4. Relative tolerance for assessing robust feasibility violation during separation phase. separation_priority_order – Optional. Default = {}. Dictionary mapping inequality constraint names to positive integer priorities for separation. Constraints not referenced in the dictionary assume a priority of 0 (lowest priority). progress_logger – Optional. Default = “pyomo.contrib.pyros”. The logger object to use for reporting. backup_local_solvers – Optional. Default = []. List of additional Solver objects to utilize as backup whenever primary local NLP solver fails to identify solution to a sub-problem. backup_global_solvers – Optional. Default = []. List of additional Solver objects to utilize as backup whenever primary global NLP solver fails to identify solution to a sub-problem. subproblem_file_directory – Optional. Path to a directory where subproblem files and logs will be written in the case that a subproblem fails to solve. bypass_local_separation – This is an advanced option. Default = False. ‘True’ to only use global solver(s) during separation; ‘False’ to use local solver(s) at intermediate separations, using global solver(s) only before termination to certify robust feasibility. bypass_global_separation – This is an advanced option. Default = False. ‘True’ to only use local solver(s) during separation; however, robustness of the final result will not be guaranteed. Use to expedite PyROS run when global solver(s) cannot (efficiently) solve separation problems. – p_robustness – This is an advanced option. Default = {}. Whether or not to add p-robustness constraints to the master problems. If the dictionary is empty (default), then p-robustness constraints are not added. See Note for how to specify arguments.

Note

Solving the master problems globally (via option solve_masters_globally=True) is one of the requirements to guarantee robust optimality; solving the master problems locally can only lead to a robust feasible solution.

Note

Selecting worst-case objective (via option objective_focus=ObjectiveType.worst_case) is one of the requirements to guarantee robust optimality; selecting nominal objective can only lead to a robust feasible solution, albeit one that has optimized the sum of first- and (nominal) second-stage objectives.

Note

To utilize option p_robustness, a dictionary of the following form must be supplied via the kwarg: There must be a key (str) called ‘rho’, which maps to a non-negative value, where ‘1+rho’ defines a bound for the ratio of the objective that any scenario may exhibit compared to the nominal objective.

## PyROS Uncertainty Sets

PyROS contains pre-implemented UncertaintySet specializations for many types of commonly used uncertainty sets. Additional capabilities for intersecting multiple PyROS UncertaintySet objects so as to create custom sets are also provided via the IntersectionSet class. Custom user-specified sets can also be defined via the base UncertaintySet class.

Mathematical representations of the sets are shown below, followed by the class descriptions.

Table 1 PyROS Uncertainty Sets
Uncertainty Set Type Set Representation
BoxSet $$Q_X = \left\{q \in \mathbb{R}^n : q^\ell \leq q \leq q^u\right\} \\ q^\ell \in \mathbb{R}^n \\ q^u \in \mathbb{R}^n : \left\{q^\ell \leq q^u\right\}$$
CardinalitySet $$Q_C = \left\{q \in \mathbb{R}^n : q = q^0 + (\hat{q} \circ \xi) \text{ for some } \xi \in \Xi_C\right\}\\ \Xi_C = \left\{\xi \in [0, 1]^n : \displaystyle\sum_{i=1}^{n} \xi_i \leq \Gamma\right\} \\ \Gamma \in [0, n] \\ \hat{q} \in \mathbb{R}^{n}_{+} \\ q^0 \in \mathbb{R}^n$$
BudgetSet $$Q_B = \left\{q \in \mathbb{R}^n_+: \displaystyle\sum_{i \in B_\ell} q_i \leq b_\ell \ \forall \ell \in \left\{1,\ldots,L\right\} \right\} \\ b_\ell \in \mathbb{R}^{L}_+$$
FactorModelSet $$Q_F = \left\{q \in \mathbb{R}^n: \displaystyle q = q^0 + \Psi \xi \text{ for some }\xi \in \Xi_F\right\} \\ \Xi_F = \left\{ \xi \in \left[-1, 1\right]^F, \left\lvert \displaystyle \sum_{f=1}^{F} \xi_f\right\rvert \leq \beta F \right\} \\ \beta \in [0,1] \\ \Psi \in \mathbb{R}^{n \times F}_+ \\ q^0 \in \mathbb{R}^n$$
PolyhedralSet $$Q_P = \left\{q \in \mathbb{R}^n: \displaystyle A q \leq b \right\} \\ A \in \mathbb{R}^{m \times n} \\ b \in \mathbb{R}^{m} \\ q^0 \in \mathbb{R}^n: {Aq^0 \leq b}$$
AxisAlignedEllipsoidalSet $$Q_A = \left\{q \in \mathbb{R}^n: \displaystyle \sum\limits_{i=1 : \atop \left\{ \alpha_i > 0 \right\} } \left(\frac{q_i - q_i^0}{\alpha_i} \right)^2 \leq 1 , \quad q_i = q^0_i \quad \forall i : \left\{\alpha_i=0\right\}\right\} \\ \alpha \in \mathbb{R}^n_+, \\ q^0 \in \mathbb{R}^n$$
EllipsoidalSet $$Q_E = \left\{q \in \mathbb{R}^n: \displaystyle q = q^0 + P^{1/2} \xi \text{ for some } \xi \in \Xi_E \right\} \\ \Xi_E = \left\{\xi \in \mathbb{R} : \xi^T\xi \leq s \right\} \\ P \in \mathbb{S}^{n\times n}_+ \\ s \in \mathbb{R}_+ \\ q^0 \in \mathbb{R}^n$$
UncertaintySet $$Q_U = \left\{q \in \mathbb{R}^n: \displaystyle g_i(q) \leq 0 \quad \forall i \in \left\{1,\ldots,m \right\}\right\} \\ m \in \mathbb{N}_+ \\ g_i : \mathbb{R}^n \mapsto \mathbb{R} \, \forall i \in \left\{1,\ldots,m\right\}, \\ q^0 \in \mathbb{R}^n : \left\{g_i(q^0) \leq 0 \ \forall i \in \left\{1,\ldots,m\right\}\right\}$$
DiscreteScenariosSet $$Q_D = \left\{q^s : s = 0,\ldots,D \right\} \\ D \in \mathbb{N} \\ q^s \in \mathbb{R}^n \forall s \in \left\{ 0,\ldots,D\right\}$$
IntersectionSet $$Q_I = \left\{q \in \mathbb{R}^n: \displaystyle q \in \bigcap_{i \in \left\{1,\ldots,m\right\}} Q_i\right\} \\ Q_i \subset \mathbb{R}^n \quad \forall i \in \left\{1,\ldots,m\right\}$$

Note

Each of the PyROS uncertainty set classes inherits from the UncertaintySet base class.

### PyROS Uncertainty Set Classes

class pyomo.contrib.pyros.uncertainty_sets.BoxSet(bounds)[source]

A hyper-rectangle (a.k.a “box”).

Parameters: bounds ((N, 2) array_like) – Lower and upper bounds for each dimension of the set.
property bounds

Lower and upper bounds for each dimension of the set.

The bounds of a BoxSet instance can be changed, such that the dimension of the set remains unchanged.

Type: (N, 2) numpy.ndarray
property dim

Dimension of the box set.

Type: int
property parameter_bounds

Uncertain parameter value bounds for the box set.

Returns: Box set bounds. list(tuple)
point_in_set(point)

Determine whether a given point lies in the uncertainty set.

Parameters: point ((N,) array-like) – Point (parameter value) of interest. is_in_set – True if the point lies in the uncertainty set, False otherwise. bool

Notes

This method is invoked at the outset of a PyROS solver call to determine whether a user-specified nominal parameter realization lies in the uncertainty set.

property type

Brief description of the type of the uncertainty set.

Type: str
class pyomo.contrib.pyros.uncertainty_sets.CardinalitySet(origin, positive_deviation, gamma)[source]

A cardinality-constrained (a.k.a. “gamma”) set.

Parameters: origin ((N,) array_like) – Origin of the set (e.g. nominal uncertain parameter values). positive_deviation ((N,) array_like) – Maximal non-negative coordinate deviation from the origin in each dimension. gamma (numeric type) – Upper bound for the number of uncertain parameters which may realize their maximal deviations from the origin simultaneously.
property dim

Dimension of the cardinality set.

Type: int
property gamma

Upper bound for the number of uncertain parameters which may maximally deviate from their respective origin values simultaneously. Must be a numerical value ranging from 0 to the set dimension.

Note that mathematically, setting gamma to 0 reduces the set to a singleton containing the center, while setting gamma to the set dimension reduces the set to a hyperrectangle with bounds [origin, origin + positive_deviation].

Type: numeric type
property origin

Origin of the cardinality set (e.g. nominal parameter values).

Type: (N,) numpy.ndarray
property parameter_bounds

Uncertain parameter value bounds for the cardinality set.

Returns: parameter_bounds – A list of 2-tuples of numerical values. Each tuple specifies the uncertain parameter bounds for the corresponding set dimension. list of tuples
point_in_set(point)[source]

Determine whether a given point lies in the cardinality set.

Parameters: point ((N,) array-like) – Point (parameter value) of interest. True if the point lies in the set, False otherwise. bool
property positive_deviation

Maximal coordinate deviations from the origin in each dimension. All entries are nonnegative.

Type: (N,) numpy.ndarray
property type

Brief description of the type of the uncertainty set.

Type: str
class pyomo.contrib.pyros.uncertainty_sets.BudgetSet(budget_membership_mat, rhs_vec, origin=None)[source]

A budget set.

Parameters: budget_membership_mat ((M, N) array_like) – Incidence matrix of the budget constraints. Each row corresponds to a single budget constraint, and defines which uncertain parameters (which dimensions) participate in that row’s constraint. rhs_vec ((M,) array_like) – Budget limits (upper bounds) with respect to the origin of the set. origin ((N,) array_like or None, optional) – Origin of the budget set. If None is provided, then the origin is set to the zero vector.
property budget_membership_mat

Incidence matrix of the budget constraints. Each row corresponds to a single budget constraint, and defines which uncertain parameters (which dimensions) participate in that row’s constraint.

Type: (M, N) numpy.ndarray
property budget_rhs_vec

Budget limits (upper bounds) with respect to the origin.

Type: (M,) numpy.ndarray
property coefficients_mat

Coefficient matrix of all polyhedral constraints defining the budget set. Composed from the incidence matrix used for defining the budget constraints and a coefficient matrix for individual uncertain parameter nonnegativity constraints.

This attribute cannot be set. The budget constraint incidence matrix may be altered through the budget_membership_mat attribute.

Type: (M + N, N) numpy.ndarray
property dim

Dimension of the budget set.

Type: int
property parameter_bounds

Uncertain parameter value bounds for the budget set.

Returns: parameter_bounds – A list of 2-tuples of numerical values. Each tuple specifies the uncertain parameter bounds for the corresponding set dimension. list of tuples
point_in_set(point)

Determine whether a given point lies in the uncertainty set.

Parameters: point ((N,) array-like) – Point (parameter value) of interest. is_in_set – True if the point lies in the uncertainty set, False otherwise. bool

Notes

This method is invoked at the outset of a PyROS solver call to determine whether a user-specified nominal parameter realization lies in the uncertainty set.

property rhs_vec

Right-hand side vector for polyhedral constraints defining the budget set. This also includes entries for nonnegativity constraints on the uncertain parameters.

This attribute cannot be set, and is automatically determined given other attributes.

Type: (M + N,) numpy.ndarray
property type

Brief description of the type of the uncertainty set.

Type: str
class pyomo.contrib.pyros.uncertainty_sets.FactorModelSet(origin, number_of_factors, psi_mat, beta)[source]

A factor model (a.k.a “net-alpha” model) set.

Parameters: origin ((N,) array_like) – Uncertain parameter values around which deviations are restrained. number_of_factors (int) – Natural number representing the dimensionality of the space to which the set projects. psi_mat ((N, number_of_factors) array_like) – Matrix designating each uncertain parameter’s contribution to each factor. Each row is associated with a separate uncertain parameter. Each column is associated with a separate factor. beta (numeric type) – Real value between 0 and 1 specifying the fraction of the independent factors that can simultaneously attain their extreme values.
property beta

Real number ranging from 0 to 1 representing the fraction of the independent factors that can simultaneously attain their extreme values.

Note that mathematically, setting beta = 0 will enforce that as many factors will be above 0 as there will be below 0 (i.e., “zero-net-alpha” model). Setting beta = 1 produces the hyper-rectangle [origin - psi @ e, origin + psi @ e], where e is a vector of ones.

Type: numeric type
property dim

Dimension of the factor model set.

Type: int
property number_of_factors

Natural number representing the dimensionality of the space to which the set projects.

This attribute is immutable, and may only be set at object construction. Typically, the number of factors is significantly less than the set dimension, but no restriction to that end is imposed here.

Type: int
property origin

Uncertain parameter values around which deviations are restrained.

Type: (N,) numpy.ndarray
property parameter_bounds

Uncertain parameter value bounds for the factor model set.

Returns: parameter_bounds – A list of 2-tuples of numerical values. Each tuple specifies the uncertain parameter bounds for the corresponding set dimension. list of tuples
point_in_set(point)[source]

Determine whether a given point lies in the factor model set.

Parameters: point ((N,) array-like) – Point (parameter value) of interest. True if the point lies in the set, False otherwise. bool
property psi_mat

Matrix designating each uncertain parameter’s contribution to each factor. Each row is associated with a separate uncertain parameter. Each column with a separate factor.

Type: (N, number_of_factors) numpy.ndarray
property type

Brief description of the type of the uncertainty set.

Type: str
class pyomo.contrib.pyros.uncertainty_sets.PolyhedralSet(lhs_coefficients_mat, rhs_vec)[source]

A bounded convex polyhedron or polytope.

Parameters: lhs_coefficients_mat ((M, N) array_like) – Left-hand side coefficients for the linear inequality constraints defining the polyhedral set. rhs_vec ((M,) array_like) – Right-hand side values for the linear inequality constraints defining the polyhedral set.
property coefficients_mat

Coefficient matrix for the (linear) inequality constraints defining the polyhedral set.

In tandem with the rhs_vec attribute, this matrix should be such that the polyhedral set is nonempty and bounded. Such a check is performed only at instance construction.

Type: (M, N) numpy.ndarray
property dim

Dimension of the polyhedral set.

Type: int
property parameter_bounds

Uncertain parameter value bounds for the polyhedral set.

Currently, an empty list, as the bounds cannot, in general, be computed without access to an optimization solver.

point_in_set(point)

Determine whether a given point lies in the uncertainty set.

Parameters: point ((N,) array-like) – Point (parameter value) of interest. is_in_set – True if the point lies in the uncertainty set, False otherwise. bool

Notes

This method is invoked at the outset of a PyROS solver call to determine whether a user-specified nominal parameter realization lies in the uncertainty set.

property rhs_vec

Right-hand side values (upper bounds) for the (linear) inequality constraints defining the polyhedral set.

Type: (M,) numpy.ndarray
property type

Brief description of the type of the uncertainty set.

Type: str
class pyomo.contrib.pyros.uncertainty_sets.AxisAlignedEllipsoidalSet(center, half_lengths)[source]

An axis-aligned ellipsoid.

Parameters: center ((N,) array_like) – Center of the ellipsoid. half_lengths ((N,) aray_like) – Semi-axis lengths of the ellipsoid.
property center

Center of the ellipsoid.

Type: (N,) numpy.ndarray
property dim

Dimension of the axis-aligned ellipsoidal set.

Type: int
property half_lengths

Semi-axis lengths.

Type: (N,) numpy.ndarray
property parameter_bounds

Uncertain parameter value bounds for the axis-aligned ellipsoidal set.

Returns: parameter_bounds – A list of 2-tuples of numerical values. Each tuple specifies the uncertain parameter bounds for the corresponding set dimension. list of tuples
point_in_set(point)

Determine whether a given point lies in the uncertainty set.

Parameters: point ((N,) array-like) – Point (parameter value) of interest. is_in_set – True if the point lies in the uncertainty set, False otherwise. bool

Notes

This method is invoked at the outset of a PyROS solver call to determine whether a user-specified nominal parameter realization lies in the uncertainty set.

property type

Brief description of the type of the uncertainty set.

Type: str
class pyomo.contrib.pyros.uncertainty_sets.EllipsoidalSet(center, shape_matrix, scale=1)[source]

A general ellipsoid.

Parameters: center ((N,) array-like) – Center of the ellipsoid. shape_matrix ((N, N) array-like) – A positive definite matrix characterizing the shape and orientation of the ellipsoid. scale (numeric type, optional) – Square of the factor by which to scale the semi-axes of the ellipsoid (i.e. the eigenvectors of the shape matrix). The default is 1.
property center

Center of the ellipsoid.

Type: (N,) numpy.ndarray
property dim

Dimension of the ellipsoidal set.

Type: int
property parameter_bounds

Uncertain parameter value bounds for the ellipsoidal set.

Returns: parameter_bounds – A list of 2-tuples of numerical values. Each tuple specifies the uncertain parameter bounds for the corresponding set dimension. list of tuples
point_in_set(point)

Determine whether a given point lies in the uncertainty set.

Parameters: point ((N,) array-like) – Point (parameter value) of interest. is_in_set – True if the point lies in the uncertainty set, False otherwise. bool

Notes

This method is invoked at the outset of a PyROS solver call to determine whether a user-specified nominal parameter realization lies in the uncertainty set.

property scale

Square of the factor by which to scale the semi-axes of the ellipsoid (i.e. the eigenvectors of the shape matrix).

Type: numeric type
property shape_matrix

A positive definite matrix characterizing the shape and orientation of the ellipsoid.

Type: (N, N) numpy.ndarray
property type

Brief description of the type of the uncertainty set.

Type: str
class pyomo.contrib.pyros.uncertainty_sets.UncertaintySet[source]

An object representing an uncertainty set for a two-stage robust optimization model. Along with a ConcreteModel object representing the corresponding deterministic model formulation, the uncertainty set object may be passed to the PyROS solver to obtain a robust model solution.

An UncertaintySet object should be viewed as merely a container for data needed to parameterize the set it represents, such that the object’s attributes do not, in general, reference the components of a Pyomo modeling object.

abstract property dim

Dimension of the uncertainty set (number of uncertain parameters in a corresponding optimization model of interest).

abstract property parameter_bounds

Bounds for the value of each uncertain parameter constrained by the set (i.e. bounds for each set dimension).

point_in_set(point)[source]

Determine whether a given point lies in the uncertainty set.

Parameters: point ((N,) array-like) – Point (parameter value) of interest. is_in_set – True if the point lies in the uncertainty set, False otherwise. bool

Notes

This method is invoked at the outset of a PyROS solver call to determine whether a user-specified nominal parameter realization lies in the uncertainty set.

class pyomo.contrib.pyros.uncertainty_sets.DiscreteScenarioSet(scenarios)[source]

A discrete set of finitely many uncertain parameter realizations (or scenarios).

Parameters: scenarios ((M, N) array_like) – A sequence of M distinct uncertain parameter realizations.
property dim

Dimension of the discrete scenario set.

Type: int
property parameter_bounds

Uncertain parameter value bounds for the discrete scenario set.

Returns: parameter_bounds – A list of 2-tuples of numerical values. Each tuple specifies the uncertain parameter bounds for the corresponding set dimension. list of tuples
point_in_set(point)[source]

Determine whether a given point lies in the discrete scenario set.

Parameters: point ((N,) array-like) – Point (parameter value) of interest. True if the point lies in the set, False otherwise. bool
property scenarios

Uncertain parameter realizations comprising the set. Each tuple is an uncertain parameter realization.

Note that the scenarios attribute may be modified, but only such that the dimension of the set remains unchanged.

Type: list(tuple)
property type

Brief description of the type of the uncertainty set.

Type: str
class pyomo.contrib.pyros.uncertainty_sets.IntersectionSet(**unc_sets)[source]

An intersection of a sequence of uncertainty sets, each of which is represented by an UncertaintySet object.

Parameters: **uncertainty_sets (dict) – PyROS UncertaintySet objects of which to construct an intersection. At least two uncertainty sets must be provided. All sets must be of the same dimension.
property all_sets

List of the uncertainty sets of which to take the intersection. Must be of minimum length 2.

This attribute may be set through any iterable of UncertaintySet objects, and exhibits similar behavior to a list.

Type: UncertaintySetList
property dim

Dimension of the intersection set.

Type: int
property parameter_bounds

Uncertain parameter value bounds for the intersection set.

Currently, an empty list, as the bounds cannot, in general, be computed without access to an optimization solver.

point_in_set(point)[source]

Determine whether a given point lies in the intersection set.

Parameters: point ((N,) array-like) – Point (parameter value) of interest. True if the point lies in the set, False otherwise. bool
property type

Brief description of the type of the uncertainty set.

Type: str

## PyROS Usage Example

We will use an example to illustrate the usage of PyROS. The problem we will use is called hydro and comes from the GAMS example problem database in The GAMS Model Library. The model was converted to Pyomo format via the GAMS Convert tool.

This model is a QCQP with 31 variables. Of these variables, 13 represent degrees of freedom, with the additional 18 being state variables. The model features 6 linear inequality constraints, 6 linear equality constraints, 6 non-linear (quadratic) equalities, and a quadratic objective. We have augmented this model by converting one objective coefficient, two constraint coefficients, and one constraint right-hand side into Param objects so that they can be considered uncertain later on.

Note

Per our analysis, the hydro problem satisfies the requirement that each value of $$\left(x, z, q \right)$$ maps to a unique value of $$y$$, which indicates a proper partition of variables between (first- or second-stage) degrees of freedom and state variables.

### Step 0: Import Pyomo and the PyROS Module

In anticipation of using the PyROS solver and building the deterministic Pyomo model:

>>> # === Required import ===
>>> import pyomo.environ as pyo
>>> import pyomo.contrib.pyros as pyros

>>> # === Instantiate the PyROS solver object ===
>>> pyros_solver = pyo.SolverFactory("pyros")


### Step 1: Define the Deterministic Problem

The deterministic Pyomo model for hydro is shown below.

Note

Primitive data (Python literals) that have been hard-coded within a deterministic model cannot be later considered uncertain, unless they are first converted to Param objects within the ConcreteModel object. Furthermore, any Param object that is to be later considered uncertain must have the property mutable=True.

Note

In case modifying the mutable property inside the deterministic model object itself is not straight-forward in your context, you may consider adding the following statement after import pyomo.environ as pyo but before defining the model object: pyo.Param.DefaultMutable = True. Note how this sets the default mutable property in all Param objects in the ensuing model instance to True; consequently, this solution will not work with Param objects for which the mutable=False property was explicitly enabled inside the model object.

>>> # === Construct the Pyomo model object ===
>>> m = pyo.ConcreteModel()
>>> m.name = "hydro"

>>> # === Define variables ===
>>> m.x1 = pyo.Var(within=pyo.Reals,bounds=(150,1500),initialize=150)
>>> m.x2 = pyo.Var(within=pyo.Reals,bounds=(150,1500),initialize=150)
>>> m.x3 = pyo.Var(within=pyo.Reals,bounds=(150,1500),initialize=150)
>>> m.x4 = pyo.Var(within=pyo.Reals,bounds=(150,1500),initialize=150)
>>> m.x5 = pyo.Var(within=pyo.Reals,bounds=(150,1500),initialize=150)
>>> m.x6 = pyo.Var(within=pyo.Reals,bounds=(150,1500),initialize=150)
>>> m.x7 = pyo.Var(within=pyo.Reals,bounds=(0,1000),initialize=0)
>>> m.x8 = pyo.Var(within=pyo.Reals,bounds=(0,1000),initialize=0)
>>> m.x9 = pyo.Var(within=pyo.Reals,bounds=(0,1000),initialize=0)
>>> m.x10 = pyo.Var(within=pyo.Reals,bounds=(0,1000),initialize=0)
>>> m.x11 = pyo.Var(within=pyo.Reals,bounds=(0,1000),initialize=0)
>>> m.x12 = pyo.Var(within=pyo.Reals,bounds=(0,1000),initialize=0)
>>> m.x13 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
>>> m.x14 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
>>> m.x15 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
>>> m.x16 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
>>> m.x17 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
>>> m.x18 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
>>> m.x19 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
>>> m.x20 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
>>> m.x21 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
>>> m.x22 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
>>> m.x23 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
>>> m.x24 = pyo.Var(within=pyo.Reals,bounds=(0,None),initialize=0)
>>> m.x25 = pyo.Var(within=pyo.Reals,bounds=(100000,100000),initialize=100000)
>>> m.x26 = pyo.Var(within=pyo.Reals,bounds=(60000,120000),initialize=60000)
>>> m.x27 = pyo.Var(within=pyo.Reals,bounds=(60000,120000),initialize=60000)
>>> m.x28 = pyo.Var(within=pyo.Reals,bounds=(60000,120000),initialize=60000)
>>> m.x29 = pyo.Var(within=pyo.Reals,bounds=(60000,120000),initialize=60000)
>>> m.x30 = pyo.Var(within=pyo.Reals,bounds=(60000,120000),initialize=60000)
>>> m.x31 = pyo.Var(within=pyo.Reals,bounds=(60000,120000),initialize=60000)

>>> # === Define parameters ===
>>> m.set_of_params = pyo.Set(initialize=[0, 1, 2, 3])
>>> nominal_values = {0:82.8*0.0016, 1:4.97, 2:4.97, 3:1800}
>>> m.p = pyo.Param(m.set_of_params, initialize=nominal_values, mutable=True)

>>> # === Specify the objective function ===
>>> m.obj = pyo.Objective(expr=m.p[0]*m.x1**2 + 82.8*8*m.x1 + 82.8*0.0016*m.x2**2 +
...                               82.8*82.8*8*m.x2 + 82.8*0.0016*m.x3**2 + 82.8*8*m.x3 +
...                               82.8*0.0016*m.x4**2 + 82.8*8*m.x4 + 82.8*0.0016*m.x5**2 +
...                               82.8*8*m.x5 + 82.8*0.0016*m.x6**2 + 82.8*8*m.x6 + 248400,
...                               sense=pyo.minimize)

>>> # === Specify the constraints ===
>>> m.c2 = pyo.Constraint(expr=-m.x1 - m.x7 + m.x13 + 1200<= 0)
>>> m.c3 = pyo.Constraint(expr=-m.x2 - m.x8 + m.x14 + 1500 <= 0)
>>> m.c4 = pyo.Constraint(expr=-m.x3 - m.x9 + m.x15 + 1100 <= 0)
>>> m.c5 = pyo.Constraint(expr=-m.x4 - m.x10 + m.x16 + m.p[3] <= 0)
>>> m.c6 = pyo.Constraint(expr=-m.x5 - m.x11 + m.x17 + 950 <= 0)
>>> m.c7 = pyo.Constraint(expr=-m.x6 - m.x12 + m.x18 + 1300 <= 0)
>>> m.c8 = pyo.Constraint(expr=12*m.x19 - m.x25 + m.x26 == 24000)
>>> m.c9 = pyo.Constraint(expr=12*m.x20 - m.x26 + m.x27 == 24000)
>>> m.c10 = pyo.Constraint(expr=12*m.x21 - m.x27 + m.x28 == 24000)
>>> m.c11 = pyo.Constraint(expr=12*m.x22 - m.x28 + m.x29 == 24000)
>>> m.c12 = pyo.Constraint(expr=12*m.x23 - m.x29 + m.x30 == 24000)
>>> m.c13 = pyo.Constraint(expr=12*m.x24 - m.x30 + m.x31 == 24000)
>>> m.c14 = pyo.Constraint(expr=-8e-5*m.x7**2 + m.x13 == 0)
>>> m.c15 = pyo.Constraint(expr=-8e-5*m.x8**2 + m.x14 == 0)
>>> m.c16 = pyo.Constraint(expr=-8e-5*m.x9**2 + m.x15 == 0)
>>> m.c17 = pyo.Constraint(expr=-8e-5*m.x10**2 + m.x16 == 0)
>>> m.c18 = pyo.Constraint(expr=-8e-5*m.x11**2 + m.x17 == 0)
>>> m.c19 = pyo.Constraint(expr=-8e-5*m.x12**2 + m.x18 == 0)
>>> m.c20 = pyo.Constraint(expr=-4.97*m.x7 + m.x19 == 330)
>>> m.c21 = pyo.Constraint(expr=-m.p[1]*m.x8 + m.x20 == 330)
>>> m.c22 = pyo.Constraint(expr=-4.97*m.x9 + m.x21 == 330)
>>> m.c23 = pyo.Constraint(expr=-4.97*m.x10 + m.x22 == 330)
>>> m.c24 = pyo.Constraint(expr=-m.p[2]*m.x11 + m.x23 == 330)
>>> m.c25 = pyo.Constraint(expr=-4.97*m.x12 + m.x24 == 330)


### Step 2: Define the Uncertainty

First, we need to collect into a list those Param objects of our model that represent potentially uncertain parameters. For purposes of our example, we shall assume uncertainty in the model parameters (m.p[0], m.p[1], m.p[2], m.p[3]), for which we can conveniently utilize the m.p object (itself an indexed Param object).

>>> # === Specify which parameters are uncertain ===
>>> uncertain_parameters = [m.p] # We can pass IndexedParams this way to PyROS, or as an expanded list per index


Note

Any Param object that is to be considered uncertain by PyROS must have the property mutable=True.

PyROS will seek to identify solutions that remain feasible for any realization of these parameters included in an uncertainty set. To that end, we need to construct an UncertaintySet object. In our example, let us utilize the BoxSet constructor to specify an uncertainty set of simple hyper-rectangular geometry. For this, we will assume each parameter value is uncertain within a percentage of its nominal value. Constructing this specific UncertaintySet object can be done as follows.

>>> # === Define the pertinent data ===
>>> relative_deviation = 0.15
>>> bounds = [(nominal_values[i] - relative_deviation*nominal_values[i],
...               nominal_values[i] + relative_deviation*nominal_values[i])
...               for i in range(4)]

>>> # === Construct the desirable uncertainty set ===
>>> box_uncertainty_set = pyros.BoxSet(bounds=bounds)


### Step 3: Solve with PyROS

PyROS requires the user to supply one local and one global NLP solver to be used for solving sub-problems. For convenience, we shall have PyROS invoke BARON as both the local and the global NLP solver.

>>> # === Designate local and global NLP solvers ===
>>> local_solver = pyo.SolverFactory('baron')
>>> global_solver = pyo.SolverFactory('baron')


Note

Additional solvers to be used as backup can be designated during the solve statement via the config options backup_local_solvers and backup_global_solvers presented above.

The final step in solving a model with PyROS is to designate the remaining required inputs, namely first_stage_variables and second_stage_variables. Below, we present two separate cases.

#### PyROS Termination Conditions

PyROS will return one of six termination conditions upon completion. These termination conditions are tabulated below.

 Termination Condition Description pyrosTerminationCondition.robust_optimal The final solution is robust optimal pyrosTerminationCondition.robust_feasible The final solution is robust feasible pyrosTerminationCondition.robust_infeasible The posed problem is robust infeasible pyrosTerminationCondition.max_iter Maximum number of GRCS iteration reached pyrosTerminationCondition.time_out Maximum number of time reached pyrosTerminationCondition.subsolver_error Unacceptable return status(es) from a user-supplied sub-solver

#### A Single-Stage Problem

If we choose to designate all variables as either design or state variables, without any control variables (i.e., all degrees of freedom are first-stage), we can use PyROS to solve the single-stage problem as shown below. In particular, let us instruct PyROS that variables m.x1 through m.x6, m.x19 through m.x24, and m.x31 correspond to first-stage degrees of freedom.

>>> # === Designate which variables correspond to first- and second-stage degrees of freedom ===
>>> first_stage_variables =[m.x1, m.x2, m.x3, m.x4, m.x5, m.x6,
...                      m.x19, m.x20, m.x21, m.x22, m.x23, m.x24, m.x31]
>>> second_stage_variables = []
>>> # The remaining variables are implicitly designated to be state variables

>>> # === Call PyROS to solve the robust optimization problem ===
>>> results_1 = pyros_solver.solve(model = m,
...                                  first_stage_variables = first_stage_variables,
...                                  second_stage_variables = second_stage_variables,
...                                  uncertain_params = uncertain_parameters,
...                                  uncertainty_set = box_uncertainty_set,
...                                  local_solver = local_solver,
...                                  global_solver= global_solver,
...                                  options = {
...                                     "objective_focus": pyros.ObjectiveType.worst_case,
...                                     "solve_master_globally": True,
...                                   })
===========================================================================================
PyROS: Pyomo Robust Optimization Solver ...
===========================================================================================
...
INFO: Robust optimal solution identified. Exiting PyROS.

>>> # === Query results ===
>>> time = results_1.time
>>> iterations = results_1.iterations
>>> termination_condition = results_1.pyros_termination_condition
>>> objective = results_1.final_objective_value
>>> # === Print some results ===
>>> single_stage_final_objective = round(objective,-1)
>>> print("Final objective value: %s" % single_stage_final_objective)
Final objective value: 48367380.0
>>> print("PyROS termination condition: %s" % termination_condition)
PyROS termination condition: pyrosTerminationCondition.robust_optimal


#### PyROS Results Object

The results object returned by PyROS allows you to query the following information from the solve call: total iterations of the algorithm iterations, CPU time time, the GRCS algorithm termination condition pyros_termination_condition, and the final objective function value final_objective_value. If the option load_solution = True (default), the variables in the model will be loaded to the solution determined by PyROS and can be obtained by querying the model variables. Note that in the results obtained above, we set load_solution = False. This is to ensure that the next set of runs shown here can utilize the original deterministic model, as the initial point can affect the performance of sub-solvers.

Note

The reported final_objective_value and final model variable values depend on the selection of the option objective_focus. The final_objective_value is the sum of first-stage and second-stage objective functions. If objective_focus = ObjectiveType.nominal, second-stage objective and variables are evaluated at the nominal realization of the uncertain parameters, $$q^0$$. If objective_focus = ObjectiveType.worst_case, second-stage objective and variables are evaluated at the worst-case realization of the uncertain parameters, $$q^{k^\ast}$$ where $$k^\ast = argmax_{k \in \mathcal{K}} f_2(x,z^k,y^k,q^k)$$ .

An example of how to query these values on the previously obtained results is shown in the code above.

#### A Two-Stage Problem

For this next set of runs, we will assume that some of the previously designated first-stage degrees of freedom are in fact second-stage ones. PyROS handles second-stage degrees of freedom via the use of decision rules, which is controlled with the config option decision_rule_order presented above. Here, we shall select affine decision rules by setting decision_rule_order to the value of 1.

>>> # === Define the variable partitioning
>>> first_stage_variables =[m.x5, m.x6, m.x19, m.x22, m.x23, m.x24, m.x31]
>>> second_stage_variables = [m.x1, m.x2, m.x3, m.x4, m.x20, m.x21]
>>> # The remaining variables are implicitly designated to be state variables

>>> # === Call PyROS to solve the robust optimization problem ===
>>> results_2 = pyros_solver.solve(model = m,
...                                   first_stage_variables = first_stage_variables,
...                                   second_stage_variables = second_stage_variables,
...                                   uncertain_params = uncertain_parameters,
...                                   uncertainty_set = box_uncertainty_set,
...                                   local_solver = local_solver,
...                                   global_solver = global_solver,
...                                   options = {
...                                      "objective_focus": pyros.ObjectiveType.worst_case,
...                                      "solve_master_globally": True,
...                                      "decision_rule_order": 1
...                                   })
===========================================================================================
PyROS: Pyomo Robust Optimization Solver ...
...
INFO: Robust optimal solution identified. Exiting PyROS.

>>> # === Compare final objective to the singe-stage solution
>>> two_stage_final_objective = round(pyo.value(results_2.final_objective_value),-1)
>>> percent_difference = 100 * (two_stage_final_objective - single_stage_final_objective)/(single_stage_final_objective)
>>> print("Percent objective change relative to constant decision rules objective: %.2f %%" % percent_difference)
Percent objective change relative to constant decision rules objective: -24...


In this example, when we compare the final objective value in the case of constant decision rules (no second-stage recourse) and affine decision rules, we see there is a ~25% decrease in total objective value.

#### The Price of Robustness

Using appropriately constructed hierarchies, PyROS allows for the facile comparison of robust optimal objectives across sets to determine the “price of robustness.” For the set we considered here, the BoxSet, we can create such a hierarchy via an array of relative_deviation parameters to define the size of these uncertainty sets. We can then loop through this array and call PyROS within a loop to identify robust solutions in light of each of the specified BoxSet objects.

>>> # This takes a long time to run and therefore is not a doctest
>>> # === An array of maximum relative deviations from the nominal uncertain parameter values to utilize in constructing box sets
>>> relative_deviation_list = [0.00, 0.10, 0.20, 0.30, 0.40]
>>> # === Final robust optimal objectives
>>> robust_optimal_objectives = []
>>> for relative_deviation in relative_deviation_list:
...   bounds = [(nominal_values[i] - relative_deviation*nominal_values[i],
...                   nominal_values[i] + relative_deviation*nominal_values[i])
...                   for i in range(4)]
...   box_uncertainty_set = pyros.BoxSet(bounds = bounds)
...   results = pyros_solver.solve(model = m,
...                                     first_stage_variables = first_stage_variables,
...                                     second_stage_variables = second_stage_variables,
...                                     uncertain_params = uncertain_parameters,
...                                     uncertainty_set = box_uncertainty_set,
...                                     local_solver = local_solver,
...                                     global_solver = global_solver,
...                                     options = {
...                                        "objective_focus": pyros.ObjectiveType.worst_case,
...                                        "solve_master_globally": True,
...                                        "decision_rule_order": 1
...                                     })
...   if results.pyros_termination_condition != pyros.pyrosTerminationCondition.robust_optimal:
...       print("This instance didn't solve to robust optimality.")
...       robust_optimal_objective.append("-----")
...   else:
...       robust_optimal_objectives.append(str(results.final_objective_value))


For this example, we obtain the following price of robustness results:

 Uncertainty Set Size (+/-) o Robust Optimal Objective % Increase x 0.00 35,837,659.18 0.00 % 0.10 36,135,191.59 0.82 % 0.20 36,437,979.81 1.64 % 0.30 43,478,190.92 17.57 % 0.40 robust_infeasible $$\text{-----}$$

Note how, in the case of the last uncertainty set, we were able to utilize PyROS to show the robust infeasibility of the problem.

o Relative Deviation from Nominal Realization

x Relative to Deterministic Optimal Objective

This clearly illustrates the impact that the uncertainty set size can have on the robust optimal objective values. Price of robustness studies like this are easily implemented using PyROS.

Warning

PyROS is still under a beta release. Please provide feedback and/or report any problems by opening an issue on the Pyomo GitHub page.