I have a linearized time-delay system as follows:
$$\frac{\mathrm d X}{\mathrm d t} = a[X(t)-X^*] + b [X(t-R) - X^*], $$ where $a$, $b$ are constants, $R$ is the constant delay, and $X^*$ is the equilibrium point.
How do I konw if, and under what condtion, this system is stable?
On
When $R=0$ or $b=0$ you can simply just look at whether $a+b<0$, but when $R\neq0$ and $b\neq0$ the stability analysis requires some more effort.
You can formulate the stability analysis by assuming a solution of the following form
$$ X(t) = C\,e^{\lambda\,t} + X^*. $$
Plugging this into the left and right hand sides of the differential equations gives
\begin{align} \frac{\mathrm d X}{\mathrm d t} &= C\,\lambda\,e^{\lambda\,t}, \\ a[X(t)-X^*] + b [X(t-R) - X^*] &= a\,C\,e^{\lambda\,t} + b\,C\,e^{\lambda\,(t - R)}. \end{align}
Equating these to each other and factoring out common terms gives
$$ \left(\lambda - a - b\,e^{-\lambda\,R}\right)\,C\,e^{\lambda\,t} = 0. $$
The trivial solution $C\,e^{\lambda\,t} = 0$ of the above equation would imply $X(t) = X^*$, so would be stable. The other possible solutions require solving
$$ \lambda = a + b\,e^{-\lambda\,R}, $$
so the whole system would be stable if all of the solutions for $\lambda$ have a negative real part. By using the substitution $\lambda = \tfrac{\mu}{R} + a$ the above equation can also be written as
$$ \mu\,e^\mu = b\,R\,e^{- a\,R} = x, $$
which has infinitely many solutions and can be found by using different branches of the Lambert W-function.
The stability condition that all of the solutions for $\lambda$ must have a negative real part can also be expressed as that all solutions for $\mu$ must have a real part larger than the real part of $-a\,R$ when $R<0$ or smaller than the real part of $-a\,R$ when $R>0$. If I remember correctly the principle branch of the Lambert W-function has the largest real part and there is no lower bound for the real part when considering all branches, but I couldn't find a source for this. So when $R>0$ the system is stable if the real part of the value of principal branch of the Lambert W-function evaluated at $x$ is smaller than the real part of $-a\,R$ and when $R<0$ it will always be unstable because there will always be a branch which will violate the previously stated condition.
For example when assuming that $R>0$ and using a couple of values for $a\,R$ and $b\,R$ gives the following image, where blue means stable and yellow means unstable:
In the above image the real part of the value of principal branch of the Lambert W-function evaluated at $x$ plus $a\,R$ is plotted. This value thus implies stability when it is negative.
The edge of stability can be found by setting the real part of $\mu$ equal to $-a\,R$, so it would be of the form $\mu=-a\,R+i\,\sigma$. Plugging this into the equation which was previously solved with the Lambert W-function gives
$$ (-a\,R+i\,\sigma)\,e^{-a\,R+i\,\sigma} = b\,R\,e^{- a\,R}. $$
Normally the real and imaginary part of the solution of the solution of the Lambert W-function is proportional to the branch number and the logarithm of the branch number respectively. So in order to give the solution of the principal branch $|\sigma|$ should be small. When $\sigma=0$ this gives $-a\,R=b\,R$, which matches with the edge in the image with $b\,R>-1$. When $\sigma\neq0$ it can be simplified to
$$ (-a\,R+i\,\sigma)\,(\cos(\sigma)+i\sin(\sigma)) = b\,R. $$
Solving this for $a\,R$ and $b\,R$ gives
\begin{align} a\,R &= \frac{\sigma}{\tan(\sigma)}, \\ b\,R &= \frac{-\sigma}{\sin(\sigma)}, \end{align}
which matches with the edge in the image with $b\,R<-1$ when $\sigma\in[-\pi,\pi]$.
Another stability analysis tool that could be used is the Nyquist Criterion. Namely one could use the following transfer function as an open loop transfer function
$$ G(s) = \frac{b\,e^{-R\,s}}{a - s}. $$
The closed loop (and thus also your system) is stable if the Nyquist contour makes one counter clockwise encirclement around the minus one point when the real part of $a$ is positive, or when it makes no counter clockwise encirclements around the minus one point when the real part of $a$ is not positive.
For a given $a$ and $R$ it is then relatively easy to read for what values of $b$ the system will be stable. Namely $b$ is just a gain, so the gain margin could be used. The Nyquist contour can be obtained by evaluating $G(s)$ at $s=j\,\omega$ with $\omega\in[-\infty,\infty]$, which gives
$$ G(j\,\omega) = \frac{b}{a^2 + \omega^2} \left(\cos(R\,\omega) - j \sin(R\,\omega)\right) \left(a + j\,\omega\right). $$
The points of interest of the stability analysis are the places where the the Nyquist contour crosses the real axis, so when the imaginary part of $G(j\,\omega)$ is zero. At these frequencies of real axis crossings the stability edge cases will be at some of the cases where the real part equals minus one. Solving for this gives
\begin{align} a\,R &= \frac{R\,\omega}{\tan(R\,\omega)}, \\ b\,R &= \frac{-R\,\omega}{\sin(R\,\omega)}. \end{align}
For $R\,\omega\in[-\pi,\pi]$ this again gives the same stability edge with $b\,R<-1$ in the image of the previous derivation. However, I did not find a way of using Nyquist to also show the same stability edge with $b\,R>-1$.
After replacing $\lambda$ with $\lambda/R$, $a$ with $a/R$ and $b$ with $b/R$, it may be assumed that $R = 1$.
Set $u(t) = X(t) - X^\ast$, write down the equation for $u$, and try $u(t) = e^{\lambda t}$. This results in the characteristic equation $$ \lambda = a + be^{-\lambda} $$ The system is stable if this transcendental equation has only solutions with negative real parts. A sufficient condition is $a + |b| < 0$.
To find a relation between $a$ and $b$ that separates stable and unstable solutions, consider a purely imaginary $\lambda = i \omega$. The real and imaginary parts of the equation then become $$ 0 = a + b \cos \omega, \quad \omega = - b \sin \omega $$
implying $\tan \omega = \frac{\omega}{a}$ and $ b = - \frac{\omega}{\sin \omega} = - a \sec \omega$.
You can now find marginally stable solutions for any $a$ by solving the first equation for $\omega$ and then finding $b$ from the second equation.