Uncertainty quantification

The Measurements extension propagates the uncertainties of the inputs to the line parameters by linear propagation, the first-order method of Measurements.jl. The extension is loaded by using Measurements.

Measurements extension

The extension defines only methods, of the functions of the package and of SpecialFunctions. The Bessel recurrences used in these methods are those of the NIST Digital Library of Mathematical Functions [6].

MeasurementsExt — Module

LineCableModelsMeasurementsExt

Compute linear propagation of uncertainty. A Measurement input gives series impedances Z [Ω/m] and shunt admittances Y [S/m] with uncertain real and imaginary parts. Define the nominal method and four methods of SpecialFunctions Bessel functions on an uncertain complex argument. The engine calls these functions for Schelkunoff surface impedances and Unified kernels.

For $z = x + iy$, write the scaled functions as $I^s_\nu(z) = \mathrm{e}^{-|x|} I_\nu(z)$ (besselix), $K^s_\nu(z) = \mathrm{e}^{z} K_\nu(z)$ (besselkx), and $J^s_\nu(z) = \mathrm{e}^{-|y|} J_\nu(z)$ (besseljx), with $i^2 = -1$. The unscaled $K_\nu(z)$ is besselk. The Cartesian partial derivatives and their order recurrences are

\[\begin{aligned} \partial_x I^s_\nu &= (I^s_{\nu-1}+I^s_{\nu+1})/2 - \mathrm{sgn}(x)I^s_\nu, \\ \partial_y I^s_\nu &= i(I^s_{\nu-1}+I^s_{\nu+1})/2, \\ \partial_x K^s_\nu &= K^s_\nu - (K^s_{\nu-1}+K^s_{\nu+1})/2, & \partial_y K^s_\nu &= i\partial_x K^s_\nu, \\ \partial_x J^s_\nu &= (J^s_{\nu-1}-J^s_{\nu+1})/2, & \partial_y J^s_\nu &= i\partial_x J^s_\nu - \mathrm{sgn}(y)J^s_\nu, \\ \partial_x K_\nu &= -(K_{\nu-1}+K_{\nu+1})/2, & \partial_y K_\nu &= i\partial_x K_\nu, \\ 2I_\nu' &= I_{\nu-1}+I_{\nu+1}, & 2K_\nu' &= -(K_{\nu-1}+K_{\nu+1}), & 2J_\nu' &= J_{\nu-1}-J_{\nu+1}. \end{aligned}\]

The modified Bessel recurrences are equation 10.29.1 of the NIST Digital Library of Mathematical Functions. The ordinary Bessel recurrence is equation 10.6.1. The derivative at the cusp of the absolute value takes sign(0) = 0.

source

Uncertain inputs

Uncertain values, of type Measurement, are accepted for any number of the data model and of the problem, such as a material property, a dimension of a cable design, a placement in the system, a property of the earth or the operating temperature. The number type T of a design is a Measurement when any of its numbers is uncertain, and the problem promotes its numbers to the same type. In this case, the computation returns Z and Y with uncertain real and imaginary parts.

The correlations between results that depend on the same uncertain input are retained by Measurements. Wires that share the same uncertain radius give correlated line parameters, and the test suite verifies that one uncertain radius shared by two wires gives √2 times the uncertainty obtained with two independent radii of the same value.

Deterministic controls

Γ and the other options are controls of the formulation or of the backend, and take deterministic values. A Measurement value given to an option throws a MethodError when the record of options is built.
The finite-element backend computes the nominal problem and returns a result without uncertainty, as described in Finite element solver.

Example

The following example computes the reference case two_insulated_wires with uncertain values of the wire radius and of the copper resistivity.

Define the materials and the cable design with the uncertain inputs:

using LineCableModels, Measurements
copper = Material(:conductor, 1.7241e-8 ± 1.0e-10, 1, 0.999994, 20, 0.00393) # [Ω·m], [°C], [1/°C]
insulation = Material(:insulator, 1.97e14, 2.3) # [Ω·m]
metal = Terminal(:core, Region(:core_metal, Circle(0.0425 ± 0.0005), copper)) # [m]
design = build(CableDesign, "insulated_wire", metal,
    Region(:core_insulation, Offset(1.0e-3), insulation)) # [m]

Build the system of two wires and compute its line parameters:

system = build(CableSystem, [design, design], [Placement(0.0, -1.0), Placement(1.0, -1.0)];
    connections = [(core = 1,), (core = 2,)], id = "two_insulated_wires", line_length = 1.0)
frequencies = 10.0 .^ range(-1, 7; length = 101) # [Hz]
problem = LineParametersProblem(system, homogeneous(rho = 0.1); frequencies) # [Ω·m]
result = compute(problem, CoaxialFormulation())
k = argmin(abs.(result.frequencies .- 50.0))
Measurements.value(result.frequencies[k]) # [Hz]

52.48074602497726

Real and imaginary parts of the self impedance of the first wire at this frequency:

(R = real(result.Z[1, 1, k]), X = imag(result.Z[1, 1, k])) # [Ω/m]

(R = 6.302e-5 ± 1.1e-7, X = 0.00043292 ± 8.6e-7)

Real and imaginary parts of the self admittance of the first wire at this frequency:

(G = real(result.Y[1, 1, k]), B = imag(result.Y[1, 1, k])) # [S/m]

(G = 4.276e-13 ± 9.3e-15, B = 1.814e-6 ± 2.1e-8)

Uncertainty of the self resistance:

Measurements.uncertainty(real(result.Z[1, 1, k])) # [Ω/m]

1.0684207588225137e-7

Tests of the propagation

The tests in test/measurements/propagation.jl compute the reference cases with uncertain inputs and verify that:

The nominal values of Z and Y equal the deterministic results within 1e-14 of the largest diagonal entry.
The derivatives with respect to each input equal the central finite differences of the deterministic computation with a step of 1e-5, within 1e-5.
The uncertainties of independent inputs add as the root sum of squares, and one uncertain input shared by two wires gives √2 times the uncertainty obtained with two independent inputs of the same value, both within 1e-12.

Computational cost

Each result evaluated inside the quadrature of the earth-return integrals stores its list of derivatives, with one entry per uncertain input, and the computation allocates memory for these lists at each evaluation of an integrand. On one workstation, for the reference case two_insulated_wires at 11 frequencies, the deterministic computation takes 6 ms, against 83 ms with one uncertain material property and 236 ms with six uncertain inputs. The example above computes 101 frequencies in about 1 s, against 0.054 s without uncertainty.