Received July 19, 2021, accepted September 6, 2021, date of publication
September 23, 2021, date of current version September 30, 2021.
10.1109/ACCESS.2021.3115082

Corresponding author: Fanjun Bu (e-mail: fb266@cornell.edu).

# Object Permanence Through Audio-Visual Representations

FANJUN BU1  CHIEN-MING HUANG 2  Department of Computer Science, Cornell
University, Ithaca, NY 14850 USA (e-mail<EMAIL_ADDRESS> Department of
Computer Science, Johns Hopkins University, Baltimore, MD 21218 USA (e-mail:
<EMAIL_ADDRESS>
###### Abstract

As robots perform manipulation tasks and interact with objects, it is probable
that they accidentally drop objects (e.g., due to an inadequate grasp of an
unfamiliar object) that subsequently bounce out of their visual fields. To
enable robots to recover from such errors, we draw upon the concept of _object
permanence_ —objects remain in existence even when they are not being sensed
(e.g., seen) directly. In particular, we developed a multimodal neural network
model—using a partial, observed bounce trajectory and the audio resulting from
drop impact as its inputs—to predict the full bounce trajectory and the end
location of a dropped object. We empirically show that: 1) our multimodal
method predicted end locations close in proximity (i.e., within the visual
field of the robot’s wrist camera) to the actual locations and 2) the robot
was able to retrieve dropped objects by applying minimal vision-based pick-up
adjustments. Additionally, we show that our method outperformed five
comparison baselines in retrieving dropped objects. Our results contribute to
enabling object permanence for robots and error recovery from object drops.

###### Index Terms:

Error Recovery, Multimodal Neural Network, Object Localization, Object
Permanence, Trajectory Prediction.

=-15pt

##  I Introduction

We all drop objects, from car keys to pens or utensils. When an object is
dropped and bounces out of sight, people retrieve the object by estimating its
landing position. This ability to retrieve objects that are outside of visual
fields is based on an understanding of _object permanence_ [1]—objects remain
in existence even when they may not be visible. As robots increasingly
interact with human-made objects, they are bound to drop objects. Drops can be
frequent not only due to the challenging nature of robot grasping [2],
imperfect sensing, complex human-made objects in unstructured human
environments, and computationally intensive pose estimations, but also because
robots may accidentally drop objects during manipulation. This paper addresses
the problem of _how robots may recover from dropping objects_ , focusing on
how to retrieve objects that bounce out of the robot’s visual field by
estimating bounce trajectories and object locations (Fig. 1).

Figure 1: We explore object permanence through audio-visual representations
and contextualize our exploration in retrieving dropped objects that bounce
out of its visual field.

Drawing inspiration from humans’ ability to locate objects using multiple
sensing modalities, we developed a multimodal deep network that takes as input
a partial, observed trajectory and impact sound to predict an object’s bounce
trajectory and end location. Our multimodal network encodes modality-specific
features from visual and auditory channels and uses observed partial
trajectories to adjust audio-driven spatial predictions enabled by a
microphone array. To illustrate our approach to object permanence, we used a
simple wooden object and compared our approach to physics-driven baselines and
baselines informed by previous state-of-the-art models. We additionally
explored the possibility of generalizing our approach to a different object
and dropping height. However, rather than focusing on developing a universal
model capable of predicting end locations for objects made of various
materials under a diversity of conditions, our goal in this paper is to
demonstrate the utility of leveraging audio input and incomplete visual data
to achieve object permanence.

Our contributions include: 1) a human-inspired multimodal approach to object
permanence, 2) a demonstration of how a robotic manipulator may use the
multimodal approach to recover from object dropping, and 3) an open dataset of
object drops for the robotics community. Next, we briefly highlight relevant
prior research that motivates this work.

##  II Background and Related Work

###  II-A Object Permanence in Robotics

The concept of object permanence originates from the fields of psychology and
child development, and has motivated research in developmental robotics [3].
Prior research on object permanence in robotics has focused on realizing the
concept of object permanence to enable robots to understand the effects of
their own actions and facilitate their interactions with people. For example,
object permanence has been explored to enable perspective taking and to
support situated dialog between humans and robots [4, 5]. It has also been
studied as a mechanism for self-monitoring, allowing robots to distinguish
self-generated movements from external movements and thus infer the effects of
their own actions. Such monitoring mechanisms can help robots form mental
representations of objects that were occluded by their own actions [6, 7, 8].
In contrast to previous research, the present work focuses on investigating
how object permanence may be realized by modeling audio-visual inputs and used
by a manipulator to retrieve objects that bounce out of its view.

###  II-B Trajectory Prediction

The problem of trajectory prediction has been widely investigated to enable
robot autonomy. For instance, a growing body of research focuses on modeling
and predicting pedestrian trajectories to enable mobile robots to navigate in
a safe and socially appropriate manner among people (e.g., [9, 10, 11, 12,
13]). Recent works have also investigated vehicle trajectory prediction in an
effort to develop autonomous vehicles (e.g., [14, 15]). The common property
shared by pedestrian trajectories and autonomous vehicle trajectories is that
both are goal-oriented. When the trajectories are not driven by intentional
agents but purely laws of physics, physics-based and dynamics models are
commonly used (e.g., [16, 17]). For example, neural network models have been
proposed to infer physical parameters from images, capturing the physical
interactions between a foam ball and different type of surfaces, and predict
post-bounce trajectories of the foam ball for a certain time interval [18,
19]. In our work, we focus on the prediction of bounce trajectories and end
locations of dropped objects with the goal to enhance robot autonomy in events
of unintended object drops.

###  II-C Audio-Visual Joint Reasoning

Visual and auditory signals are typically correlated in common daily
activities. This insight has inspired learning sound representations from
visual signals [20] and led to various explorations of object/scene
recognition through audio-visual joint reasoning. For example, cross-modal
models, leveraging audio-visual correspondence, have been used to map audio
and visual data to a shared space and locate the sound source in visual data
[21, 22]. With audio-guided visual attention, a multimodal residual network
has shown success in audio-visual event localization in video segments [23].
Audio-visual correspondence also allows audio signals to serve as supervision
for visual learning [24].

In the context of robotics and situated interactions, joint reasoning using
audio-visual inputs supports target localization, tracking, and hence
navigation (e.g., [25, 26, 27, 28]). For example, a deep neural network that
utilizes both vision and auditory inputs outperforms conventional methods that
rely on single modalities in object tracking [25]. Audio-visual fusion has
also been used to enhance speaker tracking [29, 30]. A simultaneous
localization and mapping (SLAM) framework that uses both audio and visual
inputs has been shown to localize both human speakers and the observer
effectively [26]. A variational expectation-maximization algorithm that takes
audio and visual signals as inputs is able to track the location of multiple
speakers in the scene [27]. Rather than focusing on tracking, this work
focuses on predicting a partially unobservable trajectory that resulted from
an object drop using the observed portion of the trajectory and a complete
impact sound recording.

Figure 2: Structure of our multimodal neural network.

##  III Problem Statement

Our goal is to enable robots to recover from inadvertent object drops. In this
work, we consider a robotic manipulator with a wrist camera and that does not
rely on external visual sensing capability for error recovery, and moreover,
we focus on erroneous situations in which objects bounce out of the robot’s
visual field. To contextualize our investigation, we experimentally let the
robot drop a wooden block and collected its partial trajectory, as observed
through its wrist camera, and the corresponding impact sound due to the
dropping. Our primary technical problem is the recovery of the full bounce
trajectory and end location of the dropped object given the partial trajectory
and impact sound.

##  IV Multimodal Network for Object Permanence

We developed a new network architecture111Our code is available at:
https://intuitivecomputing.jhu.edu/openscience.html (Fig. 2) that takes as
inputs an observed partial trajectory and a complete audio recording, and
outputs the object’s bounce trajectory and its end location with respect to
the robot.

###  IV-A Feature Extraction

####  IV-A1 Audio

