\section{1. Introduction}
\addcontentsline{toc}{section}{Introduction}

Crystallization from droplets is a common process in the separation and purification of compounds in solutions. The most notable example is oiling-out crystallization, a crystallization process with a secondary liquid phase appearing (liquid-liquid phase separation, LLPS) during the cooling or adding an anti-solvent\cite{kiesow_2010_solubility, lafferrre_2004_study, derdour_2010_a}. The oiling-out phenomenon was observed and reported in the crystallization of active pharmaceutical ingredients (ibuprofen\cite{codan_2012_phase}, erythromycin ethylsuccinate\cite{li_2016_process}, idebenone\cite{lu_2012_study}, ethal- propyl- and butyl- paraben\cite{yang_2010_solubility,yang_2013_nucleation}) and food components such as vanillin\cite{zhao_2012_solution, du_2016_the}, and lauric acid\cite{maeda_1997_separation}. After oiling-out, the solution becomes cloudy due to the formation of the droplets (dispersed phase) in the continuous phase but eventually separates into two layers of subphases. If the droplet phase has a high solute concentration, the oiling-out processes lead to agglomeration, lowering the quality and purity of the final crystal products\cite{deneau_2005_an,gielen_2017_agglomeration}. It is important to understand and control the phase-phase separation \cite{codan_2010_phase,sun_2018_understanding, yang_2014_sandwich} and design the oiling-out process. To control the crystal size distribution and quality of the crystals \cite{sun_2018_oilingout, yang_2014_influence, li_2016_oiling}, the crystallization processes can be designed to avoid oiling-out by lowering the initial concentration, adjusting temperature and pressure, seeding, or employing ultrasounds\cite{ogrady_2007_the, tari_2019_comparative, beck_2010_controlled, kluge_2012_emulsion}. On the other hand, some research was reported to use the oiling-out to fabricate spherical agglomeration for the convenience of formulation in downstream processing\cite{sun_2018_oilingout}.

Our latest work demonstrated that oiling-out crystallization may provide a simple approach to prepare seed crystals from a stock solution at a very low concentration of the solute\cite{zhang_2020_oilingout}. The oiling-out crystallization of beta-alanine as a model compound occurred in a process called solvent exchange, where a good solvent (high solute solubility) was continuously replaced by a poor solvent (low solute solubility).  During the solvent exchange process, the droplets of alanine solution nucleated on the substrate and the water (the good solvent) in the droplets were continuously replaced by isopropanol (the poor solvent). With more good solvent exchanged, the droplets entered into solid-liquid-liquid phase region where the crystallization occurred\cite{zhang_2020_oilingout}. We have demonstrated that the crystal morphology and the crystallization processes on the homogeneous surfaces were largely influenced by the flow rate and concentration ratio of good and poor solvents.

It is known from the formation of ternary liquid nano-droplets that the droplet size is correlated with the initial solution concentration, flow rate, or the channel height during the solvent exchange \cite{zhang_2015_formation,yu_2015_gravitational,dyett_2018_coalescence, qian_2019_surface}. On a homogeneous surface, the nucleation and coalescence of the droplets are random, leading to the heterogeneous final size of droplets by the end of solvent exchange. However, highly ordered droplet arrays can be formed on chemically micro-patterned surfaces, as the location and spatial arrangement of the droplets follow the wettable microdomains on the surface \cite{bao_2015_highly,bao_2016_controlling, devi_2017_sessile, encarnacinescobar_2019_morphology}. For the crystallization of lipid nano-droplets in arrays, the small-angle X-ray scattering (SAXS) measurement during the cooling crystallization showed faster crystallization in comparison to bulk lipid crystallization due to the droplet morphology\cite{dyett_2018_crystallization}. We expect that the size and location of oiling-out droplets are influenced by micro-patterns on the surface, which may lead to different shapes and sizes of crystals from the solvent exchange crystallization of these droplets.

In this work, we have studied the effects of micro-patterns and micro-structures on the surface on the formation of oiling-out droplets and crystals during solvent exchange. Furthermore, we will show that solvent exchange method can induce crystallization of droplets even without the oiling-out process. Our work demonstrates that the oiling-out crystallization from droplets with solvent exchange is a potential method to screen and explore the new morphology of the crystalline product. The results in this work contribute to further understanding the crystallization from droplets in order to integrate with the application of pharmaceutical production and other future implementations.

\section{2. Experimental section}
\addcontentsline{toc}{section}{Experimental section}

\subsection{2.1 Substrate preparation}
\addcontentsline{toc}{subsection}{Substrate preparation}

The bare silicon substrates were purchased from University Wafer (South Boston, MA, US). The bare silicon wafer was used as bases for forming octadecyl trichlorosilane (OTS) coated substrate and 3-aminopropyltriethoxysilane (APTES) coated substrate. To prepare the wafer, the bare silicon wafer was cleaned by sonication in deionized water (from Milli-Q Direct) and then in ethanol for 10-20 minutes each before use. The bare silicon substrate was cleansed by piranha solution for 10-20 minutes at 75 \textdegree C by a hot plate, then sonicated by deionized water for 10-20 minutes.

Bare silicon substrates were coated by OTS by following the protocol reported in the literature\cite{zhang_2008_nanobubbles}. A similar protocol was followed for APTES coating on the substrate as well. The OTS-coated wafers were used as the bases for fabricating the patterned substrate via standard photolithography followed by plasma etching as indicated in the works of Bao et al.\cite{bao_2015_highly,lei_2018_formation}. 20 \textmu m hydrophilic domain surrounded by hydrophobic background was used as hydrophilic patterned substrates, and 5 \textmu m hydrophobic domain surrounded by hydrophilic background was used as hydrophobic patterned substrates. 1, 6-hexanediol diacrylate (HDODA) micro-lenses on the substrate were prepared by following the protocol in the literature.\cite{bao_2015_highly,lei_2018_formation}

 The contact angle of the water droplet in the air on APTES-coated substrate was approximately 50-60 degrees. The contact angle from hydrophilic and hydrophobic components of the patterned substrate was approximately 20-30 and 100-110 degrees respectively.

\subsection{2.2 Solution preparation}
\addcontentsline{toc}{subsection}{Solution preparation}

Beta-alanine (ACROS organics, 99 \%), ethanol (Fisher Scientific, HPLC Grade, including 90 \% ethanol, 5 \% methanol, 5 \% isopropanol), isopropanol (IPA, Fisher Scientific, 99.9 \%) were mainly used for solution preparation. The solution was stored in a sealed container overnight before use. Solution A was 3 \% beta-alanine, 45 \% water, and 52 \% isopropanol. The solutions were prepared by mixing the beta-alanine with isopropanol and water mixture followed by sonication until the solution became homogeneous. Solution B only consisted of isopropanol. All the concentration ratios are by weight percent unless indicated otherwise.

2,5-Bis(2-(4-pyridyl)-vinylene) hydroquinone dineopentyl ether (Np-P4VB) was prepared by following the protocol reported in literature\cite{lane_2018_widegamut}. The solution of Np-P4VB was prepared similarly to the preparation of the alanine solution. The solvent ratio for solution A consisted of 15:4:0.35 (ethanol:water:mesitylene) with approximately 0.01 g of Np-P4VB per 50 mL of the solvents. A small amount of mesitylene (Sigma-Aldrich, 98 \%) was added to minimize the use of the solute as Np-P4VB was artificially created in a small quantity. The prepared solution was then stored away from the light as the solute was sensitive to light.

\subsection{2.3 Procedure of oiling-out crystallization by solvent exchange}
\addcontentsline{toc}{subsection}{Procedure of oiling-out crystallization by solvent exchange}

\begin{figure} [htp]
	\includegraphics[trim={1.5cm 4.9cm 7.7cm 2.5cm}, clip, width=0.9\columnwidth]{Figures/Updated/1.jpg}
	\caption[Experimental apparatus and ternary phase diagram]{a) Sketch of the side view of fluid chamber. The solvent exchange process is progressed as alanine solution is replaced by isopropanol. b) Ternary phase diagram of beta-alanine, isopropanol, and water at 25 \textdegree C taken from Sun et al.\cite{sun_2018_oilingout} Regions on the diagram are solid-liquid region (S-L), solid-liquid-liquid region (S-L-L), liquid-liquid region (L-L), and homogeneous liquid region (L).}
	\label{f1}
\end{figure}

The process of solvent exchange started with solution A being filled in the flow cell as shown in Figure \ref{f1} (a). Solution A was then displaced by solution B. The flow of solution B continued until all droplets crystallized and ran for 1-2 minutes extra for possible incompletion of crystallization. The syringes and flow cell were subjected to washing and sonicating by ethanol and water to avoid possible contamination or crystal leftover after each experiment. After sonication, the flow cell and syringes were dried by compressed air before assembling the flow cell for the next experiment. For the adjustment of the channel height (H), layers of double-sided tape were used to create a higher platform for the substrate to reach from 100 \textmu m to 300 \textmu m.

