Response bounds for complex systems with a localised and uncertain nonlinearity

Predicting the vibration response of complex nonlinear structures is a significant challenge: the response may involve many modes of the structure; nonlinearity precludes the use of efficient techniques developed for linear systems; and there is often uncertainty associated with the nonlinear law, even to the extent that its functional form is not always problem in a novel way. The method exploits the fact that nonlinearities are often spatially localised, and seeks the bestand worst-case system response with respect to a chosen metric by regarding the internal nonlinear force as an independent excitation to the underlying linear system. Constraints are used to capture what is thought to be known about the nonlinearity without needing to specify a particular law. This paper focuses on the case of systems with a single point nonlinearity but with arbitrarily complex underlying linear dynamics, driven by a sinusoidal force excitation. Semi-analytic upper and lower bounds are proposed for root-mean-square response metrics subject to constraints which specify that the nonlinearity should be a combination of (A) passive, (B) displacement-limited, and / or (C) force-saturating. The concept of ‘equivalent linear bounds’ is also introduced for cases where the response metric is thought to be dominated by the same frequency as the input. The bounds corresponding to a passive and displacement-limited nonlinearity are compared with Monte Carlo experimental and numerical results from an impacting beam test rig. The bounds corresponding to a passive and force-saturating nonlinearity are compared with numerical results for a frictiondamped beam. The global upper and lower bounds are satisfied for all input frequencies but are generally found to be rather conservative. The ‘equivalent linear bounds’ show remarkably good agreement for predicting the range of root-mean-square velocity responses. Finally, the principle of Maximum Entropy is used to estimate the response distributions, which was found to give surprisingly good agreement with experimental and numerical data. & 2016 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
There is a need for efficient and pragmatic tools that can predict the vibration response of complex nonlinear structures, for example in the contexts of turbine blades with frictional contacts [1]; oilwell drilling dynamics [2]; vehicle brake squeal [3]; 'buzz, squeak and rattle' problems in the automotive industry [4]; mooring line dynamics [5]; or stringed musical instruments [6]. The presence of nonlinearity precludes direct application of efficient linear methods for predicting the system response, and full time-domain simulations can be computationally intensive. In addition, nonlinear interaction laws are often highly uncertain to the extent that even the functional form of the law may not be known, for example friction contact laws [7].
For many examples of system the nonlinearities are spatially localised. The system dynamics can then be described in terms of the (few) degrees of freedom associated with the nonlinearity while still efficiently accounting for the linear dynamics of the coupled system. This property is commonly utilised with time-domain integration of Finite Element models (e.g. [8]), coupling impulse response functions to nonlinearity (e.g. [8][9][10][11]), and with the Harmonic Balance Method (e.g. [12]). This enables systems with local nonlinearities to be simulated more efficiently, but accounting for uncertainty can still be computationally demanding.
There is a large and growing body of literature investigating methods for predicting the effects of uncertainties in structural dynamics: helpful special issues include [13][14][15]. The approaches can be categorised according to the type of uncertainty to which they apply with main classifications including parametric/non-parametric; probabilistic/non-probabilistic; or epistemic/aleatory types of uncertainty. Some interesting discussion of these and other classifications can be found in [16,17], but for the purposes of this article the main distinction is made between parametric and non-parametric uncertainties. Parametric uncertainty descriptions are applicable to systems where the governing equations are known but whose specific parameter values are not precisely defined (with uncertainty specified by bounds, fuzzy numbers, or probability density functions). Non-parametric uncertainty descriptions are needed when there are unknowns associated with the governing equations of the system (sometimes referred to as 'model uncertainty', see [17]).
There are two challenges common to methods that tackle parametric uncertainty: they require multiple simulations to be carried out to predict the response of the system for different choices of system parameters; and the functional form of the governing equations of the system needs to be specified. There are methods emerging that begin to tackle the first issue. For example, in the context of probabilistic uncertainties Peherstorfer et al. [18] use importance sampling together with a combination of surrogate models and high fidelity models to obtain an efficient estimate of the response statistics. That there is a need for this kind of multi-resolution algorithm itself highlights the difficulty, and there is still a need for multiple simulations of a high resolution model. Fuzzy arithmetic is another method applicable to non-probabilistic types of uncertainty, but as described by Moens and Hanss [19] the efficiency of fuzzy arithmetic methods is still limited by the number of simulations needed to estimate response bounds for different levels of uncertainty membership. This is because the response bounds are found by optimisation, or for the upper bound 'anti-optimisation': in other words numerical optimisation is used to search the admissible set of parameters for the extreme responses. For both of these example methods the governing system equations need to be pre-specified as 'knowns' and parametric uncertainty methods intrinsically cannot account for 'model' uncertainty.
These challenges make non-parametric methods appealing: unfortunately their domain of applicability is often narrow and/or a great deal of care is needed when interpreting the results. For example Statistical Energy Analysis (SEA) enables effective prediction of the steady-state mean response and its variance for linear systems at high frequency [20], i.e. when there is significant modal overlap. The method has successfully been combined with deterministic Finite Element modelling, effectively extending its applicability to low-and mid-frequencies [21]. It remains a challenge to model nonlinear systems and only a few attempts have been made to extend SEA to nonlinear systems, e.g. [22]. Another interesting non-parametric approach is using random matrices and the principle of Maximum Entropy, or 'entropy optimisation' [23]. In this approach the system matrices (mass, stiffness and damping) are taken to be uncertain, with their mean represented by deterministic values and probability density function assigned using Maximum Entropy. Response distributions can then be calculated efficiently for linear systems. One challenge with this approach is the difficulty of interpretation: the admissible set of system matrices is not necessarily constrained by physical arguments. For example, the ensemble of possible uncertainties can include systems where remote parts of the system are directly coupled, leading to non-physical results [17]. In itself this is not a failing of Maximum Entropy, rather that the approach does not readily lend itself to including physical constraints. In addition, this approach is mainly focused on uncertainty associated with the underlying linear system rather than nonlinear interactions.
The non-parametric random matrix method has been extended to some examples of nonlinear systems, e.g. [24,25]. Ritto et al. [25] consider the nonlinear dynamics of an oilwell drillstring, with particular emphasis on developing a nonparametric description of the uncertain nonlinear bit-rock interaction. The approach uses a random matrix to represent the nonlinear interaction law: constraints are applied to the random matrix such that it satisfies three key properties: it must be positive-definite (power is always dissipated); its expectation is chosen to represent the deterministic 'mean' model; and it is not a singular matrix. The probability density function for the random matrix used for simulations is taken to have the Maximum Entropy distribution resulting from these constraints. However, in order to obtain a stochastic prediction it is necessary to solve the stochastic system equations for many realisations of the random matrix. In addition, there is a need to specify a 'dispersion parameter' that represents the degree of uncertainty. Finally the constraints on the random matrix itself are only loosely based on physical arguments: power dissipation and the notion of a mean model. This paper presents a new approach for predicting the response of complex structures when uncertainty is associated with local nonlinear interactions, by using a non-parametric approach combined with the concept of anti-optimisation to obtain response bounds. The approach also draws on Maximum Entropy as a way to estimate response distributions. The novelty of the approach described in this paper is that analytic (or semi-analytic) bounds and the distribution of the response are estimated without needing to solve for a nonlinear response. This paper builds on previous studies [26,27] which used numerical optimisation to obtain the 'worst' response, subject to constraints derived from physical arguments. This study takes a new approach and extends the work presented in [28], with the aim of seeking semi-analytic response bounds subject to a minimal set of constraints. The scope of this paper is limited to sinusoidal excitation of systems which have a single point discrete nonlinearity that is uncertain, and for which the linear dynamics may be complex (i.e. systems with many modes). Extensions to multiple discrete nonlinearities are the subject of ongoing research.
The goals of this paper are to: propose semi-analytic response bounds for complex systems with spatially local and uncertain nonlinearities; test whether the resulting bounds give useful predictions that are not overly conservative; attempt to estimate the response distribution; and to apply the approach to academic case studies representing common kinds of nonlinearity.
The adaptation of the anti-optimisation framework is described in Section 2. Hypothesised global response bounds are derived in Section 3, and the concept of 'equivalent linear bounds' is introduced in Section 4. Global lower bounds are derived in Section 5. The bounds are generalised to relax several assumptions in Section 6. The results are then applied to two specific test cases: an impacting beam (Section 7) and a friction-damped beam (Section 8). For both test cases the principle of Maximum Entropy is also used to estimate the response distribution.

