*Article* **Approximation of the Statistical Characteristics of Piecewise Linear Systems with Asymmetric Damping and Stiffness under Stationary Random Excitation**

**Tudor Sireteanu 1, Ana-Maria Mitu 1,\*, Ovidiu Solomon 1,2 and Marius Giuclea 1,2**

**Abstract:** In this paper, the dynamic response of piecewise linear systems with asymmetric damping and stiffness for random excitation is studied. In order to approximate the statistical characteristics for each significant output of piecewise linear system, a method based on transmissibility factors is applied. A stochastic linear system with the same transmissibility factor is attached, and the statistical parameters of the studied output corresponding to random excitation having rational spectral densities are determined by solving the associated Lyapunov equation. Using the attached linear systems for root mean square and for standard deviation of displacement, the shift of the sprung mass average position in a dynamic regime, due to damping or stiffness asymmetry, can be predicted with a good accuracy for stationary random input. The obtained results are compared with those determined by the Gaussian equivalent linearization method and by the numerical integration of asymmetric piecewise linear system equations. It is shown that the piecewise linear systems with asymmetrical damping and stiffness characteristics can provide a better vibration isolation (lower force transmissibility) than the linear system.

**Keywords:** asymmetric piecewise linear systems; transmissibility factors; Lyapunov equation

**MSC:** 60H35

#### **1. Introduction**

The limitations of vibration isolation systems with linear passive damping and stiffness characteristics are well known. A high damping ratio is effective in the resonance frequency range but increases the dynamic response of isolation system for higher frequencies. On the other hand, lower damping ratios could be effective above the resonance range with the cost of an unacceptable increase in the dynamic response within the resonance range. Piecewise linear (PWL) systems with asymmetric damping and stiffness characteristics can provide a lower transmissibility factor over the entire frequency range than linear systems.

Many approximate methods have been proposed for studying the vibration of systems with PWL stiffness and damping characteristics [1–6]. The dynamic behavior of PWL systems was studied in [7–9]. A piecewise linear aeroelastic system with and without a tuned vibration absorber was investigated [10]. The experimental results show that the introduction of the piecewise linear stiffness and damping significantly decreases the response amplitude at the primary resonance [11]. The beneficial effect for ride comfort of road vehicles, mainly due to the suspension damping asymmetry, which introduces a downward shift in the mean position of the sprung mass in addition to the vibratory response, has been studied [12–18]. The classical dynamics of the systems with both the statistically uncertain piecewise constant drift and diffusion were extended in [19]. Asymmetric damping forces induce the equilibrium position of the isolated body to shift

**Citation:** Sireteanu, T.; Mitu, A.-M.; Solomon, O.; Giuclea, M. Approximation of the Statistical Characteristics of Piecewise Linear Systems with Asymmetric Damping and Stiffness under Stationary Random Excitation. *Mathematics* **2022**, *10*, 4275. https://doi.org/10.3390/ math10224275

Academic Editor: Davide Valenti

Received: 17 October 2022 Accepted: 11 November 2022 Published: 15 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

downward [20]. A nonlinear interval optimization of asymmetric damper parameters for a racing car is proposed to improve road holding [21].

Various linearization methods have been developed for the analysis of nonlinear systems [22–24]. A very useful property of piecewise linear systems is the independence of their transmissibility factors with respect to the excitation amplitude [25,26]. These factors could be defined as the ratios of root mean square (rms) or standard deviation (std) of output for the same parameters of the harmonic input within the frequency range of practical interest. Therefore, a first order linear differential system can be attached to the considered piecewise linear system so that the first vector component of the attached system has the same transmissibility factor as the chosen output of the nonlinear system. This approach was employed to obtain approximate solutions of PWL systems with piecewiselinear damping with variable friction for application to semi-active control of vibration [23] and for the comparison of the on–off control strategies of vehicle suspensions [24].

In the present work, the Lyapunov equation for attached linear systems is used to approximate the first and second order statistical moments of any significant output of PWL systems with passive asymmetric damping and stiffness. In classical linearization methods, the nonlinear system is replaced by a single equivalent linear system. In the framework of the method used in the present paper, a set of attached linear systems is employed to approximate the statistical characteristics of the PWL system output. Using the attached linear systems for rms and std displacement, the shift of sprung mass average position in dynamic regime, due to damping or stiffness asymmetry, can be predicted with a good accuracy for stationary random input, as confirmed by the numerical results.

In Section 2, the asymmetrical piecewise characteristics are described. In Section 3, the mathematical model of single degree of freedom (SDOF) vibration isolation system with PWL characteristics is presented. In Section 4, the effect of asymmetry of damping and restoring characteristics on the dynamic behavior of piecewise linear systems under stationary random excitation is illustrated. In Section 5, the Gaussian equivalent linearization method for PWL systems is applied. In Section 6, the results obtained by the proposed approach are compared with those given by the Gaussian equivalent linearization method. In order to estimate the statistical characteristics for the output of asymmetric PWL systems, the corresponding attached linear systems are determined in Section 7. In the last section, the statistical characteristics of the simulated output with those calculated by solving Lyapunov equation for corresponding attached linear system are compared. The relative errors show the efficiency and applicability of this method for PWL systems.

#### **2. Modeling the Asymmetrical Piecewise Characteristics**

Figure 1 shows the plots of asymmetrical PWL stiffness in Figure 1a and damping characteristics in Figure 1b, given by

$$F\_{\mathbf{\dot{s}}}(\mathbf{x}) = \begin{cases} \begin{array}{llll} \ k\_1 \mathbf{x}, & \mathbf{x} \le \mathbf{0} \\\ k\_2 \mathbf{x}, & \mathbf{x} > \mathbf{0} \end{array}, \ F\_{\mathbf{d}}(\dot{\mathbf{x}}) = \begin{cases} \ c\_1 \dot{\mathbf{x}}, & \dot{\mathbf{x}} \le \mathbf{0} \\\ c\_2 \dot{\mathbf{x}}, & \dot{\mathbf{x}} > \mathbf{0} \end{cases} \tag{1}$$

where *k*1, *k*2 ≥ 0 are the stiffness coefficients and *c*1, *c*2 ≥ 0 are the damping coefficients and *x* is the travel of vibration isolation system.

Total hysteretic force developed by vibration isolation system for imposed harmonic motion *x*(*t*) = *X* sin *ωt*, where ω = 2π*f* , *f* is the frequency and *X* is the amplitude, is *F*h *x*, . *x* = *F*s(*x*) + *F*d . *x* . The time histories of hysteretic force *F*h *x*, . *x* , stiffness force *F*s(*x*) and damping force *F*d . *x* are illustrated in Figure 2, for parameters values shown in Table 1.

**Figure 1.** Asymmetrical PWL: (**a**) stiffness characteristics; (**b**) damping characteristics.

**Figure 2.** Time histories of forces developed by a vibration isolation system with asymmetric PWL.

**Table 1.** Values of parameters for hysteretic force.

The loops portraying the variation of damping force *F*d and total hysteretic force *F*h versus the imposed displacement *x* are shown in Figure 3.

**Figure 3.** The stiffness characteristic and hysteresis loops portraying the variation of damping force.

The enclosed area by these loops represents *E*d, the energy dissipated per cycle:

$$E\_{\rm d} = \int\_{0}^{\frac{2\pi}{w}} F\_{\rm d}(\dot{x}) \dot{x} dt = 0.5 \pi w X^{2} (c\_{1} + c\_{2}) \tag{2}$$

Figure 4 depicts the schematic model of a device with asymmetrical damping and stiffness characteristics.

**Figure 4.** Design principle of a device with PWL asymmetric damping and stiffness characteristics.

The metallic bellows, filled with hydraulic fluid, are welded at both ends, and, therefore, the fluid damper is leak proof. The asymmetry of damping force is controlled by the openings of extension and compression valves. The dimensions of valve openings and fluid viscosity must be assessed such that to have laminar flow within the range of damper operating conditions. Since the bellows geometry is identical, there is no need for any volume compensation system. The suspension springs with different stiffness are in the unloaded condition (free length) when the isolation system is in the equilibrium position. Each of them has only one end fixed on the device structure. Therefore, they work only as compression springs for both extension and compression strokes. The bellows longitudinal stiffness, being much smaller than the stiffness of springs, is neglected.

#### **3. Mathematical Model of SDOF Vibration Isolation System with PWL Characteristics**

Vibration isolation systems are widely used to reduce the dynamic forces transmitted from the base input to sprung mass (Figure 5) or from the sprung mass to the system base (Figure 6).

**Figure 5.** PWL system for mitigation the dynamic forces transmitted from the base input to the sprung mass.

**Figure 6.** PWL system for mitigation the dynamic forces transmitted from the sprung mass to the system base.

The equation of motion for both vibration isolation systems, shown in Figures 5 and 6, can be written as:

$$m\ddot{\mathbf{x}} + F\_{\mathbf{d}}(\dot{\mathbf{x}}) + F\_{\mathbf{s}}(\mathbf{x}) = P\_0(t) \tag{3}$$

where *x* = *x*1 − *x*0 is the relative displacement of sprung mass, *P*0(*t*) = −*m*.. *x*0 is the input of system shown in Figure 5, *x*1 is the absolute displacement and *x*0 is the base displacement. For the system depicted in Figure 6, *x* is the absolute displacement of sprung mass, relative to its static equilibrium position, and *P*0(*t*) = *F*0(*t*) is the force applied to the sprung mass. In both cases, *x* is the stroke (travel) of sprung mass suspension and will be called displacement (disp). The main output of interest for vibration isolation systems are the absolute accelerations of sprung mass .. *x*1, for system shown in Figure 5, and .. *x*, for system shown in Figure 6. The absolute acceleration is a measure for mitigation of dynamic forces transmitted through the sprung mass suspension. In the rest of the paper it will be called acceleration and abbreviated as acc. The analytic expressions of asymmetric damping and elastic characteristics *F*d . *x* and *F*s(*x*) can be written as

$$\begin{array}{l} F\_{\mathbf{d}}(\dot{\mathbf{x}}) = 0.5 \left[ c\_1 (1 - \text{sgn}\dot{\mathbf{x}}) + c\_2 (1 + \text{sgn}\dot{\mathbf{x}}) \right] \dot{\mathbf{x}}, \\ F\_{\mathbf{s}}(\mathbf{x}) = 0.5 \left[ k\_1 (1 - \text{sgn}\mathbf{x}) + k\_2 (1 + \text{sgn}\mathbf{x}) \right] \mathbf{x}. \end{array} \tag{4}$$

Introducing the notations

