Wrench Uncertainty Quantification and Reconfiguration Analysis in Loosely Interconnected Cooperative Systems

Cable-based parallel architectures are of significant interest due to the advantages offered over their rigid-link counterparts. Replacing articulated arms by light-weight cable elements results in simpler construction, lower cost, simpler transportation, larger workspace, ease of reconfiguration, and lower interaction of the robot structure with the environment. However, one consequence of using cable elements is the increase of system complexity due to the unilateral force constraint of the cables. Hence, in order to keep the robot fully restrained, $m=n+1$ cables are at least required for a cable-robot with n degrees-of-freedom. Nonetheless, in order to improve the wrench-closure-workspace [1], many of the cable-based systems are redundantly restrained (⁠$m>n+1$⁠). Hence, cooperation of several actuation units (e.g., pulley and cables, and gantry cranes) is a characteristic feature of such systems.

In a cooperative task, individual agents are required to either maintain a predefined position/orientation or follow a desired trajectory. However, the exact pose of each robot remains uncertain due to control error. Loose interconnections between robots (e.g., aerial towing system) and flexibility of the system structure, particularly in cable-robots with elastic cables, increase the spatial fluctuations of the agents (that cannot be fully compensated by the control system) and hence exacerbate the uncertainty in the configuration of the entire cooperative system. Moreover, the exerted force by an individual agent is not exactly known due to the inherent actuation uncertainty. The problem of configuration and actuation uncertainty becomes further critical as there exists an increasing interest in using several low-cost off-the-shelf robots in designing the multi-agent cooperative systems to reduce the total cost. Such low-cost agents are less precise and introduce more uncertainty that propagates to the output of the entire cooperative system.

While loose interconnections and base mobility in the new class of cable-robots increase the configuration uncertainty, it adds reconfigurability to the system and motivates developing reconfiguration schemes to lower the uncertainty in the system output. We show that using the redundancy, we can achieve this by reconfiguring the system, while the desired mean output is maintained for any new configuration. Hence, a systematic and efficient quantification of the system response uncertainty, given its configuration, is necessary. Efficiency in uncertainty quantification is crucial because, as we will discuss in Sec. 4, measures of the output uncertainty are constructing the objective function and constraints of the configuration optimization problem. The significance of the computational efficiency is further increased when the reconfiguration is aimed to be performed in real time.

A motivational example is discussed here to support the necessity and effectiveness of the uncertainty quantification and reconfiguration scheme in mobile cable-robots. Consider a cooperative payload transportation using flying robots (aerial towing), shown in Fig. 1(a), as an illustrative example. Cooperative systems consisting of aerial robots have received a significant amount of attention in the last decade. The low-cost off-the-shelf flying robots, such as quadrotors, are easily accessible and multiple number of them can be employed for different cooperative tasks. Several studies particularly focus on different aspects of the design and control of such systems. For example, see Refs. [2–5].

In many such applications, the system design seeks to compensate the gravitational force and move the payload in an almost quasi-static condition. Low-cost flying robots are inherently uncertain, hence, neither the load provided by individual agents nor their spatial position can be identified deterministically. Consequently, the resultant wrench delivered at the platform is random due to the uncertainty in tension and orientation of the suspending cables. The main idea is that although the uncertainty introduced by each individual agent cannot be reduced, the overall effect of these uncertainties can be minimized by reconfiguration of these flying units. Clearly, redistribution of the tension on the cables is simultaneously considered in order to maintain a prescribed mean output wrench.

Figures 1(b) and 1(c) illustrate this concept where the objective is to transport the payload from point A to point B without collision with the walls. In Fig. 1(b), the output uncertainty, corresponding to the configuration of the agents, is shown by a gray ellipse. The configuration in Fig. 1(b) can be considered unsafe for this particular task as the major axis of the uncertainty ellipse crosses the walls. Hence, reconfiguration needs to be performed such that the probability of collision is minimized as shown in Fig. 1(c), where the major axis of the uncertainty ellipse is aligned with the straight path from A to B. Note that in addition to re-orienting the ellipse, the uncertainty can also be reduced, i.e., the length of the minor and major axes can be minimized.

While the discussed example is a spatial system, the focus of the current work is on the planar systems. Besides the direct application to planar systems (e.g., see the system in Ref. [6]), the proposed approach in the current work can be used to treat the spatial systems by appropriate mapping into two-dimensional (2D) space such that the quantity of the interest can be mainly captured. For instance, in the example shown in Fig. 1, the out-of-plane uncertainty, i.e., uncertainty in the direction of the gravity, may be excluded from design criteria, and hence, a 2D formulation is adequate for the design of the reconfiguration schemes. Moreover, current study sets the ground for extension to three-dimensional (3D) wrench uncertainty analysis [7].