The ternary phase diagram of beta-alanine, water, and isopropanol is referenced in Figure \ref{f1} (b). There are five regions within the ternary phase diagram: solid-liquid region (solute lean), solid-liquid-liquid region, solid-liquid region (solute rich), liquid-liquid region, and homogeneous liquid region. The data for the phase diagram was acquired from beta-alanine experimental works from Sun et al.\cite{sun_2018_oilingout}

\subsection{2.4 Solute concentration, flow rate and channel height}
\addcontentsline{toc}{subsection}{Solute concentration, flow rate and channel height}

  The conditions for the main experiment can be viewed in Table \ref{tp2}.

\begin{table}[h!]
\centering

\begin{tabular}{ |p{1cm}|p{3cm}|p{3cm}|p{3cm}|p{3cm}|  }
 \hline
 No. & Solution A & Flow rate & channel height & Peclet number\\
 \hline
 1 & 1.5 \% alanine & 6 mL/hr & 300 \textmu m & $170$\\
 2 & 1.8 \% alanine & 6 mL/hr & 300 \textmu m & $170$\\
 3 & 3 \% alanine & 6 mL/hr & 300 \textmu m & $170$\\
 4 & 3 \% alanine & 12 mL/hr & 300 \textmu m & $330$\\
 5 & 3 \% alanine & 50 mL/hr & 300 \textmu m & $1400$\\
 6 & 3 \% alanine & 70 mL/hr & 300 \textmu m & $1900$\\
 7 & 3 \% alanine & 70 mL/hr & 200 \textmu m & $1900$\\
 8 & 3 \% alanine & 70 mL/hr & 100 \textmu m & $1900$\\
 \hline
 \multicolumn{5}{|p{13cm}|}{ Peclet number ($Pe$) is defined as $\frac{Q}{wD}$ where $Q$ is the flow rate, $w$ is the channel width, and $D$ is the diffusion constant.} \\
 \hline
\end{tabular}
\caption{List of experimental conditions.}
\label{tp2}
\end{table}

\subsection{2.5 Data collection}
\addcontentsline{toc}{subsection}{Data collection}

The videos and images for the progress of oiling-out crystallization in the flow cell during the solvent exchange were accomplished by NIKON Eclipse Ni microscope with X-Cite Series 120 Q mercury light source. The NIKON software has an auto-whiting function which corrected the color saturation and an auto-exposure correction that adjusts the brightness of the videos and images. Both of these features were utilized for our experiments. We selected video frames with a fixed interval based on the rate in droplet and crystal growth.

The initial time $t_{0}$ in our experiments is defined as the time at the start of droplet formation on the substrate. There are few droplets shown in the bulk flow during the solvent exchange before $t_{0}$ but the widespread droplet formation on the substrate occurred separately. We based our time stamps only on the observation of the substrates through an optical microscope.

\section{3. Results and Discussion}

Effects of chemical patterns on the oiling-out crystallization are coupled with the influences from other conditions including hydrophobicity/hydrophilicity of the micro-patterns, flow rate, channel geometry, and solute concentration. We will show effects from chemical micro-patterns on smooth surfaces. Finally, we will demonstrate that solvent exchange can be used to induce crystallization of a solute other than alanine.

\subsection{3.1 Oiling-out crystallization on surfaces with chemical micro-patterns}
\subsubsection{a. Hydrophobic or hydrophilic micro-patterns }

During the solvent exchange at a flow rate of 12 mL/hr, a thin layer of liquid formed in the hydrophilic surrounding area while some small droplets formed on the hydrophobic micro-patterns as shown in Figure \ref{f6a} (a)(b).
The crystals evolved to a needle-like shape or a thin film similar to crystallization on homogeneous hydrophilic surfaces by solvent exchange\cite{zhang_2020_oilingout}. In the cases where the crystallization of a thin film occurred, small holes formed at the locations where the hydrophobic domains lied.

\begin{figure} [htp]
	\includegraphics[trim={0.5cm 0.5cm 1.5cm 0cm}, clip, width=1\columnwidth]{Figures/Updated/2.jpg}
	\caption[Oiling-out crystallization on the hydrophobic patterned substrate]{Oiling-out crystallization on the hydrophobic patterned substrate. a) At 12 mL/hr flow rate, b) the zoomed in images of the red highlighted square in a). The blue dotted line indicates the boundary of thin liquid film. c) At 70 mL/hr flow rate, d) the zoomed in images of the red highlighted square in c). Channel height: 300 \textmu m. Length of the scale bar: a) and c) 100 \textmu m and  b) and d) 50 \textmu m.
}
	\label{f6a}
\end{figure}

 At a faster flow rate of 70 mL/hr, droplets formed a regular array, and the holes became more apparent where the thin crystal film propagated along the surface (Figure \ref{f6a} (c)(d)). The position of holes copied the arrangement of patterns. Peng et al.\cite{peng_2017_morphological, peng_2017_simple} showed that the formation of droplets could be controlled by the patterned substrate due to preferential nucleation of oil droplets on hydrophobic areas. The alanine-rich sub-phase from oiling-out appeared to spread on the hydrophilic areas but not on the hydrophobic domains; therefore, the local rupture of the film led to the hole formation on the crystal film. This result suggested that we may make different hole arrangements on the crystal film by changing the pattern position. In particular, the size of crystal film is large at a faster flow rate, and the hole array is more regular. Such holes may allow gases, chemicals, or particles to permeate through the film.

\begin{figure} [htp]
	\includegraphics[trim={1.2cm 6.7cm 1.0cm 3cm}, clip, width=1\columnwidth]{Figures/Updated/3.jpg}
	\caption[Oiling-out crystallization on the hydrophilic patterned substrates.]{Oiling-out crystallization on the patterned substrates. a) on hydrophilic micro-patterns b) the zoomed in images of the red highlighted square in a). Length of the scale bar: a) 100 \textmu m and b) 50 \textmu m. The condition for solvent exchange is listed as No. 6 in Table \ref{tp2}.
}
	\label{f6b}
\end{figure}

The oiling-out crystallization behaviour on the substrate with hydrophilic pattern and hydrophobic background is shown in Figure \ref{f6b} (a)(b).
The droplets in hydrophilic domains are large with uniform droplet size, possibly due to the constraints from the patterns on the surface. Smaller droplets nucleated on the hydrophobic surrounding area with less uniform size distribution. The crystallization process propagated in the path of the droplet coalescence. The crystals are irregularly shaped, and the network of crystals is related to the hydrophilic patterns. The time from the initial crystallization to completion of crystallization in the field of view is around 50 s, almost twice the time on hydrophobic patterns in Figure \ref{f6a} (a)(b). A longer time for the crystal growth may be due to more droplets on the surface going through phase separation, requiring more solvent supplied from the flow.

\begin{figure} [htp]
	\includegraphics[trim={4.6cm 7.5cm 4.3cm 2.6cm}, clip, width=1\columnwidth]{Figures/Updated/4.jpg}
	\caption[Oiling-out crystallization of droplets on homogeneous vs patterned substrates]{Oiling-out crystallization of droplets on a) homogeneous surface, b) hydrophilic patterns and hydrophobic surrounding area. Length of the scale bar: 100 \textmu m. The condition for solvent exchange is listed as No. 3 in Table \ref{tp2}.
}
	\label{f6c}
\end{figure}

At a slower flow rate of 6 mL/hr in Figure \ref{f6c} (b), the hydrophilic patterns on the substrate had an even clearer effect on the droplet formation and crystallization. The position and shape of the droplets follow the pattern of the hydrophilic domains on the surface which in turn led to numerous and smaller crystals. In contrast, on a homogeneous APTES-coated substrate, the size of the crystals was much larger on average. Further studies on pattern sizes or different chemical patterns may allow optimization of solvent exchange for a specific crystal morphology and pattern.

\begin{figure} [htp]
	\includegraphics[trim={1.0cm 10cm 1.0cm 4.3cm}, clip, width=1\columnwidth]{Figures/Updated/5.jpg}
	\caption[Oiling-out crystallization with a low alanine concentration]{Oiling-out crystallization with the low alanine concentration. The condition for the solvent exchange is No. 2 in Table \ref{tp2}. Length of the scale bar: 100 \textmu m.
}
	\label{f6}
\end{figure}

The initial concentration of the solute also has an important effect on the droplet formation in the solvent exchange system. At an even lower concentration at 1.8 \% alanine in solution A, Figure \ref{f6} shows the decrease in crystals in comparison to the concentration of 3 \% in Figure \ref{f5} (a). However, the features in the droplet formation and the crystallization are similar to each other. The main differences are the smaller surface droplet coverage and the crystallization at a lower concentration of alanine, suggesting that there is a threshold of initial concentration required for oiling-out crystallization by solvent exchange.

