Loosely interconnected cooperative systems such as cable robots are particularly susceptible to uncertainty. Such uncertainty is exacerbated by addition of the base mobility to realize reconfigurability within the system. However, it also sets the ground for predictive base reconfiguration in order to reduce the uncertainty level in system response. To this end, in this paper, we systematically quantify the output wrench uncertainty based on which a base reconfiguration scheme is proposed to reduce the uncertainty level for a given task (uncertainty manipulation). Variations in the tension and orientation of the cables are considered as the primary sources of the uncertainty responsible for nondeterministic wrench output on the platform. For nonoptimal designs/configurations, this may require complex control structures or lead to system instability. The force vector corresponding to each agent (e.g., pulley and cable) is modeled as random vector whose magnitude and orientation are modeled as random variables with Gaussian and von Mises distributions, respectively. In a probabilistic framework, we develop the closed-form expressions of the means and variances of the output force and moment given the current state (tension and orientation of the cables) of the system. This is intended to enable the designer to efficiently characterize an optimal configuration (location) of the bases in order to reduce the overall wrench fluctuations for a specific task. Numerical simulations as well as real experiments with multiple iRobots are performed to demonstrate the effectiveness of the proposed approach.

## Introduction

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.

### Motivational Examples.

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].

### Literature Review.

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.

## Preliminaries

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

### Multivariate Moments and Cumulants.

**x**,

*r*th (statistical) moment is defined as

**x**and $Dx$ is the domain where

**x**is defined. The first moment is the mean (i.e., $E[x]=x\xaf$) and the variance is the second central moment defined by $\sigma x2=E[(x\u2212x\xaf)2]$. Now, let $X=[x1,\u2026,xn]T$ be an

*n*-dimensional random vector. The components of the

*r*th moment of the random vector

**X**, denoted by $ms1,s2,\u2026,sr$, can be written as

where $sk\u2208{1,\u2026,n}$ ($\u2200k=1,\u2026,r$). From the general representation in Eq. (2), the elements of the mean (first moment) vector and second moment matrix can be obtained as $\mu s1=E[xs1]$ and $\Sigma \u0303s1,s2=E[xs1xs2]$, respectively. Similarly, the elements of the covariance matrix are $\Sigma s1,s2=E[(xs1\u2212\mu s1)(xs2\u2212\mu s2)]$.

**x**, mgf is defined as $Mx(\eta )=E[exp\u2009{\eta x}]$. The multivariate counterpart, and for random vector

**X**, is the joint moment generating function (jmgf) defined as $MX(v)=E[exp\u2009{vTX}]$. Moments of the distribution can be obtained by differentiating jmgf. However, analytical analyses are sometimes simpler when using joint cumulant generating function (jcgf). Given the jmgf, the jcgf is $KX(v)=log\u2009MX(v)$. The coefficients of the Taylor expansion of the jcgf around the origin are, in fact, the cumulants. From the Taylor expansion of the multivariate functions, it can be seen that the

*r*th cumulant, denoted by $cs1,s2,\u2026,sr$, is

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).

*n*= 2, we have $X=[x1,x2]T$ and $v=[\eta ,\zeta ]T$. The mean vector and covariance matrix can be written in terms of the cumulants as

in which $ci=((\u2202KX(\eta ,\zeta ))/\u2202xi)|(0,0)\u2009\u2009and\u2009ci,j=((\u22022KX(\eta ,\zeta ))/\u2202xi\u2202xj)|(0,0)$, where $KX(\eta ,\zeta )=log\u2009MX(\eta ,\zeta )=log\u2009E[exp\u2009{\eta x1+\zeta x2}]$. We refer the reader to Ref. [37] for a quick reference on the multivariate moments and cumulants.

### von Mises Distribution and Maximum Likelihood Parameter Estimation.

*n*-dimensional unit vector

**u**, defined on $(n\u22121)$-dimensional hypersphere (sphere when

*n*= 3 and circle when

*n*= 2), has von Mises-Fisher distribution with mean vector $u\xaf$ and dispersion parameter

*λ*, denoted by $u\u223cvMF(u\xaf,\lambda u)$, if its pdf is given by

_{u}*ν*th-order ($Re(\nu )>\u22121/2$) modified Bessel function of the first kind. For planar case when

*n*= 2, the von Mises-Fisher distribution reduces to the von Mises distribution whose density function is given by

where $\theta $ 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\xaf$ 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.

*N*observations of the random variable $\theta $ (angles of the unit vectors $u1,\u2026,uN$) that has von Mises distribution with mean $\theta \xaf$ and dispersion parameter $\lambda \theta $, denoted as $\theta \u223cM(\theta \xaf,\lambda \theta )$. It can be shown that the MLE of $\theta \xaf$ is simply the direction of the sum of the sample unit vectors, i.e., $\theta \xaf\u0302=angle(u1+\cdots +uN)$. Furthermore, it can be shown that the likelihood function is maximized by

