Skip to content

refactor: format #3710

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open
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
9 changes: 5 additions & 4 deletions docs/src/basics/Events.md
Original file line number Diff line number Diff line change
Expand Up @@ -92,8 +92,8 @@ The basic purely symbolic continuous event interface to encode *one* continuous
event is

```julia
AbstractSystem(eqs, ...; continuous_events::Vector{Equation})
AbstractSystem(eqs, ...; continuous_events::Pair{Vector{Equation}, Vector{Equation}})
AbstractSystem(eqs, _...; continuous_events::Vector{Equation})
AbstractSystem(eqs, _...; continuous_events::Pair{Vector{Equation}, Vector{Equation}})
```

In the former, equations that evaluate to 0 will represent conditions that should
Expand Down Expand Up @@ -272,7 +272,7 @@ In addition to continuous events, discrete events are also supported. The
general interface to represent a collection of discrete events is

```julia
AbstractSystem(eqs, ...; discrete_events = [condition1 => affect1, condition2 => affect2])
AbstractSystem(eqs, _...; discrete_events = [condition1 => affect1, condition2 => affect2])
```

where conditions are symbolic expressions that should evaluate to `true` when an
Expand Down Expand Up @@ -497,7 +497,8 @@ so far we aren't using anything that's not possible with the implicit interface.
You can also write

```julia
[temp ~ furnace_off_threshold] => ModelingToolkit.ImperativeAffect(modified = (;
[temp ~
furnace_off_threshold] => ModelingToolkit.ImperativeAffect(modified = (;
furnace_on)) do x, o, i, c
@set! x.furnace_on = false
end
Expand Down
2 changes: 1 addition & 1 deletion docs/src/basics/FAQ.md
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ The same principle applies to any parameter type that is not `Float64`.
@parameters p1::Int # integer-valued
@parameters p2::Bool # boolean-valued
@parameters p3::MyCustomStructType # non-numeric
@parameters p4::ComponentArray{...} # non-standard array
@parameters p4::ComponentArray{_...} # non-standard array
```

## Getting the index for a symbol
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
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
2 changes: 1 addition & 1 deletion docs/src/tutorials/linear_analysis.md
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ This is signified by the name being the middle argument to `connect`.
Of the above mentioned functions, all except for [`open_loop`](@ref) return the output of [`ModelingToolkit.linearize`](@ref), which is

```julia
matrices, simplified_sys = linearize(...)
matrices, simplified_sys = linearize(_...)
# matrices = (; A, B, C, D)
```

Expand Down
30 changes: 20 additions & 10 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 All @@ -211,7 +219,7 @@ end
function MTK.get_V_values(model::CasADiModel)
value_getter = MTK.successful_solve(model) ? CasADi.debug_value : CasADi.value
(nu, nt) = size(model.V.u)
if nu*nt != 0
if nu * nt != 0
V_vals = value_getter(model.solver_opti, model.V.u)
size(V_vals, 2) == 1 && (V_vals = V_vals')
V_vals = [[V_vals[i, j] for i in 1:nu] for j in 1:nt]
Expand All @@ -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
63 changes: 42 additions & 21 deletions ext/MTKInfiniteOptExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -48,26 +48,32 @@ 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)
@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)
@constraint(m.model, expr.lhs - expr.rhs==0)
elseif expr.relational_op === Symbolics.geq
@constraint(m.model, expr.lhs - expr.rhs0)
@constraint(m.model, expr.lhs - expr.rhs0)
else
@constraint(m.model, expr.lhs - expr.rhs0)
@constraint(m.model, expr.lhs - expr.rhs0)
end
end
MTK.set_objective!(m::InfiniteOptModel, expr) = @objective(m.model, Min, expr)
Expand All @@ -76,20 +82,27 @@ 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 +138,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 +155,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 +173,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 +210,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 +244,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
13 changes: 7 additions & 6 deletions 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 Expand Up @@ -141,17 +142,17 @@ end
function MTK.lowered_integral(m::PyomoDynamicOptModel, arg, lo, hi)
@unpack model, model_sym, t_sym, dummy_sym = m
total = 0
dt = Pyomo.pyconvert(Float64, (model.t.at(-1) - model.t.at(1))/(model.steps - 1))
dt = Pyomo.pyconvert(Float64, (model.t.at(-1) - model.t.at(1)) / (model.steps - 1))
f = Symbolics.build_function(arg, model_sym, t_sym, expression = false)
for (i, t) in enumerate(model.t)
if Bool(lo < t) && Bool(t < hi)
t_p = model.t.at(i-1)
t_p = model.t.at(i - 1)
Δt = min(t - lo, t - t_p)
total += 0.5*Δt*(f(model, t) + f(model, t_p))
total += 0.5 * Δt * (f(model, t) + f(model, t_p))
elseif Bool(t >= hi) && Bool(t - dt < hi)
t_p = model.t.at(i-1)
t_p = model.t.at(i - 1)
Δt = hi - t + dt
total += 0.5*Δt*(f(model, t) + f(model, t_p))
total += 0.5 * Δt * (f(model, t) + f(model, t_p))
end
end
PyomoVar(model.tₛ * total)
Expand Down
Loading
Loading