Skip to content

Commit

Permalink
Merge branch 'main' into incremental
Browse files Browse the repository at this point in the history
  • Loading branch information
emma58 authored Aug 19, 2024
2 parents 233b016 + ccc96bf commit 2a6213d
Show file tree
Hide file tree
Showing 123 changed files with 10,923 additions and 8,280 deletions.
107 changes: 107 additions & 0 deletions doc/OnlineDocs/contributed_packages/alternative_solutions.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,107 @@
###############################################
Generating Alternative (Near-)Optimal Solutions
###############################################

Optimization solvers are generally designed to return a feasible solution
to the user. However, there are many applications where a user needs
more context than this result. For example,

* alternative solutions can support an assessment of trade-offs between competing objectives;

* if the optimization formulation may be inaccurate or untrustworthy, then comparisons amongst alternative solutions provide additional insights into the reliability of these model predictions; or

* the user may have unexpressed objectives or constraints, which only are realized in later stages of model analysis.

The *alternative-solutions library* provides a variety of functions that
can be used to generate optimal or near-optimal solutions for a pyomo
model. Conceptually, these functions are like pyomo solvers. They can
be configured with solver names and options, and they return a list of
solutions for the pyomo model. However, these functions are independent
of pyomo's solver interface because they return a custom solution object.

The following functions are defined in the alternative-solutions library:

* ``enumerate_binary_solutions``

* Finds alternative optimal solutions for a binary problem using no-good cuts.

* ``enumerate_linear_solutions``

* Finds alternative optimal solutions for a (mixed-integer) linear program.

* ``enumerate_linear_solutions_soln_pool``

* Finds alternative optimal solutions for a (mixed-binary) linear program using Gurobi's solution pool feature.

* ``gurobi_generate_solutions``

* Finds alternative optimal solutions for discrete variables using Gurobi's built-in solution pool capability.

* ``obbt_analysis_bounds_and_solutions``

* Calculates the bounds on each variable by solving a series of min and max optimization problems where each variable is used as the objective function. This can be applied to any class of problem supported by the selected solver.


Usage Example
-------------

Many of functions in the alternative-solutions library have similar options, so we simply illustrate the ``enumerate_binary_solutions`` function. We define a simple knapsack example whose alternative solutions have integer objective values ranging from 0 to 90.

.. doctest::

>>> import pyomo.environ as pyo

>>> values = [10, 40, 30, 50]
>>> weights = [5, 4, 6, 3]
>>> capacity = 10

>>> m = pyo.ConcreteModel()
>>> m.x = pyo.Var(range(4), within=pyo.Binary)
>>> m.o = pyo.Objective(expr=sum(values[i] * m.x[i] for i in range(4)), sense=pyo.maximize)
>>> m.c = pyo.Constraint(expr=sum(weights[i] * m.x[i] for i in range(4)) <= capacity)

We can execute the ``enumerate_binary_solutions`` function to generate a list of ``Solution`` objects that represent alternative optimal solutions:

.. doctest::
:skipif: not glpk_available

>>> import pyomo.contrib.alternative_solutions as aos
>>> solns = aos.enumerate_binary_solutions(m, num_solutions=100, solver="glpk")
>>> assert len(solns) == 10

Each ``Solution`` object contains information about the objective and variables, and it includes various methods to access this information. For example:

.. doctest::
:skipif: not glpk_available

>>> print(solns[0])
{
"fixed_variables": [],
"objective": "o",
"objective_value": 90.0,
"solution": {
"x[0]": 0,
"x[1]": 1,
"x[2]": 0,
"x[3]": 1
}
}


Interface Documentation
-----------------------

.. currentmodule:: pyomo.contrib.alternative_solutions

.. autofunction:: enumerate_binary_solutions

.. autofunction:: enumerate_linear_solutions

.. autofunction:: enumerate_linear_solutions_soln_pool

.. autofunction:: gurobi_generate_solutions

.. autofunction:: obbt_analysis_bounds_and_solutions

.. autoclass:: Solution

Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
160 changes: 54 additions & 106 deletions doc/OnlineDocs/contributed_packages/doe/doe.rst
Original file line number Diff line number Diff line change
Expand Up @@ -116,64 +116,14 @@ In order to solve problems of the above, Pyomo.DoE implements the 2-stage stocha

