|
| 1 | +# Disturbance and input modeling modeling |
| 2 | + |
| 3 | +Disturbances are often seen as external factors that influence a system. Modeling and simulation of such external influences is common in order to ensure that the plant and or control system can adequately handle or suppress these disturbances. Disturbance modeling is also integral to the problem of state estimation, indeed, modeling how disturbances affect the evolution of the state of the system is crucial in order to accurately estimate this state. |
| 4 | + |
| 5 | +This tutorial will show how to model disturbances in ModelingToolkit as _disturbance inputs_. This involves demonstrating best practices that make it easy to use a single model to handle both disturbed and undisturbed systems, and making use of the model for both simulation and state estimation. |
| 6 | + |
| 7 | +## A flexible component-based workflow |
| 8 | + |
| 9 | +We will consider a simple system consisting of two inertias connected through a flexible shaft, such as a simple transmission system in a car. We start by modeling the plant _without any input signals_: |
| 10 | + |
| 11 | +```@example DISTURBANCE_MODELING |
| 12 | +using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, Test |
| 13 | +using ModelingToolkitStandardLibrary.Mechanical.Rotational |
| 14 | +using ModelingToolkitStandardLibrary.Blocks |
| 15 | +t = ModelingToolkit.t_nounits |
| 16 | +D = ModelingToolkit.D_nounits |
| 17 | +
|
| 18 | +@mtkmodel SystemModel begin |
| 19 | + @parameters begin |
| 20 | + m1 = 1 |
| 21 | + m2 = 1 |
| 22 | + k = 10 # Spring stiffness |
| 23 | + c = 3 # Damping coefficient |
| 24 | + end |
| 25 | + @components begin |
| 26 | + inertia1 = Inertia(; J = m1, phi = 0, w = 0) |
| 27 | + inertia2 = Inertia(; J = m2, phi = 0, w = 0) |
| 28 | + spring = Spring(; c = k) |
| 29 | + damper = Damper(; d = c) |
| 30 | + torque = Torque(use_support = false) |
| 31 | + end |
| 32 | + @equations begin |
| 33 | + connect(torque.flange, inertia1.flange_a) |
| 34 | + connect(inertia1.flange_b, spring.flange_a, damper.flange_a) |
| 35 | + connect(inertia2.flange_a, spring.flange_b, damper.flange_b) |
| 36 | + end |
| 37 | +end |
| 38 | +``` |
| 39 | + |
| 40 | +Here, we have added a `torque` component that allows us to add a torque input to drive the system, but we have not connected any signal to it yet. We have not yet made any attempts at modeling disturbances, and this is deliberate, we will handle this later in order to make the plant model as generically useful as possible. |
| 41 | + |
| 42 | +In order to simulate this system in the presence of disturbances, we must 1. Reason about how disturbances may affect the system, and 2. attach _disturbance inputs_ and _disturbance signals_ to the model. We distinguish between an _input_ and a _signal_ here, where we by _input_ mean an attachment point (connector) to which we may connect a _signal_, i.e., a time-varying function. |
| 43 | + |
| 44 | +We create a new model that includes disturbance inputs and signals, and attach those to the already defined plant model. We assume that each of the two inertias can be affected by a disturbance torque, such as due to friction or an unknown load on the output inertia. |
| 45 | + |
| 46 | +```@example DISTURBANCE_MODELING |
| 47 | +@mtkmodel ModelWithInputs begin |
| 48 | + @components begin |
| 49 | + input_signal = Blocks.Sine(frequency = 1, amplitude = 1) |
| 50 | + disturbance_signal1 = Blocks.Step(height = -1, start_time = 2) # We add an input signal that equals zero by default so that it has no effect during normal simulation |
| 51 | + disturbance_signal2 = Blocks.Step(height = 2, start_time = 4) |
| 52 | + disturbance_torque1 = Torque(use_support = false) |
| 53 | + disturbance_torque2 = Torque(use_support = false) |
| 54 | + system_model = SystemModel() |
| 55 | + end |
| 56 | + @equations begin |
| 57 | + connect(input_signal.output, :u, system_model.torque.tau) |
| 58 | + connect(disturbance_signal1.output, :d1, disturbance_torque1.tau) # When we connect the input _signals_, we do so through an analysis point. This allows us to easily disconnect the input signals in situations when we do not need them. |
| 59 | + connect(disturbance_signal2.output, :d2, disturbance_torque2.tau) |
| 60 | + connect(disturbance_torque1.flange, system_model.inertia1.flange_b) |
| 61 | + connect(disturbance_torque2.flange, system_model.inertia2.flange_b) |
| 62 | + end |
| 63 | +end |
| 64 | +``` |
| 65 | + |
| 66 | +This outer model, `ModelWithInputs`, contains two disturbance inputs, both of type `Torque`. It also contains three signal specifications, one for the control input and two for the corresponding disturbance inputs. Note how we added the disturbance torque inputs at this level of the model, but the control input was added inside the system model. This is a design choice that is up to the modeler, here, we consider the driving torque to be a fundamental part of the model that is always required to make use of it, while the disturbance inputs may be of interest only in certain situations, and we thus add them when needed. Since we have added not only input connectors, but also connected input signals to them, this model is complete and ready for simulation, i.e., there are no _unbound inputs_. |
| 67 | + |
| 68 | +```@example DISTURBANCE_MODELING |
| 69 | +@named model = ModelWithInputs() # Model with load disturbance |
| 70 | +ssys = structural_simplify(model) |
| 71 | +prob = ODEProblem(ssys, [], (0.0, 6.0)) |
| 72 | +sol = solve(prob, Tsit5()) |
| 73 | +using Plots |
| 74 | +plot(sol) |
| 75 | +``` |
| 76 | + |
| 77 | +A thing to note in the specification of `ModelWithInputs` is the presence of three _analysis points_, `:u`, `:d1`, and `:d2`. When signals are connected through an analysis point, we may at any time linearize the model as if the signals were not connected, i.e., as if the corresponding inputs were unbound. We may also use this to generate a julia function for the dynamics on the form ``f(x,u,p,t,w)`` where the input ``u`` and disturbance ``w`` may be provided as separate function arguments, as if the corresponding input signals were not present in the model. More details regarding this will be presented further below, here, we just demonstrate how we could linearize this system model from the inputs to the angular velocity of the inertias |
| 78 | + |
| 79 | +```@example DISTURBANCE_MODELING |
| 80 | +using ControlSystemsBase, ControlSystemsMTK # ControlSystemsMTK provides the high-level function named_ss and ControlSystemsBase provides the bodeplot function |
| 81 | +P = named_ss(model, [ssys.u, ssys.d1, ssys.d2], |
| 82 | + [ssys.system_model.inertia1.w, ssys.system_model.inertia2.w]) |
| 83 | +bodeplot(P, plotphase = false) |
| 84 | +``` |
| 85 | + |
| 86 | +It's worth noting at this point that the fact that we could connect disturbance outputs from outside of the plant-model definition was enabled by the fact that we used a component-based workflow, where the plant model had the appropriate connectors available. If the plant model had modeled the system using direct equations without connectors, this would not have been possible and the model would thus be significantly less flexible. |
| 87 | + |
| 88 | +We summarize the findings so far as a number of best practices: |
| 89 | + |
| 90 | +!!! tip "Best practices" |
| 91 | + |
| 92 | + - Use a component-based workflow to model the plant |
| 93 | + - Model the plant without disturbance inputs to make it as generic as possible |
| 94 | + - When disturbance inputs are needed, create a new model that includes the plant model and the disturbance inputs |
| 95 | + - Only add input _signals_ at the top level of the model, this applies to both control inputs and disturbance inputs. |
| 96 | + - Use analysis points to connect signals to inputs, this allows for easy disconnection of signals when needed, e.g., for linearization or function generation. |
| 97 | + |
| 98 | +## Modeling for state estimation |
| 99 | + |
| 100 | +In the example above, we constructed a model for _simulation_ of a disturbance affecting the system. When simulating, we connect an input signal of specified shape that simulates the disturbance, above, we used `Blocks.Step` as input signals. On the other hand, when performing state estimation, the exact shape of the disturbance is typically not known, we might only have some diffuse knowledge of the disturbance characteristics such as "varies smoothly", "makes sudden step changes" or "is approximately periodic with 24hr period". The encoding of such knowledge is commonly reasoned about in the frequency domain, where we specify a disturbance model as a dynamical system with a frequency response similar to the approximate spectrum of the disturbance. [For more details around this, see the in-depth tutorial notebook "How to tune a Kalman filter"](https://juliahub.com/pluto/editor.html?id=ad9ecbf9-bf83-45e7-bbe8-d2e5194f2240). Most algorithms for state estimation, such as a Kalman-filter like estimators, assume that disturbances are independent and identically distributed (i.i.d.). While seemingly restrictive at first glance, when combined with an appropriate disturbance models encoded as dynamical systems, this assumption still allows for a wide range of non i.i.d. disturbances to be modeled. |
| 101 | + |
| 102 | +When modeling a system in MTK, we essentially (without considering algebraic equations for simplicity in exposition) construct a model of a dynamical system |
| 103 | + |
| 104 | +```math |
| 105 | +\begin{aligned} |
| 106 | +\dot x &= f(x, p, t) \\ |
| 107 | +y &= g(x, p, t) |
| 108 | +``` |
| 109 | + |
| 110 | +where ``x`` is the state, ``y`` are observed variables, ``p`` are parameters, and ``t`` is time. When using MTK, which variables constitute ``x`` and which are considered part of the output, ``y``, is up to the tool rather than the user, this choice is made by MTK during the call to `@mtkbuild` or the lower-level function `structural_simplify`. |
| 111 | + |
| 112 | +If we further consider external inputs to the system, such as controlled input signals ``u`` and disturbance inputs ``w``, we can write the system as |
| 113 | + |
| 114 | +```math |
| 115 | +\begin{aligned} |
| 116 | +\dot x &= f(x, u, p, t, w) \\ |
| 117 | +y &= g(x, u, p, t) |
| 118 | +\end{aligned} |
| 119 | +``` |
| 120 | + |
| 121 | +To make use of the model defined above for state estimation, we may want to generate a Julia functions for the dynamics ``f`` and the output equations ``g`` that we can plug into, e.g., a nonlinear version of a Kalman filter or a particle filter, etc. MTK contains utilities to do this, namely, [`generate_control_function`](@ref) and [`build_explicit_observed_function`](@ref) (described in more details in ["Input output"](@ref inputoutput)). These functions take keyword arguments `disturbance_inputs` and `disturbance_argument`, that indicate which variables in the model are considered part of ``w``, and whether or not these variables are to be added as function arguments to ``f``, i.e., whether we have ``f(x, u, p, t)`` or ``f(x, u, p, t, w)``. If we do not include the disturbance inputs as function arguments, MTK will assume that the ``w`` variables are all zero, but any dynamics associated with these variables, such as disturbance models, will be included in the generated function. This allows a state estimator to estimate the state of the disturbance model, provided that this state is [observable](https://en.wikipedia.org/wiki/Observability) from the measured outputs of the system. |
| 122 | + |
| 123 | +Below, we demonstrate |
| 124 | + |
| 125 | + 1. How to add an integrating disturbance model |
| 126 | + 2. how to generate the functions ``f`` and ``g`` for a typical nonlinear state estimator with explicit disturbance inputs |
| 127 | + |
| 128 | +```@example DISTURBANCE_MODELING |
| 129 | +@mtkmodel IntegratingDisturbance begin |
| 130 | + @variables begin |
| 131 | + x(t) = 0.0 |
| 132 | + w(t) = 0.0, [disturbance = true, input = true] |
| 133 | + end |
| 134 | + @components begin |
| 135 | + input = RealInput() |
| 136 | + output = RealOutput() |
| 137 | + end |
| 138 | + @equations begin |
| 139 | + D(x) ~ w |
| 140 | + w ~ input.u |
| 141 | + output.u ~ x |
| 142 | + end |
| 143 | +end |
| 144 | +
|
| 145 | +@mtkmodel SystemModelWithDisturbanceModel begin |
| 146 | + @components begin |
| 147 | + input_signal = Blocks.Sine(frequency = 1, amplitude = 1) |
| 148 | + disturbance_signal1 = Blocks.Constant(k = 0) |
| 149 | + disturbance_signal2 = Blocks.Constant(k = 0) |
| 150 | + disturbance_torque1 = Torque(use_support = false) |
| 151 | + disturbance_torque2 = Torque(use_support = false) |
| 152 | + disturbance_model = Blocks.Integrator() |
| 153 | + system_model = SystemModel() |
| 154 | + end |
| 155 | + @equations begin |
| 156 | + connect(input_signal.output, :u, system_model.torque.tau) |
| 157 | + connect(disturbance_signal1.output, :d1, disturbance_model.input) |
| 158 | + connect(disturbance_model.output, disturbance_torque1.tau) |
| 159 | + connect(disturbance_signal2.output, :d2, disturbance_torque2.tau) |
| 160 | + connect(disturbance_torque1.flange, system_model.inertia1.flange_b) |
| 161 | + connect(disturbance_torque2.flange, system_model.inertia2.flange_b) |
| 162 | + end |
| 163 | +end |
| 164 | +
|
| 165 | +@named model_with_disturbance = SystemModelWithDisturbanceModel() |
| 166 | +``` |
| 167 | + |
| 168 | +We demonstrate that this model is complete and can be simulated: |
| 169 | + |
| 170 | +```@example DISTURBANCE_MODELING |
| 171 | +ssys = structural_simplify(model_with_disturbance) |
| 172 | +prob = ODEProblem(ssys, [], (0.0, 10.0)) |
| 173 | +sol = solve(prob, Tsit5()) |
| 174 | +using Test |
| 175 | +@test SciMLBase.successful_retcode(sol) |
| 176 | +``` |
| 177 | + |
| 178 | +but we may also generate the functions ``f`` and ``g`` for state estimation: |
| 179 | + |
| 180 | +```@example DISTURBANCE_MODELING |
| 181 | +inputs = [ssys.u] |
| 182 | +disturbance_inputs = [ssys.d1, ssys.d2] |
| 183 | +P = ssys.system_model |
| 184 | +outputs = [P.inertia1.phi, P.inertia2.phi, P.inertia1.w, P.inertia2.w] |
| 185 | +
|
| 186 | +(f_oop, f_ip), x_sym, p_sym, io_sys = ModelingToolkit.generate_control_function( |
| 187 | + model_with_disturbance, inputs, disturbance_inputs; disturbance_argument = true) |
| 188 | +
|
| 189 | +g = ModelingToolkit.build_explicit_observed_function( |
| 190 | + io_sys, outputs; inputs) |
| 191 | +
|
| 192 | +op = ModelingToolkit.inputs(io_sys) .=> 0 |
| 193 | +x0, _ = ModelingToolkit.get_u0_p(io_sys, op, op) |
| 194 | +p = MTKParameters(io_sys, op) |
| 195 | +u = zeros(1) # Control input |
| 196 | +w = zeros(length(disturbance_inputs)) # Disturbance input |
| 197 | +@test f_oop(x0, u, p, t, w) == zeros(5) |
| 198 | +@test g(x, u, p, 0.0) == [0, 0, 0, 0] |
| 199 | +
|
| 200 | +# Non-zero disturbance inputs should result in non-zero state derivatives. We call `sort` since we do not generally know the order of the state variables |
| 201 | +w = [1.0, 2.0] |
| 202 | +@test sort(f_oop(x0, u, p, t, w)) == [0, 0, 0, 1, 2] |
| 203 | +``` |
| 204 | + |
| 205 | +## Further reading |
| 206 | + |
| 207 | +To see full examples that perform state estimation with ModelingToolkit models, see the following resources: |
| 208 | + |
| 209 | + - [C codegen considered unnecessary: go directly to binary, do not pass C. Compilation of Julia code for deployment in model-based engineering](https://arxiv.org/abs/2502.01128) |
| 210 | + - [LowLevelParticleFiltersMTK.jl](https://github.com/baggepinnen/LowLevelParticleFiltersMTK.jl) |
| 211 | + |
| 212 | +## Index |
| 213 | + |
| 214 | +```@index |
| 215 | +Pages = ["disturbance_modeling.md"] |
| 216 | +``` |
| 217 | + |
| 218 | +```@autodocs |
| 219 | +Modules = [ModelingToolkit] |
| 220 | +Pages = ["systems/analysis_points.jl"] |
| 221 | +Order = [:function, :type] |
| 222 | +Private = false |
| 223 | +``` |
0 commit comments