# Artifacts and Errors in
Cross-Spectrum Phase Noise Measurements

Yannick Gruson∃, Adrian Rus⋆, Ulrich L. Rohde◇,
Alexander Roth□, and Enrico Rubiola${}^{\exists\,\nabla\,\Re}$
$\exists$ FEMTO-ST Institute, UBFC and CNRS, Besançon, France.
$\star$ Broadhurst, Bucharest, Romania. Call sign YO3HHZ.
$\Diamond$ Synergy Microwave Corp., Paterson, NJ, USA. Call sign N1UL.
$\Box$ Rohde & Schwarz, München, Germany.
$\nabla$ Istituto Nazionale di Ricerca Metrologica INRiM, Torino, Italy
$\Re$ Reference author
http://rubiola.org

(March 25, 2020)

###### Abstract

This article deals with the erratic and inconsistent phase-noise spectra often
seen in low-noise oscillators, whose floor is of the order of $-180$ dBc/Hz or
less. Such oscillators are generally measured with two-channel instruments
based on averaging two simultaneous and statistically independent measures.
Our new method consists of inserting a dissipative attenuator between the
oscillator under test and the phase-noise analyzer. The thermal noise of the
attenuator introduces a controlled amount of phase noise. We compare the phase
noise floor to the theoretical expectation with different values of the
attenuation in small steps. The analysis reveals a negative bias
(underestimation of phase noise) due to the thermal energy of the internal
power splitter at the instrument input, and an uncertainty due to crosstalk
between the two channels. In not-so-rare unfortunate cases, the bias results
in a negative phase-noise spectrum, which is an obvious nonsense. Similar
results are observed separately in three labs with instruments from the two
major brands. We give experimental evidence, full theory, and suggestions to
mitigate the problem. Our multiple-attenuators method provides quantitative
information about the correlation phenomena inside the instrument.

##  1 Introduction

Modern analyzers measure the phase noise (PN, or PM noise) by correlation and
averaging on the simultaneous measurement of the oscillator under test (DUT)
with two separate channels, each consisting on a phase detector and a
frequency reference. The DUT noise is extracted after rejecting the background
noise of the instrument, assuming that the two channels are statistically
independent. After the seminal paper [1], and the early application shown in
[2], this choice is adopted by virtually all manufacturers (Table 1). The dual
channel scheme comes in two flavors, with one or two reference oscillators. We
focus on the latter because it enables the noise rejection of the reference
oscillators, and also of the frequency synthesizers which may be interposed
between reference oscillators and phase detectors.

The correlation-and-averaging process rejects the single-channel noise
proportionally to $1/\sqrt{m}$, where $m$ is the number of averaged spectra,
that is, 5 dB per factor-of-ten. Nowadays, digital electronics provides a high
computing power and memory size for cheap, as compared to the cost and to the
complexity of RF and microwave technology. Thus, the theoretical rejection can
exceed 30 dB if the experimentalist accepts the long measurement time it
takes, ultimately limited by the time-frequency indetermination theorem.
However, such rejection cannot be achieved in practice because of fundamental
phenomena and artifacts. The thermal energy in the input power splitter [3, 4,
5, 6, 7] and impedance matching [8] first caught the attention of the
scientific community. These and other problems were addressed in three
international workshops [9, 10, 11].

Most practitioners, naively, believe that a noise analyzer always _over-
estimates_ the DUT noise because it _adds_ its own background noise. This is
not true in the case of the two-channel instruments because the cross spectrum
is the frequency-domain equivalent of the covariance. The correlation between
channels introduces systematic errors and artifacts, which can be positive or
negative. The consequence is that there is no a-priori rule to state whether
the instrument over-estimates or under-estimates the DUT noise. A problem is
that the noise rejection due to averaging, usually calculated and displayed
together with the phase noise, does not account for artifacts and systematic
errors. Another problem is that the instruments display the absolute value of
the cross spectrum, giving no warning about negative outcomes. The combination
of these facts originates erratic and misleading results.

We propose an experiment (Fig. 1) that reveals the presence systematic errors
due to unwanted correlated terms. We focus on the 100 MHz OCXOs because this
type of oscillator exhibits the lowest white PM noise floor. However trivial
the experiment may seem, nothing even broadly similar has been reported in the
literature. We provide all the details related to two specific cases, together
with the full theoretical interpretation.

Table 1: Dual-Channel Phase Noise Analyzers Brand | Type or Series
---|---
AnaPico | APPH series
Berkeley Nucleonics Corp. | 7000 series
Holzworth | HA7062 series
Jackson Labs Technologies | PhaseStation 53100A
Keysight Technologies | E5500 series
Microsemi Corporation | 3120A / 5120A
NoiseXT | DCNTS / NXA
OEwaves | HI-Q TMS
Rohde & Schwarz | FSWP series
Figure 1: Block diagram of the experiment.

##  2 Phase Noise and Thermal Energy

Let us start with a review of key facts, based on References [3, 4, 5, 6, 7].
The phase noise is described in terms the power spectral density (PSD) of the
random phase $\varphi(t)$, and denoted with $S_{\varphi}(f)$. A model that is
found useful to describe oscillators and components is the polynomial law

$S_{\varphi}(f)=\sum_{n\leq-4}^{0}\mathsf{b}_{n}f^{n}~{}~{}~{}\mathrm{[rad^{2}/Hz]}~{},$
(1)

where the term $\mathsf{b}_{0}$ is the white PM noise, $\mathsf{b}_{-1}/f$ is
the flicker PM noise, $\mathsf{b}_{-2}/f^{2}$ is the white FM noise,
$\mathsf{b}_{-3}/f^{3}$ is the flicker FM noise, $\mathsf{b}_{-4}/f^{4}$ is
the frequency random walk, and other terms can be added. The quantity
$\mathscr{L}(f)$, most often used by the manufacturers, is defined as
$\mathscr{L}(f)=(1/2)S_{\varphi}(f)$ [12]. As a matter of fact, white phase
noise is mostly of additive origin. Accordingly, it can be written as

$\mathsf{b}_{0}=\frac{N}{P}~{},$ (2)

where $P$ is the carrier power, and $N$ is the power spectral density (PSD) of
the RF noise. In this context, we prefer the unit W/Hz to J. By analogy with
the PSD $N=kT$ of the thermal noise, we associate to $\mathsf{b}_{0}$ the
equivalent temperature

$T=\frac{P\,\mathsf{b}_{0}}{k},$ (3)

where $k=1.380649{\times}10^{-23}$ J/K (exact) is the Boltzmann constant.

###  2.1 The Effect of the Attenuator

Physical insight suggests that the dissipative attenuator can only degrade the
signal-to-noise ratio (SNR), which results in increased PM noise. Focusing on
the white noise at the attenuator in and out, we use the subscripts “$i$” and
“$o$” dropping the subscript 0. For example, $\mathsf{b}_{i}$ stands for
$\mathsf{b}_{0\,i}$, and $\mathsf{b}_{o}$ for $\mathsf{b}_{0\,o}$. Thus (2)
rewrites as $\mathsf{b}_{i}=N_{i}/P_{i}$, or $\mathsf{b}_{o}=N_{o}/P_{o}$.
Assuming that everything is matched to the characteristic impedance $R_{0}$,
the RF white noise at the attenuator output is

$N_{o}=kT_{i}A^{2}+kT_{a}\left(1-A^{2}\right)~{},$ (4)

where $kT_{i}$ is the input noise, $A$ is the voltage gain of the attenuator,
$A^{2}<1$, and $T_{a}$ is the temperature of the attenuator. The term
$kT_{i}A^{2}$ means that the input noise $kT_{i}$ is attenuated by the factor
$A^{2}$, like any signal. The term $kT_{a}(1-A^{2})$ is the thermal noise
added by the attenuator. This is obvious if one replaces the oscillator with a
resistive load $R_{0}$ at the temperature $T_{a}$. In this condition the
output is equivalent to a resistor $R_{0}$ at the temperature $T_{a}$. Thus,
the total output noise is $kT_{a}$, independent of $A$. Equation (4) is well
known in radio astronomy, where it finds application in the estimation of the
effect of losses in the antenna and in the line between antenna and receiver
[13, Sec. 7-2b (Noise Temperature of an Attenuator)], and in the calibration
of the receiver [14, Sec. 4.2.4 (Receiver Calibration)].

After (2)–(4), the white PM noise at the attenuator output is

$\mathsf{b}_{o}=\frac{kT_{i}}{P_{i}}+\frac{kT_{a}\left(1-A^{2}\right)}{A^{2}P_{i}},$
(5)

which is obviously greater than $\mathsf{b}_{i}=kT_{i}/P_{i}$.

##  3 Inside the Dual-Channel Noise Analyzer

###  3.1 The Cross-Spectrum Estimator

The cross-spectrum estimator is a general tool. Here, it finds application to
two key signals inside the instrument: (i) the voltage at the outputs of the
power splitter, and (ii) the phase at the output of the detectors (Fig: 1).
Let us start with

$\displaystyle x(t)$ $\displaystyle=a(t)+c(t)$ (6) $\displaystyle y(t)$
$\displaystyle=b(t)+c(t)$ (7)

where $a(t)$ and $b(t)$ are the background noise of the channel $\mathcal{A}$
and $\mathcal{B}$. They are statistically independent, and have zero mean and
equal or similar variance. The signal $c(t)$ is the target, that is, the DUT
noise. Thus, the statistical properties of $c(t)$ are measured after averaging
out $a(t)$ and $b(t)$. As said, (6)-(7) apply to the RF signal or to the PM
noise, at choice. Thus, we can measure the RF spectrum [W/Hz], or the PM noise
spectrum [$\mathrm{rad^{2}/Hz}$]. The reader interested to know more about the
method should refer to [15], [16] and [2].

As a mathematical concept, the cross PSD is defined as the Fourier transform
of the autocorrelation function. However, under certain conditions which are
generally met in the case of physical signals digitized on a finite
acquisition time $T$, it can be evaluated using the Fourier transforms. Thus,
denoting the discrete Fourier transform with the uppercase letter, as in
$X(f)\leftrightarrow x(t)$, the one-sided cross PSD is evaluated as

$S_{yx}(f)=\frac{2}{T}\Bigl{[}Y(f)\,X^{\ast}(f)\Bigr{]}~{}.$ (8)

The denominator $T$ is the acquisition time, the superscript “$\ast$” means
complex conjugate, and the factor “2” accounts for energy conservation after
suppressing the negative frequencies. Equation (8) states a general fact, thus
it holds for one realization, for the average or for the expectation,
depending on the context. Dropping the frequency and expanding $X=A+C$ and
$Y=B+C$ we get

$S_{yx}=\frac{2}{T}\Bigl{(}B+C\Bigr{)}\>\Bigl{(}A^{\ast}+C^{\ast}\Bigr{)}~{}.$
(9)

The mathematical expectation $\mathbb{E}\left\\{S_{yx}\right\\}$ is

$\mathbb{E}\bigl{\\{}S_{yx}\bigr{\\}}=\frac{2}{T}\mathbb{E}\bigl{\\{}CC^{\ast}\bigr{\\}}=\mathbb{E}\bigl{\\{}S_{c}\bigr{\\}}$
(10)

because $\mathbb{E}\\{BA^{\ast}\\}=0$, $\mathbb{E}\\{BC^{\ast}\\}=0$, and
$\mathbb{E}\\{CA^{\ast}\\}=0$. All the useful information is in $CC^{\ast}$,
thus $S_{c}>0$. By contrast, all the background noise goes in $BA^{\ast}$,
$BC^{\ast}$ and $CA^{\ast}$, and under normal circumstances it is equally
distributed between real and imaginary part. It is therefore clear that the
optimum estimator is

$\displaystyle\widehat{S_{yx}}$
$\displaystyle=\Re\left\\{\left<S_{yx}\right>_{m}\right\\}~{}.$ (11)

This estimator has two important properties, (i) it is unbiased, and (ii) it
is the fastest because it takes in the smallest amount of background noise. A
problem with $\Re\left\\{\left<S_{yx}\right>_{m}\right\\}$ is that it is not
always positive before averaging out the background noise, or in the presence
of spurs. The negative outcomes cannot be plotted on a logarithmic scale (dB).
The FSWP [17, Equation (4)] uses the estimator

$\displaystyle\widehat{S_{yx}}$
$\displaystyle=\left|\left<S_{yx}\right>_{m}\right|~{}.$ (12)

Albeit the documentation provided by other manufacturers gives little
indication about the estimator, we believe that (12) is the most chosen
option. A reason is that it shows no negative values, thus it is always
suitable to be represented on a log scale. Another reason is that such
estimator is positively biased, and the bias decreases monotonically as $m$
increases. Thus, under normal circumstances
$\left|\right<S_{yx}\left>{}_{m}\right|$ converges to
$\mathbb{E}\left\\{S_{c}\right\\}$, after decreasing monotonically during the
measurement process. The estimator (12) matches the behavior we observe in the
regular use of the instruments, where no attenuator is inserted at the input.

