Minimal realization and approximation of commensurate linear fractional-order systems via Loewner matrix method

1.   Introduction

Creating mathematical model is an essential step for the simulation, analysis, control and design of modern systems. In the last decades, it has become ubiquitous to build model from collected input-output data by adopting data driven modeling approaches. Loewner Matrix Method (LMM) is a widely applicable data-driven modeling approach which is introduced and developed in ^[1,2,3,4]. LMM provides a systematic framework for constructing models from given noise-free measurements in the frequency or time domain ^[5]. In fact, recent investigation shows that LMM is also able to approximate the original system even for high levels of noise ^[6,7]. To make the LMM more compatible for dealing with dynamical systems with intrinsic structural properties like delay or second-order, Schulze et al. recently generalized the LMM to generate structure-preserving model ^[8].

As a natural extension of the integer-order systems, fractional-order systems received a fast increasing attention in the last years. It is revealed that fractional-order model is able to more accurately capture the basic dynamics of many real-life systems arising in electrical, electronic, mechanical, biological (e.g., cardiac tissue electrode interface). Especially in recent years, the fractional-order models were applied to the fields of chemical and bio-sciences engineering in more and more literature. Dulf et al. proposed the use of fractional-order model to represent the complex mechanisms of the biochemical processes without losing the physical meaning of gain and time constants in ^[9] and it worked better than integer-order. Toledo-Hernandez et al. extended fractional calculus to the biological reactive systems in ^[10,11], they shown that the dynamics of some reactive systems displaying atypical behavior can be represented by fractional-order differential equations. Khan proposed a fractional-order biochemical reaction model and shown that the fractional modeling has more advantage than classical integer model ^[12].

Modeling of fractional-order systems via data-driven methods has received a lot of interests. Some researchers focused on the modeling of fractional-order system (FoS) and proposed different methods to identify the parameters and order of FoS in time domain or frequency domain.The differential evolution (DE) algorithm was applied to search the optimal fractional commensurate differential order in ^[13]. Gao proposed a stable model order reduction method for fractional-order systems and achieved a great fitting effect with the original system in ^[14]. The FoS was identified by applying least squares method in ^[15]. In the paper ^[16], a subspace identification algorithm in the time-domain was proposed to identify the coefficient matrices and the order-$\alpha$ of multi-variable FoS. An algebraic approach was proposed to identify linear systems with fractional derivatives in ^[17].

In the present work, we apply and extend the LMM to construct fractional-order state space model with low computational cost from the interpolation data. It is shown that the generalized LMM is powerful to generate the desired fractional-order model with minimal realization. In particular, the unknown fractional-order $\alpha$ can be identified accurately thanks to its rank revealing property. Casagrande et al. proposed using integer-order model to approximate the FoS in the Loewner framework ^[18]. The drawback of the method is that the approximate system works badly at the high frequencies sometimes. Moreover, the order will increases as the amount of interpolation data increases. Our method can solve both of the problems effectively.

The remainder of this paper is organized as follows. Section 2 not only introduces a brief mathematical background of commensurate FoS and fractional-order time delay systems (FoTDS), but also recalls the generalized Loewner realization method. The associate data-driven realization problems about the fractional-order are introduced. Section 3 present a generalization to the commensurate FoS based on extensions of the LMM (Divided into two cases of order $\alpha$ known and unknown). An approximation method to FoDTS based on the LMM is illustrated in Section 4. To study the applicability of the proposed method, some examples are outlined in Section 5. Section 6 gives a concluding remark and discusses the future works.

2.   Preliminaries and problem statement

2.1.   Commensurate fractional-order systems

Detailed introduction to fractional-order systems is given in ^[19,20,21]. In this work, we only consider the commensurate fractional-order linear time invariant (FoLTI) continuous systems. Generally, the state space model of a FoLTI system is described by

where $x(t)\in \mathbb{R}^n, u(t)\in \mathbb{R}^m$ and $y(t) \mathbb \in \mathbb{R}^p$ are the system states, input and output vectors, respectively. $\mathcal D^{\alpha}$ is the fractional differential operator, $\boldsymbol{A}\in \mathbb{R}^{n\times n}$, $\boldsymbol{B}\in \mathbb{R}^{n\times m}$ and $\boldsymbol{C}\in \mathbb{R}^{p\times n}$ are the system matrices. The state space model described by the Eq (2.1) can be transformed into the following fractional transfer function form.

Problem 1. For the system described by the Eq (2.1), given a set of input-output frequency responses of the transfer function. we have

