\section{Introduction}
\label{sec_1}

Evidence-based program evaluation is increasingly becoming an essential tool for policymakers to fine--tune socio--economic policies (Cerulli, 2015). Generally, decision-makers run multiple rounds of similar policies at national, regional, or local level. As recurrent policies take place, they can collect massive evidence on past programs' operation and effects that increasingly enables a learning process based on past experience.

Program evaluation draws on three dimensions: ex--ante, in--itinere, and ex--post. Abstracting from in--itinere evaluation, this article focuses on how decision-makers can improve ex-ante evaluation of policy effects based on the results obtained from the (antecedent) ex-post evaluations. In tune with the recent literature, I call this process \textit{policy learning} (Athey and Wager, 2019).

How does this process take place? And more specifically, what do we mean by ex - ante program evaluation in this specific context? Answering these questions requires delving into a deeper understanding of the causal dynamic underlying a policymaker's action.

To fix the idea, consider a policymaker administrating a policy $P$ targeting a specific outcome $Y$. Either for ethical or policy constrains reasons, the policymaker hardly wants to (or can) pick up beneficiaries at random. She generally adopts selection strategies aimed at optimizing her objective function, typically by increasing the overall policy effect as much as possible. The point in question is: which effect?

It is well-known that the total effect (let's call it $\gamma$) of a policy intervention can be decomposed into two sub-effects, a direct one (let's call it $\alpha$), and an indirect one (let's call it $\beta$). The direct is the effect of the policy ``as if'' the decision--maker had selected both beneficiaries and non--beneficiaries at random. The indirect effect is that part of the total effect due to the selection process operated by the policymaker's choice of beneficiaries (the so--called \textit {policy assignment rule}). In this setting, we thus have:

\begin{equation}
\gamma(S) = \alpha + \beta(S)
\end{equation}
where the indirect effect is clearly a function of the selection process $S$ operated by the policymaker. Different policy assignment rules affect the total effect (or \textit{welfare}) of the policy, thus raising a quest for an ``optimal'' policy rule within certain classes of feasible rules (Manski, 2004).

Ex--post program evaluation techniques, either experimental (randomized control trials, RCT) or quasi--experimental (re-balancing procedures), generally aim at consistently estimating $\alpha$ by assuming $\beta$ to be a bias induced by the selection process. This is the econometrician perspective, that does not coincide however with the policymaker viewpoint. Given budget, ethical, or contextual constraints, the policymaker is in fact mainly interested in increasing the total effect of the policy, and thus finding the selection rule $S^{*}$ which maximizes this effect (Dehejia, 2005).

Recently, the literature has produced significant theoretical advances on optimal policy assignment, both in statistics, biostatistics and econometrics\footnote{Specifically, this paper refers to the econometric papers of Manski (2004) and Kitagawa and Tetenov (2012) who extended the former. In his seminal paper, Manski (2004) assumes different probability distributions of the possible experimental outcomes, and considers treatment rules assigning individuals to those treatments producing the best experimental outcomes. These ``conditional empirical success'' rules consider different subsets of the observed covariates space. The author derives bounds for the welfare regret associated to such rules showing that it converges to zero at a rate of at least $0.5$. Kitagawa and Tetenov (2012) extend Manski approach to a broader class  of rules, including threshold-based ones, and show that similar rate of convergence are achieved by these broader class of rules. Other contributions to the literature are: Dehejia, 2005; Hirano and Porter, 2009; Athey and Wager, 2019; Bhattacharya and Dupas, 2012; Zhao et al., 2012; Zhang et al., 2012.}. Following in the footsteps of this literature, this paper wants to contribute by stressing the policymaker perspective via a practical implementation of an optimal policy assignment problem (or, equivalently, of an ``empirical welfare maximizing'' problem). More specifically, by focusing on the class of \textit{threshold--based} policies, we first set up the theoretical underpinnings of the policymaker selection problem, to then offer a practical solution to this problem via an empirical illustration using the popular LaLonde (1986) training program dataset. Moreover, the paper proposes an implementation protocol of the optimal solution that can be easily programmed with standard statistical software\footnote{With regard to the application presented in this paper, I developed a Stata 16 program for both a univariate and a bivariate empirical welfare optimization threshold--rule. All codes are available upon request.}.

The structure of the paper is as follows. Section \ref{sec_2} presents the theoretical underpinnings of the optimal constrained treatment assignment. Section \ref{sec_3} proposes a protocol for carrying out optimal constrained treatment assignment. Section \ref{sec_4} illustrates an empirical application using the popular LaLonde (1986) job training dataset and discusses the results. Section \ref{sec_5} concludes the paper.

\section{Optimal constrained treatment assignment}
\label{sec_2}

Let $X$ be an individual's vector of characteristics, $Y$ an outcome of interest, $T=\{0,1\}$ a binary treatment.
A policy assignment rule $\mathcal{G}$ is a function mapping $X$ to $T$, specifying which individuals are or are not to be treated:

\[\mathcal{G}: X  \rightarrow T \]
Define the policy conditional average treatment effect as:

\begin{equation*}
\tau(X)=E(Y_{1}|X) - E(Y_{0}|X)
\end{equation*}
where $Y_{1}$ and $Y_{0}$ represent the two potential outcomes of the policy, and $E_{X}[\tau(X)]=\tau$ the average treatment effect.
The policy actual total effect (or \textit{welfare}) $W$ is defined as:

\begin{equation*}
W = \sum_{i=1}^{N}T_{i}\cdot\tau(X_{i})
\end{equation*}
and the policy \textit{unconstrained} optimal total effect (or \emph{unconstrained maximum welfare}) as:
\begin{equation*}
W^{*} = \sum_{i=1}^{N}T^{*}_{i}\cdot\tau(X_{i})
\end{equation*}
where:

\begin{equation*}
T^{*}_{i}=\mathbf{1}[\tau(X_{i})>0]
\end{equation*}
is the optimal unconstrained policy assignment.

The difference between the (unconstrained) maximum achievable welfare and the
the welfare associated to the policy actually run is called \emph{regret}, and it is defined as:

\begin{equation*}
Regret = W^{*} - W
\end{equation*}
The regret is generally positive as decision makers are unable to cherry-pick the best treatment assignment possible, given their constrains and available information.
This is particularly true for randomized control trials (RCT) as treatment is, in this case, deliberately void of any selection objectives.

Table \ref{tab:tab1} sets out an illustrative example of the previous setting with three treated and three untreated units. The conditional average treatment effect $\tau(X)$ is reported in the third column of the table, while the last row of the fourth column shows the actual policy welfare that in this numerical example is equal to $10$. By definition, the optimal (unconstrained) treatment assignment would select only those units obtaining a positive $\tau(X)$, that in this specific case are the units 1, 3, 4, and 6. This is the optimal selection assignment $T^{*}$, with the last row of the last column showing that the largest welfare achievable by the policymaker is equal to $26$. The regret is therefore equal to $26-10=16$.

Because of eligibility, budget, ethical, or institutional constrains, policymakers are in general unable to implement the optimal unconstrained policy assignment,
being therefore obliged to rely on a \emph{constrained} assignment $T'$  which selects treated units according to (some) of their characteristics.
The welfare thus obtained, $W'$, may or may not be lower than $W^{*}$. However, one can ask whether it is possible for policymakers to produce the largest possible constrained welfare.
As the set of feasible treatment assignment rules is huge, replying to the previous question requires to restrict the focus on certain \textit{classes} of policies and find, among them, the optimal assignment to treatment (\textit{optimal constrained policy assignment}).

Although there are several classes of policies, policymakers often use few of them, such as ``threshold-based'', ``linear--combination'', or ``fixed-depth decision trees'' (Kitagawa and Tetenov, 2018).

Threshold-based assignment policies are popular in policy programs as they are simple to manage and easy to interpret. For this class of policies, policymakers select one or more than one variable of outstanding importance (such as, for example, firm number of employees, bank total asset, individual age, etc.) by selecting specific values (the thresholds) discriminating between treated and untreated units. To clarify, suppose to have just one single selection variable $x$, with threshold $c$. The assignment to treatment is clearly a function of $c$:

\begin{equation*}
T_{i}(c)=T^{*}_{i} \cdot \mathbf{1}[x>=c]
\end{equation*}
with corresponding welfare:
\begin{equation*}
W(c) = \sum_{i=1}^{N}T_{i}(c)\cdot\tau(X_{i})
\end{equation*}
We define the optimal choice of the threshold $c$ as the one maximizing $W(c)$ over $c$:

\begin{equation*}
c^{*}=\mathtt{argmax}_{c}[W(c)]
\end{equation*}
If $c^{*}$ exists, the optimal constrained welfare will thus be equal to $W(c^{*})$.

Extensions to the case of two or three selection variables are straightforward. Suppose to have two selection variables, $x$ and $z$, and two corresponding thresholds, $c_{x}$ and $c_{z}$. In this case, the assignment to treatment is a function of both thresholds, and takes on this form:

\begin{equation*}
T_{i}(c_{x},c_{z})=T^{*}_{i} \cdot \mathbf{1}[x>=c_{x}] \cdot \mathbf{1}[z>=c_{z}]
\end{equation*}
\\
The previous welfare, $W(c)$, can thus be maximized over $c$ as a vector of thresholds, one for each selection variable. In two dimensions, the previous assignment rule is called \textit{quadrant--assignment}, as it selects the upper--right quadrant of the four quadrants generated by setting the thresholds.

There are some problematic aspects that can arise when seeking the optimal thresholds. It is possible that, at the optimal thresholds, the share of the units to treat could be too small, or specularly, too large. In this case, it would be meaningless to run a policy based on such too small (or too large) number of treated. To solve the problem, the policymaker can however consider either budget constraints (for example, the maximum number of units she would be able to treat given the availability of a certain amount of money), or set in advance a targeted number of units to treat (for example, a treatment share between $30\%$ and $40\%$ of the entire reference population). This procedure might lead to a sizable reduction of the welfare but would preserve the full feasibility of the policy.

A related problem arises when, for one or more selection variables, there exists the so-called \textit{angle solution} due to a monotonic effect of the selection variable(s) on the welfare. To give a practical example of this occurrence, consider as selection variable the educational attainment of individuals. In many policy contexts, the welfare (or effect) of the policy increases monotonically with the level of education. This event would lead to select only people having the highest level of educational attainment in the sample. This can be however unfeasible for two reasons: (i) the policy aims at targeting poorly educated people; (ii) the number of treated units would become too small as people with high educational attainment are generally few. The solution to this problem is similar to the one set out above, that is, the policymaker can again introduce a threshold limit, and/or a prefixed range of treatment shares, and thus run the policy according to these constraints. Again, this procedure may entail a shrinking welfare, but would preserve the possibility to eventually run the policy.

\begin{table}[hbt]
\centering
$\begin{array}{*{20}{c}}
\hline
\; \; \; \; \; \;  {ID} \; \; \; \; \; \;  &  \; \; \; \; \; \;  {T} \; \; \; \;   \; \; & \; \; \; \; \; \;   {\tau(X)} \; \; \; \;   \; \; & \; \; \; \; \; \;   {T \cdot \tau(X)} \; \; \; \;   \; \; &  \; \; \; \; \; \;   {T^{*}} \; \; \; \;  \; \;  & \; \; \; \; \; \;    {T^{*}\cdot \tau(X)}  \; \; \; \;   \; \; \\
\hline
1 & 1 & 9  & 9 & 1  & 9\\
2 & 1 & -4 & -4  & 0  & 0 \\
3 & 1 & 5  & 5  & 1  & 5 \\
4 & 0 & 6  & 0  & 1  & 6 \\
5 & 0 & -2  & 0 & 0  & 0 \\
6 & 0 & 6 & 0  & 1  & 6 \\
\hline
 &  &  & 10  &   & 26 \\
\hline
\end{array}$
\caption{Example of an optimal policy assignment rule. The regret of this policy is equal to $16=26-10$.}
\label{tab:tab1}
\end{table}

\begin{figure}[ht]
\centering
\includegraphics[width=12cm]{TT.eps}
\caption{Distribution of $\tau$(X) and $\tau_{1}$(X). Program: National Supported Work Demonstration (NSWD). Data: LaLonde (1986). Target variable: Real earnings in 1978. Estimation technique: Regression--adjustment (with observable heterogeneity).}
\label{fig:fig1}
\end{figure}

\section{The empirical welfare maximizing protocol}
\label{sec_3}

To provide an empirical strategy for carrying out a threshold-based optimal policy assignment decision,
we suggest the following procedure consisting of five steps.

\vspace{0.5cm}

\parbox{\textwidth}{

\noindent\rule{15.5cm}{0.4pt} \newline
\textbf{Procedure}. Threshold--based optimal policy assignment  \newline
\vspace{-1cm}

\noindent\rule{15.5cm}{0.4pt} \newline

\vspace{-1cm}
\begin{enumerate}

\item Suppose to have data from an RCT or from an observational study consisting of the information triple ($Y, X, T$) available for every unit involved in the program.

\item Run a quasi--experimental method with observable heterogeneity, estimate $\tau(X)$, and compute the actual total welfare of the policy $W$.

\item Identify the optimal unconstrained policy $T^{*}$, and compute $W^{*}$, i.e. the maximum total welfare achievable by the policy, and estimate the regret as $W^{*} - W $.

\item Consider a constrained selection rule $T(x,c)$ based on a given set of selection variables, $x$, and related thresholds, $c$, and define the maximum constrained welfare as $W(x,c)$.

\item Build a greed of $K$ possible values for $c \in \{c_{1},...,c_{K} \}$, compute the optimal vector of thresholds $c_{k^{*}}$ and the corresponding maximum welfare $W(x,c_{k^{*}})$ thus achieved.

\end{enumerate}
\vspace{-0.5cm}
\noindent\rule{15.5cm}{0.4pt} \newline
}

In the application proposed in the next section, we will show how to implement this procedure on a real policy dataset.

\section{Application}
\label{sec_4}

We consider the popular experimental dataset from the National Supported Work Demonstration (NSWD) used by LaLonde (1986).  As well-known, this study looked at the effectiveness of a job training program (the treatment) administrated in 1976 on the real earnings of individuals two years after the completion of the program. It includes a set of demographic, social and economic variables at individual level such as age, race, educational attainment, previous employment condition, and real earnings - as well as the treatment indicator, and the real earnings in 1978, that is our target variable.

We first estimate the average treatment effect of this program using a regression--adjustment (RA) approach and a specification of the model including as control variables: real earnings in 1974 and 1975, age, age squared, an indicator for not having a degree, an indicator for being married, one for being black, and one for being hispanic.

As the assignment to treatment of this program was randomized, the average treatment effect (ATE) is consistently estimated by the treated--control's difference--in--mean (DIM), found to be equal to 1.79 thousand dollars. For the sake of comparison, however, we also estimate ATE by a regression--adjustment obtaining a similar estimation value, 1.76, still significant at 1\%. The RA approach allows us to estimate also the average treatment effects conditional on the covariates, i.e. $\tau(X)$ and $\tau_{1}(X)$, as well as their distributions as set out in Figure \ref{fig:fig1}.

\begin{figure}[t]
\centering
\includegraphics[width=15cm]{fig2.eps}
\caption{Computation of the policy optimal selection threshold in univariate cases. Program: National Supported Work Demonstration (NSWD). Data: LaLonde (1986). Target variable: real earnings in 1978. Univariate selection variables: real earnings in 1974, age, and educational attainment.}
\label{fig:fig2}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=15cm]{fig3.eps}
\caption{Computation of the policy optimal decision boundary in the bivariate case. Program: National Supported Work Demonstration (NSWD). Data: LaLonde (1986). Target variable: real earnings in 1978. Bivariate selection variables: real earnings in 1975 and age.}
\label{fig:fig3}
\end{figure}

