Performance Improvement of Discrete-time Linear Control Systems Subject to Varying Sampling Rates Using the Tikhonov Regularization Method

Fernando Agustín Pazos Anibal Zanini Amit Bhaya

Fernando Agustín Pazos, Anibal Zanini and Amit Bhaya. Performance Improvement of Discrete-time Linear Control Systems Subject to Varying Sampling Rates Using the Tikhonov Regularization Method. International Journal of Automation and Computing. doi: 10.1007/s11633-019-1205-8
Citation: Fernando Agustín Pazos, Anibal Zanini and Amit Bhaya. Performance Improvement of Discrete-time Linear Control Systems Subject to Varying Sampling Rates Using the Tikhonov Regularization Method. International Journal of Automation and Computing. doi: 10.1007/s11633-019-1205-8

doi: 10.1007/s11633-019-1205-8

Performance Improvement of Discrete-time Linear Control Systems Subject to Varying Sampling Rates Using the Tikhonov Regularization Method

More Information
    Author Bio:

    Fernando Agustín Pazos received the B. Sc. degree in electronic engineering from the University of Buenos Aires, Argentina in 1990, the M. Sc. and D. Sc. degrees in electrical engineering from the Federal University of Rio de Janeiro, Brazil in 2000 and 2007, respectively. He is a professor in the Electronic and Telecommunications Engineering Department, the State University of Rio de Janeiro, Brazil. He is an author of “Automação de Sistemas & Robótica” (Axcel Books, Rio de Janeiro, Brazil 2002).His research interest is development of algorithms to solve several optimization problems interpreting them as closed loop control systems. E-mail: quini.coppe@gmail.com (Corresponding author) ORCID iD: 0000-0001-9951-9790

    Anibal Zanini received the B. Eng. degree in electrical engineering from Universidad Nacional de Rosario, Argentina in 1977, and recived the Ph. D. degree in electrical engineering at the Escuela Técnica Superior de Ingenieros Industriales of the Universidad Politécnica de Madri, Spain in 1983. He has worked for 20 years as project leader and senior researcher in the automation area at metallurgical industry. Since 1998, his main activity is at Universidad de Buenos Aires. He is a professor at the Engineering School of “Universidad de Buenos Aires”. His research interests include industrial control, identification and adaptive control. E-mail: anibal.zanini@gmail.com ORCID iD: 0000-0001-9641-562X

    Amit Bhaya received the B. Eng. degree in eletrical engineering from the Indian Institute of Technology, India in 1981, and received the Ph. D. degree in eletrical engineering from the University of California at Berkeley, USA in 1986. He is a professor in the Electrical Engineering Department, Graduate School of Engineering, the Federal University of Rio de Janeiro (COPPE/UFRJ), Brazil. He coauthored the books “Matrix Diagonal Stability in Systems and Computation” (Birkhäuser, 2000), and “Control Perspectives on Numerical Algorithms and Matrix Problems” (SIAM, 2006). He was an associate editor of the IEEE Transactions on Neural Networks and Learning Systems from 2010 to 2012. His research interests include systems and control theory, parallel computation, neural networks, matrix stability theory, mathematical ecology and optimization of economic dynamic systems.E-mail: amit@nacad.ufrj.brORCID iD: 0000-0002-3144-1242

    Corresponding author: email: quini.coppe@gmail.com (corresponding author), ORCID: 0000-0001-9951-9790
图(7) / 表(2)
计量
  • 文章访问数:  465
  • HTML全文浏览量:  323
  • PDF下载量:  7
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-07-16
  • 录用日期:  2019-10-08
  • 网络出版日期:  2019-12-23

Performance Improvement of Discrete-time Linear Control Systems Subject to Varying Sampling Rates Using the Tikhonov Regularization Method

doi: 10.1007/s11633-019-1205-8
    作者简介:

    Fernando Agustín Pazos received the B. Sc. degree in electronic engineering from the University of Buenos Aires, Argentina in 1990, the M. Sc. and D. Sc. degrees in electrical engineering from the Federal University of Rio de Janeiro, Brazil in 2000 and 2007, respectively. He is a professor in the Electronic and Telecommunications Engineering Department, the State University of Rio de Janeiro, Brazil. He is an author of “Automação de Sistemas & Robótica” (Axcel Books, Rio de Janeiro, Brazil 2002).His research interest is development of algorithms to solve several optimization problems interpreting them as closed loop control systems. E-mail: quini.coppe@gmail.com (Corresponding author) ORCID iD: 0000-0001-9951-9790

    Anibal Zanini received the B. Eng. degree in electrical engineering from Universidad Nacional de Rosario, Argentina in 1977, and recived the Ph. D. degree in electrical engineering at the Escuela Técnica Superior de Ingenieros Industriales of the Universidad Politécnica de Madri, Spain in 1983. He has worked for 20 years as project leader and senior researcher in the automation area at metallurgical industry. Since 1998, his main activity is at Universidad de Buenos Aires. He is a professor at the Engineering School of “Universidad de Buenos Aires”. His research interests include industrial control, identification and adaptive control. E-mail: anibal.zanini@gmail.com ORCID iD: 0000-0001-9641-562X

    Amit Bhaya received the B. Eng. degree in eletrical engineering from the Indian Institute of Technology, India in 1981, and received the Ph. D. degree in eletrical engineering from the University of California at Berkeley, USA in 1986. He is a professor in the Electrical Engineering Department, Graduate School of Engineering, the Federal University of Rio de Janeiro (COPPE/UFRJ), Brazil. He coauthored the books “Matrix Diagonal Stability in Systems and Computation” (Birkhäuser, 2000), and “Control Perspectives on Numerical Algorithms and Matrix Problems” (SIAM, 2006). He was an associate editor of the IEEE Transactions on Neural Networks and Learning Systems from 2010 to 2012. His research interests include systems and control theory, parallel computation, neural networks, matrix stability theory, mathematical ecology and optimization of economic dynamic systems.E-mail: amit@nacad.ufrj.brORCID iD: 0000-0002-3144-1242

    通讯作者: email: quini.coppe@gmail.com (corresponding author), ORCID: 0000-0001-9951-9790

English Abstract

