Let the complete joint input vector from the controller (i.e. desired reference values to the plant model) be expressed as q = l s l c and the realized joint vector (actual values at time t) be given by q ′ = l s ′ l c ′ . The moving slider bases are modeled as single degree-of-freedom (DOF) second order linear systems subject to maximum and minimum velocity constraints to mimic the motion of a screw-type linear actuator. The bases are modeled with inertia and friction and a simple proportional controller as shown in Eq. 1. For every ith actuator in the vector l _s (dropping the subscript i to indicate each individual slider joint) we model the dynamic equation as, m l ̈ s ′ + d l ̇ s ′ = g p , l s l s − l s ′ s.t.  l ̇ s , min ′ < l ̇ s ′ < l ̇ s , max ′ (1) where g _p,ls is the proportional gain and the maximum velocity that a linear screw actuator can travel in is limited to a small value ( | l ̇ s , min ′ | = | l ̇ s , max ′ | ) . The motor cable spooling l c ′ is modeled as a first order lag system, i.e. T s l ̇ c ′ + l c ′ = g p , l c l c − l c ′ (2) where T _s is the settling time and g p , l c is the proportional gain. In addition to the specified models and constraints, there are also limits placed on the position values of the actuator joints in order to set hard constraints defined by geometric limits.

The cable failure is modeled as a linear drop in the cable stiffness of the failed joint, i.e., k _i → 0 over a small but finite time interval Δt _f . Thus if there is no change in winch angle as the cable fails, the tension in the corresponding cable falls. The full dynamic model of the rCDPR can be derived via the Lagrangian approach as M x ̈ e + D x ̇ e x ̇ e = J w x e , l s ′ K q l p x e , l s ′ − l c ′ (3) where end-effector pose is given by x e = x y ϕ T and J w ( x e , l s ′ ) is the pulling map or wrench Jacobian. The prismatic length, l _p , comes from inverse kinematics. Both are dependent on the current pose, x _e , and the realized geometric attachment points, l s ′ , while the joint state stiffness matrix, K _q , is a diagonal matrix with the elements (k ₁, … , k ₄). M is the mass/inertia matrix whose coefficients are informed by the mass and shape of the end effector. This formulation ignores the effects of Coriolis forces. The damping coefficients of the damping matrix D are minor and expected to be a result of surface friction alone. An effect of low cable stiffness and damping constants is that oscillations induced due to system dynamics in operation or during cable failure do not dissipate easily. A feedforward controller for such a system would be intractable. The complete dynamic system is shown in Figure 4.

The motion models employed in the Extended Kalman Filters for the IMM (Section. 3.2) and the LQR (Section 5.4) use Eq. 4 and 5 as the motion model. Three motion models are employed to describe the seven working modes 4.1. For each working mode, the effect of the corresponding cable failure is incorporated in the dynamic model through the loss of stiffness in the respective cables. However, the model does not consider the true stiffness constant value from the plant model but instead employs static parameters corresponding to the different failure modes. In addition, it does not consider the model for actuator dynamics, but instead accepts a measurement of the current joint state, q ̃ ′ = q ′ + v ′ , where v′ is modeled as a zero mean Gaussian delta correlated noise. On the other hand, the LQR assumes that the joint input to the motion model is the desired joint state. This is further discussed in Section 5.4.

The kinematic feedforward controller does not use the dynamic models and instead relies on the computation of the forward and inverse kinetostatic solutions for the redundant system. The cable stiffness model used to compute the desired cable tensions is approximated to consider the true cable length instead of the free cable length ( k = k 0 l p ) . This was done to remove the direct dependence of the cable stiffness parameter on l _c which is a control variable.

These assumptions add unmodeled uncertainties to the control and estimation problem. The complexity introduced by this bears resemblance to similar challenges associated with real-world experimental implementations. In addition, only a single measurement source (at the end-effector) is considered. It must be emphasized that the simulated measurement noise is unbiased and Gaussian. The deteriorating effect of measurement bias (such as IMU drift) on most estimation approaches is well known and well explained and are not considered in this work. Sensors such as ultrasonic sensors, GPS devices or visual systems can help eliminate bias. The selected measurement noise characteristics in this work was not based on simulations of a real measurement device.

