Feature Request: Selective Plotting of State Variables vs Derivatives

[ChrisRackauckas-Claude]

Feature Request: Selective Plotting of State Variables vs Derivatives

Summary

Add support for plotting only state variables (without their derivatives) when calling plot(sol) on ModelingToolkit solutions. This would provide a cleaner visualization for complex models where derivative plots can be overwhelming.

Motivation

Currently, when using plot(sol) on a ModelingToolkit solution, all variables including derivatives are plotted. For models with many state variables, this can result in cluttered plots with 2x the number of series (state + derivative for each variable).

Users often want to quickly visualize just the "primary" state variables without derivatives. The current workaround requires manually specifying idxs with individual variables, which defeats the convenience of plot(myanalysis()).

Proposed Solution

Extend SymbolicIndexingInterface (SII) with functionality to distinguish between state variables and their derivatives, allowing syntax like:

# Option 1: Using SII interface
plot(sol, idxs=SII.STATE_VARIABLES)  # Only state variables
plot(sol, idxs=SII.ALL_DERIVATIVES)  # Only derivatives

# Option 2: Plotting keyword argument  
plot(sol, plot_derivs=false)  # Exclude derivatives

This would integrate with the existing plotting infrastructure while providing users more control over what gets visualized.

Minimal Working Example

using ModelingToolkit, OrdinaryDiffEqDefault, Plots
using ModelingToolkit: t_nounits as t, D_nounits as D

# Simple cartesian pendulum example
@parameters g = 9.81 L = 1.0  # gravity, length
@variables x(t) y(t)  # cartesian coordinates

# Pendulum equations with constraint x² + y² = L²
eqs = [
    D(D(x)) ~ -x/L^2 * g * y/L,
    D(D(y)) ~ -y/L^2 * g * y/L - g,
    0 ~ x^2 + y^2 - L^2  # constraint
]

@mtkbuild sys = ODESystem(eqs, t)

# Initial conditions: pendulum starts at 30° from vertical
θ₀ = π/6
u0 = [
    x => L * sin(θ₀),
    y => -L * cos(θ₀),
    D(x) => 0.0,
    D(y) => 0.0
]

# Solve the system
prob = ODEProblem(sys, u0, (0.0, 3.0), [])
sol = solve(prob)

# Current behavior: plots ALL variables including derivatives
plot(sol)  # Shows x(t), y(t), D(x)(t), D(y)(t), etc.

# Desired behavior:
plot(sol, idxs=SII.STATE_VARIABLES)  # Only x(t), y(t)
# or
plot(sol, plot_derivs=false)  # Only x(t), y(t)

Current Workaround

Users must manually specify which variables to plot:

state_variables = [x, y]  # Manually list state variables
plot(sol, idxs=state_variables)

This approach doesn't scale well for models with many state variables and requires users to know which variables are "primary" vs derivatives.

Related Issues

This relates to issue #3555 about plotting unknowns and observables, but focuses specifically on the state/derivative distinction rather than the unknown/observable distinction.

Implementation Notes

This would likely involve:

1. Extending SymbolicIndexingInterface to categorize variables by type (state vs derivative)
2. Adding plotting keywords or SII constants for selective plotting
3. Updating the default plotting behavior to optionally filter variable types

The implementation should maintain backward compatibility while providing users with more intuitive control over visualization complexity.

[ChrisRackauckas-Claude]

Implementation Details: Using System Structure to Identify Variable Types

After looking at the ModelingToolkit codebase, I found that the initialization system already has the infrastructure to distinguish between differential and algebraic variables. Here's how this could be implemented:

Current Infrastructure

ModelingToolkit has several utility functions in src/systems/systemstructure.jl and src/utils.jl for categorizing variables:

# From src/utils.jl
isdifferential(expr) = isoperator(expr, Differential)
isdiffeq(eq) = isdifferential(eq.lhs) || isoperator(eq.lhs, Shift)

# From src/systems/systemstructure.jl
isalgvar(s::SystemStructure, i)    # Check if variable i is algebraic
isdiffvar(s::SystemStructure, i)   # Check if variable i is differential 
isdervar(s::SystemStructure, i)    # Check if variable i is a derivative

# Range iterators for different variable types
algvars_range(s::SystemStructure)   # Iterator over algebraic variables
diffvars_range(s::SystemStructure)  # Iterator over differential variables  
dervars_range(s::SystemStructure)   # Iterator over derivative variables

Usage in Initialization Code

The initialization system in src/systems/nonlinear/initializesystem.jl already uses these patterns:

# Line 75: Identify which equations are differential
idxs_diff = isdiffeq.(eqs)

# Line 95-96: Separate algebraic and differential equations  
idxs_alge = .\!idxs_diff
append\!(eqs_ics, eqs[idxs_alge]) # start with algebraic equations

Proposed Implementation

For the plotting feature, we could extend SymbolicIndexingInterface with constants that use this existing infrastructure:

# In SymbolicIndexingInterface
struct StateVariables end
struct DerivativeVariables end  
struct AlgebraicVariables end

const STATE_VARIABLES = StateVariables()
const DERIVATIVE_VARIABLES = DerivativeVariables() 
const ALGEBRAIC_VARIABLES = AlgebraicVariables()

# Implementation would filter variables using the existing functions
function Base.getindex(sol::ODESolution, ::StateVariables)
    sys = sol.prob.f.sys
    # Use diffvars_range(structure) and algvars_range(structure) 
    # to get indices of non-derivative variables
end

This approach leverages the existing, well-tested infrastructure for variable categorization that's already used throughout ModelingToolkit's structural analysis and initialization systems.

[ChrisRackauckas]

@AayushSabharwal thoughts on syntax and structure here?

[AayushSabharwal]

Yeah, we could add this to SII. It already has solvedvariables and allvariables. We'll probably need to special case the plot recipe to handle this too, though.