Similar to the approaches proposed in SELDnet [31] and DOAnet [32], we
extracted magnitude and phase components from spectrograms of each channel of
the microphone array. The magnitude and phase components were then stacked
along the channel dimension and treated as a $M/2$ by $2\times C$ image, where
$M$ is the window length of the Fourier transformation and $C$ is the number
of channels (Fig. 2). After processing the whole audio sequence with 50%
overlap on window size, the shape of the audio data was $T\times M/2\times
2C$, where $T$ is the number of data points in the time dimension. The
parameters used for our feature extraction followed those in SELDnet [31] and
DOAnet [32].

####  IV-A2 Vision

Observed trajectories were represented as NumPy arrays with shape $(65,2)$
after being extracted from series of images using a color tracking algorithm
(detailed in Section V).

###  IV-B Audio Feature Encoder Module

As informed by SELDnet [31], we encoded audio features through convolutional
recurrent neural networks. To learn inter-channel and intra-channel features
across time, we used four 3D convolutional layers with kernel size $(k,k,c)$
to compute our audio representations. $(k,k)$ is the kernel size in the time-
frequency dimension, which aims to capture correlations within the channel
across time; $c$ is the total number of channels, which allows for finding
correlations across all available channels. 3D max pooling layers were added
between convolutional layers to reduce the dimension along the frequency axis.

The encoded audio representations were then passed to two bidirectional Gated
Recurrent Units (GRUs) (input dimension = 64, hidden dimension = 64, number of
layers = 1), seeking to further encode possible connections in the time
domain.

###  IV-C Trajectory Encoder Module

The trajectory encoder module is a three-layer multilayer perceptron (MLP)
that maps the $65\times 2$ partial observed trajectories from the robot’s
wrist camera to a higher-dimensional feature vector that matches the shape of
the output of the bidirectional GRUs ($T\times 128$).

###  IV-D Mutimodal Fusion & Correction Module

The audio representation ($T\times 128$) and the visual feature vector
($T\times 128$) are stacked into one 2D vector $(T\times 256)$, which is
subsequently fed into a two-layer MLP. The audio representation reflects the
complete trajectory, while the visual representation describes the observed,
partial trajectory. Consequently, the visual representation serves as a
mechanism to provide corrective adjustments to the location information
embedded in the audio representation. The output of the fusion and correction
module $(135\times 256)$ encodes each $(x,y)$ coordinate in a 256-dimensional
vector, with 135 being the time length for complete trajectories.

###  IV-E Trajectory Prediction Module

The trajectory prediction module is a three-layer MLP (hidden layer size:
$(64,16)$) that decodes the 256-dimensional vector representation from each
time step and outputs a corresponding $(x,y)$ location. We note that our
prediction includes the full trajectory so as to ensure robustness against
missing frames.

##  V Experiment

In this section, we describe an experiment that sought to evaluate the
effectiveness of our multimodal model in predicting the object bounce
trajectory and end location based on a partially observed trajectory and
impact sound.

###  V-A Task and Data Collection Setup

Our experimental task involved a robot picking up a red wooden cube (3 cm x 3
cm x 3 cm) and releasing it from the same height (0.3 meters above the table
surface) repeatedly. There was no systematic manipulation of object release
across trials. Observed differences in trajectories and impact sounds were
results of gripper friction and natural physics. Our data collection setup is
shown in Fig 1. We used a Kinova Gen3 robot, which has a built-in RGB-D wrist
camera (robot camera). A Microsoft Azure DK camera was used to collect impact
sounds; we note that we only used its embedded microphone array (seven
channels; 48k sampling rate). The microphone was positioned to reduce noise
from the robot motors. A RealSense D435 camera was mounted on the ceiling
(ceiling camera) to collect the ground truth videos, which included object
trajectories outside of the robot camera’s visual field. Both the robot and
ceiling cameras had a data rate of 30 frames per second.

Figure 3: (a) This histogram presents the distribution of travel distance
after the initial impact between the cube and surface. The red curve shows the
estimated Gaussian distribution. (b) This histogram depicts the distribution
of travel duration after the initial impact between the cube and surface. (c)
This hexbin plot shows the density of end locations in the 2D field as viewed
from the ceiling camera. Each cell represents the number of times the cube
ended up in that cell. The robot was placed at the top facing down. We only
visualize the end locations that were outside the initial visual field of the
robot’s wrist camera. Thus, the rectangular empty space in the middle
corresponds roughly to the size of the visual field.

###  V-B Experimental Data

####  V-B1 Data Collection and Processing

In each trial, the robot camera and the ceiling camera published RGB images
separately to different ROS topics. To synchronize all of the recording
devices (two cameras and a microphone array), two processes were run right
before the robot released the object; the first process listened to both ROS
topics and saved images in a rosbag file, whereas the second process recorded
audio data using the PyAudio package. All processes ran in parallel for three
seconds.

To extract object trajectories, we wrote a simple color tracking program using
OpenCV. In particular, we defined the range for the color red in HSV space and
masked out the contours of red objects. The center of the object contour was
recorded to represent the current location of the object.

Each camera, in practice, captured only a portion of the complete trajectory.
The robot camera captured the trajectory up to the moment when the cube
bounced outside of its view (observed trajectory). The ceiling camera was
initially blocked by the robot gripper and therefore was not able to observe
the release of the cube. However, the ceiling camera was able to keep track of
the cube after the cube was outside of the robot camera’s visual field (post
trajectory).

To map coordinate frames between the two cameras, four ChArUco markers were
used such that they were initially visible to both the ceiling camera and the
robot camera, allowing the post trajectory to be transformed to the coordinate
frame of the robot camera. After this transformation, the overlapping portion
between the observed trajectory and the post trajectory was removed from the
post trajectory.

We kept the last 65 coordinates recorded in the observed trajectory and the
first 70 coordinates recorded in the post trajectory to keep the model input
size fixed. The numbers 65 and 70 were chosen to ensure the complete capture
of each trajectory. If the length of an observed trajectory was greater than
65, the beginning of the observed trajectory could be removed without losing
information since the cube was still in the air while the first few
coordinates were recorded. If the length of a post trajectory was longer than
70, the last few coordinates could be removed without losing information since
the cube always stopped moving during the three-second movement window,
leaving the last few data points the same. We removed cases in which the cube
did not stop moving within 3 seconds; in most of these cases, the cube fell
off the table in our experimental setup. For similar reasons, if the length of
an observed trajectory or a post trajectory were less than 65 or 70, we could
prepend or postpend the first or last data point, respectively, to keep the
input size of our network fixed. After the above processing steps, a complete
trajectory was formed by concatenating an observed trajectory with its
corresponding post trajectory.

We note that our audio and trajectory data were not necessarily aligned. The
number of trajectory data points that was removed depended on factors such as
missing frames and travel direction, and was independent of the corresponding
audio data. We did not trim or pad audio data accordingly to keep two data
sequences aligned.

####  V-B2 Dataset

A total of 1404 trajectories222Data is available at:
https://doi.org/10.7281/T1/EP0W7Y was included in this experiment [33]. All
included data points contained full trajectories from the ceiling camera’s
point of view; instances where the experimental cube rolled out of the ceiling
camera’s view (most likely fell off the table) or never left the initial
visual field of the wrist camera were excluded. We also excluded the data
points where the cube bounces backward and collide with the robot.

Each data point consists of a three second long audio recording and two NumPy
2D arrays that contain the pixel location of the cube in the image frames from
the robot camera and the ceiling camera. An observed trajectory from the
robot’s wrist camera is a NumPy array with shape (65, 2); a complete
trajectory from the ceiling camera is a NumPy array with shape $(135,2)$. Each
row in the trajectories is the $(x,y)$ coordinate of the object in the image
plane. All trajectories were normalized for our proposed model, and the
initial locations of the trajectories were set as the origin using the
equation: $trajectory=trajectory[0]-trajectory$.

