# Generating A Crowdsourced Conversation Dataset to Combat Cybergrooming

Xinyi Zhang<EMAIL_ADDRESS>Virginia TechBlacksburgVirginiaUSA ,  Pamela J.
Wisniewski<EMAIL_ADDRESS>Vanderbilt
UniversityNashvilleTennesseeUSA ,  Jin-Hee Cho<EMAIL_ADDRESS>Virginia
TechBlacksburgVirginiaUSA ,  Lifu Huang<EMAIL_ADDRESS>Virginia
TechBlacksburgVirginiaUSA and  Sang Won Lee<EMAIL_ADDRESS>Virginia
TechBlacksburgVirginiaUSA

(2024; 20 February 2007; 12 March 2009; 5 June 2009)

###### Abstract.

##  1\. Abstract

Cybergrooming emerges as a growing threat to adolescent safety and mental
health. One way to combat cybergrooming is to leverage predictive artificial
intelligence (AI) to detect predatory behaviors in social media. However,
these methods can encounter challenges like false positives and negative
implications such as privacy concerns. Another complementary strategy involves
using generative artificial intelligence to empower adolescents by educating
them about predatory behaviors. To this end, we envision developing state-of-
the-art conversational agents to simulate the conversations between
adolescents and predators for educational purposes. Yet, one key challenge is
the lack of a dataset to train such conversational agents. In this position
paper, we present our motivation for empowering adolescents to cope with
cybergrooming. We propose to develop large-scale, authentic datasets through
an online survey targeting adolescents and parents. We discuss some initial
background behind our motivation and proposed design of the survey, such as
situating the participants in artificial cybergrooming scenarios, then
allowing participants to respond to the survey to obtain their authentic
responses. We also present several open questions related to our proposed
approach and hope to discuss them with the workshop attendees.

###### Abstract.

Cybergrooming emerges as a growing threat to adolescent safety and mental
health. One way to combat cybergrooming is to leverage predictive artificial
intelligence (AI) to detect predatory behaviors in social media. However,
these methods can encounter challenges like false positives and negative
implications such as privacy concerns. Another complementary strategy involves
using generative artificial intelligence to empower adolescents by educating
them about predatory behaviors. To this end, we envision developing state-of-
the-art conversational agents to simulate the conversations between
adolescents and predators for educational purposes. Yet, one key challenge is
the lack of a dataset to train such conversational agents. In this position
paper, we present our motivation for empowering adolescents to cope with
cybergrooming. We propose to develop large-scale, authentic datasets through
an online survey targeting adolescents and parents. We discuss some initial
background behind our motivation and proposed design of the survey, such as
situating the participants in artificial cybergrooming scenarios, then
allowing participants to respond to the survey to obtain their authentic
responses. We also present several open questions related to our proposed
approach and hope to discuss them with the workshop attendees.

Do, Not, Us, This, Code, Put, the, Correct, Terms, for, Your, Paper

††copyright: acmlicensed††journalyear: 2024††copyright:
acmlicensed††conference: The Methods for Family-Centered Design Workshop at
the CHI Conference on Human Factors in Computing Systems; May, 2024; Honolulu,
HI, USA††booktitle: The Methods for Family-Centered Design Workshop at the CHI
Conference on Human Factors in Computing Systems (CHI ’24), May, 2024,
Honolulu, HI, USA††ccs: Do Not Use This Code Generate the Correct Terms for
Your Paper††ccs: Do Not Use This Code Generate the Correct Terms for Your
Paper††ccs: Do Not Use This Code Generate the Correct Terms for Your
Paper††ccs: Do Not Use This Code Generate the Correct Terms for Your Paper

##  2\. Introduction and Background

Cybergrooming is a long-term online activity in which a predator, usually an
adult, befriends adolescents, lures them with some benefits, and eventually
solicits them for sexual encounters both online and offline (Mladenović et
al., 2021). As cybergrooming becomes more prevalent online, it has started to
make a significant impact on adolescent’s safety and health. Research has
shown that this predatory behavior negatively impacts youth mental health and
safety, leading to long-term side effects such as anxiety, depression, post-
traumatic stress, and suicide (INH, [n. d.]).

Researchers have explored different ways to combat cybergrooming. For example,
recently, with the rise of conversation agents powered by predictive
artificial intelligence (AI), researchers have also explored various
strategies to use algorithms to detect potential risks during the conversation
between predators and victims (Bours and Kulsrud, 2019). However, these
approaches primarily focus on detecting and mitigating risks and cannot be
entirely accurate due to the complex and ambiguous nature of such
conversations, and may yield misses or false positives. In addition, such
detection mechanisms come with privacy concerns for individuals, as the
algorithm needs to constantly monitor online communication. A complementary
approach is to enhance coping mechanisms and empower youth to recognize the
risks of cybergrooming by educating adolescents to be resilient.

The rapid advancement of generative AI provides new opportunities to empower
youth using conversational agents. For instance, conversational agents can
provide a simulation environment for youth to experience different kinds of
cybergrooming without exposing them to actual risks in the real world (Wang et
al., 2023). However, one challenge to developing such conversation agents is
to access datasets that provide adequate instances of authentic cybergrooming
incidents. The existing dataset that researchers (Gupta et al., 2012; Gunawan
et al., 2016; Mladenović et al., 2021; Guo et al., 2023) often rely on is the
Perverted Justice (PJ) Dataset, a public conversation corpus between predators
and victims disguised by volunteers. Yet, this dataset has some limitations.
One limitation is that adolescents in the real world may behave differently
from the ones who pose as adolescents whose goal is to lure predators. In
addition, the PJ dataset was generated between 2004 and 2016, which makes it
outdated to reflect how adolescents communicate in modern social media.
Another limitation is that the PJ dataset only has vulnerable behaviors and
without any resilient behaviors. Therefore, we do not have a dataset on how
predators respond to circumvent the resilient responses from adolescents.
These limitations make it challenging to create an authentic and comprehensive
conversation agent that can demonstrate both adolescents’ resilient and
vulnerable behaviors and how predators respond to such behaviors.

To that end, we propose leveraging both parents and adolescents in generating
data for a cybergrooming corpus. Parents naturally have a strong instinct to
protect their children from harm which could translate into proactive
behaviors. Therefore, we believe that parents can be harnessed to demonstrate
resilient behaviors. Additionally, it provides an opportunity to understand
the discrepancy between the two groups in how they perceive vulnerable and
resilient behaviors differently. In this position paper, we share how we plan
to involve both parents and adolescents to generate a dataset that can be used
to create large-scale, authentic cybergrooming conversation datasets. Such
datasets will lay a solid foundation for future research in developing
educational conversation agents or designing other educational interfaces to
empower youth to combat cybergrooming.

Figure 1. An example conversation scenario for the proposed survey

##  3\. Proposed Work: Parents/Adolescents sourcing

In our research, we propose crowdsourcing conversation responses from parents
and adolescents in simulated cybergrooming scenarios through an online survey.
We believe incorporating parents’ and adolescents’ perspectives is important
for a more comprehensive understanding of cybergrooming. Adolescents will
provide authentic responses, while parents’ responses can help us identify
effective risk-coping responses with stronger motivation to protect their
children. These data can be used to train generative NLP models capable of
providing conversations close to authentic cases. Parents and adolescents can
further provide complementary responses in many cases. For example, given the
same cybergrooming conversation scenario, parents may consider it risky, while
adolescents may treat them as safe. Such discrepancy between parents and
adolescents provides opportunities to carry out more targeted education.
Lastly, through the process of generating data, parents and adolescents will
have more awareness of how predators communicate with adolescents.

In our survey, we plan to present participants — both parents and adolescents
— conversations between victims and predators of cybergrooming. Then we will
ask participants to respond to the conversations in two different categories:
1) vulnerable behaviors that can elevate the risks of being groomed or
vulnerable, and 2) resilient behaviors that can protect themselves from the
risk. Additional follow-up questions can be used to rate how much they
perceive the risks from the given scenarios. Such survey example is shown in
Fig. 1). The team has experience involving adolescents for data collection
purposes (Razi et al., 2023, 2022), as well as co-design (Agha et al., 2023a,
b) and is well aware of potential negative implications that should be
considered when involving adolescents. For example, we do plan to recruit
parents and adolescents separately rather than jointly, as it can impact the
authenticity of their responses.

The collected data will be analyzed to identify the perceived differences
between the two groups in various components of cybergrooming: perceived
risks, resilient behaviors, and vulnerable behaviors. We will use qualitative
research methods. In addition, we plan to understand the linguistic
characteristics of four different types of data: parents’ vulnerable
behaviors, parents’ resilient behaviors, adolescents’ vulnerable behaviors,
and adolescents’ resilient behaviors Our analysis will provide insights into
developing coping skills for cybergrooming and use them to develop
conversational agents. The collected dataset will be used to develop
conversational agents.

##  4\. Open questions to discuss at the workshop

In the workshop, we would like to discuss the following questions with other
workshop attendees:

  * •

What is the effective recruitment strategy for involving adolescents and
parents in generating the dataset for cybergrooming?

  * •

Compared to traditional crowdsourcing, how should we approach our participants
(parents and adolescents) differently, considering the research context of
cybergrooming?

  * •

What ethical considerations should we take into account when collecting such a
dataset?

  * •

How should we balance the authenticity of conversation scenarios with the
protection of adolescents from exposure to sensitive content?

## References

  * (1)
  * INH ([n. d.]) [n. d.].  https://inhope.org/EN/articles/the-impact-of-online-grooming#:~:text=The%20consequences,traumatic%20stress%2C%20and%20suicidal%20thoughts.
  * Agha et al. (2023a) Zainab Agha, Karla Badillo-Urquiola, and Pamela J Wisniewski. 2023a.  ” Strike at the Root”: Co-designing Real-Time Social Media Interventions for Adolescent Online Risk Prevention.  _Proceedings of the ACM on Human-Computer Interaction_ 7, CSCW1 (2023), 1–32.
  * Agha et al. (2023b) Zainab Agha, Kelsey Miu, Sophia Piper, Jinkyung Park, and Pamela J Wisniewski. 2023b.  Co-designing user personas and risk scenarios for evaluating adolescent online safety interventions. In _Companion Publication of the 2023 Conference on Computer Supported Cooperative Work and Social Computing_. 249–253.
  * Bours and Kulsrud (2019) Patrick Bours and Halvor Kulsrud. 2019.  Detection of cyber grooming in online conversation. In _2019 IEEE International Workshop on Information Forensics and Security (WIFS)_. IEEE, 1–6.
  * Gunawan et al. (2016) Fergyanto E Gunawan, Livia Ashianti, Sevenpri Candra, and Benfano Soewito. 2016.  Detecting online child grooming conversation. In _2016 11th International Conference on Knowledge, Information and Creativity Support Systems (KICSS)_. IEEE, 1–6.
  * Guo et al. (2023) Zhen Guo, Pei Wang, Jin-Hee Cho, and Lifu Huang. 2023.  Text mining-based social-psychological vulnerability analysis of potential victims to cybergrooming: Insights and lessons learned. In _Companion Proceedings of the ACM Web Conference 2023_. 1381–1388.
  * Gupta et al. (2012) Aditi Gupta, Ponnurangam Kumaraguru, and Ashish Sureka. 2012.  Characterizing pedophile conversations on the internet using online grooming.  _arXiv preprint arXiv:1208.4324_ (2012).
  * Mladenović et al. (2021) Miljana Mladenović, Vera Ošmjanski, and Staša Vujičić Stanković. 2021.  Cyber-aggression, cyberbullying, and cyber-grooming: a survey and research challenges.  _ACM Computing Surveys (CSUR)_ 54, 1 (2021), 1–42.
  * Razi et al. (2023) Afsaneh Razi, Ashwaq AlSoubai, Seunghyun Kim, Shiza Ali, Gianluca Stringhini, Munmun De Choudhury, and Pamela J Wisniewski. 2023.  Sliding into My DMs: Detecting Uncomfortable or Unsafe Sexual Risk Experiences within Instagram Direct Messages Grounded in the Perspective of Youth.  _Proceedings of the ACM on Human-Computer Interaction_ 7, CSCW1 (2023), 1–29.
  * Razi et al. (2022) Afsaneh Razi, Ashwaq AlSoubai, Seunghyun Kim, Nurun Naher, Shiza Ali, Gianluca Stringhini, Munmun De Choudhury, and Pamela J Wisniewski. 2022.  Instagram Data Donation: A Case Study on Collecting Ecologically Valid Social Media Data for the Purpose of Adolescent Online Risk Detection. In _CHI Conference on Human Factors in Computing Systems Extended Abstracts_. 1–9.
  * Wang et al. (2023) Pei Wang, Zhen Guo, Lifu Huang, and Jin-Hee Cho. 2023.  Authentic Dialogue Generation to Improve Youth’s Awareness of Cybergrooming for Online Safety.  _Authorea Preprints_ (2023).


# On Radiation-Based Thermal Servoing:
New Models, Controls and Experiments

Luyin Hu, David Navarro-Alarcon, Andrea Cherubini and Mengying Li This work is
supported in part by the Research Grants Council (RGC) of Hong Kong under
grant number 14203917, in part by PROCORE-France/Hong Kong Joint Research
Scheme sponsored by the RGC and the Consulate General of France in Hong Kong
under grant F-PolyU503/18, in part by the Jiangsu Industrial Technology
Research Institute Collaborative Funding Scheme under grant 43-ZG9V, in part
by the Key-Area Research and Development Program of Guangdong Province 2020
under project 76 and in part by the Hong Kong PolyU under grants ZZHJ and
YBYT.Luyin Hu, David Navarro-Alarcon and Mengying Li are with The Hong Kong
Polytechnic University, Kowloon, Hong Kong. Corresponding author email:
<EMAIL_ADDRESS>Cherubini is with LIRMM, University of Montpellier CNRS,
Montpellier, France.

