Evaluating Measurement Uncertainty Using Models with Arguments Subject to a Constraint

1. Introduction

The metering of natural gas serves a variety of purposes, including billing, custody transfer, assessing conformity with contractual and legislative requirements [1,2,3]. For such assessments, it is paramount to know the uncertainty of the measurement results [4]. To obtain the measurement uncertainty for the total mass, volume or energy, it is necessary to duly propagate the measurement uncertainty [5,6,7] of, e.g., the flow rate [8,9,10,11,12,13], composition [14,15], density [16], calorific value [17,18]. It has been argued that also the dependencies between the input variables of the measurement models should be taken into account [19], which is not always done.

When calculating gas properties calculated from a composition, pressure and temperature, such as the compressibility factor, density or speed of sound, it is important to consider the covariances between the amount fractions [20]. The amount fractions of all components in a gas mixture form a composition [21] and a feature of a composition is that its amount fractions sum exactly to unity. This property of a composition has implications for the covariance matrix associated with the amount fractions and also affects the differentiation of measurement models which have a composition as one of their inputs.

The property of a composition is that its amount fractions sum to unity exactly is an example of a mathematical constraint [22,23] and it is well known from mathematics that a function with constrained arguments behaves differently from the same function while ignoring the constraint. The latter case, thus ignoring the constraint can lead to physically infeasible output. For instance, when calculating a calorific value of natural gas with a composition that adds to, say, 101 cmol mol⁻¹ has no physical meaning, yet the model of ISO 6976 [18] permits at least mathematically to do so. When using an equation of state such as GERG-2008 [24,25] or AGA8-DC92 [26,27], it is less obvious whether it is mathematically possible to evaluate these functions while ignoring the constraint. Most valid implementations of these models refuse the input when the constraint on the composition is not met.

In the industry and elsewhere, the law of propagation of uncertainty (LPU) [5,7] is the de facto standard for propagating measurement uncertainty. The LPU is part of the methods of the Guide to the Expression of Uncertainty in Measurement (GUM) [4,5,6,28,29] for evaluating and propagating uncertainty in measurement. The LPU is widely used in standards [8,9,10,11,13,14,15,18], regulatory documents [2,30] and industry guidance documents [31,32]. Measurement uncertainty is propagated using a measurement model [29].

An equation of state, when used for instance to calculate a compressibility factor under actual conditions, can be used as a measurement model. When using these equations, there are in principle two kinds of uncertainty to be considered, the model uncertainty [29], clause 11.10] and the uncertainty due to the uncertainty of the input variables, e.g., the line pressure, temperature and gas composition. For the compressibility factor computed from the GERG-2008 equation of state, the model uncertainty is given in ISO 20765-2 [25], but a method to propagate the uncertainty from the input variables using the LPU is missing [25], clause 7.6].

For applying the LPU, values for the sensitivity coefficients are necessary. These are the first-order partial derivatives of the measurement model with respect to the input variables. ISO/IEC Guide 98-3 [5], clause 5.1.3] defines sensitivity coefficients as partial derivatives that “[...] describe how the output estimate y varies with changes in the values of the input estimates $x_1,x_2,\dots ,x_N$”. The GUM [5,6,7,29] makes silently the assumption that the input variables are independent. The guidance on obtaining the sensitivity coefficients [5], clause 5.1.3] instructs to obtain these coefficients by differentiating, one variable at a time, the measurement model. Varying one variable at a time the input of a measurement model with a composition as input not possible without running the risk of infeasible output: to satisfy the constraint, if one amount fraction changes, then at least one other should change too to ensure that the sum of the amount fractions is still one. There is however a method for obtaining partial derivatives that appreciate an equality constraint [23] and this method can in principle be used with equations of state as well.

Differentiating equations of state is cumbersome when doing that analytically. Many of these are complex, and furthermore many of the calculations are iterative [33,34,35]. For instance, calculating a compressibility factor using a cubic equation of state requires solving a cubic polynomial [33], which makes the function to be differentiated implicit. A vapour-liquid equilibrium calculation with the same simple equation of state is iterative [33]. Numerically evaluating the partial derivatives only requires a valid implementation of the calculation involving the equation of state and a method to obtain the from varying the input variables the values for the sensitivity coefficients. In this work, a numerical method for using the LPU is presented for calculations involving an equation of state. The methodology does not depend on the kind of model and is capable of dealing with models with a combination of constrained and unconstrained input variables.

The LPU requires the first-order partial derivatives of the measurement model with respect to the input variables for propagating the standard uncertainties from the input variables to the output variable(s). The GUM [5,6,7,29] makes silently the assumption that the input variables are independent. The guidance on obtaining the sensitivity coefficients [5], clause 5.1.3] instructs to obtain these coefficients by differentiating, one variable at a time, the measurement model. This method does not work if the input variables are subject to a constraint. In that situation, it is no longer possible to vary one by one the input variables to evaluate the sensitivity coefficients.

This work presents a numerical method for evaluating the sensitivity coefficients of models with constrained arguments, with the emphasis on compositions, as these are relevant for the metering of natural gas, hydrogen-enriched natural gas and hydrogen. The method is based on the work from Schay [23], who presented a neat framework for differentiating functions with the input variables subject to an equality constraint. It is shown how this method can be applied for calculating a compressibility factor using the GERG-2008 equation of state [24], which is a model that has arguments subject to a constraint (the amount fractions of the components) and independent ones (temperature and pressure).