Fig. 3 visualizes characteristics of our dataset in terms of bounce duration,
distance traveled, and resulting end locations. We calculated the duration of
object bounces by computing the difference between the first peak and last
peak in audio wave data. We also computed the geodesic distance traveled by
the object (i.e., Euclidean distance between initial impact location and the
cube’s end location). The distributions of bounce duration and geodesic
distance can be modeled as Gaussian distributions through kernel density
estimation. The parameters for the distributions were dependent on factors
including object properties (e.g., shape, size, material) and surface
properties. Fig. 3 (c) shows the distribution of the object’s end location in
the ceiling camera’s visual field.

###  V-C Evaluation Baselines

We evaluated our proposed method against five baselines: two simple uni-modal
baselines, two additional uni-modal baselines inspired by state-of-the-art
(SOTA) methods, and one multimodal baseline that combined the two SOTA models.
We provide detailed descriptions below.

####  V-C1 B1: Linear Extrapolation of Visible Trajectory (simple vision-only
baseline)

One simple baseline for estimating the end location of a dropped object is to
project its trajectory linearly based on the last seen location ($P_{f}$)
before the object bouncing out of view, the corresponding velocity at that
location, and the remaining time to travel. We can compute the velocity ($v$)
at which the cube leaves the visual field by multiplying the coordinate
difference between the last two observed frames and the inverse of sampling
frequency (1/30). From impact sound, we can compute the duration of bouncing
($t$) after the cube leaves the visual field. Assuming that the acceleration
is constant and that velocity linearly decreases to zero, we can then estimate
the end location of the cube using the following equation:

$P_{f}+vt+\frac{1}{2}at^{2},\text{where }a=\frac{0-v}{t}$ (1)

####  V-C2 B2: Estimation of Audio Signal Delay (simple audio-only baseline)

Given the microphone arrays used in our data collection, localization
estimation through sound can be solved by leveraging signal delay across
microphone channels [34]. We used convolutional layers to extract delay
information across channels as our audio and visual data were not perfectly
synchronized due to noises and difference in sampling frequency. Specifically,
the output of our 7-D microphone array was treated as a $T$ by 7 image, where
$T$ denotes time. The image was fed through three convolutional 2D layers with
kernel sizes (16, 7), (8, 7), and (4 ,7). Each convolutional layer was
followed by a 2D batch normalization layer and a 2D pooling layer. The output
of the final convolutional layer was reshaped by a multi-layer perceptron.

####  V-C3 B3: Social GAN-lite (SOTA vision-only baseline)

Most trajectory prediction models utilize state-of-the-art recurrent neural
network architectures. This baseline was motivated by one such architecture:
the Social GAN [9], which utilizes long short-term memory (LSTM) [35] in
predicting human movement trajectories. In particular, we embedded each 2D
coordinate of the observed trajectory in a 64-dimensional vector. Then, an
encoder (implemented using LSTM, input dimension = 64, hidden dimension = 64,
number of layers = 1) took the 65 embedding vectors as inputs and outputted
its final hidden state. A decoder (implemented using LSTM, input dimension =
64, hidden dimension = 128, number of layers = 1) took the 64-dimensional
embedding of the last 2D location coordinate of the observed trajectory and
the final hidden state of the encoder as inputs and recursively predicted the
difference between consecutive locations sequentially for the 70 unseen
trajectory positions, which we normalized. Our choice of parameters was
consistent with that in Social GAN [9]. Note that the decoder predicted
recursively starting at the end location of the observed trajectory. Thus, we
needed to include the observed trajectory in our predicted complete
trajectory. This is different from our proposed multimodal method, where we
did not take the observed trajectory as exact in the prediction.

####  V-C4 B4: SELDnet-lite (SOTA audio-only baseline)

The modal used for this audio baseline was inspired by SELDnet [31] and
mirrored our multimodal network, except that it did not include the vision
input, vision encoder, and fusion components. The audio representation
produced by the bidirectional GRUs was passed to the trajectory prediction
network directly. Conceptually, this baseline addressed the sound localization
problem without vision-based corrective adjustments.

####  V-C5 B5: Combination of Social GAN-lite and SELDnet-lite

In addition to four uni-modal baselines described above, we created a
multimodal baseline that combined the Social GAN-lite and SELDnet-lite models.
Specifically, the output from the Social GAN-lite’s encoder and the output
from the SELDnet-lite’s bidirectional GRU were fused through a fully connected
linear layer and then reshaped to match target shape.

###  V-D Model Training

We split our 1404 data samples into training, validation, and test sets using
an 8:1:1 ratio. We trained all models for 30 epochs with the Adam optimizer
using a learning rate of $1e-4$ and a momentum parameter of 0.9. No early
stopping was applied, and model selection was performed using the validation
set.

###  V-E Evaluation Metrics

To compare model performance, we employed two evaluation metrics, focusing on
whether or not the robot was able to retrieve the dropped objects successfully
and on how accurate the predictions of the end locations were. Below, we
describe these metrics in detail.

####  V-E1 Task Success Rate

In practice, the robot is able to reposition itself using a simple vision-
based alignment method to retrieve the object when it is in the camera view
(Fig. 1). Therefore, we considered a prediction as successful if the dropped
object was visible to the robot’s wrist camera when the robot was located at
the end of the predicted trajectory. We defined task success rate as the
number of trials with a successful prediction divided by the total number of
trials.

$\text{Task Success Rate}=\frac{\text{Number of Successful
Trails}}{\text{Total Number of Trials}}$ (2)

####  V-E2 Target Displacement Error

Target displacement error represents the distance (in centimeters) between the
end location of a predicted bounce trajectory and the ground truth end
location.

###  V-F Results

We conducted independent two-tailed t-tests to compare the target displacement
errors from our multimodal model and the baseline models. Table I summarizes
the results of the model performance on the test dataset. Figure 4 shows
samples of trajectory predictions. Overall, our multimodal network
outperformed all baselines in all metrics, suggesting that the learned audio-
visual representation captured meaningful location information from both
modalities.

Linear extrapolation of visible trajectory (B1) made reasonable predictions
when the unobservable portion of the trajectory was mostly linear. However,
the majority of the collected trajectories were non-linear, influenced by the
object and surface properties (Fig. 4). The Social GAN-lite model (B3) also
tended to predict linear trajectories, although it appeared that the Social
GAN-lite model may have relied on two sources of information that were
implicitly contained in the partial observations for prediction: 1) the
distance between consecutive observed locations (capturing bounce speed) and
2) directional change in consecutive observed locations. When a bounce
trajectory was close to linear, this vision model worked reasonably well (Fig.
4 (c)).

Through estimation of audio signal delay, the B2 model seemed to capture non-
linearity and the general shape of a bounce trajectory. However, the
visualization of output trajectories (Fig. 4) indicated that the model failed
to capture the sequential property of input data, as predicted trajectories
were noisy and scattered. In contrast, the SELDnet-lite model (B4) seemed to
be able to capture non-linearity and enforce the sequential property of the
trajectory at the same time. We believe that the phase and magnitude
information extracted from the seven-channel microphone array carried more
nuanced information that the 2D visual trajectories did not represent, and
that the recurrent neural network component in B4 effectively maintained the
sequential property of trajectories. However, early errors (offsets)
propagated through the model and influenced later predictions. As shown in
Fig. 4 (d) and (e), the offsets accumulated and caused consequential
displacement error in the end.

B5, a multimodal baseline that combined B3 and B4, showed enhanced performance
in both task success rate and target displacement error. However, simply
combining the state-of-the-art uni-modal models was not sufficient. In B5,
audio and visual information were both processed sequentially and appeared to
interfere with, rather than complement, each other, resulting in discontinuous
trajectories (Fig. 4 (e)). Additionally, as we will discuss in the next
section, B5 did not generalize well to new bounce dynamics (Table II).

Our proposed multimodal model was able to utilize information from both audio
and vision modalities effectively. Specifically, instead of a recurrent neural
network (LSTM), a multi-layer perceptron was applied to partially observed
trajectory data and able to correctively adjust early prediction offsets and
reduced the error propagation.

