tdcpy.stability package#

The tdcpy.stability package contains functions for computing the characteristic roots, and spectral abscissa of time-delay systems.

Submodules#

tdcpy.stability.bounds module#

Set of functions to compute the lower bounds and upper bounds#

Implemented functions:

  1. lower_bound -> returns the lower bound at x, given (x,epsilon,gamma)

  2. upper_bound -> returns the upper bound at x, given (x,epsilon,gamma)

Original matlab files:

  • tds-control/code/+tds_control/+stability/lowerbound.m

  • tds-control/code/+tds_control/+stability/upperbound.m

tdcpy.stability.bounds.lower_bound(x, epsilon=0.15, gamma=1e-06)#

Calculates lower bound

Parameters:

  • x (float) – point to calculate lower bound at

  • epsilon (float) – small positive parameter, default 1e-2

  • gamma (float) – small positive parameter, default 1e-2

Returns:

lower bound at x

Return type:

float

Notes

  1. for |x| large enough, lower bound is (1-epsilon)*x or (1+epsilon)*x

  2. for |x| small enough, lower bound is x - gamma

Examples

>>> from tdcpy.stability.bounds import lower_bound, upper_bound
>>> lower_bound(100.0)
85.0
>>> lower_bound(-100.0)
-115.0
>>> lower_bound(0.0)
-1e-06
>>> upper_bound(100.0)
115.0
>>> upper_bound(-100.0)
-85.0
>>> upper_bound(0.0)
1e-06
tdcpy.stability.bounds.upper_bound(x, epsilon=0.01, gamma=0.01)#

Calculates upper bound

Parameters:

  • x (float) – point to calculate upper bound at

  • epsilon (float) – small positive parameter, default 1e-2

  • gamma (float) – small positive parameter, default 1e-2

Returns:

upper bound at x

Return type:

float

Notes

  1. for |x| large enough, upper bound is (1+epsilon)*x or (1-epsilon)*x

  2. for |x| small enough, upper bound is x + gamma

Examples

>>> from tdcpy.stability.bounds import lower_bound, upper_bound
>>> lower_bound(100.0)
85.0
>>> lower_bound(-100.0)
-115.0
>>> lower_bound(0.0)
-1e-06
>>> upper_bound(100.0)
115.0
>>> upper_bound(-100.0)
-85.0
>>> upper_bound(0.0)
1e-06

tdcpy.stability.characteristic_roots module#

Characteristic roots#

Set of functionalities connected to characteristic roots of DDAE

Implemented functions:

  1. roots_ddae: characteristic roots of DDAE

  2. rightmost_root: right most root and necessary things for gradient

TODO: 1. split rootd_ddae into roots_ddae_rhp and roots_ddae_region and move higher logic into high level API (tdcpy.roots)

class tdcpy.stability.characteristic_roots.RootsInfo(discretization, gamma_r_exceeds_one, index_exceeds_one, discretization_eigenvalues, max_size_evp_enforced, newton_inital_guesses, newton_final_values, newton_residuals, newton_unconverged_initial_guesses, newton_large_corrections)#

Bases: tuple

discretization#

Alias for field number 0

discretization_eigenvalues#

Alias for field number 3

gamma_r_exceeds_one#

Alias for field number 1

index_exceeds_one#

Alias for field number 2

max_size_evp_enforced#

Alias for field number 4

newton_final_values#

Alias for field number 6

newton_inital_guesses#

Alias for field number 5

newton_large_corrections#

Alias for field number 9

newton_residuals#

Alias for field number 7

newton_unconverged_initial_guesses#

Alias for field number 8

class tdcpy.stability.characteristic_roots.RightmostRootInfo(M, DM, u, v, found, max_size_evp_enforced)#

Bases: tuple

DM#

Alias for field number 1

M#

Alias for field number 0

found#

Alias for field number 4

max_size_evp_enforced#

Alias for field number 5

u#

Alias for field number 2

v#

Alias for field number 3

tdcpy.stability.characteristic_roots.roots_ddae(E, A, hA, r, **kwargs)#

