# Wrong stability results when using Padé approximation

The Padé approximation of the exponential function, $$F(s) = e^{-\tau s}$$, is used often in control theory. I wonder whether its use can lead to erroneous results regarding the stability properties of the closed-loop system.

Let's say we have a dynamical system with transfer function

$$G(s) = \frac{P(s)}{Q(s)}e^{-\tau s},$$

with $$\tau>0$$, where $$P$$ and $$Q$$ are polynomials and the degree of $$P$$ is strictly smaller than the degree of $$Q$$. The corresponding closed-loop transfer function is defined as

$$G_{CL}(s) = \frac{G(s)}{1+G(s)}.$$

If we approximate $$e^{-\tau s}$$ with the Padé approximation of order $$[1,1]$$ we obtain

$$G_{\mathrm{approx}}(s) = \frac{P(s)}{Q(s)}\frac{1-\frac{\tau}{2}s}{1+\frac{\tau}{2}s}.$$

Then, $$G_{CL, \mathrm{approx}} = \frac{G_{\mathrm{approx}}(s)}{1+G_{\mathrm{approx}}(s)}$$ becomes a rational function and we can use an algebraic criterion to tell whether it is stable. However, I guess that this does not imply that $$G_{CL}$$ is stable, but I have not been able to find a counterexample.

I tried to take $$G$$ to be a first and a second order system and tried with low and high values of $$\tau$$, but I was not able to come up with a counterexample where Routh's criterion on $$G_{CL, \mathrm{approx}}$$ and the Bode or Nyquist stability criteria for the actual transfer function lead to different conclusions regarding the stability of the closed-loop system.

Usually both the delay and the right half plane zero constrain how high you can place the closed loop bandwidth. However a delay really put a hard upper bound on the bandwidth because the phase drop increases linearly in frequency (it looks exponentially when the frequency is on a logarithmic scale). For the right half plane zero it is still possible to get an arbitrary high bandwidth, because the phase only drops at most 90 degrees, however high bandwidths do yield worse and worse margins, for example see my answer here. From this it is also possible to construct an example where the system with the Padé approximation is stable in closed-loop but with a time delay it is not, namely

$$G(s) = \frac{100\,s + 90}{s^2 + 110\,s}e^{-2\,s}, \tag{1}$$

so

$$G_\mathrm{approx}(s) = \frac{100\,s + 90}{s^2 + 110\,s}\frac{1 - s}{1 + s}. \tag{2}$$

The corresponding Nyquist diagram can be seen below, from which one can conclude that $$(2)$$ would be closed-loop stable, but $$(1)$$ would not.

However I would like to state this should only be the case when you choose a bandwidth high compared to $$1/\tau$$. Because in my example $$1/\tau=0.5$$ while the bandwidth lies close to 2 rad/s. I define the bandwidth as the frequency where the open-loop crosses the 0 dB line for the first time, this definition is sometimes also called the crossover frequency. Normally one would already design a controller which would not have a high bandwidth compared to $$1/\tau$$ if one wants to have reasonable phase and gain margins, even when using the Padé approximation for the delay in this design process. For example adding a gain of 0.848 would give $$G_\mathrm{approx}(s)$$ a phase margin of 45 degrees (however the gain margin is still only 2.27 dB). With this gain $$G(s)$$ is stable as well, however its margins a still quite a bit different, namely a phase margin of 15 degrees and a gain margin of 0.455 dB. It can be noted that the bandwidth is still 1.1 rad/s, but even if you do not have a relatively high bandwidth I would still recommend checking whether or not requirements for the margins are met while using the delay instead of the Padé approximation.

Another way of looking at this is considering the phase margin of the system without delay and look at the phase shift induced by the actual delay and the Padé approximation. The phase margin of the system is denoted with $$\phi$$ at a crossover frequency $$\omega_c$$. The phase shift at $$\omega_c$$ of the actual delay can be shown to be $$\tau\,\omega_c$$, so the maximum delay the system can actually handle can be calculated with