the right interpolation data:

and the left interpolation data:

where $\boldsymbol G(\lambda_{i}) \mathbf{r}_{i} = \mathbf{w}_{i} \text{ and } \ell_{j} \boldsymbol G(\mu_{j}) = \mathbf{v}_{j}$. Our purpose is to realize a minimal state space model ^[22] $[E, A, B, C]$ of FoS. It is divided into two parts.

a) In case that the value of $\alpha$ is pre-known, the system matrices $[E, A, B, C]$ are constructed according to the interpolation data, such that the commensurate fractional-order transfer function $H(s) = C(s^{\alpha} E-A)^{-1} B$ satisfies the interpolation data. i.e., $H(\lambda_{i})\mathbf{r}_{i} = \mathbf{w}_{i}$ and $\ell_{j}H(\mu_{j}) = \mathbf{v}_{j}$.

b) In case that the value of $\alpha$ is unknown, to find the optimal $\hat \alpha$ and construct the corresponding matrices $[E, A, B, C]$ according to the interpolation data, such that the transfer function $H(s) = C(s^{\hat \alpha} E-A)^{-1} B$ satisfies the right and left interpolation conditions.

2.2.   Fractional-order time delay systems

The transfer function of a SISO FoTDS is given as the following expression:

where $\tau$ is the time delay.

Problem 2. For a FoTDS, given a set of frequency response input-output pairs $(s_i, S_i), \ i = 1, \dots N$, where $S_i$ is obtained by sampling the transfer function Eq (2.5), i.e., $S_i = G(s_i)$. Our purpose is to find a linear integer-order model with r-order in the descriptor form:

where $x(t) \in \mathbb{R}^{r}, u(t) \in \mathbb{R}, y(t) \in \mathbb{R}$, and the associated transfer function $H(s)$ satisfies the interpolation data, i.e..

2.3.   Loewner framework realization

In order to solve the above both of problems. We resort the Loewner framework which has been widely applied to the generalized realization problem by Antoulas and his co-workers ^[1,3,23]. The Loewner framework are extended to the dynamic time delay systems in ^[24].

We briefly recall the Loewner realization in ^[1]. The Loewner matrix $\mathbb{L}$ and shifted Loewner matrix $\mathbb{L}_{\sigma}$ are defined as follows

The interpolation data is

These matrices satisfy the following Sylvester equations

Theorem 2.1. For all $x \in\left\{\lambda_{i}\right\} \cup \{ \mu_{j} \}$, let $\text{det}(x\mathbb{L} - \mathbb{L}_\sigma) \neq 0$. Then $\boldsymbol{E} = -\mathbb{L}, \ \boldsymbol{A} = -\mathbb{L}_\sigma, \ \boldsymbol{B} = \mathbf{V}, \ \boldsymbol{C} = \mathbf{W}$ is a minimal realization. The transfer function is

Theorem 2.2. Suppose that:

For some $x$, the short singular value decomposition (SVD) is computed as follow:

where $\Sigma \in \mathbb{C}^{r \times r}$ is positive definite and diagonal, $\mathbf{Y} \in \mathbb{C}^{r \times k}$ and $\mathbf{X} \in \mathbb{C}^{k \times r}$ are the orthogonal factors of the short SVD. Then a minimal realization is given as follows:

3.   Loewner framework for commensurate FoS

We extend the Loewner framework to the commensurate FoS. The fractional-order Loewner matrix $\mathbb{L}_{f}$ and shifted Loewner matrix $\sigma\mathbb{L}_{f}$ associated with the commensurate FoS are defined as following in terms of the data Eqs (2.3) and (2.4) :

where $\mu_{i}^{\alpha}-\lambda_{j}^{\alpha} \neq 0, \text{for all } i, j = 1, \cdots, k, \alpha$ is the commensurate order.

These matrices satisfy the following Sylvester equations

3.1.   Realization of FoS with pre-known $\alpha$

When we already know the commensurate order-$\alpha$, for the solution of the first part of Problem 1, There are two cases: the right amount of data and the more realistic redundant amount of data. The following theorem gives the solution for the first case.

Theorem 3.1. If det$(x^\alpha\mathbb{L}_{f}-\sigma\mathbb{L}_{f})\neq0, \ $ for all $x \in \{\lambda_{i}\}\ \cup \{\mu_{j} \}, \ i, j = 1, \cdots, k.$ Then $E = -\mathbb{L}_{f}, A = -\sigma\mathbb{L}_{f}, B = \mathbf{V} \mathit{\text{and}}C = \mathbf{W}$ is a minimal realization. The corresponding transfer function is: $H(s) = \mathbf{W}(\sigma\mathbb{L}_{f}-s^{\alpha} \mathbb{L}_{f})^{-1} \mathbf{V}.$

