Skip to content

Commit 9fd3de1

Browse files
authored
update to MTKv10 (#76)
* update to MTKv10 * bump cache version * up CI script * remove unused export * update other files * update usage of `get_u0_p` * more updates * turn off existing defaults * handle operating points * careful with state-space basis in interpolation * simplify implementation of `named_sensitivity_function` * add specification of op for original inputs does not work, since the original inputs have the wrong namespace * unnamespace * finally working * add a test * return named systems * add names to batch_ss as well * getfield -> getproperty * .sys * handle non-unique names
1 parent 78f50ad commit 9fd3de1

File tree

9 files changed

+266
-249
lines changed

9 files changed

+266
-249
lines changed

.github/workflows/CI.yml

Lines changed: 22 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -1,40 +1,45 @@
11
name: CI
22
on:
3-
- push
4-
- pull_request
3+
push:
4+
branches:
5+
- main
6+
tags: ['*']
7+
pull_request:
8+
workflow_dispatch:
9+
concurrency:
10+
# Skip intermediate builds: always.
11+
# Cancel intermediate builds: only if it is a pull request build.
12+
group: ${{ github.workflow }}-${{ github.ref }}
13+
cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }}
514
jobs:
615
test:
7-
name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
16+
name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}
817
runs-on: ${{ matrix.os }}
18+
timeout-minutes: 60
19+
permissions: # needed to allow julia-actions/cache to proactively delete old caches that it has created
20+
actions: write
21+
contents: read
922
strategy:
1023
fail-fast: false
1124
matrix:
1225
version:
13-
- '1'
26+
- '1.11'
1427
os:
1528
- ubuntu-latest
1629
arch:
1730
- x64
1831
steps:
19-
- uses: actions/checkout@v2
20-
- uses: julia-actions/setup-julia@v1
32+
- uses: actions/checkout@v4
33+
- uses: julia-actions/setup-julia@v2
2134
with:
2235
version: ${{ matrix.version }}
2336
arch: ${{ matrix.arch }}
24-
- uses: actions/cache@v4
25-
env:
26-
cache-name: cache-artifacts
27-
with:
28-
path: ~/.julia/artifacts
29-
key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
30-
restore-keys: |
31-
${{ runner.os }}-test-${{ env.cache-name }}-
32-
${{ runner.os }}-test-
33-
${{ runner.os }}-
37+
- uses: julia-actions/cache@v2
3438
- uses: julia-actions/julia-buildpkg@v1
3539
- uses: julia-actions/julia-runtest@v1
3640
- uses: julia-actions/julia-processcoverage@v1
3741
- uses: codecov/codecov-action@v4
3842
with:
39-
file: lcov.info
43+
files: lcov.info
4044
token: ${{ secrets.CODECOV_TOKEN }}
45+
fail_ci_if_error: false

Project.toml

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,8 @@
11
name = "ControlSystemsMTK"
22
uuid = "687d7614-c7e5-45fc-bfc3-9ee385575c88"
33
authors = ["Fredrik Bagge Carlson"]
4-
version = "2.3.3"
4+
version = "2.4.0"
5+
56

67
[deps]
78
ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
@@ -17,7 +18,7 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
1718
[compat]
1819
ControlSystemsBase = "1.0.1"
1920
DataInterpolations = "3, 4, 5, 6, 7"
20-
ModelingToolkit = "9.61"
21+
ModelingToolkit = "10.3"
2122
ModelingToolkitStandardLibrary = "2"
2223
MonteCarloMeasurements = "1.1"
2324
RobustAndOptimalControl = "0.4.14"

docs/src/batch_linearization.md

