tdcpy.stabopt package

Submodules

tdcpy.stabopt.controller_bfgs module

Delay controller design (stabilization via BFGS)

tdcpy.stabopt.controller_bfgs.design_bfgs(E, P, hP, K0, hK, B, C, **kwargs)

function for designing a controller via BFGS optimization

Parameters:

• E (array) – E matrix of closed-loop system

• P (array) – system matrix of open-loop system

• hP (array) – system delays of open-loop system

• K0 (array) – initial controller gain matrix

• hK (array) – controller delays

• B (array) – input matrix of closed-loop system

• C (array) – output matrix of closed-loop system

• kwargs (dict) – options: options for the optimization solver

Returns:

sol – Optimization result containing the optimal controller parameters and optimization information x: optimal controller parameters fun: optimal objective function value success: whether the optimization was successful message: description of the cause of the termination

Return type:

OptimizeResult

Notes

1. This function assumes that the system is retarded.

2. The optimization problem is non-convex, and the solution may depend on the initial controller parameters K0.

3. The optimization is performed using the BFGS algorithm, which is a quasi-Newton method for unconstrained optimization. The objective function is the spectral abscissa of the closed-loop system, which is computed using the function func_sa. The gradient of the objective function is computed using the function grad_sa.

4. The optimization options can be specified via the options key in kwargs.

Examples

>>> from tdcpy.stabopt.controller_bfgs import design_bfgs
>>> import numpy as np
>>> E = np.eye(2)
>>> P0 = np.array([[-1., 0.], [0., -2.]])
>>> P1 = np.array([[-0.5, 0.], [0., -0.5]])
>>> P = np.stack([P0, P1], axis=2)
>>> hP = np.array([0.,0.5])
>>> K0 = np.zeros((2,2,1))
>>> hK = np.array([0.])
>>> B = np.eye(2)
>>> C = np.eye(2)
>>> sol = design_bfgs(E, P, hP, K0, hK, B, C, options={"disp": True})

tdcpy.stabopt.controller_bfgs.design_granso(E, P, hP, K0, hK, B, C, **kwargs)

Parameters:

• TODO

• **kwargs – mask (array): masking gradient

• E (ndarray[tuple[Any, ...], dtype[_ScalarT]])

• P (ndarray[tuple[Any, ...], dtype[_ScalarT]])

• hP (ndarray[tuple[Any, ...], dtype[_ScalarT]])

tdcpy.stabopt.controller_bfgs.minimize_spectral_abscissa(ddae, order, **kwargs)

function to controller parameters for minimizing the spectral abscissa of a ddae

Parameters:

• ddae (DDAE) – DDAE object representing the open-loop system

• order (int) – Degree of the dynamic controller

• kwargs (dict) – initial (array): initial controller parameters, shape (nc+nu, nc+ny, nd) y_indices (array): indices of system outputs used for feedback u_indices (array): indices of system inputs used for control basis : controller structure mask (array): gradient mask options: stabilization options

Returns:

sol – Optimization result containing the optimal controller parameters and optimization information x: optimal controller parameters fun: optimal objective function value success: whether the optimization was successful message: description of the cause of the termination

Return type:

OptimizeResult

Notes

1. The optimization problem is non-convex, and the solution may depend on the initial controller parameters.

2. The default solver if L-BFGS-B.

3. The optimization options can be specified via the options key in kwargs.

Examples

>>> from tdcpy.stabopt.controller_bfgs import minimize_spectral_abscissa
>>> from tdcpy.ddae import DDAE
>>> A0 = np.array([[-1., 0.], [0., -2.]])
>>> A1 = np.array([[1., 0.], [0., 1.]])
>>> A = np.stack([A0, A1], axis=2)
>>> hA = np.array([0., 1.])
>>> Bu = np.array([ [-0.1],[-0.2]    ])
>>> B = np.stack([Bu],axis=2)
>>> hB = np.array([0.5])
>>> C = np.array(np.eye(2))
>>> C = np.stack([C],axis=2)
>>> hC = np.array([0])
>>> D = np.zeros(shape=(2,1,1), dtype=float)    # must be defined, otherwise error!
>>> hD = np.array([0.])
>>> ddae = DDAE(A=A,hA=hA,B=B,hB=hB,C=C,hC=hC,D=D,hD=hD)
>>> sol = minimize_spectral_abscissa(ddae, order=0, method="L-BFGS-B", options={"disp": True}, type = "barrier")