Representation and statistical manipulation of the uncertainty have been the subject of several studies in robotics. Early contributions in this area dates back to the work by Smith and Cheeseman [8] on the spatial uncertainty defined on the coordinate frames that was applied to the navigation problem of the mobile robots. Followed by their work and in order to adopt the convenient representation of the motion, quantification of the uncertainty in the form of homogeneous transforms and on the motion groups received special interest. Su and Lee [9] used the differential homogeneous transform to represent the uncertainty. They addressed the kinematic error propagation in a robotic manipulator used in an assembly task. Error propagation on the Euclidean motion group (with a focus on robotic manipulator applications) has been comprehensively studied by Wang and Chirikjian [10–13]. In their work, the propagation of uncertainty is first studied by convolution on SE(3) and, using the theory of Lie algebras and Lie groups, first- [10,11] and second-order [12,13] approximations to the mean and covariance were provided. This enabled the recursive approximation of the mean and covariance of the whole manipulator. The robustness of the second-order approximation as well as its superiority compared to the first-order method was verified in different numerical examples. Sovizi et al. [14–18] proposed an uncertainty characterization method in complex robotic manipulators based on the random matrix theory (RMT). It was shown that RMT-based models can appropriately capture the uncertainty only using limited information about the system variation. It is beneficial in many cases of complex architectures where the detailed statistical information are not available. Moreover, the proposed RMT-based model was able to adaptively tune the uncertainty level depending on the state/configuration of the system.

Cable-based cooperative systems, which will mostly benefit from the approach developed in this paper, have been extensively considered in many recent studies. In towing applications, the cooperative agents are commonly mobile robots such as quadrotors for aerial towing and wheeled robots for ground towing. Cheng et al. [6] proposed a planar system in which a payload was towed/carried with multiple mobile robots (Scarab) attached to the platform through cables. A comprehensive treatment to the aerial towing systems consisting of multiple flying robots has been addressed by Refs. [2,3], and [19–22]. It can be seen that the performance of the cooperative system has been degraded by inevitable uncertainty in the flying robot position and provided force. For example in Ref. [3], fluctuations in controlled position of the individual robots as well as random oscilliations of the platform around equilibrium configurations have been reported.

Base mobility and reconfiguration problem also applies to the new class of the cable-based manipulators. Although conventional cable-robots use fixed bases with varying cable lengths [23–28], recent studies [29–35] focus on reconfigurable systems where the cable lay out can be changed through mobility of the bases. The configuration optimization has been considered to improve some performance indices such as dexterity [30], tension factor [32,36], and power consumption and stiffness [34,35]. However, in this work, our main focus is to reduce the uncertainty of the generated wrench at the platform that directly impacts the platform position uncertainty.

An overview of the literature reveals that (i) new class of reconfigurable cable-based systems are increasingly evolving and extending their application to a broad range of physical problems; (ii) these systems are significantly impacted by the uncertainty especially with recent growing attitudes toward using low-cost and rapidly prototyped robotic agents in many distributed systems; and (iii) systematic treatment to the uncertainty in these systems and approaches toward reducing the system output fluctuations have not been addressed in the literature.

To this end, in this paper, we provide a probabilistic framework where the forces applied to the platform by each individual agent are treated as random vectors. The statistics of the output wrench are characterized through solving the problem of sum of random vectors. Using the closed-form expressions of the means and variances of the output wrench, we cast the reconfiguration problem into a polynomial optimization that can be efficiently solved in real-time applications. Using both numerical and real-system experiments, we show the effectiveness of the proposed approach in realizing the configurations that reduce the uncertainty in certain directions. It is noteworthy that in addition to reconfiguration planning, our quantification of the output wrench uncertainty provides useful statistical information for higher-level control and facilitates a systematic design trade-off between agents with different levels of precision.

The rest of this paper is organized as follows: Section 2 reviews the preliminaries from multivariate statistics used in this paper for our analytical analysis. The wrench uncertainty formulation and closed-form expressions of the means and variances are developed in Sec. 3. Reconfiguration scheme and the procedure to formulate and solve the optimization problem are developed in Sec. 4. Section 5 provides numerical example for optimizing the location of the bases in a static system. Furthermore, our experimental demonstration of the proposed approach is reported in Sec. 6. Finally, a brief discussion and direction of the future work are given in Sec. 7.

