*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*<sup>2</sup> ≥ 0 are the stiffness coefficients and *c*1, *c*<sup>2</sup> ≥ 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*<sup>d</sup> . *x* . The time histories of hysteretic force *F*<sup>h</sup> *x*, . *x* , stiffness force *F*s(*x*) and damping force *F*<sup>d</sup> . *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*<sup>d</sup> and total hysteretic force *F*<sup>h</sup> 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 *<sup>x</sup>* <sup>=</sup> *<sup>x</sup>*<sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>0</sup> is the relative displacement of sprung mass, *<sup>P</sup>*0(*t*) = <sup>−</sup>*m*.. *x*<sup>0</sup> is the input of system shown in Figure 5, *x*<sup>1</sup> is the absolute displacement and *x*<sup>0</sup> 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 .. *<sup>x</sup>*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*<sup>d</sup> . *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*<sup>2</sup> and bound *E*d1 = 0.5πω*X*2*c*<sup>1</sup> 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*<sup>x</sup> = 0 if *γ* = 1, β = 1; *m*<sup>x</sup> < 0 if *γ* > 1, β > 1 and *m*<sup>x</sup> > 0 if *<sup>γ</sup>* <sup>&</sup>lt; 1, <sup>β</sup> <sup>&</sup>lt; 1. It is worth noting that by assuming *<sup>m</sup>*<sup>|</sup> . x| /ω1*m*|x<sup>|</sup> ≈ 1, for all case studies considered in this work (including *γ* > 1, β < 1 or *γ* < 1, β > 1), the sign of *m*<sup>x</sup> could be predicted by determining the sign of expression *<sup>S</sup>*(*γ*, <sup>β</sup>, *<sup>ζ</sup>*1) <sup>=</sup> <sup>−</sup>-*γ*<sup>2</sup> <sup>−</sup> <sup>1</sup> + 2(β − 1)*ζ*<sup>1</sup> . , 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 <sup>0</sup>(*t*) . = σ<sup>2</sup> <sup>0</sup> 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 σ<sup>x</sup> = σx(ωe, *ζ*e) and σ . <sup>x</sup> = σ . <sup>x</sup>(ωe, *ζ*e) are the standard deviations for the solution of Equation (10).

