Fast Model Identification via Physics Engines for Data-Efficient Policy Search

1 Introduction

Reinforcement learning (RL) typically requires a lot of training data before it can provide useful skills in robotics. Learning in simulation can help reduce the dependence on realworld data but the learned policy may not work due to model inaccuracies, also known as the "reality gap". This paper presents an approach for model identification by exploiting the availability of off-the-shelf physics engines [Erez et al., 2015] for simulating robot and object dynamics in a given scene. One objective is to achieve data-efficiency for RL by reducing the need for many real-world robot experiments by providing better-simulated models. It also aims to achieve time efficiency in the context of model identification by reducing the computational effort of this process.

The accuracy of a physics engine depends on several factors. The model assumed by the engine, such as Coulomb's law of friction, may be limited or the numerical algorithm for solving the differential equations of motion may introduce errors. Moreover, the robot's and the objects' mechanical parameters, such as mass, friction, and elasticity, may be inaccurate. This work focuses on this last factor and proposes a method for identifying mechanical parameters used in physical simulations so as to assist in learning a robotic task.

Given initial real-world trajectories, potentially generated by randomly initializing a policy and performing roll-outs, the proposed approach searches for the best model parameters so that the simulated trajectories are as close as possible to the real ones. This is performed through an anytime black-box Bayesian optimization, where a belief on the optimal model is dynamically updated. Once a model with a high probability is identified, a policy search subroutine takes over the returned model and computes a policy to perform the target task. The subroutine could be a control method, such as a linear-quadratic regulator (LQR), or an RL algorithm that runs on the physics engine with the identified model. The obtained policy is then executed on the real robot and the newly observed trajectories are fed to the model identification module, to repeat the same process.

The question that arises is how accurate the identified model should be to find a successful policy. Instead of spending time searching for the most accurate model, the optimization can stop whenever a model that is sufficiently accurate to find a good policy is identified. Answering this question exactly, however, is difficult, because that would require knowing in advance the optimal policy for the real system.

The proposed solution is motivated by a key quality desired in many robot RL algorithms. To ensure safety, most robot RL algorithms constrain the changes in the policy between two iterations to be minimal and gradual. For instance, both Relative-Entropy Policy Search (REPS) [Peters et al., 2010] and Trust Region Policy Optimization (TRPO) [Schulman et al., 2015] algorithms guarantee that the KL divergence between an updated policy and the one in the previous iteration is bounded. Therefore, one can in practice use the previous best policy as a proxy to verify if there is a consensus among the most likely models on the value of the best policy in the next iteration. This is justified by the fact that the updated policy is not too different from the previous one. Thus, model identification is stopped whenever the most likely models predict almost the same value for the previously computed policy, i.e., if all the high-probability models predict a similar value for the previous policy, then any of these models could be used for searching for the next policy.

Empirical evaluation performed in simulation and on a real robot show the benefits of the framework. The initial set of experiments are performed in the OpenAI Gym [Brockman et al., 2016] with the MuJoCo simulator 1 . They demonstrate that the proposed model identification approach is more timeefficient than alternatives, and it improves the data-efficiency of policy search algorithms, such as TRPO. The second part is performed on a real-robot manipulation task. It shows that learning using a simulator with the identified model can be more data-efficient than other model-based methods, such as PILCO, and model-free ones, such as PoWER.

2 RelatedWork

This work relates to model-based RL, where a dynamical system is physically simulated. Model-based learning involves explicitly learning the unknown system dynamics and searching for an optimal policy. Variants have been applied in robotics [Dogar et al., 2012; Lynch and Mason, 1996; Scholz et al., 2014; Zhou et al., 2016], where using inaccurate models can still allow a policy to steer an RC car [Abbeel et al., 2006]. For general-purpose model-based RL, the Probabilistic Inference for Learning COntrol (PILCO) method was able to utilize a small amount of data to learn dynamical models and optimal policies [Deisenroth et al., 2011]. Replacing gradient-based optimization with a parallel, black-box algorithm (Covariance Matrix Adaptation (CMA)) [Hansen, 2006], Black-DROPS can be as data-efficient as PILCO without imposing constraints on the reward function or the policy [Chatzilygeroudis et al., 2017]. CMA was also used for automatically designing open-loop reference trajectories [Tan et al., 2016]. Trajectory optimization was performed in a physical simulator before real experiments, where robot trajectories were used to optimize simulation parameters.

1 MuJoCo: www.mujoco.org

