Exported functions and types

ModelingToolkit.System(sys::AbstractStateSpace; name::Symbol, x0 = zeros(sys.nx), x_names, u_names, y_names)

Create an System from sys::StateSpace.

Arguments:

• sys: An instance of StateSpace or NamedStateSpace.
• name: A symbol giving the system a unique name.
• x0: Initial state

The arguments below are automatically set if the system is a NamedStateSpace.

• x_names: A vector of symbols with state names.
• u_names: A vector of symbols with input names.
• y_names: A vector of symbols with output names.

GainScheduledStateSpace(systems, vt; interpolator, x = zeros((systems[1]).nx), name, u0 = zeros((systems[1]).nu), y0 = zeros((systems[1]).ny))

A linear parameter-varying (LPV) version of Blocks.StateSpace, implementing the following equations:

\[\begin{aligned} \dot{x} &= A(v) x + B(v) u \\ y &= C(v) x + D(v) u \end{aligned}\]

where v is a scalar scheduling variable.

See example usage in the gain-scheduling example.

Arguments:

• systems: A vector of ControlSystemsBase.StateSpace objects
• vt: A vector of breakpoint values for the scheduling variable v, this has the same length as systems.
• interpolator: A constructor i = interpolator(values, breakpoints) and returns an interpolator object that can be called like i(v) to get the interpolated value at v. LinearInterpolation from DataInterpolations.jl is a good choice, but a lookup table can also be used.

Connectors

• input of type RealInput connects to $u$.
• output of type RealOutput connects to $y$.
• scheduling_input of type RealInput connects to $v$.

batch_ss(sys, inputs, outputs, ops::AbstractVector{<:AbstractDict};
            t = 0.0,
            allow_input_derivatives = false,
            kwargs...)

Linearize sys in multiple operating points ops::Vector{Dict}. Returns a vector of StateSpace objects and the simplified system.

Example:

using ControlSystemsMTK, ModelingToolkit, RobustAndOptimalControl
using ModelingToolkit: getdefault
unsafe_comparisons(true)

# Create a model
@parameters t k=10 k3=2 c=1
@variables x(t)=0 [bounds = (-2, 2)]
@variables v(t)=0
@variables u(t)=0
@variables y(t)

D = Differential(t)

eqs = [D(x) ~ v
       D(v) ~ -k * x - k3 * x^3 - c * v + 10u
       y ~ x]


@named duffing = System(eqs, t)

bounds = getbounds(duffing, unknowns(duffing))
sample_within_bounds((l, u)) = (u - l) * rand() + l
# Create a vector of operating points
N = 10
ops = map(1:N) do i
    op = Dict(x => sample_within_bounds(bounds[x]) for x in keys(bounds) if isfinite(bounds[x][1]))
end

Ps, ssys = batch_ss(duffing, [u], [y], ops)
w = exp10.(LinRange(-2, 2, 200))
bodeplot(Ps, w)
P = RobustAndOptimalControl.ss2particles(Ps) # convert to a single StateSpace system with `Particles` as coefficients.
bodeplot(P, w) # Should look similar to the one above

Let's also do some tuning for the linearized models above

function batch_tune(f, Ps)
    f.(Ps)
end

Cs = batch_tune(Ps) do P
    # C, kp, ki, fig, CF = loopshapingPI(P, 6; phasemargin=45)
    C, kp, ki, kd, fig, CF = loopshapingPID(P, 6; Mt=1.3, Tf = 1/100)
    ss(CF)
end

P = RobustAndOptimalControl.ss2particles(Ps)
C = RobustAndOptimalControl.ss2particles(Cs)

nyquistplot(P * C,
            w,
            ylims = (-10, 3),
            xlims = (-5, 10),
            points = true,
            Ms_circles = [1.5, 2],
            Mt_circles = [1.5, 2])

# Fit circles that encircle the Nyquist curve for each frequency
centers, radii = fit_complex_perturbations(P * C, w; relative = false, nominal = :center)
nyquistcircles!(w, centers, radii, ylims = (-4, 1), xlims = (-3, 4))

See also trajectory_ss and fuzz.

build_quadratic_cost_matrix(linear_sys, ssys::System, costs::Vector{Pair})

For a system that has been linearized, assemble a quadratic cost matrix (for LQR or Kalman filtering) that penalizes states or outputs of simplified system ssys according to the vector of pairs costs.

The motivation for this function is that ModelingToolkit does not guarantee

• Which states are selected as states after simplification.
• The order of the states.