\begin{figure}[ht]
\centering
\includegraphics[width=15cm]{Fig4.eps}
\caption{Computation of policy optimal decision boundaries in the bivariate case, when one of the two selection variables (age) is fixed at its optimal threshold, and the threshold of the other variable (education) is varying. Program: National Supported Work Demonstration (NSWD). Data: LaLonde (1986). Target variable: real earnings in 1978. Bivariate selection variables: age and educational attainment.}
\label{fig:fig4}
\end{figure}

Now we move on by implementing the previous empirical welfare maximizing protocol, first in a univariate, and then in a bivariate selection setting.

The univariate setting considers only one selection variable. For the sake of comparison, we run the empirical welfare maximizing protocol separately on four selection variables, i.e. real earning in 1974, age, and education. For every selection variable, we compute the optimal threshold, namely, the selection criterion that would allow the policymaker to maximize the average welfare. Figure \ref{fig:fig2} illustrates the results.

The graphs herein reported show that, for real earning in 1974, the optimal threshold is $0.96$, with corresponding maximum average welfare equal to $2.65$. Results for age show that the welfare maximizing threshold is found at age 42, with a corresponding average welfare of $2.85$. This means that, if the policymaker had selected as beneficiaries people aged 42 or larger, she would have been able to obtain an increase of the welfare by 1.09 thousands dollars per beneficiary per year, which is the difference between the pure RCT effect, 1.79, and the one obtained by selecting over age.

Results for education are particularly interesting as in this case we have an angle solution due to monotonicity. It means that the empirical welfare maximization of this program would lead to select people with the highest educational attainment possible (which is around 15 years). This is clearly unfeasible, as no one would be selected in this case. As said above, this is however a minor issue, as we could adjust our empirical welfare maximization by requiring a certain ``minimum'' share of units to treat, so to make our analysis eventually feasible. Of course, adding this constraint would lead to a reduction of the obtained level of the welfare (that, in this specific case, is rather high and equal to 4.24). As the size of the welfare is huge, we expect that constraining on a minimum share of treated units can allow to obtain a larger welfare than the one obtained by randomizing the assignment to treatment.

To explore this latter aspect, i.e. to what extent introducing a constrain on the share of treated would affect the average welfare, we consider a bivariate setting, i.e. a setting characterized by two selection variables. Before moving to this, however, it seems useful to first provide an example of a bivariate optimal treatment assignment.

Figure \ref{fig:fig3} sets out an example where we carried out an empirical welfare maximization of the NSWD jointly over age and real earnings in 1975 (with real earnings in 1978 chosen again as target variable). The figure also shows the optimal estimated decision boundary\footnote{This is an estimate of the so--called \textit{Bayesian decision boundary} as defined by supervised learning classification models (Gareth et al., 2013). In our case, this boundary should be closer to the boundary of the upper-right quadrant: it is however smoother and imprecise due to the high sparseness (few observations) located in this quadrant, with the largest part of individuals placed in the lower--left quadrant.} (the curve drawn in black). Observe that in the univariate case treated above, the decision boundary collapses to a single point. The optimal thresholds are 30.5 for age, and 10.9 for real earnings in 1975.

The upper-right quadrant represents the optimal treatment zone, and it is immediate to see that only four individuals are selected as beneficiaries. Unfortunately, they represent only the $1\%$ of the entire sample. A solution to this problem would be to provide the policymaker with a \textit{menu} of optimal treatment scenarios among which she can make pondered decisions. We show this meanu using as selction variables age and education within the same dataset.

Figure \ref{fig:fig4} sets out the optimal policy assignment decision boundary over age at different education thresholds\footnote{The choice of varying the threshold of education by letting age fixed to its optimal threshold is dictated by education monotonicity as discussed above. It goes without saying that one might do the opposite as well.}. As said, this figure can be interpreted as a menu for the policymaker to choose among alternative scenarios characterized by different policy settings entailing a trade--off between the size of the policy effect and the number of units to treat. As long as the educational attainment threshold increases, we observe in fact that the average welfare increases too. The maximum of the average welfare is obtained when only one individual is treated, with an average welfare of 3.92 and an education of 14 years. Any other intermediate scenario of this menu would entail a smaller average welfare with, however, a higher number of treated units. Given a budget constrain, for example, the policymaker can cherry--pick one of these scenarios and run the policy using this menu as reference for an ex-ante optimal re-programming of the policy treatment assignment.

\section{Conclusion}
\label{sec_5}

The literature on empirical welfare maximization is growing up rapidly. Large availability of datasets on programs already carried out, either based on observational data or randomized control trials, allows researchers and policymakers to design policies for increasing social welfare by optimally fine-tuning treatment assignment. Following in the footsteps of this recent literature, this paper has stressed the policymaker's perspective by proposing a practical implementation of an optimal policy assignment problem within the class of threshold-based selection rules. Straightforward to apply in practice and to implement with standard statistical software, the proposed procedure and illustrative application can guide policymakers to improve the ex--ante design of future policies; in other words, to learn from experience.

\newpage


\section{Introduction}

Wormholes are solutions of general relativity which possess a throat.
We define a wormhole throat as a closed 2-dimensional hypersurface of minimal area.
Such solutions were first mentioned in the works~\cite{Flamm, EnstRos, Wheeler1, Wheeler2}, but the greatest interest in these objects arose after the work of Michael Morris and Kip Thorne in 1988~\cite{MorTho}.
As noted by the authors, one of the ``exotic'' properties of these objects is that the presence of a throat leads to a violation of the energy conditions for the energy-momentum tensor of matter.
An overview of research on wormholes can be found, for example, in~\cite{Visser, LoboRew}.
Wormholes have been considered in various theories of gravity, including the Brans-Dicke theory of gravity~\cite{Agnese}, in the Einstein-Born-Infeld theory of gravity~\cite{Arellano}, in the Einstein-Gauss-Bonnet theory~\cite{Bhawal}, in $f(R)$ gravity~\cite{Miguel}, in Rastall's theory of gravity~\cite{Halder} and many other theories.
Various models of matter were considered as a source of wormholes:
scalar and electromagnetic fields~\cite{Arellano, Bronnikov1, Bronnikov2}, Chaplygin gas~\cite{ChaplyginLobo}, various models of phantom energy and quintessence, including models of matter with the structure of the energy--momentum tensor of a perfect fluid~{\cite{Sushkov1, Kuhfittig2, Lobo3, Kuhfittig:2016, Sahoo}}.
Models of wormholes were also built using the thin-shell formalism, the~works of~\cite{Visser2, Visser3} being among the first.

Most of the literature is devoted to the study of static and spherically symmetric wormholes.
An~important generalization of these studies is the study of dynamic solutions.
Dynamic wormholes have been considered in various aspects. One of the ways for building dynamic models ~\cite{Kar, KarSahdev, Kim, Roman, Kuhfittig3} is to add a time-dependent scaling factor to the metric.
Dynamic models of wormholes were constructed using the thin shell formalism~\cite{Anzhong}.
General properties of an arbitrary dynamic wormhole are considered in the papers~\cite{Hayward2, Dynwh1}.
In~recent papers~{\cite{Simpson, SimpsonLobo,SimpsonVisser}}, dynamic solutions have been constructed on the basis of a family of solutions describing regular black holes.

It is worth noting that the throat of a wormhole is defined differently by different authors.
In~the static case, the~different definitions agree with each other; in the dynamic case contradictions may arise.
The~general defenition of the throat, including for time-dependent metrics, is considered, for example, in the following
works~\cite{Dynwh1, TomikawaIzumi, Bittencourt}. We will use the definition that was used in the works~\cite{Kar, Kim, Roman}.

The~aim of this work is to construct a solution describing a dynamic wormhole in the theory of gravity with dust-like matter.
The~general solution of Einstein's equations in the theory of gravity with dust-like matter for a spherically symmetric metric was obtained by
Tolman in 1934 ~\cite{Tolman34, Lemaitre33, Bondi47, Land, CosimoBambi}; the~solution will be summarized in Section~\ref{sec2}.
This solution contains three arbitrary functions.
In~Section~\ref{sec3} we construct a solution that describes a collapsing wormhole, choosing these functions in a certain way.
In~Section~\ref{sec4} we investigate the properties of the obtained solution.

\section{Gravitational Collapse of a Dust-Like Spherical Shell}
\label{sec2}
The spherically symmetric solution of the Einstein equations in the theory of gravity with dust-like matter was obtained by Tolman in 1934 \cite{Tolman34, Lemaitre33, Bondi47, Land, CosimoBambi}.  Briefly, the solution will be presented in this section.
Dust-like matter allows the choice of a frame that is both synchronous and co-moving.
In this case, the spherically symmetric metric has the form\footnote{Henceforth, we denote the speed of light as
	$c =1$.}:
\begin{equation}\label{metric-4}
ds^2=d\tau^2-e^{\lambda(\tau,R)}dR^2-r^2(\tau,R)\big[d\theta^2+\sin^2\theta d\varphi^2\big],
\end{equation}
where $\lambda(\tau,r)$ and $r(\tau,R)$ are functions of $\tau$ and $R$.
The energy--momentum tensor of dust-like matter in the co-moving frame has the form
\begin{equation}
T_{i}^{j}=\varepsilon u_i u^j, \label{fluid}
\end{equation}
where $(u^{i})=(1,0,0,0)$ is a velocity four vector, and $\varepsilon$ is the energy density.
Below, we present the solution of the Einstein equations for the metric (\ref{metric-4}) and the energy--momentum tensor (\ref{fluid}).
The function $\lambda(\tau,R)$ has the form
\begin{equation}
e^{\lambda(\tau,R)}=\dfrac{r'^2}{1+f(R)},\label{E5}
\end{equation}
where the prime means differentiation with respect to the coordinate $R$, the dot means differentiation with respect to the coordinate $\tau$,  $f(R)$ is an arbitrary function, and $f(R)$ satisfies the condition $1+f>0$. Function $r(\tau,R)$  can be represented in parametric form.
In case $f>0$, the function has the form:
\begin{eqnarray}
&r=\dfrac{F(R)}{2f(R)}\left[ \cosh\eta-1 \right],
\quad
{\tau_0}(R)-\tau=\dfrac{F(R)}{2f(R)^{3/2}}\left[ \sinh\eta -\eta\right];&\label{Tolman1}
\end{eqnarray}
in the case of $f<0$, the function has the form:
\begin{eqnarray}
&r=\dfrac{F(R)}{2|f(R)|}\left[1- \cos\eta \right],\quad
{\tau_0}(R)-\tau=\dfrac{F(R)}{2|f(R)|^{3/2}}\left[\eta- \sin\eta\right],\label{Tolman2}&
\end{eqnarray}
where $F(R)$ and ${\tau_0}(R)$ are an arbitrary functions.
In the case of $f=0$, the function has the form:
\begin{eqnarray}
r=\left(\dfrac{9F}{4}\right)^{1/3}\left[{\tau_0}(R)-\tau\right]^{2/3}.
\end{eqnarray}

The energy density $\varepsilon$ is:
\begin{equation}
8\pi k \varepsilon=\dfrac{F'}{r'r^2 }.\label{varepsilon}
\end{equation}

The scalar curvature $\mathcal{R}$ and the Kretschman scalar $\mathcal{K}$ for the Tolman solution (\ref{metric-4}--\ref{Tolman2}) read
\begin{eqnarray}
& \mathcal{R} = \dfrac{F'}{r^2r'}, \label{R} &\\
& \mathcal{K} = \dfrac{3F'^2}{r^4r'^2}+ \dfrac{16 F^2}{r^6} - \dfrac{8FF'}{r^5r'}. \label{K}&
\end{eqnarray}

The scalar curvature and the Kretschmann scalar diverge at $r=0$ or $r'=0$ (in the case $F'\neq 0$), the solution in this case has a singularity.

The general solution (\ref{metric-4})--(\ref{Tolman2})
depend on three arbitrary functions $f(R)$, $F(R)$ and ${\tau_0}(R)$.
The solution allows an arbitrary transformation of the radial coordinate $\tilde{R}=\tilde{R}(R)$, where $\tilde{R}(R)$ is a differentiable function.
Using this, we can give any of the three functions $f(R)$, $F(R)$ or ${\tau_0}(R)$ some specific form.
Thus, only two of the three functions can be considered arbitrary.
The particular cases of the solution are the Schwarzschild solution and the Friedmann model with dust.

\section{Throat Conditions}
\label{sec3}
In this section, we will obtain the conditions under which the general solution (\ref{metric-4})--(\ref{Tolman2}) will describe the wormhole geometry.
A characteristic feature of wormholes is the presence of a throat---i.e., the presence of a space-like closed two-dimensional surface of the minimum area.
Various definitions of the mouth of a wormhole are considered in the literature \cite{Dynwh1, TomikawaIzumi, Bittencourt, Kar, Kim, Roman},
we will follow the approach used, for example, in works \cite{Kar, Kim, Roman}.

