Notebook S05. Size Distribution Inversion Using Regularization
\n",
"\n",
" Notebook Version 2: This notebook has been updated to reflect changes made to the package DifferentialMobilityAnalyzers.jl to work with the Julia v1 series (tested with Julia 1.1.0). To read the original supplement published with the paper please switch to v1.0.0 of the package DifferentialMobilityAnalyzers.jl and/or download the virtual machine on zenodo.org which contains a complete installation that works with Julia 0.6.4 \n",
"\n",
"\n",
"\n",
"This notebook demonstrates how to invert a size distribution from a measured noisy response function. The notebook is a supplement to the manuscript \n",
"\n",
"Petters, M. D. (2018) A language to simplify computation of differential mobility analyzer response functions, Aerosol Science & Technology. \n",
"\n",
"### Problem Statement\n",
"\n",
"Let R be a measured response vector from a DMA. What is the measured true size distribution?\n",
"\n",
"\n",
" Figure 1. Convolution of true size distribution resulting in the response vector.\n",
"\n",
"For situations where instrument noise is neglible, $\\epsilon_i \\approx 0$ and \n",
"
$ \\rm{π£} = \\bf{A}\\rm{π}$
\n",
"In this case, the number concentration is readily obtained from a simple inversion\n",
"
$ \\rm{π} = \\bf{A^{-1}}\\rm{π£}$
\n",
"This linear matrix inversion approach can be sufficient in some cases, However, even seemingly small noise is amplified leading to unusable estimates of $\\rm{π}$ (e.g., Kandlikar and Ramachandran, 1999). Some method to dampen the noise is needed. This is a achieved via regularization. Several regularization approaches are available for this task (Kandlikar and Ramachandran, 1999). Here the Twomey inverse (Twomey, 1963) together with the L-curve method (Hansen, 2000) to find the optimal regularization parameter is applied. This Notebook demonstrates implementation of this approach using synthetic data. Application to actual data is presented in Notebook S6. \n",
"\n",
"Routines for regularization are contained in regularization.jl. The method reproduced below is taken from Hansen (2000). The inverse of π£ is found as follows. The residual norm and size of the regularized system are defined as \n",
"
\n",
"where $\\left\\lVert x \\right\\rVert_2$ denotes the 2-norm of a vector ($x$), $L_1$ and $L_2$ are the residual norm and size of the regularized system, $\\bf{A}$ is the forward matrix, $\\bf{L}$ is a weights matrix, and $π_i$ is the initial guess of the correct solution. The general inverted size distribution $π_{inv}$ is found via \n",
"
\n",
"where $\\lambda$ is the regularization parameter and $\\arg \\min$ indicates the solution from a minimization algorithm. For the special case where $\\bf{A}$ is square and $\\bf{L} = \\bf{I}$ equals the identity matrix, $π_{inv}$ can be obtained without the need for computationally expensive minimization\n",
"
\n",
"where $\\bf{A}$ is the forward matrix defined in module DifferentialMobilityAnalyzer.jl. The convolution is applied to the dN/dlnD and N fields of the size distribution. The resulting π£ distribution is what an ideal instrument would measure. In practice, counting uncertainty imposes noise on the response. Counting uncertainty is simulated using a Poisson process. The actual number of counts in each bin are \n",
"
$c_i = N_iQ_{CPC}t$
\n",
"where $N$ is the number concentration in the bin, $Q_{CPC}$ is the flow rate of the CPC and $t$ is the integration time over the bin. The integration time is the ratio of scan time and number of bins. The noisy response function is obtained by computing the expected number of counts in a bin ($c_i$) followed by drawing a random number from a Poisson distribution with expectation value $c_i$, i.e. $c_\\epsilon = P(\\mu=c_i)$, where $\\mu$ is the mean of the Poisson distribution, and $c_\\epsilon$ are the noisy counts. The noise perturbed count distribution is then converted back to concentration and denoted as \n",
"