diff --git a/doc/OnlineDocs/explanation/analysis/parmest/covariance.rst b/doc/OnlineDocs/explanation/analysis/parmest/covariance.rst index fa6c001d868..43640de6717 100644 --- a/doc/OnlineDocs/explanation/analysis/parmest/covariance.rst +++ b/doc/OnlineDocs/explanation/analysis/parmest/covariance.rst @@ -1,86 +1,106 @@ +.. _covariancesection: + Covariance Matrix Estimation ============================ -The uncertainty in the estimated parameters is quantified using the covariance matrix. -The diagonal of the covariance matrix contains the variance of the estimated parameters. -Assuming Gaussian independent and identically distributed measurement errors, the -covariance matrix of the estimated parameters can be computed using the following -methods which have been implemented in parmest. +The goal is parameter estimation (see :ref:`driversection` Section) is to estimate unknown model parameters +from experimental data. When the model parameters are estimated from the data, their accuracy is measured by +computing the covariance matrix. The diagonal of this covariance matrix contains the variance of the +estimated parameters which is used to calculate their uncertainty. Assuming Gaussian independent and identically +distributed measurement errors, the covariance matrix of the estimated parameters can be computed using the +following methods which have been implemented in parmest. 1. Reduced Hessian Method - When the objective function is the sum of squared errors (SSE): - :math:`\text{SSE} = \sum_{i = 1}^n \left(y_{i} - \hat{y}_{i}\right)^2`, - the covariance matrix is: + When the objective function is the sum of squared errors (SSE) for homogeneous data, defined as + :math:`\text{SSE} = \sum_{i = 1}^{n} \left(\boldsymbol{y}_{i} - \boldsymbol{f}(\boldsymbol{x}_{i}; + \boldsymbol{\theta})\right)^\text{T} \left(\boldsymbol{y}_{i} - \boldsymbol{f}(\boldsymbol{x}_{i}; + \boldsymbol{\theta})\right)`, the covariance matrix is: .. math:: - V_{\boldsymbol{\theta}} = 2 \sigma^2 \left(\frac{\partial^2 \text{SSE}} - {\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}}\right)^{-1}_{\boldsymbol{\theta} - = \boldsymbol{\theta}^*} - - When the objective function is the weighted SSE (WSSE): - :math:`\text{WSSE} = \frac{1}{2} \left(\mathbf{y} - f(\mathbf{x};\boldsymbol{\theta})\right)^\text{T} - \mathbf{W} \left(\mathbf{y} - f(\mathbf{x};\boldsymbol{\theta})\right)`, + \boldsymbol{V}_{\boldsymbol{\theta}} = 2 \sigma^2 \left(\frac{\partial^2 \text{SSE}} + {\partial \boldsymbol{\theta}^2}\right)^{-1}_{\boldsymbol{\theta} + = \hat{\boldsymbol{\theta}}} + + Similarly, when the objective function is the weighted SSE (WSSE) for heterogeneous data, defined as + :math:`\text{WSSE} = \frac{1}{2} \sum_{i = 1}^{n} \left(\boldsymbol{y}_{i} - + \boldsymbol{f}(\boldsymbol{x}_{i};\boldsymbol{\theta})\right)^\text{T} \boldsymbol{\Sigma}_{\boldsymbol{y}}^{-1} + \left(\boldsymbol{y}_{i} - \boldsymbol{f}(\boldsymbol{x}_{i};\boldsymbol{\theta})\right)`, the covariance matrix is: .. math:: - V_{\boldsymbol{\theta}} = \left(\frac{\partial^2 \text{WSSE}} - {\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}}\right)^{-1}_{\boldsymbol{\theta} - = \boldsymbol{\theta}^*} - - Where :math:`V_{\boldsymbol{\theta}}` is the covariance matrix of the estimated - parameters, :math:`y` are the observed measured variables, :math:`\hat{y}` are the - predicted measured variables, :math:`n` is the number of data points, - :math:`\boldsymbol{\theta}` are the unknown parameters, :math:`\boldsymbol{\theta^*}` - are the estimates of the unknown parameters, :math:`\mathbf{x}` are the decision - variables, and :math:`\mathbf{W}` is a diagonal matrix containing the inverse of the - variance of the measurement error, :math:`\sigma^2`. When the standard - deviation of the measurement error is not supplied by the user, parmest - approximates the variance of the measurement error as - :math:`\sigma^2 = \frac{1}{n-l} \sum e_i^2` where :math:`l` is the number of - fitted parameters, and :math:`e_i` is the residual for experiment :math:`i`. + \boldsymbol{V}_{\boldsymbol{\theta}} = \left(\frac{\partial^2 \text{WSSE}} + {\partial \boldsymbol{\theta}^2}\right)^{-1}_{\boldsymbol{\theta} + = \hat{\boldsymbol{\theta}}} + + Where :math:`\boldsymbol{V}_{\boldsymbol{\theta}}` is the covariance matrix of the estimated + parameters :math:`\hat{\boldsymbol{\theta}} \in \mathbb{R}^p`, :math:`\boldsymbol{y}_{i} \in \mathbb{R}^m` are + observations of the measured variables, :math:`\boldsymbol{f}` is the model function, + :math:`\boldsymbol{x}_{i} \in \mathbb{R}^{q}` are the input variables, :math:`n` is the number of experiments, + :math:`\boldsymbol{\Sigma}_{\boldsymbol{y}}` is the measurement error covariance matrix, and :math:`\sigma^2` + is the variance of the measurement error. When the standard deviation of the measurement error is not supplied + by the user, parmest approximates :math:`\sigma^2` as: + :math:`\hat{\sigma}^2 = \frac{1}{n-p} \sum_{i=1}^{n} \boldsymbol{\varepsilon}_{i}(\boldsymbol{\theta})^{\text{T}} + \boldsymbol{\varepsilon}_{i}(\boldsymbol{\theta})`, and :math:`\boldsymbol{\varepsilon}_{i} \in \mathbb{R}^m` + are the residuals between the data and model for experiment :math:`i`. In parmest, this method computes the inverse of the Hessian by scaling the - objective function (SSE or WSSE) with a constant probability factor. + objective function (SSE or WSSE) with a constant probability factor, :math:`\frac{1}{n}`. 2. Finite Difference Method - In this method, the covariance matrix, :math:`V_{\boldsymbol{\theta}}`, is - calculated by applying the Gauss-Newton approximation to the Hessian, - :math:`\frac{\partial^2 \text{SSE}}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}}` + In this method, the covariance matrix, :math:`\boldsymbol{V}_{\boldsymbol{\theta}}`, is + computed by differentiating the Hessian, + :math:`\frac{\partial^2 \text{SSE}}{\partial \boldsymbol{\theta}^2}` or - :math:`\frac{\partial^2 \text{WSSE}}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}}`, - leading to: + :math:`\frac{\partial^2 \text{WSSE}}{\partial \boldsymbol{\theta}^2}`, and + applying Gauss-Newton approximation which results in: + + .. math:: + \boldsymbol{V}_{\boldsymbol{\theta}} = \left(\sum_{i = 1}^n \boldsymbol{G}_{i}^{\text{T}} + \boldsymbol{\Sigma}_{\boldsymbol{y}}^{-1} \boldsymbol{G}_{i} \right)^{-1} + + where .. math:: - V_{\boldsymbol{\theta}} = \left(\sum_{i = 1}^n \mathbf{G}_{i}^{\mathrm{T}} \mathbf{W} - \mathbf{G}_{i} \right)^{-1} + \boldsymbol{G}_{i} = \frac{\partial \boldsymbol{f}(\boldsymbol{x}_i;\boldsymbol{\theta})} + {\partial \boldsymbol{\theta}} - This method uses central finite difference to compute the Jacobian matrix, - :math:`\mathbf{G}_{i}`, for experiment :math:`i`, which is the sensitivity of - the measured variables with respect to the parameters, :math:`\boldsymbol{\theta}`. - :math:`\mathbf{W}` is a diagonal matrix containing the inverse of the variance - of the measurement errors, :math:`\sigma^2`. + This method uses central finite difference to compute the Jacobian matrix, :math:`\boldsymbol{G}_{i}`, + for experiment :math:`i`. + + .. math:: + \boldsymbol{G}_{i}[:,\,k] \approx \frac{\boldsymbol{f}(\boldsymbol{x}_i;\theta_k + \Delta \theta_k) + \vert_{\hat{\boldsymbol{\theta}}} - \boldsymbol{f}(\boldsymbol{x}_i;\theta_k - \Delta \theta_k) + \vert_{\hat{\boldsymbol{\theta}}}}{2 \Delta \theta_k} \quad \forall \quad \theta_k \, \in \, + [\theta_1,\cdots, \theta_p] 3. Automatic Differentiation Method Similar to the finite difference method, the covariance matrix is calculated as: .. math:: - V_{\boldsymbol{\theta}} = \left( \sum_{i = 1}^n \mathbf{G}_{\text{kaug},\, i}^{\mathrm{T}} - \mathbf{W} \mathbf{G}_{\text{kaug},\, i} \right)^{-1} + \boldsymbol{V}_{\boldsymbol{\theta}} = \left( \sum_{i = 1}^n \boldsymbol{G}_{\text{kaug},\, i}^{\text{T}} + \boldsymbol{\Sigma}_{\boldsymbol{y}}^{-1} \boldsymbol{G}_{\text{kaug},\, i} \right)^{-1} + + However, this method uses implicit differentiation and the model-optimality or Karush–Kuhn–Tucker (KKT) conditions + to compute the Jacobian matrix, :math:`\boldsymbol{G}_{\text{kaug},\, i}`, for experiment :math:`i`. + + .. math:: + \boldsymbol{G}_{\text{kaug},\,i} = \frac{\partial \boldsymbol{f}(\boldsymbol{x}_i,\boldsymbol{\theta})} + {\partial \boldsymbol{\theta}} + \frac{\partial \boldsymbol{f}(\boldsymbol{x}_i,\boldsymbol{\theta})} + {\partial \boldsymbol{x}_i}\frac{\partial \boldsymbol{x}_i}{\partial \boldsymbol{\theta}} - However, this method uses the model optimality (KKT) condition to compute the - Jacobian matrix, :math:`\mathbf{G}_{\text{kaug},\, i}`, for experiment :math:`i`. +The covariance matrix calculation is only supported with the built-in objective functions "SSE" or "SSE_weighted". -The covariance matrix calculation is only supported with the built-in objective -functions "SSE" or "SSE_weighted". +In parmest, the covariance matrix can be computed after creating the +:class:`~pyomo.contrib.parmest.experiment.Experiment` class, +defining the :class:`~pyomo.contrib.parmest.parmest.Estimator` object, +and estimating the model parameters using :class:`~pyomo.contrib.parmest.parmest.Estimator.theta_est` +(all these steps were addressed in the :ref:`driversection` Section). -In parmest, the covariance matrix can be calculated after defining the -:class:`~pyomo.contrib.parmest.parmest.Estimator` object and estimating the unknown -parameters using :class:`~pyomo.contrib.parmest.parmest.Estimator.theta_est`. To -estimate the covariance matrix, with the default method being "finite_difference", call -the :class:`~pyomo.contrib.parmest.parmest.Estimator.cov_est` function, e.g., +To estimate the covariance matrix, with the default method being "finite_difference", call +the :class:`~pyomo.contrib.parmest.parmest.Estimator.cov_est` function as follows: .. testsetup:: * :skipif: not __import__('pyomo.contrib.parmest.parmest').contrib.parmest.parmest.parmest_available diff --git a/doc/OnlineDocs/explanation/analysis/parmest/driver.rst b/doc/OnlineDocs/explanation/analysis/parmest/driver.rst index 9d4d2cca1c3..50cbbb25bb7 100644 --- a/doc/OnlineDocs/explanation/analysis/parmest/driver.rst +++ b/doc/OnlineDocs/explanation/analysis/parmest/driver.rst @@ -1,7 +1,7 @@ .. _driversection: Parameter Estimation -================================== +==================== Parameter Estimation using parmest requires a Pyomo model, experimental data which defines multiple scenarios, and parameters @@ -23,8 +23,170 @@ modified within parmest to match the required specifications. The stochastic programming callback function is also defined within parmest. The callback function returns a populated and initialized model for each scenario. -To use parmest, the user creates a :class:`~pyomo.contrib.parmest.parmest.Estimator` -object which includes the following methods: +Quick Guide +=========== +We use a simple model to provide a quick guide on how to use parmest to estimate model parameters +from experimental data as well as compute their uncertainty. The model and data used in this +guide is taken from: Rooney, W. C.; Biegler, L. T. Design for Model Parameter Uncertainty Using +Nonlinear Confidence Regions. AIChE J. 2001, 47 (8), 1794–1804. + +The mathematical model of interest is: + +.. math:: + y_i(\theta_1, \theta_2, t_i) = \theta_1 \left(1 - e^{-\theta_2 t_i} \right), + \quad \forall \quad i \, \in \, {1, \ldots, n} + +Where :math:`y` is the observation of the measured variable, :math:`t` is the time, :math:`\theta_1` +is the asymptote, and :math:`\theta_2` is the rate constant. + +The experimental data is given in the table below: + +.. list-table:: Data + :header-rows: 1 + :widths: 30 30 + + * - hour + - y + * - 1 + - 8.3 + * - 2 + - 10.3 + * - 3 + - 19.0 + * - 4 + - 16.0 + * - 5 + - 15.6 + * - 7 + - 19.8 + +To use parmest to estimate :math:`\theta_1` and :math:`\theta_2` from the data, we provide the following +detailed steps: + +Step 0: Import Pyomo, parmest, Experiment Class, and Pandas +----------------------------------------------------------- + +Before solving the parameter estimation problem, the following code must be executed to import the +required packages for parameter estimation in parmest: + +.. doctest:: + :skipif: not __import__('pyomo.contrib.parmest.parmest').contrib.parmest.parmest.parmest_available + + >>> import pyomo.environ as pyo + >>> import pyomo.contrib.parmest.parmest as parmest + >>> from pyomo.contrib.parmest.experiment import Experiment + >>> import pandas as pd + +.. _ExperimentClass: + +Step 1: Create the Experiment Class for the Model +------------------------------------------------- + +parmest requires that the user creates an :class:`~pyomo.contrib.parmest.experiment.Experiment` class that +contains the annotated Pyomo model with labeled experiment outputs, unknown parameters, and measurement errors. + +A labeled Pyomo model ``m`` has the following additional suffixes (Pyomo `Suffix`): + +* ``m.experiment_outputs`` which defines experiment output (Pyomo `Param`, `Var`, or `Expression`) + and their associated data values (float, int). +* ``m.unknown_parameters`` which defines the mutable parameters or variables (Pyomo `Param` or `Var`) + to estimate along with their component unique identifier (Pyomo `ComponentUID`). + Within parmest, any parameters that are to be estimated are converted to unfixed variables. + Variables that are to be estimated are also unfixed. + +The experiment class has one required method: + +* :class:`~pyomo.contrib.parmest.experiment.Experiment.get_labeled_model` which returns the labeled Pyomo model. + Note that the model does not have to be specifically written as a + two-stage stochastic programming problem for parmest. + That is, parmest can modify the + objective, see the :ref:`EstimatorObj` Section below. + +This step shows how to create the :class:`~pyomo.contrib.parmest.experiment.Experiment` class using the +mathematical model outlined in the introduction section of this Quick Start. + +.. doctest:: + + >>> class RooneyBieglerExperiment(Experiment): + ... def __init__(self, data): + ... self.data = data + ... self.model = None + ... + ... def create_model(self): + ... # the model expects a dataframe + ... data_df = self.data.to_frame().transpose() + ... + ... # create the pyomo model + ... m = self.model = pyo.ConcreteModel() + ... + ... # add asymptote and rate constant to the model + ... m.asymptote = pyo.Var(initialize=15) + ... m.rate_constant = pyo.Var(initialize=0.5) + ... + ... # add the measured variable, y, to the model + ... m.y = pyo.Var(data_df.hour, within=pyo.PositiveReals, initialize=5) + ... + ... # add the mathematical equation for predicting y + ... def response_rule(m, h): + ... return m.y[h] == m.asymptote * (1 - pyo.exp(-m.rate_constant * h)) + ... + ... m.response_function = pyo.Constraint(data_df.hour, rule=response_rule) + ... + ... return m + ... + ... def label_model(self): + ... + ... m = self.model + ... + ... # label the experiment outputs + ... m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL) + ... m.experiment_outputs.update([(m.y[self.data['hour']], self.data['y'])]) + ... + ... # label the unknown parameters in the model + ... m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL) + ... m.unknown_parameters.update( + ... (k, pyo.ComponentUID(k)) for k in [m.asymptote, m.rate_constant] + ... ) + ... + ... # add the measurement error assumed to be constant and 0.1 + ... m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL) + ... m.measurement_error.update([(m.y[self.data['hour']], 0.1)]) + ... + ... def finalize_model(self): + ... + ... m = self.model + ... pass + ... + ... def get_labeled_model(self): + ... self.create_model() + ... self.label_model() + ... self.finalize_model() + ... + ... return self.model + + +Step 2: Load the Data and Create a List of the Model's Experiment Class +----------------------------------------------------------------------- + +After creating an :class:`~pyomo.contrib.parmest.experiment.Experiment` class instance for the model, a list of the +model's :class:`~pyomo.contrib.parmest.experiment.Experiment` class for all the experimental data points should be +created. + +.. doctest:: + + >>> data = pd.DataFrame(data=[[1, 8.3], [2, 10.3], [3, 19.0], [4, 16.0], [5, 15.6], [7, 19.8]], + ... columns=["hour", "y"]) + >>> exp_list = [] + >>> for i in range(data.shape[0]): + ... exp_list.append(RooneyBieglerExperiment(data.loc[i, :])) + +.. _EstimatorObj: + +Step 3: Create the Estimator Object +----------------------------------- + +To use parmest, the user creates an :class:`~pyomo.contrib.parmest.parmest.Estimator` object which includes +the following methods: .. autosummary:: :nosignatures: @@ -52,33 +214,20 @@ results and fit distributions to theta values. ~pyomo.contrib.parmest.graphics.fit_kde_dist A :class:`~pyomo.contrib.parmest.parmest.Estimator` object can be -created using the following code. A description of each argument is -listed below. Examples are provided in the :ref:`examplesection` -Section. +created using the following code. A description of the arguments are +listed below. -.. testsetup:: * +.. doctest:: :skipif: not __import__('pyomo.contrib.parmest.parmest').contrib.parmest.parmest.parmest_available - # Data - import pandas as pd - data = pd.DataFrame( - data=[[1, 8.3], [2, 10.3], [3, 19.0], - [4, 16.0], [5, 15.6], [7, 19.8]], - columns=['hour', 'y'], - ) - - # Create an experiment list - from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import RooneyBieglerExperiment + >>> pest = parmest.Estimator(exp_list, obj_function="SSE") - exp_list = [] - for i in range(data.shape[0]): - exp_list.append(RooneyBieglerExperiment(data.loc[i, :])) +Alternatively, the weighted sum of squared errors objective can be used. .. doctest:: :skipif: not __import__('pyomo.contrib.parmest.parmest').contrib.parmest.parmest.parmest_available - >>> import pyomo.contrib.parmest.parmest as parmest - >>> pest = parmest.Estimator(exp_list, obj_function="SSE") + >>> pest = parmest.Estimator(exp_list, obj_function="SSE_weighted") Optionally, solver options can be supplied, e.g., @@ -88,66 +237,47 @@ Optionally, solver options can be supplied, e.g., >>> solver_options = {"max_iter": 6000} >>> pest = parmest.Estimator(exp_list, obj_function="SSE", solver_options=solver_options) - -List of experiment objects --------------------------- - -The first argument is a list of experiment objects which is used to -create one labeled model for each experiment. -The template :class:`~pyomo.contrib.parmest.experiment.Experiment` -can be used to generate a list of experiment objects. - -A labeled Pyomo model ``m`` has the following additional suffixes (Pyomo `Suffix`): - -* ``m.experiment_outputs`` which defines experiment output (Pyomo `Param`, `Var`, or `Expression`) - and their associated data values (float, int). -* ``m.unknown_parameters`` which defines the mutable parameters or variables (Pyomo `Param` or `Var`) - to estimate along with their component unique identifier (Pyomo `ComponentUID`). - Within parmest, any parameters that are to be estimated are converted to unfixed variables. - Variables that are to be estimated are also unfixed. - -The experiment class has one required method: - -* :class:`~pyomo.contrib.parmest.experiment.Experiment.get_labeled_model` which returns the labeled Pyomo model. - Note that the model does not have to be specifically written as a - two-stage stochastic programming problem for parmest. - That is, parmest can modify the - objective, see :ref:`ObjFunction` below. - -Parmest comes with several :ref:`examplesection` that illustrates how to set up the list of experiment objects. -The examples commonly include additional :class:`~pyomo.contrib.parmest.experiment.Experiment` class methods to -create the model, finalize the model, and label the model. The user can customize methods to suit their needs. - -.. _ObjFunction: - Objective function ------------------- +^^^^^^^^^^^^^^^^^^ -The second argument is an optional argument which defines the -optimization objective function to use in parameter estimation. +The second argument is an optional argument if the objective function has already been included in the +Pyomo model, which defines the optimization objective function to use in parameter estimation. However, if the +objective function has not been included in the Pyomo model, like the one in the :ref:`ExperimentClass` Section +above, the user is required to supply the second argument. If no objective function is specified, the Pyomo model is used "as is" and should be defined with "FirstStageCost" and "SecondStageCost" expressions that are used to build an objective for the two-stage -stochastic programming problem. +stochastic programming problem. If the Pyomo model is not written as a two-stage stochastic programming problem in this format, the user can select the "SSE" or "SSE_weighted" built-in objective functions. If the user wants to use an objective that is different from the built-in options, a custom objective function can be defined for parameter estimation. However, -covariance matrix estimation will not support this custom objective function. The objective -function (built-in or custom) has a single argument, which is the model from a single -experiment. -The objective function returns a Pyomo -expression which is used to define "SecondStageCost". The objective -function can be used to customize data points and weights that are used -in parameter estimation. +covariance matrix estimation (see :ref:`covariancesection` Section) will not support this +custom objective function. The objective function (built-in or custom) has a single argument, +which is the model from a single experiment. The objective function returns a Pyomo expression +which is used to define "SecondStageCost". The objective function can be used to customize data +points and weights that are used in parameter estimation. Parmest includes two built-in objective functions ("SSE" and "SSE_weighted") to compute the sum of squared errors between the ``m.experiment_outputs`` model values and data values. -Suggested initialization procedure for parameter estimation problems +Step 3: Estimate the Parameters +------------------------------- + +After creating the :class:`~pyomo.contrib.parmest.parmest.Estimator` object with the desired objective function, +solve the parameter estimation problem by calling :class:`~pyomo.contrib.parmest.parmest.Estimator.theta_est`, +e.g., + +.. doctest:: + :skipif: not __import__('pyomo.contrib.parmest.parmest').contrib.parmest.parmest.parmest_available + + >>> pest = parmest.Estimator(exp_list, obj_function="SSE") + >>> obj_val, theta_val = pest.theta_est() + +Suggested Initialization Procedure for Parameter Estimation Problems -------------------------------------------------------------------- To check the quality of initial guess values provided for the fitted parameters, we suggest solving a @@ -160,11 +290,8 @@ estimation solve from the square problem solution, set optional argument ``solve argument ``(initialize_parmest_model=True)``. Different initial guess values for the fitted parameters can be provided using optional argument `theta_values` (**Pandas Dataframe**) -3. Solve parameter estimation problem by calling -:class:`~pyomo.contrib.parmest.parmest.Estimator.theta_est`, e.g., +More Examples Beyond this Quick Guide +------------------------------------- -.. doctest:: - :skipif: not parmest_available - - >>> pest = parmest.Estimator(exp_list, obj_function="SSE") - >>> obj_val, theta_val = pest.theta_est() +More detailed examples, such as parameter estimation of reaction kinetics are provided in the +:ref:`examplesection` Section. diff --git a/doc/OnlineDocs/explanation/analysis/parmest/installation.rst b/doc/OnlineDocs/explanation/analysis/parmest/installation.rst index 5cc41878719..c08527483a8 100644 --- a/doc/OnlineDocs/explanation/analysis/parmest/installation.rst +++ b/doc/OnlineDocs/explanation/analysis/parmest/installation.rst @@ -6,17 +6,30 @@ To run parmest, you will need Python version 3.x along with various Python package dependencies and the IPOPT software library for non-linear optimization. -Python package dependencies +Python Package Dependencies --------------------------- -#. numpy -#. pandas -#. pyomo -#. mpisppy (optional) -#. matplotlib (optional) -#. scipy.stats (optional) -#. seaborn (optional) -#. mpi4py.MPI (optional) +1. Install NumPy and Pandas with your preferred package manager; + both NumPy and SciPy are required dependencies of parmest. + You may install NumPy and Pandas with, for example, ``conda``: + + :: + + conda install numpy pandas + + or ``pip``: + + :: + + pip install numpy pandas +2. `Install Pyomo `_. + parmest is included in the Pyomo software package, at pyomo/contrib/parmest. +3. (*Optional*) Install ``matplotlib`` and ``scipy``: + + :: + + pip install scipy matplotlib +4. (*Optional*) Install ``mpisppy``, ``mpi4py.MPI``, and ``seaborn`` IPOPT ----- diff --git a/doc/OnlineDocs/explanation/analysis/parmest/overview.rst b/doc/OnlineDocs/explanation/analysis/parmest/overview.rst index f45c9c785b3..7d1e95c787d 100644 --- a/doc/OnlineDocs/explanation/analysis/parmest/overview.rst +++ b/doc/OnlineDocs/explanation/analysis/parmest/overview.rst @@ -10,7 +10,8 @@ for design optimization. Functionality in parmest includes: -* Model based parameter estimation using experimental data +* Model-based parameter estimation using experimental data +* Covariance matrix estimation * Bootstrap resampling for parameter estimation * Confidence regions based on single or multi-variate distributions * Likelihood ratio @@ -21,61 +22,56 @@ Background ---------- The goal of parameter estimation is to estimate values for -a vector, :math:`{\theta}`, to use in the functional form +a vector, :math:`\boldsymbol{\theta}`, to use in the functional form .. math:: - y = g(x; \theta) - -where :math:`x` is a vector containing measured data, typically in high -dimension, :math:`{\theta}` is a vector of values to estimate, in much -lower dimension, and the response vectors are given as :math:`y_{i}, -i=1,\ldots,m` with :math:`m` also much smaller than the dimension of -:math:`x`. This is done by collecting :math:`S` data points, which are -:math:`{\tilde{x}},{\tilde{y}}` pairs and then finding :math:`{\theta}` -values that minimize some function of the deviation between the values -of :math:`{\tilde{y}}` that are measured and the values of -:math:`g({\tilde{x}};{\theta})` for each corresponding -:math:`{\tilde{x}}`, which is a subvector of the vector :math:`x`. Note -that for most experiments, only small parts of :math:`x` will change -from one experiment to the next. + \boldsymbol{y}_i & = \boldsymbol{f}\left(\boldsymbol{x}_{i}, \boldsymbol{\theta}\right) + + \boldsymbol{\varepsilon}_i \quad \forall \; i \in \{1, \ldots, n\} + +where :math:`\boldsymbol{y}_{i} \in \mathbb{R}^m` are observations of the measured or output variables, +:math:`\boldsymbol{f}` is the model function, :math:`\boldsymbol{x}_{i} \in \mathbb{R}^{q}` are the decision +or input variables, :math:`\boldsymbol{\theta} \in \mathbb{R}^p` are the model parameters, +:math:`\boldsymbol{\varepsilon}_{i} \in \mathbb{R}^m` are measurement errors, and :math:`n` is the number of +experiments. The following least squares objective can be used to estimate parameter values assuming Gaussian independent and identically distributed measurement -errors, where data points are indexed by :math:`s=1,\ldots,S` +errors: .. math:: - \min_{{\theta}} Q({\theta};{\tilde{x}}, {\tilde{y}}) \equiv \sum_{s=1}^{S}q_{s}({\theta};{\tilde{x}}_{s}, {\tilde{y}}_{s}) \;\; + \min_{\boldsymbol{\theta}} \, g(\boldsymbol{x}, \boldsymbol{y};\boldsymbol{\theta}) \;\; -where :math:`q_{s}({\theta};{\tilde{x}}_{s}, {\tilde{y}}_{s})` can be: +where :math:`g(\boldsymbol{x}, \boldsymbol{y};\boldsymbol{\theta})` can be: 1. Sum of squared errors .. math:: - q_{s}({\theta};{\tilde{x}}_{s}, {\tilde{y}}_{s}) = - \sum_{i=1}^{m}\left({\tilde{y}}_{s,i} - g_{i}({\tilde{x}}_{s};{\theta})\right)^{2} + g(\boldsymbol{x}, \boldsymbol{y};\boldsymbol{\theta}) = + \sum_{i = 1}^{n} \left(\boldsymbol{y}_{i} - \boldsymbol{f}(\boldsymbol{x}_{i};\boldsymbol{\theta}) + \right)^\text{T} \left(\boldsymbol{y}_{i} - \boldsymbol{f}(\boldsymbol{x}_{i};\boldsymbol{\theta})\right) 2. Weighted sum of squared errors .. math:: - q_{s}({\theta};{\tilde{x}}_{s}, {\tilde{y}}_{s}) = - \sum_{i=1}^{m}\left(\frac{{\tilde{y}}_{s,i} - g_{i}({\tilde{x}}_{s};{\theta})}{w_i}\right)^{2} + g(\boldsymbol{x}, \boldsymbol{y};\boldsymbol{\theta}) = + \frac{1}{2} \sum_{i = 1}^{n} \left(\boldsymbol{y}_{i} - \boldsymbol{f}(\boldsymbol{x}_{i};\boldsymbol{\theta}) + \right)^\text{T} \boldsymbol{\Sigma}_{\boldsymbol{y}}^{-1} \left(\boldsymbol{y}_{i} - + \boldsymbol{f}(\boldsymbol{x}_{i};\boldsymbol{\theta})\right) -i.e., the contribution of sample :math:`s` to :math:`Q`, where :math:`w -\in \Re^{m}` is a vector containing the standard deviation of the measurement -errors of :math:`y`. Custom objectives can also be defined for parameter estimation. +where :math:`\boldsymbol{\Sigma}_{\boldsymbol{y}}` is the measurement error covariance matrix containing the +standard deviation of the measurement errors of :math:`\boldsymbol{y}`. Custom objectives can also be defined +for parameter estimation. In the applications of interest to us, the function :math:`g(\cdot)` is usually defined as an optimization problem with a large number of (perhaps constrained) optimization variables, a subset of which are -fixed at values :math:`{\tilde{x}}` when the optimization is performed. -In other applications, the values of :math:`{\theta}` are fixed +fixed at values :math:`\boldsymbol{x}` when the optimization is performed. +In other applications, the values of :math:`\boldsymbol{\theta}` are fixed parameter values, but for the problem formulation above, the values of -:math:`{\theta}` are the primary optimization variables. Note that in +:math:`\boldsymbol{\theta}` are the primary optimization variables. Note that in general, the function :math:`g(\cdot)` will have a large set of -parameters that are not included in :math:`{\theta}`. Often, the -:math:`y_{is}` will be vectors themselves, perhaps indexed by time with -index sets that vary with :math:`s`. +parameters that are not included in :math:`\boldsymbol{\theta}`.