，，2＊

1.Army Engineering University of PLA，Nanjing 211101，P.R.China；

2.State Key Laboratory of Mechanics and Control of Mechanical Structures，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，P.R.China；

0 Introduction

Time delay exists commonly in control applications，such as in digital controller where the time delay is resulted from using sampling and zero－order holder，and in human－interaction systems where the time delay is produced due to the delay of human′s response.Dynamical systems with time delays are called time－delay systems，which can be classified into two categories：Retarded type and neutral type.On one hand，the presence of time delay may deteriorate the system′s performance and even destabilizes the system，which may occur even when the delay is very short［1］.On the other hand，the effect of time delay on the system dynamics may be positive，such as in the study of sway reduction of cranes［2］.Thus，stability analysis of time－delay systems has been one of the major concerns in many control applications.

Usually，stability is in the sense of Lyapunov′s asymptotical stability.An equilibrium of a timedelay system of retarded type is asymptotically stable if all the characteristic roots of the corresponding linearized system have negative real parts only.This is true for time－delay systems of neutral type under certain conditions.Many methods and criteria have been established for the stability analysis of time－delay systems，such as the D－subdivision method［3］，the method of stability switches［4］，and the stability criteria developed on the basis of Argument Principle，including the Nyquist criterion［5］，Integral estimation criterion［6］，Stepan－Hassard theorem［7］，etc.The D－subdivision method works effectively in determining the stable regions in a parametric plane such as the feedback gain plane of a controlled system by using the critical stable conditions.The method of stability switches is preferable when only one parameter is focused.The stability criteria are used mainly for stability test of given timedelay systems，but they can also be used for stud－ying the problems of stability switches［8］.In some cases，the solutions starting from points close to a stable equilibrium may have a very long transient.

Abscissa is a real number defined as the real part of the rightmost characteristic roots of a dynamical system，and it is an index for measuring stability.An equilibrium of a time－delay system is asymptotically stable in the Lyapunov′s sense if the abscissa is negative，the smaller（the larger in absolute）the abscissa is，the better the stability is，and the solutions starting from points close to a stable equilibrium will decay to the stable equilibrium faster.Thus，in the control design of a delayed feedback control，it is important not only to determine the stable regions in the gain plane or gain space，but also to find out the pair of optimal feedback gains that minimizes the abscissa within a given stable region.Roughly speaking，all the above－mentioned stability criteria can be used for the calculation of the abscissa of a given time－delay system，among them the integral estimation criterion seems more effective in implementation.These criteria could be used directly for the calculation of the abscissa one－by－one at each gridding node，but seemingly it is not easy to obtain the optimal abscissa and the optimal feedback gains because manual intervention seems necessary in the calculation.Thus，more effective method or algorithms are needed to find the optimal abscissa and the optimal feedback gains within a stable region in the gain plane.This paper presents a simple algorithm for labeling the stable region in feedback gain plane with different abscissa on the basis of the D－subdivision method，and it is found that the sub－region labeled with the smallest abscissa is preferable in applications.For clarity in presentation，the main results are introduced for the controlled pendulum or inverted pendulum with delayed feedback.The proposed algorithm works also for labeling the stable regions of a pair of feedback gains of any controlled systems with delayed feedback.

1 Labeling of Stable Regions

Pendulum and inverted pendulum are two very popular models for many mechanical systems/structures.When a delayed acceleration－velocity－position control is used for controlling a pendulum or inverted pendulum，the controlled system is described by

wherem，c，kare the system parameters，τ1≥0，τ2≥0，τ3≥0the time delays，andka，kd，kpthe feedback gains.Eq.（1）withk＞0corresponds to apendulum and withk＜0corresponds to an inverted pendulum.Eq.（1）is called retarded （or neutral）type whenka＝0（orka≠0）.The characteristic equation is

wherep0（λ）＝mλ2＋cλ＋k.The abscissa is defined by

