A modular computational framework for the dynamic analyses of cable-driven parallel robots with different types of actuation including the effects of inertia, elasticity and damping of cables

Abstract

Dynamic simulations of the cable-driven parallel robots (CDPRs) with cable models closer to reality can predict the motions of moving platforms more accurately than those with idealisations. Hence, the present work proposes an efficient and modular computational framework for this purpose. The primary focus is on the developments required in the context of CDPRs actuated by moving the exit points of cables while the lengths are held constant. Subsequently, the framework is extended to those cases where simultaneous changes in the lengths of cables are employed. Also, the effects due to the inertia, stiffness and damping properties of the cables undergoing 3D motions are included in their dynamic models. The efficient recursive forward dynamics algorithms from the prior works are utilised to minimise the computational effort. Finally, the efficacy of the proposed framework and the need for such an inclusive dynamic model are illustrated by applying it to different application scenarios using the spatial $4$-$4$ CDPR as an example.

Nomenclature

CDPR

cable-driven parallel robot

MP

moving platform

DoF

degree(s)-of-freedom

VACTS

variable aerial cable towed system (see [1])

FD

forward dynamics

ID

inverse dynamics

EoM

equation of motion

SDE

spring and damper element

RSSLM

recursive sub-system-level Lagrange multiplier approach (see [2])

MRFE

modified rigid finite element (see, e.g., [3])

GIM

generalised inertia matrix (see, e.g., [4])

DAE

differential-algebraic equation

$n_{\textrm{k}}$ - $n_{\textrm{r}}$ CDPR

a CDPR with distinct $n_{\textrm{k}}$ exit and $n_{\textrm{r}}$ anchoring points of cables (see, e.g., [5], pp. $29$ - $30$ )

$\boldsymbol{{q}},{\dot{\boldsymbol{{q}}}},{\ddot{\boldsymbol{{q}}}} \in{\mathbb{R}^{n}}$

joint positions, velocities and accelerations, respectively

${\boldsymbol{{M}}} \in GL_{n}({\mathbb{R}})$

GIM of the robot

${\boldsymbol{{C}}} \in{\mathbb{R}^{n \times n}}$

centripetal and Coriolis matrix (see Eq. (2))

$\boldsymbol{{G}} \in{\mathbb{R}^{n}}$

generalised forces due to the gravitational potential

${\boldsymbol{\tau }^{\textrm{s}}} \in{\mathbb{R}^{n}}$

generalised forces due to the SDEs

${\boldsymbol{\tau }^{\textrm{c}}} \in{\mathbb{R}^{n}}$

generalised constraint forces (see Eq. (6))

$\boldsymbol{\eta } \in{\mathbb{R}^{m}}$

constraint functions

$\boldsymbol{\lambda } \in{\mathbb{R}^{m}}$

Lagrange multipliers appearing in the EoM (see Eq. (1))

1. Introduction

A cable-driven parallel robot (CDPR) manipulates the motion of a moving platform (MP) with the help of multiple cables connected to it. The CDPRs provide large workspaces, low inertias and high load-carrying capacities compared to their self-weights. Nevertheless, the analyses of these robots’ mechanics are challenging due to the cables’ unilaterality, that is, the ability to pull but not push, anisotropic material properties causing deflections to occur in different time scales, sagging due to their own weights and non-rigidity leading to large degrees-of-freedom (DoF). The actuation schemes employed in CDPRs also differ significantly from conventional manipulators with rigid links. Following the classification proposed in ref. [6], the two kinds of actuation schemes commonly used in CDPRs are designated in this article as Type I and Type II, respectively. In the case of Type I CDPRs, the MP is manipulated by changing the lengths of the cables hoisting it. The cables themselves are fed and retrieved via stationary guidance systems (see, e.g., Fig. 2.1 in ref. [5], p. $17$ ). In contrast, in Type II CDPRs, the cables’ lengths are held constant, whereas their terminal points (i.e., the ones not connected to the MP) are moved. These terminal points are known in CDPR parlance as the exit points of the cables (see Fig. 2).

The primary focus of the present work is to develop a modular computational framework to simulate the dynamics of Type II CDPRs with cable models closer to the physical system, given the kinematic or dynamic actuation inputs. Also, extend the analyses to those scenarios where Type I actuation is applied simultaneously by employing the existing recursive forward dynamics algorithms for computational efficiency. The motivation for the above-mentioned directions of research is presented below vis-à-vis the previous works in this field.

The modularity of the proposed framework helps accommodate diverse architectures and actuators employed for different applications of Type II CDPRs in the past. For instance, the Wirepuller-Arm-driven Redundant Parallel manipulator, an $8$ - $3$ CDPR proposed in ref. [7] for quick assembly of light payloads, employs arms driven by motors to manipulate the exit points of the cables. In ref. [8], mobile robots were employed for cooperative towing of the MPs of planar $2$ - $2$ and $3$ - $3$ CDPRs with applications in warehouse operations. Likewise, a spatial $3$ - $3$ CDPR actuated by aerial vehicles with potential applications in rescue operations and heavy payload transfers was demonstrated in ref. [9].

Further, the inclusion of simultaneous utilisation of both kinds of actuation schemes in the scope of the developed framework accommodates the following application scenarios. The Type II actuation was used for the active reconfiguration ¹ of a Type I CDPR in the Large Vessel Interface Lift-on and Lift-off (LVI Lo/Lo) crane system. It was designed to minimise the relative motions between the containers and ships while transferring cargo between two ships in the open seas with the help of a Type I actuated $8$ - $8$ CDPR attached to the end-effector of a serial crane (see, e.g., ref. [10]). A similar idea for micro-macro manipulation in ship replenishment was presented in ref. [11], where the winches of a Type I suspended $6$ - $6$ CDPR were attached to a helicopter. Further, in ref. [12], a Type I suspended $3$ - $3$ CDPR was proposed to displace its winches via trolleys over rails and attain control of the 6 DoF of the MP for cargo handling. A Variable Aerial Cable Towed System (VACTS) was introduced in ref. [1] by fastening the winches of a Type I suspended $6$ - $6$ CDPR on four aerial vehicles to transfer loads in cluttered environments. Recently, the utility of FASTKIT, a Type I $8$ - $8$ CDPR mounted on top of two mobile robots, was demonstrated in ref. [13] as a low-cost solution for logistics in small warehouses.

The existing methodologies to solve the forward kinetostatic ² problems of Type I CDPRs can be readily applied to those of the Type II CDPRs and vice versa due to the similarities in the inputs to these problems. In contrast, a clear distinction exists between the input actuator variables and the architecture parameters of these manipulators while analysing their dynamics due to the variations of the former with time. Hence, the formulations for forward dynamics (FD) analysis of Type II CDPRs are different from those of Type I robots (see, e.g., [14]), requiring individual developments to incorporate the distinctions in cables’ actuation. Nevertheless, a few common features, such as the dynamic models of cables, moving platforms, and their interactions can be utilised from the formulations developed for Type I CDPRs. The current work employs one such modular and computationally efficient formulation recently reported in ref. [15] and focuses on the required additional developments.

The formulations of dynamics of Type II CDPRs wherein manned or unmanned aerial vehicles imparting motions of the cables’ exit points were extensively studied in the literature. The primary focus had been on the dynamic models of vehicles, their formations and interactions with the MP for applications in control. Consequently, the cable models were simplified significantly, that is, the cables were usually treated as massless and inelastic force transmission elements. The constrained equations of motion (EoM) were derived using either the Newton-Euler formulation, as in ref. [16] or the Lagrangian formulation, as in ref. [17], in a coordinate-free form. Such representations are compact and avoid the effects of singularities arising from the choice of parameters used to denote the rigid bodies’ orientations. Although these features are preferred to establish theoretical foundations, they are not ideal for performing FD simulations. Hence, it is typical to resort to commercially available software, for example, SimMechanics (Simscape Multibody) in ref. [14], for FD simulations, even though inverse dynamics (ID) and control are based on custom-built dynamic models. Further, the benchmark dynamic models used for this purpose should be more accurate than those used to develop the controller. Therefore, monitoring such custom developments in the past is imperative to ensure that the framework developed in this work meets the necessary requirements.

For instance, in ref. [18], the elastic and damping properties of the cables of a planar $2$ - $1$ CDPR were introduced as massless spring and damper elements (SDEs) to demonstrate the effect of cables’ stretching on aerial vehicles connected to them. A custom framework called AuRoRA platform with virtual vehicle models was employed to perform the simulations. The cable tensions were first explicitly obtained from a separate dynamic model and were then provided as disturbances acting on the vehicles. Such a sequential treatment disregards the coupled dynamics, limiting their framework’s utility and fidelity of the results it produced. Contrary to the above-mentioned approach, in ref. [19], a commercial software, MSC ADAMS, was employed to model the cables as flexible aluminium bars, that is, massive and elastic, and then imported it to MATLAB/Simulink for realising the non-linear coupled dynamics of a $2$ - $2$ CDPR driven by quadrotors. Nevertheless, in specific situations where the MP is subjected to external disturbances or the vehicles are undergoing agile manoeuvres, the assumption that cables remain taut at all times is invalid (see, e.g., Section 4). Therefore, efforts towards incorporating the lateral and transverse flexibilities of cables into the dynamic model can also be seen in the following works.

