\section{Introduction.}
\hspace{0.5cm}
Recently the Borexino collaboration \cite{agost} published the results of analysis of the experimental energy spectrum of single events. In this work the authors fitted the experimental spectrum with the theoretical one.  As a result, the interaction rates from various sources were obtained except for CNO solar neutrinos interaction rate.

To obtain the values of the interaction rates from all sources, it is necessary to have the correct theoretical information about the shape of the measured energy spectrum from each of the possible sources. In the absence of correct information, the used procedure for fitting of experimental data with the model may not allow to obtain the interaction rates and their uncertainties from the sources used in the model.

The results of \cite{agost} most likely demonstrate the lack of complete information on possible sources of single events. The paper reports that there is no way to get the interaction rate of CNO neutrinos. Only the upper bound of $8\div10$ events per day in 100 tons of scintillator was obtained.  From 3 to 5 events per day are expected from CNO neutrino interactions depending on the model of the Sun.

\section{Energy spectrum of recoil electrons from scattering of $^{40}$K geo-antineutrino.}
\hspace{0.5cm}
The authors of the work \cite{bezr} considered an additional source of single events for the Borexino detector. This is the scattering of potassium geo-antineutrinos by electrons. Such a source was not considered in \cite{agost} because the Silicate Earth model was used in it. The concentration of potassium in this model is very small (0.024\% of the Earth mass) and the contribution to single events of the Borexino detector from potassium geo-antineutrinos is negligible.

\includegraphics[scale=0.7]{40KBorexinofig11}

\small{Fig.1. Differential energy spectra of recoil electrons from neutrino  (or antineutrino) scattering on electrons for various sources (the oscillations were not taken into account). The red line is the solar neutrino from $^{13}$N decay in CNO cycle. The green is the solar neutrino from $^{15}$O decay in CNO cycle. The blue is the solar neutrino from $^{17}$F decay in CNO cycle. The black is $^{40}$K geo-antineutrino for the concentration of potassium in modern Earth 1\% of the Earth mass}

\includegraphics[scale=0.7]{40KBorexinofig22}

\small{Fig.2. Differential energy spectra of recoil electrons from solar CNO neutrino scattering by scintillator electrons in 100 tons of scintillator (red curve) and the sum of the recoil electron spectra from the scattering of solar CNO neutrino  and from the scattering of  $^{40}$K geo-antineutrino for the concentration of potassium in the modern Earth from 1\% to 2\% of the Earth mass. The oscillations were taken into account.}
\vspace{0.5cm}

\hspace{0.5cm}
We proceeded from the prediction of the Earth elemental composition by another model of the Earth "Hydridic Earth" or "Hydrogen-rich Earth" \cite{Larin, toul}. These studies predict concentration of potassium in the original Earth up to 3.8\% of the Earth mass. In the modern Earth we can expect a larger proportion of potassium due to the loss of hydrogen by the Earth during the time of it's existence. The predictions of the Hydridic Earth model have low accuracy, so it is reasonable to search for the flux of potassium geo-antineutrinos in a wide range of values.

The work \cite{bezr}  presented the results of calculation of the energy spectrum of recoil electrons from geo-antineutrino scattering from $^{40}$K decay on scintillator electrons and the counting rate of these events in 100 tons of scintillator for the concentration of potassium in the modern Earth from 0.3\% to 3.8\% of the Earth's mass.

We present at figure 1 our calculation (using the approximations described in  \cite{bezr} and the oscillations were not taken into account) of differential energy spectra of recoil electrons produced by the scattering of neutrinos  (or antineutrinos) on electrons for different sources: red curve - solar neutrinos from the decay of $^{13}$N in the CNO cycle, the green curve corresponds to the solar neutrinos from the decay of $^{15}$O in the CNO cycle, blue curve - solar neutrinos from the decay of  $^{17}$F in the CNO cycle, black curve -  geo-antineutrinos from $^{40}$K decay for the concentration of potassium in the Earth 1\% of weight of the Earth. We can see that the spectrum of $^{40}$K geo-antineutrinos for the concentration of potassium in the Earth of 1\% of the Earth mass is similar to the spectrum of solar neutrinos from the decay of $^{13}$N in the CNO cycle. It has an upper limit larger than that of the spectrum of solar neutrinos from the decay of $^{13}$N. For our calculation we used spectra of solar neutrinos from the CNO cycle \cite{Bahcall} — original work of J. Bahcall, which are normalized to the value of the flux for large metallicity of the Sun ($4.88\cdot10^8$ cm$^{-2}$ s$^{-1}$).

Figure 2 shows the differential energy spectrum from CNO neutrino scattering by electrons in 100 tons of scintillator (red line) and the total recoil electron spectrum from the scattering of CNO neutrino and $^{40}$K geo-antineutrino for the concentration of potassium in the modern Earth from 1\% to 2\% of the Earth mass (green and blue curves respectively). The oscillations were taken into account. It means that the contribution of muon and tau neutrinos (and anty-neutrinos) scattering was taken into account.

 The intensity of events from the CNO neutrino calculated for the curve from figure 2 is 4.9 cpd/100 tons, and the additional intensities of events from $^{40}$K geo-antineutrinos is 1.6 to 3.4 cpd/100 tons.
The Borexino collaboration in \cite{agost} used a differential energy spectrum of recoil electrons from the scattering of  CNO  neutrinos on electrons, similar to the one shown in figure 2 in red, and did not obtain expected value of interaction rate ($3 \div 5$ cpd/100 tons).

We propose to perform the same study as in Borexino collaboration work \cite{agost} using the total spectrum of recoil electrons from the scattering of CNO neutrinos and from the scattering of $^{40}$K geo-antineutrino for the potassium concentration in the modern Earth from 1\% to 2\% of the Earth's mass. We expect to find a certain value for the intensity of total events from the CNO neutrinos and from the $^{40}$K geo-antineutrinos in the range of $6 \div 9$ events per day in 100 tons of scintillator. We also expect that, as a result of this analysis, several events per day will be taken away from events attributed in work  \cite{agost} to solar $^{7}$Be neutrinos interactions and $^{210}$Bi decays. The errors of the obtained intensities will become symmetrical.

We give the functions from figure 2 in table form (see table below).

The obtaining the intensity of total spectrum of recoil electrons from the scattering of CNO neutrinos and $^{40}$K geo-antineutrinos as a result of fitting analysis of Borexino single event data will be an important confirmation of the validity of the Hydridic Earth model.

