Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion docs/src/API/variables.md
Original file line number Diff line number Diff line change
Expand Up @@ -293,7 +293,7 @@ For systems that contain parameters with metadata like described above, have som
In the example below, we define a system with tunable parameters and extract bounds vectors

```@example metadata
@variables x(t)=0 u(t)=0 [input = true] y(t)=0 [output = true]
@variables x(t)=0 u(t)=0 [input=true] y(t)=0 [output=true]
@parameters T [tunable = true, bounds = (0, Inf)]
@parameters k [tunable = true, bounds = (0, Inf)]
eqs = [D(x) ~ (-x + k * u) / T # A first-order system with time constant T and gain k
Expand Down
3 changes: 2 additions & 1 deletion docs/src/basics/MTKLanguage.md
Original file line number Diff line number Diff line change
Expand Up @@ -381,7 +381,8 @@ Refer the following example for different ways to define symbolic arrays.
@parameters begin
p1[1:4]
p2[1:N]
p3[1:N, 1:M] = 10,
p3[1:N,
1:M] = 10,
[description = "A multi-dimensional array of arbitrary length with description"]
(p4[1:N, 1:M] = 10),
[description = "An alternate syntax for p3 to match the syntax of vanilla parameters macro"]
Expand Down
2 changes: 1 addition & 1 deletion docs/src/basics/Validation.md
Original file line number Diff line number Diff line change
Expand Up @@ -108,7 +108,7 @@ function ModelingToolkit.get_unit(op::typeof(dummycomplex), args)
end

sts = @variables a(t)=0 [unit = u"cm"]
ps = @parameters s=-1 [unit = u"cm"] c=c [unit = u"cm"]
ps = @parameters s=-1 [unit=u"cm"] c=c [unit=u"cm"]
eqs = [D(a) ~ dummycomplex(c, s);]
sys = System(
eqs, t, [sts...;], [ps...;], name = :sys, checks = ~ModelingToolkit.CheckUnits)
Expand Down
3 changes: 2 additions & 1 deletion docs/src/examples/perturbation.md
Original file line number Diff line number Diff line change
Expand Up @@ -99,7 +99,8 @@ sol = solve(prob)
plot(sol, idxs = substitute(x_series, ϵ => 0.1); label = "Perturbative (ϵ=0.1)")

x_exact(t, ϵ) = exp(-ϵ * t) * sin(√(1 - ϵ^2) * t) / √(1 - ϵ^2)
@assert isapprox(sol(π/2; idxs = substitute(x_series, ϵ => 0.1)), x_exact(π/2, 0.1); atol = 1e-2) # compare around 1st peak # hide
@assert isapprox(
sol(π/2; idxs = substitute(x_series, ϵ => 0.1)), x_exact(π/2, 0.1); atol = 1e-2) # compare around 1st peak # hide
plot!(sol.t, x_exact.(sol.t, 0.1); label = "Exact (ϵ=0.1)")
```

Expand Down
21 changes: 13 additions & 8 deletions docs/src/examples/sparse_jacobians.md
Original file line number Diff line number Diff line change
Expand Up @@ -23,15 +23,20 @@ function brusselator_2d_loop(du, u, p, t)
@inbounds for I in CartesianIndices((N, N))
i, j = Tuple(I)
x, y = xyd_brusselator[I[1]], xyd_brusselator[I[2]]
ip1, im1, jp1, jm1 = limit(i + 1, N), limit(i - 1, N), limit(j + 1, N),
ip1, im1, jp1,
jm1 = limit(i + 1, N), limit(i - 1, N), limit(j + 1, N),
limit(j - 1, N)
du[i, j, 1] = alpha * (u[im1, j, 1] + u[ip1, j, 1] + u[i, jp1, 1] + u[i, jm1, 1] -
4u[i, j, 1]) +
B + u[i, j, 1]^2 * u[i, j, 2] - (A + 1) * u[i, j, 1] +
brusselator_f(x, y, t)
du[i, j, 2] = alpha * (u[im1, j, 2] + u[ip1, j, 2] + u[i, jp1, 2] + u[i, jm1, 2] -
4u[i, j, 2]) +
A * u[i, j, 1] - u[i, j, 1]^2 * u[i, j, 2]
du[i,
j,
1] = alpha * (u[im1, j, 1] + u[ip1, j, 1] + u[i, jp1, 1] + u[i, jm1, 1] -
4u[i, j, 1]) +
B + u[i, j, 1]^2 * u[i, j, 2] - (A + 1) * u[i, j, 1] +
brusselator_f(x, y, t)
du[i,
j,
2] = alpha * (u[im1, j, 2] + u[ip1, j, 2] + u[i, jp1, 2] + u[i, jm1, 2] -
4u[i, j, 2]) +
A * u[i, j, 1] - u[i, j, 1]^2 * u[i, j, 2]
end
end
p = (3.4, 1.0, 10.0, step(xyd_brusselator))
Expand Down
4 changes: 2 additions & 2 deletions docs/src/examples/tearing_parallelism.md
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D

# Basic electric components
@connector function Pin(; name)
@variables v(t)=1.0 i(t)=1.0 [connect = Flow]
@variables v(t)=1.0 i(t)=1.0 [connect=Flow]
System(Equation[], t, [v, i], [], name = name)
end

Expand All @@ -36,7 +36,7 @@ function ConstantVoltage(; name, V = 1.0)
end

@connector function HeatPort(; name)
@variables T(t)=293.15 Q_flow(t)=0.0 [connect = Flow]
@variables T(t)=293.15 Q_flow(t)=0.0 [connect=Flow]
System(Equation[], t, [T, Q_flow], [], name = name)
end

Expand Down
4 changes: 3 additions & 1 deletion docs/src/tutorials/disturbance_modeling.md
Original file line number Diff line number Diff line change
Expand Up @@ -188,7 +188,9 @@ disturbance_inputs = [ssys.d1, ssys.d2]
P = ssys.system_model
outputs = [P.inertia1.phi, P.inertia2.phi, P.inertia1.w, P.inertia2.w]

(f_oop, f_ip), x_sym, p_sym, io_sys = ModelingToolkit.generate_control_function(
(f_oop, f_ip), x_sym,
p_sym,
io_sys = ModelingToolkit.generate_control_function(
model_with_disturbance, inputs, disturbance_inputs; disturbance_argument = true)

g = ModelingToolkit.build_explicit_observed_function(
Expand Down
28 changes: 19 additions & 9 deletions ext/MTKCasADiDynamicOptExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,9 @@ struct MXLinearInterpolation
t::Vector{Float64}
dt::Float64
end
Base.getindex(m::MXLinearInterpolation, i...) = length(i) == length(size(m.u)) ? m.u[i...] : m.u[i..., :]
function Base.getindex(m::MXLinearInterpolation, i...)
length(i) == length(size(m.u)) ? m.u[i...] : m.u[i..., :]
end

mutable struct CasADiModel
model::Opti
Expand Down Expand Up @@ -55,7 +57,7 @@ function (M::MXLinearInterpolation)(τ)
(i > length(M.t) || i < 1) && error("Cannot extrapolate past the tspan.")
colons = ntuple(_ -> (:), length(size(M.u)) - 1)
if i < length(M.t)
M.u[colons..., i] + Δ*(M.u[colons..., i+1] - M.u[colons..., i])
M.u[colons..., i] + Δ*(M.u[colons..., i + 1] - M.u[colons..., i])
else
M.u[colons..., i]
end
Expand All @@ -65,7 +67,9 @@ function MTK.CasADiDynamicOptProblem(sys::System, op, tspan;
dt = nothing,
steps = nothing,
guesses = Dict(), kwargs...)
prob, _ = MTK.process_DynamicOptProblem(CasADiDynamicOptProblem, CasADiModel, sys, op, tspan; dt, steps, guesses, kwargs...)
prob,
_ = MTK.process_DynamicOptProblem(
CasADiDynamicOptProblem, CasADiModel, sys, op, tspan; dt, steps, guesses, kwargs...)
prob
end

Expand Down Expand Up @@ -127,10 +131,10 @@ function MTK.lowered_integral(model::CasADiModel, expr, lo, hi)
for (i, t) in enumerate(model.U.t)
if lo < t < hi
Δt = min(dt, t - lo)
total += (0.5*Δt*(expr[i] + expr[i-1]))
total += (0.5*Δt*(expr[i] + expr[i - 1]))
elseif t >= hi && (t - dt < hi)
Δt = hi - t + dt
total += (0.5*Δt*(expr[i] + expr[i-1]))
total += (0.5*Δt*(expr[i] + expr[i - 1]))
end
end
model.tₛ * total
Expand Down Expand Up @@ -186,9 +190,13 @@ struct CasADiCollocation <: AbstractCollocation
tableau::DiffEqBase.ODERKTableau
end

MTK.CasADiCollocation(solver, tableau = MTK.constructDefault()) = CasADiCollocation(solver, tableau)
function MTK.CasADiCollocation(solver, tableau = MTK.constructDefault())
CasADiCollocation(solver, tableau)
end

function MTK.prepare_and_optimize!(prob::CasADiDynamicOptProblem, solver::CasADiCollocation; verbose = false, solver_options = Dict(), plugin_options = Dict(), kwargs...)
function MTK.prepare_and_optimize!(
prob::CasADiDynamicOptProblem, solver::CasADiCollocation; verbose = false,
solver_options = Dict(), plugin_options = Dict(), kwargs...)
solver_opti = add_solve_constraints!(prob, solver.tableau)
verbose || (solver_options["print_level"] = 0)
solver!(solver_opti, "$(solver.solver)", plugin_options, solver_options)
Expand Down Expand Up @@ -224,9 +232,11 @@ function MTK.get_t_values(model::CasADiModel)
value_getter = MTK.successful_solve(model) ? CasADi.debug_value : CasADi.value
ts = value_getter(model.solver_opti, model.tₛ) .* model.U.t
end
MTK.objective_value(model::CasADiModel) = CasADi.pyconvert(Float64, model.solver_opti.py.value(model.solver_opti.py.f))
function MTK.objective_value(model::CasADiModel)
CasADi.pyconvert(Float64, model.solver_opti.py.value(model.solver_opti.py.f))
end

function MTK.successful_solve(m::CasADiModel)
function MTK.successful_solve(m::CasADiModel)
isnothing(m.solver_opti) && return false
retcode = CasADi.return_status(m.solver_opti)
retcode == "Solve_Succeeded" || retcode == "Solved_To_Acceptable_Level"
Expand Down
54 changes: 37 additions & 17 deletions ext/MTKInfiniteOptExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -48,20 +48,26 @@ struct InfiniteOptDynamicOptProblem{uType, tType, isinplace, P, F, K} <:
end

MTK.generate_internal_model(m::Type{InfiniteOptModel}) = InfiniteModel()
MTK.generate_time_variable!(m::InfiniteModel, tspan, tsteps) = @infinite_parameter(m, t in [tspan[1], tspan[2]], num_supports = length(tsteps))
MTK.generate_state_variable!(m::InfiniteModel, u0::Vector, ns, ts) = @variable(m, U[i = 1:ns], Infinite(m[:t]), start=u0[i])
MTK.generate_input_variable!(m::InfiniteModel, c0, nc, ts) = @variable(m, V[i = 1:nc], Infinite(m[:t]), start=c0[i])
function MTK.generate_time_variable!(m::InfiniteModel, tspan, tsteps)
@infinite_parameter(m, t in [tspan[1], tspan[2]], num_supports = length(tsteps))
end
function MTK.generate_state_variable!(m::InfiniteModel, u0::Vector, ns, ts)
@variable(m, U[i = 1:ns], Infinite(m[:t]), start=u0[i])
end
function MTK.generate_input_variable!(m::InfiniteModel, c0, nc, ts)
@variable(m, V[i = 1:nc], Infinite(m[:t]), start=c0[i])
end

function MTK.generate_timescale!(m::InfiniteModel, guess, is_free_t)
@variable(m, tₛ ≥ 0, start = guess)
if !is_free_t
fix(tₛ, 1, force=true)
fix(tₛ, 1, force = true)
set_start_value(tₛ, 1)
end
tₛ
end

function MTK.add_constraint!(m::InfiniteOptModel, expr::Union{Equation, Inequality})
function MTK.add_constraint!(m::InfiniteOptModel, expr::Union{Equation, Inequality})
if expr isa Equation
@constraint(m.model, expr.lhs - expr.rhs == 0)
elseif expr.relational_op === Symbolics.geq
Expand All @@ -76,20 +82,26 @@ function MTK.JuMPDynamicOptProblem(sys::System, op, tspan;
dt = nothing,
steps = nothing,
guesses = Dict(), kwargs...)
prob, _ = MTK.process_DynamicOptProblem(JuMPDynamicOptProblem, InfiniteOptModel, sys, op, tspan; dt, steps, guesses, kwargs...)
prob,
_ = MTK.process_DynamicOptProblem(JuMPDynamicOptProblem, InfiniteOptModel, sys,
op, tspan; dt, steps, guesses, kwargs...)
prob
end

function MTK.InfiniteOptDynamicOptProblem(sys::System, op, tspan;
dt = nothing,
steps = nothing,
guesses = Dict(), kwargs...)
prob, pmap = MTK.process_DynamicOptProblem(InfiniteOptDynamicOptProblem, InfiniteOptModel, sys, op, tspan; dt, steps, guesses, kwargs...)
prob,
pmap = MTK.process_DynamicOptProblem(InfiniteOptDynamicOptProblem, InfiniteOptModel,
sys, op, tspan; dt, steps, guesses, kwargs...)
MTK.add_equational_constraints!(prob.wrapped_model, sys, pmap, tspan)
prob
end

MTK.lowered_integral(model::InfiniteOptModel, expr, lo, hi) = model.tₛ * InfiniteOpt.∫(expr, model.model[:t], lo, hi)
function MTK.lowered_integral(model::InfiniteOptModel, expr, lo, hi)
model.tₛ * InfiniteOpt.∫(expr, model.model[:t], lo, hi)
end
MTK.lowered_derivative(model::InfiniteOptModel, i) = ∂(model.U[i], model.model[:t])

function MTK.process_integral_bounds(model::InfiniteOptModel, integral_span, tspan)
Expand Down Expand Up @@ -125,7 +137,7 @@ function add_solve_constraints!(prob::JuMPDynamicOptProblem, tableau)
nᵥ = length(V)
if MTK.is_explicit(tableau)
K = Any[]
for τ in tsteps[1:end-1]
for τ in tsteps[1:(end - 1)]
for (i, h) in enumerate(c)
ΔU = sum([A[i, j] * K[j] for j in 1:(i - 1)], init = zeros(nᵤ))
Uₙ = [U[i](τ) + ΔU[i] * dt for i in 1:nᵤ]
Expand All @@ -142,14 +154,15 @@ function add_solve_constraints!(prob::JuMPDynamicOptProblem, tableau)
K = @variable(model, K[1:length(α), 1:nᵤ], Infinite(model[:t]))
ΔUs = A * K
ΔU_tot = dt * (K' * α)
for τ in tsteps[1:end-1]
for τ in tsteps[1:(end - 1)]
for (i, h) in enumerate(c)
ΔU = @view ΔUs[i, :]
Uₙ = U + ΔU * dt
@constraint(model, [j = 1:nᵤ], K[i, j]==(tₛ * f(Uₙ, V, p, τ + h * dt)[j]),
DomainRestrictions(t => τ), base_name="solve_K$i($τ)")
end
@constraint(model, [n = 1:nᵤ], U[n](τ) + ΔU_tot[n]==U[n](min(τ + dt, tsteps[end])),
@constraint(model,
[n = 1:nᵤ], U[n](τ) + ΔU_tot[n]==U[n](min(τ + dt, tsteps[end])),
DomainRestrictions(t => τ), base_name="solve_U($τ)")
end
end
Expand All @@ -159,15 +172,21 @@ struct JuMPCollocation <: AbstractCollocation
solver::Any
tableau::DiffEqBase.ODERKTableau
end
MTK.JuMPCollocation(solver, tableau = MTK.constructDefault()) = JuMPCollocation(solver, tableau)
function MTK.JuMPCollocation(solver, tableau = MTK.constructDefault())
JuMPCollocation(solver, tableau)
end

struct InfiniteOptCollocation <: AbstractCollocation
solver::Any
derivative_method::InfiniteOpt.AbstractDerivativeMethod
end
MTK.InfiniteOptCollocation(solver, derivative_method = InfiniteOpt.FiniteDifference(InfiniteOpt.Backward())) = InfiniteOptCollocation(solver, derivative_method)
function MTK.InfiniteOptCollocation(
solver, derivative_method = InfiniteOpt.FiniteDifference(InfiniteOpt.Backward()))
InfiniteOptCollocation(solver, derivative_method)
end

function MTK.prepare_and_optimize!(prob::JuMPDynamicOptProblem, solver::JuMPCollocation; verbose = false, kwargs...)
function MTK.prepare_and_optimize!(
prob::JuMPDynamicOptProblem, solver::JuMPCollocation; verbose = false, kwargs...)
model = prob.wrapped_model.model
verbose || set_silent(model)
# Unregister current solver constraints
Expand All @@ -190,7 +209,8 @@ function MTK.prepare_and_optimize!(prob::JuMPDynamicOptProblem, solver::JuMPColl
model
end

function MTK.prepare_and_optimize!(prob::InfiniteOptDynamicOptProblem, solver::InfiniteOptCollocation; verbose = false, kwargs...)
function MTK.prepare_and_optimize!(prob::InfiniteOptDynamicOptProblem,
solver::InfiniteOptCollocation; verbose = false, kwargs...)
model = prob.wrapped_model.model
verbose || set_silent(model)
set_derivative_method(model[:t], solver.derivative_method)
Expand Down Expand Up @@ -223,8 +243,8 @@ function MTK.successful_solve(model::InfiniteModel)
error("Model not solvable; please report this to github.com/SciML/ModelingToolkit.jl with a MWE.")

pstatus === FEASIBLE_POINT &&
(tstatus === OPTIMAL || tstatus === LOCALLY_SOLVED || tstatus === ALMOST_OPTIMAL ||
tstatus === ALMOST_LOCALLY_SOLVED)
(tstatus === OPTIMAL || tstatus === LOCALLY_SOLVED || tstatus === ALMOST_OPTIMAL ||
tstatus === ALMOST_LOCALLY_SOLVED)
end

import InfiniteOpt: JuMP, GeneralVariableRef
Expand Down
3 changes: 2 additions & 1 deletion ext/MTKPyomoDynamicOptExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -108,7 +108,8 @@ function MTK.add_constraint!(pmodel::PyomoDynamicOptModel, cons; n_idxs = 1)
else
cons.lhs - cons.rhs ≤ 0
end
expr = Symbolics.substitute(Symbolics.unwrap(expr), SPECIAL_FUNCTIONS_DICT, fold = false)
expr = Symbolics.substitute(
Symbolics.unwrap(expr), SPECIAL_FUNCTIONS_DICT, fold = false)

cons_sym = Symbol("cons", hash(cons))
if occursin(Symbolics.unwrap(t_sym), expr)
Expand Down
4 changes: 3 additions & 1 deletion src/bipartite_graph.jl
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,8 @@ Base.iterate(u::Unassigned, state) = nothing

Base.show(io::IO, ::Unassigned) = printstyled(io, "u"; color = :light_black)

struct Matching{U, V <: AbstractVector} <: AbstractVector{Union{U, Int}} #=> :Unassigned =#
#U=> :Unassigned =#
struct Matching{U, V <: AbstractVector} <: AbstractVector{Union{U, Int}}
match::V
inv_match::Union{Nothing, V}
end
Expand Down Expand Up @@ -634,6 +635,7 @@ function Graphs.incidence_matrix(g::BipartiteGraph, val = true)
I = Int[]
J = Int[]
for i in 𝑠vertices(g), n in 𝑠neighbors(g, i)

push!(I, i)
push!(J, n)
end
Expand Down
5 changes: 4 additions & 1 deletion src/deprecations.jl
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,7 @@ for T in [:ODEProblem, :DDEProblem, :SDEProblem, :SDDEProblem, :DAEProblem,
end

for pType in [SciMLBase.NullParameters, Nothing], uType in [Any, Nothing]

@eval function SciMLBase.$T(sys::System, u0::$uType, tspan, p::$pType; kw...)
ctor = string($T)
pT = string($(QuoteNode(pType)))
Expand Down Expand Up @@ -142,6 +143,7 @@ for T in [:NonlinearProblem, :NonlinearLeastSquaresProblem,
end
end
for pType in [SciMLBase.NullParameters, Nothing], uType in [Any, Nothing]

@eval function SciMLBase.$T(sys::System, u0::$uType, p::$pType; kw...)
ctor = string($T)
pT = string($(QuoteNode(pType)))
Expand Down Expand Up @@ -173,7 +175,8 @@ end

macro brownian(xs...)
return quote
Base.depwarn("`@brownian` is deprecated. Use `@brownians` instead", :brownian_macro)
Base.depwarn(
"`@brownian` is deprecated. Use `@brownians` instead", :brownian_macro)
$(@__MODULE__).@brownians $(xs...)
end |> esc
end
4 changes: 2 additions & 2 deletions src/inputoutput.jl
Original file line number Diff line number Diff line change
Expand Up @@ -120,7 +120,7 @@ function same_or_inner_namespace(u, var)
nu == nv || # namespaces are the same
startswith(nv, nu) || # or nv starts with nu, i.e., nv is an inner namespace to nu
occursin(NAMESPACE_SEPARATOR, string(getname(var))) &&
!occursin(NAMESPACE_SEPARATOR, string(getname(u))) # or u is top level but var is internal
!occursin(NAMESPACE_SEPARATOR, string(getname(u))) # or u is top level but var is internal
end

function inner_namespace(u, var)
Expand All @@ -129,7 +129,7 @@ function inner_namespace(u, var)
nu == nv && return false
startswith(nv, nu) || # or nv starts with nu, i.e., nv is an inner namespace to nu
occursin(NAMESPACE_SEPARATOR, string(getname(var))) &&
!occursin(NAMESPACE_SEPARATOR, string(getname(u))) # or u is top level but var is internal
!occursin(NAMESPACE_SEPARATOR, string(getname(u))) # or u is top level but var is internal
end

"""
Expand Down
6 changes: 4 additions & 2 deletions src/linearization.jl
Original file line number Diff line number Diff line change
Expand Up @@ -298,7 +298,8 @@ function (linfun::LinearizationFunction)(u, p, t)
error("Number of unknown variables ($(linfun.num_states)) does not match the number of input unknowns ($(length(u)))")
integ_cache = (linfun.caches,)
integ = MockIntegrator{true}(u, p, t, fun, integ_cache, nothing)
u, p, success = SciMLBase.get_initial_values(
u, p,
success = SciMLBase.get_initial_values(
linfun.prob, integ, fun, linfun.initializealg, Val(true);
linfun.initialize_kwargs...)
if !success
Expand Down Expand Up @@ -750,7 +751,8 @@ function linearize(sys, inputs, outputs; op = Dict(), t = 0.0,
allow_input_derivatives = false,
zero_dummy_der = false,
kwargs...)
lin_fun, ssys = linearization_function(sys,
lin_fun,
ssys = linearization_function(sys,
inputs,
outputs;
zero_dummy_der,
Expand Down
Loading
Loading