An early development in this regard is presented by Goodarzi et al. in ref. [20], wherein the cables were modelled as serial chains of rigid links connected by universal joints, and quadcopters were used for moving their exit points. However, the EoM of the CDPR was developed using a coordinate-free form of the Lagrangian formulation with the intention of its utilisation in developing the controller. Recently, in ref. [21], the cables were approximated as multiple lumped masses connected by SDEs to simulate the dynamics of a spatial $2$ - $1$ CDPR driven by quadcopters. Additionally, the effects of air resistance and contact forces from the ground were included in the model. The developed dynamic model was used to test the efficacy of the proposed controller.

As is evident from the above survey of the existing literature, the dynamic models of cables considered in previous works have evolved from a massless inelastic one to a serial chain of point masses connected by SDEs. Therefore, the current work employs one such detailed cable model in developing the computational framework to meet the upcoming demands in this field. Furthermore, it addresses the shortcomings mentioned in the existing formulations of dynamics so that the proposed framework can perform the FD simulations in an efficient manner.

A few prominent works on the CDPRs that incorporate both Type II and Type I actuation schemes are mentioned below to re-emphasise that the research gaps addressed in this work are equally relevant in such scenarios. In ref. [11], the Newton-Euler formulation was used to derive the EOM of a $6$ - $6$ CDPR with its winches mounted on a helicopter, assuming cables to be massless and inelastic force transmission elements. The FD simulations of the CDPR and helicopter, along with the proposed controllers, were performed using a framework developed in MATLAB/Simulink. The same approach was followed in the case of an $8$ - $8$ CDPR suspended through a heavy airship in ref. [22] emphasising the added complexity due to the extra cable winches. In ref. [1], such systems were referred to as VACTS and an experimental demonstration of their feasibility via a spatial $3$ - $1$ CDPR was reported for the first time. These developments justify the extension of the proposed framework’s scope for analysing such systems. Moreover, the dynamic models of cables and the associated formulation of dynamics considered in the past are still in their preliminary stages compared to the proposed work. This limitation is also pertinent to systems other than VACTS, that is, those actuated by means other than aerial vehicles.

For instance, in refs. [23] and [12], the cables were assumed to be massless and inelastic in deriving the EoM of a spatial $3$ - $3$ CDPR whose winches were moved using trolleys on rails. The incremental change over a decade is in considering the trolley’s dynamics in [12], while it was ignored in [23]. In the same category, an $8$ - $8$ CDPR with four of its eight pulleys moved over vertical rails with applications in automated masonry was recently presented in ref. [24]. The cables were modelled as massless springs, and the EOM were obtained using the Newton-Euler formulation. Despite the existing limitations, a trend towards detailed cable models can be observed. A similar tendency exists in the cases where mobile vehicles on the ground were used to move the exit points of CDPRs.

In ref. [25], the cables were treated as massless and inelastic force transmission elements, and the EoM of a $3$ - $3$ CDPR with each cable connected to a mobile crane was established. Recently, in ref. [26], the dynamic model of a spatial $4$ - $1$ CDPR attached to four Turtlebots, termed MoPICK, was developed in the V-REP environment. The cables were approximated to be serial chains of multiple massive cylindrical links connected by prismatic joints. In ref. [27], the Rayleigh–Ritz method was employed to model the non-linear longitudinal vibrations of the cables of a mobile ICaSbot, that is, a $6$ - $6$ CDPR mounted on a mobile base platform.

The significant contributions of this work can be summarised as follows. As is apparent from the reported survey of the state-of-the-art, many combinations of actuators and CDPRs were investigated in the past. However, apart from a few recent attempts, little to no effort has been spent in bringing the forward dynamic models of Type II CDPRs closer to their physical realities. The present work bridges this gap by proposing a computational framework to comprehensively model and analyse the dynamics of Type II CDPRs with kinematic and dynamic actuation inputs.

The said framework is detailed and realistic; for example, it takes into consideration the inertia, elasticity and damping of the cables while modelling their dynamics. Also, it is modular by design. In this context, “modular” implies that the components of the CDPR, such as the cables, the MP and the driving mechanisms, are treated as independent sub-systems. It does not refer to different models of links and joints used to represent these sub-systems, as in ref. [28]. For instance, as shown in Figs. 3 and 5, the spatial $4$ - $4$ CDPR was notionally decomposed to four cables, an MP, and four quadcopters to perform the analysis. As such, the dynamic model of each can be developed in isolation, and then those can be integrated to define the dynamics of the CDPR as a whole. Naturally, different combinations can be created out of the same set of modules to generate multiple CDPRs.

Further, the inputs to the analyses can be either the variations of forces imparted by the actuators or the changes in the locations of the cables’ exit points with time. The latter class of problems is addressed for the first time in the present work to the best of the knowledge of the authors.

A recursive sub-system-level Lagrange multiplier (RSSLM) FD algorithm, originally introduced in ref. [2], is employed to efficiently perform the computations in linear time. This work extends its applicability to include closed-loop rigid-flexible multi-body systems with rheonomic constraints for the first time. In addition, the generic nature of this framework for analysing CDPRs involving simultaneous application of both kinds of actuation is also demonstrated for the first time.

A brief overview of the existing dynamic models of CDPRs is presented in Section 2. The modifications necessary to incorporate the Type II actuation are elucidated in Section 3. The effectiveness of the proposed framework is demonstrated with the help of a spatial $4$ - $4$ CDPR in Section 4. In particular, its capability to incorporate different types of actuation is illustrated with the same example in Section 4.4. Finally, the discussions and conclusions of the present work are presented in Sections 5 and 6, respectively.

2. Existing dynamic models of cable-Driven parallel robot

The dynamic model of CDPRs employed in the present work is an extension of those developed in refs. [2, 15]. Therefore, a brief overview of these models is presented first for the sake of completeness.

The cables are modelled as serial chains of modified rigid finite elements (MRFEs), mentioned in refs.[3]. As depicted in Fig. 1, each MRFE consists of two rigid links connected by a prismatic (P) joint, and these elements are connected to one another by a universal (U) joint. The axial deflections, considering both the axial stiffness and damping of the cables, are lumped together in the form of SDEs at the P joints. Similarly, those due to the transverse and lateral stiffness and damping properties are lumped together as the SDEs at the U joints.

Figure 1. A typical architecture of a modified rigid finite element, adopted from [29].

The MP was modelled as a rigid body connected in parallel to multiple serial chains of MRFEs. Effectively, the dynamic model of a CDPR without the inclusion of the effects of actuation was rendered equivalent to a rigid-flexible multi-body mechanical system with closed loops. The EoM of such systems are represented as follows:

(1) \begin{align}{\boldsymbol{{M}}}_j{\ddot{\boldsymbol{{q}}}} = \boldsymbol{{f}}_j + \sum _{i=1}^{u}\boldsymbol{{J}}_{{\boldsymbol{\eta }\boldsymbol{{q}}}_{i, j}}^\top \boldsymbol{\lambda }_i,\; j=1,\dots,v, \end{align}

where $v$ denotes the number of sub-systems, that is, cables and MP, and $u$ represents the number of independent closed loops. The matrix ${\boldsymbol{{M}}}_j \in GL_{n_j}({\mathbb{R}})$ is the $j$ th sub-system’s generalised inertia matrix (GIM) (see, e.g., [4]), with $n_j$ being the number of joints in that sub-system. The vector $\boldsymbol{{q}}_j \in{\mathbb{R}^{n_j}}$ represents the joint displacements, ${\dot{\boldsymbol{{q}}}}_j \in{\mathbb{R}^{n_j}}$ their velocities and ${\ddot{\boldsymbol{{q}}}}_j \in{\mathbb{R}^{n_j}}$ their accelerations. The first term on the right side of the equation, $\boldsymbol{{f}}_j \in{\mathbb{R}^{n_j}}$ , incorporates the contributions of the generalised forces imparted by the environment ( $\boldsymbol{\tau }_j$ ), due to the SDEs ( $\boldsymbol{\tau }_j^{\textrm{s}}$ ), the centripetal and Coriolis accelerations $\left (\sum _{k=1}^u{\boldsymbol{{C}}}_{j,k}{\dot{\boldsymbol{{q}}}}_k\right )$ and the gravity load ( $\boldsymbol{{G}}_j$ ). It is computed as the following sum:

(2) \begin{align} \boldsymbol{{f}}_j = \boldsymbol{\tau }_j +{\boldsymbol{\tau }_j^{\textrm{s}}} - \sum _{k=1}^v{\boldsymbol{{C}}}_{j, k}{\dot{\boldsymbol{{q}}}}_k - \boldsymbol{{G}}_j. \end{align}

The second term in Eq. (1) indicates the generalised constraint forces necessary for maintaining the connectivity between different sub-systems. The reaction forces at such common points between a pair of sub-systems are denoted by $\boldsymbol{\lambda }_i \in{\mathbb{R}^{m_i}}$ , with $m_i$ being the minimum number of constraints required to ensure the integrity of $i$ th loop. The distribution of reaction forces at every joint was computed by transforming these via the sub-matrices of constraint Jacobian matrix, $\boldsymbol{{J}}_{{\boldsymbol{\eta }\boldsymbol{{q}}}_{i, j}} \in{\mathbb{R}^{m_i \times n_j}}$ . These transformations are computed from the following constraint functions that need to be satisfied at all times to ensure the intactness of $i$ th closed loop:

(3) \begin{align}& \boldsymbol{\eta }_i (\boldsymbol{{q}}) = \textbf{0}, \end{align}

(4) \begin{align} \Rightarrow \; & \dot{\boldsymbol{\eta }}_i = \sum _{j=1}^{v}{\frac{\partial \boldsymbol{\eta }_i}{\partial \boldsymbol{{q}}_j }} \dot{\boldsymbol{{q}}}_j =\sum _{j=1}^{v} \boldsymbol{{J}}_{{\boldsymbol{\eta }\boldsymbol{{q}}}_{i, j}} \dot{\boldsymbol{{q}}}_j=\textbf{0}, \textrm{ where } \boldsymbol{{J}}_{{\boldsymbol{\eta }\boldsymbol{{q}}}_{i, j}} ={\frac{\partial \boldsymbol{\eta }_i}{\partial \boldsymbol{{q}}_j }}. \end{align}

The additional conditions required for the computation of the unknown reaction forces, $\boldsymbol{\lambda }_i$ , are obtained as:

(5) \begin{align} \ddot{\boldsymbol{\eta }}_i ={\frac{\textrm{d} \dot{\boldsymbol{\eta }}_i}{\textrm{d} t}} = \sum _{j=1}^{v} \left ( \boldsymbol{{J}}_{{\boldsymbol{\eta }\boldsymbol{{q}}}_{i, j}} \ddot{\boldsymbol{{q}}}_j + \dot{\boldsymbol{{J}}}_{{\boldsymbol{\eta }\boldsymbol{{q}}}_{i, j}} \dot{\boldsymbol{{q}}}_j \right ) = \textbf{0}. \end{align}

A computationally efficient RSSLM approach was employed to determine the reaction forces $\boldsymbol{\lambda }_i$ and the joint accelerations ${\ddot{\boldsymbol{{q}}}}_j$ at desired instances of time using Eq. (1) and Eq. (5) for specified joint positions $\boldsymbol{{q}}_j$ and velocities ${\dot{\boldsymbol{{q}}}}_j$ .

The incorporation of the Type I actuation through time-varying lengths of the rigid links of MRFEs, in ref. [15], induced an explicit dependence of every term in Eq. (1) on time, that is:

(6) \begin{align}{{\boldsymbol{{M}}}}_j (t){{\ddot{\boldsymbol{{q}}}}}_j + \sum _{k=1}^v{\boldsymbol{{C}}}_{j, k}(t){\dot{\boldsymbol{{q}}}}_k + \boldsymbol{{G}}_j (t) ={\boldsymbol{\tau }_j^{\textrm{s}}} (t) + \boldsymbol{\tau }_j +{\boldsymbol{\tau }_j^{\textrm{c}}}(t),\;j=1,\dots,v. \end{align}

Also, the matrix $\boldsymbol{{C}}$ represents generalised forces caused by the changing inertia of cables and those associated with centripetal and Coriolis accelerations. The linear time complexity of the RSSLM approach is retained even after such necessary modifications.

The modifications required to incorporate the Type II actuation in the same framework are presented in the next section.

3. Extension of the dynamic model to include Type II actuation

A typical scenario of the Type II actuation in CDPRs is illustrated in Fig. 2. Multiple driving mechanisms, denoted by $\textrm{D}_k$ , manipulate the positions of the exit points of cables, $\boldsymbol{{b}}_k$ . As mentioned in Section 1, these mechanisms could be aerial vehicles, mobile robots, or gantry cranes. Moreover, all the cables could be connected to the same rigid body, as in the case of a crane, or connected to multiple robots for cooperative transportation. This provision is made feasible due to the modular nature of the proposed framework, which can accommodate any other variation in the components of the system as well.

Figure 2. A schematic representation of a typical type II CDPR. The $k$ th cable, labelled as $C_k$ , is attached to the moving platform at $\xi _r$ and connected to the driving mechanism $\textrm{D}_k$ at $\boldsymbol{{b}}_k$ . The point $\boldsymbol{{b}}_k$ is known as the exit point (described in section 1).

The locations of the exit points of the cables are treated as immobile in refs. [2, 15]. However, to analyse the dynamics of Type II CDPRs, the motions of these points must be considered. First, the necessary modifications to the computational model presented in ref. [2] are elucidated. Next, the formulations of dynamics are characterised depending on the kinds of inputs supplied to the analyses, that is,

1. 1. Availability of the forces applied on the exit points or the actuator efforts of the driving mechanisms (attached to these points), and

2. 2. Provision of the trajectories tracked by the exit points of cables.

Finally, the efficacy of the proposed formulations is demonstrated using the incompletely restrained spatial suspended $4$ - $4$ CDPR as an example.

3.1. Forward dynamics formulation

The dynamic model of the cable is modified to accommodate the relative motions of their exit points with respect to the global frame of reference $\boldsymbol{{o}}$ - $X_0Y_0Z_0$ . Due to these additional DoF, the configuration of cable $C_k$ is updated to:

(7) \begin{align} \boldsymbol{{q}}_j = \left [\boldsymbol{{b}}_k^\top, q_1,\dots,q_{3{n_{\textrm{e}}}}\right ]^\top \in{\mathbb{R}^{3 ({n_{\textrm{e}}}+1)}},\; j=k=1,\dots,{n_{\textrm{k}}}. \end{align}

In Eq. (7), the symbol $\boldsymbol{{b}}_k$ denotes the position of the exit point of cable $C_k$ , whereas $q_\iota, \iota =1,\dots,3{n_{\textrm{e}}},$ with $n_{\textrm{e}}$ being the number of elements, represent the configuration of the cable. The number of cables is denoted by $n_{\textrm{k}}$ . Further changes to the formulation of dynamics based on the nature of inputs provided to the analyses are outlined below.

3.1.1. Inputs in the form of forces applied at the exit points of cables

The first kind of inputs considered in the following is the forces exerted by the driving mechanisms on the cables. In such cases, assuming no other external forces are acting on the cables, the vector of external generalised forces of cable $C_k$ has the following structure:

(8) \begin{align} \boldsymbol{\tau }_k = \left [\boldsymbol{\lambda }_i^\top, 0, \dots, 0\right ]^\top \in{\mathbb{R}^{3 ({n_{\textrm{e}}}+1)}},\; k=1,\dots,{n_{\textrm{k}}},\; i={n_{\textrm{k}}}+ k. \end{align}

The force $\boldsymbol{\lambda }_{{n_{\textrm{k}}}+ k}$ in Eq. (8) acts at the exit point of cable $C_k$ .

Alternatively, the generalised forces internal to the actuators of the driving mechanisms, denoted by $\boldsymbol{\tau }_i$ , could be specified as inputs to the analyses when there is no provision for direct measurement of $\boldsymbol{\lambda }_i$ . In such scenarios, as depicted in Fig. 3, every mechanism $\textrm{D}_k,\,k=1,\dots,{n_{\textrm{k}}}$ , is treated as a separate sub-system of the CDPR to compute the required restraining forces at the corresponding exit point.

Figure 3. Driving mechanisms, $\textrm{D}_k,\,k=1,\dots,{n_{\textrm{k}}}$ , as separate sub-systems of a Type II CDPR. The symbol $\boldsymbol{\tau }_i, i={n_{\textrm{k}}}+k,$ denotes the internal actuator forces of the $k$ th driving mechanism and $\boldsymbol{\lambda }_i$ represents the reaction forces at the exit point $\boldsymbol{{b}}_k$ of the cable $C_k$ .

At first, the connectivity between the sub-systems $\textrm{D}_k$ and $C_k$ is ensured by introducing the action-reaction pair of forces, denoted by $\boldsymbol{\lambda }_i$ . Subsequently, $\boldsymbol{\lambda }_i$ are computed from the imposition of additional constraint equations, $\boldsymbol{\eta }_i=\textbf{0}$ , where:

(9) \begin{align} \boldsymbol{\eta }_i = \boldsymbol{{b}}_k -{\boldsymbol{{e}}}_k,\; k=1,\dots,{n_{\textrm{k}}},\,i={n_{\textrm{k}}}+k. \end{align}

In Eq. (9), the location of the leading end of cable $C_k$ is denoted by $\boldsymbol{{b}}_k$ and its attachment point on $\textrm{D}_k$ by ${\boldsymbol{{e}}}_k$ .

In summary, the above formulation can handle inputs either in the form of forces directly applied on the exit points of cables or those that are internal to the actuators of the driving mechanisms. However, there can be many situations where the inputs are available in the form of motions of the exit points rather than the forces acting on them. One instance of such a scenario can be where the exit points of the cables are being carried by individual aerial vehicles, such as quadrotors or helicopters – instead of incorporating their dynamics into the overall model, one can simply track their motions and use those as kinematic inputs to the overall system. Analyses of such systems are discussed next.

3.1.2. Inputs in the form of trajectories tracked by the exit points of cables

The second kind of inputs to the analyses are the trajectories of the exit points of cables. As shown in Fig. 4, the point $\boldsymbol{{b}}_k$ of cable $C_k$ is considered to follow the input trajectory ${\boldsymbol{{e}}}_k(t)$ . Consequently, the inputs at some of the joints will be kinematic while the rest are dynamic in nature. Therefore, simulation of the CDPR requires a hybrid approach, that is, both the FD and ID analyses as classified in [30], pp. $171$ – $172$ . Such problems of CDPRs are addressed for the first time in the present work, to the best of the authors’ knowledge.

Figure 4. Representative trajectories $e_k(t)$ of the exit points $\boldsymbol{{b}}_k$ of cables $C_k$ in Type II CDPRs. The reaction forces acting at the exit point $\boldsymbol{{b}}_k$ are denoted by $\boldsymbol{\lambda }_{{n_{\textrm{k}}}+k}$ .

In the present work, the hybrid dynamic analysis of the robot is circumvented by introducing the following constraint equations:

(10) \begin{align} \boldsymbol{\eta }_i(\boldsymbol{{b}}_k, t) \,:\!=\, \boldsymbol{{b}}_k -{\boldsymbol{{e}}}_k(t) = 0,\; k=1,\dots,{n_{\textrm{k}}},\, i={n_{\textrm{k}}}+k. \end{align}

Clearly, Eq. (10) represents rheonomic constraints, ³ that is, they depend explicitly on time. Therefore, the associated Lagrange multipliers, denoted by $\boldsymbol{\lambda }_i$ , are also functions of time. Specifically, differentiating $\boldsymbol{\eta }_i$ twice results in the following expression:

(11) \begin{align}{\frac{\textrm{d} \boldsymbol{\eta }_i}{\textrm{d} t}} &={\frac{\textrm{d} }{\textrm{d} t}} \left ({\frac{\partial \boldsymbol{\eta }_i}{\partial t }} + \sum _{j=1}^v \boldsymbol{{J}}_{\boldsymbol{\eta }\boldsymbol{{q}}_{i,j}}{\dot{\boldsymbol{{q}}}}_j\right ) \nonumber \\ &={\frac{\partial ^2 \boldsymbol{\eta }_i}{\partial t^2}} + \sum _{j=1}^v \left ({\frac{\partial }{\partial q_j }} \left ({\frac{\partial \boldsymbol{\eta }_i}{\partial t }}\right ){\dot{\boldsymbol{{q}}}}_j + \boldsymbol{{J}}_{\boldsymbol{\eta }\boldsymbol{{q}}_{i,j}}{\ddot{\boldsymbol{{q}}}}_j +{\frac{\textrm{d} }{\textrm{d} t}}\left (\boldsymbol{{J}}_{\boldsymbol{\eta }\boldsymbol{{q}}_{i,j}}\right ){\dot{\boldsymbol{{q}}}}_j \right ) . \end{align}

Equation (11) may be simplified significantly by taking into cognisance the following.

1. 1. The variations of $\boldsymbol{\eta }_i$ w.r.t. time and the configuration variables, $\boldsymbol{{q}}_j$ , are mutually independent, hence ${\frac{\partial }{\partial q_j }}\left ({\frac{\partial \boldsymbol{\eta }_i}{\partial t }}\right ) ={\textbf{0}}$ .

2. 2. The matrix $\boldsymbol{{J}}_{\boldsymbol{\eta }\boldsymbol{{q}}_{i,j}}$ is of constant nature, since its entries are either ones or zeros. Therefore, ${\frac{\textrm{d} }{\textrm{d} t}}\left (\boldsymbol{{J}}_{\boldsymbol{\eta }\boldsymbol{{q}}_{i,j}}\right ) ={\textbf{0}}$ .

With these, Eq. (11) reduces to:

(12) \begin{align} &{\frac{\textrm{d}^2 \boldsymbol{\eta }_i}{\textrm{d} t^2}} ={\frac{\partial ^2 \boldsymbol{\eta }_i}{\partial t^2}} + \sum _{j=1}^v \boldsymbol{{J}}_{\boldsymbol{\eta }\boldsymbol{{q}}_{i,j}}{\ddot{\boldsymbol{{q}}}}_j ={\textbf{0}},\; i={n_{\textrm{k}}}+1,\dots,2{n_{\textrm{k}}}. \end{align}

The key difference between the above form of the loop-closure equation and that associated with the scleronomic types (as in Eq. (5)) is the additional term $\frac{\partial ^2 \boldsymbol{\eta }_i}{\partial t^2}$ , accounting for the explicit dependence of $\boldsymbol{\eta }_i$ on time. Consequently, the expression of the residual constraint accelerations $\boldsymbol{\psi }_i$ (see Eq. (15)) is modified to:

(13) \begin{align} \boldsymbol{\psi }_i ={\frac{\partial ^2 \boldsymbol{\eta }_i}{\partial t^2}} + \sum _{j=1}^v \boldsymbol{{J}}_{\boldsymbol{\eta }\boldsymbol{{q}}_{i,j}} \boldsymbol{{a}}_j,\; i={n_{\textrm{k}}}+1,\dots,2{n_{\textrm{k}}}. \end{align}

In Eq. (13), the symbol $\boldsymbol{{a}}_j \in{\mathbb{R}^{n_j}}$ denotes the unconstrained joint accelerations of the $j$ th sub-system and it is given by:

(14) \begin{align} \boldsymbol{{a}}_j ={{\boldsymbol{{M}}}_j}^{-1} \boldsymbol{{f}}_j,\;j=1,\dots,v, \end{align}

where $\boldsymbol{{f}}_j$ is the cumulative effect of the generalised forces devoid of the inertial and constrained forces acting on the $j$ th sub-system, as in Eq. (2). Finally, the Lagrange multipliers can be computed using the same expressions reported in ref. [2], that is:

(15) \begin{align} \sum _{s=1}^{u} \boldsymbol{{A}}_{i, s} \boldsymbol{\lambda }_s + \boldsymbol{\psi }_i = \textbf{0},\;i=1,\dots,u,\textrm{ and }u= 2{n_{\textrm{k}}}, \end{align}

where $\boldsymbol{{A}}_{i, s} \in{\mathbb{R}^{m_i \times m_i}}$ is the sub-matrix of the generalised inertia constraint matrix, mentioned in ref. [31]. Effectively, by virtue of the generalisations mentioned above, the modified RSSLM approach can simulate the dynamics of Type II CDPRs even when the trajectories of the driving mechanisms are provided as inputs.

4. Illustrative examples

The utility of the modified computational model is demonstrated on an incompletely restrained suspended $4$ - $4$ CDPR. The details related to the architecture of the robot and the notations used in the formulation of dynamics are delineated in Section 4.1. Also, a rudimentary model used for detecting the contact between the MP and the ground so as to prevent the MP from penetrating the ground is described in Appendix B. In the first case study, the transfer of load from an initial pose to a final one based on the input trajectories of the exit points of cables is demonstrated. In the second one, four quadcopters are used to manipulate the motions of these exit points. The response of the CDPR when one of those quadcopters fails mid-flight is examined. In the final case study, the response of the robot is investigated when different types of actuation, namely, Types I & II, are simultaneously applied.

4.1. Architecture of a 4-4 cable-driven parallel robot

A $4$ - $4$ CDPR consists of four cables $C_1,\dots,C_4$ connected to the MP at $\boldsymbol{\xi }_1,\dots,\boldsymbol{\xi }_4,$ and to the driving mechanisms at $\boldsymbol{{b}}_1,\dots,\boldsymbol{{b}}_4$ . It is decomposed into four cables and an MP, as shown in Fig. 5. The length, width and height of the MP are represented by ${l_{\textrm{m}}},{w_{\textrm{m}}},{h_{\textrm{m}}},$ respectively. Since the local frame of reference $\boldsymbol{\xi }_{\textrm{c}}$ - $X_aY_aZ_a$ is considered identical to the fixed one $\boldsymbol{{o}}$ - $X_0Y_0Z_0$ at the initial configuration, the ZYZ convention of Euler angles used for representing its orientation in the previous works, for example, [2], will lead to a parametric singularity (see, e.g., [32], pp. $51$ – $52$ ). Therefore, the XYZ convention is used in this case to circumvent the singularity, and the corresponding DH parameters are listed in Table I.