Framework for predicting response bounds
The framework is conceptually straightforward: replace the actual nonlinearity with an external excitation source whose objective is to maximise the response, subject to constraints that capture some properties of the original nonlinear system. This approach can be formulated as an optimisation problem: the key advantage over other methods is that by replacing the nonlinearity with an external excitation source, it is only necessary to evaluate the linear system response when seeking response bounds. The approach differs from a traditional parametric interval analysis because rather than placing bounds on model parameters, constraints (or bounds) are placed on higher level properties of the nonlinearity such as power input, response displacement, or maximum force. The framework requires a specification of three components: a representation of the underlying system dynamics; an output response metric for which bounds are sought; and high-level constraints that describe general properties of the nonlinearity. These components and the resulting optimisation problem are presented in the following sections.

System representation
Consider a general linear system with a single localised nonlinearity, illustrated by the feedback diagram in Fig. 1. The underlying linear system can be characterised by a matrix of impulse response functions dðtÞ. This linear system is subject to a total input force fðtÞ, so that the output states are given by yðtÞ ¼ dðtÞÃfðtÞ (where n denotes convolution). The total input force f is composed of the external input force f i and an internal nonlinear force f nl , such that This paper tackles the case of a general linear system with a purely sinusoidal input force at Site (1) and with a single point nonlinearity at Site (2). It is initially assumed that the two sites are distinct: the collocated case is considered in Section 6.1. The spatially discretised input force vector can be arranged so that f i ¼ ½ f 1 ðtÞ 0 0 … T and the nonlinear internal force vector arranged so that f nl ¼ ½0 f 2 ðtÞ 0 … T . Choosing a sinusoidal input at frequency ω a gives where F 1a may be a complex amplitude, with the subscript 'a' associated with the input frequency ω a .

Response metric
A response metric M needs to be defined that captures an output of interest. This choice defines the quantity whose bounds and distribution are sought. The particular choice of response metric depends on what is of most interest for a given application: in this study, the metric is taken to be the steady-state root-mean-squared (RMS) displacement, velocity and acceleration at Site (1) where the input is applied: where α ¼ 0; 1; 2 defines whether the chosen metric is displacement, velocity or acceleration respectively. The time window t to t þ T is taken to encompass the steady-state response. RMS quantities are commonly used measures of vibration amplitude, so this represents a widely applicable choice of metric.

Constraints on the nonlinearity
The nonlinear force is the result of a localised nonlinearity which would normally be viewed as a nonlinear function of one output state y nl , illustrated by the dashed line in Fig. 1. In the proposed approach, the nonlinear force f nl will be considered to be an external input force and general properties of the nonlinearity will be represented by either equality or inequality constraints. This allows an ensemble of nonlinearities to be defined in a flexible way, without having to specify its functional form.
The choice of constraints depends on what is assumed to be known about the nonlinearity. As more information is known, more constraints can be imposed and the predicted bound becomes less conservative (as demonstrated in [26]). Some choices of constraints for particular applications can lead to discontinuous constraint functions that require numerically challenging optimisation to find solutions (e.g. [27]). The aim of this paper is to find semi-analytic bounds for common kinds of nonlinearity: this requires finding numerically convenient choices that also have a clear physical basis. With this in mind, three properties of common nonlinearities will be considered: (A) passive, i.e. the nonlinearity cannot add energy to the system; (B) displacement limited; and/or (C) force limited. These choices are motivated by: 'impacting', 'buzzing' or 'rattling' system, i.e. structures with a local, passive, displacement-limited nonlinearity (A þB); and frictiondamped systems, i.e. structures with a local, passive, force-saturating nonlinearity (A þC).
More formally, the constraint functions hðf nl ; y nl Þ ¼ ½h A h B h C T must satisfy hðf nl ; y nl Þ r 0 and have been chosen as follows: where h A is a constraint on the average power input from the nonlinearity, h B is the displacement-limit, and h C is the force saturation constraint. The time window t to t þT is taken to encompass the steady-state response, y max is the largest amplitude displacement (e.g. from a clearance nonlinearity), and f max is the saturating force (e.g. from a frictional nonlinearity). The reason that average power is chosen for the constraint h A , rather than instantaneous power, is that passive nonlinearities can readily be envisaged that have non-trivial reactive components. Note that no equality constraints are used for the present study. The physical origin of the passive constraint h A is clear: in many cases the presence of nonlinearity does not add energy to the system, e.g. friction dampers; gears with backlash; fluid damping; or nonlinearity from structural joints. The displacement constraint h B also has a clear physical basis: when systems have a gap or clearance (such as oilwell drillstrings within a borehole), then the response is limited to be within the clearance. This constraint is harder to define in practise, as when the clearance is reached then there will be some finite stiffness associated with the contact and the nominal clearance can be exceeded. Nevertheless, the mindset of this constraint is pragmatic and assumes that the contact stiffness is reasonably large such that the clearance represents a practical limit. The force saturation constraint h C is based on typical friction behaviour, where friction forces tend to saturate for a given normal preload. This limit can be hard to identify precisely, as friction characterisation is challenging [29], but the constraint nevertheless represents a pragmatic property based on typical frictional behaviour.

General optimisation problem
The problem to be solved can written as a standard optimisation problem: maximise: where M is the response metric and h is the vector of inequality constraint functions. The specific choices for the present study are defined in Eqs. (1) and (2). The input force f i is assumed to be known. The degrees of freedom for the optimisation could either be taken as the nonlinear force f nl ðtÞ or the associated state y nl ðtÞ. The solution to this optimisation problem can be sought using 'black-box' toolboxes available with standard software packages such as Matlab. This style of analysis has been investigated in previous studies [26,27]. However, the goal of this paper is to explore whether response metrics and constraints can be chosen that allow semi-analytic solutions to be found that still give useful predictions.

Global upper bounds on system response
The global solution is sought to the optimisation problem defined in Eq. (3), with response metric from Eq. (1), subject to the constraints defined by Eq. (2). The three constraints will be considered independently: in each case the upper bound response must occur when the constraint is active, and when pairs of constraints are considered together then the response bound is simply given by the more restrictive bound. This approach seeks fundamental bounds on the response of a structure given bounds on general properties of the nonlinearity: this is distinct from a traditional interval analysis which seeks bounds given possible intervals of system parameters.