2. Law of Propagation of Uncertainty

Measurement uncertainty is propagated using a measurement model. Much of the guidance on evaluating measurement uncertainty is still based on using the LPU for mutually independent quantities [5], equation (10)]. In this work, the form for dependent quantities is used [5], equation (13)], as arguments subject to a constraint are dependent. Consequently, it is assumed that the covariance matrix associated with the arguments subject to the constrained is developed in a fashion that is consistent with that constraint. A covariance matrix associated with a composition has specific features, one of them being that the elements on each row (and for that matter, each column) sum to zero. It was discussed previously that for duly propagating measurement uncertainty, the covariance matrix including the covariances should be used [20]. In this work, such a covariance matrix will be called ’properly formed’.

Using LPU requires calculating the partial derivatives of the measurement model, also known as sensitivity coefficients [5,7]. According to the ISO/IEC Guide 98-3  [5], clause 5.1.3], these coefficients “[...] describe how the output estimate y varies with changes in the values of the input estimates $x_1,x_2,\dots ,x_N$”. For variables that are independent, the partial derivatives are given by the gradient of the function f [22,23,36]

$$\nabla f=\left[{\displaystyle \frac{\partial f}{\partial x_1}},\dots ,{\displaystyle \frac{\partial f}{\partial x_N}}\right]$$

(1)

where the partial derivatives in the gradient can be obtained by analytic differentiation or evaluated numerically [37]. The numerical approach, also mentioned in ISO/IEC Guide 98-3, lends itself also for very complex models, or models in the form of an algorithm. As long as there is a valid implementation of f that returns a value for the output variable y for the input variables $x_1,\dots ,x_N$, the partial derivatives can be evaluated numerically with sufficient reliability for using the LPU [5], clause 5.1.3].

If (a subset of) the $x_1,\dots x_N$ cannot be varied independently because of a relationship between them, then usually also the way in which y responds to changes in the input quantities changes [23]. For the propagation of measurement uncertainty, the behaviour of the function including any constraints on the input quantities is relevant for propagating the measurement uncertainty. The constrained partial derivatives express this behaviour and by virtue of the description of the desired properties of the sensitivity coefficients [5], clause 5.1.3], it is that these constrained partial derivatives are needed for using the LPU.

In many cases, the expressions for the sensitivity coefficients are obtained by analytic differentiation. That works also for situations with constrained input quantities [23]. It can become however prohibitively difficult for complex measurement models, such as equations of state [25,27,33] or in instances where the measurement model takes the form of an algorithm, such as in the calculation of fluid properties at vapour-liquid equilibrium [33,38]. In those cases, a numerical method is desirable, where the variation of the input variables should respect not only any boundaries on the input quantities, but also any constraints. For the amount fractions forming a composition, the boundaries are that for any amount fraction x, $0<x<1$ (situations where $x=0$ or $x=1$ require a dedicated mathematical treatment and are not relevant for calculating properties of mixtures). The constraint on the other hand, formulated as a function g which is equal to zero is given by [23]

$$g\left(\mathbf{x}\right)=\sum _{i=1}^N{x}_i-1=0$$

(2)

where $\mathbf{x}={\left[x_1,\dots ,x_N\right]}^{\top }$ denotes the composition, a vector holding the amount fractions of the components in the mixture. In the case of an equation of state, the function f has two more arguments, namely the pressure p and temperature T. These are not constrained, but have bounds (both are non-negative).

Numerical differentiation methods use finite differences [37,39]. These differences should be formed appreciating the constraint(s) on the arguments, to ensure physically feasible output of the function f. So, the finite differences for p and T should be such that the perturbed pressure and temperature are still valid in their own right. The perturbed composition should still meet the constraint in equation (2).

3. Constrained Differentiation

Analytic differentiation of functions with arguments subject to an equality constraint has been described by Schay [23]. In this work, his method is applied to differentiating functions with a composition, but his method applies to all differentiable functions with constrained arguments, including multivariate ones.

Let $\mathrm{d}\mathbf{x}$ denote a small change (perturbation) of the constrained input quantities $\mathbf{x}$, then the change in the output of a function f, $\mathsf{\delta }y$ is given by [23]

$$\mathsf{\delta }y=\nabla fP\mathrm{d}\mathbf{x}$$

(3)

where P is the orthogonal projection matrix that depends solely on the constraint g, for a composition given in equation (2). The projection P is for a scalar function g given by [23]

$$P=I-\mathbf{n}{\mathbf{n}}^{\top }$$

(4)

where I denotes the $N\times N$ identity matrix and $\mathbf{n}$ the unit normal vector ${(\nabla g)}^{\top }/|\nabla g|$. For the constraint in equation (2),

$$\mathbf{n}={\displaystyle \frac{1}{\sqrt{N}}}{[1,\dots ,1]}^{\top }$$

(5)

and the projection matrix equals $P=I-{\displaystyle \frac{1}{N}}\mathbf{1}$ where $\mathbf{1}$ is an $N\times N$ matrix with ones. Equation (3) shows how the constrained partial derivatives are different from their unconstrained counterparts. A change in the function value equals