$$\begin{array}{l} \mathfrak{w}\_{1} = 2\pi f\_{1} = \sqrt{\frac{k\_{1}}{m}}, \; \mathfrak{w}\_{2} = 2\pi f\_{2} = \sqrt{\frac{k\_{2}}{m}}, \; \mathfrak{f}\_{1} = \frac{c\_{1}}{2\mathfrak{w}\_{1}m'}, \; \mathfrak{f}\_{2} = \frac{c\_{2}}{2\mathfrak{w}\_{2}m} \\ \mathfrak{f} = \frac{c\_{2}}{c\_{1}} = \frac{\mathfrak{f}\_{2}f\_{2}}{\mathfrak{f}\_{1}f\_{1}}, \; \gamma = \sqrt{\frac{k\_{2}}{k\_{1}}} = \frac{f\_{2}}{f\_{1}} \\ \mathfrak{f}\_{\mathsf{d}}(\dot{\mathfrak{x}}) = \frac{\mathsf{F}\_{\mathsf{d}}(\dot{\mathfrak{x}})}{m}, \; f\_{\mathsf{s}}(\mathsf{x}) = \frac{\mathsf{F}\_{\mathsf{d}}(\mathbf{x})}{m} \text{ and } p\_{0}(t) = \frac{p\_{0}(t)}{m}, \end{array} \tag{5}$$

the equation of motion (3) becomes

$$
\ddot{\mathbf{x}} + f\_{\mathbf{d}}(\dot{\mathbf{x}}) + f\_{\mathbf{s}}(\mathbf{x}) = p\_0(t) \tag{6}
$$

where

$$\begin{array}{l} f\_{\mathbf{d}}(\dot{\mathbf{x}}) = \zeta\_1 \omega\_1 \left[ (\beta + 1)\dot{\mathbf{x}} + (\beta - 1)|\dot{\mathbf{x}}| \right]\_{\mathbf{y}} \\ f\_{\mathbf{s}}(\mathbf{x}) = 0.5 \omega\_1^2 \left[ (\gamma^2 + 1)\mathbf{x} + (\gamma^2 - 1)|\mathbf{x}| \right]. \end{array} \tag{7}$$

From (2) and (5), one can see that asymmetry parameter β is the ratio of dissipated energy per rebound *E*d2 = 0.5πω*X*2*c*2 and bound *E*d1 = 0.5πω*X*2*c*1 strokes for an imposed harmonic motion.

#### **4. The Effect of Asymmetry of Damping and Restoring Characteristics on the Dynamic Behavior of Piecewise Linear Systems under Stationary Random Excitation**

In general, the asymmetry of damping or stiffness characteristics leads to a drift of sprung mass average position in dynamic regime, different from its static equilibrium position. Nevertheless, by a suitable combination of the asymmetry parameters β and *γ*, one can obtain outputs of PWL systems with almost no drift.

Suppose that *p*0(*t*) is a stationary Gaussian random process with zero mean and standard deviation σ0. If *x*(*t*) is the steady state stationary solution of Equation (6), with constant mean value *m*x= E[*x*] then E - . *x* . = E -.. *x* . = 0. Therefore, by applying the average operator corresponding to joint distribution of the output of Equation (6), *m*x is obtained as follows:

$$m\_{\mathbf{x}} = -\frac{m\_{|\mathbf{x}|}}{\gamma^2 + 1} \left[ \left( \gamma^2 - 1 \right) + 2(\beta - 1)\zeta\_1 \frac{m\_{|\mathbf{x}|}}{\omega\_1 m\_{|\mathbf{x}|}} \right] \tag{8}$$

where

$$m\_{|\mathbf{x}|} = \operatorname{E}[|\mathbf{x}|] \text{ and } m\_{|\dot{\mathbf{x}}|} = \operatorname{E}\left[|\dot{\mathbf{x}}|\right] \tag{9}$$

The relation (8) shows that *m*x = 0 if *γ* = 1, β = 1; *m*x < 0 if *γ* > 1, β > 1 and *m*x > 0 if *γ* < 1, β < 1. It is worth noting that by assuming *m*| . x| /ω1*m*|x| ≈ 1, for all case studies considered in this work (including *γ* > 1, β < 1 or *γ* < 1, β > 1), the sign of *m*x could be predicted by determining the sign of expression *S*(*γ*, β, *ζ*1) = −-*γ*2 − 1 + 2(β − 1)*ζ*1 . , without being necessary the numerical simulation values from (9).

#### **5. Gaussian Equivalent Linearization Method of PWL Systems**

The Gaussian equivalent linear system (LinEq) of system (6), where *p*0(*t*) is a stationary Gaussian process, E[*p*0(*t*)] = 0 and E*p*2 0(*t*) . = σ2 0 is written as

$$
\ddot{\mathbf{x}} + 2\mathbb{Z}\_{\mathbf{e}}\omega\_{\mathbf{e}}\dot{\mathbf{x}} + \omega\_{\mathbf{e}}^2 \mathbf{x} = p\_0(t) \tag{10}
$$

The joint probability density function of the Gaussian stationary solution of equivalent linear system is

$$\begin{array}{c} \text{g}(\boldsymbol{x}, \dot{\boldsymbol{x}}) = \operatorname{g}\_1(\boldsymbol{x}) \operatorname{g}\_2(\dot{\boldsymbol{x}})\\ \text{g}\_1(\boldsymbol{x}) = \frac{1}{\sqrt{2\pi}\sigma\_\mathbf{x}} \exp\left[-\frac{\boldsymbol{x}^2}{2\sigma\_\mathbf{x}^2}\right], \operatorname{g}\_2(\dot{\boldsymbol{x}}) = \frac{1}{\sqrt{2\pi}\sigma\_\mathbf{x}} \exp\left[-\frac{\dot{\boldsymbol{x}}^2}{2\sigma\_\mathbf{x}^2}\right] \end{array} \tag{11}$$

where σx = σx(ωe, *ζ*e) and σ . x = σ . x(ωe, *ζ*e) are the standard deviations for the solution of Equation (10).

The variance of the acceleration .. *x*1 = 2*ζ*eωe . *x* + ω2 e*x* of the equivalent linear system is

$$
\sigma\_{\ddot{\mathbf{x}}\_1}^2 = E\left[\ddot{\mathbf{x}}\_1^2\right] = 4\zeta\_\mathbf{e}^2\omega\_\mathbf{e}^2 E\left[\dot{\mathbf{x}}^2\right] + \omega\_\mathbf{e}^4 E\left[\mathbf{x}^2\right] = \omega\_\mathbf{e}^2 \left[4\zeta\_\mathbf{e}^2\sigma\_\mathbf{x}^2 + \omega\_\mathbf{e}^2\sigma\_\mathbf{x}^2\right] \tag{12}
$$

Applying the linearization criteria,

$$\begin{aligned} \varepsilon\_{\mathsf{e}}(\omega\_{\mathsf{e}}) &= E\left[\left(f\_{\mathsf{s}}(\mathsf{x}) - \omega\_{\mathsf{e}}^{2}\mathsf{x}\right)^{2}\right] = \min\_{\mathsf{e}} \qquad \frac{\partial \varepsilon\_{\mathsf{e}}(\omega\_{\mathsf{e}})}{\partial \omega\_{\mathsf{e}}} = 0\\ \varepsilon\_{\mathsf{d}}(\zeta\_{\mathsf{e}}, \omega\_{\mathsf{e}}) &= E\left[\left(f\_{\mathsf{d}}\left(\dot{\mathsf{x}}\right) - 2\zeta\_{\mathsf{e}}\omega\_{\mathsf{e}}\dot{\mathsf{x}}\right)^{2}\right] = \min\_{\mathsf{d}} \ \frac{\partial \varepsilon\_{\mathsf{d}}(\zeta\_{\mathsf{e}}, \omega\_{\mathsf{e}})}{\partial \zeta\_{\mathsf{e}}} = 0 \end{aligned} \tag{13}$$

the linear equivalent stiffness and damping coefficients are obtained using (6), (11) and (13):

$$\begin{array}{l} \omega\_{\mathbf{e}}^{2} = \frac{\mathbb{E}\left[\boldsymbol{x}f\_{\mathbf{s}}(\mathbf{x})\right]}{\mathbb{E}\left[\boldsymbol{x}^{2}\right]} = \frac{\omega\_{1}^{2} + \omega\_{2}^{2}}{2} = \frac{\omega\_{1}^{2}\left(1 + \gamma^{2}\right)}{2},\\ \zeta\_{\mathbf{e}} = \frac{\mathbb{E}\left[\boldsymbol{x}f\_{\mathbf{d}}\left(\mathbf{x}\right)\right]}{2\omega\_{\mathbf{e}}\mathbb{E}\left[\boldsymbol{x}^{2}\right]} = \frac{\zeta\_{1}\omega\_{1} + \zeta\_{2}\omega\_{2}}{2\omega\_{\mathbf{e}}} = \frac{\zeta\_{1}\omega\_{1}(1 + \beta)}{2\omega\_{\mathbf{e}}}. \end{array} \tag{14}$$

Using (14) one can write

$$f\_{\mathbf{e}} = f\_1 \sqrt{\frac{1+\gamma^2}{2}}, \ \zeta\_{\mathbf{e}} = \frac{\zeta\_1(1+\beta)}{\sqrt{2(1+\gamma^2)}}, \ f\_{\mathbf{e}} = \omega\_{\mathbf{e}}/2\pi, \ f\_1 = \omega\_1/2\pi \tag{15}$$

In order to highlight the advantage of using the vibration isolation systems with asymmetric PWL characteristics, the obtained results are compared with those of optimal linear equivalent system.

For given values of linear equivalent system *ζ*e, *f*e and chosen values of asymmetry parameters *γ*, β, relations (15) yield:

$$f\_1 = f\_{\text{e}} \sqrt{\frac{2}{1 + \gamma^2}}, \ f\_2 = \gamma f\_{1\prime} \ \zeta\_1 = \frac{\zeta\_{\text{e}} \sqrt{2(1 + \gamma^2)}}{1 + \beta} \text{ and } \zeta\_2 = \frac{\beta}{\gamma} \zeta\_1 \tag{16}$$

From (16) one can obtain the balance equation between the energy dissipated by PWL asymmetric system and its linear equivalent system per cycle for same imposed harmonic motion:

$$
\mathbb{Z}\_1 f\_1 + \mathbb{Z}\_2 f\_2 = 2f\_\text{e} \mathbb{Z}\_\text{e} \tag{17}
$$

As one can see from previous relations, there are an infinite number of PWL asymmetric systems having same linear equivalent system.

Following [27], the standard deviation of the stationary steady state acceleration of sprung mass for stochastic linear system (10) with Gaussian white noise excitation *p*0(*t*) and constant spectral density *S*0 is

$$
\sigma\_{\bar{\mathbf{x}}\_{\text{l}}\mathbf{e}} = \sqrt{2S\_0 \int\_0^\infty A\_{\bar{\mathbf{x}}\_{\text{l}}\mathbf{e}}^2(\omega) \, \mathrm{d}\omega} = \sqrt{\frac{\pi \omega \mathbf{e}\_\mathbf{e} S\_0 (1 + 4\zeta\_\mathbf{e}^2)}{2\zeta\_\mathbf{e}}},\tag{18}
$$

where *A* .. x1e(ω) is the acceleration transmissibility factor of linear equivalent system:

$$A\_{\bar{\mathbf{x}}\_{\text{l}^{\text{e}}}}(\omega) = \sqrt{\frac{4\zeta\_{\text{e}}^{2}\omega\_{\text{e}}^{2}\omega^{2} + \omega\_{\text{e}}^{4}}{\omega^{4} + 2(2\zeta\_{\text{e}}^{2} - 1)\omega\_{\text{e}}^{2}\omega^{2} + \omega\_{\text{e}}^{4}}} \tag{19}$$

The optimum value of damping ratio *ζ*e, which minimizes the std value of sprung mass acceleration is *ζ*e = 0.5, and its minimum value is σ .. x1min = √2πωe*S*0. Taking *S*0 = 1 m2s−3 and ωe = 2πrad/s, the optimum std value of acceleration is σ.. x1min ∼= 6.28 ms−2. For numerical integration, the input is a limited bandwidth white noise, and std value of acceleration is calculated as

$$
\sigma\_{\ddot{\mathbf{x}}\_1 \mathbf{e}} \cong \sqrt{2 \int\_{\omega\_{\min}}^{\omega\_{\max}} A\_{\ddot{\mathbf{x}}\_1 \mathbf{e}}^2(\omega; \zeta\_{\mathbf{e}}, \omega\_{\mathbf{e}}) \, \mathrm{d}\omega} \cong 6.21 \,\mathrm{ms}^{-2}},\tag{20}
$$