Lines changed: 93 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -25,13 +25,13 @@ eqs = [D(x) ~ v
2525
y.u ~ x]
2626
2727
28-
@named duffing = ODESystem(eqs, t, systems=[y, u])
28+
@named duffing = System(eqs, t, systems=[y, u])
2929
```
3030

3131
## Batch linearization
3232
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)
3333
```@example BATCHLIN
34-
N = 16 # Number of samples
34+
N = 5 # Number of samples
3535
xs = range(getbounds(x)[1], getbounds(x)[2], length=N)
3636
ops = Dict.(u.u => 0, x .=> xs)
3737
```
@@ -56,19 +56,24 @@ P = RobustAndOptimalControl.ss2particles(Ps) # convert to a single StateSpace sy
5656

5757
notice how some coefficients are plotted like uncertain numbers `-13.8 ± 4.3`. We can plot such models as well:
5858
```@example BATCHLIN
59-
bodeplot(P, w, legend=:bottomright) # Should look similar to the one above
59+
bodeplot(P, w, legend=:bottomright, adaptive=false) # Should look similar to the one above
6060
```
6161

6262
## Controller tuning
6363
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.
64+
65+
!!! "note" Interpolating between controllers
66+
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`.
67+
6468
```@example BATCHLIN
6569
function batch_tune(f, Ps)
6670
f.(Ps)
6771
end
6872
6973
Cs = batch_tune(Ps) do P
7074
C, kp, ki, kd, fig, CF = loopshapingPID(P, 7; Mt=1.2, Tf = 1/100)
71-
ss(CF)
75+
ss(CF, balance=false)
76+
modal_form(ss(CF, balance=true))[1]
7277
end
7378
7479
P = RobustAndOptimalControl.ss2particles(Ps)
@@ -103,22 +108,22 @@ using DataInterpolations # Required to interpolate between the controllers
103108
connect = ModelingToolkit.connect
104109
105110
closed_loop_eqs = [
106-
connect(ref.output, F.input)
111+
connect(ref.output, :r0, F.input)
107112
connect(F.output, :r, fb.input1) # Add an analysis point :r
108113
connect(duffing.y, :y, fb.input2) # Add an analysis point :y
109114
]
110115
plot(layout=2)
111116
112117
# Simulate each individual controller
113118
for C in Cs
114-
@named Ci = ODESystem(C)
119+
@named Ci = System(C)
115120
eqs = [
116121
closed_loop_eqs
117122
connect(fb.output, Ci.input)
118123
connect(Ci.output, duffing.u)
119124
]
120-
@named closed_loop = ODESystem(eqs, t, systems=[duffing, Ci, fb, ref, F])
121-
prob = ODEProblem(structural_simplify(closed_loop), [F.xd => 0], (0.0, 8.0))
125+
@named closed_loop = System(eqs, t, systems=[duffing, Ci, fb, ref, F])
126+
prob = ODEProblem(structural_simplify(closed_loop), [F.x => 0, F.xd => 0], (0.0, 8.0))
122127
sol = solve(prob, Rodas5P(), abstol=1e-8, reltol=1e-8)
123128
plot!(sol, idxs=[duffing.y.u, duffing.u.u], layout=2, lab="")
124129
end
@@ -128,39 +133,111 @@ end
128133
eqs = [
129134
closed_loop_eqs
130135
connect(fb.output, Cgs.input)
131-
connect(Cgs.output, duffing.u)
132-
connect(duffing.y, Cgs.scheduling_input) # Don't forget to connect the scheduling variable!
136+
connect(Cgs.output, :u, duffing.u)
137+
connect(duffing.y, :v, Cgs.scheduling_input) # Don't forget to connect the scheduling variable!
133138
]
134-
@named closed_loop = ODESystem(eqs, t, systems=[duffing, Cgs, fb, ref, F])
139+
@named closed_loop = System(eqs, t, systems=[duffing, Cgs, fb, ref, F])
135140
prob = ODEProblem(structural_simplify(closed_loop), [F.xd => 0], (0.0, 8.0))
136-
sol = solve(prob, Rodas5P(), abstol=1e-8, reltol=1e-8, initializealg=SciMLBase.NoInit())
141+
sol = solve(prob, Rodas5P(), abstol=1e-8, reltol=1e-8, initializealg=SciMLBase.NoInit(), dtmax=0.01)
137142
plot!(sol, idxs=[duffing.y.u, duffing.u.u], l=(2, :red), lab="Gain scheduled")
138143
plot!(sol, idxs=F.output.u, l=(1, :black, :dash, 0.5), lab="Ref")
139144
```
140145

141-
If everything worked as expected, the gain-scheduled controller should perform better than each of the included controllers individually.
146+
If everything worked as expected, the gain-scheduled controller should perform reasonably well across the entire scheduling range.
142147

143148

