A novel wavelet-based optimal linear quadratic tracker for time-varying systems with multiple delays

A new method for solving optimal tracking control of linear quadratic time-varying systems with multiple time delays in state and input variables and with combined constraints is presented in this paper. By using the relations of Chebyshev wavelets, we simulate the optimal tracking problem to a static optimization one. This alternative method is applied on different optimal tracking systems and simulation results demonstrate the effectiveness of the proposed method.


Introduction
The mathematical simulation of a time-delay system leads to a system of differential equations of delayed (retarded) type [1]. One class of these equations, the integro-differential equations, was first studied by Volterra who investigated time delay phenomena in various systems [2]. In optimal control problems, when we try to keep the output or state near a desired output or state, we are dealing with tracking problems. We find that in both state and output time-delay or regulator systems, the desired (reference) state or output is zero and in time-delay or regulator tracking system the error is to be made zero [3]. Time-delay is a common phenomenon in engineering problems [4]. Optimal time-delay tracking control as a combination of time-delay optimal control and tracking control, aims at finding the optimal control law to minimize the given performance index to make the system output track the reference in an optimal way. The tracking system is widely used in aerospace and mechanical systems, robot control. For example, consider an antenna control system to track an aircraft.
Most previous studies which have been done to solve the optimal time-delay tracking control problem, used discrete-time strategies. Ref. [5] presented an iterative method by using discretization to find the suboptimal control of a linear quadratic time-varying system with multiple delays. [6] presented a discretization approach by using the Newton center interpolation formula and the linear interpolation techniques for systems with multiple discrete and distributed time delays. In [7] an optimal tracking controller for discrete time-delay systems based on a sensitivity approximation approach was designed in which the problem under consideration was transformed into a series of difference equations without time-advance on delayed terms. In [8], the authors used a successive approximation method for delays systems with constant matrices and the resulting problems form two-point boundary value problems (TPBVPs). [9] proposed a suboptimal tracking method which obtained by finite iterations of a solution 419 of TPBVPs; for continuous-time control systems it provides good methodology but obtaining the solutions of these TPBVPs is difficult. [10] by introducing a sensitively parameter, the discrete time systems transformed into series of TPBVPs. In [11] by applying an approximation approach of differential equations, TPBVP derived from the optimal tracking problem was transformed into a sequence of linear TPBVPs without delays. Ref. [12] used a discrete-time strategy to design an optimal controller for multiple-input and multiple-output continuous time systems with multiple delays in states, inputs and outputs which the Chebyshev quadrature formula together with a linear interpolation method was employed to get an extended discrete-time model from the continuoustime multiple time-delays system. [13] converted a continuous-time input-state delayed system into an equivalent discrete-time input-state delayed model and its extended discrete-time delay-free model. Refs. [14,15] proposed a discrete variable transformation by which discrete-time linear-quadratic systems with multiple input delays were transformed into equivalent delay-free ones. [16] proposed an adaptive dynamic programming algorithm to deal with the optimal tracking problem for delay discrete-time systems. Also, there are some works which studied freedelay linear tracking optimal control, for example, [17] continues-time systems, [18] discrete-time systems and [19] time-varying systems. A discrete-time methodology has notable disadvantages. For example, the discretized system is described by using the extended high-order state-space equation. Computational difficulties with the discretetime strategy occur in the algorithm, and in the case we have to impose some constraints to the system. In this paper we introduce a novel method to solve constrained linear quadratic time-delay tracking optimal control systems. By using the parameterization method we convert the original problem to a static optimization one. This state and control parameterization method is based on Legendre and Chebyshev wavelets which consist of Legendre and Chebyshev scaling functions. The main motivations of this research are summarized as: • A continuous-time model of the optimal tracking control of a linear quadratic time-varying system with multiple delays is obtained. A major advantage of continuous-time models is that they avoid dependence on a particular timescale. Moreover, many standard numerical procedures are available to solve the simulated model.
• In the optimal tracking control, we would like on one hand, to keep the error small, but on the other hand, we must not pay higher cost to large inputs; hence, we have to try various values of the weighting matrices. From this fact, we conclude that we need a method provides good tracking in which with no concern about the algorithm of the solution we can change these matrices.
• Physical considerations imply that some constraints should be imposed on the optimal tracking control systems and unconstrained systems are less involved. An efficient method for solving a time-delay optimal tracking control problem should easily be able to resolve the problem in the cases there are constraints on the controls and states; in this case, the method can incorporate them directly into the model of the problem.
• The proposed method has a good future and a high degree of flexibility. For example, it is possible that there is a reverse time term like x(t f − t) in the state equation and/or a time-delay term like u(t − h) in a performance index, or, since fractional order dynamics appear in engineering problems [26], the state quation is fractional order; the method is capable of executing such cases. Also, the method must guarantee intersample constraint satisfaction and can be used to solve the optimal tracking problem in situations where there are multiple delays or no delays (regulator systems), the plant matrices are time-varying and/or constants.
In Section 2, the concepts of Legendre and Chebyshev wavelets required to model the problems are given. We will apply the state-control parameterization method to the tracking delay systems in Section 3. Corresponding problems will be discussed in Section 4.