The variance of the acceleration .. *x*<sup>1</sup> = 2*ζ*eω<sup>e</sup> . *x* + ω<sup>2</sup> <sup>e</sup>*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*<sup>e</sup> 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*<sup>0</sup> 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 *ζ*<sup>e</sup> = 0.5, and its minimum value is σ .. x1min <sup>=</sup> <sup>√</sup>2πωe*S*0. Taking *S*<sup>0</sup> = 1 m2s−<sup>3</sup> and ω<sup>e</sup> = 2πrad/s, the optimum std value of acceleration is σ.. x1min <sup>∼</sup><sup>=</sup> 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 σ<sup>2</sup> <sup>0</sup> =  ∞ −∞ *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*<sup>0</sup> = σ<sup>0</sup> *<sup>a</sup>* <sup>π</sup> , *b*<sup>1</sup> = σ<sup>0</sup> *a*(*a*2+*b*2) <sup>π</sup> , *<sup>a</sup>*<sup>0</sup> = 1, *<sup>a</sup>*<sup>1</sup> = <sup>2</sup>*a*, *<sup>a</sup>*<sup>2</sup> = *<sup>a</sup>*<sup>2</sup> + *<sup>b</sup>*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*<sup>0</sup> = 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*<sup>0</sup> of the applied harmonic input with variable frequency, as long as they are computed for the stationary regime. The numerical values *<sup>A</sup>*\* .. 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*<sup>1</sup> (ω) <sup>=</sup> <sup>√</sup>2σ*u*1/*P*0, where *<sup>u</sup>*<sup>1</sup> is the first component of the solution vector **u**. The system (30) is asymptotically stable if *a*1, *a*<sup>2</sup> > 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 *<sup>A</sup>*\* .. 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 (*ζ* <sup>1</sup> = *ζ*<sup>2</sup> = 0.5, *γ* = β = 1, *f*<sup>1</sup> = *f*<sup>2</sup> = 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 <sup>=</sup> 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* <sup>=</sup> *cji* <sup>=</sup> lim*t*→<sup>∞</sup> 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 <sup>∼</sup><sup>=</sup> <sup>σ</sup><sup>u</sup> <sup>1</sup> where <sup>σ</sup><sup>u</sup> <sup>1</sup> <sup>=</sup> <sup>√</sup>*c*11. The values of σ .. x1 obtained by using Lyapunov equation will be compared with those determined for linear equivalent system (10) where *ζ*<sup>e</sup> = 0.5 and ω<sup>e</sup> = 2π rad/s.

The values of transmissibility factors *<sup>A</sup>*\*xrms(ω*i*) <sup>=</sup> <sup>Ψ</sup>x(ω*i*) σp0*<sup>i</sup>* = <sup>√</sup><sup>2</sup> <sup>Ψ</sup>x(ω*i*) P0 , *<sup>i</sup>* = 1, *<sup>N</sup>* , and *<sup>A</sup>*\*xstd(ω*i*) <sup>=</sup> <sup>σ</sup>x(ω*i*) σp0*<sup>i</sup>* = <sup>√</sup>2σx(ω*i*) P0 , *<sup>i</sup>* = 1, *<sup>N</sup>* corresponding to rms <sup>Ψ</sup><sup>x</sup> and std <sup>σ</sup><sup>x</sup> 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, Ψ<sup>x</sup> and σx, can be approximated as <sup>Ψ</sup><sup>x</sup> <sup>∼</sup><sup>=</sup> <sup>√</sup>*c*11 rms and <sup>σ</sup><sup>x</sup> <sup>∼</sup><sup>=</sup> <sup>√</sup>*c*11 std. The values of *<sup>c</sup>*11 rms and *<sup>c</sup>*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 (*ζ*<sup>e</sup> = 0.5, ω<sup>e</sup> = 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 <sup>m</sup>|x|, *<sup>m</sup>*<sup>|</sup> . <sup>x</sup><sup>|</sup> 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 ω<sup>1</sup> = 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 <sup>=</sup> 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*<sup>11</sup> 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 <sup>=</sup> <sup>√</sup>*c*<sup>11</sup> <sup>=</sup> 0.515 ms−<sup>2</sup> and the mean value of displacement *<sup>m</sup>*xLyap <sup>∼</sup><sup>=</sup> <sup>−</sup>√*c*11 rms <sup>−</sup> *<sup>c</sup>*11 std<sup>=</sup> <sup>−</sup>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 <sup>1</sup>**


**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 *β<sup>i</sup>* = *β*<sup>0</sup> *<sup>i</sup>* , *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*<sup>1</sup> = *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*, *β*<sup>0</sup> <sup>1</sup>, ... , *<sup>β</sup>*<sup>0</sup> *<sup>m</sup>*) *are continuous in a rectangular parallelepiped <sup>Q</sup>* <sup>=</sup> {(*x*1, ..., *xm*) <sup>∈</sup> *<sup>R</sup><sup>m</sup>* : *<sup>x</sup>*<sup>0</sup> *<sup>i</sup>* <sup>−</sup> *ai* <sup>≤</sup> *xi* <sup>≤</sup> *<sup>x</sup>*<sup>0</sup> *<sup>i</sup>* + *ai*, *i* = 1, .., *m*} *together with their partial derivatives <sup>∂</sup>Fi ∂xi* , *i* = 1, ..., *m*. *Then there is a segment t*<sup>0</sup> − *r* ≤ *t* ≤ *t*<sup>0</sup> + *r*, *on which the system of Equation (1) has a unique solution satisfying the initial conditions xi*(*t*0) = *x*<sup>0</sup> *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*<sup>0</sup> − *r*, *t*<sup>0</sup> + *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, ... , *β<sup>m</sup>* (see for example Theorem 2).

**Theorem 2.** *If the functions Gi*, *i* = 1, ... , *m are continuously differentiable in the neighbourhood of the point M* <sup>0</sup> *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*<sup>0</sup> <sup>1</sup>, ... , *<sup>x</sup>*<sup>0</sup> *<sup>m</sup>*), (*F*<sup>0</sup> <sup>1</sup> , ... , *<sup>F</sup>*<sup>0</sup> *<sup>m</sup>*), (*β*<sup>0</sup> <sup>1</sup>, ... , *<sup>β</sup>*<sup>0</sup> *m*)*, respectively, such that the system of Equation (2) is uniquely solvable in the neighbourhood of U* × *V* × *W of the point M* <sup>0</sup> *relative to the variables β*1, ... , *βm*. *Moreover, if β<sup>i</sup>* = *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* <sup>×</sup> *V and <sup>β</sup>*<sup>0</sup> *<sup>i</sup>* = *gi*(*x*<sup>0</sup> <sup>1</sup>,..., *<sup>x</sup>*<sup>0</sup> *<sup>m</sup>*, *F*<sup>0</sup> <sup>1</sup> ,..., *<sup>F</sup>*<sup>0</sup> *<sup>m</sup>*).

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

#### *2.2. Ordinary Differential Equation*

Consider the differential equation for a fixed value of the parameter *β* = *β*<sup>0</sup>

$$\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 *<sup>M</sup>*<sup>0</sup> = (*x*0, *<sup>β</sup>*0) and *<sup>∂</sup><sup>F</sup> ∂β* % % % *M*<sup>0</sup> = 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 *β*<sup>0</sup> 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 *<sup>x</sup>*)0, *<sup>F</sup>*)<sup>0</sup> and their convergence to the estimated parameters *x*0, *F*<sup>0</sup> is investigated. Then, by analogy with the method of moments, an estimate of *β* )<sup>0</sup> is constructed and its convergence to the estimated parameter *β*<sup>0</sup> is investigated.

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

$$
\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 <sup>σ</sup>*<sup>2</sup> <sup>&</sup>lt; <sup>∞</sup> *and <sup>h</sup>* <sup>=</sup> *<sup>n</sup>*−*α*, *then, for <sup>α</sup>* <sup>&</sup>gt; <sup>1</sup>*, the estimate of <sup>x</sup>*)<sup>0</sup> *is an asymptotically unbiased and consistent estimate of the parameter <sup>x</sup>*0. *The estimate <sup>F</sup>*)<sup>0</sup> *is an asymptotically unbiased estimate of the parameter F*0. *At* <sup>1</sup> <sup>&</sup>lt; *<sup>α</sup>* <sup>&</sup>lt; 3/2*; the estimate <sup>F</sup>*)<sup>0</sup> *is a consistent estimate of F*0.

**Proof of Theorem 3.** Denote *y*˜*<sup>k</sup>* = *x*<sup>0</sup> + *F*0*kh* + *ε<sup>k</sup>* 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*˜ <sup>0</sup> are obtained by the least squares method for coefficients *x*0, *F*<sup>0</sup> 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* = *<sup>E</sup>*(*<sup>x</sup>* <sup>−</sup> *Ex*)<sup>2</sup> 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 *<sup>y</sup>*)*<sup>k</sup>* <sup>−</sup> *<sup>y</sup>*˜*<sup>k</sup>* <sup>=</sup> *xk* <sup>−</sup> *<sup>x</sup>*<sup>0</sup> <sup>−</sup> *<sup>F</sup>*0*kh*, *<sup>k</sup>* <sup>=</sup> 0, <sup>±</sup>1, ... , <sup>±</sup>*<sup>n</sup>* 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* <sup>≤</sup> 1. Then from the condition *<sup>h</sup>* <sup>=</sup> *<sup>n</sup>*−*α*, *<sup>α</sup>* <sup>&</sup>gt; 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 *<sup>x</sup>*)0, *<sup>F</sup>*)<sup>0</sup> are asymptotic unbiased estimates of *<sup>x</sup>*0, *<sup>F</sup>*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 *<sup>h</sup>* <sup>=</sup> *<sup>n</sup>*−*α*, *<sup>α</sup>* <sup>&</sup>gt; 1, estimate *<sup>x</sup>*)<sup>0</sup> is a consistent estimate of *<sup>x</sup>*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 *<sup>h</sup>* <sup>=</sup> *<sup>n</sup>*−*α*, 1 <sup>&</sup>lt; *<sup>α</sup>* <sup>&</sup>lt; 3/2, is true, the estimate *<sup>F</sup>*)<sup>0</sup> is a consistent estimate of *F*0.

**Remark 3.** *It is worth noting that Theorem 3 is true for any distribution of random variables ε<sup>k</sup> with finite variance <sup>σ</sup>*2. *Indeed it is necessary to prove limit relation Hn* <sup>=</sup> <sup>∑</sup>*<sup>n</sup> <sup>k</sup>*=−*<sup>n</sup> <sup>ε</sup>kk h* ∑*<sup>n</sup> <sup>k</sup>*=−*<sup>n</sup> <sup>k</sup>*<sup>2</sup> <sup>→</sup> 0, *<sup>n</sup>* <sup>→</sup> <sup>∞</sup>. *However, the most reasonable way to solve this question is to consider such distributions of random variables ε<sup>k</sup> as normal for σ*<sup>2</sup> < ∞/ *or stable for σ*<sup>2</sup> = ∞, *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 *<sup>M</sup>*<sup>0</sup> = (*x*0, *<sup>β</sup>*0) and *<sup>∂</sup><sup>F</sup> ∂β* % % % *M*<sup>0</sup> = 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* <sup>0</sup> = (*x*0, *F*0, *β*0), the equation is solvable with respect to *β* = *g*(*x*, *f*), while *β*<sup>0</sup> = *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*(|*β* )<sup>0</sup> <sup>−</sup> *<sup>β</sup>*0| ≤ *<sup>ε</sup>*) <sup>≥</sup> <sup>1</sup> <sup>−</sup> *<sup>ε</sup>*. Thus, for any *ε* > 0, there exists *n*0(*ε*) such that for *n* > *n*0(*ε*), the inequality holds *P*(|*β* )<sup>0</sup> <sup>−</sup> *<sup>β</sup>*0<sup>|</sup> <sup>&</sup>gt; *ε*) < *ε*, 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*<sup>0</sup> *<sup>i</sup>* , *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 *<sup>M</sup>*<sup>0</sup> and the Jacobian *<sup>∂</sup>*(*F*1,..., *Fm*) *∂*(*β*1,..., *βm*) % % % *M*<sup>0</sup> = 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 (*β*<sup>0</sup> <sup>1</sup>, ... , *<sup>β</sup>*<sup>0</sup> *<sup>m</sup>*) 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 <sup>σ</sup>*<sup>2</sup> <sup>&</sup>lt; <sup>∞</sup> *and <sup>h</sup>* <sup>=</sup> *<sup>n</sup>*−*α*, *then, for <sup>α</sup>* <sup>&</sup>gt; <sup>1</sup>*, the estimate <sup>x</sup>*)<sup>0</sup> *<sup>i</sup> is an asymptotically unbiased and consistent estimate of the parameter x*<sup>0</sup> *<sup>i</sup>* . *The estimate <sup>F</sup>*)<sup>0</sup> *<sup>i</sup> is an asymptotically unbiased estimate of the parameter F*<sup>0</sup> *<sup>i</sup>* . *For* <sup>1</sup> <sup>&</sup>lt; *<sup>α</sup>* <sup>&</sup>lt; 3/2*, the estimate <sup>F</sup>*)<sup>0</sup> *<sup>i</sup> is a consistent estimate of the value F*0 *<sup>i</sup>* , *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 <sup>1</sup>,... *β* )0 *<sup>m</sup>*), *which is a consistent estimate of the vector of parameters* (*β*<sup>0</sup> <sup>1</sup>,..., *<sup>β</sup>*<sup>0</sup> *<sup>m</sup>*).

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*<sup>t</sup>* . *We assumed that by observing the process described by this equation,* <sup>±</sup>*kh*, *<sup>k</sup>* <sup>=</sup> 0, 1, ..., *<sup>n</sup>*, *<sup>h</sup>* <sup>=</sup> *<sup>n</sup>*−5/4, *<sup>n</sup>* <sup>=</sup> 10, 000, *inaccurate observations were obtained at time points y*±*<sup>k</sup>* <sup>=</sup> *<sup>e</sup>*±*b*0*hk* <sup>+</sup> *<sup>ε</sup>*±*k*, *<sup>k</sup>* <sup>=</sup> 0, 1, ..., *<sup>n</sup>*.

Here, independent random variables *<sup>ε</sup>*±*k*, *<sup>k</sup>* = 0, 1, ..., *<sup>n</sup>*, 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 *<sup>x</sup>*0, *<sup>F</sup>*<sup>0</sup> <sup>=</sup> *<sup>F</sup>*(*x*0, *<sup>β</sup>*0) in our notation *<sup>x</sup>*)0, *<sup>F</sup>*)0, were evaluated first, then the formula for evaluating the parameter *<sup>β</sup>*<sup>0</sup> was found from the equation *<sup>F</sup>*)<sup>0</sup> <sup>=</sup> *<sup>x</sup>*)0*<sup>β</sup>* )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 *β* )<sup>0</sup> when *<sup>ε</sup>*±*<sup>k</sup>* has variance *<sup>σ</sup>*<sup>2</sup> <sup>=</sup> 1/12 left and variance *σ*<sup>2</sup> = 1/48 right.

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

Now, consider the case in which independent random variables *<sup>ε</sup>*±*k*, *<sup>k</sup>* = 0, 1, ..., *<sup>n</sup>*, 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 *β* )<sup>0</sup> when *<sup>ε</sup>*±*<sup>k</sup>* has variance *<sup>σ</sup>*<sup>2</sup> <sup>=</sup> 1/12 left and variance *σ*<sup>2</sup> = 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 *σ*<sup>0</sup> = 1, *r*<sup>0</sup> = 2, *b*<sup>0</sup> = 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*, *<sup>k</sup>* <sup>=</sup> 0, 1, ..., *<sup>n</sup>*} in increments of *<sup>h</sup>* <sup>=</sup> *<sup>n</sup>*−5/4, *<sup>n</sup>* <sup>=</sup> 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*<sup>0</sup> = 1, *y*<sup>0</sup> = 2, *z*<sup>0</sup> = 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*<sup>0</sup> *<sup>i</sup>* = *Fi*(*x*0, *y*0, *z*0, *σ*0,*r*0, *b*0), *i* = 1, 2, 3, in our notation

*<sup>x</sup>*)0, *<sup>y</sup>*)0, )*z*0, *<sup>F</sup>*)<sup>0</sup> *<sup>i</sup>* , *i* = 1, 2, 3. Further, the estimates of the parameters *σ*0,*r*0, *b*<sup>0</sup> 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, *<sup>β</sup>* )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 *σ*<sup>2</sup> decreases and so quality of obtained estimates improves. However, if the variance *σ*<sup>2</sup> = ∞ 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**

