Calculatin Uncertainty

Industrial Research Limited Report 1073

Calculating uncertainty automatically in instrumentation systems B. D. Hall

Measurement Standards Laboratory of New Zealand Lower Hutt, New Zealand 22 February 2002

Calculating uncertainty automatically in instrumentation systems B. D. Hall Industrial Research Limited Report 1073

Reference B. D. Hall; 22 February 2002. Calculating uncertainty automatically in instrumentation systems. Industrial Research Limited Report 1073; Lower Hutt, New Zealand.

Summary The report describes a simple and general algorithm for propagating uncertainty information between modules in an instrumentation system. The algorithm implements current international best-practice in the evaluation and reporting of measurement uncertainty. It allows modular systems to be designed to propagate uncertainty correctly without concern for the specific uncertainty data associated with individual modules. The approach described could lead to future systems that are: more flexible, more reliable, easier to design and easier to maintain. The technique has several important features which make it attractive for instrumentation: • It supports ‘plug-and-play’ of modular units – the ability to change one module for another without disrupting the system; • It propagates uncertainty components, as opposed to combined uncertainties – a complete set of uncertainty components is required to correctly evaluate the combined uncertainty of a measured value; • it eliminates the need to derive equations describing uncertainty for a particular system – the equations are obtained implicitly and automatically by software. The report concludes that provision of automatic uncertainty calculations in modern measurement systems is not complicated. The report contains a concise presentation of the algorithm in a mathematical form, followed by a simple application example illustrating its use. Software implementation of the technique is discussed, with pseudo C++ code annotations. A complete, proof-of-concept, C++ implementation for the example is given in an Appendix.

Contents 1

Introduction

1

2

An algorithm for propagating uncertainty

2

3

An example

4

4

Applicability and validity of the technique

6

5

A C++ implementation

6

6

Modules as intermediate steps in a calculation

8

6.1

Decomposing the power module equations . . . . . . . . . . . . . . .

8

6.2

Re-factoring ProcessorModule . . . . . . . . . . . . . . . . . . . . . . 10

6.3

Abstract modules can represent mathematical operations . . . . . . 11

7

Discussion and Conclusions

12

A

Standard practice in evaluating uncertainty

13

B

The chain rule

15

C

Software

15

C.1

Main . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

C.2

module.h . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

C.3

classes.h . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

C.4

helpers.h . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

C.5

utility functions.h . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

Legal issues The basis of this work is protected by a New Zealand Provisional Patent 1 . In addition, Industrial Research Limited retains copyright over the software presented herein.

1 Introduction This report is motivated by a problem pervading the Test & Measurement industry: there is no easy way of automatically calculating measurement uncertainty in modular measurement systems. This is an important issue for the design and realisation of reliable instrumentation systems. It is well-known that testing and validation imposes a heavy burden on performance-critical systems. Uncertainty is at the heart of the problem. Nevertheless, although many other technological issues affecting modularity in instrumentation systems are currently under industry review [1], the fundamental question of measurement uncertainty is not receiving enough attention. Modularity is an important engineering concept. A system composed of modular components is flexible because the configuration can be chosen from among sets of interchangeable parts. In addition, well-designed components generally find applications in a range of systems, resulting in economies of scale through reuse. In instrumentation systems, modularity is underpinned by microprocessor technology and standard communication interfaces between modules, which became available some thirty years ago. Since then, ‘automation’ of data acquisition and control has become ubiquitous. In essence, such automation implies that a suite of instruments, or sub-systems, is co-ordinated under computer control. This represents a modular measurement system. Measurement systems are invariably ‘value-centric’ - they return a measured value as the result of a measurement procedure. While any measurement is made with the intention of determining the value of a specific quantity (the ‘measurand’), no measurement result should be considered exact; it provides only an estimate for the measurand and so a qualifying statement of the uncertainty in this estimate is essential [2, Section 3.1]. Measurement uncertainty is important because it defines the variability of an estimate of the quantity of interest. Together, the measured value and the uncertainty determine a range of values that might reasonably be attributed to the measurand. Uncertainty is therefore intimately linked to the reliability of a value. Increasingly, reliability in measurement (more formally referred to as ‘traceability’) is demanded, because without reliable information, erroneous decisions are easily made. The high standards of reliability expected of systems in fields such as healthcare, transport and defense, etc, falls back on the accuracy of the measurements made by these systems. The value-centric design of measurement systems prevents a user from easily obtaining an objective statement of uncertainty. Measurement instruments rarely seem to report uncertainty explicitly. More often, handbooks and calibration records must be examined and interpreted in the context of a particular measurement. In complex modular systems the relationship between the uncertainty in the measurand and the uncertainty contributed by a particular module may not be obvious. 1

B. D. Hall and R. Willink, “Uncertainty Calculation System and Method”, NZ Patent Appln. No. 512212, (June, 2001)

1

