Analysis of periodic systems

pspole(psys::PeriodicStateSpace{PM}[,K]; kwargs...) -> val

Return for the periodic system psys = (A(t),B(t),C(t),D(t)) the complex vector val containing the characteristic exponents of the periodic matrix A(t) (also called the poles of the system psys).

Depending on the underlying periodic matrix type PM, the optional argument K and keyword arguments kwargs may have the following values:

• if PM = PeriodicFunctionMatrix, or PM = PeriodicSymbolicMatrix, or PM = PeriodicTimeSeriesMatrix, then K is the number of factors used to express the monodromy matrix of A(t) (default: K = 1) and kwargs are the keyword arguments of PeriodicMatrices.pseig(::PeriodicFunctionMatrix);

• if PM = HarmonicArray, then K is the number of harmonic components used to represent the Fourier series of A(t) (default: K = max(10,nh-1), nh is the number of harmonics terms of A(t)) and kwargs are the keyword arguments of PeriodicMatrices.psceig(::HarmonicArray);

• if PM = FourierFunctionMatrix, then K is the number of harmonic components used to represent the Fourier series of A(t) (default: K = max(10,nh-1), nh is the number of harmonics terms of A(t)) and kwargs are the keyword arguments of PeriodicMatrices.psceig(::FourierFunctionMatrix);

• if PM = PeriodicMatrix or PM = PeriodicArray, then K is the starting sample time (default: K = 1) and kwargs are the keyword arguments of PeriodicMatrices.psceig(::PeriodicMatrix);

pszero(psys::PeriodicStateSpace{HarmonicArray}[, N]; P, atol, rtol, fast) -> val
pszero(psys::PeriodicStateSpace{PeriodicFunctionMatrix}[, N]; P, atol, rtol, fast) -> val
pszero(psys::PeriodicStateSpace{PeriodicSymbolicMatrix}[, N]; P, atol, rtol, fast) -> val

Compute the finite and infinite zeros of a continuous-time periodic system psys = (A(t), B(t), C(t), D(t)) in val, where the periodic system matrices A(t), B(t), C(t), and D(t) are in a Harmonic Array, or Periodic Function Matrix, or Periodic Symbolic Matrix representation (the last two representation are automatically converted to a Harmonic Array representation). N is the number of selected harmonic components in the Fourier series of the system matrices (default: N = max(20,nh-1), where nh is the maximum number of harmonics terms) and the keyword parameter P is the number of full periods to be considered (default: P = 1) to build a frequency-lifted LTI representation based on truncated block Toeplitz matrices.

The computation of the zeros of the complex lifted system is performed by reducing the corresponding system pencil to an appropriate Kronecker-like form which exhibits the finite and infinite eigenvalues. The reduction is performed using orthonal similarity transformations and involves rank decisions based on rank revealing QR-decompositions with column pivoting, if fast = true, or, the more reliable, SVD-decompositions, if fast = false. For a system psys of period T, the finite zeros are determined as those eigenvalues which have imaginary parts in the interval [-ω/2, ω/2], where ω = 2π/(P*T). To eliminate possible spurious finite eigenvalues, the intersection of two finite eigenvalue sets is computed for two lifted systems obtained for N and N+2 harmonic components. The infinite zeros are determined as the infinite zeros of the LTI system (A(ti), B(ti), C(ti), D(ti)) resulting for a random time value ti. Warning: While this evaluation of the number of infinite zeros mostly provides the correct result, there is no theoretical assessment of this approach (counterexamples are welcome!).

The keyword arguments atol and rtol specify the absolute and relative tolerances for the nonzero elements of the underlying lifted system pencil, respectively. The default relative tolerance is n*ϵ, where n is the size of the smallest dimension of the pencil, and ϵ is the working machine epsilon.

 pszero(psys::PeriodicStateSpace{PeriodicMatrix}[, K]; atol, rtol, fast) -> val
 pszero(psys::PeriodicStateSpace{PeriodicArray}[, K]; atol, rtol, fast) -> val

Compute the finite and infinite zeros of a discrete-time periodic system psys = (A(t), B(t), C(t), D(t)) in val, where the periodic system matrices A(t), B(t), C(t), and D(t) are in either Periodic Matrix or Periodic Array representation. The optional argument K specifies a desired time to start the sequence of periodic matrices (default: K = 1).

The computation of zeros relies on the fast structure exploiting reduction of singular periodic pencils as described in [1], which separates a matrix pencil M-λN which contains the infinite and finite zeros and the left and right Kronecker indices. The reduction is performed using orthonal similarity transformations and involves rank decisions based on rank revealing QR-decompositions with column pivoting, if fast = true, or, the more reliable, SVD-decompositions, if fast = false.

The keyword arguments atol and rtol specify the absolute and relative tolerances for the nonzero elements of the periodic system matrices, respectively. The default relative tolerance is n*ϵ, where n is proportional with the largest dimension of the periodic matrices, and ϵ is the working machine epsilon.

References

[1] A. Varga and P. Van Dooren. Computing the zeros of periodic descriptor systems. Systems & Control Letters 50:371–381, 2003.

 isstable(psys[, K = 1]; smarg = 1, fast = false, offset = sqrt(ϵ), kwargs...) -> Bool

Return true if the continuous-time periodic system psys has only stable characteristic multipliers and false otherwise.

To assess the stability, the absolute values of the characteristic multipliers (i.e., the eigenvalues of the monodromy matrix) must be less than smarg-β, where smarg is the discrete-time stability margin (default: smarg = 1) and β is an offset specified via the keyword parameter offset = β to be used to numerically assess the stability of eigenvalues. The default value used for β is sqrt(ϵ), where ϵ is the working machine precision.

The monodromy matrix is determined as a product K state transition matrices (default: K = 1) computed by integrating numerically a homogeneous linear ODE with periodic coefficients. If fast = false (default) then the characteristic multipliers are computed using an approach based on the periodic Schur decomposition [1], while if fast = true the structure exploiting reduction [2] of an appropriate lifted pencil is employed. This later option may occasionally lead to inaccurate results for large number of matrices.

References

[1] A. Bojanczyk, G. Golub, and P. Van Dooren, The periodic Schur decomposition. Algorithms and applications, Proc. SPIE 1996.

[2] A. Varga & P. Van Dooren. Computing the zeros of periodic descriptor systems. Systems and Control Letters, 50:371-381, 2003.

isstable(psys; smarg = 1, fast = false, offset = sqrt(ϵ)) -> Bool

Return true if the discrete-time periodic system psys has only stable characteristic multipliers and false otherwise.

To assess the stability, the absolute values of the characteristic multipliers (i.e., the eigenvalues of the monodromy matrix) must be less than smarg-β, where smarg is the discrete-time stability margin (default: smarg = 1) and β is an offset specified via the keyword parameter offset = β to be used to numerically assess the stability of eigenvalues. The default value used for β is sqrt(ϵ), where ϵ is the working machine precision.

If fast = false (default) then the characteristic multipliers are computed using an approach based on the periodic Schur decomposition [1], while if fast = true the structure exploiting reduction [2] of an appropriate lifted pencil is employed. This later option may occasionally lead to inaccurate results for large number of matrices.

References

[1] A. Bojanczyk, G. Golub, and P. Van Dooren, The periodic Schur decomposition. Algorithms and applications, Proc. SPIE 1996.

[2] A. Varga & P. Van Dooren. Computing the zeros of periodic descriptor systems. Systems and Control Letters, 50:371-381, 2003.

pshanorm(psys; smarg = 1, lifting = false, offset = sqrt(ϵ)) -> nrm

Compute the Hankel-norm of a stable discrete-time periodic system psys = (A(t),B(t),C(t),D(t)). For the computation of the norm, the formulas given in [1] are employed. The Hankel norm is computed as

 nrm = maximum(sqrt(Λ(P(t)Q(t))),

where P(t) is the controllability Gramian satisfying the periodic Lyapunov equation

 P(t+1) = A(t)P(t)A(t)' + B(t)B(t)',

and Q(t) is the observability Gramian satisfying the periodic Lyapunov equation

Q(t) = A(t)'Q(t+1)A(t) + C(t)'C(t) .

To assess the stability, the absolute values of the characteristic multipliers of A(t) must be less than smarg-β, where smarg is the discrete-time stability margin (default: smarg = 1) and β is an offset specified via the keyword parameter offset = β to be used to numerically assess the stability of eigenvalues. The default value used for β is sqrt(ϵ), where ϵ is the working machine precision.

If lifting = false (default), then the norm is evaluated using an approach based on the periodic Schur decomposition of A(t), while if lifting = true the norm of the lifted cyclic system is evaluated. This later option may lead to excessive computational times for large matrices or large periods.

References

[1] S. Bittanti and P. Colaneri. Periodic Systems : Filtering and Control. Springer-Verlag London, 2009.

pshanorm(psys, K; smarg = 1, offset = sqrt(ϵ), solver = "auto", reltol = 1.e-4, abstol = 1.e-7) -> nrm

Compute the Hankel-norm of a stable continuous-time periodic system psys = (A(t),B(t),C(t),D(t)). For the computation of the norm, the approach suggested in [1] is employed, in conjunction with the multiple-shooting approach using K discretization points. For a periodic system of period T, the Hankel-norm is defined as

 nrm = sqrt(max(Λ(P(t)Q(t)))), for t ∈ [0,T]

where P(t) is the controllability Gramian satisfying the periodic differential Lyapunov equation

 .
 P(t) = A(t)P(t)A(t)' + B(t)B(t)',

and Q(t) is the observability Gramian satisfying the periodic differential Lyapunov equation

 .
-Q(t) = A(t)'Q(t)A(t) + C(t)'C(t) .

The norm is evaluated from the K time values of P(t) and Q(t) in the interval [0,T] and the precision is (usually) better for larger values of K.

To assess the stability, the absolute values of the characteristic multipliers of A(t) must be less than smarg-β, where smarg is the discrete-time stability margin (default: smarg = 1) and β is an offset specified via the keyword parameter offset = β to be used to numerically assess the stability of eigenvalues. The default value used for β is sqrt(ϵ), where ϵ is the working machine precision.

The ODE solver to be employed to convert the continuous-time problem into a discrete-time problem can be specified using the keyword argument solver, together with the required relative accuracy reltol (default: reltol = 1.e-4) and absolute accuracy abstol (default: abstol = 1.e-7). Depending on the desired relative accuracy reltol, lower order solvers are employed for reltol >= 1.e-4, which are generally very efficient, but less accurate. If reltol < 1.e-4, higher order solvers are employed able to cope with high accuracy demands.

The following solvers from the OrdinaryDiffEq.jl package can be selected:

solver = "non-stiff" - use a solver for non-stiff problems (Tsit5() or Vern9());

solver = "stiff" - use a solver for stiff problems (Rodas4() or KenCarp58());

solver = "auto" - use the default solver, which automatically detects stiff problems (AutoTsit5(Rosenbrock23()) or AutoVern9(Rodas5())).

References

[1] A. Varga, On solving periodic differential matrix equations with applications to periodic system norms computation. Proc. CDC/ECC, Seville, p.6545-6550, 2005.

psh2norm(psys; adj = false, smarg = 1, fast = false, offset = sqrt(ϵ)) -> nrm

Compute the H2-norm of a discrete-time periodic system psys = (A(t),B(t),C(t),D(t)). For the computation of the norm, the formulas given in [1] are employed. For adj = false the norm is computed as

 nrm = sqrt(sum(tr(C(t)P(t)C(t)'+D(t)*D(t)'))),

where P(t) is the controllability Gramian satisfying the periodic Lyapunov equation

 P(t+1) = A(t)P(t)A(t)' + B(t)B(t)',

while for adj = true the norm is computed as

nrm = sqrt(sum(tr(B(t)'Q(t+1)B(t)+D(t)'*D(t)))),

where Q(t) is the observability Gramian satisfying the periodic Lyapunov equation

Q(t) = A(t)'Q(t+1)A(t) + C(t)'C(t) .

The norm is set to infinity for an unstable system.

To assess the stability, the absolute values of the characteristic multipliers of A(t) must be less than smarg-β, where smarg is the discrete-time stability margin (default: smarg = 1) and β is an offset specified via the keyword parameter offset = β to be used to numerically assess the stability of eigenvalues. The default value used for β is sqrt(ϵ), where ϵ is the working machine precision.

If fast = false (default) then the norm is evaluated using an approach based on the periodic Schur decomposition of A(t), while if fast = true the norm of the lifted standard system is evaluated. This later option may occasionally lead to inaccurate results for large number of matrices.

References

[1] S. Bittanti and P. Colaneri. Periodic Systems : Filtering and Control. Springer-Verlag London, 2009.

psh2norm(psys, K; adj = false, smarg = 1, fast = false, offset = sqrt(ϵ), solver = "auto", reltol = 1.e-4, abstol = 1.e-7, quad = false) -> nrm

Compute the H2-norm of a continuous-time periodic system psys = (A(t),B(t),C(t),D(t)). For the computation of the norm, the formulas given in [1] are employed, in conjunction with the multiple-shooting approach of [2] using K discretization points. For a periodic system of period T, for adj = false the norm is computed as

 nrm = sqrt(Integral(tr(C(t)P(t)C(t)')))dt/T),

where P(t) is the controllability Gramian satisfying the periodic differential Lyapunov equation

 .
 P(t) = A(t)P(t)A(t)' + B(t)B(t)',

while for adj = true the norm is computed as

nrm = sqrt(Integral(tr(B(t)'Q(t)B(t)))dt/T),

where Q(t) is the observability Gramian satisfying the periodic differential Lyapunov equation

 .
-Q(t) = A(t)'Q(t)A(t) + C(t)'C(t) .

If quad = true, a simple quadrature formula based on the sum of time values is employed (see [2]). This option ensures a faster evaluation, but is potentially less accurate then the exact evaluation employed if quad = false (default).

The norm is set to infinity for an unstable system or for a non-zero D(t).

To assess the stability, the absolute values of the characteristic multipliers of A(t) must be less than smarg-β, where smarg is the discrete-time stability margin (default: smarg = 1) and β is an offset specified via the keyword parameter offset = β to be used to numerically assess the stability of eigenvalues. The default value used for β is sqrt(ϵ), where ϵ is the working machine precision. If fast = false (default) then the stability is checked using an approach based on the periodic Schur decomposition of A(t), while if fast = true the stability is checked using a lifting-based approach.

The ODE solver to be employed to convert the continuous-time problem into a discrete-time problem can be specified using the keyword argument solver, together with the required relative accuracy reltol (default: reltol = 1.e-4) and absolute accuracy abstol (default: abstol = 1.e-7). Depending on the desired relative accuracy reltol, lower order solvers are employed for reltol >= 1.e-4, which are generally very efficient, but less accurate. If reltol < 1.e-4, higher order solvers are employed able to cope with high accuracy demands.

The following solvers from the OrdinaryDiffEq.jl package can be selected:

solver = "non-stiff" - use a solver for non-stiff problems (Tsit5() or Vern9());

solver = "stiff" - use a solver for stiff problems (Rodas4() or KenCarp58());

solver = "auto" - use the default solver, which automatically detects stiff problems (AutoTsit5(Rosenbrock23()) or AutoVern9(Rodas5())).

References

[1] P. Colaneri. Continuous-time periodic systems in H2 and H∞: Part I: Theoretical Aspects. Kybernetika, 36:211-242, 2000.

[2] A. Varga, On solving periodic differential matrix equations with applications to periodic system norms computation. Proc. CDC/ECC, Seville, p.6545-6550, 2005.

pslinfnorm(psys, hinfnorm = false, rtolinf = 0.001, fast = true, offset = sqrt(ϵ)) -> (linfnorm, fpeak)
pshinfnorm(psys, rtolinf = 0.001, fast = true, offset = sqrt(ϵ)) -> (hinfnorm, fpeak)

Compute for a discrete-time periodic system psys = (A(t),B(t),C(t),D(t) with the lifted transfer function matrix G(λ) the L∞ norm linfnorm with pslinfnorm or the H∞ norm hinfnorm with pshinfnorm (i.e., the peak gain of G(λ)) and the corresponding peak frequency fpeak, where the peak gain is achieved. If hinfnorm = true, the H∞ norm is computed with pslinfnorm.

The L∞ norm is infinite if psys has poles on the stability domain boundary, i.e., on the unit circle. The H∞ norm is infinite if psys has unstable poles.

To check the lack of poles on the stability domain boundary, the poles of psys must not have moduli within the interval [1-β,1+β], where β is the stability domain boundary offset. The offset β to be used can be specified via the keyword parameter offset = β. The default value used for β is sqrt(ϵ), where ϵ is the working machine precision.

The keyword argument rtolinf specifies the relative accuracy for the computed infinity norm. The default value used for rtolinf is 0.001.

The computation of the L∞ norm is based on the algorithm proposed in [1]. The involved computations of characteristic multipliers are performed either with the fast reduction method of [2], if fast = true or if time-varying dimensions are present, or the generalized periodic Schur decomposition based method of [3], if fast = false.

References

[1] A. Varga. "Computation of L∞-norm of linear discrete-time periodic systems." Proc. MTNS, Kyoto, 2007.

[2] A. Varga and P. Van Dooren. Computing the zeros of periodic descriptor systems. Systems & Control Letters 50:371–381, 2003.

[3] Kressner, D. An efficient and reliable implementation of the periodic QZ algorithm. IFAC Workshop on Periodic Control Systems (PSYCO 2001), Como (Italy), August 27-28 2001. Periodic Control Systems 2001 (IFAC Proceedings Volumes), Pergamon.

pslinfnorm(psys, K = 100; hinfnorm = false, rtolinf = 0.001, offset = sqrt(ϵ), reltol, abstol, dt) -> (linfnorm, fpeak)
pshinfnorm(psys, K = 100; rtolinf = 0.001, offset = sqrt(ϵ), reltol, abstol, dt) -> (linfnorm, fpeak)

Compute for a continuous-time periodic system psys = (A(t),B(t),C(t),D(t) the L∞ norm linfnorm with pslinfnorm or the H∞ norm hinfnorm with pshinfnorm as defined in [1]. If hinfnorm = true, the H∞ norm is computed with pslinfnorm. The corresponding peak frequency fpeak, where the peak gain is achieved, is usually not determined, excepting in some limiting cases. The L∞ norm is infinite if psys has poles on the imaginary axis.

To check the lack of poles on the imaginary axis, the characteristic exponents of A(t) must not have real parts in the interval [-β,β], where β is the stability domain boundary offset. The offset β to be used can be specified via the keyword parameter offset = β. The default value used for β is sqrt(ϵ), where ϵ is the working machine precision.

A bisection based algorith, as described in [2], is employed to approximate the L∞ norm, and the keyword argument rtolinf specifies the relative accuracy for the computed infinity norm. The default value used for rtolinf is 0.001.

If hinfnorm = true, the H∞ norm is computed. In this case, the stability of the system is additionally checked and the H∞ norm is infinite for an unstable system. To check the stability, the characteristic exponents of A(t) must have real parts less than -β.

The ODE solver to be employed to compute the characteristic multipliers of the system Hamiltonian can be specified using the keyword argument solver (default: solver = "symplectic") together with the required relative accuracy reltol (default: reltol = 1.e-4), absolute accuracy abstol (default: abstol = 1.e-7) and stepsize dt (default: dt = 0). The value stepsize is relevant only if solver = "symplectic", in which case an adaptive stepsize strategy is used ifdt = 0and a fixed stepsize is used ifdt > 0. Depending on the desired relative accuracyreltol, lower order solvers are employed forreltol >= 1.e-4, which are generally very efficient, but less accurate. Ifreltol < 1.e-4`, higher order solvers are employed able to cope with high accuracy demands.

The following solvers from the OrdinaryDiffEq.jl package can be selected:

solver = "non-stiff" - use a solver for non-stiff problems (Tsit5() or Vern9());

solver = "stiff" - use a solver for stiff problems (Rodas4() or KenCarp58());

solver = "linear" - use a special solver for linear ODEs (MagnusGL6()) with fixed time step dt;

solver = "symplectic" - use a symplectic Hamiltonian structure preserving solver (IRKGL16());

solver = "auto" - use the default solver, which automatically detects stiff problems (AutoTsit5(Rosenbrock23()) or AutoVern9(Rodas5())).

References:

[1] P. Colaneri. Continuous-time periodic systems in H2 and H∞: Part I: Theoretical Aspects. Kybernetika, 36:211-242, 2000.

[2] A. Varga, On solving periodic differential matrix equations with applications to periodic system norms computation. Proc. CDC/ECC, Seville, p.6545-6550, 2005.

pstimeresp(sys, u, t, x0 = zeros(sys.nx); state_history = false, solver, reltol, abstol) -> (y, tout, x)

Compute the time response of a periodic system sys = (A(t),B(t),C(t),D(t)) to the input signals described by u and t. The time vector t consists of regularly spaced time samples. The input u can be specified as a matrix with as many columns as the number of inputs of sys, in which case its i-th row specifies the input values at time t[i]. For a discrete-time system, u should be sampled at the same sampling rate Ts as sys and t must have all time steps equal to Ts or can be set to an empty vector. For continuous-time models, the input values are interpolated between samples using zero-order hold based interpolation. The vector x0 specifies the initial state vector at time t[1] and is set to zero when omitted.

The matrix y contains the resulting time history of the outputs of sys and the vector tout contains the corresponding values of the time samples. The i-th row of y contains the output values at time tout[i]. If the keyword parameter value state_history = true is used, then the matrix x contains the resulting time history of the state vector and its i-th row x[i,:] contains the state values at time tout[i]. The column dimension of the matrix x is n, the dimension of the state vector. For a discrete-time periodic system with time-varying state dimensions, n is the maximum value of the dimensions of the state vector over one period. The components of x[i,:] have trailing zero values if the corresponding state vector dimension is less than n. By default, the state history is not saved and x = nothing.

For a continuous-time model an equivalent discretized model is determined to be used for simulation. The discretization is performed by determining the monodromy matrix as a product of state transition matrices of the extended state-space matrix [A(t) B(t); 0 0] by integrating numerically the corresponding homogeneous linear ODE. The ODE solver to be employed can be specified using the keyword argument solver, together with the required relative accuracy reltol (default: reltol = 1.e-4), absolute accuracy abstol (default: abstol = 1.e-7) and/or the fixed step length dt (default: dt = Ts/10). Depending on the desired relative accuracy reltol, lower order solvers are employed for reltol >= 1.e-4, which are generally very efficient, but less accurate. If reltol < 1.e-4, higher order solvers are employed able to cope with high accuracy demands.

The following solvers from the OrdinaryDiffEq.jl package can be selected:

solver = "non-stiff" - use a solver for non-stiff problems (Tsit5() or Vern9());

solver = "stiff" - use a solver for stiff problems (Rodas4() or KenCarp58());

solver = "linear" - use a special solver for linear ODEs (MagnusGL6()) with fixed time step dt;

solver = "symplectic" - use a symplectic Hamiltonian structure preserving solver (IRKGL16());

solver = "auto" - use the default solver, which automatically detects stiff problems (AutoTsit5(Rosenbrock23()) or AutoVern9(Rodas5())).

For large numbers of product terms, parallel computation of factors can be alternatively performed by starting Julia with several execution threads. The number of execution threads is controlled either by using the -t/--threads command line argument or by using the JULIA_NUM_THREADS environment variable.

pstimeresp(sys, u, t, x0 = zeros(sys.nx); state_history = false, solver, reltol, abstol) -> (y, tout, x)

Compute the time response of a continuous-time periodic system sys = (A(t),B(t),C(t),D(t)) to the input signals described by u and t. The time vector t consists of regularly spaced time samples. The input u is specified as a vector of time dependent signals with as many components as the number of inputs of sys. The vector x0 specifies the initial state vector at time t[1] and is set to zero when omitted.

The matrix y contains the resulting time history of the outputs of sys and the vector tout contains the corresponding values of the time samples. The i-th row of y contains the output values at time tout[i]. If the keyword parameter value state_history = true is used, then the matrix x contains the resulting time history of the state vector and its i-th row contains the state values at time tout[i]. By default, the state history is not computed and x = nothing.

The ODE solver to be employed can be specified using the keyword argument solver (see below), together with the required relative accuracy reltol (default: reltol = 1.e-4), absolute accuracy abstol (default: abstol = 1.e-7) and the fixed step length dt (default: dt = 0), only used if solver = "symplectic". Depending on the desired relative accuracy reltol, lower order solvers are employed for reltol >= 1.e-4, which are generally very efficient, but less accurate. If reltol < 1.e-4, higher order solvers are employed able to cope with high accuracy demands.

The following solvers from the OrdinaryDiffEq.jl package can be selected:

solver = "non-stiff" - use a solver for non-stiff problems (Tsit5() or Vern9());

solver = "stiff" - use a solver for stiff problems (Rodas4() or KenCarp58());

solver = "linear" - use a special solver for linear ODEs (MagnusGL6()) with fixed time step dt;

solver = "symplectic" - use a symplectic Hamiltonian structure preserving solver (IRKGL16());

solver = "auto" - use the default solver, which automatically detects stiff problems (AutoTsit5(Rosenbrock23()) or AutoVern9(Rodas5())).

psstepresp(sys[, tfinal]; ustep = ones(psys.nu), x0 = zeros(psys.nx), timesteps = 100, 
         state_history = false, abstol, reltol) -> (y, tout, x)

Compute the time response of a periodic system sys = (A(t),B(t),C(t),D(t)) to step input signals. The final time tfinal, if not specified, is set equal to the period of the periodic system sys. For a discrete-time system, the final time tfinal must be commensurate with the system sampling time Ts (otherwise it is adjusted to the nearest smaller commensurate value). The keyword argument ustep is a vector with as many components as the inputs of sys and specifies the desired amplitudes of step inputs (default: all components are set to 1). The keyword argument x0 is a vector, which specifies the initial state vector at time 0, and is set to zero when omitted. The keyword argument timesteps specifies the number of desired simulation time steps (default: timesteps = 100).

If ns is the total number of simulation values, n the number of state components, p the number of system outputs and m the number of system inputs, then the resulting ns×p×m array y contains the resulting time histories of the outputs of sys, such that y[:,:,j] is the time response for the j-th input set to ustep[j] and the rest of inputs set to zero. The vector tout contains the corresponding values of the time samples. The i-th row y[i,:,j] contains the output values at time tout[i] of the j-th step response. If the keyword parameter value state_history = true is used, then the resulting ns×n×m arrayx contains the resulting time histories of the state vector and the i-th row x[i,:,j] contains the state values at time tout[i] of the j-th step response. For a discrete-time periodic system with time-varying state dimensions, n is the maximum value of the dimensions of the state vector over one period. The components of x[i,:,j] have trailing zero values if the corresponding state vector dimension is less than n. By default, the state history is not saved and x = nothing.

The total number of simulation values ns is set as follows: for a continuous-time system ns = timesteps+1 and for a discrete-time system ns = min(timesteps,tfinal/Ts)+1.

For a continuous-time model an equivalent discretized model is determined to be used for simulation, provided the time step Δ = tfinal/timesteps is commensurate with the system period T and tfinal >= T. The discretization is performed by determining the monodromy matrix as a product of state transition matrices of the extended state-space matrix [A(t) B(t); 0 0] by integrating numerically the corresponding homogeneous linear ODE. If the time step Δ is not commensurate with the period T or tfinal < T, then numerical integrations of the underlying ODE systems are performed. The ODE solver to be employed can be specified using the keyword argument solver, together with the required relative accuracy reltol (default: reltol = 1.e-4), absolute accuracy abstol (default: abstol = 1.e-7) and/or the fixed step length dt (default: dt = Ts/10). Depending on the desired relative accuracy reltol, lower order solvers are employed for reltol >= 1.e-4, which are generally very efficient, but less accurate. If reltol < 1.e-4, higher order solvers are employed able to cope with high accuracy demands.

The following solvers from the OrdinaryDiffEq.jl package can be selected:

solver = "non-stiff" - use a solver for non-stiff problems (Tsit5() or Vern9());

solver = "stiff" - use a solver for stiff problems (Rodas4() or KenCarp58());

solver = "linear" - use a special solver for linear ODEs (MagnusGL6()) with fixed time step dt;

solver = "symplectic" - use a symplectic Hamiltonian structure preserving solver (IRKGL16());

solver = "auto" - use the default solver, which automatically detects stiff problems (AutoTsit5(Rosenbrock23()) or AutoVern9(Rodas5())).

For large numbers of product terms, parallel computation of factors can be alternatively performed by starting Julia with several execution threads. The number of execution threads is controlled either by using the -t/--threads command line argument or by using the JULIA_NUM_THREADS environment variable.

pseval(psys, val) -> Gval

Evaluate for a finite λ = val, the value Gval of the transfer function matrix G(λ) of the lifted system of the discrete-time periodic system psys. val must not be a pole of psys.