The following proof shows that the realization satisfies the interpolation data. i.e.. $H(\lambda_{i})\mathbf{r}_{i} = \mathbf{w}_{i}$ and $\ell_{j}H(\mu_{j}) = \mathbf{v}_{j}$.

Proof. Multiplying the Eq (3.2) by $s^{\alpha}$ and subtracting it from the Eq (3.2) obtain:

Multiplying Eq (3.4) by $\mathbf{e}_{i}$ on the right and setting $s = \lambda_{i}$ to imply

This proves that the Theorem 3.1is satisfied with the right interpolation data. We can also prove that it is satisfied with the left interpolation data analogously by multiplying the Eq (3.4) with $\mathbf{e}_{j}^{*}$ on the left and setting $s = \mu_{j}$.

When the Loewner pencil is regular, the minimal realization can be constructed according to the Theorem 3.1. However, the Loewner pencil is singular due to the redundant data. In case of SISO systems, the matrices $\mathbf{L} \text{ and } \mathbf{R}^{T}$ are unit vectors. The following assumption is given:

If the assumption the Eq (3.5) is satisfied, for some $x$, a short SVD is computed as following:

Where $\ Y \in \mathbb{C}^{k \times r}$ and ${X} \in \mathbb{C}^{r \times k}$, $\operatorname{rank}(x^{\alpha} \mathbb{L}_{f}-\sigma\mathbb{L}_{f}) = \operatorname{rank}(\Sigma) = \operatorname{size}(\Sigma) = r$.

Theorem 3.2. If the Eq (3.5) is satisfied, and the short SVD exists, a minimal realization $[E, A, B, C]$ is given as follows:

The following proof proposes that above theory satisfies the interpolation data.

Proof. For the system $[E, A, B, C]$, $Y \text{ and } X$ are the generalized controllability and observability matrices respectively. Detailed introduction and proof are given in ^[1] and ^[25]. That is to say:

To prove the realization satisfies $H(\lambda_{i})\mathbf{r}_{i} = \mathbf{w}_{i}, \ell_{j} H(\mu_{j}) = \mathbf{v}_{j}\ (i = 1, \cdots, k)$, the following calculations are carried out:

Therefore we demonstrate that the realization holds for the left and right interpolation conditions.

Remark 3.1. In case of MIMO systems, the assumption described by the Eq (3.5) may not be satisfied. In order to get a minimal realization, the matrix $D \in \mathbb{C}^{p \times m}$ term is considered. Then the shift Loewner matrix becomes $\sigma\mathbb{L}_f - \mathbf{L} D \mathbf{R}$, and $\mathbf{V}, \ \mathbf{W}$ are replaced by $\mathbf{V} - \mathbf{L} D, \ \mathbf{W}-D \mathbf{R}.$ In this way, a suitable $D$ can be found so that the assumption is satisfied. According to the Theorem 3.1 we can construct as follows:

The system [$E, A, B, C, D$] is a realization.

Remark 3.2. The complex conjugate terms will appear in the calculation results of the above algorithm. In order to obtain the real matrix entries, a change needs to be performed by using the matrix ${\Delta }$. Then,

Where

The complex conjugate matrix can be changed to a real matrix by the calculation of the Eq (3.8). And the Loewner matrix $\mathbb{L}_{r}$ and the shifted Loewner matrix $\sigma\mathbb{L}_{r}$ with real entries are also satisfied with the Sylvester Eqs (3.2) and (3.3). A detailed certification process is given in Appendix B of ^[26].

3.2.   Realization of FoS with unknown $\alpha$

In practical applications, the fractional-order $\alpha$ may not be accessible beforehand, only a range for the order $\alpha$ is known. Fortunately, $\alpha$ can only be a parameter satisfying $0 < \alpha < 1$ according to the underlying characterization of FoS. By sampling with allowable precision, one can find the optimal $\alpha$ from the sampled values $[\alpha_1, ....\alpha_L]$. Hereby we propose two criteria to choose the most appropriate $\alpha_l$.

Optimal $\alpha$ selection criterion 1 Minimal Order

