From 578932c484fe48a3807b88b162b229c255b05a57 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Thu, 25 Jun 2026 13:46:31 +0000 Subject: [PATCH 1/3] Add `report_initialization_nullspace` init-Jacobian diagnostic MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Underdetermined or rank-deficient initialization systems are a common source of intermittent ("flaky") initialization failures: the solve lands on different solutions depending on the solver and on the nondeterministic ordering of the assembled init unknowns/equations, and some orderings hit a singular realization and throw while others succeed. `report_initialization_nullspace(prob)` makes the offending degrees of freedom explicit. It evaluates the Jacobian of the initialization residual at the initial guess (via ForwardDiff), computes its SVD, and reports the unknowns that span the null space — the directions along which the initialization is free to move. Pinning those unknowns (or adding constraints) makes the init well-posed and deterministic. `prob` may carry initialization data (an `ODEProblem`/`DAEProblem` built from a `System`) or be an initialization `NonlinearProblem`/`NonlinearLeastSquaresProblem` directly. The per-variable weight reported is the diagonal of the null-space projector (squared row norm over an orthonormal null-space basis), which is invariant to the arbitrary basis chosen within the null space. A nullity of 0 means the init Jacobian has full column rank at the guess and is locally well-posed. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01FSDc4tFYygs26KhqZBACLS --- src/ModelingToolkit.jl | 1 + src/initialization.jl | 124 +++++++++++++++++++++++++++++++ test/group_initialization.jl | 1 + test/initialization_nullspace.jl | 36 +++++++++ 4 files changed, 162 insertions(+) create mode 100644 test/initialization_nullspace.jl diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index a2f545bce3..8a1dd61cd0 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -152,6 +152,7 @@ end @reexport using .StructuralTransformations export SemilinearODEFunction, SemilinearODEProblem +export report_initialization_nullspace export alias_elimination export linearize, linearization_function, LinearizationProblem, linearization_ap_transform diff --git a/src/initialization.jl b/src/initialization.jl index e75ec4820e..9a676acf87 100644 --- a/src/initialization.jl +++ b/src/initialization.jl @@ -47,3 +47,127 @@ function MTKBase.get_initialization_problem_type( NonlinearLeastSquaresProblem end end + +""" + report_initialization_nullspace(prob; rtol = 1e-8, atol = 0.0, var_threshold = 1e-3, verbose = true) + +Diagnose rank deficiency of a system's initialization problem by inspecting the null +space of its residual Jacobian, and report which unknowns span that null space. + +When the initialization system is underdetermined or otherwise rank deficient, the +initialization solve can converge to different solutions depending on the solver and on +the (often nondeterministic) ordering of the assembled unknowns/equations. This is a +common source of intermittent initialization failures: some orderings land on a singular +realization and throw, others succeed. This utility makes the offending degrees of freedom +explicit. + +It evaluates the Jacobian of the initialization residual at the initial guess `u0` (via +`ForwardDiff`), computes its singular value decomposition, and reports the unknowns that +participate in the null space — i.e. the directions along which the initialization is free +to move. Pinning those unknowns (or adding equations that constrain them) makes the +initialization well-posed and deterministic. + +`prob` may be a problem that carries initialization data (e.g. an `ODEProblem`/`DAEProblem` +built from a `System`), or an initialization `NonlinearProblem` / +`NonlinearLeastSquaresProblem` directly. + +# Keyword arguments + - `rtol`: a singular value `σ` is treated as zero when `σ ≤ max(atol, rtol * σmax)`, + where `σmax` is the largest singular value. Increase it to also surface + near-singular directions. + - `atol`: absolute singular-value threshold (see `rtol`). + - `var_threshold`: only unknowns whose null-space participation exceeds this value are + reported. + - `verbose`: print a human-readable report. + +# Returns + +A `NamedTuple` `(; jacobian, singular_values, nullity, variables)`. `variables` is a vector +of `unknown => weight` pairs sorted by decreasing `weight`. The weight, in `[0, 1]`, is the +squared norm of the unknown's row in an orthonormal basis of the null space (the diagonal +of the null-space projector): how strongly that unknown participates in the +underdetermined directions, independent of the arbitrary basis chosen within the null +space. A `nullity` of `0` means the initialization Jacobian has full column rank at `u0` +and the initialization is locally well-posed. + +!!! note + The Jacobian is evaluated at a single point (`u0`), so this reports the *local* rank + structure there. A structurally well-posed initialization can still be numerically + rank deficient at a particular operating point, and vice versa. +""" +function report_initialization_nullspace( + prob; rtol = 1e-8, atol = 0.0, var_threshold = 1e-3, verbose = true) + empty_result = (; jacobian = nothing, singular_values = Float64[], + nullity = 0, variables = Pair[]) + iprob = _initialization_problem(prob) + if iprob === nothing + verbose && + @info "No initialization problem to analyze: the system is fully determined by its initial conditions." + return empty_result + end + u0 = state_values(iprob) + if u0 === nothing || isempty(u0) + verbose && @info "The initialization problem has no unknowns to solve for." + return empty_result + end + u0 = collect(float.(u0)) + p = parameter_values(iprob) + f = iprob.f + nresid = f.resid_prototype !== nothing ? length(f.resid_prototype) : length(u0) + residual = if SciMLBase.isinplace(f) + u -> (du = similar(u, nresid); f(du, u, p); du) + else + u -> f(u, p) + end + J = ForwardDiff.jacobian(residual, u0) + n = length(u0) + fact = svd(J; full = true) + S = fact.S + σmax = isempty(S) ? zero(eltype(S)) : maximum(S) + threshold = max(atol, rtol * σmax) + # A right-singular vector is a null direction when its singular value is below the + # threshold; columns of `V` beyond `length(S)` (present when n > nresid) have an + # implied singular value of zero and are always null directions. + is_null = [(j <= length(S) ? S[j] : zero(σmax)) <= threshold for j in 1:n] + nullity = count(is_null) + syms = variable_symbols(iprob) + # Participation = diagonal of the null-space projector V_null * V_null', i.e. the + # squared row norm of each unknown over an orthonormal null-space basis. This is + # invariant to the arbitrary choice of basis within the null space. + weights = nullity == 0 ? zeros(n) : vec(sum(abs2, @view(fact.V[:, is_null]); dims = 2)) + variables = sort( + [syms[i] => weights[i] for i in 1:n if weights[i] > var_threshold]; + by = last, rev = true) + if verbose + println("Initialization Jacobian null-space report") + println(" residual Jacobian: $(nresid)×$(n), rank ≈ $(n - nullity), nullity ≈ $nullity") + if !isempty(S) + k = min(5, length(S)) + println(" smallest $k singular value(s): ", + round.(S[(end - k + 1):end], sigdigits = 3)) + end + if nullity == 0 + println(" Full column rank at u0 — initialization is locally well-posed.") + else + println(" Unknowns spanning the null space (participation weight ∈ [0, 1]):") + for (s, w) in variables + println(" ", lpad(string(round(w, digits = 4)), 8), " ", s) + end + isempty(variables) && + println(" (no individual unknown exceeds var_threshold = $var_threshold)") + end + end + return (; jacobian = J, singular_values = S, nullity, variables) +end + +# Return the initialization `NonlinearProblem`/`NonlinearLeastSquaresProblem` carried by +# `prob`, or `prob` itself if it is already a nonlinear problem, or `nothing` if there is +# no initialization problem to analyze. +function _initialization_problem(prob) + prob isa SciMLBase.AbstractNonlinearProblem && return prob + f = prob.f + if hasproperty(f, :initialization_data) && f.initialization_data !== nothing + return f.initialization_data.initializeprob + end + return nothing +end diff --git a/test/group_initialization.jl b/test/group_initialization.jl index 70460390c2..4d011487b8 100644 --- a/test/group_initialization.jl +++ b/test/group_initialization.jl @@ -2,5 +2,6 @@ include("shared/mtktestset.jl") @mtktestset("Guess Propagation", "guess_propagation.jl") @safetestset "Hierarchical Initialization Equations" include("hierarchical_initialization_eqs.jl") +@safetestset "Initialization Jacobian null space" include("initialization_nullspace.jl") @mtktestset("InitializationSystem Test", "initializationsystem.jl") @mtktestset("Initial Values Test", "initial_values.jl") diff --git a/test/initialization_nullspace.jl b/test/initialization_nullspace.jl new file mode 100644 index 0000000000..7aa510e365 --- /dev/null +++ b/test/initialization_nullspace.jl @@ -0,0 +1,36 @@ +using ModelingToolkit, Test +using ModelingToolkit: t_nounits as t, D_nounits as D + +# Underdetermined initialization: a point constrained to the unit circle, with no +# initial condition pinning where on the circle (or with what velocity) it starts. +# The initialization Jacobian is rank deficient, so the report must flag a nonempty +# null space. +@variables x(t) y(t) +@named sys = System([D(x) ~ y, 0 ~ x^2 + y^2 - 1], t) +sys = mtkcompile(sys) +prob = ODEProblem(sys, [], (0.0, 1.0); guesses = [x => 0.6, y => 0.8], + warn_initialize_determined = false) +res = report_initialization_nullspace(prob; verbose = false) +@test res.nullity ≥ 1 +@test !isempty(res.variables) +@test res.jacobian !== nothing +@test all(0 .≤ last.(res.variables) .≤ 1 + 1e-8) + +# Determined initialization: an algebraic unknown uniquely fixed by a pinned state has +# a full-column-rank initialization Jacobian and no null space. +@variables u(t) w(t) +@named sys2 = System([D(u) ~ -u, 0 ~ w - exp(u)], t) +sys2 = mtkcompile(sys2) +prob2 = ODEProblem(sys2, [u => 1.0], (0.0, 1.0); guesses = [w => 1.0]) +res2 = report_initialization_nullspace(prob2; verbose = false) +@test res2.nullity == 0 +@test isempty(res2.variables) + +# A system fully determined by its initial conditions has no initialization problem to +# analyze; the utility handles that gracefully rather than erroring. +@variables a(t) b(t) +@named sys3 = System([D(a) ~ -a, D(b) ~ -b], t) +sys3 = mtkcompile(sys3) +prob3 = ODEProblem(sys3, [a => 1.0, b => 2.0], (0.0, 1.0)) +res3 = report_initialization_nullspace(prob3; verbose = false) +@test res3.nullity == 0 From b237da57be1d23eeafe1e6c022316958f6a23f3c Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Fri, 26 Jun 2026 04:49:56 +0000 Subject: [PATCH 2/3] Also report redundant init equations; rename to `analyze_initialization_jacobian` MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Extend the initialization diagnostic to report both sides of the rank deficiency of the initialization residual Jacobian: - underdetermined unknowns — its (right) null space, the directions the initialization is free to move along; and - redundant equations — its left null space, combinations of initialization equations whose Jacobian rows are linearly dependent at `u0` and therefore do not locally constrain any further degree of freedom. The latter explains the common, confusing situation where an initialization has more equations than unknowns yet is still underdetermined. Renamed from `report_initialization_nullspace` to `analyze_initialization_jacobian` now that it analyzes the full rank structure (rank, nullity, redundancy) rather than only the null space. The return value gains `rank`, `redundancy`, and `redundant_equations` (a vector of `equation => weight` pairs), and the per-unknown result is renamed `underdetermined_unknowns`. Redundant equations are mapped back to the symbolic initialization equations when available. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01FSDc4tFYygs26KhqZBACLS --- src/ModelingToolkit.jl | 2 +- src/initialization.jl | 142 ++++++++++++++--------- test/group_initialization.jl | 2 +- test/initialization_jacobian_analysis.jl | 53 +++++++++ test/initialization_nullspace.jl | 36 ------ 5 files changed, 145 insertions(+), 90 deletions(-) create mode 100644 test/initialization_jacobian_analysis.jl delete mode 100644 test/initialization_nullspace.jl diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index 8a1dd61cd0..f27886e602 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -152,7 +152,7 @@ end @reexport using .StructuralTransformations export SemilinearODEFunction, SemilinearODEProblem -export report_initialization_nullspace +export analyze_initialization_jacobian export alias_elimination export linearize, linearization_function, LinearizationProblem, linearization_ap_transform diff --git a/src/initialization.jl b/src/initialization.jl index 9a676acf87..962b8dc181 100644 --- a/src/initialization.jl +++ b/src/initialization.jl @@ -49,56 +49,67 @@ function MTKBase.get_initialization_problem_type( end """ - report_initialization_nullspace(prob; rtol = 1e-8, atol = 0.0, var_threshold = 1e-3, verbose = true) + analyze_initialization_jacobian(prob; rtol = 1e-8, atol = 0.0, threshold = 1e-3, verbose = true) -Diagnose rank deficiency of a system's initialization problem by inspecting the null -space of its residual Jacobian, and report which unknowns span that null space. +Diagnose rank deficiency of a system's initialization problem by inspecting the singular +value decomposition of its residual Jacobian, and report both the **unknowns** that span +its (right) null space and the **equations** that are locally redundant (its left null +space). -When the initialization system is underdetermined or otherwise rank deficient, the -initialization solve can converge to different solutions depending on the solver and on -the (often nondeterministic) ordering of the assembled unknowns/equations. This is a -common source of intermittent initialization failures: some orderings land on a singular -realization and throw, others succeed. This utility makes the offending degrees of freedom -explicit. +When the initialization system is rank deficient the solve can converge to different +solutions, or fail, depending on the solver and on the (often nondeterministic) ordering of +the assembled unknowns/equations — a common source of intermittent initialization failures: +some orderings land on a singular realization and throw, others succeed. This utility makes +the offending degrees of freedom explicit: -It evaluates the Jacobian of the initialization residual at the initial guess `u0` (via -`ForwardDiff`), computes its singular value decomposition, and reports the unknowns that -participate in the null space — i.e. the directions along which the initialization is free -to move. Pinning those unknowns (or adding equations that constrain them) makes the -initialization well-posed and deterministic. + - **Underdetermined unknowns** (right null space): directions along which the + initialization is free to move. Pinning these unknowns (or adding equations that + constrain them) makes the initialization well-posed and deterministic. + - **Redundant equations** (left null space): combinations of initialization equations + whose Jacobian rows are linearly dependent at `u0`, i.e. equations that do not locally + constrain any additional degree of freedom. These explain why an initialization can + have more equations than unknowns yet still be underdetermined. -`prob` may be a problem that carries initialization data (e.g. an `ODEProblem`/`DAEProblem` -built from a `System`), or an initialization `NonlinearProblem` / -`NonlinearLeastSquaresProblem` directly. +It evaluates the Jacobian of the initialization residual at the initial guess `u0` (via +`ForwardDiff`) and computes its SVD. `prob` may be a problem that carries initialization +data (e.g. an `ODEProblem`/`DAEProblem` built from a `System`), or an initialization +`NonlinearProblem`/`NonlinearLeastSquaresProblem` directly. # Keyword arguments - `rtol`: a singular value `σ` is treated as zero when `σ ≤ max(atol, rtol * σmax)`, - where `σmax` is the largest singular value. Increase it to also surface - near-singular directions. + where `σmax` is the largest singular value. Increase it to also surface near-singular + directions. - `atol`: absolute singular-value threshold (see `rtol`). - - `var_threshold`: only unknowns whose null-space participation exceeds this value are + - `threshold`: only unknowns/equations whose participation exceeds this value are reported. - `verbose`: print a human-readable report. # Returns -A `NamedTuple` `(; jacobian, singular_values, nullity, variables)`. `variables` is a vector -of `unknown => weight` pairs sorted by decreasing `weight`. The weight, in `[0, 1]`, is the -squared norm of the unknown's row in an orthonormal basis of the null space (the diagonal -of the null-space projector): how strongly that unknown participates in the -underdetermined directions, independent of the arbitrary basis chosen within the null -space. A `nullity` of `0` means the initialization Jacobian has full column rank at `u0` -and the initialization is locally well-posed. +A `NamedTuple` `(; jacobian, singular_values, rank, nullity, redundancy, underdetermined_unknowns, redundant_equations)`. + + - `nullity = ncols - rank` is the number of underdetermined directions and `redundancy = + nrows - rank` the number of redundant equations. + - `underdetermined_unknowns` is a vector of `unknown => weight` pairs and + `redundant_equations` a vector of `equation => weight` pairs, each sorted by decreasing + `weight ∈ [0, 1]`. The weight is the diagonal of the corresponding null-space projector + (the squared row norm over an orthonormal basis of the right/left null space): how + strongly that unknown/equation participates in the underdetermined/redundant directions, + independent of the arbitrary basis chosen within the null space. + +A `nullity` of `0` means the Jacobian has full column rank at `u0` (no underdetermined +unknowns); a `redundancy` of `0` means full row rank (no redundant equations). !!! note The Jacobian is evaluated at a single point (`u0`), so this reports the *local* rank structure there. A structurally well-posed initialization can still be numerically rank deficient at a particular operating point, and vice versa. """ -function report_initialization_nullspace( - prob; rtol = 1e-8, atol = 0.0, var_threshold = 1e-3, verbose = true) - empty_result = (; jacobian = nothing, singular_values = Float64[], - nullity = 0, variables = Pair[]) +function analyze_initialization_jacobian( + prob; rtol = 1e-8, atol = 0.0, threshold = 1e-3, verbose = true) + empty_result = (; jacobian = nothing, singular_values = Float64[], rank = 0, + nullity = 0, redundancy = 0, underdetermined_unknowns = Pair[], + redundant_equations = Pair[]) iprob = _initialization_problem(prob) if iprob === nothing verbose && @@ -120,44 +131,71 @@ function report_initialization_nullspace( u -> f(u, p) end J = ForwardDiff.jacobian(residual, u0) - n = length(u0) + nrows, ncols = size(J) fact = svd(J; full = true) S = fact.S σmax = isempty(S) ? zero(eltype(S)) : maximum(S) - threshold = max(atol, rtol * σmax) - # A right-singular vector is a null direction when its singular value is below the - # threshold; columns of `V` beyond `length(S)` (present when n > nresid) have an - # implied singular value of zero and are always null directions. - is_null = [(j <= length(S) ? S[j] : zero(σmax)) <= threshold for j in 1:n] - nullity = count(is_null) + tol = max(atol, rtol * σmax) + rank = count(>(tol), S) + # A singular vector is a null direction when its singular value is below the + # tolerance. Right vectors (columns of `V`) beyond `length(S)` (present when + # ncols > nrows) and left vectors (columns of `U`) beyond `length(S)` (present when + # nrows > ncols) have an implied singular value of zero and are always null. + col_isnull = [(j <= length(S) ? S[j] : zero(σmax)) <= tol for j in 1:ncols] + row_isnull = [(i <= length(S) ? S[i] : zero(σmax)) <= tol for i in 1:nrows] + nullity = count(col_isnull) + redundancy = count(row_isnull) + # Participation = diagonal of the null-space projector (squared row norm over an + # orthonormal null-space basis), invariant to the arbitrary basis within the null space. + col_w = nullity == 0 ? zeros(ncols) : + vec(sum(abs2, @view(fact.V[:, col_isnull]); dims = 2)) + row_w = redundancy == 0 ? zeros(nrows) : + vec(sum(abs2, @view(fact.U[:, row_isnull]); dims = 2)) syms = variable_symbols(iprob) - # Participation = diagonal of the null-space projector V_null * V_null', i.e. the - # squared row norm of each unknown over an orthonormal null-space basis. This is - # invariant to the arbitrary choice of basis within the null space. - weights = nullity == 0 ? zeros(n) : vec(sum(abs2, @view(fact.V[:, is_null]); dims = 2)) - variables = sort( - [syms[i] => weights[i] for i in 1:n if weights[i] > var_threshold]; + eqs = _initialization_equations(iprob) + underdetermined_unknowns = sort( + [syms[i] => col_w[i] for i in 1:ncols if col_w[i] > threshold]; + by = last, rev = true) + redundant_equations = sort( + [(eqs === nothing ? i : eqs[i]) => row_w[i] for i in 1:nrows if row_w[i] > threshold]; by = last, rev = true) if verbose - println("Initialization Jacobian null-space report") - println(" residual Jacobian: $(nresid)×$(n), rank ≈ $(n - nullity), nullity ≈ $nullity") + println("Initialization Jacobian rank analysis") + println(" residual Jacobian: $(nrows)×$(ncols), rank ≈ $rank, nullity ≈ $nullity, redundancy ≈ $redundancy") if !isempty(S) k = min(5, length(S)) println(" smallest $k singular value(s): ", round.(S[(end - k + 1):end], sigdigits = 3)) end if nullity == 0 - println(" Full column rank at u0 — initialization is locally well-posed.") + println(" Full column rank at u0 — no underdetermined unknowns.") else - println(" Unknowns spanning the null space (participation weight ∈ [0, 1]):") - for (s, w) in variables + println(" Underdetermined unknowns (participation ∈ [0, 1]):") + for (s, w) in underdetermined_unknowns println(" ", lpad(string(round(w, digits = 4)), 8), " ", s) end - isempty(variables) && - println(" (no individual unknown exceeds var_threshold = $var_threshold)") + end + if redundancy == 0 + println(" Full row rank at u0 — no redundant equations.") + else + println(" Redundant equations (participation ∈ [0, 1]):") + for (e, w) in redundant_equations + println(" ", lpad(string(round(w, digits = 4)), 8), " ", e) + end end end - return (; jacobian = J, singular_values = S, nullity, variables) + return (; jacobian = J, singular_values = S, rank, nullity, redundancy, + underdetermined_unknowns, redundant_equations) +end + +# Symbolic equations of an initialization problem, in residual order, or `nothing` if not +# available (e.g. a problem not built from a `System`). +function _initialization_equations(iprob) + f = iprob.f + hasproperty(f, :sys) || return nothing + sys = f.sys + sys === nothing && return nothing + return equations(sys) end # Return the initialization `NonlinearProblem`/`NonlinearLeastSquaresProblem` carried by diff --git a/test/group_initialization.jl b/test/group_initialization.jl index 4d011487b8..a2b37d43bc 100644 --- a/test/group_initialization.jl +++ b/test/group_initialization.jl @@ -2,6 +2,6 @@ include("shared/mtktestset.jl") @mtktestset("Guess Propagation", "guess_propagation.jl") @safetestset "Hierarchical Initialization Equations" include("hierarchical_initialization_eqs.jl") -@safetestset "Initialization Jacobian null space" include("initialization_nullspace.jl") +@safetestset "Initialization Jacobian analysis" include("initialization_jacobian_analysis.jl") @mtktestset("InitializationSystem Test", "initializationsystem.jl") @mtktestset("Initial Values Test", "initial_values.jl") diff --git a/test/initialization_jacobian_analysis.jl b/test/initialization_jacobian_analysis.jl new file mode 100644 index 0000000000..c9df851bbe --- /dev/null +++ b/test/initialization_jacobian_analysis.jl @@ -0,0 +1,53 @@ +using ModelingToolkit, Test +using ModelingToolkit: t_nounits as t, D_nounits as D + +# Underdetermined initialization: a point constrained to the unit circle, with no +# initial condition pinning where on the circle (or with what velocity) it starts. +# The initialization Jacobian is column-rank deficient, so the analysis must flag a +# nonempty set of underdetermined unknowns. +@variables x(t) y(t) +@named sys = System([D(x) ~ y, 0 ~ x^2 + y^2 - 1], t) +sys = mtkcompile(sys) +prob = ODEProblem(sys, [], (0.0, 1.0); guesses = [x => 0.6, y => 0.8], + warn_initialize_determined = false) +res = analyze_initialization_jacobian(prob; verbose = false) +@test res.nullity ≥ 1 +@test !isempty(res.underdetermined_unknowns) +@test res.jacobian !== nothing +@test all(0 .≤ last.(res.underdetermined_unknowns) .≤ 1 + 1e-8) +@test res.rank + res.nullity == size(res.jacobian, 2) + +# Overdetermined initialization with a redundant equation: two consistent initial +# conditions on a single unknown. The analysis must report a redundant equation. +@variables z(t) +@named sys2 = System([D(z) ~ -z], t; initialization_eqs = [z ~ 1.0, z^2 ~ 1.0]) +sys2 = mtkcompile(sys2) +prob2 = ODEProblem(sys2, [], (0.0, 1.0); guesses = [z => 0.9], + warn_initialize_determined = false) +res2 = analyze_initialization_jacobian(prob2; verbose = false) +@test res2.redundancy ≥ 1 +@test !isempty(res2.redundant_equations) +@test all(0 .≤ last.(res2.redundant_equations) .≤ 1 + 1e-8) + +# Determined initialization: an algebraic unknown uniquely fixed by a pinned state has a +# full-rank initialization Jacobian — no underdetermined unknowns and no redundant +# equations. +@variables u(t) w(t) +@named sys3 = System([D(u) ~ -u, 0 ~ w - exp(u)], t) +sys3 = mtkcompile(sys3) +prob3 = ODEProblem(sys3, [u => 1.0], (0.0, 1.0); guesses = [w => 1.0]) +res3 = analyze_initialization_jacobian(prob3; verbose = false) +@test res3.nullity == 0 +@test res3.redundancy == 0 +@test isempty(res3.underdetermined_unknowns) +@test isempty(res3.redundant_equations) + +# A system fully determined by its initial conditions has no initialization problem to +# analyze; the utility handles that gracefully rather than erroring. +@variables a(t) b(t) +@named sys4 = System([D(a) ~ -a, D(b) ~ -b], t) +sys4 = mtkcompile(sys4) +prob4 = ODEProblem(sys4, [a => 1.0, b => 2.0], (0.0, 1.0)) +res4 = analyze_initialization_jacobian(prob4; verbose = false) +@test res4.nullity == 0 +@test res4.redundancy == 0 diff --git a/test/initialization_nullspace.jl b/test/initialization_nullspace.jl deleted file mode 100644 index 7aa510e365..0000000000 --- a/test/initialization_nullspace.jl +++ /dev/null @@ -1,36 +0,0 @@ -using ModelingToolkit, Test -using ModelingToolkit: t_nounits as t, D_nounits as D - -# Underdetermined initialization: a point constrained to the unit circle, with no -# initial condition pinning where on the circle (or with what velocity) it starts. -# The initialization Jacobian is rank deficient, so the report must flag a nonempty -# null space. -@variables x(t) y(t) -@named sys = System([D(x) ~ y, 0 ~ x^2 + y^2 - 1], t) -sys = mtkcompile(sys) -prob = ODEProblem(sys, [], (0.0, 1.0); guesses = [x => 0.6, y => 0.8], - warn_initialize_determined = false) -res = report_initialization_nullspace(prob; verbose = false) -@test res.nullity ≥ 1 -@test !isempty(res.variables) -@test res.jacobian !== nothing -@test all(0 .≤ last.(res.variables) .≤ 1 + 1e-8) - -# Determined initialization: an algebraic unknown uniquely fixed by a pinned state has -# a full-column-rank initialization Jacobian and no null space. -@variables u(t) w(t) -@named sys2 = System([D(u) ~ -u, 0 ~ w - exp(u)], t) -sys2 = mtkcompile(sys2) -prob2 = ODEProblem(sys2, [u => 1.0], (0.0, 1.0); guesses = [w => 1.0]) -res2 = report_initialization_nullspace(prob2; verbose = false) -@test res2.nullity == 0 -@test isempty(res2.variables) - -# A system fully determined by its initial conditions has no initialization problem to -# analyze; the utility handles that gracefully rather than erroring. -@variables a(t) b(t) -@named sys3 = System([D(a) ~ -a, D(b) ~ -b], t) -sys3 = mtkcompile(sys3) -prob3 = ODEProblem(sys3, [a => 1.0, b => 2.0], (0.0, 1.0)) -res3 = report_initialization_nullspace(prob3; verbose = false) -@test res3.nullity == 0 From fa819f003ccff943353f09f49b3f591d215e56bc Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Fri, 26 Jun 2026 05:21:13 +0000 Subject: [PATCH 3/3] Point over/underdetermined init warnings at `analyze_initialization_jacobian` The overdetermined/underdetermined initialization warnings now suggest calling `analyze_initialization_jacobian(prob)` on the constructed problem to see exactly which unknowns are underdetermined and which equations are redundant, rather than just reporting the equation/unknown counts. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01FSDc4tFYygs26KhqZBACLS --- .../src/problems/initializationproblem.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/lib/ModelingToolkitBase/src/problems/initializationproblem.jl b/lib/ModelingToolkitBase/src/problems/initializationproblem.jl index e24346cfea..e86ade9bf2 100644 --- a/lib/ModelingToolkitBase/src/problems/initializationproblem.jl +++ b/lib/ModelingToolkitBase/src/problems/initializationproblem.jl @@ -131,6 +131,9 @@ function overdetermined_initialization_message(neqs::Integer, nunknown::Integer, Initialization system is overdetermined. $neqs equations for $nunknown unknowns. \ Initialization will default to using least squares. $(extra) + Call `analyze_initialization_jacobian(prob)` on the constructed problem to see which \ + equations are redundant (and which unknowns, if any, remain underdetermined). + To suppress this warning, pass `warn_initialize_determined = false`. To turn this \ warning into an error, pass `fully_determined = true`. """ @@ -141,6 +144,9 @@ function underdetermined_initialization_message(neqs::Integer, nunknown::Integer Initialization system is underdetermined. $neqs equations for $nunknown unknowns. \ Initialization will default to using least squares. $(extra) + Call `analyze_initialization_jacobian(prob)` on the constructed problem to see which \ + unknowns are underdetermined (and which equations, if any, are redundant). + To suppress this warning, pass `warn_initialize_determined = false`. To turn this \ warning into an error, pass `fully_determined = true`. """