\begin{table}[ht]
\caption{Differential energy spectra of recoil electrons from solar CNO neutrino scattering by scintillator electrons in 100 tons of scintillator and the same with addition of $^{40}$K antineutrinos in units $MeV^{-1} year^{-1} (100 tons)^{-1}$. }
\label{tabl:event}
\centering
\vspace{2mm}
 \begin{tabular}{| c | c | c | c | c |}
 \hline
 T, MeV &  CNO  & CNO+1\%($^{40}$K)  & CNO+1.5\%($^{40}$K) & CNO+2\%($^{40}$K) \\
 \hline
 0.1 & 2997.64 & 5254.39 & 6334.59 & 7414.79 \\
 0.2 & 2736.43 & 4363.73 & 5142.65 & 5921.57 \\
 0.3 & 2447.10 & 3597.83 & 4148.63 & 4699.42 \\
 0.4 & 2129.82 & 2924.45 & 3304.78 & 3685.11 \\
 0.5 & 1787.78 & 2319.67 & 2574.24 & 2828.80 \\
 0.6 & 1439.26 & 1781.15 & 1944.77 & 2108.38 \\
 0.7 & 1100.80 & 1308.16 & 1407.39 & 1506.62 \\
 0.8 & 787.671 & 901.87 & 956.52 & 1011.16 \\
 0.9 & 526.847 & 579.58 & 604.81 & 630.04 \\
 1.0 & 355.736 & 371.46 & 378.98 & 386.5 \\
 1.1 & 249.9 & 249.9 & 249.9 & 249.9 \\
 1.2 & 156.71 & 156.71 & 156.71 & 156.71 \\
 1.3 & 81.17 & 81.17 & 81.17 & 81.17 \\
 1.4 & 26.75 & 26.75 & 26.75 & 26.75 \\
 1.5 & 0.303 & 0.303 & 0.303 & 0.303 \\
 \hline
\end{tabular}
\vspace{10mm}
\end{table}

\section{Reasons to use the Hydridic Earth model.}

\vspace{0.5cm}
\hspace{0.5cm} The widespread belief in the fairness of Silicate Earth model and belief in the validity of the results of work \cite{davis} that the heat flux from the Earth interior is equal to 47$\pm$2 TW do not allow to include in analysis a new source of single events in Borexino - the $^{40}$K geo-antineutrinos interaction.

The idea that the Hydridic Earth model is not correct is also widely accepted. One of the objections is that the entire Earth could have melt due to radiogenic heat and spent most part of its life in this state if the Earth contains potassium more than 1\% of the Earth weight. The current heat flux from the Earth's interior can be more than 200 TW in the frame of Hydridic Earth model which contradicts to the result of the work \cite{davis}.

We have learned how to respond to the objections to Hydridic Earth model. The answers to the main objections are published in \cite{bezr2, bezr2018, barab2019, bezr22018}. Here we note that the entire Earth could not have melt because the Hydridic Earth model contains a subsurface cooling mechanism that is activated when the subsurface is heated enough to decompose the metal hydrides. Therefore, in the Hydridic Earth model, the subsurface temperature oscillates \cite{Larin}.

It is also noted in \cite{bezr2, bezr2018, barab2019, bezr22018}  that thermal conductivity is not the main mechanism of heat transfer in the Earth, but protons and hydrogen-containing gases carry out the heat away. In these works the experimental evidences are provided that the heat flux from Earth interior can reach the several hundreds TW. These are the heating of the oceans, the temperature profile of ultra-deep wells and non-direct evidence – Moon heat flux from interior.

Moreover, the effect was found that is predicted by the Hydridic Earth model and not predicted by the Silicate Earth model.
The Hydridic Earth model predicts that the Earth's crust is positively charged and contains a large amount of positive charge in the form of protons and positive ions of various hydrogen-containing gases. We tested the validity of this prediction experimentally by detecting an excess of the concentration of positive air-ions in underground rooms over the concentration of negative air-ions. In order to make this prediction, a new model of terrestrial electricity was developed based on the Hydridic Earth model \cite{bezr2019}. The model named "Hydridic model of terrestrial electricity".  This model explains the origin of the atmospheric electric field and all the observed effects of atmospheric electricity in a single way. The model also explains the origin of telluric currents. In \cite{bezr22019}, the Hydridic model of terrestrial electricity was successfully used in the analysis of the reaction of telluric currents to earthquakes.

 We consider Hydridic Earth model as an adequate tool for the study of the Earth.

\vspace{0.5cm}

\section{Conclusion.}
\begin{enumerate}
    \item
We propose to include the scattering of $^{40}$K geo-antineutrinos by electrons in the analysis of single event energy spectrum of Borexino detector.
\item
The Hydridic Earth model predicts the potassium concentrations up to 4\% of the Earth mass.
\item
It is expected that the total interaction rate of the solar CNO neutrino and $^{40}$K geo-antineutrino will be obtained in the range of $6 \div 9$ events per day in 100 tons of scintillator.
\item
We consider Hydridic Earth model as an adequate tool for the study of the Earth.
\end{enumerate}

\section{Acknowledgment}
The authors express their gratitude to M. D. Skorokhvatov for the opportunity to discuss this note at a seminar at the NRC "Kurchatov Institute".


\section{Introduction}
\label{sec:intro}

Multi-robot systems are well suited to solve complex tasks in dynamic and dangerous environments due to their redundancy, ability to operate in parallel, and system level fault tolerance to individual failure as highlighted in \cite{brambilla2013swarm, sahin2004swarm}. Multi-Robot Task Allocation (MRTA) deals with the assignment of agents to tasks in order to achieve an overall system goal within the constraints of the deployment setting. Therefore, in order to leverage the potential multi-robot systems have to successfully operate in dynamic and dangerous environments to solve complex problems such as disaster response, search and rescue, environmental monitoring, and automated warehousing, effective methods for solving the MRTA problem are needed (\cite{gerkey2004formal}).

The MRTA problem becomes more complex when introducing morphological or behavioral heterogeneity within a deployed multi-robot system (\cite{swarmanoid2013}). However, this additional complexity comes with the benefit of improving the overall system efficiency by leveraging the strengths of individual robots within the collective. For example, in a Search and Rescue (SaR) scenario, quadrotors, which are quick and agile, are better suited for scouting and surveying while ground robots are better suited for debris clearing and resource extraction. By leveraging these strengths efficiently and allocating tasks appropriately the heterogeneous system could out perform a system comprised of only aerial or ground robots.

 Moreover, in many scenarios the MRTA problem is accompanied by timing constraints where tasks must be performed sequentially, e.g. a robot must wait for a delivery before transporting the delivered package. Problems of this type are typically referred to as scheduling problems and involve an additional complexity. This class of problems can be solved with Mixed Integer Linear Programs (MILP) that attempt to schedule all the tasks at once, however, this approach suffers from an exponential complexity as noted in \cite{gombolay2013fast}. Additionally, when deploying a multi-agent system in dynamic environments, the system must be able to respond and reschedule tasks when unavoidable and inevitable environmental disturbances occur. Lastly and specific to the SaR operation following a natural disaster (e.g. wildfires, earthquakes, hurricanes) is that the team of agents must cover the targeted area and rescue victims amongst other tasks within a short-time window. This is due to the drastic decrease in the likelihood of victims surviving after $48$ hours as highlighted in \cite{48hours}. Therefore, the SaR problem can be cast as an instance of heterogeneous multi-agent system scheduling with the objective of minimizing the time of completion of all tasks (i.e. the makespan).