The possibility of changing modules only complicates the situation. In this regard, it is important to distinguish clearly between interchangeability based on measurement function alone, as opposed to function and uncertainty. For instance, many different instruments can measure DC voltage. However, although a number may operate over a particular voltage-range, the uncertainty in their measurements is unlikely to be comparable. Clearly, accuracy in a value-centric system is difficult to predict when component instruments can be exchanged. In practice, laborious testing and re-validation may be necessary to guarantee performance. While modular measurement systems technology is continually improving, support for handling uncertainty is sorely lacking. With ‘intelligent’ modern instruments, it must be possible to do better. The uncertainty inherent in use of a particular instrument is known to its designers. If this information were built in to instruments, it could be accessed conveniently in software (essentially the equivalent of what is usually tabulated in instrument manuals) to report uncertainty, while taking into account the instrument configuration and calibration data. This report asserts that measurement systems can be designed to evaluate and propagate uncertainty automatically. It gives a simple technique that can be used to obtain a dynamic and self-consistent view of system uncertainty. If modular components are changed the resulting system immediately adapts to the change. The technique strictly adheres to the recommendations of the ‘Guide to the Expression of Uncertainty in Measurement’ (Guide) [2], published by the ISO, which is widely accepted as describing best-practice in this area (a summary is given in Appendix A). 2 This report presents the solution to a fundamental and wide-spread problem inherent in most modern instrumentation systems. It begins with a succinct presentation of an algorithm to propagate uncertainty. Then it uses a simple example to illustrate the technique. Initially, excerpts of C++ code are used to show, in principle, how software could be written to implement the technique. The report then shows that the same algorithm can be applied to the evaluation of mathematical expressions, thereby automating calculation of uncertainty in software. A complete C++ program is supplied in Appendix C showing how this technique is realised.

2 An algorithm for propagating uncertainty This section describes an algorithm for propagating values of uncertainty between modules in a measurement system. The term ‘module’ refers to a component of a measurement system. Usually a set of physically, or logically, distinct modules constitute a particular measurement system, however, the notion of a module can also be extended to abstract mathematical expressions, which can simplify software development (Section 6). Strictly, a module must have a single output but may have several inputs, or none at all. A module without inputs is referred to as a ‘leaf’ module or ‘leaf’ input. 3 For example, transducers are generally instances of leaf nodes. Modules with inputs may (as part of a more complex measurement procedure) manipulate intermediate values. In mathematical terms, a measurement can be expressed as a function of the inputs to the system xsys = fsys (x1 , · · · , xi ) , 2 3

(1)

There is an equivalent ANSI document [3] as well as guidelines prepared by NIST [4]. This alludes to the tree-like structure of interconnected modules in a system.

2

where the parameters x1 , · · · , xi can each contribute a significant component of uncertainty to the result, xsys . When a system is composed of several modules, the system measurement function can be decomposed into a set of module functions, written in the form xj = fj (Λj ) ,

(2)

where ‘j’ identifies a module, xj is the output and Λj is the set of inputs. Note that module subscripts are assigned so that j > k, where k is the subscript of any member of the set Λ j . For a system consisting of m modules, of which l are leaf modules, the system output, x m , will be obtained recursively: module m calls its input modules, which call their inputs, and so on, down to the leaf inputs (x1 , · · · , xl ). The process then reverses, as the modules called return values to their callers. The flow of data in this second phase can be regarded as an algorithm for evaluating xm . This can be expressed iteratively as for j = l + 1, · · · , m xj = fj (Λj ) .

(3)

Each step in (3) is associated with a single module and its immediate inputs, reflecting the connections in the system. Equation (3) also shows that the evaluation of a system measurement result is distributed throughout the modules of the system. We use the so-called “law of propagation of error” to evaluate the uncertainty in x m . When inputs are uncorrelated, the combined uncertainty in x m due to uncertainties in x1 , · · · , xl can be expressed as " l  2 #1/2 X ∂fm u(xi ) . (4) u(xm ) = ∂xi i=1

When inputs are correlated, the appropriate equation is      1/2 l  l X X ∂fm ∂fm u(xi ) r(xi , xj ) u(xj )  u(xm ) =  , ∂xi ∂xj

(5)

i=1 j=1

where the r(xi , xj ) are correlation coefficients and the u(xi ) are the standard uncertainties of the inputs. It is convenient to define a term for the product of a sensitivity coefficient and an input uncertainty. This will be referred to as a ‘component’ of uncertainty, so, for example, ui (xm ) ≡

∂fm u(xi ) ∂xi

(6)

is the ‘component’ of uncertainty in xm due to uncertainty in xi . Note, however, that component values take the sign of the partial derivative and are therefore not strictly-speaking uncertainties: the uncertainty associated with a measured value is always positive. The value of combined uncertainty (always positive) is obtained by adding components in quadrature if inputs are uncorrelated " l #1/2 X 2 , (7) u(xm ) = ui (xm ) i=1

or by evaluating



u(xm ) = 

l X l X i=1 j=1

1/2

ui (xm ) r(xi , xj ) uj (xm ) 

(8)

3

if there is correlation. Equations (7) and (8) are clearly equivalent to (4) and (5). Uncertainty can be propagated in a system by an algorithm similar to (3). The component ui (xm ), introduced in (6), is obtained iteratively, for any given i, by: for j = l + 1, · · · , m

X ∂fj ui (xk ) . ∂xk

ui (xj ) =

(9)

xk ∈Λj

Each step in this iteration is essentially an application of the chain rule for partial differentiation of a function (see Appendix B). The sensitivity of each module function to x i is propagated with a weighting factor of u(xi ). It is also important to note that, by definition, ui (xi ) ≡ u(xi ). Equation (9) is the main result of this report. It shows that uncertainty components can be propagated in a similar way to values. Note that the iteration order is the same as (3) and that the inputs required on the right-hand side of (9) are drawn from Λ j , the inputs to the j th module. Therefore each iteration step in (9) represents a calculation that can be performed by the module. Hence, the evaluation of system measurement uncertainty can be distributed among the modules of a system, with each module evaluating the uncertainty components pertaining to its output value. The next section will show this in the context of a simple example.

3 An example Suppose that a system is being designed to estimate the electrical power from a measurement of the potential difference across a calibrated resistor, R. If x 1 is the measured voltage then the power, y, is given by x2 y= 1. (10) R Further, the resistor is temperature dependent R = R0 (1 + α(x2 − T0 )) ,

(11)

