M ANY popular frameworks for controller design are based on the robot's model of dynamics. In the real world, however, this model can often turn out to be inaccurate, due to, for example, misspecification of the robot's physical parameters, mechanical wear and tear, and deployment-time interventions such as additional payload. While a well designed controller is robust to small inaccuracies in the dynamics, large deviations may significantly degrade its performance.
Our goal is to make corrections to the model behind the controller during deployment, through online learning using onboard sensors. Since the nature of a model is to predict the future given the past, data for supervised learning of dynamics can be collected automatically without human supervision, as time goes on and the future is revealed.
Because data are generated along the controller's trajectory that we are trying to improve, they might not contain enough information about the entire system. Nevertheless, we find it sufficient to limit the scope of learning to a local neighborhood of the current point in the current trajectory, instead of the entire system, if the learned model is updated in real time as the trajectory evolves. Fortunately, even globally complex systems, such as the highly nonlinear hybrid systems for legged locomotion, can be locally simple. Therefore, we also find it sufficient to learn with only a time-varying, locally linear model, which is computationally feasible to be updated in real time.
We first develop the intuition of online learning into a method for controllers that drive the outputs to the desired behavior based on control-affine models. We then analyze this method's theoretical properties, and evaluate it in two applications for legged locomotion.
A. Related Work 1) System identification: Given a system with known form but unknown parameters, system identification (sysID) estimates these parameters from signals given by the system ( [1]). Recent papers have applied sysID for inertial parameters of a humanoid ( [2], [3]). The parameters are assumed to be constant in time, and estimation is performed before the deployment of a controller. Thinking of identification as training and deployment as testing, sysID trains a model before deployment, and keeps the model fixed during testing. Since the goal is to model the system's behavior globally across the entire state space, sysID usually requires driving the system to diverse enough states, using diverse enough inputs. This requirement is known as persistence of excitation in control theory, and might be difficult to satisfy without many samples from the plant. In contrast, we only model the system's behavior locally, around the small neighborhood of our current state, learning a linear model even for complex systems with relatively few samples.
2) Learning dynamics: There is also a developing community in machine learning, modeling dynamics of the environment from interactions and observations ( [4], [5], [6], [7], [8], [9], [10], [11]). It has roughly the same goal as sysID, but often uses powerful tools from deep learning, and does not assume any specific form of the system; here, learning often produces a general prediction model. We diverge from this community in the global vs. local aspect (like from sysID), but embrace its philosophy of learning a general model with parameters that might not be interpretable.
3) Adaptive control: The intuition of adaptive control is to change the controller's parameters during deployment ( [12]). Online system identification ( [13]) is the most relevant subfield, since it directly concerns the model behind the controller. It has been successfully applied in manipulation ( [14], [15], [16], [17], [18]), and for the location and inertial parameters of the center of mass of a quadruped ( [19]). For online sysID, the parameters considered are very specific, and estimation relies on the physics of the model and the particular controller for the application. Our work considers parameters in a much more general sense closer to that of the machine learning community. Our parameters are functions of the state, thus are inherently time-varying and abstract. In fact, in the controlaffine form, every term of our dynamics is updated in real time as the state evolves. Furthermore, unlike sysID (online or not) whose goal is to identify the parameters, our goal is simply to give accurate predictions for the next timestep, again closer to the goal of learning. This allows our method to not rely on the specific meanings of the parameters and instead work with general model-based controllers. Another relevant subfield is L1 adaptive control ( [20], [21]), which, like our work, concerns the residual dynamics, but does not use learning.
4) Online learning: Our work performs supervised learning online, which has long been a subject of research in machine learning ( [22], [23], [24]). The two central questions are: where does the label come from, and how is learning evaluated. Traditionally ( [25]), learning has been evaluated with regret, and labels can come from a potentially adversarial oracle. Recently, the computer vision community has been using self-supervised tasks to provide labels ( [26], [27], [28], [29], [30]), and the continual learning community has been evaluating with forward and backward predictions ( [4]) c.f. Subsection II-B.
In this paper, vectors (a, α) are bold and lowercase, matrices (A, Ω) are bold and uppercase, scalars and functions (of all type signatures) are not bold. We assemble matrices and vectors like in MATLAB: [A, B] concatenates A and B horizontally with a comma, and [A; B] concatenates them vertically with a semicolon. 0 n denotes the n × n matrix of zeros, and 1 n denotes the n × n identity matrix. Also, • denotes the 2-norm for vectors (Euclidean norm) and matrices (spectral norm), unless stated otherwise. We express quantities in the nominal dynamics ᾱ with a bar, in the residual dynamics α with a tilde, and in the true (plant) dynamics α without anything on top.
Given a robotic system that is characterized by rigid-body dynamics, we denote x ∈ R n as its state, u ∈ R m its vector of control inputs, and y ∈ R d its vector of outputs. The output dynamics can almost always be written as a secondorder system of the following form ( [31]), known as controlaffine ( [32]):
We consider model-based controllers whose goal is to drive the vector of tracking errors η = [y; ẏ] to zero. The bars on top of the variables imply that they come from our assumed nominal model, which in reality can never be completely accurate. The unknown real-world dynamics are called the true (plant) model, denoted without the bars as α, β. We often use an alternative set of notations to write equation (1) simply as:
in order to emphasize the role of ᾱ and β as time-varying parameters of the output dynamics.
To make corrections to the nominal model, we incorporate two residual parameters and obtain the following form:
where α is called the weight and β is called the bias. They are written as time-varying parameters, and have the same dimensions as ᾱ and β respectively. The tildes on top of them emphasize that they are estimated from data.
To better understand these residual parameters, we manipulate equation (3) into:
Intuitively, the above equation says that the goal of learning is to make the residual model on the right-hand side account for the prediction errors of the nominal model on the left-hand side. It also reveals the role of labels vs. covariates, as we explain next in the context of learning.
For real systems, sensor data can only be collected at discrete sampling intervals. We denote each sampling timestep by an integer subscript, which converts equation ( 4) into:
Note that we are merely sampling a continuous system at discrete timesteps, so continuous-time concepts such as acceleration are still well defined. We collect a dataset of the form D := {label s , covariate s } s=t-k,...,t-1 , where s is the index of discrete timesteps, and k denotes the fixed size of the sliding time window. From equation ( 5), we have
Given a dataset, our method solves regularized least squares a.k.a. ridge regression on the labels and covariates. The weight of the solution is αt , and bias is βt . Note that in textbook-style least squares, the weight is a vector, and the label and bias are scalars; for our learning problem, the weight is a matrix in R d×m , and the label and bias are vectors in R d . But we can simply reduce this to d independent vector-scalar least squares problems. The same regularization is added independently to these d problems, since they share the same covariates; thus inversion of the covariance matrix, the most computationally costly step, is only performed once.
The solved parameters are then immediately used by the model-based controller to produce u t . In both of our later examples, the baseline controller solves for u t in an applicationspecific optimization problem with the assumed nominal parameters ᾱt and βt . We simply substitute these with ᾱt + αt and βt + βt respectively, as shown in Figure 2.
Learning is performed online, as the controller is running with the learned parameters. At the beginning, all residual parameters are initialized to zero, because there is not enough data to learn them. Once we are k steps into the trajectory, we have enough data to form D as above and solve for the residual parameters; informed by them, the controller generates an improved trajectory, which in turn generates new data that are more relevant as time goes on.
The fact that D only keeps the k most recent data points implements a natural forgetting mechanism. In reinforcement learning terms, D is called the replay buffer, which stores the off-policy data that are not generated by the current controller; in our case, data in D are generated by the old controllers using the residual parameters from previous timesteps. Because we learn small, local models, we encourage forgetting so that our model capacity can be used only for the neighborhood of our current state. This is in contrast to the vast literature in reinforcement learning [33], [24], [34], [4], where the goal is to learn a large, global model; there the replay buffer contains as much historical data as possible, and various techniques are implemented to discourage forgetting.
Our method can also be viewed as bootstrapping from a "bad" controller based on an inaccurate model to a better one. This might not be feasible, however, if the initial model deviates too much from the plant. For example, if the nominal model is so far off that the robot loses balance immediately, no useful information will be contained in the data collected. Fortunately, when deviations happen gradually over time, there will more likely be enough information for learning to maintain a controller that keeps generating useful data. We study this phenomenon empirically in Section IV.
Suppose the true (plant) output dynamics is control-affine:
We prove that our method stabilizes the tracking errors under two assumptions. The main theorem illustrates our intuition of learning in a local time window under smoothly varying dynamics, and characterizes the role of k, our window size. Denote errors in the nominal model's prediction as
with αt := α t -αt , and βt := β t -βt . Denote the prediction of the residual model as
Assumption 1: The model-based controller can stabilize the tracking errors η = [y; ẏ] if for some > 0,
In words, Assumption 1 says that the proposed model-based controller works when the proposed (nominal plus residual) model is relatively accurate; Assumption 2 says that the deviations in dynamics are relatively smooth (in the space of parameters) over time.
In addition, we denote the motor torque saturation as u < B. Denote u t = [u t ; 1] ∈ R m+1 , and
We set k ≥ m + 1, so σ min (U ) > 0, i.e. the covariance matrix of ordinary least squares (OLS) has rank m + 1. Theorem 1: Given the above assumptions, if
then the model-based controller stabilizes η.
Note that any claim of stability in Theorem 1 is completely inherited from the baseline controller, when Assumption 1 holds. Our method is agnostic to the exact type of stability e.g. exponential / asymptotic, which depends on the underlying baseline, and is orthogonal to the theory we develop.
In Theorem 1, B, d, δ α and δ β are constants determined by the application. is the model-based controller's tolerance for model inaccuracy, also independent of our method. The only quantity we tune is k, the window size, which strongly affects σ min (U ). With a large k, we pay a factor of k √ k, intuitively due to the lag in our dataset. With a small k, we pay for the decrease in σ min (U ), as α and β become more sensitive to noise. The user should tune k to find a sweet spot in the middle. In practice, we use regularized least squares instead of OLS, so σ min (U ) is always > 0 and more noise tolerant, making the balance less delicate w.r.t. choice of k. We use k = 100 in both of our applications (100 and 200 ms respectively).
Before proving Theorem 1, we state two lemmas, whose proofs are given in Subsection A of the appendix.
Lemma 1:
. Let y t = A t u t for t = 1, ..., k, and à be the OLS estimator of the dataset {(y 1 , u 1 ), ..., (y k , u k )}. If for t = 1, ..., k + 1, A t+1 -A t < δ A , and u t < B, then
where
By definition, Ãt is the least squares solution on D. We then apply Assumption 1 and Lemma 2 to finish the proof.
We now apply our method to two model-based controllers, derived from two different perspectives for different robotic platforms: a Lyapunov perspective to control the full-order dynamics of bipedal robots, and a simplified dynamics based control architecture for robust quadrupedal locomotion. We focus on identifying the components of our method in the context of each controller, without elaborating on derivations of the nominal dynamics.
Let q be the robot's configuration, and x = [q; q] be the robot's state. We define y = h(x), where h is called the virtual constraints ( [35]). For a biped, stabilizing η = [y; ẏ] means, for example, that the torso maintains a constant posture, and the legs walk in a scissor-symmetric gait.
The nominal output dynamics, whose derivation we omit, can then be written in the familiar form of equation ( 1), using Lie-derivatives of the nominal dynamics in the state space as ᾱ and β:
where D is the inverse of the mass-inertia matrix, C is the Coriolis matrix and ḡ is the gravity vector. While ᾱ(x) might not be square (d = m) in general, this particular bipedal controller has the same number of virtual constraints as actuated joints. Now the control law
We can then design v to stabilize the output dynamics using control Lyapunov functions (CLFs), a common tool in control theory for providing stability guarantees in legged locomotion ( [36]). Because η is linear in η and v, it is straightforward to find a CLF by solving the Lyapunov equation V (η) ( [37]). It is then a well known fact that V (η, v) < -cV implies exponential stability of η(t), with a constant c > 0. This motivates the following CLF-based quadratic program (CLF-QP) to solve for v:
where u min and u max are bounds of the torque saturation constraints. Since the output dynamics is already in the form of equation ( 1), it is straightforward to apply our method to obtain α and β. We can then modify the C2 in the optimization problem (17) to have
In Section IV we show that this simple modification leads to significant improvements under uncertain dynamics.
To control a quadrupedal system walking stably under large disturbance (such as heavy loads), we take the model predictive control (MPC) approach using the simplified dynamics from [38] as our baseline controller.
For quadrupedal dynamics, let p, ṗ, p ∈ R 6 be the position, velocity and acceleration of the robot's center of mass (CoM). Let f i ∈ R 3 be the ground reaction force at the robot's i th foot, with i ∈ {1, 2, 3, 4}. We also denote f = [f 1 ; f
The nominal dynamics of the CoM is given by
where ḡ ∈ R 6 is the gravity vector, D ∈ R 6×6 is the inverse mass matrix, and G ∈ R 6×12 is called the grasp map, which depends on the robot's state and is assumed to be accurate.
The goal of the model-based controller is to have p and ṗ track the desired position and velocity p d and ṗd , generated from user command. In Sec.II notations, y = pp d , we want to stabilize η = [y, ẏ] around zero. This is achieved by having p track some desired acceleration pd , generated from PD control on p d and ṗd . The model-based controller then uses equation (19) to solve for f :
stance and swing leg constraints, friction pyramid condition.
(
where more details can be found in [39], Following the outline in Section II, we modify (19) to incorporate the linear residual model:
where D is the weight, and -g is the bias. Note that the nominal dynamics in (19) has no Coriolis terms, a simplification often adopted in the literature for model-based controller design of quadrupeds with small angular velocity. While this simplification has been validated in many implementations, it is never completely accurate. Therefore, even if D = D and ḡ = g i.e. they are both accurate parameters, ( 19) is still an inaccurate description of the plant. We make no distinction, philosophically or algorithmically, between unknown dynamics e.g. payload, and unmodeled dynamics e.g. the Coriolis terms discarded by design. Our true output dynamics can take any general form. Also note that Assumption 1 is in fact not satisfied by our baseline controller due to its simplifications e.g. massless legs. In this case, stability is left to empirical validation.
Moving on, we sample equation ( 21) at discrete timesteps:
and form the dataset as
After solving for D and ḡ, we use them to modify the objective function in equation ( 20) as:
By definition, D + D must be positive definite; this is also necessary for the optimization problem above to make sense. For computational efficiency, we solve for D unconstrained, and find that our least squares solution in fact always gives D + D positive definite for our experiments.
Video of our experiments is available at https://youtu.be/Je 2Y-FQpKw ([40]). Simulations are performed in the PyBullet ( [41]) physics engine.
Our baseline controller discussed in Subsection III-A is taken from [33], which introduces its own setting and method for unknown dynamics. We perform simulation in their setting, and make comparison with their method.
The problem setting is based on RABBIT ( [42]), an underactuated planar five-link bipedal robot with seven degree-offreedom; virtual constraints and controller design are based on [37]. Model uncertainty is introduced in [33] by scaling the mass of each link by a factor of two in the real environment. The baseline CLF-QP controller falls in a few steps in this setting, due to the significant difference in dynamics between the nominal and true model.
By querying the plant, [33] uses model-free reinforcement learning (RL) to train a policy that directly adds on the original control inputs u, without reasoning about the unknown dynamics in the model space. Specifically, the commanded control inputs take the form u + u θ (x), where u is a neural network policy with parameters θ. Their reward is designed 3. Bipedal walking with mass of each link scaled by two. Both our method and that of [33] walk stably. Their RL-based method trains on 20,000 samples from the real environment before deployment. Our method trains completely online and does not sample from or anticipate the real environment, treating it as truly unknown until the robot is deployed, and results in smaller impulses of control inputs and better tracking performance. The top panel visualizes the gait generated by our method. to encourage V < -cV , where the value of V is obtained by simulating in the plant. After 20,000 samples from the plant simulated using the true dynamics, their method trains a policy which walks in the true dynamics without falling.
Our method walks stably in the same setting, training completely online without querying the plant at all before deployment. In fact, Fig. 3 shows that our method enjoys smaller impulses of control inputs and better tracking performance than the RL-based method, even though the latter had privileged access to the plant before deployment to optimize exactly for these metrics.
Online learning enables us to treat the plant as truly unknown, in terms of both data and mathematical representation, while only the latter is unknown for methods that train offline like in [33]. This philosophical difference prevents our controller from overfitting on the training environment. In particular, our controller still walks stably under the original dynamics without scaling, where the policy trained with the scaled links fails, because it overfits to the scaled dynamics.
In addition, our controller walks stably in all environments below, where the baseline and the RL-based method cannot:
1) scaling the control inputs by half, in order to simulate transmission inefficiencies and motor wear and tear; 2) scaling the mass of the torso by four, in order to simulate payload on the back of the humanoid; 3) scaling the mass of the right leg by four, as an example of asymmetric changes in dynamics. We keep the same hyper-parameters for all the experiments above, including a windows size of 100 ms (where k = 100 and each timestep is 1 ms). The robot is still able to walk The baseline has completely fallen in 2 s, but the proposed method still walks stably after 10 s (50 kg). The bottom visualizations are captured when the payload reaches the specified mass. The torque limit is reached at 25 kg.
under the scaled dynamics with a window size of 10 or 1000, but has higher norm of control inputs and tracking errors.
Our baseline controller, as discussed in Subsection III-B, is based on [38] and used subsequently in [39] and [43]. Our implementation is modified from the publicly available code of [43] on an Unitree A1 quadruped, and keeps their original parameters unless stated otherwise. The A1 weighs 12 kg and has 12 motors, three for each leg, with the stated torque limit of 35.5 Nm. We experiment in PyBullet using Unitree's URDF description, and also on a real robot. In both simulation and real world, we use a window size of k = 100 (like for the biped); the controller runs at a frequency of 500 Hz, making the dataset window 200 ms.
We command robot to walk with linear velocity of 0.5 m / s in the x-direction, while maintaining CoM height of 0.24 m. Both the baseline and the proposed method can walk stably without payload, while tracking the desired velocity and height. With 6 kg of payload, however, the baseline can barely walk at 2 / 3 the desired velocity, and sags to 2 / 3 the desired height; the robot falls with 7 kg.
The proposed method walks stably with 12 kg of payload (same as its body mass), while tracking the desired velocity up to 0.05 m / s, and the desired height up to 0.01 m; all motors torques are less than 35.5 Nm. With more than 12 kg, however, tracking becomes less accurate, and with 15 kg the robot falls. Since the payload is carried from very the beginning of simulation, the robot visibly sags for the first fifth of a second, as we collect data before we can estimate the residual parameters. With 12 kg it soon recovers from the sag, but for larger payloads it struggles to get back. Next, we experiment with gradually changing dynamics. We start with an empty payload, and increase its mass by 5 kg / s, that is, 0.001 kg per timestep, once simulation begins. The tracking errors are shown in Figure 4. The baseline falls within 2 s. We have tried to improve the baseline by tuning the PD gains for pd , but found it ineffective. This observation is reasonable, since larger gains only make pd more aggressive, but cannot help if the model-based controller fails to achieve it using the nominal dynamics. The proposed method walks stably even when the payload reaches 50 kg. Motor torques reach the specified limit at 25 kg (5 s), but the URDF allows simulation to keep running.
To facilitate hardware testing, we fit the Unitree A1 quadruped with a loading rig designed to hold up three standard 1 inch weight plates. The rig allows for incremental, discrete changes in load while the quadruped is in operation. The rig itself weighs 0.9 kg.
The experiments were designed to compare the performance of the baseline and proposed controllers under varying load conditions during operation. Two tests for each controller were performed: a step-in-place test and a 0.1 m/s forward motion test. The load conditions for the tests are shown in Table 1. Due to the manual loading process, the duration of each load varies by a small amount of transition time, typically less than 1 s. To protect the hardware from possible damage, we do not load beyond 10 kg, and limit operation at this load to 5 s.
In the transition from simulation to hardware, we had to address the problem of acceleration estimation from noisy measurements. [43] uses a Kalman filter to fuse IMU and
We have introduced a method to update model parameters through online learning, while a controller is running with past, inaccurate versions of these parameters. Unlike methods in adaptive control, our method uses learned parameters that are general functions of the state, thus inherently time-varying. To the best of our knowledge, this is the first method that applies machine learning to real-time updates (500 to 1000 Hz) in hardware walking experiments.
While the nominal models in our applications are derived from classical mechanics, our method can be applied to any black box nominal model e.g. a simulator. While our baselines are derived from classical control principles, our method can also be applied to any controller using the black box model, even a policy trained in simulation. We hope to explore these potentials in future work, under broader definitions of unknown dynamics, such as sim-to-real transfer. 2) Proof of Lemma 2: We first prove the vector version in a claim, which is used in the proof of the lemma.
Claim 1: Consider y t = u t , a t ∈ R for t = 1, . . . , k, and u t , a t ∈ R m . Suppose u t < B and a t+1 -a t < . Let ã be the OLS estimator on this dataset, then
.
Proof: Define the feasible set of weights for a dataset as A = {(a 1 , a 2 , • • • , a k ) : y t = u t , a t , a t -a t+1 ≤ , ∀t} .
Then a k can only exist in the kth component of A, denoted
Our goal is to bound max a k ∈A k ãa k . Define e t = a t -a k , where e k = 0 and e t+1 -e t < . We can rewrite A using these conditions as A = {(a 1 , • • • , a k ) : y t = u t , a t , e t = a t -a k , e k = 0, e t+1 -e t < , ∀t}.
Define E = [e T 1 ; ...; e T k ] and A = [a T 1 ; ...; a T k ]. Also, by definition of OLS, ã = (U T U) -1 U T y. Therefore
where • denotes the Hadamard operator, σ min (U) is the minimum non-zero singular value of U; the last equality follows from singular value decomposition of U. Note that
which finishes the proof of Claim 1. Now we extend the result of Claim 1 to prove Lemma 2, Note that in the context of Lemma 2, A t ∈ R d×m , and is different from the definition of A in the proof of Claim 1. We use the standard matrix norm relationship ( [44])
for any matrix X ∈ R d×m . Combining the second half of equation (25) with the lemma's assumption, we have
As explained in Subsection II-B, the matrix-vector least squares problem is solved by reducing to d independent vectorscalar sub-problems, for each dimension of y t . Each subproblem solves for one row of Ã. From equation ( 26), we already know that rows of the ground truth weight matrices satisfy the smoothness assumption in Claim 1. Therefore we can apply Claim 1 to each row of Ã, yielding
Combining this with the first half of equation ( 25) finishes the proof of Lemma 2.