Inspired by the work in \cite{rhp}, this paper proposes a Receding Horizon Planning framework to solve the heterogeneous MRTA scheduling problem and demonstrates its effectiveness in a SaR application. Inspired by Model Predictive Control, at fixed time intervals, the proposed framework detailed in Section~\ref{sec:rhp} schedules tasks for each agent up to a pre-defined time horizon and leverages a heuristic to estimate the cost to go for each schedule. As such, this framework does not suffer from the exponential complexity caused by the scheduling of the tasks and is robust to changes in the environment. Moreover, a specialized version of the framework for the SaR application is presented in Section~\ref{sec:useCase} along with a simulation environment for an abstracted SaR scenario. To validate the effectiveness of the proposed approach, Section~\ref{sec:experiments} presents extensive experimentation comparing the proposed Receding Horizon Planning based approach to a greedy scheduling scheme.

\section{Literature Review}
\label{sec:litReview}
As mentioned in Section~\ref{sec:intro}, the topic of coordination for multi-agent teams in SaR scenarios following natural disasters falls under the umbrella of Multi-Robot Task Allocation (MRTA). A comprehensive taxonomy of existing methods for MRTA can be found in \cite{gerkey2004formal}. As highlighted in \cite{khamis2015multi}, most existing approaches can be categorized as decentralized or centralized approaches.

There are two main advantages to decentralized approaches: their robustness to varying team sizes and communication failures, and their scalability with respect to the size of the agent fleet thanks to the computational burden being shared amongst the agents. Moreover, market-based approaches such as \cite{coalitions,dias2006market,vig2006market}, attempt to combine the benefits of centralized and decentralized methods by having the computational burden shared between a central entity and the remainder of the fleet. For example, in \cite{coalitions}, the authors suggest a protocol where agents communicate their respective capabilities and use this information to form the coalitions in a decentralized fashion. As such, these methods are able to generate better solutions than fully decentralized approaches while maintaining a certain level of scalability. However, since SaR scenarios typically involve a bounded number of agents and little computational constraints, scalability with respect to the size of the team is of no concern; thus making fully-centralized approaches more suitable for this application.

Moreover, the topic of task allocation specifically pertaining to SaR scenarios is well-studied, most notably, by the participants of RoboCup SaR Agent Simulation competition as highlighted in \cite{Sheh2016}. The competition setup is as follows. A heterogeneous team is to be deployed to extinguish fires and rescue victims. Specifically, there are three types of agents: ambulances which rescue victims, fire brigades which extinguish fires and police units which remove the road blockades enabling the two other type of agents to reach their desired targets faster. The state-of-the-art task allocation strategy utilized by winning teams such as MRL in the competition is K-Means clustering of the Fires/Victims followed by a cluster to agent assignment using the Hungarian Algorithm which runs in polynomial time. We refer the reader to \cite{mrl19} for details.

Inspired by the RoboCup competition, the motivation behind the development of the new simulation environment presented in Section~\ref{sec:useCase} is two-fold. First, the proposed scenario can be seen as a generalized version of the competition's scenario. Specifically, instead of having a fixed number of types of agents each associated with a single class of tasks (e.g. police units only capable of removing road blockades), through characterizing each agent through the capabilities it exhibits, the proposed simulation framework allows for the modular addition of agent types and the collaboration of a heterogeneous sub-team of agents in achieving a single task. For example, given any two agents and their potentially different capacities to transport water, in the proposed scenario, they can indeed collaborate to extinguish a target fire. Moreover, since not all SaR scenarios are identical in nature, the proposed simulation environment frames the problem as a dynamic set of pick and place tasks where victims and resources are to be delivered to target locations.

Based on this abstracted view of SaR problems, the Receding Horizon Planning framework presented in this paper aims to leverage the strengths of centralized scheduling approaches while keeping the problem size tractable and remaining robust to changes occurring in the environment. The robustness property is obtained through the repeated generation of schedules at fixed time intervals; whereas tractability of the problem size is obtained through only scheduling tasks up to a certain time horizon and leveraging a load-balancing Linear Program as a heuristic to estimate the cost-to-go. Consequently solutions produced by the proposed framework are not guaranteed global optimality since the schedules are of finite horizon. However, we present extensive empirical evidence demonstrating the effectiveness of the proposed approach in Section~\ref{sec:experiments}. In the next section, the proposed Receding Horizon Planning framework is introduced and presented in detail.

\section{The Receding Horizon Planner}
\label{sec:rhp}

In this section, we present an extended version of the Receding Horizon Planner (RHP), first presented in \cite{rhp}, aimed at solving the Single-Task robots, Single-Robot tasks, Time-extended Assignment (ST-SR-TA) problem, as defined in \cite{gerkey2004formal}, for heterogeneous teams of agents. This problem class involves building a schedule of tasks for each agent that minimizes a given cost function and is strongly $\mathcal{N}\mathcal{P}$-hard as highlighted in \cite{brucker1999scheduling}.

The brute-force approach for solving this class of problems is to enumerate all possible schedules and choose the one with the smallest associated cost. In its simplest form, the process of generating all possible schedules is done through iteratively assigning one of the remaining tasks to each agent's schedule. As such, a set of partial schedules (i.e. schedules that do not include all tasks) will be generated to which the process is applied again. Similarly to Branch and Bound (\cite{lawler1966branch}), this process can be depicted as a tree-search where each partial schedule is associated to a node $\mathcal N_i$, and partial schedules generated through subsequent assignments are depicted as children of that node. However, given the set of tasks $\mathcal T$ and the set of agents $\mathcal A$, the number of possible schedules grows with $\bigO(|\mathcal{A}|^{|\mathcal{T}|})$, rendering the enumeration of all schedules intractable. A scalable alternative is the greedy approach, which solely considers the next best option given the current partial schedule. However, this approach suffers from a reduced performance of the overall system in terms of optimality.