where x2 is the measured temperature, α is the temperature coefficient and R 0 is the resistance at temperature T0 (α, R0 and T0 have negligible uncertainty here). The measurement system is realized in three modules: a voltage sensor, temperature sensor, and a processing module (e.g., a computer), with which the sensors both communicate (Fig. 1). The sensor modules are regarded as ‘black boxes’ – their inner workings are hidden. However, both provide an ‘uncertainty’ function in their communications interface, so that values for ui (x1 ) and ui (x2 ) can be obtained by the processor module. It is assumed that the sensor measurements are independent, i.e. that u2 (x1 ) = u1 (x2 ) = 0. The discussion focuses on the processor module design. There are three physical modules in the system, so y = x3 and, from (10) and (11), y = f3 (x1 , x2 ) =

x21 . [R0 (1 + α(x2 − T0 ))]

(12)

The uncertainty component calculation is represented by the i = 3 iteration of (9), i.e.: ui (x3 ) =

X ∂f3 ui (xk ) . ∂xk

(13)

xk ∈Λ3

4

3

Processor

1

Voltage

2

Temperature

Figure 1: Two sensors are connected to a processor module to realize the power measurement system. The small empty circles represent a common communications interface, the larger circles identify the modules by number. Modules 1 and 2 are ‘leaf’ modules.

Now, Λ3 = {x1 , x2 } so (13) can be written as4 ui (x3 ) =

∂f3 ∂f3 ui (x1 ) + ui (x2 ) , ∂x1 ∂x2

(14)

and the partial derivatives are, from (12), ∂f3 ∂x1 ∂f3 ∂x2

= =

2x1 [R0 (1 + α(x2 − T0 ))] −R0 α x21 . [R0 (1 + α(x2 − T0 ))]2

(15) (16)

Equations (14) through (16) specify the requirements for an uncertainty interface function for module 3. The function should take a parameter representing the value of i, and return the desired uncertainty component. In this example, it is tempting to exploit the independence of the two leaf inputs, so the only contribution to the sum when i = 1 is xk = x1 and when i = 2 is xk = x2 . However, this has been avoided in describing the solution above to make the point that it would limit the ‘modularity’ of the processor module. To see why, suppose that x 1 and x2 are later exchanged for other devices that do have explicit dependencies on other quantities. For instance, the new sensor modules may be sensitive to ambient temperature, which would introduce another leaf input to the system (ambient temperature) and alter the connection topology of the modules (ambient temperature will be associated with values of voltage and temperature). An implementation of the more general form of (14) through (16) allows, in principle, the processor module to function correctly in spite of the substantial change to the overall uncertainty calculation. 4 By definition u1 (x1 ) and u2 (x2 ) are synonymous with u(x1 ) and u(x2 ) respectively. It is instructive to leave the subscripts on u because they correspond to a parameter passed to the uncertainty software function.

5

4 Applicability and validity of the technique The evaluation of uncertainty according to equation (9) is made possible by nature of the equations recommended in the Guide and is ideally suited for complex modular systems. The reader should note, however, that there are assumptions underlying the Guide’s approach to propagating and interpreting measurement uncertainties. First, it is assumed that the measurement function can be approximated by a Taylor series truncated beyond linear terms: in other words, that the measurement function, in the vicinity of the measurement point, can be considered linear on the scale of variations associated with the uncertainties. One would expect this assumption to be satisfied in the majority of cases. Nevertheless, it is conceivable that system software be designed to test this assumption in-situ. Second, the Guide assumes that the uncertainty associated with a measurand can be approximated by a Gaussian distribution (or Student’s t).5 If these approximations do not hold then the proposed method may not apply. The technique could be applied to complex systems that use instrument technology lacking a built-in uncertainty function. In this case, the systems designer will need to provide software to encapsulate each instrument, so that it presents a suitable module interface to the system. In other words, the system designer must write the uncertainty function, based on information provided by the manufacturer. Such designer-provided software is similar to the concept of a ‘Role Control Module’ (RCM) for instruments, as discussed in [1]. A RCM defines a generic instrument interface with the particular functions required by a system, and features specific to an instrument are not directly accessible. In this way an instrument’s role can be quite narrowly defined within the context of a system and decoupled from a particular instrument’s capabilities. The task of writing an uncertainty function therefore applies only to the instrument functionality required by the system. For a measurement result to be considered reliable, the measurement procedure should be demonstrably ‘traceable’ to primary standards and the uncertainties quantified [5]. Traceable measurements are increasingly demanded, because they carry an assurance of quality and reliability. The analysis of uncertainty contributions in a particular measurement procedure is a task requiring specialist skills. Various considerations apply when assessing uncertainties in measurement, which are covered in the Guide [2] (and [3, 4]). This paper simply assumes that systems and instrument designers will have this competency. It should be clear that by adhering to common guidelines and dealing with uncertainty calculations at module-level, rather than at system-level, the uncertainty calculations will be simplified. It will also be easier to test and validate systems, because modules can be treated independently. So, in general terms, the techniques described here will help to design and maintain complex measurement systems.

5 A C++ implementation In general, software implementation will depend on the choice of technology and on the performance and functionality required. This section presents an outline of the technique in C++. 6 To create modules that propagate uncertainty a standard interface with two function calls is used: one function returns a value, the other an uncertainty component. The value function 5 6

The Guide does allow input distributions to be of different types, e.g. uniform or triangular. More examples are available, in C++ and Basic [6, 7].

6

takes no parameters, while the uncertainty function must identify the input associated with the uncertainty component of interest. The declaration of a pure abstract C++ class for this interface looks like class Interface { public: virtual double fn() = 0; virtual double uComponent(Interface* id) = 0; };