Passive nonlinearity
It is hypothesised that the upper bound response for a passive nonlinearity is obtained if the force f 2 extracts the maximum possible power from the excitation frequency ω a and injects it into another frequency ω b which has the maximum gain G from input power P in to response metric M: This is based on the result that for a finite available power, the maximum response occurs when it is applied at a single 'optimal' frequency (see Appendix A).
Therefore, the nonlinear force f 2 ðtÞ has two frequency components: one at the input frequency ω a and the other at the injection frequency ω b : where F 2a and F 2b are the complex amplitudes of the input and injection components respectively. This abbreviated notation will be used for the subsequent analysis, where the subscript a or b denotes the corresponding frequency component of the complex amplitude. The underlying linear system dynamics can be represented as a frequency response function matrix relating the displacements and forces at Sites (1) and (2): The maximum energy is extracted at Site (2) when F 2a is equivalent to an impedance match (see Appendix B). This is achieved when: with Z ¼ ðiω a D 22a Þ Ã and ðÁÞ Ã denotes complex conjugate. Following the abbreviated notation Y 2a is the ω a component of the complex amplitude of displacement at Site (2). Physically this would correspond to the nonlinearity presenting an impedance match that optimally absorbs incoming wave energy. From Eq. (6) this leads to: which gives the peak power extracted in a cycle: The average power is half this value and it is taken that the nonlinearity construes to use all of it to excite a new frequency ω b . Therefore the peak power at this new frequency is also P h , which can be written in terms of the injected frequency ω b : The aim is to choose ω b to give the largest gain from injected power to RMS displacement (α ¼0), velocity (α ¼1) or acceleration (α ¼2): Now because F 1b ¼ 0, then Y 1b ¼ D 12b F 2b (from Eq. (6)). Combining this with Eq. (9) and rearranging to separate ω a and ω b terms gives: The worst-case frequency ω b can be found by maximising Bðω b Þ, which can be found by evaluating the function for a high resolution vector of frequencies and finding the maximum. The worst-case input force F 2b is found by combining Eqs. (9) and (10): Finally the upper bound is given by: Note that the response makes no assumption of the periodicity of the response: the choice of ω b is free. Therefore the scope of this bound is very general and is an estimate of the global limit on the metric M for any passive nonlinear system coupled at Site (2). The metric is taken to be an average property of a system under steady-state conditions, which could be periodic, quasi-periodic, or even chaotic.

Displacement-limited nonlinearity
It is hypothesised that the upper bound response for a displacement-limiting nonlinearity is obtained when the constraint is fully active, i.e. when the response y 2 ðtÞ is a square wave of amplitude y max and with a fundamental frequency that causes the maximum response at Site (1). To achieve this requires choosing the force f 2 ðtÞ to have a frequency component at the driving frequency ω a which pins the system at Site (2) (i.e. cancels the response due to the input force at Site 1), and a harmonic series of frequency components with fundamental frequency ω b : The coefficients F 2nb ðnω b Þ of the harmonic series can be chosen such that the displacement y 2 ðtÞ at Site (2) is a square wave: The upper bound response occurs when the fundamental frequency ω b of the square wave is chosen to maximise the response metric. It is straightforward to compute the response metric for a square wave displacement input over a range of fundamental frequencies and find the maximum gain from the square wave amplitude to the response metric.
For the system to be pinned at Site (2), then Y 2a ¼ 0 (where Y 2a is the complex amplitude of the ω a contribution of y 2 ).
From Eq. (6) this gives: with corresponding response at Site (1): The worst-case RMS response at Site (1) is given by: with Y 2nb given by Eq. (17). It is interesting to note that the coefficients scale according to Y 2nb p 1=n and the frequencyaveraged function D 12 =D 22 is constant (for a sufficiently broad frequency average). Therefore this upper bound only converges when α¼0, i.e. when the response metric is RMS displacement. In other words, just imposing a displacement constraint does not impose an upper bound on the output RMS velocity or acceleration.

Force-saturating nonlinearity
It is hypothesised that the upper bound response for a force-saturating nonlinearity is obtained when the constraint is fully active, i.e. by assuming that the force f 2 ðtÞ is a square wave of amplitude f max . This can be expressed as a Fourier Series with fundamental frequency ω b : with coefficients given by The upper bound response occurs when the fundamental frequency of the square wave ω b is chosen to maximise the response metric. It is straightforward to compute the response metric for a square wave input force over a range of fundamental frequencies and find the maximum gain from the square wave amplitude to the response metric. Therefore the upper bound RMS response at Site (1) is given by: with F 2nb given by Eq. (22). The coefficients scale according to F 2nb p 1=n and the frequency-averaged function D 12nb p 1=n 2 .
Therefore this upper bound converges for all α¼0,1,2, i.e. when the response metric is RMS displacement, velocity or acceleration. It is interesting that a force-saturating nonlinearity gives upper bounds for a broader set of metrics and is numerically better conditioned than a displacement-limiting nonlinearity. This is consistent with intuition: a constraint on displacement places no bound on its derivatives, while a constraint on force necessarily limits the resulting accelerations.

Equivalent linear upper bounds on system response
It is hypothesised that the bounds above are global solutions to the optimisation problem of Eq. (3): maximising the RMS response subject to energy, displacement or force constraints. As will be seen in Sections 7 and 8, these bounds can lead to overly conservative upper bound predictions. Therefore this section considers 'equivalent linear bounds': adding the constraint that the nonlinear force only has a component at the same frequency as the input frequency ω a . This might seem a severe restriction, but it is directly analogous to the 'Describing Function' approach, i.e. the Harmonic Balance method with only the fundamental frequency retained (e.g. [30]). It is often a valid assumption for displacement or velocity response to input force, rather than the acceleration response, as there is a natural roll-off with frequency for these output metrics. However, it is not expected to be valid for acceleration output metrics which are high-frequency weighted. The same is expected to hold for the equivalent linear bounds: with useful upper bound predictions expected for displacement-or velocity-based metrics, but not for acceleration-based metrics. The distinction between these equivalent linear bounds and the describing function approach is that the equivalent linear bound effectively identifies a 'worst case' describing function as a function of frequency which need not represent a physical system.

Passive nonlinearity
The worst-case passive linear response is readily derived by assuming that a passive linear system is coupled at Site The maximum right-hand side of the inequality occurs when C þD 22a j jis minimised. Both the underlying linear system and the coupled system are passive, so imag D 22a f gr0 and imag C f gr0. Therefore the maximum response occurs when In other words, the worst-case coupled compliance acts to ensure a resonance condition by adding either stiffness (real C f g40) or mass (real C f go 0). Therefore, the equivalent linear upper bound is: recalling that α¼0,1,2 corresponds to RMS displacement, velocity or acceleration.
As noted above, this bound does not constrain the sign of real C f g. Yet it may well be known that the nonlinearity is stiffness-based, only allowing compressive forces that oppose displacement, such that real C f gZ 0. This modifies the equivalent linear bound to:

Displacement-limited nonlinearity
If the displacement at Site (2) has amplitude y max , then from Eq. (6) the force F 2a is given by: where θ is the phase of the displacement response and is not specified.
The magnitude of the response at Site (1) is then given by: y max e iθ (28) and so the equivalent linear upper bound is: If the nonlinearity is constrained to be compressive, then the nonlinear force can only oppose displacement and can be expressed in terms of a compliance such that F 2a ¼ ÀY 2a =C, with real C f gZ 0. The amplitude of the response at Site (1) is given by Eq. (24), duplicated here for clarity: The bound for the system with combined displacement and compression constraints is non-trivial. Therefore, in the interest of finding semi-analytic bounds the compression constraint is considered in isolation (though the actual bound for the combined constraints could readily be found by optimisation for the complex value of the compliance C). If real D 22a f gZ0, the maximum response amplitude occurs when real C f g ¼ 0 and imag C f g ¼ Àimag D 22a f g. If real D 22a f gr0 then the response is unbounded as C ¼ ÀD 22a . Therefore the worst-case displacement-limited response allowing only compressive compliance is given by: where M 0 B is given by Eq. (29), and:

Force-saturating nonlinearity
From Eq. (6), the complex amplitude of the response at Site (1) to input forces F 1a and F 2a is given by: Constraining the nonlinear force such that F 2a j jrf max gives: Therefore the equivalent linear upper bound response for a force-saturating nonlinearity is given by:

Global lower bounds on system response
The above analysis seeks the 'worst-case' response subject to three common constraints. It is straightforward to reverse the question and seek the lower bound 'best-case' response. This can be useful when determining whether or not it is even possible to achieve a design specification. The lower bound response can answer questions such as: can a passive (or otherwise constrained) system attached at a particular location ever reduce vibration levels below some threshold?; or where do I need to attach a passive system so that it becomes possible to reduce vibration to acceptable levels? Then further thought is needed as to how to actually approach this limit. It also gives insight into what may be causing a high lowerbound by identifying which type of constraint is active.
In principle, the global lower bound response is significantly more straightforward than the upper bounds as the nonlinear force that yields the smallest response can only have a frequency contribution at the same frequency as the input. This is because any other frequency contributions can only increase an RMS-based response metric. The bounds for the passive, displacement-limited and force-saturating cases are presented in the following sections.

Passive nonlinearity
To find the best-case passive linear response it is easiest to begin by using the inverse form of Eq. (6): where K ¼ D À 1 . The force F 1 at Site (1) is given by: The best-case passive linear response is readily derived by assuming that a passive linear system is coupled at Site (2) with compliance CðωÞ such that F 2a ¼ ÀY 2a =C. The magnitude of the force F 1a j j is given by: The best-case response can be framed as finding the maximum input force required for a given displacement. The maximum right-hand side of the inequality occurs when 1=C þ K 22a is minimised. Both the underlying linear system and the coupled system are passive, so imag K 22a f gZ 0 (because K is the inverse of D) and imag C f gr0. Therefore the maximum force occurs when 1=C ¼ Àreal K 22a f gþ0i. In other words, the best-case coupled compliance acts to ensure an anti-resonance condition by adding either stiffness (real C f g4 0) or mass (real C f go 0). This is somewhat counterintuitive: it may have been expected that the best case response occurred when the added compliance absorbed most energy, but this does not yield the minimum response at a single point location.
Therefore, the global (and linear) best-case passive response is: This bound does not constrain the sign of real C f g. Yet it may well be known that the nonlinearity is stiffness-based, only allowing compressive forces that oppose displacement, such that real C f gZ0. This modifies the best-case linear response to:

Displacement-limited nonlinearity
The force F 1 at Site (1) is given by Eq. (37): When Y 2a is limited to y max , then the upper bound force is given by: Therefore the best-case response is given by: Note that if F 1a j jr K 12a j jy max then it is possible for the nonlinear force to pin the system at Site (1) without violating the displacement constraint.

Force-saturating nonlinearity
Eliminating Y 2 from Eq. (37) gives: The best-case response is found when F 1a is maximised for a given displacement Y 1a : This gives the best-case response: Note that if F 1a j jr K 12a K 22a f max then it is possible for the nonlinear force to pin the system at Site (1) without violating the force constraint.

Generalisation to different cases
The analysis above assumes that the input force and nonlinearity are located at different sites on the structure; and that the nonlinearity is 'grounded', acting between one degree of freedom and a fixed reference point. A discussion of the effect of relaxing these assumptions follows.

Collocated input force and nonlinearity
If the input force at Site (1) and the nonlinear force at Site (2) are at the same location, then the response is given by: where DðωÞ ¼ D ij ðωÞ. Following through the same working as Section 3, it is reassuring that the global bounds presented for the passive, displacement or force constraints are all still valid. Note that defining DðωÞ ¼ D ij ðωÞ causes some terms to simplify, most notably the displacement-limited bound from Eq. (20) becomes: where Y 2nb were the coefficients for a square wave. Recalling from Section 3.2 that this sum only converges when α¼0 gives the trivial result that: More care is needed to find the equivalent linear bound for the passive constraint when the input and nonlinearity are collocated. The bound M 0 A above is not invalid, but is unnecessarily conservative. Let F 2 ¼ ÀY 1 =C with passive compliance such that imag C f gr 0. Substituting into Eq. (48) gives: Therefore the maximum response is given when 1 þ D=C is at a minimum under the constraint imag C f gr 0. By a geometric argument, the minimum occurs when the complex phasor D is orthogonal to 1 þD=C, which can be expressed: Rearranging gives the solution: Substituting into Eq. (51) gives the equivalent linear passive collocated bound: The equivalent linear bounds for the displacement and force constraints are unchanged when assuming collocation.

Nonlinearity between a pair of degrees of freedom
It is often the case that a single nonlinearity occurs between a specific pair of degrees of freedom of a structure, or connecting two structures. Assume that a nonlinearity acts between a pair of degrees of freedom at Sites (2) and (3): the frequency-domain displacement vector Y and nonlinear force vector F nl can be defined to be As before, the dynamics of the underlying linear system can be written in the frequency domain in terms of a frequency response matrix: If the nonlinear interface is intrinsically massless then F 2 ¼ ÀF 3 , and a new variable defining the relative displacement can be defined: This allows a row and column of the frequency response matrix to be eliminated, giving a reduced frequency response matrix D 0 with a nonlinearity associated with a single degree of freedom Y 0 2 . This allows the analysis above to be carried out exactly as before.

Application to an example test structure I: impacting beam
The bounds derived above are now applied to an example test structure. The purpose of this study is to explore how effective the bounds are for predicting the behaviour of a system, rather than to identify artificial cases that attempt to approach or exceed the predicted bounds. Two common kinds of system are investigated: structures involving impact (passive, displacement-limited nonlinearity), and friction-damped structures (passive, force-saturating nonlinearity). The focus of this section is on a system involving impact, and forms the main comparative study. Section 8 presents a shorter summary of the equivalent comparisons for a friction-damped structure.
For both test cases the underlying linear system chosen for study is a free-free beam and is based on an experimental test rig illustrated in Fig. 2 and Fig. 3. This represents a 'complex' (i.e. multi-mode) linear system which can be coupled to a spatially localised nonlinearity with adjustable properties that can be investigated numerically and experimentally.
The nonlinearity for the impacting beam consisted of a pair of end-stops with an adjustable gap and contact material. This allowed an ensemble of data to be collected for different 'uncertain' contact properties. This system can also be simulated numerically to explore regions of the parameter space beyond practical limits imposed by the test rig. The data collected from this impacting beam experiment is compared with predicted response bounds using the passive, displacement-limited constraints.