The RHP is a task allocation scheme inspired by Model Predictive Control (MPC) which, at fixed time intervals, computes the optimal schedule for a limited number of tasks and leverages a heuristic to estimate the cost of executing the remaining tasks. Thus, the size of the optimization problem remains constant with respect to the number of tasks and can be adjusted to the computational resources available. With regard to this variable look-ahead time, the receding horizon approach is a superset of the greedy algorithms (zero look-ahead) and the full-blown optimization (infinite look-ahead).

Similarly to MPC, the number of assignments planned by the RHP is larger than the number of assignments that are executed. This creates an overlap between the consecutive optimization cycles, which reduces the loss in optimality due to the neglected future operations. The cyclic nature of this scheme allows for the incorporation of the current system state into the optimization. This feedback loop –~as in classic control~– provides robustness against disturbances and model deviations.

\begin{figure}
    \centering
    \includegraphics[width=0.7\columnwidth]{figures/DecisionTree.pdf}
    \caption{Example decision tree generated by the Receding Horizon Planner.}
    \label{fig:DecisionTree}
\end{figure}

As introduced in \cite{rhp}, The RHP utilizes a branch-and-bound method to generate the optimal schedule up to the desired time-horizon. As illustrated by Figure~\ref{fig:DecisionTree}, each node $\mathcal{N}_i$ in the search tree corresponds to a partial schedule and the addition of a new assignment of a task to an agent creates a new node. To enable an efficient exploration of the tree, the cost $J(\cdot)$  at node  $\mathcal{N}_i$ is decomposed into an accumulated cost value $g(\cdot)$ and a remaining cost estimation $h(\cdot)$
\begin{equation}
J(\mathcal N_i) = g(\mathcal N_i) + h(\mathcal N_i).
\label{eq:CostDecomposition}
\end{equation}
The accumulated cost evaluates the already assigned operations and is a measure of the consumption of resources. Since each new task assignment increases the consumption of resources, the accumulated cost increases monotonically
\begin{equation}
    g(\mathcal N_j) \geq g(\mathcal N_i),
    \label{eq:AccumulatedCost}
\end{equation}
where $\mathcal{N}_j \in Children(\mathcal{N}_i)$. The second summand in~\eqref{eq:CostDecomposition} is a lower-bound estimate of the remaining efforts to reach the overall goal. As each new task assignment reduces the outstanding efforts the cost of the remaining tasks must decay
\begin{equation}
    h(\mathcal N_j) \leq h(\mathcal N_i).
    \label{eq:RemainingCost}
\end{equation}
Moreover, the estimate $h$ must provide a lower bound of the true remaining cost at each step, therefore satisfying
\begin{equation}
    g(\mathcal N_j) - g(\mathcal N_i) \geq h(\mathcal N_i) - h(\mathcal N_j).
    \label{eq:LowerBoundCondition}
\end{equation}
We will show later, how such a lower bound estimation can be obtained by relaxation of the integer constraints. If $h$ is chosen such that~\eqref{eq:LowerBoundCondition} holds we can conclude that
\begin{equation}
        J(\mathcal N_j) \geq J(\mathcal N_i),
    \label{eq:TotalCostEstimate}
\end{equation}
the cost of each node increases as the tree grows. This allows to stop the further exploration of a branch if at any time $J(\mathcal N_i) \geq J_{\text{opt}}$ (i.e. the cost of node $\mathcal N_i$ is larger than the cost of a known solution). This strategy is guaranteed to find the optimal solution on the tree.

To eliminate symmetries and thus reduce the number of nodes to be explored in a tree, only one resource (agent) is chosen for the set of offspring-nodes that are generated from any node in the tree. The agent is chosen as
\begin{equation}
a^{\ast} = \argmin_{t \in \mathcal T , a \in \mathcal A} \ (y_{ta} + T_{ta}).
\label{eq:NextAgent}
\end{equation}
The indices $t$ and $a$ tally the available tasks $\mathcal T$ and agents $\mathcal A$ respectively. $T_{ta}$ is the duration of task $t$ if performed by agent $a$ and $y_{ta}$ is the potential start time for agent $a$ on task $t$. In other words, out of all agents, the RHP chooses the one with the earliest completion-time of all tasks. Once the agent is decided upon, all potential tasks that satisfy
\begin{equation}
y_{ta}  = \min_{t \in \mathcal T } (y_{ta^{\ast}} + T_{ta^{\ast}}),
\label{eq:NextStates}
\end{equation}
are considered as next nodes. That is, we include all tasks that can be started, before the earliest task can be finished.  This again reduces the search space in the tree exploration without affecting the optimality of the solution.

The main contribution of this paper is to extend the range of applications of the RHP to scenarios in which multiple agents, each exhibiting different capabilities, are needed to complete a single task or vice-versa. The wildfires in our SaR scenario represent the former type, where the combined effort of multiple agents is needed to extinguish a fire. The rescue operations represent the second type, where a single agent can carry multiple survivors.

The decision space in the tree search inherently includes the various agent capabilities and teaming scenarios, if this is properly described in the set of dispatching rules for each agent. These dispatching rules describe which possible tasks an agent can do, given its current location and occupation, and to what amount the agent can contribute to the overall task, i.e. its capacity. In the context of the SaR scenario, the capacity corresponds to the mundane load capacity for water or victims of each agent.

While the implementation of these dispatching rules is straightforward, the challenge for the tree search is again limiting the search space. To avoid the exploration of all possible combinations of capacities to complete a task, we make use of the heuristic $h(\cdot)$ to guide the tree search. To achieve this, we include a high level load balancing into the heuristic, which breaks the required effort for one large task down into a set of sub-tasks that can be handled by individual agents.

The general approach is to first group all agents with equivalent capabilities into different classes and then minimize the latest finishing time among all agents under the following constraints
\begin{enumerate}[(i)]
\item \label{constraints:first} The effort for each task is distributed among the different agent classes
\item \label{constraints:second}  The efforts from~(\ref{constraints:first})  for each class are distributed among the individual member agents of this class
\end{enumerate}
The detailed application to the SaR use-case is described below.

\begin{figure*}[tb]
    \centering
    \begin{minipage}{0.33\textwidth}
    \includegraphics[width=\textwidth]{figures/scenarioInitSmall.jpg}
    \end{minipage}~%
    \begin{minipage}{0.33\textwidth}
    \includegraphics[width=\textwidth]{figures/scenarioMidSmall.jpg}
    \end{minipage}~%
    \begin{minipage}{0.33\textwidth}
    \includegraphics[width=\textwidth]{figures/scenarioEndSmall.jpg}
    \end{minipage}~%
    \caption{An example scenario in the Search and Rescue Simulator. The heterogeneous team of agents is initially positioned in the starting zone (left). The orange robots denote the AGVs and the red, green and pink circles denote the ground units, drones and helicopters respectively. Hidden and identified victims are depicted as pink and green squares respectively. Fires, which can grow and spread, are depicted as red circles of radii proportional to their respective strengths. The middle and right figure depict the scenario mid and post execution respectively.}
    \label{fig:sim}
