oddsidemargin has been altered.
textheight has been altered.
marginparsep has been altered.
textwidth has been altered.
marginparwidth has been altered.
marginparpush has been altered.
The page layout violates the UAI style. Please do not change the page layout,
or include packages like geometry, savetrees, or fullpage, which change it for
you. We’re not able to reliably undo arbitrary changes to the style. Please
remove the offending package(s), or layout-changing commands and try again.

# Curiosity-Driven Multi-Criteria Hindsight Experience Replay

John B. Lanier
<EMAIL_ADDRESS>
&Stephen McAleer22footnotemark: 2
<EMAIL_ADDRESS>
&Pierre Baldi22footnotemark: 2
<EMAIL_ADDRESS>
Department of Computer Science
University of California, Irvine

###### Abstract

Dealing with sparse rewards is a longstanding challenge in reinforcement
learning. The recent use of hindsight methods have achieved success on a
variety of sparse-reward tasks, but they fail on complex tasks such as
stacking multiple blocks with a robot arm in simulation. Curiosity-driven
exploration using the prediction error of a learned dynamics model as an
intrinsic reward has been shown to be effective for exploring a number of
sparse-reward environments. We present a method that combines hindsight with
curiosity-driven exploration and curriculum learning in order to solve the
challenging sparse-reward block stacking task. We are the first to stack more
than two blocks using only sparse reward without human demonstrations.

##  1 Introduction

Goal-based reinforcement learning has become an important framework for
formulating and solving goal-based sequential decision making tasks. In goal-
based reinforcement learning, the agent’s rewards are usually dependent on
achieving a goal, and it chooses its actions using a goal-conditioned policy.
Goal-conditioned policies can enable a reinforcement learning agent to
generalize to new goals after training on a many different goals in the same
environment ([37]).

Goal-based reinforcement learning environments can be given a binary and
sparse reward that is encountered only when the goal is reached. Defining
reward in this way ensures that if the agent maximizes reward then it also
reaches the user’s intended goal, which is not necessarily true of manually-
shaped dense rewards ([41]). However, sparse rewards are also difficult to
learn from. As the length of a sparse-reward task increases, it becomes less
likely that an agent will discover how to reach its goal through random
exploration ([38]). This problem is exacerbated when a sparse reward depends
on the fulfillment of multiple goals or criteria.

Figure 1: The simulated robotic block stacking environment. A goal consists of
target positions where the blocks need to be placed, shown here as colored
spheres. In the incremental reward environment, the agent receives a reward
for each block being in its target position. In the sparse reward environment,
the agent receives a reward only when every block is in its target position.

Recently, hindsight methods have served as a popular solution to sparse-reward
goal-oriented learning by training an agent on the goals that it actually
reached in addition to those which were intended ([1]). This is done in the
hope that knowledge of how to reach randomly discovered goals will allow an
agent to generalize well enough to find its assigned goals. However, in many
environments, an agent can be asked to reach goals that are very different
from those it may discover by chance, causing such generalization to be
difficult. In these cases, the same sparse reward issues remain, making it
challenging for an agent to learn how to accomplish its given objectives.

Stacking multiple blocks in a simulated robotics environment is a sparse-
reward, goal-based task that highlights shortcomings of hindsight learning.
Multiple-block stacking is too difficult for established hindsight methods
like Deep Deterministic Policy Gradients with Hindsight Experience Replay
(DDPG+HER) to reliably solve without access to human demonstrations ([30]).
Satisfying all criteria of a block stacking goal requires learning multiple
skills to correctly place each block, and the end goals are very different
from those that the agent may discover with random exploration. Achieving
reward by correctly placing all blocks is precarious and requires long chains
of specific actions. Therefore, under a sparse reward, even with hindsight, it
is highly unlikely that an agent will discover the complex sequences of
actions required to place every block in its correct position on the stack.
Our method is the first method that is able to solve sparse-reward block
stacking for more than two blocks without access to human demonstrations.

To solve sparse-reward multi-block stacking without help from demonstration,
we use DDPG+HER combined with curiosity-driven exploration and curriculum
learning. In order to balance improved exploration with exploitation during
training, we introduce a new method of combining data from both curiosity-
based and standard policies in an off-policy fashion. Additionally, we
introduce a form of hindsight experience replay that is more sample efficient
for multi-criteria goal-based environments. We show that the advantages
introduced by each of these methods complement the others, and that the
combination of all of them is necessary to solve the hardest stacking tasks.

###  1.1 Related Work

####  1.1.1 Curiosity-Driven Exploration

We refer to curiosity-driven exploration as any method that attempts to drive
an agent to explore trajectories which it has not visited frequently before,
usually by making the agent pursue some form of exploration related objective
or reward.

Curiosity-driven exploration has been approached by training agents to
maximize information gain ([27, 21]), pursue less visited areas using state
pseudo-counts ([2, 31]), and maximize state empowerment ([18, 29, 24]).

We focus on exploration by performing actions that both challenge and improve
an agent’s ability to model the world ([22, 17, 39, 40]). We approach this by
training a dynamics model on the state transitions that our agent visits and
encouraging the agent to maximize the model’s per-sample error on those
transitions. Assuming a dynamics model is more accurate on transitions that it
has seen frequently before, such an agent seeking to challenge the dynamics
model should be inclined to visit new, rarely before seen state transitions.
Choosing actions to directly challenge an online trained dynamics model has
been shown to result in complex emergent behaviors ([19]). Using a dynamics
model’s error as an RL exploration reward can motivate an agent to seek out
novel states, sometimes solving an environment’s objective without extrinsic
rewards, and combining environmental rewards with a bonus exploration reward
has the potential to increase an agent’s learning speed and end-performance
([32, 4]). On the same note, training a model to predict the output of a
random function from state features and choosing actions to maximize its error
helped achieve state-of-the art performance on the Montezuma’s Revenge Atari
domain ([5]).

####  1.1.2 Curriculum Learning in Goal-Based Tasks

Previous applications of curriculum learning ([3]) to goal based environments
include training on a variety of tolerances for considering goals achieved
([16]), masking certain goal dimensions to allow all such values on an axis to
be sufficient for success ([13]), and generating curricula that walk backwards
from a predefined success state ([14], [28]).

Intrinsically motivated goal exploration processes (IMGEPs) have also been
used to automatically generate goals which maximize learning progress across
one ([15, 33, 25]) or multiple ([6]) tasks.

####  1.1.3 Hindsight methods

Our work builds on Hindsight Experience Replay (HER) ([1]) as a way to
effectively augment goal oriented transition samples for a replay buffer.
Hindsight has also been adapted to policy gradient settings ([37]).

Efforts have been made to increase the efficiency of HER by prioritizing the
sampling of more relevant transitions. This has been done by attributing
higher importance to transitions and trajectories in which more physical work
is done by the agent ([43]), rare goal states are achieved ([42]), or higher
temporal difference error is measured ([9]).

####  1.1.4 Block Stacking

Stacking multiple blocks with sparse rewards has been solved before using
expert demonstration in ([11]) and ([30]). Our work is a direct followup to
the latter, as we solve a similar set of environments without demonstration.

Stacking only 2-blocks with sparse rewards has been solved without
demonstration by training on an automatic curriculum which selects tasks from
a small collection, prioritizing tasks with higher changes in learning
progress ([6]) and by collecting data from multiple policies following
auxiliary objectives to accomplish predefined interesting actions ([38]).
Dense reward robotic block stacking tasks have been solved before using both a
model-based approach, PILCO, ([8, 7]) and by initializing the environment at
intermediate stages of the task ([36]).

###  1.2 Background

####  1.2.1 Reinforcement Learning

We consider the standard reinforcement learning formalism in which an agent
interacts with an environment $E$. The environment is fully observable, and
consists of a set of states $S$, a set of actions $A$, a reward function
$r:S\times A\rightarrow\mathbb{R}$, an initial state distribution $p(s_{0})$
and transition dynamics $p(s_{t+1}|s_{t},a_{t})$. At each timestep $t$, the
agent observes a state $s_{t}$, takes an action $a_{t}$, and receives a reward
$r_{t}$. The agent chooses these actions using a policy $\pi$, which is a
conditional distribution over actions given states. In this paper, we consider
deterministic policies which map directly from states to actions
$\pi:S\rightarrow A$.

The discounted sum of future rewards is defined as the return
$R_{t}=\sum_{i=t}^{T}\gamma^{(i-t)}r_{i}$ over some time horizon $T$ and with
a discounting factor $\gamma\in[0,1]$. We define $\rho^{\pi}$ as the state
visitation distribution when taking actions according to $\pi$. The goal in
reinforcement learning is to learn a policy $\pi$ to maximize the expected
return
$J=\operatorname{\mathbb{E}}_{s_{i}\sim\rho^{\pi},a_{i}\sim\pi,r_{i}\sim
E}{[R_{0}|s_{0}]}$.

The expected return when taking actions according to a specific policy $\pi$
is called the Q-function or action-value function, and is defined as:

$Q^{\pi}(s_{t},a_{t})=\operatorname{\mathbb{E}}_{s_{i>t}\sim\rho^{\pi},r_{i\geq
t}\sim E}{[R_{t}|s_{t},a_{t}]}$ (1)

which can be recursively stated as the Bellman equation:

$Q^{\pi}(s_{t},a_{t})=\operatorname{\mathbb{E}}_{r_{t},s_{t+1}\sim
E}{[r_{t}+\gamma Q^{\pi}(s_{t+1},\pi(s_{t},a_{t}))]}$ (2)

Because $\pi$ is deterministic, the expectation in equation 2 depends only on
the environment, allowing off-policy methods to learn $Q^{\pi}$ while using
transitions generated with some other stochastic policy $\beta$.

####  1.2.2 DDPG

Our work uses the Deep Deterministic Policy Gradients algorithm (DDPG) ([26]),
which is an off-policy, model-free reinforcement learning algorithm designed
for use with deep neural networks in continuous action spaces. DDPG uses an
actor-critic methodology. Two neural networks are trained: a critic $Q:S\times
A\rightarrow\mathbb{R}$ parameterized by $\theta^{Q}$, and an actor serving as
the policy $\pi:S\rightarrow A$, which is updated using the policy gradient to
directly maximize $Q^{\pi}$ with respect to the policy’s parameters
$\theta^{\pi}$:

$\nabla_{\theta^{\pi}}J=\operatorname{\mathbb{E}}_{s_{t}\sim\rho^{\beta}}[\nabla_{\theta^{\pi}}Q(s_{t},a|\theta^{Q})|_{a=\pi(s_{t}|\theta^{\pi})}]$
(3)

This quantity can be estimated with the following:

$\nabla_{\theta^{\pi}}J\approx\frac{1}{N}\sum_{i}\nabla_{a}Q(s_{i},a|\theta^{Q})|_{a=\pi(s_{i})}\nabla_{\theta^{\pi}}\pi(s_{i}|\theta^{\pi})$
(4)

The critic’s parameters $\theta^{Q}$ are updated to minimize the loss:

$L^{crit}=\frac{1}{N}\sum_{\mkern
0.0mui}{(y_{i}-Q(s_{i},a_{i}|\theta^{Q}))^{2}}$ (5)

where

$y_{t}=r_{t}+\gamma Q^{\prime}(s_{t+1},\pi^{\prime}(s_{t+1}))$ (6)

For stability, slower moving target networks $\pi^{\prime}$ and $Q^{\prime}$
are used to calculate $y_{t}$. These network’s parameters are exponential
moving averages of $\theta^{\pi}$ and $\theta^{Q}$ respectively.

DDPG maintains a replay buffer $R$ containing transition samples, which are
tuples $(s_{t},a_{t},r_{t},s_{t+1})$, and alternates between two stages. The
first stage is to gather experience for $R$ by performing rollouts on the
environment, choosing actions from a new policy $\beta=\pi+\epsilon$ where
$\epsilon$ is random. The second stage is to train $\pi$ and $Q$ on batches of
transition samples from $R$.

To efficiently gather experience, we run DDPG in parallel using multiple
workers with synchronized copies of each network, averaging parameters across
workers after each update.

Figure 2: Forward Network Connections for DDPG+HER Learner with Curiosity-
Driven Exploration. During testing, actions are taken by the exploit actor to
maximize return on environmental rewards. During training, actions can be
taken by any actor depending on which objectives we wish to emphasize.

####  1.2.3 DDPG with Goals

In our work, we follow a goal based-framework. A goal $g\in G$ is sampled each
episode, and $\pi$ and $Q$ are conditioned on these goals, making them
$\pi:S\times G\rightarrow A$ and $Q:S\times A\times G\rightarrow\mathbb{R}$.
Furthermore, the replay buffer instead stores transition samples as
$(s_{t}||g,a_{t},r_{t},s_{t+1}||g)$, where the states are each concatenated
with a goal. The environments’ reward functions $r_{t}=r_{env}(s_{t+1},g)$ are
also parameterized on whether a new state meets these goals.

####  1.2.4 Hindsight Experience Replay

In goal-based scenarios, hindsight experience replay increases the sample
efficiency of replay buffer based algorithms like DDPG by adding additional
augmented samples to the replay buffer. In doing so, HER allows the agent to
evaluate its progress not only towards the goals that it was given by the
environment, but also towards those that it actually reached in experience
gathering rollouts, thus giving the agent hindsight.

HER acts by duplicating transition samples before placing them in the replay
buffer, and in those duplicates, augmenting them by replacing the environment-
provided goals with goals that were actually reached later in the same
episode. HER requires the learning algorithm to have access to the reward
function, and the rewards in the augmented samples are updated according to
the newly replaced goals.

HER can also be implemented by expressly storing unmodified transition samples
in the replay buffer and, with a certain probability, augmenting them when
they are sampled from it. We use this method in our work.

###  1.3 Environments

The block stacking environments that we consider in this work are based on the
Fetch robot environments from the OpenAI Gym API ([34]) and are similar those
used in ([30]). We test on separate environments for $n=2$ to $4$ blocks. Each
episode, target block locations are initialized in a stack somewhere on the
surface of a table. The $n$ blocks are initialized at random locations on the
table away from the target stack location. The blocks are uniquely labeled,
and each block always goes to the same vertical position in the stack. The
agent has $25n$ timesteps before the environment resets. These environments
are fully observable, and observations include the claw state and full
position, rotation, and velocity for both the robot’s gripper and each block.
These environments’ goals specify the target positions of each block, and a
block is considered correctly placed if its position is within an error
tolerance $e$ from its target position. Actions are continuous and control the
robot gripper’s movement in 3 dimensions as well as the state of the claw.

Similar to ([30]), we consider two sparse reward formulations with these
environments: binary and incremental. We can provide a single binary reward
when the goal is fully achieved upon correctly placing all blocks:

$r_{t}^{binary}=\begin{cases}0&\text{all blocks in place}\\\
-1&\text{otherwise}\end{cases}$ (7)

We also consider incremental rewards for each block correctly placed:

$r_{t}^{incremental}=\text{no. of blocks in place}-\text{no. of blocks}$ (8)

In both cases, we also add $1$ to the reward for moving the gripper away from
the blocks once they are all correctly placed, but only the correct placement
of all the blocks determines whether a goal has been achieved.

##  2 Methods

To solve multi-block stacking with both incremental and binary rewards, we use
three methods to improve the performance of a standard multi-worker DDPG+HER
learner: curiosity-driven exploration, multi-criteria HER, and curriculum
learning.

First, we incorporate curiosity-driven exploration by training a forward
dynamics model on state transitions visited by the agent and treating the
dynamics model’s prediction error on these transitions as an exploration
reward. In order to have a certain portion of workers explore while others
exploit, we train three separate policies to maximize exploration rewards,
environmental rewards, and a weighted combination of both. Experience from
rollouts is shared among each network regardless of which policy collected it.
By doing so, we can use different policies at training time than at test time.

Second, we introduce a form of hindsight experience replay better suited for
multi-criteria goal-based environments, where a criteria in our environment is
defined as the position of a specific block. Our method randomly performs the
goal replacement operation on each independent criteria in a goal rather than
on an entire goal at once, decoupling the individual effects of each criteria
on the reward function and providing higher sample efficiency.

Third, we use curriculum learning by training the agent on two easier skill-
building environments before training on the target multi-block stacking task.

###  2.1 Curiosity Driven Exploration with Multiple Policies

We use curiosity-driven exploration to encourage an agent to visit transitions
which are novel and surprising to it. We define an auxiliary exploration
reward in addition to environmental reward, and we train separate critics for
each. The explore critic $Q_{e}$ predicts the action-value function for
exploration rewards, and the exploit critic $Q_{r}$ predicts the action-value
function for environmental rewards. We train three actor policies
$\pi_{e},\pi_{r},\pi_{c}$ which respectively maximize exploration rewards,
environmental rewards, and a weighted combination of both. By training
separate polices, we can make our agent pursue multiple and various objectives
at training time and maximize only environmental rewards at test time.

We maintain a forward dynamics neural network $D:S\times A\rightarrow S$
parameterized by $\theta^{d}$ to predict the next observation given the
current observation and action, and we train it on the same transition samples
from the replay buffer as our agent at each DDPG update step. For each
transition sample trained on, an exploration reward for the sample is defined
as the squared error between the predicted next state
$D(s_{t},a_{t}|\theta^{d})$ and the actual next state $s_{t+1}$. The minibatch
loss function for $D$ and the exploration reward is formulated as:

$L_{d}=\frac{1}{N}\sum_{\mkern 0.0mui}r_{i}^{explore}=\frac{1}{N}\sum_{\mkern
0.0mui}(s_{i+1}-D(s_{i},a_{i}|\theta^{d}))^{2}$ (9)

By passing this error $r_{i}^{explore}$ to our agent as an exploration reward
to maximize, we encourage our agent to pursue transitions that are difficult
to predict and unlike transitions currently in the replay buffer.

With two separate reward sources, we group our multiple actors and critics
into two DDPG actor-critic pairs. On exploration reward, we train $Q_{e}$ and
$\pi_{e}$ as our explore actor-critic pair. On environmental reward, we train
$Q_{r}$ and $\pi_{r}$ as our exploit actor-critic pair. In addition to these
actor-critic pairs, we also train $\pi_{c}$ as our combined actor. $\pi_{c}$
pursues both exploration and environmental reward by maximizing a weighted
average of both critics’ action-value functions.

We train multiple actors towards different objectives so that we can assign a
portion of our workers to follow an exploration related policy $\pi_{e}$ or
$\pi_{c}$ while the rest follow the exploit policy $\pi_{r}$. Doing so allows
us to diversify the experience gathered and make less sacrifices toward either
exploration or exploitation objectives than if we were to only ever choose
actions which maximize a weighted combination of the two. With multiple
actors, we can specialize our workers and maximize both environmental and
exploration rewards when gathering experience and then use $\pi_{r}$ at test
time to solely maximize environmental rewards. Below we describe these
networks in more detail.

####  2.1.1 Exploit Actor and Critic

$\pi_{r}:S\times G\rightarrow A$ and $Q_{r}:S\times G\times
A\rightarrow\mathbb{R}$ together form the exploit actor-critic pair, which is
trained on the normal DDPG+HER goal-based RL objective for maximizing return
on environmental rewards conditioned on goals. This actor-critic pair follows
the same configuration and update rules as what would be used in vanilla
DDPG+HER. The loss function for $Q_{r}$ to minimize with respect to its
parameters $\theta^{Q}_{r}$ is:

$L_{r}^{crit}=\operatorname{\mathbb{E}}_{s_{t}\sim\rho^{\beta},a_{t}\sim\beta,r_{t}\sim
E}(y_{t}^{r}-Q_{r}(s_{t},a_{t},g|\theta^{Q}_{r}))^{2}$ (10)

where $y_{t}^{r}$ is calculated using the target exploit actor and critic
$\pi^{\prime}_{r}$ and $Q^{\prime}_{r}$:

$y_{t}^{r}=r_{t}^{env}+\gamma
Q^{\prime}_{r}(s_{t+1},\pi^{\prime}_{r}(s_{t+1},g),g)$ (11)

$\pi_{r}$ is updated using the standard goal-based DDPG policy gradient to
maximize $Q_{r}$ with respect to $\pi_{r}$’s parameters $\theta^{\pi}_{r}$.

####  2.1.2 Explore Actor and Critic

$\pi_{e}:S\rightarrow A$ and $Q_{e}:S\times A\rightarrow\mathbb{R}$ together
form the explore actor-critic pair, which is trained on the objective of
maximizing return on exploration reward. Goals do not affect exploration
rewards and are not factored in these calculations. The loss function for
$Q_{e}$ to minimize with respect to its parameters $\theta^{Q}_{e}$ is:

$L_{e}^{crit}=\operatorname{\mathbb{E}}_{s_{t}\sim\rho^{\beta},a_{t}\sim\beta,r_{t}\sim
E}(y_{t}^{e}-Q_{e}(s_{t},a_{t}|\theta^{Q}_{e}))^{2}$ (12)

where $y_{t}^{e}$ is calculated using the target explore actor and critic
$\pi^{\prime}_{e}$ and $Q^{\prime}_{e}$:

$y_{t}^{e}=r_{t}^{explore}+\gamma
Q^{\prime}_{e}(s_{t+1},\pi^{\prime}_{e}(s_{t+1}))$ (13)

Likewise, $\pi_{e}$ is updated using the standard DDPG policy gradient to
maximize $Q_{e}$ with respect to $\pi_{e}$’s parameters $\theta^{\pi}_{e}$.

####  2.1.3 Combined Actor and POP-ART

Once our agent has an idea of how to find environmental rewards, it is usually
more advantageous to explore trajectories close to what actually results in
those rewards. Towards this end, we train our combined actor $\pi_{c}:S\times
G\rightarrow A$ to choose actions that maximize both the exploration and
exploitation objectives simultaneously. $\pi_{c}$ outputs actions that
maximize the weighted combination of both $Q_{e}$ and $Q_{r}$’s action-value
functions.

We intend to maintain a normalized scale at which to compare the return
estimates from $Q_{e}$ and $Q_{r}$ so that we can intuitively weight their
relative importance to $\pi_{c}$. We also need to account for the fact that
the magnitude of both action-value functions may change drastically over the
course of training. This is especially true of $Q_{e}$ which predicts the
return from the moving exploration reward function. To accomplish both of
these goals, each of the targets $y^{e}$ and $y^{r}$ for $Q_{e}$ and $Q_{r}$
are adaptively normalized such that we also maintain normalized versions
$n_{e}^{Q}$ and $n_{r}^{Q}$ of both action-value functions with the same
relative scale at all times. We can then intuitively weight the relative
importance of $n_{e}^{Q}$ and $n_{r}^{Q}$ for $\pi_{c}$ to maximize. In our
case, we weight them equally.

To do this, we use PopArt normalization ([20]), which allows us to adaptively
normalize our critics’ targets without hurting the accuracy of our
predictions. Here we only sketch PopArt informally. See ([20]) for more
details. For each critic target $y^{e}$ and $y^{r}$ we keep an online estimate
of its mean and standard deviation $\sigma_{e},\mu_{e}$ and
$\sigma_{r},\mu_{r}$. We then parameterize $Q_{e}$ and $Q_{r}$ as linear
transformations of the suitably normalized action-value functions $n_{e}^{Q}$
and $n_{r}^{Q}$:

$\begin{split}Q_{e}(s,a|\theta^{Q}_{e})&=\sigma_{e}n_{e}^{Q}(s,a|\theta^{Q}_{e})+\mu_{e}\\\
Q_{r}(s,a,g|\theta^{Q}_{r})&=\sigma_{r}n_{r}^{Q}(s,a,g|\theta^{Q}_{r})+\mu_{r}\end{split}$
(14)

$n_{e}^{Q}$ and $n_{r}^{Q}$ are the actual networks that we train
parameterized by $\theta^{Q}_{e}$ and $\theta^{Q}_{r}$, and when the
statistics $\sigma_{e},\mu_{e}$ and $\sigma_{r},\mu_{r}$ are updated, the top
layers of $n_{e}^{Q}$ and $n_{r}^{Q}$ are also adjusted to preserve equation
14. Similarly, our target critics $Q_{e}^{\prime},Q_{r}^{\prime}$ are
equivalent linear transformations of target normalized critic networks
${n_{e}^{Q}}^{\prime},{n_{r}^{Q}}^{\prime}$.

In our implementation, the scale-invariant loss functions for each of our two
critic networks $n_{e}^{Q}$ and $n_{r}^{Q}$ are:

$L_{e}^{crit}=\frac{1}{N}\sum_{\mkern
0.0mui}{(\frac{y_{i}^{e}-\mu_{j}}{\sigma_{e}}-n_{e}^{Q}(s,a|\theta^{Q}_{e}))^{2}}$
(15) $L_{r}^{crit}=\frac{1}{N}\sum_{\mkern
0.0mui}{(\frac{y_{i}^{r}-\mu_{j}}{\sigma_{r}}-n_{r}^{Q}(s,a,g|\theta^{Q}_{r}))^{2}}$
(16)

Given: Worker policies
$\pi_{0},\pi_{1},...,\pi_{w}|\pi_{i}\in\\{\pi_{e},\pi_{r},\pi_{c}\\}$

Randomly initialize networks $D,n^{Q}_{e},n^{Q}_{r},\pi_{e},\pi_{r},\pi_{c}$

Initialize target networks
${n_{e}^{Q}}^{\prime},{n_{r}^{Q}}^{\prime},\pi_{e}^{\prime},\pi_{r}^{\prime},\pi_{c}^{\prime}$

_(Execute for each parallel worker $i$):_

Initialize replay buffer $R$

for  _Epoch = $1,...,E$_ do

for  _Cycle = $1,C$_ do

for  _Episode = $1,M$_ do

Sample $\epsilon$ and set $\beta\leftarrow\pi_{i}+\epsilon$

Receive initial state $s_{0}$ and goal $g$

for  _t = $0,T$_ do

Select action $a_{t}=\beta(s_{t},g)$ with noise

Take action $a_{t}$, receive $r_{t},s_{t+1}$

Store $(s_{t}||g,a_{t},r_{t},s_{t+1}||g)$ in $R$

end for

end for

for  _Batch = $1,...,K$_ do

Sample batch $B$ from $R$ with HER augmentations

Train $D$ on $B$

foreach  _transition sample $j$ in $B$_ do

Set $r_{j}^{explore}$ and add it to sample

end foreach

Train $n^{Q}_{e},n^{Q}_{r},\pi_{e},\pi_{r},\pi_{c}$ on $B$

Update target networks

Average network parameters over workers

end for

end for

Test performance on episodes using $\pi_{r}$

end for

Algorithm 1 DDPG+HER with Curiosity

By training $n_{e}^{Q}$ and $n_{r}^{Q}$ to predict normalized action-value
functions, we can update $\pi_{c}$ to jointly maximize the evaluation from
both the explore and exploit critics with equal importance:

$\nabla_{\theta^{\pi}_{c}}J=\operatorname{\mathbb{E}}_{s_{t}\sim\rho^{\beta}}[\nabla_{\theta^{\pi}_{c}}Q_{c}^{norm}(s_{t},a,g)|_{a=\pi_{c}(s_{t},g|\theta^{\pi}_{c})}]$
(17)

where

$Q_{c}^{norm}(s_{t},a,g)=\frac{n_{e}^{Q}(s_{t},a|\theta^{Q}_{e})+n_{r}^{Q}(s_{t},a,g|\theta^{Q}_{r})}{2}$

Then for our three actors $\pi_{c},\pi_{e},\pi_{r}$ the implemented policy
gradient update rules are:

$\nabla_{\theta^{\pi}_{c}}J\approx\frac{1}{N}\sum_{i}\nabla_{a}Q_{c}^{norm}(s_{i},a,g)\nabla_{\theta^{\pi}_{c}}a|_{a=\pi_{c}(s_{i},g|\theta^{\pi}_{c})}$
(18)
$\nabla_{\theta^{\pi}_{e}}J\approx\frac{1}{N}\sum_{i}\nabla_{a}n_{e}^{Q}(s_{i},a|\theta^{Q}_{e})\nabla_{\theta^{\pi}_{e}}a|_{a=\pi_{e}(s_{i}|\theta^{\pi}_{e})}$
(19)
$\nabla_{\theta^{\pi}_{r}}J\approx\frac{1}{N}\sum_{i}\nabla_{a}n_{r}^{Q}(s_{i},a,g|\theta^{Q}_{r})\nabla_{\theta^{\pi}_{r}}a|_{a=\pi_{r}(s_{i},g|\theta^{\pi}_{r})}$
(20)

In our experiments, when we used curiosity-driven learning, we chose actions
using the combined policy $\pi_{c}$ instead of the pure explore policy
$\pi_{e}$. The pure explore policy $\pi_{e}$ is still useful to train the
explore critic $Q_{e}$ which is then used to train the combined policy
$\pi_{c}$.

###  2.2 Multi-Criteria Hindsight

We define the multiple criteria in a goal as the individual target block
positions that the goal specifies. In general, for other environments,
criteria can be elements of a goal that require learning separate skills to
accomplish. To increase the quality of data provided by hindsight experience
replay, we randomly perform the hindsight goal replacement operation
independently on each criteria in a goal that we are augmenting. This is done
instead of replacing the entire goal with one reached later in the same
episode.

Our method provides more transition samples to the agent with goals that are
only partially completed later in the same episode. With normal HER, all
hindsight augmented samples that the agent receives contain goals in which all
criteria were satisfied at a later timestep. With multi-criteria HER, the
agent will still receive a portion of goals that it later satisfied
completely, and it will also receive many goals that it later only satisfied
some criteria for.

Given:

  * •

an augmentation probability $z$

  * •

a Replay Buffer $R$

Sample a batch $B$ from $R$

foreach  _transition sample $(s_{t}||g,a_{t},r_{t},s_{t+1}||g)$ in $B$_ do

foreach  _target block position $p_{i}$ in $g$_ do

Sample $u\sim U(0,1)$

if  _$u <z$_ then

Sample a position $p^{\prime}_{i}$ that block $i$ reached later in the same
episode.

else

$p_{i}^{\prime}\leftarrow p_{i}$

end if

end foreach

$g^{\prime}\leftarrow p_{0}^{\prime}||p_{1}^{\prime}||...||p_{n}^{\prime}$

$r_{t}^{\prime}\leftarrow r(s_{t+1},g^{\prime})$

replace transition sample w/
$(s_{t}||g^{\prime},a_{t},r_{t}^{\prime},s_{t+1}||g^{\prime})$

end foreach

Pass $B$ with augmented transition samples to neural networks for training

Algorithm 2 Multi-Criteria HER Augmentation Step

In our experiments, for both binary and incremental reward formulations, using
multi-criteria HER results in significant, if not critical, improvements to
sample efficiency and inter-task generalization.

###  2.3 Curriculum

Although multi-criteria hindsight sampling allows for more sample-efficient
learning and curiosity driven exploration assists in reward discovery, it was
necessary to employ curriculum learning to successfully solve multi-block
stacking with sparse rewards. Training was broken into three stages, in which
reaching a threshold success rate in a previous stage caused the agent to
transition to the next stage. At the beginning of each stage, the DDPG
algorithm was restarted, transferring only the weights of each network from a
previous stage and reinitializing an empty replay buffer.