Experimental test rig
The beam was suspended by string rather than clamped in order to minimise boundary nonlinearity: clamping typically introduces frictional and contact nonlinearity. The strings were connected to the beam at the nodes of the first 'free-free' mode of vibration, which are asymmetric due to additional mass at Site (1), as discussed below. The strings were mounted to the beam via rubber isolation to prevent the taught string from making intermittent contact with the beam under dynamic loads. The length of the beam is L ¼1 m, the width w¼40 mm and the thickness for most of the length was d ¼3 mm. A cylindrical neodymium rod magnet (10 Â 30 mm) was mounted onto the beam at Site (1) (x 1 ¼0.005 m), which together with a coil (RS: 357-766) clamped to the bench provided the driving force. The magnet was viewed as part of the system dynamics of the beam, and there was no contact between the magnet and the coil. This arrangement was found to be particularly effective at producing a linear excitation with a large throw: single frequency input force driving the beam to a peak-peak displacement of around 10 mm with the largest harmonic in the coil current less than 40 dB below the fundamental frequency, and the largest harmonic in the acceleration response of the linear beam generally less than 40-60 dB below the fundamental frequency.
The beam was made thicker at the driving location (Site 1) such that d ¼10 mm for the first 40 mm of the beam's length. The intention was to create a mounting platform for the magnet and accelerometer (B&K Type 4369) such that the curvature of the beam at their interface was negligible under dynamic loads, in order to minimise nonlinear effects associated with their attachment. The additional mass of the mounting platform, magnet and accelerometer breaks the symmetry of mode shapes of the beam.
Deliberate nonlinearity was introduced at Site (2) (x 2 ¼0.995 m) in the form of a pair of end-stops. Hemispherical pins of either rubber or polycarbonate were used at the contact point, which were mounted on force transducers clamped to a heavy base. The nominal clearance was imposed using feeler gauges in the range 0.05-2 mm. The material choice and gap variation enabled a Monte Carlo study with adjustable contact parameters, in order to generate an ensemble of data for comparison with predictions.

Nonlinear test method
Each experimental test consisted of a sinusoidal input force at a fixed frequency for 10 s which was usually sufficient for a satisfactory steady-state to be reached. The amplitude of the input force was smoothly ramped up over 0.5 s in order to reduce starting transients. In addition, each test was repeated with and without an initial Gaussian pulse (band-limited to 100 Hz) that guaranteed the initiation of contact at Site (2) with the end-stops: it was observed that intermittent contact could sometimes be sustained even when the purely linear system response was insufficient to initiate contact. Frequencies were selected at random from a uniform distribution in the range 0.1-100 Hz.
For a given frequency, the test was repeated at increasing amplitude levels. The maximum amplitude in the sequence was chosen based on the free-free response to limit the peak displacement of the beam at Site (1) to 5 mm (based on the measured transfer function D 11 ). Some tests nevertheless resulted in larger amplitudes due to the nonlinearity and the beam occasionally impacted the excitation coil. In these cases, the data was discarded and a sequence at a new frequency was initiated. The test rig was mounted on a wooden bench which was affected by people walking past: an accelerometer on the floor detected when this occurred and repeated the test case as appropriate. A total of approximately 21,000 tests were carried out that were not affected by passers-by or coil impacts.
A quick check of the response spectra revealed that impacts with rubber end-stops excited a 40 dB bandwidth of approximately 700 Hz and polycarbonate end-stops excited approximately 5 kHz. These bandwidths were case-dependent but this provides an approximate guide. The experiments with rubber were logged using a sampling rate of f s ¼10 kHz, and for polycarbonate end-stops f s ¼30 kHz.
Several metrics were extracted from each experiment, including the root-mean-square (RMS) displacement, velocity and acceleration of the beam at Sites (1) and (2). The RMS velocity and displacement were obtained by high-pass filtering the accelerometer signal using a third-order Butterworth filter with low frequency cutoff at 10 Hz; integrating once or twice for velocity or displacement; and then subtracting the mean from the chosen steady-state time window. The accuracy was confirmed by independent measurements with a laser vibrometer (not available for all tests). The metrics were computed using the final second of data before the input force was ramped down: in most cases this represented a reasonably good approximation to a steady state. A summary of the data is shown in Fig. 5: each polygon encompasses a cloud of data points which are shaded according to the excitation amplitude, with dark shades corresponding to low amplitudes.

Numerical model
The linear dynamics of the beam can be represented as a transfer function matrix relating the displacements YðωÞ and forces FðωÞ at Sites (1) and (2)  with variables defined in Fig. 2. Note that Y 2 and F 2 are taken to be collocated: the offset shown in Fig. 2 is approximately 5 mm and is to prevent the accelerometer itself from coming into contact with the end-stops. These transfer functions have been measured and mathematically fitted using standard modal analysis procedures (e.g. [31]). By way of example, the driving point velocity transfer function at Site (2) is shown in Fig. 4(a) and (b) in the frequency ranges 0-100 Hz and 0-5 kHz respectively. The frequency range 0-100 Hz was the same range as used for excitation during the nonlinear tests (described in the next section). The solid line shows the measured transfer function and the dashed line is the modal fit: the peaks were used for the fitting algorithm hence anti-resonance troughs do not match as well. The three first bending modes are visible at 13 Hz, 38 Hz and 76 Hz. The two rigid body swinging modes are at very similar frequencies around 0.8 Hz. The transfer functions were measured using a small instrumented impulse hammer and modal parameters were estimated for the first 72 modes: a good fit was obtained up to 2.5 kHz (19 modes) and an increasingly approximate fit was obtained up to 16 kHz.
The equations of motion for the beam can be written in modal coordinates. For the nth mode: where q n is the modal contribution of the nth mode, dots denote differentiation with respect to time, ω n is the natural frequency, ζ n is the damping factor, Φ T n is a row vector of the spatially discretised mode shape, and f is a column vector of the spatially discretised forces. Recall that the force vector is ordered such that f ¼ ½ f 1 f 2 0 0… T , i.e. with the input force at Site (1) and the nonlinearity at Site (2). The nonlinear force f 2 is a function of the spatial coordinate y 2 , which is straightforward to reconstruct during time-domain simulations using y 2 ¼ Φ ð2Þ qðtÞ where Φ ð2Þ is a row vector of all the mode shape contributions to the y 2 degree of freedom. Numerical integration was carried out using Matlab's built-in ode45 algorithm.
Experimental results indicated that impacts with polycarbonate excited modes up to 5 kHz, requiring 31 modes. The computation time for a single 10 s nonlinear simulation was approximately 15 min: too slow for Monte Carlo simulations. As a trade-off all modes up to 1 kHz were included, requiring 12 modes and giving average computation times of approximately 2 min per simulation. The modal properties (Φ ð1;2Þ , ω n , and ζ n ) at y 1 and y 2 were taken directly from modal fitting of the experimentally measured transfer functions.
For the impacting beam simulations, the nonlinear law was taken to be a spring in parallel with a viscous damper: where k 0 is the stiffness coefficient, p is the exponent of the assumed power-law dependence of the contact (a linear spring corresponds to p ¼1 and Hertzian contact theory corresponds to p ¼1.5), c 0 is the viscous damping coefficient, and y 0 is the clearance. The 'min' and 'max' functions are necessary to ensure that the nonlinear force only acts in compression. To improve equivalence with experimental results, the stiffness coefficient k 0 was based on an effective linear contact stiffness k c , such that for any power-law dependence p the contact force at depth y ref would always be equivalent: giving The reference depth was taken to be y ref ¼ 0:1 mm, and the linearised contact stiffness identified from the experimental results was found to be in the range k c ¼ ½10 4 ; 10 6 . The lower and higher values corresponded to rubber and polycarbonate end-stops respectively. A numerical Monte Carlo study was carried out, varying the nonlinear contact parameters using a uniform distribution over the following ranges: 10 3 ok c o10 6 N m À 1 , 1op o1:5, and 0 oy 0 o 1 À 3 m. The input force was sinusoidal with amplitude F 1a ¼ 0:625; 1:25; 2:5; 5; 10 N. The maximum force applied in the experiments was 2.5 N so this choice deliberately overlapped and extended the range that was possible experimentally.

Results and comparisons
A summary of the experimental data is shown in Fig. 6. Approximately 21,000 data-points were obtained, which are summarised by polygons generated using Matlab's 'boundary' function. Each data-polygon represents a collection of experimentally identified steady-state RMS velocities plotted against input frequency, and the polygons are shaded according to excitation amplitude with dark shades for small amplitudes. If the system were linear and without uncertainty then the responses at a given frequency should span one decade, corresponding to the range of input force amplitudes spanning one decade. This is true for much of the input frequency range, perhaps suggesting a weakly nonlinear system. However, near 10 Hz, 25-40 Hz and 60-80 Hz the behaviour is more clearly nonlinear as the response spans two or more decades. These frequencies correspond to anti-resonances of the beam and demonstrate bistable behaviour: at these frequencies a smooth starting transient means that contact is not initiated and the beam follows the linear 'free-free' response, but when an initial Gaussian pulse is applied then contact is initiated and sustained. In addition, it can be seen that some peaks diminish with increasing input force amplitude, while others become more significant: the clearest example is the peak at 75 Hz which becomes much less significant as input force increases, while a peak at 65 Hz becomes steadily more obvious. The shift in dominance of these peaks occurs as the beam transitions from behaving more like a 'free-free' beam at low input force amplitudes, to a 'free-pinned' beam at high input force amplitudes. Fig. 6 shows the data corresponding to just the lowest amplitude input force for (a) RMS displacement; (b) RMS velocity; and (c) RMS acceleration. Also shown in these plots are the two limiting cases of the predicted linear 'free-free' or 'freepinned' beam as labelled in the figure. At first glance the labels appear to be the wrong way round: adding a constraint should raise the corresponding resonant frequencies. In fact each pinned mode does correspond to the lower free mode, and the first significant peak in the pinned response comes from the low frequency rigid body modes. It is interesting that in Fig. 6(a) and (b) many of the features of the experimental data are captured remarkably well by either the free-free or freepinned predictions, suggesting that the system is 'weakly' nonlinear. In contrast Fig. 6(c) suggests that the system is very strongly nonlinear as the linear predictions under-predict most of the experimental data by an order of magnitude. These features can be understood in terms of the significance of harmonics contributions to the RMS metric: the RMS acceleration weights high frequency contributions that are excited by impacts more strongly than the RMS velocity or displacement.
The numerical simulations can be compared with the experimental results by filtering the data for a given input force amplitude. An exact comparison is not possible because the input force amplitude for the experimental tests was scaled at a given frequency according to the predicted linear response (for practical considerations to avoid damage to the test rig); while the numerical simulations used a constant force across all frequencies and allowed larger input forces. Nevertheless comparisons are possible using F 0 ¼ 0:625 70:2 N and F 0 ¼ 1:25 7 0:2 N. Fig. 7 shows the comparison of the RMS acceleration data for these two cases: the experimental results are summarised using a boundary polygon, while the numerical simulations are plotted as crosses. It is clear that there is broadly good agreement: most of the data falls within overlapping regions, and both high amplitude peaks and anti-resonant troughs appear at similar frequencies. The most significant discrepancies in both (a) and (b) occur for input frequencies of 80-100 Hz, and is perhaps due to only including modes below 1 kHz in the numerical simulations, while impacts actually excited modes up to around 5 kHz: these high frequency modes may contribute significantly to the RMS acceleration response for this range of input frequencies.
The good agreement obtained between experiment and simulation allows either set of data to be compared with the predicted bounds. Fig. 8 shows the predicted global bounds compared with experimental and numerical results. Fig. 8 (d), and (f). The construction lines for the independent bounds are shown as thinner lines, while the limiting bounds are shown in bold. The construction lines have not been individually labelled to avoid clutter, but they can be identified as follows: bounds that slowly vary with frequency correspond to the passive constraint; lower bounds that take the shape of a typical transfer function curve and upper bounds that are mostly horizontal both correspond to the displacement constraint; and lower bounds that appear over short frequency ranges as an 'n' shape correspond to the compression constraint.
The upper bounds in Fig. 8  There are a number of interesting features of these comparisons. The global upper bounds for RMS displacement shown in Fig. 8(a) and (b) are both limited by the displacement constraint: this is because the passive constraint allows an unbounded displacement as an arbitrarily large static load could be applied without inputting any power to the system. It is clear that these bounds are extremely conservative: for most of the range of input frequencies the bound is nearly two decades larger than the experimental or simulated worst case responses. However, the bound is nearly reached at specific frequencies: 9 Hz in (a) and 9, 32, 67 Hz in (b). It is notable that at these frequencies the bound also shows a small feature, indicating that something interesting occurs here. So even though the bounds are very conservative, they nevertheless contain some useful information regarding the system response. It is also interesting that in (a) and especially in (b) the lower bound is extremely close to the upper bound at 9 Hz, indicating that a high amplitude response is guaranteed at this frequency. This prediction does not depend on detailed knowledge about the particular law associated with the nonlinearity and is true for any passive, displacement-limited compressive nonlinearity. Fig. 8(c) and (d) shows the RMS velocity comparisons. Again the predicted bounds are often two decades larger than actual worst-cases, except at particular input frequencies (e.g. 9 Hz). Fig. 8(e) and (f) shows the RMS acceleration comparisons. The upper bound response for this metric is less conservative: generally one decade larger than the simulated or experimental worst cases. It is interesting that in (e) the experimental RMS accelerations are much more scattered than the numerical simulations predict: this may be due to the high frequency modes above 1 kHz that have been neglected in the model. Nevertheless the upper bound does seem to closely encompass these experimental data points. It is reassuring that both the global lower and upper bounds are not exceeded by any of the numerical simulation data points, even though they are sometimes closely tracked, and that very few experimental data points exceed the bounds. The experimental discrepancies are due to differences in position of the anti-resonances which are sensitive to the modal fitting procedure. Also the bounds are based on including system modes below 1 kHz, but as mentioned earlier the impacting nonlinearity excites modes up to 5 kHz.