The proposed method also relates to learning from simulation. Utilizing simulation to learn a prior before learning on the robot is a common way to reduce real-world data needs for policy search and provides data-efficiency over PILCO [Cutler and How, 2015]. Policy Improvement with REsidual Model learning (PI-REM) [Saveriano et al., 2017] has shown improvement over PILCO by utilizing a simulator and only modeling the residual dynamics between the simulator and reality. Similarly, using data from simulation as the mean function of a Gaussian Process (GP) and combining with Iterative Linear Quadratic Gaussian Control (ILQG), GP-ILQG is able to improve incorrect models and generate robust controllers [Lee et al., 2017].

In addition, the proposed method utilizes Bayesian optimization (BO) to efficiently identify the model. For direct policy search, a prior can also be selected from a set of candidates using BO [Pautrat et al., 2017]. BO can also be used to optimize the policy while balancing the trade-off between simulation and real data [Marco et al., 2017], or learn a locally linear dynamical model, while aiming to maximize control performance [Bansal et al., 2017].

Traditional system identification builds a dynamics model by minimizing prediction error (e.g., using least squares) [Swevers et al., 1997]. There have been attempts to combine parametric rigid body dynamics model with nonparametric model learning for approximating the inverse dynamics [Nguyen-Tuong and Peters, 2010]. In contrast to these methods, this work uses a physics engine, and concentrates on identifying model parameters instead of learning the models from scratch. Optimizing the simulated model parameters to match actual robot data can also help model robot damage [Bongard et al., 2006]. Recent work performed in simulation proposed model identification for predicting physical parameters, such as mass and friction [Yu et al., 2017].

Different from model-based methods, model-free methods search for a policy that best solves the task without explicitly learning the system's dynamics [Sutton and Barto, 1998; Kober et al., 2013]. They have been shown effective in robotic tasks [Peters et al., 2010]. Policy learning by Weighting Exploration with the Returns (PoWER) [Kober and Peters, 2009] is a model-free policy search approach widely used for learning motor skills. The TRPO algorithm [Schulman et al., 2015] has been demonstrated in some MuJoCo simulator environments. End-to-end learning is an increasingly popular framework [Agrawal et al., 2016; Wu et al., 2015; Byravan and Fox, 2017; Finn and Levine, ], which involves successful demonstrations of physical interaction and learning a direct mapping of the sensing input to controls. These approaches usually require many physical experiments to effectively learn. Overall, model-free methods tend to require a lot of training data and can jeopardize robot safety.

A contribution of this work is to link model identification with policy search. It does not search for the most accurate model within a time budget but for an accurate enough model, given an easily computable criterion, i.e., predicting a policy's value function that is not too different from the searched policy. This idea can be used with many policy search algorithms that allow for smooth changes during learning.

3 Proposed Approach

This section provides first a system overview and then the main algorithm, focusing on model identification.

3.1 SystemOverview and Notations

The focus is on identifying physical properties of robots or objects, e.g. mass and friction, using a simulator. These physical properties are concatenated in a single vector and represented as a D-dimensional vector q ∈ Q, where Q is the space of possible physical parameters. Q is discretized with a regular grid resolution. The proposed approach returns a distribution P on a discretized Q instead of a single point q ∈ Q given the challenge of perfect model identification. In other terms, there may be multiple models that can explain an observed movement of the robot with similar accuracy. The approach is to preserve all possible explanations and their probabilities.

The online model identification takes as input a prior distribution Pt , for time-step t ≥ 0, on the discretized parameter space Q. Pt is calculated based on the initial distribution P0 and a sequence of observations (x0,m0,x1,m1, . . . ,xt-1,mt-1,xt). For instance, xt can be the state of the robot at time t and mt is a vector describing a vector of torques applied to the robot at time t. Applying mt results in changing the robot's state from xt to xt+1. The algorithm returns a distribution Pt+1 on the models Q.

The robot's task is specified by a reward function R that maps state-actions (x,m) into real numbers. A policy p returns an action m = p(x) for state x. The value V p (q) of policy p given model q is defined as V p (q) = � H t=0 R(xt , mt ), where H is a fixed horizon, x0 is a given starting state, and xt+1 = f (xt ,mt ,q) is the predicted state at time t +1 after simulating action mt in state xt given parameters q. For simplicity, the focus is on systems with deterministic dynamics.

3.2 Main Algorithm

