\section{Introduction}
\label{sec:electr-subm}

Artificial intelligence (AI) has revolutionized  maintenance operations by providing more accurate predictive tools such as failure prediction and remaining useful life (RUL) estimation in recent years \cite{zheng2017long}. However,    industries  often  rely on human judgment to use  predictive tools and  determine maintenance decisions.
\citeA{silver2016mastering}  have shown that deep reinforcement learning (RL) can achieve   a performance superior to human judgment in games such as chess and Go.  The main challenge of training deep RL for predictive maintenance decision-making is the  lack of reliable simulators in most industries \cite{khorasgani2020challenges}. Failure is often rare and complex. Therefore, it is challenging  to develop a simulator that can model these events.

Recently, several researchers have shown that deep RL  can be formulated as a supervised problem. \citeA{schmidhuber2019reinforcement} introduced Upside Down RL (UDRL). Unlike traditional
RL algorithms wherein an agent learns from  reward, UDRL learns a mapping from  desired  rewards and  observations  to  actions. The main  idea of UDRL is that supervised learning can learn a generalized  model  that can generate actions which can achieve requested rewards within requested  time. In UDRL,  learning includes two parallel algorithms: 1) A1 performs execution and exploration and 2) A2 performs  supervised learning. A1 generates actions   based on known maximum possible cumulative rewards using the current version of the model trained by A2. For the sake of exploration,  random actions are  taken   occasionally. The observations, actions and rewards are saved in a buffer. A2 uses historical data that is  extracted from the buffer  to train the model that maps the observation and maximum possible cumulative reward to optimal action.

\citeA{zha2021simplifying}  show that supervised learning can be competitive to state of the art RL  algorithms with
better  stability  and lower training  time.
\citeA{chen2021decision} introduced Decision Transformer which  expands UDRL to offline RL. Decision Transformer uses transformer architecture \cite{vaswani2017attention} to learn  a mapping from  historical observations, historical actions, and expected reward to proper actions using offline data.
\citeA{janner2021reinforcement}  also show that a supervised approach based on the transformer architecture  is competitive to  other offline RL algorithms without the requirement of learning the system dynamic model or  even the value function.
Relying purely on offline data makes offline RL  a perfect choice for maintenance decision-making wherein storing historical data to learn a policy is much more cost-effective   than developing a simulator that can model failure accurately. Using a supervised approach simplifies  the solution even more.  In this work, we use a supervised offline RL approach to learn optimal maintenance decision-making.

The rest of the paper is as follows. Section \ref{sec:Approach} presents our approach. Section \ref{sec:Experimental} presents the experimental results and Section \ref{sec:Conclusion} concludes the paper.

\section{Maintenance Decision-Making As a Supervised RL Problem}
\label{sec:Approach}

We define the problem of maintenance decision-making as learning a  maintenance policy that generates optimal maintenance actions such as  ``continuation of the operation" or ``the visitation of the repair shop" based on the observations and the previous maintenance actions. The observations can come in different formats. For example, they may include time-series sensor data such as temperature, and  event data such as  the check engine light.
A maintenance action is optimal if it leads to the maximum profit. The profit is the difference between the generated revenue and the  operating costs. The  operating costs include both maintenance  costs and  failure costs. Failure costs are often much higher than maintenance  costs because they occur without  planning. More importantly, in some cases, a failure can put the operator's safety in danger. In this paper, we  rely solely upon offline data to learn the maintenance policy. This makes our approach scalable and  practical as it does not require the development  and maintenance of a simulator.

Consider the case wherein offline data includes the observations, $o_k$, actions, $a_k$ and reward $r_k$ at each given time $k$. A trajectory $i$  includes the set of observations over time, the set of actions taken over time,  and the set of rewards gained over time. For a given window with the length of  $T$, we can extract the $T$ steps' observation history of the equipment at each given time, $k$, as:
\begin{equation}
O^T_k= \{o_{k-T+1},..., o_k\},
\label{eq:ol}
\end{equation}
using the equipment trajectory.
We can also extract $T$ steps' action history of the equipment at each given time, $k-1$, as:
\begin{equation}
A^T_{k-1}= \{a_{k-T},..., a_{k-1}\}.
\label{eq:a}
\end{equation}
Note that  unlike the observation history where we included $o_k$, here, we consider $a_{k-1}$ as the last action. This is because our goal is to use the  observation history and the action history to predict action  $a_k$.

Consider horizon $H$. We can also calculate future cumulative rewards for each trajectory from time $k$ up to the horizon, $H$,  as $r^H_k=\sum_{i=k}^{k+H} r_i$ using  offline data. Calculating $r^H$ over the past $T$ time steps, we have:
\begin{equation}
R^H_k= \{r^H_{k-T+1},..., r^H_{k}\}.
\label{eq:r}
\end{equation}

