Optimal control of linear time-varying systems with state and input delays by Chebyshev wavelets

This paper presents a novel method for finding the optimal control, state and cost of linear time-delay systems with quadratic performance indices. The basic idea here is to convert a time-delay optimal control problem into a quadratic programming one which can be easily solved using MATLAB. To implement this idea we choose a state and control parameterization method by using Chebyshev wavelets. The inverse time operational matrix of Chebyshev wavelets is introduced and applied for parameterizing state and control terms containing inverse time. The method is also applicable to linear quadratic time-delay systems with combined constraints. Illustrative examples demonstrate the validity and applicability of the approach which new expansions for initial vector functions of state and control variables are defined.


Introduction
A time-delay (TD) system is a system in which time delays occur between the application of the delayed variables to the system and their resulting effect on it.The TD system can be constructed either with inherent delays in the components or with a deliberate application of time delays into the system for some purposes [1].Delays occur frequently in chemical processes, electronic, aerospace and mechanical systems, transmission lines and industrial processes.A mathematical model of a system such as population growth, epidemic growth, economic growth and neural networks results in a delay differential equation involving the variable time [2].The extension of Pontryagin Maximum Principle on optimal control system with delays, as described by Kharatishvili in [3], constitute a twopoint boundary value problem which its solution is usually very difficult due to the coupled nature of the solutions.In some papers [8,9,10,17], orthogonal functions and polynomials, such that block-pulse and Walsh functions, Chebyshev and Legendre polynomials, were applied to find the optimal control of TD systems.We can see the properties of Chebyshev polynomials in [12] produce an approximating function which minimizes error in its application.Chebyshev polynomials form a special class of polynomials especially suited for approximation theory.These polynomials as scaled Mother functions, are basis of Chebyshev wavelets [13].
In this paper, we introduce a direct method to solve the linear quadratic optimal control problems with delays and inverse time terms.The method consists of transforming the original optimal control problem to a quadratic programming (QP) one by using Chebyshev scaling functions and then solving the converted problem.We do not need a special program to solve the QP problem and many numerical methods are available to solve it.We can use