144149
## C-Code generation
145-
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]`
150+
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]`
146151
```@example BATCHLIN
147152
using SymbolicControlSystems, ControlSystemsBase
148153
Cs_disc = c2d.(Cs, 0.05, :tustin) # Discretize the controller before generating code
149-
code = SymbolicControlSystems.print_c_array(stdout, Cs_disc[1:7:end], xs[1:7:end], "Cs")
154+
code = SymbolicControlSystems.print_c_array(stdout, Cs_disc[1:3:end], xs[1:3:end], "Cs")
150155
```
151156
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).
152157

153158
## Linearize around a trajectory
154159
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`.
155160
```@example BATCHLIN
156161
timepoints = 0:0.01:8
157-
Ps2, ssys = trajectory_ss(closed_loop, closed_loop.r, closed_loop.y, sol; t=timepoints, initialize=true, verbose=true)
162+
Ps2, ssys, ops2, resolved_ops = trajectory_ss(closed_loop, closed_loop.r, closed_loop.y, sol; t=timepoints, verbose=true);
158163
bodeplot(Ps2, w, legend=false)
159164
```
165+
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.
160166

161167
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.
162168

163169

170+
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:
171+
- Break the connection from plant output to controller input by passing `loop_openings=[closed_loop.y]`
172+
- Break the connection between the controller and the plant input by passing `loop_openings=[closed_loop.u]`
173+
- 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]`
174+
175+
We will explore these options below, starting with the first option, breaking the connection `y`:
176+
```@example BATCHLIN
177+
kwargs = (; adaptive=false, legend=false)
178+
plants, _ = trajectory_ss(closed_loop, closed_loop.u, closed_loop.y, sol; t=timepoints, verbose=true, loop_openings=[closed_loop.u]);
179+
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]);
180+
181+
closed_loopsy = feedback.(plants .* controllersy)
182+
bodeplot(closed_loopsy, w; title="Loop open at y", kwargs...)
183+
```
184+
When we open the loop at `u`, we get a slightly different result:
185+
```@example BATCHLIN
186+
controllersu, ssu = trajectory_ss(closed_loop, closed_loop.r, closed_loop.u, sol; t=timepoints, verbose=true, loop_openings=[closed_loop.u]);
187+
188+
closed_loopsu = feedback.(plants .* controllersu)
189+
bodeplot(closed_loopsu, w; title="Loop open at u", kwargs...)
190+
```
191+
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`:
192+
```@example BATCHLIN
193+
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]);
194+
195+
closed_loopsv = feedback.(plants .* controllersv)
196+
bodeplot(closed_loopsv, w; title="Loop open at v and y", kwargs...)
197+
```
198+
199+
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
200+
201+
```@example BATCHLIN
202+
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]);
203+
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]);
204+
205+
closed_loops_mimo = feedback.(controllers_mimo .* plants_mimo) # Look at complementary sensitivity function in the input, since this is a SISO system
206+
bodeplot(closed_loops_mimo, w; title="Loop open at MIMO", kwargs...)
207+
```
208+
209+
210+
211+
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.
212+
213+
```@example BATCHLIN
214+
plot(
215+
bodeplot(Cs, w, legend=false, plotphase=false, title="Designed controllers"),
216+
bodeplot(controllersy, w, legend=false, plotphase=false, title="Loop open at y"),
217+
bodeplot(controllersu, w, legend=false, plotphase=false, title="Loop open at u"),
218+
bodeplot(controllersv, w, legend=false, plotphase=false, title="Loop open at v and y"),
219+
)
220+
```
221+
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).
222+
```@example BATCHLIN
223+
using Test
224+
@test all(sminreal.(controllersv) .== sminreal.(controllersu))
225+
```
226+
227+
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
228+
```@example BATCHLIN
229+
sminreal(controllersy[end]).x
230+
```
231+
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:
232+
```@example BATCHLIN
233+
sminreal(controllersu[end]).x
234+
```
235+
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.
236+
237+
238+
239+
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`.
240+
164241
## Summary
165242
We have seen how to
166243
- 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
173250

174251

175252
```@example BATCHLIN
176-
using Test
177253
@test sol(6.99, idxs=closed_loop.duffing.y.u) ≈ 0.0 atol=0.01
178254
```

0 commit comments

Comments
 (0)