tdcpy.stabopt.controller_bfgs.find_feasible_point(E, P, hP, K0, hK, B, C, **kwargs)

function to find a feasible point for the optimization problem

for a DDAE of the form:

Parameters:

• E (array) – E matrix of closed-loop system

• P (array) – system matrix of open-loop system

• hP (array) – system delays of open-loop system

• K0 (array) – initial controller gain matrix

• hK (array) – controller delays

• B (array) – input matrix of closed-loop system

• C (array) – output matrix of closed-loop system

• kwargs – options: options for the optimization solver gamma0_threshold: threshold for gamma0 to consider a point feasible, default is 1.0

Returns:

K – feasible controller parameters, if found, otherwise None

Return type:

array

Notes

Examples

>>> import numpy as np
>>> from tdcpy.stabopt.controller_bfgs import find_feasible_point
>>> E = np.array([[1,0,0],[0,1,0],[0,0,0]])
>>> P = np.zeros((3,3,1))
>>> hP = np.array([0])
>>> K0 = np.zeros((1,3,1))
>>> hK = np.array([0])
>>> B = np.eye(3)[:,:,np.newaxis]
>>> C = np.eye(3)[:,:,np.newaxis]
>>> K_feasible = find_feasible_point(E, P, hP, K0, hK, B, C)
>>> print(K_feasible)

tdcpy.stabopt.gradients module

Functions and its gradients for stabilization

tdcpy.stabopt.gradients.func_sa(x, E, P, hP, hK, Kmask, B, C)

function for spectral abscissa and its gradient with respect to controller parameters

The system (Closed Loop) is defined as

E dxdt(t) = SUM P[i] x(t-hP[i]) + SUM B * K[j] * C x(t-hK[j]) (1)

where K = x.reshape(Kshape) are controller parameters, hK controller delays, P and hP are controlled system matrices and delays, respectively.

Parameters:

• x (array) – vectorized controller parameters

• E (array) – 2d array defining LHS of Closed Loop

• P (array) – 3d array defining the Pi matrices (controlled system)

• hP (array) – 1d array defining delays associated with P

• hK (array) – 1d array defining delays associated with K (note that x := vec(K))

• Kmask (array) – 3d array defining the gradient mask for controller parameters K, has to have same shape as K

• B (array) – left 2d matrix defining the position of controller matrices with respect to closed loop, see (1)

• C (array) – right 2d matrix defining the position of controller matrices with respect to closed loop, see (1)

Returns:

• fval (float): spectral abscissa

• grad (array): jacobian, 1d array matching shape of x

Return type:

tuple containing

tdcpy.stabopt.gradients.func_gamma(x, E, P, hP, hK, Kmask, B, C)

function for of gamma(r) and its gradient TODO

The system (Closed Loop) is defined as

E dxdt(t) = SUM P[i] x(t-hP[i]) + SUM B * K[j] * C x(t-hK[j]) (1)

where K = x.reshape(Kshape) are controller parameters, hK controller delays, P and hP are controlled system matrices and delays, respectively.

Parameters:

• x (array) – vectorized controller parameters

• TODO (r)

• E (array) – 2d array defining LHS of Closed Loop

• P (array) – 3d array defining the Pi matrices (controlled system)

• hP (array) – 1d array defining delays associated with P

• hK (array) – 1d array defining delays associated with K (note that x := vec(K))

• Kmask (array) – 3d array defining the gradient mask for controller parameters K, has to have same shape as K

• B (array) – left 2d matrix defining the position of controller matrices with respect to closed loop, see (1)

