Skip to content

Commit 8ee9a74

Browse files
authored
Merge pull request #1322 from baggepinnen/input_output
add functions to handle inputs and outputs
2 parents d130562 + a261271 commit 8ee9a74

File tree

5 files changed

+423
-0
lines changed

5 files changed

+423
-0
lines changed

src/ModelingToolkit.jl

Lines changed: 1 addition & 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")

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

src/variables.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,7 @@ struct Equality <: AbstractConnectType end # Equality connection
1616
struct Flow <: AbstractConnectType end # sum to 0
1717
struct Stream <: AbstractConnectType end # special stream connector
1818

19+
isvarkind(m, x::Num) = isvarkind(m, value(x))
1920
function isvarkind(m, x)
2021
p = getparent(x, nothing)
2122
p === nothing || (x = p)

0 commit comments

Comments
 (0)