`ControlSystemsBase.are`

`ControlSystemsBase.are`

`ControlSystemsBase.balance`

`ControlSystemsBase.balance_statespace`

`ControlSystemsBase.balreal`

`ControlSystemsBase.baltrunc`

`ControlSystemsBase.controllability`

`ControlSystemsBase.covar`

`ControlSystemsBase.ctrb`

`ControlSystemsBase.damp`

`ControlSystemsBase.dampreport`

`ControlSystemsBase.dcgain`

`ControlSystemsBase.delaymargin`

`ControlSystemsBase.gangoffour`

`ControlSystemsBase.gangofseven`

`ControlSystemsBase.gram`

`ControlSystemsBase.grampd`

`ControlSystemsBase.hinfnorm`

`ControlSystemsBase.innovation_form`

`ControlSystemsBase.innovation_form`

`ControlSystemsBase.linfnorm`

`ControlSystemsBase.margin`

`ControlSystemsBase.markovparam`

`ControlSystemsBase.observability`

`ControlSystemsBase.observer_controller`

`ControlSystemsBase.observer_filter`

`ControlSystemsBase.observer_predictor`

`ControlSystemsBase.obsv`

`ControlSystemsBase.plyap`

`ControlSystemsBase.poles`

`ControlSystemsBase.reduce_sys`

`ControlSystemsBase.relative_gain_array`

`ControlSystemsBase.relative_gain_array`

`ControlSystemsBase.similarity_transform`

`ControlSystemsBase.time_scale`

`ControlSystemsBase.tzeros`

`ControlSystemsBase.zpkdata`

`LinearAlgebra.lyap`

`LinearAlgebra.norm`

For robust analysis, see RobustAndOptimalControl.jl.

# Analysis

`ControlSystemsBase.damp`

— Method`Wn, zeta, ps = damp(sys)`

Compute the natural frequencies, `Wn`

, and damping ratios, `zeta`

, of the poles, `ps`

, of `sys`

`ControlSystemsBase.dampreport`

— Method`dampreport(sys)`

Display a report of the poles, damping ratio, natural frequency, and time constant of the system `sys`

`ControlSystemsBase.dcgain`

— Function`dcgain(sys, ϵ=0)`

Compute the dcgain of system `sys`

.

equal to G(0) for continuous-time systems and G(1) for discrete-time systems.

`ϵ`

can be provided to evaluate the dcgain with a small perturbation into the stability region of the complex plane.

`ControlSystemsBase.delaymargin`

— Method`dₘ = delaymargin(G::LTISystem)`

Return the delay margin, dₘ. For discrete-time systems, the delay margin is normalized by the sample time, i.e., the value represents the margin in number of sample times. Only supports SISO systems.

`ControlSystemsBase.gangoffour`

— Method```
S, PS, CS, T = gangoffour(P, C; minimal=true)
gangoffour(P::AbstractVector, C::AbstractVector; minimal=true)
```

Given a transfer function describing the plant `P`

and a transfer function describing the controller `C`

, computes the four transfer functions in the Gang-of-Four.

`S = 1/(1+PC)`

Sensitivity function`PS = (1+PC)\P`

Load disturbance to measurement signal`CS = (1+PC)\C`

Measurement noise to control signal`T = PC/(1+PC)`

Complementary sensitivity function

If `minimal=true`

, `minreal`

will be applied to all transfer functions.

`ControlSystemsBase.gangofseven`

— Method`S, PS, CS, T, RY, RU, RE = gangofseven(P,C,F)`

Given transfer functions describing the Plant `P`

, the controller `C`

and a feed forward block `F`

, computes the four transfer functions in the Gang-of-Four and the transferfunctions corresponding to the feed forward.

`S = 1/(1+PC)`

Sensitivity function`PS = P/(1+PC)`

`CS = C/(1+PC)`

`T = PC/(1+PC)`

Complementary sensitivity function`RY = PCF/(1+PC)`

`RU = CF/(1+P*C)`

`RE = F/(1+P*C)`

`ControlSystemsBase.margin`

— Method`wgm, gm, wpm, pm = margin(sys::LTISystem, w::Vector; full=false, allMargins=false)`

returns frequencies for gain margins, gain margins, frequencies for phase margins, phase margins

If `!allMargins`

, return only the smallest margin

