PyROS Solver

PyROS (Pyomo Robust Optimization Solver) is a Pyomo-based meta-solver for non-convex, two-stage adjustable robust optimization problems.

It was developed by Natalie M. Isenberg, Jason A. F. Sherman, 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.

PyROS is designed to operate on deterministic models of the general form

\[\begin{split}\begin{array}{clll} \displaystyle \min_{\substack{x \in \mathcal{X}, \\ z \in \mathbb{R}^{n_z}, y\in\mathbb{R}^{n_y}}} & ~~ f_1\left(x\right) + f_2(x,z,y; q^{\text{nom}}) & \\ \displaystyle \text{s.t.} & ~~ g_i(x, z, y; q^{\text{nom}}) \leq 0 & \forall\,i \in \mathcal{I} \\ & ~~ h_j(x,z,y; q^{\text{nom}}) = 0 & \forall\,j \in \mathcal{J} \\ \end{array}\end{split}\]

where:

  • \(x \in \mathcal{X}\) are the “design” variables (i.e., first-stage degrees of freedom), where \(\mathcal{X} \subseteq \mathbb{R}^{n_x}\) is the feasible space defined by the model constraints (including variable bounds specifications) referencing \(x\) only.
  • \(z \in \mathbb{R}^{n_z}\) are the “control” variables (i.e., second-stage degrees of freedom)
  • \(y \in \mathbb{R}^{n_y}\) are the “state” variables
  • \(q \in \mathbb{R}^{n_q}\) is the vector of model parameters considered uncertain, and \(q^{\text{nom}}\) 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 all variables and the uncertain parameters
  • \(g_i\left(x, z, y; q\right)\) is the \(i^\text{th}\) inequality constraint function in set \(\mathcal{I}\) (see Note)
  • \(h_j\left(x, z, y; q\right)\) is the \(j^\text{th}\) equality constraint function in set \(\mathcal{J}\) (see Note)

Note

PyROS accepts models in which bounds are directly imposed on Var objects representing components of the variables \(z\) and \(y\). These models are cast to the form above by reformulating the bounds as inequality constraints.

Note

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 the mapping is not unique, 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 of the determinstic model, we now assume that the uncertain parameters may attain any realization in a compact uncertainty set \(\mathcal{Q} \subseteq \mathbb{R}^{n_q}\) containing the nominal value \(q^{\text{nom}}\). The set \(\mathcal{Q}\) may be either continuous or discrete.

Based on the above notation, the form of the robust counterpart addressed by PyROS is

\[\begin{split}\begin{array}{ccclll} \displaystyle \min_{x \in \mathcal{X}} & \displaystyle \max_{q \in \mathcal{Q}} & \displaystyle \min_{\substack{z \in \mathbb{R}^{n_z},\\y \in \mathbb{R}^{n_y}}} \ \ & \displaystyle ~~ f_1\left(x\right) + f_2\left(x, z, y, q\right) \\ & & \text{s.t.}~ & \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{array}\end{split}\]

PyROS solves problems of this form using the Generalized Robust Cutting-Set algorithm developed in [Isenberg_et_al].

When using PyROS, please consider citing the above paper.

PyROS Required Inputs

The required inputs to the PyROS solver are:

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

These are more elaborately presented in the Solver Interface section.

Note

Any variables in the model not specified to be first-stage 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 a model.

Parameters:
  • model (ConcreteModel) – The deterministic model.
  • first_stage_variables (list of Var) – First-stage model variables (or design variables).
  • second_stage_variables (list of Var) – Second-stage model variables (or control variables).
  • uncertain_params (list of Param) – Uncertain model parameters. The mutable attribute for every uncertain parameter objects must be set to True.
  • uncertainty_set (UncertaintySet) – Uncertainty set against which the solution(s) returned will be confirmed to be robust.
  • local_solver (Solver) – Subordinate local NLP solver.
  • global_solver (Solver) – Subordinate global NLP solver.
Returns:

return_soln – Summary of PyROS termination outcome.
Return type:

ROSolveResults
Keyword Arguments:
 
  • time_limit (None or NonNegativeFloat, optional) – Wall time limit for the execution of the PyROS solver in seconds (including time spent by subsolvers). If None is provided, then no time limit is enforced.
  • keepfiles (bool, default=False) – Export subproblems with a non-acceptable termination status for debugging purposes. If True is provided, then the argument subproblem_file_directory must also be specified.
  • tee (bool, default=False) – Output subordinate solver logs for all subproblems.
  • load_solution (bool, default=True) – Load final solution(s) found by PyROS to the deterministic model provided.
  • objective_focus (ObjectiveType, optional) –

    Objective focus for the master problems:

    • ObjectiveType.nominal: Optimize the objective function subject to the nominal uncertain parameter realization.
    • ObjectiveType.worst_case: Optimize the objective function subject to the worst-case uncertain parameter realization.

    By default, ObjectiveType.nominal is chosen.

    A worst-case objective focus is required for certification of robust optimality of the final solution(s) returned by PyROS. If a nominal objective focus is chosen, then only robust feasibility is guaranteed.

  • nominal_uncertain_param_vals (list of float, default=[]) – Nominal uncertain parameter realization. Entries should be provided in an order consistent with the entries of the argument uncertain_params. If an empty list is provided, then the values of the Param objects specified through uncertain_params are chosen.
  • decision_rule_order (In[0, 1, 2], default=0) –

    Order (or degree) of the polynomial decision rule functions used for approximating the adjustability of the second stage variables with respect to the uncertain parameters.

    Choices are:

    • 0: static recourse
    • 1: affine recourse
    • 2: quadratic recourse
  • solve_master_globally (bool, default=False) – True to solve all master problems with the subordinate global solver, False to solve all master problems with the subordinate local solver. Along with a worst-case objective focus (see argument objective_focus), solving the master problems to global optimality is required for certification of robust optimality of the final solution(s) returned by PyROS. Otherwise, only robust feasibility is guaranteed.
  • max_iter (int, default=-1) – Iteration limit. If -1 is provided, then no iteration limit is enforced.
  • robust_feasibility_tolerance (NonNegativeFloat, default=0.0001) – Relative tolerance for assessing maximal inequality constraint violations during the GRCS separation step.
  • separation_priority_order (dict, default={}) – Mapping from model inequality constraint names to positive integers specifying the priorities of their corresponding separation subproblems. A higher integer value indicates a higher priority. Constraints not referenced in the dict assume a priority of 0. Separation subproblems are solved in order of decreasing priority.
  • progress_logger (str or logging.Logger, default='pyomo.contrib.pyros') – Logger (or name thereof) used for reporting PyROS solver progress. If a str is specified, then logging.getLogger(progress_logger) is used.
  • backup_local_solvers (list of Solver, default=[]) – Additional subordinate local NLP optimizers to invoke in the event the primary local NLP optimizer fails to solve a subproblem to an acceptable termination condition.
  • backup_global_solvers (list of Solver, default=[]) – Additional subordinate global NLP optimizers to invoke in the event the primary global NLP optimizer fails to solve a subproblem to an acceptable termination condition.
  • subproblem_file_directory (None, str, or path-like, default=None) – Directory to which to export subproblems not successfully solved to an acceptable termination condition. In the event keepfiles=True is specified, a str or path-like referring to an existing directory must be provided.
  • bypass_local_separation (bool, default=False) – This is an advanced option. Solve all separation subproblems with the subordinate global solver(s) only. This option is useful for expediting PyROS in the event that the subordinate global optimizer(s) provided can quickly solve separation subproblems to global optimality.
  • bypass_global_separation (bool, default=False) – This is an advanced option. Solve all separation subproblems with the subordinate local solver(s) only. If True is chosen, then robustness of the final solution(s) returned by PyROS is not guaranteed, and a warning will be issued at termination. This option is useful for expediting PyROS in the event that the subordinate global optimizer provided cannot tractably solve separation subproblems to global optimality.
  • p_robustness (dict, default={}) – This is an advanced option. Add p-robustness constraints to all master subproblems. If an empty dict is provided, then p-robustness constraints are not added. Otherwise, the dict must map a str of value 'rho' to a non-negative float. PyROS automatically specifies 1 + p_robustness['rho'] as an upper bound for the ratio of the objective function value under any PyROS-sampled uncertain parameter realization to the objective function under the nominal parameter realization.