\end{figure*}

\section{Use Case: Search And Rescue}
\label{sec:useCase}

\subsection{Problem Setup}
Inspired by the RoboCup Search and Rescue Agent competition, we setup the SaR scenario as follows. A fire breaks out in a forest near a city and starts spreading. Victims are hidden within the city and forest themselves. The objective is to locate and save all victims and extinguish all fires simultaneously using a heterogeneous team of agents as quickly as possible. Note that the fires also grow and spread, therefore solely focusing on rescuing the victims is commonly a sub-optimal strategy. Moreover, the victims are initially not visible to the agents. Therefore, exploring the map is also of paramount importance.

 Each fire is represented as a circle of radius proportional to its health and requires a specific amount of water that is also proportional to its health. Upon reaching $100\%$ health, a fire then spread to a nearby territory which is illustrated through the creation of another circle in that location. In order to keep the problem setup as general as possible, we frame the scenario as a pick and place problem where agents need to repeatedly go back to the base and deliver water/victims to the fires/hospital. Moreover, we solely consider each agent's capabilities when considering it for a given task. As such, we allow for the modular addition or removal of agents from the setup as discussed in the next subsection.

\subsection{The Heterogeneous Team}
Teams deployed for SaR are typically heterogeneous due to the variety of the tasks at hand and terrains to be navigated. Therefore, we chose to create teams composed of $4$ types of agents: Ground Units, Helicopters, Drones and Autonomous Ground Vehicles (AGVs). It is important to note that our approach is agnostic to the specific types of agents and the number existing types. In fact, the proposed approach solely considers the number of agents of each type, and each type's capabilities (as defined in \cite{gerkey2004formal}). We restrict the agent capabilities we consider to two categories. The first type deals with the mobility of the agents (e.g. what is the velocity of the agent when navigating in the forest?), and the second type considers what action the agent can perform once at the desired location (e.g. can the agent "pickup" a victim?). A tabulation of the capabilities of the $4$ types of agents is presented in Table~\ref{tab:agentCaps}. As such, each agent type's specific capabilities can be accounted for explicitly by the proposed framework in the process of generating feasible schedules. Moreover, it is worth noting that additional capabilities can be modularly added since the framework solely considers the capabilities required by each task in the process of generating new nodes. In the next subsection, we present the load-balancing Linear Program used to estimate the cost-to-go.

\subsection{Estimation of the cost-to-go}
We formulate the cost-to-go required for the RHP as a Linear Program (LP), which added to the accumulated cost provides a lower-bound on the total cost of each node. The objective of the LP is to minimize the makespan $s$, which is the time of completion of the last task.

In order to compute the makespan, one must be able to estimate the time taken to complete each task (e.g. rescuing a victim). This is a non-trivial problem, since the time of completion of a task by an agent is dependant upon the previous assignment of the agent. This difficulty also arises in the travelling salesman problem, where the time to travel to a given city depends on the last destination of the salesman. In order to overcome this difficulty, we assume that the distances between the victims/fires are negligible compared to their distances to the base. Therefore by lower-bounding all distances between the targets and the base, we can obtain a ``tight'' lower-bound on the amount of time a given agent takes to complete a trip to any target $t_{Type(a)}$ given the agent's velocity.

The decision variable for the LP are
\begin{itemize}
    \item the number of assignments of task $t$ to the agents of class $c$, denoted by $n_{c,t}$
    \item the number of assignments of task $t$ to the individual agent $j$ denoted by $m_{j,t}$
    \item the total makespan denoted by $s$
\end{itemize}
With these constraints stacked into a vector
\begin{equation}
    x^{\text{T}} = [n_{c,t}\,,\,m_{j,t}\,,\, s]
\end{equation}
the LP is formulated as follows
\begin{subequations} \label{eq:heuristicLP}
\begin{align}
\min_{x} ~& \; \; [0\, \ldots\, 0\, 1] x  \label{eq:lp:a}\\
\subjto~~& \smashoperator{\sum_{c \in \mathcal C}} C_{c}^{(t)}n_{c,t} \geq R_{t} \: \forall t \in \mathcal T \label{eq:lp:b} \\
& y_j + \smashoperator{\sum_{t \in \mathcal T}} T_t^{(c)}m_{j,t} \leq s \:\: \forall j \in \mathcal A \label{eq:lp:c}\\
& \sum_{j \in \mathcal{A}}B_j^{(c)}  m_{j,t} = n_{c,t} \:\: \forall c,t \in \mathcal C  \otimes \mathcal T \label{eq:lp:d},
\end{align}
\noeqref{eq:lp:a}\noeqref{eq:lp:b}\noeqref{eq:lp:c}
\end{subequations}
where $\mathcal T$, $\mathcal C$ and $\mathcal A$ denote the set of all tasks, agent types and individual agents respectively. The program aims to minimize the makespan (i.e. the time of completion of the last task). Moreover, constraint \eqref{eq:lp:b} ensures that the required effort for accomplishing task $t$ denoted by $R_{t}$ is matched by the agent fleet. The effort provided by the agent fleet is computed through summing the effort provided by each agent type. The capacity of agent type $c$ for task $t$ is denoted by $C_{c}^{(t)}$.  Additionally, the makespan is computed in constraint \eqref{eq:lp:c}, where $y_j$ denotes the completion of time of agent $j$'s current schedule and the execution time of task $t$ for class $c$ is denoted by $T_{t}^{(c)}$. Lastly, constraint \eqref{eq:lp:d} ensures that the sum of contributions of each individual agent of a given type matches the total contribution of that type, where $B_{j}^{(c)}$ indicates if agent $j$ is of type $c$. In the next section, we present experimental results demonstrating the effectiveness of the SaR specialized Receding Horizon Planner presented in this section compared to a greedy scheduling approach.

\begin{table}
\centering
 \begin{tabular}{|c||c c c c||}
 \hline
 Types & G. Unit & Heli. & Drone & AGV \\ [0.5ex]
 \hline\hline
 Water Cap. & 1 & 5 & 0 & 2 \\
 \hline
 Rescue Cap. & 1 & 4 & 0 & 4 \\
 \hline
 Move Forest & 0.1 & 0.5 & 0.40 & 0.20\\
 \hline
 Move City & 0.1 & 0.5 & 0.40 & 0.00 \\
 \hline
