Robust Adaptive Stabilization of Linear Time-Invariant Dynamic Systems by Using Fractional-Order Holds and Multirate Sampling Controls

Abstract

This paper presents a strategy for designing a robust discrete-time adaptive controller for stabilizing linear time-invariant (LTI) continuous-time dynamic systems. Such systems may be unstable and noninversely stable in the worst case. A reduced-order model is considered to design the adaptive controller. The control design is based on the discretization of the system with the use of a multirate sampling device with fast-sampled control signal. A suitable on-line adaptation of the multirate gains guarantees the stability of the inverse of the discretized estimated model, which is used to parameterize the adaptive controller. A dead zone is included in the parameters estimation algorithm for robustness purposes under the presence of unmodeled dynamics in the controlled dynamic system. The adaptive controller guarantees the boundedness of the system measured signal for all time. Some examples illustrate the efficacy of this control strategy.

1. Introduction

Adaptive control theory has been widely applied for stabilizing increasingly complex engineering systems with large uncertainties [1], including the incorporation of parallel multiestimation and time-delayed and hybrid models [2–6]. Such model uncertainties may come from the fact that the parameters of the dynamic system model are partially or fully unknown and/or from the presence of unmodeled dynamics [3]. On the other hand, discrete equations are useful for modeling and controlling discretized continuous-time systems in practical situations [2, 5–9] and as a tool for describing more complex nonlinear structures via discretization [10, 11]. A frequently used method to stabilize an unknown dynamic system is based on the model reference adaptive control (MRAC) problem [12]. However, the presence of unstable system zeros does not make possible the design of a controller to achieve the model-matching objective unless such unstable zeros are transmitted to the reference model [1, 2, 6, 13, 14]. Unfortunately, such a transmission cannot be available if some of the zeros of the system to be controlled are unknown as it can occur in the context of the adaptive control where the dynamic system may be completely or partially unknown. However, there are several alternative methods to circumvent this difficulty and carry out the stable adaptive control design. Two of them depend on relaxing the control objective from the MRAC to the less exigent adaptive pole-placement control (APPC) [15, 16]. In this way, the stabilization of the feedback (or closed-loop) system can be ensured although its transient behavior cannot be fixed to a predefined one. On one hand, the method in [15] includes a modification in the estimated parameters to ensure the controllability of the estimation model of the dynamic system. In this way, closed-loop unstable pole-zero cancellations are avoided, which is crucial for the controller synthesis. Such a method is applicable for both continuous-time and discrete-time dynamic systems to be controlled. On the other hand, the research [16] may be used in the case of continuous-time dynamic systems. There, a periodic piecewise constant-gain controller is added in the feedback chain. In the nonadaptive case, the gain values are those required so that the discretized system model under the fundamental sampling period and a zero-order hold (ZOH) could be stabilized. It is worth that such a control gain is piecewise constant during the sampling period in order to place the discretized poles at stable desired locations. Concretely, each sampling period is split into a certain finite number of uniform subintervals and the control gain takes a different value within each subinterval. In this way, the controller consists of a constant vector of gains. In this sense, the controller works with a sampling rate faster than that used to discretize the controlled system. In the adaptive case, the discretized dynamic system model parameters are on-line updated; then the controller gains vector is time varying and converges asymptotically to a constant one.

Another method, which does not relax the MRAC objective, to overcome the drawback of the unstable zeros of a continuous-time dynamic system is the design of a discrete-time controller with the use of a hold device combined with a multirate sampling with fast input rate in the discretization of the continuous-time system [7, 8]. In this way, an inversely stable estimated model of the discretized dynamic system can be obtained and a controller can be designed to match a stable arbitrarily chosen discrete-time reference model since all of the discretized zeros may be cancelled if suited. In this context, this paper presents a robust adaptive control scheme for stabilizing uncertain controlled continuous-time dynamic systems while matching a freely chosen discrete-time reference model, with a bounded tracking-error, to be applied when the continuous-time system is unknown and subject to the presence of unmodeled dynamics. The main novelty is that the discrete zeros may be always stabilized even if the zeros of the continuous system are not all stable. As a result, the reference model zeros can all be freely fixed. The control scheme is based on the discretization process by combining the use of a fractional-order hold (FROH) and a multirate with fast sampling control signal. In this way, the estimated discretized model can be guaranteed to be inversely stable by means of a suitable on-line updating of the multirate gains without requirements neither on the stability of the continuous-time zeros nor the size of the sampling period [17]. Such a strategy gives place to hybrid systems where continuous-time and discrete-time dynamics are mixed. Concretely, a discrete-time controller is designed to govern the behavior of a continuous-time dynamical system. Epidemic control of infectious diseases and species population control systems in Ecology are two typical examples of this class of hybrid systems [18–20]. A discrete-time control strategy depending on a pulse vaccination is designed in [18]. Such vaccination pulses work as a discrete-time control signal and the eradication of the diseases is reached provided that the vaccination rate is sufficiently large. The researches [19, 20] study, respectively, the necessary conditions which guarantee the stabilization and permanence of single-species and predator-prey systems populations in their respective habitats. Both are continuous-time dynamic systems which can be described efficiently by means of discrete-time models in order to prescribe its evolution and then to develop discrete-time control strategies to ensure the permanence of the species. All of these systems are subject to eventual changes in the continuous-time dynamics whenever the discrete-time control action takes place. In this sense, they belong to a class of switched systems whose stability and stabilization conditions have been studied in the recent literature [4, 21].

Furthermore, the presented control strategy guarantees the stabilization of the continuous-time dynamic system without any assumption about the stability of its zeros, which had been supposed in previous works [13, 22], and without requiring estimates modification in contrast with previous works on the subject [2, 15]. Such relaxations constitute the main contribution of the present work. A FROH is used since it allows a better accommodation of discrete adaptive techniques to the transient response of discrete-time controlled continuous-time dynamic systems [9]. Furthermore, the estimation algorithm includes a relative adaptation dead-zone to deal with the presence of unmodeled dynamics [14]. Such a dead-zone is crucial to ensure the estimates convergence and the stability of the adaptive control system. The stabilization is guaranteed provided that (1) the continuous-time dynamic system is stabilizable and observable   (2) the size of the unmodeled dynamics is sufficiently small, and (3) such an unmodeled dynamics can be related to the system input by means of additive and/or multiplicative stable transfer functions.

The paper is organized as follows. Section 2 presents the discretization process used to obtain an inversely stable discretized dynamic system model from a possibly inversely unstable continuous-time dynamic system. Section 3 deals with the control design to match a discrete-time reference model at sampling instants for both nonadaptive and adaptive cases. Then, the stability analysis of the designed adaptive control algorithm is presented in Section 4. Finally, simulation results, which illustrate the behavior of the adaptive control system, are shown in Section 5 and conclusions end the paper in Section 6.

2. Problem Statement

Consider a linear time invariant SISO and strictly proper continuous-time dynamic system described by the following state-space equations:𝑥̇(𝑡)=𝐴𝑥(𝑡)+𝐵𝑢(𝑡),𝑦(𝑡)=𝐶𝑥(𝑡), (2.1) where 𝑢(𝑡) and 𝑦(𝑡) are, respectively, the control (or input) and measured (or output) signals, 𝑥(𝑡)∈ℝ^𝑛 is the state vector, and 𝐴, 𝐵, and 𝐶 are constant matrices of suitable dimensions. The transfer function of (2.1) is 𝐺_𝑝(𝑠)=𝑞(𝑠)/𝑝(𝑠)=𝐶(𝑠𝐼_𝑛−𝐴)⁻¹𝐵 where 𝑛=Deg(𝑝(𝑠))≥Deg(𝑞(𝑠))=𝑚, 𝑠 denotes the Laplace transform argument [1], and 𝐼_𝑛 represents the 𝑛-order identity matrix. In the sequel, the controlled dynamic system is referred to as the “plant” to be controlled as it is commonly referred to in the Engineering Automation context. The following assumptions are made on the plant.

Assumption 1. (i) An upper-bound 𝑛 of the plant order 𝑛 is known as it is the nominal order 𝑛₀≤𝑛.
(ii) The plant realization matrices can be expressed as𝐴𝐴=⁰0𝐴₂₁𝐴₂₂𝐵,𝐵=⁰𝐵₂₁,𝐶=𝐶₀𝐶₁₂, (2.2) where {𝐴₀,𝐵₀,𝐶₀} denotes the state-space realization for the nominal model of the plant and the other matrix blocks are related to the unmodeled dynamics. In this context, the state vector is composed of two sets of state variables, namely, 𝑥≜[𝑥^𝑇₀𝑥^𝑇₁]^𝑇 where 𝑥₀∈ℝ^𝑛0 is the nominal state vector. Furthermore, the eigenvalues of the block 𝐴₂₂ are strictly stable, ‖𝐵₂₁‖≤𝜇₀, and ‖𝐶₁₂‖≤𝜇₀ for some sufficiently small real 𝜇₀>0 with ‖𝑀‖ denoting the norm of the matrix 𝑀.
(iii) The nominal pair (𝐴₀,𝐶₀) is observable.

