A two-state hysteresis model from high-dimensional friction

In prior work (Biswas & Chatterjee 2014 Proc. R. Soc. A 470, 20130817 (doi:10.1098/rspa.2013.0817)), we developed a six-state hysteresis model from a high-dimensional frictional system. Here, we use a more intuitively appealing frictional system that resembles one studied earlier by Iwan. The basis functions now have simple analytical description. The number of states required decreases further, from six to the theoretical minimum of two. The number of fitted parameters is reduced by an order of magnitude, to just six. An explicit and faster numerical solution method is developed. Parameter fitting to match different specified hysteresis loops is demonstrated. In summary, a new two-state model of hysteresis is presented that is ready for practical implementation. Essential Matlab code is provided.


Introduction
In this paper, we follow up on our recent work on lowdimensional modelling of frictional hysteresis [1]. Contributions of this paper include a different underlying frictional model with greater intuitive appeal, new analytical insights, reduction in the number of states from six to two, 1 reduction in the number of free parameters by an order of magnitude, and demonstration of fitting these parameters to several hypothetical hysteresis loops. The net result is a two-state hysteresis model that captures minor loops under small reversals within larger load paths and is ready for practical numerical implementation (simple Matlab code is provided).
For elementary background, we note that hysteresis is a largely rate-independent, irreversible phenomenon that occurs in many systems. Much research on hysteresis has been done over several decades (e.g. three volumes of Bertotti & Mayergoyz [2] and references therein). For classical papers, see, for example, Ewing [3], Rowett [4], Preisach [5], Jiles & Atherton [6]. For our 1 Two is the theoretical minimum number. Single-state models cannot capture commonly observed behaviour (see fig. 3 in [1]). present purposes, for hysteresis in mechanical systems with elastic storage and frictional dissipation, a model due to Iwan [7,8] seems promising, but is high-dimensional and deeply nonlinear with several dry friction elements. By contrast, the famous Bouc-Wen model ( [9,10]; see also [11]) is one-dimensional but fails to form minor loops under small reversals within larger load paths.
With this background, we recently studied [1] a frictional hysteretic system given by μ sgn(ẋ) + Kx = bf (t), (1.1) where x is high-dimensional; μ is diagonal; K is symmetric and positive definite; b is a column matrix; f (t) is scalar and differentiable; and the signum function 'sgn' is defined elementwise as follows: Equation (1.1) can be solved incrementally via a linear complementarity problem (LCP) [12] or, less efficiently, using other means as described later. The solution of equation (1.1) captures important aspects of hysteresis including formation of minor loops. From equation (1.1), we had developed a reduced order model with six states. The order reduction included finding the slip direction as a minimizer of a complicated function containing many signum nonlinearities, for which a convenient analytical approximation was found. However, some shortcomings remained. The choice of basis vectors involved some arbitrariness whose implications were unclear; reductions below sixth order gave poor results; and there were too many fitted parameters for practical use.
In the light of the above, this paper makes the following notable progress. A more intuitively appealing frictional system is studied here, motivated by the Iwan model [7] and yielding a simpler governing equation. The numerically obtained basis vectors are now amenable to analytical approximation, providing better analytical insight. Finally, a two-state, reduced order model is derived that allows practical parameter fitting to match a range of given data.
As far as we know, the two-state model developed here has no parallel in the literature.

