A frequency-domain method for solving linear time delay systems with constant coefficients

Mengshi Jin·Wei Chen·Hanwen Song·Jian Xu

Abstract In an active control system,time delay will occur due to processes such as signal acquisition and transmission,calculation,and actuation.Time delay systems are usually described by delay differential equations(DDEs).Since it is hard to obtain an analytical solution to a DDE,numerical solution is of necessity.This paper presents a frequency-domain method that uses a truncated transfer function to solve a class of DDEs.The theoretical transfer function is the sum of infinite items expressed in terms of poles and residues.The basic idea is to select the dominant poles and residues to truncate the transfer function,thus ensuring the validity of the solution while improving the efficiency of calculation.Meanwhile,the guideline of selecting these poles and residues is provided.Numerical simulations of both stable and unstable delayed systems are given to verify the proposed method,and the results are presented and analysed in detail.

Keywords Time delay system·Delay differential equation·Approximate solution·Transfer function

1 Introduction

Time delays exist in a wide range of engineering systems.Since time delay may improve the performance of the system,there are various applications of time delay in active vibration control and noise control[1–5].On the other hand,time delay may influence the stability and the dynamic behavior of the system[6–9],which deserves attention.Delay differential equations(DDEs)are employed to describe time-delay systems.They are a special class of differential equations,called functional differential equations.Thus,the approach to solving DDEs differs from that to solving ordinary differential equations.The main difficulty lies in the transcendental character of DDEs[10].

Extensive investigations have been carried out concerning the solution to DDE[10–16].Delay problems always bring about an infinite spectrum of frequencies.An approximation method of this type is the well-known Padé approximation,which results in a shortened repeating fraction for the approximation of the characteristic equation of the delay τ[11,12].The work of Lam[11]gave an algorithm for forming a balanced realization of the all-pass Padé approximants of e?τs.The algorithm generated the coefficients of the Padé denominators based on a scaled recurrence formula.The problem of finding reduced-order models for a class of delay systems was considered in Ref.[12].The delay element e?τswas suggested to be replaced by Padé approximants.

Another approach to solving DDE is to find the response through analyzing the entire delay spectrum instead of approximating the time delay.Several studies have been accomplished to find a solution to DDE by solving its characteristic equation under different conditions[13–15].A related study on analytic solution to linear DDEs can be found in Ref.[13].In their work,Wright[14]studied a Fourier-like analysis of the existence of the solution and its properties for the nonlinear DDEs.Similar approaches to linear DDEs were also reported by Bellman and Cooke[15].The uniqueness of the solution and its properties for the linear DDEs with varying coefficients were studied by W right[14].

Although different methods have been proposed to solve DDEs,they still have a few shortcomings.The solution to DDE can be investigated analytically by using the classical method of steps[16].However,the calculations will quickly become inconvenient without revealing essential properties of the solution.The readers can refer to Ref.[17]for the example of using the method of steps to solve DDEs.Since a few frequencies are actually dominant in the appearance of oscillations[10],it can be envisioned that an approximation in the frequency domain may be feasible.Therefore,it is of necessity to present a frequency method that guarantees the simplicity of calculation without sacrificing accuracy,which is the motivation of this paper.Jin et al.[18]studied the identification of time delay dynamic system in experiment,while this work can help to better understand the properties of delayed system.

This paper concerns solving a class of linear DDEs with constant coefficients.Firstly,the solution in the frequency domain is deduced and expressed by the transfer function,system parameters and equivalent initial conditions.The theoretical transfer function is represented by infinite number of poles and residues,which can be acquired by fitting the curve of the theoretical frequency response function(FRF).While the transfer function is the sum of infinite terms,it can be truncated by selecting the dominate poles and residues.Moreover,the guidance of selecting these dominant poles and residues is provided.