$$\mathsf{\delta }y=\nabla f\mathrm{d}\mathbf{x}$$

(6)

in the unconstrained case. Equation (6) is the fundament under the LPU [5,7], when only the first terms of the Taylor expansion of the measurement model f are considered.

To illustrate the difference between the two sets of partial derivatives, consider the calculation of the calorific value of a natural gas as described in ISO 6976. The formula for calculating the calorific value ($H_{\mathrm{mix}}$) takes the form [18]

$$H_{\mathrm{mix}}=\sum _{i=1}^N{x}_i{H}_i$$

(7)

where $H_i$ denotes the calorific value of component i ($i=1,\dots ,N$), $x_i$ its amount fraction and N the number of components in the gas mixture. The calorific values are unconstrained in this measurement model, whereas the amount fractions of all components in the mixtures are constrained as discussed previously. The constrained partial derivatives with respect to $x_i$, formally denoted by the operator ${\displaystyle \frac{{\partial }^{*}}{\partial x_i}}$, are [23]

$${\displaystyle \frac{{\partial }^{*}{H}_{\mathrm{mix}}}{\partial x_i}}=H_i-\mathrm{bar}H$$

(8)

where

$$\mathrm{bar}H={\displaystyle \frac{1}{N}}\sum _{i=1}^N{H}_i.$$

Both the constrained partial derivatives of equation (8) and their counterparts ignoring the constraint lead to exactly the same standard uncertainty $u\left(H\right)$ as the calculations in ISO 6976 [18]. The proof rests on the properties of the covariance matrix associated with the composition [20,21].

Let the Jacobian J be defined as the gradient of the multivariate, differentiable function $\mathbf{f}$, i.e., [36,37]

$$J_{ij}={\displaystyle \frac{\partial f_i}{\partial x_j}}$$

then the covariance matrix associated with the output vector $\mathbf{y}$ is obtained as [7], clause 6.2.1]

$$V_{\mathbf{y}}^{\left(0\right)}=J{V}_{\mathbf{x}}{J}^{\top }.$$

(9)

The proof that instead of J a sensitivity matrix can be used composed of constrained partial derivatives can be given as follows. From the equality constraint g, it follows that the constraint is in the form of an $(N-1)$-dimensional hyperplane $\mathcal{H}$. For compositions, the constraint is given in equation (2) and the normal vector $\mathbf{n}$ in equation (5). Let the matrix D be defined as

$$D=\left[{\mathbf{d}}_1,\dots ,{\mathbf{d}}_{N-1},\mathbf{n}\right],$$

(10)

where ${\mathbf{d}}_i$ denote $(N-1)$ independent directions in the hyperplane $\mathcal{H}$, each perpendicular to the normal vector of the constraint (equation (5)).

Now, equation (9) can be written as

$$V_{\mathbf{y}}^{\left(0\right)}=\left(JD\right)\left(D^{-1}{V}_{\mathbf{x}}{D}^{-T}\right){\left(JD\right)}^{\top }.$$

(11)

Let $V_D=D^{-1}{V}_{\mathbf{x}}{D}^{-T}$. Let L denote the Cholesky factor of $V_x$, ${\mathbf{e}}_N={[0,\dots ,0,1]}^{\top }$ and $\mathbf{0}={[0,\dots ,0,0]}^{\top }$. The constraint $\mathbf{n}\top \mathbf{x}=\kappa$, where $\kappa$ is some exactly known constant implies that

$$\begin{array}{c}\hfill 0=\mathrm{var}\left(\kappa \right)=\mathrm{var}\left({\mathbf{n}}^{\top }\mathbf{x}\right)={\mathbf{n}}^{\top }\mathrm{var}\left(\mathbf{x}\right)\mathbf{n}={\mathbf{n}}^{\top }{V}_{\mathbf{x}}\mathbf{n}={\mathbf{n}}^{\top }{L}^{\top }L\mathbf{n}={\parallel L\mathbf{n}\parallel }^2\end{array}$$

As $L\mathbf{n}=\mathbf{0}$ and from the properties of the covariance matrix is follows

$$V_{\mathbf{x}}\mathbf{n}=\mathbf{0}.$$

(12)

From the definition of D it is evident that

$$D^{\top }\mathbf{n}={\mathbf{e}}_N\mathrm{and}{D}^{-\top }{\mathbf{e}}_N=\mathbf{n}.$$

So,

$$V_D{\mathbf{e}}_N=\left(D^{-1}{V}_{\mathbf{x}}{D}^{-\top }\right){\mathbf{e}}_N=\mathbf{0}.$$

This result implies that the last column, and by symmetry the last row, of $V_D$ contain only zeros. Let ${\tilde{V}}_D$ be the matrix $V_D$ without the last column and row.

Let a directional derivative $f_{\mathbf{d}}^{\prime }$ at $\mathbf{x}$ be defined as [40]

$$f_{\mathbf{d}}^{\prime }\left(\mathbf{x}\right)=\underset{\epsilon \to 0}{lim}{\displaystyle \frac{f(\mathbf{x}+\epsilon \mathbf{d})-f\left(\mathbf{x}\right)}{\epsilon }}$$

(13)

