# Two-orbital Hubbard model vs spin $S=1$ Heisenberg model: studies on
clusters

R. Lemański, J. Matysiak

(Received June 15, 2018, in final form August 13, 2018)

###### Abstract

We perform exact numeric calculations for the two-orbital Hubbard model on the
four-site cluster. In the limit of large on-site coupling the model becomes
equivalent to the spin $S=1$ Heisenberg model. By comparing energy spectra of
these two models, we quantified the range of interaction parameters for which
the Heisenberg model satisfactorily reproduces the two-orbital Hubbard model.
Then we examined how the spectrum evolves when we are outside of this region,
focusing especially on checking of how it is modified when various ways of
interatomic hoppings of electrons between different orbitals are taken into
account. We finally show how these modifications affect the dependence of
specific heat on temperature.

Key words: multi-orbital Hubbard model, Heisenberg model, magnetic molecules

PACS: 31.15.vq, 75.10.Jm, 75.10.Pq

###### Abstract

Íàìè ïðîâåäåíî òîí èñëîâ ðîçðàõóíêè äëÿ äâîîðáòàëüíî¿ ìîäåë Õàááàðäà íà
îòèðèâóçëîâîìó êëàñòåð. Â ãðàíèö âåëèêî¿ îäíîâóçëîâî¿ âçàìîä¿ öÿ ìîäåëü ñòà
åêââàëåíòíîþ ñïí $S=1$ ìîäåë Ãàéçåíáåðãà. Ç ïîðâíÿííÿ åíåðãåòèíèõ ñïåêòðâ öèõ
äâîõ ìîäåëåé íàìè âèçíàåíî êëüêñíî äàïàçîí çíàåíü ïàðàìåòðâ âçàìîä¿, äëÿ ÿêîãî
ìîäåëü Ãàéçåíáåðãà çàäîâëüíî âäòâîðþ äâîîðáòàëüíó ìîäåëü Õàááàðäà. Ïîòì ìè
ðîçãëÿíóëè çìíè ñïåêòðó, êîëè ìè çíàõîäèìîñÿ çà ìåæàìè öüîãî äàïàçîíó,
çîñåðåäæóþèñü, çîêðåìà, íà ïåðåâðö òîãî, ÿê âí ìîäèôêóòüñÿ, êîëè âðàõîâóþòüñÿ
ðçí ñïîñîáè ìæàòîìíèõ ïåðåñêîêâ åëåêòðîíâ ìæ ðçíèìè îðáòàëÿìè. Íàñàìêíåöü,
íàìè ïîêàçàíî ÿê ö çìíè âïëèâàþòü íà òåìïåðàòóðíó çàëåæíñòü ïèòîìî¿
òåïëîìíîñò.

Ключов слова: áàãàòîîðáòàëüíà ìîäåëü Õàááàðäà, ìîäåëü Ãàéçåíáåðãà, ìàãíòí
ìîëåêóëè

##  1 Introduction

Multi-orbital Hubbard models (MOHMs) are suitable for studying correlated
electron materials with orbital degeneracies, such as transition metal
compounds [1, 2, 3]. In particular, they seem to be relevant in the
theoretical analysis of single molecular magnets (SMMs). However, magnetic
molecules have been usually described by phenomenological spin model
Hamiltonians with the dominant Heisenberg interaction term (e.g., [4, 5]).
This is because the low-energy excited states observed in these systems are
often similar to the energy spectrum of the Heisenberg model. Then, attempts
were made to describe these systems from the first principles using the DFT
method and to determine the exchange parameters entering the Heisenberg model
from the comparison of energies of various magnetic configurations [6]. The
resulting energy spectrum obtained in these calculations slightly deviates
from the spectrum of the simplest version of the Heisenberg model, but it can
be well tuned to the spectrum of the Falicov-Kimball model [6]. Unfortunately,
the calculations reported in [6] did not take into account the spin flip
processes, which made the exchange integrals too large to be compatible with
the experimental data.

It has been recently proposed to describe SMMs using the MOHM in the large
interaction limit combined with DFT calculations [7]. In this work, the values
of exchange integrals were obtained within the perturbation theory. They
turned out to be several times smaller than those obtained in [6], thus better
suited to experimental data. The point is that in [7], the DFT was not used to
calculate the energies of the system with different magnetic configurations,
but it was focused on determining the local Coulomb and Hund interactions and
the amplitudes of inter-ionic hoppings of electrons. Then, the exchange
integral was determined by putting these data into the formula obtained from
the perturbation theory in the limit of large values of the local Coulomb
interaction, when the MOHM reduces to the Heisenberg model. In fact, the
transition from the Hubbard model to the spin model by applying the
perturbation theory is a commonly used procedure in correlated electron
systems not only in relation to magnetism, but also relative to other
phenomena, such as superconductivity [8, 9, 10, 11, 12, 13, 14].

Although the approach used in [7] improved the procedure of finding the
integrals of the exchange, still the system was described only using the
Heisenberg model. Therefore, nothing was known what exactly we could gain when
instead of the Heisenberg model we examine the system using the MOHM. Indeed,
the procedure applied there did not allow one to determine the range of the
model parameters, beyond which the energy spectra of the MOHM and the
Heisenberg model clearly differ from each other. Moreover, regardless of the
orbitals between which the electrons jump from site to site, the perturbation
theory always leads to the Heisenberg model, modifying only the value of the
exchange constant. In other words, the perturbative theory does not
distinguish whether the electron jumps between the same or between different
orbitals (‘hybridization’).

To find out what is the advantage of the description of a SMM using the
microscope model MOHM over the phenomenological description using the
Heisenberg model, and in particular to see how the energy spectra of these two
models differ from each other, here we analyze the four-site ring with two
orbitals per site. Such a small ring was chosen to make it possible to perform
exact numerical calculations efficiently and rather easily. Therefore, our
analysis does not focus on a particular material, but rather on capturing a
role played by various hopping amplitudes of electrons, when the system goes
out of the range where the second order perturbation theory is justified. For
a better illustration of our results, in this work we also present a set of
curves of the specific heat versus temperature for various relationships
between amplitudes of electron hoppings between different orbitals.

Figure 1: Labelling of the hopping amplitudes $t_{mm^{\prime}}$ between the
sites $i$ and $j$ in the TOHM. $m,m^{\prime}$ indicate orbital indices.

##  2 Two-orbital Hubbard model (TOHM)

Studies of the two-orbital case require an inclusion into the single-orbital
Hubbard model of additional on-site couplings between various orbitals: the
direct and exchange Coulomb interactions that ensure the fulfillment of the
Hund’s rule. An additional complication that arises here is the possibility of
electron jumping between more than one orbital and also between different
orbitals (mixing terms or ‘hybridization’) belonging to the neighboring sites
(see figure 1). Taking into account all these factors ensures not only an
increase of the number of model parameters, but above all a significant
increase of the number of states of the system, which obviously complicates
diagonalization of the Hamiltonian.

The Hamiltonian $H_{\text{HM}}$ of the TOHM that we study in this paper is as
follows:

$\displaystyle H_{\text{HM}}$ $\displaystyle=H_{0}+H_{1}\,,$ (2.1)
$\displaystyle H_{0}$
$\displaystyle=U\sum\limits_{im}n_{im\downarrow}n_{im\uparrow}+\frac{1}{2}\sum\limits_{i,m\neq
m^{\prime},\sigma}\left(U^{\prime}n_{im\sigma}n_{im^{\prime}\bar{\sigma}}+U^{\prime\prime}n_{im\sigma}n_{im^{\prime}\sigma}\right)$
$\displaystyle\quad+\frac{1}{2}\sum\limits_{i,m\neq
m^{\prime},\sigma}\left(Jc_{im\sigma}^{\dagger}c_{im^{\prime}\bar{\sigma}}^{\dagger}c_{im\bar{\sigma}}c_{im^{\prime}\sigma}+Jc_{im\sigma}^{\dagger}c_{im\bar{\sigma}}^{\dagger}c_{im^{\prime}\bar{\sigma}}c_{im^{\prime}\sigma}\right),$
$\displaystyle H_{1}$ $\displaystyle=\sum\limits_{i\neq
j,m,m^{\prime},\sigma}t_{mm^{\prime}}c_{im\sigma}^{\dagger}c_{jm^{\prime}\sigma}\,,$

where $i$ and $j$ denote nearest-neighbour sites, $m,m^{\prime}$ label
orbitals and $\sigma,\bar{\sigma}$ label spins of electrons
($\bar{\sigma}=-\sigma$). $U$, $U^{\prime}$ and $U^{\prime\prime}$ describe
the Coulomb type on-site interactions between two electrons: $U$ — on the same
orbital and $U^{\prime}$ ($U^{\prime\prime}$) — on different orbitals with
opposite (parallel) spins, respectively. $J$ represents the on-site exchange
coupling, but it also enters the interaction constants in the TOHM. Indeed, we
take $U^{\prime}=U-2J$, $U^{\prime\prime}=U-3J$, which is valid for $T_{2g}$
orbitals in an octahedral crystal field [15].

