(usage_page)= # Usage Guide ## Basics (cpp_python_sec)= ### C++ and Python layers The core software for {{ kdflex }} is in `C++` (`C++20` standard) and uses an object-oriented design. Additionally `Python` bindings that mirror the `C++` API are available as well. The feature complete `Python` bindings provide an excellent avenue for developing and prototyping simulation software. The Python interface is especially useful for scripting, interacting with the simulation and introspection of the simulation content. This combination of languages provides a best-of-both-worlds, where users can setup and introspect their simulation using Python while obtaining the high performance of compiled `C++`. The `C++` object public interfaces are mirrored at the `Python` level and can be used for scripting and interacting with the software. The `Python` modules are included in the `Karana` Python package. Simulations are not required to use `Python` at all, and can be written entirely in `C++` if desired. When using the `C++` API, it is recommended that users build with `cmake`, as this provides the simplest way to compile/link code correctly. A sample `CMakeLists.txt` file is provided below. ```cmake # CMakeLists.txt for building script.cpp cmake_minimum_required(VERSION 3.25) project(my_project) # Find the Karana package set(Eigen3_DIR /path/to/eigen3/cmake/files) # This should be eigen 5.0 find_package(Karana REQUIRED) # Optional: Enable compile commands. This is used by LSPs. set(CMAKE_EXPORT_COMPILE_COMMANDS ON) # Example binary add_executable(my_bin main.cc) target_link_libraries(1_link Karana::kdflex) ``` It is important to note that the version of `Eigen` used is 5.0.0. Most package managers have a version of `Eigen` that is outdated by a few years and is missing critical bug fixes. Users should clone the `Eigen` package at this commit using ```bash git clone --depth 1 --branch 5.0.0 --single-branch https://gitlab.com/libeigen/eigen.git eigen_5.0.0 ``` The package can be built and installed using `cmake`. For example, by using ```bash cmake -DCMAKE_INSTALL_PREFIX=/path/to/desired/install/dir -B build . cd build make install -j 10 ``` After installation, the `Eigen` cmake files will be available at the installed location under `share/eigen3/cmake`. Of course, `C++` applications can be built without `cmake` if desired. When doing so, users will need to include the header files in `/usr/include/Karana` and link against `libkdflex.so` in `/usr/lib`. In addition, they will need to include headers for `Eigen` (again, these need to be the header files of `Eigen` 5.0) and `spdlog`, as well as link against `libspdlog` and the `libfmt` libraries. Furthermore, they will need to use the following compiler flags: `-DFMT_SHARED -DSPDLOG_COMPILED_LIB -DSPDLOG_FMT_EXTERNAL -DSPDLOG_SHARED_LIB -std=gnu++20 -msse -msse2 -mssse3 -msse4.1 -msse4.2 -mavx -mavx2`. (create_discard_sec)= ### Creating and discarding objects The general convention in {{ kdflex }} is to use `create()` static factory methods to create new instances of objects - instead of calling the constructor. Virtually all of the classes have implementations of this method. So to create an instance of class `A` call the static `A::create()` method. On the other hand, to discard an object instance you just need to call the object's `obj->discard()` method (instead of calling the destructor). Thus an instance of the {{ physbody }} can be created by calling {py:meth}`Karana.Dynamics.PhysicalBody.create` method, and can be discarded by calling the {py:meth}`Karana.Core.discard` method. ### Usage tracking of objects {{ kdflex }} supports configuration changes to the system to support the simulation of complex scenarios where multibody systems can undergo topological changes from the creation and removal of bodies, nodes, hinges etc. {{ kdflex }} uses _shared pointers_ and their associated reference counting to avoid _use after free_ errors. However, this is not sufficient since having dependencies linger past the stage where the object is no longer supposed to exist introduces vulnerabilities. {{ kdflex }} implements more rigorous usage tracking that requires all dependencies on an object to have been removed before the object can be discarded. Enforcement of such strict housekeeping can take a little get used to, but helps ensure that complex configuration changes can be carried out robustly at run-time within simulations. Such usage tracking is applied at both the `C++` and `Python` levels. The example simulations such as [2 link pendulum (basic)](generated/notebooks/example_2_link_pendulum/notebook.ipynb) include a `cleanup` section at the end where objects are explicitly discarded to ensure clean shutdown at the end of a simulation. In the event that a user runs into issues discarding objects and cannot find a way forward, assistance from a Karana employee is always available via the [forum](https://forum.karanadyn.com/) or [tickets](https://auth.karanadyn.com/) system. However, in the case that a user must move forward immediately and cannot wait for assistance, the {py:attr}`Karana.Core.DebugManager.enable_usage_tracking_map_errors` can be used to disable usage tracking: see the {ref}`DebugManager_sec` section for details. ### The dump() method Most {{ kdflex }} classes implement specializations of the {py:meth}`Karana.Core.Base.dump` method. This method can be used for introspection to display useful information specifically about an object instance. (math_sec)= ### Math layer {{ kdflex }} builds on the [Eigen](https://eigen.tuxfamily.org) math library basic classes, and extends and sub-classes them as needed. The derived classes are available under the `Karana::Math` namespace. The {cpp:type}`Karana::Math::Vec`, {cpp:type}`Karana::Math::Vec3`, {cpp:type}`Karana::Math::Mat`, {cpp:type}`Karana::Math::Mat33` etc are simple aliases around the basic `Eigen` array classes. While vectors are matrices at the `C++` level correspond to the `Eigen` classes, their corresponding `Python` equivalents are `numpy` arrays. Methods that take and return `Eigen` arrays at the `C++` level, take and return `numpy` arrays at the Python level. This by far is the largest difference between the `C++` and `Python` APIs. Otherwise, with few exceptions, the Python methods consistently replicate the C++ class and method interfaces. The Math layer also defines a special NaN value, {cpp:member}`Karana::Math::notReadyNaN`, for use in places where a numerical variable is needed but its value is has not yet been initialized. The function {cpp:func}`Karana::Math::isNotReadyNaN` checks if a given number has this special value. When not ready, the math objects in {{ kdflex }} will be populated with this special value to avoid undefined behavior or confusion with NaNs resulting from bad numerics. - - - (framestop_sec)= ## Frames Layer (attitude_sec)= ### Pose representations {{ kdflex }} supports a number of types of attitude representations including {py:class}`Karana.Math.RotationMatrix` for direction cosine rotation matrices, {py:class}`Karana.Math.RotationVector` for rotation vectors, and {py:class}`Karana.Math.UnitQuaternion` for unit quaternions. The {py:class}`Karana.Math.EulerAngles` attitude representation at the Python level and `Eigen::AngleAxis` and `Eigen::EulerAngles` attitude representations at the C++ level are also available for use. Unit quaternions are the ones used most broadly in {{ kdflex }} since they are free of singularities and computationally more efficient. Methods are available to transform between these different representations. While unit quaternions keep the vector and scalar components separately, in cases where they are combined, the convention is to keep the scalar component last after the vector elements. {{ kdflex }} uses the _passive_ approach for attitude representations. Thus the product ${}^a v = {}^a R_b {}^b v$ denotes the conversion of the the coordinates of the $v$ vector from the $b$ frame to the $a$ frame. ```{eval-rst} .. figure:: images/rotations.png :width: 500px :alt: Rotations Accumulation of rotations across frames ``` Since unit quaternions are not a minimal attitude representation, the minimal coordinate {{ rotveccls }} attitude representations are used instead as coordinates when parameterizing the relative attitude for {py:class}`Karana.Dynamics.SphericalSubhinge` instances for use with numerical integrators. Global minimal coordinate attitude representations are know to have singularities, and the {{ rotveccls }} attitude representation is no exception. Switching and re-centering charts is a way when in the neighborhood of singularities is a way to handle the singularities in practice. The {py:meth}`sanitizeCoords() ` method can be called for any {{ subhingecls }} to re-center its charts. The {py:class}`Karana.Math.HomTran` class is for homogeneous transforms that represent the relative relative pose across frames. Homogeneous transforms include both positional and orientation components. ```{eval-rst} .. _homtranDef_anchor: .. figure:: images/homtran.png :width: 400px :align: center :alt: Homogeneous transform The homogeneous transform ``` Unit quaternions are used for the attitude representation within homogeneous transforms. (spatial_notation_sec)= ### Spatial notation {{ kdflex }} uses spatial notation throughout the code base. Spatial vectors are 6 dimensional quantities with angular and linear components. The {py:class}`Karana.Math.SpatialVector` class is used for spatial vector instances. ```{eval-rst} .. _spVel: .. figure:: images/spatial_vectors.png :width: 500px :align: center :alt: Spatial vectors Spatial velocity and spatial force vectors ``` By convention, the angular component comes first within a {{ spveccls }}, followed by the the linear component. Thus spatial velocities contain angular velocities followed by linear velocities, and spatial forces have the momentum component first, followed by the linear force component. Spatial notation also extends to spatial inertias which combine the inertia tensor, with the CM offset vector and the body mass. The {py:class}`Karana.Math.SpatialInertia} class is available for spatial inertia instances. For the 3x3 rotational inertia tensors, the {{ kdflex }} convention is to use the negative integral convention for the off-diagonal terms, so that the matrix is always positive semi-definite. ```{eval-rst} .. _spinertiaDef: .. figure:: images/spinertia.png :width: 500px :align: center :alt: Spatial inertia matrix The spatial inertial matrix ``` A noteworthy fact is that {{ homtrancls }} classes have the {py:meth}`phi() ` and {py:meth}`phiStar() ` methods for carrying out rigid body transformations of spatial quantities. ```{eval-rst} .. _phiDef: .. figure:: images/phimatrix.png :width: 500px :align: center :alt: The "phi" rigid body transformation matrix The `phi` rigid body transformation matrix ``` The following summarizes expressions for rigid body transformations of spatial velocities, accelerations, forces, momenta and inertias that can be carried succinctly using these methods. ```{eval-rst} .. _phiRig: .. figure:: images/phi.png :width: 600px :align: center :alt: Rigid body transformations Rigid body transformation examples ``` - - - (frames_layer_sec)= ### Coordinate Frames {{ kdflex }} has a general purpose frames layer consisting of a directed tree of frames. The frame family of classes form a foundational layer for the multibody classes. Each {py:class}`Karana.Frame.Frame` instance represents a node in a tree of frames managed by a {py:class}`Karana.Frame.FrameContainer` instance. A {py:class}`Karana.Frame.FrameToFrame` instance is defined by a pair of from/to {{ framecls }} instances. The edges in the frame tree are defined by {py:class}`Karana.Frame.EdgeFrameToFrame` instances. Each edge contains - the relative pose across the edge defined by a {{ homtrancls }} - the relative spatial velocity across the edge defined by a {{ spveccls }} as observed from and represented in the oframe - the relative spatial acceleration across the edge defined by a {{ spveccls }} as observed from and represented in the oframe An application is responsible for setting the pose, spatial velocity and spatial acceleration for each of the edges in the frames tree. The frames layer provides support for creating {{ f2fcls }} instances - referred to as *chains* - from arbitrary *oframe/pframe* pairs of frames. Applications can query chain instances for the relative pose, spatial velocity and spatial acceleration across the oframe/pframe pair for the chain. The chain instances internally combine the contributions by the appropriate edges along the connecting path to compute the overall requested data across the path. The {py:class}`Karana.Frame.ChainedFrameToFrame` chains can be defined using arbitrary oframe/pframe {{ framecls }} pairs from the frames tree. The {py:class}`Karana.Frame.OrientedChainedFrameToFrame` chains can only be used when the oframe is an ancestor of the pframe {{ framecls }}. (create_frame_container_sec)= #### Creating a frame container We start off the creation of a frames layer by creating a {{ framectrcls }} by calling its {py:meth}`create() ` static method. This creates a {{ framectrcls }} with a single root {{ framecls }} that serves as the anchor for the rest of the frames in the system. ```{eval-rst} .. _framesimg: .. figure:: images/frame_layer.png :width: 600px :align: center :alt: frames layer The frames layer ``` Additional frames can be created by calling the {{ framecls }} {py:meth}`create() ` static method. Each new {{ framecls }} needs to be attached to a {{ framecls }} in the frames tree. This can be done by creating a {{ f2fcls }} edge by calling the {py:meth}`Karana.Frame.PrescribedFrameToFrame.create` method. This method takes a parent and unattached child {{ framecls }} arguments and attaches the child to the parent in the tree. The first (i.e. the "from") {{ framecls }} in a {{ f2fcls }} is referred to as the _oframe_ {{ framecls }}, and the second (i.e. the "to") {{ framecls }} is referred to as the _pframe_ {{ framecls }}. That `o` comes before `p` alphabetically mnemonic can be helpful to remember this naming convention. While these are the basic steps to build up the frame tree, other steps, such as adding bodies and hinges, also add frames to the frame tree. Once the frames tree has been built, the {{ framectrcls }} {py:meth}`ensureHealthy() ` method must be called to finalize the tree so that it is ready for use. A {{ framecls }} can be looked up by name using the {{ framectrcls }}'s {py:meth}`lookupFrame() ` method. Calling a {{ framecls }}'s {py:meth}`dumpFrameTree() ` methods displays the frame tree hierarchy starting with the {{ framecls }}. The following is example Python code for creating a {{ framectrcls }}, adding {{ framecls }} instances and attaching them. (create_frame_container_anchor)= ```python # create a frame container fc = FrameContainer.create() # get the root frame root_frame = fc.root() # create a new frame frmA = Frame.create("frameA", fc) # attach the new frame to the root frame f2fA = PrescribedFrameToFrame(root_frame, frmA) # create a second frame frmB = Frame.create("frameB", fc) # attach the B frame to the A frame f2fB = PrescribedFrameToFrame(frmA, frmB) # finalize the frame container once the structural changes are done fc.ensureHealthy() # display the current frame hierarchy root_frame.dumpFrameTree() ``` A frame can be detached from the frame tree by deleting its edge: ```python # get the edge e = frm.edge() # delete the edge e.discard() # verify that the edge is gone assert not frm.edge() ``` A {{ framecls }} can be reattached by creating a new {py:class}`Karana.Frame.PrescribedFrameToFrame` edge connecting it to another parent {{ framecls }}. Note that a {{ framectrcls }} cannot be made healthy(via its {py:meth}`ensureHealthy() ` method) if there are any _unattached_ frames in the system. #### Edge relative data The application is only responsible for setting the edge relative transform spatial velocity and spatial acceleration values. The relative homogeneous transform for an edge {{ f2fcls }} can be set using its {py:meth}`setRelTransform() ` method, the relative spatial velocity using its {py:meth}`setRelSpVel() ` method, and the relative spatial acceleration via its {py:meth}`setRelSpAccel() ` method. The following `Python` code illustrates how to set the edge data for the {ref}`frame creation example ` above: ```python # set the edge transform for frmA f2fA.setRelTransform(HomTran(UnitQuaternion(1, 0, 0, 1), [1,2,3])) # set the edge relative spatial velocity for frmA f2fA.setRelSpVel(SpatialVector([3, 4, 5], [1,2,3])) # set the edge relative spatial acceleration for frmA f2fA.setRelSpAccel(SpatialVector([4, 5, 1], [.1,0,.6])) ``` For some {{ f2fcls }} edges, the pose, velocity, acceleration data is fixed and needs to be set just once (e.g., for rigidly attached frames). More generally, callbacks can be registered with an {{ f2fcls }} edge to automatically update the transform and related data. This callback mechanism is used for auto-updating {{ f2fcls }} edge data for subhinges after coordinate changes, and changes to ephemeris time for Spice edge {{ f2fcls }} etc. #### Ephemerides frames {py:class}`Karana.Frame.SpiceFrame` class is a specialized {{ framecls }} class for frames associated with celestial and other bodies for ephemerides. For instance, these frames are associated with celestial bodies in the solar system such as the Sun, the Earth, the Moon, Mars etc. The {py:class}`Karana.Frame.SpiceFrameToFrame` class is for edges defined by {py:class}`SpiceFrame ` oframe/pframe pairs. The motion of a Spice {{ framecls }} can queried from ephemeris kernels from the [JPL NAIF toolkit](https://naif.jpl.nasa.gov/naif/). The {py:meth}`Karana.Frame.SpiceFrame.lookupOrCreate` method can be used to create such a frame using `NAIF` body and frame ids. The {py:class}`Karana.Frame.SpiceFrameToFrame` edges are created from pairs of Spice {{ framescls }} instances. The data for the {py:class}`Karana.Frame.SpiceFrameToFrame` edges is computed for the specified time epoch from the loaded `NAIF` kernels. Once created, these Spice {{ framecls }} instances become a part of the frames tree, and thus can be used seamlessly with any other {{ framecls }} in the frames layer. Thus for example, querying the direction of the sun, or the moon from a camera or a sun sensor attached to a terrestrial platform simply requires the creation of the Spice {{ framecls }} for the relevant celestial bodies, after which relative transform queries from any other {{ framecls }} such as the camera can be made directly. The [Spice frames example notebook](generated/notebooks/example_spice_frames/notebook.ipynb) contains a celestarium example to illustrate creating {{ framectrcls }} and {{ framecls }} instances and the use of Spice frames. #### Vector time derivatives It is important to keep in mind that there are potentially four frames involved when computing time derivatives of vectorial quantities for velocities and accelerations. These have to carefully handled in the presence of rotating frames as is common place for multibody systems. These four frames are: - the first two are the _from_ and _to_ frames that are referred to as the _oframe/pframe_ frame pair. - the next frame is the _observing_ frame, i.e. the frame in which the oframe/pframe relative quantity whose derivative is being taken is represented. The observing frame often is the same as the oframe, but that is not a requirement. - the fourth frame is the _representation_ frame, i.e. the frame in which the derivative quantity (which is itself a vector) is represented. Often this representation frame is the same as the observing frame, but this is again not a requirement. {{ kdflex }} provides the {{ f2fcls }} {py:meth}`relTransform() `, {py:meth}`relSpVel() `, and {py:meth}`relSpAccel() ` methods to return the relative pose, spatial velocity and spatial acceleration respectively of the oframe/pframe pair {{ f2fcls }} as observed from and represented in the oframe {{ framecls }}. The {{ kdflex }} frame layer convention is that the oframe is the observing and representation {{ framecls }} for the quantities returned by these methods. There are additional methods available to transformation the data into alternative observing and derivative {{ framecls }} choices (see the {ref}`observer_frame_sec` section). Another noteworthy observation is that with rotating frames, there are Coriolis acceleration terms when working with accelerations, and their specific values depends on the specific choices of the oframe, pframe, observing and representation {{ framecls }}. #### Arbitrary frame pair relative data The frames layer supports on demand computation of relative pose, spatial velocity, and spatial acceleration for any pair of frames in the system. The frame pair may be on arbitrary branches, and the connecting path may involve traversing edges in the opposite direction. Combining the edge relative transform values along the path is taken care of by the frames layer for each of the transform, velocity and acceleration levels. ```{eval-rst} .. _cSpVel: .. figure:: images/spvels.png :width: 600px :align: center :alt: combining spatial velocities Combining spatial velocities across frames ``` This combination process takes into account the transformations needed to combine data from the different oframe, pframe, observing and representation frames encountered along the path. Calling a {{ framecls }}'s (the _from/oframe_ frame) {py:meth}`frameToFrame() ` method with a _to/pframe_ {{ framecls }} argument will return a {py:class}`Karana.Frame.ChainedFrameToFrame` {{ f2fcls }} chain instance with information about the path connecting its oframe/pframe pair. The {py:meth}`relTransform() ` method can be called on the {{ f2fcls }} chain to get the relative transform for the frame pair, the {py:meth}`relSpVel() ` method for the relative spatial velocity, and {py:meth}`relSpAccel() ` method for the relative spatial acceleration of the pframe with respect to the oframe, with the oframe also serving as the observing and representation {{ framecls }}. The following code shows how to create a chain and query the chain's relative pose etc data for the above {ref}`frame creation example `: ```python # get the root/frmB frame to frame chain f2f = root.frameToFrame(frameB) # get the relative transform for the chain T = f2f.relTransform() # get the relative spatial velocity for the chain V = f2f.relSpVel() # get the relative spatial acceleration for the chain A = f2f.relSpAccel() # get the oframe and pframe for the chain of = f2f.oframe() pf = f2f.pframe() ``` The user is also free to create chains such as `frameB.frameToFrame(root)` where the _pframe_ is an ancestor of the _oframe_, or even on a different branch altogether. It is important to appreciate that this chain's relative data can well be non-zero, because even though the _root_ {{ framecls }} is the anchor (and does not move), the _root_ {{ framecls }} will appear to move when observed from another moving {{ framecls }}. The {{ framecls }}'s {py:meth}`frameToFrame() ` method first checks whether the required {{ f2fcls }} instance already exists or needs to be created, and creates one if necessary. For applications that have a recurring need for a specific {{ f2fcls }}, it is a good practice to cache the {{ f2fcls }} instance locally to avoid the costs from repeated look ups. (observer_frame_sec)= #### Changing observer and derivative frame By convention, the {{ f2fcls }} {py:meth}`relSpVel() ` and {py:meth}`relSpAccel() ` methods compute quantities relative to the oframe and with the oframe being the observer {{ framecls }}. Additional methods are available for switching to an alternative observer {{ framecls }} for vector derivatives. The {{ f2fcls }} {py:meth}`pframeObservedRelSpVel() ` and {py:meth}`pframeObservedRelSpAccel() ` methods return the oframe/pframe relative spatial velocity and acceleration respectively with the pframe being the observer frame. The {{ f2fcls }} {py:meth}`toPframeRelativeSpVel() ` and {py:meth}`toPframeRelativeSpAccel() ` methods go further and are handy methods for transforming hypothetical oframe relative spatial velocity and acceleration derivatives into pframe relative ones. The {{ f2fcls }} {py:meth}`toOframeObserved() ` method can be used to transform arbitrary pframe observed 3-vector derivative into an oframe observed vector derivative. Thus assume that we have frame `A`, a 3-vector $X$ and its `A` observed time derivative `Xdot_A`. To convert this derivative into another frame `B` observed derivative. `Xdot_B`, run ```python Xdot_B = B.frameToFrame(A).toOframeObserved(X, Xdot_A) ``` This method is general can be used for both spatial velocity and acceleration quantities. #### Remarks on frames Frames are typically created to track locations of interest in a simulation where bodies are in motion, and subhinges are articulating. The motion of some of the frames may be driven by the dynamics of bodies being simulated, while in other cases the motion may be driven by sources such as ephemeris which defines the motion of celestial bodies over time. A frame layer that uniformly includes all of the entities of interest in a simulation allows users to seamlessly make queries about the relative motion of one frame with respect to another for frames from diverse and disparate sources. Examples of specializations of frame classes are bodies and body node classes. Subhinges are specializations of edge {{ f2fcls }} classes, while body hinges are specializations of. oriented, chained {{ f2fcls }} classes. It is thus trivial for example to query the relative transform or velocity of a {{ framecls }} on one vehicle with respect to another even when there are several moving frames in between. The frames layer relieves applications from having to implement custom code for making such {{ f2fcls }} queries. – bespoke implementations that can often be fragile and error prone. It is worth remembering that the structure of a frame tree is not unique. The only requirement is that the edge transform and related values be set correctly for the specific frame tree structure. In practice, there can be hundreds if not thousands of frames in a simulation. Keeping the frame to frame data updated for all frame pair combinations in the frame tree at all times is impractical. Therefore, the frames layer uses an on demand, lazy computation model that computes relative frame pair values on demand. There is also a built in data-caching capability that avoids recomputing frame to frame data that have not changed. - - - (multibodytop_sec)= ## Multibody Layer (physical_bodies_sec)= ### Physical bodies and hinges In {{ kdflex }}, a multibody system is an instance of the {{ mbody }} class. The multibody representation consists of a tree or graph of physical bodies connected by hinges that define their permissible relative motion. Before turning to the process for creating multibody instances in the {ref}`treembody_sec` section, we first describe the bodies and hinges that populate such multibody models. {{ kdflex }} supports a number of classes for different body types derived from the {py:class}`Karana.Dynamics.BodyBase` class. The different body types are: - the {{ physbody }} class for physical rigid bodies - the {{ modalbody }} class for small deformation physical flexible bodies. This class is derived from the {{ physbody }} class. - the {py:class}`Karana.Dynamics.CompoundBody` container class for a connected {{ subtreecls }} of {{ physbodycls }} instances treated as a compound body unit. Compound bodies are not {{ physbodiescls }} themselves. Note that the {{ physbodycls }} class is sub-classed from the {py:class}`Frame ` class, and thus {{ physbodycls }} instances are also {{ framescls }} with additional attributes such as mass properties. We will use the _physical body_ terminology for concepts that only apply to physical bodies, and just _body_ terminology for ones that even apply to compound bodies. (FA_TA_models_sec)= ### Multibody representations Within {{ kdflex }}, {{ physbodycls }} instances are connected in parent/child relationship via {py:class}`Karana.Dynamics.PhysicalHinge` hinge instances. The {{ physbodycls }} members of a {{ mbodycls }} system are organized in a tree structure, with each {{ physbodycls }} having at most a single parent {{ physbodycls }}. A multibody graph system's closed-chain topology is represented as the combination of a {{ physbodycls }} tree together with an additional set of cut-joints defined as {py:class}`Karana.Dynamics.LoopConstraintCutJoint` hinge-loop-constraint instances (see {ref}`cutjoint_constraints_sec` section). For regular serial and tree topology multibody systems, the set of constraints is empty. Such graph representations are not unique. As discussed in {ref}`cutjoint_constraints_sec` section, {{ hingecls }} and {{ cutjtconstraintcls }} represent complementary ways of characterizing the permissible motion between nodes on {{ physbodiescls }}, with a {{ hingecls }} providing an _explicit_ way, and {{ cutjtconstraintcls }} providing an alternative _implicit_ way. While a {{ hingecls }} instance can be replaced with an equivalent {{ cutjtconstraintcls }} instance, the converse is not always possible without violating the tree-topology structure requirement for {{ physbodiescls }} and hinges. Among the possible graph representations for a {{ mbodycls }} system, at one extreme is a representation where all {{ physbodycls }} are independent and connected to the root via a 6-dof hinge. Thus the {{ physbodycls }} instances do not have children, and the motion constraints on each {{ physbodycls }} is defined via a set of {{ cutjtconstraintcls }} instances. We refer to such a multibody representation as a _Fully Augmented (FA)_ model. FA models are also referred to as _absolute_ coordinate multibody models in the literature. At the other end of the spectrum is a multibody graph representation consisting of a maximal spanning tree of {{ physbodiescls }} connected in a tree-topology structure via a set of inter-body {{ hingecls }} instances, together with a minimal set of {{ cutjtconstraintcls }} for cut-joints. We refer to such a multibody representation as a _Tree Augmented (TA)_ model. Note that converting all the {{ hingecls }} instances in a TA model representation into {{ cutjtconstraintcls }} instances turns it into an FA model representation. Between the span of representations book-ended by the FA and TA representations lie a whole range of intermediate representations where only a subset of the {{ hingecls }} instances are converted into {{ cutjtconstraintcls }} instances. The number of degrees of freedom (dofs) in a {{ mbodycls }} system is defined by the difference between the number of generalized velocity coordinates (see {ref}`coords_sec` section) for the {{ mbodycls }} system, and the number of constraints on the system. The number of dofs is invariant for a {{ mbodycls }} system across the different representations. This does imply that the number of coordinates and the number of constraints rise and fall together across different multibody representations since their difference is fixed. {{ kdflex }} supports both FA and TA {{ mbodycls }} representations, and arbitrary intermediate versions where some of the {{ hingescls }} in a TA model are replaced with {{ cutjtconstraintcls }} instances. Within this set of multibody representations, the TA model requires the smallest number of coordinates to describe the system configuration (and thus also the smallest number of constraints). The computational cost of the {{ kdflex }} {{ SOA }} based algorithms typically increases with the number of coordinates, and thus it is advantageous, and the default, to use TA multibody representations to the degree possible. There are methods available to transform among multibody representations. In contrast with {{ kdflex }}, most general-purpose multibody dynamics tools favor the FA multibody representations because of its simplicity. However, this simplicity comes with a large computational cost. The {{ SOA }} methodology provides tools to manage the increased complexity of TA models, and take advantage of minimal coordinate representations to develop fast computational dynamics algorithms. {{ kdflex }} uses these {{ SOA }}-based algorithms and thus favors TA model representations by default. ```{eval-rst} .. _mbodyrep: .. figure:: images/FA_TA_CE.png :width: 600px :align: center :alt: multibody representations Representation options for a multibody system ``` While TA models having a smaller number of coordinates when compared with FA models, TA models are not minimal coordinates when there are loop and coordinate constraints. There is yet another multibody representation, referred to as _Constraint Embedding (CE)_ that is truly minimal. CE models are an advanced topic, that we defer documentation for later. (body_nodes_sec)= #### Body nodes The {py:class}`Karana.Dynamics.Node` class is used for *nodes* that represent locations of interest on a {{ physbodycls }}. Once a {{ physbodycls }} has been created, {{ nodecls }} node instances can be created on the {{ physbodycls }} using the {py:meth}`Karana.Dynamics.Node.lookupOrCreate` static method. This class is derived from the {py:class}`Karana.Frame.Frame` class and the {{ nodecls }} instances are thus also {{ framecls }} instances - albeit attached to a parent {{ physbodiescls }}. One can think of a {{ nodecls }} as the special type of a {{ framecls }} that is directly attached to a {{ physbodycls }}. Node may be used to attach other {{ physbodiescls }}, or represent attachment points for sensors and actuators attached to the parent {{ physbodycls }}. A {{ nodecls }}'s {py:meth}`setBodyToNodeTransform() ` method can be used to set the (undeformed) location and orientation of a node on its parent {{ physbodycls }}. Since each {{ nodecls }} is also a {{ framecls }}, it is a simple matter to query a {{ nodecls }}'s pose, velocity etc. information relative to any other {{ framecls }} using the frames layer (see {ref}`frames_layer_sec` section). Selected {{ nodescls }} can be designated as ones that can apply spatial forces on the body, and these {{ nodescls }} are referred to as *force nodes*. Nodes that serve as attachment points for actuator devices (e.g., thrusters) that apply forces on a {{ physbodycls }} must be designated as force nodes or else the applied forces will be ignored during dynamics computations. While the location of a {{ nodecls }} on its parent {{ physbodycls }} remains fixed on a rigid {{ physbodycls }} instance, the location and orientation of nodes attached to a {{ modalbodycls }} instance (see {ref}`modal_flex_body_sec` section for more on flexible bodies) can change as the body deforms. ```{eval-rst} .. _nodes: .. figure:: images/body_nodes.png :width: 150px :align: center :alt: body nodes Body nodes ``` Before carrying out any computations, the parameters for all of the {{ physbodiescls }} and and hinges need to be set. A {{ physbodycls }}'s parameters can be set via its {py:meth}`setParams() ` and other related methods. (body_hinge_sec)= #### Connecting bodies via hinges A {{ physbodycls }} is connected to its parent {{ physbody }} via a {{ hingecls }} of {py:class}`Karana.Dynamics.PhysicalHinge` type. A {{ hingecls }} is a specialization of the {py:class}`Karana.Frame.OrientedChainedFrameToFrame` class. ```{eval-rst} .. _hinges: .. figure:: images/body_hinge.png :width: 600px :align: center :alt: hinges Body pair connected by a hinge ``` The hinge creation process automatically creates a pair of {{ nodecls }} instances, referred to as an `onode `of type {py:class}`Karana.Dynamics.HingeOnode` and a `pnode` of type {py:class}`Karana.Dynamics.HingePnode`. The {{ onodecls }} is attached to the parent {{ physbodycls }} and the {{ pnodecls }} to the child {{ physbodycls }}. These nodes serve as the oframe/pframe for the {{ hingecls }} {{ f2fcls }} chain. The {{ onodecls }} and {{ pnodecls }} define the location of the {{ hingecls }} on the pair of connected {{ physbodiescls }}. Every {{ physbodycls }} instance has a unique {{ pnodecls }} that is the connection point for its parent {{ hingecls }}, and has one or more {{ onodecls }} instances where the {{ hingecls }} for each attached child {{ physbodycls }} is located. The location and orientation of the {{ onodecls }} on the parent {{ physbodycls }} can be set using the {{ onodecls }}'s {py:meth}`setBodyToNodeTransform() ` method. The location and orientation of the {{ pnodecls }} on the child {{ physbodycls }} can be set using the {py:meth}`setBodyToJointTransform() ` method for the child {{ physbodycls }}. ```{eval-rst} .. _hnodes: .. figure:: images/body2joint.png :width: 500px :align: center :alt: hinge nodes Hinge node locations ``` The articulation of the child {{ physbodycls }} with respect to the parent {{ physbodycls }} is defined entirely by the relative motion of the hinge's {{ pnodecls }} with respect to the {{ hingecls }}'s {{ onodecls }}. The {{ onodecls }} and {{ pnodecls }} frames for a hinge coincide with the hinge coordinates are zero. New {{ physbodiescls }} can be attached to existing {{ physbodiescls }} by creating a {{ hingecls }} of a specified type. The {py:meth}`Karana.Dynamics.PhysicalHinge.create` static method creates {{ hingecls }} instances. This method takes a {py:class}`Karana.Dynamics.HingeType` argument to specify the type of hinge to be created. The {py:class}`Karana.Dynamics.PhysicalHinge` is derived from the {py:class}`Karana.Dynamics.FramePairHinge` class which implements much of the functionality. The main difference between these classes is that the latter class simply connects a oframe/pframe {{ framecls }} pair. While the {{ hingecls }} instances are used to connect {{ physbodiescls }} in the {{ mbodycls }} spanning tree. A {py:class}`FramePairHinge ` can also be used for defining {{ loopconstraintcls }} instances (see {ref}`cutjoint_constraints_sec` section). (subhinges_sec)= #### Subhinges A {{ hingecls }} by itself is a container class for a sequence of zero or more {py:class}`physical subhinges ` instances. The {py:class}`PhysicalSubhinge ` class is derived from the {py:class}`Karana.Frame.EdgeFrameToFrame` class, and hence represents an {{ f2fcls }} edge in the {{ framecls }} tree. There are a pre-defined set of subhinge types selected via the {py:class}`Karana.Dynamics.SubhingeType` available for use. Each {py:class}`Karana.Dynamics.HingeType` enum type corresponds to a specific pre-defined sequence of subhinge types. The one exception is the {py:attr}`CUSTOM ` hinge type which allows the user to specify a custom the sequence of {{ subhingecls }} types. (coords_sec)= #### Generalized Q, U etc coordinates Every {{ subhingecls }} in the {{ mbodycls }} system has a set of _generalized coordinates_ that parameterize its motion. Additionally, deformable bodies such as {{ modalbodycls }} have generalized coordinates that parameterize its deformation. All objects that contribute such generalized coordinates, including subhinges and {{ physbodiescls }}, are derived from the abstract {py:class}`Karana.Dynamics.CoordBase` class. Each {{ coordbasecls }} instance has an inherent polarity (i.e. orientation) with its coordinates parameterizing the relative motion of some frame with respect to another. Thus the coordinates for a {{ subhingecls }} that is part of a {{ hingecls }} parameterize the motion of the outboard child {{ physbodycls }} with respect to the inboard parent {{ physbodycls }} (and not the other way around!). The {{ kdflex }} conventions for generalized coordinates are described below: - The {{ coordbasecls }} {py:meth}`nU() ` method returns the number of degrees of freedom (dofs) for the {{ coordbasecls }} instance. The {py:meth}`nQ() ` method returns the number of configuration coordinates for the {{ coordbasecls }} instance. These values are generally equal, but there are exceptions such as for the {py:attr}`SPHERICAL_QUAT ` {{ subhingecls }} type where they are not. - each {{ coordbasecls }} instance has a vector of _generalized configuration coordinates_, denoted **Q** (of size `nQ()`) for the coordinates that parameterize the {{ homtrancls }} relative pose across the {{ coordbasecls }} instance. The {{ coordbasecls }} {py:meth}`getQ() ` and {py:meth}`setQ() ` methods can be used to get and set the `Q` values for the {{ coordbasecls }} instance. - each {{ coordbasecls }} instance has a vector of _generalized velocity coordinates_, denoted **U** (of size `nU()`) that parameterize the relative spatial velocity across the {{ coordbasecls }} instance. The {{ coordbasecls }} {py:meth}`getU() ` and {py:meth}`setU() ` methods can be used to get and set the `U` values for the {{ coordbasecls }} instance. In many cases, the *U* coordinates are the time derivatives of the *Q* coordinates, but this is not true when using quasi-coordinates (such as for {py:attr}`SPHERICAL ` {{ subhingecls }} type), and cases where even their sizes may differ (such as for {py:attr}`SPHERICAL_QUAT ` {{ subhingescls }} type)! - each {{ coordbasecls }} instance has a vector of _generalized acceleration coordinates_, denoted **Udot** (of size `nU()`) for the coordinates that parameterize the relative spatial acceleration across the {{ coordbasecls }} instance. The {{ coordbasecls }} {py:meth}`getUdot() ` and {py:meth}`setUdot() ` methods can be used to get and set the `Udot` values for the {{ coordbasecls }} instance. The `Udot` values are the time derivatives of the `U` velocity coordinates. In addition, - each {{ coordbasecls }} instance has a vector of _generalized coordinate rates_, denoted **Qdot** (of size `nQ()`) for the time derivatives of the `Q` coordinates. The {{ coordbasecls }} {py:meth}`getQdot() ` method can be used to get the `Qdot` values for the {{ coordbasecls }} instance. The `Qdot` values are auto-computed and thus read-only. The `Qdot` values are often the same as `U`, except for a {{ subhingecls }} of {py:attr}`SPHERICAL ` and {py:attr}`SPHERICAL_QUAT ` types which make use of angular velocity quasi-velocities for the `U` velocity coordinates. - a vector of _generalized force_ values, denoted **T** (of size `nU()`) for the forces being applied at the {{ coordbasecls }} instance. The {{ coordbasecls }} {py:meth}`getT() ` and {py:meth}`setT() ` methods can be used to get and set the `T` values for a {{ coordbasecls }} instance. - - - (treembody_sec)= ### Creating a tree multibody system The {{ mbodycls }} system is the top-level container for all the {{ physbodiescls }} belonging to the system. Each {{ mbodycls }} instance automatically comes with a *virtual root* {{ physbodycls }} that serves as the root/anchor for the actual {{ physbodiescls }} in the system. {{ kdflex }} supports {{ mbodycls }} systems with serial-chain, tree and graph topologies. Unlike serial-chains, tree systems can have branches. The distinction between tree and graph topology systems is that graphs can contain {{ loopconstraintcls }} and {{ coordconstraintcls }} between {{ physbodiescls }}. ```{eval-rst} .. _mbodytop: .. figure:: images/mbody_topologies.png :width: 600px :align: center :alt: multibody topologies Types of multibody topologies ``` A {{ mbodycls }} system can be created procedurally as well from Python DataStruct definitions, YAML, JSON, HDF5 and URDF model files. Model file exporters to these formats are also available. We begin with the steps need to create tree topology {{ mbodycls }} system. (manual_creation_sec)= #### Manual creation A {{ mbodycls }} system is created by calling the {py:meth}`Karana.Dynamics.Multibody.create` static method. This method takes as an argument the Newtonian {{ framecls }} which is used as the inertial frame for dynamics computations. This method creates a {{ mbodycls }} instance with a single virtual root body. This root body serves as an anchor for the remaining {{ physbodiescls }} that are created subsequently. New {{ physbodiescls }} can be added to the multi body system by calling the `create()` method on the appropriate physical body class. Once the {{ physbodiescls }} have been created. The {py:meth}`Karana.Dynamics.Multibody.ensureHealthy` method should be called to finalize the setup of the {{ mbodycls }} system. The {py:meth}`Karana.Dynamics.SubTree.displayModel` can be used to display the list of {{ physbodiescls }} and hinges in a {{ mbodycls }} system. ```{eval-rst} .. _displaymdl: .. figure:: images/displayModel.png :width: 500px :align: center :alt: displayModel Example output from `SubTree.displayModel()` ``` The {py:meth}`Karana.Dynamics.Multibody.dumpTree` method displays the topology of the {{ mbodycls }} system including the loop and coordinate constraints among the {{ physbodiescls }}. ```{eval-rst} .. _dumptree: .. figure:: images/dumpTree.png :width: 600px :align: center :alt: dumptree Example output from `SubTree.dumpTree()` ``` Creating {{ physbodiescls }} and hinges automatically add a number of frames and edges to the frames tree. You can examine the current state of the frame tree by calling the {py:meth}`dumpFrameTree() ` method on the root {{ framecls }}. An interesting exercise is to identify the segments of the frame tree contributed by the individual {{ physbodiescls }}, {{ subhingecls }} and {{ hingecls }} instances. ```{eval-rst} .. _bdframes: .. figure:: images/body_frames.png :width: 500px :align: center :alt: body frames Frames associated with bodies and hinges ``` Before carrying out any computations, all of the subhinge and deformation coordinates need to be initialized (see {ref}`subhinges_sec` section). The {py:meth}`Karana.Core.allReady` method can be called to run an audit to verify that there are no unset parameters. It is a good practice to call this method after {{ mbodycls }} system creation and initialization. The [2 link pendulum (basic)](generated/notebooks/example_2_link_pendulum/notebook.ipynb) example script illustrates the the manual process for creating a {{ mbodycls }} system. (coord_data_sec)= ##### The CoordData container In a {{ mbodycls }} system is created, the collection of {{ subhingescls }} and {{ modalbodiescls }} together define the system dofs. The {py:class}`Karana.Dynamics.CoordData` container class can be used to manage collections of such {{ coordbasecls }} instances. For this purpose, the {{ mbodycls }} model has - a {{ coorddatacls }} instance for all of the {{ subhingecls }} instances, and it can be accessed via the {{ mbodycls }}'s {py:meth}`subhingeCoordData() ` method - a {{ coorddatacls }} instance for all the deformation coordinates for the {{ modalbodiescls }}, and it can be accessed via the {{ mbodycls }}'s {py:meth}`bodyCoordData() ` method - a {{ coorddatacls }} instance for all the {{ cutjtconstraintcls }} coordinates, and it can be accessed via the {{ mbodycls }}'s {py:meth}`cutjointCoordData() ` method The {{ coorddatacls }} class provides methods such as {py:meth}`nQ() `, {py:meth}`getQ() `, {py:meth}`setQ() `that parallel the methods for individual subhinges, except that they work for lists of {{ coordbasecls }} instances. Thus the **Q** vector is the combined vector of the values from each of the subhinges. (See the {ref}`coords_sec` section for more information on `Q`, `U` etc coordinates.) The {py:meth}`dumpState() ` is handy for displaying the coordinates data for a {{ coorddatacls }}'s elements. ```{eval-rst} .. _coorddataimg: .. figure:: images/dumpState.png :width: 600px :align: center :alt: coord data output Example output from `CoordData.dumpState()` ``` When getting or setting the **Q** etc coordinates for a {{ coordbasecls }}, it is possible to do this directly via the {{ coordbasecls }} instance. This can be also be done indirectly via a {{ coorddatacls }} instance that contains the {{ coordbasecls }} instance. To do the so, you will need to find the correct slot (via the {{ coorddatacls }}'s {py:meth}`coordOffsets() ` method) for the {{ coordbasecls }} in the overall coordinate vector and get/set these values. While technically doable. the direct approach is more convenient to work with for individual {{ coordbasecls }} instances. (multibody_procedural_sec)= #### Procedural creation the previous section described the process for manually creating {{ physbodiescls }} one at a time. This process provides complete control over creating the right body types, hinge types and setting up the parent/child body relationships. Convenience methods are however available to add serial-chain and tree topology systems of {{ physbodiescls }} to a {{ mbodycls }} system. in bulk. The {py:meth}`Karana.Dynamics.PhysicalBody.addSerialChain` and {py:meth}`Karana.Dynamics.PhysicalBody.addTree` methods will create a serial and tree systems respectively of uniform {{ physbodiescls }}. These methods are especially useful for automating the creation of large {{ mbodycls }} systems for testing and benchmarking performance. These methods can be called multiple times to build up the {{ physbodiescls }} in the {{ mbodycls }} system. While these methods create uniform {{ physbodiescls }}, it is possible to alter the properties of {{ physbodiescls }} individually as desired. The [n link pendulum (procedural)](generated/notebooks/example_n_link_pendulum/notebook.ipynb) example illustrates the procedural creation of a {{ mbodycls }} system. (importing_model_data_files_sec)= #### Importing model data files There is also support for creating {{ mbodycls }} systems from model data files. The following model formats are currently supported: - URDF format model files widely used by the robotics community - The JPL DARTS dynamics simulation tool Python model format files - The SOADyn_types.SubGraphDS DataStruct generated `yaml`, `json`, `HDF5` and `Python` format model files. There is also support for exporting model files into the URDF and SOADyn_types.SubGraphDS DataStruct compatible formats. The [WAM arm](generated/notebooks/example_urdf/notebook.ipynb) and [M2020 rover](generated/notebooks/example_m2020/notebook.ipynb) example notebooks illustrate process for creating a {{ mbodycls }} system from model files. ```{eval-rst} .. _mbodyimp: .. figure:: images/importexport.png :width: 500px :align: center :alt: mbody import export Supported multiple multibody model file import/export formats ``` - - - (modal_flex_body_sec)= ### Flexible body systems In addition to multibody dynamics with rigid {{ physbodiescls }}, {{ kdflex }} also supports dynamics with flexible bodies, where the rigid/flex dynamics coupling is rigorously handled using the low-cost {{ SOA }} dynamics algorithms. The **FEMBridge** toolkit provides interfaces for processing finite-element-model (FEM) data from structure analysis tools such as NASTRAN to generate assumed modes models for component deformable bodies. Modeling of flexible {{ mbodycls }} system dynamics accurately and fast, allows {{ kdflex }} to support the needs for loads and structures analysis, as well as for control-system design and analysis. It can provide accurate and fast performance needed for closed-loop time-domain simulation testbeds as well as generation of linearized state-space models for frequency-domain analysis. ```{eval-rst} .. _structcontrol: .. figure:: images/loads2gnc.png :width: 600px :align: center :alt: struct control Multibody dynamics bridge between structures and controls discipline analysis needs ``` The flexible bodies are assumed to undergo small deformation. The flexibility data for such bodies is typically developed using FEA tools such as NASTRAN. Since these FEA models are typically high-dimensional, an _assumed modes_ approach, together with _component mode synthesis (CMS)_ is typically used to develop reduced order models for the flexible bodies. ```{eval-rst} .. _flexbd: .. figure:: images/flexiblebody.png :width: 600px :align: center :alt: deformable Multibody systems with deformable bodies ``` #### Creating flexible bodies A flexible body can be created in {{ kdflex }} using the {py:meth}`Karana.Dynamics.PhysicalModalBody.create` method. The number of modes is an argument for this `create()` method. The {py:meth}`nU() ` method can be used to query the number of modes for a body. Note that the returned value from this method is zero for rigid {{ physbodiescls }}. The parameters - over and beyond those for a rigid body - for a {{ modalbodycls }}, consist of - a modal stiffness vector defined by the modal frequencies. This can be set using its {py:meth}`setStiffnessVector() ` method. - a modal damping vector. This can be set using its {py:meth}`setDampingVector() ` method. - nodal vectors of size `6xnU()` for every {{ nodecls }} on the {{ modalbodycls }} (including {{ pnodecls }} and {{ onodecls }} instances). These nodal matrices are defined from the mode shapes generated by the external model reduction process. Each {{ nodecls }} on a {{ modalbodycls }} has a non-null {py:class}`Karana.Dynamics.ModalNodeDeformationProvider` instance which can be accessed via the {py:meth}`deformationProvider() ` method. The Its {py:meth}`setNodalMatrix() ` method can be used to set the nodal matrix for the {{ nodecls }}. #### Deformation coordinates Physical bodies are {{ coordbasecls }} objects. The generalized coordinates associated with a {{ modalbodycls }} are its deformation coordinates. The number of deformation coordinates can be queried using the body's {py:meth}`nU() ` method. Deformable bodies have `nU()` size `Q`, `U`, `Udot` coordinates as described in the {ref}`coords_sec` section for {{ coordbasecls }} objects. The {py:meth}`getQ() ` etc methods can be used to get and set these coordinates. The {{ mbodycls }} instance has a dedicated {{ coorddatacls }} instance with the deformation coordinates for {{ modalbodiescls }} that can be accessed via the {py:meth}`CoordData() ` method. Virtually all the {{ kdflex }} algorithms, including kinematics, Jacobians, loop constraints, dynamics etc, work for rigid as well as flexible bodies. This allows {{ kdflex }} to support arbitrary graph topology systems with arbitrary mix of rigid and flexible bodies. This general capability thus provides accurate and fast rigid/flex {{ mbodycls }} dynamics support capitalizing on the fast {{ SOA }} methods. (mbody_config_changes_sec)= ### Multibody configuration changes {{ kdflex }} supports the run-time changes to the {{ mbodycls }} system topology. The {ref}`treembody_sec` section described the steps for creating a new body. These steps can be done at any time to add a new body. Removing a body is done by simply calling the `discard()` method on the body after all dependencies have been removed as discussed in the {ref}`create_discard_sec` section. {{ kdflex }} also supports changing {{ mbodycls }} topologies from detaching and reattaching bodies as follows. - Calling the {py:meth}`detach() ` method for a {{ physbodycls }} will detach it from its current parent {{ physbodycls }} and will delete their connecting {{ hingecls }}. It will also create a new {{ hingecls }} of the specified type to attach the body to the {{ mbodycls }} virtual root body via a {py:attr}`FULL6DOF ` hinge. The new {{ hingescls }} coordinates are initialized so as to preserve the inertial pose, spatial velocity and spatial acceleration of the body across the detachment action. - The {py:meth}`reattach() ` method can also be called with a new parent {{ physbodycls }} argument. This will carry out steps similar to the {py:meth}`detach() ` method, except that the new {{ hingecls }} of the specified type is created and used to connect the body to the new parent {{ physbodycls }}. If the hinge type is {py:attr}`FULL6DOF `, then the new hinge coordinates are initialized to preserve the inertial pose, spatial velocity and spatial acceleration of the body across the reattachment. The inertial pose, velocity and acceleration of the the body are automatically preserved only for a 6-dof {{ hingecls }} type. This is based on the fact that it is the only hinge type capable of ensuring continuity in the inertial pose, velocity and accelerations of the body across the reattachment. Furthermore, using 6-dof hinges is simpler since they do not require the additional data such as the location of the new {{ onodecls }} and {{ pnodecls }}, and additional axes etc parameters. What about reattaching with a {{ hingecls }} type different than a 6 dof hinge as may be needed in practice while preserving inertial pose, velocity and accelerations of the {{ physbodycls }}. The following describes the general process for changing the {{ hingecls }} type. - to avoid the non-physical discontinuities in {{ physbodycls }} inertial pose and velocities, it is necessary that the relative pose, spatial velocities and accelerations of the body with respect to its new parent be achievable by the new desired hinge type. In practice, this may require a transition period (e.g., nulling out relative velocities during a docking scenario) to meet this requirement. - once the admissible relative pose etc conditions have been met for the new {{ hingecls }} type, record the new parent relative pose etc quantities - they will be needed later to initialize the new {{ hingecls }}'s coordinates. The {py:meth}`discard() ` method on the existing {{ hingecls }} can be called to remove it, followed by a call to {py:meth}`Karana.Dynamics.PhysicalHinge.create` method to create the new {{ hingecls }} of the right type between the new parent and the {{ physbodycls }}. - now set the {{ onodecls }}, {{ pnodecls }} locations and the axes etc parameters as needed for the new {{ hingecls }}. These steps are similar to the ones needed when attaching a newly created {{ physbodycls }} via a {{ hingecls }} as described in the {ref}`manual_creation_sec` section. - the final step is to set the {{ hingecls }} coordinates to preserve the body's original inertial pose etc. This can be done by calling the {{ hingecls }} {py:meth}`fitQ() ` method to find the best `Q` coordinates for the desired relative transform. This method will return the residual transform error, which should be zero if the desired relative transform is compatible with the new hinge type. Similarly calls to the {py:meth}`fitU() ` method should be made to initialize the `U` coordinates for the {{ hingecls }}, and {py:meth}`fitUdot() ` method to initialize the `Udot` coordinates for the {{ hingecls }}. After making such changes, it is important to make the {{ mbodycls }} system current to update it for the new configuration. Due to the structure-based nature of the {{ SOA }} algorithms, very little beyond these changes to the topology is needed for all the {{ kdflex }} algorithms to continue working. Most of the algorithms are structure-based and will simply start following the new system topology! There is no other bookkeeping updates needed for the changes to connectivity that may otherwise be required in conventional approaches. The cost of making these changes is very modest, and so easy to accommodate in complex scenarios. - - - (constraints_sec)= ### Bilateral closure constraints Bilateral constraints can be imposed to limit the admissible motions of {{ physbodiescls }}. Such constraints are can be imposed directly on the body subhinge coordinates or on the relative spatial velocity of physical body nodes. The former are referred to as *coordinate constraints* (see {ref}`coordinate_constraints_sec` section), while the latter as *loop constraints*. The presence of constraints implies that the system has non-minimal coordinates. By default this requires constraint kinematics solvers (see {ref}`constraint_kinematics_sec` section) to compute system states that satisfy the loop motion constraints. Any meaningful use of such closed-chain topology systems requires that the system states be such that the loop motion constraints are satisfied! Also, different solution techniques for the system dynamics are required since the algebraic motion constraints make the system dynamics model into a DAE system. {{ SOA }}'s constraint embedding technique can be used to get back to an ODE formulation, but that is an advanced topic for later. (cutjoint_constraints_sec)= #### Cut-Joint constraints Hinges provide an _explicit_ way to characterize the permissible constrained between a pair of {{ physbodiescls }}. An alternative, kinematically equivalent - albeit _implicit_ - way, is to define the motion constraints via *cut-joint constraints*. The {{ cutjtconstraintcls }} class defines bilateral loop constraints between {{ physbodiescls }}. This class takes advantage of the duality between hinges and constraints - in that admissible motions can be defined explicitly using a {{ hingecls }}, or implicitly via a {{ cutjtconstraintcls }}. A {{ cutjtconstraintcls }} can be created using the {py:meth}`Karana.Dynamics.LoopConstraintCutJoint.create` method. A {{ cutjtconstraintcls }} is typically created between a a pair of {{ cnodecls }} nodes of type {py:class}`Karana.Dynamics.ConstraintNode` instances attached to a pair of {{ physbodiescls }} to constrain their relative motion. The same method also can be used to create a motion constraint between a {{ physbodycls }} and the environment. The method takes a {{ f2fcls }} instance as argument to specify the {{ framecls }} pair being constrained. At least one of the {{ framescls }} is required to be a {{ cnodecls }}. The method also takes a {py:attr}`HingeType ` argument to specify the allowed motion across the {{ cutjtconstraintcls }}. A {py:class}`FramePairHinge ` instance of the specified hinge type is used internally to restrict the admissible motion for the {{ cutjtconstraintcls }}. This process allows us to use rigidly locked constraints, as well as more general cases where the motion is only partially constrained. When the motion is only partially constrained, there are additional generalized coordinates for the hinge associated with the cut-joint constraint. The {{ cutjtconstraintcls }} {py:meth}`poseError() ` etc methods can be used to compute the constraint error for individual constraints. Cut-joint constraints are required for {{ mbodycls }} systems with graph topologies, i.e. ones with topological loops. A set of {{ cutjtconstraintscls }} for the cut-joints, together with the {{ physbodiescls }} {{ subtreecls }}, are used by {{ kdflex }} to define general graph topology {{ mbodycls }} systems. While another conceivable option to convert a graph into a tree might be to split bodies rather than hinge, we do not adopt this option. There is no reasonable way to split the mass and modal properties of flexible bodies across the fracture parts. The extreme option of splitting a body while apportioning all the body properties to one part, and making the other part a dummy mass-less rigid body is also problematic since some algorithms cannot accept mass-less leaf bodies. modes bases. For rigid bodies, we can simply spit about the body frame and assign half mass properties to each half. The [Slider-crank with loop constraints](generated/notebooks/example_slider_crank/notebook.ipynb) example illustrates the creation of {{ cutjtconstraintscls }} for a graph topology {{ mbodycls }} system. (convel_loop_constraints_sec)= #### Convel loop-constraints Unlike {{ cutjtconstraintscls }}, {{ convelconstraintscls }} instances do not have a {{ hingescls }}. A {{ convelconstraintcls }} requires projections of the relative spatial velocity of a pair of {{ physbodiescls }} along an axis to track each other. Furthermore, a {{ convelconstraintcls }} does not apply at the relative pose level, and only applies at at the velocity and acceleration levels. A {{ convelconstraintcls }} is an instance of the {py:class}`Karana.Dynamics.LoopConstraintConVel` class. Since a {{ convelconstraintcls }} does not have a {{ hingecls }}, there are no coordinates associated with them either. A {{ convelconstraintcls }} can be either on the relative angular or linear velocities, and the component axis can be set via its {py:meth}`setUnitAxis() ` method. Both the {{ cutjtconstraintcls }} and {{ convelconstraintcls }} classes derive from the {py:class}`LoopConstraintBase ` class. This class's {py:meth}`hasHinge() ` method can be used to distinguish between these two loop constraint types. (coordinate_constraints_sec)= #### Coordinate constraints While {{ cutjtconstraintcls }} and {{ convelconstraintcls }} instances restrict the admissible spatial velocity across {{ cnodecls }} and {{ framecls }} instances, a {{ coordconstraintcls }} directly imposes algebraic constraints between generalized coordinates to restrict their motion. Examples of coordinate constraints are {{ physbodiescls }} coupled via gears, with their respective generalized coordinates values being subject to a fixed scalar ratios. The {py:meth}`Karana.Dynamics.CoordinateConstraint.create` method can be used to create a {{ coordconstraintcls }} instance to impose a scale ratio between the `Q` coordinates of a pair of {{ subhingecls }} (with the same number of dofs). The scale ratio can be set via the {py:meth}`setScaleRatio() ` method. This {{ coordconstraintcls }} can be used for a pair of rotational {{ subhingescls }} (e.g., for gears) or the combination of rotational and translational {{ subhingescls }}, such as for rack and pinion mechanisms. - - - (subtrees_sec)= ### Sub-Trees While most applications will carry out computations on the full {{ mbodycls }} system, it is possible to restrict computations to a **connected** set of {{ physbodiescls }} in the {{ mbodycls }} system. For instance, for simulations with multiple vehicles, the user may at times want to limit some of the kinematics, statics or dynamics computations to individual vehicles. A similar story applies to multi-limbed robotic systems where some tasks may involve just a sub-set of the limbs. An additional pertinent example would be mobile manipulator systems where the computations can be narrowed down to just the manipulator for scenarios where the mobility degrees of freedom are not active. The {py:class}`Karana.Dynamics.SubTree` class is a container class for a sub-tree of {{ physbodiescls }} for this purpose. The {{ physbodiescls }} of in this container must for a topological tree, and thus the sub-tree must have a single virtual root body that is the ancestor of all the {{ physbodiescls }} in the sub-tree. The {py:class}`Karana.Dynamics.SubGraph` is a derived class that additionally can have a list of {{ cutjtconstraintcls }} and {{ coordconstraintcls }} instances defining the loop and coordinate constraints between its {{ physbodiescls }}. The {{ mbody }} class itself is derived from the {py:class}`SubGraph ` class. All additional {{ subtreecls }} and {{ subgraphcls }} instances are children of the parent {{ mbodycls }} system, and contain subsets of the {{ physbodiescls }} in the {{ mbodycls }} system. The {ref}`computational_algorithms_sec` section describes a number of available system-level algorithms for kinematics, statics and dynamics computations that work on {{ subtreecls }} and {{ subgraphcls }} instances. These algorithms ignore the {{ physbodiescls }}, {{ hingecls }} and {{ framescls }} external to them when carrying out their computations. This generally means that the external bodies and frames are treated as kinematically frozen, and the mass properties of inboard and outboard bodies are ignored. Since a {{ mbodycls }} is also a {{ subtreecls }} and a {{ subgraphcls }}, passing in the {{ mbodycls }} allows the use of the same algorithm methods for the full system as well. A new {{ subtreecls }} is created from an existing {{ subtreecls }} instance by calling the {py:meth}`Karana.Dynamics.SubTree.create` method. This method takes the `new_root`, `use_branches` and `stop_at` arguments to tailor the {{ physbodiescls }} parent sub-tree to be included in the new sub-tree. The following describes the role of these options in more detail: (subtree_args_anchor)= - To create a {{ subtreecls }} that includes all the descendants of a specific body, specify the body as the `new_root` argument and empty lists for the `use_branches` and `stop_at` arguments. Note that the the specified body will be the virtual root of the new {{ subtreecls }}, and hence technically not belong to it. To actually include this body, specify the body's parent body (or another strict ancestor) as the `new_root` argument. You can also specify the `parent_subgraph.virtualRoot()` as the argument value if the body happens to be a base body for the parent {{ subtreecls }}. - If the `use_branches` argument is an empty list, then all branches starting from the `new_root` are included in the new {{ subtreecls }}. If a non-empty list of bodies specified as the `use_branches` argument, the included bodies are limited to only those on branches containing one of these bodies. The `use_branches` argument can thus be used to exclude unneeded branches in the parent {{ subtreecls }}. If the `new_root` argument is null, then the new-subtree uses the common ancestor body for the `use_branches` body as the root of the new {{ subtreecls }}. Thus at least one of the `new_root` or `use_branches` arguments is required to be non-null. - When the `stop_at` argument is an empty list, all bodies in the branches specified by the `new_root` and `use_branches` arguments are included in the new {{ subtreecls }}. When the `stop_at` list is non-empty, bodies that are descendants of bodies in the `stop_at` list are excluded from the new {{ subtreecls }}. Bodies in the `stop_at` list are required to be strict descendants of the new root body - since otherwise this would lead to an empty {{ subtreecls }}. - A common situation is the need to create a {{ subtreecls }} consisting of the minimal spanning tree of a set of bodies. Such a {{ subtreecls }} can be created by specifying a null `new_root` argument, and the same list of bodies as the `use_branches` and the `stop_at` argument. Arbitrary levels of nesting are allowed for {{ subtreecls }}s, and a {{ subtreecls }} can have an arbitrary number of children {{ subtreescls }}. A {{ subtreecls }}'s {py:meth}`virtualRoot() ` method can be used to get its root body. Since the set of bodies can vary across {{ subtreecls }} instances, so can the children bodies for a body. The {py:meth}`childrenBodies() ` method can be used to query the children bodies for a body in a {{ subtreecls }}. Since {{ mbodycls }} instance is also a {{ subtreecls }}, calling {py:meth}`childrenBodies() ` method on the {{ mbodycls }} returns the full list of children {{ physbodiescls }} for {{ physbodycls }}. One handy feature of any {{ subtreecls }} instance is that it automatically has a dedicated {{ framecls }} (accessed via its {py:meth}`cmFrame() ` method) to track the center of mass location of the bodies in the sub-tree. As a consequence, a {{ mbodycls }} system also has such a CM tracking frame. Since this is a {{ framecls }} instance, it is a part of the frames tree, and thus the location of a {{ subtreecls }}'s CM location can easily be queried with respect to any other {{ framecls }} in the system. (subgraphs_sec)= #### Sub-Graphs with Constraints The {py:class}`Karana.Dynamics.SubGraph` class is derived from the {py:class}`Karana.Dynamics.SubTree` class, and adds support for bilateral {{ loopconstraintcls }} and {{ coordconstraintcls }} constraints (see {ref}`constraints_sec` section) among the {{ subgraphcls }} component {{ physbodiescls }}. A {{ subgraphcls }} instance can be created using the {py:meth}`Karana.Dynamics.SubGraph.create` method. In addition to the arguments for the {py:meth}`Karana.Dynamics.SubTree.create` method - for selecting the bodies to include - the {{ subgraphcls }}, this `create()` method takes lists of {{ loopconstraintcls }} and {{ coordconstraintcls }} instances from the parent {{ subgraphcls }} to include in the new one. For convenience, a new {{ subgraphcls }} instance automatically inherits and enables all the constraints from the parent {{ subgraphcls }} that are legal for the subset of bodies in the new {{ subgraphcls }}. Such inheritance can be disabled by setting the `inherit_constraints` argument for the `create` method to `false`. ```{eval-rst} .. _subgraphs: .. figure:: images/subgraphs.png :width: 300px :align: center :alt: subgraphs SubGraph illustration ``` All constraint instances can be enabled and disabled for a {{ subgraphcls }} at run-time by calling its {py:meth}`enableConstraint() ` and {py:meth}`disableConstraint() ` methods respectively. Only the enabled constraints are taken into account in computations involving the {{ subgraphcls }}. This provides a convenient way to handle scenarios where constraints are changing at run-time such as during robot manipulation and mobility activities. The {py:class}`Karana.Dynamics.Multibody` class is a specialization of the {py:class}`Karana.Dynamics.SubGraph` class. Multiple {{ subtreecls }} and {{ subgraphcls }} classes can be created with different combinations of {{ physbodiescls }} and constraints for a {{ mbodycls }} system. Several of the computational algorithms in {{ kdflex }} act on {{ subtreecls }} and {{ subgraphcls }} instances. So while these algorithms can be invoked on the {{ mbodycls }} instance, they can also be called on {{ subtreescls }} and {{ subgraphscls }}. This allows an application to restrict the execution of these algorithms to a narrower set of the {{ physbodiescls }} in the {{ mbodycls }} system (e.g., to just the manipulator arm on a mobile platform, or to a single vehicle in a {{ mbodycls }} system with multiple vehicles). - - - (jacobian_sec)= ### Jacobians In the domain of robotics, Jacobian matrices play an important role. In the most basic case, Jacobians relate joint velocities to spatial velocities of a robot end-effector. Jacobians are configuration dependent and need to be computed on the fly in the context of robot motion control since the coordination of multiple limbs is required. In the broader context of multi-limb robots, there are multiple end-effectors. The concept of a Jacobian extends naturally to such multiple end-effectors. Again, in the context of robotics, Jacobians are needed for inverse kinematics (IK) solvers. The IK problem is one of finding a set of hinge coordinates that will position one or more end-effectors in a desired pose. This typically requires numerical solvers that often utilize these Jacobians within an iterative numerical search process. {{ kdflex }} includes support for the efficient computation of general purpose Jacobians, as well as broader constraint kinematics solver capabilities. For this, {{ kdflex }} first generalizes the notion of a Jacobian. Conventionally Jacobians are used to characterize end-effector poses and velocities with respect to the ground or the inertial frame. A more general setting is where the Jacobian is defined for a pair of frames, where both frames may be in motion. In this more general setting, the Jacobian characterized the relative motion between these pair of frames as a function of the joint coordinates. This approach subsumes the conventional case where the Jacobian is with respect to the inertial ground frame. This broader framing allows the use of the Jacobians for the more general _constraint kinematics (CK)_ problem of solving for the system state that satisfies a set of loop constraints. ```{eval-rst} .. _jacob: .. figure:: images/jacobian.png :width: 250px :align: center :alt: Jacobians Jacobians can be relative as well for multiple nodes ``` The {py:meth}`Karana.Dynamics.FrameToFrameJacobianGenerator.create` method can be used to create a Jacobian generator for a {{ framecls }} pair defined by a single {{ f2fcls }} instance for coordinates defined by a {{ coorddatacls }} instance. This object's {py:meth}`jacobian() ` method computes the relative Jacobian for the frame pair. The object does the caches bookkeeping information for reuse across multiple Jacobian calls. By default, analytical techniques are used to compute the Jacobians. The {py:class}`Karana.Dynamics.MultiJacobianGenerator` class is a more general version of the generator that can compute the Jacobian for a list of {{ f2fcls }} instances. Methods such as the {py:meth}`Karana.Dynamics.SubTree.subhingeCoordData()` can be used for the The {{ coorddatacls }} arguments when seeking Jacobians for all the subhinge coordinates in a {{ subtreecls }}. - - - (constraint_kinematics_sec)= ### Constraint kinematics The {py:class}`Karana.Dynamics.ConstraintKinematicsSolver` class is a general purpose constraint kinematics solver for solving a {{ subgraphcls }}'s coordinates that satisfy its motion constraints. The simplest way to create an instance of the solver is by calling the {{ subgraphcls }}'s {py:meth}`cks() ` method. This solver will be setup to solve for the {{ subgraphcls }} state satisfying its enabled constraints. The solver's {py:meth}`solveQ() ` etc methods can be invoked to solve for the coordinates satisfying the constraints. The returned value is the residual error. The frozen coordinates are left unchanged during the solution process. A feature of the {{ kdflex }} constraint kinematics solver is that it allows for the _freezing_ of a subset of the {{ subgraphcls }} coordinates using the {py:meth}`freezeCoord() ` method. Such freezing can be useful for {{ subgraphscls }} with redundant coordinates, where the overall strategy may be to use additional policies to manage a subset of the coordinates, and shield them from change during the solution processes. The [Slider-crank with loop constraints](generated/notebooks/example_slider_crank/notebook.ipynb) example illustrates the creation of and use of constraint kinematics solvers. (inversekin_sec)= #### Inverse kinematics In its simplest form, the _inverse kinematics (IK)_ problem is one of solving for joint coordinates required to achieve a desired pose for a frame such as the end-effector. We recast the IK problem into a constraint kinematics problem by creating a loop constraint whose satisfaction would meet the IK requirements. We can then proceed to use the constraint kinematics solver's solve methods to obtain the IK solutions. This approach is quite general, and allows us to solve IK problems involving relative pose requirements across moving frames, having multiple such requirements, and accommodate cases where the relative pose requirements are partial - such as position only, or orientation only. Consider for instance a multi-limbed robot, where the scenario calls for using one limb end-effector to anchor the robot by holding on to a handle, while planting one end-effector on the floor, and using two other limbs to grasp and move a task object from one position to another. Motion planning for this scenario requires solving the inverse kinematics for the robot across the whole motion trajectory. The following steps can be used to carry out the coordinated motion for the robot using {{ kdflex }}: - Create a {{ cutjtconstraintcls }} constraint of hinge of type {py:attr}`Karana.Dynamics.HingeType.LOCKED` between the planted end-effector and the ground, and add and enable it in the {{ subgraphcls }}. - Create a {{ cutjtconstraintcls }} of type {py:attr}`Karana.Dynamics.HingeType.LOCKED` constraint of hinge between the end-effector holding the handle, and the handle, and add and enable it in the {{ subgraphcls }}. - Create a pair of {{ cutjtconstraintcls }} constraints of hinge type {py:attr}`Karana.Dynamics.HingeType.LOCKED` between each of the end-effectors holding the task object and the point of attachment on the task object, and add and enable it in the {{ subgraphcls }}. - Enable the IK constraints in the {{ mbodycls }} instance by calling the {py:meth}`enableConstraint() ` method for each of them. - Call the {{ mbodycls }} {py:meth}`cks() ` method to get the {{ ckscls }} instance for the {{ mbodycls }}. - Select the subset of independent coordinates to use to parameterize the motion. These can be external coordinates, such as the desired position of the task object, along with redundant hinge coordinates within the {{ mbodycls }}. Use the `cks` {py:meth}`freezeCoord() ` method to freeze the coordinates that we do not want the IK solver to change (i.e. the independent coordinates). - Now set the independent coordinates to their (frozen) coordinate values. Then call the {{ ckscls }} {py:meth}`solveQ() `, {py:meth}`solveU() ` and {py:meth}`solveUdot() ` methods as needed to solve for the `Q`, `U` and `Udot` coordinates satisfying the constraints. The return value from these methods is the residual error - which should be zero if the IK solution was successful. The solution coordinates will be stored at within each {{ subhingecls }}, and the values can be queried individually, or for the full {{ mbodycls }} by calling it's `getQ()` etc methods. An important reminder is to always do the `Q`, `U` and `Udot` solutions in order, since each solution depends on the previous one. If solving for IK coordinates across a path, this step can be repeated over over with new values for the independent coordinates to get the matching set of IK solutions. - Once done with the IK solutions, disable the IK constraints from the {{ mbodycls }} by calling the {py:meth}`disableConstraint() ` method for each of them. Also, unfreeze the independent coordinates by calling the the `cks` {py:meth}`unfreezeCoord() ` method. If the IK solutions are required multiple times during the simulation scenario, an option is to create a dedicated {{ subgraphcls }} instance for IK solutions with all the bodies in the system. The user simply needs to use the {{ subgraphcls }} instead of {{ mbodycls }} in the steps above for the IK solution. The advantage of this approach is that, the enabling of constraints, and freezing of coordinates just needs to be done once - or as changes happen, without the need to undo them after use. - - - (algorithms_sec)= ### Dynamics Computations While the classes described above can be used to define the {{ mbodycls }} system model, the {py:class}`Karana.Dynamics.Algorithms` class contains several static methods for system level computational algorithms. The algorithms described here are are atomic, work-horse methods that are available for use by higher level layers in an application. These methods carry out computations for a specific state of the system, and use the `Q` and `U` generalized coordinate values set in the {{ mbodycls }} system when a method is called. Thus the kinetic energy method computes the kinetic energy at the state set in the {{ mbodycls }} instance. Before getting to the algorithms themselves, we first pause to introduce the notion of `prescribed motion` that is important for the forward dynamics related algorithms. #### Hybrid Forward Dynamics and Prescribed Motion The conventional _inverse dynamics problem_ is one of using the equations of motion to compute the `T` generalized forces for an input set of `Udot` generalized accelerations. The conventional _forward dynamics problem_ is the converse, and involves using the equations of motion to compute the `Udot` generalized accelerations for an input set of `Udot` generalized forces. However, there are times when the situation is mixed, and the inputs are a combination of `T` generalized forces at some subhinges and the `Udot` generalized accelerations at the remaining ones, and we need to compute the complimentary unknown generalized forces and accelerations. For this situation, the subhinges with input `Udot` generalized accelerations are referred to as ones undergoing `prescribed motion` since their motion problem is pre-determined. The {{ subhingecls }} {py:meth}`setPrescribed() ` method can be used to set and unset its prescribed motion property. By default, {{ subhingescls }} are non-prescribed. {{ kdflex }} includes support for the mixed situation with a _hybrid forward dynamics algorithm_. This algorithms will use the `T` generalized force values of non-prescribed {{ subhingescls }} as input and the `Udot` generalized acceleration values for prescribed ones. The hybrid algorithm will compute the unknown `Udot` values for the non-prescribed {{ subhingescls }}, and the `T` values for the prescribed ones. When none of the {{ subhingescls }} are prescribed, we fall back to the conventional forward dynamics problem. On the other hand, when all {{ subhingescls }} are prescribed, the resulting computations are those for the inverse dynamics problem. The {{ kdflex }} forward dynamics algorithm is in fact the hybrid forward dynamics algorithm so it can handle the the conventional inverse and forward dynamics cases, as well as the hybrid in-between cases. For the pure inverse dynamics case, there are however faster alternative Newton-Euler methods available for carrying out the computations. A noteworthy fact is that the prescribed property for any {{ subhingecls }} can be changed at any time during a simulation, and the hybrid algorithm will seamlessly adapt to the changes without any overhead. The equal and opposite spatial force between {{ physbodiescls }} at the {{ hingescls }} is often of interest. This quantity corresponds to Lagrange multipliers that are solved for in absolute coordinate dynamics formulations such as for the Fully-Augmented (FA) models. The {{ kdflex }} ATBI forward dynamics algorithm does not need or compute these quantities by default. However, these inter-body forces are easily computed from the ATBI algorithm intermediate quantities, and this technique is used by the {{ onodecls }} {py:meth}`getInterBodyForceTreeFwdDyn() ` method to compute the inter-body force at a {{ hingecls }}'s {{ onodecls }}. #### Gravitational acceleration The {py:meth}`setGravAccel() ` method can be used to set the gravity acceleration value for a {{ physbodycls }} in the {{ mbodycls }} system. To apply a uniform gravity field for all {{ physbodiescls }}, the {{ subtreecls }} {py:meth}`setUniformGravAccel() ` helper method can be used to set a gravity acceleration value for all the {{ physbodiescls }} in the {{ subtreecls }}. The {py:meth}`setGravityGradient() ` method is also available to set gravity gradient values for a {{ physbodycls }}. #### Dynamics in Non-Inertial Frames Normally, dynamics and equations of motion are described in a Newtonian, i.e. inertial {{ framecls }}. Inertial frames are by definition non-rotating and non-accelerating frames. The root {{ framecls }} in the frame container is by default the inertial frame in {{ kdflex }}. Any downstream frame that is non-rotating and non-accelerating is also an inertial frame and can be used for the purposes of dynamics computations. The {{ mbodycls }} virtual root body is quite often attached to the root frame so as to be a inertial frame to satisfy the requirements for the {{ kdflex }} dynamics algorithms. When using Spice {{ framescls }}, with celestial bodies, it is common practice to treat the Solar System Barycenter (SSBC) as an inertial frame (clearly an idealization, but one commonly used). Celestial bodies (such as the Earth) are typically accelerating and rotating with respect to the SSBC. Often terrestrial, and other planet referenced scenarios, are described with respect to the `planet centered rotating (PCR)` frame attached to the planetary body. As a rotating frame, the PCR frame is non-inertial. The question is how to include the effect of the PCR rotation into the dynamics computations. The conventional approach is to create a `planet centered inertial (PCI)` frame that is a non-rotating frame, but collocated with the PCR frame, and use it as the inertial frame. The PCR frame is made a child of the PCI frame. The planet is folded into the {{ mbodycls }} tree by creating a body for it and attaching it to the PCR frame. The vehicle - whose dynamics we are really interested in - is made into a child of the PCR body. To avoid mixing the dynamics of this massively heavy body with the vehicle bodies, the planetary body's parent hinge is marked as undergoing prescribed motion. With this assumption, the mass properties of the planet body have no effect on the dynamics, but the PCR rotation does get taken into account correctly in the vehicle dynamics. While technically sound, this approach does require the awkward step of including {{ physbodiescls }} for the PCI and PCR frames into the {{ physbodiescls }} tree, and adding ways to keep it in sync with the ephemerides trajectory for the planetary body. {{ kdflex }} includes a simpler and more streamlined solution for this problem that avoids the creation and inclusion of {{ physbodiescls }} for the PCI and PCR frames. This approach attaches the virtual root body to the non-inertial PCR frame. As a consequence the virtual root body has non-zero rotation and accelerations while serving as the attachment root for the vehicle bodies. The standard dynamics algorithms in {{ kdflex }} have been extended to accommodate such non-inertial virtual root {{ mbodycls }} bodies, such that they continue to capture the effects from the PCR rotation in the dynamics computations accurately. (computational_algorithms_sec)= ### Computational Algorithms {{ kdflex }} provides a large collection of fast computation algorithms to use with the {{ mbodycls }}, {{ subgraphcls }} and {{ subtreecls }} instances. The available algorithms include ones for doing inverse kinematics, inverse and forward dynamics for tree and graph systems, computing the system spatial momentum, mass properties and kinetic energy etc. The operation of the algorithms is restricted to the mini-world defined by the subset of their bodies. The algorithms thus ignore the {{ physbodiescls }}, {{ hingecls }} and {{ framescls }} external to them when carrying out their computations. This generally means that the external bodies and frames are treated as kinematically frozen, and the mass properties of inboard and outboard bodies are ignored. Since a {{ mbodycls }} is also a {{ subtreecls }} and a {{ subgraphcls }}, the methods can be used for system level computations by passing in the {{ mbodycls }} as the argument to these methods. There is no restriction on creating multiple {{ subtreecls }} and {{ subgraphcls }} instances with overlapping sets of bodies. However certain algorithms can only be run on non-overlapping {{ subtreecls }} instances. These methods can be used after the {{ subtreecls }} {py:meth}`enableAlgorithmicUse() ` has been called. This method will fail if there is another algorithmically enabled {{ subtreecls }} with overlapping set of {{ physbodiescls }}. One consequence of this is that if the {{ mbodycls }} itself has been algorithmically enabled, then no other {{ subtreecls }} can be enabled. The {{ subtreecls }}'s {py:meth}`disableAlgorithmicUse() ` method can be used to disable an enabled {{ subtreecls }}. Unless otherwise noted in the method documentation, all of the algorithm methods support rigid and flexible bodies. The available algorithms are described in the following sections. (kinematics_algorithms_sec)= #### Kinematics Algorithms The methods described here are for the system kinematics. - {py:meth}`Karana.Dynamics.Algorithms.jacobianGenerator`: Return a Jacobian generator for a {{ subtreecls }} system. This generator can be used to compute the Jacobian matrix (see {ref}`jacobian_sec` section) for a list of frames for a specified set of dofs defined by a list of {{ coorddatacls }} instances. - {py:meth}`Karana.Dynamics.Algorithms.constraintKinematicsSolver`: Return a constraint/inverse kinematics solver instance for a {{ subgraphcls }} system. Methods described in {ref}`constraint_kinematics_sec` section can be used to solve for generalized coordinates satisfying the constraints at the configuration, velocity and acceleration levels. (mass_prop_algorithms_sec)= #### System Properties Algorithms The following methods return quantities that are relevant for the analysis and design of systems. - {py:meth}`Karana.Dynamics.Algorithms.evalKineticEnergy`: Calculate the overall kinetic energy of the component bodies in a {{ subtreecls }} - {py:meth}`Karana.Dynamics.Algorithms.evalExternalSpatialForce`: Calculate the overall external spatial force on the tree from accumulating the spatial forces from all the external force nodes, active contact nodes and constraint forces from constraint nodes for the component {{ physbodiescls }} in a {{ subtreecls }} system. The returned value is at, and represented in, the {{ subtreecls }} center of mass frame. - {py:meth}`Karana.Dynamics.Algorithms.evalSpatialInertia`: Calculate the overall tree system spatial inertia from accumulating the spatial inertias for the component {{ physbodiescls }} in a {{ subtreecls }} system. - {py:meth}`Karana.Dynamics.Algorithms.evalTreeMassMatrixCRB`: Calculate the mass matrix of a {{ subtreecls }} using composite rigid body (CRB) inertia algorithm. An optional {{ coorddatacls }} can be used to customize the row/column order of the mass matrix. - {py:meth}`Karana.Dynamics.Algorithms.evalTreeMassMatrixInvDyn`: Calculate the mass matrix of a {{ subtreecls }} system using multiple invocations of the Newton-Euler inverse dynamics algorithm for computing the individual columns. An optional {{ coorddatacls }} can be used to customize the row/column order of the mass matrix. The mass matrix computed by this method is the same as the {py:meth}`Karana.Dynamics.Algorithms.evalTreeMassMatrixCRB` method, but the latter is computationally faster. - {py:meth}`Karana.Dynamics.Algorithms.evalTreeMassMatrixInvFwdDyn`: Calculate the mass matrix inverse of a {{ subtreecls }} system using forward dynamics for computing the individual columns. An optional {{ coorddatacls }} can be used to customize the row/column order of the mass matrix inverse. - {py:meth}`Karana.Dynamics.Algorithms.evalSerialChainMassMatrixInverse`: Calculate the mass matrix inverse of a serial-chain rigid-body system using {{ SOA }} based operator decomposition expressions. This method does not require the computation of the mass matrix or its inversion. An optional {{ coorddatacls }} can be used to customize the row/column order of the mass matrix inverse. - {py:meth}`Karana.Dynamics.Algorithms.evalTreeMassMatrixInverse`: Calculate the mass matrix inverse of a tree-topology rigid body system using {{ SOA }} based operator decomposition expressions. This method does not require the computation of the mass matrix or its inversion. While this method can also handle serial-chain systems, the specialized {py:meth}`Karana.Dynamics.Algorithms.evalSerialChainMassMatrixInverse()` alternative method is also available since it is much simpler for the simpler topology. An optional {{ coorddatacls }} can be used to customize the row/column order of the mass matrix inverse. - {py:meth}`Karana.Dynamics.Algorithms.evalFramesOSCM`: Calculate the `operational space compliance matrix (OSCM)` for a {{ subtreecls }} system. The `OSCM` is the reflected mass matrix inverse in the operational/task space as defined by a set of input frames. The original definition centered around the single end-effector frame for a robot arm, but the more general multiple frame version plays an important role in whole-body motion control applications for multi-limb robots. The formal definition of the OSCM is $J{\cal M}^{-1}J^*$, where $J$ denotes the Jacobian matrix for the frames, and $\cal M $ the system mass matrix. The direct computation of the OSCM by the direct evaluation of this expression requires the computation of the mass matrix and its inverse and can be computationally expensive. {{ kdflex }} instead uses a recursive {{ SOA }} algorithm that exploits the structure of the OSCM to develop a decomposition and computational algorithm that avoids even the need for computing the mass matrix. This is an important benefit since variants of the OSCM also play an important role in solving the forward dynamics for systems with constraints. (constraints_kin_algorithms_sec)= #### Constraints Related Algorithms A {{ subgraphcls }} with constraints will typically have generalized velocity coordinates that dependent on one another, as well as possibly redundant constraints. The following methods are for the computation of various constrained kinematics related properties for such {{ subgraphscls }}: - {py:meth}`Karana.Dynamics.Algorithms.evalVelocityConstraintMatrix`: Compute the velocity constraint matrix, $G_c$. The generalized velocity coordinates vector $U$ satisfies the loop and coordinate constraints on the system if $G_c * U = 0$, i.e. $U$ belongs to the null space of $G_c$. Note that the $G_c$ is configuration dependent in general. There are two forms of the $G_c$ velocity constraint matrix. When the `with_constraints` argument for the method is `true`, the number of columns in the computed matrix is the same as the total number of generalized velocity coordinates for the {{ subgraphcls }}, i.e. the value returned by its {py:meth}`nU() ` method, and $U$ contains all of the generalized velocity coordinates from the {{ subhingescls }}, {{ cutjtconstraintscls }} and {{ modalbodiescls }} velocity coordinates. The rows of the matrix correspond to the constraints on the {{ subgraphcls }}, with 6 rows for each {{ cutjtconstraintcls }}, and 1 row for each {{ convelconstraintcls }} and {{ coordconstraintcls }}. When the `with_constraints` argument for the method is `false`, the columns for the {{ cutjtconstraintscls }} coordinates are omitted from the matrix. The number of rows for each {{ cutjtconstraintcls }} in the returned matrix is equal to the number of motion constraints for the constraint and generally less than 6. Moreover, the row values also are different and map to the residuals for the constraint. Overall, the number of rows is fewer by the number of {{ cutjtconstraintcls }} velocity coordinates for the {{ subgraphcls }}. Nevertheless the $G_c * U = 0$ expression still holds, with $U$ being the smaller vector without the {{ cutjtconstraintcls }} velocity coordinates. - {py:meth}`Karana.Dynamics.Algorithms.evalIndepVelConstraintCoordIndices`: For a {{ subgraphcls }} with constraints, the number of independent generalized velocity coordinates $U$ is smaller than the overall number of {{ subgraphcls }} generalized velocity coordinates (as computed by its {py:meth}`nU() ` method). Moreover some of the constraints on the system may also be redundant. This method computes and returns a pair of lists of indices - the first being the list of indices of the best choice of _independent_ generalized velocity coordinates, denoted $U_i$, and the second being the list of indices of the constraints that are non-redundant. Note, that the results of this method are configuration dependent. Moreover, while this method returns the best choice of independent coordinates for a specific configuration, there are other possible choices for the independent coordinates. Also available is the analogous {py:meth}`Karana.Dynamics.Algorithms.evalIndepPoseConstraintIndices` and {py:meth}`Karana.Dynamics.Algorithms.evalIndepPoseCoordIndices` method for pose level constraints (useful when the pose and velocity constraints are not the same). - {py:meth}`Karana.Dynamics.Algorithms.evalIndepPoseCoordIndices`: For a {{ subgraphcls }} with constraints, the number of independent generalized coordinates $Q$ is smaller than the overall number of {{ subgraphcls }} generalized coordinates (as computed by its {py:meth}`nQ() ` method). This method computes and returns the list of indices of the best choice of _independent_ generalized coordinates, denoted $Q_i$. Note, that the results of this method are configuration dependent. Moreover, while this method returns the best choice of independent coordinates for a specific configuration, there are other possible choices for the independent coordinates. Also available is the analogous {py:meth}`Karana.Dynamics.Algorithms.evalIndepVelConstraintCoordIndices` method for velocity level constraints (useful when the pose and velocity constraints are not the same). - {py:meth}`Karana.Dynamics.Algorithms.evalIndepPoseConstraintIndices`: For a {{ subgraphcls }} with constraints, some of the constraints on the system may be redundant. This method computes and returns the list of indices of the constraints that are non-redundant. Note, that the results of this method are configuration dependent. Also available is the analogous {py:meth}`Karana.Dynamics.Algorithms.evalIndepVelConstraintCoordIndices` method for velocity level constraints (useful when the pose and velocity constraints are not the same). - {py:meth}`Karana.Dynamics.Algorithms.evalVelCoordinatePartitioning`: Since not all of the $U$ generalized velocity coordinates are independent, for constrained systems, it is possible to partition them into sets $U_i$ and $U_d$ independent and dependent coordinates respectively. This method returns a list of indices for the best choice of independent coordinates, along with the $D$ dependency matrix for computing the $U_d$ dependent coordinates from the independent coordinates using the $U_d = D\ U_i$ expression. - {py:meth}`Karana.Dynamics.Algorithms.evalDependentVelCoordinatesMatrix`: This method computes the $D$ dependency matrix for a specific choice (not necessarily the optimal one) of independent coordinate indices. When the optional `full` argument is `false`, the returned matrix is such that $U_d = D\ U_i$. However, when this argument is `true`, the returned matrix is such that $U = D\ U_i$, i.e. it maps to the full $U$ generalized velocity coordinates that include both the independent and dependent coordinates. Once again, the results of this method are configuration dependent. - {py:meth}`Karana.Dynamics.Algorithms.evalSqueezeMatrix`: For a {{ subgraphcls }} with constraints, there are generalized forces that can be applied that only effect internal forces and cause no change in motion. Such forces form a subspace whose dimension is the same as the size of the motion constraints on the system. These forces are referred to as _squeeze_ forces. This method 0 computes the $S$ squeeze force matrix such that $T_s = S\ x$ for arbitrary vector $x$. The row dimension of $S$ is {py:meth}`nU() `, while the column dimension is the number of constraints on the system. Note that this matrix may not have full rank. A corollary fact is that the squeeze force space is null when there are no constraints, because all applied generalized forces result in motion for such systems. Note that the above methods work on {{ subgraphcls }} instances. This allows the user to narrow down the constraint analysis to fragments of the overall multibody system, such as individual suspension system. a single robotic arm, just a hand as needed by simply creating and using {{ subgraphcls }} for these fragments and using the methods on them. This staged analysis capability can be very useful when working with large and complex full vehicle {{ mbodycls }} systems. (simulation_algorithms_sec)= #### Forward Dynamics Algorithms The following methods are the primarily for solving the equations of motion for simulation use. - {py:meth}`Karana.Dynamics.Algorithms.evalForwardDynamics`: This calls one of two algorithms under the hood depending on the argument passed. a. If the argument passed is a {{ subtreecls }} or a {{ subgraphcls }} without constraints, then a tree forward dynamics algorithm is used. This evaluates the articulated body inertia (ATBI) based forward dynamics. This algorithm actually implements _hybrid_ dynamics where some of the subhinges may be undergoing prescribed motion, and hence the inputs can be a combination of applied generalized forces at some subhinges, and generalized accelerations at the others. The algorithm is based on the {{ SOA }} `articulated body inertia (ATBI)` algorithm for tree systems. The algorithm has `O(N)` computational complexity, that is the cost increases linear with the number of {{ physbodiescls }} in the system. b. If the argument passed is a {{ subgraphscls }} with constraints, then method uses the {{ SOA }} based `Tree-Augmented (TA) algorithm` which involves two tree forward dynamics recursive sweeps and an additional step to compute constraint forces using the OSCM. The `Udot` generalized accelerations method computed by the method satisfy the constraints. The method is also `hybrid` in handling prescribed motion subhinges. Also, the algorithm supports the case where the hinge cut-joint constraints are actuated and thus have non-zero `T` generalized forces. The `Baumgarte` technique for constraint error stabilization can be enabled by setting the related parameters for the {{ mmcls }} instance. (embedded_control_algorithms_sec)= #### Embedded Control Algorithms The methods described here are for model based computations for use within embedded robotics control stacks. - {py:meth}`Karana.Dynamics.Algorithms.evalCmLocation`: Returns the current location of a {{ subtreecls }} center-of-mass with respect to the {{ subtreecls }}'s virtual root. This is a simple wrapper that makes use of the CM frame for the subtree to obtain its location. - {py:meth}`Karana.Dynamics.Algorithms.evalSpatialMomentum`: Calculate the overall spatial momentum for a {{ subtreecls }}. This method combines the spatial momentum contributions of the individual {{ physbodiescls }} in the subtree about a common reference frame. - {py:meth}`Karana.Dynamics.Algorithms.evalCentroidalMomentum`: Calculate the centroidal spatial momentum for a {{ subtreecls }}. This method computes the overall system spatial momentum about the system center of mass. This quantity is referred to as the `centroidal momentum` and is used in the dynamics and control of legged robotic platforms. - {py:meth}`Karana.Dynamics.Algorithms.evalCentroidalMomentumMatrix`: Calculate the centroidal momentum matrix for a {{ subtreecls }}. This method returns the matrix that maps the `U` generalized velocity coordinates to the system `centroidal momentum`. It is used in the dynamics and control of legged robotic platforms. - {py:meth}`Karana.Dynamics.Algorithms.evalGravityCompensation`: Compute the gravity compensation generalized forces for a {{ subtreecls }} system. This method computes the `T` generalized forces as a function of the system configuration that are needed to compensate for the gravitational load on the system. Ideally, these generalized forces will neutralize the effect of gravity and keep the system in a passive, steady state. - {py:meth}`Karana.Dynamics.Algorithms.evalInverseDynamics`: Evaluate the inverse dynamics for a {{ subtreecls }}. This method implements the fast Newton-Euler (NE) algorithm for the inverse dynamics computation. - {py:meth}`Karana.Dynamics.Algorithms.evalComputedTorque`: Calculate the computed torque using the Newton-Euler algorithm for a {{ subtreecls }}. This method uses the Newton-Euler inverse dynamics algorithm to compute the `T` generalized forces needed for a desired `Udot` generalized accelerations at the hinges. This method is used to compute the feed-forward motor torques for executing a motion profile. (gnc_algorithms_sec)= #### Guidance and Control Related Algorithms - {py:meth}`Karana.Dynamics.Algorithms.stateSpaceGenerator`: Create a generator that can create a linearized state space representation for a {{ mmcls }}'s {{ subtreecls }}. This method linearizes the system dynamics about the current state set point. It takes callback arguments for computing the system state derivative and output. The linearized system can be used for frequency domain control analysis. - {py:meth}`Karana.Dynamics.Algorithms.modalAnalysis`: Perform modal analysis for a {{ subtreecls }} system with deformable bodies (see {ref}`modal_flex_body_sec` section). This method creates a state space model of the flexible body system about the current configuration, and then solves an eigen-problem to compute the frequencies and mode shapes. Virtually all of these algorithms are recursive, structure-based and low-cost based on the {{ SOA }} framework. Most of these methods take a {{ subtreecls }} or {{ subgraphcls }} arguments, and hence their scope can be limited to a {{ subtreecls }} or {{ subgraphcls }} of the overall {{ mbodycls }} system if desired. We can of course use these methods on the {{ mbodycls }} instance itself (for computations on the full system) since it is itself a {{ subtreecls }} and a {{ subgraphcls }}. For example, since the {py:meth}`evalSpatialInertia() ` computes the combined spatial inertia of all the {{ physbodiescls }} in a {{ subtreecls }}, passing in the {{ mbodycls }} argument will compute the system level spatial inertia. Passing in a strict {{ subtreecls }} argument will compute the spatial inertia for just the {{ physbodiescls }} in the {{ subtreecls }}. A notable feature of {{ kdflex }} is the $O(N)$ nature of the forward dynamics algorithms - that is their computational cost scales only linearly with the number of bodies. This is in contrast with the $O(N^3)$ (cubic) computational complexity of the most conventional methods. The [Rigid body O(N) dynamics performance](generated/notebooks/example_rigid_timing/notebook.ipynb) example notebook verifies the O(N) performance of the forward dynamics algorithm by measuring computational time for a system with a varying number of bodies. The [Rigid body dynamics benchmark](generated/notebooks/example_fully_augmented_benchmark/notebook.ipynb) example compares the performance of the {{ kdflex }} forward dynamics algorithms against conventional method. Furthermore, the use of minimal coordinates avoids the need for constraint error management and allows the use of the simple ODE numerical integration techniques. - - - (simulationtop_sec)= ## Simulation Layer (system_level_models_sec)= ### System level models When modeling a physical system, the {{ mbodycls }} dynamics model is a critical, but one contributor to the overall physics of the system. A more complete model of the system requires integration of models for actuator and sensor devices attached to the system. and environmental interaction (e.g., gravity) effecting the system dynamics. {{ kdflex }} has a component based, object-oriented architecture for implementing and integrating such device and environment interaction models with the {{ mbodycls }} dynamics. ```{eval-rst} .. _clsd: .. figure:: images/closedloop.png :width: 350px :align: center :alt: closedloop Integration of multibody dynamics with component models in closed-loop simulations ``` {{ kdflex }} supports the implementation of a library of parameterized models, that can be instanced and reused across different simulations. The {py:class}`Karana.Dynamics.StatePropagator` and {py:class}`Karana.Dynamics.ModelManager` classes supports the development of system level simulations that combine the multibody dynamics with device and environment component models implemented using the {py:class}`Karana.Models.KModel` base class. The {{ kmodelcls }} classes support custom parameter classes and methods for interfacing with the multibody dynamics. (model_manager_sec)= #### The model manager A {{ mmcls }} instance can be used as the simulation manager for integrating the multibody dynamics with instances of component {{ kmodelcls }} models and other user defined classes to build up system level simulations. Such component models can be registered with the model manager for invocation and execution within the system level simulation. The {{ mmcls }} class defines the {py:attr}`MMFunctions ` struct with registries for methods of different types. The methods in these registries are invoked at different times in the system level execution. Thus the {py:attr}`pre_deriv_fns ` and the {py:attr}`post_deriv_fns ` methods are called just before and after respectively computing the {{ mbodycls }} state time derivative. Thus the `pre-deriv` methods can be used to set forces etc. that would effect the equations of motion, while the `post-deriv` methods can be used to query information (e.g., inter-body forces) that depend on the solution to the equations of motion. The {py:attr}`pre_hop_fns ` methods are called before an integration step is taken. These methods can be used to advance any discrete states within the component models. Other callbacks - including for zero-crossing condition - are available to further tailor model manager execution. While users can add methods to these registries directly, a more common (and recommended) way is to create a {{ kmodelcls }} component model subclass (see {ref}`kmodels_sec`) that has slots for compliant methods for implementation (see {ref}`kmodel_methods_sec` section). When a model is registered with the model manager, its methods are automatically registered with the model manager as well. The model manager will thus invoke the model's methods at the correct times and in the correct order during state propagation. This model manager design has several features that support the development of complex system level simulations with coupled multibody, device and environment models. Notable features are: - support for specialized {py:class}`Karana.Models.KModelParams` parameter classes to go hand-in-hand with component models sub-classed from the {py:class}`Karana.Models.KModel` class. - support for multi-rate models - support for tailoring the sequencing of method calls for the component models and the multibody dynamics (e.g., actuator force computation methods should be called before the multibody dynamics is evaluated) - ability to create multiple instances of a component model (e.g., one each for a thruster in a bank of thrusters), while allowing custom parameters for each instance. - creating a catalog of such component model classes for reuse across simulations. - run-time activation and deactivation of the component model instances - support for zero-crossing detection to terminate hops exactly when a zero-crossing condition is met. The {py:class}`MMFunctions ` with the callback registries are stored on the {py:attr}`fns ` member of the {{ mmcls }}. The creation and use of the model manager and component models can be found in several of the example scripts including - [2 link pendulum (basic)](generated/notebooks/example_2_link_pendulum/notebook.ipynb), [n link pendulum (procedural)](generated/notebooks/example_n_link_pendulum/notebook.ipynb) and [Slider-crank with loop constraints](generated/notebooks/example_slider_crank/notebook.ipynb). (system_state_sec)= ##### System state When using an integrator, the integrator takes ownership of the system state. This is necessary to accommodate the needs of adaptive step integrators which often maintain internal histories of the system state to manage error while propagating the system state. In this perspective, the {{ mbodycls }} system objects and component models serve as compute engines for computing the state derivative for the integrator to meet its state propagation needs. Simply changing the {{ mbodycls }} system coordinates is not therefore the way to update the simulation state. Such a new state has to be set by calling the {{ mmcls }} {py:meth}`setState() ` method - which also resets the internal state of the integrator. This method is used in the beginning to initialize the system state, and can also be used mid-stream should the need arise. If {py:meth}`setState() ` method is called mid-stream to set a new simulation time, the {{ mmcls }} will attempt to reschedule all {{ teventscls }} appropriately (see {ref}`timed_events_sec` section). This means, it will attempt to find the next time after the new time that a {{ teventcls }} should be executed. However, since {{ teventscls }} can have non-uniform execution times via the register_fn, it is impossible to capture all corner cases. Hence, this should be used with caution when modifying time mid-stream with non-uniform {{ teventscls }}. As a reminder, the system state vector is the concatenation of the `Q` generalized configuration coordinate values and the `U` generalized velocity coordinate values for the {{ coordbasecls }} objects such as {{ subhingescls }} and {{ modalbodycls }} deformation coordinates in the {{ mbodycls }} system. The state vector defines the current independent positions and velocities that fully describe the system's kinematic state at any given moment. A subtle issue of note in this context is the use of `global` and `local` coordinates at the {{ subhingecls }} level. As discussed earlier in the {ref}`attitude_sec` section, {{ subhingescls }} such as the {py:class}`SphericalSubhinge ` use `local coordinates` centered around `local charts` to get around the singularities associated with minimal-coordinate attitude representations. As a consequence, the numerical integrator state also makes use of the `local coordinates` for the {{ subhingescls }}. However, to shield external users from having to deal with this _under-the-hood_ issue, the {py:meth}`setState() ` method always expects the physically meaningful, and unambiguous, `global` coordinates when resetting the system state. Coordinate sanitization is automatically performed at the end of each integration step, and subhinges may re-center their local charts to convert the new _global_ coordinate values into _local_ coordinates as needed. ##### Tracing/debugging The {{ mmcls }} has a {py:attr}`trace_state_propagator ` variable that can used to enable trace messages for the {{ mmcls }}. This can be very helpful when debugging {{ mmcls }}-related and model-related issues. To get enable messages, set {py:attr}`trace_state_propagator ` to `True` and ensure the logger where you want to see the messages has a verbosity set to at least the {py:attr}`TRACE ` level. In addition, the {cpp:class}`CallbackRegistry `s in {py:attr}`fns ` have a {cpp:member}`trace_callback_registry ` variable that can be used to enable trace messages for each {cpp:class}`CallbackRegistry `. To enable the messages, set {cpp:member}`trace_callback_registry ` to `True` for each {cpp:class}`CallbackRegistry ` you want to see messages for, and ensure the logger you want to see messages on has a verbosity set to at least the {py:attr}`TRACE ` level. (state_space_sec)= #### Linearized state space models Since multibody dynamics models are nonlinear, so are the system level models. Designing controllers for such systems with nonlinear models is challenging. While there is a rich methodology for developing controllers for linear systems, this is not the case for nonlinear systems. Control system design and analysis therefore often rely on linearized state space models for control system design and frequency domain analysis (including creation of Bode plots). For nonlinear systems, linearized state space models at specific system configuration set points are utilized to support the control system design and analysis workflow. Towards this, there is built in support in {{ kdflex }} for creating linearized state space models from such {{ mmcls }} system level models for user specified sets of inputs and outputs about any configuration set point. The {py:meth}`Karana.Dynamics.Algorithms.stateSpaceGenerator` method can be used to create a {py:class}`Karana.Math.StateSpace` instance for a specified set of inputs and outputs. The {py:meth}`generate() ` can then be called on this object to linearize the system and compute the state space _A, B, C, D_ matrices. Since the state space representation is generated using linearization, its values are only valid in a small neighborhood of the configuration used to do the linearization. Additionally, such linearizations are typically meaningful only when the system is in a state of equilibrium, i.e. the accelerations are near zero. Thus, an equilibration step may be needed prior to state space model generation. Such linearized system based control-system design and analysis is complemented with more complete and thorough checkouts in closed-loop time-domain simulations with the full high-fidelity nonlinear model as described in the {ref}`system_level_models_sec` section. - - - (kmodels_sec)= ### KModel component models Simulation models provide a structured way to interact with objects in the simulation. For example, they can be used to do things like sense the position and velocity of a joint or a sensor node, apply forces via a joint or actuator node, or other tasks such as data logging. The {{ mmcls }} class provides an API for registering such component models so that their methods and functionality can be included in the overall system dynamics. At the C++ level, the component models are created by deriving from a templatized {cpp:class}`Karana::Models::KModel` class. All of the template arguments except the first indicate the types of states, parameters, etc. that the model has. The first template argument is the model class itself. The remaining template arguments are derived from specific other classes, e.g., the second template argument must be derived from {cpp:class}`Karana::Models::KModelParams`, the third from {cpp:class}`Karana::Models::KModelScratch`, etc. (kmodel_diagram_exs)= The diagram below shows the most complex KModel, one that utilizes all the template arguments. Each of these will be discussed in detail in the sections below. ```{eval-rst} .. image:: images/kmodel_full_light.svg :align: center :class: only-light .. image:: images/kmodel_full_dark.svg :align: center :class: only-dark ``` This diagram just uses {py:meth}`preHop ` and {py:meth}`preDeriv ` as some example methods `MyModel` may define. See the {ref}`kmodel_methods_sec` section below for more details. In many cases, not all of the template arguments will be needed. The default template arguments start with "No", e.g., {cpp:class}`NoParams `, to indicate they are not used. For example, suppose we had a model that just used discrete states. Then, the diagram would look like this: ```{eval-rst} .. image:: images/kmodel_discrete_states_light.svg :align: center :class: only-light .. image:: images/kmodel_discrete_states_dark.svg :align: center :class: only-dark ``` Notice that `No*` is used to indicate that the model does not have any params or scratch. The {cpp:class}`NoContinuousStates ` is not needed at the end, since not specifying this argument will pick up the default (which is {cpp:class}`NoContinuousStates `). When specified in this way, the `params`, `scratch`, and `continuous_states` members will all be `nullptr`s. (kmodel_methods_sec)= #### Model methods Their are six different methods a user can define in a model that will run at various points during the simulation. The user does not need to define all six methods, but must define at least one of them. A model can skip the implementation of unneeded methods. The model's {py:meth}`preDeriv ` and {py:meth}`postDeriv ` methods will run before and after the {{ mbodycls }} derivative, respectively. The {py:meth}`preHop ` and {py:meth}`postHop ` methods will run before and after each simulation hop, respectively. The {py:meth}`preModelStep ` and {py:meth}`postModelStep ` step methods will run before and after each of the model’s steps, respectively. The contents of these methods are model-dependent, and up to the model writer to define. The {py:meth}`preDeriv ` and {py:meth}`postDeriv ` methods are called at times dependent on the integrator type. The {py:meth}`Karana.Models.SpringDamper.preDeriv` is an example implementation of the {py:meth}`preDeriv ` method. The {py:meth}`preHop ` and {py:meth}`postHop ` methods will be called at each hop on the timeline determined by the scheduler. The {py:meth}`Karana.Models.SyncRealTime.preHop` method is an example of a {py:meth}`preHop ` method, and the {py:meth}`Karana.Models.UpdateProxyScene.postHop` method is an example of a {py:meth}`postHop ` method. Environment interaction models, such as gravity, do not have a pre-determined step size that they need to run at. However device models often do have specific rates that they need to run at based on the hardware design. The multi-rate step size can be specified for such models. The {py:meth}`preModelStep ` and {py:meth}`postModelStep ` methods will run at rates specified by the model's multi-rate step size. Model step sizes fall into one of two categories: uniform and non-uniform. Uniform model steps are the most common and are step sizes defined by a constant period. For these models, simply use the {{ kmodelcls }} {py:meth}`setPeriod() ` method to set the model’s period. Non-uniform models can use the {py:meth}`nextModelStepTime ` method to define their non-uniform rate. This method takes as input the current time and outputs the next time at which the model should step. The {py:meth}`Karana.Models.DataLogger.postModelStep` is an example of a {py:meth}`postModelStep ` step method implementation. (kmodel_params_sec)= #### Model parameters Oftentimes, models utilize parameters in their calculations. For example, a spring damper model may have the spring constant and damping value as parameters. In order to define parameters for your model, create a parameter class that is derived from {cpp:class}`Karana::Models::KModelParams` and contains the parameters you want. Then, set the {cpp:member}`params ` member variable for the model equal to a `ks_ptr` instance of this class during the model’s construction or in its `create` method. Model parameter classes should also override the {cpp:func}`isReady ` method. This method should contain logic to verify that the derived parameter class’s parameters have been initialized. For parameters that are real numbers or math objects such as vectors or matrices, see {ref}`Math ` for details on handling of not ready values. The {cpp:func}`isReady ` method will be invoked automatically for all registered models as part of the overall simulation initialization. It will be able to flag any not ready parameters for the model instance. This check will help verify for anyone creating an instance of the {{ kmodelcls }} that they have not accidentally skipped the initialization of all the necessary parameters. (kmodel_scratch_sec)= #### Model scratch Sometimes, it is useful to store intermediate results in a model for debugging purposes. The {cpp:class}`Karana::Models::KModelScratch` class gives one a place to store such information. To use it, create a class that is derived from {cpp:class}`KModelScratch ` and add class members for the scratch values the model should store. This can be accessed using the {cpp:member}`scratch ` member on the model itself. The {cpp:member}`scratch ` value should be initialized when the model is created, either in the constructor directly or in the `create` method. (kmodel_discrete_states_sec)= #### Discrete states Models may have discrete states (states that should not be integrated, but need to be stored and updated as time progresses). Examples of discrete states are the on/off state of a valve or PID integral gain, where the PID model is keeping track of the integrated error separate from the {{ mmcls }}'s integrator, e.g., the model is using an Euler integrator with a fixed step set by the model period. The {cpp:class}`Karana::Models::KModelDiscreteStates` class gives one a place to store such information. To use it, create a class that is derived from {cpp:class}`KModelDiscreteStates ` and add class members for the discrete states the model needs to store. This class can be accessed from the {cpp:member}`discrete_states ` member of the model itself. The {cpp:member}`discrete_states ` value should be initialized when the model is created, either in the constructor directly or in the `create` method. Discrete states should be updated by the model in one of the following methods: * {py:meth}`preModelStep ` * {py:meth}`preHop ` * {py:meth}`postHop ` * {py:meth}`postModelStep ` These methods run at the hop-boundaries. The method used to update the discrete states is model dependent. Models that need to update states at a rate set by the model should use either {py:meth}`preModelStep ` or {py:meth}`postModelStep `. Models that need to update their state more frequently e.g., whenever there is a "good" state, should use the {py:meth}`preHop ` or {py:meth}`postHop ` methods. A "good" state refers to one that the integrator will keep, as opposed to states in the {py:meth}`preDeriv ` and {py:meth}`postDeriv ` method that might represent states where an adaptive integrator has taken too large of a step and will be rolled back. Using the `pre*` vs. `post*` methods is also model dependent. The `pre*` methods run first, and so are a good choices when one wants to update a state immediately when a command comes in, e.g., turning on power or opening a valve. The {py:meth}`preHop ` method does not know the exact hop size. It has the intended hop size, but this may be modified by zero-crossing events, so if the true hop size is needed, then the {py:meth}`postHop ` method should be used over the {py:meth}`preHop ` method. The {py:meth}`preModelStep ` and {py:meth}`postModelStep ` do not have this same issue, as their size is set by the model and not modified by zero-crossing events. (kmodel_continuous_states_sec)= #### Continuous states Models may have continuous states (states that should be part of the {{ mmcls }}'s integrator). Examples of continuous states are values used to calculate noise in an IMU model or values used to calculate fuel burned in a thruster model. The model's continuous states must be expressible as a vector of doubles (a `km::Vec` in C\++ or `numpy.array` in Python). To use continuous states, derive a class from {cpp:class}`Karana::Models::KModelContinuousStates` and override the {cpp:func}`getX `, {cpp:func}`setX `, and {cpp:func}`getdX ` methods. These retrieve the current continuous states from the class, set the continuous states for the class, and retrieve the derivatives of the continuous states from the class, respectively. All of these methods make use of `Eigen::Ref` at the C++ level, so they should be calculated and stored on the class itself, and an `Eigen::Ref` returned from them if the model is a C++ model. The continuous state derivatives should be updated in either the {py:meth}`preDeriv ` or {py:meth}`postDeriv ` methods. A method at the derivative level must be used, as these are used by the integrator to propagate the state forward in time. The {py:meth}`preDeriv ` runs before the multibody dynamics and the {py:meth}`postDeriv ` runs after the multibody dynamics. The method used by a model often depends on whether the continuous state derivatives need updated multibody dynamics values or not. (model_register_sec)= #### Registering and unregistering models To affect the simulation, models must be registered with a {{ mmcls }}. It is recommended that when one writes a model they add a {py:meth}`ModelManager.registerModel ` call in the model's `create` method (or the `__init__` method if it is a Python model) so that this is done automatically. The {{ mmcls }} will run the model’s methods at the appropriate times. At any time, models can also be unregistered from the {{ mmcls }}. This will effectively disable the model. To register or unregister a model, simply use the {{ mmcls }} {py:meth}`registerModel() ` or {py:meth}`unregisterModel() ` methods, respectively. This provides a convenient way of activating and deactivating models at run-time in response to simulation mode changes. #### Inputs and outputs The philosophy for connecting {{ kmodelscls }} is to use shared pointers to C++ classes. This is preferable to completely encapsulating code within a KModel, as it makes it: * Easier to unit test. * Easier to update with new versions in the future. * Simplifies architecture by removing the need for a separate, shared signal buffer. This is understood easiest with an example. Suppose one wants to drive a PID KModel with the output of some profile generator. At first, it may be tempting to create an architecture similar to the one shown in the diagram below: ```{eval-rst} .. image:: images/kmodel_pid_connection_1_light.svg :align: center :class: only-light .. image:: images/kmodel_pid_connection_1_dark.svg :align: center :class: only-dark ``` However, this fully encapsulates the profile generator code into the KModel. This is undesirable for a few reasons: 1. If at a later time one wants to swap to a different profile generator, e.g., from a linear profile generator to a cubic profile generator, they need to swap out the entire KModel. 1. It makes unit-testing the profile generator code more difficult, as one needs to create a KModel and supporting infrastructure to be able to test it, rather than being able to just test the profile generator code in isolation. 1. This requires something to tie the buffers between the {{ kmodelscls }} together, or a small class that contains the buffer that the two can share. Since something like 3 above is needed anyway, is is simpler and beneficial to instead assemble the models as shown below: ```{eval-rst} .. image:: images/kmodel_pid_connection_2_light.svg :align: center :class: only-light .. image:: images/kmodel_pid_connection_2_dark.svg :align: center :class: only-dark ``` This has a few advantages: 1. If one has a base class for profile generators, and they write each of these {{ kmodelscls }} to take an instance of that base class, then one can easily swap out profile generators if need be. There is no need to keep creating a new {{ kmodelcls }} for each profile generator. This reduces the amount of code that needs to be written as new profile generators are added in the future. 1. In this architecture, one can easily write unit tests for a profile generator, as it can be created without any {{ kmodelscls }} or associated supporting infrastructure: a `ModelManager, etc. 1. The profile generator itself is the shared buffer, provided that the base class has a getter/setter for the profile value. 1. This is more extensible than the original diagram. For example, suppose that instead of updating the generator at a fixed interval with a {{ kmodelcls }}, one wants to update it at zero crossing events. One can easily do this by just wrapping the profile generator in a zero crossing event rather than a {{ kmodelcls }}, and the PID can still use the output of that profile generator as it normally would. Another example is a constant profile generator that never needs to be updated, or is being updated by the user from the Python REPL as they are debugging. This second architecture greatly expands the types of scenarios this profile generator and PID controller can used it in, as the profile generator it is no longer constrained to the update restrictions of {{ kmodelscls }}. #### C++ base class The {cpp:class}`KModel ` base class in c++ utilizes the _curiously recurring template pattern_. This pattern allows the base class to detect which methods a derived model class has defined. This provides a way to register only the methods that a user has defined for a model with the model manager when registering the model. In addition, it also provides other conveniences, such as the `params`, `scratch`, etc. member variables' types matching the derived model's class' types. ```c++ // Params class for MyModel class MyModelParams : public Karana::Models::KModelParams { public: // Override the `isReady` method bool isReady() const final; // Include parameters here float: my_param ... }; class MyModel : public Karana::Models::KModel { // Define methods for the model such as `preDeriv` or `postModelStep` here ... }; ``` The {ref}`diagrams ` at the beginning of the {ref}`kmodels_sec` section show some other examples of using the templated {cpp:class}`KModel ` class. #### Python base class The {py:class}`Karana.Models.KModel` bass class in Python compares derived methods with their base class counterparts to determine whether those methods have been overridden or not. Adding params, scratch, discrete states, and continuous states is done by populating member variables in the constructor. The typing for these members (`params`, `scratch`, `discrete_states`, and `continuous_states`) is provided by a generic in the base {py:class}`KModel ` class. Similar to C++, the `No*` versions should be used if the member will not be used. The Python versions of these classes are similar to the C++ with one exception; the {py:class}`KModelParams ` class is defined from {py:class}`DataStruct `, so the parameters can be defined with type hints (that are verified at runtime) by adding them as one would do with any other {py:class}`DataStruct `. ```python # Params class for MyModel from Karana.Models import KModelParams, KModel, NoScratch, NoDiscreteStates, NoContinuousStates class MyKModelParams(KModelParams): # Override the `isReady` method def isReady()-> bool: ... # Include parameters here my_param: float ... class MyModel(KModel[KModelParams, NoScratch, NoDiscreteStates, NoContinuousStates]): def __init__(self, name: str, sp: ModelManager): super().__init__(name, sp) self.params = MyKModelParams(f"{name}_params") sp.registerModel(self) # Define methods for the model such as `preDeriv` or `postModelStep` here ... ``` #### Debugging Every KModel has a {cpp:member}`debug_model ` boolean member variable. This variable is used to enable/disable debug messages for the model. When running a simulation, to get debug messages for a particular model, one should ensure they have a logger whose verbosity is set to at least the debug level and then set the {cpp:member}`debug_model ` member variable to `true`. To instrument models with these messages, in C++ use ```cpp if (debug_model) [[unlikely]] { double my_var = someLongFunction(); stdDebugMsg(std::format("my_var is : {}", my_var); } ``` and in Python use ```python if self.debug_model: my_var = someLongFunction() self.stdDebugMsg(f"my_var: {my_var}") ``` KModels and associated classes (params, scratch, discrete states, and continuous states) can also override `dumpString` methods. These methods will be called and combined whenever {py:meth}`KModel.dump() ` is called, which will display an overall dump for the model. In C++, these methods should be overridden like ```cpp class M1Scratch : public kmo::KModelScratch { public: using KModelScratch::KModelScratch; double t = 0.0; // Example dumpString override. The returned string should utilize the prefix in front of each line. // This can be done in a similar fashion for params, discrete states, and continuous states. std::string dumpString(std::string_view prefix = "", const kc::Base::DumpOptions * /*options*/ = nullptr) const override { return std::format("{}t: {}\n", prefix, t); } }; class M1 : public kmo::KModel { public: // Other model-related code, constructor, etc. here. // Example dumpString override for the model std::string dumpString(std::string_view prefix = "", const kc::Base::DumpOptions * /*options*/ = nullptr) const override { std::string ret = std::string(prefix); // Any custom information should go here. ret.append(std::format("Dumping model {}\n", name())); // Call the parent method. Here, we add a " " to the prefix to tab over // everything the parent adds relative to the statement we just printed above. // The parent method will take care of calling the dumpString methods for the // params, scratch, discrete states, and continuous states if appropriate. ret.append(kmo::KModel::dumpString( std::format("{} ", prefix))); return ret; } }; ``` Similarly, in Python, these methods can be overridden like ```python class M1Scratch(KModelScratch): def __init__(self): super().__init__("scratch") self.t: float = 0.0 def dumpString(self, prefix: str = "", options: DumpOptionsBase | None = None) -> str: return f"{prefix}t: {self.t}\n" class M1(KModel[KModelParams, M1Scratch, M1States, KModelContinuousStates]): # Other model-related code, constructor, etc. here. def dumpString(self, prefix: str = "", options: DumpOptionsBase | None = None) -> str: # Add our own custom dumpString line. Then, call the parent method with a tab added # to the prefix. The parent method will take care of calling the params, scratch, # discrete states, and continuous states dumpString methods if appropriate, and the # extra tab added to the prefix will tab them over relative to the custom line we # added here. return f"{prefix}Dumping mode {self.name()}\n" + super().dumpString( prefix + " ", options ) ``` (kmodel_datastructs_sec)= #### KModel DataStructs In order to serialize/deserialize a `KModel`, it must have an associated `DataStruct`: see {ref}`datastruct_sec` for details. If the model has parameters, a `DataStruct` will also need to be created for its parameters, and the parameter `DataStruct` model should be nested in the KModel's `DataStruct`, i.e., the parameter `DataStruct` should be a field of the KModel `DataStruct`. The KModel `DataStruct` should be derived from `KModelDS` and override the `toModel` method. The `toModel` method is the method that creates a `KModel` instance from the `DataStruct`. In some cases, this method will need other simulation objects, e.g., `Node`s or `PhysicalBody`s, used by the `KModel`s constructor. The `toModel` has a `SubGraphDS` as an argument that is used for this purpose. The {py:meth}`IdMixin.getObjectsCreatedById ` static method can be used to link IDs used by a `DataStruct` to object instances in the current simulation. To utilize this, the KModel's `DataStruct` should contain fields for the IDs of the objects it needs during construction. As an example, consider the portion of `SpringDamperDS` shown below. This model uses two nodes in its constructor. The IDs of these nodes are saved as fields in the `DataStruct`. Then, they are used in the `toModel` method to find the actual node objects. ```python class SpringDamperDS(KModelDS[SpringDamper]): """DataStruct for SpringDamper model. Parameters ---------- n1 : int ID of the first node of the spring damper. n2 : int ID of the second node of the spring damper. params : SpringDamperParamsDS Parameters for the SpringDamper model. """ n1: int n2: int params: SpringDamperParamsDS def toModel(self, mm: ModelManager) -> SpringDamper: """Create a SpringDamper model instance from this SpringDamperDS. Parameters ---------- mm : ModelManager The ModelManager to associate this model with. Returns ------- SpringDamper An instance of the SpringDamper model. """ # Find the nodes needed for the model by ID dark = IdMixin.getObjectsCreatedById(self.n1) if dark is None or len(dark) != 1: raise ValueError( f"Could not find unique object associated with original node 1. Found {dark} for node {self.n1}." ) else: n1 = dark[0] dark = IdMixin.getObjectsCreatedById(self.n2) if dark is None or len(dark) != 1: raise ValueError( f"Could not find unique object associated with original node 2. Found {dark} for node {self.n1}." ) else: n2 = dark[0] # Create the model model = SpringDamper(self.name, mm, n1, n2) self.params.setParams(model.params) if self.period > 0: model.setPeriod(self.period) # Keep track of objects created from this DS self.addObjectFromId(model) return model ``` In addition, the `KModel` itself should define the `toDS` method. The "Karana/SOADyn/pybind11KModel.h" file includes helper methods to create these in `pybind11` code. For example, see the excerpt from `GeneralKModels_Py.cc` below. ```c++ #include "Karana/SOADyn/pybind11KModel.h" // This macro must be called before the associated `addToDSMethod` is called. // This creates the appropriate type caster. ADD_TO_DS_METHOD_CASTER(SpringDamper, Karana.Models.GeneralKModels_types, SpringDamperDS) PYBIND11_MODULE(GeneralKModels_Py, m) { // Define SpringDamper class for pybind11 like any other class ... // Add the toDS method using the helper kd::addToDSMethod( PySpringDamper, "SpringDamper", "Karana.Models.GeneralKModels_types", "SpringDamperDS", "fromModel"); } ``` - - - (time_domain_sim_sec)= ### Time-domain simulations The primary requirement for time-domain simulations is the accurate propagation of the system state, {{ kdflex }} supports this primary objective while also accommodating several other important simulation requirements such as: - support for integrating multibody dynamics with component device and environment interaction models - support for using multi-rate models within the system simulation - support for terminating integration steps at zero-crossing boundaries - support for modeling arbitrary time delays - support for step validation and rollback - support for handling discrete, discontinuous run-time changes to system dynamics - support for activating and deactivating component models at run-time - scheduling and registering of arbitrary timed events during a simulation run - support for the use of different types of multibody dynamics solvers - support for accurate and high order numerical integrators as needed, including variable-step and stiff integrators - providing deterministic and unambiguous handling of time representations - support verification of parameter initialization across all model instances - support for periodic data logging - support for live updates to visualization and data plots at run-time - fast computational performance for closed-loop simulation use - embeddable for hardware-in-the-loop simulations The {{ spcls }} class supports these features required for the development of realistic and complex system level time-domain simulations. (integrator_sec)= #### Numerical integration For time-domain simulations, the class supports the instantiation of a {py:class}`Karana.Math.Integrator` numerical integrator from a large suite of fixed and variable step integrators for time domain state propagation and simulation. The integrator options range from the fixed step {py:class}`Karana.Math.EulerIntegrator` Euler 1-step integrator and the fourth order {py:class}`Karana.Math.RK4Integrator` integrators, to the {py:class}`Karana.Math.CVodeIntegrator` and {py:class}`Karana.Math.ArkExplicitIntegrator` large family of variable step integrators from Lawrence Livermore's [SUNDIALS](https://computing.llnl.gov/projects/sundials) numerical integrator suite. (state_propagator_sec)= #### The state propagator The {{ spcls }} class is derived from the {{ mmcls }} class and adds the functionality for propagating the {{ mmcls }} state in time using numerical integrators. ##### State propagation hops The {{ spcls }} {py:meth}`advanceTo() ` method can be used to request the advancement of the system state over time for a simulation scenario. Such a time advancement request is broken up by the state propagator into discrete state propagation time intervals referred to as time `hops`. The different call rate times of the component models are used to determine the (irregular) time `hop` sizes required to meet the overall propagation time barrier requirements. The numerical integrator propagates the system state across individual hops by using the state propagator to orchestrate calls to the methods of the registered {{ kmodelcls }} component models and the multibody dynamics in the correct order to compute the state derivatives. (zero_crossings_sec)= ##### Zero-crossings In addition to this normal state propagation process, the {{ spcls }} also has support for `zero-crossing detection`. This phrasing is used for situations where a simulation hop is required to terminate when an implicitly defined, state dependent, condition is met. ```{eval-rst} .. _zerocrs: .. figure:: images/zero_crossing.png :width: 500px :align: center :alt: zero crossing Terminating time hops at zero-crossing events ``` Examples of this are body collision events, state-dependent triggers etc. When a zero-crossing condition method is registered with the state propagator, the propagator will monitor these zero-crossing conditions, and when triggered, will use a regula-falsi iteration to find the precise time boundary when the condition triggered to terminate the hop, execute appropriate registered actions before proceeding on to the next hop. (timed_events_sec)= ##### Timed events The {{ spcls }} also allows the registering of user defined {py:class}`Karana.Dynamics.TimedEvent` instances along with associated custom handlers. The time event executions can be one time, irregular, or periodic. The {{ spcls }} will ensure that hops terminate at the registered time boundaries for such {{ teventscls }} to invoke their associated handlers. Examples of use of {{ teventscls }} are for simulation needs such as data logging, updating visualization etc. This flexible capability allows user to incorporate arbitrary events into their simulations. {{ kdflex }} uses `numpy.timedelta64` nanosecond precision class for time ({cpp:type}`Karana::Math::Ktime` at the C++ level). Its use avoids any ambiguities that can occur from double precision representations of time. - - - (scenetop_sec)= ## Geometry Layer (scene_layer_sec)= ### Scene layer {{ kdflex }} defines an abstract interface for working with 3D geometries attached to the physical {{ physbodiescls }} in the {{ mbodycls }} system. This abstract interface and its various implementations make up the {py:class}`Karana.Scene.Scene` _Scene_ layer. {{ kdflex }} provides various {{ scenecls }} implementations for different purposes and third-party backends. For example, {py:class}`Karana.Scene.WebScene` is a {{ scenecls }} implementation that provides a live 3d visualization of the simulation in the web browser using [three.js](https://threejs.org/). {py:class}`Karana.Scene.CoalScene` implements the same {{ scenecls }} interface for collision detection using [Coal](https://github.com/coal-library/coal). This architecture is extensible and allows the addition of new backend scenes by simply implementing new specializations of the {{ scenecls }} class. ```{eval-rst} .. _vizupd: .. figure:: images/simgraphics.png :width: 600px :align: center :alt: viz updates Updating of visualization scene from simulation data ``` (proxyscene_sec)= #### ProxyScene There is often a need for multiple {{ scenecls }} instances in a single simulation - for example the need for both collision detection and 3d graphics. {py:class}`Karana.Scene.ProxyScene` is a {{ scenecls }} implementation that manages other {{ scenecls }} implementations. The {{ mbodycls }} and other multibody classes by and large interface and interact with the {{ pscenecls }} manager class, and use it to mediate the interactions with the registered {{ scenecls }} backends. We refer to Scene implementations managed by {{ pscenecls }} as _client Scenes_. In addition to synchronizing state across client Scenes, {{ pscenecls }} also acts as a bridge between the {{ framecls }} and {{ scenecls }} layers. That is, {{ pscenecls }} uniquely supports attaching objects from the Scene layer to Frames from the Frame layer. The upshot is, as a {{ framecls }} moves, the {{ pscenecls }} will correspondingly move attached geometries in all of the client Scenes. ```{eval-rst} .. _proxysceneimg: .. figure:: images/proxyscene.png :width: 600px :align: center :alt: proxysceneimg ProxyScene interface to client Scenes ``` A typical setup in a full simulation will look like this: 1. Create a {{ mbody }} 2. Create a {{ pscenecls }} 3. Create Scene objects in the {{ pscenecls }} instance, and attach them to {{ physbodiescls }} and other frames 4. Create various client Scenes as needed and register them using the {{ pscenecls }} {py:meth}`registerClientScene() ` method. This order doesn't need to be strictly adhered to. It's also reasonable to populate the {{ mbody }} and {{ pscenecls }} in tandem, creating each body and its attached geometries simultaneously. A client {{ scenecls }} can also be registered at any time. {{ pscenecls }} immediately syncs up a newly registered client {{ scenecls }} and continues to update the client {{ scenecls }} with any further changes. As a simulation evolves and a {{ framecls }} moves around, {{ pscenecls }} will need to be told when to update the corresponding transforms in its client Scenes. To sync a client {{ scenecls }}'s transforms with the Frame layer, call the {{ pscenecls }} {py:meth}`update() ` with the client Scene in question as an argument. To update all client Scenes, just call {py:meth}`update() ` with no arguments. Each {py:class}`Karana.Scene.ProxySceneNode` and {py:class}`Karana.Scene.ProxyScenePart` implementation in {{ pscenecls }} has a corresponding {{ scenenodecls }} and {{ scenepartcls }} implementation in each client scene which should be used exclusively with that client Scene. In most cases, where {{ pscenecls }} is managing the Scene layer, this means creating instances of {{ pscenenodecls }} and {{ pscenepartcls }} to populate the ProxyScene (and consequently all of the client Scenes). #### Scene objects Each {{ scenecls }} is a container for instances of {py:class}`Karana.Scene.SceneNode` and {py:class}`Karana.Scene.ScenePart`. Each {{ scenenodecls }} has a {py:class}`Karana.Math.SimTran`, which consists of a position, rotation, and uniform scale factor relative to a parent. {{ scenenodecls }} instances may be attached to one another with {py:meth}`attachTo() `, giving the overall Scene a tree structure (with an implicit root node). Each {{ scenenodecls }} also has a visibility flag, which can be used to hide the node and all of its descendants. A {{ scenepartcls }} is a {{ scenenodecls }} that additionally has a 3D geometry, which can either be a primitive shape such as a sphere or cylinder, or a triangular mesh. A {{ scenepartcls }} also has a material describing its surface properties and appearance. Typically a scene will consist of a tree, where branches are {{ scenenodecls }} instances, and leaves are {{ scenepartcls }} instances. ```{eval-rst} .. image:: images/scene_tree_light.png :align: center :class: only-light .. image:: images/scene_tree_dark.png :align: center :class: only-dark ``` #### Geometries Each {{ scenepartcls }} has a {py:class}`Karana.Scene.AbstractStaticGeometry` defining its 3D shape. A geometry may be an analytical primitive shape (e.g., {py:class}`Karana.Scene.BoxGeometry`, {py:class}`Karana.Scene.SphereGeometry`), or a triangular mesh (typically {py:class}`Karana.Scene.StaticMeshGeometry`). Creating a {{ scenepartcls }} requires a `VarStaticGeometry`, which is a variant that can be a shared pointer to any of the geometry types defined in {{ kdflex }}: - {py:class}`Karana.Scene.BoxGeometry` - {py:class}`Karana.Scene.CapsuleGeometry` - {py:class}`Karana.Scene.ConeGeometry` - {py:class}`Karana.Scene.CylinderGeometry` - {py:class}`Karana.Scene.RoundFrustumGeometry` - {py:class}`Karana.Scene.SphereGeometry` - {py:class}`Karana.Scene.StaticMeshGeometry` Additional user-defined geometry types can also be used by first upcasting to {py:class}`Karana.Scene.AbstractStaticGeometry`. The same geometry instance may be passed to multiple {{ scenepartcls }} instances for easy [geometry instancing](https://en.wikipedia.org/wiki/Geometry_instancing). #### Materials Each {{ scenepartcls }} also has a material. There are currently two choices for material type: - {py:class}`Karana.Scene.PhysicalMaterial`: this is a conventional PBR material with fields such as base color, metalness, and roughness. - {py:class}`Karana.Scene.PhongMaterial`: this is commonly used material with fields for diffuse, specular, emissive colors and textures. Creating a material is a two-step process. First, create a corresponding _info_ {py:class}`Karana.Scene.PhysicalMaterialInfo` or {py:class}`Karana.Scene.PhongMaterialInfo`, and update its fields as desired. Then to create the actual material, pass the info instance into the material's `create` method. Creating a {{ scenepartcls }} requires a `VarMaterial`, which is a variant that can be a shared pointer to either material type. A part's material may be modified at any time by calling {py:meth}`setMaterial() ` to assign a new material instance to the part. #### Layers Each {{ scenepartcls }} belongs to _layers_. Layers are implemented using a `layer_t` which is a 32-bit unsigned integer, where each bit represents a different layer that a {{ scenepartcls }} may inhabit. Scene parts may then be filtered using a `layer_t` bitmask, such that a part is filtered out if it doesn't have any layers in common with the bitmask. The least significant twenty-four bits are set aside for user defined layers and will not be used internally by {{ kdflex }}. For convenience {{ kdflex }} defines the constants {py:const}`Karana.Scene.LAYER_CUSTOM0` through {py:const}`Karana.Scene.LAYER_CUSTOM23`, which are simply the first twenty-four powers of two. The most significant eight bits are reserved for layers defined by {{ kdflex }}, with constants {py:const}`Karana.Scene.LAYER_RESERVED0` through {py:const}`Karana.Scene.LAYER_RESERVED7`. {{kdflex}} currently makes use of the first four reserved layers: - {py:const}`Karana.Scene.LAYER_PHYSICAL_GRAPHICS`: reserved for visualization of real, physical objects - {py:const}`Karana.Scene.LAYER_ORNAMENTAL`: reserved for visualization of ornamental objects not suitable for sensor simulation - {py:const}`Karana.Scene.LAYER_STICK_FIGURE`: reserved for scene parts generated by {py:meth}`createStickParts() ` - {py:const}`Karana.Scene.LAYER_COLLISION`: reserved for parts used for collision detection {{ kdflex }} also defines several layer sets: - {py:const}`Karana.Scene.LAYER_NONE`: no layers, that is, the number zero. - {py:const}`Karana.Scene.LAYER_ALL`: all layers, that is, the maximum 32-bit unsigned integer. - {py:const}`Karana.Scene.LAYER_PHYSICAL`: combination of {py:const}`LAYER_PHYSICAL_GRAPHICS ` and {py:const}`LAYER_PHYSICAL_COLLISION `. This may be useful if the same geometry is suitable for both graphics and collision detection. - {py:const}`Karana.Scene.LAYER_GRAPHICS`: combination of all reserved graphical layers ({py:const}`LAYER_PHYSICAL_GRAPHICS `, {py:const}`LAYER_PHYSICAL_ORNAMENTAL `, {py:const}`LAYER_PHYSICAL_FIGURE `). This is the default bitmask used for visualization. Filtering based on layers is done in two places: 1. {py:meth}`ProxyScene.registerClientScene ` takes a `layer_t` argument, defaulting to {py:const}`LAYER_PHYSICAL_ALL `. This can be used to filter which scene parts should be implemented in a registered client scene. 2. Client scenes may filter parts, depending on the application. For example, {py:class}`Karana.Scene.GraphicalSceneCamera` has a `layer_t` bitmask controlling which layers should be rendered (defaulting to {py:const}`Karana.Scene.LAYER_GRAPHICS). #### Uniform and intrinsic scale Each {{ scenenodecls }} has a {py:class}`Karana.Math.SimTran`, which includes a scale factor relative to the node's parent. This can be used to scale entire subtrees within a scene. The scale factor is uniform, meaning it is a scalar applied in all three directions. Without this limitation it would be possible to introduce skew by chaining together rotations and non-uniform scaling transforms, which is problematic for many scene backends. To support non-uniform scaling, the {{ scenepartcls }} class has the notion of _intrinsic scaling_. The intrinsic scaling can be non-uniform, and is only applied directly to the part and is ignored when computing the world transform for any descendant parts. A part's intrinsic scale can be accessed through the {py:meth}`getIntrinsicScale ` and {py:meth}`setIntrinsicScale ` methods. #### Importing assets Scene assets such as geometries and materials are often created offline and need to be loaded from files. {{ kdflex }} defines an abstract interface {py:class}`Karana.Scene.AbstractImporter` and implementation {py:class}`Karana.Scene.AssimpImporter`, which can be used to load asset files into native {{ kdflex }} types. To load a file, first create an `AssimpImporter` instance, then invoke the {py:meth}`importFrom ` method with a filepath argument. This returns a {py:class}`Karana.Scene.ImportResult` containing any assets found in the given file. There are multitude of widely used file formats, whose contents may range from a single asset to an entire serialized scene, so the `ImportResult` may consist of any number of assets in general. #### Parts managed by a multibody Scene parts attached to a {{ physbodycls }} should preferably be managed by the body by calling {py:meth}`PhysicalBody.addScenePartSpec `. The {py:class}`Karana.Scene.ScenePartSpec` argument is a plain struct with all of the information needed to eventually create a {{ scenepartcls }}. This has a number of benefits: - The scene part's lifetime is automatically managed by the body - The scene part can be specified immediately while deferring its creation until a {{ pscenecls }} is registered with the multibody via {py:meth}`setScene `. - The scene part will be included when serializing the multibody - - - (visualization_sec)= ### Visualization {{ kdflex }} provides a Scene implementation for 3d visualization called {py:class}`Karana.Scene.WebScene`. {py:class}`Karana.Scene.WebScene` is implemented on top of [three.js](https://threejs.org/), and uses WebGL to do hardware accelerated 3d-rendering in the web browser. This makes it possible to bring up a live, interactive visualization of a simulation on a remote machine by simply opening the correct URL in a web browser. {{ kdflex }} also comes with a standalone [electron](https://www.electronjs.org/) application that can be used for local 3d visualization in lieu of the web browser option. In most cases, where {{ pscenecls }} is managing the {{ scenecls }} layer, graphics can be enabled by creating the {py:class}`Karana.Scene.WebScene` client scene and registering it via {{ pscenecls }} {py:meth}`registerClientScene() `. This will automatically populate the {py:class}`Karana.Scene.WebScene` client with whatever geometries have been set up in {{ pscenecls }}. (webui_sec)= #### WebUI To instantiate {py:class}`Karana.Scene.WebScene`, first create an instance of {py:class}`Karana.WebUI.Server`, which manages communication between the simulation process and any number of connected front-end clients. {py:class}`Karana.WebUI.Server` accepts connections on a port specified at construction time and queryable by calling {py:meth}`Karana.WebUI.Server.port`. To open a browser-based front-end, simply enter the URL to access the port over http. For example, if {py:meth}`Karana.WebUI.Server.port` returns 9090 and the simulation is running locally, enter http://localhost:9090 in your web browser's address bar. If the simulation is running on a remote machine, replace localhost with an IP address or domain that resolves to the simulation server. It's not always possible for {py:class}`Karana.WebUI.Server` to know how you will be accessing it, but {py:meth}`Karana.WebUI.Server.guessUrl` returns a best-guess URL string. To open the standalone front-end for a locally running simulation, call {py:meth}`Karana.WebUI.Server.launchLocalClient`. This launches an electron application on the same machine as the simulation that automatically connects to the simulation server. The [Visualization](generated/notebooks/example_visualization/notebook.ipynb) notebook provides an example of using the {{ kdflex }} visualization capability . - - - (collision_dynamics_sec)= ### Collision dynamics Collision detection is handled by collision engines. These detect whether objects are colliding, and provide information like the penetration of the collision and the collision normal. The {py:class}`Karana.Collision.FrameCollider` class is a generic way to capture the information coming from collision engines. Once this information is gathered, contact forces can be calculated from it. The {py:class}`Karana.Models.PenaltyContact` model passes the information from the {py:class}`FrameCollider `s to classes derived from {py:class}`Karana.Collision.ContactForceBase`: these are the classes which calculate and apply the contact forces. The diagram below illustrates the flow of data graphically: ```{eval-rst} .. image:: images/collision_diagram_light.svg :align: center :class: only-light .. image:: images/collision_diagram_dark.svg :align: center :class: only-dark ``` Each of the pieces of this pipeline are explained in detail in the subsections below. #### Collision engines Currently, {py:class}`Karana.Scene.CoalScene` is the only collision engine that comes packaged with {{ kdflex }}. However, more may be added in the future. {py:class}`Karana.Scene.CollisionScene` is an extension of {{ scenecls }} that adds a more specialized interface for querying collisions that can be fed back into the simulation for contact dynamics. {py:class}`Karana.Scene.CoalScene` is built upon the [Coal](https://github.com/coal-library/coal) project that implements {py:class}`Karana.Scene.CollisionScene`. In most cases, where {{ pscenecls }} is managing the Scene layer, collisions can be enabled by creating a {py:class}`Karana.Scene.CoalScene` client scene and registering it via {{ pscenecls }} {py:meth}`registerClientScene() `. This will automatically populate the CoalScene client with whatever geometries have been set up in {{ pscenecls }}. #### FrameCollider The {py:class}`Karana.Collision.FrameCollider` class is a helper that maps collisions from a {py:class}`Karana.Scene.CollisionScene` to the {{ framecls }} layer. {py:meth}`Karana.Collision.FrameCollider.collide` sweeps through the {py:class}`Karana.Scene.CoalScene` and invokes the registered callback for each pair of contact points between {{ framecls }}s due to their attached geometries. This process provides a pathway for returning information and invoking simulation level actions based on events triggered in the collision client scenes. To ignore collisions between frame pairs, call the {py:meth}`FrameCollider.ignoreFramePair ` method. In many cases a multibody model may have overlapping geometries between which collisions must be ignored. The convenience method {py:meth}`FrameCollider.ignoreAllCurrentlyTouchingPairs ` sweeps through and ignores all such pairs. #### PenaltyContact The {py:class}`Karana.Models.PenaltyContact` uses {py:class}`FrameCollider `(s) and applies penalty forces to {{ physbodiescls }} based on their contact penetration, relative velocity, etc. The {py:class}`PenaltyContact ` model uses a {py:class}`ContactForceBase ` class to perform the contact force calculations. This can be any class derived from {py:class}`ContactForceBase `. The {py:class}`FrameCollider `(s) and {py:class}`ContactForceBase ` class are all passed as constructor arguments to {py:class}`PenaltyContact `. For example, the {py:class}`Karana.Collision.HuntCrossley` contact force calculation uses a Hunt-Crossley contact model for normal forces and Coulomb friction for the tangential component. The {py:class}`PenaltyContact ` can handle collisions coming from multiple sources (multiple {py:class}`FrameCollider `s). #### ContactForceBase Every contact force calculation is done on classes derived from {py:class}`ContactForceBase `. Derived classes must override the {cpp:func}`computeForce ` method. This method receives the {py:class}`Karana.Collision.FrameContact` structure for the collision as well as two nodes (one for each body involved in the collision). The {cpp:func}`computeForce ` method uses this information to calculate the contact force (as a {py:class}`SpatialVector `) and returns it. #### Putting it all together The steps to add collision detection and contact forces to a simulation are as follows: 1. Create collision engines. Oftentimes, only one collision engine will be used, but there are scenarios where there may be multiple. 1. Create a {py:class}`FrameCollider ` for each collision engine. 1. Create a {py:class}`PenaltyContact ` with the {py:class}`FrameCollider `(s) and {py:class}`ContactForceBase ` as constructor arguments. 1. Set parameters for {py:class}`ContactForceBase ` class if appropriate. The [2 link pendulum (collisions)](generated/notebooks/example_pendulum_collision/notebook.ipynb) example notebook illustrates the use of collision dynamics with {{ kdflex }}. ## Simulation Core ### Prefabs {{ prefabscls }} are hierarchical, self‑contained, reusable components that can be assembled into a larger simulation. Typical use‑cases include: * A recurring subsystem, e.g., a thruster, helicopter rotors, a link in a robotic arm. * A “template” that can be instantiated many times with different parameters. * A logical grouping of bodies, joints, and associated {{ kmodelscls }}. --- #### Config, Context, and Params Each {{ prefabcls }} contains a {py:class}`Config `, {py:class}`Context `, and {py:class}`Params `. These are described in the table below. | Role | What it holds | Examples | |------|---------------|-----------------| | {py:class}`Config ` | Creation-time configuration data for the {{ prefabcls }}. This information defines items that should not change after construction, and defines the {{ prefabcls }}'s structure. These are also used to select between different "variants" of a {{ prefabcls }}.| For example, a thruster {{ prefabcls }} may have a few different thruster {{ kmodelscls }} that can be used, and which one a particular {{ prefabcls }} instance uses would be selected in the {py:class}`Config `. | | {py:class}`Context ` | Data that defines how a {{ prefabcls }} connects to the rest of the simulation. | The parent {{ physbodycls }} that the {{ prefabcls }} will connect {{ physbodiescls }} it creates to, the {{ mmcls }} to connect {{ kmodelscls }} it creates to | | {py:class}`Params ` | All other values to used to create the {{ prefabcls }}. | Mass properties, dimensions, etc. | #### Prefab initialization When a new Prefab instance is created, the following steps are executed in order: | Step | Method | Purpose | What to put in | |------|--------|---------|----------------| | 1 | {py:meth}`runBeforeInitialSetParent ` | Do any work that must be done before this {{ prefabcls }}'s parent is set. | Put in any custom logic that needs to run for the {{ prefabcls }} before its parent {{ prefabcls }} is set. | | 2 | {py:meth}`checkParentCompatible ` / {py:meth}`checkChildCompatible ` | When a {{ prefabcls }}'s parent {{ prefabcls }} is set, the parent {{ prefabcls }} runs {py:meth}`checkChildCompatible ` on the {{ prefabcls }} being added as a child, and the new child {{ prefabcls }} runs {py:meth}`checkParentCompatible ` on the new parent {{ prefabcls }}. These methods should return `True` if the parent/child is compatible, and throw an error if it is not. | Any logic to determine if the parent/child is compatible. These methods are designed to help restrict what the parent/children are allowed to be, if any restrictions exist. If the criteria are not met, then an error should be raised. These methods, if defined, will sometimes make use of {py:mod}`PrefabInterfaces `. See the {ref}`prefab_interfaces` section for details.| | 3 | {py:meth}`addChildPrefabs ` | Create any children {{ prefabscls }}. For example, a helicopter {{ prefabcls }} may create rotor {{ prefabscls }}. The initialization logic brings the new child prefabs up to the {py:meth}`addChildPrefabs ` stage. Meaning, all children {{ prefabscls }} will be added recursively before moving on to the next step, i.e., all {{ prefabscls }} will be created before any {{ prefabcls }} has multibody objects added in the next step.| Any children {{ prefabscls }} that should be created by this {{ prefabcls }}. | | 4 | {py:meth}`addMultibodyObjects ` | Create all multibody objects that belong to this {{ prefabcls }}. As in the previous step, all {{ prefabscls }} in the hierarchy will have {py:meth}`addMultibodyObjects ` called before moving on to the next step. | Calls to add multibody objects such as {{ physbodiescls }} via {py:meth}`addBody ` or {{ nodescls }} via {py:meth}`addNode `. | | 5 | {py:meth}`addKModels ` | Create all {{ kmodelscls }} that belong to this {{ prefabcls }}. As in the previous step, all {{ prefabscls }} in the hierarchy will have {py:meth}`addKModels ` called before moving on to the next step. | Logic to create any {{ kmodelscls }} that belong to this {{ prefabcls }}. | | 6 | {py:meth}`connectKModels ` | Connect {{ kmodelscls }} together. | Any logic needed to connect {{ kmodelscls }} together. This will often require finding other {{ kmodelscls }} in other {{ prefabscls }} in the hierarchy. See {ref}`prefab_hierarchy_request` for details on how to do this. | (prefab_hierarchy_request)= #### Requesting information from the Prefab hierarchy {{ prefabscls }} expose a set of helpers that allow a {{ prefabcls }} to query its ancestors or descendants for information. The {py:meth}`requestFromHierarchy ` method takes a filesystem-like path (can contain globs, wildcards, etc.), a callable, and an optional root. This method obtains any {{ prefabscls }} that match the filesystem-like path as referenced from the root, if defined, or the root of the entire {{ prefabcls }} hierarchy if not defined. Then, the callable is run on each of these matching {{ prefabscls }}. The callable should be of type `Callable[[Prefab[Any, Any, Any,]], T | None]`. In other words, it should take a {{ prefabcls }} as an argument, and return a type or None. Below is an example callable to get a {py:class}`GravityInterface ` if the {{ prefabcls }} implements the {py:class}`GravityProviderInterface ` (you can read more about {py:mod}`PrefabInterfaces ` in the {ref}`prefab_interfaces` section). ```python def getGrav(p: Prefab[Any, Any, Any]) -> GravityInterface | None: if isinstance(p, GravityProviderInterface): return p.getGravityInterface() else: return None ``` To use this to search across the entire {{ prefabcls }} hierarchy, one would make the call like this. ```python all_grav_interfaces = self.requestFromHierarchy("**", getGrav) ``` To request only gravity interfaces from children of this {{ prefabcls }}, one would make the call like this. ```python children_grav_interfaces = self.requestFromHierarchy("**", getGrav, self) ``` The other request methods follow the same pattern, but add some extra logic for convenience. {py:meth}`requestUniqueFromHierarchy ` follows the same process, but ensures that there is only one object returned after mapping the callable over all the matching {{ prefabscls }}. If there are no objects returned or more than one object returned, an error is thrown. {py:meth}`requestClosestFromHierarchy ` tracks each returned object's distance to the {{ prefabcls }} that made the request. Only the object closest is returned. If there is a tie for closest or no objects are returned, then an error is thrown. (prefab_interfaces)= #### Prefab interfaces {py:mod}`PrefabInterfaces ` provide a convenient way to ensure that a {{ prefabcls }} meets certain requirements. For example, the {py:class}`GravityProviderInterface ` ensures that {{ prefabcls }} implements a `getGravityInterface` method with the signature ```python def getGravityInterface(self) -> GravityInterface: ``` The {py:mod}`PrefabInterfaces ` are used with the `isinstance` method to ensure a {{ prefabcls }} meets that interface. For example: ```python # Let p be a Prefab, then if isinstance(p, GravityProviderInterface): # This Prefab meets the requirements of the GravityProviderInterface ... else: # This Prefab does not meet the requirements of the GravityProviderInterface ... ``` #### Adding bodies and nodes {{ prefabscls }} provide convenience methods for adding {{ physbodiescls }} and {{ nodescls }}. The {py:meth}`addBody ` method can be used to add a {{ physbodycls }}. It takes a parent {{ physbodycls }} where the new body will be attached and a child body. The child body can be either a {{ physbodycls }} or a {py:class}`BodyDS `. If it is a {{ physbodycls }}, then the method just ensures that it is indeed a child of the parent {{ physbodycls }}. If it is a {py:class}`BodyDS `, then it creates a new body based on the provided information and attaches it to the parent {{ physbodycls }}. The {py:meth}`addNode ` method functions similarly. It takes as arguments a parent {{ physbodycls }}, and the {{ nodecls }} as either a {{ nodecls }} instance or a {py:class}`NodeDS `. If it is a {{ nodecls }} instance, it ensures that the provided {{ nodecls }} is attached to the parent {{ physbodycls }}. If it is a {py:class}`NodeDS `, then it creates a new {{ nodecls }} and attaches it to the parent {{ physbodycls }}. #### Converting Prefabs to a DataStruct and vice-versa The {py:meth}`toDS ` method can be used to convert a {{ prefabcls }} to a {{ datastructcls }}. Oftentimes, this does not need to be overridden when creating a new {{ prefabcls }}, since it will automatically capture the derived {{ prefabcls }}'s {py:class}`Config `, {py:class}`Context `, and {py:class}`Params `, as well as its location in the {{ prefabcls }} hierarchy, which is typically all that is needed. When converting from a {{ datastructcls }} to a {{ prefabcls }}, the {py:meth}`fromDS ` method can be used. #### BasicPrefab The {{ basicprefabcls }} is a commonly-used {{ prefabcls }} when importing/exporting to/from various file types. For example, URDF files can be turned into a {py:class}`BasicPrefabDS `, and then into a {{ basicprefabcls }}. See the {py:class}`BasicPrefabDS `'s {py:meth}`fromFile ` method for more details on the file types supported. {{ basicprefabcls }} also has a {py:meth}`createStandalone ` method to create a standalone simulation from a {py:class}`BasicPrefabDS `. This will create and populate a {{ spcls }} and {{ mbodycls }} with the information on the {py:class}`BasicPrefabDS ` and create a {{ basicprefabcls }} from it. (datastruct_sec)= ### DataStruct The DataStruct Python class is built upon `pydantic` 's `BaseModel` and offers a useful way to define structured data in Python. It catches typing errors at runtime, which allows for early error detection and facilitates throwing better exceptions. This helps users efficiently fix problems. For fundamental usage that overlaps with `pydantic.BaseModel`, refer to the [pydantic documentation](https://docs.pydantic.dev/latest/). #### Type Checking and Validation `DataStruct` leverages `pydantic` 's robust type checking, ensuring that data conforms to defined types at creation and when fields are set. This allows for early error detection and informative exceptions. ##### Basic Validation One can define fields with specific types and add validation rules using `Annotated` and `pydantic.Field`: ```python from typing import Annotated from Karana.KUtils.DataStruct import DataStruct from pydantic import Field class A(DataStruct): """ Parameters ``````-- a : int Description for a. b : str b is a string. c : Annotated[int, Field(gt=0)] Positive integer. I'm a long description. d : float I'm a float. """ a: int b: str = "hello" c: Annotated[int, Field(gt=0)] # 'c' must be a positive integer d: float = 2.2 # This will work fine a = A(a=3, b="hello", c=4) # This will fail because 'c' must be greater than 0 # b = A(a=3, b="hi", c=-1) ``` Validation is also enforced when modifying fields after a `DataStruct` instance has been created. ```python a = A(a=3, b="hello", c=4) # Instance created successfully # This will fail as the value for 'c' is invalid # a.c = -1 ``` ##### Custom types and quantities `Karana.KUtils.Ktyping` defines convenient types that utilize the Python `pint` package. Here's an example using the `Mass` (a scalar value that must have units corresponding to mass) and `Length3` (3-vector with units that correspond to length) types: ```python from Karana.KUtils.DataStruct import DataStruct from Karana.Math.Kquantities import ureg from Karana.Math.Ktyping import Mass, Length3 import numpy as np class A(DataStruct): m: Mass l3: Length3 # These will work (can be floats or quantity objects with compatible units) a = A(m=3.0, l3=np.zeros(3)) a = A(m=3.0 * ureg.slug, l3=np.ones(3)) a = A(m=3.0 * ureg.slug, l3=np.ones(3) * ureg.mm) # This will fail due to incompatible units (ft for Mass) # a = A(m=3.0 * ureg.ft, l3=np.ones(3) * ureg.mm) ``` The system also provides custom unit-less types like `Vec3` (for 3-element `numpy` arrays) and `Mat33` (for 3x3 `numpy` arrays), which include shape and type validation. There are also convenient checks for use in Annotation, like `normCheck`, which validates against a specified norm value and tolerance: ```python import numpy as np from typing import Annotated from Karana.Math.Ktyping import Vec3, normCheck from Karana.KUtils.DataStruct import DataStruct class V3norm(DataStruct): v3: Annotated[Vec3, normCheck(2.5, 1e-5)] # Vec3 with norm ~2.5, tolerance 1e-5 # This passes because the norm (approx 2.5) is within tolerance a = V3norm(v3=np.array([0.0, 2.5 + 1e-6, 0.0])) # This fails because the norm (approx 2.5 + 1e-4) exceeds the tolerance # V3norm(v3=np.array([0.0, 2.5 + 1e-4, 0.0])) ``` #### Versioning The `DataStruct` class provides a built-in mechanism for tracking changes over time using a `version` field. This allows one to manage backward compatibility, automatically upgrade older versions if desired, and flag breaking changes. The `version` field is defined as a *2-integer tuple*. This tuple adheres to a common semantic versioning pattern to differentiate between types of changes. The *first number* in the tuple represents the *major version*. A change in the major version typically indicates *breaking changes* that are *not backward compatible*. When a major version mismatch occurs (e.g., trying to load a `1.x` version when the current is `2.x`), it usually necessitates an error, as automatic upgrades cannot be done without user intervention. For example, suppose there is a `DataStruct` with a new field added called `new_field` with no default. Further, suppose there is no good way to choose a default for this field for users. In this case, user intervention is required to choose the right value, and this is not backward compatible. The major version number should be incremented, and encountering older versions should trigger an error with a message or link to documentation that describes what the `new_field` is and how users might choose the right value and upgrade their `DataStruct`. The *second number* in the tuple represents the *minor version*. Changes to the minor version are *backward compatible*. This means that an older minor version (e.g., `1.1`) can often be automatically updated or migrated to a newer minor version (e.g., `1.2`) without data loss or functional issues. For example, suppose there is a `DataStruct` that gains a new field called `new_field` with a default value of 3.0. This is backwards compatible as older versions of this `DataStruct` can just use the new field with the new default. It is still a good idea to trigger a warning to users in this case letting them know this is happening. Furthermore, providing a way for them to upgrade, e.g., add this `new_field` and step up the version so they don't see the warning again, is helpful. This can often be accomplished by suggesting they use `toFile` to save the automatically upgraded version of the `DataStruct` to wherever they are currently loading it from. ##### Setting the current version The `DataStruct` class defines a *`_version_default`* class variable that holds the default version. To set the current, expected version for a class derived from `DataStruct`, override this `_version_default` class variable. ```python from Karana.KUtils.DataStruct import DataStruct from typing import ClassVar class MyData(DataStruct): # Set the current expected version for this data structure. # This overrides the (0, 0) default inherited from DataStruct. _version_default: ClassVar[tuple[int, int]] = (1, 0) # ... other fields ``` ##### Version validators To implement the logic for handling different versions (issuing warnings, raising errors, or performing automatic upgrades), you should use *Pydantic's `field_validator` or `model_validator`*. A `field_validator` can validate the `version` field itself and trigger warnings or errors. However, a `field_validator` *cannot modify other variables of the class*, so it cannot be used for automatic upgrades. Below is an example of a `field_validator` used for version checking and issuing warnings/errors: ```python from typing import ClassVar from pydantic import field_validator from Karana.Core import error, warn from Karana.KUtils.DataStruct import DataStruct class A(DataStruct): @field_validator("version", mode="after") @classmethod def _validateVersion(cls, version: tuple[int, int]) -> tuple[int, int]: """Validate the version. This cannot modify other variables of the class. Parameters ---------- version : tuple[int, int] The version. Returns ------- version : tuple[int, int] The validated version. """ if version == cls._version_default: # Early out for most common case where version matches the default. return version elif version < cls._version_default: # Trigger a warning for older, potentially compatible versions. warn( f"Got version {version} which is less than the current version {cls._version_default}." ) else: # Trigger an error for newer, potentially incompatible versions. raise ValueError(f"Got version {version} which is greater than the current version {cls._version_default}.") # Return the validated version return version # Example usage illustrating the field validator's behavior: a = A(version=(0, -1)) # This will trigger a warning b = A(version=(1, 0)) # This will trigger an error ``` For more complex logic, especially when there is a need to *upgrade other fields* based on a version mismatch, one should use a *`model_validator`*. A `model_validator` gets access to the full model instance (`self`), allowing it to modify other fields. This is ideal for implementing automatic backward-compatible upgrades. Below is an example of a `model_validator` that performs automatic upgrades: ```python from typing import Self, ClassVar from pydantic import model_validator from Karana.Core import error from Karana.KUtils.DataStruct import DataStruct class B(DataStruct): # Suppose this is a newly added field in 1.1 (from 1.0). f: int = 3 g: int # This class now sets its _version_default to (1, 1) overriding the (0, 0) default from DataStruct. _version_default: ClassVar[tuple[int, int]] = (1, 1) @model_validator(mode="after") def _validateVersion(self) -> Self: """Validate the version. This is a model validator, so it gets access to the full model. This can be used to upgrade other fields based on a version mismatch. Returns Self self after validation. """ if self.version == self._version_default: # Early out if the version matches the current default. return self if self.version < self._version_default: # This block handles backward-compatible upgrades (typically minor version changes) # or raises an error for major version downgrades. if self.version[0] < self._version_default[0]: # If the major version is older, this is considered a breaking change. raise ValueError("Major version change. Cannot handle! To upgrade please do .") else: # If only the minor version is older, an automatic, backward-compatible upgrade is possible. print("Automatic upgrade. Backward compatible. Add a message telling the user this has been done and what to do to avoid this warning in the future.") # Do the upgrade. Here our example upgrade sets g = f. self.g = self. if self.version > self._version_default: # Should not get a nonexistent future version. Trigger an error. raise ValueError(f"Got version {self.version}, which is a nonexistent future version.") return self ``` #### Equality and approximate equality (isApprox()) `DataStruct` instances support both strict equality checking with the usual `==` operator and approximate equality checking via the `isApprox` method. These methods recurse through the composite data types, e.g., a field that is a list of other `DataStruct` s. Consider the following example with nested `DataStruct` instances and `numpy` arrays: ```python from Karana.KUtils.DataStruct import DataStruct from Karana.KUtils.Ktyping import Vec3 from copy import deepcopy import numpy as np class A(DataStruct): a: int b: str = "hello" c: int class B(DataStruct): a: dict[str, list[list[A] | A]] a2: dict[str, dict[str, int | tuple[A, A] | Vec3]] a1 = A(a=3, b="hello", c=4) a2 = A(a=3, b="hello", c=4) b1 = B(a={"1": [[a1, a2], a1], "2": []}, a2={"3": {"3.1": 3, "3.2": (a1, a2), "3.3": np.zeros(3)}}) b2 = deepcopy(b1) b2.a2["3"]["3.3"] += 2e-8 # Introduce a small difference in a numpy array element # These two differ by 2e-8, so isApprox passes with prec=1e-7 assert b1.isApprox(b2, prec=1e-7) # These two still differ by 2e-8, so isApprox fails with prec=1e-8 assert not b1.isApprox(b2, prec=1e-8) # In their current state, b1 and b2 are not strictly equal assert b1 != b2 # After making b2 a fresh deepcopy of b1, they are strictly equal b2 = deepcopy(b1) assert b1 == b2 ``` #### Saving to and loading from files {{ datastruct }} provides convenient methods for saving to files, the {py:meth}`toFile ` method, and for loading the data back, the {py:meth}`fromFile ` method. The file type is determined by the file name's suffix. Currently, the supported formats are `json`, `yaml`, `pickle`, and `h5`. ##### Meta data When a {{ datastruct }} is saved to a file, meta data is saved with it. This is saved in the form of the {py:class}`Karana.KUtils.DataStruct.Meta` {{ datastruct }}. This data is provided via the `meta` keyword argument to the {py:meth}`toFile ` method. If this keyword argument is not supplied, then default {py:class}`Karana.KUtils.DataStruct.Meta` data will be generated automatically. One can add more information to the {py:class}`Karana.KUtils.DataStruct.Meta` class if desired by deriving a new class from {py:class}`Meta ` and adding fields said new class. When loading data using the {py:meth}`fromFile ` method, the {py:class}`Meta ` data can be returned in addition to the {{ datastruct }} being loaded by setting the `get_meta` keyword argument to `True`. See the {py:meth}`fromFile ` method for more details. ##### Custom class handling In order to save/load a `DataStruct`, all fields need to be serialized. In most cases, a `DataStruct` 's fields will be serializable without any changes. However, for custom classes, e.g., `Karana.Math.UnitQuaternion`, small modifications need to be made. Custom classes define certain methods in order to be serialized, and for `yaml` and `json`, must also be registered with the appropriate saver/loader. The classes in `kdFlex` that are often used, e.g., `UnitQuaternion` and `SpatialVector`, already have the methods defined. If saving/loading from a `pickle` or `h5` file, then nothing needs to be done, and if saving/loading from `yaml` or `json` all one must do is register them with the appropriate saver/loader. ###### YAML This example shows saving and loading a class that uses `UnitQuaternion` and `SpatialVector`. These classes already have `to_yaml` and `from_yaml` defined, so serialization and deserialization is handled automatically. ```python import numpy as np from Karana.KUtils.DataStruct import DataStruct from Karana.Math import UnitQuaternion, SpatialVector class YamlInOut(DataStruct): a1: UnitQuaternion = UnitQuaternion(0.5, 0.5, 0.5, 0.5) a2: SpatialVector = SpatialVector(np.array([1-3]), np.array([4, 6, 8])) # Example usage: dark = YamlInOut() dark.toFile("example.yaml") loaded_dark = YamlInOut.fromFile("example.yaml") assert loaded_dark == dark ``` For examples of how to write the `to_yaml` and `from_yaml` methods, see `Karana.Math._Math_Py.py`. ###### JSON This example shows saving and loading a class that uses `UnitQuaternion` and `SpatialVector`. These classes already have `to_json` and `from_json` defined, so serialization and deserialization is handled automatically. ```python import numpy as np from Karana.KUtils.DataStruct import DataStruct from Karana.Math import UnitQuaternion, SpatialVector class JsonInOut(DataStruct): a1: UnitQuaternion = UnitQuaternion(0.5, 0.5, 0.5, 0.5) a2: SpatialVector = SpatialVector(np.array([1-3]), np.array([4, 6, 8])) # Example usage: dark = JsonInOut() dark.toFile("example.json") loaded_dark = JsonInOut.fromFile("example.json") assert loaded_dark == dark ``` For examples of how to write the `to_json` and `from_json` methods, see `Karana.Math._Math_Py.py`. ###### Pickle and H5 To save/load from a `pickle` or `h5` file, the class need only be picklable. This means that its `__getstate__` and `__setstate__` must be defined. Nothing needs to be done on the `DataStruct` itself, the class just needs to define these methods. ##### File I/O Examples `DataStruct` comes pre-loaded with support for `numpy` arrays and `quantities.Quantity` s when saving to and loading from files. Custom classes can be registered as described above. It also automatically supports composite types, e.g., embedded `DataStruct` instances, lists of dictionaries of numpy arrays, etc. Here's an example demonstrating saving and loading an instance of `A` (a `DataStruct` containing `Mass` and `Length3` types) to a YAML file: ```python from pathlib import Path from Karana.KUtils.Ktyping import Mass, Length3 from Karana.KUtils.DataStruct import DataStruct import numpy as np from tempfile import NamedTemporaryFile class A(DataStruct): m: Mass l3: Length3 a = A(m=3.0, l3=np.zeros(3)) # Create an instance tf = NamedTemporaryFile(suffix=".yaml") # Create a temporary YAML file p = Path(tf.name) a.toFile(p) # Save the DataStruct to the file b = A.fromFile(p) # Load it back assert a == b # Verify equality ``` Additionally, `DataStruct` instances can be saved to and loaded from an `h5py` group within an HDF5 file. ```python import h5py from Karana.KUtils.DataStruct import DataStruct from tempfile import NamedTemporaryFile import numpy as np class A(DataStruct): x: float = 3.0 y: str = "hi" class B(DataStruct): a3: dict[str, A | int] # Create a sample B instance a1 = A(x=1.0, y="1") b = B(a3={"5": a1, "6": 3}) tf = NamedTemporaryFile(suffix=".h5") with h5py.File(tf.name, "r+") as f: g = f.create_group("new_group") # Create an HDF5 group b.toFile(g) # Save the DataStruct to the group b_new = B.fromFile(g) # Load it back from the group assert b == b_new ``` #### NestedBaseMixin The `NestedBaseMixin` provides a way to deserialize nested Pydantic classes, where the nested class is a base class and subclasses may be provided, e.g., `BodyDS` in `BodyWithContextDS` from `Karana.Dynamics.SOADyn_types`. The `NestedBaseMixin` should be added as a mixin to the base class (`BodyDS` in the example with `BodyWithContextDS`). The string provided to `NestedBaseMixin` will be the name of the field that is added to the base class and all of its derived classes by the mixin, and it is used to track the actual class type used. `pydantic` 's `SerializeAsAny` should always be used when the `NestedBaseMixin` is used, so that the correct type is serialized as well. `SerializeAsAny` is essentially duck-type serialization as opposed to strict-type serialization. Below is an example of how `NestedBaseMixin` is used: ```python from pydantic import SerializeAsAny from Karana.KUtils.DataStruct import DataStruct, NestedBaseMixin class Base(DataStruct, NestedBaseMixin("type_name")): a: int = 3 class D1(Base): b: int = 4 class D2(Base): c: int = 5 class Test(DataStruct): l: list[SerializeAsAny[Base]] = [] t = Test(l=[D1(b=0), D2(c=0), Base(a=1)]) t.toFile("test.yaml") assert Test.fromFile("test.yaml") == t ``` #### To/From DataStruct and Dictionary Converting `DataStruct` instances to and from standard Python dictionaries is straightforward. To convert a `DataStruct` instance to a dictionary, use the `model_dump` method: ```python my_dict = my_data_struct.model_dump() ``` To create a `DataStruct` instance from a dictionary, use the ``** dictionary unpacking operator: ```python my_data_struct = MyDataStruct(**my_dict) ``` #### IdMixin The {py:class}`Karana.KUtils.DataStruct.IdMixin` class is used to track objects that have been created from {{ datastructs }}. This associates the objects created via the {{ datastructs }} with the unique ID of the original object used to create the {{ datastruct }}. It is the responsibility of the class using the {py:class}`IdMixin ` to set the `_id` when a {{ datastruct }} is created from an object, and to call the {py:meth}`addObjectFromId ` when an object is created from a {{ datastruct }}. Here is a simple example: ```python from Karana.KUtils.DataStruct import DataStruct, IdMixin from Karana.Core import Base class C1(Base): def __init__(self): super().__init__("c1") self.x = 3.0 def toDS(self) -> "C1DS": c1_ds = C1DS(x=self.x) c1_ds._id = self.id() return c1_ds @staticmethod def fromDS(c1_ds: "C1DS") -> "C1": return c1_ds.toC1() class C1DS(DataStruct, IdMixin[C1]): x: float def toC1(self) -> C1: obj = C1() obj.x = self.x self.addObjectFromId(obj) return obj ``` In this example, the `toDS` method sets the `_id` member whenever a `C1DS` is created from an object, and the `toC1` method calls the {py:meth}`getObjectsCreatedById ` method whenever an object is created from a `C1DS`. The {py:class}`IdMixin ` class does the bookkeeping so that one can lookup any `C1` objects that were created from any `C1DS` with a specific `id`. For example, ```python class C2(Base): def __init__(self, c1: C1): super().__init__("c2") self.c1 = c1 self.y = 5.0 def toDS(self) -> "C2DS": c2_ds = C2DS(y=self.y, c1=self.c1.id()) c2_ds._id = self.id() return c2_ds @staticmethod def fromDS(c2_ds: "C2DS") -> "C2": return c2_ds.toC2() class C2DS(DataStruct, IdMixin[C2]): c1: int y: float def toC2(self) -> C2: """Use global lookup method.""" c1 = self.getObjectsCreatedById(self.c1) if c1 is None: raise ValueError("Cannot find C1 object") else: obj = C2(cast(C1, c1[0])) obj.y = self.y self.addObjectFromId(obj) return obj ``` The `toC2` method on `C2DS` looks for instances of `C1` that was created by any `C1DS` instance with a specific `id`. For example, they could be called like this: ```python # Code that loads c1_ds of type C1DS and c1_ds of type C2DS from a file c1 = c1_ds.toC1() c2 = c2_ds.toC2() assert c2.c1 is c1 ``` This was a simple standalone example. A more complex example using real simulation objects is shown in the {ref}`kmodel_datastructs_sec` section, where a KModel's {{ datastruct }} uses {py:class}`IdMixin ` to find {{ nodescls }} on {{ physbodiescls }} for use when creating an instance of the KModel from the {{ datastruct }}. #### Asciidoc Generation `DataStruct` instances can generate their own asciidoc documentation using the `generateAsciidoc` method. This method returns a string containing the asciidoc data, which can then be written directly to an asciidoc file. For example, if you have a `DataStruct` class named `A`, you can generate its documentation as follows: ```python # Let A be a DataStruct class and "a" be an instance of class A class A(DataStruct): """ Parameters ---------- a : int Description for a. b : str b is a string. c : Annotated[int, Field(gt=0)] Positive integer. I'm a long description. d : float I'm a float. """ a: int b: str = "hello" c: Annotated[int, Field(gt=0)] d: float = 2.2 a = A(a=3, b="hello", c=4) with open("A_doc.asciidoc", "w") as f: f.write(a.generateAsciidoc()) ``` (units_quantities_sec)= ### Units and Quantities The Kquantities module is designed to provide common quantities used throughout simulations and includes helper functions for unit conversion. It also offers functionalities to define and set a units system, with SI (International System of Units) as the default. This module acts as a simple wrapper around the Python `pint` package. #### What units and quantities are A quantity refers to a type of unit like length, force, etc. The unit is a specific instance of quantity, e.g., for the `length` quantity the units could be meters, feet, inches, etc. The `Kquantities` module defines several standard physical quantities as module-level variables. These are re-calculated internally to reflect the currently set units system. A few examples include: * length * mass * time * velocity #### check() and to_base_units() One can import these predefined quantities directly from `Karana.KUtils.Kquantities`. These imported quantities can be used to check if a value is of that quantity type using the `check` method. **Example of check() usage:** ```python from Karana.Math.Kquantities import velocity, torque, ureg v = 1.0 * ureg.m / ureg.s # Define a velocity t = 1.0 * ureg.N * ureg.m # Define a torque # Check if 'v' is a velocity quantity assert v.check(velocity) # True # Check if 't' is not a velocity quantity assert not t.check(velocity) # True # Check if 't' is a torque quantity assert t.check(torque) # True # Check if 'v' is not a torque quantity assert not v.check(torque) # True ``` All `pint` `Quantity`s also provides a `to_base_units` method to transform a given quantity into the units for the associated unit system. **Example of to_base_units usage:** ```python from Karana.Math.Kquantities import ureg import numpy as np # Define a velocity and torque v_si = 1.0 * ureg.mm / ureg.s m_si = 1.0 * ureg.lb # Convert them to the current (default SI) units system converted_v = v_si.to_base_units() converted_m = m_si.to_base_units() # In the default SI system assert converted_v.m == 0.001 # 1 mm/s = 0.001 m/s assert np.round(converted_m.m, 8) == 0.45359237 # 1 lb = 0.45359237 kg print(f"Converted velocity (SI): {converted_v}") print(f"Converted torque (SI): {converted_m}") ``` #### Changing the units system The units system can be modified using the `setUnitRegistry` function. This function allows for setting a units system using `pint`'s `UnitRegistry`. This should be done before `ureg` is ever imported from or used from `Kquantities`. **Example of setUnitRegistry usage and its effect on conversion:** ```python from pint import UnitRegistry from Karana.Math.Kquantities import setUnitRegistry, velocity, torque # Change the length unit to centimeters and mass to grams. # This should be done before importing ureg. setUnitRegistry(UnitRegistry(system="cgs")) from Karana.Math.Kquantities import ureg # Convert the original quantities again. The underlying units system has changed. v = 1.0 * ureg.m / ureg.s t = 1.0 * ureg.N * ureg.m v_converted_cm = v.to_base_units() t_converted_cm = t.to_base_units() # A velocity of 1 m/s becomes 100 cm/s. # A torque of 1 N*m (kg*m^2/s^2) becomes kg * m^2 / s^2 * (100 cm)^2 / m^2 * (1000 g) / kg = 10^7 kg * cm^2 / s^2 print(f"After setting length to centimeters:") print(f" Velocity: {v} converted to {v_converted_cm} (Expected: 1000.0)") print(f" Torque: {t} converted to {t_converted_cm} (Expected: 10,000,000.0)") assert v_converted_cm.m == 1e2 # 1 m/s = 1000 cm/s assert t_converted_cm.m == 1e7 # 1 N*m = 10^7 g*cm^2/s^2 ``` ### Var The {py:meth}`Var ` class (short for variable) can be thought of as a container that holds a value that can change over time. These containers are convenient, as they provide a common interface for logging and plotting. The value that the {py:meth}`Var ` holds is produced by a callable that is passed in at construction. This callable is invoked each time the {py:meth}`Var ` is read. In practice, think of a {py:meth}`Var ` as a live variable that can be queried at any simulation step. The callable provided can compute the value from the simulation state, read a sensor, simply return a constant, etc. All objects derived from {py:class}`Base ` have a {py:attr}`vars ` member. This member holds {py:meth}`Var `s of interest for the class. For example, the {py:attr}`mass ` {py:meth}`Var ` on {py:attr}`PhysicalBody.vars ` contains the current mass of the {{ physbody }}. {py:meth}`Var ` itself cannot be created directly. Rather, a typed version of {py:meth}`Var ` must be created. In C++, this is done via the templatized {cpp:class}`Karana::Var::Var_T` class. In Python, the common {py:meth}`Var ` types are wrapped with the type name given after the word "Var", e.g., {py:meth}`VarDouble ` for {py:meth}`Var `s with a `double` return type, or {py:meth}`VarVec ` for {py:meth}`Var `s that return a {cpp:type}`Karana::Math::Vec`. Variable sized {py:meth}`Var ` objects, e.g., for {cpp:type}`Karana::Math::Vec`, must always return the same-sized type. The size is specified as a constructor argument. See below for examples. #### Creating Var objects As mentioned above, all objects derived from {py:class}`Base ` have a {py:attr}`vars ` member. Hence, oftentimes, the {py:meth}`Var `s one is interested in have already been created and exist on the {py:attr}`vars ` member of some class. However, there are times when this is not the case and one wants to create a custom {py:meth}`Var `. In C++, a {cpp:class}`Karana::Var::Var_T` can be created as follows: ```cpp #include auto my_var = Karana::Var::Var_T::create( "temperature", []() { return 298.15; } // function that returns a constant value ); auto my_var_2 = Karana::Var::Var_T::create( "my_vec", []() { km::Vec dark(5); dark << 0, 1, 2, 3, 4; }, // function that returns a constant vector of size 5 5 ); ``` In Python, the process is similar: ```python import numpy as np from Karana.KUtils import VarDouble, VarVec # Constant, scalar value my_var = VarDouble("temperature", lambda: 298.15) # Constant, vector value my_var_2 = VarVec("my_vec", lambda: np.linspace(0,4,5), 5) ``` #### Retrieving the value In both C++ and Python, the `()` operator can be used to retrieve the value. For example, in C++ ```cpp double my_val = my_var() ``` and in Python ```python my_val = my_var() ``` #### Using Vars for logging and plotting {py:meth}`Var `s can be used for both logging and plotting. For logging, the {py:meth}`PacketTableConfig.addData ` methods have overloads that support {py:meth}`Var `s. Similarly, for plotting, the {py:class}`SinglePlotData ` constructors allow {py:meth}`Var `s. For more details on how the {py:class}`PacketTableConfig ` class is used for logging see the {ref}`data_logging_sec` section. For more details on how the {py:class}`SinglePlotData ` is used for plotting, see the {ref}`recipe_plotting_sec` recipe. ### Logging The logger in `KCore` is used to display messages to the user. Each message comes with a level called the verbosity level. The user can set verbosity level of the logger to hide messages they don't want to see. The verbosity levels, in order, are: `off`, `error`, `warn`, `info`, `debug`, and `trace`. If a logger has a verbosity level set to `info`, but a message is being sent with verbosity `debug`, then the logger won't even print the message. The `Logger` class in `KCore` has the ability to send messages to multiple loggers, e.g., the terminal output and to a file. Each logger can have its own verbosity. The verbosity levels can be modified by using `Logger.changeVerbosity`. By default, there is only one logger, the logger to the terminal. The terminal logger special as it is actually comprised of two loggers: one to `stderr` and one to `stdout`. These use some custom verbosity settings, which are covered below. The names of these loggers are `stderr` and `stdout`. The default logging level of `stdout` can be set by the environment variable `KARANA_LOG_LEVEL`. The value of this variable should be the verbosity level desired in all capital letters, e.g., `KARANA_LOG_LEVEL=INFO` for `info`. If this variable is not set or does not match any known level, then `warn` is used. A file logger can be created using `Logger.addFileLogger`. The name of the logger will be the name of the file. #### stdout/stderr We have separate loggers for `stdout` and `stderr`. These loggers are custom, they do not follow the normal verbosity rules. That is because we want to treat them as if they were one logger. We want everything `error` and above to be printed to `stderr` and we want everything `warn` and below to be printed to `stdout`. Of course, one can still set the verbosity levels of these loggers. For example, setting the `stdout` verbosity to `error`, means it will not print anything, because it can't print anything above a warn, but its verbosity level is `error`, so it can't print anything below an `error`. #### Minimizing work We want to minimize the work we do if messages will not be consumed. This issue is aimed mostly at `debug` and `trace` messages, which can be time-consuming to create and print. For example, ```cpp // Work starts here to construct the debug string here double x = my_complicated_method() std::string y = my_time_consuming_method(); std::string debug_msg = std::format("debug string for double {}. See {}.", x, y); debug(debug_msg); ``` If no sinks are consuming said messages, we don't even want to do the work to form them, i.e., we want to avoid doing the work to create `x` and `y` in the example above. As a result, the `error`, `warn`, `info`, `debug`, and `trace` functions are all overloaded to accept a `std::function` (or `Callable[[], str]` in Python). The `std::function` version first checks the minimum level of all the sinks, if no sink is going to be consuming the message, we don't even bother to run the function. If at least one sink will consume the result, then we run the `std::function` to generate the string and pass it on to the loggers like usual. Using the example above, this would look like ```cpp debug([&](){ double x = my_complicated_method() std::string y = my_time_consuming_method(); std::string debug_msg = std::format("debug string for double {}. See {}", x, y); return debug_msg; }) ``` In Python, this would look like: ```python def debug_msg(): x = my_complicated_method() y = my_time_consuming_method() return f"debug string for float {x}. See {y}" debug(debug_msg) ``` If you are writing messages in C++, and that message is just a format string, you can use the format version of these messages. For example, ```cpp error("My error with a number '{}'", 3.0) ``` as opposed to ```cpp std::string str = std::format("My error with a number '{}'", 3.0); error(str) ``` #### Deprecated There is a `deprecated` function that can be used to print deprecation messages. This prints messages at the `warn` verbosity level. It takes in extra arguments * `name` - The name of the thing being deprecated. * `date` - The date that `name` will be removed. in addition to the normal message, which should specify how users should update the thing being deprecated. In C++, things should be deprecated by using the `[[deprecated]]` attribute. This will give users a message at compile time letting them know they are using a deprecated function. If the C++ function is used by `pybind11` to expose it to Python, then the `deprecatedWrapper` from `Karana/KCore/pybind11Utils.h` should be used to wrap the method. This will print a deprecation warning on the Python side, with a Python stack trace, whenever that function is used. In Python, the same `deprecated` function exists. However, it can also be used as a decorator like this: ```python @deprecated(date(2025, 4, 25), "old_f is deprecated. Use new_f instead") def old_f(): ... ``` The `name` parameter is not required when using the decorator version, since the `name` is calculated automatically based on the name of the function, class, etc. being deprecated. Both the `deprecated` function and `decorator` will automatically append the `deprecated` message with a stack trace so users can easily update their code. (DebugManager_sec)= ### DebugManager The {py:class}`Karana.Core.DebugManager` class holds static variables used to control extra debugging behavior. These are: * `all_destroyed_verbosity` - Controls the number of objects printed when calling all destroyed, if there are undestroyed objects. A value of `nullopt` in C++ or `None` in Python will print all undestroyed objects. * `enable_locking_trace` - This will enable the `LockingBase` trace output. Enabling this will print messages whenever a `LockingBase` is made healthy or not healthy. It will also print messages when new dependents are added. This affects all `LockingBase` s. * `enable_usage_tracking_map_errors` - Setting this to false will disable the usage tracking map errors. Use this with caution, as disabling this can allow memory leaks that would otherwise be caught. This should only be used in cases where the user is stuck and cannot move forward without assistance, but cannot wait for assistance, i.e., this is a sort of escape hatch to ignore cleanup related issues, but should only be used temporarily. (data_logging_sec)= ### Data logging The `DataLogger` class provides a convenient way to record data into HDF5 files, leveraging packet tables to efficiently append data to the files throughout the simulation. This system is built around two primary components: the `PacketTableConfig` and the `H5Writer`. (data_logger_setup_sec)= #### Setting up your data logger The process of logging data involves these steps: 1. **Define your data structure** using `PacketTableConfig` instances for each table you wish to create. 2. **Initialize the logger** by creating an `H5Writer` instance. 3. **Register your data configurations** by calling `createTable` on the `H5Writer` for each `PacketTableConfig`. 4. **Log your data** by invoking the `log` method on your `H5Writer`. (packet_table_config_setup_sec)= #### Defining data with PacketTableConfig The `PacketTableConfig` class is used to register functions that will produce data for the columns of your HDF5 packet table. Each `PacketTableConfig` represents a single table in your HDF5 file. You instantiate it by providing a name, which can include forward slashes (`/`) to specify a group hierarchy within the HDF5 file; groups will be created as needed. **Example: Creating a `PacketTableConfig`** **C++:** ```cpp auto c = Karana::KUtils::PacketTableConfig::create("data"); // Or for nested groups auto c2 = Karana::KUtils::PacketTableConfig::create("mygroup/inner_group/data"); ``` **Python:** ```python from Karana.KUtils import PacketTableConfig c = PacketTableConfig("data") c2 = PacketTableConfig("mygroup/inner_group/data") ``` Once a `PacketTableConfig` is made healthy, its configuration cannot be changed, as the types become fixed at that time for the table. Attempting to modify it afterward will result in an error. ##### Adding data columns (the addData() methods) To specify what data goes into each column of your table, you use the `addData` methods of `PacketTableConfig`. These methods generally require two main arguments: * A `name` for the column. * A callable function (`f`) that takes no arguments and returns the data for that column. The `addData` methods handle various data types, including scalars, dynamic vectors/matrices, and fixed-size Eigen types. For basic scalar types (like `int`, `double`), the C++ `addData` method automatically deduces the return type of the provided function. **Example: Adding basic scalar data** **C++:** ```cpp int x = 0; double y = 3.2; double t = 0.0; c->addData("time", [&t]() { // Returns double t += 0.1; return t; }); c->addData("int", [&x]() { // Returns int x += 2; return x; }); c->addData("double", [&y]() { // Returns double y += 0.6; return y; }); ``` In Python, due to dynamic typing, separate `addData` methods are provided for different return types (e.g., `addDataInt`, `addDataFloat`). **Python:** ```python t = 0.0 def time(): global t t += 0.1 return t c.addDataFloat("time", time) # Use addDataFloat for float return type i = 1 def int_f() -> int: global i i += 1 return i c.addDataInt("int", int_f) # Use addDataInt for integer return type ``` ###### Dynamic vectors and matrices For dynamic vectors and matrices (e.g., {cpp:type}`Karana::Math::Vec`, {cpp:type}`Karana::Math::Mat` in C++ or NumPy arrays in Python), you must provide extra integer arguments to specify the **size** of the vector or the **rows and columns** of the matrix. **Example: Adding dynamic vector/matrix data** **C++:** ```cpp km::Vec v(2); // Dynamic vector of size 2 v.setZero(); km::Mat m(2, 3); // Dynamic matrix of size 2x3 m.setZero(); c->addData( "vec", [&v]() { v[0] += 1.0; v[1] += 2.0; return v; }, 2); // Specify vector_size = 2 c->addData( "mat", [&m]() { m(0, 0) += 1.0; return m; }, 2, // Specify matrix_rows = 2 3); // Specify matrix_cols = 3 ``` **Python:** ```python import numpy as np v = np.zeros(2) # Dynamic NumPy array (vector) def vec(): global v v[0] += 1.0 v[1] += 2.0 return v c.addData("vec", vec, 2) # Specify vector size = 2 m = np.zeros((2, 3)) # Dynamic NumPy array (matrix) def mat(): global m m[0,0] += 1.0 return m c.addData("mat", mat, 2, 3) # Specify matrix_rows = 3, matrix_cols = 2 ``` ###### Fixed-size eigen types For **fixed-size Eigen vectors and matrices** in C++, special `addData` overloads exist that infer dimensions directly from the Eigen type, so you don't need to provide explicit size arguments. However, these require the Eigen type to have a fixed size; attempting to use a dynamically sized Eigen type without specifying dimensions will trigger a static assertion error. **Example: Adding fixed-size eigen data** **C++:** ```cpp km::Vec3 v3{1.0, 2.0, 3.0}; // Fixed-size 3D vector km::Mat33 m3; // Fixed-size 3x3 matrix m3.setZero(); c->addData("vec3", [&v3]() { v3[0] += 0.6; v3[1] += 1.0; return v3; }); // Dimensions inferred from km::Vec3 c->addData("mat33", [&m3]() { m3(0, 0) += 0.6; return m3; }); // Dimensions inferred from km::Mat33 ``` ###### Logging vectors/matrices as scalars (the as_scalars option) For both dynamic and fixed-size vectors and matrices, an optional `as_scalars` boolean argument can be set to `true`. When `as_scalars` is true, each component of the vector or matrix is logged as a separate scalar column in the packet table, rather than the entire vector/matrix being stored in a single column. The column names are appended with `_X` for vector components (where `X` is the index) or `_X_Y` for matrix components (where `X` is the row index and `Y` is the column index). **Example: Logging as scalars** **C++:** ```cpp // Using fixed-size Eigen types as scalars c->addData( "vec3s", [&v3]() { v3[0] += 0.6; v3[1] += 1.0; return v3; }, true); // Logs as "vec3s_0", "vec3s_1", "vec3s_2" c->addData( "mat33s", [&m3]() { m3(0, 0) += 0.6; return m3; }, true); // Logs as "mat33s_0_0", "mat33s_0_1", etc. // Using dynamic types as scalars c->addData( "vecs", [&v]() { v[0] += 1.0; v[1] += 2.0; return v; }, 2, // Vector size true); // Logs as "vecs_0", "vecs_1" c->addData( "mats", [&m]() { m(0, 0) += 1.0; return m; }, 2, // Matrix rows 3, // Matrix cols true); // Logs as "mats_0_0", "mats_0_1", etc. ``` **Python:** ```python # Vector as scalars def vec(): global v v[0] += 1.0 v[1] += 2.0 return v c.addData("vec", vec, 2, as_scalars=True) # Logs as "vec_0", "vec_1" # Matrix as scalars m = np.zeros((2, 3)) def mat(): global m m[0,0] += 1.0 return m c.addData("mat", mat, 2, 3, as_scalars=True) # Logs as "mat_0_0", "mat_0_1", etc. ``` #### Logging data with H5Writer The `H5Writer` class is responsible for managing your HDF5 file and its packet tables. It handles the creation of tables based on the `PacketTableConfig`s and facilitates the logging of data. **Example: Creating an `H5Writer`** You create an `H5Writer` by providing the desired filename for your HDF5 log file. **C++:** ```cpp #include "Karana/KUtils/DataLogger.h" auto h5 = Karana::KUtils::H5Writer::create("example.h5"); ``` **Python:** ```python from Karana.KUtils import H5Writer h5 = H5Writer("example.h5") ``` ##### Creating and logging tables Once your `H5Writer` is initialized, you need to register your `PacketTableConfig` instances with it using the `createTable` method. This makes the tables "active" by default, meaning they will be logged when `log()` is called. **Example: Creating and logging tables** **C++:** ```cpp h5->createTable(c1); // Register the packet table config h5->log(); // Log data for all active tables ``` **Python:** ```python h5.createTable(c) # Register the packet table config h5.log() # Log data for all active tables ``` The `log()` method will iterate through all **active** tables and write their contents into the log file. You can also log a specific table, regardless of its active status, using `logTable`. **Example: Logging a specific table** **C++:** ```cpp h5->logTable("mygroup/data") ``` **Python:** ```python h5.logTable("mygroup/data") ``` ##### Managing active tables The `H5Writer` also provides methods to manage which tables are actively logged: * **`activateTable(name)`**: Activates a table by its name. * **`deactivateTable(name)`**: Deactivates a table by its name, preventing it from being logged by `log()`. * **`getActiveTableNames()`**: Returns a vector/list of names of all currently active tables. * **`getTableNames()`**: Returns a vector/list of names of all registered tables (active or inactive). **Example: Managing active tables** **C++:** ```cpp h5->deactivateTable("mygroup/inner_group/data") h5->log() // Only active tables will be logged h5->deactivateTable("mygroup/data") h5->logTable("mygroup/data") // Log this specific table even if inactive h5->activateTable("mygroup/data") h5.log() // Now "mygroup/data" will be logged again by log() ``` **Python:** ```python h5.deactivateTable("mygroup/inner_group/data") h5.log() # Only active tables will be logged h5.deactivateTable("mygroup/data") h5.logTable("mygroup/data") # Log this specific table even if inactive h5.activateTable("mygroup/data") h5.log() # Now "mygroup/data" will be logged again by log() ``` #### DataLogger model During a simulation, we oftentimes want data to be logged regularly. The `DataLogger` from `GeneralKModels` is designed to do just that. **Example: DataLogger model** **C++:** ```cpp dl = Karana::Models::DataLogger::create("wsm", h5_writer) dl->params.log_first_step = True dl->setPeriod(0.1) // Must register this with the StatePropagator for it to have any effect sp.registerModel(dl) ``` **Python:** ```python # Create the DataLogger model dl = DataLogger("wsm", h5_writer) dl.params.log_first_step = True dl.setPeriod(0.1) # Must register this with the StatePropagator for it to have any effect sp.registerModel(dl) ``` The example above shows creating a `DataLogger` model that will call `log` on the `h5_writer` (this is an `H5Writer` instance) every 0.1 seconds. The `log_first_step` parameter can be used to also create a log entry at the current time, e.g., when starting a simulation, this will log the data at 0 seconds, rather than the first data entry being at 0.1 seconds. #### Examples Further examples can be found below. All filepaths are given relative to the installation directory, `/usr` for Linux and `/opt/homebrew/` for macOS. * [Example data logging Python notebook](generated/notebooks/example_data_logging/notebook.ipynb) - Simple data logging example done in Python * Example data logging C++ - An example is included with your installation at `share/Karana/Tutorials/example_cc/example_data_logging` that shows the same data logging, but done entirely in C++. * Example data logging C++, setup in Python - An example is included with your installation at `/share/Karana/Tutorials/example_cc/example_data_logging_pybind11` that shows you can set up a simulation and DataLogging in Python, but have the data logging functions themselves be entirely in C++. This is the best-of-both-worlds example, where you still have the runtime performance as is everything was done in C++, but you can set things up in Python.