The second problem above, the ordering of the states, can be worked around using reorder_states, but the first problem cannot be solved by trivial reordering. This function thus accepts an array of costs for a user-selected state realization, and assembles the correct cost matrix for the state realization selected by MTK. To do this, the funciton needs the linearization (linear_sys) as well as the simplified system, both of which are outputs of linearize.

Arguments:

• linear_sys: Output of linearize, an object containing a property called C. This can be a ControlSystemsBase.StateSpace or a NamedTuple with a field C.
• ssys: Output of linearize.
• costs: A vector of pairs

build_quadratic_cost_matrix(sys::System, inputs::Vector, costs::Vector{Pair}; kwargs...)

Assemble a quadratic cost matrix (for LQR or Kalman filtering) that penalizes states or outputs of system sys according to the vector of pairs costs.

The motivation for this function is that ModelingToolkit does not guarantee

• Which states are selected as states after simplification.
• The order of the states.

The second problem above, the ordering of the states, can be worked around using reorder_states, but the first problem cannot be solved by trivial reordering. This function thus accepts an array of costs for a user-selected state realization, and assembles the correct cost matrix for the state realization selected by MTK. To do this, the funciton performs a linearization between inputs and the cost outputs. The linearization is used to determine the matrix entries belonging to states that are not part of the realization chosen by MTK.

Arguments:

• sys: The system to be linearized (not simplified).
• inputs: A vector of variables that are to be considered controlled inputs for the LQR controller.
• costs: A vector of pairs.

get_named_comp_sensitivity(sys, ap::AnalysisPoint; kwargs...)
get_named_comp_sensitivity(sys, ap_name::Symbol; kwargs...)

Call ModelingToolkitStandardLibrary.Blocks.get_comp_sensitivity while retaining signal names. Returns a NamedStateSpace object (similar to named_ss).

get_named_looptransfer(sys, ap::AnalysisPoint; kwargs...)
get_named_looptransfer(sys, ap_name::Symbol; kwargs...)

Call ModelingToolkitStandardLibrary.Blocks.get_looptransfer while retaining signal names. Returns a NamedStateSpace object (similar to named_ss).

get_named_sensitivity(sys, ap::AnalysisPoint; kwargs...)
get_named_sensitivity(sys, ap_name::Symbol; kwargs...)

Call ModelingToolkitStandardLibrary.Blocks.get_sensitivity while retaining signal names. Returns a NamedStateSpace object (similar to named_ss).

linsystems, ssys = trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), fuzzer=nothing, verbose = true, kwargs...)

Linearize sys around the trajectory sol at times t. Returns a vector of StateSpace objects and the simplified system.

Arguments:

• inputs: A vector of variables or analysis points.
• outputs: A vector of variables or analysis points.
• sol: An ODE solution object. This solution must contain the states of the simplified system, accessible through the idxs argument like sol(t, idxs=x).
• t: Time points along the solution trajectory at which to linearize. The returned array of StateSpace objects will be of the same length as t.
• fuzzer: A function that takes an operating point dictionary and returns an array of "fuzzed" operating points. This is useful for adding noise/uncertainty to the operating points along the trajectory. See ControlSystemsMTK.fuzz for such a function.
• verbose: If true, print warnings for variables that are not found in sol.
• kwargs: Are sent to the linearization functions.
• named: If true, the returned systems will be of type NamedStateSpace, otherwise they will be of type StateSpace.

RobustAndOptimalControl.named_ss(sys::ModelingToolkit.AbstractSystem, inputs, outputs; descriptor=true, simple_infeigs=true, balance=false, kwargs...)

Convert an System to a NamedStateSpace using linearization. inputs, outputs are vectors of variables determining the inputs and outputs respectively. See docstring of ModelingToolkit.linearize for more info on kwargs.

If descriptor = true (default), this method automatically converts systems that MTK has failed to produce a proper form for into a proper linear statespace system using the method described here: https://juliacontrol.github.io/ControlSystemsMTK.jl/dev/#Internals:-Transformation-of-non-proper-models-to-proper-statespace-form If descriptor = false, the system is instead converted to a statespace realization using sys[:,uinds] + sys[:,duinds]*tf('s'), which tends to result in a larger realization on which the user may want to call minreal(sys, tol) with a carefully selected tolerance.

See also ModelingToolkit.linearize which is the lower-level function called internally. The functions get_named_sensitivity, get_named_comp_sensitivity, get_named_looptransfer similarily provide convenient ways to compute sensitivity functions while retaining signal names in the same way as named_ss. The corresponding lower-level functions get_sensitivity, get_comp_sensitivity and get_looptransfer are available in ModelingToolkitStandardLibrary.Blocks and are documented in MTKstdlib: Linear analysis.

