# Production of dual species Bose-Einstein condensates of 39K and 87Rb

Cheng-Dong Mi  State Key Laboratory of Quantum Optics and Quantum Optics
Devices,
Institute of Opto-electronics, Shanxi University, Taiyuan, Shanxi 030006,
People’s Republic of China  Collaborative Innovation Center of Extreme Optics,
Shanxi University, Taiyuan, Shanxi 030006, People’s Republic of China  Khan
Sadiq Nawaz  State Key Laboratory of Quantum Optics and Quantum Optics
Devices,
Institute of Opto-electronics, Shanxi University, Taiyuan, Shanxi 030006,
People’s Republic of China  Collaborative Innovation Center of Extreme Optics,
Shanxi University, Taiyuan, Shanxi 030006, People’s Republic of China  Peng-
Jun Wang  pengjun<EMAIL_ADDRESS>State Key Laboratory of Quantum Optics and
Quantum Optics Devices,
Institute of Opto-electronics, Shanxi University, Taiyuan, Shanxi 030006,
People’s Republic of China  Collaborative Innovation Center of Extreme Optics,
Shanxi University, Taiyuan, Shanxi 030006, People’s Republic of China  Liang-
Chao Chen  State Key Laboratory of Quantum Optics and Quantum Optics Devices,
Institute of Opto-electronics, Shanxi University, Taiyuan, Shanxi 030006,
People’s Republic of China  Collaborative Innovation Center of Extreme Optics,
Shanxi University, Taiyuan, Shanxi 030006, People’s Republic of China  Zeng-
ming Meng  State Key Laboratory of Quantum Optics and Quantum Optics Devices,
Institute of Opto-electronics, Shanxi University, Taiyuan, Shanxi 030006,
People’s Republic of China  Collaborative Innovation Center of Extreme Optics,
Shanxi University, Taiyuan, Shanxi 030006, People’s Republic of China
Lianghui Huang  State Key Laboratory of Quantum Optics and Quantum Optics
Devices,
Institute of Opto-electronics, Shanxi University, Taiyuan, Shanxi 030006,
People’s Republic of China  Collaborative Innovation Center of Extreme Optics,
Shanxi University, Taiyuan, Shanxi 030006, People’s Republic of China  Jing
Zhang<EMAIL_ADDRESS>
<EMAIL_ADDRESS>State Key Laboratory of Quantum Optics and Quantum Optics
Devices,
Institute of Opto-electronics, Shanxi University, Taiyuan, Shanxi 030006,
People’s Republic of China  Collaborative Innovation Center of Extreme Optics,
Shanxi University, Taiyuan, Shanxi 030006, People’s Republic of China

###### Abstract

We report the production of 39K and 87Rb Bose-Einstein condensates (BECs) in
the lowest hyperfine states $|F=1,m_{F}=1\rangle$ simultaneously. We collect
atoms in bright/dark magneto-optical traps (MOTs) of 39K/87Rb to overcome the
light-assisted losses of 39K atoms. Gray molasses cooling on the D1 line of
the 39K is used to effectively increase the phase density, which improves the
loading efficiency of 39K into the quadrupole magnetic trap. Simultaneously,
the normal molasses are employed for 87Rb. After the microwave evaporation
cooling on 87Rb in the optically plugged magnetic trap, the atoms mixture is
transferred to a crossed optical dipole trap, where the collisional properties
of the two species in different combinations of the hyperfine states are
studied. The dual species BECs of 39K and 87Rb are obtained by further
evaporative cooling in optical dipole trap at a magnetic field of 372.6 G with
the background repulsive interspecies scattering length $a_{KRb}$ = 34 $a_{0}$
($a_{0}$ is the Bohr radius) and the intraspecies scattering length $a_{K}$ =
20.05 $a_{0}$.

##  I Introduction

Since Bose-Einstein condensates (BECs) were firstly observed in alkali atoms
bec1 ; bec2 , the ultracold atomic gases have attracted much attention to
study many interesting quantum phenomena, e.g. the BEC-BCS crossover, spin
orbit coupling crossover1 ; soc1 ; soc2 ; soc3 , Mott insulator in the optical
lattice, wangdw ; chen ; wang , and ultracold chemistry chemis1 ; sadiq2 .
Dual species ultracold gasses greatly expand the research area of this
platform to study even more complex physical phenomena that are not accessible
in a single component quantum gas, such as producing heteronuclear molecules
with long dipole-dipole interaction molecule1 ; molecule2 ; molecule3 ,
creating polarons near quantum criticality polaron1 ; polaron2 ; polaron3 and
observation of collective dynamics of a mixture of Bose and Fermi superfluids
dualsf1 . To date, several dual Bose-Bose species have been cooled to BECs in
experiments, including 23Na-87Rb fwang ,39K-87Rb jarlt , 39K-23Na schulze2018
, 41K-87Rb modugno ; Thalhammer , 88Sr-87Rb, 84Sr-87Rb pasquiou , Er-Dy
trautmann , and 133Cs-87Rb Lercher ; McCarron .

Recently, more interest is directed towards the 39K-87Rb dual-species BECs as
there is rich Feshbach resonance structure jarlt ; simoni for the precise
control of the interspecies scattering length at low magnetic field.
Nevertheless, there are many difficulties in achieving the BEC of 39K. First,
its negative background sacttering length chiara inhibit the achievement of
the BEC state. Second, the unresolved excited-state hyperfine structure of 39K
Prevedelli makes the standard sub-doppler laser cooling less efficient. Third,
the strong light-assisted losses of the dual species in the initial double
MOTs stage Marcassa result in a poor loading. All of these factors make the
achievement of 39K-87Rb dual-species BECs to be a challenging task.

In this paper, we report the experimental route to overcome these
difficulties, and finally obtain the BECs of 39K and 87Rb in the lowest
hyperfine states $|F=1,m_{F}=1\rangle$ simultaneously. Previously jarlt , the
39K-87Rb BECs were achieved with bright molasses on the D2 line in the
$|F=1,m_{F}=-1\rangle$ states at around 117.56 G of uniform magnetic Feshbach
resonance field. We follow a slightly different path to achieve the dual
species condensate in the $|F=1,m_{F}=1\rangle$ states using 372.6 G of
uniform magnetic Feshbach resonance field, and using the gray molasses on the
D1 line. At the three dimensional MOT phases, a dark spontaneous optical force
trap (dark SPOT) Ketterle is used to load the 87Rb atoms in $|F=1\rangle$
state in the central part without repumping beams, and a bright MOT loads the
39K atoms at same time. This scheme mitigates the atoms loss resulted from the
inelastic collisions and the light-induced loss of different species in the
MOT loading stage. To effectively cool the 39K, the gray molasses cooling
dipankar ; Salomon ; Sievers ; shi on D1 line is used after the loading stage
in magneto-optic trap. The obstacle of negative background sacttering length
is overcome by tuning the 39K scattering length to $a_{K}$ = 20.05 $a_{0}$ by
using an external homogeneous magnetic field at 372.6 G near a broad Feshbach
resonance centered at 403.4 G, where it is positive to facilitate sympathetic
cooling below 1$\mu$K with the background interspecies scattering length
$a_{KRb}$ = 34 $a_{0}$ during the final evaporation step in the optical dipole
trap. Using this scattering length improves the density of 39K due to the
sympathetic cooling compared to when cooling the 39K alone.

##  II The 39K and 87Rb double MOT

The experimental setup is shown in Fig. 1, which has been used for the
creation of 87Rb BEC in Ref. liangchao ; guangyu . The 87Rb and 39K atoms are
first simultaneously cooled and trapped in the 2D MOT chamber to form an
atomic beam. This beam is then pushed to the science chamber (at much lower
pressure of 10-9 Pa) with the help of a push beam containing 4 frequencies. In
the science chamber, the 39K atoms are collected in a 3D MOT while the 87Rb
atoms are collected in a dark SPOT MOT to reduce the loss of 39K atoms. We
also tried the double dark SPOT MOT to further reduce the light assisted
losses, however the 39K dark SPOT MOT did not work. This could be due to the
use of a single (unbalanced) hollow repumping beam in contrast to the use of
four hollow repumping beams (balancing each other in the opposite direction)
used in Stwalley1995 . Since the 87Rb does not need a high power repumping
beam (compared to the 39K which needs almost the same power as that of the
trapping/cooling beams due to its narrowly spaced excited state hyperfine
splitting), we do not need to balance the single hollow repumping beam for the
87Rb dark SPOT MOT. The different powers and frequencies used in the push beam
and the two MOT stages are listed in Table I while the Fig. 2 gives the
definitions of the detunings.

In the 3D MOT, there are six solid MOT laser beams (2.5 cm diameter). Each of
these six beams carries the Rb and K cooling and the repumping light only for
K, as shown in the insert image in Fig. 1. To produce the dark region in 87Rb
dark SPOT, an opaque disk with diameter of 10 mm is used to cast a shadow at
the center of the 87Rb repumping beam. The center of the hollow repumping beam
is filled with a depumping beam (shown in Fig. 2(b)). The repumping+depumping
beam is delivered to the 3D MOT separately from the six MOT beams and thus the
87Rb atoms are trapped at the trap center in $|F=1\rangle$ state out of the
cooling cycle. The depumping beam further depletes the atomic $|F=2\rangle$
state in the center of the MOT by optical pumping atoms to the $|F=1\rangle$
state. The dark SPOT MOT significantly improves the simultaneous loading of
both species with reduced light assisted losses Marcassa . The magnetic field
gradient in the 3D MOT stage is 4 G/cm along the $z$-axis which is provided by
a pair of coils in the anti-Helmholtz configuration as shown in Fig. 1 fanhao
. The same coils provide the homogeneous magnetic field for Feshbach resonance
when the current direction in one of the coils is changed.

The setup for MOT lasers is presented in Fig. 3. For 87Rb, the cooling and
repumping light for the D2 line are generated by separate lasers due to large
ground state hyperfine splitting of 6.8 GHz, and then frequency tuned by
double-passed AOMs. The cooling beam is then amplified by a tapered amplifier
(Rb-TA), as shown in Fig. 3(a). For 39K, a single master laser can generate
the cooling and repumping lights due to the smaller ground state hyperfine
splitting of 461 MHz. We use two master lasers, one for the MOT on D2 line and
the other one for the gray molasses on D1 line. The cooling beams from both
the master lasers are first combined together and then amplified by an
injection locked laser. A similar setup amplifies the 39K repumping light and
all of these are then amplified together by a single tapered amplifier (K-TA),
as shown in Fig. 3(b). This design for the 39K guarantees the perfect
collimation of the D2 and D1 line laser beams.

The compressed MOT (CMOT) step follows the MOT loading step by increasing the
magnetic field strength to 22 G/cm (duration 150 ms). The fine tuning of the
laser parameters during the CMOT reduces the radiation pressure which together
with the high magnetic field increases the density of the atomic sample. The
repumping+depumping beam for 87Rb is now blocked and the repumping light is
delivered to the atoms through the six MOT beams. The detunings and single
beam power of the lasers during the CMOT are shown in Fig. 2(c).

Figure 1:  The experimental setup for the dual 39K and 87Rb BECs showing the
two MOT cells and the ion pumps for each chamber. The Rb and K sources (at 40
and 45 0C respectively) supply atoms into the 2D MOT chamber. The science
chamber (enclosed by the quadrupole/Feshbach coils) collects the atoms in the
3D MOT from the push beams assisted atomic flux coming from the 2D MOT via a
differential pressure tube. The inset shows the top view of 3D MOT, the four
horizontal MOT beams and a hollow Rb repumping beam are visible.

After the CMOT, the magnetic field is completely turned off in 100 $\mu$s and
we perform normal molasses for the 87Rb but D1 line gray molasses for the 39K.
The D1 line gray molasses Salomon gives samples with higher density and atoms
number compared to the previously devised schemes using the D2 line laser
jarlt ; Landini . The one-photon detunings $\Delta_{GC}$ and $\Delta_{GR}$
both are blue detuned +3.33 $\Gamma$ forming a Raman $\Lambda$-type system in
which the dark and bright states formed by the coherent superposition of the
Zeeman sub-levels in the hyperfine levels $|F=1\rangle$ and $|F=2\rangle$ help
cool in the bright state and hold the atoms in the dark state thus avoiding
the chance of reheating by scattering of photons from the GM fields. The other
powers and detunings are given in Fig. 2(c). The molasses step lasts for 6.7
ms.

