# On the inversion of describing functions in nonlinear system analysis

## Huan Du1, Zhiqiang Li2, Kjell Ahlin3

1, 2School of Reliability and System Engineering, Beihang University, Beijing, China

3Xielalin Consulting, Akersberga, Sweden

2Corresponding author

Vibroengineering PROCEDIA, Vol. 12, 2017, p. 172-177. https://doi.org/10.21595/vp.2017.18538
Received 25 April 2017; accepted 29 April 2017; published 30 June 2017

Copyright © 2017 JVE International Ltd. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Views 36
Abstract.

Nonlinear mechanical vibrations are widely found in engineering structures. Traditional linear analysis methods are failing to meet the growing requirements of analysis accuracy. Because of the complexity of nonlinear system, there is no general method for nonlinear system analysis. Based on the describing function method, a polynomial inversion method may be used in nonlinear system identification. Using this method, the measured DF can be inverted into the form of force which has a clearer physical meaning and is easier to use for parameter fit. A MDOF system with two types of nonlinearities at the same location is studied to demonstrate this method.

Keywords: nonlinear system, identification, describing function, polynomial inversion method.

#### 1. Introduction

In the past, most attention was paid to linear systems, and theories and methods have been developed for linear system modeling, simulation and prediction. However, there exist lots of nonlinearities in an actual engineering system and some of them may have a significant effect in system dynamic characteristics. The study of nonlinearities is indispensable and inevitable. However, one should be forced to admit that there are no general analysis methods that can be applied to all nonlinear system in all instances .

During the past decades, lots of work has been done in identification of system nonlinearities and they can be roughly divided into two categories, time domain methods and frequency domain methods. Time domain identification of nonlinear structural models relies on processing time series while the data processing in frequency domain methods can take form of Fourier spectra, frequency response and transmissibility functions, or power spectral densities . The describing function, DF, based on the definitions of the frequency response function and the first-order harmonic balance method, was used in analyzing nonlinear structures by Tanrikulu O., et al. . A further research was made by Özer and Özgüven  in calculating the describing function’s approximate solutions by using the Sherman-Morrison Method. Aykan and et al  proposed a method to transform the describing function into the restoring force to make it more visualized.

Finding the nonlinear forces corresponding to the measured describing functions is called describing function inversion . In this paper, nonlinear system identification using describing function is introduced briefly. A polynomial inversion method of inverting the measured describing function is proposed. At last, an illustrative example is given to demonstrate this method.

#### 2.1. Describing function method

The use of non-linear describing function matrix to represent the nonlinear forces in a MDOF system was proposed and developed by Özer et al. . They established a set of feasible methods for nonlinear system identification, including early detection, nonlinear localization and parameter estimation. This method was confirmed to have an acceptable application to a MDOF system with several nonlinearities at different locations in a case study . As this method was discussed in detail in ref. , a brief introduction is given below for readers’ comprehension.

The equation of motion for a nonlinear MDOF system under excitation can be written as:

(1)
$\left[M\right]\bullet \left\{\stackrel{¨}{x}\right\}+\left[C\right]\bullet \left\{\stackrel{˙}{x}\right\}+\left[K\right]\bullet \left\{x\right\}+\left\{nlf\left(x,\stackrel{˙}{x}\right)\right\}=\left\{f\right\},$

where $\left[M\right]$, $\left[C\right]$, and $\left[K\right]$ are mass, viscous damping and stiffness matrices of the linear system, respectively. $\left\{f\right\}$ and $\left\{x\right\}$ stand for the external forces and the responses of the system. $\left\{nlf\left(x,\stackrel{˙}{x}\right)\right\}$ represents the nonlinear internal forces in the system, and it is a function of the displacement and/or velocity responses, depending on the type of nonlinearity present in the system.

A harmonic excitation with circular frequency $\omega$ may be written:

(2)
$\left\{f\right\}=\left\{F\right\}{e}^{j\omega t}.$

Using first order harmonic balance method, the response is assumed to consist of only the first fundamental harmonic, disregarding higher order harmonics. Since the response is harmonic as the assume, $\left\{x\right\}=\left\{X\right\}{e}^{j\omega t}$, the internal nonlinear forces can be harmonic too:

(3)
$\left\{nlf\left(x,\stackrel{˙}{x}\right)\right\}=\left[∆\left(x,\stackrel{˙}{x}\right)\right]\left\{X\right\}{e}^{j\omega t},$

where $\left[∆\left(x,\stackrel{˙}{x}\right)\right]$ is the ‘non-linearity matrix’ that is dependent on response amplitude. The “non-linearity matrix” can be formed by using describing function $v$ representations of internal non-linear forces as follows:

(4)
(5)
${∆}_{rj}={-v}_{rj}.$

The relationship between describing functions and nonlinear forces will be introduced in Section 2.2. From the equations above:

(6)
$\left(-{\omega }^{2}\bullet \left[M\right]+j\bullet \omega \bullet \left[C\right]+\left[K\right]\right)\bullet \left\{X\right\}+\left[∆\right]\bullet \left\{X\right\}=\left\{F\right\}.$

Then the nonlinear frequency response function (FRF) matrix $\left[{H}^{N}\right]$ can be calculated as:

(7)
$\left[{H}^{N}\right]={\left(-{\omega }^{2}\bullet \left[M\right]+j\bullet \omega \bullet \left[C\right]+\left[K\right]+\left[∆\right]\right)}^{-1}.$

The FRF of underlying linear system, $\left[{H}^{L}\right]$, can be written as:

(8)
$\left[{H}^{L}\right]={\left(-{\omega }^{2}\bullet \left[M\right]+j\bullet \omega \bullet \left[C\right]+\left[K\right]\right)}^{-1}.$

Notice that the nonlinear FRF may be different at different excitation levels due to the non-linearity matrix. The non-linearity matrix can formally be obtained by:

(9)
$\left[∆\right]={\left[{H}^{N}\right]}^{-1}-{\left[{H}^{L}\right]}^{-1}.$

Eq. (9) is the core of this method. Once we get the non-linearity matrix, we have the ability to localize and identify the nonlinearities in a MODF system. Although there still exist some challenges, such as the calculation of $\left[{H}^{L}\right]$ and $\left[{H}^{N}\right]$, they were discussed and partially solved in previous studies mentioned above. Readers who are interested can have a more comprehensive understanding in ref. [4, 7].

#### 2.2. Polynomial inversion method

By the use of DF method, the nonlinear force describing function of the response’s amplitude can be obtained. Generally, they are represented in the form of two related series, the response amplitude, $\left\{dfx\right\}$, and corresponding DF values, $\left\{df\right\}$. The describing function, which can be used for determining the type of nonlinearity as well as for parametric identification, has a quite complex form and it will be worse if there exist multiple nonlinearities in one location of the nonlinear system. This not only makes it necessary for a wealth of engineering experience but also greatly increases the difficulty for identification. An inversion method that works for general describing functions will effectively solve this problem. In the rest of this section, a polynomial inversion method will be given in detail.

The describing function for general dynamic system nonlinearity can be written as:

(10)
$v\left(x,\stackrel{˙}{x}\right)=\frac{j}{\pi A}{\int }_{0}^{2\pi }nlf\left(x,\stackrel{˙}{x}\right){e}^{-j\tau }d\tau ,$

where the response $x$ is supposed to be:

(11)
$x\left(t\right)=A\mathrm{s}\mathrm{i}\mathrm{n}\left(\omega t+\theta \right)=A\mathrm{s}\mathrm{i}\mathrm{n}\left(\tau \right),$
(12)
$\tau =\omega t+\theta .$

Eq. (10) indicates that the describing function is a complex function, and a complex function can be divided into two parts, real part and imaginary part:

(13)
$v\left(x,\stackrel{˙}{x}\right)={v}_{p}\left(x,\stackrel{˙}{x}\right)+j{\bullet v}_{q}\left(x,\stackrel{˙}{x}\right).$

The real part ${v}_{p}$ corresponds to stiffness nonlinearities while the imaginary part $j{\bullet v}_{q}$ corresponds to damping nonlinearities . For each part, the derivation for the inverse is given as follows:

An arbitrary nonlinearity can be expressed as a sum of n nonlinearities, ${N}_{1}$, ${N}_{2}$,..., ${N}_{n}$:

(14)
$nlf\left(x\right)={nlf}_{1}\left(x\right)+{nlf}_{2}\left(x\right)+\cdots +{nlf}_{n}\left(x\right)={\sum }_{i=1}^{n}{nlf}_{i}\left(x\right).$

The DF for composite nonlinearities is given by:

(15)
$v\left(A\right)=\frac{j}{\pi A}{\int }_{0}^{2\pi }{\sum }_{i=1}^{n}{nlf}_{i}\left(A\mathrm{s}\mathrm{i}\mathrm{n}\tau \right){e}^{-j\tau }d\tau =\frac{j}{\pi A}{\sum }_{i=1}^{n}{\int }_{0}^{2\pi }{nlf}_{i}\left(A\mathrm{s}\mathrm{i}\mathrm{n}\tau \right){e}^{-j\tau }d\tau ={\sum }_{i=1}^{n}{v}_{i}\left(A\right).$

