Interval sensitivity analysis of joint dynamics in complex built-up structures

The dynamic properties of joints, or inter-component connections, are a principal source of uncertainty when modelling complex built-up structures. In the present paper, we propose an interval-based sensitivity analysis (SA) to establish the contribution, or influence, of a structure’s uncertain joint dynamics towards the uncertainty of its coupled admittance. We choose an interval SA as it overcomes the limitations of conventional local and global methods; accuracy and computation efficiency, respectively. Furthermore, it avoids the need for detailed probabilistic data to describe each joint, which is often unavailable. The proposed interval SA is based on the dual formulation of the sub-structuring problem and utilises the Sherman–Morrison formula to factor out the contribution of an individual joint. A complex interval is used to represent the uncertain joint stiffness and damping, and precise bounds on the assembly’s complex admittance (also operational response) are determined following an application specific algorithm. The relative change between input (joint stiffness and damping) and output (complex admittance) interval area is chosen as a sensitivity metric and used to rank order the uncertain influence of each joint. The method is illustrated by numerical example.


Introduction
Complex engineering structures (such as vehicles, satellites, domestic products, buildings and more) are built-up from individual components that are joined with connections that rely on frictional contact or other dissipative mechanisms (such as elastomeric elements). These inter-component connections, or joints, are a principal source of uncertainty when modelling such structures.
The types of joints encountered in practical engineering structures are wide ranging, and can include bolted connections, rivets, welds and resilient elements, among many others (Ibrahim and Pettit, 2005). Likewise, their dynamics can vary greatly, from linear spring or beam-like coupling, to non-linear hysteretic, perhaps discontinuous (e.g. frictional), connections (Segalman, 2006;Marques et al., 2016). To accurately model the dynamics of a jointed structure, it is necessary to characterise the dynamics of its constituent joints; experimentally, this is a demanding task (Ren andBeards, 1995, 1998;Čelič and Boltežar, 2008;Noll et al., 2013;Tol andÖzgüven, 2015;Roettgen and Allen, 2017;De Oliveira Teloli et al., 2021), often requiring supplementary numerical models (Saeed et al., 2020b,a), and is subject to considerable experimental uncertainty.
In the development or refinement of a numerical model it is important to establish the contribution of each uncertain parameter towards the uncertainty of the model output. This is the role of Sensitivity Analysis (SA) (Saltelli et al., 2008). Methods of SA can be broadly categorised as local or global. Local methods consider linearisation of the model and are appropriate only when low level uncertainty, or weak non-linearity, is present. Based on sampling strategies, global methods consider the entire range of model behaviour and take into account non-linear effects due to large parameter variation.
When dealing with problems of uncertainty, there are two further classifications, or viewpoints: probabilistic or non-probabilistic. Conventional SA is probabilistic, describing uncertain parameters by their probability distributions. One of the principal difficulties in evaluating the uncertainty of a complex system by probabilistic means is the acquirement of detailed statistical input data. In the absence of such data, complex statistical analysis cannot be justified. It is this issue that lead to the development of nonprobabilistic, or 'possibilistic', analysis methods, which include interval analysis (Rao and Berke, 1997;Moens and Vandepitte, 2004;Muscolino et al., 2014;Łasecka-Plura and Lewandowski, 2017;Lei et al., 2019;Imholz et al., 2020;Zhao et al., 2021), convex modelling (Ben-Haim and Elishakoff, 2013) and fuzzy set theory (Chen and Rao, 1997;Yin et al., 2016). As the name suggests, possibilistic methods are not concerned with establishing the probability of an outcome occurring, rather whether or not is it possible, regardless of its likelihood.
Given the experimentally demanding nature of joint characterisation, it is often challenging to obtain the necessary statistical information to conduct a probabilistic analysis of uncertainty. For this reason, possibilistic analysis methods offer a convenient alternative. In the present paper, we consider, in particular, an interval-based analysis of joint dynamics, where uncertain joint parameters (stiffness and damping) are treated as unknown but bound by upper and lower limits.
In 2007, Moens and Vandepitte introduced the concept of interval sensitivity theory for analysing systems subject to interval-based uncertainty. In contrast to conventional local SA, interval SA is not limited to the first order behaviour of the system understudy. Rather, it considers the entire range of dynamics permissible within the parameter's interval bounds. In brief, their proposed interval SA considered the relative change in absolute interval widths on the input and the output side of the problem. Utilising a mixed optimisation and interval arithmetic strategy, a numerical example was presented whereby the frequency response function of a lumped parameter truck model was subject to an interval SA with respect to its lumped mass values.
In the present paper, we propose a (conceptually) similar interval SA methodology to identify (and rank order) the uncertainty contribution of individual joints within a complex built-up structure.
The proposed analysis is based on the dual formulation of the sub-structuring problem and utilises the Sherman-Morrison formula to factor out the contribution of an individual joint. The considered joint dynamics are limited to those of linear spring-damper like connections. Though simple in form, joints of this type have been used successfully to improve estimates of coupled assembly dynamics (Saeed et al., 2020b,a). The stiffness and damping of each joint is represented by a complex rectangular interval, and is used to establish precise interval bounds on the coupled assembly's complex admittance (also its operational response). To rank order the influence of each joint, we consider the relative change between the input (joint stiffness and damping) and output (complex admittance) interval area as a sensitivity metric. These developments are illustrated by a numerical example.
Having introduced the context of this paper, its remainder is organised as follows. Section 2 begins by introducing the governing component-based equations for complex structures with flexible interface dynamics. Section 3 provides some detail on interval analysis, including complex interval representations, and describes the concept of interval SA. Section 4 considers the interval assessment of joint variability, before section 5 demonstrates the proposed joint interval SA by numerical example. Finally, section 6 draws some concluding remarks.

