\section{INTRODUCTION}

Unlike ground robots, aerial vehicles have an advantage of being able to move freely in three-dimensional environments, so it is emerging as an alternative to overcome the limitations of ground robots. As the technological level of sensor system, flight control, and communication system has advanced, the use of UAVs has been expanded in the last few years, and many attempts have been made to apply autonomous UAVs to actual fields\cite{UAV_application1,UAV_application2,UAV_application3,UAV_application4,UAV_application5, UAV_application6}. It is essential to accurately estimate the position of the UAV in autonomous flight because the autonomous flying UAV determines a flight path based on the current position and performs a given task. In general, UAVs flying outdoors recognize their location based on GNSS. This method has the advantage that the absolute coordinates can be estimated with error from several meters to several centimeters, depending on the GNSS receiver device used. However, if the UAV localizes based only on GNSS, it has a critical limitation that it can be used only in an environment where signals are well received. GNSS signals are easily blocked by buildings and are vulnerable to disturbances. Additionally, UAVs using GNSS signals can be attacked by jamming and spoofing. Therefore, UAV localization that relies only on GNSS signals is not a robust method. In this paper, we propose an image-based localization method to estimate the global position of the UAV without using GNSS.

\subsection{Related Works}
The research for UAV localization without GNSS signals is actively being conducted. Various types of sensors have been tried to accurately estimate the location of drones in large areas. A UAV has a limited payload due to the characteristics of its platform, and if the payload is high or a sensor that consumes a lot of power is used, the flight time is shortened. For this reason, lightweight cameras with low power consumption are widely used. VIO (Visual Inertial Odometry) algorithm has been studied to reduce the error when using only a camera with an IMU\cite{delmerico2018benchmark,vins-mono,svo1,svo2,rovio,okvis}. Since this method estimates the relative position of a UAV after its initialization, there is a limitation that the absolute position cannot be known. If the UAV flies in a path without loop closing, errors will accumulate. Also, if the UAV loses its position, it is difficult to recognize the relative position from the initial position. Various methods have been carried out to solve this problem. A method using a marker such as a QR code\cite{qrcode} or an April Tag\cite{apriltag} has been proposed\cite{marker}. It is a method of calibrating the position of a UAV when a corresponding marker is recognized after pre-positioning a specific marker. This method has the advantage of being able to estimate the position very accurately once the marker is recognized. However, this localization method is not suitable for UAVs flying over large areas. There has been an attempt to develop absolute position estimation method based on the local feature on the flight path\cite{local_feature}, a method that learns landmarks such as buildings located near the flight path in advance and estimates the position of the UAV based on the learned landmark. Nevertheless, this method also has a problem that localization becomes difficult when the UAV deviates from the learned path.

Recently, methods that use free available maps have been proposed to solve this problem\cite{nassar2018deep,goforth2019gps,yol2014vision}. In particular, satellite maps, which store very large areas in the form of orthogonal images, are frequently used to estimate the location of UAVs. Many methods for localization of the UAVs by matching satellite maps with images obtained from the camera which is attached downward from the UAV have been studied\cite{goforth2019gps, nassar2018deep, yol2014vision}. In the case of satellite maps, the image was taken at a certain point of time at a height of hundreds to thousands of kilometers. To match the images taken from the UAV to satellite images, it is necessary to resolve perspective effects, lighting condition differences, seasonal differences, and regional variations. Feature-based matching methods such as SIFT\cite{SIFT}, SURF\cite{surf}, ORB\cite{ORB}, etc. fail to match as shown in Fig. 1.

   \begin{figure}[h]
      \centering
      \includegraphics[scale=0.4]{SIFT.png}
      \caption{This figure shows that SIFT algorithm, which is the widely used conventional point feature matching, fails in matching satellite images(right) and UAV images(left). The area marked in blue is the same area as the image on the left, but the SIFT algorithm fails to find the blue area.}
      \label{figurelabel}
   \end{figure}

In \cite{goforth2019gps}, they used the idea that matching satellite images taken at different times would be a similar task to matching satellite and UAV images. They proposed a localization algorithm by matching satellite maps and UAV images after learning several identical satellite maps taken at various times to consider seasonal, weather, and light disparities. This algorithm works at various environments, from gravel pits to villages. However, the limitation is that the initial position of the UAV should be exactly given and that the algorithm was only tested in a narrow area. In addition, the algorithm was verified only in a very stable situation by capturing the image along a straight path from the static 3D reconstructed map rather than from the actual flight. In \cite{nassar2018deep}, they proposed an SSM (Semantic Shape Matching) method. Some classes, which have specific shapes such as buildings and roads, are segmented from both satellite maps and UAV images and matched by comparing contours of each segmented areas. However, it is also assumed that the GPS position of the initial starting position is known. Additionally, since the shape of the segmented area is used for matching, the algorithm is highly dependent on segmentation accuracy, which makes it difficult to relocalize if the building is lost in complex areas such as the city center. In \cite{yol2014vision}, the measurement of statistical dependency between two signals (mutual information) has been proven to work more robustly than the sum of squared differences (SSD) which is widely used in vision-based matching algorithms. However, because the experiment was only done in a region as small as 130m $\times$ 100m, it is unclear whether the algorithm works well when flying over long distances or when the environment changes.

\subsection{Contributions}
In this paper, we propose an algorithm for estimating the absolute position of the UAV by matching UAV images and a numerical map without using GNSS. Unlike the previous papers, we propose a localization algorithm for UAVs in the absence of exact position coordinates at the starting point of flight. However, it is assumed that the UAV performs the mission within the designated area and exists within the map of the area. For localization, we propose the building ratio feature. Building area is segmented by a segmentation neural network. These segmented images of the buildings are divided into circular areas with various radii, and the ratio of the building area to the circular area is used as a feature. This proposed feature is rotation invariant so it has the advantage of estimating the position candidate group without knowing the position and the orientation. Additionally, it is possible to start localization again at any time even if the position is lost. This feature also works more robustly in terms of the segmentation accuracy than other shape matching algorithms. The training data sets and maps for matching are just gained from freely available online data. The previous satellite map matching-based algorithms are tested only at static environments or in narrow areas. However, we verify the feasibility of the proposed algorithm by real flight data acquired from a UAV.

\subsection{Overview}
In Section {\afterassignment\BigRoman\count255= 2}, a framework of the BRM (Building Ratio Map based) localization is introduced. Then, the proposed feature is defined. Also, a process of the building ratio map generation is described. Building ratio matching algorithm is explained in Section {\afterassignment\BigRoman\count255= 3}. We explain how the number of candidates decreases and that the global position is estimated. In Section {\afterassignment\BigRoman\count255= 4}, we analyze the results of the proposed algorithm and the VIO algorithm. Finally, we conclude with the future works in Section {\afterassignment\BigRoman\count255= 5}.

\section{Proposed Feature and Building Ratio Map}
\subsection{Proposed Algorithm Framework}
   \begin{figure}[h]
      \centering
      \includegraphics[scale=0.41]{flow_chart_new_new.png}
      \caption{Flow of the proposed localization algorithm. IMU and image data are received from a flight controller unit and a camera. From these two inputs, a local odometry is calculated. UAV images are converted to segmented building images. Then, the building ratio matching algorithm runs and the global position of the UAV is estimated.}
      \label{figurelabel}
   \end{figure}
The whole algorithm is shown in Fig. 2. IMU data is obtained from the Flight Control Unit (FCU) and RGB images are received from the monocular camera mounted downward on the UAV. The Visual Inertial Odometry algorithm is operated on two inputs to estimate the relative position from the initial position of the UAV and the relative odometry between each image frame. Images from the camera are passed through the segmentation network at regular intervals and they are converted to segmented building images. The features for matching are calculated from the segmented building images. Each image has not only feature values but also odometry constraints on the previous images. In the pre-processed numerical map, location candidates that meet these conditions are deduced. After the convergence of candidates, the global position of the UAV is estimated.

\subsection{Building Segmented Images}
$w \times h$ $(w > h)$ sized RGB images are received from the mono camera. From this image sequence, some images are extracted at regular intervals to be passed through the pre-trained segmentation network. Through the network, the building area is segmented. To create a rotation invariant feature from the original image, only the $h \times h$ square area that is obtained by cutting off both ends of the original image is used. The \textit{i}-th image created in this way is denoted as $I_{i}$ $(i=1, 2,\dots,m)$.

\subsection{Building Ratio Feature}
   \begin{figure}[h]
      \centering
      \includegraphics[scale=0.5]{feature_new.png}
      \caption{The proposed building ratio feature. Black area is the building segmented area. The ratios of building area in the certain circular areas are used as features. Each circular area is the area covered by all the region until the center.}
      \label{figurelabel}
   \end{figure}

   For the situation where there is no initial information about the location, the following features are proposed for the localization as shown in Fig. 3. The image $I_i$ is divided into \textit{n} circular regions which are denoted by \textit{$S_{k}^{i}$} \textit{$(k=1,\dots,n)$}. \textit{$S_{k}^{i}$} is an area with the radius \textit{$r_{k}$} sharing a center point with $I_i$.
   \begin{equation}
r_{k} = \cfrac{h}{2}\cdot\cfrac{n+1-k}{n}
\end{equation}
   The area of the buildings (black area in Fig. 3) in the area \textit{$S_{k}^{i}$} is defined as \textit{$B_{k}^{i}$}, where the ratio of \textit{$B_{k}^{i}$} to \textit{$S_{k}^{i}$} is used as a feature as follows:

\begin{equation}
f_{k}^{i} = \cfrac{B_{k}^{i}}{S_{k}^{i}}.
\end{equation}
Then, each $I_i$ has \textit{n} features. Since \textit{$f_{k}^{i}$} is a rotation invariant feature, it is possible to estimate the position of the UAV without the initial position based only on the RGB image taken from the UAV and the local odometry estimated based on visual inertial odometry algorithm.

