Skip to content

Commit 89954e4

Browse files
committed
Formatter (2)
1 parent aecd59b commit 89954e4

File tree

2 files changed

+95
-74
lines changed

2 files changed

+95
-74
lines changed

docs/src/basics/Events.md

Lines changed: 90 additions & 71 deletions
Original file line numberDiff line numberDiff line change
@@ -379,36 +379,39 @@ sol.ps[c] # sol[c] will error, since `c` is not a timeseries value
379379

380380
It can be seen that the timeseries for `c` is not saved.
381381

382-
383382
## [(Experimental) Imperative affects](@id imp_affects)
383+
384384
The `ImperativeAffect` can be used as an alternative to the aforementioned functional affect form. Note
385385
that `ImperativeAffect` is still experimental; to emphasize this, we do not export it and it should be
386-
included as `ModelingToolkit.ImperativeAffect`. It abstracts over how values are written back to the
386+
included as `ModelingToolkit.ImperativeAffect`. It abstracts over how values are written back to the
387387
system, simplifying the definitions and (in the future) allowing assignments back to observed values
388388
by solving the nonlinear reinitialization problem afterwards.
389389

390-
We will use two examples to describe `ImperativeAffect`: a simple heater and a quadrature encoder.
390+
We will use two examples to describe `ImperativeAffect`: a simple heater and a quadrature encoder.
391391
These examples will also demonstrate advanced usage of `ModelingToolkit.SymbolicContinousCallback`,
392392
the low-level interface that the aforementioned tuple form converts into and allows control over the
393393
exact SciMLCallbacks event that is generated for a continous event.
394394

395395
### [Heater](@id heater_events)
396+
396397
Bang-bang control of a heater connected to a leaky plant requires hysteresis in order to prevent control oscillation.
397398