Next, we use  a supervised learning method to learn a model which predicts action $a_k$ using  $O^T_k$, $A^T_{k-1}$, $R^H_k$.
\begin{equation}
\hat{a}_k=model(O^T_k, A^T_{k-1}, R^H_k)
\label{eq:model}
\end{equation}
The idea here is to solve maintenance decision-making using an offline supervised  RL approach similar to Decision Transformer
 \cite{chen2021decision} and Trajectory Transformer \cite{janner2021reinforcement}. In the training stage, we learn a mapping from  future expected  cumulative rewards to the actions at each state using offline data. In the application, we feed the model with high future expected cumulative rewards and we expect the model to generate actions to achieve such.

\begin{figure}[t]
\centering
\includegraphics[scale=.34]{Model}
\caption{Proposed architecture for smart maintenance decision-making.}
\label{fig:model}
\end{figure}

To design a reliable maintenance decision-making  model, we have to consider the following points:
\begin{itemize}
\item \textit{RUL estimation and failure prediction models:}
when enough data is available, our proposed maintenance decision-making framework can learn an optimal policy from the data  without requiring  a separate RUL estimation   or failure prediction unit. However, it is very common for  industries to  develop reliable RUL estimation models that can estimate the remaining useful life of   equipment.
When these models are already available, their output can be used as an input to our decision-making model to improve accuracy and simplify the training process as it is shown in Figure \ref{fig:model}.  The formulation is:
 \begin{equation}
\hat{a}_k=model(O^T_k, A^T_{k-1}, R^H_k,RUL_k),
\label{eq:model2}
\end{equation}
 where $RUL_k$ presents the remaining useful life at time $k$. Note that during the training, we may have access to the actual RUL. However,   during the application process, we have to use the RUL estimation.

\item \textit{Reward function:}
 the goal of predictive maintenance is to minimize the operation costs and maximize the profit. Typically,  an equipment generates profit as it operates. So to maximize the profit, the industries aim to maximize equipment utilization. On the other hand, failures can be very expensive to fix or even dangerous for the operators. Therefore, operating  equipment to the point of failure is often not a viable option. The cost of maintenance and repair is another factor that we have to consider in maintenance decision-making. In addition to the extraneous labor and parts costs, unnecessary maintenance  reduces industries' profits by lowering  utilization time. Finally, the repair cost can depend upon the state of the equipment. For example, for some equipment, early fixes may be cheaper than late-stage repairs. Designing a reward function that can capture these complexities associated with costs and profits is crucial in designing an RL algorithm for maintenance decision-making. Fortunately, in many industries, the cost and profit of equipment operations are recorded in details and therefore can be used to design the reward function.

 \item \textit{Total expected reward:} as it is shown in equation (\ref{eq:model})  and (\ref{eq:model2}), the expected reward is an input to our model. During the training, we use the actual costs to compute the future rewards for each trajectory. During the application,   we must feed high expected future rewards so that the model generates good actions. Naturally, one may ask  ``is there a limitation on how high the  expected future rewards can be during the application?"  The answer to this question is  yes. Our experimental studies  in this paper  and previous experiments done by \cite{chen2021decision} show that after a certain point, increasing the expected future rewards has an opposite effect and degrades the performance. Our supervised offline RL algorithm uses historical data to learn the mapping from expected rewards to the actions at each state. Therefore, if during the application we feed the model with an expected future reward significantly outside of the training  dataset distribution, the model may generate unreliable  actions.

 \citeA{chen2021decision} show that the total expected rewards
and the actual total returns are highly correlated as long as we feed the model with the total expected reward in the range which was observed in the offline data.  More interestingly, they show that
it is even  possible to achieve  higher returns than the maximum episodic return
in the dataset for certain specific tasks. This shows that offline supervised RL can   extrapolate accurately to some extent.
In this paper, we also observe  a high correlation between the total expected reward
and the  observed total reward for the maintenance problem.  Moreover, our experiments  demonstrate  that selecting the  total expected reward    within     the vicinity of the maximum episodic return  available
in the offline dataset is a good practical approach to set this parameter  in the maintenance application.

\item \textit{Architecture:}
\citeA{chen2021decision}  and  \citeA{janner2021reinforcement}  used the transformer architecture \cite{vaswani2017attention} to learn the model that predicts the maintenance action (see  equations (\ref{eq:model}) and (\ref{eq:model2})). However, their approach can be implemented  using any  classifier  model when we have discrete action space. Similarly,   any regression model can be used to predict actions when we have continuous  action space. The  main parameters we have to consider in selecting the model  architecture are: 1) the amount of available data,   2) the complexity of the data, and 3) the available computational resources. In the experimental section, we will demonstrate  that a simple fully connected neural network architecture  can achieve acceptable performance for our dataset.
 \end{itemize}

\section{Case Study}
\label{sec:Experimental}

In this section, we use the C-MAPSS dataset to show the application of our proposed maintenance decision-making algorithm.  This dataset includes time-series sensor measurements of jet engine under  different operational conditions and fault modes. The dataset is generated by   the Prognostic Center of Excellence at NASA Ames using the C-MAPSS simulator \cite{saxena2008turbofan}.