###### Abstract

In this paper, we introduce a new sensor-based control method that regulates
(by means of robot motions) the heat transfer between a radiative source and
an object of interest. This valuable sensorimotor capability is needed in many
industrial, dermatology and field robot applications, and it is an essential
component for creating machines with advanced thermo-motor intelligence. To
this end, we derive a geometric-thermal-motor model which describes the
relation between the robot’s active configuration and the produced dynamic
thermal response. We then use the model to guide the design of two new thermal
servoing controllers (one model-based and one adaptive), and analyze their
stability with Lyapunov theory. To validate our method, we report a detailed
experimental study with a robotic manipulator conducting autonomous thermal
servoing tasks. To the best of the authors’ knowledge, this is the first time
that temperature regulation has been formulated as a motion control problem
for robots.

###### Index Terms:

Thermoception, visual servoing, sensor-based control, robotic manipulation

Figure 1: Creatures and robots with thermo-motor intelligence: When exposed to
the sun, butterflies adjust their wings configuration to control their
temperature (a); robotic systems with thermal servoing algorithms can be used
for firefighting (b), volcano exploration (c) and industrial applications (d).

##  I INTRODUCTION

Thermal servoing is a feedback control problem that deals with the regulation
of an object’s temperature by means of motor actions of a rigid robot, which
can either manipulate the object or the heat source. It is a frontier problem
that has numerous important applications (e.g. in industrial process control,
cosmetic dermatology, fire-fighting missions, etc.) where temperature needs to
be dynamically controlled and the environment is uncertain. The quality,
performance and safety of these (otherwise open-loop) applications can be
considerably improved by incorporating thermal sensorimotor capabilities.

From a control systems perspective, the automation of this type of
temperature-critical tasks requires: (a) the computation of a geometric-
thermal-motor (GTM) model describing the relation between the robot’s motion
and the consequent thermal response, and (b) the development of a sensor-based
strategy to autonomously impose a desired heat profile onto the surface of
interest. Note that unlike other perception modalities for robot control (e.g.
vision [1], proximity [2], touch [3], audition [4] and even smell [5]),
thermoception has not been fully formalized in the literature as a bona fide
feedback signal for motion control. In the robotics community, we still lack
the framework to fully exploit it. Up to now, the overwhelming use of
thermoception in robotics has been to monitor processes but not to establish
_explicit thermal servo-loops_ [6], which are needed to accurately regulate
temperature. Our aim in this paper is precisely to develop the necessary
framework that enables the design of thermal servoing controls with radiative
heat sources.

###  I-A Related Work

Although thermal sensing is a mature technology and has a rich history in the
automation of many tasks (see e.g. [7, 8, 9, 10]), its use as a feedback
signal for robot control has not been sufficiently studied in the literature
[11], where only a few works have addressed this challenging servo-control
problem. Some representative works that deal with explicit thermal control
include: [12], where a fuzzy controller is developed to regulate the
temperature of a fuel cell actuator; [13], where the influence of temperature
in the deformation behavior of a surgical robot is investigated, and an
explicit thermal regulator is designed; [14], where a control method is
designed to maintain a constant tool temperature by adjusting the spindle
speed in a stir friction welding robot. However, in these types of methods,
temperature control is achieved by directly modulating the power of the heat-
generating components. This approach is not suitable when considering external
heat sources, e.g. wildfires [15] and sunlight [16], or when the source’s
power should not be varied, e.g. in cosmetic procedures [17].

A different strategy is to use sensor-based control, i.e. to dynamically
change the source-object geometric configuration to achieve a desired thermal
response (similar to what many organisms do [18]). This can be easily done by
rigid robots, since their basic function is motion control. Such approach
demands the development of appropriate models that can effectively capture the
system’s GTM relations. This idea has been partially demonstrated in [19],
where the optimal fixed location of multiple radiating heaters in a process is
automatically calculated to evenly imprint a desired thermal profile onto a
surface. Yet, the heater is static and the method requires exact knowledge of
all thermodynamic parameters (which are generally unknown). The proposed
approach has also the potential to be used e.g. in fire-fighting [20] or
volcano exploration robots [21] to calculate optimal trajectories that avoid
overheating or damaging the robot’s components.

The dynamic coupling between temperature and motion may seem unintuitive for
humans [22]. Yet, many organisms extensively exploit these relations. For
instance, marine animals that inhabit hydro-thermal vents manage their
energetic demands by controlling their proximity to the source [23].
Butterflies increase their body temperature by basking dorsally/laterally to
maximize sun expose [24]. Ground-dwelling insects adaptively change their body
configuration with respect to the sun-heated ground to regulate their
temperature [25]. Such advanced thermoception-based behaviors can be used to
solve many real-world problems (see Fig. 1). However, these capabilities have
not yet been fully incorporated in robot control, a discipline with good track
record of borrowing inspiration from nature [26, 27, 28], but which seems to
be lagging in this direction.

###  I-B Our Contribution

As a feasible solution to the above-mentioned issues, in this paper, we
present a rigorous formulation for robot thermal servoing with radiative
sources. The main contributions are summarized as follows:

  * •

We develop an efficient algorithm for computing in real-time the radiation-
based thermal interaction matrix which relates robot velocity and object
temperature rate.

  * •

We present a novel robot control method for automatically regulating the
temperature of grasped objects.

  * •

We report a detailed experimental study to validate the proposed theory.

To the best of the authors’ knowledge, this is the first time that temperature
regulation has been formulated in the literature as a robot servoing problem.
The proposed approach has the potential to advance the development of
multimodal motion controllers for robots.

The rest of the paper is as follows: Section II presents the mathematical
models; Section III derives the controller; Section IV reports results;
Section V gives final conclusions.

##  II MATHEMATICAL MODELING

###  II-A Notation

Throughout this manuscript, we denote all _column_ vectors by small bold
letters, e.g. $\mathbf{v}\in\mathbb{R}^{n\times 1}$, and matrices by capital
bold letters, e.g. $\mathbf{M}\in\mathbb{R}^{m\times n}$.

###  II-B Heat Transfer Model

In the following sections, we introduce basic thermodynamic concepts (we refer
the reader to [29, 30]) that are needed for developing the system’s GTM model.
To this end, consider a robot manipulator with end-effector configuration
denoted by a vector $\mathbf{x}\in\mathbb{R}^{n}$. The robot rigidly grasps
(through an adiabatic layer) a planar object whose surface temperature is to
be controlled by changing the relative pose to a heat source. The heat
transfer model is composed of three main parts (depicted in Fig. 2): (i) heat
source, (ii) heat collector (i.e. the object), and (iii) surrounding
environment. Thermophysical parameters of different parts are denoted by the
same symbol but with different subscripts. We denote the (constant)
temperature of the heat source and the (varying) temperature of the object by
${T}_{1}$ and $T_{2}$, respectively. We assume both temperatures to be
spatially uniform during the heat transfer process. The environment
temperature (assumed to be constant) is denoted by ${T}_{3}$.

######  Remark 1.

In this paper, we use the subscripts $i=1,2,3$ to denote the thermophysical
parameters of the heat source, the object and the environment, respectively (a
convention followed by many works dealing with heat transfer).

Heat transfer occurs amongst the three parts whenever ${T}_{1}$,
${T}_{2},{T}_{3}$ have different values. The direction of heat transfer is
always from a high temperature part to a low temperature part. We denote the
net energy transfer rate to the object by $Q_{2}$, where a positive value
indicates energy inflow. We introduce $q_{2}=Q_{2}/A_{2}$ to represent the
surface’s net heat flux and
$v=\mathop{}\\!\mathrm{d}T_{2}/\mathop{}\\!\mathrm{d}t$ to describe the
temporal change of the measured temperature ${T}_{2}$. According to the energy
conservation laws, these quantities satisfy the relation:

$v=\frac{1}{m_{2}c_{2}}Q_{2}$ (1)

where $m_{2}$ denotes object’s mass and $c_{2}$ denotes the material’s
specific heat. To synthesize a thermal servoing controller, it is useful to
find an expression of the following form:

$v=f(\mathbf{x},{T}_{2})$ (2)

which describes the thermal-geometric relation between the robot configuration
and the temperature rate.

###  II-C Radiation Exchange Between Planar Surfaces

In this subsection, we show how to calculate $Q_{2}$ between planar surfaces
when thermophysical properties are known. According to different mechanisms
involved in the heat transfer processes [31], the object’s net heat flux
$q_{2}$ satisfies the expression $q_{2}=q_{rad}+q_{conv}+{q}_{cond}$, for
radiative ${q}_{rad}$, convective ${q}_{conv}$ and conductive ${q}_{cond}$
fluxes. In our case of study, thermal radiation is the dominant heat transfer
mode; Note that $q_{cond}$ and $q_{conv}$ are negligible since the object is
grasped through an adiabatic layer and the source’s temperature is much higher
than those of the object and environment.

######  Assumptions 1.

We assume that the following conditions are satisfied during the task (see
Fig. 2):

  1. 1.

All surfaces have uniform temperature properties.

  2. 2.

All surfaces are gray, i.e. they are diffuse emitters with equal emittance and
absorptance.

  3. 3.

The environment/room is modeled as a black body (i.e.
$\varepsilon_{3}={\alpha}_{3}=1$, see the variables’ definition below).

The net energy transfer rate ${Q_{2}}$ has the following form:

$Q_{2}=A_{2}{q}_{rad}={A}_{2}({\alpha}_{2}{G}_{2}-{E}_{2})$ (3)

where $A_{2}$ denotes the object’s surface area, $G_{2}$ the radiative flux
incident at the surface, $\alpha_{2}\in[0,1]$ the object’s absorptance and
$E_{2}$ the heat flux emitted by a surface (i.e. emissive power) which is
approximated using the Stefan-Boltsman law [29]

$E_{2}={\varepsilon}_{2}\sigma{T}_{2}^{4}$ (4)

with $\varepsilon_{2}\in[0,1]$ the material’s emittance, and $\sigma$ is the
Stefan-Boltzmann constant.

Figure 2: Representation of the heat transfer model. A part of the object
surface is magnified to show the various heat transfer processes.

Note that for an opaque surface (i.e. with zero transmittance), its
reflectance $\rho_{i}$ and absorptance $\alpha_{i}$ satisfy
$\rho_{i}+{\alpha}_{i}=1$. Radiosity is defined as
$J_{i}={E}_{i}+{\rho}_{i}{G}_{i}$, and since for our heat source
$E_{1}\gg\rho_{1}G_{1}$, we can fairly approximate it as $J_{1}\approx E_{1}$.

######  Remark 2.

The view factor ${F}_{ij}$ represents the fraction of $J_{i}$ that is incident
on surface $j$. The view factor depends on the end-effector configuration,
i.e. $F_{ij}=F_{ij}(\mathbf{x})$. Thus, its calculation is essential for
deriving the geometric-thermal-motor model. The detailed derivation of
$F_{ij}$ for different cases and configurations is presented in later
sections. Here, we assume $F_{ij}$ in known and only focus on the derivation
of $Q_{2}$.

The radiation incident to a surface is the summation of the corresponding
portion of radiation coming from other surfaces. Thus, $G_{i}$ can be
calculated from the expression:

$A_{2}G_{2}=\sum_{j=1}^{3}F_{j2}A_{j}J_{j}=F_{12}A_{1}J_{1}+F_{22}A_{2}J_{2}+F_{32}A_{3}J_{3}$
(5)

where by using the reciprocity relation $A_{i}F_{ij}=A_{j}F_{ji}$, the
summation rule $\sum_{j=1}^{N}F_{ij}=1$, and the view factor property for
planar surfaces $F_{ii}=0$ (readers are referred to [29]), we can simplify (5)
into:

$\begin{split}{A}_{2}{G}_{2}&=F_{21}{A}_{2}J_{1}+{F}_{23}A_{2}J_{3}\\\
&=F_{21}{A}_{2}E_{1}+(1-F_{21})A_{2}E_{3}\end{split}$ (6)

Substitution of (6) into (3) yields:

$Q_{2}=A_{2}\alpha_{2}(E_{1}-E_{3})F_{21}+A_{2}\alpha_{2}E_{3}-A_{2}E_{2}$ (7)

which we substitute alongside (4) into (1) to obtain the following key
expression for the object’s temperature rate:

$v=\lambda_{1}F_{21}-\lambda_{2}{T}_{2}^{4}+\lambda_{3}$ (8)

for _constant_ scalar parameters $\lambda_{1}$, $\lambda_{2}$, and
$\lambda_{3}$ satisfying

$\lambda_{1}=\frac{A_{2}{\alpha_{2}\sigma}({\varepsilon_{1}}{T}_{1}^{4}-{T}_{3}^{4})}{m_{2}c_{2}},~{}\lambda_{2}=\frac{A_{2}{\varepsilon_{2}}\sigma}{m_{2}c_{2}},~{}\lambda_{3}=\frac{A_{2}{\alpha}_{2}\sigma{T}_{3}^{4}}{m_{2}c_{2}}$
(9)

###  II-D View Factor Analytical Definition

In this section, we provide the general expression of $F_{21}$, which we will
instantiate (in the following sections) for various configurations. To this
end, consider the elementary areas ${\mathop{}\\!\mathrm{d}}A_{1}$ and
${\mathop{}\\!\mathrm{d}}A_{2}$ on the source and object surfaces,
respectively. These areas are separated by a length $r$ that forms polar
angles $\theta_{1}$ and $\theta_{2}$ (see Fig. 3). The classical definition of
the view factor is:

$F_{21}=\frac{1}{A_{2}}\int_{A_{2}}\int_{A_{1}}\frac{\cos\theta_{2}\cos\theta_{1}}{\pi
r^{2}}\mathop{}\\!\mathrm{d}A_{2}\mathop{}\\!\mathrm{d}A_{1}$ (10)

The solution of (10) is usually complicated to derive. A variety of methods
[32, 33, 34, 35] have been proposed to calculate it. Here, we use the method
in [36], which converts the double surface integrals into double contour
integrals as follows:

$F_{21}=\frac{1}{2\pi A_{2}}\oint_{\Gamma_{1}}\oint_{\Gamma_{2}}\ln
s\mathop{}\\!\mathrm{d}\mathbf{s}_{2}\cdot\mathop{}\\!\mathrm{d}\mathbf{s}_{1},$
(11)

where $\Gamma_{i}$ denotes the contour of the $i$th surface, $\mathbf{s}_{i}$
the position vector of an arbitrary point on boundary ${\Gamma}_{i}$, and
$s=\|\mathbf{s}_{2}-\mathbf{s}_{1}\|$ the distance between two contour points.
The advantage of using this approach is its efficient computation time [37].

###  II-E Thermal Servoing with Parallel Circular Surfaces

Figure 3: Geometry of the view factor between two elementary surfaces.

In this section, we derive the thermal servoing model for two parallel source-
object surfaces. To this end, we denote the surfaces’ center and radius by
$\mathbf{c}_{i}$ and $r_{i}$, respectively. The origin of the coordinate
system $\vec{\mathbf{i}}_{1}\vec{\mathbf{j}}_{1}\vec{\mathbf{k}}_{1}$ is set
at $\mathbf{c}_{1}$, with a unit basis vector $\vec{\mathbf{k}}_{1}$ along the
normal $\vec{\mathbf{n}}_{1}$, and a unit basis vector $\vec{\mathbf{i}}_{1}$
perpendicular to the ground. We define
$\vec{\mathbf{i}}_{2}\vec{\mathbf{j}}_{2}\vec{\mathbf{k}}_{2}$ as the
translation of $\vec{\mathbf{i}}_{1}\vec{\mathbf{j}}_{1}\vec{\mathbf{k}}_{1}$,
with origin at $\mathbf{c}_{2}$. The scalars $\omega_{i}$ denote the angle
between ${\vec{\mathbf{i}}}_{i}$ and $\mathbf{s}_{i}$. We set the frames’
centers at $\mathbf{c}_{1}=[0,0,0]^{\intercal}$ and
$\mathbf{c}_{2}=[p_{1},p_{2},p_{3}]^{\intercal}$, with respect to
$\vec{\mathbf{i}}_{1}\vec{\mathbf{j}}_{1}\vec{\mathbf{k}}_{1}$. The parametric
position vectors $\mathbf{s}_{i}$ are then computed as:

$\begin{split}\mathbf{s}_{1}&=\begin{bmatrix}{r}_{1}\cos{\omega}_{1}&{r}_{1}\sin{\omega}_{1}&0\end{bmatrix}^{\intercal}\\\
\mathbf{s}_{2}&=\begin{bmatrix}{r}_{2}\cos{\omega}_{2}+p_{1}&{r}_{2}\sin{\omega}_{2}+p_{2}&p_{3}\end{bmatrix}^{\intercal}.\end{split}$
(12)

Their differential changes satisfy the following relations:

$\begin{split}&\mathop{}\\!\mathrm{d}{\mathbf{s}_{1}}=\begin{bmatrix}-{r}_{1}\sin{\omega}_{1}\mathop{}\\!\mathrm{d}{{\omega}_{1}}&{r}_{1}\cos{\omega}_{1}\mathop{}\\!\mathrm{d}{{\omega}_{1}}&0\end{bmatrix}^{\intercal}\\\
&\mathop{}\\!\mathrm{d}{{\mathbf{s}_{2}}}=\begin{bmatrix}-{r}_{2}\sin{\omega}_{2}\mathop{}\\!\mathrm{d}{{\omega}_{2}}&{r}_{2}\cos{\omega}_{2}\mathop{}\\!\mathrm{d}{{\omega}_{2}}&0\end{bmatrix}^{\intercal}\\\
&\mathop{}\\!\mathrm{d}\mathbf{s}_{1}\cdot\mathop{}\\!\mathrm{d}\mathbf{s}_{2}=r_{1}r_{2}\cos(\omega_{1}-\omega_{2})\mathop{}\\!\mathrm{d}\omega_{1}\mathop{}\\!\mathrm{d}\omega_{2}.\end{split}$
(13)

Then, the distance $s=\|\mathbf{s}_{2}-\mathbf{s}_{1}\|$ can be derived as:

$s=s(p_{1},p_{2},p_{3},\omega_{1},\omega_{2})=\left({p_{1}}^{2}+{p_{2}}^{2}+{p_{3}}^{2}+\right.\\\
2{p}_{1}(r_{2}\cos{\omega_{2}}-{r}_{1}\cos{\omega}_{1})+2{p}_{2}(r_{2}\sin{\omega_{2}}-{r}_{1}\sin{\omega}_{1})\\\
\left.+{r_{1}}^{2}+{r_{2}}^{2}-2{r_{1}}{r_{2}}\cos(\omega_{2}-\omega_{1})\right)^{\frac{1}{2}}.$
(14)

By substituting (13)–(14) into (11), $F_{21}$ can be calculated with the
expression:

$F_{21}=\int_{0}^{2\pi}\int_{2\pi}^{0}\frac{{r}_{1}{r}_{2}}{2\pi\mathrm{A}_{2}}\cos\left({\omega}_{1}-{\omega}_{2}\right)\ln{s(\mathbf{x},\omega_{1},\omega_{2})}\mathop{}\\!\mathrm{d}{\omega}_{2}\mathop{}\\!\mathrm{d}{\omega}_{1}$
(15)

where we define the end-effector position as
$\mathbf{x}=[p_{1},p_{2},p_{3}{]}^{\intercal}$. By injecting (15) into (8), we
can finally obtain the system’s thermal-geometric relation:

$v=f(\mathbf{x},T_{2})=\lambda_{1}{F}_{21}-\lambda_{2}{T_{2}}^{4}+\lambda_{3}$
(16)

where $f(\cdot)$ is the function in (2). By differentiating (16), we obtain
the following key dynamic model:

$\dot{v}=\mathbf{l}\cdot{\mathbf{u}}-4\lambda_{2}{T_{2}}^{3}v$ (17)

where $\mathbf{l}=\lambda_{1}\frac{\partial
F_{21}}{\partial\mathbf{x}}^{\intercal}$ denotes the interaction/Jacobian
matrix (vector, in this example), and
$\mathbf{u}=\dot{\mathbf{x}}\in\mathbb{R}^{n}$ is the robot’s Cartesian
velocity. The above expression is used for designing control laws for
$\dot{\mathbf{x}}$ in the following sections. By using Leibniz integral
rule[38], the interaction matrix can be expressed as:

$\mathbf{l}=\begin{bmatrix}\displaystyle\int_{0}^{2\pi}\int_{2\pi}^{0}h(p_{1}+r_{2}\cos\omega_{2}-r_{1}\cos\omega_{1})\mathop{}\\!\mathrm{d}\omega_{2}\mathop{}\\!\mathrm{d}\omega_{1}\\\
\displaystyle\int_{0}^{2\pi}\int_{2\pi}^{0}h(p_{2}+r_{2}\sin\omega_{2}-r_{1}\sin\omega_{1})\mathop{}\\!\mathrm{d}\omega_{2}\mathop{}\\!\mathrm{d}\omega_{1}\\\
\displaystyle\int_{0}^{2\pi}\int_{2\pi}^{0}h{p}_{3}\mathop{}\\!\mathrm{d}\omega_{2}\mathop{}\\!\mathrm{d}\omega_{1}\end{bmatrix}$
(18)

with the scalar $h$ defined as:

$h=\lambda_{1}\frac{r_{1}r_{2}\cos(\omega_{1}-\omega_{2})}{2\pi{A}_{2}s^{2}}\
$ (19)

Note that since it is hard to analytically compute the above double integrals,
we use a numerical method [39] to approximate $\mathbf{l}$ in real-time. See
Appendix A for details.

###  II-F Circular Surfaces in Arbitrary Configurations

Figure 4: Conceptual representation of a self-obstruction case.

In this section, we extend the parallel surfaces problem to a 6-DOF scenario,
where the end-effector configuration is now defined as
$\mathbf{x}=[p_{1},p_{2},p_{3},\theta_{x},\theta_{y},\theta_{z}{]}^{\intercal}$,
for $\theta_{i}$ as the angles around the object’s coordinate system (see Fig.
3). We denote by $\mathbf{R}$ the 3D rotation matrix corresponding to this
relative orientation. Note that in some configurations of this non-parallel
case, radiation from a source’s region cannot reach the front side of object’s
surface (hence, will not contribute to the heat inflow). We refer to this
problem as _self-obstruction_.

To model this situation, let us denote the object plane as
$D_{1}D_{2}\mathbf{c}_{2}$, for $D_{1}$ and $D_{2}$ as the intersections with
the bounded source plane. This setup is depicted in Fig. 4, where the heat
source is divided into two surfaces: the red surface composed111The symbol  ¿
$\textstyle abc$ denotes the arc that passes through the points $a$, $b$ and
$c$. of $\mathchoice{\vbox{
\hbox{\leavevmode\resizebox{34.07085pt}{0.0pt}{{\char 62\relax}}}
\nointerlineskip\hbox{\set@color$\displaystyle D_{1}D_{5}D_{2}$}}}{\vbox{
\hbox{\leavevmode\resizebox{34.07085pt}{0.0pt}{{\char 62\relax}}}
\nointerlineskip\hbox{\set@color$\textstyle D_{1}D_{5}D_{2}$}}}{\vbox{
\hbox{\leavevmode\resizebox{23.96956pt}{0.0pt}{{\char 62\relax}}}
\nointerlineskip\hbox{\set@color$\scriptstyle D_{1}D_{5}D_{2}$}}}{\vbox{
\hbox{\leavevmode\resizebox{18.83537pt}{0.0pt}{{\char 62\relax}}}
\nointerlineskip\hbox{\set@color$\scriptscriptstyle
D_{1}D_{5}D_{2}$}}}D_{2}D_{1}$ and the black surface composed of
$\mathchoice{\vbox{ \hbox{\leavevmode\resizebox{34.07085pt}{0.0pt}{{\char
62\relax}}} \nointerlineskip\hbox{\set@color$\displaystyle
D_{2}D_{6}D_{1}$}}}{\vbox{
\hbox{\leavevmode\resizebox{34.07085pt}{0.0pt}{{\char 62\relax}}}
\nointerlineskip\hbox{\set@color$\textstyle D_{2}D_{6}D_{1}$}}}{\vbox{
\hbox{\leavevmode\resizebox{23.96956pt}{0.0pt}{{\char 62\relax}}}
\nointerlineskip\hbox{\set@color$\scriptstyle D_{2}D_{6}D_{1}$}}}{\vbox{
\hbox{\leavevmode\resizebox{18.83537pt}{0.0pt}{{\char 62\relax}}}
\nointerlineskip\hbox{\set@color$\scriptscriptstyle
D_{2}D_{6}D_{1}$}}}D_{1}D_{2}$. The black surface only “sees” the object’s
backside (i.e. the robot’s gripper), thus, is omitted from the following
calculation of $F_{21}$:

$F_{21}=\frac{1}{2\pi A_{2}}\left(\oint_{\Gamma_{a}}\oint_{\Gamma_{2}}\ln
s_{a}\mathop{}\\!\mathrm{d}\mathbf{s}_{2}\mathop{}\\!\mathrm{d}\mathbf{s}_{a}+\oint_{\Gamma_{l}}\oint_{\Gamma_{2}}\ln
s_{l}\mathop{}\\!\mathrm{d}\mathbf{s}_{2}\mathop{}\\!\mathrm{d}\mathbf{s}_{l}\right)$
(20)

where $\Gamma_{a}$ denotes the arc  ¿ $\textstyle D_{1}D_{5}D_{2}$ and
$\Gamma_{l}$ the line $D_{1}D_{2}$. To derive its respective position vectors
$\mathbf{s}_{a}$ and $\mathbf{s}_{l}$, we compute the vector
$\mathbf{n}_{2}=[n_{2}^{1},n_{2}^{2},n_{2}^{3}{]}^{\intercal}$ normal to the
object plane as:

$\mathbf{n}_{2}=\mathbf{R}\begin{bmatrix}0&0&-1\end{bmatrix}^{\intercal}$ (21)

whose plane equation satisfies:

$n_{2}^{1}\left(x-p_{1}\right)+n_{2}^{2}\left(y-p_{2}\right)+n_{2}^{3}\left(z-p_{3}\right)=0$
(22)

To find the intersection with the plane, we use (21) and substitute
$x=r_{1}\cos\varphi,$ $y=r_{1}\sin\varphi$ and $z=0$ into (22), for $\varphi$
as a variable angle. Self-obstruction occurs when there exist two solutions
$\varphi_{1}$ and $\varphi_{2}$; After some algebraic operations, the arc and
line parametric vectors are obtained from:

$\mathbf{s}_{a}=\begin{bmatrix}r_{1}\cos{\varphi}\\\ r_{1}\sin{\varphi}\\\
0\end{bmatrix},~{}\mathbf{s}_{l}=\begin{bmatrix}x_{l}\\\
k_{l}(x_{l}-r_{1}\cos{\varphi_{2}})+r_{1}\sin{\varphi_{2}}\\\ 0\end{bmatrix}$
(23)

