namespace updates

Author: baggepinnen

docs/src/examples/pendulum.md

@@ -65,7 +65,7 @@ D = Differential(t)
 defs = Dict() # We may specify the initial condition here
 prob = ODEProblem(ssys, defs, (0, 3.35))
 
-sol = solve(prob, Rodas4())
+sol = solve(prob, Tsit5())
 plot(sol, idxs = joint.phi, title="Pendulum")
 ```
 The solution `sol` can be plotted directly if the Plots package is loaded. The figure indicates that the pendulum swings back and forth without any damping. To add damping as well, we could add a `Damper` from the `ModelingToolkitStandardLibrary.Mechanical.Rotational` module to the revolute joint. We do this below
@@ -74,7 +74,7 @@ The solution `sol` can be plotted directly if the Plots package is loaded. The f
 Multibody.jl supports automatic 3D rendering of mechanisms, we use this feature to illustrate the result of the simulation below:
 
 ```@example pendulum
-import GLMakie # GLMakie is another alternative, suitable for interactive plots
+import GLMakie
 Multibody.render(model, sol; filename = "pendulum.gif") # Use "pendulum.mp4" for a video file
 nothing # hide
 ```
@@ -95,12 +95,11 @@ connections = [connect(world.frame_b, joint.frame_a)
                connect(body.frame_a, joint.frame_b)]
 
 @named model = System(connections, t, systems = [world, joint, body, damper])
-model = complete(model)
 ssys = multibody(model)
 
-prob = ODEProblem(ssys, [damper.phi_rel => 1], (0, 10))
+prob = ODEProblem(ssys, [], (0, 10))
 
-sol = solve(prob, Rodas4())
+sol = solve(prob, Tsit5())
 plot(sol, idxs = joint.phi, title="Damped pendulum")
 ```
 This time we see that the pendulum loses energy and eventually comes to rest at the stable equilibrium point ``\pi / 2``.
@@ -128,12 +127,11 @@ connections = [connect(world.frame_b, joint.frame_a)
                connect(body_0.frame_a, joint.frame_b)]
 
 @named model = System(connections, t, systems = [world, joint, body_0, damper, spring])
-model = complete(model)
 ssys = multibody(model)
 
-prob = ODEProblem(ssys, [], (0, 10))
+prob = ODEProblem(ssys, [ssys.joint.s => 0, ssys.joint.v => 0], (0, 10))
 
-sol = solve(prob, Rodas4())
+sol = solve(prob, Tsit5())
 Plots.plot(sol, idxs = joint.s, title="Mass-spring-damper system")
 ```
 
@@ -155,14 +153,13 @@ connections = [connect(world.frame_b, multibody_spring.frame_a)
                 connect(root_body.frame_a, multibody_spring.frame_b)]
 
 @named model = System(connections, t, systems = [world, multibody_spring, root_body])
-model = complete(model)
 ssys = multibody(model)
 
-defs = Dict(root_body.r_0 .=> [0, 1e-3, 0]) # The spring has a singularity at zero length, so we start some distance away
+defs = Dict(collect(ssys.root_body.r_0) .=> [0, 1e-3, 0]) # The spring has a singularity at zero length, so we start some distance away
 
 prob = ODEProblem(ssys, defs, (0, 10))
 
-sol = solve(prob, Rodas4())
+sol = solve(prob, Tsit5())
 plot(sol, idxs = multibody_spring.r_rel_0[2], title="Mass-spring system without joint")
 ```
 Here, we used a [`Multibody.Spring`](@ref) instead of connecting a `Translational.Spring` to a joint. The `Translational.Spring`, alongside other components from `ModelingToolkitStandardLibrary.Mechanical`, is a 1-dimensional object, whereas multibody components are 3-dimensional objects.
@@ -173,7 +170,6 @@ push!(connections, connect(multibody_spring.spring2d.flange_a, damper.flange_a))
 push!(connections, connect(multibody_spring.spring2d.flange_b, damper.flange_b))
 
 @named model = System(connections, t, systems = [world, multibody_spring, root_body, damper])
-model = complete(model)
 ssys = multibody(model)
 prob = ODEProblem(ssys, defs, (0, 10))
 
@@ -224,18 +220,17 @@ import Multibody.Rotations
 end
 
 @named model = FurutaPendulum()
