feat: resolve symbolic loop-opening op values from the solution

lib/ModelingToolkitBase/test/analysis_points.jl

@@ -307,6 +307,65 @@
         @test length(mats_lo) == length(ts)
     end
 
+    @testset "Symbolic loop-opening op values resolve from the solution" begin
+        # A symbolic operating-point value (`opened_signal => driving_signal`) is
+        # evaluated from the solution at each time point, so the opened signal is
+        # linearized around the value it has in the loop-closed solution. The controller
+        # state equation is nonlinear in its input, so the linearization genuinely
+        # depends on the value used for the opened signal.
+        @component function NonlinPlant(; name)
+            @variables begin
+                u(t), [input = true]
+                x(t) = 1.0
+                y(t), [output = true]
+            end
+            System([D(x) ~ -x^3 + u, y ~ x], t; name)
+        end
+        @component function NonlinCtrl(; name, k = -2.0)
+            @variables begin
+                u(t), [input = true]
+                x(t) = 0.0
+                y(t), [output = true]
+            end
+            @parameters k = k
+            System([D(x) ~ -(1 + u^2) * x + k * u, y ~ x], t; name)
+        end
+        @named nlP = NonlinPlant()
+        @named nlC = NonlinCtrl()
+        nleqs = [connect(nlP.y, :nly, nlC.u), connect(nlC.y, :nlu, nlP.u)]
+        @named nlsys = System(nleqs, t; systems = [nlP, nlC])
+        ncsys = complete(nlsys)
+        nlprob = ODEProblem(mtkcompile(nlsys), [], (0.0, 5.0))
+        nlsol = solve(nlprob, Rodas5())
+        nlts = collect(0.0:0.5:5.0)
+        sym_op = Dict(ncsys.nlC.u => ncsys.nlP.y)
+
+        # Trajectory path: per-point values match the loop-closed solution.
+        mats_sym, lsys_sym, _ = linearize(
+            nlsys, :nlu, :nly; loop_openings = [:nly],
+            op = ModelingToolkit.LinearizationOpPoint(nlsol, nlts; op = sym_op)
+        )
+        iC = findfirst(v -> occursin("nlC₊x", string(v)), unknowns(lsys_sym))
+        got = [m.A[iC, iC] for m in mats_sym]
+        truth = [-(1 + nlsol(ti, idxs = ncsys.nlP.y)^2) for ti in nlts]
+        @test got ≈ truth rtol = 1.0e-6
+        @test !allequal(round.(got; digits = 6))
+
+        # Scalar path agrees.
+        mats_1, _, _ = linearize(
+            nlsys, :nlu, :nly; loop_openings = [:nly],
+            op = ModelingToolkit.LinearizationOpPoint(nlsol, 2.5; op = sym_op)
+        )
+        @test mats_1.A[iC, iC] ≈ -(1 + nlsol(2.5, idxs = ncsys.nlP.y)^2) rtol = 1.0e-6
+
+        # `missing` passes through the trajectory path without error.
+        mats_miss, _, _ = linearize(
+            nlsys, :nlu, :nly; loop_openings = [:nly],
+            op = ModelingToolkit.LinearizationOpPoint(nlsol, nlts; op = Dict(ncsys.nlC.u => missing))
+        )
+        @test length(mats_miss) == length(nlts)
+    end
+
     @testset "Complicated model" begin
         # Parameters
         m1 = 1

src/linearization.jl

@@ -13,8 +13,11 @@ extras.
 The `op` keyword argument provides additional operating-point values that are merged into
 the solution-derived operating point at every time point (taking precedence). This is how
 values for variables that are not present in `sol` are supplied — in particular the
-parameters created by `loop_openings`, e.g.
-`LinearizationOpPoint(sol, t; op = Dict(opened_signal => 0))`.
+parameters created by `loop_openings`. Values may be numbers, e.g.
+`LinearizationOpPoint(sol, t; op = Dict(opened_signal => 0))`, or symbolic expressions,
+which are evaluated from the solution at each time point. Mapping an opened signal to the
+expression that drives it, e.g. `op = Dict(opened_signal => driving_signal)`, thus
+linearizes around the value the opened signal has in the loop-closed solution.
 
 # Fields
 