whereR（λ）represents the real part of complex numberλ .D（λ）has infinite many roots，but the number of roots with positive real part must be finite.Thus the abscissaα must be finite［3］.For the case of｜ka｜＜1，the time－delay system is asymptotically stable ifα＜0.Without loss of generality，assume thatα＜0whenτ1＝0，τ2＝0，τ3＝0.Then conditions onka，kd，kpcan be obtained by using the Routh－Hurwitz criterion.Due to the continuous－dependence on the delays，the timedelay system keeps asymptotically stable if the delays are small enough.

1.1 Introduction of D－subdivision method

For given small delays，the stable regionswith respect to the gainska，kd，kpcan be obtained by using the D－subdivision method，where the boundaries of the stable regions are plotted by using the critical stable curves determined fromD（iω）＝0（i2＝－1）.For a delayed proportionalderivative （PD）feedback with equal delays，ka＝0andτ2＝τ3＝τ，for example，the boundaries of the stable regions are determined by R（D（iω））＝0and I（D（iω））＝0（where I（z）stands for the imaginary part of complex numberz），which give

whenω≠0，and forω＝0one has

The critical curves divide the gain plane（kp，kd）into many open regions，in which none，one or more are stable.For finding the stable region and labeling the stable region with abscissa，the σ－critical stable curves，determined byD（σ＋iω）＝0or equivalently R（D（σ＋iω））＝0，I（D（σ＋iω））＝0，will be used，they are plotted by｛（kp，kd）：0≤ω＜＋∞｝，where

Similar results can be obtained for a delayed PD feedback withka＝0andτ2＝τ，τ3＝2τ，and for a delayed acceleration－derivative （AD）feedback withkp＝0andτ2＝τ3＝τ，as well as the one withkp＝0andτ2＝τ，τ3＝2τ.

For any point passed by aσ－critical stable curve，the corresponding time－delay system has a characteristic root with real partσ.Theσ－critical stable curves divide the gain plane into a number of sub－regions，which can be classified into two classes：theσ－stable ones for which the correspondingD（λ）has roots with real parts less than σonly，and theσ－unstable ones for which the correspondingD（λ）has at least one root with real parts larger thanσ.The 0－stable ones are the stable regions in the Lyapunov′s sense.The stable regions as well as theσ－stable regions can be determined graphically.

Theσ－critical stable curves with different values ofσmay intersect with each other，or do not intersect with each other at all.For the case when intersection happens，the intersect point should be marked with the color of theσ－critical stable curve corresponding to the largestσ.Firstly，choose two real figuresσmin＜0，σmax＞0such that the abscissaαfor all parameter combinations in a given region of the parameter plane is in the interval［σmin，σmax］，then the process of labeling the stable region can be completed in the following major steps.

1.2 Subdivision of given region viaσ－critical stable curves with negativeσ

Starting fromσ＝σmintoσ＝0by a small step δσ，and for each node ofσ，theσ－critical stable curves are drawn，and each point on theσ－critical stable curve is marked with a designated color in the color set.Letσ2＞σ1，if theσ1－critical stable curve intersects with theσ2－critical stable curve，the intersect points should be marked with the color of theσ2－critical stable curve.Ifσminis cho－sen small enough（or equivalently，is large enough），every point in the given region of the parameter plane is labeled by the color of theσcritical stable curves withσ ∈ ［σmin，0］.

1.3 Erasure of unstable regions viaσ－critical stable curves with positiveσ

Points passed by aσ－critical stable curve with positiveσcorrespond to the case when the timedelay system has at least one pair of characteristic roots with positive real part，so they are not belong to the stable regions.Further effort is required to erase the unstable points from the given region of the parameter plane labeled by the color of theσ－critical stable curves withσ ∈ ［σmin，0］，by using white color of theσ－critical stable curvesForσfrom 0toσmax（or from σmaxto 0）by a small stepδσ，and for each node σ，plot theσ－critical stable curve by white color.Then，all the points in the given region of the parameter plane have been marked with different color characterizing the stable region，showing different level of the abscissa in［σmin，0）.

1.4 Asymptote issue

Theσ－critical stable curves can be either continuous or discontinuous.For the continuous case，the stable region can be labeled simply by using the above two steps.If there are some break points on theσ－critical stable curves，just like in the applications to time－delay systems of neutral type，the plot in these points will create some asymptotes.For generality，assume that theσ－critical stable curves are plotted by （f（σ，ω），g（σ，ω））asωvaries from 0to＋∞，and they have discontinuity.The first asymptote is defined by