Some basics from probability theory and directional statistics that are used in this paper are reviewed here.

where $sk∈{1,…,n}$ (⁠$∀k=1,…,r$⁠). From the general representation in Eq. (2), the elements of the mean (first moment) vector and second moment matrix can be obtained as $μs1=E[xs1]$ and $Σ̃s1,s2=E[xs1xs2]$⁠, respectively. Similarly, the elements of the covariance matrix are $Σs1,s2=E[(xs1−μs1)(xs2−μs2)]$⁠.

In fact, from Eqs. (4) and (5), the first cumulant is simply the first moment (mean) and second cumulant is the second central moment (covariance).

in which $ci=((∂KX(η,ζ))/∂xi)|(0,0) and ci,j=((∂2KX(η,ζ))/∂xi∂xj)|(0,0)$⁠, where $KX(η,ζ)=log MX(η,ζ)=log E[exp {ηx1+ζx2}]$⁠. We refer the reader to Ref. [37] for a quick reference on the multivariate moments and cumulants.

where $θ$ is the angle of random (unit) vector u in 2D coordinate system.

Different studies have addressed the maximum likelihood estimation (MLE) of the parameters of the von Mises-Fisher distribution, i.e., $u¯$ and λ_u [38–40]. To keep the paper as self-contained as possible, we briefly review the MLE of the parameters of the von Mises distribution for circular data.

Refer to Refs. [38] and [39] for additional details on von Mises-Fisher distribution and estimation of its parameters.

Using vector analysis, the total wrench at the platform of a planar system shown in Fig. 2 can be written as the sum of individual wrenches as

where $Fx$ and $Fy$ are the random output forces along x- and y-axes, respectively, and M is the random output moment. Also, $fix$ and $fiy$ are the random forces applied by the ith agent along x- and y-axes, respectively, and $mi$ is the random moment induced by ith agent. $θi=θ¯i+θiν$ and $Ti=T¯i+Tiν$ are the random orientation and tension of the ith cable in which $θ¯i$ and $T¯i$ are the mean orientation and tension, and $θiν$ and $Tiν$ are the angle and tension perturbations, respectively. $ri=rixê1+riyê2$ is the vector in global frame x–y pointing from the origin of the local frame u–v to the attachment point of the ith cable, as shown in Fig. 2. In this study, we characterize up to the second statistical moments of the output force vector as well as output moment at the platform.

and $DTi$ and $Dθi$ are the domains where $Ti$ and $θi$ are defined, respectively. As Gaussian models reliably represent the uncertain quantities in several physical systems, we assume the tension in each cable is normally distributed with mean $T¯i$ and standard deviation $σTi$⁠, denoted as $Ti∼N(T¯i,σTi2)$⁠. Furthermore, the cable orientation is assumed to have von Mises distribution with mean $θ¯i$ and dispersion parameter $λθi$⁠, denoted as $θi∼M(θ¯i,λθi)$⁠. The validity of these assumptions are further investigated in Sec. 6 using experimental data.

Note that the first equalities in Eqs. (21) and (22) are due to Eq. (6) and first equalities in Eqs. (23) and (24) are due to Eq. (7).

where $z(Ti)=ai2+bi2$ in which $ai=ηriTi+λθi cos(θ¯i)$ and $bi=ηviTi+λθi sin(θ¯i)$⁠.

There are many available packages that can be used to solve such polynomial optimization problems such as GloptiPoly [41] (used in this paper) and matlab optimization toolbox (e.g., fmincon function along with sqp algorithm).

We first validate the closed-form expressions of the statistical moments given by Eqs. (21)–(24), (27), and (28) by a series of the MC simulations for different number of agents (cables) in the system. Then, examples of planar cable-robots are considered to validate the effectiveness of the proposed approach in quantifying the wrench uncertainty and efficient optimization of the system state.

Three planar systems with m = 4, m = 5, and m = 6 cables are considered in order to validate the closed-form expressions of the statistical moments. The geometrical and statistical information corresponding to each system is shown in Fig. 3.

Table 1 compares the means and variances of the resultant force and moment given by Eqs. (21)–(24), (27), and (28) with the corresponding values obtained through 10,000 MC simulations.