Note

Upon successful convergence of PyROS, the solution returned is certified to be robust optimal only if:

  1. master problems are solved to global optimality (by specifying solve_master_globally=True)
  2. a worst-case objective focus is chosen (by specifying objective_focus=ObjectiveType.worst_case)

Otherwise, the solution returned is certified to only be robust feasible.

PyROS Uncertainty Sets

Uncertainty sets are represented by subclasses of the UncertaintySet abstract base class. PyROS provides a suite of pre-implemented subclasses representing commonly used uncertainty sets. Custom user-defined uncertainty set types may be implemented by subclassing the UncertaintySet class. The intersection of a sequence of concrete UncertaintySet instances can be easily constructed by instantiating the pre-implemented IntersectionSet subclass.

The table that follows provides mathematical definitions of the various abstract and pre-implemented UncertaintySet subclasses.

Table 2 Mathematical definitions of PyROS uncertainty sets of dimension \(n\).
Uncertainty Set Type Input Data Mathematical Definition
BoxSet \(\begin{array}{l} q ^{\text{L}} \in \mathbb{R}^{n}, \\ q^{\text{U}} \in \mathbb{R}^{n} \end{array}\) \(\{q \in \mathbb{R}^n \mid q^\mathrm{L} \leq q \leq q^\mathrm{U}\}\)
CardinalitySet \(\begin{array}{l} q^{0} \in \mathbb{R}^{n}, \\ \hat{q} \in \mathbb{R}_{+}^{n}, \\ \Gamma \in [0, n] \end{array}\) \(\left\{ q \in \mathbb{R}^{n} \middle| \begin{array}{l} q = q^{0} + \hat{q} \circ \xi \\ \displaystyle \sum_{i=1}^{n} \xi_{i} \leq \Gamma \\ \xi \in [0, 1]^{n} \end{array} \right\}\)
BudgetSet \(\begin{array}{l} q^{0} \in \mathbb{R}^{n}, \\ b \in \mathbb{R}_{+}^{L}, \\ B \in \{0, 1\}^{L \times n} \end{array}\) \(\left\{ q \in \mathbb{R}^{n} \middle| \begin{array}{l} \begin{pmatrix} B \\ -I \end{pmatrix} q \leq \begin{pmatrix} b + Bq^{0} \\ -q^{0} \end{pmatrix} \end{array} \right\}\)
FactorModelSet \(\begin{array}{l} q^{0} \in \mathbb{R}^{n}, \\ \Psi \in \mathbb{R}^{n \times F}, \\ \beta \in [0, 1] \end{array}\) \(\left\{ q \in \mathbb{R}^{n} \middle| \begin{array}{l} q = q^{0} + \Psi \xi \\ \displaystyle\bigg| \sum_{j=1}^{F} \xi_{j} \bigg| \leq \beta F \\ \xi \in [-1, 1]^{F} \\ \end{array} \right\}\)
PolyhedralSet \(\begin{array}{l} A \in \mathbb{R}^{m \times n}, \\ b \in \mathbb{R}^{m}\end{array}\) \(\{q \in \mathbb{R}^{n} \mid A q \leq b\}\)
AxisAlignedEllipsoidalSet \(\begin{array}{l} q^0 \in \mathbb{R}^{n}, \\ \alpha \in \mathbb{R}_{+}^{n} \end{array}\) \(\left\{ q \in \mathbb{R}^{n} \middle| \begin{array}{l} \displaystyle\sum_{\substack{i = 1: \\ \alpha_{i} > 0}}^{n} \left(\frac{q_{i} - q_{i}^{0}}{\alpha_{i}}\right)^2 \leq 1 \\ q_{i} = q_{i}^{0} \,\forall\,i : \alpha_{i} = 0 \end{array} \right\}\)
EllipsoidalSet \(\begin{array}{l} q^0 \in \mathbb{R}^n, \\ P \in \mathbb{S}_{++}^{n}, \\ s \in \mathbb{R}_{+} \end{array}\) \(\{q \in \mathbb{R}^{n} \mid (q - q^{0})^{\intercal} P^{-1} (q - q^{0}) \leq s\}\)
UncertaintySet \(g: \mathbb{R}^{n} \to \mathbb{R}^{m}\) \(\{q \in \mathbb{R}^{n} \mid g(q) \leq 0\}\)
DiscreteScenarioSet \(q^{1}, q^{2},\dots , q^{S} \in \mathbb{R}^{n}\) \(\{q^{1}, q^{2}, \dots , q^{S}\}\)
IntersectionSet \(\mathcal{Q}_{1}, \mathcal{Q}_{2}, \dots , \mathcal{Q}_{m} \subset \mathbb{R}^{n}\) \(\displaystyle \bigcap_{i=1}^{m} \mathcal{Q}_{i}\)