Computes the roots of DDAE represented via E, A, hA in region r

TDS is assumed to be defined via E, A, hA as

\[E \dot{x}(t) = \sum\limits_{k=1}^N x(t-h_{A,k}) \ldots \qquad (1)\]
region is defined either by (1) r is float:

region = {z in C: Re(z) >= r && Im(z) >= 0}

or r is list of 4 floats [Re_min, Re_max, Im_min, Im_max]
region = {z in C: Re(z) in [Re_min, Re_max] and

Im(z) in [Im_min, Im_max]}

Parameters:

  • E (array) – 2d array defining LHS of DDAE

  • A (array) – 3d array defining the Ai matrices on the RHS

  • hA (array) – 1d array defining the delays associated with A

  • r (float) – definition of region, either float (half-plane) or list of 4 floats (rectangle)

  • **kwargs –

    discretization (int): discretization for discretizing DDAE into DAE,

    optional, default None, if not specified, heuristic will be used to obtain sufficient discretization

    max_size_evp (int): maximum allowed size of eigenvalue problem (EVP)

    optional, default 600

    base_delay (float): base delay in case delays are commensurate,

    optional, default None, used in discretization heuristic

    cd (float): c_D value of provided DDAE, optional, default None, if not

    provided condition c_D < r will be checked in case of RHP region and heuristic for obtaining discretization will be envoked if condition is not satisfied, if provided, condition c_D < r will be checked and warning will be raised if not satisfied. If you want to skip this check, provide any value of cd smaller then r, setting cd=-np.inf makes sure check will never be performed.

Returns:

A tuple containing:

crarray

vector of obtained roots

roots_infoRootsInfo

RMR metadata containing:

discretizationint

discretization used for obtaining EVP

gamma_r_exceeds_onebool

flag indicating gamma(r) > 1

index_exceeds_onebool

flag that index exceeds one

discretization_eigenvaluesarray

eigenvalues of EVP

max_size_evp_enforcedbool

flag if maximum size of EVP was enforced

newton_inital_guessesarray

roots before newton corrections

newton_final_valuesarray

roots after newton corrections

newton_residualsarray

newton residuals

newton_unconverged_initial_guessesarray

mask of unconverged newton initial guesses

newton_large_correctionsarray

mask of “large” corrections

Return type:

tuple

Notes

  1. assumes non-empty, real E, A, hA

  2. assumes dimensions E.shape[:2] == A.shape[:2] and A.shape[2] == hA.shape[0]

  3. if hA[0] > 0, prepends zero-delay term to hA and A

  4. if region is half-plane and no roots found, warns user

  5. if region is rectangle and no roots found, warns user

  6. uses Newton corrections to refine roots obtained via spectral discretization

  7. if discretization is not provided, heuristic is used to obtain sufficient discretization

Examples

>>> import numpy as np
>>> from tdcpy.stability.characteristic_roots import roots_ddae
>>> E = np.array([[1,0],[0,1]])
>>> A = np.zeros(shape=(2,2,2))
>>> A[:,:,0] = np.array([[0,1],[0,0]])
>>> A[:,:,1] = np.array([[0,0],[1,0]])
>>> hA =  np.array([0,1])
>>> cr, info = roots_ddae(E,A,hA,r=-0.1)
>>> cr
array([0.61803399+0.j, -1.61803399+0.j, -0.30901699+0.95105652j,
       -0.30901699-0.95105652j])
>>> info.discretization
10
>>> info.gamma_r_exceeds_one
False
>>> info.index_exceeds_one
False
tdcpy.stability.characteristic_roots.rightmost_root(E, A, hA, r=0.0, **kwargs)#

Find the rightmost root (RMR) of DDAE represented via E, A, hA

TDS is assumed to be defined via E, A, hA as