## Wrench Uncertainty

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 *i*th agent along *x*- and *y*-axes, respectively, and $mi$ is the random moment induced by *i*th agent. $\theta i=\theta \xafi+\theta i\nu $ and $Ti=T\xafi+Ti\nu $ are the random orientation and tension of the *i*th cable in which $\theta \xafi$ and $T\xafi$ are the mean orientation and tension, and $\theta i\nu $ and $Ti\nu $ are the angle and tension perturbations, respectively. $ri=rixe\u03021+riye\u03022$ is the vector in global frame *x–y* pointing from the origin of the local frame *u–v* to the attachment point of the *i*th 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.

### Output Force Vector Statistics.

and $DTi$ and $D\theta i$ are the domains where $Ti$ and $\theta 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\xafi$ and standard deviation $\sigma Ti$, denoted as $Ti\u223cN(T\xafi,\sigma Ti2)$. Furthermore, the cable orientation is assumed to have von Mises distribution with mean $\theta \xafi$ and dispersion parameter $\lambda \theta i$, denoted as $\theta i\u223cM(\theta \xafi,\lambda \theta i)$. The validity of these assumptions are further investigated in Sec. 6 using experimental data.

### Output Moment Statistics.

where $z(Ti)=ai2+bi2$ in which $ai=\eta riTi+\lambda \theta i\u2009cos(\theta \xafi)$ and $bi=\eta viTi+\lambda \theta i\u2009sin(\theta \xafi)$.

## Base Configuration and Cable Tension Optimization

*J*is the objective function that quantifies the output wrench uncertainty. Constraints described by Eqs. (30)–(32) ensure that the mean of resultant wrench is sufficiently close to the prescribed wrench. Superscript

*d*indicates the desired (prescribed) values and $eFx,\u2009eFy,\u2009and\u2002eM$ are the maximum acceptable distances between mean and the desired values of the wrench vector elements. Inequality constraint (33) ensures the kinematic feasibility of the cable orientations (typically to avoid the collision) and (34) guarantees that optimal tensions are consistent with the force capacity of the agents. Depending on application scenario, different variant of optimization problems can be formulated by selecting various combinations of objective functions. One such formulation based on the assumption of linear combination of force and moment variances can be written as follows:

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).

## Numerical Examples

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.

### MC Verification.

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.

m = 4 | m = 5 | m = 6 | ||
---|---|---|---|---|

$E[Fx]$ | 32.05 | 39.46 | 41.17 | Eq. (21) |

32.09 | 39.55 | 41.11 | MC | |

$E[Fy]$ | −7.54 | 5.30 | −4.39 | Eq. (22) |

−7.51 | 5.23 | −4.18 | MC | |

$E[M]$ | −0.1 | −0.99 | −0.78 | Eq. (27) |

−0.11 | −0.99 | −0.77 | MC | |

$E[(Fx\u2212E[Fx])2]$ | 74.47 | 78.19 | 81.14 | Eq. (23) |

76.9 | 77.71 | 80.72 | MC | |

$E[(Fy\u2212E[Fy])2]$ | 110.65 | 112.09 | 112.75 | Eq. (24) |

112.56 | 113.35 | 112.57 | MC | |

$E[(M\u2212E[M])2]$ | 0.46 | 0.5 | 0.55 | Eq. (28) |

0.46 | 0.51 | 0.54 | MC |

m = 4 | m = 5 | m = 6 | ||
---|---|---|---|---|

$E[Fx]$ | 32.05 | 39.46 | 41.17 | Eq. (21) |

32.09 | 39.55 | 41.11 | MC | |

$E[Fy]$ | −7.54 | 5.30 | −4.39 | Eq. (22) |

−7.51 | 5.23 | −4.18 | MC | |

$E[M]$ | −0.1 | −0.99 | −0.78 | Eq. (27) |

−0.11 | −0.99 | −0.77 | MC | |

$E[(Fx\u2212E[Fx])2]$ | 74.47 | 78.19 | 81.14 | Eq. (23) |

76.9 | 77.71 | 80.72 | MC | |

$E[(Fy\u2212E[Fy])2]$ | 110.65 | 112.09 | 112.75 | Eq. (24) |

112.56 | 113.35 | 112.57 | MC | |

$E[(M\u2212E[M])2]$ | 0.46 | 0.5 | 0.55 | Eq. (28) |

0.46 | 0.51 | 0.54 | MC |

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.