Now we break the hypothesis of statistically independent channels, and we
introduce the disturbing signal $d(t)\leftrightarrow D(f)$, the same in the
two channels but for the sign, $\varsigma_{x}=\pm 1$ and $\varsigma_{y}=\pm
1$, as we did in [18]

$\displaystyle X$ $\displaystyle=A+C+\varsigma_{x}D$ (13) $\displaystyle Y$
$\displaystyle=B+C+\varsigma_{y}D~{},$ (14)

The signal $D$ is either correlated or anticorrelated. Introducing
$\varsigma=\varsigma_{x}\varsigma_{y}=\pm 1$, and expanding
$\mathbb{E}\left\\{S_{yx}\right\\}$ as above, we find

$\mathbb{E}\bigl{\\{}S_{yx}\bigr{\\}}=\mathbb{E}\bigl{\\{}S_{c}\bigr{\\}}+\varsigma\mathbb{E}\bigl{\\{}S_{d}\bigr{\\}}~{}.$
(15)

Thus, $\varsigma$ is the sign of the correlation coefficient, and the term
$\varsigma\mathbb{E}\bigl{\\{}S_{d}\bigr{\\}}$ is a systematic bias, positive
or negative. A substantially equivalent approach is found in [5, Section II],
which differs in the analysis of four separate cases (the presence or not of
the disturbing signal $d(t)$, and the sign of the correlation), instead of the
compact form (13)–(14) and (15).

The disturbing signal can be (i) the thermal energy in the input power
splitter (Sec. 3.2), (ii) the crosstalk between the two channels (Sec. 5.1),
and (iii) the AM noise pickup [19], or other effects not considered here.

###  3.2 Application to the Input Power Splitter

Figure 2: Signal and noise model of the most common power splitters. The
reactive power splitter (top) is free from dissipation, thus it has 3
insertion loss. The resistive power splitter (bottom) has 6 dB insertion loss.

Two types of power splitters are mostly used, shown on Fig 2. The loss-free
splitter is a 3 dB directional coupler terminated at one input (dark port).
The resistive splitter is a Y network which attenuates the input signal by 6
dB. Here, $x(t)$ and $y(t)$ are the RF voltages at the output of the power
splitter. Denoting with $T_{o}$ equivalent noise temperature at the power-
splitter input, and with $T_{s}$ the temperature of the splitter, trite
calculation shows that the correlated noise is

$\mathbb{E}\left\\{S_{yx}\right\\}=\frac{1}{2}k\Bigl{(}T_{o}-T_{s}\Bigr{)}$
(16)

for the 3-dB dissipation-free coupler. Interestingly, (16) is a classical
result from Johnson thermometry [20, 21], with well known application in
microwaves [22, 23].

Similarly, we find

$\mathbb{E}\left\\{S_{yx}\right\\}=\frac{1}{4}k\Bigl{(}T_{o}-T_{s}\Bigr{)}$
(17)

for the 6-dB resistive coupler. Deriving (16) and (17) from (15), $\varsigma$
does not need to appear explicitly because it always hold that $\varsigma=-1$.

Because the output power is $P_{o}/2$ for the 3-dB splitter and $P_{o}/4$ for
the 6-dB splitter, the output SNR is the same, and the white PM noise is

$\mathbb{E}\left\\{\mathsf{b}_{o}\right\\}=\frac{k(T_{o}-T_{s})}{P_{o}}~{}.$
(18)

Reference [6, Section IV] provides an extension to other less common types of
power splitter.

###  3.3 Hardware Architectures

The FSWP [24, 17] is based on the SDR (Software Defined Radio) technology
after down-converting the input to an appropriate IF (see [25] for a modern
treatise of SDR). The mixers are used in the linear region because linearity
prevents the AM noise from polluting the phase noise measurement. The use of
I-Q mixers enables to unwrap the phase, and to measure beyond the IF. Two
operating modes are used, depending on the Fourier frequency. Up to 1 MHz, the
input RF signal is down converted to 1.3 MHz. Beyond 1 MHz, the reference
synthesizers are set close to the input frequency, keeping the beat note below
10 Hz. In both cases, $I$ and $Q$ of the down-converted signal are digitized,
and the phase information is extracted in FPGA. The FSWP uses a 3-dB coupler
as the input power splitter (actually, three different couplers are switched,
for $<1$ GHz, $1{-}8$ GHz, and $8{-}50$ GHz).

The E5052B [26] is based on direct phase detection with double-balanced mixers
as the phase-to-voltage converters. The mixers are saturated at both inputs,
and driven with synchronous signals kept in quadrature. The mixer output is
digitized and processed. The power splitter is a Y resistive network.

In both cases, the signals $x(t)$ and $y(t)$ used to calculate
$\widehat{S_{yx}(f)}$ are the instantaneous phases at the detector output,
sampled and digitized. This is the case of all the spectra (Figs. 3, 5, 6 and
8) and the coefficient $\mathsf{b}_{0}$ (Fig. 4 and 7).

##  4 Experiments and Results

The experiment consists of the measurement of the white noise floor after
inserting various dissipative attenuators in the path from the oscillator
under test to the phase noise analyzer, as shown on Figure 1.

(A): all spectra overlapped

(B): Same spectra of (A), fitted with the model (see Section 5).

Figure 3: Phase noise of the Wenzel 501-04623E OCXO measured with the FSWP 26.
Figure 4: White noise floor $\mathsf{b}_{0}$ (dots) taken from Fig. 3-A,
compared to the “attenuator only” model based on (5). The “full model” plot
and the parameters are discussed in Section 5.

Two oscillators are tested, $\mathcal{A}$ a Wenzel 501-04623E, and
$\mathcal{B}$ a Wenzel 501-25900B “Golden Citrine,” both 100-MHz OCXOs
intended for the lowest-noise applications. The former dates more than 20
years ago. The latter is the top low-PM-noise oscillator by Wenzel. After
comparing to the spectra published on the web pages of several manufactures,
$\mathcal{B}$ is the OCXO that exhibits the lowest white noise we have found,
below $-190$ dBc/Hz [27].

The oscillator is clamped on a vibration-damping breadboard, of the same type
commonly seen in optical experiments. A 150-MHz low-pass filter (MiniCircuits
SLP-150) is inserted at the oscillator output. The attenuation is obtained by
stacking small-size SMA attenuators at the filter output, close to the
oscillator. The attenuators (Radiall brand) are intended for DC to 18 GHz. In
most of the tests, the phase-noise analyzer is a Rohde Schwarz FSWP 26 with
high-stability OCXO and cross-spectrum options. The phase-noise analyzers are
referenced to a T4Science Hydrogen maser, in turn monitored vs other masers of
the same type. The power is measured with a Rohde & Schwarz power meter, which
replaces temporarily the phase-noise analyzer before each measurement. The
attenuation is evaluated as the power ratio. All the experiments are done in a
Faraday cage with usual isolation transformer and EMI filters. Temperature and
humidity are stabilized to $22\pm 0.5\,\mathrm{{}^{\circ}C}$ and $50\,\%\pm
10\,\%$ by a PID control, which also guarantee a drift smaller than 0.2
K/hour. The environment control is probably overdone for PM noise
measurements, yet it helps to get conservative results.

Figure 5: Phase noise spectrum of the same oscillator of Figure 3, measured
with a Keysight E5052B phase noise analyzer.

(A): all spectra overlapped

(B): Same spectra of (A), fitted with the model (see Section 5).

Figure 6: Phase noise of a Wenzel 501-25900B “Golden Citrine” 100 MHz OCXO
measured in the same conditions and with the same instruments of Fig. 3.

Figure 3 shows the phase noise spectra of the oscillator $\mathcal{A}$,
observed with different values of the attenuation between 0 dB and 27 dB. The
experimental data (dots) on Fig. 4 are the white PM noise from Fig. 3,
averaged on a suitable region 2–3 decades wide. Surprisingly, the observed
floor does not match the “attenuator only” plot. The latter is calculated from
(5). Instead, the floor decreases monotonically from 0 dB to 15 dB
attenuation, and it increases monotonically beyond.

Measuring the oscillator $\mathcal{A}$ with a Keysight E5052B, we see that the
white PM noise decreases monotonically with the attenuation, attends a minimum
at 9 dB, and increases at higher attenuation (Fig. 5). We could not push the
attenuation beyond 15 dB because the carrier power falls below the minimum for
the E5052B.

The anomalously low white PM noise when an attenuator is introduced was first
observed by one of us (AR) in his radio amateur lab at home, measuring a
Wenzel 501-04538F 10 MHz OCXO with a FSWP 8.

Comparing Fig. 5 to Fig. 3-A, the calibration of the two instruments is
consistent within at most a small fraction of a dB. The flicker of frequency
is the same, $\mathsf{b}_{-3}=-74.5$ $\mathrm{dB\,rad^{2}Hz^{2}}$. Likewise,
the white noise floor at 0 dB attenuation is the same, $\mathsf{b}_{0}=-172.4$
$\mathrm{dB\,rad^{2}/Hz}$. The 2-dB discrepancy in the flicker PM noise is not
significant because the $\mathsf{b}_{-1}$ coefficient is hardly readable on
Fig. 5.

Figure 7: White noise floor $\mathsf{b}_{0}$ (dots) taken from Fig. 6-A,
compared to the “attenuator only” model based on (5). The “full model” plot
and the parameters are discussed in Section 5.

Figure 3-B shows the same plots of Fig. 3-A, just separated for better
readability. The most interesting fact is the appearance of dips at 1–1.5 kHz
for attenuation of $\geq 18$ dB.

Figures 6 and 7 refer to the same experiment of Fig. 3 and 4, but for the
oscillator $\mathcal{B}$. In this case the white noise floor increases
monotonically with the attenuation, but there is a significant discrepancy
between the experimental data and the “attenuator only” floor predicted by
(5). Additionally, dips are seen on Fig. 6-B at 2–20 kHz, more noticeable than
on Fig. 3-B.

Figure 8: PM noise of a “Golden Citrine” 100-MHz oscillator, measured with a
hacked FSWP. The polynomial fit gives $\mathsf{b}_{-4}=1.51{\times}10^{-7}$
$\mathrm{rad^{2}\,Hz^{3}}$ ($-68.2$ dB), $\mathsf{b}_{-3}=3.8{\times}10^{-9}$
$\mathrm{rad^{2}\,Hz^{2}}$ ($-84.2$ dB),
$\mathsf{b}_{-2}=3.39{\times}10^{-11}$ $\mathrm{rad^{2}\,Hz}$ ($-104.7$ dB),
$\mathsf{b}_{-1}=-6.31{\times}10^{-16}$ $\mathrm{rad^{2}}$
($|\mathsf{b}_{1}|=-152$ $\mathrm{dB\,rad^{2}}$), and
$\mathsf{b}_{0}=-1.26{\times}10^{-19}$ $\mathrm{rad^{2}/Hz}$
($|\mathsf{b}_{0}|=-189$ $\mathrm{dB\,rad^{2}/Hz}$). The region where
$\Re\\{S_{yx}\\}<0$ is not a valid PM-noise spectrum, but $-\Re\\{S_{yx}\\}$
provides useful information related to $\left|S_{yx}\right|$ because
$\Re^{2}\\{S_{yx}\\}\gg\Im^{2}\\{S_{yx}\\}$ almost everywhere. The ancient-
Latin expression HIC SVNT LEONES, usually translated as “here be dragons,”
refers to an unexplored land, or to a land where humans are not permitted.

Inspired by the theory (Sec. 3.1), we hacked a FSWP at the Rohde Schwarz R&D
facility in München, extracting $\Re\\{\left<S_{yx}(f)\right>\\}$ and
$\Im\\{\left<S_{yx}(f)\right>\\}$. This instrument is of the same type of that
we have in Besancon. In München we measured a third oscillator $\mathcal{C}$,
a 100-MHz Wenzel 501-25900B “Golden Citrine” OCXO, same brand and type of
$\mathcal{B}$. The DUT is connected via a 3-dB attenuator, and the FSWP had
internal 5-dB attenuation mechanically switched for better impedance matching.
Additionally, there is a 2.4 dB (typical) loss inside the FSWP, before the
power splitter. All losses accounted for, the signal level at the power
splitter input is 8.4 dBm, measured with the internal power meter. The result
is shown on Fig. 8. The phase noise is represented as
$\Re\\{\left<S_{yx}\right>\\}$. The negative, invalid outcomes are replaced
with $-\Re\\{\left<S_{yx}\right>\\}$ and shown in different color. The
quantity $|\Im\\{\left<S_{yx}\right>\\}|$ gives an indication about the
averaging limit of the instrument. Because
$|\Im\\{\left<S_{yx}\right>\\}|\ll|\Re\\{\left<S_{yx}\right>\\}|$ almost
everywhere in the spectrum, $|\Re\\{\left<S_{yx}\right>\\}|$ is a good
approximation of $|\left<S_{yx}\right>|$.

##  5 Interpretation