##  III Optical pumping and MW Evaporation

Figure 2:  Energy level diagram of 39K (a) and 87Rb (b), and the time sequence
of the detunings and powers of each laser. (a) The $\Delta_{KC}$ and
$\Delta_{KR}$ are the detunings of the trapping and repumping laser from the
atomic transitions shown in the D2 line (766.7 nm) of 39K. Similarly
$\Delta_{GC}$ and $\Delta_{GR}$ are the GM cooling and repumping laser
detunings from the transitions shown in the D1 line (770.1 nm). (b)
$\Delta_{RbC}$, $\Delta_{RbR}$ and $\Delta_{RbD}$ are the detunings of the
trapping, repumping and depumping laser respectively from the transitions
shown in the D2 line (780.2 nm) of 87Rb. (c) The time sequence of the powers
(red solid lines in the right vertical axis) and detunings (black dashed lines
in the left vertical axis) of the respective lasers is shown. The duration of
each step of the experiment is shown on the horizontal axis.

Table 1: The detunings and powers of the different laser beams during the MOT
loading. $\Gamma\simeq$6 MHz.  Figure 3:  The optical setup and frequency
shifting and locking schemes of the (a) 87Rb cooling and repumping lasers and
(b) the 39K D1 and D2 cooling and repumping lasers. The 87Rb repumping laser
does not need amplification while all of the remaining lasers are amplified by
tapered amplifiers (TAs) or injection locked lasers.

Right after the molasses step, the pumping stage follows where both the 39K
and 87Rb atoms are optically pumped to the $|F=2,m_{F}=2\rangle$ states. This
is done by applying a 2 G homogeneous magnetic field along the $x$-direction
and a $\sigma^{+}$-polarized pump light also along the same direction. The
time sequence for the optimized optical pumping of both the species is as
follows. The repumping light for 87Rb (detuning +6 $\Gamma$) remains on in the
six MOT beams while the 39K D2 line repumping light (detuning +1.6 $\Gamma$)
is turned on after the molasses step. Then after 0.5 ms the 39K D1 line pump
light (on resonance) is turned on for 0.7 ms and the 87Rb pump light (4
$\Gamma$) for 1.1 ms. The two repumping beams stay on for additional 0.1 ms.
This completes the optical pumping step.

To efficiently load the sample in the magnetic quadrupole trap, the magnetic
field gradient is ramped up to 26.5 G/cm in 1 ms and held on for 10 ms and
then to 62 G/cm in 200 ms and held at this field for 50 ms. Finally, the
magnetic trap is further compressed to 74 G/cm in 300 ms to increase the atoms
density for strong collisions in the following evaporation stage. The 20 W
green laser is turned on at full power focused ($\omega_{0}$=30 $\mu$m 1/e2
radius) at the center of the magnetic quadrupole trap to convert it to the
optically plugged magnetic quadrupole trap.

Now, the atoms of both species can be cooled by either RF evaporation or MW
evaporation of 87Rb only campbell . The RF evaporation transfers the hotter
$|F=2,m_{F}=2\rangle$ state atoms to other un-trapped Zeeman states of the
same hyperfine level but the MW evaporation transfers the
$|F=2,m_{F}=2\rangle$ state atoms to the lower hyperfine (un-trapped)
$|F=1,m_{F}=1\rangle$ state. The RF evaporation evaporates both the atomic
species but the MW only evaporates the 87Rb atoms with sympathetic cooling of
the 39K atoms. This MW evaporation technique compensates for the low number of
39K atoms in the sample at the cost of 87Rb atoms. During the MW evaporation
process, the MW frequency is scanned from 6894.7 to 6855.7 MHz in 6.2 s and
then to 6835.7 MHz in another 5 s. At the end of the evaporation stage, there
are around 1.17 $\times 10^{7}$ 39K atoms and 4.73 $\times 10^{7}$ 87Rb atoms
at about 45 $\mu$K. The magnetic field strength during the evaporation is held
constant at 74 G/cm.

##  IV Lifetime Measurement of Different State Mixtures

After the MW evaporation cooling, the atoms are transferred to the crossed
optical dipole trap, which is made by the intersection of two 1064 nm laser
beams that are frequency shifted (around 10 MHz) to avoid static interference
fringes. Power in both of the dipole trap beams is ramped up from zero to 3 W
in 500 ms. Simultaneously, the magnetic quadrupole trap strength is ramped
down from 74 G/cm to 7 G/cm in 500 ms thus the atoms are adiabatically
transferred to the optical trap. Three orthogonal pairs of Helmholtz coils are
used to cancel the effect of earth (and any other stray) magnetic fields on
the atoms and hold them stable at the center of the dipole trap. The dipole
trap center is slightly shifted from the magnetic quadrupole trap center
otherwise the green laser plugged in to the magnetic quadrupole trap center
would reduce the efficient loading of the dipole trap. The green laser is
turned off after the loading ramp is completed. After another 10 ms, the
magnetic trap strength is brought down to zero and the dipole trap alone now
holds the atoms. This completes the optical dipole trap loading.

Figure 4:  The lifetimes of the various mixture states of 39K and 87Rb. (a)
Lifetime of 39K . (b) Lifetime of 87Rb

At this stage, the number of 87Rb atoms is 8.21 $\times 10^{6}$ at 18 $\mu$K
while there are 2.15 $\times 10^{6}$ 39K atoms at 21 $\mu$K. The transfer
efficiency of the dipole trap from the quadrupole magnetic trap is around 20
percent for 39K and 17 percent for 87Rb. After the loading is complete, a
background magnetic field is also stabilized at 1 G along $z$-direction.

A forced evaporation in the optical trap then follows. The dipole trap power
is quickly reduced to 2.6 W in 30 ms to remove the hotter atoms in the wings
of the crossed dipole trap, and then the power decreases again to 1.6 W in 100
ms. These operations reduce the temperature of the sample down to 8 $\mu$K.
Subsequently, the 87Rb atoms are firstly transferred to the
$|F=1,m_{F}=1\rangle$ state adiabatically using a MW sweep of 100 ms duration.
Then a resonant flash light of 1 ms length removes the residual 87Rb atoms in
the $|2,2\rangle$ state. After another 1 ms of wait, the 39K atoms are also
transferred adiabatically to the $|F=1,m_{F}=1\rangle$ state using an RF sweep
in 100 ms and a similar flash light follows removing the remaining 39K atoms
in the $|2,2\rangle$ state. The transfer efficiency is almost 95 percent. In
this spin state preparation process, it is very important to transfer the 87Rb
atoms to the lower hyperfine level before 39K because the ground hyperfine
splitting of 87Rb is larger than that of 39K. Therefore, a mixture of
$|F=1,m_{F}=1\rangle$ state of 87Rb and $|F=2,m_{F}=2\rangle$ state of 39K is
more stable than the other way round which results in extreme losses due to
hyperfine changing collisions jarlt ; Thalhammer ; Marcassa .

To study the loss induced by this collision mechanism, we measure the lifetime
of the different spin state mixtures of the 39K and 87Rb atoms by holding
these mixtures in the same dipole trap depth for variable times and then
measuring the number of remaining atoms. In Fig. 4 the 39K (a) and 87Rb (b)
decay curves are shown for four combinations of the hyperfine states: (1)
Rb$|2,2\rangle$ $\oplus$ K$|2,2\rangle$, (2) Rb$|2,2\rangle$ $\oplus$
K$|1,1\rangle$, (3) Rb$|1,1\rangle$ $\oplus$ K$|2,2\rangle$, and (4)
Rb$|1,1\rangle$ $\oplus$ K$|1,1\rangle$.

The experimental data are fitted by the function $N_{0}$e-t/τ+$N_{r}$, with
$N_{r}$, $N_{0}$ and $t$ as free parameters ($N_{0}$ is the initial number of
atoms and $N_{r}$ is the atoms remaining after the long holding time), as the
solid lines shown in Fig. 4. The $\tau$ is the lifetime of the species in a
particular mixture. We notice that the mixture of the two species in case 2
has a very low initial 39K number $N_{0}$=0.97$\times 10^{5}$, which is one
order of magnitude lower than that in other cases. This is attributed to the
faster loss in the spin state transfer process for 39K. For case 2, we obtain
a very small lifetime constant of $\tau_{K}$=0.8 s and $\tau_{Rb}$=3.95 s. We
attribute this observed fast decay to the hyperfine changing collision, which
can be described by the equation below

$\begin{split}&{}^{39}K(F=1,m_{F}=1)+^{87}Rb(F=2,m_{F}=2)\rightarrow\\\
&{}^{39}K(F=2,m_{F}=2)+^{87}Rb(F=1,m_{F}=1)+\Delta,\end{split}$ (1)

where $\Delta$=h$\times$(6334.7-461.7) MHz $\approx$ $k_{B}$$\times$0.306 K.
The projection of total hyperfine angular momentum $m_{F,Rb}+m_{F,K}$ is
conserved in such a reaction. The released energy $\Delta$ in this process
converted to the kinetic energy of the collision atoms, is higher than the
trapping potential of the dipole trap and results in the fast atom loss. Other
hyperfine changing collisions can also happen Marcassa but the main
contribution to the losses in our system is coming from this reaction simoni ;
Thalhammer . When the spin preparation sequence is reversed (case 3), this
type of hyperfine changing collision is forbidden and longer lifetimes are
observed with $\tau_{K}$=1.7 s and $\tau_{Rb}$=7.92 s.

Figure 5:  The Feshbach resonance spectrum for the mixture 39K and 87Rb in the
$|1,1\rangle$ states. The upper panel shows the number of atoms remaining in
the trap for 39K and the lower panel shows that for 87Rb. The two intra-
species resonances at 25.9 G and 403.4 G for 39K and the single inter-species
resonance at 318.3 G for the mixture are presnted.

The measured most stable combination is in the case 4 with longest lifetimes
$\tau_{K}$=28 s and $\tau_{Rb}$=27 s. Because there is no any kind of
inelastic collision in this combination when both species are prepared in the
lowest states, the lifetime is only limited by the one-body loss induced by
the collision of background gases and the photon scattering of the dipole trap
beams. The 87Rb and 39K mixture when both in the $|2,2\rangle$ states is not
as stable (case 1) as the situation when each specie is in the $|1,1\rangle$
state (case 4) because in the former case there are still some residue atoms
in for example the $|2,1\rangle$ or other trappable state which can result in
inelastic collisional losses. Losses from other inelastic collisions reported
in Marcassa are also possible in this combination. Therefore, these losses
result in $\tau_{K}$=1.65 s and $\tau_{Rb}$=4.37 s in case 1.

All of these lifetime measurements are done while both of the species are not
condensed (at temperature of 8 $\mu$K). These decay lifetimes can be of
interest for the calculation of the interspecies scattering lengths among
different mixture states at ultra-cold temperatures.

##  V Feshbach resonances and achievement of dual BEC$\bf{s}$

Figure 6:  Absorption images of the 39K (upper row) and 87Rb (lower row) atoms
near the BEC critical temperature and the fits to the 1-D integrated optical
density along $x$-direction. The time of flight for these images is 20 ms for
39K and 30 ms for 87Rb. The field of view is 0.8 $\times$ 0.8 mm.