The Hamiltonian (2.1) preserves the total magnetisation. Therefore, in the
beginning we consider only the states with total magnetisation equal to zero,
which constitutes the largest subspace. The total number of these states is
4900. It is helpful to analyse first the spectrum of the TOHM in the atomic
limit, i.e., when all $t_{mm^{\prime}}=0$ (the lowest part of the spectrum is
shown in figure 2). We have chosen the parameters in our calculations to be
$U=10$ eV, $J=0.5$ eV. The spectrum consists of 26 degenerate energy levels.
Then, we are interested in the region of small values of the hopping constants
$t_{mm^{\prime}}$, and it is understood that the overall structure of the
spectrum in figure 2 is preserved, but each of the degenerate levels splits in
a specific way when $t_{mm^{\prime}}$ becomes non-zero. Here, we focus mainly
on finding the splitting of the lowest energy level.

Figure 2: A part of energy spectrum of the half-filled TOHM for the four-site
ring in the atomic limit. The interaction parameters are $U=10$ eV, $J=0.5$
eV, $S_{z}=0$.

To obtain the exact spectrum of the model (2.1) for our small system, we use
the exact diagonalization method.

##  3 Computation details

In the atomic limit, the eigenstates of the Hamiltonian (2.1) for a single
site together with their representations and energies are shown in table 1. As
indicated in the first column, they are also the eigenstates of the total spin
$S$ and its $z$ component $S_{z}$. The full Hilbert space base for the whole
4-site system is formed from the tensor products of the states from table 1.

For a site occupied by two electrons, the lowest energy refers to the
eigenstates 1–3 from table 1. However, obviously, after one of the electrons
jumps from one site to another, the first site will contain only one electron
and the other one three electrons. Then, the lowest energy eigenstate for the
site with one electron is one of the eigenstates 7–10 from table 1., and for
the site with three electrons, it is one of the eigenstates 11–14 from table
1.

Table 1: Eigenstates and eigenvalues of the model (2.1) for a single site with two orbitals per site. Sections separated by double horizontal lines correspond to the different numbers of electrons $n$ on the site. nr | $\ket{n,S,S_{z}}$ | State representation | Energy
---|---|---|---
1 | $\ket{2,1,1}$ | $\quantity(\uparrow|\uparrow)$ | $U-3J$
2 | $\ket{2,1,0}$ | $\dfrac{1}{\sqrt{2}}\quantity\Big(\quantity(\uparrow|\downarrow)+(\downarrow|\uparrow))$ | $U-3J$
3 | $\ket{2,1,-1}$ | $\quantity(\downarrow|\downarrow)$ | $U-3J$
4 | $\ket{2,0,0}$ | $\dfrac{1}{\sqrt{2}}\quantity\Big(\quantity(\uparrow|\downarrow)-(\downarrow|\uparrow))$ | $U-J$
5 | $\ket{2,0_{a},0}$ | $\dfrac{1}{\sqrt{2}}\quantity\Big(\quantity(\updownarrows|-)+(-|\updownarrows))$ | $U-J$
6 | $\ket{2,0_{b},0}$ | $\dfrac{1}{\sqrt{2}}\quantity\Big(\quantity(\updownarrows|-)-(-|\updownarrows))$ | $U+J$
7 | $\ket{1,1/2_{a},1/2}$ | $\quantity(\uparrow|-)$ | $0$
8 | $\ket{1,1/2_{a},-1/2}$ | $\quantity(\downarrow|-)$ | $0$
9 | $\ket{1,1/2_{b},1/2}$ | $\quantity(-|\uparrow)$ | $0$
10 | $\ket{1,1/2_{b},-1/2}$ | $\quantity(-|\downarrow)$ | $0$
11 | $\ket{3,1/2_{a},1/2}$ | $\quantity(\uparrow|\updownarrows)$ | $3U-5J$
12 | $\ket{3,1/2_{a},-1/2}$ | $\quantity(\downarrow|\updownarrows)$ | $3U-5J$
13 | $\ket{3,1/2_{b},1/2}$ | $\quantity(\updownarrows|\uparrow)$ | $3U-5J$
14 | $\ket{3,1/2_{b},-1/2}$ | $\quantity(\updownarrows|\downarrow)$ | $3U-5J$
15 | $\ket{0,0,0}$ | $\quantity(-|-)$ | $0$
16 | $\ket{4,0,0}$ | $\quantity(\updownarrows|\updownarrows)$ | $6U-10J$

In the case of a ring with 4 sites and 2 orbitals per site, we have the total
number of orbitals $r=8$. Then, at half filling, i.e., when the number of
electrons is equal to the namber of orbitals, the numbers of the states that
are needed to diagonalize the Hamiltonian are: 12870 for the TOHM, but only
$3^{4}=81$ for the Heisenberg model.

###  3.1 Perturbation theory

In the perturbation theory, $H_{1}$ is taken as the perturbation and $H_{0}$
is the Hamiltonian for which the solution is known. Contrary to the single
orbital model, in the multi-orbital case, the states with one electron per
orbital and with different spin configurations have, in general, different
energies even in the limit when $t_{mm^{\prime}}=0$ for all hopping
amplitudes. Then, in order to obtain an effective Hamiltonian with only spin
degrees of freedom, we restrict our considerations to a single subspace, where
all unperturbed states are degenerate. In particular, we focus on the lowest
energy subspace.

Here, we only provide the final result without going through the calculations.
A detailed derivation of the formulae reported below can be found in [7]. The
result, which is valid only for the subspace with the lowest energy, is that
the effective Hamiltonian turns out to be the standard Heisenberg Hamiltonian
with the effective Heisenberg constant $\Gamma_{\text{eff}}$:

$\displaystyle H_{\text{eff}}$ $\displaystyle=\Gamma_{\text{eff}}\sum_{i\neq
j}\quantity(S_{i}S_{j}-4),$ (3.1) $\displaystyle\Gamma_{\text{eff}}$
$\displaystyle=\dfrac{4\absolutevalue{t_{\text{eff}}}^{2}}{U+J}.$

$\Gamma_{\text{eff}}$ is determined by a single effective hopping constant
$t_{\text{eff}}$, which in a simple way depends on an initial hopping
constants $t_{mm^{\prime}}$:

$t_{\text{eff}}=\dfrac{\sqrt{\sum_{m,m^{\prime}}t_{mm^{\prime}}^{2}}}{2}.$
(3.2)

##  4 The lowest part of energy spectra

We start the presentation of our results by specifying the conditions for the
model parameters for which the perturbation theory is justified. Our
calculations are illustrated in figure 3 where the differences between the
energies of the ground states, of the first excited states and of the second
excited states were determined from the exact diagonalization and by using the
Heisenberg model with the effective exchange integral calculated from the
perturbation theory. From figure 3 it is clear that these differences increase
more steeply than linearly with an increase of the ratio
$t_{\text{eff}}/(U+J)$. Interestingly, the energies of the ground state differ
the most, almost twice less differ the energies of the first excited state and
even less do the ones of the second excited state.

Obviously, the limit value of $t_{\text{eff}}/(U+J)$ above which the
Heisenberg model ceases to describe this system adequately, can be fixed
arbitrarily. Here, we assumed that the difference between the state energies
should not exceed 0.03 eV. This corresponds to $t_{\text{eff}}/(U+J)$ value
equal to 0.05.

Figure 3: (Colour online) Comparison of the difference between the energies
calculated for the TOHM and for the Heisenberg model as a function of
$t_{\text{eff}}/(U+J)$ ($S_{z}=0$, $t_{12}=t_{21}=0$, $t_{11}=t_{22}$). The
data are shown for the lowest three eigenvalues. The interaction parameters
are: $U=10$ eV, $J=0.5$ eV.

Next, we present the evolution of the lowest part of the energy spectrum of
the TOHM depending on the ratio $t_{12}/t_{11}$ between the hopping amplitudes
of electrons, while the constant value of $t_{\text{eff}}$ is maintained. In
the case presented in figure 4, we have chosen $t_{\text{eff}}=0.5$. When
$t_{\text{eff}}$ is constant, then the perturbation theory gives the energy
spectrum typical of the Heisenberg model, which does not depend on the ratio
$t_{12}/t_{11}$ (see the first column in figure 4). However, the exact
solution shows that the energy spectrum of the TOHM clearly depends on the
ratio $t_{12}/t_{11}$, which is illustrated by the other columns in figure 4.
Thus, the theory of perturbation is too imprecise to distinguish the cases
with different ways of electron jumping between neighboring magnetic ions. In
fact, in figure 4, we present only the energy levels with the total $S_{z}=0$.

Figure 4: (Colour online) Comparison of the energy spectrum of the Heisenberg
model with the lowest parts of the energy spectra of the TOHM for the set of
the following relationships between the hopping amplitudes:
$t_{12}/t_{11}=0,0.2,0.4,0.6,0.8,1$. Here, $t_{21}=t_{12}$, $t_{22}=0$, $U=10$
eV and $J=0.5$ eV, $S_{z}=0$. In all these cases, the same
$t_{\text{eff}}=0.5$ is kept. For the Heisenberg model, the exchange constant
$\Gamma_{\text{eff}}=4t_{\text{eff}}^{2}/(U+J)$.

The difference between the exact solution and that based on the Heisenberg
model is clearly seen by comparing the energy of the second excited state
depending on the ratio $t_{12}/t_{11}$ for several different values of
$t_{\text{{eff}}}$, as it is shown in figure 5. For that particular state, the
difference decreases with an increasing $t_{12}/t_{11}$ ratio. The same trend
takes place for all $t_{\text{{eff}}}$ values, though the magnitude of the
difference increases exponentially.