@@ -45,6 +48,17 @@ end
 # to transfer values from the solution to the `LinearizationProblem` without having to
 # materialize and manipulate a dictionary alongisde a corresponding list of setter functions.
 
+# Values in `LinearizationOpPoint.op` may be symbolic expressions (e.g.
+# `opened_signal => driving_signal`), which are evaluated from the solution at each time
+# point rather than passed through verbatim. Resolving them here means the downstream
+# machinery (problem construction, per-point setters) only ever sees numbers: symbolic
+# values cannot be evaluated eagerly at problem construction, and resolving them against
+# the `LinearizationProblem` inside the per-point loop would read state values that are
+# only refreshed by initialization inside `solve`, i.e. those of the previous time point.
+function _is_symbolic_op_value(v)
+    return !SU.isconst(unwrap(v))
+end
+
 function _build_op_from_solution(op::LinearizationOpPoint)
     sol_sys = MTKBase.indp_to_system(op.sol)
     eqs = equations(sol_sys)
@@ -59,7 +73,7 @@ function _build_op_from_solution(op::LinearizationOpPoint)
         result[p] = getp(op.sol, p)(op.sol)
     end
     for (k, v) in op.op
-        result[unwrap(k)] = v
+        result[unwrap(k)] = _is_symbolic_op_value(v) ? op.sol(op.t; idxs = v) : v
     end
     return result
 end
@@ -76,15 +90,28 @@ function _build_op_from_solution(op::LinearizationOpPoint{S, <:AbstractVector})
     for p in parameters(sol_sys)
         param_vals[p] = getp(op.sol, p)(op.sol)
     end
+    # Split the extra op into constant entries and symbolic entries; the latter are
+    # evaluated from the solution, once per entry for all time points.
+    extra_const = Dict{SymbolicT, SymbolicT}()
+    extra_resolved = Pair{SymbolicT, eltype(eltype(sol.u))}[]
+    for (k, v) in op.op
+        if _is_symbolic_op_value(v)
+            push!(extra_resolved, unwrap(k) => op.sol(op.t; idxs = v))
+        else
+            extra_const[unwrap(k)] = v
+        end
+    end
     # Interpolate once per time point to build the per-point operating-point dict.
-    extra_op = Dict{SymbolicT, SymbolicT}(unwrap(k) => v for (k, v) in op.op)
-    return map(op.t) do ti
+    return map(enumerate(op.t)) do (i, ti)
         u = op.sol(ti)
         result = copy(param_vals)
-        for i in diff_idxs
-            result[sts[i]] = u[i]
+        for j in diff_idxs
+            result[sts[j]] = u[j]
+        end
+        merge!(result, extra_const)
+        for (k, vals) in extra_resolved
+            result[k] = vals[i]
         end
-        merge!(result, extra_op)
         result
     end
 end
@@ -1026,7 +1053,10 @@ function _check_loop_opening_op(loop_opening_params, op)
         parameters whose operating-point values are not implied by the rest of the \
         system, so they must be provided explicitly in `op` (e.g. set to zero). When \
         linearizing along a trajectory with `LinearizationOpPoint`, pass them via its \
-        `op` keyword argument: `LinearizationOpPoint(sol, t; op = Dict(signal => value))`.
+        `op` keyword argument: `LinearizationOpPoint(sol, t; op = Dict(signal => value))`, \
+        where the value may also be a symbolic expression evaluated from the solution at \
+        each time point, e.g. `Dict(opened_signal => driving_signal)` to linearize around \
+        the value the opened signal has in the loop-closed solution.
         """
     )
 end
@@ -1058,6 +1088,7 @@ function __linearize_multiple_op_barrier(ssys, lin_fun; ops, ts, allow_input_der
         for (setter, k) in zip(setters, op_keys)
             v = get(op, k, COMMON_NOTHING)
             isequal(v, COMMON_NOTHING) && continue
+            v === COMMON_MISSING && continue
             setter(prob, _resolve_op_value(prob, v))
         end
         prob.t = t