To further cool the two species, the 87Rb and 39K atoms are transferred to the
$|F=1,m_{F}=1\rangle$ state after following the exact sequence of state
preparation in the dipole trap and then we further perform the evaporation by
ramping down the dipole trap power. During these ramp down steps, the mixture
of 39K ($a_{K}$ $\sim$ -33$a_{0}$) and 87Rb ($a_{Rb}$ $\sim$ 100$a_{0}$) atoms
(interspecies $a_{RbK}$ $\sim$ 34$a_{0}$) are cooled by evaporation in the
presence of a homogeneous magnetic field of 1 G. The first ramp down takes 300
ms and the power in the dipole trap beams is reduced to 0.7 W. In another ramp
down lasting for 1 s, the power in the dipole trap beams is reduced to 0.35 W.
In the third ramp down, the dipole trap power in the two beams is reduced to
0.26 W in 500 ms. At this stage the temperature of the atoms reaches below 1
$\mu$K. At this temperature, we measure the inter- and intra-species Feshbach
resonances by scanning the external magnetic field along $z$-direction
provided by the pair of coils in the Helmholtz configuration (shown in Fig.
1). The resonances are detected as an enhancement of atoms loss wang2 ; Gerken
; zhangbei . Fig. 5 shows the results of the remaining atom number of both
species after holding time of 200 ms for different magnetic field values.

The upper panel shows two 39K intra-species resonances (25.9 G, 403.4 G) and a
single 87Rb-39K inter-species resonance (318.3 G) recorded by measuring the
number of 39K atoms after holding the two species in the magnetic field inside
the dipole trap for 200 ms (including 20 ms ramp up time of the magnetic
field). The corresponding number of 87Rb atoms is shown in the lower panel.
The single inter-species resonance is clearly visible at 318.3 G from the
simultaneous loss of both species. The other interesting feature in this
spectrum is the increase in the 87Rb atoms (encircled) at the exact 39K intra-
species resonances. This increase in the number of observed 87Rb atoms is due
to the reduced inter-species repulsion caused by the rapid loss of 39K.

When the atoms mixture is cooled to 0.8 $\mu$K, the negative background
scattering length of 39K significantly limit the further effective evaporation
cooling. Feshbach resonance can be used to tune the scattering length to avoid
this problem. Many interesting phenomena can be studied by the tuning of the
scattering length using Feshbach resonances wang3 ; yao1 ; zhangbei . To
achieve the BECs, the magnetic field of 372.6 G (which is below the 39K intra-
species resonance at 403.4 G) is ramped up in 20 ms and held 200 ms for the
field stabilization after the previously mentioned evaporation step of both
species in the $|F=1,m_{F}=1\rangle$ state. At this time, the scattering
length of 39K is $a_{K}$ = 20.05$a_{0}$, which is obtained form the well known
formula chiara ,

$\begin{split}a(B)=a_{bg}(1-\frac{\Delta}{B-B_{0}}),\end{split}$ (2)

where the background scattering length $a_{bg}$ = -29$a_{0}$, the Feshbach
resonance position $B_{0}$ = 403.4 G and the width $\Delta$ = -52 G. Then the
power in the dipole trap beams is reduced very slowly from 0.7 W to 0.2 W (in
2000 ms) and the final dual species BECs are produced.

To obtain the information of the BECs, the magnetic field is adiabatically
ramped to 354.2 G in 30 ms and held for 100 ms for tuning the 39K scattering
length close to zero. Then we perform the time of flight (TOF) imaging. The
dipole trap power is completely turned off to let the atom cloud to freely
expand for 5 ms in high magnetic field. And then we switch off the high
magnetic field and the atoms further freely expand for 20 ms for 39K and 30 ms
for 87Rb in the background field of 1 G before the absorption imaging. The
nearly pure condensates of 4.19 $\times 10^{5}$ 39 K atoms and 5.11 $\times
10^{5}$ 87 Rb atoms are produced.

Fig. 6 shows the absorption images of the optical density of the 39K and 87Rb
atoms, the Gaussian fitting for thermal component and the polynomial fitting
for the condensed fraction. The bimodal character appears in the second
columns for both species which means a BEC phase transition, and becomes
prominent in the third column of Fig. 6 for the almost pure condensates. The
first columns of Fig. 6 for both species show a single mode gaussian fit
representing a thermal gas.

##  VI Achievement of 39K BEC by single species evaporation

Figure 7:  Comparison of the 39K BEC achieved using the two different
processes. (a) Absorption image of the 39K BEC achieved using the two-species
evaporation performed at 372.6 G . (b) The same BEC produced using the single-
species evaporation performed at 391.46 G. (c, d) The respective poly-line
fits to these BECs. (e) The optical density of the BECs as a function of the
magnetic field in the two processes.

Next, we compare the evaporation efficiency of the two-species and the single-
species of 39K at different scattering lengths for 39K. Here the scattering
length $a_{Rb}$ of 87Rb is about 100$a_{0}$ and the inter-species scattering
$a_{KRb}$ is about 34$a_{0}$. For the two-species evaporation, we found that a
BEC with high number of atoms can be obtained in the broad region of positive
scattering length $a_{K}$ from about 5.7$a_{0}$ to 35.4$a_{0}$. However, for
only 39K evaporation the 39K BEC can be produced at 97.3$a_{0}$ with low atom
numbers as can be seen from Fig. 7. Here we remove the 87Rb atoms in the
$|2,2\rangle$ state using a flash light at 8 $\mu$K after the first two steps
of evaporation in the dipole trap described above.

##  VII Conclusion

In conclusion, we have achieved the dual 39K and 87Rb BECs by utilizing the
various available techniques of dark SPOT MOT, D1 line gray molasses, the
microwave evaporation and tuning the atomic scattering length using the
Feshbach resonances. We also measure the lifetimes of these various mixture
states showing interesting features especially the upper ground stretched
states lifetime is shorter than the lower hyperfine stretched states. The
difference in the lifetimes highlights the importance of the hyperfine
changing collisions in ultra-cold samples of dual species BECs. These dual
BECs with rich Feshbach resonance structure and large mass imbalance can be
used in various applications such as the formation of heteronuclear quantum
droplets. We also show how the sympathetic cooling of 39K using the 87Rb atoms
results in larger and denser BECs compared to the single species 39K
evaporative cooling.

## Acknowledgments

This research was supported by the MOST (Grants No.2016YFA0301602 and
No.2018YFA0307601), NSFC (Grants No.11974224, No.11704234, No.11804203,
No.12034011, No.12022406, No.12004229 and No.92065108), the Fund for Shanxi
“1331 Project” Key Subjects Construction, and the program of Youth Sanjin
Scholar.

## References

  * (1) Anderson M H, Ensher J R, Matthews M R, Wieman C E and Cornell E A 1995 Science 269 198
  * (2) Davis K B, Mewes M -O, Andrews M R, van Druten N J, Durfee D S, Kurn D M and Ketterle W 1995 Phys. Rev. Lett. 75 3969
  * (3) Giorgini S, Pitaevskii L P and Stringari S 2008 Rev. Mod. Phys. 80 1215
  * (4) Wang P, Yu Z, Fu Z, Miao J, Huang L, Chai S, Zhai H and Zhang J 2012 Phys. Rev. Lett. 109 095301
  * (5) Cheuk L W, Sommer A T, Hadzibabic Z, Yefsah T, Bakr W S and Zwierlein M W 2012 Phys. Rev. Lett. 109 095302
  * (6) Galitski V and Spielman I B 2013 Nature 494 49
  * (7) Wang D, Liu R, Zhu S and Scully M O 2015 Phys. Rev. Lett. 114 043602
  * (8) Chen L, Wang P, Meng Z, Huang L, Cai H, Wang D, Zhu S and Zhang J 2018 Phys. Rev. Lett. 120 193601
  * (9) Wang P, Chen L, Mi C, Meng Z, Huang L, Nawaz K S, Cai H, Wang D, Zhu S and Zhang J 2020 npj Quant. Inf. 6 18
  * (10) Liu Y, Hu M, Nichols M A, Grimes D D, Karman T, Guo H and Ni K 2020 Nat. Phys. 16 1132
  * (11) Nawaz K S, Chen L, Mi C, Chen L, Meng Z, Huang L, Wang P and Zhang J 2020 Phys. Rev. A 102 053326
  * (12) Grimm R 2019 Few-Body Syst. 51 113
  * (13) Tung S, Jiménez-García K, Johansen J, Parker C V and Chin C 2014 Phys. Rev. Lett. 113 240402
  * (14) Nawaz K S, Mi C, Chen L, Wang P and Zhang J 2019 Chin. Phys. Lett. 36 043201
  * (15) Bohn J L, Rey A M and Ye J 2017 Science 357 1002
  * (16) Ospelkaus S, Ni K -K, Wang D, de Miranda M H G, Neyenhuis B, Qu$\acute{e}$m$\acute{e}$ner G, Julienne P S, Bohn J L, Jin D S and Ye J 2010 Science 327 853
  * (17) Park J W, Will S A and Zwierlein M W 2015 Phys. Rev. Lett. 114 205302
  * (18) Schirotzek A, Wu C, Sommer A and Zwierlein M W 2009 Phys. Rev. Lett. 102 230402
  * (19) Kohstall C, Zaccanti M, Jag M, Trenkwalder A, Massignan P, Bruun G M, Schreck F and Grimm R 2012 Nature 485 615
  * (20) Yan Z Z, Ni Y, Robens C and Zwierlein M W 2020 Science 368 190
  * (21) Ferrier-Barbut I, Delehaye M, Laurent S, Grier A T, Pierce M, Rem B S, Chevy F and Salomon C 2014 Science 345 1035
  * (22) Wang F, Li X, Xiong D and Wang D 2015 J. Phys. B: At. Mol. Opt. Phys. 49 015302
  * (23) Wacker L, Jørgensen N B, Birkmose D, Horchani R, Ertmer W, Klempt C, Winter N, Sherson J and Arlt J J 2015 Phys. Rev. A 92 053602
  * (24) Schulze T A, Hartmann T, Voges K K, Gempel M W, Tiemann E, Zenesini A and Ospelkaus S 2018 Phys. Rev. A 97 023623
  * (25) Modugno G, Modugno M, Riboli F, Roati G and Inguscio M 2002 Phys. Rev. Lett. 89 190404
  * (26) Thalhammer G, Barontini G, De Sarlo L, Catani J, Minardi F and Inguscio M 2008 Phys. Rev. Lett. 100 210402
  * (27) Pasquiou B, Bayerle A, Tzanova S M, Stellmer S, Szczepkowski J, Parigger M, Grimm R and Schreck F 2013 Phys. Rev. A 88 023601
  * (28) Trautmann A, Ilzhöfer P, Durastante G, Politi C, Sohmen M, Mark M J and Ferlaino F 2018 Phys. Rev. Lett. 121 213601
  * (29) Lercher A D, Takekoshi T, Debatin M, Schuster B, Rameshan R, Ferlaino F, Grimm R and N$\ddot{a}$gerl H C 2011 Eur. Phys. J. D 65 3
  * (30) McCarron D J, Cho H W, Jenkin D L, Köppinger M P and Cornish S L 2011 Phys. Rev. A 84 011603
  * (31) Simoni A, Zaccanti M, $D^{\prime}$Errico C, Fattori M, Roati G, Inguscio M and Modugno G 2008 Phys. Rev. A 77 052705
  * (32) $D^{\prime}$Errico C, Zaccanti M, Fattori M, Roati G, Inguscio M, Modugno G and Simoni A 2007 New J. Phys. 9 223
  * (33) Prevedelli M, Cataliotti F S, Cornell E A, Ensher J R, Fort C, Ricci L, Tino G M and Inguscio M 1999 Phys. Rev. A 59 886
  * (34) Marcassa L G, Telles G D, Muniz S R and Bagnato V S 2000 Phys. Rev. A 63 013413
  * (35) Ketterle W, Davis K B, Joffe M A, Martin A and Pritchard D 1993 Phys. Rev. Lett. 70 2253
  * (36) Nath D, Easwaran R K, Rajalakshmi G and Unnikrishnan C S 2013 Phys. Rev. A 88 053407
  * (37) Salomon G, Fouche L, Wang P, Aspect A, Bouyer P and Bourdel T 2013 EPL 104 63002
  * (38) Sievers F, Kretzschmar N, Fernandes D R, Suchet D, Rabinovic M, Wu S, Parker C V, Khaykovich L, Salomon C and Chevy F 2015 Phys. Rev. A 91 023426
  * (39) Shi Z, Li Z, Wang P, Meng Z, Huang L and Zhang J 2018 Chin. Phys. Lett. 35 123701
  * (40) Chen L, Meng Z and Wang P 2017 Acta Phys. Sin. 66 083701
  * (41) Yang G, Chen L, Mi C, Wang P and Zhang J 2018 Acta Sin. Quantum Opt. 24 156
  * (42) Wang H, Gould P L and Stwalley W C 1995 Phys. Rev. A 53 R1216
  * (43) Fan H, Wang P and Zhang J 2015 Acta Sin. Quantum Opt. 21 351
  * (44) Landini M, Roy S, Carcagní L, Trypogeorgos D, Fattori M, Inguscio M and Modugno G 2011 Phys. Rev. A 84 043432
  * (45) Campbell R L D, Smith R P, Tammuz N, Beattie S, Moulder S and Hadzibabic Z 2010 Phys. Rev. A 82 063611
  * (46) Wang P, Fu Z, Chai S and Zhang J 2011 Chin. Phys. B 20 103401
  * (47) Gerken M, Tran B, Häfner S, Tiemann E, Zhu B and Weidemüler M 2019 Phys. Rev. A 100 050701(R)
  * (48) Zhang B, Li F, Deng S and Wu H 2019 Acta Sin. Quantum Opt. 25 413
  * (49) Wang J, Pan J, Cui X and Yi W 2020 Chin. Phys. Lett. 37 076701
  * (50) Yao R, Sun Z, Zhou S and Wang Y 2018 Chin. Phys. B 27 016703
  * (51)