Figure 5: (Colour online) Comparison of the difference between energies of the
second excited state calculated exactly for the MOHM and for the Heisenberg
model as a function of $t_{12}/t_{11}$, $S_{z}=0$. The interaction parameters
are: $U=10$ eV, $J=0.5$ eV. Figure 6: (Colour online) Specific heat as a
function of temperature for the Heisenberg model with
$\Gamma_{\text{eff}}=4t_{\text{eff}}/(U+J)=0.095$ eV and for various versions
of the MOHM model (for all of them $t_{\text{eff}}=0.5$ eV is kept). The
curves for the MOHM model differ by ratios of the hopping constants
$t_{12}/t_{11}$ ($t_{21}=t_{12}$, $t_{22}=0$), as indicated on the plot
legend. The parameters of the model are: $U=10$ eV, $J=0.5$ eV. Figure 7:
(Colour online) The difference between specific heat of the Heisenberg model
with $\Gamma=4t_{\text{eff}}/(U+J)=0.095$ eV and the different cases of the
MOHM model (for all of them $t_{\text{eff}}=0.5$ eV is kept). The curves for
the MOHM model differ by ratios of the hopping constants $t_{12}/t_{11}$
($t_{21}=t_{12}$, $t_{22}=0$), as indicated on the plot legend. The parameters
of the model are: $U=10$ eV, $J=0.5$ eV.

It is obvious that the distribution of the lowest energy levels determines the
dependence of specific heat on temperature as it is illustrated in figures 6
and 7. Figure 6 illustrates the $C_{v}(T)$ changes with an increase of the
$t_{12}/t_{11}$ ratio for the fixed value $t_{\text{eff}}=0.5$ for the TOHM,
as well as $C_{v}(T)$ for the Heisenberg model for comparison. In turn, in
figure 7, the differences between the values of $C_{v}(T)$ calculated for the
TOHM and the Heisenberg model with an increase of the ratio $t_{12}/t_{11}$
for the same value of $t_{\text{eff}}=0.5$ are displayed. Additionally, table
2 contains numeric values of the maxima of the curves from figure 7. Based on
the drawings given in figures 6 and 7, and especially in the table 2, it can
be supposed that from a course of the dependency $C_{v}$ you can draw
conclusions about the degree of hybridization expressed by the ratio
$t_{12}/t_{11}$.

Table 2: Maxima of the differences between the values of specific heat calculated for the TOHM and for the Heisenberg model calculated for various $t_{12}/t_{11}$ ratios. $t_{12}/t_{11}$ ratio | 0 | 0.2 | 0.4 | 0.6 | 0.8 | 1
---|---|---|---|---|---|---
maximum $\Delta C_{v}$ [meV/K] | 0.051 | 0.049 | 0.046 | 0.036 | 0.024 | 0.018

##  5 Summary and conclusions

We carried out exact numerical calculations for the four-site ring with two
orbitals and two electrons per site, using the TOHM and we compared the
results with those obtained for the respective Heisenberg model. On the basis
of this comparison, we estimated the range of the interaction parameters, and,
more specifically, the range of the $t_{\text{eff}}/(U+J)$ ratio, for which
the TOHM is well represented by the Heisenberg model.

However, the main achievement of our work is a demonstration that the energy
spectrum of the system depends in a significant way on how an electron can
jump from one site to another. Indeed, the case when an electron can only jump
between the same type of orbital is different from the case when the electron
can jump both between the same type and between different types of orbitals.
This difference manifests itself through a modification of the energy
spectrum, which evidently evolves with an increase in the degree of
hybridization between different orbitals. This, in turn, brings a change,
e.g., into the dependence of specific heat on temperature $C_{v}(T)$. It can,
therefore, be supposed that precise measurements of the specific heat would
enable us to determine the degree of hybridization between different orbitals
located on neighboring sites. However, to make this possible, the ratio of the
effective amplitude $t_{\text{eff}}$ of the electron jump to the on-site
interaction constant $U+J$ should be not very small, because only in that case
the TOHM is significantly different from the Heisenberg model.

## References

  * [1] Chao K.A., Spałek J., Oleś A.M., J. Phys. C: Solid State Phys., 1977, 10, L271,
doi:10.1088/0022-3719/10/10/002.

  * [2] Castellani C., Natoli C.R., Ranninger J., Phys. Rev. B, 1978, 18, 4945, doi:10.1103/physrevb.18.4945.
  * [3] Pizarro J.M., Calderón M.J., Liu J., Muñoz M.C., Bascones E., Phys. Rev. B, 2017, 95, 075115, doi:10.1103/PhysRevB.95.075115.
  * [4] Affronte M., Guidi T., Caciuffo R., Carretta S., Amoretti G., Hinderer J., Sheikin I., Jansen A.G.M., Smith A.A., Winpenny R.E.P., van Slageren J., Gatteschi D., Phys. Rev. B, 2003, 68, 104403, doi:10.1103/PhysRevB.68.104403.
  * [5] Kozłowski P., Kamieniarz G., J. Nanosci. Nanotechnol., 2011, 11, 9175, doi:10.1166/jnn.2011.4301.
  * [6] Brzostowski B., Lemański R., Ślusarski T., Tomecka D., Kamieniarz G., J. Nanopart. Res., 2013, 15, 1528, doi:10.1007/s11051-013-1528-2.
  * [7] Chiesa A., Carretta S., Santini P., Amoretti G., Pavarini E., Phys. Rev. Lett., 2013, 110, 157204,
doi:10.1103/physrevlett.110.157204.

  * [8] Plakida N.M., Yushankhai V.Yu., Stasyuk I.V., Physica C, 1989, 160, 80, doi:10.1016/0921-4534(89)90456-5.
  * [9] MacDonald A.H., Girvin S.M., Yoshioka D., Phys. Rev. B, 1988, 37, 9753, doi:10.1103/PhysRevB.37.9753.
  * [10] Cleveland C.L., Medina R., Am. J. Phys., 1976, 44, 44, doi:10.1119/1.10537.
  * [11] Chao K.A., Spałek J., Oleś A.M., Phys. Lett. A, 1977, 64, 163, doi:10.1016/0375-9601(77)90702-2.
  * [12] Chao K.A., Spałek J., Oleś A.M., Phys. Rev. B, 1978, 18, 3453, doi:10.1103/physrevb.18.3453.
  * [13] Micnas R., Ranninger J., Robaszkiewicz S., Rev. Mod. Phys., 1990, 62, 113, doi:10.1103/revmodphys.62.113.
  * [14] Manousakis E., Rev. Mod. Phys., 1991, 63, 1, doi:10.1103/revmodphys.63.1.
  * [15] Frésard R., Kotliar G., Phys. Rev. B, 1997, 56, 12909, doi:10.1103/physrevb.56.12909.

Ukrainian 0

Ïîðâíÿííÿ äâîîðáòàëüíî¿ ìîäåë Õàááàðäà òà ñïí $S=1$ ìîäåë Ãàéçåíáåðãà:
äîñëäæåííÿ íà êëàñòåðàõ Ð. Ëìàíüñê, ß. Ìàòèñÿê

íñòèòóò íèçüêèõ òåìïåðàòóð òà ñòðóêòóðíèõ äîñëäæåíü Ïîëüñüêî¿ àêàäåì¿ íàóê,
âóë. Îêóëüíà 2, 50-422 Âðîöëàâ, Ïîëüùà


# Torque-based Deep Reinforcement Learning for Task-and-Robot Agnostic
Learning on Bipedal Robots Using Sim-to-Real Transfer

Donghyeon Kim1, Glen Berseth2,∗, Mathew Schwartz3 and Jaeheung Park1,4,5,6,∗
*This work was supported by the National Research Foundation of Korea (NRF)
grant funded by Korea government (MSIT) (No.2021R1A2C3005914)1Donghyeon Kim
and Jaeheung Park are with the Department of Intelligence and Information,
Seoul National University, Seoul, KS013, Republic of Korea kdh0429,
<EMAIL_ADDRESS>Berseth is with the Université de Montréal, Quebec, H3C
3J7, Canada<EMAIL_ADDRESS>Schwartz is with the College of
Architecture and Design, New Jersey Institute of Technology, Newark, 07102,
USA cadop@njit.edu4,5,6Jaeheung Park is also with ASRI, RICS, Seoul National
University, Republic of Korea, and Advanced Institutes and Advanced Institutes
of Convergence Technology(AICT), Republic of Korea.∗Both are corresponding
authors.

###### Abstract

In this paper, we review the question of which action space is best suited for
controlling a real biped robot in combination with Sim2Real training. Position
control has been popular as it has been shown to be more sample efficient and
intuitive to combine with other planning algorithms. However, for position
control gain tuning is required to achieve the best possible policy
performance. We show that instead, using a torque-based action space enables
task-and-robot agnostic learning with less parameter tuning and mitigates the
sim-to-reality gap by taking advantage of torque control’s inherent
compliance. Also, we accelerate the torque-based-policy training process by
pre-training the policy to remain upright by compensating for gravity. The
paper showcases the first successful sim-to-real transfer of a torque-based
deep reinforcement learning policy on a real human-sized biped robot. The
video is available at https://youtu.be/CR6pTS39VRE.

###### Index Terms:

Reinforcement Learning, Humanoid and Bipedal Locomotion, Torque-based Control

##  I Introduction

