\documentclass[letterpaper, 10 pt, conference]{ieeeconf}
\IEEEoverridecommandlockouts  

\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{mathrsfs}
\usepackage{graphicx}


\title{\LARGE \bf
Inertia Partitioning Modular Control Framework for Reconfigurable Multibody Systems
}

\author{Mohammad Dastranj$^{1,*}$ and Jouni Mattila$^{1}$% <-this % stops a space
\thanks{\protect\rule{0.965\linewidth}{0.4pt} \indent$^{1}$The authors are with the Unit of Automation Technology and Mechanical Engineering, Faculty of Engineering and Natural Sciences, Tampere University, 33720 Tampere, Finland
        {\tt\small \{mohammad.dastranj, jouni.mattila\}@tuni.fi}}%
}

\begin{document}

\maketitle 

\begin{abstract}
    A novel modular control framework for reconfigurable rigid multibody systems is proposed, motivated by the challenges of modular control of systems with closed kinematic chains. In the framework, modularity is defined in the sense of degrees of freedom, and the inertial properties of each body are partitioned with respect to how they are reflected in the kinetic energy of the system through the motion induced by each degree of freedom. This approach inherently handles closed chains in the same manner as tree-like structures, eliminating the need for explicit constraint force calculations or formulations based on differential-algebraic equations. The proposed framework is implemented via simulation on a three-degree-of-freedom series-parallel manipulator, with the results being consistent with the expected stability and tracking performance, and indicating the framework's potential for scalability in trajectory-tracking control of multibody systems.
\end{abstract}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
Advances in the design of mechanical systems have led to greater complexity in their modeling and control. With more complex systems, monolithic control approaches tend to have higher computational complexity. The surge in computational cost imposes limits for the implementation of monolithic controllers regarding scalability and real-time performance \cite{tinoco2025review}. Therefore, the motivation arises to adopt modular control approaches to control complex dynamic systems.

\subsection{Background}
There are well-established modular control approaches in the literature, namely passivity-based control (PBC) and virtual decomposition control (VDC) frameworks. PBC leverages passivity theory to design stable modular controllers by exploiting the passivity of individual subsystems while the interconnections preserve passivity \cite{ortega2002interconnection}. The port-Hamiltonian framework \cite{vanDerSchaft2017} provides a natural and rigorous foundation for controlling physical systems by using passivity theory and energy-based methods, as this framework provides an intrinsically passive and modular representation of physical systems, including rigid \cite{berger2025port} and flexible \cite{brugnoli2021port} multibody systems.

Alternatively, VDC breaks the complex multibody system down to simpler subsystems using virtual cutting points and models the interconnection of subsystems with virtual power flows. Then, the local controllers are designed, and the concept of virtual stability is utilized to provide the stability analysis of the controller \cite{zhu2010virtual}. 

\subsection{Literature Review}
Advances in heavy-duty robotic systems, such as those in construction, have introduced architectures with extended reach and multi-ton payload capacities \cite{KOIVUMAKI201559}. These systems often employ hybrid series–parallel structures driven by hydraulic or electro-mechanical actuators. However, the inclusion of closed kinematic chains in such mechanisms introduces additional challenges to model-based control.

Regarding the port-Hamiltonian framework, this issue is addressed in \cite{vanderschaft2018generalizedporthamiltoniandaesystems} and \cite{van2020dirac}, and the result is a system of differential-algebraic equations (DAEs). The addition of constraint equations to form the system of DAEs imposes numerical challenges and stability issues in forward dynamics \cite{fernandez2022non}, including constraint drift and the need for constraint stabilization \cite{khoshnazar2024application}. Furthermore, the DAE formulation increases the computational and analytical complexity of model-based control design and stability analysis \cite{drucker2023trajectory}.

VDC also relies on modifications to its general framework to handle closed chains, as it was originally developed for multibody systems with tree-like structure by using the spatial-vector form of the Newton-Euler formalism \cite{zhu2010virtual}. Currently, no widely adopted solution exists within the VDC framework for calculating the constraint forces and torques in closed chains, and the available modifications are often system-specific, depending on the particular structure of the closed chain. For instance, \cite{petrovic2022mathematical,ding2023high,barjini2025surrogate}, and \cite{zhang2025equivalence} have provided solutions for utilizing VDC for a certain class of multibody systems including triangular closed chains comprised of several passive revolute joints and one hydraulically or electrically actuated prismatic joint, which are common in heavy-duty robotic manipulators.

These frameworks, despite the solid foundations they provide for modular control, face challenges with closed chains. These challenges motivate the development of a unified modular control framework capable of inherently handling closed chains in rigid multibody systems, while avoiding the computational burdens commonly associated with DAE formulations.

\subsection{Contributions}
While the term \textit{module} in the related literature typically indicates a physical component or a small set of physical components that is a building block of a larger system, the term is considered differently as the basis of our new approach. In this paper, the Lagrangian formalism in minimal coordinates is used, and each generalized coordinate of the minimal set is considered a module.

The Jacobian matrices of body-fixed frames of reference are central to this framework. They are necessary to calculate the generalized inertia matrices, which capture how the inertial properties of each body contribute to the kinetic energy in the minimal coordinates, encapsulating the interconnection between system components. Also, the Jacobians are used to calculate the generalized forces associated with external wrenches exerted on the system. Modularity in this sense allows the design of controllers individually for each input and facilitates controller design for multi-input systems, and the use of minimal coordinates yields a system of ordinary differential equations (ODEs) \cite{shabana2009computational} instead of DAEs. The \textit{inertia partitioning} framework alleviates the computational issues inherent to DAEs as well as the challenges in using the DAE form in model-based control and stability analysis. This new approach is also well suited for reconfigurable multibody systems, where components can be added, removed, or modified. Local generalized inertia matrices allow for handling the reconfiguration without requiring modeling from scratch. Also, the framework enables augmenting the inertial effects of actuators within the dynamic modeling and controller design.