If `full`

return also `fullPhase`

See also `delaymargin`

and `RobustAndOptimalControl.diskmargin`

`ControlSystemsBase.markovparam`

— Method`markovparam(sys, n)`

Compute the `n`

th markov parameter of discrete-time state-space system `sys`

. This is defined as the following:

`h(0) = D`

`h(n) = C*A^(n-1)*B`

`ControlSystemsBase.poles`

— Method`poles(sys)`

Compute the poles of system `sys`

.

`ControlSystemsBase.reduce_sys`

— Method`reduce_sys(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix, meps::AbstractFloat)`

Implements REDUCE in the Emami-Naeini & Van Dooren paper. Returns transformed A, B, C, D matrices. These are empty if there are no zeros.

`ControlSystemsBase.relative_gain_array`

— Method`relative_gain_array(A::AbstractMatrix; tol = 1.0e-15)`

Reference: "On the Relative Gain Array (RGA) with Singular and Rectangular Matrices" Jeffrey Uhlmann https://arxiv.org/pdf/1805.10312.pdf

`ControlSystemsBase.relative_gain_array`

— Method```
relative_gain_array(G, w::AbstractVector)
relative_gain_array(G, w::Number)
```

Calculate the relative gain array of `G`

at frequencies `w`

. G(iω) .* pinv(tranpose(G(iω)))

The RGA can be used to find input-output pairings for MIMO control using individually tuned loops. Pair the inputs and outputs such that the RGA(ωc) at the crossover frequency becomes as close to diagonal as possible. Avoid pairings such that RGA(0) contains negative diagonal elements.

- The sum of the absolute values of the entries in the RGA is a good measure of the "true condition number" of G, the best condition number that can be achieved by input/output scaling of
`G`

, -Glad, Ljung. - The RGA is invariant to input/output scaling of
`G`

. - If the RGA contains large entries, the system may be sensitive to model errors, -Skogestad, "Multivariable Feedback Control: Analysis and Design":
- Uncertainty in the input channels (diagonal input uncertainty). Plants with

- Element uncertainty. Large RGA-elements imply sensitivity to element-by-element uncertainty.

The relative gain array is computed using the The unit-consistent (UC) generalized inverse Reference: "On the Relative Gain Array (RGA) with Singular and Rectangular Matrices" Jeffrey Uhlmann https://arxiv.org/pdf/1805.10312.pdf

`ControlSystemsBase.tzeros`

— Method`tzeros(sys)`

Compute the invariant zeros of the system `sys`

. If `sys`

is a minimal realization, these are also the transmission zeros.

`ControlSystemsBase.zpkdata`

— Method`z, p, k = zpkdata(sys)`

Compute the zeros, poles, and gains of system `sys`

.

**Returns**

`z`

: Matrix{Vector{ComplexF64}}, (ny × nu)`p`

: Matrix{Vector{ComplexF64}}, (ny × nu)`k`

: Matrix{Float64}, (ny × nu)

`ControlSystemsBase.are`

— Method`are(::Continuous, A, B, Q, R)`

Compute 'X', the solution to the continuous-time algebraic Riccati equation, defined as A'X + XA - (XB)R^-1(B'X) + Q = 0, where R is non-singular.

In an LQR problem, `Q`

is associated with the state penalty $x'Qx$ while `R`

is associated with the control penalty $u'Ru$. See `lqr`

for more details.

Uses `MatrixEquations.arec`

. For keyword arguments, see the docstring of `ControlSystemsBase.MatrixEquations.arec`

, note that they define the input arguments in a different order.

`ControlSystemsBase.are`

— Method`are(::Discrete, A, B, Q, R; kwargs...)`

Compute `X`

, the solution to the discrete-time algebraic Riccati equation, defined as A'XA - X - (A'XB)(B'XB + R)^-1(B'XA) + Q = 0, where Q>=0 and R>0

In an LQR problem, `Q`

is associated with the state penalty $x'Qx$ while `R`

is associated with the control penalty $u'Ru$. See `lqr`

for more details.

Uses `MatrixEquations.ared`

. For keyword arguments, see the docstring of `ControlSystemsBase.MatrixEquations.ared`

, note that they define the input arguments in a different order.

`ControlSystemsBase.balance`

— Function`S, P, B = balance(A[, perm=true])`

