Next Article in Journal
Modified Classical Graph Algorithms for the DNA Fragment Assembly Problem
Previous Article in Journal
A Comparative Study of Modern Heuristics on the School Timetabling Problem
Article

A CS Recovery Algorithm for Model and Time Delay Identification of MISO-FIR Systems

by 1,2,* and 2
1
Key Laboratory of Advanced Process Control for Light Industry (Ministry of Education), Jiangnan University, Wuxi 214122, China
2
School of Internet of Things Engineering, Jiangnan University, Wuxi 214122, China
*
Author to whom correspondence should be addressed.
Academic Editor: Paul M. Goggans
Algorithms 2015, 8(3), 743-753; https://0-doi-org.brum.beds.ac.uk/10.3390/a8030743
Received: 29 June 2015 / Revised: 6 August 2015 / Accepted: 2 September 2015 / Published: 10 September 2015

Abstract

This paper considers identifying the multiple input single output finite impulse response (MISO-FIR) systems with unknown time delays and orders. Generally, parameters, orders and time delays of an MISO system are separately identified from different algorithms. In this paper, we aim to perform the model identification and time delay estimation simultaneously from a limited number of observations. For an MISO-FIR system with many inputs and unknown input time delays, the corresponding identification model contains a large number of parameters, requiring a great number of observations for identification and leading to a heavy computational burden. Inspired by the compressed sensing (CS) recovery theory, a threshold orthogonal matching pursuit algorithm (TH-OMP) is presented to simultaneously identify the parameters, the orders and the time delays of the MISO-FIR systems. The proposed algorithm requires only a small number of sampled data compared to the conventional identification methods, such as the least squares method. The effectiveness of the proposed algorithm is verified by simulation results.
Keywords: compressed sensing; sparse; parameter identification; gradient projection pursuit algorithm; time delay estimation compressed sensing; sparse; parameter identification; gradient projection pursuit algorithm; time delay estimation

1. Introduction

In many industry processes, time delay is an unavoidable behavior in their dynamical models. For example, in a chemical industry, analyzers need a certain amount of time to process an analysis, which is used as a part of the control loop [1]. Analyzing and designing a controller for a system with time delays depend on an effective dynamical system model. System identification is an effective way to build mathematical models of dynamical systems from observed input-output data. Conventional identification methods, such as the least squares methods [2], the stochastic gradient methods [3] and the maximum likelihood estimation methods [4], usually require a large number of sampled data and that the model structures, including time delays, are a priori known [5,6]. However, this is not the case in many situations in practice for which the time delays and systems orders have to be estimated together with the parameters. Most of the process industries, such as distillation columns, heat exchangers and reactors, can be modeled as multiple input multiple output (MIMO) or multiple input single output (MISO) systems with long time delays [1]. For a multiple input single output finite impulse response (MISO-FIR) system with many inputs and unknown input time delays, the parameterized model contains a large number of parameters, requiring a great number of observations and a heavy computational burden for identification. Furthermore, the measuring cost is an issue. Especially for chemical processes, it will take a long time to get enough data because of the long analysis time. Therefore, it is a challenging work to find new methods to identify the time delays and parameters of multivariable systems with less computational burden using a small number of sampled data.
System models of practical interest, especially the control-oriented models, often require low order, which implies that the systems contain only a few non-zero parameters. Because the time delays and orders are unknown, we can parameterize such an MISO-FIR system by a high-dimensional parameter vector, with only a few non-zero elements. These systems can be termed as sparse systems. On the other hand, a sparse MISO-FIR system contains a long impulse response, but only a few terms are non-zeros. The sparse systems are widely discussed in many fields [7,8,9,10,11]. For such sparse systems, the identification objective is to estimate the sparse parameter vector and to read the time delays and orders from the structure of the estimated parameter vector.
In this paper, we consider the identification of the sparse MISO-FIR systems by using the compressed sensing (CS) recovery theory [12,13,14,15]. CS has been widely applied in signal processing and communication systems, which enables the recovery of an unknown vector from a set of measurements under the assumption that the signal is sparse and certain conditions on the measurement matrices [16,17,18,19]. The identification of a sparse system can be viewed as a CS recovery problem. The sparse recovery problem can be solved by the l 1 -norm convex relaxation. However, the l 1 regularization schemes are always complex, requiring heavy computational burdens [20,21,22]. The greedy methods have speed advantages and are easily implemented. In this literature, a block orthogonal matching pursuit (BOMP) algorithm has been applied to identify the determined linear time invariant (LTI) and linear time variant (LTV) auto-regressive with external (ARX) models without disturbances [9]. In this paper, we will consider identifying the parameters, orders and time delays of the MISO-FIR systems from a small number of noisy measured data.
Briefly, the structure of this paper is as follows. Section 2 gives the model description of MISO-FIR systems and describes the identification problem. Section 3 presents a threshold orthogonal matching pursuit (TH-OMP) algorithm based on the compressed sensing theory. Section 4 provides a simulation example to show the effectiveness of the proposed algorithm. Finally, some concluding remarks are given in Section 5.