Cable failure in a CDPR reduces the degrees of freedom (DOFs) or the quality of the DOF (i.e. manipulability or dexterity) available at the end-effector. We control only 2 DOFs at the end-effector for the four-cable planar CDPR model in this paper. The total working modes considered can be described by 3 motion models (i) an over-constrained CDPR with 4 cables, (ii) a fully constrained CDPR with 3 cables and (iii) and under-constrained CDPR with two cables (Figure 5).

Therefore, there are 7 total working modes (mode 1) no cables have failed, (modes 2–5) only cable A, B, C, or D fail, (mode 6) cables A and B have failed, (mode 7) cables C and D have failed. If any other failure combination occurs, significant areas of the workspace become unreachable and are considered to be complete failure modes (i.e. cannot be recovered from). If the sliders were not constrained to geometric limits and free to move like mobile robots or drones, then these additional failure modes can also be considered.

An interactive multiple model filter (IMM) is a dynamic estimator (Blom and Bar-Shalom (1988)) which can be used when model changes appear suddenly or gradually over time, providing means for failure detection. The process flow of an IMM is shown in Figure 6.

In general, the bank of filters in an IMM can be realized by a variety of estimation approaches, i.e., Kalman Filters (KF), Extended Kalman Filters (EKF), Unscented Kalman Filters (UKF), Particle Filters (PF), et cetera. In this study, we use an EKF formulation in which the error dynamics are approximated by a first order Taylor series expansion about the estimated state, x ̂ , i.e. F = ∂ f x , u ∂ x , H = ∂ h x ∂ x (6) Here, F _{6 × 6} is determined numerically by a complex step derivative approximation and H _{3 × 6} is simlpy given by I 3 × 3 0 3 × 3 , since we directly measure the end-effector’s position. In this study, the EKF is applied in the discrete time domain; hence, the Euler-Maruyama discretization of the motion model yields. x k = x k − 1 + f x k − 1 , q k − 1 ′ Δ t ⏟ ϕ x k − 1 , q k − 1 ′ + ω k Δ t (7) y ̃ k = h x k + v k (8) where y ̃ ∈ R m is the measurement vector, ω k ∼ N ( 0 , Q k ) denotes the process noise with covariance Q k ∈ R n × n , and v k ∼ N ( 0 , R k ) is the measurement noise with covariance R k ∈ R m × m . The sampling time is denoted by Δt, and k is the current time step.

Let M be the number of all working modes in a CDPR at the kth time step. Then, the input to the filter is the current measurement and mixed state estimate x ̂ + 0 . The algorithm broadly consists of four stages (Gill (2019)).

1. Interaction and Mixing: The weights and the estimates from the last cycle are mixed as per their associated Markov transition probabilities given by the Markov transition matrix, A (also sometimes called the Markov switching matrix). The predicted probability for the filter to end up in mode j in each cycle, given that it was in the mode i during the previous cycle is given by

• Propagation: x ̂ k − j = x ̂ k + 0 j + f x ̂ k + 0 j , q k ′ Δ t P k − j = Φ k j P k + 0 j Φ k j T + Q k Φ j = I + F k j Δ t (11)

• And update: K k j = P k − j H k j T H k j P k − j H k j T + R k − 1 x ̂ k + j = x ̂ k − j + K k y ̃ k − h x ̂ k − j P k + j = I − K k H k P k − j (12)

Here, ‘+’ denotes the a posteriori estimate whereas ‘−’ denotes the a priori quantity before the update. The likelihood of a measurement is then given by p y ̃ k | x ̂ k − j = 1 det 2 π E k − j 1 / 2 ⋅ exp − 1 2 e k − j T E k − j − 1 e k − j with E k − j = H k j P k − j H k j T + R k e k − j = y ̃ k − y ̂ k − j (13) and where e k − j is the estimation error.

3. Model Probability Update: Now, the model probabilities are updated via the likelihood with subsequent normalization, i.e.

The weight vector, w _k , displays the importance associated with each model in the bank. For example, if there are no cable failures, the IMM will determine the first model (mode 1) to have the largest weight. This vector also informs the balance of joint inputs to the controller as described in the next section.

