Repository navigation
Add simulation option 'jacobian_mode' (auto | dd | fd) and a kink guard for the finite-difference Jacobian - #424
Draft
hubertus65 wants to merge 3 commits into
Draft
hubertus65 wants to merge 3 commits into
hubertus65 wants to merge 3 commits into
Conversation
…Jacobian When 'with_jacobian' is in effect and the FMU provides directional derivatives, PyFMI evaluates the Jacobian with one fmi2GetDirectionalDerivative call per colour group of the sparsity pattern. On FMUs where one such call costs more than an rhs evaluation (measured 1.7x on OCT exports of 140-320 state models), coloured forward differences give the same Jacobian for a third of the cost, and because the colouring is bounded by the densest Jacobian row, Jacobian evaluation was 55-85 % of the solver time on those models. Measured with CVode: absorptionplant (226 states) 164 -> 89 s; Radau5ODE 529 -> 206 s; same step counts, events and accuracy (deviation <= 4e-4 at rtol 1e-6). 'jacobian_mode' = "auto" (default) times one rhs and one directional-derivative call at initialization and picks "fd" when the derivative is the more expensive one, except for RodasODE (a Rosenbrock method needs the exact Jacobian), the sparse linear solver, and rtol < 1e-7 (forward-difference truncation error ~sqrt(eps) then stalls Newton, measured on a 318-state plant at rtol 1e-8). "dd" and "fd" force either. Without directional derivatives the mode is always "fd". The setting is scoped to the simulation: the model's force_finite_differences attribute is restored afterwards. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
A forward difference taken across a kink of the model -- a limiter, a region boundary of a property function -- returns jump/h (1e8 and more) where the derivative on this side is small. Newton-based integrators converge on the corrupted iteration matrix and reuse it for many steps. Measured on a 318- state combined-cycle plant at rtol 1e-6: two such entries in 200 Jacobian evaluations (d der(pump.s)/d s = -3.2e8 vs 0 at the end of the steam-table failure zone; d der(drum.p)/d hl = -1.3e8 vs -0.18 at t = 0) moved 91 of 308 plant states more than 1e-3 off the rtol 1e-8 reference. The guard (FMUModelME2.fd_kink_guard, default True) works per colour group of _estimate_directional_derivative: - flag: an entry that grew more than 1e3x relative to the previous Jacobian and is not negligible relative to it marks the group as suspect; the first evaluation checks every group. Kinks are transient, so this is nearly free. - confirm: a suspect group is re-differenced with a 10x smaller step; entries whose quotient changes by more than 10 % are kinks (a jump scales with 1/h, a smooth function does not). - fix: kink entries take the FMU's directional derivative for the group when it has one, else the smaller in magnitude of the forward and the backward difference (the side without the jump). Model variables are restored after each check; statistics in _fd_guard_stats. With the guard, CVode with the FD Jacobian on the plant matches the DD run's accuracy (plant median 5.0e-6 vs 9.8e-6, 4 vs 1 states beyond 1e-3, 23 events) at 5458 instead of 14353 steps and 320 instead of 511 Jacobians, each half the cost: 229 s instead of 340 s. On a 226-state plant and a 141-state vehicle model the guard changes nothing (accuracy and time within noise). FMI3's separate finite-difference routine is unchanged. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PyrxGeWuwCmAH9LishMW6p
…tion's help The default stays "auto" (as introduced); this records what the kink guard protects against and asserts the default in the option test. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PyrxGeWuwCmAH9LishMW6p
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
A new simulation option
jacobian_mode("auto"|"dd"|"fd") choosing how the PyFMIJacobian for an ME FMU is formed, and a kink guard that makes the finite-difference path
trustworthy enough to be the automatic choice. Three commits.
9823fa8"auto"times onefmi2GetDirectionalDerivativecall against one rhs evaluation at initialization and picks the cheaper; scoped to the simulation, so a later manual Jacobian call is unaffected."dd"is the previous behaviour. The mode actually used is on the algorithm object asjacobian_mode.f352f04FMUModelME2.fd_kink_guard: the kink guard on the coloured finite-difference JacobianWhy
On the 12 ME FMUs of this study the Jacobian evaluation is 55–85 % of the solver time, because
one directional-derivative call costs about 2× an rhs evaluation and the colouring is bounded by
the densest Jacobian row. Coloured finite differences cost one rhs per colour group.
"auto"declines to use finite differences when the tolerance is tight enough that theforward-difference truncation error (~
sqrt(eps)) would stall Newton (rtol < 1e-7), when thesolver needs an exact Jacobian (Rosenbrock methods), and when a directional derivative is cheap
anyway.
The kink guard, and why the option is not enough without it
Plain coloured finite differences are not safe by default, which is the interesting part of
this PR. On a 318-state combined-cycle plant, two entries per ~200 Jacobians had
jump/h ≈ 1e8— a difference taken across a steam-table region boundary — and that corrupted the Newton matrix
for many steps afterwards: 91 of 308 plant states ended more than 1e-3 off the reference, at
rtol 1e-6 and 5e-7 but not at 2e-7 or 1e-7, so no check at initialization detects it and a tighter
tolerance appears to "fix" it.
The guard: an entry that jumps more than 1e3× against the previous Jacobian flags its colour
group; a flagged group is re-differenced with
h/10; entries whose quotient then moves more than10 % are kinks and take the directional derivative if the FMU has one, otherwise the
smaller-magnitude one-sided difference.
With it, the finite-difference run matches the directional-derivative run on that model (median
deviation 5.0e-6) in 229 s instead of 340 s.
Behaviour change
FMUs whose directional derivatives cost more than an rhs now get finite differences by default.
options["jacobian_mode"] = "dd"restores the previous behaviour per run. FMUs withoutdirectional derivatives are unaffected except that the guard adds one extra Jacobian on the first
evaluation and then only re-differences flagged groups.
Review notes
f352f04touchesfmi2.pyx/fmi2.pxd; the rest isfmi_algorithm_drivers.py."dd"and then back to"auto"have been squashed out; what remains is one decision with its documentation.
test_jacobian_mode_option(mode selection under every combination of capability,solver and tolerance) and the kink-guard tests in
tests/test_fmi2.py.🤖 Generated with Claude Code
https://claude.ai/code/session_01PyrxGeWuwCmAH9LishMW6p