\end{tabular}
\vspace{3mm}
\caption{The capabilities of the $4$ types of agents: Ground Units, Helicopters, Drones and AGVs. Note for example how AGVs are unable to navigate through the forest or how the capacity for carrying water units or victims vary depending on the type of agent.}
\label{tab:agentCaps}
\end{table}

\section{Experiments}
\label{sec:experiments}
In this section, we present empirical results validating the use of a receding horizon through simulations and experiments on the Robotarium, a remotely accessible swarm accessible testbed (\cite{robotarium}). The setup of the experiments is as follows. The team is composed of 6 ground units, 3 helicopter units, 14 drones and 3 AGVs. Moreover, in each experiment $10$ initially hidden victims are randomly placed in the forest and city along with $3$ fires that were also initialized at random positions. However, the number of total fires generated in a given experiment may differ depending on the planner being used. This is due to the fact that the fires grow and spread over time and therefore require even more resources to extinguish. This is an accurate depiction of many real-world SaR scenarios and serves to emphasize that the time of completion of tasks is an important measure in such scenarios.

 The path-action planning algorithm details the execution of the generated schedules and is implemented as follows. First, it generates each agent's path to its corresponding task, then ensures that the agent takes the corresponding required action for the task if feasible. For example, once an agent tasked with rescuing a victim reaches its location, it will take the "pick-up" action if and only if the agent's maximum capacity for the number of victims is not already reached. Since path planning is not the focus of this work, we assume all agents except the AGVs possess single-integrator dynamics and use  proportional controllers to guide them to their targets. However, since the AGVs are presented as robots (GRITSBot X) in the Robotarium experiments and can indeed collide, we implemented multi-agent A* to generate way-points the agents can follow to their targets. Moreover, since we do not explicitly check for collision-avoidance in the trajectories between way-points, we also utilize Control Barrier Functions (CBFs) to instantaneously ensure collision-avoidance at all times as described in \cite{AmesBarriers, ames2014}. This is achieved through solving a Quadratic Program at each point in time that generates a minimally altered trajectory for the agents relative to their nominal trajectory to ensure collision-avoidance.

A run of $20$ simulated experiments with randomized initial conditions were run to compare several planning depths of the RHP and a greedy scheduling approach. The greedy scheduling scheme used for bench-marking the proposed RHP is a one-step look-ahead planning approach. Specifically, at each scheduling iteration, each idle agent is assigned to the task that it is closest to. The mean, median and variance of the makespans of each of the scheduling approaches over all experiments are presented in Table~\ref{tab:results}. As shown in the table, the mean makespan decreases significantly when the RHP is used. However, the rate of improvement decreases at higher planning depths, which highlights the trade-off between computing time and solution quality.

Moreover, to demonstrate the applicability of the RHP onto real systems, $10$ experiments were conducted on the Robotarium testbed comparing the RHP using a planning depth of $10$ to the greedy scheduler. The Robotarium's robots (GRITSBot X) were used instead of the three simulated AGVs as depicted in Figure~\ref{fig:robotariumExp}. To obtain the desired frequency of operation on the Robotarium ($\sim100$~Hz), a separate computing node was used to run the scheduling algorithms, and the schedules were transmitted to the agents using a publisher-subscriber protocol. The results of the experiments are presented in Table~\ref{tab:resultsR}. Indeed, the use of the RHP reduces the makespan, thus validating the applicability of the proposed scheduling approach.

\begin{figure}[t]
\centering
 \includegraphics[width=0.40\textwidth]{figures/roboExp.png}
 \caption{Search and Rescue experiment on the Robotarium using 3 GRITSBot X as the AGVs. The description of all entities are as in Figure ~\ref{fig:sim}.}
 \label{fig:robotariumExp}
\end{figure}

\begin{table}
\centering
 \begin{tabular}{|c||c c c||}
 \hline
 Types & Mean & Median & $\sigma$ \\ [0.5ex]
 \hline\hline
 Greedy & 78.78 & 53.03 & 10.63  \\
 \hline
 RHP10 & 42.54 & 40.34  & 9.54  \\
 \hline
  RHP15 & 40.09 & 38.56 &  7.95 \\
  \hline
  RHP20 & 39.87 & 36.84 &  9.39 \\
 \hline
\end{tabular}
\vspace{3mm}
\caption{Results over $20$ simulations comparing the mean, median and standard deviation of the makespans of the Receding Horizon Planner with different planning depths and greedy scheduling scheme.}
\label{tab:results}
\end{table}

\begin{table}
\centering
 \begin{tabular}{|c||c c c||}
 \hline
 Types & Mean & Median & $\sigma$ \\ [0.5ex]
 \hline\hline
 Greedy & 46.70  &  45.06  & 9.64 \\
 \hline
 RHP10 & 39.00 &  36.72 & 9.09 \\
 \hline
\end{tabular}
\vspace{3mm}
\caption{Results over $10$ Robotarium experiments comparing the mean, median and variance of the makespans of the Receding Horizon Planner with planning depth $10$ and the greedy scheduling scheme.}
\label{tab:resultsR}
\end{table}

\section{Conclusion}
\label{sec:conclusion}
This paper introduces a task allocation framework capable of scheduling tasks for heterogeneous teams of agents in a manner that is tractable and robust to changes in the environment and in the agent fleet. This was achieved through repeatedly scheduling solely up to a fixed horizon and leveraging a load-balancing Linear Program for the estimation of the cost to go. Moreover, a simulation framework for an abstracted Search and Rescue scenario inspired by the RoboCup Search and Rescue Agent Simulation competition was presented along with a specialized formulation of the Receding Horizon Planning approach. Experimental results showcase the efficacy of the proposed scheduling method in extensive multi-agent simulations and experiments on the Robotarium.


\section{Introduction}

There is a trend to forecast COVID-19 using mathematical models for the probability of moving between states from susceptible to infected, and then to a recovered state
or death (SIR models). This approach is however very sensitive to starting assumptions and tend to overestimated the virus reproductive rate. One key point that these models miss is the  individual behavioral responses and government-mandated policies that  can dramatically influence
the course of the epidemic. In Wuhan, for instance, strict social distancing was instituted on January 23rd, 2020, and by March 15th new infections were close to zero.  Taking into account this observation,
\citet{covid2020forecasting} have proposed  a statistical approach to model a empirical cumulative population death rate. However, modelling observed cumulative death numbers has the inherent problem of yielding a  highly correlated data which, if it is not taken into account, could cause misleading  inference results.