We focus on equipment F002 in this experimental study. The training dataset for this equipment includes 260 trajectories. Each trajectory ends when the equipment fails. The test dataset includes 259 trajectories. Unlike the training dataset, the trajectories in the test dataset do not include failures. This was  mainly designed  to hide the failure points from the participants in the data challenge. To show the performance of our algorithm, when the system fails, we consider the first 250 trajectories in the original training dataset as our new  training dataset and the last 10 trajectories as our first test dataset. We consider the original test dataset as the second test dataset in this paper.

\subsection{Removing the effect of operation modes}

Similar to  previous work \cite{wang2019remaining}, we first normalize sensor variables with respect to  the operation modes. The goal here is to remove the effect of different operating modes from the sensor data so that we can learn a maintenance decision-making model that works for all operating modes. Equations (\ref{eq:norm}) and (\ref{eq:norm2}) represent the normalization process.
 \begin{equation}
\hat{s}_i =\text{model}_i(\text{operation modes})
\label{eq:norm}
\end{equation}
where $\text{model}_i$ represents the regression normalization model for sensor variable $s_i$ and $\hat{s}_i$ represents the estimation of $s_i$ made by its normalization model  using three operating modes as the inputs. Note that the equipment has  21 sensor variables and we train a regression normalization model for each sensor.  We normalize each sensor variable $s_i$ using  its estimation $\hat{s}_i$ as follows:
 \begin{equation}
s_{in}=\frac{s_i}{\hat{s}_i}
\label{eq:norm2}
\end{equation}
where $s_{in}$ represents normalized $s_i$.

\begin{figure}[t]
\centering
\includegraphics[scale=.34]{sensor}
\caption{An original sensor value, $s_i$, and its normalized value, $s_{in}$, show that normalization amplifies  the effect of degradation on sensor variables by removing  the effect of operation modes.}
\label{fig:sensor}
\end{figure}
Figure \ref{fig:sensor} shows a sensor variable  before and after normalization. We can see   the degradation process in the sensor variables when the effect of operation mode is removed.

\subsection{Generating a  dataset for training   decision-making model}

Unfortunately, the C-MAPSS dataset does not include any maintenance action. In the original training dataset, during the simulation, all of the equipment has been utilized to the point of failure. In the test dataset, the simulation has stopped before the failure without taking any  maintenance action.
 To train  our models (see equations (\ref{eq:model}) and (\ref{eq:model2}) ), we must have access to historical actions. We modify the dataset by assuming two kinds of actions: 1) continue operation; 2) repair the equipment. The original dataset only includes continual operation (no actions). We consider random repairs in the dataset. We assume that after each repair, the equipment will return to its initial condition.

\subsection{Costs and rewards}

In this paper, we consider a  simple cost structure.   In our setting, each failure costs $250 + U(-50,50)$; each repair costs $25+  U(-5,5)$; and  equipment generates $1 +  U(-0.2,0.2) $ profit per operating cycle, where   $U(a,b)$ represents uniform distribution between $a$  and $b$.
 Note that the ratio between the cost of failure and repair defines how conservative the network will be in taking proactive maintenance action. In an extreme case where the repair cost is equal to or more than the cost of failure, the optimal solution would be to operate the equipment to the point of failure. In another extreme case where  repair costs nothing, the optimal solution would be to repair the equipment after each cycle.
 In practice, these numbers should be derived from facts on the ground. It is also possible that the costs change based on operation modes and the state of the equipment. For example, for a truck, it may be cheaper and faster to do repair and maintenance in a large city compared to a remote small town where  expert operators and equipment parts may be scarce.

\subsection{RUL estimation}
\begin{figure*}[t]
\centering
\includegraphics[scale=.55]{RUL_r_2}
\caption{RUL vs RUL estimation for our first test dataset (unit number 251 to 260). $R^2 = 0.53$.}
\label{fig:RUL}
\end{figure*}
As we mentioned earlier, we divided the equipment F002 dataset to train and test. We considered unit numbers  one to 250 as our new training dataset, and  unit numbers 251 to 260 as our first  test dataset. We considered the original equipment test dataset as our second test dataset.
We then used the training dataset to learn a normalization model for each sensor. We use the long short-term memory (LSTM) architecture proposed by  \cite{zheng2017long} for the RUL estimation
 using normalized sensors.

  Figure \ref{fig:RUL} shows the RUL estimation for our first test dataset. Figure \ref{fig:RUL_LSTM} shows the RUL prediction  for the first 6 trajectories in our second test dataset. We  see that the RUL estimation becomes more accurate as it gets closer to the equipment's end of life.