\[E \dot{x}(t) = \sum\limits_{k=1}^N A_k x(t-h_{A,k})\]
Parameters:

  • E (array) – 2d array defining LHS of DDAE

  • A (array) – 3d array defining the Ai matrices on the RHS

  • hA (array) – 1d array defining delays associated with A

  • r (float) – definition of complex half-plane, also discretization is perfomed around this point, optional, default 0.0

  • **kwargs –

    residual_maxfloat

    roots with norm=||M(s) @ v|| bellow this threshold will be considered as correct solutions, default 1e-6

Returns:

A tuple containing:

root_starcomplex

rightmost root

root_infoRightmostRootInfo

RMR metadata containing: M : array

characteristic matrix evaluated at root_star

DMarray

derivative of characteristic matrix evaluated at root_star

uarray

left eigenvector associated with root_star

varray

right eigenvector associated with root_star

foundbool

flag if RMR succesfully found

max_size_evp_enforcedbool

flag if maximum size of EVP was enforced

Return type:

tuple

tdcpy.stability.discretization_heuristic module#

Discretization heuristic#

References

[1] Wu, Z., & Michiels, W. Reliably computing all characteristic roots of

delay differential equations in a given right half plane using a spectral method. Journal of Computational and Applied Mathematics, 236(9), 2012, pp. 2499-2514.

tdcpy.stability.discretization_heuristic.grid_generator(m, n_grid_points=10)#

Generator for m-dimensional space grid points in [0, pi] x [-pi, pi]^{m-1}

Parameters:

  • m (int) – number of dimensions

  • n_grid_points (int) – number of grid points in each dimension (default 10)

Returns:

generator yielding m-dimensional grid points

Return type:

generator

Notes

  1. first dimension is in [0, pi]

  2. other dimensions are in [-pi, pi]

  3. total number of points generated is n_grid_points**m

Examples

>>> for point in grid_generator(2, 3):
...     print(point)
[0.         0.        ]
[0.         3.14159265]
[0.         -3.14159265]
[1.57079633 0.        ]
[1.57079633 3.14159265]
[1.57079633 -3.14159265]
[3.14159265 0.        ]
[3.14159265 3.14159265]
[3.14159265 -3.14159265]
tdcpy.stability.discretization_heuristic.commensurate_gk(E, B, C, n_k, grid_points=20)#

TODO

Return type:

ndarray[tuple[Any, …], dtype[_ScalarT]]

tdcpy.stability.discretization_heuristic.incommensurate_gk(E, B, C, tau, grid_points=20)#

TODO - in original implementation, only 3 non zero delays allowed

Return type:

ndarray[tuple[Any, …], dtype[_ScalarT]]

tdcpy.stability.discretization_heuristic.compute_n_rhp(E, B, C, tau, basic_delay=None, **kwargs)#
Parameters:
Return type:

int

kwargs:

n_minimal (int): minimal degree of discretization, default 8 n_grid (int): number of grids point in interval [0, pi] for discretizing

Psi, default 20

Returns:

number of discretization points necessary

Return type:

n (int)

Parameters:
See:
[1] Wu, Z., & Michiels, W. Reliably computing all characteristic roots of

delay differential equations in a given right half plane using a spectral method. Journal of Computational and Applied Mathematics, 236(9), 2012, pp. 2499-2514.

tdcpy.stability.discretization_heuristic.compute_n_rect(region, tau_max, **kwargs)#

Computes the deegree of the spectral discretization for rectanuglar region in complex plane

Computes the degree of the spectral discretisation and a shift of the origin such that exp(-s*taum) is sufficiently well approximated in the rectanuglar region in complex plane:

region[1]+1j*region[4] --------- region[2]+1j*region[4]
        |                                |
        |                                |
region[1]+1j*region[3] --------- region[2]+1j*region[3]
Parameters:

  • region (tuple) – tuple representing rectangular region in complex plane, example: (-10, 0, 0, 100) represents region D from complex plane D={z from C: -10<=Re(z)<=0, 0<=Im(z)<=100}

  • tau_max (float) – maximal delay, has to be >= 0

  • kwargs