Notation
The transpose of a matrix O is written O . 0 and I denote the zero and identity matrices. Kronecker product of O and I q is denoted byÔ, Kronecker product of O and I r is denoted byǑ, that is,Ô = O ⊗ I q ,Ǒ = O ⊗ I r . The superscript " * " indicates the optimal condition. C[0, t f ] denotes real-valued continuous functions over [0, t f ]. w(t) denotes the desired wavelets vector consists of Legendre or Chebyshev wavelets {w nm (t)}. The subscript "w" is for the desired wavelets. The semicolon in defining matrices is used to separate the rows.

420
WAVELET-BASED OPTIMAL LINEAR QUADRATIC TRACKER

Legendre and Chebyshev wavelets
In applications of classical Legendre and Chebyshev wavelets, we have to use large k in most cases, for example, see [24]. To overcome this drawback and others, we presented new definitions of Legendre and Chebyshev wavelets in [20], [21]. For purposes of comparison with classic types, we have used ξ k−1 in their definitions. Here we present more general definitions for Legendre and Chebyshev wavelets. We mentioned previously in [28], we see in the applications of these wavelets that we set k = 2 in their definitions. So by taking N = ξ k−1 in the definition given in [20], Legendre wavelets φ nm are defined on [0,1] as where t ∈ [0, 1] is as an independent variable, P m are the well-known Legendre polynomials, N ∈ N ≥2 is an arbitrarily selected scaling parameterand specifies the number of subintervals, n = 1, 2, . . . , N refers to the number of subinterval and specifies the location of the subinterval, m = 0, 1, . . . , M − 1 is the degree of P m and c m is c m = √ 2m + 1. Legendre wavelets are an orthonormal set with respect to the weight functions ω n (t) = 1. Again by taking N = ξ k−1 in the definition given in [21], Chebyshev wavelets ψ nm are defined on [0, 1] as where T m are the well-known Chebyshev polynomials, N, n and m are the same as before and c 0 = 1/ √ π, c m =0 = √ 2/ √ π. They form an orthogonal basis with respect to the weight function ω n (t) = 1/ 1 − (2N t − 2n + 1) 2 . We can expand a function f (t) in a series of Legendre or Chebyshev wavelets denoted by {w nm (t)} (we truncate the series with the M th term in N subintervals) as where f w is a 1 × N M vector that consists of constants, w(t) as a N M × 1 vector of any of these two wavelets and The coefficients of these scaling functions can be approximated as follows We use the definitions of Legendre and Chebyshev wavelets in (1), (2) because they have many advantages over classic definitions of them, [22][23][24].

Theorem 1
A twice differentiable function f (t), defined on [0, 1], with bounded second derivatives, say |f (t)| ≤ ρ, can be expanded as an infinite sum m of Legendre or Chebyshev wavelets, and this series converges uniformly to f (t).

Proof
The proof is the same as that was presented in [21].

