A novel feedforward–feedback suboptimal control of linear time-delay systems

Accepted Manuscript A novel feedforward-feedback suboptimal control of linear time-delay systems A. Jajarmi, M. Dehghan Nayyeri, H. Saberi Nik PII: DOI: Reference:

S0885-064X(16)00013-3 http://dx.doi.org/10.1016/j.jco.2016.02.001 YJCOM 1284

To appear in:

Journal of Complexity

Received date: 10 June 2015 Accepted date: 26 January 2016 Please cite this article as: A. Jajarmi, M.D. Nayyeri, H.S. Nik, A novel feedforward-feedback suboptimal control of linear time-delay systems, Journal of Complexity (2016), http://dx.doi.org/10.1016/j.jco.2016.02.001 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A novel feedforward-feedback suboptimal control of linear time-delay systems A. Jajarmi1,∗ , M. Dehghan Nayyeri2 , and H. Saberi Nik3 1 Department

of Electrical Engineering, University of Bojnord, Bojnord, Iran of Mathematics, University of Bojnord, Bojnord, Iran 3 Department of Mathematics, Neyshabur Branch, Islamic Azad University, Neyshabur, Iran ∗ Corresponding author: [email protected] 2 Department

Abstract This paper presents a new approach for solving the optimal control problem of linear time-delay systems with a quadratic cost functional. In this approach, a method of successive substitution is employed to convert the original time-delay optimal control problem into a sequence of linear timeinvariant ordinary differential equations (ODEs) without delay and advance terms. The obtained optimal control consists of a linear state feedback term and a forward term. The feedback term is determined by solving a matrix Riccati differential equation. The forward term is an infinite sum of adjoint vectors, which can be obtained by solving recursively the above-mentioned sequence of linear non-delay ODEs. A fast-converging iterative algorithm for this purpose is presented which provides a promising possible reduction of computational efforts. Numerical examples demonstrating the efficiency, simplicity and high accuracy of the suggested technique have been included. Simulation results reveal that just a few iterations of the proposed algorithm are required to find an accurate enough feedforward-feedback suboptimal control.

Keywords: Time-delay optimal control problems; Pontryagin’s maximum principle; successive substitution method; feedforward-feedback suboptimal control.

1

Introduction

Time-delay systems are a very important class of systems whose control and optimization have been of interest to many researchers [1]. The presence of delay makes the analysis and control of such systems much more complicated. In fact, the application of Pontryagin’s maximum principle to the optimal control problems (OCPs) with time-delay results in a system of coupled two-point boundary value problem (TPBVP) involving both delay and advance terms [2]. This problem cannot be solved analytically except in very special cases. Therefore, the main objective of all methodologies for solving the time-delay optimal control problems (TDOCPs) has been the avoidance of solving the abovementioned TPBVP. The solution of TDOCPs has been investigated in the literature from the viewpoint of approximating the infinite dimensional operators obtained by rewriting the system in the Hilbert space M2 = Rn × L2 (−τ, 0; Rn ). In this case, Banks and Burns [3] applied a numerical method based on averaging approximations to open-loop control problems for delay systems. The convergence analysis was based on approximation results for linear semigroups involving the Trotter-Kato theorem. The averaging approximation scheme has been used in [4] for both the finite and infinite horizon linear quadratic problem of delay systems. Subsequent to the development of averaging methods for delay systems, Banks and Kappel [5] developed higher order approximation schemes based on spline approximations. A new spline approximation scheme has been developed in [6] for the linear quadratic problem of delay systems with any number of pure delay terms. Theoretical convergence results were obtained in the finite horizon case, as this approximation scheme does not guarantee the uniform

2 exponential stability of the approximate semigroups. The spline approximation scheme has also been applied to the linear quadratic problem of delay systems in [7]. Although numerical simulations show better performance than the averaging scheme, no theoretical convergence results are so far available. In [8], a piecewise linear approximation theory has been developed for the finite and infinite horizon linear quadratic problem of general delay systems. Theoretical convergence results were obtained both in the finite and infinite horizon cases, as the condition of uniform exponential stability is verified. Although the approximation schemes mentioned above are attractive from theoretical aspects, they may deal with the large-scale computations when a high order accuracy is required. Concerning the solution of TDOCPs, some other alternative techniques have been developed in the literature. Eller et al. [9] have proposed a method that involves the solution of a set of successively coupled partial differential equations. This method of solution is computationally tedious in general. The sensitivity approach [10] is another method which is only applicable for the systems with small delay in the state. The method presented in [11] is also a sensitivity approach in which the original system is embedded in a class of non-delay systems using an appropriate parameter. Ref. [12] has developed a computational scheme using the technique of control parameterization. This technique involves the approximation of the control variables by piecewise constant functions with the instants of possible switching pre-assigned. The TDOCP is then converted into a finite dimensional optimization problem which can be solved by a standard quasi-Newton method. Lee [13] has used the auxiliary state variables and a P´ ade approximation to transform the problem to one without a time-delay argument. The transformed problem can be solved by a well-developed optimization algorithm. Even the approximation errors can be reduced by adding more state variables, the computational times of the numerical integrations must be increased considerably. The next method for the solution of TDOCPs is an iterative dynamic programming (IDP) [14]. Owing to the backward solution manner, however, IDP still requires interpolating the data from the last iteration to estimate the historical information of time-delay states. To avoid the data interpolation, Chen et al. [15] used the Taylor expansion to estimate time-delay state variables under the framework of IDP. However, both procedures derived from IDP often rely on a large number of time and state grids to refine the solutions. Consequently, the enlarged problem size often leads to a slow convergence. The TDOCPs have also been solved using the orthogonal functions. This technique is based on converting the underlying differential equation into integral equation through integration, approximating various signals involved in the equation by the truncated orthogonal series and utilizing the operational matrix of integration, delay and product to reduce the problem to solve a system of algebraic equations. Typical examples have used the piecewise constant basis functions (e.g. Walsh functions [16] and block-pulse functions [17]) and orthogonal polynomials (e.g. Legendre polynomials [18] and Chebyshev polynomials [19]) for their expansions. However, the computed response of the delay systems via continuous or piecewise constant basis functions is not in good agreement with the exact response of the system [20]. The main difficulty arises because of the lack of smoothness of the optimal control for time-delay systems. Recently, the hybrid system of basis functions have been used which consist of the combination of block-pulse functions with Legendre polynomials [21], Bernoulli polynomials [22] and Chebyshev polynomials [23]. One of the most important advantages of hybrid functions is the good representation of smooth and especially nonsmooth functions. Although the computational accuracy can be enhanced by evaluating more components of hybrid series, considering more terms of hybrid series enlarges the size of the resultant system of algebraic equations. The system dimension is enlarged more by increasing the number of state and control variables. Consequently, the burden of computation increases due to the increase of system dimension. As it was mentioned, all the reviewed methods above suffer from the computational difficulties especially when a high order accuracy is required. The aim of this paper is to develop a computational efficient method for solving the OCP of linear time-delay systems with a quadratic cost functional.