where 0.2 ≤ ω ≤ 128 and ω is measured in rad/s, which is a good approximation of optimum value σ .. x1min, calculated over the whole range of angular frequency [0, ∞). The results obtained by the proposed approach will be compared with those obtained by the Gaussian equivalent linearization method.

#### **6. The Response of PWL Systems to Stationary Gaussian Random Input with Rational Spectral Density (Shape Filtered White Noise)**

According to [28], the covariance function and spectral density of stationary Gaussian random input *p*0(*t*) are

$$\begin{array}{l} \mathbb{C}\_{0}(\tau) = \sigma\_{0}^{2} e^{-a|\tau|} \cos b\tau, \ a > 0, \ b \ge 0, \\\ S\_{0}(\omega) = \frac{\sigma\_{0}^{2}}{\pi} \frac{a\left(\omega^{2} + a^{2} + b^{2}\right)}{\omega^{4} + 2\left(a^{2} - b^{2}\right)\omega^{2} + \left(a^{2} + b^{2}\right)^{2}} \end{array} \tag{21}$$

where σ2 0 =  ∞ −∞ *S*0(ω)dω = 2 ∞ 0 *S*0(ω)dω . The above expression of *S*0(ω) can be viewed as the spectral density of the Gaussian stationary random process *p*0(*t*), obtained as the output of a second order shape filter to a stationary Gaussian white noise process *z*(*t*) with E[*z*(*t*)]= 0, E[*z*(*t*)*z*(*t* + τ)] = 2*πS*0*δ*(τ), where *δ*(τ) is the Dirac delta function. In order to determine the equations of the second order shape filter, the spectral density (21) is written under the form

$$S\_0(\omega) = \frac{\left|P(\text{i}\omega)\right|^2}{\left|Q(\text{i}\omega)\right|^2} = \frac{\left|b\_0(\text{i}\omega) + b\_1\right|^2}{\left|a\_0(\text{i}\omega)^2 + a\_1(\text{i}\omega) + a\_2\right|^2},\tag{22}$$

where *b*0 = σ0 *a* π , *b*1 = σ0 *a*(*a*2+*b*2) π , *a*0 = 1, *a*1 = 2*a*, *a*2 = *a*2 + *b*2.

The output *u*1(*t*) of the following first order differential system with the white noise excitation *z*(*t*) is a Gaussian stationary random process with spectral density *S*0(ω):

$$
\dot{\mathbf{u}} = \mathbf{A}\mathbf{u} + \mathbf{g}z,\tag{23}
$$

where

$$\begin{aligned} \mathbf{A} &= \begin{bmatrix} 0 & 1 \\ -a\_2 & -a\_1 \end{bmatrix}, \mathbf{u} = \begin{bmatrix} u\_1 \\ u\_2 \end{bmatrix}, \mathbf{g} = \begin{bmatrix} \mathcal{G}1 \\ \mathcal{G}2 \end{bmatrix} \\\ g\_1 &= \sigma\_0 \sqrt{\frac{\pi}{\pi}}, \ g\_2 = \sigma\_0 \sqrt{\frac{\pi}{\pi}} \left( \sqrt{a^2 + b^2} - 2a \right) \text{ and } p\_0(t) = u\_1(t). \end{aligned} \tag{24}$$

In order to study the behavior of asymmetric PWL systems excited by stationary random input with rational spectral density, a linear system of first order stochastic equations is assessed such as the first component of its solution vector has the same transmissibility factor as the chosen output of the considered piecewise linear system [23]. The statistical parameters of obtained stochastic differential equations are determined by solving the associated Lyapunov equation.

Since the mean value of PWL acceleration response system has zero mean, the transmissibility factors corresponding to standard deviation and root mean square values are identical.

The discrete values of transmissibility factor corresponding to standard deviation of acceleration .. *x*1(*t*) is defined as:

$$\tilde{A}\_{\bar{\mathbf{x}}\_{\bar{1}}}(\omega\_{i}) = \frac{\sigma\_{\bar{x}\_{1}}(\omega\_{i})}{\sigma\_{p\_{0i}}} = \frac{\sqrt{2}\,\sigma\_{\bar{x}\_{1}}(\omega\_{i})}{P\_{0}}, \; i = \overline{1,N} \tag{25}$$

These values are obtained by numerical integration of Equation (6), using Matlab Simulink, for harmonic inputs with constant amplitude and different frequencies in the twelfth octave band:

$$p\_{0i}(t) = P\_0 \sin \omega\_i t, \ \omega\_i = 2\pi f\_{i\prime} \ f\_1 = 0.3 \text{ Hz}, f\_N = 20.05 \text{ Hz}, \ f\_i = 2^{\frac{i-1}{12}} f\_1, \ i = \overline{1, N} \quad \text{(26)}$$

where *P*0 = 1 m/s2, *N* = log(*fN*/ *f*1) log(*f*2/ *f*1) + 2 = 114.

It should be mentioned that the transmissibility factors of PWL systems, with asymmetry type (affine) [29], considered in this paper, do not depend on the amplitude *P*0 of the applied harmonic input with variable frequency, as long as they are computed for the stationary regime. The numerical values *A*\* .. x1 (ω*i*), *i* = 1, *N*, can be fitted using the Least Squares Method, by analytical expressions having the form

$$A\_{\ddot{\chi}\_1}(\omega) \;= \sqrt{\frac{P\_1 \omega^2 + P\_2}{\omega^4 + Q\_1 \omega^2 + Q\_2}}. \tag{27}$$

The transmissibility factor *A* .. x1 (ω) is written as:

$$A\_{\ddot{\mathbf{x}}\_{1}}(\omega) := \left| \frac{b\_{0}(\mathbf{i}\omega) + b\_{1}}{\left(\mathbf{i}\omega\right)^{2} + a\_{1}(\mathbf{i}\omega) + a\_{2}} \right| = \sqrt{\frac{b\_{0}^{2}\omega^{2} + b\_{1}^{2}}{\omega^{4} + \left(a\_{1}^{2} - 2a\_{2}\right)\omega^{2} + a\_{2}^{2}}} \tag{28}$$

From relations (27) and (28), the following nonlinear algebraic systems of equations for unknown coefficients *b*0, *b*1, *a*1, *a*2, are obtained

$$\begin{cases} \begin{array}{c} b\_0^2 = P\_1 \\ b\_1^2 = P\_2 \end{array} \end{cases} \begin{cases} a\_1^2 - 2a\_2 = Q\_1 \\ a\_2^2 = Q\_2 \end{cases} \tag{29}$$

The attached linear system corresponding to transmissibility factor (28) can be written as

$$
\dot{\mathbf{u}} = \mathbf{A}\mathbf{u} + \mathbf{c}p\_{0\prime} \tag{30}
$$

where

$$\mathbf{A} = \begin{bmatrix} 0 & 1 \\ -a\_2 & -a\_1 \end{bmatrix}, \ \mathbf{u} = \begin{bmatrix} u\_1 \\ u\_2 \end{bmatrix}, \ \mathbf{c} = \begin{bmatrix} c\_1 \\ c\_2 \end{bmatrix}, \ \begin{cases} & c\_1 = b\_0 \\ c\_2 = b\_1 - a\_1 c\_1 \end{cases} \tag{31}$$

