Exported functions and types

Disk

Represents a perturbation disc in the complex plane. Disk(0.5, 2) represents all perturbations in the circle centered at 1.25 with radius 0.75, or in other words, a gain margin of 2 and a phase margin of 36.9 degrees.

A disk can be converted to a Nyquist exclusion disk by nyquist(disk) and plotted using plot(disk).

Arguments:

• γmin: Lower intercept
• γmax: Upper intercept
• c: Center
• r: Radius
• ϕm: Angle of tangent line through origin.

If γmax < γmin the disk is inverted. See diskmargin for disk margin computations.

Diskmargin

The notation follows "An Introduction to Disk Margins", Peter Seiler, Andrew Packard, and Pascal Gahinet

Fields:

α: The disk margin ω0: The worst-case frequency f0: The destabilizing perturbation f0 is a complex number with simultaneous gain and phase variation. This critical perturbation causes an instability with closed-loop pole on the imaginary axis at the critical frequency ω0 δ0: The uncertain element generating f0. γmin: The lower real-axis intercept of the disk (analogous to classical lower gain margin). γmax: The upper real-axis intercept of the disk (analogous to classical upper gain margin). ϕm: is analogous to the classical phase margin. σ: The skew parameter that was used to calculate the margin

Note, γmax and ϕm are in smaller than the classical gain and phase margins sicne the classical margins do not consider simultaneous perturbations in gain and phase.

The "disk" margin becomes a half plane for α = 2 and an inverted circle for α > 2. In this case, the upper gain margin is infinite. See the paper for more details, in particular figure 6.

ExtendedStateSpace{TE, T} <: AbstractStateSpace{TE}

A type that represents the two-input, two-output system

z  ┌─────┐  w
◄──┤     │◄──
   │  P  │
◄──┤     │◄──
y  └─────┘  u

where

• z denotes controlled outputs (sometimes called performance outputs)
• y denotes measured outputs
• w denotes external inputs, such as disturbances or references
• u denotes control inputs

The call lft(P, K) forms the (lower) linear fractional transform

z  ┌─────┐  w
◄──┤     │◄──
   │  P  │
┌──┤     │◄─┐
│y └─────┘ u│
│           │
│  ┌─────┐  │
│  │     │  │
└─►│  K  ├──┘
   │     │
   └─────┘

i.e., closing the lower loop around K.

An ExtendedStateSpace can be converted to a standard StateSpace by ss(P), this will keep all inputs and outputs, effectively removing the partitioning only.

When feedback is called on this type, defaults are automatically set for the feedback indices. Other functions defined for this type include

• system_mapping
• performance_mapping
• noise_mapping
• lft
• feedback has special overloads that sets defaults connections for ExtendedStateSpace.

and the following design functions expect ExtendedStateSpace as inputs

• hinfsynthesize
• h2synthesize
• LQGProblem (also accepts other types)

A video tutorial on how to use this type is available here.

se = ExtendedStateSpace(s::AbstractStateSpace; kwargs...)

The conversion from a regular statespace object to an ExtendedStateSpace creates the following system by default

\[\begin{bmatrix} A & B & B \\ C & D & D \\ C & D & D \end{bmatrix}\]

i.e., the system and performance mappings are identical, system_mapping(se) == performance_mapping(se) == s. However, all matrices B1, B2, C1, C2; D11, D12, D21, D22 are overridable by a corresponding keyword argument. In this case, the controlled outputs are the same as measured outputs.

Related: se = convert(ExtendedStateSpace, s::StateSpace) produces an ExtendedStateSpace with empty performance_mapping from w->z such that ss(se) == s.

G = LQGProblem(sys::ExtendedStateSpace, Q1, Q2, R1, R2; qQ=0, qR=0, SQ=nothing, SR=nothing)

Return an LQG object that describes the closed control loop around the process sys=ss(A,B,C,D) where the controller is of LQG-type. The controller is specified by weight matrices Q1,Q2 that penalizes state deviations and control signal variance respectively, and covariance matrices R1,R2 which specify state drift and measurement covariance respectively.

sys is an extended statespace object where the upper channel corresponds to disturbances to performance variables (w→z), and the lower channel corresponds to inputs to outputs (u→y), such that lft(sys, K) forms the closed-loop transfer function from external inputs/disturbances to performance variables.

qQ and qR can be set to incorporate loop-transfer recovery, i.e.,