-model = complete(model)
 ssys = multibody(model)
 
-prob = ODEProblem(ssys, [model.shoulder_joint.phi => 0.0, model.elbow_joint.phi => 0.1], (0, 10))
-sol = solve(prob, Rodas4())
+prob = ODEProblem(ssys, [ssys.shoulder_joint.phi => 0.0, ssys.elbow_joint.phi => 0.1], (0, 10))
+sol = solve(prob, Tsit5())
 plot(sol, layout=4)
 ```
 
 In the animation below, we visualize the path that the origin of the pendulum tip traces by providing the tip frame in a vector of frames passed to `traces`
 ```@example pendulum
 import GLMakie
-Multibody.render(model, sol, filename = "furuta.gif", traces=[model.tip.frame_a])
+Multibody.render(model, sol, filename = "furuta.gif", traces=[ssys.tip.frame_a])
 nothing # hide
 ```
 ![furuta](furuta.gif)
@@ -245,11 +240,11 @@ nothing # hide
 Let's break down how to think about directions and orientations when building 3D mechanisms. In the example above, we started with the shoulder joint, this joint rotated around the gravitational axis, `n = [0, 1, 0]`. When this joint is positioned in joint coordinate `shoulder_joint.phi = 0`, its `frame_a` and `frame_b` will coincide. When the joint rotates, `frame_b` will rotate around the axis `n` of `frame_a`. The `frame_a` of the joint is attached to the world, so the joint will rotate around the world's `y`-axis:
 
 ```@example pendulum
-get_rot(sol, model.shoulder_joint.frame_b, 0)
+get_rot(sol, ssys.shoulder_joint.frame_b, 0)
 ```
 we see that at time $t = 0$, we have no rotation of `frame_b` around the $y$ axis of the world (frames are always resolved in the world frame), but a second into the simulation, we have:
 ```@example pendulum
-R1 = get_rot(sol, model.shoulder_joint.frame_b, 1)
+R1 = get_rot(sol, ssys.shoulder_joint.frame_b, 1)
 ```
 Here, the `frame_b` has rotated around the $y$ axis of the world (if you are not familiar with rotation matrices, we can ask for the rotation axis and angle)
 ```@example pendulum
@@ -259,7 +254,7 @@ rotation_axis(R1), rotation_angle(R1)
 
 This rotation axis and angle should correspond to the joint coordinate (the orientation described by an axis and an angle is invariant to a multiplication of both by -1)
 ```@example pendulum
-sol(1, idxs=model.shoulder_joint.phi)
+sol(1, idxs=ssys.shoulder_joint.phi)
 ```
 
 !!! note "Convention"
@@ -271,25 +266,25 @@ Here, we made use of the function [`get_rot`](@ref), we will now make use of als
 
 The next body is the upper arm. This body has an extent of `0.6` in the $z$ direction, as measured in its local `frame_a`
 ```@example pendulum
-get_trans(sol, model.upper_arm.frame_b, 0)
+get_trans(sol, ssys.upper_arm.frame_b, 0)
 ```
 One second into the simulation, the upper arm has rotated around the $y$ axis of the world
 ```@example pendulum
-rb1 = get_trans(sol, model.upper_arm.frame_b, 1)
+rb1 = get_trans(sol, ssys.upper_arm.frame_b, 1)
 ```
 
-If we look at the variable `model.upper_arm.r`, we do not see this rotation!
+If we look at the variable `ssys.upper_arm.r`, we do not see this rotation!
 ```@example pendulum
-arm_r = sol(1, idxs=model.upper_arm.r)
+arm_r = sol(1, idxs=ssys.upper_arm.r)
 ```
 The reason is that this variable is resolved in the local `frame_a` and not in the world frame. To transform this variable to the world frame, we may multiply with the rotation matrix of `frame_a` which is always resolved in the world frame:
 ```@example pendulum
-get_rot(sol, model.upper_arm.frame_a, 1)*arm_r
+get_rot(sol, ssys.upper_arm.frame_a, 1)*arm_r
 ```
 We now get the same result has when we asked for the translation vector of `frame_b` above.
 ```@example pendulum
 using Test # hide
-get_rot(sol, model.upper_arm.frame_a, 1)*arm_r ≈ rb1 # hide
+get_rot(sol, ssys.upper_arm.frame_a, 1)*arm_r ≈ rb1 # hide
 nothing # hide
 ```
 