The transmissibility factor *A* .. x1 (ω) ∼= *A*σ*u*1 (ω) = √2σ*u*1/*P*0, where *u*1 is the first component of the solution vector **u**. The system (30) is asymptotically stable if *a*1, *a*2 > 0. Therefore, from the sets of real solutions of (29), one must select only the solutions that fulfill these conditions.

In what follows, the study is carried out for several asymmetric PWL systems for which the stochastic equivalent linear system is the optimal one. The parameters of PWL systems, given in Table 2, were obtained by using relations (16).

**Table 2.** The parameters of PWL systems.

In Figure 7, the transmissibility factor *A*\* .. x1 (ω*i*) obtained by numerical integration for the asymmetric PWL systems from Table 2 is compared with transmissibility factor of their stochastic equivalent linear system (*ζ* 1 = *ζ*2 = 0.5, *γ* = β = 1, *f*1 = *f*2 = 1 Hz).

**Figure 7.** The transmissibility factors for linear and asymmetric PWL systems.

Table 3 presents the standard deviation values of the acceleration obtained by using a similar relation to (20), for the transmissibility factors of asymmetric PWL systems shown in Figure 7.

**Table 3.** The transmissibility factors of asymmetric PWL systems.

The above results show that the piecewise linear systems with asymmetrical damping and stiffness characteristics can provide a better vibration isolation (lower force transmissibility) than the optimum equivalent linear system (σ .. x1 = 6.21 [m/s2]).

#### **7. Attached Linear System for Different Outputs of PWL Systems Excited by a Second Order Shape Filtered White Noise**

In order to estimate the statistical characteristics for the output of asymmetric PWL systems, the corresponding attached linear systems will be determined in the next sections. The stochastic equations of attached linear system for the piecewise linear system (6), with shape filtered white noise excitation, is obtained by combining Equations (23) and (30):

**.**

$$
\dot{\mathbf{u}} = \mathbf{A}\mathbf{u} + \mathbf{g}z,\tag{32}
$$

where

$$\mathbf{A} = \begin{bmatrix} 0 & 1 & c\_1 & 0 \\ -a\_2 & -a\_1 & c\_2 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & -a\_4 & -a\_3 \end{bmatrix}, \ \mathbf{u} = \begin{bmatrix} u\_1 \\ u\_2 \\ u\_3 \\ u\_4 \end{bmatrix}, \ \mathbf{g} = \begin{bmatrix} 0 \\ 0 \\ \mathcal{S}\_1 \\ \mathcal{S}\_2 \end{bmatrix} \tag{33}$$

$$a\_3 = 2a, \ a\_4 = a^2 + b^2, \mathcal{S}\_1 = \sigma\_0 \sqrt{\frac{\pi}{\pi'}}, \mathcal{S}\_2 = \sigma\_0 \sqrt{\frac{\pi}{\pi}} (\sqrt{a^2 + b^2} - 2a)$$

The covariance matrix **C** = *cij* , *cij* = *cji* = lim*t*→∞ E *ui*(*t*)*uj*(*t*) . , *i*, *j* = 1, 4 of the steady state stationary solution of stochastic linear system (32) satisfies [30] the Lyapunov Equation:

$$\mathbf{AC} + \mathbf{CA}^{\mathrm{T}} + 2\pi \mathbf{S}\_{0} \mathbf{gg}^{\mathrm{T}} = 0 \tag{34}$$

The standard deviation of the acceleration is estimated by σ .. x1 ∼= σu 1 where σu 1 = √*c*11. The values of σ .. x1 obtained by using Lyapunov equation will be compared with those determined for linear equivalent system (10) where *ζ*e = 0.5 and ωe = 2π rad/s.

The values of transmissibility factors *A*\*xrms(ω*i*) = Ψx(ω*i*) σp0*i* = √2 Ψx(ω*i*) P0 , *i* = 1, *N* , and *A*\*xstd(ω*i*) = σx(ω*i*) σp0*i* = √2σx(ω*i*) P0 , *i* = 1, *N* corresponding to rms Ψx and std σx of relative displacement *x*(*t*) are obtained by numerical integrations. These values can be approximated by rational expressions having the form

$$A\_{\rm X}(\omega) \;= \sqrt{\frac{P\_1 \omega^4 + P\_2 \omega^2 + P\_3}{\omega^6 + Q\_1 \omega^4 + Q\_2 \omega^2 + Q\_3}} \tag{35}$$

The transmissibility factor is written as

$$A\_{\mathbf{x}}(\omega) = \left| \frac{b\_0(\mathbf{i}\omega)^2 + b\_1(\mathbf{i}\omega) + b\_2}{(\mathbf{i}\omega)^3 + a\_1(\mathbf{i}\omega)^2 + a\_2(\mathbf{i}\omega) + a\_3} \right| = \sqrt{\frac{b\_0^2 \omega^4 + \left(b\_1^2 - 2b\_0 b\_2\right) \omega^2 + b\_2^2}{\omega^6 + \left(a\_1^2 - 2a\_2\right) \omega^4 + \left(a\_2 - 2a\_1 a\_3\right) \omega^2 + a\_3^2}}\tag{36}$$

From relations (35) and (36) one can obtain the following algebraic systems of equations for unknown coefficients *b*0, *b*1, *b*2, *a*1, *a*2, *a*3:

$$\begin{cases} b\_0^2 = P\_1\\ b\_1^2 - 2b\_0b\_2 = P\_2\\ b\_2^2 = P\_3 \end{cases}, \begin{cases} a\_1^2 - 2a\_2 = Q\_1\\ a\_2^2 - 2a\_1a\_3 = Q\_2\\ a\_3^2 = Q\_3 \end{cases} \tag{37}$$

The equations of the attached linear system having the same transmissibility factor (35) can be written as **.**

$$
\dot{\mathbf{u}} = \mathbf{A}\mathbf{u} + \mathbf{c}p\_0 \tag{38}
$$

where

$$\mathbf{A} = \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ -a\_3 & -a\_2 & -a\_1 \end{bmatrix}, \mathbf{u} = \begin{bmatrix} u\_1 \\ u\_2 \\ u\_3 \end{bmatrix}, \mathbf{c} = \begin{bmatrix} c\_1 \\ c\_2 \\ c\_3 \end{bmatrix}, \begin{cases} c\_1 = b\_0 \\ c\_2 = b\_1 - a\_1 c\_1 \\ c\_3 = b\_2 - a\_2 c\_1 - a\_1 c\_2 \end{cases} \tag{39}$$

The system (38) is asymptotically stable if *ai* > 0, for *i* = 1, 2, 3. The covariance function and spectral density of system input *p*0(*t*) are given by (21). The attached system of stochastic differential equations with white noise excitation is given by

$$
\dot{\mathbf{u}} = \mathbf{A}\mathbf{u} + \mathbf{g}z \tag{40}
$$

where

$$u\_1(t) = \mathbf{x}, \ u\_2(t) = \dot{\mathbf{x}}(t), \ u\_4 = p\_0(t), \ \mathbf{A} = \begin{bmatrix} 0 & 1 & 0 & c\_1 & 0 \\ 0 & 0 & 1 & c\_2 & 0 \\ -a\_3 & -a\_2 & -a\_1 & c\_3 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & -a\_5 & -a\_4 \end{bmatrix}, \ \mathbf{u} = \begin{bmatrix} u\_1 \\ u\_2 \\ u\_3 \\ u\_4 \\ u\_5 \end{bmatrix}, \ \mathbf{g} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ \mathcal{G}\_1 \\ \mathcal{G}\_2 \end{bmatrix} \tag{41}$$

 $a\_4 = 2a$ ,  $a\_5 = a^2 + b^2$ ,  $g\_1 = \sigma\_0 \sqrt{\frac{\pi}{\pi}}$ ,  $g\_2 = \sigma\_0 \sqrt{\frac{\pi}{\pi}} \left(\sqrt{a^2 + b^2} - 2a\right)$

The rms and std values of relative displacements of PWL system, Ψx and σx, can be approximated as Ψx ∼= √*c*11 rms and σx ∼= √*c*11 std. The values of *c*11 rms and *c*11 std are the first elements of covariance matrices **C**rms *cij*rms , **C**std *cij*std , *i*, *j* = 1, 4, obtained by solving the Lyapunov Equation (34), corresponding to attached linear systems for the transmissibility factors *A*xrms(ω) and *A*xstd(ω), respectively. The mean displacement of asymmetric PWL system is approximated by

$$m\_{\rm xLyap} \cong \sqrt{c\_{11\,\rm rms} - c\_{11\,\rm std}} \text{sgn}(\mathcal{S}(\gamma, \beta, \mathcal{J}\_1))\tag{42}$$

#### **8. Numerical Results**

In this section, the statistical characteristics of simulated output are compared with those calculated by solving the Lyapunov Equation (34) for corresponding attached linear systems. The length and sampling interval of simulated filtered white noise input *p*0(*t*) were *T* = 100 s and Δ*t* = 0.001 s. The results obtained for the study cases, given in Table 2 for PWL asymmetric systems and for their linear equivalent system (*ζ*e = 0.5, ωe = 2π rad/s), are presented in Tables 4–6.

**Table 4.** The std values and relative errors for absolute accelerations.

**Table 5.** The rms values and relative errors for relative displacements.

**Table 6.** The mean values and relative errors for relative displacements.

The last column of Table 6 shows the mean values of displacement, evaluated by using in (8) the values m|x|, *m*| . x| obtained by numerical integration of PWL equation of motion (6). It worth noting that the optimum value of damping ratio for a linear system with undamped eigenfrequency ω1 = 2π rad/s and considered random input *p*0(*t*), is *ζ*opt = 0.55. The value of standard deviation of simulated acceleration output obtained in this case is σ .. x1 opt = 0.741 m/s2 .

Table 4 shows that the simulated values σ .. x1sim are better approximated by using the proposed method than the Gaussian equivalent linearization method. Therefore, in all case studies the asymmetric PWL systems provide better vibration isolation than the optimum linear system, for both considered random inputs (band limited and shape filtered white noise).

The results presented in Tables 4 and 6 show that the relative errors of approximation between the results obtained by numerical integration of asymmetric PWL systems and those calculated by using the Lyapunov equation for linear attached systems are less than 7.5% for standard deviation of acceleration and less than 13% for mean value of displacement. As one can see, from Tables 2 and 6, as nonlinearity increases, the mean value displacement is better approximated. It should be mentioned that the Gaussian equivalent linearization method cannot provide any information about the drift of sprung mass average position in dynamic regime, as it is shown in Table 5.

In order to illustrate the application of presented method, the case 1 from Table 2, which display the strongest nonlinearity, has been chosen. In Figures 8 and 9 are plotted

the transmissibility factors, simulated and fitted for this case, as well as the values of parameters from fitting the curves given by the expressions (27) and (35).

**Figure 8.** Simulated and fitted transmissibility factors for acceleration.

**Figure 9.** Simulated and fitted transmissibility factors for std and rms displacement.

In Tables 7 and 8 are given the coefficients of attached linear systems corresponding to acceleration and displacement, obtained by solving the algebraic Equations (29), (31), (37) and (39) for parameters shown in Figures 8 and 9.

**Table 8.** Coefficients of attached linear system for std and rms displacement.

In the last column of these tables are given the values of elements *c*11 from covariance matrices obtained by solving the Lyapunov Equation (34), for the corresponding attached linear systems. Using these coefficients, are obtained the values of std acceleration σ .. x1Lyap = √*c*11 = 0.515 ms−2 and the mean value of displacement *m*xLyap ∼= −√*c*11 rms − *c*11 std= −0.18 m, according to (8). Figure 10 shows the first 30 s from the simulated time histories of input, acceleration and displacements outputs for PWL, attached linear (rms for displacement) and linear equivalent systems, obtained for case study 1.

**Figure 10.** Acceleration output of PWL, equivalent linear and attached linear systems for case 1.

In Figure 11 are plotted the spectral densities of acceleration output, determined by 1/3 octave band-pass filtering for PWL, linear equivalent and the attached linear systems.

**Figure 11.** Spectral densities of acceleration for PWL, linear equivalent and attached linear systems.

The relative errors between the areas under spectral densities that represent the variances of acceleration, given in Table 9, advocate the efficiency of the proposed method.

**Table 9.** Variances of acceleration for PWL, LinEq and LinAtt systems.

#### **9. Conclusions**

The dynamic response of piecewise linear systems with asymmetric damping and stiffness for random inputs is approximated by a method based on transmissibility factors. The application of this method does not require the numerical simulation of input and output time histories, except for obtaining the transmissibility factors by using harmonic inputs. Using these frequency characteristics, a stochastic linear system is attached for each variable of interest. The statistical parameters of the studied output corresponding to random excitations having rational spectral densities are determined by solving the associated Lyapunov equation.

The obtained results are compared with those determined by the numerical integration of asymmetric PWL response. The relative errors show the efficiency and applicability of this method for PWL systems. In addition, this approach allows the realization of vibration isolation systems with better performance than those with linear characteristics. Using the attached linear systems for rms and std displacement, the shift of sprung mass average position in dynamic regime, due to damping or stiffness asymmetry, can be predicted with a good accuracy for stationary random input.

**Author Contributions:** Conceptualization, T.S., A.-M.M., O.S. and M.G.; methodology, T.S., A.-M.M., O.S. and M.G.; writing—original draft preparation, T.S., A.-M.M., O.S. and M.G.; writing—review and editing, T.S., A.-M.M., O.S. and M.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

#### **References**

**Gurami Tsitsiashvili 1,\*,†, Marina Osipova 1,2 and Yury Kharchenko 1**

**Abstract:** In this paper, we solve the problem of estimating the parameters of a system of ordinary differential equations from observations on a short interval of argument values. By analogy with linear regression analysis, a sufficiently large number of observations are selected on this segment and the values of the functions on the right side of the system and the values of the derivatives are estimated. According to the obtained estimates, unknown parameters are determined, using the differential equations system. The consistency of the estimates obtained in this way is proved with an increase in the number of observations over a short period of argument values. Here, an algorithm for estimating parameters acts as a system. The error of the obtained estimate is an indicator of its quality. A sequence of inaccurate measurements is a random process. The method of linear regression analysis applied to an almost linear regression function is used as an optimization procedure.

**Keywords:** system of ordinary differential equations; linear regression analysis; theorem of existence and uniqueness; implicit function theorem; method of moments

**MSC:** 60J28

#### **1. Introduction**

The problem of estimating the parameters of a system of nonlinear ordinary differential equations, based on inaccurate deterministic observations, using known optimization algorithms, is solved in the papers [1–3]. An alternative approach for estimating the parameters of a deterministic recurrent sequence, observed with random additive and multiplicative errors, based on the relationships between the trajectory averages and their approximation from inaccurate observations, is proposed in [4,5].

The advantage of the first approach is the possibility of using known optimization algorithms, and the disadvantage of it is the lack of analytical estimates of the convergence rate to the estimated parameters. The advantage of the second approach is the availability of theoretical estimates of the convergence rate to the estimated parameters, and the disadvantage of it is the need to establish limit cycles or limit distributions for recurrent sequences.

Despite all the differences in these approaches, the common fact is that by increasing in the length of the observation segment, the accuracy of estimates increases and, under certain conditions, may tend to zero. At the same time, the problem of estimating parameters over a small observation interval is interesting, which is closely related to discrete optimization methods of experiment planning (see, for example, [6,7]).

In this paper, this problem is solved for a system of non-linear ordinary differential equations. At the same time, the estimation of the parameters of this system, based on inaccurate observations, is solved under the assumption that a large number of observations may be carried out over a relatively short segment. To estimate the parameters, the method

**Citation:** Tsitsiashvili, G.; Osipova, M.; Kharchenko, Y. Estimating the Coefficients of a System of Ordinary Differential Equations Based on Inaccurate Observations. *Mathematics* **2022**, *10*, 502. https://doi.org/10.3390/ math10030502

Academic Editors: Alexandru Agapie, Denis Enachescu, Vlad Stefan Barbu and Bogdan Iftimie

Received: 13 January 2022 Accepted: 2 February 2022 Published: 4 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

of linear regression analysis is used in relation to a regression function that slightly deviates from the original function in a small neighbourhood of some time moment [8–13].

This method is based on minimizing the standard deviation of a sequence of observations from a linear regression function. In this case, such a relationship is selected between the number of observations and the interval between neighbouring observations so that the resulting error in determining the parameters tends to zero when the number of observations tends to infinity.

The final stage of the parameter estimation algorithm is the substitution of estimates of the values of functions and the values of their derivatives into the original system of equations at the selected point. Further, by analogy with the method of moments, unknown parameters of the system of equations are estimated and the consistency of the estimates obtained is proved. This paper also uses the implicit function theorem, which allows us to establish that the obtained parameter estimates are consistent depending on the number of observations. Based on the results obtained, computational experiments were carried out.

Thus, elements of system analysis have been introduced into the solution of the task. Here, an algorithm for estimating parameters acts as a system. The error of the obtained estimate is an indicator of its quality. A sequence of inaccurate measurements is a random process. Furthermore, the process and the method of linear regression analysis applied to an almost linear regression function is used as an optimization procedure. It is evaluated using the theorem on the existence and uniqueness of the solution of a system of ordinary differential equations and with the help of the implicit function theorem. Additionally, known error estimates in the linear regression analysis method are used.

#### **2. Estimating the Coefficients of a System of Ordinary Differential Equations by Inaccurate Observations**

#### *2.1. Preliminaries*

Consider a system of ordinary differential equations with fixed values of parameters *βi* = *β*0 *i* , *i* = 1, . . . , *m*,

$$\frac{d\mathbf{x}\_i}{dt} = F\_i(\mathbf{x}\_1, \dots, \mathbf{x}\_{m\_i} \beta\_1^0, \dots, \beta\_m^0), \ i = 1, \dots, m\_i \tag{1}$$

where *x*1 = *x*1(*t*), ... , *xm* = *xm*(*t*) are unknown functions. In well-known monographs on the theory of ordinary differential equations (see, for example, [14,15]), the theorem of the existence and uniqueness of the solution of this system in a small neighbourhood of a certain point is formulated and proved in Theorem 1.

**Theorem 1.** *Assume that functions Fi* = *Fi*(*x*1, ... , *xm*, *β*0 1, ... , *β*0 *m*) *are continuous in a rectangular parallelepiped Q* = {(*x*1, ..., *xm*) ∈ *Rm* : *x*0 *i* − *ai* ≤ *xi* ≤ *x*0 *i* + *ai*, *i* = 1, .., *m*} *together with their partial derivatives ∂Fi ∂xi* , *i* = 1, ..., *m*. *Then there is a segment t*0 − *r* ≤ *t* ≤ *t*0 + *r*, *on which the system of Equation (1) has a unique solution satisfying the initial conditions xi*(*t*0) = *x*0 *i* , *i* = 1, . . . , *m*.

**Remark 1.** *From the Weierstrass theorem for continuous functions on a compact, it follows that the functions xi*(*t*), *i* = 1, ..., *m*, *on the segment* [*t*0 − *r*, *t*0 + *r*] *(continuity follows from differentiability) and function* % % % % *Fi* · *∂Fi ∂xi* % % % % , *i* = 1, ..., *m*, *on a set Q (due to the continuity of the multipliers) reach their highest final values Ci*.

Denote

$$M\_0 = (\mathbf{x}\_{1'}^0, \dots, \mathbf{x}\_{m'}^0 \beta\_1^0, \dots, \beta\_m^0), \ F\_i(M\_0) = F\_{i'}^0, \ i = 1, \dots, m,$$

$$M\_0' = (\mathbf{x}\_{1'}^0, \dots, \mathbf{x}\_{m'}^0 \ F\_{1'}^0, \dots, F\_{m'}^0 \ \beta\_{1'}^0, \dots, \beta\_m^0),$$

$$F\_i(\mathbf{x}\_1, \dots, \mathbf{x}\_{\mathfrak{m}}, f\_1, \dots, f\_{\mathfrak{m}}, \beta\_1, \dots, \beta\_{\mathfrak{m}}) = F\_i(\mathbf{x}\_1, \dots, \mathbf{x}\_{\mathfrak{m}}, \beta\_1, \dots, \beta\_{\mathfrak{m}}) - f\_{i\nu}$$

where *Fi* are described in Theorem 1, and consider the system of equations

$$G\_i(\mathbf{x}\_1, \dots, \mathbf{x}\_m, f\_1, \dots, f\_m, \beta\_1, \dots, \beta\_m) = 0, \ i = 1, \dots, m. \tag{2}$$

In monographs on mathematical analysis (see, for example, [16,17]), conditions are formulated, under which the system (2) may be resolved with respect to variables *β*1, ... , *βm* (see for example Theorem 2).

**Theorem 2.** *If the functions Gi*, *i* = 1, ... , *m are continuously differentiable in the neighbourhood of the point M* 0 *and the Jacobian*

$$\frac{\partial \left( G\_{1}, \ldots, G\_{m} \right)}{\partial \left( \beta\_{1}, \ldots, \beta\_{m} \right)} \Big|\_{M\_{0}'} \neq 0,\tag{3}$$

*then there are neighbourhoods U*, *V*, *W of points* (*x*0 1, ... , *x*0 *m*), (*F*0 1 , ... , *F*0 *m*), (*β*0 1, ... , *β*0 *m*)*, respectively, such that the system of Equation (2) is uniquely solvable in the neighbourhood of U* × *V* × *W of the point M* 0 *relative to the variables β*1, ... , *βm*. *Moreover, if βi* = *gi*(*x*1, ... , *xm*, *f*1, ... , *fm*), *i* = 1, ... , *m*, *is the specified solution, then all functions gi are continuously differentiable in the neighbourhood U* × *V and β*0 *i* = *gi*(*x*0 1,..., *x*0 *m*, *F*0 1 ,..., *F*0 *m*).

**Remark 2.** *When the conditions of Theorem 2 are met, the functions gi*, *i* = 1, ..., *m*, *are continuous at the point* (*x*0 1,..., *x*0 *m*, *F*0 1 ,..., *F*0 *m*).

#### *2.2. Ordinary Differential Equation*

Consider the differential equation for a fixed value of the parameter *β* = *β*0

$$\frac{d\mathbf{x}}{dt} = F(\mathbf{x}, \boldsymbol{\beta}\_0) \tag{4}$$

with the initial condition *x*(0) = *x*0, assuming that the function *F*(*x*, *β*) is continuously differentiable in the neighbourhood of a point *M*0 = (*x*0, *β*0) and *∂F ∂β* % % % *M*0 = 0. Let the inaccurate observations *y*(*t*) = *x*(*t*) + *ε*(*t*) are known for the state of *x*(*t*) at the moments *t* = *kh*, *k* = 0, ±1, . . . , ±*n*, *hn* ≤ *r*. Denote

$$\varepsilon\_k = \varepsilon(kh), \ x\_k = \mathfrak{x}(kh), \ y\_k = \mathfrak{y}(kh) = \mathfrak{x}\_k + \varepsilon\_k, \ F\_0 = F(\mathfrak{x}\_0, \beta\_0).$$

and suppose that *εk*, *k* = 0, ±1, ... , ±*n*, is a set of independent and identically distributed random variables with zero mean and variance *σ*2. The problem of estimating the parameter *β*0 of the differential Equation (4) from these observations is posed.

The solution of this problem is carried out in two stages. First, they are constructed using a modification of the least squares estimation method *x*)0, *F*)0 and their convergence to the estimated parameters *x*0, *F*0 is investigated. Then, by analogy with the method of moments, an estimate of *β* )0 is constructed and its convergence to the estimated parameter *β*0 is investigated.

**Evaluation of values** *x*0, *F*0**.** Let us introduce the notations, outlining the method for defining *x*)0, *F*)0

$$
\widehat{\mathbf{x}}\_0 = \frac{\sum\_{k=-n}^n y\_k}{2n+1}, \widehat{\mathbf{F}}\_0 = \frac{\sum\_{k=-n}^n y\_k k h}{\sum\_{k=-n}^n (kh)^2}.\tag{5}
$$

**Theorem 3.** *If σ*2 < ∞ *and h* = *n*−*α*, *then, for α* > 1*, the estimate of x*)0 *is an asymptotically unbiased and consistent estimate of the parameter x*0. *The estimate F*)0 *is an asymptotically unbiased estimate of the parameter F*0. *At* 1 < *α* < 3/2*; the estimate F*)0 *is a consistent estimate of F*0.

**Proof of Theorem 3.** Denote *y*˜*k* = *x*0 + *F*0*kh* + *εk* and put

$$
\bar{\mathfrak{x}}\_0 = \frac{\sum\_{k=-n}^n \mathfrak{y}\_k}{2n+1}, \\
\bar{F}\_0 = \frac{\sum\_{k=-n}^n \mathfrak{y}\_k k h}{\sum\_{k=-n}^n (k h)^2}.
$$

Estimates of *x*˜0, *F*˜ 0 are obtained by the least squares method for coefficients *x*0, *F*0 of linear regression [9] and satisfy the following relations

$$E\tilde{\mathbf{x}}\_0 = \mathbf{x}\_0, \ E\tilde{\mathbf{F}}\_0 = \mathbf{F}\_0, \ Var\tilde{\mathbf{x}}\_0 = \frac{\sigma^2}{2n+1}, \ Var\tilde{\mathbf{F}}\_0 = \frac{\sigma^2}{\sum\_{k=-n}^n (kh)^2}. \tag{6}$$

Here, *Ex* is mathematical expectation of arbitrary random variable *x* and *Varx* = *E*(*x* − *Ex*)2 is its variance. In turn, the following equalities are almost certainly fulfilled

$$\hat{\mathfrak{x}}\_{0} - \mathfrak{x}\_{0} = \frac{\sum\_{k=-n}^{n} (\hat{y}\_{k} - \tilde{y}\_{k})}{2n+1}, \\ \hat{F}\_{0} - \bar{F}\_{0} = \frac{\sum\_{k=-n}^{n} (\hat{y}\_{k} - \tilde{y}\_{k})kh}{\sum\_{k=-n}^{n} (kh)^{2}}. \tag{7}$$

Moreover, the differences *y*)*k* − *y*˜*k* = *xk* − *x*0 − *F*0*kh*, *k* = 0, ±1, ... , ±*n* are deterministic quantities.

The Remark 1 implies the existence of a number *C* satisfying the inequality

$$\sup\_{||t|| \le nh} |\mathfrak{x}^{\prime\prime}(t)| = \sup\_{|t| \le nh} \left| \frac{\partial F(\mathfrak{x}(t), \beta\_0)}{\partial \mathfrak{x}} F(\mathfrak{x}(t), \beta\_0) \right| = 2\mathcal{C} < \infty.$$

Then, from the Taylor formula with a residual term in the Lagrange form,

$$\mathbf{x}(kh) = \mathbf{x}(0) + F\_0kh + \frac{(kh)^2}{2} \mathbf{x}^{\prime\prime}(kh\tau\_k), \ |\tau\_k| \le 1, \ k = 0, \pm 1, \dots, \pm n\_{\text{tot}}$$

inequalities follow

$$|\mathbf{x}\_k - \mathbf{x}\_0 - F\_0kh| \le \mathcal{C}(kh)^2, \ k = 0, \pm 1, \dots, \pm n. \tag{8}$$

From the Formulas (7) and (8) for *n* → ∞, the relations follow

$$\left| \left| \hat{\mathbf{x}}\_{0} - \hat{\mathbf{x}}\_{0} \right| \leq \frac{\sum\_{k=-n}^{n} |\mathbf{x}\_{k} - \mathbf{x}\_{0} - \mathbf{F}\_{0} \mathbf{k} \mathbf{h}|}{2n+1} \right| \leq \frac{2Ch^{2} \sum\_{k=1}^{n} k^{2}}{2n+1} \sim \frac{Ch^{2} n^{2}}{3},\tag{9}$$

$$|\vec{F}\_0 - F\_0| \le \frac{\sum\_{k=-n}^{n} |(\mathbf{x}\_k - \mathbf{x}\_0 - F\_0kh)kh|}{\sum\_{k=-n}^{n} (kh)^2} \le \frac{Ch^3 \sum\_{k=1}^{n} k^3}{\sum\_{k=1}^{n} h^2 k^2} \sim \frac{Chn}{4}.\tag{10}$$

The Formulas (6), (9) and (10) lead to the relations

$$|E\hat{\mathbf{x}}\_0 - \mathbf{x}\_0| = |E\hat{\mathbf{x}}\_0 - E\hat{\mathbf{x}}\_0| \preceq \frac{Ch^2n^2}{2}, \text{ } Var\hat{\mathbf{x}}\_0 = Var\hat{\mathbf{x}}\_0. \tag{11}$$

$$|E\widehat{F}\_0 - F\_0| = |E\widehat{F}\_0 - EF\_0| \preceq \frac{3\mathcal{C}\hbar n}{4}, \; Var\widehat{F}\_0 = Var\widehat{F}\_0. \tag{12}$$

Here *an* % *bn* means that lim sup *n*→∞ *an*/*bn* ≤ 1. Then from the condition *h* = *n*−*α*, *α* > 1, and the Relations (11) and (12) we have