\subsection{Building Ratio Map Generation}
\begin{figure}
\centering
\subfigure[]{\includegraphics[width=0.48\linewidth]{numerical_map_new.png}}
\subfigure[]{\includegraphics[width=0.48\linewidth]{numerical_map_line_new.png}}
\subfigure[]{\includegraphics[width=0.48\linewidth]{numerical_map_binary_new.png}}
\subfigure[]{\includegraphics[width=0.48\linewidth]{feature_map_new.png}}
\caption{The process of building ratio map generation. (a) original numerical map. (b) numerical map which contains only building related layer. (c) binary numerical map where the building area is black part and the other part is white. (d) building ratio map which is the reference feature map.}
\end{figure}

Previously, it is described that the building ratio at predetermined areas of $I_i$ is used as a feature. In addition, we already know the building information of each location through the numerical map. Therefore, when the UAV is located at a certain position within the numerical map, it is possible to know in advance what feature values are to be calculated in the image $I_i$.
In other words, the true value of \textit{$f_{k}^{i}$} is known. We define a building ratio map \textit{$M_{k}$} which consists of the true value of \textit{$f_{k}^{i}$} which is denoted by \textit{$\tilde{f}_{k}$}.
 The example of the original numerical map is shown in Fig. 4 (a). This map includes various artifacts and natural topography. We use building information from this map as shown in Fig. 4 (b), which is composed of the exact position and outline information of buildings. Therefore, this map is used for generating the ground truth. Fig. 4 (c) is a binary map which is generated from Fig. 4 (b). It is a numerical map divided into a building area (black) and a non-building area (white). Finally, the building ratio map \textit{${M}_{k}$} is generated as shown in Fig. 4 (d).
 This building ratio map is generated as follows. First of all, the window size on the numerical map $\tilde{S}_{1}^{i}$ is calculated to be the same scale with $I_i$ as follows:
\begin{equation}
\tilde{S}_{k}^{i} = \pi(\cfrac{n+1-k}{n}\textit{$z_{l}$}\tan{\cfrac{\alpha}{2}})^{2}
\end{equation}
 where \textit{$z_{l}$} is the height of the UAV calculated from visual inertial odometry and $\alpha$ is field of view of the camera. The same number of maps as the number of features used are generated.

\section{Building Ratio Matching Algorithm}
\subsection{Preliminaries}
\begin{figure}[h]
\centering
\subfigure[]{\includegraphics[width=0.32\linewidth]{candidate_1.png}}
\subfigure[]{\includegraphics[width=0.32\linewidth]{candidate_2.png}}
\subfigure[]{\includegraphics[width=0.32\linewidth]{candidate_3.png}}
\caption{Examples of candidates. The gray areas are candidates based on the movement of the UAV. The blue areas are candidates which are calculated by matching with the building ratio map. (a) example of candidates when $\tilde{C}^{i-2}$ = $\emptyset$. (b) $C^{i+1}$ is selected by the constraints of the UAV movement. (c) $\tilde{C}^{i+1}$ is calculated by matching with the building ratio map and $C^{i+2}$ is selected by the constraints of the UAV movement.}
\end{figure}

In the following, we define some notations which will be used in the next section. First of all, ${C}^{i}$ is a set of \textit{i}-th candidates that are selected by the constraints of the UAV movement such as the travelled distance and heading angle. These candidates are depicted in gray color in Fig. 5. $\tilde{C}^{i}$ is a set of \textit{i}-th candidates that are calculated by matching the building ratio map with the feature $f_{k}^{i}$ $(k=1,2,\dots,n)$. These candidates are represented in blue color in Fig. 5. Local odometry of the UAV is calculated during the flight. So two constraints of the UAV movement between two positions, heading and traveled distance of the UAV are made. $\theta^{i}$ is the heading constraint of the UAV at $I_{i}$. It is described in Fig. 5 (b) and (c) and calculated by (4) where $\theta$ is as follows:
\begin{equation}
\theta^{i-1} = \left\{\theta \mid \theta = \arctan({\cfrac{y_{i-1}-y_{i-2}}{x_{i-1}-x_{i-2}}})\right\}
\end{equation}
where $(x_{i-1},y_{i-1})$ and $(x_{i-2},y_{i-2})$ are arbitrary points of the (\textit{i}-1)-th candidates  $\tilde{C}^{i-1}$  and (\textit{i}-2)-th candidates $\tilde{C}^{i-2}$, respectively. If $\tilde{C}^{i-2} = \emptyset$, then $\theta$ can be any degrees between 0 to 360.

\subsection{Candidate Selection}
  \begin{figure}[h]
      \centering
      \includegraphics[scale=0.35]{matching_flowchart.png}
      \caption{
Flow chart of the building ratio matching algorithm. When \textit{i} is received, the feature \textit{$f_{k}^{i}$} is calculated. Then, the previous candidates $\tilde{C}^{i-1}$ are checked for absence. If there is no candidate, new candidate $C^{i}$ is searched throughout the whole building ratio map. If candidates exist, the candidates $C^{i}$ in the current image is determined. After that, the candidates $\tilde{C}^{i}$ for the current location is recalculated by matching with the building ratio map. When the candidates are calculated, the global position of the UAV is estimated.
}
      \label{figurelabel}
   \end{figure}

The flow chart of the proposed building ratio matching algorithm is described in Fig. 6. When the segmented building image $I_i$ is generated by the method described in Section {\afterassignment\BigRoman\count255= 2.\textit{B}}, feature \textit{$f_{k}^{i}$} is calculated. Before estimating the candidates of the current position, the candidates at $I_{i-1}$, which are denoted by $\tilde{C}^{i-1}$, are checked. This is because the previous candidates impose constraints on the areas where the current candidates may exist. According to the presence or absence of the position candidates in the preceding image frame, the position estimation algorithm is divided into two.
If $I_{i}$ is the first image or all the previous candidates were eliminated, there are no candidates of the previous position. In this case, we can assume that all positions on the map can be the initial candidates of $I_{i}$ because $C^{i}$ has no movement constraint. Then
\begin{equation}
C^{i} = \{(x,y) \mid (x,y) \in (X_{M},Y_{M})\}
\end{equation}
where ($X_{M}$, $Y_{M}$) is a set of coordinates on the numerical map. And candidates are calculated on the entire map as follows:
\begin{equation}
\tilde{C}^{i} = \{(x,y) \in C^{i} \mid \sum_{k=1}^{n}|M_{k}(x,y)-f_{k}^{i}| < e_{1}\}
\end{equation}
where $M_{k}(x,y)$ is the true feature value $\tilde{f}_{k}$ at $(x,y)$ position and $e_{1}$ is the feature error threshold.
If the previous candidates exist, $C^{i}$ is estimated from the distance traveled and the heading angle of the UAV as follows:
\begin{equation}
\begin{aligned}
C^{i} & = \{(x,y) \mid |x-(x_{i-1}+d^{i}\cdot\cos{\theta^{i-1}})| < \epsilon \\
& \qquad \qquad,  |y-(y_{i-1}+d^{i}\cdot\sin{\theta^{i-1}})| < \epsilon\}
\end{aligned}
\end{equation}
where $\epsilon$ is the distance error constant and $d^{i}$ is the distance between $I_{i-1}$ and $I_{i}$ which is estimated by VIO. After the constraints-based candidates ${C}^{i}$ are selected, $\tilde{C}^{i}$ are calculated in the same way as (6). Whenever the maximum distance between the average point of $\tilde{C}^{i}$ and an arbitrary point of $\tilde{C}^{i}$ is less than the constant $d_{max}$, we can determine the convergence and the global position can be estimated by averaging $\tilde{C}^{i}$.

\section{Result}
\subsection{Implementation Details}
\textbf{Hardware configuration}
  \begin{figure}[h]
      \centering
      \includegraphics[scale=0.25]{UAV_setup_new.png}
      \caption{The UAV configuration. DJI Matrice 100 is used as the UAV frame, and Pixhawk2 is used for flight control and sensor information collection. Downward images are acquired with the D435i mono-camera. Intel NUC is used for VIO algorithm testing, images storage, and other calculations. RTK GPS is used to obtain ground truth of the UAV flight path.}
      \label{figurelabel}
   \end{figure}

The UAV configuration is shown in Fig. 7. We used DJI Matrice 100 for the drone platform. In selecting the UAV, we considered whether it can hold the wind enough, sustain enough flight time, and be loaded with the necessary sensors on top of it while being not too large or heavy. The monocular camera of Realsense D435i is used. The camera was chosen because it is light enough (72 grams) and the RGB mono camera uses a global shutter. In the case of the camera using the rolling shutter, the screen distortion occurs when the vibration occurs in the actual flight, and the visual inertial odometry algorithm does not work well. We use the Pixhawk2 for flight controller which is capable of obtaining 100Hz IMU sensor data and other sensor information. Intel NUC is used to verify the VIO algorithm result during flight and to save bag files (experimental data storage). The actual flight route is obtained via Here+ RTK GPS. Since RTK GPS has a position accuracy in the range of several cm, the actual flight path is recorded and used to obtain data for comparison with the path estimated by the proposed algorithm.

\textbf{Training network}
In the EO (Earth Observation) task that extracts the information from the satellite map, such as building segmentation and road segmentation, it was observed that U-Net\cite{u-net} showed the best performance\cite{u-net_comparison, u-net_comparison2}, thus U-Net is used for learning. The numerical map mentioned in Section {\afterassignment\BigRoman\count255= 2} is processed and used as labeled data. The size of numerical map is $4,403\times5,555$ pixels with 0.5 meter-per-pixel resolution. The resolution of satellite images is $2,048\times2,048$ pixels. Scale of the satellite images are between 0.23 and 0.47 meter-per-pixel. All the data sets are downloaded from the cities that are not the actual experimental areas and we selected areas that include buildings, fields, rivers, and roads. Numerical maps and satellite maps were resized and cropped to $512\times512$ size and used for training. We trained with a learning rate of $7 \cdot 10^{-7}$ and 124 epoch. The train accuracy is 0.940 and the test accuracy is 0.916.

