Vibrations in Single-Degrees-of-Freedom Sampled-Data Linear Mechanical Systems

This paper aims to present that the effect of sampling can result in multi-frequency vibration even in the case of a single-degree-of-freedom linear mechanical model. Even though the sampled-data systems have an infinite number of characteristic exponents due to sampling, the vibrations of these systems can still be characterized by an effective system model with a single dominant frequency. However, as this paper shows, additional harmonics become relevant, resulting in multi-frequency vibrations depending on the magnitude of applied control parameters. The vibrations presented by the time histories of vibration and their spectra resulted in numerical simulation of the sampled-data system.


Introduction
Sampled-data systems play a fundamental role in engineering due to the digital realization of the applied controller. In these systems, the plant, i.e., the controlled physical system, operates in continuous time while the controller is realized in discrete-time [1].
In sampled-data systems, the applied periodic sampling results in periodicity in the closed-loop system, and these can be analyzed in the discrete-time [2,3]. In the past decades, comprehensive studies and textbooks are available on digital control theory, just mentioned a few [4][5][6][7][8][9].
Usually, phase or gain margin [9] is used to design the controllers, while it is also common to use another design method called pole placement [4]. Based on this method, the feedback gains are determined to modify the poles or the characteristic roots of the controlled plant to achieve the desired dynamics.
In the case of discrete-time control, there is an additional transformation. The characteristic roots related to the desired dynamics are usually described in continuous time. The characteristic roots are transformed to discrete-time with a given sampling time to get the control gains for the discrete-time model. These roots are often called the characteristic multipliers [10].
The analysis of sampling effect requires the handling of infinite-dimensional mathematical models because due to periodicity, infinite number of characteristic roots called characteristic exponents [10] are obtained. However, in many cases, the effect of sampling is often neglected due to the applied sufficiently high sampling frequency or treated as a continuous-time delayed system [11][12]. Although these simplifications may be acceptable from an engineering point of view, it becomes essential to handle systems in discrete-time.
For example, the applied discrete-time proportional control in position control tasks is always unstable [13,14] even at high sampling rates. The effect of spatial discretization, i.e., quantization, can cause chaotic behavior [15,16], or the sampled-data systems can have particular concave envelope vibrations when only dry friction stabilizes the motion [17,18]. It was also highlighted in [18] that the effect of sampling could result in multi-frequency vibration even in a single-degreeoffreedom system model. Still, an explanation for a deep understanding of the phenomenon was not presented.
In recent years, the effect of sampling on system dynamics was analyzed in many cases. For example, tools and design methods were introduced in [19] for digitally controlled robots, a detailed analytical investigation was carried out related to discrete-time force control [20][21][22], or stability and bifurcation analysis was presented in [23] discrete-time position control.
The paper is organized as follows. In Section 2, the mathematical model of a single-degree-of-freedom sampled-data system is introduced. In Sections 3 and 4, the stability and the dynamic analysis corresponding sampled-data system are carried out. The frequency analysis in Section 5, and Section 6 concludes the paper.
2 Model of sampled-data system A single-degree-of-freedom mechanical system is considered where discrete-time state feedback with zero-orderhold signal recognition is used to settle in zero reference position resulted in a sampled-data system.
The corresponding equation of motion is: where q(t) is the generalized coordinate as a function of time t, and dot refers to differentiation concerning time.
The generalized mass is denoted by m and u corresponds to the piecewise constant generalized control force as the control input. In addition, t j = jτ with j   means the j th sampling instant, where τ is the sampling time. It is noted that due to the applied zero-order-hold, the generalized control force u(t j ) is constant between two consecutive sampling instants, i.e., for the time interval t  (t j , t j + τ). By the applied state feedback, the control force is determined by the linear combination of the discrete states as u(t j ) = -k p x j -k d v j , where k p and k d are the feedback gains, and x j ≔ q(t j ) and v j ≔ q̇(t j ) defines the sampled position and velocity at the beginning of the j th time interval, i.e., the discrete states.
Equation (1) forms a non-homogeneous differential equation with piecewise constant right-hand side due to the applied zero-order-hold. It follows that, between two consecutive sampling instants, Eq. (1) can be solved analytically for the discrete state variables by integrating over the given time interval.
Introducing the dimensionless time T = t/τ, the number of free parameters can be reduced, and the equation of motion can be rewritten as where T j = j means the dimensionless sampling instant, and prime denotes the differentiation with respect to the dimensionless time. With the rewritten control law u(T j ) = -k p x j -k d /τ v j the following discrete map can be derived x j+1 = Hx j , where the discrete state vector x j = [x j v j ] ┬ , and: with the definitions p k m d k m : and .