\subsubsection{b. Flow rate of the poor solvent}

\begin{figure} [htp]
	\includegraphics[trim={1.7cm 3cm 1.5cm 2.6cm}, clip, width=1\columnwidth]{Figures/Updated/6.jpg}
	\caption[Oiling-out crystallization at different flow rate]{Oiling-out crystallization at the flow rate of a) 6 mL/hr, b) 12 mL/hr, and c) 70 mL/hr. The red circle indicates a case where a alanine-lean droplets are encapsulated by the alanine-rich phase. Length of the scale bar: 100 \textmu m.
}
	\label{f5}
\end{figure}

In Figure \ref{f5}, the channel height is kept constant while the flow rate varies. The droplets at lower flow rates lasted longer before crystallization and were more likely to dissolve instead of coming to full crystallization. The completion of droplet depletion took around 290 seconds for 6 mL/hr, 200 seconds for 12 mL/hr, and 85 seconds for 70 mL/hr. All three flow rates had a single droplet that covered the hydrophilic area, suggesting that the patterns effectively controlled the droplet size regardless of the flow rate.

At the slowest flow rate of 6 mL/hr, crystals on the surface had more polygonal shapes. Some crystals may agglomerate to form clusters of crystals. For 12 mL/hr, more flat crystals grew in diamond shape while still having few polygonal shapes. At the fastest flow rate of 70 mL/hr, the droplets crystallized in a network that spreads in a radiating fashion. These crystals were more irregular in shape, and the crystallization followed the path of the droplet coalescence. Different from crystallization on homogeneous surfaces, a faster flow rate led to large quantities of crystals of the droplet shape.\cite{zhang_2020_oilingout}

\begin{figure} [htp]
	\includegraphics[trim={1.7cm 7.3cm 1.5cm 2.6cm}, clip, width=1\columnwidth]{Figures/Updated/7.jpg}
	\caption[Reproducibility of oiling-out crystallization]{Oiling-out crystallization at the flow rate of 6 mL/hr. a) and b) are the repeats of the experiments at condition No. 3 in Table \ref{tp2}. Length of the scale bar: 100 \textmu m.
}
	\label{f5a}
\end{figure}

The general features in crystallization were reproducible under the same experimental conditions. As shown in Figure \ref{f5a}, some parts of the substrate or different experimental repeats may exhibit more flat crystal growth with diamond shapes. The crystals formed on patterned surfaces are much less polydispersed compared to those on a homogeneous surface. The number of crystals per $mm^{2}$ detected by the software came out to be 391 and 401 for Figure \ref{f5a} (a) and (b), respectively. Again, this indicated that general features of oiling-out crystallization were reproducible.

Faster crystallization at a faster flow rate of solvent exchange was explained by the transport of good solvent out from the droplets\cite{zhang_2015_formation, zhang_2020_oilingout}. In brief, water diffuses out from the droplets sooner at higher flow rates, leading to the supersaturation of alanine in the droplets in a shorter time.

The flow rate can also influence the crystallization by varying the number of the crystal seeds that land on the substrates, which is significantly greater at lower flow rates as shown in Figure \ref{f5} (a). More seeds led to smaller crystals at lower flow rates compared to larger connected crystals at higher flow rates.  Two reasons may be at play for more crystal seeds at lower flow rates. One may be due to a later crystal landing time. The droplets in the flow experienced counter diffusion, which triggers the crystallization as observed for spherical alanine crystal formation from droplets in work by Sun et al.\cite{sun_2018_understanding} The longer retention time for the ternary droplets in the flow may form more crystals in the flow before reaching a certain location on the substrate. Secondly, the slower flow rate allows more crystals to settle onto the surface in competition with advection by the bulk flow. Therefore, more crystal seeds can land on the substrate during the solvent exchange process.

\subsubsection{c. Height of the channel for solvent exchange}

 The crystallization from the solvent exchange performed in the channels with three different heights is shown in Figure \ref{fpd}. The substrates with hydrophilic micro-patterns surrounded by hydrophobic areas were used in the experiment (Figure \ref{fpd} (a)-(c)). At a given flow rate of 70 mL/hr, the droplet size decreased with a decrease in the channel height. There was only one droplet per micro-domain at the channel height of 200 \textmu m or 300 \textmu m, but multiple droplets at the channel height of 100 \textmu m.

 \begin{figure} [htp]
     \centering
         \includegraphics[trim={0.5cm 0cm 2cm 1cm}, clip, width=1\columnwidth]{Figures/Updated/8.jpg}

	\caption[Decreasing droplet size with decrease in channel heights]{Decreasing droplet size with decrease in channel heights. Droplets and crystals on the patterned substrate at a) 300 \textmu m, b) 200 \textmu m, and c) 100 \textmu m and on the homogeneous substrate at channel height of d) 300 \textmu m and e) 200 \textmu m. The red outline indicates the droplet outline in the hydrophilic domain. Length of the scale bar: 50 \textmu m
}
	\label{fpd}
\end{figure}

  As a comparison, effects from channel height were also examined on homogeneous APTES-coated substrate (Figure \ref{fpd} (d)(e)). Crystallization followed the droplet formation with the channel height of 200 \textmu m and 300 \textmu m. When channel height was down to 100 \textmu m, no droplet formation was observed. The droplets decreased in size with the decrease in channel heights. The dependence of droplet sizes on the channel height is consistent with that on the micro-patterned substrate.

The previous work showed that by solvent exchange, the final volume of the droplet $V_f$ increased with the channel height $h$\cite{zhang_2015_formation,yu_2015_gravitational}. At the same flow rate (in volume per unit time), $V_{f} \sim h^{3}$ because $h$ determines how long it takes for the solution A and B to mix uniformly across the channel. The higher the channel, the longer the time for the droplets to grow \cite{zhang_2015_formation,yu_2015_gravitational}.

Above certain channel height, the effect of gravity on the mixing between two solutions may also become important. Archimedes number $Ar$ is the dimensionless number to describe the gravity effect, arising from the density difference between solution A and the displacing solution \cite{yu_2015_gravitational}.

\begin{equation}
Ar = \frac{g (h)^3}{\nu^{2}}\frac{\Delta \rho}{\rho}
\label{eq1}
\end{equation}

The density of the solvents in solution A is 0.885 $g/mL$, while solution B is 0.786 $g/mL$. The density difference $\Delta$ $\rho$ is $\sim$ 0.1 $g/mL$, the gravitational acceleration $g$ is 9.8 $m/s^2$, and viscosity $\nu$ is $2.8 \times 10^{-6} m^{2}/s$. For our channel height $h$ of 100, 200, and 300 \textmu m, the $Ar$ number is 0.16, 1.3, and 4.3, respectively. For such a large Archimedes number ($Ar > 1$) at the channel height $h$ of 300 \textmu m, the center of the parabolic mixing front in the laminar flow is shifted towards the bottom surface under the gravity. Such shift would lead to a longer growth time for droplets and crystals on the top surface than on the bottom surfaces of a horizontal channel. However, the gravity effect could be eliminated by placing the channel vertically \cite{yu_2015_gravitational}.

Another effect that may arise in a higher channel is enhanced mixing from convection \cite{zhang_2015_formation,yu_2015_gravitational}.
To estimate when the convection rolls set in due to the density difference, we calculate Rayleigh number $Ra$.
\begin{equation}
Ra=\frac{\Delta \rho g (h/2)^3}{\mu D}
\label{eq2}
\end{equation}
where $\mu$ is the dynamic viscosity of solution A and the mass diffusion coefficient $D$ is approximately $10^{-5} cm^2/s$. For the channel height $h$ of 100, 200, and 300 \textmu m, $Ra$ number is 71, 570, 1900, respectively. The convection rolls may occur for 300 \textmu m channel heights where the Rayleigh number is slightly above the critical Rayleigh number 1708\cite{faber_1995_fluid}. These convection rolls enhanced mixing conditions and led to larger droplets along with the rolls than ones away from the rolls.

\subsection{3.2 Oiling-out crystallization on micro-structured surfaces}

 A substrate with the polymeric micro-lens array was used in our experiments as shown in Figure \ref{f2a} (a). These micro-lenses could be easily prepared by polymerization of surface droplets, therefore chosen as representative surface micro-structures to reveal the influence of physical structures on oiling-out crystallization. Figure \ref{f2a} b) from $t_0$ to $t_0$ + 60 s showed progression of droplet formation where $t_0$ is defined as a timestamp before observable oiling-out. With micro-lenses on the surface, the number density of droplets is much lower than that on a homogeneous surface, possibly due to fewer nucleation sites on the micro-lenses (Figure \ref{fpd} (d)). The droplets ranged from 5.6 \textmu m to 14.2 \textmu m, much larger than droplets on a homogeneous surface. The effect of the geometric structures on oiling-out droplets appears to be similar to that on the oil droplet formation by solvent exchange\cite{dyett_2018_coalescence,peng_2016_how}.