Dual representation of flexible interface dynamics in complex built-up structures
When analysing the properties of complex structures, and the influence of joint dynamics, it is convenient to adopt a component-based approach where each structural component (s) is represented by its free interface admittance matrix Y (s) . This can be obtained numerically by inversion of the dynamic stiffness matrix Z (s) where M, K and C represent the discretised mass, stiffness and damping matrices of the component (Petyt, 2010), or equivalently, assuming proportional or structural damping, by modal summation where Φ is a truncated matrix of mass normalised mode shapes, Λ is a diagonal matrix containing the squared natural frequencies, C is a modal damping matrix that describes the component's damping characteristics (its form will depend on the type of damping considered), and I is an identity matrix (Ewins, 1984). Alternatively, Y (s) can be determined experimentally using standard measurement techniques (International Organization for Standardization, 1994). The components that comprise built-up structures may be classed as either active, or passive. Active components contain vibration generating mechanisms, and include motors, engines, pumps, etc. Passive components simply transmit or radiate energy. A generic two component activepassive assembly is depicted in Figure 1. Each component is represented by two sets of DoFs, the interface set c and the internal set a or b. This picture may be readily generalised to an arbitrary number of components.
where N denotes the number of connections, and k and r are the joint stiffness and damping. To satisfy equilibrium, the interface force g is given the form where λ are the unknown joint forces and B is a signed Boolean matrix. The presence of B T ensures that the forces applied to neighbouring components are equal and opposite, thus satisfying equilibrium. The unknown joint forces λ are found using the relaxed continuity condition as follows. The interface force g is substituted into the EoM where Y is the uncoupled admittance matrix, to give Equation (11) is substituted into the relaxed continuity condition (equation (6)) Solving for λ then yields Finally, substituting equation (13) into equation (11) leads to, which is an equation for the coupled admittance of the assembly. Note that in the presence of a rigid connection, γ n → 0.
In the present paper we are interested in the sensitivity/ uncertainty of the coupled admittance Y C due to the unknown, but bounded, interface flexibility γ n .