Pyomo.DoE Required Inputs
--------------------------------
The required inputs to the Pyomo.DoE solver are the following:
The required input to the Pyomo.DoE solver is an ``Experiment`` object. The experiment object must have a ``get_labeled_model`` function which returns a Pyomo model with four ``Suffix`` components identifying the parts of the model used in MBDoE analysis. This is in line with the convention used in the parameter estimation tool, `Parmest <https://pyomo.readthedocs.io/en/stable/contributed_packages/parmest/index.html>`_. The four ``Suffix`` components are:

* A function that creates the process model
* Dictionary of parameters and their nominal value
* A measurement object
* A design variables object
* A Numpy ``array`` containing the Prior FIM
* Optimization solver

Below is a list of arguments that Pyomo.DoE expects the user to provide.

parameter_dict : ``dictionary``
A ``dictionary`` of parameter names and values. If they are an indexed variable, put the variable name and index in a nested ``Dictionary``.

design_variables: ``DesignVariables``
A ``DesignVariables`` of design variables, provided by the DesignVariables class.
If this design var is independent of time (constant), set the time to [0]

measurement_variables : ``MeasurementVariables``
A ``MeasurementVariables`` of the measurements, provided by the MeasurementVariables class.

create_model : ``function``
A ``function`` returning a deterministic process model.

prior_FIM : ``array``
An ``array`` defining the Fisher information matrix (FIM) for prior experiments, default is a zero matrix.

Pyomo.DoE Solver Interface
---------------------------

.. figure:: uml.png
:scale: 25 %


.. autoclass:: pyomo.contrib.doe.doe.DesignOfExperiments
:members: __init__, stochastic_program, compute_FIM, run_grid_search

.. Note::
``stochastic_program()`` includes the following steps:
#. Build two-stage stochastic programming optimization model where scenarios correspond to finite difference approximations for the Jacobian of the response variables with respect to calibrated model parameters
#. Fix the experiment design decisions and solve a square (i.e., zero degrees of freedom) instance of the two-stage DOE problem. This step is for initialization.
#. Unfix the experiment design decisions and solve the two-stage DOE problem.

.. autoclass:: pyomo.contrib.doe.measurements.MeasurementVariables
:members: __init__, add_variables

.. autoclass:: pyomo.contrib.doe.measurements.DesignVariables
:members: __init__, add_variables

.. autoclass:: pyomo.contrib.doe.scenario.ScenarioGenerator
:special-members: __init__

.. autoclass:: pyomo.contrib.doe.result.FisherResults
:members: __init__, result_analysis

.. autoclass:: pyomo.contrib.doe.result.GridSearchResult
:special-members: __init__
* ``experiment_inputs`` - The experimental design decisions
* ``experiment_outputs`` - The values measured during the experiment
* ``measurement_error`` - The error associated with individual values measured during the experiment
* ``unknown_parameters`` - Those parameters in the model that are estimated using the measured values during the experiment

An example ``Experiment`` object that builds and labels the model is shown in the next few sections.