\subsection{Paper Structure}
The remainder of this paper is structured as follows. Section \ref{sec:dyn} explains the derivation of the Jacobians, the generalized inertia matrices, and the equations of motion with respect to each degree of freedom. Also, it explains how the modular aspect of the dynamic modeling allows for reconfiguration in the multibody system. The procedure for designing a controller based on the modular dynamic model of the multibody system and the stability analysis required for the controller is discussed in Section \ref{sec:control}. Section \ref{sec:simulation} models and controls a three-degree-of-freedom (3-DoF) series-parallel manipulator in a simulation by implementing the inertia partitioning framework, and the results of the simulation are presented. Finally, Section \ref{sec:conclusion} concludes the steps taken in the paper and the results that it yields with respect to the intended goals of the research.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Modular Dynamic Modeling of Multibody Systems}
\label{sec:dyn}
In the first part of this section, kinematic analysis is performed using the methodology in \cite{lynch2017modern}, whose results are required in the next part, where it is explained how the local generalized inertia matrix is obtained for each rigid body by using the spatial Jacobian corresponding to the motion of an arbitrary body-fixed frame of reference on the same body. Later, the equations of motion for each degree of freedom are obtained.

\subsection{Kinematics of the Multibody System}
For two frames of reference $\{\mathcal{A}\}$ and $\{\mathcal{B}\}$, the homogeneous transformation matrix comprised of the rotation matrix $\boldsymbol{R}_{\mathcal{A}\mathcal{B}}$ and the distance vector $\boldsymbol{p}_{\mathcal{A}\mathcal{B}}$ of frame $\{\mathcal{B}\}$ with respect to frame $\{\mathcal{A}\}$ and expressed in frame $\{\mathcal{B}\}$ is formed as
\begin{equation}
    \boldsymbol{\mathcal{T}}_{\mathcal{A}\mathcal{B}}=\begin{bmatrix}
    \boldsymbol{R}_{\mathcal{A}\mathcal{B}} &\boldsymbol{p}_{\mathcal{A}\mathcal{B}}\\
    \boldsymbol{0} & 1
    \end{bmatrix} \in SE(3)
    \label{eq,Tmat}
\end{equation}
where $SE(3)$ is the Special Euclidean group. In case of $f$ intermediate frames $\{\mathcal{F}_1\}, \, \ldots ,\,\{\mathcal{F}_f\}$ related by $\boldsymbol{\mathcal{T}}_{\mathcal{A}\mathcal{F}_1}, \, \boldsymbol{\mathcal{T}}_{\mathcal{F}_1\mathcal{F}_2}, \, \ldots \,,\,\boldsymbol{\mathcal{T}}_{\mathcal{F}_f\mathcal{B}}$ to each other, and the frames $\{\mathcal{A}\}$ and $\{\mathcal{B}\}$, the relation
\begin{equation}
    \boldsymbol{\mathcal{T}}_{\mathcal{A}\mathcal{B}}=\boldsymbol{\mathcal{T}}_{\mathcal{A}\mathcal{F}_1}\boldsymbol{\mathcal{T}}_{\mathcal{F}_1\mathcal{F}_2}\ldots\boldsymbol{\mathcal{T}}_{\mathcal{F}_f\mathcal{B}}
    \label{eq,TmatRec}
\end{equation}
holds \cite{lynch2017modern}. With $\{\mathcal{S}\}$ denoting an inertial frame and $\{\mathcal{B}\}$ indicating a body-fixed frame, the linear velocity $\boldsymbol{v}_\mathcal{B}$ and the angular velocity $\boldsymbol{\omega}_\mathcal{B}$ of frame $\{\mathcal{B}\}$ measured from the inertial frame $\{\mathcal{S}\}$ but expressed in frame $\{\mathcal{B}\}$ can be obtained using the inverse and the first time derivative of the transformation matrix $\boldsymbol{\mathcal{T}}_{\mathcal{S}\mathcal{B}}$, \cite{lynch2017modern}
\begin{equation}
    \begin{bmatrix}
        [\boldsymbol{\omega}_\mathcal{B}] & \boldsymbol{v}_\mathcal{B} \\ \boldsymbol{0} & 0
    \end{bmatrix}=\boldsymbol{\mathcal{T}}_{\mathcal{S}\mathcal{B}}^{-1} \boldsymbol{\dot{\mathcal{T}}}_{\mathcal{S}\mathcal{B}},
    \label{eq,bodytwist}
\end{equation}
considering that the notation $\left[.\right]$ denotes the skew-symmetric matrix representation of a vector. Assuming that the transformation matrix is written as a function of the minimal set of generalized coordinates $\boldsymbol{q}(t)\in \mathbb{R}^n$ and using (\ref{eq,bodytwist}), the linear Jacobian $\boldsymbol{\mathcal{J}}_\mathcal{B}$ and angular Jacobian $\boldsymbol{\mathcal{L}}_\mathcal{B}$ of frame $\{\mathcal{B}\}$ are calculated as 
\begin{equation}
\begin{split}
    \boldsymbol{v}_\mathcal{B}&=\boldsymbol{\mathcal{J}}_\mathcal{B}(\boldsymbol{q})\,\boldsymbol{\dot{q}} \\
    \boldsymbol{\omega}_\mathcal{B}&=\boldsymbol{\mathcal{L}}_\mathcal{B}(\boldsymbol{q})\,\boldsymbol{\dot{q}}.
    \label{eq,Jacob}
\end{split}
\end{equation}
Combining the linear and angular velocities in (\ref{eq,Jacob}) forms the twist vector
\begin{equation}
    \boldsymbol{\mathcal{V}}_\mathcal{B}=\begin{bmatrix}
        \boldsymbol{\omega}_\mathcal{B} \\ \boldsymbol{v}_\mathcal{B}
    \end{bmatrix}\in se(3)
    \label{eq,twist}
\end{equation}
with $se(3)$ being the Lie algebra associated with $SE(3)$. Substituting (\ref{eq,Jacob}) into (\ref{eq,twist}) yields the spatial Jacobian
\begin{equation}
    \boldsymbol{J}_\mathcal{B}=\begin{bmatrix}
        \boldsymbol{\mathcal{L}}_\mathcal{B}\\ \boldsymbol{\mathcal{J}}_\mathcal{B}
    \end{bmatrix}
    \label{eq,spatialJacobian}
\end{equation}
to be used as
\begin{equation}
    \boldsymbol{\mathcal{V}}_\mathcal{B}=\boldsymbol{J}_\mathcal{B}\boldsymbol{\dot{q}}.
    \label{eq,twistJacobian}
\end{equation}
To find the twist of frame $\{\mathcal{A}\}$ from the known twist of frame $\{\mathcal{B}\}$, the relation \cite{lynch2017modern}
\begin{equation}
    \boldsymbol{\mathcal{V}}_\mathcal{A}=\boldsymbol{Ad}_{\boldsymbol{\mathcal{T}}_{\mathcal{A}\mathcal{B}}}\boldsymbol{\mathcal{V}}_\mathcal{B}
    \label{eq,twistrelation}
