Data-driven learning and control of nonlinear system dynamics: A robust-learning approach via Sontag’s control formula

Nonlinear Dyn https://doi.org/10.1007/s11071-024-10149-1 ORIGINAL PAPER Data-driven learning and control of nonlinear system dynamics A robust-learning approach via Sontag’s control formula Yeyson A. Becerra-Mora ·José Ángel Acosta Received: 27 December 2023 / Accepted: 8 August 2024 © The Author(s) 2024 Abstract This work falls into the field of discovering the dynamic equations of stabilizable nonlinear systems, via a learning-and-control algorithm to process the data sets of trajectories previously obtained. To this end,aninterlacedmethodtolearnandcontrolnonlinear system dynamics from a set of demonstrations is proposed, under a constrained optimization framework for the unsupervised learning process. The nonlinear system is modeled as a mixture of Gaussians and Sontag’s formula together with its associated Control Lyapunov Function is proposed for learning and control. Lyapunov stability and robustness in noisy data environmentsareguaranteed,asaresult of theinclusionofcontrol in the learning-optimization problem. The performances are validated through a well-known dataset of demonstrations with handwriting complex trajectories, succeeding in all of them and outperforming previous The authors acknowledge support from the Project 2024/00000362 for emerging research group Multi-Robot And Control Systems (MACS) under the framework VII PPIT-US 2024 and by the National Program for Doctoral Formation (Minciencias-Colombia, 885-2020). Y. A. Becerra-Mora Department of Electronic Engineering, CUN, Bogota 111711, Colombia Y. A. Becerra-Mora (B)·J. Á. Acosta DepartmentofIngenieríadeSistemasyAutomática,Universidad de Sevilla, 41092 Sevilla, Spain e-mail: [email protected]; ye[email protected] J. Á. Acosta e-mail: [email protected] methods under bounded disturbances, possibly coming from inaccuracies, imperfect demonstrations, or noisy datasets. As a result, the proposed interlaced solution yields a good performance trade-off between reproductions and robustness. Therefore, this work sheds some more light on the automatic discovery of nonlinear dynamics from noisy raw data. Keywords Data-driven dynamics ·Control Lyapunov functions ·Unsupervised machine learning 1 Introduction Nowadays, the amount and variety of real-world data sets available and the ease of obtaining and recording them have caused data-driven approaches to grow exponentially. An especially interesting approach is to discover the underlying dynamics of a real system with the only information collected in the data set, either defining or discovering evolution equations. If that process is reliable, it opens the door to a tremendous number of applications ranging from unknown complex dynamics to high-dimensional dynamics. However, this entire process is useful if and only if it is automatic i.e. carried out by a dedicated computer. Currently, aided mainly by the rise of branches of artificial intelligence, there are many sophisticated methods to carry out this automated discovery of complex dynamics. Although they work reasonably well in most situations, they are not completely reliable–both at a math- 123 Y. A. Becerra-Mora, J. Á. Acosta ematical and computer level–and require some supervision or human intervention. The work presented here focuses on the mathematical level with emphasis on disturbed/noisy data sets, which are the most common real data in practice and noise is usually ignored. Many approaches have been proposed to learn and then predict complex dynamics in nonlinear systems. Machine Learning methods like supervised or unsupervised learning have aided to identify such dynamics as well as make estimates of the system’s output. Unlike mathematical models-based approaches, Machine Learning requires a dataset to learn from the system where inputs and outputs need to be identified. Someapplicationshaveaimedatdetectinganomaliesin networks [1,2], predicting the weigh category of a person [3], predicting two-dimensional trajectories from human walks [4], among many other applications. Furthermore, some studies have addressed not only the learning approach, but also the system stability as in [5,6]. Among all applications, the human-robot collaboration has taken on great relevance for the scientific community, because it has the two necessary ingredients to study novel data-driven techniques for nonlinear/complex systems: the highly non-linear dynamics of human-robot interaction while executing some complex task, and the ease of recording trajectory data sets with a controlled laboratory setup. In essence, in this application, a human user teaches a robot a complex task through many in-situ demonstrations, obtaining a dataset of trajectories followed. This data is processed with a data-driven mathematical method that attempts to discover the dynamics of a nonlinear system that reproduces these trajectories as faithfully as possible. Finally, the extracted dynamics can be programmed into the robot so that it follows these trajectories without human intervention. Furthermore, these data sets can be made available to the entire scientific community so that different researchers can test and verify new techniques, as e.g. the one used in this work. Thus, new methods to estimate the robotic paths (inherently nonlinear) have to be considered. Learning from Demonstrations (LfD) [7] is a paradigm that has been very well received recently to estimate (stable) dynamical systems (DS), hence, several works have used this approach. Some works related to LfD have been characterized by the use of regression methods like Gaussian Process Regression (GPR), Locally Weighted Projection Regression (LWPR), Gaussian Mixture Regression (GMR), and polynomial regression, among others. Another current is based on the optimization of policies as the Dynamic Movement Primitives (DMP) method. The approach presented in this work falls into the former. Modeling nonlinear trajectories of DSs from user demonstrations is possible by employing Gaussian Mixture Model (GMM) and GMR [8], with the combination of Hidden Markov Model (HMM) and GMR [9], using GMM and HMM [10], with hypernetworks and Neural Ordinary Differential Equations (NODEs) [11], or even employing GMM parameters into a Gaussian Process Regression (GPR) [12]; nevertheless, neither of them ensures stability. On the other hand, the Binary Merging method proposed in [13], ensures local stability through a numerical process and employs GMM/GMR to build a trajectory estimate. Likewise, [14] use GMM to estimate a nonlinear point-to-point trajectory and define conditions to ensure Lyapunov stability; as well as [15] propose a constrained optimization problem to learn the unknown parameters of GMM that ensures asymptotic stability. In [16] a constrained optimization problem is employed to learn a Control Lyapunov Function (CLF) from demonstrations to ensure stability on an estimated DS. Following a similar approach, an improved estimate of the DS is presented in [17] which is based on transforming the data from demonstrations with a diffeomorphic candidate. Diverse applications have been developed by using a DS estimate approach with regression methods and stability conditions such as: Obstacle avoidance, [18]; multiple attractors, [19]; capturing irregular objects in the air, [20]; rehabilitation therapy, [21]; visual serving, [22]; robotic handwriting, [23]; impedance control, [24]; virtual reality, [25]; among others. Other approaches for LfD have been addressed. Recurrent Neural Networks (RNN) have been used to learn a Lyapunov function from the user demonstrations like in [26] where the learned function is used to stabilize the system. It is worth mentioning that this work and the one in [16] have resembled approaches, both of them divide the problem into two steps: Learning a function and using it to somehow stabilize the system. Nonlinear autoregressive polynomial models and a constrained least-square estimator are employed in[27]tolearnaDSwhichislocallyasymptoticallystable. In [28] a feed-forward neural network is employed to ensure stability and to estimate the region of attrac- 123 Data-driven learning and control of nonlinear system dynamics tion with the learned Lyapunov function. Finally, based on energy tank-based controllers, the work in [29] ensures the stability of a DS dissipating its energy. Ourapproachisinspiredbytwopreviousworks,Stable Estimator of Dynamical systems (SEDS) [15] and Control Lyapunov Function-based Dynamic Movements (CLF-DM) [16]. The former learns optimal parameters from GMM to estimate a complex trajectory through GMR. In contrast, the latter learns an energy function and corrects the dataset using a virtual control in combination with a regression method, e.g. GMR. Moreover, SEDS learns the energy function and estimates the trajectory in one step, while CLF-DM do the same but in two different steps; nevertheless, CLFDM achieves a better accuracy of estimates than SEDS. Both methods have been validated in the LASA handwriting dataset [30] wherein unperturbed trajectories are presented1. In our study, we propose a hybrid learning algorithm called Control Lyapunov Function-based Stable Estimator(CLF-SE)thatcombinestheadvantagesofSEDS and CLF-DM, additionally, improves the accuracy of estimates, and deals with noisy data. i.e. our method optimizes the GMM parameters and uses a control signal to enforce global exponential stability, and robustness. As in [15,16], our work also uses LASA for the sake of having a fair comparison. Nevertheless, noise levels have been added to the dataset to validate the response of our method under external disturbances. Therefore,acontrollawisdesignedtocompensatesuch perturbations, unlike previous approaches. At this point, we would like to distinguish our line of research with the one of obtaining and/or approximating solutions of nonlinear dynamics described either by Partial Differential Equations (PDEs) or nonlinear Ordinary Differential Equations (ODEs), as e.g. in [31,32] and references therein. In that line of research, the PDEs/ODEs are known (possibly by first principles), unlike in our approach that they are completely unknown. Hence, our approach is not focused on finding numerical or analytical approximation of solutions, but rather on discovering the unknown dynamics from a known dataset. Finally, this work somewhat complements the thorough analysis carried out in [33], where in short there were two main findings about algorithms made to 1LASA dataset repository at https://bitbucket.org/khansari/ lasahandwritingdataset/src/master/. discover governing equations from data: all of them present drawbacks and reveal deficiencies; and they have a tendency to remember information that confirms preconceptions. This leads to the conclusion that, in general, their outcomes are not reliable and equivalent to the results obtained with black boxes. From our own experience we mostly agree with conclusions in [33], however we have found some other outcomes that add new constructive viewpoints to the discussion. On one hand, we first would like to highlight that there are many similarities between classical identification/estimation techniques and what nowadays is called data-driven ones. Obviously the datadriven-learning-like approaches have opened many fruitful distinct viewpoints that made them deserve ‘new labels’; however there are clear overlaps between the classical and new approaches. In this way, it has long been demonstrated that there is a deep knowledge in identifying nonlinear systems via matching data with physics first principles. Roughly speaking, there are two approaches for the discovery of nonlinear dynamics: manually and automatically; being the former the classical way and the latter the modern one in which this article is focused. Let us comment a couple of examples of the former from our own experience. Thus, we recall the intense hard work, where it certainly took us more than a year to obtain a reliable model from data, based on wellknown physical principles of Mechanics. In the end, we could construct a model to reproduce non-linear behaviors from data, such as homoclinic orbits, for a laboratory pendulum (see a summary in [34] and references therein), or more recently in real wind turbines in [35]. Of course, these are two simple examples of a large number in the field, but quite significant in terms of the nonlinearities. Therefore, the data-driven discovery can be performed with accuracy manually, but underagreatdegreeofcarefuldatafilteringandmining, i.e. supervision. Unfortunately, the amount of manual effort required was too high and would not be useful in many applications. The reliability problem becomes when we want to automate the process by ‘reducing’ it to unsupervised learning. Thus, the work proposed here sheds some more light on the automatic discovery of nonlinear dynamics from noisy raw data and demonstrates that there is an inevitable trade-off, specially in the following: 123 Y. A. Becerra-Mora, J. Á. Acosta Automation trade-off. The less knowledge of nonlinear dynamics to be discovered the poorer estimate; however, a higher level of automation. In this way, we fix a quite flexible structure for nonlinear dynamics, rather than letting the algorithm discovers this as well. Realistic trade-off. The inherent freedom in the solution provides worse estimates due to overfitting and/or overparametrization, both indeed greatly affected by noise. However, adding a compatible goal to the dynamics improves drastically the results, at the expense of performing some interaction with the dynamics to enforce the goal via control and to reject the noise disturbance. The paper is organized as follows. The problem statementanddetailedcontributionsareinSect.2;Sect. 3summarizes GMM and GMR; Sect. 4is devoted to stability and robustness and Sect. 5the learning. Sect. 6validates the performance with a dataset and Sect. 7 conclusions and future research. Notation. Collection of scalars, vectors or matrices are denotedby enclosingthem in {...}.Forx,˙x∈Rd,consider Msets of Npoints each, a pair is defined as the collection {xm,n,˙xm,n}, where the superscript integers m,ndenote a point n⊂Nof the set m⊂M. For compactness,symbolsin ‘bold’ denoteagrouped collection as e.g. pairs x:= col(x,˙x)and at a specific point as xm,n:= {xm,n,˙xm,n}.∇xV(x)and ∇2 xV(x)denote the gradient and the Hessian with respect to xof the scalar function V(x), respectively. |·|stands for vector norm, and |·| Wa weighted norm with matrix W0 where ‘’ stands for positive definite. A function α(s)∈K∞ if α(0)=0 and lim|x|→∞ α(|s|)=∞. For a matrix A wedenotewiththesuperscript’s’the symmetric matrix As=A+A. 2 Problem statement and contributions The underlying idea is to accurately estimate a ‘singular’(nonlinear)dynamics,whichisrepresentedthrough demonstrations, but where their mathematical models are completely unknown, because either are too complex to be obtained by first principles or the ultimate goal is to automate a complex task. Thus, only (noisy) input/outputdata areavailable through demonstrations. First, let us be precise with the distinction between demonstrationsandtrainingsetsusedinthelearningliterature. Demonstrations reproduce trajectories of presumably stable (nonlinear) dynamics, while training sets do not necessarily. The proposed approach learns and controls and hence it covers both. Let the trajectory to be learned the state vector x∈R2d, which encodes a set of Mdemonstrations of Ndata-points each. Denote the pairs D:= {xm,n,˙xm,n}M,N m=1,n=1which corresponds to the set of instances of an autonomous dynamical system given by ˙x=f(x)+u,(1) where f:Rd→ Rdis a continuous and smooth vector field, uthe control action of appropriate dimensions. In order to facilitate notation, let us denote in ‘bold’ each pair xm,n:= {xm,n,˙xm,n}and hence D:= {xm,n}M,N m=1,n=1. The pair estimate mechanism is basedonconditionalprobabilityusing known data xm,n as prior probabilities and posterior ones to estimate unknown ˙xm,n. Thus, for example, for path planning control the pair encodes position and velocity, respectively. Let the learned controlled system be defined as ˆ ˙x=ˆ f(x)+ˆu(x)+η(t,x), (2) where ˆ f:Rd→ Rdis a nonlinear estimate and ˆu(x) the input estimate and ηa bounded additive disturbance accounting for inaccuracies in measurements, imperfectdemonstrationsdue to boundednoiseor evenexternal disturbances coming from the original dynamics (1). Without any loss of generality, we consider the target at x∗=0, and hence f(0)=ˆ f(0)=0asin[15]. Noticethat definition(1) accountsalso foruncontrolled (open-loop) systems u=0, and hence stabilizable2 ones, but not necessarily open-loop stable. Therefore, the definition (2) allows the estimate ˆ f(x)to represent open-loop unstable dynamics, unlike the related works. This is the prelude to the following and only assumption of the original dynamics (1). Assumption 1 The system (1) is stabilizable at the target. Before enumerating the contributions let us contextualize this work in the vast field of Machine Learning (ML). We explain this with the help of the ‘big picture’ sketch shown in Fig.1. Thus, this approach falls into the automatic discovery of nonlinear system dynamics with the following ingredients: 2The system (1) is stabilizable if for any x(0)there exists u(t) such that limt→∞ x(t)=x∗,t≥0. 123 Data-driven learning and control of nonlinear system dynamics Fig. 1 Big picture of the approach (i) Problem statement Automatic discovery of system dynamics, where for ease of discussion, we explicitly show the learning parameter θ, so that (2) become ˙x(t)=ˆ f(x(t), θ)+ˆu(·,θ)   =: ˆ F(x(t),θ) +η, where, the dynamical governing equations ˆ fare unknown; ˆuare the control inputs (unlike most learning approaches); and ηare external disturbances, that to facilitate the discussion we divide intoη:= ηn+ηdwith ηnthenoisy-likedisturbance and ηdthe unexpected non-learned disturbance. (ii) Input The available information from the system is a dataset of trajectories points {x(t), ˙x(t)}for a timesequence(recallthattrainingdatasetpointsare denoted as {·}). Notice that, ηncorrupts the dataset affecting the learning (see Fig.1); (iii) Learning The parameters θof the dynamics and controller such that min θ{˙x}− ˆ f({x},θ)+ˆu({x},θ) 2 . s.t.stability and robustness guarantees This objective function is represented by J(x,θ)in Fig.1). The learning method GMM/GMR and stability/robustness conditions imposed through Lyapunov theory and control. (iv) Output The complete dynamics and control ˆ F(·) such that their solutions, namely X(t), reproduce theoriginaldatasetwithminimal error,i.e.{˙x(t)}≈ ˆ F({x(t)}), and any prediction can be generated with the explicit dynamics ˙ X(t)=ˆ F(X(t)). (v) Objective Reproduction/prediction control methods with stability and robustness guarantees, with the following requirements: I. Explicit dynamics for ˆ f(·)(analytical expression), needed for control design methods. II. Controller ˆu(·)such that guarantees stability and robustness as follows: i. Minimal error in reproductions, being able to reject noisy-like disturbances ηnin the learning procedure. ii. Robust and stable predictions, being able toreject unexpecteddisturbances ηdwhich have not been learned. Note how ηddisturbs the predictions in Fig.1. Notice that, even though (i) and (ii) are usual ingredients in most ML methods (without η), their main objective is the reproduction and/or prediction without putting effort into either the guarantees or disturbances. In the learning approach we propose, the guarantees and rejection of disturbances are imposed through the objectivefunction with control and through constraints, in (iii). Although most of ML methods have good performances in their predictions, in general, most of them do not consider either control, i.e. u≡0, or stability during learning. Moreover, most of them provide data asoutput in(iv),withoutanalytical expressionsfor ˆ f(·) limiting their use to guarantee robustness and stability methods, which is our main objective (v). In summary, the method proposed here, provides both: good repro- 123 Y. A. Becerra-Mora, J. Á. Acosta ductions of the learned dynamics (v.I) even with noisylike data inputs and the capability of rejecting disturbances, even those that have not been learned (v.II), under rigorously established robustness/stability guarantees on the discovered nonlinear dynamics. Therefore, this work complements the ML methods, going a step further than obtaining precise automatic solutions. The combined ML and control method will be more reliable and will smooth out human–machine interactions in the broadest sense. Once the problem is posed, we are in position to describe the contributions with respect to related works inthese directions: C1Stabilizability; C2Optimization and C3 Robustness. C1. The closest and seminal related works which have served as inspiration for the present one are [15] and [16]. First, with respect to the class of systems, in there the class of systems and so its estimate are restricted to be stable, which is, in fact, a hard constraint imposed in the optimization problem. In here, Assumption 1 does indeed relax that condition by requiring only stabilizability rather than stability of (1). To give a simple example, consider the dynamics ˙x=x+uas (1), which is an open-loop unstable system and hence the approach[15]or[16] cannot beused. However, it isstabilizable and thus the approach proposed in here can learn and estimate ˆ f(x)and ˆu(x), ensuring stability and convergence to the target. C2. In [16], the approach is split into two independent steps: 1) an estimate of a Lyapunov function, namely V(x), in agreement with the demonstration set, and 2) an out-the-loop correction of such set of demonstrations through u(x)to enforce the already estimated V(x)to be a Lyapunov function for ‘mislearned’ points. To avoid any misunderstanding, let us clarify that although the latter might seem a remedy to impose stability, the fact is that in [16] the learning strategy priories the estimation over the execution, unlike here. Essentially, they discard and replace the ‘mislearned’ points of ˆ f(x)while keeping the solution V(x)obtained in a decoupled optimization problem. Certainly, the main advantage is to solve two simplified and independent optimization problems for V(x) and ˆ f(x),however, thedisadvantageisthat V(x)might turn out not to be an accurate Lyapunov function for the original set of demonstrations, hence the authors propose such suitable correction. Unlike here, we propose to solve only one optimization problem in the loop for the demonstration set including V(x),ˆ f(x)and ˆu(x), which indeed is computationally lighter because it does not force the stable constraint, thus allowing to be directly used in real time for the dynamics of the form (1), i.e. without any low-level controller. C3. In this approach, we consider undesirable disturbances all collected in ηof (2), which is a major difference with related works. The definition of dynamics to be learned (1) considering the drift as the sum of f(x)and u(x)allow us to guarantee robustness. Thedisturbancescan comefrom verydifferent sources: eitherinaccuraciesinmeasurements,imperfectdemonstrations due to bounded noise and even external disturbances coming from the original dynamics. As a result, we prove that the inclusion of the appropriate u(x)in the learning-optimization problem guarantees goodperformance trade-off between reproductions and robustness in noisy data environments. 3 Regression via GMM and GMR The trajectory estimate ˆ ffrom a set of demonstrations Mis defined as a mixture modeling with a finite set of Gaussian kernels, K. Mixture modeling builds a (coarse) representation of data density through a fixed numberofmixturecomponents.MethodslikeBayesian Information Criterion (BIC) [36], Akaike Information Criterion (AIC) [37], and Deviance Information Criterion(DIC)[38]canbeusedtofindanoptimalnumberof components.GMMisdefinedasanunsupervisedlearning or clustering algorithm which is based on Expectation Maximization (EM) to optimize the parameters of the Gaussian mixture [39], likewise K-means [40] is used to initialize EM parameters. Although there are optimalregressionmethodslikeGaussianProcess [41], they suffer the curse of dimensionality, i.e. the number of data points force to a linear growing pattern in the regressor. Thus, GMM is able to construct the estimate ˆ fwith the solely unknown parameters given by priors πk, meansμkand covariance matrices kof each Gaussian kernelk,where3θk:= {πk,μ k, k}withk=1, ..., K, and μk=μx k μ˙x k,k=x kx˙x k ˙xx k˙x k.(3) 3Note that θkis defined as a collection of vectors and matrices with dimension dim(θk)=2(K×d2). 123 Data-driven learning and control of nonlinear system dynamics Hence, in GMM the probability that a data-point in the trajectory fits into a Gaussian is defined by the conditional probability density function P(xm,n|k)as a Gaussian normal distribution given by P(xm,n|k)=N(xm,n;μk, k) =(2π)2ddet(k)−1 2 exp −1 2|xm,n−μk|(k)−1, and the mixture probability density function becomes P(xm,n;θk)= K  k=1 πkP(xm,n|k)⎧ ⎪ ⎨ ⎪ ⎩ k∈1,...,K, n∈1,...,N, m∈1,...,M. An estimated velocity ˆ f(x)is described through GMR [42]as ˆ f(x)= K  k=1 γk(x)(μ˙x k+˙xx k(x k)−1(x−μx k)), (4) where we have defined the nonlinear weighting term [40] that measures the influence of different Gaussians kas γk(x):= πkP(x|k) K i=1πiP(x|i)∈R[0,1]. Even though ˆ f(x)can follow nonlinear trajectories, it is not possible to ensure asymptotic stability at the target with this method. However, Sect. 4presents the control strategy to deal with that. In Fig.2we show the GMR estimates for different curve types (C, G, W, Sine, S) named by their shapes. Even though, as Kincreases—shown downwards at the top of Fig.2—, the estimates of nonlinear trajectories are more similar to the demonstrations, simple shapes like C need a smaller Kthan complex like W to find a better fit of the trajectory. In any case, note that no trajectory reaches the target, as expected (i.e. this method is not able to ensure asymptotic stability). Furthermore, disturbances can affect the performance of the estimate. To overcome these issues a data-driven control strategy is presented in the next section. 4 Robust data-driven control with Sontag’s formula For the control action, we select the universal Sontag’s formula which is one of the main differences with relatedworks.Recall that the underlying ideais theparity between the existence of Lyapunov functions and the stabilisability with some special properties introduced in the work of Artstein [43] and subsequently in the work of Sontag [44]. Artstein proved that there exists a stabiliser that makes the origin of (2) asymptotically stable if and only if there exists a C1Lyapunov function V(x)>0, with V(0)=0, satisfying the inequality inf u∈U{∇ xV(x)( ˆ f(x)+ˆu(x))}<0(5) foratleastanuif x= 0 with U⊆Rl. Sontag named such Control Lyapunov Function (CLF). Artstein also proved that, although smooth elsewhere, such uis not necessarily continuous at the origin, and provided this necessary and sufficient condition on Vto make ucontinuous: for every >0 there is a δ>0 that whenever |x|<δ,x= 0, there is some uwith |u|<such that ∇xV(ˆ f+ˆu)<0 holds, named is the small control property later on. Defining a(x):= ∇xV(x)ˆ f(x), b(x):= ∇xV(x)and −→ b(x):= b(x)/|b(x)|2,an analytical solution of the optimal control problem (5) reads ˆu(x)=−(a(x)+ρ(|x|))−→ b(x), a(x)+ρ(|x|)>0, 0,a(x)+ρ(|x|)≤0,(6) where we dropped the arguments for compactness. Note that (6) is a simplified version of the original Sontag’s formula. Certainly, the Lyapunov function to be learned has a global minimum at x∗=0, which means that b(x∗)=0 and b(x)= 0forx= x∗, and hence ruling out the case b(x)=0 in the original formula. More importantly, Sontag proposed as optimal4 ρ(|x|)=ρ0a(x)2+|b(x)|4,ρ0>0, but in general it could be any positive one with ρ(0)=0. This is the main difference with [16] where minimum or no control effort was preferred forcing the 0-case in (6) and class Kfunction ρ(|x|)=ρ0(1−e−κ0|x|),ρ0,κ 0>0, losing the robustness that we prove in this work. In any case, as in any Lyapunov control design for nonlinear systems the finding of a CLF is quite difficult. In this work, we fix the CLF as an asymmetric biquadratic function of the error as suggested in [16,17]. It is readily seen that this is a limitation of the classes of systems that can be learned but it is kept in this work 4In [44] the control gain ρ0=1, the adapted case for ρ0>0is in [45]. 123 Y. A. Becerra-Mora, J. Á. Acosta Fig. 2 Demonstrations (dotted black lines), estimates (red lines) and target (blue star) for comparison. Moreover, we also provide the stability vs robustness requirements unlike in there. Thus, let σl:Rd→ Rbe defined as σl(x):= xPl(x−μl), σl(0)=0, for l=1,...,L, and the CLF candidate as V(x)=xP0x+ L  l=1 sign+(σl)σ l(x)2,(7) with sign+(·)a coefficient defined as 0 when its argument is negative and 1 elsewhere. Notice that (7)ispositive definite by definition upon P0>0 irrespective of Pl. Notice that the CLF candidate is not unique, other functions satisfying the stability and CLF requirements can be proposed, which are stated in the Lemma below. However, this choice strikes a good balance between having enough freedom to reproduce a wide range of different trajectories and complexity. The necessary requirements for stability are stated in the following lemma. Lemma 1 Consider V :Rd→ R+defined in (7) with Ps l0,l =1,...,L, then V (0)=0,V(x)>0for all x = 0and x =0is a global minimum. Proof The gradient and Hessian of (7) become ∇xV(x)=2P0x+2 L  l=1 sign+(σl)σl(x)(Ps lx−Plμl), ∇2 xV(x)=2P0+2 L  l=1 sign+(σl)σl(x)Ps l +2 L  l=1 sign+(σl)(Ps lx−Plμl)(Ps lx−Plμl). The necessary condition ∇xV(0)=0 comes from σl(0)=0. Noting that sign+(σl)σl≥0 and ∇xV(x)x>0 by direct calculation, hence the Hessian is positive definite everywhere, which means that (7) is globally convex.  Thefunction(7) allows to model a broad set of complexfunctions,with the asymmetryprovidedby vectors μland the number of asymmetric functions Ldefined by the user. In what follows in this section, we establish stability and robustness of the proposed approach, which optimizes while executing the feedback loop and is thus able to reject disturbances and diminish the inevitable mismatch between the reproduction and the learned dynamics. The main result is stated in the following proposition. Proposition 1 Consider the system dynamics (2) with the state feedback (6) under Assumption 1. For any reproduction the following stability (i) and robustness (ii) results hold: (i) For η=0the origin x =0is a globally exponentially stable equilibrium for ρ0>0and V (x)given by (7) is a Lyapunov function; (ii) Let η, ˙η∈K∞, with |η|2<κfor some constant κ>0. Then closed-loop trajectories are globally uniformly ultimately bounded and converge to the residual set := x∈Rd:ρ(|x|)≤2κ/ρ0. Proof Before proving the claims, the following facts are in order by construction: F1. ByLemma1,thefunction(7)satisfies∇xV(x)x= b(x)x>0forx∈Rd/{0}, and hence there exist functions α, α ∈K∞such that α(|x|)≤V(x)≤ α(|x|)). 123 Data-driven learning and control of nonlinear system dynamics F2. By construction ρ(0)=0⇔x=0 and lim|x|→∞ ρ(|x|)=∞, hence ρ∈K∞. Claim (i) [η=0]. The time derivative of (7) along the trajectories of (2) with (6) becomes ˙ V(x)=∇ xV(x)ˆ ˙x =∇ xV(x)ˆ f(x)−(a(x)+ρ(|x|))−→ b(x) =−ρ(|x|), and then F2 implies ˙ V(x)<0, for x∈Rd/{0}.The result follows, noting that from F1 V(x)is globally positive definite and radially unbounded and V(x)= 0⇔x=0,andthereforetherealwaysexistsaconstant c>0 such that ˙ V(x)<−cV(x)concluding the proof of this claim. Claim (ii) [η= 0]. Likewise, an estimate of the time derivative of (7)for(2) with (6) becomes ˙ V(x)=∇ xV(x)ˆ ˙x =∇ xV(x)ˆ f(x)+η(t,x)−(a(x)+ρ(|x|))−→ b(x) ≤−ρ(|x|)+|∇ xV(x)η(t,x)| ≤−ρ0a(x)2+|b(x)|4+ρ0|b(x)|2 2+|η(t,x)|2 2ρ0 =−ρ0 4(2a(x))2+|b(x)|4+|η(t,x)|2 2ρ0 ≤−1 2ρ(|x|) 2−κ ρ0, where in the fourth line we have made use of Young’s inequality and made some rearrangements in the fifth one. Recalling F2, |η|2<κimplies that ˙ V(x)<0for x∈Rd/. Therefore, F1 and F2 together with the set being compact, we conclude global boundedness of trajectories. Finally, the ultimately claim comes from F1 ensuring compact level-sets of V(x)and ˙ V(x)< 0forx∈Rd/, and for any initial condition, i.e. uniformly, concluding the proof.  As it is readily seen, the stability margin of (i) naturally rejects bounded disturbances, and it could be interpreted as an input–output feedback problem with additive model uncertainty. Remark 1 The class Kfunction ρ(|x|)=ρ0(1− eκ0|x|),ρ0,κ 0>0, proposed in [16] does not provide any robustness to disturbances or inaccuracies. Moreover, in [16] the authors seek for minimal or no interaction of virtual control (ρabove) to guarantee good reproductions. However, in this work we demonstrate that a good trade-off is achieved between robustness and good reproductions. To show the difference in Fig.3the control effect along the trajectories in all shapes is presented, note that demonstrations are free of noise. However, a trade-off is always needed and a discussion about this is carried out at the end of Sect. 4. A final remark about the smoothness of the trajectories is in order. Remark 2 This work focuses on smooth vector fields motivatedmainlyfrompathsduetomovementstowards a target position, which are the usual e.g. in vehicles of any type, robots and so on. Thus, the theory provided has not been stated for discontinuities in general. However, by construction the proposed controller (6) can deal paths with a discontinuity at the target, and the theory provided covers this case as well. In order to generalize the approach to discontinuous paths elsewhere, the controller must be changed, possibly to a piecewise controller. 5 Learning by optimization algorithm In this section, we define a learning algorithm to determine the optimal parameters, which is based on a constrained optimization problem. Recall that the parameters are θk:= {πk,μ k, k}for GMR, and Pland μl for the energy function. Let us define the whole set of parameters as θ:= {θk,Pl,μ l}. The parameters for GMR are initialized with GMM which in turn uses EM, while the parameters to estimate the energy function, Pland μlare initialized with identity matrices and null vectors, respectively. We propose an eager learning in which the algorithm tries to optimize all parameters at the same time, unlike [16]. Even though this is an ambitious approach, we manage to consolidate two independent steps into just one, mainly because we removed the hardest stability constraint (20b) of [15]. However, the number of parameters to be optimized can increase considerably because of the number of Gaussian kernels Kand asymmetric quadratic functions Lutilized to estimate the nonlinear trajectory as well as to try to make it asymptotically stable. The objective function is defined as the mean square error (MSE) between the data and estimates. The employed method to solve 123 Y. A. Becerra-Mora, J. Á. Acosta tation. In: 2017 IEEE International conference on systems, man, and cybernetics (SMC), pp. 2225–2230 (2017) 22. Paolillo, A., Saveriano, M.: Learning stable dynamical systems for visual servoing. In: 2022 International conference on robotics and automation (ICRA), pp. 8636–8642, (2022) 23. Göttsch, P., Olschewski, R., Werner, H.: A segmentation scheme for CLF dynamic movement control applied to robotic handwriting. IFAC-PapersOnLine, 20th IFAC World Congress. vol. 50, pp. 11459–11464, (2017) 24. Khansari-Zadeh, S.M., Khatib, O.: Learning potential functions from human demonstrations with encapsulated dynamic and compliant behaviors. Auton. Robot. 41, 45– 69 (2017) 25. Hernoux,F., Nyiri, E.,Gibaru, O.: Virtual realityfor improving safety and collaborative control of industrial robots. In: Proceedings of the 2015 virtual reality international conference, pp. 1–6 (2015) 26. Lemme, A., Neumann, K., Reinhart, R.F., Steil, J.J.: Neural learning of vector fields for encoding stable dynamical systems. Neurocomputing 141, 3–14 (2014) 27. Santos, Rafael F., Pereira, Guilherme A.S.., Aguirre, L.A.: Learning robot reaching motions by demonstration using nonlinear autoregressive models. Robot. Auton. Syst. 107, 182–195 (2018) 28. Rego, Rosana C.B.., de Araújo, Fabio, Meneghetti, U.: Lyapunov-based continuous-time nonlinear control using deep neural network applied to underactuated systems. Eng. Appl. Artif. Intell. 107, 104519 (2022) 29. Saveriano, M.: An energy-based approach to ensure the stability of learned dynamical systems. In: 2020 International conference on robotics and automation (ICRA), pp. 4407– 4413, (2020) 30. Khansari-Zadeh, S.M., Billard, A.: Handwriting dataset, (2011) 31. Khater, Mostafa M.A..: Computational simulations of propagation of a tsunami wave across the ocean. Chaos, Solitons Fractals 174, 9 (2023) 32. Khater, Mostafa M.A..: Soliton propagation under diffusive and nonlinear effects in physical systems; (1+1)- dimensional MNW integrable equation. Phys. Lett., Sect. A: Gen, Atomic Solid State Phys. 480, 8 (2023) 33. Antonelli, G., Chiaverini, S., Di Lillo, P.: On data-driven identification: is automatically discovering equations of motion from data a chimera? Nonlinear Dyn. 111(7), 6487– 6498 (2023) 34. Acosta, J.Á.: Furuta’s pendulum: a conservative nonlinear model for theory and practise. Math. Probl. Eng. 2010, 742894 (2010) 35. Zarouala, Rachad Oulad Ben., Acosta, J.Á.: Timescale separation via Rayleigh quotient in flexible wind turbines: a singularlyperturbed approach.Nonlinear Dyn. 97(4), 2723– 2738 (2019) 36. Schwarz, G.: Estimating the dimension of a model. Ann. Stat. 6(2), 461–464 (1978) 37. Akaike,H.:Anewlookatthestatisticalmodelidentification. IEEE Trans. Autom. Control 19(6), 716–723 (1974) 38. Spiegelhalter, D.J., Best, N.G., Carlin, B.P., Van Der Linde, A.: Bayesian measures of model complexity and fit. J. R. Stat. Soc.: Ser. B (Stat. Methodol.) 64(4), 583–639 (2002) 39. Dempster, A.P., Laird, N.M., Rubin, D.B.: Maximum likelihood from incomplete data via the EM algorithm. J. Roy. Stat. Soc. 39(1), 1–38 (1977) 40. Christopher, M.B.: Patter Recognition and Machine Learning. Springer, New York (2006) 41. Rasmussen, C.E., Williams, C.K.I.: Gaussian Processes for Machine Learning. The MIT Press, Cambridge (2006) 42. Cohn, D.A., Ghahramani, Z., Jordan, M.I.: Active learning with statistical models. J. Artif. Intell. Res. 4, 129–145 (1996) 43. Artstein, Z.: Stabilization with relaxed controls. Nonlinear Anal.: Theory, Methods Appl. 7(11), 1163–1173 (1983) 44. Sontag, E.D.: A ‘universal’ construction of Artstein’s theorem on nonlinear stabilization. Syst. Control Lett. 13(2), 117–123 (1989) 45. Acosta, J.Á., Dòria-Cerezo, A., Fossas, E.: Stabilisation of state-and-input constrained nonlinear systems via diffeomorphisms: a Sontag’s formula approach with an actual application. Int. J. Robust Nonlinear Control 28(13), 4032– 4044 (2018) 46. Wright, M.H.: The interior-point revolution in optimization: History, recent developments, and lasting consequences. Bull. Amer. Math. Soc. 42, 39–56 (2005) 47. Becerra-Mora, Y., Quesada Conejero, V., Rodríguez Castaño, A., Acosta, J.A.: Trajectories learning by demonstrations for unmanned superficial vehicle. XLIV Jornadas de Automática, 1–5. 2023. https://doi.org/10.17979/spudc. 9788497498609.001 Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 123