Returns:

  • tuple containing –

    • n (int): number of discretization points necessary

    • origin (complex): origin TODO

  • See –

    [1] Wu, Z., & Michiels, W. (2012) Reliably computing all

    characteristic roots of delay differential equations in a given right half plane using a spectral method. Journal of Computational and Applied Mathematics, 236(9), pp. 2499-2514.

Return type:

int

tdcpy.stability.gamma_r module#

Computation of gamma(r, DIFF)#

DIFF - delay difference equation represented via matrices (3D array) and delays (1D array)

Notes

  1. MemoizeJac decorator can be used because of (i) keep implementation as

    simillar as possible to MATLAB® and (ii) to optimize calculation of function value and jacobian as some operations are shared

class tdcpy.stability.gamma_r.GammaInfo(th, M, s, u, v)#

Bases: tuple

M#

Alias for field number 1

s#

Alias for field number 2

th#

Alias for field number 0

u#

Alias for field number 3

v#

Alias for field number 4

tdcpy.stability.gamma_r.rho_normalized_diff(DD, hDD, r, theta)#

Calculates spectral radius of matrix M

Matrix M is given as sum of rotations of normalized delay difference equation.

Parameters:

  • DD (ndarray) – Coefficient matrices packed into 3D array of shape (n, n, m). Note that coefficients for x(t) are assumed to be identity matrix and therefore omitted (see functions for converting DDAE to delay difference equation and normalizing).

  • hDD (ndarray) – Delays represented by 1D array of shape (m,). Note that delay 0 is omitted.

  • r (float) – Point from complex plane.

  • theta (ndarray) – Non-empty 1d array of rotations of shape (m-1,) in [0, 2*pi).

Returns:

Spectral radius matrix.

Return type:

ndarray

Notes

Does not check for valid inputs.

tdcpy.stability.gamma_r.func(x, DD, hDD, r, v0)#

Calculates value and jacobian of the vector function F(x)

vector x, shape=(4*ndiff + 2 + (m-1), ):

x = [Re(v), Im(v), Re(u), Im(u), Re(lambda), Im(s), th]

function F(x):