\end{equation}
is used. $\boldsymbol{Ad}_{\boldsymbol{\mathcal{T}}_{\mathcal{A}\mathcal{B}}}$ is the adjoint representation of the transform matrix $\boldsymbol{\mathcal{T}}_{\mathcal{A}\mathcal{B}}$, and it is defined as
\begin{equation}
    \boldsymbol{Ad}_{\boldsymbol{\mathcal{T}}_{\mathcal{A}\mathcal{B}}}=\begin{bmatrix}
        \boldsymbol{R}_{\mathcal{A}\mathcal{B}} & \boldsymbol{0} \\ \left[\boldsymbol{p}_{\mathcal{A}\mathcal{B}}\right] \boldsymbol{R}_{\mathcal{A}\mathcal{B}} & \boldsymbol{R}_{\mathcal{A}\mathcal{B}}
    \end{bmatrix}.
    \label{eq,adjoint}
\end{equation}
From the twist relation (\ref{eq,twistrelation}) and the spatial Jacobian (\ref{eq,spatialJacobian}), the following is obtained:
\begin{equation}
    \boldsymbol{J}_\mathcal{A}=\boldsymbol{Ad}_{\boldsymbol{\mathcal{T}}_{\mathcal{A}\mathcal{B}}}\boldsymbol{J}_\mathcal{B}.
    \label{eq,Jacobianrelation}
\end{equation}
The spatial form of kinematics has provided a path to calculate the spatial Jacobian of frames from previously known spatial Jacobians. Also, it enables the use of the twist of any arbitrary body-fixed frame to be used in the calculation of the kinetic energy.

\subsection{Generalized Inertia Matrix}
Kinetic energy of the $i^\text{th}$ rigid body in a multibody system can be calculated as
\begin{equation}
    T_i=\frac{1}{2}\boldsymbol{\mathcal{V}}_{\mathcal{B}_i}^T\boldsymbol{M}_i\boldsymbol{\mathcal{V}}_{\mathcal{B}_i}
    \label{eq,localkineticenergy}
\end{equation}
for an arbitrary body-fixed frame $\{\mathcal{B}_i\}$. The spatial inertia matrix $\boldsymbol{M}_i$ is defined as \cite{Featherstone2008}
\begin{equation}
    \boldsymbol{M}_i=\begin{bmatrix}
        \boldsymbol{I}_{\mathcal{C}_i}+m_i\left[\boldsymbol{p}_{\mathcal{B}_i\mathcal{C}_i}\right]\left[\boldsymbol{p}_{\mathcal{B}_i\mathcal{C}_i}\right]^T & m_i\left[\boldsymbol{p}_{\mathcal{B}_i\mathcal{C}_i}\right] \\ -m_i\left[\boldsymbol{p}_{\mathcal{B}_i\mathcal{C}_i}\right] & m_i \boldsymbol{\mathcal{I}}_3
    \end{bmatrix}
    \label{eq,spatialinertia}
\end{equation}
with $m_i$ being the mass of the body and $\boldsymbol{I}_{\mathcal{C}_i}$ denoting the inertia tensor of the body with respect to frame $\{\mathcal{C}_i\}$ that coincides with the center of mass. The vector $\boldsymbol{p}_{\mathcal{B}_i\mathcal{C}_i}$ is the distance vector from the origin of frame $\{\mathcal{B}_i\}$ to the origin of frame $\{\mathcal{C}_i\}$ and expressed in $\{\mathcal{B}_i\}$. $\boldsymbol{\mathcal{I}}_3$ denotes the identity matrix of order 3.
By replacing the velocities in (\ref{eq,localkineticenergy}) with (\ref{eq,twistJacobian}), the kinetic energy of the body is rewritten as
\begin{equation}
    T_i=\frac{1}{2}\boldsymbol{\dot{q}}^T \boldsymbol{\Gamma}_i(\boldsymbol{q})\boldsymbol{\dot{q}}.
    \label{eq,newlocalkinetic}
\end{equation}
The matrix $\boldsymbol{\Gamma}_i$ is the generalized inertia matrix of the $i^\text{th}$ body and can be calculated as
\begin{equation}
    \boldsymbol{\Gamma}_i=\boldsymbol{J}_{\mathcal{B}_i}^T\boldsymbol{M}_i\boldsymbol{J}_{\mathcal{B}_i}.
    \label{eq,generalizedinertialmatrix}
\end{equation}
The local generalized inertia matrices show how the inertial properties of each body are related to the motion generated by each individual state in the minimal generalized coordinates set. In other words, it represents the interconnection between bodies in the multibody system. Since the total kinetic energy of the multibody system is the sum of the kinetic energies of individual bodies, it can easily be deduced that the global generalized inertia matrix $\boldsymbol{\Gamma}$ is the sum of local generalized inertia matrices; thus,
\begin{equation}
    T=\sum_i T_i=\frac{1}{2}\boldsymbol{\dot{q}}^T \left(\sum_i\boldsymbol{\Gamma}_i(\boldsymbol{q})\right)\boldsymbol{\dot{q}}=\frac{1}{2}\boldsymbol{\dot{q}}^T \boldsymbol{\Gamma}(\boldsymbol{q})\boldsymbol{\dot{q}}.
    \label{eq,globalkinetic}
\end{equation}
Inspecting (\ref{eq,generalizedinertialmatrix}) and considering the symmetric positive definiteness of the spatial inertia matrix of each body, one observes that local generalized inertia matrices and, thus, the global one are symmetric positive definite.

\subsection{Modular Equations of Motion}
We assume that the Lagrangian of any rigid body, and thus the whole multibody system, consists of only the kinetic energy ($L=T$), and the effects of gravity and energy-conserving elements such as springs are treated as generalized forces $\boldsymbol{Q}$. The generalized forces corresponding to the external wrenches $\boldsymbol{W}_k$ applied to the origin of frame $\{\mathcal{P}_k\}$ can be calculated by \cite{Featherstone2008}
\begin{equation}
    \boldsymbol{Q}=\sum_k \boldsymbol{J}_{\mathcal{P}_k}^T \boldsymbol{W}_k
    \label{eq,generalizedforce}