We first discretize and parameterize the reference trajectory into a set of constant arc-length waypoints W s . W s = p i = x i , y i i ∈ 1 → N i (16)

where N _i is the total number of waypoints. The constant arc-length waypoints are numerically determined from the path equation. This re-parameterization ensures that regardless of the chosen trajectory equation the distance between any two waypoints is never too large. Waypoints function as an ordered set of markers without time parameterization that can be dropped from the waypoint map as the robot traverses within a threshold of it. The robot always moves toward the next waypoint in the map regardless of cable failure or any other unpredictable disturbance at the end-effector that may temporarily navigate the robot away from the desired trajectory.

Next, we use a constant-velocity carrot following controller to traverse straight line intersections between any given pair of waypoints.

The shortest distance from the current position of the robot end-effector origin, shown in Figure 7 by frame {O}, intersects with the line segment at the point p _n and forms a normal vector e whose magnitude forms the cross-track error. If the projection of the robot frame origin {O} lies outside of the line segment, the shortest distance is simply the distance to the closest end-point. Therefore, the cross-track error vector is determined by first projecting the robot origin on the line segment p through a dot product and capping the magnitude of the projection ‖ p ‖ ∈ [ 0 , l W , i ] where l W , i is the length of the current line segment. If o W is the vector between frame {O} and W i , then the cross-track error vector is always given by e = o W + ‖ p ‖ p ̂ . The reference velocity is a constant and is given by v _ref . Thus the reference velocity vector which lies along the current line segment is v r e f = v r e f p ̂ . The desired direction of the velocity vector is given by v _d = e + v _ref . The final desired velocity vector has its velocity always capped at v _ref and is therefore given by v = v r e f v ̂ d .

As the cross-track error approaches zero, the desired velocity vector and the reference velocity vector are more closely aligned. The joint states required to achieve this velocity are computed through redundancy resolution-based inverse kinematics. When the end-effector reaches within a certain threshold of the next waypoint, the waypoint map W s updates by dropping the first value from the list. This strategy allows the robot to always move towards the second waypoint in W s , while the first waypoint serves as a marker for the connecting line segment to be tracked. If the robot snaps or moves off track, it still moves towards the second waypoint ensuring that the trajectory we want to track is traversed approximately from where it deviates. The choice of arc length chosen above determines the maximum size of the missed line segment.

For a given end-effector pose, there are infinitely many solutions in the joint space of kinematically redundant CDPR. Considering that this system has two kinds of redundancies, actuation and kinematic redundancy, the solutions for determining the optimal joint states can be quite complex. In this study, a feedforward kinematic controller was applied and as a result, the redundancy resolution schemes and joint determination use kinematic models.

Given the desired joint state q _0,n of the nth working mode, the desired steady state end-effector pose x _0,n can be determined through forward kinematics which is implemented by minimizing the potential energy of the system. x 0 , n = arg min x 0 , n U = 1 2 l p − l c 0 , n T K q l p − l c 0 , n (23) s.t. l p − l c 0 , n > 0 (24) where l _p is true cable length calculated for each value of x _0,n . The cable tensions given by K q ( l p − l c 0 , n ) must always be positive. If this condition is not met, it indicates that the initial joint state provided was not feasible which is generally a result of one or more cables being slack. Fortunately, in the feedforward control setting considered in this work, the desired winch position l c 0 , n employed in the calculation of the forward kinematics is always less than the prismatic cable length l _p . This calculation is not influenced by the true input l _c that is sent to the plant. Thus, although cable slackness can exist in the plant model, it is not encountered during the calculation of the forward kinematics.

The task recovery algorithm has a parallel bank of trajectory tracking controllers, each corresponding to a working mode. Each working mode, n, runs an independent kinematic controller that propagates the desired joint inputs in time such that it is the only working mode. The controllers use simple kinematic models for the different redundant CDPR motion models.

The control loop accepts the current joint inputs, q to the plant and determines the steady-state robot pose through forward kinematics for the nth model (Figure 8). The way point tracking algorithm accepts the current state estimate, x ̂ to update the waypoint and reference trajectory information if necessary. Through the carrot following algorithm, the desired task space state is determined and the redundancy resolution scheme is applied to determine the desired joint state, q _0,n . This is propagated through forward kinematics again to determine the steady-state end-effector pose, x _0,n .