where $\mathbf{d}$ denotes the direction. Now, it is shown that the method does not rely on a particular choice of vectors ${\mathbf{d}}_1,\dots ,{\mathbf{d}}_{N-1}$ by looking at the product $J_D=JD$. From calculus, for the directional derivative $f_{\mathbf{d}}^{\prime }$ of f in direction $\mathbf{d}$ it holds that [40], section 12.1]

$$\begin{array}{cc}\hfill f_{\mathbf{d}}^{\prime }& =J\mathbf{d},\hfill \end{array}$$

(14)

so that

$$J_D=JD=(f_{{\mathbf{d}}_1},\dots ,f_{{\mathbf{d}}_{N-1}},f_{\mathbf{n}}).$$

(15)

Note that all columns of $J_D$ can be calculated using (19) or (20) except the last one if $\mathbf{x}$ is limited to values in $\mathcal{H}$. However, as the last row and column of $V_{\mathbf{d}}$ are zero, knowledge about the last column of $J_D$ is not needed. Let

$${\tilde{J}}_D=(f_{{\mathbf{d}}_1},\dots ,f_{{\mathbf{d}}_{N-1}}).$$

(16)

The result of the generalized algorithm equals

$$V_{\mathbf{y}}^{\left(1\right)}={\tilde{J}}_D{\tilde{V}}_D{\tilde{J}}_D^{\top }$$

(17)

Further below, it will be explained in more detail the correspondence with the version of the algorithm presented in Section 4. Using the fact that the last row and column of $V_D$ only contain zeros, it becomes clear that

$$V_{\mathbf{y}}^{\left(1\right)}={\tilde{J}}_D{\tilde{V}}_D{\tilde{J}}_D^{\top }=J_D{V}_D{J}_D^{\top }=V_{\mathbf{y}}^{\left(0\right)},$$

(18)

which shows that the constrained partial differentiation yields the same covariance matrix as using the LPU (see equation (9)) with the Jacobian. As no assumptions on the ${\mathbf{d}}_i$, $1\le i\le N-1$ were made except for their independence and perpendicularity to $\mathbf{n}$, the result of the algorithm is independent on the particular choice for the ${\mathbf{d}}_i$. The choice between using the unconstrained and constrained partial derivatives does however matter for ill-formed covariance matrices, see also Appendix A.

So, the issue in using the partial derivatives for constrained variables is in this instance not so much in that they could not be used when propagating measurement uncertainty using the LPU. The issues are in two other directions, viz.,

• ISO/IEC Guide 98-3  [5], clause 5.2.1] states about the sensitivity coefficients that they “[...] describe how the output estimate y varies with changes in the values of the input estimates $x_1,\dots ,x_N$"
• when using numerical differentiation, the measurement model f may not exist outside the boundaries of the constraint on the arguments.

If the MCM [6,7] is used to obtain values for the sensitivity coefficients, for example for constructing an uncertainty budget, it is affected in a similar manner as using the LPU, for it is generally not possible to vary just one quantity while keeping all others constant, as described in Annex B.1 of ISO/IEC Guide 98-3/Supplement 1  [6].

4. Numerical Constrained Differentiation

Numerical differentiation [36,37] is a practical alternative for analytic differentiation to obtain values for the sensitivity coefficients. This method is introduced here to enable propagating measurement uncertainty through involved models, such as the GERG-2008 [24] and the AGA8-DC92 [26,27] equations of state, and algorithms.

At first glance, it is not straightforward to perturb a composition and still satisfy the condition given in equation (2). Consider a normalised direction vector $\mathbf{q}$. For constrained differentiation with respect to compositions, it is required that $\mathbf{x}+\epsilon \mathbf{q}$ is also a valid composition for a sufficiently small $\epsilon$. In this case, the constraint (2) can be expressed as $g\left(\mathbf{x}\right)={\mathbf{e}}^{\top }\mathbf{x}-1=0$ where ${\mathbf{e}}^{\top }=[1,\dots ,1]$. $\mathbf{x}+\epsilon \mathbf{q}$ is a valid composition if $\mathbf{q}$ is orthogonal to $\mathbf{e}$. From the requirement that $\mathbf{x}+\epsilon \mathbf{q}$ is a valid composition, it follows that the sum of the elements of $\mathbf{q}$ is zero.

To illustrate how a composition $\mathbf{x}$ can be perturbed by a small shift, let us consider a ternary mixture ($N=3$), $\mathbf{x}={\left[x_1,x_2,x_3\right]}^{\top }$. A change in $x_1$ should be compensated by changes in $x_2$, $x_3$ or both, so that the perturbed input vector $\mathbf{x}+\mathsf{\delta }\mathbf{x}$ should again be a composition. Consequently, the sum of the elements of the perturbation $\mathsf{\delta }\mathbf{x}$ should be zero. Amongst the infinitely many options to satisfy this condition, one option is to orthogonally project the vectors ${\left[0,d{x}_2,0\right]}^{\top }$ and ${\left[0,0,d{x}_3\right]}^{\top }$ onto the hyperplane $\mathcal{H}$ containing $\mathbf{x}$ and with as normal vector $\mathbf{n}$, the gradient of the constraint. Using the resulting projection P (see also equation (4)) yields the normalized perturbation directions

