\section*{}
\vspace{-1cm}

\section{Statement of significance}

Our manuscript shows how to calibrate the thermal motion of a probe trapped in-vivo using optical tweezers. We use a new theory that makes the viscosity frequency dependent quite akin to the Jeffery's model. We find that the fits to the power spectral density are very good and we can extract the 0 Hz (DC) viscosity that match well with literature. We can make the optical tweezers calibration "at a glance". In view of the fact that calibration of optical tweezers in-vivo has remained challenging due to the variability of the intracellular viscosity and refractive index, this fit to the Jeffery's model solves the problem completely.

\section{Introduction}
Microrheology is the study of flow of matter at micron length scales and microliter sample volumes. It is particularly important when rare or precious materials are involved, like in bio-physical studies. These can be performed in an environment where conventional rheological tools cannot reach inside the  cells\cite{fabry,allan}. Rheological characterization of cells is relevant for early diagnosis of diseases like malaria \cite{malaria} and migration of cancerous cells \cite{cancer1}.  The Red Blood Cells (RBC) in malaria infected patients are known to have different elasticity than healthy ones. The Extra-Cellular Matrix (ECM) of a cell in the human body interacts and changes the stiffness of the cytoskeleton in healthy cells via mechanical adaptation, the exact mechanism of which is not very well understood. A better match of the rheological features of the ECM and the cell is required for their adherence and cell mobility. Normally, when the ECM is too stiff and the cell cannot change its cytoskeletal features, there is no adherence. In the case of cancer cells, this mechanism is suppressed thereby allowing attachment inspite of stiffness mismatch and executes metastasis. Thus, one promising strategy to address cancer could be to understand and revive this cellular mechanism to match rheological features. In view of these facets, in-vivo rheology is useful. Further, intracellular viscoelasticity has a role in diffusion of molecules and performance of chemical reactions, not to mention, such rheology also enables calibration of externally applied forces and torques on the system.

 Recently, the cell cytoplasm was shown to be a poroelastic medium, \cite{poroelastic} implying a medium with an elastic meshlike network into which a viscoelastic medium is suspended. There have been numerous attempts to ascertain both the viscosity\cite{puchkov,kuimova} and the viscoelasticity\cite{berret} of the cell cytoplasm. However, it is well known that the viscosity results are prone to the size of the probe used \cite{sizeeff}. Any probe larger than 100 nm shall yield a viscosity larger than the native one, when probed actively. It is here, that a passive detection technique that moves the probe by small amounts can be expected to avoid the elastic mesh and provide information about the cytoplasmic fluid. Thus we study the passive thermal fluctuations of a probe using optical tweezers. Attempts to use optical tweezers in-vivo \cite{oddershede,roop,goldman,reihani} have been complicated due to the problems in quantifying the optical trap stiffness accurately \cite{gross}, mainly due to the variable nature of the intracellular refractive index and the viscosity. We use a new theory to directly ascertain the unknown parameters from the fit  "at a glance".

 It is generally believed that the viscoelasticity of a medium automatically implies a segment-wise power-law behavior of the complex elastic modulus ($G^*(\omega)$ = A $\omega^{\beta}$) in the frequency domain \cite{powerlaw}.
 Recently, there has been an attempt to establish the macroscopic viscoelasticity of the medium from the microscopic Stokes Oldroyd-B model that showed a frequency-dependent viscosity quite akin to the Jeffery's model \cite{shuvojit}. We use this new strategy to attack the problem of intra-cellular viscoelasticity and find that the motional power spectral density for
 spherical polystyrene particles of radius $ a_0 = 0.5 \mu m$ fits well to the suggested model. We extract the 0 Hz (DC) viscosity and find a good match with the established intra-cellular viscosity values. We also find the relaxation time of the intracellular medium and the optical trap stiffness also emerge from the fitting parameters.

\section{Theory}

The frequency- dependent viscosity in an incompressible low Reynolds number viscoelastic medium comprising of a solvent and a polymer solute dissolved in it has been found to be given by the following expression derived from the Stokes Oldroyd-B model for linear microscopic viscoelasticity \cite{shuvojit}
\begin{equation}
     \mu(\omega)=\mu_s + \frac{\mu_p}{-i\omega\lambda+1}
    \label{oldroyd}
\end{equation}
Where $\mu_s$ is the zero frequency solvent viscosity, $\mu_p$ is the zero frequency polymer viscosity and $\lambda$ is the polymer relaxation time. This expression is very similar to the Jeffery's model of frequency dependent viscosity with the coefficients labelled differently.  Thus the viscosity of the solution at zero frequency would be, $\mu_0$=$\mu_s+\mu_p$.
Solving for the power spectral density of an optically trapped particle in a viscoelastic fluid, we get \cite{shuvojit}
\begin{equation}
 <x(\omega)x^*(\omega)>= \frac{2k_BT}{\gamma_0}\frac{(\frac{(1+\frac{\mu_p}{\mu_s})}{\lambda^2}+\omega^2)}{[(\frac{\kappa}{\gamma_0\lambda}-\omega^2)^2+\omega^2(\frac{\kappa}{\gamma_0}+\frac{1}{\lambda}(1+\frac{\mu_p}{\mu_s}))^2]}
 \label{PSD}
\end{equation}
The term $\kappa$ signifies the trap stiffness and $\gamma_0$ is the drag coefficient for only the solvent. The exact expression for the $\gamma_0$ is given by,

\begin{equation}
\gamma_0=6\pi\mu_s a_0
\label{solvent}
\end{equation}

In order to correlate with experimentally obtained power spectral density curves of translational motion, we rewrite the eq. \ref{PSD}, as

\begin{equation}
 <x(\omega)x^*(\omega)>= \beta^2 A\frac{(\frac{(1+\frac{\mu_p}{\mu_s})}{\lambda^2}+\omega^2)}{[(\frac{\kappa}{\gamma_0\lambda}-\omega^2)^2+\omega^2(\frac{\kappa}{\gamma_0}+\frac{1}{\lambda}(1+\frac{\mu_p}{\mu_s}))^2]}
 \label{fit}
\end{equation}

where the A coefficient indicates the amplitude in terms of Volts$^2$/Hz and the calibration factor is $\beta$ in (m/Volt) quite akin to the conventional calibration factor for normal media \cite{erik}. The calibration factor for the translational signal is related to temperature as

\begin{equation}
 \beta^2 A = \frac{2 k_B T}{\gamma_0}
 \label{fit1}
\end{equation}
Thus, the calibration factor $\beta$ is given as

\begin{equation}
 \beta = \sqrt{\frac{2 k_B T}{A \gamma_0}}
 \label{fit2}
\end{equation}
Fitting the power spectral density with this equation we can extract the values of relative viscosity($\frac{\mu_s+\mu_p}{\mu_s}$) of solution and polymer relaxation time constant($\lambda$).

\section{Experimental details}

To perform the experiment, we place a polystyrene particle (1  $\mu$m diameter) inside a cell attached to the glass slide (Blue Star, 75 mm length, 25 mm width and a thickness of 1.1 mm) of a sample chamber assembled, as shown in Fig. \ref{schematic}. The other side of the sample chamber is formed by a cover slip (Blue Star, number 1 size, english glass). The cell is close to the top surface of the sample chamber and illuminated in an inverted microscopy configuration using the Optical Tweezers kit (OTKB-M, Thorlabs USA) \cite{rahul}. The illumination objective is a 1.3 NA, 100x oil immersion objective from Olympus at the bottom with the illumination aperture being overfilled. The collection objective at the top is  E Plan 10x, 0.25 N.A. air-immersion objective from Nikon. The laser used for optical trapping is a diode laser from Lasever at 1064 nm wavelength which typically has a maximum power of 1.7 Watt with about 400 mW power in the sample plane. However, for the experiments performed in this manuscript, the laser was set to 200 mW power at the sample plane. An LED lamp illuminates the sample, coupled in using a dichroic mirror, as shown in Fig. \ref{schematic}. Another dichroic mirror couples the broadband visible light out of the path of the laser and illuminates a CMOS camera (Thorlabs, USA).

\begin{figure}[h]
\centering
  \includegraphics[width=\linewidth]{schematic.png}
  \caption{The schematic of the experimental set up. The MCF7 cell is attached to the top surface of a sample chamber with the probe particle located inside the cell. The Optical tweezers light traps the particle and then the forward scattered light is used to perform translational measurements.}
  \label{schematic}
\end{figure}

The translational motion is recorded by impinging the forward scattered light onto a Quadrant Photodiode (QPD) (Thorlabs, USA). This has four quadrants indicated by A, B, C, and D in Fig. \ref{schematic}. The top half minus the bottom half yields y displacement while the right half minus the left half yields the x displacement \cite{basudev}. The bandwidth of this detector is 40 KHz. The position signals emerging from the QPD are acquired by the computer using data acquisition cards (National Instruments, NI PCI 6143) which has a bandwidth of 40 KHz.

We first tested that our technique works for a solution of polymers in water. For this, a concentrated solution of polyacrylamide (PAM, 1$\%$ by weight) was prepared in a water solvent with suspended 1 $\mu$m diameter polystyrene particles to make the viscoelastic solution \cite{shuvojit2}. One such particle was trapped, and the power spectral density of translational motion recorded with the optical tweezers system. The curve has been mentioned in Fig. \ref{polymer} and fits well to the x-direction power spectral density dataset.

\begin{figure}
 \centering
 \includegraphics[width=\linewidth]{viscoelastic1.pdf}
 \caption{ This figure demonstrates a typical fit using eq. (\ref{fit}) to the case of a viscoelastic solution made from PAM in water. }
 \label{polymer}
\end{figure}

In order to create in-vivo conditions, MCF-7 (Michigan Cancer Foundation-7) cells were grown on glass slides coated with gelatin. These glass slides were first treated with piranha solution and sterilized using UV (265nm) light for 20 minutes and coated with $ 0.5\% $ gelatin solution. MCF7 cells were added on to the center of the coverslip and Dulbecco's Modified Eagle Medium (DMEM) supplemented with 10$\%$ fetal bovine serum and 1$\%$ glutamine-penicillin-streptomycin was added on top of the coverslip. 10$\mu$L of gelatin-coated 0.5$\mu$m radius polystyrene particles($1\mu$g/mL) suspended in sterile serum-free media were sonicated and added to the cells. Cells were incubated with
5\% carbondioxide and 37$^{\circ}$C for 12 hours to initiate endocytosis \cite{guo}.

 The temperature of the system was maintained using an air conditioner in the room at 26$^{\circ}$C. Care was taken to ensure that the cool air from the air-conditioner was not blowing directly onto the system.

\begin{figure}
 \centering
 \includegraphics[width=\linewidth]{in_vivo_psd.png}
 \caption{This figure represents the power spectral density (a) of x-displacements of a tracer particle optically trapped inside the cell. The data is well fitted by a viscoelastic equation of the form eq. \ref{fit}. The inset shows an image of the tracer particle of 1 $\mu$m diameter trapped inside the cell. (b) The residuals to the fit.}
 \label{invivo}
\end{figure}

The particles incubated with the cells might either be outside the cell or inside. We ascertain this by first trapping the particle and then moving the stage back and forth manually in the X direction using micrometer screws by about 10 $\mu$m. If the particle is outside the cell and not residing on the membrane, it shall readily follow the trap without too much delay. If the particle is residing on the membrane and also non-specifically bound, it cannot follow the trap much. However, if the particle is indeed inside the cell, it follows the trap within the space available to it without hitting the side membrane or organelles like the nucleus. One such video has been shown in S1. Once we ascertain the focus of the objective for which the particles can be trapped inside the cell, we mark it to find other particles at the same depth.

\section{Results and discussions}

A typical motional power spectral density for a particle inside a cell has been shown in Fig. \ref{invivo}. Every power spectra is calculated by taking a 5 second time series and averaging over 10 such spectra. The eq. (\ref{fit}) fits well to the experimental data. The fitting parameters automatically yield the A parameter in Volts$^2$/Hz and subsequently the calibration factor $\beta$, as shown in eq, (\ref{fit2}). We also extract the ratios $\frac{\kappa}{\gamma_0}$ and $1+\frac{\mu_p}{\mu_s}$. We show two more datasets indicating PSD's at two different locations of the same cell in fig. \ref{sample}. The green curve shows a typical DC viscosity of 5 times that of solvent while the blue curve shows a value of 13 times that of solvent. There can be variations by  factors of 2 to 3 inside the same cell. These spectra are also calibrated to indicate thermal motion at each frequency. We can see that the DC thermal motion is of the order of 100 nm which reduces to 3 nm at 1000 Hz. Thus, we typically achieve a frequency range of about 3000 Hz with our present configuration.

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{visco1.pdf}
    \caption{Power spectral density of the trapped particle fitted with the Oldroyd model, at different regimes of the same cell. The green curve indicates a typical PSD with DC viscosity about 5 times that of solvent while the blue curve indicates another location where it is 13 times that of solvent.  We show that there are regions showing higher DC viscosity inside the same cell by more than a factor of 2.}
    \label{sample}
\end{figure}

 If we assume that the solvent medium for the cell is water, which has been proved to be a good approximation \cite{solvent}, the $\mu_s$ is automatically that of water while the $\gamma_0$ is given by eq. (\ref{solvent}). Using this approximation, we extract values for the DC viscosity $\mu_0 + \mu_p$, indicated in fig. \ref{collated}(a) and $\lambda$, indicated in fig. \ref{collated}(b) for 58 events.

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{collated.png}
    \caption{(a) The value of DC(0Hz) viscosity for the cell cytoplasm for 58 number of events.   (b)Measurement of polymer time constant $\lambda$ for the passive motion of the probe particle in every 58 events. (c) This shows  typical trap stiffness of 1 $\mu$m diameter polystyrene particles in-vivo. }
    \label{collated}
\end{figure}

The DC (0 Hz) viscosity for the cell cytoplasm seems to vary from 2 mPa sec to 16 mPa sec, which is what one expects for typical intracellular concentrations of protein filaments like actin, microtubules and intermediate filaments dissolved in water\cite{proteinconc}. Further, the average relaxation time is 0.021 $\pm$ 0.005 sec. A typical estimate using the Maxwell's model observed a value of 0.2 sec\cite{magnetic}.

We also extract the typical trap stiffness at the same value of laser power of 200 mW at the sample plane, shown in fig. \ref{collated}(c). The trap stiffness seems to have an average value of 0.1 pN/nm.

In these measurements, we span a frequency range between 2 Hz and 3000 Hz, partly including the athermal fluctuations inside the cell ranging from 0.1 Hz to 10 Hz\cite{schmidt}. Indeed we find deviations from the thermal curve below 5 Hz, as shown in Figs. \ref{invivo} and \ref{sample}. Moreover, as indicated by Tassieri \cite{tassieri}, the tweezers can indeed probe the viscoelasticity of the cell if the measurement time is lower than active motion time, given by the Deborah number (De) being greater than 1. Since we have a frequency range extending to 3 KHz, the measurmeent time is indeed smaller than 5 Hz activity time, thus enabling this kind of approach. Once the optical tweezers is calibrated properly, newer experiments could even be designed to look specifically at the athermal fluctuations.

We also calibrate a typical time series for the x displacement of a 1 $\mu$m diameter microsphere inside the cell, while being carried by molecular motors on microtubules. It has been shown in Fig. \ref{time_series}. The red curve shows the unfiltered data while the black curve shows a median filtered data. The optical trap stiffness in this case was 0.12 pN/nm. Thus the force applied by the molecular motors on the particle is then about 10 pN, possibly indicating multiple motors at work.

 \begin{figure}
    \centering
    \includegraphics[width=\linewidth]{time_series.pdf}
    \caption{This figure indicates a calibrated x-position time series of a 1 $\mu$m diameter particle moving on a microtubule inside the cell. The red line indicates unfiltered data while the black curve indicated median filtered data.}
    \label{time_series}
\end{figure}

\section{Conclusions}
In conclusion, we could trap 1$\mu$m diameter polystyrene particle inside MCF-7 cells and obtain the power spectral density for passive motion of the particle along the X-axis inside the cell. To fit the PSD, we model the cytoplasm as a polymer network immersed in water that acts as a viscoelastic medium. This power spectrum is fitted with Jefferey's model. We see that the relative viscosity, polymer relaxation time and the trap stiffness varies inside the cell from place to place possibly due to the variation of cytoplasmic density, but maintains a good agreement with previous literature. This study of the viscoelasticity of the cell is made possible by the enhanced frequency response due to the power spectral density extending to 3 KHz, thereby ensuring a Deborah number greater than 1, when compared with activity frequency of 5 Hz or slower.

\section{Author contributions}

B.R. and S.B. designed the experiment. R. V., S.R.V. and B.R. performed the experiment and did data analysis. B. R. and S. B. wrote the manuscript.

\section*{Acknowledgements}
We thank the Indian Institute of Technology Madras, India for their seed grant.


\subsection{Simulations}
Due to the difficulty of accurately modeling the complex turbulent wind effects on a quadrotor, our simulations consider steady wind but serve to test our finite state machine, planner, controller, and estimator. For space constraints, simulation experiments are not analyzed here, but the video accompanying the paper (linked in the Supplementary Material Section) shows a landing simulation.

\begin{figure} [t]
	\begin{center}
	\vspace*{.1in}
		\includegraphics[trim =0mm 0mm 0mm 0mm, clip, width=0.35\textwidth]{figs/sq01-reduced.png}
\vspace*{-.05in}
\caption{Quadrotor vehicle used for the experiments.}
		\label{fig:sq01}
	\end{center}
\begin{center}
\includegraphics[trim =10mm 150mm 0mm 0mm, clip, width=0.35\textwidth]{figs/lbarray-reduced.jpg}
\vspace*{-.05in}
\caption{Leaf blower array used in front of the static platform showing the wind measurement procedure.}
\label{fig:lbarray}
\end{center}
\iffalse
\begin{figure}[t]
\begin{center}
\includegraphics[trim=200mm 155mm 0mm 90mm, clip, width=0.35\textwidth]{figs/jackal-dolly-1.JPG}
\vspace*{-.05in}
\caption{Landing platform on dolly towed by UGV.}
\label{fig:jackal-dolly}
\end{center}
\end{figure}
\fi
\begin{center}
\includegraphics[trim =0mm 0mm 0mm 0mm, clip, width=0.35\textwidth]{figs/windsp-eps-converted-to.pdf}
\vspace*{-.05in}
\caption{Mean wind speed (in m/s) from the leaf blowers vs.~distance to the vertical platform (with $1-\sigma$ error bars).}
\label{fig:windsp}
\end{center}
\end{figure}

\begin{figure*}[ht]
\begin{center}
\vspace*{.05in}
\includegraphics[trim =40mm 5mm 40mm 0mm, clip, width=1\textwidth]{figs/static_sept11_no-eps-converted-to.pdf}
\vspace*{-.175in}
\caption{Position and velocity tracking and estimation performance of a static platform experiment using a BLSC that does not model the turbulence. Wind blowing towards $-x$ resulting in poor tracking performance in the $x$-direction. Vertical dashed lines indicate time of the first tag detection, $t_d$.}
\label{fig:static_no}
\end{center}
\end{figure*}

\subsection{Hardware}
Static and moving platform landing tests were done (see Figs.~\ref{fig:video-screenshot} and \ref{fig:lbarray
), both of which had turbulent wind at the landing site from the leaf blowers. The quadrotor used for the hardware experiments (Fig.~\ref{fig:sq01}) weighs 0.564kg including the 1500mAh 3S battery. It is $36\times29$cm (approximately half the platform size) and its thrust-to-weight ratio is 1.75. The onboard computer is a Qualcomm Snapdragon Flight APQ8074, whose front-facing camera provides $640\times480$ black-and-white images at a rate of 30fps to the AprilTag detection module.
Hover tests were carried out to determine $\hat{b}$ in \eqref{eq:u}.
To measure the drag coefficient $\hat{c}$ and the bound on its error $\tilde{c}$, we performed tests with the quadrotor flying in front of a strong wind and the accurate pitch angle was obtained by the motion capture system.
By balancing forces, the drag could be determined, yielding $\hat{c}$.
Additionally, we used the {\tt IMU utils} package \cite{imucal} to compute the Snapdragon's IMU accelerometer and gyroscope noise density and bias random walk.
The {\tt Kalibr} visual-inertial calibration toolbox then used this IMU intrinsic information to find the camera-IMU transform \cite{furgale2015kalibr}.

\subsection{Static Platform Experiments}
\subsubsection{Experiment Setting}
The first set of hardware experiments presents a static platform in front of an array of 5 leaf blowers, as shown in Fig.~\ref{fig:lbarray}. The leaf blowers are set at two different heights and point to the negative $x$-direction.
The platform is composed of a 60$\times$60cm horizontal plate and a 60$\times$30cm vertical plate, which has attached the tag bundle described in \ref{subs:detection}.

The parameters $v_w$ and $\sigma$ of the turbulent wind were measured at distances $l$ from the vertical platform every $0.5$m until $l = 3.5$m, which is approximately the tag detection range. For every $l$, a measurement of the speed was taken every second for $60$s using a high-precision hot-wire anemometer.
Fig.~\ref{fig:windsp} shows the data obtained. Interestingly, at $l = 0.5$m, the mean speed decreases while the standard deviation increases to its maximum value. We believe this is due to the vortices caused when the flow traverses the vertical platform.

\begin{figure*}[ht]
\begin{center}
\includegraphics[trim =40mm 5mm 40mm 0mm, clip, width=1\textwidth]{figs/static_sept11_yes-eps-converted-to.pdf}
\vspace*{-.175in}
\caption{Position and velocity tracking and estimation performance of a static platform experiment using our BLSC. The tracking error remains small during the flight.}
\label{fig:static_yes}
\end{center}
	\begin{center}
		\includegraphics[trim =40mm 5mm 40mm 0mm, clip, width=1\textwidth]{figs/moving_sept12_yes-eps-converted-to.pdf}
		\vspace*{-.175in}
\caption{Position and velocity tracking and estimation performance of a moving platform experiment using our BLSC. The quadrotor quickly matches the UGV's velocity and is able to track successfully the planned trajectory.}
		\label{fig:moving_yes}
	\end{center}
\end{figure*}

\subsubsection{Results}
First, we tested a standard BLSC that does not take into account the turbulent wind generated by the leaf blowers, that is, the factors $\hat{f}$ and $K$ only consider the drag generated by the quadrotor's speed relative to the ground. Fig.~\ref{fig:static_no} shows the tracking and estimation performance obtained. The quadrotor starts at $\boldsymbol{p_q}=(-3.5, 2.0, 1.3)$m, and the platform is located at $\boldsymbol{p_p}=(1.1, -1.1, 0.8)$m. As expected, the tracking performance is poor, and the landing takes a long time: $27.2$s since the first tag detection, which occurs at $t_d=3.3$s (indicated with a vertical dashed line). Also note that the platform's estimation improves after $t_d$. In the video, the covariance ellipses for the tag's 2D position can be seen, and decrease abruptly in size at $t_d$.

Next, we used our BLSC with the same initial conditions as in Fig.~\ref{fig:static_no} to compare its improvement. The results are shown in Fig.~\ref{fig:static_yes}. It can be seen that the tracking is much better, and it only worsens slightly when the quadrotor is inside the wind field, after the time of the first tag detection $t=2.5$s. The landing time is just $3.2$s (measured since the first tag detection) even in challenging conditions, which is 8.5 times faster compared to the standard BLSC.

\subsection{Moving Platform Experiments}
\subsubsection{Experiment Setting}
To fully test our approach, we also performed landing experiments on a moving platform that is mounted on top of a dolly. The \textit{Clearpath Jackal} was used as the ground vehicle that tows the landing platform (see Fig.~\ref{fig:video-screenshot}) and carries two of the leaf blowers. This provided turbulent air at the landing pad even though the vehicles were not moving very quickly. The distance from the leaf blowers to the platform is such that this turbulent wind follows the same plot as Fig.~\ref{fig:windsp}.

\subsubsection{Results}
A standard BLSC was also tested first for this experiment setting. There was considerable tracking error and the quadrotor was not able to land on the platform before the vehicle arrived at the final point (see video linked in the Supplementary Material Section for details).
The results obtained using our BLSC are shown in Fig.~\ref{fig:moving_yes}. The quadrotor starts at $\boldsymbol{p_q}=(-1.5, 2.7, 0.9)$m, and the UGV starts at $\boldsymbol{p_p}=(-3.9, 0.1, 0.4)$m. When the maneuver begins, the UGV is commanded to move at 1m/s and rotate to its right at a rate of 4$^\circ$/s. The tag is detected at $t_d = 4.7$s, and the landing time is $6.8s$. Note that the quadrotor is at approximately $4m$ from the moving platform at $t_d$, a value larger than for the static experiment (1.5m).

\subsection{Related Work}

Autonomous UAV landing has been investigated by several researchers.
Ref.~\cite{venugopalan2012autonomous} landed a drone on a static kayak in a reservoir with mild wind and water ripple conditions, but the image processing was done off-board and the landing time was close to 1~min.
Ref.~\cite{araar2017vision} developed a system capable of landing a UAV on a moving platform with PID control, an EKF estimator, and an AprilTag visual fiducial bundle.
However, the computation was also done off-board, and there was no relative wind present: the platform's speed was just 0.18m/s and the tests were done indoors.
Ref.~\cite{falanga2017vision} demonstrated a quadrotor landing on a moving platform using only onboard sensing and computation, but the environment was not turbulent due to the tests being indoors and the platform was moving relatively slowly at 1.2m/s.
Furthermore, besides two cameras, a distance sensor was required to estimate the UAV-platform relative pose.
Ref.~\cite{borowczyk2017autonomous} demonstrated an autonomous landing of a quadrotor on a car moving at 14m/s by detecting an AprilTag on its roof.
A proportional navigation-based guidance law was used for the approach phase, and a PID controller for the landing phase, with no disturbance rejection considerations. The landing maneuver consisted of acquiring the tag while hovering, and the descent was initiated once the quadrotor stabilized over it. Additionally, the quadrotor used was large and had a broad sensor suite, including a downward-facing camera, a three-axis gimballed camera to track the target, and an inertial navigation system. Furthermore, besides the ground vehicle's GPS coordinates, the quadrotor also used the ground vehicle's IMU data to improve its estimation of the landing platform's state.
Ref.~\cite{xing2019autonomous} also demonstrated an outdoors landing, but the landing platform's speed was only 0.5m/s and the landing maneuver took 12--20~s to complete. Moreover, the ground vehicle's wheel encoders data was used to estimate its state.
Ref.~\cite{persson2019model} used model predictive control (MPC) to land a quadrotor on a moving ship. Both the UAV and the ship collaborated to reach their goal.
The waves were modeled as sinusoidal, but the wind disturbances considered were not turbulent -- this approach only compensated the effects caused by a steady-state wind. Furthermore, the vehicle needed to hover above the platform for at least 5~s, and the ship was simulated.
A simulated boat landing was also carried out in \cite{wynn2019visual}, where the platform's state was estimated fusing GPS and visual measurements.
This work incorporated a velocity feed-forward term to the controller to catch the platform, but its only consideration of external wind was an offset of the hovering position to ensure the target was inside the field of view of the downwards-facing camera, and the landing maneuver lasted 24~s.
Ref.~\cite{bhargavapuri2019vision} developed an adaptive controller to track a ground vehicle with only relative position data from ArUco tags, and tested it in outdoor experiments at 5.6m/s. But this approach only considered disturbances due to ground effect, and the landing maneuver needed 20~s for tracking and 10~s for descent.

In summary, most current approaches for quadrotor landing on moving platforms involve hovering above the platform for a period of time to visually acquire it, which is then followed by a relatively slow descent.
In contrast, our approach investigates a direct trajectory to the landing platform, which has the possible advantages of providing faster landings and enabling several UAVs to approach the platform at the same time from different directions (improving the utilization of this limited resource).
Additionally, most of current approaches do not make special considerations to reject the turbulent wind present near the target vehicle, and thus safety in challenging conditions cannot be guaranteed.

\subsection{Contributions}

 This paper demonstrates a vision-based system capable of dynamic landing (i.e., the multirotor does not need to hover above the vehicle before descending) which also accounts for the turbulent conditions that would be present near a rapidly-moving ground vehicle. The resulting framework allows thus a maneuver that will be crucial for efficient truck-drone delivery systems. The contributions of this paper are as follows:
\begin{itemize}

    \item Demonstration of a boundary layer sliding controller to incorporate and compensate for turbulence based on a model of the conditions near the landing platform.
    \item An algorithm for computing fast, vision-based dynamic landing maneuvers, is demonstrated both in simulation and hardware experiments that include challenging steady/turbulent wind conditions.
\end{itemize}

\iffalse
\subsection{Outline}
The remainder of this paper is organized as follows: Section~\ref{s:system} presents an overview of the system and describes its main components. Section~\ref{s:exp} shows the experimental results obtained both in simulation and in real life experiments. Finally, Section~\ref{s:conclusion} concludes this paper and proposes directions for future work.
\fi

\section*{Supplementary Material}

Video of the paper summary and experiments is available at \url{https://youtu.be/xKo1rY4riJQ}.

\section{INTRODUCTION}\label{s:intro}
\input{introduction.tex}

\section{SYSTEM OVERVIEW}\label{s:system}
\input{system_overview.tex}

\section{EXPERIMENTAL RESULTS}\label{s:exp}
\input{experimental_results.tex}

\section{CONCLUSION}\label{s:conclusion}
\input{conclusion.tex}

\section*{ACKNOWLEDGMENT}

Work supported by Ford Motor Company, and the authors would like to thank Michael Everett for his help configuring the UGV.

\balance

\bibliographystyle{IEEEtran}

\subsection{Finite State Machine}\label{subs:fsm}
The quadrotor's behavior is determined by a finite state machine (FSM) comprised of four states:
\begin{enumerate}
    \item \textit{Stand By}: This is the initial state of the quadrotor, which consists of taking off and hovering at a predefined altitude above the starting point. The FSM then transitions to \textit{Search} mode.
    \item \textit{Search}: The quadrotor uses simulated GPS coordinates of the unmanned ground vehicle (UGV) ---as explained in Section~\ref{subs:estimation}--- to predict a rendezvous location and flies there. When the front-facing camera detects the landing platform as described in Section~\ref{subs:detection}, the state automatically switches to \textit{Landing}.
    \item \textit{Landing}: In this mode, the quadrotor approaches the target following a direct trajectory towards it. When the distance and relative velocity UAV-UGV are below threshold values, the quadrotor switches to the \textit{End} mode. If the last detection happened more than 0.8s ago, the mode returns to \textit{Search}.
    \item \textit{End}: motors are stopped, maneuver is finished.
\end{enumerate}
\subsection{Trajectory Planning: Model Predictive Control}\label{subs:planning}

The planner solves a convex optimization problem with changing objectives depending on the state of the FSM.
In \textit{Search} mode, the UAV-UGV rendezvous point is predicted by assuming a constant linear velocity and yaw rate for the UGV, and a constant velocity for the UAV.
This position is then offset by a small distance backwards in the direction of the UGV to ensure target detection by the front-facing camera, and is used by the planner as the final position of the trajectory, while the final velocity is the UGV's.
The time taken to reach the UGV is minimized, to reduce delivery turnarounds.
In \textit{Landing} mode, the planner initially minimizes the jerk to obtain a trajectory which ensures adequate tag acquisition.
As the UAV approaches the target, the effect of disturbances increases so the planner minimizes the time spent in the turbulent area. The final position of the planned trajectory is the vertical tag's position (offset a few cm backwards) and predicted ahead (by an amount that depends on the computation time) assuming constant linear/angular velocities.
The final velocity of the trajectory is set to match the UGV's.

The minimum jerk approach produces smooth trajectories and has a long heritage for planning quadrotor paths \cite{yu2014minimum, phang2015systems, falanga2017vision}.
The trajectory is re-planned using an MPC approach every time a new estimate of the platform's position is obtained. We used CVXGEN \cite{MattingleySB12} to generate the code to solve the following convex optimization problem
\begin{equation}
\begin{aligned}\label{eq:opt_problem}
    \min_{x,u}  J &= \sum_{t=0}^{N} \ell \left( \bar{u}, d \right)[t] \\
    \text { subject to } x[t+1] &= Ax[t] + B\bar{u}[t] \\
    |a[t] |_{\infty} & \leq a_{\max} , \text{ and } |\bar{u}[t]|_{\infty} \leq j_{\max}   \\
    x[0] &=x_{0}, \, x[N]=x_f  \\
    \text {for } t &= 0\ldots N
\end{aligned}
\end{equation}
where $N$ is the timestep when the quadrotor has to reach the final state $x_f$, $\ell$ is a quadratic cost function, $\bar{u}$ is the open-loop control input (the jerk of the trajectory), $d$ is the UAV-tag distance, $x$ is the position, velocity, and acceleration of the UAV, $A$ and $B$ are, respectively, the state and input matrices for a triple integrator, and $a$ is the subvector of $x$ representing the acceleration.
Note that we can plan using this linear model for a nonlinear system because the nonlinear dynamics of the quadrotor are canceled by the ancillary controller, as explained next.

\subsection{Ancillary Controller: Boundary Layer Sliding Controller}\label{subs:controller}

While MPC has been used extensively in industry \cite{qin2003survey}, the control of systems with nonlinear dynamics requires expensive optimization. Sliding control \cite{slotine1991applied} has proven to be effective in quadrotors \cite{runcharoon2013sliding, yang2016attitude}. This control strategy guarantees bounds on the tracking error and has been combined with MPC \cite{rubagotti2010robust, lopez2018robust}. In our approach, we derive a nonlinear ancillary controller using sliding control that models the disturbances found near the landing platform.

In quadrotors, the attitude dynamics are much faster than the position dynamics, and thus control of both can be decoupled \cite{how2015linear}: the output of a controller is the setpoint for the other. The position and velocity controller is derived in this section to account for the turbulent wind present near the landing platform, and the attitude control is performed by a quaternion-based controller \cite{lopez2016low}.

The following derives the BLSC. Define the manifold $S(t)$ by
$    s=\dot{\tilde{x}}+\lambda \tilde{x} = 0,$
where $\tilde{x} = x - x_d$ and $\lambda > 0$. The objective of sliding control is to maintain $s=0$ at all times.
If the control action's frequency is high enough, zero tracking error is guaranteed \cite{slotine1991applied}.
This high control action is impractical in many applications because of actuator limits and the excitation of unmodelled dynamics.
An approach taken in \cite{runcharoon2013sliding, lopez2018robust} is BLSC, where the control discontinuity is smoothed in a thin boundary layer of thickness $\Phi$:
\begin{equation}
    \mathcal{B} := \{s :|s| \leq \Phi\}
\end{equation}

Consider a system whose dynamics can be expressed as
\begin{equation}\label{eq:dynamics}
    \ddot{x} =  f\left(x\right) + b \left( x \right) u + d
\end{equation}
where $d$ is the disturbance. Then, the BLSC strategy is
\begin{equation}\label{eq:u}
    u = \hat{b}^{-1}\left[\ddot{x}_{d}-\lambda \dot{\tilde{x}}-\hat{f}(\boldsymbol{x})-K \operatorname{sat}\left( \frac{s}{\Phi} \right)\right]
\end{equation}
where $\operatorname{sat}(\cdot)$ is the saturation function, $\hat{f}$ is the estimated acceleration caused by drag, and $K$ is determined by the uncertainty in the dynamics and the disturbance of the system. As noted before, we can plan in~\eqref{eq:opt_problem} using linear MPC because of the cancellation of $f$ in~\eqref{eq:dynamics} and~\eqref{eq:u}.

We generate turbulence using leaf blowers as shown in Figs.~\ref{fig:video-screenshot} and~\ref{fig:lbarray}.
The turbulent wind parameters are the mean $\bm{v_w}$ and standard deviation $\sigma$ of the speed. Define
$    \bm{V} = \bm{v} + \bm{v_w}$
where $\bm{v}$ is the quadrotor's speed. Then, $\bm{V}$ is the total wind speed relative to the UAV and $\hat{f}$ is
$    \hat{f} = \hat{c} \norm{\bm{V}} V,$
where $\hat{c}$ is the estimated drag coefficient of the quadrotor.

The quadrotor plans a trajectory to approach the UGV in the direction it is facing to match its speed, and thus is never outside the wind field generated by the leaf blowers during the landing maneuver. Therefore, it is reasonable to assume that this wind field is constant in the directions perpendicular to where the leaf blowers point to. The true acceleration caused by drag is
\begin{equation}
    f = c \norm{\bm{V} \pm 2\sigma \bm{u_w}} \left( V \pm 2\sigma u_w \right)
\end{equation}
where $\bm{u_w}$ is a unit vector in the direction of the wind.

The variation of $b$ is very small for a UAV with constant weight. Therefore, $\beta = \left(b_{max}/b_{min}\right)$ is approximately 1, where $b_{max}$ and $b_{min}$ are the maximum and minimum control gains respectively (or throttle gains in the context of quadrotors). Thus, $K$ is simplified as \cite{slotine1991applied}
$    K = \bar{F} + \eta,$
where $\eta > 0$ is a constant in the sliding condition
\begin{equation}
    \frac{1}{2} \frac{d}{d t} s^{2} \leq-\eta|s|.
\end{equation}
The larger the $\eta$, the faster the system will reach the sliding surface. Nevertheless, $K$ should only be as large as the disturbance magnitude requires to avoid a high-frequency control signal. $\bar{F}$ is
\begin{equation}
\begin{split}
    F =& | f - \hat{f} | \leq \bar{F}  \\
    \bar{F} =& | \left( \hat{c} + \tilde{c} \right) \norm{\bm{V} \pm 2\sigma \bm{u_w} } (V \pm 2\sigma u_w) - \hat{c} \norm{\bm{V}} V |
\end{split}
\end{equation}
where $\tilde{c} > 0$ is a bound on the absolute value of the drag coefficient error $|c - \hat{c}|$. By taking the sign that makes this coefficient larger, we have defined $K$. In our application, the quadrotor moves towards the generated wind and therefore this occurs when the $2\sigma$ is \textbf{increasing} the magnitude of $\bm{v_w}$.

\subsection{Landing Platform Estimation: Extended Kalman Filter}\label{subs:estimation}

To estimate the state of the moving platform, an extended Kalman filter (EKF) is used. This filtering algorithm minimizes the mean of the squared error and has demonstrated its efficacy in robot localization \cite{kiriy2002three, georgiev2004localization, teslic2011ekf}. The state vector of the platform is
$    \boldsymbol{x_p} = [p_x, p_y, v_p, \theta, \dot{\theta}]^\top,$
where $p_x,p_y$ are the 2D coordinates, $v_p$ is the magnitude of the velocity, $\theta$ is the orientation angle with respect to the $x$-axis, and $\dot{\theta}$ is the angular velocity.
Since real-world roads are mostly horizontal and in particular our experiments were carried on completely flat surfaces, the velocity in the $z$-direction is not estimated.
The moving platform is modeled as a unicycle with dynamics
\begin{equation}
    \boldsymbol{\dot{x}_p}(t)=f_p(\boldsymbol{x_p}(t))+\boldsymbol{w}(t),
\end{equation}
where $\boldsymbol{w}(t)$ is the process noise, assumed to be a white Gaussian noise. We consider a constant linear and angular velocity, and the UGV dynamics are
\begin{equation}
\begin{aligned}
\dot{p}_{x} &=v_{p} \cos (\theta), && \dot{v}_{p} = 0, \\
\dot{p_{y}} &=v_{p} \sin (\theta), &&
\ddot{\theta} = 0.
\end{aligned}
\end{equation}

The measurement vector is
    $\boldsymbol{z_p} = [p_x, p_y, \theta]^\top$, and
when a measurement is received, the EKF approach is used to perform an update of the estimated state.
The measurements are obtained in two ways.
First, when the quadrotor is far from the platform (that is, in the \textit{Search} state defined in Section~\ref{subs:fsm}), these measurements are obtained by adding a white Gaussian noise to the ground truth pose of the vehicle (obtained from the motion capture data) to simulate inaccurate GPS measurements that a ground vehicle could provide to the UAV for the rendezvous.
Note that receiving $\theta$ is not necessary to estimate the orientation of a moving platform because its motion could be used to infer that quantity.
Nevertheless, $\theta$ measurements are used in this work, which enables testing for static platform experiments. The update frequency is 2Hz, which is realistic for UAV applications \cite{salih2013suitability}. These first set of measurements are simply used to help the UAV locate the ground vehicle, but that information could be estimated without requiring a link between the two vehicles.

Second, and most importantly, when the quadrotor is near the ground vehicle and detects the onboard tag, we fuse both the simulated GPS and the visual detection measurements to estimate $\boldsymbol{x_p}$. The vision-based detection is explained in the next subsection, and it provides a far more accurate estimate of the position and orientation of the tag/platform. The 2D positions $p_x$ and $p_y$ and the heading angle $\theta$ are then used to update the platform's state, and the estimated velocity $v_p$ is incorporated into the vector $x_f$ in~\eqref{eq:opt_problem} to ensure the quadrotor matches the moving platform's speed at the landing point.

\subsection{Visual detection: AprilTag visual fiducial system}\label{subs:detection}
When the UAV is relatively close to the platform, visual estimation provides more accurate UAV-UGV poses than GPS. We used the AprilTag visual fiducial system \cite{olson2011tags} to obtain them, and a ROS wrapper \cite{malyuta2018guidance} based on AprilTag 2 \cite{wang2016apriltag} to interface with the core detection algorithm.
A tag bundle is a set of several coplanar tags used simultaneously by the visual fiducial system: the algorithm extracts the information of all of them to estimate a single ``bundle pose''. Thus, they are useful when accurate detection is required, which is the case in this paper.
Additionally, by using tags of different sizes, detection at various distances is ensured. We used a tag bundle comprised of a 14$\times$14cm tag on top of a 5$\times$5cm tag.
Despite the relatively small bundle size, it can be detected at a maximum distance of approximately 3.5m, and a minimum distance of about 5cm. The bundle can be seen in Fig.~\ref{fig:lbarray}.


\section{Introduction}

White dwarfs (WDs) are the final products of the evolution of average and low-mass main sequence stars. They are formed in the cores of red giant stars. Most of the stellar population will end up as a white dwarf star. Therefore, WDs are considered to be the most abundant stellar remnants. The average mass of a white dwarf is around 0.6 solar mass ($M_\odot$) and radius is roughly 10000 km. Correspondingly, their average density is approximately 10$^6$ g/cm$^3$ \citep{1971reas.book.....Z,1971tges.book.....Z,shapirobook}.

Unlike neutron stars (NSs), there are only a few equations of state (EoS) in the literature for WDs: the Chandrasekhar EoS \citep{chandrasekhar31}, the Salpeter EoS \citep{1961ApJ...134..669S}, the Relativistic Feynman-Metropolis-Teller EoS \citep{2011PhRvD..84h4007R} and other EoS which are more sophisticated and realistic \citep{2019MNRAS.490.5839B} . These equations with some modifications, including finite temperatures, magnetic field etc. are used to describe the cores of WDs and outer crusts of NSs \citep{haenselbook}.

The physical properties of WDs have been intensively studied both in Newtonian gravity (NG) and general relativity (GR) \citep{1971tges.book.....Z,shapirobook,1968Ap......4..227S, 1971Ap......7..274A,
1990RPPh...53..837K, 2011IJMPE..20..136B, 2015mgm..conf.2468B,2018GReGr..50...38C}. It has been shown that the effects of GR are crucial to analyze the stability of WDs close to the Chandrasekhar mass limit 1.44$M_\odot$ and can be neglected for low mass WDs \citep{2011PhRvD..84h4007R}.

Here the analyzes performed in Refs.~\citep{sheyse, b2016b, boshmg18} are extended to include the effects of rigid rotation in WDs employing the Chandrasekhar EoS at finite-temperatures with $\mu=2$ and $\mu=56/26$. The mass-radius ($M-R$), mass-central density ($M-\rho$), radius-central density ($R-\rho$), mass-angular velocity ($M-\Omega$), radius-angular velocity ($R-\Omega$) and angular velocity-central density ($\Omega-R$) relations for hot static and rotating WDs at Keplerian rate are constructed both in NG and GR for comparison.

The main goal of the paper is to jointly study the effects of rotation, finite temperatures and GR only, which usually are considered separately in the literature, without involving the Thomas-Fermi corrections, Coulomb interactions between electrons and ions, phase transition etc in the EoS. Though the corrections in the EoS are important in the accurate theoretical description of the $M-R$ relations of WDs, for the sake of simplicity the Chandrasekhar EoS is used throughout this work.

\section{Problem Setup and Results}\label{sec:fopres}
The physical characteristics of rigidly rotating WDs at different finite temperatures are studied here in the range of radius from 1000 km to 200 000 km. For clarity, the Chandrasekhar EoS is exploited since it is well-known and widely used to the description of WDs \citep{chandrasekhar31,2011PhRvD..84h4007R,b2016b,2019arXiv190810806C}.

The ratio of the atomic number $A$ to the number of protons $Z$ is usually denoted in the literature as $\mu=A/Z$ and all calculations in this paper are carried out by adopting $\mu=2$ for helium $^{4}_{2}$He, carbon $^{12}_{6}$C, oxygen $^{16}_{8}$O etc and $\mu=56/26\approx2.154$ for iron $^{56}_{26}$Fe for comparison.

Helium and iron are the two extreme bounds of chemical elements in WDs. All the rest elements such as carbon, oxygen, neon, magnesium etc. are between these two, correspondingly the mass-radius relations of other elements will be restricted by helium from above and by iron from below, see Ref.~\citep{2011PhRvD..84h4007R} for details.

Furthermore, WDs composed of helium were observed and their formation were analyzed in simulations with different astrophysical scenarious \citep{2001AN....322..405S, 2004ApJ...606L.147L, 2005MNRAS.362..891B, 2014A&A...571L...3I}. There are also a plethora of observational data for carbon and other WDs \citep{kepler2015, kepler2016}. The time and thermal evolution of carbon WDs has been theoretically investigated involving the nuclear burning and neutrino emission processes \citep{2019MNRAS.487..812B}. The observational support of the existence of iron-rich cores of WDs was given in Ref. \citep{1998ApJ...494..759P}. The structure and evolution of iron-core WDs has been studied in Ref. \citep{2000MNRAS.312..531P}. Some simulations demonstrate that the formation mechanism of WDs with iron-rich cores are different from the WDs composed of light elements \citep{2012ApJ...761L..23J}.

For the construction of $M-R$ relations the temperatures of the WD isothermal cores $T=T_c$ are considered here without taking into account the atmosphere of WDs. The  effective surface temperature $T_{eff}$, what is usually measured from observations, is roughly three order of magnitude less than $T_c$ according to the approximate Koester formula  $T_{eff}^4/g=2.05\times10^{-10} T_c^{2.56}$, where $g$ is the surface gravity \citep{koester2}.

To investigate static WDs in GR the Tolman-Oppenheimer-Volkoff (TOV) hydrostatic equilibrium equations are solved numerically. To construct rotating WDs the Hartle formalism is applied both in GR and NG \citep{Hartle,1968ApJ...153..807H, bosh2016} which generalizes the TOV equations by the inclusion of rotation. Note, that Ref~\citep{bosh2016} is a pedagogical paper explicitly showing the Hartle approach formulated in NG and comparing with other treatments in the literature in detail. All rotating WDs in NG are calculated at the Keplerian sequence: $\Omega_{Kep}=\sqrt{G M_{rot}/R^3_{eq}}$, where $\Omega_{Kep}=\Omega$ is the maximum angular velocity, $M_{rot}$ is the total rotating mass, and $R_{eq}$ is the corresponding equatorial radius of WDs. In GR the corrections to $\Omega_{Kep}$ owing to the angular momentum (the Lense-Thirring effect) and quadrupole moment (deformation of a star) are duly taken into account.

It should be stressed that though the Hartle formalism is designed and valid for slowly rotating stars in the approximation of $\sim \Omega^2$ (angular velocity of a star) one can extrapolate it to the Keplerian mass shedding limit for qualitative analyses. Indeed, Hartle and Thorne \citep{1968ApJ...153..807H} applied the approach for rapidly rotating stars (from supermassive stars to NSs) to study the effects of rotation on the structure. The validity and reliability of the formalism was tested in Refs.~\citep{berti2004,berti2005}. It was shown that the discrepancies between the Hartle formalism and exact computations appear close to the mass shedding limit. In addition, it has been demonstrated that the so-called Sedrakyan and Chubaryan formalism \citep{1968Ap......4..227S, 1971Ap......7..274A}, usually applied to the  study of rigidly rotating WDs and NSs, is identical with the Hartle formalism \citep{2016GrCo...22..305B}. Therefore, one can safely employ the Hartle formalism at the  mass shedding limit for rough estimates and qualitative analyses \citep{2000csnp.conf.....G,stergioulas}.

One can also study the effects of differential rotation, contribution of the atmosphere, diffent temperature profiles ect inside WDs in analogy with Refs. \citep{2019MNRAS.486.2982Y,2020MNRAS.492..978T}. However, those issues are out of the scope of the current work.

\begin{figure*}
\centering
{\includegraphics[width=3.2in]{MRforHe.eps}}
{\includegraphics[width=3.2in]{MRforFe.eps}}
\caption{Colour online. Mass-radius relations for WDs with $\mu=2$ (left panel) and for WDs with $\mu=2.154$ (right panel). In both figures the solid curves show static and the dashed curves show rotating WDs for different temperatures T=[0, 4$\times10^7$, $10^8$] K.}\label{fig:f1}
\end{figure*}

By solving the structure equations numerically,  $M-R$, $M-\rho$, $R-\rho$, $M-\Omega$, $R-\Omega$ and $\Omega-\rho$ relations for static and rotating hot WDs are obtained. The stability of hot rotating WDs has been analyzed in Refs.~\citep{bosh2018,bosh2018AR}. In Fig.~\ref{fig:f1} the $M-R$ relations are shown at different temperatures $T=[\ 0, 4\times10^7, 10^8]\ $K depicted with black, red and blue curves, respectively, for $\mu=2$ (left panel) and $\mu=2.154$ (right panel). All solid curves indicate static WDs and dashed curves correspond to rotating WDs. It is apparent that the finite-temperature effects are more prominent for larger in size (low mass) WDs and the effects of rotation are crucial for smaller in size (massive) WDs. Nevertheless both the finite temperature and rotation effects contribute to the radius and the mass of WDs. The comparison of these theoretical curves for $\mu=2$ with the estimated masses and radii of WDs from the Sloan Digital Sky Survey Data Release 4 (see Ref.~\citep{tremblay}) is given in Ref.~\citep{bosh2018AR}.

The behavior of the $M-R$ relations with $\mu=2$ and $\mu=56/26$ is similar, though the mass is larger for WDs with $\mu=2$. One can see that for a fixed mass within the range of the equatorial radius $R_{eq}=(5\times10^3-50\times10^3)$ km of Fig.~\ref{fig:f1} the higher the temperature the larger the radii. The same is true for a fixed radius: the higher the temperature the larger the mass.

\begin{figure
\centering
\includegraphics[width=1\columnwidth,clip]{MRmagnified.eps}
\caption{Colour online. Mass-radius relations of WDs close to the maximum mass.}\label{fig:f2}
\end{figure}

Fig.~\ref{fig:f2} is magnified Fig.~\ref{fig:f1} (left panel) in the range $R_{eq}=(1.1\times10^3-5\times10^3)$ km. Previously the range near the maximum mass was not analyzed thoroughly, albeit the general behavior of the $M-R$ relations was known \citep{bosh2018,bosh2018AR}. Within that range for static WDs, the higher temperatures the larger masses. This effect is natural, since due to the temperature, for a fixed central density, the pressure of partially degenerate and non-degenerate electrons increases and can sustain more mass. Instead, for rotating WDs, the higher temperatures the smaller masses. This effect is counter-intuitive, since one would expect that the effects observed in static WDs would automatically translate to rotating WDs. However, this is not the case, at least in this range of mass and radius. Here, the temperature affects the radius more than the mass. Hence, the Keplerian angular velocity is lower for hotter WDs, correspondingly the total rotating mass is less than for colder WDs.

\begin{figure*
\centering
{\includegraphics[width=3.2in]{MRmagnifiedGR.eps}}
{\includegraphics[width=3.2in]{MRmagnifiedFeGR.eps}}
\caption{Colour online. Mass-radius relations of WDs close to the maximum mass.}\label{fig:f222}
\end{figure*}

Fig.~\ref{fig:f222} shows $M-R$ relations close to the maximum mass both in NG and GR. As one can see that for the $\mu=2$ case (left panel) the role of GR is crucial. Instead GR is negligible for the $\mu=56/26$ case (right panel).

\begin{figure*
\centering
{\includegraphics[width=3.2in]{MassOmega.eps}}
{\includegraphics[width=3.2in]{ReqOmega.eps}}
\caption{Colour online. Mass-angular velocity relations (left panel) and equatorial radius-angular velocity relations (right panel) of WDs close to the maximum mass.}\label{fig:f3}
\end{figure*}

In Fig.~\ref{fig:f3} (left panel) the dependence of the total mass is constructed as a function of the angular velocity. It is clear that for a fixed mass cold WDs rotate faster. However the situation is opposite close to the maximum mass. Right panel of Fig.~\ref{fig:f3} shows the dependence of the equatorial radius on the angular velocity. Here for a fixed radius hot WDs rotate faster and close to the maximum mass the rotation rate of WDs is almost independent of the temperature.

\begin{figure*
\centering
{\includegraphics[width=3.2in]{Mrho.eps}}
{\includegraphics[width=3.2in]{Reqrho.eps}}
\caption{Colour online. Mass-central density relations (left panel) and equatorial radius-central density relations (right panel) of rotating cold and hot WDs.}\label{fig:f4}
\end{figure*}

To better understand the combined effects of rotation and finite temperatures on WDs $M-\rho$ and $R-\rho$ relations are constructed in Fig.~\ref{fig:f4}. Left panel of Fig.~\ref{fig:f4} shows that for a fixed central density rotating hot WDs in the considered density range will be less massive than colder ones. On the contrary for a fixed central density the radius of hot WDs will be always larger than the colder ones (see Fig.~\ref{fig:f4} (right panel)).

\begin{figure
\centering
\includegraphics[width=\columnwidth,clip]{Omegarho.eps}
\caption{Colour online. Angular velocity versus central density for rotating WDs.}\label{fig:f5}
\end{figure}

According to Fig.~\ref{fig:f5} for a fixed central density colder WDs will always rotate faster than hotter ones.

\begin{table*
\centering
 \caption{Static and rotating masses for hot white dwarfs for fixed central densities both in Newtonian gravity and in general relativity. The values of the mass in general relativity are given in parentheses.}
 \label{tab:example}
 \begin{tabular}{lccccccc}
  \hline
  & Configurations &$T=0$ K & $T=4\times10^7$ K & $T=10^8$ K \\
  \hline
$M/M_{\odot}$,  $\mu=2$ & Static & 1.440 (1.426) & 1.440 (1.426) & 1.440 (1.426) \\

$\rho_c=10^{10}$ g/cm$^3$& Rotating & 1.534 (1.519) & 1.532 (1.517) & 1.526 (1.511)\\
 \hline
 $M/M_{\odot}$,  $\mu=56/26$ & Static  & 1.180 (1.176) & 1.181 (1.177) & 1.181 (1.177) \\

$\rho_c=10^9$ g/cm$^3$ & Rotating  & 1.284 (1.279) & 1.278 (1.273) & 1.268 (1.263)\\
  \hline
 \end{tabular}
\end{table*}

In Table~\ref{tab:example} the static and rotating masses of WDs are presented corresponding to $\mu=2$ with the central density $\rho_c=10^{10}$ g/cm$^3$ and $\mu=56/26$ with $\rho_c=10^9$ g/cm$^3$. The values of the critical central density are chosen in accordance with the inverse $\beta$-decay and pycnonuclear reaction (for carbon WDs) instability densities, see Ref.~\citep{2013ApJ...762..117B} for details. As one can see that for static WDs the maximum mass almost does not change with the increasing temperature (see Fig.~\ref{fig:f4}), instead for rotating white dwarfs the maximum mass slightly decreases with the increasing temperature in both $\mu=2$ and $\mu=56/26$ cases.

\begin{table*
\centering
 \caption{Static and rotating masses for hot white dwarfs for a fixed radius both in Newtonian gravity and in general relativity. The values of the mass in general relativity are given in parentheses.}
 \label{tab:example2}
 \begin{tabular}{lccccccc}
  \hline
  & Configurations &$T=0$ K & $T=4\times10^7$ K & $T=10^8$ K \\
  \hline
$M/M_{\odot}$,  $\mu=2$ & Static & 1.329 (1.324) & 1.337 (1.332)& 1.347 (1.342)\\

$R_{eq}=3\times10^{3}$ km & Rotating & 1.495 (1.488) & 1.490 (1.483) & 1.484 (1.477)\\

  \hline
 \end{tabular}
\end{table*}

In Table~\ref{tab:example2} total static and rotating masses are given for a fixed equatorial radius (3000 km) at different temperatures in $\mu=2$ case only. Note that for a static case the equatorial radius reduces to the static radius. It is again evident that rotating hot WDs will possess less mass with respect to cold ones, in accordance with Fig.~\ref{fig:f2}. In $\mu=56/26$ case this effect is negligible due to the limitation in the central density caused by the inverse $\beta$-decay instability.

\section{Discussion and Conclusion}\label{sec:con}
The properties of static and rotating WDs have been investigated using the Chandrasekhar EoS at finite-temperatures. The structure equations have been solved numerically to construct $M-R$, $M-\rho$, $R-\rho$, $M-\Omega$, $R-\Omega$ and $\Omega-\rho$ relations for hot static and rotating isothermal cores of WDs. The atmosphere of WDs was not considered for simplicity.

It has been shown that temperature affects the masses of WDs at larger radii. At smaller radii the thermal effects are negligible. Rotation affects the masses of cold WDs in all density range. For hot WDs the effects of rotation are less noticeable at larger radii and more noticeable at smaller radii.

The $M-R$ relations were compared and contrasted for cold and hot, static and rotating WDs with $\mu=2$ and $\mu=56/26$. If the case with $\mu=2$ was well studied in the literature, the case with $\mu=56/26$ including the effects of both rotation and finite temperatures is investigated here for the first time. It turned out that for a fixed mass hot iron WDs are smaller in size, correspondingly denser with respect to the WDs composed of light elements.

It was found that close to the maximum mass the hotter rotating WDs possess less mass than the colder ones. This is related to the fact that unlike in the static case where the  mass of a WD is a function of the central density, radius and temperature, in the rotating case one more variable i.e. angular velocity is involved in a non-trivial way. Therefore, one observes that within the range $R_{eq}=(1\times10^3-5\times10^3)$ km  hot rotating WDs posses slightly less mass than the cold ones. To the best knowledge of the author this result is new.

In addition, to better comprehend the decrease in mass of hot rotating WDs it was useful to construct $M-\rho$, $R-\rho$, $M-\Omega$, $R-\Omega$ and $\Omega-\rho$  relations. Since in our case angular velocity is an additional degree of freedon its effects to the structure of WDs along with temperature are not straightforward.

In summary, one can see in Fig.~\ref{fig:f4} that for hot {\it static} WDs the maximum mass does not change (left panel) with increasing temperature, this is related to the Fermi temperature which is very high at this density range, instead the corresponding {\it static} radius (right panel)   changes: the hotter the WD the larger its radius. Hence, hotter WDs are fluffier than the colder ones. Consequently, larger in size hot WDs cannot rotate faster than smaller in size cold WDs with the same central density. If they cannot rotate faster, in order to fulfill the stability criteria i.e. the central density must not exceed a certain value; their rotating total mass must be lower than for the colder WDs.

Furthermore, for the sake of comparison all the computations have been performed both in NG and GR. As it was expected the effects appeared in NG are automatically translated into GR. The contribution of GR to the $M-R$ becomes relevant close to the maximum mass for $\mu=2$ cases. As for the $\mu=56/26$ case the effects of GR are negligible though are crucial for stability analysis and the correct estimation of the radius. Therefore the main focus was given to the $\mu=2$ case.

In view of the latest observational data on WDs \citep{kepler2019, 2019MNRAS.482.4570G} it would be interesting to explore WDs, taking into account the Coulomb interactions and Thomas-Fermi corrections at finite temperatures in the EoS as in Ref.~\citep{sheyse,Faussurier2017,2017JSMTE..11.3101F}, including the effects of rotation in the structure equations. That will be the issue of future studies.

\section*{Acknowledgement}
\medskip
\noindent
This work was supported by the Ministry of Education and Science of the Republic of Kazakhstan, Program IRN: BR05236494, Grants IRN: AP05135753  and IRN: AP08052311 .

\bibliographystyle{mnras}


\section{Introduction}
\label{sec:introduction}
This paper investigates methods for evaluating the interpretability of models of time series data arising  from people's interactions in immersive simulations
such as those  used for teaching in healthcare, disaster response and science education~\cite{alinier2014immersive,amir2013plan}.
In such simulations, people's interactions  engender a rich array of emergent outcomes and yield diverse opportunities for learning~\cite{smordal2012hybrid}. In the immersive simulation used in this study, Connected Worlds\footnote{Installed at the New York Hall of Science (NYSCI): \url{https://nysci.org/home/exhibits/connected-worlds/}} (CW), students interact with an ecological simulation to learn about the causal effects of their actions on environments over time~\cite{mallavarapu2019connect}.

Rich causal relationships, simultaneous participation from students and the changing dynamics of immersive simulations can make it difficult for people to determine how their interactions with the simulation caused the changes they observe in the simulated world.
Machine learning methods can be used to summarize the effects of participants' actions over various time periods.
For such methods to be effective, though, they must meet the challenge of identifying a model that is both ``true to the data'' and understandable to the target audience interested in uncovering the causal relationships.

This paper defines and solves the \emph{interpretability problem for immersive simulation settings}: determining that machine learning model, from a set of candidates, that people understand best~\cite{doshi2017roadmap,caruana2015intelligible}.
It compares the selection of a model according to a criterion that optimizes for maximum statistical information with one that optimizes for interpretability.
The ability to identify the model that is best (or among the top choices) for interpretability is essential to a system's capability to explain its conclusions~\cite{rosenfeld2019explainability}.

Our approach to addressing the interpretability problem comprises the following:
(1) select a set of machine learning models for segmenting time series data; in the domain we investigated, the segmentation is of students' interactions with CW into coherent periods of time;
(2) design tests for computing the interpretability score of a model for a given input;
(3) empirically evaluate the models with respect to their interpretability score in a user study.

To infer the boundaries of stable periods in the data of CW dynamics, we use a family of hidden Markov models (HMMs).
These HMMs are augmented with an additional ``sticky'' hyperparameter which biases the transition dynamics of the latent state-space~\cite{fox2008hdp}.
The input to each HMM is a multidimensional time series representing the response of the CW system to actions performed by students in the simulation.
The output of the HMM is a segmentation of the time series into a set of periods, which are contiguous lengths of time during which the system dynamics form a stable linear process.

We implemented two tests of interpretability for CW models: the Forward Simulation and Binary Forced Choice \cite{doshi2017roadmap}.
These tests each determine the extent to which the learnt representations are interpretable to people, albeit in different ways.
They both use a visualization of the inferred periods that shows experimental subjects snapshots of the CW system's state from the selected periods that the HMM inferred.

The results showed that the interpretability of the different models varied according to the value(s) of HMM parameters.
In particular, the HMM that optimized statistical information criteria did not optimize interpretability quality.
In addition, a fully Bayesian approach, which does not require hyperparameter tuning, offered a good balance between interpretability and performance on the theoretical statistical tests.
We argue that the Bayesian approach could be suitable for situations in which it is not possible to engage people in determining interpretability or doing so would be unethical or impractical.

This paper makes three contributions.
First, it provides an end-to-end paradigm for the design and evaluation of the interpretability of models for unsupervised learning in time series domains.
Second, it defines new interpretability tests for unsupervised time-series settings, and applies them to real-world data.
Third, in identifying the Bayesian solution, it provides an attractive alternative to model selection when human subject experimentation is not possible.
Finally, we note that the results of this investigation have been deployed in a classroom study for the purpose of assisting teachers in explaining systems thinking to students who participated in the CW simulation study.

\section{Related Work}

\citename{doshi2017roadmap} suggested three tests to evaluate how interpretable a model's representations are to people.
Forward Simulation: requires a human evaluator to predict the output of a model for a given input.
Binary Forced Choice: requires an evaluator to choose one of two plausible model explanations for a data instance.
Counterfactual Simulation: requires an evaluator to identify what must be changed in an explanation to correct it for a given data instance.

In follow-up work~\citename{lage2018human} propose a model selection process that considers both a model's accuracy and its degree of interpretability, according to one of the above tests.
They provide a framework for iteratively optimizing the interpretability of a model with a human-in-the-loop optimization procedure.
Their work applied this framework to tests in the lab in which human judgment was used to optimize supervised learning models.
Other works that studied interpretability tests for supervised learning settings include~\citename{wu2018beyond,ribeiro2016should,choi2016retain,lipton2016mythos}.
We extend this literature on interpretability by adapting the model selection process to an unsupervised learning setting, that of segmenting a multi-dimensional time series into periods.
Moreover, we implement examples of the Forward Simulation and Binary Forced Choice tests suggested by~\citename{doshi2017roadmap} and apply them to a high dimensional time series setting.

Our work was inspired by \citename{chang2009reading} who were the first to show that optimizing machine learning models in unsupervised settings using predictive log-likelihood may not induce models that are interpretable to people.
They focused on the use of topic models for finding meaningful structure in documents and they compared the models that are selected to optimize \textit{perplexity} (analogous to held-out log-likelihood) to the models that were selected by the human interpretability tests that they designed.
\citename{chang2009reading} operationalized two Forward Simulation tests for evaluating the interpretability of a topic model: word intrusion, in which the evaluator is required to identify which of several words does not belong together in one topic represented by the other words; and  topic intrusion, in which the evaluator is required to
identify which of several topics is not associated with a given document.
We extend this work to a multi-dimensional time series domain and we introduce a Binary Forced Choice test to complement the ``intrusion'' Forward Simulation test.

\section{The Connected Worlds Domain}
\label{sec:cw_desc}
Connected Worlds (CW), a multi-person ecology simulation (installed at NYSCI), aims to teach students about complex systems and systems thinking. Its immersive environment comprises four biomes (Desert, Grasslands, Jungle \& Wetlands) connected by a central water flow fed by a waterfall.
Students plant trees which flourish or die, animals arrive or depart, and rain clouds form and rain feeds the waterfall.

Students control the direction of water flows in the simulation by moving foam logs to direct water among the biomes.
Water enters the simulation through rainfall events, which are not under student control.
Figure~\ref{fig:connected_worlds_graphic} gives a snapshot of the system state, we refer to this snapshot as the session-view.
This session-view is a system-generated representation of the water flows and it directly reflects the logged water flows and levels in the simulation.
The output of a CW session is a time series recording the levels of water in the different biomes for 8 minutes at a 1Hz frequency.
The ability to model the effects of student actions on the environment were limited by two factors: The time series was the only source of information about students' interactions, and it was not possible to access the CW simulation except at NYSCI.

\begin{figure}[t]
\centering
\includegraphics[width=0.39\textwidth]{schematic_new.png}
\caption{CW session-view. Biomes are labelled on the perimeter and logs appear as thick red lines. Water (blue stream in the middle of the image) enters via the waterfall and in this image it mainly flows toward the Grasslands and the Desert.}
\label{fig:connected_worlds_graphic}
\end{figure}

The CW simulation is complex on several dimensions as a large number of students simultaneously execute actions that change the state of the simulated environment.
Each participant has a different view of what transpired, depending on the actions s/he took and the environment changes that resulted.
Students' activities are recorded as a movie (see Figure~\ref{fig:connected_worlds_graphic}) that can be shown to students and teachers.
This movie can inform discussions about the causal effects of the students' actions on simulation outcomes, but it obscures temporal dependencies in their interactions.
This limitation motivated the use of ML algorithms to better support students' understanding of the effects of their actions on the simulation's progression.

\section{Interpretability Tests for CW}
\label{sec:interpretability_tests}

\begin{figure}
\centering
\includegraphics[width=0.42\textwidth]{period_selection_representation}
\caption{
The time series is represented as a horizontal line from minute 0 to 8; red vertical lines denote sampled time points in the time series; each model is shown as a grey rectangle; models segment time series into periods delimited by white vertical lines. The forward or backward neighbour of the candidate period is selected as an intruder.}
\label{fig:selection_representation}
\end{figure}

\begin{figure*}[t]
\centering
\includegraphics[width=.85\textwidth]{test1_new.png}
\caption{Screenshot of the Forward Simulation test interface. Here $4$ of the images show water flowing towards the Desert. An intruder image, the highlighted one, comes from a different period and shows water flowing to both the Desert and the Grasslands.}
\label{fig:test1_screenshot}
\end{figure*}

Let $D$ be a time series that records the levels of water in the different biomes.
Let $M$ be a model that takes as input a time series $D$ and outputs a segmentation of $D$ into periods.
Each period aims to provide a coherent description of the water flow for a length of a time.

Importantly, a single period is insufficient for modeling the effects of students' interactions with CW, because students' sustained actions have complex effects on the system dynamics over time.
For example, when students choose to direct water to the Desert and Plains and plant trees in the Desert, the system dynamics are entirely different from the case when water is directed towards the Jungle and the Desert, and the Plains are left to dry.
We must therefore allow for multiple periods.
Each period describes a length of time where water flowed to a sufficiently stable target.
From the above example, one period can describe water that mainly flows to the Plains and to the Desert; students then move logs to re-route water flow to the Jungle, thus starting a new period.

We use an interpretability score $IS$ to measure the interpretability of a model $M$ applied to $D$.
The interpretability score is computed via an average across test instances, $T(M,D,i)$, which each take as input a model $M$, a time series $D$ and a selected point in time $i$ from the time series.
Each test instance returns True if an evaluator successfully completes a required objective.

We adapted the Forward Simulation and Binary Forced Choice tests~\cite{doshi2017roadmap} to the CW domain using the notion of \textit{candidate} and \textit{intrusion} periods.
We say that period $p$ is \textit{active} for model $M$ at time $i$ if $M$ infers the period $p$ to describe a contiguous length of time in the time series, and $p$ includes the time $i$.
Figure~\ref{fig:selection_representation} shows how the tests select candidate and intrusion periods.
First, a time point (red vertical line) is used to select a candidate period where the candidate period is the active period from model $M$ at $i$ (the active period for a model intersects with the red line).
Then, the intrusion period is selected as a direct neighbor to a candidate.
Each test is operationalized via a \emph{visualization} which presents any period as a set of images extracted from the session-view.

Figure~\ref{fig:test1_screenshot} shows an example of the Forward Simulation test on a real data instance.
As shown by this figure, the test sampled $4$ session-view images from the candidate period of model $M$ at time $i$, and a single session-view image sampled from the intrusion period.
The images were presented in a random order.
In Figure~\ref{fig:test1_screenshot}, the image that is outlined in green is the intrusion image that corresponds to the intrusion period.
A test evaluator was required to identify which image was the intrusion image.

Figure~\ref{fig:test2_screenshot} presents an example of the Binary Forced Choice test.
The test displays an unknown session-view image from a candidate period (center of screen) and additional images from two competing periods that contain this image (``Period 1'' or ``Period 2'').
Each of the two competing periods is visualized as four images sampled from the candidate or the intruder period.
The unknown image is sampled in time close to the boundary of when the candidate period transitions into the intruder period.
In Figure~\ref{fig:test2_screenshot}, Period 1, highlighted in green, is the period that correctly explains the unknown image (i.e., the images in ``Period 1'' and the ``unknown image'' are all sampled from the candidate period).
A test evaluator is required to choose between the two possible periods.

Hypothetically, the intruder period can be chosen arbitrarily, as in~\citename{chang2009reading}.
However, intrusion periods that are further away in time from the candidate period would be easier to detect due to the non-stationary evolution of the system.
We made a design decision to chose the period that is immediately adjacent to the candidate period, either forward or backward in time.
This makes it harder to distinguish between candidate and intrusion period, but provides a rigorous test for the specific choice of boundary between the two periods.

Given data set $D$ and model $M$, the interpretability score $IS$ of a model is equal to the average success of the test instances for model $M$ over multiple points $\{i\}$ in a time series $D$.
The set of time points $\{i\}$ were uniformly sampled from the time series with the additional constraint that each minute of interaction had at least one sample.
For every model we test, we hold constant the selected times $\{i\}$ in the time series (as shown in Figure~\ref{fig:selection_representation}).
In this way we control for different areas in the time series being more or less difficult to segment into coherent periods.

\begin{figure*}[t]
\centering
\includegraphics[width=0.62\textwidth]{test2_new.png}
\caption{Screenshot of the Binary Forced Choice user interface. An unknown center image needs to be associated with either ``Period 1'' or ``Period 2''. In this case, streams of water flowing to both the Grasslands and to the Jungle capture the dynamics in Period 2. Period 1 has a small amount of water reaching the Desert which is consistent with the unknown image.}
\label{fig:test2_screenshot}
\end{figure*}

\section{Modeling Students' Activities in CW}
\label{sec:model_for_segmenting_time}

In this section we describe the design of general models for segmenting students' activities into periods of time and thereafter present the specific classes of model that are used in our interpretability tests.

\subsection{Segmenting Time Series Data into Periods}
\citename{hoernle2018modeling} used a HMM to model the system responses to students' activities in CW in which the latent states of the HMM corresponded to periods.
Transitions between different states equates to the system changing between different periods, while self transitions mean the system persists within the same period.
The authors did not address the question of how to choose the number of states. To this end, we augment the HMM with a hierarchical Dirichlet process which places this non-parametric prior over the state space, following the approach detailed by \citename{teh2005sharing} and \citename{fox2008hdp}.

The ``Sticky-HMM" approach introduced by \citename{fox2008hdp} includes a hyperparameter, $\kappa$, that biases the model to persist in a state, given that it has already adopted that state.
Applied to CW, the greater the value for $\kappa$, the more the model will try to persist in any given state.
The increase in the length of periods corresponds to a decrease in the number of latent states.
The opposite is true for lower values of $\kappa$ where there is a lower bias to persist within a given state and consequently there are more periods that are inferred.
For a detailed description of the model, including the  Gibbs sampling inference scheme that is used to infer the model parameters, refer to \citename{fox2008hdp} and \citename{fox2009bayesian}.

\subsection{Model Classes}
We introduce three classes of model that segment time into periods that can be used to explain the water flows:

\begin{enumerate}
    \item  {${MK_X}$: sticky HMM with fixed $\kappa$.} We use the basic structure of the sticky HMM described by \citename{fox2008hdp} with set values for $\kappa$ to produce $10$ unique models, spanning a wide range of possible settings\footnote{$\kappa \in \{1, 5, 10, 50, 100, 150, 200, 300, 500, 700\}$.}.
    \item  {${FB}$: fully Bayesian sticky HMM with Gamma prior on $\kappa$.} This  approach places a weakly informative, conjugate Gamma prior on the hyperparameter that expresses high uncertainty over the $\kappa$ values\footnote{The \textit{(shape,rate)} parameters were chosen to be $(1,\frac{1}{4})$; empirical results were invariant to a range of these values.}.
    \item  {${Rand}$: Random baseline.} The random baseline generates periods of random length drawn from a Poisson distribution with mean set to be the mean of all other periods induced by the parametric models. The random periods are defined to include the selected time points ($\{i\}$ from Section~\ref{sec:interpretability_tests}).
\end{enumerate}

\noindent We refer to ${FB}$ as the fully Bayesian model to indicate the fact that the none of the parameters of interest are specified and consequently posterior inference is over all of the parameters in the model (including $\kappa$). This is in contrast to the ${MK_X}$ models where we explicitly set the value for the sticky parameter $\kappa$.

For models in class $1$ and $2$, we use the Gibbs sampler, described by \citename{fox2008hdp}, to perform inference over the parameters in the model, this includes inference over the state sequence and thus the period segmentation of the model. The observation distribution was chosen to be a mixture of two multivariate Gaussians with conjugate Normal-inverse-Wishart priors. This mixture model addresses the noise in the CW water flow, such as ``splashes'', which prior work has identified as a challenge in this domain~\cite{hoernle2018modeling}.

\section{Model Selection for Interpretability}
\label{sec:model_selection}
The goal of model selection is to optimize a metric such that a specific parameter setting can be chosen as the best model for use during inference.
We compare how the models from section~\ref{sec:model_for_segmenting_time} perform on both statistical tests and on the human interpretability tests outlined in section~\ref{sec:interpretability_tests}.

\subsection{Selection using Statistical Information}
\label{sec:baseline}
When  human interpretability testing is infeasible,  one could choose to optimize some proxy to interpretability~\cite{doshi2017roadmap,lage2018human}.
For example, \citename{chang2009reading} compared the proxy of held-out log-likelihood to the human interpretability score that was a result from two tests that were run on Amazon Mechanical Turk (Mturk).

Ideally, the model parameters would be optimized on held-out data using predictive log-likelihood as the objective~\cite{chang2009reading}.
However, the difficulty of collecting controlled sessions of student interaction in CW meant we had few data instances available (see limitation discussion in the next section).
To address this challenge we use statistical information criteria as a theoretical approximation to the predictive accuracy of a model~\cite{gelman2013bayesian}.

Figure~\ref{fig:dic_waic_scores} shows the two information criteria (the Deviance Information Criteria, DIC, and the Watanabe-Akaike Information Criteria, WAIC~\cite{gelman2013bayesian}) plotted as a function of the model (the random model has no notion of information criteria and so was not compared here). The data set comprised of both of the log files of students' interactions (8 minutes each). The optimal model for both DIC and WAIC is the $MK_{5}$ model but we note that $MK_{1}$, $MK_{5}$ and $MK_{10}$ all perform close to this optimal setting. Notice that the fully Bayesian model (FB) is not optimal but it is in the top $5$ models for both criteria.
\begin{figure}[t]
\centering
\includegraphics[width=0.42\textwidth]{dic_waic_scores.png}
\caption{DIC and WAIC as a function of the model (lower  is better). The $MK_5$ model is optimal, the $FB$ approach is in 5th place.}
\label{fig:dic_waic_scores}
\end{figure}

\subsection{Selection using Interpretability Test}
\label{sec:interpretability_results}
This section describes the choice of model according to interpretability quality, as measured by the interpretability tests.
The set of models used in this study includes the 12 CW models described in Section~\ref{sec:model_for_segmenting_time}.
IRB was obtained for the study.

We recruited participants from two cohorts: undergraduate engineering students in a large public university and Mturk workers (with a total of $240$ people who participated in the experimentation).
For a given time series $D$ in CW, we randomly sampled a set of 12 time points, which remained constant across all model conditions.
Each time point was used to generate a candidate and two intrusion periods (both forward and backward in time, see Figure~\ref{fig:selection_representation}), making for $2 \times 12 \times 12=288$ tests per time series.
We divided participants into two cohorts, one for Forward Simulation, and one for Binary Forced Choice tests.
Both cohorts varied the models used to generate their respective tests.
Each participant performed $20$ tests, with no more than $2$ tests generated from any given model, to ensure a representative range of models.
After making their choice, participants received brief visual feedback on whether or not their selection was in agreement with the model's choice.

All participants received a detailed tutorial about CW and the study, as well as a pre-study comprehension quiz\footnote{Tutorial pdf slides are available at \url{https://www.dropbox.com/s/pu2nxk2k0g81ql6/forijcai.pdf}}. Mturk workers were paid a base rate of $\$0.25$ for participating and a bonus structure of $\$0.1$ for each correct response.

\begin{figure}[t]
\includegraphics[width=.45\textwidth]{combined_res_data_all.png}
\caption{Effect of each model on the log-odds of a test evaluator selecting the correct response (controlling for the test evaluator, the experiment trial, log file and ordering effects).}
\label{fig:essil_test_results}
\end{figure}

We first describe results in terms of accuracy (the percent of correctly labelled test instances).
The top performing model was $MK_{200}$ with an accuracy of $83\%$ on the Forward Simulation  test and $MK_{100}$  with an accuracy of $82\%$ on the Binary Forced Choice test.
The random baseline model performed consistently poorly with an average accuracy of $53\%$ on both tests.
The fully Bayesian model achieved an accuracy of $72\%$ and $70\%$ respectively on the two tests.

To control for ordering effects, chosen time periods, data instance used, and effects of individual participants, we applied an L2 regularized logistic regression for predicting  the user specific success on the experiment trial, shown in  Figure~\ref{fig:essil_test_results}. The y-axis presents the improvement in log-odds that a model has on the expected response accuracy (higher is better).
As shown by this figure, the Forward Simulation shows a high variance with no clear maximum. In contrast, the Binary Forced Choice test has a clear maximum in the region of $MK_{100}$ and $MK_{150}$.

From Figures \ref{fig:dic_waic_scores} and \ref{fig:essil_test_results} we can infer the following four conclusions.
First, all of the models ($MK_{1},\ldots,MK_{700}, FB$) outperform the random baseline: participants are more likely to select the correct response from any of these models. This result suggests that periods of stable dynamics exist in the data and that it is possible to construct models, which describe these dynamics, that are interpretable to people.

Second, the Binary Forced Choice test is a preferable measure for interpretability to the Forward Simulation test. Figure~\ref{fig:essil_test_results} shows that the Binary Forced Choice test exhibits a clear peak (around $MK_{100}$ and $MK_{150}$) where interpretability of the model is maximized.
These models also maximized the  raw accuracy on the Binary Forced Choice test.

On the other hand, the Forward Simulation test has a greater variance across models and across data instances.
Two possible causes for this higher variance are: (1) there is more room for error in the Forward Simulation test (5 choices vs. 2 choices in Binary Forced Choice); (2) sampling a single image to represent a period (as in Forward Simulation) presents less information to the user than sampling 4 images (as in Binary Forced Choice).

Third, the best $\kappa$ settings vary for different tests and information criteria.
Model interpretability grows steadily as the value of $\kappa$ increases, with $MK_{100}$ and $MK_{150}$ being the optimal models, and then proceeds to decrease steadily.
These models are not consistent with the model $MK_5$ that optimized the information criteria.
Note that higher $\kappa$ values are ``sticky'' - they  bias the model towards longer periods, which condense too many activities to make sense to people.
On the other end of the spectrum, lower $\kappa$ values allow for more (shorter) periods that may  capture  noise in the system.   The $\kappa$ value for  models $MK_{100}$ and $MK_{150}$ represent a ``sweet spot" in between these two extremes.

Finally, the fully Bayesian model $(FB)$ performs consistently well on both information criteria and interpretability tests. It is interesting to note that while this model does not find the optimal setting (from neither the statistical information criteria nor from the human interpretability task) it does perform well across all tests, tasks and instances, and is fully automated (no human evaluation is required in order to choose an optimal parameter setting).

We conclude this section with mentioning the limitation that the user study was based on a small number $(n=2)$  of instances. This was due to the difficulty in obtaining controlled sessions of student behavior in CW. Despite this issue, the differences between the models in Figure~\ref{fig:essil_test_results} are statistically significant, having being evaluated across 12 different time points for each instance and with hundreds of evaluators.

\section{Conclusion \& Future Work}
With the growing prevalence of immersive simulations the need arises for AI systems which help people gain insight into the ways participants' activities affect the simulation outcomes.
We have studied an environmental simulation intended to teach students  about the causal effects of their actions.
Our results show that algorithms can segment time series log data into periods that are meaningful for people.
Selecting hyperparameters in these models is a challenge, especially when trying to optimize the representations they produce for their interpretability.
We have described ways to select these hyperparameters using two tests that are grounded in the literature.
We showed that the fully Bayesian method is a promising technique for implementing a model when people cannot directly assess and evaluate the models.
Our results are important for any unsupervised machine learning task for which interpretability is an important criterion, because in such cases the model selection problem will be encountered.
The work forms part of a broader project where the goal is to generate relevant summaries of the CW dynamics such that teachers can effectively engage their students in discussions about their own experiences with the simulation.

In future work we plan to explore alternative ways to measure interpretability quality in time series domains, including the design of a counterfactual simulation test~\cite{doshi2017roadmap},  and the application of our approach to additional domains.

\section*{Acknowledgments}

Nicholas Hoernle is supported by a commonwealth scholarship. This paper is based on work supported by the National Science Foundation under grant numbers \#1623091, \#1623124 and \#1623094.
 Thank you to  Michael Rovatsos and  Alex Lascarides  for their advice on this project, to the reviewers for their helpful comments and to NYSCI for assisting with the data collection.

\bibliographystyle{named}