diff --git a/lib/ModelingToolkitBase/test/analysis_points.jl b/lib/ModelingToolkitBase/test/analysis_points.jl index 0722b710d3..40fa576832 100644 --- a/lib/ModelingToolkitBase/test/analysis_points.jl +++ b/lib/ModelingToolkitBase/test/analysis_points.jl @@ -307,6 +307,65 @@ if @isdefined(ModelingToolkit) @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 diff --git a/src/linearization.jl b/src/linearization.jl index 3e7e3fd7d5..9df9e23f46 100644 --- a/src/linearization.jl +++ b/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