The dips found at 1–1.5 kHz in Fig. 3, and also at 2–20 kHz in Fig. 6, suggest
that $S_{\varphi}(f)$ changes sign at these points, being $S_{\varphi}(f)>0$
for $f<f_{\text{dip}}$, and $S_{\varphi}(f)<0$ beyond. The sign change occurs
because of $\varsigma S_{d}$ in (15), related to the fact that the displayed
$S_{\varphi}(f)$ is actually $|\left<S_{yx}(f)\right>|$, where $x$ and $y$ are
the phase of the DUT measured by the two channels inside the instrument. The
absolute value turns the sign-change into the sharp dip observed on the log
scale. This is experimentally confirmed in Fig. 8. By the way, the presence of
such dips was already predicted by a simulation in [5, Fig. 1(b) and Fig. 3].

From a theoretical standpoint, the combined effect of the attenuator (5) and
of the power splitter (18) results in

$\mathbb{E}\left\\{\mathsf{b}_{o}\right\\}=\frac{kT_{i}}{P_{i}}+\frac{k(1-A^{2})T_{a}}{A^{2}P_{i}}-\frac{kT_{s}}{A^{2}P_{i}}$
(19)

at the attenuator output. This contains two systematic effects: the attenuator
noise (positive), and the thermal energy of the power splitter (negative). At
high attenuation ($A^{2}\rightarrow 0$), the RF spectrum associated to the
noise sidebands tends to $kT_{a}$. In this condition, (19) predicts
$\mathrm{b}_{o}<0$ because the temperature $T_{s}$ inside the instrument is
obviously higher than the attenuator (and room) temperature $T_{a}$.

Let us start from the oscillator $\mathcal{A}$, the old Wenzel 501-04623E
(Fig.4). Using the absolute-value estimator, the expected $\mathsf{b}_{o}$ is

$\displaystyle\mathbb{E}\\{\widehat{\mathsf{b}_{o}}\\}$
$\displaystyle=\left|\frac{kT_{i}}{P_{i}}+\frac{k(1-A^{2})T_{a}}{A^{2}P_{i}}-\frac{kT_{s}}{A^{2}P_{i}}\right|~{}.$
(20)

Fitting the experimental points with (20) fails because there results a too
high $T_{s}$. Because the isolation between channels cannot be perfect, we
replace $T_{s}$ with $T_{s}-\varsigma T_{c}$, where $\varsigma T_{c}$
expresses the crosstalk given in terms of a temperature, and $\varsigma$ has
the same meaning as in (15). Accordingly, (20) rewrites as

$\displaystyle\mathbb{E}\\{\widehat{\mathsf{b}_{o}}\\}$
$\displaystyle=\left|\frac{kT_{i}}{P_{i}}+\frac{k(1-A^{2})T_{a}}{A^{2}P_{i}}+\frac{k(\varsigma
T_{c}-T_{s})}{A^{2}P_{i}}\right|~{}.$ (21)

Notice that there are two unknowns in (21), $T_{i}$ and $\varsigma
T_{c}-T_{s}$. The former is dominant at no attenuation ($A^{2}=1$), where the
observed PM noise is rather high. The latter is dominant at high attenuation
($A^{2}\rightarrow 0$). Because $\varsigma T_{c}-T_{s}$ appears as a single
quantity in (21), separating $\varsigma T_{c}$ from $T_{s}$ is somewhat
artificial, but it is useful in that it provides physical insight. We assume
$T_{a}=295$ K (23 $\mathrm{{}^{\circ}C}$) and $T_{s}=320$ K (47
$\mathrm{{}^{\circ}C}$) a convenient round number quite plausible for the
instrument inside. Fitting the data of Fig. 4 with (21) results in
$T_{i}=4528$ K and $T_{c}=122$ K. This is the curve labeled “full model.”
Using $\mathsf{b}_{i}=kT_{i}/P_{i}$, with $P_{i}=9.6$ mW ($+9.8$ dBm at
$A^{2}=1$), we get $\mathsf{b}_{i}=6.5{\times}10^{-10}$ $\mathrm{rad^{2}/Hz}$
($-171.9$ $\mathrm{dB\,rad^{2}/Hz}$). Comparing this value to the readout
($-172.4$ $\mathrm{dB\,rad^{2}/Hz}$ at $A^{2}=1$), the instrument introduces a
bias of $-0.5$ dB due to the combined effect of power splitter and crosstalk.

Removing the absolute value in (21) yields

$\displaystyle\mathbb{E}\\{\widehat{\mathsf{b}_{o}}\\}$
$\displaystyle=\frac{kT_{i}}{P_{i}}+\frac{k(1-A^{2})T_{a}}{A^{2}P_{i}}+\frac{k(\varsigma
T_{c}-T_{s})}{A^{2}P_{i}}~{},$ (22)

which results in $\mathsf{b}_{o}>0$ up to 15 dB attenuation, and in
$\mathsf{b}_{o}<0$ beyond.

Rewriting the polynomial model (1) for the absolute-value estimator we get

$\mathbb{E}\\{\widehat{S_{\varphi}(f)}\\}=\left|\frac{\mathsf{b}_{-3}}{f^{3}}+\frac{\mathsf{b}_{-2}}{f^{2}}+\frac{\mathsf{b}_{-1}}{f}+\frac{kT_{i}}{P_{i}}+\frac{k(1-A^{2})T_{a}}{A^{2}P_{i}}+\frac{k(\varsigma
T_{c}-T_{s})}{A^{2}P_{i}}\right|\,.$ (23)

Evaluating (23) with $\mathsf{b}_{-3}=3.5{\times}10^{-8}$
$\mathrm{rad^{2}Hz^{2}}$ ($-74.5$ $\mathrm{dB\,rad^{2}Hz^{2}}$),
$\mathsf{b}_{-2}=0$, and $\mathrm{b}_{-1}=4{\times}10^{-14}$
$\mathrm{dB\,rad^{2}}$ ($-134$ $\mathrm{dB\,rad^{2}}$), taken from Fig. 3-A,
we find the solid lines overlapped to the experimental spectra of Fig. 3-B.
The model matches the experiment, and predicts precisely the dips. These dips
occurr at $\geq 18$ dB attenuation, where $\mathsf{b}_{o}<0$.

Now we turn our attention to the oscillator $\mathcal{B}$, the Wenzel
501-25900B “Golden Citrine.”. Looking at Fig. 6-A and Fig. 7, we notice that
the white noise floor increases monotonically increasing the attenuation, and
the dips are present for all the values of the attenuation — albeit these dips
are not clear at 0 dB and 6 dB because of insufficient averaging. This
indicates that $\mathsf{b}_{o}<0$ in all cases. Evaluating (23) with the same
$T_{a}$, $T_{s}$ and $T_{c}$ as above, we find $T_{\text{eq}}=T_{i}=50$ K. The
model fits well the experimental data, as shown on Fig. 7. The temperature of
50 K is equivalent to a white noise floor of $-200$ $\mathrm{dB\,rad^{2}/Hz}$
at $+18.5$ dBm (70.5 mW) output power, with no attenuation. Finally, (23)
predicts precisely the dips seen on at Fig. (6)-B.

###  5.1 The Origin of the Crosstalk

Trying to understand the crosstalk, we look at the part of the FSWP where the
strongest and the weakest signals come close to one another, which is the
input mixer. Let us put numbers together with this idea. For linear
conversion, the LO signal should not be lower than $+20$ dBm. The phase noise
of a state-of the art synthesizer at 100 MHz carrier is of the order of $-160$
$\mathrm{dB\,rad^{2}/Hz}$. For reference, the R&S SMA100A synthesizer with the
low-phase-noise option SMA-B22 has a white floor of this order [28, data
sheet, p. 12]. At $+20$ dBm power, the white-noise sidebands are of $-140$
dBm/Hz, that is, $10^{-17}$ W/Hz. The crosstalk $kT_{c}$ we search for is of
$1.7{\times}10^{-21}$ W/Hz with $T_{c}=122$ K. This is 38 dB smaller than the
LO sidebands. A coupling of the order of $-38$ dB due to leakage is quite
plausible for a good mixer circuit. Besides, the absence of discontinuity in
the spectrum (Fig. 3 and 6) at 1 MHz indicates that the crosstalk does not
depend on the operating mode, which excludes some other parts of the
instrument. The presence of a small amount of anticorrelated flicker PM is
also possible, for the same reason. Anyway, this interpretation is just a
guess, not based on the internal design nor on specific measurements.

###  5.2 Inside the Oscillator

Figure 9: Simplified scheme of the low-noise quartz oscillators. The key point
is the interplay of filters and impedances. Otherwise, commercial oscillators
may differ from the schemes shown.

We address the question of the origin of $T_{\text{eq}}$, and why it can be
smaller than the room temperature. From our purposes, the oscillator consists
of a core (the auto-oscillator in strict sense), a buffer, and an output
filter (Fig. 9). The attenuation in the filter stopband is generally achieved
by reflecting the power back to the generator’s internal impedance.

The _conventional_ oscillators may have a lowpass or bandpass RLC filter at
the output to suppress the harmonic distortion and to solve other practical
problems (Fig. 9-A). Such filter cannot have a bandwidth smaller than a few
MHz at 100 MHz carrier because the quality factor $Q$ of these resonators is
of the order of 10–20 in practical conditions. As a consequence, the output
impedance is reasonably matched in the whole Fourier-frequency span, and the
white noise is chiefly the noise of the sustaining amplifier, where the
carrier is the weakest.

In the _thermally limited_ quartz oscillator, a quartz resonator is present
between the core and the buffer. Such filter can be the main resonator if the
carrier is extracted from the resonator’s ground pin [29], [30, Fig. 4-58 to
4-61], or a second quartz resonator (Fig. 9-B). Out of the resonator bandwidth
$\nu_{0}(1\pm 1/2Q)$, the quartz is a high impedance circuit, thus the noise
of the sustaining amplifier is not transmitted to the buffer. The noise
associated to the resonator’s motional resistance is also rejected, for the
same reason. The buffer (a common-base amplifier) has low input impedance and
low noise figure, thus the white noise is chiefly limited by the physical
temperature of the collector resistor $R_{C}$ at the output. Such oscillators
may have an additional RLC output filter of the same type discussed before.

In the _sub-thermally limited_ quartz oscillator, a quartz resonator or a
quartz filter is introduced in series to the output [31, Fig. 7], with no
further amplification. The output filter has a small cutoff frequency even
with the low $Q$ imposed by the heavy load condition. For example, taking
$Q=5000$ at 100 MHz, the cutoff frequency is $f_{c}=10$ kHz. For comparison, a
good resonator at this frequency has $Q>10^{5}$, unloaded. Out of the
bandwidth $\nu_{0}(1\pm 1/2Q)$, the output impedance is quite high
($|Z_{o}|\gg 50~{}\Omega$), which gives the appearance of a _cold_ source.
There no violation of the second principle because the filter is obviously in
thermal equilibrium with the environment. However, the _electrical access_ to
the thermal energy is open. In this condition, the input power splitter of the
noise analyzer is reasonably well matched only in the pass band, and nearly
open circuit in the stopband. In the stopband, the expected cross spectrum
relates to the thermal energy of the power splitter (and to the crosstalk, if
any), wich has negative sign in the correlation.

Simple attempts to measure the output impedance failed because the impedance
analyzers do not work in the presence of the strong carrier at the input. We
disassembled two 100 MHz oscillators, a Wenzel 501-04623E and a Wenzel
Citrine, the same type as the oscillator $\mathcal{A}$ and $\mathcal{B}$,
respectively. The oscillator $\mathcal{A}$ is of the conventional type, with a
RLC filter at the output. The white PM noise limited by the signal-to-noise
ratio in the sustaining amplifier. The oscillator $\mathcal{B}$ is of the sub-
thermally limited type, with a quartz resonator in series to the output.
Albeit we did not reverse-engineer the oscillator, the two values of
$T_{\text{eq}}$, 4528 K and 50 K, are consistent with the oscillator
architecture.

##  6 Discussion

We have seen that the (anti-)correlated noise inside the instrument can be
modeled as a temperature, which is $\varsigma T_{c}-T_{s}$. Let us look at
$T_{s}$ and $T_{c}$ separately. Because (18) is based on simple and well-
established physics, a software correction inside the instrument can
compensate for $T_{s}$ in a reliable way. An accuracy of a few kelvins is all
what is needed. Likewise, (5) can be used to correct for the effect of a
switchable dissipative attenuator, if present. By contrast, there is no
general way to compensate for $T_{c}$. We have no a priori reason to trust it
as a constant in the carrier-frequency range (4 decades), nor as reproducible
parameter across different specimens or architectures. The brute force
approach of putting the power splitter in a liquid-He cryostat [7] is not
effective because of the crosstalk. In our experiment 70% of the bias error is
due to the power splitter, 30% to the crosstalk. However, compensating for
$T_{s}$ alone is a general solution, and mitigates the problem.

The real-part estimator $\Re\left\\{\left<S_{yx}\right>_{m}\right\\}$ is
superior to the traditional estimator $\left|\left<S_{yx}\right>_{m}\right|$
in that (i) it converges faster because the background noise in
$\Im\left\\{\left<S_{yx}\right>_{m}\right\\}$ is not taken in, and (ii) it
reveals the negative, nonsensical outcomes.