$$\begin{array}{cc}\hfill {\mathbf{q}}_2& ={\displaystyle \frac{1}{\sqrt{5}}}{\left[-1,2,-1\right]}^{\top }\hfill \\ \hfill {\mathbf{q}}_3& ={\displaystyle \frac{1}{\sqrt{5}}}\left[-1,-1,2\right]{)}^{\top }\hfill \end{array}$$

such that it becomes physically meaningful and mathematically feasible to assess the change in the property of interest when the composition is changed from $\mathbf{x}$ to either of the compositions

$$\begin{array}{cc}\hfill {\mathsf{\delta }\mathbf{x}}_2& =\mathbf{x}+{\epsilon }_2{\mathbf{q}}_{\mathbf{2}}\hfill \\ \hfill {\mathsf{\delta }\mathbf{x}}_3& =\mathbf{x}+{\epsilon }_3{\mathbf{q}}_{\mathbf{3}}\hfill \end{array}$$

where ${\epsilon }_2$ and ${\epsilon }_3$ are some sufficiently small numbers1.

Given a function f having as one of its arguments a composition $\mathbf{x}$, the derivatives with respect to amount fractions of the components in the composition can be obtained numerically as follows.

For the numerical differentiation, instead of computing the limit in equation (13), a finite, small value for $\epsilon \mathbf{q}$ will be used instead

$$f_{\mathbf{q}}^{\prime }\left(\mathbf{x}\right)\approx {\mathbf{b}}_j={\displaystyle \frac{f(\mathbf{x}+\epsilon \mathbf{q})-f\left(\mathbf{x}\right)}{\epsilon }}$$

(19)

To evaluate the partial derivative exactly at the value of the estimates $\mathbf{x}$, the directional derivative can be approximated by

$$f_{\mathbf{q}}^{\prime }\left(\mathbf{x}\right)\approx {\mathbf{b}}_j={\displaystyle \frac{f(\mathbf{x}+\epsilon \mathbf{q})-f(\mathbf{x}-\epsilon \mathbf{q})}{2\epsilon }}$$

(20)

For the numerical evaluation of the constrained partial derivatives, the following algorithm can be used:

• Choose a small value for $\epsilon$ so that $0<x_i+\epsilon <1$ for all components i
• Choose an orthogonal matrix Q with columns ${\mathbf{q}}_1,\dots ,{\mathbf{q}}_{N-1}$, which are all perpendicular to the vector $\mathbf{e}$
• Compute the directional derivatives $[f_{{\mathbf{q}}_1}^{\prime },\dots ,f_{{\mathbf{q}}_{N-1}}^{\prime }]=\mathbf{b}$ using equation (19), or (20) for symmetric partial derivatives
• Now,

  $$C=\mathbf{b}{Q}^{\top }.$$

  (21)

  where C denotes the sensitivity matrix in the sense of ISO/IEC Guide 98-3  and ISO/IEC Guide 98-3/Supplement 2  [5,7]. (In the case of a scalar function fm C takes the form of a row vector, or, equivalently, an $N\times 1$ matrix.)

The algorithm can be readily further generalised to multivariate functions $\mathbf{f}$, which can be used for propagating measurement uncertainty using the LPU from ISO/IEC Guide 98-3/Supplement 2  [7].

The algorithm is based on the consideration that if the directions ${\mathbf{d}}_i$ (see equation (10)) are chosen normalized and mutually orthogonal, the matrix $D=Q_0$ becomes orthogonal so that $D^{-1}=Q_0^{\top }$. This matrix $Q_0$ enables formulating the particularly form of the algorithm, cf. equation (21). In more detail, the matrix Q from the algorithm in Section 4 corresponds to the first $N-1$ columns of $Q_0$ such that

$$\begin{array}{c}\hfill Q_0=\left(\begin{array}{cc}Q& \mathbf{n}\end{array}\right)\end{array}$$

and let

$$\mathbf{b}={\tilde{J}}_{Q_0}=J_Q=JQ.$$

From equation (9), it follows that

$$\begin{array}{cc}\hfill V_{\mathbf{y}}^{\left(0\right)}& =J{V}_{\mathbf{x}}{J}^{\top }=J\left(\begin{array}{cc}Q& \mathbf{n}\end{array}\right)\left(\begin{array}{c}{Q}^T\\ {\mathbf{n}}^{\top }\end{array}\right)V_{\mathbf{x}}\left(\begin{array}{cc}Q& \mathbf{n}\end{array}\right)\left(\begin{array}{c}{Q}^T\\ {\mathbf{n}}^{\top }\end{array}\right)J^{\top }\hfill \\ \hfill & =J\left(\begin{array}{cc}Q& \mathbf{n}\end{array}\right)\left(\begin{array}{cc}{Q}^T{V}_{\mathbf{x}}Q& \mathbf{0}\\ \mathbf{0}& 0\end{array}\right)\left(\begin{array}{c}{Q}^T\\ {\mathbf{n}}^{\top }\end{array}\right)J^{\top }\hfill \\ \hfill & =JQ{Q}^T{V}_{\mathbf{x}}Q{\left(JQ\right)}^T\hfill \\ \hfill & =\mathbf{b}{Q}^T{V}_{\mathbf{x}}Q{\mathbf{b}}^{\top }\hfill \\ \hfill & =C{V}_{\mathbf{x}}{C}^{\top }.\hfill \end{array}$$

