The invention relates to analysis of uncertainties in a system. More particularly, the invention provides a method and a system for analyzing uncertainties for a set of modules in a system in an integrated manner.
Major challenges facing industry, particularly manufacturing industries, include reducing lengthy time to market and improving the performance of existing capital assets. For example, in the case of the chemical industry, reducing the typical 5-7-year development cycle for a product may result in significant advantages in the market. In industries with relatively short cycles, enormous competitive pressures remain to accelerate the development process.
The development or improvement of a production facility generally involves several basic phases. These phases may include a technical feasibility analysis, detailed studies of the processes, pilot scale testing, detailed engineering design, building a facility, and continuous improvement of the facility. Many commercial software packages are available for various industries to assist in many of these phases. For example, for the chemical industry, computational fluid dynamics simulation packages are readily available. Further, project scheduling software packages are available for general and specific scheduling.
One concern in each phase of the development cycle is the level of uncertainties involved. The commercial packages may generally provide a point solution for a set of inputs. In order to account for uncertainties at each level, an uncertainty analysis may be required for each step or process. Such an uncertainty analysis may be required to determine the source of variations in the result of each step or process.
Uncertainty analyses may be performed using many known methods. For example, a Monte Carlo analysis may be performed for each step or process of a system. A Monte Carlo analysis may require a large number of simulations to be executed with the inputs being varied according to their underlying probability density function. The result of the Monte Carlo analysis is a distribution of the results as a function of the variations in the inputs. On a large-scale project, however, such an analysis may be cumbersome for some applications.
U.S. Pat. No. 6,173,240 discloses a method by which Monte Carlo sampling may be reduced. However, such an analysis provides results for only a single step.
The disclosed systems and methods are directed to analysis of uncertainties in a system. Uncertainties in the inputs of a system and their effect on the outputs may be efficiently analyzed by, for example, generating a simplified, yet accurate, model of the system.
Additionally, the uncertainties in several components of the system may be analyzed together, rather than individually, thereby allowing an efficient analysis of the system as a whole.
According to an aspect of the invention, a method of analyzing uncertainties in a system having at least two modules includes propagating an uncertainty distribution associated with each of a set of inputs through a module to produce a description of the uncertainty in a set of outputs of said module.
Uncertainties may be uncontrollable variations in the inputs that may cause variations in the outputs. Uncertainties may be distributed continuously or discretely over a range of Values.
A module may be any component of a system of processes, mechanisms, or algorithms. A module may include a process, a sub-process, a mechanism, an algorithm step, a calculation, or a software package simulation. Further, a module may be a part of or one or more processes, sub-processes, mechanisms, algorithm steps, calculations, simulations or other components.
Inputs are parameters that are used by one or more modules. Inputs may include, for example, internal or external parameters that may be preset, provided by a user, or provided by another module.
Outputs are parameters that are generated by one or more modules. Outputs may include parameters that are generated by a module in response to one or more inputs.
The method further includes generating a probabilistically equivalent model of the module, the equivalent model producing a model of the outputs.
The probabilistically equivalent model may be a model of a module that is less complex yet produces similar outputs for a given set of inputs. Thus, the model of the outputs generally approximates the set of outputs.
The method further includes providing the model of the outputs in a common data architecture for use as inputs by any other module in the system.
The common data architecture may be a format for presenting the data to any other module in the system in such a manner that it is readily acceptable, including any information regarding uncertainty distribution of a particular variable.
According to another aspect of the invention, a method of analyzing uncertainties in a system includes substituting at least one of a plurality modules of a system with a corresponding probabilistically equivalent module model, the equivalent module model adapted to propagate uncertainties in inputs of the module to outputs of the module. The method further includes providing outputs of each of the modules in a common data architecture for use as inputs by any other module, the architecture adapted to propagate uncertainties in the outputs to the inputs of the other module. The method further includes substituting the plurality of modules with a single probabilistically equivalent system model for propagating uncertainties in system inputs to system outputs. The single probabilistically equivalent system model may be a single, less complex module that approximates the outputs, for a given set of inputs, of a system having two or more modules.
In another aspect of the invention, a system for generating an uncertainty analysis includes a module adapted to receive a set of inputs and to produce a set of outputs as a function of the inputs. Each of the inputs has an associated uncertainty distribution. As discussed above, the uncertainty distribution may be uncontrollable variations in the input parameter. The system may further include means for propagating the uncertainty distribution of the inputs through the module to produce an uncertainty in the outputs. The means for propagating uncertainties through the module may be a process or algorithm for determining the effects of the input uncertainties on the outputs, and may include, for example, a Monte Carlo or Pattern Search analysis. The system further includes means for generating a probabilistically equivalent model of the module, the equivalent model producing model outputs. The model outputs may be a set of outputs that approximate the outputs of the module given a set of inputs. The system further includes means for providing the outputs in a common data architecture for use as inputs by any other module in the system.
In a further aspect of the invention, a system of analyzing uncertainties in a system comprises means for generating a probabilistically equivalent module model for at least one of a plurality modules of a system. The equivalent module model is adapted to propagate uncertainties in inputs of the module to outputs of the module. The system further includes two or more interacting modules and means for providing outputs of each of the modules in a common data architecture for use as inputs by any other module. The architecture is adapted to propagate uncertainties in the outputs to the inputs of the other module. The system further includes means for generating a single probabilistically equivalent system model for the plurality of modules for propagating uncertainties in system inputs to system outputs.
According to a further aspect of the invention, a system for generating an uncertainty analysis includes a modeling module adapted to receive a set of inputs and to produce a set of outputs as a function of the inputs. Each of the inputs has an associated uncertainty distribution. The system includes an uncertainty propagation module adapted to propagate the uncertainty distribution of the inputs through the modeling module to produce an uncertainty in the outputs. An equivalent model generation module is adapted to generate a probabilistically equivalent model of the modeling module, The equivalent model produces model outputs. The system further includes an output generation module adapted to provide the outputs in a common data architecture for use as inputs by any other module.
According to a still further aspect of the invention, a system of analyzing uncertainties in a system comprises an equivalent model generation module adapted to generate a probabilistically equivalent subsystem model for at least one of a plurality of subsystems, the equivalent subsystem model being adapted to propagate uncertainties in inputs of the subsystem to outputs of the subsystem. The system further includes an output generation module adapted to provide outputs of each of the subsystems in a common data architecture for use as inputs by any other subsystem. The architecture is adapted to propagate uncertainties in the outputs to the inputs of the other subsystem. The output generation module may be a module adapted to generate output in a predetermined format which, for example, includes a readily acceptable means of propagating uncertainty information. The system also includes an equivalent system generation module adapted to generate a single probabilistically equivalent system model for the plurality of subsystems for propagating uncertainties in system inputs to system outputs.
In a yet further aspect of the invention, a program product comprises machine readable program code for causing a machine to perform method steps. The program product may be, for example, a software package adapted to run on a computer, PC, laptop, mainframe or similar computing device. The program product may contain instructions to be executed. The instructions may include a list of the method steps. The method steps include propagating an uncertainty distribution associated with each of a set of inputs through a module to produce an uncertainty in a set of outputs of the module. The method steps further include generating a probabilistically equivalent model of the module, the equivalent model producing a model of the outputs. The method steps also include providing the model of the outputs in a common data architecture for use as inputs by any other module in the system.
According to another aspect of the invention, a program product comprises machine readable program code for causing a machine to perform method steps, which include substituting at least one of a plurality modules of a system with a corresponding probabilistically equivalent module model. The equivalent module model is adapted to propagate uncertainties in inputs of the module to outputs of the module. The method steps also include providing outputs of each of the modules in a common data architecture for use as inputs by any other module. The architecture is adapted to propagate uncertainties in the outputs to the inputs of the other module. The method steps further include substituting the plurality of modules with a single probabilistically equivalent system model for propagating uncertainties in system inputs to system outputs.
In a preferred embodiment, the probabilistically equivalent model is a deterministically equivalent model. Similarly, the probabilistically equivalent system model may be a deterministically equivalent system model. A deterministically equivalent model may be developed using the steps described herein. The deterministically equivalent model may be a reduced-order model, which is less complex than the actual module in that relatively few inputs may be considered in generating the model outputs.
In a preferred embodiment, propagating the uncertainty distribution includes using a Monte Carlo or Pattern Search method. Monte Carlo and Pattern Search methods are well known in the art and may include perturbing each of a plurality of variables to obtain an output uncertainty.
At least one of the set of outputs may be incorporated into at least one of the set of inputs in a feedback loop. The feedback loop allows using an output of a module to determine one or more of the inputs of the module in, for example, an iterative process.
In a preferred embodiment, an optimization module is provided for optimizing an objective function. The optimization module is adapted to receive the system outputs and to vary the system inputs. The optimization module may be a software package or a routine for either maximizing or minimizing an objective function. The objective function may be any parameter or combination of parameters whose value is desired to be either minimized or maximized. In a preferred embodiment, the objective function is a weighted function of two or more output parameters. Thus, the variable to be minimized or maximized may be a combination of several parameters.
In the following, the invention will be explained in further detail with reference to the drawings, in which:
In one embodiment of the invention, each module 16a-g may be replaced with an equivalent representation. The representation is preferably a probabilistically equivalent model. Such models may be generated according to the method described below with reference to
Now, with reference to
A wide variety of engineering and problems can be described by systems of algebraic or differential equations of the form:
where N is a model that takes as input a set of m parameters θ={θ1, θ2, . . . , θm}, that
might include, for example, reaction rate constants, initial concentrations or stoichiometric coefficients and produces as output an n-dimensional vector of state variables y={y1, y2, . . . , yn} that may be typically associated with; for example, species concentrations. There are three essential levels at which the parameter vector θ influences the model predictions y(θ). The first and easiest is the solution of the model itself given a nominal set of parameter values {circumflex over (θ)}. There are numerous tools available to accomplish this task (e.g., Kee et al. 1996). A slightly harder problem is to assess the sensitivity, S, of differential changes in y around a nominal point {circumflex over (θ)}. In this case both the model (1) and the system of adjoint sensitivity equations:
must be solved. Again, there are robust methods (e.g., Kee et al. 1996; Dunker 1984, Kramer et al. 1984) for solving (1) and (2) and, once the sensitivities have been found, they can be used to rank the relative importance of different parameters. (See, for example, Gao et al. 1995). The third, and most difficult, level, is to determine the global response of the model when the parameters are varied over a much wider range. In practice, not all values of the parameters may be equally likely, and the challenge is to combine the model response with the parameter variability.
As a way of illustrating some of the complexities associated with incorporating uncertainties consider a simple first chemical decay of the form A with a reaction rate k. The kinetics of the concentration of a species A can be described by a first order differential equation:
where y0 is the initial condition. For this very simple case the solution and the associated sensitivity are given by:
If k is an uncertain variable described by normal probability distribution with mean of k0 and standard deviation k1 i.e. k˜N[k0, k1] then the probability density function fy(k,t)[y(k,t)] of the solution for y(k,t), when k is constant throughout the solution, but uncertain, can be found analytically:
Quite clearly even though the parameter value is normally distributed the density function for the solution is lognormal. Given the probability distribution function f it is possible to characterize the uncertainty in terms of the moments. For example, the expected value or mean of y(θ) is given by (see Papoulis, 1991):
and the r-th central moments by
For the particular case (6) the expected value is given by:
There are several points that can be drawn from this example. The first is that solution using the mean value of the rate constant is not the same as the expected value i.e. y0e−k
If the analytic solution to fy(θ)[y(θ)] is not available then the key practical problem in characterization of uncertainties is evaluating the multi-dimensional integrals (7-8). A wide variety of methods have been developed and one of the simplest is the classical Monte Carlo method where the multi-dimensional integral is replaced by a finite summation of the form:
where y(θi) is the model prediction corresponding to the i-th sample point drawn from the distribution fθ(θ) and Ns is the number of sample points needed to achieve statistically stable estimates of the moments. Although Monte Carlo methods (MCM) can be used for dealing with implicit models, these methods can be prohibitively expensive, especially when the computational cost is already high. Clearly alternative approaches, which can produce results at less computational cost, are of great interest.
Some of the methods that have been developed to treat this problem include the perturbation method (Lax, 1980), the method of moments (Morgan et al., 1992), Neumann expansions (Adomian, 1980; Ghanem and Spanos, 1991), the hierarchy method (Lax, 1980), the semi-group operator method (Serrano and Unny, 1990), and the spectral-based finite element method (Ghanem and Spanos, 1991). In order to use these methods the mathematical models must be explicit functions of the parameters and the equations must be readily for manipulation. For many practical problems these constraints can be very restrictive. Some of the sampling based methods that use solutions to the models that have been developed include the stratified or pattern search methods such as the Latin Hypercube Sampling (LHS) (McKay et al., 19769; Derwent, 1987), the Fourier Amplitude Sensitivity Test (FAST) (Cukier et al., 1973, 1975, 1978; McRae et al., 1982; Koda et al. 1979), and the Walsh amplitude sensitivity procedure (WASP) (Pierce and Cukier, 1981). In practice even using the best sampling procedures described in the previous section the number of runs needed to achieve stable statistics can be prohibitively expensive.
Traditionally, the approach to the treatment of uncertainty has been to first build the model and then probe its response by varying the parameters. An alternative approach is to integrate the uncertainty at the outset. In a classic paper Wiener (1938) developed the optimal representation of a random variable in terms of a series called a “polynomial chaos” expansion (PCE):
where ω is the stochastic event, ai are constant coefficients and Hi are functionals whose m arguments are known probability density functions {ξ1(ω), ξ2(ω), . . . , ξm(ω)}. The polynomial chaos expansion, has the following four properties (Tatang, 1995): (1) Any square-integrable random variable can be approximated as closely as desired by a polynomial chaos expansion; (2) The polynomial chaos expansion is convergent in the mean-square sense; (3) The set of orthogonal polynomials is unique given the probability density function; (4) The polynomial chaos expansion is unique in representing the random variable. The probabilistic form (11) is analogous to a conventional Fourier series where a function is expanded in terms of a linear combination of sine and cosine basis functions. In practice only a finite number of terms M in (11) are used:
Given the general form (12), the next steps are to define the functionals (Hi), functions (ξ) and solve for the coefficients ai of the finite expansion. The simplest way to determine the ai's is to use the method of weighted residuals (MWR) (See, for example, Villadsen and Michelsen, 1978). The weighted residual is defined as the difference between the exact solution and the result when the series expansion is substituted into the model. For the general form (1) the j-th weighted residual is given, after a suitable change of variables from θ→ξ by
Rj(ξ)={N[ŷ(ξ)]ŷ(ξ)−f(ξ)}Wj(ξ);j=1, 2, . . . , m (13)
where Rj(ω) is the j-th residual and the Wj(ω) are weighting coefficients associated with each of the uncertain parameters in the model. If the expansion ŷ(ξ) satisfies (13) exactly then the residual is zero. Depending on the choice of weighting function and minimization method used to find the coefficients ai the method is known as a least squares, Galerkin, or collocation based MWR schemes.
In this case the coefficients ai are determined by setting the residual to be orthogonal to the space spanned by the probabilistic basis functions used in the expansion. The probabilistic form of the inner product of the residual and the weighting function, Wk(ξ), is set to zero:
The integral (14) is defined for each of the M+1 basis polynomials Hi. Once the integrals have been evaluated the system of M+1 deterministic equations can then be solved simultaneously for the coefficients ai. Two weighting functions are typically used in practice a Galerkin and a collocation formulation.:
In the Galerkin case the orthogonal trial functions are used as the weighting functions. When collocation is used δj(ξ−c) are Dirac delta functions which force the residual to vanish at the collocation points c={c1, c2, . . . , ck). While in either case the multi-dimensional integrals (14) need to be evaluated, careful choice of the functionals Hi, the weighting functions Wk and the independent functions can considerably simplify the process.
Polynomial chaos expansions are “problem specific” because of the definition of orthogonality in stochastic systems. Similar to the concept of orthogonal vectors spanning the vector space, parameter specific orthogonal polynomials are derived such that their roots are spread over the high probability region of the parameter. Two stochastic functions gi(ξ) and gj(ξ) are orthogonal when their inner product, defined using the probability distribution of the stochastic variable ξ, vanishes The definition of orthogonal polynomials is:
where gi(x) is the i-th order orthogonal polynomial. Note that the polynomials are derived solely from the probability density function of the model parameters. In general, problem-specific orthogonal polynomials can be derived by algorithms such as ORTHPOL, following the recurrence relations (Gautschi et al. 1994):
g−1(x)=0,
g0(x)=1,
gk+1(x)=(x−αk)gk(x)−βkgk−1(x),
k=0, 1, . . . , n (17)
where the coefficients αk, βk can be expressed in terms of the orthogonal polynomials following the Gram-Schmidt orthogonalization procedure:
The inner product used above is in the form of Riemann-Stieltjes integral
where the function λ(x) is the indefinite integral of the weighting function. Several different types of orthogonal expansions are summarized in Table 1.
As an illustration of the process consider the simple case
described earlier. The basic idea is to approximate y(t) using a polynomial expansion of the form:
where the gi(ξ) are the basis functionals and yi(t) are the time varying coefficients in the expansion. For the particular case of Hermite polynomials the expansion is of the form:
y(t)=y0(t)+y1(t)ξ+y2(t)(ξ2−1)+y3(t)(ξ3−3ξ)+y4(t)(ξ4−6ξ2+3)+ (21)
Applying the variational procedure described in the previous section produces a set of linear ordinary differential equations for the coefficients:
where A is the identity matrix and elements of B for the first four terms in the expansion is given by:
The key point to note about (22) is that the equations for the uncertainty coefficients are of the same structural form as the original model and so its numerical solver can be used for both forms.
In the collocation approach the residual (14) is forced to vanish at ck, the collocation points thus satisfying the model exactly at ξ=c1, ξ=c2, . . . , ξ=cM+1. For an M-th order polynomial chaos expansion, the collocation points {ck} are the roots of gM+1(ξ). Collocation points are chosen in a manner analogous to the Guassian quadrature method for evaluating integrals. In the collocation method, instead of solving once a large system like (22), the deterministic model is solved M+1 times at each of the collocation points ck. The result is a set of M+1 deterministic equations for different Ck:
After the model has been solved at each of the collocation points the set of simultaneous linear equations (24) can be solved for the coefficients y0, . . . yM. A key advantage of the collocation procedure is that it can be applied to “black box” type models where the model equations are not known explicitly because the method it requires only the solution of the model at specific values of the parameters.
This method, and the associated properties are completely generalizable to systems with many stochastic parameters. For example, if the parameters are independent:
Assuming y is a function of N independent random variables, y=y(ξ1, ξ2, . . . , ξm), an M-th order polynomial chaos approximation y(ξ1, ξ2, . . . ξm) of y is written as:
The choice of collocation points for higher order system warrants further discussion. Unless all the cross product terms are included in the expansion, only selected collocation points will be used to determine the PCE coefficients. In order to handle this situation a formal procedure has been developed to choose systematically the collocation points used in the solution procedure. Consider first a two parameter case. The collocation points for each parameter are placed in order of decreasing probability. In the case when the probability is equal (e.g., in a uniform distribution), the points are organized in increasing distance from the mean. The first pair of points, which contains the most probable values for all the parameters among the collocation points (c1, C3), is termed the anchor point (ξAnchor). For each increasing order of approximation, the corresponding variable's collocation point is perturbed. Therefore, the pairs of points (c1, C3), (C2, C3), and (c1, C4) are chosen for an approximation which has a constant term and the first order terms in ξ1 and ξ2. If the there is a bilinear term g1(ξ1)g1(ξ2) is used in the approximation, the point (C2, C4) will also be used in the coefficient evaluation process.
Given the discussion in the previous section there is a clear need for an automatic procedure to simplify the choice of the appropriate numbers of terms in the expansion of the model output variables. Using an error correction mechanism embedded into most ordinary differential equations solvers the truncation error of the response surface representation is estimated by comparing the M-th order prediction to the (M+1)-th order prediction. The model is evaluated at the collocations points corresponding to the (M+1)-th order approximation and then the model solutions are compared to the approximation obtained from the M-th order PCE at those points. The error at each of the (M+1)-th order collocation points is defined as the square of the distance between the exact solution and the M-th order approximation:
Two specific metrics are used; the sum square root (SSR) error and the relative sum square root (RSSR) error as:
The error measures in (28) can be used to guide the decision of whether more terms are needed in the PCE. The accuracy and number of terms required for the response surface approximation depends on the goal of the analysis. This procedure is implemented in a computer program that guarantees the convergence of the PCE series with increasing order. Interactions between the parameters can also be elucidated. The order of approximation is increased until the error is negligible. However, excessive number of model runs to evaluate coefficients sometimes can makes this approach computationally intensive. It is possible to analyze the error contribution from each of the variables by evaluating the individual terms, and select variables that contribute to the error as targets for higher order representation. Physical insights can also be used to guide the selection and use of cross product terms. One procedure for error control is shown in
Once the coefficients in the polynomial chaos expansion have been determined there are several other useful properties than can be determined including the probability density function of the outputs, confidence intervals, moment information, and variance apportionment to identify the critical input variables. For example, one simple way to obtain the probability distribution of a response variable from the PCE representation is by Monte Carlo (MC) sampling of the expansion itself. In essence the PCE approximation can be viewed as a reduction of the original output variable. Where MC sampling of the original complex model is prohibitively expensive, MC sampling of a linear combination of algebraic terms containing the random input variables provides a viable alternative for understanding the behavior of the random output variable. This method can als be used to derive the cumulative density function (CDF). To generate a CDF, the output samples are sorted in ascending order and the limits of each fractile are recorded. The confidence intervals can also be determined using the sorted samples. For example, a 90% confidence interval will be the range of the empirical samples after the highest and lowest 5% of the samples are discarded. If a probability density function is needed, the range of the response variables is divided into bins or intervals and the frequency of occurrence in each interval is counted based on the same procedures used to generate histograms.
One application of particular importance is the determination of the moments of the output probability distribution and their application to the analysis of variance. The moments of the distribution can be determined empirically from the MC samples; or they can be calculated directly from the PCE coefficients, using the definition of the n-th central moment (cmn). The evaluation of moments is simplified by the orthogonal properties of the polynomials. For example, if:
y=u0+u1g11(ξ1)+u2g21(ξ2), (29)
the mean value is equal to y0, and the variance of the random variable is described by
σ2=A1u12+A2u22
A1=∫g112(ξ1)pξ(ξ1)dξ1, A2=∫g212(ξ2)pξ(ξ2)dξ2. (30)
Higher moments can also be determined from the coefficients of higher order terms. The relationship between the PCE coefficients and the variance suggests the utility of the PCE approximation for variance analysis. The contribution of each input parameter can be determined from the relevant terms in the approximation. In (30), the variance contribution (VC) from ξ1 is A1 u12, while the VC from ξ2 is A2 u22. Any cross terms are apportioned among the variables involved. This kind of analysis is particularly useful for identifying input variables whose uncertainties have strong effects on the uncertain outputs.
Consider an simple series reaction mechanism of the form
where k1 and k2 are uncertain parameters described by the normal distributions k1=N[0.5,0.1] and k2=N[2.0,0.5]. The initial conditions are [A(0)]=100, [B(0)]=[C(0)]=0. Once the reactions commence, the concentrations of A, B, and C are uncertain because of the uncertain rate constants. Set out below are the steps in applying the collocation procedure for uncertainty analysis.
Step 1. Specify Uncertain Parameters. In this example, the probability distributions of k1 and k2 are assumed to be independent and Gaussian. The polynomial chaos expansions are simply:
k1=k10+k11ξ1
k2=k20+k21ξ2 (31)
where k10=0.5, k20=2.0 are the mean values of k1 and k2, and k11=0.1 and k21=0.5 are the standard deviations. Methods for developing PCE forms for other probability distributions are described in Tatang (1995).
Step 2. Generate Problem-specific polynomial chaos expansions. Since the explicit forms of distributions of k1 and k2 are known, orthogonal polynomials chaos {gi} can be generated such that the inner products, defined by ∫−∞∞gi(ξ)gj(ξ)fξ(ξ)dξ, are zero, where fξ(ξ) is the PDF of the uncertain variables ξ1 or ξ2. For standard normal distributions, the PCE are simply orthogonal Hermite polynomials defined by:
H0(ξ)=1,
H1(ξ)=ξ,
H2(ξ)=ξ2−1,
H3(ξ)=ξ3−3ξ
Etc. (32)
Step 3. Approximate Uncertain Outputs Using Polynomial Chaos Expansion. The model outputs, concentrations A (ξ1,t), B(ξ1, ξ2,t), C(ξ1, ξ2,t), are expressed as linear combinations of the orthogonal polynomials determined in Step 2. These expressions are known as the polynomial chaos expansions (PCE) for the uncertain outputs, and to first order, are given by:
The concentrations of A, B, and C, and the coefficients, A0, A1, . . . , B0, B1, . . . , C0, C1, . . . are all functions of time. At each time point, the number of coefficients, hence the number of simultaneous equations for their solution, is determined by the order of the polynomial approximation. The higher the order of the approximation, the better the approximation. In practice the procedure is to start with a low order expansion and to increase the order iteratively as needed. Linear PCE representations for A, B, and C, using Hermite polynomials, are given by:
A(ξ1,ξ2,t)=A0(t)+A1(t)ξ1+A2(t)ξ2
B(ξ1,ξ2,t)=B0(t)+B1(t)ξ1+B2(t)ξ2
C(ξ1,ξ2,t)=C0(t)+C1(t)ξ1+C2(t)ξ2. (34)
Step 4. Find the Collocation Points. Collocation points are selected to solve for the coefficients, A0, A1, A2, B0, B1, B2, C0, C1, and C2, in the approximation (33). For a linear approximation, the collocation points are determined by the roots of the second order polynomial. H2(ξ)=ξ2−1; therefore, ξ=±1. Since ξ1 and ξ2 are Gaussians and are symmetric about zero, the four pairs of collocation points (±1, ±1) are equal in probability and equal in distance to the mean (0, 0). The point (ξ1, ξ2)=(1, 1) is designated as the anchor point, the point with the highest probability. In this example, the points (−1, 1) and (1, −1) are also chosen. These correspond to (k1, k2) pairs of (k10+k11, k20+k21), (k10−k11, k20+k21), and (k10+k11, k20−k21), as listed in Table 2.
Step 5. Solve the Model at the Collocation Points. The model is formulated to take the uncertain parameters k1 and k2 as external inputs. Solutions for A(t), B(t), and C(t) are evaluated for each pair of (k1, k2) found in Step 4. The original model solver is used “as is”, since the model equations are exactly satisfied at the collocation points.
Step 6. Solve for the Coefficients of Expansion from Model Results. The model solutions for A(t), B(t), and C(t) are equated to simple algebraic equations in (27) for each of the collocation points (k1, k2) listed in Table 2. The resulting time-dependent equations for the coefficients A0, A1, A2, B0, B1, B2, C0, C1, and C2 are evaluated numerically at the selected time points. At each time point, the algebraic equations (24) are solved simultaneously for the unknowns. Since the concentration of A does not depend on the uncertain rate constant k2, the coefficient A2(t) in is exactly zero at all times.
Step 7 Estimate the Error of Approximation. The error of the linear PCE can be evaluated for each species at any given time based on the solutions at the roots of the third order Hermite polynomial (collocation points for the second order PCE). The roots to the third order approximation are ξ=0, ±{square root}{square root over (3)} and the corresponding points for each parameter combination are shown in Table 3 for a second order approximation.
An example error calculation is shown in Table 4 for species B at time 1.0 units (when the concentration of B is at its maximum).
Table 5 shows that the relative error of the linear approximation of the response surface is about 4%.
This percentage number by itself is not an absolute measure of the “goodness” of the approximation. When the expected value is close to zero, the RSSR can grow in an unbounded manner and caution should be used in interpreting the error estimates.
Step 8. Increase the Order of Approximation. One way to decrease the error of the response surface approximation, and hence of the uncertainty estimates, is to increase the order of the polynomial chaos approximation. Including higher order terms and cross product terms have the obvious utility of capturing curvature of the response surface better. There is an additional advantage. Based on the choice of collocation points, as described in Step 4, increasing the number of terms also increase the “spatial coverage” of the collocation points, making the estimate applicable over a wider range of values of the uncertain inputs. The errors associated with different orders of approximation for the concentration of B at time=1 are presented in Table 5.
Step 9: Variance Analysis. Using the formulation described in (29-30) the mean and the variance, of the spread, of the response are of particular interest.
From the individual PCE coefficients, the contribution of any particular uncertain parameter to the output variance can be calculated. The PCE coefficients give information regarding the “global sensitivity” of the response variable to the parameter.
Another application of the polynomial chaos expansion represents the distribution of the response variable as a functional of the uncertain input parameters. Monte Carlo sampling procedures can be applied to sampling the polynomial chaos expansion to obtain the probability density function of the output. With this approach, the overhead computer resource required to run a Monte Carlo analysis is small compared to the time taken to solve the model at the collocation points. For example, a 39-term PCE takes 1.8 seconds to solve and sample. If the model takes two hours to run, a overhead time of several seconds is negligible, and the time savings of using DEMM instead of Monte Carlo is the ratio of the number of model runs needed for the two methods.
In a further embodiment, after increasing the order of the approximation, the algorithm may determine whether certain inputs may be ignored due to negligible uncertainty effect. In this manner, a reduced-order, deterministically equivalent model may be achieved. As indicated in
Once a deterministically equivalent model has been created for each module, the modules may communicate with each other while maintaining proper propagation of the input uncertainties.
The geometry module 510 and the kinetics module 520 may provide inputs to a reactor model module 530 for simulating the reaction of the reactants in a chemical reactor. The reactor model may be a commercially available software package such as Chemkin. The reactor model module 530 may provide inputs to a computational fluid dynamics (CFD) module 540. The CFD module 540 may simulate the fluid dynamics of the reactants and products through a chemical reactor. A software package such as STARCD may be used to simulate the fluid dynamics. Outputs of one or more modules may be used to provide inputs to a process engineering module 550. Although
Further, an optimization module 570 may be incorporated to perform an optimization in the design of the reactor system. The optimization module 570 may also be a commercially available software package 570 and may include any of a variety of commonly known optimization algorithms, such as recursive quadratic programming or sequential quadratic programming. As indicated by the dotted lines in
As seen from the illustrated systems of
Similarly, other uncertainty distribution types may be specified for each variable. For example, the uncertainty distribution of a variable having an exponential PDF distribution may be described by providing a mean value. Other distribution types that may be described using this common data architecture may include a polynomial chaos expansion, a list of points, or a histogram. Thus, the uncertainty distribution of each variable, input or output may be associated with the variable itself in, for example, a database.
The above-described methodology has been shown to use random variables. It is contemplated within the scope of the invention to allow utilization of random processes, for example, using Karhunen-Loeve series expansions, which are well known to those skilled in the art. For details on Karhunen-Loeve series expansions, reference may be made to Papoulis, A. Probability, Random Variables, and Stochastic Processes, 3rd Edition, McGraw Hill, N.Y., 1991, as well as to Tatang, M. A., Direct Incorporation of Uncertainty in Chemical and Environmental Engineering Systems, Ph.D. Thesis, Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, Mass., 1995, each of which is hereby incorporated by reference.
With each module being able to effectively communicate its outputs to all other modules in a common data architecture, the optimization process may be automated.
In a further embodiment of the invention, the entire system of processes and subsystems may be modeled as a single system by creating a deterministically equivalent model, as described above for individual modules. In this regard, the global inputs into the system may now be treated as the inputs 10 illustrated in
Although an embodiment of the invention is described above as being applied in the chemical environment, embodiments of the invention may be employed in a broad variety of applications. Some possible areas and industries for application include, without limitation, financial analyses, oil industry, various types of networks including computer networks, transportation, circuit simulation, project scheduling, and decision analysis.
For example, in the financial arena, when alternative investment proposals are considered, there are often many uncertainties to be considered including market size, selling price, financing availability, etc. If financial risk is to be managed effectively, it is critical to be able to assess the relative contributions of different sources of uncertainties. For example, in a net present value (NPV) calculation, there are often uncertainties in the future cash flows and the discount rate that, in turn, lead to uncertainties in the project valuation. Using a method or system for uncertainty analysis according to an embodiment of the present invention, it is possible to determine the NPV probability distribution and the contributions of individual terms to the variance in the valuation. Such information is important for developing risk mitigation strategies or to determine where additional resources might be allocated to reduce the overall risk. Similar analyses may be applied to a broad spectrum of financial instruments including options pricing, portfolio management, insurance pricing, etc.
An embodiment of the application may also be applied to the area of logistics and transportation networks. A common problem in managing the logistics of moving products from factories to warehouses to markets is to manage the shipping costs when there are uncertainties in customer demands, raw material suppliers and the availability of shipping capacity. The uncertainty analysis methods and systems according to embodiments of the present invention, together the mathematical programming formulations of the logistics problem, may be used to develop robust production schedules, inventory management policies and identify optimal ways to allocate shipping.
As a further example, an embodiment of the invention may be applied to simulations of electronic circuits. Electronic circuits are typically composed of many subsystems. At the lowest level, the components might be transistors, capacitors or resistors and, at a higher level of integration, the subsystems could be amplifiers or inverters. The electrical properties of these devices can be uncertain, which in turn leads to uncertainties in the overall system performance. Current circuit simulators, such as SPICE, cannot propagate effects of component uncertainties on determining the probability distribution of predicted outputs from the simulation model. The uncertainty analysis methods and systems according to embodiments of the present invention may treat the circuit simulator as a black box and identify component uncertainties that are influencing the outputs.
Mechanical and structural analyses may also employ embodiments of the present invention. Finite elements are widely used to study the statics and dynamics of complex systems made up of simple components. The uncertainty analysis methods and systems according to embodiments of the present invention may treat uncertainties in the external loadings and physical properties and determine how they affect the predictions of the numerical model. Embodiments of the present invention, in combination with Karhunen Loeve decomposition, can also account for spatial and temporal variations in the intake parameters.
Another example of an application of an embodiment of the invention is the analysis of decisions. Structuring and analyzing a complex decision involves many uncertainties. When models are used to describe the project elements, or there are uncertainties in decision outcomes, there is a need to identify the component parts that have the most influence on the outcomes. A knowledge of the probability density function of the outcomes as described above with reference to the embodiments of the present invention enables the investor to manage the risk across a portfolio of projects.
An embodiment of the invention may be implemented on a processor such as the computer system illustrated in
A hard drive 822 may be provided to store data and instructions in a non-volatile memory, for example. Further, a random access memory 824 is provided to temporarily store instructions for the CPU. The random access memory 824 may be provided with stored instructions by, for example, an executable residing on the hard drive 822 or on an external storage medium such as a floppy disk or a CD-ROM 828.
Information on the CD-ROM 828 may be accessed by the CPU 820 via a CD-ROM drive 826. Other drives may be provided to access information from other types of external storage media.
The CPU 820 may receive instructions, commands or data from an external user through input devices such as a keyboard 830 and a mouse 840. The CPU 820 may display status, results or other information to the user on a monitor 850.
While particular embodiments of the present invention have been disclosed, it is to be understood that various different modifications and combinations are possible and are contemplated within the true spirit and scope of the appended claims. There is no intention, therefore, of limitations to the exact abstract or disclosure herein presented.