for a distance range
$x_{l}\in[r_{1}\cos{\varphi_{2}},r_{1}\cos{\varphi_{1}}]$, an angle range
$\varphi\in[\varphi_{1},\varphi_{2}]$, and a slope
$k_{l}=\frac{\sin{\varphi_{2}}-\sin{\varphi_{1}}}{\cos{\varphi_{2}}-\cos{\varphi_{1}}}$
of the line $D_{1}D_{2}$. The parametric vector on the object’s contour is
computed as:

$\mathbf{s}_{2}=\mathbf{R}\begin{bmatrix}{r}_{2}\cos{\omega}_{2}&{r}_{2}\sin{\omega}_{2}&0\end{bmatrix}^{\intercal}+\begin{bmatrix}p_{1}&p_{2}&p_{3}\end{bmatrix}^{\intercal}$
(24)

As with the parallel surface case, we compute the arc and line distances
$s_{a}=\|\mathbf{s}_{2}-\mathbf{s}_{a}\|$ and
$s_{l}=\|\mathbf{s}_{2}-\mathbf{s}_{l}\|$. The 6-DOF view factor for the self-
obstruction case is as follows:

$F_{21}=\frac{1}{2\pi\mathrm{A}_{2}}\int_{\varphi_{1}}^{\varphi_{2}}\int_{2\pi}^{0}\ln
s_{a}\mathop{}\\!\mathrm{d}{\omega}_{1}\mathop{}\\!\mathrm{d}{\varphi}\\\
+\frac{1}{2\pi\mathrm{A}_{2}}\int_{r_{1}\cos{\varphi_{2}}}^{r_{1}\cos{\varphi_{1}}}\int_{2\pi}^{0}\ln
s_{l}\mathop{}\\!\mathrm{d}{\omega}_{1}\mathop{}\\!\mathrm{d}{x_{l}.}$ (25)

With this expression, we can derive a similar GTM model
$\dot{v}=\mathbf{l}\cdot{\mathbf{u}}-4\lambda_{2}{T_{2}}^{3}v$, where
$\mathbf{u}=\dot{\mathbf{x}}\in\mathbb{R}^{6}$ and the interaction matrix is
$\mathbf{l}=\lambda_{1}\frac{\partial
F_{21}}{\partial\mathbf{x}}^{\intercal}\in\mathbb{R}^{6}$. For this 6-DOF
case, we use the following numerical differentiation method to approximate
$\mathbf{l}$:

$\mathbf{l}=\lambda_{1}\begin{bmatrix}\cfrac{F_{21}(p_{1}+\mathop{}\\!\mathrm{d}p_{1},p_{2},...,\theta_{z})-F_{21}(\mathbf{x})}{\mathop{}\\!\mathrm{d}p_{1}}\\\
\vdots\\\
\cfrac{F_{21}(p_{1},p_{2},\ldots,\theta_{z}+\mathop{}\\!\mathrm{d}\theta_{z})-F_{21}(\mathbf{x})}{\mathop{}\\!\mathrm{d}\theta_{z}}\end{bmatrix}$
(26)

Parallel programming techniques can be applied to achieve real-time
capabilities, where every element of $\mathbf{l}$ is simultaneously calculated
by an independent process.

######  Remark 3.

When the object plane and the bounded source plane do not intercept, we
classify it either as a common (non-self-obstruction) case or as a complete
self-obstruction case. This classification can be done by checking whether the
ray $\mathbf{n}_{2}$ intercepts the bounded source plane. For the common case,
the entire contour of the heat source will be used to calculate $F_{21}$. For
the complete self-obstruction case, $F_{21}=0$.

###  II-G Thermal Servoing Model with Multiple Objects

Here, we consider the case where the robot rigidly grasps $N$ objects with its
end-effector, and independently regulates the temperature of each object. For
this situation, we assume that heat exchange amongst the objects is
negligible, therefore, the derivation of the $N$-object interaction matrix
$\mathbf{L}\in\mathbb{R}^{N\times n}$ (where $n$ is the number of DOF of the
robot) is analogous to the previous sections and is simply constructed with
$N$ vectors $\mathbf{l}_{i}$ (defined for the $i$th object) as follows:

$\mathbf{L}=\begin{bmatrix}\mathbf{l}_{1}&\mathbf{l}_{2}&\cdots&\mathbf{l}_{N}\end{bmatrix}^{\intercal}.$
(27)

Figure 5: Contours of two objects approximated by a truncated Fourier series
with $5$ harmonics.

To ensure the independent control of each feedback temperature, we limit the
number of objects to be fewer than the number of robot DOF, i.e. $N\leq n$.
For this multi-object system, we construct the following structures:

$\displaystyle\boldsymbol{\tau}$
$\displaystyle=\begin{bmatrix}T_{2}^{1}&T_{2}^{2}&\cdots&T_{2}^{N}\end{bmatrix}^{\intercal}\in\mathbb{R}^{N}$
(28) $\displaystyle{\mathbf{T}}$
$\displaystyle=\operatorname{diag}\left((T_{2}^{1})^{3},(T_{2}^{2})^{3},\cdots,(T_{2}^{N})^{3}\right)\in\mathbb{R}^{N\times
N}$ (29) $\displaystyle\mathbf{v}$
$\displaystyle=\begin{bmatrix}v^{1}&v^{2}&\cdots&v^{N}\end{bmatrix}^{\intercal}\in\mathbb{R}^{N}$
(30)

The constant thermophysical parameters $\lambda_{2}$ are defined for $i$th
object as $\lambda_{2}^{i}$, and are grouped into the constant matrix:

$\boldsymbol{\Lambda}=\operatorname{diag}(\lambda_{2}^{1},\lambda_{2}^{2}\cdots,\lambda_{2}^{N})\in\mathbb{R}^{N\times
N}$ (31)

With all these terms, the geometric-thermal-motor model can be extended to a
multi-object case:

$\dot{\mathbf{v}}=\mathbf{L}{\mathbf{u}}-4{\mathbf{T}}\boldsymbol{\Lambda}\mathbf{v}$
(32)

###  II-H Irregularly Shaped Surfaces

In this section, we propose an efficient method to calculate view factors
(which, in turn, are needed for calculating interaction matrices $\mathbf{L}$)
of irregularly shaped flat surfaces. The core concept is to use truncated
Fourier series for approximating the parametric position vectors
$\mathbf{s}_{1}$ and $\mathbf{s}_{2}$. This method can be used to model both
the object and the heat source.

To see this, let us consider that the robot manipulates an object with an
irregular contour (e.g. as in Fig. 5). We can approximate the parametric
position vector of $\mathbf{s}_{2}$ with the following approach [40]:

$\mathbf{s}_{2}=\begin{bmatrix}p_{1}+\sum_{j=-F}^{F}\beta_{j}\sin(2\pi
j\phi)+b_{n}\cos(2\pi j\phi)\\\ p_{2}+\sum_{j=-F}^{N}\beta_{j}\cos(2\pi
j\phi)-b_{n}\sin(2\pi j\phi)\\\ p_{3}\end{bmatrix}$ (33)

for $\phi\in[0,1)$ as the arc-length along the contour, and $\beta_{j}$ and
$b_{j}$ constant coefficients (which can be computed as in [41]). The object’s
contour is approximated with $2F+1$ harmonic terms. The differential change
can then be computed as:

$\mathop{}\\!\mathrm{d}{\mathbf{s}}_{2}=\begin{bmatrix}\sum_{j=-F}^{F}2\pi
j\left(\beta_{j}\cos(2\pi j\phi)-b_{j}\sin(2\pi j\phi)\right)\\\
\sum_{j=-F}^{F}-2\pi j\left(\beta_{j}\sin(2\pi j\phi)+b_{j}\cos(2\pi
j\phi)\right)\\\ 0\end{bmatrix}\mathop{}\\!\mathrm{d}{\phi}$ (34)

The thermal interaction matrix $\mathbf{L}$ can now be easily derived by
following procedures similar to the ones presented in the previous sections.
In Appendix B, we present a numerical study comparing the performance of this
method with the standard discrete surface integral [42].

##  III CONTROLLER DESIGN

######  Problem statement.

Given a _constant_ temperature reference vector
$\boldsymbol{\tau}^{*}=[{T}^{*1},\ldots,{T}^{*N}{]}^{\intercal}\in\mathbb{R}^{N}$,
design a velocity-based motion controller $\mathbf{u}$ that asymptotically
minimizes the feedback error
$\Delta\boldsymbol{\tau}=\boldsymbol{\tau}-\boldsymbol{\tau}^{*}$ for all $N$
objects.

To solve this problem, in this section we propose two methods: a model-based
controller (conceptually depicted in Fig. 6) and an adaptive controller (shown
in Fig. 7). At the end of this section, we discuss target feasibility.

###  III-A Model-Based Controller

Figure 6: Schematic representation of the model-based controller.

To automatically regulate the feedback temperature of $N$ objects, we design
the following velocity control input:

$\mathbf{u}=\mathbf{L}^{+}(-D\mathbf{v}-K\Delta\boldsymbol{\tau}+4{\mathbf{T}}\boldsymbol{\Lambda}\mathbf{v})$
(35)

where
$\mathbf{L}^{+}=\mathbf{L}^{\intercal}\left(\mathbf{L}\mathbf{L}^{\intercal}\right)^{-1}$
is the Moore-Penrose pseudoinverse of $\mathbf{L}$ [43], and $D>0$ and $K>0$
are control gains.

######  Proposition 1.

Consider that thermodynamic parameters in (32) are accurately known. For this
situation, the control input (35) enforces a stable closed-loop system which
asymptotically minimizes $\|\Delta\boldsymbol{\tau}\|$.

###### Proof.

Substituting (35) into the nonlinear dynamic system (32), yields the following
closed-loop system:

$\dot{\mathbf{v}}=-D\mathbf{v}-K\Delta\boldsymbol{\tau}.$ (36)

Consider the quadratic Lyapunov function

$\mathcal{Q}(\mathbf{v},\Delta\boldsymbol{\tau})=\frac{1}{2}\|\mathbf{v}\|^{2}+\frac{1}{2}K\|\Delta\boldsymbol{\tau}\|^{2}$
(37)

whose time derivative along trajectories of (36) yields

$\begin{split}\dot{\mathcal{Q}}(\mathbf{v},\Delta\boldsymbol{\tau})&=\mathbf{v}^{\intercal}\dot{\mathbf{v}}+K\Delta\boldsymbol{\tau}^{\intercal}{\mathbf{v}}\\\
&=-D\|\mathbf{v}\|^{2}\end{split}$ (38)

which shows that the energy function is non-increasing, i.e.
$\dot{\mathcal{Q}}\leq 0$, thus, the closed-loop system is stable. By applying
the Krasovskii-LaSalle principle [44], the asymptotic minimization of
$\|\Delta\boldsymbol{\tau}\|$ can be proved. ∎

######  Remark 4.

In our proposed method, the terms $\mathbf{T}$ and $\boldsymbol{\tau}$ in the
controller (35) can be directly obtained from real-time sensor measurements.
Yet, to implement the variable $\mathbf{v}$, we use a rate estimation
algorithm based on polynomial fitting with sliding windows [45] (see Appendix
C for details).

###  III-B Adaptive Controller

In the above model-based controller, we assume that the object’s
thermophysical properties are exactly known. However, due to the differences
in material and surface conditions, it is hard to determine the true values.
In this subsection, we propose an adaptive controller which updates the
unknown parameters within the servo-loop. To this end, we start by introducing
the unknown parameters $a_{1}=\frac{1}{\lambda_{1}}$ and
${a}_{2}=\frac{\lambda_{2}}{\lambda_{1}}$, which are well-defined since
$\lambda_{1}>0$; We use the superscripts ${a}_{1}^{i}$ and ${a}_{2}^{i}$ to
distinguish them between different objects. With these parameters, we
construct the following _constant_ vector $\mathbf{a}_{1,2}\in\mathbb{R}^{N}$
and matrix $\mathbf{A}_{1,2}\in\mathbb{R}^{N\times N}$ structures:

$\displaystyle\mathbf{a}_{1}$
$\displaystyle=\begin{bmatrix}a_{1}^{1}&\cdots&a_{1}^{N}\end{bmatrix}^{\intercal},\qquad\mathbf{A}_{1}=\operatorname{diag}(\mathbf{a}_{1})>0,$
$\displaystyle\mathbf{a}_{2}$
$\displaystyle=\begin{bmatrix}a_{2}^{1}&\cdots&a_{2}^{N}\end{bmatrix}^{\intercal},\qquad\mathbf{A}_{2}=\operatorname{diag}(\mathbf{a}_{2})>0.$
(39)

By applying the dynamic expression (17) to the multi-object case and dividing
it by $\lambda_{1}^{i}$ for each $i$th object, we obtain:

$\mathbf{A}_{1}\dot{\mathbf{v}}+4\mathbf{T}{\mathbf{A}}_{2}\mathbf{v}=\mathbf{J}\mathbf{u}$
(40)

for an new interaction matrix

$\mathbf{J}=\begin{bmatrix}\mathbf{l}_{1}/\lambda_{1}^{1}&\mathbf{l}_{2}/\lambda_{1}^{2}&\cdots&\mathbf{l}_{N}/\lambda_{1}^{N}\end{bmatrix}^{\intercal}\in\mathbb{R}^{N\times
6}$ (41)

which is independent from the unknown thermophysical parameters222An analogous
result to the depth-independent interaction matrix in [46] and is entirely
computed with the $N$ gradients of the view factors. To design the adaptive
controller, it is useful to introduce the combined thermal error vector
$\boldsymbol{\zeta}=[\zeta^{1},\ldots,\zeta^{N}{]}^{\intercal}\in\mathbb{R}^{N}$
defined as:

$\boldsymbol{\zeta}=\Delta\dot{\boldsymbol{\tau}}+\mu\Delta\boldsymbol{\tau}=\mathbf{v}+\mu\Delta\boldsymbol{\tau}$
(42)

for $\mu>0$ as an arbitrary feedback gain. To control the $N$ object
temperatures, we design the following velocity input:

$\mathbf{u}=\mathbf{J}^{+}(-\mu\widehat{\mathbf{A}}_{1}\mathbf{v}-K\boldsymbol{\zeta}+4\mathbf{T}\widehat{\mathbf{A}}_{2}\mathbf{v})$
(43)

where the elements of the adaptive diagonal matrices
$\widehat{\mathbf{A}}_{i}=\operatorname{diag}(\widehat{\mathbf{a}}_{i})\in\mathbb{R}^{N\times
N}$ are computed with the update rules:

$\displaystyle\dot{\widehat{\mathbf{a}}}_{1}$
$\displaystyle=\gamma_{1}\mu\begin{bmatrix}v^{1}\zeta^{1}&\ldots&v^{N}\zeta^{N}\end{bmatrix}^{{\intercal}}\in\mathbb{R}^{N}$
(44) $\displaystyle\dot{\widehat{\mathbf{a}}}_{2}$
$\displaystyle=-4\gamma_{2}\begin{bmatrix}v^{1}\zeta^{1}(T_{2}^{1})^{3}&\ldots&v^{N}\zeta^{N}(T_{2}^{N})^{3}\end{bmatrix}^{{\intercal}}\in\mathbb{R}^{N}$
(45)

where the positive scalars $\gamma_{i}>0$ are used for tuning algorithm’s
learning rate.

Figure 7: Schematic representation of the adaptive controller.

######  Proposition 2.

The adaptive controller (43) with update rules (44)–(45) guarantees a bounded
estimation of the unknown parameters $\mathbf{a}_{1}$ and $\mathbf{a}_{2}$,
and the asymptotic minimization of the thermal error
$\|\Delta\boldsymbol{\tau}\|$.

###### Proof.

Substitution of (43) into (40) yields:

$\mathbf{A}_{1}\dot{\mathbf{v}}+4\mathbf{T}{\mathbf{A}}_{2}\mathbf{v}=-\mu\widehat{\mathbf{A}}_{1}\mathbf{v}-K\boldsymbol{\zeta}+4\mathbf{T}\widehat{\mathbf{A}}_{2}\mathbf{v}$
(46)

By adding $\mu\mathbf{A}_{1}\mathbf{v}$ to both sides of (46), noting that
$\dot{\boldsymbol{\zeta}}=\dot{\mathbf{v}}+\mu\mathbf{v}$, and performing some
algebraic operations we can obtain:

$\mathbf{A}_{1}\dot{\boldsymbol{\zeta}}+K\boldsymbol{\zeta}=-\mu\widetilde{\mathbf{A}}_{1}\mathbf{v}+4\mathbf{T}\widetilde{\mathbf{A}}_{2}\mathbf{v}.$
(47)

for error matrices
$\widetilde{\mathbf{A}}_{i}=\widehat{\mathbf{A}}_{i}-\mathbf{A}_{i}=\operatorname{diag}(\widetilde{\mathbf{a}}_{i})$,
with error vectors
$\widetilde{\mathbf{a}}_{i}=\widehat{\mathbf{a}}_{i}-\mathbf{a}_{i}$. To
analyze the stability of the closed-loop dynamical system (44)–(45) and (47),
we introduce the following Lyapunov function:

$\mathcal{H}(\boldsymbol{\zeta},\widetilde{\mathbf{a}}_{1},\widetilde{\mathbf{a}}_{2})=\frac{1}{2}\boldsymbol{\zeta}^{\intercal}\mathbf{A}_{1}\boldsymbol{\zeta}+\frac{1}{2\gamma_{1}}\|{\widetilde{\mathbf{a}}}_{1}\|^{2}+\frac{1}{2\gamma_{2}}\|{\widetilde{\mathbf{a}}}_{2}\|^{2}$
(48)

whose time derivative along (44)–(45) and (47) yields

$\begin{split}\dot{\mathcal{H}}(\boldsymbol{\zeta},\widetilde{\mathbf{a}}_{1},\widetilde{\mathbf{a}}_{2})&=\boldsymbol{\zeta}^{{\intercal}}\mathbf{A}_{1}\dot{\boldsymbol{\zeta}}+\frac{1}{\gamma_{1}}\dot{\widehat{\mathbf{a}}}_{1}^{{\intercal}}\widetilde{\mathbf{a}}_{1}+\frac{1}{\gamma_{2}}\dot{\widehat{\mathbf{a}}}_{2}^{{\intercal}}\widetilde{\mathbf{a}}_{2}\\\
&=-K\|\boldsymbol{\zeta}\|^{2}\end{split}$ (49)

which shows that the energy function is non-increasing, i.e.
$\dot{\mathcal{H}}\leq 0$, thus, the parameter estimation errors
$\widetilde{\mathbf{a}}_{i}$ are bounded. Asymptotic stability of
$\Delta\boldsymbol{\tau}$ directly follows by applying the Krasovskii-LaSalle
principle [44]. ∎

###  III-C Target Feasibility

In previous sections, we proved that $\|\Delta\boldsymbol{\tau}\|$ can be
asymptotically minimized by two automatic controllers. However, it is not
guaranteed that such error can be enforced to zero. Failure cases are caused
by the choice of unfeasible target temperatures: Intuitively, if targets are
set to too high or too low, they might be physically unachievable; In
addition, for objects fixed to the same end-effector, the difference range
between their target temperatures is constrained by the fixed spatial
relationship between the objects. In this section, we analyze two necessary
but not sufficient conditions to ensure the feasibility of the targets. Failed
experimental results are analyzed accordingly in Section IV-F.

Consider a simple case with two objects, object 1 and object 2, fixed to the
end-effector (the extension to $N$ object is straightforward). For one of the
objects, recall the thermal-geometric relation $\eqref{eq:final temperature
rate}$ and rewrite it as follows:

$v=-\lambda_{2}{T_{2}}^{4}+\lambda_{1}{F}_{21}(\mathbf{x}_{o})+\lambda_{3}$
(50)

where $\mathbf{x}_{o}$ denotes an object configuration. Let us assume there
exists a temperature $T_{2}=T_{v0}$ that makes the rate $v=0$. As temperature
is always non-negative, $T_{v0}$ can be solved as:

$T_{v0}(F_{21})=\left((\lambda_{1}{F}_{21}+\lambda_{3})/\lambda_{2}\right)^{\frac{1}{4}}$
(51)

Since the parameters $\lambda_{i}>0$ are all positive and $F_{21}\in[0,1)$,
$T_{v0}$ always exists. $T_{v0}$ represents the steady state temperature at
$\mathbf{x}_{o}$. Note that $T_{v0}$ is a function of $F_{21}$ and that
$\partial T_{v0}/\partial F_{21}>0$ is always positive. Thus, the minimum
value of $T_{v0}$ is determined when $F_{21}=0$ as:

$\operatorname{min}(T_{v0})=(\lambda_{3}/\lambda_{2})^{\frac{1}{4}}=(\alpha_{2}T_{3}^{4}/\varepsilon_{2})^{\frac{1}{4}}$
(52)

According to Kirchhoff’s law of thermal radiation [29], at thermodynamic
equilibrium, $\alpha_{2}=\varepsilon_{2}$. Thus, the minimum is:

$\operatorname{min}(T_{v0})=T_{3}$ (53)

When $F_{21}(\mathbf{x}_{o})\to 1$, the maximum value of $T_{v0}$ approaches:

$\operatorname{max}(T_{v0})\to((\lambda_{1}+\lambda_{3})/\lambda_{2})^{\frac{1}{4}}=(\alpha_{2}\varepsilon_{1}T_{1}^{4}/\varepsilon_{2})^{\frac{1}{4}}=\varepsilon_{1}^{\frac{1}{4}}T_{1}$
(54)

From (53)–(54), we derive the _first_ boundary value condition:

$T^{1*},T^{2*}\in[T_{3},\varepsilon_{1}^{\frac{1}{4}}T_{1})$ (55)

Now we discuss the limitation of the difference between target temperatures
$|\delta T^{*}|=|T^{1*}-T^{2*}|$. We denote the configuration of the object 1
and object 2 by $\mathbf{x}_{o1}=\mathbf{x}+\Delta\mathbf{x}_{1}$ and
$\mathbf{x}_{o2}=\mathbf{x}+\Delta\mathbf{x}_{2}$, respectively, where
$\Delta\mathbf{x}_{1},\Delta\mathbf{x}_{2}$ are constant displacement vectors
determined by the arrangement of objects. Their corresponding view factor are
denoted by $F_{21}(\mathbf{x}_{o1})$ and $F_{21}(\mathbf{x}_{o2})$, and its
steady-state temperatures by $T_{v0}^{1}$ and $T_{v0}^{2}$. According to
$\eqref{eq:steady state temperature}$, $|\Delta
T_{v0}|=|T_{v0}^{1}-T_{v0}^{2}|$ can be expressed as:

$|\Delta
T_{v0}(\mathbf{x}_{o1},\mathbf{x}_{o2})|=|\chi_{1}(F_{21}(\mathbf{x}_{o1}))-\chi_{2}(F_{21}(\mathbf{x}_{o2}))|$
(56)

with functions $\chi_{1}(F_{21})$ and $\chi_{2}(F_{21})$ defined as:

$\chi_{1}(F_{21})=\left(\frac{\lambda_{1}^{1}{F}_{21}+\lambda_{3}^{1}}{\lambda_{2}^{1}}\right)^{\frac{1}{4}},~{}\chi_{2}(F_{21})=\left(\frac{\lambda_{1}^{2}{F}_{21}+\lambda_{3}^{2}}{\lambda_{2}^{2}}\right)^{\frac{1}{4}}$
(57)

where $\lambda_{i}^{1}$ and $\lambda_{i}^{2}$ are the thermophysical
parameters of the two objects. Note that for the continuous function $\Delta
T_{v0}(\mathbf{x}+\Delta\mathbf{x}_{1},\mathbf{x}+\Delta\mathbf{x}_{2})$,
where $\mathbf{x}\in\mathbb{W}$ for $\mathbb{W}$ as the bounded workspace and
$\Delta\mathbf{x}_{j}$ as constant vectors, there must exist a minimum value
$\operatorname{min}(\Delta T_{v0})=\Delta T_{v0}(\mathbf{x}^{min})$ and a
maximum value $\operatorname{max}(\Delta T_{v0})=\Delta
T_{v0}(\mathbf{x}^{max})$ which encompass all possible values of $\Delta
T_{v0}$, where $\mathbf{x}^{min}$ and $\mathbf{x}^{max}$ are the end-effector
configurations corresponding to the two extreme cases. The _second_ condition
for feasible target temperatures is:

$\delta T^{*}\in[\operatorname{min}(\Delta T_{v0}),\operatorname{max}(\Delta
T_{v0})]$ (58)

A numerical (geometric) interpretation of $\mathbf{x}^{min}$ and
$\mathbf{x}^{max}$ will be discussed in Section IV-F.

##  IV RESULTS

###  IV-A Experimental Setup

Figure 8: Experimental setup. Figure 9: Different object used to test the two
controllers. Figure 10: Calibration using ArUco markers before experiments.

We conducted a series of experiments on a 4-DOF robot (3 translations and 1
rotation) to evaluate the proposed method. Fig. 8 shows the robot, whose end-
effector is replaced by a 3D printed connector fixed to an aluminum holder.
The objects are attached to the holder through an adiabatic layer to minimize
heat conduction. We prepared three different kinds of objects for temperature
control experiments (see Fig. 9): An aluminum circular sheet with
$1.5\text{\,}\mathrm{cm}$ radius and $3\text{\,}\mathrm{mm}$ thickness; A
bunny-shaped object with $1\text{\,}\mathrm{mm}$ thickness, 3D printed using
polylactic acid (PLA) material with $30\%$ infill density; A hand-shaped sheet
with $1\text{\,}\mathrm{mm}$ thickness, also 3D printed using PLA but with
$50\%$ infill density. We approximate the aluminum sheet’s and the heat
source’s thermophysical properties via standard tables [29]. The object’s
emittance, absorptance, specific heat, and density are 0.04, 0.04,
$903\text{\,}\mathrm{J}\text{$\cdot$}{\mathrm{K}}^{-1}\text{$\cdot$}{\mathrm{kg}}^{-1}$,
and $2702\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$, respectively.
The source’s emittance and absorptance are estimated as 0.25 and 0.25.

For the two 3D printed objects, different infill densities, colors, and
uncertain surface conditions make their thermophysical properties hard to be
estimated. Thus, we only consider the aluminum sheet for the experiments with
the model-based controller (adaptive control is used for the other objects). A
radiative heating platform with adjustable temperature output is used as the
heat source. The (indoor) environment temperature is assumed to be constant at
$23\text{\,}\mathrm{\SIUnitSymbolCelsius}$.

To obtain the feedback temperatures, we attach a PT100 platinum thermistor
with $0.3\text{\,}\mathrm{\SIUnitSymbolCelsius}$ accuracy and
$0.1\text{\,}\mathrm{\SIUnitSymbolCelsius}$ precision to each object. The raw
data obtained by thermistors is processed by a current-temperature
transformation module and sent to a Linux-based control computer as the
feedback signal. The motion command is calculated by the computer program and
sent to the robot under a position-stepping mode. At the beginning of the
experiments, we use an RGB camera and three ArUco markers [47] to calibrate
the configuration between the heat source and the end-effector (see Fig. 10).