$$\tau_\mathrm{max} = \frac{\phi}{\omega_c}. \tag{3}$$

The phase shift at $$\omega_c$$ of a $$[1,1]$$ Padé approximation can be shown to be $$2\tan^{-1}\left(\tau\,\omega_c/2\right)$$, so the predicted maximum delay approximation the system can handle can be calculated with

$$\tau_\mathrm{max,approx} = \frac{2\tan\left(\frac{\phi}{2}\right)}{\omega_c}. \tag{4}$$

The Taylor series of $$2\tan\left(\frac{\phi}{2}\right)$$ is

$$2\tan\left(\frac{\phi}{2}\right) = \phi + \frac{\phi^3}{12} + \frac{\phi^5}{120} + \mathcal{O}(\phi^7), \tag{5}$$

so when comparing $$(4)$$ to $$(3)$$ it can be seen that the larger $$\phi$$ the worse the predicted maximum delay according to $$[1,1]$$ Padé approximation is, since the Taylor series shown that it always over estimates the delay it can handle (only this over estimation gets really small when $$\phi$$ is small).

From this one can think of an even easier example

$$G(s) = \frac{1}{s} e^{-\tau\,s}. \tag{6}$$

This system has a phase margin of $$\phi=\pi/2$$ radians and a crossover frequency of $$\omega_c=1$$ rad/s. So according to $$(3)$$ this system can handle at most $$\pi/2\approx1.5708$$ seconds of delay, while according to $$(4)$$ the $$[1,1]$$ Padé approximation would still tell the system is stable up to a delay of $$2$$ seconds.

If one would use the $$[2,2]$$ Padé approximation instead it can be shown that similar to $$(4)$$ one gets

$$\tau_\mathrm{max,approx[2,2]} = \frac{\sqrt{9+12\tan^2\left(\frac{\phi}{2}\right)} - 3}{\tan\left(\frac{\phi}{2}\right)\omega_c}, \tag{7}$$

with

$$\frac{\sqrt{9+12\tan^2\left(\frac{\phi}{2}\right)} - 3}{\tan\left(\frac{\phi}{2}\right)} = \phi + \frac{\phi^5}{720} - \frac{\phi^7}{12096} + \mathcal{O}(\phi^9), \tag{8}$$

so when comparing $$(5)$$ with $$(8)$$ it can be seen that $$(8)$$ should stay closer to $$\phi$$ for a larger range of values for $$\phi$$. When again considering $$(6)$$ the $$[2,2]$$ Padé approximation predicts a maximum delay of $$1.5826$$ seconds, which is much closer to the actual value of $$1.5708$$ seconds compared to previous prediction of $$2$$ seconds. However often one normally has a smaller phase margin (of the open-loop without delay) so the $$[1,1]$$ Padé approximation can be sufficient in many cases. In most other cases the $$[2,2]$$ Padé approximation is probably enough. You do need to watch out if the open-loop crosses the 0 dB line again above the crossover frequency. In that case you could also use all those crossings for calculating the phase margin and crossover frequencies and plug those values into $$(3)$$ to find the smallest maximum delay that your system can handle.