\end{equation}
Using the Lagrangian formalism \cite{shabana2009computational}
\begin{equation}
    \frac{d}{dt}\left(\frac{\partial L}{\partial \boldsymbol{\dot{q}}}\right)-\frac{\partial L}{\partial \boldsymbol{q}}=\boldsymbol{Q},
    \label{eq,Lagrange}
\end{equation}
and using (\ref{eq,globalkinetic}) in accordance with the mentioned assumptions, the equations of motion for the whole system are achieved as
\begin{equation}
    \boldsymbol{\Gamma}\boldsymbol{\ddot{q}}+\boldsymbol{\dot{\Gamma}}\boldsymbol{\dot{q}}-\frac{1}{2}\nabla_{\boldsymbol{q}}\left(\boldsymbol{q}^T\boldsymbol{\Gamma}\boldsymbol{q}\right)=\boldsymbol{Q}
    \label{eq,totalEoM}
\end{equation}
with the equation of motion corresponding to the $j^\text{th}$ generalized coordinate being
\begin{equation}
    \left(\boldsymbol{\Gamma}\boldsymbol{\ddot{q}}+\boldsymbol{\dot{\Gamma}}\boldsymbol{\dot{q}}\right)_j-\frac{1}{2}\boldsymbol{\dot{q}}^T \left(\frac{\partial\boldsymbol{\Gamma}}{\partial q_j}\right)\boldsymbol{\dot{q}}=Q_j.
    \label{eq,EoM}
\end{equation}
The partial derivative $\dfrac{\partial \boldsymbol{\Gamma}}{\partial q_j}$ is element-wise, and the first time derivative of the generalized inertia matrix can be obtained using the chain rule
\begin{equation}
    \boldsymbol{\dot{\Gamma}}=\sum_j \left(\frac{\partial \boldsymbol{\Gamma}}{\partial q_j} \dot{q}_j\right).
    \label{eq,chainrule}
\end{equation}
The modular equations of motion obtained in this part will be  used in the model-based controller design in Section \ref{sec:control}.

\subsection{Multibody System Reconfiguration}
Due to the modular structure of the global generalized inertia matrix, as seen in (\ref{eq,globalkinetic}), any modification to the multibody system can be directly addressed by modifying its corresponding local generalized inertia matrix. If the physical parameters of a rigid body change, only the corresponding terms in the local generalized inertia matrix of that body need to be modified. For the case where the number of degrees of freedom is increased, there are two possibilities:
\begin{enumerate}
    \item If the new degrees of freedom have an impact on the motion of the body, the spatial Jacobian needs to be re-derived to form the new local generalized inertia matrix.
    \item If the motion of the body is independent from the newly added degrees of freedom, the former spatial Jacobian only needs to be augmented in a larger matrix with zero columns, with the same number as the number of the newly added degrees of freedom, to form the new spatial Jacobian with dimensional consistency, i.e.,
    \begin{equation}
        \boldsymbol{J}_{new}=\begin{bmatrix} \boldsymbol{J}_{old} & \boldsymbol{0}\end{bmatrix}.
        \label{eq,augmentJacob}
    \end{equation}
\end{enumerate}
In the case where degrees of freedom are excluded from the motion states of a system, a selection matrix $\boldsymbol{P}$ can be easily formed so that
\begin{equation}
    \boldsymbol{q}_{new}=\boldsymbol{P}\boldsymbol{q}_{old}
    \label{eq,selectq}
\end{equation}
holds and as a result, the new spatial Jacobian will be 
\begin{equation}
    \boldsymbol{J}_{new}= \boldsymbol{J}_{old}\boldsymbol{P}^T.
    \label{eq,selectJacob}
\end{equation}
Since the matrix $\boldsymbol{P}$ is the same for all the bodies, it can be deduced from (\ref{eq,generalizedinertialmatrix}) and (\ref{eq,globalkinetic}) that the selection matrix can be directly applied to the global generalized inertia matrix as
\begin{equation}
    \boldsymbol{\Gamma}_{new}=\boldsymbol{P}\boldsymbol{\Gamma}_{old}\boldsymbol{P}^T,
    \label{eq,selectGamma}
\end{equation}
which is computationally more efficient compared to the selection matrix being applied to the Jacobians of each body individually.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Modular Control of Multibody Systems}
\label{sec:control}
In this section, we use the modular equation of motion to design a stable model-based controller where the control law for each degree of freedom will be available individually. Then, the local control laws are merged in a single system to ensure global stability.

\subsection{Controller Design}
We can separate the generalized force $Q_j$ into two independent parts $Q_{j,e}$ and $Q_{j,a}$ where the former stands for external generalized forces applied to the system, including gravity, and the latter shows the generalized force coming from the actuator. To propose a control law for the actuator, we assume that $Q_{j,a}$ consists of two signals: $Q_{j,m}$ coming from the inverse dynamics as the model-based part, and $Q_{j,c}$ as the correction term based on errors to account for uncertainties. Therefore, we can write the model-based part of the control law from (\ref{eq,EoM}) for the desired trajectory $\boldsymbol{q}_d$ as
\begin{equation}
    Q_{j,m}=\left(\boldsymbol{\Gamma}\boldsymbol{\ddot{q}}_d+\boldsymbol{\dot{\Gamma}}\boldsymbol{\dot{q}}_d\right)_j-\frac{1}{2}\boldsymbol{\dot{q}}_d^T \left(\frac{\partial\boldsymbol{\Gamma}}{\partial q_j}\right)\boldsymbol{\dot{q}}_d-Q_{j,e}.
    \label{eq,Qm}
\end{equation}
It should be noted that in (\ref{eq,Qm}) and the remainder of this paper, the generalized inertia matrix is evaluated using actual values of the coordinates $\boldsymbol{q}$ instead of the desired values. By defining the tracking error as
\begin{equation}
        \boldsymbol{e}=\boldsymbol{q}-\boldsymbol{q}_d,
        \label{eq,error}
\end{equation}
we propose the correction term of the control law to be
\begin{equation}
    Q_{j,c}=-k_{j,p}e_j-k_{j,v}\dot{e}_j
    \label{eq,Qc}
\end{equation}
where the terms $k_{j,p}$ and $k_{j,v}$ are positive control gains. With the mentioned assumptions, the control law $Q_{j,a}$ can be obtained as
\begin{equation}
    \scalebox{0.875}{$\begin{aligned}
    &Q_{j,a}=Q_{j,m}+Q_{j,c}=\\& \left(\boldsymbol{\Gamma}\boldsymbol{\ddot{q}}_d+\boldsymbol{\dot{\Gamma}}\boldsymbol{\dot{q}}_d\right)_j-\frac{1}{2}\boldsymbol{\dot{q}}_d^T \left(\frac{\partial\boldsymbol{\Gamma}}{\partial q_j}\right)\boldsymbol{\dot{q}}_d-Q_{j,e}-k_{j,p}e_j-k_{j,v}\dot{e}_j.
    \end{aligned}$}
    \label{eq,control}
\end{equation}
For the whole system, the control law becomes
\begin{equation}
    \scalebox{0.9}{$
    \boldsymbol{Q}_a=\boldsymbol{\Gamma}\boldsymbol{\ddot{q}}_d+\boldsymbol{\dot{\Gamma}}\boldsymbol{\dot{q}}_d-\frac{1}{2}\nabla_{\boldsymbol{q}}\left(\boldsymbol{\dot{q}}_d^T\boldsymbol{\Gamma}\boldsymbol{\dot{q}}_d\right)-\boldsymbol{Q}_e-\boldsymbol{K}_p\boldsymbol{e}-\boldsymbol{K}_v\boldsymbol{\dot{e}}
    $}
    \label{eq,totalcontrol}
\end{equation}
where $\boldsymbol{K}_p$ and $\boldsymbol{K}_v$ are diagonal matrices with elements $k_{j,p}$ and $k_{j,v}$, respectively. With the assumption mentioned for control gains, these matrices are symmetric positive definite. The global form of the controller enables stability analysis for the whole system.

\subsection{Stability Analysis}
To establish the stability of the proposed controller, consider the Lyapunov function candidate 
\begin{equation}
    \scalebox{0.925}{$
    V=\frac{1}{2}\left(\boldsymbol{\dot{q}}-\boldsymbol{\dot{q}}_d\right)^T\boldsymbol{\Gamma}\left(\boldsymbol{\dot{q}}-\boldsymbol{\dot{q}}_d\right)+\frac{1}{2}\left(\boldsymbol{q}-\boldsymbol{q}_d\right)^T\boldsymbol{K}_p\left(\boldsymbol{q}-\boldsymbol{q}_d\right)$}
    \label{eq,Lyapunov}
\end{equation}
to capture both position and velocity errors. Taking the first derivative of (\ref{eq,Lyapunov}) yields
\begin{equation}
    \scalebox{0.975}{$
    \begin{aligned}
    \dot{V}=&\left(\boldsymbol{\dot{q}}-\boldsymbol{\dot{q}}_d\right)^T\boldsymbol{\Gamma}\left(\boldsymbol{\ddot{q}}-\boldsymbol{\ddot{q}}_d\right)+\frac{1}{2}\left(\boldsymbol{\dot{q}}-\boldsymbol{\dot{q}}_d\right)^T\boldsymbol{\dot{\Gamma}}\left(\boldsymbol{\dot{q}}-\boldsymbol{\dot{q}}_d\right)\\+&\left(\boldsymbol{\dot{q}}-\boldsymbol{\dot{q}}_d\right)^T\boldsymbol{K}_p\left(\boldsymbol{q}-\boldsymbol{q}_d\right).
    \end{aligned}$}
    \label{eq,Vdot}
\end{equation}
Substituting the system dynamics (\ref{eq,totalEoM}) and the control law (\ref{eq,totalcontrol}) into (\ref{eq,Vdot}) results in
\begin{equation}
    \scalebox{0.85}{$
    \begin{aligned}
    \dot{V}=&-\left(\boldsymbol{\dot{q}}-\boldsymbol{\dot{q}}_d\right)^T\boldsymbol{K}_v\left(\boldsymbol{\dot{q}}-\boldsymbol{\dot{q}}_d\right)+\left(\boldsymbol{\dot{q}}-\boldsymbol{\dot{q}}_d\right)^T\left(-\frac{1}{2}\boldsymbol{\dot{\Gamma}}\left(\boldsymbol{\dot{q}}-\boldsymbol{\dot{q}}_d\right)\right. \\&
    \left. +\frac{1}{2}\left(\nabla_{\boldsymbol{q}}\left(\boldsymbol{\dot{q}}^T\boldsymbol{\Gamma}\boldsymbol{\dot{q}}\right)-\nabla_{\boldsymbol{q}}\left(\boldsymbol{\dot{q}}_d^T\boldsymbol{\Gamma}\boldsymbol{\dot{q}}_d\right)\right)\right).
    \end{aligned}$}
    \label{eq,Vdot2}
\end{equation}
To manage the second term in (\ref{eq,Vdot2}), we refer to the equivalent form of (\ref{eq,totalEoM}) that can be written as \cite{murray1994mathematical}
\begin{equation}
    \boldsymbol{\Gamma}\boldsymbol{\ddot{q}}+\boldsymbol{C}(\boldsymbol{q},\boldsymbol{\dot{q}})\boldsymbol{\dot{q}}=\boldsymbol{Q}
    \label{eq,AltEoM}
\end{equation}
that indicates the relation
\begin{equation}
    \boldsymbol{C}(\boldsymbol{q},\boldsymbol{\dot{q}})\boldsymbol{\dot{q}}=\boldsymbol{\dot{\Gamma}}\boldsymbol{\dot{q}}-\frac{1}{2}\nabla_{\boldsymbol{q}}\left(\boldsymbol{\dot{q}}^T\boldsymbol{\Gamma}\boldsymbol{\dot{q}}\right)
    \label{eq,C}
\end{equation}
holds. Therefore, (\ref{eq,Vdot2}) with addition and subtraction of the term $\boldsymbol{\dot{\Gamma}}\left(\boldsymbol{\dot{q}}-\boldsymbol{\dot{q}}_d\right)$ can be re-written as
\begin{equation}
    \begin{split}
    \dot{V}=&-\left(\boldsymbol{\dot{q}}-\boldsymbol{\dot{q}}_d\right)^T\boldsymbol{K}_v\left(\boldsymbol{\dot{q}}-\boldsymbol{\dot{q}}_d\right)\\&-\frac{1}{2}\left(\boldsymbol{\dot{q}}-\boldsymbol{\dot{q}}_d\right)^T\left(\boldsymbol{\dot{\Gamma}}-2\boldsymbol{C}\right)\left(\boldsymbol{\dot{q}}-\boldsymbol{\dot{q}}_d\right).
    \end{split}
    \label{eq,AltVdot}
