diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 64be270..027b603 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -1,40 +1,45 @@ name: CI on: - - push - - pull_request + push: + branches: + - main + tags: ['*'] + pull_request: + workflow_dispatch: +concurrency: + # Skip intermediate builds: always. + # Cancel intermediate builds: only if it is a pull request build. + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} jobs: test: - name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} + name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} runs-on: ${{ matrix.os }} + timeout-minutes: 60 + permissions: # needed to allow julia-actions/cache to proactively delete old caches that it has created + actions: write + contents: read strategy: fail-fast: false matrix: version: - - '1' + - '1.11' os: - ubuntu-latest arch: - x64 steps: - - uses: actions/checkout@v2 - - uses: julia-actions/setup-julia@v1 + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - - uses: actions/cache@v4 - env: - cache-name: cache-artifacts - with: - path: ~/.julia/artifacts - key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }} - restore-keys: | - ${{ runner.os }}-test-${{ env.cache-name }}- - ${{ runner.os }}-test- - ${{ runner.os }}- + - uses: julia-actions/cache@v2 - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 - uses: codecov/codecov-action@v4 with: - file: lcov.info + files: lcov.info token: ${{ secrets.CODECOV_TOKEN }} + fail_ci_if_error: false diff --git a/Project.toml b/Project.toml index cdb5d3a..e66b68c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,8 @@ name = "ControlSystemsMTK" uuid = "687d7614-c7e5-45fc-bfc3-9ee385575c88" authors = ["Fredrik Bagge Carlson"] -version = "2.3.3" +version = "2.4.0" + [deps] ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e" @@ -17,7 +18,7 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" [compat] ControlSystemsBase = "1.0.1" DataInterpolations = "3, 4, 5, 6, 7" -ModelingToolkit = "9.61" +ModelingToolkit = "10.3" ModelingToolkitStandardLibrary = "2" MonteCarloMeasurements = "1.1" RobustAndOptimalControl = "0.4.14" diff --git a/docs/src/batch_linearization.md b/docs/src/batch_linearization.md index 30feb3c..0cded72 100644 --- a/docs/src/batch_linearization.md +++ b/docs/src/batch_linearization.md @@ -25,13 +25,13 @@ eqs = [D(x) ~ v y.u ~ x] -@named duffing = ODESystem(eqs, t, systems=[y, u]) +@named duffing = System(eqs, t, systems=[y, u]) ``` ## Batch linearization To perform batch linearization, we create a vector of operating points, and then linearize the model around each of these points. The function [`batch_ss`](@ref) does this for us, and returns a vector of `StateSpace` models, one for each operating point. An operating point is a `Dict` that maps variables in the MTK model to numerical values. In the example below, we simply sample the variables uniformly within their bounds specified when we created the variables (normally, we might want to linearize on stationary points) ```@example BATCHLIN -N = 16 # Number of samples +N = 5 # Number of samples xs = range(getbounds(x)[1], getbounds(x)[2], length=N) ops = Dict.(u.u => 0, x .=> xs) ``` @@ -56,11 +56,15 @@ P = RobustAndOptimalControl.ss2particles(Ps) # convert to a single StateSpace sy notice how some coefficients are plotted like uncertain numbers `-13.8 ± 4.3`. We can plot such models as well: ```@example BATCHLIN -bodeplot(P, w, legend=:bottomright) # Should look similar to the one above +bodeplot(P, w, legend=:bottomright, adaptive=false) # Should look similar to the one above ``` ## Controller tuning Let's also do some controller tuning for the linearized models above. The function `batch_tune` is not really required here, but it shows how we might go about building more sophisticated tools for batch tuning. In this example, we will tune a PID controller using the function [`loopshapingPID`](@ref). Note, this procedure is not limited to tuning a gain-scheduled PID controller, it should work for gain-scheduling of any LTI controller. + +!!! "note" Interpolating between controllers + There are multiple ways in which one could interpolate between different controllers. The two most common approaches are to interpolate their outputs, and interpolating their coefficients. When interpolating the coefficients, it is important to ensure that all controllers have the same structure for the interpolation to be meaningful. One may for example interpolate between PID coefficients, or between the coefficients of a state-space model. When interpolating state-space matrices, the systems must all share the same basis, i.e., the state variables must all have the same meaning among the interpolated systems. When converting a transfer function to state-space form, a numerical balancing is performed, this alters the meaning of the state variables and may introduce artificial dynamics to the interpolated system. We thus pass `balance=false` to the function `ss` to avoid this, or pick a form explicitly, e.g., `modal_form`. + ```@example BATCHLIN function batch_tune(f, Ps) f.(Ps) @@ -68,7 +72,8 @@ end Cs = batch_tune(Ps) do P C, kp, ki, kd, fig, CF = loopshapingPID(P, 7; Mt=1.2, Tf = 1/100) - ss(CF) + ss(CF, balance=false) + modal_form(ss(CF, balance=true))[1] end P = RobustAndOptimalControl.ss2particles(Ps) @@ -103,7 +108,7 @@ using DataInterpolations # Required to interpolate between the controllers connect = ModelingToolkit.connect closed_loop_eqs = [ - connect(ref.output, F.input) + connect(ref.output, :r0, F.input) connect(F.output, :r, fb.input1) # Add an analysis point :r connect(duffing.y, :y, fb.input2) # Add an analysis point :y ] @@ -111,14 +116,14 @@ plot(layout=2) # Simulate each individual controller for C in Cs - @named Ci = ODESystem(C) + @named Ci = System(C) eqs = [ closed_loop_eqs connect(fb.output, Ci.input) connect(Ci.output, duffing.u) ] - @named closed_loop = ODESystem(eqs, t, systems=[duffing, Ci, fb, ref, F]) - prob = ODEProblem(structural_simplify(closed_loop), [F.xd => 0], (0.0, 8.0)) + @named closed_loop = System(eqs, t, systems=[duffing, Ci, fb, ref, F]) + prob = ODEProblem(structural_simplify(closed_loop), [F.x => 0, F.xd => 0], (0.0, 8.0)) sol = solve(prob, Rodas5P(), abstol=1e-8, reltol=1e-8) plot!(sol, idxs=[duffing.y.u, duffing.u.u], layout=2, lab="") end @@ -128,25 +133,25 @@ end eqs = [ closed_loop_eqs connect(fb.output, Cgs.input) - connect(Cgs.output, duffing.u) - connect(duffing.y, Cgs.scheduling_input) # Don't forget to connect the scheduling variable! + connect(Cgs.output, :u, duffing.u) + connect(duffing.y, :v, Cgs.scheduling_input) # Don't forget to connect the scheduling variable! ] -@named closed_loop = ODESystem(eqs, t, systems=[duffing, Cgs, fb, ref, F]) +@named closed_loop = System(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=SciMLBase.NoInit()) +sol = solve(prob, Rodas5P(), abstol=1e-8, reltol=1e-8, initializealg=SciMLBase.NoInit(), dtmax=0.01) 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") ``` -If everything worked as expected, the gain-scheduled controller should perform better than each of the included controllers individually. +If everything worked as expected, the gain-scheduled controller should perform reasonably well across the entire scheduling range. ## C-Code generation -We can generate C-code to interpolate our controller using the function [`SymbolicControlSystems.print_c_array`](@ref) from [SymbolicControlSystems.jl](https://github.com/JuliaControl/SymbolicControlSystems.jl). If the controller is a standard [`ControlSystemsBase.StateSpace`](@ref) object, a function that filters the input through the controller can be generated by calling [`SymbolicControlSystems.ccode`](@ref). But if the controller is a vector of controllers representing a gain-scheduled controller, a function that creates the interpolated dynamics is written. In the code below, we shorten the vector of controllers to make the generated C-code easier to read by passing `Cs[1:7:end]` and `xs[1:7:end]` +We can generate C-code to interpolate our controller using the function [`SymbolicControlSystems.print_c_array`](@ref) from [SymbolicControlSystems.jl](https://github.com/JuliaControl/SymbolicControlSystems.jl). If the controller is a standard [`ControlSystemsBase.StateSpace`](@ref) object, a function that filters the input through the controller can be generated by calling [`SymbolicControlSystems.ccode`](@ref). But if the controller is a vector of controllers representing a gain-scheduled controller, a function that creates the interpolated dynamics is written. In the code below, we shorten the vector of controllers to make the generated C-code easier to read by passing `Cs[1:3:end]` and `xs[1:3:end]` ```@example BATCHLIN using SymbolicControlSystems, ControlSystemsBase Cs_disc = c2d.(Cs, 0.05, :tustin) # Discretize the controller before generating code -code = SymbolicControlSystems.print_c_array(stdout, Cs_disc[1:7:end], xs[1:7:end], "Cs") +code = SymbolicControlSystems.print_c_array(stdout, Cs_disc[1:3:end], xs[1:3:end], "Cs") ``` The generated code starts by defining the interpolation vector `xs`, this variable is called `Cs_interp_vect` in the generated code. The code then defines all the ``A`` matrices as a 3-dimensional array, followed by a function that performs the interpolation `interpolate_Cs_A`. This function takes the output array as the first argument, a pointer to the 3D array with interpolation matrices, the interpolation vector as well as the interpolation variable `t`, in this document called ``v``. The same code is then repeated for the matrices ``B,C,D`` as well if they require interpolation (if they are all the same, no interpolation code is written). @@ -154,13 +159,85 @@ The generated code starts by defining the interpolation vector `xs`, this variab We can linearize around a trajectory obtained from `solve` using the function [`trajectory_ss`](@ref). We provide it with a vector of time points along the trajectory at which to linearize, and in this case we specify the inputs and outputs to linearize between as analysis points `r` and `y`. ```@example BATCHLIN timepoints = 0:0.01:8 -Ps2, ssys = trajectory_ss(closed_loop, closed_loop.r, closed_loop.y, sol; t=timepoints, initialize=true, verbose=true) +Ps2, ssys, ops2, resolved_ops = trajectory_ss(closed_loop, closed_loop.r, closed_loop.y, sol; t=timepoints, verbose=true); bodeplot(Ps2, w, legend=false) ``` +Not how the closed-loop system changes very little along the trajectory, this is a good indication that the gain-scheduled controller is able to make the system appear linear. Internally, [`trajectory_ss`](@ref) works very much the same as [`batch_ss`](@ref), but constructs operating points automatically along the trajectory. This requires that the solution contains the states of the simplified system, accessible through the `idxs` argument like `sol(t, idxs=x)`. By linearizing the same system as we simulated, we ensure that this condition holds, doing so requires that we specify the inputs and outputs as analysis points rather than as variables. +We can replicate the figure above by linearizing the plant and the controller individually, by providing the `loop_openings` argument. When linearizing the plant, we disconnect the controller input by passing `loop_openings=[closed_loop.u]`, and when linearizing the controller, we have various options for disconnecting the the plant: +- Break the connection from plant output to controller input by passing `loop_openings=[closed_loop.y]` +- Break the connection between the controller and the plant input by passing `loop_openings=[closed_loop.u]` +- Break the connection `y` as well as the scheduling variable `v` (which is another form of feedback) by passing `loop_openings=[closed_loop.y, closed_loop.v]` + +We will explore these options below, starting with the first option, breaking the connection `y`: +```@example BATCHLIN +kwargs = (; adaptive=false, legend=false) +plants, _ = trajectory_ss(closed_loop, closed_loop.u, closed_loop.y, sol; t=timepoints, verbose=true, loop_openings=[closed_loop.u]); +controllersy, ssy, ops3, resolved_ops3 = trajectory_ss(closed_loop, closed_loop.r, closed_loop.u, sol; t=timepoints, verbose=true, loop_openings=[closed_loop.y]); + +closed_loopsy = feedback.(plants .* controllersy) +bodeplot(closed_loopsy, w; title="Loop open at y", kwargs...) +``` +When we open the loop at `u`, we get a slightly different result: +```@example BATCHLIN +controllersu, ssu = trajectory_ss(closed_loop, closed_loop.r, closed_loop.u, sol; t=timepoints, verbose=true, loop_openings=[closed_loop.u]); + +closed_loopsu = feedback.(plants .* controllersu) +bodeplot(closed_loopsu, w; title="Loop open at u", kwargs...) +``` +In this case, all static gains are 1. A similar result is obtained by explicitly breaking the scheduling feedback `v` in addition to an opening of `y`: +```@example BATCHLIN +controllersv, ssv = trajectory_ss(closed_loop, closed_loop.r, closed_loop.u, sol; t=timepoints, verbose=true, loop_openings=[closed_loop.y, closed_loop.v]); + +closed_loopsv = feedback.(plants .* controllersv) +bodeplot(closed_loopsv, w; title="Loop open at v and y", kwargs...) +``` + +We have thus far treated the controller as a SISO system, but we could also view it as a system with two inputs, the measurement feedback and the scheduling feedback. For completeness, we show below how to derive the corresponding MIMO systems + +```@example BATCHLIN +plants_mimo, _ = trajectory_ss(closed_loop, closed_loop.u, [closed_loop.y, closed_loop.v], sol; t=timepoints, verbose=true, loop_openings=[closed_loop.u]); +controllers_mimo, ssm = trajectory_ss(closed_loop, [closed_loop.r, closed_loop.v], closed_loop.u, sol; t=timepoints, verbose=true, loop_openings=[closed_loop.u]); + +closed_loops_mimo = feedback.(controllers_mimo .* plants_mimo) # Look at complementary sensitivity function in the input, since this is a SISO system +bodeplot(closed_loops_mimo, w; title="Loop open at MIMO", kwargs...) +``` + + + +Why are the results different depending on where we open the loop? We can understand the difference by comparing the Bode plots of the controllers. + +```@example BATCHLIN +plot( + bodeplot(Cs, w, legend=false, plotphase=false, title="Designed controllers"), + bodeplot(controllersy, w, legend=false, plotphase=false, title="Loop open at y"), + bodeplot(controllersu, w, legend=false, plotphase=false, title="Loop open at u"), + bodeplot(controllersv, w, legend=false, plotphase=false, title="Loop open at v and y"), +) +``` +if we open at both `y` and `v` or we open at `u`, we get controllers for the different values of the scheduling variable, and the corresponding measurement feedback (which is the same as the scheduling variable in this case). +```@example BATCHLIN +using Test +@test all(sminreal.(controllersv) .== sminreal.(controllersu)) +``` + +However, if we only open at `y` we get controller linearizations that _still contain the closed loop through the scheduling connection_ `v`. We can verify this by looking at what variables are present in the input-output map +```@example BATCHLIN +sminreal(controllersy[end]).x +``` +notice how the state of the plant is included in the controller, this is an indication that we didn't fully isolate the controller during the linearizaiton. If we do the same thing for the controller with the loop opened at `u`, we see that the state of the plant is not included in the controller: +```@example BATCHLIN +sminreal(controllersu[end]).x +``` +The call to `sminreal` is important here, it removes the states that are not needed to represent the input-output map of the system. The state of the full model, including the plant state, is present in the linearization before this call. + + + +The easiest way to ensure that the controller is properly disconnected from the plant while taking into account the different scheduling along the trajectory is thus to break at `u`. + ## Summary We have seen how to - Perform linearization of a nonlinear ModelingToolkit model in multiple different operating points @@ -173,6 +250,5 @@ Batch linearization in multiple different operating points is an intuitive way t ```@example BATCHLIN -using Test @test sol(6.99, idxs=closed_loop.duffing.y.u) ≈ 0.0 atol=0.01 ``` \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md index 3dfd0ad..c1285e3 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -21,7 +21,7 @@ pkg> add ControlSystemsMTK ## From ControlSystems to ModelingToolkit -Simply calling `ODESystem(sys)` converts a `StateSpace` object from ControlSystems into the corresponding [`ModelingToolkitStandardLibrary.Blocks.StateSpace`](http://mtkstdlib.sciml.ai/dev/API/blocks/#ModelingToolkitStandardLibrary.Blocks.StateSpace). If `sys` is a [named statespace object](https://juliacontrol.github.io/RobustAndOptimalControl.jl/dev/#Named-systems), the names of inputs and outputs will be retained in the `ODESystem` as connectors, that is, if `my_input` is an input variable in the named statespace object, `my_input` will be a connector of type `RealInput` in the resulting ODESystem. Names of state variables are currently ignored. +Simply calling `System(sys)` converts a `StateSpace` object from ControlSystems into the corresponding [`ModelingToolkitStandardLibrary.Blocks.StateSpace`](http://mtkstdlib.sciml.ai/dev/API/blocks/#ModelingToolkitStandardLibrary.Blocks.StateSpace). If `sys` is a [named statespace object](https://juliacontrol.github.io/RobustAndOptimalControl.jl/dev/#Named-systems), the names of inputs and outputs will be retained in the `System` as connectors, that is, if `my_input` is an input variable in the named statespace object, `my_input` will be a connector of type `RealInput` in the resulting System. Names of state variables are currently ignored. ### Example: @@ -41,7 +41,7 @@ D = Continuous-time state-space model -julia> @named P = ODESystem(P0) +julia> @named P = System(P0) Model P with 2 equations States (3): x[1](t) [defaults to 0.0] @@ -65,15 +65,15 @@ using ControlSystemsMTK, ControlSystemsBase, ModelingToolkit, RobustAndOptimalCo P = named_ss(DemoSystems.double_mass_model(outputs = [1,3]), u=:torque, y=[:motor_angle, :load_angle]) ``` -When we convert this system to an ODESystem, we get a system with connectors `P.torque` and `P.motor_angle`, in addition to the standard connectors `P.input` and `P.output`: +When we convert this system to a `System`, we get a system with connectors `P.torque` and `P.motor_angle`, in addition to the standard connectors `P.input` and `P.output`: ```@example CONNECT -@named P_ode = ODESystem(P) +@named P_ode = System(P) ``` Here, `P.torque` is equal to `P.input`, so you may choose to connect to either of them. However, since the output is multivariable, the connector `P.output` represents both outputs, while `P.motor_angle` and `P.load_angle` represent the individual scalar outputs. ## From ModelingToolkit to ControlSystems -An `ODESystem` can be converted to a named statespace object from [RobustAndOptimalControl.jl](https://github.com/JuliaControl/RobustAndOptimalControl.jl) by calling [`named_ss`](@ref) +A `System` can be converted to a named statespace object from [RobustAndOptimalControl.jl](https://github.com/JuliaControl/RobustAndOptimalControl.jl) by calling [`named_ss`](@ref) ```julia named_ss(ode_sys, inputs, outputs; op) @@ -167,9 +167,9 @@ function SystemModel(u=nothing; name=:model) ] if u !== nothing push!(eqs, connect(u.output, :u, torque.tau)) - return @named model = ODESystem(eqs, t; systems = [sens, torque, inertia1, inertia2, spring, damper, u]) + return @named model = System(eqs, t; systems = [sens, torque, inertia1, inertia2, spring, damper, u]) end - ODESystem(eqs, t; systems = [sens, torque, inertia1, inertia2, spring, damper], name) + System(eqs, t; systems = [sens, torque, inertia1, inertia2, spring, damper], name) end model = SystemModel() |> complete @@ -194,7 +194,8 @@ That's pretty cool, but even nicer is to generate some code for this symbolic sy ```@example LINEAIZE_SYMBOLIC defs = ModelingToolkit.defaults(simplified_sys) defs = merge(Dict(unknowns(model) .=> 0), defs) -x, pars = ModelingToolkit.get_u0_p(simplified_sys, defs, defs) # Extract the default state and parameter values +x = ModelingToolkit.get_u0(simplified_sys, defs) # Extract the default state and parameter values +pars = ModelingToolkit.get_p(simplified_sys, defs, split=false) fun = Symbolics.build_function(symbolic_sys, unknowns(simplified_sys), ModelingToolkit.parameters(simplified_sys); expression = Val{false}, # Generate a compiled function rather than a Julia expression diff --git a/src/ControlSystemsMTK.jl b/src/ControlSystemsMTK.jl index 418cc0f..1e987e9 100644 --- a/src/ControlSystemsMTK.jl +++ b/src/ControlSystemsMTK.jl @@ -3,7 +3,7 @@ using RobustAndOptimalControl: NamedStateSpace #= Ideas: All connections handled by ModelingToolkit. Names: -- handled either by named system, or directly in constructor to ODESystem. +- handled either by named system, or directly in constructor to System. Functions: - Give me linear system from [u1, u3] to [qm, a] If the linearization of a full system produces a named system, one could implement getindex for vectors of names and obtain the desired transfer functions. @@ -20,7 +20,7 @@ using ModelingToolkit, ControlSystemsBase using ControlSystemsBase: ssdata, AbstractStateSpace, Continuous, nstates, noutputs, ninputs # using ControlSystemIdentification using RobustAndOptimalControl, MonteCarloMeasurements -import ModelingToolkit: ODESystem, FnType, Symbolics +import ModelingToolkit: System, FnType, Symbolics using ModelingToolkit: unknowns, observed, isdifferential using Symbolics using Symbolics: jacobian, solve_for @@ -29,7 +29,7 @@ using UnPack # using SymbolicControlSystems -export sconnect, feedback, ODESystem, unknowns, observed, named_ss +export feedback, System, unknowns, observed, named_ss export batch_ss, trajectory_ss, GainScheduledStateSpace export build_quadratic_cost_matrix diff --git a/src/ode_system.jl b/src/ode_system.jl index 119f39b..29b6118 100644 --- a/src/ode_system.jl +++ b/src/ode_system.jl @@ -3,12 +3,12 @@ const AP = Union{Symbol, AnalysisPoint} import ModelingToolkitStandardLibrary.Blocks as Blocks conn = ModelingToolkit.connect t = Blocks.t -ModelingToolkit.ODESystem(sys::LTISystem; kwargs...) = ODESystem(ss(sys); kwargs...) +ModelingToolkit.System(sys::LTISystem; kwargs...) = System(ss(sys); kwargs...) """ - ModelingToolkit.ODESystem(sys::AbstractStateSpace; name::Symbol, x0 = zeros(sys.nx), x_names, u_names, y_names) + ModelingToolkit.System(sys::AbstractStateSpace; name::Symbol, x0 = zeros(sys.nx), x_names, u_names, y_names) -Create an ODESystem from `sys::StateSpace`. +Create an System from `sys::StateSpace`. # Arguments: - `sys`: An instance of `StateSpace` or `NamedStateSpace`. @@ -19,7 +19,7 @@ The arguments below are automatically set if the system is a `NamedStateSpace`. - `u_names`: A vector of symbols with input names. - `y_names`: A vector of symbols with output names. """ -function ModelingToolkit.ODESystem( +function ModelingToolkit.System( sys::AbstractStateSpace; name::Symbol, x0 = zeros(sys.nx), @@ -39,112 +39,14 @@ function ModelingToolkit.ODESystem( [uc[i].u ~ input.u[i] for i in 1:length(uc)]; [yc[i].u ~ output.u[i] for i in 1:length(yc)]; ] - extend(ODESystem(eqs, t; name, systems), ssblock) -end - -""" - sconnect(input::Function, sys::T; name) - sconnect(input::Num, sys::T; name) - -Connect a function `input(t)` to `sys.input` - -# Examples: -```julia -sconnect(sin, sys) # Connect a funciton, assumed to be a function of time -sconnect(sin(t), sys) # Connect a Num -``` -""" -function sconnect( - input::Union{Function, Num}, - sys::T; - name = Symbol("$(sys.name) with input"), -) where {T<:ModelingToolkit.AbstractTimeDependentSystem} - @named output = Blocks.RealOutput() - T( - [ - sys.input.u ~ (input isa Num ? input : input(t)) - output.u ~ sys.output.u - ], - t; - systems = [sys, output], - name, - ) -end - -""" - sconnect(sys1::T, sys2::T; name) - -Connect systems in series, equivalent to `sys2*sys1` or `series(sys1, sys2)` in ControlSystems.jl terminology -""" -function sconnect( - sys1::T, - sys2::T; - name = Symbol("$(sys1.name)*$(sys2.name)"), -) where {T<:ModelingToolkit.AbstractTimeDependentSystem} - @named output = Blocks.RealOutput() # TODO: missing size - @named input = Blocks.RealInput() # TODO: missing size - T( - [ - conn(input, sys2.input) - conn(output, sys1.output) - conn(sys2.output, sys1.input) - ], - t; - name, - systems = [sys1, sys2, output, input], - ) -end - -""" - G = ControlSystemsBase.feedback(loopgain::T; name) - -Form the feedback-interconnection -\$G = L/(1+L)\$ - -The system `G` will be a new system with `input` and `output` connectors. -""" -function ControlSystemsBase.feedback( - loopgain::T; - name = Symbol("feedback $(loopgain.name)"), -) where {T<:ModelingToolkit.AbstractTimeDependentSystem} - add = Blocks.Add(k1 = 1, k2 = -1, name = :feedback) - @named input = Blocks.RealInput() - @named output = Blocks.RealOutput() - T( - [ - input.u ~ add.input1.u - output.u ~ loopgain.output.u - conn(loopgain.output, add.input2) - conn(add.output, loopgain.input) - ], - t; - systems = [input, output, loopgain, add], - name, - ) -end - -function Base.:(*)(s1::T, s2::T) where {T<:ModelingToolkit.AbstractTimeDependentSystem} - name = Symbol(string(s1.name) * "_" * string(s2.name)) - @named input = Blocks.RealInput() - @named output = Blocks.RealOutput() - eqs = [ - conn(s1.input, s2.output) - output.u ~ s1.output.u - ] - systems = [output, s1, s2] - if any(s.name == :input for s in s2.systems) - push!(eqs, input.u ~ s2.input.u) - push!(systems, input) - end - T(eqs, t; systems, name) + extend(System(eqs, t; name, systems), ssblock) end numeric(x::Num) = x.val - function ControlSystemsBase.ss( - sys::ModelingToolkit.AbstractTimeDependentSystem, + sys::ModelingToolkit.AbstractSystem, inputs, outputs; kwargs... @@ -152,11 +54,13 @@ function ControlSystemsBase.ss( named_ss(sys, inputs, outputs; kwargs...).sys # just discard the names end +symstr(x) = Symbol(x isa AnalysisPoint ? x.name : string(x)) + """ RobustAndOptimalControl.named_ss(sys::ModelingToolkit.AbstractSystem, inputs, outputs; descriptor=true, simple_infeigs=true, kwargs...) -Convert an `ODESystem` to a `NamedStateSpace` using linearization. `inputs, outputs` are vectors of variables determining the inputs and outputs respectively. See docstring of `ModelingToolkit.linearize` for more info on `kwargs`. +Convert an `System` to a `NamedStateSpace` using linearization. `inputs, outputs` are vectors of variables determining the inputs and outputs respectively. See docstring of `ModelingToolkit.linearize` for more info on `kwargs`. If `descriptor = true` (default), this method automatically converts systems that MTK has failed to produce a proper form for into a proper linear statespace system using the method described here: https://juliacontrol.github.io/ControlSystemsMTK.jl/dev/#Internals:-Transformation-of-non-proper-models-to-proper-statespace-form @@ -166,7 +70,7 @@ If `descriptor = false`, the system is instead converted to a statespace realiza See also [`ModelingToolkit.linearize`](@ref) which is the lower-level function called internally. The functions [`get_named_sensitivity`](@ref), [`get_named_comp_sensitivity`](@ref), [`get_named_looptransfer`](@ref) similarily provide convenient ways to compute sensitivity functions while retaining signal names in the same way as `named_ss`. The corresponding lower-level functions `get_sensitivity`, `get_comp_sensitivity` and `get_looptransfer` are available in ModelingToolkitStandardLibrary.Blocks and are documented in [MTKstdlib: Linear analysis](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/linear_analysis/). """ function RobustAndOptimalControl.named_ss( - sys::ModelingToolkit.AbstractTimeDependentSystem, + sys::ModelingToolkit.AbstractSystem, inputs, outputs; descriptor = true, @@ -178,12 +82,12 @@ function RobustAndOptimalControl.named_ss( outputs = vcat(outputs) inputs = map(inputs) do inp - if inp isa ODESystem + if inp isa System @variables u(t) if u ∈ Set(unknowns(inp)) inp.u else - error("Input $(inp.name) is an ODESystem and not a variable") + error("Input $(inp.name) is an System and not a variable") end else inp @@ -192,19 +96,18 @@ function RobustAndOptimalControl.named_ss( nu = length(inputs) outputs = map(outputs) do out - if out isa ODESystem + if out isa System @variables u(t) if u ∈ Set(unknowns(out)) out.u else - error("Outut $(out.name) is an ODESystem and not a variable") + error("Outut $(out.name) is an System and not a variable") end else out end end - matrices, ssys = ModelingToolkit.linearize(sys, inputs, outputs; kwargs...) - symstr(x) = Symbol(x isa AnalysisPoint ? x.name : string(x)) + matrices, ssys, xpt = ModelingToolkit.linearize(sys, inputs, outputs; kwargs...) unames = symstr.(inputs) if nu > 0 && size(matrices.B, 2) == 2nu # This indicates that input derivatives are present @@ -214,6 +117,11 @@ function RobustAndOptimalControl.named_ss( else lsys = ss(matrices...) end + pind = [ModelingToolkit.parameter_index(ssys, i) for i in ModelingToolkit.inputs(ssys)] + x0 = xpt.x + u0 = [xpt.p[pi] for pi in pind] + xu = (; x = x0, u = u0) + extra = Dict(:operating_point => xu) # If simple_infeigs=false, the system might have been reduced and the state names might not match the original system. x_names = simple_infeigs ? symstr.(unknowns(ssys)) : [Symbol(string(nameof(sys))*"_x$i") for i in 1:lsys.nx] named_ss( @@ -222,6 +130,7 @@ function RobustAndOptimalControl.named_ss( u = unames, y = symstr.(outputs), name = string(Base.nameof(sys)), + extra, ) end @@ -292,7 +201,7 @@ get_named_looptransfer function named_sensitivity_function( fun, - sys::ModelingToolkit.AbstractTimeDependentSystem, + sys::ModelingToolkit.AbstractSystem, inputs, args...; descriptor = true, simple_infeigs = true, @@ -301,12 +210,12 @@ function named_sensitivity_function( inputs = vcat(inputs) inputs = map(inputs) do inp - if inp isa ODESystem + if inp isa System @variables u(t) if u ∈ Set(unknowns(inp)) inp.u else - error("Input $(inp.name) is an ODESystem and not a variable") + error("Input $(inp.name) is an System and not a variable") end else inp @@ -343,7 +252,7 @@ if isdefined(ModelingToolkit, :get_disturbance_system) end """ - build_quadratic_cost_matrix(linear_sys, ssys::ODESystem, costs::Vector{Pair}) + build_quadratic_cost_matrix(linear_sys, ssys::System, costs::Vector{Pair}) For a system that has been linearized, assemble a quadratic cost matrix (for LQR or Kalman filtering) that penalizes states or outputs of simplified system `ssys` according to the vector of pairs `costs`. @@ -358,7 +267,7 @@ The second problem above, the ordering of the states, can be worked around using - `ssys`: Output of [`linearize`](@ref). - `costs`: A vector of pairs """ -function build_quadratic_cost_matrix(matrices::NamedTuple, ssys::ODESystem, costs::AbstractVector{<:Pair}) +function build_quadratic_cost_matrix(matrices::NamedTuple, ssys::System, costs::AbstractVector{<:Pair}) x = ModelingToolkit.unknowns(ssys) y = ModelingToolkit.outputs(ssys) # y = getproperty.(ModelingToolkit.observed(ssys), :lhs) @@ -378,7 +287,7 @@ function build_quadratic_cost_matrix(matrices::NamedTuple, ssys::ODESystem, cost end """ - build_quadratic_cost_matrix(sys::ODESystem, inputs::Vector, costs::Vector{Pair}; kwargs...) + build_quadratic_cost_matrix(sys::System, inputs::Vector, costs::Vector{Pair}; kwargs...) Assemble a quadratic cost matrix (for LQR or Kalman filtering) that penalizes states or outputs of system `sys` according to the vector of pairs `costs`. @@ -393,8 +302,8 @@ The second problem above, the ordering of the states, can be worked around using - `inputs`: A vector of variables that are to be considered controlled inputs for the LQR controller. - `costs`: A vector of pairs. """ -function build_quadratic_cost_matrix(sys::ODESystem, inputs::AbstractVector, costs::AbstractVector{<:Pair}; kwargs...) - matrices, ssys = ModelingToolkit.linearize(sys, inputs, first.(costs); kwargs...) +function build_quadratic_cost_matrix(sys::System, inputs::AbstractVector, costs::AbstractVector{<:Pair}; kwargs...) + matrices, ssys, extras = ModelingToolkit.linearize(sys, inputs, first.(costs); kwargs...) x = ModelingToolkit.unknowns(ssys) y = ModelingToolkit.outputs(ssys) nx = length(x) @@ -417,10 +326,12 @@ function batch_linearize(sys, inputs, outputs, ops::AbstractVector{<:AbstractDic allow_input_derivatives = false, kwargs...) lin_fun, ssys = linearization_function(sys, inputs, outputs; op=ops[1], kwargs...) - lins = map(ops) do op + lins_ops = map(ops) do op linearize(ssys, lin_fun; op, t, allow_input_derivatives) end - lins, ssys + lins = first.(lins_ops) + resolved_ops = last.(lins_ops) + lins, ssys, resolved_ops end """ @@ -442,7 +353,7 @@ unsafe_comparisons(true) @variables x(t)=0 [bounds = (-2, 2)] @variables v(t)=0 @variables u(t)=0 -@variables y(t)=0 +@variables y(t) D = Differential(t) @@ -451,11 +362,12 @@ eqs = [D(x) ~ v y ~ x] -@named duffing = ODESystem(eqs, t) +@named duffing = System(eqs, t) bounds = getbounds(duffing, unknowns(duffing)) sample_within_bounds((l, u)) = (u - l) * rand() + l # Create a vector of operating points +N = 10 ops = map(1:N) do i op = Dict(x => sample_within_bounds(bounds[x]) for x in keys(bounds) if isfinite(bounds[x][1])) end @@ -498,10 +410,22 @@ nyquistcircles!(w, centers, radii, ylims = (-4, 1), xlims = (-3, 4)) See also [`trajectory_ss`](@ref) and [`fuzz`](@ref). """ function batch_ss(args...; kwargs...) - lins, ssys = batch_linearize(args...; kwargs...) - [ss(l...) for l in lins], ssys + lins, ssys, resolved_ops = batch_linearize(args...; kwargs...) + named_linsystems = map(lins) do l + # Convert to a NamedStateSpace with the same names as the original system + named_ss(ss(l.A, l.B, l.C, l.D); name = string(Base.nameof(ssys)), x = symstr.(unknowns(ssys))) + end + named_linsystems, ssys, resolved_ops end +# function unnamespace(ap) +# map(ap.outputs) do out +# ap_name = ModelingToolkit.SymbolicIndexingInterface.getname(out.u) +# new_name = join(ModelingToolkit.namespace_hierarchy(ap_name)[2:end], Symbolics.NAMESPACE_SEPARATOR) +# Symbolics.rename(ap.input.u, Symbol(new_name)) +# end +# end + """ linsystems, ssys = trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), fuzzer=nothing, verbose = true, kwargs...) @@ -515,19 +439,29 @@ Linearize `sys` around the trajectory `sol` at times `t`. Returns a vector of `S - `fuzzer`: A function that takes an operating point dictionary and returns an array of "fuzzed" operating points. This is useful for adding noise/uncertainty to the operating points along the trajectory. See [`ControlSystemsMTK.fuzz`](@ref) for such a function. - `verbose`: If `true`, print warnings for variables that are not found in `sol`. - `kwargs`: Are sent to the linearization functions. +- `named`: If `true`, the returned systems will be of type `NamedStateSpace`, otherwise they will be of type `StateSpace`. """ -function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), allow_input_derivatives = false, fuzzer = nothing, verbose = true, kwargs...) +function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), allow_input_derivatives = false, fuzzer = nothing, verbose = true, named = 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; warn_initialize_determined = false, kwargs...) - lin_fun, ssys = linearization_function(sys, inputs, outputs; warn_initialize_determined = false, kwargs...) - + lin_fun, ssys = linearization_function(sys, inputs, outputs; warn_empty_op = false, warn_initialize_determined = false, kwargs...) x = unknowns(ssys) + + # TODO: The value of the output (or input) of the input analysis points should be mapped to the perturbation vars + perturbation_vars = ModelingToolkit.inputs(ssys) + # original_inputs = reduce(vcat, unnamespace(ap) for ap in vcat(inputs)) # assuming all inputs are analysis points for now + + input_names = reduce(vcat, getproperty.(ap.outputs, :u) for ap in vcat(inputs)) + output_names = reduce(vcat, ap.input.u for ap in vcat(outputs)) + + op_nothing = Dict(unknowns(sys) .=> nothing) # Remove all defaults present in the original system defs = ModelingToolkit.defaults(sys) ops = map(t) do ti - Dict(x => robust_sol_getindex(sol, ti, x, defs; verbose) for x in x) + opsol = Dict(x => robust_sol_getindex(sol, ti, x, defs; verbose) for x in x) + # opsolu = Dict(new_u => robust_sol_getindex(sol, ti, u, defs; verbose) for (new_u, u) in zip(perturbation_vars, original_inputs)) + merge(op_nothing, opsol) end if fuzzer !== nothing opsv = map(ops) do op @@ -536,12 +470,30 @@ 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], initialize, kwargs...) - lins = map(zip(ops, t)) do (op, t) + lin_fun, ssys = linearization_function(sys, inputs, outputs; op=ops[1], kwargs...)#, initialization_abstol=1e-1, initialization_reltol=1e-1, kwargs...) # initializealg=ModelingToolkit.SciMLBase.NoInit() + # Main.lin_fun = lin_fun + # Main.op1 = ops[1] + # Main.ops = ops + # equations(lin_fun.prob.f.initialization_data.initializeprob.f.sys) + # observed(lin_fun.prob.f.initialization_data.initializeprob.f.sys) + lins_ops = 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] + # linearize(sys, inputs, outputs; op, t, allow_input_derivatives) # useful for debugging + end + lins = first.(lins_ops) + resolved_ops = last.(lins_ops) + named_linsystems = map(lins) do l + if named + # Convert to a NamedStateSpace with the same names as the original system + ynames = allunique(output_names) ? symstr.(output_names) : [Symbol(string(nameof(sys))*"_y$i") for i in 1:length(output_names)] + unames = allunique(input_names) ? symstr.(input_names) : [Symbol(string(nameof(sys))*"_u$i") for i in 1:length(input_names)] + nsys = named_ss(ss(l.A, l.B, l.C, l.D); name = string(Base.nameof(sys)), x = symstr.(unknowns(ssys)), u = unames, y = ynames) + # RobustAndOptimalControl.merge_nonunique_outputs(RobustAndOptimalControl.merge_nonunique_inputs(nsys)) + else + ss(l.A, l.B, l.C, l.D) + end end - (; linsystems = [ss(l...) for l in lins], ssys, ops) + (; linsystems = named_linsystems, ssys, ops, resolved_ops) end "_max_100(t) = length(t) > 100 ? range(extrema(t)..., 100) : t" @@ -663,10 +615,10 @@ function GainScheduledStateSpace(systems, vt; interpolator, x = zeros(systems[1] s1 = first(systems) (; nx, nu, ny) = s1 - Aint = [maybe_interp(interpolator, getindex.(getfield.(systems, :A), i, j), vt) for i = 1:nx, j = 1:nx] - Bint = [maybe_interp(interpolator, getindex.(getfield.(systems, :B), i, j), vt) for i = 1:nx, j = 1:nu] - Cint = [maybe_interp(interpolator, getindex.(getfield.(systems, :C), i, j), vt) for i = 1:ny, j = 1:nx] - Dint = [maybe_interp(interpolator, getindex.(getfield.(systems, :D), i, j), vt) for i = 1:ny, j = 1:nu] + Aint = [maybe_interp(interpolator, getindex.(getproperty.(systems, :A), i, j), vt) for i = 1:nx, j = 1:nx] + Bint = [maybe_interp(interpolator, getindex.(getproperty.(systems, :B), i, j), vt) for i = 1:nx, j = 1:nu] + Cint = [maybe_interp(interpolator, getindex.(getproperty.(systems, :C), i, j), vt) for i = 1:ny, j = 1:nx] + Dint = [maybe_interp(interpolator, getindex.(getproperty.(systems, :D), i, j), vt) for i = 1:ny, j = 1:nu] @named input = Blocks.RealInput(nin = nu) @named scheduling_input = Blocks.RealInput() @@ -678,10 +630,10 @@ function GainScheduledStateSpace(systems, vt; interpolator, x = zeros(systems[1] description = "Scheduling variable of gain-scheduled statespace system $name", ] - @variables A(v)[1:nx, 1:nx] = systems[1].A - @variables B(v)[1:nx, 1:nu] = systems[1].B - @variables C(v)[1:ny, 1:nx] = systems[1].C - @variables D(v)[1:ny, 1:nu] = systems[1].D + @variables A(t)[1:nx, 1:nx] = systems[1].A + @variables B(t)[1:nx, 1:nu] = systems[1].B + @variables C(t)[1:ny, 1:nx] = systems[1].C + @variables D(t)[1:ny, 1:nu] = systems[1].D A,B,C,D = collect.((A,B,C,D)) eqs = [ @@ -695,7 +647,7 @@ function GainScheduledStateSpace(systems, vt; interpolator, x = zeros(systems[1] for i in 1:nx]; collect(output.u .~ C * x .+ D * (input.u .- u0) .+ y0) ] - compose(ODESystem(eqs, t, name = name), [input, output, scheduling_input]) + compose(System(eqs, t, name = name), [input, output, scheduling_input]) end "LPVStateSpace is equivalent to GainScheduledStateSpace, see the docs for GainScheduledStateSpace." diff --git a/test/runtests.jl b/test/runtests.jl index a6e5c0a..bea6f84 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,8 +2,8 @@ using ControlSystemsMTK using Test @testset "ControlSystemsMTK.jl" begin - @testset "ODESystem" begin - @info "Testing ODESystem" + @testset "System" begin + @info "Testing System" include("test_ODESystem.jl") end diff --git a/test/test_ODESystem.jl b/test/test_ODESystem.jl index 3e426f0..438851f 100644 --- a/test/test_ODESystem.jl +++ b/test/test_ODESystem.jl @@ -11,34 +11,12 @@ connect = ModelingToolkit.connect P0 = tf(1.0, [1, 1]) |> ss C0 = pid(1, 1) * tf(1, [0.01, 1]) |> ss -@named P = ODESystem(P0) -@test P isa ODESystem +@named P = System(P0) +@test P isa System # @test length(ModelingToolkit.outputs(P)) == P0.ny # @test length(ModelingToolkit.inputs(P)) == P0.nu # @named nonlinear_P = sconnect(x->sign(x)*sqrt(abs(x)), P) # apply input-nonlinearity -@named C = ODESystem(C0) -@named loopgain = sconnect(C, P) -@named ref = Blocks.Sine(frequency = 1) -fb0 = feedback(loopgain, name = :fb) * ref -fb = structural_simplify(fb0) - -# @test length(unknowns(P)) == 3 # 1 + u + y -# @test length(unknowns(C)) == 4 # 2 + u + y - -x0 = Pair[ - collect(loopgain.P.x) .=> 1.0; -] - -prob = ODEProblem(fb, x0, (0.0, 10.0)) -sol = solve(prob, Rodas5()) -@test Int(sol.retcode) == 1 -isinteractive() && plot(sol) - -fb = structural_simplify(sconnect(sin, P)) -prob = ODEProblem(fb, x0, (0.0, 10.0)) -sol2 = solve(prob, Rodas5()) -@test Int(sol2.retcode) == 1 -isinteractive() && plot(sol2) +@named C = System(C0) Pc = complete(P) op = Dict(Pc.input.u => 0.0) @@ -65,7 +43,7 @@ P1 = ss(Pc, [Pc.input.u], [Pc.output.u]; op) -# === Go the other way, ODESystem -> StateSpace ================================ +# === Go the other way, System -> StateSpace ================================ x = unknowns(P) # I haven't figured out a good way to access states, so this is a bit manual and ugly @unpack input, output = P P02_named = named_ss(P, [input.u], [output.u]; op) @@ -79,11 +57,11 @@ P02 = ss(P02_named) # same for controller x = unknowns(C) @nonamespace C02 = named_ss(C, [C.input], [C.output]; op) -@test ss(C02) == C0 +@test tf(ss(C02)) ≈ tf(C0) -## Back again for a complete round trip, test that ODESystem get correct names -@named P2 = ODESystem(P02_named) +## Back again for a complete round trip, test that System get correct names +@named P2 = System(P02_named) # @test Set(unknowns(P)) ⊆ Set(unknowns(P2)) # @test Set(ModelingToolkit.inputs(P)) ⊆ Set(ModelingToolkit.inputs(P2)) # @test Set(ModelingToolkit.outputs(P)) ⊆ Set(ModelingToolkit.outputs(P2)) @@ -128,13 +106,13 @@ Cp0 = tx * pid(2, 0) * tf(1, [0.01, 1]) #|> ss P0i = deepcopy(P0) P0i.D .= 1e-8 -@named RF = ODESystem(tf(1, [0.001, 1])) # ref filter -@named RFv = ODESystem(tf(1, [0.001, 1])) -@named P = ODESystem(P0) # system model -@named Fv = ODESystem(inv(P0i[2, 1])) # feedforward vel -@named Fp = ODESystem(inv(P0i[1, 1])) # feedforward pos -@named Cv = ODESystem(Cv0) # vel controller -@named Cp = ODESystem(Cp0) # pos controller +@named RF = System(tf(1, [0.001, 1])) # ref filter +@named RFv = System(tf(1, [0.001, 1])) +@named P = System(P0) # system model +@named Fv = System(inv(P0i[2, 1])) # feedforward vel +@named Fp = System(inv(P0i[1, 1])) # feedforward pos +@named Cv = System(Cv0) # vel controller +@named Cp = System(Cp0) # pos controller fb = let ref0 = Blocks.Sine(amplitude = 0.2, frequency = 1, name = :r), disturbance = Blocks.Step(height = 100, start_time = 10, name = :d) @@ -142,7 +120,7 @@ fb = @named input = Blocks.RealInput() @named output = Blocks.RealOutput() Dₜ = Differential(t) - fb = ODESystem( + fb = System( [ conn(ref0.output, RF.input) # filter position reference expand_derivatives(Dₜ(ref0.output.u) ~ RFv.input.u) # Filter vel reference @@ -158,16 +136,15 @@ fb = name = :feedback, ) end -simplified_sys = structural_simplify(fb) +simplified_sys = mtkcompile(fb) x0 = Pair[ P.x[1] => 0.0 P.x[3] => 0.0 ] -p = Pair[] -prob = ODEProblem(simplified_sys, x0, (0.0, 20.0), p, jac = true) +prob = ODEProblem(simplified_sys, x0, (0.0, 20.0), jac = true) sol = solve(prob, Rodas5(), saveat = 0:0.01:20) if isinteractive() @show sol.retcode @@ -211,9 +188,9 @@ function SystemModel(u=nothing; name=:model) ] if u !== nothing push!(eqs, connect(u.output, :u, torque.tau)) - return @named model = ODESystem(eqs, t; systems = [sens, torque, inertia1, inertia2, spring, damper, u]) + return @named model = System(eqs, t; systems = [sens, torque, inertia1, inertia2, spring, damper, u]) end - ODESystem(eqs, t; systems = [sens, torque, inertia1, inertia2, spring, damper], name) + System(eqs, t; systems = [sens, torque, inertia1, inertia2, spring, damper], name) end model = SystemModel() |> complete @@ -237,7 +214,7 @@ sys = ss((mats...,)[1:4]...) defs = ModelingToolkit.defaults(ssys) defs = merge(Dict(unknowns(model) .=> 0), defs) -_, p = ModelingToolkit.get_u0_p(ssys, defs, defs) +p = ModelingToolkit.get_p(ssys, defs, split=false) sympars = ModelingToolkit.parameters(ssys) @@ -261,14 +238,14 @@ eqs = [connect(P.output, :plant_output, add.input2) connect(add.output, C.input) connect(C.output, :plant_input, P.input)] -sys_inner = ODESystem(eqs, t, systems = [P, C, add], name = :inner) +sys_inner = System(eqs, t, systems = [P, C, add], name = :inner) @named r = Blocks.Constant(k = 1) @named F = Blocks.FirstOrder(k = 1, T = 3) eqs = [connect(r.output, F.input) connect(F.output, sys_inner.add.input1)] -sys_outer = ODESystem(eqs, t, systems = [F, sys_inner, r], name = :outer) +sys_outer = System(eqs, t, systems = [F, sys_inner, r], name = :outer) matrices, _ = Blocks.get_sensitivity(sys_outer, [sys_outer.inner.plant_input, sys_outer.inner.plant_output]) S = ss(matrices...) @@ -282,7 +259,7 @@ Sn = get_named_sensitivity(sys_outer, [sys_outer.inner.plant_input, sys_outer.in ## Test connector names P = named_ss(ssrand(1,1,1), u=:jörgen, y=:solis) -@named Pode = ODESystem(P) +@named Pode = System(P) ModelingToolkit.isconnector(Pode.jörgen) ModelingToolkit.isconnector(Pode.solis) @@ -313,15 +290,20 @@ eqs = [connect(link1.TX1, cart.flange) connect(cart.flange, force.flange) connect(link1.TY1, fixed.flange)] -@named model = ODESystem(eqs, t, [], []; systems = [link1, cart, force, fixed]) +@named model = System(eqs, t, [], []; systems = [link1, cart, force, fixed]) lin_outputs = [cart.s, cart.v, link1.A, link1.dA] lin_inputs = [force.f.u] # => nothing to remove extra defaults op = Dict(cart.s => 10, cart.v => 0, link1.A => -pi/2, link1.dA => 0, force.f.u => 0, link1.x1 => nothing, link1.y1 => nothing, link1.x2 => nothing, link1.x_cm => nothing) +# op = Dict(cart.s => 10, cart.v => 0, link1.A => -pi/2, link1.dA => 0, force.f.u => 0) + +guesses = [link1.fy1 => 0.1, cart.f => 0.1] + G = named_ss(model, lin_inputs, lin_outputs; allow_symbolic = true, op, - allow_input_derivatives = true, zero_dummy_der = true) + allow_input_derivatives = true, zero_dummy_der = true, guesses) G = sminreal(G) +@test 10 ∈ RobustAndOptimalControl.operating_point(G).x @info "minreal" G = minreal(G) @info "poles" @@ -330,8 +312,8 @@ ps = poles(G) @test minimum(abs, ps) < 1e-6 @test minimum(abs, complex(0, 1.3777260367206716) .- ps) < 1e-10 -lsys, syss = linearize(model, lin_inputs, lin_outputs, op = op, - allow_input_derivatives = true) +lsys, syss = linearize(model, lin_inputs, lin_outputs; op = op, + allow_input_derivatives = true, guesses) lsyss, sysss = ModelingToolkit.linearize_symbolic(model, lin_inputs, lin_outputs; allow_input_derivatives = true) diff --git a/test/test_batchlin.jl b/test/test_batchlin.jl index 11804c4..e7cc035 100644 --- a/test/test_batchlin.jl +++ b/test/test_batchlin.jl @@ -19,7 +19,7 @@ eqs = [D(x) ~ v y.u ~ x] -@named duffing = ODESystem(eqs, t, systems=[y, u]) +@named duffing = System(eqs, t, systems=[y, u]) bounds = getbounds(duffing, unknowns(duffing)) sample_within_bounds((l, u)) = (u - l) * rand() + l @@ -32,14 +32,14 @@ inputs, outputs = [u.u], [y.u] Ps, ssys = batch_ss(duffing, inputs, outputs , ops) @test length(Ps) == N -@test Ps[1] == ss(linearize(duffing, inputs, outputs; op=ops[1])[1]...) -@test Ps[end] == ss(linearize(duffing, inputs, outputs; op=ops[end])[1]...) +@test Ps[1].sys == ss(linearize(duffing, inputs, outputs; op=ops[1])[1]...) +@test Ps[end].sys == ss(linearize(duffing, inputs, outputs; op=ops[end])[1]...) ## using DataInterpolations @named Cgs = GainScheduledStateSpace(Ps, xs, interpolator=LinearInterpolation) -@test Cgs isa ODESystem +@test Cgs isa System # This is tested better in the docs ## C-code generation @@ -53,7 +53,7 @@ import ModelingToolkitStandardLibrary.Blocks @named fb = Blocks.Add(k2=-1) @named ref = Blocks.Square(frequency=1/6, amplitude=0.5, offset=0.5, start_time=1) @named F = Blocks.SecondOrder(w=10, d=0.7) -@named C = ODESystem(pid(1,1,0; state_space=true, Tf=0.01)) +@named C = System(pid(1,1,0; state_space=true, Tf=0.01)) closed_loop_eqs = [ @@ -64,7 +64,7 @@ closed_loop_eqs = [ ModelingToolkit.connect(C.output, duffing.u) ] -@named closed_loop = ODESystem(closed_loop_eqs, t, systems=[duffing, C, fb, ref, F]) +@named closed_loop = System(closed_loop_eqs, t, systems=[duffing, C, fb, ref, F]) ssys = structural_simplify(closed_loop) prob = ODEProblem(ssys, [F.xd => 0], (0.0, 8.0))