$$|E\hat{x}\_0 - x\_0| \to 0, \ |E\hat{x}\_0' - F\_0| \to 0, \ n \to \infty,\tag{13}$$

that *x*)0, *F*)0 are asymptotic unbiased estimates of *x*0, *F*0.

From the Bieneme–Chebyshev inequality, the Relations (9) and (11) and the conditions *h* = *n*−*α*, *α* > 1, we get for any *δ* > 0

$$P(|\hat{\mathbf{x}}\_0 - \mathbf{x}\_0| > \delta) \le P(|\hat{\mathbf{x}}\_0 - \tilde{\mathbf{x}}\_0| + |\tilde{\mathbf{x}}\_0 - \mathbf{x}\_0|) > \delta) = P(|\tilde{\mathbf{x}}\_0 - \mathbf{x}\_0| \ge \delta - |\hat{\mathbf{x}}\_0 - \tilde{\mathbf{x}}\_0|) \le \delta$$

$$\le \frac{\sigma^2}{(2n + 1)(\delta - |\hat{\mathbf{x}}\_0 - \tilde{\mathbf{x}}\_0|)^2} \to 0, \ n \to \infty. \tag{14}$$

Thus, for *h* = *n*−*α*, *α* > 1, estimate *x*)0 is a consistent estimate of *x*0.