\end{equation}
As shown in \cite{murray1994mathematical}, the matrix $\boldsymbol{\dot{\Gamma}}-2\boldsymbol{C}$ is skew-symmetric, and as a result,
\begin{equation}
    \dot{V}=-\left(\boldsymbol{\dot{q}}-\boldsymbol{\dot{q}}_d\right)^T\boldsymbol{K}_v\left(\boldsymbol{\dot{q}}-\boldsymbol{\dot{q}}_d\right).
    \label{eq,VdotFinal}
\end{equation}
Since $V$ is positive definite and radially unbounded, and $\dot{V}$ is negative definite, almost-global asymptotic stability of the closed-loop system is proved by standard Lyapunov theory \cite{slotine_li_1991} and considerations about the stability of mechanical systems with rotational degrees of freedom \cite{bhat2000topological}. The chosen control law is used in a simulation to determine the capability of the proposed modular approach.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Simulation and Results}
\label{sec:simulation}
To assess the performance of the inertia partitioning framework for multibody systems with closed chains, a simulation is conducted on a 3-DoF series-parallel manipulator using Simscape, the physics engine of MatLab.

\subsection{Kinematic and Dynamic Modeling}
The series-parallel manipulator in the simulation is shown in Fig. \ref{iso}. The frames of reference needed for modeling this multibody system, as well as the required geometric parameters, are illustrated in Fig. \ref{frames}. The orientation of frames is chosen such that they comply with the Simscape convention for mechanical joints. Frame $\{\mathcal{S}\}$ is the inertial frame of reference, while the other frames are body-fixed. Indices $1, 2,$ and $3$ refer to body-fixed frames attached to the base, the first link, and the second link, respectively.

\begin{figure}[t]
    \centering
    \includegraphics[width=\linewidth,trim={5cm, 2cm, 5cm, 1.5cm},clip]{Img/Iso}
    \caption{The schematics for a 3-DoF series-parallel manipulator. The degrees of freedom are made up of two actuated revolute joints for the base and the first link, and an actuated prismatic joint between the first and second links.}
    \label{iso}
\end{figure}

The generalized coordinates used to model this series-parallel manipulator are $\boldsymbol{q}=\begin{bmatrix}\phi & \theta & \delta\end{bmatrix}^T$, as depicted in Fig. \ref{frames}. The intermediate coordinate $\zeta$, which is a function of $\delta$, is also used for the sake of brevity and clarity in the mathematical relations. 
Using the twist relation (\ref{eq,twistJacobian}) and the Jacobian transformation (\ref{eq,Jacobianrelation}), while considering the mechanical joints in the 3-DoF series-parallel manipulator, the Jacobians of the required frames are calculated as
\begin{equation}
    \begin{split}
        \boldsymbol{J}_{\mathcal{B}_1}&=\begin{bmatrix}
            0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0
        \end{bmatrix}^T \\
        \boldsymbol{J}_{\mathcal{C}_1}&=\boldsymbol{Ad}_{\boldsymbol{\mathcal{T}}_{\mathcal{C}_1\mathcal{B}_1}}\boldsymbol{J}_{\mathcal{B}_1}\\
        \boldsymbol{J}_{\mathcal{F}_1}&=\boldsymbol{Ad}_{\boldsymbol{\mathcal{T}}_{\mathcal{F}_1\mathcal{B}_1}}\boldsymbol{J}_{\mathcal{B}_1}\\
        \boldsymbol{J}_{\mathcal{B}_2}&=\boldsymbol{Ad}_{\boldsymbol{\mathcal{T}}_{\mathcal{B}_2\mathcal{F}_1}}\boldsymbol{J}_{\mathcal{F}_1}+\begin{bmatrix}
             0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0
        \end{bmatrix}^T \\
        \boldsymbol{J}_{\mathcal{C}_2}&=\boldsymbol{Ad}_{\boldsymbol{\mathcal{T}}_{\mathcal{C}_2\mathcal{B}_2}}\boldsymbol{J}_{\mathcal{B}_2}\\
        \boldsymbol{J}_{\mathcal{F}_2}&=\boldsymbol{Ad}_{\boldsymbol{\mathcal{T}}_{\mathcal{F}_2\mathcal{B}_2}}\boldsymbol{J}_{\mathcal{B}_2}\\
        \boldsymbol{J}_{\mathcal{B}_3}&=\boldsymbol{Ad}_{\boldsymbol{\mathcal{T}}_{\mathcal{B}_3\mathcal{F}_2}}\boldsymbol{J}_{\mathcal{F}_2}+\begin{bmatrix}
             0 & 0 & 0 & 0 & 0 & 0  \\0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & \frac{\partial \zeta}{\partial \delta} & 0 & 0 & 0
        \end{bmatrix}^T \\
        \boldsymbol{J}_{\mathcal{C}_3}&=\boldsymbol{Ad}_{\boldsymbol{\mathcal{T}}_{\mathcal{C}_3\mathcal{B}_3}}\boldsymbol{J}_{\mathcal{B}_3}\\
        \boldsymbol{J}_{\mathcal{F}_3}&=\boldsymbol{Ad}_{\boldsymbol{\mathcal{T}}_{\mathcal{F}_3\mathcal{B}_3}}\boldsymbol{J}_{\mathcal{B}_3}.
    \end{split}
    \label{eq,simJacobian}
\end{equation}

The rotation matrices between the frames are obtainable from Figs. \ref{frames} and \ref{triangle}, and the non-zero required distance vectors are available in Table \ref{tab,dist}. For the relations to be fully expressed in the minimal set of generalized coordinates, a geometric relation is needed between the intermediate coordinate $\zeta$ and the degree-of-freedom coordinate $\delta$. The triangular closed chain between the first and second links is shown more closely in Fig. \ref{triangle}. Applying the law of cosines for angle $\gamma$ yields
\begin{equation}
    \cos(\gamma)=1-\dfrac{(l_0+\delta)^2}{2L_0^2}
    \label{eq,gamma}