New frictional system
Differing somewhat from Biswas & Chatterjee [1], here we consider the intuitively simpler highdimensional frictional system sketched in figure 1. In this n-dimensional model (with n large), each spring has stiffness 1/n, and friction coefficients at the slip sites are As indicated in the figure, u(t) is a displacement input to the system, for which a force f (t) is needed. Friction forces at the slip sites are written as where the overdot denotes a time derivative and the signum function is understood to be multivalued at zero (taking any necessary value between ±1). The governing equation is In matrix form which resembles equation (1.1) but is in fact simpler because the elements of μ now have a regular variation (they are linearly increasing), K is a scalar multiple of the identity matrix, and all elements of b are identical. The output force f (t) is x n k 1 = 1 n k 2 = 1 n k n -1 = 1 n k n = 1 n Figure 1. A high-dimensional frictional system.
Incidentally, if a spring of stiffness k s is attached to the system, in parallel, being stretched by u(t), then the net output force is We will use the parameter k s later for better fitting of the model to specified hysteretic response curves; here, we note that k s in equation (2.3) has no effect on the solution of equation (2.1), which takes u(t) as its input. We solve equation (2.1) incrementally by casting it first into an LCP (as described in Biswas & Chatterjee [1]) and then by using Lemke's algorithm (as implemented by Miranda & Fackler [13]). There is in fact a large literature on solving friction problems using the LCP; readers interested in the theory may consult, e.g. Klarbring & Pang [14].
Alternatively, equation (2.1) may be regularized as follows: In the above the exponential, the absolute value within it, and the signum function are all evaluated elementwise; the fact that K is a scalar multiple of the identity and that μ is diagonal has been used to simplify the first term; and is a regularizing parameter. The justification for this regularizing method is that (i) the exponential term produces high rates of change only when the concerned absolute value exceeds unity, and (ii) the signum term outside guides that rate of change in the correct direction. Further discussion of this regularizing method is avoided to minimize the distraction from the main flow of the paper, but a numerical example is given in appendix A. Note that equation (2.4) may be solved using an ordinary differential equation (ODE) solver. For our numerical solution, the arbitrarily selected numerical parameters are as follows. We use n = 500, μ 0 = 0.002 and the two-frequency displacement input u(t) = 0.6748 sin(t) + 0.2887 sin(6.5581t). (2.5) We then solve equation (2.1) incrementally using 12 000 uniform time steps of dt = 0.001 each.
In this way, we obtain a 12 000 × 500 matrix, wherein each row is ξ T at some instant. From ξ , we find f (t) using equation (2.2). Figure 2 shows f (t) against u(t). Minor loops are seen. As emphasized by Biswas & Chatterjee [1], such minor loops are not predicted by the Bouc-Wen model or indeed any hysteresis model with a single state.
We now develop a reduced-order model from this high-dimensional hysteretic system (figure 1, equation (2.1)).

Reduced-order model
Our system is discrete. However, the largeness of n and the slow variation in μ j suggest that we might loosely think of an approximating continuous system. With this motivation, we assume In the above, m is the reduced dimension, q r (t)'s are functions of time and the φ r (x)'s are basis functions yet to be chosen.

Choice of basis functions
The singular value decomposition 2 of the 12 000 × 500 data matrix from §2 shows that the first two singular values are distinctly larger than the rest (figure 3). Figure 4a shows the first three singular vectors plotted against x 3/2 (where x = 1, 2, . . . , 500) for two different solutions using different μ 0 's and u(t)'s. Figure 4b shows that, after rescaling, the singular vectors for the two cases are similar. These observations suggest that the basis functions may be taken as functions of x 3/2 (the 3 2 power is empirical, based on the fact that the slope near x 3/2 is finite and non-zero). In order to ensure eventual decay to zero, we choose the following basis functions for lower order modelling: The free parameter α > 0 controls the decay rate of the basis functions. The actual discrete versions of these basis functions will be orthonormalized below for analytical convenience. It may be noted that the ad hoc form of equation (3.2) leads to simplification below, but also means that the possibility of a very good match with the full numerical solution has now been abandoned.

Slip criterion
Our slip criterion is this: slip cannot occur if the accompanying frictional dissipation exceeds the external work input minus the internal increase in potential energy. This criterion will help us find slip directions and rates below.
What follows in §3.2-3.4 has much in common with Biswas & Chatterjee [1], but is included for completeness. Let ξ ≈ Φq 2 Also known as the proper orthogonal decomposition (e.g. [15]).     The rate of increase in spring potential energy is d dt and the rate of work done by the external force f iṡ Substituting the above into our slip criterion yields with matrices μ, K and b as described in §2. We define and We note thatK = (1/n)I above, with I being the identity matrix; and choosing n = 500, b = 0.0326 0.0068 and 0.0285 0.0059 for α = 0.0008 and α = 0.0012, respectively. These values will be used later.
If the minimum possible value of the left-hand side of inequality (3.9) is negative, rapid slip will occur because there is no inertia; if the minimum is positive, no slip can occur; and if it is zero over a time interval, sustained slip can occur at a finite rate.
Accordingly, we will minimize G(η) + η T c at each time step with respect to η, and see if the minimum is negative or positive. Noting that G(η) + η T c is homogeneous of degree one in η, our minimization will be done subject to η T η = 1. The only difficulty is that G(η) is a complicated function. Luckily, a convenient analytical approximation for G(η) of equation (3.7) was found by Biswas & Chatterjee [1].