###  IV-B Experiments with the Model-Based Controller

Figure 11: Snapshots of a representative thermal servoing experiment: (a)
Initial position, (b)–(c) transient motion, and (d) steady-state
configuration. Figure 12: Evolution $\Delta T$ of the temperature error using
one aluminum object with the model-based controller. Figure 13: Evolution of
the temperature errors ($\Delta T^{1},\Delta T^{2}$) using two aluminum
objects with the model-based controller.

We conduct a series of thermal servoing experiments to evaluate our proposed
control methodology (see Fig. 11 for a representative experiment). Here, we
first evaluate the performance of the model-based controller with aluminum
objects (whose properties are approximately known); The experiments are
conducted with a source’s temperature of
$200\text{\,}\mathrm{\SIUnitSymbolCelsius}$.

We denote the scalar temperature error by $\Delta T=T_{2}-T^{*}$. By using the
controller (35), we enforce a closed-loop heat transfer system that resembles
a mass-spring-damper system. Therefore, the values of the stiffness/damping-
like gains $K$ and $D$ can be used to specify the system’s performance. Fig.
12 (a) demonstrates the effect of the gain $K$ on the thermal response. For
that, we set $D=0.2$ and $T^{*}=$50\text{\,}\mathrm{\SIUnitSymbolCelsius}$$
and conduct three experiments with different $K$ values. These show that when
$K=0.005$ (red curve), the error $\Delta T$ asymptotically decreases to zero
with a relatively slow speed; When $K=0.05$ (blue curve), $\Delta T$ decreases
faster and a small overshoot occurs; When $K=0.5$ (black curve), $\Delta T$
oscillates near zero with an approximate
$3\text{\,}\mathrm{\SIUnitSymbolCelsius}$ amplitude. This results shows how
the closed-loop system varies from over-damped to under-damped. Thus, the
gains should be specified according to the desired thermodynamic performance.
We further conducted experiments with the same gains ($K=0.05$, $D=0.2$) but
with different targets $T^{*}$ and found a consistent response (see Fig. 12
(b)).

Model-based experiments were also conducted to independently regulate the
temperatures of two aluminum objects, shown in Fig. 9 (b). We designed 4
experiments with different targets $\boldsymbol{\tau}^{*}$ (measured in
$\text{\,}\mathrm{\SIUnitSymbolCelsius}$). Fig. 13 depicts the minimization of
the thermal errors for these 4 experiments, with target temperatures defined
as $\boldsymbol{\tau}^{*}=[{50},{40}{]}^{{\intercal}}$,
$\boldsymbol{\tau}^{*}=[{60},{40}{]}^{{\intercal}}$,
$\boldsymbol{\tau}^{*}=[{50},{50}{]}^{{\intercal}}$, and
$\boldsymbol{\tau}^{*}=[{80},{40}{]}^{{\intercal}}$ in (a), (b), (c), and (d),
respectively. For the first three experiments where the differences between
the target temperatures $|T^{*1}-T^{*2}|$ are small (or null), the thermal
error $\|\Delta\boldsymbol{\tau}\|$ can be asymptotically minimized to zero.
However, when $|T^{*1}-T^{*2}|$ is large, as in Fig. 13 (d), the two
temperatures cannot be accurately controlled. This failure case can be
explained by the second condition for feasible targets discussed in Section
III-C.

###  IV-C Experiments with the Adaptive Controller

We designed a series of experiments to evaluate the performance of the
proposed adaptive controller. For that, we consider with three different
objects (see Fig. 9 (c)) with unknown thermophysical properties and irregular
shapes. To compute the interaction matrix, we use truncated Fourier series
with 5 harmonics terms; This approach provides a fast calculation time with a
“good enough” shape approximation. The controller’s gains are set to
$\mu=0.05$ and $K=0.15$. To initialize the parameters
$\widehat{\mathbf{a}}_{i}(0)$ at the time instance $t=0$, we use (for the
“hand” and “bunny” objects) the constant values calculated for the aluminum
object in the previous model-based controller, i.e.
$\widehat{\mathbf{a}}_{i}(0)=\mathbf{a}_{i}$; For the circular object, we
simply initialize $\widehat{\mathbf{a}}_{i}(0)$ with random values.

In this study, we report eight temperature control experiments with different
targets, objects and source conditions. Figure 14 shows the evolution of the
individual thermal errors $\Delta T^{i}$. For ease of presentation, we name
these eight experiments as exp 1, …, exp 8, and denote the corresponding
target temperature for each experiment by
$\boldsymbol{\tau}^{*1},\ldots,\boldsymbol{\tau}^{*8}$. In exp 1, we set the
three target temperatures to the same value. In exp 2 – exp 4, we only set two
targets to the same value. In exp 5 and exp 6, we set all targets to different
values, with a non-uniform thermal separated in exp 5 and a uniform one in exp
6. In exp 7 and exp 8, all targets are set to the same value, but with
different heat source conditions. The source temperature $T_{1}$ is set to
$200\text{\,}\mathrm{\SIUnitSymbolCelsius}$ in exp 1 – exp 6, to
$300\text{\,}\mathrm{\SIUnitSymbolCelsius}$ in exp 7, and varies from
$200\text{\,}\mathrm{\SIUnitSymbolCelsius}$ to
$300\text{\,}\mathrm{\SIUnitSymbolCelsius}$ in exp 8.

In all these experiments with all these different conditions, the magnitude of
the temperature error $\|\boldsymbol{\tau}\|$ asymptotically decreases to
zero. Yet, failure control experiments do happen and are reported and
discussed in Section IV-F). The results experimentally confirm that (for
_feasible_ target temperatures) the adaptive method is able to independently
regulate temperatures of various objects with different shapes and materials,
without exact knowledge of their thermophysical properties or the
source’s/environment’s temperatures.

Fig. 15 depicts the performed object trajectories during the experiments in
Fig. 14. The boundary of the circular heat source is depicted as a black
circle (and ellipse). The color of a trajectory point represents the feedback
temperature at that position; Variation from blue to red corresponds to a
change from “low” to “high”. For clarity, we depict two sets of trajectory
visualizations from different viewing angles: For Fig. 15
$(a_{1}),(b_{1}),\ldots,(h_{1})$, the trajectories are viewed in
$-\vec{\mathbf{k}}_{1}$ direction; For Fig. 15
$(a_{2}),(b_{2}),\ldots,(h_{2})$, the trajectories are viewed in
$\vec{\mathbf{i}}_{1}$ direction.

From these trajectory visualizations, we can see that when target temperatures
are set to different values, the object with a higher target temperature
usually reaches a position that is closer to the center of the heat source;
This situation will be further discussed in the Section IV-E. For the case
when target temperatures are set to the same value, the final position of the
circular aluminum sheet is always closer to the center of the heat source.
This phenomenon could be explained by the fact that absorptance of a metal is
usually much smaller than the absorptance of non-metallic materials (e.g. PLA)
[31].

###  IV-D Experiments with a Moving Heat Source

Figure 14: Evolution of the temperature errors of the three objects in the $8$
experiments with the adaptive controller, $\Delta T^{1}$, $\Delta T^{2}$, and
$\Delta T^{3}$ (measured in $\text{\,}\mathrm{\SIUnitSymbolCelsius}$). The
target temperatures are set as:
$\boldsymbol{\tau}^{*1}=\begin{bmatrix}40&40&40\end{bmatrix}^{{\intercal}}$,
$\boldsymbol{\tau}^{*2}=\begin{bmatrix}50&40&40\end{bmatrix}^{{\intercal}}$,
$\boldsymbol{\tau}^{*3}=\begin{bmatrix}50&50&35\end{bmatrix}^{{\intercal}}$,
$\boldsymbol{\tau}^{*4}=\begin{bmatrix}60&60&30\end{bmatrix}^{{\intercal}}$,
$\boldsymbol{\tau}^{*5}=\begin{bmatrix}50&45&35\end{bmatrix}^{{\intercal}}$,
$\boldsymbol{\tau}^{*6}=\begin{bmatrix}40&45&50\end{bmatrix}^{{\intercal}}$,
$\boldsymbol{\tau}^{*7}=\begin{bmatrix}45&45&45\end{bmatrix}^{{\intercal}}$,
$\boldsymbol{\tau}^{*8}=\begin{bmatrix}45&45&45\end{bmatrix}^{{\intercal}}$.
Figure 15: Spatial displacements of the three objects in the $8$ experiments
with the adaptive controller visualized from two viewing angles. We use a
blue-to-red color gradient to visualize the cold-to-hot change of temperatures
during the experiments.

In this section, we report an integrated experiment where the adaptive
controller is combined with an online ArUco tracking algorithm to achieve
temperature regulation while the heat source is changing. During setup, the
ArUco markers attached to the heat source are used to geometrically calibrate
the robot with the source, which is essential to compute the thermal
interaction matrix $\mathbf{L}$. Thus, when these relation are uncertain,
their new configurations have to be updated simultaneously.

Figure 16: Experiments with the adaptive controller with disturbances.

Here, we study the case where the robot is fixed and the heat source is
manually moved. We track marker $2$ attached to the source to obtain its
configuration. The target temperature vector is set to
$\boldsymbol{\tau}=[40,40,40{]}^{\intercal}$
$\text{\,}\mathrm{\SIUnitSymbolCelsius}$. Fig. 16 shows: (a) The initially
calibrated set up, (b) the manual movement applied to the source, (c) the
detected marker $2$ when the source is moving, and (d) the evolution of the
individual temperature errors. These results show that by continuously
updating the source-object pose, the control of the individual temperature
errors is not significantly affected and that $\|\Delta\boldsymbol{\tau}\|$
can still be asymptotically minimized. This experiment demonstrates how our
new thermal servoing method can be combined with other traditional controllers
(visual servoing in this case) to extend the sensorimotor capabilities of a
robot [48].

###  IV-E View Factor Visualization

In previous sections, we designed thermal controllers based on derived heat
transfer models. However, the models that relate $F_{21}$ and $\mathbf{x}$ are
generally complex. Therefore, part of the controlled system behaves as a
“black box” to the user. To investigate these aspects, in this section, we
introduce the visualization of the view factor $F_{21}$ with respect to the
end-effector configuration $\mathbf{x}$ as a useful tool for analyzing
radiation-based thermal servoing problems.

As an example, we take the “circular surfaces in arbitrary configurations”
case discussed in Section II-F. We implement the controlled variable method
[49] to split the 6-DOF pose $\mathbf{x}$ into two subsets: One where the
translation coordinates $p_{1}$, $p_{2}$, and $p_{3}$ (measured in
$\text{\,}\mathrm{cm}$) are the controlled variables, and another where the
rotation coordinates $\theta_{x}$, $\theta_{y}$, $\theta_{z}$ (measured in
degrees) are controlled variables. In the translation subset, rotations are
set to constant values of $\theta_{i}=0$, then, we compute $F_{21}$ for points
in a selected working range of controlled variables $p_{1},p_{2}\in[-20,20]$,
$p_{3}\in[0,30]$ with a step of $1$ (with $48,000$ points in total). In the
rotation subset, translations are similarly set to fixed constant values
$p_{1,2}=0$ and $p_{3}=$5\text{\,}\mathrm{cm}$$, then, points in the range of
$\theta_{x},\theta_{y},\theta_{z}\in[-90,90]$ are computed with an incremental
step of $2$ (with $729,000$ points in total).

We use the isosurface visualization tool provided by Plotly to visualize the
data. The translation and rotations subsets are shown in Fig. 17, where 3-DOF
end-effector configurations are represented by points in space, and the view
factor values are represented by isosurfaces with different colors (the
isosurfaces are formed by points which have the same or very close values of
$F_{21}$). This visualization method is inspired by the approximation of the
interaction matrix in (26), which reveals that $\mathbf{L}$ is positive
proportional to the directional derivative of $F_{21}(\mathbf{x})$ along
$\mathbf{x}$ as
$\mathbf{L}^{\intercal}=\lambda_{1}\nabla_{\mathbf{x}}F_{21}(\mathbf{x})$.
According to the definition of isosurface, the surface normal of every point
on the surface also points in $\nabla_{\mathbf{x}}F_{21}(\mathbf{x})$
direction. In addition, the interval distance between isosurfaces with an
equal value difference (also called “isosurface interval”) reveals the
magnitudes of the elements of $\nabla_{\mathbf{x}}F_{21}(\mathbf{x})$; A
larger distance represents a smaller magnitude.

As an example, let us analyze the translation subset shown in Fig. 17 (a). For
this single-object scenario, the normal vector at a point on the isosurface
indicates the direction of the end-effector movement (as computed from the
thermal controls (35), (43)) at that point. There are some characteristics of
these isosurfaces that can be intuitively deducted from the setup, e.g., the
symmetric spatial distribution of $F_{21}$ (due to the circular shape of the
heat source), and the proportionality of values of $F_{21}$ with respect to
the source-object separation.

However, the visualization provides two useful pieces of information. First,
that the centers of the incomplete spherical isosurfaces shift upwards when
$F_{21}$ decreases, which means that at some points, movement in the
$\vec{\mathbf{k}}_{1}$ direction will cause a decrease of $F_{21}$ (which
seems counter-intuitive). See e.g. $\mathbf{c}_{2}$ on the $F_{21}=0.1$
isosurface in Fig. 17 (a), which shows that in that configuration, the end-
effector needs to move backwards along the $\widehat{\mathbf{k}}_{1}$
direction to heat up faster. Second, the isosurface intervals at regions that
are farther from the heat source center are comparatively larger, which
indicates that the end-effector will move comparatively faster in those
regions. Similarly, Fig. 17 (b) shows the (much simpler) case where angles are
varied at a fixed position.