Note

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

PyROS Uncertainty Set Classes

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

Bases: UncertaintySet

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

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

Examples

1D box set (interval):

>>> from pyomo.contrib.pyros import BoxSet
>>> interval = BoxSet(bounds=[(1, 2)])
>>> interval.bounds
array([[1, 2]])

2D box set:

>>> box_set = BoxSet(bounds=[[1, 2], [3, 4]])
>>> box_set.bounds
array([[1, 2],
       [3, 4]])

5D hypercube with bounds 0 and 1 in each dimension:

>>> hypercube_5d = BoxSet(bounds=[[0, 1] for idx in range(5)])
>>> hypercube_5d.bounds
array([[0, 1],
       [0, 1],
       [0, 1],
       [0, 1],
       [0, 1]])
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 N of the box set.

Type:int
property parameter_bounds

Bounds in each dimension of the box set. This is numerically equivalent to the bounds attribute.

Returns:List, length N, of 2-tuples. Each tuple specifies the bounds in its corresponding dimension.
Return type: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.
Returns:is_in_set – True if the point lies in the uncertainty set, False otherwise.
Return type: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]

Bases: UncertaintySet

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.

Examples

A 3D cardinality set:

>>> from pyomo.contrib.pyros import CardinalitySet
>>> gamma_set = CardinalitySet(
...     origin=[0, 0, 0],
...     positive_deviation=[1.0, 2.0, 1.5],
...     gamma=1,
... )
>>> gamma_set.origin
array([0, 0, 0])
>>> gamma_set.positive_deviation
array([1. , 2. , 1.5])
>>> gamma_set.gamma
1
property dim

Dimension N 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 N.

Note that, mathematically, setting gamma to 0 reduces the set to a singleton containing the center, while setting gamma to the set dimension N makes the set mathematically equivalent to a BoxSet with bounds numpy.array([origin, origin + positive_deviation]).T.

Type:numeric type
property origin

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

Type:(N,) numpy.ndarray
property parameter_bounds

Bounds in each dimension of the cardinality set.

Returns:List, length N, of 2-tuples. Each tuple specifies the bounds in its corresponding dimension.
Return type: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.
Returns:True if the point lies in the set, False otherwise.
Return type: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]

Bases: UncertaintySet

A budget set.