In Fig. 4 (b), (d) and (e), the accumulation of early offsets led B4’s
predictions away from the ground truth. In contrast, our multimodal model was
able to utilize the observed information to reduce the errors from the audio
modality and produced more accurate predictions. However, the model might not
be able to recover from significant early prediction errors from the audio
input (potentially due to surrounding noise), as shown in Fig. 4 (e).

TABLE I: Results of the baselines and multimodal models on the test dataset.
TABLE II: Results of the finetuned baselines and multimodal model using
Dataset H and Dataset T on their respective test datasets. Figure 4: Sample
predictions from baselines (B3, B4, B5) and the multimodal model. Each column
shows prediction results from the same trial.

###  V-G Additional Exploration

While our goal in this paper is not to develop an omnipotent network that
predicts well in every possible scenario, we wanted to explore the
generalizability of our audio-visual representation. To this end, we collected
two additional small datasets: releasing the original object (cube) from a
different height (Dataset H, size 108) and releasing a different object from
the original height of 0.3 meters (Dataset T, size 115). Specifically, for
Dataset H, the cube was released 5 cm higher than the height used in the
original experiment. For Dataset T, we used a triangular wooden block to
replace the cube. The triangular block can be thought of as the original cube
cut in half along the diagonal. These datasets include novel bounce dynamics
that were not represented in the original dataset. In Dataset H, the cube made
the initial impact with more energy, which resulted in longer and more non-
linear trajectories. In Dataset T, the triangular block was lighter in weight
and non-symmetric. Therefore, the magnitude of impact audio was smaller, and
conversely noise was more evident. Trajectories in Dataset T also tended to be
non-linear.

To explore generalizability, we finetuned all baselines and our proposed model
using Dataset H and Dataset T separately, with a much smaller learning rate
($1e-5$) and 0 momentum value for 3 epochs. We performed a 5-fold cross-
validation on each dataset to account for small data size. Table II describes
the results of our exploratory examination. These results suggest that a small
set of data was sufficient in finetuning our multimodal modal to achieve
reasonable performance. Note that both the audio baseline and multimodal model
performed reasonably well, suggesting the importance of audio data in
providing indicative information.

##  VI Robot Demonstration

Fig. 1 illustrates a manipulator using our multimodal method to locate and
retrieve a dropped object. A full demonstration can be seen in our
supplementary video333https://youtu.be/Rj-ZZf3r4g8. To allow for retrieval
actions, we transformed model predictions in the robot camera frame to the
world frame. Though predictions were not perfect, as long as the object was
present in the robot’s camera view, the robot was able to use a simple vision-
based method to reposition itself for object grasping.

##  VII Discussion

Object permanence is crucial for autonomous robots to interact with objects
and operate in human environments robustly. In this project, we explore object
permanence through object localization and retrieval in the context of dropped
objects bouncing out of a robot’s visual field. We developed a multimodal
neural network that combines a partial, observed bounce trajectory and the
audio resulting from the drop impact to predict the full bounce trajectory and
the end location of a dropped object. Our results show that our approach
outperformed various baseline methods.

###  VII-A Object Permanence for Enhancing Robot Operations

The ability to estimate where dropped objects may be is important in enhancing
robot operations. As an example, our multimodal model can be used to let
robots recover from accidentally dropping objects efficiently. Rather than
relying on external sensing or heuristics-based search, our lightweight model
is able to provide reasonable estimations for object retrieval. This ability
to estimate object locations also provides robots with a sense of action
feasibility. If an object is too far to reach or in a position that the robot
cannot find a feasible motion plan to reach, the robot may instead ask for
human assistance, minimizing unnecessary grasping attempts. Overall, this
ability may be used to foster the fundamental skill of knowing when to ask for
help.

###  VII-B Limitations and Future Work

The task setup used in this work was experimentally controlled. Future work
should explore task settings in natural human environments, diverse dropping
scenarios (e.g., accidental drops during robot motions), and a variety of
experimental objects (e.g., everyday objects). Future work also needs to
explore different audio-visual representations that can encode richer nuances
in object properties and physical impact. For instance, both a higher release
point and a heavier object can result in an increase in magnitude. In addition
to different representations, future work should also consider other sensing
modalities that might contribute to the formation of object permanence
understanding. Lastly, more research is need to study object permanence in the
context of human-robot collaboration in which objects are jointly used,
manipulated, and shared.

## Acknowledgments

We would like to thank the Johns Hopkins University Institute for Assured
Autonomy for supporting this work.