Definition and Function approximation
Chebyshev polynomial of the first kind T m (x) defined by T m (x) = cos(m arccos(x)) is a polynomial in x of degree m which satisfies the equation (1 − x 2 )T m − xT m + m 2 T m = 0. Chebyshev polynomials form an orthogonal polynomial set on the interval [−1 , 1] with respect to the weight w(x) = (1 − x 2 ) −1/2 .Some of the useful relations are Chebyshev wavelets of the first kind ψ nm defined on the interval [0, 1], have four arguments: k ∈ N >1 specifies the number of subintervals, n = 1, 2, . . ., 2 k−1 refers to the segment number, m = 0, 1, . . ., M − 1 is the degree of T m (x) and t is as an independent variable denotes the time; where ( They form an orthogonal basis with respect to the weight function w n (t), where w n (t) = w(2 k t − 2n + 1).A function f (t) defined over the interval [0, 1] can be approximated by Chebyshev wavelets as where f and Ψ(t) are 2 k−1 M × 1 column vectors, denotes transposition and Ψ(t) = [ψ 10 (t), ψ 11 (t), . . ., ψ 1M −1 (t), ψ 20 (t), . . ., ψ 2M −1 (t), . . ., ψ 2 k−1 0 (t), . . ., OPTIMAL CONTROL OF LINEAR TD SYSTEMS BY CHEBYSHEV WAVELETS We need to to find f; when With multiplying this equation by ψ n m (t)w n (t) and then integrating both sides from n−1 2 k−1 to n 2 k−1 , we find where m = 0, 1, . . ., M − 1, n = 1, 2, . . ., 2 k−1 ; the definition of Chebyshev wavelets in (1) demands that ψ nm (t)ψ n m (t) = 0 when n = n , (6) thus by replacing 2 k t − 2n + 1 by x in the right side for n = n we obtain Theorem 1 Assume that a function f (t) is a δ-Lipschitz function on [0, 1], then f (t) can be expanded as , where this series converges toward the function f (t).Theorem 2 (Error bound) Let f (t) be a δ-Lipschitz function on [0, 1], where f (t) = Then for M > 1, we have the following error bound Hence by substituting into formula ( 7) So, after integration by parts for m ≥ 1 and some ζ, we get By using the fact that k ≥ 2 and the norm definition of a function, we can write

The product operational matrix of Chebyshev wavelets
Let f be a row vector defined in (4).Then the product operational matrix of two Chebyshev wavelet vectors f is obtained from: as a 2 k−1 M × 2 k−1 M matrix and can be found by the compact support property (6) in the form of a block diagonal matrix as f = blkdiag( f1 , f2 , . . ., f2 k−1 ).(11) Assume fn = [ fn ab ] where a = m + 1 and b = m + 1.The M × M matrices fn are determined by using (1) and the recurrence relation of Chebyshev polynomials on the left side of (10), and then equating coefficients of same wavelet on both sides of (10).So from ψ n0 (t)ψ nm (t) = 2 k/2 √ π ψ nm (t) and the following findings: where and f nξ = f nM −1 , M even 0, M odd.

The delay operational matrix of Chebyshev wavelets
The delay Chebyshev scaling function Ψ(t − t d ) is the shifted function of Ψ(t) along the t-axis and it is found to be as where D is the 2 k−1 M × 2 k−1 M delay operational matrix of Chebyshev wavelets.We use (1) to compute D; Note that if t f = 1 then we let t/t f → t.From (1) it follows that when Consequently D is a sparse matrix and has the following structure: thus this matrix can be written compactly as Throughout this paper we use 0 and I to denote the zero and identity matrices, respectively.

The inverse time operational matrix of Chebyshev wavelets
Some linear time-varying systems comprise the expressions x(t f − t) and/or u(t f − t) in their state equations.For expanding these terms by Chebyshev scaling functions we can write (by setting t/t f → t) Stat., Optim.Inf.Comput.Vol. 5, December 2017 IMAN MALMIR 307 where we denote the inverse time operational matrix of Chebyshev wavelets by Υ.To find this matrix, we can write In other words, when n = 1, 2, . . ., in which Z is a M × M diagonal matrix as

The operational matrix of Chebyshev wavelets for integration
The integration of the Chebyshev wavelet vector on the interval [0, t] can be obtained as The matrix P is called the operational matrix of Chebyshev wavelets for integration, where Y and J are M × M matrices given by OPTIMAL CONTROL OF LINEAR TD SYSTEMS BY CHEBYSHEV WAVELETS 2.6.The integration matrix of the product operational matrix of Chebyshev wavelets We need to find the integration matrix of the product of two Chebyshev scaling function vectors on the time interval [0, 1].By defining and using the properties of Chebyshev polynomials, we deduce that L is a 2 k−1 M × 2 k−1 M symmetric matrix and can be obtained by integrating the entries of Ψ(t)Ψ (t) from 0 to 1.This may be done by introducing the auxiliary row vector φ n (t) = [ψ n0 (t), ψ n1 (t), . . ., ψ nM −1 (t)] and using (1) as follows, This yields We let L n = [l n ab ] M ×M , where a and b are the same as before.When m = m = 0, we simply have l 11 = 2 π ; otherwise, from (8) we see that where and we can write Clearly one can see that when m + m is odd, then l n ab = 0.It is apparent the values of the entries l n ab of L n depend not on n but only on m and m , thus Hence the integration matrix of the product of two Chebyshev wavelet vectors takes the form where For instance if M = 8, then 3. QP representation of an optimal control problem

Problem statement and Approximation process
In this section we use the above-mentioned results to model a linear quadratic TD control problem as a QP one.Consider a linear TD system with a quadratic performance index (PI) or cost functional where x(t) ∈ R q and u(t) ∈ R r are state and control vectors, E(t), F(t), G(t), H(t), V(t) and W(t) are piecewisecontinuous matrices of compatible dimensions, t dx and t du represent state and control delays, respectively, x 0 is an initial condition vector, α(t) and β(t) are, respectively, qand r-dimensional initial state and initial control vector functions, T and Q(t) are positive semi-definite matrices and R(t) is a positive definite matrix.The problem is to find u * (t), x * (t) and J * such that the PI in ( 29) is minimized while satisfying Eqs. ( 26)-( 28), where * indicates optimal value.First we approximate the system dynamics.Since Chebyshev wavelets are defined on [0, 1], it makes sense to define a new variable by = t/t f ; replacing t by , the state equation ( 26) becomes (30) where x = t dx /t f and u = t du /t f .By integrating (30) from 0 to and rearranging, we get (31) In the approximation process we will use the symbolsˆand¯placed on top of a matrix to denote Kronecker product of the matrix and I q and Kronecker product of the matrix and I r , respectively.We parameterize the state and control vectors as x( where X and U are 2 k−1 qM × 1 and 2 k−1 rM × 1 column vectors of unknown parameters and represented by Our immediate goal is to find X and U from the optimization.Using (3), the initial state can be expanded in terms of Chebyshev wavelets as where we have Stat., Optim.Inf.Comput.Vol. 5, December 2017 310 OPTIMAL CONTROL OF LINEAR TD SYSTEMS BY CHEBYSHEV WAVELETS Also for the matrices E( ) and G( ) we can write where [E nm ] and [G nm ] are q × 2 k−1 qM and q × 2 k−1 rM matrices of constants and we take To illustrate, suppose ) cos mθ dθ, i, j = 1, 2, . . ., q.

Compatibility constraint
The proposed method follows a procedure in which the interval [0, 1] is divided into 2 k−1 subintervals.To ensure the continuity of the obtained state functions on the segment interfaces (boundary points) which are given by we enforce the following compatibility constraint: where OPTIMAL CONTROL OF LINEAR TD SYSTEMS BY CHEBYSHEV WAVELETS By using the relations of Chebyshev polynomials we find (54)

TD optimal control problem reformulation
Now we are in a position to generalize the preceding method for the optimization of linear TD systems; taken together, Eqs. ( 51), ( 49) and (52) form the following QP problem which can be easily solved by standard numerical methods such as the quadprog function in MATLAB: Our task is to determine the matrices Ξ, Λ and b for the problem and then construct the QP solver.

Remark 1
In a linear time-invariant TD system, E(t), F(t), G(t), H(t), V(t) and W(t) are matrices of constants, that is, we simply have where [ ] denotes greatest integer value (see [17]) or we can use some approximations produced by numerical manipulations.

Remark 3
For a system described by ( 26)-(28) with time-invariant weighting matrices in the PI, the entries in Q(t) and R(t) are constant and x (t)Qx(t) + u (t)Ru(t) dt; thus the matrix Ξ in (56) becomes Consider a simple example to illustrate the method.A TD system x is to be controlled to minimize the cost functional Find the optimal control, state and cost.
By constructing the matrices P and L and putting all the findings on (59), ( 60) and (61), we solve the QP problem by calling quadprog in MATLAB.Also, we solve this example for different order of approximations and a comparison of the optimal costs is made in Table 1.The time histories of the optimal control and states for k = 2 and M = 5 are depicted in Fig. 1 and it is clear that the optimal control compares very well with that was presented in [7].However, as we see in Table 1, by changing k = 3, we get a better result and this means that by increasing the number of subintervals, we can obtain more accurate results.In [18, p. 117], the author has used linear Legendre multiwavelets to solve this problem; the graph of optimal control which obtained by using k = 3 has a slight similarity with referenced graph which given by Banks and Burns in [7, p. 194] and the results do not provide a satisfactory approximation where the optimal control values are not compatible with the exact and approximated values (e.g., u * (t = 0); see [18, p. 118] and the references therein).This example shows that this method is much better than the Khellat approach to solve such time-invariant TD problems outlined in [18].For k = 2 and M = 5, we find Source Estimated cost H.T. Banks, J.A. Burns [7] 3.3991 A.Y. Lee [11] 3.4827 F. Khellat [18] 3.43254 N. Haddadi, Y. Ordokhani, M. Razzaghi [19]

Solving Example 1 with terminal inequality constraints Consider the terminal inequality constraints which
Lee in [11] imposed to this system, as To solve this problem, we must rewrite (65) in the standard form: Λ i.e z ≤ b i.e .We write the inequality as

Example 2
Let us illustrate the application of remark 2 with a simple first order system containing delays in state and control variables.Given a TD system and the PI Obtain the optimal control, state and PI.
We choose k = 6 and M = 8; according to remark 2, consider two approximation cases: the first case which based on [17] as case I and the second we propose here as case II I .We solve the problem for the two cases.From the new proposed expansion in (42) for case II, we form Constructing Ξ, Λ and b, we then obtain the solutions by using the interior-point-convex algorithm in MATLAB.
The optimal performance indices computed by the present method for each of the cases and the optimal PI which obtained in [17] are given in Table 2.The time histories of u * (t) and x * (t) for the two approximations are depicted in Fig. 2. We know from experience that the exact optimal state or control (if available) for a TD system is a piecewise-defined function, in which the boundary point(s) (or some of these points) is (are) exactly equal to the delay(s).As we see in the optimal trajectories of Fig. 2, for the first approximation, irrespective of the unexpected jump near the point t = 0.3, one obvious boundary point is t = 0.6875, while for our approximation is t = 21/32 that is closer to the control delay t = 2/3; also in the optimal controls, this point for the first and second cases is located at, in turn, 0.3125 and 11/32, which our result is closer to t = 1/3.These comparisons of two approximate solutions show that the second case is a better approximation and it provides more accurate results.

Example 3
This example was first studied in [8].We are interested in finding the optimal control and state which when applied to a TD system expressed by gives an optimal PI J * described by We approximate x(t) and u(t) by dividing the time interval [0, 1] to four subintervals (k = 3) and by mean of 4th degree polynomials (M = 5).Formula ( 14) gives n dx = 2 and n du = 1, hence according to Eqs. ( 42)-( 43) By using the expansion of t and t 2 , substituting the computed matrices in ( 56),( 57) and ( 58) and calling quadratic programming in MATLAB, we thus obtain the solutions.In Table 3, a comparison is made between the numerical results of J * obtained by this approach and those reported in [8], [20], [21] and [22].The optimal curves are plotted in Fig. 3 and it is obvious that the solutions of this approach represent smooth graphs.The simulation curves derived by this new algorithm which based on Chebyshev polynomials are much better than those which presented in [21] by using Bezier curves; the obtained graphs in [21] are hardly plausible solution curves.

Solving Example 3 with constraints
Physical considerations imply that some constraints should be imposed on the optimal control systems and unconstrained systems are less happened.The proposed method is also applicable to such systems with state, state-control and interior point constraints and final conditions.To illustrate this discussion, we impose a state constraint and a terminal condition on the present system as Case 1.
In case 1, we can write

0
, where f 1 (t) = 8(t − 0.6) 2 .To handle this constraint, we must satisfy it at discrete points t nm on every subinterval where we define: t nm = {(n − 1)M + m}/s.For satisfying the obtained inequality constraints at these points, we have to add two matrices which computed from them to the quadprog function.If the obtained state touches the constraint, we set f 1 (t) − → f 1 (t), where ∈ R >0 is an arbitrary constant (as small as possible) which enforced to prevent any violation.By adding the two matrices to the algorithm which solved the unconstrained problem, we find J * = 1.850110.But, by looking very closely at the obtained results, we see that there is a small violation.As already mentioned, by applying the latter condition, we get J * = 1.850244,where we have no violation, and x * 1 (t) does not touch the respective constraint.In case 2, we have the equality constraint; to handle this, we must add this constraint to (52).Hence we get J * = 2.099227 and x(1) = [0.831905937,0.831905937] .The optimal control and state for the two cases are depicted in Fig. 4.

Example 4
This problem was first studied in [20].The problem is to find the optimal control u * (t) which minimizes the PI subject to We choose k = 7 and M = 5; in a similar manner as explained in Example 3, we construct the required matrices and also we use the inverse time operational matrix of Chebyshev wavelets which described in (17).So, the obtained results are plotted in Fig. 5.The optimal value of the PI by the procedure is given in Table 4 and for purposes of comparison we can find the estimated values of the PI which reported in [20] and [23].As can be seen, it should be obvious that the method is so accurate whereas the derived graphs appear smooth; this fact shows the effectiveness of the proposed method.

Example 5
Consider a linear time-varying system with different delays and terminal times described by (see [24]) This system is to be controlled to minimize the PI where the delay and terminal times are as the following three cases: Case 1.
Since the delay and terminal time are changed, n dx and consequently A and D x will be changed; this affects on the optimal PI J * .In case 1, when we solve the QP problem, then by just changing A and D x according to the new value of x , we are able to solve case 3 of the problem and this is a distinct advantage of the method over some conventional methods.Selecting k = 5 and M = 8, then the values of J * by the method are reported in Table 5 and the optimal states and controls shown in Fig. 6.Because of the nature of the control cost expression u 2 (t), we have to pay higher cost for larger energy consumption.As one can see, the values of delay and final time play a significant role in reducing J * ; these facts lead us immediately to the conclusion that a TD system with a longer terminal time probably have more opportunities to be controlled with lower energy cost by using little energy.

Example 6
This example is adopted from [25].Consider a TD system which contains multiple delays as follows: We want to find the optimal control u * (t) which minimizes subject to Eqs. ( 83)-(84) with t f = 16 and different delays as Case 1.
As mentioned previously, by changing the values of the delays, the optimal value of the PI J * varies considerably regarding the optimal final state x * (t f ) (this directly affects on the optimal value of the terminal PI) and control u * (t).In case 1 we find J * = 0.135229 and in case 2, J * = 0.081048.

Example 7
This problem was studied in [1].A second order linear time-varying multi-delay system expressed by is to be controlled to minimize the cost functional Calculate the minimum value of J.
After rescaling the time interval by setting = t/3, the problem converted to finding u * ( ) and x * ( ) which minimize Using the procedure in case II of Example 2, we set x1 ≈ 21/64, x2 ≈ 17/64, and u ≈ 11/64 to obtain more accurate results; hence we choose k = 7.By selecting M = 5, 8 and using the proposed approach, the QP is solved in MATLAB within, in turn, 2.311 and 4.324 seconds.x * (t) and u * (t) are plotted in Fig. 8.As mentioned in section 4.2, for the optimal state x * 2 , one boundary point seems to be t ≈ 0.5, while we find x * 2 on the interval [0, 0.5] as OPTIMAL CONTROL OF LINEAR TD SYSTEMS BY CHEBYSHEV WAVELETS the piecewise-defined function containing eleven sub-functions.Although increasing the number of subintervals is undesirable, accurate information about the location of boundary points and the values of sub-functions is important so that by using a smooth curve fitting technique, we can replace the piecewise-defined function by an alternative function; for example, x * 2 (t) ∼ = g(t), t ∈ [0, 0.5].In Table 6 a comparison is made, so it is found that the results of applying this method are in good agreement.

Example 8
Consider an example of a typical control problem occurring in the chemical and petroleum industries with transport lag introduced by D.W. Ross [26].This problem also investigated in [27] by using Lyapunov redesign technique for linear TD systems.The system described by where x(t) = [x 1 (t), x 2 (t), x 3 (t), x 4 (t)] and u(t) = [u 1 (t), u 2 (t)] ; the performance criterion is 1 0 0 0 0 10 0 0 0 0 1 0 0 0 0 100 We first determine the value of k.Since the proposed method does not deal with infinite time problems, we consider the finite horizon version of this problem; so, assume t f = 8.We select k = 5 and M = 8.By using Remarks 1 and 3 and implementation of the algorithm in MATLAB, we solve the problem easily.The optimal states and inputs as Fig. 9 are obtained.We can see in the obtained curves that the convergence rate of the presented method is higher than the convergence rate of the method developed in [27].

Conclusion
In this paper a state and control parameterization method is introduced to solve optimal control of linear timevarying multiple delays system which containing inverse time.We have transformed the original TD problem into a new one and obtained the optimal trajectory, control and PI by solving a simple QP problem instead of using of such techniques can lead to major difficulties, such as the generalized Riccati method.The problems with combined state and control constraints can be handled by this method.Furthermore, since the method can incorporate the initial and final conditions and interior point constraints directly into the algorithm, there is no need for the separate operations of applying these conditions to the obtained solutions.In a minimum-energy system, it is interesting to see how changing the parameters like terminal time, lags and control weighted matrix affects on the PI.An efficient method for solving a TD optimal control problem should easily be able to resolve the problem for various cases, such as changing the values of delays, initial vector functions and weighting matrices; the proposed method has fairly this capability.Also, illustrative examples demonstrate that this wavelet-based method compares very well with other methods in the accuracy and order of approximation.It is worth mentioning that the convergence of the method is clearly evident in the obtained curves.The presented QP model of a TD optimal control problem has a very simple structure and implementation of it is easy and convenient.As we have seen, it is important to find accurate information about the solutions in terms of the location of boundary points and polynomial equations of sub-functions to achieve high accuracy.Since the product of the delays and the number of subintervals are not often integers, the accuracy of the method may vary with the choice of k.We have seen that the accuracy may be improved by increasing k.Of course, this enhanced accuracy is usually obtained at a cost-namely, increased computation time; hence future work includes improving the choice of k to match systems with different delays.

Figure 1 .
Figure 1.Optimal states and controls for Example 1.

Figure 3 .
Figure 3. Optimal states and control for Example 3.

Figure 4 .
Figure 4. Optimal states and controls for constrained problem.

Figure 5 .
Figure 5. Optimal states and controls for Example 4.

Figure 6 .
Figure 6.Optimal states and controls for Example 5.

Figure 7 .
Figure 7. Optimal states and controls for Example 6.

Figure 8 .
Figure 8. Optimal states and control for Example 7.

Figure 9 .
Figure 9. Optimal states and controls for Example 8.

Table 2 .
Comparison of the optimal PI (Example 2).

Table 3 .
Comparison of the optimal PI (Example 3).

Table 6 .
Comparison of estimated optimal cost with different methods (Example 7).