# High Precision Control of Tracked Field Robots in the Presence of Unknown
Traction Coefficients

Erkan Kayacan
Distributed Autonomous Systems Laboratory
Coordinated Science Laboratory
University of Illinois at Urbana-Champaign
Urbana, IL 61801
<EMAIL_ADDRESS>&Sierra N. Young
Department of Civil and Environmental Engineering
University of Illinois at Urbana-Champaign
Urbana, IL 61801
<EMAIL_ADDRESS>Joshua M. Peschel
Dept. of Agricultural and Biosystems Engineering
Iowa State University
Ames, Iowa 50011-3270
<EMAIL_ADDRESS>&Girish Chowdhary
Distributed Autonomous Systems Laboratory
Coordinated Science Laboratory
Dept. of Agricultural and Biological Engineering
University of Illinois at Urbana-Champaign
Urbana, IL 61801
<EMAIL_ADDRESS>
###### Abstract

Accurate steering through crop rows that avoids crop damage is one of the most
important tasks for agricultural robots utilized in various field operations,
such as monitoring, mechanical weeding, or spraying. In practice, varying soil
conditions can result in off-track navigation due to unknown traction
coefficients so that it can cause crop damage. To address this problem, this
paper presents the development, application, and experimental results of a
real-time receding horizon estimation and control (RHEC) framework applied to
a fully autonomous mobile robotic platform to increase its steering accuracy.
Recent advances in cheap and fast microprocessors, as well as advances in
solution methods for nonlinear optimization problems, have made nonlinear
receding horizon control (RHC) and receding horizon estimation (RHE) methods
suitable for field robots that require high frequency (milliseconds) updates.
A real-time RHEC framework is developed and applied to a fully autonomous
mobile robotic platform designed by the authors for in-field phenotyping
applications in Sorghum fields. Nonlinear RHE is used to estimate constrained
states and parameters, and nonlinear RHC is designed based on an adaptive
system model which contains time-varying parameters. The capabilities of the
real-time RHEC framework are verified experimentally, and the results show an
accurate tracking performance on a bumpy and wet soil field. The mean values
of the Euclidean error and required computation time of the RHEC framework are
respectively equal to $0.0423$ m and $0.88$ milliseconds.

##  1 Introduction

Due to the petroleum crisis and growing demands for renewable energy sources
[Tverberg, 2012], there have been increased research efforts to improve
bioenergy crop breeding for biofuel production [Wang et al., 2016, Buckeridge
et al., 2014]. The success and sustainability of biofuels are dependent on
increasing plant yields [Vega-Sanchez and Ronald, 2010, Yano and Tuberosa,
2009]. Top-yielding plant traits (phenotypes) are linked to their
corresponding genes (genotypes). Identifying useful phenotypes can help in
identifying which genes will increase the productivity of biofuel crops.
However, the process of phenotyping is cost and labor intensive. On the other
hand, stationary high throughput phenotyping platforms, or greenhouse based
phenotyping systems are bulky, costly, and require elaborate infrastructure
[Fiorani and Tuberosa, 2013]. Therefore, there is a high need for an
automated, mobile, sensory platform that could significantly increase the
throughput and accuracy of screening crops in the field.

One of the major challenges in developing robotic phenotyping platforms is
high-precision guidance and control in muddy and uneven agricultural fields
[Kayacan et al., 2014a, Kayacan et al., 2015a]. A high-precision control
scheme is necessary to avoid crop disturbance and damage. For example, Figure
1 depicts the TERRA-MEPP phenotyping robot under development at the University
of Illinois at Urbana Champaign. It is a typical tracked robot that carries
hyperspectral, infrared, visual spectrum, and lidar sensors for phenotyping.
Due to the necessity to carry large sensing payloads, the width a phenotyping
robot can be significant, for example, it is $0.48$ m for the robot shown in
this study. However, the row spacing of sorghum cultivated for biofuel
production may vary from $0.5$ m to $1$ m, with a typical row spacing measured
centerline to centerline; of $0.78$ m as used in this study. Furthermore, the
stem width of a fully mature sorghum plant and the accuracy of the global
navigation satellite system used in this study can be up to $0.03$ m so that
the effective row spacing is reduced to $0.72$ m, which leaves only $0.12$ m
(4.5 in) of space on either side of the robot. Uncertain traction coefficients
and uneven fields can make it difficult for the robot to ensure that it stays
within the $0.12$ m tolerance so as to not damage the crop. It is the goal of
this paper to describe the development of high precision guidance and motion
control system that can ensure a navigational tolerance of strictly less than
$0.12$ m in real-world field environments.

Figure 1: Field robot in sorhgum plants in Energy Farm, Urbana, IL, USA.

Tracked robots are preferred to wheeled robots for field-based phenotyping
applications due to their improved traction and large contact area with the
ground, which minimizes adverse impacts on the soil [Bekker, 1956]. Tracked
robots, however, have complex track-ground interactions and slippage due to
differential velocities between treads, which can prompt difficulties in
control. The presence of these time-varying parameters can degrade the path-
tracking performance when traditional control approaches (e.g., proportional-
integral-derivative control) are employed [Fukao et al., 2000]. The reason is
that traditional controllers are not aware of soil parameters of unknown
terrain which play a vital role in determining field robot’s speed and
steering, which are in turn utilized for developing traction control
algorithms [Kayacan et al., 2012b]. Therefore, there is a need to develop
advanced control algorithms with online parameter estimation to update key
terrain and slip parameters, which can be used to improve robot mobility and
localization.

Field robots are multi-input-multi-output systems and must be aware of
actuator limitations [Kayacan et al., 2012a, Kayacan et al., 2015b, Kayacan
and Peschel, 2016, Kayacan, 2017]. Therefore, real-time applications of linear
model predictive controllers (MPCs), i.e., RHCs, have been used for path
tracking of mobile robots, but there are significant limitations of these
linear control techniques [Kayacan et al., 2015c]. Early linear MPCs have been
designed based on a tracking-error based linear model, and the total control
input consists of the feedback and feedforward control actions [Lee et al.,
2001, Gu and Hu, 2006, Lins Barreto et al., 2014]. In this scheme, the linear
MPC is the feedback controller, and the feedforward control action is derived
by taking the reference trajectory and system model into consideration
[Klancar and Skrjanc, 2007]. However, these linear control techniques are
unable to ensure accurate tracking performance if the systems start off-track
or there exist uncertainties, e.g., unknown terrain parameters. Linear control
techniques perform well for systems that work at fixed operating-points around
which the system is linearized. However, mobile robots work at varying
operating points (e.g., tracking of a curvilinear trajectory) and are
subjected to several uncertainties (e.g., bumpy fields, variable soil
conditions) so that locally linearized models are infeasible. Additionally,
linear MPCs cause large system errors in the prediction horizon with a
possible instability of the system because the mismatch between the linearized
model and the mobile robot increases when mobile robots move away from the
target path [Falcone et al., 2007]. To overcome these limitations, a robust
tube-based MPC has been developed and tested in real-time [Kayacan et al.,
2016b]. Although the system could stay-on track, a highly accurate trajectory
tracking performance has not been obtained.

To achieve highly accurate tracking performance for constrained navigation in
row crops, nonlinear RHC approaches have been widely used, particularly in
agricultural applications as they are capable of taking constraints on states
and inputs into account and can be designed for nonlinear models [Backman et
al., 2012, Utstumo et al., 2015, Kayacan et al., 2018]. Additionally, an
accurate online estimation of traction parameters in a system model is crucial
as the tire-soil interactions change throughout operations in varying,
slippery conditions. Thus, the system model can represent all the system
dynamics and interactions in the real-time system. Although nonlinear RHC
approach is computationally intensive, availability of compact and powerful
computational packages and advances in numerical methods for solving
optimization problems have allowed for real-time applications of constrained
nonlinear optimization problems for complex systems [Naveau et al., 2017, Oh
et al., 2015]. Motivated by the previous studies, this study employs a
nonlinear receding horizon estimation and control (RHEC) framework for path
tracking with a mobile, phenotyping robot.

The main contributions of this study are as follows.

  * •

In this work, we augment the kinematic model with traction parameters to
capture soil characteristics of agricultural fields, while the system model is
specified a priori and remains unchanged during operation in traditional
implementations.

  * •

We provide robust tracking performance when uncertainty is high as uncertainty
is reduced through learning whereas traditional RHC approaches do not
typically account for model uncertainty.

  * •

Highly tracking accuracy, e.g., less than $0.12$ m, for a tracked field robot
under unknown and variable soil conditions, has been obtained in real-time by
developing and implementing a nonlinear RHEC framework with a sampling time in
the range of milliseconds.

  * •

Computation times for the nonlinear RHEC have been decreased by restricting
the number of Gauss-Newton iterations to 1 while solving the nonlinear
optimization problems.

This paper is organized as follows: Section 2 provides the system model with
traction parameters and the description of the tracked mobile robot used in
this study. Section 3 is devoted to describing the nonlinear RHEC approaches.
Implementations and solution methods are presented in Section 4. Validation
and field tests for the real-time nonlinear RHEC algorithms are presented in
Section 5. Finally, conclusions are summarized in Section 6.

##  2 Field Robot

###  2.1 System Description

The field robot as shown in Fig. 2 has been built utilizing practical, hands-
on experience with modeling, control and estimation techniques, and various
sensors and actuators. It is $0.48$ m wide, $1.02$ m long and around $45$ kgs,
allowing it to navigate $0.78$ m wide crop rows without damaging the plants
and maintaining vertical stability. The mast on the robot can be adjusted to
match the sorghum canopy height and is used as a mounting platform for both
the plant imaging sensors and a Real-Time Kinematic (RTK) differential Global
Navigation Satellite System (GNSS) to acquire positional information. A
Septentrio Altus APS-NR2 GNSS receiver (Septentrio Satellite Navigation NV,
Belgium) is used to obtain highly accurate positional information, which has a
specified position accuracy of $0.03$ m at a $5$-Hz measurements rate. The
Trimble network supplies the RTK correction signals via the 4G internet. Two
DG-158 A DC motors (SuperDroid Robots Inc., U.S.) capable of $200$ W power
output with 16.95 Newton-meter rated torque are used to drive the tracks. E2
optical encoders (US Digital, USA) on DC motors provide speed information with
an accuracy of $0.035$ $m/s$. Two 12 V, 40 Ah 4s 20c lipo batteries power the
field robot and enable over 3 hours of operation.

