# Fractional-Order LQR and State Observer for a Fractional-Order Vibratory System
## 6 Numerical Calculation Method for Fractional-Order Dynamical System  
### 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 Grunwald-Letnikov definition, which is a fractional-order differentiation, is as follows:  
\[{}_{a}D_{t}^{q}f(t)=\lim_{\Delta t\to 0}(\Delta t)^{-q}\sum_{j=0}^{[\frac{t-q}{ \Delta t}]}(-1)^{j}\bigg{(} \begin{array}{c}q\\ j\end{array} \bigg{)}f(t-j\Delta t), \tag{43}\]  
where \(\Delta t\) is the pitch width, \(q\) is the fractional-order of differentiation, \([\ ]\) is the integer part of the argument, and \(\bigg{(} \begin{array}{c}q\\ j\end{array} \bigg{)}\) indicates the generalized binomial coefficient. The revised version of the above definition for the numerical calculation is as follows:  
\[{}_{a}D_{t}^{q}f(t)=\lim_{\Delta t\to 0}\frac{1}{(\Delta t)^{q}}\sum_{j=0}^{[ \frac{t-q}{\Delta t}]}w_{j}^{(q)}f(t-j\Delta t), \tag{44}\]  
where  
\[w_{0}^{(q)}=1,\ w_{j}^{(q)}=\bigg{(}1-\frac{q+1}{j}\bigg{)}w_{j-1}^{(q)},\ j=1,\ 2,\ \cdots. \tag{45}\]  
Next, a simple scalar differential equation is considered as follows:  
\[D^{q}z(t_{k})=f(t_{k},\,z(t_{k})), \tag{46}\]  
Figure 2: Block diagram of augmented system. 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}}\sum\limits_{j=0}^{m}w_{j}^{(q)}z\Big{(}t_{ k-j}\Big{)}=\frac{1}{\left(\Delta t\right)^{q}}\bigg{[}z(t_{k})+\sum\limits_{j=1}^{m}w_{j} ^{(q)}z\Big{(}t_{k-j}\Big{)}\bigg{]}=f(t_{k},\ z(t_{k})), \tag{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(t_{k})=(\Delta t)^{q}f(t_{k},\ z(t_{k}))-\sum\limits_{j=1}^{m}w_{j}^{(q)}z \Big{(}t_{k-j}\Big{)}. \tag{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]:  
\[(\Delta t)^{q}f(t_{k},\ z(t_{k}))\ \approx(\Delta t)^{q}f(t_{k},\ z(t_{k-1})). \tag{49}\]  
By this approximation, the numerical calculation becomes possible. Moreover, if each row of the system matrix \(\mathbf{A}\) of Equation (8) is dealt with as a function \(f(t)\) in Equation (46), then the above-mentioned numerical calculation method can be used to solve the fractional-order state equation as follows [36]:  
\[\left\{ \begin{array}{l}z_{1}(t_{k})=(\Delta t)^{q}f_{1}(t_{k},\ z_{1}(t_{k-1 }),\ z_{2}(t_{k-1}),\cdots,\ z_{n}(t_{k-1}))-\sum\limits_{j=1}^{m}w_{j}^{(q)}z _{1}\Big{(}t_{k-j}\Big{)}\\ z_{2}(t_{k})=\ (\Delta t)^{q}f_{2}(t_{k},\ z_{1}(t_{k}),\ z_{2}(t_{k-1}),\cdots,\ z_{n}(t_{k-1 }))-\sum\limits_{j=1}^{m}w_{j}^{(q)}z_{2}\Big{(}t_{k-j}\Big{)}\\ \vdots\\ z_{n}(t_{k})=\ (\Delta t)^{q}f_{n}(t_{k},\ z_{1}(t_{k}),\ z_{2}(t_{k}),\cdots,\ z_{n}(t_{k-1 }))-\sum\limits_{j=1}^{m}w_{j}^{(q)}z_{n}\Big{(}t_{k-j}\Big{)}\end{array} .\right. \tag{50}\]  
In the present paper, each state \(x_{i}(t_{k})\) is applied to each \(z_{i}(t_{k})\).