Compute a similarity transform `T = S*P`

resulting in `B = T\A*T`

such that the row and column norms of `B`

are approximately equivalent. If `perm=false`

, the transformation will only scale `A`

using diagonal `S`

, and not permute `A`

(i.e., set `P=I`

).

`ControlSystemsBase.balreal`

— Method`sysr, G, T = balreal(sys::StateSpace)`

Calculates a balanced realization of the system sys, such that the observability and reachability gramians of the balanced system are equal and diagonal `diagm(G)`

. `T`

is the similarity transform between the old state `x`

and the new state `z`

such that `z = Tx`

.

Reference: Varga A., Balancing-free square-root algorithm for computing singular perturbation approximations.

`ControlSystemsBase.baltrunc`

— Method`sysr, G, T = baltrunc(sys::StateSpace; atol = √ϵ, rtol=1e-3, n = nothing, residual = false)`

Reduces the state dimension by calculating a balanced realization of the system sys, such that the observability and reachability gramians of the balanced system are equal and diagonal `diagm(G)`

, and truncating it to order `n`

. If `n`

is not provided, it's chosen such that all states corresponding to singular values less than `atol`

and less that `rtol σmax`

are removed.

`T`

is the projection matrix between the old state `x`

and the newstate `z`

such that `z = Tx`

. `T`

will in general be a non-square matrix.

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.

See also `gram`

, `balreal`

Glad, Ljung, Reglerteori: Flervariabla och Olinjära metoder.

For more advanced model reduction, see RobustAndOptimalControl.jl - Model Reduction.

**Extended help**

Note: Gramian computations are sensitive to input-output scaling. For the result of a numerical balancing, gramian computation or truncation of MIMO systems to be meaningful, the inputs and outputs of the system must thus be scaled in a meaningful way. A common (but not the only) approach is:

- The outputs are scaled such that the maximum allowed control error, the maximum expected reference variation, or the maximum expected variation, is unity.
- The input variables are scaled to have magnitude one. This is done by dividing each variable by its maximum expected or allowed change, i.e., $u_{scaled} = u / u_{max}$

Without such scaling, the result of balancing will depend on the units used to measure the input and output signals, e.g., a change of unit for one output from meter to millimeter will make this output 1000x more important.

`ControlSystemsBase.controllability`

— Method```
controllability(A, B; atol, rtol)
controllability(sys; atol, rtol)
```

Check for controllability of the pair `(A, B)`

or `sys`

using the PHB test.

The return value contains the field `iscontrollable`

which is `true`

if the rank condition is met at all eigenvalues of `A`

, and `false`

otherwise. The returned structure also contains the rank and smallest singular value at each individual eigenvalue of `A`

in the fields `ranks`

and `sigma_min`

.

Technically, this function checks for controllability from the origin, also called reachability.

`ControlSystemsBase.covar`

— Method`P = covar(sys, W)`

Calculate the stationary covariance `P = E[y(t)y(t)']`

of the output `y`

of a `StateSpace`

model `sys`

driven by white Gaussian noise `w`

with covariance `E[w(t)w(τ)]=W*δ(t-τ)`

(δ is the Dirac delta).

Remark: If `sys`

is unstable then the resulting covariance is a matrix of `Inf`