## References

  * [1] J. Piaget and M. Cook, _The origins of intelligence in children_. International Universities Press New York, 1952, vol. 8, no. 5.
  * [2] A. Saxena, J. Driemeyer, and A. Y. Ng, “Robotic grasping of novel objects using vision,” _The International Journal of Robotics Research_ , vol. 27, no. 2, pp. 157–173, 2008.
  * [3] J. Luo, “Rethinking piaget for a developmental robotics of object permanence,” in _2007 IEEE 6th International Conference on Development and Learning_ , 2007, pp. 235–240.
  * [4] D. Roy, K.-Y. Hsiao, and N. Mavridis, “Mental imagery for a conversational robot,” _IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics)_ , vol. 34, no. 3, pp. 1374–1383, 2004.
  * [5] K.-y. Hsiao, N. Mavridis, and D. Roy, “Coupling perception and simulation: Steps towards conversational robotics,” in _Proceedings 2003 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS 2003)(Cat. No. 03CH37453)_ , vol. 1, 2003, pp. 928–933.
  * [6] S. Bechtle, G. Schillaci, and V. V. Hafner, “First steps towards the development of the sense of object permanence in robots,” in _2015 Joint IEEE International Conference on Development and Learning and Epigenetic Robotics (ICDL-EpiRob)_ , 2015, pp. 283–284.
  * [7] ——, “On the sense of agency and of object permanence in robots,” in _2016 Joint IEEE International Conference on Development and Learning and Epigenetic Robotics (ICDL-EpiRob)_ , 2016, pp. 166–171.
  * [8] C. Lang, G. Schillaci, and V. V. Hafner, “A deep convolutional neural network model for sense of agency and object permanence in robots,” in _2018 Joint IEEE 8th International Conference on Development and Learning and Epigenetic Robotics (ICDL-EpiRob)_ , 2018, pp. 257–262.
  * [9] A. Gupta, J. Johnson, L. Fei-Fei, S. Savarese, and A. Alahi, “Social gan: Socially acceptable trajectories with generative adversarial networks,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2018, pp. 2255–2264.
  * [10] K. D. Katyal, G. D. Hager, and C.-M. Huang, “Intent-aware pedestrian prediction for adaptive crowd navigation,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_ , 2020, pp. 3277–3283.
  * [11] A. Alahi, K. Goel, V. Ramanathan, A. Robicquet, L. Fei-Fei, and S. Savarese, “Social lstm: Human trajectory prediction in crowded spaces,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 961–971.
  * [12] Y. Tanaka, J. Kinugawa, Y. Sugahara, and K. Kosuge, “Motion planning with worker’s trajectory prediction for assembly task partner robot,” in _2012 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2012, pp. 1525–1532.
  * [13] A. Rudenko, L. Palmieri, M. Herman, K. M. Kitani, D. M. Gavrila, and K. O. Arras, “Human motion trajectory prediction: A survey,” _The International Journal of Robotics Research_ , vol. 39, no. 8, pp. 895–935, 2020\.
  * [14] N. Deo and M. M. Trivedi, “Convolutional social pooling for vehicle trajectory prediction,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops_ , 2018, pp. 1468–1476.
  * [15] A. Houenou, P. Bonnifait, V. Cherfaoui, and W. Yao, “Vehicle trajectory prediction based on motion model and maneuver recognition,” in _2013 IEEE/RSJ international conference on intelligent robots and systems_ , 2013, pp. 4363–4369.
  * [16] X. Chen, Y. Tian, Q. Huang, W. Zhang, and Z. Yu, “Dynamic model based ball trajectory prediction for a robot ping-pong player,” in _2010 IEEE International Conference on Robotics and Biomimetics_ , 2010, pp. 603–608.
  * [17] Y. Huang, D. Xu, M. Tan, and H. Su, “Trajectory prediction of spinning ball for ping-pong player robot,” in _2011 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2011, pp. 3434–3439.
  * [18] S. Purushwalkam, A. Gupta, D. M. Kaufman, and B. Russell, “Bounce and learn: Modeling scene dynamics with real-world bounces,” _arXiv preprint arXiv:1904.06827_ , 2019.
  * [19] C. Innamorati, B. Russell, D. Kaufman, and N. J. Mitra, “Neural re-simulation for generating bounces in single images,” _ICCV_ , 2019.
  * [20] Y. Aytar, C. Vondrick, and A. Torralba, “Soundnet: Learning sound representations from unlabeled video,” _Advances in neural information processing systems_ , vol. 29, pp. 892–900, 2016.
  * [21] R. Arandjelovic and A. Zisserman, “Look, listen and learn,” in _Proceedings of the IEEE International Conference on Computer Vision_ , 2017, pp. 609–617.
  * [22] ——, “Objects that sound,” in _Proceedings of the European Conference on Computer Vision (ECCV)_ , 2018, pp. 435–451.
  * [23] Y. Tian, J. Shi, B. Li, Z. Duan, and C. Xu, “Audio-visual event localization in unconstrained videos,” in _Proceedings of the European Conference on Computer Vision (ECCV)_ , 2018, pp. 247–263.
  * [24] A. Owens, J. Wu, J. H. McDermott, W. T. Freeman, and A. Torralba, “Ambient sound provides supervision for visual learning,” in _European conference on computer vision_ , 2016, pp. 801–816.
  * [25] J. Wilson and M. C. Lin, “Avot: Audio-visual object tracking of multiple objects for robotics,” 2020.
  * [26] A. Chau, K. Sekiguchi, A. A. Nugraha, K. Yoshii, and K. Funakoshi, “Audio-visual slam towards human tracking and human-robot interaction in indoor environments,” in _2019 28th IEEE International Conference on Robot and Human Interactive Communication (RO-MAN)_ , 2019, pp. 1–8.
  * [27] Y. Ban, X. Alameda-Pineda, L. Girin, and R. Horaud, “Variational bayesian inference for audio-visual tracking of multiple speakers,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2019.
  * [28] S. Pleshkova and A. Bekiarski, “Audio visual attention models in the mobile robots navigation,” in _New Approaches in Intelligent Image Analysis_ , 2016, pp. 253–294.
  * [29] S. T. Shivappa, B. D. Rao, and M. M. Trivedi, “Audio-visual fusion and tracking with multilevel iterative decoding: Framework and experimental evaluation,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 4, no. 5, pp. 882–894, 2010.
  * [30] X. Qian, A. Brutti, O. Lanz, M. Omologo, and A. Cavallaro, “Multi-speaker tracking from an audio–visual sensing device,” _IEEE Transactions on Multimedia_ , vol. 21, no. 10, pp. 2576–2588, 2019.
  * [31] S. Adavanne, A. Politis, and T. Virtanen, “Localization, detection and tracking of multiple moving sound sources with a convolutional recurrent neural network,” _arXiv preprint arXiv:1904.12769_ , 2019.
  * [32] ——, “Direction of arrival estimation for multiple sound sources using convolutional recurrent neural network,” in _2018 26th European Signal Processing Conference (EUSIPCO)_ , 2018, pp. 1462–1466.
  * [33] F. Bu and C.-M. Huang, “Dataset: Audio-visual representations of object drops,” 2020. [Online]. Available: https://doi.org/10.7281/T1/EP0W7Y
  * [34] J. Fan, Q. Luo, and D. Ma, “Localization estimation of sound source by microphones array,” _Procedia Engineering_ , vol. 7, pp. 312–317, 2010.
  * [35] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” _Neural computation_ , vol. 9, no. 8, pp. 1735–1780, 1997.

|  Fanjun Bu is a Ph.D. candidate in the Department of Computer Science at the
Cornell University. He received his B.S. degree in Applied Math and Computer
Science from the Johns Hopkins University. His research interest includes
human-robot interaction, human-computer interaction, and cognitive robotics.
---|---
|  Chien-Ming Huang is a John C. Malone Assistant Professor in the Department
of Computer Science and the Malone Center for Engineering in Healthcare at the
Johns Hopkins University. His research focuses on designing and developing
interactive systems that provide physical, social, and cognitive support to
people with diverse needs. Dr. Huang completed his postdoctoral training at
Yale University, his PhD in Computer Science at the University of
Wisconsin–Madison, and his MS in Computer Science at the Georgia Institute of
Technology. His research has received media coverage from MIT Technology
Review, Tech Insider, and Science Nation.
---|---


# Operational improvements for an algorithm to noninvasively measure the orbit
response matrix in storage rings

Volker Ziemann, Uppsala University, 75120 Uppsala, Sweden

###### Abstract

We improve the algorithm to noninvasively update the response matrix using
information from the orbit-feedback system, described in [1]. The new version
is capable of adapting to slow changes of the lattice, albeit at the expense
of limiting the accuracy.

##  1 Introduction

The orbit-response matrix relates the changes in the excitation of steering
magnets to observed position changes on the beam position monitor system [2,
3]. It is the workhorse needed to correct the beam positions and to analyze
discrepancies between an idealized model of the accelerator to the “real” one
using codes like LOCO [4, 5]. Usually, the response matrix is either derived
from a computer model or it is measured, which normally requires some
dedicated beam time. In [1] we presented a method to improve the response
matrix by exploiting correlations between the position changes and the
excitations of steering magnets caused by an orbit-feedback system. The system
runs quasi “on the side” and does not perturb the running accelerator.

Unfortunately the rate of convergence of this system, expecially with very
accurate position monitors, is very slow. Moreover, in [1] we assumed that the
model is stationary, which real accelerators, however, often are not; for
example, correcting the tunes or closing the gap of an undulator in
synchrotron light sources slightly affects the beam optics and thereby the
response matrix of the accelerator. To account for these effects, we describe
a modification of the algorithm from [1] to make it much more agile to respond
to slow and small changes of the underlying system. We base our discussion on
well-known methods from the theory of system identification described in [6].

In the next section we briefly review the model and the algorithm as well as
the improvements, before we simulate its performance in the next one. The
tradeoff between speed and accuracy of the algorithm are explored in Section 4
before we come to the conclusions.

##  2 Model

As in [1], we model the dependence of readings from $n$ position monitors by
the $\ket{x}$ on $m$ steering magnet excitations by a dynamical system

$\ket{x_{t+1}}=\ket{x_{t}}+B\ket{u_{t}}+\ket{w_{t}}\qquad\mathrm{with}\qquad\ket{u_{t}}=-K\ket{x_{t}}\
,$ (1)

where the subscript $t$ denotes a discrete time step from one iteration to the
next, $B$ is the $n\times m$ dimensional orbit response matrix, $K$ is the
$m\times n$-dimensional correction matrix of the orbit correction system, and
$\ket{w_{t}}$ describes noise in the system, characterized by the expectation
value
$\mathbf{E}\\{\ket{w_{s}}\bra{w_{t}}\\}=\sigma_{w}^{2}\delta_{st}\mathbf{1}$.
Here $\mathbf{1}$ is the $n\times n$ unit matrix, $\delta_{st}$ is the
Kronecker symbol, and $\sigma_{w}$ is the rms magnitude of the noise. We
borrow the notation with bra and ket vectors from quantum mechanics, because
keeping track of many inner and outer products becomes transparent. Here a ket
denotes a column vector and a bra denotes a row vector. Throughout this
report, the notation is consistent with [1].

Our task is now to determine an estimate $\hat{B}^{ij}_{T}$ of the matrix
elements $B^{ij}$ from recordings of all monitor readings $x_{t}^{i}$ with
$1\leq i\leq n$ and steerer excitations $u_{t}^{j}$ with $1\leq j\leq m$. Here
the subscripts denote times steps and superscripts label monitors and
steerers. We point out that the estimated matrix $\hat{B}$ depends on the time
step $T$ and typically improves as more samples are included when $T$ grows.
Note the caret to indicate that $\hat{B}$ is an estimate.