Measuring the oscillator $\mathcal{A}$, the experimentalist may be satisfied
of the spectra taken with no attenuation ($A^{2}=1$) because

  * •

Two instruments from the major brands, with similar correlation algorithm but
radically different in the RF architecture and in the detection principle, are
in perfect agreement.

  * •

The systematic error in the white noise, revealed by our rather complex
experiment, is of a mere $-0.5$ dB, not alarming.

Conversely, the white noise floor measured on the oscillator $\mathcal{B}$ is
a complete nonsense because $\Re\left\\{\left<S_{yx}(f)\right>_{m}\right\\}<0$
inside the instrument.

Unlike most domains of metrology (mass, length, etc.), a PM noise spectrum
consists of hundreds or thousands of points on the $S_{\varphi}(f)$ plot. The
common ditto _too much information is no information_ rises the question of
the nature of the measurand. General experience indicates that the polynomial
law (1) describes well the PM noise spectrum of quartz ad dielectric
oscillators, thus a small number (4–5) of parameters $\mathsf{b}_{n}$ tell the
whole story. In optics, some additional terms appear, like bumps and blue
noise, which call for a small number of additional coefficients [32].
Regardless of the model we choose, a small number of missing points, like the
negative spurs we have seen, is not a real nuisance and can be ignored. The
experienced scientist does this after visual inspection. By contrast, an
irregular behavior over a wide frequency range has to be taken seriously.

Ultimately, some concepts found in the International Vocabulary of Metrology
(VIM) [33] and in the Guide to the Expression of Uncertainty in Measurement
(GUM) [34] (see also [35, 36]) should be introduced in phase noise
measurements. Going through the VIM, the following definitions are relevant to
our experiments: _Type A_ and _Type B_ evaluation of uncertainty (2.28 and
2.29), the _influence quantities_ (2.52), the _definitional uncertainty_
(2.27), and the _null measurement uncertainty_ (4.29). Because none of us is a
true expert of uncertainty in metrology, subtleties may escape from our
attention. However, this article shows that the assessment of uncertainty in
PM noise is still at a too rudimentary stage. The _following digression_ is
intended to _stimulate a discussion_ , with no intention of stating rules.

The type A uncertainty $u_{A}$ can be processed by a statistical analysis of
the time series, or in our case of a series of spectra. Conversely, the type B
uncertainty $u_{B}$ can be determined by other means, chiefly the analysis of
the system. The combined uncertainty (GUM Sec. 2.3.4 and Sec. 5) is
$u_{C}=\sqrt{u_{A}^{2}+u_{B}^{2}}$. In engineering, implicit reference is
often made to the _expanded uncertainty_ (VIM 2.35), with a coverage
probability of 95 %.

Most of our practical knowledge about uncertainty comes from Fig. 3 to 7. By
contrast, Fig. 8 provides the strong evidence of negative outcomes in a large
portion of spectrum. In turn, Fig. 8 supports our conclusion that some
portions of spectra are correctly interpreted as negative outcomes
($\mathsf{b}_{0}<0$ on Fig. 3 to 7), made positive by the absolute-value
estimator.

The single-channel background of the instrument is chiefly zero-average
Gaussian noise with white or colored spectral distribution, thus it falls in
the type A uncertainty. Averaging on $m$ spectra, the single-channel is
reduced by a factor of $1/\sqrt{n}$ for the absolute value estimator, and
$1/\sqrt{2m}$ for the real-part estimator. The practical impact on the
estimation of the white noise can be easily made negligible by averaging on a
large $m$, say $10^{6}$, ultimately limited by the time-frequency uncertainty
theorem.

Two contributions to the type B uncertainty are obviously identified, a
calibration factor and the correlated effects discussed. This suggests a
minimalist model like $u_{B}(S_{\varphi})=a_{1}S_{\varphi}+a_{0}$. The
calibration factor $a_{1}$ is dominant on the left-hand side of the spectrum,
where $S_{\varphi}(f)$ is higher. Based on the $1/f^{3}$ noise seen on Fig. 3
and 5 (same oscillator, measured with two instruments of different brand, RF
architecture, and principles), we can infer that the uncertainty is not
greater than 0.2–0.3 dB, i.e., $\lesssim 2\,\%$ in the phase-to-number
conversion. The same is expected on Fig. 8. Anyway such small value is not an
issue in the laboratory practice. By contrast, the noise $k(\varsigma
T_{c}-T_{s})/P$ show up in the white noise region, which is our main concern.
Accepting to go through a complex process, we can measure $k(\varsigma
T_{c}-T_{s})/P$ and corrected for it, with a residual uncertainty. Otherwise,
we can use the knowledge gathered to infer the uncertainty in a similar
experiment.

Figure 10: Uncertainty concepts, adapted from the VIM [33].

Here we see the importance of the null measurement uncertainty. The quantity
$k(\varsigma T_{c}-T_{s})/P$, or its uncertainty after correction, sets the
detection threshold, i.e., theminimum value of $S_{\varphi}(f)$ that can be
measured. This concept is illustrated in Fig. 10. The cases (A) and (B) comply
with rule that $S_{\varphi}(f)\pm u>0$, with (B) having smaller $u$ as a
result of a larger number of averaged spectra (same $u_{B}$ but smaller
$u_{A}$). By contrast, in (C) the uncertainty bar hits the negative values. In
this case, the correct way to express the result result is $S_{\varphi}(f)=0$,
with a detection threshold $u$. The case (D) should be discarded at sight,
giving the same result $S_{\varphi}(f)=0$ as in (C).

Let us show the above concepts on the spectrum of Fig. 8. In this case we did
only one experiment, thus we have no measure of $\varsigma T_{c}-T_{s}$, and
no correction for $T_{s}$ is implemented inside the FSWP. Thus, we take
$|\varsigma T_{c}-T_{s}|$ rounded to 500 K as a convenient estimate of the
instrument limit. With $P=6.9$ mW ($+8.4$ dBm), $\mathsf{b}_{0}=kT/P$ gives
$u(S_{\varphi})=10^{-18}$ $\mathrm{rad^{2}/Hz}$. Accordingly, the detection
threshold is of $-180$ $\mathrm{dB\,rad^{2}/Hz}$. The forbidden region from 15
kHz to 3 MHz, where $\Re\\{S_{\varphi}(f)\\}<0$, corresponds to the case (D)
of Fig. 10. Beyond 3 MHz the spectrum is not trusted, being affected by other
limitations. Anyway, the latter region has little or no practical importance.

Finally, we have seen that the output impedance $Z_{o}(f)$ produces erratic
results if it changes significantly in the analysis bandwidth. This opens the
question of whether $Z_{o}(f)$ goes in the definitional uncertainty (it is
inside the DUT), it goes in the B-type uncertainty, or if it is an influence
quantity. The role of impedance mismatch is well known in microwave noise
measurements [37, 38], but these concepts have not been transposed to PM
noise.

##  7 Conclusions

Our method consists of introducing various values of dissipative attenuation
between the oscillator under test and the phase-noise analyzer. This method is
new. It provides quantitative information about the unwanted correlated
effects inside the analyzer, and helps to assess the _null measurement
uncertainty_ , i.e., the minimum amount of phase noise that can be detected.
In some circumstances, inserting an attenuator results in lower white PM noise
floor. When this happens, gross errors are around the corner. The idea that
the (anti-)correlated noise inside the instrument can be described in terms of
the thermal energy $k(\varsigma T_{c}-T_{s})$ is also new. This parameter
accounts for the temperature of the power splitter at the instrument input,
and the crosstalk between the two channels.

The experiments described provide the evidence that pushing the noise
rejection too far by averaging on a large number of data may result in
misleading or grossly wrong results. The reason is in residual correlated
effects, not under control. In general terms, under-estimating the DUT noise
is obviously worse than over-estimating it.

Impedance matching in the whole analysis bandwidth is a critical issue. Sub-
thermally limited oscillators make use of a narrowband reactive filter at the
output, which exploits impedance mismatch in the stopband to deliver the
lowest noise floor. However, such filter results in anticorrelated noise due
to the thermal energy in the power splitter at the instrument input. From a
different standpoint, the benefit of a sub-thermally limited oscillator is
unclear to us if the oscillator is intended to be a part of a system at room
temperature.

## Disclaimer

Our strong statements require an equally strong disclaimer about the
commercial products we refer to. We experimented on them because they were on
hand at the right time, as opposite to gathering parts with this research in
mind. By no means we criticize these products, nor we endorse them. The
problems and the inconsistencies we describe relate to unintended, strange, or
weird use of these products. Driven by the genuine scientific curiosity, we
share our knowledge with the ultimate intent to contribute to better
understanding the physics and the technology of phase noise metrology. We hope
that no misunderstanding will arise, and we apologize if this will happen.

## Acknowledgments

This work is funded by the ANR Programme d’Investissement d’Avenir under the
Oscillator IMP project (contract ANR-11-EQPX-0033-OSC-IMP) and the FIRST-TF
network (contract ANR-10-LABX-48-01), and by the Région Bourgogne Franche
Comté. We thank AR Electronique (the similarity between this brand and the
initials of one of the authors is a random outcome), Besancon, for hosting us
in their R&D lab for some additional measurements, and Cristian Bolovan and
Radu Ohreac, Rohde & Schwarz Romania, for generously lending a FSWP 8.

