Exported functions and types
In addition to the documentation for this package, we encourage the reader to explore the documentation for ControlSystems.jl and RobustAndOptimalControl.jl that contains functionality and types for basic control analysis and design, as well as the documentation of ModelingToolkit for modeling and simulation.
Index
ControlSystems.Simulator
ControlSystems.Simulator
ControlSystemsBase.BodemagWorkspace
ControlSystemsBase.DelayLtiSystem
ControlSystemsBase.DelayLtiSystem
ControlSystemsBase.LsimWorkspace
ControlSystemsBase.StateSpace
ControlSystemsBase.StaticStateSpace
ControlSystemsBase.StepInfo
ControlSystemsBase.TransferFunction
JuliaSimControl.AutoTuningProblem
JuliaSimControl.AutoTuningResult
JuliaSimControl.Eq29
JuliaSimControl.Eq32
JuliaSimControl.FixedGainObserver
JuliaSimControl.FunctionSystem
JuliaSimControl.FunctionSystem
JuliaSimControl.FunctionSystem
JuliaSimControl.GainMarginObjective
JuliaSimControl.LMI
JuliaSimControl.LinearSlidingModeController
JuliaSimControl.MPC.BoundsConstraint
JuliaSimControl.MPC.CollocationFinE
JuliaSimControl.MPC.CollocationFinE
JuliaSimControl.MPC.ControllerInput
JuliaSimControl.MPC.ControllerOutput
JuliaSimControl.MPC.DifferenceCost
JuliaSimControl.MPC.GenericMPCProblem
JuliaSimControl.MPC.LQMPCProblem
JuliaSimControl.MPC.LinearMPCConstraints
JuliaSimControl.MPC.LinearMPCModel
JuliaSimControl.MPC.MPCConstraints
JuliaSimControl.MPC.MPCIntegrator
JuliaSimControl.MPC.MPCParameters
JuliaSimControl.MPC.MultipleShooting
JuliaSimControl.MPC.MultipleShooting
JuliaSimControl.MPC.NonlinearMPCConstraints
JuliaSimControl.MPC.OSQPSolver
JuliaSimControl.MPC.Objective
JuliaSimControl.MPC.Objective
JuliaSimControl.MPC.ObjectiveInput
JuliaSimControl.MPC.ObserverInput
JuliaSimControl.MPC.QMPCProblem
JuliaSimControl.MPC.RobustMPCModel
JuliaSimControl.MPC.StageConstraint
JuliaSimControl.MPC.StageCost
JuliaSimControl.MPC.TerminalCost
JuliaSimControl.MPC.TerminalStateConstraint
JuliaSimControl.MPC.Trapezoidal
JuliaSimControl.MPC.Trapezoidal
JuliaSimControl.MaximumSensitivityObjective
JuliaSimControl.MaximumTransferObjective
JuliaSimControl.OperatingPoint
JuliaSimControl.OvershootObjective
JuliaSimControl.PhaseMarginObjective
JuliaSimControl.RiseTimeObjective
JuliaSimControl.SettlingTimeObjective
JuliaSimControl.SimulationObjective
JuliaSimControl.SlidingModeController
JuliaSimControl.StateFeedback
JuliaSimControl.Stateful
JuliaSimControl.StepTrackingObjective
JuliaSimControl.StructuredAutoTuningProblem
JuliaSimControl.StructuredAutoTuningResult
JuliaSimControl.SuperTwistingSMC
JuliaSimControl.SuperTwistingSMC
LowLevelParticleFilters.AdvancedParticleFilter
LowLevelParticleFilters.AuxiliaryParticleFilter
LowLevelParticleFilters.DAEUnscentedKalmanFilter
LowLevelParticleFilters.ExtendedKalmanFilter
LowLevelParticleFilters.KalmanFilter
LowLevelParticleFilters.ParticleFilter
LowLevelParticleFilters.SqKalmanFilter
LowLevelParticleFilters.TupleProduct
LowLevelParticleFilters.UnscentedKalmanFilter
LowLevelParticleFilters.UnscentedKalmanFilter
RobustAndOptimalControl.Disk
RobustAndOptimalControl.Diskmargin
RobustAndOptimalControl.ExtendedStateSpace
RobustAndOptimalControl.ExtendedStateSpace
RobustAndOptimalControl.LQGProblem
RobustAndOptimalControl.LQGProblem
RobustAndOptimalControl.NamedStateSpace
RobustAndOptimalControl.UncertainSS
RobustAndOptimalControl.nyquistcircles
RobustAndOptimalControl.δ
Base.step
CommonSolve.solve
CommonSolve.solve
CommonSolve.solve
CommonSolve.solve
CommonSolve.solve
ControlSystems.rlocus
ControlSystemsBase.G_CS
ControlSystemsBase.G_CS
ControlSystemsBase.G_PS
ControlSystemsBase.G_PS
ControlSystemsBase.add_input
ControlSystemsBase.add_output
ControlSystemsBase.append
ControlSystemsBase.are
ControlSystemsBase.are
ControlSystemsBase.array2mimo
ControlSystemsBase.balance
ControlSystemsBase.balance_statespace
ControlSystemsBase.balreal
ControlSystemsBase.baltrunc
ControlSystemsBase.bode
ControlSystemsBase.bodemag!
ControlSystemsBase.bodeplot
ControlSystemsBase.bodev
ControlSystemsBase.bodev
ControlSystemsBase.c2d
ControlSystemsBase.c2d
ControlSystemsBase.c2d
ControlSystemsBase.c2d_poly2poly
ControlSystemsBase.c2d_roots2poly
ControlSystemsBase.c2d_x0map
ControlSystemsBase.comp_sensitivity
ControlSystemsBase.controllability
ControlSystemsBase.covar
ControlSystemsBase.ctrb
ControlSystemsBase.d2c
ControlSystemsBase.d2c
ControlSystemsBase.d2c_exact
ControlSystemsBase.dab
ControlSystemsBase.damp
ControlSystemsBase.dampreport
ControlSystemsBase.dcgain
ControlSystemsBase.delay
ControlSystemsBase.delaymargin
ControlSystemsBase.diagonalize
ControlSystemsBase.evalfr
ControlSystemsBase.extended_gangoffour
ControlSystemsBase.feedback
ControlSystemsBase.feedback
ControlSystemsBase.feedback2dof
ControlSystemsBase.feedback2dof
ControlSystemsBase.freqresp
ControlSystemsBase.freqresp!
ControlSystemsBase.freqrespv
ControlSystemsBase.freqrespv
ControlSystemsBase.freqrespv
ControlSystemsBase.gangoffour
ControlSystemsBase.gangoffourplot
ControlSystemsBase.gangofseven
ControlSystemsBase.gram
ControlSystemsBase.grampd
ControlSystemsBase.hinfnorm
ControlSystemsBase.impulse
ControlSystemsBase.innovation_form
ControlSystemsBase.innovation_form
ControlSystemsBase.input_comp_sensitivity
ControlSystemsBase.input_comp_sensitivity
ControlSystemsBase.input_names
ControlSystemsBase.input_sensitivity
ControlSystemsBase.input_sensitivity
ControlSystemsBase.iscontinuous
ControlSystemsBase.isdiscrete
ControlSystemsBase.isproper
ControlSystemsBase.isstable
ControlSystemsBase.kalman
ControlSystemsBase.laglink
ControlSystemsBase.leadlink
ControlSystemsBase.leadlinkat
ControlSystemsBase.leadlinkcurve
ControlSystemsBase.lft
ControlSystemsBase.linfnorm
ControlSystemsBase.loopshapingPI
ControlSystemsBase.loopshapingPID
ControlSystemsBase.lqr
ControlSystemsBase.lsim
ControlSystemsBase.lsim
ControlSystemsBase.lsim
ControlSystemsBase.lsim!
ControlSystemsBase.margin
ControlSystemsBase.marginplot
ControlSystemsBase.markovparam
ControlSystemsBase.minreal
ControlSystemsBase.minreal
ControlSystemsBase.nicholsplot
ControlSystemsBase.nonlinearity
ControlSystemsBase.nyquist
ControlSystemsBase.nyquistplot
ControlSystemsBase.nyquistv
ControlSystemsBase.nyquistv
ControlSystemsBase.observability
ControlSystemsBase.observer_controller
ControlSystemsBase.observer_filter
ControlSystemsBase.observer_predictor
ControlSystemsBase.obsv
ControlSystemsBase.output_comp_sensitivity
ControlSystemsBase.output_comp_sensitivity
ControlSystemsBase.output_names
ControlSystemsBase.output_sensitivity
ControlSystemsBase.output_sensitivity
ControlSystemsBase.pade
ControlSystemsBase.pade
ControlSystemsBase.parallel
ControlSystemsBase.pid
ControlSystemsBase.pidplots
ControlSystemsBase.place
ControlSystemsBase.placePI
ControlSystemsBase.place_knvd
ControlSystemsBase.plyap
ControlSystemsBase.poles
ControlSystemsBase.pzmap
ControlSystemsBase.reduce_sys
ControlSystemsBase.relative_gain_array
ControlSystemsBase.relative_gain_array
ControlSystemsBase.rgaplot
ControlSystemsBase.rstc
ControlSystemsBase.rstd
ControlSystemsBase.sensitivity
ControlSystemsBase.series
ControlSystemsBase.setPlotScale
ControlSystemsBase.sigma
ControlSystemsBase.sigmaplot
ControlSystemsBase.sigmav
ControlSystemsBase.sigmav
ControlSystemsBase.similarity_transform
ControlSystemsBase.sminreal
ControlSystemsBase.ss
ControlSystemsBase.ss
ControlSystemsBase.ssdata
ControlSystemsBase.ssrand
ControlSystemsBase.stabregionPID
ControlSystemsBase.starprod
ControlSystemsBase.state_names
ControlSystemsBase.stepinfo
ControlSystemsBase.system_name
ControlSystemsBase.tf
ControlSystemsBase.thiran
ControlSystemsBase.time_scale
ControlSystemsBase.to_sized
ControlSystemsBase.tzeros
ControlSystemsBase.zpconv
ControlSystemsBase.zpk
ControlSystemsBase.zpkdata
DescriptorSystems.dss
JuliaSimControl.ESC
JuliaSimControl.MPC.IpoptSolver
JuliaSimControl.MPC.get_xu
JuliaSimControl.MPC.modelfit
JuliaSimControl.MPC.optimize!
JuliaSimControl.MPC.optimize!
JuliaSimControl.MPC.optimize!
JuliaSimControl.MPC.rms
JuliaSimControl.MPC.rollout
JuliaSimControl.MPC.step!
JuliaSimControl.MPC.step!
JuliaSimControl.MPC.step!
JuliaSimControl.OptimizedPID
JuliaSimControl.PIESC
JuliaSimControl.app_autotuning
JuliaSimControl.app_modelreduction
JuliaSimControl.build_K_function
JuliaSimControl.build_controlled_dynamics
JuliaSimControl.frequency_response_analysis
JuliaSimControl.frequency_response_analysis
JuliaSimControl.hinfsyn_lmi
JuliaSimControl.inverse_lqr
JuliaSimControl.inverse_lqr
JuliaSimControl.inverse_lqr
JuliaSimControl.inverse_lqr
JuliaSimControl.ispassive_lmi
JuliaSimControl.mussv
JuliaSimControl.mussv
JuliaSimControl.mussv
JuliaSimControl.mussv_DG
JuliaSimControl.mussv_tv
JuliaSimControl.mussv_tv
JuliaSimControl.poly_approx
JuliaSimControl.poly_approx
JuliaSimControl.pqr
JuliaSimControl.predicted_cost
JuliaSimControl.rk4
JuliaSimControl.simplified_system
JuliaSimControl.spr_synthesize
LinearAlgebra.lyap
LinearAlgebra.norm
LowLevelParticleFilters.commandplot
LowLevelParticleFilters.correct!
LowLevelParticleFilters.correct!
LowLevelParticleFilters.correct!
LowLevelParticleFilters.debugplot
LowLevelParticleFilters.densityplot
LowLevelParticleFilters.forward_trajectory
LowLevelParticleFilters.forward_trajectory
LowLevelParticleFilters.log_likelihood_fun
LowLevelParticleFilters.loglik
LowLevelParticleFilters.logsumexp!
LowLevelParticleFilters.mean_trajectory
LowLevelParticleFilters.mean_trajectory
LowLevelParticleFilters.metropolis
LowLevelParticleFilters.reset!
LowLevelParticleFilters.reset!
LowLevelParticleFilters.reset!
LowLevelParticleFilters.simulate
LowLevelParticleFilters.smooth
LowLevelParticleFilters.smooth
LowLevelParticleFilters.smoothed_cov
LowLevelParticleFilters.smoothed_mean
LowLevelParticleFilters.smoothed_trajs
LowLevelParticleFilters.update!
LowLevelParticleFilters.weighted_cov
LowLevelParticleFilters.weighted_mean
LowLevelParticleFilters.weighted_mean
RobustAndOptimalControl.add_disturbance
RobustAndOptimalControl.add_input_differentiator
RobustAndOptimalControl.add_input_integrator
RobustAndOptimalControl.add_low_frequency_disturbance
RobustAndOptimalControl.add_low_frequency_disturbance
RobustAndOptimalControl.add_measurement_disturbance
RobustAndOptimalControl.add_output_differentiator
RobustAndOptimalControl.add_output_integrator
RobustAndOptimalControl.add_resonant_disturbance
RobustAndOptimalControl.add_resonant_disturbance
RobustAndOptimalControl.baltrunc2
RobustAndOptimalControl.baltrunc_coprime
RobustAndOptimalControl.baltrunc_unstab
RobustAndOptimalControl.bilinearc2d
RobustAndOptimalControl.bilinearc2d
RobustAndOptimalControl.bilinearc2d
RobustAndOptimalControl.bilineard2c
RobustAndOptimalControl.bilineard2c
RobustAndOptimalControl.bilineard2c
RobustAndOptimalControl.blocksort
RobustAndOptimalControl.broken_feedback
RobustAndOptimalControl.closedloop
RobustAndOptimalControl.connect
RobustAndOptimalControl.controller_reduction
RobustAndOptimalControl.controller_reduction_plot
RobustAndOptimalControl.controller_reduction_weight
RobustAndOptimalControl.dare3
RobustAndOptimalControl.diskmargin
RobustAndOptimalControl.diskmargin
RobustAndOptimalControl.diskmargin
RobustAndOptimalControl.expand_symbol
RobustAndOptimalControl.extended_controller
RobustAndOptimalControl.extended_controller
RobustAndOptimalControl.feedback_control
RobustAndOptimalControl.ff_controller
RobustAndOptimalControl.find_lft
RobustAndOptimalControl.fit_complex_perturbations
RobustAndOptimalControl.frequency_separation
RobustAndOptimalControl.frequency_weighted_reduction
RobustAndOptimalControl.fudge_inv
RobustAndOptimalControl.gain_and_delay_uncertainty
RobustAndOptimalControl.gainphaseplot
RobustAndOptimalControl.glover_mcfarlane
RobustAndOptimalControl.glover_mcfarlane
RobustAndOptimalControl.glover_mcfarlane_2dof
RobustAndOptimalControl.h2norm
RobustAndOptimalControl.h2synthesize
RobustAndOptimalControl.hankelnorm
RobustAndOptimalControl.hanus
RobustAndOptimalControl.hess_form
RobustAndOptimalControl.hinfassumptions
RobustAndOptimalControl.hinfgrad
RobustAndOptimalControl.hinfnorm2
RobustAndOptimalControl.hinfpartition
RobustAndOptimalControl.hinfsignals
RobustAndOptimalControl.hinfsynthesize
RobustAndOptimalControl.hsvd
RobustAndOptimalControl.ispassive
RobustAndOptimalControl.loop_diskmargin
RobustAndOptimalControl.loop_diskmargin
RobustAndOptimalControl.loop_scale
RobustAndOptimalControl.loop_scaling
RobustAndOptimalControl.lqr3
RobustAndOptimalControl.makeweight
RobustAndOptimalControl.measure
RobustAndOptimalControl.modal_form
RobustAndOptimalControl.muplot
RobustAndOptimalControl.mvnyquistplot
RobustAndOptimalControl.named_ss
RobustAndOptimalControl.named_ss
RobustAndOptimalControl.named_ss
RobustAndOptimalControl.ncfmargin
RobustAndOptimalControl.neglected_delay
RobustAndOptimalControl.neglected_lag
RobustAndOptimalControl.noise_mapping
RobustAndOptimalControl.nu_reduction
RobustAndOptimalControl.nu_reduction_recursive
RobustAndOptimalControl.nugap
RobustAndOptimalControl.nugap
RobustAndOptimalControl.partition
RobustAndOptimalControl.partition
RobustAndOptimalControl.passivity_index
RobustAndOptimalControl.passivityplot
RobustAndOptimalControl.performance_mapping
RobustAndOptimalControl.robstab
RobustAndOptimalControl.schur_form
RobustAndOptimalControl.show_construction
RobustAndOptimalControl.sim_diskmargin
RobustAndOptimalControl.sim_diskmargin
RobustAndOptimalControl.sim_diskmargin
RobustAndOptimalControl.specificationplot
RobustAndOptimalControl.splitter
RobustAndOptimalControl.ss2particles
RobustAndOptimalControl.ssdata_e
RobustAndOptimalControl.stab_unstab
RobustAndOptimalControl.static_gain_compensation
RobustAndOptimalControl.structured_singular_value
RobustAndOptimalControl.structured_singular_value
RobustAndOptimalControl.sumblock
RobustAndOptimalControl.sys_from_particles
RobustAndOptimalControl.system_mapping
RobustAndOptimalControl.uss
RobustAndOptimalControl.uss
RobustAndOptimalControl.uss
RobustAndOptimalControl.uss
RobustAndOptimalControl.vec2sys
RobustAndOptimalControl.δc
RobustAndOptimalControl.δr
StatsAPI.predict!
StatsAPI.predict!
StatsAPI.predict!
Docstrings
JuliaSimControl
Docstrings of the MPC submodule are located under MPC.
JuliaSimControl.Eq29
— TypeEq29 <: InverseLQRMethod
R
is a weighting matrix, determining the relative importance of matching each input to that provided by the original controller.
JuliaSimControl.Eq32
— TypeEq32 <: InverseLQRMethod
Eq32(q1, r1)
r1
penalizes deviation from the original control signal while q1
penalizes deviation from the original effect the cotrol signal had on the state. In other words, the r1
term tries to use the same actuator configuration as the original controller, while the q1
term ensures that the effec of the controller is the same.
JuliaSimControl.LMI
— TypeLMI <: InverseLQRMethod
R
is a weighting matrix, determining the relative importance of matching each input to that provided by the original controller.
JuliaSimControl.inverse_lqr
— MethodQ,R,S,P,L2 = inverse_lqr(method::Eq29, G::AbstractStateSpace, L)
Solve the inverse optimal control problem that finds the LQR cost function that leads to a controller approximating state feedback controller with gain matrix L
. P
is the solution to the Riccati equation and L2
is the recreated feedback gain. Note: S
will in general not be zero, and including this cross term in the cost function may be important. If including S
is not possible, use the LMI
method to find a cost function with as small S
as possible.
Creates the stage-cost matrix
\[\begin{bmatrix} L^T R L & -L^T R\\ -R L & R \end{bmatrix} = \begin{bmatrix} Q & -S\\ -S^T & R \end{bmatrix} \]
Ref: "Designing MPC controllers by reverse-engineering existing LTI controllers", E. N. Hartley, J. M. Maciejowski
JuliaSimControl.inverse_lqr
— MethodQ,R,S,P,L2 = inverse_lqr(method::Eq32, G::AbstractStateSpace, L)
Solve the inverse optimal control problem that finds the LQR cost function that leads to a controller approximating state feedback controller with gain matrix L
. P
is the solution to the Riccati equation and L2
is the recreated feedback gain. Note: S
will in general not be zero, and including this cross term in the cost function may be important. If including S
is not possible, use the LMI
method to find a cost function with as small S
as possible.
Creates the cost function
\[||Bu - BLu||^2_{q_1} + ||u - Lu||^2_{r_1}\]
Ref: "Designing MPC controllers by reverse-engineering existing LTI controllers", E. N. Hartley, J. M. Maciejowski
JuliaSimControl.inverse_lqr
— MethodQ,R,S,P,L2 = inverse_lqr(method::GMF)
Solve the inverse optimal control problem that finds the LQR cost function on the form
\[x'Qx + u'Ru\]
that leads to an LQR controller approximating a Glover McFarlane controller. Using this method, S = 0
.
Ref: Rowe and Maciejowski, "Tuning MPC using H∞ Loop Shaping" .
Example
disc = (x) -> c2d(ss(x), 0.01)
G = tf([1, 5], [1, 2, 10]) |> disc # Plant model
W1 = tf(1,[1, 0]) |> disc # Loop shaping weight
gmf2 = glover_mcfarlane(G; W1) # Design Glover McFarlane controller
Q,R,S,P,L = inverse_lqr(GMF(gmf2)) # Get cost function
using Test
L3 = lqr(G*W1, Q, R) # Test equivalence to LQR
@test L3 ≈ -L2
JuliaSimControl.inverse_lqr
— MethodQ,R,S,P,L2 = inverse_lqr(method::LMI, G::AbstractStateSpace, L; optimizer = Hypatia.Optimizer)
Solve the inverse optimal control problem that finds the LQR cost function that leads to a controller approximating state feedback controller with gain matrix L
. P
is the solution to the Riccati equation and L2
is the recreated feedback gain. Note, L
is supposed to be used with negative feedback, i.e., it's designed such that u = -Lx
.
Solves a convex LMI problem minimizing the cross term S
. Note: there is no guarantee that the corss term will be driven to 0 exactly. If S
remains large in relation to Q
and R
, S
must be included in the cost function for high fidelity reproduction of the original controller.
Ref: "Designing MPC controllers by reverse-engineering existing LTI controllers", E. N. Hartley, J. M. Maciejowski
JuliaSimControl.FixedGainObserver
— TypeFixedGainObserver{F <: Function, x} <: AbstractFilter
FixedGainObserver(sys::AbstractStateSpace, x0, K)
A linear observer, similar to a Kalman filer, but with a fixed measurement feedback gain. The gain can be designed using, e.g., pole placement or solving a Riccati equation. For a robust observer, consider using glover_mcfarlane
followed by inverse_lqr
.
JuliaSimControl.OperatingPoint
— TypeOperatingPoint(x, u, y)
OperatingPoint()
Structure representing an operating point around which a system is linearized. If no arguments are supplied, an empty operating point is created.
Arguments:
x
: Stateu
: Control inputy
: Output
JuliaSimControl.StateFeedback
— TypeStateFeedback{F <: Function, x} <: AbstractFilter
StateFeedback(discrete_dynamics, x0, nu, ny)
StateFeedback(sys::FunctionSystem, x0)
An observer that uses the dynamics model without any measurement feedback. This observer can be used as an oracle that has full knowledge of the state. Note, this is often an unrealistic assumption in real-world contexts and open-loop observers can not account for load disturbances. Use of this observer in a closed-loop context creates a false closed loop.
LowLevelParticleFilters.UnscentedKalmanFilter
— MethodUnscentedKalmanFilter(sys::FunctionSystem, R1, R2, d0=MvNormal(Matrix(R1)); p = SciMLBase.NullParameters())
Convencience constructor for systems of type FunctionSystem
.
JuliaSimControl.mussv
— Methodmussv(M::AbstractMatrix; optimizer = Hypatia.Optimizer{eltype(M)}, bu = opnorm(M), tol = 0.001)
Compute (an upper bound of) the structured singular value μ for diagonal Δ of complex perturbations (other structures of Δ are handled by supplying the block structure mussv(M, blocks)). M
is assumed to be an (n × n × N_freq) array or a matrix.
Solves a convex LMI.
Arguments:
bu
: Upper bound for bisection.tol
: tolerance.
See also mussv
, mussv_tv
, mussv_DG
.
Extended help
By default, the Hypatia solver (native Julia) is used to solve the problem. Other solvers that handle the this problem formulation is
- Mosek (requires license)
- SCS
The solver can be selected using the keyword argument optimizer
.
JuliaSimControl.mussv
— Methodmussv(M::AbstractStateSpace, blocks; optimizer = Hypatia.Optimizer, bu0 = 20, tol = 0.001)
Compute (an upper bound of) the structured singular value μ of statespace model M
interconnected with uncertainty structure described by blocks
. Reference: MAE509 (LMIs in Control): Lecture 14, part C - LMIs for Robust Control with Structured Uncertainty
Example:
The following example illustrates how to use the structured singular value to determine how large diagonal complex uncertainty can be added at the input of a plant P
before the closed-loop system becomes unstable
julia> Δ = uss([δc(), δc()]); # Diagonal uncertainty element
julia> a = 1;
julia> P = ss([0 a; -a -1], I(2), [1 a; 0 1], 0) * (I(2) + Δ);
julia> K = ss(I(2));
julia> G = lft(P, -K);
julia> stabmargin = 1/mussv(G) # We can handle 134% of the modeled uncertainty
1.3429508196721311
julia> # Compare with the input diskmargin
diskmargin(K*system_mapping(P), -1)
Disk margin with:
Margin: 1.3469378397689664
Frequency: -0.40280561122244496
Gain margins: [-0.3469378397689664, 2.3469378397689664]
Phase margin: 84.67073122411068
Skew: -1
Worst-case perturbation: missing
Extended help
By default, the Hypatia solver is used to solve the problem. Other solvers that handle the this problem formulation is
- Mosek (requires license)
- Hypatia.jl (native Julia)
- Clarabel.jl (native Julia)
- SCS (typically performs poorly for this problem)
The solver can be selected using the keyword argument optimizer
.
JuliaSimControl.mussv
— Methodmussv(G::UncertainSS; kwargs...)
Compute (an upper bound of) the structured singular value μ of uncertain system model G
.
JuliaSimControl.mussv_DG
— Methodmussv_DG(M::AbstractMatrix; optimizer = Hypatia.Optimizer{eltype(M)}, bu = opnorm(M), tol = 0.001)
Compute (an upper bound of) the structured singular value μ for diagonal Δ of real perturbations (other structures of Δ are not yet supported). M
is assumed to be an (n × n × N_freq) array or a matrix.
See mussv
for more details. See also mussv
, mussv_tv
, mussv_DG
.
JuliaSimControl.mussv_tv
— Methodmussv_tv(G::AbstractStateSpace, blocks; optimizer = Hypatia.Optimizer, bu = 20, tol = 0.001)
Compute (an upper bound of) the structured singular value margin when G
is interconnected with time-varying uncertainty structure described by blocks
. This value will in general be larger than the one returned by mussv
, but can be used to guarantee stability for infinitely fast time-varying perturbations, i.e., if the return value is < 1, the system is stable no matter how fast the dynamics of the perturbations change.
The result will in general be more accurate if G
is passed rather than a matrix M
, unless a very dense grid around the critical frequency is used for to calculate M
.
Solves a convex LMI.
Arguments:
bu
: Upper bound for bisection.tol
: tolerance.
JuliaSimControl.mussv_tv
— Methodmussv_tv(M::AbstractArray{<:Any, 3}; optimizer = Hypatia.Optimizer, bu = 20, tol = 0.001)
mussv_tv(G::UncertainSS; optimizer = Hypatia.Optimizer, bu0 = 20, tol = 0.001)
Compute (an upper bound of) the structured singular value μ for diagonal, complex and time-varying Δ(t) using constant (over frequency) matrix scalings. This value will in general be larger than the one returned by mussv
, but can be used to guarantee stability for infinitely fast time-varying perturbations, i.e., if the return value is < 1, the system is stable no matter how fast the dynamics of the perturbations change.
M
is assumed to be an (n × n × N_freq) array or a matrix. G
is an UncertainSS
. The result will in general be more accurate if G
is passed rather than M
, unless a very dense grid around the critical frequency is used for to calculate M
.
Solves a convex LMI.
Arguments:
bu
: Upper bound for bisection.tol
: tolerance.
See also mussv
, mussv_tv
, mussv_DG
.
Extended help
The normal μ is calculated by minimizing
\[\operatorname{min}_D ||D(\omega) M(\omega) D(\omega)^{-1}||\]
where a unique $D(\omega)$ is allowed for each frequency. However, in this problem, the $D$ scalings are constant over frequency, yielding a more conservative answer, with the additional bonus of being applicable for time-varying perturbations.
This strong guarantee can be used to prove stability of nonlinear systems by formulating them as linear systems with time-varying, norm-bounded perturbations. For such systems, mussv_tv < 1
is a sufficient condition for stability. See Boyd et al., "Linear Matrix Inequalities in System and Control Theory" for more details.
JuliaSimControl.hinfsyn_lmi
— MethodK, γ = hinfsyn_lmi(P::ExtendedStateSpace;
opt = Hypatia.Optimizer(), γrel = 1.01, ϵ = 1e-3, balance = true, perm = false,)
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. This implementation solves a convex LMI problem.
Arguments:
opt
: A MathOptInterface solver constructor.γ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
.ϵ
: A small offset to enforce strict LMI inequalities. This can be tuned if numerical issues arise.balance
: Perform a balancing transformation onP
usingControlSystemsBase.balance_statespace
.perm
: Ifbalance=true, perm=true
, the balancing transform is allowed to permute the state vector. This is not allowed by default, but can improve the numerics slightly if allowed.
The Hypatia solver takes the following arguments https://github.com/chriscoey/Hypatia.jl/blob/42e4b10318570ea22adb39fec1c27d8684161cec/src/Solvers/Solvers.jl#L73
JuliaSimControl.ispassive_lmi
— Methodispassive_lmi(P::AbstractStateSpace{ControlSystemsBase.Continuous}; ftype = Float64, opt = Hypatia.Optimizer{ftype}(), ϵ = 0.001, balance = true, verbose = true, silent_solver = true)
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.
This functions solves a convex LMI related to the KYP (positive real) lemma.
Arguments:
balance
: Balance the system before calculations?verbose
: Print status messages
JuliaSimControl.spr_synthesize
— MethodK, Gcl = spr_synthesize(P0::ExtendedStateSpace{Continuous};, opt = Hypatia.Optimizer, balance = true, verbose = true, silent_solver = true, ϵ = 1e-6)
Design a strictly positive real controller (passive) that optimizes the closed-loop H₂-norm subject to being passive.
For plants that are known to be passive, control using a passive controller is guaranteed to be stable.
The returned controller is supposed to be used with positive feedback, so ispassive(-K)
should hold. The resulting closed-loop system from disturbances to performance outputs is also returned, Gcl = lft(P0, K)
.
Implements the algorithm labeled as "Pseudocode 1" in "Synthesis of strictly positive real H2 controllers using dilated LMI", Forbes 2018
Arguments:
P0
: AnExtendedStateSpace
object. This object can be designed using H₂ or H∞ methods. See, e.g.,hinfpartition
.opt
: A JuMP compatible solver.balance
: Perform balancing of the statespace system before solving.verbose
: Print info?silent_solver
:ϵ
: A small numerical constant to enforce strict positive definiteness.
See also h2synthesize
, hinfsynthesize
.
JuliaSimControl.FunctionSystem
— MethodFunctionSystem(sys::ODESystem, u, y; z=nothing, w=nothing, kwargs...)
Generate code for the dynamics and output of sys
by calling build_controlled_dynamics
. See build_controlled_dynamics
for more details.
JuliaSimControl.build_controlled_dynamics
— Functionf, xout, pout = build_controlled_dynamics(sys, u; kwargs...)
f, obsf, xout, pout = build_controlled_dynamics(sys, u, y; z=nothing, w=nothing, kwargs...)
Build a function on the form (x,u,p,t) -> ẋ where
x
are the differential statesu
are control inputsp
are parameterst
is timekwargs
are sent toModelingToolkit.build_function
f
is a tuple of functions, one out of palce and one in place(x,u,p,t) -> ẋ
and(ẋ,x,u,p,t) -> nothing
xout
contains the order of the states included in the dynamicspout
contains the order of the parameters included in the dynamics
If in addition to u
, outputs y
are also specified, an additional observed function tuple is returned.
Example
This example forms a feedback system and builds a function of the dynamics from the reference r
to the output y
.
using ModelingToolkit, Test
@variables t x(t)=0 y(t)=0 u(t)=0 r(t)=0
@parameters kp=1
D = Differential(t)
eqs = [
u ~ kp * (r - y) # P controller
D(x) ~ -x + u # Plant dynamics
y ~ x # Output equation
]
@named sys = ODESystem(eqs, t)
funsys = JuliaSimControl.build_controlled_dynamics(sys, r, y; checkbounds=true)
x = zeros(funsys.nx) # sys.x
u = [1] # r
p = [1] # kp
xd = funsys(x,u,p,1)
varmap = Dict(
funsys.x[] => 1,
kp => 1,
)
@test xd == ModelingToolkit.varmap_to_vars(varmap, funsys.x)
JuliaSimControl.FunctionSystem
— TypeFunctionSystem{TE <: ControlSystemsBase.TimeEvolution, F, G}
A structure representing the dynamical system
\[�egin{aligned} x′ &= f(x,u,p,t)\ y &= g(x,u,p,t) �nd{aligned}\]
To build a FunctionSystem for a DAE on mass-matrix form, use an ODEFunction
as f
f = ODEFunction(dae_dyn, mass_matrix = M)
To obtain the simplified system and default values for the initial condition and parameters, see simplified_system
.
Fields:
dynamics::F
measurement::G
timeevol::TE
: ControlSystemsBase.TimeEvolutionx
: statesu
: controlled inputsy
: measured outputsw
: disturbance inputsz
: performance outputsp
: parameters
JuliaSimControl.FunctionSystem
— MethodFunctionSystem(f, g; kwargs...)
FunctionSystem(f, g, Ts::Real; kwargs...)
Constructor for FunctionSystem
.
Arguments:
f
: Discrete dynamics with signature (x,u,p,t)g
: Measurement function with signature (x,u,p,t)Ts
: If the sample timeTs
is provided, the system represents a discrete-time system, otherwise the dynamics is assumed to be continuous.kwargs
: Signal names
JuliaSimControl.simplified_system
— Methodsimplified_system(funcsys::FunctionSystem)
Obtain the result from structural_simplify
that was obtained after input-output processing. This system is often of a lower order than the original system. To obtain the default initial condition and parameters of the simplified system, call
ssys = simplified_system(funcsys)
defs = ModelingToolkit.defaults(ssys)
x0, p0 = ModelingToolkit.get_u0_p(ssys, defs, defs)
ControlSystems and RobustAndOptimalControl
ControlSystems.Simulator
— TypeSimulator
Fields:
P::StateSpace
f = (x,p,t) -> x
y = (x,t) -> y
ControlSystems.Simulator
— MethodSimulator(P::StateSpace, u = (x,t) -> 0)
Used to simulate continuous-time systems. See function ?solve
for additional info.
Usage:
using OrdinaryDiffEq, Plots
dt = 0.1
tfinal = 20
t = 0:dt:tfinal
P = ss(tf(1,[2,1])^2)
K = 5
reference(x,t) = [1.]
s = Simulator(P, reference)
x0 = [0.,0]
tspan = (0.0,tfinal)
sol = solve(s, x0, tspan, Tsit5())
plot(t, s.y(sol, t)[:], lab="Open loop step response")
ControlSystems.rlocus
— Methodroots, Z, K = rlocus(P::LTISystem, K = 500)
Compute the root locus of the SISO LTISystem P
with a negative feedback loop and feedback gains between 0 and K
. rlocus
will use an adaptive step-size algorithm to determine the values of the feedback gains used to generate the plot.
roots
is a complex matrix containing the poles trajectories of the closed-loop 1+k⋅G(s)
as a function of k
, Z
contains the zeros of the open-loop system G(s)
and K
the values of the feedback gain.
ControlSystemsBase.lsim
— Methodres = lsim(sys::DelayLtiSystem, u, t::AbstractArray{<:Real}; x0=fill(0.0, nstates(sys)), alg=MethodOfSteps(Tsit5()), abstol=1e-6, reltol=1e-6, force_dtmin=true, kwargs...)
Simulate system sys
, over time t
, using input signal u
, with initial state x0
, using method alg
.
Arguments:
t
: Has to be an AbstractVector
with equidistant time samples (t[i] - t[i-1]
constant) u
: Function to determine control signal ut
at a time t
, on any of the following forms:
u
: Function to determine control signaluₜ
at a timet
, on any of the following forms:- A constant
Number
orVector
, interpreted as a constant input. - Function
u(x, t)
that takes the internal state and time, note, the state representation for delay systems is not the same as for rational systems. - In-place function
u(uₜ, x, t)
. (Slightly more efficient)
- A constant
alg, abstol, reltol
and kwargs...
: are sent to DelayDiffEq.solve
.
This methods sets force_dtmin=true
by default to handle the discontinuity implied by, e.g., step inputs. This may lead to the solver taking a long time to solve ill-conditioned problems rather than exiting early with a warning.
Returns an instance of SimResult
which can be plotted directly or destructured into y, t, x, u = res
.
ControlSystemsBase.lsim
— Methodlsim(sys::HammersteinWienerSystem, u, t::AbstractArray{<:Real}; x0=fill(0.0, nstates(sys)), alg=Tsit5(), abstol=1e-6, reltol=1e-6, kwargs...)
Simulate system sys
, over time t
, using input signal u
, with initial state x0
, using method alg
.
Arguments:
t
: Has to be anAbstractVector
with equidistant time samples (t[i] - t[i-1]
constant)u
: Function to determine control signaluₜ
at a timet
, on any of the following forms: Can be a constantNumber
orVector
, interpreted asuₜ .= u
, or Functionuₜ .= u(x, t)
, or In-place functionu(uₜ, x, t)
. (Slightly more efficient)alg, abstol, reltol
andkwargs...
: are sent toOrdinaryDiffEq.solve
.
Returns an instance of SimResult
.
ControlSystemsBase.BodemagWorkspace
— MethodBodemagWorkspace(sys::LTISystem, N::Int)
BodemagWorkspace(sys::LTISystem, ω::AbstractVector)
BodemagWorkspace(R::Array{Complex{T}, 3}, mag::Array{T, 3})
BodemagWorkspace{T}(ny, nu, N)
Generate a workspace object for use with the in-place function bodemag!
. N
is the number of frequency points, alternatively, the input ω
can be provided instead of N
. Note: for threaded applications, create one workspace object per thread.
Arguments:
mag
: The output array ∈ 𝐑(ny, nu, nω)R
: Frequency-response array ∈ 𝐂(ny, nu, nω)
ControlSystemsBase.DelayLtiSystem
— Typestruct DelayLtiSystem{T, S <: Real} <: LTISystem
Represents an LTISystem with internal time-delay. See ?delay
for a convenience constructor.
ControlSystemsBase.DelayLtiSystem
— MethodDelayLtiSystem{T, S}(sys::StateSpace, Tau::AbstractVector{S}=Float64[]) where {T <: Number, S <: Real}
Create a delayed system by specifying both the system and time-delay vector. NOTE: if you want to create a system with simple input or output delays, use the Function delay(τ)
.
ControlSystemsBase.LsimWorkspace
— MethodLsimWorkspace(sys::AbstractStateSpace, N::Int)
LsimWorkspace(sys::AbstractStateSpace, u::AbstractMatrix)
LsimWorkspace{T}(ny, nu, nx, N)
Generate a workspace object for use with the in-place function lsim!
. sys
is the discrete-time system to be simulated and N
is the number of time steps, alternatively, the input u
can be provided instead of N
. Note: for threaded applications, create one workspace object per thread.
ControlSystemsBase.StateSpace
— TypeStateSpace{TE, T} <: AbstractStateSpace{TE}
An object representing a standard state space system.
See the function ss
for a user facing constructor as well as the documentation page creating systems.
Fields:
A::Matrix{T}
B::Matrix{T}
C::Matrix{T}
D::Matrix{T}
timeevol::TE
ControlSystemsBase.StaticStateSpace
— MethodStaticStateSpace(sys::AbstractStateSpace)
Return a HeteroStateSpace
where the system matrices are of type SMatrix.
NOTE: This function is fundamentally type unstable. For maximum performance, create the static system manually, or make use of the function-barrier technique.
ControlSystemsBase.StepInfo
— TypeStepInfo
Computed using stepinfo
Fields:
y0
: The initial value of the step response.yf
: The final value of the step response.stepsize
: The size of the step.peak
: The peak value of the step response.peaktime
: The time at which the peak occurs.overshoot
: The overshoot of the step response.settlingtime
: The time at which the step response has settled to withinsettling_th
of the final value.settlingtimeind::Int
: The index at which the step response has settled to withinsettling_th
of the final value.risetime
: The time at which the response rises fromrisetime_th[1]
torisetime_th[2]
of the final valuei10::Int
: The index at which the response reachesrisetime_th[1]
i90::Int
: The index at which the response reachesrisetime_th[2]
res::SimResult{SR}
: The simulation result used to compute the step response characteristics.settling_th
: The threshold used to computesettlingtime
andsettlingtimeind
.risetime_th
: The thresholds used to computerisetime
,i10
, andi90
.
ControlSystemsBase.TransferFunction
— MethodF(s)
, F(omega, true)
, F(z, false)
Notation for frequency response evaluation.
- F(s) evaluates the continuous-time transfer function F at s.
- F(omega,true) evaluates the discrete-time transfer function F at exp(imTsomega)
- F(z,false) evaluates the discrete-time transfer function F at z
Base.step
— Methody, t, x = step(sys[, tfinal])
y, t, x = step(sys[, t])
Calculate the response of the system sys
to a unit step at time t = 0
. If the final time tfinal
or time vector t
is not provided, one is calculated based on the system pole locations.
The return value is a structure of type SimResult
. A SimResul
can be plotted by plot(result)
, or destructured as y, t, x = result
.
y
has size (ny, length(t), nu)
, x
has size (nx, length(t), nu)
ControlSystemsBase.G_CS
— MethodG_CS(P, C)
The closed-loop transfer function from (-) measurement noise or (+) reference to control signal. Technically, the transfer function is given by (1 + CP)⁻¹C
so SC
would be a better, but nonstandard name.
▲
│e₁
│ ┌─────┐
d₁────+──┴──► P ├─────┬──►e₄
│ └─────┘ │
│ │
│ ┌─────┐ -│
e₂◄──┴─────┤ C ◄──┬──+───d₂
└─────┘ │
│e₃
▼
input_sensitivity
is the transfer function from d₁ to e₁, (I + CP)⁻¹output_sensitivity
is the transfer function from d₂ to e₃, (I + PC)⁻¹input_comp_sensitivity
is the transfer function from d₁ to e₂, (I + CP)⁻¹CPoutput_comp_sensitivity
is the transfer function from d₂ to e₄, (I + PC)⁻¹PCG_PS
is the transfer function from d₁ to e₄, (1 + PC)⁻¹PG_CS
is the transfer function from d₂ to e₂, (1 + CP)⁻¹C
ControlSystemsBase.G_PS
— MethodG_PS(P, C)
The closed-loop transfer function from load disturbance to plant output. Technically, the transfer function is given by (1 + PC)⁻¹P
so SP
would be a better, but nonstandard name.
▲
│e₁
│ ┌─────┐
d₁────+──┴──► P ├─────┬──►e₄
│ └─────┘ │
│ │
│ ┌─────┐ -│
e₂◄──┴─────┤ C ◄──┬──+───d₂
└─────┘ │
│e₃
▼
input_sensitivity
is the transfer function from d₁ to e₁, (I + CP)⁻¹output_sensitivity
is the transfer function from d₂ to e₃, (I + PC)⁻¹input_comp_sensitivity
is the transfer function from d₁ to e₂, (I + CP)⁻¹CPoutput_comp_sensitivity
is the transfer function from d₂ to e₄, (I + PC)⁻¹PCG_PS
is the transfer function from d₁ to e₄, (1 + PC)⁻¹PG_CS
is the transfer function from d₂ to e₂, (1 + CP)⁻¹C
ControlSystemsBase.add_input
— Functionadd_input(sys::AbstractStateSpace, B2::AbstractArray, D2 = 0)
Add inputs to sys
by forming
\[\begin{aligned} x' &= Ax + [B \; B_2]u \\ y &= Cx + [D \; D_2]u \\ \end{aligned}\]
If B2
is an integer it will be interpreted as an index and an input matrix containing a single 1 at the specified index will be used.
Example: The following example forms an innovation model that takes innovations as inputs
G = ssrand(2,2,3, Ts=1)
K = kalman(G, I(G.nx), I(G.ny))
sys = add_input(G, K)
ControlSystemsBase.add_output
— Functionadd_output(sys::AbstractStateSpace, C2::AbstractArray, D2 = 0)
Add outputs to sys
by forming
\[\begin{aligned} x' &= Ax + Bu \\ y &= [C; C_2]x + [D; D_2]u \\ \end{aligned}\]
If C2
is an integer it will be interpreted as an index and an output matrix containing a single 1 at the specified index will be used.
ControlSystemsBase.append
— Methodappend(systems::StateSpace...), append(systems::TransferFunction...)
Append systems in block diagonal form
ControlSystemsBase.are
— Methodare(::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
— Methodare(::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.array2mimo
— Methodarray2mimo(M::AbstractArray{<:LTISystem})
Take an array of LTISystem
s and create a single MIMO system.
ControlSystemsBase.balance
— FunctionS, 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.balance_statespace
— FunctionA, 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 [TA/T TB; 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
ControlSystemsBase.balreal
— Methodsysr, 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
— Methodsysr, 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.bode
— Methodmag, phase, w = bode(sys[, w]; unwrap=true)
Compute the magnitude and phase parts of the frequency response of system sys
at frequencies w
. See also bodeplot
mag
and phase
has size (ny, nu, length(w))
. If unwrap
is true (default), the function unwrap!
will be applied to the phase angles. This procedure is costly and can be avoided if the unwrapping is not required.
For higher performance, see the function bodemag!
that computes the magnitude only.
ControlSystemsBase.bodemag!
— Methodmag = bodemag!(ws::BodemagWorkspace, sys::LTISystem, w::AbstractVector)
Compute the Bode magnitude operating in-place on an instance of BodemagWorkspace
. Note that the returned magnitude array is aliased with ws.mag
. The output array mag
is ∈ 𝐑(ny, nu, nω).
ControlSystemsBase.bodeplot
— Functionfig = bodeplot(sys, args...)
bodeplot(LTISystem[sys1, sys2...], args...; plotphase=true, balance = true, kwargs...)
Create a Bode plot of the LTISystem
(s). A frequency vector w
can be optionally provided. To change the Magnitude scale see setPlotScale
. The default magnitude scale is "log10" (absolute scale).
- If
hz=true
, the plot x-axis will be displayed in Hertz, the input frequency vector is still treated as rad/s. balance
: Callbalance_statespace
on the system before plotting.
kwargs
is sent as argument to RecipesBase.plot.
ControlSystemsBase.bodev
— Methodbodev(sys::LTISystem, w::AbstractVector; $(Expr(:kw, :unwrap, true)))
For use with SISO systems where it acts the same as bode
but with the extra dimensions removed in the returned values.
ControlSystemsBase.bodev
— Methodbodev(sys::LTISystem; )
For use with SISO systems where it acts the same as bode
but with the extra dimensions removed in the returned values.
ControlSystemsBase.c2d
— Functionc2d(G::DelayLtiSystem, Ts, method=:zoh)
ControlSystemsBase.c2d
— Functionsysd = c2d(sys::AbstractStateSpace{<:Continuous}, Ts, method=:zoh; w_prewarp=0)
Gd = c2d(G::TransferFunction{<:Continuous}, Ts, method=:zoh)
Convert the continuous-time system sys
into a discrete-time system with sample time Ts
, using the specified method
(:zoh
, :foh
, :fwdeuler
or :tustin
).
method = :tustin
performs a bilinear transform with prewarp frequency w_prewarp
.
w_prewarp
: Frequency (rad/s) for pre-warping when using the Tustin method, has no effect for other methods.
See also c2d_x0map
Extended help
ZoH sampling is exact for linear systems with piece-wise constant inputs (step invariant), i.e., the solution obtained using lsim
is not approximative (modulu machine precision). ZoH sampling is commonly used to discretize continuous-time plant models that are to be controlled using a discrete-time controller.
FoH sampling is exact for linear systems with piece-wise linear inputs (ramp invariant), this is a good choice for simulation of systems with smooth continuous inputs.
To approximate the behavior of a continuous-time system well in the frequency domain, the :tustin
(trapezoidal / bilinear) method may be most appropriate. In this case, the pre-warping argument can be used to ensure that the frequency response of the discrete-time system matches the continuous-time system at a given frequency. The tustin transformation alters the meaning of the state components, while ZoH and FoH preserve the meaning of the state components. The Tustin method is commonly used to discretize a continuous-time controller.
The forward-Euler method generally requires the sample time to be very small relative to the time constants of the system, and its use is generally discouraged.
Classical rules-of-thumb for selecting the sample time for control design dictate that Ts
should be chosen as $0.2 ≤ ωgc⋅Ts ≤ 0.6$ where $ωgc$ is the gain-crossover frequency (rad/s).
ControlSystemsBase.c2d
— MethodQd = c2d(sys::StateSpace{Continuous}, Qc::Matrix, Ts; opt=:o)
Qd, Rd = c2d(sys::StateSpace{Continuous}, Qc::Matrix, Rc::Matrix, Ts; opt=:o)
Qd = c2d(sys::StateSpace{Discrete}, Qc::Matrix; opt=:o)
Qd, Rd = c2d(sys::StateSpace{Discrete}, Qc::Matrix, Rc::Matrix; opt=:o)
Sample a continuous-time covariance or LQR cost matrix to fit the provided discrete-time system.
If opt = :o
(default), the matrix is assumed to be a covariance matrix. The measurement covariance R
may also be provided. If opt = :c
, the matrix is instead assumed to be a cost matrix for an LQR problem.
Measurement covariance (here called Rc
) is usually estimated in discrete time, and is in this case not dependent on the sample rate. Discretization of the measurement covariance only makes sense when a continuous-time controller has been designed and the closest corresponding discrete-time controller is desired.
The method used comes from theorem 5 in the reference below.
Ref: "Discrete-time Solutions to the Continuous-time Differential Lyapunov Equation With Applications to Kalman Filtering", Patrik Axelsson and Fredrik Gustafsson
On singular covariance matrices: The traditional double integrator with covariance matrix Q = diagm([0,σ²])
can not be sampled with this method. Instead, the input matrix ("Cholesky factor") of Q
must be manually kept track of, e.g., the noise of variance σ²
enters like N = [0, 1]
which is sampled using ZoH and becomes Nd = [1/2 Ts^2; Ts]
which results in the covariance matrix σ² * Nd * Nd'
.
Example:
The following example designs a continuous-time LQR controller for a resonant system. This is simulated with OrdinaryDiffEq to allow the ODE integrator to also integrate the continuous-time LQR cost (the cost is added as an additional state variable). We then discretize both the system and the cost matrices and simulate the same thing. The discretization of an LQR contorller in this way is sometimes refered to as lqrd
.
using ControlSystemsBase, LinearAlgebra, OrdinaryDiffEq, Test
sysc = DemoSystems.resonant()
x0 = ones(sysc.nx)
Qc = [1 0.01; 0.01 2] # Continuous-time cost matrix for the state
Rc = I(1) # Continuous-time cost matrix for the input
L = lqr(sysc, Qc, Rc)
dynamics = function (xc, p, t)
x = xc[1:sysc.nx]
u = -L*x
dx = sysc.A*x + sysc.B*u
dc = dot(x, Qc, x) + dot(u, Rc, u)
return [dx; dc]
end
prob = ODEProblem(dynamics, [x0; 0], (0.0, 10.0))
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
cc = sol.u[end][end] # Continuous-time cost
# Discrete-time version
Ts = 0.01
sysd = c2d(sysc, Ts)
Qd, Rd = c2d(sysd, Qc, Rc, opt=:c)
Ld = lqr(sysd, Qd, Rd)
sold = lsim(sysd, (x, t) -> -Ld*x, 0:Ts:10, x0 = x0)
function cost(x, u, Q, R)
dot(x, Q, x) + dot(u, R, u)
end
cd = cost(sold.x, sold.u, Qd, Rd) # Discrete-time cost
@test cc ≈ cd rtol=0.01 # These should be similar
ControlSystemsBase.c2d_poly2poly
— Methodc2d_poly2poly(ro, Ts)
returns the polynomial coefficients in discrete time given polynomial coefficients in continuous time
ControlSystemsBase.c2d_roots2poly
— Methodc2d_roots2poly(ro, Ts)
returns the polynomial coefficients in discrete time given a vector of roots in continuous time
ControlSystemsBase.c2d_x0map
— Functionsysd, x0map = c2d_x0map(sys::AbstractStateSpace{<:Continuous}, Ts, method=:zoh; w_prewarp=0)
Returns the discretization sysd
of the system sys
and a matrix x0map
that transforms the initial conditions to the discrete domain by x0_discrete = x0map*[x0; u0]
See c2d
for further details.
ControlSystemsBase.comp_sensitivity
— Method ▲
│e₁
│ ┌─────┐
d₁────+──┴──► P ├─────┬──►e₄
│ └─────┘ │
│ │
│ ┌─────┐ -│
e₂◄──┴─────┤ C ◄──┬──+───d₂
└─────┘ │
│e₃
▼
input_sensitivity
is the transfer function from d₁ to e₁, (I + CP)⁻¹output_sensitivity
is the transfer function from d₂ to e₃, (I + PC)⁻¹input_comp_sensitivity
is the transfer function from d₁ to e₂, (I + CP)⁻¹CPoutput_comp_sensitivity
is the transfer function from d₂ to e₄, (I + PC)⁻¹PCG_PS
is the transfer function from d₁ to e₄, (1 + PC)⁻¹PG_CS
is the transfer function from d₂ to e₂, (1 + CP)⁻¹C
ControlSystemsBase.controllability
— Methodcontrollability(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
— MethodP = 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 (DWD' .!= 0) will equal Inf
for continuous-time systems.
ControlSystemsBase.ctrb
— Methodctrb(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.d2c
— Functiond2c(sys::AbstractStateSpace{<:Discrete}, method::Symbol = :zoh; w_prewarp=0)
Convert discrete-time system to a continuous time system, assuming that the discrete-time system was discretized using method
. Available methods are `:zoh, :fwdeuler´.
w_prewarp
: Frequency for pre-warping when using the Tustin method, has no effect for other methods.
See also d2c_exact
.
ControlSystemsBase.d2c
— FunctionQc = d2c(sys::AbstractStateSpace{<:Discrete}, Qd::AbstractMatrix; opt=:o)
Resample discrete-time covariance matrix belonging to sys
to the equivalent continuous-time matrix.
The method used comes from theorem 5 in the reference below.
If opt = :c
, the matrix is instead assumed to be a cost matrix for an LQR problem.
Ref: Discrete-time Solutions to the Continuous-time Differential Lyapunov Equation With Applications to Kalman Filtering Patrik Axelsson and Fredrik Gustafsson
ControlSystemsBase.d2c_exact
— Functiond2c_exact(sys::AbstractStateSpace{<:Discrete}, method = :causal)
Translate a discrete-time system to a continuous-time system by one of the substitutions
- $z^{-1} = e^{-sT_s}$ if
method = :causal
(default) - $z = e^{sT_s}$ if
method = :acausal
The translation is exact in the frequency domain, i.e., the frequency response of the resulting continuous-time system is identical to the frequency response of the discrete-time system.
This method of translation is useful when analyzing hybrid continuous/discrete systems in the frequency domain and high accuracy is required.
The resulting system will be be a static system in feedback with pure delays. When method = :causal
, the delays will be positive, resulting in a causal system that can be simulated in the time domain. When method = :acausal
, the delays will be negative, resulting in an acausal system that can not be simulated in the time domain. The acausal translation results in a smaller system with half as many delay elements in the feedback path.
ControlSystemsBase.dab
— MethodX,Y = dab(A,B,C)
Solves the Diophantine-Aryabhatta-Bezout identity
$AX + BY = C$, where $A, B, C, X$ and $Y$ are polynomials and $deg Y = deg A - 1$.
See Computer-Controlled Systems: Theory and Design, Third Edition Karl Johan Åström, Björn Wittenmark
ControlSystemsBase.damp
— MethodWn, zeta, ps = damp(sys)
Compute the natural frequencies, Wn
, and damping ratios, zeta
, of the poles, ps
, of sys
ControlSystemsBase.dampreport
— Methoddampreport(sys)
Display a report of the poles, damping ratio, natural frequency, and time constant of the system sys
ControlSystemsBase.dcgain
— Functiondcgain(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.delay
— Methoddelay(tau)
delay(tau, Ts)
delay(T::Type{<:Number}, tau)
delay(T::Type{<:Number}, tau, Ts)
Create a pure time delay of length τ
of type T
.
The type T
defaults to promote_type(Float64, typeof(tau))
.
If Ts
is given, the delay is discretized with sampling time Ts
and a discrete-time StateSpace object is returned.
Example:
Create a LTI system with an input delay of L
L = 1
tf(1, [1, 1])*delay(L)
s = tf("s")
tf(1, [1, 1])*exp(-s*L) # Equivalent to the version above
ControlSystemsBase.delaymargin
— Methoddₘ = 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.diagonalize
— Methoddsys = diagonalize(s::StateSpace, digits=12)
Diagonalizes the system such that the A-matrix is diagonal.
ControlSystemsBase.evalfr
— Methodevalfr(sys, x)
Evaluate the transfer function of the LTI system sys at the complex number s=x (continuous-time) or z=x (discrete-time).
For many values of x
, use freqresp
instead.
ControlSystemsBase.extended_gangoffour
— Functionextended_gangoffour(P, C, pos=true)
Returns a single statespace system that maps
w1
reference or measurement noisew2
load disturbance
to
z1
control errorz2
control input
z1 z2
▲ ┌─────┐ ▲ ┌─────┐
│ │ │ │ │ │
w1──+─┴─►│ C ├──┴───+─►│ P ├─┐
│ │ │ │ │ │ │
│ └─────┘ │ └─────┘ │
│ w2 │
└────────────────────────────┘
The returned system has the transfer-function matrix
\[\begin{bmatrix} I \\ C \end{bmatrix} (I + PC)^{-1} \begin{bmatrix} I & P \end{bmatrix}\]
or in code
# For SISO P
S = G[1, 1]
PS = G[1, 2]
CS = G[2, 1]
T = G[2, 2]
# For MIMO P
S = G[1:P.ny, 1:P.nu]
PS = G[1:P.ny, P.ny+1:end]
CS = G[P.ny+1:end, 1:P.ny]
T = G[P.ny+1:end, P.ny+1:end] # Input complimentary sensitivity function
The gang of four can be plotted like so
Gcl = extended_gangoffour(G, C) # Form closed-loop system
bodeplot(Gcl, lab=["S" "CS" "PS" "T"], plotphase=false) |> display # Plot gang of four
Note, the last input of Gcl is the negative of the PS
and T
transfer functions from gangoffour2
. To get a transfer matrix with the same sign as G_PS
and input_comp_sensitivity
, call extended_gangoffour(P, C, pos=false)
. See glover_mcfarlane
from RobustAndOptimalControl.jl for an extended example. See also ncfmargin
and feedback_control
from RobustAndOptimalControl.jl.
ControlSystemsBase.feedback
— Methodfeedback(sys1::AbstractStateSpace, sys2::AbstractStateSpace;
U1=:, Y1=:, U2=:, Y2=:, W1=:, Z1=:, W2=Int[], Z2=Int[],
Wperm=:, Zperm=:, pos_feedback::Bool=false)
Basic usefeedback(sys1, sys2)
forms the (negative) feedback interconnection
┌──────────────┐
◄──────────┤ sys1 │◄──── Σ ◄──────
│ │ │ │
│ └──────────────┘ -1
│ |
│ ┌──────────────┐ │
└─────►│ sys2 ├──────┘
│ │
└──────────────┘
If no second system sys2
is given, negative identity feedback (sys2 = 1
) is assumed. The returned closed-loop system will have a state vector comprised of the state of sys1
followed by the state of sys2
.
Advanced usefeedback
also supports more flexible use according to the figure below
┌──────────────┐
z1◄─────┤ sys1 │◄──────w1
┌─── y1◄─────┤ │◄──────u1 ◄─┐
│ └──────────────┘ │
│ α
│ ┌──────────────┐ │
└──► u2─────►│ sys2 ├───────►y2──┘
w2─────►│ ├───────►z2
└──────────────┘
U1
, W1
specifies the indices of the input signals of sys1
corresponding to u1
and w1
Y1
, Z1
specifies the indices of the output signals of sys1
corresponding to y1
and z1
U2
, W2
, Y2
, Z2
specifies the corresponding signals of sys2
Specify Wperm
and Zperm
to reorder the inputs (corresponding to [w1; w2]) and outputs (corresponding to [z1; z2]) in the resulting statespace model.
Negative feedback (α = -1) is the default. Specify pos_feedback=true
for positive feedback (α = 1).
See also lft
, starprod
, sensitivity
, input_sensitivity
, output_sensitivity
, comp_sensitivity
, input_comp_sensitivity
, output_comp_sensitivity
, G_PS
, G_CS
.
The manual section From block diagrams to code contains higher-level instructions on how to use this function.
See Zhou, Doyle, Glover (1996) for similar (somewhat less symmetric) formulas.
ControlSystemsBase.feedback
— Methodfeedback(sys)
feedback(sys1, sys2)
For a general LTI-system, feedback
forms the negative feedback interconnection
>-+ sys1 +-->
| |
(-)sys2 +
If no second system is given, negative identity feedback is assumed
ControlSystemsBase.feedback2dof
— Methodfeedback2dof(P,R,S,T)
feedback2dof(B,A,R,S,T)
- Return
BT/(AR+ST)
where B and A are the numerator and denominator polynomials ofP
respectively - Return
BT/(AR+ST)
ControlSystemsBase.feedback2dof
— Methodfeedback2dof(P::TransferFunction, C::TransferFunction, F::TransferFunction)
Return the transfer function P(F+C)/(1+PC)
which is the closed-loop system with process P
, controller C
and feedforward filter F
from reference to control signal (by-passing C
).
+-------+
| |
+-----> F +----+
| | | |
| +-------+ |
| +-------+ | +-------+
r | - | | | | | y
+--+-----> C +----+----> P +---+-->
| | | | | |
| +-------+ +-------+ |
| |
+--------------------------------+
ControlSystemsBase.freqresp!
— Methodfreqresp!(R::Array{T, 3}, sys::LTISystem, w_vec::AbstractVector{<:Real})
In-place version of freqresp
that takes a pre-allocated array R
of size (ny, nu, nw)`
ControlSystemsBase.freqresp
— Methodsys_fr = freqresp(sys, w)
Evaluate the frequency response of a linear system
w -> C*((iw*im*I - A)^-1)*B + D
of system sys
over the frequency vector w
.
ControlSystemsBase.freqrespv
— Methodfreqrespv(G::AbstractMatrix, w_vec::AbstractVector{<:Real}; )
For use with SISO systems where it acts the same as freqresp
but with the extra dimensions removed in the returned values.
ControlSystemsBase.freqrespv
— Methodfreqrespv(G::Number, w_vec::AbstractVector{<:Real}; )
For use with SISO systems where it acts the same as freqresp
but with the extra dimensions removed in the returned values.
ControlSystemsBase.freqrespv
— Methodfreqrespv(sys::LTISystem, w_vec::AbstractVector{W}; )
For use with SISO systems where it acts the same as freqresp
but with the extra dimensions removed in the returned values.
ControlSystemsBase.gangoffour
— MethodS, 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 functionPS = (1+PC)\P
Load disturbance to measurement signalCS = (1+PC)\C
Measurement noise to control signalT = PC/(1+PC)
Complementary sensitivity function
If minimal=true
, minreal
will be applied to all transfer functions.
ControlSystemsBase.gangoffourplot
— Methodfig = gangoffourplot(P::LTISystem, C::LTISystem; minimal=true, plotphase=false, Ms_lines = [1.0, 1.25, 1.5], Mt_lines = [], sigma = true, kwargs...)
Gang-of-Four plot.
sigma
determines whether a sigmaplot
is used instead of a bodeplot
for MIMO S
and T
. kwargs
are sent as argument to RecipesBase.plot.
ControlSystemsBase.gangofseven
— MethodS, 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 functionPS = P/(1+PC)
CS = C/(1+PC)
T = PC/(1+PC)
Complementary sensitivity functionRY = PCF/(1+PC)
RU = CF/(1+P*C)
RE = F/(1+P*C)
ControlSystemsBase.gram
— Methodgram(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
— MethodU = 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
— MethodNinf, ω_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.impulse
— Methody, t, x = impulse(sys[, tfinal])
y, t, x = impulse(sys[, t])
Calculate the response of the system sys
to an impulse at time t = 0
. For continous-time systems, the impulse is a unit Dirac impulse. For discrete-time systems, the impulse lasts one sample and has magnitude 1/Ts
. If the final time tfinal
or time vector t
is not provided, one is calculated based on the system pole locations.
The return value is a structure of type SimResult
. A SimResul
can be plotted by plot(result)
, or destructured as y, t, x = result
.
y
has size (ny, length(t), nu)
, x
has size (nx, length(t), nu)
See also lsim
.
ControlSystemsBase.innovation_form
— Methodsysi = 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
— Methodsysi = 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.input_comp_sensitivity
— Methodinput_comp_sensitivity(P,C)
Transfer function from load disturbance to control signal.
- "Input" signifies that the transfer function is from the input of the plant.
- "Complimentary" signifies that the transfer function is to an output (in this case controller output)
▲
│e₁
│ ┌─────┐
d₁────+──┴──► P ├─────┬──►e₄
│ └─────┘ │
│ │
│ ┌─────┐ -│
e₂◄──┴─────┤ C ◄──┬──+───d₂
└─────┘ │
│e₃
▼
input_sensitivity
is the transfer function from d₁ to e₁, (I + CP)⁻¹output_sensitivity
is the transfer function from d₂ to e₃, (I + PC)⁻¹input_comp_sensitivity
is the transfer function from d₁ to e₂, (I + CP)⁻¹CPoutput_comp_sensitivity
is the transfer function from d₂ to e₄, (I + PC)⁻¹PCG_PS
is the transfer function from d₁ to e₄, (1 + PC)⁻¹PG_CS
is the transfer function from d₂ to e₂, (1 + CP)⁻¹C
ControlSystemsBase.input_names
— Methodinput_names(P)
input_names(P, i)
Get a vector of strings with the names of the inputs of P
, or the i
:th name if an index is given.
ControlSystemsBase.input_sensitivity
— Methodinput_sensitivity(P, C)
Transfer function from load disturbance to total plant input.
- "Input" signifies that the transfer function is from the input of the plant.
▲
│e₁
│ ┌─────┐
d₁────+──┴──► P ├─────┬──►e₄
│ └─────┘ │
│ │
│ ┌─────┐ -│
e₂◄──┴─────┤ C ◄──┬──+───d₂
└─────┘ │
│e₃
▼
input_sensitivity
is the transfer function from d₁ to e₁, (I + CP)⁻¹output_sensitivity
is the transfer function from d₂ to e₃, (I + PC)⁻¹input_comp_sensitivity
is the transfer function from d₁ to e₂, (I + CP)⁻¹CPoutput_comp_sensitivity
is the transfer function from d₂ to e₄, (I + PC)⁻¹PCG_PS
is the transfer function from d₁ to e₄, (1 + PC)⁻¹PG_CS
is the transfer function from d₂ to e₂, (1 + CP)⁻¹C
ControlSystemsBase.iscontinuous
— Methodiscontinuous(sys)
Returns true
for a continuous-time system sys
, else returns false
.
ControlSystemsBase.isdiscrete
— Methodisdiscrete(sys)
Returns true
for a discrete-time system sys
, else returns false
.
ControlSystemsBase.isproper
— Methodisproper(tf)
Returns true
if the TransferFunction
is proper. This means that order(den)
= order(num))
ControlSystemsBase.isstable
— Methodisstable(sys)
Returns true
if sys
is stable, else returns false
.
ControlSystemsBase.kalman
— Methodkalman(Continuous, A, C, R1, R2)
kalman(Discrete, A, C, R1, R2; direct = false)
kalman(sys, R1, R2; direct = false)
Calculate the optimal asymptotic Kalman gain.
If direct = true
, the observer gain is computed for the pair (A, CA)
instead of (A,C)
. This option is intended to be used together with the option direct = true
to observer_controller
. Ref: "Computer-Controlled Systems" pp 140. direct = false
is sometimes referred to as a "delayed" estimator, while direct = true
is a "current" estimator.
To obtain a discrete-time approximation to a continuous-time LQG problem, the function c2d
can be used to obtain corresponding discrete-time covariance matrices.
To obtain an LTISystem that represents the Kalman filter, pass the obtained Kalman feedback gain into observer_filter
. To obtain an LQG controller, pass the obtained Kalman feedback gain as well as a state-feedback gain computed using lqr
into observer_controller
.
The args...; kwargs...
are sent to the Riccati solver, allowing specification of cross-covariance etc. See ?MatrixEquations.arec/ared
for more help.
ControlSystemsBase.laglink
— Methodlaglink(a, M; [Ts])
Returns a phase retarding link, the rule of thumb a = 0.1ωc
guarantees less than 6 degrees phase margin loss. The bode curve will go from M
, bend down at a/M
and level out at 1 for frequencies > a
\[\dfrac{s + a}{s + a/M} = M \dfrac{1 + s/a}{1 + sM/a}\]
ControlSystemsBase.leadlink
— Functionleadlink(b, N, K=1; [Ts])
Returns a phase advancing link, the top of the phase curve is located at ω = b√(N)
where the link amplification is K√(N)
The bode curve will go from K
, bend up at b
and level out at KN
for frequencies > bN
The phase advance at ω = b√(N)
can be plotted as a function of N
with leadlinkcurve()
Values of N < 1
will give a phase retarding link.
\[KN \dfrac{s + b}{s + bN} = K \dfrac{1 + s/b}{1 + s/(bN)}\]
See also leadlinkat
laglink
ControlSystemsBase.leadlinkat
— Functionleadlinkat(ω, N, K=1; [Ts])
Returns a phase advancing link, the top of the phase curve is located at ω
where the link amplification is K√(N)
The bode curve will go from K
, bend up at ω/√(N)
and level out at KN
for frequencies > ω√(N)
The phase advance at ω
can be plotted as a function of N
with leadlinkcurve()
Values of N < 1
will give a phase retarding link.
See also leadlink
laglink
ControlSystemsBase.leadlinkcurve
— Functionleadlinkcurve(start=1)
Plot the phase advance as a function of N
for a lead link (phase advance link) If an input argument start
is given, the curve is plotted from start
to 10, else from 1 to 10.
See also leadlink, leadlinkat
ControlSystemsBase.lft
— Functionlft(G, Δ, type=:l)
Lower and upper linear fractional transformation between systems G
and Δ
.
Specify :l
lor lower LFT, and :u
for upper LFT.
G
must have more inputs and outputs than Δ
has outputs and inputs.
For details, see Chapter 9.1 in Zhou, K. and JC Doyle. Essentials of robust control, Prentice hall (NJ), 1998
ControlSystemsBase.linfnorm
— MethodNinf, ω_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.loopshapingPI
— MethodC, kp, ki, fig, CF = loopshapingPI(P, ωp; ϕl, rl, phasemargin, form=:standard, doplot=false, Tf, F)
Selects the parameters of a PI-controller (on parallel form) such that the Nyquist curve of P
at the frequency ωp
is moved to rl exp(i ϕl)
The parameters can be returned as one of several common representations chosen by form
, the options are
:standard
- $K_p(1 + 1/(T_i s) + T_d s)$:series
- $K_c(1 + 1/(τ_i s))(τ_d s + 1)$:parallel
- $K_p + K_i/s + K_d s$
If phasemargin
is supplied (in degrees), ϕl
is selected such that the curve is moved to an angle of phasemargin - 180
degrees
If no rl
is given, the magnitude of the curve at ωp
is kept the same and only the phase is affected, the same goes for ϕl
if no phasemargin is given.
Tf
: An optional time constant for second-order measurement noise filter on the formtf(1, [Tf^2, 2*Tf/sqrt(2), 1])
to make the controller strictly proper.F
: A pre-designed filter to use instead of the default second-order filter that is used ifTf
is given.doplot
plot thegangoffourplot
andnyquistplot
of the system.
See also loopshapingPID
, pidplots
, stabregionPID
and placePI
.
ControlSystemsBase.loopshapingPID
— MethodC, kp, ki, kd, fig, CF = loopshapingPID(P, ω; Mt = 1.3, ϕt=75, form=:standard, doplot=false, lb=-10, ub=10, Tf = 1/1000ω, F = nothing)
Selects the parameters of a PID-controller such that the Nyquist curve of the loop-transfer function $L = PC$ at the frequency ω
is tangent to the circle where the magnitude of $T = PC / (1+PC)$ equals Mt
. ϕt
denotes the positive angle in degrees between the real axis and the tangent point.
The default values for Mt
and ϕt
are chosen to give a good design for processes with inertia, and may need tuning for simpler processes.
The gain of the resulting controller is generally increasing with increasing ω
and Mt
.
Arguments:
P
: A SISO plant.ω
: The specification frequency.Mt
: The magnitude of the complementary sensitivity function at the specification frequency, $|T(iω)|$.ϕt
: The positive angle in degrees between the real axis and the tangent point.doplot
: If true, gang of four and Nyquist plots will be returned infig
.lb
: log10 of lower bound forkd
.ub
: log10 of upper bound forkd
.Tf
: Time constant for second-order measurement noise filter on the formtf(1, [Tf^2, 2*Tf/sqrt(2), 1])
to make the controller strictly proper. A practical controller typically sets this time constant slower than the default, e.g.,Tf = 1/100ω
orTf = 1/10ω
F
: A pre-designed filter to use instead of the default second-order filter.
The parameters can be returned as one of several common representations chosen by form
, the options are
:standard
- $K_p(1 + 1/(T_i s) + T_ds)$:series
- $K_c(1 + 1/(τ_i s))(τ_d s + 1)$:parallel
- $K_p + K_i/s + K_d s$
See also loopshapingPI
, pidplots
, stabregionPID
and placePI
.
Example:
P = tf(1, [1,0,0]) # A double integrator
Mt = 1.3 # Maximum magnitude of complementary sensitivity
ω = 1 # Frequency at which the specification holds
C, kp, ki, kd, fig, CF = loopshapingPID(P, ω; Mt, ϕt = 75, doplot=true)
ControlSystemsBase.lqr
— Methodlqr(sys, Q, R)
lqr(Continuous, A, B, Q, R, args...; kwargs...)
lqr(Discrete, A, B, Q, R, args...; kwargs...)
Calculate the optimal gain matrix K
for the state-feedback law u = -K*x
that minimizes the cost function:
J = integral(x'Qx + u'Ru, 0, inf) for the continuous-time model dx = Ax + Bu
. J = sum(x'Qx + u'Ru, 0, inf) for the discrete-time model x[k+1] = Ax[k] + Bu[k]
.
Solve the LQR problem for state-space system sys
. Works for both discrete and continuous time systems.
The args...; kwargs...
are sent to the Riccati solver, allowing specification of cross-covariance etc. See ?MatrixEquations.arec / ared
for more help.
To obtain also the solution to the Riccati equation and the eigenvalues of the closed-loop system as well, call ControlSystemsBase.MatrixEquations.arec / ared
instead (note the different order of the arguments to these functions).
To obtain a discrete-time approximation to a continuous-time LQR problem, the function c2d
can be used to obtain corresponding discrete-time cost matrices.
Examples
Continuous time
using LinearAlgebra # For identity matrix I
using Plots
A = [0 1; 0 0]
B = [0; 1]
C = [1 0]
sys = ss(A,B,C,0)
Q = I
R = I
L = lqr(sys,Q,R) # lqr(Continuous,A,B,Q,R) can also be used
u(x,t) = -L*x # Form control law,
t=0:0.1:5
x0 = [1,0]
y, t, x, uout = lsim(sys,u,t,x0=x0)
plot(t,x', lab=["Position" "Velocity"], xlabel="Time [s]")
Discrete time
using LinearAlgebra # For identity matrix I
using Plots
Ts = 0.1
A = [1 Ts; 0 1]
B = [0;1]
C = [1 0]
sys = ss(A, B, C, 0, Ts)
Q = I
R = I
L = lqr(Discrete, A,B,Q,R) # lqr(sys,Q,R) can also be used
u(x,t) = -L*x # Form control law,
t=0:Ts:5
x0 = [1,0]
y, t, x, uout = lsim(sys,u,t,x0=x0)
plot(t,x', lab=["Position" "Velocity"], xlabel="Time [s]")
ControlSystemsBase.lsim!
— Methodres = lsim!(ws::LsimWorkspace, sys::AbstractStateSpace{<:Discrete}, u, [t]; x0)
In-place version of lsim
that takes a workspace object created by calling LsimWorkspace
. Notice, if u
is a function, res.u === ws.u
. If u
is an array, res.u === u
.
ControlSystemsBase.lsim
— Methodresult = lsim(sys, u[, t]; x0, method])
result = lsim(sys, u::Function, t; x0, method)
Calculate the time response of system sys
to input u
. If x0
is omitted, a zero vector is used.
The result structure contains the fields y, t, x, u
and can be destructured automatically by iteration, e.g.,
y, t, x, u = result
result::SimResult
can also be plotted directly:
plot(result, plotu=true, plotx=false)
y
, x
, u
have time in the second dimension. Initial state x0
defaults to zero.
Continuous-time systems are simulated using an ODE solver if u
is a function (requires using ControlSystems). If u
is an array, the system is discretized (with method=:zoh
by default) before simulation. For a lower-level interface, see ?Simulator
and ?solve
. For continuous-time systems, keyword arguments are forwarded to the ODE solver. By default, the option dtmax = t[2]-t[1]
is used to prevent the solver from stepping over discontinuities in u(x, t)
. This prevents the solver from taking too large steps, but may also slow down the simulation when u
is smooth. To disable this behavior, set dtmax = Inf
.
u
can be a function or a matrix of precalculated control signals and must have dimensions (nu, length(t))
. If u
is a function, then u(x,i)
(for discrete systems) or u(x,t)
(for continuous ones) is called to calculate the control signal at every iteration (time instance used by solver). This can be used to provide a control law such as state feedback u(x,t) = -L*x
calculated by lqr
. To simulate a unit step at t=t₀
, use (x,t)-> t ≥ t₀
, for a ramp, use (x,t)-> t
, for a step at t=5
, use (x,t)-> (t >= 5)
etc.
Note: The function u
will be called once before simulating to verify that it returns an array of the correct dimensions. This can cause problems if u
is stateful. You can disable this check by passing check_u = false
.
For maximum performance, see function lsim!
, available for discrete-time systems only.
Usage example:
using ControlSystems
using LinearAlgebra: I
using Plots
A = [0 1; 0 0]
B = [0;1]
C = [1 0]
sys = ss(A,B,C,0)
Q = I
R = I
L = lqr(sys,Q,R)
u(x,t) = -L*x # Form control law
t = 0:0.1:5
x0 = [1,0]
y, t, x, uout = lsim(sys,u,t,x0=x0)
plot(t,x', lab=["Position" "Velocity"], xlabel="Time [s]")
# Alternative way of plotting
res = lsim(sys,u,t,x0=x0)
plot(res)
ControlSystemsBase.margin
— Methodwgm, 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.marginplot
— Functionfig = marginplot(sys::LTISystem [,w::AbstractVector]; balance=true, kwargs...)
marginplot(sys::Vector{LTISystem}, w::AbstractVector; balance=true, kwargs...)
Plot all the amplitude and phase margins of the system(s) sys
.
- A frequency vector
w
can be optionally provided. balance
: Callbalance_statespace
on the system before plotting.
kwargs
is sent as argument to RecipesBase.plot.
ControlSystemsBase.markovparam
— Methodmarkovparam(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.minreal
— Functionminreal(tf::TransferFunction, eps=sqrt(eps()))
Create a minimal representation of each transfer function in tf
by cancelling poles and zeros will promote system to an appropriate numeric type
ControlSystemsBase.minreal
— Methodminreal(sys::T; fast=false, kwargs...)
Minimal realisation algorithm from P. Van Dooreen, The generalized eigenstructure problem in linear system theory, IEEE Transactions on Automatic Control
For information about the options, see ?ControlSystemsBase.MatrixPencils.lsminreal
See also sminreal
, which is both numerically exact and substantially faster than minreal
, but with a much more limited potential in removing non-minimal dynamics.
ControlSystemsBase.nicholsplot
— Functionfig = nicholsplot{T<:LTISystem}(systems::Vector{T}, w::AbstractVector; kwargs...)
Create a Nichols plot of the LTISystem
(s). A frequency vector w
can be optionally provided.
Keyword arguments:
text = true
Gains = [12, 6, 3, 1, 0.5, -0.5, -1, -3, -6, -10, -20, -40, -60]
pInc = 30
sat = 0.4
val = 0.85
fontsize = 10
pInc
determines the increment in degrees between phase lines.
sat
∈ [0,1] determines the saturation of the gain lines
val
∈ [0,1] determines the brightness of the gain lines
Additional keyword arguments are sent to the function plotting the systems and can be used to specify colors, line styles etc. using regular RecipesBase.jl syntax
This function is based on code subject to the two-clause BSD licence Copyright 2011 Will Robertson Copyright 2011 Philipp Allgeuer
ControlSystemsBase.nonlinearity
— Methodnonlinearity(f)
nonlinearity(T, f)
Create a pure nonlinearity. f
is assumed to be a static (no memory) nonlinear function from $f : R -> R$.
The type T
defaults to Float64
.
NOTE: The nonlinear functionality in ControlSystemsBase.jl is currently experimental and subject to breaking changes not respecting semantic versioning. Use at your own risk.
Example:
Create a LTI system with a static input nonlinearity that saturates the input to [-1,1].
tf(1, [1, 1])*nonlinearity(x->clamp(x, -1, 1))
See also predefined nonlinearities saturation
, offset
.
Note: when composing linear systems with nonlinearities, it's often important to handle operating points correctly. See ControlSystemsBase.offset
for handling operating points.
ControlSystemsBase.nyquist
— Methodre, im, w = nyquist(sys[, w])
Compute the real and imaginary parts of the frequency response of system sys
at frequencies w
. See also nyquistplot
re
and im
has size (ny, nu, length(w))
ControlSystemsBase.nyquistplot
— Functionfig = nyquistplot(sys; Ms_circles=Float64[], Mt_circles=Float64[], unit_circle=false, hz=false, critical_point=-1, kwargs...)
nyquistplot(LTISystem[sys1, sys2...]; Ms_circles=Float64[], Mt_circles=Float64[], unit_circle=false, hz=false, critical_point=-1, kwargs...)
Create a Nyquist plot of the LTISystem
(s). A frequency vector w
can be optionally provided.
unit_circle
: if the unit circle should be displayed. The Nyquist curve crosses the unit circle at the gain crossover frequency.Ms_circles
: draw circles corresponding to given levels of sensitivity (circles around -1 with radii1/Ms
).Ms_circles
can be supplied as a number or a vector of numbers. A design staying outside such a circle has a phase margin of at least2asin(1/(2Ms))
rad and a gain margin of at leastMs/(Ms-1)
.Mt_circles
: draw circles corresponding to given levels of complementary sensitivity.Mt_circles
can be supplied as a number or a vector of numbers.critical_point
: point on real axis to mark as critical for encirclements- If
hz=true
, the hover information will be displayed in Hertz, the input frequency vector is still treated as rad/s. balance
: Callbalance_statespace
on the system before plotting.
kwargs
is sent as argument to plot.
ControlSystemsBase.nyquistv
— Methodnyquistv(sys::LTISystem, w::AbstractVector; )
For use with SISO systems where it acts the same as nyquist
but with the extra dimensions removed in the returned values.
ControlSystemsBase.nyquistv
— Methodnyquistv(sys::LTISystem; )
For use with SISO systems where it acts the same as nyquist
but with the extra dimensions removed in the returned values.
ControlSystemsBase.observability
— Methodobservability(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
— Methodcont = 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 systemL
: State-feedback gainu = -Lx
K
: Observer gain
See also observer_predictor
and innovation_form
.
ControlSystemsBase.observer_filter
— Methodobserver_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
— Methodobserver_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
— Functionobsv(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.output_comp_sensitivity
— Methodoutput_comp_sensitivity(P,C)
Transfer function from measurement noise / reference to plant output.
- "output" signifies that the transfer function is from the output of the plant.
- "Complimentary" signifies that the transfer function is to an output (in this case plant output)
▲
│e₁
│ ┌─────┐
d₁────+──┴──► P ├─────┬──►e₄
│ └─────┘ │
│ │
│ ┌─────┐ -│
e₂◄──┴─────┤ C ◄──┬──+───d₂
└─────┘ │
│e₃
▼
input_sensitivity
is the transfer function from d₁ to e₁, (I + CP)⁻¹output_sensitivity
is the transfer function from d₂ to e₃, (I + PC)⁻¹input_comp_sensitivity
is the transfer function from d₁ to e₂, (I + CP)⁻¹CPoutput_comp_sensitivity
is the transfer function from d₂ to e₄, (I + PC)⁻¹PCG_PS
is the transfer function from d₁ to e₄, (1 + PC)⁻¹PG_CS
is the transfer function from d₂ to e₂, (1 + CP)⁻¹C
ControlSystemsBase.output_names
— Methodoutput_names(P)
output_names(P, i)
Get a vector of strings with the names of the outputs of P
, or the i
:th name if an index is given.
ControlSystemsBase.output_sensitivity
— Methodoutput_sensitivity(P, C)
Transfer function from measurement noise / reference to control error.
- "output" signifies that the transfer function is from the output of the plant.
▲
│e₁
│ ┌─────┐
d₁────+──┴──► P ├─────┬──►e₄
│ └─────┘ │
│ │
│ ┌─────┐ -│
e₂◄──┴─────┤ C ◄──┬──+───d₂
└─────┘ │
│e₃
▼
input_sensitivity
is the transfer function from d₁ to e₁, (I + CP)⁻¹output_sensitivity
is the transfer function from d₂ to e₃, (I + PC)⁻¹input_comp_sensitivity
is the transfer function from d₁ to e₂, (I + CP)⁻¹CPoutput_comp_sensitivity
is the transfer function from d₂ to e₄, (I + PC)⁻¹PCG_PS
is the transfer function from d₁ to e₄, (1 + PC)⁻¹PG_CS
is the transfer function from d₂ to e₂, (1 + CP)⁻¹C
ControlSystemsBase.pade
— Methodpade(G::DelayLtiSystem, N)
Approximate all time-delays in G
by Padé approximations of degree N
.
ControlSystemsBase.pade
— Methodpade(τ::Real, N::Int)
Compute the N
th order Padé approximation of a time-delay of length τ
.
See also thiran
for discretization of delays.
ControlSystemsBase.parallel
— Methodparallel(sys1::LTISystem, sys2::LTISystem)
Connect systems in parallel, equivalent to sys2+sys1
ControlSystemsBase.pid
— FunctionC = pid(param_p, param_i, [param_d]; form=:standard, state_space=false, [Tf], [Ts])
Calculates and returns a PID controller.
The form
can be chosen as one of the following
:standard
-Kp*(1 + 1/(Ti*s) + Td*s)
:series
-Kc*(1 + 1/(τi*s))*(τd*s + 1)
:parallel
-Kp + Ki/s + Kd*s
If state_space
is set to true
, either Kd
has to be zero or a positive Tf
has to be provided for creating a filter on the input to allow for a state space realization. The filter used is 1 / (1 + s*Tf + (s*Tf)^2/2)
, where Tf
can typically be chosen as Ti/N
for a PI controller and Td/N
for a PID controller, and N
is commonly in the range 2 to 20. The state space will be returned on controllable canonical form.
For a discrete controller a positive Ts
can be supplied. In this case, the continuous-time controller is discretized using the Tustin method.
Examples
C1 = pid(3.3, 1, 2) # Kd≠0 works without filter in tf form
C2 = pid(3.3, 1, 2; Tf=0.3, state_space=true) # In statespace a filter is needed
C3 = pid(2., 3, 0; Ts=0.4, state_space=true) # Discrete
The functions pid_tf
and pid_ss
are also exported. They take the same parameters and is what is actually called in pid
based on the state_space
parameter.
ControlSystemsBase.pidplots
— Methodpidplots(P, args...; params_p, params_i, params_d=0, form=:standard, ω=0, grid=false, kwargs...)
Display the relevant plots related to closing the loop around process P
with a PID controller supplied in params
on one of the following forms:
:standard
-Kp*(1 + 1/(Ti*s) + Td*s)
:series
-Kc*(1 + 1/(τi*s))*(τd*s + 1)
:parallel
-Kp + Ki/s + Kd*s
The sent in values can be arrays to evaluate multiple different controllers, and if grid=true
it will be a grid search over all possible combinations of the values.
Available plots are :gof
for Gang of four, :nyquist
, :controller
for a bode plot of the controller TF and :pz
for pole-zero maps and should be supplied as additional arguments to the function.
One can also supply a frequency vector ω
to be used in Bode and Nyquist plots.
See also loopshapingPI
, stabregionPID
ControlSystemsBase.place
— Functionplace(A, B, p, opt=:c; direct = false)
place(sys::StateSpace, p, opt=:c; direct = false)
Calculate the gain matrix K
such that A - BK
has eigenvalues p
.
place(A, C, p, opt=:o)
place(sys::StateSpace, p, opt=:o)
Calculate the observer gain matrix L
such that A - LC
has eigenvalues p
.
If direct = true
and opt = :o
, the the observer gain K
is calculated such that A - KCA
has eigenvalues p
, this option is to be used together with direct = true
in observer_controller
.
Note: only apply direct = true
to discrete-time systems.
Ref: "Computer-Controlled Systems" pp 140.
Uses Ackermann's formula for SISO systems and place_knvd
for MIMO systems.
Please note that this function can be numerically sensitive, solving the placement problem in extended precision might be beneficial.
ControlSystemsBase.placePI
— MethodC, kp, ki = placePI(P, ω₀, ζ; form=:standard)
Selects the parameters of a PI-controller such that the poles of closed loop between P
and C
are placed to match the poles of s^2 + 2ζω₀s + ω₀^2
.
The parameters can be returned as one of several common representations chose by form
, the options are
:standard
- $K_p(1 + 1/(T_i s))$:series
- $K_c(1 + 1/(τ_i s))$ (equivalent to above for PI controllers):parallel
- $K_p + K_i/s$
C
is the returned transfer function of the controller and params
is a named tuple containing the parameters. The parameters can be accessed as params.Kp
or params["Kp"]
from the named tuple, or they can be unpacked using Kp, Ti, Td = values(params)
.
See also loopshapingPI
ControlSystemsBase.place_knvd
— Methodplace_knvd(A::AbstractMatrix, B, λ; verbose = false, init = :s)
Robust pole placement using the algorithm from
"Robust Pole Assignment in Linear State Feedback", Kautsky, Nichols, Van Dooren
This implementation uses "method 0" for the X-step and the QR factorization for all factorizations.
This function will be called automatically when place
is called with a MIMO system.
Arguments:
init
: Determines the initialization strategy for the iterations for find theX
matrix. Possible choices are:id
(default),:rand
,:s
.
ControlSystemsBase.plyap
— MethodXc = 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.poles
— Methodpoles(sys)
Compute the poles of system sys
.
ControlSystemsBase.pzmap
— Functionfig = pzmap(fig, system, args...; hz = false, kwargs...)
Create a pole-zero map of the LTISystem
(s) in figure fig
, args
and kwargs
will be sent to the scatter
plot command.
To customize the unit-circle drawn for discrete systems, modify the line attributes, e.g., linecolor=:red
.
If hz
is true, all poles and zeros are scaled by 1/2π.
ControlSystemsBase.reduce_sys
— Methodreduce_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
— Methodrelative_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
— Methodrelative_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.rgaplot
— Functionrgaplot(sys, args...; hz=false)
rgaplot(LTISystem[sys1, sys2...], args...; hz=false, balance=true)
Plot the relative-gain array entries of the LTISystem
(s). 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. balance
: Callbalance_statespace
on the system before plotting.
kwargs
is sent as argument to Plots.plot.
ControlSystemsBase.rstc
— MethodSee ?rstd
for the discrete case
ControlSystemsBase.rstd
— MethodR,S,T = rstd(BPLUS,BMINUS,A,BM1,AM,AO,AR,AS)
R,S,T = rstd(BPLUS,BMINUS,A,BM1,AM,AO,AR)
R,S,T = rstd(BPLUS,BMINUS,A,BM1,AM,AO)
Polynomial synthesis in discrete time.
Polynomial synthesis according to "Computer-Controlled Systems" ch 10 to design a controller $R(q) u(k) = T(q) r(k) - S(q) y(k)$
Inputs:
BPLUS
: Part of open loop numeratorBMINUS
: Part of open loop numeratorA
: Open loop denominatorBM1
: Additional zerosAM
: Closed loop denominatorAO
: Observer polynomialAR
: Pre-specified factor of R,
e.g integral part [1, -1]^k
AS
: Pre-specified factor of S,
e.g notch filter [1, 0, w^2]
Outputs: R,S,T
: Polynomials in controller
See function dab
how the solution to the Diophantine- Aryabhatta-Bezout identity is chosen.
See Computer-Controlled Systems: Theory and Design, Third Edition Karl Johan Åström, Björn Wittenmark
ControlSystemsBase.sensitivity
— MethodThe output sensitivity function$S_o = (I + PC)^{-1}$ is the transfer function from a reference input to control error, while the input sensitivity function $S_i = (I + CP)^{-1}$ is the transfer function from a disturbance at the plant input to the total plant input. For SISO systems, input and output sensitivity functions are equal. In general, we want to minimize the sensitivity function to improve robustness and performance, but practical constraints always cause the sensitivity function to tend to 1 for high frequencies. A robust design minimizes the peak of the sensitivity function, $M_S$. The peak magnitude of $S$ is the inverse of the distance between the open-loop Nyquist curve and the critical point -1. Upper bounding the sensitivity peak $M_S$ gives lower-bounds on phase and gain margins according to
\[ϕ_m ≥ 2\text{sin}^{-1}(\frac{1}{2M_S}), g_m ≥ \frac{M_S}{M_S-1}\]
Generally, bounding $M_S$ is a better objective than looking at gain and phase margins due to the possibility of combined gain and pahse variations, which may lead to poor robustness despite large gain and pahse margins.
▲
│e₁
│ ┌─────┐
d₁────+──┴──► P ├─────┬──►e₄
│ └─────┘ │
│ │
│ ┌─────┐ -│
e₂◄──┴─────┤ C ◄──┬──+───d₂
└─────┘ │
│e₃
▼
input_sensitivity
is the transfer function from d₁ to e₁, (I + CP)⁻¹output_sensitivity
is the transfer function from d₂ to e₃, (I + PC)⁻¹input_comp_sensitivity
is the transfer function from d₁ to e₂, (I + CP)⁻¹CPoutput_comp_sensitivity
is the transfer function from d₂ to e₄, (I + PC)⁻¹PCG_PS
is the transfer function from d₁ to e₄, (1 + PC)⁻¹PG_CS
is the transfer function from d₂ to e₂, (1 + CP)⁻¹C
ControlSystemsBase.series
— Methodseries(sys1::LTISystem, sys2::LTISystem)
Connect systems in series, equivalent to sys2*sys1
ControlSystemsBase.setPlotScale
— MethodsetPlotScale(str)
Set the default scale of magnitude in bodeplot
and sigmaplot
. str
should be either "dB"
or "log10"
. The default scale if none is chosen is "log10"
.
ControlSystemsBase.sigma
— Methodsv, w = sigma(sys[, w])
Compute the singular values sv
of the frequency response of system sys
at frequencies w
. See also sigmaplot
sv
has size (min(ny, nu), length(w))
ControlSystemsBase.sigmaplot
— Functionsigmaplot(sys, args...; hz=false balance=true, extrema)
sigmaplot(LTISystem[sys1, sys2...], args...; hz=false, balance=true, extrema)
Plot the singular values of the frequency response of the LTISystem
(s). 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. balance
: Callbalance_statespace
on the system before plotting.extrema
: Only plot the largest and smallest singular values.
kwargs
is sent as argument to Plots.plot.
ControlSystemsBase.sigmav
— Methodsigmav(sys::LTISystem, w::AbstractVector; )
For use with SISO systems where it acts the same as sigma
but with the extra dimensions removed in the returned values.
ControlSystemsBase.sigmav
— Methodsigmav(sys::LTISystem; )
For use with SISO systems where it acts the same as sigma
but with the extra dimensions removed in the returned values.
ControlSystemsBase.similarity_transform
— Methodsyst = 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.sminreal
— Methodsminreal(sys)
Compute the structurally minimal realization of the state-space system sys
. A structurally minimal realization is one where only states that can be determined to be uncontrollable and unobservable based on the location of 0s in sys
are removed.
Systems with numerical noise in the coefficients, e.g., noise on the order of eps
require truncation to zero to be affected by structural simplification, e.g.,
trunc_zero!(A) = A[abs.(A) .< 10eps(maximum(abs, A))] .= 0
trunc_zero!(sys.A); trunc_zero!(sys.B); trunc_zero!(sys.C)
sminreal(sys)
In contrast to minreal
, which performs pole-zero cancellation using linear-algebra operations, has an 𝑂(nₓ^3) complexity and is subject to numerical tolerances, sminreal
is computationally very cheap and numerically exact (operates on integers). However, the ability of sminreal
to reduce the order of the model is much less powerful.
See also minreal
.
ControlSystemsBase.ss
— Methodsys = ss(A, B, C, D) # Continuous
sys = ss(A, B, C, D, Ts) # Discrete
Create a state-space model sys::StateSpace{TE, T}
with matrix element type T
and TE is Continuous
or <:Discrete
.
This is a continuous-time model if Ts
is omitted. Otherwise, this is a discrete-time model with sampling period Ts
.
D
may be specified as 0
in which case a zero matrix of appropriate size is constructed automatically. sys = ss(D [, Ts])
specifies a static gain matrix D
.
To associate names with states, inputs and outputs, see named_ss
.
ControlSystemsBase.ssdata
— MethodA, B, C, D = ssdata(sys)
A destructor that outputs the statespace matrices.
ControlSystemsBase.ssrand
— Methodssrand(T::Type, ny::Int, nu::Int, nstates::Int; proper=false, stable=true, Ts=nothing)
Returns a random StateSpace
model with ny
outputs, nu
inputs, and nstates
states, whose matrix elements are normally distributed.
It is possible to specify if the system should be proper
or stable
.
Specify a sample time Ts
to obtain a discrete-time system.
ControlSystemsBase.stabregionPID
— Functionkp, ki, fig = stabregionPID(P, [ω]; kd=0, doplot=false, form=:standard)
Segments of the curve generated by this program is the boundary of the stability region for a process with transfer function P(s) The provided derivative gain is expected on parallel form, i.e., the form kp + ki/s + kd s, but the result can be transformed to any form given by the form
keyword. The curve is found by analyzing
\[P(s)C(s) = -1 ⟹ \\ |PC| = |P| |C| = 1 \\ arg(P) + arg(C) = -π\]
If P
is a function (e.g. s -> exp(-sqrt(s)) ), the stability of feedback loops using PI-controllers can be analyzed for processes with models with arbitrary analytic functions See also loopshapingPI
, loopshapingPID
, pidplots
ControlSystemsBase.starprod
— Methodstarprod(sys1, sys2, dimu, dimy)
Compute the Redheffer star product.
length(U1) = length(Y2) = dimu
and length(Y1) = length(U2) = dimy
For details, see Chapter 9.3 in Zhou, K. and JC Doyle. Essentials of robust control, Prentice hall (NJ), 1998
ControlSystemsBase.state_names
— Methodstate_names(P)
state_names(P, i)
Get a vector of strings with the names of the states of P
, or the i
:th name if an index is given.
ControlSystemsBase.stepinfo
— Methodstepinfo(res::SimResult; y0 = nothing, yf = nothing, settling_th = 0.02, risetime_th = (0.1, 0.9))
Compute the step response characteristics for a simulation result. The following information is computed and stored in a StepInfo
struct:
y0
: The initial value of the responseyf
: The final value of the responsestepsize
: The size of the steppeak
: The peak value of the responsepeaktime
: The time at which the peak occursovershoot
: The percentage overshoot of the responseundershoot
: The percentage undershoot of the response. If the step response never reaches below the initial value, the undershoot is zero.settlingtime
: The time at which the response settles withinsettling_th
of the final valuesettlingtimeind
: The index at which the response settles withinsettling_th
of the final valuerisetime
: The time at which the response rises fromrisetime_th[1]
torisetime_th[2]
of the final value
Arguments:
res
: The result from a simulation usingstep
(orlsim
)y0
: The initial value, if not provided, the first value of the response is used.yf
: The final value, if not provided, the last value of the response is used. The simulation must have reached steady-state for an automatically computed value to make sense. If the simulation has not reached steady state, you may provide the final value manually.settling_th
: The threshold for computing the settling time. The settling time is the time at which the response settles withinsettling_th
of the final value.risetime_th
: The lower and upper threshold for computing the rise time. The rise time is the time at which the response rises fromrisetime_th[1]
torisetime_th[2]
of the final value.
Example:
G = tf([1], [1, 1, 1])
res = step(G, 15)
si = stepinfo(res)
plot(si)
ControlSystemsBase.system_name
— Methodsystem_name(nothing::LTISystem)
Return the name of the system. If the system does not have a name, an empty string is returned.
ControlSystemsBase.tf
— Methodsys = tf(num, den[, Ts])
sys = tf(gain[, Ts])
Create as a fraction of polynomials:
sys::TransferFunction{SisoRational{T,TR}} = numerator/denominator
where T is the type of the coefficients in the polynomial.
num
: the coefficients of the numerator polynomial. Either scalar or vector to create SISO systems
or an array of vectors to create MIMO system.
den
: the coefficients of the denominator polynomial. Either vector to create SISO systems
or an array of vectors to create MIMO system.
Ts
: Sample time if discrete time system.
The polynomial coefficients are ordered starting from the highest order term.
Other uses:
tf(sys)
: Convertsys
totf
form.tf("s")
,tf("z")
: Create the continuous-time transfer functions
, or the discrete-time transfer functionz
.numpoly(sys)
,denpoly(sys)
: Get the numerator and denominator polynomials ofsys
as a matrix of vectors, where the outer matrix is of sizen_output × n_inputs
.
ControlSystemsBase.thiran
— Methodthiran(τ::Real, Ts)
Discretize a potentially fractional delay $τ$ as a Thiran all-pass filter with sample time Ts
.
The Thiran all-pass filter gives an a maximally flat group delay.
If $τ$ is an integer multiple of $Ts$, the Thiran all-pass filter reduces to $z^{-τ/Ts}$.
Ref: T. I. Laakso, V. Valimaki, M. Karjalainen and U. K. Laine, "Splitting the unit delay [FIR/all pass filters design]," in IEEE Signal Processing Magazine, vol. 13, no. 1, 1996.
ControlSystemsBase.time_scale
— Methodtime_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)
ControlSystemsBase.to_sized
— Methodto_sized(sys::AbstractStateSpace)
Return a HeteroStateSpace
where the system matrices are of type SizedMatrix.
NOTE: This function is fundamentally type unstable. For maximum performance, create the sized system manually, or make use of the function-barrier technique.
ControlSystemsBase.tzeros
— Methodtzeros(sys)
Compute the invariant zeros of the system sys
. If sys
is a minimal realization, these are also the transmission zeros.
ControlSystemsBase.zpconv
— Methodzpc(a,r,b,s)
form conv(a,r) + conv(b,s)
where the lengths of the polynomials are equalized by zero-padding such that the addition can be carried out
ControlSystemsBase.zpk
— Methodzpk(gain[, Ts])
zpk(num, den, k[, Ts])
zpk(sys)
Create transfer function on zero pole gain form. The numerator and denominator are represented by their poles and zeros.
sys::TransferFunction{SisoZpk{T,TR}} = k*numerator/denominator
where T
is the type of k
and TR
the type of the zeros/poles, usually Float64 and Complex{Float64}.
num
: the roots of the numerator polynomial. Either scalar or vector to create SISO systems
or an array of vectors to create MIMO system.
den
: the roots of the denominator polynomial. Either vector to create SISO systems
or an array of vectors to create MIMO system.
k
: The gain of the system. Obs, this is not the same asdcgain
.Ts
: Sample time if discrete time system.
Other uses:
zpk(sys)
: Convertsys
tozpk
form.zpk("s")
: Create the transferfunctions
.
ControlSystemsBase.zpkdata
— Methodz, 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)
LinearAlgebra.lyap
— Methodlyap(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
— Functionnorm(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.
RobustAndOptimalControl.Disk
— TypeDisk
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 pahse 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 interceptc
: Centerr
: Radiusϕm
: Angle of tangent line through origin.
If γmax < γmin the disk is inverted. See diskmargin
for disk margin computations.
RobustAndOptimalControl.Diskmargin
— TypeDiskmargin
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.
RobustAndOptimalControl.ExtendedStateSpace
— TypeExtendedStateSpace{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 outputsw
denotes external inputs, such as disturbances or referencesu
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 forExtendedStateSpace
.
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.
RobustAndOptimalControl.ExtendedStateSpace
— Methodse = 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
.
RobustAndOptimalControl.LQGProblem
— TypeG = 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
RobustAndOptimalControl.LQGProblem
— MethodLQGProblem(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
RobustAndOptimalControl.NamedStateSpace
— TypeSee named_ss
for a convenient constructor.
RobustAndOptimalControl.UncertainSS
— TypeUncertainSS{TE} <: AbstractStateSpace{TE}
Represents LFT_u(M, Diagonal(Δ))
RobustAndOptimalControl.nyquistcircles
— Typenyquistcircles(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
RobustAndOptimalControl.δ
— Typeδ(N=32)
Create an uncertain element of N
uniformly distributed samples ∈ [-1, 1]
ControlSystemsBase.G_CS
— MethodG_CS(l::LQGProblem)
ControlSystemsBase.G_PS
— MethodG_PS(l::LQGProblem)
ControlSystemsBase.input_comp_sensitivity
— Methodinput_comp_sensitivity(l::LQGProblem)
ControlSystemsBase.input_sensitivity
— Methodinput_sensitivity(l::LQGProblem)
ControlSystemsBase.output_comp_sensitivity
— Methodoutput_comp_sensitivity(l::LQGProblem)
ControlSystemsBase.output_sensitivity
— Methodoutput_sensitivity(l::LQGProblem)
ControlSystemsBase.ss
— Functionss(A, B1, B2, C1, C2, D11, D12, D21, D22 [, Ts])
Create an ExtendedStateSpace
.
DescriptorSystems.dss
— MethodDescriptorSystems.dss(sys::AbstractStateSpace)
Convert sys
to a descriptor statespace system from DescriptorSystems.jl
RobustAndOptimalControl.add_disturbance
— Methodadd_disturbance(sys::StateSpace, Ad::Matrix, Cd::Matrix)
See CCS pp. 144
Arguments:
sys
: System to augmentAd
: The dynamics of the disturbanceCd
: How the disturbance states affect the states ofsys
. This matrix has the shape (sys.nx, size(Ad, 1))
See also add_low_frequency_disturbance
, add_resonant_disturbance
RobustAndOptimalControl.add_input_differentiator
— Functionadd_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.
RobustAndOptimalControl.add_input_integrator
— Functionadd_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
RobustAndOptimalControl.add_low_frequency_disturbance
— Functionadd_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 ofsys
, and defaults tosys.B
. Ifmeasurement=true
, this matrix indicates how the disturbance states affect the outputs ofsys
, and defaults toI(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.
RobustAndOptimalControl.add_low_frequency_disturbance
— Methodadd_low_frequency_disturbance(sys::StateSpace, Ai::Integer; ϵ = 0)
A disturbance affecting only state Ai
.
RobustAndOptimalControl.add_measurement_disturbance
— Methodadd_measurement_disturbance(sys::StateSpace{Continuous}, Ad::Matrix, Cd::Matrix)
Create the system
Ae = [A 0; 0 Ad]
Be = [B; 0]
Ce = [C Cd]
RobustAndOptimalControl.add_output_differentiator
— Functionadd_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)
RobustAndOptimalControl.add_output_integrator
— Functionadd_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.
RobustAndOptimalControl.add_resonant_disturbance
— Methodadd_resonant_disturbance(sys::AbstractStateSpace, ω, ζ, Bd::AbstractArray)
Bd
: The disturbance input matrix.
RobustAndOptimalControl.add_resonant_disturbance
— Methodadd_resonant_disturbance(sys::StateSpace{Continuous}, ω, ζ, Ai::Int; measurement = false)
Augment sys
with a resonant disturbance model.
Arguments:
ω
: Frequencyζ
: Relative damping.Ai
: The affected statemeasurement
: 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 + ω²).
RobustAndOptimalControl.baltrunc2
— Methodsysr, 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.
RobustAndOptimalControl.baltrunc_coprime
— Methodsysr, 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 passingRobustAndOptimalControl.DescriptorSystems.glcf
.kwargs
: Are passed toDescriptorSystems.gbalmr
RobustAndOptimalControl.baltrunc_unstab
— Functionbaltrunc_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.
RobustAndOptimalControl.bilinearc2d
— Methodbilinearc2d(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[]
RobustAndOptimalControl.bilinearc2d
— Methodbilinearc2d(sys::ExtendedStateSpace, Ts::Number)
Applies a Balanced Bilinear transformation to a discrete-time extended statespace object
RobustAndOptimalControl.bilinearc2d
— Methodbilinearc2d(sys::StateSpace, Ts::Number)
Applies a Balanced Bilinear transformation to a discrete-time statespace object
RobustAndOptimalControl.bilineard2c
— Methodbilineard2c(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[]
RobustAndOptimalControl.bilineard2c
— Methodbilineard2c(sys::ExtendedStateSpace)
Applies a Balanced Bilinear transformation to continuous-time extended statespace object
RobustAndOptimalControl.bilineard2c
— Methodbilineard2c(sys::StateSpace)
Applies a Balanced Bilinear transformation to continuous-time statespace object
RobustAndOptimalControl.blocksort
— Methodblocks, 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 sizeN
[N, 0]
denotes a repeated complex block of sizeN
[ny, nu]
denotes a full complex block of sizeny × nu
RobustAndOptimalControl.broken_feedback
— Methodbroken_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
RobustAndOptimalControl.closedloop
— Functionclosedloop(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 reference only, other disturbance signals (B1) are ignored. See feedback
for a more advanced option.
Use static_gain_compensation
to adjust the gain from references acting on the input B2, dcgain(closedloop(l))*static_gain_compensation(l) ≈ I
RobustAndOptimalControl.connect
— Methodconnect(systems, connections; w1, z1 = (:), 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 connectedconnections
: a vector of pairs output => input, where each pair maps an output to an input. Each output must appear as an output in one ofsystems
, and similarly each input must appear as an input in one ofsystems
. 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 theaddP
block toP
's input:uP
w1
: external signals to be used as inputs in the constructed system. Use(:)
to indicate all signalsz1
: outputs of the constructed system. Use(:)
to indicate all signalsverbose
: Issue warnings for signals that have no connectionunique
: Iftrue
, all input names must be unique. Iffalse
, 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
]
w1 = [:uF] # External inputs
G = connect([F, R, C, P, addP, addC], connections; w1)
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.
RobustAndOptimalControl.controller_reduction
— Functioncontroller_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
.
RobustAndOptimalControl.controller_reduction_plot
— Functioncontroller_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.
RobustAndOptimalControl.controller_reduction_weight
— Methodcontroller_reduction_weight(P::ExtendedStateSpace, K)
Lemma 19.1 See Robust and Optimal Control Ch 19.1
RobustAndOptimalControl.dare3
— Methoddare3(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
.
RobustAndOptimalControl.diskmargin
— Functiondiskmargin(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 thehinfnorm
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
.
RobustAndOptimalControl.diskmargin
— Methoddiskmargin(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
.
RobustAndOptimalControl.diskmargin
— Methoddiskmargin(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
.
RobustAndOptimalControl.expand_symbol
— Methodexpand_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.
RobustAndOptimalControl.extended_controller
— Functionextended_controller(l::LQGProblem, L = lqr(l), K = kalman(l))
Returns an expression for 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 e = ry - y -> u
is given by
Ce = extended_controller(l)
Ce = named_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 = [:ry^l.ny, :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
RobustAndOptimalControl.extended_controller
— Methodextended_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).
RobustAndOptimalControl.feedback_control
— MethodG = 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
.
RobustAndOptimalControl.ff_controller
— Functionff_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$
See also observer_controller
.
RobustAndOptimalControl.find_lft
— Methodl, 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]
RobustAndOptimalControl.fit_complex_perturbations
— Methodcenters, 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).
RobustAndOptimalControl.frequency_separation
— Methodfrequency_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 ω
.
RobustAndOptimalControl.frequency_weighted_reduction
— Functionsysr, 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.
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)
RobustAndOptimalControl.fudge_inv
— Functionfudge_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.
RobustAndOptimalControl.gain_and_delay_uncertainty
— Methodgain_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)
RobustAndOptimalControl.gainphaseplot
— Functiongainphaseplot(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.
RobustAndOptimalControl.glover_mcfarlane
— FunctionK, γ, 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:
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.
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$
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$.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.
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).Analyze the design and if all the specifications are not met make further modifications to the weights.
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
RobustAndOptimalControl.glover_mcfarlane
— FunctionK, γ, 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"
RobustAndOptimalControl.glover_mcfarlane_2dof
— FunctionK, γ, 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 modelTref
: 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 tohinfsynthesize
.
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"])
RobustAndOptimalControl.h2norm
— Methodn = 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 pencilA-λ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
.
RobustAndOptimalControl.h2synthesize
— FunctionK, 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.
RobustAndOptimalControl.hankelnorm
— Methodn, 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
.
RobustAndOptimalControl.hanus
— MethodWh = 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"
RobustAndOptimalControl.hess_form
— Methodsysm, 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
RobustAndOptimalControl.hinfassumptions
— Methodflag = hinfassumptions(P::ExtendedStateSpace; verbose=true)
Check the assumptions for using the γ-iteration synthesis in Theorem 1.
RobustAndOptimalControl.hinfgrad
— Method∇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 ω
).
RobustAndOptimalControl.hinfnorm2
— Methodn, ω = 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
.
RobustAndOptimalControl.hinfpartition
— MethodP = 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 LTISystem
s.
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
RobustAndOptimalControl.hinfsignals
— Methodhinfsignals(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 functionsS : w → e
: From input to errorCS : w → u
: From input to controlT : w → y
: From input to output
RobustAndOptimalControl.hinfsynthesize
— MethodK, γ, 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 achievescheck
: 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 likeftype = BigFloat
.
See the example folder for example usage.
RobustAndOptimalControl.hsvd
— Methodhsvd(sys::AbstractStateSpace)
Return the Hankel singular values of sys
, computed as the eigenvalues of QP
Where Q
and P
are the Gramians of sys
.
RobustAndOptimalControl.ispassive
— Methodispassive(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
.
RobustAndOptimalControl.loop_diskmargin
— Methodloop_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
RobustAndOptimalControl.loop_diskmargin
— Methodloop_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
RobustAndOptimalControl.loop_scale
— Functionloop_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
.
RobustAndOptimalControl.loop_scaling
— Functionloop_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.
RobustAndOptimalControl.lqr3
— Methodlqr3(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)\]
RobustAndOptimalControl.makeweight
— Methodmakeweight(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 gainmid
: A number specifying the frequency at which the gain is 1, or a tuple(freq, gain)
. Ifgain_mid
is not specified, the geometric mean ofhigh
andlow
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)
RobustAndOptimalControl.measure
— Methodmeasure(s::NamedStateSpace, names)
Return a system with specified states as measurement outputs.
RobustAndOptimalControl.modal_form
— Methodsysm, 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
RobustAndOptimalControl.muplot
— Functionmuplot(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.
RobustAndOptimalControl.mvnyquistplot
— Functionfig = 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
RobustAndOptimalControl.named_ss
— Methodnamed_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.
RobustAndOptimalControl.named_ss
— Methodnamed_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) #
RobustAndOptimalControl.named_ss
— Methodnamed_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
RobustAndOptimalControl.ncfmargin
— Methodm, ω = 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
RobustAndOptimalControl.neglected_delay
— Methodneglected_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)
RobustAndOptimalControl.neglected_lag
— Methodneglected_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)
RobustAndOptimalControl.noise_mapping
— Functionnoise_mapping(P::ExtendedStateSpace)
Return the system from w -> y See also performance_mapping
, system_mapping
, noise_mapping
RobustAndOptimalControl.nu_reduction
— Functionnu_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.
RobustAndOptimalControl.nu_reduction_recursive
— Functionnu_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 aStateSpace{TE, Particles}
(a multi-model).g
: The radius of the balls in the νgap cover.gap
: An optional precomputed gapverbose
: Print progress
RobustAndOptimalControl.nugap
— Functionnugap(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
.
RobustAndOptimalControl.nugap
— Methodnugap(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.
RobustAndOptimalControl.partition
— Methodpartition(P::AbstractStateSpace, nw::Int, nz::Int)
RobustAndOptimalControl.partition
— Methodpartition(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
is the complement ofu
.z
is the complement ofy
.
RobustAndOptimalControl.passivity_index
— Methodpassivity_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
.
RobustAndOptimalControl.passivityplot
— Functionpassivityplot(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
.
RobustAndOptimalControl.performance_mapping
— Functionperformance_mapping(P::ExtendedStateSpace)
Return the system from w -> z See also performance_mapping
, system_mapping
, noise_mapping
RobustAndOptimalControl.robstab
— Methodrobstab(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.
RobustAndOptimalControl.schur_form
— Methodsysm, 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
RobustAndOptimalControl.show_construction
— Methodshow_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
RobustAndOptimalControl.sim_diskmargin
— Functionsim_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
.
RobustAndOptimalControl.sim_diskmargin
— Functionsim_diskmargin(L, σ::Real = 0, l=1e-3, u=1e3)
Return the smallest simultaneous diskmargin over the grid l:u See also ncfmargin
.
RobustAndOptimalControl.sim_diskmargin
— Methodsim_diskmargin(L, σ::Real, w::AbstractVector)
sim_diskmargin(L, σ::Real = 0)
Simultaneuous diskmargin at the outputs of L
. Users should consider using diskmargin
.
RobustAndOptimalControl.specificationplot
— Functionspecificationplot([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 pointshz
: truensigma
: typemax(Int) number of singular values to shows_labels
: `[ "σ(S)", "σ(CS)", "σ(T)",
]`
w_labels
: `[ "γ σ(Wₛ⁻¹)", "γ σ(Wᵤ⁻¹)", "γ σ(Wₜ⁻¹)",
]`
colors
:[:blue, :red, :green]
colors for $S$, $CS$ and $T$
RobustAndOptimalControl.splitter
— Functionsplitter(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 splitn
: Number of splits
RobustAndOptimalControl.ss2particles
— Methodss2particles(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
.
RobustAndOptimalControl.ssdata_e
— MethodA, B1, B2, C1, C2, D11, D12, D21, D22 = ssdata_e(sys)
RobustAndOptimalControl.stab_unstab
— Methodstab, 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).
RobustAndOptimalControl.static_gain_compensation
— Functionstatic_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
RobustAndOptimalControl.structured_singular_value
— Methodstructured_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.
RobustAndOptimalControl.structured_singular_value
— Methodμ = 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)
.
RobustAndOptimalControl.sumblock
— Methodsumblock(ex::String; Ts = 0, n = 1)
Create a summation node (named statespace system) that sums (or subtracts) vectors of length n
.
Arguments:
Ts
: Sample timen
: The length of the input and output vectors. Setn=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
RobustAndOptimalControl.sys_from_particles
— Methodsys_from_particles(P, i)
sys_from_particles(P)
Return the i
th 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
.
RobustAndOptimalControl.system_mapping
— Functionsystem_mapping(P::ExtendedStateSpace)
Return the system from u -> y See also performance_mapping
, system_mapping
, noise_mapping
RobustAndOptimalControl.uss
— Functionuss(d::AbstractVector{<:δ}, Ts = nothing)
Create a diagonal uncertain statespace object with the uncertain elements d
on the diagonal.
RobustAndOptimalControl.uss
— Functionuss(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.
RobustAndOptimalControl.uss
— Functionuss(D11, D12, D21, D22, Δ, Ts = nothing)
Create an uncertain statespace object with only gin matrices.
RobustAndOptimalControl.uss
— Methoduss(d::δ{C, F}, Ts = nothing)
Convert a δ object to an UncertainSS
RobustAndOptimalControl.vec2sys
— Functionvec2sys(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
RobustAndOptimalControl.δc
— Functionδ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.
RobustAndOptimalControl.δr
— Functionδ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.
LowLevelParticleFilters.AdvancedParticleFilter
— MethodAdvancedParticleFilter(N::Integer, dynamics::Function, measurement::Function, measurement_likelihood, dynamics_density, initial_density; p = SciMLBase.NullParameters(), threads = false, kwargs...)
This type represents a standard particle filter but affords extra flexibility compared to the ParticleFilter
type, e.g., non-additive noise in the dynamics and measurement functions.
See the docs for more information: https://baggepinnen.github.io/LowLevelParticleFilters.jl/stable/#AdvancedParticleFilter-1
Arguments:
N
: Number of particlesdynamics
: A discrete-time dynamics function(x, u, p, t, noise=false) -> x⁺
. It's important that thenoise
argument defaults tofalse
.measurement
: A measurement function(x, u, p, t, noise=false) -> y
. It's important that thenoise
argument defaults tofalse
.measurement_likelihood
: A function(x, u, y, p, t)->logl
to evaluate the log-likelihood of a measurement.dynamics_density
: This field is not used by the advanced filter and can be set tonothing
.initial_density
: The distribution of the initial state.threads
: use threads to propagate particles in parallel. Only activate this if your dynamics is thread-safe.SeeToDee.SimpleColloc
is not thread-safe by default due to the use of internal caches, butSeeToDee.Rk4
is.
LowLevelParticleFilters.AuxiliaryParticleFilter
— MethodAuxiliaryParticleFilter(args...; kwargs...)
Takes exactly the same arguments as ParticleFilter
, or an instance of ParticleFilter
.
LowLevelParticleFilters.DAEUnscentedKalmanFilter
— MethodDAEUnscentedKalmanFilter(ukf; g, get_x_z, build_xz, xz0, threads=false)
An Unscented Kalman filter for differential-algebraic systems (DAE).
Ref: "Nonlinear State Estimation of Differential Algebraic Systems", Mandela, Rengaswamy, Narasimhan
This filter is still considered experimental and subject to change without respecting semantic versioning. Use at your own risk.
Arguments
ukf
is a regularUnscentedKalmanFilter
that containsdynamics(xz, u, p, t)
that propagates the combined statexz(k)
toxz(k+1)
and a measurement function with signature(xz, u, p, t)
g(x, z, u, p, t)
is a function that should fulfillg(x, z, u, p, t) = 0
get_x_z(xz) -> x, z
is a function that decomposesxz
intox
andz
build_xz(x, z)
is the inverse ofget_x_z
xz0
the initial full state.threads
: If true, evaluates dynamics on sigma points in parallel. This typically requires the dynamics to be non-allocating (use StaticArrays) to improve performance.
Assumptions
- The DAE dynamics is index 1 and can be written on the form
\[\begin{aligned} ẋ &= f(x, z, u, p, t) \quad &\text{Differential equations}\ 0 &= g(x, z, u, p, t) \quad &\text{Algebraic equations}\ y &= h(x, z, u, p, t) \quad &\text{Measurements} \begin{aligned}\]
the measurements may be functions of both differential state variables x
and algebraic variables z
. Please note, the actual dynamics and measurement functions stored in the internal ukf
should have signatures (xz, u, p, t)
, i.e., they take the combined state (descriptor) containing both x
and z
in a single vector as dictated by the function build_xz
. It is only the function g
that is assumed to actually have the signature g(x,z,u,p,t)
.
LowLevelParticleFilters.ExtendedKalmanFilter
— TypeExtendedKalmanFilter(kf, dynamics, measurement)
ExtendedKalmanFilter(dynamics, measurement, R1,R2,d0=MvNormal(Matrix(R1)); nu::Int, p = SciMLBase.NullParameters(), α = 1.0, check = true)
A nonlinear state estimator propagating uncertainty using linearization.
The constructor to the extended Kalman filter takes dynamics and measurement functions, and either covariance matrices, or a KalmanFilter
. If the former constructor is used, the number of inputs to the system dynamics, nu
, must be explicitly provided with a keyword argument.
The filter will internally linearize the dynamics using ForwardDiff.
The dynamics and measurement function are on the following form
x(t+1) = dynamics(x, u, p, t) + w
y = measurement(x, u, p, t) + e
where w ~ N(0, R1)
, e ~ N(0, R2)
and x(0) ~ d0
See also UnscentedKalmanFilter
which is typically more accurate than ExtendedKalmanFilter
. See KalmanFilter
for detailed instructions on how to set up a Kalman filter kf
.
LowLevelParticleFilters.KalmanFilter
— TypeKalmanFilter(A,B,C,D,R1,R2,d0=MvNormal(R1); p = SciMLBase.NullParameters(), α=1)
The matrices A,B,C,D
define the dynamics
x' = Ax + Bu + w
y = Cx + Du + e
where w ~ N(0, R1)
, e ~ N(0, R2)
and x(0) ~ d0
The matrices can be time varying such that, e.g., A[:, :, t]
contains the $A$ matrix at time index t
. They can also be given as functions on the form
Afun(x, u, p, t) -> A
For maximum performance, provide statically sized matrices from StaticArrays.jl
α is an optional "forgetting factor", if this is set to a value > 1, such as 1.01-1.2, the filter will, in addition to the covariance inflation due to $R_1$, exhibit "exponential forgetting" similar to a Recursive Least-Squares (RLS) estimator. It is thus possible to get a RLS-like algorithm by setting $R_1=0, R_2 = 1/α$ and $α > 1$ ($α$ is the inverse of the traditional RLS parameter $α = 1/λ$). The exact form of the covariance update is
\[R(t+1|t) = α AR(t)A^T + R_1\]
Tutorials on Kalman filtering
The tutorial "How to tune a Kalman filter" details how to figure out appropriate covariance matrices for the Kalman filter, as well as how to add disturbance models to the system model. See also the tutorial in the documentation
LowLevelParticleFilters.ParticleFilter
— MethodParticleFilter(N::Integer, dynamics, measurement, dynamics_density, measurement_density, initial_density; threads = false, p = SciMLBase.NullParameters(), kwargs...)
See the docs for more information: https://baggepinnen.github.io/LowLevelParticleFilters.jl/stable/#Particle-filter-1
Arguments:
N
: Number of particlesdynamics
: A discrete-time dynamics function(x, u, p, t) -> x⁺
measurement
: A measurement function(x, u, p, t) -> y
dynamics_density
: A probability-density function for additive noise in the dynamics. UseAdvancedParticleFilter
for non-additive noise.measurement_density
: A probability-density function for additive measurement noise. UseAdvancedParticleFilter
for non-additive noise.initial_density
: Distribution of the initial state.
LowLevelParticleFilters.SqKalmanFilter
— TypeSqKalmanFilter(A,B,C,D,R1,R2,d0=MvNormal(R1); p = SciMLBase.NullParameters(), α=1)
A standard Kalman filter on square-root form. This filter may have better numerical performance when the covariance matrices are ill-conditioned.
The matrices A,B,C,D
define the dynamics
x' = Ax + Bu + w
y = Cx + Du + e
where w ~ N(0, R1)
, e ~ N(0, R2)
and x(0) ~ d0
The matrices can be time varying such that, e.g., A[:, :, t]
contains the $A$ matrix at time index t
. They can also be given as functions on the form
Afun(x, u, p, t) -> A
The internal fields storing covariance matrices are for this filter storing the upper-triangular Cholesky factor.
α is an optional "forgetting factor", if this is set to a value > 1, such as 1.01-1.2, the filter will, in addition to the covariance inflation due to $R_1$, exhibit "exponential forgetting" similar to a Recursive Least-Squares (RLS) estimator. It is thus possible to get a RLS-like algorithm by setting $R_1=0, R_2 = 1/α$ and $α > 1$ ($α$ is the inverse of the traditional RLS parameter $α = 1/λ$). The form of the covariance update is
\[R(t+1|t) = α AR(t)A^T + R_1\]
Ref: "A Square-Root Kalman Filter Using Only QR Decompositions", Kevin Tracy https://arxiv.org/abs/2208.06452
LowLevelParticleFilters.TupleProduct
— TypeTupleProduct(v::NTuple{N,UnivariateDistribution})
Create a product distribution where the individual distributions are stored in a tuple. Supports mixed/hybrid Continuous and Discrete distributions
LowLevelParticleFilters.UnscentedKalmanFilter
— TypeUnscentedKalmanFilter(dynamics, measurement, R1, R2, d0=MvNormal(Matrix(R1)); p = SciMLBase.NullParameters(), ny, nu)
A nonlinear state estimator propagating uncertainty using the unscented transform.
The dynamics and measurement function are on the following form
x' = dynamics(x, u, p, t) + w
y = measurement(x, u, p, t) + e
where w ~ N(0, R1)
, e ~ N(0, R2)
and x(0) ~ d0
The matrices R1, R2
can be time varying such that, e.g., R1[:, :, t]
contains the $R_1$ matrix at time index t
. They can also be given as functions on the form
Rfun(x, u, p, t) -> R
For maximum performance, provide statically sized matrices from StaticArrays.jl
ny, nu
indicate the number of outputs and inputs.
LowLevelParticleFilters.commandplot
— Functioncommandplot(pf, u, y, p=parameters(pf); kwargs...)
Produce a helpful plot. For customization options (kwargs...
), see ?pplot
. After each time step, a command from the user is requested.
- q: quit
- s n: step
n
steps
LowLevelParticleFilters.correct!
— Functioncorrect!(kf::SqKalmanFilter, u, y, p = parameters(kf), t::Real = index(kf); R2 = get_mat(kf.R2, kf.x, u, p, t))
For the square-root Kalman filter, a custom provided R2
must be the upper triangular Cholesky factor of the covariance matrix of the measurement noise.
LowLevelParticleFilters.correct!
— Functionll, e = correct!(f, u, y, p = parameters(f), t = index(f))
Update state/covariance/weights based on measurement y
, returns loglikelihood and prediction error (the error is always 0 for particle filters).
LowLevelParticleFilters.correct!
— Function(; ll, e, S, Sᵪ, K) = correct!(kf::AbstractKalmanFilter, u, y, p = parameters(kf), t::Integer = index(kf), R2)
The correct step for a Kalman filter returns not only the log likelihood ll
and the prediction error e
, but also the covariance of the output S
, its Cholesky factor Sᵪ
and the Kalman gain K
.
If R2
stored in kf
is a function R2(x, u, p, t)
, this function is evaluated at the state before the correction is performed. The measurement noise covariance matrix R2
stored in the filter object can optionally be overridden by passing the argument R2
, in this case R2
must be a matrix.
LowLevelParticleFilters.debugplot
— Functiondebugplot(pf, u, y, p=parameters(pf); runall=false, kwargs...)
Produce a helpful plot. For customization options (kwargs...
), see ?pplot
.
runall=false:
if true, runs all time steps befor displaying (faster), if false, displays the plot after each time step.
The generated plot becomes quite heavy. Initially, try limiting your input to 100 time steps to verify that it doesn't crash.
LowLevelParticleFilters.densityplot
— Functiondensityplot(x,[w])
Plot (weighted) particles densities
LowLevelParticleFilters.forward_trajectory
— Functionsol = forward_trajectory(kf::AbstractKalmanFilter, u::Vector, y::Vector, p=parameters(kf))
Run a Kalman filter forward
Returns a KalmanFilteringSolution: with the following
x
: predictionsxt
: filtered estimatesR
: predicted covariance matricesRt
: filter covariancesll
: loglik
sol
can be plotted
plot(sol::KalmanFilteringSolution; plotx = true, plotxt=true, plotu=true, ploty=true)
LowLevelParticleFilters.forward_trajectory
— Functionsol = forward_trajectory(pf, u::AbstractVector, y::AbstractVector, p=parameters(pf))
Run the particle filter for a sequence of inputs and measurements. Return a solution with x,w,we,ll = particles, weights, expweights and loglikelihood
If MonteCarloMeasurements.jl is loaded, you may transform the output particles to Matrix{MonteCarloMeasurements.Particles}
using Particles(x,we)
. Internally, the particles are then resampled such that they all have unit weight. This is conventient for making use of the plotting facilities of MonteCarloMeasurements.jl.
sol
can be plotted
plot(sol::ParticleFilteringSolution; nbinsy=30, xreal=nothing, dim=nothing)
LowLevelParticleFilters.log_likelihood_fun
— Methodll(θ) = log_likelihood_fun(filter_from_parameters(θ::Vector)::Function, priors::Vector{Distribution}, u, y, p)
returns function θ -> p(y|θ)p(θ)
LowLevelParticleFilters.loglik
— Functionll = loglik(filter, u, y, p=parameters(filter))
Calculate loglikelihood for entire sequences u,y
LowLevelParticleFilters.logsumexp!
— Functionll = logsumexp!(w, we [, maxw])
Normalizes the weight vector w
and returns the weighted log-likelihood
https://arxiv.org/pdf/1412.8695.pdf eq 3.8 for p(y) https://discourse.julialang.org/t/fast-logsumexp/22827/7?u=baggepinnen for stable logsumexp
LowLevelParticleFilters.mean_trajectory
— Methodx,ll = mean_trajectory(pf, u::Vector{Vector}, y::Vector{Vector}, p=parameters(pf))
This method resets the particle filter to the initial state distribution upon start
LowLevelParticleFilters.mean_trajectory
— Methodmean_trajectory(sol::ParticleFilteringSolution)
mean_trajectory(x::AbstractMatrix, we::AbstractMatrix)
Compute the weighted mean along the trajectory of a particle-filter solution. Returns a matrix of size T × nx
. If x
and we
are supplied, the weights are expected to be in the original space (not log space).
LowLevelParticleFilters.metropolis
— Functionmetropolis(ll::Function(θ), R::Int, θ₀::Vector, draw::Function(θ) = naive_sampler(θ₀))
Performs MCMC sampling using the marginal Metropolis (-Hastings) algorithm draw = θ -> θ'
samples a new parameter vector given an old parameter vector. The distribution must be symmetric, e.g., a Gaussian. R
is the number of iterations. See log_likelihood_fun
Example:
filter_from_parameters(θ) = ParticleFilter(N, dynamics, measurement, MvNormal(n,exp(θ[1])), MvNormal(p,exp(θ[2])), d0)
priors = [Normal(0,0.1),Normal(0,0.1)]
ll = log_likelihood_fun(filter_from_parameters,priors,u,y,1)
θ₀ = log.([1.,1.]) # Initial point
draw = θ -> θ .+ rand(MvNormal(0.1ones(2))) # Function that proposes new parameters (has to be symmetric)
burnin = 200 # If using threaded call, provide number of burnin iterations
# @time theta, lls = metropolis(ll, 2000, θ₀, draw) # Run single threaded
# thetam = reduce(hcat, theta)'
@time thetalls = LowLevelParticleFilters.metropolis_threaded(burnin, ll, 5000, θ₀, draw) # run on all threads, will provide (2000-burnin)*nthreads() samples
histogram(exp.(thetalls[:,1:2]), layout=3)
plot!(thetalls[:,3], subplot=3) # if threaded call, log likelihoods are in the last column
LowLevelParticleFilters.reset!
— Methodreset!(kf::AbstractKalmanFilter; x0)
Reset the initial distribution of the state. Optionally, a new mean vector x0
can be provided.
LowLevelParticleFilters.reset!
— MethodReset the filter to initial state and covariance/distribution
LowLevelParticleFilters.reset!
— Methodreset!(kf::SqKalmanFilter; x0)
Reset the initial distribution of the state. Optionally, a new mean vector x0
can be provided.
LowLevelParticleFilters.simulate
— Functionx,u,y = simulate(f::AbstractFilter, T::Int, du::Distribution, p=parameters(f), [N]; dynamics_noise=true, measurement_noise=true)
x,u,y = simulate(f::AbstractFilter, u, p=parameters(f); dynamics_noise=true, measurement_noise=true)
Simulate dynamical system forward in time T
steps, or for the duration of u
, returns state sequence, inputs and measurements du
is a distribution of random inputs.
A simulation can be considered a draw from the prior distribution over the evolution of the system implied by the selected noise models. Such a simulation is useful in order to evaluate whether or not the noise models are reasonable.
If MonteCarloMeasurements.jl is loaded, the argument N::Int
can be supplied, in which case N
simulations are done and the result is returned in the form of Vector{MonteCarloMeasurements.Particles}
.
LowLevelParticleFilters.smooth
— Functionxb,ll = smooth(pf, M, u, y, p=parameters(pf))
xb,ll = smooth(pf, xf, wf, wef, ll, M, u, y, p=parameters(pf))
Perform particle smoothing using forward-filtering, backward simulation. Return smoothed particles and loglikelihood. See also smoothed_trajs
, smoothed_mean
, smoothed_cov
LowLevelParticleFilters.smooth
— FunctionxT,RT,ll = smooth(kf::KalmanFilter, u::Vector, y::Vector, p=parameters(kf))
xT,RT,ll = smooth(kf::ExtendedKalmanFilter, u::Vector, y::Vector, p=parameters(kf))
Returns smoothed estimates of state x
and covariance R
given all input output data u,y
LowLevelParticleFilters.smoothed_cov
— Methodsmoothed_cov(xb)
Helper function to calculate the covariance of smoothed particle trajectories
LowLevelParticleFilters.smoothed_mean
— Methodsmoothed_mean(xb)
Helper function to calculate the mean of smoothed particle trajectories
LowLevelParticleFilters.smoothed_trajs
— Methodsmoothed_trajs(xb)
Helper function to get particle trajectories as a 3-dimensions array (N,M,T) instead of matrix of vectors.
LowLevelParticleFilters.update!
— Functionll, e = update!(f::AbstractFilter, u, y, p = parameters(f), t = index(f))
Perform one step of predict!
and correct!
, returns loglikelihood and prediction error
LowLevelParticleFilters.weighted_cov
— Methodweighted_cov(x,we)
Similar to weighted_mean
, but returns covariances
LowLevelParticleFilters.weighted_mean
— Methodx̂ = weighted_mean(x,we)
Calculated weighted mean of particle trajectories. we
are expweights.
LowLevelParticleFilters.weighted_mean
— Methodx̂ = weighted_mean(pf)
x̂ = weighted_mean(s::PFstate)
StatsAPI.predict!
— Functionpredict!(f, u, p = parameters(f), t = index(f))
Move filter state forward in time using dynamics equation and input vector u
.
StatsAPI.predict!
— Functionpredict!(kf::SqKalmanFilter, u, p = parameters(kf), t::Real = index(kf); R1 = get_mat(kf.R1, kf.x, u, p, t))
For the square-root Kalman filter, a custom provided R1
must be the upper triangular Cholesky factor of the covariance matrix of the process noise.
StatsAPI.predict!
— Functionpredict!(kf::AbstractKalmanFilter, u, p = parameters(kf), t::Integer = index(kf); R1)
Perform the prediction step (updating the state estimate to $x(t+1|t)$). If R1
stored in kf
is a function R1(x, u, p, t)
, this function is evaluated at the state before the prediciton is performed. The dynamics noise covariance matrix R1
stored in kf
can optionally be overridden by passing the argument R1
, in this case R1
must be a matrix.