At the same time, from the Relations (10), (12) and (13) for *h* = *n*−*α*, 1 < *α* < 3/2, we get for any *δ* > 0

$$P(|\hat{F}\_0 - F\_0| > \delta) \le P(|\hat{F}\_0 - \bar{F}\_0| + |\bar{F}\_0 - F\_0| > \delta)) = $$

$$P(|\bar{F}\_0 - F\_0| > \delta - |\hat{F}\_0 - \bar{F}\_0|) \le \frac{3\sigma^2}{h^2 n^3 (\delta - |\hat{F}\_0 - \bar{F}\_0|)^2} \to 0, \ n \to \infty. \tag{15}$$

Therefore, if the condition *h* = *n*−*α*, 1 < *α* < 3/2, is true, the estimate *F*)0 is a consistent estimate of *F*0.

**Remark 3.** *It is worth noting that Theorem 3 is true for any distribution of random variables εk with finite variance σ*2. *Indeed it is necessary to prove limit relation Hn* = ∑*n k*=−*n εkk h* ∑*n k*=−*n k*2 → 0, *n* → ∞. *However, the most reasonable way to solve this question is to consider such distributions of random variables εk as normal for σ*2 < ∞/ *or stable for σ*2 = ∞, *because Hn has normal/stable distribution also.*

**Evaluation of parameter** *β*0. Consider the equation

$$F(\widehat{x}\_0, \boldsymbol{\beta}) = \widehat{F}\_0. \tag{16}$$

**Theorem 4.** *In conditions of Theorem 3, Equation (16) has a unique solution β* )0, *which is a consistent estimate of the parameter β*0.

**Proof of Theorem 4.** Since the function *F*(*x*, *β*) is continuously differentiable in the neighbourhood of the point *M*0 = (*x*0, *β*0) and *∂F ∂β* % % % *M*0 = 0,, then the conditions of the theorem for the function *G*(*x*, *f* , *β*) = *F*(*x*, *β*) − *f* are fulfilled. So, in some neighbourhood of the point *M* 0 = (*x*0, *F*0, *β*0), the equation is solvable with respect to *β* = *g*(*x*, *f*), while *β*0 = *g*(*x*0, *F*0). Then, from the Remark 2, we get that for any *ε* > 0 there exists *δ*(*ε*) > 0 such that in the neighbourhood {(*x*, *f*) : |*x* − *x*0| ≤ *δ*(*ε*), | *f* − *F*0| ≤ *δ*(*ε*)} of the point (*x*0, *F*0) the inequality |*β* − *β*0| ≤ *ε* is executed.

It follows that with the specified choice of *δ*(*ε*), the relation is fulfilled

$$|\hat{x}\_0 - x\_0| \le \delta(\varepsilon), \ |\hat{F}\_0 - F\_0| \le \delta(\varepsilon) \implies |\hat{\beta}\_0 - \beta\_0| \le \varepsilon. \tag{17}$$

In turn, from the Relations (14) and (15) it follows that for any *ε* and *δ*(*ε*) there is such a *n*0(*ε*, *δ*(*ε*)), that for any *n* > *n*0(*ε*, *δ*(*ε*)) inequalities are fair

$$P(|\hat{x}\_0 - x\_0| \le \delta(\varepsilon)) \ge 1 - \frac{\varepsilon}{2},\ P(|\hat{F}\_0 - F\_0| \le \delta(\varepsilon)) \ge 1 - \frac{\varepsilon}{2}.\tag{18}$$

Therefore, from the Relations (17) and (18) we have *P*(|*β* )0 − *β*0| ≤ *ε*) ≥ 1 − *ε*. Thus, for any *ε* > 0, there exists *n*0(*ε*) such that for *n* > *n*0(*ε*), the inequality holds *P*(|*β* )0 − *β*0| > *ε*) < *ε*, which means consistency (convergence in probability at *n* → ∞) of the constructed estimate.

#### *2.3. System of Differential Equations*

Consider a system (1) with initial conditions *xi*(0) = *x*0 *i* , *i* = 1, ... , *m*. We assume that the functions *Fi*(*x*1, ... , *xm*, *β*1, ... , *βm*), *i* = 1, ... , *m*, are continuously differentiable in the neighbourhood of the point *M*0 and the Jacobian *∂*(*F*1,..., *Fm*) *∂*(*β*1,..., *βm*) % % % *M*0 = 0. Inaccurate observations are known *yi*(*t*) = *xi*(*t*) + *εi*(*t*) for the state *xi*(*t*), *i* = 1, ... , *m*, at moments *t* = *kh*, *k* = 0, ±1, ... , ±*n*, *hn* ≤ *r*. Let *εi*(*kh*), *k* = 0, ±1, ... , ±*nh*, *i* = 1, ... , *m*, is a set of independent and identically distributed random variables with zero mean and variance *σ*2. The task is to estimate the vector of parameters (*β*0 1, ... , *β*0 *m*) of a system of differential Equation (1) based on these observations.

Denote

$$\hat{x}\_i^0 = \frac{\sum\_{k=-n}^n y\_i(kh)}{2n+1}, \ \hat{F}\_i^0 = \frac{\sum\_{k=-n}^n y\_i(kh)}{\sum\_{k=-n}^n (kh)^2}, \ i = 1, \ldots, m. \tag{19}$$

**Theorem 5.** *If σ*2 < ∞ *and h* = *n*−*α*, *then, for α* > 1*, the estimate x*)0 *i is an asymptotically unbiased and consistent estimate of the parameter x*0 *i* . *The estimate F*)0 *i is an asymptotically unbiased estimate of the parameter F*0 *i* . *For* 1 < *α* < 3/2*, the estimate F*)0 *i is a consistent estimate of the value F*0 *i* , *i* = 1, . . . , *m*.

Consider a system of equations

$$F\_{\hat{\mathbf{i}}}(\hat{\mathbf{x}}\_1, \dots, \hat{\mathbf{x}}\_m, \beta\_1, \dots, \beta\_m) = F\_{\hat{\mathbf{i}}}^0, \mathbf{i} = 1, \dots, m. \tag{20}$$

**Theorem 6.** *In conditions of Theorem 5 the system of Equation (20) has a unique solution* (*β* )0 1,... *β* )0 *m*), *which is a consistent estimate of the vector of parameters* (*β*0 1,..., *β*0 *m*).

The proofs of the Theorems 5 and 6 almost verbatim repeat the proofs of the Theorems 3 and 4.

**Remark 4.** *Theorems 3–6 are devoted to ordinary differential equations of the first order and their systems. However, it is possible to spread them to ordinary differential equations and their systems of arbitrary order. For this purpose it is possible to use for examples results of [12,13].*

#### *2.4. Computational Experiment*

**Example 1.** *The computational experiment was conducted first for the Cauchy problem*

$$\frac{d\mathbf{x}}{dt} = F(\mathbf{x}, \boldsymbol{\beta}\_0) = \boldsymbol{\beta}\_0 \mathbf{x}, \; \mathbf{x}(0) = 1, \; \boldsymbol{\beta}\_0 = 0.5.$$

*The solution of this equation has the form x* = *eb*0*t* . *We assumed that by observing the process described by this equation,* ±*kh*, *k* = 0, 1, ..., *n*, *h* = *n*−5/4, *n* = 10, 000, *inaccurate observations were obtained at time points y*±*k* = *e*±*b*0*hk* + *ε*±*k*, *k* = 0, 1, ..., *n*.

Here, independent random variables *ε*±*k*, *k* = 0, 1, ..., *n*, are distributed uniformly on the segment [−1/2, 1/2] left and on the segment [−1/4, 1/4] right. According to the Formula (5), the parameters *x*0, *F*0 = *F*(*x*0, *β*0) in our notation *x*)0, *F*)0, were evaluated first, then the formula for evaluating the parameter *β*0 was found from the equation *F*)0 = *x*)0*β* )0. Table 1 shows the results of a computational experiment conducted 1000 times, namely, the interval distribution (5 intervals) of relative frequencies *β* )0.

**Table 1.** Interval distribution of estimate *β* )0 when *ε*±*k* has variance *σ*2 = 1/12 left and variance *σ*2 = 1/48 right.

