\section{Introduction}

Numerical design of microalgae production technologies has been for one decade a source of many interesting challenges not only in engineering but also in the area of scientific computing~\cite{Eilers1993,Masojidek2003,Yoo2012,Lamare2018}.
The potential of these emerging photosynthetic organisms finds interests for the cosmetics, pharmaceutical fields, feed, food and - in the longer term - green chemistry and energy applications~\cite{Wijffels2010}. Outdoor production is mainly carried out in open bioreactors with a raceway shape. Algae grow while exposed to solar radiation in such circular basins, where the water is set in motion by a paddle wheel~\cite{Chiaramonti2013} which homogenizes the medium for ensuring an equidistribution of the nutrients
 and guarantees that each cell will have regularly access to light~\cite{Demory2018}.
The algae are periodically harvested, and their concentration is maintained around an optimal value~\cite{Munoz-Tamayo2013,Posten2016}. Light penetration is strongly reduced by the algal biomass, resulting in limited algal concentrations, which stay  below 1\%~\cite{Bernard2015} of the liquid mass. Above this value, the light extinction is so high that a large fraction of the population evolves in the dark and does not grow anymore. At low biomass density, a fraction of the solar light is not used by the algae and the productivity is suboptimal. Theoretical works have determined the optimal biomass for maximising productivity \cite{Masci2010,Grognard2014}. Remark that other possibilities exist to improve the entire photoproduction process in order to optimize its potential. For instance, a non flat topography shape has been observed to improve the growth rate of the algae in a specific case in~\cite{Bernard02970776}. Here, we are going to prove this result, considering different cases. For this we develop a coupled model to describe the growth of algae in a raceway pond, accounting for the light that they receive. More precisely, extending the study of~\cite{Bernard2013}
 by combining the Han photosynthesis equations with an hydrodynamic law based on the shallow water equations.

This approach enables us to formulate an optimization problem where the raceway topography is designed to maximize the productivity. For this problem, we present an adjoint-based optimization scheme which includes the constraints associated to the shallow water regime. On the contrary to a widespread belief, we prove that the flat topography is optimal in a periodic case for productivity in laminar regime. However, non-trivial topographies can be obtained in other contexts, e.g., when an extra mixing strategy is included in the model. Note that in the examples considered in our numerical tests such topographies only slightly improved the biomass production.

The outline of the paper is as follows. In Section 2, we present the biological and hydrodynamic models underlying our coupled model. In Section 3, we describe the optimization problem and a corresponding numerical optimization procedure. Section 4 is devoted to the numerical results obtained with our approach. We conclude in Section 5 with some perspectives opened by this work.

\section{Coupling hydrodynamic and biological models}

Our approach is based on a coupling between the hydrodynamic behavior of the particles and the evolution of the photosystems driven by the light intensity they receive when traveling across the raceway pond.

\subsection{Modeling the photosystems dynamics}

We consider the Han model~\cite{Han2001} which describes the dynamics of the reaction centers. These subunits of the photosynthetic process harvest photons and transfer their energy to the cell to fix CO$_2$.  In this compartmental model, the photosystems can be described by three different states: open and ready to harvest a photon ($A$), closed while processing the absorbed photon energy ($B$), or inhibited if several photons have been absorbed simultaneously ($C$). The relation of these three states are presented in Fig.~\ref{fig:han}.
\begin{figure}[thpb]
\begin{center}
\usetikzlibrary{shapes.geometric}
\tikzset{block/.style = {draw,rectangle,minimum height = 2.5em,minimum width = 2.5em}}
\begin{tikzpicture}
\node at (0,0)[circle,draw,scale=1.5] (A){$A$};
\node at (3,0)[circle,draw,scale=1.5] (B){$B$};
\node at (6,0)[circle,draw,scale=1.5] (C){$C$};
\node at (1.5,0) [block,draw,scale=0.8](sigI){$\sigma I$};
\node at (4.5,0) [block,draw,scale=0.8](kdsigI){$k_d\sigma I$};
\node at (1.5,1.5) [block,draw,scale=0.8](tau1){$\tau^{-1}$};
\node at (4.5,1.5) [block,draw,scale=0.8](kr){$k_r$};
\node at (0,2)(Pho1){Photon $I$};
\node at (3,2)(Pho2){Photon $I$};
\draw (A)--(sigI)
[->](sigI)--(B);
\draw (B)--(kdsigI)
[->](kdsigI)--(C);
\draw (C) edge[bend right=40] (kr)
(B) edge[<-,bend left=40] (kr);
\draw (B) edge[bend right=40] (tau1)
(A) edge[<-,bend left=40] (tau1);
\draw [red,->](Pho1)--(A);
\draw [red,->](Pho2)--(B);
\end{tikzpicture}
\end{center}
\caption{Scheme of the Han model, representing the probability of state transition, as a function of the photon flux density.}
\label{fig:han}
\end{figure}