wherexandyare the gain values of the feedback control.

Except the first asymptote，stability switches do not occur at the both sides of the asymptotes，which are not a part of critical stable curves，so we have to avoid these asymptotes when plotting theσ－critical stable curves.For a fixedσ，let 0＜ωc1，ωc2，…，ωck，… ，be the roots of the denominator off（σ，ω）andg（σ，ω），then all the points（f（σ，ωci），g（σ，ωci））withi＝1，2，…should be avoided in plotting theσ－critical stable curves with thisσ.

2 Algorithm

Below is the algorithm for labeling the stable region in feedback gain plane of a time－delay system with abscissa.N，Kare two integers satisfyingσmin＝N＊δσandσmax＝K＊δσ.

INPUT：The color setC；σ－critical stable equations：x＝f（σ，ω），y＝g（σ，ω）；the range of considered regionS；frequency lengthM；frequency step lengthδω ；negative integerN；positive integerK；step lengthδσ.

OUTPUT：The labeled stable regionS.

Step 1Forn＝N，N＋1，…，Kdo Steps 2—3.

Step 2σ←n＊δσ.Ifσ＜0，choose the color identified byσinC，else choose white color.

Step 3m＝0，1，…，M，do Steps 4—5.

Step 4ω←m＊δω.

Step 5Plot point （f（σ，ω），g（σ，ω））with the designated color.

Step 6Output the plot ofS.

Terminate.

If the denominators off（σ，ω）andg（σ，ω）have nonzero real rootωσforωwith a fixedσ，Step 4will be replaced with

Step 4ω←m＊δω.to Step 3.

When plotting theσ－critical stable curves with the available mathematical software，we draw lines rather than isolated points.To avoid plotting the asymptotes，we calculate the break points off（σ，ω）andg（σ，ω）with the index ofω firstly，and then draw the piecewise curves.A successful application of the proposed algorithm requires a suitable estimation ofσminandσmax.In many applications，the delays are small，thus the estimated values ofσminandσmaxcan be chosen based on the abscissa when all the delays equal zero.When the delays are not small，the estimation ofσminandσmaxis left for further investigation.

3 Examples

In the following two case studies，only the simply connected stable region close to the origin of the gain plane is considered.

3.1 Example 1

Consider the following controlled system in dimensionless form

wherepandqare the gain values of position and velocity respectively，andris the ratio coefficient of position delay and velocity delay.The delay values are assumed small，and only the stable region containing the origin of the parameter plane is considered.The characteristic function of system （Eq.（3））is

Separating the real and imaginary parts ofD（σ＋iω）＝0，and solving the gainsp，dfrom linear equations R（D（σ＋iω））＝0，I（D（σ＋iω））＝0gives

Fig.1shows the stable regions close the origin in the gain plane of system （Eq.（3））withr＝1and different parameter combinations，labeled with abscissa within ［－6，0］by using the proposed algorithm.The algorithm can be implemented with Matlab.Fig.2 presents the labeled stable regions of system （Eq.（3））withr＝2，labeled with abscissa within ［－8，0］.Both cases show that the increase ofτnot only shrinks the stable region but also decreases the abscissa，and on the contrary，the increase ofξnot only enlarges the stable region but also increases the abscissa.In addition，F(xiàn)igs.1—2show that the stable regions of system （Eq.（3））withr＝2are much larger than the corresponding ones withr＝1.This means that from the viewpoint of stable region，a delayed PD feedback with the delay in position feedback double of that in velocity is preferable in applications.

Fig.1 Labeled stable regions of system （Eq.（3））with r＝1，σmin＝－6，σmax＝15in［－4，80］×［－3，17］

Fig.2 Labeled stable regions of system （Eq.（3））with r＝2，σmin＝－8，σmax＝15in［－8，300］×［－4，40］

3.2 Example 2

Consider the following delayed system in dimensionless form