Remark 2.1. (i) Given any state-space realization, there always exists an infinite number of state-space realizations in triangular form as (2.2) which have the same transfer function that the former has. Each one of such state-space realizations may be obtained by applying an appropriate coordinates transformation on the original realization. Then, whenever the original state-space realization of the plant is not a triangular form, an appropriate coordinates transformation is required to obtain a triangular form state-space realization as (2.2).
(ii) The internal representation (2.2) gives place to an input-output relation defined by a transfer function as𝐺_𝑝(𝑠)=𝐺₀(𝑠)(1+Δ_𝑚(𝑠))+Δ_𝑎(𝑠), (2.3) where 𝐺₀(𝑠) denotes the transfer function of the plant nominal model and Δ_𝑚(𝑠) and Δ_𝑎(𝑠) are two rational transfer functions related to the multiplicative and additive unmodeled dynamics, respectively. The poles of Δ_𝑚(𝑠) and Δ_𝑎(𝑠) are the eigenvalues of the block 𝐴₂₂ and their gains depend on the norms of the blocks 𝐵₂₁ and 𝐶₁₂.
(iii) If the nominal triple {𝐴₀,𝐵₀,𝐶₀} is known, then a classical pole-placement controller may be synthesized without using estimation. However, this knowledge is not necessary to synthesize adaptive control with the less restrictive knowledge of 𝑛₀, which is the nominal plant order.

The plant can be unstable and of nonstable inverse. Then, the use of the model-matching technique, with a free-chosen reference model, for the controller synthesis can be used with a discrete-time controller. In such a case, the reconstruction process of the continuous-time plant input from the discrete-time control output gives the possibility of obtaining an inversely stable discretized plant model. Such a reconstruction has to be made with the use of a hold device, a FROH in the most general case, combined with a multirate with fast input rate. This multirate provides free-design parameters, related to multirate gains, which can be adjusted so that the discretized plant model could be of stable inverse. In this way, the model-matching technique can be used to synthesize a discrete-time controller to stabilize the continuous-time plant while matching a freely chosen reference model at sampling instants. The plant input obtained from such a reconstruction method is given by𝑢(𝑡)=𝛼_𝑗{𝑢(𝑘)+𝛽𝑢(𝑘)−𝑢(𝑘−1)𝑇(𝑡−𝑘𝑇)} (2.4) for 𝑡∈[𝑘𝑇+(𝑗−1)𝑇^′,𝑘𝑇+𝑗𝑇^′), 𝑗∈{1,2,…,𝑁}, where 𝛽∈[−1,1]∩ℝ is the FROH correcting gain, T is the sampling period for the state and output (slow sampling) which is uniformly divided in 𝑁 subperiods of length 𝑇^′=𝑇/𝑁 (fast sampling) to generate the fast sampling plant input, 𝑢(𝑘) denotes the value of the controller output sequence at the instant kT, for all 𝑘∈ℤ⁺₀≜ℤ⁺∪{0}, and 𝛼_𝑗∈ℝ are the multirate gains. The technique of using 𝑇≠𝑇^′<𝑇, with 𝑁 exceeding some prescribed lower bound, is the key feature for always achieving a stable discrete-time transfer function numerator even if the continuous-time plant transfer function 𝐺_𝑝(𝑠) has some critically stable or unstable zero. In this sense, the FROH device operates on the sequence {𝑢(𝑘)} defined at the slow sampling instants kT and then the input 𝑢(𝑡) is generated over each subperiod 𝑇^′ with the corresponding gain 𝛼_𝑗. Such gains have to be suitably selected to ensure the stability of the zeros of the discretized plant model which relates the sequences {𝑢(𝑘)} and {𝑦(𝑘)} (plant output sequence) defined over the sampling period 𝑇.

The state-space representation corresponding to the discrete-time plant obtained from the discretization of the continuous-time plant by applying the FROH with the multirate is given by𝑥₀(𝑘+1)=𝐹(𝑇)𝑥₀(𝑘)+𝐻₁(𝑇)𝑢(𝑘)+𝐻₂(𝑇)𝑢(𝑘−1),𝑦(𝑘)=𝐶₀𝑥₀(𝑘)+𝜂(𝑘), (2.5) where 𝜂(𝑘) denotes the contribution of the unmodeled dynamics to the discretized plant output, 𝐹(𝑇)=𝜓^𝑁(𝑇)=𝜙(𝑇)=𝑒^𝐴0𝑇∈ℝ^𝑛0×𝑛0 is the state transition matrix of the continuous-time nominal dynamic system valued during a sampling period, and𝐻₁(𝑇)=^𝑁_ℓ=1𝛼_ℓ𝜓^𝑁−ℓ(𝑇)[(1+ℓ−1𝑁𝛽)Γ(𝑇^′𝛽)+𝑇Γ^′(𝑇^′)]=𝐶_Δ(𝑇)𝑔∈ℝ^𝑛0×1,𝐻₂(𝑇)=−𝛽^𝑁_ℓ=1𝛼_ℓ𝜓^𝑁−ℓ(𝑇)[ℓ−1𝑁Γ(𝑇^′1)+𝑇Γ^′(𝑇^′)]=−𝛽𝐶^′_Δ(𝑇)𝑔∈ℝ^𝑛0×1, (2.6) withΓ(𝑇^′)=^𝑇′₀𝜙(𝑇^′−𝑠)𝐵₀𝑑𝑠∈ℝ^𝑛0×1,𝐶_Δ(𝑇)=[𝜓^𝑁−1(𝑇)Δ₁(𝑇)⋯𝜓(𝑇)Δ_𝑁−1(𝑇)Δ_𝑁(𝑇)]∈ℝ^𝑛0×𝑁,Γ^′(𝑇^′)=^𝑇′₀𝜙(𝑇^′−𝑠)𝐵₀𝑠𝑑𝑠∈ℝ^𝑛0×1,𝐶^′_Δ(𝑇)=[𝜓^𝑁−1(𝑇)Δ^′₁(𝑇)⋯𝜓(𝑇)Δ^′_𝑁−1(𝑇)Δ^′_𝑁(𝑇)]∈ℝ^𝑛0×𝑁,Δ_𝑗(𝑇)=(1+𝑗−1𝑁𝛽)Γ(𝑇^′𝛽)+𝑇Γ^′(𝑇^′)∈ℝ^𝑛0×1,Δ^′_𝑗(𝑇)=𝑗−1𝑁Γ(𝑇^′1)+𝑇Γ^′(𝑇^′)∈ℝ^𝑛0×1,𝑔=[𝛼₁⋯𝛼_𝑁]^𝑇. (2.7)

2.1. Input-Output Relation for the Discretized Plant

From the output equation of (2.5), substituting the state-space equation and iterating 𝑛₀ times, it follows that𝑦(𝑘)=𝐶₀{𝐹^𝑛0𝑥₀(𝑘−𝑛₀)+𝐻₁𝑢(𝑘−1)+^𝑛0−1_𝑖=1𝐹^𝑖−1(𝐹𝐻₁+𝐻₂)𝑢(𝑘−𝑖−1)+𝐹^𝑛0−1𝐻₂𝑢(𝑘−𝑛₀−1)}+𝜂(𝑘), (2.8) where the argument 𝑇 in 𝐹(𝑇), 𝐻₁(𝑇), and 𝐻₂(𝑇) has been omitted for simplicity. In a similar way, one can obtain that𝑦(𝑘−ℓ)=𝐶₀{𝐹^𝑛0−ℓ𝑥₀(𝑘−𝑛₀)+𝐻₁𝑢(𝑘−ℓ−1)+^𝑛0−1_𝑖=ℓ+1𝐹^{𝑖−ℓ−1}(𝐹𝐻₁+𝐻₂)𝑢(𝑘−𝑖−1)+𝐹^{𝑛0−ℓ−1}𝐻₂𝑢(𝑘−𝑛₀−1)}+𝜂(𝑘−ℓ) (2.9) for ℓ∈{1,2,…,𝑛₀−1}. (2.9) together with 𝑦(𝑘−𝑛₀)=𝐶₀𝑥₀(𝑘−𝑛₀)+𝜂(𝑘−𝑛₀) can be rewritten as𝑌_𝑣(𝑘−1)=𝑉𝑥₀(𝑘−𝑛₀)+Π_𝑣(𝑘−1)+Φ𝑈_𝑣(𝑘−2), (2.10) where𝑌_𝑣(𝑘−1)=[𝑦(𝑘−𝑛₀)𝑦(𝑘−𝑛₀+1)⋯𝑦(𝑘−2)𝑦(𝑘−1)]^𝑇,𝑉=𝐶^𝑇₀(𝐶₀𝐹)^𝑇⋯(𝐶₀𝐹^𝑛0−2)^𝑇(𝐶₀𝐹^𝑛0−1)^𝑇𝑇,Π_𝑣(𝑘−1)=[𝜂(𝑘−𝑛₀)𝜂(𝑘−𝑛₀+1)⋯𝜂(𝑘−2)𝜂(𝑘−1)]^𝑇,𝑈_𝑣(𝑘−2)=[𝑢(𝑘−𝑛₀−1)𝑢(𝑘−𝑛₀)⋯𝑢(𝑘−3)𝑢(𝑘−2)]^𝑇,𝐶Φ=00⋯00₀𝐻₂𝐶₀𝐻₁𝐶⋯00₀𝐹𝐻₂𝐶₀(𝐹𝐻₁+𝐻₂𝐶)⋯00⋮⋮⋱⋮⋮₀𝐹^𝑛0−3𝐻₂𝐶₀𝐹^𝑛0−4(𝐹𝐻₁+𝐻₂)⋯𝐶₀𝐻₁0𝐶₀𝐹^𝑛0−2𝐻₂𝐶₀𝐹^𝑛0−3(𝐹𝐻₁+𝐻₂)⋯𝐶₀(𝐹𝐻₁+𝐻₂)𝐶₀𝐻₁ (2.11) with 𝑉 being the observability matrix for the discretized plant nominal model. By substituting the expression for 𝑥₀(𝑘−𝑛₀), obtained from (2.10), in (2.8) it follows that𝑦(𝑘)=−^𝑛0_𝑖=1𝑎_𝑖𝑦(𝑘−𝑖)+^𝑛𝑑_𝑖=1𝑏_𝑖𝑢(𝑘−𝑖)+Ω(𝑘)=−^𝑛0_𝑖=1𝑎_𝑖𝑦(𝑘−𝑖)+^𝑛𝑑𝑁_𝑖=1𝑗=1𝛼_𝑗𝑏_𝑖,𝑗𝑢(𝑘−𝑖)+Ω(𝑘) (2.12) where 𝑛_𝑑=𝑛₀ if 𝛽=0 or 𝑛_𝑑=𝑛₀+1 if 𝛽≠0, and𝑎_𝑖=−[𝐶₀𝐹^𝑛0𝑉⁻¹]_{𝑛0−𝑖+1},𝑏_𝑖=^𝑁_𝑗=1𝛼_𝑗𝑏_𝑖,𝑗=𝐶₀𝐹^𝑖−2(𝐹𝐻₁+𝐻₂)−[𝐶₀𝐹^𝑛0𝑉⁻¹Φ]_{𝑛0−𝑖+2}for𝑖∈{2,3,…,𝑛₀𝑏},₁=^𝑁_𝑗=1𝛼_𝑗𝑏_1,𝑗=𝐶₀𝐻₁,𝑏_𝑛0+1=^𝑁_𝑗=1𝛼_𝑗𝑏_𝑛0+1,𝑗=𝐶₀𝐹^𝑛0−1𝐻₂−[𝐶₀𝐹^𝑛0𝑉⁻¹Φ]₁,Ω(𝑘)=𝜂(𝑘)−^𝑛0_𝑖=1[𝐶₀𝐹^𝑛0𝑉⁻¹]_{𝑛0−𝑖+1}𝜂(𝑘−𝑖) (2.13) with [𝑣]_𝑖 denoting the 𝑖th component of the vector 𝑣 and having into account the expressions (2.6) for 𝐻₁(𝑇) and 𝐻₂(𝑇). Note that Ω(𝑘) contains the contribution of the unmodeled dynamics to the discretized plant model output.