Interval analysis
Interval analysis falls within a family of non-probabilistic methods for evaluating uncertainty in complex systems (Langley, 2000;Faes and Moens, 2020). One of the main difficulties when evaluating the uncertainty of a complex system by probabilistic means is the acquirement of detailed statistical input data. This is particularly so for joint dynamics. In the absence of such data, complex statistical analysis cannot be justified. It is this issue that has lead to the development of non-probabilistic, or 'possibilistic', analysis methods. With an interval-based analysis, the uncertain parameter x is not treated statistically. Rather, it is taken to be bounded, such that where □ and □ denote lower and upper bounds, respectively. As an interval, {x} is characterised by a pair of numbers representing its lower and upper bounds where we define the interval mid-point and half width, respectively, as Unlike a probability distribution, the interval bounds provide no detailed information on the likelihood of a particular value of x, only that its true value must occur within the stated bounds.
There are two common ways to extend the interval representation to complex variables (required to represent joint stiffness and damping). In one, the complex interval is represented by a circular disk in the complex plane, characterised by its centre point and radius (Petkovic and Ljiljana, 1998;Rump, 1999). In the other, the real and imaginary parts of the complex number are represented by independent real intervals. This leads to a rectangular interval in the complex plane (Petkovic and Ljiljana, 1998;Rokne and Lancaster, 1971;Lohner and Gudenberg, 1985). In the present paper, we will adopt the rectangular complex interval representation as it enables the independent assignment of joint stiffness and damping.
Adopting the centre point and half width notation above, a complex rectangular interval takes the form an example of which is shown in Figure 2.
Having defined the notion of a real/complex interval, it is the role of interval analysis to compute the upper and lower bounds of functions whose inputs are intervals For simple functions, this can be achieved by redefining the basic arithmetic operations + À × ÷ as interval operations Å . Ä / (Petkovic and Ljiljana, 1998). This leads to the so-called natural interval extension of the function. For more complex problems, global optimisation (Neumaier and Hansen, 1994) or vertex (Dong and Shah, 1987) methods can be used, among other problem specific algorithms.
Care should be taken when using the natural interval extension of a function. Whilst for real intervals the operations Å . Ä / are optimal (i.e. they provide the smallest possible bounds that enclose the true range of outcomes), this is not so for complex intervals. Specifically, complex division. The natural extension of complex interval division can be expressed as where {x} and {y} are each taken to represent a complex interval. The complex interval obtained as per equation (20) is not optimal. If standard interval arithmetic operations are applied to complex division, the numerator and denominator intervals are treated independently. For example, the instance of {y r } in the numerator is considered independent to that being squared in the denominator fy 2 r g. The same issue occurs with {y i }. This so-called 'dependency problem' can lead to a huge overestimation of the interval. Similar dependency issues are often encountered for more complex functions where specific intervals occur more than once and cannot be factored appropriately. This represents one of the main limiting factors of conventional interval arithmetic analysis. For complex division, special algorithms are required to compute the optimal rectangular interval (Rokne and Lancaster, 1971;Lohner and Gudenberg, 1985;Mahmood and Soylu, 2020).
As a final remark on interval analysis, it is noted that non-interval quantities can readily be incorporated within interval arithmetic by considering an interval whose upper and lower bounds are of equal value,