Equivalent linear bounds
It was observed in Section 7.4 that key features of the RMS velocity response could be predicted using either the 'freefree' or 'free-pinned' linear response of the beam. This suggests that the RMS velocity is dominated by the frequency contribution at the drive frequency, with other harmonics being less significant. This metric is a strong candidate for applying the 'equivalent linear bounds' as derived in Section 4. Fig. 9(a) shows the linear bounds on RMS velocity compared with both experimental results (shaded polygons) and numerical simulations (crosses) for an input force amplitude of F 0 ¼0.625 N. Fig. 9(b) shows the same comparison using a higher input force of F 0 ¼10 N, but without experimental data as this was beyond the practical range. Note that the linearity assumption only affects the upper bound: the global lower bound only requires finding the best case linear response (see Section 5). The upper bounds in Fig. 9 are based on Eqs. (26) (passive bound) and (30) (displacement bound).
From Fig. 9 it is striking that for both input force amplitudes the linear bound encompasses the ensemble of experimental and numerical data: it is not an overly conservative bound and captures many key features of the data. The level of agreement is particularly apparent for (b) where both upper and lower bounds closely encompass the ensemble of data across the frequency range. It also gives good agreement in (a) despite being somewhat conservative in the range of input frequencies 40-70 Hz.
These 'equivalent-linear' bounds appear to be a very useful tool when the response is dominated by the input driving force. They are less useful for the RMS displacement and RMS acceleration data, which are significantly affected by sub-and super-harmonics respectively. However, to eliminate artificial drift from integrating accelerometer signals the experimental  RMS displacement data was high-pass filtered which removed the very low frequency response. This had the effect that the resulting RMS filtered displacement response was in fact dominated by the input frequency. Fig. 10 shows the comparison between the equivalent linear bounds and the experimental data (but not the simulation data) which shows that the equivalent linear bound is useful for predicting the ensemble of experimental RMS displacement data.
The original global bounds are so conservative because they allow extreme transfer of energy from the input frequency to another frequency (or combination of frequencies). In practice it would be hard to design a system that actually does this, though that is not the purpose of finding global bounds. The modified bounds are so much better because it is much closer to what actually happens in this case: most energy remains at the input frequency and the nonlinearity is encapsulated by the way in which the equivalent linear impedance satisfies energy and displacement constraints.

Estimating confidence bounds using Maximum Entropy
Identifying bounds is useful for quantifying the possible range of responses, and for identifying features of interest such as input frequencies that are guaranteed to result in a large response. The bounds derived above also give insight into the relevant physical processes that govern the response: whether that is the available power that can be transferred to other frequencies, or the gap between end-stops. It is also desirable to estimate something about the probability distribution of responses across an ensemble. This information is not available from constraints alone, and the factors that govern the distribution of outputs is both complicated and unknown. However, it is possible to ask: is there a rational choice of distribution that is only based on the information contained by the constraints without introducing any additional assumptions? One possible answer is to apply the principle of Maximum Entropy [32].
This section makes an initial attempt to see whether Maximum Entropy could yield useful predictions, but a full treatment is beyond the scope of this paper. It is recognised that there are philosophical and practical challenges, but these are beyond the scope of this study and are the subject of ongoing work.
The Maximum Entropy principle for a continuous set of variables seeks to maximise the 'information entropy' H: where p(x) is the probability density function that is sought, R is the region of admissible values of x, and m(x) is a 'measure function' that is required to satisfy invariance of H under a change of variables x-f ðxÞ. If only bounds on x are known then the Maximum Entropy probability distribution is pðxÞ ¼ mðxÞ. This is problematic because it effectively necessitates defining m(x), which is a probability density function that describes 'complete ignorance'. By the principle of indifference [32] it is natural to take m(x) to be a uniform distribution. This is not to claim that it is a correct or unique choice of distribution, only that it is a suitable choice that results from knowledge only of the bounds. The difficulty is that x can be parameterised in different ways, with uniform distributions on those parameters mapping to non-uniform distributions for other parameterisations. In this section a first choice is made for m on pragmatic grounds as a first exploration of whether the general approach has value, and a rigorous analysis is left for future work. It is easiest to apply Maximum Entropy to the equivalent linear force F 2a subject to the constraints on energy, displacement, and compression. These constraints lead to an admissible region of values for the complex variable F 2a , which can also be mapped to an admissible region for Y 2a using Eq. (6) with a given input force F 1a . The complex variable Y 2a is most readily parameterised in terms of magnitude and phase, or real and imaginary parts. The measure function mðY 2a Þ is chosen to be uniform over the admissible set of real and imaginary parts of Y 2a as these are dimensionally consistent (rather than uniform over magnitude and phase). Therefore the Maximum Entropy distribution is a uniform distribution for Y 2a over a region determined by the constraints. The choice of using F 2a or Y 2a is arbitrary as a uniform distribution for Y 2a maps to a uniform distribution for F 2a .
The most straightforward way to carry out calculations is to generate a uniform grid of possible values for real Y 2a f g and imag Y 2a f g, then retain only those points that satisfy the imposed constraints (passivity, displacement, and compression). More efficient schemes could be envisaged but as long as the number of points that remains is 'enough' to obtain an estimate of the probability density function, then this is deemed sufficient for the present study. One might ask why a uniform grid rather than randomised sampling from a uniform distribution: both methods were tested and found to give converging results for a large number of samples, but the uniform grid approach converged with fewer samples (approximately 10 4 samples gave smooth estimates for the 10, 50 and 90 percentiles of the Maximum Entropy distribution). Despite this high number of samples needed, the computation is still highly efficient because it involves calculating the linear response to external input forces. The total computation time for 10 4 samples over 10 3 input frequencies in the range 0-100 Hz was approximately 10 s.
In order to generate an efficient grid of values, bounds for the maximum value of jY 2a j are needed so that the grid resolution is reasonable after applying the constraints. This is achieved by applying the same principles as in Section 4. For passivity: for displacement: and for compression, when real D 22a 4 0 f g : It should be noted that the discrepancies between Maximum Entropy and numerical data in Fig. 11(a) do not invalidate the use of Maximum Entropy. It has been commented by Jaynes [32] that Maximum Entropy gives a logically self-consistent subjective description of uncertainty: it cannot be validated or invalidated by comparison with data obtained by randomised experiments. In particular, Jaynes notes that '…the principle of maximum entropy is most useful to us in just those cases where it fails to predict the correct experimental results.'  The difference between Maximum Entropy and the experimental or simulated data is simply because insufficient data has been used as an input to Maximum Entropy. In other words the experimental and simulated ensembles of data do not span the whole set of uncertainties described by either the bounds ( Fig. 9(a)) or Maximum Entropy ( Fig. 11(a)). What additional information is needed to constrain the predicted bounds or Maximum Entropy solution? Or conversely what additional uncertainty needs to be included in the simulation to 'fill in the gaps'? The most likely candidate is the distribution of nonlinear stiffness values used in the Monte Carlo experiment and numerical simulations: only two materials were used for the end-stops (rubber and polycarbonate) and their effective stiffnesses were used for numerical simulation. A linear (rather than logarithmic) distribution for k c was used so only 10 percent of simulations used values in the range 10 3 o k c o10 5 while 90 percent of simulations were in the range 10 5 ok c o 10 6 . In addition, the rubber end-stop in the experimental tests was only used for the larger clearances. Reconstructing the response from a Maximum Entropy solution shows that high amplitude response in the 40-60 Hz range occurs when 10 3 o k c o10 5 , and only when the clearance is sufficiently small that these are the effective stiffness values for the nonlinearity. Fig. 12 shows the results of a few additional numerical simulations which use a logarithmic distribution for k c in the range 10 3 o k c o10 5 . It can be seen that these results do indeed begin to fill the 'gap'. The new simulations also introduce features in the response above 80 Hz that raise additional questions: however these are left for another study.
It can be seen from Fig. 11 that the Maximum Entropy distribution matches the data extremely well in (b) and less well in (a). This was also true for the equivalent linear bound predictions, but what makes the Maximum Entropy solution interesting is that it gives some indication of the predicted spread of data. Fig. 13 shows the relative standard deviation σ rel (defined as the standard deviation divided by the mean) as a function of equally spaced frequency bands for the numerical simulations, experimental results, and Maximum Entropy distribution for the low amplitude input force corresponding to Fig. 11. Even though the mean response differs from the experimental and numerical data, the predicted variation is still in reasonably good agreement.

Application to an example test structure II: friction-damped beam
The bounds are now applied to a friction-damped structure to test the usefulness of the passive and force-saturating bounds. The results are somewhat similar, so this section presents a deliberately shorter summary. Experimental tests were not carried out with a friction damper, but numerical simulations were carried out using the same suspended 'free-free' beam, in this case coupled at Site (2) to an ensemble of 'uncertain' friction laws.
For the nonlinear time-domain simulations, the nonlinear law was taken to be a continuous approximation to Coulomb's law (based on [33]): where μ s is approximately the static coefficient of friction, μ d is the high-velocity asymptote of the dynamic coefficient of friction, β 1 governs the sharpness of transition when the direction of sliding is reversed, and β 2 governs the sharpness of transition from static to dynamic friction limits. For the purposes of demonstrating the applicability of the bounds to friction-damped systems it is not necessary to know whether or not this law is a good approximation to reality: rather it demonstrates typical features of friction in a numerically convenient form. The equations of motion (Eqs. (56) and (64)) can be written in first-order form and solved using the Matlab function ode45. A numerical Monte Carlo study was carried out, varying the nonlinear contact parameters using a uniform distribution over the following ranges: 0:1 oμ s N 0 o 1, 0:5 oμ d =μ s o 1, 10 À 3 oβ 1 o 10 À 2 ms À 1 , and 1 o β 2 =β 1 o10. The input force was sinusoidal with amplitudes F 1a ¼ 0:625; 1:25; 2:5; 5, and 10 N. It is reassuring that the numerical solutions are fully bounded by the global bounds in Fig. 14(a), albeit the upper bound is again rather conservative at most frequencies. Nevertheless the upper bound is approached near resonant frequencies and both upper and lower bounds give some indication of frequencies of particular interest. In (b) it can be seen that the equivalent linear bound gives an extremely tight bound compared with the simulations, indicating that the RMS velocity is dominated by the same frequency as the input. This turned out to be true for the RMS acceleration and displacement for most simulations (not shown), indicating that the friction-damped system is more weakly nonlinear than the impacting beam. The Maximum Entropy distribution in (c) gives a very good match to the data in terms of both the average response (50 percentile) and the distribution (10 and 90 percentiles). Some details in the response features are not predicted by Maximum Entropy, for example the small peaks near 35 Hz and 70 Hz, but they are adequately encompassed by the predicted variability. This level of agreement is consistent across the input force amplitudes tested for cases with both large and small variability in the response. It is surprising that Maximum Entropy seems to perform so well: it is a subjective measure of uncertainty, so comparisons with ensembles of simulations would not necessarily be expected to yield a close match. The good level of agreement suggests that perhaps the simulations are in some sense 'random enough', and there may be potential for Maximum Entropy to be a useful tool in broader classes of systems with 'sufficiently uncertain' nonlinearities. This is speculative and is left for further investigation.

Conclusions
This paper presents a novel approach for predicting steady-state response metrics for complex multi-degree of freedom systems with a single nonlinearity. The approach considers the nonlinear internal force as an independent excitation, and defines an optimisation (or 'anti-optimisation') problem that seeks the best or worst case with respect to an output metric of interest for the system (e.g. root-mean-square velocity). The 'nonlinear' force is subject to constraints that describe what is known about the nonlinearity (e.g. that it is passive). Semi-analytic upper and lower bounds can be found by careful choice of the output response metric and nonlinearity constraints.
The output metrics considered were the root-mean-square displacement, velocity and acceleration. The constraints defining the nonlinearity were chosen to be a combination of (A) passive; (B) displacement-limited; and (C) forcesaturating. Semi-analytic bounds were derived for each constraint independently, and these bounds can be combined depending on the kind of nonlinearity. Under some conditions bounds could be derived for nonlinearities that were constrained to be compressive (i.e. force opposing displacement).
The bounds were compared with results from experimental and numerical tests. The experimental test rig consisted of a suspended beam, driven sinusoidally at one end with a double end-stop nonlinearity at the other. A Monte Carlo experiment was carried out by varying the nonlinear contact properties, generating an ensemble of approximately 21,000 tests. A numerical model was developed based on the same underlying linear system and with similarly variable contact properties. An equivalent numerical Monte Carlo experiment was also carried out, and the steady-state metrics were found to be in very good agreement with the experimental results.
The bounds associated with a passive and displacement-limited nonlinearity were compared with the impacting beam experimental and numerical results. The global upper and lower bounds encompassed all numerical data points and nearly all experimental data points. The global bounds were found to be generally very conservative compared with results from the specific impacting beam test system. Nevertheless they could still be useful, revealing frequencies of high amplitude response or identifying frequencies of interest.
It was found that the RMS velocity metric from the impacting beam results appeared to be dominated by the input excitation frequency. This led to the development of 'equivalent linear bounds': constraining the nonlinear force to be at the same frequency as the input frequency. The linear bounds were much less conservative than the global bounds, and gave excellent agreement with the steady-state RMS velocity data from experiment and simulations. This approach is directly analogous to the 'Describing Function' approach (or Harmonic Balance Method with just the fundamental frequency retained) and would be expected to be valid under similar conditions, i.e. when the response metric is dominated by the excitation frequency.
A numerical experiment was carried out using the same underlying linear system (the suspended beam), but with a frictional nonlinearity. Monte Carlo simulation results were compared with bounds corresponding to the passive and forcesaturating constraints. Again the global bounds were in agreement but were generally rather conservative, and the equivalent linear bounds gave excellent agreement.
The principle of Maximum Entropy was used for both test systems to predict the response distribution associated with the equivalent linear bounds. For the impacting beam the predicted average and variability in the response gave good agreement with experimental and numerical results, particularly for high input force amplitudes. At low input force amplitudes, the Maximum Entropy solution differed significantly from the data: this revealed that the ensemble of nonlinearities encompassed by Maximum Entropy was much broader than encompassed by the numerical and experimental tests, demonstrating that care is needed when interpreting Maximum Entropy predictions. The Maximum Entropy distribution for the friction-damped system gave excellent agreement across all metrics considered and for all input force amplitudes simulated.
The results obtained show that the simple constraints and metrics considered in this paper can yield useful predictions of the response of a complex system with a single, uncertain and harsh nonlinearity. Generalisation to systems with multiple nonlinearities is left for future work.

Data access
Additional data related to this publication is available at the University of Cambridge data repository: http://dx.doi.org/ 10.17863/CAM.679 response U 2 is given by: Writing real and imaginary parts separately, let F 2 ¼ F R þiF I , so that: This is maximised when giving and giving Combining these such that F 2 ¼ F R þiF I gives: To show that this is equivalent to an impedance match, let F 2 ¼ ÀiωU 2 =ðiωD 22 Þ Ã . This gives: Rearranging for F 2 gives: which is the same result as Eq. (80).