## References

  * [1] F. L. Walls, S. R. Stain, J. E. Gray, and D. J. Glaze, “Design considerations in state-of-the-art signal processing and phase noise measurement systems,” in _Proc. Freq. Control Symp._ , Atlantic City, NJ, USA, Jun. 1976, pp. 269–274.
  * [2] W. F. Walls, “Cross-correlation phase noise measurements,” in _Proc. Freq. Control Symp._ , Hershey, PA, May 1992, pp. 257–261.
  * [3] E. Rubiola and V. Giordano, “Correlation-based phase noise measurements,” _Rev. Sci. Instrum._ , vol. 71, no. 8, pp. 3085–3091, Aug. 2000.
  * [4] E. N. Ivanov and F. L. Walls, “Interpreting anomalously low voltage noise floors in two-channel measurement systems,” _IEEE Trans. Ultras. Ferroelec. Freq. Contr._ , vol. 49, no. 1, pp. 11–19, Jan. 2002.
  * [5] C. W. Nelson, A. Hati, and D. A. Howe, “A collapse of the cross-spectral function in phase noise metrology,” _Rev. Sci. Instrum._ , vol. 85, no. 3, pp. 024 705, 1–7, Mar. 2014.
  * [6] A. Hati, C. W. Nelson, and D. A. Howe, “Cross-spectrum measurement of thermal-noise limited oscillators,” _Rev. Sci. Instrum._ , vol. 87, p. 034708, Mar. 2016.
  * [7] A. Hati, C. W. Nelson, D. P. Pappas, and D. A. Howe, “Phase noise measurements with a cryogenic power-splitter to minimize the cross-spectral collapse effect,” _Rev. Sci. Instrum._ , vol. 88, p. 114707, Nov. 2017.
  * [8] J. Breitbarth, “Source impedance influence on cross-correlation phase noise measurements,” in _Proc. Intl. Freq. Control Symp. and Europ. Freq. Time Forum Joint Meeting_ , Prague, Czech Republic, Jul. 2013, pp. 434–437.
  * [9] “Cross Spectrum $\mathscr{L}(f)$ Workshop,” Besançon, France, Jul. 2017, side event at the 2017 Europ. Freq. Time Forum and Freq. Control Symp. Joint Meeting, July 10–13.
  * [10] “Cross Spectrum $\mathscr{L}(f)$ Workshop,” Denver, CO, USA, Apr. 2015, side event at the 2015 Europ. Freq. Time Forum and Freq. Control Symp. Joint Meeting, April 13–16.
  * [11] “European Workshop on Cross-Spectrum Phase Noise Measurements,” Laboratoire National de métrologie et d’Essays (LNE), Paris, France, Dec. 2014\.
  * [12] E. S. Ferre-Pikal, _IEEE Standard Definitions of Physical Quantities for Fundamental Frequency and Time Metrology–Random Instabilities (IEEE Standard 1139-2008)_ , IEEE, New York, Feb. 2009.
  * [13] M. E. Tiuri, _Radio-Telescope Receivers,_Chap. 7 of J. D. Kraus,_ Radio Astronomy_. New York: McGraw Hill, 1966.
  * [14] K. Rohlfs and T. L. Wilson, _Tools of Radio Astronomy_ , 3rd ed. Berlin, Heidelberg, New York: Springer Verlag, 2000\.
  * [15] E. Rubiola and F. Vernotte, “The cross-spectrum experimental method,” arXiv:1004.5539 [physics.ins-det], Apr. 2010.
  * [16] M. Sampietro, L. Fasoli, and G. Ferrari, “Spectrum analyzer with noise reduction by cross-correlation technique on two channels,” _Rev. Sci. Instrum._ , vol. 70, no. 5, pp. 2520–2525, May 1999.
  * [17] G. Feldhaus and A. Roth, “A 1 MHz to 50 GHz direct down-conversion phase noise analyzer with cross-correlation,” in _Proc. Europ. Freq. Time Forum_ , York, UK, Apr. 4–7, 2016, Conference Proceedings.
  * [18] Y. Gruson, V. Giordano, U. L. Rohde, A. K. Poddar, and E. Rubiola, “Cross-spectrum pm noise measurement, thermal energy, and metamaterial filters,” _IEEE Trans. Ultras. Ferroelec. Freq. Contr._ , vol. 64, no. 3, pp. 634–642, Mar. 2017.
  * [19] E. Rubiola and R. Boudot, “The effect of AM noise on correlation phase noise measurements,” _IEEE Trans. Ultras. Ferroelec. Freq. Contr._ , vol. 54, no. 5, pp. 926–932, May 2007, also arXiv:physics/0609147.
  * [20] J. Fischer, “Progress towards a new definition of the kelvin,” _Metrologia_ , vol. 52, no. 5, pp. S364–S375, Aug. 2015.
  * [21] D. R. White, R. Galleano, A. Actis, H. Brixy, M. De Groot, J. Dubbeldam, A. L. Reesink, F. Edler, H. Sakurai, R. L. Shepard, and G. J. C., “The status of Johnson noise thermometry,” _Metrologia_ , vol. 33, no. 4, pp. 325–335, Aug. 1996.
  * [22] C. M. Allred, “A precision noise spectral density comparator,” _J. Res. NBS_ , vol. 66C, pp. 323–330, Oct-Dec 1962.
  * [23] G. J. Halford, “Noise comparators for s and x bands,” _IEEE Trans. Instrum. Meas._ , vol. 15, no. 4, pp. 310–317, Dec. 1966.
  * [24] Rohde & Schwarz, “R&S FSWP Phase Noise Analyzer and VCO Tester,” Brochures and Data Sheets available on https://www.rohde-schwarz.com/us/brochure-datasheet/fswp/, accessed December 21, 2019.
  * [25] U. L. Rohde, J. C. Whitaker, and H. Zahnd, _Communications Receivers_ , 4th ed. New York: Mc Graw Hill Education, 2017.
  * [26] Keysight Technologies, Inc., “E5500 series, phase noise measurement solutions,” available on https://www.keysight.com/us/en/assets/7018-01214/data-sheets/5989-0851.pdf, accessed December 21, 2019.
  * [27] Wenzel Associates, Inc., “Low Noise Crystal Oscillators > VHF Citrine,” brochure available on http://www.wenzel.com/wp-content/uploads/VHF-Citrine.pdf, accessed December 21, 2019.
  * [28] Rohde & Schwarz, “R&S SMA100A Signal Generator,” Brochures and Data Sheets available on https://www.rohde-schwarz.com/us/brochure-datasheet/sma100a/, accessed December 21, 2019.
  * [29] U. L. Rohde, “Crystal oscillator provides low noise,” _Electronic Design_ , vol. 21, p. 11 and 14, Oct. 11, 1975.
  * [30] ——, _Microwave and Wireless Synthesizers_. New York: John Wiley & Sons, 1997.
  * [31] ——, “Mathematical analysis and design of ultra stable low noise 100 MHz crystal oscillator with differential limiter and its possibilities in frequency standards,” in _Proc. Freq. Control Symp._ , May 31 – Jun. 2 1978, pp. 409–425.
  * [32] C. E. Calosso, C. Clivati, and S. Micalizio, “Avoiding aliasing in allan variance: An application to fiber link data analysis,” _IEEE Trans. Ultras. Ferroelec. Freq. Contr._ , vol. 63, no. 4, pp. 646–655, Apr. 2016.
  * [33] Joint Committee for Guides in Metrology (JCGM), _International Vocabulary of Metrology — Basic and General Concepts and Associated Terms (VIM), document JCGM 200:2012_ , 3rd ed., 2012, available on https://www.bipm.org/en/publications/guides/.
  * [34] ——, _Evaluation of measurement data –– Guide to the expression of uncertainty in measurement (GUM), document JCGM 100:2008_ , 1st ed., Sep. 2008, available on https://www.bipm.org/en/publications/guides/.
  * [35] B. D. Hall, “Evaluating the measurement uncertainty of complex quantities: a selective review,” _Metrologia_ , vol. 53, pp. S25–S31, 2016.
  * [36] A. Possolo and H. K. Iyer, “Concepts and tools for the evaluation of measurement uncertainty,” _Rev. Sci. Instrum._ , vol. 88, p. 011301, 2017\.
  * [37] T. Y. Otoshi, “The effect of mismatched components on microwave noise-temperature calibrations,” _IEEE Trans. Microw. Theory Tech._ , vol. 16, no. 9, pp. 675–686, Sep. 1968.
  * [38] S. M. Wedge and D. B. Rutledge, “Wave techniques for noise modeling and measurement,” _IEEE Trans. Microw. Theory Tech._ , vol. 40, no. 11, pp. 2004–2012, Nov. 1992.


Domain randomization (DR) cannot provide optimal policies for adapting the learning agent to the dynamics of the environment, although it can generalize sub-optimal policies to work in a transferred domain. In this paper, we present Universal Policy with Embedding System Identification (UPESI) as an implicit system identification (SI) approach with universal policies (UPs), as a learning-based control method to execute optimal actions adaptively in environments with various dynamic properties. Previous approaches of SI for adaptive policies either conduct explicit SI, which is testified to be an ill-posed problem, or suffer from low efficiency without leveraging the simulation oracle. We propose to conduct SI in the embedding space of system dynamics by leveraging a learned forward dynamics model, and use Bayesian optimization for the SI process given transition data in a new environment. The identified embeddings are applied as additional input to the UP to enable its dynamics adaptability. Experiments demonstrate the advantageous performances of our proposed UP with embedding SI over standard DR and conventional SI approaches on both low-dimensional and high-dimensional simulation tasks.

§ INTRODUCTION

Robot learning [1, 2, 3] has been testified to successfully work not only in simulation but also on real-world robotic control. In this domain, reinforcement learning (RL) [4, 5, 6] methods are typically applied for robotic control via sim-to-real transfer [2, 7, 8]. As one of the most widely applied methods in sim-to-real research literature, domain randomization (DR) can help to learn policies that are more general and robust to be applied in various environments which do not permit direct access to their underlying dynamics. System identification (SI) is another main category of methods that can help to bridge the domain gap. SI methods usually configure the physical parameters from historical transition data, either in an explicit or implicit manner. Both approaches have been proven to be feasible for some control tasks [3, 1], when encountering the sim-to-real transfer problems, or more generally, domain transfer problems. However, a successful execution of the task does not necessarily indicate that an optimal control strategy is achieved, which leaves the space for further improvement based on solving the existing defects of above approaches.

Problems of existing methods are investigated previously. Conventional DR methods [9] are able to provide a policy with better generalization in various environments, but not the optimal policy for the settings in a specific environment. The optimization objective is additionally taken with respect to an expectation of the distribution of system dynamic properties, rather than directly on the true dynamic settings of the testing environment. Therefore, being aware of the true dynamics of a new environment is necessary for accomplishing optimal control on a task, which is lacking for DR in present works.

SI can help with bridging this gap via an active inferring of system dynamics with the necessary data collected.  In [1], an online system identification (OSI) module is used for explicitly configuring the dynamics parameters of environments, with a subsequent universal policy (UP) conditioned on the identified parameters to generate adaptive actions in new environments. However, both the UP and OSI modules can be infeasible in practice for general robot learning tasks according to the discussions in [8]. OSI may not be feasible since the trajectories may not provide enough information for accurately predicting the dynamics parameters in a point-estimation manner, therefore distributions can better represent the uncertainty in the dynamics identification. The problem can be ill-posed to directly predict dynamics parameters, according to [10]. Multiple potential dynamics models can lead to the same trajectories as a non-injective function, due to the entanglement of dynamics parameters. For UP, it can be hard to train when the dimension of dynamics parameters is high.

Present methods for system dynamics identification usually require aligned trajectories for constructing the contrastive loss [10, 11], which we find not to be a necessity. In previous methods, the trajectories alignment can usually be achieved with setting same initial states and executing the same policy, or rollout the simulation environment on a give state and certain action for a single timestep at each time. However, not all simulators can support a direct setup to be any state in the environments. Setting up the initial states only will lead to compounding errors on the trajectories due to the dynamics difference or different random seeds, or requires undesirable engineering efforts. Our proposed method solves the trajectories alignment procedure by leveraging on a well-trained forward dynamics prediction model.

This paper makes the following contributions:
1) We demonstrate the limited performances of typical DR methods when the dynamics of testing environment is uncertain or various, compared against our proposed method;
2) We propose a method to handle dynamics uncertainty and unmodeled effects [3] separately, with SI module for implicitly identifying the dynamics parameters, as well as DR for unmodeled effects involving observation noise, observation delay and action noise;
3) For system dynamics identification, our method does not depend on aligned trajectories collected in the source and target domains, but with randomly sampled trajectories carried out by the same control policy;

4) Our method learns a regularized dynamics embedding rather than applying the oracle dynamics parameters as the representation of system dynamics, to handle the entangled effects or redundant information within dynamics parameters. A universal policy conditioned on the learned dynamics embedding is trained with DR for both dynamics and unmodeled effects.

§ RELATED WORK

Existing methods for bridging the domain transfer gap (e.g., sim-to-real) can be broadly categorized in to the following classes: 1). Domain Randomization [9] randomizes either visual features [12] or dynamics parameters [2] in a source domain to train policies with better adaptability to the target domain. The branch of works also include the recent process of structured DR [13] and active DR [14]. 2). System Identification usually requires an extensive data collection and calibration process to mitigate the gap between the source and target environments. Explicit SI can be incorporated with a universal policy [1] for achieving an adaptive control of the real robot for various environment settings. Implicit SI usually uses recurrent units like long short-term memory (LSTM) networks with sequential inputs to preserves information about environment dynamics [2, 3]. 3). Domain Adaptation (DA) [15, 16] applies transfer learning techniques to match the distribution of source domain data with the target domain data, often applied in visuomotor control with images.  4). Strategy Optimization (SO) [17, 7, 18] requires evaluating a family of policies (called strategy) in the target domain and selecting the one with the best performance, usually with sampling-based methods like Bayesian optimization (BO) [19] and CMA-ES [20].

Our works focus on the SI approach with universal policies trained for domain transfer. DR in general cannot provide optimal policy due to the induced noise for increasing model generality. The optimized objective for the control policy is additionally taken with respect to an expectation of randomized parameters in environments. Thus only sub-optimal policies can be performed in the target domain even for the case that target sample data falls in the domain of randomized distributions of the source. Traditional explicit SI methods [1] directly configure the system parameter, which is testified to be an ill-posed problem [10, 8] due to that the entangled effects of multiple parameters will lead to non-unique identification results. This problem no longer exists for implicit SI [2, 3]. However, the true system parameters in the source domain are not leveraged for implicit SI, which hinders the learning efficiency of this approach. To this end, we found the embedding of system dynamics parameters can be achieved with encoders to improve the above approaches.

§ PRELIMINARIES

§.§ Notations

A typical formulation of RL problems follows a standard Markov decision process (MDP), which can be represented as $(\mathcal{S}, \mathcal{A}, R, \mathcal{T}, \gamma)$, where $\mathcal{S}$ and $\mathcal{A}$ are feasible sets of state and action, and $R$ is reward function $R(s,a)$: $\mathcal{S}\times \mathcal{A}\rightarrow \mathbb{R}$. $\mathcal{T}$ defines the transition probability from current state $s$ and action $a$ to a next state $s^\prime$ based on a fixed dynamics setting: $\mathcal{T}(s^\prime|s,a)$, and $\gamma\in(0,1)$ is a discount factor. However, it is not sufficient for a learned policy to be applied in a transferred domain. We consider a partially observable MDP (POMDP) with randomized dynamics as $(\mathcal{S}, \mathcal{A}, \Theta, \mathcal{T}, \mathcal{O}, \mathcal{P}_{\mathcal{O}}, R, \gamma)$, with additional $\Theta$ as the $\mathbb{R}^d$ space of dynamics parameters ($d$ is the number of dynamics parameters), $\mathcal{O}$ as the observation space and $\mathcal{P}_{\mathcal{O}}$ as the emission probability distribution: $\mathcal{P}_{\mathcal{O}}(o|s)$. Moreover, the transition probability $\mathcal{T}$ will then be different from the standard MDP, and it further depends on the dynamics parameters $\theta$ as $\mathcal{T}(s^\prime|s,a,\theta), \theta\in\Theta$. Note that here we define $\mathcal{P}_{\mathcal{O}}$ to be independent on dynamics parameters $\theta$, since we disentangle the observed transition process into the transition $\mathcal{T}$ in underlying state space $\mathcal{S}$ and the omission $\mathcal{P}_{\mathcal{O}}$ of observations from states, so only the transition of underlying states $\mathcal{T}$ depends on the dynamics parameters $\theta$.