As pointed out in ^[1,25], one of the main advantages of the Loewner framework is that the minimal order of the interpolating model can be obtained by evaluating the rank of the Loewner Matrix Pencil. It shows that $r \ (\text{ i.e., rank}(x\mathbb{L} - \sigma\mathbb{L}))$ is the order of constructed system in ^[7].

For fractional-order system, given $\alpha_l$ from $[\alpha_1, ....\alpha_L]$, the minimal order $n(\alpha_l)$ is obtained by the Equation (3.5).

The commensurate order $\alpha$ that minimizes the system highest order $n$ is optimal.

Optimal $\alpha$ selection criterion 2 Minimal Interpolation Error

In order to research the optimal commensurate order-$\alpha$, the interpolation data is used to research the optimal order in a least-squares method.We turn this problem into an optimization problem. Suppose that we have constructed a system state space model $[{E}_{l}, {A}_{l}, {B}_{l}, {C}_{l}]$ with order-$\alpha_{l}$ in terms of the interpolation data Eqs (2.3) and (2.4), the transfer function is:

Here a new set of interpolation data like the Eqs (2.3) and (2.4) is obtained by sampling the transfer function Eq (2.2). The right and left interpolation data:

The data described the Eq (3.10) is used to fit the transfer function in least-squares method. The following minimization problem can be solved. The error can be derived in least squares,

For each given commensurate differential order $\alpha_{l} \in\left(0, 1\right)$, the coefficient matrix of the fractional-order system is identified by the method proposed in the Section 3.1. And the function value of $J(\alpha_{l})$ can be calculated when taking different values of $\alpha_{l}$ separately according to the Eq (3.11). Then we can look for the $\alpha_{l}$ that minimizes the $J(\alpha_{l})$ as an estimate of the fractional differential order.

4.   Rational approximation of FoTD systems

In this section, the solution of the problem 2 is proposed. Likewise, we apply the Loewner framework to find a linear integer-order model to approximate FoTD systems. The frequency response pairs $S_i = G(s_i)$ are divided into the right and left interpolation data. To simplify the exposition, let the number of input-output pairs be even. i.e., $N = 2 k$. Thus the interpolation data is given as follows according to the Eqs (2.9) and (2.10) with $\mathbf{r}_{i} = \ell_{j} = 1, \ (i, j = 1, \cdots, k)$.

the right interpolation data

the left interpolation data

where $\mathbf{L} \text{ and } \mathbf{R}^{T}$ are unit vectors, $\boldsymbol G(\lambda_{i}) = \mathbf{w}_{i} \text{ and } \boldsymbol G(\mu_{j}) = \mathbf{v}_{j}$. Therefore the Loewner matrix $\mathbb{L}$ and shifted Loewner matrix $\mathbb{L}_{\sigma}$ are computed according to the Eq (2.8).

Then a linear integer-order model with $r$-order can be constructed in terms of Theorem 2.1 or Theorem 2.2. The detailed proof is given in ^[1].

5.   Examples

5.1.   Example 1

Consider a simple fractional system is described by the following fractional differential equations.

where the commensurate order $\alpha = 0.5$, thus the associated transfer function is

We chose the right input frequencies $\Lambda = \text{diag} (2j, -2j, 4j, -4j)$, and the left input frequencies $\boldsymbol {M} = \text{diag} (j, -j, 3j, -3j)$.

Since the system $\Sigma$ is a single-input single-output system, we set the the right and left input directions as $\mathbf{R} = \mathbf{L}^{T} = \left[1 \quad 1 \quad 1 \quad 1\right]$, then the right and left responses $\mathbf{V} \text{ and } \mathbf{W}$ can be computed by sampling $G(s)$ at the right and left frequencies.

To generate state-space realization with real entries, all the above complex matrices are transformed into real matrices according to the Eq (3.8). The real matrices are computed as follows with $\hat{\Delta } = \frac{1}{\sqrt{2}}\ {\mathbf{blkdiag}} \left(\left[\begin{array}{rr}{1} & {-j} \\ {1} & {j}\end{array}\right], \left[\begin{array}{rr}{1} & {-j} \\ {1} & {j}\end{array}\right] \ \right)$.

with $\mathbf{R}_{r} = \mathbf{L}_{r}^{T} = \left(1393/985 \quad 0 \quad 1393/985 \quad 0 \right)$ and the Loewner pencil $(\mathbb{L}_{r}, \; \sigma\mathbb{L}_{r})$ as:

We compute the rank of Loewner matrix and shifted Loewner matrix is equal to 2, so the Loewner pencil $(\mathbb{L}_{r}, \; \sigma\mathbb{L}_{r})$ is not regular.

While ${D} = 0$, according to Theorem 3.2, we check the assumption Eq (3.5) with $x \in\left\{2, -2, 4, -4\right\} \cup\left\{ 1, -1, 3, -3 \right\}$.

we compute the short $SVD$ of $\left(\mathbb{L}_{r} - \sigma\mathbb{L}_{r}\right)$ Then,

Thus, we can construct the minimal realization of the fractional-order system.

and the transfer function is :

5.2.   Example 2

In this example, the optimal commensurate order-$\alpha$ will be identified by the least squares method. The system and the interpolation data is the same as that in the Example 1. We take $\alpha_{l} = [0.1, \ 0.2, \ \cdots, \ 0.9]$. For each given the order $\alpha_{l}$, the system matrices can be constructed. Another right and left input frequencies $\Lambda = (6j, -6j, 8j, -8j)$, and $\boldsymbol {M} = (5j, -5j, 7j, -7j)$ are chosen. Then the function value of the Eq (3.11) can be calculated and shown as the Table 1. It shows that the $J(\alpha_{l})$ is minimum while $\alpha_{l} = 0.5$. That is to say the commensurate order of interpolation model is the same as that of the original system $\Sigma$.

On the other hand, the highest order of the system model constructed can be observe as the Table 1. We can find the order is minimum only while $\alpha_{l} = 0.5$ via comparing the highest orders under different commensurate order $\alpha_{l}$.

Then we increased the data with $k = q = 50$, the error function $J(\alpha_{l})$ and highest order of the system model under $\alpha_{l} = 0.1, 0.2, \cdots, 1$ are shown as Table 2. It indicates that the highest order increases with the increase in the amount of data if the commensurate order $\alpha$ is not equal to 0.5, a optimal commensurate order that minimizes the highest order and the error function $J(\alpha_{l})$.

5.3.   Example 3

Consider a linear fractional-order time delay system whose transfer function is given as

and consider the following two disjoint sets of input frequencies $\boldsymbol \Lambda = \text{diag} (0.2, 0.4, 0.6, 0.8)$, and $\boldsymbol {M} = \text{diag} (0.1, 0.3, 0.5, 0.7)$. And the corresponding frequency responses of the transfer function are:

Then the Loewner matrix and shifted Loewner matrix can be calculated as:

Loewner Matrix Pencil is full rank, a integer order model with 4-th order can be constructed according to the Theorem 2.1. The transfer function of approximation model is

The step responses of the original system and approximate system are shown in Figure 1, the red solid line represents original system. It shows that this approximate model works well for step responses.

Figure 2 is the bode diagram of the approximate system and the original system. It shows that the approximate model is slightly better at low frequencies and worse at high frequencies. Moreover the phase advances the original system $2\pi$. We can solve this problem by performing phase correction on the approximate system.

6.   Conclusions and future work

We present a generalization of the LMM for realization of commensurate fractional-order model. The concept and definition of the generalized fractional-order Loewner Matrix are proposed, which enriches the Loewner framework. More importantly, generalized Loewner framework derives the most simplest model (in the sense of McMillan order) while the commensurate fractional power $\alpha$ coincides the true value, which provides an insightful interpretation on the minimal realization both in fractional-order setting and in integer-order setting. An illustrative example is included for demonstration. For fractional-order time delay systems, time and frequency responses of the approximate model based on the Loewner framework can match the original system.

In this paper, commensurate fractional-order model is considered, i.e., the fractional powers in the model are integer multiples of a single real number $\alpha$. Further adopting LMM for building non-commensurate models and fractional time delay models would be worth exploring in the future. It should be pointed out that extending the Loewner framework to non-commensurate models is not straight forward as the formulation of Loewner Matrix becomes more complicated.

Acknowledgments

This work was supported in part by National Natural Science Foundation of China under Grant (No. 61873336, 61873335), in part by the National Key Research and Development Program (No. 2020YFB1708200), in part by the Foreign Expert Program (No. 20WZ2501100) granted by the Shanghai Science and Technology Commission of Shanghai Municipality (Shanghai Administration of Foreign Experts Affairs), in part by 111 Project (No. D18003) granted by the State Administration of Foreign Experts Affairs, in part by the Fundamental Research Funds for the Central Universities under Grant (No. FRF-BD-19-002A).

Conflict of interest

All authors declare no conflicts of interest in this paper.