### Reconfiguration Scheme.

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\xaf)/max(T\xaf)$. Maximizing the TF provides an improved tension balance among the cables and it has been shown [36] that it is equivalent to minimizing $JTF=\u2211i=1mT\xafi$. So, in Eq. (36), we have $J=JTF=\u2211i=1mzi$. Second and third objective functions are $J=Var(Fx)$ (by choosing $\gamma 1=1,\u2009\gamma 2=0$, and $\gamma 3=0$) and $J=Var(Fy)$ (by choosing $\gamma 1=0,\u2009\gamma 2=1$, and $\gamma 3=0$), respectively. For all optimization problems, we set $Wd=[45\u2009N,\u200960\u2009N,\u20090\u2009N\u22c5m]T,\u2009\lambda \theta i=100$, and $\sigma Ti=2\u2009N\u2009\u2200i=1,\u2026,m$. The lower and upper bounds on the mean tension of the cables are assumed to be $T\xafLi=0\u2009N$ and $T\xafUi=100\u2009N\u2009\u2200i=1,\u2026,m$. Assuming a minimum tension of $0\u2009N$ allows the algorithm to remove a cable if the optimal configuration can be achieved by fewer number of cables. Moreover, $\theta \xafi$'s can be chosen anywhere on the circle, i.e., $\theta \xafi\u2208[0,\u20092\pi ]$. 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, $\theta \xafi$'s are restricted by certain bounds.

For the optimal $T\xafi$'s and $\theta \xafi$'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 $\theta 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 $\theta \xafi$'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\u2009N2$ in Fig. 4(a) to $9.95\u2009N2$ 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\u2009N2$ (Fig. 4(a)) to $7.6\u2009N2$ (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\u2009m$ 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\u2009N,\u200960\u2009N,120\u2009N\u22c5m]T,\u2009\lambda \theta i=100$, and $\sigma Ti=2\u2009N\u2009\u2200i=1,\u2026,m$. Additionally, the mean tension of the cables are restricted to $T\xafLi=0\u2009N$ and $T\xafUi=40\u2009N\u2009\u2200i=1,\u2026,m$. The domain from which $\theta \xafi$'s can be chosen are shown by dashed arcs whose lengths are $\pi \u2009rad$. As shown in Fig. 5, the variance of the output moment can be reduced from $56.04\u2009N\u22c5m2$ (when TF is maximized) to $39.02\u2009N\u22c5m2$ (when $Var(M)$ is minimized) in a system with *m* = 3 cables. Similarly, when *m* = 4, $Var(M)$ is reduced from $61.91\u2009(N\u22c5m)2$ to $37.09\u2009(N\u22c5m)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.

## Experimental Study

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.