\textbf{Numerical map and satellite map}
Only freely available maps are used for localization and training data set generation. Satellite images are downloaded from Kakao Map\cite{kakao}. We use numerical map which is freely available at National Geographic Information Institute of Korea\cite{numerical}. The layers from B0014111 to B0014118 at the dxf file of the numerical map are used.

\subsection{Acquisition of the Test Data}
  \begin{figure}[h]
\centering
\subfigure[]{\includegraphics[width=0.6\linewidth]{path.png}}
\subfigure[]{\includegraphics[width=0.13\linewidth]{uav_image.png}}
\subfigure[]{\includegraphics[width=0.13\linewidth]{u-net_image.png}}
\caption{(a) The flight path. UAV took off at the start point and landed at the end point after flying along the square path. (b) Some examples of the image frames taken by the RGB camera during the UAV flight. (c) The images corresponding to (b) in which the buildings are segmented through the segmentation network.}
\end{figure}

\begin{table}[b]
\caption{Parameters used in the experiment}
\label{parameters}
\begin{center}
\begin{tabular}{|c|c|c|}
\hline
Parameter & Value & Meaning\\
\hline\hline
\textit{n} & 3 & 3 features are used in the experiment\\
\hline
$\alpha$ & $43{^\circ}$ & Field of view of the monocular camera\\
\hline
$z_{l}$ & 150m & the flight height is fixed to \textit{$z_{l}$}\\
\hline
$e_{1}$ & 0.3 & feature error threshold in (5)\\
\hline
$\epsilon$ & 25m & distance error constant in (6)\\
\hline
$d_{max}$ & 75m & convergence determination constant\\
\hline
\end{tabular}
\end{center}
\end{table}

  The experimental area is near Daegwallyeong, Pyeongchang, Gangwon-do, South Korea. Fig. 8 (a) shows the satellite images of the area where experiment was conducted. The blue line in Fig.8 (a) is the flight path. The flight was performed in a rectangular path and, to simulate the long distance flight, the experimental data was used under the condition in which loop closing algorithm in the VIO is not used. Fig.8 (b) shows some examples of the UAV images and Fig.8 (c) shows the corresponding building segmented images. The UAV flied at an average speed of 5m/s at a height of 150m. The total flight distance is 1,075m. The UAV images were acquired with 30Hz. From this image sequence, the building ratio matching algorithm was run every 5 sec. This is because the aircraft flied at an altitude of 150m and the aircraft had to move about 25m to make a meaningful change in the input images. Other parameters used are shown in Table {\afterassignment\BigRoman\count255= 1}.

\begin{figure*}[h]
\centering
\subfigure[]{\includegraphics[width=0.31\linewidth]{vins_result_new.png}}
\subfigure[]{\includegraphics[width=0.31\linewidth]{proposed_result1_new.png}}
\subfigure[]{\includegraphics[width=0.31\linewidth]{proposed_result2_new.png}}
\caption{Results of the localization. Green line is the ground truth (trajectory of RTK GPS) and blue line is the estimated trajectory. (a) The trajectory estimated by VINS-mono algorithm. (b) The trajectory estimated by building ratio matching and VINS-mono. The UAV is localized at the red circle (First converged point). To compare with (c), only VINS-mono is used after the first converged point. (c) The trajectory estimated by building ratio matching and VINS-mono. The UAV is localized twice (two red circles).}
\end{figure*}

   \begin{figure}[h]
      \centering
      \includegraphics[scale=0.5]{searching.png}
      \caption{Result of the first candidates calculation (blue). The image on the left shows the candidates estimated at the entire map. The image on the right is an enlarged image of the actual flight area on the entire map.}
      \label{figurelabel}
   \end{figure}

\subsection{Experimental Results}

\begin{table}[b]
\caption{Root-Mean-Square-Error of each result}
\begin{center}
\begin{tabular}{|c c c c|}
\hline
RMSE & VINS-mono & Proposed 1 & Proposed 2 \\
\hline\hline
Whole path & 48.27m & 31.25m & \textbf{7.53m} \\
\hline
After first convergence & 79.58m & 77.12m & \textbf{12.01m} \\
\hline
\end{tabular}
\end{center}
\end{table}

Results of the VIO algorithm and proposed algorithm are shown in Fig. 9. The ground truth trajectory which is obtained from RTK GPS is represented in green color. The trajectories which are estimated by VIO and the proposed algorithm is represented in blue color. VINS-mono\cite{vins-mono} is used as a comparison algorithm. VINS-mono algorithm can not estimate the initial point, thus the start point was given to VINS-mono. Fig. 9 (a) is the trajectory estimated by only VINS-mono. A trajectory which uses only VINS-mono after the first convergence is generated as shown in Fig. 9 (b) for the fair comparison when both VINS-mono and the proposed algorithm when the initial point is known (TABLE {\afterassignment\BigRoman\count255= 2} Proposed 1). Fig. 9 (c) is the trajectory of the proposed algorithm (TABLE {\afterassignment\BigRoman\count255= 2} Proposed 2). In this case, the UAV calculated the candidates until it landed.

The RMSE values are shown in TABLE {\afterassignment\BigRoman\count255= 2}. The whole length of the path is 1,075m. The RMSE of the whole path of Proposed2 which uses the proposed algorithm until it landed is significantly lower than VINS-mono result which are 48.27m and 7.53m. In the case of the initial points of both VINS-mono and the proposed algorithm are known, the RMSE of VINS-mono was 77.12m while RMSE of our proposed method was 12.01m.

      \begin{figure}[t]
      \centering
      \includegraphics[scale=0.55]{candidates_table.png}
      \caption{Some examples of the candidates calculation. This figure starts with the candidate estimated by the first image frame and shows how the candidates are converged.}
      \label{figurelabel}
   \end{figure}

   Fig. 10 shows the first candidates (blue) of the proposed algorithm. The image on the left is the entire numerical map with the candidates. The image on the right is an enlarged image of the actual flight area as shown in Fig. 8 (a). Some examples of the candidates at the actual flight area are described on Fig. 11. The candidates are initially scattered and converge at frame 27. The global position of the UAV is firstly estimated at this frame.

\section{CONCLUSIONS}
In this paper, we proposed a novel localization algorithm, named as BRM (Building Ratio Map based) localization, for a UAV based on building ratio matching using numerical map and the UAV images in GNSS denied environments. The conventional approaches need the initial position to estimate the absolute position of the UAV. However, our algorithm estimates the position of the UAV when the initial position is given within wide area e.g. \textit{$2.2km\times2.6km$}. Since the UAV is used to accomplish a given goal in a particular region, the assumption that the initial location is known in the range of \textit{$km^{2}$} is a resonable assumption and it can be extended to wider areas in the proposed algorithm. The proposed algorithm is verified using real UAV flight data which include UAV images, IMU data, and GPS data.

The experiment has been carried out in Pyeongchang, Gangwon-do. The error is analyzed by comparing with the actual flight path measured by RTK GPS. The candidates have converged twice during the flight with 6.1m and 7.8m errors, respectively, which are small enough since these are after 382m and 891m long flights at a height of 150m. It is confirmed that the RMSE (Root Mean Square Error) of the proposed algorithm is smaller than that of VINS-mono. In the absence of the initial position, the position is better estimated than VINS-mono, and after two convergences, the error about the entire path is much smaller than that of VINS-mono. For the fair comparison when the initial position is known, we analyze the position error after the convergence of candidates. This means that the errors of the proposed algorithm and VINS-mono are compared when both are provided with the starting position. In this case, the result shows that the RMSE is significantly reduced from that of VINS-mono.

In this paper, only building ratio information is used as features. Because of this, it takes time to estimate the initial position. For the same reason, the proposed algorithm is suitable for UAVs flying long distances outdoors. So, it can be applied to long-distance goods delivery systems that use UAVs, and reconnaissance in the air in areas where attacks such as GPS spoofing can occur. However, due to the characteristics of the features used, the current algorithm has limitations in indoor navigation. If the proposed algorithm is improved by learning various classes that are likely to be indoors, it can be extended to indoor navigation in the future. For the future work, various classes such as roads which are mostly present in various environments will be used as feature information as well. Learning of various classes may improve the accuracy and the convergence speed of the proposed algorithm. Additionally, we expect that the improved algorithm will operate robustly in various environments. In such settings, the algorithm can be verified to be more robust with additional experiments in complex urban environments.

\bibliographystyle{IEEEtran}


\section{Data sets} \label{sec:datasets}

\subsection{Artificial data generation} \label{subsec:data-generation}

Both artificial data sets were generated using the GluonTS toolkit. More specifically,
\begin{itemize}
	\item
	$D_2$ was generated using `ComplexSeasonalTimeSeries` with the following parameters: blabla
	\item
	$D_3$ was generated through the general `RecipeDataset` class with the following recipe details: blabla
\end{itemize}

\subsection{More plots}

\begin{figure}[!htb]
  \begin{subfigure}[t]{1\textwidth}
    \centering
	\includegraphics[width=0.32\linewidth]{./figures/data/datasets/artificial/easy/test/_1.pdf}
	\includegraphics[width=0.32\linewidth]{./figures/data/datasets/artificial/easy/test/_2.pdf}
	\includegraphics[width=0.32\linewidth]{./figures/data/datasets/artificial/easy/test/_3.pdf}
    \caption{Easy artificial data samples}
  \end{subfigure}

  \medskip

  \begin{subfigure}[t]{1\textwidth}
    \centering
	\includegraphics[width=0.32\linewidth]{./figures/data/datasets/artificial/hard/test/_4.pdf}
	\includegraphics[width=0.32\linewidth]{./figures/data/datasets/artificial/hard/test/_6.pdf}
	\includegraphics[width=0.32\linewidth]{./figures/data/datasets/artificial/hard/test/_8.pdf}
    \caption{Hard artificial data samples}
  \end{subfigure}

  \medskip

  \begin{subfigure}[t]{1\textwidth}
    \centering
	\includegraphics[width=0.32\linewidth]{./figures/data/datasets/icc/test/_1.pdf}
	\includegraphics[width=0.32\linewidth]{./figures/data/datasets/icc/test/_2.pdf}
	\includegraphics[width=0.32\linewidth]{./figures/data/datasets/icc/test/_3.pdf}
    \caption{Real traffic data samples}
  \end{subfigure}
  \caption{title}
\end{figure}

\newpage

\section{GluonTS forecasting models details}

We provide here some details about the predictive models that were used to forecast the traffic.

\paragraph{Hyper-parameters}

blabla

\newpage

\section{Results}

blabla

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

Modern cloud computing providers automate the scaling of hardware resources to meet
the demand of the hosted application. The key consideration when designing scaling policies is
the trade-off between reducing costs by minimizing the allocated resources
and satisfying customers by supplying sufficient resources for the application to run nominally.
Auto-scaling is a well-studied topic for which comprehensive
surveys and reviews are available \cite{Qu2018Survey, makridakis2018m4}.

We consider complex, large-scale applications with hundreds or thousands of
servers, for which reactive auto-scaling policies often exhibit shortcomings
stemming from unrealistic assumptions.
For example, while releasing hosts is usually fast, there
are limits to how quickly the cloud provider can fulfill requests to add more
instances to a fleet, and in particular for large volumes.
Reactive scaling assuming instantaneous resource addition will then lag behind the true demand and
the fleet may be at risk of being under-provisioned before a traffic peak and over-provisioned
afterwards.
Predictive scaling can circumvent these limitations by using a
forecast of the application's workload to make scaling decisions ahead of time while
taking throughput constraints into account.
Using a probabilistic forecast, we can design policies that are optimal for the
level of risk aversion of the application owner, in the sense that any quantile
of the random workload can be estimated.
This approach is particularly well suited for applications with
strongly seasonal (daily) traffic patterns and for large fleet sizes, for
which the time to scale the fleet up is significant.

Our contributions to predictive auto-scaling for cloud applications are as
follows.
\textit{\textbf{(i)}} We develop an approach for analyzing and evaluating auto
scaling policies for cloud applications. Our approach includes a risk aware cost
function and a realistic model for the scaling behavior of the application
taking into account real world constraints on throughput and latency.
\textit{\textbf{(ii)}} We show how \emph{probabilistic} forecasts naturally lend
themselves to satisfy the risk aversion of the owner of the application.
\textit{\textbf{(iii)}} We show through experiments that a simple heuristic becomes
optimal in the case of high risk aversion which is the practically most relevant
case.