To further analyze the effect of robot motion on heat transfer process, we now
present a 4-DOF visualization result in Fig. 18. Based on the previous 3-DOF
translation subset, we add one more controlled variable $\theta_{x}$ to extend
it to 4-DOF. We vary the rotation $\theta_{x}$ from $0$ to
$90$\text{\,}\mathrm{\SIUnitSymbolCelsius}$$ and depict the change of
translation subset isosurfaces in Fig. 18 (a), (b), (c). For the previous
rotation subset, we add one more controlled variable $p_{3}$ and vary it from
$1$ to $10$ $\text{\,}\mathrm{cm}$, as shown in Fig. 18 (d), (e), (f).
Similarly, we vary $p_{1}$ from $0$ to $10$ $\text{\,}\mathrm{cm}$ and depict
the change of rotation subset isosurfaces in Fig. 18 (g), (h), (i). An
animation of this experiment is included in the supplementary video.

Figure 17: Isosurfaces visualization of the two view factor subsets: (a)
Translation subset and (b) rotation subset. Figure 18: Controlled variable
visualization of the view factor isosurfaces.

###  IV-F Unfeasible Thermal Targets

In Section III-C, we discussed two necessary but not sufficient conditions for
feasible targets. When one of the two conditions is not fulfilled, the
temperature error cannot be minimized to zero. This section reports and
analyzes two failed experiments with the proposed adaptive controller where
the target temperatures are set to
$\boldsymbol{\tau}^{*}=[80,80,80{]}^{\intercal}$ and
$\boldsymbol{\tau}^{*}=[70,35,35{]}^{\intercal}$
$\text{\,}\mathrm{\SIUnitSymbolCelsius}$; The temperature errors of each
coordinate $\Delta T^{i}$ are depicted in Fig. 19 (a), (c). The evolution of
the error $\|\Delta\boldsymbol{\tau}\|$ for the two experiments is shown in
Fig. 19 (b), (d). In this experimental study, we found that when all the
individual target temperatures are set relatively “high” ($80$
$\text{\,}\mathrm{\SIUnitSymbolCelsius}$ in this example), its corresponding
errors converge to a local minimum. Also, when the difference between
individual target temperatures is too large, one of the objects might more
closely reach its target, while the other will present steady-state errors.

In Section III-C, we prove that the steady-state temperature of an object
heated by a radiative source is directly proportional to $F_{21}$. Thus, the
geometry of view factor isosurfaces is a useful tool for analyzing the such
reachability conditions. Here, we discuss a simple but representative case
where two aluminum circular sheets with radius
$r_{o1}=$1.5\text{\,}\mathrm{cm}$$, $r_{o2}=$4.5\text{\,}\mathrm{cm}$$ are
attached to the end-effector at $\mathbf{o}_{1}$ and $\mathbf{o}_{2}$ (see
Fig. 20) and heated by a source with
$T_{1}=$200\text{\,}\mathrm{\SIUnitSymbolCelsius}$$. The center of the end-
effector is at $\mathbf{o}_{e}$, and $l_{e1}=l_{e2}=$2\text{\,}\mathrm{cm}$$
are the distances between the centers $\mathbf{o}_{1}$ and $\mathbf{o}_{2}$
and the end-effector $\mathbf{o}_{e}$. The view factors of the objects are
calculated based on the same setup as in previous sections. We use the
visualization method where three translations are the controlled variables,
for a parallel object and source surfaces.

By using the expression (51) and assuming that the thermophysical properties
are the same as mentioned in Section IV-A, the view factor values
corresponding to steady-state temperatures
$$30\text{\,}\mathrm{\SIUnitSymbolCelsius}$,$40\text{\,}\mathrm{\SIUnitSymbolCelsius}$,$50\text{\,}\mathrm{\SIUnitSymbolCelsius}$$
are calculated as $0.12$, $0.37$, and $0.65$. According to this one-to-one
correspondence between the isosurface and the steady-state temperature, to
automatically reach the target temperature $T^{*}$ can be geometrically
interpreted as positioning the object center over the isosurface that
corresponds to $T^{*}$. Similarly, determining the feasibility of target
temperatures $T^{*1}$ and $T^{*2}$ of two objects attached to the same end-
effector is identical to finding whether there exists an end-effector pose
that places both objects onto their “desired isosurfaces”.

An example is shown in Fig. 21, where we denote the steady-state temperatures
of objects 1 and 2 by $T_{ss}^{1}$ and $T_{ss}^{2}$. Fig. 21 (a)–(b) show the
steady-state temperature isosurfaces of objects 1 and 2 where
$T_{ss}^{1}=T_{ss}^{12}=$30\text{\,}\mathrm{\SIUnitSymbolCelsius}$,$40\text{\,}\mathrm{\SIUnitSymbolCelsius}$,$50\text{\,}\mathrm{\SIUnitSymbolCelsius}$$.
Since the two objects are circular plates with different radii, the shapes of
their isosurfaces are slightly different; We use red and blue color to
differentiate them, and are jointly depicted in the same coordinate system in
Fig. 21 (c).

Figure 19: Evolution of the temperature errors with unfeasible target
temperatures: (top) $\boldsymbol{\tau}^{*}=[80,80,80{]}^{\intercal}$
$\text{\,}\mathrm{\SIUnitSymbolCelsius}$, and (bottom)
$\boldsymbol{\tau}^{*}=[70,35,35{]}^{\intercal}$.

In Fig. 21 (d), (e), (f), different combinations of target temperatures
$T^{1*},T^{2*}$ and their corresponding isosurfaces are depicted. These
figures graphically demonstrate how for thermal targets
$T^{*1}=T^{*2}=$30\text{\,}\mathrm{\SIUnitSymbolCelsius}$$ (depicted in Fig.
21 (d)), and $T^{*1}=$30\text{\,}\mathrm{\SIUnitSymbolCelsius}$$,
$T^{*2}=$40\text{\,}\mathrm{\SIUnitSymbolCelsius}$$ (depicted in Fig. 21 (e)),
the end-effector can position the objects into their desired final isosurfaces
(corresponding to their target steady-state temperatures); The initial
position of this trajectory is colored in blue, the final in red. However, for
the case where $T^{*1}=$50\text{\,}\mathrm{\SIUnitSymbolCelsius}$$ and
$T^{*2}=$30\text{\,}\mathrm{\SIUnitSymbolCelsius}$$, Fig. 21 (f) graphically
demonstrates that the minimum distance between the two isosurfaces is larger
than $l_{e1}+l_{e2}$; Therefore, these target temperatures $T^{*1}$ and
$T^{*2}$ are _not feasible_. Similarly, if $l_{e1}+l_{e2}$ is larger than the
maximum distance between two target isosurfaces, that combination of $T^{*1}$
and $T^{*2}$ is also unfeasible.

In general, thermophysical properties, view factors, and fixed spatial
relationships between objects are the main three factors that determine the
feasibility of thermal targets. The geometric interpretation of the
feasibility problem might also be useful for path planning-like algorithms
dealing with thermal servoing problems.

The accompanying multimedia file demonstrates the performance of our new
control methodology with multiple experimental results.

Figure 20: Conceptual illustration of two objects fixed to an end-effector for
analyzing unfeasible target temperatures. Figure 21: Geometric explanation of
the target temperature feasibility using the steady-state temperature
isosurfaces.

##  V CONCLUSION

In this paper, we present a new robotic temperature control technique based on
heat radiation to automatically regulate temperatures of multiple objects. For
that, we provide a comprehensive formulation of different scenarios of thermal
servoing problems. Two asymptotically stable controllers, one model-based and
one adaptive, are designed and validated by a series of experiments where
temperatures of three different objects are independently regulated. We also
discussed potential applications of the isosurface visualization, such as
analyzing the geometry of the seemingly invisible heat transfer process.

The key concept of the proposed method is to exploit the geometric-thermal-
motor relations between the heat source and the surface for automatically
computing motion controls. This advanced feedback control capability is needed
to improve the performance of many economically-important applications.
However, from the point of view of generality, the proposed algorithm has many
limitations, since we only consider cases where heat radiation is dominant.
For cases where other heat transfer modes are dominant or comparable (objects
in contact with non-adiabatic surfaces, electrical equipment cooled by high
speed air flow, human skin treated by laser thermal excitation, food heated up
in a pan, etc.), different heat transfer models need to be analyzed. Another
possible solution is to implement model-free control algorithms that primarily
rely on collected data instead of analytical models (which in our case were
formulated based on fundamental physical principles).

For future work, we would like to integrate thermal servoing with existing
visual and proximity servoing algorithms; This multimodal perceptual and
control capability is essential for developing advanced robotic temperature
control systems in complex scenarios, such as service tasks in human
environments and and intelligent industrial manufacturing. Our team is
currently working towards developing algorithms which simultaneously consider
the three basic modes of heat transfer. For this situation, thermal images
(which provide detailed temperature profiles of an object surface) may be used
as a sensing system. We encourage interested readers to work along these open
research directions.

##  VI APPENDIX

###  VI-A Online Estimation of the Interaction Matrix

We use the numerical definite integral solver provided by SciPy [39] to
approximate the interaction matrix in real time. The essence of numerical
integral is to divide the complex integrand into small subsections, and
approximate each subsection with a polynomial that is easy to integrate. A
commonly used method is the composite Simpson’s rule[50], which approximates
the subsection by quadratic polynomials. The general form of composite
Simpson’s rule is as follows:

$\begin{split}&\int_{c}^{d}g(x)\mathop{}\\!\mathrm{d}x\approx\\\
&\frac{h}{3}\left[g\left(x_{0}\right)+2\sum_{j=1}^{n/2-1}g\left(x_{2j}\right)+4\sum_{j=1}^{n/2}g\left(x_{2j-1}\right)+g\left(x_{n}\right)\right]\end{split}$
(59)

where $n$ is the number of subintervals, $x_{j}=c+jh$ for
$j=\mathrm{0,1,...,n-1,n}$ with $h=(d-c)/n$; in particular, $x_{0}=c$ and
$x_{n}=d$.

###  VI-B Evaluation of the Estimated Interaction Matrix

In Section II-H, we propose a new method that uses truncated Fourier series to
approximate the view factor between irregularly shaped objects. Since the
interaction matrices are estimated according to the view factor, the accuracy
and computation time of the approximation directly affects the performance of
the system. In this section, we conduct a case study where the standard
discrete surface integral (DSI) method[42] and our proposed method are used to
calculate the view factor between the bunny object and the circular heat
source used in experiments. The center of the heat source is
$\mathbf{c}_{1}=\begin{bmatrix}0&0&0\end{bmatrix}^{{\intercal}}$ and the
center of the bunny object is
$\mathbf{c}_{2}=\begin{bmatrix}0&0&5\end{bmatrix}^{{\intercal}}$ (measured in
$\text{\,}\mathrm{cm}$). We denote the number of discrete surfaces by $n$ and
the number of harmonics by $h$.

Since there is no explicit formula to calculate the view factor of irregularly
shaped objects, we take the estimated value $F_{21}=0.764602$ obtained by the
DSI method when $n$ is set to $20000$ as the ground truth. The ground truth is
then used to compute the error of the estimation when using different values
of $n$ and $h$. The comparison result in terms of accuracy and computation
time is depicted in Fig. 22. We find that the result obtained by the proposed
method converges to the ground truth in a short period of time. Thus, a
desirable estimation of the interaction matrix can be obtained in real time.

Figure 22: Comparison between the discrete surface integral method and the
proposed truncated Fourier series algorithm.

###  VI-C Polynomial Fitting with Sliding Window

Due to the limitation of the sensor accuracy and sampling rate, the raw
feedback temperatures are discrete and noisy. This is undesirable for directly
estimating the true values of $\mathbf{v}$. To this end, we implement
polynomial fitting with a sliding window method, to estimate the values of
$\boldsymbol{\tau}$ and $\mathbf{v}$ from the raw feedback data. For a single
object, we denote its temperature as $T_{t}$ where $t$ is a variable of
sampling time. A series of sampling times is denoted by
$t_{1},t_{2},\cdots,t_{n}$. The size of the sliding window is set to be 10
data samples. When $n<10$, the estimation is at the initialization stage and
the robot will not move. When $n\geq 10$, we denote the collected 10
temperature values by a vector
$\mathbf{Y}_{s}=\begin{bmatrix}T_{t_{n-9}}&T_{t_{n-8}}&\cdots&T_{t_{n}}\end{bmatrix}^{{\intercal}}$.
We then fit 10 sample points in the sliding window to the polynomial of order
3 as:

$\widehat{\mathbf{Y}}_{s}=\mathbf{P}\mathbf{c}^{t_{n}}=\begin{bmatrix}(t_{n-9})^{3}&(t_{n-9})^{2}&(t_{n-9})&1\\\
(t_{n-8})^{3}&(t_{n-8})^{2}&(t_{n-8})&1\\\ \vdots&\vdots&\vdots\\\
(t_{n})^{3}&(t_{n})^{2}&(t_{n})&1\end{bmatrix}\begin{bmatrix}c_{3}^{t_{n}}\\\
c_{2}^{t_{n}}\\\ c_{1}^{t_{n}}\\\ c_{0}^{t_{n}}\\\ \end{bmatrix}$ (60)

where $\mathbf{c}^{t_{n}}$ is the coefficient vector of the polynomial at
sample time $t_{n}$. To minimize
$\|\mathbf{Y}_{s}-\widehat{\mathbf{Y}}_{s}\|^{2}$, the coefficients are
computed as $\mathbf{c}^{t_{n}}={\mathbf{P}}^{+}\mathbf{Y}_{s}$. Then, the
temperature and temperature rate of a single object can be estimated as:

$\displaystyle
T_{t_{n}}=c^{t_{n}}_{3}(t_{n})^{3}+c^{t_{n}}_{2}(t_{n})^{2}+c^{t_{n}}_{1}(t_{n})+c^{t_{n}}_{0}$
(61) $\displaystyle
v_{t_{n}}=3c^{t_{n}}_{3}(t_{n})^{2}+2c^{t_{n}}_{2}(t_{n})+c^{t_{n}}_{1}$ (62)

When the new data is obtained, we update $\widehat{\mathbf{Y}}_{s}$ and
$\mathbf{P}$. By following the same procedure, object temperature rate can be
estimated accordingly.

## References

  * [1] D. Navarro-Alarcon, Y.-H. Liu, J. G. Romero, and P. Li, “On the visual deformation servoing of compliant objects: Uncalibrated control methods and experiments,” _Int. Journal of Robotics Research_ , vol. 33, no. 11, pp. 1462–1480, 2014.
  * [2] A. Cherubini and F. Chaumette, “Visual navigation of a mobile robot with laser-based collision avoidance,” _Int. Journal of Robotics Research_ , vol. 32, no. 2, pp. 189–205, 2013.
  * [3] M. Tirindelli, M. Victorova, J. Esteban, S. T. Kim, D. Navarro-Alarcon, _et al._ , “Force-ultrasound fusion: Bringing spine robotic-us to the next “level”,” _IEEE Robot. Autom. Lett._ , vol. 5, no. 4, pp. 5661–5668, 2020.
  * [4] A. Magassouba, N. Bertin, and F. Chaumette, “Audio-based robot control from interchannel level difference and absolute sound energy,” in _Proc. IEEE Int. Conf. Intelligent Robots and Systems_ , 2016, pp. 1992–1999.
  * [5] F. Rahbar, A. Marjovi, P. Kibleur, and A. Martinoli, “A 3-d bio-inspired odor source localization and its validation in realistic environmental conditions,” in _IEEE/RSJ Int. Conf. on Intelligent Robots and Systems_ , 2017, pp. 3983–3989.
  * [6] R. Gade and T. B. Moeslund, “Thermal cameras and applications: A survey,” _Machine Vision and Applications_ , vol. 25, pp. 245–262, 2014.
  * [7] E. Benli, R. L. Spidalieri, and Y. Motai, “Thermal multisensor fusion for collaborative robotics,” _IEEE Transactions on Industrial Informatics_ , vol. 15, no. 7, pp. 3784–3795, 2019.
  * [8] Y. Cao, B. Xu, Z. Ye, J. Yang, Y. Cao, _et al._ , “Depth and thermal sensor fusion to enhance 3d thermographic reconstruction,” _Opt. Express_ , vol. 26, no. 7, pp. 8179–8193, Apr 2018.
  * [9] Y. He and R. Yang, “Eddy current volume heating thermography and phase analysis for imaging characterization of interface delamination in cfrp,” _IEEE Trans. on Industrial Informatics_ , vol. 11, no. 6, pp. 1287–1297, 2015\.
  * [10] W. W.-L. Lai, K.-K. Lee, and C.-S. Poon, “Validation of size estimation of debonds in external wall’s composite finishes via passive infrared thermography and a gradient algorithm,” _Construction and Building Materials_ , vol. 87, pp. 113 – 124, 2015.
  * [11] A. Cherubini and D. Navarro-Alarcon, “Sensor-based control for human-robot collaboration: Fundamentals, challenges and opportunities,” _Front. in Neurorobotics (in press)_ , vol. 1, no. 1, pp. 1–21, 2020.
  * [12] W. Binrui, J. Yinglian, X. Hong, and W. Ling, “Temperature control of pem fuel cell stack application on robot using fuzzy incremental pid,” in _2009 Chinese Control and Decision Conference_. IEEE, 2009, pp. 3293–3297.
  * [13] M. Ho and J. P. Desai, “Towards a mri-compatible meso-scale sma-actuated robot using pwm control,” in _2010 3rd IEEE RAS & EMBS International Conference on Biomedical Robotics and Biomechatronics_. IEEE, 2010, pp. 361–366.
  * [14] J. De Backer, G. Bolmsjö, and A.-K. Christiansson, “Temperature control of robotic friction stir welding using the thermoelectric effect,” _The International Journal of Advanced Manufacturing Technology_ , vol. 70, no. 1-4, pp. 375–383, 2014.
  * [15] J.-H. Kim and B. Y. Lattimer, “Real-time probabilistic classification of fire and smoke using thermal imagery for intelligent firefighting robot,” _Fire Safety Journal_ , vol. 72, pp. 40–49, 2015.
  * [16] A. Cammarata, “Optimized design of a large-workspace 2-dof parallel robot for solar tracking systems,” _Mechanism and machine theory_ , vol. 83, pp. 175–186, 2015.
  * [17] M. Muddassir, D. Gomez, L. Hu, S. Chen, and D. Navarro-Alarcon, “Robotics meets cosmetic dermatology: Development of a novel vision-guided system for skin photo-rejuvenation,” _IEEE/ASME Trans. Mechatronics (under review)_ , vol. 1, no. 1, pp. 1–15, 2020.
  * [18] D. McKemy, “Temperature sensing across species,” _Eur. J. Physiol._ , vol. 454, pp. 777–791, 2007.
  * [19] J. Mlynek, R. Knobloch, and R. Srb, “Optimization of a heat radiation intensity and temperature field on the mould surface.” in _ECMS_ , 2016, pp. 425–431.
  * [20] A. Imdoukh, A. Shaker, A. Al-Toukhy, D. Kablaoui, and M. El-Abd, “Semi-autonomous indoor firefighting uav,” in _Int. Conf. on Advanced Robotics_ , 2017, pp. 310–315.
  * [21] G. Muscato, F. Bonaccorso, L. Cantelli, D. Longo, and C. D. Melita, “Volcanic environments: Robots for exploration and measurement,” _IEEE Robotics & Automation Magazine_, vol. 19, no. 1, pp. 40–49, 2012.
  * [22] M. Fu, W. Weng, W. Chen, and N. Luo, “Review on modeling heat transfer and thermoregulatory responses in human body,” _Journal of Thermal Biology_ , vol. 62, pp. 189 – 200, 2016, modeling bioheat transfer processes and thermoregulatory responses.
  * [23] R. Lee, “Thermal tolerances of deep-sea hydrothermal vent animals from the northeast pacific,” _Biological Bulletin_ , vol. 205, no. 2, pp. 98–101, 2003.
  * [24] C.-C. Tsai, R. A. Childers, N. N. Shi, C. Ren, J. N. Pelaez, _et al._ , “Physical and behavioral adaptations to prevent overheating of the living wings of butterflies,” _Nature communications_ , vol. 11, no. 1, pp. 1–14, 2020.
  * [25] M. L. May, “Insect thermoregulation,” _Annual review of entomology_ , vol. 24, no. 1, pp. 313–349, 1979.
  * [26] L. Natale, G. Metta, and G. Sandini, “Development of auditory-evoked reflexes: Visuo-acoustic cues integration in a binocular head,” _Rob. Auton. Syst._ , vol. 39, no. 2, pp. 87 – 106, 2002.
  * [27] G. Arechavaleta, J. Laumond, H. Hicheur, and A. Berthoz, “An optimality principle governing human walking,” _IEEE Trans. Robot._ , vol. 24, no. 1, pp. 5–14, 2008.
  * [28] S. Na, Y. Qiu, A. E. Turgut, J. Ulrich, T. Krajník, _et al._ , “Bio-inspired artificial pheromone system for swarm robotics applications,” _Adaptive Behavior_ , pp. 1–21, 2020.
  * [29] T. L. Bergman, F. P. Incropera, D. P. DeWitt, and A. S. Lavine, _Fundamentals of heat and mass transfer_. John Wiley & Sons, 2011.
  * [30] M. F. Modest, _Radiative heat transfer_. Academic press, 2013.
  * [31] I. Lienhard and H. John, _A heat transfer textbook_. phlogiston press, 2005.
  * [32] M. Vujičić, N. Lavery, and S. Brown, “View factor calculation using the monte carlo method and numerical sensitivity,” _Communications in numerical methods in Engineering_ , vol. 22, no. 3, pp. 197–203, 2006.
  * [33] J. C. Chai, J. P. Moder, and K. C. Karki, “A procedure for view factor calculation using the finite-volume method,” _Numerical Heat Transfer: Part B: Fundamentals_ , vol. 40, no. 1, pp. 23–35, 2001.
  * [34] S. C. Mishra, A. Shukla, and V. Yadav, “View factor calculation in the 2-d geometries using the collapsed dimension method,” _International communications in heat and mass transfer_ , vol. 35, no. 5, pp. 630–636, 2008\.
  * [35] M. Vujičić, N. Lavery, and S. Brown, “Numerical sensitivity and view factor calculation using the monte carlo method,” _Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science_ , vol. 220, no. 5, pp. 697–702, 2006.
  * [36] E. M. Sparrow, _Radiation heat transfer_. Routledge, 2018.
  * [37] V. R. Rao and V. Sastri, “Efficient evaluation of diffuse view factors for radiation,” _International journal of heat and mass transfer_ , vol. 39, no. 6, pp. 1281–1286, 1996.
  * [38] H. Flanders, “Differentiation under the integral sign,” _The American Mathematical Monthly_ , vol. 80, no. 6, pp. 615–627, 1973.
  * [39] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, _et al._ , “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,” _Nature Methods_ , vol. 17, pp. 261–272, 2020.
  * [40] G. P. Tolstov, _Fourier series_. Courier Corporation, 2012.
  * [41] D. Navarro-Alarcon and Y.-H. Liu, “Fourier-based shape servoing: a new feedback method to actively deform soft objects into desired 2-d image contours,” _IEEE Transactions on Robotics_ , vol. 34, no. 1, pp. 272–279, 2018.
  * [42] T. Muneer, S. Ivanova, Y. Kotak, and M. Gul, “Finite-element view-factor computations for radiant energy exchanges,” _Journal of Renewable and Sustainable Energy_ , vol. 7, no. 3, p. 033108, 2015.
  * [43] Y. Nakamura, _Advanced robotics: redundancy and optimization_. Boston, MA: Addison-Wesley Longman, 1991.
  * [44] M. Vidyasagar, _Nonlinear systems analysis_. SIAM, 2002.
  * [45] N. Titchener, S. Colliss, and H. Babinsky, “On the calculation of boundary-layer parameters from discrete data,” _Experiments in Fluids_ , vol. 56, no. 8, p. 159, 2015.
  * [46] Y.-H. Liu, H. Wang, C. Wang, and K. K. Lam, “Uncalibrated visual servoing of robots using a depth-independent interaction matrix,” _IEEE Trans. Robot._ , vol. 22, no. 4, pp. 804–817, Aug. 2006.
  * [47] S. Garrido-Jurado, R. Muñoz-Salinas, F. J. Madrid-Cuevas, and M. J. Marín-Jiménez, “Automatic generation and detection of highly reliable fiducial markers under occlusion,” _Pattern Recognition_ , vol. 47, no. 6, pp. 2280–2292, 2014.
  * [48] D. Navarro-Alarcon, J. Qi, J. Zhu, and A. Cherubini, “A Lyapunov-stable adaptive method to approximate sensorimotor models for sensor-based control,” _Front. in Neurorobotics_ , vol. 14, no. 59, pp. 1–12, 2020\.
  * [49] C. Williams, “Research methods,” _Journal of Business & Economics Research (JBER)_, vol. 5, no. 3, 2007.
  * [50] D. Hai and R. Smith, “An elementary proof of the error estimates in simpson’s rule,” _Mathematics Magazine_ , vol. 81, no. 4, pp. 295–300, 2008.

|  Luyin Hu received his BEng. degree in Mechanical Engineering from The Hong
Kong Polytechnic University, KLN, Hong Kong in 2019. Currently, he is pursuing
an MPhil degree in Mechanical Engineering at the same university, where he
also works as a Research Assistant. His research interests include multimodal
robot perception, servomechanisms, and control system design.
---|---
|  David Navarro-Alarcon (GS’06–M’14–SM’19) received the Ph.D. degree in
Mechanical and Automation Engineering from The Chinese University of Hong Kong
(CUHK), Shatin, Hong Kong, in 2014. He was an Assistant (Research) Professor
at the CUHK T Stone Robotics Institute, from 2015 to 2017. Since 2017, he has
been with The Hong Kong Polytechnic University, KLN, Hong Kong, where he is an
Assistant Professor at the Department of Mechanical Engineering. His current
research interests include perceptual robotics and control theory.
---|---
|  Andrea Cherubini received the M.Sc. in Mechanical Engineering in 2001 from
the University of Rome La Sapienza and a second M.Sc. in Control Systems in
2003 from the University of Sheffield, U.K. In 2008, he received the Ph.D. in
Control Systems from La Sapienza. From 2008 to 2011, he was postdoc at INRIA
Rennes. Since 2011, he is Associate Professor at Université de Montpellier.
His research interests include: physical human–robot interaction, and
manipulation of soft objects.
---|---
|  Mengying Li received her PhD degree in Mechanical and Aerospace Engineering
from University of California San Diego. From 2018 to 2020, She was a
Postdoctoral Scholar in the Center for Energy Research of UC San Diego. She
joined The Hong Kong Polytechnic University as an Assistant Professor in the
Department of Mechanical Engineering since 2020. Her research interests
include atmospheric radiative heat transfer for renewable energy integration
and design of multi-generation systems.
---|---