§ METHODOLOGY

Overview of our method, including four steps: (a) training the dynamics prediction model with inputs of observation $o$, action $a$, and embedding $\alpha$ of system dynamics parameter $\theta$; (b) training the universal policy with a learned encoder for providing the embedding of system dynamics; (c) optimization in the embedding space with data from a target domain using Bayesian optimization and the learned forward prediction model; (d) inference in the target domain environment with identified dynamics embedding and the universal policy. The blue shaded arrows indicate the gradient flows during model training. The red arrows indicate a trained/fitted models or parameters are applied (with fixed weights/values) somewhere else.

In this paper, we propose the method for domain transfer with not only DR but also SI in an embedding space. The overview of our method is shown in Fig. <ref>. Three key components are involved in the framework:

Dynamics Encoder: $\mathbf{E}(\theta) \mapsto \alpha$

Forward Dynamics Predictor: $\mathbf{F}(o,a,\alpha)\mapsto o^\prime$

Universal Policy: $\pi(o, \alpha) \mapsto a$

The Dynamics Encoder $\mathbf{E}$ is an embedding network to generate low-dimensional embedding $\alpha$ from the system parameters $\theta$. In practical robotic learning tasks, dozens of system parameters can be involved in simulation [8, 3], which would severely increase the difficulty of universal policy learning. Moreover, explicit identification of each parameter value is neither necessary nor impossible for general cases [10]. For example, an increment of friction coefficients on a specific joint may be counteracted by a decrease of mass of some link bodies when we simply look at the trajectories executed by the robot with the same torques applied. The entangled effects of system parameters are considerable in SI process. The Forward Dynamics Predictor $\mathbf{F}$ is learned to mirror both the transition function $\mathcal{T}(s^\prime|s,a,\theta)$ of the environment with randomized dynamics $\theta\in\Theta$ and the observation emission function $\mathcal{P}_{\mathcal{O}}(o|s)$ for transferring states $s,s^\prime$ to observations. By leveraging the conditional dependency on embedding $\alpha$ in function $\mathbf{F}$, we propose to directly optimize $\alpha$ in the embedding space using Bayesian optimization when given the data from the target domain, for achieving the embedding SI. The Universal Policy $\pi(o, \alpha)$ basically follows the standard settings in [1], where a concatenation of observations and identified parameter embedding is taken as policy inputs. By conditioning the policy on the system properties, it is capable of adaptively selecting optimal actions in different system settings and randomized parameters $\theta$. We will show details in later sections.

§.§ Theoretical Insights

In the following, we provide the theoretical justification of our method. Since we assume the observation emission process has no impact on the dynamics of underlying state space, the analysis below will be based on a fully observable MDP with randomized dynamics but no observations. Apart from this, another difference of our implementation in practice from the theoretical analysis is the usage of dynamics embedding, i.e., using $\alpha$ instead of $\theta$. These two differences will not cause the loss of generality for our analysis.

Following a similar setup in [21], we extend the definitions of different occupancy measures for MDP with randomized dynamics as in Table <ref>.

tableDifferent occupancy measure for MDP with randomized dynamics

Occupancy Measure     State-Action      State Transitions     Joint w/o Dynamics     Joint w/ Dynamics

Denotation     $\rho_\pi(s,a, \theta)$      $\rho_\pi(s,s^\prime)$     $\rho_\pi(s,a,s^\prime)$     $\rho_\pi(s,a,\theta,s^\prime)$
Support     $\mathcal{S}\times\mathcal{A}\times\Theta$     $\mathcal{S}\times\mathcal{S}$     $\mathcal{S}\times\mathcal{A}\times\mathcal{S}$     $\mathcal{S}\times\mathcal{A}\times\Theta\times\mathcal{S}$
Definition     $\rho_\pi(s,\theta)\pi(a|s,\theta)$     $\int_{\mathcal{A}\times\Theta}\rho_\pi(s,a,\theta)\mathcal{T}(s^\prime|s,a,\theta)\text{d} a \text{d}\theta$
$\int_{\Theta}\rho_\pi(s,a,\theta)\mathcal{T}(s^\prime|s,a,\theta)\text{d}\theta$     $\rho_\pi(s,a,\theta)\mathcal{T}(s^\prime|s,a,\theta)$