Control algorithms for robots can be extremely complex and difficult to
implement. They change widely based on the robot themselves and tasks, and a
continuing line of robotics research has shown the difficulties faced when
high-dof multi-link robots, such as biped robots, are controlled [1, 2].
Recently, reinforcement learning (RL) has become a popular method for
implementing robotic controllers, in part due to its ability to be applied to
challenging tasks and robots. For example, by introducing RL control methods
to biped robots, challenging tasks could be performed while effectively
utilizing the advantages of biped robots.

Figure 1: Diagram of position-based deepRL and torque-based deepRL. For
position-based deepRL, the PD gain which is proper for task 1, robot M, and
sim-to-real transfer is tuned in advance. In torque-based deepRL, we propose
to pre-train the policy with the gravity compensation torque.

Biped robots (and other legged robots) make discontinuous contact with the
ground, enabling greater mobility compared to wheeled robots. However, the
discontinuous contact requires switching the contact discretely at proper
timing and location [2, 3], making the biped robot challenging to control.
Additionally, the control challenge varies depending on the task and robot,
and it has resulted in designing solutions that work for more specific robots
and tasks. For example, although a target pose control (position-based
deepRL), which is a widely-used control space in RL, has many advantages, such
as it is more explainable and it promotes fast learning [4], it comes at the
cost of introducing what is essentially a mass-spring-damper system with gain
parameters into the control design. When controlling biped robots, the gain
parameters are often tuned for each robot and every task the robot performs
[5, 6, 7, 8] as shown in Fig. 1. This tuning complicates the creation of
robot-and-task-agnostic algorithms and requires manually searching PD gains,
and adding the gains to the optimization space increases the size and
difficulty of the optimization [8].

Additionally, since the reality gap exists between simulation and the real-
world due to various factors, such as contact dynamics and state estimation
gap [7, 9], we should promote control methods that can be used more broadly
whether in a simulation or real robot. While transferring the trained policy
to the robot hardware, the reality gap would result in a tracking performance
difference and the timing of contact with the ground cannot exactly match
between the trained environment and the real-world. However, when early
contact occurs in the position-based deepRL policy, the PD controller will try
to track the target position no matter how much torque is generated, and the
robot can lose balance due to high impact force [10]. In this sense,
compliance is advantageous to overcome the reality gap because compliance
implicitly controls the energy transfer to the environment and reduces the
impact force [11].

Therefore, it would be better to develop control methods that require less
parameter tuning and reduce the reality gap, because they will scale better,
be more reusable, and be independent of the robot platform, task, and
environment (simulation or real). One option that fits nicely into this puzzle
is to control the motor torques directly as shown in Fig. 1 because there is
no PD gain that predetermines the behavior and the torque-based control is
known to be compliant [10, 11]. Early work on deep RL [12] has compared
position-based and torque-based deepRL and reported that the torque-based
deepRL generally shows inferior performance with slow training. However, the
performance was evaluated only on a simulator, and the Sim2Real problem was
not taken into account. Therefore, it would be worthy to analyze the action
space on the real robot from a robotics perspective.

In this paper, we revisit the question of which action-space is the best in
general. First, from the learning perspective, we show that the torque-based
deepRL is a general algorithm that can be applied out of the box without
manually tuning gains for each robot (Atlas, TOCABI) or task (squat, walk,
run) compared to position-based deepRL. Second, from the sim-to-real
perspective, we include an analysis of real robot hardware and show that the
position-based deepRL policy is more vulnerable to the reality gap due to its
limited compliance. We also show that the torque-based deepRL policy is
compliant by its nature, and effectively overcomes the reality gap. Third,
since torque-based policies have been described as exhibiting slow training,
we propose to pre-train the torque-based deepRL policy with a gravity
compensation torque to accelerate initial training and reduce the sample
inefficiency. Lastly, this is the first time that a torque-based deepRL policy
is successful on a heavy human-sized biped robot ($\approx$100kg) with stable
walking. It is demonstrated that a biped robot can be torque-controlled at
various control frequencies with RL, ranging from 250 Hz to as low as 62 Hz,
which could suggest a new method to the torque control society. We provide
evidence for these claims by comparing position-based and torque-based deepRL
policies to each other, and this results in a single control algorithm that
achieves a larger space of tasks without additional manual tuning and is
adequate for real-world legged robots.

##  II Related Work

###  II-A RL on Legged Robot

RL is being actively applied on legged robots because it can exploit full-
model dynamics, leading to higher performance with less domain knowledge on
complex floating-base robot dynamics with discrete contact [4, 13, 14, 15]
while showing robust behavior compared to model-based controllers [4, 16].
Many researchers have demonstrated RL on legged robots not only in simulations
but also on real robots [17, 18], although there are not enough hardware
verifications on humanoid robots yet [16, 19, 9, 20]. However, when moving
research interests from simulation to the real robot, early works [15, 4]
adopted popular action space used in the graphics [12, 6] and did not fully
consider how the action space would affect on the real robot, and the chosen
action space became a common standard.

###  II-B Action Space on Legged Robots

When controlling a legged robot with RL, various action spaces can be
designed. Not only action spaces that were compared in early work on RL by
[12], but also any other action space such as task space [21] can be adopted
by using an additional low-level controller to track the target output by the
policy. The most widely used approach is a method that involves using a policy
to generate a target joint angle $q_{target,\pi_{\theta}}$ and tracking it
with a PD controller using gain parameters $K_{p},K_{d}$ (position-based
deepRL policy). Plenty of works using the position-based deepRL emphasized the
importance of selecting a proper gain. However, how to determine the PD gain
in RL is not unified and decided by the rule of thumb. Since each task and
robot requires a different amount of torque range, the PD controller requires
some insight to set properly for each individual robot morphology and task [5,
6]. For example, [22] mentioned that they used a relatively very low PD gain
on a humanoid robot that has heavy legs and a high-ratio gearbox. . On the
other side, [9] emphasized tuning the PD controller to behave similarly on the
real robot and the simulation is critical because the gap of the PD controller
response can result in non-observed states and increase instabilities.

Alternatively, the action space can be designed to generate the torque command
directly from the policy. To the best of the authors’ knowledge, the only work
in this framework on a legged robot has been implemented by [23] on a
quadruped robot. They implemented the torque-based deepRL policy with a
relatively high control rate and showed successful sim-to-real outdoor
walking. However, their work was implemented on a quadruped robot whose impact
force could be low, and it has not been discussed which aspects caused the
advantages of torque control. In this paper, we demonstrate that the torque-
based deepRL policy can be applied on various robot platforms and tasks with
minimal tuning along with a comparison of behavioral differences arising from
intrinsic compliance.

###  II-C Compliance of Position Control and Torque Control

When a robot works in a known and structured environment, position control can
be adopted to track the target with a relatively simple setup. However, when
the task involves contact with an unstructured environment and the tracking
error unexpectedly increases, the PD controller tries to reduce the error with
all available torque [10]. Due to this stiff disturbance rejection behavior,
compliance is lacking in the position control, and a perturbation is created
and transmitted into the robot. Meanwhile, when using torque control, the
torque does not abruptly change, and the robot shows robust and safe
interaction with the environment without sacrificing tracking accuracy [11,
24, 25].

###  II-D Sim-to-Real

A fundamental solution to reduce the discrepancies between the real-world and
the simulation is to bring the simulator closer to reality by developing a
high-fidelity simulator [4], or identifying the robot model accurately [15].
However, since it is hard to completely imitate reality in the simulator,
dynamics randomization is applied to train the robot robustly in a certain
range of randomized parameters [26, 27]. In addition, implicitly identifying
the randomized environment parameter in the latent vector and using the
information as part of a state [14, 13] is also being adopted. In this paper,
we show a new way to reduce the reality gap by exploiting the intrinsic
compliance of the controller in the environment where contact is involved.

##  III Background

###  III-A Reinforcement Learning

We model a robot control problem as a discrete-time Markov Decision Process
(MDP) with a discounted expected return objective. At each time step $t$, the
robot control policy $\pi_{\theta}$ observes a state $s_{t}$, and the policy
samples an action $a_{t}\sim\pi_{\theta}(\cdot|s_{t})$ based on the state. The
agent applies the action, transitioning the robot state $s_{t}$ to a new state
$s_{t+1}\sim p(\cdot|s_{t},a_{t})$, and a reward
$r_{t}=r(s_{t},a_{t},s_{t+1})$ is calculated accordingly. The agent tries to
learn a policy through interaction with the environment that maximizes the
expected return $J(\pi_{\theta})$, which is the expected cumulative discounted
reward over a finite-horizon T

$J(\pi_{\theta})=E_{\tau\sim
p(\tau|\pi_{\theta})}[\sum_{t=0}^{T-1}\gamma^{t}r_{t}],$ (1)

where $\tau$ is the trajectory when executing the policy $\pi_{\theta}$.

###  III-B Biped Robot TOCABI

In this paper, we use a lab-made biped robot TOCABI for hardware verification.
TOCABI is a human-sized humanoid that is 1.8 meters tall and weighs 100
kilograms [28]. It is designed according to the proportions of the human body
and has 12 actuated joints in the lower body and 21 actuated joints in the
upper body. Gears with a high reduction ratio of 100:1 are used due to heavy
weight, and the current control-based servo drives, which communicate with the
real-time control computer through the EtherCAT interface, control each joint
motor with a sampling frequency of 2kHz.

##  IV Method