Stability analysis
Stépán et al. [19] show methods that handle the difficulties of the digital effects, and based on this method the stability analysis can be carried out. Later, this method was used in the study of various engineering tasks, e.g., in force control [20][21][22], position control [23][24][25][26], or balancing [27] tasks. Using the discrete map x j+1 = Hx j , the time series of discrete states can be determined for any initial states x 0 as x n = H n x 0 . This series forms a multidimensional geometric series called Neumann series, which convergence is equivalent to ρ(H) < 1 where ρ means the spectral radius [7]. It follows, the eigenvalues of matrix H called characteristic multipliers determine the asymptotic stability of the system, which is stable if, and only if, the eigenvalues are inside the complex unit disk [7], i.e., |z 1,2 | < 1.
Using the characteristic equation of matrix H: the asymptotic stability of the system can be determined. By applying the Mobius transformation z = (w + 1)/(w -1) with w  C, the stability can be determined by the RouthHurwitz stability test for w which results in the stability conditions as p > 0, d > p/2 and 2 > d. The corresponding stable domain is illustrated in Fig. 1 in the frame of the dimensionless control gains p and d as the shaded triangular area. The discriminant curve marked by the dashed line shows that within the stable domain, two distinct real roots or a complex conjugate pair are formed depending on the selected parameters.
Based on the characteristic equation (Eq. 4), the characteristic multipliers can be expressed as and the types of stability loss at the stability boundaries, can be determined according to Eq. (5) collected in Table 1 [10].
Although, in the case of sampled-data systems, there are cases where vibrations are still observed in the time domain even though the roots are real. Thus, further investigation is needed regarding the location of the characteristic multipliers.

Dynamic analysis
The location of the characteristic multipliers can be classified by the discriminant curve ∆ and also by the curves given by the scalar invariants of matrix H, i.e., T D : : det p/ shows the location of the real part of roots, and  = ⇔ = + 0 1 2 d p/ separates the two cases when the roots have a different sign. The discriminant curve is given by ∆ : These curves divide the stable domain into six regions, as is collected in Table 2.
Based on Eq. (5), the different locations of characteristic multipliers can be determined, and these are shown by blue dots in Fig. 2 where the red boxes represent /2.
The location of the characteristic exponents determines the different types of possible motion. Thus, applying the inverse transformation as s = ln(z)/τ in general, the characteristic multipliers are transformed back to s domain to obtain the characteristic exponents. It is noted that due to the introduced dimensionless time, the transformation rule simplifies to s = ln(z).
To apply the inverse transformation in an easier way, the characteristic multipliers are reformulated to the complex exponential form z = ρ exp(iϑ), where ρ is the magnitude and ϑ is the argument of z. By taking into consideration that ϑ is periodic, i.e., ϑ = ϑ 0 + 2kπ with k   where ϑ 0 is the principal value which corresponds to k = 0, the characteristic exponents' yield where ϑ 0 can be determined by the two-argument variant of arctangent function.