\begin{figure} [htp]
	\includegraphics[trim={1.7cm 4.3cm 1.5cm 2.6cm}, clip, width=1\columnwidth]{Figures/Updated/9.jpg}
	\caption[Oiling-out crystallization on the lens array]{Oiling-out crystallization on the lens array. a) an image of the lens array on the substrate b) the snapshot of the solvent exchange process. Length of the scale bar: a) 20 \textmu m and b) 40 \textmu m. Channel height: 300 \textmu m.
}
	\label{f2a}
\end{figure}

More importantly, the micro-lenses have a significant effect on the crystallization process shown in Figure \ref{f2a} (b). The crystals were easily detached from the substrate, and the crystal shapes were more irregular in comparison to those on the smooth substrates.\cite{zhang_2020_oilingout} The easier detachment of crystals may be attributed to the weak adhesion between the crystal plates and the spherical-cap shaped lenses on the substrate. In addition, the lenses caused the crystals to grow in a slightly elevated manner, which also favours crystal detachment. The self-detachment of crystals during the solvent exchange makes it possible for the collection of the crystals from the flow at the exit of the chamber. These collected crystals may be useful as seeds to trigger crystallization in the bulk crystallization\cite{zhang_2020_oilingout}.

\subsection{3.3 Fiber crystal formation from droplets without oiling-out}

\begin{figure} [htp]
	\includegraphics[trim={0cm 0cm 0cm 0cm}, clip, width=1\columnwidth]{Figures/Updated/10.jpg}
	\caption[Ternary phase of mesitylene, ethanol, and water and a diagram of fiber formation in surface droplet]{a) Ternary phase diagram of mesitylene, water, and ethanol from IUPAC-NIST solubility database\cite{skrzecz_1997_solubility} b) Sketch of Np-P4VB fiber formation during solvent exchange.
}
	\label{fdrawfib}
\end{figure}

Finally, we will show that the crystal formation by solvent exchange can be applied beyond the oiling-out systems. Np-P4VB was used as the solute for crystallization from droplets consisting of mesitylene, ethanol, and water. The solubility of Np-P4VB is the highest in ethanol and is the lowest in water. During the solvent exchange, solution A contained Np-P4VB, mesitylene, ethanol, and water, and solution B was water as presented in Figure \ref{fdrawfib} (a). Initially, surface droplets were mainly mesitylene with some ethanol and Np-P4VB. As the solvent exchange continued, the droplets became saturated with Np-P4VB as ethanol diffused out into the flow of water and crystallization of Np-P4VB was induced inside the droplets (Figure \ref{fdrawfib} (b)). The crystallization continued until the solute in the droplets was depleted.

\begin{figure} [htp]
\centering
	\includegraphics[trim={1.1cm 0.5cm 5cm 3cm}, clip, width=1\columnwidth]{Figures/Updated/11.jpg}
	\caption[Fiber formation in droplet on OTS substrate through solvent exchange]{Fiber formation in droplet at flow rate of a) 100 mL/hr on OTS substrate through solvent exchange. $t_0$ is defined as the time when the droplet formation started. Crystal fiber on the substrate after solvent exchange at flow rate of b) 12 mL/hr c) 75 mL/hr d) 100 mL/hr. Length of the scale bar: a) 100 \textmu m b-d) 250 \textmu m. e) Stitched images of crystal fiber formed for b). Total of 25 images were stitched together with 10 \% overlapping to manually stitch the images.
}
	\label{ffg}
\end{figure}

The crystal of Np-P4VB formed fibers during the solvent exchange. At a flow rate of 100 mL/hr, the fibers formed inside the droplets and extended with the solvent exchange process as shown in Figure \ref{ffg} (a). Different from the oiling-out crystallization of alanine, the crystallization initiated from the solute-rich droplets instead of the landed seed crystals, evident from the location of the crystal fibers. The onset of the crystals suggests that supersaturation may be from the reduced ethanol concentration in the droplets during the solvent exchange.

 The flow rate affected the crystal fiber formation, and longer fibers were formed at lower flow rates (Figure \ref{ffg} (b)-(d)). For a faster flow rate of 100 mL/hr, there was a significant decrease in fiber length (Figure \ref{ffg} (d)). To capture the bigger picture of the fiber formed on the substrate, the images of the fiber after solvent exchange at 12 mL/hr flow rate have been stitched as a larger image (Figure \ref{ffg} (e)). The results demonstrate that solvent exchange can be applied to control crystallization from droplets containing the solute, not limited to oiling-out systems.

\section{Conclusion}

We show that surface micro-patterns can be an effective way to control oiling-out crystallization during solvent exchange. A thin film of crystals can be created with holes on the substrates with hydrophobic micro-patterns. More uniform crystals form on the surfaces with hydrophilic micro-domains, as compared to a homogeneous surfaces. The channel height and flow rate of solvent exchange can be varied to control the oiling-out droplets and crystallization. By using representative micro-structures of polymeric micro-lenses on a substrate, the crystals detached easily from the surface, which may be used for collecting seed crystals to trigger crystallization in the bulk crystallization. Beyond oiling-out systems, solvent exchange can be used to induce crystallization by forming droplets containing the solute that does not exhibit oiling-out behaviour.

The results in the study show great potential in using solvent exchange for crystallization from droplets on the solid surface. Further studies in other chemicals will be valuable for using solvent exchange as a way to control, separate, and purify the crystal product for pharmaceutical and many other applications.

\section{Acknowledgements}

The project is supported by the Natural Science and Engineering Research Council of Canada (NSERC) and Future Energy Systems (Canada First Research Excellence Fund). This research was undertaken, in part, thanks to funding from the Canada Research Chairs program. We also thank Al Meldrum and Sergey Vagin for generously supplying us with Np-P4VB for our experiments.
\cite{codan_2010_phase,codan_2012_phase}


\subsubsection*{References}}
\usepackage{mathtools}
\usepackage{booktabs}
\usepackage{tikz}