It is worth noting that the method presented by Qin[6]may lead to the opposite judgement of the stability of the delayed system when the term with time delay was expanded by Maclaurin series.To be more specific,after the delayed term was expanded,one may get the conclusion that the system was unstable while the system was actually stable.However,the method proposed in this paper will not result in the wrong judge of the stability.In other words,the truncation of the transfer function does not affect the stability of the original delayed system so that one will not make mistakes when judging the stability of the system.In addition,since the method is developed in the frequency domain,there will be fewer cumulative errors appearing in other iteration methods,which contributes to the accuracy of calculation.Another feature of the method is that usually a few truncation terms are sufficient for getting a reasonably accurate solution,leading to high efficiency.

The rest of this paper is organized as follows.Section 2 describes the category of the DDEs concerned in this paper.The frequency method to solve the DDE is proposed in Sect.3.Numerical examples of both stable and unstable systems are presented to illustrate the method in Sect.4.Conclusions are provided in Sect.5.

2 Problem formulation

The equation of motion of a time delay system with n degrees of-freedom(DOF)is given as

with initial conditions

where M,C,K ∈ Rn×nare the mass matrix,damping matrix,stiffness matrix,respectively.u(t)∈Rnis the vector of external excitation.G ∈ Rn×nis the feedback gain matrix,and τ is the time-delay.A ll the system parameters and feedback parameters are real constants.x(t)∈Rnis the vector of the displacement response of the time delay system.

The characteristic equation of Eq.(1)is

where s represents the pole of the system.

In this paper,we discuss the systems where poles are complex conjugate pairs as

The solution to Eq.(1)consists of two parts.The first part is the response of the system over the time interval of t ∈ [0,τ),while the second part is the response over the time interval of t∈ [τ,∞).In this paper,the system over the time interval of t∈ [0,τ)and t∈ [τ,∞)w ill be termed as open-loop system and closed-loop system,respectively.

The derivation of the proposed method can be developed for a general class of inputs.In this research we specifically focus on transient response to highlight the essential aspect of the method proposed.Response under external excitation can be derived in a similar manner.Therefore,this paper presents the derivation of the following DDE

3 Solution to time delay system

The equation of motion of the open-loop system is given as

with initial conditions

where xop(t)∈Rnis the vector of the displacement response of the open-loop system.The subscript?oprefers to “open loop”.

Since it is straightforward to obtain the solution to Eq.(6),the process is omitted here[19].The displacement and velocity responses can be acquired and denoted as xop(t)and,respectively.Substituting t=τ into xop(t)andyields the responses at the instant t=τ,which are given as

3.1 Solution to closed-loop system

This section presents the process of solving for the solution to the closed-loop system in detail.The equation of motion of the closed-loop system is a DDE,which can be expressed as

where xcl(t)∈Rnis the vector of the displacement response of the closed-loop system.The subscript?clrefers to“closed loop”.

The equivalent initial conditions are denoted as

It is worth noting that xcl0andare two unknown vectors and will be determined with continuity conditions later.Equivalent initial conditions are the responses of the closed loop system at the instant t=0.However,the domain of Eq.(9)is t ∈ [τ,∞),which means there is no response of the closed-loop system over the time interval t∈ [0,τ).Therefore,the equivalent initial conditions have no physical meaning.This is the reason why they are termed as equivalent initial conditions.

The solution to the closed-loop system can be expressed as

where L?1(·)is inverse Laplace transform and

is the transfer function of the closed-loop system.

3.2 Equivalent initial condition

As stated in the Sect.2,the solution to the time delay system is calculated over two time intervals,which means the responses of the open-loop and closed-loop systems are respectively acquired.However,since the displacement and velocity of the two systems are continuous at the instant t=τ,we have the following continuity conditions

The left-hand side of Eq.(13)is the responses of the closed-loop system containing the unknown equivalent initial conditions xcl0and,while the right-hand side is the responses of the open-loop system and have already been obtained in Eq.(8).

Once Eq.(13)is solved,complete solution to the closed loop system can be completely obtained by substituting the specific values of the equivalent initial conditions into Eq.(11).

3.3 Truncation of transfer function

The transfer function given in Eq.(12)can be written in the form of a matrix as

Substituting Eq.(14)into Eq.(11)yields the solution to the closed-loop system,given by

where mjand cqjare the elements of the mass matrix and damping matrix,respectively.

The element of the transfer function matrix can be expressed as the sum of terms involving the poles and their residues as

