Basic conversions of descriptor system models

  gprescale(sys; withB = true, withC = true) -> (sysbal, D1, D2)

Balance a descriptor system sys = (A-λE,B,C,D) by reducing the 1-norm of the matrix

         S =  ( abs(A)+abs(E)  abs(B) )
              (    abs(C)        0    )

by row and column balancing using diagonal matrices D1 and D2 to make the rows and columns of

              diag(D1,I)  * S * diag(D2,I)

as close in norm as possible.

The balancing can be performed optionally on the following particular matrices:

    S = abs(A)+abs(E),             if withB = false and withC = false ,
    S = ( abs(A)+abs(E)  abs(B) ), if withC = false,    
    S = ( abs(A)+abs(E) ),         if withB = false .
        (   abs(C)     )

The diagonal elements of the resulting D1 and D2 are the nearest integer powers of 2 resulting from the optimal diagonal scaling determined from a modified version of the algorithm of [1]. For a standard system with E = I, D1 = inv(D2) and D2 is computed using an extension of the scaling approach implemented in the LAPACK family *GEBAL.f.

This function is merely an interface to the functions lsbalance! of the MatrixPencils package.

[1] F.M.Dopico, M.C.Quintana and P. van Dooren, "Diagonal scalings for the eigenstructure of arbitrary pencils", SIMAX, 43:1213-1237, 2022.

c2d(sysc, Ts, meth = "zoh"; x0, u0, standard = true, fast = true, prewarp_freq = freq, 
           state_mapping = false, simple_infeigs = true, 
           atol = 0, atol1 = atol, atol2 = atol, rtol = n*ϵ) -> (sysd, xd0, Mx, Mu)

Compute for the continuous-time descriptor system sysc = (A-sE,B,C,D) with the proper transfer function matrix Gc(λ) and for a sampling time Ts, the corresponding discretized descriptor system sysd = (Ad-zEd,Bd,Cd,Dd) with the transfer function matrix Gd(z) according to the selected discretization method specified by meth. The keyword argument standard specifies the option to compute a standard state-space realization of sysd (i.e., with Ed = I), if standard = true (default), or a descriptor system realization if standard = false. The keyword argument simple_infeigs = true indicates that only simple infinite eigenvalues of the pair (A,E) are to be expected (default). The setting simple_infeigs = false indicates that possible uncontrollable or unobservable higher order infinite generalized eigenvalues of the pair (A,E) are present and have to be removed. xd0 is the mapped initial condition of the state of the discrete-time system sysd determined from the initial conditions of the state x0 and input u0 of the continuous-time system sysc. The keyword argument state_mapping = true specifies the option to compute the state mapping matrices Mx and Mu such that the values xc(t) and xd(t) of the system state vectors of the continuous-time system sysc and of the discrete-time system sysd, respectively, and of the input vector u(t) are related as xc(t) = Mx*xd(t)+Mu*u(t). If state_mapping = false (the default option), then Mx = nothing and Mu = nothing.

The following discretization methods can be performed by appropriately selecting meth:

"zoh"     - zero-order hold on the inputs (default); 

"foh"     - linear interpolation of inputs (also known as first-order hold);

"impulse" - impulse-invariant discretization; 

"Tustin"  - Tustin transformation (also known as trapezoidal integration): a nonzero prewarping frequency
            `freq` can be specified using the keyword parameter `prewarp_freq = freq` to ensure 
            `Gd(exp(im*freq*Ts)) = Gc(im*freq)`.

The underlying pencil manipulation algorithms employ rank determinations based on either the use of rank revealing QR-decomposition with column pivoting, if fast = true (default), or the SVD-decomposition, if fast = false. The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

The keyword arguments atol1, atol2 and rtol specify, respectively, the absolute tolerance for the nonzero elements of A, B, C, D, the absolute tolerance for the nonzero elements of E, and the relative tolerance for the nonzero elements of A, B, C, D and E. The default relative tolerance is n*ϵ, where n is the order of the square matrices A and E, and ϵ is the working machine epsilon. The keyword argument atol can be used to simultaneously set atol1 = atol and atol2 = atol.

rd = c2d(rc, Ts, meth = "zoh"; prewarp_freq = freq, atol = 0, rtol = n*ϵ)

Compute for the continuous-time rational transfer function rc and for a sampling time Ts, the corresponding discretized rational transfer function rd according to the selected discretization method specified by meth.

The following discretization methods can be performed by appropriately selecting meth:

"zoh"     - zero-order hold on the inputs (default); 

"foh"     - linear interpolation of inputs (also known as first-order hold);

"impulse" - impulse-invariant discretization; 

"matched" - matched pole-zero method (see [1]); 

"Tustin"  - Tustin transformation (also known as trapezoidal integration): a nonzero prewarping frequency
            `freq` can be specified using the keyword parameter `prewarp_freq = freq` to ensure 
            `rd(exp(im*freq*Ts)) = rc(im*freq)`.

The keyword arguments atol and rtol specify, respectively, the absolute and relative tolerances, respectively, for the nonzero coefficients of the numerator and denominator polynomials of rc. The default relative tolerance is n*ϵ, where n is the order of rc (i.e., the maximum degree of the numerator and denominator polynomials) and ϵ is the working machine epsilon.

References:

[1] G.F. Franklin, D.J. Powell and M.L. Workman, Digital Control of Dynamic Systems (3rd Edition), Prentice Hall, 1997.

Rd = c2d(Rc, Ts, meth = "zoh"; prewarp_freq = freq, atol = 0, rtol = n*ϵ)

