Skip to content

Commit f1032d7

Browse files
committed
Merge branch 'master' into myb/error_check
2 parents 798e6e7 + ee8890b commit f1032d7

File tree

13 files changed

+803
-6
lines changed

13 files changed

+803
-6
lines changed

Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "ModelingToolkit"
22
uuid = "961ee093-0014-501f-94e3-6117800e7a78"
33
authors = ["Chris Rackauckas <[email protected]>"]
4-
version = "8.9.0"
4+
version = "8.10.0"
55

66
[deps]
77
AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"

docs/make.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,7 @@ makedocs(
3030
"Basics" => Any[
3131
"basics/AbstractSystem.md",
3232
"basics/ContextualVariables.md",
33+
"basics/Variable_metadata.md",
3334
"basics/Composition.md",
3435
"basics/Validation.md",
3536
"basics/DependencyGraphs.md",

docs/src/basics/Variable_metadata.md

Lines changed: 95 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,95 @@
1+
# Symbolic metadata
2+
It is possible to add metadata to symbolic variables. The following
3+
information can be added (note, it's possible to extend this to user-defined metadata as well)
4+
5+
## Input or output
6+
Designate a variable as either an input or an output using the following
7+
```@example metadata
8+
using ModelingToolkit
9+
@variables u [input=true]
10+
isinput(u)
11+
```
12+
```@example metadata
13+
@variables y [output=true]
14+
isoutput(y)
15+
```
16+
17+
## Bounds
18+
Bounds are useful when parameters are to be optimized, or to express intervals of uncertainty.
19+
20+
```@example metadata
21+
@variables u [bounds=(-1,1)]
22+
hasbounds(u)
23+
```
24+
```@example metadata
25+
getbounds(u)
26+
```
27+
28+
## Mark input as a disturbance
29+
Indicate that an input is not available for control, i.e., it's a disturbance input.
30+
31+
```@example metadata
32+
@variables u [input=true, disturbance=true]
33+
isdisturbance(u)
34+
```
35+
36+
## Mark parameter as tunable
37+
Indicate that a parameter can be automatically tuned by automatic control tuning apps.
38+
39+
```@example metadata
40+
@parameters Kp [tunable=true]
41+
istunable(Kp)
42+
```
43+
44+
## Probability distributions
45+
A probability distribution may be associated with a parameter to indicate either
46+
uncertainty about it's value, or as a prior distribution for Bayesian optimization.
47+
48+
```julia
49+
using Distributions
50+
d = Normal(10, 1)
51+
@parameters m [dist=d]
52+
hasdist(m)
53+
```
54+
```julia
55+
getdist(m)
56+
```
57+
58+
## Additional functions
59+
For systems that contain parameters with metadata like described above have some additional functions defined for convenience.
60+
In the example below, we define a system with tunable parameters and extract bounds vectors
61+
62+
```@example metadata
63+
@parameters t
64+
Dₜ = Differential(t)
65+
@variables x(t)=0 u(t)=0 [input=true] y(t)=0 [output=true]
66+
@parameters T [tunable = true, bounds = (0, Inf)]
67+
@parameters k [tunable = true, bounds = (0, Inf)]
68+
eqs = [
69+
Dₜ(x) ~ (-x + k*u) / T # A first-order system with time constant T and gain k
70+
y ~ x
71+
]
72+
sys = ODESystem(eqs, t, name=:tunable_first_order)
73+
```
74+
```@example metadata
75+
p = tunable_parameters(sys) # extract all parameters marked as tunable
76+
```
77+
```@example metadata
78+
lb, ub = getbounds(p) # operating on a vector, we get lower and upper bound vectors
79+
```
80+
```@example metadata
81+
b = getbounds(sys) # Operating on the system, we get a dict
82+
```
83+
84+
85+
## Index
86+
```@index
87+
Pages = ["Variable_metadata.md"]
88+
```
89+
90+
## Docstrings
91+
```@autodocs
92+
Modules = [ModelingToolkit]
93+
Pages = ["variables.jl"]
94+
Private = false
95+
```

src/ModelingToolkit.jl

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -115,6 +115,7 @@ using .BipartiteGraphs
115115

116116
include("variables.jl")
117117
include("parameters.jl")
118+
include("inputoutput.jl")
118119

119120
include("utils.jl")
120121
include("domains.jl")
@@ -172,6 +173,7 @@ export NonlinearSystem, OptimizationSystem
172173
export ControlSystem
173174
export alias_elimination, flatten
174175
export connect, @connector, Connection, Flow, Stream, instream
176+
export isinput, isoutput, getbounds, hasbounds, isdisturbance, istunable, getdist, hasdist, tunable_parameters
175177
export ode_order_lowering, liouville_transform
176178
export runge_kutta_discretize
177179
export PDESystem

src/inputoutput.jl

Lines changed: 248 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,248 @@
1+
using Symbolics: get_variables
2+
"""
3+
inputs(sys)
4+
5+
Return all variables that mare marked as inputs. See also [`unbound_inputs`](@ref)
6+
See also [`bound_inputs`](@ref), [`unbound_inputs`](@ref)
7+
"""
8+
inputs(sys) = filter(isinput, states(sys))
9+
10+
"""
11+
outputs(sys)
12+
13+
Return all variables that mare marked as outputs. See also [`unbound_outputs`](@ref)
14+
See also [`bound_outputs`](@ref), [`unbound_outputs`](@ref)
15+
"""
16+
function outputs(sys)
17+
o = observed(sys)
18+
rhss = [eq.rhs for eq in o]
19+
lhss = [eq.lhs for eq in o]
20+
unique([
21+
filter(isoutput, states(sys))
22+
filter(x -> x isa Term && isoutput(x), rhss) # observed can return equations with complicated expressions, we are only looking for single Terms
23+
filter(x -> x isa Term && isoutput(x), lhss)
24+
])
25+
end
26+
27+
"""
28+
bound_inputs(sys)
29+
30+
Return inputs that are bound within the system, i.e., internal inputs
31+
See also [`bound_inputs`](@ref), [`unbound_inputs`](@ref), [`bound_outputs`](@ref), [`unbound_outputs`](@ref)
32+
"""
33+
bound_inputs(sys) = filter(x->is_bound(sys, x), inputs(sys))
34+
35+
"""
36+
unbound_inputs(sys)
37+
38+
Return inputs that are not bound within the system, i.e., external inputs
39+
See also [`bound_inputs`](@ref), [`unbound_inputs`](@ref), [`bound_outputs`](@ref), [`unbound_outputs`](@ref)
40+
"""
41+
unbound_inputs(sys) = filter(x->!is_bound(sys, x), inputs(sys))
42+
43+
"""
44+
bound_outputs(sys)
45+
46+
Return outputs that are bound within the system, i.e., internal outputs
47+
See also [`bound_inputs`](@ref), [`unbound_inputs`](@ref), [`bound_outputs`](@ref), [`unbound_outputs`](@ref)
48+
"""
49+
bound_outputs(sys) = filter(x->is_bound(sys, x), outputs(sys))
50+
51+
"""
52+
unbound_outputs(sys)
53+
54+
Return outputs that are not bound within the system, i.e., external outputs
55+
See also [`bound_inputs`](@ref), [`unbound_inputs`](@ref), [`bound_outputs`](@ref), [`unbound_outputs`](@ref)
56+
"""
57+
unbound_outputs(sys) = filter(x->!is_bound(sys, x), outputs(sys))
58+
59+
"""
60+
is_bound(sys, u)
61+
62+
Determine whether or not input/output variable `u` is "bound" within the system, i.e., if it's to be considered internal to `sys`.
63+
A variable/signal is considered bound if it appears in an equation together with variables from other subsystems.
64+
The typical usecase for this function is to determine whether the input to an IO component is connected to another component,
65+
or if it remains an external input that the user has to supply before simulating the system.
66+
67+
See also [`bound_inputs`](@ref), [`unbound_inputs`](@ref), [`bound_outputs`](@ref), [`unbound_outputs`](@ref)
68+
"""
69+
function is_bound(sys, u, stack=[])
70+
#=
71+
For observed quantities, we check if a variable is connected to something that is bound to something further out.
72+
In the following scenario
73+
julia> observed(syss)
74+
2-element Vector{Equation}:
75+
sys₊y(tv) ~ sys₊x(tv)
76+
y(tv) ~ sys₊x(tv)
77+
sys₊y(t) is bound to the outer y(t) through the variable sys₊x(t) and should thus return is_bound(sys₊y(t)) = true.
78+
When asking is_bound(sys₊y(t)), we know that we are looking through observed equations and can thus ask
79+
if var is bound, if it is, then sys₊y(t) is also bound. This can lead to an infinite recursion, so we maintain a stack of variables we have previously asked about to be able to break cycles
80+
=#
81+
u Set(stack) && return false # Cycle detected
82+
eqs = equations(sys)
83+
eqs = filter(eq->has_var(eq, u), eqs) # Only look at equations that contain u
84+
# isout = isoutput(u)
85+
for eq in eqs
86+
vars = [get_variables(eq.rhs); get_variables(eq.lhs)]
87+
for var in vars
88+
var === u && continue
89+
if !same_or_inner_namespace(u, var)
90+
return true
91+
end
92+
end
93+
end
94+
# Look through observed equations as well
95+
oeqs = observed(sys)
96+
oeqs = filter(eq->has_var(eq, u), oeqs) # Only look at equations that contain u
97+
for eq in oeqs
98+
vars = [get_variables(eq.rhs); get_variables(eq.lhs)]
99+
for var in vars
100+
var === u && continue
101+
if !same_or_inner_namespace(u, var)
102+
return true
103+
end
104+
if is_bound(sys, var, [stack; u]) && !inner_namespace(u, var) # The variable we are comparing to can not come from an inner namespace, binding only counts outwards
105+
return true
106+
end
107+
end
108+
end
109+
false
110+
end
111+
112+
113+
114+
"""
115+
same_or_inner_namespace(u, var)
116+
117+
Determine whether or not `var` is in the same namespace as `u`, or a namespace internal to the namespace of `u`.
118+
Example: `sys.u ~ sys.inner.u` will bind `sys.inner.u`, but `sys.u` remains an unbound, external signal. The namepsaced signal `sys.inner.u` lives in a namspace internal to `sys`.
119+
"""
120+
function same_or_inner_namespace(u, var)
121+
nu = get_namespace(u)
122+
nv = get_namespace(var)
123+
nu == nv || # namespaces are the same
124+
startswith(nv, nu) || # or nv starts with nu, i.e., nv is an inner namepsace to nu
125+
occursin('', var) && !occursin('', u) # or u is top level but var is internal
126+
end
127+
128+
function inner_namespace(u, var)
129+
nu = get_namespace(u)
130+
nv = get_namespace(var)
131+
nu == nv && return false
132+
startswith(nv, nu) || # or nv starts with nu, i.e., nv is an inner namepsace to nu
133+
occursin('', var) && !occursin('', u) # or u is top level but var is internal
134+
end
135+
136+
"""
137+
get_namespace(x)
138+
139+
Return the namespace of a variable as a string. If the variable is not namespaced, the string is empty.
140+
"""
141+
function get_namespace(x)
142+
sname = string(x)
143+
parts = split(sname, '')
144+
if length(parts) == 1
145+
return ""
146+
end
147+
join(parts[1:end-1], '')
148+
end
149+
150+
"""
151+
has_var(eq, x)
152+
153+
Determine whether or not an equation or expression contains variable `x`.
154+
"""
155+
function has_var(eq::Equation, x)
156+
has_var(eq.rhs, x) || has_var(eq.lhs, x)
157+
end
158+
159+
has_var(ex, x) = x Set(get_variables(ex))
160+
161+
# Build control function
162+
163+
"""
164+
(f_oop, f_ip), dvs, p = generate_control_function(sys::AbstractODESystem, dvs = states(sys), ps = parameters(sys); implicit_dae = false, ddvs = if implicit_dae
165+
166+
For a system `sys` that has unbound inputs (as determined by [`unbound_inputs`](@ref)), generate a function with additional input argument `in`
167+
```
168+
f_oop : (u,in,p,t) -> rhs
169+
f_ip : (uout,u,in,p,t) -> nothing
170+
```
171+
The return values also include the remaining states and parameters, in the order they appear as arguments to `f`.
172+
173+
# Example
174+
```
175+
using ModelingToolkit: generate_control_function, varmap_to_vars, defaults
176+
f, dvs, ps = generate_control_function(sys, expression=Val{false}, simplify=true)
177+
p = varmap_to_vars(defaults(sys), ps)
178+
x = varmap_to_vars(defaults(sys), dvs)
179+
t = 0
180+
f[1](x, inputs, p, t)
181+
```
182+
"""
183+
function generate_control_function(
184+
sys::AbstractODESystem;
185+
implicit_dae=false,
186+
has_difference=false,
187+
simplify=true,
188+
kwargs...
189+
)
190+
191+
ctrls = unbound_inputs(sys)
192+
if isempty(ctrls)
193+
error("No unbound inputs were found in system.")
194+
end
195+
196+
# One can either connect unbound inputs to new parameters and allow structural_simplify, but then the unbound inputs appear as states :( .
197+
# One can also just remove them from the states and parameters for the purposes of code generation, but then structural_simplify fails :(
198+
# To have the best of both worlds, all unbound inputs must be converted to `@parameters` in which case structural_simplify handles them correctly :)
199+
sys = toparam(sys, ctrls)
200+
201+
if simplify
202+
sys = structural_simplify(sys)
203+
end
204+
205+
dvs = states(sys)
206+
ps = parameters(sys)
207+
208+
dvs = setdiff(dvs, ctrls)
209+
ps = setdiff(ps, ctrls)
210+
inputs = map(x->time_varying_as_func(value(x), sys), ctrls)
211+
212+
eqs = [eq for eq in equations(sys) if !isdifferenceeq(eq)]
213+
foreach(check_derivative_variables, eqs)
214+
# substitute x(t) by just x
215+
rhss = implicit_dae ? [_iszero(eq.lhs) ? eq.rhs : eq.rhs - eq.lhs for eq in eqs] :
216+
[eq.rhs for eq in eqs]
217+
218+
219+
# TODO: add an optional check on the ordering of observed equations
220+
u = map(x->time_varying_as_func(value(x), sys), dvs)
221+
p = map(x->time_varying_as_func(value(x), sys), ps)
222+
t = get_iv(sys)
223+
224+
# pre = has_difference ? (ex -> ex) : get_postprocess_fbody(sys)
225+
226+
args = (u, inputs, p, t)
227+
if implicit_dae
228+
ddvs = map(Differential(get_iv(sys)), dvs)
229+
args = (ddvs, args...)
230+
end
231+
pre, sol_states = get_substitutions_and_solved_states(sys)
232+
f = build_function(rhss, args...; postprocess_fbody=pre, states=sol_states, kwargs...)
233+
f, dvs, ps
234+
end
235+
236+
"""
237+
toparam(sys, ctrls::AbstractVector)
238+
239+
Transform all instances of `@varibales` in `ctrls` appearing as states and in equations of `sys` with similarly named `@parameters`. This allows [`structural_simplify`](@ref)(sys) in the presence unbound inputs.
240+
"""
241+
function toparam(sys, ctrls::AbstractVector)
242+
eqs = equations(sys)
243+
subs = Dict(ctrls .=> toparam.(ctrls))
244+
eqs = map(eqs) do eq
245+
substitute(eq.lhs, subs) ~ substitute(eq.rhs, subs)
246+
end
247+
ODESystem(eqs, name=sys.name)
248+
end

0 commit comments

Comments
 (0)