@@ -300,23 +295,23 @@ The next joint, the elbow joint, has the rotational axis `n = [0, 0, 1]`. This i
 
 The lower arm is finally having an extent along the $y$-axis. At the final time when the pendulum motion has been fully damped, we see that the second frame of this body ends up with an $y$-coordinate of `-0.6`:
 ```@example pendulum
-get_trans(sol, model.lower_arm.frame_b, 12)
+get_trans(sol, ssys.lower_arm.frame_b, 12)
 ```
 
 If we rotate the vector of extent of the lower arm to the world frame, we indeed see that the only coordinate that is nonzero is the $y$ coordinate:
 ```@example pendulum
-get_rot(sol, model.lower_arm.frame_a, 12)*sol(12, idxs=model.lower_arm.r)
+get_rot(sol, ssys.lower_arm.frame_a, 12)*sol(12, idxs=ssys.lower_arm.r)
 ```
 
-The reason that the latter vector differs from `get_trans(sol, model.lower_arm.frame_b, 12)` above is that `get_trans(sol, model.lower_arm.frame_b, 12)` has been _translated_ as well. To both translate and rotate `model.lower_arm.r` into the world frame, we must use the full transformation matrix $T_W^A \in SE(3)$:
+The reason that the latter vector differs from `get_trans(sol, ssys.lower_arm.frame_b, 12)` above is that `get_trans(sol, ssys.lower_arm.frame_b, 12)` has been _translated_ as well. To both translate and rotate `ssys.lower_arm.r` into the world frame, we must use the full transformation matrix $T_W^A \in SE(3)$:
 
 ```@example pendulum
-r_A = sol(12, idxs=model.lower_arm.r)
+r_A = sol(12, idxs=ssys.lower_arm.r)
 r_A = [r_A; 1] # Homogeneous coordinates
 
-get_frame(sol, model.lower_arm.frame_a, 12)*r_A
+get_frame(sol, ssys.lower_arm.frame_a, 12)*r_A
 ```
-the vector is now coinciding with `get_trans(sol, model.lower_arm.frame_b, 12)`.
+the vector is now coinciding with `get_trans(sol, ssys.lower_arm.frame_b, 12)`.
 
 
 ## Control-design example: Pendulum on cart
@@ -336,12 +331,7 @@ using Plots
 gray = [0.5, 0.5, 0.5, 1]
 @component function Cartpole(; name, use_world = false)
     systems = @named begin
-        if use_world
-            fixed = World()
-        else
-            # In case we wrap this model in an outer model below, we place the world there instead
-            fixed = Fixed()
-        end
+        fixed = Fixed()
         cart = BodyShape(m = 1, r = [0.2, 0, 0], color=[0.2, 0.2, 0.2, 1], shape="box")
         mounting_point = FixedTranslation(r = [0.1, 0, 0])
         prismatic = Prismatic(n = [1, 0, 0], axisflange = true, color=gray, state_priority=100)
@@ -352,7 +342,7 @@ gray = [0.5, 0.5, 0.5, 1]
     end
 
     vars = @variables begin
-        u(t) = 0
+        u(t)
         x(t)
         v(t)
         phi(t)
@@ -374,7 +364,7 @@ gray = [0.5, 0.5, 0.5, 1]
         w ~ revolute.w
     ]
 
-    return System(equations, t; name, systems)
+    return System(equations, t, vars, []; name, systems)
 end
 @component function CartWithInput(; name)
     systems = @named begin
@@ -390,16 +380,15 @@ end
     return System(equations, t; name, systems)
 end
 @named model = CartWithInput()
-model = complete(model)
 ssys = multibody(model)
-prob = ODEProblem(ssys, [model.cartpole.prismatic.s => 0.0, model.cartpole.revolute.phi => 0.1], (0, 10))
+prob = ODEProblem(ssys, [ssys.cartpole.prismatic.s => 0.0, ssys.cartpole.prismatic.v => 0.0, ssys.cartpole.revolute.phi => 0.1], (0, 10))
 sol = solve(prob, Tsit5())
 plot(sol, layout=4)
 ```
 As usual, we render the simulation in 3D to get a better feel for the system:
 ```@example pendulum
 import GLMakie