Instead of considering a mathematical model such as SIR, that rely heavily on parameters assumptions, we will directly work as in \citet{covid2020forecasting} or \citep[]{zhou2020forecasting} with  empirical data. To this end,  we propose to model the daily number of deaths using a Poisson distribution with a rate parameter that is proportional to a  Skew Normal density. Using a Bayesian approach and a prior epidemic China COVID-19 history, we forecast the total number of deaths in Peru for the next seventy days.

\section{Methods}

\subsection{Model}

Let  $Y_{t_1}, Y_{t_2}, \ldots, Y_{t_n}$ be the number of COVID-19 deaths at times $t_1, t_2, \ldots, t_n$, where time is measured from the first reported death due to COVID-19, and let us suppose that the death rate will hit a platoon due to the government intervention. We propose then the model

\begin{eqnarray}
Y(t_i) &\sim& Poisson(\lambda(t_i)), \ i =1,2,\ldots,n
\label{eq1}
\end{eqnarray}
with  death rate
\begin{eqnarray}
\lambda(t_i) &=& p g(t_i \mid  \alpha, \beta, \eta)K,
\label{eq2}
\end{eqnarray}

\noindent where $g(t_i \mid  \alpha, \beta, \eta)$ denotes the density function of a Skew Normal distribution with location, scale, and shape parameters,  $\alpha$, $\beta$, and $\eta$, respectively; $p$ is a maximum asymptotic level parameter  and $K$  is the population size.  The choice of the  Skew Normal distribution  is motivated by its  flexibility on the tails and their asymmetry, which  can force, as one could  expect, a rapidly increase rate at the first stages of the pandemic and a slower decrease of this rate at the last stages. A  negative binomial distribution can also be considered for the deaths numbers, but we found empirically a better fit with the Poisson model.

Our model differs of the approach taken by  \citet{covid2020forecasting}, who directly models the cumulative death rate $\Lambda(t) = \int_0^t \lambda(s)ds$ with a term proportional to a symmetric cumulative normal distribution. This difference is not only found in the formulation, but also in the estimation procedure. While these authors incorporated first the Chinese data--more concretely the time from
when the initial death rate exceeds 1e-15 to the implementation of social distancing-- into their model throughout a location-specific inflection point parameter or a maximum death rate and then used a sort of credibility model between short-range and long-rate variants, we followed a Bayesian approach that incorporates the Chinese data as a prior distribution. The posterior predictive distribution, which is the goal of this model, is then a dynamic object that can be updated with new data   about the number of deaths or any other reliable information  that may be incorporated as covariates in the model.

Apart from the total number of deaths, three main quantities of interest  can be easily derived from our model. First,  a time to threshold death rate, which  provides  the time after which only 0.01 $\%$ of deaths will be observed on the population. This will be defined as the 0.99 quantile of the $g$ distribution. Another characteristic of interest is the inflection point, defined as the time at which the death rate reaches its maximum level.

\subsection{Priors}