### Agents Force Capacity and Uncertainty Characterization.

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 $\sigma Ti,\u2009\lambda \theta i,\u2009T\xafiL$,and $T\xafiU$ for $i=1,\u2026,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\xafiL$ and $T\xafiU$. 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 $\sigma 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\u22121={u\u2208\mathbb{R}n:||u||=1}$, where in planar system, *n* = 2, and $S1$ is a unit circle. Let $Uj={u1j,\u2026,uNj}$ be the set of *N* samples of **u** corresponding to the *j*th motion command. It can be seen that the data sets *U ^{j}* ($j=1,\u2026,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 $\lambda \theta 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\xafnp=[0,\u20091]T$. Using Rodrigues' formula [42], we calculate a rotation matrix

*R*that rotates the sample mean vector $u\xafj$ (corresponding to data set

_{j}*U*) to $u\xafnp$, i.e., $u\xafnp=Rju\xafj$. Then, the data set $U$ is constructed as follows:

^{j}The samples in $U$ have common mean vector $u\xafnp$ and dispersion parameter *λ _{u}*. Now, we can use Eq. (11) to estimate $\lambda \theta i$.

Table 2 summarizes the estimated parameters corresponding to three iRobots used in this study.

iRobot 1 | iRobot 2 | iRobot 3 | |
---|---|---|---|

$T\u0302\xafiL\u2009(N)$ | 3 | 3 | 4 |

$T\u0302\xafiU\u2009(N)$ | 5 | 6 | 6 |

$\sigma \u0302Ti\u2009(N)$ | 0.95 | 1.13 | 1.13 |

$\lambda \u0302\theta i$ | 1.208 × 10^{3} | 338.7 | 778.7 |

iRobot 1 | iRobot 2 | iRobot 3 | |
---|---|---|---|

$T\u0302\xafiL\u2009(N)$ | 3 | 3 | 4 |

$T\u0302\xafiU\u2009(N)$ | 5 | 6 | 6 |

$\sigma \u0302Ti\u2009(N)$ | 0.95 | 1.13 | 1.13 |

$\lambda \u0302\theta i$ | 1.208 × 10^{3} | 338.7 | 778.7 |

### Reconfiguration Scheme: Minimizing the Uncertainty Level.

We implement different optimization scenarios simulated in Sec. 5 on the cooperative system shown in Fig. 6. The desired task is determined to provide an output wrench $Wd=[3\u2009N,\u20095\u2009N,0\u2009N\u22c5m]T$. However, as the system output is uncertain, we can only achieve the desired wrench in an average sense. We set the acceptable distance from the mean wrench vector to be $1%$ of the desired wrench. Figure 9 shows both the simulation and real-system experimental results corresponding to three optimization scenarios: TF maximization (Figs. 9(a) and 9(b)), minimization of $Var(Fx)$ (Figs. 9(c) and 9(d)), and minimization of the $Var(Fy)$ (Figs. 9(e) and 9(f)).

Referring to Fig. 9, the results obtained through experimental studies are in a good agreement with simulation results. The first row shows the result of our attempt to maximize the tension factor and simultaneously maintaining a mean output wrench close to the desired wrench. The simulation results predict $Var(Fx)$ and $Var(Fy)$ to be $2.36\u2009N2$ and $1.16\u2009N2$, respectively. The sample estimate of these variances, obtained from the experimental measurements, are $2.68\u2009N2$ and $1.43\u2009N2$, respectively. The KS joint density plot for experimental data is also closely resembling the one obtained through the numerical simulation. In the second row of Fig. 9, we are aiming to rotate the uncertainty ellipse such that its minor axis is aligned with the *x*-direction, i.e., the variance of $Fx$ is minimized. Figure 9(d) verifies that reconfiguration of the iRobots according to the optimal solution can successfully minimize the variance of $Fx$. It can be seen that $Var(Fx)$ reduces from $2.68\u2009N2$, when we maximize TF (Fig. 9(b)) to $0.35\u2009N2$, when the objective is minimizing $Var(Fx)$ (Fig. 9(d)). Similarly, when the objective is to minimize $Var(Fy)$, the simulation predicts the variance to be $0.49\u2009N2$ and $Var(Fy)$ is $0.43\u2009N2$ in the actual system that is a significant reduction compared with $Var(Fy)=1.43\u2009N2$, when TF is maximized. Different factors are responsible for the small difference between simulation and experimental results. In our simulation, we ignored the dependency of the $\sigma Ti$ and $\lambda \theta i$ on $T\xafi$. Moreover, depending on the position of the iRobots and the traction provided by each robot, cables became slack for a very short time during the experiment that resulted in unexpected drop in the output force (e.g., small cluster of points in Fig. 9(f)). There is also an unavoidable error in adjusting the orientation of the cables (configuring the iRobots) that is not taken into account in our formulation. In addition, we have not considered the sensing uncertainty associated with the force transducer. However, neglecting these factors enables the analytical treatment to the problem as developed in Secs. 3 and 4, and as presented in Fig. 9, the proposed model can still provide high fidelity, hence, justifying these assumptions.

## Discussion

The quantification of wrench uncertainty in cable-robot systems was considered in this paper. By taking advantage of the special lay out of the parallel (cable) architectures, we were able to formulate the resultant wrench in an uncertain system through sum of the random vectors. The uncertainty induced by each agent was considered due to the variation in the tension and orientation of the cables that were modeled by Gaussian and von Mises distributions, respectively. The closed-form expressions of the means and variances of the resultant force and moment were provided. This enabled us to efficiently and accurately quantify the wrench uncertainty for a given state of the system. This, in turn, facilitates off-line or real-time optimization procedures that adaptively optimize the system state (configuration and cable tensions) in order to improve precision indices (functions of force and moment variances), i.e., reducing the uncertainty level. We provided examples of the static cable-based systems to examine the effectiveness of the proposed approach. Further, in an experimental study, we implemented the developed method on a real cooperative system where analysis of the system output verified the applicability and acceptable accuracy of the proposed approach in reducing the uncertainty level.

One of the key assumptions in this work was the statistical independency of the random variables. While such assumptions are legitimate for loosely interconnected systems, further experimental investigations are needed for systems with more mechanical interconnections. Moreover, our development was based on the assumption that the cooperative task is performed in a quasi-static condition. More complex uncertainty analysis is required for (highly) dynamic manipulations. Furthermore, spatial systems can be partially analyzed by the approach proposed in this paper; however, our formulation for planar systems can be extended to spatial systems by using spherical statistics [7]. The 3D extension of this study will apply to a variety of the cooperative systems such as aerial towing with flying robots (e.g., quadcopters).

## Funding Data

National Science Foundation (Award Nos. IIS-1319084 and CNS-1314484).