To this end we employ standard methods from the theory of system
identification [6, 7] and write Equation 1 for one monitor labeled $i$

$x^{i}_{s+1}-x^{i}_{s}=\left(\begin{array}[]{ccc}u^{1}_{s}&\dots&u^{m}_{s}\\\
\end{array}\right)\left(\begin{array}[]{c}\hat{B}^{i1}\\\ \vdots\\\
\hat{B}^{im}\end{array}\right)\ ,$ (2)

which provides us with information about row $i$ of $\hat{B}$. Stacking many
copies of this equation for successive time steps $1\leq s\leq T$ on top of
each other leads to

$\left(\begin{array}[]{c}x^{i}_{2}-x^{i}_{1}\\\ \vdots\\\
x^{i}_{T+1}-x^{i}_{T}\end{array}\right)=U_{T}\left(\begin{array}[]{c}\hat{B}^{i1}_{T}\\\
\vdots\\\ \hat{B}^{im}_{T}\end{array}\right)\qquad\mathrm{with}\qquad
U_{T}=\left(\begin{array}[]{ccc}u^{1}_{1}&\dots&u^{m}_{1}\\\ &\vdots&\\\
u^{1}_{T}&\dots&u^{m}_{T}\\\ \end{array}\right)\ .$ (3)

As $T$ increases the matrix $U_{T}$ grows by one line in each time step and we
gather more and more information about the row $i$ of $\hat{B}_{T}$ after time
step $T$. In this way Equation 3 becomes a highly overdetermined linear system
that can be solved in the least-squares sense by the pseudo inverse [8]

$\left(\begin{array}[]{c}\hat{B}^{i1}_{T}\\\ \vdots\\\
\hat{B}^{im}_{T}\end{array}\right)=\left(U_{T}^{\top}U_{T}\right)^{-1}U_{T}^{\top}\left(\begin{array}[]{c}x^{i}_{2}-x^{i}_{1}\\\
\vdots\\\ x^{i}_{T+1}-x^{i}_{T}\end{array}\right)\ .$ (4)

Of course we have to repeat the same procedure for all other rows of
$\hat{B}_{T}$ to obtain the complete estimate of the response matrix after $T$
time steps. Equation 4 describes a linear map from the vector with the
position differences on the right-hand side onto the vector with row $i$ of
$\hat{B}_{T}$. Therefore [8] $P_{T}=\left(U_{T}^{\top}U_{T}\right)^{-1}$ is
the empirical (data-driven) covariance matrix of the $\hat{B}_{T}$ after
multiplying with the error bars of the positions, which is $\sigma^{2}_{w}$.
The error bars $\sigma(\hat{B})$ of the fitted $\hat{B}_{T}$ are therefore
approximately given by the square root of the diagonal elements of
$\sigma_{w}^{2}\left(U_{T}^{\top}U_{T}\right)^{-1}$ up to a factor of order
unity.

Instead of storing and inverting $U_{T}$ after each times step, we employ the
Sherman-Morrison formula [9] to iteratively update $\hat{B}_{T}$ and the
empirical covariance matrix $P_{T}=\left(U_{T}^{\top}U_{T}\right)^{-1}$ that
appears in Equation 4. In each time step the row vector
$\bra{u_{T+1}}=(u^{1}_{T+1},\dots,u^{m}_{T+1})$ is added to the bottom of
$U_{T}$ which allows us to write
$P_{T+1}^{-1}=P_{T}^{-1}+\ket{u_{T+1}}\bra{u_{T+1}}$.

In contrast to [1], here we introduce a factor $\alpha=1-1/N_{f}$ that weighs
down the older samples [10], where $N_{f}$ is the exponential time constant
(in units of iterations) that controls this “forgetting.” We therefore write

$P_{T+1}^{-1}=\alpha P_{T}^{-1}+\ket{u_{T+1}}\bra{u_{T+1}}$ (5)

whereas in [1] we had $\alpha=1.$ Note that $0\leq\alpha\leq 1$ assigns a
weight to all rows of $U_{T}$, except the most recent one.

All derivations from [1] to invert Equation 5 are still valid, provided we
substitute $P_{T}\to P_{T}/\alpha$ in Equations 6 and 7 from [1]. After some
straightforward algebra we obtain for the updated empirical covariance matrix
$P_{T+1}$

$P_{T+1}=\frac{1}{\alpha}\left[P_{T}-\frac{P_{T}\ket{u_{T+1}}\bra{u_{T+1}}P_{T}}{\alpha+\braket{u_{T+1}}{P_{T}}{u_{T+1}}}\right]\
.$ (6)

and for the updated response matrix $\hat{B}_{T+1}$

$\hat{B}_{T+1}=\hat{B}_{T}+\frac{\left(\ket{x_{T+2}}-\ket{x_{T+1}}-\hat{B}_{T}\ket{u_{T+1}}\right)\bra{u_{T+1}}P_{T}}{\alpha+\braket{u_{T+1}}{P_{T}}{u_{T+1}}}\
.$ (7)

We refer to Appendix A and B in [1], as well as [10] for more details. We
point out the enormous advantage of the iterative procedure to update $P_{T}$
and $\hat{B}_{T}$ over repeatedly solving Equation 4 for $\hat{B}_{T+1}$. Here
we only have to store $P_{T}$ and $\hat{B}_{T}$ in memory and update them with
Equations 6 and 7 as new information represented by the monitor readings and
corresponding steerer excitations becomes available.

The following code snippet illustrates how to implement the algorithm in
MATLAB

      function [Bhatnew,Pnew,xnew]=one_iteration(Bhat,P,x,alpha)
      global sig Btilde Breal K
      u=-K*x;                                    % eq. 1, second part
      xnew=x+Breal*u+sig*randn(size(x));         % eq. 1, first part
      tmp=u’*P;                                  % <u|P
      denominv=1/(alpha+tmp*u);                  % 1/(alpha+<u|P|u>)
      Pnew=(P-tmp’*tmp*denominv)/alpha;          % eq. 6
      Bhatnew=Bhat+(xnew-x-Bhat*u)*tmp*denominv; % eq. 7

where all variables are consistently named to those used in the text. In the
next section we explore the algorithm with numerical simulations.

##  3 Simulation

We test the updated algorithm with the same model used in [1]; a FODO ring
with ten cells having phase advances of $\mu_{x}/2\pi=0.228$ in the horizontal
plane and $\mu_{y}/2\pi=0.238$ in the vertical. A position monitor and steerer
are placed at the same location as the (thin quad) focusing quadrupoles. We
use the “ideal” response matrix $\tilde{B}$ to calculate the correction matrix
$K=\left(\tilde{B}^{\top}\tilde{B}\right)^{-1}\tilde{B}^{\top}$ that appears
in Equation 1. We then randomly vary the focal lengths of all quadrupoles by 5
% to create a “real” response matrix $B$. One of the quadrupoles is varied by
an additional 5 %, which results in a second “real” response matrix $\bar{B}$
that we will use as an example to model changes to the beam optics. In all
simulations we use $\sigma_{w}=0.1\,$mm to quantify the monitor errors.

Figure 1: $|P_{T}|_{rms}$ (top) and $|b_{T}|_{rms}$ (bottom) as a function of
the iterations for $N_{f}=\infty$ (solid black), $N_{f}=200\,000$ (dashed
red), and $N_{f}=50\,000$ (dot dashed blue).

In order to assess the performance of our algorithm we introduce the
estimation error $b_{T}=\hat{B}_{T}-B$ as the difference between estimate
$\hat{B}_{T}$ and the “real” response matrix $B$ (or $\hat{B}$). The rms value
of all its matrix elements $|b_{T}|_{rms}$ can be calculated from

$|b_{T}|_{rms}=\sqrt{\frac{\mathop{\mathrm{Trace}}\left((\hat{B}_{T}-B)^{\top}(\hat{B}_{T}-B)\right)}{nm}}\
.$ (8)