whereaandpare the gain value of acceleration and velocity respectively，＜1，andris the ra－tio coefficient of acceleration delay and velocity delay.The corresponding characteristic function is

D（λ）＝（1＋ae－rλτ）λ2＋ （2ξ＋de－λτ）λ＋1

Fig.3shows the stable regions close the ori－gin in the gain plane of system （Eq.（4））withr＝1and different parameter combinations，labeled with abscissa within［－19，0］by using the proposed algorithm.Fig.4presents the labeled stable regions of system （Eq.（4））withr＝2，labeled with abscissa within［－13，0］.Again，both cases show that the increase ofτnot only shrinks the stable region but also decreases the abscissa，and on the contrary，the increase ofξnot only enlarges the stable region but also increases the abscissa.In addition，F(xiàn)igs.3—4show that the stable regions of system （Eq.（4））withr＝2are much larger than the corresponding ones withr＝1.This means that from the viewpoint of stable region，a delayed AD feedback with the delay in acceleration feedback double of that in velocity is preferable in applications，which is agreement with the result given in Ref.［9］.

Fig.3 Labeled stable regions of system （Eq.（4））with r＝1，σmin＝－19，σmax＝15in［－1，1］×［－4，20］

Fig.4 Labeled stable regions of system （Eq.（4））with r＝2，σmin＝－13，σmax＝15in［－1，1］×［－4，40］

4 Conclusions

An algorithm based on the D－subdivision method is proposed for labeling the stable region in the plane of feedback gains of dynamical systems under a delayed feedback control with different abscissa.Two main steps in the labeling process are required，one is subdivision of the stable region，and the other is erasure of the unstable regions.The labeling simply uses a color in a designated color set to plot theσ－critical stable curves，and can be easily implemented with computer codes.A successful application of the proposed algorithm requires a suitable estimation of the abscissa.The two case studies show that for the controlled pendulum with a delayed feedback，the stable region can be substantially enlarged if the delays are properly chosen.The algorithm works for labeling the stable regions of a pair of feedback gains for any controlled systems with a delayed feedback.

Acknowledgements

This work was supported by the National Natural Science Foundation of China （No.11372354）.The authors thank Dr.Zhang Li of Nanjing University of Aeronautics and Astronautics for bringing a labeling algorithm based on DDE－BIFTOOL into their attention.

［1］ HUANG R，HU H Y，ZHAO Y H.Designing ac－tive flutter suppression for high－dimensional aeroelastic systems involving a control delay［J］.Journal of Fluids and Structures，2012，34（10）：33－50.

［2］ NAYFEH A，MASOUD Z，NAYFEH N.A smart Sway controller for cranes from theory to laboratory to industry［C］//Vibration Problems ICoVP 2011 Supplement the 10th International Conference on Vibration Problems.Liberec，Czech：Technical University of Liberec，2011：14－29.

［3］ QIN Y X，LIU Y Q，WANG L，ZHENG Z X.Stability of dynamic systems with time delays［M］.Beijing：Science Press，1989.（in Chinese）

［4］ KUANG Y.Delay differential equations with applications to population dynamics［M］.New York：Academic Press，1993.

［5］ FU M Y，OLBROT A W，POLIS M P.Robust stability for time－delay systems：The edge theorem and graphical tests［J］.IEEE Transactions on Automatic Control，1989，34（8）：813－820.

［6］ XU Q，STEPAN G，WANG Z.Delay－dependent stability analysis by using delay－independent integral evaluation［J］.Automatica，2016，70（8）：153－157.

［7］ STEPAN G.Retarded dynamical systems：stability and characteristic functions［M］.Essex：Longman Scientific，1989.

［8］ XU Q，STEPAN G，WANG Z.Numerical stability test of linear time－delay systems of neutral type［C］//Time Delay Systems，Advances in Delays and Dynamics 7.Berlin Heidelberg：Springer，2017：77－91.

［9］ WANG Z，HU H，XU Q，et al.Effect of delay combinations on stability and Hopf bifurcation of an oscillator with acceleration－derivative feedback［J］.International Journal of Non－Linear Mechanics，2017，94（9）：392－399.