Compute for the continuous-time rational transfer function matrix Rc and for a sampling time Ts, the corresponding discretized rational transfer function matrix Rd according to the selected discretization method specified by meth.

The following discretization methods can be performed by appropriately selecting meth:

  "zoh"     - zero-order hold on the inputs (default); 

  "foh"     - linear interpolation of inputs (also known as first-order hold);

  "impulse" - impulse-invariant discretization; 
  
  "matched" - matched pole-zero method (see [1]); 
  
  "Tustin"  - Tustin transformation (also known as trapezoidal integration): a nonzero prewarping frequency
              `freq` can be specified using the keyword parameter `prewarp_freq = freq` to ensure 
              `Rd(exp(im*freq*Ts)) = Rc(im*freq)`.

The keyword arguments atol and rtol specify, respectively, the absolute and relative tolerances, respectively, for the nonzero coefficients of the numerator and denominator polynomials of the elements of Rc. The default relative tolerance is 10*ϵ, where ϵ is the working machine epsilon.

References:

[1] G.F. Franklin, D.J. Powell and M.L. Workman, Digital Control of Dynamic Systems (3rd Edition), Prentice Hall, 1997.

R = dss2rm(sys; fast = true, atol = 0, atol1 = atol, atol2 = atol, gaintol = atol, rtol = min(atol1,atol2) > 0 ? 0 : n*ϵ, val)

Build for the descriptor system sys = (A-λE,B,C,D) the rational matrix R(λ) = C*inv(λE-A)*B+D representing the transfer function matrix of the system sys.

The keyword arguments atol1 and atol2 specify the absolute tolerances for the elements of A, B, C, D, and, respectively, of E, and rtol specifies the relative tolerances for the nonzero elements of A, B, C, D and E. The default relative tolerance is n*ϵ, where n is the maximal dimension of state, input and output vectors, and ϵ is the working machine epsilon. The keyword argument atol can be used to simultaneously set atol1 = atol, atol2 = atol and gaintol = atol.

The keyword argument gaintol specifies the threshold for the magnitude of the nonzero elements of the gain matrix C*inv(γE-A)*B+D, where γ = val if val is a number or γ is a randomly chosen complex value of unit magnitude, if val = missing. Generally, val should not be a zero of any of entries of R.

Method: Each rational entry of R(λ) is constructed from its numerator and denominator polynomials corresponding to its finite zeros, finite poles and gain using the method of [1].

References:

[1] A. Varga Computation of transfer function matrices of generalized state-space models. Int. J. Control, 50:2543–2561, 1989.

P = dss2pm(sys; fast = true, atol = 0, atol1 = atol, atol2 = atol, gaintol = 0, rtol = min(atol1,atol2) > 0 ? 0 : n*ϵ, val)

Build for the descriptor system sys = (A-λE,B,C,D) the polynomial matrix P(λ) = C*inv(λE-A)*B+D representing the transfer function matrix of the system sys.

The keyword arguments atol1 and atol2 specify the absolute tolerances for the elements of A, B, C, D, and, respectively, of E, and rtol specifies the relative tolerances for the nonzero elements of A, B, C, D and E. The default relative tolerance is n*ϵ, where n is the maximal dimension of state, input and output vectors, and ϵ is the working machine epsilon. The keyword argument atol can be used to simultaneously set atol1 = atol, atol2 = atol and gaintol = atol.

The keyword argument gaintol specifies the threshold for the magnitude of the nonzero elements of the gain matrix C*inv(γE-A)*B+D, where γ = val if val is a number or γ is a randomly chosen complex value of unit magnitude, if val = missing. Generally, val should not be a zero of any of entries of P.

Method: Each entry of P(λ) is constructed from the polynomial corresponding to its finite zeros and gain using the method of [1].

References:

[1] A. Varga Computation of transfer function matrices of generalized state-space models. Int. J. Control, 50:2543–2561, 1989.

gbilin(sys, g; compact = true, minimal = false, standard = true, atol = 0, atol1 = atol, atol2 = atol, rtol = n*ϵ) -> (syst, ginv)

Compute for the descriptor system sys = (A-λE,B,C,D) with the transfer function matrix G(λ) and a first degree real rational transfer function g = g(δ), the descriptor system realization syst = (At-δEt,Bt,Ct,Dt) of G(g(δ)) corresponding to the bilinear transformation λ = g(δ) = (aδ+b)/(cδ+d). For a continuous-time transfer function g(δ), δ = s, the complex variable in the Laplace transform, while for a discrete-time transfer function, δ = z, the complex variable in the Z-transform. syst inherits the sampling-time of sys1. sysi1 is the transfer function ginv(λ) = (d*λ-b)/(-c*λ+a) representing the inverse of the bilinear transformation g(δ) (i.e., g(ginv(λ)) = 1).

The keyword argument compact can be used to specify the option to compute a compact descriptor realization without non-dynamic modes, if compact = true (the default option) or to disable the ellimination of non-dynamic modes if compact = false.

The keyword argument minimal specifies the option to compute minimal descriptor realization, if minimal = true, or a nonminimal realization if minimal = false (the default option).

The keyword argument standard specifies the option to compute a standard state-space (if possible) realizations of syst, if standard = true (default), or a descriptor system realization if standard = false.

The keyword arguments atol1, atol2 and rtol specify, respectively, the absolute tolerance for the nonzero elements of A, B, C, D, the absolute tolerance for the nonzero elements of E, and the relative tolerance for the nonzero elements of A, B, C, D and E. The default relative tolerance is n*ϵ, where n is the order of the square matrices A and E, and ϵ is the working machine epsilon. The keyword argument atol can be used to simultaneously set atol1 = atol and atol2 = atol.