A Lyapunov-Stable Adaptive Method to Approximate Sensorimotor Models for Sensor-Based Control

In this article, we present a new scheme that approximates unknown sensorimotor models of robots by using feedback signals only. The formulation of the uncalibrated sensor-based regulation problem is first formulated, then, we develop a computational method that distributes the model estimation problem amongst multiple adaptive units that specialize in a local sensorimotor map. Different from traditional estimation algorithms, the proposed method requires little data to train and constrain it (the number of required data points can be analytically determined) and has rigorous stability properties (the conditions to satisfy Lyapunov stability are derived). Numerical simulations and experimental results are presented to validate the proposed method.


INTRODUCTION
Robots are widely used in industry to perform a myriad of sensor-based applications ranging from visually servoed pick-and-place tasks to force-regulated workpiece assemblies (Nof, 1999). Their accurate operation is largely due to the fact that industrial robots rely on fixed settings that enable the exact characterization of the tasks' sensorimotor model. Although this full characterization requirement is fairly acceptable in industrial environments, it is too stringent for many service applications where the mechanical, perceptual and environment conditions are not exactly known or might suddenly change (Navarro-Alarcon et al., 2019), e.g., in domestic robotics (where environments are highly dynamic), field robotics (where variable morphologies are needed to navigate complex workspaces), autonomous systems (where robots must adapt and operate after malfunctions), to name a few cases.
In contrast to industrial robots, the human brain has a high degree of adaptability that allows it to continuously learn sensorimotor relations. The brain can seemingly coordinate the body (whose morphology persistently changes throughout life) under multiple circumstances: severe injuries, amputations, manipulating tools, using prosthetics, etc. It can also recalibrate corrupted or modified perceptual systems: a classical example is the manipulation experiment performed in Kohler (1962) with image inverting goggles that altered a subject's visual system. In infants, motor babbling is used for obtaining (partly from scratch and partly innate) a coarse sensorimotor model that is gradually refined with repetitions (Von Hofsten, 1982). Providing robots with similar incremental and life-long adaptation capabilities is precisely our goal in this paper.
From an automatic control point of view, a sensorimotor model is needed for coordinating input motions of a mechanism with output sensor signals (Huang and Lin, 1994), e.g., controlling the shape of a manipulated soft object based on vision  or controlling the balance of a walking machine based on a gyroscope (Yu et al., 2018). In the visual servoing literature, the model is typically represented by the so-called interaction matrix (Hutchinson et al., 1996;Cherubini et al., 2015), which is computed based on kinematic relations between the robot's configuration and the camera's image projections. In the general case, sensorimotor models depend on the physics involved in constructing the output sensory signal; If this information is uncertain (e.g., due to bending of robot links, repositioning of external sensors, deformation of objects), the robot may no longer properly coordinate actions with perception. Therefore, it is important to develop methods that can efficiently provide robots with the capability to adapt to unforeseen changes of the sensorimotor conditions.
Classical methods in robotics to compute this model (see Sigaud et al., 2011 for a review) can be roughly classified into structure-based and structure-free approaches . The former category represents "calibrationlike" techniques [e.g., off-line (Wei et al., 1986) or adaptive (Wang et al., 2008;Liu et al., 2013;Navarro-Alarcon et al., 2015)] that aim to identify the unknown model parameters. These approaches are easy to implement, however, they require exact knowledge of the analytical structure of the sensory signal (which might not be available or subject to large uncertainties). Also, since the resulting model is fixed to the mechanical/perceptual/environmental setup that was used for computing it, these methods are not robust to unforeseen changes.
For the latter (structure-free) category, we can further distinguish between two main types (Navarro-Alarcon et al., 2019): instantaneous and distributed estimation. The first type performs online numerical approximations of the unknown model (whose structure does not need to be known); Some common implementations include e.g., Broyden-like methods (Hosoda and Asada, 1994;Jagersand et al., 1997;Alambeigi et al., 2018) and iterative gradient descent rules (Navarro-Alarcon et al., 2015;Yip et al., 2017). These methods are robust to sudden configuration changes, yet, as the sensorimotor mappings are continuously updated, they do not preserve knowledge of previous estimations (i.e., it's model is only valid for the current local configuration). The second type distributes the estimation problem amongst multiple computing units; The most common implementation is based on (highly nonlinear) connectionists architectures (Li and Cheah, 2014;Lyu and Cheah, 2018;Hu et al., 2019). These approaches require very large amounts of training data to properly constrain the learning algorithm, which is impractical in many situations. Other distributed implementations (based on SOM-like sensorimotor "patches, " Kohonen, 2013) are reported e.g., in Zahra and Navarro-Alarcon (2019), Pierris andDahl (2017), andEscobar-Juarez et al. (2016), yet, the stability properties of its algorithms are not rigorously analyzed.
As a solution to these issues, in this paper we propose a new approach that approximates unknown sensorimotor models based on local data observations only. In contrast to previous state-of-the-art methods, our adaptive algorithm has the following original features: • It requires few data observations to train and constrain the algorithm (which allows to implement it in real-time). • The number of minimum data points to train it can be analytically obtained (which makes data collection more effective). • The stability of its update rule can be rigorously proved (which enables to deterministically predict its performance).
The proposed method is general enough to be used with different types of sensor signals and robot mechanisms. The rest of the manuscript is organized as follows: section 2 presents preliminaries, section 3 describes the proposed method, section 4 reports the conducted numerical study, and section 5 gives final conclusions.

Notation
Along this note we use very standard notation. Column vectors are denoted with bold small letters m and matrices with bold capital letters M. Time evolving variables are represented as m t , where the subscript * t denotes the discrete time instant. Gradients of functions b = β(m) : M → B are denoted as ∇β(m) = (∂β/∂m) ⊺ .

Configuration Dependant Feedback
Consider a fully-actuated robotic system whose instantaneous configuration vector (modeling e.g., end-effector positions in a manipulator, orientation in a robot head, etc.) is denoted by the vector x t ∈ R n . Such model can only be used to represent traditional rigid systems, thus, it excludes soft/continuum mechanisms (Falkenhahn et al., 2015) or robots driven by elastic actuators . Without loss of generality, we assume that its coordinates are all represented using the same unitless range 1 . To perform a task, the robot is equipped with a sensing system that continuously measure a physical quantity whose instantaneous values depend on x t . Some examples of these types of configuration-dependent feedback signals are: geometric features in an image (Tirindelli et al., 2020), forces applied onto a compliant surface (Navarro-Alarcon et al., 2014;Bouyarmane et al., 2019), proximity to an object (Cherubini and Chaumette, 2013), intensity of an audio source (Magassouba et al., 2016), attitude of a balancing body (Defoort and Murakami, 2009), shape of a manipulated object (Navarro-Alarcon and , temperature from a heat source (Saponaro et al., 2015), etc.
Let y t ∈ R m denote the vector of feedback features that quantify the task; Its coordinates might be constructed with raw measurements or be the result of some processing. We model the instantaneous relation between this sensor signal and the robot's configuration as (Chaumette and Hutchinson, 2006): Remark 1. Along this paper, we assume that the feedback feature functional f (x t ) is smooth (at least twice differentiable) and its Jacobian matrix has a full row/column rank (which guarantees the existence of its (pseudo-)inverse).

Uncalibrated Sensorimotor Control
In our formulation of the problem, it is assumed that the robotic system is controlled via a standard position/velocity interface (as in e.g., Whitney, 1969;Siciliano, 1990), a situation that closely models the majority of commercial robots. With position interfaces, the motor action u t ∈ R n represents the following displacement difference: Such kinematic control interface renders the typical stiff behavior present in industrial robots (for this model, external forces do not affect the robot's trajectories). The methods in this paper are formulated using position commands, however, these can be easily transformed into robot velocities v t ∈ R n by dividing u t by the servo controller's time step dt as follows u t / dt = v t . The expression that describes how the motor actions result in changes of feedback features is represented by the first-order difference model 2 : where the configuration-dependent matrix A(x t ) = ∂f /∂x t ∈ R m×n represents the traditional sensor Jacobian matrix of the system (also known as the interaction matrix in the visual servoing literature; Hutchinson et al., 1996). To simplify notation, throughout this paper we shall omit its dependency on x t and denote it as A t = A(x t ). The flow vector δ t = A t u t ∈ R m represents the sensor changes that result from the action u t . Figure 1 conceptually depicts these quantities. The sensorimotor control problem consists in computing the necessary motor actions for the robot to achieve a desired sensor configuration. Without loss of generality, in this note, such configuration is characterized as the regulation of the feature vector y t toward a constant target y * . The necessary motor action to reach the target can be computed by minimizing the following quadratic cost function: where λ > 0 is a gain and sat(·) a standard saturation function (defined as in e.g., Chang et al., 2018). The rationale behind the minimization of the cost (4) is to find an incremental motor command u t that forward-projects into the sensory space (via the interaction matrix A t ) as a vector pointing toward the target y * . By iteratively commanding these motions, the distance y t − y * is expected to be asymptotically minimized.
To obtain u t , let us first compute the extremum ∇J(u t ) = 0, which yields the normal equation FIGURE 1 | Representation of a configuration trajectory x t , its associated transformation matrices A t and motor actions u t , that produce the measurements y t and sensory changes δ t .
Solving (5) for u t , gives rise to the motor command that minimizes J: Nakamura, 1991), whose existence is guaranteed as A t has a full column/row rank (depending on whichever is larger n or m). Yet, note that for the case where m > n, the cost function J can only be locally minimized.
Note that the computation of (6) requires exact knowledge of A t . To analytically calculate this matrix, we need to fully calibrate the system, which is too restrictive for applications where the sensorimotor model is unavailable or might suddenly change. This situation may happen if the mechanical structure of the robot is altered (e.g., due to bendings or damage of links), or the configuration of the perceptual system is changed (e.g., due to relocating external sensors), or the geometry of a manipulated object changes (e.g., due to grasping forces deforming a soft body), to name a few cases. Without this information, the robot may not properly coordinate actions with perception. In the following section, we describe our proposed solution.

Discrete Configuration Space
Since the (generally non-linear) feature functional (1) is smooth, the Jacobian matrix A t = ∂f /∂x t is also expected to smoothly change along the robot's configuration space. This situation means that a local estimation A of the true matrix A t around a configuration point x i is also valid around the surrounding neighborhood (Sang and Tao, 2012). We exploit this simple yet powerful idea to develop a computational method that distributes the model estimation problem amongst various units that specialize in a local sensorimotor map.
It has been proved in the sensor-based control community (Cheah et al., 2003) that rough estimations of A t (combined with the rectifying action of feedback) are sufficient for guiding the robot with sensory signals. However, note that large deviations from such configuration point x i may result in model inaccuracies. Therefore, the local neighborhoods cannot be too large.
Consider a system with N computing units distributed around the robot's configuration space (see Figure 2). The location of these units can be defined with many approaches, e.g., with selforganization (Kohonen, 2001), random distributions, uniform distributions, etc. (Haykin, 2009). To each unit, we associate the following 3-tuple: The weight vector w l ∈ R n represents a configuration x t of the robot where w l = x t . The matrix A l t ∈ R n×m stands for a local approximation of A t (w l ) evaluated at the point w l . The purpose of the structure D l is to store sensor and motor observations d t = {x t , u t , δ t }, that are collected around the vicinity of w l through babbling-like motions (Saegusa et al., 2009). The structure D l is constructed as follows: for τ > 0 as the total number of observations, which once collected, they remain constant during the learning stage. Note that x i and x i+1 are typically not consecutive time instances. The total number τ of observations is assumed to satisfy τ > mn.

Initial Learning Stage
We propose an adaptive method to iteratively compute the local transformation matrix from data observations. To this end, consider the following quadratic cost function for the lth unit: for F(u k ) ∈ R m×mn as a regression-like matrix defined as FIGURE 2 | Representation of the lth computing unit and the neighboring data used to approximate the local sensorimotor model. The black and red dashed depict the Gaussian and its square approximation. and a vector of adaptive parameters a l t ∈ R nm constructed as: where the scalarâ lij t denotes the ith row jth column element of the matrix A l t . The scalar h lk represents a Gaussian neighborhood function centered at the lth unit and computed as: where σ > 0 (representing the standard deviation) is used to control the width of the neighborhood. By using h lk , the observations' contribution to the cost (9) proportionally decreases with the distance to w l . The dimension of the neighborhood is defined such that h ≈ 0 is never satisfied for any of its observations x k . In practice, it is common to approximate the Gaussian shape with a simple "square" region, which presents the highest approximation error around its corners (see e.g., Figure 2 where the sampling point d τ +1 is within its boundary).
To compute an accurate sensorimotor model, the data points in (8) should be as distinctive as possible (i.e., the motor observations u t should not be collinear). This requirement can be fairly achieved by covering the uncertain configuration with curved/random motions.
The following gradient descent rule is used for approximating the transformation matrix A t at the lth unit: where u j k and δ i k denote the jth and ith components of the vectors u k and δ k , respectively.

Remark 2.
There are other estimation methods in the literature that also make use of Gaussian functions, e.g., radial basis functions (RBF) (Li and Cheah, 2014) to name an instance. However, RBF (in its standard formulation) use configurationdependent Gaussians to modulate a set of weights (which provide non-linear approximation capabilities), whereas in our case, the Gaussians are used but within the weights' adaptation law to proportionally scale the contribution of the collected sensorymotor data (our method provides a linear approximation within the neighborhood). Our Gaussian weighted approach most closely resembles the one used in self organizing maps (SOM) (Kohonen, 2013) to combine surrounding data observations.

Lyapunov Stability
In this section, we analyse the stability properties of the proposed update rule by using discrete-time Lyapunov theory (Bof et al., 2018). To this end, let us first assume that the transformation matrix satisfies: for any configuration x j around the neighborhood defined by D l (this situation implies that A(·) is constant around the vicinity of w l ). Therefore, we can locally express around w l the sensor changes as: where a l = [a l11 , a l12 , . . . , a lmn ] ⊺ ∈ R mn denotes the vector of constant parameters, for a lij as the ith row jth column of the unknown matrix A(w l ). To simplify notation, we shall denote F k = F(u k ).
Proposition 1. For a number mn of linearly independent vectors u k , the adaptive update rule (13) asymptotically minimizes the magnitude of the parameter estimation error a l t − a l .
Proof: Consider the following quadratic (energy-like) function: Computing the forward difference of V l t yields: for a symmetric matrix ∈ R mn×mn defined as follows: with H = diag(h l1 I τ , . . . , h lτ I τ ) ∈ R mτ ×mτ as a positive-definite diagonal matrix, I τ ∈ R τ ×τ as an identity matrix and ∈ R mτ ×mn constructed with τ matrices F k as follows: To prove the asymptotic stability of (13), we must first prove the positive-definiteness of the dissipation-like matrix (van der Algorithm 1: Compute a suitable γ . 1: γ ← initial value < 1, µ ← small step 2: repeat 3: γ ← γ − µ 4: until C > 0 Schaft, 2000). To this end, note that since the "tall" observations' matrix is exactly known and H is diagonal and positive (hence full-rank), we can always find a gain γ > 0 to guarantee that the symmetric matrix is also positive-definite, and therefore, full-rank. Next, let us rearrange mn linearly independent row vectors from as follows: which shows that has a full column rank, hence, the matrix = γ ⊺ C > 0 is positive-definite. This condition implies that V l t+1 − V l t < 0 for any a l t − a l = 0. Asymptotic stability of the parameter's estimation error directly follows by invoking Lyapunov's direct method (Bof et al., 2018).

Remark 3.
There are two conditions that need to be satisfied to ensure the algorithm's stability. The first condition is related to the magnitude of the learning gain γ . Large gain values may lead to numerical instabilities, which is a common situation in discretetime adaptive systems. To find a "small enough" gain γ > 0, we can conduct the simple 1D search shown in Algorithm 1. An eigenvalue test on C can be used to verify (20). The second condition is related to the linear independence (i.e., the non-collinearity) of the motor actions u t . Such independent vectors are needed for providing a sufficient number of constraints to the estimation algorithm (this condition can be easily satisfied by performing random babblinglike motions).

Localized Adaptation
Once the cost function (9) has been minimized, the computed transformation matrix A t locally approximates the robot's sensorimotor model around the lth unit. Note that the stability of the total N units is analogous the analysis shown in the previous section; A global analysis is out of the scope of this work.
The associated local training data (8) must then be released from memory to allow for new relations to be learnt-if needed. However, for the case where changes in the sensorimotor conditions occur, the model may contain inaccuracies in some or all computing units, and thus, its transformation matrices cannot be used for controlling the robot's motion. To cope with this issue, we need to first quantitatively assess such errors. For that, the following weighted distortion metric is introduced: where B > 0 denotes a positive-definite diagonal weight matrix to homogenize different scales in the approximation error e t = A s u t − δ t ∈ R m . The scalar index s is found by solving the search problem: To enable adaptation of problematic units, we evaluate the magnitude of the metric U t , and if found to be larger than an arbitrary threshold U t > |ε|, new motion and sensor data must be collected around the sth computing unit to construct the revised structure D s by using a push approach: that updates the topmost observation and discards the oldest (bottom) data, so as to keep a constant number τ of data points. The transformation matrices are then computed with the new data.

Motion Controller
The update rule (13) computes an adaptive transformation matrix A l t for each of the N units in the system. To provide a smooth transition between different units, let us introduce the matrix L t ∈ R m×n which is updated as follows 3 : where η > 0 is a tuning gain. The above matrix represents a filtered version of A s t , where s denotes the index of the active unit, as defined in (23). With this approach, the transformation matrix smoothly changes between adjacent neighborhoods, while providing stable values in the vicinity of the active unit; It can be seen as a continuous interpolation between adjacent neighborhoods.
The motor command with adaptive model is implemented as follows: The stability of this kinematic control method can be analyzed with its resulting closed-loop first-order system (a practice also commonly adopted with visual servoing controllers; Chaumette and Hutchinson, 2006). To this end, we use a small displacement approach (motivated by the local target provided by the saturation function), where we introduce the increment vector i = − sat(y t − y * ) and define the local reference position y = y t + i ∈ R m . Let us consider the case when the N units have minimized the cost functions (9). Note that the asymptotic Proposition 2. For n ≥ m (i.e., more/equal motor actions than feedback features), the "stiff " kinematic control input (26) provides the local feedback error y t − y with asymptotic stability.
Proof: Substitution of the controller (26) into the difference model (3) yields the closed-loop system: Adding ±y to (27) and after some algebraic operation, we obtain: which for a gain satisfying 0 < λ < 1, it implies local asymptotic stability of the small displacement error (y t − y) (Kuo, 1992).
Remark 4. Note that the above stability analysis assumes that robot's trajectories are not perturbed by external forces and that the estimated interaction matrix locally satisfies A t L # t A t ≈ A t around the active neighborhood.

CASE OF STUDY
In this section, we validate the performance of the proposed method with numerical simulations and experiments. A visionbased manipulation task with a deformable cable is used as our case of study (Bretl and McCarthy, 2014): It consists in the robot actively deforming the object into a desired shape by using visual feedback of the cable's contour (see e.g., Zhu et al., 2018). Soft object manipulation tasks are challengingand relevant to the fundamental problem addressed here-since the sensorimotor models of deformable objects are typically unknown or subject to large uncertainties (Sanchez et al., 2018). Therefore, the transformation matrix relating the shape feature functional and the robot motions is difficult to compute. The proposed algorithm will be used to adaptively approximate the unknown model. Figure 3 conceptually depicts the setup of this sensorimotor control problem.

Simulation Setup
For this study, we consider a planar robot arm that rigidly grasps one end of an elastic cable, whose other end is static; We assume that the total motion of this composed cablerobot system remains on the plane. A monocular vision sensor observes the manipulated cable and measures its 2D contour in real-time. The dynamic behavior of the elastic cable is simulated as in Wakamatsu and Hirai (2004) by using the minimum energy principle (Hamill, 2014), whose solution is computed using the CasADi framework (Andersson et al., 2019). The cable is assumed to have negligible plastic behavior. All numerical simulation algorithms are implemented in MATLAB.
FIGURE 3 | Representation of the cable manipulation case of study, where a vision sensor continuously measures the cable's feedback shape y t , which must be actively deformed toward y * .

FIGURE 4 | Various configurations of the visually measured cable profile (black solid line) and its approximation with Fourier series (red dashed line).
The cable simulation code is publicly available at https://github. com/Jihong-Zhu/cableModelling2D.
Let the long vector s t ∈ R 2α represents the 2D profile of the cable, which is simulated using a resolution of α = 100 data points. To perform the task, we must compute a vector of feedback features y t that characterizes the object's configuration. For that, we use the approach described in Digumarti et al. (2019) and Navarro-Alarcon and Liu (2018) that approximates s t with truncated Fourier series (in our case, we used four harmonics), and then constructs y t with the respective Fourier coefficients (Collewet and Chaumette, 2000). The use of these coefficients as feedback signals enable us to obtain a compact representation of the object's configuration, however, it complicates the analytical derivation of the matrix A t .

Approximation of the Matrix A t
To construct the data structure (8), we collect τ = 40 data observations d t at random locations around the manipulation workspace. Next, we define local neighborhoods centered at the configuration points w 1 = [0.3, 0.5], w 2 = [0.5, 0.5], w 3 = [0.5, 0.3], and w 4 = [0.5, 0.5]. These neighborhoods are defined with a standard deviation of σ = 1.3. With the collected observations, l = 1, . . . , 4 matrices A l t are computed using the update rule (14). Figure 4 depicts the measured shape (black solid line) of the cable at the four points w l and the shape that is approximated (red FIGURE 5 | Profile of the function G that is computed along the circular trajectory passing through the points in Figure 4; The "switch" label indicates the instant when A l t switches to different one. dashed line) with the feedback feature vector y t (i.e the Fourier coefficients). It shows that four harmonics provide sufficient accuracy for representing the object's configuration. To evaluate the accuracy of the computed discrete configuration space and its associated matrices A l t , we conduct the following test: The robot is commanded to move the cable along a circular trajectory that passes through the four points w l . The following energy function is computed throughout this trajectory: which quantifies the accuracy of the local differential mapping (3). The index l switches (based on the solution of 23) as the robot enters a different neighborhood. Figure 5 depicts the profile of the function G along the trajectory. We can see that this error function increases as the robot approaches the neighborhood's boundary. The "switch" label indicates the time instant when A l t switches to different (more accurate) matrix, an action that decreases the magnitude of G. This result confirms that the proposed adaptive algorithm provides local directional information on how the motor actions transform into sensor changes.

Sensor-Guided Motion
In this section, we make use of the approximated sensorimotor model to guide the motion of a robotic system based on feedback features. To this end, various cable shapes are defined as target configurations y * (to provide physically feasible targets, these shapes are collected from previous sensor observations). The target configurations are then given to the motion controller (26) to automatically perform the task. The controller implemented with saturation bounds of | sat(·)| ≤ 2 and a feedback gain λ = 0.1. Figure 6 depicts the progression of the cable shapes obtained during these numerical simulations. The initial y 0 and the intermediate configurations are represented with solid black curves, whereas the final shape y * is represented with red dashed curves. To assess the accuracy of the controller, the following cost function is computed throughout the shaping motions: For these four shaping actions, Figure 7 depicts the time evolution of the function E. This figure clearly shows that the feedback error is asymptotically minimized. Now, consider the setup depicted in Figure 8, which has two 3-DOF robots jointly manipulating the deformable cable. For this more complex scenario, the total configuration vector x t must be constructed with the 3-DOF pose (position and orientation) vectors of both robot manipulators as x t = Training of the sensorimotor model is done similarly as with the single-robot case described above; The same feedback gains and controller parameters are also used in this test. Figure 9 depicts the initial shape y 0 and intermediate configurations (black solid curves), as well as the respective final shape y * (red dashed curve) of the cable. Note that as more input DOF can be controlled by the robotic system, the object can be actively deformed into more complex configurations (cf. the achieved S-shape curve with the profiles in Figure 6). The result demonstrates that the approximated sensorimotor model provides sufficient directional information to the controller to properly "steer" the feature vector y t toward the target y * .
FIGURE 8 | Representation of a two-robot setup where both systems must jointly shape the cable into a desired form.  We now compare the performance of our method (using the same manipulation task shown in Figures 8, 9) with two state-of-the-art approaches commonly used for guiding robots with unknown sensorimotor models. To this end, we consider the classical Broyden update rule (Broyden, 1965) and the recursive least-squares (RLS) (Hosoda and Asada, 1994). These two methods are used for estimating the matrix A that is needed to compute the control input (6). To compare their performance,   the cost function E is evaluated throughout their respective trajectories; The same feedback gain λ = 0.1 is used for these three methods. Figure 10 depicts the time evolution of E computed with the three methods. This result demonstrates that the performance of our method is comparable to the other two classical approaches.

Experiments
To validate the proposed theory, we developed an experimental platform composed of a three degrees-of-freedom serial robotic manipulator (DOBOT Magician), a Linux-based motion control system (Ubuntu 16.04), and a USB Webcam (Logitech C270); Image processing is performed by using the OpenCV libraries (Bradski, 2000). A sampling time of dt ≈ 0.04 s is used in our Linux-based control system. In this setup, the robot rigidly grasps an elastic piece of pneumatic air tubing, whose other end is attached to the ground. The 3-DOF mechanism has a double parallelogram structure that enables to control the gripper's x-y-z position while keeping a constant orientation. For this experimental study, we only control 2-DOF of the robot such it manipulates the tubing with plane motions. Figure 11 depicts the setup.
We conduct similar vision-guided experiments with the platform as the ones described in the previous section. For these tasks, the elastic tubing must be automatically positioned into a desired contour. The configuration dependant feedback for this task is computed with the observed contour of the object by using two harmonic terms (Navarro-Alarcon and ). The sensorimotor model is similarly approximated around four configuration points (as in Figure 4), by performing random motions and collecting sensor data. Figure 12 depicts snapshots of the conducted experiments, where we can see the initial and final configurations of the system. The red curves represent the (static) target configuration y * . For these two targets, Figure 13 depicts the respective time evolution profiles of the energy function E, where we can clearly see that the feedback error is asymptotically minimized. The control inputs u t used during the experiments are depicted in Figures 14, 15. These motion commands are computed from raw vision measurements and a saturation threshold of ±1 is applied to its values. This results demonstrate that the approximated model can be used to locally guide motions of the robot with sensor feedback.

CONCLUSION
In this paper, we describe a method to estimate sensorimotor relations of robotic systems. For that, we present a novel adaptive rule that computes local sensorimotor relations in real-time; The stability of this algorithm is rigorously analyzed and its convergence conditions are derived. A motion controller to coordinate sensor measurements and robot motions is proposed. Simulation and experimental results with a cable manipulation case of study are reported to validate the theory.
The main idea behind the proposed method is to divide the robot's configuration workspace into discrete nodes, and then, locally approximate at each node the mappings between robot motions and sensor changes. This approach resembles the estimation of piecewise linear systems, except that in our case, the computed model represents a differential Jacobian-like relation.
The key to guarantee the stability of the algorithm lies in collecting sufficient linear independent motor actions (such condition can be achieved by performing random babbling motions).
The main limitation of the proposed algorithm is the local nature of its model, which can be improved by increasing the density of the distributed computing units. Another issue is related to the scalability of its discretized configuration space. Note that for 3D spaces, the method can fairly well approximate the sensorimotor model, yet for multiple DOF (e.g., more than 6) the data is difficult to manage and visualize.
As future work, we would like to implement our adaptive method with other sensing modalities and mechanical configurations, e.g., with an eye-in-hand visual servoing (where the camera orientation is arbitrary) and with variable morphology manipulators (where the link's length and joint's configuration are not known).

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
DN-A conceived the algorithm and drafted the manuscript. JQ and JZ performed the numerical simulation results. AC analyzed the theory and revised the paper. All authors contributed to the article and approved the submitted version.