The operational matrices of Legendre and Chebyshev wavelets
The following properties of the desired wavelets will be used for applying the wavelet-based method: where P w is the integration operational matrix, the Riemann-Liouville fractional integral operator and the Caputo fractional derivative of order α denoted by I α and D α , P α w is the fractional integration operational matrix of the desired wavelets in the Riemann-Liouville sense, D w is the derivative operational matrix, D α w is the Caputo fractional derivative operational matrix,f w is the product operational matrix of the desired wavelets for f w , D ι is the delay operational matrix for a time-delay h ι , Γ w is the integration matrix of the product of two desired wavelets vectors on [0, 1], and D t is the piecewise delay operational matrix for a piecewise delay h(t). D ι and D t are the same for the both wavelets, so we do not use the subscript "w" for them. For more details, see [20], [21], [27], [28]. In [21], the error generated by using the wavelets approximation has been studied. All the properties mentioned above have been presented in the general form in these texts and we just set N = ξ k−1 in their formulas.

Optimal tracker for linear quadratic time-delay system
Consider a linear time-varying system with multiple time delays described bẏ x(0) = x 0 (14) and a quadratic performance index as where x(t) and u(t) are qand r-dimensional state and control vectors, respectively, A(t), A h (t), B(t) and B h (t) are piecewise-continuous matrices of compatible dimensions, ∆A(t) and ∆A h (t) denote the parameter uncertainties, h x and h u denote time delays, d(t) as a q-dimensional vector represents disturbances, f (t) is a qdimensional initial state vector function, g(t) is an r-dimensional initial control vector function, x 0 is an initial condition vector, Q is a positive semi-definite matrix, R is a positive definite matrix and r(t) is a q-dimensional desired or reference state vector. The main purpose of the matrix T is to ensure that the error at the terminal time is as small as possible. So, this matrix should be positive semi-definite. Our objective is to control this system in such a way that the state x(t) tracks the desired state r(t) as close as possible during the time interval [0, t f ]. The optimal tracking problem is to find u * (t), x * (t) and J * for the time-delay system (12)- (14) such that the performance index in (15) is minimized.

422
WAVELET-BASED OPTIMAL LINEAR QUADRATIC TRACKER

Open-Loop Optimal Controller (OLOC)
We find OLOC from solving the optimization problem which is presented in the following.

Theorem 2
The optimal tracking problem as minimizing the performance index (15) for the plant described by (12)- (14), by using the properties of the desired wavelets can be solved by forming a quadratic programming (QP) problem as where ℵ w , Λ w , and b w are known constant matrices and χ w as a combination of parameters of the state and control parameterization is the solution of the problem in terms of the desired wavelets. The function J(χ w ) in this QP is convex.