• C (array) – right 2d matrix defining the position of controller matrices with respect to closed loop, see (1)

Returns:

• fval (float): TODO

• grad (array): TODO

Return type:

tuple containing

tdcpy.stabopt.gradients.func_cd(x, E, P, hP, Kmask, hK, B, C, uE, vE, **kwargs)

TODO

Assumptions:

1. cd is a function of x (and can be changed via changing p)

2. hP[0] == 0.0, hP are unique and sorted len > 1

Parameters:

• x (ndarray[tuple[Any, ...], dtype[_ScalarT]])

• E (ndarray[tuple[Any, ...], dtype[_ScalarT]])

• P (ndarray[tuple[Any, ...], dtype[_ScalarT]])

• hP (ndarray[tuple[Any, ...], dtype[_ScalarT]])

• Kmask (ndarray[tuple[Any, ...], dtype[_ScalarT]])

• hK (ndarray[tuple[Any, ...], dtype[_ScalarT]])

• B (ndarray[tuple[Any, ...], dtype[_ScalarT]])

• C (ndarray[tuple[Any, ...], dtype[_ScalarT]])

• uE (ndarray[tuple[Any, ...], dtype[_ScalarT]])

• vE (ndarray[tuple[Any, ...], dtype[_ScalarT]])

tdcpy.stabopt.gradients.grad_gamma0(x, E, P, hP, Kmask, hK, B, C)

gradient of gamma0 function for delay-difference equations

Computes the gradient of the function \(\gamma_0(p)\) of the delay-difference equation defined by matrices D and delays hD with respect to controller parameters.

Parameters:

• x (npt.NDArray) – optimization variables

• E (npt.NDArray) – matrix defining LHS of closed loop

• P (npt.NDArray) – system matrix of open-loop system

• hP (npt.NDArray) – system delays of open-loop system

• Kmask (npt.NDArray) – mask for the controller entries

• hK (npt.NDArray) – controller delays

• B (npt.NDArray) – input matrix of closed-loop system

• C (npt.NDArray) – output matrix of closed-loop system

Returns:

• fvalfloat

  spectral abscissa

• gradarray

  gradient, 1d array matching shape of x

Return type:

A tuple containing

Notes

The function computes the gradient of the function \(\gamma_0\) of the delay-difference equation defined by matrices D and delays hD with respect to controller parameters. For the normalized DDE defined as .. math:

0 = I x(t) + H_1 x(t-h_1) + ... + H_m x(t-h_m),

with H_i = linalg.solve(D[:, :, 0], D[:, :, i]), hH_i = hD[i].

• assumes that the controller parameters affect the DDE

Examples

>>> import numpy as np
>>> from tdcpy.stabopt.gradients import grad_gamma0, gradient_test
>>> np.random.seed(0)
>>> E = np.array([[1., 0., 0.], [0., 1., 0.], [0., 0., 0.]])
>>> P = np.random.rand(3,3,2)
>>> hP = np.array([0., 1.])
>>> hK = np.array([0., 1.])
>>> Kmask = np.ones((1,3,2), dtype=bool)
>>> B = np.random.rand(3,1)
>>> C = np.random.rand(3,3)
>>> fval, grad = grad_gamma0(x=np.random.rand(6), E=E, P=P, hP=hP, Kmask=Kmask, hK=hK, B=B, C=C)
>>> print(f"gamma0: {fval}, grad: {grad}")
gamma0: 95.72387379041314, grad: [37.77946154 37.77946154  6.93863878  6.93863878 20.07062041 20.07062041]
>>> # test gradient
>>> grad_analytical, grad_numerical = gradient_test(grad_gamma0,x=np.random.rand(6),args=(E, P, hP, Kmask, hK, B, C))
Gradient test passed norm=5.389855703635268e-11 < 1e-06

tdcpy.stabopt.gradients.gradient_test(func, x, args, h=0.0001, tolerance=1e-06)

function to test the numerical accuracy of the computed gradient using central differences The method uses central differences for testing the numerical gradient.