Their evolution satisfy the following dynamical system
\begin{equation*}
\left\{
\begin{array}{lr}
\dot{A} = -\sigma I A + \frac B{\tau},\\
\dot{B} =  \sigma I A - \frac B{\tau} + k_rC - k_d\sigma I B,\\
\dot{C} = -k_r C + k_d \sigma I B.
\end{array}
\right.
\end{equation*}
Here $A, B$ and $C$ are the relative frequencies of the three possible states with
\begin{equation}\label{eq:abc}
A+B+C=1,
\end{equation}
and $I$ is the photon flux density, a continuous time-varying signal. The other parameters are $\sigma$, that stands for the specific photon absorption, $\tau$ which is the turnover rate, $k_r$ which represents the photosystem repair rate and $k_d$ which is the damage rate. As shown in~\cite{Lamare2018}, one can use~\eqref{eq:abc}, a fast-slow approximation and singular perturbation theory to reduce this system to a single evolution equation:
\begin{equation}\label{eq:ct}
\dot{C} = -\alpha(I) C + \beta(I),
\end{equation}
where
$\alpha(I) = \beta(I) + k_r,\text{with } \beta(I) =  k_d\tau \frac{(\sigma I)^2}{\tau \sigma I+1}.$
The net specific growth rate is obtained by balancing photosynthesis and respiration, which gives
\begin{equation}\label{eq:mu}
\mu(C,I) = -\gamma(I)C + \zeta(I),
\end{equation}
where
$\zeta(I) = \gamma(I)-R,\text{ with } \gamma(I) = \frac{k\sigma I}{\tau \sigma I+1}.$
Here, $k$ is a factor that relates received energy with growth rate. The term $R$ represents the respiration rate.

In this framework, the dynamics of the biomass $X$ in an open system where a fraction of the biomass is continuously sampled and is derived from $\mu$:
\begin{equation}\label{eq:Xt}
\dot{X} = \mu(C,I)X-D X,
\end{equation}
where $D$ is the dilution rate. The growth rate~\eqref{eq:mu} associated with the steady state, for a constant $I$ is given by
\begin{equation}\label{eq:muss}
\mu(I) = -\gamma(I)C^*(I) + \zeta(I),
\end{equation}
where
\begin{equation}\label{eq:css}
C^*(I) =\frac{\beta(I)}{\alpha(I)}.
\end{equation}

\subsection{Steady 1D shallow water equations}\label{sec:model_hydro}

The shallow water equations are one of the most popular model for describing geophysical flows, which is derived from the free surface incompressible Navier-Stokes equations (see for instance~\cite{Gerbeau2001}). In the current study, we focus on the smooth steady state solutions of the shallow water equations in a laminar regime. Such steady states are governed by the following partial differential equations:
\begin{align}
&\partial_x(h u) = 0, \label{eq:sv1} \\
&\partial_x(h u^2+g\frac{h^2}2) = -g h \partial_x z_b, \label{eq:sv2}
\end{align}
where $h$ is the water elevation, $u$ is the horizontal averaged velocity of the water, the constant $g$ stands for the gravitational acceleration, and $z_b$ defines the topography. The free surface $\eta$ is given by $\eta = h + z_b$ and the averaged discharge $Q=hu$. This system is presented in Fig.~\ref{fig:shallow}.
\begin{figure}[thpb]
\begin{center}
\begin{tikzpicture}
\draw [thick](0,-4) -- (6,-4);
\draw (0,-4) node[anchor=north] {0};
\draw (6,-4) node[anchor=north] {$L$};
\draw[thick, ->] (0,-4) -- (0,1) node[anchor=east] {$z$};
\draw (0,0) node[anchor=east] {$0$};
\draw[dashed,->] (0,0) -- (6,0) node[anchor=west]  {$x$};
\draw [dashed, <-](1.5,-0.48) -- (1.5,0);
\draw (1.5,0.3) node{$\eta(x)$};
\draw[thick,red,<-] (2.5,-0.25) -- (2.5,1);
\draw[thick,red,<-] (3,0) -- (3,1);
\draw[thick,red,<-] (3.5,0.2) -- (3.5,1);
\draw[thick,red,<-] (4,0.35) -- (4,1);
\draw[thick,red,<-] (4.5,0.4) -- (4.5,1);
\draw[thick,red,<-] (5,0.3) -- (5,1);
\draw[thick,red,<-] (5.5,0.1) -- (5.5,1);
\draw (4,1.5) node {$I_s$};
\draw (0,0.2) sin (1.5,-0.5) cos (3,0) sin (4.5,0.4) cos (6,-0.2);
\draw[fill, gray] (0,-4) -- (0,-3.6) sin (1.5,-3.2) cos (3,-3.4) sin (4.5,-3.6) cos (6,-3.2) -- (6,-4);
\draw [dashed,->](5.3,0) -- (5.3,-3.5);
\draw (5.4,-1.8) node[anchor=west] {$z_b(x)$};
\draw[dashed,<-] (3.5,0.15) -- (3.5,-3.5);
\draw (3.6,-1.8) node[anchor=west]{$h(x)$};
\draw [->](1.3,-2) -- (1.6,-2) ;
\node at (2,-1.8){$u(x)$};
\draw [->](1.3,-1.5) -- (1.6,-1.5);
\draw [->](1.3,-2.5) -- (1.6,-2.5);
\draw [->](1.3,-1) -- (1.6,-1);
\draw [->](1.3,-3) -- (1.6,-3);
\end{tikzpicture}
\end{center}
\caption{Representation of the hydrodynamic model.}
\label{fig:shallow}
\end{figure}
The $z$ axis represents the vertical direction and the $x$ axis represents the horizontal direction. Besides, $I_s$ represents the light intensity at the free surface (assumed to be constant).

Integrating~\eqref{eq:sv1}, we get
\begin{equation}\label{eq:Q0}
h u=Q_0,
\end{equation}
for a fixed positive constant $Q_0$, which implies a constant discharge in space. Then~\eqref{eq:sv2} can be rewritten as
\begin{equation}\label{eq:sv2mod}
h u\partial_x u + h\partial_x g h + h\partial_x g z_b = 0.
\end{equation}
Let us assume that $h$ is strictly positive. Dividing~\eqref{eq:sv2mod} by $h$ and using~\eqref{eq:Q0} to eliminate $u$, we get
\begin{equation*}
\partial_x \Big( \frac{Q_0^2}{2h^2} + g (h + z_b)\Big) = 0.
\end{equation*}
This equation corresponds to the Bernoulli's principle. Now let us consider two fixed constants $h(0), z_b(0)\in \R$. For all $x\in [0,L]$, we obtain
\begin{equation*}
\frac{Q_0^2}{2h^2} + g (h + z_b) = \frac{Q_0^2}{2h^2(0)} + g (h(0) + z_b(0)) =:M_0,
\end{equation*}
meaning that the topography $z_b$ satisfies
\begin{equation}\label{eq:zb}
z_b = \frac{M_0}g - \frac{Q_0^2}{2g h^2} - h.
\end{equation}
\begin{rem}\label{rem:flu}
Define the Froude number for the steady state by $Fr=u/\sqrt{g h}$. The situation $Fr<1$ corresponds to the subcritical case (i.e. the flow regime is fluvial) while $Fr>1$ is to the supercritical case (i.e. the flow regime is torrential). In particular, the threshold value of $h$ for $Fr=1$ is given by
\begin{equation*}
h_c := (\frac{Q^2_0}{g})^{\frac 13}.
\end{equation*}
\end{rem}

Because of~\eqref{eq:zb}, $h$ is the solution of a third order polynomial equation. Given a smooth topography $z_b$, if
\begin{equation*}
h_c + z_b + \frac{Q_0^2}{2g h^2_c} - \frac{M_0}g < 0.
\end{equation*}
there exists a unique positive smooth solution of~\eqref{eq:zb} which satisfies the subcritical flow condition (see~\cite[Lemma 1]{Michel-Dansac2016}). The velocity $u$ can then be obtained thanks to~\eqref{eq:Q0}.

\subsection{Lagrangian trajectories of the algae and captured light intensity}

Let $z(t)$ be the depth of a particle at time $t$ in the raceway pond. We first determine the Lagrangian trajectory of an algal cell that starts at a given position $z(0)$ at time 0.

From the incompressibility of the flow, we have $\nabla \cdot \underline{\textbf{u}} = 0$ with $\underline{\textbf{u}} = (u(x),w(x,z))$. Here,
$w(x,z)$ is the vertical velocity. This implies that
\begin{equation}\label{eq:incompress}
\partial_x u + \partial_z w = 0.
\end{equation}
Integrating~\eqref{eq:incompress} from $z_b$ to $z$ gives:
\begin{equation*}
\begin{split}
0 &= \int_{z_b}^z \big(\partial_x u(x) + \partial_z w(x, z) \big)\D z \\
&= \partial_x \int_{z_b}^z  u(x) \D z +\int_{z_b}^z \partial_z w(x, z) \D z\\
& =  \partial_x \big((z-z_b)u(x)\big) + w(x,z) - w(x, z_b) \\
& = (z - z_b) \partial_x u(x) - u(x)\partial_x z_b + w(x, z),
\end{split}
\end{equation*}
where we have used the kinematic condition at the bottom (i.e. $w(x,z_b) = u(x)\partial_x z_b$). It then follows from~\eqref{eq:zb} that
\begin{equation*}
w(x , z) = (\frac{M_0}g - \frac{3u^2(x)}{2g} - z)u'(x).
\end{equation*}
The Lagrangian trajectory is characterized by the system
\begin{equation}\label{eq:xtzt}
\begin{pmatrix}
\dot{x}(t)\\
\dot{z}(t)
\end{pmatrix}
=
\begin{pmatrix}
u(x(t))\\
w(x(t),z(t))
\end{pmatrix}.
\end{equation}
\begin{rem}
Note that a Brownian motion can be included in these equations to take into account the fluid viscosity. In this case, the optimisation procedure required a large set of simulations and an averaging strategy. Such questions are out of the range of this paper.
\end{rem}
Up to now  and for the sake of simplicity, we denote by $z(x)$ the depth of a particle at $x$. From~\eqref{eq:xtzt}, we get
\begin{equation}\label{eq:z'x}
z' := \frac{\dot{z}}{\dot{x}} = (\frac {M_0}g -\frac{3u^2}{2g} - z)\frac{u'}{u}.
\end{equation}
Note that from~\eqref{eq:Q0} and~\eqref{eq:zb}, we have $\eta=h+z_b = \frac{M_0}g - \frac{u^2}{2g}$, which implies $\eta'=-u u'/g$. Multiplying~\eqref{eq:z'x} both sides by $u$ and using the formulation of $\eta, \eta'$ gives
\begin{equation*}
\begin{split}
z' u + z u'&= (\eta - \frac{u^2}g)u'\\
&=\eta u' + \eta' u,
\end{split}
\end{equation*}
which implies that $\big(u(z-\eta)\big)'=0$. We then obtain

\begin{equation}\label{eq:zx}
z(x) = \eta(x) + \frac{u(0)}{u(x)}(z(0) - \eta(0)).
\end{equation}
The computation of trajectories in the shallow water system can be carried out with this formula, where $u(0)$, $z(0)$ and $\eta(0)=h(0)+z_b(0)$ are given as boundary conditions, whereas $u(x)$ can be obtained by~\eqref{eq:Q0} and~\eqref{eq:zb}. On the contrary, $h(x)$, hence $u(x)$ will be considered  in section~\ref{sec:opti} as a parameter to be optimized, and~\eqref{eq:zb} will be used to (post-)compute $z_b(x)$.

\begin{rem}
Since $Q_0$ is chosen to be positive, $h$ is necessarily positive and so does $u$ from~\eqref{eq:Q0}. Moreover, if $z(0)$ belongs to  $[z_b(0),\eta(0)]$, then $z(x)$ belongs to $[z_b(x),\eta(x)]$. In particular, choosing $z(0)=z_b(0)$ in~\eqref{eq:zx} and using~\eqref{eq:Q0} gives $z(x) = z_b(x)$. In the same way, we find that $z(x)=\eta(x)$ when $z(0) = \eta(0)$.
\end{rem}

To evaluate the light intensity observed on the trajectory $z$, we assume that the growth process occurs at a much slower time scale than those of hydrodynamics and is, as such, negligible for one lap over the raceway. As a consequence, the turbidity is supposed to be constant over the considered time scale. In this framework, the Beer-Lambert law describes how light is attenuated with depth:
\begin{equation}\label{eq:Beer}
I(x,z) = I_s\exp\Big(-\varepsilon (\eta(x) - z)\Big).
\end{equation}
Here $\varepsilon$ is the light extinction coefficient. We assume that the system is perfectly mixed so that the concentration of the biomass $X$ in~\eqref{eq:Xt} is homogeneous and $\varepsilon$ is constant. Combining~\eqref{eq:Beer} with~\eqref{eq:zx}, we get the following  expression for the captured light intensity along the trajectory $z$
\begin{equation}\label{eq:Ixz0}
I(x,z(x)) =I_s \exp\Big(-\varepsilon  \frac{u(0)}{u(x)} (\eta(0) -z(0) ) \Big).
\end{equation}
It follows that in this approach, the light intensity $I$ couples the hydrodynamic model and Han model: the trajectories of the algae define the received light intensity, which is used in the photosystems dynamics.

We can then derive the equation satisfied by $C$. Indeed, repeating the reasoning done to get~\eqref{eq:z'x} with \eqref{eq:ct}, we find a time-free reformulation, namely
\begin{equation}\label{eq:c'x}
C':=\frac{\dot{C}}{\dot{x}} = -\frac{\alpha(I)}{u}C +  \frac{\beta(I)}{u},
\end{equation}
where all the functions on the right-hand side only depend on $x$.

\section{Optimization problem}\label{sec:opti}

In this section, we define the optimization problem associated with our biological-hydrodynamic model. We first introduce our procedure for a simple problem, and then extend this to other variant cases.

\subsection{Optimization functional}

The average net specific growth rate over the domain is defined by
\begin{equation}\label{eq:mucontinous}
\bar{\mu} := \frac 1L\int_0^L\frac 1{h(x)}\int_{z_b(x)}^{\eta(x)} \mu\big(C(x, z), I(x, z)\big) \D z \D x.
\end{equation}
In order to tackle numerically this optimization problem, let us consider a vertical discretization. Let $N_z$ denotes the number of particles, and consider a uniform vertical discretization of their initial position:
\begin{equation*}
z_i(0) = \eta(0) - \frac {i-\frac12}{N_z}h(0), \quad i=1,\ldots,N_z.
\end{equation*}
From~\eqref{eq:zx}, we obtain
\begin{equation*}
z_{i}(x) - z_{i+1}(x) = \frac 1{N_z}h(x), \quad i=1,\ldots,N_z,
\end{equation*}
meaning that the distribution of particles remains uniform along the trajectories. To simplify notations, we write $I_i(x)$ instead of $I(x,z_i(x))$ hereafter.

Let $C_i(x)$ (resp. $I_i(x)$) the photo-inhibition state (resp. the light intensity) associated with the trajectories $z_i(x)$. Then the semi-discrete average net specific growth rate in the raceway pond can be defined from~\eqref{eq:mucontinous} by
\begin{equation}\label{eq:mudiscret}
\begin{split}
\bar \mu_{\Delta} :&= \frac 1L \int_0^L \frac 1{h(x)}\frac 1{N_z}h(x)\sum_{i=1}^{N_z}\mu( C_i(x),I_i(x) ) \D x\\
&= \frac 1{L N_z}\sum_{i=1}^{N_z} \int_0^L\mu( C_i(x),I_i(x) ) \D x.
\end{split}
\end{equation}

From now on, we focus only on the subcritical case, i.e. $Fr<1$, see Remark~\ref{rem:flu}. As we have mentioned in the previous section, in this regime, a given topography $z_b$ corresponds to a unique water height $h$ which verifies this assumption. Note that the volume of our system is simply given by
\begin{equation}\label{eq:vol}
V = \int_0^L h(x) \D x.
\end{equation}
Therefore, we choose to parameterize $h$ by a vector $a:=[a_1,\cdots,a_N]\in \R^{N}$ in order to handle the volume of our system. This vector will be the variable to be optimized. Once we find the optimum parameter $a^*$, we can determinate the associated topography by means of~\eqref{eq:zb}.
\begin{rem}
In usual shallow water solver, equations of type~\eqref{eq:zb}  are usually consider to compute $h$ in the simulations. Here, we use this equation in the opposite way, i.e., to recover the topography $z_b$ from $h$.
\end{rem}

Our goal is to optimize the topography to maximize $\bar \mu_{\Delta}$. Up to now, we denote explicitly that the functions depend on $a$.
In the case of one trajectory, we consider the functional
\begin{equation*}
J(a) = \frac1L\int_0^L \frac{-\gamma(I(x;a))C(x) + \zeta(I(x;a))}{u(x;a)} \D x ,
\end{equation*}
where $C$ satisfies the following parameterized version of~\eqref{eq:c'x}
\begin{equation}\label{eq:cond}
C' + \frac{\alpha(I(x;a))}{u(x;a)}C =  \frac{\beta(I(x;a))}{u(x;a)}.
\end{equation}
The optimization problem then reads:

\textit{Find $a^*$ solving the maximization problem:}
\begin{equation}\label{eq:Jopt}
\max_{a\in \R^{N}} J(a).
\end{equation}

\subsection{Optimality System for maximizing growth rate}

In this section, we introduce the optimality system derived from the optimization problem~\eqref{eq:Jopt} and some of its variants. For the sake of simplicity, we first remain in the setting where one single trajectory is considered, and then extend it to a multiple trajectories system. Next, we analyse the case where $C$ is periodic, and finally, we present a more general case where the volume can also be optimized. For simplicity of notation, we omit $x$ in the notation.

\subsubsection{One trajectory problem}

Define the Lagrangian of Problem~\eqref{eq:Jopt} by
\begin{equation*}
\begin{split}
\mathcal{L}(C,p,a) = \frac1L&\int_0^L \frac{-\gamma(I(a))C +\zeta(I(a))}{u(a)} \D x\\
- &\int_0^L p\big(C' + \frac{\alpha(I(a))C-\beta(I(a))}{u(a)}\big) \D x,
\end{split}
\end{equation*}
where $p$ is the Lagrange multiplier associated with the constraint~\eqref{eq:cond}.

The optimality system is obtained by cancelling all the partial derivatives of $\mathcal{L}$. Differentiating $\mathcal{L}$ with respect to $p$ and equating the resulting expression to zero gives~\eqref{eq:cond}. Integrating the terms $\int p C'\D x$ on the interval $[0,L]$ by parts and differentiating $\L$ with respect to $C$ and $C(L)$ gives rise to
\begin{equation}\label{eq:adjoint}
\left\{
\begin{array}{lr}
\partial_{C}\mathcal{L} = p' - p \frac{\alpha(I(a))}{u(a)} - \frac1L\frac{\gamma(I(a))}{u(a)} \\
\partial_{C(L)}\mathcal{L}  = p(L).\\
\end{array}
\right.
\end{equation}
Given a vector $a$, let us still denote by $C, p$ the corresponding solutions of~\eqref{eq:cond} and~\eqref{eq:adjoint}. The gradient $\nabla J(a)$ is obtained by
\begin{equation*}
\nabla J(a) = \partial_{a}\mathcal{L},
\end{equation*}
where
\begin{equation}\label{eq:grad1}
\begin{split}
\partial_{a}\mathcal{L} = \frac1L&\int_0^L \frac{-\gamma'(I(a))C+\zeta'(I(a))}{u(a)}\partial_a I(a)\D x \\
-\frac1L&\int_0^L \frac{-\gamma (I(a))C + \zeta(I(a))}{u^2(a)}\partial_a u(a)\D x \\
+&\int_0^L p\frac{-\alpha'(I(a))C + \beta'(I(a))}{u(a)}\partial_a I(a)\D x \\
-&\int_0^L p\frac{-\alpha(I(a))C + \beta(I(a))}{u^2(a)}\partial_a u(a) \D x.
\end{split}
\end{equation}

\subsubsection{Multiple trajectories problem}

We now extend the previous framework to deal with multiple trajectories. From~\eqref{eq:mudiscret}, the objective functional is defined by:
\begin{equation}\label{eq:muNza}
\bar \mu_{\Delta}(a) = \frac 1{L N_z}\sum_{i=1}^{N_z}\int_0^L\frac{-\gamma(I_i(a))C_i + \zeta(I_i(a))}{u(a)}   \D x,
\end{equation}
where $C_i$ satisfies~\eqref{eq:cond} for the light intensity $I_i$ with $i=1,\ldots,N_z$. Let us still denote by $\L$ the associated Lagrangian. Similar computations as in the previous section give rise to $N_z$ systems similar to~\eqref{eq:adjoint}, where $C$ is replaced by $C_i$. Denoting by $p_{i}$ the associated Lagrange multipliers, the partial derivative $\partial_a\mathcal{L}$ (hence the gradient $\nabla \bar \mu_{\Delta}(a)$) is then given by
\begin{equation*}
\begin{split}
\partial_{a}\mathcal{L} = \frac 1{L N_z}&\sum_{i=1}^{N_z}\int_0^L \frac{-\gamma'(I_i(a))C_i+\zeta'(I_i(a))}{u(a)}\partial_a I_i(a)\D x \\
-\frac 1{L N_z}&\sum_{i=1}^{N_z}\int_0^L \frac{-\gamma (I_i(a))C_i + \zeta(I_i(a))}{u^2(a)}\partial_a u(a)\D x \\
+&\sum_{i=1}^{N_z}\int_0^L p_{i}\frac{-\alpha'(I_i(a))C_i + \beta'(I_i(a))}{u(a)}\partial_a I_i(a)\D x \\
-&\sum_{i=1}^{N_z}\int_0^L p_{i}\frac{-\alpha(I_i(a))C_i + \beta(I_i(a))}{u^2(a)}\partial_a u(a) \D x.
\end{split}
\end{equation*}

\subsubsection{Periodic problem}

We now consider a variant of our problem, where the photo-inhibition state $C$ is periodic, with a period corresponding to a raceway lap. This situation
occurs, e.g., when an appropriate harvest is performed after each lap. To describe the corresponding model, let us first consider a variant of the usual Cauchy problem for~\eqref{eq:c'x}:

\textit{Given $I\in L^2(0,L;R)$, find $(C_0,C)\in \R\times \C(0,L;\R)$ such that}
\begin{equation}\label{eq:periodiceq}
\left\{
\begin{array}{lr}
C'(x) = \frac{-\alpha(I(x))C(x) + \beta(I(x))}{u(x)},\quad x\in[0,L]&\\
C(L) = C(0)= C_0.&
\end{array}
\right.
\end{equation}

The following theorem gives the existence and uniqueness of a (weak) solution of the problem~\eqref{eq:periodiceq}
\begin{thm}\label{thm:perC}
Given $I\in L^2(0,L;R)$, there exists a unique couple $(C_0,C) \in \R\times \C(0,L;\R)$ satisfying
\begin{equation}\label{periodicC}
\left\{
\begin{array}{lr}
C(x) = C_0 + \int_0^x\frac{-\alpha(I(s)) C(s) + \beta(I(s))}{u(s)}\D s&,\\
C(L)  = C_0&
\end{array}
\right.
\end{equation}
for all $x\in [0,L]$.
\end{thm}

The proof is given in Appendix A. The optimality system is obtained as in the previous sections.

\begin{rem}\label{rem:pper}
Note that the periodicity of $C$ implies that $p$ is also periodic. Indeed, since $C(L)=C(0)$, then by differentiating $\mathcal{L}$ with respect to $C(L)$, we have
\begin{equation*}
\partial_{C(L)} \mathcal{L} = p(L) - p(0).
\end{equation*}
so that equating the above equation to zero gives the periodicity for $p$. Similar computations apply for multiple trajectories problem.
\end{rem}

Next, we state a result for Problem~\eqref{eq:Jopt} in the case where $C$ is periodic.

\begin{prop}\label{prop:flattopo}
Assume $h$ is of the form
\begin{equation}\label{eq:Fourier}
h(x;a) = a_0 + \sum_{n=1}^N a_n \sin(2n\pi \frac x L),
\end{equation}
with $a:=[a_1,\ldots,a_N]$ chosen such that the subcritical condition is satisfied (see Remark~\ref{rem:flu}). Then $\nabla J(0)=0$.
\end{prop}
As a consequence, the flat topography is critical point for Problem~\eqref{eq:Jopt}. Note that the larger $N$ is, the less valid is our hydrodynamic model valid, see Section~\ref{sec:model_hydro}, where a smooth topography is assumed. Hence, limit situations where $N\rightarrow+\infty$ are not considered in what follows.

\begin{pf}
Let $q \in[0,1]$, consider an initial position by
\begin{equation}\label{eq:z0}
z(0;a) =  \eta(0;a) - q h(0;a).
\end{equation}
Replacing $z(0;a)$ by the previous expression in~\eqref{eq:Ixz0} and using~\eqref{eq:Q0}, gives
\begin{equation}\label{eq:Ia}
I(a) = I_s \exp(-\varepsilon q\frac{Q_0}{u(a)}) = I_s \exp(-\varepsilon q h(a)).
\end{equation}
Introduce $a^* = [0, \ldots, 0]$. Replacing $a$ by $a^*$ in~\eqref{eq:Q0} and~\eqref{eq:Ia}, we find that $u(a^*)$ and $I(a^*)$ are both constants with respect to $x$. For a given $C(0)$,~\eqref{eq:cond} gives
\begin{equation}\label{eq:Cx}
C(x) = e^{-\frac{\alpha(I(a^*))}{u(a^*)}x}C(0) + \frac{\beta(I(a^*))}{\alpha(I(a^*))}(1 - e^{-\frac{\alpha(I(a^*))}{u(a^*)}x}).
\end{equation}
Since $C$ is periodic (i.e. $C(L)=C(0)$), we get from the previous equation that $ C(0) = \frac{\beta(I(a^*))}{\alpha(I(a^*))}$. Inserting this value in~\eqref{eq:Cx}, we find
\begin{equation*}
C(x) = \frac{\beta(I(a^*))}{\alpha(I(a^*))}, \quad \forall x\in[0,L],
\end{equation*}
which corresponds to the steady state~\eqref{eq:css}, that we denote by $C^*$ hereafter. We then compute all the partial derivatives appearing in~\eqref{eq:grad1}. We get:
\begin{equation*}
\begin{split}
&\partial_{a_n} u(a) = -\frac{Q_0}{h^2(a)}\partial_{a_n} h(a),\\
&\partial_{a_n} I(a) = -\varepsilon q I(a)\partial_{a_n} h(a),
\end{split}
\end{equation*}
but $\partial_{a_n} h(a) = \sin(2n\pi \frac x L)$, so that the integral from $0$ to $L$ of $\partial_{a_n} h(a)$ is zero. Replacing $a$ by $a^*$ in~\eqref{eq:grad1} , the first integral can be computed as
\begin{equation*}
\begin{split}
&\int_0^L \frac {-\gamma'(I(a^*))C^* + \zeta'(I(a^*))}{u(a^*)} \partial_{a_n} I(a^*) \D x\\
= &\frac {-\gamma'(I(a^*))C^* + \zeta'(I(a^*))}{u(a^*)} \int_0^L\partial_{a_n} I(a^*) \D x \\
= &0.
\end{split}
\end{equation*}
and the second integral can be computed as
\begin{equation*}
\begin{split}
 &\int_0^L\frac {-\gamma(I(a^*))C^*+\zeta(I(a^*))}{u^2(a^*)} \partial_{a_n} u(a^*)\D x \\
=&\frac {-\gamma(I(a^*))C^*+\zeta(I(a^*))}{u^2(a^*)} \int_0^L\partial_{a_n} u(a^*)\D x \\
= &0.
\end{split}
\end{equation*}
Since $p$ is also periodic (see Remark~\ref{rem:pper}), it can be proved in much the same way that the two last terms in~\eqref{eq:grad1} also cancels when $a=a^*$. The result follows.
\end{pf}

\begin{rem}
Numerically, we observe that the flat topography is actually optimal in the periodic case.
\end{rem}

\begin{rem}
Note that if $C$ is not periodic,~\eqref{eq:Cx} implies that $C$ depends also on the space variable $x$ and the computations for the first-two integrals in the proof above will no longer hold. For the same reason, $p$ shall be not constant and the last two integrals in~\eqref{eq:grad1}, hence the gradient will not cancel in general. In other words, the flat topography is not an optimum for the case $C$ non periodic, which is confirmed by our numerical tests  (see Section~\ref{sec:testC0}).
\end{rem}

\subsubsection{Non-constant volume problem for maximizing areal productivity}\label{sec:nonconstant}

In general, the extinction coefficient is an affine function with respect to the biomass $X$~\cite{Martinez201811}
\begin{equation}\label{eq:eps}
\varepsilon(X) = \alpha_0 X + \alpha_1,
\end{equation}
where $\alpha_0>0$ the specific light extinction coefficient of the microalgae specie and $\alpha_1$ the background turbidity that summarizes the light absorption and diffusion due to all non-microalgae components. Recall that according to~\eqref{eq:Xt}, the biomass concentration $X$ depends on the growth rate $\mu$ and on the harvesting of the raceway $D$.

\begin{rem}\label{rem:optXcond}
A standard criterion to determine a relevant value of $X$~\cite{Masci2010,Grognard2014}, and then an harvesting frequency, consists in regulating it to a value such that the steady-state value of the growth rate $\mu$ at the (average) bottom depth $\bar z_b$ is 0, i.e.,
\begin{equation}\label{eq:muIzb}
\mu(I_{\bar z_b})=0,
\end{equation}
where $\mu(I)$ is given by~\eqref{eq:muss}. Solving~\eqref{eq:muIzb} provides a value of $I_{\bar z_b}$ thus of the extinction coefficient $\varepsilon(X)$, and finally of the biomass concentration $X$. This assumption is usually called the \textit{compensation condition}, for which the growth at the bottom compensates exactly the respiration.
\end{rem}

In this framework, maximizing areal productivity is a relevant target. Productivity per unit of surface for a given biomass concentration $X$ is given by:
\begin{equation}\label{eq:Pi}
\Pi = \bar \mu X \frac{V}{S},
\end{equation}
where $S$ is the ground surface of the raceway pond and $\bar \mu$ is defined in~\eqref{eq:mucontinous}. To optimize this objective function, a more general setting where the reactor volume varies can be considered. The total biomass $X$ in a given volume is supposed to satisfy~\eqref{eq:muIzb} and we assume that it is maintained constant to achieve the
compensation conditions, see Remark~\ref{rem:optXcond}.

Let us keep the parameterization $h$ given in~\eqref{eq:Fourier}. The volume of the system satisfies $V=a_0L$ (see~\eqref{eq:vol}). Moreover, since we consider a 1D framework, we have $S=L$. The biomass concentration $X$ is chosen to be solution of~\eqref{eq:muIzb}. We have
\begin{equation*}
X(a_0)\frac{V(a_0)}{S}= \alpha_2 - \alpha_3 a_0,
\end{equation*}
where $\alpha_2=\frac{1}{\alpha_0}\ln(\frac{I_s}{I_{\bar z_b}})$ and $\alpha_3=\frac{\alpha_1}{\alpha_0}$, and $\alpha_0, \alpha_1$ are given in~\eqref{eq:eps}. The detail of the computation is given in Appendix C. Let us then consider the extended parameter vector $\tilde a:=[a_0,a]$. For simplicity, we again consider one single trajectory and optimize the objective function associated with~\eqref{eq:Pi} defined by
\begin{equation*}
\tilde J(\tilde a) = \frac{\alpha_2-\alpha_3a_0}{L}\int_0^L \frac{-\gamma(I(\tilde a))C +\zeta(I(\tilde a)) }{u(\tilde a)} \D x,
\end{equation*}
where $C$ satisfies~\eqref{eq:cond} with $a$ replacing by $\tilde a$.
The corresponding optimization problem reads:

\textit{Find $\tilde a^*$ solving the maximization problem:}
\begin{equation}\label{eq:Joptilde}
\max_{\tilde a\in \R^{N+1}} \tilde J(\tilde a).
\end{equation}
The Lagrangian of Problem~\eqref{eq:Joptilde} is
\begin{equation*}
\begin{split}
\tilde{\mathcal{L}}(C,p,\tilde a) = &\frac{\alpha_2-\alpha_3a_0}{L}\int_0^L \frac{-\gamma(I(\tilde a))C+\zeta(I( \tilde a))}{u(\tilde a)} \D x\\
- &\int_0^L p\big(C' + \frac{\alpha(I(\tilde a))C-\beta(I(\tilde a))}{u(\tilde a)}\big) \D x,
\end{split}
\end{equation*}
where $p$ is the Lagrangian multiplier associated with the constraint~\eqref{eq:cond} for $C$. By computations similar to that of the previous section, we find as optimality system
\begin{equation*}
\left\{
\begin{array}{lr}
\partial_{C} \tilde{\mathcal{L}} = p' - p \frac{\alpha(I(\tilde a))}{u(\tilde a)} -\frac{\gamma(I(\tilde a))}{u(\tilde a)}\frac{\alpha_2-\alpha_3a_0}{L} \\
\partial_{C(L)}\tilde{\mathcal{L}}  = p(L).
\end{array}
\right.
\end{equation*}
However, there is an extra element in the gradient $\tilde J(\tilde a)$, namely the derivative with respect to $a_0$, so that $\nabla \tilde J(\tilde a): = [\partial_{a_0}\tilde{\mathcal{L}},  \partial_{a}\tilde{\mathcal{L}}]$, where
\begin{equation}\label{eq:grad1tilde1}
\begin{split}
\partial_{a_0}\tilde{\mathcal{L}} =  \frac{\alpha_2-\alpha_3a_0}{L}&\int_0^L \frac{-\gamma'(I(\tilde a))C+\zeta'(I(\tilde a))}{u(\tilde a)}\partial_{a_0} I(\tilde a)\D x\\
-\frac{\alpha_2-\alpha_3a_0}{L}&\int_0^L \frac{-\gamma (I(\tilde a))C + \zeta(I(\tilde a))}{u^2(\tilde a)}\partial_{a_0} u(\tilde a)\D x \\
-\frac{\alpha_3}{L}&\int_0^L \frac{-\gamma (I(\tilde a))C + \zeta(I(\tilde a))}{u(\tilde a)}\D x \\
+&\int_0^L p\frac{-\alpha'(I(\tilde a))C + \beta'(I(\tilde a))}{u(\tilde a)}\partial_{a_0} I(\tilde a)\D x \\
-&\int_0^L p\frac{-\alpha(I(\tilde a))C + \beta(I(\tilde a))}{u^2(\tilde a)}\partial_{a_0} u(\tilde a) \D x,
\end{split}
\end{equation}
\begin{equation}\label{eq:grad1tilde2}
\begin{split}
\partial_{a}\tilde{\mathcal{L}} =  \frac{\alpha_2-\alpha_3a_0}{L}&\int_0^L \frac{-\gamma'(I(\tilde a))C+\zeta'(I(\tilde a))}{u(\tilde a)}\partial_a I(\tilde a)\D x \\
-\frac{\alpha_2-\alpha_3a_0}{L}&\int_0^L \frac{-\gamma (I(\tilde a))C + \zeta(I(\tilde a))}{u^2(\tilde a)}\partial_a u(\tilde a)\D x \\
+&\int_0^L p\frac{-\alpha'(I(\tilde a))C + \beta'(I(\tilde a))}{u(\tilde a)}\partial_a I(\tilde a)\D x \\
-&\int_0^L p\frac{-\alpha(I(\tilde a))C + \beta(I(\tilde a))}{u^2(\tilde a)}\partial_a u(\tilde a) \D x.
\end{split}
\end{equation}
However, we still find that the flat topography is optimal in that case.
\begin{prop}
Assume $h$ is of the form~\eqref{eq:Fourier} and $C$ is periodic, then $\tilde J(\tilde a^*)=0$ with $\tilde a^*=[a^*_0,0,\cdots,0]$.
\end{prop}

\begin{pf}
The proof is similar to that of Proposition~\ref{prop:flattopo}. For a given $z(0;\tilde a)$ defined in~\eqref{eq:z0}, using the result in Appendix C to replace $\varepsilon$ in~\eqref{eq:Ia}, we have
\begin{equation*}
I(\tilde a) = I_s \exp\Big(- \frac 1{a_0}\ln(\frac{I_s}{I_{\bar z_b}})q h(\tilde a)\Big).
\end{equation*}
The two extra partial derivatives appearing in~\eqref{eq:grad1tilde1} are computed by
\begin{equation*}
\begin{split}
&\partial_{a_0} u(\tilde a) = -\frac{Q_0}{h^2(\tilde a)}\partial_{a_0} h(\tilde a),\\
&\partial_{a_0} I(\tilde a) = -( -\frac 1{a^2_0}h(\tilde a) + \frac 1{a_0}\partial_{a_0} h(\tilde a) )\ln(\frac{I_s}{I_{\bar z_b}})q I(\tilde a),
\end{split}
\end{equation*}
where $\partial_{a_0} h(\tilde a)=1$. The flat topography corresponds to the parameter $\tilde a^* = [a^*_0,0,\ldots,0]$. Then $\partial_{a_0} u(\tilde a^*)=-\frac{Q_0}{{a^*_0}^2}$ and $\partial_{a_0} I(\tilde a^*)=0$. Hence the first integral and the fourth integral in $\partial_{a_0}\tilde{\mathcal{L}}$ in~\eqref{eq:grad1tilde1} are zero. Following a reasoning similar to that of Proposition~\ref{prop:flattopo}, we find $C^*$ is constant. Hence, the last integral in $\partial_{a_0}\tilde{\mathcal{L}}$ is zero. Finally for
\begin{equation}\label{eq:a0}
a^*_0=\frac{\alpha_2}{2\alpha_3},
\end{equation}
the sum of the second and the third integral in $\partial_{a_0}\tilde{\mathcal{L}}$ is zero.

As for~\eqref{eq:grad1tilde2}, it follows from computations similar to that of the proof of Proposition~\ref{prop:flattopo} that this term equals zero, which completes the proof.
\end{pf}

\section{Numerical Experiments}

In this section, we show some optimal topographies obtained in the various previous frameworks.

\subsection{Numerical Methods}

We start with an algorithm to solve the optimization problem and the numerical solver that we use for our numerical tests.

\subsubsection{Gradient Algorithm}

We detail now a gradient-based optimization procedure and apply it to Problem~\eqref{eq:Jopt}. The procedure is given in Algorithm~\ref{algo}.
\begin{algorithm}
\begin{algorithmic}[1]
\State \textbf{Input}: Tol$>0$, $\rho>0$.
\State \textbf{Initial guess}: $a$.
\State \textbf{Output}: $a$
\State Set $err :=$ Tol$+1$ and define $h$ by~\eqref{eq:Fourier} using the input data.
\While{$err>$Tol and $\norm{h}_{\infty} > h_{c}$}
\State Compute $u$ by~\eqref{eq:Q0} and $I$ by~\eqref{eq:Ixz0}.
\State Set $C$ as the solution of~\eqref{eq:cond}.
\State Set $p$ as the solution of~\eqref{eq:adjoint}.
\State Compute the gradient $\nabla J$ by~\eqref{eq:grad1}.
\State $a = a + \rho \nabla J$,
\State Set $err := \| \nabla J\|$.
\EndWhile
\end{algorithmic}
\caption{Gradient-based optimization algorithm}
\label{algo}
\end{algorithm}

Note that $\norm{h}_{\infty} := \max_{x\in[0,L]}h(x)$. In addition to a numerical tolerance criterion on the magnitude of the gradient, we need to take into account a constraint on the height $h$, to guarantee that the simulated flow remains in a subcritical regime, see Remark~\ref{rem:flu} (and in the range of industrial constraints, see~\cite{Chiaramonti2013}).

A similar algorithm can be considered to tackle the multiple trajectories Problem~\eqref{eq:muNza}. Remark that since no interaction between trajectories is considered, the gradient computation can be partially parallelized when computing $I_i$ (hence $C_i$) and $p_i$.

Finally, the optimality system presented in Section~\ref{sec:nonconstant} can still be solved with Algorithm~\ref{algo} by adapting the objective functional $\tilde J$ and the gradient $\nabla \tilde J$.

\subsubsection{Numerical Solvers}

To solve our optimization problem numerically, we introduce a supplementary space discretization with respect to $x$. In this way, let us take a space increment $\Delta x$, set $N_x =[L/\Delta x]$ and $x^{n_x}=n_x\Delta x$ for $n_x=0,\ldots,N_x$. We choose to use the Heun's method for computing $C$ via~\eqref{eq:cond}. Following a first-discretize-then-optimize strategy, we get that the Lagrange multiplier $p$ is also computed by a
Heun's type scheme.
Note that this scheme is still explicit since it solves a backward dynamics starting from $p(L) = 0$.

\subsection{Parameter settings}

We now detail the parameters used in our simulations.

\subsubsection{Parameterization of $h$}

As mentioned in Proposition~\ref{prop:flattopo}, we  parameterize $h$ by a truncated Fourier series~\eqref{eq:Fourier}.
The parameters to be optimized will be the Fourier coefficients $a:=[a_1,\ldots,a_N]$. Note that we choose to fix $a_0$, since it is related to the volume $V$ of the raceway with $V=a_0L$

\subsubsection{Parameter for the models}

The spatial increment is set to $\Delta x=\SI{0.01}{m}$ so that the convergence of the numerical scheme has been ensured, and we set the raceway length $L=\SI{100}{m}$, the averaged discharge $Q_0 = \SI{0.04}{m^2.s^{-1}}$, the initial value $a_0(=h(0;a))=\SI{0.4}{m}$ and $z_b(0)=-\SI{0.4}{m}$ to stay in standard ranges for a raceway. The free-fall acceleration $g=\SI{9.81}{m.s^{-2}}$. All the numerical parameters values for Han's model are taken from~\cite{Grenier2020} and given in Table~\ref{tab:parameter}.
\begin{table}[htpb]
\caption{Parameter values for Han Model.}
\begin{center}
\begin{tabular}{|c|c|c|}
\hline
$k_r$  & 6.8 $10^{-3}$ & s$^{-1}$\\
\hline
$k_d$ & 2.99 $10^{-4}$  & -\\
\hline
$\tau$ & 0.25 & s\\
\hline
$\sigma$ & 0.047 & m$^2\cdot(\mu$mol)$^{-1}$\\
\hline
$k$  & 8.7 $10^{-6}$ & -\\
\hline
$R$ &  1.389 $10^{-7}$ & s$^{-1}$\\
\hline
\end{tabular}
\end{center}
\label{tab:parameter}
\end{table}

In order to determinate the light extinction coefficient $\varepsilon$, let us assume that only $1\%$ of light can be captured by the cells at the bottom of the raceway, i.e. for $z=z_b$, meaning that $I_b=0.01I_s$, we choose $I_s=\SI{2000}{\mu mol.m^{-2}.s^{-1}}$ which approximates the maximum light intensity, e.g., at summer in the south of France. Then $\varepsilon$ can be computed by
\begin{equation*}
\varepsilon = (1/h(0;a))\ln(I_s/I_b).
\end{equation*}

Note that in the case where $a_0$ is also a parameter to be optimized, we take from~\cite{Martinez201811} the specific light extinction coefficient of the microalgae specie $\alpha_0=\SI{0.2}{m^2\cdot g}$C the background turbidity and $\alpha_1=\SI{10}{m^{-1}}$, under these settings, $a^*_0$ defined in~\eqref{eq:a0} is around $\SI{0.316}{m}$.

\subsection{Numerical results}

We now test the influence of various parameters on optimal topographies. In all our experiments, we always observe that the topographies satisfy $\norm{h}_\infty> h_c$.

\subsubsection{Influence of vertical discretization}

The first test consists in studying the influence of the vertical discretization number $N_z$. We choose $N=5$ and take 100 random initial guesses of $a$. Note that the choice of $a$ should respect the subcritical condition. Let $N_z$ varies from $1$ to $100$, and we compute the average value of $\bar \mu_{\Delta}$ for each $N_z$. The results are shown in Fig.~\ref{fig:cvNz}.
\begin{figure}[thpb]
\centering
\includegraphics[scale=0.4]{Is_2000_q_001_N_5_L_100_dx_001_convergence_Nz_optim_h_x}
\caption{The value of the functional $\bar \mu_{\Delta}$ for $N_z=[1,100]$.}
\label{fig:cvNz}
\end{figure}
We observe numerical convergence when $N_z$ grows, showing the convergence towards the model continuous in space. In view of these results, we take hereafter $N_z=50$.

\subsubsection{Influence of the initial condition}\label{sec:testC0}

The second test consists in studying the influence of the initial condition $C_0$ on the optimal shape of the raceway pond. We set the numerical tolerance Tol$=10^{-10}$, and consider $N=5$ terms in the truncated Fourier series as an example to show the shape of the optimal topography. As an initial guess, we consider the flat topography, meaning that $a$ is set to 0. Besides, we compare the optimal topographies obtained with $C_0=0.1$ and with $C_0=0.9$. The result is shown in Fig.~\ref{fig:C0depend}.
\begin{figure}[thpb]
\centering
\includegraphics[scale=0.4]{Is_2000_q_001_N_5_Nz_50_L_100_dx_001_C0_01_optim_topo_optim_h_x.pdf}
\includegraphics[scale=0.4]{Is_2000_q_001_N_5_Nz_50_L_100_dx_001_C0_09_optim_topo_optim_h_x.pdf}
\caption{The optimal topography for $C_0=0.1$ (top) and $C_0=0.9$ (bottom). The red thick line represents the topography ($z_b$), the blue thick line represents the free surface ($\eta$), and all the other curves between represent the different trajectories.}
\label{fig:C0depend}
\end{figure}
A slight difference between the two optimal topographies is observed. Note that this difference remains when the spatial increment $\Delta x$ goes to zero.

\subsubsection{Influence of Fourier series truncation}

The next test is dedicated to the influence of the order of the truncation of the Fourier series used to parameterize the water elevation $h$.

Set $N=[0,5,10,15,20]$ and keep all the other parameter values. Table~\ref{tab:N} shows the optimal value of the objective function $\bar \mu_{\Delta}(a^*)$ for different order of truncation $N$ of the Fourier series.
\begin{table}[htpb]
\caption{The objective function for different order of truncation of the Fourier series.}
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
$N$ & Iter & $\bar \mu_{\Delta}(\si{d^{-1}})$ & $\log_{10}(\|\nabla J\|)$\\
\hline
0 & 0 & 0.94764 & $-$ \\
\hline
5 & 16 & 0.96942 & -10.368944 \\
\hline
10 & 17 & 0.97509 & -10.262891 \\
\hline
15 & 17 & 0.97734 & -10.010322 \\
\hline
20 & 18 & 0.97851 & -10.196963\\
\hline
\end{tabular}
\end{center}
\label{tab:N}
\end{table}
As expected, the result shows a slight increase of the optimal value of the objective function $\bar \mu_{\Delta}(a^*)$ when $N$ becomes larger. However, the corresponding values of $\bar \mu_{\Delta}(a^*)$ remain close to the one associated with a flat topography.

\subsubsection{Simulation with paddle wheel}

The last test aims to simulate a more realistic raceway pond situation, where a paddle wheel considered in the system. More precisely, we simulate several laps, with a paddle wheel that mix up the algae after each laps. The paddle wheel is modelled by a permutation matrix $P$ that rearrange the trajectories at each lap. In our test, $P$ is an anti-diagonal matrix with the entries one. This choice actually corresponds to an optimum and has been shown in~\cite{Bernard02970756}, where other choices will be investigated.

Meanwhile, this permutation matrix $P$ corresponds to the permutation
\begin{equation*}
\pi = (1 \ N_z)(2 \ N_z-1)(3 \ N_z-2)\cdots,
\end{equation*}
where we use the standard notation of cycles in the symmetric group. Note that $\pi$ is of order two. The photo-inhibition state $C$ is then set to be $2$-periodic (i.e $C^1(0)=PC^2(L)$). The details of the optimization procedure are given in Appendix B.

We choose a truncation of order $N=5$ in the Fourier series. The initial guess $a$ is still zero. Fig.~\ref{fig:2laps} presents the shape of the optimal topography and the evolution of the photo-inhibition state $C$ for two laps. In particular, we observe that the state $C$ is actually periodic for each laps. This result is actually proved in~\cite{Bernard02970756} in the case of a flat topography.
\begin{figure}[thpb]
\centering
\includegraphics[scale=0.4]{Is_2000_q_001_N_5_Nz_50_Nl_2_L_100_dx_001_C0_per_optim_topo_optim_h_matrix.pdf}
\includegraphics[scale=0.4]{Is_2000_q_001_N_5_Nz_50_Nl_2_L_100_dx_001_C0_per_C_optim_h_matrix.pdf}
\caption{The optimal topography (top) and the evolution of the photo-inhibition state $C$ (bottom) for two laps.}
\label{fig:2laps}
\end{figure}

The resulting optimal topography in this case is not flat. However, the increase in the optimal value of the objective function $\bar \mu_{\Delta}$ compared to a flat topography is around $0.217\%$, and compare to a flat topography and non permutation case is around $0.265\%$, in both cases, the increase remains small.

\section{CONCLUSIONS AND FUTURE WORKS}

A non flat topography slightly enhances the average growth rate. However the gain remains very limited and it is not clear if the difficulty to design such a pattern could be compensated by the increase in the process productivity.

A flat topography cancels the gradient of the objective functional in many situations where $C$ is assumed to be periodic. However, including a mixing device gives rise to an optimal non flat topography with a slight gain of the average growth rate.

\bibliographystyle{plain}


\section*{\bibname}}
\usepackage{apalike}

\usepackage{geometry}
\geometry{verbose,tmargin=1in,bmargin=1in,lmargin=1in,rmargin=1in}
\hypersetup{
	colorlinks,linkcolor=red,anchorcolor=blue,citecolor=blue}
\setlength{\parindent}{0pt}
\usepackage{multirow}

\usepackage{authblk}

\makeatletter
\newcommand{\printfnsymbol}[1]{%
	\textsuperscript{\@fnsymbol{#1}}%
}
\makeatother

\begin{document}

\title{Sample Complexity Bounds for Two Timescale Value-based Reinforcement Learning Algorithms}

\author[ ]{Tengyu Xu, Yingbin Liang}
\affil[ ]{Department of Electrical and Computer Engineering, The Ohio State University}
\affil[ ]{\{xu.3260, liang.889\}@osu.edu}

\date{}
\maketitle

\begin{abstract}
Two timescale stochastic approximation (SA) has been widely used in value-based reinforcement learning algorithms. In the policy evaluation setting, it can model the linear and nonlinear temporal difference learning with gradient correction (TDC) algorithms as linear SA and nonlinear SA, respectively. In the policy optimization setting, two timescale nonlinear SA can also model the greedy gradient-Q (Greedy-GQ) algorithm. In previous studies, the non-asymptotic analysis of linear TDC and Greedy-GQ has been studied in the Markovian setting, with diminishing or accuracy-dependent stepsize. For the nonlinear TDC algorithm, only the asymptotic convergence has been established. In this paper, we study the non-asymptotic convergence rate of two timescale linear and nonlinear TDC and Greedy-GQ under Markovian sampling and with accuracy-independent constant stepsize. For linear TDC, we provide a novel non-asymptotic analysis and show that it attains an $\epsilon$-accurate solution with the optimal sample complexity of $\mathcal{O}(\epsilon^{-1}\log(1/\epsilon))$ under a constant stepsize. For nonlinear TDC and Greedy-GQ, we show that both algorithms attain $\epsilon$-accurate stationary solution with sample complexity $\mathcal{O}(\epsilon^{-2})$. It is the first non-asymptotic convergence result established for nonlinear TDC under Markovian sampling and our result for Greedy-GQ outperforms the previous result orderwisely by a factor of $\mathcal{O}(\epsilon^{-1}\log(1/\epsilon))$.
\end{abstract}

\section{Introduction}

Two timescale stochastic approximation (SA) algorithms have wide applications in reinforcement learning (RL) \cite{sutton2018reinforcement}. Typically, two timescale SA algorithms involve iterations of two types of variables updated at different speeds, i.e., the stepsizes for two iterates are chosen differently so that one iterate runs much faster than the other \cite{borkar1997stochastic,borkar2009stochastic}. Such algorithms are widely used to solve both policy evaluation and policy optimization problems in RL, in which the goal of policy evaluation is to estimate the expected total reward (i.e. value function) of a target policy, and the goal of policy optimization is to search for a policy with the optimal expected total reward.

In the policy evaluation problem, temporal difference (TD) learning \cite{sutton1988learning} is one of the most widely used algorithms when a linear function class is utilized to approximate the value function. However, in the off-policy setting, in which the target policy to be evaluated is different from the behavior policy that generates samples, TD learning may diverge to infinity. To overcome such an issue, \cite{sutton2009fast} proposed the two timescale linear TD with gradient correction (TDC) algorithm, which has convergence guarantee in the off-policy setting. The two timescale linear TDC is a special case of two timescale linear SA, whose asymptotic convergence has been established in \cite{sutton2009fast,borkar2009stochastic} and \cite{yu2017convergence,tadic2004almost,yaji2016stochastic} for the i.i.d. and Markovian settings, respectively. The non-asymptotic convergence rate of two timescale linear TDC/SA has also been studied. In the i.i.d. setting, under diminishing stepsize, \cite{dalal2017finite} established the sample complexity of $\mathcal{O}(\epsilon^{-1.5})$, and an improved complexity of $\mathcal{O}(\epsilon^{-1})$ was later established in \cite{dalal2019tale}. In the Markovian setting, \cite{xu2019two} established the complexity of $\mathcal{O}(\epsilon^{-1.5}\log^2(1/\epsilon))$ under a diminishing stepsize, and \cite{gupta2019finite} established the complexity of $\mathcal{O}(\epsilon^{-1-\zeta}\log^2(1/\epsilon))$ under a $\epsilon$-dependent stepsize, where $\zeta$ can be an arbitrarily small positive constant. Recently, \cite{kaledin2020finite} provides a tighter complexity bound of $\mathcal{O}(\epsilon^{-1})$ for two timescale linear SA under a diminishing stepsize. Although having progressed significantly, existing convergence guarantee were established either under a {\em diminishing} stepsize or a drastically small {\em $\epsilon$-level} stepsize, which yield very slow convergence and are rarely used in practice.
\begin{list}{$\bullet$}{\topsep=0.ex \leftmargin=0.15in \rightmargin=0.in \itemsep =0.01in}
	\item {\em Thus, the first goal of this paper is to investigate the two timescale linear TDC {\bf under a constant stepsize (not $\epsilon$-dependent)}, which is commonly adopted in practice, and to provide the finite-sample convergence guarantee for such a case. This necessarily requires a new approach differently from the existing ones. }

\end{list}

When a nonlinear function is utilized to approximate the value function, TD learning still suffers from the divergence issue \cite{tsitsiklis1996analysis}. To address that, \cite{bhatnagar2009convergent} proposed the two timescale nonlinear TDC, which can be modeled as a two timescale nonlinear SA.
The asymptotic convergence of two timescale nonlinear SA has been well established in \cite{borkar1997stochastic,tadic2004almost,karmakar2018two}. However, the non-asymptotic convergence of two timescale nonlinear SA has only been established in the i.i.d.\ setting under some restrict assumptions such as global (local) stability and local linearizion \cite{borkar2018concentration,mokkadem2006convergence}.
So far, the non-asymptotic convergence performance of two timescale nonlinear TDC has not been studied under the general Markovian sampling.
\begin{list}{$\bullet$}{\topsep=0.ex \leftmargin=0.15in \rightmargin=0.in \itemsep =0.01in}
	\item {\em The second goal of this paper is to provide the first non-asymptotic convergence analysis for two timescale nonlinear TDC with a constant stepsize, under {\bf Markovian sampling}, and without restricted assumptions.}
\end{list}

Moreover, in the policy optimization problem, Q-learning \cite{watkins1992q} has been widely used and has achieved significant success in practice. However, in the function approximation setting, Q-learning does not have convergence guarantee \cite{baird1995residual} unless under some restricted regularity assumptions \cite{melo2008analysis, zou2019finite, cai2019neural}. In corresponding to this, \cite{maei2010gq} proposed the Greedy-GQ algorithm in the linear function approximation setting, in which the algorithm is guaranteed to converge to a locally optimal policy without restricted assumptions. Similarly to nonlinear TDC algorithms, Greedy-GQ also adopts a two timescale update scheme, and is a special case of two timescale nonlinear SA. Under single-sample update and Markovian sampling, \cite{wang2020finite} provided the non-asymptotic convergence rate of Greedy-GQ with diminishing stepsize, which achieves the complexity of $\mathcal{O}(\epsilon^{-3}\log(\epsilon^{-1}))$. However, such a rate does not attain the typical complexity order of nonconvex optimization, and can be potentially improved with a larger stepsize.
\begin{list}{$\bullet$}{\topsep=0.ex \leftmargin=0.15in \rightmargin=0.in \itemsep =0.01in}
	\item {\em The last focus of this paper is to provide an improved non-asymptotic convergence rate for two timescale Greedy-GQ under a constant stepsize.}
\end{list}

\subsection{Our Contributions}

For two timescale linear TDC, we show that it achieves the sample complexity of $\mathcal{O}(\epsilon^{-1}\log(\epsilon^{-1}))$, which has the optimal dependence on $\epsilon$ due to the lower bound given in \cite{dalal2019tale}. Such a rate has been established in \cite{kaledin2020finite}, but only under a {\em diminishing} stepsize, which is rarely used in practice due to the slow empirical performance. In contrast, our guarantee is established under a {\bf constant (not $\epsilon$-dependent) stepsize}, which is commonly used in practice. Our analysis approach leverages the mini-batch sampling for each iteration to control the convergence error, which is significantly different from that in \cite{kaledin2020finite}, and can be of independent interest.

For two timescale nonlinear TDC, we establish the first non-asymptotic convergence rate under {\bf Markovian sampling}. We show that the mini-batch two timescale nonlinear TDC algorithm achieves the sample complexity of $\mathcal{O}(\epsilon^{-2})$.

For two timescale Greedy-GQ, we
show that mini-batch two timescale Greedy-GQ with a constant stepsize and under Markovian sampling achieves the sample complexity of $\mathcal{O}(\epsilon^{-2})$. Our result orderwisely outperforms the previous result of Greedy-GQ with diminishing stepsize in \cite{wang2020finite} by a factor of $\mathcal{O}(\epsilon^{-1}\log(1/\epsilon))$.

\subsection{Related Work}
Due to the vast amount of studies on SA and value-based RL algorithms, we include here only the studies that are highly related to our work.

\textbf{Policy evaluation with {\em linear} function approximation.} In the on-policy setting, TD learning \cite{sutton1988learning} has been proposed to solve the policy evaluation problem in the linear function approximation setting. The non-asymptotic convergence rate of TD learning has been established in \cite{dalal2018finite} for the i.i.d. setting and in \cite{bhandari2018finite,srikant2019finite,hu2019characterizing} for the Markovian setting. \cite{cai2019neural} explored the linearizable structure of neural networks in the overparameterized regime, and studied the non-asymptotic convergence rate of TD learning with neural network approximation.
\cite{zou2019finite} studied the convergence rate of SARSA with linear function approximation in the Markovian setting, which can been viewed as a policy evaluation with dynamic changing transition kernel.

In the off-policy setting, GTD, GTD2 and TDC have been proposed to solve the divergence issue of TD learning \cite{sutton2008convergent,sutton2009fast,maei2011gradient}. The convergence rate of one timescale GTD and GTD2 algorithms has been established in \cite{liu2015finite} by converting the objective into a convex-concave saddle problem in the i.i.d. setting, and was further generalized to the Markovian setting in \cite{wang2017finite}. For two timescale linear TDC, in the i.i.d.\ setting, the non-asymptotic analysis was provided in \cite{dalal2017finite,dalal2019tale}. In the Markovian setting, the non-asymptotic convergence rate was first established in \cite{xu2019two} under diminishing stepsize and in \cite{gupta2019finite} under constant stepsize. The result in \cite{xu2019two} was later improved by \cite{kaledin2020finite} to achieve the optimal convergence rate.

\textbf{Policy evaluation with {\em nonlinear} function approximation.}
Two timescale nonlinear TDC is proposed by \cite{bhatnagar2009convergent}, in which a smooth nonlinear function is utilized to approximate the value function. Nonlinear TDC with i.i.d.\ samples is a special case of two time-scale nonlinear SA with martingale noise, whose asymptotic convergence has been established in \cite{bhatnagar2009convergent,maei2011gradient} by using asymptotic convergence results in nonlinear SA \cite{borkar1997stochastic,borkar2009stochastic,tadic2004almost}. Under the global/local asymptotic stability assumptions or local linearizion assumption, the non-asymptotic convergence of two timescale nonlinear SA with martingale noise has been studied in \cite{borkar2018concentration}. Under certain stability assumptions, the asymptotic convergence of two timescale nonlinear SA with Markov noise was established in \cite{karmakar2016asymptotic,karmakar2018two}.
A concurrent study \cite{qiu2020single} also investigated nonlinear TDC and obtained the same sample complexity of $\mathcal{O}(\epsilon^{-2})$ as our result. However, \cite{qiu2019finite} only considered the i.i.d. setting, whereas we considered the more general Markovian setting.

\textbf{Policy optimization with {\em linear} function approximation.}
Q-learning \cite{watkins1992q} is one of the most widely used value-based policy optimization algorithms. The asymptotic and non-asymptotic convergence have been established for Q-learning with linear function approximation in \cite{melo2008analysis} and \cite{zou2019finite}, respectively, under certain regularity assumption. Under a similar regularity assumption, \cite{cai2019neural} established the convergence rate of Q-Learning in the neural network approximation setting. However, without regularity assumptions, Q-Learning does not have convergence guarantee in the function approximation setting.
\cite{maei2010toward} proposed two timescale Greedy-GQ to solve the divergence issue of Q-Learning with linear function approximation, and the asymptotic convergence of Greedy-GQ was also established therein. Recently, \cite{wang2020finite} studied the non-asymptotic convergence rate of Greedy-GQ under diminishing stepsize in the Markovian setting. In this paper, we provide an orderwisely better convergence rate than that in \cite{wang2020finite}.

\section{Markov Decision Process}\label{sc: mdp}

Consider a Markov decision process (MDP) denoted $(\mathcal S, \mathcal A, \mathsf{P},r,\gamma)$. Here, $\mathcal S\subset \mathbb{R}^d$ is a state space, $\mathcal A$ is an action set, $\mathsf{P}=\mathsf{P}(s^\prime|s,a)$ is the transition kernel,  $r(s, a, s^\prime)$ is the reward function bounded by $r_{\max}$, and $\gamma\in(0,1)$ is the discount factor. A stationary policy $\pi$ maps a state $s\in \mathcal S$ to a probability distribution $\pi(\cdot|s)$ over the action space $\mathcal A$. At time-step $t$, suppose the process is in some state $s_t\in \mathcal S$. Then an action $a_t\in \mathcal A$ is taken based on the distribution $\pi(\cdot|s_t)$, the system transitions to a next state $s_{t+1}\in \mathcal S$ governed by the transition kernel $\mathsf{P}(\cdot|s_t,a_t)$, and a reward $r_t=r(s_t, a_t, s_{t+1})$ is received. We assume the associated Markov chain $p(s^\prime|s)=\sum_{a\in\mathcal A}p(s^\prime|s,a)\pi(a|s)$ is ergodic, and let $\mu_\pi$ be the induced stationary distribution of this MDP, i.e., $\sum_{s}p(s^\prime|s)\mu_{\pi}(s)=\mu_{\pi}(s^\prime)$. The state value function for policy $\pi$ is defined as: $V^\pi\left(s\right)=\mathbb E[\sum_{t=0}^{\infty}\gamma^t r(s_t,a_t, s_{t+1})|s_0=s,\pi]$, and the state-action value function is defined as: $Q^\pi(s,a)=\mathbb E[\sum_{t=0}^{\infty}\gamma^t r(s_t,a_t, s_{t+1})|s_0=s, a_0=a, \pi]$. It is known that $V^\pi(s)$ is the unique fixed point of the Bellman operator $T^\pi$, i.e., $V^\pi(s) = T^\pi V^\pi(s)\coloneqq r^\pi(s)+\gamma\mathbb E_{s^\prime |s} V^\pi(s^\prime)$, where $r^\pi(s)=\mathbb E_{a, s^\prime|s}r(s,a,s^\prime)$ is the expected reward of the Markov chain induced by the policy $\pi$. We take the following standard assumption for the MDP in this paper, which has also been adopted in previous works \cite{bhandari2018finite,zou2019finite,srikant2019finite,xu2019two,xu2020reanalysis}.

\begin{assumption}[Geometric ergodicity]\label{ass1}
	There exist constants $\kappa>0$ and $\rho\in(0,1)$ such that
	\begin{flalign*}
	\sup_{s\in\mathcal S}d_{TV}(\mathbb P(s_t|s_0=s),\mu_{\pi_b})\leq \kappa\rho^t, \forall t\geq 0,
	\end{flalign*}
	where $\mathbb P(s_t|s_0=s)$ is the distribution of $s_t$ conditioned on $s_0=s$ and $d_{TV}(P,Q)$ denotes the total-variation distance between the probability measures $P$ and $Q$.
\end{assumption}
\Cref{ass1} holds for any time-homogeneous Markov chain with finite state space and any uniformly ergodic Markov chain with general state space \cite{bhandari2018finite,zou2019finite,xu2019two}.

\section{Two Timescale TDC with Linear Function Approximation}\label{sc: linearTDC}
In this section we first introduce the two timescale linear TDC algorithm to solve the policy evaluation problem, and then present our convergence rate result.

\subsection{Algorithm}

When $\mathcal S$ is large or infinite, a linear function $\hat{v}(s,\theta)=\phi(s)^\top\theta$ is often used to approximate the value function $V^\pi(s)$, where $\phi(s)\in\mathbb{R}^d$ is a fixed feature vector for state $s$ and $\theta\in\mathbb{R}^d$ is a parameter vector. We can also write the linear approximation in the vector form as $\hat{v}(\theta)={\rm \Phi} \theta$, where ${\rm \Phi}$ is the $|\mathcal S|\times d$ feature matrix. Without loss of generality, we assume that the feature vector $\ltwo{\phi(s)}\leq 1$ for all $s\in\mathcal S$ and the columns of the feature matrix $\Phi$ are linearly independent.
Here we consider policy evaluation problem in the off-policy setting. Namely, a sample path $\{ (s_t,a_t,s_{t+1}) \}_{t\geq 0}$ is generated by the Markov chain according to a behavior policy $\pi_b$, but our goal is to obtain the value function of a target policy $\pi$, which is different from $\pi_b$.

To find a parameter $\theta^*\in\mathbb{R}^d$ with $ \mathbb E_{\mu_{\pi_b}}\hat{v}(s,\theta^*)= \mathbb E_{\mu_{\pi_b}}T^\pi\hat{v}(s,\theta^*)$. The linear TDC algorithm \cite{sutton2009fast} updates the parameter by minimizing the mean-square projected Bellman error (MSPBE) objective, defined as
\begin{flalign*}
J(\theta)&=\mathbb E_{\mu_{\pi_b}}[\hat{v}(s,\theta)-{\rm \Pi} T^\pi\hat{v}(s,\theta)]^2,
\end{flalign*}
where ${\rm \Pi}$ is the orthogonal projection operation onto the function space $\hat{\mathcal V}=\{\hat{v}(\theta)\  |\  \theta\in\mathbb{R}^d\  \text{and}\ \hat{v}(\cdot, \theta)=\phi(\cdot)^\top \theta   \}$. When the columns of the feature matrix $\Phi$ are linearly independent, \cite{sutton2009fast} shows that $J(\theta)$ is strongly convex and has $\theta^*=-A^{-1}b$ as its global minimum, i.e., $J(\theta^*)=0$, where $A=\mathbb E_{\mu_{\pi_b}}[(\gamma\mathbb E_\pi[\phi(s^\prime)|s]-\phi(s))\phi(s)]$ and $b=\mathbb E_{\mu_{\pi_b}}[\mathbb E_\pi[r(s,a,s^\prime)|s]\phi(s)]$. A convenient way to find $\theta^*$ is to minimize the MSPBE objective function $J(\theta)$ using the gradient descent method: $\theta_{t+1}=\theta_t - \frac{\alpha}{2} \nabla J(\theta_t)$, where $\alpha>0$ is the stepsize and the gradient $\nabla J(\theta)$ was derived by \cite{bhatnagar2009convergent} as follows:
{\begin{flalign}\label{linear-TDCupdate}
	-\frac{1}{2}\nabla J(\theta) = \mathbb E_{\mu_{\pi_b}}[\mathbb E_\pi[\delta(\theta)|s]\phi(s)]-\gamma\mathbb E_{\mu_{\pi_b}}[ \mathbb E_\pi[\phi(s^\prime)|s]\phi(s)^\top ]w(\theta),
	\end{flalign}}
where $\delta(\theta)=r(s,a,s^\prime) + \gamma \hat{v}(s^\prime,\theta) - \hat{v}(s,\theta)$ is the temporal difference error, $w(\theta)\coloneqq \mathbb E_{\mu_{\pi_b}}[\phi(s)\phi(s)^\top]^{-1}\mathbb E_{\mu_{\pi_b}}[ \mathbb E_\pi[\delta(\theta)|s]\phi(s)]$.
In practice, stochastic gradient descent (SGD) method is usually adopted to perform the update in \cref{linear-TDCupdate} approximately. However, directly sampling is not applicable to $w(\theta)$. To solve such an issue, an auxiliary parameter $w_t$ can be introduced to estimate the vector $w(\theta_t)$, i.e., $w_t\approx w(\theta_t)$, by solving a linear SA with the following corresponding ODE:
\begin{flalign*}
\dot{w}=-\mathbb E_{\mu_{\pi_b}}[\phi(s)\phi(s)^\top]w + \mathbb E_{\mu_{\pi_b}}[\mathbb E_\pi[\delta(\theta)|s]\phi(s)].
\end{flalign*}

Given $w_t$, the parameter $\theta_t$ can then be updated with a stochastic approximation of $\nabla J(\theta_t)$ obtained via directly sampling:
\begin{flalign}
\theta_{t+1}=\theta_t + \alpha \frac{1}{\lone{\mathcal B_t}}\sum_{j\in \mathcal B_t}g(\theta_t, w_t,x_j),\label{eq: linearTDC_stoc}
\end{flalign}
where $\mathcal B_t$ is the mini-batch sampled from the MDP, $g(\theta_t, w_t,x_j)=\rho(s_j, a_j)(\delta_j(\theta_t)\phi(s_j)-\gamma\phi(s_{j+1})\phi(s_j)^\top w_t)$, $\rho(s,a)=\pi(a|s)/\pi_b(a|s)$ is the importance weighting factor with $\rho_{\max}$ being its maximum value, and $x_j$ denotes the sample $(s_j,a_j,s_{j+1})$.

\Cref{al: linearTDC} is an {\bf online }algorithm based on a {\bf single sample path}. \Cref{al: linearTDC} adopts a two timescale update scheme, in which parameters $\theta_t$ and $w_t$ are updated simultaneously but with different stepsizes. Specifically, the main parameter $\theta_t$ iterates at a slow timescale with a smaller stepsize, and the auxiliary parameter $w_t$ iterates at a fast timescale with a larger stepsize. By doing so, $w_t$ can be close to $w(\theta_t)$ asymptotically, so that $\theta_t$ is updated approximately in the direction of $-\nabla J(\theta)$.
\Cref{al: linearTDC} utilizes an accuracy-independent constant stepsize, i.e., $\alpha, \beta=\mathcal{O}(1)$ for both the updates of $\theta_t$ and $w_t$, and a mini-batch of samples $\{ (s_j, a_j, s_{j+1}) \}_{i_t\leq j\leq i_t+M-1}$ are taken sequentially from the trajectory at each iteration to perform the update. As we will show later, linear TDC in this setting is guaranteed to converge to the global optimal with an arbitrary accuracy level.

\begin{algorithm}[H]
	\null
	\caption{Two Timescale Linear TDC}
	\label{al: linearTDC}
	\begin{algorithmic}[1]
		\STATE {\bfseries Input:} batch size $M$, learning rate $\alpha$ and $\beta$
		\STATE {\bfseries Sampling:} A trajectory $\{ s_j, a_j\}_{j\geq 0}$ is sampled by following the behaviour policy $\pi_b$
		\STATE \textbf{Initialization: } $\theta_0$ and $w_0$
		\FOR{$t=0,\cdots,T-1$}
		\STATE $i_t=tM$
		\STATE $w_{t+1} = w_t + \beta\frac{1}{M}\sum_{j=i_t}^{i_t+M-1}(-\phi(s_j)\phi(s_j)^\top w_t + \rho(s_j, a_j)\delta_j(\theta_t)\phi(s_j))$
		\STATE $\theta_{t+1}=\theta_t + \alpha\frac{1}{M}\sum_{j=i_t}^{i_t+M-1} \rho(s_j, a_j)(\delta_j(\theta_t)\phi(s_j)-\gamma \phi(s_{j+1})\phi(s_j)^\top w_t)$
		\ENDFOR
		\STATE {\bfseries Output:} $\theta_T$
	\end{algorithmic}
\end{algorithm}

\subsection{Convergence Analysis}
We define matrix $C=-\mathbb E_{\mu_{\pi_b}}[\phi(s)\phi(s)^\top]$. Let $\lambda_1=\lone{\lambda_{\max}(A^\top C^{-1} A)}$, $\lambda_2=\lone{\lambda_{\max}(C)}$ and $R_\theta=\ltwo{\theta^*}$. The following theorem provides the convergence rate and sample complexity of \Cref{al: linearTDC}.
\begin{theorem}\label{thm1}
	Suppose \Cref{ass1} hold. Consider \Cref{al: linearTDC} of two timescale linear TDC update. Let the stepsize $\alpha\leq \min\left\{ \frac{1}{8\lambda_1}, \frac{\lambda_1\lambda_2}{12}, \frac{\sqrt{\lambda_2\beta}}{4\sqrt{6}\rho_{\max}} \frac{\lambda_2\sqrt{\lambda_2}\beta}{16\rho^2_{\max}}, \frac{\lambda_1\lambda_2\beta}{64\rho^2_{\max}}, \frac{\lambda_1\lambda^2_2\beta}{768}  \right\}$, $\beta \leq  \min\left\{ \frac{1}{8\lambda_2}, \frac{\lambda_2}{4}\right\}$ and the batch size $M\geq 128\left( \rho^2_{\max} + \frac{1}{\lambda^2_2} \right)\frac{1+(\kappa-1)\rho}{1-\rho}\max\left\{1,\frac{8\beta+8\lambda_2\beta^2}{\lambda_1\lambda_2\alpha}, \frac{8+12\lambda_1\alpha}{\lambda_1} \right\}$. Then we have
	{\small \begin{flalign}
		&\mathbb E[\ltwo{\theta_{T}-\theta^*}^2] \leq \left(1- \frac{\min\{  \lambda_1\alpha, \lambda_2\beta \}}{8} \right)^{T}\Delta_0 + \frac{A_1}{M},\label{eq: thm1_ub}
		\end{flalign}}where $\Delta_0 = \ltwo{w_0 - w^*(\theta_0)}^2 + \ltwo{\theta_0 - \theta^*}^2$, where $A_1$ is a constant defined in \cref{eq: 32} in \Cref{sc: proflemma1}. Furthermore, let $M\geq \frac{2A_1}{\epsilon}$ and $T\geq \frac{8}{\min\{  \lambda_1\alpha, \lambda_2\beta \}}\ln\left( \frac{2\Delta_0}{\epsilon} \right)$. The total sample complexity for \Cref{al: linearTDC} to achieve an $\epsilon$-accurate optimal solution $\theta^*$, i.e., $\mathbb E[\ltwo{\theta_T-\theta^*}^2]\leq \epsilon$, is given by
	{\small \begin{flalign*}
		TM=\Theta\left( \frac{1}{\epsilon}\log\left( \frac{1}{\epsilon} \right) \right).
		\end{flalign*}}
\end{theorem}
\Cref{thm1} shows that the convergence error of \Cref{al: linearTDC} consists of two terms: the first term is the transient error decreasing at an exponential rate, and the second term is the variance error that diminishes as the batch size $M$ increases. This is in contrast to the single-sample TDC under constant stepsizes, which suffers from the variance and bias errors with order $\mathcal{O}(\beta^2/\alpha)$ \cite{gupta2019finite}. Thus, $\epsilon$-level small stepsizes $\alpha$ and $\beta$ are required in single-sample TDC to reduce the variance error to achieve the required $\epsilon$-accurate optimal solution, which can slow down the practical convergence speed significantly. In contrast, mini-batch TDC can attain high accuracy with a large constant (not $\epsilon$-level) stepsize. Our result of $\mathcal{O}(\epsilon^{-1}\log(1/\epsilon))$ achieves the optimal complexity order due to the lower bound given in \cite{dalal2019tale}. In contrast to the same sample complexity established in \cite{kaledin2020finite}, which is applicable only under a diminishing stepsize, our result given in \Cref{thm1} is applicable under the constant stepsize, which is practically preferred due to the much better performance.

We next provide a sketch of the proof for \Cref{thm1}.
\begin{proof}[{\bf Proof Sketch of \Cref{thm1}}]
	The proof of \Cref{thm1} consists of the following three steps. At $t$-th step, we call $\ltwo{\theta_t-\theta^*}^2$ as the training error and $\ltwo{w_t - w(\theta_t)}^2$ as the tracking error.

	\noindent\textbf{Step 1:} We establish the following induction relationships for the tracking error:
	\begin{flalign}\label{eq: sketch1}
	\mathbb E\left[ \ltwo{w_{t+1}-w(\theta_{t+1})}^2 \right]&\leq (1 - \Theta(\lambda_2\beta) +  \Theta(\alpha^2/\beta) )\mathbb E\left[ \ltwo{w_{t}-w(\theta_{t})}^2 \right] \nonumber\\
	&\quad + \Theta(\alpha^2/\beta + \lambda_1\alpha)\mathbb E[\ltwo{\theta_t - \theta^*}^2] + \Theta(1/M).
	\end{flalign}

	\noindent\textbf{Step 2:} We then establish the induction relationships for the training error:
	\begin{flalign}\label{eq: sketch2}
	\mathbb E\left[\ltwo{\theta_{t+1}-\theta^*}^2\right]&\leq ( 1 - \Theta(\lambda_1\alpha) + \Theta( \alpha^2))\mathbb E\left[\ltwo{\theta_{t}-\theta^*}^2\right] \nonumber\\
	&\quad + \Theta(\alpha + \alpha^2)\mathbb E\left[ \ltwo{w_{t}-w(\theta_{t})}^2 \right] + \Theta(1/M).
	\end{flalign}

	\noindent\textbf{Step 3:} Combing \cref{eq: sketch1} and \cref{eq: sketch2} and letting the stepsize $\alpha$ and $\beta$ and batch size $M$ satisfy the requirement specified in \Cref{thm1}, we establish the induction relationship of $\Delta_t =\mathbb E[\ltwo{\theta_{t}-\theta^*}^2] + \mathbb E[ \ltwo{w_{t}-w(\theta_{t})}^2]$ as follows:
	\begin{flalign}\label{eq: sketch3}
	\Delta_{t+1} \leq \left(1- \Theta(\min\{  \lambda_1\alpha, \lambda_2\beta \}) \right) \Delta_t + \Theta( 1/M).
	\end{flalign}
	Applying \cref{eq: sketch3} recursively from $t=T-1$ to $0$ yields the desired convergence result.
\end{proof}

\section{Two Timescale TDC with Nonlinear Function Approximation}\label{sc: nonlinearTDC}

In this section we first introduce the nonlinear two timescale TDC algorithm to solve the policy evaluation problem, then we provide our non-asymptotic convergence rate result.
\subsection{Algorithm}\label{sc: nonlineartd}

In this section we consider policy evaluation problem with nonlinear function approximation, in which a parameterized smooth {\em nonlinear} function $\hat{v}(s,\theta)$ is used to approximate the value function $V^\pi(s)$. \cite{bhatnagar2009convergent} proposed an algorithm to find a parameter for the approximator $\hat{v}(s,\theta)$, named nonlinear TDC.
The nonlinear TDC updates the parameter by minimizing the following mean-square projected Bellman error objective defined as:
\begin{flalign}
J(\theta)&=\mathbb E_{\mu_{\pi}}[\hat{v}(s,\theta)-{\rm \Pi}_\theta T^\pi \hat{v}(s,\theta)]^2,\label{mspbe1}
\end{flalign}
where ${\rm \Pi}_\theta$ is the orthogonal projection operation into the function space $\bar{\mathcal V}=\{\bar{V}(s,\zeta)\  |\  \zeta\in\mathbb{R}^d\  \text{and}\ \bar{v}(s, \zeta)=\phi_\theta(s)^\top \zeta \,\text{with}\, (\phi_\theta(s))_i=\nabla_{\theta_i}\hat{v}(s,\theta)   \}$.
In general, since $J(\theta)$ defined in \cref{mspbe1} is nonconvex with respect to the parameter $\theta$, finding the global minimum of $J(\theta)$ is NP-hard. However, we can still apply gradient descent method to find a local optimum (i.e., first-order stationary point) of $J(\theta)$, via updating the parameter $\theta$ iteratively as $\theta_{t+1}=\theta_t - \frac{\alpha_t}{2} \nabla J(\theta_t)$, where $\alpha_t>0$ is the stepsize and the gradient $\nabla J(\theta)$ was derived by \cite{bhatnagar2009convergent} as follows:
{\begin{flalign}\label{nonlinear-TDCupdate}
	-\frac{1}{2}\nabla J(\theta) = \mathbb E[\delta(\theta)\phi_\theta(s)]-\gamma\mathbb E[\phi_\theta(s^\prime)\phi_\theta(s)^\top ]w(\theta)-h(\theta,w(\theta)),
	\end{flalign}}
where $\delta(\theta)=r(s,a,s^\prime) + \gamma \hat{v}(s^\prime,\theta)-\hat{v}(s,\theta)$ is the temporal difference and
{\begin{flalign*}
	&w(\theta)\coloneqq \mathbb E[\phi_\theta(s)\phi_\theta(s)^\top]^{-1}\mathbb E[\delta(\theta)\phi_\theta(s)],\nonumber\\ &h(\theta,u)\coloneqq \mathbb E[(\delta(\theta)-\phi_\theta(s)^\top u)\nabla^2_\theta V_\theta(s)u].
	\end{flalign*}}
Similarly to linear TDC studied in \Cref{sc: linearTDC}, in order to estimate the gradient in \cref{nonlinear-TDCupdate}, an auxiliary parameter $w_t$ can be used to estimate the vector $w(\theta_t)$, i.e., $w_t\approx w(\theta_t)$, by solving a linear SA with the following corresponding ODE:
\begin{flalign}\label{eq: ODE1}
\dot{w}=-\mathbb E[\phi_\theta(s)\phi_\theta(s)^\top]w + \mathbb E[\delta(\theta)\phi_\theta(s)].
\end{flalign}
Given $w_t$, the parameter $\theta_t$ can then be updated with a stochastic approximation of $\nabla J(\theta_t)$ obtained via directly sampling:
\begin{flalign}
\theta_{t+1}&=\theta_t + \alpha_t\frac{1}{\lone{\mathcal B_t}}\sum_{j\in \mathcal B_t} g(\theta_t, w_t,x_j), \label{eq: nonlinearTDC_stoc}
\end{flalign}
where $\mathcal B_t$ is the minibatch sampled from the MDP, $x_j$ denotes the sample $(s_j,a_j,s_{j+1})$ and we define $g(\theta_t, w_t,x_j)=\delta_j(\theta_t)\phi_{\theta_t}(s_j)-\gamma\phi_{\theta_t}(s_{j+1})\phi_{\theta_t}(s_j)^\top w_t-h_j(\theta_t,w_t)$, where $h_j(\theta_t,w_t)=(\delta_j(\theta_t)-\phi_{\theta_t}(s_j)^\top w_t)\nabla^2_\theta V_{\theta_t}(s_j)w_t$.
The nonlinear TDC algorithm is shown in \Cref{al: nonlinearTDC}. Similarly to \Cref{al: linearTDC}, here we also use a mini-batch of samples for each update.

\begin{algorithm}[H]
	\null
	\caption{Two Time-scale Nonlinear TDC}
	\label{al: nonlinearTDC}
	\begin{algorithmic}[1]
		\STATE \textbf{Input} batch size $M$, learning rate $\alpha$ and $\beta$
		\STATE {\bfseries Sampling:} A trajectory $\{ s_j, a_j\}_{j\geq 0}$ is sampled by following the policy $\pi$
		\STATE \textbf{Initialization: } $\theta_0$ and $w_0$
		\FOR{$t=0, 1, ..., T-1$}
		\STATE $i_t=tM$

		\STATE $w_{t+1} = w_t + \beta\frac{1}{M}\sum_{j=i_t}^{i_t+M-1}(-\phi_{\theta_t}(s_j)\phi_{\theta_t}(s_j)^\top w_t + \delta_j(\theta_t)\phi_{\theta_t}(s_j))$
		\STATE $\theta_{t+1}=\theta_t + \alpha\frac{1}{M}\sum_{j=i_t}^{i_t+M-1}(\delta_j(\theta_t)\phi_{\theta_t}(s_j)-\gamma\phi_{\theta_t}(s_{j+1})\phi_{\theta_t}(s_j)^\top w_t-h_j(\theta_t,w_t))$
		\ENDFOR
		\OUTPUT $\tilde{\theta}_{\hat{T}}$ with $\hat{T}$ chosen uniformly from $\{1,\cdots,T\}$
	\end{algorithmic}
\end{algorithm}

\subsection{Convergence Analysis}

Our analysis of \Cref{al: nonlinearTDC} will be based on the following assumptions.
\begin{assumption}[Bounded feature]\label{ass2}
	For any state $s\in \mathcal{S}$ and any vector $\theta\in \mathbb{R}^d$, we have $\ltwo{\phi_\theta(s)}\leq C_\phi$, $\lone{V(s,\theta)}\leq C_v$ and $\lF{\nabla^2_\theta V(s,\theta)}\leq D_v$, where $C_\phi$, $C_v$ and $D_v$ are positive constants.
\end{assumption}
\begin{assumption}[Smoothness]\label{ass3}
	For any state $s\in \mathcal{S}$ and any vector $\theta, \theta^\prime\in \mathbb{R}^d$, we have $\lone{V(s,\theta)-V(s,\theta^\prime)}\leq L_v\ltwo{\theta-\theta^\prime}$, $\ltwo{\phi_\theta(s)-\phi_{\theta^\prime}(s)}\leq L_\phi \ltwo{\theta-\theta^\prime}$, and $\ltwo{\nabla^2_\theta V(s,\theta) - \nabla^2_\theta V(s,\theta^\prime) } \leq L_h\ltwo{\theta-\theta^\prime} $, where $L_v$, $L_\phi$, and $L_h$ are positive constants.
\end{assumption}
\begin{assumption}[Non-singularity]\label{ass4}
	For any vector $\theta\in \mathbb{R}^d$, we have $\text{eig}\{\mathbb E[\phi_\theta(s) \phi_\theta(s)^\top]\} \geq \lambda_v$, where $\lambda_v$ is a positive constant.
\end{assumption}
\begin{assumption}[Lipschitz gradient]\label{ass5}
	For any vector $\theta, \theta^\prime$ and $w,w^\prime\in \mathbb{R}^d$, and any sample $x$, we have $\ltwo{\nabla J(\theta)- \nabla J(\theta^\prime)}\leq L_J\ltwo{\theta-\theta^\prime}$, and $\ltwo{g(\theta, w,x)-g(\theta, w^\prime,x)}\leq L_e\ltwo{w-w^\prime}$ where $L_J$ and $L_e$ are positive constants.
\end{assumption}
Assumptions \ref{ass2}-\ref{ass5} are equivalent to the assumptions adopted in the original nonlinear TDC analysis \cite{bhatnagar2009convergent}, and can be satisfied by appropriately choosing the approximation function class $\hat{v}(s,\theta)$. The following theorem characterizes the converge rate and sample complexity of \Cref{al: nonlinearTDC}.

\begin{theorem}\label{thm2}
	Consider the two timescale nonlinear TDC algorithm in \Cref{al: nonlinearTDC}. Suppose Assumptions \ref{ass1}-\ref{ass5} hold. Let the stepsize $\beta\leq \min\{\frac{\lambda_v}{8C^4_\phi}, \frac{8}{\lambda_v} \}$ and $\alpha\leq \min\{\frac{1}{2L_J}, \frac{\lambda_v \beta}{8\sqrt{2}L_w L_e}, \frac{L_J \lambda^2_v\beta^2}{384 L^2_wL^2_e} \}$. We have
	\begin{flalign*}
	\mathbb E\left[\ltwo{\nabla J(\theta_{\hat{T}})}^2\right]\leq \frac{8(J(\theta_0)-\mathbb E[J(\theta_T)])}{\alpha T} + \frac{B_1\ltwo{w_0-w(\theta_0)}^2}{T} + \frac{B_2}{M},
	\end{flalign*}

	where $B_1$ and $B_2$ are constants defined in \Cref{sc: proflemma2} in \cref{eq: 33}. Furthermore, let $M\geq \frac{2	B_2}{\epsilon}$ and $T\geq \frac{2}{\epsilon}\left[ \frac{8J(\theta_0)}{\alpha} + B_1\ltwo{w_0-w(\theta_0)}^2 \right]$. The total sample complexity for \Cref{al: nonlinearTDC} to achieve an $\epsilon$-accurate stationary point, i.e., $\mathbb E\big[ \ltwo{\nabla J(\theta_{\hat{T}})}^2 \big]\leq \epsilon$, is given by
	\begin{flalign*}
	TM=\Theta\left( \frac{1}{\epsilon^2} \right).
	\end{flalign*}
\end{theorem}

\Cref{thm2} shows that the convergence error of \Cref{al: nonlinearTDC} consists of three terms: the first two terms are the transient error decreasing at a sublinear rate as $T$ increases, and the third term contains the variance and bias errors that diminish as the batch size $M$ increases.
We next provide a sketch of the proof for \Cref{thm2}.
\begin{proof}[{\bf Proof Sketch of \Cref{thm2}}]
	The proof of \Cref{thm2} consists of the following four steps.

	\noindent\textbf{Step 1:} We first provide \Cref{lemma: fixpointlip} to show that $w(\theta)$ is $L_w$-Lipschitz:
	\begin{flalign*}
	\ltwo{w(\theta)-w(\theta^\prime)}\leq L_w\ltwo{\theta - \theta^\prime},\quad \text{for all }\theta,\theta^\prime\in \mathbb{R}^d.
	\end{flalign*}
	This property is crucial for the convergence analysis of two time-scale nonlinear TDC. It indicates that if $\theta_t$ changes slowly, then $w(\theta_t)$ also changes slowly. This allows our finite time analysis to be over a slowly changing linear SA with corresponding ODE defined in \cref{eq: ODE1}, guaranteeing that $\ltwo{w_t - w(\theta_t)}^2$ is small in an amortized sense.

	\noindent\textbf{Step 2:} We then establish the induction relationships for the tracking error $\ltwo{w_t - w(\theta_t)}^2$:
	\begin{flalign}\label{eq: sketch4}
	\mathbb E\left[ \ltwo{w_{t+1}-w(\theta_{t+1})}^2 \right] &\leq \left(1 - \Theta(\lambda_v\beta)  \right)\mathbb E\left[ \ltwo{w_{t}-w(\theta_{t})}^2 \right] \nonumber\\
	&\quad + \Theta(\alpha^2/\beta)\mathbb E\left[\ltwo{\nabla J(\theta_t)}^2\right] + \Theta(1/M).
	\end{flalign}

	\noindent\textbf{Step 3:} We then establish the induction relationships for the gradient norm $\ltwo{\nabla J(\theta_t)}^2$:
	\begin{flalign}\label{eq: sketch5}
	\left(\Theta(\alpha) - \Theta(\alpha^2) \right) \mathbb E\left[\ltwo{\nabla J(\theta_t)}^2\right]&\leq \mathbb E[J(\theta_t)] - \mathbb E[J(\theta_{t+1})] + \Theta(\alpha + \alpha^2)\mathbb E\left[ \ltwo{w_{t}-w(\theta_{t})}^2 \right] +  \Theta(1/M).
	\end{flalign}

	\noindent\textbf{Step 4:} Applying \cref{eq: sketch4} and \cref{eq: sketch5} recursively from $t=T-1$ to $0$ and combing those two results together yield
	\begin{flalign*}
	\left( \Theta(\alpha) - \Theta(\alpha^2) - \Theta(\alpha^3) \right)\sum_{t=0}^{T-1}\mathbb E\left[\ltwo{\nabla J(\theta_t)}^2\right]\leq J(\theta_0) - \mathbb E[J(\theta_T)] + \Theta((\alpha + \alpha^2)/\beta)\ltwo{w_0 - w(\theta_0)}^2 + \Theta(1/M).
	\end{flalign*}
	Letting the stepsize $\alpha$ and $\beta$ and the batch size $M$ satisfies the requirement specified in \Cref{thm2}, we can then obtain the desired convergence result.
\end{proof}

\section{Policy Optimization: Greedy-GQ Algorithm}

In this section, we will provide the non-asymptotic convergence result of Greedy-GQ \cite{maei2010toward}, which is also a two timescale nonlinear SA algorithm.

Greedy-GQ was proposed in \cite{maei2010toward} to solve the divergence issue of Q-Learning in the linear function approximation setting.
In Greedy-GQ, the goal of the agent is to learn an optimal policy for the MDP with respect to the total expected discounted reward.
In the linear function approximation setting, a linear function $\hat{Q}(s,a,\theta)=\phi(s,a)^\top \theta$ is used to approximate the state-action value function $Q(s,a)$, where $\phi(s,a)\in \mathbb{R}^d$ is a fixed feature vector for state-action pair $(s,a)$ and $\theta\in \mathbb{R}^d$ is a parameter vector. Without loss of generality, we assume that the feature vector $\ltwo{\phi(s,a)}\leq 1$ for all $(s,a)\in\mathcal S\times \mathcal A$ and the columns of the feature matrix $\Phi$ are linearly independent.  In this setting, we hope to find a solution $\theta$ that satisfies
\begin{flalign}\label{eq: gq_goal}
\Pi T^{\pi_{\theta}}\hat{Q}(s,a,\theta)=\hat{Q}(s,a,\theta),\quad \text{for all } (s,a)\in \mathcal S\times\mathcal A,
\end{flalign}
where $\pi_{\theta}$ is the soft-max greedy policy with respect to the state-action value function $\hat{Q}(s,a,\theta)$, i.e., $\pi_\theta(a|s)=\frac{\exp(\tau\hat{Q}(s,a,\theta))}{\sum_{a^\prime \in \mathcal A}\exp(\tau\hat{Q}(s,a^\prime,\theta))}$, where $\tau>0$ is the temperature parameter, and $T^{\pi_{\theta}}$ denotes the Bellman operator with policy $\pi_{\theta}$.
Similarly to the TDC algorithms, Gready-GQ searches a parameter that satisfies \cref{eq: gq_goal} by minimizing a projected Bellman error objective function defined as:
\begin{flalign}\label{mspbe2}
J(\theta)=\mathbb E_{\mu_{\pi_b}}[\hat{Q}(s,a,\theta)-{\rm \Pi} T^{\pi_\theta} \hat{Q}(s,a,\theta)]^2,
\end{flalign}
where $\delta(\theta)=r(s,a,s^\prime) + \gamma \hat{Q}(s^\prime,b,\theta)-\hat{Q}(s,a,\theta)$ is the temporal difference error, with $a \sim \pi_\theta(\cdot|s)$ and $b \sim \pi_\theta(\cdot|s^\prime)$. Since $J(\theta)$ is nonconvex and smooth everywhere, we can apply gradient descent method to find a local optimal (stationary point) of the objective $J(\theta)$ via applying the update $\theta_{t+1}=\theta_t - \frac{\alpha}{2}\nabla J(\theta_t)$ iteratively, in which
\begin{flalign*}
\frac{1}{2}\nabla J(\theta_t)&=-\mathbb E_{\mu_{\pi_b}}[\mathbb E_{\pi_\theta}[\delta(\theta)|s,a]\phi(s,a)] + \gamma \mathbb E_{\mu_{\pi_b}}[\mathbb E_{\pi_\theta}[\phi(s^\prime,b)|s,a]\phi(s,a)^\top] w(\theta)
\end{flalign*}
where
{\small \begin{flalign*}
	w(\theta)=\mathbb E_{\mu_{\pi_b}}[\phi(s,a)\phi(s,a)^\top]^{-1}\mathbb E_{\mu_{\pi_b}}[\mathbb E_{\pi_\theta}[\delta(\theta)|s,a]\phi(s,a)].
	\end{flalign*} }
Similarly to the nonlinear TDC algorithms in \cref{sc: nonlinearTDC}, here an auxiliary parameter $w_t$ is adopted to estimate the vector $w(\theta_t)$ by solving a linear SA with the following corresponding ODE:
{\small \begin{flalign*}
	\dot{w}=-\mathbb E_{\mu_{\pi_b}}[\phi(s,a)\phi(s,a)^\top]w + \mathbb E_{\mu_{\pi_b}}[\mathbb E_{\pi_\theta}[\delta(\theta)|s,a]\phi(s,a)].
	\end{flalign*}}
Then, $\theta_t$ can be updated via direct sampling
\begin{flalign*}
\theta_{t+1}=\theta_t + \alpha_t \frac{1}{\lone{\mathcal B_t}}\sum_{j\in \mathcal B_t} g(\theta_t, w_t,x_j),
\end{flalign*}
where $\rho_{\theta}(s,a)=\pi_\theta(a|s)/\pi_b(a|s)$ is the importance weighting factor bounded by $\rho_{\max}$ and we define $g(\theta_t, w_t,x_j)=\rho_{\theta_t}(s_{j+1}, a_{j+1})(\delta_j(\theta_t)\phi(s_j,a_j) - \gamma \phi(s_{j+1},a_{j+1}) \phi(s_j,a_j)^\top w_t)$. The two timescale Greedy-GQ algorithm is shown below.

\begin{algorithm}[H]
	\null
	\caption{Two Timescale Greedy-GQ}
	\label{al: GreedyGQ}
	\begin{algorithmic}[1]
		\STATE {\bfseries Input:} batch size $M$, learning rate $\alpha$ and $\beta$
		\STATE {\bfseries Sampling:} A trajectory $\{ s_j, a_j\}_{j\geq 0}$ is sampled by following the behaviour policy $\pi_b$
		\STATE \textbf{Initialization: } $\theta_0$ and $w_0$
		\FOR{$t=0,\cdots,T-1$}
		\STATE $i_t=tM$
		\STATE {\small $w_{t+1} = w_t + \beta\frac{1}{M}\sum_{j=i_t}^{i_t+M-1}(-\phi(s_j, a_j)\phi(s_j, a_j)^\top w_t + \rho_{\theta_j}(s_j, a_j)\delta_j(\theta_t)\phi(s_j, a_j))$}
		\STATE {\small $\theta_{t+1}=\theta_t + \alpha\frac{1}{M}\sum_{j=i_t}^{i_t+M-1} \rho_{\theta_t}(s_{j+1}, a_{j+1})(\delta_j(\theta_t)\phi(s_j)-\gamma \phi(s_{j+1}, a_{j+1})\phi(s_j, a_j)^\top w_t)$}
		\ENDFOR
		\OUTPUT $\tilde{\theta}_{\hat{T}}$ with $\hat{T}$ chosen uniformly from $\{1,\cdots,T\}$
	\end{algorithmic}
\end{algorithm}

By slightly abusing notations in \Cref{sc: nonlinearTDC}, we make the follow standard assumptions.
\begin{assumption}[Non-singularity]\label{ass6}
	We have $(\max_{\theta\in \mathbb{R}^d}\lone{\lambda_{\max}\{A^\top_\theta C^{-1} A_\theta\}})^{-1}=\lambda_1$ and $\lone{\lambda_{\max}\{C\}} = \lambda_2$, where $A_\theta=\mathbb E_{\mu_{\pi_b}}[(\gamma\mathbb E_{\pi_\theta}[\phi(s^\prime)|s]-\phi(s))\phi(s)^\top]$ and $C=-\mathbb E_{\mu_{\pi_b}}[\phi(s)\phi(s)^\top]$ and $\lambda_1$ and $\lambda_2$ are positive constants.
\end{assumption}
\begin{assumption}[Bounded importance factor]\label{ass7}
	For any state-action pair $(s,a)\in \mathcal{S}\times\mathcal{A}$ and any $\theta\in \mathbb{R}^d$, we have $\rho_\theta(s,a)\leq \rho_{\max}$, where $\rho_{\max}$ is a positive constant.
\end{assumption}
Note that \Cref{ass7} can be satisfied when the behaviour policy is non-degenerated for all states.
Moreover, \cite{wang2020finite} provide the following Lipschitz property of the gradient $\nabla J(\theta)$.
\begin{lemma}
	Suppose \Cref{ass1} and \Cref{ass6} hold, for any $\theta, \theta^\prime\in \mathbb{R}^d$, we have $\ltwo{\nabla J(\theta)-\nabla J(\theta^\prime)}\leq L_J\ltwo{\theta-\theta^\prime}$, where $L_J$ is a positive constant.
\end{lemma}
Note that the Greedy-GQ algorithm in \Cref{al: GreedyGQ} and nonlinear TDC algorithm in \Cref{al: nonlinearTDC} share similar structures. Both objectives are nonconvex and both algorithms adopt a two timescale update scheme, in which the fast timescale iteration corresponds to a linear SA and the slow time-scale iteration corresponds to a nonlinear SA.
Thus, the analysis of two time-scale nonlinear TDC in \Cref{sc: nonlinearTDC} can be extended to study the convergence rate of Greedy-GQ algorithm. The following theorem characterizes the convergence rate and sample complexity of \Cref{al: GreedyGQ}.
\begin{theorem}\label{thm3}
	Consider the two timescale Greedy-GQ algorithm in \Cref{al: GreedyGQ}. Suppose Assumptions \ref{ass1}, \ref{ass6} and \ref{ass7} hold. Let the stepsize $\beta\leq \min\{ \frac{\lambda_2}{4}, \frac{8}{\lambda_2} \}$ and $\alpha\leq \min\{ \frac{1}{8L_J},\frac{\lambda_2\sqrt{\lambda_2}}{8\sqrt{2}\rho_{\max}}\beta, \frac{L_J\lambda^3_2}{5312\rho^2_{\max}\lambda^2_1}\beta^2 \}$, and batch size $M\geq \frac{1+(\kappa-1)\rho}{1-\rho} \max\{128\left( \rho^2_{\max} + \frac{1}{\lambda^2_2} \right)[1+ \frac{\lambda^2_2\beta}{4\alpha^2}(\frac{2\beta}{\lambda_2} + 2\beta^2)], \frac{\beta^2\lambda^3_2(\rho_{\max}+1)^4}{\rho^2_{\max}\alpha^2}\}$. We have
	{\begin{flalign*}
		\mathbb E\left[\ltwo{\nabla J(\theta_{\hat{T}})}^2\right]&\leq \frac{8(J(\theta_0)-\mathbb E[J(\theta_T)])}{\alpha T} + \frac{192\rho^2_{\max}}{\lambda_2\beta}\frac{\ltwo{w_0 - w^*(\theta_0)}^2}{T} + \frac{32C_1[1+(\kappa-1)\rho]}{M(1-\rho)},
		\end{flalign*}}

	where $C_1$ is a positive constant defined in \cref{eq: 30} in \Cref{sc: appGQ}. Furthermore, let $M\geq \frac{64C_2[1+(\kappa-1)\rho]}{(1-\rho)\epsilon}$ and $T\geq \frac{2}{\epsilon}\left[ \frac{8J(\theta_0)}{\alpha} + \frac{192\rho^2_{\max}\ltwo{w_0-w(\theta_0)}^2}{\lambda_2\beta} \right]$. The total sample complexity for \Cref{al: nonlinearTDC} to achieve an $\epsilon$-accurate stationary point, i.e., $\mathbb E\big[ \ltwo{\nabla J(\theta_{\hat{T}})}^2 \big]\leq \epsilon$, is given by
	\begin{flalign*}
	TM=\Theta\left( \frac{1}{\epsilon^2} \right).
	\end{flalign*}
\end{theorem}

Similarly to \Cref{thm2}, in \Cref{thm3} we show that \Cref{al: GreedyGQ} converges to an $\epsilon$-accurate stationary point with sample complexity $\mathcal{O}(\epsilon^{-2})$. Note that \cite{wang2020finite} studied the convergence rate of two timescale Greedy-GQ with diminishing stepsize, which achieves the complexity of $\mathcal{O}(\epsilon^{-3}\log(1/\epsilon))$. \Cref{thm3} for two timescale Greedy-GQ with constant stepsize outperforms the result in \cite{wang2020finite} by a factor of $\mathcal{O}(\epsilon^{-1}\log(1/\epsilon))$, indicating that the constant stepsize can significantly improve the convergence rate of two timescale Greedy-GQ algorithm.

\section{Conclusion}

In this paper, we study the convergence rate for two timescale linear and nonlinear TDC and Greedy-GQ under Markovian sampling and constant stepsize. Specifically, we show that the complexity result of linear TDC orderwisely achieves the optimal convergence rate under a constant stepsize. Our result for nonlinear TDC is the first under Markovian sampling. Moreover, our sample complexity result of Greedy-GQ outperforms the previous result orderwisely.
For future work, it is interesting to apply more advance optimization techniques, e.g., acceleration, variance reduction, to further improve the convergence performance of the value-based RL algorithms studied in this paper.

\newpage
\bibliographystyle{apalike}