so the algorithm provides a sensitivity matrix C that meets the descriptions in the GUM [5,7]. In the case that J cannot be evaluated because the function f is not defined outside the region defined by the constraint, the covariance matrix associated with $\mathbf{y}$ cannot be defined and evaluated by means of equation (9). In this case however, the sensitivity matrix defined in equation (21) can be used instead. This case is currently not covered by the GUM [5,7].

A complementary viewpoint is obtained by considering the orthogonal projection P which was introduced in equation (4). Noting that an orthogonal projection fulfils $P^{\top }=P$ and $P^2=P$, it follows from equation (12) that for a well-formed covariance matrix [20,21] it holds that

$$V_{\mathbf{x}}=P{V}_{\mathbf{x}}P.$$

As a consequence,

$$\begin{array}{c}\hfill V_{\mathbf{y}}=J{V}_{\mathbf{x}}{J}^T=JP{V}_{\mathbf{x}}P{J}^T=\tilde{C}{V}_{\mathbf{x}}{\tilde{C}}^{\top },\end{array}$$

where

$$\tilde{C}=JP={\left(P{J}^{\top }\right)}^{\top }.$$

(22)

This matrix $\tilde{C}$ as orthogonal projection of J on $\mathcal{H}$ offers an alternative definition to the matrix $C=\mathbf{b}{Q}^{\top }$ as constructed in the algorithm. In the next section it is shown that $\tilde{K}$ is actually identical to K if $V_{\mathbf{x}}$ is of rank $N-1$.

If the degeneracy of $V_{\mathbf{x}}$ is solely due to the constraint, then $V_{\mathbf{x}}$ has rank $N-1$ and its unique Cholesky factor L has also rank $N-1$. Assume that by specific choices of the matrix Q, the projection operator of equation (22) or any other method two constrained sensitivity matrices $C_1$ and $C_2$ would have been obtained, i.e.,

$$\begin{array}{cc}\hfill V_{\mathbf{y}}& =C_1{V}_{\mathbf{x}}{C}_1^{\top }=C_2{V}_{\mathbf{x}}{C}_2^{\top }\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill C_1\mathbf{n}& =\mathbf{0}\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill C_2\mathbf{n}& =\mathbf{0}.\hfill \end{array}$$

(25)

Then $(C_1-C_2)L^{\top }=\mathbf{0}$, and as L has rank $N-1$ and $L\mathbf{n}=\mathbf{0}$, it must hold that $C_1-C_2=\alpha {\mathbf{n}}^{\top }$ for some scalar $\alpha$. In view of equations () and () this is only possible if $\alpha =0$, i.e., if $C_1=C_2$. As the factors $C_1$ and $C_2$ occur twice in (23), the uniqueness is actually up to a sign, so, $C_1=\pm C_2$.

Therefore, if $V_{\mathbf{x}}$ has rank $N-1$, the constrained sensitivity matrix is essentially unique and equal to $\tilde{C}$ from equation (22), i.e., the orthogonal projection of the Jacobian applied from the right side or, equivalently, the transpose of the orthogonal projection of the transposed Jacobian.

The matrix Q can be obtained in different ways. One such way is using a matrix derived from the Helmert matrix [41], which is available as the function ilrBase in the R package ’compositions’ [42]. Alternatively, the (modified) Gram-Schmidt algorithm can be used [43]. For a composition, the following orthogonal matrix can be used. For each column $j=1,\dots ,N$, the first j elements are equal to

$${\displaystyle \frac{-1}{\sqrt{j(j+1)}}}$$

and the following element is equal to

$${\displaystyle \frac{j}{\sqrt{j(j+1)}}}$$

The remaining elements are zero. For $N=5$, using these formulæ,

$$Q=\left[\begin{array}{cccc}-0.70711\hfill & \hfill -0.40825& -0.28868\hfill & \hfill -0.22361\\ 0.70711\hfill & \hfill -0.40825& -0.28868\hfill & \hfill -0.22361\\ 0\hfill & \hfill 0.81650& -0.28868\hfill & \hfill -0.22361\\ 0\hfill & \hfill 0& 0.86603\hfill & \hfill -0.22361\\ 0\hfill & \hfill 0& 0\hfill & \hfill 0.89443\end{array}\right]$$

(26)

Using the matrix Q from (26) and the composition of a five-component natural gas in Table 1, the sensitivity coefficients are evaluated for the molar mass of the natural gas using the model of ISO 6976 [18]. The molar mass is given by [18]

$$M_{\mathrm{mix}}=\sum _{i=1}^N{x}_i{M}_i$$

(27)

where $x_i$ denotes the (normalised) amount fraction of component i, $M_i$ its molar mass and $M_{\mathrm{mix}}$ the molar mass of the mixture. The numerical method given in equations (21) and (26) are compared with the values obtained using the analytic formula, viz. [23],

$${\displaystyle \frac{{\partial }^{*}\mathrm{bar}M}{\partial x_i}}=M_i-{\displaystyle \frac{1}{N}}\sum _{j=1}^N{M}_j$$

In Table 1, the sensitivity coefficients are given using the numerical method ($c_{\mathrm{num}}$) versus analytic differentiation including the constraint on the $x_i$ ($c_{\mathrm{an}}$). The sensitivity coefficients are obtained using equation (19).

