diff --git a/Project.toml b/Project.toml index 7b0997823c..57a64a0ce7 100644 --- a/Project.toml +++ b/Project.toml @@ -64,6 +64,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6" +FMI = "14a09403-18e3-468f-ad8a-74f8dda2d9ac" HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327" InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57" LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800" @@ -72,9 +73,10 @@ LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800" MTKBifurcationKitExt = "BifurcationKit" MTKChainRulesCoreExt = "ChainRulesCore" MTKDeepDiffsExt = "DeepDiffs" +MTKFMIExt = "FMI" MTKHomotopyContinuationExt = "HomotopyContinuation" -MTKLabelledArraysExt = "LabelledArrays" MTKInfiniteOptExt = "InfiniteOpt" +MTKLabelledArraysExt = "LabelledArrays" [compat] AbstractTrees = "0.3, 0.4" @@ -106,6 +108,7 @@ FindFirstFunctions = "1" ForwardDiff = "0.10.3" FunctionWrappers = "1.1" FunctionWrappersWrappers = "0.1" +FMI = "0.14" Graphs = "1.5.2" HomotopyContinuation = "2.11" InfiniteOpt = "0.5" @@ -141,7 +144,7 @@ SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2" StaticArrays = "0.10, 0.11, 0.12, 1.0" StochasticDiffEq = "6.72.1" StochasticDelayDiffEq = "1.8.1" -SymbolicIndexingInterface = "0.3.36" +SymbolicIndexingInterface = "0.3.37" SymbolicUtils = "3.7" Symbolics = "6.19" URIs = "1" diff --git a/ext/MTKFMIExt.jl b/ext/MTKFMIExt.jl new file mode 100644 index 0000000000..603ec8906c --- /dev/null +++ b/ext/MTKFMIExt.jl @@ -0,0 +1,930 @@ +module MTKFMIExt + +using ModelingToolkit +using SymbolicIndexingInterface +using ModelingToolkit: t_nounits as t, D_nounits as D +using DocStringExtensions +import ModelingToolkit as MTK +import SciMLBase +import FMI + +""" + $(TYPEDSIGNATURES) + +A utility macro for FMI.jl functions that return a status. Will terminate on +fatal statuses. Must be used as `@statuscheck FMI.fmiXFunction(...)` where +`X` should be `2` or `3`. Has an edge case for handling tuples for +`FMI.fmi2CompletedIntegratorStep`. +""" +macro statuscheck(expr) + @assert Meta.isexpr(expr, :call) + fn = expr.args[1] + @assert Meta.isexpr(fn, :.) + @assert fn.args[1] == :FMI + fnname = fn.args[2] + + instance = expr.args[2] + is_v2 = startswith("fmi2", string(fnname)) + + fmiTrue = is_v2 ? FMI.fmi2True : FMI.fmi3True + fmiStatusOK = is_v2 ? FMI.fmi2StatusOK : FMI.fmi3StatusOK + fmiStatusWarning = is_v2 ? FMI.fmi2StatusWarning : FMI.fmi3StatusWarning + fmiStatusFatal = is_v2 ? FMI.fmi2StatusFatal : FMI.fmi3StatusFatal + fmiTerminate = is_v2 ? FMI.fmi2Terminate : FMI.fmi3Terminate + fmiFreeInstance! = is_v2 ? FMI.fmi2FreeInstance! : FMI.fmi3FreeInstance! + return quote + status = $expr + fnname = $fnname + if status !== nothing && ((status isa Tuple && status[1] == $fmiTrue) || + (!(status isa Tuple) && status != $fmiStatusOK && + status != $fmiStatusWarning)) + if status != $fmiStatusFatal + $fmiTerminate(wrapper.instance) + end + $fmiFreeInstance!(wrapper.instance) + wrapper.instance = nothing + error("FMU Error in $fnname: status $status") + end + end |> esc +end + +@static if !hasmethod(FMI.getValueReferencesAndNames, Tuple{FMI.fmi3ModelDescription}) + """ + $(TYPEDSIGNATURES) + + This is type piracy, but FMI.jl is missing this implementation. It allows + `FMI.getStateValueReferencesAndNames` to work. + """ + function FMI.getValueReferencesAndNames( + md::FMI.fmi3ModelDescription; vrs = md.valueReferences) + dict = Dict{FMI.fmi3ValueReference, Array{String}}() + for vr in vrs + dict[vr] = FMI.valueReferenceToString(md, vr) + end + return dict + end +end + +""" + $(TYPEDSIGNATURES) + +A component that wraps an FMU loaded via FMI.jl. The FMI version (2 or 3) should be +provided as a `Val` to the function. Supports Model Exchange and CoSimulation FMUs. +All inputs, continuous variables and outputs must be `FMI.fmi2Real` or `FMI.fmi3Float64`. +Does not support events or discrete variables in the FMU. Does not support automatic +differentiation. Parameters of the FMU will have defaults corresponding to their initial +values in the FMU specification. All other variables will not have a default. Hierarchical +names in the FMU of the form `namespace.variable` are transformed into symbolic variables +with the name `namespace__variable`. + +# Keyword Arguments + +- `fmu`: The FMU loaded via `FMI.loadFMU`. +- `tolerance`: The tolerance to provide to the FMU. Not used for v3 FMUs since it is not + supported by FMI.jl. +- `communication_step_size`: The periodic interval at which communication with CoSimulation + FMUs will occur. Must be provided for CoSimulation FMU components. +- `reinitializealg`: The DAE initialization algorithm to use for the callback managing the + FMU. For CoSimulation FMUs whose states/outputs are used in algebraic equations of the + system, this needs to be an algorithm that will solve for the new algebraic variables. + For example, `OrdinaryDiffEqCore.BrownFullBasicInit()`. +- `type`: Either `:ME` or `:CS` depending on whether `fmu` is a Model Exchange or + CoSimulation FMU respectively. +- `name`: The name of the system. +""" +function MTK.FMIComponent(::Val{Ver}; fmu = nothing, tolerance = 1e-6, + communication_step_size = nothing, reinitializealg = SciMLBase.NoInit(), type, name) where {Ver} + if Ver != 2 && Ver != 3 + throw(ArgumentError("FMI Version must be `2` or `3`")) + end + if type == :CS && communication_step_size === nothing + throw(ArgumentError("`communication_step_size` must be specified for Co-Simulation FMUs.")) + end + # mapping from MTK variable to value reference + value_references = Dict() + # defaults + defs = Dict() + # unknowns of the system + states = [] + # differential variables of the system + # this is a subset of `states` in the case where the FMU has multiple names for + # the same value reference. + diffvars = [] + # variables that are derivatives of diffvars + dervars = [] + # observed equations + observed = Equation[] + # need to separate observed equations for duplicate derivative variables + # since they aren't included in CS FMUs + der_observed = Equation[] + + # parse states + fmi_variables_to_mtk_variables!(fmu, FMI.getStateValueReferencesAndNames(fmu), + value_references, diffvars, states, observed) + # create a symbolic variable __mtk_internal_u to pass to the relevant registered + # functions as the state vector + if isempty(diffvars) + # no differential variables + __mtk_internal_u = Float64[] + elseif type == :ME + # to avoid running into `structural_simplify` warnings about array variables + # and some unfortunate circular dependency issues, ME FMUs use an array of + # symbolics instead. This is also not worse off in performance + # because the former approach would allocate anyway. + # TODO: Can we avoid an allocation here using static arrays? + __mtk_internal_u = copy(diffvars) + elseif type == :CS + # CS FMUs do their own independent integration in a periodic callback, so their + # unknowns are discrete variables in the `ODESystem`. A default of `missing` allows + # them to be solved for during initialization. + @parameters __mtk_internal_u(t)[1:length(diffvars)]=missing [guess = diffvars] + push!(observed, __mtk_internal_u ~ copy(diffvars)) + end + + # parse derivatives of states + # the variables passed to `postprocess_variable` haven't been differentiated yet, so they + # should match one variable in states. That's the one this is the derivative of, and we + # keep track of this ordering + derivative_order = [] + function derivative_order_postprocess(var) + idx = findfirst(isequal(var), states) + idx === nothing || push!(derivative_order, states[idx]) + return var + end + fmi_variables_to_mtk_variables!( + fmu, FMI.getDerivateValueReferencesAndNames(fmu), value_references, dervars, + states, der_observed; postprocess_variable = derivative_order_postprocess) + @assert length(derivative_order) == length(dervars) + + # parse the inputs to the FMU + inputs = [] + fmi_variables_to_mtk_variables!(fmu, FMI.getInputValueReferencesAndNames(fmu), + value_references, inputs, states, observed; postprocess_variable = v -> MTK.setinput( + v, true)) + # create a symbolic variable for the input buffer + __mtk_internal_x = copy(inputs) + if isempty(__mtk_internal_x) + __mtk_internal_x = Float64[] + end + + # parse the outputs of the FMU + outputs = [] + fmi_variables_to_mtk_variables!(fmu, FMI.getOutputValueReferencesAndNames(fmu), + value_references, outputs, states, observed; postprocess_variable = v -> MTK.setoutput( + v, true)) + # create the output buffer. This is only required for CoSimulation to pass it to + # the callback affect + if type == :CS + if isempty(outputs) + __mtk_internal_o = Float64[] + else + @parameters __mtk_internal_o(t)[1:length(outputs)]=missing [guess = zeros(length(outputs))] + push!(observed, __mtk_internal_o ~ outputs) + end + end + + # parse the parameters + params = [] + # multiple names for the same parameter are treated as parameter dependencies. + parameter_dependencies = Equation[] + fmi_variables_to_mtk_variables!( + fmu, FMI.getParameterValueReferencesAndNames(fmu), value_references, + params, [], parameter_dependencies, defs; parameters = true) + # create a symbolic variable for the parameter buffer + __mtk_internal_p = copy(params) + if isempty(__mtk_internal_p) + __mtk_internal_p = Float64[] + end + + derivative_value_references = UInt32[value_references[var] for var in dervars] + state_value_references = UInt32[value_references[var] for var in diffvars] + output_value_references = UInt32[value_references[var] for var in outputs] + input_value_references = UInt32[value_references[var] for var in inputs] + param_value_references = UInt32[value_references[var] for var in params] + + # create a parameter for the instance wrapper + # this manages the creation and deallocation of FMU instances + buffer_length = length(diffvars) + length(outputs) + if Ver == 2 + @parameters (wrapper::FMI2InstanceWrapper)(..)[1:buffer_length] = FMI2InstanceWrapper( + fmu, derivative_value_references, state_value_references, output_value_references, + param_value_references, input_value_references, tolerance) + else + @parameters (wrapper::FMI3InstanceWrapper)(..)[1:buffer_length] = FMI3InstanceWrapper( + fmu, derivative_value_references, state_value_references, + output_value_references, param_value_references, input_value_references) + end + + # any additional initialization equations for the system + initialization_eqs = Equation[] + + if type == :ME + # the wrapper is a callable struct which returns the state derivative and + # output values + # symbolic expression for calling the wrapper + call_expr = wrapper(__mtk_internal_u, __mtk_internal_x, __mtk_internal_p, t) + + # differential and observed equations + diffeqs = Equation[] + for (i, var) in enumerate([dervars; outputs]) + push!(diffeqs, var ~ call_expr[i]) + end + for (var, dervar) in zip(derivative_order, dervars) + push!(diffeqs, D(var) ~ dervar) + end + + # instance management callback which deallocates the instance when + # necessary and notifies the FMU of completed integrator steps + finalize_affect = MTK.FunctionalAffect(fmiFinalize!, [], [wrapper], []) + step_affect = MTK.FunctionalAffect(Returns(nothing), [], [], []) + instance_management_callback = MTK.SymbolicDiscreteCallback( + (t != t - 1), step_affect; finalize = finalize_affect, reinitializealg = reinitializealg) + + push!(params, wrapper) + append!(observed, der_observed) + elseif type == :CS + _functor = if Ver == 2 + FMI2CSFunctor(state_value_references, output_value_references) + else + FMI3CSFunctor(state_value_references, output_value_references) + end + @parameters (functor::(typeof(_functor)))(..)[1:(length(__mtk_internal_u) + length(__mtk_internal_o))] = _functor + # for co-simulation, we need to ensure the output buffer is solved for + # during initialization + for (i, x) in enumerate(collect(__mtk_internal_o)) + push!(initialization_eqs, + x ~ functor( + wrapper, __mtk_internal_u, __mtk_internal_x, __mtk_internal_p, t)[i]) + end + + diffeqs = Equation[] + + # use `ImperativeAffect` for instance management here + cb_observed = (; inputs = __mtk_internal_x, params = copy(params), + t, wrapper, dt = communication_step_size) + cb_modified = (;) + # modify the outputs if present + if symbolic_type(__mtk_internal_o) != NotSymbolic() + cb_modified = (cb_modified..., outputs = __mtk_internal_o) + end + # modify the continuous state if present + if symbolic_type(__mtk_internal_u) != NotSymbolic() + cb_modified = (cb_modified..., states = __mtk_internal_u) + end + initialize_affect = MTK.ImperativeAffect(fmiCSInitialize!; observed = cb_observed, + modified = cb_modified, ctx = _functor) + finalize_affect = MTK.FunctionalAffect(fmiFinalize!, [], [wrapper], []) + # the callback affect performs the stepping + step_affect = MTK.ImperativeAffect( + fmiCSStep!; observed = cb_observed, modified = cb_modified, ctx = _functor) + instance_management_callback = MTK.SymbolicDiscreteCallback( + communication_step_size, step_affect; initialize = initialize_affect, + finalize = finalize_affect, reinitializealg = reinitializealg + ) + + # guarded in case there are no outputs/states and the variable is `[]`. + symbolic_type(__mtk_internal_o) == NotSymbolic() || push!(params, __mtk_internal_o) + symbolic_type(__mtk_internal_u) == NotSymbolic() || push!(params, __mtk_internal_u) + + push!(params, wrapper, functor) + end + + eqs = [observed; diffeqs] + return ODESystem(eqs, t, states, params; parameter_dependencies, defaults = defs, + discrete_events = [instance_management_callback], name, initialization_eqs) +end + +""" + $(TYPEDSIGNATURES) + +A utility function which accepts an FMU `fmu` and a mapping from value reference to a +list of associated names `varmap`. A symbolic variable is created for each name. The +associated value reference is kept track of in `value_references`. In case there are +multiple names for a value reference, the symbolic variable for the first name is pushed +to `truevars`. All of the created symbolic variables are pushed to `allvars`. Observed +equations equating identical variables are pushed to `obseqs`. `defs` is a dictionary of +defaults. + +# Keyword Arguments +- `parameters`: A boolean indicating whether to use `@parameters` for the symbolic + variables instead of `@variables`. +- `postprocess_variable`: A function applied to each created variable that should + return the updated variable. This is useful to add metadata to variables. +""" +function fmi_variables_to_mtk_variables!( + fmu::Union{FMI.FMU2, FMI.FMU3}, varmap::AbstractDict, + value_references::AbstractDict, truevars, allvars, + obseqs, defs = Dict(); parameters = false, postprocess_variable = identity) + for (valRef, varnames) in varmap + stateT = FMI.dataTypeForValueReference(fmu, valRef) + snames = Symbol[] + ders = Int[] + for name in varnames + sname, der = parseFMIVariableName(name) + push!(snames, sname) + push!(ders, der) + end + if parameters + vars = [postprocess_variable(MTK.unwrap(only(@parameters $sname::stateT))) + for sname in snames] + else + vars = [postprocess_variable(MTK.unwrap(only(@variables $sname(t)::stateT))) + for sname in snames] + end + for i in eachindex(vars) + der = ders[i] + vars[i] = MTK.unwrap(vars[i]) + for j in 1:der + vars[i] = D(vars[i]) + end + vars[i] = MTK.default_toterm(vars[i]) + end + for i in eachindex(vars) + if i == 1 + push!(truevars, vars[i]) + else + push!(obseqs, vars[i] ~ vars[1]) + end + value_references[vars[i]] = valRef + end + append!(allvars, vars) + defval = FMI.getStartValue(fmu, valRef) + defs[vars[1]] = defval + end +end + +""" + $(TYPEDSIGNATURES) + +Parse the string name of an FMI variable into a `Symbol` name for the corresponding +MTK variable. Return the `Symbol` name and the number of times it is differentiated. +""" +function parseFMIVariableName(name::AbstractString) + name = replace(name, "." => "__") + der = 0 + if startswith(name, "der(") + idx = findfirst(',', name) + if idx === nothing + name = @view name[5:(end - 1)] + der = 1 + else + der = parse(Int, @view name[(idx + 1):(end - 1)]) + name = @view name[5:(idx - 1)] + end + end + return Symbol(name), der +end + +""" + $(TYPEDEF) + +A struct which manages instance creation and deallocation for v2 FMUs. + +# Fields + +$(TYPEDFIELDS) +""" +mutable struct FMI2InstanceWrapper + """ + The FMU from `FMI.loadFMU`. + """ + const fmu::FMI.FMU2 + """ + The value references for derivatives of states of the FMU, in the order that the + caller expects them to be returned when calling this struct. + """ + const derivative_value_references::Vector{FMI.fmi2ValueReference} + """ + The value references for the states of the FMU. + """ + const state_value_references::Vector{FMI.fmi2ValueReference} + """ + The value references for outputs of the FMU, in the order that the caller expects + them to be returned when calling this struct. + """ + const output_value_references::Vector{FMI.fmi2ValueReference} + """ + The parameter value references. These should be in the same order as the parameter + vector passed to functions involving this wrapper. + """ + const param_value_references::Vector{FMI.fmi2ValueReference} + """ + The input value references. These should be in the same order as the inputs passed + to functions involving this wrapper. + """ + const input_value_references::Vector{FMI.fmi2ValueReference} + """ + The tolerance with which to setup the FMU instance. + """ + const tolerance::FMI.fmi2Real + """ + The FMU instance, if present, and `nothing` otherwise. + """ + instance::Union{FMI.FMU2Component{FMI.FMU2}, Nothing} +end + +""" + $(TYPEDSIGNATURES) + +Create an `FMI2InstanceWrapper` with no instance. +""" +function FMI2InstanceWrapper(fmu, ders, states, outputs, params, inputs, tolerance) + FMI2InstanceWrapper(fmu, ders, states, outputs, params, inputs, tolerance, nothing) +end + +Base.nameof(::FMI2InstanceWrapper) = :FMI2InstanceWrapper + +""" + $(TYPEDSIGNATURES) + +Common functionality for creating an instance of a v2 FMU. Does not check if +`wrapper.instance` is already present, and overwrites the existing value with +a new instance. `inputs` should be in the order of `wrapper.input_value_references`. +`params` should be in the order of `wrapper.param_value_references`. `t` is the current +time. Returns the created instance, which is also stored in `wrapper.instance`. +""" +function get_instance_common!(wrapper::FMI2InstanceWrapper, inputs, params, t) + wrapper.instance = FMI.fmi2Instantiate!(wrapper.fmu)::FMI.FMU2Component + if !isempty(inputs) + @statuscheck FMI.fmi2SetReal(wrapper.instance, wrapper.input_value_references, + Csize_t(length(wrapper.param_value_references)), inputs) + end + if !isempty(params) + @statuscheck FMI.fmi2SetReal(wrapper.instance, wrapper.param_value_references, + Csize_t(length(wrapper.param_value_references)), params) + end + @statuscheck FMI.fmi2SetupExperiment( + wrapper.instance, FMI.fmi2True, wrapper.tolerance, t, FMI.fmi2False, t) + @statuscheck FMI.fmi2EnterInitializationMode(wrapper.instance) + return wrapper.instance +end + +""" + $(TYPEDSIGNATURES) + +Create an instance of a Model Exchange FMU. Use the existing instance in `wrapper` if +present and create a new one otherwise. Return the instance. + +See `get_instance_common!` for a description of the arguments. +""" +function get_instance_ME!(wrapper::FMI2InstanceWrapper, inputs, params, t) + if wrapper.instance === nothing + get_instance_common!(wrapper, inputs, params, t) + @statuscheck FMI.fmi2ExitInitializationMode(wrapper.instance) + eventInfo = FMI.fmi2NewDiscreteStates(wrapper.instance) + @assert eventInfo.newDiscreteStatesNeeded == FMI.fmi2False + # TODO: Support FMU events + @statuscheck FMI.fmi2EnterContinuousTimeMode(wrapper.instance) + end + + return wrapper.instance +end + +""" + $(TYPEDSIGNATURES) + +Create an instance of a CoSimulation FMU. Use the existing instance in `wrapper` if +present and create a new one otherwise. Return the instance. + +See `get_instance_common!` for a description of the arguments. +""" +function get_instance_CS!(wrapper::FMI2InstanceWrapper, states, inputs, params, t) + if wrapper.instance === nothing + get_instance_common!(wrapper, inputs, params, t) + if !isempty(states) + @statuscheck FMI.fmi2SetReal(wrapper.instance, wrapper.state_value_references, + Csize_t(length(wrapper.state_value_references)), states) + end + @statuscheck FMI.fmi2ExitInitializationMode(wrapper.instance) + end + return wrapper.instance +end + +""" + $(TYPEDSIGNATURES) + +Call `fmiXCompletedIntegratorStep` with `noSetFMUStatePriorToCurrentPoint` as false. +""" +function partiallyCompleteIntegratorStep(wrapper::FMI2InstanceWrapper) + @statuscheck FMI.fmi2CompletedIntegratorStep(wrapper.instance, FMI.fmi2False) +end + +""" + $(TYPEDSIGNATURES) + +If `wrapper.instance !== nothing`, terminate and free the instance. Also set +`wrapper.instance` to `nothing`. +""" +function reset_instance!(wrapper::FMI2InstanceWrapper) + wrapper.instance === nothing && return + FMI.fmi2Terminate(wrapper.instance) + FMI.fmi2FreeInstance!(wrapper.instance) + wrapper.instance = nothing +end + +""" + $(TYPEDEF) + +A struct which manages instance creation and deallocation for v3 FMUs. + +# Fields + +$(TYPEDFIELDS) +""" +mutable struct FMI3InstanceWrapper + """ + The FMU from `FMI.loadFMU`. + """ + const fmu::FMI.FMU3 + """ + The value references for derivatives of states of the FMU, in the order that the + caller expects them to be returned when calling this struct. + """ + const derivative_value_references::Vector{FMI.fmi3ValueReference} + const state_value_references::Vector{FMI.fmi3ValueReference} + """ + The value references for outputs of the FMU, in the order that the caller expects + them to be returned when calling this struct. + """ + const output_value_references::Vector{FMI.fmi3ValueReference} + """ + The parameter value references. These should be in the same order as the parameter + vector passed to functions involving this wrapper. + """ + const param_value_references::Vector{FMI.fmi3ValueReference} + """ + The input value references. These should be in the same order as the inputs passed + to functions involving this wrapper. + """ + const input_value_references::Vector{FMI.fmi3ValueReference} + """ + The FMU instance, if present, and `nothing` otherwise. + """ + instance::Union{FMI.FMU3Instance{FMI.FMU3}, Nothing} +end + +""" + $(TYPEDSIGNATURES) + +Create an `FMI3InstanceWrapper` with no instance. +""" +function FMI3InstanceWrapper(fmu, ders, states, outputs, params, inputs) + FMI3InstanceWrapper(fmu, ders, states, outputs, params, inputs, nothing) +end + +Base.nameof(::FMI3InstanceWrapper) = :FMI3InstanceWrapper + +""" + $(TYPEDSIGNATURES) + +Common functionality for creating an instance of a v3 FMU. Since v3 FMUs need to be +instantiated differently depending on the type, this assumes `wrapper.instance` is a +freshly instantiated FMU which needs to be initialized. `inputs` should be in the order +of `wrapper.input_value_references`. `params` should be in the order of +`wrapper.param_value_references`. `t` is the current time. Returns `wrapper.instance`. +""" +function get_instance_common!(wrapper::FMI3InstanceWrapper, inputs, params, t) + if !isempty(params) + @statuscheck FMI.fmi3SetFloat64(wrapper.instance, wrapper.param_value_references, + params) + end + @statuscheck FMI.fmi3EnterInitializationMode( + wrapper.instance, FMI.fmi3False, zero(FMI.fmi3Float64), t, FMI.fmi3False, t) + if !isempty(inputs) + @statuscheck FMI.fmi3SetFloat64( + wrapper.instance, wrapper.input_value_references, inputs) + end + + return wrapper.instance +end + +""" + $(TYPEDSIGNATURES) + +Create an instance of a Model Exchange FMU. Use the existing instance in `wrapper` if +present and create a new one otherwise. Return the instance. + +See `get_instance_common!` for a description of the arguments. +""" +function get_instance_ME!(wrapper::FMI3InstanceWrapper, inputs, params, t) + if wrapper.instance === nothing + wrapper.instance = FMI.fmi3InstantiateModelExchange!(wrapper.fmu)::FMI.FMU3Instance + get_instance_common!(wrapper, inputs, params, t) + @statuscheck FMI.fmi3ExitInitializationMode(wrapper.instance) + eventInfo = FMI.fmi3UpdateDiscreteStates(wrapper.instance) + @assert eventInfo[1] == FMI.fmi2False + # TODO: Support FMU events + @statuscheck FMI.fmi3EnterContinuousTimeMode(wrapper.instance) + end + + return wrapper.instance +end + +""" + $(TYPEDSIGNATURES) + +Create an instance of a CoSimulation FMU. Use the existing instance in `wrapper` if +present and create a new one otherwise. Return the instance. + +See `get_instance_common!` for a description of the arguments. +""" +function get_instance_CS!(wrapper::FMI3InstanceWrapper, states, inputs, params, t) + if wrapper.instance === nothing + wrapper.instance = FMI.fmi3InstantiateCoSimulation!( + wrapper.fmu; eventModeUsed = false)::FMI.FMU3Instance + get_instance_common!(wrapper, inputs, params, t) + if !isempty(states) + @statuscheck FMI.fmi3SetFloat64( + wrapper.instance, wrapper.state_value_references, states) + end + @statuscheck FMI.fmi3ExitInitializationMode(wrapper.instance) + end + return wrapper.instance +end + +""" + $(TYPEDSIGNATURES) +""" +function partiallyCompleteIntegratorStep(wrapper::FMI3InstanceWrapper) + enterEventMode = Ref(FMI.fmi3False) + terminateSimulation = Ref(FMI.fmi3False) + @statuscheck FMI.fmi3CompletedIntegratorStep!( + wrapper.instance, FMI.fmi3False, enterEventMode, terminateSimulation) + @assert enterEventMode[] == FMI.fmi3False + @assert terminateSimulation[] == FMI.fmi3False +end + +""" + $(TYPEDSIGNATURES) +""" +function reset_instance!(wrapper::FMI3InstanceWrapper) + wrapper.instance === nothing && return + FMI.fmi3Terminate(wrapper.instance) + FMI.fmi3FreeInstance!(wrapper.instance) + wrapper.instance = nothing +end + +@register_array_symbolic (fn::FMI2InstanceWrapper)( + states::Vector{<:Real}, inputs::Vector{<:Real}, params::Vector{<:Real}, t::Real) begin + size = (length(states) + length(fn.output_value_references),) + eltype = eltype(states) + ndims = 1 +end + +@register_array_symbolic (fn::FMI3InstanceWrapper)( + wrapper::FMI3InstanceWrapper, states::Vector{<:Real}, + inputs::Vector{<:Real}, params::Vector{<:Real}, t::Real) begin + size = (length(states) + length(fn.output_value_references),) + eltype = eltype(states) + ndims = 1 +end + +""" + $(TYPEDSIGNATURES) + +Update the internal state of the ME FMU and return a vector of updated values +for continuous state derivatives and output variables respectively. Needs to be a +callable struct to enable symbolic registration with an inferred return size. +""" +function (wrapper::Union{FMI2InstanceWrapper, FMI3InstanceWrapper})( + states, inputs, params, t) + instance = get_instance_ME!(wrapper, inputs, params, t) + + # TODO: Find a way to do this without allocating. We can't pass a view to these + # functions. + states_buffer = zeros(length(states)) + outputs_buffer = zeros(length(wrapper.output_value_references)) + # Defined in FMIBase.jl/src/eval.jl + # Doesn't seem to be documented, but somehow this is the only way to + # propagate inputs to the FMU consistently. I have no idea why. + instance(; x = states, u = inputs, u_refs = wrapper.input_value_references, + p = params, p_refs = wrapper.param_value_references, t = t) + # the spec requires completing the step before getting updated derivative/output values + partiallyCompleteIntegratorStep(wrapper) + instance(; dx = states_buffer, dx_refs = wrapper.derivative_value_references, + y = outputs_buffer, y_refs = wrapper.output_value_references) + return [states_buffer; outputs_buffer] +end + +""" + $(TYPEDSIGNATURES) + +An affect function for use inside a `FunctionalAffect`. This should be triggered at the +end of the solve, regardless of whether it succeeded or failed. Expects `p` to be a +1-length array containing the index of the instance wrapper (`FMI2InstanceWrapper` or +`FMI3InstanceWrapper`) in the parameter object. +""" +function fmiFinalize!(integrator, u, p, ctx) + wrapper_idx = p[1] + wrapper = integrator.ps[wrapper_idx] + reset_instance!(wrapper) +end + +""" + $(TYPEDEF) + +A callable struct useful for initializing v2 CoSimulation FMUs. When called, updates the +internal state of the FMU and gets updated values for output variables. + +# Fields + +$(TYPEDFIELDS) +""" +struct FMI2CSFunctor + """ + The value references of state variables in the FMU. + """ + state_value_references::Vector{FMI.fmi2ValueReference} + """ + The value references of output variables in the FMU. + """ + output_value_references::Vector{FMI.fmi2ValueReference} +end + +function (fn::FMI2CSFunctor)(wrapper::FMI2InstanceWrapper, states, inputs, params, t) + states = states isa SubArray ? copy(states) : states + inputs = inputs isa SubArray ? copy(inputs) : inputs + params = params isa SubArray ? copy(params) : params + instance = get_instance_CS!(wrapper, states, inputs, params, t) + if isempty(fn.output_value_references) + return eltype(states)[] + else + return FMI.fmi2GetReal(instance, fn.output_value_references) + end +end + +@register_array_symbolic (fn::FMI2CSFunctor)( + wrapper::FMI2InstanceWrapper, states::Vector{<:Real}, + inputs::Vector{<:Real}, params::Vector{<:Real}, t::Real) begin + size = (length(states) + length(fn.output_value_references),) + eltype = eltype(states) + ndims = 1 +end + +""" + $(TYPEDSIGNATURES) + +An affect function designed for use with `ImperativeAffect`. Should be triggered during +callback initialization. `m` should contain the key `:states` with the value being the +state vector if the FMU has continuous states. `m` should contain the key `:outputs` with +the value being the output vector if the FMU has output variables. `o` should contain the +`:inputs`, `:params`, `:t` and `:wrapper` where the latter contains the `FMI2InstanceWrapper`. + +Initializes the FMU. Only for use with CoSimulation FMUs. +""" +function fmiCSInitialize!(m, o, ctx::FMI2CSFunctor, integrator) + states = isdefined(m, :states) ? m.states : () + inputs = o.inputs + params = o.params + t = o.t + wrapper = o.wrapper + if wrapper.instance !== nothing + reset_instance!(wrapper) + end + + instance = get_instance_CS!(wrapper, states, inputs, params, t) + if isdefined(m, :states) + @statuscheck FMI.fmi2GetReal!(instance, ctx.state_value_references, m.states) + end + if isdefined(m, :outputs) + @statuscheck FMI.fmi2GetReal!(instance, ctx.output_value_references, m.outputs) + end + + return m +end + +""" + $(TYPEDSIGNATURES) + +An affect function designed for use with `ImperativeAffect`. Should be triggered +periodically to communicte with the CoSimulation FMU. Has the same requirements as +`fmiCSInitialize!` for `m` and `o`, with the addition that `o` should have a key +`:dt` with the value being the communication step size. +""" +function fmiCSStep!(m, o, ctx::FMI2CSFunctor, integrator) + wrapper = o.wrapper + states = isdefined(m, :states) ? m.states : () + inputs = o.inputs + params = o.params + t = o.t + dt = o.dt + + instance = get_instance_CS!(wrapper, states, inputs, params, integrator.t) + if !isempty(inputs) + FMI.fmi2SetReal( + instance, wrapper.input_value_references, Csize_t(length(inputs)), inputs) + end + @statuscheck FMI.fmi2DoStep(instance, integrator.t - dt, dt, FMI.fmi2True) + + if isdefined(m, :states) + @statuscheck FMI.fmi2GetReal!(instance, ctx.state_value_references, m.states) + end + if isdefined(m, :outputs) + @statuscheck FMI.fmi2GetReal!(instance, ctx.output_value_references, m.outputs) + end + + return m +end + +""" + $(TYPEDEF) + +A callable struct useful for initializing v3 CoSimulation FMUs. When called, updates the +internal state of the FMU and gets updated values for output variables. + +# Fields + +$(TYPEDFIELDS) +""" +struct FMI3CSFunctor + """ + The value references of state variables in the FMU. + """ + state_value_references::Vector{FMI.fmi3ValueReference} + """ + The value references of output variables in the FMU. + """ + output_value_references::Vector{FMI.fmi3ValueReference} +end + +function (fn::FMI3CSFunctor)(wrapper::FMI3InstanceWrapper, states, inputs, params, t) + states = states isa SubArray ? copy(states) : states + inputs = inputs isa SubArray ? copy(inputs) : inputs + params = params isa SubArray ? copy(params) : params + instance = get_instance_CS!(wrapper, states, inputs, params, t) + + if isempty(fn.output_value_references) + return eltype(states)[] + else + return FMI.fmi3GetFloat64(instance, fn.output_value_references) + end +end + +@register_array_symbolic (fn::FMI3CSFunctor)( + wrapper::FMI3InstanceWrapper, states::Vector{<:Real}, + inputs::Vector{<:Real}, params::Vector{<:Real}, t::Real) begin + size = (length(states) + length(fn.output_value_references),) + eltype = eltype(states) + ndims = 1 +end + +""" + $(TYPEDSIGNATURES) +""" +function fmiCSInitialize!(m, o, ctx::FMI3CSFunctor, integrator) + states = isdefined(m, :states) ? m.states : () + inputs = o.inputs + params = o.params + t = o.t + wrapper = o.wrapper + if wrapper.instance !== nothing + reset_instance!(wrapper) + end + instance = get_instance_CS!(wrapper, states, inputs, params, t) + if isdefined(m, :states) + @statuscheck FMI.fmi3GetFloat64!(instance, ctx.state_value_references, m.states) + end + if isdefined(m, :outputs) + @statuscheck FMI.fmi3GetFloat64!(instance, ctx.output_value_references, m.outputs) + end + + return m +end + +""" + $(TYPEDSIGNATURES) +""" +function fmiCSStep!(m, o, ctx::FMI3CSFunctor, integrator) + wrapper = o.wrapper + states = isdefined(m, :states) ? m.states : () + inputs = o.inputs + params = o.params + t = o.t + dt = o.dt + + instance = get_instance_CS!(wrapper, states, inputs, params, integrator.t) + if !isempty(inputs) + FMI.fmi3SetFloat64(instance, wrapper.input_value_references, inputs) + end + eventEncountered = Ref(FMI.fmi3False) + terminateSimulation = Ref(FMI.fmi3False) + earlyReturn = Ref(FMI.fmi3False) + lastSuccessfulTime = Ref(zero(FMI.fmi3Float64)) + @statuscheck FMI.fmi3DoStep!( + instance, integrator.t - dt, dt, FMI.fmi3True, eventEncountered, + terminateSimulation, earlyReturn, lastSuccessfulTime) + @assert eventEncountered[] == FMI.fmi3False + @assert terminateSimulation[] == FMI.fmi3False + @assert earlyReturn[] == FMI.fmi3False + + if isdefined(m, :states) + @statuscheck FMI.fmi3GetFloat64!(instance, ctx.state_value_references, m.states) + end + if isdefined(m, :outputs) + @statuscheck FMI.fmi3GetFloat64!(instance, ctx.output_value_references, m.outputs) + end + + return m +end + +end # module diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index 10aba4d8a9..cde6261665 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -291,4 +291,6 @@ export MTKParameters, reorder_dimension_by_tunables!, reorder_dimension_by_tunab export HomotopyContinuationProblem +function FMIComponent end + end # module diff --git a/src/parameters.jl b/src/parameters.jl index ced7767718..91121b7cbb 100644 --- a/src/parameters.jl +++ b/src/parameters.jl @@ -33,7 +33,12 @@ end function getcalledparameter(x) x = unwrap(x) - return getmetadata(x, CallWithParent) + # `parent` is a `CallWithMetadata` with the correct metadata, + # but no namespacing. `operation(x)` has the correct namespacing, + # but is not a `CallWithMetadata` and doesn't have any metadata. + # This approach combines both. + parent = getmetadata(x, CallWithParent) + return CallWithMetadata(operation(x), metadata(parent)) end """ diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl index 8161a9572b..f1a152b8e7 100644 --- a/src/structural_transformation/symbolics_tearing.jl +++ b/src/structural_transformation/symbolics_tearing.jl @@ -593,7 +593,7 @@ function tearing_reassemble(state::TearingState, var_eq_matching, @set! sys.unknowns = unknowns obs, subeqs, deps = cse_and_array_hacks( - obs, subeqs, unknowns, neweqs; cse = cse_hack, array = array_hack) + sys, obs, subeqs, unknowns, neweqs; cse = cse_hack, array = array_hack) @set! sys.eqs = neweqs @set! sys.observed = obs @@ -627,7 +627,7 @@ if all `p[i]` are present and the unscalarized form is used in any equation (obs not) we first count the number of times the scalarized form of each observed variable occurs in observed equations (and unknowns if it's split). """ -function cse_and_array_hacks(obs, subeqs, unknowns, neweqs; cse = true, array = true) +function cse_and_array_hacks(sys, obs, subeqs, unknowns, neweqs; cse = true, array = true) # HACK 1 # mapping of rhs to temporary CSE variable # `f(...) => tmpvar` in above example @@ -696,6 +696,7 @@ function cse_and_array_hacks(obs, subeqs, unknowns, neweqs; cse = true, array = tempvar; T = Symbolics.symtype(rhs_arr))) tempvar = setmetadata( tempvar, Symbolics.ArrayShapeCtx, Symbolics.shape(rhs_arr)) + vars!(all_vars, rhs_arr) tempeq = tempvar ~ rhs_arr rhs_to_tempvar[rhs_arr] = tempvar push!(obs, tempeq) @@ -718,12 +719,16 @@ function cse_and_array_hacks(obs, subeqs, unknowns, neweqs; cse = true, array = Symbolics.shape(sym) != Symbolics.Unknown() || continue arg1 = arguments(sym)[1] cnt = get(arr_obs_occurrences, arg1, 0) - cnt == 0 && continue arr_obs_occurrences[arg1] = cnt + 1 end for eq in neweqs vars!(all_vars, eq.rhs) end + + # also count unscalarized variables used in callbacks + for ev in Iterators.flatten((continuous_events(sys), discrete_events(sys))) + vars!(all_vars, ev) + end obs_arr_eqs = Equation[] for (arrvar, cnt) in arr_obs_occurrences cnt == length(arrvar) || continue @@ -737,7 +742,9 @@ function cse_and_array_hacks(obs, subeqs, unknowns, neweqs; cse = true, array = # try to `create_array(OffsetArray{...}, ...)` which errors. # `term(Origin(firstind), scal)` doesn't retain the `symtype` and `size` # of `scal`. - push!(obs_arr_eqs, arrvar ~ change_origin(Origin(firstind), scal)) + rhs = scal + rhs = change_origin(firstind, rhs) + push!(obs_arr_eqs, arrvar ~ rhs) end append!(obs, obs_arr_eqs) append!(subeqs, obs_arr_eqs) @@ -764,10 +771,10 @@ getindex_wrapper(x, i) = x[i...] # PART OF HACK 2 function change_origin(origin, arr) - return origin(arr) + return Origin(origin)(arr) end -@register_array_symbolic change_origin(origin::Origin, arr::AbstractArray) begin +@register_array_symbolic change_origin(origin::Any, arr::AbstractArray) begin size = size(arr) eltype = eltype(arr) ndims = ndims(arr) diff --git a/src/structural_transformation/utils.jl b/src/structural_transformation/utils.jl index bd24a1d017..0de7199a23 100644 --- a/src/structural_transformation/utils.jl +++ b/src/structural_transformation/utils.jl @@ -218,11 +218,18 @@ function find_eq_solvables!(state::TearingState, ieq, to_rm = Int[], coeffs = no all_int_vars = true coeffs === nothing || empty!(coeffs) empty!(to_rm) + + vars_buffer = Set() for j in 𝑠neighbors(graph, ieq) var = fullvars[j] isirreducible(var) && (all_int_vars = false; continue) a, b, islinear = linear_expansion(term, var) a, b = unwrap(a), unwrap(b) + vars!(vars_buffer, b) + if islinear && isequal(a, 0) && var in vars_buffer + islinear = false + end + empty!(vars_buffer) islinear || (all_int_vars = false; continue) a = ModelingToolkit.fold_constants(a) b = ModelingToolkit.fold_constants(b) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index eaf31a9c5f..492b5dac8c 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -77,6 +77,13 @@ function has_functional_affect(cb) (affects(cb) isa FunctionalAffect || affects(cb) isa ImperativeAffect) end +function vars!(vars, aff::FunctionalAffect; op = Differential) + for var in Iterators.flatten((unknowns(aff), parameters(aff), discretes(aff))) + vars!(vars, var) + end + return vars +end + #################################### continuous events ##################################### const NULL_AFFECT = Equation[] @@ -333,6 +340,22 @@ function continuous_events(sys::AbstractSystem) filter(!isempty, cbs) end +function vars!(vars, cb::SymbolicContinuousCallback; op = Differential) + for eq in equations(cb) + vars!(vars, eq; op) + end + for aff in (affects(cb), affect_negs(cb), initialize_affects(cb), finalize_affects(cb)) + if aff isa Vector{Equation} + for eq in aff + vars!(vars, eq; op) + end + elseif aff !== nothing + vars!(vars, aff; op) + end + end + return vars +end + #################################### discrete events ##################################### struct SymbolicDiscreteCallback @@ -469,6 +492,28 @@ function discrete_events(sys::AbstractSystem) cbs end +function vars!(vars, cb::SymbolicDiscreteCallback; op = Differential) + if symbolic_type(cb.condition) == NotSymbolic + if cb.condition isa AbstractArray + for eq in cb.condition + vars!(vars, eq; op) + end + end + else + vars!(vars, cb.condition; op) + end + for aff in (cb.affects, cb.initialize, cb.finalize) + if aff isa Vector{Equation} + for eq in aff + vars!(vars, eq; op) + end + elseif aff !== nothing + vars!(vars, aff; op) + end + end + return vars +end + ################################# compilation functions #################################### # handles ensuring that affect! functions work with integrator arguments diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 20293068a9..7d1c00a6a6 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -1369,11 +1369,13 @@ function InitializationProblem{iip, specialize}(sys::AbstractSystem, for sym in unknowns(isys) haskey(fullmap, sym) || continue symbolic_type(fullmap[sym]) == NotSymbolic() || continue + is_array_of_symbolics(fullmap[sym]) && continue u0T = promote_type(u0T, typeof(fullmap[sym])) end for eq in observed(isys) haskey(fullmap, eq.lhs) || continue symbolic_type(fullmap[eq.lhs]) == NotSymbolic() || continue + is_array_of_symbolics(fullmap[eq.lhs]) && continue u0T = promote_type(u0T, typeof(fullmap[eq.lhs])) end if u0T != Union{} diff --git a/src/systems/imperative_affect.jl b/src/systems/imperative_affect.jl index 2f489913ad..aee95fd051 100644 --- a/src/systems/imperative_affect.jl +++ b/src/systems/imperative_affect.jl @@ -38,7 +38,7 @@ in the returned tuple, in which case the associated field will not be updated. skip_checks::Bool end -function ImperativeAffect(f::Function; +function ImperativeAffect(f; observed::NamedTuple = NamedTuple{()}(()), modified::NamedTuple = NamedTuple{()}(()), ctx = nothing, @@ -48,18 +48,18 @@ function ImperativeAffect(f::Function; collect(values(modified)), collect(keys(modified)), ctx, skip_checks) end -function ImperativeAffect(f::Function, modified::NamedTuple; +function ImperativeAffect(f, modified::NamedTuple; observed::NamedTuple = NamedTuple{()}(()), ctx = nothing, skip_checks = false) ImperativeAffect( f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) end function ImperativeAffect( - f::Function, modified::NamedTuple, observed::NamedTuple; ctx = nothing, skip_checks = false) + f, modified::NamedTuple, observed::NamedTuple; ctx = nothing, skip_checks = false) ImperativeAffect( f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) end function ImperativeAffect( - f::Function, modified::NamedTuple, observed::NamedTuple, ctx; skip_checks = false) + f, modified::NamedTuple, observed::NamedTuple, ctx; skip_checks = false) ImperativeAffect( f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) end @@ -216,3 +216,20 @@ function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs. end scalarize_affects(affects::ImperativeAffect) = affects + +function vars!(vars, aff::ImperativeAffect; op = Differential) + for var in Iterators.flatten((observed(aff), modified(aff))) + if symbolic_type(var) == NotSymbolic() + if var isa AbstractArray + for v in var + v = unwrap(v) + vars!(vars, v) + end + end + else + var = unwrap(var) + vars!(vars, var) + end + end + return vars +end diff --git a/src/systems/index_cache.jl b/src/systems/index_cache.jl index 623623da42..e4ca0b6b48 100644 --- a/src/systems/index_cache.jl +++ b/src/systems/index_cache.jl @@ -299,6 +299,9 @@ function IndexCache(sys::AbstractSystem) for v in vs if (idx = get(disc_idxs, v, nothing)) !== nothing push!(timeseries, idx.clock_idx) + elseif iscall(v) && operation(v) === getindex && + (idx = get(disc_idxs, arguments(v)[1], nothing)) !== nothing + push!(timeseries, idx.clock_idx) elseif haskey(observed_syms_to_timeseries, v) union!(timeseries, observed_syms_to_timeseries[v]) elseif haskey(dependent_pars_to_timeseries, v) diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl index 602fc7cac9..f71be80189 100644 --- a/src/systems/nonlinear/initializesystem.jl +++ b/src/systems/nonlinear/initializesystem.jl @@ -186,6 +186,14 @@ function generate_initializesystem(sys::AbstractSystem; defs[eq.lhs] = eq.rhs end + # even if `p => tovar(p)` is in `paramsubs`, `isparameter(p[1]) === true` after substitution + # so add scalarized versions as well + for k in collect(keys(paramsubs)) + symbolic_type(k) == ArraySymbolic() || continue + for i in eachindex(k) + paramsubs[k[i]] = paramsubs[k][i] + end + end eqs_ics = Symbolics.substitute.([eqs_ics; trueobs], (paramsubs,)) vars = [vars; collect(values(paramsubs))] for k in keys(defs) diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl index 6ce5704daa..55cbc32317 100644 --- a/src/systems/problem_utils.jl +++ b/src/systems/problem_utils.jl @@ -88,6 +88,18 @@ function symbols_to_symbolics!(sys::AbstractSystem, varmap::AbstractDict) end end +""" + $(TYPEDSIGNATURES) + +Utility function to get the value `val` corresponding to key `var` in `varmap`, and +return `getindex(val, idx)` if it exists or `nothing` otherwise. +""" +function get_and_getindex(varmap, var, idx) + val = get(varmap, var, nothing) + val === nothing && return nothing + return val[idx] +end + """ $(TYPEDSIGNATURES) @@ -115,8 +127,9 @@ function add_fallbacks!( val = map(eachindex(var)) do idx # @something is lazy and saves from writing a massive if-elseif-else @something(get(varmap, var[idx], nothing), - get(varmap, ttvar[idx], nothing), get(fallbacks, var, nothing)[idx], - get(fallbacks, ttvar, nothing)[idx], get(fallbacks, var[idx], nothing), + get(varmap, ttvar[idx], nothing), get_and_getindex(fallbacks, var, idx), + get_and_getindex(fallbacks, ttvar, idx), get( + fallbacks, var[idx], nothing), get(fallbacks, ttvar[idx], nothing), Some(nothing)) end # only push the missing entries @@ -561,6 +574,10 @@ function maybe_build_initialization_problem( p = unwrap(p) stype = symtype(p) op[p] = get_temporary_value(p) + if iscall(p) && operation(p) === getindex + arrp = arguments(p)[1] + op[arrp] = collect(arrp) + end end for v in missing_unknowns diff --git a/src/systems/systems.jl b/src/systems/systems.jl index 04c50bc766..9c8c272c5c 100644 --- a/src/systems/systems.jl +++ b/src/systems/systems.jl @@ -57,6 +57,7 @@ function structural_simplify( ks = collect(keys(defs)) # take copy to avoid mutating defs while iterating. for k in ks if Symbolics.isarraysymbolic(k) && Symbolics.shape(k) !== Symbolics.Unknown() + defs[k] === missing && continue for i in eachindex(k) defs[k[i]] = defs[k][i] end diff --git a/src/utils.jl b/src/utils.jl index c3011c2a79..3d9dae3ba9 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -403,6 +403,9 @@ function vars!(vars, eq::Equation; op = Differential) end function vars!(vars, O; op = Differential) if isvariable(O) + if iscall(O) && operation(O) === getindex && iscalledparameter(first(arguments(O))) + O = first(arguments(O)) + end if iscalledparameter(O) f = getcalledparameter(O) push!(vars, f) @@ -1114,6 +1117,23 @@ function subexpressions_not_involving_vars!(expr, vars, state::Dict{Any, Any}) symtype(expr) <: Union{Real, AbstractArray{<:Real}} || return expr Symbolics.shape(expr) == Symbolics.Unknown() && return expr haskey(state, expr) && return state[expr] + + op = operation(expr) + args = arguments(expr) + if op == StructuralTransformations.change_origin + newargs = map(args) do arg + if symbolic_type(arg) == NotSymbolic() + if is_array_of_symbolics(arg) + return map(arg) do _arg + subexpressions_not_involving_vars!(_arg, vars, state) + end + end + return arg + end + subexpressions_not_involving_vars!(arg, vars, state) + end + return maketerm(typeof(expr), op, newargs, metadata(expr)) + end vs = ModelingToolkit.vars(expr) intersect!(vs, vars) if isempty(vs) @@ -1123,8 +1143,6 @@ function subexpressions_not_involving_vars!(expr, vars, state::Dict{Any, Any}) state[expr] = var return var end - op = operation(expr) - args = arguments(expr) if (op == (+) || op == (*)) && symbolic_type(expr) !== ArraySymbolic() indep_args = [] dep_args = [] @@ -1143,7 +1161,14 @@ function subexpressions_not_involving_vars!(expr, vars, state::Dict{Any, Any}) return op(indep_term, dep_term) end newargs = map(args) do arg - symbolic_type(arg) != NotSymbolic() || is_array_of_symbolics(arg) || return arg + if symbolic_type(arg) == NotSymbolic() + if is_array_of_symbolics(arg) + return map(arg) do _arg + subexpressions_not_involving_vars!(_arg, vars, state) + end + end + return arg + end subexpressions_not_involving_vars!(arg, vars, state) end return maketerm(typeof(expr), op, newargs, metadata(expr)) diff --git a/test/extensions/Project.toml b/test/extensions/Project.toml index 5f7afe222a..72501554b5 100644 --- a/test/extensions/Project.toml +++ b/test/extensions/Project.toml @@ -2,6 +2,8 @@ BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a" +FMI = "14a09403-18e3-468f-ad8a-74f8dda2d9ac" +FMIZoo = "724179cf-c260-40a9-bd27-cccc6fe2f195" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327" InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57" diff --git a/test/extensions/fmi.jl b/test/extensions/fmi.jl new file mode 100644 index 0000000000..145989ab29 --- /dev/null +++ b/test/extensions/fmi.jl @@ -0,0 +1,309 @@ +using ModelingToolkit, FMI, FMIZoo, OrdinaryDiffEq +using ModelingToolkit: t_nounits as t, D_nounits as D +import ModelingToolkit as MTK + +const FMU_DIR = joinpath(@__DIR__, "fmus") + +@testset "Standalone pendulum model" begin + fmu = loadFMU("SpringPendulum1D", "Dymola", "2022x"; type = :ME) + truesol = FMI.simulate( + fmu, (0.0, 8.0); saveat = 0.0:0.1:8.0, recordValues = ["mass.s", "mass.v"]) + + function test_no_inputs_outputs(sys) + for var in unknowns(sys) + @test !MTK.isinput(var) + @test !MTK.isoutput(var) + end + end + @testset "v2, ME" begin + fmu = loadFMU("SpringPendulum1D", "Dymola", "2022x"; type = :ME) + @mtkbuild sys = MTK.FMIComponent(Val(2); fmu, type = :ME) + test_no_inputs_outputs(sys) + prob = ODEProblem{true, SciMLBase.FullSpecialize}( + sys, [sys.mass__s => 0.5, sys.mass__v => 0.0], (0.0, 8.0)) + sol = solve(prob, Tsit5(); reltol = 1e-8, abstol = 1e-8) + @test SciMLBase.successful_retcode(sol) + + @test sol(0.0:0.1:8.0; + idxs = [sys.mass__s, sys.mass__v]).u≈collect.(truesol.values.saveval) atol=1e-4 + # repeated solve works + @test_nowarn solve(prob, Tsit5()) + end + @testset "v2, CS" begin + fmu = loadFMU("SpringPendulum1D", "Dymola", "2022x"; type = :CS) + @named inner = MTK.FMIComponent( + Val(2); fmu, communication_step_size = 0.001, type = :CS) + @variables x(t) = 1.0 + @mtkbuild sys = ODESystem([D(x) ~ x], t; systems = [inner]) + test_no_inputs_outputs(sys) + + prob = ODEProblem{true, SciMLBase.FullSpecialize}( + sys, [sys.inner.mass__s => 0.5, sys.inner.mass__v => 0.0], (0.0, 8.0)) + sol = solve(prob, Tsit5(); reltol = 1e-8, abstol = 1e-8) + @test SciMLBase.successful_retcode(sol) + + @test sol(0.0:0.1:8.0; + idxs = [sys.inner.mass__s, sys.inner.mass__v]).u≈collect.(truesol.values.saveval) rtol=1e-2 + end + + fmu = loadFMU("SpringPendulum1D", "Dymola", "2023x", "3.0"; type = :ME) + truesol = FMI.simulate( + fmu, (0.0, 8.0); saveat = 0.0:0.1:8.0, recordValues = ["mass.s", "mass.v"]) + @testset "v3, ME" begin + fmu = loadFMU("SpringPendulum1D", "Dymola", "2023x", "3.0"; type = :ME) + @mtkbuild sys = MTK.FMIComponent(Val(3); fmu, type = :ME) + test_no_inputs_outputs(sys) + prob = ODEProblem{true, SciMLBase.FullSpecialize}( + sys, [sys.mass__s => 0.5, sys.mass__v => 0.0], (0.0, 8.0)) + sol = solve(prob, Tsit5(); reltol = 1e-8, abstol = 1e-8) + @test SciMLBase.successful_retcode(sol) + + @test sol(0.0:0.1:8.0; + idxs = [sys.mass__s, sys.mass__v]).u≈collect.(truesol.values.saveval) atol=1e-4 + # repeated solve works + @test_nowarn solve(prob, Tsit5()) + end + @testset "v3, CS" begin + fmu = loadFMU("SpringPendulum1D", "Dymola", "2023x", "3.0"; type = :CS) + @named inner = MTK.FMIComponent( + Val(3); fmu, communication_step_size = 0.001, type = :CS) + @variables x(t) = 1.0 + @mtkbuild sys = ODESystem([D(x) ~ x], t; systems = [inner]) + test_no_inputs_outputs(sys) + + prob = ODEProblem{true, SciMLBase.FullSpecialize}( + sys, [sys.inner.mass__s => 0.5, sys.inner.mass__v => 0.0], (0.0, 8.0)) + sol = solve(prob, Tsit5(); reltol = 1e-8, abstol = 1e-8) + @test SciMLBase.successful_retcode(sol) + + @test sol(0.0:0.1:8.0; + idxs = [sys.inner.mass__s, sys.inner.mass__v]).u≈collect.(truesol.values.saveval) rtol=1e-2 + end +end + +@mtkmodel SimpleAdder begin + @variables begin + a(t) + b(t) + c(t) + out(t) + out2(t) + end + @parameters begin + value = 1.0 + end + @equations begin + out ~ a + b + value + D(c) ~ out + out2 ~ 2c + end +end + +@mtkmodel StateSpace begin + @variables begin + x(t) + y(t) + u(t) + end + @parameters begin + A = 1.0 + B = 1.0 + C = 1.0 + _D = 1.0 + end + @equations begin + D(x) ~ A * x + B * u + y ~ C * x + _D * u + end +end + +@testset "IO Model" begin + function build_simple_adder(adder) + @variables a(t) b(t) c(t) [guess = 1.0] + @mtkbuild sys = ODESystem( + [adder.a ~ a, adder.b ~ b, D(a) ~ t, + D(b) ~ adder.out + adder.c, c^2 ~ adder.out + adder.value], + t; + systems = [adder]) + # c will be solved for by initialization + # this tests that initialization also works with FMUs + prob = ODEProblem(sys, [sys.adder.c => 2.0, sys.a => 1.0, sys.b => 1.0], + (0.0, 1.0), [sys.adder.value => 2.0]) + return sys, prob + end + + @named adder = SimpleAdder() + truesys, trueprob = build_simple_adder(adder) + truesol = solve(trueprob, abstol = 1e-8, reltol = 1e-8) + @test SciMLBase.successful_retcode(truesol) + + @testset "v2, ME" begin + fmu = loadFMU(joinpath(FMU_DIR, "SimpleAdder.fmu"); type = :ME) + @named adder = MTK.FMIComponent(Val(2); fmu, type = :ME) + @test MTK.isinput(adder.a) + @test MTK.isinput(adder.b) + @test MTK.isoutput(adder.out) + @test MTK.isoutput(adder.out2) + @test !MTK.isinput(adder.c) && !MTK.isoutput(adder.c) + + sys, prob = build_simple_adder(adder) + sol = solve(prob, Rodas5P(autodiff = false), abstol = 1e-8, reltol = 1e-8) + @test SciMLBase.successful_retcode(sol) + + @test truesol(sol.t; + idxs = [truesys.a, truesys.b, truesys.c, truesys.adder.c]).u≈sol[[ + sys.a, sys.b, sys.c, sys.adder.c]] rtol=1e-7 + end + @testset "v2, CS" begin + fmu = loadFMU(joinpath(FMU_DIR, "SimpleAdder.fmu"); type = :CS) + @named adder = MTK.FMIComponent( + Val(2); fmu, type = :CS, communication_step_size = 1e-3, + reinitializealg = BrownFullBasicInit()) + @test MTK.isinput(adder.a) + @test MTK.isinput(adder.b) + @test MTK.isoutput(adder.out) + @test MTK.isoutput(adder.out2) + @test !MTK.isinput(adder.c) && !MTK.isoutput(adder.c) + + sys, prob = build_simple_adder(adder) + sol = solve(prob, Rodas5P(autodiff = false), abstol = 1e-8, reltol = 1e-8) + @test SciMLBase.successful_retcode(sol) + + @test truesol(sol.t; idxs = [truesys.a, truesys.b, truesys.c]).u≈sol.u rtol=1e-2 + # sys.adder.c is a discrete variable + @test truesol(sol.t; idxs = truesys.adder.c)≈sol(sol.t; idxs = sys.adder.c) rtol=1e-3 + end + + function build_sspace_model(sspace) + @variables u(t)=1.0 x(t)=1.0 y(t) [guess = 1.0] + @mtkbuild sys = ODESystem( + [sspace.u ~ u, D(u) ~ t, D(x) ~ sspace.x + sspace.y, y^2 ~ sspace.y + sspace.x], t; + systems = [sspace] + ) + + prob = ODEProblem( + sys, [sys.sspace.x => 1.0], (0.0, 1.0), [sys.sspace.A => 2.0]; use_scc = false) + return sys, prob + end + + @named sspace = StateSpace() + truesys, trueprob = build_sspace_model(sspace) + truesol = solve(trueprob, abstol = 1e-8, reltol = 1e-8) + @test SciMLBase.successful_retcode(truesol) + + @testset "v3, ME" begin + fmu = loadFMU(joinpath(FMU_DIR, "StateSpace.fmu"); type = :ME) + @named sspace = MTK.FMIComponent(Val(3); fmu, type = :ME) + @test MTK.isinput(sspace.u) + @test MTK.isoutput(sspace.y) + @test !MTK.isinput(sspace.x) && !MTK.isoutput(sspace.x) + + sys, prob = build_sspace_model(sspace) + sol = solve(prob, Rodas5P(autodiff = false); abstol = 1e-8, reltol = 1e-8) + @test SciMLBase.successful_retcode(sol) + + @test truesol(sol.t; + idxs = [truesys.u, truesys.x, truesys.y, truesys.sspace.x]).u≈sol[[ + sys.u, sys.x, sys.y, sys.sspace.x]] rtol=1e-7 + end + + @testset "v3, CS" begin + fmu = loadFMU(joinpath(FMU_DIR, "StateSpace.fmu"); type = :CS) + @named sspace = MTK.FMIComponent( + Val(3); fmu, communication_step_size = 1e-4, type = :CS, + reinitializealg = BrownFullBasicInit()) + @test MTK.isinput(sspace.u) + @test MTK.isoutput(sspace.y) + @test !MTK.isinput(sspace.x) && !MTK.isoutput(sspace.x) + + sys, prob = build_sspace_model(sspace) + sol = solve(prob, Rodas5P(autodiff = false); abstol = 1e-8, reltol = 1e-8) + @test SciMLBase.successful_retcode(sol) + + @test truesol( + sol.t; idxs = [truesys.u, truesys.x, truesys.y]).u≈sol[[sys.u, sys.x, sys.y]] rtol=1e-2 + @test truesol(sol.t; idxs = truesys.sspace.x).u≈sol(sol.t; idxs = sys.sspace.x).u rtol=1e-2 + end +end + +@testset "FMUs in a loop" begin + function build_looped_adders(adder1, adder2) + @variables x(t) = 1 + @mtkbuild sys = ODESystem( + [D(x) ~ x, adder1.a ~ adder2.out2, + adder2.a ~ adder1.out2, adder1.b ~ 1.0, adder2.b ~ 2.0], + t; + systems = [adder1, adder2]) + prob = ODEProblem( + sys, [adder1.c => 1.0, adder2.c => 1.0, adder1.a => 2.0], (0.0, 1.0)) + return sys, prob + end + @named adder1 = SimpleAdder() + @named adder2 = SimpleAdder() + truesys, trueprob = build_looped_adders(adder1, adder2) + truesol = solve(trueprob, Tsit5(), reltol = 1e-8) + @test SciMLBase.successful_retcode(truesol) + + @testset "v2, ME" begin + fmu = loadFMU(joinpath(FMU_DIR, "SimpleAdder.fmu"); type = :ME) + @named adder1 = MTK.FMIComponent(Val(2); fmu, type = :ME) + @named adder2 = MTK.FMIComponent(Val(2); fmu, type = :ME) + sys, prob = build_looped_adders(adder1, adder2) + sol = solve(prob, Rodas5P(autodiff = false); reltol = 1e-8) + @test SciMLBase.successful_retcode(sol) + @test truesol(sol.t; + idxs = [truesys.adder1.c, truesys.adder2.c]).u≈sol( + sol.t; idxs = [sys.adder1.c, sys.adder2.c]).u rtol=1e-7 + end + @testset "v2, CS" begin + fmu = loadFMU(joinpath(FMU_DIR, "SimpleAdder.fmu"); type = :CS) + @named adder1 = MTK.FMIComponent( + Val(2); fmu, type = :CS, communication_step_size = 1e-3) + @named adder2 = MTK.FMIComponent( + Val(2); fmu, type = :CS, communication_step_size = 1e-3) + sys, prob = build_looped_adders(adder1, adder2) + sol = solve(prob, Tsit5(); reltol = 1e-8) + @test truesol(sol.t; + idxs = [truesys.adder1.c, truesys.adder2.c]).u≈sol( + sol.t; idxs = [sys.adder1.c, sys.adder2.c]).u rtol=1e-3 + end + + function build_looped_sspace(sspace1, sspace2) + @variables x(t) = 1 + @mtkbuild sys = ODESystem([D(x) ~ x, sspace1.u ~ sspace2.x, sspace2.u ~ sspace1.y], + t; systems = [sspace1, sspace2]) + prob = ODEProblem(sys, [sspace1.x => 1.0, sspace2.x => 1.0], (0.0, 1.0)) + return sys, prob + end + @named sspace1 = StateSpace() + @named sspace2 = StateSpace() + truesys, trueprob = build_looped_sspace(sspace1, sspace2) + truesol = solve(trueprob, Rodas5P(), reltol = 1e-8) + @test SciMLBase.successful_retcode(truesol) + + @testset "v3, ME" begin + fmu = loadFMU(joinpath(FMU_DIR, "StateSpace.fmu"); type = :ME) + @named sspace1 = MTK.FMIComponent(Val(3); fmu, type = :ME) + @named sspace2 = MTK.FMIComponent(Val(3); fmu, type = :ME) + sys, prob = build_looped_sspace(sspace1, sspace2) + sol = solve(prob, Rodas5P(autodiff = false); reltol = 1e-8) + @test SciMLBase.successful_retcode(sol) + @test truesol(sol.t; + idxs = [truesys.sspace1.x, truesys.sspace2.x]).u≈sol( + sol.t; idxs = [sys.sspace1.x, sys.sspace2.x]).u rtol=1e-7 + end + + @testset "v3, CS" begin + fmu = loadFMU(joinpath(FMU_DIR, "StateSpace.fmu"); type = :CS) + @named sspace1 = MTK.FMIComponent( + Val(3); fmu, type = :CS, communication_step_size = 1e-4) + @named sspace2 = MTK.FMIComponent( + Val(3); fmu, type = :CS, communication_step_size = 1e-4) + sys, prob = build_looped_sspace(sspace1, sspace2) + sol = solve(prob, Rodas5P(autodiff = false); reltol = 1e-8) + @test SciMLBase.successful_retcode(sol) + @test truesol(sol.t; + idxs = [truesys.sspace1.x, truesys.sspace2.x]).u≈sol( + sol.t; idxs = [sys.sspace1.x, sys.sspace2.x]).u rtol=1e-2 + end +end diff --git a/test/extensions/fmus/SimpleAdder.fmu b/test/extensions/fmus/SimpleAdder.fmu new file mode 100644 index 0000000000..ed12961bfb Binary files /dev/null and b/test/extensions/fmus/SimpleAdder.fmu differ diff --git a/test/extensions/fmus/SimpleAdder/README.md b/test/extensions/fmus/SimpleAdder/README.md new file mode 100644 index 0000000000..b26c247bce --- /dev/null +++ b/test/extensions/fmus/SimpleAdder/README.md @@ -0,0 +1,3 @@ +SimpleAdder.fmu is a v2 FMU built using OpenModelica. The file `SimpleAdder.mo` contains the +modelica source code for the model. The file `buildFMU.mos` is the OpenModelica script used +to build the FMU. diff --git a/test/extensions/fmus/SimpleAdder/SimpleAdder.mo b/test/extensions/fmus/SimpleAdder/SimpleAdder.mo new file mode 100644 index 0000000000..ce97f09a86 --- /dev/null +++ b/test/extensions/fmus/SimpleAdder/SimpleAdder.mo @@ -0,0 +1,12 @@ +block SimpleAdder + parameter Real value = 1.0; + input Real a; + input Real b; + Real c(start = 1.0, fixed = true); + output Real out; + output Real out2; +equation + der(c) = out; + out = a + b + value; + out2 = 2 * c; +end SimpleAdder; diff --git a/test/extensions/fmus/SimpleAdder/buildFmu.mos b/test/extensions/fmus/SimpleAdder/buildFmu.mos new file mode 100644 index 0000000000..df2c9bd58e --- /dev/null +++ b/test/extensions/fmus/SimpleAdder/buildFmu.mos @@ -0,0 +1,12 @@ +OpenModelica.Scripting.loadFile("SimpleAdder.mo"); getErrorString(); + +installPackage(Modelica); + +setCommandLineOptions("-d=newInst"); getErrorString(); +setCommandLineOptions("-d=initialization"); getErrorString(); +setCommandLineOptions("-d=-disableDirectionalDerivatives"); getErrorString(); + +cd("output"); getErrorString(); +buildModelFMU(SimpleAdder, version = "2.0", fmuType = "me_cs"); getErrorString(); +system("unzip -l SimpleAdder.fmu | egrep -v 'sources|files' | tail -n+3 | grep -o '[A-Za-z._0-9/]*$' > BB.log") + diff --git a/test/extensions/fmus/StateSpace.fmu b/test/extensions/fmus/StateSpace.fmu new file mode 100644 index 0000000000..6138e1fdb2 Binary files /dev/null and b/test/extensions/fmus/StateSpace.fmu differ diff --git a/test/extensions/fmus/StateSpace/0001-tmp-commit.patch b/test/extensions/fmus/StateSpace/0001-tmp-commit.patch new file mode 100644 index 0000000000..951af2d4fd --- /dev/null +++ b/test/extensions/fmus/StateSpace/0001-tmp-commit.patch @@ -0,0 +1,425 @@ +From 10ac73996a7c09772ff31ccd6e030112d3bb5c43 Mon Sep 17 00:00:00 2001 +From: Aayush Sabharwal +Date: Wed, 8 Jan 2025 11:04:01 +0000 +Subject: tmp commit + +--- + StateSpace/FMI3.xml | 31 ++---- + StateSpace/config.h | 18 ++-- + StateSpace/model.c | 230 +++++++++++++------------------------------- + 3 files changed, 83 insertions(+), 196 deletions(-) + +diff --git a/StateSpace/FMI3.xml b/StateSpace/FMI3.xml +index 8d2e4d9..002fbab 100644 +--- a/StateSpace/FMI3.xml ++++ b/StateSpace/FMI3.xml +@@ -33,36 +33,23 @@ + + + +- +- +- ++ + +- +- +- ++ + +- +- +- ++ + +- +- +- ++ + +- +- ++ + +- +- ++ + +- +- ++ + +- +- ++ + +- +- ++ + + + +diff --git a/StateSpace/config.h b/StateSpace/config.h +index 707e500..8b422e2 100644 +--- a/StateSpace/config.h ++++ b/StateSpace/config.h +@@ -44,15 +44,15 @@ typedef struct { + uint64_t m; + uint64_t n; + uint64_t r; +- double A[M_MAX][N_MAX]; +- double B[M_MAX][N_MAX]; +- double C[M_MAX][N_MAX]; +- double D[M_MAX][N_MAX]; +- double x0[N_MAX]; +- double u[N_MAX]; +- double y[N_MAX]; +- double x[N_MAX]; +- double der_x[N_MAX]; ++ double A; ++ double B; ++ double C; ++ double D; ++ double x0; ++ double u; ++ double y; ++ double x; ++ double der_x; + } ModelData; + + #endif /* config_h */ +diff --git a/StateSpace/model.c b/StateSpace/model.c +index 8a47e74..9c170d8 100644 +--- a/StateSpace/model.c ++++ b/StateSpace/model.c +@@ -3,77 +3,29 @@ + + + void setStartValues(ModelInstance *comp) { +- +- M(m) = 3; +- M(n) = 3; +- M(r) = 3; +- + // identity matrix +- for (int i = 0; i < M_MAX; i++) +- for (int j = 0; j < N_MAX; j++) { +- M(A)[i][j] = i == j ? 1 : 0; +- M(B)[i][j] = i == j ? 1 : 0; +- M(C)[i][j] = i == j ? 1 : 0; +- M(D)[i][j] = i == j ? 1 : 0; +- } +- +- for (int i = 0; i < M_MAX; i++) { +- M(u)[i] = i + 1; +- } +- +- for (int i = 0; i < N_MAX; i++) { +- M(y)[i] = 0; +- } +- +- for (int i = 0; i < N_MAX; i++) { +- M(x)[i] = M(x0)[i]; +- M(x)[i] = 0; +- } +- ++ M(m) = 1; ++ M(n) = 1; ++ M(r) = 1; ++ ++ M(A) = 1; ++ M(B) = 1; ++ M(C) = 1; ++ M(D) = 1; ++ ++ M(u) = 1; ++ M(y) = 0; ++ M(x) = 0; ++ M(x0) = 0; + } + + Status calculateValues(ModelInstance *comp) { +- +- // der(x) = Ax + Bu +- for (size_t i = 0; i < M(n); i++) { +- +- M(der_x)[i] = 0; +- +- for (size_t j = 0; j < M(n); j++) { +- M(der_x)[i] += M(A)[i][j] * M(x)[j]; +- } +- } +- +- for (size_t i = 0; i < M(n); i++) { +- +- for (size_t j = 0; j < M(r); j++) { +- M(der_x)[i] += M(B)[i][j] * M(u)[j]; +- } +- } +- +- +- // y = Cx + Du +- for (size_t i = 0; i < M(r); i++) { +- +- M(y)[i] = 0; +- +- for (size_t j = 0; j < M(n); j++) { +- M(y)[i] += M(C)[i][j] * M(x)[j]; +- } +- } +- +- for (size_t i = 0; i < M(r); i++) { +- +- for (size_t j = 0; j < M(m); j++) { +- M(y)[i] += M(D)[i][j] * M(u)[j]; +- } +- } +- ++ M(der_x) = M(A) * M(x) + M(B) * M(u); ++ M(y) = M(C) * M(x) + M(D) * M(u); + return OK; + } + + Status getFloat64(ModelInstance* comp, ValueReference vr, double values[], size_t nValues, size_t* index) { +- + calculateValues(comp); + + switch (vr) { +@@ -82,66 +34,40 @@ Status getFloat64(ModelInstance* comp, ValueReference vr, double values[], size_ + values[(*index)++] = comp->time; + return OK; + case vr_A: +- ASSERT_NVALUES((size_t)(M(n) * M(n))); +- for (size_t i = 0; i < M(n); i++) { +- for (size_t j = 0; j < M(n); j++) { +- values[(*index)++] = M(A)[i][j]; +- } +- } ++ ASSERT_NVALUES(1); ++ values[(*index)++] = M(A); + return OK; + case vr_B: +- ASSERT_NVALUES((size_t)(M(m) * M(n))); +- for (size_t i = 0; i < M(m); i++) { +- for (size_t j = 0; j < M(n); j++) { +- values[(*index)++] = M(B)[i][j]; +- } +- } ++ ASSERT_NVALUES(1); ++ values[(*index)++] = M(B); + return OK; + case vr_C: +- ASSERT_NVALUES((size_t)(M(r) * M(n))); +- for (size_t i = 0; i < M(r); i++) { +- for (size_t j = 0; j < M(n); j++) { +- values[(*index)++] = M(C)[i][j]; +- } +- } ++ ASSERT_NVALUES(1); ++ values[(*index)++] = M(C); + return OK; + case vr_D: +- ASSERT_NVALUES((size_t)(M(r) * M(m))); +- for (size_t i = 0; i < M(r); i++) { +- for (size_t j = 0; j < M(m); j++) { +- values[(*index)++] = M(D)[i][j]; +- } +- } ++ ASSERT_NVALUES(1); ++ values[(*index)++] = M(D); + return OK; + case vr_x0: +- ASSERT_NVALUES((size_t)M(n)); +- for (size_t i = 0; i < M(n); i++) { +- values[(*index)++] = M(x0)[i]; +- } ++ ASSERT_NVALUES(1); ++ values[(*index)++] = M(x0); + return OK; + case vr_u: +- ASSERT_NVALUES((size_t)M(m)); +- for (size_t i = 0; i < M(m); i++) { +- values[(*index)++] = M(u)[i]; +- } ++ ASSERT_NVALUES(1); ++ values[(*index)++] = M(u); + return OK; + case vr_y: +- ASSERT_NVALUES((size_t)M(r)); +- for (size_t i = 0; i < M(r); i++) { +- values[(*index)++] = M(y)[i]; +- } ++ ASSERT_NVALUES(1); ++ values[(*index)++] = M(y); + return OK; + case vr_x: +- ASSERT_NVALUES((size_t)M(n)); +- for (size_t i = 0; i < M(n); i++) { +- values[(*index)++] = M(x)[i]; +- } ++ ASSERT_NVALUES(1); ++ values[(*index)++] = M(x); + return OK; + case vr_der_x: +- ASSERT_NVALUES((size_t)M(n)); +- for (size_t i = 0; i < M(n); i++) { +- values[(*index)++] = M(der_x)[i]; +- } ++ ASSERT_NVALUES(1); ++ values[(*index)++] = M(der_x); + return OK; + default: + logError(comp, "Get Float64 is not allowed for value reference %u.", vr); +@@ -153,58 +79,40 @@ Status setFloat64(ModelInstance* comp, ValueReference vr, const double values[], + + switch (vr) { + case vr_A: +- ASSERT_NVALUES((size_t)(M(n) * M(n))); +- for (size_t i = 0; i < M(n); i++) { +- for (size_t j = 0; j < M(n); j++) { +- M(A)[i][j] = values[(*index)++]; +- } +- } ++ ASSERT_NVALUES(1); ++ M(A) = values[(*index)++]; + break; + case vr_B: +- ASSERT_NVALUES((size_t)(M(n) * M(m))); +- for (size_t i = 0; i < M(n); i++) { +- for (size_t j = 0; j < M(m); j++) { +- M(B)[i][j] = values[(*index)++]; +- } +- } ++ ASSERT_NVALUES(1); ++ M(B) = values[(*index)++]; + break; + case vr_C: +- ASSERT_NVALUES((size_t)(M(r) * M(n))); +- for (size_t i = 0; i < M(r); i++) { +- for (size_t j = 0; j < M(n); j++) { +- M(C)[i][j] = values[(*index)++]; +- } +- } ++ ASSERT_NVALUES(1); ++ M(C) = values[(*index)++]; + break; + case vr_D: +- ASSERT_NVALUES((size_t)(M(r) * M(m))); +- for (size_t i = 0; i < M(r); i++) { +- for (size_t j = 0; j < M(m); j++) { +- M(D)[i][j] = values[(*index)++]; +- } +- } ++ ASSERT_NVALUES(1); ++ M(D) = values[(*index)++]; + break; + case vr_x0: +- ASSERT_NVALUES((size_t)M(n)); +- for (size_t i = 0; i < M(n); i++) { +- M(x0)[i] = values[(*index)++]; +- } ++ ASSERT_NVALUES(1); ++ M(x0) = values[(*index)++]; + break; + case vr_u: +- ASSERT_NVALUES((size_t)M(m)); +- for (size_t i = 0; i < M(m); i++) { +- M(u)[i] = values[(*index)++]; +- } ++ ASSERT_NVALUES(1); ++ M(u) = values[(*index)++]; ++ break; ++ case vr_y: ++ ASSERT_NVALUES(1); ++ M(y) = values[(*index)++]; + break; + case vr_x: +- if (comp->state != ContinuousTimeMode && comp->state != EventMode) { +- logError(comp, "Variable \"x\" can only be set in Continuous Time Mode and Event Mode."); +- return Error; +- } +- ASSERT_NVALUES((size_t)M(n)); +- for (size_t i = 0; i < M(n); i++) { +- M(x)[i] = values[(*index)++]; +- } ++ ASSERT_NVALUES(1); ++ M(x) = values[(*index)++]; ++ break; ++ case vr_der_x: ++ ASSERT_NVALUES(1); ++ M(der_x) = values[(*index)++]; + break; + default: + logError(comp, "Set Float64 is not allowed for value reference %u.", vr); +@@ -217,8 +125,6 @@ Status setFloat64(ModelInstance* comp, ValueReference vr, const double values[], + } + + Status getUInt64(ModelInstance* comp, ValueReference vr, uint64_t values[], size_t nValues, size_t* index) { +- +- + switch (vr) { + case vr_m: + ASSERT_NVALUES(1); +@@ -289,49 +195,43 @@ Status eventUpdate(ModelInstance *comp) { + } + + size_t getNumberOfContinuousStates(ModelInstance* comp) { +- return (size_t)M(n); ++ return M(n) == 1 ? M(n) : 1; + } + + Status getContinuousStates(ModelInstance* comp, double x[], size_t nx) { + +- if (nx != M(n)) { +- logError(comp, "Expected nx=%zu but was %zu.", M(n), nx); ++ if (nx != 1) { ++ logError(comp, "Expected nx=%zu but was %zu.", 1, nx); + return Error; + } + +- for (size_t i = 0; i < M(n); i++) { +- x[i] = M(x)[i]; +- } ++ x[0] = M(x); + + return OK; + } + + Status setContinuousStates(ModelInstance* comp, const double x[], size_t nx) { + +- if (nx != M(n)) { +- logError(comp, "Expected nx=%zu but was %zu.", M(n), nx); ++ if (nx != 1) { ++ logError(comp, "Expected nx=%zu but was %zu.", 1, nx); + return Error; + } + +- for (size_t i = 0; i < M(n); i++) { +- M(x)[i] = x[i]; +- } ++ M(x) = x[0]; + + return OK; + } + + Status getDerivatives(ModelInstance* comp, double dx[], size_t nx) { + +- if (nx != M(n)) { +- logError(comp, "Expected nx=%zu but was %zu.", M(n), nx); ++ if (nx != 1) { ++ logError(comp, "Expected nx=%zu but was %zu.", 1, nx); + return Error; + } + + calculateValues(comp); + +- for (size_t i = 0; i < M(n); i++) { +- dx[i] = M(der_x)[i]; +- } ++ dx[0] = M(der_x); + + return OK; + } +-- +2.34.1 + diff --git a/test/extensions/fmus/StateSpace/README.md b/test/extensions/fmus/StateSpace/README.md new file mode 100644 index 0000000000..108147443d --- /dev/null +++ b/test/extensions/fmus/StateSpace/README.md @@ -0,0 +1,7 @@ +StateSpace.fmu is a v3 FMU built using a modified version of the StateSpace FMU in Modelica's +Reference-FMUs. The link to the specific commit hash used is + +[https://github.com/modelica/Reference-FMUs/tree/0e1374cad0a7f0583a583ce189060ba7a67c27a7]() + +The git patch containing the changes is provided as `0001-tmp-commit.patch`. The FMU can be +built using the instructions provided in the repository's README. diff --git a/test/runtests.jl b/test/runtests.jl index 4a2ec80e6a..1cf53402fd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -114,6 +114,7 @@ end if GROUP == "All" || GROUP == "Extensions" activate_extensions_env() + @safetestset "FMI Extension Test" include("extensions/fmi.jl") @safetestset "HomotopyContinuation Extension Test" include("extensions/homotopy_continuation.jl") @safetestset "Auto Differentiation Test" include("extensions/ad.jl") @safetestset "LabelledArrays Test" include("labelledarrays.jl")