s. Entries corresponding to direct feedthrough (D*W*D' .!= 0) will equal `Inf`

for continuous-time systems.

`ControlSystemsBase.ctrb`

— Method```
ctrb(A, B)
ctrb(sys)
```

Compute the controllability matrix for the system described by `(A, B)`

or `sys`

.

Note that checking for controllability by computing the rank from `ctrb`

is not the most numerically accurate way, a better method is checking if `gram(sys, :c)`

is positive definite or to call the function `controllability`

.

`ControlSystemsBase.gram`

— Method`gram(sys, opt; kwargs...)`

Compute the grammian of system `sys`

. If `opt`

is `:c`

, computes the controllability grammian. If `opt`

is `:o`

, computes the observability grammian.

See also `grampd`

For keyword arguments, see `grampd`

.

**Extended help**

Note: Gramian computations are sensitive to input-output scaling. For the result of a numerical balancing, gramian computation or truncation of MIMO systems to be meaningful, the inputs and outputs of the system must thus be scaled in a meaningful way. A common (but not the only) approach is:

- The outputs are scaled such that the maximum allowed control error, the maximum expected reference variation, or the maximum expected variation, is unity.
- The input variables are scaled to have magnitude one. This is done by dividing each variable by its maximum expected or allowed change, i.e., $u_{scaled} = u / u_{max}$

Without such scaling, the result of balancing will depend on the units used to measure the input and output signals, e.g., a change of unit for one output from meter to millimeter will make this output 1000x more important.

`ControlSystemsBase.grampd`

— Method`U = grampd(sys, opt; kwargs...)`

Return a Cholesky factor `U`

of the grammian of system `sys`

. If `opt`

is `:c`

, computes the controllability grammian `G = U*U'`

. If `opt`

is `:o`

, computes the observability grammian `G = U'U`

.

Obtain a `Cholesky`

object by `Cholesky(U)`

for observability grammian

Uses `MatrixEquations.plyapc/plyapd`

. For keyword arguments, see the docstring of `ControlSystemsBase.MatrixEquations.plyapc/plyapd`

`ControlSystemsBase.hinfnorm`

— Method`Ninf, ω_peak = hinfnorm(sys; tol=1e-6)`

Compute the H∞ norm `Ninf`

of the LTI system `sys`

, together with a frequency `ω_peak`

at which the gain Ninf is achieved.

`Ninf := sup_ω σ_max[sys(iω)]`

if `G`

is stable (σ_max = largest singular value) := `Inf' if`

G` is unstable

`tol`

is an optional keyword argument for the desired relative accuracy for the computed H∞ norm (not an absolute certificate).

`sys`

is first converted to a state space model if needed.

The continuous-time L∞ norm computation implements the 'two-step algorithm' in:**N.A. Bruinsma and M. Steinbuch**, 'A fast algorithm to compute the H∞-norm of a transfer function matrix', Systems and Control Letters (1990), pp. 287-293.

For the discrete-time version, see:**P. Bongers, O. Bosgra, M. Steinbuch**, 'L∞-norm calculation for generalized state space systems in continuous and discrete time', American Control Conference, 1991.

See also `linfnorm`

.

`ControlSystemsBase.innovation_form`

— Method```
sysi = innovation_form(sys, R1, R2[, R12])
sysi = innovation_form(sys; sysw=I, syse=I, R1=I, R2=I)
```

Takes a system

```
x' = Ax + Bu + w ~ R1
y = Cx + Du + e ~ R2
```

and returns the system

```
x' = Ax + Kv
y = Cx + v
```

where `v`

is the innovation sequence.

If `sysw`

(`syse`

) is given, the covariance resulting in filtering noise with `R1`

(`R2`

) through `sysw`

(`syse`

) is used as covariance.

See Stochastic Control, Chapter 4, Åström

`ControlSystemsBase.innovation_form`

— Method`sysi = innovation_form(sys, K)`

Takes a system

```
x' = Ax + Bu + Kv
y = Cx + Du + v
```

and returns the system

```
x' = Ax + Kv
y = Cx + v
```

where `v`

is the innovation sequence.

See Stochastic Control, Chapter 4, Åström

`ControlSystemsBase.linfnorm`

— Method`Ninf, ω_peak = linfnorm(sys; tol=1e-6)`

Compute the L∞ norm `Ninf`

of the LTI system `sys`

, together with a frequency `ω_peak`

at which the gain `Ninf`

is achieved.

`Ninf := sup_ω σ_max[sys(iω)]`

(σ_max denotes the largest singular value)

`tol`

is an optional keyword argument representing the desired relative accuracy for the computed L∞ norm (this is not an absolute certificate however).

`sys`

is first converted to a state space model if needed.

The continuous-time L∞ norm computation implements the 'two-step algorithm' in:**N.A. Bruinsma and M. Steinbuch**, 'A fast algorithm to compute the H∞-norm of a transfer function matrix', Systems and Control Letters (1990), pp. 287-293.

For the discrete-time version, see:**P. Bongers, O. Bosgra, M. Steinbuch**, 'L∞-norm calculation for generalized state space systems in continuous and discrete time', American Control Conference, 1991.

See also `hinfnorm`

.

`ControlSystemsBase.observability`

— Method`observability(A, C; atol, rtol)`

Check for observability of the pair `(A, C)`

or `sys`

using the PHB test.

The return value contains the field `isobservable`

which is `true`

if the rank condition is met at all eigenvalues of `A`

, and `false`

otherwise. The returned structure also contains the rank and smallest singular value at each individual eigenvalue of `A`

in the fields `ranks`

and `sigma_min`

.

`ControlSystemsBase.observer_controller`

— Method`cont = observer_controller(sys, L::AbstractMatrix, K::AbstractMatrix; direct=false)`

**If direct = false**

Return the observer_controller `cont`

that is given by `ss(A - B*L - K*C + K*D*L, K, L, 0)`

such that `feedback(sys, cont)`

produces a closed-loop system with eigenvalues given by `A-KC`

and `A-BL`

.

This controller does not have a direct term, and corresponds to state feedback operating on state estimated by `observer_predictor`

. Use this form if the computed control signal is applied at the next sampling instant, or with an otherwise large delay in relation to the measurement fed into the controller.

Ref: "Computer-Controlled Systems" Eq 4.37

**If direct = true**

Return the observer controller `cont`

that is given by `ss((I-KC)(A-BL), (I-KC)(A-BL)K, L, LK)`

such that `feedback(sys, cont)`

produces a closed-loop system with eigenvalues given by `A-BL`

and `A-BL-KC`

. This controller has a direct term, and corresponds to state feedback operating on state estimated by `observer_filter`

. Use this form if the computed control signal is applied immediately after receiveing a measurement. This version typically has better performance than the one without a direct term.

Ref: Ref: "Computer-Controlled Systems" pp 140 and "Computer-Controlled Systems" pp 162 prob 4.7

**Arguments:**

`sys`

: Model of system`L`

: State-feedback gain`u = -Lx`

`K`

: Observer gain

See also `observer_predictor`

and `innovation_form`

.

`ControlSystemsBase.observer_filter`

— Method`observer_filter(sys, K; output_state = false)`

Return the observer filter

\[\begin{aligned} x̂(k|k) &= (I - KC)Ax̂(k-1|k-1) + (I - KC)Bu(k-1) + Ky(k) \\ \end{aligned}\]

with the input equation `[(I - KC)B K] * [u(k-1); y(k)]`

.

Note the time indices in the equations, the filter assumes that the user passes the *current* $y(k)$, but the *past* $u(k-1)$, that is, this filter is used to estimate the state *before* the current control input has been applied. This causes a state-feedback controller acting on the estimate produced by this observer to have a direct term.

This is similar to `observer_predictor`

, but in contrast to the predictor, the filter output depends on the current measurement, whereas the predictor output only depend on past measurements.

The observer filter is equivalent to the `observer_predictor`

for continuous-time systems.

Ref: "Computer-Controlled Systems" Eq 4.32

`ControlSystemsBase.observer_predictor`

— Method```
observer_predictor(sys::AbstractStateSpace, K; h::Int = 1, output_state = false)
observer_predictor(sys::AbstractStateSpace, R1, R2[, R12]; output_state = false)
```

If `sys`

is continuous, return the observer predictor system

\[\begin{aligned} x̂' &= (A - KC)x̂ + (B-KD)u + Ky \\ ŷ &= Cx + Du \end{aligned}\]

with the input equation `[B-KD K] * [u; y]`

If `sys`

is discrete, the prediction horizon `h`

may be specified, in which case measurements up to and including time `t-h`

and inputs up to and including time `t`

are used to predict `y(t)`

.

If covariance matrices `R1, R2`

are given, the kalman gain `K`

is calculated using `kalman`

.

If `output_state`

is true, the output is the state estimate `x̂`

instead of the output estimate `ŷ`

.

See also `innovation_form`

, `observer_controller`

and `observer_filter`

.

`ControlSystemsBase.obsv`

— Function```
obsv(A, C, n=size(A,1))
obsv(sys, n=sys.nx)
```

Compute the observability matrix with `n`

rows for the system described by `(A, C)`

or `sys`

. Providing the optional `n > sys.nx`

returns an extended observability matrix.

Note that checking for observability by computing the rank from `obsv`

is not the most numerically accurate way, a better method is checking if `gram(sys, :o)`

is positive definite or to call the function `observability`

.

`ControlSystemsBase.plyap`

— Method`Xc = plyap(sys::AbstractStateSpace, Ql; kwargs...)`

Lyapunov solver that takes the `L`

Cholesky factor of `Q`

and returns a triangular matrix `Xc`

such that `Xc*Xc' = X`

.

`ControlSystemsBase.similarity_transform`

— Method`syst = similarity_transform(sys, T; unitary=false)`

Perform a similarity transform `T : Tx̃ = x`

on `sys`

such that

```
Ã = T⁻¹AT
B̃ = T⁻¹ B
C̃ = CT
D̃ = D
```

If `unitary=true`

, `T`

is assumed unitary and the matrix adjoint is used instead of the inverse. See also `balance_statespace`

.

`ControlSystemsBase.time_scale`

— Method```
time_scale(sys::AbstractStateSpace{Continuous}, a; balanced = false)
time_scale(G::TransferFunction{Continuous}, a; balanced = true)
```

Rescale the time axis (change time unit) of `sys`

.

For systems where the dominant time constants are very far from 1, e.g., in electronics, rescaling the time axis may be beneficial for numerical performance, in particular for continuous-time simulations.

Scaling of time for a function $f(t)$ with Laplace transform $F(s)$ can be stated as

\[f(at) \leftrightarrow \dfrac{1}{a} F\big(\dfrac{s}{a}\big)\]

The keyword argument `balanced`

indicates whether or not to apply a balanced scaling on the `B`

and `C`

matrices. For statespace systems, this defaults to false since it changes the state representation, only `B`

will be scaled. For transfer functions, it defaults to true.

**Example:**

The following example show how a system with a time constant on the order of one micro-second is rescaled such that the time constant becomes 1, i.e., the time unit is changed from seconds to micro-seconds.

```
Gs = tf(1, [1e-6, 1]) # micro-second time scale modeled in seconds
Gms = time_scale(Gs, 1e-6) # Change to micro-second time scale
Gms == tf(1, [1, 1]) # Gms now has micro-seconds as time unit
```

The next example illustrates how the time axis of a time-domain simulation changes by time scaling

```
t = 0:0.1:50 # original time axis
a = 10 # Scaling factor
sys1 = ssrand(1,1,5)
res1 = step(sys1, t) # Perform original simulation
sys2 = time_scale(sys, a) # Scale time
res2 = step(sys2, t ./ a) # Simulate on scaled time axis, note the `1/a`
isapprox(res1.y, res2.y, rtol=1e-3, atol=1e-3)
```

`LinearAlgebra.lyap`

— Method`lyap(A, Q; kwargs...)`

Compute the solution `X`

to the discrete Lyapunov equation `AXA' - X + Q = 0`

.

Uses `MatrixEquations.lyapc / MatrixEquations.lyapd`

. For keyword arguments, see the docstring of `ControlSystemsBase.MatrixEquations.lyapc / ControlSystemsBase.MatrixEquations.lyapd`

`LinearAlgebra.norm`

— Function`norm(sys, p=2; tol=1e-6)`

`norm(sys)`

or `norm(sys,2)`

computes the H2 norm of the LTI system `sys`

.

`norm(sys, Inf)`

computes the H∞ norm of the LTI system `sys`

. The H∞ norm is the same as the L∞ for stable systems, and Inf for unstable systems. If the peak gain frequency is required as well, use the function `hinfnorm`

instead. See `hinfnorm`

for further documentation.

`tol`

is an optional keyword argument, used only for the computation of L∞ norms. It represents the desired relative accuracy for the computed L∞ norm (this is not an absolute certificate however).

`sys`

is first converted to a `StateSpace`

model if needed.

`ControlSystemsBase.balance_statespace`

— Function```
A, B, C, T = balance_statespace{S}(A::Matrix{S}, B::Matrix{S}, C::Matrix{S}, perm::Bool=false)
sys, T = balance_statespace(sys::StateSpace, perm::Bool=false)
```

Computes a balancing transformation `T`

that attempts to scale the system so that the row and column norms of [T*A/T T*B; C/T 0] are approximately equal. If `perm=true`

, the states in `A`

are allowed to be reordered.

The inverse of `sysb, T = balance_statespace(sys)`

is given by `similarity_transform(sysb, T)`

This is not the same as finding a balanced realization with equal and diagonal observability and reachability gramians, see `balreal`

## Videos

Basic usage of robustness analysis with JuliaControl