The sensitivity coefficients in the ultimate two columns of Table 1 are identical up to four decimal figures, which is sufficient for the propagation of measurement uncertainty using the LPU. In the numerical evaluation of the sensitivity coefficients, $\epsilon$ was chosen to be 1 % of the smallest amount fraction in the composition. This choice ensures that $\mathbf{x}\pm \mathsf{\delta }\mathbf{x}$ are valid compositions. According to ISO/IEC Guide 98-3 [44], clause 5.1.3 NOTE 2], a perturbation in the order of a standard uncertainty is adequate for the purpose of propagating measurement uncertainty. This choice of $\epsilon$ just does that. The only concern can be that for compositions involving components at the μmol mol⁻¹-level, this choice may lead to unacceptable cancellation in the calculation of the differences.

Using a constant value of $\epsilon$ and the first order formula (19), there are in total N calls to the function f needed. For equations like equation (27), this feature is not important, but the more involved the evaluation of f is, the more important this deliberation gets. The differences between the sensitivity coefficients using equation (19) and the ones obtained using analytic differentiation are in the order of 10⁻¹² g mol⁻¹. Using equation (20) instead, thus obtaining the symmetric partial derivatives, the errors are slightly smaller in absolute value but still in the order of 10⁻¹² g mol⁻¹. This computation requires $2N$ calls to the function f. Table 2 shows the same comparison, but now for an 11-component natural gas mixture. The differences between the numerically evaluated sensitivity coefficients and the ones obtained through analytic differentiation are in the order of 10⁻⁹ g mol⁻¹. This performance is also sufficient for propagating the measurement uncertainty when computing the molar calorific value H or the compressibility factor Z.

The calculations can be readily performed in any kind of software, from mainstream spreadsheet software to environments like R [45] or Python [46]. The arguments of a function f that do not belong to the composition can be handled as usual. They can be evaluated numerically directly using the approach described in GUM:1995 [5], clause 5.1.3 NOTE 2].

5. Propagation of Measurement Uncertainty Through the GERG-2008 Equation of State

To illustrate the method, the computation of the compressibility factor Z using the GERG-2008 equation of state [24,25] is revisited. Apart from the model uncertainty stated in ISO 20765-2 [25], the uncertainty associated with the measured composition $\mathbf{x}$, pressure p and temperature T is propagated using the LPU for dependent input quantities from ISO/IEC Guide 98-3 [5]. The normalised natural gas composition is given in Table 3 [20]. The correlation matrix associated with the composition is

$$R_{\mathbf{x}}=\left[\begin{array}{ccccc}1\hfill & \hfill 0.0635& -0.0703\hfill & \hfill 0.0367& -0.1543\hfill \\ 0.0635\hfill & \hfill 1& -0.0605\hfill & \hfill 0.0320& -0.1341\hfill \\ -0.0703\hfill & \hfill -0.0605& 1\hfill & \hfill -0.2531& -0.8782\hfill \\ 0.0367\hfill & \hfill 0.0320& -0.2531\hfill & \hfill 1& -0.1609\hfill \\ -0.1543\hfill & \hfill -0.1341& -0.8782\hfill & \hfill -0.1609& 1\hfill \end{array}\right]$$

(28)

The actual conditions are $T=300.00\mathrm{K}$ and $p=6.20\mathrm{M}\mathrm{Pa}$. The standard uncertainty associated with the temperature is $u\left(T\right)=0.20\mathrm{K}$ and that associated with the pressure is $u\left(p\right)=0.05\mathrm{M}\mathrm{Pa}$. The normalised composition $\mathbf{x}$, expressed in amount fractions, is given in Table 3. The associated correlation matrix ${\mathbf{R}}_{\mathbf{x}}$ is shown in equation (28). The covariance matrix ${\mathbf{V}}_{\mathbf{x}}$ can obtained as [7], definition 3.21]

$${\mathbf{V}}_{\mathbf{x}}=\mathbf{D}{\mathbf{R}}_{\mathbf{x}}\mathbf{D}$$

(29)

where $\mathbf{D}=\mathrm{diag}\{u\left(x_1\right),\dots ,u\left(x_N\right)\}$. To propagate the measurement uncertainty associated with the composition, actual temperature and pressure using the GERG-2008 equation of state, equations (19) and (21) are used to obtain the sensitivity coefficients for the composition. For the temperature and pressure, the process from ISO/IEC Guide 98-3 [5], clause 5.1.3 NOTE 2] can be used. So, evaluating the equation of state for $\mathbf{x},T+\mathsf{\delta }T,p$ and a second time for $\mathbf{x},T,p+\mathsf{\delta }p$, enables evaluating these two partial derivatives, which are approximated by

$${\displaystyle \frac{\partial f}{\partial T}}={\displaystyle \frac{f(\mathbf{x},T+\mathsf{\delta }T,p)-f(\mathbf{x},T,p)}{\mathsf{\delta }T}},$$

(30)

and

$${\displaystyle \frac{\partial f}{\partial p}}={\displaystyle \frac{f(\mathbf{x},T,p+\mathsf{\delta }p)-f(\mathbf{x},T,p)}{\mathsf{\delta }p}}.$$

(31)