fig = 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: Call balance_statespace on the system before plotting.
• adjust_phase_start: If true, the phase will be adjusted so that it starts at -90*intexcess degrees, where intexcess is the integrator excess of the system.
• adaptive: If true, an adaptive frequency grid is used in order to keep the number of plotted points low, while resolving features in the frequency response well. If a manually provided frequency vector is used, this may be downsampled before plotting.

kwargs is sent as argument to RecipesBase.plot.

(; A, B, C, D), simplified_sys, extras = linearize(sys, inputs, outputs;    t=0.0, op = Dict(), allow_input_derivatives = false, zero_dummy_der=false, kwargs...)
(; A, B, C, D), extras                 = linearize(simplified_sys, lin_fun; t=0.0, op = Dict(), allow_input_derivatives = false, zero_dummy_der=false)

Linearize sys between inputs and outputs, both vectors of variables. Return a NamedTuple with the matrices of a linear statespace representation on the form

\[\begin{aligned} ẋ &= Ax + Bu\\ y &= Cx + Du \end{aligned}\]

The first signature automatically calls linearization_function internally, while the second signature expects the outputs of linearization_function as input.

op denotes the operating point around which to linearize. If none is provided, the default values of sys are used.

If allow_input_derivatives = false, an error will be thrown if input derivatives ($u̇$) appear as inputs in the linearized equations. If input derivatives are allowed, the returned B matrix will be of double width, corresponding to the input [u; u̇].

zero_dummy_der can be set to automatically set the operating point to zero for all dummy derivatives.

The return value extras is a NamedTuple (; x, p, t) containing the result of the initialization problem that was solved to determine the operating point.

See also linearization_function which provides a lower-level interface, linearize_symbolic and ModelingToolkit.reorder_unknowns.

See extended help for an example.

The implementation and notation follows that of "Linear Analysis Approach for Modelica Models", Allain et al. 2009

Extended help

This example builds the following feedback interconnection and linearizes it from the input of F to the output of P.

  r ┌─────┐       ┌─────┐     ┌─────┐
───►│     ├──────►│     │  u  │     │
    │  F  │       │  C  ├────►│  P  │ y
    └─────┘     ┌►│     │     │     ├─┬─►
                │ └─────┘     └─────┘ │
                │                     │
                └─────────────────────┘

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
function plant(; name)
    @variables x(t) = 1
    @variables u(t)=0 y(t)=0
    eqs = [D(x) ~ -x + u
           y ~ x]
    System(eqs, t; name = name)
end

function ref_filt(; name)
    @variables x(t)=0 y(t)=0
    @variables u(t)=0 [input = true]
    eqs = [D(x) ~ -2 * x + u
           y ~ x]
    System(eqs, t, name = name)
end

function controller(kp; name)
    @variables y(t)=0 r(t)=0 u(t)=0
    @parameters kp = kp
    eqs = [
        u ~ kp * (r - y),
    ]
    System(eqs, t; name = name)
end

@named f = ref_filt()
@named c = controller(1)
@named p = plant()

connections = [f.y ~ c.r # filtered reference to controller reference
               c.u ~ p.u # controller output to plant input
               p.y ~ c.y]

@named cl = System(connections, t, systems = [f, c, p])

lsys0, ssys = linearize(cl, [f.u], [p.x])
desired_order = [f.x, p.x]
lsys = ModelingToolkit.reorder_unknowns(lsys0, unknowns(ssys), desired_order)

@assert lsys.A == [-2 0; 1 -2]
@assert lsys.B == [1; 0;;]
@assert lsys.C == [0 1]
@assert lsys.D[] == 0

## Symbolic linearization
lsys_sym, _ = ModelingToolkit.linearize_symbolic(cl, [f.u], [p.x])

@assert substitute(lsys_sym.A, ModelingToolkit.defaults(cl)) == lsys.A

lin_fun, simplified_sys = linearization_function(sys::AbstractSystem, inputs, outputs; simplify = false, initialize = true, initialization_solver_alg = TrustRegion(), kwargs...)

Return a function that linearizes the system sys. The function linearize provides a higher-level and easier to use interface.

lin_fun is a function (variables, p, t) -> (; f_x, f_z, g_x, g_z, f_u, g_u, h_x, h_z, h_u), i.e., it returns a NamedTuple with the Jacobians of f,g,h for the nonlinear sys (technically for simplified_sys) on the form

