Propagating bias and precision errors using the perturbation method

PROPAGATING BIAS A N D PRECISION ERRORS USING THE PERTURBATION METHOD

/ /

Philip A. Jones Bechtel National, Inc.

Michael A. Friedman Electric Power Research Institute

Calculations based on measured variables are subject to some uncertainty due to errors in measurements. These measurement errors, while unavoidable, can be quantified. The usual method for determining the uncertainty in a calculated variable, given the uncertainties lop expected errors] in the measured variables that define it, involves the use of the Taylor series method. With this method, the partial derivative of the result with re* spect to each of the measured variables must be determined. This can be tedious and time consuming. When results are calculated by computer, another method is available that does not require partial differentiation. This alternative is known as the perturbation method. This paper explains how the perturbation method can be used in conjunction with the methodology for propagating measurement errors set forth in ANSI/ASME Performance Test Code 19.1-1985. The FORTRAN coding required to implement the perturbation ISSN0019-0578/90/04/0071/07/$2.50© ISA 1990

method is also described. The program, which propagates the measurement uncertainties, is easily combined with the results calculation code and requires little or no maintenance. The code was developed to determine uncerteinties associated with boiler performance calculations that were carried out as part of a test program for circulating fluidized bed combustion sponsored by the Electric Power Research Institute.

INTRODUCTION In general, when a test involving measured parameters is performed, these independent variables are used to calculate values for one or more dependent variables. In addition, it is often useful to calculate the uncertainty of a dependent variable, given the uncertainties of the independent variables used in its determination. The method specified for propagating measurement uncertainties in ANSI/ASME Power Test Code 19.11984 (PTC 19.1 [1]) requires that the first term of a Taylor series expansion be used to estimate the sensiISA Transactions

•

Vol. 29, No. 4

71

tivity of the dependent variable to each independent variable. To do this, the partial derivative of the dependent variable with respect to each of its independent variables must be determined analytically. When using a computer to calculate test results, much less code development is required to calculate the effect of the measurement uncertainties on the result using the perturbation method. The perturbation method has been introduced in previous papers by Moffat (1982 [2] and 1985 [3]). In this paper the perturbation method and its use with the PTC 19.1 methodology is explained, the FORTRAN coding required to propagate bias and precision errors by this method is also described.

each P~, and is found by evaluating the partial derivative of R with respect to P~: 0 i = dR/dP i

In Eqs. (2) and (3), the terms 0; x SPi and 0~ x BP~ are first term approximations from the Taylor series representations for SRPi

=f(Pl, P2. . . .

Pi + SPi . . . . PJ) -f(P1, P1 . . . . Pi . . . . PJ)

R = f(Pl, P2. . . .

Pg)

(1)

where P1 through Pj are the average values for the measured variables. Associated with the average value of each measured variable is a precision index, SP~, and a bias error, BP r The precision index quantifies the random errors observed in repeated measurements of a given independent variable. The bias error represents the systematic or fixed error associated with the measurement of a given independent variable. The bias error for each measurement of a given independent variable is constant. According to the methodology set forth in PTC 19.1, the precision indices and bias errors of the independent variables must be propagated separately to obtain a precision index and a bias error for the dependent variable. The resulting precision index and bias can then be combined to calculate the total error of the dependent variable. The propagation equations given in PTC 19.1 are as follows:

s

B =

(0 i x SPi)2

(0; x Bpi)2

(2)

(3)

(5)

and BRPi = f(P1, P2 . . . . Pi + BPi . . . . P J) -f(Pp P1 . . . . Pi . . . . PJ)

THE P£wrlJlzgI~ATIOIII Mk'THO£1 A dependent variable (or result), R, is a function of J independent (measured) variables:

(4)

(6)

where SRP~ and BRP i represent quantities by which the result changes when the ith input parameter is changed by its precision index and bias error, respectively. The terms 0; x SP i and 0~ x BP~ are, therefore, approximately equal to SRP~ and BRP~ as represented in Eqs. (5) and (6). To use the perturbation method for the propagation of precision and bias errors, a computer code must be developed to calculate the result, R, a total of 2 x J + 1 times. The first evaluation of R is performed using the average values of the measured parameters. R' is then evaluated J times, perturbing a different measured parameter by its precision index and subtracting the original, unperturbed R from each of the perturbed Rs. This process yields SRP 1 through SRPj. The process is repeated an additional J times, perturbing the averaging values of the measured parameters by their bias errors to yield BRP 1 through BRP~. The precision index and bias error of the result are calculated as follows:

S = ~i~-~l SRPi2

(7)

B =

(8)

BRPi 2

The term 0; x SP i also appears in the Welch-Satterthwaite formula for V, the degrees of freedom of S: 4

The terms S and B are the precision index of the result and the bias error of the result, respectively. The term 0; is the sensitivity coefficient associated with 72

ISA Transactions • Vol. 29, No. 4

V =

(9) J i=1

((o, x SC)'/v C)

where VPi is the number of degrees of freedom for the ith measured variable. SRPi may be used once in lieu of the term 0. x SP~, yielding the following equation for V: Sr4

(10)

J ~,(SRpi4/Vp i) i=1 Once V has been found, the appropriate Student t value for a given number of degrees of freedom can be obtained from a two-tailed Student t table. The precision index multiplied by the appropriate value of t is equal to the precision error. With the precision error known, the uncertainty in the result U can be calculated with the following equation: U =4B 2+(txS) 2

(11)

THE FORTRAN CODING A flow chart representing the FORTRAN code required to implement the Perturbation Method is shown in Figure 1. In block (a) of the flowchart, the arrays required for the uncertainty analysis program are dimensioned. These arrays, P, BP, SP, and VP, must be dimensioned for the number of independent variables. For each dependent (measured) variable, these arrays contain values for the average of the measurements (P), for the bias error of the measurements (BP), for the precision index of the measurements (SP), and the number of degrees of freedom (VP). These values comprise the set of inputs to this program and are initialized in block (b). Also initialized are J, which is the number of independent variables; I, which is a counter and is initially set to zero; and SR and BR, the precision index and bias error of the result, which should both initially be set to zero. In block (c) a call is made to the subroutine CALC. The P array (containing the set of measured values) is passed as an argument to the CALC subroutine. The function defining the dependent variable is evaluated by the CALC subroutine, and the result, R, is passed back as an argument. Blocks (d) through (q) comprise a loop which is executed J times. In block (d), the counter I is incremented by one, so that a different one of the J independent variables is perturbed each pass through the loop. In block (e), the Ith independent variable is per-

turbed by its corresponding precision index, so that the subroutine CALC can be called again in block (f). The difference between the perturbed result, R' (returned by subroutine CALC), and the original, unperturbed R determined in block (c) is found in block (g). This difference, SRP, is squared and then accumulated in the variable SR in block (h). In block (i), the Ith element of the P array is reset to its original value. If the number of degrees of freedom associated with the precision index of the lth input variable are greater than zero, the command represented by block (j) results in the execution of block (k). Block (k) represents the accumulation of the value that will be used as the denominator term (DENV) in the calculation of the number of degrees of freedom for the result (performed in block (t)). If the Ith value of VP is equal to zero, this step is not executed, in accordance with block (j). The logic behind blocks (1) through (p) is identical to that of blocks (e) through (i). However, in blocks (1) through (p), the bias error is being propagated rather than the precision index. In block (q), the value of I is checked against J, If the value of I is not yet equal to J, the loop will be executed again. Blocks (r), (s), and (t) use the values accumulated in the variables SR, BR, and DENV to calculate SR, BR, and VR--the precision index, bias error, and degrees of freedom of the result, respectively. Once the degrees of freedom for the result are known, a call to subroutine STUDENT in block (u) can be used to determine the appropriate Student t value, T. The value of VR is passed as an argument to subroutine STUDENT, and the corresponding value for T is passed back as an argument. Subroutine STUDENT is simply a look-up table. With BR, T, and SR all known, UR, the uncertainty of the result, can be found as shown in block (v).

A PERTURBATION METHOD EXAMPLE To illustrate the use of the perturbation method, consider the example problem that is worked out using the Taylor series method in Section 5.2 of PTC 19.1. In this example, the object is to calculate the uncertainty associated with a water flow measurement made using a 6 in. by 4 in. venturi. The equation specified for determining the mass flow of the water is:

MDOT =

0.099702 x C x Yx DT2 x F x ~ x / = x / 1 - (DT[DI)4

ISA Transactions • Vol. 29, No. 4

H (12)

73

/

[

DIMENSION P(J}, ~ (a) sp(J) ,sP (J), vp (d~j

I INITIALIZE P,BP, I (b) SP,VP,J,SR,BR, I

P (1)=P(1)+BP(1)

(i)

CALL CALC(P,R')

(m)

1 C~L CAT.C(P,R) I (c)

C I=l+l

BR=BR+BRp2

(d)

l (o)

P(1)=P(1) -BP(1) l (P)

CALL CALC(P,R')

}

(f) FALSE

I

SRP=R'- R

(g) 4 ~ TR~

l

SR=SR~+SRP2

(h)

SR= ~-ff

1

FALSE

(q) (r)

BR=B/~ "

(s)

VR=S R4/DENV

(t)

(j) {u)

I DENV=DENV +s~4/vp (I)

©

(k)

I ur~.~/,3,~2+crxsRI2

@

Figure 1. Flowchart for Perturbation Method Coding 74

ISA Transactions •

Vol. 29, No. 4

(v)

where:

MDOT

= mass flow rate of water lb m/s

C

= discharge coefficient, dimensionless

Y

= expansion factor for water, dimensionless

DT

=

throat diameter, in.

DI

=

inlet diameter, in.

F

= thermal expansion number for water, dimensionless

G

=

H

= differential pressure head, in. H20

water density, lb m/ft3

Table 1 shows nominal values for the independent variables along with the absolute bias errors, absolute precision indices, and degrees of freedom associated with each independent variable. These numbers comprise the sets of inputs to the perturbation method program. The value of MDOT for this set of inputs is calculated to be 138.386 lb m/s. In Table 2, the intermediate results for the uncertainty calculation by the Taylor series method are shown. The values for the relative sensitivities, 0'(/), are taken from PTC 19.1. The relative sensitivity is defined by the following equation: o ' ( 0 = (dR/g)/[de(t)/e(r)]

case is equal to MDOT. The terms 0(/) x BP(I) and 0(/) x SP(I) displayed in Table 2 represent the contributions to the bias error of the result and the precision index of the result from each of the input variables. The term (0(I) x SP(1)4/VP(1) is used in the denominator of the Welch-Satterthwaite formula for determining the number of degrees of freedom of the result. Table 3 displays the values calculated by the perturbation method that correspond to the values shown in the last 3 columns of Table 2. The values for BRP(I), SRP(1), and DENV(1) that appear in Table 3 correspond to the values in Table 2 that appear in the 0(/) x BP(I), 0(I) x SP(1), and (0(/) x SP(I))4/VP(!) columns, respectively. Comparisons of the values in the two tables demonstrate that for this problem, the perturbation and the Taylor series methods yield the same values to the third significant figure. The values for the bias error and precision index of the result can be found by root-sum-squaring the values in the bias error and precision index contribution columns of both tables. The values obtained are displayed in Table 4, which also displays the degrees of freedom for the result calculated by both methods. The total uncertainty, calculated by root-sum-squaring the bias error of the result and the product of the precision index of the result and the appropriate Student t value (t = 2 for 50 degrees of freedom) is also shown in Table 4. As can be seen, the uncertainty of MDOT as calculated by the perturbation method is consistent with the uncertainty calculated by the Taylor series method to the third significant figure.

(13)

Therefore, the absolute sensitivities (0(/)) are calculated by multiplying 0'(/) times R/P(1), where R in this

CONCLUSIONS The perturbation method for propagating bias and precision errors has been shown to be quite simple to

Table 1 Uncertainty Calculation Input Set

Independent Variable

Array Element # (1)

Nominal Value (P(I))

Absolute Bias Error (BP(1))

Absolute Precision Index (SP(I))

Degrees of Freedom (VP(1)) 0

C

1

0.984

0.0075

0

Y

2

I

0

0

0

DT

3

4

0.001

0

0

DI

4

6

0.002

0

0

F

5

I

0

0

0

G

6

62.37

0.004

0.002

40

H

7

100

0.3

0.4

50

I S A Transactions •

Vol. 29, No. 4

75

Table 2 Uncertainty Calculation Intermediate Results (Taylor Series Method)

Variable Name

Relative Sensitivity (0'(1)= (dR/R)/ (dP(I)/P(1)))

Absolute Sensitivity (0(I)= dRldP(1))

0(I) x BP(I)

O(I) x SP(I)

O(I) x SPO))4/ VP(I) 0

C

1

1.406E+02

1.0553+00

0

Y

0

0

0

0

0

DT

2.4923

8.622E+01

8.622E-02

0

0

DI

0.4923

I. 135E+01

2.271E-02

0

0

F

0

0

0

0

0

G

0.5

1. 109E+00

4.438E--03

2.21E-03

6.059E-13

H

0.5

6.919E-01

2.076E-01

2.768E-01

1.174E--04

Table 3 Uncertainty Calculation Intermediate Results (Perturbation Method) Array Variable Name

Element # (I)

BRP(I)

SRP(I)

DENV(1)

C

I

1.055E+00

0

0

Y

2

0

0

0

DT

3

8.625E-02

0

0

DI

4

F

5

G H

-2.269E--02

0

0

0

0

0

6

4.438E--03

2.219E--03

6.059E-13

7

2.074E--01

2.765E--01

I.169315--04

Table 4 Uncertainty Calculation Results (Taylor Series and Perturbation Methods) Taylor

Series Method

Perturbation Method

1.079

1.079

0.2768

0.2765

Degrees of Freedom:

50.01

50.01

Absolute Uncertainty of Result (lb m/S):

1.212

1.212

0.8761

0.8759

Absolute Bias Error Ib m/S): Absolute Precision Index (Ib m/S):

RelativeUncertaintyof Result(%):

76

ISA Transactions

•

Vol. 29. No. 4

implement when compared to the Taylor series method. Its incorporation into the data analysis computer code at the fluidized bed demonstration program at Nucla, Colorado, has enabled the test team to quickly and easily determine the uncertainties in 200 dependent variables. The FORTRAN code, which propagates the uncertainties into these 200 dependent variables, consists of only 30 lines. The use of the perturbation method is highly recommended whenever measurement uncertainties are to be propagated into a result and the results calculations are being performed on a computer.

NOMENCLATURE BP~ Absolute bias error associated with the ith in-

S

r

SRP~

Contribution to the absolute precision index of the result due to SPi; see Eq. (5)

0i

Absolute sensitivity - the ratio of the change in a result to. a unit change in the ith parameter; see Eq. (4)

0i'

Relative sensitivity - the ratio of the fractional change in a result to the fractional change in the ith parameter; see Eq. (13)

Ur

Uncertainty in the result; see Eq. 11

VPi

Number of degrees of freedom for the ith input parameter

v,

Number of degrees of freedom for the result; see Eqs. (9) and (10)

put parameter B

Absolute bias error of the result; see Eqs. (3) and (8)

BRP.

Contribution to the absolute bias error of the result due to BP~; see Eq. (6)

REFERENCES 1.

J

Number of independent parameter contributing to a result

Pi

Average of the ith input parameter

R

Result of a calculation based upon one or more input parameters

R,

Result of a calculation in which one input parameter has been incremented by the amount of its bias error or precision index

SP~

Absolute precision index associated with the ith input parameter

Absolute precision index of the result

.

3.

ASME Performance Test Code Committee No. 19.1 on Measurement Uncertainty, "Measurement Uncertainty," PTC 19.1-1985, American Society of Mechanical Engineers. Moffat, R. J., "Contributions to the Theory of Single Sample Uncertainty Analysis," ASME Journal of Fluids Engineering, Vol. 104, June 1982, pp 250-260. Moffat, R. J., "Using Uncertainty Analysis in the Planning of an Experiment," ASME Journal of Fluids Engineering, Vol. 107, June 1985, pp 173-178.

ISA Transactions • Vol. 29, No. 4

77