In these expressions, $\mathbf{x}$ denotes the composition in Table 3. If desired, they can also be obtained exactly at the value of the estimate like shown in equations (19) and (20) for differentiating f with respect to the amount fractions. For the composition in Table 3 and setting $\epsilon$ equal to 0.01 times the smallest amount fraction in $\mathbf{x}$ ($=0.0002421\mathrm{mol}/\mathrm{mol}$, see Table 4), and using the perturbation directions as defined in (26), the following compositions and compressibility factors are obtained (see Table 4). From these, the compressibility factors at 300.0 $\mathrm{K}$ and 6.20 $\mathrm{M}$$\mathrm{Pa}$ were calculated (see Table 5).

From the quotients $\mathsf{\delta }Z/\epsilon$ in the ultimate row of Table 5, using equation (21), the sensitivity coefficients can be computed (see Table 6).

Using equations (30) and (31), the sensitivity coefficients for temperature and pressure can be obtained. Setting $\mathsf{\delta }T=0.2\mathrm{K}$ and $\mathsf{\delta }p=0.00372\mathrm{M}\mathrm{Pa}$, the compressibility factor at a temperature $T+\mathsf{\delta }T$ is 0.870025. The compressibility factor at a pressure $p+\mathsf{\delta }p$ is 0.869601. From these and the compressibility factor at actual conditions (see Table 5, column labelled ${\mathbf{x}}_{\mathbf{0}}$), the sensitivity coefficients can be obtained, whose values are -1.9326e-8 /$\mathrm{Pa}$ and that with respect to temperature 1.7609 /$\mathrm{K}$. The value of the compressibility factor Z is 0.869672 with a standard uncertainty $u\left(Z\right)$ of 0.000375. This standard uncertainty does not include the model uncertainty, which is stated to be 0.1 % relative. Considering that most predicted compressibility factors will lie within $\pm 0.1\%$ of the predicted value, this relative uncertainty is interpreted as an expanded uncertainty with a coverage probability of 95 %. Assuming a normal distribution, the relative standard uncertainty is 0.05 %, so 0.000435. Combined with the uncertainty contribution from propagating the uncertainty associated with the composition, pressure and temperature provides a standard uncertainty $u\left(Z\right)$ of 0.000574. Expressed as a relative standard uncertainty, it is 0.066 %. So, the magnitude of the model uncertainty is about the same as that arising due to the composition, pressure and temperature measurements.

6. Conclusions

In this work, it was shown how functions with constrained arguments can be numerically differentiated to obtain values for the sensitivity coefficients needed for propagating measurement uncertainty. The proposed procedure was demonstrated on the calculation of the compressibility factor of natural gas at metering conditions using the GERG-2008 equation of state, for which duly propagating the measurement uncertainty from the composition, pressure and temperature is challenging. The numerical differentiation method does not depend on a particular choice of the directions chosen for the differentiation and provides a sensitivity matrix that provides the same covariance matrix for the output quantities of the measurement model as the Jacobian would, if that matrix can be obtained. The method is applicable to all differentiable functions of a composition and can thus be applied to all models used in, e.g., physical chemistry that rely on the composition and is thereby a leap step forward not only for custody transfer and billing, but also for engineering as a whole where safety margins need to be established. The numerical method for obtaining the sensitivity coefficients extends the guidance of the GUM, which (silently) assumes that the partial derivatives can be formed by varying one by one the input quantities of the measurement model. This way of differentiating functions is neither possible in models having a composition as argument, nor for any other measurement models where the arguments are subject to a mathematical constraint. The GUM should be revised in the sense that either its scope should be narrowed to measurement models where the input quantities are unconstrained, or, better, it should duly deal with this kind of measurement models. The utility of the developed numerical method goes way beyond using it for models with chemical compositions only. It was proved that the constrained and unconstrained partial derivatives with respect to the amount fractions of models in standards like ISO 6976 for calculating the calorific value, molar mass and compressibility factor of a mixture provide the same output when propagating measurement uncertainty, provided that a properly formed covariance matrix is used for the input composition. As many of the standards like ISO 6976 do not explicitly require using such a covariance matrix, it matters what partial derivatives are used with the LPU. The constrained partial derivatives provide uncertainties that are closer to the ones obtained with a properly formed covariance matrix, so these are to be preferred. The (unconstrained) partial derivatives in these standards do not live up to the explanations in ISO/IEC Guide 98-3 regarding what they are supposed to express and hence are technically flawed. The numerical procedure is easy to implement on any platform that supports matrix multiplication, which includes mainstream spreadsheet software. The formulæ for the orthogonal matrix needed to perturb the composition can be readily implemented too. It is founded on the same footing as the guidance provided in ISO/IEC Guide 98-3 [5], and can be used with the LPU from ISO/IEC Guide 98-3 and ISO/IEC Guide 98-3/Supplement 2 [7]. In the examples developed, the asymmetric partial derivatives were very close to the symmetric ones, so the former can be used as well. The GUM [5,7] presumes that the symmetric ones are used though. In the example given for the compressibility factor from the GERG-2008 equation of state, the uncertainty contribution from propagating the uncertainty associated with the temperature, pressure and composition measurements was of the same magnitude as the model uncertainty for density and the compression factor as given in ISO 20765-2. The example shows how to meet the requirement in clause 7.6 of this standard that requires that the uncertainty of the input quantities to the equation of state should be duly propagated.