Interval sensitivity analysis
Sensitivity analysis (SA) is a tool used to allocate, or apportion, the uncertainty in a system's output to that of its input parameters (Saltelli et al., 2008). Methods of SA can be broadly categorised as local or global. Local methods assume a small variance on the input parameters and consider a linearisation of the system's dynamics. Such methods are typically based on partial derivatives (obtained analytically or by numerical approximation). Local methods are advantageous in terms of computational efficiency, but limited in terms of accuracy, especially for non-linear systems with reasonable levels of input variance. Global methods overcome this limitation by sampling the model output over the space of input parameters. Whilst computationally more expensive, this approach ensures that model non-linearity and input parameter distributions are accounted for. Hence, more a robust analysis is achieved.
The local and global methods described above are limited, respectively, by the assumption of linearity and computational efficiency. An alternative approach, overcoming these limitations, is obtained by considering an interval-based SA (Moens and Vandepitte, 2007). Whilst interval SA does not provide a detailed analysis like the global approach, it avoids the need to sample the entire input space; instead of representing each input as a probability distribution, intervals are considered and the bounds of the output are determined.
To illustrate the concept of interval SA we consider the function f( ), whose inputs and output are treated as complex rectangular intervals. We consider the interval output resulting from each input individually Based on equation (21), a measure of sensitivity can be obtained in several ways. Moens and Vandepitte (2007) considered the change in absolute interval widths on the input (z) and the output (w) side of the problem and proposed the sensitivity metric (assuming real valued intervals) For reasons discussed later, in the present paper we choose to define a complex interval sensitivity metric instead as the ratio of output interval area to input interval area S n ¼ u n À u n v n À v n x n À x n y n À y n : In this way, S n describes the relative increase in uncertainty considering the entire range of possible outputs in the complex plane; given the same range of possible inputs, a larger value of S n indicates a greater range of possible outputs and a greater sensitivity of the function f( ).
The reason for this particular choice of sensitivity metric is that, for the problem considered here (assembly admittance in the presence of joint uncertainty), the output interval takes a very similar form to the Jacobian matrix that describes the partial derivatives of the assembly admittance with respect to joint flexibility. Hence, with this definition of interval sensitivity there is a clear connection to linear/ gradient-based SA.

Interval assessment of joint variability
In this section we will consider the interval-based assessment of uncertain joint dynamics, and its application in determining the maximum and minimum bounds of an assembled structure's admittance matrix.
As shown in section 2, the dual formulation, in the presence of a flexible interface, leads to the following equation for the assembled admittance matrix where the interface flexibility matrix Γ is defined as With attention focused towards the influence of a single joint, it is convenient to separate its contribution by means of the Sherman-Morrison formula (i.e. a rank-one matrix update) (Hager, 1989) where A is any invertible matrix, tr(□) is the trace operator, and X is a matrix of rank one. To apply the Sherman-Morrison formula, the interface flexibility matrix Γ is split into two parts; Γ (n) , which contains the nth joint flexibility of interest, and Γ ðñÞ , which contains all remaining terms Note that Γ (n) is a rank-one matrix, so can take on the role of X in equation (26). Hence, the matrix inverse in equation (24) can be rewritten in the form where A ¼ BYB T þ Γ ðñÞ . Note that Γ (n) is a single entry matrix. Its only non-zero value, the interface flexibility γ n , can be brought out as a scalar term, and Γ (n) replaced by P n , which is a zero matrix bar the nth diagonal entry whose value is one. Equation (28) then reduces to is a complex scalar value, and Z n = k n + ir n is the dynamic stiffness of the nth interface connection. Note that g n is a non-linear function of γ n ; the joint, although linear in its dynamics, has a non-linear effect on the assembled structure's admittance matrix. This non-linearity limits the application of local/gradient-based SA. Note that this nonlinearity should not be confused with response nonlinearity, whereby the structural response is dependent upon the amplitude of said response. Rather, it is system's coupled admittance matrix that is related non-linearly to the dynamic properties of the joint.
Herein we propose that the interface stiffness Z n be treated as a complex rectangular interval {Z n }, such that its bounds are supposed to contain all possible values of joint stiffness and damping. The complex interval {Z} = {k} + i {r} is characterised by a pair of real intervals, namely where x and x denote the lower and upper bounds of the respective interval, k c and r c are the stiffness and damping mid-points, and k Δ and r Δ are the associated half widths. Treating the interface stiffness as a complex interval {Z n }, and substituting equations (29) and (30) into equation (24), leads to the interval expression or where fz n g ¼ fZ n gÅtrðP n A À1 Þ, and {Y C } is a complex interval admittance matrix whose bounds enclose all possible values obtainable given the range of permissible joint dynamics specified by the interval {Z n }. It is important that the complex interval division of equation (33) is performed optimally using an appropriate algorithm. Unlike complex addition, subtraction and multiplication, the natural interval extension of complex division is not optimal and can lead to a large overestimation of the resulting bound (Rokne and Lancaster, 1971;Lohner and Gudenberg, 1985;Mahmood and Soylu, 2020). In the present paper, we propose an application specific algorithm (see Appendix B) to obtain the precise bounds of the complex interval {Y C } (also the magnitude interval {|Y C |}).
Note that Y ðñÞ C contains no intervals, and represents the structure's admittance when γ n = 0, that is, when the nth connection is made rigid. Hence, fΔY ðnÞ C g represents an interval matrix that describes the effect, or influence, of the nth joint's relaxation on the coupled admittance.
It is further noted that the matrix product within fΔY ðnÞ C g represents the assembly's Jacobian matrix, J n , with respect to the interface flexibility γ n , evaluated for a rigid connection (see Appendix A) This Jacobian matrix would form the basis of a local/ gradient-based SA (Meggitt, 2022). From the above, it is clear that fΔY ðnÞ C g is proportional to the system Jacobian, but includes the interval contribution of the uncertain joint through the interval term {z}. It is for this reason that we choose to define interval sensitivity as per equation (23), as opposed to the derivative of the interval width as per equation (22) (Moens and Vandepitte, 2007).