2. Problem Description

Consider an MISO-FIR system with time delays:
y ( t ) = i = 1 r z d i B i ( z ) u i ( t ) + v ( t )
where u i ( t ) is the i-th input, y ( t ) the output, d i the time delay of the i-th input channel and v ( t ) a white noise with zero mean and variance σ 2 , B i ( z ) , i = 1 , 2 , , r are the polynomials in the unit backward shift operator z 1 (i.e., z 1 y ( t ) = y ( t 1 ) ) and:
B i ( z ) : = b i 1 z 1 + b i 2 z 2 + + b i n b i z n b i
Assume that y ( t ) = 0 , u ( t ) = 0 and v ( t ) = 0 for t 0 and the orders n b i , the time delays d i , as well as the parameters b i j are unknown. Thus, the identification goal is to estimate the unknown parameters b i j , orders n b i and the time delays d i from the observations.
Define the information vector:
φ ( t ) : = [ u 1 ( t 1 ) , , u 1 ( t d 1 ) , u 1 ( t d 1 1 ) , , u 1 ( t d 1 n b 1 ) , , u 1 ( t l ) , , u r ( t 1 ) , , u r ( t d r ) , u r ( t d r 1 ) , , u r ( t d r n b r ) , , u r ( t l ) ] T R N
where l is the input maximum regression length. l must be chosen to capture all of the time delays and the degrees of the polynomials B i ( z ) . Thus, l max ( d i + n b i ) . The dimension of φ ( t ) is N = l r . The corresponding parameter vector θ is:
θ : = [ 0 , , 0 , b 11 , , b 1 n b 1 , 0 , , 0 , b i 1 , , b i n b i , 0 , , 0 , b r 1 , , b r n b r , 0 , , 0 ] T R N
Then, Equation (1) can be written as:
y ( t ) = φ T ( t ) θ + v ( t )
Because of the unknown orders and time delays, the parameter vector θ contains many zeros, but only a few non-zeros. According to the CS theory, the parameter vector θ can be viewed as a sparse signal, and the number of non-zero entries K = i = 1 r n b i is the sparsity level of θ [8]. Let θ 0 denote the number of non-zero entries in θ. Then, the identification problem can be described as:
θ ^ = arg min θ 0 , s . t . y φ T ( t ) θ 2 < ε ,
where θ ^ is the estimate of θ and ε > 0 is the error tolerance.
For m observations, define the stacked matrix and vectors as:
Φ m : = φ T ( 1 ) φ T ( 2 ) φ T ( m ) R m × N , y m : = y ( 1 ) y ( 2 ) y ( m ) R m , v m : = v ( 1 ) v ( 2 ) v ( m ) R m
Then, we have:
y m = Φ m θ + v m
If there are enough measurements, i.e., m N , according to the least squares principle, we can get the least squares estimate of θ,
θ ^ LS = ( Φ m T Φ m ) 1 Φ m T y m
However, from Equation (2), it is obvious that N is a large number, especially for systems with a large number of inputs. Therefore, it will take a lot of time and effort to get enough measurements. In order to improve the identification efficiency, in this paper, we aim to identify the parameter vector θ using less observations ( K < m < N ) based on the CS recovery theory.

3. Identification Algorithm