Parameters:
  • budget_membership_mat ((L, 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 ((L,) 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.

Examples

3D budget set with one budget constraint and no origin chosen (hence origin defaults to 3D zero vector):

>>> from pyomo.contrib.pyros import BudgetSet
>>> budget_set = BudgetSet(
...     budget_membership_mat=[[1, 1, 1]],
...     rhs_vec=[2],
... )
>>> budget_set.budget_membership_mat
array([[1, 1, 1]])
>>> budget_set.budget_rhs_vec
array([2])
>>> budget_set.origin
array([0., 0., 0.])

3D budget set with two budget constraints and custom origin:

>>> budget_custom = BudgetSet(
...     budget_membership_mat=[[1, 0, 1], [0, 1, 0]],
...     rhs_vec=[1, 1],
...     origin=[2, 2, 2],
... )
>>> budget_custom.budget_membership_mat
array([[1, 0, 1],
       [0, 1, 0]])
>>> budget_custom.budget_rhs_vec
array([1, 1])
>>> budget_custom.origin
array([2, 2, 2])
property budget_membership_mat

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

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

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

Type:(L,) 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:(L + N, N) numpy.ndarray
property dim

Dimension N of the budget set.

Type:int
property origin

Origin of the budget set.

Type:(N,) numpy.ndarray
property parameter_bounds

Bounds in each dimension of the budget set.

Returns:List, length N, of 2-tuples. Each tuple specifies the bounds in its corresponding dimension.
Return type: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.
Returns:is_in_set – True if the point lies in the uncertainty set, False otherwise.
Return type: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:(L + 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]

Bases: UncertaintySet

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, F) 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. Number of columns F of psi_mat should be equal to number_of_factors.
  • beta (numeric type) – Real value between 0 and 1 specifying the fraction of the independent factors that can simultaneously attain their extreme values.

Examples

A 4D factor model set with a 2D factor space:

>>> from pyomo.contrib.pyros import FactorModelSet
>>> import numpy as np
>>> fset = FactorModelSet(
...     origin=np.zeros(4),
...     number_of_factors=2,
...     psi_mat=np.full(shape=(4, 2), fill_value=0.1),
...     beta=0.5,
... )
>>> fset.origin
array([0., 0., 0., 0.])
>>> fset.number_of_factors
2
>>> fset.psi_mat
array([[0.1, 0.1],
       [0.1, 0.1],
       [0.1, 0.1],
       [0.1, 0.1]])
>>> fset.beta
0.5
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). If beta = 1, then the set is numerically equivalent to a BoxSet with bounds [origin - psi @ np.ones(F), origin + psi @ np.ones(F)].T.

Type:numeric type
property dim

Dimension N of the factor model set.

Type:int
property number_of_factors

Natural number representing the dimensionality F 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

Bounds in each dimension of the factor model set.

Returns:List, length N, of 2-tuples. Each tuple specifies the bounds in its corresponding dimension.
Return type: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.
Returns:True if the point lies in the set, False otherwise.
Return type: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, F) 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]

Bases: UncertaintySet

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. Each entry is an upper bound for the quantity lhs_coefficients_mat @ x, where x is an (N,) array representing any point in the polyhedral set.

Examples

2D polyhedral set with 4 defining inequalities:

>>> from pyomo.contrib.pyros import PolyhedralSet
>>> pset = PolyhedralSet(
...     lhs_coefficients_mat=[[-1, 0], [0, -1], [-1, 1], [1, 0]],
...     rhs_vec=[0, 0, 0, 1],
... )
>>> pset.coefficients_mat
array([[-1,  0],
       [ 0, -1],
       [-1,  1],
       [ 1,  0]])
>>> pset.rhs_vec
array([0, 0, 0, 1])
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 N of the polyhedral set.

Type:int
property parameter_bounds

Bounds in each dimension of the polyhedral set.

Currently, an empty list is returned, 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.
Returns:is_in_set – True if the point lies in the uncertainty set, False otherwise.
Return type: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]

Bases: UncertaintySet

An axis-aligned ellipsoid.

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

Examples

3D origin-centered unit hypersphere:

>>> from pyomo.contrib.pyros import AxisAlignedEllipsoidalSet
>>> sphere = AxisAlignedEllipsoidalSet(
...     center=[0, 0, 0],
...     half_lengths=[1, 1, 1]
... )
>>> sphere.center
array([0, 0, 0])
>>> sphere.half_lengths
array([1, 1, 1])
property center

Center of the ellipsoid.

Type:(N,) numpy.ndarray
property dim

Dimension N of the axis-aligned ellipsoidal set.

Type:int
property half_lengths

Semi-axis lengths.

Type:(N,) numpy.ndarray
property parameter_bounds

Bounds in each dimension of the axis-aligned ellipsoidal set.

