Fractional-Order LQR and State Observer for a Fractional-Order Vibratory System

Abstract

The present study uses linear quadratic regulator (LQR) theory to control a vibratory system modeled by a fractional-order differential equation. First, as an example of such a vibratory system, a viscoelastically damped structure is selected. Second, a fractional-order LQR is designed for a system in which fractional-order differential terms are contained in the equation of motion. An iteration-based method for solving the algebraic Riccati equation is proposed in order to obtain the feedback gains for the fractional-order LQR. Third, a fractional-order state observer is constructed in order to estimate the states originating from the fractional-order derivative term. Fourth, numerical simulations are presented using a numerical calculation method corresponding to a fractional-order state equation. Finally, the numerical simulation results demonstrate that the fractional-order LQR control can suppress vibrations occurring in the vibratory system with viscoelastic damping.

1. Introduction

Fractional calculus is a form of calculus in which the order of differentiation and integration is expanded or generalized from an integer number to a non-integer number. However, fractional calculus is not novel or strange and has a long history that traces back to Leibniz, who was one of the founders of calculus [1]. During these past 30 years, fractional calculus has been applied in engineering fields, especially in control engineering. In control engineering, dynamical systems, which cannot be described sufficiently by conventional integer-order calculus, have been modeled with fractional calculus, and controllers that use fractional calculus have been designed. In other words, a revolution in the construction methods of control system has occurred.

Fractional calculus has been introduced into proportional-integral-differential (PID) control systems as a typical classical control method. Several examples can be cited: PI^λD^μ control [2], robust fractional-order PI control using the auto-tuning method [3], the application of this fractional-order PI control method to a manipulator robot [4], fractional-order PD control of a smart beam [5,6], a tuning rule for a fractional-order PD controller for a second-order system [7], fractional-order PID control for a magnetically levitated system [8], fractional-order PD control to suppress vibrations that appear in a three-story structure [9], a graphical design method for tuning a fractional-order PI or PD controller [10], a fractional-order PID controller to suppress vibration that occurs in a cantilever beam using electromagnetic actuators [11], the meaning of the parameters used in fractional-order PID controllers [12], fractional-order PID control for a pumped storage unit [13], fractional-order PID load-frequency control for a power system [14], tuning rules for fractional-order PID controllers [15], PI^λD^μ control for trajectory tracking by a 3-DOF parallel manipulator [16], fractional-order PID control for an automatic voltage regulating system [17], and fractional-order PID control for a hydraulic turbine regulating system [18]. Thus, research examples are too numerous to mention.

Fractional calculus has also been used for linear quadratic regulator (LQR) control as a typical modern control theory. In the sense that the controller is applicable to controlled systems with fractional-order dynamics, we use the terminology “fractional-order LQR” in this paper. Similarly, some examples can be cited: a state-feedback LQR controller for fractional-order systems [19], an LQR-based PID controller design using a fractional-order integral performance index [20], a fractional-order LQR [21], a fractional-order LQR control to suppress vibrations appearing in a three-story building [22], and fractional-order PID control for time-delay systems [23]. Again, research examples are too numerous to count. However, it is known that the procedure to derive a fractional-order LQR, as described in a previous report [21], is too complicated to be practical.

Therefore, in the present paper, a novel iteration-based method is first proposed to solve the algebraic Riccati equation to derive the fractional-order LQR. Since all states of the controlled system are fed back in the LQR control scheme, all of the states are assumed to be detected. However, in reality, it is very difficult to measure the states originating from the fractional-order derivative terms. A state observer is then established in order to estimate all of the states, including the fractional-order derivative states. The fractional-order LQR control system is then constructed using the estimated states. Next, by way of example, a viscoelastic damper system modeled with a fractional-order differential equation is considered, and fractional-order LQR control is conducted for the viscoelastic damper system. The effectiveness of the proposed design method is then confirmed. Finally, the proposed fractional-order LQR and the fractional-order LQR derived using the method proposed by Sierociuk and Vinagre [19] are compared with respect to control effects.

The remainder of the present paper is organized as follows. Some basic concepts, such as stability, controllability, and observability for a fractional-order dynamical system, are explained in Section 2. After a viscoelastic damper system is selected as a controlled object in Section 3, the proposed design method for fractional-order LQR control is explained in Section 4. In Section 5, a fractional-order state observer is constructed in order to estimate the fractional-order derivative states, which is required for the realization of fractional-order LQR. A numerical calculation method corresponding to a fractional-order dynamical system is introduced in Section 6. In Section 7, numerical simulation results reveal that the fractional-order LQR, from which feedback gains are derived with the proposed design method, can effectively suppress vibrations that occur in the viscoelastic damper system. Finally, concluding remarks are summarized in Section 8.

2. Background

2.1. Stability of Fractional-Order Dynamical System

Denoting the differential operator as D, a linear fractional-order autonomous system can be given as follows:

$$D^{\alpha }\mathit{x}\left(t\right)=\mathit{A}\mathit{x}\left(t\right), \mathit{x}\left(0\right)={\mathit{x}}_0,$$

(1)

where $\alpha$ is the non-integer order of differentiation, $\mathit{x}\left(t\right)$ is the state vector and $\mathit{A}$ is the system matrix.

The dynamical behavior and stability of a fractional-order dynamical system can be analyzed by investigating the poles of the system on the complex plane (s-plane). However, in the case of a fractional-order dynamical system, such as Equation (1), the poles need to be investigated on a new s^α-plane. In the present paper, the s^α-plane is called the ω-plane. Here, ω can be defined as follows:

$$\omega =s^{\alpha },$$

(2)

where s is given as:

$$s=r{e}^{i\theta }.$$

(3)

Then, the definitions given by Equations (2) and (3) yield the following:

$$\omega =s^{\alpha }={\left(r{e}^{i\theta }\right)}^{\alpha }=r^{\alpha }{e}^{i\alpha \theta }.$$

(4)

Using Equation (4), lines on the s-plane can be mapped to corresponding lines on the ω-plane. Practically, let us consider the case in which $\alpha =0.5$. The imaginary axis on the s-plane is expressed as $\theta =\pm \pi /2$, and $0<r<+\infty$. Then, the imaginary axis on the ω-plane can be expressed as follows:

$$\omega =r^{\frac{1}{2}}{e}^{\pm i\frac{\pi }{4}}.$$

(5)