398-
```@example events
399+
```@example events
399400
@variables temp(t)
400401
params = @parameters furnace_on_threshold=0.5 furnace_off_threshold=0.7 furnace_power=1.0 leakage=0.1 furnace_on(t)::Bool=false
401402
eqs = [
402403
D(temp) ~ furnace_on * furnace_power - temp^2 * leakage
403404
]
404405
```
406+
405407
Our plant is simple. We have a heater that's turned on and off by the clocked parameter `furnace_on`
406408
which adds `furnace_power` forcing to the system when enabled. We then leak heat porportional to `leakage`
407-
as a function of the square of the current temperature.
409+
as a function of the square of the current temperature.
408410

409411
We need a controller with hysteresis to conol the plant. We wish the furnace to turn on when the temperature
410412
is below `furnace_on_threshold` and off when above `furnace_off_threshold`, while maintaining its current state
411413
in between. To do this, we create two continous callbacks:
414+
412415
```@example events
413416
using Setfield
414417
furnace_disable = ModelingToolkit.SymbolicContinuousCallback(
@@ -422,42 +425,49 @@ furnace_enable = ModelingToolkit.SymbolicContinuousCallback(
422425
@set! x.furnace_on = true
423426
end)
424427
```
428+
425429
We're using the explicit form of `SymbolicContinuousCallback` here, though
426430
so far we aren't using anything that's not possible with the implicit interface.
427-
You can also write
431+
You can also write
432+
428433
```julia
429-
[temp ~ furnace_off_threshold] => ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, i, c
434+
[temp ~ furnace_off_threshold] => ModelingToolkit.ImperativeAffect(modified = (;
435+
furnace_on)) do x, o, i, c
430436
@set! x.furnace_on = false
431437
end
432438
```
439+
433440
and it would work the same.
434441

435442
The `ImperativeAffect` is the larger change in this example. `ImperativeAffect` has the constructor signature
443+
436444
```julia
437-
ImperativeAffect(f::Function; modified::NamedTuple, observed::NamedTuple, ctx)
445+
ImperativeAffect(f::Function; modified::NamedTuple, observed::NamedTuple, ctx)
438446
```
447+
439448
that accepts the function to call, a named tuple of both the names of and symbolic values representing
440449
values in the system to be modified, a named tuple of the values that are merely observed (that is, used from
441450
the system but not modified), and a context that's passed to the affect function.
442451

443452
In our example, each event merely changes whether the furnace is on or off. Accordingly, we pass a `modified` tuple
444-
`(; furnace_on)` (creating a `NamedTuple` equivalent to `(furnace_on = furnace_on)`). `ImperativeAffect` will then
453+
`(; furnace_on)` (creating a `NamedTuple` equivalent to `(furnace_on = furnace_on)`). `ImperativeAffect` will then
445454
evaluate this before calling our function to fill out all of the numerical values, then apply them back to the system
446455
once our affect function returns. Furthermore, it will check that it is possible to do this assignment.
447456

448457
The function given to `ImperativeAffect` needs to have one of four signatures, checked in this order:
449-
* `f(modified::NamedTuple, observed::NamedTuple, ctx, integrator)::NamedTuple` if the function needs the low-level integrator,
450-
* `f(modified::NamedTuple, observed::NamedTuple, ctx)::NamedTuple` if the function needs the user-defined context,
451-
* `f(modified::NamedTuple, observed::NamedTuple)::NamedTuple` if the function also reads observed values from the system,
452-
* `f(modified::NamedTuple)::NamedTuple` if the function only writes values (unknowns or parameters) to the system.
453-
The `do` block in the example implicitly constructs said function inline. For exposition, we use the full version (e.g. `x, o, i, c`) but this could be simplified to merely `x`.
458+
459+
- `f(modified::NamedTuple, observed::NamedTuple, ctx, integrator)::NamedTuple` if the function needs the low-level integrator,
460+
- `f(modified::NamedTuple, observed::NamedTuple, ctx)::NamedTuple` if the function needs the user-defined context,
461+
- `f(modified::NamedTuple, observed::NamedTuple)::NamedTuple` if the function also reads observed values from the system,
462+
- `f(modified::NamedTuple)::NamedTuple` if the function only writes values (unknowns or parameters) to the system.
463+
The `do` block in the example implicitly constructs said function inline. For exposition, we use the full version (e.g. `x, o, i, c`) but this could be simplified to merely `x`.
454464

455465
The function `f` will be called with `observed` and `modified` `NamedTuple`s that are derived from their respective `NamedTuple` definitions.
456-
In our example, if `furnace_on` is `false`, then the value of the `x` that's passed in as `modified` will be `(furnace_on = false)`.
466+
In our example, if `furnace_on` is `false`, then the value of the `x` that's passed in as `modified` will be `(furnace_on = false)`.
457467
The modified values should be passed out in the same format: to set `furnace_on` to `true` we need to return a tuple `(furnace_on = true)`.
458468
We use Setfield to do this in the example, recreating the result tuple before returning it.
459469

460-
Accordingly, we can now interpret the `ImperativeAffect` definitions to mean that when `temp = furnace_off_threshold` we
470+
Accordingly, we can now interpret the `ImperativeAffect` definitions to mean that when `temp = furnace_off_threshold` we
461471
will write `furnace_on = false` back to the system, and when `temp = furnace_on_threshold` we will write `furnace_on = true` back
462472
to the system.
463473

@@ -468,7 +478,8 @@ ss = structural_simplify(sys)
468478
prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 10.0))
469479
sol = solve(prob, Tsit5())
470480
plot(sol)
471-
hline!([sol.ps[furnace_off_threshold], sol.ps[furnace_on_threshold]], l = (:black, 1), primary = false)
481+
hline!([sol.ps[furnace_off_threshold], sol.ps[furnace_on_threshold]],
482+
l = (:black, 1), primary = false)
472483
```
473484

474485
Here we see exactly the desired hysteresis. The heater starts on until the temperature hits
@@ -477,71 +488,76 @@ point the furnace turns on again until `furnace_off_threshold` and so on and so
477488
is effectively regulating the temperature of the plant.
478489

479490
### [Quadrature Encoder](@id quadrature)
491+
480492
For a more complex application we'll look at modeling a quadrature encoder attached to a shaft spinning at a constant speed.
481493
Traditionally, a quadrature encoder is built out of a code wheel that interrupts the sensors at constant intervals and two sensors slightly out of phase with one another.
482494
A state machine can take the pattern of pulses produced by the two sensors and determine the number of steps that the shaft has spun. The state machine takes the new value
483495
from each sensor and the old values and decodes them into the direction that the wheel has spun in this step.
484496

485497
```@example events
486-
@variables theta(t) omega(t)
487-
params = @parameters qA=0 qB=0 hA=0 hB=0 cnt::Int=0
488-
eqs = [D(theta) ~ omega
489-
omega ~ 1.0]
498+
@variables theta(t) omega(t)
499+
params = @parameters qA=0 qB=0 hA=0 hB=0 cnt::Int=0
500+
eqs = [D(theta) ~ omega
501+
omega ~ 1.0]
490502
```
503+
491504
Our continous-time system is extremely simple. We have two states, `theta` for the angle of the shaft
492-
and `omega` for the rate at which it's spinning. We then have parameters for the state machine `qA, qB, hA, hB`
505+
and `omega` for the rate at which it's spinning. We then have parameters for the state machine `qA, qB, hA, hB`
493506
and a step count `cnt`.
494507

495508
We'll then implement the decoder as a simple Julia function.
509+
496510
```@example events
497-
function decoder(oldA, oldB, newA, newB)
498-
state = (oldA, oldB, newA, newB)
499-
if state == (0, 0, 1, 0) || state == (1, 0, 1, 1) || state == (1, 1, 0, 1) ||
500-
state == (0, 1, 0, 0)
501-
return 1
502-
elseif state == (0, 0, 0, 1) || state == (0, 1, 1, 1) || state == (1, 1, 1, 0) ||
503-
state == (1, 0, 0, 0)
504-
return -1
505-
elseif state == (0, 0, 0, 0) || state == (0, 1, 0, 1) || state == (1, 0, 1, 0) ||
506-
state == (1, 1, 1, 1)
507-
return 0
508-
else
509-
return 0 # err is interpreted as no movement
510-
end
511+
function decoder(oldA, oldB, newA, newB)
512+
state = (oldA, oldB, newA, newB)
513+
if state == (0, 0, 1, 0) || state == (1, 0, 1, 1) || state == (1, 1, 0, 1) ||
514+
state == (0, 1, 0, 0)
515+
return 1
516+
elseif state == (0, 0, 0, 1) || state == (0, 1, 1, 1) || state == (1, 1, 1, 0) ||
517+
state == (1, 0, 0, 0)
518+
return -1
519+
elseif state == (0, 0, 0, 0) || state == (0, 1, 0, 1) || state == (1, 0, 1, 0) ||
520+
state == (1, 1, 1, 1)
521+
return 0
522+
else
523+
return 0 # err is interpreted as no movement
511524
end
525+
end
512526
```
527+
513528
Based on the current and old state, this function will return 1 if the wheel spun in the positive direction,
514529
-1 if in the negative, and 0 otherwise.
515530

516-
The encoder state advances when the occlusion begins or ends. We model the
531+
The encoder state advances when the occlusion begins or ends. We model the
517532
code wheel as simply detecting when `cos(100*theta)` is 0; if we're at a positive
518533
edge of the 0 crossing, then we interpret that as occlusion (so the discrete `qA` goes to 1). Otherwise, if `cos` is
519534
going negative, we interpret that as lack of occlusion (so the discrete goes to 0). The decoder function is
520535
then invoked to update the count with this new information.
521536

522537
We can implement this in one of two ways: using edge sign detection or right root finding. For exposition, we
523-
will implement each sensor differently.
538+
will implement each sensor differently.
524539

525540
For sensor A, we're using the edge detction method. By providing a different affect to `SymbolicContinuousCallback`'s
526541
`affect_neg` argument, we can specify different behaviour for the negative crossing vs. the positive crossing of the root.
527542
In our encoder, we interpret this as occlusion or nonocclusion of the sensor, update the internal state, and tick the decoder.
543+
528544
```@example events
529-
qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta) ~ 0],
530-
ModelingToolkit.ImperativeAffect((; qA, hA, hB, cnt), (; qB)) do x, o, i, c
531-
@set! x.hA = x.qA
532-
@set! x.hB = o.qB
533-
@set! x.qA = 1
534-
@set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB)
535-
x
536-
end,
537-
affect_neg = ModelingToolkit.ImperativeAffect(
538-
(; qA, hA, hB, cnt), (; qB)) do x, o, c, i
539-
@set! x.hA = x.qA
540-
@set! x.hB = o.qB
541-
@set! x.qA = 0
542-
@set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB)
543-
x
544-
end)
545+
qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta) ~ 0],
546+
ModelingToolkit.ImperativeAffect((; qA, hA, hB, cnt), (; qB)) do x, o, i, c
547+
@set! x.hA = x.qA
548+
@set! x.hB = o.qB
549+
@set! x.qA = 1
550+
@set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB)
551+
x
552+
end,
553+
affect_neg = ModelingToolkit.ImperativeAffect(
554+
(; qA, hA, hB, cnt), (; qB)) do x, o, c, i
555+
@set! x.hA = x.qA
556+
@set! x.hB = o.qB
557+
@set! x.qA = 0
558+
@set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB)
559+
x
560+
end)
545561
```
546562

547563
The other way we can implement a sensor is by changing the root find.
@@ -550,29 +566,32 @@ the root is crossed. This makes it trickier to figure out what the new state is.
550566
Instead, we can use right root finding:
551567

552568
```@example events
553-
qBevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta - π / 2) ~ 0],
554-
ModelingToolkit.ImperativeAffect((; qB, hA, hB, cnt), (; qA, theta)) do x, o, i, c
555-
@set! x.hA = o.qA
556-
@set! x.hB = x.qB
557-
@set! x.qB = clamp(sign(cos(100 * o.theta - π / 2)), 0.0, 1.0)
558-
@set! x.cnt += decoder(x.hA, x.hB, o.qA, x.qB)
559-
x
560-
end; rootfind = SciMLBase.RightRootFind)
569+
qBevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta - π / 2) ~ 0],
570+
ModelingToolkit.ImperativeAffect((; qB, hA, hB, cnt), (; qA, theta)) do x, o, i, c
571+
@set! x.hA = o.qA
572+
@set! x.hB = x.qB
573+
@set! x.qB = clamp(sign(cos(100 * o.theta - π / 2)), 0.0, 1.0)
574+
@set! x.cnt += decoder(x.hA, x.hB, o.qA, x.qB)
575+
x
576+
end; rootfind = SciMLBase.RightRootFind)
561577
```
562-
Here, sensor B is located `π / 2` behind sensor A in angular space, so we're adjusting our
578+
579+
Here, sensor B is located `π / 2` behind sensor A in angular space, so we're adjusting our
563580
trigger function accordingly. We here ask for right root finding on the callback, so we know
564-
that the value of said function will have the "new" sign rather than the old one. Thus, we can
581+
that the value of said function will have the "new" sign rather than the old one. Thus, we can
565582
determine the new state of the sensor from the sign of the indicator function evaluated at the
566583
affect activation point, with -1 mapped to 0.
567584

568585
We can now simulate the encoder.
586+
569587
```@example events
570-
@named sys = ODESystem(
571-
eqs, t, [theta, omega], params; continuous_events = [qAevt, qBevt])
572-
ss = structural_simplify(sys)
573-
prob = ODEProblem(ss, [theta => 0.0], (0.0, pi))
574-
sol = solve(prob, Tsit5(); dtmax = 0.01)
575-
sol.ps[cnt]
588+
@named sys = ODESystem(
589+
eqs, t, [theta, omega], params; continuous_events = [qAevt, qBevt])
590+
ss = structural_simplify(sys)
591+
prob = ODEProblem(ss, [theta => 0.0], (0.0, pi))
592+
sol = solve(prob, Tsit5(); dtmax = 0.01)
593+
sol.ps[cnt]
576594
```
595+
577596
`cos(100*theta)` will have 200 crossings in the half rotation we've gone through, so the encoder would notionally count 200 steps.
578-
Our encoder counts 198 steps (it loses one step to initialization and one step due to the final state falling squarely on an edge).
597+
Our encoder counts 198 steps (it loses one step to initialization and one step due to the final state falling squarely on an edge).

src/systems/diffeqs/odesystem.jl

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -444,7 +444,7 @@ function build_explicit_observed_function(sys, ts;
444444
param_only = false,
445445
op = Operator,
446446
throw = true,
447-
array_type=:array)
447+
array_type = :array)
448448
if (isscalar = symbolic_type(ts) !== NotSymbolic())
449449
ts = [ts]
450450
end
@@ -589,10 +589,12 @@ function build_explicit_observed_function(sys, ts;
589589
oop_mtkp_wrapper = mtkparams_wrapper
590590
end
591591

592-
output_expr = isscalar ? ts[1] : (array_type == :array ? MakeArray(ts, output_type) : MakeTuple(ts))
592+
output_expr = isscalar ? ts[1] :
593+
(array_type == :array ? MakeArray(ts, output_type) : MakeTuple(ts))
593594
# Need to keep old method of building the function since it uses `output_type`,
594595
# which can't be provided to `build_function`
595-
oop_fn = Func(args, [], pre(Let(obsexprs, output_expr, false))) |> array_wrapper[1] |> oop_mtkp_wrapper |> toexpr
596+
oop_fn = Func(args, [], pre(Let(obsexprs, output_expr, false))) |> array_wrapper[1] |>
597+
oop_mtkp_wrapper |> toexpr
596598
oop_fn = expression ? oop_fn : eval_or_rgf(oop_fn; eval_expression, eval_module)
597599

598600
if !isscalar

0 commit comments

Comments
 (0)