In the same fashion, we introduce
$|P_{T}|_{rms}=\sqrt{\mathop{\mathrm{Trace}}\left(P_{T}^{\top}P_{T}\right)/m^{2}}$.

In a first simulation, we initialize $\hat{B}$ with the “ideal” matrix from
the computer model, while we use the real matrix $B$ to model the response of
the “real” system with equation 1. The upper panel in Figure 1 shows
$|P_{T}|_{rms}$ for one million iterations and the lower panel shows
$|b_{T}|_{rms}$. Each panel shows curves for three values of the forgetting
parameter $N_{f}$. The black curves correspond to $N_{f}=\infty$ or
$\alpha=1$, the case already covered in [1]. The red dashed curves correspond
to $N_{f}=200\,000$ and the blue dot-dashed curves to $N_{f}=50\,000$. From
the lower panel we observe that decreasing values of $N_{f}$ indeed cause
$|b_{T}|_{rms}$ to decrease more quickly, albeit at the expense of an
deteriorated asymptotic behavior. The read and blue curves no longer approach
zero, as the black one was shown to do in [1]. This observation is consistent
with the evolution of $|P_{T}|_{rms}$ shown in the upper panel. Instead of
decreasing to zero, as the black curve does, the read and blue curves
asymptotically approach finite limiting values. Considering that $P_{T}$ is
the empirical covariance matrix that describes the error bars of the
$\hat{B}_{T}$ we cannot expect them to approach $B$ arbitrarily close, as they
do with $N_{f}=\infty$.

Figure 2: $|b_{T}|_{rms}$ for the four values of $N_{f}$ indicated in the top
right corner of the respective plots. After 500 000 iterations the “real”
system matrix $B$ is replaced by $\bar{B}$ which causes the spike in the
middle of the plots.

Figure 2 shows $|b_{T}|_{rms}$ for one million iterations for six values of
$N_{f}$ where, after 500 000 iterations we replace the “real” response matrix
$B$ by $\bar{B}$, which is derived from a lattice with one quadrupole value
changed by 5 %, as mentioned before. At the same time, we also replace $B$ by
$\bar{B}$ when calculating $|b_{T}|_{rms}$, because after the new optics is in
place, we expect the algorithm to converge towards $\bar{B}$, rather than $B$.
The values of $N_{f}$ are indicated on the top left corner of the six plots.
We observe that the algorithm in all cases uses the first half of the plot to
approach $B$, already seen in Figure 1. But after the new optics is in place
after 500 000 iterations, the best current approximation $\hat{B}_{T}$ differs
significantly from new reference $\bar{B}$, which causes the spike in the
middle part of all the plots. Forgetting the old configuration and approaching
the new reference happens on the time scale given by $N_{f}$. Smaller values
lead to a faster approach. Again, at the expense of a faster approach being
paid for by an elevated asymptotic level, which we note in the title bar of
each plot to vary from around 40 mm/rad for $N_{f}=100\,000$ to more than 120
mm/rad for $N_{f}=10\,000$.

Figure 3: $|b_{T}|_{rms}$ for one million iterations while the “real” system
matrix $B$ is replaced by $\bar{B}$ after $5\times 10^{5}$ iterations. For the
first $10^{5}$ iterations we use $N_{f}=10\,000$, then use $N_{f}=200\,000$
between $10^{5}$ and $4\times 10^{5}$ iterations, where we again set
$N_{f}=10\,000$ while the change of the system matrix occurs. After $6\times
10^{5}$ iterations we set $N_{f}=200\,000$.

In order to remedy the deteriorated asymptotic level we explored whether it is
possible to temporarily adjust $N_{f}$ to allow the system to quickly react to
anticipated changes in the lattice, for example, to accommodate an undulator
gap to be closed. Figure 3 shows the configuration leading to the bottom right
plot in Figure 2 with $N_{f}=10\,000$, only here temporarily increase $N_{f}$
to 200 000 between iteration 100 000 and 400 000 and again after iteration 600
000. We clearly see that the approximation gets better during the windows with
$N_{f}=200\,000$. Once $N_{f}$ is decreased to 10 000 less information is
available and the approximation gets worse, even before the change that causes
the spike. But the system is much more agile to react and quickly adapts to
the new system, albeit with bad precision until the larger values of $N_{f}$
after iteration 600 000 improves the estimate significantly.

Based on the discussion in this section, we suggest to adapt $N_{f}$, and
thereby $\alpha$ to the anticipated running mode of the accelerator. If there
is along period of tranquility, a large value is beneficial, only to be
changed once more activity, for example, tune corrections or changes of
undulator gaps are imminent.

In the next section we will theoretically analyze the time-dependent behavior
of the system.

##  4 Convergence

Figure 4: Simulations for $|P_{T}|_{rms}$ (top) and $|b_{T}|_{rms}$ (bottom)
as well as analytic results as a function of the iterations for $N_{f}=2\times
10^{5}$ (left) and $N_{f}=5\times 10^{4}$ (right). Note how $|b_{T}|_{rms}$
initially follows and exponential form until it becomes comparable to the
noise floor, indicated by the blue dot-dashed line.

Equations 6 and 7 describe the time dependence and thus also the convergence
of the response matrix towards the “real” one. We first consider Equation 6,
because it only depends on the most recent steerer excitation through
$\ket{u_{T}}\bra{u_{T}}$, where $\ket{u_{T}}=-K\ket{x_{T}}$. Orbit correction
systems are always configured such that the correction matrix $K$ is closely
related to the inverse of the response matrix $B$, such that the largest
eigenvalue of $1-BK$ is small. By virtue of
$\ket{x_{T}}=(1-BK)\ket{x_{T-1}}+\ket{w_{T-1}}$ this implies that the most
recent position $\ket{x_{T}}$ holds no or very little memory of all previous
position and is dominated by noise $\ket{w_{T-1}}$. As a consequence we find
$Q=\mathbf{E}\left\\{\ket{u_{T}}\bra{u_{T}}\right\\}\approx\sigma_{w}^{2}KK^{\top}$
is a constant matrix. See [1] for a more detailed discussion and how to
include small additional variations of the steerers—so-called dithering—in the
analysis.

We now insert this averaged matrix into Equation 6 and find

$\hat{P}_{T+1}=\frac{1}{\alpha}\left[\hat{P}_{T}-\frac{\hat{P}_{T}Q\hat{P}_{T}}{\alpha+\mathop{\mathrm{Trace}}(Q\hat{P}_{T})}\right]\
,$ (9)

where we rewrite the expectation value in the denominator as a trace. We
placed a caret over $P_{T}$ to distinguish it from the solution of Equation 6.
In order to solve this system we now neglect the denominator, which is close
to unity and we are left with
$\hat{P}_{T+1}=(\hat{P}_{T}-\hat{P}_{T}Q\hat{P}_{T})/\alpha$ and by
subtracting $\hat{P}_{T}$ on both sides we obtain
$d\hat{P}_{T}/dT=\hat{P}_{T+1}-\hat{P}_{T}=(\hat{P}_{T}-\alpha\hat{P}_{t}-\hat{P}_{T}Q\hat{P}_{T})/\alpha$.
Here we also introduced the difference between $\hat{P}_{T}$ in two times
steps as a differential. We observe that the matrix $Q$ is by construction
symmetric and we can therefore diagonalize it, which leads to $Q=ODO^{\top}$
with $D=\mathop{\mathrm{diag}}(\lambda_{1},\dots,\lambda_{m})$ and an
orthogonal matrix $O$. The starting matrix $\hat{P}_{0}$ is the unit matrix
and always diagonal. Therefore Equation 9 can be written as $m$ independent
equations for each of the diagonal elements $p_{j,T}$ of $\hat{P}_{T}$. The
differential equation for $\hat{P}_{T}$ thus defines the corresponding one for
each of the modes with its particular eigenvalue $\lambda_{j}$

$\frac{dp_{j,T}}{dT}=\left(\frac{1-\alpha}{\alpha}\right)p_{j,T}-\left(\frac{\lambda_{j}}{\alpha}\right)p_{j,T}^{2}\
,$ (10)