Figure 2: Field robot in sorghum breeding plots in Energy farm, IL, USA.
Figure 3: Block diagram of hardware

The real-time estimation and control framework is implemented and executed on
an on-board computer (Raspberry Pi 3 Model B), which is equipped with $1.2$
GHz Quad Core Cortex-A53 64-Bit CPU, $1$ GB RAM. It acquires the GNSS data,
the speed of the wheels, and controls the tracked robot by sending signals to
the Kangaroo x2 motion controller (Dimension Engineering, USA). The generated
control input is sent to the motion controller via the on-board computer, and
speed measurements of the wheels are received with the same way. Positional
measurements are obtained via Bluetooth. The sampling frequency is set to
$5$-Hz due to the fact that the maximum update rate of the GNSS is equal to
$5$-Hz.

The block diagram of the hardware is illustrated in Fig. 3. The nonlinear
receding horizon estimation (RHE) receives the measurements (12) from onboard
sensors, and estimates the full states and parameters (6)-(8). The states and
parameters estimated by the nonlinear RHE and a reference trajectory are used
as inputs into the nonlinear RHC, which generates the control signal, i.e.,
the desired yaw rate. The desired yaw rate and speed as reference command
signals are then sent to the Kangaroo x2 Motion Controller that adds self-
tuning motion control to a Sabertooth dual 12A motor driver (Dimension
Engineering, USA), which allows us to control the speed of the DC motors. The
Kangaroo x2 Motion Controller, which is a two channel self-tuning PID
controller, functions as the field robot’s low-level controller by using
feedback from the encoders attached to the DC motors to determine the required
control signals. The Kangaroo x2 motion controller outputs the modified
command signals to Sabertooth dual 12A motor driver, which correlates the
given control signals to the necessary output voltages needed by the DC
motors. In other words, the inner loop, the speed of the DC motors, is
controlled by the Kangaroo x2 motion controller, which eliminates the need of
manually tune the low-level controller, and executed at a rate of $50$-Hz.

###  2.2 System Model with Traction Coefficients

The schematic diagram of the field robot, which is a typical tracked mobile
robot, is illustrated in Fig. 4. The velocities of two driven tracks result in
linear velocity $\nu=(\nu_{l}+\nu_{r})/2$ and angular velocity
$\omega=(\nu_{r}-\nu_{l})/b$ with the distance between tracks $b$. The
traditional kinematic model is formulated as follows [Kanayama et al., 1990]:

$\displaystyle\dot{x}$ $\displaystyle=$ $\displaystyle\nu\cos{\theta}$
$\displaystyle\dot{y}$ $\displaystyle=$ $\displaystyle\nu\sin{\theta}$
$\displaystyle\dot{\theta}$ $\displaystyle=$ $\displaystyle\omega$ (1)

where $x$ and $y$ denote the position of the field robot, $\theta$ denotes the
yaw angle, $\nu$ denotes the speed, $\omega$ denotes the yaw rate.

The effect of terrain characteristics on robot performance and track should be
well covered because the tracks are the unique connections between the ground
and robot, and nearly all forces and moments applied to the robot are
transmitted through the tracks [Iagnemma et al., 2004, Ray, 2009]. From track-
terrain interaction dynamics, soil characteristic plays a vital role in
determining vehicle speed and steering, which in turn are utilized for
developing traction control algorithms. The knowledge of soil parameters of
unknown terrain is then advantageous for improving vehicle performance [Lee
and Iagnemma, 2016, Lee et al., 2016]. Therefore, instead of using the
traditional kinematic model of a mobile robot (2.2), an adaptive nonlinear
kinematic model that is an extension of the traditional model is derived by
adding two traction parameters ($\mu,\kappa$) to minimize deviations between
the real-time system and system model of the field robot in this study. The
field robot can be formulated with the following equations:

$\displaystyle\dot{x}$ $\displaystyle=$ $\displaystyle\mu\nu\cos{\theta}$
$\displaystyle\dot{y}$ $\displaystyle=$ $\displaystyle\mu\nu\sin{\theta}$
$\displaystyle\dot{\theta}$ $\displaystyle=$ $\displaystyle\kappa\omega$ (2)

where $\mu$ and $\kappa$ denote the traction parameters. It is noted that they
must be between zero and one, and it is inherently arduous to measure them. If
the traction parameters are equal to 1, i.e., $\mu=\kappa=1$, there exist no
slip. The percentages of longitudinal and side slips are respectively found as
$1-\mu$ and $1-\kappa$. These parameters show the effective speed $\mu v$ and
steering $\kappa\omega$ of the field robot. To avoid bias, these traction
parameters must be estimated along with the full system state in each
iteration based on a number of past measurements.

Figure 4: Schematic illustration of the field robot

In the rest of the paper, we denote a nonlinear system model as

$\dot{\xi}(t)=f\Big{(}\xi(t),u(t),p(t)\Big{)}$ (3)

where $\xi$ $\in$ $\mathbb{R}^{n_{\xi}}$ is the state vector, $u$ $\in$
$\mathbb{R}^{n_{u}}$ is the control input, $p$ $\in$ $\mathbb{R}^{n_{p}}$ is
the system parameter vector,
$f(\cdot,\cdot,\cdot):\mathbb{R}^{n_{\xi}+n_{u}+n_{p}}\longrightarrow\mathbb{R}^{n_{\xi}}$
is the continuously differentiable state update function and
$f(0,0,p)=0\;\forall t$. The derivative of $\xi$ with respect to $t$ is
denoted by $\dot{\xi}$ $\in$ $\mathbb{R}^{n_{\xi}}$.

Similarly, a nonlinear measurement model denoted $z(t)$ can be described with
the following equation:

$z(t)=h\Big{(}\xi(t),u(t),p(t)\Big{)}$ (4)

where $h:R^{n_{\xi}+n_{u}+n_{p}}\longrightarrow R^{n_{z}}$ is the measurement
function which describes the relation between the variables of the system
model and the measured outputs of the real-time system.

The state, parameters, input and output vectors are respectively denoted as
follows:

$\displaystyle\xi$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{ccc}x&y&\theta\end{array}\right]^{T}$ (6)
$\displaystyle p$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{ccc}\nu&\mu&\kappa\end{array}\right]^{T}$
(8) $\displaystyle u$ $\displaystyle=$
$\displaystyle\begin{array}[]{c}\omega\end{array}$ (10) $\displaystyle z$
$\displaystyle=$
$\displaystyle\left[\begin{array}[]{cccc}x&y&\nu&\omega\end{array}\right]^{T}$
(12)

##  3 Nonlinear Receding Horizon Estimation and Control

###  3.1 Nonlinear Receding Horizon Estimation

The least squares estimation problem uses all available past measurements to
estimate the current values of the system states. In this estimation scheme,
the number of measurements increases and the optimization problem is enlarged
since a new measurement is acquired at every time instant. Therefore, solving
the relevant optimization problem becomes infeasible in real-time as time
approaches infinity. However, in practice, the computational burden of the
estimation problem must remain constant so that the number of past
measurements considered in the estimation scheme must also be constant.
Furthermore, the relationship between the system model with current estimated
parameters and measurements in the past can be inconsistent. To deal with this
problem, past measurements are deemed less influential or left out. The idea
of the RHE approach is to use only a constant number of the past measurements
so that the oldest data are removed in the moving window when new data is
acquired [Robertson et al., 1996]. Thus, a fixed size window will be moving
over the increasing measurement data so that the optimization problem will be
fixed size. Then, the required method to incorporate the pre-information from
the old measurements, which are not inside the moving window, and to prevent
the heavy influence of the old measurements is discussed in subsection 3.1.1.

State estimation alone is not enough to know the system behavior in real-time
when uncertain systems are considered; therefore, parameter estimation is
required to determine uncertain parameters. A parametric least square
estimation subject to the system model an/or boundary conditions was studied
extensively in the literature, and there are many software packages to solve
optimization problems for offline parameter estimation [Houska et al., 2011].
Moreover, two approximations have been proposed for online parameter
estimation, which is necessary to determine system behavior in conjunction
with state estimation accurately. In the first choice, model parameters are
assumed as so-called random constant by a differential equation
$\dot{\xi}_{p}=0$ with initial value $\xi_{p}(t_{k})=p_{k}$. This approach
results in time-invariant parameters over the estimation horizon. If jumps or
drifts for parameters are expected, which is the case in practice under
varying working conditions, a model bias will occur. To solve parameter jumps
and drifts problem, the model parameters must be assumed as time-varying.
Model parameters are assumed as so-called ”random walk” by a differential
equation $\dot{\xi}_{p}=\frac{d^{p}}{\triangle t}$ with sampling time
$\triangle t$ and initial value $\xi_{p}(t_{k})=p_{k}$ [Hughes, 1995].

In the nonlinear RHE method, the parameters are assumed to be time-invariant
in estimation horizon and not subject to process noise in the estimation
horizon. However, it is also assumed that the parameters are time-varying
Gaussian random variables in the arrival cost. For this reason, additional
weighting factors must be added as the variance of the parameters noise and
the parameters will appear only in the arrival cost. The extended weighting
matrix $W$ $\in$ $\mathbb{R}^{(n_{xi}+n_{p})\times(n_{\xi}+n_{p})}$ is written
as follows:

$\displaystyle W=diag(W_{\xi},W_{p})$ (13)

where $W_{\xi}$ $\in$ $\mathbb{R}^{n_{\xi}\times n_{\xi}}$ and $W_{p}$ $\in$
$\mathbb{R}^{n_{p}\times n_{p}}$ denote the weighting matrices for the state
noise and parameter pseudo-variance. Large values in $W$ result in large
variations in the corresponding state and parameter estimates, and vice-versa.

The nonlinear RHE is formulated as follows:

$\displaystyle\underset{\xi(t),p,u(t)}{\text{min}}$
$\displaystyle\frac{1}{2}\bigg{\\{}\left\|\begin{array}[]{c}\hat{\xi}-\xi(t_{k-N+1})\\\
\hat{p}-p\end{array}\right\|^{2}_{H_{N}}+\sum_{i=k-N+1}^{k}\|z_{m}(t_{i})-z(t_{i})\|^{2}_{H_{k}}\bigg{\\}}$
(14) s. t. $\displaystyle\dot{\xi}(t)=f\big{(}\xi(t),u(t),p\big{)}$
$\displaystyle z(t)=h\big{(}\xi(t),u(t),p\big{)}$
$\displaystyle\xi_{min}\leq\xi(t)\leq\xi_{max}$ $\displaystyle p_{min}\leq
p\leq p_{max}\qquad\qquad\forall t\in[t_{k-N},t_{k}]$

where the deviations of the oldest states in the estimation horizon and the
parameters from the prior estimates $\hat{\xi}$ and $\hat{p}$ are minimized by
a symmetric positive semi-definite weighting matrix $H_{N}$. Deviations of the
measured and system outputs in the estimation horizon are minimized by a
symmetric positive semi-definite weighting matrix $H_{k}$, which is equal to
the inverse of the measurement noise covariance matrix $V_{k}$ [Ferreau et
al., 2012]. As seen from the formulation, the objective function consists of
two parts: the arrival and quadratic costs. The arrival cost stands for the
early measurements $t=[t_{0,k-N+1}]$, and summarizes the deviations of the
system outputs with respect to the measured outputs before the beginning of
the estimation horizon. The quadratic cost stands for the recent measurements
$t=[t_{t-N+1,k}]$, and corresponds to deviations of the system outputs with
respect to the measured outputs within the estimation window. The estimation
horizon is denoted by $N$, and lower and upper constraints on the model states
and parameters are respectively denoted by $\xi_{min}$, $\xi_{max}$, $p_{min}$
and $p_{max}$.

####  3.1.1 Arrival Cost