Proof
We can use two procedures based on (7a) or (7c). First we use the integration procedure. For using the wavelet method, we must change the range of t such that 0 ≤ t ≤ 1. We set τ = t/t f to rescale (12), this implieṡ where Assume that U is a class of such admissible controls. To each u ∈ U, the system has a solution as the corresponding state function x(·|u). Let us define a new state vector as Using Theorem 1, extending (3)-(5) for vector functions, we parameterize this new state and the control vectors as followsx wherex w and u w are N qM × 1 and N rM × 1 column vectors of unknown parameters and Using the latter point, we expand the initial condition and the desired state as where x 0 w and r w are known N qM × 1 column vectors given by In (23), a 0 n0 := [x 0 , 0 1×q(M −1) ] and in (24), the constant coefficients {r a nm }, for a = 1, 2, . . . , q, are obtained by (6).
Using this and (9) in (22), we can write , so according to (13) we have (3), where by taking n x = N τ x and n u = N τ u , f xw and g uw are, respectively, N qM × 1 and N rM × 1 column vectors of constants defined by f a nm and for b = 1, 2, . . . , r, g b nm can be calculated using formula (6). Thus from (18), (26), (9), we find Now we express the time-varying matrices in (16) in terms of the desired scaling functions. Using (3) for A(τ ), we can write (32) Similar to r(τ ), the disturbance can be expressed as First, we substitute (17) in (16). Then, as was mentioned in [24], we integrate the resulting equation from 0 to τ , substitute (18), (21), (22), (25), (26), (29)-(33) and then we use (7a) and (8), hencê Thus, by factoringx w and u w representing the new state and control vector of the problem, we have where we let s = N M .
From the definition of this wavelet, we see the fact that the time interval [0, 1] is divided into N subintervals. In order to ensure continuity in the obtained states across these subintervals, the following compatibility constraint [25] is enfotced at the interface points (τ ι ) of two consecutive subintervals: We assume that r(t) is in C[0, t f ], so for all a, it is necessary that , the compatibility constraint for the new state is expressed aŝ where by ρ ι := [ where By setting (17) in the (rescaled) performance index (15), using (18) and (10), we can write As a result Taking (37), (36) and (34) together, the optimal tracking control problem is transformed into a QP problem as: where we set where the vector χ and the matrix ℵ are the same as (38) and (39).
For convexity, we must show that for all λ ∈ [0, 1], z 1 , z 2 , we have Since Γ w and w(1)w (1) are symmetric matrices, it is clear that ℵ w is symmetric; this yields z 1 ℵ w z 2 = z 2 ℵ w z 1 , for z 1 , z 2 ∈ R N M (q+r) . Thus by assuming J(z . ) ≥ 0, we can write so we see that J is a convex function, if and only if ℵ w is (symmetric) positive semidefinite. For showing ℵ w is positive semidefinite, we know a block diagonal matrix is positive semidefinite if and only if its diagonal blocks are positive semidefinite. Assume that ε 1 , ε 2 , . . . , ε N M are the eigenvalues of Γ w and υ 1 , υ 2 , . . . , υ q are the eigenvalues of Q, then the eigenvalues of Γ w ⊗ Q are ε i υ j , 1 ≤ i ≤ N M , 1 ≤ j ≤ q. Also assume that {ε i } are corresponding eigenvectors of Γ w and {υ j } are corresponding eigenvectors of Q. Then η := {ε 1 ⊗ υ 1 , ε 1 ⊗ υ 2 , . . . , ε 2 ⊗ υ 2 , . . . , ε N M ⊗ υ q } are corresponding corresponding eigenvectors of Γ w ⊗ Q. It follows directly that since Γ w is positive definite and Q is positive semidefinite, Γ w ⊗ Q is positive semidefinite. We can obtain similar statements for other blocks of ℵ w ; hence ℵ w is positive semidefinite.

Corollary 1
Let u * exact (t) be the OLOP of the original problem and u * (t) be the OLOP of the QP problem. Then for n x , n u ∈ N, lim M →∞ u * (t) = u * exact (t). Also by assuming similar statement for J, lim The first goal is to find χ from solving the optimization problem which is static in nature. Standard numerical methods are available to solve this QP problem. Hence we do not need a special program. We can use the quadprog function provided by the optimization toolbox in MATLAB. This toolbox widely has been used to solve constrained and unconstrained optimization problems and like the case we have simple delay or regulator systems, [21], [29], in which the value of optimal cost is one of its default outputs. In this work we use the interior-point-convex algorithm of the quadprog function in MATLAB R2013b; the proposed algorithm provides accurate solutions, and is fast and stable. Detailed information about handling various cases of plant matrices, constraints and ..., can be found in our previous works.
By doing so we obtain the OLOC and its corresponding state from (18)- (20). In the next section, we define a closed-loop suboptimal controller which obtained from the results of this section. But before doing this, we present some important remarks.

Remark 1 (Systems with multiple state and input delays)
For a system with multiple delays described over 0 ≤ t ≤ t f aṡ we replace (Ã hw + ∆Ã hw )D x and (Ã hw + ∆Ã hw )f xw in (40)-(43) by V µ=1 (Ã hµw + ∆Ã hµw )D µ and V µ=1 (Ã hµw + ∆Ã hµw )f µw , respectively, in which every D µ is obtained from (9) and every f µw is constructed according to the value of h µ from (27), where n µ = N h µ /t f . Also for input delay terms, we do the same and replacẽ B hwĎ u andB hw g uw by W ν=1B hν wĎ ν and W ν=1B hν w g νw , respectively, in which every g νw is constructed according to h ν from (28) and n ν = N h ν /t f .

Remark 2 (Systems with piecewise constant delays or time-varying delays)
For a system with a piecewise constant state delay h(t) described over 0 ≤ t ≤ t f bẏ we replace D x and f xw by D t and f t w which D t is the piecewise delay operational matrix of Legendre or Chebyshev wavelets given in (11) and f t w is the piecewise wavelets expansion of f (t − h(t)). For a piecewise constant input delay, we do the same. Moreover, for a system with a time-varying time-delay h(t) by using the time-partition technique, as we did in Ref. [21], after selecting N , in each subinterval [t i−1 , t i ], where i = 1, 2, . . . , N , we take h i ≈ h(ti−1)+h(ti) 2 and construct the matrix D t . Hence we can apply this operational matrix on such systems in which delays are time-varying.

Remark 3 (Fractional tracking delay systems)
If the state equation in (12) is given with the Caputo derivative of order α as where 0 < α < 1, by applying the α-integral (the Riemann-Liouville fractional integral) to both sides of (44) and using (7b), we just replace t f andP w in (40), (41) of the first QP model by t α f andP α w . Using (7d), we do the same for the second QP model and replaceD w in (42), (43) byD α w . So, we extend the QP models to such systems.

Closed-Loop Optimal Controller (CLOC)
One of main motivations of defining new wavelets is (by eliminating inaccuracies of the previous definitions) to use them to obtain the closed-loop solutions [21]. Engineers prefer the closed-loop optimal control to an optimal control problem. Open-loop solutions are not of great interest in practice. Because in practical engineering problems, we need a closed-loop solution in the major of cases and open-loop solutions have not the practical significance.
We propose a suboptimal controller which has been applied on some real world optimization problem previously as Example 4 in [30] and Example 6 in [28]. The proposed CLOC of the system described in (12) where we use the results of OLOC to find the optimal feedback matrices K and K h and the vector l(t). Using (17) in (45), we can write Using the relations of these wavelets, we can compute the unknowns in (46). The proposed method is summarized in Algorithm 1. First, we use our results to the interval [h x , t f ] to compute K and K h . Then by using K and from the properties of the desired wavelets, we compute l(t).
Depending on the other types of the system (systems with multiple delays or a piecewise constant delay or a time-varying delay or with a piecewise reference state), we can propose another types of CLOCs.

Algorithm 1
Input: Λ w , ℵ w , b w Output: χ w for x * (t) and u * (t), and consequently u * CLOC (t) Initialization: choose N from the problem and set M = 8; select the method and form the QP model 1. Use the quadprog solver to find x * (t) and u * (t) 2. Design the CLOC for the tracking problem 3. Find K, K h and l(t) 4. Plot u * CLOC (t) and compare with u * (t) 5. If the results are in good agreement, go to 6 otherwise go back to 2 6. END

Numerical examples
We shall examine the proposed method in this section for solving optimal tracking control of some systems. We first show applicability of the method on a simple problem and then on some complex problems.

Example 1
The problem involves the minimization of subjected to the system of delayed differential equations and initial conditions such aṡ for the case which x(t f = 15) is free and admissible optimal control and states are unbounded, see [5].
After rescaling (47)-(50), we choose N = 30; then we use the method with the both wavelets. The solution curves are presented in Fig. 1. A comparison of the optimal value of the cost functional (47) as J * is given in Table  1. It is clear that the obtained result is in good agreement. We used two QP models for both wavelets. Although we see that there are some significant differences between elements of these QP models, but the results are very close.
For performing a better tracking by the system we increase the value of weighting matrix Q such that: 0 0 ]. Then x * (t) and u * (t) for this weighting matrix are plotted in Fig. 2; also we get J * = 60.882278290621493 by Legendre wavelet method and J * = 60.882278290621464 by Chebyshev wavelet method. These results mean that when we increase the values of the weighting matrix Q, the state of the new system is able to track the reference state better with lower error, but we have to pay higher cost for larger control effort. Suppose, to achieve a better tracking, that instead of increasing the value of the error weighted matrix Q, we decrease the value of the control weighted matrix R, such as: R new = 0.1R = 0.005. This gives J * = 6.088227829062145, thus we get a lower cost while the graphs of x * (t) and u * (t) are the same as those obtained with new Q (in Fig. 2). We want on one hand, to keep the new state small and on the other hand, we must not pay higher cost to large controls; this leads us immediately to the conclusion that we have to try various values of the weighting matrix R. The computed optimal feedback gain for the two error weighted matrices are given in Table 2.