Args: TODO

func: cost function h: step size E, P, hP, hK, Kmask, B, C: arguments

Returns:

tuple containing:

• fgradnpt.NDArray

  analytical gradient

• fgrad_numnpt.NDArray

  numerical gradient

Notes

The function computes the numerical gradient using central differences and compares it to the analytical gradient provided by the function func. If the norm of the difference between the two gradients is less than the specified tolerance, the test is considered passed.

Examples

>>> import numpy as np
>>> from tdcpy.stabopt.gradients import gradient_test
>>> def func_example(x, E, P, hP, hK, Kmask, B, C):
...     # example function returning cost and gradient
...     cost = np.sum(x**2)  # simple quadratic cost
...     grad = 2*x           # gradient of the cost
...     return cost, grad
>>> x0 = np.array([1.0, 2.0, 3.0])
>>> E = np.eye(2)
>>> P = np.random.rand(2, 2, 2)
>>> hP = np.array([0.1, 0.2])
>>> hK = np.array([0.1, 0.2])
>>> Kmask = np.ones((2, 2, 2), dtype=bool)
>>> B = np.random.rand(2, 2)
>>> C = np.random.rand(2, 2)
>>> fgrad, fgrad_num = gradient_test(func_example, x0, (E, P, hP, hK, Kmask, B, C))
Gradient test passed norm=0.0 < 1e-06

Parameters:

• func (Callable)

• x (ndarray[tuple[Any, ...], dtype[_ScalarT]])

• args (tuple)

• h (float)

• tolerance (float)

tdcpy.stabopt.utils module

Stabopt utils

tdcpy.stabopt.utils.diff_dependency_mask(Kmask, uE, vE, B, C, **kwargs)

Creates dependency mask for coefficients of delay difference eqations on controller parameters

Parameters:

• Kmask (npt.NDArray) – 3d boolean array of shape (p, q, hK) where p is the number of inputs, q the number of outputs and hK the number of controller delays. Kmask[:,:,i] is the mask for the i-th controller term.

• uE (npt.NDArray) – left null space of E matrix of the delay difference equation

• vE (npt.NDArray) – right null space of E matrix of the delay difference equation

• B (npt.NDArray) – input matrix of the delay difference equation

• C (npt.NDArray) – output matrix of the delay difference equation

• **kwargs (additional arguments, currently not used)

Returns:

3d boolean array of shape (p, q, hK) where p is the number of inputs, q the number of outputs and hK the number of controller delays. Kmask[:,:,i] is the mask for the i-th controller term, indicating which coefficients of the delay difference equation depend on the controller parameters.

Return type:

npt.NDArray

Notes

The dependency mask is computed by evaluating the expression

uE.T @ B @ K[:,:,i] @ C @ vE

for each controller term i, where K[:,:,i] is the mask for the i-th controller term. If the result is non-zero, the corresponding coefficients of the delay difference equation depend on the controller parameters.

Examples

>>> import numpy as np
>>> from tdcpy.stabopt.utils import diff_dependency_mask
>>> Kmask = np.array([[[1, 0], [0, 1]], [[0, 1], [1, 0]]], dtype=bool)
>>> uE = np.array([[1], [0]])
>>> vE = np.array([[1], [0]])
>>> B = np.array([[1, 0], [0, 1]])
>>> C = np.array([[1, 0], [0, 1]])
>>> diff_dependency_mask(Kmask, uE, vE, B, C)
array([[[ True, False],
        [False,  True]]])

tdcpy.stabopt.utils.check_diff(E, B, C, **kwargs)

Checks if delay difference equation is dependent on controller parameters

Parameters:

• TODO

• **kwargs

• E (ndarray[tuple[Any, ...], dtype[_ScalarT]])

• B (ndarray[tuple[Any, ...], dtype[_ScalarT]])

• C (ndarray[tuple[Any, ...], dtype[_ScalarT]])

Returns:

TODO

Return type:

bool

Module contents