Returns:List, length N, of 2-tuples. Each tuple specifies the bounds in its corresponding dimension.
Return type: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.
Returns:is_in_set – True if the point lies in the uncertainty set, False otherwise.
Return type: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]

Bases: UncertaintySet

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.

Examples

3D origin-centered unit hypersphere:

>>> from pyomo.contrib.pyros import EllipsoidalSet
>>> import numpy as np
>>> hypersphere = EllipsoidalSet(
...     center=[0, 0, 0],
...     shape_matrix=np.eye(3),
...     scale=1,
... )
>>> hypersphere.center
array([0, 0, 0])
>>> hypersphere.shape_matrix
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])
>>> hypersphere.scale
1

A 2D ellipsoid with custom rotation and scaling:

>>> rotated_ellipsoid = EllipsoidalSet(
...     center=[1, 1],
...     shape_matrix=[[4, 2], [2, 4]],
...     scale=0.5,
... )
>>> rotated_ellipsoid.center
array([1, 1])
>>> rotated_ellipsoid.shape_matrix
array([[4, 2],
       [2, 4]])
>>> rotated_ellipsoid.scale
0.5
property center

Center of the ellipsoid.

Type:(N,) numpy.ndarray
property dim

Dimension N of the ellipsoidal set.

Type:int
property parameter_bounds

Bounds in each dimension of the ellipsoidal set.

Returns:List, length N, of 2-tuples. Each tuple specifies the bounds in its corresponding dimension.
Return type: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.
Returns:is_in_set – True if the point lies in the uncertainty set, False otherwise.
Return type: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]

Bases: object

An object representing an uncertainty set to be passed to the PyROS solver.

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 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.
Returns:is_in_set – True if the point lies in the uncertainty set, False otherwise.
Return type: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]

Bases: UncertaintySet

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.

Examples

2D set with three scenarios:

>>> from pyomo.contrib.pyros import DiscreteScenarioSet
>>> discrete_set = DiscreteScenarioSet(
...     scenarios=[[1, 1], [2, 1], [1, 2]],
... )
>>> discrete_set.scenarios
[(1, 1), (2, 1), (1, 2)]
property dim

Dimension N of the discrete scenario set.

Type:int
property parameter_bounds

Bounds in each dimension of the discrete scenario set.

Returns:List, length N, of 2-tuples. Each tuple specifies the bounds in its corresponding dimension.
Return type: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.
Returns:True if the point lies in the set, False otherwise.
Return type: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 of tuples
property type

Brief description of the type of the uncertainty set.

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

Bases: UncertaintySet

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

Parameters:**unc_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.

Examples

Intersection of origin-centered 2D box (square) and 2D hypersphere (circle):

>>> from pyomo.contrib.pyros import (
...     BoxSet, AxisAlignedEllipsoidalSet, IntersectionSet,
... )
>>> square = BoxSet(bounds=[[-1.5, 1.5], [-1.5, 1.5]])
>>> circle = AxisAlignedEllipsoidalSet(
...     center=[0, 0],
...     half_lengths=[2, 2],
... )
>>> # to construct intersection, pass sets as keyword arguments
>>> intersection = IntersectionSet(set1=square, set2=circle)
>>> intersection.all_sets
UncertaintySetList([...])
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.
Returns:True if the point lies in the set, False otherwise.
Return type:bool
property type

Brief description of the type of the uncertainty set.

Type:str

PyROS Usage Example

In this section, we illustrate the usage of PyROS with a modeling example. The deterministic problem of interest is called hydro (available here), a QCQP taken from the GAMS Model Library. We have converted the model to Pyomo format using the GAMS Convert tool.