Remark 2.2. (i) The observability of the pair (𝐴₀,𝐶₀) implies the nonsingularity of 𝑉 whenever 𝑇≥𝑇₀>0, for some sufficiently small real 𝑇₀, and conversely. Note that 𝑉 tends to be singular as the sampling period 𝑇 tends to zero since 𝐹(𝑇) would tend to the identity matrix.
(ii) The modeled part of the discretized plant (2.12) can be expressed equivalently as the discrete transfer function𝐺_𝑑𝐵(𝑧)=^𝑑(𝑧)𝐴_𝑑(𝑧), (2.14) where 𝐴_𝑑(𝑧)=𝑧^𝑛0+1+∑^𝑛0_𝑖=1𝑎_𝑖𝑧^{𝑛0−𝑖+1} and 𝐵_𝑑(𝑧)=∑^𝑛0+1_𝑖=1𝑏_𝑖𝑧^{𝑛0−𝑖+1} with 𝑧 being the Z-transform argument used in discrete-time transfer functions [23]. In the particular case that 𝛽=0 it follows from (2.6) that 𝐻₂(𝑇)=0, then 𝑏_𝑛0+1=0 since the first column of Φ is zero and then the transfer function 𝐺_𝑑(𝑧) possesses a zero-pole cancellation at 𝑧=0; that is, its order is 𝑛₀ instead of 𝑛₀+1. In the rest of the paper, 𝑛_𝑑∈ℤ⁺ is used to denote the order of such a transfer function being 𝑛_𝑑=𝑛₀ if 𝛽=0 (i.e., a ZOH is used in the discretization process) or 𝑛_𝑑=𝑛₀+1 if 𝛽≠0. The parameter 𝑛_𝑑 allows a unified definition of polynomials and discrete transfer functions for different degrees associated with 𝛽=0 and 𝛽≠0.

Note that the coefficients 𝑏_𝑖, for 𝑖∈{1,2,…,𝑛_𝑑}, of the polynomial 𝐵_𝑑(𝑧) depend on the multirate gains 𝛼_𝑗, for 𝑗∈{1,…,𝑁}, included in 𝐻₁(𝑇) and 𝐻₂(𝑇); that is,𝑣=𝑀𝑔, (2.15) where 𝑣=[𝑏₁𝑏₂⋯𝑏_𝑛𝑑]^𝑇 and 𝑀=[𝑏_𝑖,𝑗]∈ℝ^{𝑛𝑑×𝑁}. The components 𝑏_𝑖,𝑗 depend on the sampling period 𝑇, the correcting gain 𝛽∈[−1,1] of the FROH, and the matrices 𝐴₀,𝐵₀ and 𝐶₀ which define the plant nominal model. In this context, if the multirate gain vector is suitably chosen, one may fix the coefficients 𝑏_𝑖 at desired values and, in this way, one may place the zeros of the discretized plant nominal model at desired locations, namely, within the stability domain. This is the strategy to get the stabilization of the closed-loop system by means of a model-matching controller.

Assumption 2. The correcting gain 𝛽 of the FROH and the sampling period 𝑇 are chosen such that 𝑀 is a full-rank matrix.

Remark 2.3. In the non-adaptive case (known plant parameters), Assumption 2 is crucial to calculate the multirate parameterization 𝑔 from (2.15) provided that 𝑁≥𝑛_𝑑. If 𝑁=𝑛_𝑑, 𝑔=𝑀⁻¹𝑣 is the unique solution for the multirate gains which places the discretized plant zeros at prefixed locations, those linked to the roots of the polynomial whose coefficients are the components of 𝑣. In this way, if the zeros of such a polynomial are within the stability domain, then the discretized plant transfer function is inversely stable. On the contrary, if 𝑁>𝑛_𝑑, different solutions can be obtained for 𝑔. However, in the adaptive case, such an assumption can be relaxed if the parameters estimation algorithm guarantees that the rank of the matrix 𝑀(𝑘), composed with the estimates of 𝑏_𝑖,𝑗 (namely, 𝑏̂_𝑖,𝑗), is 𝑛_𝑑 for all 𝑘∈ℤ⁺₀.

The discretized plant model (2.12) can be written as𝑦(𝑘)=𝜃^𝑇_𝑎𝜑_𝑦(𝑘−1)+^𝑛𝑑_𝑖=1𝜃^𝑇_𝑏,𝑖𝑢(𝑘−𝑖)+Ω(𝑘)=𝜃^𝑇𝜑(𝑘−1)+Ω(𝑘), (2.16) where𝜃=[𝜃^𝑇_𝑎𝜃^𝑇_𝑏,1𝜃^𝑇_𝑏,2⋯𝜃^𝑇_{𝑏,𝑛𝑑}]^𝑇,𝜃_𝑎=[−𝑎₁−𝑎₂⋯−𝑎_𝑛0]^𝑇,𝜃_𝑏,𝑖=[𝑏_𝑖,1𝑏_𝑖,2⋯𝑏_𝑖,𝑁]^𝑇,𝜑(𝑘−1)=[𝜑^𝑇_𝑦(𝑘−1)𝑢^𝑇(𝑘−1)𝑢^𝑇(𝑘−2)⋯𝑢^𝑇(𝑘−𝑛_𝑑)]^𝑇,𝜑_𝑦(𝑘−1)=[𝑦(𝑘−1)𝑦(𝑘−2)⋯𝑦(𝑘−𝑛₀)]^𝑇,𝑢(𝑘−𝑖)=[𝛼₁𝑢(𝑘−𝑖)𝛼₂𝑢(𝑘−𝑖)⋯𝛼_𝑁𝑢(𝑘−𝑖)]^𝑇 (2.17) for 𝑖∈{1,2,…,𝑛_𝑑} and 𝑗∈{1,…,𝑁}. In the following, the case 𝑁=𝑛_𝑑 is considered.

3. Adaptive Control Design

The control objective is the adaptive stabilization of the continuous-time plant while matching, with a bounded tracking-error, a stable discrete-time free-design reference model 𝐺_𝑚(𝑧)=𝐵_𝑚(𝑧)/𝐴_𝑚(𝑧) at the sampling instants. The perfect tracking is not achievable due to the presence of unknown unmodeled dynamics. A self-tuning regulator scheme is used to meet the control objective [2, 15]. The control law structure is firstly presented for the non-adaptive case, that is, when the plant to be controlled is known. Then, an extension to the adaptive case is developed, which is the main interest of the paper. In such a case, a recursive algorithm of least-squares type with an adaptation dead-zone is used to obtain an estimation of the unknown parameters included in the vector 𝜃 of (2.17) at each sampling instant. Then, the multirate gains are updated in order to guarantee the inverse stability of the transfer-like function associated to the discretized plant estimated model. Such a model is used to parameterize the adaptive controller.

3.1. Known Plant

