Active damping for vibration control based on the switched stiffness technique
Rafael Ayala^{1} , Luis Moreno-Ahedo^{2}
^{1, 2}Faculty of Engineering, Universidad Autónoma de Baja California, Mexicali, México
^{1}Corresponding author
Journal of Vibroengineering, Vol. 23, Issue 4, 2021, p. 891-907.
https://doi.org/10.21595/jve.2021.21729
Received 6 October 2020; received in revised form 27 January 2021; accepted 25 February 2021; published 14 March 2021
An active damping controller based on the switched stiffness technique is developed and applied to vibration mitigation in a lightly damped structure. The controller either increases or decreases the system stiffness according to a state-dependent rule. A novel stability analysis based on the Floquet theory is proposed and employed to analyze the mathematical model of a one-degree-of-freedom system where the friction forces are taken into account. This novel analysis allows us to prove the exponential stability of the system origin, to establish a tuning procedure for the controller gain, to solve an optimization problem and to show the controller robustness against parameter uncertainties. Experimental verification is conducted to validate the effectiveness of the controller, and it is shown that the controller is feasible for vibration control problems.
- State-dependent control law produces a periodic response
- The period of the system response can be calculated analytically
- A lightly damped system is stabilized and optimal values can be obtained for the controller
Keywords: switchable stiffness, periodic control, Floquet theory.
1. Introduction
The stiffness modification is a well-known semi-active control technique for the vibration mitigation in structures and it was first introduced by Chen [1] and studied extensively by Onoda [2, 3]. This technique is achieved by semi-active devices whose stiffness is changed according to a state-dependent control rule. The most employed rule corresponds to switch the device between a maximum and minimum stiffness value, i.e., the semi-active device performs a switched stiffness control strategy [4].
Over the years, an enormous amount of research has been devoted to the implementation and experimental verification of the switched stiffness technique through developing semi-active stiffness devices. For instance, in [2] a piezoelectric actuator is implemented, while in [5] a spring with a mechanical arrangement is used to switch the stiffness. In [6] a valve is used as a semi-active device, while in [7] an electromagnetic device is developed, and recently a magnetostrictive transducer is developed in [8], to name a few.
Although a considerable body of research has been done in developing semi-active stiffness devices, few attempts have been made to investigate the relation between the device parameters and their performance for vibration control in structures. Pioneer work was done in [9] where it is performed a characterization of the parameters of the semi-active devices by introducing a dimensionless parameter that relates the maximum and minimum stiffness. In the same work, it is found the optimal parameter value for which the damped response is improved. However, in practical terms, such optimal value cannot be handled by the switched stiffness devices. Consequently, the limitations of the switched stiffness technique are increasingly apparent.
This study aims to extend the applicability of the switched stiffness technique by considering that the control action is applied actively via an actuator. For illustrative purposes, the active damping of a one-degree-of-freedom system is analyzed considering the free vibration case. That is, consider the physical realization of the one-degree-of-freedom system given by the mass-spring system shown in Fig. 1, then the switched stiffness technique is applied actively by the force $u$ defined by $u=-{k}_{c}\mathrm{s}\mathrm{g}\mathrm{n}\left(\stackrel{-}{x}{\stackrel{-}{x}}^{\text{'}}\right)\stackrel{-}{x}$, where ${k}_{c}$ is the controller gain. The control action $u$ increases or decreases the system stiffness at the specific instants.
Thereupon, our main goal is to determine the effectiveness of active control for vibration mitigation in structures considering a switched stiffness technique. Also, we conduct experimental verification of the approach.
Fig. 1. Actuated mass-spring system with friction
This paper is organized as follows. The Section 2 is dedicated to problem formulation and some fundamental ideas underlying the main problem are presented, while Section 3 details our theoretical findings based on the application of the Floquet theorem, to clarify our results, background information about the theory is briefed in the appendix. The Section 4 provides numerical and experimental results that support our findings and some final comments on the results are given in the conclusion Section 5.
2. Problem formulation
Consider the one-degree-of-freedom system given by the lightly damped mass-spring system shown in Fig. 1. The body slides on a surface with friction, and we assume that the dissipative phenomenon due to the friction force ${f}_{s}$ can be approximated by a linear viscous damping force $b\dot{x}$, where $b$ is the damping coefficient. In practice, this assumption provides mathematical tractability and simplicity, besides the damping coefficient can be estimated by well-known techniques [10, 11]. The characterization of the involved dissipative forces as a linear viscous damping force is a proven approximation and is feasible for a neighborhood around the origin. Accordingly, the mathematical model of the system is given by:
where $m$ is the mass, $k$ is the stiffness, $u$ is the control law, the prime denotes differentiation with respect to time $\tau $ and $\left(\stackrel{-}{x}\right(0),{\stackrel{-}{x}}^{\text{'}}(0\left)\right)=({\stackrel{-}{x}}_{0},{\stackrel{-}{x}}_{0}^{\text{'}})$ are non-zero initial conditions. The control law $u$ is defined by:
where ${k}_{c}$ is the controller gain and $sgn(\cdot )$ denotes the standard sign function.
Let us define the non-dimensional parameters: the damping ratio $\xi =b/\left[2\sqrt{km}\right]$, the natural frequency ${\omega}_{0}^{2}=k/m$ and ${\omega}_{c}^{2}={k}_{c}/m$. Thus, the time scale $t={\omega}_{0}\tau $ and the variable change $x=\stackrel{-}{x}/L$, where $L$ is a constant, transform Eq. (1) into:
where the dots denote differentiation with respect to time $t$, the initial conditions are reduced to $({x}_{0},{\dot{x}}_{0})=({\stackrel{-}{x}}_{0}/L,\mathrm{}{\stackrel{-}{x}}_{0}^{\text{'}}/\left[{\omega}_{0}L\right])$ and the tuning parameter $\epsilon $ is defined by $\epsilon ={\omega}_{c}^{2}/{\omega}_{0}^{2}={k}_{c}/k$, hence, the controller gain is given by:
When the viscous damping is not present, Eq. (3) is known as the Reid equation which has been used to describe the hysteretic damping occurred intrinsically in materials [12, 13]. Therefore, it is said that the switched stiffness technique induces a nonlinear dissipative phenomenon related to the hysteretic damping [14, 15]. For the analysis of this equation, several nonlinear analysis techniques have been broadly applied. For instance, in [14] the solution is approximated by the harmonic balance technique and its stability is determined by the variational equation, in [16] an extended analysis is done using perturbation techniques and the multi-degree-of-freedom case is studied in [17], while in [15, 18] a piecewise analytic approach is carried out, in [19] the phase plane method is developed and applied and [9] applies an approximation technique based on method of variation of parameters. A novel approach using linear techniques was also developed in [20].
From Eq. (3) written as state space model:
we deduce that the origin is an isolated equilibrium point, and the existence and uniqueness of the solution of Eq. (5) are guaranteed because the function $\mathbf{f}$ satisfies the Lipschitz condition, it is bounded and piecewise continuous at $\mathbf{x}$ and it admits a limited number of finite discontinuities [21].
In general, the semi-active control technique posses inherent stability as it is pointed out in [22]. Notwithstanding, for Eq. (3), the main question is about under what circumstances of the parameters its origin is stable. The stability of the origin can be established by the Lyapunov theorem as follows.
Consider the candidate Lyapunov function $V\left(\mathbf{x}\right)=\frac{1}{2}{x}^{2}+\frac{1}{2}{\dot{x}}^{2}$, then, its derivative along the trajectories of Eq. (5) yields:
which is negative semi-defined, then by the Lyapunov theorem the trivial solution of Eq. (3) is stable. However, applying the Young inequality^{}$xy\le \frac{{\u03f5}^{p}}{p}{x}^{p}+\frac{1}{{\u03f5}^{q}q}{y}^{q}$, $\forall \u03f5>0$, $\frac{1}{q}+\frac{1}{p}=1$ to Eq. (6) yields:
where for $\epsilon >0$ and $\xi \ge 0$, the function $\dot{V}<0$ is negative defined and the origin is asymptotically stable.
From the above result, it is shown that the Lyapunov analysis delimits the parameter range to ensure asymptotic stability. However, in general, there is not a procedure to choose an appropriated tuning parameter $\epsilon $ neither its optimum value, given a prescribed performance index. To illustrate this issue, consider the numerical solutions of Eq. (3) plotted in Fig. 2 for different values of the tuning parameter $\epsilon $, along with the corresponding switching functions. From Fig. 2 we deduce several facts: first, the trajectories converge to zero as it was proved by Lyapunov analysis; second, the settling time of the response in the second case is shorter than the first case. Accordingly, a tuning methodology must be devised.
In Fig. 2, one more remarkable fact takes place, the switching function $p(x,\dot{x})=\mathrm{s}\mathrm{g}\mathrm{n}\left(x\dot{x}\right)$ seems to behave periodically in time, at this point it is important to point out that we claim that the sign function of the product $x\dot{x}$ behaves periodically in time. This fact allows us to perform a stability analysis treating Eq. (3) as a linear periodic equation instead of a nonlinear differential equation. The principal advantage of doing this is that for linear systems several well-known mathematical tools can be employed instead of nonlinear approximation techniques. The main mathematical tool for periodic linear systems is the Floquet theory which for reference is briefed in the appendix to clarify our findings. In the following sections, the application of the Floquet theory shall allow us to prove the exponential stability of the origin, to establish a tuning procedure for the controller gain, to solve an optimization problem and to prove the controller robustness.
Fig. 2. Numerical solutions of Eq. (3) with initial conditions $({x}_{0},{\dot{x}}_{0})=\left(\mathrm{1.25,0.25}\right)$: a) when $\xi =0.15$ and $\epsilon =0.1$ and b) when $\xi =0.15$ and $\epsilon =0.45$
a)
b)
3. Main results
The damped Reid equation Eq. (3), illustrates the case when the dissipative forces, which are ubiquitous in real-world systems, are taken into account. In this work, they are characterized and quantified by a viscous damping force, this assumption is feasible for a neighborhood around the origin [23]. We recall that in the absence of damping, Eq. (3) is reduced to the well-known Reid equation which has been broadly studied employing nonlinear analysis techniques, such studies in some cases are cumbersome and approximated.
We apply the Floquet theory to analyze the damped Reid equation. This approach is novel and simple, and it provides analytic results without employing approximations.
To apply the Floquet theory, it is necessary to rewrite the nonlinear differential equation Eq. (3) as a linear periodic differential equation, that is, we must prove that the function $\mathrm{s}\mathrm{g}\mathrm{n}\left(x\dot{x}\right)$ has a fixed period $T$ which must be independent of the initial conditions.
Recall that Fig. 2 suggests that switching function $\mathrm{s}\mathrm{g}\mathrm{n}\left(x\dot{x}\right)$ behaves periodically in time. In order to corroborate this fact, the set of arbitrary initial conditions $({x}_{0},{\dot{x}}_{0})=({\stackrel{-}{x}}_{0}/L,{\stackrel{-}{x}}_{0}^{\text{'}}/\left[{\omega}_{0}L\right])$ of Eq. (3) are nondimensionalized and split into two sets in the following way: if $L={\stackrel{-}{x}}_{0}$ and ${\stackrel{-}{x}}_{0}^{\text{'}}=0$ then the nondimensionalized initial conditions $({x}_{0},{\dot{x}}_{0})=\left(\mathrm{1,0}\right)$ along with Eq. (3) allow us to define the initial position problem, while if $L={\stackrel{-}{x}}_{0}^{\text{'}}/{\omega}_{0}$ and ${\stackrel{-}{x}}_{0}=0$ then the nondimensionalized initial conditions $({x}_{0},{\dot{x}}_{0})=\left(\mathrm{0,1}\right)$ along with Eq. (3) defined the initial velocity problem. In this manner, we define a set of linearly independent initial conditions.
For the initial position and velocity problem above defined, their solutions ${x}_{p}$ and ${x}_{v}$, along with the corresponding switching functions are plotted in Fig. 3(a) and Fig. 3(c) respectively. The plots reveal two remarkable facts: in both cases, the switching function $\mathrm{s}\mathrm{g}\mathrm{n}\left(x\dot{x}\right)$ behaves periodically in time as we expected and its period $T$ is fixed, and the second is that solutions approach asymptotically to the trivial solution with the same and fixed period. In this case the system described by Eq. (3) is said to be an asymptotically isochronous system, for a rigorous definition and examples see [24]. To prove the periodicity of the function, we must provide an analytic expression for the period $T$ by calculating the switching instants of the solutions as follows.
Fig. 3. Numerical solutions Eq. (3) when $\xi =0.25$ and $\epsilon =0.15$: a) solution for the initial velocity problem and its phase portrait b), c) solution for the initial position problem and its phase portrait d)
a)
b)
c)
d)
3.1. Periodicity
The phase portraits of the corresponding solutions ${x}_{v}$ and ${x}_{p}$ are plotted in Fig. 3(b) and Fig. 3(d) respectively. To calculate the switching instants ${t}_{i}$, we examine both solutions at each quadrant of their respective phase portraits, where the function $\mathrm{s}\mathrm{g}\mathrm{n}\left(x\dot{x}\right)$ is well-defined, and hence, Eq. (3) is solvable and has analytic solutions.
3.1.1. Initial velocity problem
For the initial velocity problem the quadrants are named as in Fig. 3(b). At each quadrant we have the following:
I. At the first quadrant $\mathrm{s}\mathrm{g}\mathrm{n}\left(x\dot{x}\right)=1$, then Eq. (3) is reduced to:
under the initial conditions $({x}_{0},{\dot{x}}_{0})=(0,\mathrm{}1)$. Its characteristic polynomial is ${s}^{2}+2\xi s+1+\epsilon =0$ whose roots are $s=-\xi \pm \sqrt{-(1-{\xi}^{2}+\epsilon )}$ for a oscillatory response $1-{\xi}^{2}+\epsilon >0$, then, $s=-\xi \pm i{\omega}_{h}$, where ${\omega}_{h}=\sqrt{1-{\xi}^{2}+\epsilon}$. Thus, the solution is given by ${x}_{1}\left(t\right)={c}_{1}{e}^{(-\xi +i{\omega}_{h})t}+{c}_{2}{e}^{(-\xi -i{\omega}_{h})t}$, where ${c}_{1}$ and ${c}_{2}$ are found applying the initial conditions, that is, the solution is ${x}_{1}\left(t\right)=\frac{1}{{\omega}_{h}}{e}^{-\xi t}\mathrm{s}\mathrm{i}\mathrm{n}{\omega}_{h}t$ and its derivative ${\dot{x}}_{1}\left(t\right)={e}^{-\xi t}\left(\mathrm{c}\mathrm{o}\mathrm{s}{\omega}_{h}t-\frac{\xi}{{\omega}_{h}}\mathrm{s}\mathrm{i}\mathrm{n}{\omega}_{h}t\right)$. At the first switching instant $t={t}_{1}$ we have ${\dot{x}}_{1}\left({t}_{1}\right)=0\text{,}$ i.e, $\mathrm{c}\mathrm{o}\mathrm{s}{\omega}_{h}{t}_{1}-\frac{\xi}{{\omega}_{h}}\mathrm{s}\mathrm{i}\mathrm{n}{\omega}_{h}{t}_{1}=0$ then ${t}_{1}=\frac{1}{{\omega}_{h}}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\frac{{\omega}_{h}}{\xi}$ and ${x}_{1}\left({t}_{1}\right)$ yields ${x}_{1}^{0}={x}_{1}\left({t}_{1}\right)={e}^{-\xi {t}_{1}}\frac{1}{\sqrt{{\xi}^{2}+{\omega}_{h}^{2}}}$.
II. At second quadrant the relation $\mathrm{s}\mathrm{g}\mathrm{n}\left(x\dot{x}\right)=-1$ holds, then Eq. (3) yields:
where $({x}_{1}^{0},{\dot{x}}_{1}^{0}=0)$. The roots of the characteristic polynomial are complex $s=-\xi {\omega}_{0}\pm \sqrt{-(1-{\xi}^{2}-\epsilon )}$ if $(1-{\xi}^{2})-\epsilon >0$, thus $s=-\xi {\omega}_{0}\pm i{\omega}_{l}$ where ${\omega}_{l}=\sqrt{1-{\xi}^{2}-\epsilon}$. Then, the solution is ${x}_{2}\left(t\right)={c}_{1}{e}^{(-\xi +i{\omega}_{l})t}+{c}_{2}{e}^{(-\xi -i{\omega}_{l})t}$, where ${c}_{1}$ and ${c}_{2}$ are integration constants applying the initial conditions, the solution is ${x}_{2}\left(t\right)=\frac{{x}_{1}^{0}{e}^{-\xi t}}{{\omega}_{l}}({\omega}_{l}\mathrm{c}\mathrm{o}\mathrm{s}t{\omega}_{l}+\xi \mathrm{s}\mathrm{i}\mathrm{n}t{\omega}_{l})$ and its derivative yields ${\dot{x}}_{2}\left(t\right)=-\frac{{x}_{1}^{0}{e}^{-\xi t}({\xi}^{2}+{\omega}_{l}^{2})}{{\omega}_{l}}\mathrm{s}\mathrm{i}\mathrm{n}t{\omega}_{l}$. At the switching instant ${t}_{2}$ we have ${x}_{2}\left({t}_{2}\right)=0$, that is, ${\omega}_{l}\mathrm{c}\mathrm{o}\mathrm{s}{t}_{2}{\omega}_{l}+\xi \mathrm{s}\mathrm{i}\mathrm{n}{t}_{2}{\omega}_{l}=0$ then ${t}_{2}=\frac{1}{{\omega}_{l}}\left(\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\right(-\frac{{\omega}_{l}}{\xi})+\pi )\text{,}$ and ${\dot{x}}_{2}\left({t}_{2}\right)$ yields ${\dot{x}}_{2}^{0}={\dot{x}}_{2}\left({t}_{2}\right)=-\frac{{x}_{1}^{0}{e}^{-\xi {t}_{2}}}{{\omega}_{l}}({\xi}^{2}+{\omega}_{l}^{2})\mathrm{s}\mathrm{i}\mathrm{n}{t}_{2}{\omega}_{l}={e}^{-\xi ({t}_{1}+{t}_{2})}\frac{\sqrt{{\xi}^{2}+{\omega}_{l}^{2}}}{\sqrt{{\xi}^{2}+{\omega}_{h}^{2}}}$.
III. At the third quadrant $\mathrm{s}\mathrm{g}\mathrm{n}\left(x\dot{x}\right)=1$, then Eq. (3) is reduced to ${\ddot{x}}_{3}+2\xi {\dot{x}}_{3}+(1+\epsilon ){x}_{3}=0$ with the initial conditions $({x}_{2}^{0}=0,{\dot{x}}_{2}^{0})$. The solution ${x}_{3}$ and its derivative are given by ${x}_{3}\left(t\right)=\frac{{\dot{x}}_{2}^{0}}{{\omega}_{h}}{e}^{-\xi t}\mathrm{s}\mathrm{i}\mathrm{n}{\omega}_{h}t$ and ${\dot{x}}_{3}\left(t\right)={\dot{x}}_{2}^{0}{e}^{-\xi t}\left(\mathrm{c}\mathrm{o}\mathrm{s}{\omega}_{h}t-\frac{\xi}{{\omega}_{h}}\mathrm{s}\mathrm{i}\mathrm{n}{\omega}_{h}t\right)$ respectively. At ${t}_{3}$ we have ${\dot{x}}_{3}\left({t}_{3}\right)=0$, that is, ${t}_{3}=\frac{1}{{\omega}_{h}}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\frac{{\omega}_{h}}{\xi}$, and ${x}_{3}\left({t}_{3}\right)$ yields ${x}_{3}^{0}={e}^{-\xi ({t}_{1}+{t}_{2}+{t}_{3})}\frac{\sqrt{{\xi}^{2}+{\omega}_{l}^{2}}}{{\xi}^{2}+{\omega}_{h}^{2}}$.
IV. In the last quadrant $\mathrm{s}\mathrm{g}\mathrm{n}\left(x\dot{x}\right)=-1$, then ${\ddot{x}}_{4}+2\xi {\dot{x}}_{4}+(1-\epsilon ){x}_{4}=0$ under the initial conditions $({x}_{3}^{0},{\dot{x}}_{3}^{0}=0)$. The solution is ${x}_{4}\left(t\right)=\frac{{x}_{3}^{0}}{{\omega}_{l}}{e}^{-\xi t}({\omega}_{l}\mathrm{c}\mathrm{o}\mathrm{s}{\omega}_{l}t+\xi \mathrm{s}\mathrm{i}\mathrm{n}{\omega}_{l}t)$ and its derivative is ${\dot{x}}_{4}\left(t\right)=-{x}_{3}^{0}\frac{{e}^{-\xi t}}{{\omega}_{l}}({\xi}^{2}+{\omega}_{l}^{2})\mathrm{s}\mathrm{i}\mathrm{n}t{\omega}_{l}$. At ${t}_{4}$ the solution obeys ${x}_{4}\left({t}_{4}\right)=0\text{,}$ that is, ${t}_{4}=\frac{1}{{\omega}_{l}}\left(\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\right(-\frac{{\omega}_{l}}{\xi})+\pi )$, and ${\dot{x}}_{4}^{0}={\dot{x}}_{4}\left({t}_{4}\right)=-\frac{{x}_{3}^{0}{e}^{-{t}_{4}\xi}}{{\omega}_{l}}(\xi +{\omega}_{l}^{2})\mathrm{s}\mathrm{i}\mathrm{n}{t}_{4}{\omega}_{l}={e}^{-\xi ({t}_{1}+{t}_{2}+{t}_{3}+{t}_{4})}\frac{{\xi}^{2}+{\omega}_{l}^{2}}{{\xi}^{2}+{\omega}_{h}^{2}}$.
From the above analysis, the switching function $\mathrm{s}\mathrm{g}\mathrm{n}\left(x\dot{x}\right)$ can be written as the periodic function ${q}_{v}\left(t\right)$:
where the period $T$ is the sum of all the switching instants deduced above, i.e., $T={\sum}_{i=1}^{4}\mathrm{\u200d}{t}_{i}$, then:
where ${\omega}_{l}=\sqrt{1-{\xi}^{2}-\epsilon}$ and ${\omega}_{h}=\sqrt{1-{\xi}^{2}+\epsilon}$.
The value of the trajectory at the period $T$ is given by:
where the values of ${\omega}_{h}$, ${\omega}_{l}$ and $T$ were substituted.
Notice that conditions $1+\epsilon -{\xi}^{2}>0$ and $1-\epsilon -{\xi}^{2}>0$ are imposed to guarantee the periodicity of the switching function.
3.1.2. Initial position problem
Regarding the initial position problem, Fig. 3(c) shows that switching instants are different from those of the initial velocity problem, that is, the switching function $\mathrm{s}\mathrm{g}\mathrm{n}\left(x\dot{x}\right)$ for the initial position problem commutes at different instants. Following the above procedure can be shown that the period of the function remains as in Eq. (11). Thus, the switching function is written as the periodic function ${q}_{p}\left(t\right)={q}_{p}(t+T)$, and the value of the trajectory at $T$ is:
From the previous analysis, we deduce that for the initial velocity and initial position problems, the corresponding switching function can be written by either the $T$-periodic function ${q}_{v}\left(t\right)$ or ${q}_{p}\left(t\right)$ respectively, where the period is equal to Eq. (11). This result is formalized in the following theorem.
Theorem 1. The system Eq. (3) can be written as a linear periodic system:
where $q\left(t\right)=q(t+T)$ is a periodic function with period $T$ given by the relation Eq. (11).
To clarify the above result, the function $q\left(t\right)$ is provided by either the functions ${q}_{v}\left(t\right)$ or ${q}_{p}\left(t\right)$ respectively, hence, its definition depends on which set of standard initial conditions is chosen. However, the Floquet theorem requires only the knowledge of the period $T$, thus, the definition of $q\left(t\right)$ is not necessary to be provided explicitly. Therefore, the period and the value of the trajectories ${x}_{p}$-${x}_{v}$ at $T$ given by Eq. (12-13) allow us to determine the stability through the monodromy matrix.
3.2. Stability analysis
The stability of periodic system desribed by Eq. (14) is determined by the monodromy matrix $\mathbf{\Phi}\left(T\right)=\mathbf{W}\left(T\right){\mathbf{W}}^{-1}\left(0\right)$, where $\mathbf{W}$ is the Wronskian matrix. The matrix $\mathbf{W}$ is given by the linearly independent solutions ${x}_{v}$ and ${x}_{p}$. Therefore, $\mathbf{\Phi}\left(T\right)$ yields:
where column vectors are given by Eqs. (12-13). The multipliers $\mu $ are the eigenvalues of the monodromy matrix:
the multiplier corresponds to a double semi-simple eigenvalue. While the relation Eq. (24) provides the Floquet exponent $\rho $:
From the above results allow us to state the following stability criterion
Corollary 1. The origin of the system Eq. (3) is asymptotically stable if $\epsilon <1.$
Proof 1. The Floquet multiplier Eq. (16) satisfies $\mu <1$ when $\epsilon <1$. Therefore, by the Floquet theory the origin is asymptotically stable.
The representation of Floquet exponent Eq. (17) is unique since $\mu <1$ and by corollary 4 the Floquet factors are real.
The Floquet theorem guarantees asymptotic stability. However, exponential stability can actually be proved. For that aim, Eq. (3) is reduced to a linear system with constant coefficients
Corollary 2. The system Eq. (3) is reducible to the linear system with constant coefficients:
where the solution is given by $y\left(t\right)={e}^{\rho t}{y}_{0}$, hence, the origin of system Eq. (3) is exponentially stable.
Proof 2. The Lyapunov transformation $\mathbf{x}=\mathbf{F}\left(t\right)\mathbf{y}$ transforms the system Eq. (3) into $\dot{\mathbf{y}}=\mathbf{K}\mathbf{y}$, where the matrix $\mathbf{K}=\left(\frac{1}{T}\right)\mathrm{l}\mathrm{n}\mathbf{\Phi}\left(T\right)$ is real and the monodromy matrix $\mathbf{\Phi}\left(T\right)$ is given by Eq. (15). Applying the lemma 1 (see Appedix) the trivial solution of Eq. (3) is exponentially stable.
Notice that the present analysis provides analytic closed expressions in terms of the system parameters for determining the stability, also, to establish a controller design procedure and to solve an optimization problem based on a performance index, as follows.
3.3. Controller design
By the Floquet theorem the exponent $\rho $ given in Eq. (24) provides a stability criterion, that is, the exponent meets $\left|\rho \right|<1$ when the tuning parameter and the damping ratio satisfy $0<\xi <1$ and $0<\epsilon <1$ respectively. The damping ratio is bounded because we consider a lightly damped structure and $\epsilon $ by its definition and the stability analysis. Consider the parameter space $(\epsilon ,\xi )$ restricted to these bounds, then Floquet exponent in terms of the parameters is numerically calculated and plotted as the contour plot given in Fig. 4. Notice that for some values $(\epsilon ,\xi )$ the period $T$ given by Eq. (11) is complex, such complex values do not have a physical meaning.
Fig. 4. Contour plot of the Floquet exponent $\rho (\epsilon ,\xi )$
The plot provides the stability chart in the parameter space for which $\left|\rho \right|\le 1$ avoiding complex values. The white region represents the values $(\epsilon ,\xi )$ for which the period $T$ is complex, while the region in color represents those real values for which the trivial solution of the system is exponentially stable.
The gain ${k}_{c}=\epsilon k$ of the controller is determined using the stability chart by considering the estimated value of the equivalent damping ratio $\xi $ and choosing an appropriated value for the tuning parameter $\epsilon $. For instance, consider the estimated damping ratio $\xi \approx 0.0932$ and the value $\epsilon =0.4$. In Fig. 4 the corresponding point $(\epsilon ,\xi )$ label as ${P}_{1}$.
A greater distance from the poles to the imaginary axis allows a faster convergence to zero of the system response. For this reason, it is important to find parameter setting that maximizes this distance. Regarding the optimal $\epsilon $, it is possible to state an optimization problem considering the root loci of the poles $\rho (\epsilon ,\xi )$ of the reduced linear system described in Eq. (18), that is, an appropriate stability margin can be provided by maximizing the distance between the poles and the imaginary axis.
The maximum distance is accomplished when the derivative of the pole is equal to zero, that is when ${\left.d\rho /d\epsilon \right|}_{\xi}=0$. Notice that such derivative exits and corresponds to a Frechet derivative because the Floquet exponent is semi-simple [25].
Therefore, for a damping ratio $\xi \in \left[\mathrm{0,1}\right]$, the derivative ${\left.d\rho /d\epsilon \right|}_{\xi}=0$ yields a transcendental equation which is solved for the optimal value ${\epsilon}_{op}$ applying the Newton-Rapson method. The result of this iterative procedure is shown in the contour plot by the set of points denoted by $d\rho /d\epsilon =0$. Therefore, the set of points $({\epsilon}_{op},\xi )$ are the values where the distance between the poles and the imaginary axis is maximized.
On the other hand, the tuning procedure employs the stability chart to choose the tuning parameter $\epsilon $ for a given damping ratio $\xi $. As expected, this design procedure entails some source of uncertainty, mainly in the identification process of the damping ratio. Therefore, it is natural to ask what about the uncertainty of the parameters; hence, we shall perform a sensitivity analysis to study the effects of parameter variations on the Floquet exponent.
3.4. Sensitivity analysis
Let consider that ${\xi}_{0}$ and ${\epsilon}_{0}$ are the nominal values, then the nominal Floquet exponent is given by ${\rho}_{0}=\rho ({\xi}_{0},{\epsilon}_{0})$, we wish to measure the change of ${\rho}_{0}$ as the parameters $\xi $ and $\epsilon $ change from their nominal values. For that aim, let $\mathrm{\Delta}\epsilon =\epsilon \frac{\left[c\right]}{100}\mathrm{\%}$ denotes a change in $\epsilon $, where $\left[c\right]$ represents percent changes, for example, $\left[c\right]=$ 1, 5, 10 %, etc. For a change in $\epsilon $ the Floquet exponent yields ${\rho}_{\mathrm{\Delta}\epsilon}=\rho ({\xi}_{0},{\epsilon}_{0}+\mathrm{\Delta}\epsilon )$, then, $\mathrm{\Delta}\rho ={\rho}_{0}-{\rho}_{\mathrm{\Delta}\epsilon}$ measure the variation with respect to $\epsilon $ and let $\left[\frac{\mathrm{\Delta}\rho}{{\rho}_{0}}\right]=100\frac{\mathrm{\Delta}\rho}{{\rho}_{0}}$ denotes the percent change.
Consider $d\rho =\frac{d\rho}{d\epsilon}d\epsilon ={S}_{\epsilon}^{\rho}d\epsilon $, where ${S}_{\epsilon}^{\rho}=\frac{d\rho}{d\epsilon}$ is known as the sensitivity function of $\rho $ with respect to $\epsilon $, see [26]. Then, $\mathrm{\Delta}\rho \approx {S}_{\epsilon}^{\rho}\mathrm{\Delta}\epsilon $ and the relation $\left[\frac{\mathrm{\Delta}\rho}{{\rho}_{0}}\right]$ yields:
Eq. (19) measures the change in $\rho $ as $\epsilon $ varies. However, the parameter $\xi $ may also vary then:
where ${S}_{\xi}^{\rho}=\frac{d\rho}{d\xi}$. The above expression measures the change in $\rho $ as $\xi $ varies. In addition, $\rho $ can also vary when both parameters $\rho $ and $\xi $ have variations then:
Notice that the sensitivity functions ${S}_{\epsilon}^{\rho}$ and ${S}_{\xi}^{\rho}$ are well-defined since the Floquet exponent is semi-simple, and hence their respective derivative corresponds to a Frechet derivative [25].
When the tuning parameter corresponds to the optimal value ${\epsilon}_{op}$ the sensitivity function ${S}_{\epsilon}^{\rho}$ is approximately equal to zero since ${\left.d\rho /d\epsilon \right|}_{{\epsilon}_{op}}\approx 0$. Then, we can deduce that the controller is robust against small changes in the parameters.
4. Numerical and experimental results
In this section, we carry out the simulation and experimental verification of the switched stiffness active controller. The numerical simulation of the mass-spring system is based on the system parameters of the experimental model depicted in Fig. 5.
Fig. 5. Mass-spring car system: 1 – actuator, 2 – optical encoder, 3 – mass-spring car, 4 – rack-and-pinion gear-set, 5 – scale of ±3.0 cm
4.1. Numerical results
The dynamics of the experimental mass-spring system can be described by Eq. (1), whose parameters correspond to: $m=1.2\mathrm{}$kg, $b=4.5\mathrm{}$N/m and $k=487.8$ Ns/m respectively, where the estimated damping ratio is equal to $\xi =0.0932$. The methodology to estimate the damping ratio is shown in a later section. Based on these parameter we carry out the simulation of controller.
4.1.1. Simulation
The controller gain is ${k}_{c}=k\epsilon =487.8\epsilon $ and the tuning parameter $\epsilon $ is chosen arbitrarily as example as in the stability chart of Fig. 4. As illustrative example, the point ${P}_{1}(\xi ,\epsilon )=\left(\mathrm{0.0932,0.4}\right)$ in the stability chart is selected, thus, the controller gain is ${k}_{c}\approx 195$.
The numerical solution ${\widehat{x}}_{p}$ of the controlled system Eq. (1) for the initial condition ${\widehat{x}}_{p}\left(0\right)=2.0$ cm is plotted in Fig. 6, along with the control law Eq. (2) and the switching function.
The mass-spring system under the action of the controller becomes a linear periodic system, whose multiplier $\mu $ and Floquet exponent $\rho $ are calculated by Eq. (16-17), that is, $\mu =0.224682$ and $\rho =-0.215485$ respectively. To validate our predictions, the Floquet exponent may roughly be estimated from the numerical solution using the logarithmic decrement technique. From Fig. 6(a) the logarithmic decrement $\stackrel{~}{\rho}\approx \frac{1}{2\pi}\mathrm{l}\mathrm{o}\mathrm{g}\frac{2}{0.4474}=0.2383$ approximates the Floquet exponent. The error in the approximation is because the actual system response can not the approximated by the oscillatory response of the second-order linear system as the logarithmic decrement technique intends to do.
Notice the amount of stiffness introduced into the system is given by relation ${k}_{c}\mathrm{s}\mathrm{g}\mathrm{n}\left(\stackrel{-}{x}{\stackrel{-}{x}}^{\text{'}}\right)$. On the other hand, the design procedure of the controller entails some sources of uncertainty, due mainly to the damping ratio identification. To overcome this disadvantage and to show the robustness of the controller against parameter uncertainty, a sensitivity analysis is carried out in the following.
4.1.2. Robustness analysis
The nominal parameters are $({\xi}_{0},{\epsilon}_{0})=\left(\mathrm{0.0932,0.4}\right)$ and the corresponding Floquet exponent is ${\rho}_{0}=0.215485$. We consider that both parameters vary, then Eq. (19) is used to tabulate, in Table 1, the sensitivity analysis results. The table relates the percent variations of the Floquet exponent ${\rho}_{0}$ due to variations of the parameters $({\xi}_{0},{\epsilon}_{0})$.
Fig. 6. Numerical solution of the control system Eq. (1) where $m=1.2$ kg, $b=4.5$ N/m, $k=487.8$ Ns/m and ${k}_{c}=195$: a) position ${\widehat{x}}_{p}$, b) corresponding control law $u$, c) switching function
The data in Table 1 is interpreted as follows. For instance, consider a variation of 25 % of the parameter $\xi $ with no variations in $\epsilon $, then from the table 1 the percent variation is ${\left[\mathrm{\Delta}\rho /{\rho}_{0}\right]}_{\epsilon}=$ 14 %. As second example, now consider a variation in both nominal parameters, for instance, $\xi $ with 20 % and $\epsilon $ with 10 %, then Floquet exponent has a variation of ${\left[\mathrm{\Delta}\rho /{\rho}_{0}\right]}_{\epsilon ,\xi}=15\mathrm{}\mathrm{\%}$. In both examples, the corresponding results are pointed out in table 1 by the color text.
Table 1. Percent variation $\left[\mathrm{\Delta}\rho /{\rho}_{0}\right]$ of the Floquet exponent ${\rho}_{0}=0.215485$ as parameters vary from their nominal values $({\xi}_{0},{\epsilon}_{0})=\left(\mathrm{0.0932,0.4}\right)$
$\u2206\xi $, % | $\u2206\epsilon $, % | ||||
0 | 5 | 10 | 25 | 30 | |
0 | 0 | 3 | 6 | 14 | 17 |
10 | 5 | 8 | 11 | 19 | 22 |
20 | 9 | 12 | 15 | 23 | 26 |
30 | 13 | 16 | 19 | 27 | 30 |
50 | 21 | 24 | 27 | 35 | 38 |
To visualize the effect of parameters variation on the system response, two examples are plotted in Fig. 6(a). The corresponding solutions of Eq. (1) that shows $\mathrm{\Delta}\xi =25\mathrm{}\mathrm{\%}$ with no variation in $\epsilon $ (bolded) and a variation in both parameters as $\mathrm{\Delta}\epsilon =10\mathrm{}\mathrm{\%}$, $\mathrm{\Delta}\xi =20\mathrm{}\mathrm{\%}$ (blue).
In addition, the Euclidean norm of the solution may be used as a quantitative measure for validate the robustness of the controller, i.e., for the nominal solution ${\Vert {\widehat{x}}_{p}\Vert}_{2}=0.1823$, and for the present examples ${\Vert [{\widehat{x}}_{p}{]}_{\mathrm{\Delta}\xi =25\mathrm{\%}}\Vert}_{2}=0.1725$ and ${\Vert [{\widehat{x}}_{p}{]}_{\mathrm{\Delta}\epsilon =10\mathrm{\%},\mathrm{\Delta}\xi =20\mathrm{\%}}\Vert}_{2}=0.1721$ respectively. From the above sensitivity analysis, one deduces that the controller is robust against parameter uncertainty.
4.2. Experimental model
The experimental model is comprised of the mass-spring car which is constrained to move horizontally on a rail, an optical encoder which measures mass motion, a DC servomotor with a rack-and-pinion gear-set adaptation which provides the control action; notice that the distance traveled by the mass is ±3 cm. Besides, Fig. 7 shows the fully assembled experimental setup. The dSPACE data acquisition platform is used to register the experimental data and to implement the control law in real-time by using the block diagram formulation through the Matlab/Simulink software, where the sample time of the data acquisition is ${T}_{s}=1$ ms.
Fig. 7. Experimental rig: 1 – power supply, 2 – dSPACE computer, 3 – PC with Matlab/Simulink, 4 – ECP rectilinear plant model 201, 5 – dSPACE panel connector
The mass-spring car experiments a friction force during its motion. For practical purposes, we assume that such force can be quantified by a linear viscous damping force, where the logarithmic decrement technique is applied to estimate the damping coefficient.
To apply the technique, the free vibration response ${\stackrel{~}{x}}_{p}$ is experimentally determined considering the initial position ${\stackrel{~}{x}}_{p}\left(0\right)=2.0$ cm. and zero initial velocity. The period of the free vibration response is found as $T\approx 0.3130$ seg. and a second peak amplitude $A\approx 1.05$. From the logarithmic decrement we have $\frac{1}{2\pi}\mathrm{l}\mathrm{o}\mathrm{g}\frac{2}{1.05}=\xi \sqrt{1-{\xi}^{2}}$, then, solving for the damping ratio yields $\xi \approx 0.0932$.
The stiffness and damping parameters are calculated by the relations ${\omega}_{0}=2\pi /\left[T\sqrt{1-{\xi}^{2}}\right]$, $k={\omega}_{0}^{2}m$ and $b=2\xi \sqrt{mk}$, where the mass is equal to $m=1.2$ Kg. Thus, the parameters are approximately equal to $k=487.7952$ N/m and $b=4.51$ Ns/m respectively.
4.3. Experimental verification
A crucial step in the controller implementation is to calculate the switching function $z=\mathrm{s}\mathrm{g}\mathrm{n}\left(x\dot{x}\right)$ since the velocity $\dot{x}$ must be obtained. Therefore, the controller performance is based on which technique is used to calculate the velocity. One common technique is the numerical differentiation of the position, but this approach may yield a velocity contaminated by noise. A second approach is the use of state observers, but this approach depends on an appropriate choose of the observer and the tuning of its gains. Nevertheless, in this paper, we developed a novel and heuristic approach to calculate the switching function without velocity measurement, where we take into account the dynamics of the mass during its motion.
4.3.1. Switching function calculation
To explain how we calculate the switching function without velocity measurement, we compare our approach to the position differentiation approximation scheme. The state observer formulation is not considered in the following comparison.
Fig. 8 shows the Matlab/Simulink block diagram for the experimental implementation of both approaches. The orange block calculates the switching function by differentiating the position while the cyan block executes the algorithm to calculate the switching function without the velocity measurement. The yellow block explained in the appendix corresponds to the input-output data acquisition relations of the experiment setup. The block diagram is converted to computer code and loaded into the dSPACE real-time platform which runs at a sample time of ${T}_{s}=1.0$ ms.
Fig. 8. Block diagram for calculating the switching function $\mathrm{s}\mathrm{g}\mathrm{n}\left(x\dot{x}\right)$ without measure the velocity
The experiment results of both approaches are shown in Fig. 9. Where Fig. 9(a) shows the position ${\stackrel{~}{x}}_{p}$ of the mass and its velocity, which has been calculated by numerical differentiation, for illustrative purposes, the velocity is scaled by the natural frequency, that is, ${\dot{\stackrel{~}{x}}}_{p}{\omega}_{0}$. The corresponding switching function based on numerical differentiation is plotted in Fig. 9(b).
Even when, in this particular case, the numerical differentiation scheme provides quite good agreement with our approach, it is well-known that it may be susceptible to be contaminated by noise and; and, in general, it is avoided.
We base our approach on the motion of the mass. The approach determines when the position ${\stackrel{~}{x}}_{p}$ crosses the zero and when it is decreasing, these operations yield the two logic signals plotted in Fig. 9(c) and Fig. 9(d) respectively. The exclusive OR logic operation between these signals yields the signal plotted in Fig. 9(e), and finally the signal conditioning yields the switching function plotted in Fig. 9(f) this signal is the same of Fig. 9(b). Therefore, our approach does not require velocity measurement, in this way the control law is implemented.
Fig. 9. Time history of the signals used for calculating the switching function
4.3.2. Experimental results
Fig. 10 shows the block diagram to implement the controller, the blue block calculate the switching function. The experimental results of the control law 2 are shown in Fig. 11, for control gain ${k}_{c}=195$ used in the numerical simulation. For comparative purposes, in Fig. 11 we plot the experimental response ${\stackrel{~}{x}}_{p}$ and the numerical simulation ${\widehat{x}}_{p}$.
Fig. 10. Block diagram of the controller implementation for the mass-spring system
Notice that there is a discrepancy between the experimental result and the solution. This fact is due mainly to the action of additional dynamics not taken into account in the model like the dry friction, the backlash of the rack-and-pinion gear-set, and the induced dynamic by the DC servomotor. Nevertheless, the controller is robust and can handle these drawbacks.
Fig. 11. Experimental result of the control system with ${k}_{c}=195$
5. Conclusions
The switched stiffness technique is a semi-active vibration control technique achieved by semi-active stiffness devices to change the total system stiffness to mitigate the vibration in structures. The system response depends on the parameters of the involved semi-active device, however, in real-world situations, the semi-active device cannot handle the optimal values that improve the system response. Consequently, the limitations of the switched stiffness technique are increasingly apparent.
To tackle this issue, the authors developed a feasible controller to vibration mitigation of a single-degree-of-freedom system where the dissipative phenomenon is considered. We show that the vibration control system behaves periodically in time, this remarkable fact allows us to analyze the system as a periodic differential equation, where the Floquet theory is applied to establish closed analytic expression for stability conditions in terms of parameters system instead of linear approximations. This allows us to prove the exponential stability of the system, to provide a tuning scheme for the controller gain, to perform an optimization process and to show the robustness of the controller by carrying out a sensitivity analysis.
Besides, the conducted experimental verification shows good agreement with the theoretical analysis and supports and validates our approach. The controller can be extended to be used for shock isolation.
References
- Chen J. C. Response of large space structures with stiffness control. Journal of Spacecraft and Rockets, Vol. 21, Issue 5, 1984, p. 463-467. [Publisher]
- Onoda J., Endot T., Tamaoki H., Watanabe N. Vibration suppression by variable-stiffness members. AIAA Journal, Vol. 29, Issue 6, 1991, p. 977-983. [Publisher]
- Onoda J., Sanot T., Kamiyama K. Active, passive, and semiactive vibration suppression by stiffness variation. AIAA Journal, Vol. 30, Issue 12, 1992, p. 2922-2929. [Publisher]
- Goh C. J., Caughey T. K. A quasi-linear vibration suppression technique for large space structures via stiffness modification. International Journal of Control, Vol. 41, Issue 3, 1985, p. 803-812. [Publisher]
- Ramaratnam A., Jalili N. A switched stiffness approach for structural vibration control: Theory and real-time implementation. Journal of Sound and Vibration, Vol. 291, Issues 1-2, 2006, p. 258-274. [Publisher]
- Leavitt J., Jabbar F., Bobrow J. E. Optimal performance of variable stiffness devices for structural control. Journal of Dynamic Systems, Measurement and Control, Transactions of the ASME, Vol. 129, Issue 2, 2007, p. 171-177. [Publisher]
- Ledezma Ramirez D., Ferguson N., Brennan M. An experimental switchable stiffness device for shock isolation. Journal of Sound and Vibration, Vol. 331, Issue 23, 2012, p. 4987-5001. [Publisher]
- Scheidler J. J., Asnani V. M., Dapino M. J. Dynamically tuned magnetostrictive spring with electrically controlled stiffness. Smart Materials and Structures, Vol. 25, Issue 3, 2016, p. 035007. [Publisher]
- Winthrop M. F., Baker W. P., Cobb R. G. A variable stiffness device selection and design tool for lightly damped structures. Journal of Sound and Vibration, Vol. 287, Issues 4-5, 2005, p. 667-682. [Publisher]
- Adhikari S., Woodhouse J. Identification of damping. Part I: Viscous damping. Journal of Sound and Vibration, Vol. 251, Issue 3, 2002, p. 477-490. [Publisher]
- Liang J. W., Feeny B. F. Identifying coulomb and viscous friction from free-vibration decrements. Nonlinear Dynamics, Vol. 16, Issue 4, 1998, p. 337-347. [Publisher]
- Reid T. J. Free vibration and hysteretic damping. The Journal of the Royal Aeronautical Society, Vol. 60, 1956, p. 544-283. [Publisher]
- Bert C. Material damping. Journal of Sound and Vibration, Vol. 29, Issue 2, 1973, p. 129-153. [Publisher]
- Caughey T. K., Vijayaraghavan A. Free and forced oscillations of a dynamic system with “linear hysteretic damping” (non-linear theory). International Journal of Non-Linear Mechanics, Vol. 5, Issue 3, 1970, p. 533-555. [Publisher]
- Zhang Y., Iwan W. D. Some observations on two piecewise-linear dynamic systems with induced hysteretic damping. International Journal of Non-Linear Mechanics, Vol. 38, Issue 5, 2003, p. 753-765. [Publisher]
- Caughey T., Vijayaraghavan A. Stability analysis of the periodic solution of a piecewise-linear non-linear dynamic system. International Journal of Non-Linear Mechanics, Vol. 11, Issue 2, 1976, p. 127-134. [Publisher]
- Caughey T., Vijayaraghavan A. Forced oscillations of a piecewise-linear non-linear dynamic system with several degrees of freedom. International Journal of Non-Linear Mechanics, Vol. 12, Issue 6, 1977, p. 339-353. [Publisher]
- Inaudi J. A. Modulated homogeneous friction: a semi–active damping strategy. Earthquake Engineering and Structural Dynamics, Vol. 26, Issue 3, 1997, p. 361-376. [Publisher]
- Liu C. S. Reid’s passive and semi-active hysteretic oscillators with friction force dependence on displacement. International Journal of Non-Linear Mechanics, Vol. 41, Issues 6-7, 2006, p. 775-786. [Publisher]
- Moreno Ahedo L., Diarte Acosta S. Stability analysis of linear systems with switchable stiffness using the Floquet theory. Journal of Vibration and Control, Vol. 25, Issue 5, 2019, p. 963-976. [Publisher]
- Khalil H. K. Nonlinear Systems. Prentice Hall, New York, 2nd edition, 1996. [Search CrossRef]
- Garrido H., Curadelli O., Ambrosini D. On the assumed inherent stability of semi-active control systems. Engineering Structures, Vol. 159, 2018, p. 286-298. [Publisher]
- Liang J. W. Damping estimation via energy-dissipation method. Journal of Sound and Vibration, Vol. 307, Issue 1, 2007, p. 349-364. [Publisher]
- Calogero F. Isochronous Systems. Oxford University Press, Oxford, 2008. [Publisher]
- Sun J. G. Sensitivity analysis of multiple eigenvalues (I). Journal of Computational Mathematics, Vol. 1, Issue 1, 1988, p. 28-38. [Search CrossRef]
- Frank P. M. Introduction to System Sensitivity Theory. Academic Press, New York, 1978. [Search CrossRef]
- Yakubovich V. A., Starzhinskii V. M. Linear Differential Equations with Periodic Coefficients. Wiley, New York, 1975. [Search CrossRef]
- Montagnier P., Paige C. C., Spiteri R. J. Real Floquet factors of linear time-periodic systems. Systems and Control Letters, Vol. 50, Issue 4, 2003, p. 251-262. [Publisher]