The processor module must implement this interface. Assuming that the sensors can be accessed through a pointer to Interface, a processor module class would look something like class ProcessorModule : public Interface { private: Interface* x1; Interface* x2; double R0,T0,alpha; public: ProcessorModule(Interface* v,Interface* t) : x1(v), x2(t) {} virtual double fn() { return x1.fn() * x1.fn() / resistance( x2.fn() ); } virtual double uComponent(Interface* id) { double t1 = 2.0 * x1.fn() / resistance( x2.fn() ); t1 *= x1.uComponent(id); double t2 = - x1.fn() * x1.fn() * R0 * alpha / resistance( x2.fn() ); t2 *= x2.uComponent(id); return t1 + t2; } } private: double resistance(double T){ return R0 * (1 + alpha * (T - T0)); } };

Note that fn() evaluates equation (12), and uComponent() evaluates equation (14) with (15) and (16). In general, when a module calculates uncertainty it performs a single iteration of (9). Calls are made to the module’s inputs to obtain values for ui (xk ). One parameter is required in each uncertainty function call, corresponding to i – here a pointer to Interface is used. By using Interface, the processor module is made flexible. Either the voltage or temperature sensors could be replaced without requiring changes to the definition. ProcessorModule always obtains information derived from connected sensors, so when a change occurs the appropriate uncertainty or value is used. In addition, because the ProcessorModule implements Interface, it can be used as a component in another system. The implementation can nevertheless be improved upon. In particular, it is possible to automate the uncertainty calculations, so that there is no longer a need to derive equations (15) and (16). The next two sections develop this approach.

7

6 Modules as intermediate steps in a calculation The development of software can be simplified by extending the concept of a module to abstract elements of an instrument measurement function. Indeed, the need for an explicit expression of the module function derivatives can be eliminated entirely; in the example there is no need to derive equations (15) and (16). This is achieved by assembling a set of software objects to represent the module function. The resulting structure can evaluate both function values and uncertainties.

6.1 Decomposing the power module equations In this section, the equations for power and its uncertainty will be decomposed into a number of abstract modules which can be represented as software objects. This software will be outlined in section 6.2, and developed further in Appendix C. The reader should note that more modules are used to analyse the example, so y no longer refers to x3 , and x3 no longer refers to power. If a simple pocket-calculator is used to evaluate equation (12), a number of intermediate results are obtained before the value for y. Leaving aside steps associated with the constants 7 R0 , α, T0 and the initial values for x2 and x1 , there are two intermediate steps required, associated here with x3 and x4 : x3 = x21 , x4 = R0 [1 + α (x2 − T0 )] , x3 y= . x4

(17)

These steps are represented in Fig. 2. Nodes terminating a branch in this diagram are either input quantities or the final result; those within a branch represent intermediate steps. The evaluation of this set of equations for y (now equivalent to x 5 ) is of the form of (3). In fact, the nodes can be considered modules. y = x3 / x 4

x4 = R0[1 + α (x2 - T0) ]

x3 = ( x1 )2

x1

x2

Figure 2: The schematic representation of the decomposition of the power measurement model equation (12) into equations (17). R0 , α and T0 are considered exact.

Corresponding sets of equations can be found for the uncertainty components associated with equation set (17) (formally, these are obtained by applying equation (21) from Appendix A). 7

Skipped here because uncertainties are not involved.

8

For the component u1 (y) the set is u1 (x3 ) = 2 x1 u(x1 ) , u1 (x4 ) = 0 , u1 (x3 ) ; u1 (y) = x4

(18)

and for u2 (y) u2 (x3 ) = 0 , u2 (x4 ) = R0 α u(x2 ) , x3 u2 (y) = − 2 u2 (x4 ) . x4

(19)

These are clearly the equivalent of steps in (9). The topology associated with these equations, shown in Figs. 3 and 4, is identical to Fig. 2, although nodes can now represent values or uncertainties, or both, according to the form of the equations. u1(y) = u1( x3 ) / x4 x4 = R0[1 + α (x2 - T0) ]

u1( x3 ) = 2 x1 u( x1 )

x1

x2

u( x1 ) Figure 3: The uncertainty calculation steps of equations (18).

u2(y) = -u2( x4 ) x3 / ( x4 )2

x4 = R0[1 + α (x2 - T0) ] x3 = ( x 1 ) 2

u2( x4 ) = R0 α u( x2 )

x1

x2 u( x2 )

Figure 4: The uncertainty calculation steps of equations (19).

The figures emphasize that the evaluation of a measurement value or an uncertainty is hierarchical: each step depends only on previous results. Hence, as the evaluation proceeds, uncertainty and value can be calculated together. Furthermore, expressions for the intermediate uncertainty steps are derived from the first partial derivatives of the expressions for the measurement value evaluation steps.

9

6.2 Re-factoring ProcessorModule Simple mathematical expressions can be evaluated efficiently by decomposing them into tree-like structures [8]8 and applying a top-down recursive algorithm like (3). This section shows how a tree of software objects, representing mathematical operations and behaving as modules, can be generated automatically, thereby simplifying the job of writing code for a particular measurement system. While the ProcessorModule class represents the physical modularity of the system, it does not decompose the problem to the degree shown in figures 2-4. The design can be improved by re-factoring ProcessorModule using the following three classes: Square, Resistance and Divide: Square describes an object that simply squares its input (module). class Square : public Interface { private: Interface* x1; public: Square(Interface* x) : x1(x) {} virtual double fn() { return x1->fn() * x1->fn(); } virtual double uComponent(Interface* id) { return 2.0 * x1->fn() * x1->uComponent(id); } };

Resistance describes the temperature dependent resistance and its uncertainty. class Resistance : public Interface { private: double R0,T0,alpha; Interface* x2; public: Resistance(Interface* x) : x2(x) {} virtual double fn(){ return R0 * ( 1 + alpha * ( x2->fn() - T0 ) ); } virtual double uComponent(Interface* id) { return R0 * alpha * x2->uComponent(id); } };