The CS theory, which was introduced by Cand e ` s, Romberg and Tao [12,13] and Donoho [14], has been widely used in signal processing. It indicates that an unknown signal vector can be recovered from a small set of measurements under the assumptions that the signal vector is sparse and some conditions on the measurement matrix. The identification problem can be viewed as a CS recovery problem in that the parameter vector θ is sparse.
From Equation (4), we can see that the output vector y m can be written as a linear combination of the columns of Φ m plus the noise vector v m , i.e.,
y m = ϕ 1 θ 1 + ϕ 2 θ 2 + + ϕ i θ i + + ϕ N θ N + v m
where θ i is the i-th element of θ and ϕ i is the i-th column of Φ m . Because there are only K non-zero parameters in θ, the main idea is to find the K non-zero items at the right-hand side of Equation (5).
In this paper, we propose a TH-OMP algorithm to do the identification. First, we give some notations. Let k = 1 , 2 , be the iterative number and λ k be the solution support of the k-th iterative number; Λ k is a set composed of λ i , i = 1 , 2 , , k ; r k denotes the residual of the k-th iterative number; Φ Λ k is the sub-information matrix composed by the k columns of Φ m indexed by Λ k . θ ^ k is the parameter estimation of the k-th iterative number. The TH-OMP algorithm is initialized as r 0 = y m , Λ 0 = and θ ^ Λ 0 = 0 . At the k-th iteration, minimizing the criteria function:
ε ( i ) = min r k 1 ϕ i θ i 2 , i = 1 , 2 , , N
with respect to θ i yields:
θ i = ϕ i T r k 1 θ i 2
Plugging it back into ε ( i ) , we have:
ε ( i ) = min θ i ϕ i θ i r k 1 2 = ϕ i T r k 1 ϕ i 2 ϕ i r k 1 2 = r k 1 2 2 ( ϕ i T r k 1 ) 2 ϕ i 2 + ( ϕ i T r k 1 ) 2 ϕ i 2 = r k 1 2 ( ϕ i T r k 1 ) 2 ϕ i 2 = r k 1 2 ( 1 ϕ i ϕ i T r k 1 ) 2
From the above equation, we can see that the quest for the smallest error is equivalent to the quest for the largest inner product between the residual r k 1 and the normalized column vectors of Φ m . Thus, the k-th solution support can be obtained by:
λ k = arg max i = 1 , 2 , , N | r k 1 , ϕ i ϕ i |
Update the support set Λ k and the sub-information matrix Φ Λ k by:
Λ k : = Λ k 1 λ k , Φ Λ k : = Φ Λ k 1 ϕ λ k
The next step is to find the k-th provisional parameter estimate from the sub-information matrix Φ Λ k and measurement vector y m . Minimizing the cost function:
J ( θ Λ k ) = y m Φ Λ k θ Λ k 2
leads to the k-th parameter estimate:
θ ^ Λ k = ( Φ Λ k T Φ Λ k ) 1 Φ Λ k T y m
Then, the k-th residual can be computed by:
r k = y m Φ Λ k θ ^ Λ k
If the system is noise free, then the residual will be zero after K iterations, and the non-zero parameters, as well as their positions can be exactly determined. This is the so-called orthogonal matching pursuit (OMP) algorithm. However, for systems disturbed by noises, the OMP algorithm cannot give the sparsest solution, because the residual r K is non-zero after K iterations, and the iteration continues. In order to solve this problem, an appropriate small threshold ϵ can be set to filter the parameter estimate θ ^ Λ k . If | θ ^ Λ ϵ | < ϵ , let θ ^ Λ ϵ = 0 , where θ ^ Λ ϵ is the i-th element of θ ^ Λ k . The choice of ϵ will be discussed in the simulation part. Denote θ ^ Λ k ϵ as the parameter estimate after the filtering. Then, the residual can be computed by:
r k = y m Φ Λ k θ ^ Λ k ϵ
if r k < ε , then stop the iteration; otherwise, go to apply another iteration. If the iteration stops at k, we can recover the parameter estimate θ ^ R N from θ ^ Λ k ϵ R k based on the support set Λ k .
It is worth noting that taking the derivative of J ( θ Λ k ) with respect to θ ^ Λ k gives:
J ( θ ^ Λ k ) θ ^ Λ k = Φ Λ k T ( y m Φ Λ k θ ^ Λ k ) = Φ Λ k T r k = 0
The above equation shows that the residual r k is orthogonal to the sub-information matrix Φ Λ k . Therefore, it does not require computing the inner products of r k and the columns in Φ Λ k during the ( k + 1 ) -th iteration. Modify Equation (7) as:
λ k = arg max i Ω \ Λ k 1 | r k 1 , ϕ i ϕ i | , Ω = { 1 , 2 , , N }
Then, the computational burden can be reduced.
From the recovered parameter estimation θ ^ , the orders and time delays of each input channel can be estimated. The orders n b i can be determined by the number of elements of each non-zero block. Then, the time delay estimates can be computed from the orders n b i and the input maximum regression length l. It can be seen from Equation (2) that θ contains ( r + 1 ) zero blocks. Denote the number of zeros of each zero-block as n i , i = 1 , 2 , , r + 1 . Then, the time delay estimates can be computed by:
d ^ 1 = n 1 , d ^ i = n i ( l d ^ i 1 n b i ) , i = 2 , 3 , , r
The proposed TH-OMP algorithm can be implemented as the following steps.
  • Set the initial values: let k = 0 , r 0 = y m , Λ 0 = and θ ^ Λ 0 = 0 . Give the error tolerant ε.
  • Increment k by one; compute the solution support λ k using Equation (14).
  • Update the support set Λ k and the sub-information matrix Φ Λ k using Equation (8).
  • Compute the parameter estimate θ ^ Λ k by using Equation (10).
  • Using the threshold ϵ to filter θ ^ Λ k to get θ ^ Λ k ϵ .
  • Compute the residual using Equation (12). If r k < ε , stop the iteration; otherwise, go back to Step 2 to apply another iteration.
  • Recover the parameter estimate θ ^ R N from θ ^ Λ k ϵ R k based on the support set Λ k .
  • Read the orders n b i from θ ^ R N and estimate the time delays using Equation (15).

4. Simulation Example

Consider the steam water heat exchanger system in Figure 1, where the steam flow u 1 and the water flow u 2 are two manipulated variables to control the water temperature y. The discrete input-output representation of the process is expressed as the following equation.
Figure 1. The steam water heat exchanger.
Figure 1. The steam water heat exchanger.
Algorithms 08 00743 g001
y ( t ) = z d 1 B 1 ( z ) u 1 ( t ) + z d 2 B 2 ( z ) u 2 ( t ) + v ( t )
B 1 ( z ) = 0 . 95 z 1 + 0 . 72 z 2 + 0 . 40 z 3 ,
B 2 ( z ) = 0 . 63 z 1 0 . 47 z 2 0 . 25 z 3 ,
d 1 = 42 , d 2 = 20
Let l = 50 ; then the parameter vector to be identified is:
θ = [ 0 42 , 0 . 95 , 0 . 72 , 0 . 40 , 0 25 , 0 . 63 , 0 . 47 , 0 . 25 , 0 27 ] T R N
It can be seen from the above equation that the dimension of the parameter vector θ is N = 100 , but only six of the parameters are non-zeros.
In the identification, the inputs { u i ( t ) } , i = 1 , 2 are taken as independent persistent excitation signal sequences with zero mean and unit variances and { v ( t ) } as white noise sequences with zero mean and variance σ 2 = 0 . 10 2 . Let ϵ = 0 . 05 , applying the proposed TH-OMP algorithm and the least squares (LS) algorithm to estimate the system with m measurements. The parameter estimation errors δ : = ( θ θ ^ ) / θ versus m are shown in Figure 2.
From Figure 2, we can see that under a limited dataset ( m < 100 ) , the parameter estimation error of the TH-OMP algorithm approaches zero, and the estimation accuracy is much higher than that of the LS algorithm.
Let m = 80 ; using the TH-OMP algorithm to estimate θ, the parameter estimation errors δ versus iteration k are shown in Figure 3; the non-zero parameter estimates versus iteration k are shown in Table 1 and Figure 4, where the dash dot lines are true parameters and the solid lines are parameter estimates.
Figure 2. The estimation error δ versus m measurements.
Figure 2. The estimation error δ versus m measurements.
Algorithms 08 00743 g002
Figure 3. The estimation error δ versus iteration k (m = 80, σ 2 = 0 . 10 2 ).
Figure 3. The estimation error δ versus iteration k (m = 80, σ 2 = 0 . 10 2 ).
Algorithms 08 00743 g003
Figure 4. The parameters estimates versus iteration k.
Figure 4. The parameters estimates versus iteration k.
Algorithms 08 00743 g004
Table 1. The estimation error δ versus iteration k (m=80, σ 2 = 0 . 10 2 ).
Table 1. The estimation error δ versus iteration k (m=80, σ 2 = 0 . 10 2 ).
k b 11 b 12 b 13 b 21 b 22 b 23 δ ( % )
11.05670.00000.00000.00000.00000.000077.8404
21.00560.73570.00000.00000.00000.000061.0813
31.02740.73490.0000−0.69820.00000.000044.8210
40.94220.76350.0000−0.7006−0.46040.000031.8601
50.97480.74290.3734−0.6551−0.46050.000016.9639
60.94360.70410.3936−0.6397−0.4727−0.25221.3977
70.94950.70710.3973−0.6408−0.4715−0.25371.1650
80.94660.71270.3986−0.6422−0.4707−0.25090.9755
90.94870.71200.4022−0.6386−0.4722−0.24900.8137
100.94580.71290.3988−0.6380−0.4766−0.24980.8841
True values0.95000.72000.4000−0.6300−0.4700−0.2500
After 10 iterations, the estimation error is δ = 0 . 8841 % , and the estimated parameter vector is:
θ ^ ϵ = 0 . 05 = [ 0 42 , 0 . 9472 , 0 . 7105 , 0 . 3921 , 0 25 , 0 . 6246 , 0 . 4751 , 0 . 2524 , 0 27 ] T R N
If we choose ϵ = 0 . 01 , then we can get the estimated parameter vector as:
θ ^ ϵ = 0 . 01 = [ 0 4 , 0 . 0217 , 0 25 , 0 . 0303 , 0 11 , 0 . 9458 , 0 . 7129 , 0 . 3988 , 0 25 0 . 6220 , 0 . 4634 , 0 . 2502 , 0 13 , 0 . 0244 , 0 13 ] T R N
It is obvious that the parameter vector θ has only three zero blocks in this example. However, we can see from the above equation that there exist three undesirable parameter estimates θ ^ 5 , θ ^ 31 and θ ^ 87 . Moreover, the three parameter estimates are much smaller than other parameters. Therefore, according to the structure of θ, the parameters θ ^ 5 , θ ^ 31 and θ ^ 87 can be set to zeros. Then, the estimation result is the same as the case when ϵ = 0 . 05 . This implies that even if a too small ϵ is chosen, we can still obtain the effective parameter estimates based on the model structures.
From Equations (15) and (16), we can obtain the orders n b 1 = n b 2 = 3 and get the time delay estimates as:
d ^ 1 = 42 , d ^ 2 = n 2 ( l d 1 n b 1 ) = 25 ( 50 42 3 ) = 20
From Figure 2, Figure 3 and Figure 4, Table 1 and Equations (16) and (17), we can conclude that the TH-OMP algorithm converges faster than the LS algorithm, can achieve a high estimation accuracy by K iterations with finite measurements ( m < N ) and can effectively estimate the time delay of each channel.

5. Conclusions

This paper presents a TH-OMP algorithm for MISO-FIR systems with unknown orders and time delays. The parameter estimates, orders and time delays can be simultaneously estimated from a small number of observation data. The proposed method can be simply implemented and can reduce the measuring cost, as well as improve the identification efficiency. The simulation results are given to demonstrate the performance of the proposed method.

Acknowledgments

This paper is supported by the National Natural Science Foundation of China (Nos. 61203111, 61304138) and the Natural Science Foundation of Jiangsu Province (China; BK20130163).

Author Contributions

Yanjun Liu worked on the main details of the algorithm and was responsible for writing the manuscript. Taiyang Tao helped with the MATLAB simulation.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Normey-Rico, J.E.; Camacho, E.F. Control of Dead-Time Processes; Spinger: Berlin, Germany, 2007. [Google Scholar]
  2. Liu, Y.J.; Ding, F. Convergence properties of the least squares estimation algorithm for multivariable systems. Appl. Math. Model. 2013, 37, 476–483. [Google Scholar] [CrossRef]
  3. Liu, Y.J.; Ding, R. Consistency of the extended gradient identification algorithm for multi-input multi-output systems with moving average noises. Int. J. Comput. Math. 2013, 90, 1840–1852. [Google Scholar] [CrossRef]
  4. Li, J.H.; Ding, F.; Hua, L. Maximum likelihood newton recursive and the newton iterative estimation algorithms for hammerstein CARAR systems. Nonlinear Dyn. 2014, 75, 235–245. [Google Scholar] [CrossRef]
  5. Ljung, L. System Identification, Theory for the User, 2nd ed.; Prentince Hall: Englewood Cliffs, NJ, USA, 1999. [Google Scholar]
  6. Ding, F. System Identification-New Theory and Methods; Science Press: Beijing, China, 2013. [Google Scholar]
  7. Haupt, J.; Bajwa, W.U.; Raz, G.; Nowak, R. Toeplitz compressed sensing matrices with applications to sparse channel estimation. IEEE Trans. Inf. Theory 2010, 56, 5862–5875. [Google Scholar] [CrossRef]
  8. Sanandaji, B.M.; Vincent, T.L.; Wakin, M.B. Exact topology identification of large-scale interconnected dynamical systems from compressive observations. In Proceedings of the 2011 American Control Conference, San Francisco, CA, USA, 29 June–1 July 2011; pp. 649–656.
  9. Sanandaji, B.M.; Vincent, T.L.; Wakin, M.B.; Toth, R. Compressive system identification of LTI and LTV ARX models. In Proceedings of the 50th IEEE Conference on Decision and Control and European Control Conference, Orlando, FL, USA, 12–15 December 2011; pp. 791–798.
  10. Gui, G.; Fumiyuki, A. Stable adaptive sparse filtering algorithms for estimating multiple-input multiple-output channels. IET Commun. 2014, 8, 1032–1040. [Google Scholar] [CrossRef]
  11. Gui, G.; Xu, L.; Fumiyuki, A. Normalized least mean square-based adaptive sparse filtering algorithms for estimating multiple-input multiple-output channels. Wirel. Commun. Mob. Comput. 2015, 15, 1079–1088. [Google Scholar] [CrossRef]
  12. Candès, E.J.; Romberg, J.; Tao, T. Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information. IEEE Trans. Inf. Theory 2006, 52, 489–509. [Google Scholar] [CrossRef]
  13. Candès, E.J. Compressive sampling. In Proceedings of the 2006 International Congress of Mathematics; The European Mathematical Society: Madrid, Spain, 2006; pp. 1433–1452. [Google Scholar]
  14. Donoho, D.L. Compressed sensing. IEEE Trans. Inf. Theory 2006, 52, 1289–1306. [Google Scholar] [CrossRef]
  15. Elad, M. Optimized projections for compressed sensing. IEEE Trans. Signal Process. 2007, 55, 5695–5702. [Google Scholar] [CrossRef]
  16. Baraniuk, R. A lecture on compressive sensing. IEEE Signal Process. Mag. 2007, 24, 118–121. [Google Scholar] [CrossRef]
  17. Candès, E.J.; Wakin, M.B. An introduction to compressive sampling. IEEE Signal Process. Mag. 2008, 25, 21–30. [Google Scholar] [CrossRef]
  18. Tropp, A.J. Just relax: Convex programming methods for ifentifying sparse signals in noise. IEEE Trans. Inf. Theory 2006, 52, 1030–1051. [Google Scholar] [CrossRef]
  19. Candès, E.J. Robust uncertainty principles: Extra signal reconstruction from highly frequency information. IEEE Trans. Inf. Theory 2006, 52, 489–509. [Google Scholar] [CrossRef]
  20. Tóth, R.; Sanandaji, B.M.; Poolla, K.; Vincent, T.L. Compressive system identification in the linear time-invariant framework. In Proceedings of the 50th IEEE Conference on Decision and Control and European Control Conference, Orlando, FL, USA, 12–15 December 2011; pp. 783–790.
  21. Yang, Z.; Zhang, C.S.; Deng, J.; Lu, W.M. Orthonormal expansion 1 minimization algorithms for compressed sensing. IEEE Trans. Signal Process. 2011, 59, 6285–6290. [Google Scholar] [CrossRef]
  22. Le, V.L.; Lauer, F.; Bloch, G. Selective 1 minimization for sparse recovery. IEEE Trans. Autom. Control 2014, 59, 3008–3013. [Google Scholar] [CrossRef]
Back to TopTop