The arrival cost in (14) is necessary for the stability of the nonlinear RHE.
The estimation horizon for the arrival cost is equal to $t_{k-N}$ and
approaches to infinity as time goes to infinity. Therefore, a real-time
solution for the RHE formulation will be infeasible after a finite time. In
order to forestall this problem, an iterative approach, a smoothed EKF update,
is proposed for the approximation of the arrival cost. It is run with
measurements at time $t_{k-N}$ and the inverted Kalman covariance matrix
$H_{N}$ in every iteration to update Kalman estimates $\hat{\xi},\hat{p}$. The
main advantage of using the inverted Kalman covariance matrix is to have an
adaptive weighting matrix. In case of using a constant weighting matrix, if
the deviations in the arrival cost are small and the weighting matrix is
large, then the chattering effect could be observed on the estimates. If the
deviations in the arrival cost are large and the weighting matrix is small,
then the estimates could not reach the true values in a reasonable time.

The reference estimated values for the states and parameters $\hat{\xi}$ and
$\hat{p}$ are acquired from the previous solution of RHE. Bounded arrival cost
is a requirement because if the inverted Kalman covariance is too heavy, the
arrival cost may approach infinity. Therefore, the contributions of the past
measurements to the inverted Kalman covariance $H_{N}$ are down-weighted by a
process noise covariance matrix $W$ in (13). It is to be noted that $H_{N}$ is
upper bounded by the inverse of $W$ since $W^{-1}-H_{N}$ is positive semi-
definite. The calculation of $H_{N}$ can be found in [Robertson, 1996].

###  3.2 Nonlinear Receding Horizon Control

RHC approach has the ability to anticipate the system behavior on a finite
horizon by minimizing a cost function consisting of the references, states,
and inputs and deal with hard constraints on states and inputs. A nonlinear
system model is described in (3). At every time step, the states and inputs
have to fulfill the following:

$\xi\in\mathbb{\Xi},\;\;\;u\in\mathbb{U}$ (15)

where $\mathbb{U}$ is compact subset of $\mathbb{R}^{n_{u}}$, i.e.,
$\mathbb{U}\subseteq\mathbb{R}^{n_{u}}$, $\mathbb{\Xi}$ is closed subset of
$\mathbb{R}^{n_{\xi}}$, i.e., $\mathbb{\Xi}\subseteq\mathbb{R}^{n_{\xi}}$, and
each set contains the origin in its interior. The constraints on each input
are uncoupled in the sense that the feasible regions of the inputs do not
affect each other.

The following formulation is solved at each sampling time $t$:

$\displaystyle\underset{\xi(t),u(t)}{\text{min}}$
$\displaystyle\frac{1}{2}\bigg{\\{}\Big{\\{}\sum_{i=k+1}^{k+N-1}\|\xi_{r}(t_{i})-\xi(t_{i})\|^{2}_{Q_{k}}+\|u_{r}(t_{i})-u(t_{i})\|^{2}_{R}\Big{\\}}$
(16)
$\displaystyle\qquad+\|\xi_{r}(t_{k+N})-\xi(t_{k+N})\|^{2}_{Q_{N}}\bigg{\\}}$
s. t. $\displaystyle\xi(t_{k})=\hat{\xi}(t_{k})$
$\displaystyle\dot{\xi}(t)=f\big{(}\xi(t),u(t),p\big{)}$
$\displaystyle\xi_{min}\leq\xi(t)\leq\xi_{max}\quad t\in[t_{k+1},t_{k+N}]$
$\displaystyle u_{min}\leq u(t)\leq u_{max}\quad t\in[t_{k+1},t_{k+N-1}]$