In _stage 1_ , the agent trains on a non-stacking version of the block
environment to help it learn fundamental skills that are transferable to the
target block stacking task. The stage 1 environment is initialized with the
same number of randomly placed blocks as the target stacking task. Each
episode, rather than in a stack, the blocks’ target positions are randomly
placed on the surface of the table. A single block’s target position may also
be in the air instead. This stage is designed to provide less challenging
tasks in which the agent can more easily discover the basic block manipulation
mechanics necessary for completing the harder stacking task.

In _stage 2_ , the agent trains on actual block stacking with the environment
initialized at various intermediate stages of completion. At each episode, a
random number of the n blocks between $0$ and $n-1$ are initialized already in
the correct position on the stack. Some targets may also still be on the table
rather than on the stack.

Finally, in _stage 3_ , the agent trains on the target block stacking task, in
which all blocks were consistently initialized on the table, away from their
target locations on the stack.

Figure 3: Success rates and per-episode rewards for block stacking with binary
rewards. Success rates and per-episode reward values shown here are for the
respective curriculum stage’s task in which they are measured.

##  3 Experiments

Figure 4: Success rates and per-episode rewards for block stacking with
incremental rewards. Success rates and per-episode reward values shown here
are for the respective curriculum stage’s task in which they are measured.

In this section, we show our method’s performance on the block stacking tasks
using both binary and incremental rewards. Stacking 2, 3, and 4 blocks were
tested. Ablations are also shown to demonstrate the effectiveness of each our
methods. We performed tests using the following configurations:

_All 3_ : Multi-criteria HER, curiosity-driven exploration, and curriculum
learning are all used with our DDPG+HER learner.

_No Curiosity_ : Multi-criteria HER and curriculum learning are used, but all
actions are chosen using $\pi_{r}$. An explore actor-critic and combined actor
are not trained.

_No Multi-Criteria_ : Curiosity-driven exploration and curriculum learning are
used, however HER is done in the original way as defined in ([1]).

_No Curriculum_ : Multi-criteria HER and curiosity-driven exploration are
used, however the agent only trains on stage 3 of the curriculum, which is the
actual target task of multi-block stacking.

_Vanilla DDPG+HER_ : None of the three methods introduced in section 2 are
used. This is the original DDPG+HER algorithm as in ([1]). All actions are
chosen using $\pi_{r}$.

We trained our agent using 8 to 32 parallel workers depending on the
difficulty of the task. When curiosity driven-exploration was used, during
experience gathering rollouts, we assigned half of the workers to take actions
using $\pi_{c}$, and the other half using $\pi_{r}$. Also during experience
gathering, we applied parameter-space noise ([35]) to the actor networks used
and gaussian noise to the actions chosen. Comprehensive hyper-parameter
details can be found in the supplementary materials associated with this
paper.

Success rates and per-episode reward were measured during discrete testing
phases in every epoch of training. During testing, actions were always chosen
using $\pi_{r}$. An episode was considered successful if its goal $g$ was
achieved during the episode’s final state $s_{T}$.

Success rate and per-episode reward statistics were a moving average over the
last 100 episodes tested on. These two statistics are shown as a function of
total environment interaction timesteps for binary reward tasks in Figure 3
and for incremental reward tasks in Figure 4.

Table 1: Highest Success Rates with Binary Rewards over 100 Episode Sliding Window Method | Stack-2 | Stack-3 | Stack-4
---|---|---|---
All-3 | 1.00 | 0.95 | 0.00
No Curiosity | 1.00 | 0.00 | -
No Multi-Criteria | 1.00 | 0.00 | -
No Curriculum | 0.00 | - | -
Curriculum Only | 0.00 | - | -
Vanilla DDPG+HER | 0.00 | - | -
Table 2: Highest Success Rates with Incremental Rewards over 100 Episode Window Method | Stack-2 | Stack-3 | Stack-4
---|---|---|---
All-3 | 1.00 | 0.98 | 0.79
No Curiosity | 0.99 | 0.94 | -
No Multi-Criteria | 0.00 | 0.00 | -
No Curriculum | 0.00 | 0.00 | -
Curriculum Only | 0.00 | - | -
Vanilla DDPG+HER | 0.00 | - | -

Tables 2 and 1 show the highest success rates for each method on the target
bock stacking tasks with binary and incremental reward formulations. For
methods that used a curriculum but did not reach the target task in the third
stage, the final network weights were used to test performance at the target
block stacking task anyways.

Vanilla DDPG+HER was unable to solve block-stacking with any number of blocks
and either reward formulation.

Stacking 2 blocks with either reward formulation was solvable as long as the
agent trained on the curriculum and used multi-criteria HER. Using curiosity-
driven exploration without multi-criteria HER allowed the agent to make
progress on stage 1 of the curriculum, but when incremental rewards were
given, it failed to generalize between the stage 1 task and the stage 2 task
well enough to continue learning.

Stacking 3 blocks with incremental rewards required the use of both curriculum
learning and multi-criteria HER to solve. With binary rewards, stacking 3
blocks required the use of all three methods, as curiosity-driven exploration
was necessary to find a reward signal.

Due to limits on computational resources, stacking 4 blocks was only tested
with all three methods to measure the best possible performance. No progress
was made on the binary reward environment, and in the incremental reward
environment, a max success rate of $0.79$ was reached on the target block
stacking task.

Multi-criteria HER provided clear improvements to sample efficiency, and was
necessary for stacking three or more blocks.

Agents with curiosity-driven exploration learned to solve tasks with less
environment interactions than those without. With incremental rewards, block
stacking was easy enough to be solved without curiosity-driven exploration,
however with binary rewards, curiosity was required to solve stacking 3
blocks.

Finally, curriculum learning was necessary for any of the stacking tasks, as
no method could progress on the target stacking task without first training on
stages 1 and 2.

##  4 Conclusion

By combining curiosity-based exploration with curriculum learning and multi-
criteria HER, we are the first to solve sparse reward multi-block stacking
without demonstrations. This work shows that even very challenging sparse
reward environments can be solved through a combination of existing
techniques. In future work, other methods of intrinsic exploration such as Go-
Explore ([12]) might prove more effective than curiosity-driven exploration
when combined with HER. In our work, we generate curricula in a hand-designed
way based on domain knowledge. This might not be possible in more complex
domains such as real-world robotics. Because of this, further research in
automatically generating curricula is likely to be a fruitful direction when
combined with HER.

#### References