L = lqr(A, B, Q1+qQ*C'C, Q2)
K = kalman(A, C, R1+qR*B*B', R2)

Increasing qQ will add more cost in output direction, e.g., encouraging the use of cheap control, while increasing qR adds fictious dynamics noise, makes the observer faster in the direction we control.

Example

In this example we will control a MIMO system with one unstable pole and one unstable zero. When the system has both unstable zeros and poles, there are fundamental limitations on performance. The unstable zero is in this case faster than the unstable pole, so the system is controllable. For good performance, we want as large separation between the unstable zero dynamics and the unstable poles as possible.

s = tf("s")
P = [1/(s+1) 2/(s+2); 1/(s+3) 1/(s-1)]
sys = ExtendedStateSpace(ss(P)) # Controlled outputs same as measured outputs and state noise affects at inputs only. 
eye(n) = Matrix{Float64}(I,n,n) # For convinience

qQ = 0
qR = 0
Q1 = 10eye(2)
Q2 = 1eye(2)
R1 = 1eye(2)
R2 = 1eye(2)

G = LQGProblem(sys, Q1, Q2, R1, R2, qQ=qQ, qR=qR)

T = comp_sensitivity(G)
S = sensitivity(G)
Gcl = closedloop(G)*static_gain_compensation(G)
plot(
    sigmaplot([S,T, Gcl],exp10.(range(-3, stop=3, length=1000)), lab=["S" "T" "Gry"]),
    plot(step(Gcl, 5))
)

Extended help

Several functions are defined for instances of LQGProblem

• closedloop
• extended_controller
• ff_controller
• gangoffour
• G_CS
• G_PS
• input_comp_sensitivity
• input_sensitivity
• output_comp_sensitivity
• output_sensitivity
• system_mapping
• performance_mapping
• static_gain_compensation
• gangoffourplot
• kalman
• lft
• lqr
• observer_controller

A video tutorial on how to use the LQG interface is available here

Introduction of references

The most principled way of introducing references is to add references as measured inputs to the extended statespace model, and to let the performance output z be the differences between the references and the outputs for which references are provided.

A less cumbersome way is not not consider references when constructing the LQGProblem, and instead pass the z keyword arugment to extended_controller in order to obtain a closed-loop system from state references to controlled outputs, and use some form of inverse of the DC gain of this system (or one of its subsystems) to pre-compensate the reference input.

LQGProblem(P::ExtendedStateSpace)

If only an ExtendedStateSpace system is provided, e.g. from hinfpartition, the system P is assumed to correspond to the H₂ optimal control problem with

C1'C1    = Q1
D12'D12  = Q2
SQ       = C1'D12 # Cross term

B1*B1'   = R1
D21*D21' = R2
SR       = B1*D21' # Cross term

and an LQGProblem with the above covariance matrices is returned. The system description in the returned LQGProblem will have B1 = C1 = I. See Ch. 14 in Robust and optimal control for reference.

Example:

All the following ways of obtaining the H2 optimal controller are (almost) equivalent

using Test
G = ss(tf(1, [10, 1]))
WS = tf(1, [1, 1e-6]) 
WU = makeweight(1e-2, 0.1, 100) 
Gd = hinfpartition(G, WS, WU, [])

K, Gcl = h2synthesize(Gd)              # First option, using H2 formulas
K2, Gcl2 = h2synthesize(Gd, 1000)      # Second option, using H∞ formulas with large γ

lqg = LQGProblem(Gd)                   # Third option, convert to an LQGProblem and obtain controller
K3 = -observer_controller(lqg)

@test h2norm(lft(Gd, K )) ≈ 3.0568 atol=1e-3
@test h2norm(lft(Gd, K2)) ≈ 3.0568 atol=1e-3
@test h2norm(lft(Gd, K3)) ≈ 3.0568 atol=1e-3

See named_ss for a convenient constructor.

UncertainSS{TE} <: AbstractStateSpace{TE}

Represents LFT_u(M, Diagonal(Δ))

nyquistcircles(w, centers, radii)

Plot the nyquist curve with circles. It only makes sense to call this function if the circles represent additive uncertainty, i.e., were calculated with relative=false.

See also fit_complex_perturbations

δ(N=32)

Create an uncertain element of N uniformly distributed samples ∈ [-1, 1]

G_CS(l::LQGProblem)

G_PS(l::LQGProblem)

input_comp_sensitivity(l::LQGProblem)

input_sensitivity(l::LQGProblem)

output_comp_sensitivity(l::LQGProblem)

output_sensitivity(l::LQGProblem)

ss(A, B1, B2, C1, C2, D11, D12, D21, D22 [, Ts])

Create an ExtendedStateSpace.

DescriptorSystems.dss(sys::AbstractStateSpace)

Convert sys to a descriptor statespace system from DescriptorSystems.jl

add_disturbance(sys::StateSpace, Ad::Matrix, Cd::Matrix)

See CCS pp. 144

Arguments:

• sys: System to augment
• Ad: The dynamics of the disturbance
• Cd: How the disturbance states affect the states of sys. This matrix has the shape (sys.nx, size(Ad, 1))

See also add_low_frequency_disturbance, add_resonant_disturbance

add_input_differentiator(sys::StateSpace, ui = 1:sys.nu; goodwin=false)

Augment the output of sys with the difference u(k+1)-u(k)

Arguments:

• ui: An index or vector of indices indicating which inputs to differentiate.
• goodwin: If true, the difference operator will use the Goodwin δ operator, i.e., (u(k+1)-u(k)) / sys.Ts.

The augmented system will have the matrices

[A 0; 0 0]  [B; I]  [C 0; 0 -I]  [D; I]

with length(ui) added states and outputs.

add_input_integrator(sys::StateSpace, ui = 1, ϵ = 0)

Augment the output of sys with the integral of input at index ui, i.e., y_aug = [y; ∫u[ui]] See also add_low_frequency_disturbance

add_low_frequency_disturbance(sys::StateSpace; ϵ = 0, measurement = false)
add_low_frequency_disturbance(sys::StateSpace, Cd; ϵ = 0, measurement = false)

Augment sys with a low-frequency (integrating if ϵ=0) disturbance model. If an integrating input disturbance is used together with an observer, the controller will have integral action.

• Cd: If adding an input disturbance. this matrix indicates how the disturbance states affect the states of sys, and defaults to sys.B. If measurement=true, this matrix indicates how the disturbance states affect the outputs of sys, and defaults to I(sys.ny).

Arguments:

• ϵ: Move the integrator pole ϵ into the stable region.
• measurement: If true, the disturbance is a measurement disturbance, otherwise it's an input diturbance.

add_low_frequency_disturbance(sys::StateSpace, Ai::Integer; ϵ = 0)

A disturbance affecting only state Ai.

add_measurement_disturbance(sys::StateSpace{Continuous}, Ad::Matrix, Cd::Matrix)

Create the system

Ae = [A 0; 0 Ad]
Be = [B; 0]
Ce = [C Cd]

add_differentiator(sys::StateSpace{<:Discrete})

Augment the output of sys with the numerical difference (discrete-time derivative) of output, i.e., y_aug = [y; (y-y_prev)/sys.Ts] To add both an integrator and a differentiator to a SISO system, use

Gd = add_output_integrator(add_output_differentiator(G), 1)

add_output_integrator(sys::StateSpace{<:Discrete}, ind = 1; ϵ = 0)

Augment the output of sys with the integral of output at index ind, i.e., y_aug = [y; ∫y[ind]] To add both an integrator and a differentiator to a SISO system, use

Gd = add_output_integrator(add_output_differentiator(G), 1)

Note: numerical integration is subject to numerical drift. If the output of the system corresponds to, e.g., a velocity reference and the integral to position reference, consider methods for mitigating this drift.

add_resonant_disturbance(sys::AbstractStateSpace, ω, ζ, Bd::AbstractArray)

• Bd: The disturbance input matrix.

add_resonant_disturbance(sys::StateSpace{Continuous}, ω, ζ, Ai::Int; measurement = false)

Augment sys with a resonant disturbance model.

Arguments:

• ω: Frequency
• ζ: Relative damping.
• Ai: The affected state
• measurement: If true, the disturbace is acting on the output, this will cause the controller to have zeros at ω (roots of poly s² + 2ζωs + ω²). If false, the disturbance is acting on the input, this will cause the controller to have poles at ω (roots of poly s² + 2ζωs + ω²).

sysr, hs = baltrunc2(sys::LTISystem; residual=false, n=missing, kwargs...)

Compute the a balanced truncation of order n and the hankel singular values

For keyword arguments, see the docstring of DescriptorSystems.gbalmr, reproduced below

gbalmr(sys, balance = false, matchdc = false, ord = missing, offset = √ϵ,
       atolhsv = 0, rtolhsv = nϵ, atolmin = atolhsv, rtolmin = rtolhsv, 
       atol = 0, atol1 = atol, atol2 = atol, rtol, fast = true) -> (sysr, hs)

Compute for a proper and stable descriptor system sys = (A-λE,B,C,D) with the transfer function matrix G(λ), a reduced order realization sysr = (Ar-λEr,Br,Cr,Dr) and the vector hs of decreasingly ordered Hankel singular values of the system sys. If balance = true, a balancing-based approach is used to determine a reduced order minimal realization of the form sysr = (Ar-λI,Br,Cr,Dr). For a continuous-time system sys, the resulting realization sysr is balanced, i.e., the controllability and observability grammians are equal and diagonal. If additonally matchdc = true, the resulting sysr is computed using state rezidualization formulas (also known as singular perturbation approximation) which additionally preserves the DC-gain of sys. In this case, the resulting realization sysr is balanced (for both continuous- and discrete-time systems). If balance = false, an enhanced accuracy balancing-free approach is used to determine the reduced order system sysr.

If ord = nr, the resulting order of sysr is min(nr,nrmin), where nrmin is the order of a minimal realization of sys determined as the number of Hankel singular values exceeding max(atolmin,rtolmin*HN), with HN, the Hankel norm of G(λ). If ord = missing, the resulting order is chosen as the number of Hankel singular values exceeding max(atolhsv,rtolhsv*HN).

To check the stability of the eigenvalues of the pencil A-λE, the real parts of eigenvalues must be less than -β for a continuous-time system or the moduli of eigenvalues must be less than 1-β for a discrete-time system, 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 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 ϵ is the working machine epsilon and n is the order of the system sys. The keyword argument atol can be used to simultaneously set atol1 = atol and atol2 = atol.

If E is singular, the uncontrollable infinite eigenvalues of the pair (A,E) and the non-dynamic modes are elliminated using minimal realization techniques. The rank determinations in the performed reductions are based on rank revealing QR-decompositions with column pivoting if fast = true or the more reliable SVD-decompositions if fast = false.

Method: For the order reduction of a standard system, the balancing-free method of [1] or the balancing-based method of [2] are used. For a descriptor system the balancing related order reduction methods of [3] are used. To preserve the DC-gain of the original system, the singular perturbation approximation method of [4] is used in conjunction with the balancing-based or balancing-free approach of [5].

References

[1] A. Varga. Efficient minimal realization procedure based on balancing. In A. El Moudni, P. Borne, and S.G. Tzafestas (Eds.), Prepr. of the IMACS Symp. on Modelling and Control of Technological Systems, Lille, France, vol. 2, pp.42-47, 1991.

[2] M. S. Tombs and I. Postlethwaite. Truncated balanced realization of a stable non-minimal state-space system. Int. J. Control, vol. 46, pp. 1319–1330, 1987.

[3] T. Stykel. Gramian based model reduction for descriptor systems. Mathematics of Control, Signals, and Systems, 16:297–319, 2004.

[4] Y. Liu Y. and B.D.O. Anderson Singular Perturbation Approximation of Balanced Systems, Int. J. Control, Vol. 50, pp. 1379-1405, 1989.

[5] Varga A. Balancing-free square-root algorithm for computing singular perturbation approximations. Proc. 30-th IEEE CDC, Brighton, Dec. 11-13, 1991, Vol. 2, pp. 1062-1065.

sysr, hs, info = baltrunc_coprime(sys; residual = false, n = missing, factorization::F = DescriptorSystems.gnlcf, kwargs...)

Compute a balanced truncation of the left coprime factorization of sys. See baltrunc2 for additional keyword-argument help.

Coprime-factor reduction performs a coprime factorization of the model into $P(s) = M(s)^{-1}N(s)$ where $M$ and $N$ are stable factors even if $P$ contains unstable modes. After this, the system $NM = \begin{bmatrix}N & M \end{bmatrix}$ is reduced using balanced truncation and the final reduced-order model is formed as $P_r(s) = M_r(s)^{-1}N_r(s)$. For this method, the Hankel signular values of $NM$ are reported and the reported errors are $||NM - N_rM_r||_\infty$. This method is of particular interest in closed-loop situations, where a model-reduction error $||NM - N_rM_r||_\infty$ no greater than the normalized-coprime margin of the plant and the controller, guaratees that the closed loop remains stable when either $P$ or $K$ are reduced. The normalized-coprime margin can be computed with ncfmargin(P, K) (ncfmargin).

Arguments:

• factorization: The function to perform the coprime factorization. A non-normalized factorization may be used by passing RobustAndOptimalControl.DescriptorSystems.glcf.
• kwargs: Are passed to DescriptorSystems.gbalmr

baltrunc_unstab(sys::LTISystem; residual = false, n = missing, kwargs...)

Balanced truncation for unstable models. An additive decomposition of sys into sys = sys_stable + sys_unstable is performed after which sys_stable is reduced. The order n must not be less than the number of unstable poles.

See baltrunc2 for other keyword arguments.

bilinearc2d(Ac::AbstractArray, Bc::AbstractArray, Cc::AbstractArray, Dc::AbstractArray, Ts::Number; tolerance=1e-12)

Balanced Bilinear transformation in State-Space. This method computes a discrete time equivalent of a continuous-time system, such that

\[G_d(z) = s2z[G_c(s)]\]

in a manner which accomplishes the following (i) Preserves the infinity L-infinity norm over the transformation (ii) Finds a system which balances B and C, in the sense that $||B||_2=||C||_2$ (iii) Satisfies $G_c(s) = z2s[s2z[G_c(s)]]$ for some map z2s[]

bilinearc2d(sys::ExtendedStateSpace, Ts::Number)

Applies a Balanced Bilinear transformation to a discrete-time extended statespace object

bilinearc2d(sys::StateSpace, Ts::Number)

Applies a Balanced Bilinear transformation to a discrete-time statespace object

bilineard2c(Ad::AbstractArray, Bd::AbstractArray, Cd::AbstractArray, Dd::AbstractArray, Ts::Number; tolerance=1e-12)

Balanced Bilinear transformation in State-Space. This method computes a continuous time equivalent of a discrete time system, such that

G_c(z) = z2s[G_d(z)]

in a manner which accomplishes the following (i) Preserves the infinity L-infinity norm over the transformation (ii) Finds a system which balances B and C, in the sense that ||B||2=||C||2 (iii) Satisfies Gd(z) = s2z[z2s[Gd(z)]] for some map s2z[]

bilineard2c(sys::ExtendedStateSpace)

Applies a Balanced Bilinear transformation to continuous-time extended statespace object

bilineard2c(sys::StateSpace)

Applies a Balanced Bilinear transformation to continuous-time statespace object

blocks, M = blocksort(P::UncertainSS)

Returns the block structure of P.Δ as well as P.M permuted according to the sorted block structure. blocks is a vector of vectors with the block structure of perturbation blocks as described by μ-tools, i.e.

• [-N, 0] denotes a repeated real block of size N
• [N, 0] denotes a repeated complex block of size N
• [ny, nu] denotes a full complex block of size ny × nu

broken_feedback(L, i)

Closes all loops in square MIMO system L except for loops i. Forms L1 in fig 14. of "An Introduction to Disk Margins", Peter Seiler, Andrew Packard, and Pascal Gahinet

closedloop(l::LQGProblem, L = lqr(l), K = kalman(l))

Closed-loop system as defined in Glad and Ljung eq. 8.28. Note, this definition of closed loop is not the same as lft(P, K), which has B1 instead of B2 as input matrix. Use lft(l) to get the system from disturbances to controlled variables w -> z.

The return value will be the closed loop from filtred reference only, other disturbance signals (B1) are ignored. See feedback for a more advanced option. This function assumes that the control signal is computed as u = r̃ - Lx̂ (not u = L(xᵣ - x̂)), i.e., the feedforward signal r̃ is added directly to the plant input. r̃ must thus be produced by an inverse-like model that takes state references and output the feedforward signal.

Use static_gain_compensation to adjust the gain from references acting on the input B2, dcgain(closedloop(l))*static_gain_compensation(l) ≈ I

connect(systems, connections; external_inputs, external_outputs = (:), verbose = true, unique = true, kwargs...)

Create block connections using named inputs and outputs.

Addition and subtraction nodes are achieved by creating a linear combination node, i.e., a system with a D matrix only.

Arguments:

• systems: A vector of named systems to be connected
• connections: a vector of pairs output => input, where each pair maps an output to an input. Each output must appear as an output in one of systems, and similarly each input must appear as an input in one of systems. All inputs must have unique names and so must all outputs, but an input may have the same name as an output. In the example below the connection :uP => :uP connects the output :uP of the addP block to P's input :uP
• external_inputs: external signals to be used as inputs in the constructed system. Use (:) to indicate all signals
• external_outputs: outputs of the constructed system. Use (:) to indicate all signals
• verbose: Issue warnings for signals that have no connection
• unique: If true, all input names must be unique. If false, a single external input signal may be connected to multiple input ports with the same name.

Note: Positive feedback is used, controllers that are intended to be connected with negative feedback must thus be negated.

Example: The following complicated feedback interconnection

                 yF
               ┌────────────────────────────────┐
               │                                │
    ┌───────┐  │  ┌───────┐        ┌───────┐    │    ┌───────┐
uF  │       │  │  │       | yR     │       │ yC │ uP │       │   yP
────►   F   ├──┴──►   R   │────+───►   C   ├────+────►   P   ├───┬────►
    │       │     │       │    │   │       │         │       │   │
    └───────┘     └───────┘    │   └───────┘         └───────┘   │
                               │                                 │
                               └─────────────────────────────────┘

can be created by

F = named_ss(ssrand(1, 1, 2, proper=true), x=:xF, u=:uF, y=:yF)
R = named_ss(ssrand(1, 1, 2, proper=true), x=:xR, u=:uR, y=:yR)
C = named_ss(ssrand(1, 1, 2, proper=true), x=:xC, u=:uC, y=:yC)
P = named_ss(ssrand(1, 1, 3, proper=true), x=:xP, u=:uP, y=:yP)

addP = sumblock("uP = yF + yC") # Sum node before P
addC = sumblock("uC = yR - yP") # Sum node before C

connections = [
    :yP => :yP # Output to input
    :uP => :uP
    :yC => :yC
    :yF => :yF
    :yF => :uR
    :uC => :uC
    :yR => :yR
]
external_inputs = [:uF] # External inputs

G = connect([F, R, C, P, addP, addC], connections; external_inputs)

If an external input is to be connected to multiple points, use a splitter to split up the signal into a set of unique names which are then used in the connections.

controller_reduction(P::ExtendedStateSpace, K, r, out=true; kwargs...)

Minimize ||(K-Kᵣ) W||∞ if out=false ||W (K-Kᵣ)||∞ if out=true See Robust and Optimal Control Ch 19.1 out indicates if the weight will be applied as output or input weight.

This function expects a *positive feedback controller K.

This method corresponds to the methods labelled SW1/SW2(SPA) in Andreas Varga, "Controller Reduction Using Accuracy-Enhancing Methods" SW1 is the default method, corresponding to out=true.

This method does not support unstable controllers. See the reference above for alternatives. See also stab_unstab and baltrunc_unstab.

controller_reduction_plot(G, K)

Plot the normalized-coprime margin (ncfmargin) as a function of controller order when baltrunc_coprime is used to reduce the controller. Red, orange and green bands correspond to rules of thumb for bad, okay and good coprime uncertainty margins. A value of 0 indicate an unstable closed loop.

If G is an ExtendedStateSpace system, a second plot will be shown indicating the $H_∞$ norm between inputs and performance outputs $||T_{zw}||_\infty$ when the function controller_reduction is used to reduce the controller.

The order of the controller can safely be reduced as long as the normalized coprime margin remains sufficiently large. If the controller contains integrators, it may be advicable to protect the integrators from the reduction, e.g., if the controller is obtained using glover_mcfarlane, perform the reduction on info.Gs, info.Ks rather than on K, and form Kr using the reduced Ks.

See glover_mcfarlane or the docs for an example.

controller_reduction_weight(P::ExtendedStateSpace, K)

Lemma 19.1 See Robust and Optimal Control Ch 19.1

dare3(P::AbstractStateSpace, Q1::AbstractMatrix, Q2::AbstractMatrix, Q3::AbstractMatrix; full=false)

Solve the discrete-time algebraic Riccati equation for a discrete LQR cost augmented with control differences

\[x^{T} Q_1 x + u^{T} Q_2 u + Δu^{T} Q_3 Δu, \quad Δu = u(k) - u(k-1)\]

If full, the returned matrix will include the state u(k-1), otherwise the returned matrix will be of the same size as Q1.

diskmargin(L, σ = 0)
diskmargin(L, σ::Real, ω)

Calculate the disk margin of LTI system L. L is supposed to be a loop-transfer function, i.e., it should be square. If L = PC the disk margin for output perturbations is computed, whereas if L = CP, input perturbations are considered. If the method diskmargin(P, C, args...) is used, both are computed. Note, if L is MIMO, a simultaneous margin is computed, see loop_diskmargin for single loop margins of MIMO systems.

The implementation and notation follows "An Introduction to Disk Margins", Peter Seiler, Andrew Packard, and Pascal Gahinet.

The margins are aviable as fields of the returned objects, see Diskmargin.

Arguments:

• L: A loop-transfer function.
• σ: If little is known about the distribution of gain variations then σ = 0 is a reasonable choice as it allows for a gain increase or decrease by the same relative amount. The choice σ < 0 is justified if the gain can decrease by a larger factor than it can increase. Similarly, the choice σ > 0 is justified when the gain can increase by a larger factor than it can decrease. If σ = −1 then the disk margin condition is αmax = inv(MT). This margin is related to the robust stability condition for models with multiplicative uncertainty of the form P (1 + δ). If σ = +1 then the disk margin condition is αmax = inv(MS)
• kwargs: Are sent to the hinfnorm calculation
• ω: If a vector of frequencies is supplied, the frequency-dependent disk margin will be computed, see example below.

Example:

L = tf(25, [1,10,10,10])
dm = diskmargin(L, 0)
plot(dm) # Plot the disk margin to illustrate maximum allowed simultaneous gain and phase variations.

nyquistplot(L)
plot!(dm, nyquist=true) # plot a nyquist exclusion disk. The Nyquist curve will be tangent to this disk at `dm.ω0`
nyquistplot!(dm.f0*L) # If we perturb the system with the worst-case perturbation `f0`, the curve will pass through the critical point -1.

## Frequency-dependent margin
w = exp10.(LinRange(-2, 2, 500))
dms = diskmargin(L, 0, w)
plot(dms; lower=true, phase=true)

Example: relation to Ms and Mt

Ms, wMs = hinfnorm(input_sensitivity(P, C)) # Input Ms
dm = diskmargin(C*P, 1) # Input diskmargin, skew = +1
isapprox(Ms/(Ms-1), dm.gainmargin[2], rtol=1e-2) # Guaranteed gain margin based on Ms
isapprox(inv(Ms), dm.margin, rtol=1e-2)
isapprox(dm.ω0, wMs, rtol=1e-1)


Mt, wMt = hinfnorm(input_comp_sensitivity(P, C)) # Input Mt
dm = diskmargin(C*P, -1) # Input diskmargin, skew = -1
isapprox(inv(Mt), dm.margin, rtol=1e-2)
isapprox(dm.ω0, wMt, rtol=1e-1)

See also ncfmargin and loop_diskmargin.

diskmargin(P::LTISystem, C::LTISystem, σ, w::AbstractVector, args...; kwargs...)

Simultaneuous diskmargin at outputs, inputs and input/output simultaneously of P. Returns a named tuple with the fields input, output, simultaneous_input, simultaneous_output, simultaneous where input and output represent loop-at-a-time margins, simultaneous_input is the margin for simultaneous perturbations on all inputs and simultaneous is the margin for perturbations on all inputs and outputs simultaneously.

Note: simultaneous margins are more conservative than single-loop margins and are likely to be much lower than the single-loop margins. Indeed, with several simultaneous perturbations, it's in general easier to make the system unstable. It's not uncommon for a simultaneous margin involving two signals to be on the order of half the size of the single-loop margins.

See also ncfmargin and loop_diskmargin.

diskmargin(L::LTISystem, σ::Real, ω)

Calculate the diskmargin at a particular frequency or vector of frequencies. If ω is a vector, you get a frequency-dependent diskmargin plot if you plot the returned value. See also ncfmargin.

expand_symbol(s::Symbol, n::Int)

Takes a symbol and an integer and returns a vector of symbols with increasing numbers appended to the end. E.g., (:x, 3) -> [:x1, :x2, :x3]

The short-hand syntax s^n is also available, e.g., :x^3 == expand_symbol(:x, 3).

Useful to create signal names for named systems.

extended_controller(l::LQGProblem, L = lqr(l), K = kalman(l); z = nothing)

Returns a statespace system representing the controller that is obtained when state-feedback u = L(xᵣ-x̂) is combined with a Kalman filter with gain K that produces state estimates x̂. The controller is an instance of ExtendedStateSpace where C2 = -L, D21 = L and B2 = K.

The returned system has inputs [xᵣ; y] and outputs the control signal u. If a reference model R is used to generate state references xᵣ, the controller from (ry, y) -> u where ry - y = e is given by

Ce = extended_controller(l)
Ce = named_ss(ss(Ce); x = :xC, y = :u, u = [R.y; :y^l.ny]) # Name the inputs of Ce the same as the outputs of `R`.
connect([R, Ce]; u1 = R.y, y1 = R.y, w1 = [R.u; :y^l.ny], z1=[:u])

Since the negative part of the feedback is built into the returned system, we have

C = observer_controller(l)
Ce = extended_controller(l)
system_mapping(Ce) == -C

Please note, without the reference pre-filter, the DC gain from references to controlled outputs may not be identity. If a vector of output indices is provided through the keyword argument z, the closed-loop system from state reference xᵣ to outputs z is returned as a second return argument. The inverse of the DC-gain of this closed-loop system may be useful to compensate for the DC-gain of the controller.

extended_controller(K::AbstractStateSpace)

Takes a controller and returns an ExtendedStateSpace version which has augmented input [r; y] and output y (z output is 0-dim).

G = feedback_control(P, K)

Return the (negative feedback) closed-loop system from input of K to output of P while outputing also the control signal (output of K), i.e., G maps references to [y; u]

Example:

The following are two equivalent ways of achieving the same thing

G = ssrand(3,4,2)
K = ssrand(4,3,2)

Gcl1 = feedback_control(G, K) # First option

# Second option using named systems and connect
G = named_ss(G, :G)
K = named_ss(K, :K)
S = sumblock("Ku = r - Gy", n=3) # Create a sumblock that computes r - Gy for vectors of length 3

z1 = [G.y; K.y] # Output both plant and controller outputs
w1 = :r^3       # Extenal inputs are the three references into the sum block
connections = [K.y .=> G.u; G.y .=> G.y; K.u .=> K.u] # Since the sumblock uses the same names as the IO signals of G,K, we can reuse these names here
Gcl2 = connect([G, K, S], connections; z1, w1)

@test linfnorm(minreal(Gcl1 - Gcl2.sys))[1] < 1e-10 # They are the same

To include also an input disturbance, use

Gcl = feedback(K, P, W2=:, Z2=:, Zperm=[(1:ny).+nu; 1:nu]) # y,u from r,d

See also extended_gangoffour.

ff_controller(l::LQGProblem, L = lqr(l), K = kalman(l))

Return the feedforward controller $C_{ff}$ that maps references to plant inputs: $u = C_{fb}y + C_{ff}r$

The following should hold

Cff = RobustAndOptimalControl.ff_controller(l)
Cfb = observer_controller(l)
Gcl = feedback(system_mapping(l), Cfb) * Cff # Note the comma in feedback, P/(I + PC) * Cff
dcgain(Gcl) ≈ I # Or some I-like non-square matrix 

Note, if extended_controller is used, the DC-gain compensation above cannot be used. The extended_controller assumes that the references enter like u = L(xᵣ - x̂).

See also observer_controller.

l, res = find_lft(sys::StateSpace{<:Any, <:StaticParticles{<:Any, N}}, δ) where N

NOTE: This function is experimental.

Given an systems sys with uncertain coefficients in the form of StaticParticles, search for a lower linear fractional transformation M such that lft(M, δ) ≈ sys.

δ can be either the source of uncertainty in sys, i.e., a vector of the unique uncertain parameters that were used to create sys. These should be constructed as uniform randomly distributed particles for most robust-control theory to be applicable. δ can also be an integer, in which case a numer of δ sources of uncertainty are automatically created. This could be used for order reduction if the number of uncertainty sources in sys is large.

Note, uncertainty in sys is only supported in A and B, C and D must be deterministic.

Returns l::LFT that internaly contains all four blocks of M as well as δ. Call ss(l,sys) do obtain lft(M, δ) ≈ sys.

Call Matrix(l) to obtain M = [M11 M12; M21 M22]

centers, radii = fit_complex_perturbations(P, w; relative=true, nominal=:mean)

For each frequency in w, fit a circle in the complex plane that contains all models in the model set P, represented as an LTISystem with Particles coefficients. Note, the resulting radii can be quite unstable if the number of particles is small, in particular if the particles are normally distributed instead of uniformly.

If realtive = true, circles encompassing |(P - Pn)/Pn| will be returned (multiplicative/relative uncertainty). If realtive = false, circles encompassing |P - Pn| will be returned (additive uncertainty).

If nominal = :mean, the mean of P will be used as nominal model. If nominal = :first, the first particle will be used. See MonteCarloMeasurements.with_nominal to set the nominal value in the first particle. If nominal = :center, the middle point (pmaximum(ri)+pminimum(ri))/2 will be used. This typically gives the smallest circles.

See also nyquistcircles to plot circles (only if relative=false).

frequency_separation(sys, ω)

Decomponse sys into sys = sys_slow + sys_fast, where sys_slow contain all modes with eigenvalues with absolute value less than ω and sys_fast contain all modes with eigenvalues with absolute value greater than or equal to ω.

sysr, hs = frequency_weighted_reduction(G, Wo, Wi, r=nothing; residual=true)

Find Gr such that $||Wₒ(G-Gr)Wᵢ||∞$ is minimized. For a realtive reduction, set Wo = inv(G) and Wi = I.

If residual = true, matched static gain is achieved through "residualization", i.e., setting

\[0 = A_{21}x_{1} + A_{22}x_{2} + B_{2}u\]

where indices 1/2 correspond to the remaining/truncated states respectively. This choice typically results in a better match in the low-frequency region and a smaller overall error.

Ref: Andras Varga and Brian D.O. Anderson, "Accuracy enhancing methods for the frequency-weighted balancing related model reduction"

Note: This function only handles exponentially stable models. To reduce unstable and marginally stable models, decompose the system into stable and unstable parts using stab_unstab, reduce the stable part and then add the unstable part back.

Example:

The following example performs reduction with a frequency focus between frequencies w1 and w2.

using DSP
w1 = 1e-4
w2 = 1e-1
wmax = 1
fc = DSP.analogfilter(DSP.Bandpass(w1, w2, fs=wmax), DSP.Butterworth(2))
tfc = DSP.PolynomialRatio(fc)
W = tf(DSP.coefb(tfc), DSP.coefa(tfc))
rsys, hs = frequency_weighted_reduction(sys, W, 1)

fudge_inv(s::AbstractStateSpace, ε = 0.001)

Allow inverting a proper statespace system by adding a tiny (ε) feedthrough term to the D matrix. The system must still be square.

gain_and_delay_uncertainty(kmin, kmax, Lmax)

Return a multiplicative weight to represent the uncertainty coming from neglecting the dynamics k*exp(-s*L) where k ∈ [kmin, kmax] and L ≤ Lmax. This weight is slightly optimistic, an expression for a more exact weight appears in eq (7.27), "Multivariable Feedback Control: Analysis and Design"

See also neglected_lag and neglected_delay.

Example:

a = 10
P = ss([0 a; -a 0], I(2), [1 a; -a 1], 0) # Plant
W0 = gain_and_delay_uncertainty(0.5, 2, 0.005) |> ss # Weight
W = I(2) + W0*I(2) * uss([δc(), δc()]) # Create a diagonal real uncertainty weighted in frequency by W0
Ps = P*W # Uncertain plant
Psamples = rand(Ps, 500) # Sample the uncertain plant for plotting
w = exp10.(LinRange(-1, 3, 300)) # Frequency vector
bodeplot(Psamples, w)

gainphaseplot(P)
gainphaseplot(P, re, im)

Plot complex perturbantions to the plant P and indicate whether or not the closed-loop system is stable. The diskmargin is the largest disk that can be fit inside the green region that only contains stable variations.

K, γ, info = glover_mcfarlane(G::AbstractStateSpace, γ = 1.1; W1=1, W2=1)

Design a controller for G that maximizes the stability margin ϵ = 1/γ with normalized coprime factor uncertainty using the method of Glover and McFarlane

γ = 1/ϵ = ||[K;I] inv(I-G*K)*inv(M)||∞
G = inv(M + ΔM)*(N + ΔN)

γ is given as a relative factor above γmin and must be greater than 1, i.e., if γ = 1.1, the controller will be designed for γ = 1.1*γmin.

We want γmin (which is always ≥ 1) as small as possible, and we usually require that γmin is less than 4, corresponding to 25% allowed coprime uncertainty.

Performance modeling is incorporated in the design by calling glover_mcfarlane on the shaped system Gs = W2*G*W1 and then forming the controller as K = W1*Ks*W2. Using this formulation, traditional loop shaping can be done on Gs = W2*G*W1. The plant shaping is handled internally if keyword arguments W1, W2 are used and the returned controller is already scaled. In this case, Gs and Ks are included in the info named tuple for inspection.

See also glover_mcfarlane_2dof to design a feedforward filter as well and baltrunc_coprime for controller order reduction. When reducing the order of the calculated controller, reduce the order of info.Ks and form Kr=W1*Ksred*W2. Verify the robustness using ncfmargin(info.Gs, Ksred) as well as ncfmargin(G, Kr).

Example:

Example 9.3 from the reference below.

using RobustAndOptimalControl, ControlSystemsBase, Plots, Test
G = tf(200, [10, 1])*tf(1, [0.05, 1])^2     |> ss
Gd = tf(100, [10, 1])                       |> ss
W1 = tf([1, 2], [1, 1e-6])                  |> ss
K, γ, info = glover_mcfarlane(G, 1.1; W1)
@test info.γmin ≈ 2.34 atol=0.005
Gcl = extended_gangoffour(G, K) # Form closed-loop system

fig1 = bodeplot([G, info.Gs, G*K], lab=["G" "" "G scaled" "" "Loop transfer"])
fig2 = bodeplot(Gcl, lab=["S" "KS" "PS" "T"], plotphase=false) # Plot gang of four

fig3 = plot(step(Gd*feedback(1, info.Gs), 3), lab="Initial controller")
plot!(step(Gd*feedback(1, G*K), 3), lab="Robustified")
fig4 = nyquistplot([info.Gs, G*K], ylims=(-2,1), xlims=(-2, 1),
    Ms_circles = 1.5,
    lab = ["Initial controller" "Robustified"],
    title = "Loop transfers with and without robustified controller"
)
plot(fig1, fig2, fig3, fig4)

Example of controller reduction: The order of the controller designed above can be reduced maintaining at least 2/3 of the robustness margin like this

e,_ = ncfmargin(info.Gs, info.Ks)
Kr, hs, infor = baltrunc_coprime(info.Ks, n=info.Ks.nx)
n = findlast(RobustAndOptimalControl.error_bound(hs) .> 2e/3) # 2/3 e sets the robustness margin
Ksr, hs, infor = baltrunc_coprime(info.Ks; n)
@test ncfmargin(info.Gs, Ksr)[1] >= 2/3 * e
Kr = W1*Ksr
bodeplot([G*K, G*Kr], lab=["L original" "" "L Reduced" ""])

This gives a final controller Kr of order 3 instead of order 5, but a very similar robustness margin. You may also call

controller_reduction_plot(info.Gs, info.Ks)

to help you select the controller order.

Ref: Sec 9.4.1 of Skogestad, "Multivariable Feedback Control: Analysis and Design"

Extended help

Skogestad gives the following general advice:

1. Scale the plant outputs and inputs. This is very important for most design procedures and is sometimes forgotten. In general, scaling improves the conditioning of the design problem, it enables meaningful analysis to be made of the robustness properties of the feedback system in the frequency domain, and for loop-shaping it can simplify the selection of weights. There are a variety of methods available including normalization with respect to the magnitude of the maximum or average value of the signal in question. If one is to go straight to a design the following variation has proved useful in practice:

   • The outputs are scaled such that equal magnitudes of cross-coupling into each of the outputs is equally undesirable.
   • Each input is scaled by a given percentage (say 10%) of its expected range of operation. That is, the inputs are scaled to reflect the relative actuator capabilities.

2. Order the inputs and outputs so that the plant is as diagonal as possible. The relative gain array relative_gain_array can be useful here. The purpose of this pseudo-diagonalization is to ease the design of the pre- and post-compensators which, for simplicity, will be chosen to be diagonal.

   Next, we discuss the selection of weights to obtain the shaped plant $G_s = W_2 G W_1$ where $W_1 = W_p W_a W_g$

3. Select the elements of diagonal pre- and post-compensators $W_p$ and $W_2$ so that the singular values of $W_2 G W_p$ are desirable. This would normally mean high gain at low frequencies, roll-off rates of approximately 20 dB/decade (a slope of about 1) at the desired bandwidth(s), with higher rates at high frequencies. Some trial and error is involved here. $W_2$ is usually chosen as a constant, reflecting the relative importance of the outputs to be controlled and the other measurements being fed back to the controller. For example, if there are feedback measurements of two outputs to be controlled and a velocity signal, then $W_2$ might be chosen to be diag([1, 1, 0.1]), where 0.1 is in the velocity signal channel. $W_p$ contains the dynamic shaping. Integral action, for low frequency performance; phase-advance for reducing the roll-off rates at crossover, and phase-lag to increase the roll-off rates at high frequencies should all be placed in $W_p$ if desired. The weights should be chosen so that no unstable hidden modes are created in $G_s$.

4. Optional: Introduce an additional gain matrix $W_g$ cascaded with $W_a$ to provide control over actuator usage. $W_g$ is diagonal and is adjusted so that actuator rate limits are not exceeded for reference demands and typical disturbances on the scaled plant outputs. This requires some trial and error.

5. Robustly stabilize the shaped plant $G_s = W_2 G W_1$ , where $W_1 = W_p W_a W_g$, using glover_mcfarlane. First, the maximum stability margin $ϵ_{max} = 1/γ_{min}$ is calculated. If the margin is too small, $ϵmax < 0.25$, then go back and modify the weights. Otherwise, a γ-suboptimal controller is synthesized. There is usually no advantage to be gained by using the optimal controller. When $ϵ_{max}$ > 0.25 (respectively $γ_{min}$ < 4) the design is usually successful. In this case, at least 25% coprime factor uncertainty is allowed, and we also find that the shape of the open-loop singular values will not have changed much after robust stabilization. A small value of ϵmax indicates that the chosen singular value loop-shapes are incompatible with robust stability requirements. That the loop-shapes do not change much following robust stabilization if γ is small (ϵ large), is justified theoretically in McFarlane and Glover (1990).

6. Analyze the design and if all the specifications are not met make further modifications to the weights.

7. Implement the controller. The configuration shown in below has been found useful when compared with the conventional set up. This is because the references do not directly excite the dynamics of $K$, which can result in large amounts of overshoot (classical derivative kick). The constant prefilter ensures a steady-state gain of 1 between r and y, assuming integral action in $W_1$ or $G$ (note, the K returned by this function has opposite sign compared to that of Skogestad, so we use negative feedback here).

Anti-windup can be added to $W_1$ but putting $W_1$ on Hanus form after the synthesis, see hanus.

       ┌─────────┐      ┌────────┐      ┌────────┐
    r  │         │    us│        │  u   │        │  y
   ───►│(K*W2)(0)├──+──►│   W1   ├─────►│   G    ├────┬──►
       │         │  │-  │        │      │        │    │
       └─────────┘  │   └────────┘      └────────┘    │
                    │                                 │
                    │                                 │
                    │   ┌────────┐      ┌────────┐    │
                    │   │        │  ys  │        │    │
                    └───┤   K    │◄─────┤   W2   │◄───┘
                        │        │      │        │
                        └────────┘      └────────┘

Keywords: nfcsyn, coprimeunc

K, γ, info = glover_mcfarlane(G::AbstractStateSpace{<:Discrete}, γ = 1.1; W1=1, W2=1, strictly_proper=false)

For discrete systems, the info tuple contains also feedback gains F, L and observer gain Hkf such that the controller on observer form is given by

\[x^+ = Ax + Bu + H_{kf} (Cx - y)\\ u = Fx + L (Cx - y)\]

Note, this controller is not strictly proper, i.e., it has a non-zero D matrix. The controller can be transformed to observer form for the scaled plant (info.Gs) by Ko = observer_controller(info), in which case the following holds G*K == info.Gs*Ko (realizations are different).

If strictly_proper = true, the returned controller K will have D == 0. This can be advantageous in implementations where computational delays are present. In this case, info.L == 0 as well.

Ref discrete version: Iglesias, "The Strictly Proper Discrete-Time Controller for the Normalized Left-Coprime Factorization Robust Stabilization Problem"

K, γ, info = glover_mcfarlane_2dof(G::AbstractStateSpace{Continuous}, Tref::AbstractStateSpace{Continuous}, γ = 1.1, ρ = 1.1;
W1 = 1, Wo = I, match_dc = true, kwargs...)

Joint design of feedback and feedforward compensators

\[K = \begin{bmatrix} K_1 & K_2 \end{bmatrix}\]

   ┌──────┐   ┌──────┐        ┌──────┐    ┌─────┐
r  │      │   │      │        │      │    │     │
──►│  Wi  ├──►│  K1  ├───+───►│  W1  ├───►│  G  ├─┐y
   │      │   │      │   │    │      │    │     │ │
   └──────┘   └──────┘   │    └──────┘    └─────┘ │
                         │                        │
                         │    ┌──────┐            │
                         │    │      │            │
                         └────┤  K2  ◄────────────┘
                              │      │
                              └──────┘

Where the returned controller K takes the measurement vector [r; y] (positive feedback), i.e., it includes all blocks Wi, K1, K2, W1. If match_dc = true, Wi is automatically computed to make sure the static gain matches Tref exactly, otherwise Wi is set to I. The info named tuple contains the feedforward filter for inspection (info.K1 = K1*Wi).

Arguments:

• G: Plant model
• Tref: Reference model
• γ: Relative γ
• ρ: Design parameter, typically 1 < ρ < 3. Increase to emphasize model matching at the expense of robustness.
• W1: Pre-compensator for loop shaping.
• Wo: Output selction matrix. If there are more measurements than controlled variables, this matrix let's you select which measurements are to be controlled.
• kwargs: Are sent to hinfsynthesize.

Ref: Sec. 9.4.3 of Skogestad, "Multivariable Feedback Control: Analysis and Design". The reference contains valuable pointers regarding gain-scheduling implementation of the designed controller as an observer with feedback from estimated states. In order to get anti-windup protection when W1 contains an integrator, transform W1 to self-conditioned Hanus form (using hanus) and implement the controller like this

W1h = hanus(W1)             # Perform outside loop

# Each iteration
us = filter(Ks, [r; y])     # filter inputs through info.Ks (filter is a fictive function that applies the transfer function)
u  = filter(W1h, [us; ua])  # filter us and u-actual (after input saturation) through W1h
ua = clamp(u, lower, upper) # Calculate ua for next iteration as the saturated value of u

Example:

using RobustAndOptimalControl, Plots
P = tf([1, 5], [1, 2, 10]) # Plant
W1 = tf(1,[1, 0]) |> ss    # Loop shaping controller

Tref = tf(1, [1, 1])^2 |> ss # Reference model (preferably of same order as P)

K1dof, γ1, info1 = glover_mcfarlane(ss(P), 1.1; W1)
K2dof, γ2, info2 = glover_mcfarlane_2dof(ss(P), Tref, 1.1, 1.1; W1)

G1 = feedback(P*K1dof)
G2 = info2.Gcl

w = exp10.(LinRange(-2, 2, 200))
bodeplot(info2.K1, w, lab="Feedforward filter")
plot([step(G1, 15), step(G2, 15), step(Tref, 15)], lab=["1-DOF" "2-DOF" "Tref"])

n = h2norm(sys::LTISystem; kwargs...)

A numerically robust version of norm using DescriptorSystems.jl

For keyword arguments, see the docstring of DescriptorSystems.gh2norm, reproduced below

gh2norm(sys, fast = true, offset = sqrt(ϵ), atol = 0, atol1 = atol, atol2 = atol, atolinf = atol, rtol = n*ϵ)

Compute for a descriptor system sys = (A-λE,B,C,D) the H2 norm of its transfer function matrix G(λ). The H2 norm is infinite, if G(λ) has unstable poles, or, for a continuous-time, the system has nonzero gain at infinity. If the pencil A-λE has uncontrollable and/or unobservable unstable eigenvalues on the boundary of the stability domain, then a reduced order realization is determined first (see below) to eliminate these eigenvalues.

To check the stability, the eigenvalues of the pole pencil A-λE must have real parts less than -β for a continuous-time system or have moduli less than 1-β for a discrete-time system, 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.

For a continuous-time system sys with E singular, a reduced order realization is determined first, without uncontrollable and unobservable finite and infinite eigenvalues of the pencil A-λE. For a discrete-time system or for a system with invertible E, a reduced order realization is determined first, without uncontrollable and unobservable finite eigenvalues of the pencil A-λE. The rank determinations in the performed reductions are based on rank revealing QR-decompositions with column pivoting if fast = true or the more reliable SVD-decompositions if fast = 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 keyword argument atolinf is the absolute tolerance for the gain of G(λ) at λ = ∞. The used default value is atolinf = 0. The default relative tolerance is n*ϵ, where ϵ is the working machine epsilon and n is the order of the system sys. The keyword argument atol can be used to simultaneously set atol1 = atol and atol2 = atol.

K, Cl = h2synthesize(P::ExtendedStateSpace, γ = nothing)

Synthesize H₂-optimal controller K and calculate the closed-loop transfer function from w to z. Ref: Cha. 14.5 in Robust and Optimal Control.

If γ = nothing, use the formulas for H₂ in Ch 14.5. If γ is a large value, the H∞ formulas are used. As γ → ∞, these two are equivalent. The h∞ formulas do a coordinate transfromation that handles slightly more general systems so if you run into an error, it might be worth trying setting γ to something large, e.g., 1000.

n, hsv = hankelnorm(sys::LTISystem; kwargs...)

Compute the hankelnorm and the hankel singular values

For keyword arguments, see the docstring of DescriptorSystems.ghanorm, reproduced below

ghanorm(sys, fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol = n*ϵ) -> (hanorm, hs)

Compute for a proper and stable descriptor system sys = (A-λE,B,C,D) with the transfer function matrix G(λ), the Hankel norm hanorm = $\small ||G(\lambda)||_H$ and the vector of Hankel singular values hs of the minimal realizatioj of the system.

For a non-minimal system, the uncontrollable and unobservable finite and infinite eigenvalues of the pair (A,E) and the non-dynamic modes are elliminated using minimal realization techniques. The rank determinations in the performed reductions are based on rank revealing QR-decompositions with column pivoting if fast = true or the more reliable SVD-decompositions if fast = 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 ϵ is the working machine epsilon and n is the order of the system sys. The keyword argument atol can be used to simultaneously set atol1 = atol and atol2 = atol.

Wh = hanus(W)

Return Wh on Hanus form. Wh has twice the number of inputs, where the second half of the inputs are "actual inputs", e.g., potentially saturated. This is used to endow W with anti-windup protection. W must have an invertable D matrix and be minimum phase.

Ref: Sec 9.4.5 of Skogestad, "Multivariable Feedback Control: Analysis and Design"

sysm, T, HF = hess_form(sys)

Bring sys to Hessenberg form form.

The Hessenberg form is characterized by A having upper Hessenberg structure. T is the similarity transform applied to the system such that

sysm ≈ similarity_transform(sys, T)

HF is the Hessenberg-factorization of A.

See also modal_form and schur_form

flag = hinfassumptions(P::ExtendedStateSpace; verbose=true)

Check the assumptions for using the γ-iteration synthesis in Theorem 1.

∇A, ∇B, ∇C, ∇D, hn, ω = hinfgrad(sys; rtolinf=1e-8, kwargs...)
∇A, ∇B, ∇C, ∇D        = hinfgrad(sys, hn, ω)

Compute the gradient of the H∞ norm w.r.t. the statespace matrices A,B,C,D. If only a system is provided, the norm hn and the peak frequency ω are automatically calculated. kwargs are sent to hinfnorm2. Note, the default tolerance to which the norm is calculated is set smaller than default for hinfnorm2, gradients will be discontinuous with any non-finite tolerance, and sensitive optimization algorithms may require even tighter tolerance.

In cases where the maximum singular value is reached at more than one frequency, a random frequency is used.

If the system is unstable, the gradients are NaN. Strategies to find an initial stabilizing controllers are outlined in Apkarian and D. Noll, "Nonsmooth H∞ Synthesis" in IEEE Transactions on Automatic Control.

An rrule for ChainRules is defined using this function, so hn is differentiable with any AD package that derives its rules from ChainRules (only applies to the hn return value, not ω).

n, ω = hinfnorm2(sys::LTISystem; kwargs...)

A numerically robust version of hinfnorm using DescriptorSystems.jl

For keyword arguments, see the docstring of DescriptorSystems.ghinfnorm, reproduced below

ghinfnorm(sys, rtolinf = 0.001, fast = true, offset = sqrt(ϵ), atol = 0, atol1 = atol, atol2 = atol, rtol = n*ϵ) -> (hinfnorm, fpeak)

Compute for a descriptor system sys = (A-λE,B,C,D) with the transfer function matrix G(λ) the H∞ norm hinfnorm (i.e., the peak gain of G(λ)) and the corresponding peak frequency fpeak, where the peak gain is achieved. The H∞ norm is infinite if G(λ) has unstable poles. If the pencil A-λE has uncontrollable and/or unobservable unstable eigenvalues, then a reduced order realization is determined first (see below) to eliminate these eigenvalues.

To check the stability, the eigenvalues of the pencil A-λE must have real parts less than -β for a continuous-time system or have moduli less than 1-β for a discrete-time system, 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.

For a continuous-time system sys with E singular, a reduced order realization is determined first, without uncontrollable and unobservable finite and infinite eigenvalues of the pencil A-λE. For a discrete-time system or for a system with invertible E, a reduced order realization is determined first, without uncontrollable and unobservable finite eigenvalues of the pencil A-λE. The rank determinations in the performed reductions are based on rank revealing QR-decompositions with column pivoting if fast = true or the more reliable SVD-decompositions if fast = false.

The keyword arguments atol1, atol2, and rtol, specify, respectively, the absolute tolerance for the nonzero elements of matrices 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 ϵ is the working machine epsilon and n is the order of the system sys. The keyword argument atol can be used to simultaneously set atol1 = atol and atol2 = atol.

P = hinfpartition(G, WS, WU, WT)

Transform a SISO or MIMO system $G$, with weighting functions $W_S, W_U, W_T$ into an LFT with an isolated controller, and write the resulting system, $P(s)$, on a state-space form. Valid inputs for $G$ are transfer functions (with dynamics, can be both MIMO and SISO, both in tf and ss forms). Valid inputs for the weighting functions are empty arrays, numbers (static gains), and LTISystems.

Note, system_mapping(P) is equal to -G.

Extended help

For ill-conditioned MIMO plants, the $S, CS, T$ weighting may result in controllers that "invert" the plant, which may result in poor robustness. For such systems, penalizing $GS$ and $T$ may be more appropriate. Ref: "Inverting and noninverting H∞ controllers", Urs Christen, Hans Geering

hinfsignals(P::ExtendedStateSpace, G::LTISystem, C::LTISystem)

Use the extended state-space model, a plant and the found controller to extract the closed loop transfer functions.

• Pcl : w → z : From input to the weighted functions
• S : w → e : From input to error
• CS : w → u : From input to control
• T : w → y : From input to output

K, γ, mats = hinfsynthesize(P::ExtendedStateSpace; gtol = 1e-4, interval = (0, 20), verbose = false, tolerance = 1.0e-10, γrel = 1.01, transform = true, ftype = Float64, check = true)

Computes an H-infinity optimal controller K for an extended plant P such that $||F_l(P, K)||∞ < γ$(lft(P, K)) for the smallest possible γ given P. The routine is known as the γ-iteration, and is based on the paper "State-space formulae for all stabilizing controllers that satisfy an H∞-norm bound and relations to risk sensitivity" by Glover and Doyle.

Arguments:

• gtol: Tolerance for γ.
• interval: The starting interval for the bisection.
• verbose: Print progress?
• tolerance: For detecting eigenvalues on the imaginary axis.
• γrel: If γrel > 1, the optimal γ will be found by γ iteration after which a controller will be designed for γ = γopt * γrel. It is often a good idea to design a slightly suboptimal controller, both for numerical reasons, but also since the optimal controller may contain very fast dynamics. If γrel → ∞, the computed controller will approach the 𝑯₂ optimal controller. Getting a mix between 𝑯∞ and 𝑯₂ properties is another reason to choose γrel > 1.
• transform: Apply coordiante transform in order to tolerate a wider range or problem specifications.
• ftype: construct problem matrices in higher precision for increased numerical robustness. If the calculated controller achieves
• check: Perform a post-design check of the γ value achieved by the calculated controller. A warning is issued if the achieved γ differs from the γ calculated during design. If this warning is issued, consider using a higher-precision number type like ftype = BigFloat.

See the example folder for example usage.

hsvd(sys::AbstractStateSpace)

Return the Hankel singular values of sys, computed as the eigenvalues of QP Where Q and P are the Gramians of sys.

ispassive(P; kwargs...)

Determine if square system P is passive, i.e., $P(s) + Pᴴ(s) > 0$.

A passive system has a Nyquist curve that lies completely in the right half plane, and satisfies the following inequality (dissipation of energy)

\[\int_0^T y^T u dt > 0 ∀ T\]

The negative feedback-interconnection of two passive systems is stable and parallel connections of two passive systems as well as the inverse of a passive system are also passive. A passive controller will thus always yeild a stable feedback loop for a passive system. A series connection of two passive systems is not always passive.

See also passivityplot, passivity_index.

loop_diskmargin(P, C, args...; kwargs...)

Calculate the loop-at-a-time diskmargin for each output and input of P. See also diskmargin, sim_diskmargin. Ref: "An Introduction to Disk Margins", Peter Seiler, Andrew Packard, and Pascal Gahinet

loop_diskmargin(L, args...; kwargs...)

Calculate the loop-at-a-time diskmargin for each output of L.

See also diskmargin, sim_diskmargin. Ref: "An Introduction to Disk Margins", Peter Seiler, Andrew Packard, and Pascal Gahinet

loop_scale(L::LTISystem, w = 0)

Find the optimal diagonal scaling matrix D such that D\L(iw)*D has a minimized condition number at frequency w. Applicable to square L only. Use loop_scaling to obtain D.

loop_scaling(M0::Matrix, tol = 0.0001)

Find the optimal diagonal scaling matrix D such that D\M0*D has a minimized condition number. Applicable to square M0 only. See also structured_singular_value with option dynamic=true. Use loop_scale to find and apply the scaling to a loop-transfer function.

lqr3(P::AbstractStateSpace, Q1::AbstractMatrix, Q2::AbstractMatrix, Q3::AbstractMatrix)

Calculate the feedback gain of the discrete LQR cost function augmented with control differences

\[x^{T} Q_1 x + u^{T} Q_2 u + Δu^{T} Q_3 Δu, \quad Δu = u(k) - u(k-1)\]

makeweight(low, f_mid, high)
makeweight(low, (f_mid, gain_mid), high)

Create a weighting function that goes from gain low at zero frequency, through gain gain_mid to gain high at ∞

Arguments:

• low: A number specifying the DC gain
• mid: A number specifying the frequency at which the gain is 1, or a tuple (freq, gain). If gain_mid is not specified, the geometric mean of high and low is used.
• high: A number specifying the gain at ∞

using ControlSystemsBase, Plots
W = makeweight(10, (5,2), 1/10)
bodeplot(W)
hline!([10, 2, 1/10], l=(:black, :dash), primary=false)
vline!([5], l=(:black, :dash), primary=false)

measure(s::NamedStateSpace, names)

Return a system with specified states as measurement outputs.

sysm, T, E = modal_form(sys; C1 = false)

Bring sys to modal form.

The modal form is characterized by being tridiagonal with the real values of eigenvalues of A on the main diagonal and the complex parts on the first sub and super diagonals. T is the similarity transform applied to the system such that

sysm ≈ similarity_transform(sys, T)

If C1, then an additional convention for SISO systems is used, that the C-matrix coefficient of real eigenvalues is 1. If C1 = false, the B and C coefficients are chosen in a balanced fashion.

E is an eigen factorization of A.

See also hess_form and schur_form

muplot(sys, args...; hz=false)
muplot(LTISystem[sys1, sys2...], args...; hz=false)

Plot the structured singular values (assuming time-varying diagonal complex uncertainty) of the frequency response of the LTISystem(s). This plot is similar to sigmaplot, but scales the loop-transfer function to minimize the maximum singular value. Only applicable to square systems. A frequency vector w can be optionally provided.

If hz=true, the plot x-axis will be displayed in Hertz, the input frequency vector is still treated as rad/s.

kwargs is sent as argument to Plots.plot.

fig = mvnyquistplot(sys, w;  unit_circle=true, hz = false, kwargs...)

Create a Nyquist plot of the LTISystem. A frequency vector w must be provided.

• unit_circle: if the unit circle should be displayed

If hz=true, the hover information will be displayed in Hertz, the input frequency vector is still treated as rad/s.

kwargs is sent as argument to plot.

Example

w = 2π .* exp10.(LinRange(-2, 2, 500))
W = makeweight(0.40, 15, 3) # frequency weight for uncertain dynamics
Pn = tf(1, [1/60, 1]) |> ss # nominal plant
d = δss(1,1)                # Uncertain dynamics

Pd = Pn*(I(1) + W*d)        # weighted dynamic uncertainty on the input of Pn
Pp = rand(Pd, 200)          # sample the uncertain plant
Gcl = lft(Pd, ss(-1))       # closed loop system
structured_singular_value(Gcl) # larger than 1 => not robustly stable
unsafe_comparisons(true)
mvnyquistplot(Pp, w, points=true) # MV Nyquist plot encircles origin for some samples => not robustly stable

named_ss(sys::AbstractStateSpace, name; x, y, u)

If a single name of the system is provided, the outputs, inputs and states will be automatically named y,u,x with name as prefix.

named_ss(sys::AbstractStateSpace{T}; x, u, y)

Create a NamedStateSpace system. This kind of system uses names rather than integer indices to refer to states, inputs and outputs.

• If a single name is provided but a vector of names is expected, this name will be used as prefix followed by a numerical index.
• If no name is provided, default names (x,y,u) will be used.

Arguments:

• sys: A system to add names to.
• x: A list of symbols with names of the states.
• u: A list of symbols with names of the inputs.
• y: A list of symbols with names of the outputs.

Example

G1 = ss(1,1,1,0)
G2 = ss(1,1,1,0)
s1 = named_ss(G1, x = :x, u = :u1, y=:y1)
s2 = named_ss(G2, x = :z, u = :u2, y=:y2)

s1[:y1, :u1] # Index using symbols. Uses prefix matching if no exact match is found.

fb = feedback(s1, s2, r = :r) # 

named_ss(sys::ExtendedStateSpace;       kwargs...)
named_ss(sys::ExtendedStateSpace, name; kwargs...)

Assign names to an ExtendedStateSpace. If no specific names are provided for signals z,y,w,u and statesx, names will be generated automatically.

Arguments:

• name: Prefix to add to all automatically generated names.
• x
• u
• y
• w
• z

m, ω = ncfmargin(P, K)

Normalized coprime factor margin, defined has the inverse of

\[\begin{Vmatrix} \begin{bmatrix} I \\ K \end{bmatrix} (I + PK)^{-1} \begin{bmatrix} I & P \end{bmatrix} \end{Vmatrix}_\infty\]

A margin ≥ 0.25-0.3 is a reasonable for robustness.

If controller K stabilizes P with margin m, then K will also stabilize P̃ if nugap(P, P̃) < m.

See also extended_gangoffour, diskmargin, controller_reduction_plot.

Extended help

• Robustness with respect to coprime factor uncertainty does not necessarily imply robustness with respect to input uncertainty. Skogestad p. 96 remark 4

neglected_delay(Lmax)

Return a multiplicative weight to represent the uncertainty coming from neglecting the dynamics exp(-s*L) where L ≤ Lmax. "Multivariable Feedback Control: Analysis and Design" Ch 7.4.5

See also gain_and_delay_uncertainty and neglected_lag.

Example:

a = 10
P = ss([0 a; -a 0], I(2), [1 a; -a 1], 0) # Plant
W0 = neglected_delay(0.005) |> ss # Weight
W = I(2) + W0*I(2) * uss([δc(), δc()]) # Create a diagonal real uncertainty weighted in frequency by W0
Ps = P*W # Uncertain plant
Psamples = rand(Ps, 500) # Sample the uncertain plant for plotting
w = exp10.(LinRange(-1, 3, 300)) # Frequency vector
bodeplot(Psamples, w)

neglected_lag(τmax)

Return a multiplicative weight to represent the uncertainty coming from neglecting the dynamics 1/(s*τ + 1) where τ ≤ τmax. "Multivariable Feedback Control: Analysis and Design" Ch 7.4.5

See also gain_and_delay_uncertainty and neglected_delay.

Example:

a = 10
P = ss([0 a; -a 0], I(2), [1 a; -a 1], 0) # Plant
W0 = neglected_lag(0.05) |> ss # Weight
W = I(2) + W0*I(2) * uss([δc(), δc()]) # Create a diagonal real uncertainty weighted in frequency by W0
Ps = P*W # Uncertain plant
Psamples = rand(Ps, 100) # Sample the uncertain plant for plotting
w = exp10.(LinRange(-1, 3, 300)) # Frequency vector
sigmaplot(Psamples, w)

noise_mapping(P::ExtendedStateSpace)

Return the system from w -> y See also performance_mapping, system_mapping, noise_mapping

nu_reduction(G, g=0.1; gap = nugap(G))

Reduce the number of particles in an uncertain system G by removing all particles that are within the νgap g of the nominal system Gₙ.

Note: If G has a stochastic interpretation, i.e., the coefficients come from some distribution, this interpretation will be lost after reduction, mean values and standard deviations will not be preserved. The reduced system should instead be interpreted as preserving worst-case uncertainty.

If the gap = nugap(G) has already been precomputed, it can be supplied as an argument to avoid potentially costly recomputaiton.

nu_reduction_recursive(G, g = 0.1; gap = nugap(G), keepinds = Set{Int}(1), verbose = false)

Find a νgap cover of balls of radius g (in the νgap metric) that contain all realizations in G.

If the gap = nugap(G) has already been precomputed, it can be supplied as an argument to avoid potentially costly recomputaiton. If a manually computed gap is supplied, you must also supply keepinds=Set{Int}(index) where index is the index of the nominal system in G used to compute gap.

The returned cover Gr is of the same type as G, but with a smaller number of particles. A controller designed for Gr that achieves a ncfmargin of at least g for all realizations in Gr will stabilize all realizations in the original G. The extreme case cover where Gr = Gnominal is a single realization only can be computed by calling g = nugap(G, i) where i is the index of the nominal system in G.

Arguments:

• G: An uncertain model in the form of a StateSpace{TE, Particles} (a multi-model).
• g: The radius of the balls in the νgap cover.
• gap: An optional precomputed gap
• verbose: Print progress

nugap(G; map = map)

Compute the νgap between the nominal system Gₙ represented by the first particle index in G, and all other systems in G. Returns a Particles object with the νgap for each system in G.

See with_nominal to endow uncertain values with a nominal value, and nominal to extract the nominal value.

The value returned by this function, νᵧ is useful for robust synthesis, by designing a controller for the nominal system Gₙ, that achieves an ncfmargin of at least νᵧ is guaranteed to stabilize all realizations within G.

To speed up computation for large systems, a threaded or distributed map function can be supplied, e.g., ThreadTools.tmap or Distributed.pmap.

nugap(sys0::LTISystem, sys1::LTISystem; kwargs...)

Compute the ν-gap metric between two systems. See also ncfmargin.

For keyword arguments, see the docstring of DescriptorSystems.gnugap, reproduced below

gnugap(sys1, sys2; freq = ω, rtolinf = 0.00001, fast = true, offset = sqrt(ϵ), 
       atol = 0, atol1 = atol, atol2 = atol, rtol = n*ϵ) -> (nugapdist, fpeak)

Compute the ν-gap distance nugapdist between two descriptor systems sys1 = (A1-λE1,B1,C1,D1) and sys2 = (A2-λE2,B2,C2,D2) and the corresponding frequency fpeak (in rad/TimeUnit), where the ν-gap distance achieves its peak value.

If freq = missing, the resulting nugapdist satisfies 0 <= nugapdist <= 1. The value nugapdist = 1 results, if the winding number is different of zero in which case fpeak = [].

If freq = ω, where ω is a given vector of real frequency values, the resulting nugapdist is a vector of pointwise ν-gap distances of the dimension of ω, whose components satisfies 0 <= maximum(nugapdist) <= 1. In this case, fpeak is the frequency for which the pointwise distance achieves its peak value. All components of nugapdist are set to 1 if the winding number is different of zero in which case fpeak = [].

The stability boundary offset, β, to be used to assess the finite zeros which belong to the boundary of the stability domain can be specified via the keyword parameter offset = β. Accordingly, for a continuous-time system, these are the finite zeros having real parts within the interval [-β,β], while for a discrete-time system, these are the finite zeros having moduli within the interval [1-β,1+β]. The default value used for β is sqrt(ϵ), where ϵ is the working machine precision.

Pencil reduction algorithms are employed to compute range and coimage spaces which perform 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 atol1, atol2 and rtol, specify, respectively, the absolute tolerance for the nonzero elements of A1, A2, B1, B2, C1, C2, D1 and D2, the absolute tolerance for the nonzero elements of E1 and E2, and the relative tolerance for the nonzero elements of all above matrices. The default relative tolerance is n*ϵ, where ϵ is the working machine epsilon and n is the maximum of the orders of the systems sys1 and sys2. The keyword argument atol can be used to simultaneously set atol1 = atol, atol2 = atol.

The keyword argument rtolinf specifies the relative accuracy to be used to compute the ν-gap as the infinity norm of the relevant system according to [1]. The default value used for rtolinf is 0.00001.

Method: The evaluation of ν-gap uses the definition proposed in [1], extended to generalized LTI (descriptor) systems. The computation of winding number is based on enhancements covering zeros on the boundary of the stability domain and infinite zeros.

References:

[1] G. Vinnicombe. Uncertainty and feedback: H∞ loop-shaping and the ν-gap metric. Imperial College Press, London, 2001.

partition(P::AbstractStateSpace, nw::Int, nz::Int)

partition(P::AbstractStateSpace; u, y, w=!u, z=!y)

Partition P into an ExtendedStateSpace.

• u indicates the indices of the controllable inputs.
• y indicates the indices of the measurable outputs.
• w indicates the indices of the disturbance inputs (uncontrollable), by default w is the complement of u.
• z indicates the indices of the performance outputs (not neccesarily measurable), by default z is the complement of y.

passivity_index(P; kwargs...)

Return

\[γ = \begin{Vmatrix} (I-P)(I+P)^{-1} \end{Vmatrix}_∞\]

If $γ ≤ 1$, the system is passive. If the system has unstable zeros, $γ = ∞$

The negative feedback interconnection of two systems with passivity indices γ₁ and γ₂ is stable if $γ₁γ₂ < 1$.

A passive system has a Nyquist curve that lies completely in the right half plane, and satisfies the following inequality (dissipation of energy)

\[\int_0^T y^T u dt > 0 ∀ T\]

The negative feedback-interconnection of two passive systems is stable and parallel connections of two passive systems as well as the inverse of a passive system are also passive. A passive controller will thus always yeild a stable feedback loop for a passive system. A series connection of two passive systems is not always passive.

See also ispassive, passivityplot.

passivityplot(sys, args...; hz=false)
passivityplot(LTISystem[sys1, sys2...], args...; hz=false)

Plot the passivity index of a LTISystem(s). The system is passive for frequencies where the index is < 0.

A frequency vector w can be optionally provided.

If hz=true, the plot x-axis will be displayed in Hertz, the input frequency vector is still treated as rad/s.

kwargs is sent as argument to Plots.plot.

See passivity_index for additional details. See also ispassive, passivity_index.

performance_mapping(P::ExtendedStateSpace)

Return the system from w -> z See also performance_mapping, system_mapping, noise_mapping

robstab(M0::UncertainSS, w=exp10.(LinRange(-3, 3, 1500)); kwargs...)

Return the robust stability margin of an uncertain model, defined as the inverse of the structured singular value. Currently, only diagonal complex perturbations supported.

sysm, T, SF = schur_form(sys)

Bring sys to Schur form.

The Schur form is characterized by A being Schur with the real values of eigenvalues of A on the main diagonal. T is the similarity transform applied to the system such that

sysm ≈ similarity_transform(sys, T)

SF is the Schur-factorization of A.

See also modal_form and hess_form

show_construction([io::IO,] sys::LTISystem; name = "temp", letb = true)

Print code to io that reconstructs sys.

• letb: If true, the code is surrounded by a let block.

julia> sys = ss(tf(1, [1, 1]))
StateSpace{Continuous, Float64}
A = 
 -1.0
B = 
 1.0
C = 
 1.0
D = 
 0.0

Continuous-time state-space model

julia> show_construction(sys, name="Jörgen")
Jörgen = let
    JörgenA = [-1.0;;]
    JörgenB = [1.0;;]
    JörgenC = [1.0;;]
    JörgenD = [0.0;;]
    ss(JörgenA, JörgenB, JörgenC, JörgenD)
end

sim_diskmargin(P::LTISystem, C::LTISystem, σ::Real = 0)

Simultaneuous diskmargin at both outputs and inputs of P. Ref: "An Introduction to Disk Margins", Peter Seiler, Andrew Packard, and Pascal Gahinet https://arxiv.org/abs/2003.04771 See also ncfmargin.

sim_diskmargin(L, σ::Real = 0, l=1e-3, u=1e3)

Return the smallest simultaneous diskmargin over the grid l:u See also ncfmargin.

sim_diskmargin(L, σ::Real, w::AbstractVector)
sim_diskmargin(L, σ::Real = 0)

Simultaneuous diskmargin at the outputs of L. Users should consider using diskmargin.

specificationplot([S,CS,T], [WS,WU,WT])

This function visualizes the control synthesis using the hinfsynthesize with the three weighting functions $W_S(s), W_U(s), W_T(s)$ inverted and scaled by γ, against the corresponding transfer functions $S(s), C(s)S(s), T(s)$, to verify visually that the specifications are met. This may be run using both MIMO and SISO systems.

Keyword args

• wint: (-3, 5) frequency range (log10)
• wnum: 201 number of frequency points
• hz: true
• nsigma: typemax(Int) number of singular values to show
• s_labels: `[ "σ(S)", "σ(CS)", "σ(T)",

]`

• w_labels: `[ "γ σ(Wₛ⁻¹)", "γ σ(Wᵤ⁻¹)", "γ σ(Wₜ⁻¹)",

]`

• colors: [:blue, :red, :green] colors for $S$, $CS$ and $T$

splitter(u::Symbol, n::Int, timeevol = Continuous())

Return a named system that splits an input signal into n signals. This is useful when an external signal entering a block diagram is to be connected to multiple inputs. See the tutorial https://juliacontrol.github.io/RobustAndOptimalControl.jl/dev/hinf_connection/ for example usage. An alternative way of connecting an external input to several input ports with the same name is to pass connect(..., unique=false).

Arguments:

• u: Named of the signal to split
• n: Number of splits

ss2particles(G::Vector{<:AbstractStateSpace})

Converts a vector of state space models to a single state space model with coefficient type MonteCarloMeasurements.Particles.

See also sys_from_particles.

A, B1, B2, C1, C2, D11, D12, D21, D22 = ssdata_e(sys)

stab, unstab = stab_unstab(sys; kwargs...)

Decompose sys into sys = stab + unstab where stab contains all stable poles and unstab contains unstable poles. See gsdec(sys; job = "finite", prescale, smarg, fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol = nϵ) -> (sys1, sys2)

Compute for the descriptor system sys = (A-λE,B,C,D) with the transfer function matrix G(λ), the additive spectral decomposition G(λ) = G1(λ) + G2(λ) such that G1(λ), the transfer function matrix of the descriptor system sys1 = (A1-λE1,B1,C1,D1), has only poles in a certain domain of interest Cg of the complex plane and G2(λ), the transfer function matrix of the descriptor system sys2 = (A2-λE2,B2,C2,0), has only poles outside of Cg.

If prescale = true, a preliminary balancing of the descriptor system pair (A,E) is performed. The default setting is prescale = MatrixPencils.balqual(sys.A,sys.E) > 10000, where the function pbalqual from the MatrixPencils package evaluates the scaling quality of the linear pencil A-λE.

The keyword argument smarg, if provided, specifies the stability margin for the stable eigenvalues of A-λE, such that, in the continuous-time case, the stable eigenvalues have real parts less than or equal to smarg, and in the discrete-time case, the stable eigenvalues have moduli less than or equal to smarg. If smarg = missing, the used default values are: smarg = -sqrt(ϵ), for a continuous-time system, and smarg = 1-sqrt(ϵ), for a discrete-time system), where ϵ is the machine precision of the working accuracy.

The keyword argument job, in conjunction with smarg, defines the domain of interest Cg, as follows:

for job = "finite", Cg is the whole complex plane without the point at infinity, and sys1 has only finite poles and sys2 has only infinite poles (default); the resulting A2 is nonsingular and upper triangular, while the resulting E2 is nilpotent and upper triangular;

for job = "infinite", Cg is the point at infinity, and sys1 has only infinite poles and sys2 has only finite poles and is the strictly proper part of sys; the resulting A1 is nonsingular and upper triangular, while the resulting E1 is nilpotent and upper triangular;

for job = "stable", Cg is the stability domain of eigenvalues defined by smarg, and sys1 has only stable poles and sys2 has only unstable and infinite poles; the resulting pairs (A1,E1) and (A2,E2) are in generalized Schur form with E1 upper triangular and nonsingular and E2 upper triangular;

for job = "unstable", Cg is the complement of the stability domain of the eigenvalues defined by smarg, and sys1 has only unstable and infinite poles and sys2 has only stable poles; the resulting pairs (A1,E1) and (A2,E2) are in generalized Schur form with E1 upper triangular and E2 upper triangular and nonsingular.

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

The separation of the finite and infinite eigenvalues is performed using rank decisions based on rank revealing QR-decompositions with column pivoting if fast = true or the more reliable SVD-decompositions if fast = false. for keyword arguments (argument job is set to "stable" in this function).

static_gain_compensation(l::LQGProblem, L = lqr(l))
static_gain_compensation(A, B, C, D, L)

Find $L_r$ such that

dcgain(closedloop(G)*Lr) ≈ I

structured_singular_value(M0::UncertainSS, [w::AbstractVector]; kwargs...)

• w: Frequency vector, if none is provided, the maximum μ over a grid 1e-3 : 1e3 will be returned.

An example is provided in the documentation

μ = structured_singular_value(M; tol=1e-4, scalings=false, dynamic=false)

Compute (an upper bound of) the structured singular value μ for diagonal Δ of complex perturbations (other structures of Δ are not yet supported). M is assumed to be an (n × n × N_freq) array or a matrix.

We currently don't have any methods to compute a lower bound, but if all perturbations are complex the spectral radius ρ(M) is always a lower bound (usually not a good one).

If scalings = true, return also a n × nf matrix Dm with the diagonal scalings D such that

D = Diagonal(Dm[:, i])
σ̄(D\M[:,:,i]*D)

is minimized.

If dynamic = true, the perturbations are assumed to be time-varying Δ(t). In this case, the same scaling is used for all frequencies and the returned D if scalings=true is a vector d such that D = Diagonal(d).

sumblock(ex::String; Ts = 0, n = 1)

Create a summation node (named statespace system) that sums (or subtracts) vectors of length n.

Arguments:

• Ts: Sample time
• n: The length of the input and output vectors. Set n=1 for scalars.

When using sumblock to form block diagrams, note how the system returned from sumblock has input names corresponding to the right-hand side of the expression and output names corresponding to the variable on the left-hand side. You will thus typically list connections like :y => :y in the connection list to the connect function. See the tutorials

• https://juliacontrol.github.io/RobustAndOptimalControl.jl/dev/hinf_connection/
• https://juliacontrol.github.io/RobustAndOptimalControl.jl/dev/api/#RobustAndOptimalControl.connect-Tuple{Any}

for example usage

Examples:

julia> sumblock("uP = vf + yL")
NamedStateSpace{Continuous, Int64}
D = 
 1  1

With state  names: 
     input  names: vf yL
     output names: uP


julia> sumblock("x_diff = xr - xh"; n=3)
NamedStateSpace{Continuous, Int64}
D = 
 1  0  0  -1   0   0
 0  1  0   0  -1   0
 0  0  1   0   0  -1

With state  names: 
     input  names: xr1 xr2 xr3 xh1 xh2 xh3
     output names: x_diff1 x_diff2 x_diff3
     

julia> sumblock("a = b + c - d")
NamedStateSpace{Continuous, Int64}
D = 
 1  1  -1

With state  names: 
     input  names: b c d
     output names: a

sys_from_particles(P, i)
sys_from_particles(P)

Return the ith system from a system P with Particles coefficients.

If called without an index, return a vector of systems, one for each possibly i.

This function is used to convert from an uncertain representation using Particles to a "multi-model" representation using multiple StateSpace models.

See also ss2particles and MonteCarloMeasurements.nominal.

system_mapping(P::ExtendedStateSpace)

Return the system from u -> y See also performance_mapping, system_mapping, noise_mapping

uss(d::AbstractVector{<:δ}, Ts = nothing)

Create a diagonal uncertain statespace object with the uncertain elements d on the diagonal.

uss(D::AbstractArray, Δ, Ts = nothing)

If only a single D matrix is provided, it's treated as D11 if Δ is given, and as D22 if no Δ is provided.

uss(D11, D12, D21, D22, Δ, Ts = nothing)

Create an uncertain statespace object with only gin matrices.

uss(d::δ{C, F}, Ts = nothing)

Convert a δ object to an UncertainSS

vec2sys(v::AbstractArray, ny::Int, nu::Int, ts = nothing)

Create a statespace system from the parameters

v = vec(sys) = [vec(sys.A); vec(sys.B); vec(sys.C); vec(sys.D)]

Use vec(sys) to create v.

This can be useful in order to convert to and from vectors for, e.g., optimization.

julia> sys  = ss(tf(1, [1, 1]))
StateSpace{Continuous, Float64}
A = 
 -1.0
B = 
 1.0
C = 
 1.0
D = 
 0.0

Continuous-time state-space model

julia> v    = vec(sys)
4-element Vector{Float64}:
 -1.0
  1.0
  1.0
  0.0

julia> sys2 = vec2sys(v, sys.ny, sys.nu)
StateSpace{Continuous, Float64}
A = 
 -1.0
B = 
 1.0
C = 
 1.0
D = 
 0.0

Continuous-time state-space model

δc(val::Complex = complex(0.0), radius::Real = 1.0, name)

Create a complex, uncertain parameter. If no name is given, a boring name will be generated automatically.

δr(val::Real = 0.0, radius::Real = 1.0, name)

Create a real, uncertain parameter. If no name is given, a boring name will be generated automatically.