Skip to content

Commit 4d9cc94

Browse files
RXGottliebRXGottlieb
committed
McCormick rule update; general cleanup
--Added McCormick transformation rules --Cleaned up deprecated Expr handling --Removed some extraneous comments --Added functionality for extracting model terms and setting bounds
1 parent 694dcfa commit 4d9cc94

File tree

9 files changed

+666
-238
lines changed

9 files changed

+666
-238
lines changed

‎robert_test_stuff.jl‎

Lines changed: 172 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,172 @@
1+
using ModelingToolkit
2+
3+
@variables t x1(t) x2(t) x3(t) x4(t)
4+
@parameters p1 p2 p3 p4 p5
5+
D = Differential(t)
6+
7+
@named testcase = ODESystem(
8+
[D(x1) ~ -(p1+p2+p3)*x1 + p5*x4,
9+
D(x2) ~ (p2*x1) + -(p4*x2),
10+
D(x3) ~ (p3*x1 + p4*x2),
11+
D(x4) ~ p1*x1 - p5*x4])
12+
13+
using DifferentialEquations: solve
14+
using Plots: plot
15+
prob = ODEProblem(structural_simplify(testcase), [x1 => 1, x2 => 0, x3 => 0, x4 => 0],
16+
(0.0, 120.0), [p1 => 0.02, p2 => 0.03, p3 => 0.04, p4 => 0.05, p5 => 0.06])
17+
sol = solve(prob)
18+
plot(sol, vars=[x1, x2, x3, x4])
19+
20+
new = toexpr(testcase.eqs[1].rhs)
21+
binarize!(new)
22+
factor!(new)
23+
24+
@named test2 = ODESystem([D(x1) ~ (p1+p2+p3),
25+
D(x2) ~ p1*x1,
26+
D(x3) ~ x1])
27+
t = toexpr(test2.eqs[1].rhs)
28+
factor!(toexpr(test2.eqs[1].rhs))
29+
30+
@parameters t1
31+
32+
33+
apply_transform(IntervalTransform(), testcase)
34+
35+
36+
# THIS IS MY TEST CASE, THESE SIX LINES
37+
using ModelingToolkit, OrdinaryDiffEq
38+
@parameters t
39+
@variables x(t) y(t)
40+
D = Differential(t)
41+
@named test = ODESystem([D(x) ~ exp(x+y^2)])
42+
apply_transform(IntervalTransform(), test)
43+
44+
45+
# THIS TEST DOES NOT WORK, BECAUSE THE VECTOR-VALUED PARAMETERS
46+
# AND VARIABLES MEAN I NEED TO PRE-EMPTIVELY MAKE/USE OVER/UNDERESTIMATES
47+
# (This is possible but would need a little bit of a rewrite/change. and
48+
# some thinking about how to handle cases like that)
49+
using ModelingToolkit, OrdinaryDiffEq
50+
@parameters t
51+
@variables x[1:2,1:3,1:4](t)
52+
D = Differential(t)
53+
@named test2 = ODESystem([D(x[1]) ~ x[1]*x[2]])
54+
apply_transform(IntervalTransform(), test2)
55+
56+
57+
# Trying out getting the thing in a variable form
58+
a = gensym(:aux)
59+
b = Symbol(string(a)[3:5] * string(a)[7:end])
60+
c = genvar(b)[1]
61+
62+
63+
64+
@named test3 = ODESystem([D(x1) ~ 5,
65+
D(x2) ~ x1*p1])
66+
apply_transform(IntervalTransform(), test3)
67+
68+
69+
70+
a = gensym(:test)
71+
# b = :(Symbolics._parse_vars(:variables, Real, $a))
72+
# eval(b)
73+
c = Symbolics._parse_vars(:variables, Real, a)
74+
75+
76+
using ModelingToolkit, OrdinaryDiffEq
77+
a = gensym(:var)
78+
@parameters t σ ρ β
79+
@variables x(t) y(t) z(t) $a(t)
80+
D = Differential(t)
81+
82+
eqs = [D(D(x)) ~ σ*(y-x),
83+
D(y) ~ x*(ρ-z)-y,
84+
D(z) ~ x*y - β*z]
85+
86+
eqs = [D(x) ~ σ+x]
87+
88+
@named sys = ODESystem(eqs)
89+
sys = ode_order_lowering(sys)
90+
91+
u0 = [D(x) => 2.0,
92+
x => 1.0,
93+
y => 0.0,
94+
z => 0.0]
95+
96+
p = [σ => 28.0,
97+
ρ => 10.0,
98+
β => 8/3]
99+
100+
tspan = (0.0,100.0)
101+
prob = ODEProblem(sys,u0,tspan,p,jac=true)
102+
103+
println(typeof(eqs))
104+
println(length(eqs))
105+
apply_transform(IntervalTransform(), sys)
106+
107+
108+
## New test!
109+
using ModelingToolkit, OrdinaryDiffEq
110+
@parameters t p
111+
@variables x(t)
112+
eqs = [D(x) ~ p[1]*x[1]^2 + p[1]]
113+
114+
# Try to do what @variables does, but with a symbol of choice
115+
Symbolics._parse_vars(:variables, Real, xs)
116+
# where xs is xs... is the names of the symbols to make into variables
117+
118+
a = gensym(:new)
119+
b = Symbol(string(a)[3:5] * string(a)[7:end])
120+
121+
genvar(b)
122+
genvar(b, [:t, :m])
123+
124+
125+
a = gensym(:var)
126+
println("The new variable a is $a")
127+
128+
for i = 1:1000
129+
gensym(:tes)
130+
end
131+
132+
133+
# function Base.iterate(a::Symbol)
134+
# return (a, nothing)
135+
# end
136+
# function Base.iterate(a::Symbol, ::Nothing)
137+
# return (a, nothing)
138+
# end
139+
Base.iterate(a::Symbol) = (a, nothing)
140+
Base.iterate(a::Symbol, state) = nothing
141+
142+
143+
# Original system from HDO Blackbox
144+
# dx[1] = -(p[1]+p[2]+p[3])*x[1] #Guaiacol
145+
# dx[2] = (p[2]*x[1]) + -(p[4]*x[2]) #Methoxycyclohexanone
146+
# dx[3] = (p[3]*x[1] + p[4]*x[2]) #Cyclohexanol
147+
# dx[4] = p[1]*x[1] #Fouled amount
148+
149+
# @variables t x(t) RHS(t) # independent and dependent variables
150+
# @parameters τ # parameters
151+
# D = Differential(t) # define an operator for the differentiation w.r.t. time
152+
153+
# # your first ODE, consisting of a single equation, indicated by ~
154+
# @named fol_separate = ODESystem([ RHS ~ (1 - x)/τ,
155+
# D(x) ~ RHS ])
156+
157+
158+
# using DifferentialEquations: solve
159+
# using Plots: plot
160+
161+
# prob = ODEProblem(structural_simplify(fol_separate), [x => 0.0], (0.0,10.0), [τ => 3.0])
162+
# sol = solve(prob)
163+
# plot(sol, vars=[x,RHS])
164+
165+
166+
xs = [1,2,3,4,5]
167+
function test!(vals::Vector{Int64})
168+
p = collect(1:length(vals))
169+
deleteat!(p, 3)
170+
push!(p, 3)
171+
vals[:] = vals[p]
172+
end