The proposed control law is obtained from𝑅(𝑧)𝑢(𝑘)=𝑇(𝑧)𝑐(𝑘)−𝑆(𝑧)𝑦(𝑘) (3.1) for all 𝑘∈ℤ⁺₀ where {𝑐(𝑘)} is the reference input sequence. The reconstruction of the plant input 𝑢(𝑡) is made by using (2.4), with the control sequence {𝑢(𝑘)} obtained from (3.1) and the multirate gains 𝛼_𝑗, for all 𝑗∈{1,…,𝑁}, obtained from (2.15) with an appropriate choice of 𝑣 to guarantee the inverse stability of the discretized plant nominal model; that is, such gains fix the numerator 𝐵_𝑑(𝑧) of (2.14) to a prefixed one 𝐵^′(𝑧), whose coefficients are the components of 𝑣, with the roots within the stability domain. An important mathematical issue in the current context is that the proposed method allows the stabilization of the numerator polynomial of the discretized transfer function by an appropriate choice of the multirate gains even if its continuous-time counterpart is unstable or critically stable.

The discrete-time transfer function of the closed-loop system obtained from the application of the control law (3.1) to the discretized plant model of transfer function (2.14) is given by𝑌(𝑧)=𝐵𝐶(𝑧)^′(𝑧)𝑇(𝑧)𝐴(𝑧)𝑅(𝑧)+𝐵^′=(𝑧)𝑆(𝑧)𝑇(𝑧)𝐴(𝑧)+𝑆(𝑧) (3.2) if 𝑅(𝑧)=𝐵^′(𝑧) is chosen. In this way, the polynomial 𝐵^′(𝑧) is cancelled in the closed-loop system so that it does not generate discrete plant zeros. It should be stable (i.e., with zeros in |𝑧|<1) to cancel it if the usual methods without multirate techniques are used, [12, 14, 22]. Otherwise, it cannot be cancelled by the controller and it has then to be transmitted as a factor of the numerator of the closed-loop discrete transfer function (3.2). If it is transmitted, then the reference model is not of complete free choice since it has to contain this polynomial factor. The method proposed in this manuscript allows the stabilization of the discretized plant zeros. As a result, the reference model transfer function is always freely chosen with no zero transmission constraints by using the multirate technique with appropriate choice of the multirate gains. The control polynomials 𝑇(𝑧) and 𝑆(𝑧) to meet the model-matching objective are obtained from𝑇(𝑧)=𝐵_𝑚(𝑧)𝐴_𝑠(𝑧)𝐴(𝑧)+𝑆(𝑧)=𝐴_𝑚(𝑧)𝐴_𝑠(𝑧) (3.3) with the following degree constraints, required for controller realizability, in the controller synthesis:Deg[𝐴_𝑚(𝑧)]=Deg[𝐴(𝑧)]−Deg[𝐴_𝑠(𝑧)],Deg[𝑆(𝑧)]=Deg[𝐴(𝑧)]−1=𝑛_𝑑−1=𝑁−1,Deg[𝑇(𝑧)]=Deg[𝐵_𝑚(𝑧)]+Deg[𝐴_𝑠(𝑧)]≤𝑁−1, (3.4) where 𝐴_𝑠(𝑧) is a stable monic polynomial of zero-pole cancellations of the closed-loop system. In this way, the nominal closed-loop system matches the reference model at the sampling instants, but not the true closed-loop system due to the presence of unmodeled dynamics. However, the tracking-error is guaranteed to be bounded at all sampling times subject to Assumption 1(ii).

3.2. Unknown Plant

If the continuous-time plant parameters are unknown, then the vector 𝜃 in (2.17), which is composed of the discretized plant model parameters, is also unknown. However, all of the above control design in Section 3.1 remains valid if such a parameter vector is updated by an estimation algorithm. Such an algorithm provides an adaptation of each parameter 𝑏_𝑖,𝑗, namely, 𝑏̂_𝑖,𝑗(𝑘), for 𝑖,𝑗∈{1,…,𝑁} and all 𝑘∈ℤ⁺₀. Then, the multirate gains 𝛼_𝑗, now their estimates being 𝛼_𝑗(𝑘), are calculated from an equation similar to (2.15) by replacing the matrix 𝑀 by its estimated 𝑀(𝑘). In this way, the numerator of the corresponding discretized plant estimated model is fixed to 𝐵^′(𝑧). Note that such a numerator is time invariant although the estimates 𝜃_𝑏,𝑖(𝑘), for 𝑖∈{1,…,𝑁}, are time varying. Furthermore, the controller parameterization can be obtained from 𝑅(𝑧)=𝐵^′(𝑧) and equations similar to (3.3) by replacing the discretized plant polynomial 𝐴(𝑧) by its corresponding estimated one at the current sampling instant, namely, 𝐴(𝑧,𝑘) [2]. In this context, the polynomials 𝑇(𝑧) and 𝑅(𝑧) have to be calculated once for all since 𝐵_𝑚(𝑧), 𝐴_𝑠(𝑧) and 𝐵^′(𝑧), are time invariant. On the contrary, 𝑆(𝑧), now being 𝑆(𝑧,𝑘), is updated at each running sampling instant since the polynomial 𝐴(𝑧,𝑘) is time varying.

The time-varying multirate gains 𝛼_𝑗(𝑘) are used to calculate the plant input within the inter sample period via (2.4). Then, a time-varying input-output relation for the discretized plant is derived by following similar steps to those in Section 2.1. In this sense, one obtains𝑦(𝑘−ℓ)=𝐶₀{𝐹^𝑛0−ℓ𝑥₀(𝑘−𝑛₀)+𝐻₁(𝑘−ℓ−1)𝑢(𝑘−ℓ−1)+^𝑛0−1_𝑖=ℓ+1𝐹^{𝑖−ℓ−1}(𝐹𝐻₁(𝑘−𝑖−1)+𝐻₂(𝑘−𝑖))×𝑢(𝑘−𝑖−1)+𝐹^{𝑛0−ℓ−1}𝐻₂(𝑘−𝑛₀)𝑢(𝑘−𝑛₀−1)}+𝜂(𝑘−ℓ) (3.5) for ℓ∈{0,1,…,𝑛₀−1}, where 𝐻₁(𝑖) and 𝐻₂(𝑖) are obtained from equations similar to those in (2.6) by replacing 𝛼_𝑗 by 𝛼_𝑗(𝑖), for 𝑖∈{𝑘−1,…,𝑘−𝑛₀} and 𝑗∈{1,…,𝑁}. Then𝑌_𝑣(𝑘−1)=𝑉𝑥₀(𝑘−𝑛₀)+Π_𝑣(𝑘−1)+Φ(𝑘−2,…,𝑘−𝑛₀)𝑈_𝑣(𝑘−2), (3.6) where 𝑌_𝑣(𝑘−1), 𝑉, Π_𝑣(𝑘−1), and 𝑈_𝑣(𝑘−2) are defined in (2.11) whileΦ(𝑘−2,…,𝑘−𝑛₀Φ)=00⋯00_2,1Φ_2,2Φ⋯00_3,1Φ_3,2Φ⋯00⋮⋮⋱⋮⋮_𝑛0−1,1Φ_𝑛0−1,2⋯Φ_{𝑛0−1,𝑛0−1}0Φ_𝑛0,1Φ_𝑛0,2⋯Φ_{𝑛0,𝑛0−1}Φ_𝑛0,𝑛0 (3.7) with Φ_𝑖,𝑗=𝐶₀𝐹^{𝑖−𝑗−1}(𝐹𝐻₁(𝑘−𝑛₀+𝑗−2)+𝐻₂(𝑘−𝑛₀+𝑗−1)), Φ_𝑖,𝑖=𝐶₀𝐻₁(𝑘−𝑛₀+𝑖−2), and Φ_𝑖,1=𝐶₀𝐹^𝑖−2𝐻₂(𝑘−𝑛₀) for 𝑖,𝑗∈{2,3,…,𝑛₀} and 𝑖>𝑗. By substituting in (3.5) the expression for 𝑥₀(𝑘−𝑛₀), obtained from (3.6), for ℓ=0 it follows that𝑦(𝑘)=−^𝑛0_𝑖=1𝑎_𝑖𝑦(𝑘−𝑖)+𝐶₀𝐻₁+(𝑘−1)𝑢(𝑘−1)_{𝑛0−1𝑖=1}{𝐶₀𝐹^𝑖−1(𝐹𝐻₁(𝑘−𝑖−1)+𝐻₂(𝑘−𝑖))−[𝐶₀𝐹^𝑛0𝑉⁻¹Φ(𝑘−2,…,𝑘−𝑛₀)]_{𝑛0−𝑖+1}}×𝑢(𝑘−𝑖−1)+(𝐶₀𝐹^𝑛0−1𝐻₂(𝑘−𝑛₀)−[𝐶₀𝐹^𝑛0𝑉⁻¹Φ(𝑘−2,…,𝑘−𝑛₀)]₁)𝑢(𝑘−𝑛₀+−1)_𝑛0𝑖=1𝑎_𝑖𝜂(𝑘−𝑖)+𝜂(𝑘), (3.8) or𝑦(𝑘)=−^𝑛0_𝑖=1𝑎_𝑖𝑦(𝑘−𝑖)+𝐶₀𝐻₁+(𝑘−1)𝑢(𝑘−1)_{𝑛0−1𝑖=1}{𝐶₀𝐹^𝑖−1(𝐹𝐻₁(𝑘−1)+𝐻₂(𝑘−1))−[𝐶₀𝐹^𝑛0𝑉⁻¹Φ(𝑘−1)]_{𝑛0−𝑖+1}}𝑢(𝑘−𝑖−1)+(𝐶₀𝐹^𝑛0−1𝐻₂(𝑘−1)−[𝐶₀𝐹^𝑛0𝑉⁻¹Φ(𝑘−1)]₁)𝑢(𝑘−𝑛₀−1)+^𝑛0_𝑖=1𝑎_𝑖𝜂(𝑘−𝑖)+𝜂(𝑘)+Λ_𝛼(𝑘−1), (3.9) whereΛ_𝛼=(𝑘−1)_{𝑛0−1𝑖=1}{𝐶₀𝐹^𝑖−1(𝐹Δ𝐻₁(𝑘−𝑖−1)+Δ𝐻₂(𝑘−𝑖))−[𝐶₀𝐹^𝑛0𝑉⁻¹ΔΦ]_{𝑛0−𝑖+1}}𝑢(𝑘−𝑖−1)+(𝐶₀𝐹^𝑛0−1Δ𝐻₂(𝑘−𝑛₀)−[𝐶₀𝐹^𝑛0𝑉⁻¹ΔΦ]₁)𝑢(𝑘−𝑛₀−1) (3.10) with Δ𝐻_ℓ(𝑘−𝑖−1)=𝐻_ℓ(𝑘−𝑖−1)−𝐻_ℓ(𝑘−1), for ℓ=1,2, and ΔΦ=Φ(𝑘−2,…,𝑘−𝑛₀)−Φ(𝑘−1). Note that Ω(𝑘)=∑^𝑛0_𝑖=1𝑎_𝑖𝜂(𝑘−𝑖)+𝜂(𝑘) arises from the unmodeled dynamics of the continuous-time plant while Λ_𝛼(𝑘−1) comes from the fact that the multirate gains are time varying. Both terms constitute the unmodeled dynamics for the discretized plant model. On one hand, the latter tends to zero if the multirate gains converge to constant values. On the other hand, the former is such that|Ω(𝑘)|≤Ω(𝑘)=𝜇₁𝜌(𝑘)+𝜇₂ (3.11) for some known real constants 𝜇₁≥0 and 𝜇₂≥0, where𝜌(𝑘)=Sup_{0≤𝑘′≤𝑘}{|𝑤^𝑇𝑥(𝑘^′)|𝜎^{𝑘−𝑘′}} (3.12) for all 𝑘∈ℤ⁺₀, some known constant vector 𝑤, and some known real constant 𝜎∈(0,1) with𝑥(𝑘)=[𝑦(𝑘−1)⋯𝑦(𝑘−𝑛₀)𝑢(𝑘−1)⋯𝑢(𝑘−𝑁)]^𝑇 (3.13) in view of Assumption 1(ii) [14]. Finally, one can express the discretized plant model as𝑦(𝑘)=−^𝑛0_𝑖=1𝑎_𝑖𝑦(𝑘−𝑖)+^𝑛𝑑_𝑖=1𝑏_𝑖(𝑘−1)𝑢(𝑘−𝑖)+Ω(𝑘)+Λ_𝛼(𝑘−1)=−^𝑛0_𝑖=1𝑎_𝑖𝑦(𝑘−𝑖)+^𝑛𝑑𝑁_𝑖=1𝑗=1𝑏_𝑖,𝑗𝛼_𝑗(𝑘−1)𝑢(𝑘−𝑖)+Ω(𝑘)+Λ_𝛼(𝑘−1)=𝜃^𝑇𝜑(𝑘−1)+Ω(𝑘)+Λ_𝛼(𝑘−1), (3.14) where 𝜑(𝑘−1) is built like 𝜑(𝑘−1) in (2.17) by replacing the components of 𝑢(𝑘−𝑖), namely, 𝛼_𝑗𝑢(𝑘−𝑖), by the corresponding 𝛼_𝑗(𝑘−1)𝑢(𝑘−𝑖), for 𝑖,𝑗∈{1,…,𝑁}.

