Example 2.12 - Stability analysis of a Smith predictor with delay mismatch

We will follow example 2.12 from [Appeltans and Michiels, 2023] Section 2.6.2. In there, for given plant and controller

\[H(s) = \frac{1}{s + 1}e^{-s\tau} , C(s) = \frac{s}{2} + 2,\]

the classical smith predictor is constructed. Assuming the delay mismatch \(\delta\) we arrive at the quasipolynomial

\[D_{cl}(s) = \frac{3}{2} s + 3 + (\frac{s}{2} + 2) e^{-s\tau} + (\frac{-s}{2} - 2) e^{-s(\tau + \delta)},\]

which we will analyse for stability in \((\tau, \delta)\)-parameter space.

/home/runner/work/tdcpy/tdcpy/src/tdcpy/common/quasipoly.py:136: SyntaxWarning: invalid escape sequence '\d'
  \dot{x}(t) = A_0*x(t - hA_0) + ... + A_{mA}*x(t - hA_{mA})

import numpy as np
import tdcpy
import tdcpy.plot
import matplotlib.pyplot as plt
from tdcpy.common.quasipoly import qp_to_ndde
from tdcpy.common.quasipoly import compress_qp, qp_to_ndde


tau, delta = 1., 0.5

coeffs = np.array([[3.,1.5],[2.0,0.5],[-2.0,-0.5]])
delays = np.array([0.,tau,tau+delta])

A, hA, H, hH = qp_to_ndde(coeffs,delays,ascending=True)

ndde = tdcpy.NDDE(A=A,hA=hA,H=H,hH=hH)

tau_grid = np.linspace(0, 8, 201)
delta_grid = np.linspace(-8, 10, 451)
Z = np.zeros((len(delta_grid), len(tau_grid)))

# for i2 in range(0,len(tau_grid)-1):
#     tau = tau_grid[i2]
#     for i1 in range(0,len(delta_grid)):
#         delta = delta_grid[i1]
#         if np.abs(tau+delta)<1e-8:
#             # case: tau+delta = 0
#             if np.abs(tau)<1e-8:
#                 tau = 1e-8
#             hH[0],hH[1] = tau, 1e-8
#             hA[1],hA[2] = tau, 1e-8
#             ndde = tdcpy.NDDE(H=H,hH=hH,A=A,hA=hA)
#             Z[i1,i2] = tdcpy.strong_spectral_abscissa(ndde)
#         elif tau + delta < 0:
#             # case: "real" delay cannot be negative
#             Z[i1,i2] = -np.inf
#         else:
#             hH[0], hH[1] = tau, tau+delta
#             hA[1], hA[2] = tau, tau+delta
#             ndde = tdcpy.NDDE(H=H,hH=hH,A=A,hA=hA)
#             Z[i1,i2] = tdcpy.strong_spectral_abscissa(ndde)

# X, Y = np.meshgrid(tau_grid,delta_grid)
# plt.contour(X,Y,Z,[0])
# plt.plot(plt.xlim(),[0,0],'k-.')
# plt.show()