\begin{figure*}[t]
\centering
\includegraphics[scale=.55]{RUL_LSTM}
\caption{RUL vs RUL estimation for second test dataset using LSTM neural networks. $R^2 = 0.42$}
\label{fig:RUL_LSTM}
\end{figure*}
Our first test  dataset (trajectories 251 to 260 of the original training dataset) has a more similar distribution to our training dataset (trajectories one to 250 of the original training dataset) and therefore our model predicts RUL more accurately for this dataset compared to the second test dataset ($R^2 = 0.53$ vs $R^2 = 0.42$).

\subsection{Offline RL}

In this subsection, we train a   multilayer perceptron (MLP)  neural network  to predict one of these two actions: 1) continue operation or 2) repair the equipment, based on the expected future rewards. Our network has one hidden layer with 100 neurons. We use the rectified linear activation function (ReLU),  we apply the Adam optimization and set the  learning rate  equal to $0.001$.
We train two models: in the first model, we only use sensor data and the  expected rewards to predict the action. In the second model, we add RUL as an additional feature. For the model with RUL, we use actual RUL during the training and the estimation of RUL during the test time.

\subsection{Application}

We consider each experiment in the test dataset and episode. An episode  ends when 1) its trajectory ends, 2) the equipment  fails, or 3) the equipment goes for  repair. Each cycle of operation brings 1 unit of profit, and each repair costs 25 units. An early repair  lowers the profits by ending the episode. A failure costs 250 units. We assume that  repair decisions should be made at least 10 cycles before  failure occurs. This is equal to the time that the equipment requires to get to the repair shop or to stop  operation safely.   Figure \ref{fig:results1} shows that no action would cost -135.7 on average (see the purple line) for our first  test dataset. This is because the first test dataset includes failures and each episode without action would end with a failure.
Obviously, allowing the equipment to operate  to the point of failure    is not an acceptable solution for any business.
\begin{figure}[t]
\centering
\includegraphics[scale=.28]{dataset_1.pdf}
\caption{Average cumulative rewards as a function of expected reward for the first test dataset. The maximum average cumulative rewards for the  model with RUL achieves at the total expected reward =  79, and the maximum average cumulative rewards for the model without RUL achieves at the total expected reward =  100.5.}
\label{fig:results1}
\end{figure}
Figure \ref{fig:TestResults} shows that no action would earn  96 units of profit  on average in the second test dataset (see the purple line). This is because in the second dataset, the episode often ends before reaching its  failure point.

The optimal solution is to take the equipment for repair
 10 cycles before reaching  failure. In this case, we maximize  profit as we utilize the equipment to its full potential while  avoiding the high cost of failure. The blue line in Figure \ref{fig:results1} shows this scenario. We achieved 88.3 average cumulative reward using this approach for the first test dataset.   Of course, this is not a practical solution because it requires the knowledge of exact RUL. An alternative option would be to use  the RUL estimation instead  of the actual RUL. The orange line in Figure \ref{fig:results1} shows this scenario for the first test dataset. We achieved 82.7 average cumulative reward using this approach.

\begin{figure}[t]
\centering
\includegraphics[scale=.28]{dataset_2.pdf}
\caption{Average cumulative rewards as a function of expected reward for the second test dataset.The maximum average cumulative rewards for the  model with RUL achieves at the total expected reward =  75, and the maximum average cumulative rewards for the model without RUL achieves at the total expected reward =  96.}
\label{fig:TestResults}
\end{figure}

 Note that this approach is not fully data-driven as it requires the decision-making  logic of taking the equipment for repair 10 cycles before  failure. This logic is trivial for this simple example. However, when we have more complicated cases, it may not be that simple to come up with an optimal solution. For example, the repair cost may change based on the system  conditions. This may make our simple logic suboptimal as the repair 10 cycles  before  failure can be more expensive than earlier repairs.  It is also possible that we  have more than one part to repair. For these cases, offline RL presents a better alternative as it learns the optimal policy from data.

 We can see in Figure \ref{fig:TestResults} that the difference between the perfect action and the action based on RUL estimation is wider in the second dataset compared to the first test dataset. This is because the RUL estimation in the second dataset is  less accurate and thus makes decisions based on RUL estimation less reliable.

 Figure \ref{fig:results1} and Figure \ref{fig:TestResults} show the offline RL performance for two scenarios: 1) in the first one, we used real RUL during the training, and RUL estimation during the test as an input to our offline RL model, 2) in the second case, we trained the offline RL purely relying on sensor data.  The results show that the model achieves similar performance with or without  RUL estimation.
  We believe that this is because in this dataset, the normalized sensors capture the degradation process  and having access to RUL estimation does not provide significant additional information for the model.

Note that the  offline RL outperforms decisions  based on RUL estimation in the second test dataset (see Figure \ref{fig:TestResults}). Figure \ref{fig:TestResults} shows that the performance of the model with the RUL estimation peaked  at total requested rewards equal to 75 and the performance of the model without RUL estimation peaked at total requested rewards equal to 96. However, both models perform  better than the  decision based on RUL estimation for a fairly wide range of total requested rewards.
This presents a huge potential for offline RL in maintenance decision-making, especially for real-life problems when estimating RUL is not a trivial task.