Although the full system is highly non-linear, local stabilization (vibration dampening) is performed utilizing a Linear Quadratic Regulator. The LQR assumes the system model in Eq. 5 and ignores the actuator dynamics. This does not cause unwelcome destabilizing effects as the update frequency of the LQR is low compared to the dynamic joint controllers. The system model is linearized around the desired steady-state operating point (x _0,n , q _0,n ) that results from a combination of forward and inverse kinematics in the feedforward controller described in Section 5.2. The local system dynamics then arise as δ x n , k + 1 = A n δ x n , k + B n δ q n , k A n = ∂ f n x , q ∂ x x 0 , n , q 0 , n , B n = ∂ f n x , q ∂ q x 0 , n , q 0 , n (25) where k is the time-step index within the LQR feedback control loop. Please note that the desired steady-state operation point (x _0,n , q _0,n ) is also time-varying, but at a much slower update rate than the LQR loop. The slow update rate of the feedforward controller that provides the desired steady-state operating point is useful here as the redundancy resolution algorithm is the most computationally expensive piece of the framework. Therefore, a zero-order hold assumption without loss of generality is implied in Eq. 25. The dynamic correction δ q _n,k for the kinematically optimal joint input q _n,k is now determined via a quadratic cost function, i.e. δ q n , k = arg min q n ∑ k = 1 ∞ δ x n , k T Q δ x n , k + δ q n , k T R δ q n , k (26) subject to the dynamics in Eq. 25 and where Q and R are positive semi-definite and positive definite weighting matrices, respectively. The tuning of these terms depends on the desired overall impact of the corresponding states on the cost function as well as limits on the control effort. Here, the structure for Q and R has been chosen as diagonal. Furthermore, the diagonal terms for the joint inputs ( R ) are weighted equally. For Q , the entries corresponding to the velocities in x and y outweigh the other values, as it is desired to minimize velocity offsets dominantly. The resulting controller is then given by q n , k = q 0 , n − K n δ x ̂ n , k = q 0 , n − K n x ̂ k − x 0 , n (27) where K n is the LQR gain resulting for the nth model while the estimated state x ̂ k from the IMM framework has been employed as the feedback variable for the controller. When tuned well, the LQR stabilization not only mitigates dynamic disturbances due to the elasticity of the cables during normal operation, but also dampens significant oscillations from cable breakage quickly. An additional beneficial effect from the LQR inclusion is observed, as even slackness after failure is more efficiently eliminated. Here, the kinematic solution provides a new equilibrium point, and the LQR inherently rushes the system to arrive at this new stable kinematic solution, thus significantly aiding in maintaining positive cable tensions. For the over-constrained (4-cable) and fully constrained (3-cable) rCDPR models, only θ _0,n is subject to LQR stabilization as this is enough to dampen perturbations in both x- and y-directions. For the under-constrained (2-cable) model, the redundant joints l _0,n are also utilized to dampen disturbance components orthogonal to the cable direction.

The bank of controllers independently determine the next steady-state input for each failure mode. The bank of LQRs generate optimal joint inputs that are mixed with the weighting matrix to the plant model. The final mixed input gives precedence to the dominant working mode identified by the IMM. Due to this mixing, the final mixed input from this stage will always be subject to some perturbation. However, this perturbation is considered a reasonable trade-off to accommodate a solution that is robust and reliable during cable failure conditions.

The mixed joint input is applied to the dynamic model (plant) and the measurements of the end-effector state inform the state estimation and weight vector through the IMM. The complete algorithmic framework is described in Figure 9.

The plant model described in Section 3 is perturbed with a process and measurement noise to mimic unmodeled dynamics. The cables are modeled as linear springs (where stiffness is a function of free cable length) with no mass and no sag, while all other model parameters are constants. Cable failure is modeled as a drop in cable stiffness to zero over a finite time. If the stiffness is zero, the cable can have no tension and thus has no impact on the platform dynamics. In addition, if the dynamic input to the system causes a situation where l _c > l _p , the tension in the dynamic equation is set to zero to mimic cable sag. The plant dynamics and estimation run at a frequency of 1,000 Hz while the closed loop LQR feedback is run at 50 Hz. The feedforward kinematic model, which is the most computationally-intensive system, is run at 2 Hz. The entire framework is modeled in MATLAB. The model parameters are given in Table 1.