\end{equation}
with $l_0$ being the length of the cylinder-piston pair when the piston is at its lower end, i.e., $\delta=0$. Because $\gamma$ and $\zeta$ are supplementary angles, and considering the fact that the value of rotation $\zeta$ with respect to the $z$-axis of frame $\{\mathcal{B}_3\}$ lies in a subset of the interval $(-\pi, 0)$, the relationship between $\zeta$ and $\delta$ is obtained as
\begin{equation}
    \zeta=-\cos^{-1}\left(\dfrac{(\delta+l_0)^2}{2L_0^2}-1\right),
    \label{eq,zeta}
\end{equation}
and as a result
\begin{equation}
    \frac{\partial\zeta}{\partial\delta}=\dfrac{\delta+l_0}{L_0^2\sqrt{1-\left(\dfrac{(\delta+l_0)^2}{2L_0^2}-1\right)^2}}\;.
    \label{eq,zetadelta}
\end{equation}

\begin{table}[b!]
    \centering
    \caption{Required non-zero distance vectors between the frames in the 3-DoF series-parallel manipulator. All values are in meters.}
    \begin{tabular}{|r@{\,=\,}l|}
        \hline
        \rule{0pt}{10pt}$\boldsymbol{p}_{\mathcal{B}_1\mathcal{C}_1}$ & $\begin{bmatrix} 0 & 0 & 0.103 \end{bmatrix}^T$ \\[2pt] \hline
        \rule{0pt}{10pt}$\boldsymbol{p}_{\mathcal{B}_1\mathcal{F}_1}$ & $\begin{bmatrix} 0 & 0.206 & 0.075 \end{bmatrix}^T$ \\[2pt] \hline
        \rule{0pt}{10pt}$\boldsymbol{p}_{\mathcal{B}_2\mathcal{C}_2}$ & $\begin{bmatrix} 0.959 & 0.001 & -0.077 \end{bmatrix}^T$ \\[2pt] \hline
        \rule{0pt}{10pt}$\boldsymbol{p}_{\mathcal{B}_2\mathcal{F}_2}$ & $\begin{bmatrix} 2 & 0 & 0 \end{bmatrix}^T$ \\[2pt] \hline
        \rule{0pt}{10pt}$\boldsymbol{p}_{\mathcal{B}_3\mathcal{C}_3}$ & $\begin{bmatrix} 1.041 & 0.001 & -0.077 \end{bmatrix}^T$ \\[2pt] \hline
        \rule{0pt}{10pt}$\boldsymbol{p}_{\mathcal{B}_3\mathcal{F}_3}$ & $\begin{bmatrix} 2 & 0 & 0 \end{bmatrix}^T$ \\[2pt] \hline
    \end{tabular}
    \label{tab,dist}
\end{table}

\begin{figure}[b!]
    \centering
    \begin{minipage}[h]{0.45\textwidth}
        \includegraphics[width=\linewidth,trim={5cm, 1cm, 2cm, 0cm},clip]{Img/Side}
        \label{side}
    \end{minipage}
    \hfill
    \begin{minipage}[h]{0.45\textwidth}
        \includegraphics[width=\linewidth,trim={7cm, 4cm, 7cm, 3cm},clip]{Img/Top}
        \label{top}
    \end{minipage}
    \hfill
    \caption{The required frames of reference and geometric parameters for modeling the 3-DoF series-parallel manipulator. Some frames are shown only in one of the figures for clarity. \textit{Top} - Side view perpendicular to the links' body-fixed frames $x$-$y$ plane. \textit{Bottom} - Top view of the base, perpendicular to its frames' $x$-$y$ planes.}
    \label{frames}
\end{figure}



With the Jacobians obtained from (\ref{eq,simJacobian}) and (\ref{eq,zetadelta}), the generalized inertia matrix for each body is obtained using (\ref{eq,generalizedinertialmatrix}) with their summation yielding the global generalized inertia matrix $\boldsymbol{\Gamma}$. To avoid extra complexity, it is assumed that the cylinder-piston pair has negligible mass compared to the rest of the system. To calculate the generalized forces corresponding to the weight of the bodies, the Jacobians $\boldsymbol{J}_{\mathcal{C}_i}$ and the rotation matrices $\boldsymbol{R}_{\mathcal{C}_i\mathcal{S}}$ of body-fixed frames $\{\mathcal{C}_i\}$ coincident with the centers of mass are needed. With those matrices available and assuming the weight to be an external wrench, the corresponding generalized forces can be calculated by using (\ref{eq,generalizedforce})
\begin{equation}
    \boldsymbol{Q}=\sum_i \boldsymbol{J}_{\mathcal{C}_i}^T \boldsymbol{W}_i=\sum_i\boldsymbol{J}_{\mathcal{C}_i}^T\begin{bmatrix}
        \boldsymbol{R}_{\mathcal{C}_i\mathcal{S}} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{R}_{\mathcal{C}_i\mathcal{S}}
    \end{bmatrix}\begin{bmatrix}
        \boldsymbol{0}\\-m_ig
    \end{bmatrix}
    \label{eq,weight}
\end{equation}
with $g$ denoting the gravitational acceleration of the Earth, and zero matrices $\boldsymbol{0}$ having the appropriate dimensions. With all the required components of the control law (\ref{eq,control}) available either by calculation or by selection, the chosen 3-DoF series-parallel manipulator can be controlled for trajectory tracking.

\begin{figure}[t]
    \centering
    \includegraphics[width=0.75\linewidth,trim={5cm, 1cm, 5cm, 1.5cm},clip]{Img/Tri}
    \caption{The geometry of the triangular closed chain. This closed chain consists of one actuated prismatic joint and three passive revolute joints.}
    \label{triangle}
\end{figure}

\subsection{Trajectory Tracking}
The values corresponding to the physical parameters of the 3-DoF series-parallel manipulator are expressed in Table \ref{tab,param}. Using the geometry available in Fig. \ref{frames}, the relationship between the Cartesian coordinates of the end effector and the generalized coordinates can be expressed as
\begin{equation}
    \begin{bmatrix}
        x \\ y \\ z
    \end{bmatrix}=\begin{bmatrix}
        L \cos(\phi)(\cos(\theta)+\cos(\theta-\zeta)) \\ L \sin(\phi)(\cos(\theta)+\cos(\theta-\zeta)) \\ L(\sin(\theta)+\sin(\theta-\zeta))
    \end{bmatrix}.
    \label{eq,IK}
\end{equation}