Response interval
In the presence of an active component (i.e. a vibration source) the operational response u C is often of greater importance than the assembly admittance. Hence, we look to extend the interval treatment above to account for an assembly's operational activity.
The operational response of a complex assembly can be expressed in the form where Y C is the coupled admittance matrix, and f S is a vector of so-called blocked forces (Moorhouse et al., 2009). The blocked force is an independent source quantity that describes its operational activity. It has become an internationally standardised method for experimentally characterising vibratory sources (International Organization for Standardization, 2019), and is used extensively across several engineering sectors, including automotive ( As per equation (36), the complex interval of an operational response can be obtained by where it is noted that Ä represents an interval matrix multiplication. Like the assembly admittance, the interval response is made up of two parts; the first is that of the assembly when the nth joint is considered rigid u ðñÞ C , and the second is an interval contribution describing a relaxation of the nth joint fΔu ðnÞ C g.

Numerical example
In this section we will demonstrate the interval SA for joint dynamics through a numerical example. The system considered is illustrated in Figure 3. It is an 11 DoF massspring-damper system (numeric values are available via the supplementary MATLAB codes). It comprises two subsystems (green and blue) connected by three joint elements (red), whose properties are considered uncertain.
In this numerical example we will (a) for a single joint, demonstrate the interval computation for the real, imaginary and magnitude admittance bounds, and compare the results against those obtained by a sample-based estimation, and (b) rank order the uncertainty contribution of each joint towards the assembly admittance using the proposed sensitivity metric.
In what follows, each joint is represented by its centre point stiffness and damping, and parametrised by an uncertainty factor α □ = □ Δ /□ c which controls the width of the interval bound fk n g ¼ ½ð1 À α k, n Þk c, n ð1 þ α k, n Þk c, n fr n g ¼ ½ð1 À α r, n Þr c, n ð1 þ α r, n Þr c, n , Unless otherwise stated, the centre point stiffness and damping values are set at 1 × 10 4 (N/m) and 100 (Ns/m), respectively, with uncertainty factors α k,n = 0.5 and α r,n = 0.25.
To validate the interval computation a sample-based estimate is used. The complex joint interval fZ n g ¼ k k Â Ã þ i r r ½ is logarithmically sampled between the lower and upper bounds of its real and imaginary parts, generating a total of n = 100 2 sample points per frequency. These samples are used to compute n complex admittance matrices, the max/min values of which are taken as estimates of the complex interval {Y C }. The magnitude of each complex admittance is also taken to estimate {|Y C |}.
We begin by considering the interval bounds due to a single joint; joint 1 is represented by the complex interval {Z 1 }, whilst joints 2 and 3 are considered deterministic (α k,2 = α k,3 = 0, α r,2 = α r,3 = 0). Shown in Figure 4 are the interval bounds of the real and imaginary parts of the complex admittance fY C 1;11 g. Orange markers denote the sample-based estimates and solid lines the interval computations obtained as per the algorithm described in Appendix B. To plot the intervals on a logarithmic scale, magnitudes are taken for both the real and imaginary parts; positive values are shown in blue, and negative values in green.
It is clear from Figure 4 that the interval computation (see Appendix B) correctly determines the bounds on both the real and imaginary admittance. These bounds will be used shortly to compute the interval sensitivity metric, as per equation (40).
Shown in Figure 5 are the magnitude bounds of the complex admittance fjY C 1;11 jg. Also shown, in black, is the mid-point response, obtained by considering k 1 = k c,1 and r 1 = r c,1 . As expected, we see that the interval computation successfully captures the shifts in resonant and antiresonant frequencies due to the unknown interface flexibility; this is particularly evident around 19 Hz where the second system resonance varies between 17-20 Hz.
It is interesting to note that the lower bound of the ≈70 Hz anti-resonance, goes to zero (hence cannot be plotted). This indicates that there exists a combination of joint stiffness and damping that yields zero transmission through the interface. Figures 4 and 5 provide a visual indication of the sensitivity of the assembly admittance with respect to joint 1. It is clear that the certain modes have a greater sensitivity to joint uncertainty than others.

Interval sensitivity analysis
With attention focused towards the influence of joint uncertainty, our definition of complex interval sensitivity (see equation (23)) takes the form  where Z n = k n + i r n is the joint dynamic stiffness, and ΔY ðnÞ C i, j describes the dynamic contribution of the nth joint.
The sensitivity S n describes the relative amplification of uncertainty in terms of the complex interval areas that describe the joint and coupled admittance. The greater the metric, the greater the influence that joint has on the assembly admittance. By comparing the sensitivity metric of each joint it is possible to identify the dominant source of uncertainty in the system. For clarity, we present all sensitivity metrics in a normalised form where a value of approx. 1 indicates that the uncertainty of the assembly admittance is almost entirely due to the nth joint. Similarly, a value of approx. 0 indicates that the joint contributes negligibly to the admittance uncertainty.
It should be noted that the sensitivity metric considers the influence of each joint independently (remaining joints are treated as deterministic with centre point stiffness and damping); interaction effects between multiple joints are not considered.
Shown in Figure 6 are the magnitude interval bounds with respect to each uncertain joint (taking α k,n = 0.5, α r,n = 0.25), alongside their corresponding (normalised) sensitivity metrics. For clarity, sample-based estimates have been omitted. Figure 7 similarly presents the magnitude interval and sensitivity metrics for the reduced uncertainty levels of α k,n = 0.05 and α r,n = 0.025.
The relative widths of the magnitude bounds provide a visual indication of the influence of each uncertain joint, though their respective sensitivities depend also on the relative size of the input interval; this cannot be interpreted from the magnitude bounds alone. Rather, the sensitivity metric b S n provides an indication of the relative increase in uncertainty by comparing the input and output interval areas for each joint. Though a rank ordering, this result can be used to aid model refinement, by targeting those joints with the greatest sensitivity. From Figure 6, we might conclude that above approx 50 Hz, joints 1 and 3 have the greatest sensitivity, whilst below 50 Hz joints 1 and 2 appear to dominate. Careful comparison of the sensitivity metrics presented in Figures 6 and 7 uncover an interesting result. The sensitivity with respect to a particular joint depends on the size of the input interval; the sensitivity metric captures the nonlinear contribution of the uncertain joint's dynamics towards the system response (which arises due to equation (30)). Though only a minor effect in the present case study, more complex systems can exhibit much greater effects (Meggitt, 2022). This ability to capture non-linear behaviour is a principal advantage of interval-based SA over conventional gradient-based methods, which rely on a linearisation of the system dynamics.
For more complex systems, the sensitivities of several joints can be combined to form a single sensitivity value. For example, the stiffness and damping in each coordinate direction (x, y, z) for a single connection point might be combined to yield a sensitivity value for the connection as a whole (Meggitt, 2022).

Conclusions
Joints are a principal source of uncertainty when modelling the dynamics of complex built-up structures. To better understand their influence, it is of interest to apportion the uncertainty of a system's output (FRF or operational response) to that of its uncertain joint dynamics by means of sensitivity analysis (SA). To overcome the limitations of conventional 'probabilistic' SA, we have proposed an interval-based SA for joint dynamics in complex built-up structures.
We introduce an application-specific algorithm that enables precise bounds on the real, imaginary and magnitude admittance to be determined, based on the interval uncertainty of each joint. To characterise the sensitivity of the system to a particular joint, we propose a complex interval sensitivity metric, representing a ratio of complex interval areas on the input and output of the system. This metric indicates the relative increase or decrease in uncertainty, in an interval sense, due to the dynamics of the system. This metric can be used to identify those joints that contribute unacceptably to the total uncertainty of the system, for example, for model refinement purposes.
The advantages of the proposed interval SA, over conventional (probabilistic) local and global methods are many: 1) As per the interval paradigm, there is no need for a detailed probabilistic description of each joint, which is often unavailable. Rather, only the upper and lower limits of joint stiffness and damping are required.
2) The interval SA is not limited to the linear behaviour of the system (as per local gradient-based methods); all system non-linearities contained within the input interval bounds are accounted for; hence, the method can be termed 'pseudo-global'. The ability of interval SA to capture the system's non-linear dependence on the joint dynamics was demonstrated by numerical example, where it was shown that the sensitivity metric provided a different rank ordering when subject to different levels of input uncertainty. 3) Unlike standard global methods, which sample the input parameter space over their probabilistic representations, the proposed interval SA is computed analytically, and so avoids the computational burden of sampling over a multi-dimension input space.
Though, a notable limitation of the proposed interval SA is its inability to treat multiple joints at a time, both in terms of uncertainty propagation (obtaining the admittance bounds) and SA. Consequently, the interaction between multiple uncertain joints cannot be accounted for, nor can the admittance bounds due to several uncertain joints be determined. Extension of the proposed interval SA to the multi-joint case was considered beyond the scope of this paper, and is left as further work.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) received no financial support for the research, authorship, and/or publication of this article.

Supplementary material
All MATLAB codes are available on request.

Jacobian of the coupled admittance
In this section we derive the Jacobian matrix for the coupled admittance with respect to a specific joint flexibility. We begin by considering the complex differential of the assembly admittance Noting that interest lies in the differential interface flexibility d(Γ), the above reduces to Application of the product rule, dðABCÞ ¼ dðAÞBC þ AdðBÞC þ ABdðCÞ, then leads to Recalling the differential of a matrix inverse (Hjørungnes, 2013), dðAÞ À1 ¼ ÀA À1 dðAÞA À1 , yields which reduces to In equation (46), we can interpret dY C as the small (complex) change in admittance, due to a small (complex) change in the interface flexibility dðΓÞ. Importantly, dðΓÞ represents a small change in each matrix entry. Focussing attention on the nth interface connection in particular, such that dðΓÞ → P n dγ n , where P n is a zero matrix bar the nth diagonal entry, whose value is one, equation (46) becomes where dγ n is a scalar differential of the nth interface flexibility, and J n is a complex Jacobian matrix containing the partial derivatives of the coupled assembly admittance, with respect to the nth interface flexibility which represent open, closed and open circular domains, respectively. Note that the direction of the inequalities (i.e. whether they are open or closed) depends on the relative sign of x, x, y, and y. The intersection of the four domains describes a conformal mapping from the complex plane z = x + iy (specifically, the complex rectangular interval {z}), to the domain of w = u + iv. An example is illustrated in Figure 8. The complex interval {w} represents the smallest rectangle that encloses the domain w. The subsequent addition of Y ðñÞ C in equation (50) acts to translate each domain in the complex plane. This is equivalent to translating the centre points of each circular domain, for example, Hence, evaluation of the complex interval admittance {Y C } reduces to finding the smallest enclosing rectangle that surrounds the intersection of the translated circular domains described above.
An estimate of magnitude interval admittance {|Y C |} can subsequently be found by considering each vertex of the complex rectangularinterval {Y C }, taking smallest/largest absolute distance from the origin as the minimum/ maximum interval bound. However, for certain combinations of J and z this can lead to a large over estimation of the true magnitude bound, as the complex (rectangular) interval {Y C } naturally encloses regions that lie outside the true domain of Y C . An alternative approach considers the relative size of each circular domain and finds the precise magnitude bounds. This approach is detailed in the following subsection. The interval bounds of the domain w (or the translated domain Y C ) will correspond to the intersection points of the four circular domains, or the extremum values of these domains, depending on their relative sizes.

B.1 Interval algorithm
As an example, from Figure 8 we see that if the radius of the mapped x > x domain (in blue) exceeds that of the y > y domain (in green), the maximum value of u will be the extremum of the y > y domain (its radius). In contrast, if the radius of the y > y domain exceeds that of the x > x domain, then the maximum value of u will correspond to their intersection point.
The intersection points of the four circular domains can be found by solving the following simultaneous equations for each domain pair ðu int À x 1 Þ 2 þ ðv int À y 1 Þ 2 ¼ x 2 1 þ y 2 1 (60) ðu int À x 2 Þ 2 þ ðv int À y 2 Þ 2 ¼ x 2 2 þ y 2 2 , where x n and y n describe the centre point of each circular domain's boundary, and x 2 n þ y 2 n its radius squared. Solving the above for u int and v int we obtain the coordinates of the two intersection points. To account for the addition of Y ðñÞ C in equation (50), the intersection points are simply translated according to ðu int þ realðY ðñÞ C Þ, v int þ imagðY ðñÞ C ÞÞ. The extremum points of each circular domain's boundary depend on the complex numerator J, which acts to rotate the domain about the origin, and can be obtained by considering the problem geometrically. As an example, the mapped x > x domain (closed, blue) has the extremum u points where realðY ðñÞ C Þ accounts for the subsequent addition of Y ðñÞ C in equation (50), and denotes the angle of rotation. Similarly, the extremum v points are Similar expressions can be obtained for the remaining domains.
The bounds of the complex interval fY C g ¼ Y ðñÞ C ÅðJ n /fz n gÞ are obtained as follows: 1. The coordinates for all critical (intersection and extremum) points are computed using equations (60)-(66), etc. 2. Coordinates are checked against inequalities 55-58; those that satisfy all inequalities are retained. Note, care should be taken to avoid rounding errors which may cause critical points to fall outside an inequality. 3. Finally, the complex interval bounds are obtained by taking the max/min of the remaining coordinates. Figure 9 is an example of the interval computation, including the mapped domains and critical points.

Shown in
To determine the optimum magnitude interval {|Y C |} care should be taken to avoid coordinates that lie within the complex rectangular bound of {Y C }, but outside the true bound of Y C . To obtain the upper magnitude bound, the intersection and extremum points are each are checked against inequalities 55-58. After discarding the coordinates that do not satisfy all inequalities, the upper bound is found by calculating the magnitude of the remaining coordinates and taking the maximum. For the lower magnitude bound, a similar procedure is followed, however, to avoid underestimating the lower bound a further set of points must be included; those on the domain boundaries that lie closest to the origin.
Note that care should be taken when apply the interval computation in the presence large background noise. Based on physical arguments, the complex interval {z n } should not contain the origin. This would imply that some combination of joint stiffness and damping would yield an infinite mobility. However, in the presence of background noise, this scenario is possible. Interval computations whereby 0 + i0 2 {z n } should be discarded.