-Multibody.render(model, sol, filename = "cartpole.gif", traces=[model.cartpole.pendulum.frame_b])
+Multibody.render(model, sol, filename = "cartpole.gif", traces=[ssys.cartpole.pendulum.frame_b])
 nothing # hide
 ```
 ![cartpole](cartpole.gif)
@@ -413,7 +402,7 @@ We start by linearizing the model in the upward equilibrium position using the f
 ```@example pendulum
 import ModelingToolkit: D_nounits as D
 using LinearAlgebra
-@named cp = Cartpole(use_world = true)
+@named cp = Cartpole()
 cp = complete(cp)
 inputs = [cp.u] # Input to the linearized system
 outputs = [cp.x, cp.phi, cp.v, cp.w] # These are the outputs of the linearized system
@@ -422,13 +411,13 @@ op = Dict([ # Operating point to linearize in
     cp.revolute.phi => 0 # Pendulum pointing upwards
 ]
 )
-matrices, simplified_sys = linearize(multibody(cp), inputs, outputs; op)
+matrices, simplified_sys = linearize(cp, inputs, outputs; op)
 matrices
 ```
 This gives us the matrices $A,B,C,D$ in a linearized statespace representation of the system. To make these easier to work with, we load the control packages and call `named_ss` instead of `linearize` to get a named statespace object instead:
 ```@example pendulum
 using ControlSystemsMTK
-lsys = named_ss(multibody(cp), inputs, outputs; op) # identical to linearize, but packages the resulting matrices in a named statespace object for convenience
+lsys = named_ss(cp, inputs, outputs; op) # identical to linearize, but packages the resulting matrices in a named statespace object for convenience
 ```
 
 ### LQR and LQG Control design
@@ -486,13 +475,12 @@ LQGSystem(args...; kwargs...) = System(Ce; kwargs...)
     return System(equations, t; name, systems)
 end
 @named model = CartWithFeedback()
-model = complete(model)
 ssys = multibody(model)
-prob = ODEProblem(ssys, [model.cartpole.prismatic.s => 0.1, model.cartpole.revolute.phi => 0.35], (0, 12))
+prob = ODEProblem(ssys, [ssys.cartpole.prismatic.s => 0.1, ssys.cartpole.revolute.phi => 0.35], (0, 12))
 sol = solve(prob, Tsit5())
-cp = model.cartpole
+cp = ssys.cartpole
 plot(sol, idxs=[cp.prismatic.s, cp.revolute.phi, cp.motor.f.u], layout=3)
-plot!(sol, idxs=model.reference.output.u, sp=1, l=(:black, :dash), legend=:bottomright)
+plot!(sol, idxs=ssys.reference.output.u, sp=1, l=(:black, :dash), legend=:bottomright)
 ```
 
 ```@example pendulum
@@ -557,14 +545,13 @@ normalize_angle(x::Number) = mod(x+3.1415, 2pi)-3.1415
     return System(equations, t; name, systems)
 end
 @named model = CartWithSwingup()
-model = complete(model)
 ssys = multibody(model)
-cp = model.cartpole
+cp = ssys.cartpole
 prob = ODEProblem(ssys, [cp.prismatic.s => 0.0, cp.revolute.phi => 0.99pi], (0, 5))
 sol = solve(prob, Tsit5(), dt = 1e-2, adaptive=false)
-plot(sol, idxs=[cp.prismatic.s, cp.revolute.phi, cp.motor.f.u, model.E], layout=4)
+plot(sol, idxs=[cp.prismatic.s, cp.revolute.phi, cp.motor.f.u, ssys.E], layout=4)
 hline!([0, 2pi], sp=2, l=(:black, :dash), primary=false)
-plot!(sol, idxs=[model.switching_condition], sp=2)
+plot!(sol, idxs=[ssys.switching_condition], sp=2)
 ```
 
 

src/orientation.jl

@@ -372,12 +372,12 @@ function resolve_dyade2(R, D1)
     R*D1*R'
 end
 
-function is_frame(x)
+function is_frame(frame)
     md = get_metadata(frame)
     get(md, Multibody.ModelingToolkit.IsFrame, false)
 end
 
-function is_frame_2d(x)
+function is_frame_2d(frame)
     md = get_metadata(frame)
     get(md, Multibody.PlanarMechanics.IsFrame2D, false)
 end