Consequently, a decrease in variance *σ*2 improves the quality of the obtained estimates sufficiently clearly.

Now, consider the case in which independent random variables *ε*±*k*, *k* = 0, 1, ..., *n*, are distributed normally with mean 0 and variance *σ*2. Table 2 shows the results of a computational experiment conducted 1000 times, namely, the interval distribution (five intervals) of relative frequencies *β* )0.

**Table 2.** Interval distribution of estimate *β* )0 when *ε*±*k* has variance *σ*2 = 1/12 left and variance *σ*2 = 1/48 right.

Consequently, the quality of obtained results for disturbances distributed normally behaves like in a case of uniform distribution.

**Example 2.** A computational experiment was also carried out for the system of Lorentz equations

$$\begin{cases} \frac{dx}{dt} = F\_1(x, y, z, \sigma\_0, r\_0, b\_0) = \sigma\_0(y - x),\\ \frac{dy}{dt} = F\_2(x, y, z, \sigma\_0, r\_0, b\_0) = x(r\_0 - z) - y,\\ \frac{dz}{dt} = F\_3(x, y, z, \sigma\_0, r\_0, b\_0) = xy - b\_0 z, \end{cases} \tag{21}$$

with the given initial conditions *x*(0) = 1, *y*(0) = 2, *z*(0) = 1, in the case of *σ*0 = 1, *r*0 = 2, *b*0 = 3. The solution of this system is not written out explicitly, but it is solved by the finite difference method. We write out the corresponding equations for the grid {±*kh*, *k* = 0, 1, ..., *n*} in increments of *h* = *n*−5/4, *n* = 10,000:

$$\begin{cases} \mathbf{x}\_{\pm(k+1)} = \mathbf{x}\_{\pm k} \pm \sigma\_0 h(y\_{\pm k} - \mathbf{x}\_{\pm k}),\\ \mathbf{y}\_{\pm(k+1)} = y\_{\pm k} \pm h(\mathbf{x}\_{\pm k}(r\_0 - z\_{\pm k}) - y\_{\pm k}),\\ z\_{\pm(k+1)} = z\_{\pm k} \pm h(\mathbf{x}\_{\pm k} y\_{\pm k} - b\_0 z\_{\pm k}), \end{cases} \tag{22}$$

*x*0 = 1, *y*0 = 2, *z*0 = 1. We assumed that by observing the process described by these equations, inaccurate observations were obtained

$$X\_{\pm k} = \mathfrak{x}\_{\pm k} + \varepsilon\_1(\pm hk), \ \ Y\_{\pm k} = \mathfrak{y}\_{\pm k} + \varepsilon\_2(\pm hk), \ Z\_{\pm k} = z\_{\pm k} + \varepsilon\_3(\pm hk), \ k = 0, 1, \ldots, n\_r$$

where *εi*(±*hk*), *i* = 1, 2, 3, *k* = 0, 1, ... , *n*, are independent random variables, distributed uniformly over a segment [−1/2, 1/2]. According to the Formula (19), the parameters were evaluated first *x*0, *y*0, *z*0, *F*0 *i* = *Fi*(*x*0, *y*0, *z*0, *σ*0,*r*0, *b*0), *i* = 1, 2, 3, in our notation

*x*)0, *y*)0, )*z*0, *F*)0 *i* , *i* = 1, 2, 3. Further, the estimates of the parameters *σ*0,*r*0, *b*0 were found from the relations

$$\begin{cases} \tilde{F}\_1^0 = \sigma(\hat{y}\_0 - \hat{x}\_0), \\ \tilde{F}\_2^0 = \hat{x}\_0(r - \hat{z}\_0) - \hat{y}\_0, \\ \tilde{F}\_3^0 = \hat{x}\_0 \hat{y}\_0 - b\hat{z}\_0. \end{cases} \tag{23}$$

Table 3 shows the results of a computational experiment conducted 1000 times, namely, the interval distribution of relative frequencies )*σ*0, )*r*0, *β* )0.

**Table 3.** Interval distribution of estimates )*σ*0, )*r*0, )*b*0.

#### **3. Conclusions**

Remarks 3 and 4 indicate the following possible generalizations of the results obtained in Theorems 3–6. First, we should consider the case when the variance of random perturbations *σ*2 decreases and so quality of obtained estimates improves. However, if the variance *σ*2 = ∞ like in a case of heavy tails of disturbances distributions, then it is necessary to consider stable distribution of random variables *εk*. Secondly, we should consider the case of ordinary differential Equations (and systems) of higher than the unit order.

Furthermore, at last, along with systems of ordinary differential equations, the proposed method for estimating parameters may be applied to equations or systems of partial differential equations. As a basis for the development of this method of parameter estimation, the theorem of the existence of a solution of partial differential equations system in the vicinity of a certain point may be taken (see, for example, [18]).

#### **4. Discussion**

The solution of the considered problem involves the choice of an experimental plan, the use of the theorem of existence and uniqueness for a system of ordinary differential equations, the implicit function theorem and the method of linear regression analysis. Linear regression analysis is based on minimizing of the root-mean-square deviation of the sequence of observations from the linear regression function.

Practically, all the considered generalizations of the results obtained in the paper arise at the junction of several scientific directions. These include probability theory and mathematical statistics, ordinary differential and partial differential equations and their systems, and mathematical analysis. Such tasks arising at the junction of several research directions are usually considered in the system analysis, management and information processing. This circumstance determines the choice of this research topic and the ways to solve the task and an application of optimization procedures.

**Author Contributions:** Conceptualization, G.T. and Y.K.; methodology and formal analysis, G.T.; analysis of instrument capabilities in parameter estimation, Y.K.; checking the received formulas and numerical experiments, M.O. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data supporting reported results were obtained by Yu. Kharchenko and M. Osipova.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


*Edited by Davy Vanden Broeck*

Cervical cancer is the second most prevalent cancer among women worldwide, and infection with Human Papilloma Virus (HPV) has been identified as the causal agent for this condition. The natural history of cervical cancer is characterized by slow disease progression, rendering the condition, in essence, preventable and even treatable when diagnosed in early stages. Pap smear and the recently introduced prophylactic vaccines are the most prominent prevention options, but despite the availability of these primary and secondary screening tools, the global burden of disease is unfortunately still very high. This book will focus on the clinical aspects of HPV and related disease, highlighting the latest developments in this field.

ISBN 978-953-307-860-1

From Bench to Bedside - A Clinical Perspective

Human Papillomavirus and Related Diseases -

Photo by Dr\_Microbe / iStock


Gold Catalyzed Asymmetric Transformations

*Susana Porcel García*

## **Abstract**

In this chapter, the strategies developed to attain asymmetric reactions with gold are disclosed. Because of its preferred linear arrangement, to induce asymmetry, gold(I) needs to fulfill one of the following requirements: a) the use of bulky chiral ligands, that create a chiral pocket around the active site, b) the coordination to bifunctional ligands capable to establish secondary interactions with substrates, or c) tight ion pairing with chiral counteranions. On the other hand, gold(III) profits of a square-planar coordination mode, which approaches chiral ligands to substrates. However, its tendency to be reduced leads to difficulties for its applications in catalytic asymmetric transformations. Pioneering works using cyclometaled structures, have found the balance between stability and activity, showing its potential in asymmetric transformations.

**Keywords:** gold, catalysis, bulky ligands, bifunctional ligands, ion pairing


**3. Neuromodulation applications for urological conditions**

### **3.1 Overactive bladder, urgency urinary incontinence and urgency-frequency syndromes**

Bladder overactivity manifests in a number of urinary conditions, depending on the pathophysiology and associated conditions and symptoms. Overactive bladder (OAB), defined by a compelling frequent urge to void, is not a precession of urgency urinary incontinence (UUI), nor is it a more defined form urgency-frequency syndromes: these are all an overlapping number of conditions where evidence of overactivity of the detrusor muscle may or may not be demonstrable, but is subjectively reported by patients and often objectively measurable.

The treatment for these conditions is mainly conservative and medical, be it targeting the bladder muscle or the other offending factors that lead to the overactivity, followed by intradetrusor botulinum injections, which has attained a more defined role in the OAB treatment scheme. SNM is an established mode of treatment for cases of OAB, UUI and urgency-frequency syndromes that are refractory to medical treatment, and despite arguments and established results and testaments, is yet to be designated a more primary or first line place in the treatment of these conditions [95].

### *3.1.1 Mode of effect in OAB*

It has been shown that SNM has an established modulatory effect both on micturition reflexes and higher brain centers. The SNM electrical charging of sacral roots alters neural activity, stimulating somatic afferents that signal to higher brain centers and in part restore normal control over the bladder while also inhibiting certain sensory pathways to suppress reflex bladder hyperactivity. From animal models, evidence suggests this effect is achieved through SNM's inhibition of abnormal sensory input from the pudendal nerve and neuropathological C-fibers, affecting release of μ-opiods and glutamate and suppressing bladder reflexes [16].

### *3.1.2 SNM efficacy in OAB*

Efficacy of SNM is perhaps most studied and evidently reported in refractory OAB [96]. Analysis of five trials have analytically shown significantly higher success rates for SNM in treatment of OAB compared to standard medical treatment, and equally as efficacious as intradetrusor botulinum injections with less side effects

**229**

phases [9, 104].

*Neuromodulation in Urology: Current Trends and Future Applications*

associated with the latter including risk of post-injection urinary retention and

follow-up, 83% of implants were found to have sustained efficacy [68].

**3.2 Non-obstructive urinary retention and Fowler's syndrome**

urinary retention in a "neurologically-intact" patient is.

been shown to restore normal voiding activity [98].

**3.3 Neurogenic lower urinary tract dysfunction**

*3.2.1 Mode of effect*

effect of the micturition reflex [99].

In one prospectively conducted multicenter trial on OAB patients, the 5-year success rate of SNM was 67%, with the most common adverse event or reason for failure demonstrated to be an undesirable change in stimulation, followed by site pain and ineffectiveness of treatment [97]. The InSite trial reported on one of the longest prospective follow-ups for SNM implants for refractory OAB. At 36-months

Non-obstructive urinary retention (NOUR) is one of the main indications for SNM therapy. It denotes an unidentifiable mechanical cause that may obstruct urinary outflow from the urinary bladder, resulting in urinary retention. It may be the result of an established neurological disease, as is the case in the acute phase of spinal shock after spinal cord trauma, or in a minority of MS patients. Neither of these conditions are indicated for SNM treatment. However, chronic or recurrent

One form of NOUR is termed Fowler's syndrome after the neurophysiologist Professor Clare J. Fowler who first described it in 1985. It is a cluster of symptoms and findings identified in a typically young woman with unexplained urinary retention, increased electromyographic activity of the external urinary sphincter and its failure to relax, and some associations to other female syndromes have been described including polycystic ovaries. Application of SNM in these patients has

Researchers have used a number of animal models to establish the mode of effect SNM exerts in NOUR and Fowler's. Basic science evidence suggests that by blocking the inhibitory effect that abnormal afferent activity from the external urethral and anal sphincters has on micturition, restoration of the ability of the patient to void occurs. This stimulation is through blockade of the pudendal nerve's stimulatory

Lower urinary tract symptoms resulting from neurological disease are varied, and thus, determination of these symptoms and assessment is necessary before consideration for neuromodulation as not all symptoms would be ideally treated using this modality. Neurological diseases that have documented voiding dysfunction elements include SCI, MS, Parkinson's disease, cerebrovascular accidents, and diabetic neuropathy. Congenital neurologic disorders such myelomeningoceles are becoming apparent causes of voiding dysfunction in adults and SNM candidates, as management of these pediatric disorders improves, and these patients grow into the adult population [9]. Previously thought to lack efficacy in neurogenic LUTD because of lack of an intact nervous system, SNM is emerging as an efficacious therapeutic modality for this population of patients especially in reducing incontinence episodes [9, 100–104]. The concept of neural remodeling as a hypothesized effect of SNM has also been visited as a potential role in neurogenic LUTD, particularly in acute spinal shock