\section{Conclusion}
\label{sec:Conclusion}
In this paper, we presented a framework to use offline reinforcement learning for maintenance decision-making. The results show  that offline RL can provide  decision-making  competitive  to using RUL estimation. Moreover, offline RL can generate acceptable solutions even  without requiring learning the RUL estimation model. Our approach provides a framework that can be applied to more complicated maintenance decision-making challenges.

\begin{singlespace}
\bibliographystyle{apacite}


\section{Introduction}

Simple random walk SRW is the cumulative sum of i.i.d. $\pm 1$ increments, with some initial integer value. Reflected simple random walk RSRW is non-negative, with increments like SRW whenever positive, and increment $1$ when zero. The local time $L$ at zero up to time $n$ of either process is the number of its visits to zero up to time $n$. A formula will be derived for the joint distribution of the value of the process at time $n$ and the local time at zero up to time $n$, for SRW and RSRW. Corresponding formulas in continuous time exist for Brownian Motion (\cite{BorodinSalminen},\cite{JainPruitt}).

\bigskip

The motivation for this study is an M/M/$c$/$n$ service system observed when all servers are busy. The reflected random walk is the waiting room occupancy and local time is the number of arriving customers that fill up the queue to capacity. The joint distribution of RSRW and $L$ determine the distribution of the number of waiting customers an arriving customer will meet until entering into service. At the time of a pandemic, this distribution may provide a tool to determine a safe maximal queue size. This is a companion paper of Hassin, Meilijson and Perlman \cite{HMP} that considers Nash equilibrium, safe and socially optimal maximal queue sizes, under a utility function that take into account the negative network effects of waiting time, total time spent with waiting customers and the number of these customers, the subject matter that motivates the current study.

\bigskip

Formally, let $\epsilon_t \ ; \ t=1,2,\dots$ be i.i.d. with $P(\epsilon_t=1)=1-P(\epsilon_t=-1)=q \ , \ q \in (0,1)$, be the increments of simple random walk (SRW) $S_T=s+\sum_{t=1}^T \epsilon_t$ starting at $S_0=s$. Reflected simple random walk RSRW $R$ is defined as $R_T=s+\sum_{t=1}^T \delta_t$ (starting at $R_0=s \ge 0$), where $\delta_t=\epsilon_t$ if $R_{t-1}>0$ and $\delta_t=1$ otherwise (i.e., if $R_{t-1}=0$). The trajectories of $R$ are non-negative, increasing to $1$ whenever at zero and behaving as $S$ otherwise. $R$ can be identified with $|S|$ if $q={1 \over 2}$, but not otherwise. Let $L_0=0$ and for $t>0, L_t = L_{t-1} + I_{\{R_t=0\}}$ be the {\em local time} of RSRW $R$ at zero, the cumulative number of visits to zero in positive time.

\section{Reflected simple random walk}

Let $NT(T,s,r,l)$ be the number of $R$-trajectories of length $T \ge 1$ with $R_0=s \ge 0$, $R_T=r \ge 0$ and $L_T=l$. The ease of identifying the binomial distribution of $S_T$ (w.l.o.g $s=0$) stems from two distinct elements: All trajectories with $S_T=r$ have the same probability $q^{{T+r} \over 2} (1-q)^{{T-r} \over 2}$ and the number of such trajectories is straightforward to identify as $\binom{T}{{{T+r} \over 2}}$. This is not the case for the distribution of $R_T$, but as the next theorem shows, it is for the joint distribution of $R_T$ and $L_T$.

\begin{theorem} \label{localtimereflected}
Let $T \ge 1$ and $s, r, l$ be non-negative integers. Let $R$ be reflected simple random walk whose increments are positive with probability $q \in (0,1)$. Then

\begin{equation} \label{jointRL}
P(R_T=r, L_T=l | R_0 = s)  =  NT(T,s,r,l)q^{{{T+r-s} \over 2}-l+I_{\{r=0\}}-I_{\{s=0\}}}(1-q)^{{{T-r+s} \over 2}} \\
\end{equation}
where the number of trajectories $NT(T,s,r,l)$ is zero except for the feasible cases: $s$ any non-negative integer, $r \in \{s-T, s-T+2, s-T+4 \dots, s+T\} \cap [0,s+T], 0 \le l \le max(0,{{T-s-r} \over 2}+I_{\{s > 0\}})$. For these feasible cases the coefficients $MT(T,s,r,l)$ are as follows.