where srandare the pole and its complex conjugate pair.rRijandrRij?are the residue and its complex conjugate pair.

Thus,Eq.(16)can be further rewritten in detail as

It can be seen from Eq.(17)that the transfer function is the sum of infinite terms.From the view of feasibility,the transfer function should be truncated and written in the form of finite items as

where N is the truncation order.

The truncated transfer functionin Eq.(19)is the sum of finite terms expressed in terms of the poles and residues.Theoretically,there are numerous choices we can make to select the N pairs of poles to truncate the transfer function.Therefore,it is necessary to provide a rule of selecting the poles.

As can be seen from the expressions of the solution(Eqs.(15)and(18)),the real parts of the poles are the powers of the exponents.This means the larger a real part is,the greater influence its corresponding pole will ultimately have on the solution.In other words,the poles with larger real parts dominate the solution.Consequently,the rule of selecting poles is based on the values of their real parts.

Specifically,if the system is stable,we first sort the poles according to the values of their real parts(see Eq.(4)).The sequence rule is shown in the following inequality as

It is easy to understand that the approximate solution with higher truncation order will have higher accuracy.Nevertheless,in practical problems usually a small number of terms would be sufficient.

Considering an unstable system with U unstable poles,we still sort the poles according to the magnitudes of their real parts.The sequence rule is presented as

It can be seen that even if there is only one pair of poles chosen in the truncated transfer function,the stability of the original system will not be changed since the selected poles have positive real parts.Therefore,this method of truncation will not lead to the wrong judgment of the stability of the delayed system.

When the truncation order N is determined,the poles that satisfy the inequality(20)or(21)and their corresponding residues will be chosen to truncate the transfer function.Then,substituting the poles and the residues into Eqs.(15)and(18),the approximate solution to the closed-loop system can be obtained with the infinite signs∞r(nóng)eplaced by the truncation order N,given by

where

The velocity can be obtained by differentiating the displacement.Then,the equivalent initial conditions can be determined by continuity conditions given in Eq.(13).Lastly,substituting the specific values of the equivalent initial conditions into Eq.(22)yields the approximation solution to the closed-loop system.

So far,the solution to the original time delay system has been completely obtained since the solution to the open loop system and the one to the closed-loop system have been acquired.

4 Numerical simulation

Numerical solutions to DDEs can be obtained by XPPAUT[20],a common software used in the calculation of differential equations.Since there is no analytical solution to DDE,the solutions obtained by XPPAUT will be regarded as“accurate solutions”,which the approximate solutions will be compared with.

In this section,numerical simulations of both stable and unstable time delay systems are given to show the efficiency and accuracy of the proposed method.The equations of motion of both stable and unstable time delay systems are the same as Eq.(1).

4.1 Stable time delay system

In this case,the system parameters are

The feedback parameters are set as

The initial conditions are given as

According to Eq.(12),the theoretical FRF matrix of the stable delayed system and that of the corresponding uncontrolled system can be obtained.Take one FRF element,say H1,1,as an example.The two FRF curves of the stable delayed system and the uncontrolled system are shown in Fig.1.

After using the vector fitting method[21]to fit the FRF curve of the delayed system,we can obtain as many poles as we want.The first 30 poles(15 pairs)are shown on the complex plane in Fig.2.

It can be seen from Fig.2 that the real parts of the poles are all negative since the delayed system is stable.Also,the poles are sorted by their imaginary parts.However,according to sequence(20),the poles should be sorted by their real parts.Table 1 shows the first 10 pairs of the poles sorted according to sequence(20).The circles with squares in Fig.2 represented the 10 pairs of poles listed in Table 1.

Fig.1 Curves of FRF H1,1 of the stable delayed and uncontrolled system

Fig.2 Distribution of the poles of the stable delayed system

Table1 The sorted poles according to real parts

It can be seen from Table 1 that the poles at the beginning of the sequence have real parts with larger values.As is shown in Eq.(23),the real parts of the poles are the powers of the exponents.Thus,the first three frequencies dominate in the approximate solution will be around 10.5,3.8,and 9.0Hz.