The results tabulated in Table 1 verify the validity of the developed analytical solutions to the statistical moments of the output wrench. Given the fact that Eqs. (21)–(24), (27), and (28) are the exact solutions, the sample means and variances obtained from MC simulations will get closer to those given by Eqs. (21)–(24), (27), and (28) by increasing the number of realizations.

Using a laptop computer with 2.2 GHz central processing unit and 16 GB memory, the running time for MC-based calculation of wrench statistics (with 10,000 samples) for m = 6 was 1.92 s. This indicates that such MC-based simulations are not applicable to the real-time optimization problems where wrench statistics need to be computed at each iteration (thousands of iterations may be needed before finding a minimum). However, using Eqs. (21)–(24), (27), and (28), the running time reduced to 0.0094 s enabling a real-time reconfiguration of the system to reduce the uncertainty of the output wrench.

Here, we inspect the proposed reconfiguration strategy in a static cable-based architecture that provides a static wrench at the platform (end effector). The bases are kept fixed during the task; however, they can adaptively reconfigure to improve the performance index (i.e., uncertainty level here) for each specific task. The real-system experimental proof of the concept and validation of the simulation results reported in this section is presented in Sec. 6.

In the first set of our examples, the desired output moment is set to be zero. So, all cables are attached to one point on the platform. In optimization of the tension and orientation of the cables, we examine three objective functions. First is the tension factor given by $TF=min(T¯)/max(T¯)$⁠. Maximizing the TF provides an improved tension balance among the cables and it has been shown [36] that it is equivalent to minimizing $JTF=∑i=1mT¯i$⁠. So, in Eq. (36), we have $J=JTF=∑i=1mzi$⁠. Second and third objective functions are $J=Var(Fx)$ (by choosing $γ1=1, γ2=0$⁠, and $γ3=0$⁠) and $J=Var(Fy)$ (by choosing $γ1=0, γ2=1$⁠, and $γ3=0$⁠), respectively. For all optimization problems, we set $Wd=[45 N, 60 N, 0 N⋅m]T, λθi=100$⁠, and $σTi=2 N ∀i=1,…,m$⁠. The lower and upper bounds on the mean tension of the cables are assumed to be $T¯Li=0 N$ and $T¯Ui=100 N ∀i=1,…,m$⁠. Assuming a minimum tension of $0 N$ allows the algorithm to remove a cable if the optimal configuration can be achieved by fewer number of cables. Moreover, $θ¯i$'s can be chosen anywhere on the circle, i.e., $θ¯i∈[0, 2π]$⁠. This may be not feasible in real systems (e.g., iRobot cooperative system described in Sec. 6), where different agents must maintain a certain distance to avoid collision. However, for simulation purpose and to identify minimum possible uncertainty level, a bound is not prescribed for the cable angles in this example. For the next example where the output moment uncertainty is considered and for the real-system experiment, described in Sec. 6, $θ¯i$'s are restricted by certain bounds.

For the optimal $T¯i$'s and $θ¯i$'s corresponding to each objective function, we generate 10,000 samples of the output force vector by sampling $Ti$'s from Gaussian distribution, given by Eq. (17), and $θi$'s from von Mises distribution, given by Eq. (9), and substituting the samples into Eq. (12). A pdf is approximated from the generated samples using (bivariate) kernel smoothing (KS) density approximation technique. Figure 4 shows the optimal configuration and cable tensions as well as the resulting KS joint pdf of the output force corresponding to each optimization problem. In the configuration plots, the acceptable range of the $θ¯i$'s are shown by dashed circular arcs and the (optimal) location of the bases are shown by diamond markers. In KS joint density plots, hot colors show higher densities and lower densities are shown by cold colors. Figure 4(a) corresponds to a system with m = 2 cables and the objective is to maximize TF. Second row (Figs. 4(b) and 4(c)) shows the results when the objective is solely minimizing $Var(Fx)$⁠. Comparison of Fig. 4(b) with Fig. 4(a) indicates that the uncertainty in $Fx$ cannot be notably reduced, when m = 2. However, increasing the number of cables facilitates a significant reduction in the $Var(Fx)$⁠, where it reduces from $20.68 N2$ in Fig. 4(a) to $9.95 N2$ in Fig. 4(c). Third row of Fig. 4 illustrates the results when the objective is to reduce $Var(Fy)$⁠. Similarly, reduction in uncertainty level of $Fy$ is not significant when only two cables are used; however, a major reduction from $15.16 N2$ (Fig. 4(a)) to $7.6 N2$ (Fig. 4(e)) can be achieved in a system with six cables. Roughly speaking, and referring to Figs. 4(c) and 4(e), by optimizing the tension and orientation of the cables, we are, in fact, rotating the uncertainty ellipse to align its minor axis with a direction along which the maximum precision is required. Furthermore, the results show that in addition to reorienting the ellipse, the length of axes can also be reduced.