‎src/interval/interval.jl‎

Lines changed: 51 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -5,33 +5,23 @@ Rules for constructing interval bounding expressions
55
# Structure used to indicate an overload with intervals is preferable
66
struct IntervalTransform <: AbstractTransform end
77

8-
# Creates names for interval state variables [DEPRECATED]
9-
# function var_names(::IntervalTransform, s::Num)
10-
# if s.val.metadata.value[1]==:variables
11-
# arg_list = Symbol[]
12-
# for i in s.val.arguments
13-
# push!(arg_list, get_name(i))
14-
# end
15-
# sL = genvar(Symbol(string(s.val.f)*"_lo"), arg_list)
16-
# sU = genvar(Symbol(string(s.val.f)*"_hi"), arg_list)
17-
# elseif s.val.metadata.parent.value[1]==:parameters
18-
# sL = genparam(Symbol(string(s.val.name)*"_lo"))
19-
# sU = genparam(Symbol(string(s.val.name)*"_hi"))
20-
# else
21-
# error("Not a variable or a parameter, check type of s")
22-
# end
23-
# return sL, sU
24-
# end
25-
268
function var_names(::IntervalTransform, s::Term{Real, Base.ImmutableDict{DataType, Any}}) #The variables
279
arg_list = Symbol[]
28-
for i in s.arguments
29-
push!(arg_list, get_name(i))
10+
if haskey(s.metadata, ModelingToolkit.MTKParameterCtx)
11+
sL = genparam(Symbol(string(get_name(s))*"_lo"))
12+
sU = genparam(Symbol(string(get_name(s))*"_hi"))
13+
else
14+
for i in s.arguments
15+
push!(arg_list, get_name(i))
16+
end
17+
sL = genvar(Symbol(string(get_name(s))*"_lo"), arg_list)
18+
sU = genvar(Symbol(string(get_name(s))*"_hi"), arg_list)
3019
end
31-
sL = genvar(Symbol(string(s.f)*"_lo"), arg_list)
32-
sU = genvar(Symbol(string(s.f)*"_hi"), arg_list)
3320
return sL, sU
3421
end
22+
function var_names(::IntervalTransform, s::Real)
23+
return s, s
24+
end
3525
function var_names(::IntervalTransform, s::Term{Real, Nothing}) #Any terms like "Differential"
3626
if length(s.arguments)>1
3727
error("Multiple arguments not supported.")
@@ -56,11 +46,33 @@ function var_names(::IntervalTransform, s::Term{Real, Nothing}) #Any terms like
5646
return sL, sU
5747
end
5848
function var_names(::IntervalTransform, s::Sym) #The parameters
59-
sL = genparam(Symbol(string(s.name)*"_lo"))
60-
sU = genparam(Symbol(string(s.name)*"_hi"))
49+
sL = genparam(Symbol(string(get_name(s))*"_lo"))
50+
sU = genparam(Symbol(string(get_name(s))*"_hi"))
6151
return sL, sU
6252
end
6353

54+
function translate_initial_conditions(::IntervalTransform, prob::ODESystem, new_eqs::Vector{Equation})
55+
vars, params = extract_terms(new_eqs)
56+
var_defaults = Dict{Any, Any}()
57+
param_defaults = Dict{Any, Any}()
58+
59+
for (key, val) in prob.defaults
60+
name_lo = String(get_name(key))*"_"*"lo"
61+
name_hi = String(get_name(key))*"_"*"hi"
62+
for i in vars
63+
if in(String(i.f.name), (name_lo, name_hi))
64+
var_defaults[i] = val
65+
end
66+
end
67+
for i in params
68+
if in(String(i.name), (name_lo, name_hi))
69+
param_defaults[i] = val
70+
end
71+
end
72+
end
73+
return var_defaults, param_defaults
74+
end
75+
6476

6577

6678
# function var_names(::IntervalTransform, s::Number)
@@ -75,17 +87,26 @@ get_name(x::Sym{SymbolicUtils.FnType{Tuple{Any}, Real}, Nothing}) = x.name
7587
"""
7688
Takes x[1,1] returns :x_1_1
7789
"""
78-
function get_name(z::Term{SymbolicUtils.FnType{Tuple, Real}, Nothing})
79-
d = value(z).val.f.arguments
80-
x = string(d[1])
90+
function get_name(s::Term{SymbolicUtils.FnType{Tuple, Real}, Nothing})
91+
d = s.arguments
92+
new_var = string(d[1])
8193
for i in 2:length(d)
82-
x = x*"_"*string(i)
94+
new_var = new_var*"_"*string(d[i])
8395
end
84-
Symbol(x)
96+
return Symbol(new_var)
8597
end
8698

8799
function get_name(s::Term)
88-
return get_name(s.f)
100+
if haskey(s.metadata, ModelingToolkit.MTKParameterCtx)
101+
d = s.arguments
102+
new_param = string(d[1])
103+
for i in 2:length(d)
104+
new_param = new_param*"_"*string(d[i])
105+
end
106+
return Symbol(new_param)
107+
else
108+
return get_name(s.f)
109+
end
89110
end
90111
function get_name(s::Sym)
91112
return s.name

‎src/interval/rules.jl‎

Lines changed: 34 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -6,17 +6,17 @@ using IfElse
66
# to do correctly rounded stuff on GPUs but otherwise the operators should work out)
77

88
#=
9-
Unitary rules
9+
Unitary Rules
1010
=#
11-
function transform_rule(::IntervalTransform, ::Symbol, yL, yU, xL, xU)
11+
function transform_rule(::IntervalTransform, nothing, yL, yU, xL, xU)
1212
rl = Assignment(yL, xL)
1313
ru = Assignment(yU, xU)
1414
AssignmentPair(rl, ru)
1515
end
1616
function transform_rule(::IntervalTransform, ::typeof(exp), yL, yU, xL, xU)
1717
rl = Assignment(yL, exp(xL))
1818
ru = Assignment(yU, exp(xU))
19-
AssignmentPair(rl, ru)
19+
return rl, ru
2020
end
2121

2222
#=
@@ -25,31 +25,43 @@ Binary Rules
2525
function transform_rule(::IntervalTransform, ::typeof(+), zL, zU, xL, xU, yL, yU)
2626
rl = Assignment(zL, xL + yL)
2727
ru = Assignment(zU, xU + yU)
28-
AssignmentPair(rl, ru)
28+
return rl, ru
2929
end
3030
function transform_rule(::IntervalTransform, ::typeof(-), zL, zU, xL, xU, yL, yU)
3131
rl = Assignment(zL, xL - yU)
3232
ru = Assignment(zU, xU - yL)
33-
AssignmentPair(rl, ru)
33+
return rl, ru
3434
end
3535
function transform_rule(::IntervalTransform, ::typeof(*), zL, zU, xL, xU, yL, yU)
36-
rl = Assignment(zL, :(ifelse($yL >= 0.0,
37-
ifelse($xL >= 0.0, $xL*$yL,
38-
ifelse($xU <= 0.0, $xL*$yU, $xL*$yU)),
39-
ifelse($yU <= 0.0,
40-
ifelse($xL >= 0.0, $xU*$yL,
41-
ifelse($xU <= 0.0, $xU*$yU, $xU*$yL)),
42-
ifelse($xL > 0.0, $xU*$yL,
43-
ifelse($xU < 0.0, $xL*$yU, min($xU*$yL, $xU*$yU)))))))
44-
ru = Assignment(zU, :(ifelse($yL >= 0.0,
45-
ifelse($xL >= 0.0, $xU*$yU,
46-
ifelse($xU <= 0.0, $xU*$yL, $xU*$yU)),
47-
ifelse($yU <= 0.0,
48-
ifelse($xL >= 0.0, $xL*$yU,
49-
ifelse($xU <= 0.0, $xL*$yL, $xL*$yL)),
50-
ifelse($xL > 0.0, $xU*$yU,
51-
ifelse($xU < 0.0, $xL*$yL, max($xL*$yL, $xU*$yU)))))))
52-
AssignmentPair(rl, ru)
36+
rl = Assignment(zL, IfElse.ifelse(yL >= 0.0,
37+
IfElse.ifelse(xL >= 0.0, xL*yL,
38+
IfElse.ifelse(xU <= 0.0, xL*yU, xL*yU)),
39+
IfElse.ifelse(yU <= 0.0,
40+
IfElse.ifelse(xL >= 0.0, xU*yL,
41+
IfElse.ifelse(xU <= 0.0, xU*yU, xU*yL)),
42+
IfElse.ifelse(xL > 0.0, xU*yL,
43+
IfElse.ifelse(xU < 0.0, xL*yU, min(xU*yL, xU*yU))))))
44+
ru = Assignment(zU, IfElse.ifelse(yL >= 0.0,
45+
IfElse.ifelse(xL >= 0.0, xU*yU,
46+
IfElse.ifelse(xU <= 0.0, xU*yL, xU*yU)),
47+
IfElse.ifelse(yU <= 0.0,
48+
IfElse.ifelse(xL >= 0.0, xL*yU,
49+
IfElse.ifelse(xU <= 0.0, xL*yL, xL*yL)),
50+
IfElse.ifelse(xL > 0.0, xU*yU,
51+
IfElse.ifelse(xU < 0.0, xL*yL, max(xL*yL, xU*yU))))))
52+
return rl, ru
53+
end
54+
55+
56+
function transform_rule(::IntervalTransform, ::typeof(min), zL, zU, xL, xU, yL, yU)
57+
rl = Assignment(zL, min(xL, yL))
58+
ru = Assignment(zU, min(xU, yU))
59+
return rl, ru
60+
end
61+
function transform_rule(::IntervalTransform, ::typeof(max), zL, zU, xL, xU, yL, yU)
62+
rl = Assignment(zL, max(xL, yL))
63+
ru = Assignment(zU, max(xU, yU))
64+
return rl, ru
5365
end
5466

5567
# TODO: /, ^

0 commit comments

Comments
 (0)