diff --git a/doc/OnlineDocs/_static/theme_overrides.css b/doc/OnlineDocs/_static/theme_overrides.css
index 936b3f95b13..d0eddc5faf9 100644
--- a/doc/OnlineDocs/_static/theme_overrides.css
+++ b/doc/OnlineDocs/_static/theme_overrides.css
@@ -72,6 +72,10 @@ dl.py.method dt em span.n {
padding-right: 8px !important;
}
+.rst-content table.scrollwide {
+ overflow: scroll !important;
+ display: block;
+}
/* OLD theme overrides
diff --git a/doc/OnlineDocs/explanation/solvers/index.rst b/doc/OnlineDocs/explanation/solvers/index.rst
index a50a604f93e..d416603404d 100644
--- a/doc/OnlineDocs/explanation/solvers/index.rst
+++ b/doc/OnlineDocs/explanation/solvers/index.rst
@@ -6,7 +6,7 @@ Solvers
persistent
gdpopt
- pyros
+ pyrosdir/index
mindtpy
mcpp
multistart
diff --git a/doc/OnlineDocs/explanation/solvers/pyrosdir/getting_started.rst b/doc/OnlineDocs/explanation/solvers/pyrosdir/getting_started.rst
new file mode 100644
index 00000000000..09b7f71cb1b
--- /dev/null
+++ b/doc/OnlineDocs/explanation/solvers/pyrosdir/getting_started.rst
@@ -0,0 +1,547 @@
+==========================
+Getting Started with PyROS
+==========================
+
+.. contents:: Table of Contents
+ :depth: 3
+ :local:
+
+
+.. _pyros_installation:
+
+Installation
+============
+In advance of using PyROS to solve robust optimization problems,
+you will need (at least) one local nonlinear programming (NLP) solver
+(e.g.,
+`CONOPT `_,
+`IPOPT `_,
+`Knitro `_)
+and (at least) one global NLP solver
+(e.g.,
+`BARON `_,
+`COUENNE `_,
+`SCIP `_)
+installed and licensed on your system.
+
+PyROS can be installed as follows:
+
+1. :ref:`Install Pyomo `.
+ PyROS is included in the Pyomo software package, at pyomo/contrib/pyros.
+2. Install NumPy and SciPy with your preferred package manager;
+ both NumPy and SciPy are required dependencies of PyROS.
+ You may install NumPy and SciPy with, for example, ``conda``:
+
+ ::
+
+ conda install numpy scipy
+
+ or ``pip``:
+
+ ::
+
+ pip install numpy scipy
+3. (*Optional*) Test your installation:
+ install ``pytest`` and ``parameterized``
+ with your preferred package manager (as in the previous step):
+
+ ::
+
+ pip install pytest parameterized
+
+ You may then run the PyROS tests as follows:
+
+ ::
+
+ python -c 'import os, pytest, pyomo.contrib.pyros as p; pytest.main([os.path.dirname(p.__file__)])'
+
+ Some tests involving deterministic NLP solvers may be skipped
+ if
+ `IPOPT `_,
+ `BARON `_,
+ or
+ `SCIP `_
+ is not
+ pre-installed and licensed on your system.
+
+
+Quickstart
+==========
+We now provide a quick overview of how to use PyROS
+to solve a robust optimization problem.
+
+Consider the nonconvex deterministic QCQP
+
+.. math::
+ :nowrap:
+
+ \[\begin{array}{clll}
+ \displaystyle \min_{\substack{x \in [-100, 100], \\ z \in [-100, 100], \\ (y_1, y_2) \in [-100, 100]^2}}
+ & ~~ x^2 - y_1 z + y_2 & \\
+ \displaystyle \text{s.t.}
+ & ~~ xy_1 \geq 150 (q_1 + 1)^2 \\
+ & ~~ x + y_2^2 \leq 600 \\
+ & ~~ xz - q_2 y_1 = 2 \\
+ & ~~ y_1^2 - 2y_2 = 15
+ \end{array}\]
+
+in which
+:math:`x` is the sole first-stage variable,
+:math:`z` is the sole second-stage variable,
+:math:`y_1, y_2` are the state variables,
+and :math:`q_1, q_2` are the uncertain parameters.
+
+The uncertain parameters :math:`q_1, q_2`
+each have a nominal value of 1.
+We assume that :math:`q_1, q_2`
+can independently deviate from their
+nominal values by up to :math:`\pm 10\%`,
+so that :math:`(q_1, q_2)` is constrained in value to the
+interval uncertainty set :math:`\mathcal{Q} = [0.9, 1.1]^2`.
+
+.. note::
+ Per our analysis, our selections of first-stage variables
+ and second-stage variables in the present example
+ satisfy our
+ :ref:`assumption that the state variable values are
+ uniquely defined `.
+
+
+Step 0: Import Pyomo and the PyROS Module
+-----------------------------------------
+
+In anticipation of using the PyROS solver and building the deterministic Pyomo
+model:
+
+.. _pyros_quickstart_module_imports:
+
+.. doctest::
+
+ >>> import pyomo.environ as pyo
+ >>> import pyomo.contrib.pyros as pyros
+
+Step 1: Define the Solver Inputs
+--------------------------------
+
+Deterministic Model
+^^^^^^^^^^^^^^^^^^^
+
+The model can be implemented as follows:
+
+.. _pyros_quickstart_model_construct:
+
+.. doctest::
+
+ >>> m = pyo.ConcreteModel()
+ >>> # parameters
+ >>> m.q1 = pyo.Param(initialize=1, mutable=True)
+ >>> m.q2 = pyo.Param(initialize=1, mutable=True)
+ >>> # variables
+ >>> m.x = pyo.Var(bounds=[-100, 100])
+ >>> m.z = pyo.Var(bounds=[-100, 100])
+ >>> m.y1 = pyo.Var(bounds=[-100, 100])
+ >>> m.y2 = pyo.Var(bounds=[-100, 100])
+ >>> # objective
+ >>> m.obj = pyo.Objective(expr=m.x ** 2 - m.y1 * m.z + m.y2)
+ >>> # constraints
+ >>> m.ineq1 = pyo.Constraint(expr=m.x * m.y1 >= 150 * (m.q1 + 1) ** 2)
+ >>> m.ineq2 = pyo.Constraint(expr=m.x + m.y2 ** 2 <= 600)
+ >>> m.eq1 = pyo.Constraint(expr=m.x * m.z - m.y1 * m.q2 == 2)
+ >>> m.eq2 = pyo.Constraint(expr=m.y1 ** 2 - 2 * m.y2 == 15)
+
+
+Observe that the uncertain parameters :math:`q_1, q_2` are implemented
+as mutable :class:`~pyomo.core.base.param.Param` objects.
+See the
+:ref:`Uncertain parameters section of the
+Solver Interface documentation `
+for further guidance.
+
+
+First-Stage and Second-Stage Variables
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+We take ``m.x`` to be the sole first-stage variable and ``m.z``
+to be the sole second-stage variable:
+
+.. doctest::
+
+ >>> first_stage_variables = [m.x]
+ >>> second_stage_variables = [m.z]
+
+
+Uncertain Parameters
+^^^^^^^^^^^^^^^^^^^^
+The uncertain parameters are represented by ``m.q1`` and ``m.q2``:
+
+.. doctest::
+
+ >>> uncertain_params = [m.q1, m.q2]
+
+Uncertainty Set
+^^^^^^^^^^^^^^^
+As previously discussed, we take the uncertainty set to be
+the interval :math:`[0.9, 1.1]^2`,
+which we can implement as a
+:class:`~pyomo.contrib.pyros.uncertainty_sets.BoxSet` object:
+
+.. doctest::
+
+ >>> box_uncertainty_set = pyros.BoxSet(bounds=[(0.9, 1.1)] * 2)
+
+Further information on PyROS uncertainty sets is presented in the
+:ref:`Uncertainty Sets documentation `.
+
+Subordinate NLP Solvers
+^^^^^^^^^^^^^^^^^^^^^^^
+We will use IPOPT as the subordinate local NLP solver
+and BARON as the subordinate global NLP solver:
+
+.. doctest::
+ :skipif: not (baron_available and baron.license_is_valid() and ipopt_available)
+
+ >>> local_solver = pyo.SolverFactory("ipopt")
+ >>> 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 :meth:`~pyomo.contrib.pyros.pyros.PyROS.solve` method
+ does not successfully solve a subproblem to an appropriate termination
+ condition. These alternative solvers can be provided through the optional
+ keyword arguments ``backup_local_solvers`` and ``backup_global_solvers``
+ to the PyROS :meth:`~pyomo.contrib.pyros.pyros.PyROS.solve` method.
+
+In advance of using PyROS, we check that the model can be solved
+to optimality with the subordinate global solver:
+
+.. _pyros_quickstart_solve_deterministic:
+
+.. doctest::
+ :skipif: not (baron_available and baron.license_is_valid() and ipopt_available)
+
+ >>> pyo.assert_optimal_termination(global_solver.solve(m))
+ >>> deterministic_obj = pyo.value(m.obj)
+ >>> print(f"Optimal deterministic objective value: {deterministic_obj:.2f}")
+ Optimal deterministic objective value: 5407.94
+
+
+Step 2: Solve With PyROS
+------------------------
+PyROS can be instantiated through the Pyomo
+:class:`~pyomo.opt.base.solvers.SolverFactory`:
+
+.. doctest::
+
+ >>> pyros_solver = pyo.SolverFactory("pyros")
+
+Invoke PyROS
+^^^^^^^^^^^^^^^^^
+We now use PyROS to solve the model to robust optimality
+by invoking the :meth:`~pyomo.contrib.pyros.pyros.PyROS.solve`
+method of the PyROS solver object:
+
+.. _pyros_quickstart_single-stage-problem:
+
+.. doctest::
+ :skipif: not (baron_available and baron.license_is_valid() and ipopt_available)
+
+ >>> results_1 = pyros_solver.solve(
+ ... # required arguments
+ ... model=m,
+ ... first_stage_variables=first_stage_variables,
+ ... second_stage_variables=second_stage_variables,
+ ... uncertain_params=uncertain_params,
+ ... uncertainty_set=box_uncertainty_set,
+ ... local_solver=local_solver,
+ ... global_solver=global_solver,
+ ... # optional arguments: passed directly to
+ ... # solve to robust optimality
+ ... objective_focus="worst_case",
+ ... solve_master_globally=True,
+ ... ) # doctest: +ELLIPSIS
+ ==============================================================================
+ PyROS: The Pyomo Robust Optimization Solver...
+ ...
+ Robust optimal solution identified.
+ ...
+ All done. Exiting PyROS.
+ ==============================================================================
+
+
+PyROS, by default, logs to the output console the progress of the optimization
+and, upon termination, a summary of the final result.
+The summary includes the iteration and solve time requirements,
+the final objective function value, and the termination condition.
+For further information on the output log,
+see the :ref:`Solver Output Log documentation `.
+
+
+.. note::
+
+ PyROS, like other Pyomo solvers, accepts optional arguments
+ passed indirectly through the keyword argument ``options``.
+ This is discussed further in the
+ :ref:`Optional Arguments section of the
+ Solver Interface documentation `.
+ Thus, the PyROS solver invocation in the
+ :ref:`preceding code snippet `
+ is equivalent to:
+
+ .. code-block::
+
+ results_1 = pyros_solver.solve(
+ model=m,
+ first_stage_variables=first_stage_variables,
+ second_stage_variables=second_stage_variables,
+ uncertain_params=uncertain_params,
+ uncertainty_set=box_uncertainty_set,
+ local_solver=local_solver,
+ global_solver=global_solver,
+ # optional arguments: passed indirectly to
+ # solve to robust optimality
+ options={
+ "objective_focus": "worst_case",
+ "solve_master_globally": True,
+ },
+ )
+
+
+Inspect the Results
+^^^^^^^^^^^^^^^^^^^
+The PyROS :meth:`~pyomo.contrib.pyros.pyros.PyROS.solve` method
+returns a results object,
+of type :class:`~pyomo.contrib.pyros.solve_data.ROSolveResults`,
+that summarizes the outcome of invoking PyROS on a robust optimization problem.
+By default, a printout of the results object is included at the end of the solver
+output log.
+Alternatively, we can display the results object ourselves using:
+
+.. code::
+
+ >>> print(results_1) # output may vary # doctest: +SKIP
+ Termination stats:
+ Iterations : 3
+ Solve time (wall s) : 0.917
+ Final objective value : 9.6616e+03
+ Termination condition : pyrosTerminationCondition.robust_optimal
+
+
+We can also query the results object's individual attributes:
+
+.. code::
+
+ >>> results_1.iterations # total number of iterations
+ 3
+ >>> results_1.time # total wallclock time; may vary # doctest: +SKIP
+ 0.917
+ >>> results_1.final_objective_value # final objective value; may vary # doctest: +ELLIPSIS
+ 9661.6...
+ >>> results_1.pyros_termination_condition # termination condition
+
+
+Since PyROS has successfully solved our problem,
+the final solution has been automatically loaded to the model.
+We can inspect the resulting state of the model
+by invoking, for example, ``m.display()`` or ``m.pprint()``.
+
+For a general discussion of the PyROS solver outputs,
+see the
+:ref:`Overview of Outputs section of the
+Solver Interface documentation `.
+
+
+Try Higher-Order Decision Rules
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+PyROS uses polynomial decision rules to approximate the adjustability
+of the second-stage variables to the uncertain parameters.
+The degree of the decision rule polynomials is
+specified through the optional keyword argument
+``decision_rule_order`` to the PyROS
+:meth:`~pyomo.contrib.pyros.pyros.PyROS.solve` method.
+By default, ``decision_rule_order`` is set to 0,
+so that static decision rules are used.
+Increasing the decision rule order
+may yield a solution with better quality:
+
+.. _pyros_quickstart_example-two-stg:
+
+.. doctest::
+ :skipif: not (baron_available and baron.license_is_valid() and ipopt_available)
+
+ >>> results_2 = pyros_solver.solve(
+ ... model=m,
+ ... first_stage_variables=first_stage_variables,
+ ... second_stage_variables=second_stage_variables,
+ ... uncertain_params=uncertain_params,
+ ... uncertainty_set=box_uncertainty_set,
+ ... local_solver=local_solver,
+ ... global_solver=global_solver,
+ ... objective_focus="worst_case",
+ ... solve_master_globally=True,
+ ... decision_rule_order=1, # use affine decision rules
+ ... ) # doctest: +ELLIPSIS
+ ==============================================================================
+ PyROS: The Pyomo Robust Optimization Solver...
+ ...
+ Robust optimal solution identified.
+ ...
+ All done. Exiting PyROS.
+ ==============================================================================
+
+
+Inspecting the results:
+
+.. code::
+
+ >>> print(results_2) # output may vary # doctest: +SKIP
+ Termination stats:
+ Iterations : 5
+ Solve time (wall s) : 2.730
+ Final objective value : 6.5403e+03
+ Termination condition : pyrosTerminationCondition.robust_optimal
+
+
+Notice that when we switch from optimizing over static decision rules
+to optimizing over affine decision rules,
+there is a ~32% decrease in the final objective
+value, albeit at some additional computational expense.
+
+
+Analyzing the Price of Robustness
+---------------------------------
+In conjunction with standard Pyomo control flow tools,
+PyROS facilitates an analysis of the "price of robustness",
+which we define to be the increase in the robust optimal objective value
+relative to the deterministically optimal objective value.
+
+Let us, for example, consider optimizing robustly against
+an interval uncertainty set :math:`[1 - p, 1 + p]^2`,
+where :math:`p` is the half-length of the interval.
+We can optimize against intervals of increasing half-length :math:`p`
+by iterating over select values for :math:`p` in a ``for`` loop,
+and in each iteration, solving a robust optimization problem
+subject to a corresponding
+:class:`~pyomo.contrib.pyros.uncertainty_sets.BoxSet` instance:
+
+.. doctest::
+ :skipif: not (baron_available and baron.license_is_valid() and ipopt_available)
+
+ >>> results_dict = dict()
+ >>> for half_length in [0.0, 0.1, 0.2, 0.3, 0.4]:
+ ... print(f"Solving problem for {half_length=}:")
+ ... results_dict[half_length] = pyros_solver.solve(
+ ... model=m,
+ ... first_stage_variables=first_stage_variables,
+ ... second_stage_variables=second_stage_variables,
+ ... uncertain_params=uncertain_params,
+ ... uncertainty_set=pyros.BoxSet(
+ ... bounds=[(1 - half_length, 1 + half_length)] * 2
+ ... ),
+ ... local_solver=local_solver,
+ ... global_solver=global_solver,
+ ... objective_focus="worst_case",
+ ... solve_master_globally=True,
+ ... decision_rule_order=1,
+ ... ) # doctest: +ELLIPSIS
+ ...
+ Solving problem for half_length=0.0:
+ ...
+ Solving problem for half_length=0.1:
+ ...
+ Solving problem for half_length=0.2:
+ ...
+ Solving problem for half_length=0.3:
+ ...
+ Solving problem for half_length=0.4:
+ ...
+ All done. Exiting PyROS.
+ ==============================================================================
+
+Using the :py:obj:`dict` populated in the loop,
+and the
+:ref:`previously evaluated deterministically optimal
+objective value `,
+we can print a tabular summary of the results:
+
+.. doctest::
+ :skipif: not (baron_available and baron.license_is_valid() and ipopt_available)
+
+ >>> for idx, (half_length, res) in enumerate(results_dict.items()):
+ ... if idx == 0:
+ ... # print table header
+ ... print("=" * 71)
+ ... print(
+ ... f"{'Half-Length':15s}"
+ ... f"{'Termination Cond.':21s}"
+ ... f"{'Objective Value':18s}"
+ ... f"{'Price of Rob. (%)':17s}"
+ ... )
+ ... print("-" * 71)
+ ... # print table row
+ ... obj_value, percent_obj_increase = float("nan"), float("nan")
+ ... is_robust_optimal = (
+ ... res.pyros_termination_condition
+ ... == pyros.pyrosTerminationCondition.robust_optimal
+ ... )
+ ... is_robust_infeasible = (
+ ... res.pyros_termination_condition
+ ... == pyros.pyrosTerminationCondition.robust_infeasible
+ ... )
+ ... if is_robust_optimal:
+ ... # compute the price of robustness
+ ... obj_value = res.final_objective_value
+ ... price_of_robustness = (
+ ... (res.final_objective_value - deterministic_obj)
+ ... / deterministic_obj
+ ... )
+ ... elif is_robust_infeasible:
+ ... # infinite objective
+ ... obj_value, price_of_robustness = float("inf"), float("inf")
+ ... print(
+ ... f"{half_length:<15.1f}"
+ ... f"{res.pyros_termination_condition.name:21s}"
+ ... f"{obj_value:<18.2f}"
+ ... f"{100 * price_of_robustness:<.2f}"
+ ... )
+ ... print("-" * 71)
+ ...
+ =======================================================================
+ Half-Length Termination Cond. Objective Value Price of Rob. (%)
+ -----------------------------------------------------------------------
+ 0.0 robust_optimal 5407.94 -0.00
+ -----------------------------------------------------------------------
+ 0.1 robust_optimal 6540.31 20.94
+ -----------------------------------------------------------------------
+ 0.2 robust_optimal 7838.50 44.94
+ -----------------------------------------------------------------------
+ 0.3 robust_optimal 9316.88 72.28
+ -----------------------------------------------------------------------
+ 0.4 robust_infeasible inf inf
+ -----------------------------------------------------------------------
+
+
+The table shows the response of the PyROS termination condition,
+final objective value, and price of robustness
+to the half-length :math:`p`.
+Observe that:
+
+* The optimal objective value for the interval of half-length
+ :math:`p=0` is equal to the optimal deterministic objective value
+* The objective value (and thus, the price of robustness)
+ increases with the half-length
+* For large enough half-length (:math:`p=0.4`), the problem
+ is robust infeasible
+
+Therefore, this example clearly illustrates the potential
+impact of the uncertainty set size on the robust optimal
+objective function value
+and the ease of analyzing the price of robustness
+for a given optimization problem under uncertainty.
+
+
+Beyond the Basics
+=================
+A more in-depth guide to incorporating PyROS into a
+Pyomo optimization workflow is given
+in the :ref:`Usage Tutorial `.
diff --git a/doc/OnlineDocs/explanation/solvers/pyrosdir/index.rst b/doc/OnlineDocs/explanation/solvers/pyrosdir/index.rst
new file mode 100644
index 00000000000..494bd83eb4d
--- /dev/null
+++ b/doc/OnlineDocs/explanation/solvers/pyrosdir/index.rst
@@ -0,0 +1,38 @@
+============
+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) `_
+and `Carbon Capture Simulation for Industry Impact (CCSI2)
+`_
+projects.
+
+.. toctree::
+ :maxdepth: 1
+ :caption: Index of PyROS Documentation
+
+ Methodology Overview
+ Getting Started
+ Solver Interface
+ Uncertainty Sets
+ Solver Output Log
+ Usage Tutorial
+
+
+Citing PyROS
+============
+If you use PyROS in your research,
+please acknowledge PyROS by citing [IAE+21]_.
+
+
+Feedback and Reporting Issues
+=============================
+Please provide feedback and/or report any problems by
+`opening an issue on the Pyomo GitHub page `_.
diff --git a/doc/OnlineDocs/explanation/solvers/pyrosdir/overview.rst b/doc/OnlineDocs/explanation/solvers/pyrosdir/overview.rst
new file mode 100644
index 00000000000..75a24ca8a43
--- /dev/null
+++ b/doc/OnlineDocs/explanation/solvers/pyrosdir/overview.rst
@@ -0,0 +1,110 @@
+.. _pyros_overview:
+
+==========================
+PyROS Methodology Overview
+==========================
+
+PyROS can accommodate optimization models with:
+
+* **Continuous variables** only
+* **Nonlinearities** (including **nonconvexities**) in both the
+ variables and uncertain parameters
+* **First-stage degrees of freedom** and **second-stage degrees of freedom**
+* **Equality constraints** defining state variables,
+ including implicitly defined state variables that cannot be
+ eliminated from the model via reformulation
+* **Uncertain parameters** participating in the inequality constraints,
+ equality constraints, and/or objective function
+
+Supported deterministic models are nonlinear programs (NLPs)
+of the general form
+
+.. _deterministic-model:
+
+.. math::
+ :nowrap:
+
+ \[\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}\]
+
+where:
+
+* :math:`x \in \mathcal{X}` denotes the first-stage degree of freedom variables
+ (or design variables),
+ of which the feasible space :math:`\mathcal{X} \subseteq \mathbb{R}^{n_x}`
+ is defined by the model constraints
+ (including variable bounds specifications) referencing :math:`x` only
+* :math:`z \in \mathbb{R}^{n_z}` denotes the second-stage degree of freedom
+ variables (or control variables)
+* :math:`y \in \mathbb{R}^{n_y}` denotes the state variables
+* :math:`q \in \mathbb{R}^{n_q}` denotes the model parameters considered
+ uncertain, and :math:`q^{\text{nom}}` is the vector of nominal values
+ associated with those
+* :math:`f_1\left(x\right)` is the summand of the objective function that depends
+ only on the first-stage degree of freedom variables
+* :math:`f_2\left(x, z, y; q\right)` is the summand of the objective function
+ that depends on all variables and the uncertain parameters
+* :math:`g_i\left(x, z, y; q\right)` is the :math:`i^\text{th}`
+ inequality constraint function in set :math:`\mathcal{I}`
+ (see :ref:`first Note `)
+* :math:`h_j\left(x, z, y; q\right)` is the :math:`j^\text{th}`
+ equality constraint function in set :math:`\mathcal{J}`
+ (see :ref:`second Note `)
+
+.. _var-bounds-to-ineqs:
+
+.. note::
+
+ PyROS accepts and automatically reformulates models with:
+
+ 1. Interval bounds on components of :math:`(x, z, y)`
+ 2. Ranged inequality constraints
+
+
+.. _pyros_unique_state_vars:
+
+.. note::
+ A key assumption of PyROS is that
+ for every
+ :math:`x \in \mathcal{X}`,
+ :math:`z \in \mathbb{R}^{n_z}`,
+ :math:`q \in \mathcal{Q}`,
+ there exists a unique :math:`y \in \mathbb{R}^{n_y}`
+ for which :math:`(x, z, y, q)`
+ satisfies the equality constraints
+ :math:`h_j(x, z, y, q) = 0\,\,\forall\, j \in \mathcal{J}`.
+ If this assumption is not met,
+ then the selection of state
+ (i.e., not degree of freedom) variables :math:`y` is incorrect,
+ and one or more entries of :math:`y` should be appropriately
+ redesignated to be part of either :math:`x` or :math:`z`.
+
+In order to cast the robust optimization counterpart of the
+:ref:`deterministic model `,
+we now assume that the uncertain parameters :math:`q` may attain
+any realization in a compact uncertainty set
+:math:`\mathcal{Q} \subseteq \mathbb{R}^{n_q}` containing
+the nominal value :math:`q^{\text{nom}}`.
+The set :math:`\mathcal{Q}` may be **either continuous or discrete**.
+
+Based on the above notation,
+the form of the robust counterpart addressed by PyROS is
+
+.. math::
+ :nowrap:
+
+ \[\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}\]
+
+PyROS accepts a deterministic model and accompanying uncertainty set
+and then, using the Generalized Robust Cutting-Set algorithm developed
+in [IAE+21]_, seeks a solution to the robust counterpart.
diff --git a/doc/OnlineDocs/explanation/solvers/pyrosdir/solver_interface.rst b/doc/OnlineDocs/explanation/solvers/pyrosdir/solver_interface.rst
new file mode 100644
index 00000000000..85b28829386
--- /dev/null
+++ b/doc/OnlineDocs/explanation/solvers/pyrosdir/solver_interface.rst
@@ -0,0 +1,273 @@
+.. _pyros_solver_interface:
+
+======================
+PyROS Solver Interface
+======================
+
+.. contents:: Table of Contents
+ :depth: 2
+ :local:
+
+Instantiation
+=============
+
+The PyROS solver is invoked through the
+:meth:`~pyomo.contrib.pyros.pyros.PyROS.solve` method
+of an instance of the PyROS solver class, which can be
+instantiated as follows:
+
+.. code-block::
+
+ import pyomo.environ as pyo
+ import pyomo.contrib.pyros as pyros # register the PyROS solver
+ pyros_solver = pyo.SolverFactory("pyros")
+
+
+Overview of Inputs
+==================
+
+Deterministic Model
+-------------------
+PyROS is designed to operate on a single-objective deterministic model
+(implemented as a :class:`~pyomo.core.base.PyomoModel.ConcreteModel`),
+from which the robust optimization counterpart is automatically inferred.
+All variables of the model should be continuous, as
+mixed-integer problems are not supported.
+
+First-Stage and Second-Stage Variables
+--------------------------------------
+A model may have either first-stage variables,
+second-stage variables, or both.
+PyROS automatically considers all other variables participating
+in the active model components to be state variables.
+Further, PyROS assumes that the state variables are
+:ref:`uniquely defined by the equality constraints `.
+
+
+.. _pyros_uncertain_params:
+
+Uncertain Parameters
+--------------------
+Uncertain parameters can be represented by either
+mutable :class:`~pyomo.core.base.param.Param`
+or fixed :class:`~pyomo.core.base.var.Var` objects.
+Uncertain parameters *cannot* be directly
+represented by Python literals that have been hard-coded into the
+deterministic model.
+
+A :class:`~pyomo.core.base.param.Param` object can be made mutable
+at construction by passing the argument ``mutable=True`` to the
+:class:`~pyomo.core.base.param.Param` constructor.
+If specifying/modifying the ``mutable`` argument
+is not straightforward in your context,
+then add the following lines of code to your script
+before setting up your deterministic model:
+
+
+.. code-block::
+
+ import pyomo.environ as pyo
+ pyo.Param.DefaultMutable = True
+
+All :class:`~pyomo.core.base.param.Param` objects declared
+after the preceding code statements will be made mutable by default.
+
+
+Uncertainty Set
+---------------
+The uncertainty set is represented by an
+:class:`~pyomo.contrib.pyros.uncertainty_sets.UncertaintySet`
+object.
+See the :ref:`Uncertainty Sets documentation `
+for more information.
+
+Subordinate NLP Solvers
+-----------------------
+PyROS requires at least one subordinate
+local nonlinear programming (NLP) solver (e.g., Ipopt or CONOPT)
+and subordinate global NLP solver (e.g., BARON or SCIP)
+to solve subproblems.
+
+.. note::
+
+ In advance of invoking the PyROS solver,
+ check that your deterministic model can be solved
+ to optimality by either your subordinate local or global
+ NLP solver.
+
+
+.. _pyros_optional_arguments:
+
+Optional Arguments
+------------------
+The optional arguments are enumerated in the documentation of the
+:meth:`~pyomo.contrib.pyros.pyros.PyROS.solve` method.
+
+Like other Pyomo solver interface methods,
+the PyROS
+:meth:`~pyomo.contrib.pyros.pyros.PyROS.solve`
+method
+accepts the keyword argument ``options``,
+which must be a :class:`dict`
+mapping names of optional arguments to
+:meth:`~pyomo.contrib.pyros.pyros.PyROS.solve`
+to their desired values.
+If an argument is passed directly by keyword and
+indirectly through ``options``,
+then the value passed directly takes precedence over the
+value passed through ``options``.
+
+.. warning::
+
+ All required arguments to the PyROS
+ :meth:`~pyomo.contrib.pyros.pyros.PyROS.solve` method
+ must be passed directly by position or keyword,
+ or else an exception is raised.
+ Required arguments passed indirectly through the ``options``
+ setting are ignored.
+
+
+Separation Priority Ordering
+----------------------------
+The PyROS solver supports custom prioritization of
+the separation subproblems (and, thus, the constraints)
+that are automatically derived from
+a given model for robust optimization.
+Users may specify separation priorities through:
+
+- (Recommended) :class:`~pyomo.core.base.suffix.Suffix` components
+ with local name ``pyros_separation_priority``,
+ declared on the model or any of its sub-blocks.
+ Each entry of every such
+ :class:`~pyomo.core.base.suffix.Suffix`
+ should map a
+ :class:`~pyomo.core.base.var.Var`
+ or :class:`~pyomo.core.base.constraint.Constraint`
+ component to a value that specifies the separation
+ priority of all constraints derived from that component
+- The optional argument ``separation_priority_order``
+ to the PyROS :py:meth:`~pyomo.contrib.pyros.pyros.PyROS.solve`
+ method. The argument should be castable to a :py:obj:`dict`,
+ of which each entry maps the full name of a
+ :class:`~pyomo.core.base.var.Var`
+ or :class:`~pyomo.core.base.constraint.Constraint`
+ component to a value that specifies the
+ separation priority of all constraints
+ derived from that component
+
+Specification via :class:`~pyomo.core.base.suffix.Suffix` components
+takes precedence over specification via the solver argument
+``separation_priority_order``.
+Moreover, the precedence ordering among
+:class:`~pyomo.core.base.suffix.Suffix`
+components is handled by the Pyomo
+:class:`~pyomo.core.base.suffix.SuffixFinder` utility.
+
+A separation priority can be either
+a (real) number (i.e., of type :py:class:`int`, :py:class:`float`, etc.)
+or :py:obj:`None`.
+A higher number indicates a higher priority.
+The default priority for all constraints is 0.
+Therefore a constraint can be prioritized [or deprioritized]
+over the default by mapping the constraint to a positive [or negative] number.
+In practice, critical or dominant constraints are often
+prioritized over algorithmic or implied constraints.
+
+Constraints that have been assigned a priority of :py:obj:`None`
+are enforced subject to only the nominal uncertain parameter realization
+provided by the user. Therefore, these constraints are not imposed robustly
+and, in particular, are excluded from the separation problems.
+
+
+.. _pyros_solver_outputs:
+
+Overview of Outputs
+===================
+
+.. _pyros_output_results_object:
+
+Results Object
+--------------
+The :meth:`~pyomo.contrib.pyros.pyros.PyROS.solve` method returns
+an :class:`~pyomo.contrib.pyros.solve_data.ROSolveResults` object.
+
+When the PyROS :meth:`~pyomo.contrib.pyros.pyros.PyROS.solve` method
+has successfully solved a given robust optimization problem,
+the
+:attr:`~pyomo.contrib.pyros.solve_data.ROSolveResults.pyros_termination_condition`
+attribute of the returned
+:attr:`~pyomo.contrib.pyros.solve_data.ROSolveResults`
+object is set to
+:attr:`~pyomo.contrib.pyros.util.pyrosTerminationCondition.robust_optimal`
+if and only if:
+
+1. Master problems are solved to global optimality
+ (by passing ``solve_master_globally=True``)
+2. A worst-case objective focus is chosen
+ (by setting ``objective_focus``
+ to :attr:`~pyomo.contrib.pyros.util.ObjectiveType.worst_case`)
+
+Otherwise, the termination condition is set to
+:attr:`~pyomo.contrib.pyros.util.pyrosTerminationCondition.robust_feasible`.
+
+The
+:attr:`~pyomo.contrib.pyros.solve_data.ROSolveResults.final_objective_value`
+attribute of the results object depends on
+the value of the optional ``objective_focus`` argument to the
+:meth:`~pyomo.contrib.pyros.pyros.PyROS.solve` method:
+
+* If ``objective_focus`` is set to
+ :attr:`~pyomo.contrib.pyros.util.ObjectiveType.nominal`,
+ then the objective is evaluated subject to
+ the nominal uncertain parameter realization
+* If ``objective_focus`` is set to
+ :attr:`~pyomo.contrib.pyros.util.ObjectiveType.worst_case`,
+ then the objective is evaluated subject
+ the uncertain parameter realization that induces the worst-case
+ objective value
+
+The second-stage variable and state variable values in the
+:ref:`solution loaded to the model `
+are evaluated similarly.
+
+.. _pyros_output_final_solution:
+
+Final Solution
+--------------
+PyROS automatically loads the final solution found to the model
+(i.e., updates the values of the variables of the deterministic model)
+if and only if:
+
+1. The argument ``load_solution=True`` has been passed to PyROS
+ (occurs by default)
+2. The
+ :attr:`~pyomo.contrib.pyros.solve_data.ROSolveResults.pyros_termination_condition`
+ attribute of the returned
+ :attr:`~pyomo.contrib.pyros.solve_data.ROSolveResults` object
+ is either
+ :attr:`~pyomo.contrib.pyros.util.pyrosTerminationCondition.robust_optimal`
+ or
+ :attr:`~pyomo.contrib.pyros.util.pyrosTerminationCondition.robust_feasible`
+
+Otherwise, the solution is lost.
+
+If a solution is loaded to the model,
+then,
+as mentioned in our discussion of the
+:ref:`results object `,
+the second-stage variables and state variables
+of the model are updated according to
+the value of the optional ``objective_focus`` argument to
+the :meth:`~pyomo.contrib.pyros.pyros.PyROS.solve` method.
+The uncertain parameter objects are left unchanged.
+
+
+Solver Output Log
+-----------------
+When the PyROS
+:meth:`~pyomo.contrib.pyros.pyros.PyROS.solve` method
+is invoked to solve an RO problem,
+the progress and final result are reported through a highly
+configurable logging system.
+See the :ref:`Solver Output Log documentation `
+documentation for more information.
diff --git a/doc/OnlineDocs/explanation/solvers/pyrosdir/solver_log.rst b/doc/OnlineDocs/explanation/solvers/pyrosdir/solver_log.rst
new file mode 100644
index 00000000000..95d3d39b4fa
--- /dev/null
+++ b/doc/OnlineDocs/explanation/solvers/pyrosdir/solver_log.rst
@@ -0,0 +1,282 @@
+.. _pyros_solver_log:
+
+=======================
+PyROS Solver Output Log
+=======================
+
+.. contents:: Table of Contents
+ :depth: 1
+ :local:
+
+
+.. _pyros_solver_log_appearance:
+
+Default Format
+==============
+
+When the PyROS
+:meth:`~pyomo.contrib.pyros.pyros.PyROS.solve` method
+is called to solve a robust optimization problem,
+your console output will, by default, look like this:
+
+
+.. _solver-log-snippet:
+
+.. code-block:: text
+ :caption: PyROS solver output log for the :ref:`Quickstart example `.
+ :linenos:
+
+ ==============================================================================
+ PyROS: The Pyomo Robust Optimization Solver, v1.3.9.
+ Pyomo version: 6.9.5.dev0 (devel {main})
+ Commit hash: unknown
+ Invoked at UTC 2025-09-21T00:00:00.000000+00:00
+
+ Developed by: Natalie M. Isenberg (1), Jason A. F. Sherman (1),
+ John D. Siirola (2), Chrysanthos E. Gounaris (1)
+ (1) Carnegie Mellon University, Department of Chemical Engineering
+ (2) Sandia National Laboratories, Center for Computing Research
+
+ The developers gratefully acknowledge support from the U.S. Department
+ of Energy's Institute for the Design of Advanced Energy Systems (IDAES).
+ ==============================================================================
+ ================================= DISCLAIMER =================================
+ PyROS is still under development.
+ Please provide feedback and/or report any issues by creating a ticket at
+ https://github.com/Pyomo/pyomo/issues/new/choose
+ ==============================================================================
+ Solver options:
+ time_limit=None
+ keepfiles=False
+ tee=False
+ load_solution=True
+ symbolic_solver_labels=False
+ objective_focus=
+ nominal_uncertain_param_vals=[1, 1]
+ decision_rule_order=1
+ solve_master_globally=True
+ max_iter=-1
+ robust_feasibility_tolerance=0.0001
+ separation_priority_order={}
+ progress_logger=
+ backup_local_solvers=[]
+ backup_global_solvers=[]
+ subproblem_file_directory=None
+ subproblem_format_options={'bar': {'symbolic_solver_labels': True}}
+ bypass_local_separation=False
+ bypass_global_separation=False
+ p_robustness={}
+ ------------------------------------------------------------------------------
+ Preprocessing...
+ Done preprocessing; required wall time of 0.587s.
+ ------------------------------------------------------------------------------
+ Model Statistics:
+ Number of variables : 8
+ Epigraph variable : 1
+ First-stage variables : 1
+ Second-stage variables : 1 (1 adj.)
+ State variables : 2 (2 adj.)
+ Decision rule variables : 3
+ Number of uncertain parameters : 2 (2 eff.)
+ Number of constraints : 12
+ Equality constraints : 3
+ Coefficient matching constraints : 0
+ Other first-stage equations : 0
+ Second-stage equations : 2
+ Decision rule equations : 1
+ Inequality constraints : 9
+ First-stage inequalities : 0
+ Second-stage inequalities : 9
+ ------------------------------------------------------------------------------
+ Itn Objective 1-Stg Shift 2-Stg Shift #CViol Max Viol Wall Time (s)
+ ------------------------------------------------------------------------------
+ 0 5.4079e+03 - - 3 7.9226e+00 0.815
+ 1 5.4079e+03 6.0451e-10 1.0717e-10 2 1.0250e-01 1.194
+ 2 6.5403e+03 1.0018e-01 7.4564e-03 1 1.7142e-03 1.604
+ 3 6.5403e+03 1.9372e-16 3.6853e-06 2 1.6673e-03 1.993
+ 4 6.5403e+03 0.0000e+00 2.9067e-06 0 9.8487e-05g 2.969
+ ------------------------------------------------------------------------------
+ Robust optimal solution identified.
+ ------------------------------------------------------------------------------
+ Timing breakdown:
+
+ Identifier ncalls cumtime percall %
+ -----------------------------------------------------------
+ main 1 2.970 2.970 100.0
+ ------------------------------------------------------
+ dr_polishing 4 0.227 0.057 7.6
+ global_separation 9 0.486 0.054 16.4
+ local_separation 45 0.739 0.016 24.9
+ master 5 0.672 0.134 22.6
+ master_feasibility 4 0.095 0.024 3.2
+ preprocessing 1 0.587 0.587 19.8
+ other n/a 0.164 n/a 5.5
+ ======================================================
+ ===========================================================
+
+ ------------------------------------------------------------------------------
+ Termination stats:
+ Iterations : 5
+ Solve time (wall s) : 2.970
+ Final objective value : 6.5403e+03
+ Termination condition : pyrosTerminationCondition.robust_optimal
+ ------------------------------------------------------------------------------
+ All done. Exiting PyROS.
+ ==============================================================================
+
+
+Observe that the log contains the following information
+(listed in order of appearance):
+
+
+* **Introductory information and disclaimer** (lines 1--19):
+ Includes the version number, author
+ information, (UTC) time at which the solver was invoked,
+ and, if available, information on the local Git branch and
+ commit hash.
+* **Summary of solver options** (lines 20--41): Enumeration of
+ specifications for optional arguments to the solver.
+* **Preprocessing information** (lines 42--44):
+ Wall time required for preprocessing
+ the deterministic model and associated components,
+ i.e., standardizing model components and adding the decision rule
+ variables and equations.
+* **Model component statistics** (lines 45--62):
+ Breakdown of model component statistics.
+ Includes components added by PyROS, such as the decision rule variables
+ and equations.
+ The preprocessor may find that some second-stage variables
+ and state variables are mathematically
+ not adjustable to the uncertain parameters.
+ To this end, in the logs, the numbers of
+ adjustable second-stage variables and state variables
+ are included in parentheses, next to the total numbers
+ of second-stage variables and state variables, respectively;
+ note that "adjustable" has been abbreviated as "adj."
+* **Iteration log table** (lines 63--70):
+ Summary information on the problem iterates and subproblem outcomes.
+ The constituent columns are defined in detail in
+ :ref:`the table that follows `.
+* **Termination message** (lines 71--72): One-line message briefly summarizing
+ the reason the solver has terminated.
+* **Timing statistics** (lines 73--89):
+ Tabulated breakdown of the solver timing statistics, based on a
+ :class:`pyomo.common.timing.HierarchicalTimer` printout.
+ The identifiers are as follows:
+
+ * ``main``: Time elapsed by the solver.
+ * ``main.dr_polishing``: Time spent by the subordinate solvers
+ on polishing of the decision rules.
+ * ``main.global_separation``: Time spent by the subordinate solvers
+ on global separation subproblems.
+ * ``main.local_separation``: Time spent by the subordinate solvers
+ on local separation subproblems.
+ * ``main.master``: Time spent by the subordinate solvers on
+ the master problems.
+ * ``main.master_feasibility``: Time spent by the subordinate solvers
+ on the master feasibility problems.
+ * ``main.preprocessing``: Preprocessing time.
+ * ``main.other``: Overhead time.
+
+* **Final result** (lines 90--95):
+ A printout of the
+ :class:`~pyomo.contrib.pyros.solve_data.ROSolveResults`
+ object that is finally returned.
+* **Exit message** (lines 96--97): Confirmation that the
+ solver has been exited properly.
+
+The iteration log table is designed to provide, in a concise manner,
+important information about the progress of the iterative algorithm for
+the problem of interest.
+The constituent columns are defined in the
+table below.
+
+.. _pyros_table-iteration-log-columns:
+
+.. list-table:: PyROS iteration log table columns.
+ :widths: 10 50
+ :header-rows: 1
+
+ * - Column Name
+ - Definition
+ * - Itn
+ - Iteration number, equal to one less than the total number of elapsed
+ iterations.
+ * - Objective
+ - Master solution objective function value.
+ If the objective of the deterministic model provided
+ has a maximization sense,
+ then the negative of the objective function value is displayed.
+ Expect this value to trend upward as the iteration number
+ increases.
+ A dash ("-") is produced in lieu of a value if the master
+ problem of the current iteration is not solved successfully.
+ * - 1-Stg Shift
+ - Infinity norm of the relative difference between the first-stage
+ variable vectors of the master solutions of the current
+ and previous iterations. Expect this value to trend
+ downward as the iteration number increases.
+ A dash ("-") is produced in lieu of a value
+ if the current iteration number is 0,
+ there are no first-stage variables,
+ or the master problem of the current iteration is not solved successfully.
+ * - 2-Stg Shift
+ - Infinity norm of the relative difference between the second-stage
+ variable vectors (evaluated subject to the nominal uncertain
+ parameter realization) of the master solutions of the current
+ and previous iterations. Expect this value to trend
+ downward as the iteration number increases.
+ A dash ("-") is produced in lieu of a value
+ if the current iteration number is 0,
+ there are no second-stage variables,
+ or the master problem of the current iteration is not solved successfully.
+ An asterisk ("*") is appended to the value if decision rule
+ polishing was unsuccessful.
+ * - #CViol
+ - Number of second-stage inequality constraints found to be violated during
+ the separation step of the current iteration.
+ Unless a custom prioritization of the model's second-stage inequality
+ constraints is specified (through the ``separation_priority_order`` argument),
+ expect this number to trend downward as the iteration number increases.
+ A "+" is appended if not all of the separation problems
+ were solved successfully, either due to custom prioritization, a time out,
+ or an issue encountered by the subordinate optimizers.
+ A dash ("-") is produced in lieu of a value if the separation
+ routine is not invoked during the current iteration.
+ * - Max Viol
+ - Maximum scaled second-stage inequality constraint violation.
+ Expect this value to trend downward as the iteration number increases.
+ A 'g' is appended to the value if the separation problems were solved
+ globally during the current iteration.
+ A dash ("-") is produced in lieu of a value if the separation
+ routine is not invoked during the current iteration, or if there are
+ no second-stage inequality constraints.
+ * - Wall time (s)
+ - Total time elapsed by the solver, in seconds, up to the end of the
+ current iteration.
+
+.. _pyros_solver_log_verbosity:
+
+Configuring the Output Log
+==========================
+
+The PyROS solver output log is produced by the
+Python logger (:py:class:`logging.Logger`) object
+derived from the optional argument ``progress_logger``
+to the PyROS :meth:`~pyomo.contrib.pyros.pyros.PyROS.solve` method.
+By default, the PyROS solver argument ``progress_logger``
+is taken to be the :py:obj:`logging.INFO`-level
+logger with name ``"pyomo.contrib.pyros"``.
+The verbosity level of the output log can be adjusted by setting the
+:py:mod:`logging` level of the progress logger.
+For example, the level of the default logger can be set to
+:py:obj:`logging.DEBUG` with:
+
+.. code-block::
+
+ import logging
+ logging.getLogger("pyomo.contrib.pyros").setLevel(logging.DEBUG)
+
+We refer the reader to the
+:doc:`official Python logging library documentation `
+for further guidance on (customization of) Python logger objects.
diff --git a/doc/OnlineDocs/explanation/solvers/pyrosdir/tutorial/deterministic_heatmap.png b/doc/OnlineDocs/explanation/solvers/pyrosdir/tutorial/deterministic_heatmap.png
new file mode 100644
index 00000000000..b5d844bcc7a
Binary files /dev/null and b/doc/OnlineDocs/explanation/solvers/pyrosdir/tutorial/deterministic_heatmap.png differ
diff --git a/doc/OnlineDocs/explanation/solvers/pyrosdir/tutorial/dr0_heatmap.png b/doc/OnlineDocs/explanation/solvers/pyrosdir/tutorial/dr0_heatmap.png
new file mode 100644
index 00000000000..c589cd65c83
Binary files /dev/null and b/doc/OnlineDocs/explanation/solvers/pyrosdir/tutorial/dr0_heatmap.png differ
diff --git a/doc/OnlineDocs/explanation/solvers/pyrosdir/tutorial/dr1_heatmap.png b/doc/OnlineDocs/explanation/solvers/pyrosdir/tutorial/dr1_heatmap.png
new file mode 100644
index 00000000000..e7fd1d99b5b
Binary files /dev/null and b/doc/OnlineDocs/explanation/solvers/pyrosdir/tutorial/dr1_heatmap.png differ
diff --git a/doc/OnlineDocs/explanation/solvers/pyrosdir/tutorial/por_heatmaps.png b/doc/OnlineDocs/explanation/solvers/pyrosdir/tutorial/por_heatmaps.png
new file mode 100644
index 00000000000..927439b1236
Binary files /dev/null and b/doc/OnlineDocs/explanation/solvers/pyrosdir/tutorial/por_heatmaps.png differ
diff --git a/doc/OnlineDocs/explanation/solvers/pyrosdir/tutorial/por_sensitivity.png b/doc/OnlineDocs/explanation/solvers/pyrosdir/tutorial/por_sensitivity.png
new file mode 100644
index 00000000000..805bd035bcc
Binary files /dev/null and b/doc/OnlineDocs/explanation/solvers/pyrosdir/tutorial/por_sensitivity.png differ
diff --git a/doc/OnlineDocs/explanation/solvers/pyrosdir/tutorial/reactor_cooler.png b/doc/OnlineDocs/explanation/solvers/pyrosdir/tutorial/reactor_cooler.png
new file mode 100644
index 00000000000..f4c86412763
Binary files /dev/null and b/doc/OnlineDocs/explanation/solvers/pyrosdir/tutorial/reactor_cooler.png differ
diff --git a/doc/OnlineDocs/explanation/solvers/pyrosdir/tutorial/tutorial.rst b/doc/OnlineDocs/explanation/solvers/pyrosdir/tutorial/tutorial.rst
new file mode 100644
index 00000000000..8bce5f79a94
--- /dev/null
+++ b/doc/OnlineDocs/explanation/solvers/pyrosdir/tutorial/tutorial.rst
@@ -0,0 +1,942 @@
+.. _pyros_tutorial:
+
+====================
+PyROS Usage Tutorial
+====================
+
+.. contents:: Table of Contents
+ :depth: 3
+ :local:
+
+This tutorial is an in-depth guide on how to
+use PyROS to solve a two-stage robust optimization problem.
+The problem is taken from the area of chemical process systems design.
+
+Setup
+-----
+To successfully run this tutorial, you will need to
+:ref:`install PyROS `
+along with at least one local
+nonlinear programming
+(NLP) solver and at least one global NLP solver.
+In particular, this tutorial uses
+`IPOPT `__
+as the local solver,
+`BARON `__
+as the global solver,
+and `COUENNE `__
+as a backup global solver.
+
+You will also need to
+`install Matplotlib `__,
+which is used to generate plots in this tutorial.
+Further, we recommend installing
+an `interactive Matplotlib backend
+`_
+for quick and easy viewing of plots.
+
+
+Prepare the Deterministic Model
+-------------------------------
+
+PyROS is designed to operate on a user-supplied deterministic NLP. We
+now set out to prepare a deterministic NLP that can be solved tractably
+with subordinate NLP optimizers.
+
+Formulate the Model
+~~~~~~~~~~~~~~~~~~~
+
+Consider the reactor-cooler system below.
+
+.. figure:: reactor_cooler.png
+ :alt: Reactor-cooler system.
+
+ Reactor-cooler system process flowsheet,
+ adapted from `Djelassi (2020) `__.
+ Constants are set in boldface.
+
+A stream of chemical species :math:`E` enters the
+reactor with a molar flow rate :math:`F_0 = 45.36\,\text{kmol/h}`,
+absolute temperature :math:`T_0 = 333\,\text{K}`,
+concentration :math:`\boldsymbol{c}_{E0} = 32.04\, \text{kmol}/\text{m}^3`,
+and heat capacity
+:math:`\boldsymbol{c}_p = 167.4\,\text{kJ}/ \text{kmol}\,\text{K}`.
+Inside the reactor, the exothermic reaction :math:`E \to F` occurs
+at temperature :math:`T_1` and with a conversion of 90%,
+so that :math:`\boldsymbol{c}_{E1} = 0.1\boldsymbol{c}_{E0}`.
+We assume that the reaction
+follows first-order kinetics, with a rate constant
+:math:`\boldsymbol{k}_\text{R}` of nominal value
+:math:`10\,\text{hr}^{-1}` and normalized activation energy
+:math:`\boldsymbol{E/R} = 555.6\,\text{K}`.
+Further, the molar heat of reaction is
+:math:`\boldsymbol{\Delta H}_\text{R}=-23260\,\text{kJ/kmol}`.
+
+A portion of the product is cooled to a temperature :math:`T_2`
+then recycled to the reactor.
+Cooling water, with heat capacity
+:math:`\boldsymbol{c}_{w,p} = 4.184\,\text{kJ}/\text{kg}\,\text{K}`
+and inlet
+temperature :math:`\boldsymbol{T}_{w1} = 300\,\text{K}`, is used as the
+cooling medium. The heat transfer coefficient of the cooling unit is of
+nominal value
+:math:`\boldsymbol{U} = 1635\,\text{kJ}/\text{m}^2\,\text{kg}\,\text{h}`.
+
+We are interested in optimizing the design and operation of the system:
+
+- The design is specified through the reactor volume :math:`\hat{V}` and
+ the area :math:`A` of the heat exchanger used to cool the recycle stream.
+- The operation of the system can be adjusted via the reactor outlet
+ temperature :math:`T_1` and recycle stream flow rate :math:`F_1`.
+
+Modified from the formulations presented in `Halemane and Grossmann
+(1983) `__, `Yuan et
+al. (2017) `__, `Djelassi
+(2020) `__, and `Isenberg et
+al. (2021) `__, the deterministic
+optimization model for the system of interest can be written as
+
+.. math::
+
+
+ \begin{array}{cl}
+ \displaystyle\min_{\substack{
+ (\hat{V},A, c_{E1}) \in [0, 10] \times [0, 20]\\
+ (F_1, T_1) \in \mathbb{R}^{3}\\
+ (T_2, T_{w2}, \dot{Q}_{\text{HE}}, F_w, V) \in \mathbb{R}^{5}
+ }}
+ &
+ 691.2 \hat{V}^{0.7}
+ + 873.6 A^{0.6}
+ + 7.056 F_{1}
+ + 1.760 F_{w}
+ \\
+ \text{s.t.}
+ &
+ \displaystyle\boldsymbol{F}_{0} \left(\frac{\boldsymbol{c}_{E0} - \boldsymbol{c}_{E1}}{\boldsymbol{c}_{E0}}\right)
+ = V \boldsymbol{k}_{\text{R}} \exp{\left(-\frac{\boldsymbol{E}}{\boldsymbol{R}T_1}\right)}\boldsymbol{c}_{E1}
+ \\
+ & \displaystyle
+ -\boldsymbol{\Delta}\boldsymbol{H}_{\text{R}}\boldsymbol{F}_{0}
+ \left(\frac{\boldsymbol{c}_{E0} - \boldsymbol{c}_{E1}}{\boldsymbol{c}_{E0}}\right)
+ =
+ \boldsymbol{F}_{0}\boldsymbol{c}_{p}(T_1 - \boldsymbol{T}_0)
+ + \dot{Q}_{\text{HE}}
+ \\
+ &
+ \dot{Q}_{\text{HE}}
+ =
+ F_1 \boldsymbol{c}_{p}(T_1 - T_2)
+ \\
+ &
+ \dot{Q}_{\text{HE}}
+ =
+ F_w \boldsymbol{c}_{w,p}(T_{w2} - \boldsymbol{T}_{w1})
+ \\
+ & \displaystyle
+ \dot{Q}_{\text{HE}}
+ =
+ \boldsymbol{U}A \left(
+ \frac{(T_1 - T_{w2}) - (T_2 - \boldsymbol{T}_{w1})}{
+ \ln{(T_1 - T_{w2})} - \ln{(T_{2} - \boldsymbol{T}_{w1})}
+ }
+ \right)
+ \\
+ &
+ 0 \leq V \leq \hat{V}
+ \\
+ &
+ 311 \leq T_1, T_2 \leq 389
+ \\
+ &
+ 301 \leq T_{w2} \leq 355
+ \\
+ &
+ 1 \leq T_1 - T_2
+ \\
+ &
+ 1 \leq T_{w2} - \boldsymbol{T}_{w1}
+ \\
+ &
+ 11.1 \leq T_1 - T_{w2}
+ \\
+ &
+ 11.1 \leq T_2 - \boldsymbol{T}_{w1}
+ \end{array}
+
+in which:
+
+- :math:`V` is the reactor holdup volume during operation (in
+ :math:`\text{m}^3`)
+- :math:`F_w` is the mass flow rate of cooling water (in
+ :math:`\text{kg/hr}`)
+- :math:`T_2` is the temperature to which the recycle stream is cooled
+ (in :math:`\text{K}`)
+- :math:`T_{w2}` is the cooling water return temperature (in
+ :math:`\text{K}`)
+- :math:`\dot{Q}_\text{HE}` is the cooling duty of the heat exchanger
+ (in :math:`\text{kW}`)
+
+The objective function yields the total annualized cost,
+with units of $/yr.
+
+Once the design :math:`(\hat{V}, A)` and operation :math:`(F_1, T_1)` of
+the system are specified, the state variables
+:math:`(V, \dot{Q}_\text{HE}, F_w, T_{w2}, T_2)` are calculated using
+the equality constraints, which comprise a square system of nonlinear
+equations.
+
+
+Implement the Model
+~~~~~~~~~~~~~~~~~~~
+
+We now implement the deterministic model for the reactor-cooler system.
+First, we import Pyomo:
+
+.. code::
+
+ >>> import pyomo.environ as pyo
+
+and write a function for building the
+model (with the variables uninitialized):
+
+.. code::
+
+ >>> def build_model():
+ ... m = pyo.ConcreteModel()
+ ...
+ ... # certain parameters
+ ... m.cA0 = pyo.Param(initialize=32.040)
+ ... m.cA1 = pyo.Param(initialize=0.1 * m.cA0)
+ ... m.EovR = pyo.Param(initialize=555.6)
+ ... m.delHr = pyo.Param(initialize=-23260)
+ ... m.cp = pyo.Param(initialize=167.400)
+ ... m.cwp = pyo.Param(initialize=4.184)
+ ... m.F0 = pyo.Param(initialize=45.36)
+ ... m.T0 = pyo.Param(initialize=333)
+ ... m.Tw1 = pyo.Param(initialize=300)
+ ...
+ ... # uncertain parameters
+ ... m.kR = pyo.Param(initialize=10, mutable=True)
+ ... m.U = pyo.Param(initialize=1635, mutable=True)
+ ...
+ ... # first-stage variables
+ ... m.Vhat = pyo.Var(bounds=(0, 20))
+ ... m.A = pyo.Var(bounds=(0, 10))
+ ...
+ ... # second-stage variables
+ ... m.F1 = pyo.Var()
+ ... m.T1 = pyo.Var(bounds=(311, 389))
+ ...
+ ... # state variables
+ ... m.V = pyo.Var(bounds=(0, None))
+ ... m.Qhe = pyo.Var()
+ ... m.T2 = pyo.Var(bounds=(311, 389))
+ ... m.Tw2 = pyo.Var(bounds=(301, 355))
+ ... m.Fw = pyo.Var()
+ ...
+ ... # objective and constituent expressions
+ ... m.capex = pyo.Expression(expr=691.2 * m.Vhat ** 0.7 + 873.6 * m.A ** 0.6)
+ ... m.opex = pyo.Expression(expr=1.76 * m.Fw + 7.056 * m.F1)
+ ... m.obj = pyo.Objective(expr=m.capex + m.opex)
+ ...
+ ... # equality constraints
+ ... m.reactant_mol_bal = pyo.Constraint(
+ ... expr=(
+ ... m.F0 * ((m.cA0 - m.cA1) / m.cA0)
+ ... == m.V * m.kR * pyo.exp(-m.EovR / m.T1) * m.cA1
+ ... ),
+ ... )
+ ... m.reactant_heat_bal = pyo.Constraint(
+ ... expr=(
+ ... -m.delHr * m.F0 * ((m.cA0 - m.cA1) / m.cA0)
+ ... == m.F0 * m.cp * (m.T1 - m.T0)
+ ... + m.Qhe
+ ... )
+ ... )
+ ... m.heat_bal_process = pyo.Constraint(
+ ... expr=m.Qhe == m.F1 * m.cp * (m.T1 - m.T2)
+ ... )
+ ... m.heat_bal_util = pyo.Constraint(
+ ... expr=m.Qhe == m.Fw * m.cwp * (m.Tw2 - m.Tw1)
+ ... )
+ ...
+ ... @m.Constraint()
+ ... def hex_design_eq(mdl):
+ ... dt1 = mdl.T1 - mdl.Tw2
+ ... dt2 = mdl.T2 - mdl.Tw1
+ ... lmtd_expr = (dt1 - dt2) / (pyo.log(dt1) - pyo.log(dt2))
+ ... return m.Qhe == m.A * m.U * lmtd_expr
+ ...
+ ... # inequality constraints
+ ... m.V_con = pyo.Constraint(expr=(m.V <= m.Vhat))
+ ... m.T1T2_con = pyo.Constraint(expr=(1 <= m.T1 - m.T2))
+ ... m.Tw1Tw2_con = pyo.Constraint(expr=(1 <= m.Tw2 - m.Tw1))
+ ... m.T1Tw2_con = pyo.Constraint(expr=(11.1 <= m.T1 - m.Tw2))
+ ... m.T2Tw1_con = pyo.Constraint(expr=(11.1 <= m.T2 - m.Tw1))
+ ...
+ ... return m
+ ...
+
+
+.. note::
+
+ Observe that the :class:`~pyomo.core.base.param.Param`
+ objects representing the potentially uncertain parameters
+ :math:`\boldsymbol{k}_\text{R}` and :math:`\boldsymbol{U}`
+ are declared with the argument ``mutable=True``,
+ as PyROS requires that :class:`~pyomo.core.base.param.Param`
+ objects representing uncertain parameters
+ be mutable.
+ Alternatively,
+ :math:`\boldsymbol{k}_\text{R}` and :math:`\boldsymbol{U}`
+ may have instead been implemented as fixed
+ :class:`~pyomo.core.base.var.Var` objects,
+ as follows:
+
+ .. code-block::
+
+ m.kR = pyo.Var(initialize=10)
+ m.U = pyo.Var(initialize=1635)
+ m.kR.fix(); m.U.fix()
+
+ For more information on implementing uncertain parameters for PyROS,
+ see the
+ :ref:`Uncertain Parameters section of the Solver Interface
+ documentation `.
+
+
+For convenience, we also write a function to initialize the model's
+variables values:
+
+.. code::
+
+ >>> from pyomo.util.calc_var_value import calculate_variable_from_constraint
+ >>>
+ >>> def initialize_model(m, Vhat=20, A=10, F1=50, T1=389):
+ ... # set first-stage and second-stage variable values
+ ... m.Vhat.set_value(Vhat)
+ ... m.A.set_value(A)
+ ... m.F1.set_value(F1)
+ ... m.T1.set_value(T1)
+ ...
+ ... # solve equations for state variables
+ ... calculate_variable_from_constraint(
+ ... variable=m.V,
+ ... constraint=m.reactant_mol_bal,
+ ... )
+ ... calculate_variable_from_constraint(
+ ... variable=m.Qhe,
+ ... constraint=m.reactant_heat_bal,
+ ... )
+ ... calculate_variable_from_constraint(
+ ... variable=m.T2,
+ ... constraint=m.heat_bal_process,
+ ... )
+ ... calculate_variable_from_constraint(
+ ... variable=m.Tw2,
+ ... constraint=m.hex_design_eq,
+ ... )
+ ... calculate_variable_from_constraint(
+ ... variable=m.Fw,
+ ... constraint=m.heat_bal_util,
+ ... )
+ ...
+
+And finally, a function to build the model and
+and initialize the variable values:
+
+.. code::
+
+ >>> def build_and_initialize_model(**init_kwargs):
+ ... m = build_model()
+ ... initialize_model(m, **init_kwargs)
+ ... return m
+ ...
+
+We may now instantiate and initialize the model as follows:
+
+.. code::
+
+ >>> m = build_and_initialize_model()
+
+The following helper function will be useful for inspecting
+the current solution:
+
+.. code::
+
+ >>> def print_solution(model):
+ ... print(f"Objective ($/yr) : {pyo.value(model.obj):.2f}")
+ ... print(f"Reactor volume (m^3) : {model.Vhat.value:.2f}")
+ ... print(f"Cooler area (m^2) : {model.A.value:.2f}")
+ ... print(f"F1 (kmol/hr) : {model.F1.value:.2f}")
+ ... print(f"T1 (K) : {model.T1.value:.2f}")
+ ... print(f"Fw (kg/hr) : {model.Fw.value:.2f}")
+ ...
+
+Inspecting the initial model solution:
+
+.. code::
+
+ >>> print_solution(m) # may vary
+ Objective ($/yr) : 13830.89
+ Reactor volume (m^3) : 20.00
+ Cooler area (m^2) : 10.00
+ F1 (kmol/hr) : 50.00
+ T1 (K) : 389.00
+ Fw (kg/hr) : 2484.43
+
+
+Solve the Model Deterministically
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+We use IPOPT to solve the model to local optimality:
+
+.. _pyros_tutorial_nominal_solve:
+
+.. code::
+
+ >>> ipopt = pyo.SolverFactory("ipopt")
+ >>> pyo.assert_optimal_termination(ipopt.solve(m, tee=True, load_solutions=True))
+ Ipopt ...
+ ...
+ EXIT: Optimal Solution Found.
+
+
+We are able to solve the model to local optimality. Inspecting the
+solution, we notice reductions in the objective and the main variables
+of interest compared to the initial point used:
+
+.. _pyros_tutorial_inspect_nominal:
+
+.. code::
+
+ >>> print_solution(m) # may vary
+ Objective ($/yr) : 9774.58
+ Reactor volume (m^3) : 5.32
+ Cooler area (m^2) : 7.45
+ F1 (kmol/hr) : 88.32
+ T1 (K) : 389.00
+ Fw (kg/hr) : 2278.57
+
+
+Assess Impact of Parametric Uncertainty
+---------------------------------------
+
+Suppose the reaction rate constant :math:`\boldsymbol{k}_\text{R}` and
+heat transfer coefficient :math:`\boldsymbol{U}` are uncertain. We
+assume that each parameter may deviate from its nominal value by up to
+5%, and that the deviations are independent. Thus, the joint
+realizations of the uncertain parameters are confined to a rectangular
+region, that is, a box.
+
+Given a *fixed* design :math:`(\hat{V}, A)`, we wish to assess whether
+we can guarantee that the operational variables :math:`(F_1, T_1)`, and
+concomitantly, the state
+:math:`(V, \dot{Q}_\text{HE}, T_2, T_{w2}, F_w)`, can be adjusted to a
+feasible solution under any plausible realization of the uncertain
+parameters. This assessment can be carried out with the following
+function:
+
+.. code::
+
+ >>> # module imports needed
+ >>> import numpy as np
+ >>> import matplotlib.pyplot as plt
+ >>> import matplotlib.patches as patches
+ >>>
+ >>> def plot_feasibility(solutions, solver, samples=200, test_set_size=10):
+ ... # seed the random number generator for deterministic sampling
+ ... rng = np.random.default_rng(123456)
+ ...
+ ... # nominal uncertain parameter realizations
+ ... nom_vals = np.array([10, 1635])
+ ...
+ ... # sample points from box uncertainty set of specified test size
+ ... point_samples = np.empty((samples, 2))
+ ... point_samples[0] = nom_vals
+ ... point_samples[1:] = rng.uniform(
+ ... low=nom_vals * (1 - test_set_size / 100),
+ ... high=nom_vals * (1 + test_set_size / 100),
+ ... size=(samples - 1, 2),
+ ... )
+ ...
+ ... costs = np.empty((len(solutions), point_samples.shape[0]), dtype=float)
+ ... mdl = build_model()
+ ... for sol_idx, (size, sol) in enumerate(solutions.items()):
+ ... # fix the first-stage variables
+ ... mdl.Vhat.fix(sol[0])
+ ... mdl.A.fix(sol[1])
+ ...
+ ... for pt_idx, pt in enumerate(point_samples):
+ ... # update parameter realization to sampled point
+ ... mdl.kR.set_value(pt[0])
+ ... mdl.U.set_value(pt[1])
+ ...
+ ... # update the values of the operational variables
+ ... initialize_model(mdl, Vhat=sol[0], A=sol[1])
+ ...
+ ... # try solving the model to inspect for feasibility
+ ... res = solver.solve(mdl, load_solutions=False)
+ ... if pyo.check_optimal_termination(res):
+ ... mdl.solutions.load_from(res)
+ ... costs[sol_idx, pt_idx] = pyo.value(mdl.obj)
+ ... else:
+ ... costs[sol_idx, pt_idx] = np.nan
+ ...
+ ... # now generate the plot(s)
+ ... fig, axs = plt.subplots(
+ ... figsize=(0.5 * (len(solutions) - 1) + 5 * len(solutions), 4),
+ ... ncols=len(solutions),
+ ... squeeze=False,
+ ... sharey=True,
+ ... layout="constrained",
+ ... )
+ ... for sol_idx, (size, ax) in enumerate(zip(solutions, axs[0])):
+ ... # track realizations for which solution feasible
+ ... is_feas = ~np.isnan(costs[sol_idx])
+ ...
+ ... # realizations under which design is feasible
+ ... # (colored by objective)
+ ... plot = ax.scatter(
+ ... point_samples[is_feas][:, 0],
+ ... point_samples[is_feas][:, 1],
+ ... c=costs[sol_idx, is_feas],
+ ... vmin=np.nanmin(costs),
+ ... vmax=np.nanmax(costs),
+ ... cmap="plasma_r",
+ ... marker="o",
+ ... )
+ ... # realizations under which design is infeasible
+ ... ax.scatter(
+ ... point_samples[~is_feas][:, 0],
+ ... point_samples[~is_feas][:, 1],
+ ... color="none",
+ ... edgecolors="black",
+ ... label="infeasible",
+ ... marker="^",
+ ... )
+ ... if size != 0:
+ ... # boundary of the box uncertainty set mapped to the design
+ ... rect = patches.Rectangle(
+ ... nom_vals * (1 - size / 100),
+ ... *tuple(nom_vals * 2 * size / 100),
+ ... facecolor="none",
+ ... edgecolor="black",
+ ... linestyle="dashed",
+ ... label=f"{size}% box set",
+ ... )
+ ... ax.add_patch(rect)
+ ...
+ ... ax.legend(bbox_to_anchor=(1, -0.15), loc="upper right")
+ ... ax.set_xlabel(r"$k_\mathrm{R}$ (per hr)")
+ ... ax.set_ylabel("$U$ (kJ/sqm-h-K)")
+ ...
+ ... is_in_set = np.logical_and(
+ ... np.all(nom_vals * (1 - size / 100) <= point_samples, axis=1),
+ ... np.all(point_samples <= nom_vals * (1 + size / 100), axis=1),
+ ... )
+ ... feas_in_set = np.logical_and(is_feas, is_in_set)
+ ...
+ ... # add plot title summarizing statistics of the results
+ ... ax.set_title(
+ ... f"Solution for {size}% box set\n"
+ ... "Avg ± SD objective "
+ ... f"{costs[sol_idx, is_feas].mean():.2f} ± {costs[sol_idx, is_feas].std():.2f}\n"
+ ... f"Feas. for {feas_in_set.sum()}/{is_in_set.sum()} samples in set\n"
+ ... f"Feas. for {is_feas.sum()}/{len(point_samples)} samples overall"
+ ... )
+ ...
+ ... cbar = fig.colorbar(plot, ax=axs.ravel().tolist(), pad=0.03)
+ ... cbar.ax.set_ylabel("Objective ($/yr)")
+ ...
+ ... plt.show()
+ ...
+
+Applying this function to the design that was deterministically
+optimized subject to the nominal realization of the uncertain
+parameters:
+
+.. code::
+
+ >>> plot_feasibility(
+ ... # design variable values
+ ... solutions={0: (m.Vhat.value, m.A.value)},
+ ... # solver to use for feasibility testing
+ ... solver=ipopt,
+ ... # size of the uncertainty set (percent maximum deviation from nominal)
+ ... test_set_size=5,
+ ... )
+
+.. image:: deterministic_heatmap.png
+ :alt: Reactor-cooler system.
+
+
+Clearly, the nominally optimal :math:`(\hat{V}, A)` is robust
+infeasible, as the operation of the system cannot be feasibly adjusted
+subject to approximately half of the tested scenarios. Observe that
+infeasibility occurs subject to parameter realizations in which the rate
+constant :math:`\boldsymbol{k}_\text{R}` is below its nominal value.
+This suggests that for such realizations, the design
+:math:`(\hat{V}, A)` is not sufficiently large to allow for the 90%
+reactor conversion requirement to be met.
+
+Use PyROS to Obtain Robust Solutions
+------------------------------------
+
+We have just confirmed that the nominally optimal design for the reactor
+cooler system is robust infeasible. Thus, we now use PyROS to optimize
+the design while explicitly accounting for the impact of parametric
+uncertainty.
+
+Import PyROS
+~~~~~~~~~~~~
+
+We will need to import the PyROS module in order to instantiate the
+solver and required arguments:
+
+.. code::
+
+ >>> import pyomo.contrib.pyros as pyros
+
+Construct PyROS Solver Arguments
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+We now construct the required arguments to the PyROS solver.
+A general discussion on all PyROS solver arguments is given in the
+:ref:`Solver Interface documentation `.
+
+Deterministic Model
+^^^^^^^^^^^^^^^^^^^
+
+We have already constructed the deterministic model.
+
+First-stage Variables and Second-Stage Variables
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+As previously discussed, the first-stage variables comprise the
+design variables :math:`(\hat{V}, A)`, whereas the second-stage
+variables comprise the operational decision variables
+:math:`(F_1, T_1)`. PyROS automatically infers the state variables of
+the model by inspecting the active objective and constraint components.
+
+.. code::
+
+ >>> first_stage_variables = [m.A, m.Vhat]
+ >>> second_stage_variables = [m.F1, m.T1]
+
+Uncertain Parameters and Uncertainty Set
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Following from our prior feasibility analysis, we take
+:math:`\boldsymbol{k}_\text{R}` and :math:`\boldsymbol{U}` to be the
+uncertain parameters, confined in value to a box set, such that each
+parameter may deviate from its nominal value by up to 5%.
+Thus, we compile the :class:`~pyomo.core.base.param.Param`
+objects representing
+:math:`\boldsymbol{k}_\text{R}` and :math:`\boldsymbol{U}` into a list
+and represent the uncertainty set with an instance of the PyROS
+:class:`~pyomo.contrib.pyros.uncertainty_sets.BoxSet` class:
+
+.. code::
+
+ >>> uncertain_params = [m.kR, m.U]
+ >>> uncertainty_set = pyros.BoxSet(bounds=[
+ ... [param.value * (1 - 0.05), param.value * (1 + 0.05)] for param in uncertain_params
+ ... ])
+
+Subsolvers
+^^^^^^^^^^
+
+PyROS requires subordinate deterministic NLP optimizers to solve the
+subproblems of its underlying algorithm. At least one local NLP solver
+and one global NLP solver are required. We will use IPOPT (already
+constructed) as the local NLP subsolver and BARON as the global NLP
+subsolver. For subproblems not solved successfully by BARON, we use
+COUENNE as a backup.
+
+.. code::
+
+ >>> # already constructed local subsolver IPOPT.
+ >>> # global subsolvers:
+ >>> baron = pyo.SolverFactory("baron", options={"MaxTime": 10})
+ >>> couenne = pyo.SolverFactory("couenne", options={"max_cpu_time": 10})
+
+
+Invoke PyROS
+~~~~~~~~~~~~
+
+We are now ready to invoke PyROS on our model.
+We do so by instantiating the PyROS solver interface:
+
+.. code::
+
+ >>> pyros_solver = pyo.SolverFactory("pyros")
+
+and invoking the :meth:`~pyomo.contrib.pyros.pyros.PyROS.solve` method:
+
+.. _pyros_tutorial_static_ro_solve:
+
+.. code::
+
+ >>> pyros_solver.solve(
+ ... # mandatory arguments
+ ... model=m,
+ ... first_stage_variables=first_stage_variables,
+ ... second_stage_variables=second_stage_variables,
+ ... uncertain_params=uncertain_params,
+ ... uncertainty_set=uncertainty_set,
+ ... local_solver=ipopt,
+ ... global_solver=baron,
+ ... # optional arguments
+ ... backup_global_solvers=[couenne],
+ ... )
+ ==============================================================================
+ PyROS: The Pyomo Robust Optimization Solver, ...
+ ...
+ Robust feasible solution identified.
+ ...
+ All done. Exiting PyROS.
+ ==============================================================================
+
+
+
+By default, the progress and final result of the PyROS solver
+is logged to the console.
+The :ref:`Solver Output Log documentation `
+provides guidance on how the solver log is to be interpreted.
+The :meth:`~pyomo.contrib.pyros.pyros.PyROS.solve` method
+returns an :class:`~pyomo.contrib.pyros.solve_data.ROSolveResults`
+object summarizing the results.
+
+
+Inspect the Solution
+~~~~~~~~~~~~~~~~~~~~
+
+Inspecting the solution, we see that the overall objective is increased
+compared to when the model is
+:ref:`solved deterministically `.
+The cooler area
+:math:`A` and recycle stream flow :math:`F_1` are reduced,
+but the reactor volume :math:`\hat{V}`
+and cooling water utility flow rate :math:`F_w`
+are increased:
+
+.. _pyros_tutorial_inspect_static:
+
+.. code::
+
+ >>> print_solution(m) # may vary
+ Objective ($/yr) : 10064.11
+ Reactor volume (m^3) : 5.59
+ Cooler area (m^2) : 7.19
+ F1 (kmol/hr) : 85.39
+ T1 (K) : 389.00
+ Fw (kg/hr) : 2444.42
+
+
+We can also confirm the robust feasibility of the solution empirically:
+
+.. code::
+
+ >>> plot_feasibility({5: (m.Vhat.value, m.A.value)}, solver=ipopt, test_set_size=5)
+
+
+.. image:: dr0_heatmap.png
+
+
+
+Try Higher-Order Decision Rules to Improve Solution Quality
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+For tractability purposes, the underlying algorithm of PyROS uses
+polynomial decision rules to approximate (restrict) the adjustability of
+the second-stage decision variables (that is, :math:`F_1` and :math:`T_1`
+for the present model) to the uncertain parameters. By default, a static
+approximation is used, such that the second-stage decisions are modeled
+as nonadjustable. A less restrictive approximation can be used by
+increasing the order of the decision rules to 1, through the optional
+argument ``decision_rule_order``:
+
+.. code::
+
+ >>> pyros_solver.solve(
+ ... # mandatory arguments
+ ... model=m,
+ ... first_stage_variables=first_stage_variables,
+ ... second_stage_variables=second_stage_variables,
+ ... uncertain_params=uncertain_params,
+ ... uncertainty_set=uncertainty_set,
+ ... local_solver=ipopt,
+ ... global_solver=baron,
+ ... # optional arguments
+ ... backup_global_solvers=[couenne],
+ ... decision_rule_order=1,
+ ... )
+ ==============================================================================
+ PyROS: The Pyomo Robust Optimization Solver, ...
+ ...
+ Robust feasible solution identified.
+ ...
+ All done. Exiting PyROS.
+ ==============================================================================
+
+
+
+Inspecting the solution, we see that the cost is reduced compared to when
+:ref:`a static decision rule approximation is used `,
+as a smaller cooling water flow rate :math:`F_w` is required
+since the cooler area :math:`A` is increased:
+
+.. code::
+
+ >>> print_solution(m) # may vary
+ Objective ($/yr) : 9855.95
+ Reactor volume (m^3) : 5.59
+ Cooler area (m^2) : 7.45
+ F1 (kmol/hr) : 88.32
+ T1 (K) : 389.00
+ Fw (kg/hr) : 2278.57
+
+
+Further, our empirical check confirms that the solution
+is robust:
+
+.. code::
+
+ >>> plot_feasibility({5: (m.Vhat.value, m.A.value)}, solver=ipopt, test_set_size=5)
+
+
+.. image:: dr1_heatmap.png
+
+
+Assess Impact of Uncertainty Set on the Solution Obtained
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+We now perform a price-of-robustness study, in which
+we assess the response of the solution obtained to the
+size of the uncertainty set.
+This study can be easily performed by placing the PyROS solver
+invocation in a ``for`` loop and recording the result at each iteration.
+
+Invoke PyROS in a ``for`` Loop
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The PyROS solver invocation can easily be made in a ``for`` loop. At
+each iteration of the loop, we use PyROS to solve the RO problem subject
+to the uncertainty set of the corresponding size:
+
+.. code::
+
+ >>> res_dict = dict()
+ >>> obj_vals = dict()
+ >>> capex_vals = dict()
+ >>> opex_vals = dict()
+ >>> vhat_vals = dict()
+ >>> area_vals = dict()
+ >>> for percent_size in [0, 2.5, 5, 7.5, 10]:
+ ... mdl = build_and_initialize_model()
+ ... unc_set = pyros.BoxSet(bounds=[
+ ... [param.value * (1 - percent_size / 100), param.value * (1 + percent_size / 100)]
+ ... for param in [mdl.kR, mdl.U]
+ ... ])
+ ... print(f"Solving RO problem for uncertainty set size {percent_size}:")
+ ... res_dict[percent_size] = res = pyros_solver.solve(
+ ... model=mdl,
+ ... first_stage_variables=[mdl.Vhat, mdl.A],
+ ... second_stage_variables=[mdl.F1, mdl.T1],
+ ... uncertain_params=[mdl.kR, mdl.U],
+ ... uncertainty_set=unc_set,
+ ... local_solver=ipopt,
+ ... global_solver=baron,
+ ... decision_rule_order=1,
+ ... backup_global_solvers=[couenne],
+ ... )
+ ... if res.pyros_termination_condition == pyros.pyrosTerminationCondition.robust_feasible:
+ ... obj_vals[percent_size] = pyo.value(mdl.obj)
+ ... capex_vals[percent_size] = pyo.value(mdl.capex)
+ ... opex_vals[percent_size] = pyo.value(mdl.opex)
+ ... vhat_vals[percent_size] = pyo.value(mdl.Vhat)
+ ... area_vals[percent_size] = pyo.value(mdl.A)
+ ...
+ Solving RO problem for uncertainty set size 0:
+ ...
+ Solving RO problem for uncertainty set size 2.5:
+ ...
+ Solving RO problem for uncertainty set size 5:
+ ...
+ Solving RO problem for uncertainty set size 7.5:
+ ...
+ Solving RO problem for uncertainty set size 10:
+ ...
+ All done. Exiting PyROS.
+ ==============================================================================
+
+
+Visualize the Results
+^^^^^^^^^^^^^^^^^^^^^
+
+We can visualize the results of our price-of-robustness analysis, as follows:
+
+.. code::
+
+ >>> fig, (obj_ax, vhat_ax, area_ax) = plt.subplots(
+ ... ncols=3,
+ ... figsize=(19, 4),
+ ... layout="constrained",
+ ... )
+ >>> plt.subplots_adjust(wspace=0.4, hspace=0.6)
+ >>>
+ >>> # plot costs
+ >>> obj_ax.plot(obj_vals.keys(), obj_vals.values(), label="total", marker="o")
+ >>> obj_ax.plot(capex_vals.keys(), capex_vals.values(), label="CAPEX", marker="s")
+ >>> obj_ax.plot(opex_vals.keys(), opex_vals.values(), label="OPEX", marker="^")
+ >>> obj_ax.set_xlabel("Deviation from Nominal Value (%)")
+ >>> obj_ax.set_ylabel("Cost ($/yr)")
+ >>> obj_ax.legend()
+ >>>
+ >>> # plot reactor volume
+ >>> vhat_ax.plot(vhat_vals.keys(), vhat_vals.values(), marker="o")
+ >>> vhat_ax.set_xlabel("Deviation from Nominal Value (%)")
+ >>> vhat_ax.set_ylabel(r"Reactor Volume ($\mathrm{m}^3$)")
+ >>>
+ >>> # plot cooler area
+ >>> area_ax.plot(area_vals.keys(), area_vals.values(), marker="o")
+ >>> area_ax.set_xlabel("Deviation from Nominal Value (%)")
+ >>> area_ax.set_ylabel(r"Cooler Heat Transfer Area ($\mathrm{m}^2$)")
+ >>> area_ax.set_ylim([7.45, 7.46])
+ >>>
+ >>> plt.show()
+
+
+.. image:: por_sensitivity.png
+
+Notice that the costs and reactor volume increase with the size
+of the uncertainty set, whereas the heat transfer area of the
+cooler does not vary.
+
+
+Assess Robust Feasibility of the Solutions
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+We can also visualize the robustness of each solution:
+
+.. code::
+
+ >>> plot_feasibility(
+ ... {key: (vhat_vals[key], area_vals[key]) for key in vhat_vals},
+ ... solver=ipopt,
+ ... test_set_size=10,
+ ... )
+
+
+
+.. image:: por_heatmaps.png
+
+
+Notice that:
+
+- Every solution is found to be robust feasible subject to its corresponding
+ uncertainty set, but robust infeasible subject to strict supersets.
+- As the size of uncertainty set is increased, so is the average objective
+ value.
diff --git a/doc/OnlineDocs/explanation/solvers/pyrosdir/uncertainty_sets.rst b/doc/OnlineDocs/explanation/solvers/pyrosdir/uncertainty_sets.rst
new file mode 100644
index 00000000000..1cba49f065e
--- /dev/null
+++ b/doc/OnlineDocs/explanation/solvers/pyrosdir/uncertainty_sets.rst
@@ -0,0 +1,108 @@
+.. _pyros_uncertainty_sets:
+
+======================
+PyROS Uncertainty Sets
+======================
+
+.. contents:: Table of Contents
+ :depth: 1
+ :local:
+
+
+Overview
+========
+In PyROS, the uncertainty set of a robust optimization problem
+is represented by an instance of a subclass of the
+:class:`~pyomo.contrib.pyros.uncertainty_sets.UncertaintySet`
+abstract base class.
+PyROS provides a suite of
+:ref:`pre-implemented concrete subclasses `
+to facilitate instantiation of uncertainty sets
+that are commonly used in the optimization literature.
+:ref:`Custom uncertainty set types `
+can be implemented by subclassing
+:class:`~pyomo.contrib.pyros.uncertainty_sets.UncertaintySet`.
+
+
+.. note::
+ The :class:`~pyomo.contrib.pyros.uncertainty_sets.UncertaintySet`
+ class is an abstract class and therefore cannot be
+ directly instantiated.
+
+
+.. _pyros_pre_implemented_types:
+
+Pre-Implemented Uncertainty Set Types
+=====================================
+The pre-implemented
+:class:`~pyomo.contrib.pyros.uncertainty_sets.UncertaintySet`
+subclasses are enumerated below:
+
+.. autosummary::
+
+ ~pyomo.contrib.pyros.uncertainty_sets.AxisAlignedEllipsoidalSet
+ ~pyomo.contrib.pyros.uncertainty_sets.BoxSet
+ ~pyomo.contrib.pyros.uncertainty_sets.BudgetSet
+ ~pyomo.contrib.pyros.uncertainty_sets.CardinalitySet
+ ~pyomo.contrib.pyros.uncertainty_sets.DiscreteScenarioSet
+ ~pyomo.contrib.pyros.uncertainty_sets.EllipsoidalSet
+ ~pyomo.contrib.pyros.uncertainty_sets.FactorModelSet
+ ~pyomo.contrib.pyros.uncertainty_sets.IntersectionSet
+ ~pyomo.contrib.pyros.uncertainty_sets.PolyhedralSet
+
+
+Mathematical definitions of the pre-implemented
+:class:`~pyomo.contrib.pyros.uncertainty_sets.UncertaintySet`
+subclasses are provided below.
+
+.. _pyros_uncertainty_sets_math_defs:
+
+.. list-table:: Mathematical definitions of PyROS uncertainty sets of dimension :math:`n`.
+ :header-rows: 1
+ :class: scrollwide
+
+ * - Uncertainty Set Type
+ - Input Data
+ - Mathematical Definition
+ * - :class:`~pyomo.contrib.pyros.uncertainty_sets.AxisAlignedEllipsoidalSet`
+ - :math:`\begin{array}{l} q^0 \in \mathbb{R}^{n}, \\ \alpha \in \mathbb{R}_{+}^{n} \end{array}`
+ - :math:`\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\}`
+ * - :class:`~pyomo.contrib.pyros.uncertainty_sets.BoxSet`
+ - :math:`\begin{array}{l} q ^{\text{L}} \in \mathbb{R}^{n}, \\ q^{\text{U}} \in \mathbb{R}^{n} \end{array}`
+ - :math:`\{q \in \mathbb{R}^n \mid q^\mathrm{L} \leq q \leq q^\mathrm{U}\}`
+ * - :class:`~pyomo.contrib.pyros.uncertainty_sets.BudgetSet`
+ - :math:`\begin{array}{l} q^{0} \in \mathbb{R}^{n}, \\ b \in \mathbb{R}_{+}^{L}, \\ B \in \{0, 1\}^{L \times n} \end{array}`
+ - :math:`\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\}`
+ * - :class:`~pyomo.contrib.pyros.uncertainty_sets.CardinalitySet`
+ - :math:`\begin{array}{l} q^{0} \in \mathbb{R}^{n}, \\ \hat{q} \in \mathbb{R}_{+}^{n}, \\ \Gamma \in [0, n] \end{array}`
+ - :math:`\left\{ q \in \mathbb{R}^{n} \middle| \begin{array}{l} \exists\,\xi \in [0, 1]^n\,:\\ \quad \,q = q^{0} + \hat{q} \circ \xi \\ \quad \displaystyle \sum_{i=1}^{n} \xi_{i} \leq \Gamma \end{array} \right\}`
+ * - :class:`~pyomo.contrib.pyros.uncertainty_sets.DiscreteScenarioSet`
+ - :math:`q^{1}, q^{2},\dots , q^{S} \in \mathbb{R}^{n}`
+ - :math:`\{q^{1}, q^{2}, \dots , q^{S}\}`
+ * - :class:`~pyomo.contrib.pyros.uncertainty_sets.EllipsoidalSet`
+ - :math:`\begin{array}{l} q^0 \in \mathbb{R}^n, \\ P \in \mathbb{S}_{++}^{n}, \\ s \in \mathbb{R}_{+} \end{array}`
+ - :math:`\{q \in \mathbb{R}^{n} \mid (q - q^{0})^{\intercal} P^{-1} (q - q^{0}) \leq s\}`
+ * - :class:`~pyomo.contrib.pyros.uncertainty_sets.FactorModelSet`
+ - :math:`\begin{array}{l} q^{0} \in \mathbb{R}^{n}, \\ \Psi \in \mathbb{R}^{n \times F}, \\ \beta \in [0, 1] \end{array}`
+ - :math:`\left\{ q \in \mathbb{R}^{n} \middle| \begin{array}{l} \exists\,\xi \in [-1, 1]^F \,:\\ \quad q = q^{0} + \Psi \xi \\ \quad \displaystyle\bigg| \sum_{j=1}^{F} \xi_{j} \bigg| \leq \beta F \end{array} \right\}`
+ * - :class:`~pyomo.contrib.pyros.uncertainty_sets.IntersectionSet`
+ - :math:`\mathcal{Q}_{1}, \mathcal{Q}_{2}, \dots , \mathcal{Q}_{m} \subset \mathbb{R}^{n}`
+ - :math:`\displaystyle \bigcap_{i=1}^{m} \mathcal{Q}_{i}`
+ * - :class:`~pyomo.contrib.pyros.uncertainty_sets.PolyhedralSet`
+ - :math:`\begin{array}{l} A \in \mathbb{R}^{m \times n}, \\ b \in \mathbb{R}^{m}\end{array}`
+ - :math:`\{q \in \mathbb{R}^{n} \mid A q \leq b\}`
+
+
+.. _pyros_custom_sets:
+
+Custom Uncertainty Set Types
+============================
+A custom uncertainty set type
+in which all uncertain parameters
+are modeled as continuous quantities
+can be implemented by subclassing
+:class:`~pyomo.contrib.pyros.uncertainty_sets.UncertaintySet`.
+For discrete sets, we recommend using the pre-implemented
+:class:`~pyomo.contrib.pyros.uncertainty_sets.DiscreteScenarioSet`
+subclass instead of implementing a new set type.
+PyROS does not support mixed-integer uncertainty set types.
diff --git a/pyomo/contrib/pyros/pyros.py b/pyomo/contrib/pyros/pyros.py
index 0205d54a00e..11c3213b249 100644
--- a/pyomo/contrib/pyros/pyros.py
+++ b/pyomo/contrib/pyros/pyros.py
@@ -78,15 +78,25 @@ def _get_pyomo_version_info():
@SolverFactory.register(
"pyros",
- doc="Robust optimization (RO) solver implementing "
- "the generalized robust cutting-set algorithm (GRCS)",
+ doc="Pyomo Robust Optimization Solver (PyROS): "
+ "implementation of a generalized robust cutting-set algorithm (GRCS)",
)
class PyROS(object):
- '''
- PyROS (Pyomo Robust Optimization Solver) implementing a
+ """
+ Pyomo Robust Optimization Solver (PyROS): implementation of a
generalized robust cutting-set algorithm (GRCS)
- to solve two-stage NLP optimization models under uncertainty.
- '''
+ for the solution of two-stage nonlinear programs
+ under uncertainty.
+
+ We recommend instantiating this class as follows:
+
+ .. code::
+
+ >>> import pyomo.environ as pyo
+ >>> import pyomo.contrib.pyros as pyros
+ >>> pyros_solver = pyo.SolverFactory("pyros")
+
+ """
CONFIG = pyros_config()
_LOG_LINE_LENGTH = 78
diff --git a/pyomo/contrib/pyros/uncertainty_sets.py b/pyomo/contrib/pyros/uncertainty_sets.py
index 8ec3cb13486..06a6fdeb816 100644
--- a/pyomo/contrib/pyros/uncertainty_sets.py
+++ b/pyomo/contrib/pyros/uncertainty_sets.py
@@ -477,8 +477,12 @@ class UncertaintySet(object, metaclass=abc.ABCMeta):
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,
+ `UncertaintySet` is an abstract base class for constructing
+ specific subclasses and instances of uncertainty sets.
+ Therefore, `UncertaintySet` cannot be instantiated directly.
+
+ A concrete `UncertaintySet` instance 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.
"""
@@ -1186,13 +1190,31 @@ def dim(self):
class BoxSet(UncertaintySet):
"""
- A hyper-rectangle (i.e., "box").
+ A hyperrectangle (or box).
Parameters
----------
bounds : (N, 2) array_like
Lower and upper bounds for each dimension of the set.
+ Notes
+ -----
+ The :math:`n`-dimensional box set is defined by
+
+ .. math::
+
+ \\left\\{
+ q \\in \\mathbb{R}^n\\,|
+ \\,q^\\text{L} \\leq q \\leq q^\\text{U}
+ \\right\\}
+
+ in which
+ :math:`q^\\text{L} \\in \\mathbb{R}^n` refers to
+ ``numpy.array(bounds)[:, 0]``,
+ and
+ :math:`q^\\text{U} \\in \\mathbb{R}^n` refers to
+ ``numpy.array(bounds)[:, 1]``.
+
Examples
--------
1D box set (interval):
@@ -1211,7 +1233,7 @@ class BoxSet(UncertaintySet):
5D hypercube with bounds 0 and 1 in each dimension:
- >>> hypercube_5d = BoxSet(bounds=[[0, 1] for idx in range(5)])
+ >>> hypercube_5d = BoxSet(bounds=[[0, 1]] * 5)
>>> hypercube_5d.bounds
array([[0, 1],
[0, 1],
@@ -1363,6 +1385,29 @@ class CardinalitySet(UncertaintySet):
may realize their maximal deviations from the origin
simultaneously.
+ Notes
+ -----
+ The :math:`n`-dimensional cardinality set is defined by
+
+ .. math::
+
+ \\left\\{ q \\in \\mathbb{R}^n\\,\\middle|
+ \\,\\exists\\, \\xi \\in [0, 1]^n \\,:\\,
+ \\left[
+ \\begin{array}{l}
+ q = q^0 + \\hat{q} \\circ \\xi \\\\
+ \\displaystyle \\sum_{i=1}^n \\xi_i \\leq \\Gamma
+ \\end{array}
+ \\right]
+ \\right\\}
+
+ in which
+ :math:`q^\\text{0} \\in \\mathbb{R}^n` refers to ``origin``,
+ the quantity :math:`\\hat{q} \\in \\mathbb{R}_{+}^n`
+ refers to ``positive_deviation``,
+ and :math:`\\Gamma \\in [0, n]` refers to ``gamma``.
+ The operator ":math:`\\circ`" denotes the element-wise product.
+
Examples
--------
A 3D cardinality set:
@@ -1430,7 +1475,7 @@ def origin(self, val):
def positive_deviation(self):
"""
(N,) numpy.ndarray : Maximal coordinate deviations from the
- origin in each dimension. All entries are nonnegative.
+ origin in each dimension. All entries should be nonnegative.
"""
return self._positive_deviation
@@ -1461,15 +1506,16 @@ def positive_deviation(self, val):
def gamma(self):
"""
numeric type : Upper bound for the number of uncertain
- parameters which may maximally deviate from their respective
+ parameters that 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
+ to a singleton containing the point represented by
+ ``self.origin``, 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``.
+ ``numpy.array([self.origin, self.origin + self.positive_deviation]).T``.
"""
return self._gamma
@@ -1662,6 +1708,23 @@ class PolyhedralSet(UncertaintySet):
``lhs_coefficients_mat @ x``, where `x` is an (N,)
array representing any point in the polyhedral set.
+ Notes
+ -----
+ The :math:`n`-dimensional polyhedral set is defined by
+
+ .. math::
+
+ \\left\\{
+ q \\in \\mathbb{R}^n\\,
+ \\middle| \\, A q \\leq b
+ \\right\\}
+
+ in which
+ :math:`A \\in \\mathbb{R}^{m \\times n}` refers to
+ ``lhs_coefficients_mat``,
+ and
+ :math:`b \\in \\mathbb{R}^m` refers to ``rhs_vec``.
+
Examples
--------
2D polyhedral set with 4 defining inequalities:
@@ -1875,6 +1938,7 @@ class BudgetSet(UncertaintySet):
Each row corresponds to a single budget constraint,
and defines which uncertain parameters
(which dimensions) participate in that row's constraint.
+ All entries should be of value 0 or 1.
rhs_vec : (L,) array_like
Budget limits (upper bounds) with respect to
the origin of the set.
@@ -1882,6 +1946,28 @@ class BudgetSet(UncertaintySet):
Origin of the budget set. If `None` is provided, then
the origin is set to the zero vector.
+ Notes
+ -----
+ The :math:`n`-dimensional budget set is defined by
+
+ .. math::
+
+ \\left\\{
+ q \\in \\mathbb{R}^n\\,\\middle|
+ \\begin{pmatrix} B \\\\ -I \\end{pmatrix} q
+ \\leq \\begin{pmatrix} b + Bq^0 \\\\ -q^0 \\end{pmatrix}
+ \\right\\}
+
+ in which
+ :math:`B \\in \\{0, 1\\}^{\\ell \\times n}` refers to
+ ``budget_membership_mat``,
+ the quantity
+ :math:`I` denotes the :math:`n \\times n` identity matrix,
+ the quantity
+ :math:`b \\in \\mathbb{R}_{+}^\\ell` refers to ``rhs_vec``,
+ and
+ :math:`q^0 \\in \\mathbb{R}^n` refers to ``origin``.
+
Examples
--------
3D budget set with one budget constraint and
@@ -1966,6 +2052,7 @@ def budget_membership_mat(self):
constraints. Each row corresponds to a single budget
constraint and defines which uncertain parameters
participate in that row's constraint.
+ All entries should be of value 0 or 1.
"""
return self._budget_membership_mat
@@ -2220,6 +2307,34 @@ class FactorModelSet(UncertaintySet):
independent factors that can simultaneously attain
their extreme values.
+ Notes
+ -----
+ The :math:`n`-dimensional factor model set is defined by
+
+ .. math::
+
+ \\left\\{ q \\in \\mathbb{R}^n\\,
+ \\middle|
+ \\,
+ \\exists\\, \\xi \\in [-1, 1]^F \\,:\\,
+ \\left[
+ \\begin{array}{l}
+ q = q^0 + \\Psi \\xi \\\\
+ \\displaystyle
+ \\bigg| \\sum_{i=1}^n \\xi_i \\bigg|
+ \\leq \\beta F
+ \\end{array}
+ \\right]
+ \\right\\}
+
+ in which
+ :math:`q^\\text{0} \\in \\mathbb{R}^n` refers to ``origin``,
+ the quantity :math:`F` refers to ``number_of_factors``,
+ the quantity :math:`\\Psi \\in \\mathbb{R}^{n \\times F}`
+ refers to ``psi_mat``,
+ and :math:`\\beta \\in [0, 1]` refers to ``beta``.
+
+
Examples
--------
A 4D factor model set with a 2D factor space:
@@ -2367,8 +2482,8 @@ def beta(self):
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
- ``[self.origin - psi @ np.ones(F), self.origin + psi @ np.ones(F)].T``.
+ then any number of factors can be above 0 or below 0
+ simultaneously.
"""
return self._beta
@@ -2582,7 +2697,7 @@ def validate(self, config):
class AxisAlignedEllipsoidalSet(UncertaintySet):
"""
- An axis-aligned ellipsoid.
+ An axis-aligned ellipsoidal region.
Parameters
----------
@@ -2591,18 +2706,46 @@ class AxisAlignedEllipsoidalSet(UncertaintySet):
half_lengths : (N,) array_like
Semi-axis lengths of the ellipsoid.
+ See Also
+ --------
+ EllipsoidalSet : A general ellipsoidal region.
+
+ Notes
+ -----
+ The :math:`n`-dimensional axis-aligned ellipsoidal set is defined by
+
+ .. math::
+
+ \\left\\{
+ q \\in \\mathbb{R}^n\\,
+ \\middle|\\,
+ \\begin{array}{l}
+ \\displaystyle
+ \\sum_{\\substack{i = 1 \\\\ \\alpha_i > 0}}^n
+ \\bigg( \\frac{q_i - q_i^0}{\\alpha_i}\\bigg)^2
+ \\leq 1
+ \\\\
+ q_i = q_i^0\\quad\\forall\\,i \\,:\\,\\alpha_i = 0
+ \\end{array}
+ \\right\\}
+
+ in which
+ :math:`q^0 \\in \\mathbb{R}^n` refers to ``center``,
+ and :math:`\\alpha \\in \\mathbb{R}_{+}^n` refers to
+ ``half_lengths``.
+
Examples
--------
- 3D origin-centered unit hypersphere:
+ 3D origin-centered unit ball:
>>> from pyomo.contrib.pyros import AxisAlignedEllipsoidalSet
- >>> sphere = AxisAlignedEllipsoidalSet(
+ >>> ball = AxisAlignedEllipsoidalSet(
... center=[0, 0, 0],
- ... half_lengths=[1, 1, 1]
+ ... half_lengths=[1, 1, 1],
... )
- >>> sphere.center
+ >>> ball.center
array([0, 0, 0])
- >>> sphere.half_lengths
+ >>> ball.half_lengths
array([1, 1, 1])
"""
@@ -2794,7 +2937,7 @@ def validate(self, config):
class EllipsoidalSet(UncertaintySet):
"""
- A general ellipsoid.
+ A general ellipsoidal region.
Parameters
----------
@@ -2814,6 +2957,34 @@ class EllipsoidalSet(UncertaintySet):
Exactly one of `scale` and `gaussian_conf_lvl` should be
None; otherwise, an exception is raised.
+ See Also
+ --------
+ AxisAlignedEllipsoidalSet : An axis-aligned ellipsoidal region.
+
+ Notes
+ -----
+ The :math:`n`-dimensional ellipsoidal set is defined by
+
+ .. math::
+
+ \\left\\{
+ q \\in \\mathbb{R}^n\\,|
+ \\,(q - q^0)^\\intercal \\Sigma^{-1}(q - q^0) \\leq s
+ \\right\\}
+
+ in which
+ :math:`q^0 \\in \\mathbb{R}^n` refers to ``center``,
+ the quantity
+ :math:`\\Sigma \\in \\mathbb{R}^{n \\times n}`
+ refers to ``shape_matrix``,
+ and :math:`s \\geq 0` refers to ``scale``.
+
+ The quantity :math:`s` is related to the Gaussian confidence level
+ (``gaussian_conf_lvl``) :math:`p \\in [0, 1)`
+ by :math:`s = \\chi_{n}^2(p)`, in which
+ :math:`\\chi_{n}^2(\\cdot)` is the quantile function
+ of the chi-squared distribution with :math:`n` degrees of freedom.
+
Examples
--------
A 3D origin-centered unit ball:
@@ -2834,7 +3005,7 @@ class EllipsoidalSet(UncertaintySet):
>>> ball.scale
1
- A 2D ellipsoid with custom rotation and scaling:
+ A 2D ellipsoidal region with custom rotation and scaling:
>>> rotated_ellipsoid = EllipsoidalSet(
... center=[1, 1],
@@ -2849,7 +3020,7 @@ class EllipsoidalSet(UncertaintySet):
>>> rotated_ellipsoid.scale
0.5
- A 4D 95% confidence ellipsoid:
+ A 4D 95% confidence ellipsoidal region:
>>> conf_ellipsoid = EllipsoidalSet(
... center=np.zeros(4),
@@ -3197,14 +3368,26 @@ def validate(self, config):
class DiscreteScenarioSet(UncertaintySet):
"""
- A discrete set of finitely many uncertain parameter realizations
- (or scenarios).
+ A set of finitely many distinct points (or scenarios).
Parameters
----------
scenarios : (M, N) array_like
A sequence of `M` distinct uncertain parameter realizations.
+ Notes
+ -----
+ The :math:`n`-dimensional discrete set is defined by
+
+ .. math::
+
+ \\left\\{
+ q^1, q^2, \\dots, q^m
+ \\right\\}
+
+ in which :math:`q^i \\in \\mathbb{R}^n`
+ refers to ``scenarios[i - 1]`` for :math:`i = 1, 2, \\dots, m`.
+
Examples
--------
2D set with three scenarios:
@@ -3388,7 +3571,7 @@ def validate(self, config):
class IntersectionSet(UncertaintySet):
"""
- An intersection of a sequence of uncertainty sets, each of which
+ An intersection of two or more uncertainty sets, each of which
is represented by an `UncertaintySet` object.
Parameters
@@ -3398,6 +3581,19 @@ class IntersectionSet(UncertaintySet):
an intersection. At least two uncertainty sets must
be provided. All sets must be of the same dimension.
+ Notes
+ -----
+ The :math:`n`-dimensional intersection set is defined by
+
+ .. math::
+
+ \\mathcal{Q}_1 \\cap \\mathcal{Q}_2 \\cap \\cdots
+ \\cap \\mathcal{Q}_m
+
+ in which :math:`\\mathcal{Q}_i \\subset \\mathbb{R}^n`
+ refers to the uncertainty set ``list(unc_sets.values())[i - 1]``
+ for :math:`i = 1, 2, \\dots, m`.
+
Examples
--------
Intersection of origin-centered 2D box (square) and 2D
@@ -3411,7 +3607,8 @@ class IntersectionSet(UncertaintySet):
... center=[0, 0],
... half_lengths=[2, 2],
... )
- >>> # to construct intersection, pass sets as keyword arguments
+ >>> # to construct intersection, pass sets as keyword arguments.
+ >>> # keywords are arbitrary
>>> intersection = IntersectionSet(set1=square, set2=circle)
>>> intersection.all_sets # doctest: +ELLIPSIS
UncertaintySetList([...])