Furthermore, the negative part of the real axis on the s-plane is expressed as $\theta =\pm \pi$, and $0<r<+\infty$. Then, the negative part of the real axis on the ω-plane can be expressed as follows:

$$\omega =r^{\frac{1}{2}}{e}^{\pm i\frac{\pi }{2}}.$$

(6)

Thus, the imaginary axis and the negative part of the real axis can be mapped on the ω-plane as illustrated in Figure 1.

Classification of the dynamical behavior and stability in each region in Figure 1 is as follows [24]: unstable in region ①, damped oscillation in region ②, and overdamping in region ③. Therefore, the necessary and sufficient condition for the linear autonomous system given in Equation (1) to be asymptotically stable for an arbitrary initial value ${\mathit{x}}_0$ is that the arguments of the eigenvalues ${\lambda }_i$ of the system matrix $\mathit{A}$ satisfy the following inequality [25,26]:

$$\left|\arg \left({\lambda }_i\right)\right|>\frac{\alpha \pi }{2}, (i=1, 2, \cdots n),$$

(7)

where n is the order of the system.

2.2. Controllability of Fractional-Order Dynamical System

Next, a linear fractional-order non-autonomous system can be given as follows:

$$\{\begin{array}{c}\begin{array}{cc}{D}^{\alpha }\mathit{x}\left(t\right)=\mathit{A}\mathit{x}\left(t\right)+\mathit{B}u\left(t\right),& \mathit{x}\left(0\right)={\mathit{x}}_0\end{array}\\ y\left(t\right)=\mathit{C}\mathit{x}\left(t\right)\end{array},$$

(8)

where $\alpha$ designates the non-integer order of differentiation, $\mathit{x}\left(t\right)$ is the state vector, $u\left(t\right)$ is the input, $y\left(t\right)$ is the output, $\mathit{A}$ is the system matrix, $\mathit{B}$ is the input matrix, and $\mathit{C}$ is the output matrix.

For an arbitrary initial state $\mathit{x}\left(0\right)={\mathit{x}}_0$, a time $t_f>0$, and a desired value ${\mathit{x}}_f$, if an input $u\left(t\right)$ that enables a solution to satisfy $\mathit{x}\left(t_f\right)={\mathit{x}}_f$ exists, then the system described by Equation (8) is controllable. In addition, the propositions shown below are equivalent [25,26].

• A pair of matrices $\left(\mathit{A}, \mathit{B}\right)$ is controllable.
• The rank of the following controllability matrix ${\mathit{N}}_c$ is n:

  $${\mathit{N}}_c=\left[\mathit{B} \mathit{A}\mathit{B} \cdots {\mathit{A}}^{n-1}\mathit{B}\right].$$

  (9)

2.3. Observability of Fractional-Order Dynamical System

If the initial state $\mathit{x}\left(0\right)={\mathit{x}}_0$ can be determined based on the time responses of the input $u\left(t\right)$ and the output $y\left(t\right)$ at an arbitrary time $t_1>0$, then the system defined by Equation (8) is observable. In addition, the propositions given below are equivalent [25,26].

• A pair of matrices $\left(\mathit{C}, \mathit{A}\right)$ is observable.
• The rank of the following observability matrix ${\mathit{N}}_o$ is n:

  $${\mathit{N}}_o={\left[{\mathit{C}}^T {\left(\mathit{C}\mathit{A}\right)}^T \cdots {\left(\mathit{C}{\mathit{A}}^{n-1}\right)}^T\right]}^T.$$

  (10)

3. Fractional-Order Vibratory System with Viscoelasticity

First, a viscoelastically damped structure is considered as an example of a fractional-order vibratory system in order to explain LQR control. The equation of motion for forced vibration of a one degree-of-freedom viscoelastically damped system with a fractional derivative is:

$$\left[m{D}^2+c{D}^q+k{D}^0\right]x\left(t\right)=f\left(t\right),$$

(11)

where $x\left(t\right)$ is the displacement, m is the mass, c is the viscoelastic damping coefficient, k is the spring constant, t is time, and $f\left(t\right)$ is the external force. Dividing both sides of Equation (11) by m leads to the following normal form:

$$\left[D^2+2\zeta {\omega }_n{}^{2-q}{D}^q+{\omega }_n{}^2{D}^0\right]x\left(t\right)=u\left(t\right),$$

(12)

where ${\omega }_n{}^2= \frac{k}{m}$, $\zeta = \frac{c}{2m{\omega }_n^{2-q}}$ and $u\left(t\right)=\frac{1}{m}f\left(t\right)$.

It is known that q has a value of approximately 0.5 for most viscoelastic materials [27], and is therefore fixed at 0.5 in this paper. Consequently, the viscoelastic damping model can be expressed using a one degree-of-freedom equation of motion known as the generalized Voigt model [28]:

$$\left[D^2+2\zeta {\omega }_n{}^{\frac{3}{2}}{D}^{\frac{1}{2}}+{\omega }_n{}^2{D}^0\right]x\left(t\right)=u\left(t\right),$$

(13)

where $u\left(t\right)$ is the control input.

Next, based on the above equation, a fractional-order state equation is derived. Using a sequence of fractional differentiations with a difference in order of 1/2, the state vector $\mathit{x}\left(t\right)$ can be defined as follows:

$$\mathit{x}\left(\mathrm{t}\right)={\left[\begin{array}{cc}x\left(t\right)& D^{\frac{1}{2}}x\left(t\right)\end{array}\hspace{1em}\begin{array}{cc}{D}^{\frac{2}{2}}x\left(t\right)& D^{\frac{3}{2}}x\left(t\right)\end{array}\right]}^T={\left[\begin{array}{cc}\begin{array}{ccc}{x}_1\left(t\right)& x_2\left(t\right)& x_3\left(t\right)\end{array}& x_4\left(t\right)\end{array}\right]}^T.$$

(14)

Consequently, Equations (13) and (14) lead to the following fractional-order state equation:

$$D^{\frac{1}{2}}\mathit{x}\left(t\right)=\mathit{A}\mathit{x}\left(t\right)+\mathit{B}u\left(t\right),$$

(15)

where

$$\mathit{A}=\left[\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\\ -{\omega }_n^2& -2\zeta {\omega }_n^{\frac{3}{2}}& 0& 0\end{array}\right], \mathit{B}=\left[\begin{array}{c}0\\ 0\\ 0\\ 1\end{array}\right].$$

(16)

4. Fractional-Order LQR Method for State Feedback Control

4.1. Design Method for Fractional-Order LQR

In this section, a fractional-order LQR is designed for the fractional-order vibratory system described in the previous section. The control input is assumed to be given as follows:

$$u\left(t\right)=-\mathit{F}\mathit{x}\left(t\right),$$

(17)

where F designates the feedback-gain matrix.

First, the fractional-order state equation is forced to transform into an integer-order state equation. Setting the closed-loop system matrix as ${\mathit{A}}_{cl}=\mathit{A}-\mathit{B}\mathit{F}$, Equation (17) is substituted into Equation (15), and the following relation is obtained:

$$D^{\frac{1}{2}}\mathit{x}\left(t\right)={\mathit{A}}_{cl}\mathit{x}\left(t\right).$$

(18)

The operator $D^{\frac{1}{2}}$ is applied to both sides of Equation (18) once, and then the following equation is obtained [29]:

$$D\mathit{x}\left(t\right)={\mathit{A}}_{cl}^2\mathit{x}\left(t\right).$$

(19)

Second, a quadratic cost function J is considered for the controlled object modeled by Equation (15):

$$J={{\displaystyle \int }}_0^{\infty }\left[\mathit{x}{\left(t\right)}^T\mathit{Q}\mathit{x}\left(t\right)+u{\left(t\right)}^T\mathit{R}u\left(t\right)\right]dt,$$

(20)

where the weighting matrices $\mathit{Q}\in {\mathbb{R}}^{n\times n}$ and $\mathit{R}\in {\mathbb{R}}^{r\times r}$ are positive semi-definite and positive definite, respectively, and the pair of matrices (${\mathit{Q}}^{\frac{1}{2}}$, A) is observable. Moreover, Equation (17) is substituted into Equation (20) in order to obtain the following:

$$J={{\displaystyle \int }}_0^{\infty }\mathit{x}{\left(t\right)}^T\left[\mathit{Q}+{\mathit{F}}^T\mathit{R}\mathit{F}\right]\mathit{x}\left(t\right)dt.$$

(21)

Modifying Equation (21) using both Equation (19) and the Lyapunov equation yields the following equations:

$$J=\mathit{x}{\left(0\right)}^T\mathit{P}\mathit{x}\left(0\right),$$

(22)

$$\mathit{P}{\mathit{A}}_{cl}^2+{\left({\mathit{A}}_{cl}^2\right)}^T\mathit{P}=-\mathit{Q}-{\mathit{F}}^T\mathit{R}\mathit{F},$$

(23)

where P is a positive-definite symmetric solution that satisfies the algebraic Riccati equation. Here, using Lagrange’s method for undetermined multipliers, matrices P and F, which minimize Equation (22), can be obtained under the constraint condition of Equation (23). If Equation (23) is assumed to be expressed as $\mathit{G}\left(n\times n\right)=\left\{g_{ij}\right\}=\mathbf{0}$, then the number of undetermined multipliers is required to be $n^2$. Therefore, the undetermined multipliers are arranged in the form of $\mathit{K}\left(n\times n\right)=\left\{k_{ij}\right\}$. As a result, the Lagrangian function can be expressed as follows [29]:

$$L=J+{{\displaystyle \sum }}_{i=1}^n{{\displaystyle \sum }}_{j=1}^n{k}_{ij}{g}_{ij}=J+tr\left[{\mathit{K}}^T\mathit{G}\right]=\mathit{x}{\left(0\right)}^T\mathit{P}\mathit{x}\left(0\right)+tr\left[{\mathit{K}}^T\left\{\mathit{P}{\mathit{A}}_{cl}^2+{\left({\mathit{A}}_{cl}^2\right)}^T\mathit{P}+\mathit{Q}+{\mathit{F}}^T\mathit{R}\mathit{F}\right\}\right],$$

(24)

where $tr$ designates the trace of the matrix.

Third, the Lagrangian function $L$ is differentiated with respect to the matrices $\mathit{F}$, $\mathit{P}$, and $\mathit{K}$, respectively, in order to obtain each extreme value. Finally, the equations to obtain the matrices $\mathit{F}$, $\mathit{P}$, and $\mathit{K}$ are derived.

However, it has been clarified that the feedback-gain matrix $\mathit{F}$ cannot be obtained based on the above-mentioned method reported in a previous study [29]. This is because the equation resulting from differentiation directly with respect to F becomes too complicated to be appropriate for obtaining the optimal feedback-gain matrix ${\mathit{F}}_{opt}$.

4.2. Iteration-Based Method for Obtaining Optimal Feedback Gains

Therefore, in order to overcome this difficulty, the approach must be altered to an iteration-based method in order to obtain the optimal feedback-gain matrix ${\mathit{F}}_{opt}$. To do so, ${\mathit{A}}_{cl}^2$ in Equation (24) is expressed as follows:

$${\mathit{A}}_{cl}^2={\left(\mathit{A}-\mathit{B}\mathit{F}\right)}^2=\left(\mathit{A}-\mathit{B}{\mathit{F}}_a\right)\left(\mathit{A}-\mathit{B}{\mathit{F}}_b\right).$$

(25)

Here, ${\mathit{F}}_a$ is properly determined so that all the real parts of the eigenvalues of ${\mathit{A}}_{cl}^2$ are negative, and ${\mathit{F}}_b$ is treated as a variable. By setting ${\mathit{A}}_{cl\_a}=\mathit{A}-\mathit{B}{\mathit{F}}_a$, Equation (24) can be rewritten using Equation (25) as follows:

$$L=\mathit{x}{\left(0\right)}^T\mathit{P}\mathit{x}\left(0\right)+tr\left[{\mathit{K}}^T\left\{\mathit{P}{\mathit{A}}_{cl\_a}\left(\mathit{A}-\mathit{B}{\mathit{F}}_b\right)+{\left[{\mathit{A}}_{cl\_a}\left(\mathit{A}-\mathit{B}{\mathit{F}}_b\right)\right]}^T\mathit{P}+\mathit{Q}+{\mathit{F}}_b^T\mathit{R}{\mathit{F}}_b\right\}\right].$$

(26)

As a result, differentiating the above equation with respect to F_b, P, and K, respectively, produces the following relationships [30]:

$$\frac{\partial L}{\partial \mathit{P}}=\mathit{K}{\left(\mathit{A}-\mathit{B}{\mathit{F}}_b\right)}^T{\mathit{A}}_{cl\_a}^T+{\mathit{A}}_{cl\_a}\left(\mathit{A}-\mathit{B}{\mathit{F}}_b\right)\mathit{K}+\mathit{x}\left(0\right)\mathit{x}{\left(0\right)}^T=\mathbf{0},$$

(27)

$$\frac{\partial L}{\partial {\mathit{F}}_b}=-{\left(\mathit{P}{\mathit{A}}_{cl\_a}\mathit{B}\right)}^T\mathit{K}-{\left({\mathit{A}}_{cl\_a}\mathit{B}\right)}^T\mathit{P}{\mathit{K}}^T+{\mathit{R}}^T{\mathit{F}}_b\mathit{K}+\mathit{R}{\mathit{F}}_b{\mathit{K}}^T=\mathbf{0},$$

(28)

$$\frac{\partial L}{\partial \mathit{K}}=\mathit{P}{\mathit{A}}_{cl\_a}\left(\mathit{A}-\mathit{B}{\mathit{F}}_b\right)+{\left(\mathit{A}-\mathit{B}{\mathit{F}}_b\right)}^T{\mathit{A}}_{cl\_a}^T\mathit{P}+\mathit{Q}+{\mathit{F}}_b^T\mathit{R}{\mathit{F}}_b=\mathbf{0}.$$

(29)

From Equations (27) and (28), the following equation is obtained because the matrices R, P, and K are symmetric:

$${\mathit{F}}_b={\mathit{R}}^{-1}{\mathit{B}}^T{\mathit{A}}_{cl\_a}^T\mathit{P}.$$

(30)

Next, Equations (29) and (30) generate the following algebraic Riccati equation:

$$\mathit{P}{\mathit{A}}_{cl\_a}\mathit{A}+{\left({\mathit{A}}_{cl\_a}\mathit{A}\right)}^T\mathit{P}-\mathit{P}{\mathit{A}}_{cl\_a}\mathit{B}{\mathit{R}}^{-1}{\mathit{B}}^T{\mathit{A}}_{cl\_a}^T\mathit{P}+\mathit{Q}=\mathbf{0}.$$

(31)

In the end, an iterative method to obtain the optimal feedback-gain matrix ${\mathit{F}}_{opt}$ can be explained as follows. Here, P can be obtained from Equation (31), and F_b can be determined by substituting the obtained P into Equation (30). Next, the determined F_b is treated as a new F_a, and Equations (25)–(31) are repeatedly calculated until the feedback-gain matrix converges to ${\mathit{F}}_b={\mathit{F}}_a$. Eventually, the converged feedback-gain matrix can be considered to be the optimal feedback-gain matrix ${\mathit{F}}_{opt}$ to minimize the quadratic cost function J.

5. Fractional-Order LQR Method for Output Feedback Control

5.1. Necessity of Fractional-Order State Observer

The fractional-order LQR control is a kind of state feedback control. Therefore, estimation of the states is an important factor in state feedback control. Numerical simulation of the control may be conducted under the condition that all of the states of the system are assumed to be detected. However, such a condition is not effective in many real cases. Even though sensors are present in the controlled system, obtaining all of the states of the system in most cases is difficult. Furthermore, the viscoelastic damper system, as the controlled object in the present study, contains a fractional-order derivative term that cannot be detected with ordinary sensors in general. Therefore, in order to carry out fractional-order LQR control in practice, it is necessary to estimate the fractional-order derivative states based on information about the observable states.

Accordingly, a scheme that enables estimation of the remaining states based on the input signals and the output signals measured by the sensors is called a state observer. The state observer can be designed as long as the system is observable. If all of the states of the controlled system are estimated using the state observer, then the state feedback control can be conducted using the estimated states. Differences between the estimated states and the actual states are called estimation errors. When the estimation errors become zero, the estimated states match the actual states. Therefore, the state observer must be constructed so that the estimation errors converge to zero.

The estimation errors for the fractional-order system asymptotically go to zero by means of the following design of the state observer:

$$D^{\alpha }\widehat{\mathit{x}}\left(t\right)=\mathit{A}\widehat{\mathit{x}}\left(t\right)+\mathit{B}u\left(t\right)+\mathit{H}\left(y\left(t\right)-\mathit{C}\widehat{\mathit{x}}\left(t\right)\right),$$

(32)

where $\widehat{\mathit{x}}\left(t\right)$, $\mathit{H}$, and $\mathit{C}$ are the estimated state vector, the state observer gain, and the output matrix for the controlled system, respectively. The state observer gain $\mathit{H}$ is selected in order to satisfy the following inequality [31]:

$$\left|\arg \left({\lambda }_{so\_i}\right)\right|>\frac{\alpha \pi }{2}, (i=1, 2, \cdots n),$$

(33)

where ${\lambda }_{so\_i}$ designates an eigenvalue of the matrix ${\mathit{A}}_{so}=\mathit{A}-\mathit{H}\mathit{C}$.

5.2. Configuration of Fractional-Order State Observer

Next, using Equation (32), a state observer is configured for the viscoelastic damper system. First, the viscoelastic damper, as the controlled object, is assumed to be given by the state equation of Equations (15) and (16). In addition, setting $y\left(t\right)$ as the observable, the output equation for the viscoelastic damper system is assumed to be expressed as follows:

$$y\left(t\right)=\mathit{C}\mathit{x}\left(t\right),$$

(34)

where

$$\mathit{C}=\left[\begin{array}{cc}{c}_1& c_2\end{array}\hspace{1em}\begin{array}{cc}{c}_3& c_4\end{array}\right]$$

(35)

in which the values of $c_1$, $c_2$, $c_3$, and $c_4$ are determined by what is observable. Here, if the pair of matrices $\left(\mathit{C}, \mathit{A}\right)$ is observable, then a fractional-order state observer takes the following form:

$$D^{0.5}\widehat{\mathit{x}}\left(t\right)=\mathit{A}\widehat{\mathit{x}}\left(t\right)+\mathit{B}u\left(t\right)+\mathit{H}\left(y\left(t\right)-\mathit{C}\widehat{\mathit{x}}\left(t\right)\right),$$

(36)

$$u\left(t\right)=-\mathit{F}\widehat{\mathit{x}}\left(t\right).$$

(37)

The subtraction of Equation (15) from Equation (36) yields the fractional-order differential equation for the estimation error $\mathit{e}\left(t\right)$:

$$D^{0.5}\mathit{e}\left(t\right)=\left(\mathit{A}-\mathit{H}\mathit{C}\right)\mathit{e}\left(t\right),$$

(38)

where

$$\mathit{e}\left(t\right)=\widehat{\mathit{x}}\left(t\right)-\mathit{x}\left(t\right).$$

(39)

From Equation (38), the poles of the matrix $\left(\mathit{A}-\mathit{H}\mathit{C}\right)$ can be moved freely by the state observer gain $\mathit{H}$, and the convergence characteristics of the estimation error $\mathit{e}\left(t\right)$ can be determined freely using the pole assignment method.

Moreover, Equations (15), (34), (36), (37), and (39) produce the following relationships:

$$\begin{gathered} D^{\frac{1}{2}}\mathit{x}\left(t\right)=\mathit{A}\mathit{x}\left(t\right)+\mathit{B}u\left(t\right)=\mathit{A}\mathit{x}\left(t\right)+\mathit{B}\left(-\mathit{F}\widehat{\mathit{x}}\right)=\mathit{A}\mathit{x}\left(t\right)-\mathit{B}\mathit{F}\left(\mathit{e}\left(t\right)+\mathit{x}\left(t\right)\right) \\ =\left(\mathit{A}-\mathit{B}\mathit{F}\right)\mathit{x}\left(t\right)-\mathit{B}\mathit{F}\mathit{e}\left(t\right), \end{gathered}$$

(40)

$$\begin{gathered} D^{\frac{1}{2}}\mathit{e}\left(t\right)=D^{\frac{1}{2}}\widehat{\mathit{x}}\left(t\right)-D^{\frac{1}{2}}\mathit{x}\left(t\right)=\mathit{A}\widehat{\mathit{x}}\left(t\right)+\mathit{B}u\left(t\right)+\mathit{H}\left(y\left(t\right)-\mathit{C}\widehat{\mathit{x}}\right)-\mathit{A}\mathit{x}\left(t\right)-\mathit{B}u\left(t\right) \\ =\mathit{A}\left(\widehat{\mathit{x}}\left(t\right)-\mathit{x}\left(t\right)\right)-\mathit{H}\left(\mathit{C}\widehat{\mathit{x}}\left(t\right)-\mathit{C}\mathit{x}\left(t\right)\right)=\left(\mathit{A}-\mathit{H}\mathit{C}\right)\left(\widehat{\mathit{x}}\left(t\right)-\mathit{x}\left(t\right)\right) \\ =\left(\mathit{A}-\mathit{H}\mathit{C}\right)\mathit{e}\left(t\right). \end{gathered}$$

(41)

Equations (40) and (41) generate the state equation for the augmented system that includes the state observer as follows [32]:

$$\left[\begin{array}{c}{D}^{\frac{1}{2}}\mathit{x}\left(t\right)\\ D^{\frac{1}{2}}\mathit{e}\left(t\right)\end{array}\right]=\left[\begin{array}{cc}\mathit{A}-\mathit{B}\mathit{F}& -\mathit{B}\mathit{F}\\ \mathbf{0}& \mathit{A}-\mathit{H}\mathit{C}\end{array}\right]\left[\begin{array}{c}\mathit{x}\left(t\right)\\ \mathit{e}\left(t\right)\end{array}\right].$$

(42)

Figure 2 is a block diagram of the augmented system, in which the state observer is incorporated with the feedback control system.

6. Numerical Calculation Method for Fractional-Order Dynamical System

6.1. Numerical Solution of Fractional-Order State Equation

Since a fractional-order state equation contains fractional-order derivative terms, normal numerical calculation methods, such as the Euler method, cannot be used to simulate this equation. Therefore, several numerical calculation approaches corresponding to a differential equation and a state equation that include fractional-order derivative terms have already been proposed [33,34,35,36]. One of the numerical calculation methods corresponding to a fractional-order state equation is introduced in this section [35,36].

First, the Grünwald–Letnikov definition, which is a fractional-order differentiation, is as follows:

$${}_a{}^{}D{}_t^{q}f\left(t\right)=\underset{\Delta t\to 0}{\lim }{\left(\Delta t\right)}^{-q}{{\displaystyle \sum }}_{j=0}^{\left[\frac{t-a}{\Delta t}\right]}{\left(-1\right)}^j\left(\begin{array}{c}q\\ j\end{array}\right)f\left(t-j\Delta t\right),$$

(43)

where $\Delta t$ is the pitch width, $q$ is the fractional-order of differentiation, $[\hspace{1em}]$ is the integer part of the argument, and $\left(\begin{array}{c}q\\ j\end{array}\right)$ indicates the generalized binomial coefficient. The revised version of the above definition for the numerical calculation is as follows:

$${}_a{}^{}D{}_t^{q}f\left(t\right)=\underset{\Delta t\to 0}{\lim }\frac{1}{{\left(\Delta t\right)}^q}{{\displaystyle \sum }}_{j=0}^{\left[\frac{t-a}{\Delta t}\right]}{w}_j^{\left(q\right)}f\left(t-j\Delta t\right),$$

(44)

where

$$w_0^{\left(q\right)}=1, w_j^{\left(q\right)}=\left(1-\frac{q+1}{j}\right)w_{j-1}^{\left(q\right)}, j=1, 2, \cdots .$$

(45)

Next, a simple scalar differential equation is considered as follows:

$$D^{q}z\left(t_k\right)=f\left(t_k, z\left(t_k\right)\right),$$

(46)

where the number of calculations is $k$, and the time is set to $t_k=k\Delta t$. Applying Equation (44) to the left-hand side of the above equation yields [36]:

$$\frac{1}{{\left(\Delta t\right)}^q}{{\displaystyle \sum }}_{j=0}^m{w}_j^{\left(q\right)}z\left(t_{k-j}\right)=\frac{1}{{\left(\Delta t\right)}^q}\left[z\left(t_k\right)+{{\displaystyle \sum }}_{j=1}^m{w}_j^{\left(q\right)}z\left(t_{k-j}\right)\right]=f\left(t_k, z\left(t_k\right)\right),$$

(47)

where $m=\left[\frac{t_k}{\Delta t}\right]=k$.

Furthermore, Equation (47) is transformed so that this equation can be solved in terms of $z$, and then the following relation is obtained:

$$z\left(t_k\right)={\left(\Delta t\right)}^{q}f\left(t_k, z\left(t_k\right)\right)-{{\displaystyle \sum }}_{j=1}^m{w}_j^{\left(q\right)}z\left(t_{k-j}\right).$$

(48)

If the pitch width $\Delta t$ is sufficiently small, then the first term on the right-hand side of Equation (48) can be approximated as follows [36]:

$${\left(\Delta t\right)}^{q}f\left(t_k, z\left(t_k\right)\right) \approx {\left(\Delta t\right)}^{q}f\left(t_k, z\left(t_{k-1}\right)\right).$$

(49)

By this approximation, the numerical calculation becomes possible.

Moreover, if each row of the system matrix $\mathit{A}$ of Equation (8) is dealt with as a function $f\left(t\right)$ in Equation (46), then the above-mentioned numerical calculation method can be used to solve the fractional-order state equation as follows [36]:

$$\{\begin{array}{c}{z}_1\left(t_k\right)={\left(\Delta t\right)}^q{f}_1\left(t_k, z_1\left(t_{k-1}\right), z_2\left(t_{k-1}\right),\cdots , z_n\left(t_{k-1}\right)\right)-{\displaystyle \sum }_{j=1}^m{w}_j^{\left(q\right)}{z}_1\left(t_{k-j}\right)\\ z_2\left(t_k\right)= {\left(\Delta t\right)}^q{f}_2\left(t_k, z_1\left(t_k\right), z_2\left(t_{k-1}\right),\cdots , z_n\left(t_{k-1}\right)\right)- {\displaystyle \sum }_{j=1}^m{w}_j^{\left(q\right)}{z}_2\left(t_{k-j}\right)\\ ⋮\\ z_n\left(t_k\right)= {\left(\Delta t\right)}^q{f}_n\left(t_k, z_1\left(t_k\right), z_2\left(t_k\right),\cdots , z_n\left(t_{k-1}\right)\right)- {\displaystyle \sum }_{j=1}^m{w}_j^{\left(q\right)}{z}_n\left(t_{k-j}\right)\end{array}.$$

(50)

In the present paper, each state $x_i\left(t_k\right)$ is applied to each $z_i\left(t_k\right).$ In the case that the Caputo definition is used, because of the relationship between the Caputo definition and the Grünwald–Letnikov definition, the initial value for each state must be dealt with as follows:

$$z_i\left(t_k\right)=x_i\left(t_k\right)-x_i\left(0\right).$$

(51)

Equation (51) is substituted into Equation (50), and then Equation (50) can be rewritten as follows:

$$\{\begin{array}{c}{x}_1\left(t_k\right)={\left(\Delta t\right)}^q{f}_1\left(t_k, x_1\left(t_{k-1}\right), x_2\left(t_{k-1}\right),\cdots , x_n\left(t_{k-1}\right)\right)-{\displaystyle \sum }_{j=1}^m{w}_j^{\left(q\right)}\left\{x_1\left(t_{k-j}\right)-x_1\left(0\right)\right\}+x_1\left(0\right)\\ x_2\left(t_k\right)={\left(\Delta t\right)}^q{f}_2\left(t_k, x_1\left(t_k\right), x_2\left(t_{k-1}\right),\cdots , x_n\left(t_{k-1}\right)\right)-{\displaystyle \sum }_{j=1}^m{w}_j^{\left(q\right)}\left\{x_2\left(t_{k-j}\right)-x_2\left(0\right)\right\}+x_2\left(0\right)\\ ⋮\\ x_n\left(t_k\right)={\left(\Delta t\right)}^q{f}_n\left(t_k, x_1\left(t_k\right), x_2\left(t_k\right),\cdots , x_n\left(t_{k-1}\right)\right)-{\displaystyle \sum }_{j=1}^m{w}_j^{\left(q\right)}\left\{x_n\left(t_{k-j}\right)-x_n\left(0\right)\right\}+x_n\left(0\right)\end{array}.$$

(52)

Moreover, when the amount of past calculation data is large, the calculation of the second term on the right-hand side of each row in Equation (52) becomes large. Therefore, the parameter $L_0$, which determines the truncated term in the past calculation data, is introduced so that the second term on the right-hand side of each row in Equation (52) is made to be approximated as follows [36]:

$${{\displaystyle \sum }}_{j=1}^m{w}_j^{\left(q\right)}\left\{x_i\left(t_{k-j}\right)-x_i\left(0\right)\right\} \approx {{\displaystyle \sum }}_{j=1}^{\min \left(L_0, k\right)}{w}_j^{\left(q\right)}\left\{x_i\left(t_{k-j}\right)-x_i\left(0\right)\right\}.$$

(53)

Using the numerical calculation method described above, the fractional-order state equation and the fractional-order state observer can be simulated numerically.

6.2. Comparison between Numerical and Exact Solutions

Next, the solution obtained by the above-mentioned numerical calculation method is compared with the exact solution. The exact solution of Equation (15) can be given as follows using the Mittag–Leffler function ${\mathit{E}}_{a_1, a_2}$ [29,37]:

$$\mathit{x}\left(t\right)={\mathit{E}}_{\frac{1}{2}, 1}\left(\mathit{A}{t}^{\frac{1}{2}}\right)\mathit{x}\left(0\right)+{{\displaystyle \int }}_0^t{\left(t-\tau \right)}^{\frac{1}{2}-1}{\mathit{E}}_{\frac{1}{2}, \frac{1}{2}}\left(\mathit{A}{\left(t-\tau \right)}^{\frac{1}{2}}\right)\mathit{B}u\left(\tau \right)d\tau ,$$

(54)

$${\mathit{E}}_{a_1, a_2}\left(z\right)={{\displaystyle \sum }}_{k=0}^{\infty }\frac{z^k}{\mathrm{\Gamma }\left(a_{1}k+a_2\right)}, \left(a_1>0, a_2>0\right).$$

(55)

In addition, the parameters for the viscoelastic damper are set to $\zeta =0.1$ and ${\omega }_n=3.0$ rad/s. The initial conditions for the states are set to $\mathit{x}\left(0\right)={\left[\begin{array}{ccc}0& 0& \begin{array}{cc}1& 0\end{array}\end{array}\right]}^T$. For each state, the numerical calculation results and the exact solutions obtained under the free vibration condition are compared in the following figures.

First, Figure 3 shows calculation results for simulation conditions with a larger $\Delta t$ than those in Figure 4. The numerical calculations in Figure 4 are more accurate than those in Figure 3, because the gap between the numerical calculation results and the exact solution in Figure 4 is smaller than that in Figure 3. Based on this result, it is understood that a smaller $\Delta t$ improves the accuracy of numerical calculations. A numerical calculation method that enables a highly accurate simulation of a state equation that includes fractional-order derivative terms has thus been clarified.

7. Illustrative Examples

7.1. Fractional-Order LQR Control for Viscoelastic Damper System

In this section, the viscoelastic damper described in Section 3 is set to be the controlled object. Moreover, the feedback gains of the fractional-order LQR are obtained using the method explained in Section 4, and the states are estimated with the fractional-order state observer described in Section 5. In addition, each state is simulated using the numerical calculation method introduced in Section 6. Another derivation method for the feedback gains of the fractional-order LQR, which is different from the method proposed herein, was proposed by Sierociuk and Vinagre [19]. Accordingly, vibration waveforms for the states controlled by the LQR using the feedback gains obtained with the proposed method are compared with those obtained using the Sierociuk and Vinagre method.

The parameters for the viscoelastic damper system are set as $\zeta =0.1$ and ${\omega }_n=3.0$ rad/s. The observable is chosen to be the displacement. Therefore, the output matrix $\mathit{C}$ in Equation (35) is as follows:

$$\mathit{C}=\left[\begin{array}{ccc}1& 0& \begin{array}{cc}0& 0\end{array}\end{array}\right].$$

(56)

First, the stability of the viscoelastic damper system expressed by Equation (15) is investigated. The eigenvalues of the system matrix $\mathit{A}$ are given as:

$$\lambda =1.2263\pm 1.3098i, -1.2263\pm 1.1366i.$$

(57)

Here, the stability of the system is determined using Equation (7):

$$\frac{2}{\pi }\left|\arg \left({\lambda }_i\right)\right|=0.5210, 0.5210, 1.5241, 1.5241>0.5.$$

(58)

As a result of this inequality, this system can be said to be asymptotically stable.

Next, the controllability and the observability of the viscoelastic damper system are examined. The ranks of the matrices ${\mathit{N}}_c$ and ${\mathit{N}}_o$ in Equations (9) and (10) can be calculated as follows:

$$rank\left({\mathit{N}}_c\right)=4, rank\left({\mathit{N}}_o\right)=4.$$

(59)

The pair of matrices $\left(\mathit{A}, \mathit{B}\right)$ for the viscoelastic damper system is controllable, and the pair of matrices $\left(\mathit{C}, \mathit{A}\right)$ for the viscoelastic damper system is observable. Based on these results, the system can be said to be controllable and observable. Therefore, a state feedback controller and a state observer can be designed for this system.

If the weighting matrices $\mathit{Q}$ and $\mathit{R}$ in the quadratic cost function $J$ are designed as follows:

$$\mathit{Q}=\left[\begin{array}{cccc}10& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 10& 0\\ 0& 0& 0& 0\end{array}\right], \mathit{R}=1,$$

(60)

then the optimal feedback-gain matrix ${\mathit{F}}_{opt}$ is obtained as follows:

$${\mathit{F}}_{opt}=\left[\begin{array}{ccc}0.5394& -0.5520& \begin{array}{cc}3.3285& -0.0098\end{array}\end{array}\right].$$

(61)

In this case, the weights in $\mathit{Q}$ are placed only on the integer-order derivative states. With these feedback gains, the eigenvalues for the closed-loop system-matrix ${\mathit{A}}_{cl}=\mathit{A}-\mathit{B}\mathit{F}$ are shown below:

$${\lambda }_{cl}=0.8475\pm 1.5897i, -0.8426\pm 1.4931i.$$

(62)

From Equation (62), since the poles of the closed-loop system are located in region ② in Figure 1, the system is asymptotically stable.

Next, a fractional-order state observer is designed using Equation (36) by the pole-assignment method. Furthermore, generally, when a state observer and a regulator are used simultaneously, the state observer must be designed to be able to estimate the states of the system faster than the regulator makes the states converge to zeros. Here, the eigenvalues of the matrix ${\mathit{A}}_{so}=\mathit{A}-\mathit{H}\mathit{C}$ are designed as follows:

$${\lambda }_{so}=-10.0, -9.0, -8.0, -7.0.$$

(63)

The stability of these eigenvalues can be determined using Equation (7) as:

$$\frac{2}{\pi }\left|\arg \left({\lambda }_{cl\_i}\right)\right|=0.6882, 0.6882, 1.3271, 1.3271>0.5,$$

(64)

$$\frac{2}{\pi }\left|\arg \left({\lambda }_{so\_i}\right)\right|=2, 2, 2, 2>0.5.$$

(65)

The closed-loop system ${\mathit{A}}_{cl}$ and the state observer ${\mathit{A}}_{so}$ are asymptotically stable. Judging from Equations (62) and (63), the eigenvalues of the matrix ${\mathit{A}}_{so}$ are more stable than those of the matrix ${\mathit{A}}_{cl}$. This proves that state estimation by the state observer can be designed to be faster than state convergence by the regulator. Consequently, the state observer gain matrix $\mathit{H}$ is obtained as follows:

$$\mathit{H}={10}^3\times \left[\begin{array}{ccc}0.0340& 0.4310& \begin{array}{cc}2.4130& 4.9957\end{array}\end{array}\right].$$

(66)

In practice, in the simulation, the initial conditions for the states and those for the estimated values are assumed to be as follows:

$$\mathit{x}\left(0\right)={\left[\begin{array}{ccc}0& 0& \begin{array}{cc}1& 0\end{array}\end{array}\right]}^T,$$

(67)

$$\widehat{\mathit{x}}\left(0\right)={\left[\begin{array}{ccc}0& 0& \begin{array}{cc}1.2& 0\end{array}\end{array}\right]}^T.$$

(68)

In addition, the parameter in the numerical calculation method is set as $\Delta t=0.0005$.

Using the above-mentioned values, whether the fractional-order state observer can estimate all of the states is first investigated under the free vibration condition. Figure 5 shows each estimation error, and Figure 6 compares each actual state and each estimated state. As shown in Figure 5, each estimation error converges to zero, although the estimation error about $x_4$ converges more slowly than other errors. Figure 6 demonstrates that each state can be estimated successfully by the fractional-order state observer.

Next, fractional-order LQR control is conducted. In comparing the states with control and those without control, the gaps between the actual states and the estimated values are shown in Figure 7, and the state evolution is illustrated in Figure 8. The results with control are obtained with output feedback control using a fractional-order LQR.

Figure 7 confirms that the estimation of the states by the fractional-order state observer is performed successively. Judging from the results in Figure 8, compared with the no-control case, good damping effects are confirmed in the case that fractional-order LQR control is applied.

Consequently, using the feedback gains obtained with the iteration-based method and the fractional-order state observer, it is proven that fractional-order LQR control can be carried out for a controlled object that includes fractional-order derivative states in its state equation.

7.2. Comparison between Proposed Method and Sierociuk and Vinagre Method

Moreover, the control results for the fractional-order LQR using the feedback gains obtained using the proposed method and the fractional-order LQR using the feedback gains obtained with the Sierociuk and Vinagre method [19] are compared. Using the same weighting matrices, $\mathit{Q}$ and $\mathit{R}$ in Equation (60), the feedback gains are calculated with the Sierociuk and Vinagre method as follows:

$${\mathit{F}}_{SV}=\left[\begin{array}{ccc}0.5394& 14.7824& \begin{array}{cc}13.0640& 5.1116\end{array}\end{array}\right].$$

(69)

The poles of the closed-loop system-matrix $\mathit{A}-\mathit{B}{\mathit{F}}_{SV}$ are as follows:

$$Poles=\{\begin{array}{c}-1.5649\pm 1.5733i\\ -0.9908\pm 0.9774i\end{array}.$$

(70)

These poles are located in region ③ in Figure 1. As a result, the controlled states with the feedback gains in Equation (61) and the controlled states with the feedback gains in Equation (69) can be compared, as shown in Figure 9. The results with control are obtained with output feedback control using a fractional-order LQR.

From Figure 9, it is clear that both feedback gains exhibit good control effects. In particular, in the vibration waveforms controlled using the feedback gains obtained with the Sierociuk and Vinagre method, all of the states are confirmed to be overdamped. However, it is obvious that the convergence of $x_1$ (displacement) to the target value is considerably slow.

Next, the design of the weighting matrices $\mathit{Q}$ and $\mathit{R}$ in the quadratic cost function J is modified as follows:

$$\mathit{Q}=\left[\begin{array}{cccc}10& 0& 0& 0\\ 0& 10& 0& 0\\ 0& 0& 10& 0\\ 0& 0& 0& 10\end{array}\right], \mathit{R}=1.$$

(71)

In this case, the weights in $\mathit{Q}$ are placed on not only the integer-order derivative states, but also the fractional-order derivative states. Using weighting matrices $\mathit{Q}$ and $\mathit{R}$, the feedback gain matrix derived using the proposed method is as follows:

$$\mathit{F}=\left[\begin{array}{ccc}0.5394& -2.3077& \begin{array}{cc}3.7258& -0.9831\end{array}\end{array}\right].$$

(72)

The poles of the closed-loop system matrix $\mathit{A}-\mathit{B}\mathit{F}$ are as follows:

$$Poles=\{\begin{array}{c}-1.0806\pm 1.5848i\\ -0.5891\pm 1.4985i\end{array}.$$

(73)

These poles are located in region ② in Figure 1. On the other hand, the feedback gain matrix obtained by the Sierociuk and Vinagre method is as follows:

$${\mathit{F}}_{SV}=\left[\begin{array}{ccc}0.5394& 16.2648& \begin{array}{cc}15.1137& 6.3425\end{array}\end{array}\right].$$

(74)

The poles of the closed-loop system matrix $\mathit{A}-\mathit{B}{\mathit{F}}_{SV}$ are as follows:

$$Poles=\{\begin{array}{c}-2.3380\pm 0.3989i\\ -0.8332\pm 1.0007i\end{array}.$$

(75)

These poles are located in region ③ in Figure 1. The states controlled with the above-mentioned feedback-gain matrices are compared in Figure 10. The results with control are obtained with output feedback control using a fractional-order LQR.

Figure 10 demonstrates that the results of the states are oscillations damped by the feedback-gain matrix $\mathit{F}$ and are overdamped by the feedback-gain matrix ${\mathit{F}}_{SV}$. Furthermore, the convergence property of $x_1$ in Figure 10 is similar to that in Figure 9. However, in the proposed method, the vibration suppression effect becomes worse when the weights in $\mathit{Q}$ are placed on the fractional-order derivative states. On the other hand, in the Sierociuk and Vinagre method, the weights on the fractional-order derivative states in $\mathit{Q}$ can play an important role for better vibration control effects.

8. Conclusions

The present study investigated whether the LQR control method could be applied to a state equation with fractional-order derivative states. A design method for a fractional-order LQR control system was proposed in order to deal with a controlled object that included fractional-order derivative terms in the equation of motion.

First, in order to explain the design method, a viscoelastic damper was modeled by an equation of motion including a fractional-derivative term. The equation of motion was able to be transformed into a fractional-order state equation.

Second, a fractional-order LQR was applied to control the dynamics of the viscoelastically damped structure. A design method for fractional-order LQR control was proposed using the Lyapunov equation and Lagrange’s method for undetermined multipliers. A new iteration-based approach to solve the algebraic Riccati equation was established in order to obtain the feedback gains for the fractional-order LQR.

Third, a fractional-order state observer was constructed in order to estimate the fractional-order derivative states. All states for the viscoelastic damper system were successfully estimated using the fractional-order state observer.

Fourth, a numerical calculation algorithm corresponding to a differential equation containing fractional derivative terms was introduced. The numerical calculation algorithm was based on the Grünwald–Letnikov definition of fractional differentiation. Numerical calculation results using the introduced method and the exact solution of the viscoelastic damper model were compared in order to investigate its calculation accuracy.

Finally, using the introduced numerical simulation method, LQR control effect simulations were carried out. The numerical simulations revealed that the iterative method was valid for obtaining the feedback gains for the fractional-order LQR. The fractional-order LQR controller was confirmed to be capable of successfully suppressing vibrations appearing in a vibratory system with viscoelasticity. Moreover, the fractional-order LQR obtained by the Sierociuk and Vinagre method and that derived by the proposed method were compared with respect to their control effects.

However, the fractional-order state observer in the present paper is not optimal, and therefore needs to be improved. In addition, control effects by a fractional-order LQR to suppress vibrations appearing in a vibratory system with viscoelasticity must be compared with control effects by an integer-order LQR.