3 In the proposed method, the application of Pontryagin’s maximum principle to the TDOCPs results in a system of coupled TPBVP involving both delay and advance terms [2]. Applying a successive substitution scheme, the coupled TPBVP is transformed into a sequence of linear time-invariant ordinary differential equations (ODEs) without delay and advance terms. The obtained optimal control consists of both feedback and forward portions. The feedback portion is obtained by solving a matrix Riccati differential equation. The forward portion is an infinite sum of adjoint vectors which can be obtained by solving recursively the above-mentioned sequence of linear non-delay ODEs. In this study, a fast-converging iterative algorithm is presented which provides a promising possible reduction of computational efforts. Numerical results indicate that just a few iterations are required to find an accurate enough feedforward-feedback suboptimal control. However, when a high order iteration is needed, using longer CPU time will become unavoidable. The paper is organized as follows. Section 2 describes the problem statement. In Section 3, the method of successive substitution is employed to propose an open-loop optimal control design strategy. The convergence analysis of the proposed approach is also provided in this section. Section 4 suggests a novel feedforward-feedback optimal control design strategy. In this section, an iterative algorithm and a discussion concerning the main feature of this algorithm are also given. The effectiveness of the proposed approach is verified by solving four numerical examples in Section 5. Finally, conclusions and future works are given in the last section.

2

Problem statement

