Sensitivity matrix, computation

If we consider the limiting case where p=0 and q O, i.e., the case where there are no unknown parameters and only some of the initial states are to be estimated, the previously outlined procedure represents a quadratically convergent method for the solution of two-point boundary value problems. Obviously in this case, we need to compute only the sensitivity matrix P(t). It can be shown that under these conditions the Gauss-Newton method is a typical quadratically convergent "shooting method." As such it can be used to solve optimal control problems using the Boundary Condition Iteration approach (Kalogerakis, 1983). 

Step 2. Given the current estimate of the parameters, k , compute the parameter sensitivity matrix, G the response variables f(x k[Pg.161]

At this point we can summarize the steps required to implement the Gauss-Newton method for PDE models. At each iteration, given the current estimate of the parameters, ky we obtain w(t,z) and G(t,z) by solving numerically the state and sensitivity partial differential equations. Using these values we compute the model output, y(t k(i)), and the output sensitivity matrix, (5yr/5k)T for each data point i=l,...,N. Subsequently, these are used to set up matrix A and vector b. Solution of the linear equation yields Ak(jH) and hence k°M) is obtained. The bisection rule to yield an acceptable step-size at each iteration of the Gauss-Newton method should also be used. 

The solution of Equation 10.28 is obtained in one step by performing a simple matrix multiplication since the inverse of the matrix on the left hand side of Equation 10.28 is already available from the integration of the state equations. Equation 10.28 is solved for r=l,...,p and thus the whole sensitivity matrix G(tr,) is obtained as [gi(tHt), g2(t,+1),- - , gP(t,+i)]. The computational savings that are realized by the above procedure are substantial, especially when the number of unknown parameters is large (Tan and Kalogerakis, 1991). With this modification the computational requirements of the Gauss-Newton method for PDE models become reasonable and hence, the estimation method becomes implementable. 

If the Newton-Raphson method is used to solve Eq. (1), the Jacobian matrix (df/3x)u is already available. The computation of the sensitivity matrix amounts to solving the same Eq. (59) with m different right-hand side vectors which form the columns — (3f/<5u)x. Notice that only the partial derivatives with respect to those external variables subject to actual changes in values need be included in the m right-hand sides. 

Computing the sensitivities is time consuming. Fortunately the direct integral approximation of the sensitivity matrix and its principal component analysis can offer almost the same information whenever the direct integral method of parameter estimation applies. 

Dougherty and Rabitz [8] point out that for many applications it is not necessary to compute the entire sensitivity matrix, but only those columns for species considered to be of interest, such as those susceptible to measurement. There are however, certain advantages to computing the entire Green s Function matrix, principally the ability of time scaling in cases where sensitivities are required at several points in time. [6] For the purpose of this paper, the entire matrix was computed. 

The sensitivity trajectory W t), when requested, is computed by the method of Caracotsios and Stewart (1985), extended here to include Eq. (B.l-4b). The sensitivities are computed at to and after each step, via a direct linear algorithm that utilizes the current matrix G of equation (B.l-5b) or (B.l-5a) to solve equation (B.l-4b) or a backward-difference form of (B.l-4a). 

The main task of DDAPLUS is to compute the state vector u (and the sensitivity matrix W if requested) at the next output point. The status of the solution upon return to the calling program is described by the parameter Idid, whose possible values are explained below. 

The method of the principal component analysis of the rate sensitivity matrix with a previous preselection of necessary species is a relatively simple and effective way for finding a subset of a large reaction mechanism that produces very similar simulation results for the important concentration profiles and reaction features. This method has an advantage over concentration sensitivity methods, in that the log-normalized rate sensitivity matrix depends algebraically on reaction rates and can be easily computed. For large mechanisms this could provide considerable time savings for the reduction process. This method has been applied for mechanism reduction to several reaction schemes [96-102]. 

The methods for solving multi-dimensional EM inverse problems are usually based on optimization of the model parameters by applying different inversion schemes. The key problem in the optimization technique, as we demonstrated in Chapter 5, is the calculation of the Prechet derivative (sensitivity matrix), which usually requires much computational time. 

This computer program carries out calculations for analysing the sensitivity, and also carries out the analysis in main components of the sensitivity matrix, the analysis of the reaction fluxes, the simplification of the mechanisms from the above information, and the determination of the quasi-stationary species. KINALC is a post-processor of the CHEMKIN computer programs. 

An overview of the methods used previously in mechanism reduction is presented in Tomlin et al. (1997). The present work uses a combination of existing methods to produce a carbon monoxide-hydrogen oxidation scheme with fewer reactions and species variables, but which accurately reproduces the dynamics of the full scheme. Local concentration sensitivity analysis was used to identify necessary species from the full scheme, and a principle component analysis of the rate sensitivity matrix employed to identify redundant reactions. This was followed by application of the quasi-steady state approximation (QSSA) for the fast intermediate species, based on species lifetimes and quasi-steady state errors, and finally, the use of intrinsic low dimensional manifold (ILDM) methods to calculate the mechanisms underlying dimension and to verify the choice of QSSA species. The origin of the full mechanism and its relevance to existing experimental data is described first, followed by descriptions of the reduction methods used. The errors introduced by the reduction and approximation methods are also discussed. Finally, conclusions are drawn about the results, and suggestions made as to how further reductions in computer run times can be achieved. 

The weighting factors in each case play decisive roles in determining which regions of the wavefimction converge first, and to what extent. They also determine how sensitive the computed K matrix element would be, to large emn S in the unconverged regions of the wavefunction. 

For performing the Gaufi-Newton iteration we have to integrate Eq. (7.3.2) and to compute the sensitivity matrix with respect to the initial values. In both subtasks we make use of the solution properties by applying coordinate projection. 

---