Given a reward function R and a simulator with model parameters q, there are many ways to search for a policy p that maximizes value V p (q). For example, one can use Differential Dynamic Programming (DDP), Monte Carlo (MC) methods, or if a good policy cannot be found with alternatives, execute a model-free RL algorithm on the simulator. The choice of a particular policy search method is open and depends on the task. The main loop of the system is presented in Algorithm 1. This meta-algorithm consists in repeating three main steps: (1) data collection using the real robot, (2) model identification using a simulator, and (3) policy search in simulation using the best-identified model.

3.3 Value-GuidedModel Identification

The process, explained in Algorithm 2, consists of simulating the effects of actions mi on the robot in states xi under various values of parameters q and observing the resulting states ˆ xi+1, for i = 0, . . . ,t. The goal is to identify the model parameters that make the outcomes of the simulation as close as possible to the real observed outcomes. In other terms, the following black-box optimization problem is solved:

t ← 0;

Initialize distribution P over Q to a uniform distribution; Initialize policy p;

repeat

Execute policy p for H iterations on the real robot, and collect new state-action-state data {(xi, mi, xi+1)} for i = t, . . . , t +H -1;

Run Algorithm 2 with collected state-action-state data and reference policy p for updating distribution P; Initialize a policy search algorithm (e.g. TRPO) with p and run the algorithm in the simulator with the model argmaxq∈QP(q) to find an improved policy p ′ ; p ←p ′ ;

until Timeout;

Algorithm 1: Main Loop

wherein xi and xi+1 are the observed states of the robot at times i and i+1, mi is the action that moved the robot from xi to xi+1, and f (xi,mi,q) = ˆ xi+1, the predicted state at time i+1 after simulating action mi in state xi using q.

High-fidelity simulations are computationally expensive. Thus, when searching for the optimal parameters in Eq. 1, it is important to minimize the number of evaluations of function E. This is achieved through Entropy Search [Hennig and Schuler, 2012], which explicitly maintains a belief on the optimal parameters. This is unlike other Bayesian optimization methods, such as Expected Improvement, which only maintains a belief on the objective function. The following description explains how this technique is adapted so as to provide an automated stopping criterion.

The error function E does not have an analytical form but is learned from a sequence of simulations with a small number of parameters qk ∈ Q. To choose these parameters efficiently in a way that quickly leads to accurate parameter estimation, a belief about the actual error function is maintained. This is a probability measure over the space of all functions E : RD →R, represented by a Gaussian Process (GP) [Rasmussen and Williams, 2005] with mean vector m and covariance matrix K. The mean m and covariance K of the GP are learned from data: {(q0,E(q0)), . . . , (qk,E(qk))}, where qk is a vector of physical robot properties, and E(qk) is the accumulated distance between observed states and simulated states given qk.

The probability distribution P on the identity of the best physical model q ∗ , returned by the algorithm, is computed from the learned GP as

where H is the Heaviside step function, i.e., H(E(q ′ )- E(q)) = 1 if E(q ′ ) ≥ E(q) and H(E(q ′ )-E(q)) = 0 otherwise, and pm,K(E) is the probability of E according to the learned GP mean m and covariance K. Intuitively, P(q) is the expected number of times that q happens to be the minimizer of E when E is a function distributed according to the GP.

Input: state-action-state data {(xi, mi, xi+1)} for i = 0, . . . , t a discretized space of possible values of physical properties Q, a reference policy p, minimum and maximum number of evaluated models kmin, kmax, model confidence threshold h, value error threshold e ; Output: probability distribution P over Q; Sample q0 ∼Uniform(Q); L← /0; k ←0; stop← f alse; repeat // Calculating the accuracy of model qk lk ← 0; for i = 0 to t do Simulate {(xi, mi)} using a physics engine with physical parameters qk and get the predicted next state ˆ xi+1 = f (xi,mi,qk) ; lk ← lk + ‖ ˆ xi+1 - xi+1‖2; end L ← L ∪ {(qk , lk)}; Calculate GP(m,K) on error function E, where E(q ) = l, using data (q , l) ∈ L; // Monte Carlo sampling Sample E1, E2, . . . , En ∼ GP(m, K) in Q; foreach q ∈ Q do P(q) ≈ 1 n n � j=0 1q=argmin q ′∈Q E j (q ′ ) (2) end // Selecting the next model to evaluate qk+1 = argminq ∈Q P(q ) log (P(q )) ; k ← k + 1; // Checking the stopping condition if k ≥ kmin then q ∗ ← argmaxq ∈Q P(q ); Calculate the values V p (q) with all models q that have a probability P(q) ≥ h by using the physics engine for simulating trajectories with models q ; if �q,∀P(q)≥h P(q)|V p (q)-V p (q ∗ )| ≤ e then stop ← t rue; end end if k = kmax then stop ← t rue; end until stop = true;