where $Q_{k}\in\mathbb{R}^{n_{\xi}\times n_{\xi}}$,
$R\in\mathbb{R}^{n_{u}\times n_{u}}$ and $Q_{N}\in\mathbb{R}^{n_{\xi}\times
n_{\xi}}$ are symmetric and positive semi-definite weighting matrices,
$\xi_{r}$ and $u_{r}$ are the state and input references, $\xi$ and $u$ are
the states and inputs, $t_{k}$ is the current time, $N$ is the prediction
horizon, $\hat{\xi}(t_{k})$ is the estimated state vector by the RHE,
$\xi_{min}$, $\xi_{max}$, $u_{min}$ and $u_{max}$ denote respectively the
upper and lower constraints on the state and input, $i=k+1,\ldots,k+N-1$,
$u(t)=[u(t_{k+1}),\ldots,u(t_{k+N-1}]$ is the input sequence over prediction
horizon $N$, and $\xi(t)$ is the state trajectory obtained by applying the
control sequence $u(t)$ to the system. The first element of the input sequence
$u(t)$ is applied to the systems:

$u(t_{k+1},\xi_{k+1})=u^{*}(t_{k+1})$ (17)

and the aforementioned nonlinear RHC formulation is solved again over a
shifted prediction horizon for the subsequent time instant. It is important to
point out that the control input $u^{*}(t_{k+1})$ is precisely the same as it
would be if all immeasurable states and parameters acquire values equal to
their estimates based on the estimation up to current time $t_{k}$ due to the
certainty equivalence principle.

The first cost is the stage cost, which is the cost throughout the prediction
horizon, and the second cost is the terminal penalty, which is the cost at the
end of the prediction horizon. The terminal penalty is stated for stability
reasons [Mayne et al., 2000, Rawlings and Mayne, 2009].

The formulation in (16) is a nonlinear and nonconvex optimization problem, and
the computational complexity depends on the order of the system, the
nonlinearity of the system, the horizon length and the nonlinear optimization
solver.

##  4 Implementations

###  4.1 Implementation of Nonlinear RHE

The inputs of the nonlinear RHE are the variables in the output vector (12),
which includes the position, yaw rate and speed of the field robot. The
outputs of the nonlinear RHE, the position, yaw angle, speed and traction
parameters, are the full state and parameter vectors (6)-(8). The nonlinear
RHC requires full state and parameter as an input; therefore, the estimated
values by the nonlinear RHE are fed to the nonlinear RHC to generate the
desired yaw rate applied to the field robot.

The nonlinear RHE formulation is solved at every sampling instant with the
following constraints on the traction parameters:

$\displaystyle 0\leq$ $\displaystyle\mu$ $\displaystyle\leq 1$ $\displaystyle
0\leq$ $\displaystyle\kappa$ $\displaystyle\leq 1$ (18)

The weighting matrix $H_{k}$ represents noise characteristics of sensors and
is determined using standard deviation of measurements. The measurements have
been perturbed by Gaussian noise with standard deviation of
$\sigma_{x}=\sigma_{y}=0.03$ m, $\sigma_{\omega}=0.0175$ $rad/s$,
$\sigma_{\nu}=0.05$ $m/s$ based on experimental analysis. Additionally, the
weighting matrix $W$ represents process noise, which is chosen based on the
objective. Low gain in the process noise results in better estimation
accuracy; however, it causes time-lag between true and estimated values.
Therefore, in this paper, the weighting coefficients for the measured states
and parameters (e.g., x an y positions, speed) are selected large while the
weighting coefficients for the immeasurable states and parameters (e.g., yaw
angle and traction parameters) are selected small. The following weighting
matrices $H_{k}$ and $W$ are used in nonlinear RHE:

$\displaystyle H_{k}=V_{k}^{-1}$ $\displaystyle=$ $\displaystyle
diag(\sigma_{x}^{2},\sigma_{y}^{2},\sigma_{\nu}^{2},\sigma_{\omega}^{2})^{-1}$
$\displaystyle=$ $\displaystyle
diag(0.03^{2},0.03^{2},0.05^{2},0.0175^{2})^{-1}$ $\displaystyle W$
$\displaystyle=$ $\displaystyle
diag(x^{2},y^{2},\theta^{2},\nu^{2},\mu^{2},\kappa^{2})$ (19) $\displaystyle=$
$\displaystyle diag(10.0^{2},10.0^{2},0.1^{2},1.0^{2},0.25^{2},0.25^{2})$

###  4.2 Implementation of Nonlinear RHC

The nonlinear RHC formulation is solved at every sampling instant with the
following constraint on the input:

$$-5.73$$\mathrm{\SIUnitSymbolDegree}\text{\,}{\mathrm{s}}^{-1}$\leq\omega(t)\leq$5.73$$\mathrm{\SIUnitSymbolDegree}\text{\,}{\mathrm{s}}^{-1}$$
(20)

The state and input references for the field robot are changed online and
defined as follows:

$\displaystyle\xi_{r}=[x_{r},y_{r},\theta_{r}]^{T}\quad\textrm{and}\quad
u_{r}=\omega_{r}$ (21)

where $x_{r}$ and $y_{r}$ are the position references, $\omega_{r}$ is the yaw
rate reference, and the yaw angle reference is calculated from the position
references as follows:

$\displaystyle\theta_{r}=\operatorname{atan2}{(\dot{y}_{r},\dot{x}_{r})}+\lambda\pi$
(22)

where $\lambda$ describes the desired direction of the field robot (
$\lambda=0$ for forward and $\lambda=1$ for backward). If the yaw rate
reference calculated from the reference trajectory is used as the input
reference, steady state error might occur in case of a mismatch between the
system model and real system. Therefore, the measured yaw rate is used as the
input reference to penalize the input rate in the objective function.

Three different weighting matrices $Q_{k}$ are selected to evaluate the yaw
angle effect on the trajectory tracking in Section 5.1.

$\displaystyle Q_{k}^{1}$ $\displaystyle=$ $\displaystyle diag(1,1,10)$ (23)
$\displaystyle Q_{k}^{2}$ $\displaystyle=$ $\displaystyle diag(1,1,0)$ (24)
$\displaystyle Q_{k}^{3}$ $\displaystyle=$ $\displaystyle diag(1,1,1)$ (25)

and the third one $Q_{k}^{3}$ is used in Section 5.2. The weighting matrices
$R$ and $Q_{N}$ for all experimental studies are selected as follows:

$\displaystyle R=10\quad\textrm{and}\quad Q_{N}=10\times Q_{k}$ (26)

The weighting matrix for the input $R$ is set larger than the one for the
states $Q_{k}$ to ensure well damped closed-loop system behavior. The
weighting matrix for the terminal penalty $Q_{N}$ is set to $10$ times larger
than the one for the states $Q_{k}$. Therefore, the last deviations between
the predicted states and their references in the prediction horizon are
minimized in the objective function $10$ times more than the previous points
in the prediction horizon. The reason is that the error at the end of the
prediction horizon plays a very critical role for the stability issue of the
control algorithm.

###  4.3 Extended Kalman Filter

An extended Kalman filter (EKF) is designed based on the traditional kinematic
model in (2.2). The discrete-time model used by the EKF is written with a
sampling interval $T_{s}$ as follows:

$\displaystyle x_{k+1}$ $\displaystyle=$ $\displaystyle
x_{k}+T_{s}\nu_{k}\cos{\theta_{k}}$ $\displaystyle y_{k+1}$ $\displaystyle=$
$\displaystyle y_{k}+T_{s}\nu_{k}\sin{\theta_{k}}$ $\displaystyle\theta_{k+1}$
$\displaystyle=$ $\displaystyle\theta_{k}+T_{s}\omega_{k}$ (27)

The general form of the estimated system model is:

$\displaystyle\widehat{\xi}_{k+1}$ $\displaystyle=$ $\displaystyle
f(\widehat{\xi}_{k},u_{k})+w_{k}$ $\displaystyle\widehat{z}_{k+1}$
$\displaystyle=$ $\displaystyle h(\widehat{\xi}_{k},u_{k})+v_{k}$ (28)

The differences between the kinematic model and the estimation model are the
process noise $w_{k}$ and the observation noise $v_{k}$ both in the state and
the measurement equations. They are both assumed to be independent with zero
mean multivariate Gaussian noises with covariance matrices $W_{k}$ and
$V_{k}$, respectively:

$\displaystyle w_{k}\backsim N(0,W_{k})$ $\displaystyle v_{k}\backsim
N(0,V_{k})$ (29)

The covariance matrix $V_{k}$ is formulated in (4.1), and the covariance
matrix $W_{k}$ is formulated as follows:

$\displaystyle W_{k}$ $\displaystyle=$ $\displaystyle
diag(x^{2},y^{2},\theta^{2},\nu^{2})$ (30) $\displaystyle=$ $\displaystyle
diag(10.0^{2},10.0^{2},0.1^{2},1.0^{2})$

The convariance matrix $W_{k}$ for the EKF is different than the covariance
matrix for the nonlinear RHE (4.1) because the traction parameters are not
estimated by the EKF. The EKF is designed based on the traditional kinematic
model (2.2) because it cannot deal with constraints on parameters. If the
traction parameters are estimated, their estimates might be lower than zero
[Haseltine and Rawlings, 2005], which results in instability of the control
algorithm.

###  4.4 Solution Methods

RHEC methods for systems require online solutions of nonlinear least square
optimization problems at each sampling time. In this study, single solution
method, consisting of a fusion between multiple shooting and generalized
Gauss-Newton methods, has been used to solve both nonlinear RHEC optimization
problems. This approach is valid because the formulation of nonlinear RHC
problem is akin to that of nonlinear RHE problem. The generalized Gauss-Newton
method, derived from the classical Newton method, was developed for least-
squared problems. This method is advantageous because it does not require
difficult computations of the second derivatives; however, it is challenging
to foreknow the number of iterations to reach a solution of the desired
accuracy. To overcome this challenge, the solution proposed in [Diehl et al.,
2002, KRA,] has been used where the number of Gauss-Newton iterations is
restricted to $1$, and the initial value of each optimization problem takes on
the value of the previous one intelligently. Hence this improves the
convergence of the Gauss-Newton method.

Figure 5: Illustration of the solution method

Unlike [Vukov et al., 2015], the Gauss-Newton iteration is divided into two
parts: preparation and feedback parts. As illustrated in Fig. 5, the
preparation part is executed prior to the feedback part, and the feedback part
is executed after measurements for RHE and estimates for RHC are available. In
the preparation part, the system dynamics are integrated with the previous
solution, and objectives, constraints, and corresponding sensitives are
evaluated. In the feedback part, a single quadratic programming is solved with
the current measurements for the RHE and the current estimates for the RHC.
Thus, the new estimates for the RHE and the new control signal for the RHC are
obtained. Compared to the classical method, this method minimizes feedback
delay and produces similar results with higher computational efficiency [Diehl
et al., 2002].

To solve the constrained nonlinear optimization problems in the nonlinear
RHEC, the _ACADO_ code generation tool was used. _ACADO_ is an open source
software package for optimization problems that generates C-code [Houska et
al., 2011]. The obtained dense quadratic problem sub-problems are solved by
the online quadratic problem solver [Ferreau et al., 2008].

##  5 Experimental Results

The aim of the real-time experiments is to track a predefined trajectory with
a field robot exposed to varying soil conditions characterized by a wet and
bumpy field. Highly accurate guidance is required to avoid crop, e.g.,
sorghum, damage during planting, weeding, and growth; however, slippage makes
this guidance problem difficult. In this study, a space-based trajectory
approach is employed rather than a time-based trajectory approach. The reason
is that the speed of the field robot is featured with large uncertainty
arising from the variation of terrain, load, turning resistance and
surrounding obstacles. These uncertainties can impose large time deviations,
which challenge the tracking of the reference time-based trajectory. However,
it is to be pointed out that the field robot is not forced to be at a
particular point on the trajectory at a certain sampling instant in space-
based trajectory approach.

The reference generation method in this paper is as follows. If the field
robot starts off-track, the closest point on the target trajectory is computed
first, and then the desired points in the horizon are determined. A trajectory
consisting of straight and curved lines is tracked so that the performance of
the designed framework for both scenarios is investigated. The prediction and
estimation horizons for the nonlinear RHEC are set to 3 seconds.

###  5.1 Scenario A: Effect of Yaw Angle Weight on RHC

In this subsection, the coefficient of the yaw angle in the weighting matrix
for nonlinear RHC is evaluated. Trajectory tracking performance of nonlinear
RHCs based on nonlinear RHE is shown in Fig. 6. It is evident that if the
coefficient of the yaw angle is larger than the coefficients for the position
as in (23), the field robot cannot reach the desired trajectory after a
reasonable time. However, if the coefficient of the yaw angle is set to zero
as in (24), in other words, the yaw angle error between the reference yaw
angle and actual yaw angle is not minimized, the field robot oscillates around
the reference trajectory. To overcome this problem, the beginning of the
desired points is defined as a fixed forward distance from the closest point
on the reference trajectory at every sampling instant, and the optimum
distance was selected regarding the speed of vehicles [Kayacan et al., 2016a].

Figure 6: Reference and actual trajectories for different weighting of the yaw
angle. These results clearly demonstrate that the penalization of yaw angle
error is extremely important in RHC formulation. The presented approach in the
paper should be contrasted with the typical approach of choosing small
coefficients for R matrix, which essentially leads to low-gain control. By
minimizing yaw angle error in the cost function in RHC formulation explicitly,
large coefficients for R matrix can be employed, which leads to better
tracking performance.

In this paper, the minimization of the yaw angle error in the nonlinear RHC
design is implemented to solve the aforementioned limitations. The coefficient
of the yaw angle is set to equal to the coefficients for the position as in
(25). It can be seen from Fig. 6 that the field robot reaches the reference
trajectory after a reasonable time and stay on-track, which highlights the
importance of correctly weighting of the yaw angle in the cost function in the
space-based trajectory approach.

###  5.2 Scenario B: Comparison of RHCs based on RHE and EKF

Motivated by results from the previous section, the yaw angle error is
minimized in the nonlinear RHC formulations based on the nonlinear RHE and the
extended Kalman filter (EKF). As can be seen from Fig. 7, the field robot
reaches the reference trajectory after it is started off-track and stays-on
track for nonlinear RHCs based on both the nonlinear RHE and the EKF.

Figure 7: Reference and actual trajectories. The Euclidian errors are plotted
in Fig. 8.

The Euclidean errors calculated using raw GNSS data to the space-based
reference trajectory for nonlinear RHCs based on both the nonlinear RHE and
EKF are shown in Fig. 8. The mean values of the Euclidian error for the
nonlinear RHCs based on the nonlinear RHE and the EKF are $0.0423$ m and
$0.0514$ m, respectively. The nonlinear RHEC framework benefits from traction
parameter estimates and results in less error compared to the RHC based on the
EKF. As previously mentioned in Section 1, the available space on either side
of the robot is limited to only $0.12$ m; therefore, the error must be smaller
than this limit to avoid crop damage and keep the robot-centered in the row.
The number of violations is shown in Fig. 9. The results of multiple
experiments indicate that the RHEC framework does not violate this error
constraint, while the RHC based on the EKF violates it 17 times during the
path tracking of straight lines. This demonstrates the capability of the RHEC
framework.

Figure 8: Euclidian error calculated using raw GNSS data
Figure 9: Number of violations for straight lines tracking. The available
space on either side of the field robot is restricted to $0.12$ m. The field
robot controlled by the RHEC framework does not violate this limit; therefore,
it does not cause crop damage.

The nonlinear RHE performance for the yaw angle and traction coefficients is
shown in Fig. 10. It is difficult to measure the yaw angle in practice, but it
plays a vital role in the trajectory tracking problem. If it is not estimated,
then the system model is deteriorated so that model-based controllers cannot
be employed. As can be seen, the estimated values for the traction parameters
are within the constraints specified in (4.1) and consistently stay close to
the upper bound. It is observed that model parameter estimates stabilize at
certain values which secure stable trajectory-tracking. Therefore, these
traction parameters’ estimates play a significant role in the nonlinear RHC
performance, and consequently the trajectory tracking performance. It is to be
noted that hard-coding can result in more accurate trajectory tracking
performance. However, the soil conditions must be in a steady-state, and the
identification of parameters is required, which is a time-consuming task
[Kayacan et al., 2014b, Kayacan et al., 2013]. For this reason, the nonlinear
RHE is beneficial because it provides the estimates of the traction
parameters, and the framework ensures an accurate trajectory tracking.

Figure 10: Nonlinear RHE output: The yaw angle and traction parameters
estimates

The measured and estimated speeds are shown in Fig. 11. The speed of the field
robot is constant because a space-based trajectory approach is employed in the
paper, and the nonlinear RHE can cope with noise on the speed measurements.

Figure 11: Actual and estimated speeds

The output of the nonlinear RHC, the reference yaw rate, is shown in Fig. 12.
The output of the nonlinear RHC reaches the bounds, but never violates them.
This demonstrates the capability of dealing with constraints on inputs.
Moreover, the performance of the low-level controller for the speed of the DC
motors is sufficient.

Figure 12: Reference and actual yaw rates

Real-time algorithms for RHE and RHC comprise respectively three states, two
parameters, one control input and 15 prediction intervals, and three states,
one control input, and 15 prediction intervals. The computations times for the
nonlinear RHEC are summarized in Table 1. As seen from the table, the average
computation times for the nonlinear RHE and RHC are equal to $0.48$ ms and
$0.40$ ms, respectively. Thus, the overall computation time for the nonlinear
RHEC framework is equal to $0.88$ ms with a reasonable maximum computation
time of $2.8534$ ms for real-time applications.

Table 1: Computation times of the nonlinear RHEC. | Minimum | Average | Maximum
---|---|---|---
| (ms) | (ms) | (ms)
Nonlinear RHE |  |  |
Preparation | 0.2921 | 0.3523 | 1.0075
Feedback | 0.0766 | 0.1309 | 0.4925
Overall | 0.3687 | 0.4832 | 1.5000
Nonlinear RHC |  |  |
Preparation | 0.2725 | 0.3751 | 1.2077
Feedback | 0.0149 | 0.0215 | 0.1457
Overall | 0.2874 | 0.3966 | 1.3534

A single quadratic programming iteration at each time instant may result in a
sub-optimal solution; therefore, it is necessary to observe the Karush-Kuhn-
Tucker (KTT) tolerance, which measures the performance of RHEC methods in
terms of optimality. The KKT tolerances for the RHEC are shown in Fig. 13. It
is to be noted that the KKT tolerance is equal to zero for linear systems due
to the fact that a quadratic programming at each iteration is solved
precisely. Therefore, non-zero values demonstrate the non-linearity of the
RHEC problems. The largest KKT tolerances for the RHC arise at the beginning
of the experiments because the field robot initially started off-track in the
sense that error values are large. The reason is that since a single iteration
is completed, the obtained solution is only a rough approximation to the
optimal solution. It can be seen that the KKT tolerance decreases with time
until the field robot stays on-track. The trajectory references continuously
change throughout the experiment; therefore, the KKT tolerance cannot further
decrease. Furthermore, although the KKT tolerance for the RHE is large at the
beginning due to poor initialization, it decreases after several time periods.

Figure 13: KKT tolerances

##  6 Conclusion

A computationally efficient, nonlinear RHEC framework based on an adaptive
nonlinear model has been investigated for the guidance of the field robot. The
experimental results in the field trials have shown that the nonlinear RHE can
perform precise online estimation of the immeasurable states and parameters.
Also, the nonlinear RHC has the capability of accurately controlling the
field, phenotyping robot with reasonable accuracy. The importance of correctly
weighting of the yaw angle in the cost function is highlighted for the space-
based trajectory approach. The mean value of the Euclidean error to the
reference trajectory is $0.0423$ m. Additionally, the RHEC framework does not
violate the upper bound of the required tracking accuracy, i.e., $0.12$ m,
while the RHC based on the EKF violates it 17 times for the tracking of
straight lines. Although RHEC framework is computationally intensive, the
proposed solution method is computationally efficient, with a mean total
required computation time of $0.88$ ms and the worst-case computation time of
$2.85$ ms.

The high accuracy of this nonlinear RHEC framework is of importance to plant
breeders and farmers because there is a growing need for rapid, mobile, field-
based phenotyping technologies that can navigate safely and autonomously
through row based crops. By developing a framework that can handle unknown and
immeasurable slip parameters, field robots will be able to navigate a variety
of agricultural field conditions without damaging the crop.

Recent developments in microprocessors technology and fast solution tools for
RHEC framework have changed the well-known paradigm in a way that the belief
of using RHEC for only relatively slow dynamic systems is no longer true. The
comparative results presented in this paper also show that RHEC
implementations for agile robotic systems are eligible to obtain highly
accurate tracking performance.

#### Acknowledgments

The information, data, or work presented herein was funded in part by the
Advanced Research Projects Agency-Energy (ARPA-E), U.S. Department of Energy,
under Award Number DE-AR0000598. The views and opinions of authors expressed
herein do not necessarily state or reflect those of the United States
Government or any agency thereof.

Dr. Erkan Kayacan is a postdoctoral researcher at Massachusetts Institute of
Technology since February 2018.

## References

  * [KRA, ]
  * [Backman et al., 2012] Backman, J., Oksanen, T., and Visala, A. (2012).  Navigation system for agricultural machines: Nonlinear model predictive path tracking.  Computers and Electronics in Agriculture, 82:32 – 43.
  * [Bekker, 1956] Bekker, M. (1956).  Theory of Land Locomotion.  University of Michigan Press, Ann Arbor, MI.
  * [Buckeridge et al., 2014] Buckeridge, M. S., Carpita, N. C., McCann, M. C., et al. (2014).  Plants and bioenergy.  Springer,.
  * [Diehl et al., 2002] Diehl, M., Bock, H., Schlöder, J. P., Findeisen, R., Nagy, Z., and Allgöwer, F. (2002).  Real-time optimization and nonlinear model predictive control of processes governed by differential-algebraic equations.  Journal of Process Control, 12(4):577 – 585.
  * [Falcone et al., 2007] Falcone, P., Borrelli, F., Asgari, J., Tseng, H. E., and Hrovat, D. (2007).  Predictive active steering control for autonomous vehicle systems.  Control Systems Technology, IEEE Transactions on, 15(3):566–580.
  * [Ferreau et al., 2012] Ferreau, H., Kraus, T., Vukov, M., Saeys, W., and Diehl, M. (2012).  High-speed moving horizon estimation based on automatic code generation.  In Decision and Control (CDC), 2012 IEEE 51st Annual Conference on, pages 687–692.
  * [Ferreau et al., 2008] Ferreau, H. J., Bock, H. G., and Diehl, M. (2008).  An online active set strategy to overcome the limitations of explicit mpc.  International Journal of Robust and Nonlinear Control, 18(8):816–830.
  * [Fiorani and Tuberosa, 2013] Fiorani, F. and Tuberosa, R. (2013).  Future scenarios for plant phenotyping.  Annu Rev Plant Biol, 64:267–291.
  * [Fukao et al., 2000] Fukao, T., Nakagawa, H., and Adachi, N. (2000).  Adaptive tracking control of a nonholonomic mobile robot.  IEEE Transactions on Robotics and Automation, 16(5):609–615.
  * [Gu and Hu, 2006] Gu, D. and Hu, H. (2006).  Receding horizon tracking control of wheeled mobile robots.  Control Systems Technology, IEEE Transactions on, 14(4):743–749.
  * [Haseltine and Rawlings, 2005] Haseltine, E. L. and Rawlings, J. B. (2005).  Critical evaluation of extended kalman filtering and moving-horizon estimation.  Industrial & Engineering Chemistry Research, 44(8):2451–2460.
  * [Houska et al., 2011] Houska, B., Ferreau, H. J., and Diehl, M. (2011).  Acado toolkit—an open-source framework for automatic control and dynamic optimization.  Optimal Control Applications and Methods, 32(3):298 – 312.
  * [Hughes, 1995] Hughes, B. (1995).  Random walks and random environments.  Oxford, 1.
  * [Iagnemma et al., 2004] Iagnemma, K., Kang, S., Shibly, H., and Dubowsky, S. (2004).  Online terrain parameter estimation for wheeled mobile robots with application to planetary rovers.  IEEE Transactions on Robotics, 20(5):921–927.
  * [Kanayama et al., 1990] Kanayama, Y., Kimura, Y., Miyazaki, F., and Noguchi, T. (1990).  A stable tracking control method for an autonomous mobile robot.  In Proceedings., IEEE International Conference on Robotics and Automation, pages 384–389.
  * [Kayacan, 2017] Kayacan, E. (2017).  Multiobjective $h_{\infty}$ control for string stability of cooperative adaptive cruise control systems.  IEEE Transactions on Intelligent Vehicles, 2(1):52–61.
  * [Kayacan et al., 2012a] Kayacan, E., Bayraktaroglu, Z. Y., and Saeys, W. (2012a).  Modeling and control of a spherical rolling robot: a decoupled dynamics approach.  Robotica, 30(4):671–680.
  * [Kayacan et al., 2018] Kayacan, E., Kayacan, E., Chen, I.-M., Ramon, H., and Saeys, W. (2018).  On the Comparison of Model-Based and Model-Free Controllers in Guidance, Navigation and Control of Agricultural Vehicles, pages 49–73.  Springer International Publishing, Cham.
  * [Kayacan et al., 2015a] Kayacan, E., Kayacan, E., Ramon, H., Kaynak, O., and Saeys, W. (2015a).  Towards agrobots: Trajectory control of an autonomous tractor using type-2 fuzzy logic controllers.  IEEE/ASME Transactions on Mechatronics, 20(1):287–298.
  * [Kayacan et al., 2012b] Kayacan, E., Kayacan, E., Ramon, H., and Saeys, W. (2012b).  Velocity control of a spherical rolling robot using a grey-pid type fuzzy controller with an adaptive step size.  IFAC Proceedings Volumes, 45(22):863 – 868.  10th IFAC Symposium on Robot Control.
  * [Kayacan et al., 2013] Kayacan, E., Kayacan, E., Ramon, H., and Saeys, W. (2013).  Modeling and identification of the yaw dynamics of an autonomous tractor.  In 2013 9th Asian Control Conference (ASCC), pages 1–6.
  * [Kayacan et al., 2014a] Kayacan, E., Kayacan, E., Ramon, H., and Saeys, W. (2014a).  Distributed nonlinear model predictive control of an autonomous tractor-trailer system.  Mechatronics, 24(8):926 – 933.
  * [Kayacan et al., 2014b] Kayacan, E., Kayacan, E., Ramon, H., and Saeys, W. (2014b).  Nonlinear modeling and identification of an autonomous tractor-trailer system.  Computers and Electronics in Agriculture, 106(Supplement C):1 – 10.
  * [Kayacan et al., 2015b] Kayacan, E., Kayacan, E., Ramon, H., and Saeys, W. (2015b).  Robust tube-based decentralized nonlinear model predictive control of an autonomous tractor-trailer system.  IEEE/ASME Transactions on Mechatronics, 20(1):447–456.
  * [Kayacan et al., 2015c] Kayacan, E., Kayacan, E., Ramon, H., and Saeys, W. (2015c).  Towards agrobots: Identification of the yaw dynamics and trajectory tracking of an autonomous tractor.  Computers and Electronics in Agriculture, 115:78 – 87.
  * [Kayacan and Peschel, 2016] Kayacan, E. and Peschel, J. (2016).  Robust model predictive control of systems by modeling mismatched uncertainty.  IFAC-PapersOnLine, 49(18):265 – 269.  10th IFAC Symposium on Nonlinear Control Systems NOLCOS 2016.
  * [Kayacan et al., 2016a] Kayacan, E., Peschel, J. M., and Kayacan, E. (2016a).  Centralized, decentralized and distributed nonlinear model predictive control of a tractor-trailer system: A comparative study.  In 2016 American Control Conference (ACC), pages 4403–4408.
  * [Kayacan et al., 2016b] Kayacan, E., Ramon, H., and Saeys, W. (2016b).  Robust trajectory tracking error model-based predictive control for unmanned ground vehicles.  IEEE/ASME Transactions on Mechatronics, 21(2):806–814.
  * [Klancar and Skrjanc, 2007] Klancar, G. and Skrjanc, I. (2007).  Tracking-error model-based predictive control for mobile robots in real time.  Robotics and Autonomous Systems, 55(6):460 – 469.
  * [Lee et al., 2016] Lee, S. U., Gonzalez, R., and Iagnemma, K. (2016).  Robust sampling-based motion planning for autonomous tracked vehicles in deformable high slip terrain.  In 2016 IEEE International Conference on Robotics and Automation (ICRA), pages 2569–2574.
  * [Lee and Iagnemma, 2016] Lee, S. U. and Iagnemma, K. (2016).  Robust motion planning methodology for autonomous tracked vehicles in rough environment using online slip estimation.  In 2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pages 3589–3594.
  * [Lee et al., 2001] Lee, T.-C., Song, K.-T., Lee, C.-H., and Teng, C.-C. (2001).  Tracking control of unicycle-modeled mobile robots using a saturation feedback controller.  Control Systems Technology, IEEE Transactions on, 9(2):305–318.
  * [Lins Barreto et al., 2014] Lins Barreto, J. C., Scolari Conceicao, A. G., Dorea, C. E. T., Martinez, L., and De Pieri, E. R. (2014).  Design and implementation of model-predictive control with friction compensation on an omnidirectional mobile robot.  IEEE/ASME Transactions on Mechatronics, 19(2):467–476.
  * [Mayne et al., 2000] Mayne, D., Rawlings, J., Rao, C., and Scokaert, P. (2000).  Constrained model predictive control: Stability and optimality.  Automatica, 36(6):789 – 814.
  * [Naveau et al., 2017] Naveau, M., Kudruss, M., Stasse, O., Kirches, C., Mombaur, K., and Souères, P. (2017).  A reactive walking pattern generator based on nonlinear model predictive control.  IEEE Robotics and Automation Letters, 2(1):10–17.
  * [Oh et al., 2015] Oh, H., Kim, S., and Tsourdos, A. (2015).  Road-map assisted standoff tracking of moving ground vehicle using nonlinear model predictive control.  IEEE Transactions on Aerospace and Electronic Systems, 51(2):975–986.
  * [Rawlings and Mayne, 2009] Rawlings, J. B. and Mayne, D. Q. (2009).  Model predictive control: Theory and design.  Nob Hill Pub.
  * [Ray, 2009] Ray, L. E. (2009).  Estimation of terrain forces and parameters for rigid-wheeled vehicles.  IEEE Transactions on Robotics, 25(3):717–726.
  * [Robertson, 1996] Robertson, D. (1996).  Development and Statistical Interpretation of Tools for Nonlinear Estimation.  Auburn University.
  * [Robertson et al., 1996] Robertson, D. G., Lee, J. H., and Rawlings, J. B. (1996).  A moving horizon-based approach for least-squares estimation.  AIChE Journal, 42(8):2209–2224.
  * [Tverberg, 2012] Tverberg, G. (2012).  Oil supply limits and the continuing financial crisis.  Energy, 37(1):27–34.
  * [Utstumo et al., 2015] Utstumo, T., Berge, T. W., and Gravdahl, J. T. (2015).  Non-linear model predictive control for constrained robot navigation in row crops.  In 2015 IEEE International Conference on Industrial Technology (ICIT), pages 357–362.
  * [Vega-Sanchez and Ronald, 2010] Vega-Sanchez, M. and Ronald, P. (2010).  Genetic and biotechnological approaches for biofuel crop improvement.  Curr Opin Plant Biol, 21(2):218–224.
  * [Vukov et al., 2015] Vukov, M., Gros, S., Horn, G., Frison, G., Geebelen, K., Jørgensen, J., Swevers, J., and Diehl, M. (2015).  Real-time nonlinear mpc and mhe for a large-scale mechatronic application.  Control Engineering Practice, 45:64 – 78.
  * [Wang et al., 2016] Wang, Y., Fan, C., Hu, H., Li, Y., Sun, D., Wang, Y., and Peng, L. (2016).  Genetic modification of plant cell walls to enhance biomass yield and biofuel production in bioenergy crops.  Biotechnology Advances, 34(5):997 – 1017.
  * [Yano and Tuberosa, 2009] Yano, M. and Tuberosa, R. (2009).  Genome studies and molecular genetics-from sequence to crops: genomics comes of age.  Curr Opin Plant Biol, 12(2):103–106.