From 815138d4bfac3092686711c512f418d39e039c42 Mon Sep 17 00:00:00 2001 From: ThummeTo <83663542+ThummeTo@users.noreply.github.com> Date: Mon, 3 Jul 2023 08:48:54 +0200 Subject: [PATCH 1/3] v0.15.6 (#91) * FMU parameter sensitivities during optimization * minor bug fix * changed default initialization values * support for zero state FMUs * zero state FMU simulation * minor adjustments * PrepareValue AbstractVector fix (see #92) (#93) * added fmi2GetSolutionDerivative to fmi2_library.md (#90) * PrepareValue AbstractVector fix (see #92) --------- Co-authored-by: adribrune <94375638+adribrune@users.noreply.github.com> * minor fix --------- Co-authored-by: CasBex <123587431+CasBex@users.noreply.github.com> Co-authored-by: adribrune <94375638+adribrune@users.noreply.github.com> --- Project.toml | 4 +- src/FMI2/c.jl | 7 +- src/FMI2/convert.jl | 19 + src/FMI2/ext.jl | 7 +- src/FMI2/prep.jl | 10 +- src/FMI2/sens.jl | 718 ++++++++++++++++--------------------- src/FMIImport.jl | 13 +- test/FMI2/sensitivities.jl | 49 ++- 8 files changed, 403 insertions(+), 424 deletions(-) diff --git a/Project.toml b/Project.toml index 6f59c27..9d2555c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "FMIImport" uuid = "9fcbc62e-52a0-44e9-a616-1359a0008194" authors = ["TT ", "LM ", "JK "] -version = "0.15.5" +version = "0.15.6" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" @@ -16,7 +16,7 @@ ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" [compat] ChainRulesCore = "1.15.0" EzXML = "1.1.0" -FMICore = "0.17.0" +FMICore = "0.17.2" ForwardDiffChainRules = "0.1.1" SciMLSensitivity = "7.27.0" ZipFile = "0.10.0" diff --git a/src/FMI2/c.jl b/src/FMI2/c.jl index 68d83ad..5215672 100644 --- a/src/FMI2/c.jl +++ b/src/FMI2/c.jl @@ -599,7 +599,7 @@ function fmi2SetReal(c::FMU2Component, if track if status == fmi2StatusOK - for j in (c.A, c.B, c.C, c.D) + for j in (c.A, c.B, c.C, c.D, c.E, c.F) if any(collect(v in j.∂f_refs for v in vr)) FMICore.invalidate!(j) end @@ -1562,6 +1562,8 @@ function fmi2SetTime(c::FMU2Component, time::fmi2Real; soft::Bool=false, track:: FMICore.invalidate!(c.B) FMICore.invalidate!(c.C) FMICore.invalidate!(c.D) + FMICore.invalidate!(c.E) + FMICore.invalidate!(c.F) end end @@ -1620,6 +1622,9 @@ function fmi2SetContinuousStates(c::FMU2Component, FMICore.invalidate!(c.A) FMICore.invalidate!(c.C) + + FMICore.invalidate!(c.E) + FMICore.invalidate!(c.F) end end diff --git a/src/FMI2/convert.jl b/src/FMI2/convert.jl index b64be77..51c013b 100644 --- a/src/FMI2/convert.jl +++ b/src/FMI2/convert.jl @@ -107,6 +107,25 @@ function fmi2ModelVariablesForValueReference(md::fmi2ModelDescription, vr::fmi2V ar end +""" +ToDo. +""" +function fmi2DataTypeForValueReference(md::fmi2ModelDescription, vr::fmi2ValueReference) + mv = fmi2ModelVariablesForValueReference(md, vr)[1] + if !isnothing(mv.Real) + return fmi2Real + elseif !isnothing(mv.Integer) || !isnothing(mv.Enumeration) + return fmi2Integer + elseif !isnothing(mv.Boolean) + return fmi2Boolean + elseif !isnothing(mv.String) + return fmi2String + else + @assert false "fmi2TypeForValueReference(...): Unknown data type for value reference `$(vr)`." + end + return nothing +end + """ fmi2StringToValueReference(md::fmi2ModelDescription, name::String) diff --git a/src/FMI2/ext.jl b/src/FMI2/ext.jl index 565c533..728341c 100644 --- a/src/FMI2/ext.jl +++ b/src/FMI2/ext.jl @@ -149,6 +149,7 @@ function fmi2Load(pathToFMU::String; unpackPath::Union{String, Nothing}=nothing, # parse modelDescription.xml fmu.modelDescription = fmi2LoadModelDescription(pathToModelDescription) fmu.modelName = fmu.modelDescription.modelName + fmu.isZeroState = (length(fmu.modelDescription.stateValueReferences) == 0) if isa(type, fmi2Type) fmu.type = type @@ -455,6 +456,7 @@ function fmi2Instantiate!(fmu::FMU2; component.visible = visible component.jacobianUpdate! = fmi2SampleJacobian! + # setting a jacobian update function dependent on DIrectionalDerivatives-Functionality is present in the FMU updFct = nothing if fmi2ProvidesDirectionalDerivative(fmu) updFct = (jac, ∂f_refs, ∂x_refs) -> fmi2GetJacobian!(jac.mtx, component, ∂f_refs, ∂x_refs) @@ -466,6 +468,8 @@ function fmi2Instantiate!(fmu::FMU2; component.B = FMICore.FMUJacobian{fmi2Real, fmi2ValueReference}(fmu.modelDescription.derivativeValueReferences, fmu.modelDescription.inputValueReferences, updFct) component.C = FMICore.FMUJacobian{fmi2Real, fmi2ValueReference}(fmu.modelDescription.outputValueReferences, fmu.modelDescription.stateValueReferences, updFct) component.D = FMICore.FMUJacobian{fmi2Real, fmi2ValueReference}(fmu.modelDescription.outputValueReferences, fmu.modelDescription.inputValueReferences, updFct) + component.E = FMICore.FMUJacobian{fmi2Real, fmi2ValueReference}(fmu.modelDescription.derivativeValueReferences, isnothing(fmu.optim_p_refs) ? Array{fmi2ValueReference,1}() : fmu.optim_p_refs, updFct) + component.F = FMICore.FMUJacobian{fmi2Real, fmi2ValueReference}(fmu.modelDescription.outputValueReferences, isnothing(fmu.optim_p_refs) ? Array{fmi2ValueReference,1}() : fmu.optim_p_refs, updFct) # register component for current thread fmu.threadComponents[Threads.threadid()] = component @@ -1022,7 +1026,8 @@ function fmi2Set(comp::FMU2Component, vrs::fmi2ValueReferenceFormat, srcArray::A @assert isa(srcArray[i], String) "fmi2Set(...): Unknown data type for value reference `$(vr)` at index $(i), should be `String`, is `$(typeof(srcArray[i]))`." retcodes[i] = fmi2SetString(comp, vr, srcArray[i]) elseif mv.Enumeration != nothing - @warn "fmi2Set(...): Currently not implemented for fmi2Enum." + @assert isa(srcArray[i], Union{Real, Integer}) "fmi2Set(...): Unknown data type for value reference `$(vr)` at index $(i), should be `Enumeration` (`Integer`), is `$(typeof(srcArray[i]))`." + retcodes[i] = fmi2SetInteger(comp, vr, Integer(srcArray[i])) else @assert false "fmi2Set(...): Unknown data type for value reference `$(vr)` at index $(i), is `$(mv.datatype.datatype)`." end diff --git a/src/FMI2/prep.jl b/src/FMI2/prep.jl index 6d0049e..2839b55 100644 --- a/src/FMI2/prep.jl +++ b/src/FMI2/prep.jl @@ -113,7 +113,7 @@ function prepareSolveFMU(fmu::FMU2, #retcode = fmi2SetContinuousStates(c, x0) #@assert retcode == fmi2StatusOK "fmi2Simulate(...): Setting initial state failed with return code $(retcode)." retcodes = fmi2Set(c, fmu.modelDescription.stateValueReferences, x0; filter=setBeforeInitialization) - @assert all(retcodes .== fmi2StatusOK) "fmi2Simulate(...): Setting initial inputs failed with return code $(retcode)." + @assert all(retcodes .== fmi2StatusOK) "fmi2Simulate(...): Setting initial states failed with return code $(retcode)." end # enter (hard) @@ -125,13 +125,13 @@ function prepareSolveFMU(fmu::FMU2, # parameters if parameters !== nothing retcodes = fmi2Set(c, collect(keys(parameters)), collect(values(parameters)); filter=setInInitialization) - @assert all(retcodes .== fmi2StatusOK) "fmi2Simulate(...): Setting initial parameters failed with return code $(retcode)." + @assert all(retcodes .== fmi2StatusOK) "fmi2Simulate(...): Setting initial parameters failed with return code $(retcodes)." end # inputs if inputs !== nothing retcodes = fmi2Set(c, collect(keys(inputs)), collect(values(inputs)); filter=setInInitialization) - @assert all(retcodes .== fmi2StatusOK) "fmi2Simulate(...): Setting initial inputs failed with return code $(retcode)." + @assert all(retcodes .== fmi2StatusOK) "fmi2Simulate(...): Setting initial inputs failed with return code $(retcodes)." end # start state @@ -139,7 +139,7 @@ function prepareSolveFMU(fmu::FMU2, #retcode = fmi2SetContinuousStates(c, x0) #@assert retcode == fmi2StatusOK "fmi2Simulate(...): Setting initial state failed with return code $(retcode)." retcodes = fmi2Set(c, fmu.modelDescription.stateValueReferences, x0; filter=setInInitialization) - @assert all(retcodes .== fmi2StatusOK) "fmi2Simulate(...): Setting initial inputs failed with return code $(retcode)." + @assert all(retcodes .== fmi2StatusOK) "fmi2Simulate(...): Setting initial inputs failed with return code $(retcodes)." end # exit setup (hard) @@ -153,7 +153,7 @@ function prepareSolveFMU(fmu::FMU2, # ME specific if type == fmi2TypeModelExchange - if x0 == nothing + if x0 == nothing && !c.fmu.isZeroState x0 = fmi2GetContinuousStates(c) end diff --git a/src/FMI2/sens.jl b/src/FMI2/sens.jl index 1ba4474..39d0b73 100644 --- a/src/FMI2/sens.jl +++ b/src/FMI2/sens.jl @@ -62,19 +62,23 @@ Not all options are available for any FMU type, e.g. setting state is not suppor - `x`: An array containing the states to be set. Not supported by CS-FMUs. - `u`: An array containing the inputs to be set. - `u_refs`: An array of value references to indicate which system inputs want to be set. +- `p`: An array of FMU parameters to be set. +- `p_refs`: An array of parameter references to indicate which system parameter sensitivities need to be determined. - `t`: A scalar value holding the system time to be set. # Returns (as Tuple) - `y::Union{AbstractVector{<:Real}, Nothing}`: The system output `y` (if requested, otherwise `nothing`). - `dx::Union{AbstractVector{<:Real}, Nothing}`: The system state-derivaitve (if ME-FMU, otherwise `nothing`). """ -function (fmu::FMU2)(;dx::Union{AbstractVector{<:Real}, Nothing}=nothing, - y::Union{AbstractVector{<:Real}, Nothing}=nothing, - y_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}=nothing, - x::Union{AbstractVector{<:Real}, Nothing}=nothing, - u::Union{AbstractVector{<:Real}, Nothing}=nothing, - u_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}=nothing, - t::Union{Real, Nothing}=nothing) +function (fmu::FMU2)(;dx::AbstractVector{<:Real}=Vector{fmi2Real}(), + y::AbstractVector{<:Real}=Vector{fmi2Real}(), + y_refs::AbstractVector{<:fmi2ValueReference}=Vector{fmi2ValueReference}(), + x::AbstractVector{<:Real}=Vector{fmi2Real}(), + u::AbstractVector{<:Real}=Vector{fmi2Real}(), + u_refs::AbstractVector{<:fmi2ValueReference}=Vector{fmi2ValueReference}(), + p::AbstractVector{<:Real}=fmu.optim_p, + p_refs::AbstractVector{<:fmi2ValueReference}=fmu.optim_p_refs, + t::Real=-1.0) c = nothing @@ -87,7 +91,11 @@ function (fmu::FMU2)(;dx::Union{AbstractVector{<:Real}, Nothing}=nothing, fmi2ExitInitializationMode(c) end - c(;dx=dx, y=y, y_refs=y_refs, x=x, u=u, u_refs=u_refs, t=t) + if t == -1.0 + t = c.next_t + end + + c(;dx=dx, y=y, y_refs=y_refs, x=x, u=u, u_refs=u_refs, p=p, p_refs=p_refs, t=t) end """ @@ -110,33 +118,38 @@ Not all options are available for any FMU type, e.g. setting state is not suppor - `x`: An array containing the states to be set. Not supported by CS-FMUs. - `u`: An array containing the inputs to be set. - `u_refs`: An array of value references to indicate which system inputs want to be set. +- `p`: An array of FMU parameters to be set. +- `p_refs`: An array of parameter references to indicate which system parameter sensitivities need to be determined. - `t`: A scalar value holding the system time to be set. # Returns (as Tuple) - `y::Union{AbstractVector{<:Real}, Nothing}`: The system output `y` (if requested, otherwise `nothing`). - `dx::Union{AbstractVector{<:Real}, Nothing}`: The system state-derivaitve (if ME-FMU, otherwise `nothing`). """ -function (c::FMU2Component)(;dx::Union{AbstractVector{<:Real}, Nothing}=nothing, - y::Union{AbstractVector{<:Real}, Nothing}=nothing, - y_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}=nothing, - x::Union{AbstractVector{<:Real}, Nothing}=nothing, - u::Union{AbstractVector{<:Real}, Nothing}=nothing, - u_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}=nothing, - t::Union{Real, Nothing}=nothing) - - if y_refs != nothing && length(y_refs) > 0 - if y == nothing +function (c::FMU2Component)(;dx::AbstractVector{<:Real}=Vector{fmi2Real}(), + y::AbstractVector{<:Real}=Vector{fmi2Real}(), + y_refs::AbstractVector{<:fmi2ValueReference}=Vector{fmi2ValueReference}(), + x::AbstractVector{<:Real}=Vector{fmi2Real}(), + u::AbstractVector{<:Real}=Vector{fmi2Real}(), + u_refs::AbstractVector{<:fmi2ValueReference}=Vector{fmi2ValueReference}(), + p::AbstractVector{<:Real}=c.fmu.optim_p, + p_refs::AbstractVector{<:fmi2ValueReference}=c.fmu.optim_p_refs, + t::Real=c.next_t) + + if length(y_refs) > 0 + if length(y) <= 0 y = zeros(fmi2Real, length(y_refs)) end end - @assert y == nothing || (length(y) == length(y_refs)) "Length of `y` must match length of `y_refs`." - @assert u == nothing || (length(u) == length(u_refs)) "Length of `u` must match length of `u_refs`." + @assert (length(y) == length(y_refs)) "Length of `y` must match length of `y_refs`." + @assert (length(u) == length(u_refs)) "Length of `u` must match length of `u_refs`." + @assert (length(p) == length(p_refs)) "Length of `p` must match length of `p_refs`." if fmi2IsModelExchange(c.fmu) if c.type == fmi2TypeModelExchange::fmi2Type - if dx == nothing + if length(dx) <= 0 dx = zeros(fmi2Real, fmi2GetNumberOfStates(c.fmu.modelDescription)) end end @@ -144,85 +157,86 @@ function (c::FMU2Component)(;dx::Union{AbstractVector{<:Real}, Nothing}=nothing, if fmi2IsCoSimulation(c.fmu) if c.type == fmi2TypeCoSimulation::fmi2Type - @assert dx == nothing "Keyword `dx != nothing` is invalid for CS-FMUs. Setting a state-derivative is not possible in CS." - @assert x == nothing "Keyword `x != nothing` is invalid for CS-FMUs. Setting a state is not possible in CS." - @assert t == nothing "Keyword `t != nothing` is invalid for CS-FMUs. Setting explicit time is not possible in CS." + @assert length(dx) <= 0 "Keyword `dx != nothing` is invalid for CS-FMUs. Setting a state-derivative is not possible in CS." + @assert length(x) <= 0 "Keyword `x != nothing` is invalid for CS-FMUs. Setting a state is not possible in CS." + @assert t < 0.0 "Keyword `t != nothing` is invalid for CS-FMUs. Setting explicit time is not possible in CS." end end - # ToDo: This is necessary, because NonconvexUtils/ForwardDiff can't handle arguments with type `Nothing`. - if t == nothing - t = -1.0 - end - - # ToDo: This is necessary, because NonconvexUtils/ForwardDiff can't handle arguments with type `Ptr{Nothing}`. + # ToDo: This is necessary, because ForwardDiffChainRules.jl can't handle arguments with type `Ptr{Nothing}`. cRef = nothing ignore_derivatives() do cRef = pointer_from_objref(c) cRef = UInt64(cRef) end - # ToDo: This is necessary, because NonconvexUtils/ForwardDiff can't handle arguments with type `Nothing`. - if u == nothing || length(u) <= 0 - return eval!(cRef, dx, y, y_refs, x, t) - - elseif x == nothing || length(x) <= 0 - return eval!(cRef, dx, y, y_refs, u, u_refs, t) - - end - - return eval!(cRef, dx, y, y_refs, x, u, u_refs, t) + return eval!(cRef, dx, y, y_refs, x, u, u_refs, p, p_refs, t) end -function _eval!(cRef::UInt64, - dx::Union{AbstractVector{<:Real}, Nothing}, - y::Union{AbstractVector{<:Real}, Nothing}, - y_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, - x::Union{AbstractVector{<:Real}, Nothing}, - u::Union{AbstractVector{<:Real}, Nothing}, - u_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, - t::Union{Real, Nothing}) +function eval!(cRef::UInt64, + dx::AbstractVector{<:Real}, + y::AbstractVector{<:Real}, + y_refs::AbstractVector{<:fmi2ValueReference}, + x::AbstractVector{<:Real}, + u::AbstractVector{<:Real}, + u_refs::AbstractVector{<:fmi2ValueReference}, + p::AbstractVector{<:Real}, + p_refs::AbstractVector{<:fmi2ValueReference}, + t::Real) c = unsafe_pointer_to_objref(Ptr{Nothing}(cRef)) - @assert x == nothing || !isdual(x) "_eval!(...): Wrong dispatched: `x` is ForwardDiff.Dual, please open an issue with MWE." - @assert u == nothing || !isdual(u) "_eval!(...): Wrong dispatched: `u` is ForwardDiff.Dual, please open an issue with MWE." - @assert t == nothing || !isdual(t) "_eval!(...): Wrong dispatched: `t` is ForwardDiff.Dual, please open an issue with MWE." + @assert (!isdual(x) && !istracked(x)) "eval!(...): Wrong dispatched: `x` is ForwardDiff.Dual/ReverseDiff.TrackedReal, please open an issue with MWE." + @assert (!isdual(u) && !istracked(u)) "eval!(...): Wrong dispatched: `u` is ForwardDiff.Dual/ReverseDiff.TrackedReal, please open an issue with MWE." + @assert (!isdual(t) && !istracked(t)) "eval!(...): Wrong dispatched: `t` is ForwardDiff.Dual/ReverseDiff.TrackedReal, please open an issue with MWE." + @assert (!isdual(p) && !istracked(p)) "eval!(...): Wrong dispatched: `p` is ForwardDiff.Dual/ReverseDiff.TrackedReal, please open an issue with MWE." x = unsense(x) t = unsense(t) u = unsense(u) + # p = unsense(p) # no need to unsense `p` because it is not beeing used further # set state - if x != nothing + if length(x) > 0 && !c.fmu.isZeroState fmi2SetContinuousStates(c, x) end # set time - if t != nothing && t >= 0.0 + if t >= 0.0 fmi2SetTime(c, t) end # set input - if u != nothing + if length(u) > 0 fmi2SetReal(c, u_refs, u) end # get derivative - if dx != nothing + if length(dx) > 0 if isdual(dx) - #@info "dx is dual!" - dx_tmp = collect(ForwardDiff.value(e) for e in dx) - fmi2GetDerivatives!(c, dx_tmp) + + dx_tmp = nothing + + if c.fmu.isZeroState + dx_tmp = [1.0] + else + dx_tmp = collect(ForwardDiff.value(e) for e in dx) + fmi2GetDerivatives!(c, dx_tmp) + end + T, V, N = fd_eltypes(dx) dx[:] = collect(ForwardDiff.Dual{T, V, N}(dx_tmp[i], ForwardDiff.partials(dx[i]) ) for i in 1:length(dx)) else - fmi2GetDerivatives!(c, dx) + if c.fmu.isZeroState + dx[:] = [1.0] + else + fmi2GetDerivatives!(c, dx) + end end end # get output - if y != nothing + if length(y) > 0 if isdual(y) #@info "y is dual!" y_tmp = collect(ForwardDiff.value(e) for e in y) @@ -237,14 +251,6 @@ function _eval!(cRef::UInt64, end end - if isnothing(y) - y = Array{Float64, 1}() - end - - if isnothing(dx) - dx = Array{Float64, 1}() - end - if c.fmu.executionConfig.concat_y_dx return vcat(y, dx) # [y..., dx...] else @@ -252,7 +258,8 @@ function _eval!(cRef::UInt64, end end -function _frule(Δtuple, +function ChainRulesCore.frule(Δtuple, + ::typeof(eval!), cRef, dx, y, @@ -260,10 +267,12 @@ function _frule(Δtuple, x, u, u_refs, + p, + p_refs, t) Δtuple = undual(Δtuple) - Δself, ΔcRef, Δdx, Δy, Δy_refs, Δx, Δu, Δu_refs, Δt = Δtuple + Δself, ΔcRef, Δdx, Δy, Δy_refs, Δx, Δu, Δu_refs, Δp, Δp_refs, Δt = Δtuple ### ToDo: Somehow, ForwardDiff enters with all types beeing Float64, this needs to be corrected. @@ -275,30 +284,32 @@ function _frule(Δtuple, t = undual(t) u = undual(u) x = undual(x) + + p = undual(p) y_refs = undual(y_refs) - if y_refs != nothing - y_refs = convert(Array{UInt32,1}, y_refs) - end - + y_refs = convert(Array{UInt32,1}, y_refs) + u_refs = undual(u_refs) - if u_refs != nothing - u_refs = convert(Array{UInt32,1}, u_refs) - end - + u_refs = convert(Array{UInt32,1}, u_refs) + + p_refs = undual(p_refs) + p_refs = convert(Array{UInt32,1}, p_refs) + ### c = unsafe_pointer_to_objref(Ptr{Nothing}(cRef)) - outputs = (y != nothing && length(y_refs) > 0) - inputs = (u != nothing && length(u) > 0) - derivatives = (dx != nothing && length(dx) > 0) - states = (x != nothing && length(x) > 0) - times = (t != nothing && t >= 0.0) + outputs = (length(y_refs) > 0) + inputs = (length(u_refs) > 0) + derivatives = (length(dx) > 0) + states = (length(x) > 0) + times = (t >= 0.0) + parameters = (length(p_refs) > 0) - Ω = _eval!(cRef, dx, y, y_refs, x, u, u_refs, t) + Ω = eval!(cRef, dx, y, y_refs, x, u, u_refs, p, p_refs, t) - # time, states and inputs where already set in `_eval!`, no need to repeat it here + # time, states and inputs where already set in `eval!`, no need to repeat it here ∂y = ZeroTangent() ∂dx = ZeroTangent() @@ -309,15 +320,17 @@ function _frule(Δtuple, Δx = convert(Vector{fmi2Real}, Δx) end - if derivatives && states - ∂dx += fmi2JVP!(c, :A, c.fmu.modelDescription.derivativeValueReferences, c.fmu.modelDescription.stateValueReferences, Δx) - c.solution.evals_∂ẋ_∂x += 1 - #@info "$(Δx)" - end + if states + if derivatives + ∂dx += fmi2JVP!(c, :A, c.fmu.modelDescription.derivativeValueReferences, c.fmu.modelDescription.stateValueReferences, Δx) + c.solution.evals_∂ẋ_∂x += 1 + #@info "$(Δx)" + end - if outputs && states - ∂y += fmi2JVP!(c, :C, y_refs, c.fmu.modelDescription.stateValueReferences, Δx) - c.solution.evals_∂y_∂x += 1 + if outputs + ∂y += fmi2JVP!(c, :C, y_refs, c.fmu.modelDescription.stateValueReferences, Δx) + c.solution.evals_∂y_∂x += 1 + end end end @@ -328,20 +341,41 @@ function _frule(Δtuple, Δu = convert(Vector{fmi2Real}, Δu) end - if derivatives && inputs - ∂dx += fmi2JVP!(c, :B, c.fmu.modelDescription.derivativeValueReferences, u_refs, Δu) - c.solution.evals_∂ẋ_∂u += 1 + if inputs + if derivatives + ∂dx += fmi2JVP!(c, :B, c.fmu.modelDescription.derivativeValueReferences, u_refs, Δu) + c.solution.evals_∂ẋ_∂u += 1 + end + + if outputs + ∂y += fmi2JVP!(c, :D, y_refs, u_refs, Δu) + c.solution.evals_∂y_∂u += 1 + end end + end + + if Δp != NoTangent() && length(Δp) > 0 - if outputs && inputs - ∂y += fmi2JVP!(c, :D, y_refs, u_refs, Δu) - c.solution.evals_∂y_∂u += 1 + if !isa(Δp, AbstractVector{fmi2Real}) + Δp = convert(Vector{fmi2Real}, Δp) + end + + if parameters + if derivatives + ∂dx += fmi2JVP!(c, :E, c.fmu.modelDescription.derivativeValueReferences, p_refs, Δp) + c.solution.evals_∂ẋ_∂p += 1 + end + + if outputs + ∂y += fmi2JVP!(c, :F, y_refs, p_refs, Δp) + c.solution.evals_∂y_∂p += 1 + end end end if c.fmu.executionConfig.eval_t_gradients # partial time derivatives are not part of the FMI standard, so must be sampled in any case - if Δt != NoTangent() && t != nothing && times && (derivatives || outputs) + if Δt != NoTangent() && times && (derivatives || outputs) dt = 1e-6 # ToDo: Find a better value, e.g. based on the current solver step size @@ -386,7 +420,7 @@ function _frule(Δtuple, end end - #@info "frule: ∂y=$(∂y) ∂dx=$(∂dx)" + @debug "frule: ∂y=$(∂y) ∂dx=$(∂dx)" ∂Ω = nothing if c.fmu.executionConfig.concat_y_dx @@ -410,26 +444,30 @@ function isZeroTangent(d::AbstractArray{<:ZeroTangent}) return true end -function _rrule(cRef, +function ChainRulesCore.rrule(::typeof(eval!), + cRef, dx, y, y_refs, x, u, u_refs, + p, + p_refs, t) @assert !isa(cRef, FMU2Component) "Wrong dispatched!" c = unsafe_pointer_to_objref(Ptr{Nothing}(cRef)) - outputs = (!isnothing(y) && length(y_refs) > 0) - inputs = (!isnothing(u) && length(u) > 0) - derivatives = (!isnothing(dx) && length(dx) > 0) - states = (!isnothing(x) && length(x) > 0) - times = (!isnothing(t) && t >= 0.0) + outputs = (length(y_refs) > 0) + inputs = (length(u_refs) > 0) + derivatives = (length(dx) > 0) + states = (length(x) > 0) + times = (t >= 0.0) + parameters = (length(p_refs) > 0) - Ω = _eval!(cRef, dx, y, y_refs, x, u, u_refs, t) + Ω = eval!(cRef, dx, y, y_refs, x, u, u_refs, p, p_refs, t) # if !inputs # Ω = eval!(cRef, dx, y, y_refs, x, t) @@ -481,34 +519,57 @@ function _rrule(cRef, n_dx_x = ZeroTangent() n_dx_u = ZeroTangent() + n_dx_p = ZeroTangent() n_dx_t = ZeroTangent() + n_y_x = ZeroTangent() n_y_u = ZeroTangent() + n_y_p = ZeroTangent() n_y_t = ZeroTangent() - #@info "rrule pullback ȳ, d̄x = $(ȳ), $(d̄x)" + @debug "rrule pullback ȳ, d̄x = $(ȳ), $(d̄x)" dx_refs = c.fmu.modelDescription.derivativeValueReferences x_refs = c.fmu.modelDescription.stateValueReferences - if derivatives && states - n_dx_x = fmi2VJP!(c, :A, dx_refs, x_refs, d̄x) - c.solution.evals_∂ẋ_∂x += 1 - end + if derivatives + if states + n_dx_x = fmi2VJP!(c, :A, dx_refs, x_refs, d̄x) + c.solution.evals_∂ẋ_∂x += 1 + end - if derivatives && inputs - n_dx_u = fmi2VJP!(c, :B, dx_refs, u_refs, d̄x) - c.solution.evals_∂ẋ_∂u += 1 - end + if inputs + n_dx_u = fmi2VJP!(c, :B, dx_refs, u_refs, d̄x) + c.solution.evals_∂ẋ_∂u += 1 + end - if outputs && states - n_y_x = fmi2VJP!(c, :C, y_refs, x_refs, ȳ) - c.solution.evals_∂y_∂x += 1 + if parameters + n_dx_p = fmi2VJP!(c, :E, dx_refs, p_refs, d̄x) + c.solution.evals_∂ẋ_∂p += 1 + + # if rand(1:100) == 1 + # #@info "$(c.E.mtx)" + # #@info "$(p_refs)" + # @info "$(fmi2GetJacobian(c, dx_refs, p_refs))" + # end + end end - if outputs && inputs - n_y_u = fmi2VJP!(c, :D, y_refs, u_refs, ȳ) - c.solution.evals_∂y_∂u += 1 + if outputs + if states + n_y_x = fmi2VJP!(c, :C, y_refs, x_refs, ȳ) + c.solution.evals_∂y_∂x += 1 + end + + if inputs + n_y_u = fmi2VJP!(c, :D, y_refs, u_refs, ȳ) + c.solution.evals_∂y_∂u += 1 + end + + if parameters + n_y_p = fmi2VJP!(c, :F, y_refs, p_refs, ȳ) + c.solution.evals_∂y_∂p += 1 + end end if c.fmu.executionConfig.eval_t_gradients @@ -572,305 +633,156 @@ function _rrule(cRef, ȳ_refs = ZeroTangent() t̄ = n_y_t + n_dx_t - opt = Array{Any, 1}() - if !isnothing(x) - x̄ = n_y_x + n_dx_x - push!(opt, x̄) - end - if !isnothing(u) - ū = n_y_u + n_dx_u - ū_refs = ZeroTangent() - push!(opt, ū) - push!(opt, ū_refs) - end + x̄ = n_y_x + n_dx_x - #@info "rrule: $((ret...,))" + ū = n_y_u + n_dx_u + ū_refs = ZeroTangent() - return (f̄, c̄Ref, d̄x, ȳ, ȳ_refs, opt..., t̄) + p̄ = n_y_p + n_dx_p + p̄_refs = ZeroTangent() + + @debug "rrule: $((f̄, c̄Ref, d̄x, ȳ, ȳ_refs, x̄, ū, ū_refs, p̄, p̄_refs, t̄))" + + return (f̄, c̄Ref, d̄x, ȳ, ȳ_refs, x̄, ū, ū_refs, p̄, p̄_refs, t̄) end return (Ω, eval_pullback) end -# EVAL! WITH `x` and `u` - -function eval!(cRef::UInt64, - dx::Union{AbstractVector{<:Real}, Nothing}, - y::Union{AbstractVector{<:Real}, Nothing}, - y_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, - x::Union{AbstractVector{<:Real}, Nothing}, - u::Union{AbstractVector{<:Real}, Nothing}, - u_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, - t::Union{Real, Nothing}) - - return _eval!(cRef, dx, y, y_refs, x, u, u_refs, t) -end - -function ChainRulesCore.frule((Δself, ΔcRef, Δdx, Δy, Δy_refs, Δx, Δu, Δu_refs, Δt), - ::typeof(eval!), - cRef, - dx, - y, - y_refs, - x, - u, - u_refs, - t) - - return _frule((Δself, ΔcRef, Δdx, Δy, Δy_refs, Δx, Δu, Δu_refs, Δt), - cRef, dx, y, y_refs, x, u, u_refs, t) -end - -function ChainRulesCore.rrule(::typeof(eval!), - cRef, - dx, - y, - y_refs, - x, - u, - u_refs, - t) - - return _rrule(cRef, dx, y, y_refs, x, u, u_refs, t) -end - +# dx, y, x, u, t @ForwardDiff_frule eval!(cRef::UInt64, - dx::Union{AbstractVector{<:ForwardDiff.Dual}, Nothing}, - y::Union{AbstractVector{<:ForwardDiff.Dual}, Nothing}, - y_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, - x::Union{AbstractVector{<:ForwardDiff.Dual}, Nothing}, - u::Union{AbstractVector{<:ForwardDiff.Dual}, Nothing}, - u_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, - t::Union{ForwardDiff.Dual, Nothing}) - -@ForwardDiff_frule eval!(cRef::UInt64, - dx::Union{AbstractVector{<:Real}, Nothing}, - y::Union{AbstractVector{<:Real}, Nothing}, - y_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, - x::Union{AbstractVector{<:Real}, Nothing}, - u::Union{AbstractVector{<:Real}, Nothing}, - u_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, - t::Union{ForwardDiff.Dual, Nothing}) - -@ForwardDiff_frule eval!(cRef::UInt64, - dx::Union{AbstractVector{<:Real}, Nothing}, - y::Union{AbstractVector{<:Real}, Nothing}, - y_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, - x::Union{AbstractVector{<:ForwardDiff.Dual}, Nothing}, - u::Union{AbstractVector{<:Real}, Nothing}, - u_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, - t::Union{Real, Nothing}) - -@ForwardDiff_frule eval!(cRef::UInt64, - dx::Union{AbstractVector{<:Real}, Nothing}, - y::Union{AbstractVector{<:Real}, Nothing}, - y_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, - x::Union{AbstractVector{<:Real}, Nothing}, - u::Union{AbstractVector{<:ForwardDiff.Dual}, Nothing}, - u_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, - t::Union{Real, Nothing}) + dx ::AbstractVector{<:ForwardDiff.Dual}, + y ::AbstractVector{<:ForwardDiff.Dual}, + y_refs::AbstractVector{<:fmi2ValueReference}, + x ::AbstractVector{<:ForwardDiff.Dual}, + u ::AbstractVector{<:ForwardDiff.Dual}, + u_refs::AbstractVector{<:fmi2ValueReference}, + p ::AbstractVector{<:Real}, + p_refs::AbstractVector{<:fmi2ValueReference}, + t ::ForwardDiff.Dual) @grad_from_chainrules eval!(cRef::UInt64, - dx::Union{AbstractVector{<:ReverseDiff.TrackedReal}, Nothing}, - y::Union{AbstractVector{<:ReverseDiff.TrackedReal}, Nothing}, - y_refs::Union{AbstractVector{<:UInt32}, Nothing}, - x::Union{AbstractVector{<:ReverseDiff.TrackedReal}, Nothing}, - u::Union{AbstractVector{<:ReverseDiff.TrackedReal}, Nothing}, - u_refs::Union{AbstractVector{<:UInt32}, Nothing}, - t::Union{ReverseDiff.TrackedReal, Nothing}) - -@grad_from_chainrules eval!(cRef::UInt64, - dx::Union{AbstractVector{<:Real}, Nothing}, - y::Union{AbstractVector{<:Real}, Nothing}, - y_refs::Union{AbstractVector{<:UInt32}, Nothing}, - x::Union{AbstractVector{<:Real}, Nothing}, - u::Union{AbstractVector{<:Real}, Nothing}, - u_refs::Union{AbstractVector{<:UInt32}, Nothing}, - t::Union{ReverseDiff.TrackedReal, Nothing}) - -@grad_from_chainrules eval!(cRef::UInt64, - dx::Union{AbstractVector{<:Real}, Nothing}, - y::Union{AbstractVector{<:Real}, Nothing}, - y_refs::Union{AbstractVector{<:UInt32}, Nothing}, - x::Union{AbstractVector{<:ReverseDiff.TrackedReal}, Nothing}, - u::Union{AbstractVector{<:Real}, Nothing}, - u_refs::Union{AbstractVector{<:UInt32}, Nothing}, - t::Union{Real, Nothing}) - -@grad_from_chainrules eval!(cRef::UInt64, - dx::Union{AbstractVector{<:Real}, Nothing}, - y::Union{AbstractVector{<:Real}, Nothing}, - y_refs::Union{AbstractVector{<:UInt32}, Nothing}, - x::Union{AbstractVector{<:Real}, Nothing}, - u::Union{AbstractVector{<:ReverseDiff.TrackedReal}, Nothing}, - u_refs::Union{AbstractVector{<:UInt32}, Nothing}, - t::Union{Real, Nothing}) - -# EVAL! WITH `x`, WITHOUT `u` - -function eval!(cRef::UInt64, - dx::Union{AbstractVector{<:Real}, Nothing}, - y::Union{AbstractVector{<:Real}, Nothing}, - y_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, - x::Union{AbstractVector{<:Real}, Nothing}, - t::Union{Real, Nothing}) - - t = unsense(t) - x = unsense(x) - - return _eval!(cRef, dx, y, y_refs, x, nothing, nothing, t) -end - -function ChainRulesCore.frule((Δself, ΔcRef, Δdx, Δy, Δy_refs, Δx, Δt), - ::typeof(eval!), - cRef, - dx, - y, - y_refs, - x, - t) - - return _frule((Δself, ΔcRef, Δdx, Δy, Δy_refs, Δx, NoTangent(), NoTangent(), Δt), - cRef, dx, y, y_refs, x, nothing, nothing, t) -end - -function ChainRulesCore.rrule(::typeof(eval!), - cRef, - dx, - y, - y_refs, - x, - t) - - return _rrule(cRef, dx, y, y_refs, x, nothing, nothing, t) -end - -@ForwardDiff_frule eval!(cRef::UInt64, -dx::Union{AbstractVector{<:ForwardDiff.Dual}, Nothing}, -y::Union{AbstractVector{<:ForwardDiff.Dual}, Nothing}, -y_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, -x::Union{AbstractVector{<:ForwardDiff.Dual}, Nothing}, -t::Union{ForwardDiff.Dual, Nothing}) - -@ForwardDiff_frule eval!(cRef::UInt64, -dx::Union{AbstractVector{<:Real}, Nothing}, -y::Union{AbstractVector{<:Real}, Nothing}, -y_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, -x::Union{AbstractVector{<:Real}, Nothing}, -t::Union{ForwardDiff.Dual, Nothing}) - + dx ::AbstractVector{<:ReverseDiff.TrackedReal}, + y ::AbstractVector{<:ReverseDiff.TrackedReal}, + y_refs::AbstractVector{<:UInt32}, + x ::AbstractVector{<:ReverseDiff.TrackedReal}, + u ::AbstractVector{<:ReverseDiff.TrackedReal}, + u_refs::AbstractVector{<:UInt32}, + p ::AbstractVector{<:Real}, + p_refs::AbstractVector{<:UInt32}, + t ::ReverseDiff.TrackedReal) + +# x, p @ForwardDiff_frule eval!(cRef::UInt64, -dx::Union{AbstractVector{<:Real}, Nothing}, -y::Union{AbstractVector{<:Real}, Nothing}, -y_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, -x::Union{AbstractVector{<:ForwardDiff.Dual}, Nothing}, -t::Union{Real, Nothing}) - -@grad_from_chainrules eval!(cRef::UInt64, -dx::Union{AbstractVector{<:ReverseDiff.TrackedReal}, Nothing}, -y::Union{AbstractVector{<:ReverseDiff.TrackedReal}, Nothing}, -y_refs::Union{AbstractVector{<:UInt32}, Nothing}, -x::Union{AbstractVector{<:ReverseDiff.TrackedReal}, Nothing}, -t::Union{ReverseDiff.TrackedReal, Nothing}) + dx ::AbstractVector{<:Real}, + y ::AbstractVector{<:Real}, + y_refs::AbstractVector{<:fmi2ValueReference}, + x ::AbstractVector{<:ForwardDiff.Dual}, + u ::AbstractVector{<:Real}, + u_refs::AbstractVector{<:fmi2ValueReference}, + p ::AbstractVector{<:ForwardDiff.Dual}, + p_refs::AbstractVector{<:fmi2ValueReference}, + t ::Real) @grad_from_chainrules eval!(cRef::UInt64, -dx::Union{AbstractVector{<:Real}, Nothing}, -y::Union{AbstractVector{<:Real}, Nothing}, -y_refs::Union{AbstractVector{<:UInt32}, Nothing}, -x::Union{AbstractVector{<:Real}, Nothing}, -t::Union{ReverseDiff.TrackedReal, Nothing}) + dx ::AbstractVector{<:Real}, + y ::AbstractVector{<:Real}, + y_refs::AbstractVector{<:UInt32}, + x ::AbstractVector{<:ReverseDiff.TrackedReal}, + u ::AbstractVector{<:Real}, + u_refs::AbstractVector{<:UInt32}, + p ::AbstractVector{<:ReverseDiff.TrackedReal}, + p_refs::AbstractVector{<:UInt32}, + t ::Real) + +# t +@ForwardDiff_frule eval!(cRef::UInt64, + dx ::AbstractVector{<:Real}, + y ::AbstractVector{<:Real}, + y_refs::AbstractVector{<:fmi2ValueReference}, + x ::AbstractVector{<:Real}, + u ::AbstractVector{<:Real}, + u_refs::AbstractVector{<:fmi2ValueReference}, + p ::AbstractVector{<:Real}, + p_refs::AbstractVector{<:fmi2ValueReference}, + t ::ForwardDiff.Dual) @grad_from_chainrules eval!(cRef::UInt64, -dx::Union{AbstractVector{<:Real}, Nothing}, -y::Union{AbstractVector{<:Real}, Nothing}, -y_refs::Union{AbstractVector{<:UInt32}, Nothing}, -x::Union{AbstractVector{<:ReverseDiff.TrackedReal}, Nothing}, -t::Union{Real, Nothing}) - -# EVAL! WITH `u` WITHOUT `x` - -function eval!(cRef::UInt64, - dx::Union{AbstractVector{<:Real}, Nothing}, - y::Union{AbstractVector{<:Real}, Nothing}, - y_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, - u::Union{AbstractVector{<:Real}, Nothing}, - u_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, - t::Union{Real, Nothing}) - - return _eval!(cRef, dx, y, y_refs, nothing, u, u_refs, t) -end - -function ChainRulesCore.frule((Δself, ΔcRef, Δdx, Δy, Δy_refs, Δu, Δu_refs, Δt), - ::typeof(eval!), - cRef, - dx, - y, - y_refs, - u, - u_refs, - t) - - return _frule((Δself, ΔcRef, Δdx, Δy, Δy_refs, NoTangent(), Δu, Δu_refs, Δt), - cRef, dx, y, y_refs, nothing, u, u_refs, t) -end - -function ChainRulesCore.rrule(::typeof(eval!), - cRef, - dx, - y, - y_refs, - u, - u_refs, - t) - - return _rrule(cRef, dx, y, y_refs, nothing, u, u_refs, t) -end - -@ForwardDiff_frule eval!(cRef::UInt64, -dx::Union{AbstractVector{<:ForwardDiff.Dual}, Nothing}, -y::Union{AbstractVector{<:ForwardDiff.Dual}, Nothing}, -y_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, -u::Union{AbstractVector{<:ForwardDiff.Dual}, Nothing}, -u_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, -t::Union{ForwardDiff.Dual, Nothing}) - + dx ::AbstractVector{<:Real}, + y ::AbstractVector{<:Real}, + y_refs::AbstractVector{<:UInt32}, + x ::AbstractVector{<:Real}, + u ::AbstractVector{<:Real}, + u_refs::AbstractVector{<:UInt32}, + p ::AbstractVector{<:Real}, + p_refs::AbstractVector{<:UInt32}, + t ::ReverseDiff.TrackedReal) + +# x @ForwardDiff_frule eval!(cRef::UInt64, -dx::Union{AbstractVector{<:Real}, Nothing}, -y::Union{AbstractVector{<:Real}, Nothing}, -y_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, -u::Union{AbstractVector{<:Real}, Nothing}, -u_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, -t::Union{ForwardDiff.Dual, Nothing}) + dx ::AbstractVector{<:Real}, + y ::AbstractVector{<:Real}, + y_refs::AbstractVector{<:fmi2ValueReference}, + x ::AbstractVector{<:ForwardDiff.Dual}, + u ::AbstractVector{<:Real}, + u_refs::AbstractVector{<:fmi2ValueReference}, + p ::AbstractVector{<:Real}, + p_refs::AbstractVector{<:fmi2ValueReference}, + t ::Real) +@grad_from_chainrules eval!(cRef::UInt64, + dx ::AbstractVector{<:Real}, + y ::AbstractVector{<:Real}, + y_refs::AbstractVector{<:UInt32}, + x ::AbstractVector{<:ReverseDiff.TrackedReal}, + u ::AbstractVector{<:Real}, + u_refs::AbstractVector{<:UInt32}, + p ::AbstractVector{<:Real}, + p_refs::AbstractVector{<:UInt32}, + t ::Real) + +# u @ForwardDiff_frule eval!(cRef::UInt64, -dx::Union{AbstractVector{<:Real}, Nothing}, -y::Union{AbstractVector{<:Real}, Nothing}, -y_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, -u::Union{AbstractVector{<:ForwardDiff.Dual}, Nothing}, -u_refs::Union{AbstractVector{<:fmi2ValueReference}, Nothing}, -t::Union{Real, Nothing}) - -@grad_from_chainrules eval!(cRef::UInt64, -dx::Union{AbstractVector{<:ReverseDiff.TrackedReal}, Nothing}, -y::Union{AbstractVector{<:ReverseDiff.TrackedReal}, Nothing}, -y_refs::Union{AbstractVector{<:UInt32}, Nothing}, -u::Union{AbstractVector{<:ReverseDiff.TrackedReal}, Nothing}, -u_refs::Union{AbstractVector{<:UInt32}, Nothing}, -t::Union{ReverseDiff.TrackedReal, Nothing}) + dx ::AbstractVector{<:Real}, + y ::AbstractVector{<:Real}, + y_refs::AbstractVector{<:fmi2ValueReference}, + x ::AbstractVector{<:Real}, + u ::AbstractVector{<:ForwardDiff.Dual}, + u_refs::AbstractVector{<:fmi2ValueReference}, + p ::AbstractVector{<:Real}, + p_refs::AbstractVector{<:fmi2ValueReference}, + t ::Real) @grad_from_chainrules eval!(cRef::UInt64, -dx::Union{AbstractVector{<:Real}, Nothing}, -y::Union{AbstractVector{<:Real}, Nothing}, -y_refs::Union{AbstractVector{<:UInt32}, Nothing}, -u::Union{AbstractVector{<:Real}, Nothing}, -u_refs::Union{AbstractVector{<:UInt32}, Nothing}, -t::Union{ReverseDiff.TrackedReal, Nothing}) + dx ::AbstractVector{<:Real}, + y ::AbstractVector{<:Real}, + y_refs::AbstractVector{<:UInt32}, + x ::AbstractVector{<:Real}, + u ::AbstractVector{<:ReverseDiff.TrackedReal}, + u_refs::AbstractVector{<:UInt32}, + p ::AbstractVector{<:Real}, + p_refs::AbstractVector{<:UInt32}, + t ::Real) + +# p +@ForwardDiff_frule eval!(cRef::UInt64, + dx ::AbstractVector{<:Real}, + y ::AbstractVector{<:Real}, + y_refs::AbstractVector{<:fmi2ValueReference}, + x ::AbstractVector{<:Real}, + u ::AbstractVector{<:Real}, + u_refs::AbstractVector{<:fmi2ValueReference}, + p ::AbstractVector{<:ForwardDiff.Dual}, + p_refs::AbstractVector{<:fmi2ValueReference}, + t ::Real) @grad_from_chainrules eval!(cRef::UInt64, -dx::Union{AbstractVector{<:Real}, Nothing}, -y::Union{AbstractVector{<:Real}, Nothing}, -y_refs::Union{AbstractVector{<:UInt32}, Nothing}, -u::Union{AbstractVector{<:ReverseDiff.TrackedReal}, Nothing}, -u_refs::Union{AbstractVector{<:UInt32}, Nothing}, -t::Union{Real, Nothing}) - + dx ::AbstractVector{<:Real}, + y ::AbstractVector{<:Real}, + y_refs::AbstractVector{<:UInt32}, + x ::AbstractVector{<:Real}, + u ::AbstractVector{<:Real}, + u_refs::AbstractVector{<:UInt32}, + p ::AbstractVector{<:ReverseDiff.TrackedReal}, + p_refs::AbstractVector{<:UInt32}, + t ::Real) \ No newline at end of file diff --git a/src/FMIImport.jl b/src/FMIImport.jl index 8ce1754..d274598 100644 --- a/src/FMIImport.jl +++ b/src/FMIImport.jl @@ -13,15 +13,8 @@ using FMICore: fmi2Component, fmi3Instance # functions that have (currently) no better place # Receives one or an array of values and converts it into an Array{typeof(value)} (if not already). -function prepareValue(value) - if isa(value, Array) && length(size(value)) == 1 - return value - else - return [value] - end - - @assert false "prepareValue(...): Unknown dimension of structure `$dim`." -end +prepareValue(value) = [value] +prepareValue(value::AbstractVector) = value export prepareValue, prepareValueReference # wildcards for how a user can pass a fmi[X]ValueReference @@ -59,7 +52,7 @@ export fmi2SetTime, fmi2SetContinuousStates, fmi2EnterEventMode, fmi2NewDiscrete export fmi2GetDerivatives!, fmi2GetEventIndicators!, fmi2GetContinuousStates!, fmi2GetNominalsOfContinuousStates!, fmi2GetRealOutputDerivatives! # FMI2_convert.jl -export fmi2StringToValueReference, fmi2ValueReferenceToString, fmi2ModelVariablesForValueReference +export fmi2StringToValueReference, fmi2ValueReferenceToString, fmi2ModelVariablesForValueReference, fmi2DataTypeForValueReference export fmi2GetSolutionState, fmi2GetSolutionTime, fmi2GetSolutionValue, fmi2GetSolutionDerivative # FMI2_int.jl diff --git a/test/FMI2/sensitivities.jl b/test/FMI2/sensitivities.jl index e95ea15..1b132e3 100644 --- a/test/FMI2/sensitivities.jl +++ b/test/FMI2/sensitivities.jl @@ -22,6 +22,8 @@ u = [2.0] u_refs = fmu.modelDescription.inputValueReferences y = [0.0, 0.0] y_refs = fmu.modelDescription.outputValueReferences +p_refs = fmu.modelDescription.parameterValueReferences +p = zeros(length(p_refs)) dx = [0.0, 0.0] t = 0.0 @@ -32,9 +34,11 @@ function reset!(c::FMIImport.FMU2Component) c.solution.evals_∂y_∂u = 0 c.solution.evals_∂ẋ_∂t = 0 c.solution.evals_∂y_∂t = 0 + c.solution.evals_∂ẋ_∂p = 0 + c.solution.evals_∂y_∂p = 0 - @test length(dx) == 2 - @test length(y) == 2 + @test length(dx) == length(fmu.modelDescription.derivativeValueReferences) + @test length(y) == length(fmu.modelDescription.outputValueReferences) end # evaluation: set state, get state derivative (out-of-place) @@ -52,12 +56,19 @@ ydx = fmu(;x=x, u=u, u_refs=u_refs, y=y, y_refs=y_refs) # evaluation: set state and inputs, get state derivative (in-place) and outputs (in-place) ydx = fmu(;x=x, u=u, u_refs=u_refs, y=y, y_refs=y_refs, dx=dx) +# evaluation: set state and inputs, parameters, get state derivative (in-place) and outputs (in-place) +ydx = fmu(;x=x, u=u, u_refs=u_refs, y=y, y_refs=y_refs, dx=dx, p=p, p_refs=p_refs) + # known results atol= 1e-8 A = [0.0 1.0; -10.0 0.0] B = [0.0; 1.0] C = [0.0 1.0; -10.0 0.0] D = [0.0; 1.0] +E = [0.0 0.0 0.0 0.0 0.0 0.0; + 0.0 10.0 0.1 10.0 -3.0 5.0] +F = [0.0 0.0 0.0 0.0 0.0 0.0; + 0.0 10.0 0.1 10.0 -3.0 5.0] dx_t = [0.0, 0.0] y_t = [0.0, 0.0] @@ -256,5 +267,39 @@ j_rwd = ReverseDiff.jacobian(_f, [t]) @test c.solution.evals_∂y_∂t == 3 reset!(c) +# Jacobian E=∂dx/∂p +_f = _p -> fmu(;x=x, p=_p, p_refs=p_refs) +_f(p) +j_fwd = ForwardDiff.jacobian(_f, p) +j_zyg = Zygote.jacobian(_f, p)[1] +j_rwd = ReverseDiff.jacobian(_f, p) +j_smp = fmi2SampleJacobian(c, fmu.modelDescription.derivativeValueReferences, fmu.modelDescription.parameterValueReferences) +j_get = fmi2GetJacobian(c, fmu.modelDescription.derivativeValueReferences, fmu.modelDescription.parameterValueReferences) + +@test isapprox(j_fwd, E; atol=atol) +@test isapprox(j_zyg, E; atol=atol) +@test isapprox(j_rwd, E; atol=atol) +@test isapprox(j_smp, E; atol=atol) +@test isapprox(j_get, E; atol=atol) + +reset!(c) + +# Jacobian F=∂y/∂p +_f = _p -> fmu(;p=_p, p_refs=p_refs, y_refs=y_refs)[1:length(y_refs)] +_f(p) +j_fwd = ForwardDiff.jacobian(_f, p) +j_zyg = Zygote.jacobian(_f, p)[1] +j_rwd = ReverseDiff.jacobian(_f, p) +j_smp = fmi2SampleJacobian(c, fmu.modelDescription.derivativeValueReferences, fmu.modelDescription.parameterValueReferences) +j_get = fmi2GetJacobian(c, fmu.modelDescription.derivativeValueReferences, fmu.modelDescription.parameterValueReferences) + +@test isapprox(j_fwd, F; atol=atol) +@test isapprox(j_zyg, F; atol=atol) +@test isapprox(j_rwd, F; atol=atol) +@test isapprox(j_smp, F; atol=atol) +@test isapprox(j_get, F; atol=atol) + +reset!(c) + # clean up fmi2Unload(fmu) \ No newline at end of file From 3dfb218608ce6e6f47e20db9072c02198e42d645 Mon Sep 17 00:00:00 2001 From: ThummeTo <83663542+ThummeTo@users.noreply.github.com> Date: Fri, 21 Jul 2023 15:13:48 +0200 Subject: [PATCH 2/3] V0.15.7 (#94) * dependencies inc * adjustments for eigenvalue recording --- Project.toml | 8 ++++---- src/FMI2/convert.jl | 1 + src/FMI2/ext.jl | 14 +++++++++++++- src/FMI2/sens.jl | 12 ++++++++++-- test/FMI2/sensitivities.jl | 4 ++-- 5 files changed, 30 insertions(+), 9 deletions(-) diff --git a/Project.toml b/Project.toml index 9d2555c..96b6f30 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "FMIImport" uuid = "9fcbc62e-52a0-44e9-a616-1359a0008194" authors = ["TT ", "LM ", "JK "] -version = "0.15.6" +version = "0.15.7" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" @@ -14,10 +14,10 @@ SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" [compat] -ChainRulesCore = "1.15.0" +ChainRulesCore = "0.10.7, 1" EzXML = "1.1.0" -FMICore = "0.17.2" +FMICore = "0.17.3" ForwardDiffChainRules = "0.1.1" -SciMLSensitivity = "7.27.0" +SciMLSensitivity = "7.31.0" ZipFile = "0.10.0" julia = "1.6" diff --git a/src/FMI2/convert.jl b/src/FMI2/convert.jl index 51c013b..87db3a8 100644 --- a/src/FMI2/convert.jl +++ b/src/FMI2/convert.jl @@ -6,6 +6,7 @@ using ChainRulesCore: ignore_derivatives import SciMLSensitivity.ForwardDiff +# ToDo: Replace by multiple dispatch version ... # Receives one or an array of value references in an arbitrary format (see fmi2ValueReferenceFormat) and converts it into an Array{fmi2ValueReference} (if not already). function prepareValueReference(md::fmi2ModelDescription, vr::fmi2ValueReferenceFormat) tvr = typeof(vr) diff --git a/src/FMI2/ext.jl b/src/FMI2/ext.jl index 728341c..0d64955 100644 --- a/src/FMI2/ext.jl +++ b/src/FMI2/ext.jl @@ -128,8 +128,20 @@ Retrieves all the pointers of binary functions. - FMISpec2.0.2: 2.2.7 Definition of Model Variables (ModelVariables) """ -function fmi2Load(pathToFMU::String; unpackPath::Union{String, Nothing}=nothing, type::Union{Symbol, fmi2Type, Nothing}=nothing, cleanup::Bool=true, logLevel::FMULogLevel=FMULogLevelWarn) +function fmi2Load(pathToFMU::String; unpackPath::Union{String, Nothing}=nothing, type::Union{Symbol, fmi2Type, Nothing}=nothing, cleanup::Bool=true, logLevel::Union{FMULogLevel, Symbol}=FMULogLevelWarn) # Create uninitialized FMU + + if isa(logLevel, Symbol) + if logLevel == :info + logLevel = FMULogLevelInfo + elseif logLevel == :warn + logLevel = FMULogLevelWarn + elseif logLevel == :error + logLevel = FMULogLevelError + else + @assert false "Unknown logLevel symbol: `$(logLevel)`, supported are `:info`, `:warn` and `:error`." + end + end fmu = FMU2(logLevel) if startswith(pathToFMU, "http") diff --git a/src/FMI2/sens.jl b/src/FMI2/sens.jl index 39d0b73..be38e7a 100644 --- a/src/FMI2/sens.jl +++ b/src/FMI2/sens.jl @@ -13,7 +13,6 @@ import SciMLSensitivity.ReverseDiff using ChainRulesCore import ForwardDiffChainRules: @ForwardDiff_frule import SciMLSensitivity.ReverseDiff: @grad_from_chainrules -#import NonconvexUtils: @ForwardDiff_frule import ChainRulesCore: ZeroTangent, NoTangent, @thunk import FMICore: fmi2ValueReference @@ -225,7 +224,16 @@ function eval!(cRef::UInt64, end T, V, N = fd_eltypes(dx) - dx[:] = collect(ForwardDiff.Dual{T, V, N}(dx_tmp[i], ForwardDiff.partials(dx[i]) ) for i in 1:length(dx)) + dx[:] = collect(ForwardDiff.Dual{T, V, N}(dx_tmp[i], ForwardDiff.partials(dx[i]) ) for i in 1:length(dx)) + + elseif istracked(dx) + + dx_tmp = zeros(fmi2Real, length(dx)) + fmi2GetDerivatives!(c, dx_tmp) + + for e in 1:length(dx) + dx[e].value = dx_tmp[e] + end else if c.fmu.isZeroState dx[:] = [1.0] diff --git a/test/FMI2/sensitivities.jl b/test/FMI2/sensitivities.jl index 1b132e3..d537ff0 100644 --- a/test/FMI2/sensitivities.jl +++ b/test/FMI2/sensitivities.jl @@ -290,8 +290,8 @@ _f(p) j_fwd = ForwardDiff.jacobian(_f, p) j_zyg = Zygote.jacobian(_f, p)[1] j_rwd = ReverseDiff.jacobian(_f, p) -j_smp = fmi2SampleJacobian(c, fmu.modelDescription.derivativeValueReferences, fmu.modelDescription.parameterValueReferences) -j_get = fmi2GetJacobian(c, fmu.modelDescription.derivativeValueReferences, fmu.modelDescription.parameterValueReferences) +j_smp = fmi2SampleJacobian(c, fmu.modelDescription.outputValueReferences, fmu.modelDescription.parameterValueReferences) +j_get = fmi2GetJacobian(c, fmu.modelDescription.outputValueReferences, fmu.modelDescription.parameterValueReferences) @test isapprox(j_fwd, F; atol=atol) @test isapprox(j_zyg, F; atol=atol) From 7f75ca11fdfda2ac82cb6f9f504fc2d66a0d5748 Mon Sep 17 00:00:00 2001 From: ThummeTo <83663542+ThummeTo@users.noreply.github.com> Date: Wed, 2 Aug 2023 16:21:48 +0200 Subject: [PATCH 3/3] V0.15.9 (v0.15.8) (#95) * SciMLSensitivity bump version * version inc * version dec --- Project.toml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 96b6f30..061405f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "FMIImport" uuid = "9fcbc62e-52a0-44e9-a616-1359a0008194" authors = ["TT ", "LM ", "JK "] -version = "0.15.7" +version = "0.15.8" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" @@ -14,10 +14,10 @@ SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" [compat] -ChainRulesCore = "0.10.7, 1" +ChainRulesCore = "1.16" EzXML = "1.1.0" FMICore = "0.17.3" ForwardDiffChainRules = "0.1.1" -SciMLSensitivity = "7.31.0" +SciMLSensitivity = "7.35, 7.36" ZipFile = "0.10.0" julia = "1.6"