## References

  * [1] Marcin Andrychowicz, Filip Wolski, Alex Ray, Jonas Schneider, Rachel Fong, Peter Welinder, Bob McGrew, Josh Tobin, OpenAI Pieter Abbeel and Wojciech Zaremba  “Hindsight experience replay”  In _Advances in Neural Information Processing Systems_ , 2017, pp. 5048–5058
  * [2] Marc Bellemare, Sriram Srinivasan, Georg Ostrovski, Tom Schaul, David Saxton and Remi Munos  “Unifying count-based exploration and intrinsic motivation”  In _Advances in Neural Information Processing Systems_ , 2016, pp. 1471–1479
  * [3] Yoshua Bengio, Jérôme Louradour, Ronan Collobert and Jason Weston  “Curriculum learning”  In _Proceedings of the 26th annual international conference on machine learning_ , 2009, pp. 41–48  ACM
  * [4] Yuri Burda, Harri Edwards, Deepak Pathak, Amos Storkey, Trevor Darrell and Alexei A Efros  “Large-scale study of curiosity-driven learning”  In _arXiv preprint arXiv:1808.04355_ , 2018
  * [5] Yuri Burda, Harri Edwards, Amos Storkey and Oleg Klimov  “Exploration by random network distillation”  In _arXiv preprint arXiv:1810.12894_ , 2018
  * [6] Cédric Colas, Olivier Sigaud and Pierre-Yves Oudeyer  “CURIOUS: Intrinsically Motivated Multi-Task, Multi-Goal Reinforcement Learning”  In _arXiv preprint arXiv:1810.06284_ , 2018
  * [7] Marc Peter Deisenroth and Carl Edward Rasmussen  “PILCO: A model-based and data-efficient approach to policy search”  In _Proceedings of the 28th International Conference on machine learning (ICML-11)_ , 2011, pp. 465–472
  * [8] Marc Peter Deisenroth, Carl Edward Rasmussen and Dieter Fox  “Learning to control a low-cost manipulator using data-efficient reinforcement learning”  In _Robotics: Science and Systems VII_ MIT Press, 2011, pp. 78–85
  * [9] Ameet Deshpande and Srikanth Sarma  “Improvements on Hindsight Learning”  In _arXiv preprint arXiv:1809.06719_ , 2018
  * [10] Prafulla Dhariwal, Christopher Hesse, Oleg Klimov, Alex Nichol, Matthias Plappert, Alec Radford, John Schulman, Szymon Sidor, Yuhuai Wu and Peter Zhokhov  “Openai baselines” https://github.com/openai/baselines, 2018
  * [11] Yan Duan, Marcin Andrychowicz, Bradly Stadie, OpenAI Jonathan Ho, Jonas Schneider, Ilya Sutskever, Pieter Abbeel and Wojciech Zaremba  “One-shot imitation learning”  In _Advances in neural information processing systems_ , 2017, pp. 1087–1098
  * [12] Adrien Ecoffet, Joost Huizinga, Joel Lehman, Kenneth O Stanley and Jeff Clune  “Go-Explore: a New Approach for Hard-Exploration Problems”  In _arXiv preprint arXiv:1901.10995_ , 2019
  * [13] Manfred Eppe, Sven Magg and Stefan Wermter  “Curriculum goal masking for continuous deep reinforcement learning”  In _arXiv preprint arXiv:1809.06146_ , 2018
  * [14] Carlos Florensa, David Held, Markus Wulfmeier, Michael Zhang and Pieter Abbeel  “Reverse curriculum generation for reinforcement learning”  In _arXiv preprint arXiv:1707.05300_ , 2017
  * [15] Sébastien Forestier, Yoan Mollard and Pierre-Yves Oudeyer  “Intrinsically motivated goal exploration processes with automatic curriculum learning”  In _arXiv preprint arXiv:1708.02190_ , 2017
  * [16] Pierre Fournier, Olivier Sigaud, Mohamed Chetouani and Pierre-Yves Oudeyer  “Accuracy-based Curriculum Learning in Deep Reinforcement Learning”  In _arXiv preprint arXiv:1806.09614_ , 2018
  * [17] Goren Gordon and Ehud Ahissar  “Hierarchical curiosity loops and active sensing”  In _Neural Networks_ 32 Elsevier, 2012, pp. 119–129
  * [18] Karol Gregor, Danilo Jimenez Rezende and Daan Wierstra  “Variational intrinsic control”  In _arXiv preprint arXiv:1611.07507_ , 2016
  * [19] Nick Haber, Damian Mrowca, Li Fei-Fei and Daniel LK Yamins  “Learning to Play with Intrinsically-Motivated Self-Aware Agents”  In _arXiv preprint arXiv:1802.07442_ , 2018
  * [20] Hado P Hasselt, Arthur Guez, Matteo Hessel, Volodymyr Mnih and David Silver  “Learning values across many orders of magnitude”  In _Advances in Neural Information Processing Systems_ , 2016, pp. 4287–4295
  * [21] Rein Houthooft, Xi Chen, Yan Duan, John Schulman, Filip De Turck and Pieter Abbeel  “Vime: Variational information maximizing exploration”  In _Advances in Neural Information Processing Systems_ , 2016, pp. 1109–1117
  * [22] Rituraj Kaushik, Konstantinos Chatzilygeroudis and Jean-Baptiste Mouret  “Multi-objective Model-based Policy Search for Data-efficient Learning with Sparse Rewards”  In _arXiv preprint arXiv:1806.09351_ , 2018
  * [23] Diederik P Kingma and Jimmy Ba  “Adam: A method for stochastic optimization”  In _arXiv preprint arXiv:1412.6980_ , 2014
  * [24] Alexander S Klyubin, Daniel Polani and Chrystopher L Nehaniv  “Empowerment: A universal agent-centric measure of control”  In _Evolutionary Computation, 2005. The 2005 IEEE Congress on_ 1, 2005, pp. 128–135  IEEE
  * [25] Adrien Laversanne-Finot, Alexandre Péré and Pierre-Yves Oudeyer  “Curiosity driven exploration of learned disentangled goal spaces”  In _arXiv preprint arXiv:1807.01521_ , 2018
  * [26] Timothy P Lillicrap, Jonathan J Hunt, Alexander Pritzel, Nicolas Heess, Tom Erez, Yuval Tassa, David Silver and Daan Wierstra  “Continuous control with deep reinforcement learning”  In _arXiv preprint arXiv:1509.02971_ , 2015
  * [27] Daniel Ying-Jeh Little and Friedrich Tobias Sommer  “Learning and exploration in action-perception loops”  In _Frontiers in neural circuits_ 7 Frontiers, 2013, pp. 37
  * [28] Stephen McAleer, Forest Agostinelli, Alexander Shmakov and Pierre Baldi  “Solving the Rubik’s Cube with Approximate Policy Iteration”  In _International Conference on Learning Representations_ , 2019
  * [29] Shakir Mohamed and Danilo Jimenez Rezende  “Variational information maximisation for intrinsically motivated reinforcement learning”  In _Advances in neural information processing systems_ , 2015, pp. 2125–2133
  * [30] Ashvin Nair, Bob McGrew, Marcin Andrychowicz, Wojciech Zaremba and Pieter Abbeel  “Overcoming exploration in reinforcement learning with demonstrations”  In _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , 2018, pp. 6292–6299  IEEE
  * [31] Georg Ostrovski, Marc G Bellemare, Aaron van den Oord and Rémi Munos  “Count-based exploration with neural density models”  In _arXiv preprint arXiv:1703.01310_ , 2017
  * [32] Deepak Pathak, Pulkit Agrawal, Alexei A Efros and Trevor Darrell  “Curiosity-driven exploration by self-supervised prediction”  In _International Conference on Machine Learning (ICML)_ 2017, 2017
  * [33] Alexandre Péré, Sébastien Forestier, Olivier Sigaud and Pierre-Yves Oudeyer  “Unsupervised Learning of Goal Spaces for Intrinsically Motivated Goal Exploration”  In _arXiv preprint arXiv:1803.00781_ , 2018
  * [34] Matthias Plappert, Marcin Andrychowicz, Alex Ray, Bob McGrew, Bowen Baker, Glenn Powell, Jonas Schneider, Josh Tobin, Maciek Chociej and Peter Welinder  “Multi-goal reinforcement learning: Challenging robotics environments and request for research”  In _arXiv preprint arXiv:1802.09464_ , 2018
  * [35] Matthias Plappert, Rein Houthooft, Prafulla Dhariwal, Szymon Sidor, Richard Y Chen, Xi Chen, Tamim Asfour, Pieter Abbeel and Marcin Andrychowicz  “Parameter space noise for exploration”  In _arXiv preprint arXiv:1706.01905_ , 2017
  * [36] Ivaylo Popov, Nicolas Heess, Timothy Lillicrap, Roland Hafner, Gabriel Barth-Maron, Matej Vecerik, Thomas Lampe, Yuval Tassa, Tom Erez and Martin Riedmiller  “Data-efficient deep reinforcement learning for dexterous manipulation”  In _arXiv preprint arXiv:1704.03073_ , 2017
  * [37] Paulo Rauber, Avinash Ummadisingu, Filipe Mutz and Juergen Schmidhuber  “Hindsight policy gradients”  In _arXiv preprint arXiv:1711.06006_ , 2017
  * [38] Martin Riedmiller, Roland Hafner, Thomas Lampe, Michael Neunert, Jonas Degrave, Tom Van de Wiele, Volodymyr Mnih, Nicolas Heess and Jost Tobias Springenberg  “Learning by Playing-Solving Sparse Reward Tasks from Scratch”  In _arXiv preprint arXiv:1802.10567_ , 2018
  * [39] Jürgen Schmidhuber  “A possibility for implementing curiosity and boredom in model-building neural controllers”  In _Proc. of the international conference on simulation of adaptive behavior: From animals to animats_ , 1991, pp. 222–227
  * [40] Jürgen Schmidhuber  “Formal theory of creativity, fun, and intrinsic motivation (1990–2010)”  In _IEEE Transactions on Autonomous Mental Development_ 2.3 IEEE, 2010, pp. 230–247
  * [41] Matej Večerik, Todd Hester, Jonathan Scholz, Fumin Wang, Olivier Pietquin, Bilal Piot, Nicolas Heess, Thomas Rothörl, Thomas Lampe and Martin Riedmiller  “Leveraging demonstrations for deep reinforcement learning on robotics problems with sparse rewards”  In _arXiv preprint arXiv:1707.08817_ , 2017
  * [42] Rui Zhao and Volker Tresp  “Curiosity-Driven Experience Prioritization via Density Estimation”, 2018
  * [43] Rui Zhao and Volker Tresp  “Energy-Based Hindsight Experience Prioritization”  In _arXiv preprint arXiv:1810.01363_ , 2018

##  Appendix A Links

A video showcasing this project is available at
https://youtu.be/stZX4o0H8Ro

Code for our modified DDPG Learner is available at
https://github.com/CDMCH/ddpg-with-curiosity-and-multi-criteria-her
and code for our block stacking environments is available at
https://github.com/CDMCH/gym-fetch-stack

Our DDPG learner uses code modified from the OpenAI baselines repository
([10]).

##  Appendix B Experiment Details

Observation and goal network inputs were normalized to have a mean of zero and
standard deviation of one. Once normalized, they were also clipped to the
range [-5, 5].

All networks were fully connected with 3 hidden layers and 256 hidden units in
each layer. Hidden layers used ReLU activations, while the output layers of
actor networks used tanh. The action space was re-scaled to the fit the tanh
range of [-1, 1], and to prevent vanishing gradients, the preactivations of
the actor output layers were penalized by the square of their magnitude with a
coefficient of 0.001.

The DDPG algorithm was run in parallel using multiple message passing
interface (MPI)-based workers. Network parameters and normalization statistics
were averaged across workers during update steps. The actor policy, $\pi_{e}$,
$\pi_{r}$, or $\pi_{c}$ that each worker used during experience gathering was
set as a hyperparameter. All workers used $\pi_{r}$ during performance
testing. Different worker amounts were used depending on the difficulty of the
task:

Table 3: Parallel Worker Amounts by Task Task | Number of MPI Workers
---|---
Stack 2, Sparse Rewards | 8
Stack 3, Sparse Rewards | 32
Stack 4, Sparse Rewards | 32
Stack 2, Incremental Rewards | 8
Stack 3, Incremental Rewards | 8
Stack 4, Incremental Rewards | 32

The following hyperparameters were used in our experiments:

Table 4: Hyperparameters for Block Stacking Tasks Hyperparameter | Value
---|---
Optimizer | Adam ([23])
$n_{r}$ Learning Rate | 0.001
$n_{r}$ L2 Regularization Coefficient | 0
$\pi_{r}$ Learning Rate | 0.001
Target Exploit Actor-Critic Polyak-averaging Coefficient | 0.001
$n_{e}$ Learning Rate | 0.001
$n_{e}$ L2 Regularization Coefficient | 0.01
$\pi_{e}$ Learning Rate | 0.001
Target Explore Actor-Critic Polyak-averaging Coefficient: | 0.05
$\pi_{c}$ Learning Rate | 0.001
$\pi_{c}$ Explore vs Exploit Critic Weighting | $0.5,0.5$
$D$ Learning Rate | 0.007
Episode Time Horizon | $50*\text{num blocks}$
$\gamma$ | $1-1/{\text{episode time horizon}}$
MPI Worker Replay Buffer Size |  $10^{6}$ transitions
Parameter Space Noise $\sigma$ Target | 0.1
Guassian Action Noise $\sigma$ | 0.04
Traditional HER Augmentation Probability (when used) | 0.8
Multi-Criteria HER Augmentation Probability (when used) | 0.8
Cycles per Epoch | 50
Experience Gathering Episodes per Cycle | 8 (per MPI worker)
Training Batches per cycle | 8
Network Update Batch Size | 1024 transitions (per MPI worker)
Test Episode Rollouts Per Epoch | 50


# Determination of geopotential by hydrogen masers based on global navigation
satellite system

###### Abstract.

According to the general relativity theory, a clock runs faster at a position
with higher geopotential. Thus, the geopotential difference can be determined
by comparing the clock’s vibration frequencies. Here we report on the
experiments to determine the geopotential difference between two remote sites
based on precise point positioning time-frequency transfer technique. The
experiments include the remote clock comparison and the local clock comparison
using two CH1-95 active hydrogen masers acted as the reference clocks of
global navigation satellite system time-frequency receivers. Considering the
local clock comparison as calibration, the determined geopotential difference
between two sites by our experiments is 12,142.3 (112.4) m2/s2, quite close to
the value 12,153.3 (2.3) m2/s2 computed by the EIGEN-6C4 model. Results show
that the proposed approach here for determining geopotential difference is
feasible, operable, and promising.