The algorithm used to obtain an on-line adaptation 𝜃(𝑘) of the unknown parameters vector 𝜃 is described below.

3.2.1. Estimation Algorithm

An “a priori” estimated parameters vector is obtained by using a recursive least-squares algorithm defined by𝜃⁰(𝑘)=𝜃⁰𝜑(𝑘−1)+𝑠(𝑘)𝑃(𝑘−1)_𝑛(𝑘−1)𝑒⁰_𝑛(𝑘)𝛾(𝑘)+𝑠(𝑘)𝜑^𝑇_𝑛(𝑘−1)𝑃(𝑘−1)𝜑_𝑛,𝜑(𝑘−1)𝑃(𝑘)=𝑃(𝑘−1)−𝑠(𝑘)𝑃(𝑘−1)_𝑛(𝑘−1)𝜑^𝑇_𝑛(𝑘−1)𝑃(𝑘−1)𝛾(𝑘)+𝑠(𝑘)𝜑^𝑇_𝑛(𝑘−1)𝑃(𝑘−1)𝜑_𝑛(𝑘−1) (3.15) for all 𝑘∈ℤ⁺, where 𝜑_𝑛(𝑘−1)=𝜑(𝑘−1)/(1+‖𝜑(𝑘−1)‖), 𝑃(𝑘−1) is initialized such that 𝑃(0)=𝑃^𝑇(0)>0 (denoting positive definiteness), 𝛾(𝑘)>0, 𝑒⁰_𝑛(𝑘)=𝑒⁰(𝑘)/(1+‖𝜑(𝑘−1)‖) with𝑒⁰(𝑘)=−𝜃^0𝑇(𝑘−1)𝜑(𝑘−1)+Ω(𝑘)+Λ_𝛼(𝑘−1) (3.16) being the “a priori” estimation error as well as 𝜃⁰(𝑘−1)=𝜃⁰(𝑘−1)−𝜃 the “a priori” parametrical error, and 𝑠(𝑘) is a relative adaptation dead-zone defined as𝑠(𝑘)=0if𝑒⁰_𝑎𝑛(𝑘)≤𝜇𝜂_𝑎𝑛(𝑘),𝑓(𝑘)𝑒⁰_𝑎𝑛(𝑘)otherwise (3.17) for some 𝜇>1, where 𝑒⁰_𝑎𝑛(𝑘)=((𝑒⁰_𝑛(𝑘))²+𝜑^𝑇_𝑛(𝑘−1)𝑃²(𝑘−1)𝜑_𝑛(𝑘−1))^1/2 is an augmented normalized error, 𝜂_𝑎𝑛(𝑘)=(1+𝛾⁻¹(𝑘)𝜑^𝑇_𝑛(𝑘−1)𝑃(𝑘−1)𝜑_𝑛(𝑘−1))^1/2𝜂_𝑇𝑛(𝑘) with 𝜂_𝑇𝑛(𝑘)=𝜂_𝑇(𝑘)/(1+‖𝜑(𝑘−1)‖) and 𝜂_𝑇(𝑘) being an upper bound for |𝜂_𝑇(𝑘)|=|Ω(𝑘)+Λ_𝛼(𝑘−1)|, and𝑓(𝑘)=0if𝑒⁰_𝑎𝑛(𝑘)≤𝜇𝜂_𝑎𝑛𝑒(𝑘),_0𝑎𝑛(𝑘)−𝜇𝜂_𝑎𝑛(𝑘)otherwise. (3.18) This algorithm provides an estimation 𝜃⁰(𝑘) of the parameters vector. Then, an “a posteriori” estimates vector is obtained as follows.

Estimates Modification
This algorithm consists of three main steps as follows.Step 1. Build the matrix 𝑀⁰(𝑘)=[𝑏̂⁰_𝑖,𝑗(𝑘)]∈ℝ^𝑁×𝑁, for 𝑖,𝑗∈{1,2,…,𝑁}, from the “a priori” estimates 𝜃⁰_𝑏,𝑖(𝑘), included in 𝜃⁰(𝑘), of the corresponding 𝜃_𝑏,𝑖 defined in (2.17).Step 2. 𝑀(𝑘)=𝑀⁰(𝑘)∶If  |Det(𝑀(𝑘))|≥𝛿₀  then  𝜃_𝑏,𝑖(𝑘)=𝜃⁰_𝑏,𝑖(𝑘)         else while  |Det(𝑀(𝑘))|<𝛿₀                𝑀(𝑘)=𝑀(𝑘)+𝛿𝐼_𝑁           end;           for  𝑖=1  to  𝑁                𝜃_𝑏,𝑖(𝑘)=𝑀_𝑖(𝑘)           end,end.Step 3. 𝜃(𝑘)=[𝜃_𝑎^0𝑇(𝑘)𝜃^𝑇_𝑏,1(𝑘)𝜃^𝑇_𝑏,2(𝑘)⋯𝜃^𝑇_𝑏,𝑁(𝑘)]^𝑇,for some positive real constants 𝛿≪1 and 𝛿₀≪1, and where 𝑀_𝑖(𝑘) denotes the 𝑖-th row of 𝑀(𝑘) and 𝐼_𝑁 the identity matrix of dimension 𝑁×𝑁.

Once the estimated parameters are updated, the multirate gains vector 𝑔(𝑘)=[𝛼₁(𝑘)⋯𝛼_𝑁(𝑘)]^𝑇 is calculated from an equation similar to (2.15) by replacing the matrix 𝑀 by 𝑀(𝑘)=[𝑏̂_𝑖,𝑗(𝑘)]∈ℝ^𝑁×𝑁.

