From 5f37c34c689c66ce172737a99a5085f9915d55ff Mon Sep 17 00:00:00 2001 From: TT Date: Fri, 21 Jul 2023 14:34:18 +0200 Subject: [PATCH] adjustments for new FMIImport/FMICore versions --- Project.toml | 5 +- docs/src/faq.md | 37 ++++++++++-- src/FMI2/sim.jl | 122 +++++++++++++++++++++++++++++++++++++--- src/extensions/Plots.jl | 10 ++-- 4 files changed, 154 insertions(+), 20 deletions(-) diff --git a/Project.toml b/Project.toml index aae166a0..2d91a739 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "FMI" uuid = "14a09403-18e3-468f-ad8a-74f8dda2d9ac" authors = ["TT ", "LM ", "JK "] -version = "0.12.3" +version = "0.12.4" [deps] DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" @@ -9,6 +9,7 @@ DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" FMIExport = "31b88311-cab6-44ed-ba9c-fe5a9abbd67a" FMIImport = "9fcbc62e-52a0-44e9-a616-1359a0008194" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Requires = "ae029012-a4dd-5104-9daa-d747884805df" ThreadPools = "b189fb0b-2eb5-4ed4-bc0c-d34c51242431" @@ -17,7 +18,7 @@ ThreadPools = "b189fb0b-2eb5-4ed4-bc0c-d34c51242431" DiffEqCallbacks = "2.26.0" DifferentialEquations = "7.7.0" FMIExport = "0.2.0" -FMIImport = "0.15.6" +FMIImport = "0.15.7" ProgressMeter = "1.7.0" Requires = "1.3.0" ThreadPools = "2.1.1" diff --git a/docs/src/faq.md b/docs/src/faq.md index e4579c6a..0bef5d48 100644 --- a/docs/src/faq.md +++ b/docs/src/faq.md @@ -3,14 +3,16 @@ This list some common - often numerical - errors, that can be fixed by better understanding the ODE-Problem inside your FMU. +---------- - -## Solving non-linear system +## Solving non-linear system failes ### Description Error message or warning, that solving of a non-linear system failed, close to the simulation start time. ### Example -- `Solving non-linear system 101 failed at time=3e-05.` +``` +Solving non-linear system 101 failed at time=3e-05. +``` ### Reason This could be, because the first step of the integration is accepted by the solver's error estimation, but shouldn't. This is usually, if the first step is picked to large by the solver's start step size heuristics. @@ -18,7 +20,7 @@ This could be, because the first step of the integration is accepted by the solv ### Fix - Try a small start value for the integration with keyword `dt`. - +---------- ## Access denied ### Description @@ -35,4 +37,29 @@ Access denied This is because your OS doesn't allow to interact with the binaries shipped with *FMI.jl*. ### Fix -This can easily be solved by fixing the binary's permission options or is automatically fixed if Julia runs with admin privileges. \ No newline at end of file +This can easily be solved by fixing the binary's permission options or is automatically fixed if Julia runs with admin privileges. + +---------- + +## Double Callback Crossing +### Description +Error message, that solving failed because of `double callback crossing`. + +### Example +``` +ERROR: +Double callback crossing floating pointer reducer errored. Report this issue. +``` + +### Reason +This is because the event instant (time point) of an FMU event indicator can't be found precisely. + +### Fix +This can be solved by allowing for more interpolation points during searching of the zero crossing: +``` +fmu.executionConfig.rootSearchInterpolationPoints = 1000 # default value is 10 +``` +This will have negative performance impact on systems with extreme amount of events (thousands per second). +For systems with only a few events there won't be a noticable slow down. + +---------- \ No newline at end of file diff --git a/src/FMI2/sim.jl b/src/FMI2/sim.jl index 1ac0068c..49495692 100644 --- a/src/FMI2/sim.jl +++ b/src/FMI2/sim.jl @@ -15,7 +15,9 @@ import FMIImport: prepareSolveFMU, finishSolveFMU, handleEvents import FMIImport: undual using FMIImport.ChainRulesCore -import FMIImport.ForwardDiff + +import FMIImport.ReverseDiff +import LinearAlgebra: eigvals import ProgressMeter import ThreadPools @@ -168,6 +170,64 @@ function saveValues(c::FMU2Component, recordValues, x, t, integrator, inputFunct return (fmiGet(c, recordValues)...,) end +function saveEventIndicators(c::FMU2Component, recordEventIndicators, x, t, integrator, inputFunction, inputValues) + + @assert c.state == fmi2ComponentStateContinuousTimeMode "saveEventIndicators(...): Must be in continuous time mode." + + c.solution.evals_saveeventindicators += 1 + + #x_old = fmi2GetContinuousStates(c) + #t_old = c.t + + if !c.fmu.isZeroState + fmi2SetContinuousStates(c, x) + end + fmi2SetTime(c, t) + if inputFunction != nothing + fmi2SetReal(c, inputValues, inputFunction(c, x, t)) + end + + #fmi2SetContinuousStates(c, x_old) + #fmi2SetTime(c, t_old) + + out = zeros(fmi2Real, c.fmu.modelDescription.numberOfEventIndicators) + fmi2GetEventIndicators!(c, out) + + return (out[recordEventIndicators]...,) +end + +function saveEigenvalues(c::FMU2Component, x, t, integrator, inputFunction, inputValues) + + @assert c.state == fmi2ComponentStateContinuousTimeMode "saveEigenvalues(...): Must be in continuous time mode." + + c.solution.evals_saveeigenvalues += 1 + + #x_old = fmi2GetContinuousStates(c) + #t_old = c.t + + if !c.fmu.isZeroState + fmi2SetContinuousStates(c, x) + end + fmi2SetTime(c, t) + if inputFunction != nothing + fmi2SetReal(c, inputValues, inputFunction(c, x, t)) + end + + #fmi2SetContinuousStates(c, x_old) + #fmi2SetTime(c, t_old) + + A = ReverseDiff.jacobian(_x -> FMI.fx(c, _x, [], t), x) + eigs = eigvals(A) + + vals = [] + for e in eigs + push!(vals, real(e)) + push!(vals, imag(e)) + end + + return (vals...,) +end + function fx(c::FMU2Component, dx::AbstractArray{<:Real}, x::AbstractArray{<:Real}, @@ -246,6 +306,8 @@ Keywords: - solver: Any Julia-supported ODE-solver (default is the DifferentialEquations.jl default solver, currently `AutoTsit5(Rosenbrock23())`) - customFx: [deprecated] Ability to give a custom state derivative function ẋ=f(x,t) - recordValues: Array of variables (strings or variableIdentifiers) to record. Results are returned as `DiffEqCallbacks.SavedValues` + - recordEventIndicators: Array or Range of event indicators (identified by integers) to record + - recordEigenvalues: Boolean value, if eigenvalues shall be computed and recorded - saveat: Time points to save values at (interpolated) - setup: Boolean, if FMU should be setup (default=true) - reset: Union{Bool, :auto}, if FMU should be reset before simulation (default reset=:auto) @@ -264,6 +326,8 @@ function fmi2SimulateME(fmu::FMU2, c::Union{FMU2Component, Nothing}=nothing, tsp solver = nothing, customFx = nothing, recordValues::fmi2ValueReferenceFormat = nothing, + recordEventIndicators::Union{AbstractArray{<:Integer, 1}, UnitRange{<:Integer}, Nothing} = nothing, + recordEigenvalues::Bool=false, saveat = nothing, x0::Union{AbstractArray{<:Real}, Nothing} = nothing, setup::Union{Bool, Nothing} = nothing, @@ -275,7 +339,8 @@ function fmi2SimulateME(fmu::FMU2, c::Union{FMU2Component, Nothing}=nothing, tsp inputFunction = nothing, parameters::Union{Dict{<:Any, <:Any}, Nothing} = nothing, dtmax::Union{Real, Nothing} = nothing, - callbacks = [], + callbacksBefore = [], + callbacksAfter = [], showProgress::Bool = true, kwargs...) @@ -287,11 +352,11 @@ function fmi2SimulateME(fmu::FMU2, c::Union{FMU2Component, Nothing}=nothing, tsp if inputFunction != nothing if hasmethod(inputFunction, Tuple{fmi2Real}) _inputFunction = (c, u, t) -> inputFunction(t) - elseif hasmethod(inputFunction, Tuple{Union{FMU2Component, Nothing}, fmi2Real}) + elseif hasmethod(inputFunction, Tuple{Union{FMU2Component, Nothing}, fmi2Real}) _inputFunction = (c, u, t) -> inputFunction(c, t) - elseif hasmethod(inputFunction, Tuple{Union{FMU2Component, Nothing}, AbstractArray{fmi2Real,1}}) + elseif hasmethod(inputFunction, Tuple{Union{FMU2Component, Nothing}, AbstractArray{fmi2Real,1}}) _inputFunction = (c, u, t) -> inputFunction(c, u) - elseif hasmethod(inputFunction, Tuple{AbstractArray{fmi2Real,1}, fmi2Real}) + elseif hasmethod(inputFunction, Tuple{AbstractArray{fmi2Real,1}, fmi2Real}) _inputFunction = (c, u, t) -> inputFunction(u, t) else _inputFunction = inputFunction @@ -303,6 +368,7 @@ function fmi2SimulateME(fmu::FMU2, c::Union{FMU2Component, Nothing}=nothing, tsp inputValueReferences = prepareValueReference(fmu, inputValueReferences) savingValues = (length(recordValues) > 0) + savingEventIndicators = !isnothing(recordEventIndicators) hasInputs = (length(inputValueReferences) > 0) hasParameters = (parameters !== nothing) hasStartState = (x0 !== nothing) @@ -311,7 +377,7 @@ function fmi2SimulateME(fmu::FMU2, c::Union{FMU2Component, Nothing}=nothing, tsp cbs = [] - for cb in callbacks + for cb in callbacksBefore push!(cbs, cb) end @@ -416,11 +482,12 @@ function fmi2SimulateME(fmu::FMU2, c::Union{FMU2Component, Nothing}=nothing, tsp push!(cbs, stepCb) end - if savingValues + if savingValues dtypes = collect(fmi2DataTypeForValueReference(c.fmu.modelDescription, vr) for vr in recordValues) - fmusol.values = SavedValues(Float64, Tuple{dtypes...}) + fmusol.values = SavedValues(fmi2Real, Tuple{dtypes...}) fmusol.valueReferences = copy(recordValues) + savingCB = nothing if saveat === nothing savingCB = SavingCallback((u,t,integrator) -> saveValues(c, recordValues, u, t, integrator, _inputFunction, inputValueReferences), fmusol.values) @@ -433,6 +500,45 @@ function fmi2SimulateME(fmu::FMU2, c::Union{FMU2Component, Nothing}=nothing, tsp push!(cbs, savingCB) end + if savingEventIndicators + dtypes = collect(fmi2Real for ei in recordEventIndicators) + fmusol.eventIndicators = SavedValues(fmi2Real, Tuple{dtypes...}) + fmusol.recordEventIndicators = copy(recordEventIndicators) + + savingCB = nothing + if saveat === nothing + savingCB = SavingCallback((u,t,integrator) -> saveEventIndicators(c, recordEventIndicators, u, t, integrator, _inputFunction, inputValueReferences), + fmusol.eventIndicators) + else + savingCB = SavingCallback((u,t,integrator) -> saveEventIndicators(c, recordEventIndicators, u, t, integrator, _inputFunction, inputValueReferences), + fmusol.eventIndicators, + saveat=saveat) + end + + push!(cbs, savingCB) + end + + if recordEigenvalues + dtypes = collect(Float64 for _ in 1:2*length(c.fmu.modelDescription.stateValueReferences)) + fmusol.eigenvalues = SavedValues(fmi2Real, Tuple{dtypes...}) + + savingCB = nothing + if saveat === nothing + savingCB = SavingCallback((u,t,integrator) -> saveEigenvalues(c, u, t, integrator, _inputFunction, inputValueReferences), + fmusol.eigenvalues) + else + savingCB = SavingCallback((u,t,integrator) -> saveEigenvalues(c, u, t, integrator, _inputFunction, inputValueReferences), + fmusol.eigenvalues, + saveat=saveat) + end + + push!(cbs, savingCB) + end + + for cb in callbacksAfter + push!(cbs, cb) + end + # if auto_dt == true # @assert solver !== nothing "fmi2SimulateME(...): `auto_dt=true` but no solver specified, this is not allowed." # tmpIntegrator = init(c.problem, solver) diff --git a/src/extensions/Plots.jl b/src/extensions/Plots.jl index bfc0e2c1..8e970255 100644 --- a/src/extensions/Plots.jl +++ b/src/extensions/Plots.jl @@ -41,7 +41,7 @@ function fmiPlot(solution::FMUSolution; kwargs...) fmiPlot!(fig, solution; kwargs...) return fig end -function fmiPlot!(fig, solution::FMUSolution; +function fmiPlot!(fig::Plots.Plot, solution::FMUSolution; states::Union{Bool, Nothing}=nothing, values::Union{Bool, Nothing}=nothing, stateEvents::Union{Bool, Nothing}=nothing, @@ -209,9 +209,9 @@ Extended the original plot-command by plotting FMUs. For further information seek `?fmiPlot`. """ -function Plots.plot(solution::FMUSolution, args...; kwargs...) - fmiPlot(solution, args...; kwargs...) +function Plots.plot(solution::FMUSolution; kwargs...) + fmiPlot(solution; kwargs...) end -function Plots.plot!(fig, solution::FMUSolution, args...; kwargs...) - fmiPlot!(fig, solution, args...; kwargs...) +function Plots.plot!(fig::Plots.Plot, solution::FMUSolution; kwargs...) + fmiPlot!(fig, solution; kwargs...) end