Divide describes a class of object that takes the ratio of its two inputs. class Divide : public Interface { private: Interface* x3; Interface* x4; public: Divide(Interface* n,Interface* d) : x3(n), x4(d) {} virtual double fn() { return x3->fn() / x4->fn(); } virtual double uComponent(Interface* id) { double den = x4->fn() * x4->fn(); double num = x4->fn() * x3->uComponent(id) - x3->fn() * x4->uComponent(id); return num / den; } 8

Strictly, they are directed acyclic graphs because nodes can link to several parts of the structure. The term ‘parse-tree’ is also used.

10

};

In terms of these three classes, a rather more elegant class can be defined for the processor module, along the lines of: class FactoredProcessorModule : public Interface { private: Interface* power; double R0,T0,alpha; public: ProcessorModule(Interface* v,Interface* t) { Interface* R( new Resistance(t) ); Interface* v2( new Square(v) ); power = new Divide(v2,R); } virtual double fn() { return power.fn(); } virtual double uComponent(Interface* id) { return power.uComponent(id); } };

When constructed, objects of this class create three objects: one for the resistor, one for the operation of squaring the voltage and one for taking the ratio of voltage squared to resistance. They are linked together to form a network of modules that can evaluate the power and its uncertainty. Note that calls to the functions fn() and uComponent() are now simply forwarded to the pointer, power.

6.3 Abstract modules can represent mathematical operations The two classes Square and Divide are generic. They perform elementary mathematical operations, and calculate uncertainty, on their inputs. A library of such classes, representing all the arithmetic operations and a set of standard mathematical functions (i.e., add, subtract, multiply, sin, log, etc), will allow most measurement functions to be coded directly (i.e., no expressions for uncertainty components are required). This approach is very similar to a technique known as ‘Automatic Differentiation’ (AD), used in numerical analysis, which can automatically evaluate the derivatives of an arbitrary expression [9]. AD is distinct from both symbolic differentiation and numerical finite-difference methods; it is exact (to within numerical round-off) and efficient. For measurement systems, an adaptation of the AD technique to propagate uncertainty is straightforward. In Appendix C a small library of elementary classes is given and applied to the example. Using this class library, the Resistance class can be re-written as the following Listing 1: main.cpp class adResistance : public Interface { private: ModulePtr rValue; public: adResistance( const ModulePtr& x2 , double R0 , double T0 , double alpha ) { ➥

11

Listing 1: main.cpp (continued) rValue = R0 * ( 1.0 + alpha * (x2 - T0 ) ); } virtual double fn() const { return rValue->fn(); } virtual double uComponent (const Interface* id) const { return rValue->uComponent(id); } virtual Dependents dependsOn() const { return rValue->dependsOn(); } };

Note, how an expression for resistance is assembled during the construction of an adResistance object, by calls to overloaded arithmetic operators. These are ‘helper functions’ associated with particular classes implementing arithmetic operations. For example, the helper functions to create an instance of Addition, an object that adds its inputs are: Listing 2: helpers.h static inline ModulePtr operator+(const ModulePtr& l,const ModulePtr& r) { return new Addition(l,r); } static inline ModulePtr operator+(double l,const ModulePtr& r) { ModulePtr tmp = leafConstant(l); return new Addition(tmp,r); } static inline ModulePtr operator+(const ModulePtr& l,double r) { ModulePtr tmp = leafConstant(r); return new Addition(l,tmp); }

There are three variants to allow for the two combinations of a pure number and an object (i.e. module) and a pair of objects. The class Addition is simply: Listing 3: classes.h class Addition : public Interface { private: ModulePtr lhs; ModulePtr rhs; public: Addition(const ModulePtr& l,const ModulePtr& r) : lhs(l), rhs(r) {} virtual double fn() const { return lhs->fn() + rhs->fn(); } virtual double uComponent(const Interface* id) const { return lhs->uComponent(id) + rhs->uComponent(id); } virtual Dependents dependsOn() const { return bivariateDependents(lhs,rhs); } };

7 Discussion and Conclusions This report has shown that a simple algorithm can enable modular measurement systems to propagate uncertainty information supplied by their component parts. The evaluation of

12

uncertainty according to equation (9) is distributed among the modules of a system, just as calculations of value are normally handled (i.e. equation (3)). The technique provides a system framework for combining and propagating uncertainty information that is independent of the specifics of module uncertainty data. Each module will evaluate the uncertainty of its output, which may depend on module-inputs. In short, the technique has re-factored the uncertainty problem for modular systems and isolated the system design from module properties. This then provides for independent testing and validation of modules and systems. The technique adheres to current best-practice for evaluating uncertainty. It does so with a minimum of complexity, both to the external specification of software interfaces and with regard to internal implementation details. Emerging standards for future generations of measurement systems and sensors have acknowledged the importance of uncertainty in measurement but do not provide ways to automate it [10]. Our method provides a general solution to this problem. The technique has several important features which make it attractive for instrumentation. First, it supports ‘plug-n-play’ of modular units – the ability to change one module for another without disrupting the system. Thus, in the example, one can change the thermometer, or the voltmeter, without interfering with the processor module – if a replacement sensor has different uncertainty characteristics, these are automatically propagated through the system. The plug-and-play property will provide huge benefits for systems that must operate under strict performance requirements. In these systems, if a part is changed there is often a need to re-validate the entire system’s performance – a process that is both time-consuming and costly. With our technique, there is a simple way for modules to transmit a value with an associated uncertainty. Thus the evaluation of uncertainty is localised on the scale of individual modules, which makes it possible to test and validate their performance individually. When a system is designed to use our technique, the impact on measurement accuracy of any change is always explicit. It will be much easier to validate such systems and they will be safer to operate. Secondly, our technique has the useful property of propagating uncertainty components, as opposed to combined uncertainties. Uncertainty components reflect the different influence factors in a final combined uncertainty: a complete set of components is needed to describe the combined uncertainty of the measurand. Having access to this information identifies dominant uncertainties and may be useful where the dominant term changes according to operating conditions. It also allows correlations between input quantities to be handled correctly. A third important feature of the technique is that the derivation of the equations for uncertainty can be automated. This simplifies the development of measurement software as well as making it more robust, because there is no need to separately analyze the measurement function. Furthermore, a software library required to implement this feature can be independently tested and validated and then reused many times. In conclusion, the provision of automatic uncertainty calculations in modern measurement systems is not complicated. The techniques described in this report will lead to systems that are: more flexible, more reliable, easier to design and easier to maintain.