Approximation of G(η)
G(η) is homogeneous of degree one in η. In Biswas & Chatterjee [1], a similar function was encountered and the following approximation was considered: with A a fitted symmetric and positive definite matrix; also, β = 0.5 was found to be near-optimal and selected due to analytical convenience. We use the same approximation here (again with β = 0.5).
Fitting of the matrix A was described in detail in the earlier paper. Here we fix n = 500, let μ 0 vary, and fit A for α = 0.0008 and α = 0.0012. We find that to an excellent approximation

Slip direction
The slip direction η minimizes y = G(η) + η T c for given Introducing a Lagrange multiplier for the constraint, using the approximation G(η) ≈ η T Aη, taking the gradient with respect to η, and letting 2λ =λ, we obtain The minimizing 2 × 1 matrix η is found by solving a 4 × 4 eigenvalue problem (see appendix B). If the corresponding then slip occurs in the direction of η.

Reduced-order model using incremental map
The unit vector η (computed as outlined in appendix B) gives the direction of slip, but the actual rate of slipq remains to be found. In Biswas & Chatterjee [1], we used a stiff system of ODEs with a large ad hoc gain parameter. Here, we use a faster explicit incremental formulation as follows.
As mentioned earlier, a direct comparison of these results with the earlier high-dimensional simulation is not meaningful because we have adopted analytical expressions for the basis functions (equation (3.2)) instead of numerically computed proper orthogonal modes from actual solutions ( figure 4). The high-dimensional model has thus motivated the structure of the lower dimensional model, and now we work directly with the latter.

Final reduced-order model using a differential equation
Although the numerical results that follow were obtained using the incremental map given above, some users may prefer to have a differential equation for the hysteresis model. We now present one. We also present the entire hysteresis model compactly and algorithmically below. A detailed example calculation is given in appendix C. Since A can be diagonalized by an orthogonal coordinate transformation, for compactness we henceforth assume that A is diagonal, with elements increasing order. In the m = 2 case, we introduce a scalar factorμ > 0 to write A =μ σ 0 0 1 , 0< σ < 1.
(ii) A matrixK, which is a scalar multiple of the identity. Herē 2. System inputs: u is the system input, andu is known at each instant. 3. The state vector: q is an m × 1 state vector (here m = 2). Given the above system matrices and inputs, we first compute c =Kq −bu. Subsequently, the possible slip direction η is computed as a function of A and c as described in appendix B, using straightforward matrix calculations of order 2m.
Using the above, we define the intermediate quantity 4  Figure 6. Fitting a hysteresis curve. Dotted points represent given data. Vertical distances between data and fitted curve are squared and added. That sum is minimized with respect to the fitted parameters.
The above system of ODEs does not need a large ad hoc 'gain' parameter as was used in Biswas & Chatterjee [1].