Figure 5. Notional decomposition of the $4$ - $4$ cable-driven parallel robot into four cables and a moving platform.

The action-reaction forces at the disassembled joints are determined via the computation of the Lagrange multipliers $\boldsymbol{\lambda }_1,\dots,\boldsymbol{\lambda }_4$ associated with the loop-closure constraint functions in Eq. (3). The dimension of the configuration space of every cable is given by ${n_{\textrm{k}}}=3({n_{\textrm{e}}}+1), k=1,\dots,4,$ and hence, that of the $4$ - $4$ CDPR is $n=12{n_{\textrm{e}}}+18$ . The number of constraint functions is $m=24$ , that is, three per each end of the cable.

4.2. Pick-and-place operation by the 4-4 cable-driven parallel robot

The numerical values of the architecture and inertia parameters of the robot are listed in Table II. Further, the values of the various tolerances used in the simulation are the same as the ones reported in Table C.6 of [2]. The unstrained lengths of the cables are:

(16) \begin{align} l_k=0.75 \textrm{ m}, k=1,\dots,4. \end{align}

Since these unstrained lengths are small and do not vary with time, only five MRFEs are used for modelling each cable; thereby, ${n_{\textrm{e}}}=5$ . Subsequently, the coefficients of stiffness of the SDEs are determined using the expressions reported in ref. [2]. Their numerical values are given by:

(17) \begin{align} &{s_{\textrm{a}}} = 1.16 \times 10^7\,\textrm{N/m},\;{s_{\textrm{t}}} ={s_{\textrm{l}}}= 11.64\,\textrm{Nm/rad}. \end{align}

In Eq. (17), the subscripts a, t and l denote the values associated with deflections in the axial, transverse and lateral directions, respectively. Likewise, the calculated damping coefficients of cables are as follows:

(18) \begin{align}{d_{\textrm{a}}} = 2.40 \times 10^{5}\,\textrm{Ns/m},\;{d_{\textrm{t}}}={d_{\textrm{l}}}= 0.24\,\textrm{Nms/rad}. \end{align}

The initial pose of the robot is shown in Fig. 6. As mentioned before, the initial configuration of the MP is selected such that it rests on the ground with its local frame of reference $\boldsymbol{\xi }_{\textrm{c}}$ - $X_aY_aZ_a$ oriented in the same direction as the fixed frame of reference $\boldsymbol{{o}}$ - $X_0Y_0Z_0$ . Therefore, its position and orientation are given by:

(19) \begin{align} \boldsymbol{{q}}_5 = \left [0.70, 1.00, 0.10, 0, 1.57, 0\right ]^\top . \end{align}

Next, the initial configurations of the cables are chosen to be taut and oriented along the vertical direction. Accordingly, their numerical values are obtained as follows:

(20) \begin{align} &\boldsymbol{{q}}_1=\left [0.40, 0.60, 0.95, 0, 3.14, 0, \dots \right ]^\top, \boldsymbol{{q}}_2=\left [1.00, 0.60, 0.95, 0, 3.14, 0, \dots \right ]^\top, \nonumber \\ &\boldsymbol{{q}}_3=\left [0.40, 1.40, 0.95, 0, 3.14, 0, \dots \right ]^\top, \boldsymbol{{q}}_4=\left [1.00, 1.40, 0.95, 0, 3.14, 0, \dots \right ]^\top . \end{align}

Finally, the trajectories of the exit points of cables are planned as described below.

Figure 6. Initial configuration of the $4$ - $4$ cable-driven parallel robot with its moving platform resting on the ground, that is, the plane $X_0Y_0$ .

At first, the path followed by every point $\boldsymbol{{b}}_k, k=1,\dots,4,$ is divided into four line segments. As depicted in Fig. 7, the first segment, $L_1$ , is meant to vertically lift the MP from its initial location; the second, $L_2$ , is meant to spatially ascend it; the third, $L_3$ , is meant to spatially descend it; and the last, $L_4$ , is meant to vertically lower it to the desired destination on the ground. Next, the trajectories along each segment are defined via cubic polynomials of time, assuming that the points $\boldsymbol{{b}}_k$ start and end with zero velocities.

Figure 7. Specified identical path followed by the exit points of the cables $\boldsymbol{{b}}_k$ , $k=1,\dots,4$ . It comprises four line segments, $L_1$ to $L_4$ .

The input trajectories of the exit points are represented using the difference $\Delta{\boldsymbol{{e}}}_k(t) ={\boldsymbol{{e}}}_k(t)-{\boldsymbol{{e}}}_k(0),\,k=1,\dots,4$ . The components of the vector $\Delta{\boldsymbol{{e}}}_k=\left [\Delta e_{k\textrm{x}}, \Delta e_{k\textrm{y}}, \Delta e_{k\textrm{z}}\right ]^\top$ are defined as the following piece-wise functions:

(21) \begin{align} \Delta e_{k\textrm{x}} & = \Delta e_{k\textrm{y}} = \left \{ \begin{array}{l@{\quad}l} 0, & 0 \leq t \leq 2,\\[5pt] \dfrac{(t-2)^2(7-2t)}{2\sqrt{2}}, & 2 \lt t \leq 3,\\[15pt] \dfrac{82 + t (t(21 - 2 t) -72)}{2\sqrt{2}}, & 3 \lt t \leq 4,\\[9pt] \dfrac{1}{\sqrt{2}}, & 4 \lt t \leq 6. \end{array}\right ., \nonumber \\[10pt] \Delta e_{k\textrm{z}} &= \left \{ \begin{array}{l@{\quad}l} \dfrac{t^2(3-t)}{4}, & 0 \leq t \leq 2,\\[9pt] \dfrac{2-\sqrt{3} (t-2)^2 (2t-7)}{2}, & 2 \lt t \leq 3,\\[9pt] \dfrac{2+\sqrt{3}+\sqrt{3} (t-3)^2 (2t-9)}{2}, & 3 \lt t \leq 4,\\[9pt] \dfrac{(t-6)^2 (t-3)}{4}, & 4 \lt t \leq 6; \end{array}\right .\;k=1,\dots,4. \end{align}

The variations of $\Delta{\boldsymbol{{e}}}_k$ with time are depicted in Fig. 8. A video file named “motionOf44 cdprN5_ipTrjExPts.mp4” depicts the motion of the robot for the specified input.

Figure 8. Input trajectory of the exit points of the cables $\Delta{\boldsymbol{{e}}}_k ={\boldsymbol{{e}}}_k(t) -{\boldsymbol{{e}}}_k(0)$ , $k=1,\dots,4$ . The legends $\Delta e_{k\textrm{x}}$ , $\Delta e_{k\textrm{y}}$ , $\Delta e_{k\textrm{z}}$ represent the vector components of $\Delta{\boldsymbol{{e}}}_k$ . Due to the symmetry in the chosen path, $\Delta e_{k\textrm{x}}$ is identical to $\Delta e_{k\textrm{y}}$ .

The displacements of the MP for the input trajectories of $\boldsymbol{{b}}_k$ are shown in Fig. 9, using the difference $\Delta \boldsymbol{{q}}_5 (t)= \boldsymbol{{q}}_5 (t)-\boldsymbol{{q}}_5(0)$ . In addition, the variations of its instantaneous velocities and accelerations with time are depicted in Figs. 10a, 10b and Figs. 10c, 10d, respectively. Apart from the small yaw motion of the MP initiated at the end of vertical lift, that is, at $t=2$ s, there is no significant change in the orientation of the MP (see Fig. 9b).

Figure 9. Variation in the configuration of the moving platform of the $4$ - $4$ cable-driven parallel robot, $\Delta \boldsymbol{{q}}_5(t) = \boldsymbol{{q}}_5(t)-\boldsymbol{{q}}_5(0)$ , for the inputs given in Eq. (21).

Figure 10. Variations in the linear, angular velocities $\dot{\boldsymbol{{q}}}_5(t)$ and accelerations $\ddot{\boldsymbol{{q}}}_5(t)$ of the moving platform of the $4$ - $4$ cable-driven parallel robot, corresponding to the inputs given in Eq. (21).

Further, the variations of the linear displacements of the MP (shown in Fig. 9a) are similar to those of the exit points of the cables (see Fig. 8). Ideally, these variations would be the same if the dynamics of the robot were ignored. Therefore, an error measure $\delta{\boldsymbol{\xi }_{\textrm{c}}}(t)$ used to quantify such effects is defined as:

(22) \begin{align} \delta{\boldsymbol{\xi }_{\textrm{c}}}(t) = \Delta{\boldsymbol{\xi }_{\textrm{c}}}(t)- \Delta{\boldsymbol{{e}}}_k(t). \end{align}

The changes in $\delta{\boldsymbol{\xi }_{\textrm{c}}} (t)$ with time are depicted in Fig. 11. As is evident from the results, a kinematic estimation of the displacements of the MP will be erroneous. In this case, the maximum deviations are given by ${\| \delta{\boldsymbol{\xi }_{\textrm{c}}}\|}_\infty = \left [4.46, 3.30, 0.16\right ]^\top \times 10^{-2}$ m, which are approximately $\left [10, 8, 0\right ]^\top \%$ relative error in the respective directions. Similarly, the errors in estimating the linear velocities are ${\| \delta{\dot{\boldsymbol{\xi }}_{\textrm{c}}}\|}_\infty =\left [0.34, 0.32, 0.02\right ]^\top$ m/s and the linear accelerations are ${\| \delta{\ddot{\boldsymbol{\xi }}_{\textrm{c}}}\|}_\infty =\left [3.37, 4.90, 0.32\right ]^\top$ m/s ${}^2$ . Therefore, the dynamics of the cables and the MP should be included in the analyses for better accuracy.

Figure 11. Variations in the error $\delta{\boldsymbol{\xi }_{\textrm{c}}} = \Delta{\boldsymbol{\xi }_{\textrm{c}}}- \Delta \boldsymbol{{b}}_k(t)$ with time. The scalar components of $\delta{\boldsymbol{\xi }_{\textrm{c}}}$ are denoted by $\delta \xi _{\textrm{x}}$ , $\delta \xi _{\textrm{y}}$ and $\delta \xi _{\textrm{z}}$ .

The magnitudes of the axial extensions of the cables are minimal ( $O\left (10^{-6}\right )$ m) throughout the simulation of the robot. The changes in the constraint forces at the anchoring points of the cables $\boldsymbol{{c}}_k$ are shown in Fig. 12, and those corresponding to the exit points $\boldsymbol{{b}}_k$ are depicted in Fig. 13. All the cables pull the MP vertically for the entire duration of the simulation because ${\lambda }_{i_{\textrm{z}}}(t) \lt 0,\;\forall i=1,\dots,4,$ and ${\lambda }_{i_{\textrm{z}}}(t) \gt 0,\;\forall i=5,\dots,8$ . The former indicates that the cables counter the weight and vertical inertial forces of the MP, while the latter implies that the driving mechanisms are being pulled during the operation. Hence, these reaction forces complement each other, as seen by the change in sign of the respective plots. Also, the abrupt change in $|{\lambda }_{i_{\textrm{z}}}|$ in the vicinity of $t=6$ s is due to the impact of the MP on the ground.

Figure 12. Variations in the values of components of the reaction forces $\boldsymbol{\lambda }_i,\,i=1,\dots,4,$ with time, corresponding to the simulation of the $4$ - $4$ cable-driven parallel robot for the inputs given in Eq. (21).

Figure 13. Variations in the values of components of the reaction forces $\boldsymbol{\lambda }_i,\,i=5,\dots,8,$ with time, corresponding to the simulation of the $4$ - $4$ cable-driven parallel robot for the inputs given in Eq. (21).

The computational time ⁴ taken to perform the simulation with the input trajectories ${\boldsymbol{{e}}}_k(t)$ of the points $\boldsymbol{{b}}_k$ in Eq. (21) is $2261$ s (approximately, $37$ minutes). Similar to the results presented in ref. [15], the solver increases the temporal resolution to accurately determine the robot’s response to the transitions, as the exit points switch segments in their paths. However, as depicted in Fig. 14, the majority of the time is spent in capturing the dynamics of the robot while spatially ascending and descending the MP. This could be attributed to the initiation and the rapid variations in transverse and lateral deflections of cables, which is also evident from the MP’s displacements shown in Fig. 9 from $t=2$ s to $5$ s.

Figure 14. Variations in the step size $\Delta t$ used by the solver ode15s at every instance of the simulation of the $4$ - $4$ CDPR for the input trajectories of the exit points of the cables $\boldsymbol{{b}}_k$ in Eq. (21). The sustenance of smaller steps can be seen during the spatial ascend and descend of the MP, that is, segments $L_2$ and $L_3$ in Fig. 7.

4.2.1. Validation of the simulation results

The satisfaction of the imposed kinematic constraints is verified by employing the error measures $e_1$ and $e_2$ defined in Appendix E. The variations of errors $e_1,\,e_2$ of the constraints responsible for the connectivity between the cables and the MP are depicted in Fig. 15. Similarly, those responsible for tracking the input trajectory of the exit points of cables are shown in Fig. 16. In both cases, the residues are within respective desirable limits for the entire duration of the simulation. Hence, the obtained results are confirmed to be valid.

Figure 15. Variations in the errors $e_1$ and $e_2$ with time, corresponding to the constraint functions associated with the connectivity between the cables and the MP. The definitions of $e_1$ and $e_2$ are given in Eqs. (E1) and (E2), respectively.

Figure 16. Variations in the errors $e_1$ and $e_2$ with time, corresponding to the constraint functions responsible for incorporating the input trajectories of the exit points of cables. The definitions of $e_1$ and $e_2$ are given in Eqs. (E1) and (E2), respectively.

4.3. Dynamic response of the robot in the event of mid-flight failure in the actuation

In this example, the driving mechanisms connected to the $4$ - $4$ CDPR are chosen to be four identical quadcopters. First, the dynamic model of the quadcopter is described. Next, the temporal changes in the forces internal to the quadcopters are devised. Finally, the response of the robot starting from the configuration shown in Fig. 6 is presented.

The architecture of the quadcopter that drives the cable $C_k$ is shown in Fig. 17. A local frame of reference, ${\boldsymbol{{e}}}_k$ - $X_{\textrm{m}_k}Y_{\textrm{m}_k}Z_{\textrm{m}_k}$ , is attached at its centre of mass. The orientation of this frame with respect to the inertial frame of reference $\boldsymbol{{o}}$ - $X_0Y_0Z_0$ is represented using the XYZ convention of Euler angles ${\alpha }_k,\beta _k,\gamma _k$ . Therefore, the configuration of the $k$ th quadcopter is defined by $\boldsymbol{{q}}_j=\left [\boldsymbol{{e}}_k^\top,{\alpha }_k,\beta _k,\gamma _k\right ]^\top$ , $j=k+5$ . The DH parameters required for modelling these sub-systems are listed in Table I.

Figure 17. Schematic of the model of quadcopter. The thrust forces and reactive moments of the $k$ th quadcopter are denoted by $f_{1_k},\dots,f_{4_k}$ and $m_{1_k},\dots,m_{4_k}$ , respectively. The reaction forces from the $k$ th cable is denoted by $\boldsymbol{\lambda }_{{n_{\textrm{k}}}+k}$ .

Further, the thrust forces and reactive moments of the rotors are directed along the $Z_{\textrm{m}_k}$ -axis at all times. The magnitudes of the former are denoted by $f_{1_k}, \dots, f_{4_k}$ , while those of the latter are represented by $m_{1_k}, \dots, m_{4_k}$ . Moreover, the ratio of these moments to the forces, denoted by $c_{\textrm{tm}}$ , is constant for specific propeller designs (see, e.g., [33]). Therefore, only four of these eight magnitudes can be varied independently. Furthermore, drag forces are neglected ⁵ in the present work.

Accordingly, the equations of motion of the quadcopters can be expressed in the form of Eq. (1). Subsequently, the external forces, $\boldsymbol{\tau }_j$ , are computed using the expression:

(23) \begin{align} &\boldsymbol{\tau }_j = \left [f_{\textrm{t}_k}\sin{\beta _k}, -f_{\textrm{t}_k}\sin{\alpha _k}\cos{\beta _k}, f_{\textrm{t}_k}\cos{\alpha _k}\cos{\beta _k}, \left(f_{2_k}-f_{4_{k}}\right)r_k, \left(f_{3_k}-f_{1_{k}}\right)r_k,m_{\textrm{t}_k}\right ]^\top, \nonumber \\[3pt] &\textrm{where } f_{\textrm{t}_k}= f_{1_k}+\dots +f_{4_k},\,m_{\textrm{t}_k} = m_{1_k}+\dots +m_{4_k},\;k=1,\dots,4,\;j=k+5. \end{align}

In Eq. (23), the symbol $r_k$ denotes the offset of the rotors’ axes of rotation from the point $\boldsymbol{{e}}_k$ . Since there are no SDEs associated with these sub-systems, the vector $\boldsymbol{\tau }_j^{\textrm{s}}$ is null. In addition, the sub-matrices of the constraint Jacobian matrix ${\boldsymbol{{J}}}_{\boldsymbol{\eta }\boldsymbol{{q}}}$ associated with the constraint functions in Eq. (9) are constant. Therefore, there is no need for numerical estimation of these matrices, and their entries are either ones or zeroes. Furthermore, the sparsity of the GIMs ${\boldsymbol{{M}}}_j,\,j=6,\dots,9,$ is taken into account while computing their inverses. The numerical values of the architecture and inertia parameters of the quadcopters used in the simulation are listed in Table III.

The input actuator forces of the quadcopters are planned as follows. First, the yaw motion of quadcopters is prevented by imposing the following conditions on the reactive moments of the rotors:

(24) \begin{align} m_{1_k}=-m_{3_k},\,m_{2_k}=-m_{4_k},\;k=1,\dots,4. \end{align}

Next, the changes in the thrust forces for the four quadcopters are considered as:

(25) \begin{align} &f_{1_k} = f_{2_k} = f_{3_k} = f_{4_k} = f_k,\;k=1,\dots,4,\nonumber \\[10pt] &f_k = \left \{ \begin{array}{l l} \dfrac{5}{2}\sin{\left (\dfrac{\pi }{2} t\right )}, & 0 \lt t \leq \dfrac{1}{2},\\[10pt] \dfrac{5}{2}, & \dfrac{1}{2}\lt t \leq 3,\\[7pt] 0, & k=1, t \gt 3, \\[5pt] \dfrac{5}{2}, & k \neq 1, t \gt 3. \end{array}\right . \end{align}

Their variations with time are shown in Fig. 18. Finally, the corresponding values of the reactive moments are computed using the constant $c_{\textrm{tf}}$ . A video file named “motionOf44cdprN5_ipForces_vertInitCond.mp4” depicts the motion of the robot for the specified input when started from the configuration shown in Fig. 6.

Figure 18. Temporal variations of the input thrust forces of the quadcopters, $f_{1_k} = f_{2_k} = f_{3_k} = f_{4_k} = f_k$ , $k=1,\dots,4$ .

The difference $\Delta \boldsymbol{{q}}_5 = \boldsymbol{{q}}_5 (t)-\boldsymbol{{q}}_5(0)$ is used for studying the linear and angular displacements of the MP due to the input forces mentioned above. The variations of $\Delta \boldsymbol{{q}}_5$ with time are depicted in Fig. 19. Initially, the collective forces exerted by the quadcopters are not sufficient to lift the MP. However, after $t=0.48$ s, the MP starts to move up vertically, that is, $\Delta \xi _{\textrm{z}} \gt 0$ (see Fig. 19a). Subsequently, the MP accelerates in that direction with no significant change in its orientation. Even though the actuators of the first quadcopter are shut off at $t=3$ s to simulate its failure, the MP continues to ascend due to inertia of motion. However, its orientation no longer remains the same as the initial one (see Fig. 19b). Eventually, due to insufficient support from the remaining three quadcopters, the MP starts to descend at $t=4.21$ s. The corresponding changes in the velocities and accelerations of the MP are depicted in Figs. 20a, 20b and Figs. 20c, 20d, respectively.

Figure 19. Variation in the configuration of the moving platform of the $4$ - $4$ cable-driven parallel robot, $\Delta \boldsymbol{{q}}_5(t) = \boldsymbol{{q}}_5(t)-\boldsymbol{{q}}_5(0)$ , for the inputs given in Eq. (25).

Figure 20. Variations in the linear, angular velocities $\dot{\boldsymbol{{q}}}_5(t)$ and accelerations $\ddot{\boldsymbol{{q}}}_5(t)$ of the moving platform of the $4$ - $4$ cable-driven parallel robot, corresponding to the inputs given in Eq. (25).

Further, the variations in the reaction forces at the trailing ends of the cables are shown in Fig. 21. Similarly, the changes in those at the exit points of the cables are depicted in Fig. 22. As is evident from the plots, the variations of reaction forces in different cables remain the same until $t=3$ s. Moreover, the cables start to pull the MP vertically, that is, ${\lambda }_{i_z}\lt 0,\,i=1,\dots,4,$ only after the applied actuator forces reach a specific limit at $t=0.31$ s. Also, due to the inertia of rest of the MP, slightly larger forces are required to initiate the lift. Hence, an abrupt decrease in the magnitude of reaction forces $\left |{\lambda }_{i_z}\right |$ can be seen after $t=0.48$ s. Furthermore, the first cable becomes non-supportive and acts as additional load on the remaining quadcopters after $t=3$ s.

Figure 21. Variations in the values of components of the reaction forces $\boldsymbol{\lambda }_i,\,i=1,\dots,4,$ with time, corresponding to the simulation of the $4$ - $4$ cable-driven parallel robot for the inputs given in Eq. (25).

Figure 22. Variations in the values of components of the reaction forces $\boldsymbol{\lambda }_i,\,i=5,\dots,8,$ with time, corresponding to the simulation of the $4$ - $4$ cable-driven parallel robot for the inputs given in Eq. (25).

Figure 23. Variations in the time steps $\Delta t$ of the solver ode15s at every instance of the simulation of the $4$ - $4$ CDPR for the changes in the actuator forces of the quadcopters in Eq. (25). The sustenance of smaller steps can be seen after the failure of the first quadcopter after $t=3$ s. Such a failure renders the cable to be non-supportive and leads to disruptions in the motions of the remaining cables and quadcopters.

There is no observable difference in the orientation of the quadcopters. However, the changes in their linear displacements are as depicted in Fig. 24. Although the actuator inputs of the first quadcopter cease to exist after $t=3$ s, its motion is governed by the reaction forces $-\boldsymbol{\lambda }_5$ , shown in Fig. 22a.

Figure 24. Variations in the linear displacements, $\Delta \boldsymbol{{b}}_k = \boldsymbol{{b}}_k(t)-\boldsymbol{{b}}_k(0)$ , of the quadcopters, $\textrm{D}_k$ , $k=1,\dots,4,$ for the inputs specified in Eq. (25).

Similar to the results presented in Section 4.2, the deflections in the cables due to their axial compliance are minimal throughout the simulation. The computational time taken to obtain the response of $4$ - $4$ CDPR for the inputs given in Eq. (25) is $1.38 \times 10^{4}$ s. It takes only $7.39$ s to obtain results for $t\in [0, 3)$ s. Therefore, as seen in Fig. 23, most of the time is spent obtaining the results for $t \in [3, 6]$ s, where the cables undergo transverse and lateral displacements.

Furthermore, even with the same input forces as in the previous one, when the simulation is initiated using a different configuration of cables, such as the one shown in Fig. 25, it took only $121$ s of computational time. This difference can be attributed to gravity assisting the cables to deflect in the transverse direction when the quadcopters fail to apply sufficient pulling forces on the cables. In the former, due to the alignment of cables with the direction of gravity, that is, the vertical axis $\boldsymbol{{b}}_o$ - $Z_0$ (shown in Fig. 6), both gravity and the quadcopters induce compressive loads. A video file named “motionOf44cdprN5_ipForces_incInitCond.mp4” depicts the motion of the robot for the specified input when started from the configuration shown in Fig. 25.

Figure 25. Initial configuration of the $4$ - $4$ cable-driven parallel robot with its moving platform resting on the ground, that is, the plane $X_0Y_0$ . Also, in contrast to the configuration depicted in Fig. 6, the cables are not aligned along the direction of gravity, the axis $Z_0$ .

4.3.1. Validation of the simulation results

The values of the errors $e_1$ and $e_2$ defined in Appendix E, are computed for the entire duration of the simulation. The variations in these errors associated with the constraint functions $\boldsymbol{\eta }_i,\,i=1,\dots,4,$ are shown in Fig. 26, and those related to $\boldsymbol{\eta }_i,\,i=5,\dots,8,$ are depicted in Fig. 27. In both cases, the errors are within desirable limits. Therefore, the connectivity of the cables with the MP and the quadcopters is verified to be intact during the simulation of the robot.

Figure 26. Variations in the errors $e_1$ and $e_2$ with time, corresponding to the constraint functions associated with the connectivity between the cables and the MP. The definitions of $e_1$ and $e_2$ are given in Eqs. (E1) and (E2), respectively.

Figure 27. Variations in the errors $e_1$ and $e_2$ with time, corresponding to the constraint functions responsible for incorporating the input trajectories of the forces of the quadcopters. The definitions of $e_1$ and $e_2$ are given in Eqs. (E1) and (E2), respectively.

4.4. Pick-and-place operation with sinusoidal feed rates of the cables

The simulation is initiated from the configuration of the robot depicted in Fig. 25. In addition to the input trajectories of the exit points defined in Eq. (21), the input feed rates of the cables are provided using the following expression:

(26) \begin{align} \dot{l}_k = \frac{{\omega _{\textrm{c}}}}{10} \cos{\left ({\omega _{\textrm{c}}}t + (k-1)\frac{\pi }{10}\right )},\;k=1,\dots,4,\;{\omega _{\textrm{c}}}=\frac{\pi }{5}. \end{align}

The associated input trajectories of the cables’ lengths are depicted in Fig. 28. Accordingly, the masses, moments of inertia and locations of the centres-of-mass of rigid links of MRFEs are updated at every time instance. The stiffness and damping coefficients are recomputed. The necessary changes to the formulation of dynamics noted in ref. [15] are introduced.

Figure 28. Variations in the input unstrained lengths of cables of the 4-4 CDPR, $l_k,\,k=1,\dots,4$ , associated with the feed rates given in Eq. (26).

The changes in the position and orientation of the MP are represented using the difference $\Delta \boldsymbol{{q}}_5$ in Figs. 29a, 29b, respectively. The associated variations in the velocity and acceleration of the MP are shown in Figs. 30a, 30b and Figs. 30c, 30d, respectively. As opposed to the results presented in Section 4.2, significant changes in the MP’s angular displacements can be observed in this case due to the time-varying lengths of cables.

Figure 29. Variations in the configuration of the moving platform of the 4-4 cable-driven parallel robot, $\Delta \boldsymbol{{q}}_5= \boldsymbol{{q}}_5(t)- \boldsymbol{{q}}_5(0)$ , with time for the inputs given in Eqs. (21) and (26).

Figure 30. Variations in the linear, angular velocities $\dot{\boldsymbol{{q}}}_5(t)$ and accelerations $\ddot{\boldsymbol{{q}}}_5(t)$ of the moving platform of the $4$ - $4$ cable-driven parallel robot, corresponding to the inputs given in Eqs. (21) and (26).

Figure 31. Variations in the errors $e_1$ and $e_2$ with time, corresponding to the constraint functions associated with the connectivity between the cables and the MP. The definitions of $e_1$ and $e_2$ are mentioned in Eq. (E1) and Eq. (E2), respectively.

As is the case with the previous case studies, the axial extensions of the cables are minimal for the entire duration of the simulation. The variations in the forces exerted by the MP on the cables at their anchoring points are shown in Fig. 32. The forces responsible for the desired motion of the exit points of cables are depicted in Fig. 34. Evidently, all the cables pull the MP during the operation of the robot. However, the load is not shared uniformly amongst the cables due to the temporal changes in their lengths. A depiction of the response of the CDPR to simultaneous application of feeding of cables and the movement of their exit points is included as a video file named “motionOf44cdprN5_diffMoA.mp4.”

Figure 32. Variations in the values of components of the reaction forces $\boldsymbol{\lambda }_i,\,i=1,\dots,4,$ with time, corresponding to the simulation of the $4$ - $4$ cable-driven parallel robot for the inputs given in Eqs. (25) and (26).

Figure 33. Variations in the errors $e_1$ and $e_2$ with time, corresponding to the constraint functions that ensure the trajectories of the exit points of the cables to be the same as the inputs. The definitions of $e_1$ and $e_2$ are mentioned in Eq. (E1) and Eq. (E2), respectively.

Figure 34. Variations in the values of components of the reaction forces $\boldsymbol{\lambda }_i,\,i=5,\dots,8,$ with time, corresponding to the simulation of the $4$ - $4$ cable-driven parallel robot for the inputs given in Eqs. (25) and (26).

4.4.1. Validation of the simulation results

The dynamic evolution of the robot is in accordance with the imposed constraints. The residues of the constraint functions $\boldsymbol{\eta }_i,\,i=1,\dots,4,$ responsible for the intactness of the connectivity between the cables and the MP are shown in Fig. 31 and those $\boldsymbol{\eta }_i,\,i=5,\dots,8,$ indicating the accuracy in tracking the desired trajectory of the exit points are depicted in Fig. 33.

5. Discussions

A summary of the present work and its key contributions vis-à-vis the existing literature are presented in the following.

The present work is primarily distinguished from the previous ones in its focus on the modelling and simulation of the dynamics of the Type II CDPRs, since the control-related aspects seem to dominate the literature. To that end, a new approach (reported in ref. [2]) with efficient recursive algorithms is utilised. Also, two different kinds of inputs were considered in the analyses. In cases where the trajectories of the exit points of cables are specified, the hybrid dynamic analyses of Type II CDPRs are addressed by incorporating additional rheonomic constraints in the formulation of dynamics.

The hybrid dynamics of serial and tree-type architectures have been extensively studied in the past (see, e.g., [34]). Such approaches may not be easily extended to parallel architectures due to the additional constraint forces associated with the combination of different sub-systems. In addition, the ID analyses of Type II CDPRs are not straightforward due to the associated underactuation in cables and actuators, as explained in ref. [35]. Therefore, additional rheonomic constraints are used in the present study to avoid the ID analysis. As mentioned before, such a class of problems is addressed for the first time, to the best of the authors’ knowledge, in the present work.

The second distinction of the proposed computational framework is in the employed model of the cables. As mentioned in Section 1, most of the previous studies on the Type II CDPRs assume the cables to be massless, inelastic and taut. It was justified based on the scale of applications being generally small. However, such approximations do not hold true in the case of aggressive manoeuvres or when the mass of the cable is comparable to that of the MP or when there are low tensions in the cables. Therefore, a few studies have considered slack models of cables, as outlined in Section 1.

Nonetheless, the slack models utilised in the prior works are incomplete, that is, the stiffness and damping properties of the cables were ignored while treating these as serial chains of rigid bodies connected by U joints in refs. [20, 35]. Recently, lumped masses connected by SDEs were employed to closely approximate the physics of cables; see, for example, [21]. Nevertheless, the stiffness and damping properties in the transverse and lateral directions were not incorporated in these works. Naturally, these distinctions continue to hold good even when the application domain is extended to include simultaneous actuation of Type I and II since no prior work has been found to employ a comprehensive dynamic model of the cables in such a scenario.

The scope of application of the RSSLM approach is also extended to the FD analyses of rigid-flexible multi-body autonomous systems with rheonomic constraints, as opposed to those with only scleronomic constraints in the previous works (see, e.g., [2, 15]). This aspect of the proposed framework, aided by the modular nature of the associated mathematical formulation, affords a fast expansion of the scope of such analyses. In this article, the CDPRs with Type I and Type II actuations are discussed. However, since the internal dynamics of the mechanisms responsible for the motion of the exit points of the cables may be completely abstracted out in favour of their final kinematic outputs (i.e., their motion trajectories), the very nature of such mechanisms becomes completely irrelevant to the analyses. Consequently, in addition to the FD analyses of several variants of Type I and Type II CDPRs, the proposed framework can also be applied to many situations in different domains, such as air-to-air refuelling of aircrafts, underway replenishments and the feeding and reeling-in of towed (individual) cables and ropes from a ship or a submarine.

In summary, the present work brings together a number of theoretical and computational elements, such as rheonomic constraints; a comprehensive MRFE model including inertia, elasticity and damping; and the computationally efficient RSSLM approach, in the context of analysing the dynamics of Type II CDPRs. Furthermore, it demonstrates simultaneous application of both types of actuation for the first time, to the best of the authors’ knowledge. These developments constitute key enablers for broadening the scope of analysis for similar complex systems including rigid and flexible elements, without having to make too many simplifying assumptions.

The proposed framework can be further improved in the following manner. It is assumed that there are no collisions between the cables or the MP. Such idealisations may not hold true in the case of redundantly restrained CDPRs or cluttered environments. Therefore, suitable collision models are required to analyse the contact dynamics of the cables and the MP for close proximities among these or with the surrounding obstacles. As investigated in ref. [8], the contact models become a necessity in analysing the cooperative towing of a payload using mobile robots, an example of the planar Type II CDPR. In such scenarios, the payload is always in contact with the ground and slides on it.

6. Conclusions

The computational model presented in the prior works has been updated to include the Type II actuation, wherein the locations of the exit points of cables are varied with time. Hereof, input trajectories tracked by the exit points or the external forces applied at these points are incorporated into the model. The former is achieved via the addition of rheonomic constraint equations. In comparison, the latter is demonstrated to involve minimal changes in the dynamic model when the forces acting at the exit points are provided. However, if such reaction forces are not available, then the dynamics of the mechanisms driving these points are analysed along with the cables and the MP. Subsequently, the efficacy of the improved computational model is exemplified by simulations of the $4$ - $4$ CDPR. The modularity and generality of the proposed work to simulate the dynamics of CDPRs when both types of actuation are applied simultaneously is demonstrated using the same example. Finally, the results are validated by checking the connectivity of the cables with the MP and the actuator mechanisms for the entire duration of the simulations. The codes developed for this purpose are available for download via the link: https://github.com/TejaKrishnaMamidi/rsslm.