A Standard practice in evaluating uncertainty Guidelines for the evaluation of measurement uncertainties have been published in the Guide to the Expression of Uncertainty in Measurement [2]. The recommendations of Guide are accepted as current best-practice and have been adopted by national measurement institutes and accredited measurement laboratories world-wide.

13

The measurement of a quantity Xm , the measurand9 , is described by a function Xm = f (X1 , X2 , . . . , Xl ), interpreted as “...that function which contains every quantity, including all corrections and correction factors, that can contribute a significant component of uncertainty to the measurement result”[2, 4.1.2]. The quantities X 1 , X2 , . . . , Xl and Xm are random variables, however, the parameters describing their distribution (mean and standard deviation) are not known exactly. Estimates of the means, x 1 , x2 , . . . , xl , must be used to estimate the measurand as xm = f (x1 , x2 , . . . , xl ) . (20) A ‘standard uncertainty’, u(xi ), is associated with each xi and is understood to be an estimate of the standard deviation of Xi . The quality of this estimate can be characterized by a number of ‘degrees-of-freedom’, νi . If νi is infinite the estimate is considered exact, otherwise νi is related to the relative uncertainty of u(xi ) [2, G.4.2]. Degrees of freedom has its conventional meaning when associated with a normally distributed random variable. The uncertainty in xm depends on the uncertainties in x1 , x2 , . . . , xl and on the form of f . This report use the notation10 ∂f u(xi ) , (21) ui (xm ) ≡ ∂xi where the partial derivative is more formally expressed as ∂f ∂f . (22) ≡ ∂xi ∂Xi x1 ,x2 ,...,xl

The combined standard uncertainty in xm is the root-sum-square of these components11 " l #1/2 X 2 uc (xm ) = ui (xm ) i=1

=

"

l  X ∂f i=1

∂xi

u(xi )

2 #1/2

,

(23)

provided input quantities are uncorrelated. With correlation the combined standard uncertainty should be evaluated as  1/2 l X l X uc (xm ) =  ui (xm ) r(xi , xj ) uj (xm ) i=1 j=1

1/2 l X l X ∂f ∂f u(xi ) r(xi , xj ) u(xj ) =  ∂xi ∂xj 

(24)

i=1 j=1

where r(xi , xj ) ≡ u(xi , xj )/(u(xi ) u(xj )) is the correlation coefficient and u(xi , xj ) is the estimated covariance of Xi and Xj . If one or more of the input-quantity uncertainties has finite degrees of freedom, a value for the ‘effective degrees-of-freedom’, νeff , of u(y) should be calculated using the Welch-Satterthwaite formula [2, G.4] l

u4c (xm ) X u4i (y) = . νeff νi

(25)

i=1

9

In the Guide, ‘Y ’ denotes the measurand (and y as its estimate). The Guide has a similar notation for uncertainty components in the measurand but defines them as the modulus of the ui (xm ) used in this report. We prefer to retain the sign because the implementation of our algorithm is simplified. Note that, by definition, ui (xi ) ≡ u(xi ). 11 The notation for combined uncertainty used by the Guide is uc (xi ). The subscript ‘c’ is omitted (i.e. u(xi )) when the meaning is clear from the context, for example with input quantities. 10

14

The Welch-Satterthwaite formula does not apply when correlations are present – the Guide makes no specific recommendations in this case.

B The chain rule The steps of equation (9) are closely related to the chain rule for partial differentiation. This can be seen more clearly by writing the final iteration (j = m) of (9) in the equivalent form ui (xm ) =

X ∂xm ∂xk ∂xm u(xi ) = u(xi ) . ∂xi ∂xk ∂xi

(26)

xk ∈Λm

The common factor u(xi ) can be canceled, leaving X ∂xm ∂xk ∂xm = . ∂xi ∂xk ∂xi

(27)

xk ∈Λm

k Because, in general, the xk terms in this equation may not explicitly depend on xi , the ∂x ∂xi terms may be expanded further by applying the chain rule. That is the reason for iterating in equation (9); each step provides intermediate results needed in subsequent iterations.

C Software This appendix consists of a C++ implementation of the automated uncertainty technique, that was outlined in section 6, applied to the example. The listings in this section have been integrated with the LATEX source using a literate programming tool called ProgDOC [11]. The code represents a complete working C++ example. The Microsoft Visual C++ 6.0 compiler was used together with STLPort release 4.5 [12] and the shared pointer library from Boost 1 23 0 [13]. Aside from showing how the technique can be implemented, this appendix extends the technique presented so far by implementing the calculation of ‘combined standard uncertainty’ (23). A similar approach could be applied to evaluate the degrees-of-freedom (25). To evaluate (23)-(25), these calculations have to iterate over the set of leaf-inputs to a system. This requires a ‘system-wide’ view, rather than the single-module one. Nevertheless, implementation is not difficult and does not compromise the flexibility and modularity claimed. Code to do this is given in section C.5. It relies on an extension to the definition of the Interface class, in C.2