Consider a linear system with state time-delay described by: ( x(t) ˙ = Ax(t) + A1 x(t − τ ) + Bu(t), x(t) = φ(t), −τ ≤ t ≤ 0,

0 ≤ t ≤ tf ,

(1)

where x(t) ∈ Rn and u(t) ∈ Rm are the state and control vectors, respectively, τ > 0 is the constant time-delay, A, A1 and B are real constant matrices of appropriate dimensions, and φ(t) ∈ Rn is a continuous initial state function on [−τ, 0]. In addition, it is assumed that the pair (A + A1 , B) is controllable [24]. Given φ ∈ C ([−τ, 0]; Rn ) and u ∈ L2 (0, tf ; Rm ), there exists a unique solution x(t; φ, u) for Eq. (1) that is absolutely continuous on [0, tf ] [25]. The objective is to find the optimal control law u∗ (t) over t ∈ [0, tf ] which minimizes the following quadratic cost functional: 1 1 J = xT (tf )Qf x(tf ) + 2 2

Z

tf 0


xT (t)Qx(t) + uT (t)Ru(t) dt,

(2)

subject to the system (1) where Qf ∈ Rn×n and Q ∈ Rn×n are positive semi-definite matrices and R ∈ Rm×m is a positive definite matrix. According to the Pontryagin’s maximum principle of OCPs with time-delay [2], the necessary conditions of optimality can be written as:   x(t) ˙ = Ax(t) + A1 x(t − τ ) − BR−1 B T λ(t), 0 ≤ t ≤ tf ,    λ(t) ˙ = −Qx(t) − AT λ(t) + A2 (t)λ(t + τ ), 0 ≤ t ≤ tf , (3)  x(t) = φ(t), −τ ≤ t ≤ 0,    λ(t ) = Q x(t ), f f f

where λ(t) ∈ Rn is the co-state vector, A2 (t) = −χ[t0 ,tf −τ ] (t)AT1 and χ[t0 ,tf −τ ] (t) is the characteristic

4 function on [t0 , tf − τ ] defined by: χ[t0 ,tf −τ ] (t) =

(

1, 0,

t ∈ [t0 , tf − τ ], otherwise.

(4)

Note that in the case of a quadratic functional, the Pontryagin’s maximum principle is also a sufficient condition. Also, the optimal control law is given by: u∗ (t) = −R−1 B T λ(t),

0 ≤ t ≤ tf .

(5)

As it is shown the problem (3) is a system of coupled TPBVP involving both delay and advance terms. Unfortunately, this problem cannot be solved analytically except in very special cases. To overcome this difficulty, a successive substitution method will be presented in the next section.

3

Open-loop optimal control design

In this section, first we apply a method of successive substitution in order to construct a sequence of approximations that converges to the solution of TPBVP (3). Then we prove the convergence of the proposed approach. For this purpose, consider the following initial value problem (IVP) corresponding to the original TPBVP (3): " # " # " #  x(t) A x(t − τ ) x(t) ˙  1  =M + , 0 ≤ t ≤ tf ,   λ(t) ˙ λ(t) A2 (t)λ(t + τ ) (6)  x(t) = φ(t), −τ ≤ t ≤ 0,    λ(0) = ζ,

where M is a 2n × 2n real constant matrix defined by:   A −BR−1 B T M, , −Q −AT

(7)

and ζ ∈ Rn is an unknown parameter. This parameter will be identified using the boundary condition λ(tf ) = Qf x(tf ). The solution of IVP (6) (and equivalently TPBVP (3)) can be expressed in the series P P∞ ˜ ˜ i (t) will be determined expansion forms x(t) = ∞ x ˜ (t) and λ(t) = λ (t) in which x ˜i (t) and λ i=0 i i=0 i recursively, as discussed hereinafter. To this end, we construct the following sequence of IVPs by applying a successive substitution scheme: " # # " ˙ 0 (t)  x ˜ x ˜ (t)  0  , 0 ≤ t ≤ tf , =M   ˜˙ ˜ 0 (t) λ λ0 (t) (8)  x ˜ (t) = φ(t), −τ ≤ t ≤ 0,  0   λ ˜ 0 (0) = ζ, " # " # # " ˙ i (t)  x ˜ A x ˜ (t − τ ) x ˜ (t)  1 i−1 i  + , 0 ≤ t ≤ tf , =M  ˙  ˜ i (t) ˜ i−1 (t + τ ) ˜ i (t) λ A2 (t)λ λ (9)  x ˜ (t) = 0, −τ ≤ t ≤ 0,  i   λ ˜ i (0) = 0.

Note that, Eq. (8) denotes the zeroth problem which has no delay and advance terms. Also, Eq. (9) denotes the ith problem, i ≥ 1, whose delay and advance terms are no longer considered as unknowns

5 since they have been obtained from the (i − 1)th problem. Therefore, solving the presented sequence is a recursive process. From Eqs. (8)-(9) we can construct the following solutions sequence: # " # "  φ(0) ˜0 (t)  x = S(t) , 0 ≤ t ≤ tf , ˜ 0 (t) (10) λ ζ   x ˜0 (t) = φ(t), −τ ≤ t ≤ 0, # " # Z " t  A x ˜ (θ − τ ) x ˜ (t) 1 i−1 i  dθ, 0 ≤ t ≤ tf , S(t − θ) = ˜ i−1 (θ + τ ) ˜ i (t) (11) A2 (θ)λ λ 0   x ˜i (t) = 0, −τ ≤ t ≤ 0,

where S(t) is the state transient matrix of M defined by S(t) , eM t .

Theorem 3.1. (Convergence analysis) Consider the vector functions sequences {xk (t)}∞ k=0 and {λk (t)}∞ as follows: k=0 k k X X ˜ i (t), λ (12) x ˜i (t), λk (t) , xk (t) , i=0

i=0

˜ i (t) are determined by Eqs. (10)-(11), recursively. Then sequences {xk (t)}∞ and where x ˜i (t) and λ k=0 {λk (t)}∞ k=0 converge uniformly to the solutions of IVP (6).

Proof. By substituting Eqs. (10)-(11) into expressions in (12), we obtain the following sequence of approximations # " # "   x0 (t) = S(t) φ(0) , 0 ≤ t ≤ t , f (13) λ0 (t) ζ   x0 (t) = φ(t), −τ ≤ t ≤ 0, " # # " # Z " t  A x (θ − τ ) x (t) x (t) 1 k−1 0 k  S(t − θ) dθ, 0 ≤ t ≤ tf , = + (14) A (θ)λ λk (t) λ0 (t) 0 2 k−1 (θ + τ )   xk (t) = φ(t), −τ ≤ t ≤ 0.

Define α ,

sup kS(t − θ)k, β , sup kφ(t)k, γ , kζk, η , max(β, γ) and κ , kA1 k. Notice that

t,θ∈[0,tf ]

t∈[−τ,0]

α ≥ 1, and sup kA2 (θ)k = κ since A2 (θ) = −χ[t0 ,tf −τ ] (θ)AT1 . Therefore, Eq. (13) implies that θ∈[0,tf ]

kx0 (t)k ≤ αη, kλ0 (t)k ≤ αη,

−τ ≤ t ≤ tf ,

(15)

0 ≤ t ≤ tf + τ.

(16)

From Eqs. (14) we have       Z t A1 x0 (θ − τ ) x1 (t) x0 (t) S(t − θ) dθ, − = A2 (θ)λ0 (θ + τ ) λ1 (t) λ0 (t) 0

0 ≤ t ≤ tf .

Then, with the aid of Eqs. (15)-(17), we obtain

      Z t  

x1 (t)

x0 (θ − τ ) x0 (t)



λ1 (t) − λ0 (t) ≤ ακ

λ0 (θ + τ ) dθ 0

    Z t  

x0 (θ − τ ) 0

dθ

+

≤ ακ

λ0 (θ + τ ) 0 0 ≤ 2α2 κηt, 0 ≤ t ≤ tf .

(17)

(18)

6 Also, Eq. (14) yields       Z t x2 (t) A1 (x1 (θ − τ ) − x0 (θ − τ )) x1 (t) S(t − θ) − dθ, = λ2 (t) A2 (θ)(λ1 (θ + τ ) − λ0 (θ + τ )) λ1 (t) 0

0 ≤ t ≤ tf ,

and then, with the aid of Eq. (18), we obtain

      Z t  

x2 (t)

x1 (θ − τ ) − x0 (θ − τ ) x1 (t)



λ2 (t) − λ1 (t) ≤ ακ

λ1 (θ + τ ) − λ0 (θ + τ ) dθ 0   Z t  

x1 (θ − τ ) − x0 (θ − τ ) 0

+ ≤ ακ

0 λ1 (θ + τ ) − λ0 (θ + τ ) 0 2 t ≤ 22 α3 κ2 η , 0 ≤ t ≤ tf . 2! Hence, by the mathematical induction

  

xk (t) xk−1 (t)

λk (t) − λk−1 (t)

 k

≤ 2k αk+1 κk η t ,

k!

 

dθ

0 ≤ t ≤ tf .

(19)

(20)

(21)

Applying the triangle inequality, for any r, we have

  

xk+r (t) xk (t)

λk+r (t) − λk (t) 

 k+r i k+1 X

i i+1 i t k+1 t

≤ 2 α κ η ≤ αη(2ακ) e2ακt .

i! (k + 1)!

(22)

i=k+1

∞ xk (t) Therefore the sequence is uniformly convergent. Taking the limit in Eq. (14) when λk (t) k=0 k −→ ∞ shows that the limit of this sequence is the solution of IVP (6).  As it is shown in Theorem 3.1, for a given co-state initial value ζ, the solution sequences {xk (t)}∞ k=0 ∗ denotes the value of ζ correand {λk (t)}∞ converge uniformly to the solutions of IVP (6). Let ζ k=0 sponding to the boundary condition λ(tf ) = Qf x(tf ). Once the value of ζ ∗ is known, the open-loop optimal control u∗ (t) and the corresponding state trajectory x∗ (t) of the TDOCP (1)-(2) are determined as: ∞ X ˜ i (t), 0 ≤ t ≤ tf , λ (23) u∗ (t) = −R−1 B T i=0

x∗ (t) =

∞ X i=0

x ˜i (t),

0 ≤ t ≤ tf ,

(24)

˜ i (t) for i ≥ 0 are determined recursively from Eqs. (10)-(11). However, the value of where x ˜i (t) and λ ∗ ˜ i (t) in Eqs. (10)-(11) ζ is not known in general and thus the vector functions sequences x ˜i (t) and λ depend on both t and the unknown parameter ζ, i.e. we have x ˜i (t, ζ) and x ˜i (t, ζ). To overcome this difficulty in practical applications, after sufficient iterations of the proposed approach we consider the boundary condition λM (tf , ζˆM ) = Qf xM (tf , ζˆM ) where M is a finite positive integer and ζˆM is an appropriate approximation for ζ ∗ . Thus ζˆM should be the real root of the following system of linear algebraic equations M M X X ˜ i (tf , ζˆM ) = Qf x ˜i (tf , ζˆM ). (25) λ i=0

i=0

Following this procedure we have ζˆM → ζ ∗ as M → ∞.

7 Remark 3.1. (Extension) The results of this section can be extended directly to the case of multipledelays when the state equation (1) is of the form  p X  x(t) Ai x(t − τi ) + Bu(t), 0 ≤ t ≤ tf , ˙ = Ax(t) + (26) i=1   x(t) = φ(t), −τmax ≤ t ≤ 0, where τi is a constant positive scalar, τmax = max(τi ) and Ai is a real constant matrix of appropriate i

dimensions. In this case, the proof of Theorem 3.1 proceed as before, except that the parameter κ must be replaced by pµ where µ , max(kAi k). The results can also be extended to the distributed delay i

systems when the state equation (1) is of the form  Z 0  x(t) A(s)x(t + s)ds + Bu(t), ˙ = Ax(t) +  x(t) = φ(t),

−τ

0 ≤ t ≤ tf ,

(27)

−τ ≤ t ≤ 0,

where A1 ∈ L2 (−τ, 0; Rn×n ). In this case, the proof of Theorem 3.1 proceed as before, except that the parameter κ must be replaced by τ µ where µ , sup (kA(s)k). s∈[−τ,0]

4

Feedforward-feedback optimal control design

In this section, the results of previous sections are modified in order to design an optimal control with both feedback and forward portions. The optimal control in Eq. (5) can be implemented as a closed-loop optimal control if the co-state vector obtained consists of a linear function of the state vector and a nonlinear term in the form λ(t) = K(t)x(t) + g(t),

(28)

where K(t) is a n × n real symmetric positive semi-definite matrix and g ∈ Rn is the adjoint vector. To show the existence of this relation, we express the optimality conditions given by Eq. (3) in the following form " " # # " # Z tf  A x(θ − τ ) x(t ) x(t)  1 f  S(tf − θ) dθ, 0 ≤ t ≤ tf ,   λ(t ) = S(tf − t) λ(t) + A2 (θ)λ(θ + τ ) t f (29)  x(t) = φ(t), −τ ≤ t ≤ 0,    λ(t ) = Q x(t ). f

f

f

If the transition matrix is partitioned

S(t) =



S11 (t) S12 (t) S21 (t) S22 (t)



(30)

,

where Sij is a n × n matrix and the integral is replaced by 2n × 1 vector

written

x(tf ) = S11 (tf − t)x(t) + S12 (tf − t)λ(t) + f1 (t),

λ(tf ) = S21 (tf − t)x(t) + S22 (tf − t)λ(t) + f2 (t).



 f1 (t) , Eq. (29) can be f2 (t) (31) (32)

8 Replacing x(tf ) and λ(tf ) from Eqs. (31)-(32) into the boundary condition λ(tf ) = Qf x(tf ) we obtain S21 (tf − t)x(t) + S22 (tf − t)λ(t) + f2 (t) = Qf (S11 (tf − t)x(t) + S12 (tf − t)λ(t) + f1 (t))

(33)

λ(t) = [S22 (tf − t) − Qf S12 (tf − t)]−1 [Qf S11 (tf − t) − S21 (tf − t)] x(t) + [S22 (tf − t) − Qf S12 (tf − t)]−1 [Qf f1 (t) − f2 (t)] , K(t)x(t) + g(t).

(34)

Solving for λ(t) yields

The definition of K(t) and g(t) are aparent by inspection of Eq. (34). Kalman [26] has shown that the required inverse exists for all t ∈ [0, tf ]. However, there is an alternative computational route to obtain K(t) and g(t). Let us begin with Eq. (28). Differentiating both sides with respect to t we get ˙ ˙ λ(t) = K(t)x(t) ˙ + K(t)x(t) + g(t). ˙

(35)

˙ Substituting x(t) ˙ and λ(t) from Eq. (3) and using Eq. (28) to eliminate λ(t) we obtain h i ˙ K(t) + K(t)A + AT K(t) − K(t)BR−1 B T K(t) + Q x(t)  + g(t) ˙ + AT g(t) − K(t)BR−1 B T g(t) −A2 (t)g(t + τ ) − A2 (t)K(t + τ )x(t + τ ) + K(t)A1 x(t − τ )] = 0.

(36)

Because this must be satisfied for all x(t), we obtain the following matrix Riccati differential equation ( ˙ K(t) = −K(t)A − AT K(t) + K(t)BR−1 B T K(t) − Q, (37) K(tf ) = Qf , and the adjoint vector differential equation ( g(t) ˙ = −(A − BR−1 B T K(t))T g(t) + A2 (t)g(t + τ ) + A2 (t)K(t + τ )x(t + τ ) − K(t)A1 x(t − τ ), g(tf ) = 0.

(38) It is worth mentioning that g(t) can also be computed by using the recursive procedure in this paper. ˜ i (t, ζ) for i ≥ 0 have been determined from Eqs. (10)-(11). Then g(t, ζ) can Assume that x ˜i (t, ζ) and λ be easily computed as: g(t, ζ) =

∞ X i=0

˜ i (t, ζ) − K(t) λ

∞ X

x ˜i (t, ζ)

∞ ∞   X X ˜ i (t, ζ) − K(t)˜ g˜i (t, ζ), λ xi (t, ζ) , = i=0

i=0

(39)

i=0

i ≥ 0,

˜ i (t, ζ) − K(t)˜ where g˜i (t, ζ) = λ xi (t, ζ) is the ith -order adjoint vector and ζ should be evaluated such that the boundary condition λ(tf ) = Qf x(tf ) is satisfied. Finally, the closed-loop optimal control law becomes: ∞ X u∗ (t, ζ) = −R−1 B T K(t)x(t) −R−1 B T g˜i (t, ζ), 0 < t < tf . (40) | {z } i=0 feedback | {z } feedforward

u∗ (t, ζ)

From Eq. (40) it is observed that consists both feedback and forward portions. Note that, the feedback term does not need to be computed by the proposed approach since x(t) is the available state trajectory of the closed-loop system. However, the forward term depends on the unknown parameter ζ which should be evaluated using the boundary condition λ(tf ) = Qf x(tf ).

9 Remark 4.1. (Infinite-time horizon case) The procedure assumes the final time to be finite while in practical control systems we may consider the case of tf → ∞. In this case, the following matrix Riccati algebraic equation is solved instead of the matrix Riccati differential equations (37) −KA − AT K + KBR−1 B T K − Q = 0,

(41)

where K is a n × n symmetric positive semi-definite matrix. However, we were not able to prove the convergence of the solutions sequence in Eqs. (13)-(14) for the infinite-time horizon problem. The reason is that the time interval tends to infinity and the supremum-norm of the transition matrix may become unbounded. Future work can be focused on extending the method presented in this paper to construct a uniform convergent solutions sequence for the infinite-time linear quadratic problem of time-delay systems.

4.1

Feedforward-feedback suboptimal control design

It is almost impossible to obtain the optimal control law as in Eq. (40) since Eq. (40) contains infinite series. Moreover, the parameter ζ corresponding to the boundary condition λ(tf ) = Qf x(tf ), i.e. ζ ∗ , is not known in general. Therefore, in practical applications, the feedback portion remains exact, and only the forward portion is approximated by replacing ∞ with a finite positive integer M . In addition, the unknown parameter ζ ∗ is approximated by ζˆM where ζˆM is the solution of the system of linear algebraic equations (25). Thus, the M th order suboptimal control is obtained as: uM (t, ζˆM ) =

 −1 T  −R B K(t)x(t),

M = 0, M X −1 B T K(t)x(t) − R−1 B T  g˜i (t, ζˆM ), −R  i=1

M ≥ 1.

(42)

The integer M in Eq. (42) is generally determined according to a concrete control precision. For example, the M th order suboptimal control in Eq. (42) has the desirable accuracy if for a given small enough constant ǫ > 0, the following condition holds:

    ˆM )

xM (t, ζˆM ) x (t, ζ M −1

< ǫ, eM , sup (43)

λ (t, ζˆ ) − λ

ˆ M M M −1 (t, ζM ) t∈[0,tf ] where eM is a supremum-norm error. Then the M th -order cost functional value JM is computed as Z  1 tf  T 1 T x (t)Qx(t) + uTM (t, ζˆM )RuM (t, ζˆM ) dt, (44) JM = x (tf )Qf x(tf ) + 2 2 t0

where x(t) is the closed-loop state trajectory obtained by applying uM (t, ζˆM ) to the original time-delay system (1). Notice that, if the values of M − 1 and M provide similar results, thus generating a stop, then according to Eq. (14), we have xk−1 = xk and λk−1 = λk for all k > M . Therefore, increasing M the value of JM does not change. Here, we give an iterative algorithm to design an accurate enough feedforward-feedback suboptimal control. The main feature of the algorithm is then discussed. 4.1.1

Algorithm: Designing the feedforward-feedback suboptimal control ˜ 0 (t, ζ) from Step 1. Obtain K(t) from Eq. (37) and also the zeroth-order terms x ˜0 (t, ζ) and λ Eq. (8). Then, compute the zeroth-order adjoint vector g˜0 (t, ζ) from Eq. (39).

10 Step 2. Let i = 1. ˜ i (t, ζ) from Eq. (9). Then, compute the ith -order Step 3. Obtain the ith -order terms x ˜i (t, ζ) and λ adjoint vector g˜i (t, ζ) from Eq. (39). Step 4. Let M = i. Compute ζˆM from Eq. (25). If the condition (43) holds for the given small enough constant ǫ > 0, go to step 5; else, replace i by i + 1 and go to step 3. Step 5. uM (t, ζˆM ) is the desirable feedforward-feedback suboptimal control. 4.1.2

Discussion

In this section, a discussion concerning the main feature of the proposed algorithm is given. This algorithm enables us to efficiently solve the TDOCP (1)-(2) by using a method of successive substitution. Through the algorithm, with the aid of Eqs. (8)-(9), a sequence of nonhomogeneous linear time-invariant non-delay ordinary differential equations, which have the same size as the original TPBVP (3), is solved in a recursive manner. Moreover, the homogeneous parts of these equations at different iterations are the same and only the nonhomogeneous parts are updated, recursively. These facts provide a promising possible reduction of computational efforts. In addition, the proposed algorithm converges extremely fast because, as seen from Eq. (21), the sequences kxk (t) − xk−1 (t)k and rk . Furtherkλk (t) − λk−1 (t)k are dominated by a fast uniformly convergent sequence of the form k! more, if the tolerance error bound ǫ is chosen small enough, the approximation ζˆM of the co-state initial value will be very close to its optimal value ζ ∗ and thus the boundary condition at tf will be satisfied tightly.

5

Numerical examples

In this section, four numerical examples are given to demonstrate the effectiveness of the proposed approach. In these examples, the results of our suggested technique are compared with the existing results. The numerical results confirm the theoretical results in the case of finite-time horizon problem. However, our experience has indicated that the scheme performs also well in the case of infinite-time horizon problem. This is shown by an example which already has been considered in [6–8]. Example 5.1. This example was first introduced by Eller et al. [9], and has been subsequently considered by Jamshidi and Malek-Zavarei [11], Oh and Luus [27], Dadebo and Luus [14], and more recently by Santos and S´ anchez-D´ıaz [28], to illustrate their methods. Consider a first order linear time-delay system described by: ( x(t) ˙ = x(t) + x(t − 1) + u(t), 0 ≤ t ≤ 2, (45) x(t) = 1, −1 ≤ t ≤ 0, with the quadratic cost functional: J=

Z

0

2


x2 (t) + u2 (t) dt.

(46)

The problem is to find the optimal control u∗ (t) which minimizes expression (46) subject to the timedelay system (45).

11 Here, we apply the proposed approach to solve this problem. From Eqs. (45)-(46) we have A = 1, A1 = 1, B = 1, Qf = 0 and Q = R = 2. Therefore, the matrix Riccati differential equation (37) becomes: ( ˙ K(t) + 2K(t) − 21 K 2 (t) + 2 = 0, (47) K(2) = 0, whose unique solution is: K(t) = −

√

0.1182114931e

2t

0.1426938951e

2t

√

− 33.83765733e− − 7.008008292e−

√

2t

√

2t

(48)

.

Following the recursive procedure described in Section 3, the zeroth IVP (8) becomes:  ˜ 0 (t), 0 ≤ t ≤ 2,  x ˜˙ 0 (t) = x ˜0 (t) − 0.5λ     ˜˙ ˜ 0 (t), 0 ≤ t ≤ 2, λ0 (t) = −2˜ x(0) (t) − λ  x ˜0 (0) = 1,    λ ˜ 0 (0) = ζ.

(49)

From Eq. (49) we have:

√

x ˜0 (t) = 0.125(1.171572876 + 1.414213562ζ)e− √2t −0.125(−6.828427124 + 1.414213562ζ)e 2t ,

(50)

0 ≤ t ≤ 2,

√

˜ 0 (t) = 0.25(3.414213562ζ + 2.828427124)e− 2t λ √ −0.25(−0.585786438ζ + 2.828427124)e 2t ,

(51)

0 ≤ t ≤ 2.

Substituting Eqs. (50)-(51) into Eq. (9) for i = 1 and noting that x ˜0 (t) = 1, −1 ≤ t ≤ 0, the first IVP is obtained whose solution is in the form: √

− 2t x ˜1 (t) = −0.5 + ((0.03668348235t − 0.02467865367)ζ + 0.08897265484 + 0.03038959181t)e √ (52) ((0.02467865365 − 0.1064852556t)ζ + 0.4110273452 + 0.5141562972t)e 2t , 0 ≤ t ≤ 1, √

x ˜1 (t) = ((0.11042785 + 0.02588834764t)ζ − 0.7061417822 + 0.02144660953t)e− √2t ((−0.4440003117 − 0.1508883476t)ζ + 4.040672708 + 0.7285533906t)e 2t ,

1 ≤ t ≤ 2, √

˜ 1 (t) = 1 + ((0.1771235212t − 0.1925260453)ζ + 0.1467339294t + 0.3688187959)e− 2t λ √ ((0.1925260453 + 0.088215274t)ζ − 1.368818796 − 0.4259410231t)e 2t , 0 ≤ t ≤ 1, √

˜ 1 (t) = ((0.8349695209 + 0.125t)ζ − 3.159554134 + 0.1035533906t)e− 2t λ √ ((0.3160452081 + 0.125t)ζ − 3.097402872 − 0.6035533905t)e 2t ,

1 ≤ t ≤ 2,

(53)

(54)

(55)

and so on. If we select M = 2, then the co-state initial value and the minimum value of the cost functional (46) are determined as ζˆ2 = 8.1684022784 and J2 = 6.2196152583, respectively, after the second step of the iterative algorithm. This value of J2 compares well with those given by Eller et al. [9], Jamshidi and Malek-Zavarei [11], Oh and Luus [27], Dadebo and Luus [14] and Santos and S´ anchez-D´ıaz [28]. Comparison results are listed in Table 1. Simulation curves of u2 (t) and the corresponding closed-loop state trajectory are shown in figure 1.

12

Table 1: The cost functional values for different schemes (Example 5.1). Method Eller et al. [9] Jamshidi and Malek-Zavarei [11] Oh and Luus [27] Dadebo and Luus [14] Santos and S´ anchez-D´ıaz [28] The proposed method (M = 2)

Cost functional value J 6.45 6.5 6.23711 6.26775 6.97 6.2196152583

0

1

−0.5

0.9

Closed−loop trajectory x(t)

−1

0.8

−1.5

x (t)

u(t)

0.7 −2 −2.5

0.6 0.5

−3 0.4

−3.5 −4 −4.5

0.3

Suboptimal control u2(t) 0

0.5

1

1.5

2

0.2

0

0.5

t (sec.)

(a) The control variable u(t).

1

1.5

2

t (sec.)

(b) The state variable x(t).

Figure 1: Second-order suboptimal control and the corresponding closed-loop state trajectory for for Exam-

ple 5.1.

Example 5.2. This example was considered in [3] and also studied in [13, 22, 29]. Consider the optimal control of a harmonic oscillator with retarded damping. It involves the minimization of the cost functional: Z 1 2 2 u (t)dt, (56) J = 5x21 (2) + 2 0 with the second order linear time-delay system:   x˙ 1 (t) = x2 (t), 0 ≤ t ≤ 2,    x˙ (t) = −x (t) − x (t − 1) + u(t), 2 1 2 x1 (t) = 10, −1 ≤ t ≤ 0,    x (t) = 0, −1 ≤ t ≤ 0. 2

0 ≤ t ≤ 2,

(57)

The problem is to find the optimal control u∗ (t) which minimizes expression (56) subject to the timedelay system (57). The exact solution for this problem is given by [3]: ( δ sin(2 − t) + ( 2δ )(1 − t) sin(t − 1), u∗ (t) = δ sin(2 − t), 1 ≤ t ≤ 2,

0 ≤ t < 1,

(58)

where δ = 2.5599 and the optimal cost is J ∗ = 3.3991 [3]. Here, we solve this problem by using the suggested algorithm with the tolerance error bound ǫ = 0.02. Simulation results at different

13 iteration times are summarized in Table 2. From Table 2 it is observed that the convergence is achieved after only four iterations, i.e. e2 = 0.0121918086 < ǫ. In addition, with the used precision, the co-state initial value is determined as ζˆ4 = [0.6794209611, −1.2497991708] T and a minimum value of J4 = 3.3991221016 is obtained in the fourth iteration. This value of J4 compares very well with those achieved by using a semi-group approximation scheme (when the problem is considered as an infinite dimensional OCP in the space M2 ) [3] (J = 3.4827), combined parameter and function optimization algorithm [13] (J = 3.4827) and eight basis functions of linear Legendre multi-wavelets [29] (J = 3.43254). Notice that, a minimum value of J = 3.21663 has also been reported in [22]. In that work, first the original OCP is approximated by using a hybrid function approximation method. Then, the cost is calculated for the approximate problem. Therefore, the reported value is the optimal cost for the approximate problem and is not a true value for the original problem. Figure 2 shows the fourth order suboptimal control u4 (t) and the corresponding closed-loop state trajectories, obtained by four iterations of our suggested algorithm, which compare very well with the optimal solutions. Table 2: Simulation results of Example 5.2 at different iteration times. i (iteration time) 0 1 2 3 4

Cost functional value Ji 6.4755994156 3.4747645300 3.3998363195 3.3991653665 3.3991221016

2.5

Supremum-norm error ei 2.7908489582 0.3721008386 0.0547036367 0.0121918086

10

0

Closed−loop trajectory x (t) 2

−1

Optimal trajectory x*(t)

8

2

2

−2 6

x 1 (t)

u(t) 1

x 2 (t)

−3

1.5

4

−4 −5

2 0.5

Suboptimal control u4(t) 0

0.5

1

−7

*

Optimal trajectory x1(t)

Optimal control u*(t) 0

−6

Closed−loop trajectory x1(t)

0

1.5

2

t (sec.)

(a) The control variable u(t).

−2

0

0.5

1

1.5

2

t (sec.)

(b) The state variable x1 (t).

−8

0

0.5

1

1.5

2

t (sec.)

(c) The state variable x2 (t).

Figure 2: Fourth-order suboptimal control and the corresponding closed-loop state trajectories for Example 5.2.

Example 5.3. This example was already considered in [8] where a piecewise linear approximation scheme was used for computing the infinite dimensional Riccati operator in regulator problems for time-delay systems. Minimize: Z
1 2 2
1 2 2 x1 (2) + x2 (2) + J= u1 (t) + u22 (t) dt, (59) 2 2 0

14

Table 3: Simulation results of Example 5.3 at different iteration times. i (iteration time) 0 1 2 3 4 5 6 7 8 9

subject to:

Cost functional value Ji 0.3333333333 1.5187087597 1.5395777363 1.3959807371 1.3905067209 1.4022417517 1.4026747456 1.4016881515 1.4016516448 1.4017343910

Supremum-norm error ei 0.6111120521 0.5112948346 0.0539492160 0.0428904816 0.0045255861 0.0035979142 0.0003796343 0.0003018154 0.0000318460

# # " "  1 0 0 0    u(t), x(t − 1) + x(t) ˙ =   0 1 1 0   " #    1    , x(t) = 1

0 ≤ t ≤ 2, (60)

−1 ≤ t ≤ 0.

The true solution of this problem was given in [30]. Here, we applied the iterative algorithm expressed in Section 4.1.1 to solve this problem. The resulting values for the cost functional Ji and the supremum-norm error ei computed for different values of i are given in Table 3. As it is shown in Table 3, for the tolerance error bound ǫ = 0.0001 the convergence is achieved after nine iterations, i.e. e9 = 0.0000318460 < ǫ. In this case, the co-state initial value is determined as ζˆ9 = [1.0623381790, 0.8713547895] T . Table 3 also indicates that the error is reduced further by increasing the iteration time of the algorithm. Moreover, a minimum value of J9 = 1.4017343910 is obtained after the ninth step of the iterative algorithm. This value of J9 compares very well with the true solution [30] (J ∗ = 1.4017) and the value computed by the piecewise linear approximation scheme [8] (J = 1.4017). Simulation curves of the ninth-order suboptimal controls and the corresponding closed-loop state trajectories are depicted in figure 2 which are in good agreement with the available literature [8, 30]. Example 5.4. This example is an infinite-time horizon linear quadratic problem for time-delay systems which has been investigated in the literature from the viewpoint of approximating the infinite dimensional Riccati operator obtained by rewriting the system in the Hilbert space M2 [6–8]. This is the problem of minimizing: Z ∞
J= x2 (t) + u2 (t) dt, (61) 0

subject to

( x(t) ˙ = x(t) + x(t − 1) + u(t), x(t) = sin(πt), −1 ≤ t ≤ 0.

t ≥ 0,

(62)

In order to apply the suggested algorithm in Section 4.1.1, we select tf = 20 and ǫ = 0.005. In this case, the used precision is satisfied after the seventh step of the iterative algorithm, i.e. e7 = 0.0036413610 < ǫ. In addition, the resulting value for the cost functional is J7 = 0.3217880340 which

15

−0.1

0.5

Ninth−order suboptimal control

−0.2 −0.3

0

−0.5

u 2 (t)

u 1 (t)

−0.4

−0.6 −0.7

−0.5

−1

−0.8 −0.9

−1.5

−1 −1.1

Ninth−order suboptimal control 0

0.5

1

1.5

−2

2

0

0.5

t (sec.)

1

1.5

2

t (sec.)

(a) The control variable u1 (t).

(b) The control variable u2 (t).

1

1.25

Closed−loop trajectory x1(t)

0.9

Closed−loop trajectory x2(t)

1.2

0.8

1.15

0.7

x 2 (t)

x 1 (t)

1.1 0.6 0.5

1.05 1

0.4 0.95

0.3 0.2

0.9

0.1

0.85

0

0.5

1

1.5

2

0

0.5

t (sec.)

1

1.5

2

t (sec.)

(c) The state variable x1 (t).

(d) The state variable x2 (t).

Figure 3: Ninth-order suboptimal control and the corresponding closed-loop state trajectories for Example 5.3. compares favorably with those computed by the approximation schemes in [6–8]. Despite the fact that the new scheme converges numerically in this example, we were not able to establish the convergence of the proposed scheme for the infinite-time horizon problems. Future work will involve the infinite-time horizon linear quadratic problem for time-delay systems, which in this paper has been only sketched. Example 5.5. In this example, we consider a simplified model for a wind tunnel at the NASA Langley Research Center in the form [31]       −a 0 0 0 ka 0 0 x(t) ˙ = 0 0 1  x(t) +  0 0 0  x(t − 0.33) +  0  u(t), t ≥ 0, (63) 0 −ω 2 −2ξω 0 0 0 ω2 

 −0.1 x(t) =  8.547  , −0.33 ≤ t ≤ 0, 0

(64)

1 where = 1.964, k = −0.0117, ξ = 0.8 and ω = 6. The problem is to find the optimal control u∗ (t) a which minimizes the following infinite-time horizon quadratic cost functional Z ∞
104 x21 (t) + u2 (t) dt, (65) J= 0

subject to the dynamic constraints (63)-(64).

16

Table 4: Simulation results of Example 5.5 at different iteration times. i (iteration time) 0 1 2 3 .. . 8 9 10 11 .. .

26 27 28 29 .. .

48 49 50 51

Cost functional value Ji 98.284846 159.502388 140.745792 130.592298 .. . 136.264027 137.114451 136.491533 136.032396 .. .

136.405982 136.397176 136.404168 136.409833 .. .

136.404857 136.404939 136.404873 136.404818

Supremum-norm error ei 0.38380329 0.13326576 0.07774956 .. . 0.00816955 0.00620229 0.00456378 0.00337520 .. .

0.00008006 0.00006456 0.00005125 0.00004153 .. .

0.00000073 0.00000060 0.00000048 0.00000039

CPU time (sec.) 0.498477 2.251606 4.046088 5.772680 .. . 15.256115 16.951540 18.611821 20.209162 .. .

40.244765 41.391497 42.527167 43.676102 .. .

56.912967 57.584903 58.254324 58.935752

This example has been investigated based on abstract operators approximation schemes provided in [6–8, 33]. The true solution for this problem has also been given in [32]. Here we select tf = 20 and apply the suggested algorithm in Section 4.1.1 to solve this problem. Simulation results including the cost functional value, supremum-norm error and the elapsed CPU time at different iteration times are listed in Table 4. From Table 4 it is observed that for M = 51 a minimum value of J51 = 136.404818 is achieved within a 58.935752 sec. of CPU time. This value of J51 compares very well with the true value (J = 136.40490) [32] and those computed by the approximation schemes in [6–8, 33]. However, comparing the CPU time with the methods presented in [6–8,33] cannot be done because such papers do not report the elapsed CPU time for this problem. Simulation curves of the suboptimal control u51 (t) and the corresponding closed-loop state trajectories are depicted in figure 4 which are in good agreement with the available literature. Remark 5.1. As it is shown in Tables 2-4, the supremum-norm error ei given by Eq. (43) decreases as the iteration time increases. This is due to the fact that, as seen from Eq. (21), the sequences kxk (t) − xk−1 (t)k and kλk (t) − λk−1 (t)k are dominated by a fast uniformly convergent sequence of rk . However, the proposed algorithm does not guarantee a monotone decreasing of the the form k! cost functional values Ji . This may provide the approximate cost functional assumes values at some iterations smaller than the optimal one (see for example the iteration times i = 0, 3, 8, 27 in Table 4). However, simulation results reveal that the proposed algorithm follows a decreasing manner for the values of |Ji − Ji−1 | as the iteration time increases.

17

0.5

0 −0.01

0 −0.02 −0.03

x 1 (t)

u(t)

−0.5

−1

−1.5

−0.04 −0.05 −0.06 −0.07

−2

−0.08 −2.5 −0.09

Suboptimal control u51(t) −3

0

5

10

15

−0.1

20

Closed−loop trajectory x1(t) 0

5

t (sec.)

(a) The control variable u(t).

8

0

6

−5

x 3 (t)

5

x 2 (t)

15

20

(b) The state variable x1 (t).

10

4

2

0

−10

−15

−20

−2

−25

Closed−loop trajectory x2(t) −4

10

t (sec.)

0

5

10

15

t (sec.)

(c) The state variable x2 (t).

Closed−loop trajectory x3(t) 20

−30

0

5

10

15

20

t (sec.)

(d) The state variable x3 (t).

Figure 4: Suboptimal control u51 (t) and the corresponding closed-loop state trajectories for Example 5.5.

6

Conclusion

In this paper, a new successive substitution scheme has been developed for the finite-time linear quadratic control of time-delay systems. The obtained optimal control consists of both feedback and forward portions. The feedback portion is obtained by solving a matrix Riccati differential equation. The forward portion is an infinite sum of adjoint vectors, which can be obtained by solving recursively a sequence of nonhomogeneous linear time-invariant non-delay ODEs. The results not only demonstrate the efficiency, simplicity, and high accuracy of the suggested technique, but also indicate its effectiveness in practical use. Future work can be focused on extending the present method for solving the infinite-time horizon linear quadratic problem for time-delay systems.

Acknowledgements The authors gratefully acknowledge the helpful comments and suggestions of the reviewers, which have improved the manuscript.

References [1] Malek-Zavarei, M., Jamshidi, M.: Time Delay Systems: Analysis, Opimization and Applications. North-Holland, Amesterdam (1987)

18 [2] Kharatishvili, G.L.: The maximum principle in the theory of optimal process with time-lags. Dokl. Akad. Nauk SSSR 136, 39-42 (1961) [3] Banks, H.T., Burns, J.A.: Hereditary control problem: numerical methods based on averaging approximations. SIAM J. Control Optim. 16, 169-208 (1978) [4] Gibson, J.S.: Linear-quadratic optimal control of hereditary differential systems: Infinite dimensional Riccati equations and numerical approximations. SIAM J. Control Optim. 21, 95-139 (1983) [5] Banks, H.T., Kappel, F.: Spline approximations for functional differential equations. J. Differ. Equations 34, 496-522 (1979) [6] Kappel, F., Salamon, D.: Spline approximation for retarded systems and the Riccati equation. SIAM J. Control Optim. 25, 1082-1117 (1987) [7] Banks, H.T., Rosen, G.I., Ito, K.: A spline based technique for computing Riccati operators and feedback controls in regulator problems for delay equations. SIAM J. Sci. Statist. Comput. 5, 830-855 (1984) [8] Propst, G.: Piecewise linear approximation for hereditary control problems. SIAM J. Control Optim. 28, 70-96 (1990) [9] Eller, D.H., Aggarwal, J.K., Banks, H.T.: Optimal control of linear time delay systems. IEEE Trans. Autom. Control 14, 678-687 (1969) [10] Inoue, K., Akashi, H., Ogino, K., Sawaragi, Y.: Sensitivity approaches to optimization of linear systems with time-delay. Automatica 17, 671-676 (1971) [11] Jamshidi, M., Malek-Zavarei, M.: Suboptimal design of linear control systems with time delay. Proc. IEE 12, 1743-1746 (1972) [12] Wong, K.H., Clements, D.J., Teo, K.L.: Optimal control computation for nonlinear time-Lag systems. J. Optim. Theory Appl. 47, 91-107 (1985) [13] Lee, A.Y.: Numerical solution of time-delayed optimal control problems with terminal inequality constraints. Optim. Control Appl. Meth. 14, 203-210 (1993) [14] Dadebo, S. A., Luus, R.: Optimal control of time-delay systems by dynamic programming. Optim. Control Appl. Meth. 13, 29-41 (1992) [15] Chen, C.L., Sun, D.Y., Chang, C.Y.: Numerical solution of time-delayed optimal control problems by iterative dynamic programming. Optim. Control Appl. Meth. 21, 91-105 (2000) [16] Pananisamy, K.R., Rao, G.P.: Optimal control of linear systems with delays in state and control via Wash functions. IEEE Proc. 130, 300-312 (1983) [17] Hwang, C., Shih, Y.P.: Optimal control of delay systems via block pulse functions. J. Optim. Theory Appl. 45, 101-112 (1985) [18] Perng, M.H.: Direct approach for the optimal control of linear time-delay systems via shifted Legendre polynomials. Int. J. Control 43, 1897-1904 (1986) [19] Horng, I.R., Chou, J.H.: Analysis, parameter estimation and optimal control of time-delay systems via Chebyshev series. Int. J. Control 41, 1221-1234 (1985)

19 [20] Datta, K.B., Mohan, B.M.: Orthogonal Functions in Systems and Control. World Scientific, Singapore (1995) [21] Marzban, H.R., Razzaghi, M.: Optimal control of linear delay systems via hybrid of block-pulse and Legendre polynomials. J. Franklin I. 341, 279-293 (2004) [22] Haddadi, N.,Ordokhani, Y., Razzaghi, M.: Optimal control of delay systems by using a hybrid functions approximation. J Optim. Theory Appl. 153, 338-356 (2012) [23] Marzban, H.R., Hoseini, S.M.: Solution of linear optimal control problems with time delay using a composite Chebyshev finite difference method. Optim. Control Appl. Meth. 34, 253-274 (2013) [24] Hewer, G.A.: A note on controllability of linear systems with time delay. IEEE Trans. Autom. Control 17, 733-734 (1972) [25] O˘ guzt¨ oreli, M.N.: Time-lag Control Systems. Academic Press, New York, 1966. [26] Kalman, R.E.: Contributions to the theory of optimal control. Bol. Soc. Mat. Mex. 5, 102-119 (1960) [27] Oh, S.E., Luus, R.: Optimal feedback control of time-delay systems. AIChE J. 22, 140-147 (1976) [28] Santos, O., S´ anchez-D´ıaz, G.: Suboptimal control based on hill-climbing method for time delay systems. IET Control Theory Appl. 5, 1441-1450 (2007) [29] Khellat, F.: Optimal control of linear time-delayed systems by linear Legendre multiwavelets. J. Optim. Theory Appl. 143, 107-121 (2009) [30] Banks, H.T., Burns, J.A., Cliff, E.M., Thrift, P.R.: Numerical Solutions of Hereditary Control Problems via an Approximation Technique, LCDS Technical Report 75-6, Brown University (1975) [31] Armstrong, E.S., Tripp, J.S.: An application of multivariable design techniques to the control of the National Transonic Facility, NASA Tech. Paper 1887, NASA Langley Research Center, Hampton, VA (1981) [32] Manitius, A., Tran, H.: Numerical simulation of a nonlinearfeedback controllerfor a wind tunnel model involving a time delay, Optimal Control Appl. Meth. 7, 19-39 (1986) [33] Germani, A., Manes, C., Pepe, P.: A Twofold Spline Approximation for Finite Horizon LQG Control of Hereditary Systems, SIAM J. Control Optim. 39, 1233-1295 (2000)