The distribution P from Eq. 3 does not have a closedform expression. Therefore, Monte Carlo (MC) sampling is employed for estimating P. The process samples vectors [E(q ′ )]q ′ ∈Q containing values that E could take, according to the learned GP, in the discretized space Q. Then P(q) is estimated by counting the ratio of sampled vectors of the values of simulation error E where q happens to make the lowest error, as indicated in Eq. 2 in Alg. 2.

Finally, the computed distribution P is used to select the next vector qk+1 to use as a physical model in the simulator. This process is repeated until the entropy of P drops below a certain threshold, or until the algorithm runs out of the allocated time budget. The entropy of P is given as �q∈Q-Pmin(q) log (Pmin(q)). When the entropy of P is close to zero, the mass of distribution P is concentrated around a single vector q , corresponding to the physical model that best explains the observations. Hence, next qk+1 should be selected so that the entropy of P would decrease after adding the data point (qk+1,E(qk+1))} to train the GP and re-estimate P using the new mean m and covariance K in Eq. 3.

Entropy Search methods follow this reasoning and use MC again to sample, for each potential choice of qk+1, a number of values that E(qk+1) could take according to the GP in order to estimate the expected change in the entropy of P and choose the parameter vector qk+1 that is expected to decrease the entropy of P the most. The existence of a secondary nested process of MC sampling makes this method impractical for online optimization. Instead, this work presents a simple heuristic for choosing the next qk+1. In this method, called Greedy Entropy Search, the next qk+1 is chosen as the point that contributes the most to the entropy of P,

This selection criterion is greedy because it does not anticipate how the output of the simulation using qk+1 would affect the entropy of P. Nevertheless, this criterion selects the point that is causing the entropy of P to be high. That is a point qk+1 with a good chance P(qk+1) of being the real model, but also with a high uncertainty P(qk+1)log ( 1 P(qk+1) ). Initial experiments suggested that this heuristic version of Entropy Search is more practical than the original Entropy Search method because of the computationally expensive nested MC sampling loops used in the original method. The actual sampled values are not restricted to lie on the discretized grid. Once a grid value is selected as the candidate with the highest entropy, a local optimization process using L-BFGS takes place to further optimize the parameter value.

The stopping condition of Alg. 2 depends on the predicted value of a reference policy p. The reference policy is one that will be used in the main algorithm (Alg. 1) as a starting point in the policy search with the identified model. That is also the policy executed in the previous round of the main algorithm. Many policy search algorithms (such as REPS and TRPO) guarantee that the KL divergence between consecutive policies p and p ′ is minimal. Therefore, if the difference |V p (q)-V p (q ∗ )| for two given models q and q ∗ is smaller than a threshold e, then the difference |V p ′ (q)-V p ′ (q ∗ )| should also be smaller than a threshold that is a function of e and KL(p‖p ′ ). A full proof of this conjecture is the subject of an upcoming work. In practice, this means that if q and q ∗ are two models with high probabilities, and |V p (q)-V p (q ∗ )| ≤ e then there is no point in continuing the optimization to find out which one of the two models is actually the most accurate because both models will result in similar policies. The same argument could be used when there are more than two models with high probabilities.

4 Experimental Results

VGMI is evaluated both in simulation and on a real robot.

4.1 Experiments on RL Benchmarks in Simulation

Setup: The simulation experiments are performed in OpenAI Gym with the MuJoCo simulator. The space of unknown physical models q is:

Inverted Pendulum (IP): A pendulum is connected to a cart, which moves linearly. The dimensionality of Q is two, one for the mass of the pendulum and one for the cart.

Swimmer: The swimmer is a 3-link planar robot. Q has three dimensions, one for the mass of each link.

Hopper: The hopper is a 4-link planar mono-pod robot. Thus, the dimensionality of Q is four.

Walker2D: The walker is a 7-link planar biped robot. Thus, the dimensionality of Q is seven.

HalfCheetah: The halfcheetah is a 7-link planar cheetah robot. The dimensionality of the space Q is seven.

The simulator starts the ground truth mass, which is only used for rollouts to generate data for model identification. For inaccurate simulators as priors, the mass was randomly increased or decreased by 10 to 15%. All the policies are trained with TRPO using rllab [Duan et al., 2016]. The policy network has 2 hidden layers with 32 neurons each. VGMI is compared against a) Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [Tan et al., 2016] and b) Least Square (LS) optimization using L-BFGS-B [Swevers et al., 1997]. All of them optimize the objective function of Eq. 1.