Sub-domain #1
In this sub-domain, as it is shown in Fig. 2, the characteristic multipliers z 1,2  R are z 1 > z 2 > 0. Applying the inverse transformation, the characteristic exponents are s k k i     . 2 Location of the characteristic multipliers in the complex plain z As an example, the typical motion of this sub-domain is presented as black in Fig. 4. As it can be seen, the system shows aperiodic motion, i.e., there are no vibrations.
The dominant motion can be characterized by k = 0. Thus the dominant characteristic exponents become s 1 = ln(ρ 1 ) and s 2 = ln(ρ 2 ). If these are taken into account, the generalized position q(T) can be given as where the constants c 1 and c 2 come from the initial conditions q(0) = x 0 and q'(0) = v 0 . The corresponding motion is presented as dashed red in Fig. 4.

Sub-domains #2 and #3
In this sub-domain, as it is shown in Fig. 2, the characteristic multipliers z 1,2  R are z 1 > 0 and z 2 < 0. Applying the inverse transformation, the characteristic exponents become s k k i The difference between these sub-domains can be found in the sign of /2. In sub-domain #2, /2 > 0 Þ 0 > ρ 1 > ρ 2 , while in sub-domain #3, /2 < 0 Þ 0 > ρ 2 > ρ 1 . In terms of motion, it means that in sub-domain #2, the first-order dynamics, while in sub-domain #3, the vibration becomes dominant.
In both cases, the dominant motion cannot be characterized by k = 0 because it results in a single complex number without its conjugate pair. It can be concluded that if one of the characteristic multipliers is real and negative, additional harmonics are needed corresponding to k = -1. As an example, the typical motion of this sub-domain is presented as black in Figs. 7-8.
If the dominant characteristic exponents s 1 = ln(ρ 1 ) and s 2,3 = ln(ρ 2 ) ± iπ are taken into account, the generalized position q(T) can be given as where the constants c 1 , c 2 and c 3 come from the initial conditions. But, based on the original system, only two conditions q(0) = x 0 and q'(0) = v 0 are known. Therefore, a further condition is needed to get c 3 . For example, it can be given from an additional property of the original system, i.e., due to the applied zero-order-hold, the initial acceleration a 0 = -px 0 -dv 0 is constant.
Nevertheless, after some algebraic manipulation, Eq. (8) can be reformulated to where the constants c 1 , c 2 can be determined without any other condition. But, the specific solution results in a complex-valued function, and its real part means the position function. The corresponding motion is presented as red in Figs. 7-8.
Similar to the sub-domains #2 and #3, the dominant motion cannot be characterized without the effect of higher harmonics corresponding to k = -1. As an example, the typical motion of this sub-domain is presented as black in Fig. 10.
If the dominant characteristic exponents s 1,2 = ln(ρ 1 ) ± iπ and s 2,3 = ln(ρ 2 ) ± iπ are taken into account, the generalized position q(T) can be given as where the constants c 1 and c 2 come from the initial conditions. The corresponding motion is presented as dashed red in Fig. 10. In Fig. 10, a slight difference can be observed between the simulation result and the vibrations determined by dominant characteristic exponents. However, a qualitative agreement can be observed between the curves.

Sub-domain #5
In this sub-domain, as it is shown in Fig. 2, the characteristic multipliers are z i k k  location of the characteristic multipliers and the corresponding characteristic exponents are shown in Fig. 11.
Similar to the previous sub-domains, the dominant motion cannot be characterized without the effect of higher harmonics corresponding to k = -1. As an example, the typical motion of this sub-domain is presented as black in Fig. 12.
If the dominant characteristic exponents s 1,2 = ln(ρ) ± iϑ 0 and s 3,4 = ln(ρ) ± i(2π -ϑ 0 ) are taken into account, the generalized position q(T) can be given as where ω 0 = ϑ 0 and ω 1 = 2π -ϑ 0 and the constants c n with n = 1..4 come from the initial conditions. But, based on the original system, only two conditions x 0 and v 0 are known. The further condition can be given as the initial acceleration a 0 = −px 0 -dv 0 is constant, and the initial jerk is j 0 = 0. The corresponding motion is presented as red in Fig. 12.