\[\begin{aligned} ẋ &= f(x, z, u) \\ 0 &= g(x, z, u) \\ y &= h(x, z, u) \end{aligned}\]

where x are differential unknown variables, z algebraic variables, u inputs and y outputs. To obtain a linear statespace representation, see linearize. The input argument variables is a vector defining the operating point, corresponding to unknowns(simplified_sys) and p is a vector corresponding to the parameters of simplified_sys. Note: all variables in inputs have been converted to parameters in simplified_sys.

The simplified_sys has undergone mtkcompile and had any occurring input or output variables replaced with the variables provided in arguments inputs and outputs. The unknowns of this system also indicate the order of the unknowns that holds for the linearized matrices.

Arguments:

• sys: A System of ODEs. This function will automatically apply simplification passes on sys and return the resulting simplified_sys.
• inputs: A vector of variables that indicate the inputs of the linearized input-output model.
• outputs: A vector of variables that indicate the outputs of the linearized input-output model.
• simplify: Apply simplification in tearing.
• initialize: If true, a check is performed to ensure that the operating point is consistent (satisfies algebraic equations). If the op is not consistent, initialization is performed.
• initialization_solver_alg: A NonlinearSolve algorithm to use for solving for a feasible set of state and algebraic variables that satisfies the specified operating point.
• autodiff: An ADType supported by DifferentiationInterface.jl to use for calculating the necessary jacobians. Defaults to using AutoForwardDiff()
• kwargs: Are passed on to find_solvables!

See also linearize which provides a higher-level interface.

C, 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 in fig.
• lb: log10 of lower bound for kd.
• ub: log10 of upper bound for kd.
• Tf: Time constant for second-order measurement noise filter on the form tf(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ω or Tf = 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)

RobustAndOptimalControl.named_ss(sys::ModelingToolkit.AbstractSystem, inputs, outputs; descriptor=true, simple_infeigs=true, balance=false, kwargs...)

Convert an System to a NamedStateSpace using linearization. inputs, outputs are vectors of variables determining the inputs and outputs respectively. See docstring of ModelingToolkit.linearize for more info on kwargs.

If descriptor = true (default), this method automatically converts systems that MTK has failed to produce a proper form for into a proper linear statespace system using the method described here: https://juliacontrol.github.io/ControlSystemsMTK.jl/dev/#Internals:-Transformation-of-non-proper-models-to-proper-statespace-form If descriptor = false, the system is instead converted to a statespace realization using sys[:,uinds] + sys[:,duinds]*tf('s'), which tends to result in a larger realization on which the user may want to call minreal(sys, tol) with a carefully selected tolerance.

See also ModelingToolkit.linearize which is the lower-level function called internally. The functions get_named_sensitivity, get_named_comp_sensitivity, get_named_looptransfer similarily provide convenient ways to compute sensitivity functions while retaining signal names in the same way as named_ss. The corresponding lower-level functions get_sensitivity, get_comp_sensitivity and get_looptransfer are available in ModelingToolkitStandardLibrary.Blocks and are documented in MTKstdlib: Linear analysis.

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

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

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

Arguments:

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

Example

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

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

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

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

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