Since China was the first country to have experienced  a drastic drop in infections and deaths, we are proposing to incorporate this data, into our Peruvian death rate predictions, through a prior distribution. Figure \ref{fig:lfig0} shows the empirical distribution of number of reported deaths in China.  One can notice that the reported numbers on February 13 and 14 (red points) were 254 and 13 deaths, respectively. There is some controversy\footnote{https://www.bbc.com/news/world-51495484}  about the reported numbers from China on these days, the reason why we are considering  the average number for these days.

Figure \ref{fig:lfig1} shows the observed and predicted death rates (A) and the daily number of deaths (B) in China. The predicted rates and their associated 95\% prediction credible intervals were obtained, under a Bayesian approach, by considering a non-informative prior and the Chinese official death reports. Table \ref{tab:tab1} summarizes the information given in Figure \ref{fig:lfig1}. This information will be considered as a prior for modelling the number of deaths in Per\'u with the exception of the $p$  parameter, where a weakly informative prior, $N(0,10^2)$, will be considered for $\log{(p)}$.

\begin{figure}[!htb]
\begin{center}	\includegraphics[width=15.0cm,height=10.0cm]{Figure0.pdf}
\caption{Empirical distribution of number of deaths reported in China. Points in red represents the information on February 13 and 14, respectively.}\label{fig:lfig0}
\end{center}
\end{figure}

\begin{table}[ht]
\caption{China: Mean, standard deviation (SD) and quantiles 0.25,0.5 and 0.975 of the posterior distribution}
\begin{center}
\begin{tabular}{lccccc}
  \hline
 & Mean & SD & 2.5\% & 50\% & 97.5\% \\
 \hline
 $\beta$ & 16.52 & 0.37 & 15.81 & 16.52 & 17.24 \\
  $\log{(\alpha)}$ & 2.91 & 0.02 & 2.87 & 2.91 & 2.95 \\
  $\eta$ & 2.34 & 0.14 & 2.07 & 2.34 & 2.64 \\
  $\log{(p)}$ & 0.87 & 0.02 & 0.84 & 0.87 & 0.91 \\
   \hline
\end{tabular}
\end{center}
\label{tab:tab1}
\end{table}

\begin{figure}[!htb]
\begin{center}	\includegraphics[width=15.0cm,height=10.0cm]{Figure1.pdf}
\caption{Figure A: Estimated cumulative mortality rate for China and its associated 95\% predictive interval where points represent the empirical data. Figure B: Estimated number of infected for China and its associated 95\% predictive interval where points represent the empirical data.}\label{fig:lfig1}
\end{center}
\end{figure}

Initial values for the estimation process will be sampled from the prior distribution that was constructed with the Chinese data.

\subsection{Data}

The Chinese data to be used in this work was obtained from  the European Centre for Disease Prevention and Control, institution  that  daily publishes statistics on the COVID-19 pandemic. The daily death  reports in Per\'u, on the other hand,  were obtained from the local authorities.

\subsection{Inference}

Taking into account the observed likelihood function, easily derived from (\ref{eq1}) and (\ref{eq2}), the posterior distribution of the vector parameter $(p,\alpha,\beta,\eta)$ has a complex expression. To obtain samples from it, we will make use of Hamiltonian MCMC methods as the ones implemented in the R package RStan \citep{rstan}. This package implements an adaptive Hamiltonian Monte Carlo algorithm (also known as HMC) using a No-U-Turn sampler (NUTS) for the stepsize parameter in order to generate efficient transitions to the posterior distribution \citep*{carpenter2017stan}.

\section{Results: Forecasting for Per\'u}

We now briefly describe our predictions for the spread of COVID-19 in Per\'u. The time here will be understood to be measured in days after the first reported death in the country.

\subsection{Number of deaths}

Figure \ref{fig:lfig2}.A shows the estimated cumulative number of deaths for Per\'u and its associated 95\% predictive credibility interval. In particular, the expected total number of deaths is 141.8 (95\%CI: 99 - 192) and 339.7 (95\% CI: 239 - 462) after 20 and 30 days, respectively.

In addition, Figure \ref{fig:lfig2}.B shows the expected number of deaths per day and its associated 95\% predictive credibility interval. In particular, the expected number of deaths is 17.5 (95\%CI: 9 - 28) and 19.3 (95\% CI: 10 - 31) at 20 and 30 days, respectively. Overall, the total number of deaths is expected to be 611.6 (95\%CI: 604.3 - 833.7) persons.

\begin{figure}[!htb]
	\begin{center}	\includegraphics[width=15.0cm,height=10.0cm]{Figure2.pdf}
		\caption{A: Estimated cumulative number of deaths for Per\'u and its associated 95\% predictive interval. B: Estimated number of deaths for Per\'u and its associated 95\% predictive  interval.}\label{fig:lfig2}
	\end{center}
\end{figure}

\subsection{Time to threshold death rate}

The model estimates that the number of days, since the first death, to achieve the threshold death rate  will be 63.8 (95\%CI: 52.9 - 65.84) days. At this time, the 99\% of expected deaths will have been observed.

\subsection{Inflection point}

The estimated inflection point is 25.9 (95\%CI: 25.1 - 26.8), which means that we expect to spend approximately 26 days before the COVID-19 death rate starts to decline after the first reported death case.

\section{Sensitivity analysis}

Three different scenarios were considered, each keeping the weakly informative prior log-normal distribution for $p$. In all cases,  the posterior mean for the fitted Chinese model  was taken as the mean prior for the modelling of data from Per\'u. The first scenario (I) takes also the posterior Chinese variances as  the prior Peruvian variances. The second scenario (II) induces flexibility in the prior distribution by increasing their corresponding  standard deviations by a factor of five. Finally, the third scenario (III)  increases the standard deviations by a factor of ten.

Table \ref{tab:tab2} and Figure \ref{fig:lfig4} shows that the point estimates are not heavily affected, but the precision pays the price of not having yet enough data from Per\'u.

\begin{table}[ht]
\caption{Mean (95\% credible interval) of the posterior distribution under three scenarios. Scenario I: Prior distribution using the Chinese information and $p$ free, Scenario II: Same as Scenario I, but multiplying by 5 the standard deviations, and Scenario III: Same as Scenario I, but multiplying by 10 the standard deviations.}
\begin{center}
\begin{tabular}{cccc}
Scenario & Time to threshold & Inflection point & Total number of deaths \\
\hline
I        & 63.8 (61.9 - 65.8) & 26.0 (25.1 - 26.8)
& 611.6 (437.2 - 833.7)\\
II & 64.7 (55.4 - 74.9) & 26.2 (22.0 - 30.3)
& 568.7 (299.4 - 1037.5) \\
III & 66.6 (48.9 - 88.5) & 26.6 (19.0 - 34.9)
& 629.6 (225.8 - 1588.4 )\\
\hline
\end{tabular}
\end{center}
\label{tab:tab2}
\end{table}

\begin{figure}[!htb]
	\begin{center}		\includegraphics[width=15.0cm,height=10.0cm]{Figure4.pdf}
		\caption{Predictive distribution under three scenarios. Scenario I: Prior distribution using Chinese information and $p$ free (black). Scenario II: Prior standard deviations of Scenario A were multiply by 5 (red). Scenario III: Prior standard deviations of Scenario A were multiply by 10 (green).}\label{fig:lfig4}
	\end{center}
\end{figure}

\section{Discussion}

Considering the information on daily number of deaths from China and from Per\'u, our estimation of total number of deaths will be 611.6 (95\% CI: 437.2 - 833.7) and 99\% of those deaths will occur 63.8 (95\% CI: 52.9 - 65.8) days after the first death reported case.

The present study has some limitations. Although Per\'u has not followed all the measures taken  in China, we assumed that Peruvian interventions will have similar effects as the ones observed in that country. However, we expect our model is flexible enough to be driven by the Peruvian data. Furthermore, the model will increase its precision as more data from Per\'u becomes available.

Several extensions are possible for the model. For instance, covariates as daily social mobility indicators or country age distributions can be included in the model, in particular on the $p$ parameter that measures the total number of deaths.

We expect our model can be useful to guide some policies that need to be taken by the Peruvian government in order to overcome the COVID-19 pandemic. For example, the proposed model can be useful to measure the impact of the COVID-19 pandemic on the Peruvian health system.

\section{Acknowledgement}

We would like to thanks Rodrigo Carrillo Larco for providing us the detailed information of Per\'u based on reports by the Peruvian Ministry of Health and Dr. Daniel Manrique-Vallier for earlier discussions and helpful suggestions on the topic.


\section{Introduction}

Throughout $G$ is a graph of order $n$ with adjacency matrix $A$.
Let $I$, $J$ and $D$ be the identity matrix, the all-ones matrix and the diagonal matrix with the degrees of $G$ on the diagonal,
respectively.
Any matrix of the form $U=\a A+\b I + \c J + \d D$ with $\a,\b,\c,\d\in\Re$, $\a\neq 0$ (as usuall, $\Re$ is the set of real numbers)
is called a {\em universal adjacency matrix} of $G$.
Several well-studied types of matrices associated to $G$ are a special case of a universal matrix.
For example, if $(\alpha,\beta,\c,\delta)=(1,0,0,0)$, $(-1,0,0,1)$, $(1,0,0,1)$, $(-2,-1,1,0)$ we obtain the adjacency,
the Laplacian, the signless Laplacian and the Seidel matrix, respectively,
and for $(\a,\b,\c,\d)=(-1,-1,1,0)$, $(1,n,-1,-1)$, $(-1,n-2,1,-1)$, $(2,1,-1,0)$
we get the matrices of the complement $\overline{G}$ of $G$ of the mentioned types.

If $G$ is regular of degree $d$, then there is an easy relation between the spectra of the various universal adjacency matrices.
Indeed, if $d=\x_1,\ \x_2,\ \ldots,\ \x_n$ are the eigenvalues of $A$, then
$\a d +\b+\c n+\d d,\ \a \x_2+\b+\d d,\ \ldots,\ \a \x_n+\b+\d d$ are the eigenvalues of $U$.
Let us denote the characteristic polynomial of a matrix $M$, $\det(\x I-M)$, by $\f(M,\x)$.
Then the relation between the spectrum of $A$ and the spectrum of $U$ can be formulated as:
\begin{equation}\label{UA}
\frac{\f(U,\a \x+\b+\d d)}{\a\x-\a d - \c n}
= \a^{n-1} \frac{\f(A,\x)}{\x-d}.
\end{equation}
In this note we consider the case that