\section{Introduction}
Temporal logic is a promising tool for robotics applications and explainable AI in that it can be used to represent rich, complex task objectives in the form of human-readable logical specifications. Robot systems equipped with the capability to perform STL synthesis can implement powerful, intuitive command interfaces using STL. We are interested in developing STL synthesis for practical, real-world applications by viewing the synthesis problem as a form of probabilistic inference.

More widespread adoption of formal methods in robot systems, in our view, is limited by two main factors. First, the computational requirements of existing synthesis methods are seen as prohibitive. Second, achieving computational efficiency is seen to compromise semantic expressivity.

The approach we adopt in this paper is to address challenges in both computational efficiency and expressivity by viewing STL synthesis as probabilistic inference.
This probabilistic approach fits naturally with robotics perception-action pipelines for important information gathering tasks such as search, target tracking, and mapping.

We first present a probabilistic extension of STL, called random STL (RSTL), that is designed to support synthesis of robot trajectories that satisfy specifications defined over uncertain events in the environment.
Formally, RSTL extends STL to include \emph{uncertain} semantic labels.
The gain in expressivity is the capacity to specify tasks in a way that facilitates robust behaviour in real-world settings; uncertainty that is inherent to practical environments can be anticipated explicitly instead of handled reactively.
Synthesising trajectories is a differentiable problem that largely resembles previously proposed differentiable measures of robustness for deterministic STL synthesis.

We then present three approximate methods for computing the probability of satisfaction of RSTL formulae. Incorporating these methods, we implement synthesis using GPU-accelerated gradient ascent and output the most probable sequence of control actions to satisfy the given specification.

To evaluate our method, we report empirical results that illustrate correctness, convergence, and scalability properties. Further, our method inherits the benefits of gradient-based MPC, including the anytime property and predictable computation time per iteration.

To demonstrate practical use in common robotics scenarios, we provide case studies involving target search and occupancy grids. These examples show that desirable behaviour naturally arises, such as prioritising targets whose location uncertainty is increasing over targets that are physically proximate. They also demonstrate the use of complex predicates such occupancy grids which are prevalent robotics applications. Computational efficiency is shown to parallel recent advances in approximate synthesis for deterministic STL and, importantly, can be further improved through additional GPU hardware to enable development of highly capable robot systems.

We view the main contribution of this work as a step towards the feasibility of temporal logic for robotics in practice through the introduction of a new method for synthesis as probabilistic inference that can accommodate powerful task specifications and that has useful performance characteristics. Our work also provides the basis for interesting theoretical extensions that would allow temporal logic specifications to be integrated with estimation methods and belief-based planning.

\section{Related Work}
To model tasks in uncertain environments, one approach is to model the robots' dynamics as a discrete Markov decision process~(MDP), where each state is assigned semantic labels with corresponding uncertainty~\cite{uncertain_ltl_1,takashi_tanaka,yoo2012probabilistic,yoo2013provably}.
A natural task specification tool is linear temporal logic~(LTL), given which a product MDP~\cite{uncertain_ltl_1} is constructed from an automaton and a sequence of actions are found that maximize the probability of satisfaction.
These ideas can be extended to the continuous case by judiciously partitioning the environment~\cite{tulip,yoo2016online,lee2020hierarchical}.
However, we find that partitioning is computationally prohibitive for online operations in uncertain environments, because any change in belief about the environment would lead to invalidation and expensive recomputation.
Moreover, the construction of a product MDP is a computationally expensive operation, and improving its scalability remains an open problem.

Signal temporal logic~(STL)~\cite{donze_stl} is defined over continuous signals. Unlike in LTL, satisfaction is determined using continuous-valued \emph{robustness}~\cite{donze_monitoring}. Existing work uses STL to specify a task defined over a set of deterministic classes of conditions on the environment uncertainty, such as chance constraints or variance limits~\cite{Sadigh-RSS-16,stream,belief_space}.
However, since the robustness evaluation is not differentiable, a common approach is to synthesise solutions using \emph{mixed integer linear programming} (MILP) which scales exponentially with the size of mission horizon~\cite{Sadigh-RSS-16,raman_mpc}. To address the inherent complexity, the robustness metric is approximated to smooth the search space and find a solution using gradient ascent.

\section{Problem Formulation}

Suppose we have a robot with $N$-dimensional state~$\mathbf{x}_{t} \in \mathbb{R}^N$ and  control actions $\mathbf{u}_{t} \in \mathbb{U}$, where $\mathbb{U}$ is a continuous set of admissible control actions.
The robot's dynamics is uncertain, and is modelled by a discrete-time, continuous-space MDP~$\mathcal{P}(\mathbf{x}_{t+1} \mid \mathbf{x}_{t}, \mathbf{u}_{t})$ between $t$ and $t+1$, so that the trajectory distribution over a horizon $T$ is given by:
\begin{equation}\label{eq:dynamic_model}
    \mathcal{P}(\mathbf{X} \mid \mathbf{U} ) = \mathcal{P}(\mathbf{x}_{1}) \prod_{t=1}^{T} \mathcal{P}(\mathbf{x}_{t+1} \mid \mathbf{x}_{t}, \mathbf{u}_{t}),
\end{equation}
where $\mathbf{X}\equiv\mathbf{x}_{1}...\mathbf{x}_{T}$ and $\mathbf{U}\equiv\mathbf{u}_{1}...\mathbf{u}_{T}$.

The robot encounters a finite set of random events~$\mathcal{E}=\{E^{1}, ..., E^{M}\}$~(e.g., `object detected'),
whose probability of occurrence depends on robot's state $\mathbf{x}_{t}$ and time $t$.
We are interested in finding control actions~$\mathbf{U}^{*}$ that maximises the probability of satisfying a task~$\Phi$ defined over~$\mathcal{E}$~(e.g. `detect all objects').
Namely, this is a synthesis problem:
\begin{problem}[Synthesis] \label{problem:synthesis}
    Given the uncertain dynamics~\eqref{eq:dynamic_model}, and a task specification~$\Phi$ over a set of random events~$\mathcal{E}$ with probability of satisfaction $\mathcal{P}\left( (\mathbf{X}, t) \models \Phi \right) $, find an optimal sequence of controls~$\mathbf{U}^*$ such that the trajectory $\mathbf{X}$ maximises the probability of satisfying~$\Phi$ over time horizon~$T$:
    \begin{equation}
        \mathbf{U}^{*} = \argmax_{\mathbf{U} \in \mathbb{U}^{T}}~\mathcal{P}\left( (\mathbf{X}, t) \models \Phi \right)
        ,
    \end{equation}
    with respect to time~$t = 1$.
\end{problem}

\section{STL Synthesis as Probabilistic Inference}
In this section, we first introduce \emph{random signal temporal logic}~(RSTL), a probabilistic extension of STL that allows specification of tasks~$\Phi$ over random events~$\mathcal{E}$.
We then present a probabilistic inference formulation of Problem~\ref{problem:synthesis}.

\subsection{Random STL Formulae}\label{sec:meat:rstl}
We model the random events $\mathcal{E}=\{E^{1}, ..., E^{M}\}$ as (not necessarily independent) Bernoulli random variables that are dependent on robot's state and time, with conditional probability of occurrence:
\begin{equation}\label{eq:cond_bern}
    \mathcal{P}(E^{i} = 1 \mid \mathbf{x}_{t}, t) = \mathcal{P}^{i}(\mathbf{x}_{t}, t).
\end{equation}
In other words, each $E^{i}$ is a Bernoulli \emph{random field} over $\mathbb{R}^{N} \times \mathbb{R}^{+}$.
Given a set of random events $\mathcal{E}$, the syntax of an RSTL formula~$\Phi$ is given by:
\begin{equation}
    \Phi := E \mid \neg \Phi \mid \Phi \wedge \Psi \mid \Phi \mathcal{U}_{[t_{1}, t_{2}]} \Psi
    ,
\end{equation}
where $E \in \mathcal{E}$, and $\Psi, \Phi$ are RSTL formulae.
$\neg$ is logical negation, $\wedge$ is logical conjunction.
$\mathcal{U}$ is the temporal operator `Until', and $\Phi \mathcal{U}_{[t_1,t_2]} \Psi$ means $\Phi$ must hold true between time~$[t_1,t_2]$ until~$\Psi$.
Other operators such as~$\vee$ (disjunction), $\mathcal{F}_{[t_1,t_2]}$ (`in \textbf{F}uture', i.e., eventually) and~$\mathcal{G}_{[t_1,t_2]}$ (`\textbf{G}lobally', i.e., always) can be derived from the syntax the same way as deterministic STL~\cite{donze_stl}.
The events $E$ will be referred to as `event predicates'.

RSTL is \emph{random} in the sense that, for a given trajectory $\mathbf{X}$, the satisfaction of an RSTL formula~$\Phi$ is a Bernoulli random event.
The probability of satisfaction is the \emph{expected} rate of satisfaction computed over sampled instances of the event predicates $\hat{\mathcal{E}} = \{ e_{1}, ..., e_{M} \} \sim \mathcal{E}$:
\begin{equation}\label{eq:prob_sat_monte_carlo}

    \mathcal{P}( (\mathbf{X}, t) \models \Phi ) = \underset{\hat{\mathcal{E}} \sim \mathcal{E}}{\mathbb{E}}\textbf{Sat}(\mathbf{X}, t, \hat{\mathcal{E}}, \Phi)
    ,
\end{equation}
where $\textbf{Sat}(\mathbf{X}, t, \hat{\mathcal{E}}, \Phi)$ is a deterministic function evaluated recursively as:
\begin{equation} \begin{aligned}
    \textbf{Sat}(\mathbf{X}, t, \hat{\mathcal{E}}, E^{i})            &\equiv e^{i}(\mathbf{x}_{t}, t) \sim \mathcal{P}^{i}(\mathbf{x}_{t}, t)\\
    \textbf{Sat}(\mathbf{X}, t, \hat{\mathcal{E}}, \neg \Phi)        &\equiv \neg \textbf{Sat}(\mathbf{X}, t, \hat{\mathcal{E}}, \Phi) \\
    \textbf{Sat}(\mathbf{X}, t, \hat{\mathcal{E}}, \Phi \wedge \Psi) &\equiv \textbf{Sat}(\mathbf{X}, t, \hat{\mathcal{E}}, \Phi) \wedge \textbf{Sat}(\mathbf{X}, t, \hat{\mathcal{E}}, \Psi) \\

    \textbf{Sat}(\mathbf{X}, t, \hat{\mathcal{E}}, \Phi \mathcal{U}_{[t_{1}, t_{2}]} \Psi) &\equiv\\
    \bigvee_{\tau_{1} \in t + [t_{1}, t_{2}]}
    \bigwedge_{\tau_{2} \in t_{1} + [0, \tau_{1}]} &\textbf{Sat}(\mathbf{X}, \tau_{2}, \hat{\mathcal{E}}, \Phi) \vee \textbf{Sat}(\mathbf{X}, \tau_{1}, \hat{\mathcal{E}}, \Psi).
\end{aligned} \end{equation}
Note that $e^{i}(\mathbf{x}_{t}, t) \in \{0, 1\}$ is a sample from $E^{i}$ at robot state $\mathbf{x}_{t}$ and time $t$.

\subsection{STL Synthesis as Inference}
Since both task satisfaction and robot dynamics are probabilistic, it is natural to ask if Problem~\ref{problem:synthesis} can be solved solely within the realm of probability theory.
This is achieved by the control-as-inference paradigm~\cite{kappen,levine}, which has been shown to not only encompass existing optimal control problems, but also to enable new approaches.
We follow a similar development and present an inference formulation of Problem~\ref{problem:synthesis}.

In this formulation, the problem is modelled by the joint distribution among task satisfaction, robot trajectory, and control actions:
\begin{equation}
    \mathcal{P}( \Phi_{t}, \mathbf{X}, \mathbf{U}) = \mathcal{P}( \Phi_{t} \mid \mathbf{X}) \mathcal{P}(\mathbf{X} \mid \mathbf{U}) \mathcal{P}(\mathbf{U} ),
\end{equation}
where $\mathcal{P}( \Phi_{t} \mid \mathbf{X} ) \equiv \mathcal{P}( (\mathbf{X}, t) \models \Phi)$ denotes the probability satisfaction of $\Phi$ given $\mathbf{X}$ with respect to time $t$.

Here, $\mathcal{P}(\mathbf{U})$ is our prior belief on what the control actions should be, and is representative of the admissible control space $\mathbb{U}$ in the synthesis formulation.
For example, if $\mathcal{P}(\mathbf{U})$ is a zero-mean Gaussian prior, it is equivalent to penalising quadratic control cost.
The prior can derive from other knowledge, e.g., an imitation-learnt prior as~\cite{rowan} does for optimal control.

A balance between admissibility and probability of satisfaction is captured by the \emph{posterior} probability of control actions given that the task is satisfied:
\begin{equation}\begin{aligned}
    \mathcal{P}(\mathbf{U} \mid \Phi_{t}) &\propto \mathcal{P}(\Phi_{t} \mid \mathbf{U} ) \mathcal{P}(\mathbf{U} ) \\
                                 &= \mathcal{P}(\mathbf{U}) \int \mathcal{P}(\Phi_{t} \mid \mathbf{X}) \mathcal{P}(\mathbf{X} \mid \mathbf{U})d\mathbf{X},
\end{aligned}
\end{equation}
We thus pose Problem~\ref{problem:synthesis} as a maximum a posteriori~(MAP) inference problem:
\begin{problem}[Inference]
Given the dynamic model~\eqref{eq:dynamic_model} and an RSTL task specification~$\Phi$, find the MAP control actions $\mathbf{U}^{*}$ given the robot's trajectory satisfies $\Phi$ with respect to $t$:
\begin{equation}\label{eq:stl_map}
    \mathbf{U}^{*} = \argmax_{\mathbf{U}} \mathcal{P}(\mathbf{U} \mid \Phi_{t}).

\end{equation}
\end{problem}

\section{Approximate Gradient Ascent}
In this section, we first present approximate methods that allow analytical evaluation of~\eqref{eq:prob_sat_monte_carlo}.
We then present a gradient-ascent scheme on these approximate evaluations.
\subsection{Conditional Independence Approximation}
To compute~\eqref{eq:prob_sat_monte_carlo} analytically, we observe that~\eqref{eq:prob_sat_monte_carlo} applies logical operations to samples from Bernoulli random variables.
A convenient approximation is the product relation for independent Bernoulli random variables~$A$ and~$B$:
\begin{equation}\label{eq:independence}
    \mathcal{P}(A \wedge B) = \mathcal{P}(A)\mathcal{P}(B).
\end{equation}

Technically, the product relation holds true if the operands are \emph{conditionally independent} given $\mathbf{X}$.
The conditionally independent~(CI)-approximation is defined by one of the authors~\cite{chanyeol_probabilistic} as follows:
\begin{definition}[CI-approximation~\cite{chanyeol_probabilistic}]
Given an RSTL formula $\Phi$, the CI-approximation of $\mathcal{P}(\Phi_{t} \mid x)$ is defined by:
\begin{equation}
\begin{aligned}\label{eq:prob_sat}
\mathcal{P}(E^{i}_{t} \mid \mathbf{X}) &\equiv \mathcal{P}^{i}(\mathbf{X}, t)\\
\mathcal{P}(\neg \Phi_{t} \mid \mathbf{X}) &\equiv 1 - \mathcal{P}(\Phi_{t} \mid \mathbf{X})\\
\mathcal{P}(\bigwedge_{i} \Phi_{t}^{i} \mid \mathbf{X}) &\equiv \prod_{i} \mathcal{P}(\Phi_{t}^{i} \mid \mathbf{X})\\
\mathcal{P}( (\mathcal{G}_{[t_{1}, t_{2}]}\Phi)_{t} \mid \mathbf{X}) &\equiv \prod_{\tau \in t+[t_{1}, t_{2}]}\mathcal{P}(\Phi_{\tau} \mid \mathbf{X})\\
\end{aligned}
\end{equation}
\end{definition}
\subsection{Log-odds Transform}

The output range for CI-approximation of~$\mathcal{P}$~\eqref{eq:prob_sat} is~$[0, 1]$ since it computes probability.
This can lead to numerical instability, and gradient ascent often leads to poor convergence.
A natural re-parameterisation for Bernoulli random variables is the \emph{log-odds}:
\begin{equation}\label{eq:logit_sat:pred}
    \mathcal{L}(A) = \log \frac{ \mathcal{P}(A) }{ \mathcal{P}(\neg A) }
    ,
\end{equation}

It can be shown with some algebraic manipulations that re-writing the CI rule for pairwise disjunction~$\vee$ in~\eqref{eq:prob_sat} in terms of log-odds leads to:
\begin{equation}\begin{aligned}
    \mathcal{L}(A \vee B) &= \log \frac{\mathcal{P}(A)\mathcal{P}(B) + \mathcal{P}(\neg A)\mathcal{P}(B) + \mathcal{P}(A)\mathcal{P}(\neg B)}{\mathcal{P}(\neg A)\mathcal{P}(\neg B)}  \\
                           &= \mathfrak{lse}(\mathcal{L}(A), \mathcal{L}(B), \mathcal{L}(A) + \mathcal{L}(B)) \label{eq:logit_sat:or},

\end{aligned} \end{equation}
where $A$ and $B$ are independent Bernoulli random variables, and $\mathfrak{lse}$ is the \emph{log-sum-exp} function:
\begin{equation}
    \mathfrak{lse}(\mathcal{L}_{1}, ..., \mathcal{L}_{N})  = \log \sum_{i} \exp{\mathcal{L}_{i}}
    .
\end{equation}
A series of disjunction operations~\eqref{eq:logit_sat:or} is then:
\begin{equation}\label{eq:sum_over_all_subsets} \begin{aligned}
    \mathcal{L} \left( \bigvee_{i \in I} A_{i} \right) &= \log \sum_{J \in 2^{I}} \exp \sum_{j \in J} \mathcal{L}(A_{j})
    ,
\end{aligned} \end{equation}
where $2^{I}$ denotes the power set of $I$.

Computing log-sum-exp over sum of all subsets is clearly cumbersome.
We avoid such computation by using the relationship between elementary symmetric polynomials and monic polynomials.
Observe that the summations over $j \in J$ can be taken out of the exponential as products.
Then, we have all elementary symmetric polynomials over $A_{i}$ less 1.
We thus arrive at a more compact expression:
\begin{equation}\label{eq:monic_polynomial}
    \mathcal{L} \left( \bigvee_{i} A_{i} \right) = \log \left( \prod_{i} \big(1 + \exp\mathcal{L}(A_{i})\big) - 1 \right).
\end{equation}

Now, the CI computation rules in the log-odds domain are given as follows:
\begin{definition}[CI-approximate log-odds]
Given an RSTL formula $\Phi$, the CI-approximation of log-odds of satisfaction $\mathcal{L}(\Phi_{t} \mid \mathbf{X})$ is calculated by:
\begin{equation}\label{eq:logits_sat}
\begin{aligned}
\mathcal{L}(E_{t}^{i} \mid \mathbf{X}) &\equiv \log \mathcal{P}^{i}(\mathbf{X}, t) -  \log (1 - \mathcal{P}^{i}(\mathbf{X}, t)) \\
\mathcal{L}(\neg \Phi_{t} \mid \mathbf{X}) &\equiv - \mathcal{L}(\Phi_{t} \mid \mathbf{X})\\
\mathcal{L} \left( \bigvee_{i}\Phi_{t}^{i} \mid \mathbf{X} \right) &\equiv \log \left( \prod_{i} (1 + \exp\mathcal{L}(\Phi_{t}^{i} \mid \mathbf{X})) - 1 ) \right) \\% \mathcal{L}(\Phi_{t} \wedge \Psi_{t} \mid x) &\equiv -\mathfrak{lse}(-\mathcal{L}( \Phi_{t} \mid x), -\mathcal{L}(\Psi_{t} \mid x), -\mathcal{L}( \Phi_{t} \mid x) - \mathcal{L}(\Psi_{t} \mid x))\\
\mathcal{L}( (\mathcal{F}_{[t_{1}, t_{2}]}\Phi)_{t} \mid \mathbf{X}) &\equiv \log \Big( \prod_{\mathclap{\tau \in t + [t_{1}, t_{2}]}} \big(1 + \exp\mathcal{L}(\Phi_{\tau} \mid \mathbf{X}\big) - 1 \Big). \\
\end{aligned}
\end{equation}
\end{definition}

It is interesting to note that the proposed computation rules for probability of satisfaction exhibits strong similarities to existing work on deterministic STL synthesis and model checking~\cite{donze_stl,lindermann,smooth_cumulative,belta_amgm}. In the log-odds domain, certain satisfaction (i.e. probability of 1) translates to $\infty$, certain dissatisfaction is $-\infty$, and absolute uncertainty (i.e. probability of 0.5) is $0$, which are the behaviours of spatial robustness measure for deterministic STL introduced in~\cite{donze_stl}.
Further, the log-sum-exp function has been used in deterministic STL synthesis~\cite{lindermann,smooth_cumulative} as a smooth approximation of the maximum function.
Finally, A similar expression to~\eqref{eq:monic_polynomial} was presented in~\cite{belta_amgm} as an alternative robustness measure for deterministic STL. The authors report encouragement of repeated satisfaction, which is consistent with the probability of disjunction increasing with increasing probability of the disjuncts.

Given that the approaches that do not encourage repeated satisfaction~\cite{lindermann,smooth_cumulative} still report acceptable results, we consider the following approximation for disjunction:
\begin{equation}
\begin{aligned}
        \mathcal{L}(A \vee B) &\approx \mathfrak{lse}(\mathcal{L}(A), \mathcal{L}(B))) \\
        &= \log{ \frac{\mathcal{P}(A) \mathcal{P}(\neg B) + \mathcal{P}(\neg A) \mathcal{P}( B) }{\mathcal{P}(\neg A) \mathcal{P}(\neg B)} }
  .

\end{aligned}
\end{equation}
This ignores repeated satisfaction by omitting the $\mathcal{P}(A)\mathcal{P}(B)$ term from the numerator of~\eqref{eq:logit_sat:or}.
Meanwhile, there is a potential numerical benefit that log-sum-exp can be computed numerically stably with the so-called \emph{log-sum-exp trick}, while the product in~\eqref{eq:monic_polynomial} may underflow.
We thus define the \emph{mutually exclusive}~(ME) approximation as follows.
\begin{definition}[ME approximation]
Given an RSTL formula $\Phi$, the ME approximation of $\mathcal{L}(\Phi_{t} \mid \mathbf{X})$ is calculated by:
\begin{equation}\label{eq:logits_approx_sat}
\begin{aligned}
\mathcal{L}(E_{t}^{i} \mid \mathbf{X}) &\equiv \log \mathcal{P}^{i}(\mathbf{X}, t) -  \log (1 - \mathcal{P}^{i}(\mathbf{X}, t)) \\
\mathcal{L}(\neg \Phi_{t} \mid \mathbf{X}) &\equiv - \mathcal{L}(\Phi_{t} \mid \mathbf{X})\\
\mathcal{L}(\Phi_{t} \vee \Psi_{t} \mid \mathbf{X}) &\equiv \mathfrak{lse}(\mathcal{L}( \Phi_{t} \mid x), \mathcal{L}(\Psi_{t} \mid \mathbf{X}))\\
\mathcal{L}( (\mathcal{F}_{I}\Phi)_{t} \mid \mathbf{X}) &\equiv \log \sum_{\tau \in t + I} \exp\mathcal{L}(\Phi_{\tau} \mid \mathbf{X}).
\end{aligned}
\end{equation}
\end{definition}
\begin{figure}[t]
    \centering
    \subfloat[$\mathcal{F}(A) \wedge \mathcal{F}(B)$\label{fig:soundness:FaFb}]{\includegraphics[width=0.45\columnwidth]{Figures/soundness/soundness_original.png}}
    \subfloat[$\mathcal{F}(A \wedge \mathcal{F}(B))$\label{fig:soundness:sequencing}]{\includegraphics[width=0.45\columnwidth]{Figures/soundness/soundness_sequencing.png}}

    \caption{Comparison of the CI~(green) and ME~(blue) approximation results against MC estimates~(`Empirical'). CI is exact for $\mathcal{F}(A) \wedge \mathcal{F}(B)$, while ME underestimates. For $\mathcal{F}(A \wedge \mathcal{F}(B))$, the error increases, but not significantly. Naive method's result showed numerically insignificant difference to CI, and is omitted.}
    \label{fig:soundness}
\end{figure}
\subsection{Synthesis with Gradient-Ascent}
With the probability or log-odds of satisfaction computed, we synthesise a MAP control sequence~$\mathbf{U}^*$ that maximises the posterior probability.
We use Jensen's inequality to bound the log of posterior probability~\eqref{eq:stl_map}:
\begin{equation}
\begin{aligned}
    \log \mathcal{P}( \mathbf{U} \mid \Phi_{t} ) \geq &\underset{\mathbf{X}  \sim \mathcal{P}(\mathbf{X}  \mid \mathbf{U} )}{\mathbb{E}}[\log \mathcal{P}( \Phi_{t} \mid \mathbf{X}) ) ]
                                                            +\log \mathcal{P}(\mathbf{U} ).
\end{aligned}
\end{equation}

Subsequently, we maximize the lower bound:
\begin{equation}\label{eq:expected_log_prior}
    \mathbf{U}^{*} = \argmax_{\mathbf{U}} \underset{\mathbf{X} \sim \mathcal{P}(\mathbf{X} \mid \mathbf{U})}{\mathbb{E}} \left[ \log \mathcal{P}( \Phi_{t} \mid \mathbf{U}) \right]  + \log \mathcal{P}(\mathbf{U}).
\end{equation}
\begin{figure*} [t]
    \centering
    \subfloat[$N_{s}=1$]{\includegraphics[width=0.32\textwidth]{Figures/convergence/convergence-Ns-1-upto-2000.png}}\label{fig:convergence:ns-1}

    \subfloat[$N_{s}=50$]{\includegraphics[width=0.32\textwidth]{Figures/convergence/convergence-Ns-50-upto-2000.png}}\label{fig:convergence:ns-50}
    \subfloat[$N_{s}=100$]{\includegraphics[width=0.32\textwidth]{Figures/convergence/convergence-Ns-100-upto-2000.png}}\label{fig:convergence:ns-100}
    \caption{Comparison of convergence with  number of trajectory samples. Solid lines are medians.}
    \label{fig:convergence}
\end{figure*}

The maximisation is done by gradient ascent on~\eqref{eq:expected_log_prior}.
Because the expectation in~\eqref{eq:expected_log_prior} is intractable, we replace it with an empirical mean over a $N_{s}$ number of trajectory samples, so that the $i$-th gradient ascent step is:
\begin{equation}\label{eq:gradient_update}
    \hat{\mathbf{U}}^{i+1} = \hat{\mathbf{U}}^{i} +  \frac{1}{N_{s}} \sum_{j} \frac{\partial}{\partial \hat{\mathbf{U}}^{i} } [ \log \mathcal{P}( \Phi_{t} \mid \mathbf{X}_{1:T}^{j}(\hat{\mathbf{U}}^{i} )) + \log \mathcal{P}(\hat{\mathbf{U}}^{i} ) ]
    ,
\end{equation}
where $\hat{\mathbf{U}}^{i+1}=\begin{bmatrix}\hat{\mathbf{u}}_{1}^{i}...\hat{\mathbf{u}}_{T}^{i}\end{bmatrix}$.
Each trajectory sample $\mathbf{X}^{j}(\hat{\mathbf{U}}^{i})$ is obtained by propagating the dynamic model~\eqref{eq:dynamic_model} forward in time including actuation uncertainty.
Note that, as long as the predicates' distributions and the dynamic model are differentiable, so is~\eqref{eq:gradient_update}.
For both CI and ME approximations, analytical gradients can be computed easily using autograd frameworks such as Tensorflow~\cite{tensorflow} or PyTorch~\cite{pytorch,karen_leung}.
Therefore, we do not present the expressions here.

\section{Empirical Analysis}

We evaluate the practical performance characteristics of the proposed method.
We first demonstrate that the proposed CI~\eqref{eq:logits_sat} and ME~\eqref{eq:logits_approx_sat} methods reasonably approximate the ground truth~\eqref{eq:prob_sat_monte_carlo}.
We then examine the convergence characteristics of the gradient ascent~\eqref{eq:gradient_update} solution, and demonstrate the computational benefits of GPU acceleration.

\begin{figure} [t]
    \centering
    \includegraphics[width=0.6\columnwidth]{Figures/computation_time/compute_time_new.png}
    \caption{Comparison of average computation time per gradient ascent step between combinations of CPU~(red) and GPU~(green) with CI~(upward triangle and ME~(downward triangle). GPU shows 4-fold improvement in scalability. Variance was in the order of $10^{-4}$ for all configurations.}
    \label{fig:computation}
\end{figure}

\subsection{Quality of Approximation}

In this section, we evaluate whether the CI and ME approximations compute the probability of satisfaction accurately.
Because the CI computation rule assumes independence amongst operands, we expect it to be exact if 1) all predicates are conditionally independent; 2) each predicate is independent across time and space; and 3) only one operator uses each predicate.
For example, assuming 1) and 2) hold, we expect the CI rule to be exact on $\mathcal{F}A \wedge \mathcal{F}B$, but not $\mathcal{F}(A \wedge \mathcal{F}B)$, because the disjuncts of the outer $\mathcal{F}$ operator are not independent.
The ME computation rules~\eqref{eq:logits_approx_sat} will not be exact in any case.

We validate these hypotheses by comparing against a 1000-sample Monte Carlo~(MC) approximation of RSTL probability of satisfaction~\eqref{eq:prob_sat_monte_carlo}.
We used the trajectories from the first 2000 gradient ascent steps generated from the target search scenario (Fig.~\ref{fig:two_targets}).
For simplicity, we evaluated each predicates' marginal probability independently before sampling, so that the first two conditions of exactness hold.

In Fig.~\ref{fig:soundness:FaFb}, CI (green) and ME (blue) results are compared against the MC estimate for $\mathcal{F}A \wedge \mathcal{F}B$.
It can be seen that the CI method matches the MC result as expected, while ME consistently underestimates.
This is expected, because ME does not account for multiple satisfaction.

Figure~\ref{fig:soundness:sequencing} shows comparison for $\mathcal{F}(A \wedge \mathcal{F}B)$.
As $\mathcal{F}B$ is double-counted by the outer $\mathcal{F}$, CI and ME tend to over-estimate, but not by much.
ME continues to underestimate, and CI matches the MC closely,
showing that CI and ME are reasonable approximations for practical applications.

\subsection{Convergence}\label{sec:convergence}

Gradient-based methods cannot guarantee globally optimal solutions unless the objective is convex.
We analyse the convergence characteristics of the proposed computation rules in the target search scenario~(Fig.~\ref{fig:two_targets}).
We randomly generated 100 initial conditions from the control prior $\mathcal{P}(\mathbf{U})$, and ran the gradient ascent step for 2000 iterations, with $N_{s} = 1, 50, 100$ number of trajectory samples.

Figure~\ref{fig:convergence} shows the probability of satisfaction over gradient ascent steps for naive CI~(\eqref{eq:prob_sat}, red), log-odds CI~(\eqref{eq:logits_sat}, green) and log-odds ME~(\eqref{eq:logits_approx_sat}, blue) methods with varying number of trajectory samples $N_{s}$.
It can be seen that while log-odds CI and ME methods find global and local optima, the naive CI method does not find any.
The naive CI method's failure is attributed to the computation rules~\eqref{eq:prob_sat} being bounded to $[0, 1]$ range, which causes numerical errors to build up.
Note that log-odds ME reports lower probability of satisfaction due to its underestimation property, but we found that the resulting trajectories were still similar.

With increasing $N_{s}$, the variance in probability of satisfaction is decreased.
In practice, this means the generated plan will more reliably account for state uncertainty, which is useful for, e.g., the collision avoidance scenario in Fig.~\ref{fig:indoor}.
\subsection{Computation Time}

The results in Sec.~\ref{sec:convergence} illustrate that it is important to use multiple initial conditions and more trajectory samples to ensure reliable operation.
However, this would inevitably increase the computation time as well.

GPU acceleration is a prominent means to circumvent the issue of computation time, but not all algorithms benefit from GPU acceleration.
To determine if our proposed methods benefit from GPU acceleration, we compare the computation time per gradient ascent step of our Tensorflow~\cite{tensorflow} implementation between GPU and CPU.
We used all combinations between $N_{s} = 1, 10, 50, 100$ and $N_{u} = 1, 10, 50, 100$, and computed the mean over 1000 gradient ascent steps.
We used a desktop with CPU~(Intel i5-9500) and a GPU~(NVIDIA RTX2060) to conduct the experiment.

Figure~\ref{fig:computation} shows the computation time with varying number of initial conditions~$N_u$ and the number of state samples~$N_s$.
We found that the total number of samples $N_{u} \times N_{s}$ explains all changes.
The result shows that the computation time with GPU is significantly lower than that with CPU, and that using a GPU leads to $4$-fold improvement in scalability.
This demonstrates that the proposed gradient ascent method benefits from GPU acceleration.
\section{Case Studies}

We demonstrate two example cases that illustrate the benefit of using RSTL for task specification in uncertain environments.
The proposed gradient ascent method was implemented in Tensorflow~\cite{tensorflow}.
For all examples, we consider a robot described by a bicycle dynamic model:
\begin{equation} \begin{aligned}
    \dot{\mathbf{x}}_t = \begin{bmatrix}
        \dot{x}_{t} \\
        \dot{y}_{t} \\
        \dot{\theta}_{t}
    \end{bmatrix} =
    \begin{bmatrix}
        V_{t} \cos \theta_{t} \\
        V_{t} \sin \theta_{t} \\
        \omega_{t} + \epsilon_{t}
    \end{bmatrix}

    ,
\end{aligned} \end{equation}
where $\epsilon_{t} \sim \mathcal{N}(0, \sigma_{u})$ is white Gaussian noise.
The control inputs are $\mathbf{u}_{t} = \begin{bmatrix}V_{t} & \omega_{t} \end{bmatrix}^{\top}$.

\subsection{2D Target Search}\label{sec:application:evaders}

\begin{figure}[t]
\centering

    \subfloat[$t=15$\label{fig:two_targets_good_25}]{\includegraphics[width=0.45\columnwidth]{Figures/two_targets/two-targets-good-t-15-crop.png}}
    \subfloat[$t=10$\label{fig:two_targets_bad_25}]{\includegraphics[width=0.45\columnwidth]{Figures/two_targets/two-targets-bad-t-10-crop.png}} \\

    \subfloat[$t=45$\label{fig:two_targets_good_45}]{\includegraphics[width=0.45\columnwidth]{Figures/two_targets/two-targets-good-t-45-crop.png}}
    \subfloat[$t=45$\label{fig:two_targets_bad_45}]{\includegraphics[width=0.45\columnwidth]{Figures/two_targets/two-targets-bad-t-45-crop.png}} \\
    \caption{A 2D target search scenario. The robot~(blue) is tasked with detecting both Tom~(Green) and Jerry~(Red). The global optimum with $\mathcal{P}(\Phi \mid \mathbf{x}) \approx 0.7$ ~(left column) is to detect Jerry before uncertainty grows. A local optimum with $\mathcal{P}(\Phi \mid \mathbf{x}) \approx 0.5$~(right column) prefers Tom, who is closer. Circles show 1-covariance bound. \label{fig:two_targets}}
\end{figure}

\begin{figure}[t]
    \centering

    \subfloat[$\sigma_{u}=10^{-4}$\label{fig:indoor:complex_no_noise}]{\includegraphics[width=0.48\columnwidth]{Figures/indoor/targets_not_green_until.png}}
    \subfloat[$\sigma_{u}=0.1$\label{fig:indoor:complex_noise}]{\includegraphics[width=0.48\columnwidth]{Figures/indoor/targets_not_green_until_noised_2.png}}
    \caption{Results for complex indoor mission. The robot's~(blue) task is a conjunction of `never visit Rob~(red) or Bob~(blue) before visiting sanitising station~(green), and avoid obstacles~(grey colormap)', and 'visit Rob~(red) and Bob~(blue)'. With higher actuation uncertainty~$\sigma_{u}$, the trajectory becomes further from the walls. Solid lines are the robot's nominal trajectory in the absence of noise. Transparent blue lines are the noised samples used during synthesis. The start location is top-left corner.\label{fig:indoor}}

\end{figure}

We consider a 2D target search scenario, where a robot is tasked with detecting possibly moving targets in the environment: \emph{Tom} and \emph{Jerry}. Jerry, as usual, is moving with increasing uncertainty, while Tom is stationary with high certainty. Figure~\ref{fig:two_targets} depicts an example where the mean paths for Tom and Jerry are shown in green and red. The growing uncertainty over time is shown around the mean. Note that since Tom (in green) is known to be stationary, its uncertainty does not grow over time. The robot starts at~$[0, 0]^\top$.

The task of finding Tom and Jerry can be expressed using RSTL as $\Phi = \mathcal{F}(D^{\text{Tom}}) \wedge \mathcal{F}(D^{\text{Jerry}})$~(i.e., `eventually detect Tom and eventually detect Jerry').

We model the events $D^{\text{Tom}}$ and $D^{\text{Jerry}}$ as follows.
If the location is known, the robot detects Tom and Jerry with likelihood modelled by:
\begin{equation}\label{eq:detection_likelihood}
    \mathcal{P}(D_{t} \mid \mathbf{x}_{t}, \mathbf{z}_{t}) = P_{D} \exp \left( \frac{|| \mathbf{x}_{t} - \mathbf{z}_{t} ||^{2}} {2 r_{D}^{2}} \right)
    ,
\end{equation}
where $\mathbf{z}_{t}$ is the location of the target, $r_{D}$ is the radius of detection, and $P_{D}$ controls the peak.

Tom and Jerry are modelled by a constant acceleration model, which is a linear Gaussian system.
The mean $\bar{\mathbf{z}}^{a,b}_{t}$ and $\Sigma^{a,b}_{t}$ are propagated given the robots' belief at $t=0$ using the standard prediction equations.
The marginal probability of detection accounting for their uncertainty is given by:
\begin{equation}\label{eq:pred:detection}
\begin{aligned}
    \mathcal{P}(D_{t} \mid \mathbf{x}_{t}) &= \int \exp \left( \frac{|| \mathbf{x}_{t} - \bar{\mathbf{z}}_{t} ||^{2}} {2 r_{d}^{2}} \right) \mathcal{N}(\bar{\mathbf{z}}_{t}, \Sigma^{z}_{t}) d\mathbf{z}_{t} \\
                                 &= \sqrt{2\pi^{N}} r_{d}^{N} \mathcal{N}(\bar{\mathbf{z}}_{t}, \Sigma^{z}_{t} + r_{d}^{2} I),
\end{aligned}
\end{equation}
where $\mathcal{N}(\bar{\mathbf{z}}_{t}, \Sigma^{z}_{t})$ is the multivariate Gaussian PDF.

Figure~\ref{fig:two_targets} shows global and local optima found using the log-odds CI computation rule~\eqref{eq:logits_sat}.
It can be seen that the global optimum~(Figs.~\ref{fig:two_targets_good_25} and~\ref{fig:two_targets_good_45}) is to detect Jerry first at $t=15$~(Fig.~\ref{fig:two_targets_good_25}), and to return to Tom at $t=45$~(Fig.~\ref{fig:two_targets_good_45}), while the local optimum is to detect Tom first at $t=10$~(Fig.~\ref{fig:two_targets_bad_25}) and then Jerry later.
This is because the uncertainty of Jerry grows unlike Tom, and the optimal trajectory should detect Jerry first before its uncertainty grows.
This demonstrates that RSTL naturally reasons over uncertainty.
\subsection{Complex Missions in an Indoor Environment}\label{sec:application:indoor}
Consider a nursing robot in an indoor environment, modelled as an occupancy grid $O$ such that $O_{i,j}$ is the probability of obstacle occupancy.
Collision with an obstacle is modelled by an interpolation:
\begin{equation}
    \mathcal{P}(O \mid \mathbf{x}_{t}) = \sum_{i,j} K_{ij}(\mathbf{x}_{t})O_{i, j},
\end{equation}
where $K_{ij}(\mathbf{x})$ denotes the interpolant.

The robot cares for two patients, Rob and Bob.
The doctor asks the robot to avoid obstacles, and to never visit any of the patients before visiting the sanitising station, which can be written as an RSTL formula:
\begin{equation}
    \Phi_{1} = (\neg( D^{Rob} \vee D^{Bob} ) \mathcal{U} D^{San} ) \wedge \mathcal{G}(\neg O),
\end{equation}
where $D^{Rob}$, $D^{Bob}$, and $D^{San}$ are distributed as per~\eqref{eq:pred:detection}.

Now, in addition to the previous command, the doctor asks the robot to visit the two patients:
\begin{equation}
    \Phi_{2} = \mathcal{F}(D^{Rob}) \wedge \mathcal{F}(D^{Bob}) \wedge \Phi_{1}.
\end{equation}

We created an occupancy grid from a realistic dataset commonly used in perception research~\cite{Bhoram,lan_and_me}, and compared the results with low~($10^-4 rads^{-1}$) and high~($\sigma_{u} = 0.1 \text{rad}s^{-1}$) actuation uncertainty.
The results are shown in Fig.~\ref{fig:indoor}.
In both cases, the generated trajectory is correct, visiting the sanitising station first, and then the two patients.
Interestingly, the trajectory changes drastically when control noise increases.
The path with small control noise in Fig.~\ref{fig:indoor:complex_no_noise} is aggressively close to the wall, whereas the path with control noise in Fig.~\ref{fig:indoor:complex_noise} is more conservative in that the robot keeps distance away from the wall by manoeuvring around the obstacle.
This demonstrates that the proposed probabilistic formulation enables risk-averse behaviour in STL synthesis, a crucial property for practical applications.
\section{Conclusion}
We have presented a probabilistic inference perspective on STL synthesis based on RSTL and corresponding algorithms.
Our method exhibits appealing computational and expressivity characteristics that suit practical robotics applications.
We anticipate that the inference formulation of STL synthesis presented in this paper will accelerate the development of explainable AI techniques through seamless integration between formal methods and machine learning techniques as did the optimal control-as-inference paradigm~\cite{kappen,levine}. Many avenues of future work arise from our results; one of the most exciting is integration with estimation methods~\cite{lee2019online,best2015bayesian,rhinehart2019precog} that would allow multi-robot systems to augment or replace explicit communication for behaviour coordination with trajectory predictions derived from specifications~\cite{disr,decmcts}.


\section{Introduction}
\quad\quad In geometric analysis, gradient estimates play an important role in the studying
elliptic and parabolic equations on Riemannian manifolds. In \cite{Yau}, Yau proved
that if the Ricci curvature of manifold $(M^n,g)$ satisfies ${\rm Ric}_M\geq-(n-1)K$
for some constant $K\geq 0$, then any positive harmonic function $u$ on $(M^n,g)$ satisfies
\[
\lim\limits_{B_R(x_0)}\frac{|\nabla u|}{u}\le c_n\left(\frac{1}{R}+\sqrt{K}\right),
\]
where $B_R(x_0)$ is a geodesic ball with radius $R$ and center $x_0\in M^n$. As a
consequence, any positive harmonic function on manifolds with non-negative Ricci
curvature must be constant. Later, in \cite{LY86}, Li and Yau derived gradient
estimates for parabolic equations on Riemannian manifolds with or without boundary.
Moreover, they applied parabolic gradient estimates to study upper and lower bounds
of the heat kernel, eigenvalue estimates and Betti numbers estimates on manifolds.
Motivated by Li-Yau's gradient estimate technique, in \cite{Ham93}, Hamilton
established a new gradient estimate for the heat equation. His result may allow one
to compare the temperature of two different points at the same time on compact manifolds.
In \cite{SZ06}, Souplet and Zhang extended Hamilton's estimate to a localized version.
Their result enables the comparison of temperature distribution instantaneously,
without any lag in time, even for noncompact manifolds.

On the other hand, there exist some interesting works regarding to Riemannian manifolds
with boundary. In \cite{Roger}, Chen generalized Li-Yau gradient estimates to compact
Riemannian manifolds with possible non-convex boundary. As an application, he gave a
lower bound of the first Neumann eigenvalue in terms of geometric invariants with the
boundary. In \cite{Wang97}, Wang developed Li-Yau gradient estimates for the heat
equation under Neumann boundary condition. As an application, he obtained upper and
lower bounds for the heat kernel satisfying Neumann boundary conditions on compact
manifolds with non-convex boundary. In \cite{Hsu}, Hsu proved global Li-Yau gradient
estimates for the conjugate heat equation on compact manifolds with boundary, whose
metric evolving by the Ricci flow. In \cite{Olive},  Oliv\'e proved Li-Yau gradient
estimates for the heat equation, with Neumann boundary conditions, on a compact
Riemannian submanifold with boundary,  satisfying some integral Ricci curvature
assumption. For more related works, we refer to \cite{DKN, DLT, Wu15} and the
references therein.

Besides the above results related to Neumann boundary conditions, Kunikawa and Sakurai
\cite{KS20} recently proved Yau type gradient estimates for harmonic functions with
some Dirichlet boundary condition. More precisely, they proved that
\begin{theorem}\label{keita1}
Let $(M, g)$ be an $n$-dimensional, complete Riemannian manifold with compact boundary.
Assume ${\rm Ric}_M\ge-(n-1)K$ for some constant $K\ge0$ and the mean curvature of
$\partial M$ satisfies $H_{\partial M}\ge-(n-1)\sqrt{K}$.
Let $u:B_R(\partial M)\to(0,\infty)$ be a positive harmonic function with Dirichlet
boundary condition (i.e., it is constant on the boundary). If the derivative $u_\nu$
in the direction of the outward unit normal vector $\nu$ is non-negative over $\partial M$,
then
\[
\sup\limits_{B_{R/2}(\partial M)}\frac{|\nabla u|}{u}\leq c_n\left(\frac{1}{R}+\ \sqrt[]{K}\right).
\]
Here $H_{\partial M}$ denotes the mean curvature of $\partial M$ and the Dirichlet boundary
condition means that $u$ is constant on $\partial M$.
\end{theorem}

Kunikawa and Sakurai \cite{KS20} further proved local Souplet-Zhang gradient
estimates for the heat equations with Dirichlet boundary condition.
\begin{theorem}\label{keita2}
Let $(M, g)$ be an $n$-dimensional complete Riemannian manifold with compact boundary.
Assume ${\rm Ric}_M\ge-(n-1)K$ and $H_{\partial M}\ge0$. Let $0<u<A$ for
some constant $A>0$, be a solution to the heat equation on
$Q_{R,T}(\partial M):=B_R(\partial M)\times[-T,0]$. If $u$ satisfies the Dirichlet boundary
condition, and $u_\nu\ge0$ and $\partial_tu\le0$ over $\partial M\times[-T, 0]$, then
\[
\sup\limits_{Q_{R/2, T/4}(\partial M)}\frac{|\nabla u|}{u}\le c_n\left(\frac{1}{R}+\frac{1}{\sqrt{T}}+\sqrt{K}\right)\left(1+\log\frac{A}{u}\right).
\]
\end{theorem}

Inspired by the Kunikawa-Sakurai work, in this paper we will study gradient
estimates for weighted harmonic functions and weighted heat equation with Dirichlet boundary
condition. In particular, our results improve Theorems \ref{keita1} and \ref{keita2}
on the manifold case.

Before stating our results, we fix some notations. A smooth metric measure space (also called
weighted manifold) is a triple $(M,g,e^{-f}dv)$, where $(M,g)$ is a complete $n$-dimensional
Riemannian manifold, $f$ is a smooth real-valued function on $M$ and $dv$ is the volume form
with respect to the metric $g$. Denote by $\nabla$ and $\Hess$ the gradient and Hessian operators.
The weighted Laplacian is defined by $\Delta_f:=\Delta-\nabla f\cdot\nabla$. On $(M,g,e^{-f}dv)$,
a natural generalization of the Ricci curvature is called $m$-Bakry-\'Emery Ricci curvature, defined by
\[
{\rm Ric}_f^m:={\rm Ric}+\Hess f-\frac{{\nabla f \otimes \nabla f}}{{m-n}},\quad n\le m\le\infty.
\]
When $m=n$, we regard $f$ to be constant and ${\rm Ric}_f^m={\rm Ric}$. When $m=\infty$,
we have $\infty$-Bakry-\'Emery Ricci curvature ${\rm Ric}_f={\rm Ric}_f^{\infty}$. For a complete
manifold $M$ with compact boundary $\partial M$, let $H_f:=H-\nabla f\cdot\nu$ be weighted
mean curvature on $\partial M$, where $\nu$ is the outer unit normal vector to $\partial M$
(see e.g. \cite{XMZ, LiWei}) and the mean curvature $H$ is defined with respect to $\nu$.
But if $\nu$ is the inner unit normal vector to $\partial M$, then $H_f:=H+\nabla f\cdot\nu$.
We refer the reader to \cite{Saku} for further discussions on $H_f$ and weighted manifolds
with boundary.

Let us first state gradient estimates for the $f$-harmonic function with Dirichlet boundary
condition on $(M,g,e^{-f}dv)$, which is a mild generalization of Theorem \ref{keita1}.

\begin{theorem}\label{GE1}
Let $(M,g,e^{-f}dv)$ be an $n$-dimensional smooth metric measure space with compact boundary.
Assume ${\rm Ric}_f^m\ge-(m-1)K$ and $H_f\ge-L$ for some constants $K,L\ge0$. Let
$u: B_R(\partial M)\to (0, \infty)$ be a positive $f$-harmonic function (i.e. $\Delta_fu=0$)
with Dirichlet boundary condition. If the derivative $u_\nu$ in the direction of the outward
unit normal vector $\nu$ is non-negative over $\partial M$, then
\[
\sup\limits_{B_{R/2}(\partial M)}\frac{|\nabla u|}{u}\le
c_m\left(\frac{1}{R}+L+\sqrt{K}\right).
\]
\end{theorem}
We do not know if there exists a gradient estimate under a lower bound of ${\rm Ric}_f$.
Our present proof strongly depends on a refined Kato inequality in Lemma \ref{lixd},
which seems to be not true on weighted manifolds (see \cite{Kato}). When $L=(m-1)\sqrt{K}$
and $f$ is constant, our result returns to Theorem \ref{keita1}.

Second, we prove elliptic gradient estimates for the $f$-heat equation with Dirichlet
boundary condition on $(M,g,e^{-f}dv)$ by improving the argument of \cite{DungDung19},
which seems to be new even for manifolds.
\begin{theorem}\label{GE2}
Let $(M, g, e^{-f}dv)$ be a complete smooth metric measure space with compact boundary.
Assume ${\rm Ric}_{f}\ge-(n-1)K$ for some constant $K\ge0$ and $H_f\ge 0$.
Let $0<u\le A$ for some constant $A>0$ be a solution to the $f$-heat
equation $u_t=\Delta_fu$ on $Q_{R, T}(\partial M):=B_R(\partial M)\times[-T, 0]$.
If $u$ satisfies the Dirichlet boundary condition (i.e. $u(\cdot, t)\mid_{\partial M}$
is constant for each fixed $t\in[-T, 0]$), $u_\nu\geq0$ and $\partial_tu\le0$ over
$\partial M\times[-T, 0]$, then
\[
\sup\limits_{Q_{R/2, T/2}(\partial M)}\frac{|\nabla u|}{u}
\le c_n\left(\frac{\sqrt{D}+1}{R}+\frac{1}{\sqrt{T}}+\sqrt{K}\right)\sqrt{1+\log\frac{A}{u}},
\]
where $D=1+\log A-\log(\inf\limits_{Q_{R,T}(\partial M)}u)$.
\end{theorem}

As a consequence, we obtain the following Liouville theorems.

\begin{corollary}\label{Liouhar}
Let $(M, g, e^{-f}dv)$ be a complete smooth metric measure space with compact boundary.
Assume ${\rm Ric}^m_f\ge 0$ and $H_f\ge0$.  Let
$u: M\to (0, \infty)$ be a positive $f$-harmonic function
with Dirichlet boundary condition. If $u_\nu\ge 0$ over $\partial M$, then
$u$ is constant.
\end{corollary}

\begin{corollary}\label{liouville}
Let $(M, g, e^{-f}dv)$ be a complete smooth metric measure space with compact boundary.
Assume ${\rm Ric}_f\ge 0$ and $H_f\geq0$. Let $u$ be an ancient solution to
the $f$-heat equation.
\begin{enumerate}
\item Suppose that $u>0$. If $u_\nu\geq0$ and $\partial_tu\leq0$ on $\partial M\times(-\infty, 0]$ and
$u(x,t)=e^{o(\rho(x)+{|t|})}$ near infinity, then $u$ is constant.
\item If $u_\nu\geq0$ and $\partial_tu\leq0$ on $\partial M\times(-\infty, 0]$ and
$u(x,t)={o(\rho(x)+{|t|})}$ near infinity, then $u$ is constant.
\end{enumerate}
Here, $\rho(x)=\rho_{\partial M}(x)$ denotes the Riemannian distance from the boundary.
\end{corollary}

\begin{remark}
We would like to emphasize that Corollary \ref{liouville} is better than those in Corollary 1.7 in
\cite{KS20}. Because gradient estimates in Theorem \ref{GE2} (also for manifolds)
control the growth at infinity better than those in Theorem \ref{keita2}. Indeed,
our gradient estimate is sharp in both spatial and time directions due to an example
that $M=(-\infty, 0]$, $f=2x$, $u=e^{x-t}$ is an ancient positive solution to the
$f$-heat equation, $u_\nu\ge0$, $\partial_tu\le0$ on $\partial M\times(-\infty, 0]$
and its growth near infinity is $e^{|x|+|t|}$.
\end{remark}

The rest of this paper is organized as follows. In Section \ref{sec2}, we will
introduce some basic facts on smooth metric measure space with boundary. In
particular, we will give the Laplacian comparison with boundary. In Section \ref{sec3},
we will prove gradient estimates for $f$-harmonic functions with Dirichlet boundary
condition, i.e. Theorem \ref{GE1}. In Section \ref{sec3}, we will prove Theorem
\ref{GE2} and Corollaries \ref{Liouhar} and \ref{liouville}.

\section{Preliminaries}\label{sec2}
	\setcounter{equation}{0}
\quad\quad We start to introduce some basic facts. For a smooth metric measure space
$(M^n, g, e^{-f})dv$ with boundary $\partial M$, the distance function from
the boundary is defined by
\[
\rho:=\rho_{\partial M}:=d(\cdot, \partial M).
\]
Note that $\rho$ is smooth outside of the cut locus for the boundary (see \cite{Sak17}).
The Laplacian comparison for the distance function on smooth metric measure
spaces with boundary was first proved in \cite{WZZ}; see also \cite{Saku}.
Here, we give a general statement due to \cite{Sakurai17}.
We assume that $x\in B_{R}(\partial M)$, let $B(x, d(x, \partial M))\subset M$ be the largest geodesic ball centered at $x$ and $z=\partial B(x, d(x, \partial M))\cap \partial M$. Suppose that $\gamma_{z, x}(t)$ is the geodesic curve starting from $z$ and connecting $z$ and $x$.  We see that $\gamma'_{z,x}(0)$ is the unit inner normal vector for $\partial M$ at $z$. Since $\gamma_{z,x}(d(x, \partial M))=x$ and $d(x, \partial M)\leq R$, Lemma 6.1 in \cite{Sakurai17} implies the
Laplacian comparison theorem for $\infty$-Bakry-\'{E}mery Ricci curvature, which will be used in our proof of Theorems \ref{GE1} and \ref{GE2}.

\begin{theorem}[\cite{Sakurai17}]\label{Lapcom}
Let $(M^n, g, e^{-f}dv)$ be a complete smooth metric measure space with compact boundary $\partial M$. Assume that ${\rm Ric}_f\geq-K$ and $H_f\ge-L$ for some constants $K\ge0$ and $L\in\mathbb{R}$. Then
\[
\Delta_f\rho(x)\leq KR+L
\]
for all $x\in B_R(\partial M)$.
\end{theorem}
Next, we recall the so-called Reilly formula (see \cite{Lima, PeterLi, XMZ}).
\begin{proposition}\label{Reilly}
For all $\varphi\in\mathcal{C}^\infty(M)$,
\begin{equation*}
\begin{aligned}
(|\varphi|^2)_\nu=&2\varphi_\nu(\Delta_f\varphi-\Delta_{\partial M, f}(\varphi\mid_{\partial M}-\varphi_\nu H_f))
+2g_{\partial M}(\nabla_{\partial M}\varphi\mid_{\partial M}, \nabla_{\partial M}\varphi_\nu)\\
&+{\rm II}(\nabla_{\partial M}\varphi\mid_{\partial M}, \nabla_{\partial M}\varphi_{\partial M}).
\end{aligned}
\end{equation*}
\end{proposition}

For the $f$-harmonic function, we have a useful inequality, which was
proved in Theorem 2.2 of \cite{Li05}.
\begin{lemma}[see \cite{Li05}]\label{lixd}
Let $(M, g, e^{-f}dv)$ be a smooth metric measure space with ${\rm Ric}_f^m\geq-(m-1)K$. If $u$ is a positive $f$-harmonic function and $m\geq n$, then $\phi:=|\nabla\log u|$ satisfies
\[
\Delta_f\phi\geq-(m-1)K\phi-2\frac{m-2}{m-1}\frac{\langle\nabla\phi, \nabla u\rangle}{u}+\frac{\phi^3}{n-1}.
\]
\end{lemma}

For the $f$-heat equation, we have another useful inequality, which was established in
\cite{DungDung19} for the manifold case.
	\begin{lemma}\label{eq:l1}
Under the same assumption as in Theorem \ref{GE2}, let $w={{\left| \nabla h \right|}^2}$,
where $h=\sqrt{1+\log(A/u)}$.  For any $\left( x,t \right)\in Q_{R, T}(\partial M)$,
\begin{equation}\label{eq:2.3}
\Delta_f w-{w_t}\ge-2(n-1)K w
+2\left(2h-\frac{1}{h} \right)\left\langle \nabla w,\nabla h\right\rangle+2\left(2+\frac{1}{h^2} \right){w^2}.
\end{equation}
\end{lemma}
\begin{proof}
Let $u(x,t)$ be a solution to the $f$-heat equation and $0<u\le A$ for some constant $A$  in
$Q_{R, T}(\partial M)$. Let $B=Ae$ and $h=\sqrt{\log(B/u)}=\sqrt{1+\log(A/u)}\ge 1$.
Then $u=Be^{-h^2}$ and $\log u=\log B-{{h}^{2}}$. By the $f$-heat equation, we directly
compute that
\begin{equation}\label{eq:2.2}
{{h}_{t}}=\Delta_f h+|\nabla h|^2\left( \frac{1}{h}-2h \right).
\end{equation}

On the other hand, by the Bochner formula and $\Ric_f\ge-(n-1)K$, we deduce that
\begin{equation*}
\begin{aligned}
\Delta_fw-w_t&= 2|\nabla^2h|^2+2\Ric_f(\nabla h, \nabla h)+2\langle\nabla\Delta_f h,\nabla h\rangle-w_t\\
			&\ge -2(n-1)Kw+2\langle \nabla \Delta_f h,\nabla h\rangle-w_t.
\end{aligned}
\end{equation*}
		By the equality \eqref{eq:2.2}, we obtain
\begin{equation}\label{eq:2.4}
\begin{aligned}
			\Delta_fw-{{w}_{t}}
			&\ge -2(n-1)Kw+2\left\langle \nabla \left( {{h}_{t}}+{{\left| \nabla h \right|}^{2}}\left( 2h-\frac{1}{h} \right) \right),\nabla h \right\rangle -{{w}_{t}}\\
			&\ge -2(n-1)Kw+2\left\langle \nabla \left( {{h}_{t}} \right),\nabla h \right\rangle +2\left( 2h-\frac{1}{h} \right)\left\langle \nabla \left( {{\left| \nabla h \right|}^{2}} \right),\nabla h \right\rangle\\
			&\quad+2{{\left| \nabla h \right|}^{2}}\left\langle \nabla \left( 2h-\frac{1}{h} \right),\nabla h \right\rangle  -{{w}_{t}}.
\end{aligned}
\end{equation}
Observe that $2\langle\nabla(f_t),\nabla h\rangle =(|\nabla h|^2)_t=w_t$,
and
$$\nabla \left( 2h-\frac{1}{h} \right)=2\nabla h+\frac{\nabla h}{{{h}^{2}}}=\left( 2+\frac{1}{{{h}^{2}}} \right)\nabla h.$$
Hence, \eqref{eq:2.4} implies the desired inequality.
\end{proof}

\section{Gradient estimates for harmonic functions}\label{sec3}
\quad\quad In this section we will follow an argument of \cite{KS20}
to prove Theorem \ref{GE1}.
\begin{proof}[Proof of Theorem \ref{GE1}]
Define $G:=(R^2-\rho^2)\phi$, where $\phi$ is as in Lemma \ref{lixd}. Assume $G$
obtains its maximal at $x_1\in B_R(\partial M)$.

\textbf{Case 1: }If $x_1\in B_R(\partial M)\setminus\partial M$, we
may assume that $x_1\notin {\rm Cut}(\partial M)$ by the Calabi's argument.
At $x_1$, we have
$\nabla G=0$ and $\Delta_fG\le0$. Consequently,
$$\frac{\nabla \rho^2}{R^2-\rho^2}=\frac{\nabla \phi}{\phi}, \quad -\frac{\Delta_f\rho^2}{R^2-\rho^2}+\frac{\Delta_f\phi}{\phi}-\frac{2\left\langle\nabla\rho^2, \nabla\phi\right\rangle}{(R^2-\rho^2)\phi}\leq0.$$
Hence
\begin{equation}\label{kse1}
 \frac{\Delta_f\phi}{\phi}-\frac{\Delta_f\rho^2}{R^2-\rho^2}-2\frac{|\nabla\rho^2|^2}{(R^2-\rho^2)^2}\le 0.
\end{equation}
Since  ${\rm Ric}_f^m\ge-(m-1)K$ (and hence ${\rm Ric}_f\ge-(m-1)K$),
by the Laplacian comparison of Theorem \ref{Lapcom}, we have
\begin{equation}\label{kse2}
\Delta_f\rho^2=2|\nabla\rho|^2+2\rho\Delta_f\rho\le 2+2\rho((m-1)KR+L)\le2+2LR+2(m-1)KR^2,
\end{equation}
where we used $|\nabla\rho|=1$ and $\rho\le R$.
Since $|\nabla\rho^2|=4\rho^2|\nabla\rho|^2=4\rho^2$, Lemma \ref{lixd}
together with \eqref{kse1} and \eqref{kse2} implies
\begin{equation}
\begin{aligned}\label{goodterm}
0&\ge\frac{\Delta_f\phi}{\phi}-\frac{2+2LR+2(m-1)KR^2}{R^2-\rho^2}-\frac{8\rho^2}{(R^2-\rho^2)^2}\\
&\ge-(m-1)K-\frac{2(m-2)}{m-1}\frac{\langle\nabla\phi, \nabla u\rangle}{\phi u}+\frac{\phi^2}{m-1}\\
&\quad-\frac{2+2LR+2(m-1)KR^2}{R^2-\rho^2}-\frac{8\rho^2}{(R^2-\rho^2)^2}.
\end{aligned}
\end{equation}
at $x_1$. By the definition of $G$ and
\[
\frac{\langle\nabla\phi, \nabla u\rangle}{\phi u}=\frac{2\rho\langle\nabla\phi, \nabla u\rangle}{(R^2-\rho^2)u}
\le\frac{2\rho\phi}{R^2-\rho^2},
\]
\eqref{goodterm} can be written as
\begin{equation*}
\begin{aligned}
0&\ge\frac{G^2}{m-1}-\frac{4(m-2)}{m-1}\rho G-[2+2LR+2(m-1)KR^2](R^2-\rho^2)-8\rho^2-(m-1)K(R^2-\rho^2)^2\\
&\ge\frac{G^2}{m-1}-\frac{4(m-2)}{m-1}R G-[2+2LR+2(m-1)KR^2]R^2-8R^2-(m-1)KR^4\\
&\ge\frac{G^2}{2(m-1)}-[(c_m+2LR)]R^2-3(m-1)KR^4
\end{aligned}
\end{equation*}
at $x_1$, where we used the Cauchy-Schwarz inequality. This gives
\[
G(x_1)\le c_m\left(\sqrt{1+LR}\,R+\sqrt{K}R^2\right).
\]
By the definition of $G$, we conclude
\[
\frac{3R^2}{4}\sup\limits_{B_{R/2}(\partial M)}\frac{|\nabla u|}{u}\le
c_m\left(\sqrt{1+LR}\,R+\sqrt{K}R^2\right)
\]
and the desired result follows.

\textbf{Case 2: }If $x_1\in\partial M$, the proof follows by the argument of \cite{KS20}
which is originated in \cite{LiYau2}. We include it for the completeness.
Indeed, at $x_1$, $G_\nu\geq0$, $\phi_\nu=\frac{G_\nu}{R^2}\geq0$ and hence
$(\phi^2)_\nu\geq0$. Since $u$ satisfies the Dirichlet boundary condition and
$u_\nu\geq0$, then $|\nabla u|=u_\nu$. We note that
$\phi=|\nabla u|$ and $\log u$ also satisfies Dirichlet boundary condition. So
Proposition \ref{Reilly} gives
\begin{align*}
0\le(\phi^2)_\nu&=(|\nabla\log u|^2)_\nu=2(\log u)_\nu(\Delta_f\log u-(\log u)_\nu H_f)\\
&=2\frac{u_\nu}{u}\left(-\frac{|\nabla u|^2}{u^2}-\frac{u_\nu}{u}H_f\right)=2\frac{(u_\nu)^2}{u^2}\left(-\frac{u_{\nu}}{u}-H_f\right),
\end{align*}
which implies
$$\phi(x_1)=\frac{u_\nu}{u}\le - H_f\le L.$$
Hence $G(x_1)=R^2\phi(x_1)\le LR^2$. This gives
$$\frac{3R^2}{4}\sup\limits_{B_{R/2}(\partial M)}\frac{|\nabla u|}{u}\le LR^2.$$
The proof is complete.
\end{proof}

\section{Gradient estimates for heat equations}\label{sec4}
To prove Theorem \ref{GE2}, we introduce a smooth cut-off function as in \cite{KS20}.
This cut-off function is originally used in \cite{LY86} (see also \cite{SZ06}).
	\begin{lemma}\label{eq:l2}
	 There exists a smooth cut-off function $\psi=\psi(x,t)$ supported  in
$Q_{R,T}(\partial M)$ such that
		\begin{itemize}
			\item[(i)] $\psi =\psi \left( \rho_{\partial M}(x),t \right)\equiv \psi \left( r,t \right);\psi \left( r,t \right)=1$ in ${{Q}_{\frac{R}{2},\frac{T}{2}}(\partial M)},\quad 0\le \psi \le 1.$
			\item[(ii)] $\psi$ is decreasing as a radial function in the spatial variables, and $\frac{\partial \psi }{\partial r}=0$ in ${{Q}_{{R}/{2}\;,T}(\partial M)}.$
			\item[(iii)] $\left| \frac{\partial \psi }{\partial t} \right|\frac{1}{{{\psi }^{{1}/{2}\;}}}\le \frac{C}{T}$,
\quad $\left| \frac{\partial \psi }{\partial r} \right|\le \frac{{{C}_{\varepsilon }}{{\psi }^{\varepsilon }}}{{{R}}}\quad \text{and}\quad \left| \frac{{{\partial }^{2}}\psi }{{{\partial }^{2}}{{r}^{2}}} \right|\le \frac{{{C}_{\varepsilon }}{{\psi }^{\varepsilon }}}{{{R}^{2}}}, \quad 0<\varepsilon <1.$
		\end{itemize}
	\end{lemma}

Using the cut-off function, we have
	\begin{lemma}\label{mainestimate}
If $\Phi:=\psi(\Delta_f w-{w_t})+w\Delta_f\psi-w{{\psi}_t}-2\frac{|\nabla\psi|^2}{\psi}w$, then
\[
\psi {w^2}\le c\left(\frac{D^2+1}{R^4}+\frac{1}{T^2}+K^2\right)+\frac{\Phi}{4},
\]
where $D:=1+\log A-\log(\underset{Q_{R,T}(\partial M)}{\mathop{\inf}}\,u)$. Here $c$ denotes a constant
depending only on $n$ whose value may change from line to line in the following.
	\end{lemma}
	\begin{proof}
Plugging \eqref{eq:2.3} into $\Phi$,
		\begin{align}
	\Phi&\ge -2(n-1)K\psi w-2\left( 2h-\frac{1}{h} \right)\left\langle \nabla h,\nabla \psi  \right\rangle w\nonumber\\
		&\quad+2\left( 2+\frac{1}{{{h}^{2}}} \right)\psi {{w}^{2}}+w{{\Delta_f }}\psi -w{{\psi }_{t}}-2\frac{{{\left| \nabla \psi  \right|}^{2}}}{\psi }w.\nonumber
		\end{align}
This is equivalent to
		\begin{align}
		4\psi {{w}^{2}}\le& \frac{2{{h}^{2}}}{1+2{{h}^{2}}}2(n-1)K\psi w-\frac{4h\left( 1-2{{h}^{2}} \right)}{1+2{{h}^{2}}}\left\langle \nabla h,\nabla \psi  \right\rangle w\nonumber\\
		&-\frac{2{{h}^{2}}}{1+2{{h}^{2}}}w{{\Delta_f }}\psi +\frac{2{{h}^{2}}}{1+2{{h}^{2}}}w{{\psi }_{t}}+\frac{4{{h}^{2}}}{1+2{{h}^{2}}}\frac{{{\left| \nabla \psi  \right|}^{2}}}{\psi }w+\Phi.\nonumber
		\end{align}
Since $0<\frac{2{{h}^{2}}}{1+2{{h}^{2}}}\le 1$ and $0<\frac{2}{1+2{{h}^{2}}}\le 2$,  we get
\begin{equation}\label{eq:2.5}
	4\psi {w^2}\le 2(n-1)K\psi w-\frac{4h\left( 1-2{h^2} \right)}{1+2{{h}^{2}}}\langle \nabla h,\nabla \psi\rangle w-\frac{2h^2}{1+2h^2}w\Delta_f\psi+w|{{\psi}_t}|+\frac{2{{| \nabla \psi|}^2}}{\psi }w+\Phi.
\end{equation}
In the following we will estimate each term of the right hand side of \eqref{eq:2.5}.
Since $\Ric_f\ge-(n-1)K$ and $H_f\ge 0$, by the Laplacian comparison in Theorem \ref{Lapcom},
$\Delta_f\rho\le(n-1)KR$. Using this and Lemma \ref{eq:l2}, we have
\begin{equation}\label{eq:2.6}
\begin{aligned}
-\frac{2{{h}^{2}}}{1+2{{h}^{2}}}w\Delta_f \psi &=-\frac{2{{h}^{2}}}{1+2{{h}^{2}}}w\left( {{\psi }_r}\Delta_f\rho+{{\psi }_{rr}}|\nabla \rho|^2\right)\\
&\le \frac{2{{h}^{2}}}{1+2{{h}^{2}}}w\left(|{\psi}_r|(n-1)KR+|{\psi}_{rr}|\right)\\
&\le {{\psi }^{1/2\ }}w\frac{\left| {{\psi }_{rr}} \right|}{{{\psi }^{1/2\ }}}+(n-1)KR{{\psi }^{1/2}}w\frac{\left| {{\psi }_r} \right|}{{{\psi }^{1/2}}}\\
&\le \frac{3}{5}\psi {{w}^{2}}+c\left[ {{\left( \frac{\left| {{\psi }_{rr}} \right|}{{{\psi }^{1/2\ }}} \right)}^2}
+\left((n-1)KR\frac{\left| {{\psi }_r} \right|}{{{\psi }^{1/2\ }}} \right)^2\right]\\
&\le \frac{3}{5}\psi {{w}^{2}}+\frac{c}{{{R}^4}}+cK^2.
\end{aligned}
\end{equation}
By the Young inequality, we also have
\begin{equation}\label{eq:2.7}
\begin{aligned}
-\frac{4h\left( 1-2{{h}^{2}} \right)}{1+2{{h}^{2}}}\left\langle \nabla h,\nabla \psi  \right\rangle w
&\le 4h\frac{\left| 1-2{{h}^{2}} \right|}{1+2{{h}^{2}}}\left| \nabla \psi  \right|\left| \nabla h \right|w
\\&\le 4h\left| \nabla \psi  \right|{{w}^{3/2}}=4h\left| \nabla \psi  \right|{{\psi }^{-3/4\ }}{{\left( \psi {{w}^{2}} \right)}^{3/4\ }}\\
&\le \frac{3}{5}\psi {{w}^{2}}+ch^4\frac{{{\left| \nabla \psi  \right|}^{4}}}{{{\psi }^{3}}}\\
&\le \frac{3}{5}\psi {{w}^{2}}+\frac{c{{D}^{2}}}{{{R}^{4}}},
\end{aligned}
\end{equation}
where $D:=\log B-\log( \underset{Q_{R,T}(\partial M)}{\mathop{\inf }}\,u)$.
By the Cauchy-Schwarz inequality, it is not hard to see that the following estimates hold.
We first estimate that
\begin{equation}\label{eq:2.8}
2(n-1)K\psi w\le \frac{3}{5}\psi {{w}^{2}}+{4cK^{2}};
\end{equation}
for $w\left| {{\psi }_{t}} \right|$, we estimate that
\begin{equation}\label{eq:2.9}
w\left| {{\psi }_{t}} \right|={{\psi }^{{1}/{2}\;}}w\frac{\left| {{\psi }_{t}} \right|}{{{\psi }^{{1}/{2}\;}}}
\le \frac{3}{5}{{\left( {{\psi }^{{1}/{2}\;}}w \right)}^{2}}+c{{\left( \frac{\left| {{\psi }_{t}} \right|}{{{\psi }^{{1}/{2}\;}}} \right)}^{2}}\le \frac{3}{5}\psi {{w}^{2}}+\frac{c}{T^2};
\end{equation}
we also estimate that
\begin{equation}\label{eq:2.10}
\frac{2|\nabla\psi|^2}{\psi}w=2\left( {{| \nabla \psi|}^2}{{\psi }^{-3/2}} \right)\left( {{\psi }^{1/2\;}}w \right)
\le \frac{3}{5}\psi {w^2}+c\frac{{{| \nabla \psi|}^4}}{{{\psi }^3}}\le \frac{3}{5}\psi {w^2}+\frac{c}{{R^4}}.
\end{equation}
We substitute \eqref{eq:2.6}-\eqref{eq:2.10} into the right hand side of \eqref{eq:2.5}, and get that
\[
\psi w^2\le c\left( \frac{D^2+1}{R^4}+\frac{1}{T^2}+K^2\right)+\frac{\Phi}{4}.
\]
and the result follows.
\end{proof}
Now we are ready to prove Theorem \ref{GE2} by improving the argument of \cite{KS20}.
\begin{proof}[Proof of Theorem \ref{GE2}]
We may assume $u$ is non-constant in $Q_{R/2, T/2}$. If $u$ is constant,
it is trivial. Hence we let $\psi w$ be a positive function at a point in $Q_{R/2, T/2}(\partial M)$
and achieves its maximal value at  $\left(x_1,t_1\right)$ in $Q_{R/2, T/2}(\partial M)$.
We first claim $x_1\not\in\partial M$. Indeed, by contradiction, assume that $x_1\in\partial M$.
Then at $(x_1, t_1)$, $(\psi w)_\nu\geq0$, consequently $\psi_\nu w+\psi w_\nu=\psi w_\nu\geq0$;
in particular, $w_\nu\geq0$. Since $w=|\nabla\sqrt{\log(B/u)}|^2$, and $h=\sqrt{\log(B/u)}$ also
satisfies the Dirichlet boundary condition, Proposition \ref{Reilly} implies
$0\le w_\nu=(|\nabla{h}|^2)_\nu=2{h}_\nu(\Delta_f h-h_\nu H_f)$.
Since $u$ satisfies the Dirichlet boundary condition, then $|\nabla u|=u_{\nu}$. Hence,
\[
h_\nu=-\frac{u_\nu}{2u\ \sqrt{\log(B/u)}}=-\frac{|\nabla u|}{2u\ \sqrt{\log(B/u)}}=-w^{1/2}.
\]
We now compute
\begin{equation*}
\begin{aligned}
	\Delta_fh&={\rm div}\left(-\frac{\nabla u}{2u\ \sqrt{\log(B/u)}}\right)-\left\langle\nabla f, \nabla h\right\rangle\\
	&=-\frac{\Delta u}{2u\ \sqrt{\log(B/u)}}-\frac{1}{2}\left\langle\nabla u, \nabla\left(\frac{1}{u\ \sqrt[]{\log(B/u)}}\right)\right\rangle+\left\langle\nabla f, \frac{\nabla u}{2u\ \sqrt[]{\log(B/u)}}\right\rangle\\
	&=-\frac{\Delta_f u}{2u\ \sqrt{\log(B/u)}}-\frac{1}{2}\left(-\frac{|\nabla u|^2}{u^2\sqrt{\log(B/u)}}+\frac{1}{2}\frac{|\nabla u|^2}{u^2\left(\log(B/u)\right)^{3/2}}\right)\\
	&=-\frac{u_t}{2uh}+\left(\frac{2h^2-1}{h}\right)w.
\end{aligned}
\end{equation*}
	Plugging these identities into the above inequality, we obtain
\begin{equation*}
\begin{aligned}
	0&\leq -2w^{1/2}\left(-\frac{u_t}{2uh}+\left(\frac{2h^2-1}{h}\right)w+w^{1/2}H_f\right)\\
	&\leq-2w\left(\left(\frac{2h^2-1}{h}\right)w^{1/2}+H_f\right)
\end{aligned}
\end{equation*}
where we used $u_t\leq0$ in the last inequality. This implies
\[
\left(\frac{2h^2-1}{h}\right)w^{1/2}\leq-H_f\leq0.
\]
Since $h\ge1$, $w=0$ at $(x_1, t_1)$. This means $\psi w\equiv0$ on $Q_{R,T}(\partial M)$. This is a contracdiction.

We have shown that $x_1\not\in\partial M$. By the standard argument of Calabi (see \cite{Cal57}), we may assume that $x_1\not\in\partial M\cup{\rm Cut}(\partial M)$. Now at $\left(x_1,t_1\right)$, Lemma \ref{mainestimate} gives
\begin{equation}\label{econclu}
w^2\le c\left(\frac{D^2+1}{R^4}+\frac{1}{T^2}+K^2 \right)+\frac{\Phi}{4}.
\end{equation}
On the other hand, since $(x_1, t_1)$ is a maximal point, we have
$\nabla(\psi w)=0$, $\Delta_f(\psi w)\le 0$ and $(\psi w)_t\ge 0$
at $(x_1, t_1)$. Thus, at $(x_1,t_1)$,
		$$0\ge {{\Delta_f }}\left( \psi w \right)-{{\left( \psi w \right)}_{t}}=\psi \left( {{\Delta_f }}w-{{w}_{t}} \right)+w\left( {{\Delta_f }}\psi -{{\psi }_{t}} \right)+2\left\langle \nabla w,\nabla \psi  \right\rangle,$$
that is, $\Phi(x_1, t_1)\leq0$. Hence, \eqref{econclu} implies
		$$(\psi w)(x, t)\leq (\psi w)(x_1, t_1)\leq c^{1/2}\left( \frac{D+1}{R^2}+\frac{1}{T}+K \right)$$
for all $(x, t)\in Q_{R,T}(\partial M)$. Since $\psi\equiv1$ in $Q_{R/2, T/2}(\partial M)$, by the definition of $w$, we get
$$\frac{|\nabla u|}{u}\leq 2c^{1/4}\left(\frac{\sqrt{D}+1}{R}+\frac{1}{\sqrt{T}}+\sqrt{K}\right)\sqrt{1+\log\frac{A}{u}}.$$
The proof is complete.
\end{proof}
Corollary \ref{Liouhar} follows by letting $R\to\infty$. So we only give a proof of Corollary \ref{liouville}.
\begin{proof}[Proof of Corollary \ref{liouville}]
Now $K=L=0$ and $u_t=\Delta_f u$. Let $v=u+1$ an then $v_t=\Delta_f v$. Moreover, $u$ and $v$ have
the same growth at infinity. Hence, without loss of generality, we may assume $u\geq 1$.
Fix $(x_0, t_0)$. We apply Theorem \ref{GE2} to $Q_{R,R}(\partial M)=B_R(\partial M)\times [t_0-R, t_0]$
and get
\begin{equation*}
\begin{aligned}
\frac{|\nabla u|}{u}(x_0, t_0)
&\le c_m\left( \frac{\sqrt{1+\log A}}{R}+\frac{1}{\sqrt{R}}\right)\sqrt{1+\log \frac{A}{u(x_0, t_0)}}\\
&\le c_m\left( \frac{\sqrt{o(R+|R|)}}{R}+\frac{1}{\sqrt{R}}\right)\sqrt{o(R+|R|)-\log(u(x_0, t_0))}.
\end{aligned}
\end{equation*}
Letting $R\to\infty$, we get $|\nabla u(x_0, t_0)|=0$ and $u$ is constant because $(x_0, t_0)$ is arbitrary.
The proof of (1) is complete.

The proof of (2) is similar as in \cite{SZ06}, we omit the details.
\end{proof}
\begin{remark}
In Theorem 4.2 of \cite{WZZ}, the authors proved that $\rho$ is finite
if some extra conditions on $K$ and $L$ are given. We point out that our gradient
estimates and Liouville theorems work well in this case. In fact, if $\rho$ is finite,
then all boundary balls $B_{\partial M}(R)$ are the same if $R$ is large. The proof
still work without any change. Therefore, we can let $R\to\infty$ even when $\rho$ is finite.
\end{remark}

	\vskip 0.3cm
	\noindent
\section*{Acknowledgment: }The second named author thanks Prof. Keita Kunikawa for useful comments and suggestions on weighted manifolds with $\infty$-Bakry-\'{E}mery curvature. In particular, he thanks Prof. Keita Kunikawa for pointing out Lemma 6.1 in \cite{Sakurai17} that leads to an improvement of this paper.