To find the best algorithm that can be used to train various tasks on various
robot platforms and can be transferred from simulation to the real-robot with
minimal tuning, we set two environments. In these environments, the agent
generates either target joint angles (position-based deepRL policy) or torque
commands (torque-based deepRL policy) following the MDP formulation in Sec.
III-A as follows. A diagram of each control method is also shown in Fig. 1,
and in this paper, we focus on the lower body, and the upper body is PD
controlled along with a gravity compensation to keep its default pose.

###  IV-A State Space

The state space $S\in R^{30}$ consists of the robot’s base orientation in
Euler angles $\alpha_{t}\in R,\beta_{t}\in R,\gamma_{t}\in R$, the joint
position $q_{t}\in R^{12}$, the low-pass filtered joint velocity
$\dot{q}_{lpf,t}\in R^{12}$, the phase information represented by sine and
cosine $sin(2\pi\phi_{t})\in R,cos(2\pi\phi_{t})\in R$, and the command
velocity $v_{t}^{cmd}\in R$.

$s_{t}=[\alpha_{t},\beta_{t},\gamma_{t},q_{t},\dot{q}_{lpf,t},sin(2\pi\phi_{t}),cos(2\pi\phi_{t}),v_{t}^{cmd}].$
(2)

The phase variable $\phi_{t}$ increments with time from 0 to 1 and then
returns to 0 when it reaches 1 to represent the current phase of the cyclic
motion. However, the phase $\phi_{t}$ can also be modulated with the phase
modulation action $a_{\delta\phi,t}$ as it will be described in Sec. IV-B3.
The command velocity $v_{t}^{cmd}$ is given to the robot when the task
involves tracking a user-input velocity (e.g. walk, run). To simulate the
sensor noise, we inject a Gaussian noise $w_{t}\sim N(0,\sigma)$ to the joint
position according to the encoder resolution of the robot hardware
($\sigma=1e^{-4}rad$). A joint velocity $\dot{q}_{noise,t}$ is computed from
the noise-injected joint position and then low-pass filtered with a cut-off
frequency of 4Hz to be used as a state $\dot{q}_{lpf}$.

###  IV-B Action Space

The action space $A\in R^{13}$ is composed of $12$ actions to generate joint
commands and one action to modulate the phase $\phi_{t}$. The action is
sampled from a Gaussian distribution whose mean is the output of the policy
($\mu_{\theta}$), and the standard deviation ($\Sigma$) is fixed to a
predetermined value $\pi_{\theta}\sim N(\mu_{\theta},\Sigma^{2})$. Depending
on the control method, either position-based or torque-based deep RL policy
generates commands as follows.

####  IV-B1 Position-based DeepRL

In the position-based deepRL method, the policy outputs the target joint
angle, and the target joint angle is tracked by a low-level PD controller.

$\tau_{cmd}=K_{p}(q_{\text{target},\pi_{\theta}}-q)+K_{d}(-\dot{q}).$ (3)

In general, the actor-network updates the target joint angle with 30-100Hz,
and the PD controller runs faster with 500-2000Hz.

In our implementation, the actor updates the action with 250Hz and the PD
controller runs at 2000Hz. The default pose $q_{default}$ is set to a knee-
bent standing configuration. Lastly, the standard deviation $\Sigma$ is set to
a scaled value of each joint position limit ($q_{limit}/s_{q}$) rather than
the same value across all joints. This will respect the hardware design by
allowing larger exploration at the joint whose joint limit is designed to be
high to move a wider range of motion.

####  IV-B2 Torque-based DeepRL

A torque-based deepRL policy directly outputs the torque,

$\tau_{cmd}=\pi_{\theta}.$ (4)

In contrast to the position-based deepRL policy, whose behavior is largely
predetermined by the PD gain, there is no inductive bias on the torque-based
deepRL. The bias introduced from the gains has made learning policies for
specific tasks easier but has also limited the reuse of control methods.
Therefore, torque-based deepRL, as we will show, can be used as a more robot-
and-task agnostic method for learning control policies. Also, because the
generated control input is not restricted to a mass-spring-damper-like system,
there is a higher possibility of achieving improved performance than the
position-based deepRL policy. Additionally, as opposed to the position-based
deepRL policy whose compliance is restricted by the PD controller, the torque
does not abruptly change. This compliance results in robust and safe
interaction with the environment, especially for unexpected contacts.

In our implementation, the actor updates the command torque with 250Hz. Also,
the standard deviation is set to a scaled value of each joint torque limit
($\Sigma=\tau_{limit}/s_{\tau}$), and this will allow larger torque
exploration to the joints whose torque capacity is large, respecting the
hardware design.

####  IV-B3 Phase Modulation

The action is also composed with a phase modulation action $a_{\delta\phi,t}$.
The phase variable $\phi_{t}$ represents the current phase in cyclic motion
and the reference motion is updated accordingly. However, if the phase
variable advances linearly with time by a fixed value, the period of the
learned motion is also determined accordingly and cannot be modulated [6].
This would generate sub-optimal behavior when the robot is commanded to walk
at high speed. The robot could generate a wide step with a fixed period
$T_{ref}$ to track the target speed, but it would be better if the robot could
also reduce the walking period. Therefore, in our implementation, a phase
modulation action $a_{\delta\phi,t}$ is additionally added to adjust the
period and timing of the motion as follows.

$\phi_{t+1}=fmod(\phi_{t}+{\Delta t\over T_{ref}}+a_{\delta\phi,t},1.0)$ (5)

The action space of phase modulation action $a_{\delta\phi,t}$ is designed to
be within $[0,\Delta t]$.

###  IV-C Reward

The reward is formulated to execute a given task while imitating a reference
motion. In addition, several regularization rewards are used for real-robot
implementation. The overall reward is composed as follows.

$\begin{split}r_{t}&=r^{base,imitate}_{t}+r^{q,imitate}_{t}+r^{c,imitate}_{t}+r^{v,task}_{t}\\\
&+r^{\dot{q},reg}_{t}+r^{\ddot{q},reg}_{t}+r^{F,reg}_{t}+r^{\Delta
F,reg}_{t}+r^{\tau,reg}_{t}+r^{\Delta\tau,reg}_{t}\end{split}$ (6)

The reward can be categorized into imitation reward, task reward, and
regularization reward, and each term is defined in Table I. Imitation rewards
guide the robot to follow the reference base orientation expressed in
quaternion $q^{ref}_{base,t}$, reference joint angle $q^{ref}_{t}$, and foot
contact status. Specifically, the foot contact imitation reward
$r^{c,imitate}_{t}$ is given if the current foot contact status
$c_{r,t},c_{l,t}\in\\{0,1\\}$ matches with a predetermined desired foot
contact status (double-support phase $\phi_{DSP}$ or single-support phase
$\phi_{SSP,r},\phi_{SSP,l}$). The task reward encourages the robot to track
the command velocity $v_{t}^{cmd}$ on locomotion tasks. Lastly, regularization
rewards penalize joint velocity, joint acceleration, contact force, contact
force difference, joint torque, and joint torque difference, respectively for
real-world implementation.