Sub-domain #6
In this sub-domain, as it is shown in Fig. 2, the characteristic multipliers are z i k ( ) = ( )± + ( ) ρ ϑ π where 0 < ϑ 0 < π/2. The location of the characteristic multipliers and the corresponding characteristic exponents are shown in Fig. 13.
As an example, the typical motion of this sub-domain is presented as black in Fig. 14.
In this sub-domain, the dominant motion can be characterized by k = 0. If the dominant characteristic exponents s 1,2 = ln(ρ) ± iϑ 0 are taken into account, the generalized position q(T) can be given as where ω 0 = ϑ 0 and the constants c 1 and c 2 come from the initial conditions The corresponding motion is presented as red in Fig. 14.

Frequency analysis
According to the Nyquist-Shannon sampling theorem [4], a continuous signal containing finite frequency components can be fully recovered from a finite number of samples if the sampling frequency f S is at least twice the maximum frequency f max present in the signal, i.e., the following condition has to be held f max ≤ f S /2 = f N , where f N is the folding, or Nyquist frequency [4]. It is noted that due to the introduced dimensionless time, the sampling frequency is also dimensionless, and it becomes f S = 1. When this condition 2f max ≤ f S is not satisfied, aliasing [4] (frequency folding) may occur. However, this phenomenon does not influence the stability of the system. In sampled-data systems, applying the Nyquist-Shannon sampling theorem may be difficult to interpret because the response of these systems contains an infinite number of To avoid this problem, the dominant frequency component is usually used instead of the maximum frequency because the power of the upper harmonics is relatively tiny. For example, in [28], the natural frequency of the uncontrolled system, while in [29][30], the natural frequency related to the virtual stiffness parameter was used.
In the case of the Nyquist-Shannon sampling theorem for the dominant frequency f 0 , the following condition is obtained 2f 0 ≤ f S . It means that the characteristic multipliers can have the principal phase angle ϑ 0 ≤ π.
Simulation results are presented as black in Figs. 15-17 where the control parameters were selected from the stability boundary where Neimark-Sacker bifurcation [10] occurs. The left panels of these figures present the position signal in time, while the right panels show their spectra.
The selected control parameters and the obtained characteristic multipliers are collected in Table 3.
Based on these, the natural angular frequencies ω 0 , ω 1 = 2π -ω 0 and the related natural frequencies f 0 and f 1 = f S -f 0 were determined and collected in Table 3.
Figs. 15 to 17 correspond to the stability boundary of sub-domains #5 and #6 while Fig. 16 is related to the point where these sub-domains meet each other. In these figures, the position signal of the dominant motion is presented as dashed red.
It is noted, at the point of the boundary of sub-domains #5 and #6, the dominant motion can be determined with and without the upper harmonics corresponding to k = -1.
Based on the spectra presented in Figs. 15-17, it can be seen that when the further harmonics corresponding to k = -1 has to be taken into account, the difference of the frequency components called beating frequency f beat = f 1 -f 0 = f S -2f 0 is smaller (or equal) to the folding frequency, i.e., In can be concluded that the condition in Eq. (13) gives an upper limit to the sampling frequency.

Conclusions
In the present paper, the main characteristics of sampleddata system were investigated through an example of a position control task with discrete-time state-feedback. It was highlighted that even a single degrees-of-freedom system could have unexpected multi-frequency vibrations due to the periodicity of the characteristic multipliers.
In summary, in the case where the characteristic multipliers are positive real numbers or complex conjugate pairs with positive real-part, the higher harmonics become omitted. In all other cases, they must be taken into account.
In addition, the paper explores the possibility of creating a dominant model from a fundamentally infinite-dimensional system and formulates a condition for deciding whether it is necessary to consider higher harmonics.