Remark 3.1. (i) The proposed estimates modification process avoids that the time-varying matrix 𝑀(𝑘) be close to a singular one. In this way, a bounded vector of multirate gains is obtained at all sampling instants.
(ii) Note that the estimate 𝜃⁰_𝑎(𝑘) corresponding to the parameters of 𝜃_𝑎 is not affected by the modification algorithm. In fact, such a modification only affects the entries in the main diagonal of 𝑀⁰(𝑘). Also, note that the instruction while of the second step is executed a finite number of times since there exists a finite integer number ℓ such that |Det[𝑀(𝑘)]|=|Det[𝑀⁰(𝑘)+ℓ𝛿𝐼_𝑁]|=|(ℓ𝛿)^𝑁+𝑓₀(𝛿,𝜃⁰_𝑏,𝑖(𝑘))|≥𝛿₀ for 𝑖∈{1,…,𝑁} and some function 𝑓₀(⋅,⋅).
(iii) From (3.10) to (3.13) and the construction of 𝜑(𝑘−1), there exist some real constants 𝜐_𝑖≥0, for 𝑖∈{1,2,3,4} such that the sequences {𝜂₁(𝑘)}, with 𝜂₁(𝑘)=𝜐₁‖𝜑(𝑘−1)‖+𝜐₂, and {𝜂₂(𝑘)}, with 𝜂₂(𝑘)=𝜐₃‖𝑥(𝑘−1)‖+𝜐₄, are both upper bounds for {|𝜂_𝑇(𝑘)|}.

The estimation algorithm together with such a multirate gains adaptation possesses the following properties.

Lemma 3.2 (main properties of the estimation algorithm). (i) 𝑃(𝑘) is uniformly bounded for all 𝑘∈ℤ⁺₀, and it asymptotically converges to a finite, at least semidefinite positive, limit as 𝑘→∞.
(ii) 𝜃⁰(𝑘) and 𝜃(𝑘) are uniformly bounded for all samples and converge to finite limits.
(iii) 𝑓(𝑘)<∞ for all 𝑘∈ℤ⁺₀ and lim_𝑘→∞𝑓(𝑘)=0.
(iv) The vector 𝑔(𝑘) is bounded for all samples and converges to a finite limit.

The proof is made in Appendix A.

4. Stability Analysis

The plant discretized model can be written as follows:𝑦(𝑘)=𝑦(𝑘)+𝑒(𝑘)=𝜃^𝑇(𝑘−1)𝜑(𝑘−1)+𝑒(𝑘)=−^𝑛0_𝑖=1𝑎_𝑖(𝑘−1)𝑦(𝑘−𝑖)+^𝑁𝑁_𝑖=1𝑗=1𝑏̂_𝑖,𝑗(𝑘−1)𝛼_𝑗(𝑘−1)𝑢(𝑘−𝑖)+𝑒(𝑘)=−^𝑛0_𝑖=1𝑎_𝑖(𝑘−1)𝑦(𝑘−𝑖)+^𝑁_𝑖=1𝑏^′_𝑖𝑢(𝑘−𝑖)+𝑒(𝑘) (4.1) and the adaptive control law (3.1), replacing 𝑆(𝑧) by 𝑆(𝑘,𝑧), as1𝑢(𝑘)=𝑏^′₁{^𝑁−1_𝑖=1(𝑠₁̂(𝑘−1)𝑎_𝑖(𝑘−1)−𝑠̂_𝑖+1𝑠̂(𝑘−1))𝑦(𝑘−𝑖)+𝜎(𝛽)₁(𝑘−1)𝑎_𝑁(𝑘−1)𝑦(𝑘−𝑁)−^𝑁−1_𝑖=1(𝑠₁̂(𝑘−1)𝑏^′_𝑖+𝑏^′_𝑖+1−𝑠̂)𝑢(𝑘−𝑖)₁(𝑘−1)𝑏^′_𝑁𝑢(𝑘−𝑁)+^𝑁_𝑖=1𝑡_𝑖𝑐(𝑘−𝑖+1)−𝑠₁̂(𝑘−1)𝑒(𝑘)}, (4.2) where 𝑆(𝑧,𝑘−1)=∑^𝑁_𝑖=1𝑠_𝑖̂(𝑘−1)𝑧^𝑁−𝑖, 𝑇(𝑧)=∑^𝑁_𝑖=1𝑡_𝑖𝑧^𝑁−𝑖, 𝑅(𝑧)=𝐵^′(𝑧)=∑^𝑁_𝑖=1𝑏^′_𝑖𝑧^𝑁−𝑖, and the binary-valued function𝜎(𝛽)=1if𝛽=0,0otherwise (4.3) have been used. By combining (4.1) and (4.2), the discrete-time closed-loop system can be written as𝑥(𝑘)=Λ(𝑘−1)𝑥(𝑘−1)+Ψ₁1(𝑘−1)𝑒(𝑘)+𝑏^′₁Ψ₂^𝑁_𝑖=1𝑡_𝑖𝑐(𝑘−𝑖+1), (4.4) where𝑥(𝑘)=[𝑦(𝑘)𝑦(𝑘−1)⋯𝑦(𝑘−𝑛₀+1)𝑢(𝑘)𝑢(𝑘−1)⋯𝑢(𝑘−𝑁+1)]^𝑇,Ψ₁(𝑘−1)=[10⋯0−𝑠₁̂(𝑘−1)/𝑏^′_1𝑛0+10⋯0]^𝑇∈ℝ^{(𝑛0+𝑁)×1},Ψ₂=[00⋯01^}_𝑛0+10⋯0]^𝑇∈ℝ^{(𝑛0+𝑁)×1},=−𝑎Λ(𝑘−1)₁(𝑘−1)−𝑎₂(𝑘−1)⋯−𝑎_𝑛0−1(𝑘−1)−𝑎_𝑛0(𝑘−1)𝑏^′₁𝑏^′₂⋯𝑏^′_𝑁−1𝑏^′_𝑁10⋯0000⋯0001⋯0000…00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00⋯10𝑓00⋯00₁(𝑘−1)𝑓₂(𝑘−1)⋯𝑓_𝑛0−1(𝑘−1)𝑓_𝑛0(𝑘−1)−ℎ₁(𝑘−1)−ℎ₂(𝑘−1)⋯−ℎ_𝑁−1(𝑘−1)−ℎ_𝑁00⋯0010⋯0000⋯0001⋯00⋮⋮⋱⋮⋮⋮⋮⋱⋮⋮00⋯0000⋯10 (𝑘−1)(4.5) with 𝑓_𝑛0(𝑘−1)=(1/𝑏^′₁)[𝑠₁̂(𝑘−1)𝑎_𝑛0(𝑘−1)−(1−𝜎(𝛽))𝑠̂_𝑛0+1(𝑘−1)], 𝑓_𝑖(𝑘−1)=(1/𝑏^′₁)[𝑠₁̂(𝑘−1)𝑎_𝑖(𝑘−1)−𝑠̂_𝑖+1(𝑘−1)], for 𝑖∈{1,2,…,𝑛₀−1}, ℎ_𝑁(𝑘−1)=(𝑏^′_𝑁/𝑏^′₁)𝑠₁̂(𝑘−1), and ℎ_𝑖(𝑘−1)=(1/𝑏^′₁)[𝑠₁̂(𝑘−1)𝑏^′_𝑖+𝑏^′_𝑖+1] for 𝑖∈{1,2,…,𝑁−1}. Note that 𝑎_𝑖(𝑘−1) are uniformly bounded from Lemma 3.2. Also, 𝑠_𝑖̂(𝑘−1) are uniformly bounded from the resolution of an equation similar to (3.3) by replacing the polynomials 𝐴(𝑧) and 𝑆(𝑧) by 𝐴(𝑘−1,𝑧) and 𝑆(𝑘−1,𝑧), respectively. The following theorem, whose proof is made in Appendix B, establishes the main stability result of the adaptive control system.

Theorem 4.1 (main stability result). (i) The adaptive control law stabilizes the plant model (3.14) in the sense that {𝑢(𝑘)} and {𝑦(𝑘)} are bounded for all finite initial states and any bounded reference input sequence {𝑐(𝑘)} subject to Assumption 1.
(ii) The control and measured signals of the continuous-time dynamic system, 𝑢(𝑡) and 𝑦(𝑡), are bounded for all t.

5. Simulations