Example 2
Consider a non-square multi-input multi-output controllable and observable system (see [13]) with the performance index Also, x(t) ∈ R 5 is the state vector, y(t) ∈ R 2 is the output vector, and u(t) ∈ R 2 is the control vector. The problem is to find the optimal states and controls for the time-delay system (51), which minimizes (52). We set: Q = 10 3 I 2 and h x = 0.2, h u = 0.1, and t f = 6. Since the performance index is not in the form that we presented, we define x 6 (t) = 2x 1 (t) + x 3 (t) and x 7 (t) = 1.5x 2 (t) + 1.2x 4 (t) + x 5 (t). Then we reformulate the plant and the performance index matrices, for example, we have Q = 2[0 5×7 ; 0 2×5 , 10 3 I 2 ]. Solving the new system, our results are shown in Figs. 3 and 4.
As one can see in the given tracking systems, we sometimes have to change the control weighted matrix of the system in order to enable it to better track the reference state. By changing the weighting matrix R, in using the method of [13], we have to repeat the steps of the modeling process but in our method we can change it without any concern and even we can replace it by a time-varying matrix (as we shall see later).

Example 3
To see the applicability of our method on time-varying systems, consider a system described by (12) or (44), where in which U is the unit step function or the Heaviside function over 0 ≤ t ≤ t f . This system is to be controlled to minimize the performance index in order that the state x 1 (t) tracks the desired trajectory r(t), where The terminal time is t f = 4 and delays are h x = 0.5, h u = 1.5. In the following, we consider this optimal control problem with different orders and constraints as Cases 1 and 2; we take: Case 1. α = 0.9 and the controls and states of the fractional tracking delay system are unconstrained.
Case 2. α = 1 and a. the controls and states are unconstrained. b.
By the proposed method for Chebyshev wavelet, we find the optimal solutions of the problem for all cases. For modeling the problem with the isoperimetric constraint as Case 2(b), we use the method given in [30], but one should be a little bit careful here. From (17), the term t f 0 r(t)dt must be considered in the reformulated constraint. For Case 3(c), we use the method given in [29]. The solution curves for all cases are given in Figs. 5-8. The optimal values of the performance indices are reported in Table 3. As claimed, we see in Case 1 the fact that the method has the ability to solve fractional time-delay linear-quadratic optimal control problems. For Case 1, we find From the experience of studying optimal control of some time-delay systems in [24], [21], we find that behavioral changes may arise in time-delay systems on each side of delays or times which are associated with delays (for example at h x , h u , t f − h u ). Hence we have to redesign the CLOCs. In Case 2(b) which the isoperimetric constraint is imposed, we redesign the controller The initial designs are shown in Fig. 7(b). Using the given concepts, the redesigned controllers are shown in Fig. 7(c). From Corollary 1, with reference to the OLOCs, it is obvious that the redesigned controllers are in good agreement. For Case 2(b), we find         From the nature of the constraint given in Case 2(c), we have to calculate the feedback gains for three subintervals separately. Comparing Figs. 8(b) and 8(c), we see that the redesigned CLOCs are in good agreement with OLOCs.

