diff --git a/doc/OnlineDocs/explanation/analysis/doe/FIM_sensitivity.png b/doc/OnlineDocs/explanation/analysis/doe/FIM_sensitivity.png deleted file mode 100644 index af6b75cbbea..00000000000 Binary files a/doc/OnlineDocs/explanation/analysis/doe/FIM_sensitivity.png and /dev/null differ diff --git a/doc/OnlineDocs/explanation/analysis/doe/doe.rst b/doc/OnlineDocs/explanation/analysis/doe/doe.rst index 6e77e55dfd6..ebf12f8ecfa 100644 --- a/doc/OnlineDocs/explanation/analysis/doe/doe.rst +++ b/doc/OnlineDocs/explanation/analysis/doe/doe.rst @@ -1,10 +1,19 @@ +.. _pyomo.doe: + Pyomo.DoE ========= - -**Pyomo.DoE** (Pyomo Design of Experiments) is a Python library for model-based design of experiments using science-based models. - -Pyomo.DoE was developed by **Jialu Wang** and **Alexander W. Dowling** at the University of Notre Dame as part of the `Carbon Capture Simulation for Industry Impact (CCSI2) `_. -project, funded through the U.S. Department Of Energy Office of Fossil Energy. +**Pyomo.DoE** (Pyomo Design of Experiments) is a Python library for model-based design +of experiments using science-based models. + +Pyomo.DoE was originally created by **Jialu Wang** and **Alexander W. Dowling** at the +University of Notre Dame as part of the `Carbon Capture Simulation for Industry Impact +(CCSI2) `_ +project, funded through the U.S. Department Of Energy Office of Fossil Energy with +assistance from **John Siirola**, **Bethany Nicholson**, **Miranda Mundt**, and **Hailey Lynch**. +Significant improvements and extensions were contributed by **Dan Laky**, and +**Shuvashish Mondal** with funding from the +`Process Optimization & Modeling for Minerals Sustainability (PrOMMiS) `_ project +and the `University of Notre Dame `_. If you use Pyomo.DoE, please cite: @@ -15,28 +24,40 @@ AIChE Journal 68.12 (2022): e17813. `https://doi.org/10.1002/aic.17813` Methodology Overview --------------------- -Model-based Design of Experiments (MBDoE) is a technique to maximize the information gain of experiments by directly using science-based models with physically meaningful parameters. It is one key component in the model calibration and uncertainty quantification workflow shown below: +Model-based Design of Experiments (MBDoE) is a technique to maximize the information +gain from experiments by directly using science-based models with physically meaningful +parameters. It is one key component in the model calibration and uncertainty +quantification workflow shown below: -.. figure:: flowchart.png - :scale: 25 % +.. figure:: pyomo_workflow_new.png + :align: center + :scale: 90 % - The exploratory analysis, parameter estimation, uncertainty analysis, and MBDoE are combined into an iterative framework to select, refine, and calibrate science-based mathematical models with quantified uncertainty. Currently, Pyomo.DoE focuses on increasing parameter precision. + The parameter estimation, uncertainty analysis, and MBDoE are + combined into an iterative framework to select, refine, and calibrate science-based + mathematical models with quantified uncertainty. Currently, Pyomo.DoE focuses on + increasing parameter precision. -Pyomo.DoE provides the exploratory analysis and MBDoE capabilities to the Pyomo ecosystem. The user provides one Pyomo model, a set of parameter nominal values, +Pyomo.DoE provides the exploratory analysis and MBDoE capabilities to the +Pyomo ecosystem. The user provides one Pyomo model, a set of parameter nominal values, the allowable design spaces for design variables, and the assumed observation error model. -During exploratory analysis, Pyomo.DoE checks if the model parameters can be inferred from the postulated measurements or preliminary data. +During exploratory analysis, Pyomo.DoE checks if the model parameters can be +inferred from the postulated measurements or preliminary data. MBDoE then recommends optimized experimental conditions for collecting more data. -Parameter estimation packages such as :ref:`Parmest ` can perform parameter estimation using the available data to infer values for parameters, +Parameter estimation packages such as :ref:`Parmest ` can perform +parameter estimation using the available data to infer values for parameters, and facilitate an uncertainty analysis to approximate the parameter covariance matrix. -If the parameter uncertainties are sufficiently small, the workflow terminates and returns the final model with quantified parametric uncertainty. -If not, MBDoE recommends optimized experimental conditions to generate new data. +If the parameter uncertainties are sufficiently small, the workflow terminates +and returns the final model with quantified parametric uncertainty. +If not, MBDoE recommends optimized experimental conditions to generate new data +that will maximize information gain and eventually reduce parameter uncertainty. Below is an overview of the type of optimization models Pyomo.DoE can accommodate: * Pyomo.DoE is suitable for optimization models of **continuous** variables * Pyomo.DoE can handle **equality constraints** defining state variables -* Pyomo.DoE supports (Partial) Differential-Algebraic Equations (PDAE) models via Pyomo.DAE -* Pyomo.DoE also supports models with only algebraic constraints +* Pyomo.DoE supports (Partial) Differential-Algebraic Equations (PDAE) models via :ref:`Pyomo.DAE ` +* Pyomo.DoE also supports models with only algebraic equations The general form of a DAE problem that can be passed into Pyomo.DoE is shown below: @@ -47,7 +68,7 @@ The general form of a DAE problem that can be passed into Pyomo.DoE is shown bel \dot{\mathbf{x}}(t) = \mathbf{f}(\mathbf{x}(t), \mathbf{z}(t), \mathbf{y}(t), \mathbf{u}(t), \overline{\mathbf{w}}, \boldsymbol{\theta}) \\ \mathbf{g}(\mathbf{x}(t), \mathbf{z}(t), \mathbf{y}(t), \mathbf{u}(t), \overline{\mathbf{w}},\boldsymbol{\theta})=\mathbf{0} \\ \mathbf{y} =\mathbf{h}(\mathbf{x}(t), \mathbf{z}(t), \mathbf{u}(t), \overline{\mathbf{w}},\boldsymbol{\theta}) \\ - \mathbf{f}^{\mathbf{0}}\left(\dot{\mathbf{x}}\left(t_{0}\right), \mathbf{x}\left(t_{0}\right), \mathbf{z}(t_0), \mathbf{y}(t_0), \mathbf{u}\left(t_{0}\right), \overline{\mathbf{w}}, \boldsymbol{\theta})\right)=\mathbf{0} \\ + \mathbf{f}^{\mathbf{0}}\left(\dot{\mathbf{x}}\left(t_{0}\right), \mathbf{x}\left(t_{0}\right), \mathbf{z}(t_0), \mathbf{y}(t_0), \mathbf{u}\left(t_{0}\right), \overline{\mathbf{w}}, \boldsymbol{\theta}\right)=\mathbf{0} \\ \mathbf{g}^{\mathbf{0}}\left( \mathbf{x}\left(t_{0}\right),\mathbf{z}(t_0), \mathbf{y}(t_0), \mathbf{u}\left(t_{0}\right), \overline{\mathbf{w}}, \boldsymbol{\theta}\right)=\mathbf{0}\\ \mathbf{y}^{\mathbf{0}}\left(t_{0}\right)=\mathbf{h}\left(\mathbf{x}\left(t_{0}\right),\mathbf{z}(t_0), \mathbf{u}\left(t_{0}\right), \overline{\mathbf{w}}, \boldsymbol{\theta}\right) \end{array}\] @@ -66,8 +87,9 @@ where: * :math:`\mathbf{t} \in \mathbb{R}^{N_t \times 1}` is a union of all time sets. .. note:: - * Parameters and design variables should be defined as Pyomo ``Var`` components on the model to use ``direct_kaug`` mode, and can be defined as Pyomo ``Param`` object if not using ``direct_kaug``. - + * Parameters and design variables should be defined as Pyomo ``Var`` components + when building the model using the ``Experiment`` class so that users can use both + ``Parmest`` and ``Pyomo.DoE`` seamlessly. Based on the above notation, the form of the MBDoE problem addressed in Pyomo.DoE is shown below: .. math:: @@ -75,26 +97,29 @@ Based on the above notation, the form of the MBDoE problem addressed in Pyomo.Do \begin{equation} \begin{aligned} - \underset{\boldsymbol{\varphi}}{\max} \quad & \Psi (\mathbf{M}(\mathbf{\hat{y}}, \boldsymbol{\varphi})) \\ - \text{s.t.} \quad & \mathbf{M}(\boldsymbol{\hat{\theta}}, \boldsymbol{\varphi}) = \sum_r^{N_r} \sum_{r'}^{N_r} \tilde{\sigma}_{(r,r')}\mathbf{Q}_r^\mathbf{T} \mathbf{Q}_{r'} + \mathbf{V}^{-1}_{\boldsymbol{\theta}}(\boldsymbol{\hat{\theta}}) \\ - & \dot{\mathbf{x}}(t) = \mathbf{f}(\mathbf{x}(t), \mathbf{z}(t), \mathbf{y}(t), \mathbf{u}(t), \overline{\mathbf{w}}, \boldsymbol{\theta}) \\ - & \mathbf{g}(\mathbf{x}(t), \mathbf{z}(t), \mathbf{y}(t), \mathbf{u}(t), \overline{\mathbf{w}},\boldsymbol{\theta})=\mathbf{0} \\ - & \mathbf{y} =\mathbf{h}(\mathbf{x}(t), \mathbf{z}(t), \mathbf{u}(t), \overline{\mathbf{w}},\boldsymbol{\theta}) \\ - & \mathbf{f}^{\mathbf{0}}\left(\dot{\mathbf{x}}\left(t_{0}\right), \mathbf{x}\left(t_{0}\right), \mathbf{z}(t_0), \mathbf{y}(t_0), \mathbf{u}\left(t_{0}\right), \overline{\mathbf{w}}, \boldsymbol{\theta})\right)=\mathbf{0} \\ - & \mathbf{g}^{\mathbf{0}}\left( \mathbf{x}\left(t_{0}\right),\mathbf{z}(t_0), \mathbf{y}(t_0), \mathbf{u}\left(t_{0}\right), \overline{\mathbf{w}}, \boldsymbol{\theta}\right)=\mathbf{0}\\ - &\mathbf{y}^{\mathbf{0}}\left(t_{0}\right)=\mathbf{h}\left(\mathbf{x}\left(t_{0}\right),\mathbf{z}(t_0), \mathbf{u}\left(t_{0}\right), \overline{\mathbf{w}}, \boldsymbol{\theta}\right) + \underset{\boldsymbol{\varphi}}{\max} \quad & \Psi (\mathbf{M}(\boldsymbol{\hat{\theta}}, \boldsymbol{\varphi})) \\ + \text{s.t.} \quad & \mathbf{M}(\boldsymbol{\hat{\theta}}, \boldsymbol{\varphi}) = \sum_r^{N_r} \sum_{r'}^{N_r} \tilde{\sigma}_{(r,r')}\mathbf{Q}_r^\mathbf{T} \mathbf{Q}_{r'} + \mathbf{M}_0 \\ + & \dot{\mathbf{x}}(t) = \mathbf{f}(\mathbf{x}(t), \mathbf{z}(t), \mathbf{y}(t), \mathbf{u}(t), \overline{\mathbf{w}}, \boldsymbol{\hat{\theta}}) \\ + & \mathbf{g}(\mathbf{x}(t), \mathbf{z}(t), \mathbf{y}(t), \mathbf{u}(t), \overline{\mathbf{w}},\boldsymbol{\hat{\theta}})=\mathbf{0} \\ + & \mathbf{y} =\mathbf{h}(\mathbf{x}(t), \mathbf{z}(t), \mathbf{u}(t), \overline{\mathbf{w}},\boldsymbol{\hat{\theta}}) \\ + & \mathbf{f}^{\mathbf{0}}\left(\dot{\mathbf{x}}\left(t_{0}\right), \mathbf{x}\left(t_{0}\right), \mathbf{z}(t_0), \mathbf{y}(t_0), \mathbf{u}\left(t_{0}\right), \overline{\mathbf{w}}, \boldsymbol{\hat{\theta}})\right)=\mathbf{0} \\ + & \mathbf{g}^{\mathbf{0}}\left( \mathbf{x}\left(t_{0}\right),\mathbf{z}(t_0), \mathbf{y}(t_0), \mathbf{u}\left(t_{0}\right), \overline{\mathbf{w}}, \boldsymbol{\hat{\theta}}\right)=\mathbf{0}\\ + &\mathbf{y}^{\mathbf{0}}\left(t_{0}\right)=\mathbf{h}\left(\mathbf{x}\left(t_{0}\right),\mathbf{z}(t_0), \mathbf{u}\left(t_{0}\right), \overline{\mathbf{w}}, \boldsymbol{\hat{\theta}}\right) \end{aligned} \end{equation} where: * :math:`\boldsymbol{\varphi}` are design variables, which are manipulated to maximize the information content of experiments. It should consist of one or more of :math:`\mathbf{u}(t), \mathbf{y}^{\mathbf{0}}({t_0}),\overline{\mathbf{w}}`. With a proper model formulation, the timepoints for control or measurements :math:`\mathbf{t}` can also be degrees of freedom. -* :math:`\mathbf{M}` is the Fisher information matrix (FIM), estimated as the inverse of the covariance matrix of parameter estimates :math:`\boldsymbol{\hat{\theta}}`. A large FIM indicates more information contained in the experiment for parameter estimation. +* :math:`\mathbf{M}` is the Fisher information matrix (FIM), approximated as the inverse of the covariance matrix of parameter estimates :math:`\boldsymbol{\hat{\theta}}`. A large FIM indicates more information contained in the experiment for parameter estimation. * :math:`\mathbf{Q}` is the dynamic sensitivity matrix, containing the partial derivatives of :math:`\mathbf{y}` with respect to :math:`\boldsymbol{\theta}`. -* :math:`\Psi` is the design criteria to measure FIM. -* :math:`\mathbf{V}_{\boldsymbol{\theta}}(\boldsymbol{\hat{\theta}})^{-1}` is the FIM of previous experiments. +* :math:`\Psi` is the scalar design criteria to measure the information content in the FIM. +* :math:`\mathbf{M}_0` is the sum of all the FIMs from previous experiments. -Pyomo.DoE provides four design criteria :math:`\Psi` to measure the size of FIM: +Pyomo.DoE provides five design criteria :math:`\Psi` to measure the information in the FIM. +The covariance matrix of parameter estimates is approximated as the inverse of the FIM, +i.e., :math:`\mathbf{V} \approx \mathbf{M}^{-1}`. +We can use the FIM or the covariance matrix to define the design criteria. .. list-table:: Pyomo.DoE design criteria :header-rows: 1 @@ -104,38 +129,57 @@ Pyomo.DoE provides four design criteria :math:`\Psi` to measure the size of FIM - Computation - Geometrical meaning * - A-optimality - - :math:`\text{trace}({\mathbf{M}})` - - Dimensions of the enclosing box of the confidence ellipse + - :math:`\text{trace}(\mathbf{V}) = \text{trace}(\mathbf{M}^{-1})` + - Minimizing this is equivalent to minimizing enclosing box of the confidence ellipse + * - Pseudo A-optimality + - :math:`\text{trace}(\mathbf{M})` + - Maximizing this is equivalent to maximizing the dimensions of the enclosing box of the Fisher Information Matrix * - D-optimality - - :math:`\text{det}({\mathbf{M}})` - - Volume of the confidence ellipse + - :math:`\det(\mathbf{M}) = 1/\det(\mathbf{V})` + - Maximizing this is equivalent to minimizing confidence-ellipsoid volume * - E-optimality - - :math:`\text{min eig}({\mathbf{M}})` - - Size of the longest axis of the confidence ellipse + - :math:`\lambda_{\min}(\mathbf{M}) = 1/\lambda_{\max}(\mathbf{V})` + - Maximizing this is equivalent to minimizing the longest axis of the confidence ellipse * - Modified E-optimality - - :math:`\text{cond}({\mathbf{M}})` - - Ratio of the longest axis to the shortest axis of the confidence ellipse + - :math:`\text{cond}(\mathbf{M}) = \text{cond}(\mathbf{V})` + - Minimizing this is equivalent to minimizing the ratio of the longest axis to the shortest axis of the confidence ellipse + +.. note:: + A confidence ellipse is a geometric representation of the uncertainty in parameter + estimates. It is derived from the covariance matrix :math:`\mathbf{V}`. In order to solve problems of the above, Pyomo.DoE implements the 2-stage stochastic program. Please see Wang and Dowling (2022) for details. Pyomo.DoE Required Inputs -------------------------------- -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, :ref:`Parmest `. The four ``Suffix`` components are: +The required input to the Pyomo.DoE is a subclass of the :ref:`Parmest ` ``Experiment`` class. +The subclass must have a ``get_labeled_model`` method which returns a Pyomo `ConcreteModel` +containing four Pyomo ``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, +:ref:`Parmest `. The four Pyomo ``Suffix`` components are: * ``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 +* ``measurement_error`` - The error associated with individual values measured during the experiment. It is passed as a standard deviation or square root of the diagonal elements of the observation error covariance matrix. Pyomo.DoE currently assumes that the observation errors are Gaussain and independent both in time and across measurements. * ``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. +An example of the subclassed ``Experiment`` object that builds and labels the model is shown in the next few sections. Pyomo.DoE Usage Example ----------------------- -We illustrate the use of Pyomo.DoE using a reaction kinetics example (Wang and Dowling, 2022). -The Arrhenius equations model the temperature dependence of the reaction rate coefficient :math:`k_1, k_2`. Assuming a first-order reaction mechanism gives the reaction rate model. Further, we assume only species A is fed to the reactor. +We illustrate the use of Pyomo.DoE using a reaction kinetics example (Wang and Dowling, 2022). + +.. math:: + :nowrap: + + \begin{equation} + A \xrightarrow{k_1} B \xrightarrow{k_2} C + \end{equation} +The Arrhenius equations model the temperature dependence of the reaction rate coefficient :math:`k_1, k_2`. Assuming a first-order reaction mechanism gives the reaction rate model. Further, we assume only species A is fed to the reactor. + .. math:: :nowrap: @@ -154,21 +198,33 @@ The Arrhenius equations model the temperature dependence of the reaction rate co :math:`C_A(t), C_B(t), C_C(t)` are the time-varying concentrations of the species A, B, C, respectively. -:math:`k_1, k_2` are the rates for the two chemical reactions using an Arrhenius equation with activation energies :math:`E_1, E_2` and pre-exponential factors :math:`A_1, A_2`. +:math:`k_1, k_2` are the rate constants for the two chemical reactions using an Arrhenius equation with activation energies :math:`E_1, E_2` and pre-exponential factors :math:`A_1, A_2`. The goal of MBDoE is to optimize the experiment design variables :math:`\boldsymbol{\varphi} = (C_{A0}, T(t))`, where :math:`C_{A0},T(t)` are the initial concentration of species A and the time-varying reactor temperature, to maximize the precision of unknown model parameters :math:`\boldsymbol{\theta} = (A_1, E_1, A_2, E_2)` by measuring :math:`\mathbf{y}(t)=(C_A(t), C_B(t), C_C(t))`. 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 and create an ``Experiment`` class ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +.. note:: + + This example uses the data file ``result.json``, located in the Pyomo repository at: + ``pyomo/contrib/doe/examples/result.json``, which contains the nominal parameter + values, and measurements for the reaction kinetics experiment. .. doctest:: - >>> # === Required import === + # === Required import === >>> import pyomo.environ as pyo + >>> from pyomo.dae import ContinuousSet, DerivativeVar >>> from pyomo.contrib.doe import DesignOfExperiments + >>> from pyomo.contrib.parmest.experiment import Experiment + >>> import pathlib >>> import numpy as np + >>> import json +Subclass the :ref:`Parmest ` ``Experiment`` class to define the reaction +kinetics experiment and build the Pyomo ConcreteModel. + .. literalinclude:: /../../pyomo/contrib/doe/examples/reactor_experiment.py :start-after: ======================== :end-before: End constructor definition @@ -176,7 +232,8 @@ Step 0: Import Pyomo and the Pyomo.DoE module and create an ``Experiment`` class Step 1: Define the Pyomo process model ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -The process model for the reaction kinetics problem is shown below. We build the model without any data or discretization. +The process model for the reaction kinetics problem is shown below. Here, we build +the model without any data or discretization. .. literalinclude:: /../../pyomo/contrib/doe/examples/reactor_experiment.py :start-after: Create flexible model without data @@ -185,7 +242,8 @@ The process model for the reaction kinetics problem is shown below. We build the Step 2: Finalize the Pyomo process model ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -Here we add data to the model and finalize the discretization. This step is required before the model can be labeled. +Here, we add data to the model and finalize the discretization using a new method to +the class. This step is required before the model can be labeled. .. literalinclude:: /../../pyomo/contrib/doe/examples/reactor_experiment.py :start-after: End equation definition @@ -194,7 +252,8 @@ Here we add data to the model and finalize the discretization. This step is requ Step 3: Label the information needed for DoE analysis ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -We label the four important groups as defined before. +We label the four important groups as Pyomo Suffix components as mentioned before by +adding a ``label_experiment`` method. .. literalinclude:: /../../pyomo/contrib/doe/examples/reactor_experiment.py :start-after: End model finalization @@ -203,7 +262,8 @@ We label the four important groups as defined before. Step 4: Implement the ``get_labeled_model`` method ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -This method utilizes the previous 3 steps and is used by `Pyomo.DoE` to build the model to perform optimal experimental design. +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_experiment.py :start-after: End constructor definition @@ -212,34 +272,117 @@ This method utilizes the previous 3 steps and is used by `Pyomo.DoE` to build th 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. +After creating the subclass of the ``Experiment`` class, 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. +Additionally, it helps to initialize the model for the optimal experimental design step. -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. +Pyomo.DoE can perform exploratory sensitivity analysis with the ``compute_FIM_full_factorial`` method. +The ``compute_FIM_full_factorial`` method 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. +Pyomo.DoE supports plotting the results from ``compute_FIM_full_factorial`` method +with the ``draw_factorial_figure`` method. -The following code executes the above problem description: +The following code defines the ``run_reactor_doe`` function. This function encapsulates +the workflow for both sensitivity analysis (Step 5) and optimal design (Step 6). .. literalinclude:: /../../pyomo/contrib/doe/examples/reactor_example.py - :start-after: Read in file - :end-before: End sensitivity analysis + :language: python + :start-after: # This software is distributed under the 3-clause BSD License. + :end-before: if __name__ == "__main__": + :linenos: -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: +After defining the function, we will call it to perform the exploratory analysis and +the optimal experimental design. -.. figure:: FIM_sensitivity.png - :scale: 50 % +.. literalinclude:: /../../pyomo/contrib/doe/examples/reactor_example.py + :language: python + :start-after: if __name__ == "__main__": + :dedent: 4 -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. +A design exploration for the initial concentration and temperature as experimental +design variables with 9 values for each, produces the the five figures for +five optimality criteria using ``compute_FIM_full_factorial`` and +``draw_factorial_figure`` methods as shown below: -Step 6: Performing an optimal experimental design -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +|plot1| |plot2| -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. +|plot3| |plot4| -.. literalinclude:: /../../pyomo/contrib/doe/examples/reactor_example.py - :start-after: Begin optimal DoE - :end-before: Print out a results summary +|plot5| + +.. |plot1| image:: example_reactor_compute_FIM_D_opt.png + :width: 48 % + +.. |plot2| image:: example_reactor_compute_FIM_A_opt.png + :width: 48 % + +.. |plot3| image:: example_reactor_compute_FIM_pseudo_A_opt.png + :width: 48 % + +.. |plot4| image:: example_reactor_compute_FIM_E_opt.png + :width: 48 % + +.. |plot5| image:: example_reactor_compute_FIM_ME_opt.png + :width: 48 % + +The heatmaps show the values of the objective functions, 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 example, the D-optimality (upper left subplot) heatmap figure shows that the +most informative region is around :math:`C_{A0}=5.0` M, :math:`T=500.0` K with +a :math:`\log_{10}` determinant of FIM being around 19, +while the least informative region is around :math:`C_{A0}=1.0` M, :math:`T=300.0` K, +with a :math:`\log_{10}` determinant of FIM being around -5. For D-, Pseudo A-, and +E-optimality we want to maximize the objective function, while for A- and Modified +E-optimality we want to minimize the objective function. + +However, in this sensitivity analysis plot (heatmap), we only varied the initial +concentration and the initial temperature, while the temperature at other time +points is fixed at 300 K. + +.. math:: + :nowrap: + + \[ + T(t) = \begin{cases} + T_0, & t \le 0.125 \\ + 300\ \text{K}, & t > 0.125 + \end{cases} + \] + +If :math:`T_0 = 300\ \text{K}`, the reaction is conducted under strictly isothermal +conditions. Because the temperature is constant, the sensitivities of the species +concentrations with respect to the Arrhenius parameters (:math:`A_i` and :math:`E_i`) +become linearly dependent. This high correlation means the effects of the +pre-exponential factor and the activation energy cannot be uniquely distinguished +from the measurements. Consequently, the Fisher Information Matrix (FIM) becomes +ill-conditioned, resulting in a near-zero determinant and a very large condition number. + +To break this correlation and make the parameters identifiable, introducing a time- +varying temperature profile (for example, a temperature step or a ramp) is required. +As shown in the heatmap, when the initial temperature :math:`T_0` differs from the +subsequent 300 K baseline, such a temperature change breaks the linear dependence, +yielding a well-conditioned FIM and identifiable parameters. -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 +Step 6: Performing an optimal experimental design +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +In Step 5, we defined the ``run_reactor_doe`` function. This function constructs +the DoE object and performs the exploratory sensitivity analysis. The way the function +is defined, it also proceeds immediately to the optimal experimental design step +(applying ``run_doe`` on the ``DesignOfExperiments`` object). +We can initialize the model with the result we obtained from the exploratory +analysis (optimal point from the heatmaps) to help the optimal design step to speed +up convergence. However, implementation of this initialization is not shown here. + +After applying ``run_doe`` on the ``DesignOfExperiments`` object, +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 :math:`\log_{10}` determinant of the FIM is 19.32. + diff --git a/doc/OnlineDocs/explanation/analysis/doe/example_reactor_compute_FIM_A_opt.png b/doc/OnlineDocs/explanation/analysis/doe/example_reactor_compute_FIM_A_opt.png new file mode 100644 index 00000000000..b642439bdaf Binary files /dev/null and b/doc/OnlineDocs/explanation/analysis/doe/example_reactor_compute_FIM_A_opt.png differ diff --git a/doc/OnlineDocs/explanation/analysis/doe/example_reactor_compute_FIM_D_opt.png b/doc/OnlineDocs/explanation/analysis/doe/example_reactor_compute_FIM_D_opt.png new file mode 100644 index 00000000000..33d9ccc00c4 Binary files /dev/null and b/doc/OnlineDocs/explanation/analysis/doe/example_reactor_compute_FIM_D_opt.png differ diff --git a/doc/OnlineDocs/explanation/analysis/doe/example_reactor_compute_FIM_E_opt.png b/doc/OnlineDocs/explanation/analysis/doe/example_reactor_compute_FIM_E_opt.png new file mode 100644 index 00000000000..fd6d2ea4d4f Binary files /dev/null and b/doc/OnlineDocs/explanation/analysis/doe/example_reactor_compute_FIM_E_opt.png differ diff --git a/doc/OnlineDocs/explanation/analysis/doe/example_reactor_compute_FIM_ME_opt.png b/doc/OnlineDocs/explanation/analysis/doe/example_reactor_compute_FIM_ME_opt.png new file mode 100644 index 00000000000..c5be372803d Binary files /dev/null and b/doc/OnlineDocs/explanation/analysis/doe/example_reactor_compute_FIM_ME_opt.png differ diff --git a/doc/OnlineDocs/explanation/analysis/doe/example_reactor_compute_FIM_pseudo_A_opt.png b/doc/OnlineDocs/explanation/analysis/doe/example_reactor_compute_FIM_pseudo_A_opt.png new file mode 100644 index 00000000000..3465244d52b Binary files /dev/null and b/doc/OnlineDocs/explanation/analysis/doe/example_reactor_compute_FIM_pseudo_A_opt.png differ diff --git a/doc/OnlineDocs/explanation/analysis/doe/flowchart.png b/doc/OnlineDocs/explanation/analysis/doe/flowchart.png deleted file mode 100644 index 2e66566d2f6..00000000000 Binary files a/doc/OnlineDocs/explanation/analysis/doe/flowchart.png and /dev/null differ diff --git a/doc/OnlineDocs/explanation/analysis/doe/grid-1.png b/doc/OnlineDocs/explanation/analysis/doe/grid-1.png deleted file mode 100644 index 6c96bbab7e0..00000000000 Binary files a/doc/OnlineDocs/explanation/analysis/doe/grid-1.png and /dev/null differ diff --git a/doc/OnlineDocs/explanation/analysis/doe/pyomo_workflow_new.png b/doc/OnlineDocs/explanation/analysis/doe/pyomo_workflow_new.png new file mode 100644 index 00000000000..ce27285afb0 Binary files /dev/null and b/doc/OnlineDocs/explanation/analysis/doe/pyomo_workflow_new.png differ diff --git a/doc/OnlineDocs/explanation/analysis/doe/uml.png b/doc/OnlineDocs/explanation/analysis/doe/uml.png deleted file mode 100644 index e5280987722..00000000000 Binary files a/doc/OnlineDocs/explanation/analysis/doe/uml.png and /dev/null differ diff --git a/doc/OnlineDocs/explanation/modeling/dae.rst b/doc/OnlineDocs/explanation/modeling/dae.rst index ac43ed35005..130813d61ff 100644 --- a/doc/OnlineDocs/explanation/modeling/dae.rst +++ b/doc/OnlineDocs/explanation/modeling/dae.rst @@ -1,3 +1,5 @@ +.. _pyomo.dae: + Dynamic Optimization with pyomo.DAE =================================== diff --git a/pyomo/contrib/doe/examples/reactor_example.py b/pyomo/contrib/doe/examples/reactor_example.py index c40799183bf..c480d954686 100644 --- a/pyomo/contrib/doe/examples/reactor_example.py +++ b/pyomo/contrib/doe/examples/reactor_example.py @@ -19,7 +19,10 @@ # Example for sensitivity analysis on the reactor experiment # After sensitivity analysis is done, we perform optimal DoE def run_reactor_doe( - n_points_for_design=9, + n_points_for_C0: int = 9, + n_points_for_T0: int = 9, + C0_bounds_for_factorial: tuple = (1, 5), + T0_bounds_for_factorial: tuple = (300, 700), compute_FIM_full_factorial=True, plot_factorial_results=True, figure_file_name="example_reactor_compute_FIM", @@ -31,8 +34,10 @@ def run_reactor_doe( Parameters ---------- - n_points_for_design : int, optional - number of points to use for the design ranges, by default 9 + n_points_for_C0 : int, optional + number of points to use for the C0 design range, by default 9 + n_points_for_T0 : int, optional + number of points to use for the T0 design range, by default 9 compute_FIM_full_factorial : bool, optional whether to compute the full factorial design, by default True plot_factorial_results : bool, optional @@ -90,8 +95,8 @@ def run_reactor_doe( if compute_FIM_full_factorial: # Make design ranges to compute the full factorial design design_ranges = { - "CA[0]": [1, 5, n_points_for_design], - "T[0]": [300, 700, n_points_for_design], + "CA[0]": [*C0_bounds_for_factorial, n_points_for_C0], + "T[0]": [*T0_bounds_for_factorial, n_points_for_T0], } # Compute the full factorial design with the sequential FIM calculation @@ -155,4 +160,4 @@ def run_reactor_doe( if __name__ == "__main__": - run_reactor_doe() + run_reactor_doe(figure_file_name=None) diff --git a/pyomo/contrib/doe/tests/test_doe_solve.py b/pyomo/contrib/doe/tests/test_doe_solve.py index 80a77b5f012..9e31accf894 100644 --- a/pyomo/contrib/doe/tests/test_doe_solve.py +++ b/pyomo/contrib/doe/tests/test_doe_solve.py @@ -616,7 +616,8 @@ def test_doe_full_factorial(self): 923636.598578955, ] ff = run_reactor_doe( - n_points_for_design=2, + n_points_for_C0=2, + n_points_for_T0=2, compute_FIM_full_factorial=False, plot_factorial_results=False, run_optimal_doe=False, @@ -884,7 +885,8 @@ def cleanup_files(): # Run the reactor example run_reactor_doe( - n_points_for_design=1, + n_points_for_C0=1, + n_points_for_T0=1, compute_FIM_full_factorial=True, plot_factorial_results=True, figure_file_name=prefix_linear, @@ -902,7 +904,8 @@ def cleanup_files(): # Run the reactor example with log scale run_reactor_doe( - n_points_for_design=1, + n_points_for_C0=1, + n_points_for_T0=1, compute_FIM_full_factorial=True, plot_factorial_results=True, figure_file_name=prefix_log,