Similar heuristics are used to auto-scale over 40\,000 of Amazon's internal
auto-scaling groups, as well as external applications making use of AWS
Auto-Scaling\footnote{https://aws.amazon.com/autoscaling/}.

We discuss related work in Sec.~\ref{sec:related-work}. We then introduce and
formalize the predictive auto-scaling problem and state our modeling assumptions
in Sec.~\ref{sec:problem}. We present the scaling policies and evalute them
empirically in Sec.~\ref{sec:scaling_policies} using synthetic and real-world
data.

\section{Related Work} \label{sec:related-work}
Scaling resources in cloud environments is a mature and active area of research,
see e.g.~\cite{Qu2018Survey, makridakis2018m4}
for reviews.
Predictive scaling has been studied alongside reactive scaling
(\cite{Wu2019Anna} is a recent example in the context of storage systems) and
products such as AWS Auto-Scaling now offer both predictive and reactive scaling
techniques. Hybrid methods taking advantage of the strengths of both schemes
have also been studied~\cite{coordinated}.

Approaches that rely on \emph{forecasting} techniques
(like~\cite{Messias2016,Kirchhoff19,Lang2016,Ma2018,taft2018}) are the most
closely related to our work
but state of the art probabilistic forecasts are rarely considered.
They are, however, key for optimal decision
making~\cite{Faloutsos2018,Faloutsos2019}.
This is well-known in the supply chain literature~\cite{simchi13logistics}.
In many settings, neural network architectures for forecasting have been shown
to deliver superior predictive accuracy than traditional
techniques~\cite{smyl2018m4,deepAR,rangapuram2018deep}. They are used in AWS
auto-scaling for this reason.

An exception is~\cite{Borkowski2019} which relies on a probabilistic time series
model to infer the true state of the application (but not to forecast workload).
The discovery of this state allows them to handle the cost of scaling
operations, a problem which is more relevant in a streaming scenario than in our
scenario.
Their use of multiple metrics in a multivariate time series models is relevant
for the subject matter of this paper,
it should be explored in future work
(e.g., through using modern multivariate forecasting models~\cite{multivariate}).

\section{Problem statement} \label{sec:problem}

In this section we formalize the predictive scaling problem and discuss our modeling assumptions.
We consider a fleet of servers that run a horizontally scalable application.
Our goal is to scale the number of servers in an \emph{optimal} way to match the actual demand.

We assume a work-load estimation model to be known, outputting an approximation
of the demand (ideal number of hosts) for the application,
given an observable metric.\footnote{
  This is no small assumption. Depending on the application,
  modeling the work-load may require significant work and it is a research area
  in its own right e.g.,~\cite{Lang2015,SchaffnerJ13,Duggan2011,Mozafari2013,Ma2018,das2016automated}.
}
The selected target metric must not be derived from the scaling decisions for
the forecasting model to be an unbiased predictor of the demand.
The sum of the number of bytes received and sent by every server in the fleet is one such measure, which we found to be pertinent as it is strongly correlated with the
CPU and memory usage of the servers.

We denote this variable by $v_t$ and a simple linear workload model yields that
the number of desired hosts $z_t$ is proportional to the measured load $z_t = \xi v_t$.
In practice the function linking $v_t$ to $z_t$ may be more complex and
estimated with more sophisticated models.

\subsection{Probabilistic Forecast and cost function}\label{subsec:forecast_and_cost}
Given the workload representation $v_t$, a forecasting model can be trained to
generate estimates of the probability distribution for the future
  $P(v_{t+1}, v_{t+2}, \dots, v_{t+T} \mid \mbox{past})$.\footnote{
  We assume discrete time steps. Real cloud events and scaling actions occur in continuous time,  however a discretization in steps that
  match real scaling time-frames avoids requiring complex time series models
  such as~\cite{turkmen19}.
  Delays for the host provider are a few minutes for releasing or a few dozens
  of minutes to acquire new hosts.
  Therefore, a discretization into minute steps seems a non-restrictive assumption.
}
In practice, workload trends are dominated by daily and weekly patterns.

In order to make optimal scaling decisions we need to define the cost function
that describes the cost for over-provisioning and under-provisioning the fleet.
We assume that the cost is composed of two terms.
The first one is hardware cost, which
scales linearly with the number of provisioned hosts.
The second term captures the service's performance. This may be, for instance, loss of revenue or customers
caused by the increase in latency when the service is under-provisioned.

The latter part is harder to quantify and measure. We make the following
assumption: at each time step $t$ there is a critical capacity $r^\star_t$ such
that the service will perform optimally if the provisioned capacity $r_t$ is
above $r^\star_t$.
Adding more capacity than $r^\star_t$ does not improve the performance of the
service while increasing the operational cost. If $r_t < r^\star_t$, the
application's performance deteriorates. The cost associated with this is
proportional to $r^\star_t - r_t$, typically with a large proportionality constant.
In practice $r^\star_t$ is never observed, but $z_t$ should approximate it
closely if the workload model is good.

This translates into an asymmetric cost function which is, up to a constant
factor, given by the $\alpha$-quantile loss function (sometimes also called the pinball loss,
see~\cite{gneiting2007probabilistic} for a thorough discussion) described in
Equation~\eqref{eq:pinball}.
\begin{equation}\label{eq:pinball}
  \Lambda_\alpha(r_t, z_t) = \begin{cases}
    (1 - \alpha) \cdot (r_t - z_t) & \mbox{if }r_t > z_t \;,\\
    \alpha \cdot (z_t - r_t) & \mbox{otherwise.}
  \end{cases}
\end{equation}
The quantile $\alpha\in (0,1)$ corresponds to the risk aversion factor.
Note that $z_t$ is a random variable as the forecast of $v_t$ is probabilistic.
The expected cost $\mathbb{E}[ \Lambda(r_t, z_t) ]$ is minimized when $r_t =
q_\alpha(z_t)$, i.e., the $\alpha$-quantile of $z_t$.

\subsection{Model for cloud provider behavior}

A central entity in our problem statement is the cloud or host provider that
provisions the instances for the fleet (e.g.\ the AWS EC2 service).
The application owner can request new instances from this provider or release
superfluous instances.

While the exact behavior of the host provider may be complicated, we make the
following simplifying assumptions:
\textit{(i)} the provider has an infinite reservoir of instances,
\textit{(ii)} releasing an instance back to the provider is instantaneous and
\textit{(iii)} requests for new instances are handled as depicted in
Figure~\ref{fig:providermodel}.
That is, the provider has a fixed number of slots for provisioning,
which can be empty (white) or occupied (gray).
The number of hosts requested by the user are summed to $R \ge 0$.
When slots are empty they are immediately filled and $R$ is decreased.
When a slot is filled, the corresponding instance is being provisioned and it
takes a (random) time $\tau$ until the host is available.
At this point the slot becomes empty again.
For most applications, these assumptions were found not to be overly
restrictive.
\begin{figure}[htbp]
	\centering
	  \begin{tikzpicture}[scale=0.4]
	  \node[draw, rounded corners=4pt,align=center] at (4,3) {$R$ outstanding requests};
	  \draw [rounded corners=4pt] (0,0) rectangle ++ (8,1);
	  \draw[fill=lightgray] (0.5,0.5) circle (0.3cm);
	  \draw[fill=lightgray] (1.5,0.5) circle (0.3cm);
	  \draw[fill=white] (2.5,0.5) circle (0.3cm);
	  \draw[fill=lightgray] (3.5,0.5) circle (0.3cm);
	  \draw[fill=white] (4.5,0.5) circle (0.3cm);
	  \draw[fill=white] (5.5,0.5) circle (0.3cm);
	  \draw[fill=white] (6.5,0.5) circle (0.3cm);
	  \draw[fill=lightgray] (7.5,0.5) circle (0.3cm);
	  \draw [->,line width=0.8] (4,2.3) --  node [text width=4cm,midway,right] {wait for free slot} (4,1.2);
	  \draw [->,line width=0.8] (8,0.5) --  node [] {} (10.8,0.5);
	  \node[draw, rounded corners=4pt,align=center] at (15,0.5) {alive after time $\tau \sim P(\tau)$};
	  \end{tikzpicture}
	 \caption{\label{fig:providermodel}Process for provisioning new hosts.}
\end{figure}

We denote by $\tau$ the delay between asking for additional resources and those
resources becoming available.
Its distribution $P(\tau)$ depends not only on how quickly the
provider can spin up a new instance, but also on the time it takes
to install and deploy the necessary packages, the data (e.g.\ cache), and to start
the application.
Typical values for real-world applications are between
a few minutes and half an hour.
This can vary significantly with the application and the deployment approach.
For these reasons, the resource type can exhibit a complex structure as
illustrated by the histogram of delays in Figure~\ref{fig:providerdelay}.
\begin{figure}[htpb]
  \centering
  \includegraphics[width=0.7\linewidth, height=3cm]{./amper_delay.pdf}
  \\[-3.5ex]\caption{
    \label{fig:providerdelay}Histogram of $P(\tau)$ from individual
    scaling events of a real application (count on vertical axis).
  }
\end{figure}

When the fleet size (and the corresponding scaling amounts) is large, the
limited number of slots in Figure~\ref{fig:providermodel} translates into an
effective rate limit $\rho$ for the number of hosts that can be added per unit
of time.
Our assumptions are simplications: the delays and the throughput $\rho$ are likely to be coupled, to vary over time, and to depend on the current and past request values. However, we found
these assumptions to approximate the real world behavior well.

\subsection{Scaling Cost Model}
Next, we combine the aforementioned models to calculate the \emph{expected cost} of the scaling policies.

Let us fix a time horizon $n > 0$ and consider a probabilistic forecast of length
$n$ generated at time $0$.
Using Equation~\eqref{eq:pinball}, the expected cost for a sequence of future resource values $r = (r_1, \dots, r_n)$ is
\begin{align}
  \mathcal{L}_\alpha(r) &= \sum_{i=1}^{n} \int \Lambda_\alpha(r_i, z_i) P(z_1, \dots, z_n)  \,d z_1 \dots d z_{n} \nonumber\\
                       &= \sum_{i=1}^{n} \mathbb{E}_{z_i\sim P(z_i)} \; \Lambda_\alpha(r_i, z_i) \;.
\end{align}
We denote the actual cost by $\mathcal{L}_\alpha^\star(r) = \sum_{i=1}^{n}
\Lambda_\alpha(r_i, r^\star_i)$, which is dependent on the true demand $r^\star$
and cannot be directly minimized.
If the workload model is adequate (that is, $z$ is a good estimate of
$r^\star$), minimizing $\mathcal{L}_\alpha(r)$ should provide a reasonable
resources pattern for the loss $\mathcal{L}_\alpha^\star(r)$. The choice of
$\alpha$ depends on the relative costs of the service being under- or
over-provisioned as discussed in Section~\ref{subsec:forecast_and_cost}.

The optimal solution of $\mathcal{L}_\alpha(r)$
is $r = \left(q_\alpha(z_1), \dots, q_\alpha(z_n)\right)$ where $q_\alpha(z_i)$
is the $\alpha$-quantile of $z_i$.
However, we cannot act on $r_t$ directly; instead, at each time step $t \in \mZ$
we are able to request $\posreq_t \in \mathbb{N}$ and release $\negreq_t \in
\mathbb{N}$ hosts.
As pointed out above, releases are instantaneous while every single requested
instance will arrive $\tau$ time steps later.
As we have a forecast of the demand for the next $n$ time steps,
we want to optimally pick the values $\posreq_0, \dots, \posreq_n$ and
$\negreq_0, \dots, \negreq_n$,
given the past requests $\posreq_{t' < 0}$ and $\negreq_{t' < 0}$ and the past
resources values $r_{t' \le 0}$.
The larger the forecast horizon $n$, the earlier scaling decisions can be taken.
But forecasting accuracy decreases as the horizon increases, so the choice of
$n$ is a trade-off that depends on the periodicity of the time-series and on the
effectiveness of the forecasting model.

\header{Resources estimation} Denote by $r_{-\infty}$ the host count at the start of
the auto-scaling process, before taking any decision.
If we are able to estimate accurately the host provider's behavior, we can
compute for each $t$ an approximation $\hat r_t$ of $r_t$ as a function of $\posreq_{t' \le t}$ and $\negreq_{t' \le t}$.
Let $\hat\delay$ be an estimation of the true host provider random delay
for acquiring new instances.
Then the random resource estimation at time $t$ can be written as
\begin{equation}
\hat{r}_t(q, f) = r_{-\infty} + \sum_{i \le t} \bigg( \sum_{j = 1}^{\posreq_i}
  \mathbbm{1}\{i + \hat\delay_{i, j} \le t\} - \negreq_i \bigg)
\end{equation}
where $(\hat\delay_{i, j})_{i, j \in \mZ} \overset{i.i.d.}{\sim} \hat\delay$
is a collection of random variables representing the delays of the potential
positive requests.
When the number of hosts is large, it may be approximated in
expectation in the following way
\begin{equation}
  \hat{r}_t(q, f) = r_{-\infty} + \sum_{i \le t} \Big( \posreq_i \cdot
  \mathbbm{P}\{i + \hat\delay \le t\} - \negreq_i \Big)
\end{equation}
by using the fact that, in expectation, a request $q$ done at time $t$ will yield $q \cdot
\mathbbm{P}\{\hdelay = \delta\}$ hosts $\delta$ time steps later.

\header{Optimization Problem Formulation} In summary, we wish to solve the following optimization problem,
given that $\posreq_{t' < 0}$ and $\negreq_{t' < 0}$ are fixed.
\begin{equation}\label{eq:optimization_problem}
\begin{aligned}
	& \underset{\substack{\posreq_0, \dots, \posreq_n, \\ \negreq_0, \dots, \negreq_n}} {\text{minimize}}
	& & \mathbb{E}_{z} \, \mathbb{E}_{\hat{r}(\posreq, \negreq)} \, \mathcal{L}\left( \hat{r}(\posreq, \negreq), z \right)  \\
	& \text{subject to}
	& & \posreq_i \le \hat\rho, \negreq_i \leq r_i, \posreq_i \in \mathbb{N}, \negreq_i \in \mathbb{N}
\end{aligned}
\end{equation}
where $\hat\rho$ is an estimation of the throughput $\rho$.
The objective function in Equation~\eqref{eq:optimization_problem} is convex
as the quantile loss is convex and $\hat{r}$ is linear.
The integrality constraints make the optimization challenging but, assuming a
large-enough number of hosts, a relaxation followed by a randomized rounding produce
near-optimal solutions.

The double expectation in Equation~\eqref{eq:optimization_problem} may be
approximated with Monte-Carlo averaging; it only requires to sample from $z$ and
from $\hat r(q, f)$ (the latter being done with samples of the delay $\hat\tau$).
From here, the optimization of Equation\eqref{eq:optimization_problem} results in minimizing
a convex piecewise-linear function. This problem is equivalent to solving a
linear-program, which can be done efficiently.
Once the objective is optimized, we can request the $k$ first values at the
right time and repeat the process $k$ time steps later in what amounts to
lookahead optimization.

The variant of the auto-scaling problem that we consider in this paper is, to
the best of our knowledge, new to the literature.
Prior work either assumes the cloud to be fully elastic and resources to come
available immediately or with a delay (e.g., ~\cite{Wu2019Anna}),
but the constraints introduced by the limited throughput of the cloud provider have not yet been studied.

\section{Policies \& Experiments}\label{sec:scaling_policies}

In the following, we describe the scaling policies that we compare empirically.
Forecasts are made every hour, over two days, and with 5 minutes granularity.
Every forecast overrides the previous one as it is considered to be more accurate.
Whenever a new forecast is made, new decisions are taken for the next hour,
where requests and releases are planned with a granularity of 5 minutes.
The choice of these parameters depends greatly on the workload patterns that
must be addressed;
we found these to be adapted to the daily trends that are very common at AWS.

\header{Policies}
The following two policies are frequently used in industrial predictive
auto-scaling applications~\cite{Qu2018Survey}.
The \emph{maximum observed needed capacity} policy keeps the host count at the
maximum that was needed in a recent past (for example one day or one week).
That kind of policy is very conservative but commonly employed due to its robustness.
Another baseline, \emph{reacts periodically to the current need}.
That is, every 5 minutes, it estimates the desired host count $z = \xi v$ and requests
resources accordingly.
The estimation can be adjusted by some factor in order to be more or less conservative.
This is better known as \emph{reactive} scaling and usually performs well if the
resource needs do not vary too rapidly.

We introduce the simple \emph{forecast shifting} policy based on forecasts that
is described in Algorithm~\ref{alg:forecast_shifting}. It works in three steps:
\emph{(i)} computing the $\alpha$-quantile of the forecast $z$ and adjusting it
backward to take into account the limiting throughput $\rho$,
\emph{(ii)} splitting the positive and negative requests, and
\emph{(iii)} shifting the positive ones by the mean of the delay $\hat\tau$ to
call for them earlier.
\begin{algorithm}
    \caption{Forecast shifting}\label{alg:forecast_shifting}
    \begin{algorithmic}[1]
        \State \textbf{Input:} Probabilistic forecast $z = (z_1, \dots, z_n)$,
        guessed random delay $\hat{\tau}$ and throughput $\hat{\rho}$, current host count $r_0$
        \State $z'_t \gets \text{quantile}_\alpha(z_t) \; \forall t$ \algorithmiccomment{\emph{(i)}}
        \For{$t = n - 1,\dots, 1$}
            \State $z'_t := \max(z'_t, z'_{t + 1} - \hat{\rho})$
        \EndFor
        \State $\Delta_t \gets z'_t - z'_{t - 1}$ for $t = 2,\dots, n$ and $\Delta_1 = z'_1 - r_0$
        \algorithmiccomment{\emph{(ii)}}
        \State $q,\,f \gets \max(\Delta, 0),\,\min(\Delta, 0)$
        \State Shift $q$ in the past according to $\mathbb{E}[\hat\delay]$ \algorithmiccomment{\emph{(iii)}}
        \State Request $q$ and release $f$
    \end{algorithmic}
\end{algorithm}
Figure~\ref{fig:qshift} (bottom) illustrates the backward adjustment performed in step
\emph{(i)}. The forecasts (orange) must be compensated to take into account the
throughput $\rho$.
This algorithm can be adapted to take into account the full distributions of the
delay rather than only the expectation. However, all our experiments showed that
considering the full distribution was not leading to any loss improvement.

For large risk aversion the cost of underpredicting is much higher than over
prediction. Without constraints the $\alpha$-quantile is the optimal solution
for this asymmetric trade-off. In the limit $\alpha \to 1$ the optimal capacity
with constraints is equal to the $\alpha$-quantile whenever possible but never
below it. This is what the forecast shifting method achieves.

The \emph{optimal policy} is to solve the problem as posed in
Equation~\eqref{eq:optimization_problem} with standard mathematical programming
frameworks like SCIP~\cite{GleixnerBastubbeEifleretal.2018}.
\begin{figure}[!h]
  \begin{subfigure}[t]{0.5\textwidth}
    \begin{tikzpicture}
      \centering
      \node (img1) {
        \includegraphics[width=\linewidth, height=3cm]{./optim_risk09-edited.pdf}
      };
      \node[left=of img1, node distance=0cm, rotate=90, anchor=center,yshift=-1.1cm,font=\color{black}] {count};
    \end{tikzpicture}
  \end{subfigure}\\[-2.5ex]
  \begin{subfigure}[t]{0.5\textwidth}
    \begin{tikzpicture}
      \centering
      \node (img1) {
        \includegraphics[width=\linewidth, height=3cm]{./quantile_risk09.pdf}
      };
      \node[below=of img1, node distance=0cm, yshift=1.15cm,font=\color{black}] {time [5 mins]};
      \node[left=of img1, node distance=0cm, rotate=90, anchor=center,yshift=-1.1cm,font=\color{black}] {count};
    \end{tikzpicture}
  \end{subfigure}\\[-4ex]
  \caption{
    \label{fig:qshift}
      Host counts for optimization (top) and shifting (bottom) policies
      simulated on $\boldsymbol{D}_1$. The quantile $\alpha = 0.9$ of the
      forecasts cannot be always matched because of the limiting throughput
      $\rho$.
  }
\end{figure}

\header{Datasets and predictors} In order to compare the 4 scaling methods
mentioned above, we run experiments in a simulated environment (using the
framework SimPy~\cite{simpy}) on one real-world ($\boldsymbol{R_1}$) and two
artificial ($\boldsymbol{D_1}$ and $\boldsymbol{D_2}$) datasets.
The real-world dataset $\boldsymbol{R_1}$ is composed of $1\,000$ randomly
selected Amazon auto-scaling groups over a period of 6 hours.
The datasets $\boldsymbol{D_1}$ (low noise) and $\boldsymbol{D_2}$ (high noise)
are generated with strong daily/weekly patterns, some linear trends, and
Gaussian noise. For all of them, a DeepAR~\cite{deepAR} forecasting model is
trained and makes new predictions every hour.

\begin{figure}[!h]
  \begin{subfigure}[t]{0.5\textwidth}
    \begin{tikzpicture}
	    \centering
      \node (img1) {
        \includegraphics[width=\linewidth, height=4cm]{./box_plot_D1_losses_comparison.pdf}
      };
      \node[left=of img1, node distance=0cm, rotate=90, anchor=center,yshift=-1.2cm,font=\color{black}] {loss};
    \end{tikzpicture}
  \end{subfigure}\\[-3ex]
  \begin{subfigure}[t]{0.5\textwidth}
    \begin{tikzpicture}
	    \centering
      \node (img1) {
        \includegraphics[width=\linewidth, height=4cm]{./box_plot_D2_losses_comparison.pdf}
      };
      \node[left=of img1, node distance=0cm, rotate=90, anchor=center,yshift=-1.2cm,font=\color{black}] {loss};
    \end{tikzpicture}
  \end{subfigure}\\[-3ex]
  \begin{subfigure}[t]{0.5\textwidth}
    \begin{tikzpicture}
	    \centering
      \node (img1) {
        \includegraphics[width=\linewidth, height=4cm]{./box_plot_ICC_losses_comparison.pdf}
      };
      \node[below=of img1, node distance=0cm, yshift=1.35cm,font=\color{black}] {policy};
      \node[left=of img1, node distance=0cm, rotate=90, anchor=center,yshift=-1.2cm,font=\color{black}] {loss};
    \end{tikzpicture}
  \end{subfigure}\\[-4ex]
	\caption{
    \label{fig:los_comparison}
    Losses obtained on the datasets $\boldsymbol{D}_1$ (top), $\boldsymbol{D}_2$
    (middle), and $\boldsymbol{R}_1$ (bottom) for the presented policies.
  }
\end{figure}

\header{Results} Figure~\ref{fig:los_comparison} provides the costs
associated to the scaling policies.
The first two columns, \emph{Max week} and \emph{Max day} refer to the first heuristic,
where the maximum is taken over an entire week and a day respectively, leading to over-capacity.
\emph{Instant} and \emph{Instant (tuned)} refer to the reactive policy.
\emph{Shift} refers to the forecast shifting procedure depicted in
Algorithm~\ref{alg:forecast_shifting} and \emph{Optim 500} is the optimization
scaling policy for which the solver completes 500 iterations.

From all panels of Figures~\ref{fig:los_comparison} it is apparent that the
forecasting shifting policy compares favorably to other approaches in terms of
costs, in particular with respect to the optimal policies. This is particularly
visible on the two simulated datasets. On real-world data the median loss using
the shifting policy is larger than that of the Max day and Max week policies,
however the latter two result in larger extreme losses.

Figure \ref{fig:qshift} shows the outcome of the optimization and the shift
policies on a time-series from $\boldsymbol{D}_1$.
It turns out that when $\alpha$ is close to 1, the scaling patterns of the
forecast shifting method converge to the ones obtained by optimizing the
objective function~\eqref{eq:optimization_problem}.
We also observe that the forecasting shifting policy comes at a fraction of the
optimization computational cost. Therefore, a practical predictive auto-scaling
solution should adopt this policy.
The datasets increase in strength of seasonality in relation to the noise level
as $\boldsymbol{R}_1 \to \boldsymbol{D}_2 \to \boldsymbol{D}_1$.
Intuitively, the stronger this signal to noise ratio is, the more important
seasonal scaling becomes.
For weak seasonality the simple max scaling baselines already work well -- no
scaling is necessary in the limit of very weak seasonality.
For large noise, the forecasting and optimization problems also become more
challenging, which increases the variance for the optimization method in this
case.

\section{Conclusion} \label{sec:conclusion}

We introduced a predictive auto-scaling problem formulation which takes a random
delay of resource availability and a limiting throughput into account.
Incorporating probabilistic forecasts, as opposed to merely point forecasts, we
obtain a rigorous mathematical optimization problem formulation for which we
provide a heuristics that works as well as optimal solutions in practically
relevant scenarios at a fraction of the compute cost.

\part*{Appendix}
\end{document}


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

As a byproduct of domestic (sewage) and industrial wastewater treatment plants, waste activated sludge (WAS) is pervasive.  WAS is comprised of different groups of microorganisms, organic and inorganic matter agglomerated together in a polymeric network formed by microbial extracellular polymeric substances (EPS) and cations~\cite{Li:1990}.  Treating and then disposing of WAS is costly and can negatively impact the environment. Anaerobic digestion, in which microorganisms are used to breakdown organic matter in the absence of oxygen, is commonly used to stabilize and reduce the volume of the WAS.  In the final stage of anaerobic digestion, a methane-rich biogas is produced which can be collected and later used as fuel.  The time it takes for anaerobic digestion to reach completion can be as high as 20~days and can be a bottleneck in the treatment of wastewater.  For a comprehensive review of anaerobic digestion of WAS see Ref~\cite{Appels:2008}.

Efforts to pretreat WAS before anaerobic digestion to enhance the rate of methane production have shown promise.  The pretreatment methods include thermal cycling, sonication, application of high pressure, chemical and enzymatic treatment, microwave irradiation, and exposure to high-intensity pulsed electric fields~\cite{Appels:2008, Eskicioglu:2009, Saha:2011, Salerno:2009, Lee:2010, Sheng:2011}.  The goal of all pretreatment methods is to enhance the bioavailability of the WAS to the methane-forming microorganisms in the anaerobic digester.  These enhancements can be achieved by puncturing or otherwise damaging cell membranes and/or cell walls, breaking down or reducing the characteristic size of the floc structure, and converting complex organic molecules to simpler forms~\cite{Salerno:2009}.

To this point, the microwave treatments have been restricted to frequencies of 2.54~GHz and 915~MHz due to the availability of inexpensive high-power microwave tubes (magnetrons) used in conventional and industrial microwave ovens~\cite{Saha:2011,Park:2004,Eskicioglu:2007,Hong:2006}.  However, there are no fundamental reasons to suppose that these are the optimal frequencies with which to pretreat WAS.  Additionally, the method of irradiating the samples must be carefully evaluated.  A load impedance presented by the sample that is poorly matched to the source impedance of the generator will result in an inefficient system in which not all of the incident power is absorbed by the sample.  The objective of this work was to characterize the complex permittivity of WAS over a broad frequency range and to identify the mechanisms leading to the observed dispersions.  The measurements and theoretical insights into the active microscopic mechanisms allow for the selection of candidate frequencies with which to optimize the electromagnetic pretreatment of WAS.  Furthermore, custom applicators can be designed and constructed such that the incident signal power is very efficiently transferred to the sample.

\section{Experimental Methods}
\label{sec:methods}

When a material is exposed to a time-harmonic electric field $\mathbf{E}e^{j\omega t}$ of angular frequency $\omega$, the total current density $\mathbf{J}$ is the sum of the conduction and displacement current densities:
\begin{equation}
\mathbf{J}=\left(\sigma_\mathrm{dc}+j\omega\varepsilon_0\varepsilon_\mathrm{r}\right)\mathbf{E}=\left[\left(\sigma_\mathrm{dc}+\omega\varepsilon_0\varepsilon^{\prime\!\prime}\right)+j\omega\varepsilon_0\varepsilon^\prime\right]\mathbf{E},\label{eq:J}
\end{equation}
where $\varepsilon_0$ is the permittivity of free space and $\sigma_\mathrm{dc}$ and \mbox{$\varepsilon_\mathrm{r}=\varepsilon^\prime-j\varepsilon^{\prime\!\prime}$} are the dc conductivity and relative permittivity of the material.  In this work, the electrical properties ($\varepsilon^\prime$, $\varepsilon^{\prime\!\prime}$, and $\sigma_\mathrm{dc}$) of WAS were measured using open-ended coaxial probes.  The open-ended coaxial probe was chosen because it is a broadband ($10^5$--$10^{10}$~Hz) measurement suitable for liquid or semisolid samples and requires no sample preparation.  There exist numerous other measurement techniques.  Parallel plate geometries are useful for frequencies below 30~MHz.  Transmission measurements are broadband but, in their simplest configuration, require samples machined to fit tightly inside a waveguide or coaxial transmission line.  Resonant cavities are very precise, but operate only at discrete frequencies~\cite{Kaatze:2006}.

The probes were fabricated from short ($\approx 20$~cm) sections of commercially available semirigid coaxial transmission line.  For frequencies less than 1~GHz, UT-141 coaxial cable was used (shield inside diameter \mbox{$D=2.99$~mm}, center conductor diameter \mbox{$d=0.92$~mm}).  At higher frequencies, evanescent modes at the probe tip can be excited resulting in a breakdown of the measurement technique.  These unwanted modes can be suppressed by switching to a smaller diameter probe.  Measurements up to 40~GHz were made using UT-085 coaxial cable (\mbox{$D=1.67$~mm}, \mbox{$d=0.51$~mm}).  One end of the probe was fitted with a connector and the opposite end was polished flat creating an open circuit condition.  During the measurements, the open end of the probe was submerged into the material under test (MUT) resulting in a mismatched load.  Signals incident from a vector network analyzer (VNA) are partially reflected at the MUT-probe interface.  By analyzing the amplitude and phase of the reflection coefficient as a function of frequency, the complex permittivity and dc conductivity of the MUT were determined.  The experimental geometry is shown in Fig.~\ref{fig:fringe}.
\begin{figure}[t]
\centerline{\includegraphics[width=0.7\columnwidth]{Fig1.pdf}}
\caption{A schematic of the experimental geometry showing the open end of the coaxial probe submerged in the MUT.  The calibration plane of the VNA is established at the opposite end of the probe.  The reflection coefficient $\rho_\mathrm{m}(\omega)$ is measured directly by the VNA and can be related to the reflection coefficient $\Gamma(\omega)$ at the probe tip.  The dotted lines represent electric field lines in the vicinity of the interface between the probe tip and the MUT.  The impedance of the interface is modeled as parallel shunt capacitances $C_\mathrm{f}$ and $\left(\varepsilon_\mathrm{r}-j\sigma_\mathrm{dc}/\omega\varepsilon_0\right)C_0$.}
\label{fig:fringe}
\end{figure}

The reflection coefficient $\Gamma$ at the probe tip is determined from the characteristic impedance $Z_0=50~\Omega$ of the coaxial transmission line and the effective load impedance $Z_\mathrm{L}$ present at the interface between the tip and the MUT:
\begin{equation}
\Gamma(\omega)=\frac{Z_\mathrm{L}-Z_0}{Z_\mathrm{L}+Z_0}.\label{eq:Gamma}
\end{equation}
This load impedance is typically modeled as two parallel shunt capacitors $C_\mathrm{f}$ and $C(\varepsilon_\mathrm{r},\sigma_\mathrm{dc})$. Here, $C_\mathrm{f}$ is purely reactive and accounts for fringing of the electromagnetic fields that occurs within the dielectric of the coaxial probe whereas $C(\varepsilon_\mathrm{r},\sigma_\mathrm{dc})=\left(\varepsilon_\mathrm{r}-j\sigma_\mathrm{dc}/\omega\varepsilon_0\right)C_0$ accounts for fringing that occurs within the MUT~\cite{Stuchly:1980}.  The load admittance at the probe tip can then be written:
\begin{equation}
Z_\mathrm{L}^{-1}=C_0\left(\omega\varepsilon^{\prime\!\prime}+\frac{\sigma_\mathrm{dc}}{\varepsilon_0}\right)+j\omega\left(C_\mathrm{f}+\varepsilon^\prime C_0\right).\label{eq:ZL}
\end{equation}
The measured reflection coefficient $\rho_\mathrm{m}$ at the calibration plane of the VNA can be related to $\Gamma$ by modeling the coaxial probe as a two-port network.  The relevant scattering parameters, which account for losses along the length of the probe and for spurious reflections, were determined using three calibration terminations: short-circuit, open-circuit, and a standard liquid of known permittivity~\cite{Kraszewski:1983}.  A measurement of the real and imaginary components of $\Gamma$ allows Eqns.~\ref{eq:Gamma} and \ref{eq:ZL} to be used to determine both $\varepsilon^\prime$ and $\varepsilon^{\prime\!\prime}+\sigma_\mathrm{dc}/\omega\varepsilon_0$ at each measurement frequency.

Prior to characterizing the electrical properties of WAS, the experimental technique was verified using two control samples.  In both cases, distilled water at \mbox{$25.0\pm0.1^\circ$C} was used as the standard liquid termination.  The first control sample was 99.9\% methyl alcohol which has negligible $\sigma_\mathrm{dc}$.  The relative permittivity extracted from $\Gamma$ was in excellent quantitative agreement with earlier measurements made by Bao and coworkers~\cite{Bao:1996}.  Our measurements of methyl alcohol extended from 10~MHz to 40~GHz.  The lower end of the experimental bandwidth is set by the resolution of the reflection coefficient measurement while the upper end is set by the onset of higher-order evanescent modes at the probe-MUT junction.

The second control sample was a 0.03~molarity solution of 99.8\% NaCl dissolved in distilled water.  The real and imaginary parts of \mbox{$\varepsilon_\mathrm{r}-j\sigma_\mathrm{dc}/\omega\varepsilon_0$} of the solution were extracted from the reflection coefficient measurements.  As expected, at high frequencies we achieved an accurate measurement of $\varepsilon_\mathrm{r}$ of pure water~\cite{Buchner:1999}.  At sufficiently low frequencies, $\varepsilon^{\prime\!\prime}$ is dominated by the $\sigma_\mathrm{dc}$ term such that the imaginary component of \mbox{$\varepsilon_\mathrm{r}-j\sigma_\mathrm{dc}/\omega\varepsilon_0$} varies as the inverse of frequency.  The intrinsic behavior of the real component, on the other hand, is independent of $\sigma_\mathrm{dc}$.  However, when submerged in an electrolyte, the coaxial probe acquires a surface charge resulting in a so-called electric double layer at the probe-MUT interface such that its effective impedance is modified.  This electrode polarization effect can be reliably modeled and corrected for by including a complex polarization impedance $Z_\mathrm{p}$ in series with the parallel combination of $C_\mathrm{f}$ and $C(\varepsilon_\mathrm{r},\sigma_\mathrm{dc})$ shown in Fig.~\ref{fig:fringe}~\cite{Schwan:1992, Raicu:1998}.  In our measurements of the NaCl solution, the effects of electrode polarization were clearly noticeable below 10~MHz.  The correction scheme reliably removed these effects allowing us to extend the bottom end of our measurement bandwidth by more than an order of magnitude.

\section{Waste-Activated Sludge Permittivity}\label{sec:WAS}

The open-ended coaxial probe was next used to determine the unknown permittivity of thickened WAS (TWAS) obtained from the Wastewater Treatment Facility (WWTF) in Kelowna, British Columbia, Canada.  As routine practice, before final disposal, TWAS samples are further dewatered at the WWTF using a centrifuge such that the solid content of the resulting material (called ``sludge cake'') is $18\pm 2\%$ by weight.  In this study, the permittivity and conductivity of both TWAS ($4.5\%$ solids) and sludge cake (18\% solids) were measured.  First, $\varepsilon_\mathrm{r}$ and $\sigma_\mathrm{dc}$ of the 18\% sample were determined from reflection coefficient measurements made using open-ended coaxial probes.  Next, distilled water was added to the 18\% sludge cake to produce a sample representative of TWAS, with a solid content of $4.5\pm0.5$\%, and the measurements of the electrical properties were repeated. We chose to dilute the sludge cake with distilled water so that we could directly compare data obtained from 18\%- and 4.5\%-solid solutions prepared from a single sample.  To be sure that the 4.5\% sample was a good representation of TWAS, we also measured a 4.5\% TWAS sample obtained directly from the Kelowna WWTF.  For all measurements, a 0.03~molarity NaCl solution was used as the liquid calibration standard as it mimics some of the essential electrical properties expected from the WAS samples.

When interpreting the results of the permittivity and conductivity of the WAS samples, comparisons were made to the results of previous studies of various biological tissues.  For a thorough review see Ref.~\cite{Pethig:1984} and references therein.  Figure~\ref{fig:WASpermRe} shows the real component of \mbox{$\varepsilon_\mathrm{r}-j\sigma_\mathrm{dc}/\omega\varepsilon_0$} for both the 4.5\% and 18\% TWAS samples after correcting for electrode polarization effects.
\begin{figure}[t]
\centerline{\includegraphics[width=0.7\columnwidth]{Fig2.pdf}}
\caption{Plot of $\varepsilon^\prime$ of WAS as a function of frequency at 25$^\circ$C.  Samples with both 4.5\% (black) and 18\% (blue) solid concentrations were measured.  The electrode polarization artifacts have been removed and the low-frequency upturns are features intrinsic to the samples.  At 4.5\%, results for both sludge cake diluted with distilled water (black points) and TWAS from the WWTF (open squares) are shown. The solid red curve is the real component of the permittivity of pure water.}
\label{fig:WASpermRe}
\end{figure}
For all TWAS samples, the reflection coefficient measurements were repeated at least five times.  In between measurements the probe tip was removed from the sample, cleaned, and then reinserted into the sample at a new location.  The reflection coefficients were very insensitive to the probe location indicating that the bulk material properties of the samples were measured.   Moreover, multiple TWAS samples were measured at both the 4.5\% and 18\% concentrations.  Consistent results were obtained among samples of the same concentration.  Note that $\varepsilon^\prime$ of the 4.5\% TWAS sample obtained from the WWTF is very similar to that of the sludge cake sample that has been diluted with distilled water.

Above 100~MHz, $\varepsilon^\prime$ of the 4.5\% sample is nearly identical to that of pure water.  This characteristic, referred to as $\gamma$-dispersion, is very common in biological samples and is due to the relaxation of polar water molecules.  Below 100~MHz, $\varepsilon_\mathrm{r}$ of the 4.5\% sample shows an upturn marking the onset of a distinct dispersion mechanism.  This effect, termed $\beta$-dispersion, is ubiquitous in samples of biological tissue and is typically active from 100~kHz to 10~MHz~\cite{Pethig:1984}.  The dispersion arises from the Maxwell-Wagner effect in which charge accumulates at the junction between materials of different permittivities and conductivities.  In biological samples, cell membranes form barriers between the intra- and extracellular fluids that in general have distinct electrical properties. The 18\% sample has a higher concentration of suspended cells and, as expected, the extracted $\varepsilon^\prime$ shown in Fig.~\ref{fig:WASpermRe} exhibits an enhanced $\beta$-dispersion below 100~MHz.

It is worth noting that $\beta$-dispersion may be a useful quantitative probe of cell viability in  WAS samples.  If a WAS sample is treated by a process that induces lysis, the number of intact cell membranes dividing the intra- and extracellular fluids would be reduced; presumably resulting in a suppression of $\beta$-dispersion.  It may be possible to use permittivity measurements in the range from 0.1 to 100~MHz as a quick and accurate initial indicator of treatment effectiveness.

In contrast to the 4.5\% sample, $\varepsilon^\prime$ of the 18\% sample shows no plateau marking the transition from one dispersion mechanism to the next.  Rather, the broad dispersive feature from 100~MHz to 10~GHz indicates that (i) the $\beta$-dispersion associated the charging of cell membranes extends to higher frequencies and (ii) there are other dispersion mechanisms active in this frequency range.  These features, generally termed $\delta$-dispersion, are attributed to a variety of mechanisms.  Most commonly $\delta$-dispersion is thought to arise from dielectric relaxation of water molecules bound to adjacent proteins.  The rotational relaxation time of bound water molecules is longer than that of bulk water causing the dielectric dispersion to be shifted to lower frequency.  Much of the bound water will be associated with proteins present in cell membranes~\cite{Pethig:1984}.  It is plausible that exposing WAS to electromagnetic radiation in the UHF range can bring about physiological changes that may enhance anaerobic digestion.  For example, if tuned to the relaxation rate of water molecules bound to membrane proteins, the electromagnetic energy will be, in part, absorbed directly by the intracellular fluid which may lead to cell lysis after sufficient exposure.

Focusing now on Fig.~\ref{fig:WASpermIm}, $\varepsilon^{\prime\!\prime}+\sigma_\mathrm{dc}/\omega\varepsilon_0$ of both samples is dominated by the conductivity term below 1~GHz.
\begin{figure}[t]
\centerline{\includegraphics[width=0.7\columnwidth]{Fig3.pdf}}
\caption{Plot of $\varepsilon^{\prime\!\prime}+\sigma_\mathrm{dc}/\omega\varepsilon_0$ of WAS as a function of frequency at 25$^\circ$C.
Electrode polarization corrections have negligible effect on these data.  The dashed lines have slope -1 and represent $\sigma_\mathrm{dc}/\omega\varepsilon_0$ contributions.  At 4.5\%, results for both sludge cake diluted with distilled water (black points) and TWAS from the WWTF (open squares) are shown.  The solid curve is $\varepsilon^{\prime\!\prime}$ of pure water.}
\label{fig:WASpermIm}
\end{figure}  The dashed lines in the figure follow a $\omega^{-1}$ frequency dependence and are drawn using $\sigma_\mathrm{dc}=0.34~\Omega^{-1}\mathrm{m}^{-1}$ and $\sigma_\mathrm{dc}=0.68~\Omega^{-1}\mathrm{m}^{-1}$ for the 4.5\% and 18\% samples respectively.  For frequencies below 100~MHz, linear fits to $\log\left(\varepsilon^{\prime\!\prime}+\sigma_\mathrm{dc}/\omega\varepsilon_0\right)$ as a function of $\log\omega$ were performed to extract $q$ from the observed $\omega^{-q}$ behavior.  For the 4.5\% sample, $q$ was $0.9961\pm 0.003$.  A value so close to one implies that $\varepsilon_\mathrm{r}$ exhibits very little dispersion over this frequency range, consistent with the observed $\varepsilon^\prime$ behavior. For the 18\% sample, on the other hand, with a value of $0.9728\pm 0.0003$, $q$ was noticeably less than one.  The only plausible explanation being that, in the higher concentration sample, dispersive mechanisms are active over this frequency range.  Above a few gigahertz, the conductivity term is negligible and the data are dominated by $\varepsilon^{\prime\!\prime}$.  As expected, the 4.5\% sample tracks the $\gamma$-dispersion of pure water very closely, whereas the 18\% sample shows significant deviations due to an active $\delta$-dispersion.  Again, the sludge cake sample that has been diluted with distilled water mimics the behavior of the TWAS sample from the WWTF.

\section{Summary \& Discussion}\label{sec:conclusions}

To the best of our knowledge, the electrical properties of WAS have not been previously reported.  We have preformed this characterization at two different concentrations of solid content over a broad range of frequencies and have identified the mechanisms leading to the observed dispersions.  The experimental and theoretical insights gained can be used to guide optimizations of the electromagnetic pretreatment of WAS.

Enhancements of the methane production rate are most likely to arise from punctured or otherwise compromised cell membranes that increase the bioavailability of the WAS to methane-forming microorganisms present during digestion.  At very high frequencies (tens of gigahertz) the WAS samples absorb energy from the electromagnetic fields via rotations of unbound dipolar water molecules.  This energy, dissipated as heat, will diffuse throughout the rest of the sample by conduction.  It is difficult to imagine how this microwave treatment would differ from a thermal cycling of WAS achieved, for example, using a hotplate.

The $\delta$-dispersion observed in the 18\% sample will be sensitive to electromagnetic treatments, between 100~MHz and 10~GHz.  Relaxation of protein-bound water molecules attached to cell membranes could be effective at initially concentrating the dissipated energy at the sites of the cells and could potentially induce cell lysis.  The reported benefits seen in anaerobic digestion when using WAS pretreated with electromagnetic irradiation at 915~MHz and 2.54~GHz may be in part a consequence of this effect.  Radio frequency treatments between hundreds of kilohertz and tens of megahertz are particularly intriguing as they will activate the Maxwell-Wagner effect causing the cell membranes to be repeatedly charged and discharged.  This process will also concentrate the initial dissipated energy at the cell sites and may be more effective at inducing cell lysis as the observed $\beta$-dispersion is significantly stronger than the $\delta$-dispersion.  Due to the higher concentration of suspended cells, $\beta$-dispersion is enhanced in the 18\% sample.  This observation suggests that electromagnetic pretreatments at radio frequencies are likely to be most effective when applied to samples having a high solid content.

As an outcome of this study, we are currently designing and building a series of custom high-power electromagnetic applicators that span a broad range of frequencies.  For each applicator, our permittivity and conductivity measurements allow matching circuits to be designed so as to achieve maximum power transfer from the source to the WAS load.  The overall efficiency of any practical large-scale applicator will be of critical importance.  The effectiveness of each applicator will be evaluated by quantifying cell lysis that results from the electromagnetic pretreatments and by measuring enhancements of the methane production rates from pretreated WAS during anaerobic digestion.

\ack

The financial support of the Natural Science and Engineering Research Council of
Canada Strategic Project Grant (\#396519-10) is gratefully acknowledged.  The Agilent 8722ES VNA used in this work was provided by CMC Microsystems and their support is gratefully acknowledged.