\begin{eqnarray}
NT(T,0,r,l) & = & \binom{T-l-1}{{{T+r} \over 2}-1}-\binom{T-l-1}{{{T-r} \over 2}-l-1} \nonumber \\
NT(T,s>0,r,0)& = & \binom{T}{{{T+r-s} \over 2}}-\binom{T}{{{T-r-s} \over 2}} \nonumber \\
NT(T,s>0,r,l>0) & = & \binom{T-l}{{{T-r-s} \over 2}-l+1}-\binom{T-l}{{{T-r-s} \over 2}-l} \label{jointRLtraj}
\end{eqnarray}
\end{theorem}

\begin{proof}
For initial state $s=0$ the first increment is $+1$ a.s. Hence, the probability of arriving from $s=0$ to $(r,l)$ in time $T$ is the same as the probability of arriving from $s=1$ to $(r,l)$ in time $T-1$. It is straightforward to check that the top formula in (\ref{jointRL}) derives accordingly from any of the other two, that are the same formula for $s=1$ (but not for $s>1$).

The single-trajectory probabilities are immediate to ascertain, as is the fact that all trajectories with the same triple $(s,r,l)$ are equally likely. What remains is the combinatorial part.

For $s>0, r>0, l=0$, $NT(T,s,r,0)$ is the number of non-negative trajectories of length $T$ starting at $s$ and ending at $r>0$ without ever visiting zero. Simple random walk has $\binom{T}{{{T+r-s} \over 2}}$ as many trajectories, from which those that visit zero have to be subtracted. By the usual reflection principle, the latter are just as many as the number of SRW paths from $s$ to $-r$, totalling $\binom{T}{{{T-r-s} \over 2}}$. Hence,
$NT(T,s,r,0)=\binom{T}{{{T+r-s} \over 2}}-\binom{T}{{{T-r-s} \over 2}}$, as claimed.

For $s>0,l>0$, the set to be counted (with $NT(T,s,r,l)$ elements) is the set $A(T,s,r,l)$ of non-negative trajectories of length $T$ with ${{T+r-s} \over 2}$ positive increments and $l$ visits to zero in positive time. Let $l>0$. For every element in $A(T,s,r,l)$ there is a well defined time $t$, the time of last visit to zero. Remove the increment $\delta_t=-1$ from the string of increments to obtain a shorter string $(\delta_1, \delta_2,\cdots,\delta_{t-1},\delta_{t+1},\cdots,\delta_T)$ with partial sums $s=R_0, R_1, \cdots,R_{T-1}$. This shorter RW trajectory $R$ has length $T-1$, the same number ${{T+r-s} \over 2}$ of positive increments, and last visit to $1$ at time $t-1$. Furthermore, this mapping from the set of $A(T,s,r,l)$ with last visit to zero at time $t$ to the set of $A(T-1,s,r+1,l-1)$ with last visit to $1$ at time $t-1$ is one-to-one. As long as $l-1>0$,  $t$ partitions $A(T,s,r,l)$ (differently but) just as it partitions $A(T,s,r+1,l-1)$, so this mapping is onto and these two sets have the same number of elements. Since the number of trajectories is  invariant when $T$ and $l$ are decreased by $1$ while at the same time increasing $r$ by $1$. Hence, the answer is the same as $NT(T-l+1,s,r+l-1,1)$. The plan is to carry out one more transition from $l=1$ to $l=0$. The mapping described above to $NT(T-l,s,r+l,0)$ is correct only if $s=1$ (when the existence of a visit to $1$ is guaranteed) and needs a different treatment for $s>1$, in which case the trajectories counted by $NT(T-l,s,r+l,0)$ include those that don't visit $1$ at all, whereas the one-to-one-onto mapping was only to those that visit $1$. The trajectories going from $s>1$ to $r+l$ in time $T-1$ without visiting $1$ are the same as those going from $s-1$ to $r+l-1$ in time $T-l$ without visiting $0$, i.e., $NT(T-l,s-1,r+l-1,0)$ in number.

Summarizing and performing some book-keeping,

\bigskip

$NT(T,s,r,l) = NT(T-l,s,r+l,0)=\binom{T-l}{{{T+r-1} \over 2}}-\binom{T-l}{{{T+r+1} \over 2}}$ for the case $s=1$, while for $s>1$:
\begin{eqnarray}
 & & NT(T-l,s,r+l,0) - NT(T-l,s-1,r+l-1,0) \nonumber \\
 & = & \Big[\binom{T-l}{{{T-l+r+l-s} \over 2}}-\binom{T-l}{{{T-l-r-l-s} \over 2}}\Big]-\Big[\binom{T-l}{{{T-l+r+l-1-(s-1)} \over 2}}-\binom{T-l}{{{T-l-r-l+1-(s-1)} \over 2}}\Big] \nonumber \\
 & = & \Big[\binom{T-l}{{{T+r-s} \over 2}}-\binom{T-l}{{{T-r-s} \over 2}-l}\Big]-\Big[\binom{T-l}{{{T+r-s} \over 2}}-\binom{T-l}{{{T-r-s} \over 2}-l+1}\Big] \nonumber \\
 & = & \binom{T-l}{{{T-r-s} \over 2}-l+1}-\binom{T-l}{{{T-r-s} \over 2}-l} \nonumber