Lei Wang1, Wei Xu1,2, Li-Hong Li1, Peng-Fei Zhang1, An Ning1, Rui Xu1, Peng
Cheng1 Zi-Yu Shen3, Wen-Bin Shen1,4,*
1 Time and frequency Geodesy Center, School of Geodesy and Geomatics, Wuhan
University, Wuhan 430079, China.
2 School of Geographic Information and Tourism, Chuzhou University, Chuzhou
239000, China.
3 School of Resource, Environmental Science and Engineering, Hubei University
of Science and Technology, Xianning 437100, China.
4 State Key Laboratory of Information Engineering in Surveying, Mapping and
Remote Sensing, Wuhan University, Wuhan 430079, China.
E-mail<EMAIL_ADDRESS>

## Introduction

The general relativity theory (GRT)[1] states that a clock at a position with
higher geopotential ticks quicker than the identical clock at a position with
lower geopotential[2]. Previous studies have shown that one can use the
gravity frequency shift approach[3] to determine the geopotential difference
between arbitrary two sites where precise clocks are set by comparing the
clock’s vibration frequencies[4, 5, 6, 7]. This approach provides advantages
in some aspects to overcome the drawbacks of the conventional approach that
needs combining leveling with additional gravimetry[8], which is limited by
the error accumulation with the increase of the length of the measurement-line
and the impediment of mountains, rivers, and oceans when two points are
located at separated islands or continents[3]. However, actual application of
the new approach requires high-performance clocks with a fractional frequency
accuracy of 1×10-17 to 1×10-18 at least, corresponding to a resolution of
about 1 m2/s2 to 0.1 m2/s2 in geopotential[7].

With quick development of science and technology of atomic clocks, the
frequency instability for state-of-the-art optical clocks has reached 10-19
level[9, 10, 11], and the systematic uncertainties has achieved 10-18 and
below 10-18 level[10, 12, 13, 14, 15], corresponding to an error equivalent to
less than 1 s over the lifetime of the Universe[16]. Several groups in the
world completed relevant experiments to determine the geopotential difference
using optical clocks based on the optical fibre frequency transfer (OFFT)
approach[4, 5, 7]. These experiments have demonstrated that the results by
adopting the OFFT approach are consistent well with the independent
measurements by the conventional approach, which indicates that OFFT is
prospective and could be practically applied to determining the geopotential.
By OFFT approach, it is necessary to compare two clocks located at separated
sites connected with optical fibres. If the two sites are far away or located
at different continents separated by oceans, it is inconvenient or difficult
to connect clocks with optical fibres. Additionally, it is almost unrealistic
to construct an optical fibre network to connect two arbitrary sites due to
the high expense, especially anywhere in a country with a vast territory.
Therefore, the satellite-based time-frequency transfer plays significant
role[17, 18, 19], including global navigation satellite system (GNSS),
communication satellites, and other satellite constellations that thanks to
the fast development of space science and technology.

In recent decades, the GNSS techniques have been applied extensively to high-
accuracy time-frequency transfer since the International GNSS Service (IGS)
has provided precise GNSS orbit and clock products[20]. With the rapid
development of the precise point positioning (PPP) technique in recent
years[21, 22, 23], comparing with the PPP technique solving the carrier phase
ambiguities as floating numbers that limits the performance of the time-
frequency transfer to the low 10-16, the PPP ambiguity resolution (PPP-AR)
technique has been proposed to recover the integer nature of the carrier phase
ambiguities[24]. According to the latest results, the PPP-AR technique can
provide time-frequency transfer instability with an improved long-term
performance of order 7×10-16/$T$ for clock comparison, where $T$ is the
duration in days of continuous phase measurements, thus reaching a sub 10-16
level after one week of averaging[25] and can be applied for determination of
geopotential. Several different kinds of experiments related to satellite-
based time-frequency transfer techniques have been carried out to determine
geopotential difference, including for instance, GNSS common-view (CV)[26],
two-way satellite time and frequency transfer (TWSTFT)[27]. Here, we report
new experiments for determining the geopotential difference between two sites
using two hydrogen masers based on the PPP time-frequency transfer (PTFT)
using PPP-AR technique.

## Experiments

The experiments are carried out between Jiugongshan Time-Frequency Station
(JTFS) in Xianning and Luojiashan Time-Frequency Station (LTFS) in Wuhan,
China. The straight-line distance between the two stations is approximately
129 km, and the height difference between them is about 1,245 m. We employed
two parallel independent sets of instruments, denoted as equipment set I and
equipment set II. Both equipment sets have the same specifications and
performances to avoid or weaken the system errors induced by the instruments.
Each equipment set includes a CH1-95 active hydrogen maser (AHM) (C1/C2), a
GNSS time-frequency receiver, a choke-ring antenna and relevant cables
(FIGURE. 1). The hydrogen masers were manufactured by the factory KVARZ in
Russia, and the production (nominal) frequency instability of either hydrogen
maser is less than 5×10-16 @ day. After being transported from KVARZ to
Beijing in China, the performances of the two hydrogen masers, officially
inspected by the National Institute of Metrology, China, are consistent with
the nominal ones (FIGURE. 2e). The consistency demonstrates that the
transportation process has very little or even no effect on the performances
of the hydrogen masers, which provides strong supports for the experiments.
The hydrogen masers provided time and frequency standards with one pulse per
second (1PPS) and 10 MHz for GNSS time-frequency receivers, and the receivers
can obtain the pseudorange and carrier phase observations through the choke-
ring antennas that were applied to reducing the multipath effects while
receiving navigation signals from visible GNSS satellites. These equipment
sets were installed at JTFS or/and LTFS, respectively, to form the GNSS
observation stations for time-frequency transfer. For convenience, when
equipment set N (N=1, 2) was installed at JTFS/LTFS, the GNSS station is
denoted as J00N/L00N.

We implemented both remote clock comparison (RCC) and local clock comparison
(LCC) denoted as AB experiment and AA experiment, respectively. The schematic
outline of experiments is shown in FIGURE. 1. In both AB and AA experiments,
the CH1-95 AHM and GNSS time-frequency receiver worked in well-conditioned
indoor with relative humidity of 60$\%$ and stable temperature of about 25 °C,
the main environment requirements for the CH1-95 AHM working stably, to
largely keep from potential environment disturbances, and the choke-ring
antenna was set on the roof for good observation sight.

Figure 1. Schematic of the experimental setups and process. Two sets of
instruments, each including a CH1-95 AHM, denoted as C1 or C2, a GNSS time-
frequency receiver, a choke-ring antenna and relevant cables to form a GNSS
observation station for time-frequency transfer. Control systems are used to
avoid potential environmental influences. (1) In AB experiment, C1 and C2 are
used as external clocks for GNSS stations J001 at JTFS and L002 at LTFS,
respectively. Clock comparison is conducted by PTFT using GNSS pseudorange and
carrier phase observations. (2) In AA experiment, the whole set of instruments
at LTFS is transported to JTFS with an air-conditioned vehicle (along the red
arrow). C2 is placed next to C1 and forms a new GNSS station J002 at JTFS.

In AB experiment, the equipment set I was installed at JTFS and the equipment
set II was installed at LTFS to construct GNSS stations J001 and L002 to form
the remote time-frequency transfer link, and the observations lasted for 14
days. For calibration, we need AA experiment to determine the inherent
frequency offset between C1 and C2 without any adjustments. The equipment set
II was transported to JTFS from LTFS using an air-conditioned vehicle for
about four-hours journey, during which the C2 always ran normally powered by
an uninterruptible power supply (UPS). Then, the equipment set II was
installed at JTFS, where C2 was placed next to and at the same height as C1
(FIGURE. 1). A new GNSS station J002 at JTFS was constructed and formed the
local time-frequency transfer link with J001 that had been constructed, and
the observations also lasted for 14 days in AA experiment. The two hydrogen
masers ran continuously and freely without any artificial adjustments
throughout AB and AA experiments, which is a critical step to determine
geopotential difference using atomic clocks.

## Data Analysis

In AA experiment, at beginning the C2 ran less reliably due to the vibration
and shaking during the four-hours transportation, and about two days later it
ran with good performance. Here, we selected continuous and stable
observations data with quality control (enough visible satellites, less signal
disturbance, etc., see Supplementary Information (SI)) to resolve the clock
offset series of GNSS stations, and then determine the time difference series
between two clocks (see Methods). In AB and AA experiments, we selected the
GPS data segments from days of year (DOY) 59 to 69 and 76 to 85 in 2021,
respectively. The sampling interval of all GPS observations used in this study
is 30 s. Previous studies and experiments have shown that more visible
satellites and smaller time dilution of precision (TDOP) values can improve
the accuracy of the GNSS positioning and timing and thus improve the frequency
stability of time-frequency transfer link between two stations[28] . Each GNSS
station in AB or AA experiment can observe more than 7 GPS satellites every
day, and the TDOP values of all stations are less than 1.7 (see SI). These
demonstrate that all GNSS stations have excellent observation conditions.

Figure 2. Time and frequency difference series and frequency instability of
clock comparisons. a, Time difference series (blue dots) and the fitting line
(red line) in AB experiment from DOY 59 to 69 in 2021. Inset shows the time
difference series of 0.1 day in DOY 59, 2021 and green circles mark jumps. b,
Time difference series (blue dots) and the fitting line (red line) in AA
experiment, whose y-axis subtracting 400 ns shares the same y-axis with a.
Inset shows the time difference series of 0.1 day in DOY 76, 2021 and green
circles mark jumps. c, Frequency difference series (blue curve) and the
average value of frequency differences (red line) of the link in AB experiment
where the frequency drift has been removed. d, Frequency difference series
(blue curve) and the average value of frequency differences (red line) of the
link in AA experiment where the frequency drift has been removed. e, Frequency
instability of the CH1-95 AHM. The nominal frequency instability of either
hydrogen maser (red curve with dots) and the frequency instabilities inspected
authoritatively (green curve with triangles and blue curve with squares)
expressed as Allan deviation (ADEV) present excellent consistency. f,
Frequency instability of links in AB (red curve with dots) and AA (blue curve
with triangles) experiments expressed as ADEV. The error bars show the
uncertainty.