which has the solution

$p_{j,T}=\frac{\beta}{\lambda_{j}/\alpha+(\beta/x_{0}-\lambda_{j}/\alpha)e^{-\beta
T}}$ (11)

with the abbreviation $\beta=(1-\alpha)/\alpha=1/(N_{f}-1)$. We clearly see
that asymptotically $p_{T}$ approaches the finite limit
$p_{j,\infty}=\alpha\beta/\lambda_{j}=1/(N_{f}\lambda_{j})$. Considering that
the error bars of reconstructed response are given in terms of $P_{T}$ we can
expect that increasing $N_{f}$ improves the approximation. Note that, apart
from $N_{f}$, only the eigenvalues $\lambda_{j}$ of
$Q=\sigma_{w}^{2}KK^{\top}$ enter. In particular, the asymptotic values are
therefore inversely proportional to $\sigma_{w}$; the algorithm works better
with noisy monitors, because it “learns from noise.” From the $p_{j,T}$ we can
reconstruct $\hat{P}_{T}$ from

$\hat{P}_{t}=O\mathop{\mathrm{diag}}(p_{1,T},\dots,p_{m,T})O^{\top}$ (12)

from which we derive $|\hat{P}_{T}|_{rms}$ in the same was as for
$|P_{T}|_{rms}$ that comes from the numerical simulation. The upper panels in
Figure 4 show them for $N_{f}=2\times 10^{5}$ on the left and $N_{f}=5\times
10^{4}$ on the right. The agreement between simulation and Equation 12 in both
cases is very good. Also the approach to finite asymptotic values is clearly
visible.

Substituting $\ket{x_{T+2}}-\ket{x_{T+1}}=B\ket{u_{T+1}}$ in Equation 7 allows
us to analyze the convergence of $\hat{B}_{T}$ towards $B$ from

$\left(\hat{B}_{T+1}-B\right)=\left(\hat{B}_{T}-B\right)-\left(\hat{B}_{T}-B\right)\frac{\ket{u_{T+1}}\bra{u_{T+1}}P_{T}}{\alpha+\braket{u_{T+1}}{P_{T}}{u_{T+1}}}\
.$ (13)

We now introduce $\hat{b}_{T}=\hat{B}_{T}-B$ to simplify writing and omit the
denominator with the trace, as before. Moreover, we replace $P_{T}$ by its
approximation $\hat{P}_{T}$ and replace $\ket{u_{T+1}}\bra{u_{T+1}}$ by its
expectation value $Q$, which brings us to
$\hat{b}_{T+1}=\hat{b}_{T}(1-Q\hat{P}_{T})$ and by turning the difference
equation into a differential equation with $\hat{b}_{T+1}-\hat{b}_{T}\approx
d\hat{b}_{T}/dT$ to $d\hat{b}_{T}/dT=-\hat{b}_{T}Q\hat{P}_{T}$. Despite
$\hat{b}_{T}$ not being simultaneously diagonal with $\hat{P}_{T}$ and $Q$, we
make the daring assumption that there are corresponding modes with eigenvalues
$\kappa_{j,T}$, such that we can write

$\frac{d\kappa_{j,T}}{\kappa_{j,T}}=-\lambda_{j}p_{j,T}dT=-\frac{\lambda_{j}\beta
dT}{\lambda_{j}/\alpha+(\beta/p_{0}-\lambda_{j}/\alpha)e^{-\beta T}}\ ,$ (14)

where we substituted $p_{j,T}$ from Equation 11. Integrating both sides, where
we note that the integral on the right-hand side is elementary, we find

$\log\left(\frac{\kappa_{j,T}}{\kappa_{j,0}}\right)=-\lambda_{j}\beta\left[\frac{\alpha
T}{\lambda_{j}}+\frac{\alpha}{\beta\lambda_{j}}\log\left(\frac{\lambda_{j}}{\alpha}+\left(\frac{\beta}{p_{0}}-\frac{\lambda_{j}}{\alpha}\right)e^{-\beta
T}\right)\right]\approx-\alpha\beta T$ (15)

where we only kept the term linear in $T$ as the leading contribution.
Replacing $\alpha\beta=1/N_{f}$ we see that the time scale on which the
difference between $\hat{B}_{T}$ and the real response matrix $B$ vanishes is
given by $e^{-T/N_{f}}$, at least in the dominant order. Since this applies to
all modes $\kappa_{j}$ we feel that the daring assumption is acceptable.

On the bottom panels in Figure 4 we show the evolution of $|b_{T}|_{rms}$
coming from a simulation as black lines and $|b_{0}|_{rms}e^{-T/N_{f}}$ as the
dot-dashed red line which shows a reasonable agreement. We also observe that
the exponential reduction only works during the initial phase until
$|b_{T}|_{rms}$ becomes comparable to the error bars that are proportional to
the magnitude of the matrix elements of the empirical covariance matrix
$|P_{T}|_{rms}$. We therefore also show the $\sqrt{|P_{T}|_{rms}}\sigma_{w}$
as an indication of these error bars. Once the exponential part of the
convergence comes to a point where $|b_{T}|_{rms}$ becomes comparable to
$\sqrt{|P_{T}|_{rms}}\sigma_{w}$ it no longer improves. The only way at this
point is to increase $N_{f}$ to reduce the asymptotic values of
$|P_{T}|_{rms}$ and thus lowers the noise floor, which allows $|b_{T}|_{rms}$
to decrease further. But this is just what Figure 3 shows.

In the simulations we used the same number of monitors and steerers ($n=m$),
but this restriction can be overcome using the methods discussed in Section
VII in [1] that we do not repeat here. Other aspects discussed there, such as
delays in the system, remain equally valid.

##  5 Conclusions

We presented an improved version of the algorithm to noninvasively measure the
orbit response matrix in storage rings. Following [10] we introduce a time
horizon $N_{f}$ after which the algorithm “forgets” old information, which
makes it much more agile to respond to new information, for example, due to a
changed “real” response matrix. We found the time constant (in numbers of
iterations) of convergence, at least initially, is given by $N_{f}$, smaller
values are favorable. On the other hand, we also found that the asymptotically
achievable accuracy is proportional to $1/N_{f}$, thus favoring large values
of $N_{f}$.

It is however, possible, to dynamically adjust $N_{f}$ to the prevailing
conditions of operation. In long periods of tranquility $N_{f}$ can be
increased, only to reduce it, once changes to the accelerator configuration
are expected.

We gratefully acknowledge fruitful discussions with Ingvar Ziemann, University
of Pennsylvania in Philadelphia.

## References

  * [1] I. Ziemann, V. Ziemann, Noninvasively improving the orbit-response matrix while continuously correcting the orbit, Physical Review Accelerators and Beams 24 (2021) 072804
  * [2] M. Minty, F. Zimmermann, Measurement and Control of Charged Particle Beams, Springer, Heidelberg, 2003.
  * [3] X. Huang, Beam-based correction and optimization for accelerators, CRC press, Boca Raton, 2020.
  * [4] J. Corbett et al., A Fast Model Calibration Procedure for Storage Rings, Proceedings of the Particle Accelerator Conference PAC93, Washington, 1993, p. 108.
  * [5] J. Safranek, Experimental determination of storage ring optics using orbit response measurements, Nuclear Instruments and Methods A 388 (1997) 27.
  * [6] G. Goodwin, R. Payne, Dynamic System Identification, Academic Press, London, 1977.
  * [7] L. Ljung, System Identification; theory for the user, 2nd ed., Prentice Hall, New Jersey, 1999.
  * [8] V. Ziemann, Regression Models and Hypothesis Testing. In: Physics and Finance. Undergraduate Lecture Notes in Physics. Springer, Cham. https://doi.org/10.1007/978-3-030-63643-2_7.
  * [9] W. Press et al., Numerical Recipes, 2nd ed., Cambridge University Press, Cambridge, 1992.
  * [10] Section 7.3 in [6].