The uncertainty in output moment can be similarly treated as the output force vector. Figure 5 shows the results corresponding to TF maximization (first row) and minimization of the $Var(M)$ (second row). The platform (gray circles) is chosen to be circular for simplicity in the presentation of the results; however, it can be of any general geometrical shape. There are ten possible attachment points on a circular arc with a radius $2 m$ whose center is matching with the center of the circular platform. We examine two systems with m = 3 and m = 4 cables and the optimization problem is solved for $(10m)$ choices of the attachment points to find the one that provides the least uncertainty level in output moment.

For the results shown in Fig. 5, we have $Wd=[45 N, 60 N,120 N⋅m]T, λθi=100$⁠, and $σTi=2 N ∀i=1,…,m$⁠. Additionally, the mean tension of the cables are restricted to $T¯Li=0 N$ and $T¯Ui=40 N ∀i=1,…,m$⁠. The domain from which $θ¯i$'s can be chosen are shown by dashed arcs whose lengths are $π rad$⁠. As shown in Fig. 5, the variance of the output moment can be reduced from $56.04 N⋅m2$ (when TF is maximized) to $39.02 N⋅m2$ (when $Var(M)$ is minimized) in a system with m = 3 cables. Similarly, when m = 4, $Var(M)$ is reduced from $61.91 (N⋅m)2$ to $37.09 (N⋅m)2$⁠. It is also noteworthy that the attachment points can be chosen depending on the objective/application of the system. For example, referring to Figs. 5(b) and 5(d), the algorithm chooses a different set of attachment points when the objective is changed to minimize the output moment variance.

We apply our method to a cooperative system shown in Fig. 6. In our experiment, the desired output moment is considered to be 0 and, therefore, all cables are attached to one point on the force transducer. Three mobile robots used in this study are iRobot Create kits that communicate (using a serial port) with a laptop computer to receive the velocity command which determines the tension in the cables. An ATI Industrial Automation six-axis force/torque transducer is used to measure the tension provided by each robot and when all three robots are in use. OptiTrack cameras and reflective markers (not shown in Fig. 6) are used for initialization and reconfiguration of the robots.

We choose (nominally) identical robots; however, their characteristics are different and hence need to be identified individually. The quantities of interest to construct the optimization problem (described by Eqs. (36)–(43)) are $σTi, λθi, T¯iL$⁠,and $T¯iU$ for $i=1,…,3$⁠. To find an estimate of these quantities, we send a 50 s long motion command to the iRobot that is tied to the force transducer by a cable. We measure the force in both x- and y-directions. This process is repeated for two more motion commands that provide more tractions. There is a pause between the motion commands in order to facilitate segmenting the data points.

Figure 7 shows the magnitude of the force vector versus time for three iRobots used in our experiment. The mean of the traction samples for each motion command is shown above the corresponding data points. We use these mean values to estimate $T¯iL$ and $T¯iU$⁠. Moreover, we center the entire data points (traction values) around zero traction by subtracting the corresponding mean from each set of the data points. Then, a sample estimate of $σTi$ is obtained to be used in our optimization problem.

The force vector orientation can be represented by the random unit vector u with von Mises-Fisher distribution given by Eq. (8) defined on the unit sphere $Sn−1={u∈ℝn:||u||=1}$⁠, where in planar system, n = 2, and $S1$ is a unit circle. Let $Uj={u1j,…,uNj}$ be the set of N samples of u corresponding to the jth motion command. It can be seen that the data sets U^j (⁠$j=1,…,3$⁠) have different sample mean vectors, as shown in Fig. 8. However, in order to use the MLE approach presented in Sec. 2.2 for estimating $λθi$⁠, all data points must share the same mean and dispersion parameter. We use a theorem discussed in Ref. [40] to address this problem as following. An arbitrary mean vector is first chosen. Here, we choose the north pole, i.e., $u¯np=[0, 1]T$⁠. Using Rodrigues' formula [42], we calculate a rotation matrix R_j that rotates the sample mean vector $u¯j$ (corresponding to data set U^j) to $u¯np$⁠, i.e., $u¯np=Rju¯j$⁠. Then, the data set $U$ is constructed as follows: