Grgur Kovač Flowers

Adrien Laversanne-Finot
email: adrien.laversanne-finot@inria.fr

Pierre-Yves Oudeyer
email: pierre-yves.oudeyer@inria.fr

GRIMGEP: Learning Progress for Robust Goal Sampling in Visual Deep Reinforcement Learning

Designing agents, capable of learning autonomously a wide range of skills is critical in order to increase the scope of reinforcement learning. It will both increase the diversity of learned skills and reduce the burden of manually designing reward functions for each skill. Self-supervised agents, setting their own goals, and trying to maximize the diversity of those goals have shown great promise towards this end. However, a currently known limitation of agents trying to maximize the diversity of sampled goals is that they tend to get attracted to noise or more generally to parts of the environments that cannot be controlled (distractors). When agents have access to predefined goal features or expert knowledge, absolute Learning Progress (ALP) provides a way to distinguish between regions that can be controlled and those that cannot. However, those methods often fall short when the agents are only provided with raw sensory inputs such as images. In this work we extend those concepts to unsupervised image-based goal exploration. We propose a framework that allows agents to autonomously identify and ignore noisy distracting regions while searching for novelty in the learnable regions to both improve overall performance and avoid catastrophic forgetting. Our framework can be combined with any stateof-the-art novelty seeking goal exploration approaches. We construct a rich 3D image based environment with distractors. Experiments on this environment show that agents using our framework successfully identify interesting regions of the environment, resulting in drastically improved performances. The source code is available at https://sites.google.com/view/grimgep.

Introduction

Recent work in reinforcement learning has shown that robots can learn complex individual skills such as grasping [START_REF] Kalashnikov | Qt-opt: Scalable deep reinforcement learning for vision-based robotic manipulation[END_REF], locomotion [Haarnoja et al., 2018a[START_REF] Yang | Data efficient reinforcement learning for legged robots[END_REF], and manipulation tasks [START_REF] Gu | Deep reinforcement learning for robotic manipulation with asynchronous off-policy updates[END_REF]. However, the standard paradigm in reinforcement learning often involves learning policies specific to each task. This imposes an engineering burden as a reward function has to be manually designed for each task. On the contrary, humans are capable to autonomously learn a wide variety of skills by setting their own goal and practicing them, with little to no supervision. They naturally move from simple tasks in the early infancy to more complex tasks as they grow. Designing autonomous learning agents, capable of setting their own goals, similarly to how human learn is thus a promising avenue of research. Such agents would autonomously discover general purpose policies, that could later be used for specific purpose, reducing the need for outside intervention and manual engineering.

In the absence of a specific objective, an efficient strategy in order to learn a general purpose policy is to learn to reach as many states as possible [START_REF] Ashvin V Nair | Visual reinforcement learning with imagined goals[END_REF][START_REF] Vitchyr | Skew-fit: State-covering self-supervised reinforcement learning[END_REF][START_REF] Nasiriany | Disco rl: Distribution-conditioned reinforcement learning for general-purpose policies[END_REF][START_REF] Laversanne-Finot | Curiosity driven exploration of learned disentangled goal spaces[END_REF][START_REF] Colas | Curious: intrinsically motivated modular multi-goal reinforcement learning[END_REF]. Combining a goal reaching policy with an exploration mechanism, even in its simplest form, can allow solving otherwise impossible environments. For example, GO-Explore [START_REF] Ecoffet | First return, then explore[END_REF] showed that when goal-reaching is assumed solved by expert knowledge (manually setting the environment state to the goal state), performing random actions from interesting goals states allows to solve hard exploration problems. This form of goal exploration involves two challenges: how to learn a goal-reaching policy and how to explore as many novel goals as possible. This first challenge is solved by using recent reinforcement learning algorithms [START_REF] Andrychowicz | Hindsight experience replay[END_REF], Haarnoja et al., 2018b] together with good goal representations [START_REF] Ashvin V Nair | Visual reinforcement learning with imagined goals[END_REF][START_REF] Warde-Farley | Unsupervised control through non-parametric discriminative rewards[END_REF][START_REF] Laversanne-Finot | Curiosity driven exploration of learned disentangled goal spaces[END_REF]. In order to solve the second issue, many approaches try to maximize the diversity of sampled goals, for example by using importance sampling techniques [START_REF] Vitchyr | Skew-fit: State-covering self-supervised reinforcement learning[END_REF] in order to sample states uniformly from the underlying distribution.

While optimizing the diversity of sampled goals has been shown to enable robots to solve complex manipulation tasks [START_REF] Vitchyr | Skew-fit: State-covering self-supervised reinforcement learning[END_REF], it falls short in environments containing unpredictable distractors. This problem is often referred to as the noisy TV problem as an agent aiming at maximizing novelty would get attracted by a TV displaying a sequence of continually novel and unpredictable images [Kaplan andOudeyer, 2007, Burda et al., 2019].

One way to solve this problem is to explore parts of the environment with high learning progress [START_REF] Schmidhuber | Curious model-building control systems[END_REF], Oudeyer et al., 2007], which is called competence progress when it measures progress towards self-generated goals [START_REF] Baranes | Active learning of inverse models with intrinsically motivated goal exploration in robots[END_REF]Oudeyer, 2013, Forestier et al., 2017]. For example, SAGG-RIAC [START_REF] Baranes | Active learning of inverse models with intrinsically motivated goal exploration in robots[END_REF] gradually clusters the environment into sub-regions of different learning progress and samples goals in the more promising regions (e.g. the regions with high learning progress). This approach has been shown to enable learning of inverse models in high-dimensional robots. In other approaches, the goal space is organized along several modules encoding for different types of goals: in that case goal sampling is hierarchical [START_REF] Forestier | Intrinsically motivated goal exploration processes with automatic curriculum learning[END_REF][START_REF] Colas | Curious: intrinsically motivated modular multi-goal reinforcement learning[END_REF]. Another interesting aspect of learning progress is that it allows avoiding catastrophic forgetting [START_REF] Colas | Curious: intrinsically motivated modular multi-goal reinforcement learning[END_REF] by identifying skills that were previously mastered, but are not mastered anymore. Using the learning progress to sample goals in image based environments is challenging due to the high dimensionality of the observations and as only been done in simple settings [START_REF] Laversanne-Finot | Curiosity driven exploration of learned disentangled goal spaces[END_REF].

In this paper we provide a solution to action-induced distractors in image-based environments for RL-based goal exploration algorithms. The general idea is to cluster the environment into different parts for which the agent can then estimate the associated learning progress. From a high level perspective our method can be viewed as a way to learn a prior over possible goals. The goal of this prior is to detect which parts of the environment should be explored. This prior can then be combined with any goal exploration algorithm in order to guide exploration. We experimentally demonstrate that our framework improves the overall performances of various novelty-based goal exploration algorithms and autonomously focuses the exploration on interesting learnable regions of the environment.

Contributions

The main contribution of this work is to study how learning progress can be estimated in image-based environments and how it can be used to guide the exploration of existing novelty-based goal exploration algorithms in Deep RL. The main contributions are:

• We introduce the GRIMGEP framework that allows estimating learning progress in imagebased environments for goal sampling in goal-conditioned Deep RL

• We combine it with existing goal exploration algorithms, and show that it improves the performances and avoids catastrophic forgetting of these algorithms in the presence of distractors.

• We provide an easily customizable 3D image-based environment constructed for studying noisy distractors in goal driven exploration tasks.

Related work

Goal exploration algorithms Intrinsically Motivated Goal Exploration Processes (IMGEPs) are a class of exploration algorithms that explore by repeatedly setting goals for themselves that they then try to achieve. They exist in two forms: population-based IMGEPs leveraging non-parametric models (e.g. [START_REF] Baranes | Active learning of inverse models with intrinsically motivated goal exploration in robots[END_REF], [START_REF] Forestier | Intrinsically motivated goal exploration processes with automatic curriculum learning[END_REF]), and, as we study in this paper, goal-conditioned RL IMGEPs leveraging Deep RL techniques (e.g. Nair et al. [2018], Pong et al.

[2020], [START_REF] Colas | Curious: intrinsically motivated modular multi-goal reinforcement learning[END_REF], [START_REF] Nasiriany | Disco rl: Distribution-conditioned reinforcement learning for general-purpose policies[END_REF]). Combined with an efficient curiosity mechanism for sampling goals, this approach has been shown to enable high-dimensional robots to learn very efficiently locomotion skills [START_REF] Baranes | Active learning of inverse models with intrinsically motivated goal exploration in robots[END_REF], manipulation of objects [START_REF] Rolf | Goal babbling permits direct learning of inverse kinematics[END_REF][START_REF] Sao | Socially guided intrinsic motivation for robot learning of motor skills[END_REF][START_REF] Ashvin V Nair | Visual reinforcement learning with imagined goals[END_REF][START_REF] Vitchyr | Skew-fit: State-covering self-supervised reinforcement learning[END_REF] or tool use [START_REF] Forestier | Intrinsically motivated goal exploration processes with automatic curriculum learning[END_REF]. Methods using absolute learning progress [START_REF] Baranes | Active learning of inverse models with intrinsically motivated goal exploration in robots[END_REF] to drive goal sampling were shown to scale up to real world environments with potentially many forms of distractors, including action-induced distractors [START_REF] Forestier | Intrinsically motivated goal exploration processes with automatic curriculum learning[END_REF][START_REF] Colas | Curious: intrinsically motivated modular multi-goal reinforcement learning[END_REF]. However, these works relied on abstract hand-defined goal and state spaces (e.g. based on object positions and velocities).

An exception is Laversanne-Finot et al. [2018], which used learning progress to sample goals in learned latent spaces, but this relied on population-based learners and required offline pre-training.

Here, we adapt the learning progress approach to Deep RL goal-exploring agents that perceive their environments through pixels, and study their properties in this context.

Recently some approaches studied goal exploration in the image-based goal-conditioned Deep RL framework. For instance, Nair et al.

[2018] learns a goal-conditioned policy on top of a learned embedding of the environment. [START_REF] Nasiriany | Disco rl: Distribution-conditioned reinforcement learning for general-purpose policies[END_REF] proposes a richer task representation in the form of a goal distribution. [START_REF] Warde-Farley | Unsupervised control through non-parametric discriminative rewards[END_REF] learns a goal conditioned policy by maximizing the mutual information between the goal state and the achieved state. While the goal policy achieves a good performance, the heuristic used for goal sampling is very simple and was not shown to scale to large environments with distractors. [START_REF] Vitchyr | Skew-fit: State-covering self-supervised reinforcement learning[END_REF] improves upon [START_REF] Ashvin V Nair | Visual reinforcement learning with imagined goals[END_REF] by designing a mechanism that incentivizes goal sampling to focus on exploring the frontier of the distribution of known goals, implementing a form of novelty search [START_REF] Lehman | Abandoning objectives: Evolution through the search for novelty alone[END_REF] (see also [START_REF] Eysenbach | Diversity is all you need: Learning skills without a reward function[END_REF] for a related approach). However, such methods can be naturally attracted by distractors that generate unpredictable or uncontrollable novel observations (either due to randomness, lack of cognitive capacity or lack of access to hidden information). Here, we experimentally identify this limit, and show that the learning progress based goal sampling mechanism we introduce can be used as a wrapper around these algorithms to enable them to avoid irrelevant or unlearnable regions of the environment, while keeping their efficient novelty-based exploration in learnable regions.