A continuous-time dynamic system defined by the matrices𝐴=017.50011.50000−1000110−150,𝐵=−10.05,𝐶=0110.05 (5.1) in the state-space and by the transfer function𝐺(𝑠)=𝑠−1(𝑠−5)(𝑠+3.5)(1+0.05𝑠+15)+0.05𝑠+10 (5.2) is considered. This plant is supposed unknown and the adaptive control strategy described in the paper will be used to stabilize it. Such a strategy is based on a discretization process using a FROH with 𝛽=0.3 for a slow sampling time 𝑇=0.3 and a multirate device with 𝑁=3 to place the zeros of the estimated discretized plant model within the stability domain. The time-varying discrete transfer-like function corresponding to such an estimated model is 𝐺_𝑑(𝑧,𝑘)=(𝐵_𝑑(𝑧,𝑘)/𝐴_𝑑(𝑧,𝑘))=((𝑏̂₁(𝑘)𝑧²+𝑏̂₂(𝑘)𝑧+𝑏̂₃(𝑘))/𝑧(𝑧²+𝑎₁(𝑘)𝑧+𝑎₂(𝑘))) with 𝑏̂_𝑖(𝑘)=∑³_𝑗=1𝑏̂_𝑖𝑗(𝑘)𝛼_𝑗(𝑘), for all 𝑖∈{1,2,3} and all 𝑘∈ℤ⁺₀. The estimates 𝑏̂_𝑖𝑗(𝑘) are provided by the estimation algorithm described in Section 3.2.1 with 𝛾(𝑘)=0.01 for all 𝑘∈ℤ⁺₀ and initialized with 𝑃(0)=55×𝐼₁₁, 𝜃_𝑎(0)=[−12.0793.92]^𝑇, 𝜃_𝑏,1(0)=[0.4920.3790.33]^𝑇, 𝜃_𝑏,2(0)=[−0.435−0.556−0.759]^𝑇, and 𝜃_𝑏,3(0)=[0.0190.0660.139]^𝑇. Also, 𝜂_𝑇(𝑘)=𝜐₁‖𝜑(𝑘−1)‖+𝜐₂, with 𝜐₁=0.0055 and 𝜐₂=0.01, is used as upper bound for |𝜂_𝑇(𝑘)| and 𝜇=1.1 for the dead-zone included in such an algorithm. The values 𝛿=𝛿₀=10⁻⁶ are taken in Step 2 of the estimates modification procedure. The gains 𝛼_𝑗(𝑘), for 𝑗∈{1,2,3}, are on-line updated in order to fix 𝐵(𝑧,𝑘) to the polynomial 𝐵^′(𝑧)=𝑧²−0.25. The control objective is to match the reference model given by 𝐺_𝑚(𝑧)=((𝑧²+0.1𝑧+0.083)/(𝑧³+0.3𝑧²−0.09𝑧−0.027)). Figures 1 and 2 display, respectively, the continuous-time dynamic system measured and control signals under a unitary step external input sequence {𝑐(𝑘)}. Note that both signals are bounded for all time; that is, the stabilization of the closed-loop system is reached. Furthermore, the plant output sequence {𝑦(𝑘)} converges asymptotically to {𝑦_𝑚(𝑘)}. Figures 3, 4, 5, and 6 show the evolution of the estimates, included as components of 𝜃(𝑘)=[𝜃^𝑇_𝑎(𝑘)𝜃^𝑇_𝑏,1(𝑘)𝜃^𝑇_𝑏,2(𝑘)𝜃^𝑇_𝑏,3(𝑘)]^𝑇, during the simulation while Figure 7 displays the evolution of the multirate gains 𝑔(𝑘)=[𝛼₁(𝑘)𝛼₂(𝑘)𝛼₃(𝑘)]^𝑇. Note that therefore the estimated parameters as the multirate gains converge to a set of constant values as 𝑡 tends to infinite. Finally, Figures 8(a) and 8(b) show the evolution of the continuous-time dynamic system measured signal, each one in a different time interval, if the same estimation algorithm without the dead-zone is used to stabilize the system. The behavior of the adaptive control system is improved with the inclusion of the relative adaptation dead-zone in this particular example as one can see by comparing the signal 𝑦(𝑡) in Figure 1 with those displayed in Figures 8(a) and 8(b). This conclusion cannot be generalized for all cases since one can search examples where the inclusion of the dead-zone does not improve the performance of the control system. However, the inclusion of the relative dead-zone guarantees the stabilization of the adaptive control system under the presence of unmodeled dynamics, which is what justifies its inclusion in the parameters estimation algorithm.

Figure 1 

Continuous-time dynamic system and reference model measured signals.

Figure 2 

Continuous-time dynamic system control signal.

Figure 3 

Estimates of the parameters 𝑎1 and 𝑎2.

Figure 4 

Estimates of the parameters 𝑏1,1, 𝑏1,2 and 𝑏1,3.

Figure 5 

Estimates of the parameters 𝑏2,1, 𝑏2,2 and 𝑏2,3.

Figure 6 

Estimates of the parameters 𝑏3,1, 𝑏3,2 and 𝑏3,3.

Figure 7 

Multirate gains evolution.

(a)

(b)

Figure 8 

Continuous-time dynamic system measured signal without dead-zone in the parameters estimation algorithm.

6. Conclusions

An adaptive control strategy for stabilizing linear time-invariant continuous-time dynamic systems, being possibly of inverse nonstable, and subject to the presence of unmodeled dynamics has been presented. The control design is based on the discrete-time model reference adaptive (MRAC) control method. Therefore, an inversely stable discretized model of the continuous-time dynamic system is required to achieve the stabilization objective while matching a freely chosen discrete-time reference model. Such a requirement is guaranteed by using a fractional-order hold (FROH) combined with a multirate device with fast sampling input in the discretization process of the continuous-time system. In this context, an estimation algorithm is used to on-line update the multirate gains such that the zeros of the transfer-like function associated to the estimated model of the discretized system are within the open-unit complex circle at all sampling instants. The parameters of such an estimated model are used to on-line parameterize the adaptive control law. The estimation algorithm includes a relative adaptation dead-zone for dealing with the presence of unmodeled dynamics. The stability of the adaptive control system is proved under the assumption that the nominal model of the continuous-time dynamic system is observable, an upper-bound of its order is known, and the contribution of the unmodeled dynamics to the system output is sufficiently small. Finally, the performance of the adaptive control system is shown by means of simulation results. The stabilization of the system is manifested although the intersample behavior of the measured signal could be improved. A future potential research may be the use of the same control technique in a multiestimation scheme for improving such an inter-sample behavior.

Appendices

A. Proof of Lemma 3.2

(i) Equation (32b) and the matrix inversion lemma lead to 𝑃⁻¹(𝑘)=𝑃⁻¹(𝑘−1)+𝛾⁻¹(𝑘)𝑠(𝑘)𝜑_𝑛(𝑘−1)𝜑^𝑇_𝑛(𝑘−1)>0 for all 𝑘∈ℤ⁺ provided that 𝑃(0)=𝑃^𝑇(0)>0. Then, {𝑃(𝑘)} is a nonnegative and monotonic nonincreasing matrix sequence. Thus, 0≤𝑃(𝑘)≤𝑃(0) and 𝑃(𝑘) asymptotically converges to a finite limit as 𝑘→∞.

(ii) By considering the nonnegative sequence 𝑉(𝑘)=𝜃^0𝑇(𝑘)𝑃⁻¹(𝑘)𝜃⁰(𝑘)+Tr{𝑃(𝑘)} and using the matrix inversion lemma in (32b), it follows that𝑉(𝑘)−𝑉(𝑘−1)=−𝑠(𝑘)((𝑒⁰_𝑎𝑛(𝑘))²−(𝜂_𝑎𝑛(𝑘))²)𝛾(𝑘)+𝑠(𝑘)𝜑^𝑇_𝑛(𝑘−1)𝑃(𝑘−1)𝜑_𝑛(𝑘−1)≤−((𝜇²−1)/𝜇²)𝑠(𝑘)(𝑒⁰_𝑎𝑛(𝑘))²𝛾(𝑘)+𝑠(𝑘)𝜑^𝑇_𝑛(𝑘−1)𝑃(𝑘−1)𝜑_𝑛(𝑘−1)≤−((𝜇²−1)/𝜇²)(𝑓(𝑘))²𝛾(𝑘)+𝑠(𝑘)𝜑^𝑇_𝑛(𝑘−1)𝑃(𝑘−1)𝜑_𝑛(𝑘−1)≤0, (A.1) where (32a) and the definition of the “a priori” estimation error have been used. Then, 𝑉(𝑘)≤𝑉(0)<∞ and ‖𝜃⁰(𝑘)‖≤(𝜆_max{𝑃(0)}/𝜆_min{𝑃(0)})‖𝜃⁰(0)‖+𝜆_max{𝑃(0)}Tr{𝑃(0)}<∞ where 𝜆_max(𝑀) and 𝜆_min(𝑀) denote the maximum and the minimum eigenvalues of the matrix 𝑀, respectively. It implies that 𝜃⁰(𝑘) and then also 𝜃⁰(𝑘) are uniformly bounded. Then, 𝜃(𝑘) is also bounded since the modification algorithm guarantees the boundedness of 𝑀(𝑘) provided that 𝜃⁰(𝑘) is bounded. Moreover, 𝑉(𝑘) asymptotically converges to a finite limit as 𝑘→∞ from its definition and the fact that such a sequence is nonnegative and monotonic non-increasing. Then, 𝜃⁰(𝑘), and also 𝜃⁰(𝑘), converges to a finite limit as 𝑘→∞ since 𝑃(𝑘) also converges as it has been proved in (i). Then, 𝑀(𝑘) and 𝜃(𝑘) also converge to finite limits as 𝑘→∞.

(iii) From (A.1), it follows that𝜇²−1𝜇^2𝑘_𝑖=1(𝑓(𝑖))²𝛾(𝑖)+𝑠(𝑖)𝜑^𝑇_𝑛(𝑖−1)𝑃(𝑖−1)𝜑_𝑛(𝑖−1)≤𝑉(0)−𝑉(𝑘)≤𝑉(0)<∞ (A.2) for all 𝑘∈ℤ⁺. Then 𝑓(𝑘)<∞ for all 𝑘∈ℤ⁺ and also lim_𝑘→∞𝑓(𝑘)=0.

(iv) The boundedness and convergence of the estimation model parameters vector together with the nonsingularity of matrix 𝑀(𝑘) (guaranteed by the modification algorithm) imply the boundedness and convergence of the vector 𝑔(𝑘) obtained by resolution of (2.15) replacing 𝑀 by 𝑀(𝑘).

B. Proof of Theorem 4.1