Based on the PTFT, the time difference between the hydrogen maser (as the
external clock for receiver) at the station of interest and the reference time
of precise products used in resolving is called receiver clock offset in PPP,
and it can be estimated as a parameter every epoch to form the clock offset
series (see Methods). Then the time difference series between two hydrogen
masers can be obtained by differencing the clock offset series of two stations
(see Methods) and consequently the time difference series can be transformed
to the frequency difference series (see SI). To determine the frequency shift
between C1 and C2 with high accuracy and stability, we applied the latest
open-source software, PRIDE PPP-AR[29], to processing the GPS observations
(see Methods) with excellent ambiguity resolution (see SI). Then, the time
difference series between C1 and C2 in AB and AA experiments were obtained
(FIGURE. 2a-b), from which we determined the change rate of time difference
after removing the quadratic term, equivalent to the frequency drift between
C1 and C2, by directly linear fitting with least-squares (see Methods). The
gravity frequency shift caused by geopotential difference between JTFS and
LTFS comes from the frequency shift in AB experiment subtracting the frequency
offset in AA experiment. For this resolution, some jumps occurred in time
difference series (such as the insets in FIGURE. 2a-b) can influence the
accuracy of the determined frequency shift. Hence we take the following
measures: 1) transform time difference series to frequency difference series
(see SI); 2) replace frequency difference jump with linearly interpolated
values determined by the two values before and after any interior jumps; 3)
remove the frequency drift; 4) compute the average value of frequency
differences (FIGURE. 2c-d). Thus, the average values of frequency differences
demonstrate respectively the frequency shift in AB experiment and the
frequency offset in AA experiment between C1 and C2, and the gravity frequency
shift between A and B can be obtained by removing the frequency offset in AA
experiment from the frequency shift in AB experiment. Based on the determined
gravity frequency shift, we may determine the geopotential difference between
JTFS and LTFS.

## Result