\end{eqnarray}
\end{proof}

The marginal distribution of $R_T$ can be alternatively obtained by Markovian matrix methods. Let $Q$ be the transition matrix of RSRW, with all entries zero except $Q(0,1)=1$ and $Q(i,i+1)=1-Q(i,i-1)=q$ for $i \ge 1$. Now let $Q_1=Q(0:N,0:N)$ for any $N$ exceeding $T+s$. Then the marginal distribution of $R_T$ is row $s$ of $Q_1^T$.

Expected local time at zero up to time $T$, starting at $s$, can be obtained as follows. Since it is the sum over n (from $0$ to $T$) of $Q^n(s,0)$, it is the $(s,0)$ entry of $(I-Q_1)^{-1}(I-Q_1^{T+1})$. But expectations of nonlinear functions of local time, such as its distribution, seem to be outside the standard range of algebraic methods.

\begin{figure}[hbt!]
\begin{center}
{\includegraphics[scale=0.40]{T80S3Q40.pdf}}
\qquad
{\includegraphics[scale=0.40]{T80S3Q45.pdf}}
\qquad
{\includegraphics[scale=0.40]{T80S3Q50.pdf}}
\qquad
{\includegraphics[scale=0.40]{T80S3Q55.pdf}}
\end{center}
\caption{Joint distribution of value and local time at zero for reflected random walk starting at $s=3$ and ending at time $T=80$. From top left to bottom right, $q=0.40, 0.45, 0.50, 0.55$, displaying a transition from downdrift and large local time at zero to diffusion with positive drift and few returns to zero. Value has been displayed up to $35$ and local time at zero up to $25$, switching coordinates from top to bottom.}
\label{fig:jointdist}
\end{figure}

\begin{figure}[hbt!]
\begin{center}
{\includegraphics[scale=0.40]{expectedvalue.pdf}}
\qquad
{\includegraphics[scale=0.40]{expectedlocaltime.pdf}}
\end{center}
\caption{Joint distribution of value and local time at zero for reflected random walk starting at $s=3$ and ending at time $T=80$, for $q=0.40, 0.45, 0.50, 0.55$. Left, conditional expectation of local time given value ($q=0.40$ is highest). Right, conditional expectation value given local time ($q=0.40$ is lowest).}
\label{fig:condexp}
\end{figure}

\section{Simple random walk}

Rather than counting trajectories directly, their statistics will be derived from those for reflected simple random walk. In essence, each {\em solid bridge} of RSRW gives rise to two solid bridges of SRW, where a solid bridge is a sub-trajectory starting and ending at zero without zero values in between.