Fig.3 Accurate and approximate solutions of DOF 1(N=1)

4.1.1 Free decay

We consider the free decay of the stable time delay system,which means there is no external excitation.Choosing different truncation orders,we can get several groups of approximate solutions by the mentioned method.Three groups of results with truncation order N being 1,3,and 10 are picked out and analysed to show the validity of the method.

With the truncation order N=1,Fig.3 shows the comparison between the accurate solution and the approximate solution of DOF 1.

As can be seen from the enveloping lines,the accurate solution(blue line)includes multiple frequencies,while the approximate solution(red line)includes only one frequency,which is especially obvious in the zoom-in view of the responses from 0.6 to 1.2s.During the last second,the tendency the red curve,which represents approximate solution coincides well with the blue one standing for accurate solution.However,the amplitudes of the red curve are larger than those of the blue one.Similar phenomena can be observed for DOF 2 and DOF 3.Therefore,the figures of comparisons are omitted here.

Now the truncation order is chosen as N=3 and Fig.4 shows the comparison between the accurate solution and the approximate solution of DOF 1.

According to Fig.4,the tendency of the approximate solution(red line)is consistent with that of the accurate solution(blue line).The zooms of the beginning and the last second show that the curve of the approximate solution coincides quite well with that of the accurate solution.The consistency of the two curves indicates approximation solution is capable of describing the dominant components in the accurate solution.Combining Figs.3and4,it can be obviously noticed that the approximate solution with N=3 is far more satisfactory than that with N=1.

Fig.4 Accurate and approximate solutions of DOF 1(N=3)

Now the truncation order is selected as N=10.The approximate solutions of DOF 1,2,and 3 are obtained and presented in the following.Besides,the fast Fourier transforms(FFTs)of the solutions are conducted to indicate the frequency components,as is shown in Figs.5,6,and 7.

On one hand,Figs.5a,6a,and 7a indicate that the approximation solutions of the three degrees of freedom are highly close to the accurate solutions.On the other hand,Figs.5b,6b,and 7b show the curves of the FFTs of the approximate solutions are in good agreement with those of the accurate ones.Take Fig.7b as an example.The three peaks in Fig.7b have the frequencies of 3.784,9.094,and 10.437Hz,respectively,which are consistent with the first three frequencies shown in Table 1.

Choosing truncation orders ranging from 1 to 10 and repeating the above procedures,we can obtain the corresponding approximate solutions.Here,the figures of the results are omitted due to the limit of space.Figure 8 shows the errors of the FFTs of these approximate solutions with different truncation orders(Figs.9,10).

It is obvious that with the increase of truncation order,the accuracy of the approximate solution is improved.Also,the error reaches a plateau when the truncation order increases to 3,which means that approximate results with the truncation orders higher than 3 are quite similar to those with N=3.The convergence of the curves reveals that the truncation order does not have to be very high to get an approximate solution with satisfactory accuracy.Based on this,the proposed method is proved to have convenience and efficiency in calculation without the sacrifice of accuracy.

4.1.2 Forced response

Now we consider the forced response of the stable time delay system.A sinusoidal excitation u1(t)=200sin(6πt)is applied on DOF 1 and no external excitation is applied on the other two degrees of freedom.Therefore,the frequency of the only external excitation is 3 Hz.

Fig.5 Accurate and approximate solutions of DOF 1(N=10).a Solutions to the DDE.b FFTs of the solutions

Fig.6 Accurate and approximate solutions of DOF 2(N=10).a Solutions to the DDE.b FFTs of the solutions

Fig.7 Accurate and approximate solutions of DOF 3(N=10).a Solutions to the DDE.b FFTs of the solutions

Fig.8 Errors of the FFTs with different truncation orders

The process of obtaining approximate solution of the forced response is similar to that of the free response.The following shows the approximation results with the truncation order N=10.

It can be seen from the above figures that the approximation solutions of the three degrees of freedom are highly close to the accurate solutions both in time domain and in frequency domain.The three peaks in Fig.11b with the frequencies of 3.792,9.071,and 10.445Hz are consistent with the first three frequencies shown in Table 1,respectively.The peak with the frequency of 3.0Hz is caused by the only external sinusoidal excitation.