C.1 Main The main function takes the form Listing 4: main.cpp int main() { // Simulate sensors with leaf modules ModulePtr x1 = leafVariable(5,0.5); ModulePtr x2 = leafVariable(25,1.0);

// 5 +/- 0.5 V // 25 +/- 1 degrees C ➥

15

Listing 4: main.cpp (continued)

// R0=1000 ohms, T0=20 degrees C, alpha=0.01 ModulePtr R(new adResistance(x2,1000,20,0.01) ); // Power equation ModulePtr power = x1 * x1 / R; cout

<< "Power value = " << value(power) << " +/- " << uCombined(power) << endl << endl; cout << "Power uncertainty wrt voltage=" << uComponent(power,x1) << endl; cout << "Power uncertainty wrt temperature=" << uComponent(power,x2) << endl; return 0;

}

A separate class for the processor module has not been defined here, the main program instead treats the module power as a measurement result and reports its value and uncertainty. The sensor leaf inputs are simulated by creating Leaf objects, x1 and x2, using the helper function leafModule(). The temperature dependence of the resistance has been encapsulated by the adResistance class, which was presented in Listing 1 on page 11. Here, an instance of that class is referred to by R. This enables the power module equation to be expressed as in (10).

C.2 module.h This file specifies the module interface functions. Note that Interface now has a third function, dependsOn(). This is needed to perform the calculation of combined uncertainty (see section C.5). The ModulePtr type is a smart pointer that uses reference counting to determine when its pointer is no longer required. This eases memory management (i.e. when to delete modules that have been created by a new command). Another useful aspect of ModulePtr is that its use allows us to override the arithmetic operators (pointer arguments are not sufficient to resolve alternative operator definitions). Listing 5: module.h #ifndef module h #define module h // Allow definition of the conversion operator from shared ptr to T*. // This is not recommended by Boost! #define BOOST SMART PTR CONVERSION #include "boost/smart ptr.hpp" #include <set> namespace UN { // BEGIN Module typedef std::set Dependents; struct Interface { virtual double fn() const = 0; ➥

16

Listing 5: module.h (continued) virtual double uComponent(const Interface*) const = 0; virtual Dependents dependsOn() const = 0; }; typedef boost::shared ptr ModulePtr; // END Module } #endif

C.3 classes.h The classes.h file defines the types of modules that can be created. A more complete implementation would add mathematical functions (e.g. sine, cosine, etc). There are really two types of leaf module, one representing a known constant (with no uncertainty) the other representing an uncertain value. The Leaf class is used for both. It is possible to evaluate pure derivatives by Automatic Differentiation by setting uncertainty value of a leaf to unity (the default value). Listing 6: classes.h #ifndef classes h #define classes h #include #include "module.h" namespace UN { // BEGIN Leaf class Leaf : public Interface { private: double value ; double uncert ; public: Leaf(double v,double u=1.0) : value (v), uncert (u) {} virtual double fn() const { return value ; } virtual double uComponent(const Interface* id) const { return ( id == this ) ? uncert : 0.0 ; } virtual Dependents dependsOn() const { Dependents tmp; if(uncert != 0.0) tmp.insert(this); return tmp; } }; // END Leaf // BEGIN bivariantDependents static inline Dependents bivariateDependents(const ModulePtr& lhs,const ModulePtr& rhs) { Dependents ldep( lhs->dependsOn() ); Dependents rdep = rhs->dependsOn(); ldep.insert(rdep.begin(),rdep.end()); ➥

17

Listing 6: classes.h (continued)

return ldep; } // END bivariantDependents // BEGIN Addition class Addition : public Interface { private: ModulePtr lhs; ModulePtr rhs; public: Addition(const ModulePtr& l,const ModulePtr& r) : lhs(l), rhs(r) {} virtual double fn() const { return lhs->fn() + rhs->fn(); } virtual double uComponent(const Interface* id) const { return lhs->uComponent(id) + rhs->uComponent(id); } virtual Dependents dependsOn() const { return bivariateDependents(lhs,rhs); } }; // END Addition // BEGIN Subtraction class Subtraction : public Interface { private: ModulePtr lhs; ModulePtr rhs; public: Subtraction(const ModulePtr& l,const ModulePtr& r) : lhs(l), rhs(r) {} virtual double fn() const { return lhs->fn() - rhs->fn(); } virtual double uComponent(const Interface* id) const { return lhs->uComponent(id) - rhs->uComponent(id); } virtual Dependents dependsOn() const { return bivariateDependents(lhs,rhs); } }; // END Subtraction // BEGIN Multiplication class Multiplication : public Interface { private: ModulePtr lhs; ModulePtr rhs; public: Multiplication(const ModulePtr& l,const ModulePtr& r) : lhs(l), rhs(r) {} virtual double fn() const { return lhs->fn() * rhs->fn(); } virtual double uComponent(const Interface* id) const { return lhs->uComponent(id) * rhs->fn() + lhs->fn() * rhs->uComponent(id); } virtual Dependents dependsOn() const { return bivariateDependents(lhs,rhs); } }; // END Multiplication // BEGIN Division class Division : public Interface { private: ModulePtr lhs; ModulePtr rhs; public: ➥

18

Listing 6: classes.h (continued) Division(const ModulePtr& l,const ModulePtr& r) : lhs(l), rhs(r) {} virtual double fn() const { const double den = rhs->fn(); assert(den != 0.0); return lhs->fn() / den; } virtual double uComponent(const Interface* id) const { const double den = rhs->fn(); assert(den != 0.0); const double num = lhs->uComponent(id) * rhs->fn() lhs->fn() * rhs->uComponent(id); return num / den; } virtual Dependents dependsOn() const { return bivariateDependents(lhs,rhs); } }; // END Division } #endif