Adhering to the notation above, $S$ stands for SRW and $R$ stands for RSRW. Unlike RSRW where the probability of a trajectory from $s$ to $r$ depends on $l$, for SRW or simple random walk bridge $l$ can determine whether the trajectory is possible or impossible, but no internal information on the trajectory can influence the positive value of the probability because, in statistical terms, the sum of the increments is a sufficient statistic for $q$. Each trajectory from $s$ to $r$ in time $T$ has probability $q^{{{T+r-s} \over 2}}(1-q)^{{{T-r+s} \over 2}}$ and the only question is combinatorial, how to split the number $\binom{T}{{{T+r-s} \over 2}}$ of trajectories into the various subsets determined by $l$. Taking care of special cases first, the number of trajectories of length $T$ from $s \neq 0$ to $r$ with $(l=0,r \neq 0)$ or $(l=1,r = 0)$ is the same as the corresponding number of trajectories for RSRW from $|s|$ to $|r|$. The number of trajectories of length $T$ with $l=1$ from $s \neq 0$ to $r \neq 0$ is the same as the number of reflected random walk trajectories with $l=1$ from $|s|$ to $|r|$. For all other cases, $l'=l+I_{\{s=0\}}$ (augmented $l$) is at least $2$ and the trajectory has $l'-1$ solid bridges. Hence, the number of SRW trajectories is the corresponding number for RSRW (from $|s|$ to $|r|$) multiplied by $2^{l'}$.

Summarizing,

\begin{theorem} \label{localtimesimple}
Let $T \ge 2$ and $s, r, l$ be integers (with $l \ge 0$). Let $S$ be simple random walk whose increments are positive with probability $q \in (0,1)$. Then
\begin{equation} \label{jointSL}
P(S_T=r, L_T=l | S_0 = s) = MT(T,s,r,l)q^{{{T+r-s} \over 2}}(1-q)^{{{T-r+s} \over 2}}
\end{equation}
where the number of trajectories $MT(T,s,r,l)$ is zero except for the feasible cases: $s$ any integer, $r \in \{s-T, s-T+2, s-T+4, \dots, s+T\}, 0 \le l \le max(0,{{T-|s|-|r|} \over 2}+I_{\{s \neq 0\}})$. For these feasible cases the coefficients $MT(T,s,r,l)$ are expressed in terms of $l'=l+I_{\{s=0\}}$ as

\begin{equation}
\label{jointRLtrajsimple}
MT(T,s,r,l) = \begin{cases}
\begin{aligned}
    2^l \bigg [\binom{T-l-1}{{{T+|r|} \over 2}-1}-\binom{T-l-1}{{{T+|r|} \over 2}}\bigg] & ~s=0 \\
    \binom{T}{{{T+|r|-|s|} \over 2}}-\binom{T}{{{T-|r|-|s|} \over 2}} & ~s \neq 0, l=0 \\
    2^{l-1} \bigg [\binom{T-l}{{{T+|r|+|s|} \over 2}-1}-\binom{T-l}{{{T+|r|+|s|} \over 2}}\bigg] & ~s \neq 0, l>0
\end{aligned}
\end{cases}
\end{equation}

\end{theorem}

\bigskip

\section{Mixture distributions}

As expressed by (\ref{jointRLtraj}) and (\ref{jointRLtrajsimple}), and illustrated in Figures \ref{fig:jointdist} and \ref{fig:bothcases}, each of the two joint distributions (\ref{jointRL}) and (\ref{jointSL}) of value at time $T$ and local time at zero up to time $T$, starting at $s>0$, is a mixture of two (discrete but) smooth distributions, one modelling paths with no visits to zero ($L=0$, middle term in (\ref{jointRLtraj})) and the other covering all paths that visit zero ($L>0$, last term in (\ref{jointRLtraj})). Figure \ref{fig:bothcases} displays RSRW (left) and SRW (right) for $T=160, s=5,q=0.45$.

\bigskip

{\bf The univariate mixture component $\{L=0\}$ is log-concave in $r$}. This distribution is analyzed in detail for even $r$. The case of odd $r$ only requires a minor adjustment in the auxiliary terms $W$ and $x$ to be defined. Consider the logarithm of the probability (\ref{jointRL}) as a function of $r$, for fixed $T, q, l=0$ and $s>0$. It is the sum of $\log(NT)$ and a linear function of $r$. Hence, if $\log(NT)$ is a concave function of $r$, the $\{L=0\}$-component of (\ref{jointRL}) is log-concave for each $q$. To show concavity of $\log(NT)$ means to show that the increment $\log(NT(T,s,r+2,0))-\log(NT(T,s,r,0))$ is a decreasing function of $r$. In other words, that ${{NT(T,s,r+2,0)} \over {NT(T,s,r,0)}} \le {{NT(T,s,r,0)} \over {NT(T,s,r-2,0)}}$. The proof will be finished by expressing the cross-product $NT(T,s,r,0)^2-NT(T,s,r-2,0)NT(T,s,r+2,0)$ as the sum of three non-negative terms, where $W={{T-s} \over 2}<{{T} \over 2}$ and $x={r \over 2}$. These are formal equalities. Whenever the RHS is not well defined (non-positive terms in the denominator), the LHS is non-negative, as needed.

\begin{eqnarray}
\binom{T}{W+x}^2-\binom{T}{W+x-1} \binom{T}{W+x+1}& = &\binom{T}{W+x}^2 {{T+1} \over {(W+x+1)(T-W-x+1)}} \nonumber \\
\binom{T}{W-x}^2-\binom{T}{W-x-1} \binom{T}{W-x+1}& = &\binom{T}{W-x}^2 {{T+1} \over {(W-x+1)(T-W+x+1)}} \nonumber
\end{eqnarray}
\begin{eqnarray}
 & & \binom{T}{W-x+1} \binom{T}{W+x+1}+\binom{T}{W+x-1} \binom{T}{W-x-1}-2\binom{T}{W+x} \binom{T}{W-x} \nonumber \\
 & = & (T+1)(T-1-2 W)[{1 \over {(W+1)^2-x^2}}+{1 \over {(T-W+1)^2-x^2}}]  \label{sausage}
\end{eqnarray}

Log concavity of the bivariate mixture component $\{L>0\}$ in $(r,l)$ is under study, and the report will be updated accordingly. Concavity on a finite subset of $R_2$ requires special treatment.Kannai \cite{kannai} is a good start.

\begin{figure}[hbt!]
\begin{center}
{\includegraphics[scale=0.40]{RSRWpair.pdf}}
\qquad
{\includegraphics[scale=0.40]{SRWpair.pdf}}
\end{center}
\caption{Joint distribution of value and local time at zero for reflected (left) and simple (right) random walk starting at $s=5$ and ending at time $T=160$, for $q=0.45$. Both graphs display mixtures of two distributions, one with zero local time and the other with visits to zero. The value scale of SRW has zero at coordinate 60.}
\label{fig:bothcases}
\end{figure}

\section*{Funding} This research was supported by the Israel Science Foundation, grant No. 1898/21.

\bibliographystyle{plain}