Pyomo.DoE Usage Example
-----------------------
Expand Down Expand Up @@ -203,89 +153,87 @@ The goal of MBDoE is to optimize the experiment design variables :math:`\boldsym
The observation errors are assumed to be independent both in time and across measurements with a constant standard deviation of 1 M for each species.


Step 0: Import Pyomo and the Pyomo.DoE module
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Step 0: Import Pyomo and the Pyomo.DoE module and create an ``Experiment`` class
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. doctest::

>>> # === Required import ===
>>> import pyomo.environ as pyo
>>> from pyomo.dae import ContinuousSet, DerivativeVar
>>> from pyomo.contrib.doe import DesignOfExperiments, MeasurementVariables, DesignVariables
>>> from pyomo.contrib.doe import DesignOfExperiments
>>> import numpy as np

.. literalinclude:: ../../../../pyomo/contrib/doe/examples/reactor_experiment.py
:start-after: ========================
:end-before: End constructor definition

Step 1: Define the Pyomo process model
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

The process model for the reaction kinetics problem is shown below.
The process model for the reaction kinetics problem is shown below. We build the model without any data or discretization.

.. literalinclude:: ../../../../pyomo/contrib/doe/examples/reactor_kinetics.py
:language: python
:pyobject: create_model
.. literalinclude:: ../../../../pyomo/contrib/doe/examples/reactor_experiment.py
:start-after: Create flexible model without data
:end-before: End equation definition

.. literalinclude:: ../../../../pyomo/contrib/doe/examples/reactor_kinetics.py
:language: python
:pyobject: disc_for_measure
Step 2: Finalize the Pyomo process model
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. note::
The model requires at least two options: "block" and "global". Both options requires the pass of a created empty Pyomo model.
With "global" option, only design variables and their time sets need to be defined;
With "block" option, a full model needs to be defined.
Here we add data to the model and finalize the discretization. This step is required before the model can be labeled.

.. literalinclude:: ../../../../pyomo/contrib/doe/examples/reactor_experiment.py
:start-after: End equation definition
:end-before: End model finalization

Step 2: Define the inputs for Pyomo.DoE
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. literalinclude:: ../../../../pyomo/contrib/doe/examples/reactor_compute_FIM.py
:language: python
:start-at: # Control time set
:end-before: ### Compute
Step 3: Label the information needed for DoE analysis
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

We label the four important groups as defined before.

Step 3: Compute the FIM of a square MBDoE problem
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. literalinclude:: ../../../../pyomo/contrib/doe/examples/reactor_experiment.py
:start-after: End model finalization
:end-before: End model labeling

This method computes an MBDoE optimization problem with no degree of freedom.
Step 4: Implement the ``get_labeled_model`` method
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

This method can be accomplished by two modes, ``direct_kaug`` and ``sequential_finite``.
``direct_kaug`` mode requires the installation of the solver `k_aug <https://github.com/dthierry/k_aug>`_.
This method utilizes the previous 3 steps and is used by `Pyomo.DoE` to build the model to perform optimal experimental design.

.. literalinclude:: ../../../../pyomo/contrib/doe/examples/reactor_compute_FIM.py
:language: python
:start-after: ### Compute the FIM
:end-before: # test result
.. literalinclude:: ../../../../pyomo/contrib/doe/examples/reactor_experiment.py
:start-after: End constructor definition
:end-before: Create flexible model without data

Step 4: Exploratory analysis (Enumeration)
Step 5: Exploratory analysis (Enumeration)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Exploratory analysis is suggested to enumerate the design space to check if the problem is identifiable,
i.e., ensure that D-, E-optimality metrics are not small numbers near zero, and Modified E-optimality is not a big number.

Pyomo.DoE accomplishes the exploratory analysis with the ``run_grid_search`` function.
It allows users to define any number of design decisions. Heatmaps can be drawn by two design variables, fixing other design variables.
1D curve can be drawn by one design variable, fixing all other variables.
The function ``run_grid_search`` enumerates over the design space, each MBDoE problem accomplished by ``compute_FIM`` method.
Therefore, ``run_grid_search`` supports only two modes: ``sequential_finite`` and ``direct_kaug``.
Pyomo.DoE can perform exploratory sensitivity analysis with the ``compute_FIM_full_factorial`` function.
The ``compute_FIM_full_factorial`` function generates a grid over the design space as specified by the user. Each grid point represents an MBDoE problem solved using ``compute_FIM`` method. In this way, sensitivity of the FIM over the design space can be evaluated.

.. literalinclude:: ../../../../pyomo/contrib/doe/examples/reactor_grid_search.py
:language: python
:pyobject: main
The following code executes the above problem description:

Successful run of the above code shows the following figure:
.. literalinclude:: ../../../../pyomo/contrib/doe/examples/reactor_example.py
:start-after: Read in file
:end-before: End sensitivity analysis

.. figure:: grid-1.png
:scale: 35 %
An example output of the code above, a design exploration for the initial concentration and temperature as experimental design variables with 9 values, produces the four figures summarized below:

A heatmap shows the change of the objective function, a.k.a. the experimental information content, in the design region. Horizontal and vertical axes are two design variables, while the color of each grid shows the experimental information content. Taking the Fig. Reactor case - A optimality as example, A-optimality shows that the most informative region is around $C_{A0}=5.0$ M, $T=300.0$ K, while the least informative region is around $C_{A0}=1.0$ M, $T=700.0$ K.
.. figure:: FIM_sensitivity.png
:scale: 50 %

Step 5: Gradient-based optimization
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
A heatmap shows the change of the objective function, a.k.a. the experimental information content, in the design space. Horizontal and vertical axes are the two experimental design variables, while the color of each grid shows the experimental information content. For A optimality (top left subfigure), the figure shows that the most informative region is around :math:`C_{A0}=5.0` M, :math:`T=300.0` K, while the least informative region is around :math:`C_{A0}=1.0` M, :math:`T=700.0` K.