As shown in FIGURE. 2a-b, the time differences between C1 and C2 in AB and AA
experiments are reduced by about 235 ns for 11 days and about 97 ns for 10
days, respectively. The AA experiment demonstrates that there exists the
inherent frequency offset between C1 and C2, and therefore AA experiment is
meaningful and necessary as calibration for determining the geopotential
difference. As shown in FIGURE. 2e-f, the frequency instabilities of links in
AB and AA experiments expressed as ADEVs computed by Stable32 software
(http://www.stable32.com/) can respectively reach 9.2468 ×10-16 and 8.4241
×10-16 at the averaging time of 120000 s, which are roughly consistent with
the nominal frequency instability of the hydrogen masers (FIGURE. 2e). These
results demonstrate that the PTFT has taken full advantages of the
performances of CH1-95 AHM. The resolved frequency differences between C1 and
C2 in AB and AA experiments (FIGURE. 2c-d), derived from the average values of
frequency differences transformed from time difference series, are
-2.4729×10-13 and -1.1219×10-13, respectively. Thus, the gravity frequency
shift caused by the geopotential difference is -1.3510×10-13 with an
uncertainty of 1.2501×10-15 derived from the ADEVs of links in AB and AA
experiments based on the error propagation law. Hence, the determined
geopotential difference between A and B is 12142.3 (112.4) m2/s2 (the numbers
in the parentheses are the 1 $\sigma$ uncertainties referred to the
corresponding last digits of the quoted results[5]). Compared with the
corresponding model value 12153.3 (2.3) m2/s2 computed by the EIGEN-6C4
model[30] (see Methods), the deviation is 11.0 (112.4) m2/s2 (equivalent
approximately to 1.1 (11.5) m in height).

## Conclusion

Our experiments have shown an approach to directly determine the geopotential
difference between two remote sites based on the PTFT using active hydrogen
masers. Although the accuracy at present is not high enough compared with the
conventional approach, the experiments here have demonstrated a feasible,
convenient, and economic approach to determine geopotential difference between
arbitrary two points. The accuracy of the results is mainly influenced by the
performances of the hydrogen masers, the observation noises induced by GNSS
satellite signals, observation environments and instruments. Thanks to quick
development of time and frequency science and technology in recent years,
high-performance clocks (including portable optical clocks) with frequency
instability at 10-19 level have been successively generated[9, 10, 11]. The
free-space dissemination of time and frequency with 10-19 instability over 113
km has been reported that could lay the groundwork for satellite
time–frequency dissemination[31] and the fractional frequency measurement
uncertainty of 7.6 × 10-21 has been reached[16]. Hence, in the future the
conventional approach could be supplemented or even replaced by PTFT using
high-performance clocks and advanced techniques, greatly impacting future the
way of determining geopotential and promoting the developments of time-
frequency science and geodesy.

## Methods

### Geopotential difference determination based on gravity frequency shift.

Suppose two atomic clocks, and , have been a priori calibrated, and they are
placed at two points 1 and 2, respectively. After a standard time duration ,
accurate to the level of , the geopotential difference between these two
clocks can be expressed as[2, 3]:

(1) $\Delta W_{12}=\left(W_{2}-W_{1}\right)\approx-\frac{\Delta
t_{12}}{T}c^{2}\approx-\frac{\Delta f_{12}}{f_{0}}c^{2}$

where $W_{1}$ and $W_{2}$ denote the geopotentials at points 1 and 2,
respectively; $\Delta t_{12}=t_{2}-t_{1}$ is the time difference between C1
and C2 , where $t_{1}$ and $t_{2}$ are respectively the time durations
recorded by clocks C1 and C2 after a time duration $T$; $\frac{\Delta
f_{12}}{f_{0}}$ is the (relative) gravity frequency shift between C1 and C2,
where $\Delta f_{12}=f_{2}-f_{1}$ and $f_{0}$ is the nominal frequency of the
atomic clock; $c$ is the speed of light in the vacuum.

### Geopotential difference determination by global gravity field model
EIGEN-6C4.

The geodetic coordinates of the GNSS station $i\left(i=1,2\right)$ can be
resolved along with clock offset while estimating parameters. However, that is
the coordinates at the center of GNSS receiver antenna, where the geopotential
$W_{ant}$ and gravity $g_{ant}$ at antenna can be computed online applying
global gravity field model EIGEN-6C4 from ICGEM website (https://icgem.gfz-
potsdam.de/). The vertical height difference between GNSS receiver antenna and
the maser cavity of CH1-95 AHM, $\Delta H=H_{ant}-H_{m}$, is measured by laser
range finder, where $H_{ant}$ and $H_{m}$ are the vertical heights above the
geoid of antenna and maser cavity, respectively. Then, the geopotential at the
maser cavity of CH1-95 AHM C1/C2 can be obtained as $W^{i}_{m}\approx
W^{i}_{ant}+g^{i}_{ant}\cdot\Delta H^{i}$, and the geopotential difference
between hydrogen masers is $\Delta W^{12}_{m}=W^{2}_{m}-W^{1}_{m}$.

### Precise point positioning time-frequency transfer (PTFT).

As the method of determining geopotential difference based on GRT described,
the key is to compare two atomic clock’s vibration frequencies. The PTFT can
be applied to completing this task. With this technique, the clock offset
between the atomic clock time $\left(t_{A}\right)$ and the reference time of
precise products $\left(t_{r}\right)$ can be estimated as a parameter using
GNSS pseudorange and carrier phase observations[22, 32]. The atomic clock can
provide time and frequency standards for the GNSS time-frequency receiver, and
$t_{r}$ is a very stable ensemble time scale released by data analysis center,
here used is IGS data canter of Wuhan university (http://www.igs.gnsswhu.cn/).
The clock offset between $t_{A}$ and $t_{r}$ can be expressed as
$t^{P}_{A}-t_{r}$ for any GNSS station, where $t^{P}_{A}$ denotes the local
time recorded by the atomic clock at point $P$. Thus, the time difference
$\Delta t^{P}_{Q}$between two atomic clocks at any two points $P$ and $Q$ can
be computed by differencing $t^{P}_{A}-t_{r}$ and $t^{Q}_{A}-t_{r}$, expressed
as[32]:

(2) $\Delta
t_{PQ}=t^{Q}_{A}-t^{P}_{A}=\left(t^{Q}_{A}-t_{r}\right)-\left(t^{P}_{A}-t_{r}\right)$

After a standard time duration $T$, the change rate of time difference,
$\Delta\frac{t_{PQ}}{T}$, can be determined by linear fitting with least
square. According to equation (2), $t_{r}$ is used only as a reference or
“bridge”. Therefore, arbitrary two clocks can be compared based on the PTFT
technique, no matter where the clocks are located. In addition, the PTFT is
easy to put into operation that is not required for synchronous observations
and can provide excellent short-term and medium-term frequency stability.
Therefore, this technique can be prospectively adopted and will be widely
applied for determining geopotential difference in the future.

## References

  * [1] Einstein, A. Die Feldgleichungen der Gravitation. Sitzungsberichte der Königlich Preu$\beta$ischen Akad. der Wissenschaften 1, 844–847 (1915).
  * [2] Bjerhammar, A. On a relativistic geodesy. Bull. Géodésique 59, 207–220 (1985).
  * [3] Shen, W., Ning, J., Liu, J., Li, J. & Chao, D. Determination of the geopotential and orthometric height based on frequency shift equation. Nat. Sci. 03, 388–396 (2011).
  * [4] Takamoto, M. et al. Test of general relativity by a pair of transportable optical lattice clocks. Nat. Photonics 14, 411–415 (2020).
  * [5] Takano, T. et al. Geopotential measurements with synchronously linked optical lattice clocks. Nat. Photonics 10, 662–666 (2016).
  * [6] Lisdat, C. et al. A clock network for geodesy and fundamental science. Nat. Commun. 7, 1–7 (2016).
  * [7] Grotti, J. et al. Geodesy and metrology with a transportable optical clock. Nat. Phys. 14, 437–441 (2018).
  * [8] Hofmann-Wellenhof, B. & Moritz, H. Physical Geodesy. (Physical Geodesy, by B. Hofmann-Wellenhof and H. Moritz. Approx. 500 p. 100 illus. 3-211-23584-1. Berlin: Springer, 2005., 2005).
  * [9] Campbell, S. L. et al. A Fermi-degenerate three-dimensional optical lattice clock. Science (80-. ). 358, 90–94 (2017).
  * [10] McGrew, W. F. et al. Atomic clock performance enabling geodesy below the centimetre level. Nature 564, 87–90 (2018).
  * [11] Oelker, E. et al. Demonstration of 4.8 × 10-17 stability at 1 s for two independent optical clocks. Nat. Photonics 13, 714–719 (2019).
  * [12] Brewer, S. M. et al. 27Al+ Quantum-Logic clock with a systematic uncertainty below 10-18. Phys. Rev. Lett. 123, 33201 (2019).
  * [13] Bloom, B. J. et al. An optical lattice clock with accuracy and stability at the 10-18 level. Nature 506, 71–75 (2014).
  * [14] Huang, Y. et al. Liquid-Nitrogen-Cooled Ca+ optical clock with systematic uncertainty of 3×10-18. Phys. Rev. Appl. 17, 1 (2022).
  * [15] Nicholson, T. L. et al. Systematic evaluation of an atomic clock at 2 × 10-18 total uncertainty. Nat. Commun. 6, (2015).
  * [16] Bothwell, T. et al. Resolving the gravitational redshift across a millimetre-scale atomic sample. 602, (2022).
  * [17] Allan, D. W. et al. Accuracy of international time and frequency comparisons via global positioning system satellites in common-view. IEEE Trans. Instrum. Meas. 34, 118–125 (1985).
  * [18] Fujieda, M. et al. Carrier-phase-based two-way satellite time and frequency transfer. IEEE Trans. Ultrason. Ferroelectr. Freq. Control 59, 2625–2630 (2012).
  * [19] Petit, G. The TAIPPP pilot experiment. 2009 IEEE Int. Freq. Control Symp. Jt. with 22nd Eur. Freq. Time Forum 116–119 (2009) doi:10.1109/FREQ.2009.5168153.
  * [20] Kouba, J. A Guide to using International GNSS Service ( IGS ) products. Geod. Surv. Div. Nat. Resour. Canada Ottawa 6, 34 (2009).
  * [21] Geng, J., Chen, X., Pan, Y. & Zhao, Q. A modified phase clock/bias model to improve PPP ambiguity resolution at Wuhan University. J. Geod. 93, 2053–2067 (2019).
  * [22] Laurichesse, D., Mercier, F., Berthias, J. P., Broca, P. & Cerri, L. Integer ambiguity resolution on undifferenced GPS phase measurements and its application to PPP and satellite precise orbit determination. Navig. J. Inst. Navig. 56, 135–149 (2009).
  * [23] Ge, M., Gendt, G., Rothacher, M., Shi, C. & Liu, J. Resolution of GPS carrier-phase ambiguities in precise point positioning (PPP) with daily observations. J. Geod. 82, 389–399 (2008).
  * [24] Petit, G. et al. 1 × 10-16 frequency transfer by GPS PPP with integer ambiguity resolution. Metrologia 52, 301–309 (2015).
  * [25] Petit, G. Sub-10–16 accuracy GNSS frequency transfer with IPPP. GPS Solut. 25, 1–9 (2021).
  * [26] Kopeikin, S. M. et al. Chronometric measurement of orthometric height differences by means of atomic clocks. Gravit. Cosmol. 22, 234–244 (2016).
  * [27] Cheng, P. et al. Measuring height difference using two-way satellite time and frequency transfer. Remote Sens. 14, 1–15 (2022).
  * [28] Su, K. & Jin, S. Triple-frequency carrier phase precise time and frequency transfer models for BDS-3. GPS Solut. 23, (2019).
  * [29] Geng, J. et al. PRIDE PPP-AR: an open-source software for GPS PPP ambiguity resolution. GPS Solut. 23, 1–10 (2019).
  * [30] Ince, E. S. et al. ICGEM-15 years of successful collection and distribution of global gravitational models, associated services and future plans. Earth System Science Data (2019) doi:10.5194/essd-2019-17.
  * [31] Physics, Q., National, H., Astronomical, X., Standards, F. P. & Time, N. Free-space dissemination of time and frequency with 10-19 instability over 113 km. (2022) doi:10.1038/s41586-022-05228-5.
  * [32] Petit, G. & Arias, E. F. Use of IGS products in TAI applications. J. Geod. 83, 327–334 (2009).
  * [33] Landskron, D. & Böhm, J. VMF3/GPT3: refined discrete and empirical troposphere mapping functions. J. Geod. 92, 349–360 (2018).
  * [34] Pan, L., Xiaohong, Z. & Fei, G. Ambiguity resolved precise point positioning with GPS and BeiDou. J. Geod. 91, 25–40 (2017).
  * [35] Barnes, J. A. et al. Characterization of frequency stability. IEEE Trans. Instrum. Meas. IM–20, 105–120 (1971).
  * [36] Rutman, J. Characterization of phase and frequency instabilities in precision frequency sources: fifteen years of progress. Proc. IEEE 66, 1048–1075 (1978).

## Supplementary information

### Visibility of GPS satellites.

FIGURE. S1 presents the daily average values of the number of visible GPS
satellites and TDOP values of GNSS stations in AB and AA experiments. The
visibility of the GPS satellites at station L002 in AB experiment is worse
than that at station J002 in AA experiment. Meanwhile, the TDOP values at L002
are larger than that at J002, which means that the time-frequency transfer
accuracy of L002 is worse than J002. The visibility of the GPS satellites and
TDOP values at J001 and J002 are nearly the same, respectively. Hence, the
frequency stability of the link in AA experiment is better than that in AB
experiment, which is demonstrated in FIGURE. 2f.

Figure S1. Daily average values of the number of visible satellites and TDOP
values of GNSS stations. a, Number of visible GPS satellites (histograms) and
TDOP values (curves with points) of GNSS stations J001 (red) at JTFS and L002
(blue) at LTFS in AB experiment from DOY 59 to 69 in 2021. b, Number of
visible GPS satellites (histograms) and TDOP values (curves with points) of
GNSS stations J001 (red) and J002 (green) at JTFS in AA experiment from DOY 76
to 85 in 2021.

### Precise point positioning ambiguity resolution (PPP-AR) technique.

In our experiments, to process the GPS pseudorange and carrier phase
observations reliably and accurately using the PPP-AR technique, we applied
the latest open-source software, PRIDE PPP-AR[29]. This software has been
developed by the PRIDE Lab at the GNSS Research Center of Wuhan University. It
can provide float solutions with wide-lane and narrow-lane ambiguity estimates
and integer ambiguity resolution recovering the integer nature of single-
station ambiguities by using the phase clock/bias products. PRIDE PPP-AR has
been demonstrated that the ambiguity-fixed positions outperform those by float
solutions. According to the PPP parameters estimation, the accuracy of the
ambiguity-fixed clock offset can also be improved and performed well in time-
frequency transfer links. The strategies and options of GPS observation data
processing using PRIDE PPP-AR are presented in TABLE S1.

Table S1. Summary of the GPS observation data processing strategies and options for PRIDE PPP-AR Items |  Strategies or options
---|---
Solutions |  PPP-AR
Observables |  Ionosphere-free combination of GPS pseudorange and carrier phase
Cutoff elevation |  7°
Sampling interval |  30 seconds
Strict editing |  YES
Zenith tropospheric delay |  PWC: 60 minutes
Horizontal tropospheric gradients |  PWC: 720 minutes
Troposphere mapping function |  Vienna Mapping Function 3 (VMF3) [33]
Correcting second-order ionospheric delays |  YES
Ambiguity fixing |  ROUNDING
Ambiguity duration |  600 seconds
Cutoff elevation for PPP-AR |  15°
Widelane decision |  Deviation: 0.2; Sigma: 0.15 (Cycle)
Narrowlane decision |  Deviation: 0.15; Sigma: 0.15 (Cycle)
Figure S2. Fixing rates of ambiguity resolution for GNSS stations. In AB
experiment (RCC), GNSS stations are J001 (red) at JTFS and L002 (blue) at
LTFS. In AA experiment (LCC), GNSS stations are J001 (red) and J002 (green) at
JTFS.

### Fixing rates of ambiguity resolution.

The fixing rate of ambiguity resolution is a useful indicator to evaluate the
performance of PPP-AR, expressed as[34]

(S1) $FR=\frac{N_{fixed}}{N_{total}}\times 100\%$

where $N_{fixed}$ denotes the number of both wide lane and narrow lane
ambiguities-fixed, $N_{total}$ represents the number of total ambiguities-
float. FIGURE. S2 shows the fixing rates of ambiguity resolution in our
experiments. The fixing rates of all stations are larger than 93.4% and the
fixing rate of J001 at JTFS in AB experiment can reach 97.6%. These
demonstrate that the carrier phase ambiguities in our experiments have been
fixed effectively and the results resolved by PPP-AR are reliable and
accurate.

### Transformation from time difference to frequency difference.

Consider a signal generator whose instantaneous output $V(t)$ can be written
as[35]

(S2) $V(t)=\left[V_{0}+\varepsilon(t)\right]sin\left[2\pi
f_{0}t+\varphi(t)\right]$

where $V_{0}$ and $f_{0}$ are the nominal amplitude and frequency, $t$ is the
time variable, $\varepsilon(t)$ and $\varphi(t)$ are the instantaneous
fluctuations of amplitude and phase, respectively. The instantaneous
fractional frequency deviation from the nominal frequency, $y(t)$, is defined
as[35]

(S3) $y(t)=\frac{f(t)-f_{0}}{f_{0}}=\frac{\dot{\varphi}(t)}{2\pi f_{0}}$

where $\dot{\varphi}(t)=\frac{d\varphi(t)}{dt}$ is sufficiently small for $t$.
And the instantaneous fractional phase deviation is defined as[35]

(S4) $x(t)=\frac{\varphi(t)}{2\pi f_{0}}$

From equations (S3) and (S4) it reads

(S5) $y(t)=\frac{dx(t)}{dt}$

Combining equations (S2) and (S4), one obtains

(S6) $V(t)=\left[V_{0}+\varepsilon(t)\right]sin\left[2\pi
f_{0}\left(t+x(t)\right)\right]$

where $x(t)$ presents the reading deviation of a clock running from the signal
given by equation (S2), and here we call it time deviation.

Since any frequency-measurement technique does involve a finite time interval
over which the measurement is performed, the instantaneous fractional
frequency deviation is not an observable. Therefore, the average value of
fractional frequency deviation $\bar{y}_{\tau_{k}+\tau}$over a time interval
$\tau$ beginning at $\tau_{k}$ provides a more useful quantity directly
related to an experimental result, and it can be expressed as[36]

(S7)
$\bar{y}_{\tau_{k}+\tau}=\frac{1}{\tau}\int^{\tau_{k}+\tau}_{\tau_{k}}y(t)dt=\frac{1}{\tau}\int^{\tau_{k}+\tau}_{\tau_{k}}dx(t)=\frac{x(\tau_{k}+\tau)-x(\tau_{k})}{\tau}$

where the average value of fractional frequency deviation denotes the change
rate of time deviation from $\tau_{k}$ to $\tau_{k}+\tau$.

To compare two clocks, we determine the frequency difference between two
clocks by subtracting one fractional frequency deviation from another one. In
the time interval $\left[0,n\tau\right]$, the average value of frequency
difference between two clocks can be expressed as:

(S8)
$\bar{y}^{12}_{dif}=\frac{1}{n}\sum^{n}_{i=1}\left[\bar{y}^{2}_{i\tau}-\bar{y}^{1}_{i\tau}\right]=\frac{{\left[x^{2}(n\tau)-x^{1}(n\tau)\right]}-{\left[x^{2}(0)-x^{1}(0)\right]}}{n\tau}$