\begin{table}[b]
\caption{Physical parameter values for the 3-DoF series-parallel manipulator}
\centering
\begin{tabular}{|l|c|}
\hline
\rule{0pt}{10pt}Parameter & Value  \\[2pt]
\hline
\rule{0pt}{10pt}$L$ &  $2$   $m$ \\[2pt] \hline
\rule{0pt}{10pt}$L_0$ &  $0.35$   $m$ \\[2pt] \hline
\rule{0pt}{10pt}$l_0$&  $0.425$  $m$  \\[2pt] \hline
\rule{0pt}{10pt}$m_1$&  $20$  $kg$  \\[2pt] \hline
\rule{0pt}{10pt}$m_2$&  $60$  $kg$  \\[2pt] \hline
\rule{0pt}{10pt}$m_3$&  $60$  $kg$  \\[2pt] \hline
\rule{0pt}{15pt}$\boldsymbol{I}_{\mathcal{C}_1}$&  $\begin{bmatrix}
    0.536&0&0& \\ 0&0.554&0& \\ 0&0&0.789&
\end{bmatrix}$  $kg.m^2$  \\[10pt] \hline
\rule{0pt}{15pt}$\boldsymbol{I}_{\mathcal{C}_2}$&  $\begin{bmatrix}
    0.311&-0.065&0.098 \\ -0.065&22.7&-0.003 \\ 0.098&-0.003&22.8
\end{bmatrix}$  $kg.m^2$  \\[10pt] \hline
\rule{0pt}{15pt}$\boldsymbol{I}_{\mathcal{C}_3}$&  $\begin{bmatrix}
    0.311&0.065&-0.098 \\ 0.065&22.7&-0.003 \\ -0.098&-0.003&22.8
\end{bmatrix}$  $kg.m^2$  \\[10pt] \hline
\rule{0pt}{15pt}$g$&  $9.81$  $\dfrac{m}{s^2}$  \\[2pt] \hline
\rule{0pt}{10pt}$k_{p,i}$ and $k_{v,i}$&  $20$   \\[2pt] \hline
\end{tabular}
\label{tab,param}
\end{table}

The desired trajectory for the end effector in the simulation is illustrated in Fig. \ref{cart}, both in three-dimensional Cartesian space and with respect to time. Using (\ref{eq,IK}), the desired trajectory is expressed in the generalized coordinates with respect to time in Fig. \ref{conf}. This trajectory is considered in the control law (\ref{eq,control}) for the trajectory tracking task.

\begin{figure}[b!]
    \centering
    \begin{minipage}[h]{0.45\textwidth}
        \includegraphics[width=\linewidth]{Img/traj_xyz_a.pdf}
    \end{minipage}
    \hfill
    \begin{minipage}[h]{0.45\textwidth}
        \includegraphics[width=\linewidth]{Img/traj_xyz_b.pdf}
    \end{minipage}
    \hfill
    \caption{Desired trajectory for the end effector of the 3-DoF series-parallel manipulator. \textit{Top} - The desired trajectory in Cartesian space. \textit{Bottom} - The desired trajectory with respect to time.}
    \label{cart}
\end{figure}

\begin{figure}[t]
    \centering
    \includegraphics[width=\linewidth]{Img/traj_q}
    \caption{Desired trajectory for the end effector of the 3-DoF series-parallel manipulator expressed in generalized coordinates.}
    \label{conf}
\end{figure}

With the desired trajectory and the necessary physical parameters, and the previously calculated variables, the control law (\ref{eq,control}) is applied to the manipulator. The results of the trajectory tracking are presented in Figs. \ref{track} and \ref{error}. Fig. \ref{track} shows the followed trajectory against the desired trajectory in Cartesian space, and also with respect to time. Fig. \ref{error} illustrates the tracking error with respect to time. The root mean squared error (RMSE) values for the tracking error in Cartesian space are reported in Table \ref{tab,rmse}.

\begin{figure}[t]
    \centering
    \begin{minipage}[b]{0.45\textwidth}
        \includegraphics[width=\linewidth]{Img/track_xyz_a.pdf}
    \end{minipage}
    \hfill
    \begin{minipage}[b]{0.45\textwidth}
        \includegraphics[width=\linewidth]{Img/track_xyz_b.pdf}
    \end{minipage}
    \hfill
    \caption{Trajectory tracking results for the end effector of the 3-DoF series-parallel manipulator. \textit{Top} - Desired and actual trajectories in Cartesian space. \textit{Bottom} - Comparison of desired and actual trajectories with respect to time.}
    \label{track}
\end{figure}

\begin{figure}[t!]
    \centering
    \includegraphics[width=\linewidth]{Img/error}
    \caption{Trajectory tracking error in Cartesian space with respect to time.}
    \label{error}
\end{figure}

\begin{table}[t]
\caption{Root mean squared error values for trajectory tracking in Cartesian space}
\centering
\begin{tabular}{|c|c|c|c|}
\hline
\rule{0pt}{10pt}  & X & Y & Z  \\[2pt]
\hline
\rule{0pt}{10pt}RMSE & 0.002 m & 0.003 m & 0.014 m \\[2pt] \hline
\end{tabular}
\label{tab,rmse}
\end{table}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusion}
\label{sec:conclusion}
The challenges associated with handling closed kinematic chains in modular control of multibody systems with existing methods motivated us to propose a modular control framework that inherently handles closed kinematic chains. The inertia partitioning framework presented in this paper benefits from a different look at the concept of modularity and provides a straightforward procedure to design controllers for trajectory tracking of multibody systems. Choosing the degree-of-freedom states as the modules instead of individual rigid bodies bypasses the need to calculate the constraint forces, which is the main challenge in existing modular control approaches to closed kinematic chains. 

In order to demonstrate the controller design procedure and the performance of the inertia partitioning framework, a 3-DoF series-parallel manipulator was chosen as a multibody system containing a closed kinematic chain. The controller was designed in a scalable step-by-step process, and it was applied on the series-parallel manipulator.Trajectory tracking results indicate desirable performance in tracking accuracy. The simulation is a testament to the achievement of the goals of this research.

The proposed approach presents the first step toward a generic modular control framework. Further studies will strengthen the approach by addressing reconfiguration, uncertainty compensation, and experimental implementation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Acknowledgement}
We acknowledge the financial support of the Finnish Ministry of Education and Culture through the Intelligent Work Machines Doctoral Education Pilot Program (IWM VN/3137/2024-OKM-4).

\bibliography{ref.bib}
\bibliographystyle{ieeetr}
\end{document}