Repository navigation
Add linearize_hybrid for partition-wise linearization of clocked systems - #5173
baggepinnen wants to merge 8 commits into
Conversation
Benchmark Results (Julia vlts)Time benchmarks
Memory benchmarks
|
Benchmark Results (Julia v1)Time benchmarks
Memory benchmarks
|
|
I tried this branch out in three settings: its own test file, a set of small models exercising On the small models the partitioning, the boundary bookkeeping and the connection indices come out Findings below, roughly by severity. 1. A partition inherits the assertions of the whole model
The continuous partition inheriting an assertion on a clocked variable: using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
dt = 0.1
@variables x(t) y(t) u(t) yd(t) ud(t)
eqs = [yd ~ Sample(dt)(y), ud ~ -2yd, u ~ Hold(ud), D(x) ~ -x + u, y ~ x]
@named sys = System(eqs, t; assertions = Dict(ud < 10.0 => "ud out of range"))
linearize_hybrid(sys, [], [x]; op = Dict(x => 0.0, ud => 0.0, y => 0.0))
# UndefVarError: `ud` not defined in `ModelingToolkitBase`and a clocked partition inheriting an assertion on a continuous variable, with k = ShiftIndex(Clock(dt))
eqs = [yd ~ Sample(dt)(y), ud ~ 0.5ud(k - 1) - 2yd, u ~ Hold(ud), D(x) ~ -x + u, y ~ x]
@named sys = System(eqs, t; assertions = Dict(x < 10.0 => "x out of range"))
linearize_hybrid(sys, [], [x]; op = Dict(x => 0.0, ud => 0.0, y => 0.0), warn_missing_op = false)
# UndefVarError: `x` not defined in `ModelingToolkit`A clocked partition with no state variable of its own hides the problem, since the assertion is This is what the heat-pump model trips on, see section 3b. 3. The heat-pump loop
(a) (b) Giving all fifty a value in (c) Filtering the assertions down to the symbols of the partition locally, both clocked from 4. The initialization system of the continuous partition is overdeterminedLinearizing the continuous partition of a hybrid model emits It fired five times in one run of so the surplus equation appears once the operating point has fixed both 5. QAFive qualified accesses to non-public ModelingToolkitBase names, all from the new file, fail the QA 6. Smaller observations
|
|
Thank you for the thorough review. All findings are addressed in the latest push; details per item below. 1 and 3b, assertions. 3a, parameters of other partitions. Each partition system now keeps only the parameters that occur in its equations, observed equations, the initial values and guesses of its unknowns, its boundary and input parameters, and, transitively, the bindings and initial values of those. A parameter of the plant bound to 3c and 4, initialization constraints. Three changes:
On the question whether the 5, QA. The five qualified accesses are added to 6.
CI. The downstream test file is no longer registered in the test groups and its DyadRegistry-only dependencies are removed from |
|
Retried everything on 10af7f8, same environment as before (Julia 1.13.0 for the in-repo tests, What is fixed
On the heat-pump loop, run with no workarounds at all, the three failures I reported are gone: no What still failsOne failure remains, in the continuous partition only, after both clocked partitions are through. raised from Giving all fifty a value in at the same place, this time a time-dependent parameter. What I ruled out, in case it saves someone the same detour:
I also have no first-hand comparison against the ordinary compilation path on this particular model: Since the model is internal and does not reduce, I am happy to run further probes on it if that |
|
Thank you for the second round. The remaining failure is addressed in the latest push. Cause. The partition systems cleared the initialization equations of the model, and Fix. Tests in Related fix. Observation, not addressed here. A parameter that occurs only in an initialization equation is not a parameter of the compiled system: The heat-pump loop was not run here; the fix is verified on the reductions above. If the rerun on your side still fails, the remaining failure is likely unrelated to the partitioning. |
|
I ran Five model instances are excluded because the ordinary path fails on them as well The first two groups also break this pull request's own downstream file. Running 1. The history of an array-valued discrete state reaches the generated code unscalarized14 model instances. An observed equation that refers to an array-valued clocked variable as a using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
dt = 0.1
k = ShiftIndex(Clock(dt))
@variables uc(t) u(t) (x(t))[1:2] y(t)
eqs = [uc ~ sin(t), u ~ Sample(dt)(uc),
x[1](k) ~ u(k), x[2](k) ~ x[1](k - 1), y ~ sum(x)]
@named sys = System(eqs, t; initialization_eqs = [x[1](k - 1) ~ 0.0])
linearize_hybrid(sys, [], [y]; warn_missing_op = false)
# ERROR: UndefVarError: `xₜ₋₁` not defined in `ModelingToolkit`Writing the same output as In the package this affects every component that keeps a shift register or a coefficient product 2. An array-valued output of a discrete partition is not found in the partition system13 model instances, all of A = [0.5 0.1; -0.2 0.8]
B = reshape([1.0, 0.5], 2, 1)
C = reshape([1.0, 0.0], 1, 2)
@variables uc(t) u(t) (x(t))[1:2] (y(t))[1:1]
eqs = [uc ~ sin(t), u ~ Sample(dt)(uc), x(k) ~ A * x(k - 1) + B * [u(k)], y ~ C * x]
@named sys = System(eqs, t; initialization_eqs = [x[1](k - 1) ~ 0.0, x[2](k - 1) ~ 0.0])
linearize_hybrid(sys, [], [y[1]]; warn_missing_op = false)
# ERROR: ArgumentError: Symbol (y(t))[1] is not present in the system.
linearize_hybrid(sys, [], [y]; warn_missing_op = false)
# ERROR: ArgumentError: Symbol (y(t))[1] is not present in the system.Raised from 3. A one-sample delay whose output is read by a second equation
@variables uc(t) u(t) y(t) z(t)
eqs = [uc ~ sin(t), u ~ Sample(dt)(uc), y ~ u(k - 1), z ~ y]
@named sys = System(eqs, t; initialization_eqs = [u(k - 1) ~ 0.0])
mtkcompile(sys; additional_passes = [SynchToolkit.compile_lustre]) # compiles and simulates
linearize_hybrid(sys, [], []; warn_missing_op = false)
# ERROR: ArgumentError: Equations (1), unknowns (2), and initial conditions (2) are of
# different lengths.Raised from
4. A boolean-valued clocked variable fails with an opaque error5 model instances: @variables uc(t) u(t) y(t)
@variables s(t)::Bool
eqs = [uc ~ sin(t), u ~ Sample(dt)(uc),
s(k) ~ ifelse(s(k - 1), u(k) > -0.5, u(k) > 0.5),
y(k) ~ ifelse(s(k), 1.0, -1.0)]
@named sys = System(eqs, t; initialization_eqs = [s(k - 1) ~ false])
mtkcompile(sys; additional_passes = [SynchToolkit.compile_lustre]) # compiles and simulates
linearize_hybrid(sys, [], [y]; warn_missing_op = false)
# ERROR: TypeError: non-boolean (Float64) used in boolean contextThe same with 5. An operating point taken from a solution reintroduces the overdetermined initialization18 model instances, every hybrid closed loop in the package. This is finding 4 of my first review @variables x(t)=0.0 y(t) u(t) yd(t) ud(t)
@parameters kp = 2.0
eqs = [yd ~ Sample(dt)(y), ud ~ -kp * yd, u ~ Hold(ud), D(x) ~ -x + u, y ~ x]
@named sys = System(eqs, t)
ssys = mtkcompile(sys; additional_passes = [SynchToolkit.compile_lustre])
prob = ODEProblem(ssys, [], (0.0, 1.0))
sol = solve(prob, Tsit5())
linearize_hybrid(sys, [], [y]; op = Dict(x => 0.0), warn_missing_op = false) # quiet
linearize_hybrid(sys, [], [y]; op = LinearizationOpPoint(sol, 0.5), warn_missing_op = false)
# Warning: Initialization system is overdetermined. 1 equations for 0 unknowns.Both give What is rightWhere a partition is produced, it is correct. I checked this numerically rather than structurally. For each open-loop component the wrapper For four closed loops around a The multi-rate constructs are handled: Smaller observations
|
|
I ran this branch over the models of JuliaComputing/QuanserComponents.jl, a Furuta-pendulum package whose models have a few properties the existing coverage does not: a multibody continuous plant compiled through Environment: Julia 1.13.0, ModelingToolkit 11.45.0 at 227da52 with ModelingToolkitBase from What worksContinuous-only models are identical to A sampled-data loop around the multibody plant is right.
The multirate structure of The findings below are ordered by how much they block. 1. A partition in which a boundary signal occurs only under a shift
whenever a discrete partition reads a boundary signal only at the previous tick. The discriminating case is one term: This blocks using ModelingToolkit, DiscreteComponents, SynchToolkit, OrdinaryDiffEqTsit5
using ModelingToolkit: t_nounits as t, Shift, toggle_namespacing
using BlockComponents.Continuous: FirstOrder
using BlockComponents.Math: Gain
const dt = 0.1
# `DiscreteComponents.UnitDelay` written out: the input is read at the previous tick and
# never at the current one.
@component function Delay(; name, initial_condition = 0.0, also_current = false)
@variables u(t) [input = true] y(t) [output = true]
@parameters ic = initial_condition
rhs = also_current ? Shift(t, -1)(u) + 1.0e-9 * u : Shift(t, -1)(u)
return System([y ~ rhs], t, [u, y], [ic]; name,
initialization_eqs = [Shift(t, -1)(u) ~ ic])
end
function delay_loop(; name, block)
@named plant = FirstOrder(k = 1, T = 1)
@named smp = Sampler()
@named clock = PeriodicClock(; dt)
@named zoh = ZeroOrderHold(initial_condition = 0)
@named gain = Gain(k = -0.5)
c = block()
eqs = [connect(plant.y, smp.u)
connect(smp.y, c.u, clock.y)
connect(c.y, gain.u)
connect(gain.y, zoh.u)
connect(zoh.y, plant.u)]
return System(eqs, t, [], []; name, systems = [plant, smp, clock, zoh, gain, c])
end
# The loop compiles and simulates.
@named msim = delay_loop(block = () -> UnitDelay(; name = :c, initial_condition = 0))
ssim = mtkcompile(msim; additional_passes = [SynchToolkit.compile_lustre])
solve(ODEProblem(ssim, [ssim.plant.x => 1.0], (0.0, 1.0)), Tsit5())(1.0; idxs = ssim.plant.y)
# 0.2077879112376569
function lin(block)
@named m = delay_loop(; block)
nns = toggle_namespacing(m, false)
linearize_hybrid(m, [], [nns.plant.y]; op = Dict(nns.plant.x => 0.0),
warn_missing_op = false)
end
lin(() -> UnitDelay(; name = :c, initial_condition = 0)) # ArgumentError, as above
lin(() -> Delay(; name = :c)) # the same ArgumentError
lin(() -> Delay(; name = :c, also_current = true)) # OK
lin(() -> DiscreteDerivative(; name = :c)) # OK
lin(() -> ExponentialFilter(; name = :c)) # OKOn the package model directly: using QuanserComponents, MultibodyComponents, SynchToolkit
import DyadCompilerPasses
@named model = QuanserComponents.FurutaSwingup()
nns = ModelingToolkit.toggle_namespacing(model, false)
linearize_hybrid(model, :u_plant,
[nns.qubependulum.shoulder_angle, nns.qubependulum.elbow_angle];
op = Dict(nns.qubependulum.elbow_joint.phi => pi,
nns.qubependulum.shoulder_joint.phi => 0.0,
nns.qubependulum.elbow_joint.w => 0.0,
nns.qubependulum.shoulder_joint.w => 0.0,
nns.qubependulum.voltage => 0.0),
additional_passes = [SynchToolkit.compile_lustre],
MultibodyComponents.linsys...,
optimize = [DyadCompilerPasses.LDIV_RULE])
# ArgumentError: Equations (10), unknowns (13), and initial conditions (13) are of different lengths.2. A discrete partition with no inputs fails under
|
| Model | Result |
|---|---|
QubePendulum in a World |
equals linearize (finding 3 on the AD backend) |
FurutaPredictionModel |
equals linearize (finding 3) |
FrictionAndBackEMF, Energy |
equals linearize |
SwingupCatch loop (FurutaSwingup less GoHome) |
correct, validated against the simulation |
FurutaSwingup |
finding 1 |
FurutaHardware |
finding 1 |
FurutaFriction |
finding 2 |
FurutaIdentification |
finding 2 |
multirate loop with AlphaBetaGammaFilter and Latest |
correct three-partition split |
FurutaMPCSwingup, FurutaMPCMultirateSwingup |
finding 5 |
FurutaMPCHardware, FurutaMPCMultirateHardware |
finding 5 |
|
Thank you for the sweep. The five findings fall into three groups. 1 and 2: array-valued history and array-valued outputsBoth reproduce only without JuliaComputing/StateSelection.jl#161. With the registry 3: a one-sample delay whose output is consumedThis is a bug in the alias elimination of discrete systems and is independent of @parameters p
@variables u(t) y(t) z(t)
@named sys = System([u ~ p, y ~ u(k - 1), z ~ y], t)
ss = mtkcompile(sys)
equations(ss) # Shift(t, 1)(y(t)) ~ yₜ₊₁(t)
observed(ss) # Differential(t, 1)(y(t)) ~ p, z(t) ~ y(t), uₜ₋₁(t) ~ y(t)The alias The fix is in ModelingToolkitTearing: 4: Boolean-valued variables
This is an error rather than a warning under 5: operating point from a solutionThe values of a Smaller observationsA non-finite entry in any of |
|
Thank you for the QuanserComponents sweep. The findings are addressed in the latest push, with 1. A boundary signal used only under a shiftThis is the alias-elimination bug of the DiscreteComponents comment (finding 3 there), now 2. A discrete partition with no inputs under
|
|
Retried the whole QuanserComponents sweep on 3acd78a. Same environment as before, in two variants: one with the released StateSelection 1.11.1, one with JuliaComputing/StateSelection.jl#163 at 70685b78a9 and its ModelingToolkitTearing developed in. Everything is confirmed fixed, and the models that were blocked now linearize and give the right answers. Findings 1 to 41, the boundary signal used only under a shift. Confirmed as an upstream bug and confirmed fixed by #163. The same script both ways:
With #163, 2, the discrete partition without inputs. Fixed. 3, where it previously fell back and carried 2.9e-8. The same holds one level up: the continuous partition of 4, the input-derivative diagnostic. Fixed; the message names the bound parameter: One nit: raised from 5, the MPC models. Unchanged and out of scope, as agreed: The operating points. The new check does exactly what was needed. On That names the symbol and both values, and is what I had to reconstruct by hand from New checks the fixes made possibleWith Analysis points on the real model. Operating point from a solution of the swing-up run. Taking the 10 s swing-up of
Two poles of that loop sit at exactly The rest of the sweepUnchanged and still right, now exactly rather than nearly:
Summary by model on 3acd78a with #163:
|
|
Rerun of the full sweep on First, my error on findings 1 and 2. Your diagnosis is exactly right and the cause is worth Corrected environment, used for everything below: Julia 1.13.0, ModelingToolkit at The sweep on
|
Verdict: changes needed before merge (reviewed at Adds Risk assessment
Findings
Push a fix and the PR is reviewed again automatically at the new head. 🤖 Posted by an AI agent — harness: Claude Code · model: claude-opus-5-5[1m] (fleet master); review by Devin CLI (Mac) / fusion-claude-opus-5-5-high-sidekick-swe-2-medium |
|
This is a draft pr for a reason, can you please stop the bot spam? |
|
A review was ran on every PR of the month, it's no big deal and just one posting? Every PR needs a risk assessment and your posts didn't have one. |
…stems `linearize` accounts for the continuous equations only, so a model with a synchronous controller linearizes to the open-loop plant. `linearize_hybrid` splits the model into its clock partitions with the clock inference used by `mtkcompile`, compiles every partition with the standard structural simplification and differentiates the generated functions of each partition numerically (ForwardDiff through DifferentiationInterface, with an automatic `AutoFiniteDiff` retry). The result holds one linear system per partition with its state, input and output variables and sample interval, and the connections between partitions across clock boundaries, so that the sampled-data loop can be assembled with `ControlSystemsBase.feedback` or `RobustAndOptimalControl.connect`. User inputs are bound to time-independent parameters before clock inference. Boundary terms (`Sample`, `Hold`, clock changes and `SampleTime`) become plain parameters of their partition. `get_sensitivity`, `get_comp_sensitivity` and `get_looptransfer` accept `hybrid = true`, and `linearize_hybrid` accepts analysis points. `linearization_function` is split into the compilation step and `_linearization_function_compiled` without behaviour change. ModelingToolkitBase: `find_all_parameter_equations` treats operator terms that are unknowns or parameters of the system as atomic symbols, the incomplete initialization branch of `InitializationProblem` allows the parameter-only equations it appends, and shifts are no longer distributed into operators other than `Shift` and `Differential`. Closes #5166. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01QBgBrXDidWG1G2KF64gAte
Partition systems no longer inherit the assertions of the model, keep only the parameters they use together with their boundary and input parameters, and keep the initial values of their unknowns and parameters only. Values of signals crossing a clock boundary that the operating point does not provide are evaluated from the partition they originate from, so an operating point only needs the state variables; operating-point values of observed variables that are expressions of parameters only are not turned into initialization constraints. The finite-difference fallback is limited to errors of automatic differentiation, the preprocessing mirrors `mtkcompile`, and `show` reports the state dimension. The downstream test file is not registered in the test groups, since its dependencies are registered in DyadRegistry only, and the non-public ModelingToolkitBase names reached from the new file are listed in the QA allowlist. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01QBgBrXDidWG1G2KF64gAte
The partition systems cleared the initialization equations of the model, and the initialization problem of the continuous partition was built with `algebraic_only = true`, which excludes initialization equations. Parameters bound to `missing` and constant states (`der(x) = 0`, compiled to a parameter) that the model determines through initialization equations were therefore left undetermined: the underdetermined initialization problem returned guessed values, or failed with `Could not evaluate value of parameter` when the structurally singular problem dropped the undetermined symbols to parameters. `_compile_partition!` now keeps the initialization equations whose symbols are all known to the partition, and keeps the parameters they refer to. Before the continuous partition is linearized, the equations that determine at least one unknown or `missing`-bound parameter whose value the operating point does not fix are passed to the initialization problem through a new `initialization_eqs` keyword of `_linearization_function_compiled`; observed variables are expanded to the symbols they depend on for this selection. An operating-point value takes precedence over such an equation. For a discrete partition, equations of the form `x(k-1) ~ expr` and `x ~ expr` provide the history values. Parameters bound to `missing` that nothing determines take their guess, or zero with the missing-value warning. `LinearizationFunction` no longer differentiates with respect to an empty input vector, so `linearize(sys, [], outputs)` works and returns `B` and `D` with zero columns. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01QBgBrXDidWG1G2KF64gAte
`clock_boundary(hl, from, to)` collects the connections from one partition of a `HybridLinearization` to another as the index vectors `outputs` (into the outputs of the partition the signals leave) and `inputs` (into the inputs of the partition they enter), which are the index sets required by the advanced interface of `ControlSystemsBase.feedback`. The documentation and the tests assemble the sampled-data loops from these index sets instead of hard-coded indices or filtered connection lists, and the documentation states the relation to `starprod` and `lft`. The documentation also describes how to check whether the signals entering a sampler carry content above the Nyquist frequency of the sampler: the gain from all inputs of the continuous partition, the external inputs and the held signals, to the sampled outputs is evaluated above the Nyquist frequency and compared with the gain below it. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01QBgBrXDidWG1G2KF64gAte
For a signal leaving a discrete partition for a partition on a slower clock, the gain of the fast partition alone in the band between the two Nyquist frequencies is a sufficient but conservative check, since it ignores the attenuation of the continuous partition. The documentation now describes the check on the composite of the continuous partition, discretized at the fast rate, and the fast partition, from all inputs of the composite, including the held outputs of the slow partition, to the boundary signal, with a multi-rate example. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01QBgBrXDidWG1G2KF64gAte
…tries in `linearize_hybrid` A partition whose unknowns, inputs or outputs are Boolean-valued now produces an error naming the variables instead of a `TypeError` raised from the generated code. An operating point given as a `LinearizationOpPoint` provides the values of the unknowns and parameters of each partition only; the values a solution holds for observed variables repeat the constraints the unknowns satisfy and made the initialization of the continuous partition overdetermined. Non-finite entries in the matrices of a partition, which arise when an equation is not differentiable at the operating point, are reported by a warning. The downstream tests gain a loop with a `UnitDelay` whose output is consumed. Compiling its discrete partition requires `StateSelection.var_derivative!` to build the shifted form of an alias target as a `Shift` for discrete systems, a fix in ModelingToolkitTearing. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01QBgBrXDidWG1G2KF64gAte
…t-free discrete partitions The inline solution of linear subsystems derives the element type of its buffers from the unknowns and the first registered input of the system. The input and boundary parameters of a clock partition were not registered, so the buffers stayed `Float64` when only the parameters carried dual numbers and the Jacobians with respect to the inputs of a multibody plant fell back to finite differences. `_mtkcompile!` gains the keyword argument `input_parameters`, parameters that represent inputs and are registered as such, and the initialization system accepts inputs that are time-independent parameters. `LinearizationFunction` stores its input symbols, so that the input-derivative diagnostic names them also when the inputs are parameters rather than variables marked as inputs. Discrete partitions without inputs no longer evaluate Jacobians with respect to an empty input vector, which finite differences reject. A warning reports when the initialization of the continuous partition does not retain the requested value of a differential unknown, which happens when the operating point also fixes signals the model determines from the state. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01QBgBrXDidWG1G2KF64gAte
…entries and drop two unneeded changes An operating point that fixes the output of the continuous partition alongside its state makes the initialization overdetermined, which ModelingToolkit reports; the test expected no warning, which the daemon-based test runs could not verify since they forward stdout only. The new public names are added to the canonical `@docs` block of the API reference, which the documentation build requires. The changes to `find_all_parameter_equations` and to the incomplete-initialization branch of `InitializationProblem` are removed: they were needed when the boundary terms of a partition were parameters of the partition system, and are not since those terms are replaced by plain parameters. The hybrid and downstream tests pass without them. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01QBgBrXDidWG1G2KF64gAte
3acd78a to
9011757
Compare
|
Addressed at the new head, which is rebased onto
🤖 Generated with Claude Code |
This branch has not been deployed
2 participants
Closes #5166 (second approach: separate linearization of all clock partitions).
Summary
linearizeaccounts for the continuous equations of a model only, so a model with a synchronous (clocked) controller linearizes to the open-loop plant. This PR addslinearize_hybrid(sys, inputs, outputs; op, t, autodiff, ...), which splits the model into its clock partitions with the clock inference already used bymtkcompile, compiles every partition separately with the standard structural simplification, and differentiates the generated functions of each partition numerically (ForwardDiff through DifferentiationInterface, with an automaticAutoFiniteDiffretry). No symbolic differentiation is involved, so models calling arbitrary Julia functions are supported.The result is a
HybridLinearizationholding oneClockPartitionLinearizationper partition (continuous partition first) and a list of connections. Each partition carries the state-space matrices, the state, input and output variable lists, the clock and the sample interval. The inputs and outputs of a partition consist of the user-specified signals belonging to that partition (first, in the user's order;nu_user,ny_user) followed by the signals crossing clock boundaries into and out of the partition. The connections(; from, output, to, input)are the index pairs required by the advanced interface ofControlSystemsBase.feedbackor the signal names forRobustAndOptimalControl.connect.clock_boundary(hl, from, to)collects the connections from one partition to another as the index vectorsoutputsandinputsfor the advanced interface ofControlSystemsBase.feedback.For a discrete partition, the state consists of the values at the previous tick and the matrices describe the update at a tick,
x(k) = A x(k-1) + B u(k),y(k) = C x(k-1) + D u(k), whose transfer function isC(zI - A)⁻¹B + D. This is the form produced by the discrete-system compilation of ModelingToolkit (history variablesxₜ₋₁as unknowns, observed equations shifted back to the current tick).Implementation notes
u ~ u_inputbefore clock inference. This keeps the input a variable of its partition and avoids the back-shifting of time-dependent parameters that the discrete reassembly performs.Sample(...),Hold(...),SynchToolkit.Latest(...),SampleTime()) are turned into parameters the waymtkcompile!does for the continuous partition of a hybrid system, and then replaced by plain parameters named after the signal and the operator (u_hold,y_sample,lat₊u_latest,sample_time), so that the compilation and code generation of a discrete partition need not treat operators defined outside ModelingToolkit as atomic symbols. The variable a boundary term is applied to, after stripping nestedSample/Hold, identifies the source partition. Boundary terms that no equation of their partition reads (theHold(w)of a nested clock changeSample(clk)(Hold(w))) are compiled but not exposed. The operating-point value of theSampleTime()parameter is the sample interval of the partition's clock.ModelingToolkitTearing.substitute_sample_timeis not used: itssampletimehelper matches everyAbstractClockin its first@matcharm and therefore never returns a sample interval (worth fixing in ModelingToolkitTearing separately).linearization_functionis split into the compilation step and_linearization_function_compiled, which operates on a compiled system and is reused for the continuous partition (initialization, DAE handling andLinearizationOpPointsupport included). Behaviour is unchanged.get_sensitivity,get_comp_sensitivityandget_looptransferaccepthybrid = trueand return aHybridLinearization;linearize_hybridaccepts analysis points as inputs and outputs. The transformation part ofget_linear_analysis_functionis factored intolinear_analysis_transform.LinearizationOpPoint(sol, t); timeseries (callback-saved discrete) parameters are evaluated att. The values of boundary signals that are not inopare evaluated from the operating point of the partition they originate from (continuous partition first), so an operating point only needs the state variables of the partitions; history variables take the value of the variable they are the history of. Remaining missing values default to zero with a warning.missing-bound parameter whose value the operating point does not fix, are passed to the initialization problem of the continuous partition. Parameters bound tomissingand constant states (der(x) = 0, compiled to a parameter) that the model determines through initialization equations therefore take their values from them, and an operating-point value for such a symbol takes precedence. For a discrete partition, initialization equations of the formx(k-1) ~ exprorx ~ exprprovide the history values ofx. Parameters bound tomissingthat nothing determines take their guess, or zero with the missing-value warning.LinearizationOpPoint, the values of the unknowns and parameters of each partition are taken; the values the solution holds for observed variables are not used, since they would repeat the constraints the unknowns already satisfy and made the initialization of the continuous partition overdetermined. Values given explicitly in the dictionary of theLinearizationOpPointare treated like a dictionary operating point.x0field of the partition holds the values used._mtkcompile!gains the keyword argumentinput_parameters). The inline solution of linear subsystems (inline_linear_sccs, used by MultibodyComponents) derives the element type of its buffers from the unknowns and the first registered input, so without the registration the buffers stayedFloat64when only the parameters carried dual numbers and the input Jacobians of a multibody plant fell back to finite differences. The initialization system accepts inputs that are time-independent parameters.Prerequisites in ModelingToolkitTearing
DiscreteComponents.DiscreteStateSpace,MovingAverageFilter,FIRFilter,MPCComponents.ACADOSMPC, ...) require Shift unscalarized array variables with their scalarized elements JuliaComputing/StateSelection.jl#161.y ~ u(k - 1)is followed by an aliasz ~ y(DiscreteComponents.UnitDelaywith its output connected) require Build shifted variables invar_derivative!for discrete systems JuliaComputing/StateSelection.jl#163:StateSelection.var_derivative!built the shifted form of an alias target as aDifferentialalso for discrete systems, which left the compiled partition with two unknowns and one equation. This is a bug ofmtkcompileon such discrete systems independent of this PR.Related pull requests
Three fixes found while developing this feature are submitted separately, each with a regression test that fails on
master; this branch carries copies of them until they are merged, after which it is rebased and the ModelingToolkitBase floor is raised to a release that contains them:linearizewithout inputs and name differentiated inputs from the linearization function #5225:linearize(sys, [], outputs)failed in ForwardDiff with an invalid tag, and the input-derivative diagnostic indexedinputs(sys), which is empty when the inputs of a compiled system are parameters not marked as inputs.LinearizationFunctionnow stores its input symbols. Discrete partitions without inputs are guarded in the same way, since finite differences reject an empty vector of differentiation variables.ShiftandDifferential#5226:distribute_shiftno longer distributes shifts into operator terms other thanShiftandDifferential. An operator term such as a clock change defined by another package denotes the value at the current tick; distributing the back-shift of the discrete reassembly intoLatest(x, init)produced a reference to a nonexistentxₜ₋₁.Tests
test/linearize_hybrid.jl(InterfaceII): continuous-only model equal tolinearize; single-clock loop assembled withfeedbackand compared with the closed-loop recursion; discrete-only model; second-order difference equation and an array-valued discrete state space with array parameters; all-linear multi-rate model with a nested clock change; operating points as a dictionary and from a solution;AutoFiniteDiffand the fallback for a registered function that rejects dual numbers, including a discrete partition without inputs; warnings for events; the error for Boolean variables, the input-derivative diagnostic, the warning for non-finite entries and the warning for requested state values that the initialization does not retain; analysis points on continuous and discrete connections includingget_sensitivity(...; hybrid = true).test/downstream/linearize_hybrid_components.jl(Downstream, uses DiscreteComponents, BlockComponents and SynchToolkit):DiscretePIDStandardloop whose assembled sampled-data closed loop matches the SynchToolkit simulation;DiscreteStateSpaceSISO and MIMO with array inputs and outputs compared withss(A, B, C, D, Ts);MovingAverageFilter; aUnitDelayfollowed by a gain in a loop; a multi-rate loop withLatest; operating point from a SynchToolkit solution of a nonlinear loop, without initialization warnings; analysis points.The downstream test file is not registered in the test groups, since
DiscreteComponentsandSynchToolkitare registered in DyadRegistry only; its header describes the environment it runs in. ItsDiscreteStateSpacecases shift unscalarized array variables and were run with JuliaComputing/StateSelection.jl#161; the in-repo array tests pass with the released ModelingToolkitTearing.Documentation
New section "Hybrid systems with clock partitions" in
docs/src/basics/Linearization.md, including the assembly of the sampled-data loop from the boundary index sets, its relation tostarprod/lft, and checks of the content of the signals entering a sampler above its Nyquist frequency, for the continuous partition and, through the composite of the continuous and the faster partitions, for signals entering a slower clock partition; a note in the linear analysis tutorial.🤖 Generated with Claude Code
https://claude.ai/code/session_01QBgBrXDidWG1G2KF64gAte