In order to determine the throat conditions, following the approach used in \cite{Roman}, we constructed an embedding diagram for the metric (\ref{metric-4}).
For convenience, we rewrote the function $f(R)$ in the form $f(R)=-b(R)/R$,  where $ b (R) $ is an arbitrary function.
Taking into account the relation (\ref{E5}), the metric (\ref{metric-4}) can be represented as
\begin{equation}\label{metric-41}
ds^2=d\tau^2-\dfrac{r'^2 dR^2}{1-b/R}-r^2\big[d\theta^2+\sin^2\theta d\varphi^2\big].
\end{equation}

Metric (\ref{metric-41}) is defined for $1-b(R)/R>0$, and for $b(R)=R$ has a singularity.
Consider the metric (\ref{metric-41}) on the spatial slice $\tau=const$ and $\theta=\pi/2$
\begin{equation}\label{3_d}
dl^2=\dfrac{r'^2 dR^2}{1-b/R}+r^2d\varphi^2,
\end{equation}
as a metric on a surface of revolution $\rho=\rho(z)$ embedded in a three-dimensional space with an Euclidean metric
\begin{equation}\label{R3d}
dl^2=dz^2+d\rho^2+\rho^2d\varphi^2,
\end{equation}
where $z$, $\rho$ and $\varphi$ are cylindrical coordinates (see Figure \ref{fig7}).
Comparing (\ref{3_d}) and (\ref{R3d}), we get
\begin{eqnarray}
&\rho^2=r^2,&\\
&dz^2+d\rho^2=\dfrac{r'^2dR^2}{1-b/R}.&
\end{eqnarray}
Taking into account that for constant $\tau=const$
\begin{eqnarray}
&d\rho|_{\tau_0}=r'_R dR,&
\end{eqnarray}
we have
\begin{eqnarray}
&\dfrac{d\rho}{dz}=\left(\dfrac{R}{b}-1\right)^{1/2},\quad \dfrac{d^2\rho}{dz^2}=\dfrac{b-R b'}{2 r' b^2}\label{ee2}.&
\end{eqnarray}
\begin{figure}[h]
\begin{center}\begin{tabular}{ccc}
\pgfplotsset{width=7cm,height=6.8cm}
\begin{tikzpicture}
\begin{axis}[

	axis y line = left,
    axis x line = middle,
		 xlabel     = $\rho$,

	 xtick=\empty, ytick=\empty,
		samples     = 160,
    domain      = 0.5:3,
    xmin = 0.5, xmax = 3,
    ymin = -1.8, ymax = 1.8,
]
\addplot[blue,thick] { sqrt(x-1)};
\addplot[blue,thick] { -sqrt(x-1)};
 \coordinate [label=right:$z$] (A) at (0cm, 5cm);
\end{axis}
\end{tikzpicture}
& \,\,\,\,\,\,\,\,\,\,&
\begin{tikzpicture}
\begin{axis}[view={60}{30},axis lines=none,colormap/cool, fill opacity=.9
	\addplot3[surf,
		samples=30,
        domain=-2:2,y domain=0:2*pi,
        z buffer=sort]
        ({(x^2+1) * cos(deg(y))},
		 {(x^2+1) * sin(deg(y))},
		 x);
\end{axis}
\end{tikzpicture}
\end{tabular}\end{center}
\caption{The embedding diagram. The function $\rho(z)$ is shown on the left panel, and the surface obtained by rotating the curve $\rho(z)$ about the axis $Oz$ is shown on the right panel.\label{fig7} }
\end{figure}

The throat of the wormhole will have the shape of a sphere, which is located at a certain value of the radial coordinates $R=R_{th}$.
On the embedding diagram, the sphere $R=R_{th}$ corresponds to a circle of radius $\rho$ on the surface of revolution;
at the throat the radius of the circle $\rho(z)$ has a minimum.
Conditions for the minimum of the function $\rho(z)$ at $R=R_{th}$ have the form:
\begin{eqnarray}
&\dfrac{d\rho}{dz}\Big|_{R_{th}}=0,\quad \dfrac{d^2\rho}{dz^2}\Big|_{R_{th}}>0.&\label{ee4}
\end{eqnarray}
Comparing (\ref{ee2}) and (\ref{ee4}), we get {throat conditions for the metric} (\ref{metric-41}):
\begin{eqnarray}
&b(R_{th})=R_{th},\label{ee3}\\
&\dfrac{1-b'(R_{th})}{r'b(R_{th})}>0\label{ee5}.&
\end{eqnarray}

We will assume that the radial coordinate takes the following positive range of values:
$0<R_{th}<R<+\infty$, where the value $R=R_{th}$ corresponds to the throat.
The condition (\ref{ee3}) implies that the function $b(R)$ is positive at the throat: $b(R_{th})>0$.
Therefore,  $b(R)$ is positive in some neighborhood of the throat by continuity.
Therefore, the function $f=-b/R$ is negative in the vicinity of the throat,
and to construct a wormhole model  we will consider the solution (\ref{Tolman2}) with a negative value of $f$.

\section{Collapsing Wormhole Model with Dust-like Matter}
\label{sec4}
The general solution (\ref{metric-4})--(\ref{Tolman2}) of the Einstein equations for a spherically symmetric metric in the theory of gravity with dust-like matter depend on three arbitrary functions $f(R)$, $F(R)$ and ${\tau_0}(R)$. In this section, we will consider a special case describing a wormhole.
We choose arbitrary functions as follows\footnote{This case is not the only one possible. Further we will restrict ourselves to considering this example as one of the simplest.}:
\begin{eqnarray}
&F=R,\quad {\tau_0}(R)={\tau_0}, \quad f(R)=-\dfrac{b}{R},\label{FT}&
\end{eqnarray}
where ${\tau_0}$, $b$ are constants, $b>0$ and $f<0$. The solution takes the form:
\begin{eqnarray}
&ds^2=d\tau^2-\dfrac{r'^2dR^2}{1-b/R}-r^2(t,R)\big[d\theta^2+\sin^2\theta d\varphi^2\big],\label{Twh1}&\\
&r=\dfrac{R^2}{2b}\left[1- \cos\eta \right],\quad {\tau_0}-\tau=\dfrac{R^{5/2}}{2b^{3/2}}\left[\eta- \sin\eta\right],\label{Twh3}&
\end{eqnarray}
where $R\in[b,+\infty)$. In this case, the first throat condition (\ref {ee3}) is satisfied for $R=b$.
The second condition (\ref{ee5}) will be satisfied when $\left.r'\right|_{R = b}>0$.
We calculate the derivative of the function $r(\tau,R)$ (\ref{Twh3}):
\begin{eqnarray}
&r'=\dfrac{R}{b(1-\cos\eta)}\left[
(1-\cos\eta)^2-\dfrac{5}{4}\sin\eta(\eta-\sin\eta) \right].\label{r'}
\end{eqnarray}

For  $\eta\in(0,2\pi)$ and $R\in[b,+\infty)$ the expression  (\ref{r'})  is positive: $r'>0$.
Therefore, for the solution (\ref{Twh1}, \ref{Twh3}), the throat conditions (\ref{ee3}, \ref{ee5}) are satisfied for $R=b$.

To build a model of a wormhole, two identical sets $M_{+}$ and $M_{-}$ are required, each of which has a metric (\ref{Twh1}, \ref{Twh3}) with the same parameter $b$
($b_{+}=b_{-}$):
$$
M_{\pm}=\left\{ (\tau,R,\theta,\varphi)| \; R\in [b,+\infty) \right\}.
$$
$M_{+}$ and  $M_{-}$ are manifolds with boundaries. The boundary is a timelike hypersurface
$$
\Sigma_{\pm}=\left\{ (\tau,R,\theta,\varphi)| \; R=b \right\}.
$$

We identify $M_{+}$ and $M_{-}$ by boundaries. We get a new manifold $M=M_{+}\cup M_{-}$.
The resulting space $M$ has a wormhole geometry that connects two space--times, in each of which the radial coordinate takes the values $R\in[b,+\infty)$.
The throat of the wormhole corresponds to the value of the radial coordinate $R=b$.

As $R\to b$, one of the metric coefficients in (\ref{Twh1}), approaches infinity.
However, the values of $r'$ (\ref{r'}) and $r$ (\ref{Twh3}) are nonzero for $\eta\in(0,2\pi)$, and the scalar curvature (\ref{R}) and the Kretschmann scalar (\ref{K}) are regular in this case. Thus, the singularity of the metric at the throat has a purely coordinate character.

It is known \cite{Visser2, Visser3, Lanczos, Darmois, Lichnerovich, Israel, Mars},
that the energy--momentum tensor on the shell $S_{ij}$ is proportinal to the Dirac delta function:
$$
T_{ij}=T_{ij}^{+}+T_{ij}^{-}+\delta(R-b)S_{ij},
$$
where $T_{ij}^{\pm}$ is the energy--momentum tensor in the corresponding region $M_{\pm}$.
$S_{ij}$ is calculated as follows:
$$
S_{ij}=\dfrac{1}{8\pi G}\left( [k_{ij}]-[k]h_{ij} \right)
$$
where $k_{ij}$ is the second fundamental form of $\Sigma$,
$[k_{ij}]$ is the difference of limiting values of the second fundamental form on both sides of the surface $\Sigma$ in $M^{\pm}$
$$
[k_{ij}]=k_{ij}|_+-k_{ij}|_-=\left[h_i^{\alpha}h_j^{\beta}\nabla_{\beta}n_{\alpha}\right],
$$
$h_{\alpha\beta}$ is a projection tensor on $\Sigma$,
$$
h_{\alpha\beta}=g_{\alpha\beta}-n_{\alpha}n_{\beta},
$$
and $n_{\alpha}$ is the normal vector to $\Sigma$.
In our case, the unit normal vector to the surface has the form
$
(n^i)=\left( 0, \pm e^{-\lambda/2},0,0 \right).
$
The tensor $[k_{ij}]$ on the surface $R=b$
$$
[k_{ij}]=
\left(
\begin{array}{cccc}
0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 \\
0 & 0 & 2r\sqrt{1-b/R} & 0 \\
0 & 0 & 0 & 2r\sqrt{1-b/R}\sin^2\theta
\end{array}\right)
$$
equals zero. Consequently, the shell momentum energy tensor
$
S_{ij}=0
$
is zero. Thus, there is no material shell on the surface $\Sigma$, and the energy--momentum tensor of the model contains only dust-like matter.

\section{Analysis of the Model}

In this section, we will consider the dynamics of the model from the point of view of an external observer.
We will assume that the dusty matter is distributed inside some spherical layer.
Outside this layer, the geometry is described by the spherically symmetric Schwarzschild solution in the vacuum \cite{Land}:
\begin{eqnarray}
&ds^2=\left(1-\dfrac{r_g}{r}\right)dt^2-\left(1-\dfrac{r_g}{r}\right)^{-1}dr^2-r^2d\Omega^2,\label{Sh}&
\end{eqnarray}
where  $r_g=2M$ is a gravitational radius, $M$ is a mass that creates a gravitational field.
In the absence of pressure, each element of matter moves along a geodesic path. All geodesics are radial due to spherical symmetry.
From the point of view of an external observer, the radial coordinate of the surface of the star is a function of time $r=r(\tau)$, where $r(\tau)$ is the radial geodesic in the Schwarzschild space--time.
It is convenient to represent the geodesic equation in parametric form \cite{Lightman}:
\begin{eqnarray}
&r=\dfrac{r_0}{2}(1+\cos\chi),\label{geod_Sh1}&\\
&\tau=\left(\dfrac{r_0^3}{8M}\right)^{1/2}(\chi+\sin\chi),\label{geod_Sh2}&\\
&u^t=\dfrac{dt}{d\tau}=\dfrac{(1-2M/r_0)^{1/2}}{1-2M/R},\label{geod_Sh3}&
\end{eqnarray}
where the parameter $\chi\in[0,\pi]$, $\tau$ is a time coordinate in the co-moving coordinate system, $r_0$ is the value of the radial coordinate at the initial moment of time $\tau=0$ ($\chi=0$).
Expressions (\ref{geod_Sh1}--\ref{geod_Sh3}) are valid provided that at the initial moment of time $\tau=0$ the shell was at rest $dr/d\tau=0$.
Substituting expressions (\ref{geod_Sh1}--\ref{geod_Sh3}) into the metric (\ref{Sh}), we find the 3-geometry of the surface of the star in the external metric:
\begin{eqnarray}
&{}^{(3)}ds^2=d\tau^2-R^2(\tau)d\Omega^2&\label{ds_out}\\
&=\dfrac{r_0^3}{8M}(1+\cos\chi)^2d\chi^2-\dfrac{r_0^2}{4}(1+\cos\chi)^2d\Omega^2.\nonumber&
\end{eqnarray}

Now we will consider the internal geometry of the dusty layer (\ref{Twh1}),   (\ref{Twh3}).
Let us assume that dust-like  matter is enclosed in a spherical layer $R\in[b,R_{0}]$, where $R_0$ is its outer radius.
The internal 3-geometry of the surface of the star is found by substituting expression (\ref{Twh3}) into the metric (\ref{Twh1}) at $R=R_0$:
\begin{eqnarray}
&{}^{(3)}ds^2=d\tau^2-r^2(R_0,t)d\Omega^2&\label{ds_in}\\
&=\dfrac{R_0^5}{4b^3}(1-\cos\eta)^2d\eta^2-\dfrac{R_0^4}{4b^2}(1-\cos\eta)^2d\Omega^2.\nonumber&
\end{eqnarray}

On the surface of the star the Schwarzschild metric should be smoothly matched with the intrinsic metric.
Comparing (\ref{ds_out}) and (\ref{ds_in}), we see that these 3-geometries are smoothly sewn together if we identify the following parameters
\begin{eqnarray}
&\eta=\pi-\chi,\quad r_0=\dfrac{R_0^2}{b},\quad M=\dfrac{R_0}{2}. \label{junctioncons}&
\end{eqnarray}

At the initial moment of time $\tau=0$ ($\chi=0$), the dust-like layer with the outer radius $r_0$ begins to collapse from the rest position towards the center.
As follows from formulas (\ref{geod_Sh1}, \ref{geod_Sh2}), at the moment of collapse, the outer radius vanishes $r=0$,
and the watch of the co-moving observer shows  $\tau_0=\dfrac{\pi r_0^{5/2}}{2b^{3/2}}$ ($\chi=\pi$).
The time $\tau_0 $ is the lifetime of the wormhole.
At the moment of collapse, solutions (\ref{Twh1}, \ref{Twh3}) are singular, since the scalar curvature and the Kretschmann scalar (\ref{R}, \ref{K}) approach infinity.
By construction, $R_0>b$, which means that due to the relations (\ref{junctioncons}) the initial outer radius of the shell will be greater than the gravitational radius---i.e., $r_0>r_g$. From the point of view of an infinitely distant observer, gravitational collapse will last an infinitely long time.

\begin{figure}[h]
\centerline{
\includegraphics[scale=0.6]{Shells.pdf}
}
\caption{The figure shows the graphs of the function $r(R,\tau)/b$ depending on $\tau/b$ for different values of $R$: $R=b$, $R=1.1b$ and $R=R_0=1.2b$.
Red-blue lines are the null radial geodesics. (See the text for more details.)}
\label{fig2}
\end{figure}

{
It is very important to answer the question: Is a photon able to cross the collapsing wormhole? To solve this problem we need to find null geodesics inside the collapsing dust ball. An equation for radial geodesics can be directly found from the condition $ds^2=0$. The metric (\ref{Twh1}) yields
\begin{equation}\label{geod_R}
\frac{dR}{d\tau}=\mp \frac{1}{r'}\,\sqrt{1-\frac{b}{R}},
\end{equation}
where $r'$ is given by Eq. (\ref{r'}). Using the relations (\ref{Twh3}), one can rewrite Eq. (\ref{geod_R}) as follows
\begin{equation}\label{geod_r}
\frac{dr}{d\tau}=\left(\frac{b}{R}\right)^{1/2}\,
\left[-\frac{\sin\eta}{1-\cos\eta} \mp\sqrt{\frac{R}{b}-1}\right].
\end{equation}
Integrating this equation gives null geodesics in the coordinates $(r,\tau)$. Note that the sign $"-"$ corresponds to ingoing geodesics moving towards the wormhole throat, while $"+"$ corresponds to outgoing geodesics moving away from the throat. At the throat $R=b$ Eq. (\ref{geod_r}) reduces to
\begin{equation}
\left.\frac{dr}{d\tau}\right|_{R=b}=-\frac{\sin\eta}{1-\cos\eta},
\end{equation}
therefore null geodesics are tangent to the cycloid $R=b$.

In the figure \ref{fig2} we illustrate the main features of the obtained solution. In particular, the dependence of the radius $r(R,t)$ of various layers of the dust-like shell $R\in[b,R_0]$ on time $\tau$ for $R_0=1.2b$ (see Eq. (\ref{Twh3})) is shown.
Region \textit{I}: $ R <b $ is not included in space--time and is not considered.
Region \textit{II}: $R\in[b,R_0]$ is filled with dust-like matter.
Region \textit{III} is the empty area outside the dust-like shell.
The three solid curves show the dynamics of the dust layers over time: curve at $R=R_0$ is the outer shell;
curve at $R=b$ is the wormhole throat; intermediate curve for the inner layer $b<R<R_0$, $R=1.1b$.
The dash-dotted line corresponds to the gravitational radius $r=r_g$.
In a finite proper time $\tau_0$, the shell reaches the gravitational radius and collapses towards the singularity $r=0$.
Red-blue lines are the null radial geodesics. The red curves are ingoing geodesic lines of photons who started their travel towards the wormhole throat from the outer shell. These photons are reaching the throat in a finite time; at the point where the photon is crossing the throat, its geodesic line becomes tangent to the cycloid $R=b$. After crossing the throat, the photon continues its movement on the other side of the wormhole space--time. In the figure \ref{fig2} the blue curves correspond to outgoing geodesic lines of photons who started their travel towards the outer shell from the throat.
}

\begin{figure}[h]
\centerline{
\includegraphics[scale=0.6]{Epsilon.pdf}}
\caption{The figure shows the graphs of the energy density $8\pi k\varepsilon b^2$ depending on $\tau/b$ for different values of $R$: $R=b$, $R=1.1b$ and $R=R_0$, $R_0=1.2b$.  As $\tau\to\tau_0$ the energy tends to infinity, the vertical asymptote is shown by dash-dot on the right in the figure.}\label{fig3}
\end{figure}
Figure \ref{fig3} shows the dependence (\ref{varepsilon}) of the energy density $\varepsilon$ on the time $\tau$ for different values of $R\in[b,R_0]$.
The energy density is positive, and at the moment of collapse $\tau\to\tau_0$, the energy tends to infinity.
The energy--momentum tensor of ``ordinary'' matter satisfies certain physical requirements, one of which is the null energy condition:
$$
T_{\mu\nu}V^{\mu}V^{\nu}\leqslant 0
$$
for any null vector $V^{\mu}$. It is known \cite{Visser}, that the null energy condition can be violated in wormhole space--time.
Consider a region in the vicinity of the wormhole throat $R\in[b,R_0]$, where the model is described by solutions (\ref{Twh1}, \ref{Twh3}).
We choose a null vector $V^{\ mu}=\left( 1,e^{-\lambda/2},0,0 \right)$ and calculate the following contraction:
\begin{equation}
T_{\mu\nu}V^{\mu}V^{\nu}=\varepsilon,\label{NEC}
\end{equation}
which is positive for $R\in[b,R_0]$, $\tau\in [0,\tau_0]$. Thus, in this model, the null energy condition is fulfilled.

We also present this metric in the limiting case $\tau\to\tau_0$ ($\eta\to+0$, $\chi\to\pi-0$):
\begin{eqnarray}
&
r(\tau,R)=R^{1/3}\left[\dfrac{3}{2}\left(\tau_0-\tau\right)\right]^{2/3},
&\\
&
ds^2=d\tau^2-\left(\dfrac{3}{2}\dfrac{\tau_0-\tau}{R}\right)^{4/3}\left[\dfrac{1}{12^{2/3}}\dfrac{dR^2}{1-b/R}+R^2\left(d\theta^2+\sin^2\theta d\varphi^2\right)\right]
&
\end{eqnarray}

\section{Conclusion}
It is well known that static wormhole configurations in general relativity (GR) are possible only if matter threading the wormhole throat is ``exotic''---i.e., violates a number of energy conditions. For this reason, it is impossible to construct {static} wormholes supported only by dust-like matter which satisfies all usual energy conditions. However, this is not the case for non-static configurations.
In 1934, Tolman found a general solution describing an evolution of a spherical dust shell in GR. In this particular case, Tolman's solution describes the collapsing dust ball; the inner space--time structure of the ball corresponds to the Friedmann universe filled by a dust.
In the present work we have used Tolman's solution in order to construct a dynamic spherically symmetric wormhole solution in GR with dust-like matter.
For this aim, we have considered the specific subclass of arbitrary functions $f(R)$, $F(R)$ and $t_0(R)$ constituting Tolman's solution,
and written solutions (\ref{Twh1})--(\ref{Twh3}), which satisfy the throat conditions.
The solution constructed represents the collapsing dust ball with the inner wormhole space--time structure.

To construct a wormhole model we considered two identical solutions $M_{+}$ and $M_{-}$ with metrics (\ref{Twh1}), (\ref{Twh3}), in each of which the radial coordinate $R$ takes the range of values $R\in[b,+\infty)$. The sets $M_{+}$ and $M_{-}$ are identified by their boundaries $\Sigma$. The result is a new manifold $M=M_{+}\cup M_{-}$ that has a wormhole geometry.
The throat of the wormhole corresponds to the junction surface $\Sigma$: $R=b$.
Due to the smoothness of the functions, the identity of $M_{+}$ and $M_{-}$ and to the throat conditions, the energy--momentum tensor $S_{ij}$ of the shell $\Sigma$ is equal to zero. Thus, the energy--momentum tensor of the constructed model contains only dust-like matter, and does not contain a contribution of the Dirac delta function on the junction surface.
For $\eta\in(0,2\pi)$, the scalar curvature (\ref{R}) and the Kretschmann scalar (\ref{K}) are regular, and the solution does not have singularities.

We assumed that the dusty matter is enclosed inside some spherical layer and considered the dynamics of a wormhole from the point of view of an external observer.
Gluing solutions (\ref{Twh1}), (\ref{Twh3}) and the Schwarzschild solution, we obtained the condition (\ref{junctioncons}) that connects the model parameters, the gravitational radius and the mass of the dust ball.
The dynamics of dust layers over time is depicted in Figure \ref{fig2}.
At the initial moment of time, a dusty layer begins to collapse towards the center, at the moment of collapse the outer radius of the shell vanishes, and the clock of the co-moving observer shows $\tau_0={\pi r_0^{5/2}}/{2b^{3/2}}$. The $\tau_0$ is the wormhole's lifetime.
From the point of view of an infinitely distant observer, gravitational collapse will last indefinitely.
At the moment of collapse, the scalar curvature and the Kretschmann scalar (\ref{R}), (\ref{K}) become singular and the solution has a singularity.
Figure \ref{fig3} shows the dependence of the energy density of the dust layer on time, the energy density is positive, and tends to infinity at the moment of collapse. It is shown that the null energy condition for this model is satisfied (\ref{NEC}). Therefore, the dust-like matter that the ball is made of satisfies the usual energy conditions and cannot prevent the collapse.

\acknowledgments{We acknowledge the contribution of Evdokim Isanaev on the very beginning of this work. P.E.K. and S.V.S. are supported by RSF grant No. 16-12-10401. Partially, this work was done in the framework of the Russian Government Program of Competitive Growth of the Kazan Federal University.}


\section{Introduction}
\label{sec:introduction}

The always-increasing capacity of wireless systems is promoting the design of applications and services with increasingly challenging demands, such as video streaming, teleconference, telepresence, \acrfull{xr}, among others~\cite{tgayUsageModel}.
In order to meet the demand in terms of data rate of such applications, the latest versions of the \gls{wifi} standard, i.e., IEEE~802.11ad and 802.11ay, also known as \gls{wigig}, offer the possibility to communicate over the mmWave band at 60~GHz, where multiple 2.16~GHz channels are available.
By taking advantage of techniques such as channel bonding and \gls{mimo}, and by introducing novel features to the protocol stack, these standards can provide data rates over 100~Gbps~\cite{11ayMagazine}.

However, many applications also have very stringent \gls{qos} requirements, in particular in terms of delay and jitter, which may be incompatible with the stochastic nature of contention-based channel access mechanisms generally supported by legacy \glspl{wlan}.
To address this problem, the \gls{wigig} standards introduced a contention-free access mechanism that allows a \gls{sta} to reserve radio resources at regular time intervals.
These resources are organized in blocks, called \glspl{sp}, and the standards only specify the basic procedures to ensure inter-vendor compatibility, leaving the design and the implementation of the scheduler to the manufacturers.

Regarding the practical design, handling multiple periodic traffic streams can be problematic, especially when traffic flows with different periodicities coexist.
In this case, it is necessary to anticipate collisions among different periodic allocations and either adjust them or, in the worst case, reject new incompatible requests.
Furthermore, even if a collection of requests with identical traffic requirements is considered, upon receiving a new request, the scheduler needs to decide whether to rearrange the previously allocated resources to improve fairness and efficiency, or to maintain the original schedule and then best accommodate the new request, in order not to perturb the pre-existing streams but potentially reaching suboptimal resource allocation.
Moreover, \glspl{sp} are subject to a number of constraints, described in \cref{sec:framework_description}, which need to be accounted for when designing and optimizing scheduling algorithms.

With these challenges in mind, in this work we address both admission control and resource allocation for multiple periodic traffic sources, following the constraints given by the \gls{wigig} standards.
Specifically, we cast the periodic scheduling problem within the \gls{wigig} allocation framework and design a simple and efficient algorithm to check for the feasibility of a new request.
We then propose a simple admission control algorithm with limited scheduling capabilities, as well as a more elaborate and optimized strategy to increase the admission rate and, possibly, the fairness among independent flows.
Finally, we compare these two policies and shed some light on basic trade-offs.

The rest of the paper is organized as follows.
The resource allocation framework is described in \cref{sec:framework_description}, while \cref{sec:state_of_the_art} provides an overview of the State of the Art on related problems, and motivates our need to fill the gap of the current literature in this topic.
Then, we present the proposed algorithms in \cref{sec:scheduling_algorithms}. Performance analysis is presented in \cref{sec:results}.
Finally, in \cref{sec:conclusions} we draw our conclusions and propose possible extensions of this work.

\section{Framework Description}
\label{sec:framework_description}
Based on~\cite{standard802.11_2016}, \glspl{sta} can request the \gls{ap} to reserve periodic transmission intervals by sending a control frame containing the required periodicity ($\p{}$) and the minimum and maximum duration of each allocation ($[\Tmin{}, \Tmax{}]$).
The \gls{ap} advertises the allocated \glspl{sp} to the \glspl{sta} at each \gls{bi}, specifying the starting time, duration, and periodicity of each block.
The allocation needs to comply with a number of constraints:
\begin{enumerate}
 \item Periodicity ($\p{}$) can only be an integer multiple ($\p{} \in \mathbb{N}$) or an integer fraction ($\p{}^{-1} \in \mathbb{N}$) of a \gls{bi} (\ensuremath{T_{\rm BI}}), thus the block periodicity interval will be $\Tp{} = \p{} \, \ensuremath{T_{\rm BI}} $;
 \item Allocation blocks cannot be scheduled across the \gls{bi} boundaries;
 \item The allocated block duration $\Tblk{}$ should fall in the range $[\Tmin{}, \Tmax{}]$ specified in the resource request.
\end{enumerate}
A more detailed description of the constraints imposed by the standard can be found in~\cite{standard802.11_2016} and~\cite{mohebi2020challenges}.

Since this work is focused on  allocation algorithms for periodic traffic sources, we neglect the \glspl{cbap}, which is present in each \gls{bi} for asynchronous traffic.
In addition, to compare the scheduling algorithms in challenging conditions, we assume that the allocated resources will be maintained indefinitely, so that the channel load increases progressively as new resource reservations are accepted.
For the sake of simplicity and clarity, we also assume that the allocation blocks of a given accepted request are not fractioned into multiple disjoint intervals (i.e., each \gls{sp} will consist of a single time interval of duration $\Tblk{}$).
Furthermore, we consider a \textit{strict periodicity} constraint, which prevents the scheduler from changing the starting time of already allocated blocks, while the block duration $\Tblk{}$ can be freely changed within the interval $\qty[\Tmin{}, \Tmax{}]$.

\begin{figure}[tp]
 \centering

 \newcommand\fscale{1}
 \input{img/notation_joint_alloc.tex}
 \caption{Example of allocations $A_1$ (orange) and $A_2$ (blue), with $\p{2}=2\p{1}$.}
 \label{fig:joint_alloc}
\end{figure}

\section{State of the Art}
\label{sec:state_of_the_art}

Many works analyzing the \gls{mac} layer of the \gls{wigig} standards focus mostly on \glspl{cbap}, which extend the traditional WiFi access to cope with directional communications, either neglecting \gls{sp} allocations or considering extremely simple allocation schemes.
In \cite{babich2009throughput}, the authors proposed a mathematical framework for the analysis of \gls{e2e} metrics in 802.11-based systems, comparing throughput and average packet delay in scenarios where the nodes are equipped with advanced antenna systems.
The characteristics of the \gls{dcf} were taken into account, for which a theoretical performance analysis was carried out in \cite{bianchi2000performance}.
Instead, the authors in~\cite{pielli2019analytical} present a model to assess the performance of \glspl{cbap} for the IEEE 802.11ad standard, taking into account a directional channel model and the presence of scheduled \glspl{sp}, but they do not focus on how to assign such \glspl{sp}.

To the best of our knowledge, little work has been done on contention-free scheduling for \gls{wigig} networks.
In~\cite{hemanth2013performance,rajan2016saturation}, the authors analyze the case where all contention-free allocations occupy the beginning of each \gls{bi}, while the rest of the interval is left for a single \gls{cbap}.
This allocation strategy, however, cannot support requests for periodic resource allocations with time period shorter than $\ensuremath{T_{\rm BI}}{}$.
The authors of~\cite{khorov2016mathematical}, instead, propose an accurate mathematical analysis of the performance of a realistic \gls{vbr} traffic source in the presence of channel errors, when using a periodic resource allocation scheme, but do not tackle the problem of scheduling multiple periodic allocations at once.

On the other hand, the problem of periodic scheduling has been widely studied in other areas, such as real-time computation and task scheduling, where the objective is to complete tasks within a given time, while minimizing the resource utilization.
For example, the authors of~\cite{liu73scheduling} develop and compare heuristic algorithms for scheduling tasks with hard periodic deadlines and constant resource utilization, showing that a deadline-first approach ensures maximum resource utilization.
In~\cite{ramamritham95allocation}, the authors try to schedule safety-critical periodic tasks with precedence constraints, distributed over multi-processor systems using an adapted deadline-first approach, while the authors of~\cite{cheng95allocation} use simulated annealing to optimize a similar problem.
Finally, \cite{zhu03multiple} finds a low-overhead optimal solution (from a resource utilization point of view) assuming that tasks have a fixed resource requirement.

All these approaches, however, cannot be directly used in \gls{wigig} systems, either because they are not compliant with the constraints imposed by the resource allocation procedures (i.e., granularity of the allocation periods, \gls{bi} boundaries), or because they cannot exploit the \gls{wigig} standards' flexibility (e.g., the dynamic allocation of $\Tblk{}$). This work contributes to fill the gap by proposing admission control and scheduling algorithms that account for the specific features of \gls{mmw} \glspl{wlan}.

\section{Scheduling Algorithms}
\label{sec:scheduling_algorithms}
\begin{figure}[t!]
 \newcommand\fscale{1}
 \centering
 \begin{subfigure}[b]{\columnwidth}
  \centering

  \input{img/scheduling_alg_fail_0.tex}
  \caption{Scheduling step 0.}
  \label{fig:feasibility_check_fail_0}
 \end{subfigure}
 \\
 \begin{subfigure}[b]{\columnwidth}
  \centering

  \input{img/scheduling_alg_fail_1.tex}
  \caption{Scheduling step 1.}
  \label{fig:feasibility_check_fail_1}
 \end{subfigure}
 \\
 \begin{subfigure}[b]{\columnwidth}
  \centering

  \input{img/scheduling_alg_fail_2.tex}
  \caption{Scheduling step 2.}
  \label{fig:feasibility_check_fail_2}
 \end{subfigure}
 \caption{Feasibility check for an \textit{infeasible} pair of allocations, where $A_1$ (blue) was a pre-existing allocation with $\p{1}=\frac{1}{2}$, and the algorithm is checking whether a new allocation $A_2$ (orange) with $\p{2}=\frac{1}{3}$ is compatible.}
 \label{fig:feasibility_check_fail}
\end{figure}

\begin{algorithm}[t!]
 \footnotesize
 \caption{Feasibility check under strong periodicity conditions (see \cref{fig:feasibility_check_fail,fig:feasibility_check_succ}).}
 \label{alg:strict_periodicity_feasibility}
 \begin{algorithmic}[1]
  \REQUIRE $\qty{A_1,\ldots, A_{N-1}}$ (fixed), $A_N$ (new allocation), $\qty[\ensuremath{t_{\rm min}}, \ensuremath{t_{\rm max}}]$

  \STATE Compute $\Tp{1,\ldots,N}$
  \STATE $\tstart{N} = t_{0, \rm start}^N \leftarrow \ensuremath{t_{\rm min}}$ \COMMENT{\cref{fig:feasibility_check_fail_0}}

  \WHILE{$t_{0, \rm end}^{N} < \ensuremath{t_{\rm max}}$}
  \STATE Check for collisions in $\left[\ensuremath{t_{\rm min}}, \ensuremath{t_{\rm min}} + \Tp{1,\ldots,N} \right)$
  \IF{no collisions}
  \STATE $\ensuremath{t_{\rm feas}} \leftarrow \tstart{N}$
  \RETURN $A_N$ is feasible with starting time $\ensuremath{t_{\rm feas}}$ \COMMENT{\cref{fig:feasibility_check_succ}}
  \ELSE
  \STATE Allocation block $h \in A_N$ collides with allocation block $k \in A_i$, for some $i \in {1,\ldots,N-1}$
  \STATE $\Delta t \leftarrow t_{k, \rm end}^i - t_{h, \rm start}^N$
  \STATE $\tstart{N} \leftarrow \tstart{N} + \Delta t$
  \ENDIF
  \ENDWHILE

  \RETURN $A_N$ is not a feasible allocation \COMMENT{\cref{fig:feasibility_check_fail}}
 \end{algorithmic}
\end{algorithm}

We denote by $A_n= (\tstart{n}, \Tp{n}, \Tblk{n})$ the allocation for the $n$-th traffic stream, where $\tstart{n}$ is the starting epoch, $\Tp{n}$ is the period, and $ \Tblk{n}$ is the allocated duration of each individual block.
Therefore, the allocation consists of a sequence of blocks, where the $k$-th block of the $n$-th traffic stream takes the interval $b_k^n = \qty(t_{k, \rm start}^n, t_{k, \rm end}^{n})$, where
\begin{equation}
 \begin{aligned}\label{eq:alloc_blk_start_end}
  t_{k, \rm start}^n & = \tstart{n} + k \Tp{n}           \;;  \\
  t_{k, \rm end}^{n} & = \tstart{n} + k \Tp{n} + \Tblk{n} \;;
 \end{aligned}
\end{equation}
for $k=0,1,2,\ldots\;$.
A graphical example is shown in \cref{fig:joint_alloc}.

Following this definition we can say that, given $N$ distinct allocations $A_1,\ldots,A_N$, they are jointly periodic over a period
\begin{equation}
  T_p^{1,\ldots,N} = \mathrm{lcm}\left(T_p^i, \ldots, T_p^N \right),
\end{equation}
where $\mathrm{lcm}$ indicates the \textit{least common multiple} of the periods.
Note that, since all block periods are integer multiples or fractions of \ensuremath{T_{\rm BI}}, the \gls{lcm} can always be properly defined~\cite{zazkis2015lcm} as
\begin{equation}
  \mathrm{lcm}\left( \frac{a}{b}, \frac{c}{d} \right) = \frac{\mathrm{lcm}(a,c)}{\mathrm{gcf}(b,d)},
\end{equation}
where $\mathrm{gfc}$ is the \textit{greatest common factor}.

Given the periodicity of the allocation patterns, a new allocation $A_N$ should start within a time interval $\Tp{N}$ since the beginning of the \gls{bi}.
Moreover, a necessary requirement for admission is that in an interval of duration $T_P^{1,\ldots,N}$, no block $b_h^N$ overlaps with any block $b_k^n$, $n \in \{1,,\ldots,N-1\}$, $\forall h,k\geq 0$.

The remainder of this section is structured as follows: in \cref{sub:feasibility_check_algorithm} we will illustrate an algorithm to efficiently check whether a new allocation is compatible with a pre-existing schedule, in \cref{sub:simple_scheduler} we will present a simple scheduling algorithm, and finally in \cref{sub:max_min_fair_scheduler} we will describe in detail a more complex algorithm that aims at minimizing the rejection of new allocations under the \textit{strict periodicity} assumption.

\subsection{Feasibility Check Algorithm}
\label{sub:feasibility_check_algorithm}

\begin{figure}[t!]
 \newcommand\fscale{1}
 \centering
 \begin{subfigure}[b]{\columnwidth}
  \centering

  \input{img/scheduling_alg_succ_0.tex}
  \caption{Scheduling step 0.}
  \label{fig:feasibility_check_succ_0}
 \end{subfigure}
 \\
 \begin{subfigure}[b]{\columnwidth}
  \centering

  \input{img/scheduling_alg_succ_1.tex}
  \caption{Scheduling step 1.}
  \label{fig:feasibility_check_succ_1}
 \end{subfigure}
 \vspace{0.1pt}
 \\
 \begin{subfigure}[b]{\columnwidth}
  \centering

  \input{img/scheduling_alg_succ_2.tex}
  \caption{Feasible intervals.}
  \label{fig:feasibility_check_succ_2}
 \end{subfigure}
 \caption{Feasibility check for a \textit{feasible} pair of allocations, where $A_1$ (blue) was a pre-existing allocation with $\p{1}=\frac{1}{4}$, and the algorithm is checking whether a new allocation $A_2$ (orange) with $\p{2}=\frac{1}{2}$ is compatible.}
 \label{fig:feasibility_check_succ}
\end{figure}

This feasibility check can be performed as described in \cref{alg:strict_periodicity_feasibility}, whose arguments consist of the existing allocations $A_1,\ldots,A_{N-1}$, the new request $A_N$ and a search interval $\qty[\ensuremath{t_{\rm min}}, \ensuremath{t_{\rm max}}]$.
For reasons that will be clear later, we assume that the existing allocations cannot be changed, while for $A_N$ only $\tstart{N}$ can be modified, keeping the period $T_p^N$ and the block duration $\Tblk{N}$ fixed.
Based on these input values, the aim of the procedure is to find the earliest feasible starting time $\ensuremath{t_{\rm feas}}^N$ such that a block of duration $\Tblk{N}$ fits in the search interval.

To do so, starting from $\ensuremath{t_{\rm min}}$, the algorithm progressively shifts the starting time by an interval $\Delta t$ (described in \cref{alg:strict_periodicity_feasibility}) until either all feasibility conditions are met, or $t_{0, \rm end}^{N} > \ensuremath{t_{\rm max}}$, in which case the allocation $A_N$ with block duration $\Tblk{N}$ is rejected.

A trivial example involves $A_1$, i.e., the first received allocation request from a \gls{sta}.
In this case, $\ensuremath{t_{\rm min}}$ will be set to the start time of the first \gls{bi} following the reception of the request, while $\ensuremath{t_{\rm max}} = \Tp{N}$ to guarantee the periodicity.
Since no previous allocated \glspl{sp} exist, $A_1$ is immediately accepted with $\ensuremath{t_{\rm feas}} = \ensuremath{t_{\rm min}}$.
It is important to highlight that, however, by choosing specific combinations of input parameters, \cref{alg:strict_periodicity_feasibility} can be used also by more advanced scheduling schemes, as explained later.

Given any feasible starting time $\ensuremath{t_{\rm feas}}$, it is useful to compute the rightmost \textit{boundary} of the allocation, i.e., the largest interval $[\ensuremath{t_{\rm feas}},\ensuremath{t_{\rm lim}}]$ that would still make $b_0^N \in [\ensuremath{t_{\rm feas}},\ensuremath{t_{\rm lim}}]$ and, in turn, $A_N$ feasible, even for larger values $\tstart{N}$ and $\Tblk{N}$.
This boundary can be computed by finding the minimum distance between each $b_h \in A_N$ and each $b_k \in A_n$, $n \neq N$.
The final results will be the minimum measured distance.
A graphical illustration of how the algorithm behaves when the new request is infeasible is shown in \cref{fig:feasibility_check_fail}, while a new feasible request is shown in \cref{fig:feasibility_check_succ}.

Following this definition and the above numerical example, the first allocation request to be generated, i.e., $A_1$, will find itself in the optimal condition where $\ensuremath{t_{\rm feas}}=\ensuremath{t_{\rm min}}$ and $\ensuremath{t_{\rm lim}}=\ensuremath{t_{\rm max}}$.

In general, multiple feasible intervals exist.
To find an exhaustive list, we can iterate \cref{alg:strict_periodicity_feasibility} with $\tstart{N}$ initialized to the start time of the \gls{bi}, and progressively updated at each iteration with the value of $\ensuremath{t_{\rm lim}}$ found in the previous execution.
This procedure continues until the shift of $\tstart{N}$ leads to an infeasible allocation.
We define the list of feasible intervals (which depend on $\Tblk{N}$) as $\mathcal{I}_N = \qty{I_1^N,\ldots,I_M^N}$, where $I_m^N = \qty[\tfeasmn{m}{N}, \tlimmn{m}{N}]$, $m=1,\ldots,M$ (see \cref{fig:feasibility_check_succ_2}).
Hence, each of these intervals delimits the finite number of intervals in which the new allocation $A_N$ can be fitted, considering all previous allocations.
A good scheduling algorithm should then assess which interval yields the best overall performance, possibly trying to optimize a target \gls{kpi}.

\subsection{Simple Scheduler}
\label{sub:simple_scheduler}
The first scheduler that we propose assumes that the block duration and periodicity of already accepted traffic streams cannot be varied.
Then, a new request $A_N$ with a block duration of $\Tblk{N} \in [\Tmin{N}, \Tmax{N}]$ can be accepted only if there exists a feasible interval in $\mathcal{I}_N $ with a duration of at least $\Tmin{N}$.
Therefore, the maximum amount of available resources that can be allocated to $A_N$ is determined by the longest feasible interval, or by $\Tmax{N}$, whichever is smaller; $\tstart{N}$ and $\Tblk{N}$ need to be set accordingly.
We can already notice that, using this simple first-come-first-served approach, the latest requests are highly disadvantaged if the first ones require big slices of time resources.
In the long term, as we will see in \cref{sec:results}, this could lead not only to poor performance in terms of fairness, but also to a very low admission rate.

\subsection{Max-Min Fair Scheduler}
\label{sub:max_min_fair_scheduler}
A more flexible approach consists in dynamically adapting the duration of the allocated intervals within the admissible range, $\Tblk{n} \in [\Tmin{n}, \Tmax{n}]$, $\forall n=1,\ldots,N$, in order to distribute time resources among all traffic streams in a fairer manner.

Consider the following parameterized block duration:
\begin{equation}\label{eq:parameterized_tblk}
 \Tblk{n}(r) = \Tmin{n} + r_n(\Tmax{n} - \Tmin{n}), \quad r_n \in [0,1] \;.
\end{equation}
We consider a scheduler to be fair if $r_{\rm min} = \min_{n} \{r_n\}$ cannot be increased without breaking the limits imposed by some allocation under the \textit{strict periodicity} constraint (see \cref{sec:framework_description}).
The scheduling algorithm, then, should assign the largest possible \gls{sp} to each allocation, while respecting all the constraints.\footnote{Note that if $\Tmin{n}=\Tmax{n}$, $r_n$ has no meaning. For simplicity, this case has not been included in this study.}

To fit a new traffic stream, the pre-existing allocations will thus have to either maintain or reduce their block duration, depending on whether and how the new allocation collides with them.
This will lead to a lower rejection rate with respect to the \textit{\nameref{sub:simple_scheduler}} (\cref{sub:simple_scheduler}), and more fairness among requests distributed in time.

The proposed algorithm is here presented in two parts: the first part describes how the allocation scheme works (\cref{ssub:allocation_algorithm}), while the second part describes the fairness paradigm (\cref{ssub:optimally_fair_allocation}).

\subsubsection{Allocation Algorithm}
\label{ssub:allocation_algorithm}

Differently from the simple scheduler, this scheduler can change the block duration within the range imposed by the requesting \gls{sta}, i.e., $\Tblk{n} \in [\Tmin{n}, \Tmax{n}]$.
To reduce the rejection rate, we check the feasibility of a new allocation $A_N$ (\cref{sub:feasibility_check_algorithm}) by assuming all existing allocations are shrunk to their minimum, i.e., $\Tblk{n} = \Tmin{n}$ for $n=1,\ldots,N$.
If $A_N$ is infeasible even under these conditions, then the allocation cannot be granted without disattending the requests of some previously accepted flow.
Therefore, $A_N$ is rejected.
Conversely, if $A_N$ is feasible, it gets accepted, and in a later step the algorithm will try to increase the resource utilization of all allocations fairly.

From now on, we use the symbol $^*$ to represent the parameter values at the end of the execution of the algorithm.
We recall that, based on the strict periodicity assumption, the starting times of the already allocated blocks cannot change

Note that, given a set of feasible allocations, reducing any $r_n$ (and, in turn, the $\Tblk{n}$) still yields a valid configuration.
Similarly, a valid configuration for $A_N$ with $\tstart{N}$ and $r_N \geq 0$ will remain valid if $\tstart{N*} \geq \tstart{N}$ and $\tend{N*} = \tstart{N*} + \Tblk{N}(r_N^*) \leq \tend{N}$.
We thus consider the following constraints:
\begin{subequations}\label{eq:optimally_fair_alloc_constraints}
 \begin{align}
  \label{eq:rn_bound}
  r_n^*       & \leq r_n, \, \forall n \leq N;        \\
  \label{eq:tstart_bound}
  \tstart{N*} & \geq \tstart{N}; \\
  \label{eq:tend_bound}
  \tend{N*}   & \leq \ensuremath{t_{\rm lim}}^N.
 \end{align}
\end{subequations}

\floatstyle{spaceruled
\restylefloat{algorithm
\begin{algorithm}[t!]
 \footnotesize
 \caption{Max-min fair scheduling.}
 \label{alg:complex_allocation}
 \begin{algorithmic}[1]
  \REQUIRE $A_1,\ldots, A_N$, $\Tp{1,\ldots,N}$

  \STATE Compute $\mathcal{I}_N$ considering $\Tblk{n}=\Tmin{n} \forall n=1,\ldots,N$
  \FORALL{$I_m^N = \qty[\ensuremath{t_{\rm feas}}, \ensuremath{t_{\rm lim}}]^N_m \in \mathcal{I}_N$}
  \STATE $\tstart{N} \leftarrow \ensuremath{t_{\rm feas}}$
  \STATE Set $r_N$ such that $\Tblk{N} = \min\qty{\Tmax{N}, \ensuremath{t_{\rm lim}} - \ensuremath{t_{\rm feas}}}$ \COMMENT{\cref{eq:parameterized_tblk}}

  \FORALL{$A_n, \; n=1,\ldots,N-1$}

  \FORALL{(block $k \in A_n$) $\in \Tp{1,\ldots,N}$}
  \IF{block $k$ collides with $A_N$}
  \STATE Update $r_n^*$, $r_N^*$, $\tstart{N*}$ \COMMENT{\cref{ssub:optimally_fair_allocation}}

  \IF{$r_n^* < r_n$}
  \STATE Add/update $A_n$ to a list $\mathcal{C}$ of colliding allocations
  \STATE Memorize $r_{n,\rm prev} \leftarrow r_n$
  \ENDIF

  \STATE Update $r_n$, $r_N$, $\tstart{N}$
  \ENDIF

  \ENDFOR
  \ENDFOR

  \FORALL{$A_n \in \mathcal{C}$}
  \STATE Compute $\Delta t= \ensuremath{t_{\rm lim}} - \tstart{n}$ for $A_n$ given $A_N$ \COMMENT{see \cref{sec:scheduling_algorithms}}
  \STATE $\Tblk{n} \leftarrow \min\qty{\Tblk{n}(r_{n,\rm prev}), \Delta t}$ \COMMENT{Try to improve the allocation duration if $A_N$ has been further reduced}
  \ENDFOR

  \STATE Compute allocation score $s_m \leftarrow \min_{n=1,\ldots,N} r_n$
  \ENDFOR

  \RETURN The configuration which maximizes the allocation score $\qty{s_m}$
 \end{algorithmic}
\end{algorithm}

The algorithm starts by considering the first feasible interval $I_1^N$, which ensures a valid configuration when $r_n=0$, $\forall n \leq N$.
The new request is temporarily accepted with $\tstart{N}=\tfeasmn{1}{N}$ and maximum possible $r_N$, such that $\Tblk{N} = \min\qty{\Tmax{N}, \ensuremath{t_{\rm lim}} - \ensuremath{t_{\rm feas}}}$.

Then, the algorithms try to re-balance the resource allocation by increasing all $\{ r_n, \, \forall n \leq N \}$ to their previous values.
Given that feasible intervals $\mathcal{I}_N$ were computed considering all allocations with minimum duration, though, setting $\tstart{N}=\tfeasmn{1}{N}$ may (or may not) create a collision with a generic $A_i$ when setting $r_i\geq 0$ back to its previous value.

On the other hand, thanks to the information given by $\ensuremath{t_{\rm lim}}$, we can always choose $r_N$ such that the new allocation does not collide with a previous one, on the right.

Collisions can be found iteratively over each block of each previous allocation in a joint period.

If for a certain block $b_k^n \in A_n$ and a block $b_h^N \in A_N$ we have
\begin{equation}
 t_{k, \rm start}^{n} + \Tblk{n}(r_n) \geq t_{h, \rm start}^{N}.
\end{equation}
then the two allocations are in conflict, as shown in \cref{fig:collision_baseline}.
In this case, $\tstart{N*}$, $r_N^*$, and $r_n^*$ have to be updated following the constraints in \eqref{eq:optimally_fair_alloc_constraints}, as described in \cref{ssub:optimally_fair_allocation}.

The constraints from \eqref{eq:optimally_fair_alloc_constraints}, the existence of a non-empty set of feasible intervals, and the iterative nature of the problem ensure that the algorithm will stop in a finite time with a valid configuration.
Since each feasible interval $I_m^N \in \mathcal{I}_N$ has one locally fairest configuration, the exhaustive search described in \cref{alg:complex_allocation} is able to find the globally fairest configuration by exhaustive search.

\begin{figure}[tbp]
  \centering

  \newcommand\fscale{1}
  \input{img/collision_baseline.tex}
  \caption{Representation of a collision between $A_n$ and $A_N$.}
  \label{fig:collision_baseline}
 \end{figure}

\subsubsection{Optimally fair allocation}
\label{ssub:optimally_fair_allocation}

In this section, we will discuss how fairness can be achieved given a pair of colliding allocations $A_n$, $n\in \{1,\ldots,N-1\}$, and $A_N$.
In \cref{ssub:allocation_algorithm} we explained how such a collision can be found, e.g, between blocks $b_k^n$ and $b_h^N$.
For the sake of clarity, in this section we will drop the notation for the specific colliding blocks.

In order to fully exploit the available resources, looking at \cref{fig:collision_baseline}, we force $\tstart{N} = \tend{n}$ and $\tend{N} \leq \ensuremath{t_{\rm lim}}$.
In this way we make $A_N$ start right after $A_n$, still respecting the limits imposed by $\ensuremath{t_{\rm lim}}$.

While possibly not being optimal, this is still a sensible choice for a greedy approach that tries to maximize the fairness of the current configuration.
By doing so and imposing $r_n = r_N = r^*$, what we call the \textit{fairness equation}, and by noting that $r^*\leq 1$ should hold, we have that
\begin{equation}\label{eq:bounded_rstar}
 r^* = \min \qty{1, \frac{\ensuremath{t_{\rm lim}} - \tstart{n} - \Tmin{n} - \Tmin{N}}{\qty(\Tmax{n} - \Tmin{n}) + \qty(\Tmax{N} - \Tmin{N})}}.
\end{equation}
We call $r^*$ the \textit{fair allocation ratio}, and note that if $r^*<1$, it must be that $\tend{N}=\ensuremath{t_{\rm lim}}$, whereas if $r^*=1$ in general $\tend{N} \leq \ensuremath{t_{\rm lim}}$ by construction.

Depending on the initial conditions of the problem, there is a number of different cases which have to be properly managed in order to obtain a fair distribution of resources.

First of all, if $r_n \leq r^*$, following \cref{eq:rn_bound}, it means that previous adjustments do not make it possible for $A_n$ to obtain more resources while still ensuring a valid configuration, and thus $r_n^* = r_n$.
Furthermore, since we assume that a collision happens between $A_n$ and $A_N$ with this configuration, $A_N$ has to be delayed setting $\tstart{N*} = \tstart{n} + \Tblk{n}(r_n^*) > \tstart{N}$.

In case also $r_N \leq r^*$, allocation $A_N$ cannot be extended either.
Since both allocations have $r_n, r_N \leq r^*$, they will both surely fit in the feasible interval.
If, instead, $r_N > r^*$, $A_N$ can obtain $r_N^* \geq r_n^*$, i.e., $\Tblk{N*} = \min\qty{\Tblk{N}(r_N), \ensuremath{t_{\rm lim}} - \tstart{N*}}$.

On the other hand, if $1 \geq r_n > r^*$, the block duration must be reduced so that $r_n^* = r^*$.
Then, if also $r_N > r^*$, both allocations must be trimmed and are fairly allocated, i.e., $r_n^*=r_N^*=r^* < 1$.
It follows from the properties of \cref{eq:bounded_rstar} that $\tend{N*} = \ensuremath{t_{\rm lim}}$ and $\tstart{N*} > \tstart{N}$.

Finally, if $r_N \leq r^* < r_n \leq 1$, and therefore $r_N < 1$, the properties of \cref{eq:bounded_rstar} imply that $\tend{N}=\ensuremath{t_{\rm lim}}$.
Since $A_N$ cannot be extended without possibly reducing the allocation ratio of other allocations, $\tstart{N*}=\tstart{N}$ and $r_N^* = r_N$.
Since, by assumption, $\tend{n} > \tstart{N}$, the duration of $A_n$ needs to be reduced so that $\Tblk{n*} = \tstart{N*} - \tstart{n} < \Tblk{n}$.

\section{Results}
\label{sec:results}

\begin{figure*}[t!]
\newcommand\fheight{0.5\columnwidth}
\newcommand\fwidth{0.9\columnwidth}

 \begin{subfigure}[b]{\textwidth}
  \centering
  \input{img/scenario1_legend.tex}
 \end{subfigure}
 \\
 \vspace{-2ex}
 \\
 \hspace*{\fill}%
 \begin{subfigure}[t]{0.45\textwidth}
  \centering
  \input{img/accrate_eta_bound_0.1.tex}
  \vspace*{-2ex}
  \caption{Acceptance Rate vs. $\rho$.}
  \label{fig:accrate_rho_bound}
 \end{subfigure}
 \hspace*{\fill}%
 \begin{subfigure}[t]{0.45\textwidth}
  \centering
  \input{img/jain_eta_0.1.tex}
  \caption{Jain's Index vs. $\rho$.}
  \label{fig:jain_rho}
 \end{subfigure}
 \hspace*{\fill}%
 \\
 \vspace*{2ex}
 \\
 \hspace*{-2pt}
 \begin{subfigure}[t]{0.45\textwidth}
  \centering
  \input{img/tblk_eta_0.1.tex}
  \caption{Normalized average $\Tblk{}$ vs. $\rho$.}
  \label{fig:tblk_rho}
 \end{subfigure}
 \hspace*{6ex}%
 \begin{subfigure}[t]{0.45\textwidth}
  \centering
  \input{img/tblk_accrate_0.1.tex}
  \caption{Normalized average $\Tblk{}$ vs. Acceptance Rate.}
  \label{fig:tblk_accrate}
 \end{subfigure}
 \hspace*{\fill}%

 \caption{Results for Scenario~1.}
 \label{fig:results_scenario_1}
\end{figure*}

In this section, we evaluate the algorithms described in \cref{sec:scheduling_algorithms}.
The proposed schedulers have been implemented in Python, only focusing on their capabilities of allocating communication resources to the different traffic streams.
The simulations proposed here do not include full stack behaviors, which will be investigated in future works, and their aim is thus to highlight the fundamental characteristics of each algorithm.

Based on their design criteria, we expected the two algorithms to differ mainly with respect to three \glspl{kpi}, namely the acceptance rate of new requests, the fairness among accepted allocations, and the average scheduled block duration.
To highlight these characteristics we shape the offered traffic based on three parameters, namely: the \textit{average allocation request}
\begin{equation}\label{eq:avg_allocation_request}
T_{\rm avg} = \frac{\Tmin{} + \Tmax{}}{2};
\end{equation}
the \textit{interval ratio}
\begin{equation}\label{eq:interval_ratio}
\rho = \frac{\Tmin{}}{\Tmax{}}  \in [0, 1];
\end{equation}
and the \textit{load factor}
\begin{equation}\label{eq:load_factor}
\lambda = \frac{T_{\rm avg}}{\Tp{}}.
\end{equation}
Note that, for a given average allocation request $T_{\rm avg}$, a low interval ratio $\rho$ corresponds to very flexible allocations, while $\rho=1$ corresponds to rigid allocations where $\Tmin{}=\Tmax{}$.

The proposed algorithms are compared in two different simulation scenarios:
\begin{itemize}
 \item \textit{Scenario~1}: all traffic streams are homogeneous, i.e., all requests have the same parameters.
       Specifically, we consider the case with periodicity $\Tp{}=\frac{\ensuremath{T_{\rm BI}}}{3}$, load factor $\lambda=0.1$, and $\rho \in (0,1)$.
       The impact of different periodicities and load factors is also discussed.
 \item \textit{Scenario~2}: multiple non-homogeneous applications coexist on the same network, thus generating traffic streams with different characteristics.
      We analyze a scenario where traffic streams can be of class $C_1$ or $C_2$, with periodicity $\Tp{C_1}=\frac{\ensuremath{T_{\rm BI}}}{3}$ and $\Tp{C_2}=\frac{\ensuremath{T_{\rm BI}}}{5}$, respectively.
       Both classes have load factor $\lambda=0.1$ and interval ratio $\rho=0.1$.
\end{itemize}

To evaluate the performance of the algorithms, we propose three \glspl{kpi} for \textit{Scenario~1}, shown in \cref{fig:results_scenario_1}.
Note that, given the discrete behavior of the problem and the lack of randomness in the proposed algorithms, the plots cannot be smoothed by running multiple repetitions.

\begin{figure*}[t!]
 \newcommand\fheight{0.5\columnwidth}
 \newcommand\fwidth{0.9\columnwidth}
 %
 \begin{subfigure}[b]{\textwidth}
  \centering
  \input{img/scenario_2_2_legend.tex}
 \end{subfigure}
 \\
 \hspace*{\fill}%
 \begin{subfigure}[b]{0.45\textwidth}
  \centering

  \input{img/conditional_acceptance_c1_3000run_21pt.tex}
  \caption{First allocation is of class $C_1$.}
  \label{fig:conditional_c1}
 \end{subfigure}
 \hspace*{\fill}%
 \begin{subfigure}[b]{0.45\textwidth}
  \centering

  \input{img/conditional_acceptance_c2_3000run_21pt.tex}
  \caption{First allocation is of class $C_2$.}
  \label{fig:conditional_c2}
 \end{subfigure}
 \hspace*{\fill}%
 \\
 \vspace{1ex}%
 \begin{subfigure}[b]{\textwidth}
  \centering
  \input{img/scenario2_legend.tex}
 \end{subfigure}
 \\
 \hspace*{\fill}%
 \begin{subfigure}[b]{0.45\textwidth}
  \centering

  \input{img/variability_3000run_21pt.tex}
  \caption{Variability of scheduled allocations.}
  \label{fig:variability}
 \end{subfigure}
 \hspace*{\fill}%
 \begin{subfigure}[b]{0.45\textwidth}
  \centering

  \input{img/occupancy_3000run_21pt.tex}
  \caption{BI occupancy ratio.}
  \label{fig:occupancy}
 \end{subfigure}
 \hspace*{\fill}%

 \caption{Results for Scenario 2.

  }
 \label{fig:results_scenario_2}
\end{figure*}

The first metric is the acceptance rate (\cref{fig:accrate_rho_bound}), defined as the ratio between the number of accepted allocation requests and the maximum number of acceptable requests.
To compute this achievable upper bound, since all allocations share the same parameters we ignore the strict periodicity assumption and calculate how many allocations with minimum duration $\Tblk{}=\Tmin{}$ can fit in a period $\Tp{}$, which is equal to $N_{\rm max}(\rho) = \left\lfloor \frac{\Tp{}}{\Tmin{}(\rho)} \right\rfloor$ and shown as a red, dashed line.
The acceptance rates can thus be normalized in the interval $[0,1]$, where $1$ means that the scheduler reaches the peak acceptance rate.
Since as $\rho \rightarrow 0$, $\Tmin{} \rightarrow 0$ and thus $N_{\rm max} (\rho) \rightarrow \infty$, we consider at most 100 allocations.

As expected, the \textit{simple scheduler} suffers from a lower acceptance rate than the \textit{max-min fair} one, even though, starting from $\rho=0.5$, the two algorithms tend to behave similarly.
In fact, more rigid allocations do not give enough flexibility to the \textit{max-min fair scheduler} to perform its optimization, thus yielding similar performance to the much \textit{simple scheduler}.

The second metric is \textit{Jain's Fairness Index} (\cref{fig:jain_rho}), defined as:
\begin{equation}
 \mathcal{J} = \frac{\left( \sum_n x_n \right)^2}{n \cdot \sum_n (x_n)^2},
\end{equation}
where only accepted allocations are counted and $x_n$ can either be the block duration $\Tblk{n}$ or the block duration ratio $r_n$.
If the $\{ x_n \}$ are all equal, then $\mathcal{J}(x)=1$.
On the other hand, the more unequal the values of $\{ x_n \}$, the closer the metric to its minimum $\mathcal{J}(x) = \frac{1}{N}$.

Based on the results plotted in \cref{fig:jain_rho}, both algorithms behave fairly with respect to the accepted allocations.

The third metric, shown in \cref{fig:tblk_rho}, offers a different perspective considering the average normalized block duration $\overline{\Tblk{}}/\Tmax{}(\rho)$.
As expected, the \textit{simple scheduler} shows an oscillating trend, due to the discrete behavior of the allocations.
In fact, the last scheduled allocation will only reduce \Tblk{} down to $\Tmin{}(\rho)$, thus, if the portion of DTI left by the previous allocations is less than $\Tmin{}(\rho)$, no additional allocations can be fitted.
On the other hand, the \textit{max-min fair scheduler} will try to reduce all allocations up to their minimum duration in order to avoid rejecting new ones, granting more accepted allocations at the cost of an overall lower block duration.

Finally, the two most discriminating metrics, namely the average normalized block duration and the acceptance rate, are plotted against each other in \cref{fig:tblk_accrate}.
In general, the \textit{simple scheduler} tends to favor higher average block duration for a lower acceptance rate, while the \textit{max-min fair scheduler} tends to favor acceptance rate at the cost of a lower average block duration, as expected.
Both algorithms are able to ensure high fairness to the accepted allocations, generally well above 0.85.

Similar behaviors were also observed for load factors $\lambda \in \{0.025, 0.4\}$, not shown here.
As expected, higher loads tend to have more pronounced variability in both the average block duration and the fairness granted to the accepted allocations.
Curiously, regardless of the load factor, for values of the interval ratio larger than $\rho\approx 0.5$, the two algorithms tend to have very similar performance due to the more rigid allocation requests that do not allow the \textit{max-min fair scheduler} to exploit its agility.

\ml{what about $\Tp{}>\ensuremath{T_{\rm BI}}$? does anything change?} \md{if I recall correctly, we obtained similar results and so we avoid further studies on this case}

\textit{Scenario~2} allows us to analyze the impact of allocations with different periodicities on the overall network performance, as a function of the probability $P(C_1)$ that a request $C_1$ is offered to the system.

Since allocations with different periods coexist, it is mandatory to decide how many allocations should be offered to the schedulers, as this will affect how the acceptance rate is normalized.
We define the \textit{minimum occupancy} as the minimum BI occupancy ratio of the allocation of category $C_i$, namely $O_{\rm min}^i = \Tmin{i}/\Tp{i}$.
Allocations are offered to the schedulers as long as the cumulative minimum offered occupancy does not exceed the value of 1.

In \cref{fig:conditional_c1,fig:conditional_c2} we show the biasing effect of the first accepted allocation on the proposed schedulers.
Clearly, the \textit{simple scheduler} suffers a strong and symmetric effect, meaning that once the first allocation is scheduled with maximum duration, it will be harder for subsequent allocations with a different period to fit the constrained BI, making the scheduler favor allocations with the same period.
On the other hand, it is significantly harder to interpret the behavior of the more complex \textit{max-min fair scheduler}.
From further results, not shown here for lack of space, it is possible to notice that allocations with a lower average BI occupancy are favored, with a slight preference towards those with lower values of \Tp{} and \Tmin{}.
As also shown here, in fact, allocations with lower values of \Tp{}, such as $C_2$ with respect to $C_1$, tend to fragment the BI more, making it harder to then fit allocations with different periodicity and higher \Tmin{}.

To further confirm this biasing behavior, we show the variability $\nu$ among the scheduled allocations, defined as
\begin{equation}
  \nu_{1,2}=\frac{\min\{| C_1 |, | C_2 |\}}{\max\{| C_1 |, | C_2 |\}},
\end{equation}
where $| C_i |$ represents the number of accepted allocations of category $C_i$.
This metric takes values in $[0,1]$, where a value of 0 means that only allocations of a single type have been accepted, while a value of 1 means that the same number of allocations of both categories have been accepted.
Results are shown in \cref{fig:variability}, which confirms that the \textit{simple scheduler} favors a more homogeneous BI allocation, while the \textit{max-min fair scheduler} shows once again more flexibility, being able to accommodate fairly requests from different classes, as shown by the higher variability of the accepted allocations.

Finally, we studied how efficiently the two algorithms are able to use the radio resources by measuring the BI occupancy ratio, i.e., the ratio between scheduled and unscheduled air time.
It can be noticed that, while the \textit{simple scheduler} accepts fewer requests, it is able to use almost all available resources.
This is due to the fact that the scheduler tends to accept homogeneous allocations, allowing them to be packed more efficiently in the BI.
On the other hand, the \textit{max-min fair scheduler} successfully fits multiple allocations of both types, but the constraints on the periodicity and the minimum block duration $\Tmin{}$ prevent it from fully utilizing the whole \gls{bi} when a mixture of the two types of sources is presented.
Nonetheless, it ensures very high occupancy ratios, always above 95\% for the shown example.

\section{Conclusions}
\label{sec:conclusions}

In this paper, we presented a framework for periodic scheduling in \gls{wigig}-compatible devices.
We proposed two heuristic algorithms, \textit{simple} and \textit{max-min fair schedulers}, and accurately described their inner workings.
Finally, we assessed their performance in two different scenarios, showing that the \textit{max-min fair scheduler} tends to trade resource availability for a much higher acceptance rate, contrary to the \textit{simple scheduler}'s behavior, while both schedulers obtained a high Jain's fairness index for the accepted allocations.

Even working in a simplified settings without considering further sources of complexity from other parts of the communication stack, it was possible to notice that both the design and the evaluation of \gls{wigig}-specific scheduling algorithms for periodic sources is highly non-trivial and can show surprising results.
In fact, while the formalization of the problem is straightforward, scheduling algorithms often have to deal with many hard-to-predict edge cases, which greatly increases the difficulty of designing a general algorithm.

Future works will focus on multiple objectives.
A first objective is to implement these algorithms in a full-stack simulator, allowing the study of APP-layer performance metrics for a range of possible applications.
A second objective is to extend the framework by relaxing some of the assumptions made in \cref{sec:framework_description} and test the impact of each one of them on the overall performance of a \gls{wigig} system.
A third objective is to extend the study for scheduling multiple allocations at once, a possibility given by \gls{mimo} techniques.
To further achieve a realistic evaluation, the impact of accurate wireless channel simulations, physical layer design, as well as user mobility will also be studied.

\bibliographystyle{IEEEtran}