Step 6: Performing an optimal experimental design
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Pyomo.DoE accomplishes gradient-based optimization with the ``stochastic_program`` function for A- and D-optimality design.
In step 5, the DoE object was constructed to perform an exploratory sensitivity analysis. The same object can be used to design an optimal experiment with a single line of code.

This function solves twice: It solves the square version of the MBDoE problem first, and then unfixes the design variables as degree of freedoms and solves again. In this way the optimization problem can be well initialized.
.. literalinclude:: ../../../../pyomo/contrib/doe/examples/reactor_example.py
:start-after: Begin optimal DoE
:end-before: Print out a results summary

.. literalinclude:: ../../../../pyomo/contrib/doe/examples/reactor_optimize_doe.py
:language: python
:pyobject: main
When run, the optimal design is an initial concentration of 5.0 mol/L and an initial temperature of 494 K with all other temperatures being 300 K. The corresponding log-10 determinant of the FIM is 13.75


6 changes: 3 additions & 3 deletions doc/OnlineDocs/contributed_packages/gdpopt.rst
Original file line number Diff line number Diff line change
Expand Up @@ -93,10 +93,10 @@ An example that includes the modeling approach may be found below.
Variables:
x : Size=1, Index=None
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : -1.2 : 0.0 : 2 : False : False : Reals
None : -1.2 : 0 : 2 : False : False : Reals
y : Size=1, Index=None
Key : Lower : Value : Upper : Fixed : Stale : Domain
None : -10 : 1.0 : 10 : False : False : Reals
None : -10 : 1 : 10 : False : False : Reals
<BLANKLINE>
Objectives:
objective : Size=1, Index=None, Active=True
Expand All @@ -106,7 +106,7 @@ An example that includes the modeling approach may be found below.
Constraints:
c : Size=1
Key : Lower : Body : Upper
None : 1.0 : 1.0 : 1.0
None : 1.0 : 1 : 1.0

.. note::

Expand Down
1 change: 1 addition & 0 deletions doc/OnlineDocs/contributed_packages/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ Contributed packages distributed with Pyomo:
.. toctree::
:maxdepth: 1

alternative_solutions.rst
community.rst
doe/doe.rst
gdpopt.rst
Expand Down
9 changes: 8 additions & 1 deletion pyomo/common/fileutils.py
Original file line number Diff line number Diff line change
Expand Up @@ -286,10 +286,17 @@ def find_dir(
)


_exeExt = {'linux': None, 'windows': '.exe', 'cygwin': '.exe', 'darwin': None}
_exeExt = {
'linux': None,
'freebsd': None,
'windows': '.exe',
'cygwin': '.exe',
'darwin': None,
}

_libExt = {
'linux': ('.so', '.so.*'),
'freebsd': ('.so', '.so.*'),
'windows': ('.dll', '.pyd'),
'cygwin': ('.dll', '.so', '.so.*'),
'darwin': ('.dylib', '.so', '.so.*'),
Expand Down
7 changes: 4 additions & 3 deletions pyomo/common/tests/test_download.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@
import pyomo.common.envvar as envvar

from pyomo.common import DeveloperError
from pyomo.common.fileutils import this_file
from pyomo.common.fileutils import this_file, Executable
from pyomo.common.download import FileDownloader, distro_available
from pyomo.common.log import LoggingIntercept
from pyomo.common.tee import capture_output
Expand Down Expand Up @@ -173,7 +173,8 @@ def test_get_os_version(self):
self.assertTrue(v.replace('.', '').startswith(dist_ver))

if (
subprocess.run(
Executable('lsb_release').available()
and subprocess.run(
['lsb_release'],
stdout=subprocess.DEVNULL,
stderr=subprocess.DEVNULL,
Expand Down Expand Up @@ -206,7 +207,7 @@ def test_get_os_version(self):
self.assertEqual(_os, 'win')
self.assertEqual(_norm, _os + ''.join(_ver.split('.')[:2]))
else:
self.assertEqual(ans, '')
self.assertEqual(_os, '')

self.assertEqual((_os, _ver), FileDownloader._os_version)
# Exercise the fetch from CACHE
Expand Down
Loading

0 comments on commit 2a6213d

Please sign in to comment.