Curiosity-driven Deep RL A common trend to improve exploration in the classical reinforcement learning setting with sparse external rewards has been to supplement the task reward with intrinsic reward [START_REF] Bellemare | Unifying count-based exploration and intrinsic motivation[END_REF][START_REF] Pathak | Curiosity-driven exploration by self-supervised prediction[END_REF], Burda et al., 2018a[START_REF] Raileanu | RIDE: rewarding impact-driven exploration for procedurallygenerated environments[END_REF]. Approaches considering image-based low-level perception have used intrinsic rewards measuring various forms of novelty, based on counts [START_REF] Bellemare | Unifying count-based exploration and intrinsic motivation[END_REF] or prediction errors [START_REF] Pathak | Curiosity-driven exploration by self-supervised prediction[END_REF], Burda et al., 2018b]. While these approaches can give impressive results, they have fallen short on simple environments in the presence of a distractor that is partially controlled by one of the agent's actions. [START_REF] Kim | Active world model learning with progress curiosity[END_REF] studied the use of learning progress as an intrinsic reward to enable Deep RL agents to be robust to distractors. However, this approach relied on high-level disentangled state representations and did not include the notion of goal exploration.

Problem definition

Open-ended unsupervised exploration in image based environments

In the problem of open-ended unsupervised exploration the agent is tasked with acquiring a diversity of skills while having no access to the internal dynamics of the environment or expert knowledge. The only information the agent receives from the environment is in the form of observations that it collects. The agent must discover new states and learn how to reach them at the same time, often in the form of a goal policy.

This problem is particularly challenging when the agent has only access to observations in the form of images for two main reasons. First, since no measure of similarity between a state image and the goal image is given, the agent has to autonomously train the goal-conditioned reward function in an unsupervised manner. Second, since no information is given about which images represent feasible states in the environment, a mechanism must be constructed to autonomously create and sample plausible goals.

Intrinsically motivated goal exploration

Intrinsically motivated goal exploration processes (IMGEP) are a natural way to address this problem.

IMGEPs are a family of approaches that explore by repeatedly setting goals and then trying to achieve them. Those approaches usually consist of two parts: a goal reaching policy, and a goal sampling mechanism. Having a good policy ensures that goals are reached efficiently, and the goal sampling mechanism serves to discover new goals. The combination of these two elements ensures that the agent explores the environment efficiently. In this work we focus on the goal sampling mechanism.

In the case of image-based environments, goals cannot be sampled directly from the image space as such images are unlikely to represent images feasible in the environment. A simple solution, is to sample goals from the history of observations. However, in many cases, this history will not be filled uniformly from observations of the environment. For example, it will often be biased towards observations that are close to the starting position of the agent. Thus, sampling uniformly goals from that history will not provide very diverse goals. Therefore, in order to find interesting goals, current state-of-the-art methods rely on sampling methods based on the novelty of the visited states to maximize the diversity of sampled goals. A pseudocode of such an algorithm in shown by black lines in Algorithm 1.

However, those novelty-seeking approaches fall short in environments with distractors, i.e. where unpredictable (possibly structured) images can be perceived in certain situations (like a TV showing a flow of continuously novel images). Such distractors are problematic as a goal depicting a specific randomly occurring state is infeasible. Furthermore, as they produce a large amount of novelty, current approaches are attracted to them and oversample infeasible goals resulting in reduced performance and catastrophic forgetting.

IMGEPs using Absolute Learning Progress (ALP) instead of novelty have been shown to address both of those issues in simpler (not image-based) environments. In this paper, we present the GRIMGEP approach that extends the ALP based exploration to unsupervised image-based exploration. GRIMGEP uses ALP estimates to detect learnable regions of the goal space and then uses current novelty-based approaches only inside these regions. The extensions introduced by GRIMGEP are shown by blue lines in Algorithm 1.

Evaluation

In order to evaluate such goal exploration algorithm, the evaluation objective used by the experimenter consists in measuring the performance of the learned goal-reaching policy over a test set of goals which represents well the diversity of learnable skills in the chosen environment (note the distribution of goals in the test set is unknown to the learner):

argmax θ g∼G f (g, τ π θ ,g ) dg (1)
Where G is the goal space, τ π θ ,g the trajectory resulting from following the goal-conditioned policy π, parameterized by θ, and aiming for goal g. f is the goal fulfillment evaluation function specifying to what extent was the goal g reached in episode τ π θ ,g . This function is also defined by the experimenter for the purpose of evaluation, and the agent has no access to it during training.

In our case of image-based goal exploration, f is a measure of similarity between the goal and the final state i.e. to what extent are the same objects visible on the two images. Since it is unfeasible to integrate over the whole image space, a fixed set of testing goals (test set) is used for evaluation. For further details about the evaluation and the evaluation metric refer to sections 6.1, and appendix C.

Current novelty-based exploration approaches

In this work, we consider two representative novelty seeking goal exploration approaches: Skewfit [START_REF] Vitchyr | Skew-fit: State-covering self-supervised reinforcement learning[END_REF] and CountBased. Both Skewfit and Countbased train a generative model (β-VAE) online on the encountered states. They use the negative L2 distance between the goal and the state in this latent space as the reward for training the goal-conditioned policy. Both sample goals from a set of all the encountered states, however the way in which they prioritize goal sampling is different. This same prioritization method is also used for sampling the data to train the generative model.

Skewfit [START_REF] Vitchyr | Skew-fit: State-covering self-supervised reinforcement learning[END_REF] was shown to outperform a number of baselines for unsupervised goal exploration in environments without noisy distractors. It samples goals prioritized by how novel they appear to the current generative model (VAE) i.e. the probability of sampling a goal is inversely proportional to the probability of that goal under the generative model's distribution. For a more detailed description of Skewfit refer to Appendix A.1.

CountBased's goal sampling is inspired by how Go-Explore samples cells [START_REF] Ecoffet | First return, then explore[END_REF].

The states are downsampled to a 3x3 image, and quantized to 4 values (4 3 different colors). All the encountered states are counted as their quantized versions and the goal sampling biased accordingly, i.e. the probability of sampling a goal is inversely proportional to the number times its quantized version was encountered. For a more detailed description of CountBased refer to appendix A.2.

GRIMGEP

We present the Goal Regions guided Intrinsically Motivated Goal Exploration Process (GRIMGEP) framework. The GRIMGEP framework can be combined with any goal exploration algorithm, improving it by allowing the agent to identify and to sample goals from the interesting and learnable regions of the goal space.

GRIMGEP extends the ALP approach to unsupervised image-based exploration by fitting a clustering function on the storage of encountered states, applying this clustering function to a separate history of (attempted goal, achieved last state) pairs, estimating ALP for each of the clusters using this history, and sampling a cluster based on these ALP estimates. The wrapped IMGEP (which can be novelty based) is then used to sample a goal from the storage off all the encountered states, but constrained to those goals that the clustering function would assign to the sampled cluster. This pipeline is repeated for every sampled goal.

The motivation behind this scheme is that clusters corresponding to distracting unlearnable regions will have small ALP (performance doesn't change) and will be sampled rarely, while regions that are currently being learned or forgotten (performance rising or falling) will be sampled more. 2) The absolute learning progress (ALP) of each cluster is then computed by the ALP Estimation component using the history of all attempted goals and the corresponding outcomes. 3) In the Prior Construction component a cluster is sampled using the ALP estimates. The goal sampling prior is then constructed as the masking distribution assigning uniform probability over goals inside the sampled cluster (over all the goals in the history of encountered states that the clustering function would assign to this cluster) and 0 probability to goals outside the cluster. 4) A goal is then sampled, inside The Underlying IMGEP, from the distribution formed by combining the goal prior and the underlying IMGEP's novelty-based goal sampling distribution i.e. a novel looking goal is sampled from the sampled cluster. For further details see section 5 and the appendix B.

The GRIMGEP framework consists of the following four components as shown in Figure 1:

I. Clustering component.
It creates a clustering function that maps any goal to a cluster. Every time when a goal is to be sampled all the encountered states are passed through the backbone of a pretrained YOLO-V3 [START_REF] Redmon | Yolov3: An incremental improvement[END_REF] object detector (taken from [START_REF] Guo | Gluoncv and gluonnlp: Deep learning in computer vision and natural language processing[END_REF], Apache License), followed by a global average pooling layer [START_REF] Lin | Network in network[END_REF]. A PCA dimensionality reduction[F.R.S., 1901], and a GMM that follows it are both fit on those state embeddings. We leverage a pretrained object detector as it provides a latent space suitable for clustering different realistic objects, however this constrains our approach to realistic 3D environments. The clustering component is depicted in figure 5 of the appendix.

II. ALP Estimation component.

It assigns an absolute learning progress estimate to each cluster.

To estimate ALP a memory is kept containing all the goals that were sampled, the last state that was achieved while aiming for that goal, and the epoch in which this goal was sampled.

At the beginning of every epoch, the performances (rewards) corresponding to each of the previously sampled goals and the achieved last states are recomputed. We do so using the reward function from the underlying IMGEP.

By arranging the goals from the (goal, last_state) pairs history into clusters (using the clustering function fit in the clustering component) we construct a history of performances for each cluster. The performance estimates from the same epoch in the same cluster are averaged and only the last l epochs are kept in each history. The ALP is then estimated as the distance between the mean of the first and the last halves of the history as in the following equation:

ALP = | 1 l/2 l/2 i=0 h i -1 l/2 l i=l/2 h i | (2)
, where l is the length of the cluster's history, h is the cluster's history. This component is further explained in algorithm 3 in the appendix.

III. Prior Construction component. It creates the prior distribution. A cluster is sampled according to the per cluster ALP estimates. This sampling is achieved by a Multi-Armed Bandit with ALP estimates as utilities. The goal sampling prior is then defined as a masking distribution i.e. uniform over all goals from the encountered states history that the clustering function would assign to the sampled cluster, and zero elsewhere. This distribution is the means by which the underlying IMGEP is constrained to sample a goal only from the sampled cluster.

IV. Underlying IMGEP. It can be any goal exploration process that constructs a distribution for sampling goals. The underlying IMGEPs we study (Skewfit and CountBased) both sample goals from a set of encountered image states. When used by themselves (outside the GRIMGEP framework) this sampling is performed according to a novelty maximizing distribution p imgep . However, when used inside the GRIMGEP framework, the goal sampling distribution p imgep is multiplied with the prior distribution to define the final goal sampling distribution p as in the following equation:

p(g) = p prior (g)p imgep (g) g∈R p prior (g)p imgep (g) , ( 3
)
where R is the set of all encountered states. This multiplication has a masking effect on the novelty prioritizing distribution p imgep . Distribution p has non-zero probabilities only inside the cluster sampled by the Prior Construction component where they are proportional to those of p imgep .

Using p instead of p imgep for sampling a goal has the effect of using the underlying IMGEP only in the cluster sampled by the Prior Construction component.

Experiments

In this section, we present the experiments demonstrating the effectiveness of the GRIMGEP framework. We consider two novelty-seeking goal exploration processes: Skewfit [START_REF] Vitchyr | Skew-fit: State-covering self-supervised reinforcement learning[END_REF] and CountBased (inspired by GO-Explore [START_REF] Ecoffet | First return, then explore[END_REF]). We study their behavior when used alone and when wrapped inside the GRIMGEP framework on a 3D first-person environment containing a distractor. Additional experiments on a simpler 2D environment are explained in the appendix E.

1

More precisely, we aim to address the following questions:

• How do current approaches behave in the presence of action-induced distractors?

• How does the GRIMGEP framework change the behavior of current approaches in the presence of noisy distractors? • How important are the ALP estimates for the performance of the GRIMGEP framework?

Explore3D environment

Explore3D is a 3D first-person environment implemented using the MiniWorld [Chevalier-Boisvert, 2018] environment simulator (released under the Apache license). The environment consists of 14 different objects located in three separate rooms (TV room, Gallery room, Office room). The topology of the environment and the different rooms are depicted in Fig. 2. The agent's observation space includes RGB images depicting the first-person view of the environment as shown in figures 2b-2d. The agent can interact with the environment through 7 continuous actions: four for movement, two for rotating the camera and one for turning the TV on or off. The starting position is in the TV room facing the clock as depicted by the red triangle in Fig. 8a.

At the beginning of each episode the TV is in the OFF position (a black screen). When the agent turns on the TV a random image from the ImageNet [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF] validation set is sampled and shown on the screen. The turned on TV then changes the image every step with a probability of 0.1.

The environment provides a challenge for exploration in the following two aspects: 1) if no explicit exploration incentive is used the agent will not be able to explore farther from the starting TV room, and 2) novelty based exploration will be attracted to the diversity of images on the TV. Furthermore, the reaching of goals depicting a particular image on the TV is almost impossible as the agent has no control over which image, from the 50000 ImageNet validation images, will be shown on the screen.

For evaluation, we construct a test set consisting of 20 goal images depicting all the objects from the environment except the TV. The metric used for evaluation is the f1 score between the objects visible in the last state and in the goal. For further details about this metric, refer to appendix C.

Results

In the following section we present the experiments with Skewfit and CountBased. We refer to the GRIMGEP wrapped versions of those algorithms as GRIM-Skewfit and GRIM-CountBased. As Skewfit tends to behave unstably if exploring only with the novelty bonuses, the authors introduced a regularization hyperparameter α. α enables interpolation between a uniform (α = 0) and a completely exploratory distribution (α = -1) making Skewfit more robust.

Table 1 shows the performances at epoch 1000, and the percentage of all goals sampled up to epoch 1000 that depict the active TV. This is shown for CountBased, Skewfit with two different values α, and for those three approaches wrapped using the GRIMGEP framework. We can see that, for all three approaches, GRIMGEP improves the performance by reducing the sampling of TV goals.

In figure 3 the performance and the sampling of goals depicting the active TV are shown. More precisely, in figures 3c and 3d Skewfit with α = -0.25 is compared to GRIM-Skewfit with α = -0.75. This choice is made to compare the best (in terms of peak performance) Skewfit with the best GRIM-Skewfit. The remaining training curves are shown in the figure 6 of the appendix.

In figure 3 we can see that GRIM-CountBased and GRIM-Skewfit do not suffer from catastrophic forgetting and are able to achieve greater performances compared to Skewfit and CountBased. This is explained by observing that, unlike their unwrapped versions, neither GRIM-Skewfit nor GRIM-CountBased are drawn to the distractor.

These experiments show that wrapping the novelty seeking approaches with the GRIMGEP framework alleviates the problems introduced by noisy distractors.

How important are the ALP estimates for the performance of the GRIMGEP framework?

We study this question by doing an ablation study on the GRIMGEP wrapped CountBased. Figure 4 shows how GRIM-CountBased behaves when the clusters are sampled uniformly (GRIM-UNI-CountBased) as compared to when using learning progress (GRIM-LP-CountBased). The figure shows performance and percentages off all goals sampled for the distractor(active TV), the object visible from the starting position (clock), an object from the Gallery room (Impression), and an object from the office room (office chair), the figures for the rest of the objects are shown in the appendix 7. We can see that even when sampling clusters uniformly (GRIM-UNI-CountBased) the TV is not being oversampled. Since the agent's starting position is that of looking at the clock, at the beginning a lot of clusters consist mostly of goals depicting the clock. This results in both low sampling of the goals with the active TV, and high sampling of the goals with the clock. By looking at the sampling of goals with Impression and the office desk, we can see that GRIM-UNI-Countbased is not able to explore further away from the clock as opposed to GRIM-LP-CountBased.

These experiments show that, for CountBased, LP plays a significant role in the final performance of GRIMGEP by pushing the exploration farther from the starting position while at the same time not sampling the TV.

Conclusion

We studied the problem of unsupervised image-based goal exploration in first-person 3D environments with action-induced noisy distractors. Having demonstrated that current novelty-based approaches fail in such environments, we introduced the GRIMGEP framework. GRIMGEP can be wrapped around any novelty seeking goal exploration process enabling it to ignore the distractor. It detects the noisy distracting regions of the environment by estimating absolute learning progress, and uses the novelty seeking exploration only in the interesting regions. We have shown that wrapping state-of-the-art exploration algorithms inside the GRIMGEP framework allows them to ignore the distracting regions to ultimately improve their performance and alleviate catastrophic forgetting.

Limitations The limitation of this work is that the experiments are constrained to a simulated environment.

It remains an open question to see how the approach would behave in the real physical world with more sophisticated kinds of distractors.

Negative societal impact The autonomous exploration approach presented here focuses on discovering diverse behaviors and do not consider human preferences. That means that such algorithms when run in the real world would not have any bias for moral or socially acceptable behavior. It is, therefore, important that in the future such methods are coupled with mechanisms enabling alignment with human values and preferences. The code. Containing a new environment. Refer to the abstract for the URL. (d) Did you discuss whether and how consent was obtained from people whose data you're using/curating? [N/A] (e) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [No] As the only data we used was ImageNet [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF], we didn't find it nesseccary.

A Additional details on the studied underlying IMGEPs

A.1 Skewfit algorithm Skew-fit was proposed in [START_REF] Vitchyr | Skew-fit: State-covering self-supervised reinforcement learning[END_REF] primarily as an iterative algorithm that aims to train a generative model modeling the uniform distribution over the feasible goal space. To ensure this feasibility, Skew-fit constructs a dataset for training the generative model only from the already observed states. They can do this because they, like us, assume the goal space to be equivalent to the state space. This algorithm is then applied to unsupervised goal-driven exploration. The generative model is a VAE, and for its training dataset the replay buffer is used. They propose two mechanisms for goal sampling: sampling directly from the current VAE (which they call sampling from the generative model q G ψ ), and sampling from the skewed replay buffer (which they call sampling from the skewed distribution p skewed . In this work, though both mechanisms could be used with the GRIMGEP fremework, we focus on the latter as we found it worked better. This is aligned with the conclusion from [START_REF] Vitchyr | Skew-fit: State-covering self-supervised reinforcement learning[END_REF] as they use p skewed for more complex environments.

They train a goal-conditioned policy using SAC [Haarnoja et al., 2018b] coupled with HER [START_REF] Andrychowicz | Hindsight experience replay[END_REF] with the reward being the current VAE's negative L2 latent distance between the goal and the state.

In every epoch, goals are sampled, episodes run, new data added to the replay buffer, and both the VAE and the agent are trained. It should be noted that for training the SAC agent, data is sampled uniformly from the replay buffer, but for training the VAE and sampling of replacement goals for HER, the sample sampling mechanism as for goal sampling is used (a skewed replay buffer).

They denote this algorithm Skewfit + RIG though for simplicity we denote it Skewfit.

About the α hyperparameter It is relevant to note that Skewfit has a regularization hyperparameter α which interpolates between an uniform and a skewed distribution (-1 being completely skewed and 0 being uniform). This is done by exponentiating the original probability with α and then normalizing. In their work they experiment with 4 different values of this parameter (-0.25, -0.5, -0.75, -1.0) [START_REF] Vitchyr | Skew-fit: State-covering self-supervised reinforcement learning[END_REF].

In our experiments we found that -0.25 achieves better peak performance for Skewfit, but when wrapped inside the GRIMGEP frames -0.75 works better. As can be seen in figure 6.

That demonstrates another benefit of the GRIMGEP framework. Using Skewfit inside the GRIMGEP framework enables us to use more aggressive exploration bonuses (α = -0.75) without training becoming unstable.

A.2 CountBased approach

This approach is very similar to Skewfit, but it prioritizes goal sampling using a different intrinsic reward. This intrinsic reward is computed using a count-based novelty measure. As shown in equation 4, we downsize an image (to 3x3) and quantize each channel into 4 different values (4 3 different colors), then we count how many times the representation was observed. The goal's intrinsic reward is then computed by exponentiating this count by a hyperparmeter α ∈ [-1, 0], as in [START_REF] Vitchyr | Skew-fit: State-covering self-supervised reinforcement learning[END_REF]. As in Skewfit, we use this reward to construct the distribution to prioritize data for training the VAE, sampling exploration goals and sampling HER replacement goals. When using this approach as part of the GRIMGEP framework we multiply this distribution with the prior as in equation 3.

R CB (image) = count(quantize(downsize(image))) α (4)

B Details on the GRIMGEP approach

In this section we will present additional details of the GRIMGEP framework.

The biased distribution, p which was constructed by Prior construction component, is used everywhere where the underlying IMGEP would normally use p imgep . Therefore, for three purposes: 1) sampling goals for exploration, 2) sampling replacement goals for HER, and 3) since all of our underlying IMGEPS train their own VAE, biasing the training of this VAE.

A more detailed description of shown in the algorithm 2. First, we do a warmup phase to fill the replay buffer by doing random actions. Then, for the first start_exploration epochs, we do not use any kind of intrinsic rewards but just resample goals uniformly from the replay buffer. After this, we start using intrinsic rewards to prioritize the goal sampling. In practice, we make use of parallel computing by sampling ten goals per epoch. 5. The pretrained YOLO-v3 object detector whose backbone was used was taken from the GluonCV [START_REF] Guo | Gluoncv and gluonnlp: Deep learning in computer vision and natural language processing[END_REF] library ('yolo3_darknet53_coco'). This model was pretrained on the COCO [START_REF] Lin | Microsoft COCO: common objects in context[END_REF] object detection dataset (Creative Commons 4.0).

II. ALP Estimation component: This component's pseudocode is shown in the algorithm 3. The whole history of (goal, last_state) pairs is separated by the cluster in which the goal was assigned to. The performance estimates are recomputed as the reward function is the L2 1. Using per-cluster ALP estimates, sample a cluster according to the probability defined in the following equation:

p(c) = 4 5 LP T c C i=1 LP T c + 1 5 1 C (5)
,where C is the number of clusters and T is a hyperparameter. This is a Multi-Armed Bandit with a 20% uniform prior. The hyperparmaeter T enables us to set how much priority we want to give to the high ALP clusters, which is usefull if we have a lot of clusters with small LPs. 2. Construct the prior according to the following equation:

prior(g) = 1 nc , cl(g) == c 0, else (6)
,where c is the sampled cluster, n c is the number of goals from the replay buffer that would be assigned to cluster c, and cl is the clustering function assigning any state to a cluster. This is essentially a masking distribution giving uniform probabilities for goals inside the sampled cluster and zero for the rest. Another way of looking at this prior is saying what we sample a cluster and allow the underlying IMGEP to sample the goals only from this cluster.

B.1 Additional details about the experiments

In this section, we will show some additional figures explaining the results discussed in the main text 6.

Figure 6 shows the success rates and the percentage of all sampled goals that depict the active TV for CountBased, Skewfit with two different values of hyperparmaeter α, compared to those three approaches wrapped inside the GRIMGEP framework. We can see that GRIMGEP improves the performance for all three underlying IMGEPs and doesn't suffer from catastrophic forgetting. We can also see that, even though Skewfit with α = -0.75 eventually stops oversampling goals with the active TV, its performance is still diminished as a consequence of oversampling in the beginning. Shown for Countbased, Skewfit with two different values of α, compared to the three approaches wrapped inside the GRIMGEP framework. We can see that the performance of all the three approaches is improved, and catastrophic forgetting avoided, when wrapped inside the GRIMGEP framework.

The dots depict statistical significance (Welch t-test, p=0.05, 15 seeds)

Figure 7 shows the success rate and the goal sampling statistics for every object in the environment.

It compares the CountBased approach to its GRIMGEP wrapped counterpart. The first row shows objects from the TV room, the second row objects from the GALLERY room, and the third row objects from the OFFICE room. We can see that if ALP is not used, i.e. the clusters are sampled uniformly, the agent continues to sample many goals from the TV room, especially goals depicting the clock. This demonstrates the inability of the agent to explore further from the starting position (the starting position is in the TV room looking at the clock). Furthermore, we can see that the objects from the other rooms are sampled more if LP is used.

C Details on the evaluation metric

To evaluate the fulfillment of a goal, we construct the visible entities f1 score metric. Access to this metric is only available in evaluation. Using the environment simulator, a list of visible objects (entities) is extracted for both the goal and the last state from the episode. By treating objects visible in the goal as ground truth, precision and recall, and finally f1 score are calculated.

This metric, although it captures well the task of "seeing the goal image", is not perfect. We can notice that is insensitive both to the angle and the distance at which the agent is looking at the objects. Furthermore, in some cases an object may be visible on the state image with only a pixel of two. This difference, even though it's very hard for the agent to fix, will be treated as if the object is fully present. Nonetheless, we found the metric sufficiently accurate for comparing different algorithms.

D Hyperparametres

For the underlying IMGEPs, all the hyperparameters, including the ones for the training VAE, are the same as in [START_REF] Vitchyr | Skew-fit: State-covering self-supervised reinforcement learning[END_REF] in the "Visual Door" experiments, except the α hyperparameter (see Appendix A.1).

Other hyperparameters are shown in table 2.

E Preliminary experiments in the 2D PlaygroundRGB environment

The preliminary experiments used in constructing the GRIMGEP framework were done on a simpler 2D environment called PlaygroundRGB. Since some changes were made in the architecture for extending the agent to the 3D environment, the GRIMGEP architecture used in these experiments on PlaygroundRGB environment is not the same as the one described in the main paper.

In this section we will explain the environment those differences and the experiments. Clustering component Since the 2D environment does not have realistic objects and has a simpler underlying latent space an additional VAE was trained to create the clustering latent space on which to train the GMM i.e. the YOLO backbone, global average pooling, and the PCA were replaced with this VAE. This additional VAE is not the same VAE as the one used inside the underlying IMGEP to compute the rewards for training the policy. It is important that this VAE's latent space creates the features which are relevant for separating the regions and not on the specific details relevant for training the agent. For this purpose, we reduce the size of the VAE and its latent space. For the hyperparameters, see table 3. We train this VAE online after each epoch on the data uniformly sampled from the replay buffer.

In the experiments on the PlaygroundRGB environment ten different GMM [START_REF] Edward | The infinite gaussian mixture model[END_REF] models each having a different number of clusters (1, 3, . . . , 19) was trained. Then we the best one was chosen by their AIC [START_REF] Bozdogan | Model selection and akaike's information criterion (aic): The general theory and its analytical extensions[END_REF] score. This mechanism was inspired by the one used in [START_REF] Portelas | Teacher algorithms for curriculum learning of deep rl in continuously parameterized environments[END_REF]. The best GMM was set as the clustering function. Each epoch the process is repeated and a new GMM selected.

In the experiments on the Explore3D environment, the number of clusters was simply treated as the hyperparameter (k) as the procedure described above introduced unnecessary computation.

Figure 1 :

 1 Figure 1: Goal sampling procedure in the GRIMGEP framework. 1) First, the Clustering component clusters the goal space i.e. fits a clustering function.2) The absolute learning progress (ALP) of each cluster is then computed by the ALP Estimation component using the history of all attempted goals and the corresponding outcomes. 3) In the Prior Construction component a cluster is sampled using the ALP estimates. The goal sampling prior is then constructed as the masking distribution assigning uniform probability over goals inside the sampled cluster (over all the goals in the history of encountered states that the clustering function would assign to this cluster) and 0 probability to goals outside the cluster. 4) A goal is then sampled, inside The Underlying IMGEP, from the distribution formed by combining the goal prior and the underlying IMGEP's novelty-based goal sampling distribution i.e. a novel looking goal is sampled from the sampled cluster. For further details see section 5 and the appendix B.

Figure 2 :

 2 Figure 2: Different rooms of the 3D environment

Figure 3 :

 3 Figure 3: The performance and the percentage of sampled goals depicting the active TV (± std.err.). The dots depict statistical significance (Welch t-test, p=0.05, 15 seeds)

Figure 4 :

 4 Figure 4: The performance and the percentage of sampled goals depicting the various objects (± std. err.) The dots depict statistical significance (Welch t-test, p=0.05, 15 seeds)

Checklist 1 .

 1 For all authors... (a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [Yes] (b) Did you describe the limitations of your work? [Yes] refer to paragraph "Limitations' in the conclusion (c) Did you discuss any potential negative societal impacts of your work? -understand more autonomous learning agents if they are put in real world open environment [Yes] refer to paragraph "Negative societal impact" in the conclusion (d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes] 2. If you are including theoretical results... (a) Did you state the full set of assumptions of all theoretical results? [N/A] (b) Did you include complete proofs of all theoretical results? [N/A] 3. If you ran experiments... (a) Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)? [Yes] As an URL in the abstract (b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] Refer to the appendix D (c) Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)? [Yes] Refer to figure captions. (d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [No] 4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets... (a) If your work uses existing assets, did you cite the creators? [Yes] Pong et al. [2020], Chevalier-Boisvert [2018], Krizhevsky et al. [2012] and Guo et al. [2020] are all cited in the paper. (b) Did you mention the license of the assets? [Yes] Refer to the footnote on page 1 for Pong et al. [2020], and section 6.1 for Chevalier-Boisvert [2018], section 5 under "clustering component" for Guo et al. [2020] (c) Did you include any new assets either in the supplemental material or as a URL? [Yes]

Figure 5 :

 5 Figure5: The clustering component. A goal image is passed through the backbone of a pretrained yolo detector followed by a global average pooling operation. The dimensionality of this vector is then reduced using a PCA transformation and a GMM is used to assign the goal to a cluster.

Figure 6 :

 6 Figure6: The performance and the percentage of sampled goals depicting the active TV(± std. err.). Shown for Countbased, Skewfit with two different values of α, compared to the three approaches wrapped inside the GRIMGEP framework. We can see that the performance of all the three approaches is improved, and catastrophic forgetting avoided, when wrapped inside the GRIMGEP framework. The dots depict statistical significance (Welch t-test, p=0.05, 15 seeds)

  Figure 7: The Performance and the sampling distribution of all the objects in the environment. Shown for GRIMGEP wrapped CountBased, as compared to GRIMGEP wrapped CountBased which doesn't use absolute learning progress(ALP), i.e. samples the cluster uniformly. We can see that ALP enables the agent to move from the starting TV room (figures b-e) to the other rooms of the environment(f-n), thereby reaching superior performances (figure a). The dots depict statistical significance (Welch t-test, p=0.05, 15 seeds).

E. 1

 1 Changes to the GRIMGEP between the experiments in the PlaygorunrRGB(2D environment) an the Explore3D

  Figure 10: Comparison of the GRIMGEP framework with the sampling of clusters according to ALP and uniformly on the preliminary 2D PlaygroundRGB environment. It is visible that GRIMGEP works better when LP is used to select the most interesting cluster. The dots depict statistically significant (p < 0.05, Welch's t-test) results when compared to the LP version (GRIM-LP-imgep_name). The shaded areas correspond to standard errors and the bold line to the mean (smoothed over 25 epochs).

Figure 11 :

 11 Figure 11: Comparison of Online-RIG alone and in combination with the GRIMGEP framework.On the preliminary 2D PlaygroundRGB environment. We can see that GRIMGEP improves the performance of OnlineRIG by focusing on the Object room. Ten seeds were used, and the dots depict statistically significant (p < 0.05, Welch's t-test) results. The shaded areas correspond to standard errors and the bold line to the mean (smoothed over 25 epochs).

  After each epoch of the generative model training, a new Skewed dataset is constructed. It is constructed by setting each state's weight (probability of being sampled) as inversely proportional to the current generative model's probability of that state. This skewed dataset is then used to train the generative model in the next epoch. This dataset's entropy increases in each epoch, so the model converges to modeling a uniform distribution.

  For each of those goals, we sample a new cluster, construct a new prior and then sample a new goal using the probability distribution obtained by combining the prior and the goal sampling probability distribution of the IMGEP. This is depicted by the algorithm 2 where the blue colors depict GRIMGEP addition.

	Algorithm 2: Intrinsically Motivated Goal Exploration Process in the Deep RL setting (blue lines
	correspond to GRIMGEP additions)
	π.init() R.init()	// goal-conditioned policy and the reward function (VAE)
	encoutered_states = []
	performance_history = []	// A history of (sampled_goal, last_state) pairs
	// random acting rollouts to fill the replay buffer
	for N_warmup do
	traj = env.random_rollout()
	encountered_states.add_states(traj)
	// exploration
	for ep in N_epochs do
	for N_parallel_processes do
	if ep >start_exploration then
	// novelty-based exploration of the underlying IMGEP
	goal_intrinsic_rewards = compute_intrinsic_rewards(encountered_states)
	// Assign a cluster_id to every goal from the perfromance_history
	cluster_ids = clustering_fn(performance_history)
	// Recompute the performances and estimate per-cluster ALP
	per_cluster_alp = estimate_ALPs(cluster_ids, performance_history, R)	// by algo. 3
	sampled_cluster=sample(clusters, per_cluster_alp)	// by eq. 5
	// Create a masking distribution by eq. 6
	prior=construct_prior(sampled_cluster, encountered_states)
	// Integrate the bias by eq. 3
	goal_intrinsic_rewards *= prior
	else
	goal_intrinsic_rewards = [1,1. . . 1]	// at the beginning sample uniformly
	goal = prioritized_sample_goal(encountered_states, goal_intrinsic_rewards)
	trajectory = env.policy_rollout(π, goal)
	encountered_states.add_states(trajectory)
	performance_history.add_pair((goal, trajectory.last_state))
	R.train() π.train()		// fit reward function (VAE) and policy
	clustering_fn.train(encountered_states)	// fit PCA and GMM
	I. Clustering component: The pipeline of this component is shown in figure

  The clustering VAE hyperparameterss are shown in table3.

	representation size	3
	batch_size	128
	beta	1
	lr	0.001
	Encoder
	kernel sizes	[5, 3]
	num of channels [4, 4]
	strides	[3, 2]
	Decoder
	kernel sizes	[3, 3]
	num of channels [4, 4]
	strides	[2, 2]

Table 3 :

 3 The Clustering VAE hyperparmaeters. A reduced version of the training VAE used in the underlying IMGEPs.

We reused the follwing codebases in our experiments: https://github.com/vitchyr/rlkit (MIT licence), GluonCV[Guo et al.,

2020] (Apache License), sklearn[START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] (BSD licence)

E.2 PlaygroundRGB

The PlaygroundRGB environment consists of three rooms through which the agent can move. The agent's actions space consists of 3 continuous actions which control the position of the gripper and can close/open /the arm. It observes the environment as an image of a top down view of the current room.

The available rooms are depicted in Fig. 8 and the topology of the environment in Fig. 8a. The agent always starts in the Start room. All the possible goals inside this room are very easy and require only moving the gripper to the correct location. The Object room represents the interesting part of the environment as it contains a movable object and is the only non-distracting part of the environment. The TV room plays the role of an action induced noisy distractor. This room contains a TV that can be turned on by closing the gripper. When the TV is turned ON the location of the TV and the background color are randomized (a random color from a set of 5).

In short, to solve the task, exploration should guide the agent out of the Start room, away from the TV room, and into the Object room.

Since goals from the Start room are easy, we expect any goal exploration algorithm to learn goals inside this room. However, only algorithms exploring well should master goals inside the Object room. This is why, for evaluation, we construct a static test set of 25 goals from the Object room. Goals are completed if in the last state both objects are in the correct location. The performance of the agent is the average success over this evaluation set.

For evaluation, the 25 test goals form the OBJECT room were constructed by selecting 5 possible locations (center, NW, NE, SW, SE) and doing the Cartesian product of these locations for both objects. We evaluate the location of an object as correct if its L-∞ distance from the goal's location is less then 0.2.

E.2.1 Experiments

In these experiments we address the same questions as in the Explore3D. Those questions are:

• How do current approaches behave in the presence of action-induced distractors?

• How does the GRIMGEP framework change the behavior of current approaches in the presence of noisy distractors? • How important are the ALP estimates for the performance of the GRIMGEP framework?

We study this questions using two novelty seeking IMGEPs (Skewfit and CountBased) and one which doesn't have exploration incentives and samples goals uniformly from the history of encountered states (OnlineRIG).

How do current approaches behave in the presence of noisy distractors?

We can see, in Fig. 9, that Skewfit is heavily drawn to the noisy part of the TV room resulting in very low sampling of other parts of the environment, notably the Object room. As a result of not exploring the Object room enough, the final performance diminishes.

CountBased also samples a lot of goals in the TV room. However, in comparison to Skewfit, the focus is separated between both the TV-on and TV-off goals (see Fig. 9d and Fig. 9e).

OnlineRIG is also not able to achieve high performances. The reason is that, due to the lack of exploration incentive, it samples goals mostly from the Start room as can be seen in Figure 11).

Overall, this experiment demonstrates that this environment requires exploration incentives but, that these incentives should not be novelty based.

How does the GRIMGEP framework change the behavior of current approaches in the presence of noisy distractors?

To answer this question we wrap the Skewfit and CountBased with the GRIMGEP framework. As can be seen in Fig. 9, both GRIM-Skewfit and GRIM-CountBased, focus more on the Object room and less on the TV room than their unwrapped counterparts. This results in much better performances at the end of the training. When using OnlineRIG inside the GRIMGEP framework we can again observe a strong focus on the Object room (see Fig. 11b), but we can also see that this focus alone is not sufficient to greatly improve performance. This leads us to the conclusion that GRIMGEP successfully does two things: 1) it detects the relevant part of the goal space, and 2) successfully uses the novelty seeking exploration in this relevant region.

How important are the ALP estimates for the performance of the GRIMGEP framework?

We study this question by doing an ablation study on the cluster sampling technique. In the experiments in Figure 10a we test how cluster sampling based on LP (GRIM-LP-imgep_name) differs in comparison to uniform cluster sampling (GRIM-UNI-imgep_name). As shown in Fig 10a, GRIMGEPs that use LP outperform the GRIMGEPs that don't. Furthermore, we can see that, when sampling the clusters uniformly, the Object room is sampled considerably less. We can see that GRIMGEP improves the performance of both Skewfit, CountBased and OnlineRIG. Furthermore, we can see that the performance of GRIM-OnlineRIG is inferior to that of GRIM-Skewfit and GRIM-CountBased. Since OnlineRIG doesn't use exploration bonuses (samples goals uniformly), this shows that GRIMGEP is able to both detect the interesting region of the goal space and reuse the novelty-based exploration of CountBased and Skewfit in this region. Ten seeds were used, and the dots depict statistically significant (p < 0.05, Welch's t-test) results. The shaded areas correspond to standard errors and the bold line to the mean (smoothed over 25 epochs).


Denise Aregba-Driollet

Afaf Bouharguane
email: afaf.bouharguane@math.u-bordeaux.fr

S T Éphane Brull
email: stephane.brull@math.u-bordeaux.fr

A KINETIC DISCONTINUOUS GALERKIN METHOD FOR THE NONCONSERVATIVE BITEMPERATURE EULER MODEL

Keywords: nonconservative hyperbolic system, discrete BGK approximation, discontinuous Galerkin methods, Runge-Kutta methods AMS subject classification 65M08, 35L60, 76X05, 35Q31

This paper is devoted to the construction of a discontinuous Galerkin discretisation for the nonconservative bitemperature Euler system via a a discrete BGK formulation. This formulation is compatible with the entropy properties of the system and thus provides admissible solutions. The DG method is used to approximate the linear transport part of the BGK model while the force and source-terms are treated implicitly but with explicit expressions. High order in time has also been investigated using SSP Runge-Kutta methods. We numerically show the good agreement of our results with the ones provided by other schemes, including solutions with shocks.

Introduction

This paper is devoted to the approximation of the bitemperature Euler system from a Discontinuous Galerkin (DG hereafter) method applied to a kinetic formulation proposed by Aregba and Natalini [START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF] and that is based on the discrete BGK model. The present fluid model is able to treat out of equilibrium regimes. The standard strategy in plasma physics, for simulating such regimes is to develop PIC methods, that directly solve the kinetic equations. But these methods are computationally very expensive. Therefore, the bitemperature model is a compromise between the precision of the kinetic models and the lower numerical cost of fluid models. This system is made of two conservation equations for mass and momentum and of two non conservative equations for ionic and electronic energies. It describes the interaction of a mixture of thermalized ions and electrons in a quasi-neutral regime. This system is nonconservative because of the presence of a relaxation source term and of products between velocity and pressure gradients. Those products make difficult the definition of weak solutions. In [START_REF] Maso | Definition and weak stability of nonconservative products[END_REF], a general framework has been developed in order to define shocks in this context by using families of paths. The generalization of this approach to a numerical setting has been considered in [START_REF] Parés | Path-conservative numerical methods for nonconservative hyperbolic systems[END_REF]. However, even if the path is known theoretically, the numerical determination of the path is delicate [START_REF] Abgrall | A comment on the computation of non-conservative products[END_REF]. In [START_REF] Coquel | Numerical methods for weakly ionized gas[END_REF], the authors consider the bitemperature Euler system with diffusive terms and by assuming that the electrons are isentropic. In that case, the system is transformed into a system of conservation laws. This approach has also been used in [START_REF] Sangam | Derivation and numerical approximation of two-temperature Euler plasma model[END_REF] where magnetic fields are considered and a conservation equation on electronic entropy is derived. In [START_REF] Brull | Lhébrard Modelling and entropy satisfying relaxation scheme for the nonconservative bitemperature Euler system with transverse magnetic field[END_REF], magnetic fields are also considered in a transverse magnetic configuration. Moreover, a Suliciu scheme is derived and proved to be entropic. In [START_REF] Wargnier | Numerical treatment of the nonconservative product in a multiscale fluid model for plasmas in thermal nonequilibrium: application to solar physics[END_REF], the authors perform a Chapman-Enskog expansion by introducing a small parameter representing the ratio between electronic and ionic molecular masses. At the end, they obtain a system with an hyperbolic part for ions and a parabolic part for electrons.

Discrete BGK schemes have been introduced in [START_REF] Natalini | A discrete kinetic approximation of entropy solutions to multidimensional scalar conservation laws[END_REF][START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF]. The method is performed by a transport projection approach. The advantage is to put the nonlinearity inside the relaxation term whereas the transport term is linear with constant advection velocities.

More precisely, the bitemperature model under consideration is obtained in [START_REF] Aregba-Driollet | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF], by using an hydrodynamic limit starting from a BGK model coupled with Poisson and Ampère equations in a quasi-neutral regime. In particular, it is shown that the nonconservative terms come from the Ohm's law defining the electric field. By using this asymptotic a finite volume BGK scheme has been designed in [START_REF] Aregba-Driollet | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF] by generalizing the Aregba-Natalini method [START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF] to a nonconservative setting. In particular, a force term is incorporated in the discrete BGK formulation in order to deal with the nonconservative terms. This formulation leads to the resolution of an hyberbolic system from a BGK relaxation process. Hence we have to deal with a linear advection term and a source term. Moreover, a Suliciu approach is developed in [START_REF] Aregba-Driollet | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF] and several comparisons are performed with different schemes. Next, in [START_REF] Aregba-Driollet | Modelling and numerical study of the polyatomic bitemperature euler system[END_REF] the kinetic approach has been generalized to a polyatomic setting by using a kinetic model with a continuous energy variable. The bi-dimensional case has been considered in [START_REF] Aregba-Driollet | A discrete velocity numerical scheme for the two-dimensional bitemperature Euler system[END_REF] and a second order finite volume scheme has been obtained. In [START_REF] Aregba-Driollet | About viscous approximations of the bitemperature euler system[END_REF], a Navier-Stokes system has been derived from a Chapman-Enskog expansion and by computing viscous terms, generalizing the model considered in [START_REF] Chalons | Coquel Navier-Stokes equations with several independant pressure laws and explicit predictorcorrector schemes Numerische Math[END_REF]. However, as far as we know, there is no work on the implementation of high order methods for non-conservative bitemperature model. For similar works, recently, in [START_REF] Abgrall | High order preserving deferred correction implicit-explicit schemes for kinetic models[END_REF][START_REF] Abgrall | Some preliminary results on high order asymptotic preserving computationally explicit kinetic scheme[END_REF], the authors developed a discrete BGK formulation for the compressible Euler system and next perform a RD scheme in the space variable. The time discretization is obtained by using a DeC method [START_REF] Abgrall | High order schemes for hyperbolic problems using globally continous approximation and avoiding mass matrices[END_REF]. In [START_REF] Coulette | High-order implicit palindromic discontinuous Galerkin method for kinetic-relaxation approximation[END_REF] the authors constructed a high-order implicit palindromic discontinuous Galerkin method from kinetic-relaxation approximation for solving general hyperbolic systems of conservation laws. In [START_REF] Gerhard | Michel-Dansac Unconditionally stable and parallel Discontinuous Galerkin solver[END_REF], this approach is used to solve Maxwell's equations. The aim of the present work is to provide high order methods for the bitemperature Euler system which use a DG scheme on the discrete BGK formulation with k-th degree basis for the spacial discretization and with an order (k + 1)-SSP-Runge-Kutta method for the time integration. This paper is organised as follows. Section 2 introduces the bitemperature Euler system. Section 3 deals with the presentation of the discrete BGK models in the conservative and in the nonconservative case. In particular the obtention of the bitemperature model from the kinetic formulation developed in [START_REF] Aregba-Driollet | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF] is explained. In section 4, the DG method applied to the discrete BGK approximation is described. In particular, the order in time is increased in order to be consistent with the order in the space variable. Section 5 is devoted to numerical experiments. We investigate the accuracy of the proposed schemes on some examples. We also compare our methods with other existing numerical schemes. An important result is that they converge to the same solutions, even in the presence of shocks.

Bitemperature Euler system

Subscripts e and i respectively denote electronic and ionic quantities. We denote by ρ e and ρ i the electronic and ionic densities, ρ = ρ e + ρ i the total density, m e and m i the related masses, c e and c i the mass fractions. These variables satisfy [START_REF] Abgrall | High order schemes for hyperbolic problems using globally continous approximation and avoiding mass matrices[END_REF] ρ e = m e n e = c e ρ,

ρ i = m i n i = c i ρ, m e > 0, m i > 0, c e + c i = 1.
Quasineutrality is assumed, so that the ionization ratio Z = n e /n i is a constant. This implies that the electronic and ionic mass fractions are constant and given by ( 2)

c e = Zm e m i + Zm e , c i = m i m i + Zm e .
Electronic and ionic velocities u e , u i are assumed to be in equilibrium in the model. Hence, u e = u i = u, where u denotes mixture velocity. The pressure of each species satisfies a gammalaw with its own γ exponent :

(3) p e = (γ e -1)ρ e ε e = n e k B T e , p i = (γ i -1)

ρ i ε i = n i k B T i , γ e > 1, γ i > 1,
where k B is the Boltzmann constant (k B > 0), ε α and T α represent respectively the internal specific energy and the temperature of species α for α = e, i.

Denoting by | • | the euclidean norm in R D , the total energies for the particles are defined by ( 4)

E α = ρ α ε α + 1 2 ρ α |u| 2 , α = e, i.
We denote by ν ei ≥ 0 the interaction coefficient between the electronic and ionic temperatures. The model consists of two conservative equations for mass and momentum and two non-conservative equations for the energies:

(5)

         ∂ t ρ + div(ρu) = 0, ∂ t (ρu) + div(ρu ⊗ u + (p e + p i )I) = 0, ∂ t E e + div(u(E e + p e )) -u • ∇ (c i p e -c e p i ) = ν ei (T i -T e ), ∂ t E i + div(u(E i + p i )) + u • ∇ (c i p e -c e p i ) = -ν ei (T i -T e ),
where I represents the identity matrix in R 3 . In the following we denote

(6) U = (ρ, ρu, E e , E i ), U α = (c α ρ, c α ρu, E α ).
The system (5) is hyperbolic, diagonalisable and owns 3 eigenvalues λ -, λ 0 (with multiplicity D + 1 where D is the space dimension), λ + : for any ω ∈ S D-1

λ -= u • ω -a, λ 0 = u • ω, λ + = u • ω + a where (7) a = α=e,i γ α p α ρ
is the sound velocity. The fields related to λ ± are genuinely nonlinear, while the field related to λ 0 is linearly degenerate. Defining the total energy E = E e + E i and the total pressure p = p e + p i , one can note that if U is a solution of system (5) then (ρ, ρu, E) satisfies the following conservative system:

(8)      ∂ t ρ + div(ρu) = 0, ∂ t (ρu) + div(ρu ⊗ u + pI) = 0, ∂ t E + div(u(E + p)) = 0.
If γ e = γ i this is the wellknown monotemperature Euler system. But even in this case, one has to deal with one more equation to determine electronic and ionic temperatures. If γ e ̸ = γ i system (8) is not closed. We want to underline the fact that in both cases, the solutions of system [START_REF] Aregba-Driollet | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF] are to be defined in the context of non-conservative equations were the product of a possibly discontinuous function with a Dirac measure appears. To give a sense to such solutions, one has to bring more physical information. In [START_REF] Aregba-Driollet | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF] we obtained solutions of (5) as hydrodynamic limits of solutions of an underlying, physically realistic BGK model. The entropy-entropy flux of species α being defined as

(9) η α (U α ) = - ρ α m α (γ α -1) ln (γ α -1)ρ α ε α (ρ α ) γα + C , Q α (U α ) = η α (U α )u,
the total entropy-entropy flux pair for ( 5) is ( 10)

η(U) = η e (U e ) + η i (U i ), Q(U) = η(U)u
and we proved the following entropy inequality for these hydrodynamic limits:

(11) ∂ t η(U) + divQ(U) ≤ - ν ei k B T i T e (T i -T e ) 2 .
We then defined an admissible solution of (5) as a solution satisfying this inequality.

In the following we consider the 1D version of system (5):

(12)            ∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρu) + ∂ x (ρu 2 + p e + p i ) = 0, ∂ t E e + ∂ x (u(E e + p e )) -u∂ x (c i p e -c e p i ) = ν ei (T i -T e ), ∂ t E i + ∂ x (u(E i + p i )) + u∂ x (c i p e -c e p i ) = -ν ei (T i -T e ).

BGK models

This section is devoted to the presentation of discrete BGK models that have been introduced for system of conservation laws in [START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF] and then generalized to the non-conservative case in [START_REF] Aregba-Driollet | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF], [START_REF] Aregba-Driollet | Nonconservative hyperbolic systems in fluid mechanics[END_REF].

3.1. Underlying kinetic (BGK) models for the conservative compressible Euler system. We start from BGK models for the Euler monotemperature equations. Denoting

U = (ρ, ρu, E) ∈ Ω ⊂ R 3 , F (U ) = (ρu, ρu 2 + p, u(E + p)),
the Euler system is a system of conservation laws:

(13) ∂ t U + ∂ x F (U ) = 0.
We assume that p = (γ -1)(E -1 2 ρu 2 ). We follow the framework proposed by F. Bouchut in [START_REF] Bouchut | Construction of BGK Models with a Family of Kinetic Entropies for a Given System of Conservation Laws[END_REF]. We define a measure space (X, dξ), a real valued function a defined on X, a "maxwellian function" M from R 3 × X onto R k , and a "moment operator" P from X to L(R k , R 3 ) such that for all U ∈ Ω:

X P (ξ)(M (U, ξ))dξ = U, X P (ξ) (a(ξ)M (U, ξ)) dξ = F (U ). Let f ε (x, t, ξ) ∈ R k be a solution of ∂ t f ε + a(ξ)∂ x f ε = 1 ε (M (U (f ε ), ξ) -f ε ), with U (f ε )(x, t) = X P (ξ)(f ε (x, t, ξ))dξ. Formally if lim ε→0 f ε = f , then f (x, t, ξ) = M (U (f )(x, t), ξ
) and U (f ) is a solution of [START_REF] Brull | Lhébrard Modelling and entropy satisfying relaxation scheme for the nonconservative bitemperature Euler system with transverse magnetic field[END_REF].

In [START_REF] Bouchut | Construction of BGK Models with a Family of Kinetic Entropies for a Given System of Conservation Laws[END_REF], conditions are given for the existence of microscopic entropies compatible with all the entropies of the macroscopic limit.

3.1.1. Example 1: a physically realistic BGK model. Here we set

X = R 3 , ξ = v, dξ = dv, a(ξ) = v 1 and M (U, v) ∈ R is given by (14) M (U ) = n (2πk B T /m) 3/2 exp(- |v -u| 2 2k B T /m
).

The moment operator is defined as

P (ξ)(M ) = (m, mv 1 , m|v 2 | 2 )M. As f ε (x, t, v) ∈ R,
X = {1, 2}, a(ξ) = λ ξ with λ 2 > λ 1 , P (ξ) = I d , k = 3 and M (U, 1) = λ 2 U -F (U ) λ 2 -λ 1 , M (U, 2) = -λ 1 U + F (U ) λ 2 -λ 1 ,
see [START_REF] Aregba-Driollet | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF]. For any Euler entropy, the existence of related microscopic entropies is ensured under Liu's subcharacteristic condition, see [START_REF] Bouchut | Construction of BGK Models with a Family of Kinetic Entropies for a Given System of Conservation Laws[END_REF]:

σ(F ′ (U )) ⊂]λ 1 , λ 2 [.

3.2.

BGK models for the nonconservative bitemperature Euler equations. We take a BGK model for the monotemperature Euler system (13) with γ = γ e and γ = γ i . We choose

X e = X i = X, a e (ξ) = a i (ξ) = a(ξ), P e (ξ) = P i (ξ) = P (ξ).
In order to approximate the nonconservative products we define a linear operator N such that

X P (ξ)(N M α (U α , ξ))d α ξ = -(0, ρ α , ρ α u α ).
In the case of example 1:

N f = v 1 ∂ v1 f. In the case of example 2: N f (ξ) =   0 0 0 -1 0 0 0 -1 0     f 1 (ξ) f 2 (ξ) f 3 (ξ)   .
The equations for f ε e and f ε i are coupled with the ones for the electric field E:

(15)                          ∂ t f ε e + a(ξ)∂ x f ε e + q e m e E ε N f ε e = 1 ε (M e -f ε e ) + B ei (f ε e , f ε i ), ∂ t f ε i + a(ξ)∂ x f ε i + q i m i E ε N f ε i = 1 ε (M i -f ε i ) + B ie (f ε e , f ε i ), ∂ t E ε = - 1 ε 2 q e m e ρ ε e u ε e + q i m i ρ ε i u ε i , ∂ x E ε = 1 ε 2 q e m e ρ ε e + q i m i ρ ε i .
The source-term B αβ is such that

X P (ξ)B αβ dξ = (0, 0, ν ei (T β -T α )).
When ε tends to 0, we have formally

u e = u i = u, q e m e ρ e + q i m i ρ i = 0, M α (U α ) = f α .
Quasineutrality holds: ρ = ρc e = ρc i and c e , c i are the constants defined in [START_REF] Abgrall | A comment on the computation of non-conservative products[END_REF]. By taking the moments, it comes that ( 16)

                 ∂ t ρ α + ∂ x (ρ α u) = 0, α = e, i, ∂ t (ρ α u) + ∂ x (ρ α u 2 + p α ) - q α m α Eρ α = 0, α = e, i, ∂ t E e + ∂ x (u(E e + p e )) - q e m e Eρ e u = ν ei (T i -T e ), ∂ t E i + ∂ x (u(E i + p i )) - q i m i Eρ i u = -ν ei (T i -T e ). Considering      ∂ t (ρc e u) + ∂ x (ρc e u 2 + p e ) - ρ e q e m e E = 0, ∂ t (ρc i u) + ∂ x (ρc i u 2 + p i ) - ρ i q i m i E = 0 leads to the expression of E ρ e q e m e E = - ρ i q i m i E = c i ∂ x p e -c e ∂ x p i and ∂ t (ρu) + ∂ x (ρu 2 + p e + p i ) = 0. Hence U = (ρ, ρu, E e , E i
) is a solution of system (5). Theorem 3.1. Suppose that there exists microscopic entropies for the kinetic model ( 15) related to the entropy η. Let U be a solution of the Euler bitemperature model ( 5) obtained by passing to the limit in [START_REF] Collela | The piecewse parabolic method (PPM) for gas dynamical simulations[END_REF]. Then one has the entropy inequality

∂ t η(U) + ∂ x Q(U) ≤ - ν ei k B T i T e (T i -T e ) 2 .
We define such a solution U as an admissible solution.

In the case of examples 1 and 2 above, the microscopic entropies exist, see [START_REF] Aregba-Driollet | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF] for details.

A Discontinuous Galerkin (DG) scheme

This section is devoted to discretisation of the discrete BGK model developped in section 3. We denote ∆x and ∆t the space and time steps and we mesh the real line by cells

C K = [x K-1 2 , x K+ 1 2 ] with x K+ 1 2 -x K-1 2 = ∆x.
In practice an interval is considered with appropriate boundary conditions : K ∈ {1, . . . , N }. We also denote by P k the space of all polynomials of degree at most k.

4.1. Preliminary. Let us consider a transport equation ( 17)

∂ t f + v∂ x f = 0
where v ∈ R and f (t, x, v) ∈ R. We look for an approximation of f (t, •, v) under the form

N K=1 f K (t, x, v) where each f K (t, •, v) has C K as support.
Let {Φ K j , j = 0, . . . , k -1} be a basis of polynomial functions defined on C K . Multiplying equation ( 17) by Φ K i and integrating over

C K yields C K ∂ t f K (t, x, v)Φ K i (x)dx-v C K f K (t, x, v)(Φ K i ) ′ (x)dx+v fK+ 1 2 Φ K i (x K+ 1 2 )-v fK-1 2 Φ K i (x K-1 2 ) = 0.
We have to choose a flux v fK+ 1 2 = h(f (x -

K+ 1 2 ), f (x + K+ 1 2
), v). We can set the following

(18) h(f, g, v) = v λ 2 λ 2 -λ 1 f - λ 1 λ 2 -λ 1 g + λ 1 λ 2 λ 2 -λ 1 (g -f )
where

λ 1 ≤ 0 ≤ λ 2 .
Another choice is the upwind flux:

denoting v + = max(v, 0), v -= max(-v, 0) (19) h(f, g, v) = v + f -v -g. Now we approximate f K (t, x, v) as (20) f K (t, x, v) = k-1 j=0 f K j (t, v)Φ K j (x).
We obtain

k-1 j=0 M K ij ∂ t f K j (t, v) -v k-1 j=0 S K ij f K j (t, v) + v fK+ 1 2 Φ K i (x K+ 1 2 ) -v fK-1 2 Φ K i (x K-1 2 ) = 0
where we have defined two matrices M K and S K :

M K ij = C K Φ K j Φ K i dx, S K ij = C K Φ K j Φ K i ′ dx. Denoting f K = (f K 0 , . . . , f K k-1 )
we have an ordinary differential system (21)

M K ∂ t f K (t, v) + S K f K (t, v) = F (t, v)
We denote f n+1 (x, v) = X ∆t f n (x, v) the obtained numerical scheme.

It is well-known that DG schemes may oscillate when sharp discontinuities are present in the solution. Hence in order to control these instabilities we consider a generalized slope limiters [START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications[END_REF]. We first define the interface fluxes as

v - j+ 1 2 = f K j -m f K j -f K j (x j+ 1 2 ), f K j -f K j-1 , f K j+1 -f K j v + j-1 2 = f K j + m f K j -f K j (x j-1 2 ), f K j -f K j-1 , f K j+1 -f K j
where f K j is the average of f K on C K and where m is the minmod function limiter

m(a 1 , a 2 , a 3 ) = s • min j |a j | if s = sign(a 1 ) = sign(a 2 ) = sign(a 3 ) 0 otherwise
Then the generalized slope limiter technique consists in replacing f K j on each cell C K with ΛΠ h defined by

ΛΠ h (f K j ) = f K j if v + j-1 2 = f K j (x j-1 2 ) and v + j+ 1 2 = f K j (x j+ 1 2 ) f K j + (x-xj ) ∆x/2 m(v (1) j , f K j+1 -f K j , f K j -f K j-1 ) otherwise. Remark 1. Define U (t, x) = X P (ξ)(f (t, x, ξ))dξ. By linearity if f (t, x, v) = N K=1 k-1 j=0 f K j (t, v)Φ K j (x) then U (t, x) = N K=1 k-1 j=0 U K j (t)Φ K j (x) with U K j (t) = X P (ξ)(f K j (t, ξ))dξ.

4.2.

A DG scheme for the bitemperature Euler system. Suppose that n ≥ 0 being fixed, we have U n = (ρ, ρu, E e , E i ). We define U n α = (ρ α , ρ α u, E α ) with ρ α = ρc α .

Step 1: projection onto equilibrium We define

f n α = M α (U n α ), α = e, i.
Step 2: transport by DG method. For α = e, i f

n+ 1 2 α = X ∆t f n α .
In this case, a simple Euler explicit scheme is used to discretize the time variable. According to remark 1 we then have

U n+ 1 2 α (x) = X P (ξ)(f n+ 1 2 α (x, ξ))dξ = N K=1 k-1 j=0 U K,n+ 1 2 α,j Φ K j (x) with U K,n+ 1 2 α,j = X P (ξ)(f K,n+ 1 2 j (ξ))dξ.
Step 3: force and source terms. For α = e, i f

n+ 3 4 α = f n+ 1 2 α -∆t q α m α E n+1 N f n+ 3 4 α + ∆tB αβ (f n+ 3 4 e , f n+ 3 4 i ), β ̸ = α and U n+1 α = X P (ξ)(f n+ 3 4 α (ξ))dξ = (ρ n+1 α , ρ n+1 α u n+1 α , E n+1 α ).
We obtain

U n+1 α = U n+ 1 2 α -∆t q α m α E n+1 N U n+1 α + S(T n+1 e , T n+1 i ).
Step 4: coupling with Maxwell-Ampère and Poisson equations.

     q e m e ρ n+1 e + q i m i ρ n+1 i = 0, q e m e ρ n+1 e u n+1 e + q i m i ρ n+1 i u n+1 i = 0.
Hence

u n+1 i = u n+1 e . We set ρ n+1 = ρ n+1 e + ρ n+1 i , u n+1 = u n+1 i = u n+1 e .
As q e = -e and q i = Ze, we get as in the continuous case ρ n+1 e = c e ρ n+1 , ρ n+1 i = c i ρ n+1 . We have analogously to ( 16)

(22)                            ρ n+1 α = ρ n+ 1 2 α , α = e, i, ρ n+1 e u n+1 = ρ n+ 1 2 e u n+ 1 2 e + ∆t q e m e E n+1 ρ n+1 e , ρ n+1 i u n+1 = ρ n+ 1 2 i u n+ 1 2 i + ∆t q i m i E n+1 ρ n+1 i , E n+1 e = E n+ 1 2 e + ∆t q e m e E n+1 ρ n+1 e u n+1 + ∆tν ei (T n+1 i -T n+1 e
),

E n+1 i = E n+ 1 2 i + ∆t q i m i E n+1 ρ n+1 i u n+1 -∆tν ei (T n+1 i -T n+1 e
).

The equations on mass and momentum give ρ n+1 and u n+1 :

(23)    ρ n+1 = ρ n+ 1 2 , ρ n+1 u n+1 = ρ n+ 1 2 e u n+ 1 2 e + ρ n+ 1 2 i u n+ 1 2 i .
We then compute E n+1 as in the continuous case:

c e ρ n+1 u n+1 = ρ n+ 1 2 e u n+ 1 2 e + ∆t q e m e E n+1 ρ n+1 c e , c i ρ n+1 u n+1 = ρ n+ 1 2 i u n+ 1 2 i + ∆t q i m i E n+1 ρ n+1 c i , hence ∆t q e m e E n+1 ρ n+1 e = -∆t q i m i E n+1 ρ n+1 i = -c i ρ n+ 1 2 e u n+ 1 2 e + c e ρ n+ 1 2 i u n+ 1 2 i and (24)
     E n+1 e = E n+ 1 2 e + u n+1 -c i ρ n+ 1 2 e u n+ 1 2 e + c e ρ n+ 1 2 i u n+ 1 2 i + ∆tν ei (T n+1 i -T n+1 e ) E n+1 i = E n+ 1 2 i + u n+1 c i ρ n+ 1 2 e u n+ 1 2 e -c e ρ n+ 1 2 i u n+ 1 2 i + ∆tν ei (T n+1 e -T n+1 i ).
As 24) can be solved explicitly. Finally we obtain U n+1 (x), an approximation of U(t n+1 , x).

T n+1 α = E n+1 α ρ n+1 α - 1 2 (u n+1 ) 2 (γ α -1)m α k B system (

4.3.

Higher order in time. The scheme described in subsection 4.2 can be viewed as a fractional step method which computes the approximate solution U n+1 h at time t n+1 as a function of U n h : (25)

U n+1 = E ∆t (U n ).
This method is only first order in time, even if each step is high order, basically because in the underlying Trotter formula, when two operators A and B do not commute, exp(∆t(A + B)) = exp(∆tA)exp(∆tB) + O(∆t 2 ). We can view this procedure as an explicit RK1 Euler scheme applied to the semi-discretized system obtained by performing the DG spacial discretization of system ( 5):

∂ t U h = G(U h ).
It is important to increase the order in time when one increases the order in space, otherwise you will not observe any significant improvement of the numerical results, see section 5 below.

For this purpose, we use higher order N steps explicit Runge Kutta schemes. The approximate solution U n h being known, we set ( 26)

Y 1 = U n h , V 1 = Y 1 , Y i+1 = E ∆t (V i ), V i+1 = d (N ) i1 U n h + d (N ) i2 Y i+1 , i = 1, . . . , N, and (27)
U n+1 h = V N +1 .
The scheme is defined by a N × 2 matrix

D (N ) = d (N ) ij
. Here we set:

D (1) = (0, 1) , D (2) = 0 1 1 2 1 2 , D (3) =   0 1 3 4 1 4 1 3 2 3   .
The case N = 1 is the explicit Euler method. It is used with a P 0 -DG discretization. The case N = 2 is the RK2 Heun method, used with a P 1 DG discretization. This method is Strong Stability Preserving (SSP) [START_REF] Isherwood | Sigal Strong stability preserving integrating factor Runge-Kutta methods[END_REF]. For the third-order N = 3, we consider the Shu-Osher RK3-SSP scheme [START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications[END_REF], used with P 2 -DG discretization. It is to be noted that the case N = 2 can also be viewed as a DeC method [START_REF] Abgrall | High order schemes for hyperbolic problems using globally continous approximation and avoiding mass matrices[END_REF].

We consider here the following CFL condition:

CF L = ∆t ∆x |v| ≤ 1 2k + 1
where k is the degree of the polynomial [START_REF] Cockburn | Runge-Kutta Discontinuous Galerkin Methods for Convection-Dominated Problems[END_REF].

Numerical results

This section is devoted to the numerical validation of the numerical scheme that is constructed in the previous section. The numerical method is firstly applied to the compressible Euler system (8) and next to the bitemperature Euler system (5).

Compressible Euler system.

5.1.1. Euler isentropic. In this subsection we test the convergence of the scheme for Euler equations (8) by considering the case of isentropic flow i.e when γ = 3 and p = ρ γ , with the initial conditions  

ρ 0 u 0 p 0   =   1 + 0.5 sin(πx) 0 ρ γ 0  
where the domain is Ω = [-1, 1], the final time T = 0.1 and we take CF L = 0.1. We plot in Figure 1 the numerical order of convergence of the method. The considered error is the L 2 -norm.

We observe that we get k + 1 order when we choose a P k -DG discretisation the space with a time discretisation of order k + 1.

  if x ∈ [-5, -4],   ρ 0 u 0 p 0   =   1 + 0.2 sin(5x) 0 1   if x ∈ [-4, 5],
on the domain [-5, 5] and the final time of the problem is T = 1.8. The reference solution represented in Figure 2 is obtained with the P 2 -DG method with 5000 points. We compare in Figure 2 the results obtained with the P 0 , the P 1 and the P 2 method for 512 points with the reference solution. We observe that the oscillations are well captured by the P 2 reconstruction and that the precision increases with the order of the scheme. 5.1.3. Blast waves. In order to highlight the advantages for using a high order method, we consider the test case proposed by Collela and Woodward [START_REF] Collela | The piecewse parabolic method (PPM) for gas dynamical simulations[END_REF] devoted to the compressible Euler system for γ = 1.4. We consider the initial conditions

ρ 0 = 1, u 0 = 0, p 0 = 10 3 1 [0,0.1[ + 10 -2 1 [0.1,0.9[ + 10 2 1 [0.9,1]
. The density and the energy are displayed in Figure 3 for 1000 points in space. As observed in Figure 3, the P 2 reconstruction is able to catch correctly the second pick.

All the previous numerical tests have been tested in [START_REF] Abgrall | High order preserving deferred correction implicit-explicit schemes for kinetic models[END_REF] for a RD scheme and analogous results are obtained. 5.2. Bitemperature Euler system. Next the numerical method is applied to the bitemperature Euler system [START_REF] Bouchut | Construction of BGK Models with a Family of Kinetic Entropies for a Given System of Conservation Laws[END_REF]. In 5.2.1, we perform a convergence study for the DG kinetic scheme for the P 1 and the P 2 reconstruction. In the next test cases, the DG kinetic scheme for the P 2 reconstruction is compared with the kinetic scheme for a finite volume space discretization (FV kinetic scheme) and the Suliciu method developed in [START_REF] Aregba-Driollet | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF]. 5.2.1. Analytical test case. We look for a smooth analytical solution of system [START_REF] Bouchut | Construction of BGK Models with a Family of Kinetic Entropies for a Given System of Conservation Laws[END_REF]. We assume that ρ and u are constant, so that the system reduces to (28)

     ∂ x (ZT e + T i ) = 0, ∂ t T e + u∂ x T e = β e (T i -T e ), ∂ t T i + u∂ x T i = β i (T e -T i ) with β α = ν ei m α (γ α -1) ρc α k B , α ∈ {e, i}.
We choose initial temperatures T e (x, 0) = T e,0 (x), T i (x, 0) = T i,0 (x) such that (29) ∂ x (ZT e,0 + T i,0 ) = 0.

It is easy to compute the solution (T e , T i ) of the last two equations: denoting

β 1 + β 2 = β, µ = γe-1 Z + γ i -1: (30)        T e (x + ut, t) = 1 µ γ e -1 Z
(T e,0 (x) -T i,0 (x))e -βt + (γ i -1)T e,0 (x) + γ e -1 Z T i,0 (x) , T i (x + ut, t) = 1 µ (γ i -1)(T i,0 (x) -T e,0 (x))e -βt + (γ i -1)T e,0 (x) + γ e -1 Z T i,0 (x) .

Then one can observe that if γ e ̸ = γ i the first equation of ( 28) cannot be satisfied unless ∂ x T e,0 = ∂ x T i,0 = 0. Therefore we set γ e = γ i . In that case, if (29) is satisfied then (30) gives a solution of the bitemperature Euler system.

We choose ν ei = 1, T = 0.1, Ω = [0, 1] and CF L = 0.1. We plot in Figure 4 the convergence of the method for the bitemperature Euler system [START_REF] Bouchut | Construction of BGK Models with a Family of Kinetic Entropies for a Given System of Conservation Laws[END_REF]. In contrast to compressible Euler system, the numerical simulation show a convergence of order k when using a P k -DG for the space discretisation and with a (k + 1)-SSP-RK time approximation. The results for electronic and ionic temperatures are displayed in fig. 5 and 6 at time t = 4.09 10 -7 s for 10000 points in space. They are obtained with the Suliciu method, the finite volume kinetic method constructed in [START_REF] Aregba-Driollet | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF] and the DG method with a P 2 reconstruction developed in this paper. The three methods show analogous results. We take γ e = 5/3 and γ i = 7/5. This test case is a rarefaction wave computed at time t = 4.09 10 -7 s for 10000 points in space. In this test case, an analytical solution can be computed. The results displayed in Figure 7 and in Figure 8 compare for electronic and ionic temperatures the P 2 reconstruction with the finite volume kinetic scheme and the Suliciu scheme developed in [START_REF] Aregba-Driollet | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF]. All the scheme show a good agreement with the exact solution.

5.2.4. Stationary shock. We consider here the test case of the stationary shock presented in [START_REF] Aregba-Driollet | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF] in the case ν ei = 100 and with the following parameters:

(31) k B = 1.0, m e = 10 -3 , m i = 1.0, Z = 1.0, γ e = γ i = 5/3

In that case, the left and the right states of the Riemann problem are the following

ρ L = 1.001, u L = 10, T e,L = 1, T i,R = 1, ρ R = 3.640330609, u R = 2.749750250, T e,R = 3, T i,R = 17.5060240977.
The results displayed in Figure 9 and in Figure 10 are computed for 10000 points in space and at time t = 0.005s. 5.2.5. Sod test case. We finally consider the test case of the sod test case in with the parameters (31), ν ei = 0 and a final time 0.05 s

ρ L = 1, u L = 0, T e,L = 1, T i,R = 1, ρ R = 0.125, u R = 0, T e,R = 2, T i,R = 3.
Figures 12 and 11 represent ionic and electronic temperatures for 10000 points in space with Suliciu, kinetic and discontinuous Galerkin method with a P 2 reconstruction. The three schemes show the same results. In particular, the "plateaux" after the shocks have the same amplitude.

The DG kinetic scheme shows some oscillations for the electronic temperatures. However, these oscillations do not propagate.

Conclusion and perspectives

In this paper, we have developed a DG-kinetic scheme for the bitemperature Euler system for different order of space discretization. The principle is to consider a discrete BGK model as in [START_REF] Aregba-Driollet | Modelling and numerical approximation for the nonconservative bitemperature Euler model[END_REF] and to construct a DG discretisation with k-th degree basis for the space discretization and with an (k + 1)-SSP-Runge-Kutta method for the time discretization. Due to the kinetic model, a special treatment has been used to implement the order in time. This method has been illustrated on several test cases and numerical order have been investigated on both conservative Euler equations and on nonconservative bitemperature Euler model. It is to be noted that in the nonconservative case, when shocks occur, the electronic and ionic temperatures cannot be predicted analytically, even in the framework of a Riemann problem as for a conservative system. They can depend on the numerical viscosity. A crucial fact here is that even in the presence of shocks we observe that whatever the order, the DG method converges to the same solutions as the ones obtained in previous articles.

The generalization of the present work to a two dimensional framework as done in [START_REF] Aregba-Driollet | A discrete velocity numerical scheme for the two-dimensional bitemperature Euler system[END_REF] for finite volumes is postponed to a future paper. Moreover in the present case, only the electric field has been taken into account. The question of the presence of the magnetic fields can also be considered as in [START_REF] Brull | A kinetic approach of the bi-temperature Euler model[END_REF] for a transverse magnetic field.

Figure 1 .

 1 Figure 1. Convergence of u for Euler equations

Figure 2 .

 2 Figure 2. Density for the Shu-Osher test case.

Figure 3 .

 3 Figure 3. Density (left) and Energy (right) for the blast waves.

Figure 4 . 3 , γ i = 7 5 .

 435 Figure 4. Convergence of electronic temperature for the bitemperature Euler system

5. 2 . 3 .

 23 Double rarefaction wave. Using similar physical parameters as in subsection 5.2.2 we consider here a rarefaction wave for ν ei = 0 which is a Rieman problem ρ L = 1, u L = -10 5 , T e,L = 2.3 10 7 , T i,L = 2.3 10 6 , ρ R = 1, u R = 10 5 , T e,R = 2.3 10 7 , T i,R = 2.3 10 6 .

Figure 5 .Figure 6 .Figure 7 .

 567 Figure 5. Electronic temperature for the double shock computed with the DG kinetic scheme, the finite volume kinetic scheme, the Suliciu scheme, with 10000 points in space

Figure 8 .Figure 9 .

 89 Figure 8. Ionic temperature for the double rarefaction computed with the DG kinetic scheme, the finite volume kinetic scheme, the Suliciu scheme and the exact solution

Figure 10 .Figure 11 .

 1011 Figure 10. Electronic temperature for the stationary choc computed with the DG kinetic scheme, the finite volume kinetic scheme, the Suliciu scheme

Figure 12 .

 12 Figure 12. Ionic temperature for sod test case computed with the DG kinetic scheme, the finite volume kinetic scheme, the Suliciu scheme computed with 10000 points

  it is a rank one model. This model is compatible with the physical entropy of Euler system.

3.1.2. Example 2: a discrete velocity BGK model. Here