SciML
Internal Error: Variables Incorrectly Marked as Zero During Structural Simplification in ODE System
Describe the bug 🐞 When running a differential equation system using ModelingToolkit, an internal error occurs where variables (N(t))[1] and (N(t))[2] are unexpectedly marked as zero in the equation 0 ~ -H(t) + myfun(T(t), N(t)). This issue arises during the structural transformation call.
Expected behavior The system should solve the ODE without internal errors, accurately computing the values of the variables (N(t))[1] and (N(t))[2] without them being incorrectly marked as zero.
Minimal Reproducible Example 👇
using ModelingToolkit
using Symbolics
using ModelingToolkit: t_nounits as t, D_nounits as D, scalarize
using OrdinaryDiffEq
myfun(T, x::Vector) = 100.0*T*x[1] + 50.0*T*x[2]
@register_symbolic myfun(T, x::Vector)
function myc(; name)
vars = @variables begin
(N(t))[1:2]
H(t)
T(t)
end
eqs = [scalarize(D(N) .~ -N), D(H) ~ H/sum(scalarize(N)),
H - myfun(T, N) ~ 0.0]
ODESystem([eqs...;], t, collect(Iterators.flatten(vars)), []; name)
end
component = myc(; name = :myc)
simple_prob = structural_simplify(component; check_consistency = true)
Error & Stacktrace ⚠️ It is a warning
┌ Warning: Internal error: Variable (N(t))[1] was marked as being in 0 ~ -H(t) + myfun(T(t), N(t)), but was actually zero
└ @ ModelingToolkit.StructuralTransformations C:\Users\vinicius\.julia\packages\ModelingToolkit\Jc09X\src\structural_transformation\utils.jl:233
┌ Warning: Internal error: Variable (N(t))[2] was marked as being in 0 ~ -H(t) + myfun(T(t), N(t)), but was actually zero
└ @ ModelingToolkit.StructuralTransformations C:\Users\vinicius\.julia\packages\ModelingToolkit\Jc09X\src\structural_transformation\utils.jl:233
Environment (please complete the following information):
- Output of
using Pkg; Pkg.status()
Status `C:\Users\vinicius\OneDrive - NTNU\Research\issuemtk\Project.toml`
[961ee093] ModelingToolkit v9.35.1
[1dea7af3] OrdinaryDiffEq v6.89.0
[0c5d862f] Symbolics v6.5.0
- Output of
using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)
Status `C:\Users\vinicius\OneDrive - NTNU\Research\issuemtk\Manifest.toml`
[47edcb42] ADTypes v1.7.1
[1520ce14] AbstractTrees v0.4.5
[7d9f7c33] Accessors v0.1.37
[79e6a3ab] Adapt v4.0.4
[66dad0bd] AliasTables v1.1.3
[ec485272] ArnoldiMethod v0.4.0
[4fba245c] ArrayInterface v7.16.0
[4c555306] ArrayLayouts v1.10.3
[e2ed5e7c] Bijections v0.1.9
[62783981] BitTwiddlingConvenienceFunctions v0.1.6
[8e7c35d0] BlockArrays v1.1.0
[2a0fbf3d] CPUSummary v0.2.6
[00ebfdb7] CSTParser v3.4.3
[d360d2e6] ChainRulesCore v1.24.0
[fb6a15b2] CloseOpenIntervals v0.1.13
[861a8166] Combinatorics v1.0.2
[a80b9123] CommonMark v0.8.12
[38540f10] CommonSolve v0.2.4
[bbf7d656] CommonSubexpressions v0.3.1
[f70d9fcc] CommonWorldInvalidations v1.0.0
[34da2185] Compat v4.16.0
[b152e2b5] CompositeTypes v0.1.4
[a33af91c] CompositionsBase v0.1.2
[2569d6c7] ConcreteStructs v0.2.3
[187b0558] ConstructionBase v1.5.7
[adafc99b] CpuId v0.3.1
[a8cc5b0e] Crayons v4.1.1
[9a962f9c] DataAPI v1.16.0
[864edb3b] DataStructures v0.18.20
[e2d170a0] DataValueInterfaces v1.0.0
[2b5f629d] DiffEqBase v6.154.0
[459566f4] DiffEqCallbacks v3.9.0
[77a26b50] DiffEqNoiseProcess v5.23.0
[163ba53b] DiffResults v1.1.0
[b552c78f] DiffRules v1.15.1
[a0c0ee7d] DifferentiationInterface v0.5.15
[31c24e10] Distributions v0.25.111
[ffbed154] DocStringExtensions v0.9.3
[5b8099bc] DomainSets v0.7.14
[7c1d4256] DynamicPolynomials v0.6.0
⌅ [06fc5a27] DynamicQuantities v0.13.2
[4e289a0a] EnumX v1.0.4
[f151be2c] EnzymeCore v0.7.8
[d4d017d3] ExponentialUtilities v1.26.1
[e2ba6199] ExprTools v0.1.10
⌅ [6b7a57c9] Expronicon v0.8.5
[7034ab61] FastBroadcast v0.3.5
[9aa1b823] FastClosures v0.3.2
[29a986be] FastLapackInterface v2.0.4
[1a297f60] FillArrays v1.13.0
[64ca27bc] FindFirstFunctions v1.3.0
[6a86dc24] FiniteDiff v2.24.0
[1fa38f19] Format v1.3.7
[f6369f11] ForwardDiff v0.10.36
[069b7b12] FunctionWrappers v1.1.3
[77dc65aa] FunctionWrappersWrappers v0.1.3
[d9f16b24] Functors v0.4.12
[46192b85] GPUArraysCore v0.1.6
[c145ed77] GenericSchur v0.5.4
[c27321d9] Glob v1.3.1
[86223c79] Graphs v1.11.2
[3e5b6fbb] HostCPUFeatures v0.1.17
[34004b35] HypergeometricFunctions v0.3.24
[615f187c] IfElse v0.1.1
[d25df0c9] Inflate v0.1.5
[18e54dd8] IntegerMathUtils v0.1.2
[8197267c] IntervalSets v0.7.10
[3587e190] InverseFunctions v0.1.16
[92d709cd] IrrationalConstants v0.2.2
[82899510] IteratorInterfaceExtensions v1.0.0
[692b3bcd] JLLWrappers v1.6.0
[682c06a0] JSON v0.21.4
[98e50ef6] JuliaFormatter v1.0.60
[ccbc3e58] JumpProcesses v9.13.5
[ef3ab10e] KLU v0.6.0
[ba0b0d4f] Krylov v0.9.6
[b964fa9f] LaTeXStrings v1.3.1
[984bce1d] LambertW v0.4.6
[23fbe1c1] Latexify v0.16.5
[10f19ff3] LayoutPointers v0.1.17
[5078a376] LazyArrays v2.2.0
[d3d80556] LineSearches v7.3.0
[7ed4a6bd] LinearSolve v2.34.0
[2ab3a3ac] LogExpFunctions v0.3.28
[bdcacae8] LoopVectorization v0.12.171
[d8e11817] MLStyle v0.4.17
[1914dd2f] MacroTools v0.5.13
[d125e4d3] ManualMemory v0.1.8
[bb5d69b7] MaybeInplace v0.1.4
[e1d29d7a] Missings v1.2.0
[961ee093] ModelingToolkit v9.35.1
[46d2c3a1] MuladdMacro v0.2.4
[102ac46a] MultivariatePolynomials v0.5.6
[d8a4904e] MutableArithmetics v1.4.6
[d41bc354] NLSolversBase v7.8.3
[77ba4419] NaNMath v1.0.2
[8913a72c] NonlinearSolve v3.14.0
[6fe1bfb0] OffsetArrays v1.14.1
[429524aa] Optim v1.9.4
[bac558e1] OrderedCollections v1.6.3
[1dea7af3] OrdinaryDiffEq v6.89.0
[89bda076] OrdinaryDiffEqAdamsBashforthMoulton v1.1.0
[6ad6398a] OrdinaryDiffEqBDF v1.1.1
[bbf590c4] OrdinaryDiffEqCore v1.4.0
[50262376] OrdinaryDiffEqDefault v1.1.0
[4302a76b] OrdinaryDiffEqDifferentiation v1.1.0
[9286f039] OrdinaryDiffEqExplicitRK v1.1.0
[e0540318] OrdinaryDiffEqExponentialRK v1.1.0
[becaefa8] OrdinaryDiffEqExtrapolation v1.1.0
[5960d6e9] OrdinaryDiffEqFIRK v1.1.0
[101fe9f7] OrdinaryDiffEqFeagin v1.1.0
[d3585ca7] OrdinaryDiffEqFunctionMap v1.1.1
[d28bc4f8] OrdinaryDiffEqHighOrderRK v1.1.0
[9f002381] OrdinaryDiffEqIMEXMultistep v1.1.0
[521117fe] OrdinaryDiffEqLinear v1.1.0
[1344f307] OrdinaryDiffEqLowOrderRK v1.1.0
[b0944070] OrdinaryDiffEqLowStorageRK v1.2.1
[127b3ac7] OrdinaryDiffEqNonlinearSolve v1.1.0
[c9986a66] OrdinaryDiffEqNordsieck v1.1.0
[5dd0a6cf] OrdinaryDiffEqPDIRK v1.1.0
[5b33eab2] OrdinaryDiffEqPRK v1.1.0
[04162be5] OrdinaryDiffEqQPRK v1.1.0
[af6ede74] OrdinaryDiffEqRKN v1.1.0
[43230ef6] OrdinaryDiffEqRosenbrock v1.1.1
[2d112036] OrdinaryDiffEqSDIRK v1.1.0
[669c94d9] OrdinaryDiffEqSSPRK v1.2.0
[e3e12d00] OrdinaryDiffEqStabilizedIRK v1.1.0
[358294b1] OrdinaryDiffEqStabilizedRK v1.1.0
[fa646aed] OrdinaryDiffEqSymplecticRK v1.1.0
[b1df2697] OrdinaryDiffEqTsit5 v1.1.0
[79d7bb75] OrdinaryDiffEqVerner v1.1.1
[90014a1f] PDMats v0.11.31
[65ce6f38] PackageExtensionCompat v1.0.2
[d96e819e] Parameters v0.12.3
[69de0a69] Parsers v2.8.1
[e409e4f3] PoissonRandom v0.4.4
[f517fe37] Polyester v0.7.16
[1d0040c9] PolyesterWeave v0.2.2
[85a6dd25] PositiveFactorizations v0.2.4
[d236fae5] PreallocationTools v0.4.24
[aea7be01] PrecompileTools v1.2.1
[21216c6a] Preferences v1.4.3
[27ebfcd6] Primes v0.5.6
[43287f4e] PtrArrays v1.2.1
[1fd47b50] QuadGK v2.11.0
[74087812] Random123 v1.7.0
[e6cf234a] RandomNumbers v1.6.0
[3cdcf5f2] RecipesBase v1.3.4
[731186ca] RecursiveArrayTools v3.27.0
[f2c3362d] RecursiveFactorization v0.2.23
[189a3867] Reexport v1.2.2
[ae029012] Requires v1.3.0
[ae5879a3] ResettableStacks v1.1.1
[79098fc4] Rmath v0.7.1
[7e49a35a] RuntimeGeneratedFunctions v0.5.13
[94e857df] SIMDTypes v0.1.0
[476501e8] SLEEFPirates v0.6.43
[0bca4576] SciMLBase v2.51.0
[c0aeaf25] SciMLOperators v0.3.10
[53ae85a6] SciMLStructures v1.5.0
[efcf1570] Setfield v1.1.1
[727e6d20] SimpleNonlinearSolve v1.12.0
[699a6c99] SimpleTraits v0.9.4
[ce78b400] SimpleUnPack v1.1.0
[a2af1166] SortingAlgorithms v1.2.1
[47a9eef4] SparseDiffTools v2.20.0
[0a514795] SparseMatrixColorings v0.4.0
[e56a9233] Sparspak v0.3.9
[276daf66] SpecialFunctions v2.4.0
[aedffcd0] Static v1.1.1
[0d7ed370] StaticArrayInterface v1.8.0
[90137ffa] StaticArrays v1.9.7
[1e83bf80] StaticArraysCore v1.4.3
[82ae8749] StatsAPI v1.7.0
[2913bbd2] StatsBase v0.34.3
[4c63d2b9] StatsFuns v1.3.1
[7792a7ef] StrideArraysCore v0.5.7
[2efcf032] SymbolicIndexingInterface v0.3.28
[19f23fe9] SymbolicLimits v0.2.2
[d1185830] SymbolicUtils v3.5.0
[0c5d862f] Symbolics v6.5.0
[3783bdb8] TableTraits v1.0.1
[bd369af6] Tables v1.12.0
[8ea1fca8] TermInterface v2.0.0
[8290d209] ThreadingUtilities v0.5.2
[a759f4b9] TimerOutputs v0.5.24
[0796e94c] Tokenize v0.5.29
[d5829a12] TriangularSolve v0.2.1
[410a4b4d] Tricks v0.1.9
[781d530d] TruncatedStacktraces v1.4.0
[5c2747f8] URIs v1.5.1
[3a884ed6] UnPack v1.0.2
[1986cc42] Unitful v1.21.0
[a7c27f48] Unityper v0.1.6
[3d5dd08c] VectorizationBase v0.21.70
[19fa3120] VertexSafeGraphs v0.2.0
[1d5cc7b8] IntelOpenMP_jll v2024.2.1+0
[856f044c] MKL_jll v2024.2.0+0
[efe28fd5] OpenSpecFun_jll v0.5.5+0
⌅ [f50d1b31] Rmath_jll v0.4.3+0
[1317d2d5] oneTBB_jll v2021.12.0+0
[0dad84c5] ArgTools v1.1.1
[56f22d72] Artifacts
[2a0f44e3] Base64
[ade2ca70] Dates
[8ba89e20] Distributed
[f43a241f] Downloads v1.6.0
[7b1f6079] FileWatching
[9fa8497b] Future
[b77e0a4c] InteractiveUtils
[4af54fe1] LazyArtifacts
[b27032c2] LibCURL v0.6.4
[76f85450] LibGit2
[8f399da3] Libdl
[37e2e46d] LinearAlgebra
[56ddb016] Logging
[d6f4376e] Markdown
[a63ad114] Mmap
[ca575930] NetworkOptions v1.2.0
[44cfe95a] Pkg v1.10.0
[de0858da] Printf
[3fa0cd96] REPL
[9a3f8284] Random
[ea8e919c] SHA v0.7.0
[9e88b42a] Serialization
[1a1011a3] SharedArrays
[6462fe0b] Sockets
[2f01184e] SparseArrays v1.10.0
[10745b16] Statistics v1.10.0
[4607b0f0] SuiteSparse
[fa267f1f] TOML v1.0.3
[a4e569a6] Tar v1.10.0
[8dfed614] Test
[cf7118a7] UUIDs
[4ec0a83e] Unicode
[e66e0078] CompilerSupportLibraries_jll v1.1.1+0
[deac9b47] LibCURL_jll v8.4.0+0
[e37daf67] LibGit2_jll v1.6.4+0
[29816b5a] LibSSH2_jll v1.11.0+1
[c8ffd9c3] MbedTLS_jll v2.28.2+1
[14a3606d] MozillaCACerts_jll v2023.1.10
[4536629a] OpenBLAS_jll v0.3.23+4
[05823500] OpenLibm_jll v0.8.1+2
[bea87d4a] SuiteSparse_jll v7.2.1+1
[83775a58] Zlib_jll v1.2.13+1
[8e850b90] libblastrampoline_jll v5.8.0+1
[8e850ede] nghttp2_jll v1.52.0+1
[3f19e933] p7zip_jll v17.4.0+2
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated -m`
- Output of
versioninfo()
Julia Version 1.10.4
Commit 48d4fd4843 (2024-06-04 10:41 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: Windows (x86_64-w64-mingw32)
CPU: 16 × 12th Gen Intel(R) Core(TM) i5-1250P
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, alderlake)
Threads: 1 default, 0 interactive, 1 GC (on 16 virtual cores)
Environment:
JULIA_EDITOR = code
JULIA_NUM_THREADS =
@YingboMa do you know what this warning is about? It's something to do with the symbolic array usage but it seems benign?
seeing this as well. the equations it warns against include a registered array function which generates a ton of screen output, which shows up during testing.
I have also noticed huge structural_simplify performance hits with this warning. I don't know if it's just due to the potentially large amount of printing to STERROR it does but simplifications can take seconds then I invoke this error and its close to an hour or sometimes longer. I have tested it where I change my vectors from length 5 to length 500 and it changes the simplify time from half a second to 90 mins. Is there a temporary way to disable the warning to fix performance
I am getting the same issue with a @register_array_symbolic function. Any updates on this issue?
@register_array_symbolic won't solve it, but likely exacerbate the problem. The workaround is to scalarize the problematic function call. If @variables x(t)[1:3] is an array unknown, function calls should have f([x[1], x[2], x[3]]) instead of f(x) if they can't be traced into. Otherwise, it's better to not register the function. This is of course very general advice, it's easier to "solve" this problem for specific examples.
The true solution for this requires better support for array symbolics in Symbolics.jl/SymbolicUtils.jl, and then some changes to structural_simplify.
@AayushSabharwal Back again, with the same error on a new problem :) Now I have this function
function calc_pos(ps, vel, pulley_l0)
prob, getter, setter, s_idxs = ps
setter(prob, [vel, pulley_l0])
sol = solve(prob, NewtonRaphson(autodiff=AutoFiniteDiff(), verbose=false); abstol=1e-5, reltol=1e-5, verbose=false)
return getter(sol)
end
@register_array_symbolic calc_pos(ps::Tuple, vel::AbstractMatrix, pulley_l0::AbstractVector) begin
size = (3, length(ps[4]))
eltype = Float64
end
Where vel is a 3xN matrix, and pulley_l0 is a vector. Not registering the function is not an option, because I understandably get this error:
ERROR: LoadError: MethodError: no method matching Float64(::Num)
But registering the function causes the huge printouts from structural_simplify. Scalarizing the problem isn't possible either, as that would mean having 3xN arguments for vel...
I want to add the calc_pos function, without it being called with Num. What should I try to solve this?
When registering the function: it prints out a giant equation first which floods the REPL, so I can only see the last part of the warning:
giant equation ... , but was actually zero
└ @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/CvDvM/src/structural_transformation/utils.jl:262
Without registering:
ERROR: LoadError: MethodError: no method matching Float64(::Num)
The type `Float64` exists, but no method is defined for this combination of argument types when trying to construct it.
Closest candidates are:
(::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
@ Base rounding.jl:265
(::Type{T})(::T) where T<:Number
@ Core boot.jl:900
Float64(::IrrationalConstants.Sqrthalfπ)
@ IrrationalConstants ~/.julia/packages/IrrationalConstants/lWTip/src/macro.jl:131
...
Stacktrace:
[1] convert(::Type{Float64}, x::Num)
@ Base ./number.jl:7
[2] setindex!(A::Vector{Float64}, x::Num, i::Int64)
@ Base ./array.jl:987
[3] macro expansion
@ ./multidimensional.jl:981 [inlined]
[4] macro expansion
@ ./cartesian.jl:64 [inlined]
[5] _unsafe_setindex!(::IndexLinear, A::Vector{Float64}, x::Symbolics.Arr{Num, 2}, I::Base.ReshapedArray{Int64, 2, UnitRange{Int64}, Tuple{}})
@ Base ./multidimensional.jl:979
[6] _setindex!
@ ./multidimensional.jl:967 [inlined]
[7] setindex!
@ ./abstractarray.jl:1413 [inlined]
[8] set_parameter!(p::MTKParameters{…}, val::Symbolics.Arr{…}, pidx::ModelingToolkit.ParameterIndex{…})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/CvDvM/src/systems/parameter_buffer.jl:357
[9] set_parameter!
@ ~/.julia/packages/SymbolicIndexingInterface/gQU1g/src/value_provider_interface.jl:67 [inlined]
[10] SetParameterIndex
@ ~/.julia/packages/SymbolicIndexingInterface/gQU1g/src/parameter_indexing.jl:678 [inlined]
[11] #44
@ ~/.julia/packages/SymbolicIndexingInterface/gQU1g/src/parameter_indexing.jl:695 [inlined]
[12] (::Base.var"#4#5"{SymbolicIndexingInterface.var"#44#45"{…}})(a::Tuple{SymbolicIndexingInterface.SetParameterIndex{…}, Symbolics.Arr{…}})
@ Base ./generator.jl:37
[13] iterate
@ ./generator.jl:48 [inlined]
[14] collect(itr::Base.Generator{Base.Iterators.Zip{Tuple{Vector{…}, Vector{…}}}, Base.var"#4#5"{SymbolicIndexingInterface.var"#44#45"{NonlinearProblem{…}}}})
@ Base ./array.jl:791
[15] map
@ ./abstractarray.jl:3495 [inlined]
[16] MultipleSetters
@ ~/.julia/packages/SymbolicIndexingInterface/gQU1g/src/parameter_indexing.jl:695 [inlined]
[17] ParameterHookWrapper
@ ~/.julia/packages/SymbolicIndexingInterface/gQU1g/src/parameter_indexing.jl:646 [inlined]
[18] calc_pos(ps::Tuple{…}, vel::Symbolics.Arr{…}, pulley_l0::Symbolics.Arr{…})
@ Main ~/Code/MWE/mwe/pulley-tethers.jl:233
[19] model(se::Settings3)
@ Main ~/Code/MWE/mwe/pulley-tethers.jl:310
[20] main()
@ Main ~/Code/MWE/mwe/pulley-tethers.jl:376
[21] top-level scope
@ ~/Code/MWE/mwe/pulley-tethers.jl:387
[22] include(fname::String)
@ Main ./sysimg.jl:38
[23] top-level scope
@ REPL[1]:1
in expression starting at /home/bart/Code/MWE/mwe/pulley-tethers.jl:386
Some type information was truncated. Use `show(err)` to see complete types.
I need the variable in the first part of that message, unfortunately. Assuming it is vel, the only real solution is to pass collect(vel) to calc_pos.
Unfortunately in that case I'm not really sure how to fix this. It's basically a limitation of our compiler.
Here is the code if you want to try:
using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, Parameters
using ModelingToolkit: t_nounits as t, D_nounits as D, setp, getu
using NonlinearSolve
struct Point
type::Symbol
position::Union{Vector{Float64}, Nothing}
velocity::Union{Vector{Float64}, Nothing}
mass::Union{Float64, Nothing}
force::Vector{Float64}
end
struct Tether
points::Tuple{Int, Int}
l0::Union{Float64, Nothing}
stiffness::Float64
damping::Float64
end
struct Pulley
tethers::Tuple{Int, Int}
sum_length::Float64
end
# protune-speed-system.jpg
points = [
Point(:fixed, [0, 0, 2], zeros(3), nothing, zeros(3)), # Fixed point
Point(:fixed, [1, 0, 2], zeros(3), nothing, zeros(3)), # Fixed point
Point(:fixed, [2, 0, 2], zeros(3), nothing, zeros(3)), # Fixed point
Point(:fixed, [5.5, 0, 2], zeros(3), nothing, zeros(3)), # Fixed point
Point(:quasi_static, [1, 0, -1], zeros(3), 1.0, zeros(3)),
Point(:quasi_static, [0.5, 0, -2], zeros(3), 1.0, zeros(3)),
Point(:dynamic, [1.5, 0, -2], zeros(3), 1.0, zeros(3)),
Point(:dynamic, [1, 0, -3], zeros(3), 1.0, zeros(3)),
Point(:dynamic, [2, 0, -3], zeros(3), 1.0, zeros(3)),
Point(:dynamic, [1, 0, -10], zeros(3), 1.0, [0., 0., -50]),
Point(:dynamic, [2, 0, -10], zeros(3), 1.0, [0., 0., -50]),
]
stiffness = 614600
damping = 4730
tethers = [
Tether((1, 6), norm(points[1].position - points[6].position), stiffness, damping),
Tether((2, 5), norm(points[2].position - points[5].position), stiffness, damping),
Tether((3, 7), norm(points[3].position - points[7].position), stiffness, damping),
Tether((4, 9), norm(points[4].position - points[9].position) - 0.5, stiffness, damping),
Tether((5, 6), norm(points[5].position - points[6].position) - 0.1, stiffness, damping),
Tether((5, 7), norm(points[5].position - points[7].position), stiffness, damping),
Tether((6, 8), norm(points[6].position - points[8].position), stiffness, damping),
Tether((7, 8), norm(points[7].position - points[8].position), stiffness, damping),
Tether((7, 9), norm(points[7].position - points[9].position), stiffness, damping),
Tether((8, 10), norm(points[8].position - points[10].position), stiffness, damping),
Tether((9, 11), norm(points[9].position - points[11].position), stiffness, damping),
]
pulleys = [
Pulley((5, 6), (tethers[5].l0 + tethers[6].l0))
Pulley((8, 9), (tethers[8].l0 + tethers[9].l0))
]
@with_kw mutable struct Settings3 @deftype Float64
g_earth::Vector{Float64} = [0.0, 0.0, -9.81] # gravitational acceleration [m/s²]
v_wind_tether::Vector{Float64} = [2, 0.0, 0.0]
rho = 1.225
cd_tether = 0.958
l0 = 50 # initial tether length [m]
v_ro = 2 # reel-out speed [m/s]
d_tether = 4 # tether diameter [mm]
rho_tether = 724 # density of Dyneema [kg/m³]
c_spring = 614600 # unit spring constant [N]
rel_compression_stiffness = 0.01 # relative compression stiffness [-]
damping = 473 # unit damping constant [Ns]
segments::Int64 = 2 # number of tether segments [-]
α0 = π/10 # initial tether angle [rad]
duration = 0.2 # duration of the simulation [s]
save::Bool = false # save png files in folder video
end
function set_tether_diameter!(se, d; c_spring_4mm = 614600, damping_4mm = 473)
se.d_tether = d
se.c_spring = c_spring_4mm * (d/4.0)^2
se.damping = damping_4mm * (d/4.0)^2
end
function calc_initial_state(points, tethers, pulleys)
POS0 = zeros(3, length(points))
VEL0 = zeros(3, length(points))
L0 = zeros(length(pulleys))
V0 = zeros(length(pulleys))
for i in eachindex(points)
POS0[:, i] .= points[i].position
VEL0[:, i] .= points[i].velocity
end
for i in eachindex(pulleys)
L0[i] = tethers[pulleys[i].tethers[1]].l0
V0[i] = 0.0
end
POS0, VEL0, L0, V0
end
function calc_spring_forces(pos::AbstractMatrix{T}, vel, pulley_l0) where T
# loop over all tethers to calculate spring forces
spring_force = zeros(T, length(tethers))
spring_force_vec = zeros(T, 3, length(tethers))
segment = zeros(T, 3)
unit_vector = zeros(T, 3)
rel_vel = zeros(T, 3)
for (tether_idx, tether) in enumerate(tethers)
found = false
for (pulley_idx, pulley) in enumerate(pulleys)
if tether_idx == pulley.tethers[1] # each tether should only be part of one pulley
l0 = pulley_l0[pulley_idx]
found = true
break
elseif tether_idx == pulley.tethers[2]
l0 = pulley.sum_length - pulley_l0[pulley_idx]
found = true
break
end
end
if !found
l0 = tether.l0
end
p1, p2 = tether.points[1], tether.points[2]
segment .= pos[:, p2] .- pos[:, p1]
len = norm(segment)
unit_vector .= segment ./ len
rel_vel .= vel[:, p1] .- vel[:, p2]
spring_vel = rel_vel ⋅ unit_vector
spring_force[tether_idx] = (tether.stiffness * tether.l0 * (len - l0) - tether.damping * tether.l0 * spring_vel)
spring_force_vec[:, tether_idx] .= spring_force[tether_idx] .* unit_vector
end
return spring_force_vec, spring_force
end
function calc_acc(se::Settings3, pos::AbstractMatrix{T}, vel, pulley_l0) where T
spring_force_vec, spring_force = calc_spring_forces(pos, vel, pulley_l0)
pulley_acc = zeros(T, length(pulleys))
for (pulley_idx, pulley) in enumerate(pulleys)
M = 3.1
pulley_force = spring_force[pulley.tethers[1]] - spring_force[pulley.tethers[2]]
pulley_acc[pulley_idx] = pulley_force / M
end
acc = zeros(T, 3, length(points))
force = zeros(T, 3)
for (point_idx, point) in enumerate(points)
if point.type === :fixed
acc[:, point_idx] .= 0.0
else
force .= 0.0
for (j, tether) in enumerate(tethers)
if point_idx in tether.points
inverted = tether.points[2] == point_idx
if inverted
force .-= spring_force_vec[:, j]
else
force .+= spring_force_vec[:, j]
end
end
end
force .+= point.force
acc[:, point_idx] .= force ./ point.mass .+ se.g_earth
end
end
return acc, pulley_acc
end
@register_symbolic calc_acc(se::Settings3, pos, vel, pulley_l0)
function create_pos_prob(se::Settings3, s_idxs, d_idxs)
POS0, VEL0, L0, V0 = calc_initial_state(points, tethers, pulleys)
@parameters begin
dynamic_pos[1:3, eachindex(d_idxs)]
vel[1:3, eachindex(points)]
pulley_l0[eachindex(pulleys)]
end
@variables begin
static_pos(t)[1:3, eachindex(s_idxs)]
static_acc(t)[1:3, eachindex(s_idxs)]
end
pos = zeros(Num, 3, length(points))
for (j, i) in enumerate(d_idxs)
for k in 1:3
pos[k, i] = dynamic_pos[k, j]
end
end
for (j, i) in enumerate(s_idxs)
for k in 1:3
pos[k, i] = static_pos[k, j]
end
end
eqs = []
for (j, i) in enumerate(s_idxs)
eqs = [
eqs
static_acc[:, j] ~ zeros(3)
static_acc[:, j] ~ calc_acc(se, pos, vel, pulley_l0)[1][:, i]
]
end
eqs = reduce(vcat, Symbolics.scalarize.(eqs))
@mtkbuild ns = NonlinearSystem(eqs, [static_pos, static_acc], [dynamic_pos, vel, pulley_l0])
u0 = [
static_pos => POS0[:, s_idxs]
]
ps = [
dynamic_pos => POS0[:, d_idxs]
vel => VEL0
pulley_l0 => L0
]
prob = NonlinearProblem(ns, u0, ps; verbose=false)
getter = getu(prob, static_pos)
setter = setp(prob, [vel, pulley_l0])
return prob, getter, setter, s_idxs
end
function calc_pos(ps, vel, pulley_l0)
prob, getter, setter, s_idxs = ps
setter(prob, [vel, pulley_l0])
sol = solve(prob, NewtonRaphson(autodiff=AutoFiniteDiff(), verbose=false); abstol=1e-5, reltol=1e-5, verbose=false)
return getter(sol)
end
@register_array_symbolic calc_pos(ps::Tuple, vel::AbstractMatrix, pulley_l0::AbstractVector) begin
size = (3, length(ps[4]))
eltype = Float64
end
function model(se::Settings3)
@parameters c_spring0=se.c_spring/(se.l0/se.segments) l_seg=se.l0/se.segments
@parameters rel_compression_stiffness = se.rel_compression_stiffness
@variables begin
pos(t)[1:3, eachindex(points)]
vel(t)[1:3, eachindex(points)]
acc(t)[1:3, eachindex(points)]
force(t)[1:3, eachindex(points)]
pulley_force(t)[eachindex(pulleys)]
pulley_l0(t)[eachindex(pulleys)] # first tether length in pulley
pulley_vel(t)[eachindex(pulleys)]
pulley_acc(t)[eachindex(pulleys)]
segment(t)[1:3, eachindex(tethers)]
unit_vector(t)[1:3, eachindex(tethers)]
l_spring(t), c_spring(t), damping(t), m_tether_particle(t)
len(t)[eachindex(tethers)]
l0(t)[eachindex(tethers)]
rel_vel(t)[1:3, eachindex(tethers)]
spring_vel(t)[eachindex(tethers)]
spring_force(t)[eachindex(tethers)]
spring_force_vec(t)[1:3, eachindex(tethers)] # spring force from spring p1 to spring p2
end
POS0, VEL0, L0, V0 = calc_initial_state(points, tethers, pulleys)
defaults = Pair{Num, Real}[]
guesses = Pair{Num, Real}[]
eqs = [
D(pulley_l0) ~ pulley_vel
D(pulley_vel) ~ pulley_acc - 10pulley_vel
]
s_idxs = Int[]
for (point_idx, point) in enumerate(points)
if point.type === :fixed
eqs = [
eqs
pos[:, point_idx] ~ point.position
vel[:, point_idx] ~ zeros(3)
]
elseif point.type === :dynamic
eqs = [
eqs
D(pos[:, point_idx]) ~ vel[:, point_idx]
D(vel[:, point_idx]) ~ acc[:, point_idx]
]
defaults = [
defaults
[pos[j, point_idx] => POS0[j, point_idx] for j in 1:3]
[vel[j, point_idx] => 0 for j in 1:3]
]
elseif point.type === :quasi_static
push!(s_idxs, point_idx)
guesses = [
guesses
[pos[j, point_idx] => POS0[j, point_idx] for j in 1:3]
[vel[j, point_idx] => 0 for j in 1:3]
]
else
println("wrong type")
end
end
d_idxs = setdiff(eachindex(points), s_idxs)
ps = create_pos_prob(se, s_idxs, d_idxs)
pos = collect(pos)
vel = collect(vel)
pulley_l0 = collect(pulley_l0)
for (j, i) in enumerate(s_idxs)
eqs = [
eqs
pos[:, i] ~ calc_pos(ps, vel, pulley_l0)[:, j]
vel[:, i] ~ zeros(3)
]
end
defaults = [
defaults
pulley_l0 => L0
pulley_vel => V0
]
eqs = [
eqs
vec(acc) ~ vec(calc_acc(se, pos, vel, pulley_l0)[1])
vec(pulley_acc) ~ vec(calc_acc(se, pos, vel, pulley_l0)[2])
]
eqs = reduce(vcat, Symbolics.scalarize.(eqs))
@named sys = ODESystem(eqs, t)
println("struct simplify")
@time sys = structural_simplify(sys; simplify=false)
sys, pos, vel, defaults, guesses
end
se = Settings3()
sys, pos, vel, defaults, guesses = model(se)
nothing
Are there any available alternatives I can try for a nested solver in a ModelingToolkit ODE, as this method does not work? Or is there a possibility that this will be fixed?
Edit: I figured out now that I can just set the acceleration to zero... Problem fixed / I don't need the registered function for now.
Hi I have a problem with a function, vol ~ vol_pTz(model,p,T,zᵢ), that has a vector input and has the issue of warning "Variable (vol₊zᵢ(t))[4] was marked as being in ... but was actually zero". I can't (don't want to) change my function to a list of scalars as my vector can change by the user input and can be non trival in size. This looks benign as my model runs corrrectly with this warning is there a way to supress the warning like we have for warn_initialize_determined = false
BTW the function is defined as:
function vol_pTz(model,p,T,z) return volume(model,p,T,z) end @register_symbolic vol_pTz(model::EoSModel,p::Float64,T::Float64,z::Array{Float64, 1})
where volume(model,p,T,z) is a Calpeyron function
Hi has any progress been made with this warning.
Could we get a way of switch off the warning or writing to the terminal. At the moment I happy to live with the knowledge but now want to ignore as it’s not affecting my solution.
I have functions that are registered as symbolic arrays. All looks correct in the model definition and solution. It’s a pretty simple steady state. However. I have 100s of warnings that will only grow as I make my problems more complex
Part of the issue I have is my functions need to be registered symbolic otherwise MTK complains. I then am forced to use scalarize on the vector return to set up an equation for each item in the vector. This is then the root of all the warnings. Each one of the equations made by scalerize has this warning.
ok I have wrapped my MTK code in @suppress_err local …. MTK code end
this stops the annoying warnings.