1[Estimated Dynamics]
The estimated dynamics distribution is defined as the distribution of dynamic parameters estimated from given the transition tuples $\{(s,a,s')\}$, and it can be calculated as:

\begin{equation}
p(\theta|s,a,s') = \frac{\rho(s,a,\theta,s')}{\rho(s,a,s')} = \frac{\mathcal{T}(s'|s,a,\theta)\rho(s,a,\theta)}{\int_{\Theta}\mathcal{T}(s'|s,a,\bar{\theta})\rho(s,a,\bar{\theta})\text{d}\bar{\theta}}
\end{equation}

In practice, the estimated dynamics can be represented as a parameterized model with learnable parameters $\psi$ and optimized with dataset $\mathcal{D}^s$ from the source domain. During training or evaluation, the predicted distribution and true distribution are denoted as $p^s_\psi(\theta|s,a,s')$ and $p^t(\theta|s,a,s')$ respectively. The optimization objective for SI can be formulated as minimizing the discrepancy of estimated and true dynamics distributions, e.g., $\min_\psi \mathbb{D}_\text{KL}[p^s_\psi(\theta|s,a,s')||p^t(\theta|s,a,s')]$. We will prove that the method of leveraging the dynamics prediction model to achieve SI (step (c). in the scheme as in Fig. <ref>) is a valid approach for this objective.

Given the same dataset distribution $\rho^s(s,a,s')=\rho^t(s,a,s')$, the difference between the estimated dynamics distribution $p^s$ and the true distribution $p^t$ can be characterized as:

\begin{align}
&\mathbb{D}_\text{KL}[p^s(\theta|s,a,s')||p^t(\theta|s,a,s')] \\
% &=\mathbb{D}_\text{KL}[\rho^s(s,a,\theta,s')||\rho^t(s,a,\theta,s')] - \mathbb{D}_\text{KL}[\rho^s(s,a,s')||\rho^t(s,a,s')] \\
\end{align}

3[Forward Dynamics Prediction]
The forward dynamics prediction model is defined as the distribution of next state $s'$ given the current tuple $(s,a,\theta)$ in a MDP with randomized dynamics, as an approximation of the true underlying dynamics in the source/target domain:

\begin{equation}
f(s'|s,a,\theta)\approx \mathcal{T}(s'|s,a,\theta)=\frac{\rho(s,a,\theta,s')}{\rho(s,a,\theta)}
\end{equation}

In practice, this model approximation process can be achieved with optimization on collected dataset in the source domain. We denote the parameters of forward dynamics prediction as $\Phi$, which are learned in the approximation process with objective $\min_\Phi \mathbb{D}_\text{KL}[f^s_\Phi(s'|s,a, \theta)||f^t(s'|s,a, \theta)]$ (step (a).), where $\{(s,a,\theta)\}$ are given and $f^t(s'|s,a, \theta)$ are the true distributions of the forward prediction. After training, $f(s'|s,a,\theta)$ is evaluated on the testing dataste (e.g., the target domain dataset), with its predicted distribution denoted as $f^s_\Phi(s'|s,a,\theta)\approx\frac{\rho^s(s,a,\theta,s')}{\rho^s(s,a,\theta)}$ and true distribution as $f^t(s'|s,a,\theta)\approx\frac{\rho^t(s,a,\theta,s')}{\rho^t(s,a,\theta)}$. In our proposed method leveraging the forward dynamics prediction, the SI process is characterized as $\min_\theta \mathbb{D}_\text{KL}[f^s_\Phi(s'|s,a, \theta)||f^t(s'|s,a, \theta)]$ (step (c).), where $\{(s,a)\}$ and $\Phi$ are given. In conventional SI, without $f(s'|s,a,\theta)$, this optimization objective actually becomes $\mathbb{D}_\text{KL}[\mathcal{T}^s(s'|s,a, \theta)||\mathcal{T}^t(s'|s,a, \theta)]$, where $\mathcal{T}^s$ and $\mathcal{T}^t$ are the transition functions of the source and target domains respectively, and it requires the simulator to set the same $(s,a)$ as the real world does all the time for a sim-to-real transfer example.

Now we have another lemma of discrepancy between $f^s(s'|s,a,\theta)$ and $f^t(s'|s,a,\theta)$ as follows.
4 The distance of the distribution from the forward dynamics prediction and the true distribution can be formulated with the KL-divergence, it thus follows:

\begin{align}
% &\mathbb{D}_\text{KL}[f^s(s'|s,a, \theta)||f^t(s'|s,a, \theta)]\nonumber \\
% &\approx \mathbb{D}_\text{KL}[\rho^s(s,a,\theta,s')||\rho^t(s,a,\theta,s')]

\mathbb{D}_\text{KL}[f^s||f^t]\approx \mathbb{D}_\text{KL}[\rho^s||\rho^t] -\mathbb{D}_\text{KL}[\rho^s(s,a,\theta)||\rho^t(s,a,\theta)]
\end{align}

where $f^{s,t}$ are shorten for $f^{s,t}(s'|s,a, \theta)$, and $\rho^{s,t}$ are shorten for $\rho^{s,t}(s,a,\theta,s')$.

The optimization of forward dynamics prediction is increasing the lower bound of the optimization objective for improving the estimated dynamics, i.e.,

\begin{equation}
\mathbb{D}_\text{KL}[p^s||p^t]\geq  \mathbb{D}_\text{KL}[f^s||f^t]
\end{equation}

where $p^{s,t}$ are $p^{s,t}(\theta|s,a,s')$ and $f^{s,t}$ are $f^{s,t}(s'|s,a,\theta)$.

The above theorem tells us that a direct optimization of estimated dynamics for parameters $\phi$ can be approximately achieved with optimization of parameters $\Phi$ for forward dynamics prediction then leverage it to optimize $\theta$.
Proofs of above theorems and lemmas are all provided in App. <ref>.

In further analysis, we note that step (a) and (c) in our method are actually a SI process, not in original parameter space but an embedding space. As proposed in [16], the domain-invariant representation learning can be achieved by both minimizing the marginal discrepancy (called marginal distributions alignment) and minimizing the conditional discrepancy (called conditional distributions alignment) for an object classification task. Here we try to achieve a domain-invariant forward prediction model by optimizing the embedding of system dynamics. We will show how our method satisfies the process of minimizing conditional discrepancy. The reason that marginal distribution alignment is not achieved in current settings is due to the lack of true system dynamics parameters in the target domain, which are required in marginal distribution alignment.

For conditional distributions alignment, we first modify the definition in [16] with additional conditional variable. Then we narrow down the definition for general domain-transfer models to the specific forward dynamics prediction model in our method.
6[Conditional Distributions Alignment (modified from [16])]
Given two domains $D^s=\{x^s_i, y^s_i, z^s_i\}^{N^s}_{i=1}$ and $D^T=\{x^t_i, y^t_i, z^t_i\}^{N^t}_{i=1}$ drawn from random variables $(X^s\times Y^s\times Z^s)$ and $(X^t\times Y^t\times Z^t)$ with different output conditional probability distributions $Pr(Y^s|X^s, Z^s)\neq Pr(Y^t|X^t, Z^t)$, conditional alignment corresponds to finding the transformation $X\times Z \overset{g(X)}{\longrightarrow} E\times Z \overset{f(E,Z)}{\longrightarrow} Y$ such that the discrepancy between the transformed conditional distributions is minimized, i.e., $Pr(Y^s|g(X^s), Z^s)=Pr(Y^t|g(X^t), Z^t)$. Note that the number of additional variables $Z$ is not limited to be one.

In our case with forward dynamics prediction model $\mathbf{F}(o,a,\alpha)$ and dynamics encoder $\mathbf{E}(\theta)$, we have $Pr(o^{t \prime}|o^t, a^t, \theta^t)\neq Pr(o^{s \prime}|o^s, a^s, \theta^s)$, and conditional alignment for the forward dynamics prediction model optimizes the dynamics embedding $\alpha$ such that, based on the optimized dynamics embedding, the resulting prediction distributions are minimized, i.e., $Pr(o^{s \prime}|o^s,a^s;\mathbf{E}(\theta^s))=Pr(o^{t \prime}|o^t,a^t;\alpha^t)$. In the source domain the embedding $\alpha^s=\mathbf{E}(\theta^s)$ is optimized via back-propagating through the dynamics encoder $\mathbf{E}$, while in the target domain $\alpha^t$ is directly optimized with methods like BO due to the lack of true dynamics parameters. From above analysis, we show that our method (specifically step (a) and (c)) accomplishes an extended conditional distribution alignment for identifying the dynamics embedding.

§.§ Universal Policy with Embedding System Identification

In previous methods, SI is usually achieved with aligned trajectories [10, 11], which is found to be unnecessary in our method. We will detail the formulation of two approaches as follows. Suppose the datasets in the source and target domains are as $\mathcal{D}^s=\{(o^s, a^s, o^s{}^\prime)\}$ and $\mathcal{D}^t=\{(o^t, a^t, o^t{}^\prime)\}$, $\mathcal{D}^s$ is collected with dynamics randomization (without observation noise and action noise). Apart from deriving actions from the same policy as $a^s=\pi(o^s)$ and $a^t=\pi(o^t)$, most present methods assume the underlying states to be the same in the source and target domains, as $o^t\sim\mathcal{P_O}(\cdot|s)$ and $o^s\sim\mathcal{P_O}(\cdot|s)$ for $\forall o^t\in \mathcal{D}^t, o^s\in\mathcal{D}^s$, so as to define a loss function as $L(\{o^s{}^\prime, o^t{}^\prime\})$  for optimizing the identified system dynamics parameters. The trajectory in underlying state space $\tau=(s_{T=0}, s_{T=1}, s_{T=2}, ..., s_{T=n})$ is required to be aligned in the source domain and the target domain to form a valid loss function, especially for the initial state $s^{T=0}$. In our experiments, we find this assumption in dataset collection both unnecessary and inconvenient. On the one hand, due to the potential dynamics differences, the underlying states may not be well aligned for the dataset collected in source and target, e.g., $\tau^s=(s_{T=0}^s, s_{T=1}^s, s_{T=2}^s, ..., s_{T=n}^s)\neq \tau^t=(s_{T=0}^t, s_{T=1}^t, s_{T=2}^t, ..., s_{T=n}^t)$ especially for the latest states in trajectories. On the other hand, manually set the states for the source domain to match with the target domain may not be feasible or inconvenient for some practical cases, e.g., hard to set a simulator to a certain state, etc. Therefore, in our proposed method, the trajectory alignment is no longer required for the SI process. Specifically, we assume that a forward dynamics prediction model $\mathbf{F}$: $\hat{o^\prime}=\mathbf{F}(o,a;\theta)$ trained with source domain dataset $\mathcal{D}^s$ will accurately predict the next observations in the target domain dataset $\mathcal{D}^t$ only when the input dynamics embedding $\alpha$ of the target environment is accurate. Therefore we can construct a loss/score function for optimizing the dynamics embedding based on the forward prediction results in the target domain.

In our method, the forward dynamics prediction model $\mathbf{F}$: $\hat{o}^\prime=\mathbf{F}(o,a;\alpha)$, is trained with $\mathcal{D}^s$ and further applied on $\mathcal{D}^t$ to fit $\alpha$, assuming that the target data is within the distribution of source data. This can be achieved by increasing the randomization ranges in the source domain until satisfactory. For learning the forward prediction model $\mathbf{F}_\Phi$ and dynamics encoder $\mathbf{E}_\psi$, we have the following objectives:

\begin{align}
&\min_\Phi \mathbb{E}_{(o,a,o^\prime, \theta)\sim\mathcal{D}^s}[L_1(\mathbf{F}_\Phi(o,a;\mathbf{E}(\theta)), o^\prime)] \label{eq:forward_prediction} \\
&\min_\psi \mathbb{E}_{(o,a,o^\prime, \theta)\sim\mathcal{D}^s}[L_1(\mathbf{F}(o,a;\mathbf{E}_\psi(\theta)), o^\prime)+\lambda L_2(\mathbf{D}(\mathbf{E}_\psi(\theta)))] \label{eq:dynamics_encoder}
\end{align}

where both $L_1$ and $L_2$ are mean squared error (MSE) loss function in our experiments, and $\lambda$ is a trade-off coefficient for balancing the dynamics prediction performance and the reconstruction of dynamics parameters through encoding and decoding. This is corresponding to step (a) in Fig. <ref>.

The objective for optimizing the embedding $\alpha$ with the learned forward dynamics prediction function $\mathbf{F}_{\Phi^*}$ is:

\begin{equation}
\alpha^* = \arg \min_\alpha \mathbb{E}_{(o,a,o^\prime)\sim\mathcal{D}^t}[L_3(\mathbf{F}_{\Phi^*}(o,a;\alpha), o^\prime)]
\label{equ:embed}
\end{equation}

where $L_3$ is also MSE in our experiments. This is corresponding to step (c) in Fig. <ref>.

The universal policy $\pi$: $a=\pi(o;\alpha)$ in our method is trained in the source domain and applied for inference in the target domain. During the training process of the universal policy, observation noise and action noise are applied with randomized parameters, together with the randomization of dynamics parameters. The objective for learning the universal policy is:

\begin{equation}
\phi^* = \arg \max_\phi \mathbb{E}_{\theta}[J_\theta(\pi_\phi(\cdot|\mathbf{E}(\theta)))]
\label{equ:up}
\end{equation}

where $\phi$ are parameters of the policy $\pi$. This is corresponding to step (b) in Fig. <ref>.

Combining above objectives as in Eq. (<ref>) and (<ref>), our method is capable of achieving the optimal policy inference for any dynamics setting via optimized embedding SI:

\begin{equation}
\pi^* = \pi_{\phi^*}(\cdot|\alpha^*)
\end{equation}

which corresponds to step (d) in Fig. <ref>. The pseudo-code for the entire UPESI algorithm is shown as Alg. <ref>.

# (a). TRAIN DYNAMICS PREDICTION MODEL
Initialize data buffer $\mathcal{D}^s$, prediction model $\mathbf{F}_\Phi$ and encoder $\mathbf{E}_\psi$, and a pretrained policy $\pi_\text{ini}$
Execute policy $\pi_\text{ini}$ on robot to collect samples $(o,a,o^\prime)$ under dynamics $\theta$;
$\mathcal{D}^s=\mathcal{D}^s\bigcup (o,a,o^\prime, \theta)$

Sample a batch of data $\mathcal{D}^b$ from $\mathcal{D}^s$
Optimize $\mathbf{F}_\Phi$ and $\mathbf{E}_\psi$ with Eq. (<ref>) and Eq. (<ref>) respectively using $\mathcal{D}^b$

# (b). TRAIN UNIVERSAL POLICY
Fix encoder parameters $\mathbf{E}_\psi$
Initialize universal policy $\pi_\phi$
$o^\prime, r$=env.step($a$;$\theta$)
Update $\pi_\phi$ with Eq. (<ref>)

# (c). FIT DYNAMICS EMBEDDING
Fix dynamics prediction model $\mathbf{F}_\Phi$
Collect data buffer $\mathcal{D}^t$ in the target domain following the same way as the dynamics training data buffer $\mathcal{D}^s$
Bayesian optimize the embedding $\alpha$ with $\mathcal{D}^t$ as Eq. (<ref>);

# (d). INFERENCE IN A TARGET DOMAIN
Rollout the universal policy $\pi$ to generate action: $a=\pi(o;\alpha)$.
Universal Policy with Embedding System Identification (UPESI)

§ EXPERIMENT

§.§ Comparison Methods

The comparison involves the following methods:

(1). No DR. A conservative policy is trained in a certain environment with DR, as a comparison baseline.

(2). DR only. A general policy is trained with randomized environments on specified parameters and distributions.

(3). DR+UP (True). An adaptive universal policy is trained with true system parameters as additional inputs in randomized environments.

(4). DR+UP+SI. The policy is trained in the same manner as (3), but tested with a learned SI module to predict dynamics parameters from historical transitions.

(5). DR+UP+Encoding (BO). As the proposed method, an adaptive universal policy is trained with the embedding of system parameters as additional inputs in randomized environments, also with a BO process for configuring the embedding in new environment using the learned forward prediction model.

(6). DR+UP+Encoding (True). As an oracle for the proposed method, the universal policy is trained in the same way as (5), but the embedding is given by the true system parameters going through the learned encoder.

Method (1) works as a baseline for all other methods, which basically represents the optimal policy for a certain training environment and also the most conservative policy when testing in various environments. The comparison of (3) and (4) will imply the potential effects caused by the deficiency in SI module. The comparison of (5) and (6) will show the effects of embedding configuration based on BO when there are no true system parameters but only with samples of transition provided.

§.§ Experimental Setup

Different configurations for the environment InvertedDoublePendulum-v2.

Environments InvertedDoublePendulum-v2, HalfCheetah-v3 in OpenAI Gym MuJoCo are used for testing  the above methods.

For InvertedDoublePendulum-v2, the environment is randomized for its five important dynamics parameters, including damping, gravity, two geometry lengths and density, which are detailed in App. <ref> Tab. <ref>. The encoded latent space is chosen to have a dimension of 2 as the disentangled compact representation of system dynamics. Fig. <ref> shows visualization of three different configurations for the task scene with randomized parameters. For Halfcheetah-v3, 13 dynamics parameters are randomized as detailed in App. <ref> Tab. <ref>, with an embedding space of dimension 4. For training both the conservative or adaptive universal policies, we use the twin delayed deep deterministic policy gradient (TD3) algorithm [22], with 4 MLP layers for both the policy and the value networks. The dynamics networks and SI networks have 4 MLP layers. SI model predicts the system dynamics parameters with a stack of 5 frames of transitions, using the same training data as in the training of dynamics prediction models. In our experiments, each method is tested with three runs with different random seeds.

§.§ Experimental Results

Average learning performances of different methods for two environments. Methods include no DR (black dashed lines), DR only, DR with UP, DR with UP and encoding.

Fig. <ref> shows the average learning performances in terms of episode rewards in training for all methods on two environments. The policies trained without DR are evaluated in their specific training environments as optimal baselines. The shaded areas indicating the standard deviations. Each run takes 2000 or 25000 episodes of training the policies with TD3 algorithm for InvertedDoublePendulum-v2 and HalfCheetah-v3 respectively. The universal policies with either embedding SI or true system parameters significantly outperforms the DR baselines in two environments, although not as optimal as the no-DR policies since the episode reward is directly evaluated in the training environments (i.e. non-randomized for no-DR case but randomized for others).

Bayesian optimization on two-dimensional dynamics embedding space for the environment InvertedDoublePendulum-v2. Two cases with different dynamics parameters are displayed. The blue "crosses" are the query points in BO process, the red "stars" are the best points after BO, and the green "triangles" are the true values for dynamics embedding derived from the true system parameters. Note that the “triangles" and “stars" are almost overlapped.

For the embedding SI process, we take the InvertedDoublePendulum-v2 as an example. Fig. <ref> shows the BO process for configuring the embeddings of two sets of randomly sampled system parameters following the step (c) in our proposed method (as in Fig. <ref>), with the dynamics prediction model trained with 10000 episodes of policy rollouts and BO for each parameter configuration with 1000 episodes of data. After 500 iterations of BO for embedding configuration, the best query point is already very close to the true embedding values as shown in Fig. <ref>. Results are similar for HalfCheetah-v3 just with a higher dimensional embedding. 2000 episodes of data are leveraged for training the dynamics prediction model and 100 episodes of data are used for BO in identifying each new testing environment.

Test performances of different methods for two environments. The white dots are the mean values, and the black bold lines indicate the 25%-75% quantiles.

tableTest performances: mean and standard deviations of episode rewards for different methods on two environments.

Env     Method      Episode Reward

6*pendulum     No DR     $5168.9\pm{4324.4}$
DR only     $2791.1\pm{2960.8}$
DR+UP (True)     $7586.8\pm{3435.7}$
DR+UP+SI     $2153.9\pm{3741.0}$
DR+UP+Encoding (BO)     $\mathbf{7303.9\pm{3574.2}}$
DR+UP+Encoding(True)     $7474.0\pm{3454.7}$
6*halfcheetah     No DR     $3213.6\pm{2917.9}$
DR only     $7594.2\pm{2571.1}$
DR+UP (True)    $8956.7\pm{2279.9}$
DR+UP+SI     $5679.7\pm{2994.4}$
DR+UP+Encoding (BO)     $\mathbf{8051.5\pm{2646.6}}$
DR+UP+Encoding(True)     $8194.5\pm{2385.5}$

Fig. <ref> and Tab. <ref> display the test performances of different methods on randomly sampled system parameters within the same distributions of DR, as the step (d) in Fig. <ref>. We show the episode reward distributions of different methods as a violin plot in Fig. <ref>. Policies for each method are tested on 10 randomly sampled environment dynamics for 10 episodes each, as a total of 100 episodes for each method. We can see that although the method without DR can achieve good rewards in training since there is no randomized dynamics but a fixed one, in test cases it has a heavy allocation for both the head and tail on the performance distribution. The reason is that the lack of randomized dynamics show some unseen cases incapable of being handled well by the policy. The method with only DR also does not work well because of the hardness in training an optimal policy in randomized environments, which is testified to be an ill-posed problem. For both environments, the heavy tails in no-DR method are greatly alleviated for methods with UP, due to its awareness of system dynamics. However, the true dynamics parameters of the testing environment are usually not accessible, which requires a SI module to configure. Our experiments show that the SI process with a stack of 5 frames of historical transitions (see similar settings in [1]) is not capable of providing an accurate estimation of the system dynamics,  which severely degrades the performances for UP methods. For both environments, our proposed method with embedding SI using BO process shows advantageous performances over other methods, even as good as the one with true system parameters.

§ CONCLUSIONS AND DISCUSSIONS

We propose to optimize the dynamics embedding rather than directly configuring the system parameters from historical transitions, and demonstrate its advantageous performances over standard DR and UP methods for a general domain transfer setting, with some primary tests on both a low-dimensional and a high-dimensional simulated environments. The deficiency for DR policies of not being able to achieve optimal actions and the difficulties within normal SI methods for directly configuring the system parameters are revealed. Our future work involves extending the current framework to more complex robot learning tasks, as well as its application on sim-to-real transfer problem, which is a subset of the domain transfer in general. The difficulty of UP learning due to the expanded input spaces with high-dimensional system parameters are not exposed in current experiments, therefore some higher-dimensional tasks will be investigated as well.

§ APPENDIX

§.§ Theoretical Proofs

Given the same dataset distribution $\rho^s(s,a,s')=\rho^t(s,a,s')$, the difference between the estimated dynamics distribution $p^s$ and the true distribution $p^t$ can be characterized as:

\begin{align}
&\mathbb{D}_\text{KL}[p^s(\theta|s,a,s')||p^t(\theta|s,a,s')] \\
% &=\mathbb{D}_\text{KL}[\rho^s(s,a,\theta,s')||\rho^t(s,a,\theta,s')] - \mathbb{D}_\text{KL}[\rho^s(s,a,s')||\rho^t(s,a,s')] \\
\end{align}

\begin{align}
\text{RHS}&=\mathbb{D}_\text{KL}[\rho^s(s,a,\theta,s')||\rho^t(s,a,\theta,s')]\\
&-\mathbb{D}_\text{KL}[\rho^s(s,a,s')||\rho^t(s,a,s')] \\
&=\int_{\mathcal{S}\times\mathcal{A}\times\Theta\times\mathcal{S}}\rho^s(s,a,\theta,s')\log(\frac{\rho^s(s,a,\theta,s')}{\rho^t(s,a,\theta,s')} \\ &\times\frac{\rho^t(s,a,s')}{\rho^s(s,a,s')})\text{d}s\text{d}a\text{d}\theta\text{d}s' \\
&=\int_{\mathcal{S}\times\mathcal{A}\times\Theta\times\mathcal{S}}\rho^s(s,a,\theta,s')\log\frac{p^s(\theta|s,a,s')}{p^t(\theta|s,a,s')}\text{d}s\text{d}a\text{d}\theta\text{d}s' \\
&=\int_{\Theta}p^s(\theta|s,a,s')\log\frac{p^s(\theta|s,a,s')}{p^t(\theta|s,a,s')}\text{d}\theta \\
&=\mathbb{D}_\text{KL}[p^s(\theta|s,a,s')||p^t(\theta|s,a,s')] = \text{LHS}.
\end{align}

4 The distance of the distribution from the forward dynamics prediction and the true distribution can be formulated with the KL-divergence, it thus follows:

\begin{align}
% &\mathbb{D}_\text{KL}[f^s(s'|s,a, \theta)||f^t(s'|s,a, \theta)]\nonumber \\
% &\approx \mathbb{D}_\text{KL}[\rho^s(s,a,\theta,s')||\rho^t(s,a,\theta,s')]

\mathbb{D}_\text{KL}[f^s||f^t]\approx \mathbb{D}_\text{KL}[\rho^s||\rho^t] -\mathbb{D}_\text{KL}[\rho^s(s,a,\theta)||\rho^t(s,a,\theta)]
\end{align}

where $f^{s,t}$ are shorten for $f^{s,t}(s'|s,a, \theta)$, and $\rho^{s,t}$ are shorten for $\rho^{s,t}(s,a,\theta,s')$.

Similar as the proof of Lemma <ref>, we have,

\begin{align}
\text{RHS}&=\mathbb{D}_\text{KL}[\rho^s(s,a,\theta,s')||\rho^t(s,a,\theta,s')] \\
&- \mathbb{D}_\text{KL}[\rho^s(s,a,\theta)||\rho^t(s,a,\theta)] \\
&=\int_{\mathcal{S}\times\mathcal{A}\times\Theta\times\mathcal{S}}\rho^s(s,a,\theta,s')\log(\frac{\rho^s(s,a,\theta,s')}{\rho^t(s,a,\theta,s')} \\
&\times\frac{\rho^t(s,a,\theta)}{\rho^s(s,a,\theta)})\text{d}s\text{d}a\text{d}\theta\text{d}s' \\
&\approx\int_{\mathcal{S}\times\mathcal{A}\times\Theta\times\mathcal{S}}\rho^s(s,a,\theta,s')\log\frac{f^s(s'|s,a,\theta)}{f^t(s'|s,a,\theta)}\text{d}s\text{d}a\text{d}\theta\text{d}s' \\
&=\int_{\Theta}f^s(s'|s,a,\theta)\log\frac{f^s(s'|s,a,\theta)}{f^t(s'|s,a,\theta)}\text{d}s' \\
&=\mathbb{D}_\text{KL}[f^s(s'|s,a,\theta)||f^t(s'|s,a,\theta)] = \text{LHS}
\end{align}

The optimization of forward dynamics prediction is increasing the lower bound of the optimization objective for improving the estimated dynamics, i.e.,

\begin{equation}
\mathbb{D}_\text{KL}[p^s||p^t]\geq  \mathbb{D}_\text{KL}[f^s||f^t]
\end{equation}

where $p^{s,t}$ are $p^{s,t}(\theta|s,a,s')$ and $f^{s,t}$ are $f^{s,t}(s'|s,a,\theta)$.

With Lemma <ref> and <ref>, we have,

\begin{align}
\text{LHS}&=\mathbb{D}_\text{KL}[\rho^s(s,a,\theta,s')||\rho^t(s,a,\theta,s')]\\
&=\mathbb{D}_\text{KL}[\rho^s(s,a,\theta,s')||\rho^t(s,a,\theta,s')] \\
&- \mathbb{D}_\text{KL}[\rho^s(s,a,\theta)||\rho^t(s,a,\theta)] + \mathbb{D}_\text{KL}[\rho^s(s,a,\theta)||\rho^t(s,a,\theta)]\\
&\approx \mathbb{D}_\text{KL}[f^s(s'|s,a,\theta)||f^t(s'|s,a,\theta)] \\
&+ \mathbb{D}_\text{KL}[\rho^s(s,a,\theta)||\rho^t(s,a,\theta)]\\
&\geq \mathbb{D}_\text{KL}[f^s(s'|s,a,\theta)||f^t(s'|s,a,\theta)]
\end{align}

The last line of above proof is due to that $\mathbb{D}_\text{KL}[\rho^s(s,a,\theta)||\rho^t(s,a,\theta)]$ is non-negative.

§ RANDOMIZED PARAMETERS

tableRandomized dynamics parameters (5-dimensional) for InvertedDoublePendulum-v2

Variable     Range

damping     [0.02, 0.3]
gravity     [8.5, 11.0]
length1    [0.3, 0.9]
length2     [0.3, 0.9]
density     [0.5, 1.5]

tableRandomized dynamics parameters (13-dimensional) for HalfCheetah-v3

Variable     Range

gravity     [5.5, 14.0]
bthigh damping     [3.0, 9.0]
bshin damping     [1.5, 7.5]
bfoot damping     [1.0, 5.0]
fthigh damping     [1.5, 7.5]
fshin damping     [1.0, 5.0]
ffoot damping     [0.2, 2.8]
bthigh stiffness     [100, 380]
bshin stiffness     [20, 340]
bfoot stiffness     [10, 230]
fthigh stiffness     [20, 340]
fshin stiffness     [20, 220]
ffoot stiffness     [10, 110]

[1]
W. Yu, J. Tan, C. K. Liu, and G. Turk, “Preparing for the unknown: Learning a
universal policy with online system identification,” arXiv preprint
arXiv:1702.02453, 2017.

[2]
X. B. Peng, M. Andrychowicz, W. Zaremba, and P. Abbeel, “Sim-to-real transfer
of robotic control with dynamics randomization,” in 2018 IEEE
international conference on robotics and automation (ICRA).1em plus
0.5em minus 0.4emIEEE, 2018, pp. 3803–3810.

[3]
O. M. Andrychowicz, B. Baker, M. Chociej, R. Jozefowicz, B. McGrew,
J. Pachocki, A. Petron, M. Plappert, G. Powell, A. Ray et al.,
“Learning dexterous in-hand manipulation,” The International Journal
of Robotics Research, vol. 39, no. 1, pp. 3–20, 2020.

[4]
R. S. Sutton, A. G. Barto et al., Introduction to reinforcement
learning.1em plus 0.5em minus 0.4emMIT press Cambridge, 1998,
vol. 135.

[5]
Y. Duan, X. Chen, R. Houthooft, J. Schulman, and P. Abbeel, “Benchmarking deep
reinforcement learning for continuous control,” in International
conference on machine learning.1em plus 0.5em minus 0.4emPMLR,
2016, pp. 1329–1338.

[6]
H. Dong, H. Dong, Z. Ding, S. Zhang, and Chang, Deep Reinforcement
Learning.1em plus 0.5em minus 0.4emSpringer, 2020.

[7]
W. Yu, V. C. Kumar, G. Turk, and C. K. Liu, “Sim-to-real transfer for biped
locomotion,” arXiv preprint arXiv:1903.01390, 2019.

[8]
E. Valassakis, Z. Ding, and E. Johns, “Crossing the gap: A deep dive into
zero-shot sim-to-real transfer for dynamics,” arXiv preprint
arXiv:2008.06686, 2020.

[9]
J. Tobin, R. Fong, A. Ray, J. Schneider, W. Zaremba, and P. Abbeel, “Domain
randomization for transferring deep neural networks from simulation to the
real world,” in 2017 IEEE/RSJ international conference on intelligent
robots and systems (IROS).1em plus 0.5em minus 0.4emIEEE,
2017, pp. 23–30.

[10]
S. Zhu, D. Surovik, K. Bekris, and A. Boularias, “Efficient model
identification for tensegrity locomotion,” in 2018 IEEE/RSJ
International Conference on Intelligent Robots and Systems (IROS).1em
plus 0.5em minus 0.4emIEEE, 2018, pp. 2985–2990.

[11]
Y. Chebotar, A. Handa, V. Makoviychuk, M. Macklin, J. Issac, N. Ratliff, and
D. Fox, “Closing the sim-to-real loop: Adapting simulation randomization
with real world experience,” in 2019 International Conference on
Robotics and Automation (ICRA).1em plus 0.5em minus 0.4emIEEE,
2019, pp. 8973–8979.

[12]
S. James, A. J. Davison, and E. Johns, “Transferring end-to-end visuomotor
control from simulation to real world for a multi-stage task,” in
Conference on Robot Learning.1em plus 0.5em minus 0.4emPMLR, 2017, pp. 334–343.

[13]
A. Prakash, S. Boochoon, M. Brophy, D. Acuna, E. Cameracci, G. State,
O. Shapira, and S. Birchfield, “Structured domain randomization: Bridging
the reality gap by context-aware synthetic data,” in 2019
International Conference on Robotics and Automation (ICRA).1em plus
0.5em minus 0.4emIEEE, 2019, pp. 7249–7255.

[14]
B. Mehta, M. Diaz, F. Golemo, C. J. Pal, and L. Paull, “Active domain
randomization,” in Conference on Robot Learning.1em plus 0.5em
minus 0.4emPMLR, 2020, pp. 1162–1176.

[15]
M. Wang and W. Deng, “Deep visual domain adaptation: A survey,”
Neurocomputing, vol. 312, pp. 135–153, 2018.

[16]
A. K. Tanwani, “Domain-invariant representation learning for sim-to-real
transfer,” arXiv preprint arXiv:2011.07589, 2020.

[17]
W. Yu, C. K. Liu, and G. Turk, “Policy transfer with strategy optimization,”
arXiv preprint arXiv:1810.05751, 2018.

[18]
W. Yu, J. Tan, Y. Bai, E. Coumans, and S. Ha, “Learning fast adaptation with
meta strategy optimization,” IEEE Robotics and Automation Letters,
vol. 5, no. 2, pp. 2950–2957, 2020.

[19]
J. Mockus, Bayesian approach to global optimization: theory and
applications.1em plus 0.5em minus 0.4emSpringer Science &
Business Media, 2012, vol. 37.

[20]
N. Hansen, A. Ostermeier, and A. Gawelczyk, “On the adaptation of arbitrary
normal mutation distributions in evolution strategies: The generating set
adaptation.” in ICGA.1em plus 0.5em minus 0.4emCiteseer, 1995, pp. 57–64.

[21]
C. Yang, X. Ma, W. Huang, F. Sun, H. Liu, J. Huang, and C. Gan, “Imitation
learning from observations by minimizing inverse dynamics disagreement,” in
Advances in Neural Information Processing Systems, 2019, pp. 239–249.

[22]
S. Fujimoto, H. Hoof, and D. Meger, “Addressing function approximation error
in actor-critic methods,” in International Conference on Machine
Learning.1em plus 0.5em minus 0.4emPMLR, 2018, pp. 1587–1596.