fixes in batchlin

Author: baggepinnen

docs/Project.toml

@@ -6,7 +6,8 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739"
 MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
+OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 RobustAndOptimalControl = "21fd56a4-db03-40ee-82ee-a87907bee541"
 SymbolicControlSystems = "886cb795-8fd3-4b11-92f6-8071e46f71c5"

docs/src/batch_linearization.md

@@ -96,7 +96,7 @@ If you plot the Nyquist curve using the `plotly()` backend rather than the defau
 Above, we tuned one controller for each operating point, wouldn't it be nice if we had some features to simulate a [gain-scheduled controller](https://en.wikipedia.org/wiki/Gain_scheduling) that interpolates between the different controllers depending on the operating pont? [`GainScheduledStateSpace`](@ref) is such a thing, we show how to use it below. For fun, we simulate some reference step responses for each individual controller in the array `Cs` and end with simulating the gain-scheduled controller.
 
 ```@example BATCHLIN
-using OrdinaryDiffEq
+using OrdinaryDiffEqNonlinearSolve, OrdinaryDiffEqRosenbrock
 using DataInterpolations # Required to interpolate between the controllers
 @named fb  = Blocks.Add(k2=-1)
 @named ref = Blocks.Square(frequency=1/6, amplitude=0.5, offset=0.5, start_time=1)
@@ -134,7 +134,7 @@ eqs = [
 ]
 @named closed_loop = ODESystem(eqs, t, systems=[duffing, Cgs, fb, ref, F])
 prob = ODEProblem(structural_simplify(closed_loop), [F.xd => 0], (0.0, 8.0))
-sol = solve(prob, Rodas5P(), abstol=1e-8, reltol=1e-8, initializealg=NoInit())
+sol = solve(prob, Rodas5P(), abstol=1e-8, reltol=1e-8, initializealg=SciMLBase.NoInit())
 plot!(sol, idxs=[duffing.y.u, duffing.u.u], l=(2, :red), lab="Gain scheduled")
 plot!(sol, idxs=F.output.u, l=(1, :black, :dash, 0.5), lab="Ref")
 ```

src/ode_system.jl

@@ -523,6 +523,8 @@ Linearize `sys` around the trajectory `sol` at times `t`. Returns a vector of `S
 function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), allow_input_derivatives = false, fuzzer = nothing, verbose = true, kwargs...)
     maximum(t) > maximum(sol.t) && @warn("The maximum time in `t`: $(maximum(t)), is larger than the maximum time in `sol.t`: $(maximum(sol.t)).")
     minimum(t) < minimum(sol.t) && @warn("The minimum time in `t`: $(minimum(t)), is smaller than the minimum time in `sol.t`: $(minimum(sol.t)).")
+
+    # NOTE: we call linearization_funciton twice :( The first call is to get x=unknowns(ssys), the second call provides the operating points.
     lin_fun, ssys = linearization_function(sys, inputs, outputs; kwargs...)
     x = unknowns(ssys)
     defs = ModelingToolkit.defaults(sys)
@@ -536,8 +538,10 @@ function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), allow_inp
         ops = reduce(vcat, opsv)
         t = repeat(t, inner = length(ops) ÷ length(t))
     end
+    lin_fun, ssys = linearization_function(sys, inputs, outputs; op=ops[1], kwargs...)
     lins = map(zip(ops, t)) do (op, t)
         linearize(ssys, lin_fun; op, t, allow_input_derivatives)
+        # linearize(sys, inputs, outputs; op, t, allow_input_derivatives, initialize=false)[1]
     end
     (; linsystems = [ss(l...) for l in lins], ssys, ops)
 end

test/test_batchlin.jl

@@ -1,6 +1,8 @@
 using ControlSystemsMTK, ModelingToolkit, RobustAndOptimalControl, MonteCarloMeasurements
 using OrdinaryDiffEqNonlinearSolve, OrdinaryDiffEqRosenbrock
+using ModelingToolkitStandardLibrary.Blocks
 using ModelingToolkit: getdefault
+using Test
 unsafe_comparisons(true)
 
 # Create a model
@@ -24,7 +26,7 @@ sample_within_bounds((l, u)) = (u - l) * rand() + l
 # Create a vector of operating points
 N = 10
 xs = range(getbounds(x)[1], getbounds(x)[2], length=N)
-ops = Dict.(x .=> xs)
+ops = Dict.(u.u => 0, x .=> xs)
 
 inputs, outputs = [u.u], [y.u]
 Ps, ssys = batch_ss(duffing, inputs, outputs , ops)
@@ -55,11 +57,11 @@ import ModelingToolkitStandardLibrary.Blocks
 
 
 closed_loop_eqs = [
-    connect(ref.output, :r, F.input)
-    connect(F.output, fb.input1)
-    connect(duffing.y, :y, fb.input2)
-    connect(fb.output, C.input)
-    connect(C.output, duffing.u)
+    ModelingToolkit.connect(ref.output, :r, F.input)
+    ModelingToolkit.connect(F.output, fb.input1)
+    ModelingToolkit.connect(duffing.y, :y, fb.input2)
+    ModelingToolkit.connect(fb.output, C.input)
+    ModelingToolkit.connect(C.output, duffing.u)
 ]
 
 @named closed_loop = ODESystem(closed_loop_eqs, t, systems=[duffing, C, fb, ref, F])
@@ -73,6 +75,7 @@ sol = solve(prob, Rodas5P(), abstol=1e-8, reltol=1e-8)
 
 time = 0:0.1:8
 inputs, outputs = [duffing.u.u], [duffing.y.u]
+op = Dict(u.u => 0)
 Ps2, ssys = trajectory_ss(closed_loop, closed_loop.r, closed_loop.y, sol; t=time)
 @test length(Ps2) == length(time)
 # bodeplot(Ps2)