TABLE I: Reward Definition Reward | Definition
---|---
$r^{base,imitate}_{t}$ | $0.3\cdot exp(-13.2||q^{ref}_{base,t}-q_{base,t}||$
$r^{q,imitate}_{t}$ | $0.35\cdot exp(-4.0||q^{ref}_{t}-q_{t}||_{2}^{2})$
$r^{c,imitate}_{t}$ | $\begin{cases}0.2\text{ if}\begin{cases}c_{r,t}=1,c_{r,t}=1,\phi_{t}\in\Phi_{DSP}\\\ c_{r,t}=1,c_{r,t}=0,\phi_{t}\in\Phi_{SSP,r}\\\ c_{r,t}=0,c_{r,t}=1,\phi_{t}\in\Phi_{SSP,l}\\\ \end{cases}\\\ 0\text{ else }\end{cases}$
$r^{v,task}_{t}$ | $0.3\cdot exp(-3.0||v_{t}^{cmd}-v_{t}||_{2}^{2})$
$r^{\dot{q},reg}_{t}$ | $0.05\cdot exp(-0.01||\dot{q}_{t}||_{2}^{2})$
$r^{\ddot{q},reg}_{t}$ | $0.05\cdot exp(-20.0||\ddot{q}_{t}||_{2}^{2})$
$r^{F,reg}_{t}$ | $0.1\cdot exp(-0.0005(||F_{r,t}||_{2}+||F_{l,t}||_{2}))$
$r^{\Delta F,reg}_{t}$ | $0.1\cdot exp(-0.0005(||\Delta F_{r,t}||_{2}+||\Delta F_{l,t}||_{2}))$
$r^{\tau,reg}_{t}$ | $0.05\cdot exp(-0.01(||\tau_{cmd,t}||_{2}))$
$r^{\Delta\tau,reg}_{t}$ | $0.2\cdot exp(-0.01(||\tau_{cmd,t}-\tau_{cmd,t-1}||_{2}))$

###  IV-D Gravity Pre-training

Generally, since $\pi_{\theta}$ is initialized to a small value, the robot
must first learn to support its weight in the early stages of learning, which
deteriorates sample efficiency. To accelerate the initial training and
increase sample efficiency, we propose a relatively weak inductive bias on the
torque-based deepRL policy without damaging compliance. This method only
requires pre-training a torque-based policy to maintain the initial pose of
the robot (gravity compensation).

To collect the pre-training data, the gravity compensation torque-based policy
is obtained based on the contact-consistent whole-body model [29]. The gravity
torque is computed based on the base orientation
$\alpha_{t},\beta_{t},\gamma_{t}$, current joint configuration $q_{t}$, and
contact state of each leg $c_{r,t},c_{l,t}$, and experience is collected
across various configurations by randomly initializing the robot and
perturbing the robot. Except for state variables that affect the gravity
torque ($\alpha_{t},\beta_{t},\gamma_{t},q_{t},c_{r,t},c_{l,t}$), other state
variables are randomly set at each time step. Likewise, the phase modulation
action $a_{\delta\phi,t}$ is also set to a random value. $200,000$ experience
samples are collected on 8 parallel environments for 5 hours and required less
than 10 minutes to pre-train the policy. Note that once the policy is pre-
trained, it can be reused as long as the state space and action space are not
modified.

###  IV-E Training Setup

#####  IV-E0a Dynamics Randomization

To narrow the reality-gap and improve the robustness in the real-world
implementation, the dynamics parameters are randomized during training, as
shown in Table II. The default values of
$m_{default},I_{default},p_{CoM,default}$ are acquired from the robot CAD
model, while the default values of joint friction $v_{joint,default}$ and
joint damping $f_{joint,default}$ follow the gear specification sheet. Since
our real robot actually controls not a torque but a motor current, as opposed
to the simulation, it is also important to randomize the motor constant.
Additionally, the delay is randomized between 2ms and 6ms, considering the
delay of mechanical response.

TABLE II: Dynamics Randomization Parameters Parameter | Range | Unit
---|---|---
Link Mass | ${[}0.6,1.4{]}\times m_{default}$ | kg
Link Inertia | ${[}0.6,1.4{]}\times I_{default}$ | kg$\cdot m^{2}$
Link Center of Mass | ${[}0.6,1.4{]}\times p_{CoM,default}$ | m
Joint Damping | ${[}0.6,1.4{]}\times v_{joint,default}$ | Nm$\cdot$s/rad
Joint Friction | ${[}0.6,1.4{]}\times f_{joint,default}$ | Nm
Motor Constant | ${[}1.0,1.1{]}\times c_{motor,default}$ | -
Delay | ${[}0.5,1.5{]}\times t_{delay,default}$ | ms

#####  IV-E0b RL Algorithm

As an RL algorithm, Proximal Policy Optimization (PPO) is used to train the
policy. Both the actor and the critic are modeled as multi-layer perceptions
with 2 hidden layers of 256 ReLU units. The agent is updated with a batch size
of $128$ after $16,384$ samples are collected, and $80,000,000$ samples are
collected in total. The learning rate linearly decreases from $5e^{-5}$ to
$1e^{-6}$ as the training progresses. The maximum episode length is set to
16s, which is $8,000$ simulation steps, and the early termination occurs when
the robot’s links, other than both feet, contact the ground.

#####  IV-E0c Tasks

In this paper, we selected three tasks (squat, walk, run) so that the
difficulty level of each task differs. In the squat motion, the robot bends
its knee for 2s, maintains its pose for 2s, and then stretches the knee for 2s
without requiring foot contact switching, and the reference motion for the
squat is handcrafted. For the walking task, the robot is commanded to follow
the target velocity by switching the contact, and a single clip of walking
motion with a $1.8$s cycle is acquired from the model-based controller [30]
for reference motion. For the running task, the robot follows a wider range of
target velocity, and the motion capture data from DeepMimic [6] is used as a
reference motion.

Figure 2: Trained squat motion Figure 3: Trained running motion

##  V Result

We aim to answer four questions about torque-based policies: (1) Does pre-
training with gravity compensation torque accelerate the learning of the
torque-based policy (Sec. V-A), (2) Can the torque-based policy be applied to
a wider range of tasks and robots without additional tuning (Sec. V-B), (3)
Does the torque-based policy exhibit compliant behavior (Sec. V-C1), and (4)
Does the compliant behavior of the torque-based policy reduce the reality gap
in sim-to-real transfer (Sec. V-C2).

###  V-A Effect of gravity pre-training

When using the torque-based deepRL policy, the sample efficiency in the early
stages of learning is normally low. First, we verified that the pre-trained
policy is capable of maintaining its pose on various configurations,
indicating that the gravity torque is properly learned and compensated. We
then train the torque-based deep RL policy on squat and walking tasks, using
the pre-trained model as a starting point, to ensure that pre-training works
regardless of whether the task requires foot contact transitions. The learning
curves of the position-based and the torque-based deepRL policy with/without
gravity pre-training are shown in Fig. 4 with mean and variance of three runs
for each task and policy. The pre-trained torque-based policy learns at a
comparable or faster speed than the position-based deepRL policy depending on
the task because the pre-trained torque policy can stand still from the
beginning of training, similar to the position-based deep RL policy. Also, as
training progresses, the pre-trained torque policy converges faster than the
torque policy trained from scratch, demonstrating that gravity pre-training
accelerates the training procedure.

Figure 4: Effect of gravity pre-training

###  V-B Torque-based RL for Task-and-Robot Agnostic Learning

For the position-based deepRL policy, PD gains should be tuned depending on
the task and robot to achieve optimal policies, unlike the torque-based deepRL
policy. In this subsection, we aim to verify this by showing that tuning gains
is critical for the training(Experiment 1), and tuning gains that work across
tasks (Experiment 2) and robot platforms (Experiment 3) is difficult for the
position-based deepRL policy compared to the torque-based deepRL policy in
simulation.

Experiment 1: The behavior of the robot in position-based deepRL is
significantly impacted by the PD gain, which in turn affects the training
process. In this experiment, we investigate the effect of the gain and
standard deviation of the action distribution on the training of a single task
(walk) on a single robot (TOCABI). We examine 16 different combinations by
testing four different gains and four standard deviations. The gain starts
from $K_{p}=K_{default}$, which is the same gain used by the model-based
controller[30], and then it is halved $K_{p}=K_{default}/s_{p}$, where
$s_{p}=\\{1,2,4,8\\}$. The standard deviation of the action distribution is
set to a scaled value of the joint limit $\Sigma=q_{limit}/s_{q}$, where
$s_{q}=\\{100,200,400,800\\}$ respectively. Note that $K_{p}=K_{default}$ with
a standard deviation of $\Sigma=q_{limit}$/100 would generate the largest
exploration torque during training.

In Fig. 5, we display the learning curves for each gain and action standard
deviation. It is observed that the robot cannot learn to walk when a high gain
of $K_{p}=K_{default}$ is used. The gain is a critical parameter, and finding
an appropriate value that balances compliance and tracking performance for the
given task is necessary. However, upon examining four runs when
$K_{p}=K_{default}$/2, we discover that not only the gain but also the
standard deviation of the action distribution affects the success of the
training. The robot begins to walk when the standard deviation is less than
$q_{limit}$/200. This demonstrates that setting a suitable gain for the
position-based deepRL policy is crucial, and that tuning the gain becomes even
more complex when combined with the standard deviation of the action
distribution.

Figure 5: Effect of gain and action noise on training

Experiment 2: To apply a control method on diverse tasks, the controller
should have the potential to satisfy the required specifications (e.g.
tracking performance, compliance) of multiple tasks. To understand the
limitations of tuning position-controlled policies compared to a torque-based
learning method, our second experiment compares position-based deepRL policies
with PD gains to a torque-based deepRL policy across multiple tasks on a
single robot. We attempted to train position-based deepRL policies across a
set of tasks (squat, walk, run) using a _single_ set of PD gain across all
tasks. However, it was difficult to find a single set of gains that could be
applied to all tasks simultaneously because each task required different
levels of tracking performance and compliance. To train the position-based
deepRL policy successfully, not only the PD gain but also the action standard
deviation should be tuned for each task. Specifically, by examining various PD
gains $K_{p}=K_{default}/s_{p}$, $s_{p}=\\{1,2,4,8\\}$, and action standard
deviations $\Sigma=q_{limit}/s_{q}$, where $s_{q}=\\{100,200,400,800\\}$, the
squat motion and walking motion could be learned by fixing the gain to
$K_{p}=K_{default}/4,K_{default}/8$ and tuning the action standard deviation
for each task ($s_{q}=100$ for squat, and $s_{q}=100,200,400$ for walking).
However, the running motion could not be learned at all with gains of
$K_{p}=K_{default}/4,K_{default}/8$, and $K_{p}=K_{default}/2$ should be
adopted. In contrast, the torque-based control method did not require further
parameter tuning across tasks, and it was relatively robust to the action
standard deviation by learning all tasks with $s_{\tau}=5,10,20$. This is
because the torque-based method can utilize a wider range of control space and
has a higher potential to satisfy the requirements of multiple tasks since it
is not pre-restricted and not inductively biased, unlike the position-based
deepRL policy that uses the PD controller.

Experiment 3: In the third experiment, we investigate whether the same gains
can be used as the robot platform changes and adopt Atlas and TOCABI as the
comparison candidates. These robots are chosen to demonstrate that, although
they are both human-sized robots with a mass of 90 kg and 104 kg,
respectively, and have the same order of lower-body joint direction (YRPPPR),
the gain is restricted to a specific robot. Specifically, when using the gains
of TOCABI $K_{P}=K_{default}/4,K_{default}/8$ in the position-based deepRL
policy, which is verified in Experiment 2, on Atlas for the walking task, the
simulation became unstable and the training could not progress. This is due to
the kinematics and dynamics differences between the two robots, making it
extremely difficult to find a single set of gains that would fit both robots.
Instead, we had to tune a separate gain for Atlas. In addition, the gain had
to be tuned for each joint by allocating a larger gain for joints that would
sustain gravity. In contrast, the same parameter (action standard deviation)
used in TOCABI $s_{\tau}=5,10,20$ could be simply transferred to Atlas, and no
trial-and-error was required in the torque-based deepRL policy. Moreover,
setting the action standard deviation $\Sigma$ to a scaled value of the torque
limit of each joint as proposed in Sec. IV-B2 was successful without
exception.

###  V-C Torque-based RL for Sim2Real Transfer

In this subsection, we evaluate the analysis of simulation and real robot
hardware to determine the generality of torque-based control and how its
built-in compliance can assist in Sim2Real transfer.

####  V-C1 Compliance Analysis

We start by demonstrating how the intrinsic compliance of position-based and
torque-based deepRL policies differ by examining the reaction of both policies
to unexpected contact (Experiment 4) and how the PD gain affects compliance in
the position-based deepRL policy (Experiment 5) on the TOCABI robot while
performing a walking task.

Experiment 4: In this experiment, we demonstrate how the position-based deepRL
policy and torque-based deepRL policy respond to unexpected contact and
analyze compliance behavior. To simulate the unexpected contact, the robot is
trained on flat terrain, and after training is completed, an obstacle is
placed in the environment in the simulation. The obstacle is modelled as a box
with a one cm height lying on the ground 60 cm from the robot’s starting
position and results in early contact with the ground. In Fig. 6, the command
torques of the ankle pitch joint of a foot stepping on the obstacle, which is
the most influenced joint by perturbation, are plotted. In Fig. 6(a), when the
robot steps on an unexpected obstacle at around 4s with position-based deepRL
policy, the command torque abruptly increases compared to when the robot walks
on flat terrain. This is because the low-level PD controller generates large
torque when the tracking error increases due to early contact. The sudden
impact is propagated from the foot to other links and the robot loses balance
eventually. However, as seen in Fig. 6(b), the torque-based deepRL policy does
not show any peak torque although the robot steps on an unexpected obstacle
and smoothly passes the obstacle.

We conducted the same experiment on the real robot and the time-lapse of both
policies encountering an obstacle is shown in Fig. 7. By examining the contact
force of both policies when stepping on the obstacle, the result showed that
the position-based deepRL policy produced a larger impact force when stepping
on an unexpected obstacle and fails to overcome the obstacle, eventually
losing its balance as in the simulation. Conversely, the torque-based deepRL
policy did not show any sudden torque-change when the robot encounters an
unexpected obstacle. This compliant behavior is consistent with the simulation
result and demonstrates that the torque-based deepRL policy shows greater
compliance than the position-based deepRL policy.

Figure 6: Command torque with/without obstacle in simulation Figure 7: Time-
lapse of position-based deepRL policy and torque-based deepRL policy stepping
on an unexpected obstacle

Experiment 5: In the position-based deepRL policy, the low-level PD controller
restricts compliance, and the level of compliance is predetermined by the PD
gain. In this experiment, we examine how compliance in the position-based
deepRL policy varies as we adjust the gain in the same simulation environment
as Experiment 4. To achieve this, we evaluate the success rate as the target
velocity changes when implementing position-based deep RL policies trained
with various gains and the torque-based deep RL policy. Success is defined as
not falling for 16 seconds, which is the maximum episode duration of the
training environment. The target velocity increases by 0.02 m/s in each
episode, with 25 trials ranging from 0.02 m/s to 0.5 m/s, which fall within
the range of target velocities used during training.

The position-based deepRL policy with $K_{p}=K_{default}/8$ ($s_{p}=8$), which
was verified in Sec. V-B, resulted in the robot losing balance 10 times out of
25 trials. Then, we inspected various PD gains of the position-based deepRL
policy by continuously decreasing the P gain by half, starting from the
proportional gain $K_{p}=K_{default}/8$. Each policy is trained from scratch
with the corresponding gain, and the results in Table III indicate that the
success rate increases as more compliance is introduced by using a low gain
requiring trial and error to find a proper gain for the desired level of
compliance. In contrast, the torque-based deepRL policy never fell and
remained stable at all target velocities regardless of which action standard
deviation $s_{\tau}$ was used.

Gain Scale Factor $s_{p}$ | 8 | 16 | 32 | 64 | 128
---|---|---|---|---|---
Success Rate | 15/25 | 18/25 | 22/25 | 22/25 | 25/25
TABLE III: Success rate of position-based deepRL policy according to the gain

####  V-C2 Benefits of Compliance for Overcoming Reality Gap

Since we have verified the inherent compliance of the torque-based deepRL
policy in Sec. V-C1, we demonstrate its benefits when deploying the trained
policy to environments other than those used for training. To this end, we
show that the torque-based deepRL policy can quickly adapt to an unstructured
environment with fine-tuning(Experiment 6) and advantageous for Sim2Real
transfer by exploiting its intrinsic compliance (Experiment 7, 8).

Experiment 6: In order to demonstrate the effectiveness of compliant
controllers in handling and adapting to unstructured environments, both
position-based and torque-based deep RL policies are initially trained on flat
terrain and then fine-tuned on uneven terrain. The height of the uneven
terrain is randomly generated using the ’hfield’ asset of MuJoCo, as depicted
in Fig. 8. The experiments conducted in Sec. V-C1 suggest that the torque-
based deepRL policy is less susceptible to unexpected uneven ground due to its
inherent compliance, and hence has a greater likelihood of survival. The fine-
tuning curve in Fig. 9 supports this, with the torque-based deepRL policy
exhibiting a higher episode reward at the beginning of fine-tuning by
surviving longer than the position-based deepRL policy. Consequently, the
torque-based deepRL policy is able to explore longer on various terrains and
adapts more quickly, ultimately resulting in a higher final performance than
the position-based deep RL policy. This reveals that, while the position-based
deepRL policy can learn how to handle unexpected environments by experiencing
various terrains, as shown in several studies [7, 13, 31], the inherent
compliance of the torque-based deepRL policy leads to more effective and
stable locomotion in unstructured environments.

Figure 8: Uneven terrain generated in MuJoCo Figure 9: Mean episode reward
during fine-tuning on uneven terrain

Experiment 7: To examine how compliance contributes to robust sim-to-real
transfer, we first _simulate_ the reality gap by training both policies
without dynamics randomization and then randomizing the parameters in
simulation after training is completed. After training both policies with
default parameters in Table. II, the parameters are randomly scaled between
$[0.7,1.3]$ from their default values. Additionally, a parameter that scales
leg length is also randomized to effectively simulate early or late contact
with the ground, resulting in a total of $n_{parameter}=8$ randomized
parameters. At the start of each episode, the parameters are uniformly
randomized, and the sampled parameters are logged with the episode reward and
a success indicator at the end of the episode. The success indicator is set to
True when the robot walks for the maximum episode length (16s), and the
episode reward is a cumulative reward during an episode to measure motion
quality.

Both policies encountered $4,870$ sets of randomized parameters, and the
position-based deepRL policy could withstand $n_{alive}=3,477$ sets of
parameters without falling, while the torque-based deepRL policy exhibited
greater robustness by successfully handling $n_{alive}=3,604$ sets of
parameters. When examining the episode rewards of the successful runs, as
shown in Fig. 10, the torque-based deepRL policy exhibited a higher mean
reward with less variance than the position-based deepRL policy. This
indicates that the torque-based deepRL policy produces higher-quality motion
by more accurately tracking the reference motion and generating smaller
contact forces as a result of its compliance.

To analyze which parameters each controller is sensitive to, we applied
principal component analysis (PCA) to the failure data ($D_{n_{fail}\times
n_{parameter}}$) for each policy. With PCA, we were able to identify an axis
that effectively represents the failure sets and determine the sensitive
parameters by identifying the dominant elements of the principal axis. The
results of the analysis revealed that both the position-based and torque-based
deepRL policies are primarily affected by delay, with dominant values of 0.708
and 0.625, respectively. Additionally, the principal axis showed an
interesting result indicating that the position-based deepRL policy is more
vulnerable to leg length randomization compared to the torque-based deepRL
policy. This finding supports the argument that the position-based deepRL
policy is susceptible to unexpected contact timing due to a lack of
compliance. However, the position-based deepRL policy demonstrated greater
resilience to motor constant and mass randomization.

Figure 10: Histogram of episode reward of alive runs

Experiment 8:

Figure 11: Contact force comparison in simulation and real robot

In this experiment, we demonstrate how compliance is helpful for Sim2Real
transfer. We transfer the trained policy to the real robot and commanded it to
walk on flat ground with a target velocity of 0.2m/s. Fig. 11 shows the
contact force of both policies on the real robot along with the contact force
in the simulation. In Fig. 11(a), the impact force of the position-based
deepRL policy on the real robot ($\approx$ 1800N) is much larger than in
simulation (<1400N), while in Fig. 11(b), the impact force of the torque
policy is very similar both on the real robot and in simulation due to the
compliance of the torque policy. Although the robot could walk on flat terrain
with both policies, the reality gap on the position-based deepRL policy
continuously resulted in a stamping motion of the foot due to high impact
force, and the motion was less stable. Also, the contact timing gradually
mismatched with the simulation in the position-based policy as the walking
progressed. Additionally, as demonstrated in Fig. 7, when a larger reality gap
is introduced by an unexpected obstacle, the position-based deepRL policy
could not withstand the reality gap. These comparisons indicate that torque-
based policies can more smoothly handle differences between the simulation and
the real world.

###  V-D Torque-based DeepRL Policy Control Frequency

Lastly, we argue that the torque-based deepRL policy can be applied with
various frequencies. In traditional humanoid torque control, it is common to
use high control frequencies of up to 1-4 kHz. However, the torque-based
deepRL policy was implemented on the real robot with control frequencies of
62.5 Hz, 125 Hz, and 250 Hz, and it was not difficult to reduce the control
rate without episode reward drop. The main difference was that the training
time was longer at lower control rates to collect the same amount of training
samples.

##  VI Conclusion

In this work, we investigated which action space is suitable not only for
task-and-robot agnostic learning but also for reducing the reality gap on
biped robots. By analyzing the proposed torque-based deepRL policy alongside
the widely-used position-based deepRL policy, it is demonstrated that the
torque-based deepRL policy can learn to squat, walk, and run with minimal
tuning. Additionally, the torque-based deepRL policy did not require further
parameter tuning when the robot platform changes from TOCABI to Atlas, making
it suitable for task-and-robot agnostic learning. Furthermore, it is shown
that the torque-based deepRL policy is inherently compliant, and this
compliance is beneficial when the trained policy encounters an environment
different from the one it was trained on by reducing the unexpected impact
force. Lastly, we have accelerated the training of the torque-based deepRL
policy by pre-training it with a gravity torque. This is the first successful
attempt to implement a torque control method with deep RL on a human-sized
humanoid, and we believe this could suggest a new way to actively take
advantage of torque-based control methods in RL.

## References

  * [1] J. Englsberger, G. Mesesan, A. Werner, and C. Ott, “Torque-based dynamic walking-a long way from simulation to experiment,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2018, pp. 440–447.
  * [2] A. W. Winkler, C. D. Bellicoso, M. Hutter, and J. Buchli, “Gait and trajectory optimization for legged systems through phase-based end-effector parameterization,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 3, pp. 1560–1567, 2018.
  * [3] R. Deits and R. Tedrake, “Footstep planning on uneven terrain with mixed-integer convex optimization,” in _2014 IEEE-RAS international conference on humanoid robots_. IEEE, 2014, pp. 279–286.
  * [4] J. Hwangbo, J. Lee, A. Dosovitskiy, D. Bellicoso, V. Tsounis, V. Koltun, and M. Hutter, “Learning agile and dynamic motor skills for legged robots,” _Science Robotics_ , vol. 4, no. 26, p. eaau5872, 2019.
  * [5] G. Feng, H. Zhang, Z. Li, X. B. Peng, B. Basireddy, L. Yue, Z. Song, L. Yang, Y. Liu, K. Sreenath, _et al._ , “Genloco: Generalized locomotion controllers for quadrupedal robots,” _arXiv preprint arXiv:2209.05309_ , 2022\.
  * [6] X. B. Peng, P. Abbeel, S. Levine, and M. Van de Panne, “Deepmimic: Example-guided deep reinforcement learning of physics-based character skills,” _ACM Transactions On Graphics (TOG)_ , vol. 37, no. 4, pp. 1–14, 2018.
  * [7] S. Gangapurwala, M. Geisert, R. Orsolino, M. Fallon, and I. Havoutis, “Rloc: Terrain-aware legged locomotion using reinforcement learning and optimal control,” _IEEE Transactions on Robotics_ , vol. 38, no. 5, pp. 2908–2927, 2022.
  * [8] X. B. Peng, G. Berseth, and M. Van de Panne, “Dynamic terrain traversal skills using reinforcement learning,” _ACM Transactions on Graphics (TOG)_ , vol. 34, no. 4, pp. 1–11, 2015.
  * [9] D. Rodriguez and S. Behnke, “Deepwalk: Omnidirectional bipedal gait by deep reinforcement learning,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 3033–3039.
  * [10] J. Buchli, M. Kalakrishnan, M. Mistry, P. Pastor, and S. Schaal, “Compliant quadruped locomotion over rough terrain,” in _2009 IEEE/RSJ international conference on Intelligent robots and systems_. IEEE, 2009, pp. 814–820.
  * [11] A. Calanca, R. Muradore, and P. Fiorini, “A review of algorithms for compliant control of stiff and fixed-compliance robots,” _IEEE/ASME transactions on mechatronics_ , vol. 21, no. 2, pp. 613–624, 2015.
  * [12] X. B. Peng and M. van de Panne, “Learning locomotion skills using deeprl: Does the choice of action space matter?” in _Proceedings of the ACM SIGGRAPH/Eurographics Symposium on Computer Animation_ , 2017, pp. 1–13.
  * [13] J. Lee, J. Hwangbo, L. Wellhausen, V. Koltun, and M. Hutter, “Learning quadrupedal locomotion over challenging terrain,” _Science robotics_ , vol. 5, no. 47, p. eabc5986, 2020.
  * [14] A. Kumar, Z. Fu, D. Pathak, and J. Malik, “Rma: Rapid motor adaptation for legged robots,” _arXiv preprint arXiv:2107.04034_ , 2021.
  * [15] J. Tan, T. Zhang, E. Coumans, A. Iscen, Y. Bai, D. Hafner, S. Bohez, and V. Vanhoucke, “Sim-to-real: Learning agile locomotion for quadruped robots,” _arXiv preprint arXiv:1804.10332_ , 2018.
  * [16] Z. Li, X. Cheng, X. B. Peng, P. Abbeel, S. Levine, G. Berseth, and K. Sreenath, “Reinforcement learning for robust parameterized locomotion control of bipedal robots,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 2811–2817.
  * [17] L. Smith, I. Kostrikov, and S. Levine, “A walk in the park: Learning to walk in 20 minutes with model-free reinforcement learning,” _arXiv preprint arXiv:2208.07860_ , 2022.
  * [18] T. Miki, J. Lee, J. Hwangbo, L. Wellhausen, V. Koltun, and M. Hutter, “Learning robust perceptive locomotion for quadrupedal robots in the wild,” _Science Robotics_ , vol. 7, no. 62, p. eabk2822, 2022.
  * [19] S. Masuda and K. Takahashi, “Sim-to-real learning of robust compliant bipedal locomotion on torque sensor-less gear-driven humanoid,” _arXiv preprint arXiv:2204.03897_ , 2022.
  * [20] Z. Xie, G. Berseth, P. Clary, J. Hurst, and M. van de Panne, “Feedback control for cassie with deep reinforcement learning,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2018, pp. 1241–1246.
  * [21] H. Duan, J. Dao, K. Green, T. Apgar, A. Fern, and J. Hurst, “Learning task space actions for bipedal locomotion,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 1276–1282.
  * [22] R. P. Singh, M. Benallegue, M. Morisawa, R. Cisneros, and F. Kanehiro, “Learning bipedal walking on planned footsteps for humanoid robots,” in _2022 IEEE-RAS 21st International Conference on Humanoid Robots (Humanoids)_. IEEE, 2022, pp. 686–693.
  * [23] S. Chen, B. Zhang, M. W. Mueller, A. Rai, and K. Sreenath, “Learning torque control for quadrupedal locomotion,” _arXiv preprint arXiv:2203.05194_ , 2022\.
  * [24] J. Englsberger, A. Werner, C. Ott, B. Henze, M. A. Roa, G. Garofalo, R. Burger, A. Beyer, O. Eiberger, K. Schmid, _et al._ , “Overview of the torque-controlled humanoid robot toro,” in _2014 IEEE-RAS International Conference on Humanoid Robots_. IEEE, 2014, pp. 916–923.
  * [25] M. Focchi, T. Boaventura, C. Semini, M. Frigerio, J. Buchli, and D. G. Caldwell, “Torque-control based compliant actuation of a quadruped robot,” in _2012 12th IEEE international workshop on advanced motion control (AMC)_. IEEE, 2012, pp. 1–6.
  * [26] 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)_. IEEE, 2018, pp. 3803–3810.
  * [27] Z. Xie, P. Clary, J. Dao, P. Morais, J. Hurst, and M. Panne, “Learning locomotion skills for cassie: Iterative design and sim-to-real,” in _Conference on Robot Learning_. PMLR, 2020, pp. 317–329.
  * [28] M. Schwartz, J. Sim, J. Ahn, S. Hwang, Y. Lee, and J. Park, “Design of the humanoid robot tocabi,” in _2022 IEEE-RAS 21st International Conference on Humanoid Robots (Humanoids)_. IEEE, 2022, pp. 322–329.
  * [29] J. Park and O. Khatib, “Contact consistent control framework for humanoid robots,” in _Proceedings 2006 IEEE International Conference on Robotics and Automation, 2006. ICRA 2006._ IEEE, 2006, pp. 1963–1969.
  * [30] M.-J. Kim, D. Lim, G. Park, and J. Park, “Humanoid balance control using centroidal angular momentum based on hierarchical quadratic programming,” in _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2022, pp. 6753–6760.
  * [31] X. B. Peng, G. Berseth, and M. Van de Panne, “Terrain-adaptive locomotion skills using deep reinforcement learning,” _ACM Transactions on Graphics (TOG)_ , vol. 35, no. 4, pp. 1–12, 2016.