The number of iterations for policy search varies on problem difficulty. The main loop in Alg. 1 is executed for 20 iterations for IP, 100 iterations for Swimmer, 100 iterations for Hopper, 200 iterations for Walker2D and 400 iterations for HalfCheetah. VGMI follows Alg. 2 and is executed every 10 iterations, i.e., H = 10 in Alg. 1. The same iterations are used for CMA and LS and they are also given the same time budget for model identification. All the results are the mean and stand deviation of 20 independent trials.

Results: Performance is reported for total training time, i.e., including both model identification and policy search. Fig. 2 shows the cumulative reward per trajectory on the ground truth system as a function of total time. VGMI always achieves the highest reward across benchmarks and earlier than CMA and LS. CMA performs better than LS, though relies on a large number of physically simulated trajectories.

Fig. 3 assesses the automated stopping of VGMI. The red line is the error threshold e in the predicted value. The green line is when VGMI stops based on e. The same algorithm was also executed with increasing pre-set time budgets (xaxis). The reality gap in predicting the value of the best policy was empirically estimated using the best-identified and the ground truth model. This ensures that VGMI does not stop prematurely. When it stops (green line), the actual error (blue curve) is indeed below the predefined threshold e.

Fig. 4 shows that the adaptive stopping criterion of VGMI performs better than any fixed time budget for model identification. Thus, the proposed stopping criterion effectively balances the trade-off between model identification and policy search. Fig. 5 demonstrates that Greedy Entropy Search (GES) in the context of VGMI results in improved convergence relative to the original Entropy Search.

4.2 Real Robot Pushing Experiments

Setup: Consider Figs. 6 and 7, where Motoman assists (Baxter) to pick up a bottle. The object is known but can be reached only by Motoman. The intersection of the robots' reachable workspace is empty and Motoman must learn an action to push the bottle 1m away within Baxter's reachability. If Motoman simply executes a maximum velocity push, the object falls off the table. Similarly, if the object is rolled too slowly, it can get stuck in the region between the two robots. Both outcomes are undesirable as they require human intervention to reset the scene.

The goal is to find an optimal policy with parameter h representing the pushing velocity. The pushing direction is towards the target and the hand pushes the object at its geometric center. No human effort was needed to reset the scene. A pushing velocity limit and a pushing direction were set so the object is always in Motoman's workspace. The approach iteratively searches for the best h by uniformly sampling 20 different velocities in simulation, and identifies the object model parameters q ∗ (the mass and friction coefficient) using trajectories from rollouts by running VGMI as in Alg. 2. VGMI is executed after each rollout, i.e., H = 1 in Alg. 1. The method is compared to PoWER [Kober and Peters, 2009] and PILCO [Deisenroth et al., 2011]. For PoWER, the reward function is r = e-dist , where dist is the distance between the object position after pushing and the desired target. For PILCO, the state space is the 3D object position.

Results: Two metrics are used for evaluating performance: 1) The distance between the final object location after being pushed and the desired goal location; 2) The number of times the object falls off the table. Fig. 8 shows that in the realworld experiments the method achieves both lower final object location error and fewer number of object drops, which is important for robot learning that minimizes human effort. The model-free PoWER results in higher location error and more object drops. PILCO performs better than PoWER as it learns a dynamical model but the model is not as accurate as the VGMI one. Since simple policy search is used for VGMI, the performance is expected to be better in more advanced policy search methods, such as combining PoWER with VGMI. A video with the experiments can be found on https://goo.gl/Rv4CDa.

5 Conclusion

This paper presents a practical approach that integrates physics engines and Bayesian Optimization for model identification to increase RL data efficiency. It identifies a model that is good enough to predict a policy's value function that is similar to the current optimal policy. The approach can be used with any policy search algorithm that guarantees smooth changes in the learned policy. It can also help the real-world applicability of motion planners that reason over dynamics and operate with physics engines [Bekris and Kavraki, 2008; Littlefield et al., 2017]. Both simulated and real experiments show that VGMI can decrease the number of rollouts needed for an optimal policy. Future work includes analyzing the method's properties, such as expressing the conditions under which the inclusion of model identification reduces the need for physical rollouts and the speed-up in convergence.

Acknowledgments

This work was sponsored by NSF IIS-1734492, IIS-1723869, a NASA ECF award to Dr. Bekris and U.S. Army Research Lab Collaborative Agreement # W911NF-10-2-0016.