julia> prob = ODEProblem(ssys, x0, (0.0, 5.0); guesses)
ERROR: BoundsError: attempt to access 0-element Vector{Float64} at index [1]
Stacktrace:
[1] throw_boundserror(A::Vector{Float64}, I::Tuple{Int64})
@ Base ./essentials.jl:15
[2] getindex(A::Vector{Float64}, i::Int64)
@ Base ./essentials.jl:919
[3] scalar_index(val::Any, idx::Int64)
@ SymbolicUtils ~/.julia/packages/SymbolicUtils/QB5bd/src/symbolic_ops/getindex.jl:197
[4] __stable_getindex(arr::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{…}, sidxs::SymbolicUtils.StableIndex{…})
@ SymbolicUtils ~/.julia/packages/SymbolicUtils/QB5bd/src/symbolic_ops/getindex.jl:227
[5] var"##_stable_getindex_1#774"(arr::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{…}, sidxs::SymbolicUtils.StableIndex{…})
@ SymbolicUtils ~/.julia/packages/SymbolicUtils/QB5bd/src/symbolic_ops/getindex.jl:208
[6] _stable_getindex_1(arr::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{…}, sidxs::SymbolicUtils.StableIndex{…})
@ SymbolicUtils ~/.julia/packages/SymbolicUtils/QB5bd/src/cache.jl:432
[7] getindex
@ ~/.julia/packages/SymbolicUtils/QB5bd/src/symbolic_ops/getindex.jl:201 [inlined]
[8] Fix
@ ./operators.jl:1193 [inlined]
[9] call_composed
@ ./operators.jl:1100 [inlined]
[10] call_composed
@ ./operators.jl:1099 [inlined]
[11] ComposedFunction
@ ./operators.jl:1096 [inlined]
[12] _any(f::ComposedFunction{Base.Fix2{…}, Base.Fix1{…}}, itr::SymbolicUtils.StableIndices, ::Colon)
@ Base ./anyall.jl:124
[13] any
@ ./anyall.jl:115 [inlined]
[14] (::ModelingToolkitBase.var"#795#796")(v::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal})
@ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/x7JKa/src/systems/problem_utils.jl:1307
[15] map!(f::ModelingToolkitBase.var"#795#796", iter::Base.ValueIterator{ModelingToolkitBase.AtomicArrayDict{…}})
@ Base ./abstractdict.jl:680
[16] build_operating_point(sys::System, op::Dict{Any, Any}; fast_path::Bool)
@ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/x7JKa/src/systems/problem_utils.jl:1304
[17] build_operating_point
@ ~/.julia/packages/ModelingToolkitBase/x7JKa/src/systems/problem_utils.jl:1293 [inlined]
[18] process_SciMLProblem(constructor::Type, sys::System, op::Vector{…}; build_initializeprob::Bool, implicit_dae::Bool, t::Float64, guesses::Vector{…}, warn_initialize_determined::Bool, initialization_eqs::Vector{…}, eval_expression::Bool, eval_module::Module, fully_determined::Nothing, check_initialization_units::Bool, u0_eltype::Nothing, tofloat::Bool, u0_constructor::typeof(identity), p_constructor::typeof(identity), check_length::Bool, symbolic_u0::Bool, warn_cyclic_dependency::Bool, circular_dependency_max_cycle_length::Int64, circular_dependency_max_cycles::Int64, initsys_mtkcompile_kwargs::@NamedTuple{}, substitution_limit::Int64, use_scc::Bool, time_dependent_init::Bool, algebraic_only::Bool, missing_guess_value::ModelingToolkitBase.MissingGuessValue.var"typeof(MissingGuessValue)", allow_incomplete::Bool, is_initializeprob::Bool, kwargs::@Kwargs{…})
@ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/x7JKa/src/systems/problem_utils.jl:1369
[19] (ODEProblem{…})(sys::System, op::Vector{…}, tspan::Tuple{…}; callback::Nothing, check_length::Bool, eval_expression::Bool, expression::Type, eval_module::Module, check_compatibility::Bool, kwargs::@Kwargs{…})
@ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/x7JKa/src/problems/odeproblem.jl:95
[20] ODEProblem
@ ~/.julia/packages/ModelingToolkitBase/x7JKa/src/problems/odeproblem.jl:87 [inlined]
[21] ODEProblem
@ ./none:-1 [inlined]
[22] #ODEProblem#822
@ ./none:-1 [inlined]using ModelingToolkit
using Multibody
import ModelingToolkit: t_nounits as t, D_nounits as D
import ModelingToolkitStandardLibrary.Mechanical.Rotational
using OrdinaryDiffEq
@component function DyadBot3DWheels(; name)
track = 0.13#, [description = "Distance between wheels"]
body_height = 0.1#, [description = "Height of body above wheel axle"]
pars = @parameters begin
wheel_radius = 0.04
body_mass = 0.1
wheel_mass = 0.05
wheel_I_axis = 5e-5
wheel_I_long = 1e-5
d_wheel = 0.1, [description = "Wheel rotational damping coefficient"]
end
systems = @named begin
world = World()
# Two individual rolling wheels (allows tilting)
wheel_left = SlippingWheel(
radius = wheel_radius,
m = wheel_mass,
I_axis = wheel_I_axis,
I_long = wheel_I_long,
state = false,
)
wheel_right = SlippingWheel(
radius = wheel_radius,
m = wheel_mass,
I_axis = wheel_I_axis,
I_long = wheel_I_long,
# iscut = true, # Avoid over-constraining
state = true,
)
# Revolute joints for wheel spin
revolute_left = Revolute(n = [0, 0, 1], axisflange = true, phi0=0, w0=0)
revolute_right = Revolute(n = [0, 0, 1], axisflange = true, phi0=0, w0=0)
# Axis connecting wheels
rod_left = FixedTranslation(r = [0, 0, track/2])
rod_right = FixedTranslation(r = [0, 0, -track/2])
# Axis body (small mass at center)
axis_body = Body(m = 0.01, r_cm = [0, 0, 0])
# Main body extending upward
body = BodyShape(
m = body_mass,
I_22 = 0.01*0.03^2,
I_11 = 0.01*0.05^2,
I_33 = 0.01*0.05^2,
r = [0, body_height, 0],
)
# Rotational damping for wheel joints
damper_left = Rotational.Damper(d = d_wheel)
damper_right = Rotational.Damper(d = d_wheel)
end
eqs = [
# Connect wheels to revolute joints
connect(wheel_left.frame_a, revolute_left.frame_b)
connect(wheel_right.frame_a, revolute_right.frame_b)
# Connect revolute joints to axis via rods
connect(revolute_left.frame_a, rod_left.frame_b)
connect(revolute_right.frame_a, rod_right.frame_b)
# Connect rods to axis center
connect(rod_left.frame_a, axis_body.frame_a)
connect(rod_right.frame_a, axis_body.frame_a)
# Connect body to axis center
connect(axis_body.frame_a, body.frame_a)
# Connect dampers to wheel revolute joints
connect(damper_left.flange_a, revolute_left.axis)
connect(damper_left.flange_b, revolute_left.support)
connect(damper_right.flange_a, revolute_right.axis)
connect(damper_right.flange_b, revolute_right.support)
]
System(eqs, t, [], pars; systems, name)
end
##
@named model_wheels = DyadBot3DWheels()
model_wheels = complete(model_wheels)
ssys = multibody(model_wheels)
x0 = [
# collect(ssys.axis_body.Q̂) .=> [1, 0, 0, 0];
]
guesses = unknowns(ssys) .=> 0.0
prob = ODEProblem(ssys, x0, (0.0, 5.0); guesses)