M * v - s * v = 0 u^{H} * M - s * u^{H} = 0 u^{H} * v - 1 = 0 v0^{H} * v - 1 = 0 Im(conj(lambda)*(u^{H}*DD{k}*v)*exp(-r*hDD(k))*exp(1j*theta(k)) = 0 for k = 1,…,m

with th = [0; th_v], v0 a normalization vector and

M = DD[0]*exp(-r*hDD[0])*exp(1j*th[1]) + … + DD[m]*exp(-r*hDD[m])*exp(1j*th[m])

using fsolve (#optim. variables: 4*ndiff + 2 + (m-1), #constraints: 4*ndiff + 2 + (m-1) ).

Parameters:

  • x (ndarray) – 1D array of shape (4*ndiff + 2 + (m-1), ) representing optimization variables

  • DD (ndarray) – Coefficient matrices packed into 3D array of shape (n, n, m).

  • hDD (ndarray) – Delays represented by 1D array of shape (m,).

  • r (float) – Point in complex plane.

  • v0 (ndarray) – Normalization vector of shape (n, ).

Returns:

  • tuple containing

  • y (ndarray) – 1D array representing function F evaluated at x

  • jac (ndarray) – Jacobian of F evaluated at x, shape (2*(2*ndiff+2) + n_opt, 2*(2*ndiff+1) + n_opt)

Notes

  1. shape of x (n_opt, ), n_opt = 4*ndiff + 2 + (m-1)

  2. shape of jac (2*(2*n_diff+2) + n_opt, 2*(2*n_diff+1) + n_opt), n_diff … dimension of the state vector of the delay-difference equation

tdcpy.stability.gamma_r.l2_func(*args)#

returns L2 of F(x) and its jacobian

tdcpy.stability.gamma_r.theta_generator(n_opt, n_theta=10, is_real=False)#

Generator for theta grid points in [0, 2*pi)^{m}

Parameters:
tdcpy.stability.gamma_r.gamma_normalized_diff(DD, hDD, r, **kwargs)#

Computes gamma(r) of the normalized delay difference equation

Normalized delay difference equation takes form

0 = x(t) + DD[0] * x(t-hDD[0]) + … + DD[m-1]*x(t-hDD[m-1]), (1)

where m == len(DD) == len(hDD).

The gamma(r) of (1) is given by the following optimization problem:
find maximum theta from [0, 2*pi)^m of expression:

rho( SUM for all k DD[k]*exp(-r*hDD[k])*exp(1j*theta[k]) ) (2)

where rho(.) is spectral radius of its matrix argument.

Parameters:

  • DD (ndarray) – Coefficient matrices packed into 3D array shaped (n,n,m), note that coefficients for x(t) are assumed to be identity matrix and therefore omitted (see functions for converting DDAE to delay difference equation and normalizing).

  • hDD (ndarray) – Delays represented by 1D array shaped (m,), note that delay 0 is omitted.

  • r (float) – Point from complex plane.

  • kwargs –

    n_theta (int): theta discretization, has to be > 0, default 10 correction (bool): if correction is applied, default True scipy_root_method (str): scipy.optimize.root method, default ‘lm’,

    i.e. Levenberg-Marquardt algorithm, note: carefull, not all methods attempt to solve problem

    scipy_root_tol (float): Tolerance for termination. For detailed

    control, use scipy_root_options, default None

    scipy_root_callback (function): Optional callback function. It is

    called on every iteration as callback(x, f) where x is the current solution and f the corresponding residual. For all methods but ‘hybr’ and ‘lm’.

    scipy_root_options (dict): a dictionary of solver options (method),

    default None

Returns:

  • gamma (float): quantity gamma(r, DD, hDD)

  • info (GammaInfo): gamma metadata

Return type:

tuple containing

Notes

1. for all kwargs starting with ‘scipy_*’ check the following documentation https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.root.html 2. uses grid search + optional correction via scipy.optimize.root

Examples

>>> import numpy as np
>>> from tdcpy.stability.gamma_r import gamma_normalized_diff
>>> DD = np.zeros((2,2,2), dtype=complex)
>>> DD[:,:,0] = np.array([[0.7, 0], [0, 0.7]])
>>> DD[:,:,1] = np.array([[0.1, 0], [0, 0.1]])
>>> hDD = np.array([1.0, 2.0])
>>> r = 0.5 + 0.5j
>>> gamma, info = gamma_normalized_diff(DD, hDD, r, n_theta=5, correction=True)
>>> print(f"gamma: {gamma}")
gamma: 0.4613594059159876
>>> print(f"theta: {info.th}")
theta: [0.  0.5]
>>> print(f"dominant eigenvalue: {info.s}")
dominant eigenvalue: (0.40488096939597285-0.22118748167138752j)
>>> print(f"right eigenvector: {info.v}")
right eigenvector: [0.-0.j 1.-0.j]
>>> print(f"left eigenvector: {info.u}")
left eigenvector: [ 0.        +0.j         -0.68163876-0.73168887j]
tdcpy.stability.gamma_r.gamma_diff(D, hD, r, **kwargs)#

Computes gamma(r) of delay difference equation (DIFF)

Delay difference equation takes form

0 = D[0]*x(t) + D[1] * x(t-hD[1]) + … + DD[m-1]*x(t-hDD[m-1]), (1)

where m == len(DD) == len(hDD).

quantity gamma(r; D, hD) is then: 1. gamma(r; D, hD) = 0 IF number of delays (vector hD) is less then 2 2. obtained via predictor corrector approach, i.e.

2a. normalize DIFF (multiply equation (1) by inverse of D[0]) and

omit first delay = 0 and first normalized matrix = identity

2b. call gamma_diff_normalized

Parameters:

  • D (array) – Coefficient matrices packed into 3D array shaped (n, n, m).

  • hD (array) – Delays represented by 1D array shaped (m,).

  • r (float) – Point from complex plane.

  • kwargs – kwargs passed into gamma_diff_normalized function.

Returns:

  • gamma (float): quantity gamma(r, D, hD)

  • info (GammaInfo): gamma metadata

Return type:

tuple containing

Notes

  1. if compressed version of DIFF contains 2 or more delays, invertibility of D[0] is assumed.

  2. DIFF representation (D, hD) can be emtpy, result will be gamma(r; D, hD) = 0.0. Test for emptyness is hD.size == 0.

  3. for r = 0.0, quantity gamma(r; D, hD) DOES NOT depend on the delays, see implementation of gamma_diff_normalized

Examples

>>> import numpy as np
>>> from tdcpy.stability.gamma_r import gamma_diff
>>> D = np.zeros(shape=(2,2,3))
>>> D[:,:,0] = np.array([[1.0, 0.0], [0.0, 1.0]])
>>> D[:,:,1] = np.array([[0.5, 0.0], [0.0, 0.5]])
>>> D[:,:,2] = np.array([[0.2, 0.0], [0.0, 0.2]])
>>> hD = np.array([0.0, 1.0, 2.0])
>>> r = 0.0
>>> gamma_val, gamma_info = gamma_diff(D, hD, r)
>>> print(gamma_val)
0.7
>>> print(gamma_info)
GammaInfo(th=array([0., 0.]), M=array([[0.7.+0.j , 0. +0.j ],
        [0. +0.j , 0.7+0.j ]]), s=(0.7+0j), u=array([0.+0.j, 1.+0.j]),
        v=array([0.-0.j, 1.-0.j]))

tdcpy.stability.newton module#

Newton method for increasing precision of roots#

tdcpy.stability.newton.newton_correction(roots0, E, A, hA, **kwargs)#

Newton corrections applied initial guess roots0

root s (element of roots0) is an initial guess of the eigenvalue of characteristic matrix M(s) and its derivative dM(s) .. math:

M(s) = E*s - A_0*e^{-s*h_{A,0}} - ... - A_{mA}*e^{-s*hA_{mA}}
dM(s) = E + hA_0*A_0*e^{-s*hA_0} + ... + hA_{mA}*A_{mA}*e^{-s*hA_{mA}}

and we look for a solution of system n+1 non-linear equations .. math:

M(s) * v      ==  0121
v0.H * v - 1  ==  0

with jacobian .. math:

J = egin{bmatrix} M(s) & dM(s)*v \
     v0.H &0      \end{bmatrix}

Update rule for newton corrections stands as .. math:

[dv.T, ds].T = PINV( J ) @ [M(s)*v, v0.H * v - 1].T

v = v - dv
s = s - ds
Parameters:

  • roots0 (npt.NDArray) – 1D vector of initial guesses for roots.

  • E (npt.NDArray) – E matrix from TDS representation shaped (n, n).

  • A (npt.NDArray) – A matrices from TDS representation shaped (n, n, mA).

  • hA (npt.NDArray) – hA vector of delays from TDS representation shaped (mA,).

  • **kwargs (dict, optional) –

    Additional arguments:

    • return_residuals (bool): Whether to also return residuals or not. Default False.

    • tol (float): Absolute tolerance. Default 1e-10.

    • max_iterations (int): Maximum number of newton iterations. Default 20.

Returns:

  • roots (npt.NDArray): 1D array of roots with improved precision.

  • residuals (npt.NDArray): 1D array of residuals.

  • converged_mask (npt.NDArray): 1D bool array indicating if root converged.

  • large_correction_mask (npt.NDArray): 1D bool array indicating if correction is large.

Return type:

tuple containing

Notes

  1. The expected shapes of input arrays:

    • E … (n,n)

    • A … (n,n,mA)

    • hA … (mA,)

    shapes are not checked to speed up computation.

  2. It is possible some corrections won’t converge (residual check)

  3. It is possible some corrections converge, hoewever corrections are

    “large” (see “newton method - basin of attraction”), these corrections should be taken with a grain of salt because root can converge to root which is already present in solution.

tdcpy.stability.spectral_abscissa module#

Spectral abscissa#

We distinguish between three types:

  1. spectral abscissa (denote alpha)

  2. spectral abscissa of delay-difference equation (denote cd)

  3. strong spectral abscissa := MAX(alpha, cd)

Functions#

  • spectral_abscissa: placeholder for classic spectral abscissa routines.

  • func: evaluate f(r)=gamma(r)-1 and its derivative at a scalar r.

  • CdRootProblem: root-problem wrapper object used by scipy.optimize.root.

  • spectral_abscissa_diff: compute the delay-difference contribution cd.

References

[1] Michiels, W. and Niculescu, S.I. (2014). Stability, control, and computation for time-delay systems: an eigenvalue-based approach. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA

tdcpy.stability.spectral_abscissa.spectral_abscissa()#
tdcpy.stability.spectral_abscissa.func(r, DD, hDD, gamma_kwargs)#

Evaluate f(r) = gamma(r) - 1 and its derivative.

Parameters:

  • r (float) – Real number defining right-half complex plane for characteristic roots computation

  • DD (npt.NDArray) – Normalized delay-difference matrices defining the operator.

  • hDD (npt.NDArray) – Delays corresponding to the slices of DD (shape must match DD.shape[2]).

  • gamma_kwargs (dict) – Additional keyword arguments forwarded to gamma_normalized_diff.

Returns:

  • fval (float) – The scalar value gamma(r) - 1.

  • df (float) – The derivative df/dr evaluated at r.

Notes

The implementation calls gamma_normalized_diff to obtain the normalized spectral radius and associated metadata, then computes the derivative using the expression derived from the numerator involving the left/right eigenvectors returned by the helper.

class tdcpy.stability.spectral_abscissa.CdRootProblem(DD, hDD)#

Bases: object

root_options = {'jac': True, 'scipy_root_callback': None, 'scipy_root_method': 'hybr', 'scipy_root_options': None, 'scipy_root_tol': None}#
gamma_kwargs = {}#
fun(r)#
Return type:

tuple[float, float]

tdcpy.stability.spectral_abscissa.spectral_abscissa_diff(DD, hDD, **kwargs)#

Compute the strong spectral abscissa of a normalized delay-difference equation.

This routine follows the formulation in Michiels & Niculescu (2014) and locates roots of f(r)=gamma(r)-1 using scipy.optimize.root to determine the delay-difference contribution to the strong spectral abscissa.

Parameters:

  • DD (npt.NDArray) – Normalized delay-difference matrices of shape (n, n, m).

  • hDD (npt.NDArray) – Delays corresponding to the third dimension of DD (shape (m,)).

  • **kwargs (dict) –

    Optional keyword arguments:

    • cd0 (float): initial guess for the strong spectral abscissa, default 0.0.

    • scipy_root_method (str): method for scipy.optimize.root, default 'hybr'.

    • scipy_root_tol (float): tolerance for the root solver.

    • scipy_root_callback (callable): optional callback callback(x, f).

    • scipy_root_options (dict): solver-specific options.

    • gamma_kwargs (dict): kwargs forwarded to gamma_normalized_diff.

Returns:

  • cd (float) – The strong spectral abscissa of the associated delay-difference equation. May be -np.inf or np.inf in degenerate cases.

  • info (CDInfo) – Metadata returned from gamma_normalized_diff corresponding to the root with smallest residual (or None when unavailable).

Return type:

tuple[float, GammaInfo]

Notes

  • If DD.shape[1] == 0 the function returns -np.inf to indicate no meaningful contribution from the delay-difference operator.

  • If the root-finding fails but the zero-frequency gamma satisfies gamma(0) >= 1, the routine treats the strong spectral abscissa as np.inf; otherwise -np.inf is used. In the current implementation a failed root solve raises NotImplementedError instead of returning a metadata-rich result.

References

[1] Michiels, W. and Niculescu, S.I. (2014). Stability, control, and

computation for time-delay systems: an eigenvalue-based approach. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA

Module contents#