RFC: Debugging utilities for initialization

[AayushSabharwal]

MTK's initialization can be useful, but a simple warning telling the user their initialization system is over/under-determined is not very useful for debugging. This issue is to discuss tools MTK should provide to aid users in debugging issues with initialization.

Better error messages for missing guesses

If the user doesn't provide required guesses, we throw a MissingVariablesError which is identical to what we would throw had they not specified an initial condition for the variable at all. This doesn't convey the difference between a missing initial condition and a missing guess. The easiest solution here would be a try-catch block around the call to InitializationProblem in problem_utils.jl which would intercept specifically MissingVariablesError or MissingParametersError and throw a MissingGuessesError instead. The overhead of a try-catch is negligible compared to the time we take for codegen, or structural_simplify of the initialization system.

get_initialization_problem, get_initialization_system

This is a simple function which would take a problem or solution and return the initialization problem/system stored inside that is used during solve.

Logging when generating initialization system

Basically an optional @info every time generate_initializesystem does something interesting (mainly add an equation). This could be a verbose keyword, but generate_initializesystem can add a lot of equations and the user may thus want to only log e.g. equations added from initial conditions, or ones added for parameter initialization. Being able to add some sort of tags to log messages and allowing users to filter them easily would be helpful. Alternatively, we could use ScopedValues (ScopedValues.jl for 1.10 compatibility) to allow the user to toggle flags that determine which messages are logged. This is less extensible than custom filter functions for logs, but easier to setup and may be sufficient.

solve_initialization

A function which takes a problem, a nonlinear solver and tolerance (same keyword arguments as SciMLBase.get_initial_values) and returns the result of solving the initialization problem (a NonlinearSolution). This allows users to see which equations have a high residual and are failing initialization.

substitute_defaults

A function which takes a system and returns a new one with all defaults substituted into the equations/observed. Could also have a filter function to subset the substituted defaults.

substitute_parameters_from_problem

A function which takes a system and corresponding problem, and substitutes parameter values from the problem into the system's equations/observed. Useful to see "what's the relation between the variables I'm solving for" without having to parse parameter symbols. Could also have a filter function to subset the substituted parameters.

[AayushSabharwal]

I also think that we should explicitly document at least some of these as utility functions that are public API but not stable, designed to be used for debugging models but not as part of a package/finished product. That way we can change/add/remove them at will.

[SebastianM-C]

How can one identify which equations violate the initialization constraints? It would be possible to tell if a default or an initialization equation is the reason for the failure or we just have the equation(s)?

[AayushSabharwal]

So the idea is that with the logging it should be possibly to identify where each equation in the initialization system came from (with some effort, since structural simplification likes to shuffle terms around). Now suppose you have the initialization system isys from get_initialization_system and the solution nlsol from solve_initialization. By observing the indexes of nlsol.resid, it is possible to identify which of equations(isys) are not satisfied. Whether it is not satisfied due to the equation being incorrect, or because a default is incorrect, or otherwise is very model specific and I'm not aware of a general method of identifying such an issue.

From the logs, you know where that equation came from (a symbolic default, initialization_eqs, an observed equation, etc.). Presumably observed equations are good, so if it comes from one of the other two sources that could be something to look into. Alternatively, suppose the index of the unsatisfied equation is i. Looking at full_equations(isys)[i] can be helpful. If there are only parameters involved, then it's either an incorrect parameter value or a bad Initial value (Initial(x) is also a parameter). The involved Initial variables are visible in the equation, so their values can be checked. substitute_parameters_from_problem with a filter can be used to substitute parameters whose value is known to be correct, leaving behind the remaining parameters whose value needs to be checked.

The approach above is what I would do. With these debugging utilities, I also plan to add a doc page highlighting such a workflow.