(i) Λ(𝑘−1) is bounded since the estimated plant parameters 𝑎_𝑖(𝑘−1), for 𝑖∈{1,…,𝑛₀}, and the controller parameters 𝑠_𝑗̂(𝑘−1), for 𝑗∈{1,…,𝑁}, are bounded thanks to 𝜃(𝑘−1) and 𝑔(𝑘−1) are bounded for all 𝑘∈ℤ⁺, see Lemma 3.2. The eigenvalues of Λ(𝑘−1) are in |𝑧|<1 since they are the roots of 𝐴_𝑚(𝑧), 𝐴_𝑠(𝑧) and 𝐵^′(𝑧), which are stable. Furthermore,^𝑘_{𝑘′=𝑘0+1}‖Λ(𝑘^′)−Λ(𝑘^′−1)‖²≤𝛾₀+𝛾₁(𝑘−𝑘₀) (B.1) for all integers 𝑘>𝑘₀≥0 and some positive real constants 𝛾₀ and 𝛾₁, with 𝛾₁ being sufficiently small by using slow enough estimation rates via a suitable 𝑃(0) in (32b). Thus, the unforced time-varying system 𝑥(𝑘)=Λ(𝑘−1)𝑥(𝑘−1) is exponentially stable and its transition matrix 𝜙(𝑘,𝑘^′)=∏^𝑘−1_𝑗=𝑘′Λ(𝑗) satisfies ‖𝜙(𝑘,𝑘^′)‖≤𝜌₁𝜎^{𝑘−𝑘′}₀ for all integer 𝑘≥𝑘^′ where 𝜎₀∈(0,1) is an upper bound for the absolute value of the closed-loop stability abscissa and 𝜌₁ is a non-dependent constant [2]. Note that 𝜎₀ depends on the freely chosen reference model. From (4.4), 𝑥(𝑘)=𝜙(𝑘,𝑘₀)𝑥(𝑘₀)+^𝑘_{𝑘′=𝑘0}𝜙(𝑘,𝑘^′)(Ψ₁(𝑘^′−1)𝑒(𝑘^′1)+𝑏^′₁Ψ₂^𝑁_𝑖=1𝑡_𝑖𝑐(𝑘^′−𝑖+1)). (B.2) Then,‖𝑥(𝑘)‖≤𝜌₁𝜎^𝑘−𝑘0₀‖𝑥(𝑘₀)‖+^𝑘_{𝑘′=𝑘0}𝜌₁𝜎^{𝑘−𝑘′}₀(𝜌₂+𝜌₃|𝑒(𝑘^′)|), (B.3) for some positive real constants 𝜌_𝑖, 𝑖∈{1,2,3}, since Ψ₁(𝑘) is uniformly bounded and provided that the sequence {𝑐(𝑘)} is bounded. From (3.16) and (4.1), it follows that 𝑒(𝑘)=𝑒⁰(𝑘)+(𝜃⁰(𝑘−1)−𝜃(𝑘−1))^𝑇𝜑(𝑘−1), and then|𝑒(𝑘)|≤|𝑒⁰(𝑘)|+ℓ(𝑘)𝛿‖𝜑(𝑘−1)‖ (B.4) where ℓ(𝑘) denotes the finite number of times that the instruction while in the Step 2 of the estimates modification algorithm is executed at the current sampling time 𝑘𝑇. By substituting (B.4) in (B.3), one obtains‖𝑥(𝑘)‖≤𝜌₄+^𝑘_{𝑘′=𝑘0}𝜌₅𝜎^{𝑘−𝑘′}₀(|𝑒⁰(𝑘^′)|+ℓ(𝑘^′)𝛿‖𝜑(𝑘^′−1)‖), (B.5) for some positive real constants 𝜌₄ and 𝜌₅ from the boundedness of 𝑥(𝑘₀). By using that 𝜑_𝑛(𝑘−1)=𝜑(𝑘−1)/(1+‖𝜑(𝑘−1)‖) and 𝑒⁰_𝑛(𝑘)=𝑒⁰(𝑘)/(1+‖𝜑(𝑘−1)‖), it follows that‖𝑥(𝑘)‖≤𝜌₄+^𝑘_{𝑘′=𝑘0}𝜌₅𝜎^{𝑘−𝑘′}₀(|𝑒⁰_𝑛(𝑘^′)|+ℓ(𝑘^′)𝛿‖𝜑_𝑛(𝑘^′−1)‖)(1+‖𝜑(𝑘^′−1)‖)≤𝜌₄+^𝑘_{𝑘′=𝑘0}𝜌₆𝜎^{𝑘−𝑘′}₀𝑒⁰_𝑎𝑛(𝑘^′)(1+‖𝜑(𝑘^′−1)‖) (B.6) for some positive real constant 𝜌₆ and, also, by taking into account that (|𝑒⁰_𝑛(𝑘)|+ℓ(𝑘)𝛿‖𝜑_𝑛(𝑘−1)‖)≤𝜌^′((𝑒⁰_𝑛(𝑘))²+𝜑^𝑇_𝑛(𝑘−1)𝑃²(𝑘−1)𝜑_𝑛(𝑘−1))^1/2=𝜌^′𝑒⁰_𝑎𝑛(𝑘) for some positive real constant 𝜌^′. Moreover, from (B.6),‖𝑥(𝑘)‖≤𝜌₄+^𝑘_{𝑘′=𝑘0}𝜌₆𝜎^{𝑘−𝑘′}₀𝜑𝜇(1+^𝑇𝑛(𝑘^′−1)𝑃(𝑘^′−1)𝜑_𝑛(𝑘^′−1)𝛾(𝑘^′))^1/2𝜂_𝑇(𝑘^′)+^𝑘_{𝑘′=𝑘0}𝜌₆𝜎^{𝑘−𝑘′}₀𝑓(𝑘^′)(1+‖𝜑(𝑘^′−1)‖) (B.7) where 𝑒⁰_𝑎𝑛(𝑘^′) has been split into the two additive terms 𝑓(𝑘^′) and 𝑒⁰_𝑎𝑛(𝑘^′)−𝑓(𝑘^′)≤𝜇𝜂_𝑎𝑛(𝑘^′). Then,‖𝑥(𝑘)‖≤𝜌₇+𝜌₆𝜇𝜐₃1−𝜎₀Sup_{𝑘0≤𝑘′≤𝑘}𝜑{(1+^𝑇𝑛(𝑘^′−1)𝑃(𝑘^′−1)𝜑_𝑛(𝑘^′−1)𝛾(𝑘^′))^1/2}Sup_{𝑘0≤𝑘′≤𝑘}{‖𝑥(𝑘^′+)‖}_{𝑘𝑘′=𝑘0}𝜌₆𝜎^{𝑘−𝑘′}₀𝑓(𝑘^′)(1+‖𝜑(𝑘^′−1)‖) (B.8) for some positive constant 𝜌₇ by taking into account that 𝜂_𝑇(𝑘)=𝜐₃‖𝑥(𝑘−1)‖+𝜐₄ (see Remark 3.1(iii)). Furthermore, the right-hand side of (B.8) is monotonic nondecreasing in 𝑘. Then, for 𝑘≫𝑘₀, it follows that‖𝑥(𝑘)‖≤𝜌₈+^𝑘_{𝑘′=𝑘0}𝜌₉𝜎^{𝑘−𝑘′}₀𝑓(𝑘^′)(1+‖𝜑(𝑘^′−1)‖) (B.9) provided that 𝜐₃<((1−𝜎₀)/𝜌₆𝜇)(Sup_{𝑘0≤𝑘′≤𝑘}{(1+𝜑^𝑇_𝑛(𝑘^′−1)𝑃(𝑘^′−1)𝜑_𝑛(𝑘^′−1)/𝛾(𝑘^′))^1/2})⁻¹, for some positive constants 𝜌₈ and 𝜌₉. By taking into account that ‖𝑔(𝑘)‖ is bounded for all 𝑘∈ℤ⁺₀, then‖𝜑(𝑘−1)‖≤𝜌₁₀+𝜌₁₁Sup_{𝑘0≤𝑘′≤𝑘}{‖𝑥(𝑘^′−1)‖} (B.10) for some positive constants 𝜌₁₀ and 𝜌₁₁. By substituting (B.10) in (B.9) and by taking into account that the right-hand side of (B.9) is monotonic nondecreasing in 𝑘, it follows that,Sup_{𝑘0≤𝑘′≤𝑘}{‖𝑥(𝑘^′)‖}≤𝜌₁₂+^𝑘_{𝑘′=𝑘0}𝜌₁₃𝑓(𝑘^′)Sup_{𝑘0≤𝑘′′≤𝑘′}{‖𝑥(𝑘^′′−1)‖} (B.11) for some positive constants 𝜌₁₂ and 𝜌₁₃, where the boundedness of 𝑓(𝑘) has been used, see Lemma 3.2. The use of Gronwall's Lemma [24] in (B.11) leads toSup_{𝑘0≤𝑘′≤𝑘}{‖𝑥(𝑘^′)‖²}≤𝜌₁₄+^𝑘_𝑖=𝑘0(_{𝑖<𝑗<𝑘}(1+𝜌₁₅𝑓²(𝑗)))𝜌₁₄𝜌₁₅𝑓²(𝑖)<∞ (B.12) for some positive constants 𝜌₁₄ and 𝜌₁₅. It implies that the sequence ‖𝑥(𝑘)‖ is bounded for all 𝑘∈ℤ⁺ provided that the initial condition is bounded. Then, the sequences {𝑢(𝑘)} and {𝑦(𝑘)} are also bounded.

(ii) The adaptive control algorithm ensures that there are not finite escape times, and then the boundedness of the sequences {𝑢(𝑘)} and {𝑦(𝑘)} guarantees that of the continuous-time signals from continuity arguments of the solutions of differential equations.

Acknowledgments

The authors are very grateful to MCYT for its support through Grants nos. DPI2006-00714 and DPI 2009-07197. They are also very grateful to the reviewers for their useful comments, which have helped to improve the previous versions of the manuscript.