C.4 helpers.h This file defines the overloaded operator functions that enable mathematical expressions to be written with module arguments. There are also functions to help create leaf modules. Listing 7: helpers.h #ifndef helpers #define helpers #include "classes.h" namespace UN { // BEGIN leaf helpers static inline ModulePtr leafVariable(double v,double u = 1.0) { return new Leaf(v,u); } static inline ModulePtr leafConstant(double v) { return new Leaf(v,0.0); } // END leaf helpers // BEGIN Addition helpers static inline ModulePtr operator+(const ModulePtr& l,const ModulePtr& r) { return new Addition(l,r); } ➥

19

Listing 7: helpers.h (continued) static inline ModulePtr operator+(double l,const ModulePtr& r) { ModulePtr tmp = leafConstant(l); return new Addition(tmp,r); } static inline ModulePtr operator+(const ModulePtr& l,double r) { ModulePtr tmp = leafConstant(r); return new Addition(l,tmp); } // END Addition helpers // BEGIN Subtraction helpers static inline ModulePtr operator-(const ModulePtr& l,const ModulePtr& r) { return new Subtraction(l,r); } static inline ModulePtr operator-(double l,const ModulePtr& r) { ModulePtr tmp = leafConstant(l); return new Subtraction(tmp,r); } static inline ModulePtr operator-(const ModulePtr& l,double r) { ModulePtr tmp = leafConstant(r); return new Subtraction(l,tmp); } // END Subtraction helpers // BEGIN Multiplication helpers static inline ModulePtr operator*(const ModulePtr& l,const ModulePtr& r) { return new Multiplication(l,r); } static inline ModulePtr operator*(double l,const ModulePtr& r) { ModulePtr tmp = leafConstant(l); return new Multiplication(tmp,r); } static inline ModulePtr operator*(const ModulePtr& l,double r) { ModulePtr tmp = leafConstant(r); return new Multiplication(l,tmp); } // END Multiplication helpers // BEGIN Division helpers static inline ModulePtr operator/(const ModulePtr& l,const ModulePtr& r) { return new Division(l,r); } static inline ModulePtr operator/(double l,const ModulePtr& r) { ModulePtr tmp = leafConstant(l); return new Division(tmp,r); } static inline ModulePtr operator/(const ModulePtr& l,double r) { ➥

20

Listing 7: helpers.h (continued) ModulePtr tmp = leafConstant(r); return new Division(l,tmp); } // END Division helpers } #endif

C.5 utility functions.h This file defines functions that return various properties of a measurement result, namely: value, combined uncertainty and uncertainty components. The implementation of combined uncertainty relies on the dependsOn() member of Interface. This provides a set of leaf inputs to the system. Note that uCombined implements (23), which applies only when inputs are uncorrelated. Listing 8: utility functions.h #ifndef utility fns h #define utility fns h #include #include "module.h" namespace UN { // BEGIN uCombined static inline double uCombined(ModulePtr m) { Dependents dep( m->dependsOn() ); const Dependents::iterator last = dep.end(); Dependents::iterator it; double uc = 0.0; for(it = dep.begin(); it != last; ++it) { uc += m->uComponent(*it) * m->uComponent(*it); } return sqrt( uc ); } // END uCombined static inline double value(ModulePtr m) { return m->fn(); } static inline double uComponent(ModulePtr m,ModulePtr n) { return m->uComponent(n); } } #endif

21

Acknowledgement The author is very grateful to R. Willink for many very interesting and fruitful discussions regarding this work and for suggestions that have improved this manuscript.

References [1] Mueller J and Oblad R 2000 Architecture drives test system standards IEEE Spectrum September 68 [2] ISO, Guide to the expression of uncertainty in measurement, International Organisation for Standardization, Geneva, 2nd edition, 1995. [3] ANSI/NCSL Z540-2-1997, American National Standard for Expressing Uncertainty – U.S. Guide to the Expression of Uncertainty in Measurement, American National Standards Institute, 1997. [4] B. N. Taylor and C. E. Kuyatt, “Guidelines for evaluating and expressing the uncertainty of nist measurement results,” Technical Note 1297, Nat. Inst. Stand. Tech., 1994. [5] J. V. Nicholas and D. R. White “Traceable temperatures” 2nd ed. (John Wiley & Sons, Chichester, 1995) [6] Hall B D 2001 Calculating measurement uncertainty with reverse automatic differentiation Industrial Research Report 1041 [7] Perera S, Hall B D 2002 Measurement uncertainty software tools for PCs Industrial Research Report 1072 [8] R. Sedgewick, Algorithms in C++, Parts 1-4: Fundamentals, Data Structure, Sorting, Searching, Addison-Wesley, 1999. [9] L. B. Rall and G. F. Corliss, “An introduction to automatic differentiation,” in Computational Differentiation: Techniques Applications, and Tools, M. Berz, C. H. Bischof, G. F. Corliss, and A. Griewankpp, Eds., pp. 1–17. SIAM, Philadelphia, September 1996. [10] IEEE Std 1451.1-1999 ‘IEEE Standard for a Smart Transducer Interface for Sensors and Actuators – Network Capable Application Processor Information model’ (IEEE, New York, 2000); K. L. Lee and R. D. Schneeman, “Internet-based distributed measurement and control applications,” IEEE Instrum. Meas. Mag., pp. 23–27, June 1999. [11] http://www.progdoc.org/ [12] http://www.stlport.org/ [13] http://www.boost.org/

22