named_ss(sys::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

(; A, B, C, D), simplified_sys = linearize_symbolic(sys::AbstractSystem, inputs, outputs; simplify = false, allow_input_derivatives = false, kwargs...)

Similar to linearize, but returns symbolic matrices A,B,C,D rather than numeric. While linearize uses ForwardDiff to perform the linearization, this function uses Symbolics.jacobian.

See linearize for a description of the arguments.

Extended help

The named tuple returned as the first argument additionally contains the jacobians f_x, f_z, g_x, g_z, f_u, g_u, h_x, h_z, h_u of

\[\begin{aligned} ẋ &= f(x, z, u) \\ 0 &= g(x, z, u) \\ y &= h(x, z, u) \end{aligned}\]

where x are differential unknown variables, z algebraic variables, u inputs and y outputs.

    get_sensitivity(sys, ap::AnalysisPoint; kwargs)
    get_sensitivity(sys, ap_name::Symbol; kwargs)

Compute the sensitivity function in analysis point ap. The sensitivity function is obtained by introducing an infinitesimal perturbation d at the input of ap, linearizing the system and computing the transfer function between d and the output of ap.

Arguments:

• kwargs: Are sent to ModelingToolkit.linearize

See also get_comp_sensitivity, get_looptransfer.

get_comp_sensitivity(sys, ap::AnalysisPoint; kwargs)
get_comp_sensitivity(sys, ap_name::Symbol; kwargs)

Compute the complementary sensitivity function in analysis point ap. The complementary sensitivity function is obtained by introducing an infinitesimal perturbation d at the output of ap, linearizing the system and computing the transfer function between d and the input of ap.

Arguments:

• kwargs: Are sent to ModelingToolkit.linearize

See also get_sensitivity, get_looptransfer.

get_looptransfer(sys, ap::AnalysisPoint; kwargs)
get_looptransfer(sys, ap_name::Symbol; kwargs)

Compute the (linearized) loop-transfer function in analysis point ap, from ap.out to ap.in.

Arguments:

• kwargs: Are sent to ModelingToolkit.linearize

See also get_sensitivity, get_comp_sensitivity, open_loop.

StateSpace(A, B, C, D = 0; x = zeros(size(A,1)), u0 = zeros(size(B,2)), y0 = zeros(size(C,1)), name)

A linear, time-invariant state-space system on the form.

\[\begin{aligned} ẋ &= Ax + Bu \\ y &= Cx + Du \end{aligned}\]

Transfer functions can also be simulated by converting them to a StateSpace form.

y0 and u0 can be used to set an operating point, providing them changes the dynamics from an LTI system to the affine system

\[\begin{aligned} ẋ &= Ax + B(u - u0) \\ y &= Cx + D(u - u0) + y0 \end{aligned}\]

For a nonlinear system

\[\begin{aligned} ẋ &= f(x, u) \\ y &= h(x, u) \end{aligned}\]

linearized around the operating point x₀, u₀, we have y0, u0 = h(x₀, u₀), u₀.

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

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

See also sys_from_particles.

StateSpace{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

ccode(G; simplify = identity, cse = true)

Return a string with C-code for filtering a signal u through G.

If G is a transfer function, the system must be SISO, for MIMO systems, use a StateSpace model instead.

With a transfer function as input, the code will return a double corresponding to the single output. With a StateSpace model as input, the code will produce a function that takes a double pointer double *y as the first input argument. Make sure that y points to an array of length G.ny before calling the function.

The state is internally handled by C static variables, so the generated code is thus stateful.

Arguments:

• G: A linear system
• simplify: A function for symbolic simplification. You may try Sympy.simplify, but for large systems, this will take a long time to compute.
• cse: Perform common subexpression elimination. This generally improves the performance of the generated code.

print_c_array(io, a::Vector{<:AbstractArray}, t::AbstractVector, name = "mat"; cse = false, s = "", print_vector = true, print_logic = true, struct_name::Union{Nothing, String} = nothing, struct_type = nothing, ivecname = name * "_interp_vect")

Write C-code for interpolating between arrays a. The array t contains the interpolation points.

print_c_array(io, sys::Vector{<:AbstractStateSpace}, t::AbstractVector, name = "sys"; cse = false, s = "", en = "", struct_name::Union{Nothing, String} = nothing, struct_type = nothing)

Write C-code for an interpolated linear system. The interpolation vector t defines the interpolation points, this vector is expected to be of the same length as the vector of linear systems sys.

• s, en: are strings that are appended at the start and end of variables names in the C-code.
• struct_name: If provided, the interpolation matrices will be placed inside a struct with this name.
• struct_type: If the struct name is used, provide also the C type of the struct.

fuzz(op, p; N = 10, parameters = true, variables = true)

"Fuzz" an operating point op::Dict by changing each non-zero value to an uncertain number with multiplicative uncertainty p, represented by N samples, i.e., p = 0.1 means that the value is multiplied by a N numbers between 0.9 and 1.1.

parameters and variables indicate whether to fuzz parameters and state variables, respectively.

This function modifies all variables the same way. For more fine-grained control, load the MonteCarloMeasurements package and use the Particles type directly, followed by MonteCarloMeasurements.particle_dict2dict_vec(op), i.e., the following makes uncertain_var uncertain with a 10% uncertainty:

using MonteCarloMeasurements
op = ModelingToolkit.defaults(sys)
op[uncertain_var] = op[uncertain_var] * Particles(10, Uniform(0.9, 1.1))
ops = MonteCarloMeasurements.particle_dict2dict_vec(op)
batch_ss(model, inputs, outputs, ops)

If you have more than one uncertain parameter, it's important to use the same number of particles for all of them (10 in the example above).

To make use of this function in trajectory_ss, pass something like

fuzzer = op -> ControlSystemsMTK.fuzz(op, 0.02; N=10)

to fuzz each operating point 10 times with a 2% uncertainty. The resulting number of operating points will increase by 10x.