Eq. (14) and Eq. (15) indicate that the DF and the inverse DF transformations satisfy the superposition principle. By polynomial fitting, an arbitrary function can be approximately represented in the form of a sum of polynomials. Theoretically, the higher the polynomial order we use, the more accurate result we will get.

Because the DF yields zero when the corresponding force is even, we represent the input-output characteristic for a nonlinearity as an odd nth-order polynomial:

(16)
$nlf\left(x\right)\approx {c}_{1}x{+c}_{2}x\left|x\right|+\cdots +{c}_{n-1}{x}^{n-2}\left|x\right|+{c}_{n}{x}^{n}.$

The $m$th item of DF can be computed by integrating over the interval $0<\tau <\pi /2$, as:

(17)
${v}_{m}\left(A\right)=\frac{4}{\pi A}{\int }_{0}^{\frac{\pi }{2}}{nlf}_{m}\left(A\mathrm{s}\mathrm{i}\mathrm{n}\tau \right)\mathrm{s}\mathrm{i}\mathrm{n}\tau d\tau =\frac{4}{\pi A}{\int }_{0}^{\frac{\pi }{2}}{c}_{m}{\left(A\mathrm{s}\mathrm{i}\mathrm{n}\tau \right)}^{m}\mathrm{s}\mathrm{i}\mathrm{n}\tau d\tau$

Thus, the DF is written as:

(18)
$v\left(A\right)={c}_{1}+{\frac{8}{3\pi }c}_{2}A+\frac{3}{4}{c}_{3}{A}^{2}+\frac{32}{15\pi }{c}_{4}{A}^{3}+\frac{5}{8}{c}_{5}{A}^{4}+\cdots .$

The polynomial inversion method is composed by the following major steps: 1) Get the approximate describing function by the polynomial fitting of DF values calculated; 2) Transform the polynomial coefficients by Eq. (17) respectively and get the approximate nonlinear force function (NLFa); 3) Plot the NLFa, make a guess of the type of nonlinearity and give an initial value of its parameters from the graphic information (taking the saturation for example, $NLF=nlf\left(k,M\right)$, $k$ is the gain of linear part and $M$ is the saturation value).

In order to improve the accuracy of identification, the $fm$ in search function of MATLAB is used. Calculate the DF from the NLF we get in Step 3. Comparing this DF with the original DF we can finally get the parameters’ values with accuracy through iterations. An illustrative example is given in the next section.

#### 3. An illustrative example

A four DOFs system with two nonlinear elements in the same DOF is studied here.

Fig. 1. 4 DOFs nonlinear system used in case System parameters:

$M=$ [1, 0, 0, 0; 0, 2, 0, 0; 0, 0, 3, 0; 0, 0, 0, 5];

$C=$ [10, –5, 0, 0; –5, 10, –5, 0; 0, –5, 10, –5; 0, 0, –5, 10];

$K=$ [1000, –500, 0, 0; –500, 1000, –500, 0; 0, –500, 1000, –500; 0, 0, –500, 1000].

The nonlinear elements between m1 and ground are ‘cubic stiffness’, ${nlf}_{1}$ and ‘gap with end stiffness’ ${nlf}_{2}$:

${nlf}_{1}=1.7e5×{{x}_{1}}^{3},$
${nlf}_{2}=\left\{\begin{array}{ll}0,& \left|{x}_{1}\right|<0.008,\\ 2.3e2×\left(\left|{x}_{1}\right|-0.008\right)\bullet \mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left({x}_{1}\right),& \left|{x}_{1}\right|\ge 0.008.\end{array}\right\$

The sinusoidal excitation is applied at DOF1. The system is linear under small excitation, while the nonlinearities appear when the excitation gets big enough.

The simulation tests are performed using harmonic force with frequency 2 Hz, Fig. 2. The FRFs obtained from 0.1 N amplitude harmonic force are taken as linear FRFs of system. For the nonlinear FRFs the force amplitude is stepped up to10 N. The simulation is performed with functions in a nonlinear MATLAB toolbox .

Fig. 2. Linear and Nonlinear FRF ${H}_{11}$ Fig. 3. Sum of NLN of different DOFs Fig. 4. Comparison of calculated and fitted DF Fig. 5. Comparison of theoretical and fitted $nlf$ From Fig. 3 it is found that the nonlinearities only exist between DOF1 and ground. By using measured FRF the describing function of system nonlinearities can be calculated by Eq. (9) and is shown in Fig. 4 (black). By observing, a turning point can obviously be found at the displacement 0.008 [m]. Polynomial fitting in interval [0, 0.008] and [0.008, $+\infty$] is made with degrees 3 and 5, respectively. The polynomial fitted DF curve is obtained and shown in Fig. 4 (red). For the convenience of display, all values of displacement are multiplied by 1e+3. The polynomial coefficients, , for each part are:

Dividing ${coef}_{df1}$ and ${coef}_{df2}$ by the corresponding transform coefficients from Eq. (17), the polynomial coefficients for nonlinear force expressions are found. The comparison of fitted (red) and theoretical (black) nonlinear force functions is shown in Fig. 5:

By observing Fig. 5, when $0, the expression of force should be in the form of $nlf\left(x\right)={a}_{3}×{x}^{3}$ and the corresponding theoretical describing function, with parameter ${a}_{3}$, is $df\left(A\right)={a}_{3}×{0.75×A}^{2}$. Setting the difference between DF with ${a}_{3}$ and DF calculated as the FUN in $fm$ in search (FUN, ${a}_{3}$) of MATLAB and giving ${a}_{3}$ an initial value 1.706$e5$, we can get the value of ${a}_{3}$ with accuracy, ${a}_{3}=$ 169687. The change of parameter ${a}_{3}$ with the number of iterations is shown in Fig. 6. As it’s a pure cubic stiffness in the first interval, the nonlinearities could be a combination of a cubic stiffness and a nonlinearity with gap. Compare the values of fitted $nlf$ in Fig. 5 and $169687×{x}^{3}$$\left(x\ge 0.008\right)$. The guess for the nonlinear force is $nlf\left(x\right)=169687×{x}^{3}+k\left(x-0.008\right)$. Similar to parameter ${a}_{3}$, we do the parameter fitting using fminsearch. Give an initial value $1e3\bullet 2.0955$ and parameter $k$ stays at the value $k=$ 228.0690 after iterations.

Fig. 6. Iteration of parameter ${a}_{3}$ Fig. 7. Iteration of parameter $k$ #### 4. Conclusions

Identification of nonlinear system is considered by the inversion of describing functions. The method in this paper reduces the difficulty but improves the accuracy of nonlinear system identification. Thanks to the broad applicability of polynomial fitting, the polynomial inversion method should work out well for most situations in which the nonlinearities are in-dependent of frequency. The result of the illustrative example shows that this method is capable to identify the nonlinear systems with acceptable accuracy.

1. Kerschen G., et al. Review: past, present and future of nonlinear system identification in structural dynamics. Mechanical Systems and Signal Processing, Vol. 20, Issue 3, 2006, p. 505-593. [Search CrossRef]
2. Noël J. P., Kerschen G. Nonlinear system identification in structural dynamics: 10 more years of progress. Mechanical Systems and Signal Processing, Vol. 83, 2017, p. 2-35. [Search CrossRef]
3. Tanrikulu O., et al. Forced harmonic response analysis of nonlinear structures using describing functions. AIAA Journal, Vol. 31, Issue 7, 1993, p. 1313-1320. [Search CrossRef]
4. Özer M. B., Özgüven H. N., Royston T. J. Identification of structural nonlinearities using describing functions and Sherman-Morrison method. Mechanical Systems and Signal Processing, Vol. 23, Issue 1, 2009, p. 30-44. [Search CrossRef]
5. Aykan M., Özgüven H. N. Parametric identification of nonlinearity in structural systems using describing function inversion. Mechanical Systems and Signal Processing, Vol. 40, Issue 1, 2013, p. 356-376. [Search CrossRef]
6. Jalali H., Bonab B. T., Ahmadian H. Identification of weakly nonlinear systems using describing function inversion. Experimental Mechanics, Vol. 51, Issue 5, 2011, p. 739-747. [Search CrossRef]
7. Zhongge Zhao, Chuanri Li, Kjell Ahlin, Huan Du Nonlinear system identification with the use of describing functions – a case study. Vibroengineering Procedia, Vol. 8, 2016, p. 33‑38. [Search CrossRef]
8. Gelb A., Vander Velde W. E. Multiple-Input Describing Functions and Nonlinear System Design. McGraw-Hill, 1968. [Search CrossRef]
9. Ahlin K. Toolbox for simulation and parameter identification of nonlinear mechanical systems. Proceedings of IMAC 25, 2007. [Search CrossRef]

#### Cited By

 Advances in Energy and Power Engineering 君 张 2020