Choosing truncation orders ranging from 1 to 10 and repeating the above procedures,we can obtain the approximate solutions to the DDE.Figure 12 shows the errors of the FFTs of these approximate solutions with different truncation orders.

It is obvious that with the increase of the truncation order,the accuracy of the approximate solution is improved.The convergence of the curves shown in Fig.12 reveals the truncation order does not have to be very high to get a quite satisfactory approximate solution.Therefore,the proposed method makes the calculation efficient without sacrificing accuracy.

4.2 Unstable time-delay system

In this case,we study an unstable time-delay system whose system parameters and initial conditions are the same as those of the stable delayed system(see Sect.4.1).The differences are the feedback parameters,which are set as

Since the system is unstable,no external force is applied.

According to Eq.(12),the theoretical FRF matrix of the unstable delayed system and that of the corresponding uncontrolled system can be obtained.Take one FRF element,say H1,1,as an example.The two FRF curves are shown in Fig.13.

After using the vector fitting method to fit the FRF curve of the delayed system,we can obtain as many poles as we want.The first 30 poles are shown on the complex plane in Fig.14.

It can be seen from Fig.14 that four pairs of poles(marked by crosses)have positive real parts which means the delayed system is unstable.According to sequence(21),the poles should be sorted by their real parts.Table 2 shows the first 10 pairs of the poles sorted according to sequence(21).The poles with squares in Fig.14 represented the10 pairs of poles listed in Table 2.

Fig.9 Accurate and approximate solutions of DOF 1(N=10).a Solutions to the DDE.b FFTs of the solutions

Fig.10 Accurate and approximate solutions of DOF 2(N=10).a Solutions to the DDE.b FFTs of the solutions

Fig.11 Accurate and approximate solutions of DOF 3(N=10).a Solutions to the DDE.b FFTs of the solutions

Fig.12 Errors of the FFTs with different truncation orders

Fig.13 The curves of FRF H1,1 of the unstable delayed and uncontrolled system

Fig.14 The distribution of the poles of the unstable delayed system

Table2 Sorted poles according to real parts

Table 2 shows that the poles in the front of the sequence have the real parts with larger values.After choosing the finite poles with larger real parts,the frequencies of these poles will dominant in the appearance of oscillations.

The process of obtaining approximate solution of the unstable delayed system similar to that of the stable delayed system.Using the proposed method and choosing the truncation order as N=15,we can get the approximate solutions,shown in Fig.15.

Since this delayed system is unstable and will never be put into practice,it is unnecessary to compare the specific values of the solutions.Figure 15 illustrates that the proposed method will not mistakenly treat an unstable delayed system as a stable one.This means the stability of the delayed system will not be changed due to the truncation of the transfer function,enabling the method to ensure the correctness of stability.

Fig.15 Comparisons between accurate and approximate solutions with N=15.a Solutions of DOF 1.b Solutions of DOF 2.c Solutions of DOF 3

5 Conclusions

The aim of this paper is to provide a frequency method to obtain the solutions of a class of linear DDEs with constant coefficients.

Firstly,the equation of motion is given,which clearly describes the form of the DDE concerned in this paper.The solution to the DDE consists of two parts.One is the solution over t ∈ [0,τ),which is called an open-loop system.The other is the solution to the system over t ∈ [τ,∞),called a closed-loop system.

Secondly,the solution to the closed-loop system is expressed in terms of the transfer function.Since the theoretical transfer function consists of infinite poles,it can be truncated,which is the basic idea of this method.While there are infinite poles,only a few of them dominate the solution.Therefore,the guideline of selecting poles is provided according to their real parts.Poles with larger real parts and their corresponding residues are selected to truncate the transfer function.

Finally,numerical simulations are given where the approximate solutions obtained by the presented method are compared with the solutions obtained by a well-known software.As is shown in the results,only a few truncation items are needed to get a satisfactory approximate solution,which makes the method quite efficient.

AcknowledgementsThis work was supported by the National Natural Science Foundation of China(11272235).