\newcommand{\swap}[3][-]{#3#1#2}

\title{Adaptive Belief Discretization for POMDP Planning}

\author[1]{\href{mailto:Harry Q. Bovik <harryq@example.edu>?Subject=Your UAI 2021 paper}{Divya Grover}{}}
\author[1,2]{Christos Dimitrakakis}
\affil[1]{%
    Department of Computer Science and Engineering\\
    Chalmers University, Sweden
}
\affil[2]{%
    Department of Informatics\\
    University of Oslo, Norway
}

\begin{document}
\maketitle

\begin{abstract}
Partially Observable Markov Decision Processes
(POMDP) is a widely used model to represent the
interaction of an environment and an agent, under
state uncertainty. Since the agent does not observe
the environment state, its uncertainty is typically
represented through a probabilistic belief. While
the set of possible beliefs is infinite, making exact
planning intractable, the belief space’s complexity
(and hence planning complexity) is characterised
by its covering number. Many POMDP solvers
uniformly discretize the belief space and give the
planning error in terms of the (typically unknown)
covering number. We instead propose an adaptive
belief discretization scheme, and give its associated
planning error. We furthermore characterise the
covering number with respect to the POMDP parameters. This allows us to specify the exact memory requirements on the planner, needed to bound
the value function error. We then propose a novel,
computationally efficient solver using this scheme.
We demonstrate that our algorithm is highly competitive with the state of the art in a variety of
scenarios.
\end{abstract}

\section{Introduction}

We are interested in sequential decision making problems, where an agents is interacting with a partially observable environment, meaning that it is unaware of the actual state of the environment. This interaction can be formalised through a Partially Observable Markov Decision Process (POMDP)~\cite{smallwood1973optimal}. Since the agent does not know the true state of the environment, it must instead maintain a subjective belief representing its uncertainty about the current environment state. This is typically expressed as a probability distribution over the states.

Online POMDP algorithms plan ahead by starting from the current belief and building a tree-like structure in the space of possible future beliefs, by enumerating different environment responses to the agent's actions. However, exact solutions to some lookahead horizon require an exponential size tree with respect to that horizon. Typical approximations include action selection heuristics that only consider promising  actions and leaf node approximations, that truncate the tree and replace the value of the leaf nodes with heuristics, or upper or lower bounds. While both these techniques are effective when building shallow trees, we instead propose an online POMDP algorithm that builds deep trees. In particular, we propose an adaptive discretization scheme that generates a cover with a small size, thus reducing memory and computational complexity.

\subsection{Related work}
Current state-of-the-art online POMDP solvers include POMCP~\cite{silver2010monte} and DESPOT~\cite{ye2017despot}. Like many older~\cite{ross2007bayes} and newer~\cite{kurniawati2016online} works they are based on inference using particle filtering. In building the lookahead tree, they keep its size in check by only considering ‘promising’ sequence of actions, by maintaining a statistical estimate of the action values (known as Upper Confidence Bound). DESPOT in addition, uses importance sampling to not miss important beliefs. But neither algorithm gives a structured way to check the approximate inference's effect on the planning error. While the belief discretization we propose results in a form of approximate inference as well, the focus here is on how the approximations made in the inference and in the planning process \emph{interact}.

A different method for approximate inference in POMDPs includes  belief representation methods like ~\cite{roy2005finding}. These implicitly define a kernel function over the belief space, but give no bound on the resultant solver's planning error.

One concept we use in our theoretical analysis is the covering number of the belief space. It is a measure of the density of a set of points in any continuous space. Only recently has it been used as a complexity measure for the POMDP. In particular,~\citet{zhang2012covering} performs an experimental estimation of the covering number for various environments, by using breadth first search (BFS). However, unlike this paper, the authors do not provide a POMDP solver or derive any dependence between the underlying POMDP parameters and the covering number of the belief space.

Covering number was also used in the analysis of SARSOP~\cite{kurniawati2008sarsop} by ~\cite{lee2007makes}. In particular, it shows how the planner's error depends on the (unknown) covering number. Such solvers still depends exponentially on the lookahead, as noted by~\cite{zhang2012covering} (refer Lemma 1). Finally,~\citet{zhang2014covering} also propose PGVI, which like SARSOP, uses depth first search (DFS) to estimate the optimal action. Analytically, they are the same, since PGVI's complexity is identical to that of SARSOP.\footnote{Verify that the tuned parameters are set to identical values in both works.}

\paragraph{Contributions:}
\begin{enumerate}
    \item We propose explicit lookahead dependent covers for the belief space. This significantly reduces the memory requirement on the planner.
    \item We give the exact dependence of the planning error on this adaptive discretization scheme.
    \item We give an upper bound on the memory requirement of the planner by giving bounds on the covering number of the belief space.
    \item We propose a novel, computationally efficient planner with direct control on the planning error with respect to its tuning parameters.
\end{enumerate}

\paragraph{Organization.} In Section 2, we define the POMDP model, state the used equations and introduce the concept of cover sets. In Section 3, we follow up with analysis, giving a proof of sketch first, then stating our main result, and finally the proofs. In Section 4, we describe our proposed algorithm and Section 5 is dedicated to experiments.

\section{Setting}
We divide the section into three parts: First giving the definitions of POMDP model, second giving the equations used to compute various quantities of interest in a POMDP and third part introducing the concept of $\epsilon$-cover.
\subsection{Model Definition}
\begin{definition}[MDP]
A Markov Decision Process (MDP) is a discrete-time stochastic process that provides a formal framework for decision making agents. An MDP $\mu \defn (S,A,T,R,\disc)$ is composed of a state space $S$, an action space $A$, a reward distribution $R$ and a transition function $T$. The transition function $T \defn \Pr(s_{t+1}|s_t ,a_t)$ dictates the distribution over next states $s_{t+1}$ given the present state-action pair $(s_t, a_t)$. The reward distribution $R \defn \Pr(r_{t+1}|s_t ,a_t)$  dictates the obtained reward with $r \in [0,\reward_{\max}]$.
\end{definition}
The discount $\disc \in [0,1)$ signifies the effective horizon $\horizon \defn \frac{1}{1 - \disc}$ of the agent's interaction with its environment. \emph{Henceforth, we consider only finite state, action and observation spaces}.
\newline

In many interesting settings, the agent can only indirectly estimate the state $s_t$ via a set of observations $o_t\in O$ and their (state-dependent) emission probability. Formally,
\begin{definition}[POMDP]
A Partially Observable MDP, or POMDP, denoted by $(S,A,T,R,,\disc,O,Z)$ tuple is an MDP model augmented with an observation $o_t\in O$ set and their corresponding emission probabilities $Z\defn\Pr(o_t|s_t)$.
\end{definition}
In such cases, the best an agent could hope to achieve is to act optimally in an expected sense (averaging over its estimation of the current state). Formally, the objective is to find action maximizing
\begin{equation}
 V^*(b_0) \defn \lim_{\horizon \rightarrow \infty} V_{\horizon}(b_0)
\end{equation}
where $V_{\horizon}(b_0)$, called the belief value function, is defined as
\begin{align}

	V_{\horizon}(b_0)
	=
	\max_a \mathbb{E} ( \sum_{k=0}^{H-1} \disc^{k} r_{k+1}|b_0,a) \label{eqn:val}

	\\
	b_0 \in \Delta^{|S|} \notag \\
	\Delta^{|S|} \defn \{ b \in {[0,1]}^{|S|} \; | \; \sum_s b(s) = 1 \}
\end{align}
Where $b_0$ denotes its current state estimation, known as belief. Formally, belief is a probability distribution over the state space whose element $b(s)$ gives the probability that the system’s true state is $s$. It is a compact representation of the agent's complete history of interaction with its environment, i.e, the action-observation sequence it followed.
\newline

Figure~\ref{fig:tree} shows the interaction of an agent with its environment. At each time step, it evaluates the best possible action by anticipating the effect of a sequence of actions on future belief and reward. This process of anticipation is called planning and the corresponding tree-like structure is called the belief tree. Although in Figure~\ref{fig:tree}, we only calculate effect of a fixed action, you can still note the exponential growth of the number of nodes with depth.
\subsection{Equations}
Here we describe the \textbf{standard} computational equations of the POMDP model.
To compute a new belief $b_{k+1}$, given the current belief $b_k$, action $a$ and next observation $o_{k+1}$, one uses the following set of equations:
\begin{align}
   \Pr(s_{k+1}) = \sum_{s'} T(s_{k+1}|s',a) b_k(s') \label{eq:marginal}\\
    b_{k+1}(s_{k+1}) = \frac{1}{\eta} \; Z(o_{k+1}|s_{k+1}) \Pr(s_{k+1}) \label{eq:posterior}\\
    \eta \defn \sum_{s_{k+1}} Z(o_{k+1}|s_{k+1}) \Pr(s_{k+1}) \notag
\end{align}
collectively referred to as the \emph{Bayes filter}. Equation~\ref{eq:marginal} is state marginal distribution after transition has occurred under the given action, while equation~\ref{eq:posterior} gives the Bayesian posterior (called the next belief) of the states.
\newline
Solving a POMDP involves first generating the belief tree, then doing backward induction~\cite{bellman1952theory} on it to calculate the optimal action at the root node. The corresponding Bellman equation for POMDP is:
\begin{equation}\label{eq:backward-induction}
 V(b_k)
	=
	\max_a \mathbb{E}_{\Pr(r_{k+1},o_{k+1}|a,b_k)} ( \; r_{k+1} + V(b_{k+1}) \; )
\end{equation}
Where
\begin{align*}
    \Pr(r_{k+1},o_{k+1}|a,b_k) &= \Pr(o_{k+1}|a,b_k) \Pr(r_{k+1}|a,b_k) \\
    \Pr(r_{k+1}|a,b_k) &= \sum_{s'} R(r_{k+1}|s',a) b_k(s') \\
    \Pr(o_{k+1}|a,b_k) &= \eta \\
\end{align*}
Note that $b_{k+1}$ is uniquely determined by $o_{k+1}$ using Bayes filter.

\begin{figure*}[th!]
		\centering
			\begin{tikzpicture}
			\node[RV] at (-5,0) (bt) {$b_t$};
			\node[RV] at (5,0) (bt1) {$b_{t+1}$};
			\node[select] at (-5,-1) (ah) {$a$};
			\node[select] at (5,-1) (ah1) {$a$};
			\draw[->] (bt) to (ah);
			\draw[->] (bt1) to (ah1);

			\node[RV] at (-4,-2) (bo1) {$b_{o_2}$};
			\node[RV] at (-6,-2) (bo2) {$b_{o_1}$};
			\draw[->] (ah) to (bo1);
			\draw[->] (ah) to (bo2);

			\node[RV] at (6,-2) (bo11) {$b_{o_2}$};
			\node[RV] at (4,-2) (bo12) {$b_{o_1}$};
			\draw[->] (ah1) to (bo11);
			\draw[->] (ah1) to (bo12);

			\node[select] at (-4,-3) (ahp1) {$a$};
			\node[select] at (-6,-3) (ahp2) {$a$};
			\draw[->] (bo1) to (ahp1);
			\draw[->] (bo2) to (ahp2);
			\node[RV] at (-5,-4) (bo3) {$b_{o_4}$};
			\node[RV] at (-7,-4) (bo4) {$b_{o_3}$};
			\draw[->] (ahp1) to (bo3);
			\draw[->] (ahp2) to (bo4);
			\draw[dashed,->] (ahp1) to (-3,-4);
			\draw[dashed,->] (ahp2) to (-5.5,-4);

			\node[select] at (6,-3) (ahp11) {$a$};
			\node[select] at (4,-3) (ahp12) {$a$};
			\draw[->] (bo11) to (ahp11);
			\draw[->] (bo12) to (ahp12);
			\node[RV] at (7,-4) (bo13) {$b_{o_4}$};
			\node[RV] at (3,-4) (bo14) {$b_{o_3}$};
			\draw[->] (ahp11) to (bo13);
			\draw[->] (ahp12) to (bo14);
			\draw[dashed,->] (ahp11) to (5.5,-4);
			\draw[dashed,->] (ahp12) to (4.5,-4);

			\draw[thick,->] (bt) -- (bt1) node[midway,above] {$t$};
			\node[RV] at (0,-2) (h1) {$d = 1$};
			\draw[->] (0,0) to (h1);
			\draw[dashed,-] (h1) to (-2,-2);
			\draw[dashed,-] (h1) to (2,-2);

			\node[RV] at (0,-4) (h2) {$d = 2$};
			\draw[->] (h1) to (h2);
			\draw[dashed,-] (h2) to (-2,-4);
			\draw[dashed,-] (h2) to (2,-4);
			\end{tikzpicture}
			\caption{Online POMDP planning}\label{fig:tree}
\end{figure*}

\subsection{Cover sets and covering number}
For a set $\mathcal{X} \subset \mathbb{R}^n$, define metric $d_p: x \times y \rightarrow \mathbb{R} \; | \; x,y\in \mathcal{X}$:
$$ d_p(x,y) \defn {\left( \sum^n_{i=1} |x_i - y_i|^p \right)}^{\frac{1}{p}} $$
\begin{definition}
Given a metric space $(\mathcal{X},d_p)$. Define $\epsilon$-ball of a point $x_0 \in \mathcal{X}$ as $$c(x_0,\epsilon) \defn \{ x \in \mathcal{X} \; | \; d_p(x_0,x) < \epsilon \}$$
Then for a given set $B \subset \mathcal{X}$, define an $\epsilon$-cover of $B$ as
$$C_B(\epsilon) \defn \{x \in B \; | \; B \subset \cup \; c(x,\epsilon) \}$$
Similarly an improper cover of $B$ is defined as
\begin{align*}
    \hat{C}_B(\epsilon) \defn \{x \in \mathcal{X} \; | \; B \subset \cup \; c(x,\epsilon) \}
\end{align*}
Therefore the elements of an improper cover do not necessarily have to come from set $B$ itself.\\

The covering number is defined as $$\mathcal{N}(\epsilon,B,d_p) \defn \min  |C_B(\epsilon)| $$ i.e, the minimum cardinality sets amongst such sets.\\

An improper covering number is similarly defined
$$\hat{\mathcal{N}}(\epsilon,B,d_p) \defn \min  |\hat{C}_B(\epsilon)|$$
A standard result used here is:
\begin{equation}
\mathcal{N}(2\epsilon,B,d_p) \leq \hat{\mathcal{N}}(\epsilon,B,d_p) \leq \mathcal{N}(\epsilon,B,d_p) \label{eq:improperCover}
\end{equation}
\end{definition}

Note that we refer to an $\epsilon$-cover as a cover set when the value $\epsilon$ is either implied or not important in the context.

\section{Analysis}
Our analysis relies on creating a sequence of covers for the belief, so that covers further away in the future are sparser. We show that we can uniformly bound the value function error by an appropriate increase in the sparsity that depends on the depth and discount factor.

Our main theorem, the \emph{planning error theorem}, gives a bound on the number of representative beliefs we need to keep track of, to bound the value function error at the root. This is accompanied by a lemma bounding the memory required by our discretization scheme and algorithm.

In particular, for exact $\horizon$-step lookahead planning, the corresponding belief tree would have order $O({|A|}^\horizon{|O|}^\horizon)$ nodes. This is prohibitive to store in memory even for a moderate horizon. A remedy to this problem comes from the fact that belief values functions are Lipschitz continuous.
While previous work~\cite[c.f.][]{lee2007makes} used this property to obtain a uniform cover, in our work we make the cover sparser as we go deeper in the tree.

The Lipschitz property allows us to control tree size by expanding only a certain representative node's subtree and substituting its value for all beliefs represented by that particular belief node. We can then relate the approximation errors at different depths through discounting~\cite[c.f.][]{kearns2002sparse}.
Combined together these ideas gives the required planning error theorem.

\subsection{Main Result}
\begin{definition}\label{def:discretization}
A $\epsilon$-discretization of a belief space $\Delta^{|S|}$ is defined\footnote{Ignoring the normalization constant.} as:
$$A(\epsilon) = \{ \sum_{i=1}^{S} \alpha_i {\mathbb{1}}_i | \; \alpha_i \in \{0,\epsilon,2\epsilon,\ldots,1 - \epsilon,1\}  \}$$
i.e, a vector space with ${\mathbb{1}}_i$ as basis and coefficients $\alpha_i$. ${\mathbb{1}}_i$ denotes a unit vector with ${\mathbb{1}}_i(j) = 0 \; | j\neq i$.
\end{definition}
The error between root belief's value function to its true value function is given by the planning error theorem:
\begin{theorem}\label{theorem-1}
For any belief $b_0$, and target error $\epsilon_v$, there exists a belief tree $\mathcal{T}$ of height $h_0$ rooted at $b_0$, with separate $\epsilon_d$-discretization of the belief spaces at each depth $d$, given by the sequence
$\epsilon_d \defn \frac{\epsilon_{b_0}}{\disc^d}$, such that $b_0$'s value $\hat{V}(b_0)$ at the root satisfies
$$|\hat{V}(b_0) - V^*(b_o)| \leq \epsilon_v$$
for parameter values:
\begin{align*}
h_0 = \log_{\disc} \frac{ {(1-\disc)} \epsilon_v}{2 \reward_{\max}} \\
\epsilon_{b_0} = \frac{ {(1-\disc)} \epsilon_v}{2 \disc \reward_{\max} h_0}
\end{align*}
\end{theorem}
Contrast $\epsilon_{b_0}$ with the value of parameter $\delta (= \frac{ {{(1-\disc)}^2} \epsilon_v}{2 \disc \reward_{\max} } )$ in ~\cite{lee2007makes}. Both denote the level of discretization required. While our scheme demands a finer initial discretization ($\epsilon_{b_0}$) at the root (by a factor of $\frac{1}{h_0 (1-\disc)}$), it quickly (because $\epsilon_d \defn \frac{\epsilon_{b_0}}{\disc^h}$) becomes coarser (better) as we expand the tree deeper.
\newline
Now we state our Lemma giving total memory required to build the discretized belief tree $\mathcal{T}$.
\begin{lemma}\label{lemma:memory}
For a tree $\mathcal{T}$ of total height $h_0$ with separate $\epsilon_d$- discretized belief spaces at each level $d$, as defined in Theorem~\ref{theorem-1}. The total memory required $M_{h_0}$ is bounded by:
$$M_{h_0} \leq {\ceil{\frac{2}{\epsilon_{b_0} (1-\disc)}}}^S$$
\end{lemma}
Note that the memory requirement of SARSOP and PGVI is the covering number of the minimal $\delta$-cover for their specified tolerance value $\delta$~\cite{lee2007makes}. But they don't give an upper bound on its value, because the corresponding cover sets are only implicit.
\subsection{Proofs}
Since most proofs depend on backward induction equation~\ref{eq:backward-induction}, we use the following indexing convention: For a belief tree of total height $h_0$, we define its height by index $i$, hence $i=0$ denotes the leaf nodes level, while $i=h_0$ denotes the root level. Also for depth $d$, the relation $d = h_0 - i$ holds. This convention is also followed in the previous literature.
\newline
We first give a Lemma on Lipschitz continuity of value function at different heights:
\begin{lemma}[Lipschitz continuity]\label{lemma:continuity}
For any two belief nodes $b$ and $b'$ at height $i$, we have
$$ {||V_i (b) - V_i (b')||}_1 \leq L_i \; \delta_i$$
where
\begin{align*}
L_i \leq \mathcal{R}_{\max} \frac{1-\disc^i}{1-\disc} + \disc^i L_0 \\
\delta_i = {||b-b'||}_1
\end{align*}
\end{lemma}
Proof is a simple consequence of H$\ddot{o}$lder's inequality.
\begin{align*}
    {||V_i (b) - V_i (b')||}_1 &\defn \\
    &= \sum_s V_i(s) (b(s)-b'(s)) \\
    &\leq {||V_i(s)||}_{\infty} \times {||b(s)-b'(s)||}_1 \\
    &= \max_s V_i(s) \times \delta_i \\
    &= L_i \delta_i
\end{align*}
where second line follows from H$\ddot{o}$lder's inequality.
\newline
Hence by using equation~\ref{eq:backward-induction} we get
\begin{align*}
L_i &\defn \max_s V_i \\
\max_s V_i &= \max_s \; ( r + \disc V_{i-1} ) \\
\max_s V_i &=  \mathcal{R}_{\max} + \disc \max_s V_{i-1}
\end{align*}
which applied recursively gives the result $\square$
\newline

We next give a Lemma relating value function error across different heights:
\begin{lemma}\label{lemma:error}
Maximum error for any node at height $h$, denoted by $\epsilon_h$ is given by:
$$ \epsilon_h \leq \sum^{h-1}_{i=0} \disc^{h-i} L_i \delta_i + \disc^h \epsilon_0$$
where $\epsilon_0$ denotes the error at leaf nodes.
\end{lemma}
Proof. For a belief $b$ at level $i=0$, consider its error in value function $\epsilon(b)$ as:
\begin{align*}
\epsilon(b) &= |\hat{V}(b) - V^*(b)| \\
&\leq |\hat{V}(b) - V(b')| + |V(b') - V^*(b)| \\
&\leq L_0 \delta_0 + \epsilon_0
\end{align*}
Where line 2 comes from the fact that it maybe the case that value of $b$ is substituted by value of another belief $b'$, constituting the $\delta_0$-cover at level $i=0$. Moreover, for any parent of $b$ (say $b_i$), at level $i=1$, its error $\epsilon_1 \leq \disc \epsilon(b)$~\cite{kearns2002sparse}. Hence by recursion, we get:
\begin{align*}
    \epsilon_1 &\leq \disc (L_0 \delta_0 + \epsilon_0) \\
    \epsilon_2 &\leq \disc (L_1\delta_1 + \disc L_0 \delta_0 + \disc \epsilon_0) \\
    &\ldots \\
    \epsilon_h &\leq \sum^{h-1}_{i=0} \disc^{h-i} L_i \delta_i + \disc^h \epsilon_0
\end{align*}
Where $\epsilon_i$ denotes the error in value function for any belief at level $i$.
\paragraph{Proof of Theorem~\ref{theorem-1}:} Equating root error in Lemma~\ref{lemma:error} to target error $\epsilon_v$ gives:
$$ \sum^{h-1}_{i=0} \disc^{h-i} L_i \delta_i + \disc^h \epsilon_0 = \epsilon_v $$
Combining this with Lemma~\ref{lemma:continuity} gives:
$$ \sum^{h-1}_{i=0} \disc^{h-i} \left( \reward_{\max} \frac{1-\disc^i}{1-\disc} + \disc^i L_0 \right) \delta_i + \disc^h \epsilon_0  = \epsilon_v$$
Assuming $L_0 = \frac{\reward_{\max}}{1-\disc}$, and leaf error $\epsilon_0 = \frac{\reward_{\max}}{1-\disc}$, we get
\begin{align*}
\sum^{h-1}_{i=0} \disc^{h-i} \epsilon_0 (1-\disc^i + \disc^i) \delta_i + \gamma^h \epsilon_0 &= \epsilon_v\\
\disc^h \epsilon_0 \left( \sum^{h-1}_{i=0} \frac{\delta_i}{\disc^i} + 1 \right) &= \epsilon_v \\
\left( \sum^{h-1}_{i=0} \frac{\epsilon_{b_0}}{\disc^{h-i-1} \disc^i} + 1\right)&= \frac{\epsilon_v}{\disc^h \epsilon_0}\\
\frac{h \epsilon_{b_0}}{\disc^{h-1}} + 1 = \frac{\epsilon_v}{\disc^h \epsilon_0}
\end{align*}
Where we use the fact that
\begin{align*}
\delta_i &\defn \epsilon_d = \frac{\epsilon_{b_0}}{\disc^d} = \frac{\epsilon_{b_0}}{\disc^{h-i-1}}
\end{align*}
Finally, substituting
$$h=h_0=\log_{\disc} \frac{\epsilon_v}{2\epsilon_0} = \log_{\disc} \frac{\epsilon_v (1-\disc)}{2 \reward_{\max}}$$
gives $\epsilon_{b_0} = \frac{\disc^{h_0-1}}{h_0} \; \square$
\newline

We now calculate the total memory required to store a discretized tree $\mathcal{T}$ by proving Lemma~\ref{lemma:memory}. Note that $A(\epsilon)$ acts as $\epsilon$-cover for the set ${[0,1]}^{|S|}$. Since belief space $\Delta^{|S|} \subset {[0,1]}^{|S|}$, $A(\epsilon)$ also forms an improper cover for it. We bound the covering number for each depth $d$, by bounding the cardinality of the corresponding sets $A(\epsilon_d)$. We consider the $\epsilon_d \defn \frac{\epsilon_b}{\disc^d}$-discretization of any arbitrary sequence of ${[0,1]}^{|S|}$ spaces.
\paragraph{Proof of Lemma~\ref{lemma:memory}:} For any $h_0$ length sequence of discretized spaces, total memory required:
\begin{align*}
\hat{M}_{h_0} &\defn \sum_d \mathcal{N}(\epsilon_d,A(\epsilon_d),d_1) \\
&= \sum_d |A(\frac{\epsilon_b}{\disc^d})| \\
log \sum_d |A(\frac{\epsilon_b}{\disc^d})| &\leq |S| \sum_d \log(1 + \floor{\frac{\disc^d}{\epsilon_b}}) \\
&\leq |S| \sum_{d=0}^{h_0} \floor{\frac{\disc^d}{\epsilon_b}} \:\: | \; \disc^{h_0} \geq \epsilon_b \\
&\leq |S| \sum_{d=0}^{h_0} \frac{\disc^d}{\epsilon_b} \\
\hat{M}_{h_0} &\leq {\ceil{\frac{1}{\epsilon_b (1-\disc)}}}^{|S|}
\end{align*}
From equation~\ref{eq:improperCover}
$$\mathcal{N}(\epsilon,\Delta^{|S|},d_1) \leq \mathcal{N}(\epsilon/2,{[0,1]}^{|S|},d_1)$$
Combing the two results proves the Lemma $\square$

\section{Finite Memory Planner}
We now introduce our main algorithm, Finite Memory Planner (FMP). This is an approximate online POMDP algorithm, with the following two main steps:
\begin{enumerate}
    \item Building an approximation $\mathcal{T}$ to the full belief tree, achieved here by recursively calling Algorithm~\ref{algo1}.
    \item Calculating the optimal action at the root and its value function using equation~\ref{eq:backward-induction-modified} on $\mathcal{T}$.
\end{enumerate}
The lookahead algorithm for building the belief tree $\mathcal{T}$ is described in Algorithm~\ref{algo1}.

\begin{algorithm}[t!]
		\caption{FMP (Finite Memory Planner) \\ {\bfseries Parameters:} Horizon $h_0$, discretization levels $\epsilon_d$}
		\begin{algorithmic}[1]
			\STATE {\bfseries Input:} Nodes with possible beliefs $b_d\in B_d$ at depth $d$,\\ tolerances $\epsilon_d$, current height $h$
			\IF{$h=h_0$}
			\STATE return []
			\ENDIF
			\STATE children = []
			\FOR{all $b_d \in B_d$}
			\FOR{all actions a}
			\FOR{all next observations $o_{d+1}$}
			\STATE $b_{d+1},o_{prob},r =$ BayesFilter$(b_{d},o_{d+1},a)$
			\label{inference-step}
			\STATE $\hat{b}_{d+1} = NN(b_{d+1},A(\epsilon_{d+1}))$\label{nn-step}
			\STATE children.append$([\hat{b}_{d+1},o_{prob},a,r])$
			\ENDFOR
			\ENDFOR
			\ENDFOR
			\STATE FMP(children,$\epsilon_{d+1}$,$h+1$)
		\end{algorithmic}
		\label{algo1}
	\end{algorithm}

\paragraph{Building the tree.}
The first stage is building a tree to a fixed depth, as follows.
\begin{enumerate}
    \item For each belief $b_d$ at depth $d$, action $a$ and observation $o$, we call ‘BayesFilter’ function at line~\ref{inference-step}, which uses equations~\ref{eq:posterior} and \ref{eq:marginal} to compute its true posterior and other attributes. They are later needed to compute the root value and optimal action using backward induction (equation~\ref{eq:backward-induction-modified}).
    \item At line~\ref{nn-step}, we replace the true posterior $b_{d+1}$ by its nearest neighbour in the $A(\epsilon_d)$ set which serves as its $\epsilon_d$-cover. This is computationally cheap because of uniform discretization (refer Definition~\ref{def:discretization}).
\end{enumerate}
We remark that the computational cost of Bayesian inference (line~\ref{inference-step}) is a shared burden for all POMDP solvers, and is not any more expensive for FMP, than for other solvers. One possible advantage of FMP is that inference can be quite flexible, e.g., one may use exact, variational inference or sampling based inference at line~\ref{inference-step}.\footnote{Contrast this to POMCP, which is solely based on Monte Carlo (particle filter) inference.}
Finally, we make a note that even $\epsilon$-cover based algorithms SARSOP and PGVI are significantly different. In contrast to FMP, both perform DFS to build the tree, thereby forgoing any chance of building a minimal $\epsilon$-cover at each depth. Computationally, FMP is cheaper as it does not need to cluster points at any level, while both SARSOP and PGVI do.
\paragraph{Calculating the value.}
In our scheme, the subsequent beliefs $b_{k+1}$ of any belief $b_k$ are replaced by some approximations $\hat{b}_{k+1}$ in the $\epsilon_{k+1}$-cover. Hence, the backwards induction equation is now
\begin{equation}\label{eq:backward-induction-modified}
 V(b_k)
	=
	\max_a \mathbb{E}_{\Pr(r_{k+1},o_{k+1}|a,b_k)} ( \; r_{k+1} + V(\hat{b}_{k+1}) \; )
\end{equation}
When the observation space is large we approximate the above by sampling $o_{k+1} \sim \Pr(o_{k+1}|a,b_k)$ and $r_{k+1} \sim \Pr(r_{k+1}|a,b_k)$.\footnote{In particular,  we sample $(o_{k+1},r_{k+1})$ in \emph{RandomPOMDP} environment but only $(r_{k+1})$ in \emph{RockSample} environment since it has binary observation space.}
Note that the algorithm returns the complete approximate value function for every depth $k+1$ to the previous depth $k$.

\section{Experiment}
We compare with State-of-the-art planner POMCP~\cite{silver2010monte} mainly, on RandomPOMDP and RockSample environments. We also compare with PBVI~\cite{pineau2003point} on RandomPOMDP and DESPOT~\cite{ye2017despot} on RockSample environment. The environment description is as follows:
\begin{enumerate}
    \item \textbf{RandomPOMDP ($RP[s,sp]$):} A random POMDP of given state, action and observation space. The reward, transition and observation distributions is generated using a sparsity parameter $sp \in \{0.3,0.6,0.9\}$ with $\mathcal{R}_{\max} = 1$. We experiment with different state space sizes $|S|$ and sparsity parameters. We fix $|Z|=|S|$ and $|A|=4$. Denote a particular configuration by $RP[s,sp]$.
    \item \textbf{RockSample ($RS[n,k]$):} Environment proposed in ~\cite{smith2012heuristic}, its state space size ranges from 12 thousand to 7 million. Although the belief space is very sparse, consisting of only 8 to 15 partially observed quantities. RS[n,k] denotes a map of size $n \times n$, searched by a agent, looking to collect good rocks between k possible total rocks. Its has a binary observation space, with distance dependent observation likelihood, giving more accurate observations as the agent moves closer to the rocks. Its goal is to collect the good rocks and exit the map. The agent knows its location and that of rocks, but is unaware of the state (good/bad) of each rock.
\end{enumerate}
\subsection{setup}
We take discount $\disc=0.95$ (standard) and run each experiment for $90$ steps\footnote{Inline with POMCP and DESPOT paper for RockSample, and sufficient for RandomPOMDP (since $0.95^{90} \approx 0.01$)}. We perform 50 experiments for each (Planner,POMDP) combination. We use POMDPy library~\cite{emami2015pomdpy} for RockSample simulator and PyPOMDP library~\cite{PyPOMDP} to solve our generated RandomPOMDP environment with POMCP and PBVI. While for the RockSample environment, we simply report the published (and tuned) results for DESPOT and POMCP (Table.1 in  ~\cite{ye2017despot}). Like \cite{smith2012point}, we mask the known states in RockSample environment. Finally, we also use GNU parallel~\cite{tange2018gnu} to run parallel experiments.
\paragraph{Parameter tuning:} It is worth noting that there are minimal tuning parameters to begin with, just $h_0$ and the sequences $\{\epsilon_d\}$. The solution quality is fairly linear with fineness of discretization, and inversely proportional to the error tolerance levels. In the experiments, the values are selected such that computation remains within satisfactory online planning limit (5 second for all environments up to RS[11,11], while 30 seconds for RS[15,15]). In particular, the following values were selected: $\epsilon_d = \{\epsilon_0 = \frac{1}{5},\ldots,\epsilon_H = \frac{1}{5}\}$ or $\{\epsilon_0 = \frac{1}{5},\epsilon_1 = \frac{1}{4},\ldots,\epsilon_H = \frac{1}{2}\}$ (depending on planning limit) and tree depth $H=5$ or $H=4$ (depending on planning limit). Note that the $\epsilon_d$-sequence and $H$ chosen are arbitrary and gave sufficiently good results without finetuning. For RandomPOMDP, we the allocate maximum possible budget (simulation time) to POMCP\footnote{We tuned for exploration constant $c\in \{3,10,100\}$.} and PBVI.
\newline

The source code is available online and with supplementary material. The general FMP planner is made available in a pseudo-code style implementation in Python language.

\subsection{Results and discussion}
\begin{table}[!th]
    \centering
    \begin{tabular}{|c|c|c|c|}
 \hline
 RS[n,k] & RS[7,8] & RS[11,11] & RS[15,15]\\
 \hline
 |S| & 12,544 & 247,808 & 7,372,800 \\
 \hline
 |A| & 13 & 16 & 20 \\
 \hline
 FMP & 20.24$\pm$0.85 & 21.49$\pm$0.63 & 15.67$\pm$0.37 \\
 \hline
 POMCP$^*$ & 20.71$\pm$0.21 & 20.01$\pm$0.23 & 15.32$\pm$0.28 \\
 \hline
 DESPOT$^*$ & 20.93$\pm$0.30 & 21.74$\pm$0.30 & 18.64$\pm$0.28 \\
 \hline
    \end{tabular}
    \caption{Total Discounted Reward. $^*$Reported values.}
    \label{table:rocksample}
\end{table}

\begin{table}[!ht]
    \centering
    \begin{tabular}{|c|c|c|c|}
    \hline
    RP[s,sp] & FMP & POMCP & PBVI \\
         \hline
    RP[100,0.9] & 6.02 & 5.22 & 5.82 \\
         \hline
    RP[100,0.6] & 4.38 & 3.99 & 4.65 \\
         \hline
    RP[100,0.3] & 2.42 & 1.57 & 2.40 \\
         \hline
    RP[60,0.9] & 6.04 & 5.90 & 5.01 \\
         \hline
    RP[60,0.6] & 4.65 & 4.94 & 4.23 \\
         \hline
    RP[60,0.3] & 2.42 & 2.67 & 2.27 \\
         \hline
    RP[30,0.9] & 6.36 & 6.37 & 6.26 \\
         \hline
    RP[30,0.6] & 4.12 & 5.33 & 4.26 \\
         \hline
    RP[30,0.3] & 3.40 & 3.53 & 3.74 \\
         \hline
    \end{tabular}
    \caption{Total Discounted Reward.}
    \label{table:RandomPOMDP}
\end{table}

The results are given in Table~\ref{table:rocksample}  and \ref{table:RandomPOMDP} for RockSample and RandomPOMDP\footnote{For RandomPOMDP, the variance was low and we do not report it simply for clarity.} respectively. It is clear that the performance of FMP is competitive in all settings, even with relatively lax parameters. For RockSample, FMP's python implementation is significantly handicapped compared to reported values of POMCP and DESPOT, obtained by parameter tuning and with much faster C++ implementation. We believe FMP further would benefit from such speed optimization. The RandomPOMDP performance gives some hints to FMP behaviour. It is evident that FMP performs best when there are many beliefs to begin with, i.e, it performs best in low sparsity setting. For high sparsity, it performed competitively to the rest (for larger state sizes), and we guess it misses out in other cases simple due to less lookahead (only $H=5$, again due to planning limit and Python language implementation).
\newline

Finally, we contrast FMP with other planners by noting their tuning and implementation complexity. For example, POMCP requires depth of tree, number of particles used for inference, simulation budget, UCB constant, as parameters. DESPOT in addition to all that, requires maintaining upper/lower bounds on belief nodes. Most other solvers apply such heuristics as well.
\section{Conclusion}

We have a presented a simple algorithm for online planning in POMDPs that is highly competitive with the state-of-the-art despite its apparent simplicity. The main idea behind FMP was to design an algorithm where the inference and planning are tightly coupled to minimise the estimation error. In particular, the algorithm performs approximate inference with adaptive belief discretization so as to upper bound the planning error, with the coarseness of the approximation carefully tuned so as not to overly increase the size of the planning tree. This makes the algorithm both theoretically attractive and practical.
\newline

In future work, there is possibility to extend in multiple directions. we would like to extend our method to continuous spaces, which can be done as the proof strategy only needs kernel function over beliefs. For example, kernels for the sufficient statistics of the belief over continuous state space could be explored. The analysis could also be extended to get regret bounds in Reinforcement Learning setting. Both directions may lead to potential solvers that can handle a wide variety of problems.