432
WAVELET-BASED OPTIMAL LINEAR QUADRATIC TRACKER Moreover,, by modifying the performance index (53) and the desired trajectory of the system given in Case 2(a), we consider another case as: Case 3. For the system in Case 2(a), the new performance index is and a. r(t) = cos t −1 + cos t . b. r(t) = cos t sin t . By making some changes in the model of Case 2(a) without doing any further work, we solve the new problems and this is a distinct advantage. As before, the results are presented in Figs. 9 and 10, and Table 3.
In order to get a better tracking in Cases 3(a) and . By using the proposed method, the optimal states and controls for these cases are shown in Figs. 11 and 12; as can be seen, we reach a better tracking. The optimal values of (54) become J * = 3.876154 and J * = 2.460061, respectively. So, the system with the new time-varying control weighted matrix, tracks the reference state better with lower cost. The errors for the two weighting matrices are shown in Fig. 13. As we did previously, we redesigned the CLOCs in Case 3(a).

Conclusion
An alternative method is introduced to find the optimal control, state and performance index of different types of linear time-varying tracking systems with delays. In the proposed procedure, we can easily change the weighting matrices and impose the combined constraints. When we increase the value of the error weighted matrix, then the output is able to track the reference input better with lower output error, but we have to pay higher cost for larger control effort of the designed system. As can be seen, to better tracking we must try various values of the control weighted matrix. The new optimal tracker presented by this paper can be applied to the tracking system regardless of the system stability, minimum phase properties, the dimension and order of the system, equal number of input and output, the number of delays, the types of desired states and initial functions, and fractional order dynamics.