The hydro model features 31 variables, of which 13 are degrees of freedom and 18 are state variables. Moreover, there are 6 linear inequality constraints, 12 linear equality constraints, 6 non-linear (quadratic) equality constraints, and a quadratic objective. We have extended 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, in accordance with our earlier note, indicates a proper partitioning of the model variables into (first-stage and 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 straightforward 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. For all Param objects declared after this statement, the attribute mutable is set to True by default. Hence, non-mutable Param objects are now declared by explicitly passing the argument mutable=False to the Param constructor.

>>> # === 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 the 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 object m.p (itself an indexed Param object).

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

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 use 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 NLP optimizers can be automatically used in the event the primary subordinate local or global optimizer passed to the PyROS solve() method does not successfully solve a subproblem to an appropriate termination condition. These alternative solvers are provided through the optional keyword arguments backup_local_solvers and backup_global_solvers.

The final step in solving a model with PyROS is to construct 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 defined through the pyrosTerminationCondition enumeration and tabulated below.

Table 3 PyROS termination conditions.
Termination Condition Description
robust_optimal The final solution is robust optimal
robust_feasible The final solution is robust feasible
robust_infeasible The posed problem is robust infeasible
max_iter Maximum number of GRCS iteration reached
time_out Maximum number of time reached
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-stage
>>> #     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,
...     objective_focus=pyros.ObjectiveType.worst_case,
...     solve_master_globally=True,
...     load_solution=False,
... )
===========================================================================================
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(f"Final objective value: {single_stage_final_objective}")
Final objective value: 48367380.0
>>> print(f"PyROS termination condition: {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:

  • iterations: total iterations of the algorithm
  • time: total wallclock time (or elapsed time) in seconds
  • pyros_termination_condition: the GRCS algorithm termination condition
  • final_objective_value: the final objective function value.

The preceding code snippet demonstrates how to retrieve this information.

If we pass load_solution=True (the default setting) to the solve() method, then the solution at which PyROS terminates will be loaded to the variables of the original deterministic model. Note that in the preceding code snippet, we set load_solution=False to ensure the next set of runs shown here can utilize the initial point loaded to the original deterministic model, as the initial point may 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^{\text{nom}}\). 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 = \mathrm{argmax}_{k \in \mathcal{K}}~f_2(x,z^k,y^k,q^k)\).

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 degrees of freedom. PyROS handles second-stage degrees of freedom via the use of polynomial decision rules, of which the degree is controlled through the optional keyword argument decision_rule_order to the PyROS solve() method. In this example, we select affine decision rules by setting decision_rule_order=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,
...     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 single-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 "
...       f"objective: {percent_difference:.2f}")
Percent objective change relative to constant decision rules objective: -24...

For this example, we notice a ~25% decrease in the final objective value when switching from a static decision rule (no second-stage recourse) to an affine decision rule.

The Price of Robustness

In conjunction with standard Python control flow tools, PyROS facilitates a “price of robustness” analysis for a model of interest through the evaluation and comparison of the robust optimal objective function value across any appropriately constructed hierarchy of uncertainty sets. In this example, we consider a sequence of box uncertainty sets centered on the nominal uncertain parameter realization, such that each box is parameterized by a real value specifying a relative box size. To this end, we construct an iterable called relative_deviation_list whose entries are float values representing the relative sizes. We then loop through relative_deviation_list so that for each relative size, the corresponding robust optimal objective value can be evaluated by creating an appropriate BoxSet instance and invoking the PyROS solver:

>>> # 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,
...         objective_focus=pyros.ObjectiveType.worst_case,
...         solve_master_globally=True,
...         decision_rule_order=1,
...     )
...     is_robust_optimal = (
...         results.pyros_termination_condition
...         == pyros.pyrosTerminationCondition.robust_optimal
...     )
...     if not is_robust_optimal:
...         print(f"Instance for relative deviation: {relative_deviation} "
...               "not solved to robust optimality.")
...         robust_optimal_objectives.append("-----")
...     else:
...         robust_optimal_objectives.append(str(results.final_objective_value))

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

Table 4 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{-----}\)

Notice that PyROS was successfully able to determine the robust infeasibility of the problem under the largest uncertainty set.

o Relative Deviation from Nominal Realization

x Relative to Deterministic Optimal Objective

This example clearly illustrates the potential impact of the uncertainty set size on the robust optimal objective function value and demonstrates the ease of implementing a price of robustness study for a given optimization problem under uncertainty.

Note

Please provide feedback and/or report any problems by opening an issue on the Pyomo GitHub page.