Fernando Agustín Pazos, Anibal Zanini and Amit Bhaya. Performance Improvement of Discrete-time Linear Control Systems Subject to Varying Sampling Rates Using the Tikhonov Regularization Method. International Journal of Automation and Computing. doi: 10.1007/s11633-019-1205-8
Citation: Fernando Agustín Pazos, Anibal Zanini and Amit Bhaya. Performance Improvement of Discrete-time Linear Control Systems Subject to Varying Sampling Rates Using the Tikhonov Regularization Method. International Journal of Automation and Computing. doi: 10.1007/s11633-019-1205-8
    • Discrete-time control systems are usually subject to several kinds of disturbances in practical applications, such as inaccuracies, measurements errors or timing problems. With respect to timing problems, ideally all the system components (plant, sensors, actuators, controllers and interfaces) are synchronous, but in reality sensing delays, excessive computational load, communication errors or overload in the communication channels occur frequently and are responsible for the differences in the moments in which each device acts[1]. In many practical applications, small time-variations or asynchronicities are well tolerated, but when these time-variations are not negligible compared to the dynamics of the controlled process, or when synchronicity between the components is imperative (as in many communication processes), they must be considered in the control design, otherwise the control action may be severely impaired.

      Networked control systems (NCS) are systems where the exchange of information between the components is performed through a network. A NCS reduces installation and maintenance costs. However, this sharing may produce asynchronism in the information processed by the controllers. Such asynchronism can be characterized by network-induced delays, data packet dropout, multi-rate sampling, varying sampling rates, among other less frequent sources of asynchronism[14]. An important challenge in the context of digital control systems and in particular for NCS is the modeling of such asynchronisms. When ignored, they compromise the system performance and even its stability[5]. In [1], several types of asynchronisms are described and a model for network-induced delays is presented.

      Several types of asynchronisms in dynamical discrete-time linear systems and commonly found in NCS lead to the so called discrete-time switched systems (DTSS), where the state vectors are given by the vector sequence defined by the discrete linear inclusion:

      $$ { x}_{k+1}\in\{A_{i,k}{ x}_k\,|\,A_{i,k}\in{\cal M}\},\quad k = 0,1,\cdots $ $ (1)

      where $ { x}_k\in {\bf R}^n $, $ A_{i,k}\in {\bf R}^{n\times n} $ and $ {\cal M}\subset {\bf R}^{n\times n} $ is a set of transition matrices[6].

      A discrete-time switched system given by (1) is said to be stable if all its trajectories converge to the origin. This condition is equivalent to requiring all infinite products of matrices taken from the set $ {\cal M} $ converge to zero[7]. An important feature of the DTSS is that all the matrices $ A_i\in{\cal M} $ can be stable (which implies that their eigenvalues are within the unit disc), but it does not mean that the DTSS is stable[8].

      A DTSS can arise from a dynamical continuous-time linear system sampled at variable time intervals, in this case the transition matrices $ A_{i,k} $ depend on the sampling period $ h_k $. A number of works (see [712] and references cited in [5]), consider linear systems sampled at time rates belonging to a finite set of allowable and known sampling periods, i.e., $ h_k\in\{h_1;\cdots;h_p\} $. In these works, state feedback controllers are designed in order to stabilize the DTSS for all the possible sampling periods. Pazos et al.[13] extend this approach by allowing the sampling period to belong to a continuous time interval, i.e., $ h_k\in[h_{\min},h_{\max}] $; this case implies that the set $ {\cal M} $ has infinite elements. Liberzon and Tempo[14] address the stability of continuous-time switched systems (CTSS) with an infinite number of transition matrices analyzing the existence of a common quadratic Liapunov function. Hetel et al.[5] present a discrete-time model that can simultaneously treat non-uniform sampling and time varying delays; this representation allows the application of several classical stability/stabilizability results from the domain of robust control. In [15], a discrete-time control system with nonuniform sampling is considered, suficient conditions for the closed-loop exponential stability are established and a controller design method, formulated as a convex optimization problem with linear matrix inequalities (LMI) constraints, is given. Hua et al.[16] address the finite-time stabilization problem for a class of NCS with short time-varying delays and sampling jitter; conditions for the design of robust controllers are expressed in terms of LMI. Sun[17] considers the stabilization of large-scale DTSS subject to parameter uncertainties, which include uncertain sampling rates. In [18, 19], sliding mode controllers capable of stabilizing stochastic discrete-time switched systems subject to time-varying delays are addressed.

      In this paper, stabilization and performance improvement of dynamical discrete-time linear systems sampled at varying sampling rates are addressed by means of a state feedback controller designed using the Tikhonov regularization method. The proposed controller is easy to implement. In addition, it reduces an error function so that the state vector is closer to the target state (noise-free state) than it would be if no regularization method were applied. In addition, at each iteration a new optimal regularization parameter is calculated, which results in an adaptive gain controller more efficient than the use of static feedback controllers, as proposed in [810, 13], among many other references.

      The paper is organized as follows. Section 2 presents some preliminaries and summmarizes the results presented in recently published works that address the use of the Tikhonov regularization method in control systems. Section 3 describes the method proposed here and proves its efficiency by reducing an error function. Section 4 shows the results of computational experiments that illustrate the performance of the method. Finally, Section 5 presents some general conclusions.

    • The design of robust feedback controllers capable of attenuating several types of disturbances that affect discrete-time control systems is well known in the literature. However, not much attention is given to the use of the Tikhonov regularization method applied in stabilization and performance improvement of control systems. Some exceptions can be found in [20-22]. Next, we summarize the approach addressed in these works.

      Let the dynamical continuous-time linear system (CTLS) be given by

      $$ \begin{split} & \dot{{ x}} = A(t){ x}(t)+B(t){ u}(t)\\ &{ x}(t_0) = { x}_0 \end{split} $ $ (2)

      where $ t\in J: = [t_0,t_1]\subset[0,\infty) $, $ { x}\in L_2(J,V) $ is the state variable, $ V $ is a Hilbert space possibly in infinite dimensional setting (see [23]), $ { x}_0\in V $, $ A\,:J\times V\rightarrow V $ is a densely defined closed linear operator on $ V $ which generates a strongly continuous semi-group $ {\cal T}: = \{T(t)\,|\,t\geq 0\} $ (about semi-group theory see [24, 25] and references therein), $ B\,:J\times U\rightarrow V $ is a bounded linear operator, $ U $ is a Hilbert space and $ { u}\in L_2(J,U) $ is called a control signal. Note that the control signal is bounded on the interval $ J $. The inner product and the corresponding norm on a Hilbert space are denoted by $ \langle\cdot,\cdot\rangle $ and $ \|\cdot\| $, respectively. The mild solution of (2) is given by

      $$ { x}(t) = T(t-t_0){ x}_0+\int_{t_0}^t T(t-s)B(s){ u}(s)\partial s $ $ (3)

      where $ t\in J $, $ t>t_0 $.

      In the case that $ V $ is a finite dimensional Hilbert space ${\bf C}^n $, the solution of (2) from a given $ t_0 $ is given by [26]:

      $$ { x}(t) = \varPhi(t,t_0){ x}_0+\int_{t_0}^t\varPhi(t,s)B(s){ u}(s)\partial s $ $ (4)

      where $ \varPhi(t,t_0):J\rightarrow {\bf C}^{n\times n} $ is the state transition matrix of $ \dot{{ x}} = A(t){ x} $ (see [26]). In the case that (2) is a linear time invariant (LTI) system, with $ A\in {\bf C}^{n\times n} $ being a constant matrix, $ \varPhi(t,s) = {\rm e}^{A(t-s)} $.

      The control system (2) is said to be exactly controllable over $ J $ if for every $ { x}_\tau\in V $, $ \tau\in J $, $ \tau>t_0 $, there exists $ { u}\in L_2(J,U) $ such that the mild solution (3) satisfies the condition $ { x}(\tau) = { x}_\tau $[24, 26]. When (2) is a LTI system with A and B being constant matrices, it is controllable if and only if the Kalman′s condition:

      $$ {{\rm{rank}}} [B\;AB\;A^2B\;\cdots\;A^{n-1}B] = n $ $

      is satisfied.

      The control system (2) is said to be approximately controllable over $ J $ if for every $ \epsilon>0 $, $ { x}_\tau\in V $, $ \tau\in J $, $ \tau>t_0 $, there exists $ { u}\in L_2(J,U) $ such that the mild solution (3) satisfies the condition[20, 24, 25]:

      $$ \|{ x}(\tau)-{ x}_\tau\|<\epsilon . $ $

      For a given $ t_0,\,\tau\in J $, $ \tau>t_0 $, let us define a control linear operator $ {\cal C}:\,L_2(J,U)\rightarrow V $ as

      $$ {\cal C}{ u}: = \int_{t_0}^\tau T(\tau-s)B(s){ u}(s)\partial s . $ $ (5)

      If $ V $ and $ U $ are finite dimensional Hilbert spaces $ {\bf C}^n $ and $ {\bf C}^m $, respectively, the $ C_0 $ semi-group $ T(\tau-s) $ is substituted by the matrix $ \varPhi(\tau,s) $[24]. The system (2) is also said to be controllable if and only if there exists $ { u}(t)\in L_2(J,U) $ such that

      $$ {\cal C}{ u} = { v} $ $ (6)

      where $ { v}: = { x}_\tau-T(\tau-t_0){ x}_0\;\in V $ (equal to ${ x}_\tau- $$ \varPhi(\tau,t_0){ x}_0 $ in the finite dimensional case), i.e., $ \mbox{ran}\; ({\cal C}) = V $[25].

      In the case of (2) to be controllable, the problem to find a control signal $ { u}(t) $ able to lead the state vector $ { x} $ from $ { x}_0 $ to $ { x}_\tau $ in a finite time interval $ [t_0,\tau]\subseteq J $ is reduced to solve the operator linear (6). The function $ { u}(t) $ that satisfies (6) can be expressed as[24]

      $$ { u}: = {\cal C}^*({\cal C}{\cal C}^*)^{-1}{ v} $ $ (7)

      where $ {\cal C}^* $ is the adjoint operator of $ {\cal C} $. An adjoint operator between Hilbert spaces is defined as follows[23, 25]: let $ X,\,Y $ be Hilbert spaces and $ T:X\rightarrow Y $ be a bounded linear operator, then there exists a unique bounded linear operator $ T^*:Y\rightarrow X $ such that for all $ { x}\in X,\,{ y}\in Y $,

      $$ \langle T{ x},{ y}\rangle = \langle{ x},T^*{ y}\rangle .$ $

      At this point, let us see how $ {\cal C}^* $ looks in the finite dimensional case. For every $ { r}\in V $,

      $$ \begin{split} \langle{\cal C}{ u},{ r}\rangle=& \int_{t_0}^\tau\langle \varPhi(\tau,s)B(s){ u}(s),{ r}\rangle\partial s =\\ & \int_{t_0}^\tau\langle{ u}(s),B^*(s)\varPhi^*(\tau,s){ r}\rangle\partial s =\\ & \langle{ u},{\cal C}^*{ r}\rangle . \end{split} $ $

      Thus,

      $$ ({\cal C}^*{ r})(s) = B^*(s)\varPhi^*(\tau,s){ r} $ $

      and the non-negative linear operator $ {\cal C}{\cal C}^*:V\rightarrow V $ can be expressed as

      $$ {\cal C}{\cal C}^*{ r} = \int_{t_0}^\tau\varPhi(\tau,s)B(s)B^*(s)\varPhi^*(\tau,s){ r}\partial s . $ $

      The operator $ {\cal C}{\cal C}^* $ (also denoted as controllability Gramian[26]) is positive definite if and only if the system (2) is controllable over $ J $[24]. In the case that (2) is LTI with $ A\in {\bf C}^{n\times n} $ and $ B\in {\bf C}^{n\times m} $, this operator can be expressed as the nonsingular matrix

      $$ W_c: = \int_{t_0}^\tau {\rm e}^{A(\tau-s)}BB^*{\rm e}^{A^* (\tau-s)}\partial s $ $ (8)

      and the control signal (7) as

      $$ { u}(t): = B^* {\rm e}^{A^*(\tau-t)}W_{\rm c}^{-1}{ v} $ $ (9)

      where $ { v} = { x}_\tau-{\rm e}^{A(\tau-t_0)}{ x}_0 $.

      When the system (2) is not exactly controllable, the control signal that minimizes the error $ \|{ x}(\tau)-{ x}_\tau\| $ can be found solving the unconstrained optimization problem

      $$ \min\|{\cal C}{ u}-{ v}\| = \mathop {\inf }\limits_{{{g}} \in {L_2}(J,U)} \|{\cal C}{ g}-{ v}\|.$ $

      The argument that minimizes this problem $ { u}(t)\in L_2(J,U) $ is called the least squares control signal of the system (2) and this function satisfies[24]

      $$ {\cal C^*C}{ u} = {\cal C}^*{ v} . $ $ (10)
    • According to [21], approximately controllable linear systems are ill-posed in the sense of Tikhonov. Roughly speaking, ill-posed problems can be defined as those where small perturbations in the observed data (due to noise or inaccurate measurements, for example) may produce large effects on the computed solutions. In these cases, a regularization method should be applied in order to dampen or filter out the noise effects. Regularization means to modify the original ill-posed problem by a nearby better-posed one that is less sensitive to disturbances, in such a way that the solution of the regularized problem is closer to the noise-free solution than the exact solution found solving the original ill-posed problem[27, 28]. Regularization methods require the choice of a scalar regularization parameter, and their effectiveness depends strongly on this choice. The choice of the regularization parameter searches to balance the error due to noise with the error due to regularization. The most used selection methods of the regularization parameter are based on the Morozov discrepancy principle[29], the generalized cross validation function[30], and the L curve[31, 32]. Pazos and Bhaya[33] present an adaptive choice of the Tikhonov regularization parameter based on the minimization of a conveniently chosen Liapunov function using gradient descent algorithms. In [34], a survey of different methods to calculate the regularization parameter is presented. Other parameter choice methods can be found in [3539] among many other references.

      A usual need for most methods of parameter selection is knowledge of the singular values of the operator, which requires a computationally expensive procedure when applied to large-scale problems. However, finite dimensional dynamical control systems like (2) usually have small dimensions, and therefore the aforementioned parameter choice methods are applicable in these systems without much computational effort.

      In the following, we describe the most usual form of the Tikhonov regularization method applied to linear systems.

      In the case of the discrete ill-posed linear problem described by the unconstrained minimization problem $ \min\|A{ x}-{ b}\|^2 $, where $ A\in {\bf R}^{m\times n} $, $ { b}\in {\bf R}^m $ and $ { x}\in {\bf R}^n $ is the vector of unknowns, its minimizing argument is given by the vector that satisfies the normal equation $ A ^{\rm{T}} A{ x} = A ^{\rm{T}}{ b} $. The ill-posedness of the problem is determined by the condition number of the matrix $ A ^{\rm{T}} A $: perturbing the vector $ { b} $ could result in a solution $ { x} $ far from the expected solution. The most usual form of the Tikhonov regularization method applied to this problem seeks to minimize the functional $ \|A{ x}-{ b}\|^2+\lambda\|{ x}\|^2 $, or equivalently, to solve

      $$ (A ^{\rm{T}} A+\lambda I){ x} = A ^{\rm{T}}{ b} $ $

      where $ \lambda\in {\bf R}_+ $ is the scalar regularization parameter. It is easy to prove that the condition number of $ (A ^{\rm{T}} A+\lambda I) $ is smaller than the condition number of ATA, so the regularized problem is better conditioned than the original ill-posed normal equation.

      In the case of the least squares control of an ill-posed continuous-time linear system described by (10), its regularization by the most usual form of the Tikhonov method consists in the substitution of (10) by

      $$ ({\cal C^*C}+\lambda I){ u} = {\cal C}^*{ v} $ $ (11)

      where $ I:U\rightarrow U $ is the identity operator in $ U $[23].

      Nair et al.[20] apply the Tikhonov regularization method to an approximately controllable infinite dimensional linear system exemplified by the following controlled heat equation:

      $$ \begin{split} & \frac{\partial z}{\partial t}(s,t) = \frac{\partial^2 z}{\partial s^2}(s,t)+u(s,t),\quad s\in[0,\ell]\\ &z(0,t) = z(\ell,t) = 0,\quad\forall t>0\\ &z(s,0) = z_0(s) \end{split} $ $ (12)

      where $ z(s,t) $ represents the temperature along a rod of length $ \ell $ at position $ s $ at time $ t $, $ z_0 \in L_2([0, \ell], {\bf R}) $ is the initial temperature profile and $ u(s,t) $ is the heat input along the rod.

      The upper bound of the error committed by regularization of the problem is calculated in [20]. The regularization parameter $ \lambda $ is calculated in order to keep the error due to regularization smaller than a predetermined limit derived from a tolerance $ \epsilon $ (see details in [20]).

      In [21], a weighted Tikhonov regularization method is proposed to solve approximately controllable linear systems. The resulting well-posed regularized operator equation is given by

      $$ \left(({\cal C}^*{\cal C})^{\frac{\alpha+1}{2}}+\lambda I\right){ u} = \left({\cal C}^*{\cal C}\right)^{\frac{\alpha-1}{2}}{\cal C}^*{ v} $ $ (13)

      where $ \alpha\in(0,1] $ is an user-defined parameter and $ \lambda\in {\bf R}_+ $ is the regularization parameter. Note that if $ \alpha = 1 $, (13) coincides with the classical form of the Tikhonov regularization method (11). The authors apply the Morozov discrepancy principle[29] to calculate the regularization parameter $ \lambda $. The regularization method (13) is also applied to the approximately controllable heat control system (12).

      Katta and Sukavanam[22] apply the classical form of the Tikhonov method to regularize an approximately controllable nonlinear system. This system is described by (2) plus a nonlinear term $ { f}({ x}(t),t) $ which is Lipschitz continuous on $ { x} $. The method is applied to solve the heat control system (12) plus a term $ \sin(z(s,t)) $ in the right-hand side of the partial differential equation.

    • In this section, we consider the application of the Tikhonov regularization method in the design of a feedback controller for a finite dimensional discrete-time linear control system. The sampling of the continuous-time linear system (2), considering $ A\in {\bf R}^{n\times n} $ and $ B\in {\bf R}^{n\times m} $ as constant matrices, and the application of the control signal to the plant through a zero-order hold circuit[40] yield

      $$ { x}_{k+1} = \varPhi(h){ x}_k+\varGamma(h){ u}_k $ $ (14)

      where $ \varPhi(h) = {\rm e}^{Ah} $, $ \varGamma(h) = \int_0^h{\rm e}^{A\tau}\partial\tau B $ and $ h $ is the sampling period. If $ A $ is a non singular matrix: $\varGamma(h) = $$ A^{-1}({\rm e}^{Ah}-I)B $. From (14):

      $$ \varGamma ^{\rm{T}}({ x}_{k+1}-\varPhi{ x}_k) = \varGamma ^{\rm{T}}\varGamma{ u}_k $ $ (15)

      where the dependence on the sampling period was omitted. Equation (15) is (10) in the finite dimensional discrete-time case. If the linear system expressed as (14) is ill-posed, a small perturbation in its right-hand side (i.e., due to asynchronism or noisy measurements) leads to a state $ { x}_{k+1} $ far from the expected state which will be denoted as $ { x}_{{ d},k+1} $. In this case, a regularization method should be applied in order to attenuate the perturbation effects. From (15), the control signal calculated using the classical form of the Tikhonov regularization method is given by

      $$ { u}_{k,\lambda}: = (\varGamma ^{\rm{T}}\varGamma+\lambda I)^{-1}\varGamma ^{\rm{T}} ({ x}_{{ d},k+1}-\varPhi{ x}_k) $ $ (16)

      in such a way that the application of this control signal to the system (14) leads to the regularized state

      $$ { x}_{k+1,\lambda} = \varPhi{ x}_k+\varGamma{ u}_{k,\lambda}.$ $ (17)

      The goal of the regularization method is to achieve that the regularized state is closer to the target state to be attained than the unregularized state obtained from (14), i.e., $ \|{ x}_{k+1,\lambda}-{ x}_{{ d},k+1}\|<\|{ x}_{k+1}-{ x}_{{ d},k+1}\| $, for an adequately chosen parameter $ \lambda $.

      Fig. 1 shows a block diagram of the linear system where the transference function of the controller H is given by $ (\varGamma ^{\rm{T}} \varGamma+\lambda I)^{-1}\varGamma ^{\rm{T}} $ in the case of the Tikhonov regularization method, and by $ (\varGamma ^{\rm{T}}\varGamma)^{-1}\Gamma ^{\rm{T}} $ in the case of the least squares control.

      Figure 1.  Block diagram of the discrete-time linear system

      In the following, we consider the system (14) with a perturbation given by a variable sampling interval, such that $ h = \bar{h}+\Delta h $, where $ \bar{h} $ is the nominal sampling interval and $ \Delta h $ is the disturbance. We assume $ |\Delta h|\ll\bar{h} $. Of course, $ \Delta h $ can be positive or negative. Note that by considering a varying sampling interval $ h $, (14) becomes a time variant system.

      Note also that, if a control signal calculated using a static feedback controller $ { u}_k = L{ x}_k $, $ L\in {\bf R}^{m\times n} $ is applied, system (14) becomes a DTSS as described in (1). As noted in [8], the transition matrix $ \varPhi(h)+\varGamma(h)L $ can have all its eigenvalues within the unit disc for all allowable sampling interval h, but it does not mean that the system is stable. Variations in the sampling rate can potentially unstabilize the controlled system.

      The desired state $ { x}_{{ d},k+1} $ is given by the left hand side of (14) when $ \Delta h = 0 $. As $ \Delta h $ is usually unknown, the matrices in (16) are calculated using the nominal sampling interval $ \bar{h} $.

      Note that if $ \lambda = 0 $ then $ { u}_{k,\lambda} = { u}_k $ and $ { x}_{k+1,\lambda} = { x}_{k+1} $.

      Next, we show that for an adequately chosen regularization parameter $ \lambda $, the error of the regularized system (16) and (17) is smaller than that presented by the unregularized one given by (14).

      Let us introduce the following nomenclature:

      $$ \begin{split} \varGamma(h) =& \int_0^{\bar{h}+\Delta h}{\rm e}^{A\tau}\partial\tau B = \int_0^{\bar{h}} {\rm e}^{A\tau}\partial\tau B+\\ & \int_{\bar{h}}^{\bar{h}+\Delta h}{\rm e}^{A\tau}\partial\tau B := \varGamma(\bar{h})+\Delta\varGamma . \end{split} $ $

      The assumption $ |\Delta h|\ll\bar{h} $ also implies that $\|\Delta\varGamma\|\ll $$ \|\varGamma(\bar{h})\| $, whatever the matrix norm is used. At the k-th iteration, the error of the regularized system is given by

      $$ \begin{split} &\epsilon_k(\lambda): = \|{ x}_{k+1,\lambda}-{ x}_{{ d},k+1}\| =\\ & \quad\|\varPhi(h){ x}_k+\varGamma(h){ u}_{k,\lambda}-(\varPhi(\bar{h}){ x}_k+\varGamma(\bar{h}){ u}_k)\| =\\ & \quad \|({\rm e}^{Ah}-{\rm e}^{A\bar{h}}){ x}_k +\varGamma(h) \times\\ &\quad \left((\varGamma(\bar{h}) ^{\rm{T}}\varGamma(\bar{h})+\lambda I)^{-1}\varGamma(\bar{h}) ^{\rm{T}}({ x}_{{ d},k+1}-\varPhi(\bar{h}){ x}_k)\right)-\\ &\quad \varGamma(\bar{h}){ u}_k\|= \|({\rm e}^{Ah}-{\rm e}^{A\bar{h}}){ x}_k+\\ & \quad \left(\varGamma(\bar{h})+\Delta\varGamma\right)\left((\varGamma(\bar{h}) ^{\rm{T}}\varGamma(\bar{h}) +\lambda I)^{-1}\varGamma(\bar{h}) ^{\rm{T}}\right)\times \\ & \quad \varGamma(\bar{h}){ u}_k-\varGamma(\bar{h}){ u}_k\| = \|({\rm e}^{Ah}-{\rm e}^{A\bar{h}}){ x}_k+\\ & \quad \left[(\varGamma+\Delta\varGamma)(\varGamma ^{\rm{T}}\varGamma+\lambda I)^{-1}\varGamma ^{\rm{T}}\varGamma-\varGamma\right]{ u}_k\| \end{split} $ $ (18)

      where the dependence on $ \bar{h} $ was omitted in the last row. For $ \lambda = 0 $, this error is equal to $ \|{ x}_{k+1}-{ x}_{{ d},k+1}\| $, i.e., the error introduced by disturbances in the sampling rate. If there is no error in the sampling rate, i.e., $ \Delta h = 0 $, the error is equal to that introduced by the use of regularization $ \epsilon_k(\lambda) = \|\left[\varGamma(\varGamma ^{\rm{T}}\varGamma+\lambda I)^{-1}\varGamma ^{\rm{T}}\varGamma-\varGamma\right]{ u}_k\| $, which, in turn, is obviously zero when $ \lambda = 0 $. The optimal value of $ \lambda $ is the one that balances these two sources of error, minimizing (18).

      The error $ \epsilon_k(\lambda) $ is a quasiconvex scalar function for all $ \lambda\geq 0 $. For small values of $ \lambda $, the error decreases monotonically to a minimum value as $ \lambda $ increases from $ 0 $.

      The optimal value of $ \lambda $ which minimizes (18), calculated by zeroing its derivative, is given by

      $$ \begin{split} \lambda_k^* =& \frac{{ u}_k ^{\rm{T}}\varGamma ^{\rm{T}}\varGamma\left[(\varGamma+\Delta\varGamma) ^{\rm{T}}(\varGamma+\Delta\varGamma)\right]^{-1}(\varGamma+\Delta\varGamma) ^{\rm{T}}(\varGamma{ u}_k-({\rm e}^{Ah}-{\rm e}^{A\bar{h}}){ x}_k)}{\|\left[(\varGamma+\Delta\varGamma) ^{\rm{T}}(\varGamma+\Delta\varGamma)\right]^{-1}(\varGamma+\Delta\varGamma) ^{\rm{T}}(\varGamma{ u}_k-({\rm e}^{Ah}-{\rm e}^{A\bar{h}}){ x}_k)\|^2}-\\ &\frac{\|\varGamma\left[(\varGamma+\Delta\varGamma) ^{\rm{T}}(\varGamma+\Delta\varGamma)\right]^{-1}(\varGamma+\Delta\varGamma) ^{\rm{T}}(\varGamma{ u}_k-({\rm e}^{Ah}-{\rm e}^{A\bar{h}}){ x}_k)\|^2}{\|\left[(\varGamma+\Delta\varGamma) ^{\rm{T}}(\varGamma+\Delta\varGamma)\right]^{-1}(\varGamma+\Delta\varGamma) ^{\rm{T}}(\varGamma{ u}_k-({\rm e}^{Ah}-{\rm e}^{A\bar{h}}){ x}_k)\|^2} . \end{split} $ $ (19)

      Note that if $ \Delta h = 0 $, i.e., there is no uncertainty in the sampling interval, $ h = \bar{h} $, $ \Delta\varGamma = { 0} $, and the optimal value of $ \lambda_k^* $ in (19) is $ 0 $ which implies $ \epsilon_k(0) = 0 $.

      Of course, as $ \Delta h $ is not previously known, it is not possible to calculate the optimal value of $ \lambda_k^* $ given by (19) in order to calculate the optimal control signal $ { u}_{k,\lambda} $ that minimizes the error due to uncertainties in the sampling rate. Other parameter choice methods as those mentioned in Section 2.1 must been used instead.

      It is easy to see from (19) that in the scalar case ($ m = n = 1 $), the value of $ \lambda $ that minimizes the error is given by

      $$ \lambda_k^* = \gamma^2\,\frac{\Delta\gamma\, u_k+({\rm e}^{a h}-{\rm e}^{a\bar{h}})x_k}{\gamma\, u_k-({\rm e}^{a h}-{\rm e}^{a\bar{h}})x_k} $ $

      where we denote as $ \gamma $ and a to be the matrices $ \varGamma $ and A, respectively, in the scalar case. In this case, $\Delta\gamma = $$ \dfrac{b}{a}({\rm e}^{a h}-{\rm e}^{a\bar{h}}) $ and $ \gamma = \dfrac{b}{a}({\rm e}^{a\bar{h}}-1) $, where b is the matrix B in the scalar case.

      Note also that the ill-posedness of the linear system (15) is characterized by a large condition number of the matrix $ \varGamma ^{\rm{T}}\varGamma $. However, the application of the Tikhonov regularization method with an adequately chosen parameter $ \lambda $ reduces the error produced by uncertainties in the sampling rate despite this condition number.

      With the result presented above, we enunciate Lemma 1:

      Lemma 1. Consider the discrete-time linear system (14) with a varying sampling rate $ \bar{h}+\Delta h $. When a control signal $ { u}_k $ is applied at an iteration $ k $, after an uncertain sampling period, the system produces a state vector $ { x}_{k+1} $ which is different from the desired state $ { x}_{{ d},k+1} $ (which is the state produced by the same control signal when $ \Delta h = 0 $). The error committed is given by $ \epsilon_k(0) $ in (18). If a control signal $ { u}_{k,\lambda} $ given by (16) is applied instead, with a value of $ \lambda $ between $ 0 $ and that expressed in (19), the error $ \epsilon_k(\lambda) $ decreases, reaching a minimum value when $ \lambda $ is equal to that calculated in (19).

    • In this section, we present the results of the regularized control signal applied to a discrete-time linear systems with varying sampling rates. The simulations were carried out using Matlab 2017.

      The following example is a paper machine head box model; it was extracted from [41] and it was used in [13]. Let a dynamical continuous-time linear system be given by

      $$\begin{split} &\dot{{ x}} = A{ x}+B{ u}= \qquad\\ & \qquad\left[\begin{array}{ccc}-0.2\,&0.1\,&1\\-0.05&0&0\\0&0&-1\end{array}\right]{ x}+\left[\begin{array}{cc}0\,&1\\0&0.7\\1&0\end{array}\right]{ u} .\end{split} $ $ (20)

      System (20) is a stable system with eigenvalues $ \sigma\in\{ -0.170\,7,\, -0.029\,3,\, -1\} $. System (20) is sampled at a varying sampling rate with Gaussian distribution of mean $ \bar{h} = 1\,{\rm s}$ and standard deviation equal to $ 0.4\,{\rm s} $. The resulting discrete time system sampled at the nominal time interval $ \bar{h} $ is stable, the matrix $ \varPhi(\bar{h}) $ has eigenvalues $ \sigma\in\{ 0.843\,1,\, 0.971\,1,\, 0.367\,9\} $. The simulations were carried out during $ 10\,{\rm s}$. The sampling intervals were previously and randomly chosen as $h_k\in\{1.437\,3\,{\rm s},\, 1.443\,7\,{\rm s}, $$ 0.654\,5\,{\rm s},\,1.030\,9\,{\rm s},\, 0.514\,4\,{\rm s},\, 0.554\,6\,{\rm s},\, 0.997\,3\,{\rm s},\, 1.613\,1\,{\rm s},\, 0.692\,1\,{\rm s},$$1.148\,6\,{\rm s},\, 0.909\,8\,{\rm s},\, 1.446\,9\,{\rm s}\} $. With these sampling intervals, the instants of time when the state vector was sampled were $ t_k\in\{0,\, 1.437\,3\,{\rm s},\, 2.881\,0\,{\rm s},\, 3.535\,5\,{\rm s},\, 4.566\,4\,{\rm s},\, $ $5.080\,8\,{\rm s},\, 5.635\,4\,{\rm s},\, 6.632\,7\,{\rm s},\, 8.245\,8\,{\rm s},\, 8.937\,9\,{\rm s},\, 10.086\,5\,{\rm s}\} $.

      The initial state was chosen as $ { x}_0 = [2\;-3\;-1] ^{\rm{T}} $.

      In a first series of experiment, the control signal was chosen as $ { u}_k = \left[\sin\left(\pi \dfrac{t_k}{5}\right)\;\;\cos\left(\pi \dfrac{t_k}{3}\right)\right] ^{\rm{T}} $. This control signal is applied to system (20) through a zero-order hold circuit in order to produce the unregularized state vector

      $$ { x}_{k+1} = \varPhi(h_k){ x}_k+\varGamma(h_k){ u}_k. $ $

      The desired state vector is calculated as if there were no uncertainties in the sampling interval, i.e.,

      $$ { x}_{{ d},k+1} = \varPhi(\bar{h}){ x}_k+\varGamma(\bar{h}){ u}_k $ $

      which is used to calculate the regularized control signal $ { u}_{k,\lambda} $ in (16). This control signal applied to system (20) through a zero-order hold circuit produced the regularized state vector $ { x}_{k+1,\lambda} $ given by (17).

      In this first series of experiments, the regularization parameter $ \lambda $ has been chosen by hand tuning as a constant equal to $ 0.03 $.

      Fig. 2 shows the components of the regularized control signal $ { u}_{k,\lambda} $ versus time.

      Figure 2.  Components of the regularized control vector $ { u}_{k,\lambda} $ versus time

      Fig. 3 shows the components of the desired state vector $ { x}_{{ d},k} $, the unregularized state $ { x}_k $ and the regularized state $ { x}_{k,\lambda} $ versus time. The last two vectors are the output of the continuous-time system (20) if the inputs $ { u}_k $ and $ { u}_{k,\lambda} $ were applied, respectively.

      Figure 3.  State vectors versus time. The solid lines are the components of the discrete target state. The dotted lines are the output of the continuous-time plant if the unregularized control signal $ { u}_k $ is applied. The dashed lines are the output of the continuous-time plant if the regularized control signal $ { u}_{k,\lambda} $ is applied.(Color versions of the figures in this paper are available online.)

      At this point, let us define the mean error produced throughout all the iterations as

      $$ \bar{\epsilon}(\lambda): = \frac{\sum\limits_{k = 1}^N\epsilon_k(\lambda)}{N} = \frac{\sum\limits_{k = 1}^N\|{ x}_{k,\lambda}-{ x}_{{ d},k}\|}{N} $ $ (21)

      where N is the final iteration. In this experiment, carried out during $ 10\,{\rm s} $, $ N = 10 $.

      The mean error depends on the value of the regularization parameter. Fig. 4 shows the mean error (21) versus the regularization parameter $ \lambda $. The curve presents a minimum value at $ \lambda = 0.03 $. As mentioned in Section 3, when $ \lambda = 0 $, $ { u}_{k,\lambda} = { u}_k $ and $ { x}_{k+1,\lambda} = { x}_{k+1} $, so the error $ \epsilon_k(0) $ calculated in (18), as well as the mean error $ \bar{\epsilon}(0) $ are equal to those presented by the unregularized state vector.

      Figure 4.  Mean error given by (21) versus the regularization parameter $ \lambda $

      An open question is how to choose a suitable value of the regularization parameter. As mentioned in Section 3, (19) cannot be used in practical applications because the value of $ \Delta h $ is not previously known. The methods mentioned in Section 2.1 were also tested (see a description about these methods in [2931], respectively). The minimum of the generalized cross validation function was calculated using the function gcv, the corner of the L curve was calculated using the function l_curve, and the value of $ \lambda $ that satisfies the discrepancy principle was calculated using the function discrep. All these functions are available in Hansen′s regularization toolbox for Matlab[42]. As the left hand side of (15) changes at each iteration, a different value of $ \lambda $, instead a constant one as used in our test, is calculated by these functions at each iteration. This results in an adaptive gain controller that is more efficient than a static state feedback controller. Table 1 shows the mean error for different constant values of $ \lambda $ as well as the mean error produced when this parameter is calculated by the three functions aforementioned and by choosing $ \lambda_k^* $ calculated as in (19) at each iteration. Note that although the unregularized state vector does not depend on $ \lambda $, the desired state vector does, because it is calculated using the regularized state at the former iteration, so the mean error of the unregularized system also depends on $ \lambda $.

      Table 1.  Mean error of the regularized state vector $\sum_{k = 1}^N \dfrac{\|{ x}_{k,\lambda}-{ x}_{{ d},k}\|}{N}$ and mean error of the unregularized state vector $\sum_{k = 1}^N\dfrac{\|{ x}_{k}-{ x}_{{ d},k}\|}{N}$ for different values of the parameter $\lambda$ in the first series of experiments

      $\lambda$00.030.06$\lambda_k^*$gcvl_curvediscrep
      Regularized state vector0.210 60.207 10.208 40.064 50.311 20.311 20.327 7
      Unregularized state vector0.210 60.224 40.246 60.128 70.503 80.503 80.604 0

      The difference between the mean error of the regularized state vector (produced when $ { u}_{k,\lambda} $ is applied), and the mean error of the unregularized state vector (produced when $ { u}_k $ is applied) quantifies the performance improvement of the proposed method. In Table 1, the largest difference between the mean errors is produced when the function discrep is used.

      Note also that to use a parameter $ \lambda_k^* $ that minimizes $ \epsilon_k(\lambda) $ different at each iteration using (19) produces the minimum mean error.

      In a second series of experiments, we use as control signal that presented in [13]: $ { u}_k = L{ x}_k $, where the state feedback controller is chosen as

      $$ L = \left[ {\begin{array}{*{20}{c}} {0\;\;}&{0\;\;}&{ - 1}\\ {0.5}&{ - 1.707\,1}&1 \end{array}} \right] . $ $

      Thus, the target state vector is given by

      $$ { x}_{{ d},k+1} = \left(\varPhi(\bar{h})+\varGamma(\bar{h})L\right){ x}_{k}. $ $

      The unregularized state vector is given by

      $$ { x}_{k+1} = \left(\varPhi(h)+\varGamma(h)L\right){ x}_k. $ $

      The regularized state vector $ { x}_{k+1,\lambda} $ is given by (17) and the regularized control signal $ { u}_{k,\lambda} $ is

      $$ \begin{split} { u}_{k,\lambda} = &\left(\varGamma ^{\rm{T}}\varGamma+\lambda I\right)^{-1}\varGamma ^{\rm{T}}({ x}_{{ d},k+1}-\varPhi{ x}_k) =\\ &\left(\varGamma ^{\rm{T}}\varGamma+\lambda I\right)^{-1}\varGamma ^{\rm{T}}\varGamma L{ x}_k \end{split} $ $

      where the matrix $ \varGamma $ depends on $ \bar{h} $.

      The initial state and the sampling times were the same that those used in the former test. The test was also carried out during $ 10\,{\rm s} $, so $ N = 10 $ in (21).

      The regularized parameter was chosen by hand tunning as a constant equal to $ \lambda = 0.5 $.

      Fig. 5 shows the components of the regularized control signal versus time.

      Figure 5.  Components of the regularized control vector versus time

      Fig. 6 shows the components of the target state vector, the unregularized state vector, and the regularized state vector versus time.

      Figure 6.  State vectors versus time. The solid lines are the components of the discrete target state, the dotted lines are the output of the continuous-time plant if the unregularized control signal $ { u}_k $ is applied, the dashed lines are the output of the continuous-time plant if the regularized control signal $ { u}_{k,\lambda} $ is applied.

      Fig. 7 shows the mean error of the regularized system versus the regularization parameter $ \lambda $. The plot shows that the value of $ \lambda $ that minimizes the mean error is $ 0.5 $.

      Figure 7.  Mean error of the regularized state vector $ \bar{\epsilon} $ versus the regularization parameter $ \lambda $

      Here again, we test different constant values of $ \lambda $, as well as the values of $ \lambda $ calculated by the functions gcv, l_curve and discrep and the value of $ \lambda_k^* $ calculated as in (19) at each iteration. The mean error presented by the regularized state vector and by the unregularized state vector are reported in Table 2.

      Table 2.  Mean error of the regularized state vector $\sum_{k = 1}^N\dfrac{\|{ x}_{k,\lambda}-{ x}_{{ d},k}\|}{N}$ and mean error of the unregularized state vector $\sum_{k = 1}^N\dfrac{\|{ x}_{k}-{ x}_{{ d},k}\|}{N}$ for different values of the parameter $\lambda$ in the second series of experiments

      $\lambda$00.250.5$\lambda_k^*$gcvl_curvediscrep
      Regularized state vector0.461 80.301 10.214 60.160 50.232 70.232 70.318 3
      Unregularized state vector0.461 80.484 70.686 00.572 10.618 20.618 20.787 4

      Once again, as it happened in the first series of experiments, the choice of an optimal regularization parameter $ \lambda_k^* $ calculated as in (19) produces the minimum mean error. The second largest difference between the mean errors produced by the regularized state vector and by the unregularized state vector is produced when the regularization parameter $ \lambda $ is calculated using the function discrep. Note also that the choice of a constant regularization parameter close to the one which minimizes the curve plotted in Fig. 7 produces a mean error of the regularized state vector smaller than that produced when a variable regularization parameter is calculated using any of the three standard functions.

      Algorithm 1 shows how to implement the proposed method in detail.

      Algorithm 1. Tikhonov regularization method applied to discrete-time linear control systems subject to varying sampling rates.

      1) for k:=0 to N do

      2)  Given a state vector $ { x}_k $ and a control signal $ { u}_k $

      3)  Calculate $ { x}_{{ d},k+1}=\varPhi(\bar{h}){ x}_k+\varGamma(\bar{h}){ u}_k $

      4)  Calculate $ \lambda $ using a standard parameter choice method

      5)  Calculate the regularized control signal $ { u}_{k,\lambda}=\left(\varGamma(\bar{h}) ^{\rm{T}}\varGamma(\bar{h})+\lambda I\right)^{-1}\varGamma(\bar{h}) ^{\rm{T}}\left({ x}_{{ d},k+1}-\varPhi(\bar{h}){ x}_k\right) $ (16)

      6)  Apply the regularized control signal to the system in order to obtain after an uncertain sampling period $ h $ the regularized state ${ x}_{k+1,\lambda}= $$ \varPhi(h){ x}_k+\varGamma(h){ u}_{k,\lambda} $

      7) end for

    • The use of the Tikhonov regularization method to design a state feedback controller to be applied in a discrete-time linear control system in order to attenuate the error produced by uncertainties in the sampling interval is successful for an adequately chosen regularization parameter. Lemma 1 proves that at each iteration there exists a value of $ \lambda $ that minimizes the error. The tests presented in Section 4 illustrate the efficiency of the proposed method. Since the previous knowledge of the regularization parameter $ \lambda $ which minimizes the error $ \epsilon_k $ at each iteration is not possible, as well as the value of $ \lambda $ which minimizes the mean error $ \bar{\epsilon} $, a standard parameter choice method must be applied. Among the three well known methods tested in the experiments reported here, based on the minimization of the cross validation function, the search for the corner of the L curve, and the calculation of the parameter which satisfies the discrepancy principle, the first two methods seem to present a smaller mean error of the regularized state vector than the last method (in both of the tests the gcv and the l_curve methods produced the same results), whereas the last method presents the largest difference in relation to the mean error produced by the application of the unregularized control signal, as reported in Tables 1 and 2. However, more tests and an analytical study must be done to confirm this preliminary conclusion.

      The error reduction achieved by using regularization also reduces the possibility of having a sequence of transition matrices which makes the DTSS unstable, as noted in [8].

      As future work, we propose to investigate if other type of disturbances, i.e., noise in the measured state $ { x}_k $, could also have its effects attenuated through the use of a regularized control signal $ { u}_{k,\lambda} $. Another future work is to investigate if the regularization scheme proposed here is applicable to discrete-time nonlinear control systems subject to several types of disturbances.

    • The authors would like to thank the editor and the anonymous reviewers for their valuable comments and helpful suggestions that have contributed to the improvement of the paper.

参考文献 (42)

目录

    /

    返回文章
    返回