Figure 10 demonstrates the application of the complete task recovery algorithm. The red trajectory shows the true path of the center of the end-effector while the green trajectory shows the desired path. The blue dots are the predefined waypoints. The failure of the first cable causes the end-effector to jump away from its current path, but the algorithm immediately responds by rearranging the remaining three winches and pulling the end-effector back. This breakaway from the trajectory and subsequent recovery with three cables is shown by the first kink in the red line part way through the trajectory in Figure 10B. The failure of the second cable and subsequent task recovery is seen further down the trajectory in Figure 10C.

The IMM correctly identifies the failure, mixing and switching the inputs as needed to automatically recover the trajectory tracking. In Figure 11, the weights vector correctly assigns the largest values to the current working mode. At the forty second mark, we can see that the weight vector has correctly identified the cable A failure, followed by the correct identification of the cable C failure at the 10 s mark. The estimation error plots are in Figure 12.

The robot is maintained at constant velocity throughout its trajectory (Figure 13). The desired velocity during trajectory tracking is 0.005 m/s or 5 mm/s which is a reasonable speed for the dimensions of the workspace. During recovery, the end-effector velocity is not a constraint, so the controller, mainly through the efforts of the LQR, attempts to move the end-effector as fast as dynamically feasible.

The trajectory error is minimized by the combined efforts of feedforward controller and the LQR. The banded regions in Figure 14 represent sections where the robot is in recovery mode, attempting to minimize the trajectory error while compensating for the changes in slider and winch position.

The robot does not enter singularities in the feedforward controller. Figure 15 shows the plots of a single slider position (slider C) for the second motion model. In the example trajectory shown in Figure 10 this model is weighted highest between the 40 s and 80 s timestamps. Thus the mixed output from the LQR, l begins to weight heavily in the favour of the control input l _0,2 from model 2, while the true slider position l′ begins to catch up. When cable C breaks at the 80th second, this model begins to fall out of favor and the slider position deviates from this value.

Figures 16, 17 depict the change in slider position and position angle of winch respectively. The shaded regions represent the failed state of the respective cable. The values of the slider and winch positions in these regions are perturbed by the combined attempts of the optimization to continue modulating this actuator and the weighted input mixing but has no effect on the end-effector as the corresponding cable is broken. This is seen clearly in the plots observing the true tension in the cables in Figure 18. The position values are not automatically settled to zero to account for cases where the cable failure may have been misidentified by the estimator.

A minimum value of pre-tension in the cables is maintained at 0.5 N as a good practice to eliminate the chances of encountering cable slackness. This tension is seen to be maintained throughout the run with the minimum required tensions in the remaining cables to propagate the end-effector along the trajectory. This demonstrates that the tension minimization scheme is working well. When a cable snaps, there are brief moments of cable slackness while the tensions re-distribute as the sliders reposition themselves and the winches wind within their designed dynamics. These dynamics do not affect the steady-state input from the feedforward controller and thus the LQR bears the burden of recovery.

The redundancy resolution optimization for the slider positions maximizes the manipulability and wrench quality. In Figure 19 we see these properties hold high values in the first two-thirds of the plot. These are the regions where the optimization-determined slider positions are actively applied. This optimization is no longer relevant for the two cable cases, where the slider position instead functions to replace the lost degree of freedom.

For completeness, this section compares the results of the previous section with a feedforward controller applied directly on the system without the LQR. The update rate of the feedforward controller is 50 Hz and all other parameters are kept the same. The trajectory results are presented in Figure 20.

The trajectory error plot in Figure 21 indicates that the oscillations do eventually dampen and the kinematic controller is generally successful after this due to the slow velocities at the end-effector.

The cable breakage causes the end-effector to swing causing periodic slackness in the cable tensions (Figure 22).

The estimation depends solely on the end-effector measurements, so it is unaffected by the controller dynamics as seen in Figure 23.