Fitting parameters to given data
As described in §3.5, for a two-state model we have five fitted parameters (μ, σ ,ᾱ,b 1 andb 2 ). If we used a three-state model, we would add one diagonal entry in A, no new parameters toK, and one element tō b, obtaining a model with seven fitted parameters.
Introduction of an added spring in parallel with constant k s , as in equation (2.3), would make it six fitted parameters for the two-state model and eight fitted parameters for the three-state model. For clarifying that a spring in parallel is implied, we will refer to these latter two as 5 + 1 and 7 + 1, respectively.
We now fit some hypothetical hysteresis loops using our two state, 5 + 1 parameter model, using nonlinear least squares as depicted schematically in figure 6 (the minimization was done using Matlab's built-in function fminsearch).

Results and discussion
We now show results of fitting the hysteresis model to some arbitrary input data. Results are depicted graphically here; numerical values of fitted parameters are given in appendix D. Figure 7 shows fitting of three hysteresis loops, numbered 1 to 3. Each case is depicted in three parts, namely (a), (b) and (c). Parts 1(a) to 3(a) show the prescribed or desired loop shapes (half the cycle). Parts 1(b) to 3(b) show the corresponding hysteresis loops obtained by fitting parameters. These fitted parameters are then used to plot hysteresis loops using smaller (85 and 70%) input amplitudes, and parts 1(c) to 3(c) show these loops corresponding to 100% (blue solid), 85% (black dotted) and 70% (red dashed) amplitude.
The model fits the above given data (figure 7) well. However, the model does not work well for hysteresis curves with two distinct changes of slope. For example, figure 8 shows how a hysteresis curve made of three straight lines is not captured very accurately by the two-state model (or even, in attempts not documented here, by models with three or four states). However, overall, our model fits a reasonable range of data usefully well.
In figure 9, we consider fitting of minor loops. In case 1(a), we specify two minor loops within a major loop. The corresponding fit is shown in 1(b). In case 2(a), we try to thicken one of the minor loops of case 1(a). The model captures the loop thickening in the third quadrant of the plot, but similar changes occur in the first quadrant as well, because our hysteresis loops are symmetrical about the origin.
The model developed in this paper, as explained and demonstrated above, has several advantages over our earlier work [1]. These advantages include a more intuitive underlying frictional system, analytical insights into basis functions, a minimal number of states (two), a small number of fitted parameters, and the ability to match a reasonable range of hysteretic behaviours.
The computational complexity of our model exceeds that of the Bouc-Wen model but compensates by capturing minor loops. Direct computation with the Iwan model for arbitrary forcing, in high dimensions, would be significantly more complex than for our model.  Further study may clarify the precise advantages of including more than two states in the hysteresis model; why (apparently) two distinct changes in slope are difficult for the model to capture; and how parameter fitting can be done more reliably and efficiently than using general purpose minimization routines with random initial guesses.
New research questions might also now be addressed somewhat more easily; these include control strategies for such hysteretic systems, as well as the nonlinear dynamics of systems that include elements with such hysteretic behaviour.

Closing note
An anonymous reviewer of this article has brought to our attention the recent work by Scerrato et al. [16,17]. These works are interesting and complementary to our approach, possibly opening up new lines of research for both.
(2) In particular, in these works, a micro-structural model for dissipation in concrete is considered, while our starting model is abstract although physical, and not motivated by any specific material. Both models allow complex stress histories, but their work considers multiaxial stress states while ours so far does not. In that work, the experimental hysteresis loops fitted are asymmetric, while ours are symmetric. We have fewer state variables than they do.
Thus, their work motivates our approach to try to incorporate triaxiality in stress and asymmetry in hysteresis loops, while our work suggests new methods of approximation and model-order reduction that may lead to improvements in their approach. an orthogonal transformation we can take A to be diagonal: We reproduce equation (3.12) below: From equation (B 2) and Since η T Aη =μ(η 2 1 σ + η 2 2 ), by equations (B 3) and (B 4), we havē Now, somewhat fortuitously, we consider the 4 × 4 matrix It turns out thatλ of equation (B 5) is an eigenvalue of the above B. Setting We now use a matrix determinant lemma (e.g. [18]) which says that for a square matrix H and column matrices g and h of appropriate sizes, By this lemma, equation (B 7) gives which, upon multiplying the terms out, gives equation (B 5). Thus,λ is an eigenvalue of B. Now consider the corresponding eigenvectors Since ψ is part of an eigenvector it can be scaled such that ψ T ψ = 1, and the above equation remains unaffected. However, ψ remains indeterminate up to a '±' sign. Comparing equations (B 15) and (B 1), we find η = ±ψ. Thus, the algorithm for computing η is as follows. First, we construct matrix B as in equation (B 6). We find its eigenvaluesλ. For every realλ, we find the corresponding portion of its eigenvector, ψ, normalized to ψ T ψ = 1, and with sign chosen such that ψ T c ≤ 0. Searching through all such ψ (numbering at least 2 as shown by Biswas & Chatterjee [1], and at most 2m which is the size of B), we select the one that minimizes ψ T Aψ + ψ T c. That ψ is the slip direction η.

B.2. Matlab code for computing η
The Matlab code below finds the slip direction η given a diagonal A and 2 × 1 vector c.  Appendix C. A two-mass system with a hysteretic damper Figure 11 shows a two degree of freedom oscillator. There are two unit masses attached by springs of unit stiffness as shown in the figure. A two-state hysteretic damper is attached to the first mass. The displacements of the masses are taken as u 1 and u 2 , with u 1 being the input to the hysteretic damper.
We arbitrarily choose the parameter values

Equations of motion of the two masses arë
where f is the hysteretic damper force. In addition, there are two first-order differential equations for the two-dimensional state q, as given earlier in §3. 5. The hysteretic force f is calculated at each instant following equation (3.18): f = u 1 − q Tb . Figure 12a shows the system response for an arbitrary initial condition. The blue solid line is displacement u 1 , the black dashed line is the displacement u 2 and the red dotted line is the hysteretic force f. For this two degree of freedom system, the response has more than one frequency. Consequently, u 1 shows short reversals within larger oscillations. Consequently, the hysteretic force shows minor reversals (figure 12b).