• I haven't come across this definition of bandwidth before. Is there a reference where I can read more about it? What are its properties? I guess that is the input frequency above which the system dissipates some energy and before it, it adds energy (roughly speaking). Why does this happen if the bandwidth is larger than $1/\tau$? – Nicolas Bert Johnson May 17 at 22:30
• I used this definition for the bandwidth since this way was taught at my university. Here it is called the crossover frequency instead. Defining the bandwidth as the frequency where the closed-loop goes below -3 dB often yields bigger values then the crossover frequency. Namely crossover looks at $|G(j\omega)|=1$ while the bandwidth looks at $$\left|\frac{G(j\omega)}{1+G(j\omega)}\right|=0.71$$ so at the crossover one usually has that $|1+G(j\omega)|<1$ (true when you have small margins) so the condition for the bandwidth is not met yet. – Kwin van der Veen May 17 at 23:56
• Often both yield similar values, they only start to differ a lot when you have poor margins, for example the 2 and 1.1 rad/s in my answer turn into 142 and 118 rad/s. So you could still use as a rule of thumb that if the (other) bandwidth is below $1/\tau$ the first order Padé approximation can be used to approximate a delay. The bandwidth might be a little bit conservative compared to the crossover though, I mainly mentioned it in my answer as an order of magnitude comparison, not a hard cut-off. – Kwin van der Veen May 18 at 0:16
• @NicolasBertJohnson I added quite a bit to my answer. In the section where I use the phase margin it can also be seen that the crossover frequency is useful there as well. One can also always check closedloop stability by looking at the phase margin of the openloop (without delay) and use equation $(3)$ to find the maximum delay the system could handle. – Kwin van der Veen May 20 at 3:25

You can find a counterexample by using a BIBO-stable open-loop system with poles close to the imaginary axis and an intermediate value of $$\tau$$ (not too low, not too high). Here is an example:

$$G(s) = \frac{2}{(s+0.5)^2 + 1}e^{-0.6s}.$$

The closed-loop function $$G_c(s) = G(s)/(1+G(s))$$ is approximated by $$\tilde{G}_{c} = -\frac{80-24s}{12s^3 + 52s^2 + 31s + 130}.$$

Its poles are $$-0.0085 \pm j1.5842$$ and $$-4.3163$$ (I determined these using MATLAB), so the approximate closed-loop system is BIBO stable (you can use an algebraic criterion to verify this), but the Nyquist plot of the actual open-loop transfer function is

Indeed, for $$\tau > 0.5758$$, the closed-loop system becomes BIBO unstable, but it is not until $$\tau>0.6189$$ that you will be able to detect this using a $$[1,1]$$-order Padé approximation.

Here is some MATLAB code to determine these values:

xt = 0;
for td = 0.5:0.0001:0.65,
G = tf(2, [1 1 1.25]);
H = G*tf(1,1,'InputDelay', td);
Hpade =  G*tf([-0.5*td 1], [0.5*td 1]);

Gm = margin(H);
x = (abs(isstable(T)- (Gm > 1)));
if x~=xt,
disp(td);
end
xt = x;
end


Some notes:

1. The Padé approximation of order $$[1,1]$$ is not too good. Use a higher order approximation (but not too high as this will lead to numerical problems)
2. The the best of my knowledge, you cannot use the Padé approximation to guarantee that the actual system is stable (at least not without additional requirements)
3. Don't use Padé approximations if you want to study the stability properties of your system. Use the Bode criterion (if applicable) or the Nyquist criterion (recommended).
• I think you might have made a typo, because the Nyquist diagram seems to correspond to $$G(s)=\frac{2}{(s+0.5)^2+1}e^{-0.7\,s}$$ and not $$G(s)=\frac{2\,s}{(s+0.5)^2+1}e^{-0.7\,s}.$$ However for $\bar{G}_c(s)$ you did seemed to have used the $G(s)$ with the extra $s$ in the numerator. – Kwin van der Veen May 18 at 2:08
• @KwinvanderVeen that's a lot for noticing! I'll fix that. – Pantelis Sopasakis May 18 at 2:15
• I think your answer is still not completely correct, namely I have added to my answer a way to calculate the maximum delays and for your example system this would give $\tau_\mathrm{max,approx}=0.6188$ (close enough) and $\tau_\mathrm{max}=0.5758$ (this one does differ significantly). – Kwin van der Veen May 20 at 1:42
• @KwinvanderVeen You might be right. I determined these values by trial and error (maybe there was an error in my code). I'll double-check tomorrow). Thanks a lot. – Pantelis Sopasakis May 20 at 1:44
• @KwinvanderVeen I updated my answer with a MATLAB script to compute these values. Thanks again. – Pantelis Sopasakis May 20 at 2:29