The ICS recommends SNM as an option for control of urinary symptoms in patients with stable neurological conditions who are at a low risk of developing

*DOI: http://dx.doi.org/10.5772/intechopen.92287*

urinary tract infections [41, 74].

### *Neuromodulation in Urology: Current Trends and Future Applications DOI: http://dx.doi.org/10.5772/intechopen.92287*

associated with the latter including risk of post-injection urinary retention and urinary tract infections [41, 74].

In one prospectively conducted multicenter trial on OAB patients, the 5-year success rate of SNM was 67%, with the most common adverse event or reason for failure demonstrated to be an undesirable change in stimulation, followed by site pain and ineffectiveness of treatment [97]. The InSite trial reported on one of the longest prospective follow-ups for SNM implants for refractory OAB. At 36-months follow-up, 83% of implants were found to have sustained efficacy [68].

### **3.2 Non-obstructive urinary retention and Fowler's syndrome**

Non-obstructive urinary retention (NOUR) is one of the main indications for SNM therapy. It denotes an unidentifiable mechanical cause that may obstruct urinary outflow from the urinary bladder, resulting in urinary retention. It may be the result of an established neurological disease, as is the case in the acute phase of spinal shock after spinal cord trauma, or in a minority of MS patients. Neither of these conditions are indicated for SNM treatment. However, chronic or recurrent urinary retention in a "neurologically-intact" patient is.

One form of NOUR is termed Fowler's syndrome after the neurophysiologist Professor Clare J. Fowler who first described it in 1985. It is a cluster of symptoms and findings identified in a typically young woman with unexplained urinary retention, increased electromyographic activity of the external urinary sphincter and its failure to relax, and some associations to other female syndromes have been described including polycystic ovaries. Application of SNM in these patients has been shown to restore normal voiding activity [98].

### *3.2.1 Mode of effect*

*Neurostimulation and Neuromodulation in Contemporary Therapeutic Practice*

As the name suggests, this modality is applied to areas in close proximity to target internal nerves. These include the pudendal nerve, be it through transcutaneous stimulation in the vagina in a female or in the perineal region in the male, or both the pudendal and sacral nerves when applied to the sacral skin. DGN is also a form of transcutaneous electrical nerve stimulation (TENS). It is advocated as a less invasive and low-cost neuromodulation system that can also be taught to patients

Multiple small-sized trials have demonstrated improvements in symptom scores and efficacy in patients with refractory OAB or MS with bladder hyperactivity. However, although it is safe, the durability of its effect has been called into

**3.1 Overactive bladder, urgency urinary incontinence and urgency-frequency**

Bladder overactivity manifests in a number of urinary conditions, depending on the pathophysiology and associated conditions and symptoms. Overactive bladder (OAB), defined by a compelling frequent urge to void, is not a precession of urgency urinary incontinence (UUI), nor is it a more defined form urgency-frequency syndromes: these are all an overlapping number of conditions where evidence of overactivity of the detrusor muscle may or may not be demonstrable, but is subjectively

The treatment for these conditions is mainly conservative and medical, be it targeting the bladder muscle or the other offending factors that lead to the overactivity, followed by intradetrusor botulinum injections, which has attained a more defined role in the OAB treatment scheme. SNM is an established mode of treatment for cases of OAB, UUI and urgency-frequency syndromes that are refractory to medical treatment, and despite arguments and established results and testaments, is yet to be designated a more primary or first line place in the treatment of these

It has been shown that SNM has an established modulatory effect both on micturition reflexes and higher brain centers. The SNM electrical charging of sacral roots alters neural activity, stimulating somatic afferents that signal to higher brain centers and in part restore normal control over the bladder while also inhibiting certain sensory pathways to suppress reflex bladder hyperactivity. From animal models, evidence suggests this effect is achieved through SNM's inhibition of abnormal sensory input from the pudendal nerve and neuropathological C-fibers, affecting release of μ-opiods and glutamate and suppressing bladder reflexes [16].

Efficacy of SNM is perhaps most studied and evidently reported in refractory OAB [96]. Analysis of five trials have analytically shown significantly higher success rates for SNM in treatment of OAB compared to standard medical treatment, and equally as efficacious as intradetrusor botulinum injections with less side effects

**3. Neuromodulation applications for urological conditions**

reported by patients and often objectively measurable.

*2.4.3 Transcutaneous electrical nerve stimulation*

for self-application [5].

question [9, 94].

**syndromes**

conditions [95].

*3.1.1 Mode of effect in OAB*

*3.1.2 SNM efficacy in OAB*

**228**

Researchers have used a number of animal models to establish the mode of effect SNM exerts in NOUR and Fowler's. Basic science evidence suggests that by blocking the inhibitory effect that abnormal afferent activity from the external urethral and anal sphincters has on micturition, restoration of the ability of the patient to void occurs. This stimulation is through blockade of the pudendal nerve's stimulatory effect of the micturition reflex [99].

### **3.3 Neurogenic lower urinary tract dysfunction**

Lower urinary tract symptoms resulting from neurological disease are varied, and thus, determination of these symptoms and assessment is necessary before consideration for neuromodulation as not all symptoms would be ideally treated using this modality. Neurological diseases that have documented voiding dysfunction elements include SCI, MS, Parkinson's disease, cerebrovascular accidents, and diabetic neuropathy. Congenital neurologic disorders such myelomeningoceles are becoming apparent causes of voiding dysfunction in adults and SNM candidates, as management of these pediatric disorders improves, and these patients grow into the adult population [9].

Previously thought to lack efficacy in neurogenic LUTD because of lack of an intact nervous system, SNM is emerging as an efficacious therapeutic modality for this population of patients especially in reducing incontinence episodes [9, 100–104]. The concept of neural remodeling as a hypothesized effect of SNM has also been visited as a potential role in neurogenic LUTD, particularly in acute spinal shock phases [9, 104].

The ICS recommends SNM as an option for control of urinary symptoms in patients with stable neurological conditions who are at a low risk of developing upper tract deterioration from controlled voiding [4]. It is thus important to stress the need for proper assessment and continued evaluation of these patients as urinary retention, acute or chronic, could have consequences including urinary tract infection and renal failure [9].

### *3.3.1 Spinal cord injury*

SCI, especially complete transection, has long been accepted as a contraindication for sacral neuromodulation on the basis of a disturbed neural circuit. However, numerous reports have been reviewed that show promising results for SNM in the management of neurologically-stable SCI patients, even those with complete disruption [105]. In the acute phase of spinal shock where the bladder is atonic, SNM has been found to facilitate neurogenic remodeling as researchers theorize and demonstrate sustained SNM effects and remodeling in the brain [9, 104].

In a review of eight studies where SNM was employed in the management of lower urinary tract dysfunction in SCI patients, the success rate of the test phase was a shy 45%, but that later translated into a 75% success rate once the screened patients proceeded with IPG implantation. The treatment was well-tolerated and safe without any unexpected adverse events [106].

### *3.3.2 Multiple sclerosis*

MS is of special interest to neuro-urologists as the disease manifests with a spectrum of urinary symptoms and progresses with different patterns in this spectrum along the course of the disorder as well. Demyelination, the pathological hallmark of MS, eventually affects lower urinary tract nerves, resulting in dysfunction. Up to 80% of patients show neuro-urological symptoms within 10 years of diagnosis, most frequently bladder overactivity. As a matter of fact, voiding dysfunction is the first sign of the disease in up to 10% of patients [107].

Though not FDA approved, neuromodulation has been applied in MS patients for years, and its efficacy has been repeatedly demonstrated. SNM and PTNS have been shown in a number of series to decrease urinary symptoms and improve the quality of life of MS patients who demonstrate bladder overactivity; however, although SNM is approved for NOUR, it has not shown any benefit for MS patients demonstrating "hypoactive" urinary bladders with retention [107, 108].

What remains an important issue for MS patients being considered for SNM is appropriate patient counseling and communication with their treating physician or neurologist to assess the need for MRI examination in the future as well as stability of the disease, as disease progression and relapse would negatively affect the SNM outcomes [9, 105, 107].

### *3.3.3 Diabetic cystopathy*

Diabetic cystopathy is a condition that describes the neuromuscular effect long-standing diabetes has on the urinary bladder. Part of the condition stems from diabetic neuropathy, while another part may stem from vasculopathy affecting the detrusor muscle itself. In the application of neuromodulation to the control of overactivity symptoms resulting from diabetic cystopathy, promising results from series were overshadowed by a substantially higher than average rate of infections (17%) compared to the accepted average, as would be expected from any foreign body implantation in diabetic patients especially those with poor glycemic control [105].

**231**

*Neuromodulation in Urology: Current Trends and Future Applications*

Chronic pelvic pain syndrome in males and its predominantly female counterpart interstitial cystitis/bladder pain syndrome are chronic conditions of pelvic pain and voiding dysfunction with a poorly understood etiology [109]. Off-label use of SNM in the treatment of these disorders is established with significant results, and similar to its unknown etiology, the way SNM provides subjective and objective improvements in bladder pain syndrome for example is yet to be clearly defined, with obvious differences in outcomes between IC/BPS and non-IC/BPS CPPS [110]. Many theories have been suggested for this mode of effect, from restoration of balance between excitatory and inhibitory signals in the pelvic plexus at different spinal levels as well as SNM's modulatory effect on bladder function and in turn pain. Another issue for consideration is the bilateral or multiple sacral root involvement in bladder and pelvic pain disorders, thus S3 stimulation may be insufficient to providing symptomatic relief, and some researchers have demonstrated efficacy

A multitude of studies and researchers are reporting on promising results for SNM in symptomatic management of CPP disorders, demonstrating improvements in pain indices and quality of life measures particularly relating to improvements in sleep, social life and sexual activity [110, 111]. With 10% of patients of IC/BPS reaching a severe stage refractory to conservative and other modes of management, SNM has found an emerging role in the therapeutic void for this condition. Success rates of SNM in IC/BPS have been reported to be high, north of 80% in some series, with apparent and significant objective improvements in pelvic pain and specific interstitial cystitis symptom scores as well as improvements in daytime frequency,

The ICS based on grade C evidence released a level III recommendation that designates SNM as an option for patients who are deemed non-responsive to conservative treatment measures of IC/BPS and non-IC CPPS [4]. However, large randomized controlled trials are lacking, perhaps in part due to the mixed spectrum of CPP disorders, both pathologically and symptomatically, heterogenous patient population, and unclear etiologies, and poorly understood differences in outcomes

The effects of neuromodulation, particularly SNM, on improving sexual function among female patients, and male patients to an extent, are becoming more evident in the literature [114]. Dysfunction of the pudendal nerve, an important nerve in sexual stimulation, has been demonstrated in both refractory OAB and

In a cohort of female patients who received SNM implants for OAB, urgencyfrequency syndrome or NOUR, improvements in both female sexual function index and quality of life indices were reported, though they were not correlated [116].

*DOI: http://dx.doi.org/10.5772/intechopen.92287*

**3.4 Special populations and effects**

*3.4.1 CPPS and IC/BPS*

of bilateral stimulation [14].

between the disorders [109, 110, 113].

*3.4.2 Sexual function*

NOUR [115].

tion therapy [117].