Estimating the Probability of Meeting a Deadline \\ in Hierarchical Plans
Liat Cohen
Solomon Eyal Shimony
Gera Weiss
Liat Cohen
Solomon Eyal Shimony
Gera Weiss
Computer Science Department
Ben Gurion University of The Negev
Beer-Sheva
Israel 84105
\{liati
shimony
geraw\
Given a hierarchical plan (or schedule) with uncertain task times, we may need to determine the probability that a given plan will satisfy a given deadline. This problem is shown to be NP-hard for series-parallel hierarchies. We provide a polynomial-time approximation algorithm for it. Computing the expected makespan of an hierarchical plan is also shown to be NP-hard. We examine the approximation bounds empirically and demonstrate where our scheme is superior to sampling and to exact computation.%, using some task networks.
deadline_arxiv.tex



\documentclass{article}


% Use the postscript times font!
\usepackage{times}

% the following package is optional:
\usepackage{latexsym} 

\usepackage{multirow}



%   used packages
%----------------------
\usepackage{todonotes}

\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{mathrsfs}

\usepackage{amsfonts}
\usepackage{graphicx}
\usepackage{framed}
\usepackage{tikz}
\usetikzlibrary{arrows}
\usepackage{amsthm} 
\usepackage[ruled,vlined,linesnumbered]{algorithm2e} 
\usepackage{makecell}
\usepackage{tikz}
\newcommand\diag[4]{%
  \multicolumn{1}{p{#2}|}{\hskip-\tabcolsep
  $\vcenter{\begin{tikzpicture}[baseline=0,anchor=south west,inner sep=#1]
  \path[use as bounding box] (0,0) rectangle (#2+2\tabcolsep,\baselineskip);
  \node[minimum width={#2+2\tabcolsep},minimum height=\baselineskip+\extrarowheight] (box) {};
  \draw (box.north west) -- (box.south east);
  \node[anchor=south west] at (box.south west) {#3};
  \node[anchor=north east] at (box.north east) {#4};
 \end{tikzpicture}}$\hskip-\tabcolsep}}
 
\newtheorem{theorem}{Theorem}
\newtheorem{claim}{Claim}
\newtheorem{observation}{Observation}
\newtheorem{lemma}{Lemma}
\newtheorem{proposition}{Proposition}
\newtheorem{definition}{Definition}
\newtheorem{example}{Example}

\newcommand\Tau{\tau}

\DeclareMathOperator{\support}{support}
\DeclareMathOperator{\Trim}{Trim}

\pdfinfo{
/Title{Estimating Utility of Hierarchical Plans}
/Author{}
}
\begin{document}




%\title{Hierarchical plans: Computing the Probability of Meeting a Deadline}
\title{Estimating the Probability of Meeting a Deadline \\ in Hierarchical Plans}

%\author{Liat Cohen \and Solomon Eyal Shimony \and Gera Weiss}

\author{Liat Cohen \and Solomon Eyal Shimony \and Gera Weiss \\
Computer Science Department\\
Ben Gurion University of The Negev\\
Beer-Sheva, Israel 84105\\
\{liati,shimony,geraw\}@cs.bgu.ac.il
}



\maketitle

\begin{abstract}
Given a hierarchical plan (or schedule) with uncertain task times, we may need to determine
the probability that a given plan will satisfy a given deadline.
This problem is shown to be NP-hard for series-parallel hierarchies. We provide a polynomial-time approximation algorithm for it. 
Computing the expected makespan of an hierarchical plan is also shown to be NP-hard.
We examine the approximation bounds empirically and demonstrate 
where our scheme is superior to sampling and to exact computation.%, using some task networks.

\end{abstract}  

\section{Introduction}

Numerous planning tools produce plans that call for executing tasks non-linearly.
Usually, such plans are represented as a tree, where the
leaves indicate primitive tasks, and other nodes represent compound tasks consisting of executing
their sub-tasks either in parallel (also called ``concurrent'' tasks \cite{gabaldon2002programming}) or in sequence.
\cite{erol1994htn,nau1998control,nau2003shop2,kelly2008offline}.

Given such a hierarchical plan representation, it is frequently of interest to evaluate its desirability in terms of
resource consumption, such as fuel, cost, or time. The answer to such questions can be used to
decide which of a set of plans, all valid as far as achieving the goal(s) are
concerned, is better given a user-specified utility function. Another reason to compute these
distributions is to support runtime monitoring of resources, generating alerts to the execution
software or human operator if resource consumption in practice
has a high probability of surpassing a given threshold.

While most tools aim at good average performance of the plan, in which case one may ignore the 
full distribution and consider only the expected resource consumption~\cite{bonfietti2014disregarding}, our paper focuses
on providing guarantees for the probability of meeting deadlines. This type of analysis is needed, e.g., in
Service-Level-Agreements (SLA) where guarantees of the form: ``response time less than 1mSec in at least 95\% of the cases'' are common~\cite{buyya2011sla}.

We assume that an hierarchical plan is given in the form of a tree, with uncertain resource consumption of the primitive actions
in the network, provided as a probability distribution. The problem is to compute
a property of interest of the distribution for the entire task network.
In this paper, we focus mainly on the issue of computing the probability $P(t<T)$ of satisfying a deadline $T$
(i.e. that the {\em makespan} $t$ of the plan is less than a given value).
Since in the above-mentioned applications for these computations, one needs results in real-time (for monitoring) or
multiple such computations (in comparing candidate plans), efficient computation here is crucial, and is more
important than in, e.g., off-line planning.

We show that computing $P(t<T)$ is NP-hard even for the
simple sum of independent random variables (r.v.s), the first contribution of this paper. A deterministic polynomial-time approximation scheme for
this problem is proposed, the second contribution of this paper. Error bounds are analyzed and are shown to be tight.
For discrete r.v.s with finite support, finding the distribution of the maximum can be done in low-order polynomial time.
However, when compounded with errors generated due to approximation in subtrees, handling this case requires careful
analysis of the resulting error.
The approximations developed for both sequence and parallel nodes are combined into an overall
algorithm for task trees, with an analysis of the resulting error bounds, yielding a polynomial-time (additive error) approximation scheme for computing the
probability of satisfying a deadline for the complete network, another contribution of this paper. 

We also briefly consider computing {\em expected} makespan. Since for discrete r.v.s, in parallel nodes one can compute an
exact distribution efficiently, it is easy to compute an {\em expected makespan} in this case as well as for sequence nodes. Despite that, 
we show that for trees with both parallel
and sequence nodes, computing the expected makespan is NP-hard.

Experiments are provided in order to examine the quality of approximation in practice when compared to the theoretical
error bounds. A simple sampling scheme is also provided as a yardstick, even though the sampling does not come with
error guarantees, but only bounds in probability. 
%This is done for
%randomly generated task networks, as well as some task networks taken from the (blinded for review) team task descriptions implementation for
%some of the DARPA robotics (simulation phase) challenge scenarios, and some task networks generated by the JSHOP2 planner
%for the package delivery problem.
Finally, we examine our results in light of related work in the fields of planning and scheduling, as well as
probabilistic reasoning. %The related work suggests some directions in which this work can be extended.


\section{Problem statement}\label{sec:formal}

We are given a hierarchical plan represented as a task tree consisting of three types of nodes: primitive actions as leaves, sequence nodes, and parallel nodes.
Primitive action nodes contain distributions over their resource consumption. 
Although any other resource can be represented, we will assume henceforth,
in order to be more concrete, that the only resource of interest is time.
A sequence node $v_s$ denotes a task that has been decomposed
into subtasks, represented by the children of $v_s$, which must be executed in sequence 
in order to execute $v_s$.
We assume that a subtask of $v_s$ begins as soon as its predecessor in the sequence terminates.
Task $v_s$ terminates when its last subtask terminates.
A parallel node $v_p$ also denotes a decomposed task, but subtasks begin execution in parallel
immediately when task $v_p$ begins execution;  $v_p$ terminates
as soon as all of the children of $v_p$ terminate.

Resource consumption is uncertain, and described as probability distributions in the leaf nodes.
We assume  that the distributions are  independent (but {\em not} necessarily identically distributed).
We also assume initially that the r.v.s are discrete and have finite support (i.e. number of values
for which the probability is non-zero).
As the resource of interest is assumed to be completion time, let each leaf node $v$ have
a completion-time distribution $P_v$, in some cases represented as a cumulative distribution function form (CDF) $F_v$. 

The main computational problem tackled in this paper is: Given a task tree $\tau $ and a deadline $T$, 
compute the probability that $\tau $ satisfies the deadline $T$ (i.e. terminate in time $t \leq T$). 
We show that this problem is NP-hard and provide an approximation algorithm.
The above {\em deadline problem} reflects a step utility function: a constant positive utility $U$ for all $t$ less than or equal to a deadline 
time $T$, and $0$ for all $t>T$. We also briefly consider a linear utility function, requiring
computation of the expected completion time of $\tau$ and show that this {\em expectation problem} is also NP-hard.


\section{Sequence nodes}\label{sec:seq}

Since the makespan of a sequence node is the sum of durations of its components,
the deadline problem on sequence nodes entails computation of (part of) the CDF of a sum of r.v.s which is an NP-hard problem (shown in Section \ref{sec:complexity}). Thus, there is a need for an approximation algorithm for sequence nodes. The main ingredient in our approximation scheme is the $\Trim$ operator specified as follows:

\begin{definition}[The $\Trim$ operator]
For a discrete r.v. $X$ and a parameter $\varepsilon > 0$, 
consider the sequence of elements in the support
of $X$ defined recursively by:
$x_1 = \min \operatorname{support}(X)$
and, if the set $\{ x >x_i\colon Pr( x_i < X \leq x) > \varepsilon  \}$ is not empty, 
$x_{i+1} = \min \{ x>x_i \colon Pr( x_i < X \leq x) > \varepsilon  \}.$
Let $l$ be the length of this sequence, i.e., let $l$ be the first index for which the set above is empty. For notational convenience, define $x_{l+1}=\infty$. Now, define $X'=\Trim(X,\varepsilon)$ to be the random variable specified  by:
$$Pr(X'{=}x):= \begin{cases}
 Pr(x_{i}{\leq} X {<} x_{i+1}) & \text{if } x{=}x_i{\in}\{x_1,\dots,x_l\}; \\
 0 & \text{otherwise.}
 \end{cases}
$$
\end{definition}
For example, if $X$ is a r.v. such that $Pr(X{=}1){=}0.1, Pr(X{=}2){=}0.1$
and $Pr(X{=}4){=}0.8$, the r.v. $X'{=}\Trim(X,0.5)$ is given by
$Pr(X'{=}1){=}0.2$ and $Pr(X'{=}4){=}0.8$. In words, the $\Trim$ operator removes 
consecutive domain values whose accumulated probability is
less than $\varepsilon$ and adds their probability mass to the element in the support that precedes them.

Using the $\Trim$ operator, we are now ready to introduce the main operator of this section:
%\begin{definition}%[The $\operatorname{Sequence}$ operator]
%\begin{equation*}
%\begin{split}
%\operatorname{Sequence}(X_1&,\dots,X_n,\varepsilon) :=\\ &\Trim(\operatorname{Sequence}(X_1,\dots,X_{n-1},\varepsilon)+X_n, \varepsilon )
%\end{split}
%\end{equation*}
%and $\operatorname{Sequence}(X_1,\varepsilon)=\Trim(X_1,\varepsilon)$.
%\end{definition}

\begin{definition}
Let $\operatorname{Sequence}(X_1,\dots,X_n,\varepsilon)$ be $\Trim(\operatorname{Sequence}(X_1,\dots,X_{n-1},\varepsilon){+}X_n, \varepsilon )$
and $\operatorname{Sequence}(X_1,\varepsilon){=}\Trim(X_1,\varepsilon)$.
\end{definition}


$\operatorname{Sequence}$ takes a set of r.v.s and computes a r.v. that represents an approximation of
their sum by applying the $\Trim$ operator after adding each of the variables. The parameter $\varepsilon$, (see below)
specifies the accuracy of approximation.


%As the dealine problem, computing $Pr(\sum_{i=0}^{n} Y_{i}<T)$, is NP-hard, we next present an approximation algorithm for this problem.
%In this paper we represent random variables as an explicit mapping of values to probabilities, such that the sum of the probabilities over all values is 1.
%Each of the $Y_i$, with $1\leq i \leq n$, is a nonegative discrete random variable with finite support assumed to be bounded by $m$.

%For the specification of the algorithm we will use the following two notations:
%\begin{itemize}
%\item For a (discrete) random variable $X$, we denote $\support(X)=\{ x \in \mathbb{R} \colon Pr(X=x) \neq 0\}$. 
%\item A distribution of a random variable $X$ is represented by a list $D_X$ of $(x, p)$ pairs, where $x\in \support(X)$
%and $p$ is the probability $Pr(X=x)$. We assume that $D_X$ is always sorted in increasing order of $x$.
%$D_X[x]$ is used to denote $Pr(X=x)$, and $\support(D_X)$ denotes a sorted list of the values in $\support(X)$.
%\end{itemize}


The $\operatorname{Sequence}$ operator can be implemented by the procedure outlined in Algorithm~\ref{alg:sequence}. The algorithm computes
the distribution $\sum_{i=0}^{n} X_{i}$ using convolution (the Convolve() operator in line 3) in a straightforward
manner. Computing the convolution is itself straightforward and not further discussed here.
However, since the support of the resulting distribution may be exponential in the number of convolution operations, the algorithm
must trim the distribution representation to avoid this exponential blow-up. This decreases the support size, while introducing error. The trick is
to keep the support size under control, while making sure that the error does not increase beyond a desired tolerance.
Note that the size of the support can also be decreased by simple ``binning'' schemes, but these do not provide the desired guarantees.
In the algorithm, the PDF of a r.v. $X_i$ is represented by the list $D_{X_i}$, which consists of $(x, p')$ pairs, where $x\in \support(X_i)$
and $p'$ is the probability $Pr(X_i{=}x)$, the latter denoted by $D_{X_i}[x]$.
We assume that $D_{X_i}$ is kept sorted in increasing order of $x$. 

%%%%%%%%%%%%%%%%%%%%%%%%NEED TO ADD
%
\begin{algorithm}
  \DontPrintSemicolon
  \SetKwFunction{Sequence}{Sequence} 
  \SetKwFunction{Convolv}{Conv}
  \SetKwFunction{Trim}{Trim}
  %\SetKwProg{myalg}{Algorithm}{}{}
  %-----------------------------------%
  %\myalg{\Sequence{$D_{Y_1},\dots,D_{Y_n}$ , $T$, $\varepsilon$}}{
  %$k=n/\varepsilon$\;
   $D=((0,1))$ //  Dummy random var.: $0$ with prob. $1$ \;
  
  \For{$i=1$ \emph{\KwTo} $n$} {
	$D$=Convolve($D, D_{X_i}$)\; 
	$D$=\Trim{$D, \varepsilon$}\;
	}
  
  \Return $D$
  
  %-------------------------------%
  \SetKwProg{myproc}{Procedure}{}{}
  \myproc{
  \Trim{$D$,$\varepsilon$}
  }{ 
             $D'=()$, $p=0$ \;	  
  		$d_{0} = d_{prev} = \min \support(D)$\;
		\ForEach{$d \in  \support(D) \setminus \{d_{0}\}$ in ascending order}{			 	
			\eIf{$p+D[d]\leq \varepsilon$}{
				$p=p+D[d]$\;
			}
			{ Append $(d_{prev}, D[d_{prev}]+p)$ to $D'$\;
			  $d_{prev} = d$, $p=0$\;
			}
                 }
       		 Append $(d_{prev}, D[d_{prev}]+p)$ to  $D'$\;
  \Return $D'$}
     
\caption{Sequence ($X_1,\dots,X_n$ , $\varepsilon$)}  
\label{alg:sequence}
\end{algorithm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

We proceed to show that Algorithm~\ref{alg:sequence} indeed approximates the sum of the r.v.s, and
analyze its accuracy/efficiency trade-off.
A notion of approximation relevant to deadlines is:

\begin{definition}
For r.v.s $X'$, $X$, and $\varepsilon\in[0,1]$  we write $X'\approx_\varepsilon X$  if  $0 \leq Pr(X' \leq T) - Pr(X \leq T) \leq \varepsilon$ for all $T>0$.
\end{definition}

Note that this definition is asymmetric, because, as shown below,
our algorithm never underestimates the exact value.
%This notion of approximation is useful in the context of estimating the probabilities of missing deadlines in the following sense. 
%If, for example, $X'\approx_\varepsilon X$ and we use the approximation $X'$ instead of an accurate variable $X$ to compute 
%the probability that the value of the variable exceeds some threshold $T$, we miss the true answer (the one that we would have got if we used the accurate variable) 
%by at most $\varepsilon$. The motivation for this is, of course, that it may sometimes be easier (in terms of computation times or memory) 
%to compute these probabilities for $X'$ than it is for $X$.
For the proof of the Lemma~\ref{Convolv} below, we establish the following technical claim (can be proven by induction on $n$): 
\begin{claim}
Let $(a_i)_{i=1}^n$ and $(b_i)_{i=1}^n$ be sequences of real numbers such that $\sum_{i=1}^k a_i \geq 0$ for all $1
\leq k \leq n$ and $b_1 \geq b_2 \geq \cdots \geq b_n \geq 0$, then $\sum_{i=1}^n a_i b_i \geq 0$
\label{pos-conv}
\end{claim}

We now bound the approximation error of sums of r.v.s:

\begin{lemma} \label{Convolv}
For discrete r.v.s $X_1,X_1',X_2,X_2'$ and $\varepsilon_1,\varepsilon_2 \in [0,1]$, if $X_1'\approx_{\varepsilon_1} X_1$ and  $X_2'\approx_{\varepsilon_2} X_2$,
then $X_1'+X_2'\approx_{\varepsilon_1+\varepsilon_2} X_1+X_2$.
\end{lemma}

\begin{proof}
{\small 
$Pr(X_1'+X_2' \leq T)  -  Pr(X_1+X_2 \leq T)$
\begin{align*}
%&Pr(X_1'+X_2' \leq T)  -  Pr(X_1+X_2 \leq T)  \\
&{=}\sum_{j=1}^{T}\!Pr(X_1'{=}j)\!\!\! \underbrace{Pr(X_2' {\leq} T{-}j)}_{\leq Pr(X_2 {\leq} T{-}j) + 
\varepsilon_2} \!\!\!- \sum_{j=1}^{T} Pr(X_1{=}j) Pr(X_2 {\leq} T{-}j) \\
%
%&{\leq}\! \sum_{j=1}^{T} Pr(X_1'{=}j)\!(Pr(X_2{\leq} T{-}j){+}\varepsilon_2){-} Pr(X_1{=}j)\!Pr(X_2{\leq} T{-}j) \\
%&{=}
&{\leq}\! \sum_{j=1}^{T}\! (Pr(X_1'{=}j){-}Pr(X_1{=}j))\!\underbrace{Pr(X_2{\leq} T{-}j)}_{\in [0,1]} {+} \varepsilon_2\!\sum_{j=1}^{T}\! Pr(X_1'{=}j)  \\
&{\leq}\! \underbrace{\sum_{j=1}^{T} (Pr(X_1'{=}j)-Pr(X_1{=}j))}_{\leq \varepsilon_1} + \varepsilon_2 \underbrace{\sum_{j{=}1}^{T} Pr(X_1'{=}j)}_{\in [0,1]} 
\leq\varepsilon_1 + \varepsilon_2.
\end{align*}}

Finally, we show that the  difference is also nonnegative:%between $Pr(X_1'+X_2' \leq T)$ and $Pr(X_1+X_2 \leq T)$ is nonnegative:
{ \small
\begin{align*}
&Pr(X_1'+X_2' \leq T)  -  Pr(X_1+X_2 \leq T)  \\
%
&{=}\sum_{j=1}^{T}\!Pr(X_1'{=}j)Pr(X_2' {\leq} T{-}j) {-} \sum_{j{=}1}^{T} Pr(X_1{=}j) Pr(X_2 {\leq} T{-}j) \\
%
&{=}\sum_{j=1}^{T}\!(Pr(X_1'{=}j) - Pr(X_1{=}j))Pr(X_2{\leq} T{-}j) \\[-.35cm]
      &\hspace{1.5cm} {+}  \sum_{j=1}^{T} Pr(X_1'{=}j) (Pr(X_2' {\leq} T{-}j) - Pr(X_2 {\leq} T{-}j))  
\end{align*}}
The first term here is non-negative by Claim~\ref{pos-conv}, the second is nonnegative because it
is a sum of nonnegative numbers. 
\end{proof}
\noindent We now show that $\Trim(X,\varepsilon)$ is an $\varepsilon$-approximation of $X$:
\begin{lemma} \label{Trim}
$\Trim(X,\varepsilon)\approx_\varepsilon X$

\end{lemma}
\begin{proof}
Let $X'=\Trim(X,\varepsilon)$. 
Let $x_1{<}\cdots{<}x_m$ be the support of $X'$ and let $l=\max\{ i \colon x_i\leq T\}$. We have, 
\begin{equation}
Pr(X'=x_i) = Pr(x_i \leq X < x_{i+1})
\label{collect}
\end{equation}
because, after $\Trim$, the probabilities of elements that were removed from the support are assigned to the element that precedes them.  
From Equation~\eqref{collect} we get:
{\small
\begin{align*}
&Pr(X' \leq T ) - Pr(X \leq T)  \\%[.25cm]
 &{=} \sum_{i=0}^{l-1} (Pr(X'{=}x_i)- Pr(x_i\leq X <x_{i+1})) \\[-.2cm]
& \hspace{2.7cm}+(Pr(X'{=}{x_l})-Pr(x_l \leq X \leq T)) \\%[.25cm]
&{=}Pr(X'{=}{x_l}) {-} (Pr(x_l {\leq} X {<} x_{l+1}) {-} Pr(T {<} X {<} x_{l+1})) \\
&{=}Pr(T < X < x_{l+1}) \in (0,\varepsilon]
\end{align*}}
The inequality $Pr(T < X < x_{l+1}) \leq \varepsilon$ follows from the observation that, 
for all $i$, $Pr(x_i < X < x_{i+1}) < \varepsilon$, because $p$ is never greater than $\varepsilon$
in Algorithm~\ref{alg:sequence}.
\end{proof}



To bound the amount of memory needed for our approximation algorithm,
the next lemma bounds the size of the support of the trimmed r.v.:
\vspace*{-0.1cm}

\begin{lemma} \label{SizeD}
$|\operatorname{support}(\Trim(X,\varepsilon))| \leq 1/\varepsilon$
\end{lemma}

\vspace*{-0.1cm} 
\begin{proof}
Let $X'=\Trim(X,\varepsilon)$ and let $x_1<\cdots<x_m$ be the support of $X'$. And, for notational convenience, let $x_{m+1}=\infty$.
%From the definition of random variable,  $\sum_{x=x_1}^{x_m} Pr(X'{=}x)=1$. 
% Gera: removed the above line because it is trivial.
Let $p_i =\sum_{x_i<x<x_{i+1}} Pr(X{=}x)$. Then, $1 = \sum_{i=1}^{m} Pr(X'{=}x_i) =Pr(X{=}x_1)+\sum_{i=1}^{m-1} (p_i+ Pr(X{=}x_{i+1}))+p_m$.
According to algorithm~\ref{alg:sequence}, lines 11-12, $p_i \leq \varepsilon $ and $p_i+Pr(X{=}x_{i+1})>\varepsilon$ for all $0 \leq i < m$. Therefore, $1=\sum_{i=1}^{m} Pr(X'{=}x_i) > Pr(X'{=}x_1) + (m-1)\cdot \varepsilon + p_m $. 
%By direct algebraic manipulation, 
Using the fact $Pr(X'{=}x_1)>0$, we get: $(m-1)\cdot \varepsilon<1\Rightarrow m<1/\varepsilon+1 \Rightarrow m\leq 1/\varepsilon$. 
\end{proof}

\vspace*{-0.2cm} 
These lemmas highlight the main idea behind our approximation algorithm: the $\Trim$ operator trades off approximation error 
for a reduced size of the support. 
The fact that this trade-off is linear allows us to get a linear approximation error in polynomial time, as shown below:

\vspace*{-0.1cm} 

\begin{theorem}
If $X_i' \approx_{\varepsilon_i} X_i$ for all $i \in\{1,\dots,n\}$ and $\hat{X}= Sequence(X_1',\dots,X_n', \varepsilon)$ then
$\hat{X} \approx_e \sum_{i=1}^{n} X_{i}$, where $e={\sum_{i=1}^n\varepsilon_i + n \varepsilon}$. 
\label{appSeqTheorem}
\end{theorem} 
%\begin{theorem}
%If $X_i' \approx_{\varepsilon_i} X_i$ for all $i \in\{1,\dots,n\}$ and $\hat{X}= Sequence(X_1',\dots,X_n', \varepsilon/n)$ then
%$\hat{X} \approx_e \sum_{i=1}^{n} X_{i}$, where $e={\sum_{i=1}^n\varepsilon_i + \varepsilon}$. 
%\label{appSeqTheorem}
%\end{theorem} 
\vspace*{-0.2cm} 

\begin{proof} (outline) For $n$ iterations, from Lemma~\ref{Convolv}, we get an accumulated error of 
$\varepsilon_1 +\dots+ \varepsilon_n$. From Lemma~\ref{Trim}, we get an additional error of at most $n\varepsilon$ due to trimming. 
% From Lemma~\ref{SizeD}, $|D|\leq 2n/\varepsilon$ we get overall run-time of $O(\frac{1}{\varepsilon} \cdot n^2 \cdot m)$.
\end{proof}

\vspace*{-0.3cm} 

\begin{theorem} \label{appSeqComplexTheorem}
Assuming that $m \leq 1/\varepsilon$, the procedure
$Sequence(X_1',\dots,X_n', \varepsilon)$ can be computed
in time $O((nm/\varepsilon)\log(m/\varepsilon))$ using $O(m/\varepsilon)$ memory, where $m$ is the size of the largest support of any of the $X'_i$s.
\end{theorem} 

\vspace*{-0.3cm} 
\begin{proof}
From Lemma~\ref{SizeD}, the size of list $D$ in Algorithm \ref{alg:sequence} is at most $m/\varepsilon$
just after the convolution, after which it is trimmed, so the space complexity is $O(m/\varepsilon)$.
$\operatorname{Convolve}$ thus takes time $O((m/\varepsilon)\log(m/\varepsilon))$, where the logarithmic factor is required internally
for sorting. Since the runtime of  the $\Trim$ operator is linear, and the
outer loop iterates $n$ times, the overall 
run-time of the algorithm is $O((n m/\varepsilon) \log(m/\varepsilon))$.
\end{proof}

%To show that the error bound provided in Theorem~\ref{appSeqTheorem} is tight, we
%demonstrate by the following example 
%that there are $n$ random variables $X_1,\dots,X_n$ for which $Sequence(X_1,\dots,X_n,\varepsilon/n)$ 
%indeed results in error $\varepsilon$.
\vspace*{-0.2cm} 


\begin{example}
The error bound provided in Theorem~\ref{appSeqTheorem} is tight, i.e.
$Sequence(X_1,\dots,X_n,\varepsilon/n)$ may
results in error $\varepsilon$:
Let $0 {\leq} \varepsilon {<}1$ and $n {\in} \mathbb{N}$ such that $1{-}\varepsilon {>} {\varepsilon }/{n}$.
% i.e.,  $\varepsilon$ is small or $n$ is large. 
Consider, for very small $\delta>0$, the r.v. $X_1$ defined by:
$$
Pr(X_1{=}x) {=} \begin{cases}
\delta  & x=0, \\
{\varepsilon }/{(n (1-\delta )^{x})} & x\in\{1,\dots,n\}, \\
1{-}\delta{-}\sum_{x=1}^n \frac{\varepsilon }{n(1-\delta )^{x}} & x=n+1, \\
0 & \text{otherwise}
\end{cases}
$$
and, for $i=2,\dots,n$, let the r.v.s $X_i$ be such that $Pr(X_i{=}0){=}1{-}\delta$, $Pr(X_i{=}n^2){=}\delta$, and zero otherwise.
%$$
%Pr(X_i{=}x) =\begin{cases}
%1-\delta  & x=0, \\
%\delta     & x=n^2, \\
%0 & \text{otherwise}
%\end{cases}
%$$
%The distribution of $X=X_1+X_2$ is
%$$
%Pr(X{=}x) {=} 
%\begin{cases}
%\delta(1-\delta)  & x=0, \\
%{\varepsilon }/{n} & x=1, \\
%{\varepsilon }/{(n (1-\delta )^{x-1})} & x\in\{2,\dots,n\}, \\
%(1{-}\delta)Pr(X_1{=}n{+}1)& x=n+1, \\
%\delta Pr(X_1{=}x{-}n^2)& n^2 {\leq} x {\leq} n^2{+}n{+}1 \\
%%\delta Pr(X_1{=}x{-}n^2)& x{\in}\{n^2, \dots,n^2{+}n{+}1\}, \\
%0 & \text{otherwise}
%\end{cases}
%$$
%The idea here is that the convolution with $X_2$ results in a random variable that is similar in ``shape'' to $X_1$, 
%if we ignore numbers that tend to zero as $\delta$ approaches zero. The convolution also
%modifies the probability $Pr(X{=}1)$ from slightly greater than  $\varepsilon/n$ to precisely $\varepsilon/n$, which will then allow it
%to be trimmed.
%
%Then, if we apply $\Trim(X_1+X_2,\varepsilon/n)$, when $\delta$ is sufficiently small, we get the random variable $X'$ whose probability distribution is:
%$$
%Pr(X'{=}x) {=} 
%\begin{cases}
%\delta(1-\delta)+{\varepsilon }/{n}  & x=0, \\
%{\varepsilon }/{(n (1-\delta )^{x-1})} & x\in\{2,\dots,n\}, \\
%1-Pr(X'{<}n+1)& x=n+1, \\
%0 & \text{otherwise.}
%\end{cases}
%$$ 
%
%Note that indeed trimming shifts the mass from $Pr(X=1)=\varepsilon/n$ to $Pr(X'=0)$.
%This repeats in all steps so, after $n$ steps, we get a random variable $X'$ such that $Pr(X'{=}0) \xrightarrow{\,\,\delta \to 0\,\,} \varepsilon$. Therefore, 
%$Pr(Sequence(X_1,\dots,X_n,\varepsilon/n) {\leq} 0) {-} Pr(X_1{+}\dots{+}X_n {\leq} 0)$
%%\begin{align*}
%%&Pr(Sequence(X_1,\dots,X_n,\varepsilon/n) {\leq} 0)  \\
%%&             \hspace{4cm} - Pr(X_1+\dots+X_n \leq 0)
%%\end{align*}
%approaches $\varepsilon$ as $\delta$ approaches zero
%which means that there exist no $\varepsilon' {<} \varepsilon$
%such that $Sequence(X_1,\dots,X_n,\varepsilon/n) \approx_{\varepsilon'} X_1 + \cdots + X_n$ for all $\delta >0$.



\end{example}
\vspace*{-.4cm} 

\section{Parallel nodes}\label{sec:par}

\SetKwFunction{Parallel}{Parallel}

Unlike sequence composition, the deadline problem for parallel composition is easy to compute, since the execution time of
a parallel composition is the maximum of the durations:
\vspace*{-0.2cm}
{\small
\begin{align}
%&Pr(\max\{X_1,\dots,X_n\}\leq T)= \nonumber \\
Pr(\max_{i\in[1:n]}X_i\leq T)
%&=Pr(X_1\leq T \wedge\dots\wedge X_n\leq T)
{=}Pr(\bigwedge_{i=1}^n X_i\leq T) 
{=}\prod_{i=1}^n Pr(X_i\leq T)
\label{eq:par}
\end{align}}

\noindent where the last equality follows from independence of the r.v.s.
We denote the construction of the CDF using Equation~\eqref{eq:par} by $\Parallel(X_1,\dots , X_n)$.
If the r.v.s are all discrete with finite support, $\Parallel(X_1,\dots , X_n)$
incurs linear space, and computation time $O(nm log(n))$.

If the task tree consists only of parallel nodes, one can
compute the exact CDF, with the same overall runtime.
However, when the task tree contain both sequence and parallel nodes we may get
only approximate CDFs as input, and now the above straightforward computation can compound the errors.
When the input CDFs are themselves approximations, we bound the resulting error:

\begin{lemma} \label{appPalTheorem}
For discrete r.v.s $X_1', \dots X_n'$, $X_1, \dots,X_n$, if for all $i=1,\dots,n$,  $X_i' \approx_{\varepsilon_i} X_i$ and $0\leq\varepsilon_i\leq  \frac{1}{n (K n+1)}$ for some $K>0$,
then, for any $\varepsilon \geq \varepsilon_i$, we have: $\max_{i\in[1:n]}X_i' \approx_{e} \max_{i\in[1:n]}X_i$ where $e=\sum_{i=1}^n \varepsilon_i + \varepsilon/K$.
\end{lemma}
\begin{proof}
%$Pr(\max\{X_1',\dots, X_n'\} {\leq} T) {-} Pr(\max\{X_1,\dots, X_n\} {\leq} T)$
$Pr(\max_{i\in[1:n]}X'_i
 \leq T) {-} Pr(\max_{i\in[1:n]}X_i
 \leq T)$
{\small
\begin{align*} 
%&Pr(\max\{X_1',\dots, X_n'\} {\leq} T) {-} Pr(\max\{X_1,\dots, X_n\} {\leq} T)\\
%&=\prod_{i=1}^n Pr(X_i'\leq T)-\prod_{i=1}^n Pr(X_i \leq T)\\
&\leq \prod_{i=1}^n (Pr(X_i\leq T)+\varepsilon_i) {-} \prod_{i=1}^n Pr(X_i \leq T) \\
&\leq \prod_{i=1}^n (1+\varepsilon_i) - 1 
\leq 1+\sum_{i=1}^n \varepsilon_i + \sum_{k=2}^n {n \choose k}  \varepsilon^k - 1 \\
&\leq \sum_{i=1}^n \varepsilon_i + \!\!\!\!\underbrace{\sum_{k=2}^n n^k  \varepsilon^k}_{\text{sum of a geo. series}}\!\!\!\!
\leq \sum_{i=1}^n \varepsilon_i + \frac{n^2 \varepsilon^2}{1-n \varepsilon}
\leq \sum_{i=1}^n \varepsilon_i + \varepsilon/K\\
%\leq (1+\varepsilon)^n  {-} 1 \\
%&= \sum_{k=0}^{n}{n\choose k}\varepsilon ^k - 1  
%\leq \!\!\!\underbrace{\sum_{k=1}^{n} n^k \varepsilon ^k}_{\text{sum of geo. series}} 
%\!\!\!< \dfrac{n\varepsilon}{1-n\varepsilon}
%\leq (1+n)\varepsilon
\end{align*}}
\vspace*{-0.7cm} 

Since $Pr(X_i'\leq T) > Pr(X_i \leq T)$ for each $i$, this expression is nonnegative.
\end{proof}

%Lemma~\ref{appPalTheorem} will allow us to use, for parallel nodes, the efficient computation indicated in Theorem~\ref{TPar}
%provided only that the approximation of the children are good enough, as elaborated in the next section.
\vspace*{-0.5cm} 
\section{Task trees: mixed sequence/parallel}

Given a task tree $\tau$ and a accuracy requirement $0<\varepsilon<1$, we generate a distribution for a r.v. $X_{\tau}' $
approximating the true duration distribution  $X_{\tau}$ for the task tree. 
We introduce the algorithm and prove that the algorithm indeed returns an $\varepsilon$-approximation of the completion time of the plan. 
For a node $v$, let $\tau_v $ be the sub tree with $v$ as root and let $child_v$ be the set of children of $v$. 
We use the notation $|\tau|$ to denote the total number of nodes in $\tau$.

\begin{algorithm}
  \DontPrintSemicolon
  \SetKwFunction{Network}{Network} 
  \SetKwFunction{Sequence}{Sequence}
  \SetKwFunction{Parallel}{Parallel}
  
  %-----------------------------------%
  
  {

 Let $v$ be the root of $\tau$  // Hence, $\tau_v=\tau$ \;
	$n_v = |child_v|$ \;
  \If{$v$ is a $\operatorname{Primitive}$ node}{
	\Return the distribution of $v$\;
			}
	
  \If{$v$  is a $\operatorname{Sequence}$ node}{
  	\For{$c \in child_v$} {
	    $X_c$ = \Network($\tau_c$, $\frac{|\tau_c|\varepsilon}{|\tau_v|}$)\;
	 }
	 \Return \Sequence($\{X_c\}_{c \in child_v}$, $\frac{\varepsilon}{n_v|\tau_v|}$)	
	 
	}
			
	\If{$v$ is a $\operatorname{Parallel}$ node}{
	 	\For{$c \in child_v$} {
			$X_c$ = \Network($\tau_c$, $\min(\frac{|\tau_{c}|\varepsilon }{|\tau_v|}   ,\frac{1}{n_v (|\tau_v| n_v+1)})$)\;
		}
		\Return \Parallel($\{X_c\}_{c \in child_v}$)
	}	
  }
%----------------------------------5
  
\caption{Network($\tau$, $\varepsilon$) }\label{alg:approx}
  
\end{algorithm}

%We introduce the algorithm and prove that the algorithm indeed returns an $\varepsilon$-approximation of the completion time of the plan. 
%For a node $v$, let $\tau_v $ be the sub tree with $v$ as root and let $child_v$ be the set of children of $v$. 
%We use the notation $|\tau|$ to denote the total number of nodes in $\tau$.

%Note that given a task tree $\tau$, we can assume, without loss of generality, that there are no parallel nodes as children of parallel nodes, and no sequence nodes as children of sequence nodes, otherwise, we can easily merge them.
Algorithm~2, that implements the operator $\Network$, is a straightforward postorder traversal of the task tree. The only remaining issue is handling the error,  in an amortized approach, as seen in the proof of the following theorem.

\begin{theorem}
Given a task tree $\tau$, let $X_{\tau}$ be a r.v. representing the true
distribution of the completion time for the network. Then $\Network (\tau , \varepsilon)\approx_\varepsilon X_{\tau}$.
\end{theorem}

\begin{proof} 
By induction on $|\tau|$. 
Base: $|\tau|=1$, the node must be primitive, and $\Network$ will just return the distribution unchanged which is obviously an $\varepsilon$-approximation
of itself. Suppose the claim is true for $1 \leq |\tau| < n$. Let $\tau$ be a task tree of size $n$ and
let $v$ be the root of $\tau$. If $v$ is a $\operatorname{Sequence}$ node, by the induction hypothesis that $X_c\approx_{|\tau_c|\varepsilon/|\tau_v|} X_{\tau_c}$, and by Theorem~\ref{appSeqTheorem}, the maximum accumulated error is 
$\sum_{c \in child_v}
{|\tau_c|\varepsilon}/|\tau_v| + {\varepsilon}/|\tau_v|$ = $(n-1)\varepsilon/|\tau_v|+{\varepsilon}/|\tau_v|= \varepsilon$ for $v$, therefore, $\Sequence(\{X_c\}_{c\in child_v}, {\varepsilon}/{n}) \approx_\varepsilon X_\tau$ as required. 
If $v$ is a $\operatorname{Parallel}$ node, by the induction hypothesis that $X_c  \approx_{e_c} X_{\tau_c}$, where $e_c = \min(\frac{|\tau_{c}|\varepsilon}{|\tau_v|}, \frac{1}{n_v (|\tau_v| n_v+1)})$ 
so $\sum_{c\in child_v} e_c \leq   \sum_{c\in child_v} \frac{|\tau_{c}|\varepsilon}{|\tau_v|}  \leq \varepsilon  - \varepsilon/ |\tau_v|  $. Then, by Lemma~\ref{appPalTheorem}, using $K=|\tau_v|$ and $n=n_v$, we get that $\Parallel(\{X_c\}_{c\in child_v})\approx_\varepsilon  X_\tau$ as required. 
\end{proof}

\begin{theorem}\label{th:TTalgcomplexity}
Let  $N$ be the size of the  task tree  $\tau$,  and $M$ the size of the maximal support of each of the primitive tasks.
If $0 \leq \varepsilon \leq \frac{1}{N (N^2+1)}$ and $M < N/\varepsilon$, the  $\Network$ approximation algorithm 
runs in time $O((N^5 /\varepsilon^2)\log(N^3/\varepsilon^2))$, using
$O(N^3/\varepsilon^2) $ memory.
\end{theorem}

\begin{proof}
The run-time and space bounds can be derived from the bounds on $\Sequence$ and on $\Parallel$, as follows. In the $\Network$ algorithm, the trim accuracy parameter is less than or equal to ${\varepsilon}/{N}$. 
The support size (called $m$ in Theorem~\ref{appSeqComplexTheorem}) of the variables input to $\Sequence$ are $O(N^2/\varepsilon)$.
Therefore, the complexity of the $\Sequence$ algorithm is 
$O((N^4 /\varepsilon^2)\log(N^3/\varepsilon^2))$ and the complexity of the $\Parallel$ operator is $O((N^3/\varepsilon)\log(N))$. The time and space for sequence dominate, so the total time complexity is $N$ times the complexity of $\Sequence$ and the space complexity is that of $\Sequence$.
\end{proof}

If the constraining assumptions on $M$ and $\varepsilon$ in Theorem \ref{th:TTalgcomplexity} are lifted, the complexity is still polynomial:
replace one instance of $1/\varepsilon$ by $max(m,1/\varepsilon)$, and the other by $max(1/\varepsilon,N (N^2+1))$
in the runtime complexity expression.

\section{Complexity results}\label{sec:complexity}

We show that the deadline problem is NP-hard, 
even for a task tree consisting only of primitive tasks and one sequence node, i.e.
{\em linear plans}.
%We begin with a lemma on sums of random variables.

  %(we call the following problem SumDiscreteRV).
\begin{lemma} \label{SumDiscreteRV}
Let $Y=\{Y_1,\dots,Y_n\}$ be a set of discrete real-valued r.v.s specified by probability mass functions 
with finite supports, $T \in \mathbb{Z}$, and $p\in[0,1]$. Then, deciding  whether $Pr(\sum_{i=0}^{n} Y_{i}<T)>p $ is NP-Hard.
\end{lemma}

\begin{proof} By reduction from \textit{SubsetSum}~\cite[problem number SP13]{Garey:1990:CIG:574848}. Recall that
$\textit{SubsetSum}$ is: given a set $S=\{s_{1},\dots,s_{n}\}$ of integers, and integer target value $T$, is there a subset of $S$ whose sum is exactly $T$? 
Given an instance of \textit{SubsetSum}, create
the two-valued r.v.s $Y_{1},\dots,Y_{n} $ with $Pr(Y_i=s_i)=1/2$ and $Pr(Y_i=0)=1/2$. 
By construction, there exists a subset of $S$ summing to $T$ if and only if $Pr(\sum_{i=0}^{n} Y_{i}=T)>0$.

Suppose that algorithm $ A(Y,T,p) $ can decide $Pr(\sum_{i=0}^{n} Y_{i}<T)>p $ 
in polynomial time. Then, since the r.v.s $Y_i$ are two-valued uniform r.v.s, the only possible values of $p$ are integer multiples of ${1}/{2^n}$,
and we can to compute $p=Pr(\sum_{i=0}^{n} Y_{i}<T)$ using a binary search on $p$ using $n$ calls to $A$.
To determine whether $Pr(\sum_{i=0}^{n} Y_{i}=T)>0$, simply use this scheme twice, since $Pr(\sum_{i=0}^{n} Y_{i}=T)>0$
is true if and only if $Pr(\sum_{i=0}^{n} Y_{i}<T) < Pr(\sum_{i=0}^{n} Y_{i}<T+1)$.
\end{proof}
 
\begin{theorem}
Finding the probability that a task tree satisfies a deadline $T$ is NP-hard.
\end{theorem}

\begin{proof} Given a task tree consisting of  leaf nodes, all being children of a single sequence node, its makespan
is the sum of the completion times of the leaves. The theorem follows immediately from Lemma ~\ref{SumDiscreteRV}. 
\end{proof}


Finally, we consider the linear utility function, i.e. the problem
of computing an expected makespan of a task network.
Note that although for linear plans the {\em deadline problem} is NP-hard, the {\em expectation problem} is trivial because the expectation
of the sum of r.v.s $X_i$ is equal to the sum of the expectations of the $X_i$s.
For {\em parallel nodes}, it is easy to compute the CDF and therefore also easy to
compute the expected value.
Despite that, for task networks consisting of {\em  both} sequence nodes and
parallel nodes, these methods cannot be effectively combined, and in fact, we have:
\begin{theorem}
Computing the expected completion time of a task network is NP-hard.
\end{theorem}

\begin{proof}
By reduction from subset sum. Construct r.v.s (``primitive tasks'') $Y_i$ as in the proof of Lemma~\ref{SumDiscreteRV}, and denote
by $X$ the r.v. $\sum_{i=1}^{n} Y_{i}$. Construct one parallel node
with two children, one being the a sequence node having the
completion time distribution
defined by $X$, the other being a primitive task that has a completion time
$T_j$ with probability 1. (We will use more than one such case,
which differ only in the value of $T_j$, hence the subscript $j$).
Denote by $M_j$ the r.v. that represents the completion time
distribution of the parallel node, using this
construction, with the respective $T_j$. Now consider computing the
expectation of the $M_j$ for the following cases: $T_1 = T+{1}/{2}$
and $T_2 = T+{1}/{4}$.
Thus we have, for $j\in\{ 1, 2\}$, by construction and the definition of
expectation:
{\small
\begin{eqnarray*}
E[M_i] &=& T_j Pr(X\leq T_j) + \sum_{x > T_j} x~Pr(X = x)  \\
          &=& T_j Pr(X\leq T)+ \sum_{x \geq T+1} x~Pr(X = x)
\end{eqnarray*}}
where the second equality follows from the $Y_i$ all being
integer-valued r.v.s (and therefore $X$ is also integer valued).
Subtracting these expectations, we have $E[M_1]-E[M_2]=\frac{1}{4}P(X\leq T)$.
%\begin{eqnarray*}
%(T+\frac{1}{2}) Pr(X\leq T) + \sum_{x \geq T+1} x Pr(X=x) - \\
%(T+\frac{1}{4}) Pr(X\leq T) + \sum_{x \geq T+1} x Pr(X=x) 
%\end{eqnarray*}
%which reduces to $\frac{1}{4}P(X\leq T)$. Therefore, using the computed expected values, we can compute
Therefore, using the computed expected values, we can compute
$P(X\leq T)$, and thus also $P(X=T)$, in polynomial time.
%Using two such computations,  we can also compute $P(X=T)$,
%and the latter is non-zero if and only if there is a subset of $S$
%summing to $T$.
\end{proof}

\section{Empirical Evaluation}\label{sec:exp}


%\begin{table}[tbh]
%{\tiny
%\begin{tabular}{|l|l|l|| p{1.5cm}|p{1.5cm}|p{2.4cm}||l|l|l|l|l|}
%\hline
%\multicolumn{1}{|c|}{\multirow{2}{*}{Task Tree}} & \multicolumn{1}{c|}{\multirow{2}{*}{$N$}} & \multicolumn{1}{c|}{\multirow{2}{*}{$M$}} & \multicolumn{3}{c|}{Approx. alg. err., given $\varepsilon$} & \multicolumn{4}{c|}{Sampling alg. err., given \# samples} \\ \cline{4-10} 
%\multicolumn{1}{|c|}{} & \multicolumn{1}{c|}{} & \multicolumn{1}{c|}{} & \multicolumn{1}{c|}{0.1} & \multicolumn{1}{c|}{0.01} & \multicolumn{1}{c|}{0.001} & \multicolumn{1}{c|}{$10^{3}$} & \multicolumn{1}{c|}{$10^{4}$} & \multicolumn{1}{c|}{$10^{5}$} & \multicolumn{1}{c|}{$10^{6}$} \\ \hline \hline
%\multirow{3}{*}{Drive(DRC)} & 47 & 2 & [-0.0052, 0.0086] & [-0.0004, 0.0004] & [-$3.2 {\cdot} 10^{-5}$, $3.4 {\cdot} 10^{-5}$] & 0.0206 & 0.0072 & 0.0031 & 0.0009 \\ \cline{2-10} 
% & 47 & 4 & [-0.0096, 0.019] & [-0.0009, 0.0013] & [-$9.2 {\cdot} 10^{-5}$, $1.3 {\cdot} 10^{-4}$] & 0.0476 & 0.0075 & 0.0046 & 0.0011 \\\cline{2-10} 
% & 47 & 10 & [-0.014, 0.028] & [-0.0014, 0.0025] & [-$9.5 {\cdot} 10^{-5}$, $1.4 {\cdot} 10^{-4}$] & 0.0236 & 0.0083 & 0.0024 & 0.0015 \\ \hline
%\multirow{2}{*}{Seq} & 10 & 4 & [-0.027, 0.041] & [-0.0027, 0.0041] & [-$2.2 {\cdot} 10^{-4}$, $2.5 {\cdot} 10^{-4}$] & 0.0224 & 0.008 & 0.0017 & 0.0016  \\ \cline{2-10} 
% & 10 & 10 &[-0.0316, 0.0615]  & [0.0033, 0.0067] & [$-2.6 {\cdot} 10^{-4}$, $5.2 {\cdot} 10^{-4}$] & 0.027 & 0.0117 & 0.0038 & - \\ \hline 
%\multirow{1}{*}{Random} & 50 & 4 & 0.007 & 0.0007 &  - & 0.0243 & 0.0084 & 0.0024& - \\   \hline
%  \multirow{2}{*}{Logistics14} & 45 & 4  & 0.004 & 0.0004 & $3.3\cdot 10^{-5}$& 0.036 & 0.008 & 0.002 & 0.0006 \\ \cline{2-10}  
% & 45 & 10  &0.005 & 0.0005 & 0 & 0.03 & 0.013 & 0.002 &0.001  \\ \hline
%\end{tabular}
%\caption{Estimation errors}
%\label{tab:errors}
%}
%\vspace*{-0.5cm} 
%\end{table}

We examine our approximation bounds in practice, and compare the results to
exact computation of the CDF and to a simple stochastic sampling scheme. Three types of task trees
are used in this evaluation:
task trees used as execution plans for the (blinded for review) team entry in the DARPA
robotics challenge (DRC simulation phase),  linear plans (seq), and plans for
for the Logistics domain (from IPC2 http://ipc.icaps-conference.org/).
The primitive task distributions were uniform
distributions discretized to $M$ values. 

In the Logistics domain, packages are to be transported by trucks or airplanes.
Hierarchical plans were generated by the JSHOP2 planner \cite{nau2003shop2} for this domain and consisted of one parallel node (packages delivered in parallel),
with children all being sequential plans.
The duration distribution of all primitive tasks is uniform but the support parameters were determined by the type of the task, 
in some tasks the distribution is fixed (such as for load and unload) and in others the distribution depends on the velocity of the vehicle
and on the distance to be travelled.

\begin{table}[h]
{
\begin{tabular}{|l|l|l||l|l||l|l|}
\hline
\multirow{2}{*}{Task Tree} & \multirow{2}{*}{$N$} & \multirow{2}{*}{$M$} & \multicolumn{2}{l|}{Approx. alg.,  $\varepsilon$} & \multicolumn{2}{l|}{Sampling,\# samples} \\ \cline{4-7} 
 &  &  & 0.1 & 0.01 & $10^{4}$ & $10^{6}$ \\ \hline
 \hline
\multirow{3}{*}{Drive(DRC)} & 47 & 2 & \!\!\![-0.005, 0.009]\!\!\! & \!\!\![-0.0004, 0.0004]\!\!\! & 0.0072 & 0.0009 \\ \cline{2-7} 
 & 47 & 4 & \!\!\![-0.01, 0.02]\!\!\! & \!\!\![-0.0009, 0.001]\!\!\! & 0.0075 & 0.0011  \\ \cline{2-7} 
 & 47 & 10 & \!\!\![-0.01, 0.03]\!\!\! & \!\!\![-0.001, 0.003]\!\!\! & 0.0083 & 0.0015  \\ \hline
\multirow{2}{*}{Seq} & 10 & 4 & \!\!\![-0.03, 0.04]\!\!\! & \!\!\![-0.003, 0.004]\!\!\! & 0.008 & 0.0016 \\ \cline{2-7} 
 & 10 & 10 & \!\!\![-0.03, 0.06]\!\!\!  & \!\!\![0.003, 0.007]\!\!\! & 0.0117  & 0.001  \\ \hline % Added expected result as last number, to fix
\multirow{2}{*}{Logistics} & 45 &  4& \!\!\![-0.004,0.004]\!\!\! & \!\!\![-0.0004,0.0004]\!\!\! & 0.008 & 0.0006  \\ \cline{2-7} 
 & 45 & 10 & \!\!\![-0.005,0.006]\!\!\! & \!\!\![-0.0004,0.0006]\!\!\!& 0.013 & 0.001 \\ \hline
\end{tabular}
\vspace*{-0.2cm} 
\caption{Estimation errors}
\vspace*{-0.4cm} 
\label{tab:errors}
}
%\vspace*{-0.5cm} 
\end{table}

%\vspace*{-0.3cm} 
After running our approximation algorithm we also ran a variant
that uses an inverted version of the $\Trim$ operator, providing a {\em lower} bound of the CDF, as well as the upper bound
generated by Algorithm~\ref{alg:approx}. Running both variants allows us to bound the actual error, costing
only a doubling of the run-time. Despite the fact that our error bound is theoretically tight, in practice
and with actual distributions, according to Table~\ref{tab:errors}, the resulting error
in the algorithm was usually much better than the theoretical $\varepsilon$ bound.

%\begin{figure}[h!]
%\includegraphics[width=.49\textwidth]{plan31}
%\caption{Task network for pick-up challenge task}
%\label{fig:pickup}
%\vspace*{-0.35cm}
%\end{figure}


%
%\begin{table}[htb]
%{\tiny
%\begin{tabular}{|l|l|l|l||l|l|l||l|l|l|l|l|}
%\hline
%\multicolumn{1}{|c|}{\multirow{2}{*}{Task Tree}} & \multicolumn{1}{c|}{\multirow{2}{*}{$N$}} & \multicolumn{1}{c|}{\multirow{2}{*}{$M$}} & \multirow{2}{*}{Exact} & \multicolumn{3}{c|}{Approx. alg., with $\varepsilon$} & \multicolumn{4}{c|}{Sampling alg., with \# samples} \\ \cline{5-11} 
%\multicolumn{1}{|c|}{} & \multicolumn{1}{c|}{} & \multicolumn{1}{c|}{} &  & \multicolumn{1}{c|}{0.1} & \multicolumn{1}{c|}{0.01} & \multicolumn{1}{c|}{0.001} & \multicolumn{1}{c|}{$10^{3}$} & \multicolumn{1}{c|}{$10^{4}$} & \multicolumn{1}{c|}{$10^{5}$} & \multicolumn{1}{c|}{$10^{6}$}  \\ \hline \hline
%\multirow{2}{*}{Drive (DRC)} & 47 & 2 & 1.49 & 0.141 & 1.14 & 1.49 & 0.187 & 1.92 & 19.11 & 190.4  \\ \cline{2-11} 
%& 47 & 4 & 18.9 & 0.34 & 7.91 & 16.11 & 0.21 & 2.1 & 20.95 &  211.5 \\ \cline{2-11}
% & 47 & 10 & $> 2$h & 1.036 & 32.94 & 390.5 & 0.28 & 2.81 & 28.6 &279.1   \\ \hline
% \multirow{2}{*}{Seq 10} & 10 & 4 & 0.23 & 0.003 & 0.02 & 0.148  & 0.054 & 0.545 &5.336 & 54.22 \\ \cline{2-11} 
% & 10 & 10 & 10.22 & 0.008 & 0.073 & 0.692 & 0.071 & 0.724 & 7.18 & - \\ \hline
%\multirow{1}{*}{Random} & 50 & 4 & $>2$h & 1.1544 & 19.77 & 390.58 & 5.676 & 55.021 & 590.17& - \\   \hline
% \multirow{2}{*}{Logistics} & 45 & 4 & 373.3& 0.2 & 7 & 82.9 & 0.25 & 2.5 & 25.6 & 256 \\ \cline{2-11}  
% & 45 & 10 & $> 4$h & 2.19 & 120 & 6101 & 0.31 & 3.12 & 31.3 &314 \\ \hline
%\end{tabular}
%\caption{Runtime comparison (run times in seconds)}
%\label{tab:runtimes}
%}
%\end{table}
%
\begin{table}[h]
{
\begin{tabular}{|l|l|l|l||l|l||l|l|}
\hline
\multirow{2}{*}{Task Tree} & \multirow{2}{*}{$N$} & \multirow{2}{*}{$M$} & \multirow{2}{*}{Exact} & \multicolumn{2}{l|}{Approx., $\varepsilon$} & \multicolumn{2}{l|}{Sampling, \# samples} \\ \cline{5-8} 
 &  &  &  & 0.1 & 0.01 & $10^{4}$ & $10^{6}$ \\ \hline
 \hline
\multirow{3}{*}{Drive (DRC)} & 47 & 2 &1.49 & 0.141 & 1.14 &  1.92 &  190.4  \\ \cline{2-8} 
 & 47 & 4 &18.9 & 0.34 & 7.91  & 2.1 &  211.5  \\ \cline{2-8} 
  & 47 & 10 & $> 2$h & 1.036 & 32.94 & 2.81  &279.1 \\ \hline
\multirow{2}{*}{Seq} & 10 & 4 &  0.23 & 0.003 & 0.02  & 0.545 & 54.22  \\ \cline{2-8} 
 & 10 & 10 & 10.22 & 0.008 & 0.073  & 0.724 &  72.4 \\ \hline
\multirow{2}{*}{Logistics} & 45 & 4 & 373.3& 0.2 & 7 & 2.5 & 256  \\ \cline{2-8} 
 & 45 & 10 & $> 4$h & 2.19 & 120 &  3.12  &314  \\ \hline
\end{tabular}
\vspace*{-0.2cm} 
\caption{Runtime comparison (run times in seconds)}
\vspace*{-0.3cm} 
\label{tab:runtimes}
}
\end{table}
%\vspace*{-0.3cm} 

We ran the exact algorithm, our approximation algorithm with $\varepsilon \in \{ 0.1, 0.01, 0.001\}$, and a simple simulation with 
$10^3$ to $10^7$ samples, on networks from the DRC
implementation, sequence nodes with 10, 20, and 50 children, and 20 Logistics domain plans, and
several values of $M$. Results for a typical indicative subset (regretfully reduced due to space limits) 
are shown in tables \ref{tab:errors} (error comparison) and \ref{tab:runtimes} (runtime comparison).
%measured from the true result when available,
% and from the bounds generated by the approximation algorithm using $\varepsilon = 0.0001$
%when the exact algorithm timed out (over 2 hours). 
The exact algorithm times out in some cases.
Both our approximation algorithm and the sampling algorithm handle all these cases, as our algorithm's runtime is polynomial in $N$, $M$, and ${1}/{\varepsilon}$
as is the sampling algorithm's (time linear in number of samples). 

The advantage of the approximation algorithm is mainly in providing bounds with
certainty as opposed to the bounds in-probability provided by sampling. 
Additionally, as predicted by theory, accuracy of the approximation algorithm
improves linearly with ${1}/{\varepsilon}$ (and almost linear in runtime), whereas accuracy of sampling improves only as a square root of the number of
samples. Thus, even in cases where sampling initially outperformed the approximation algorithm, increasing the required accuracy for both algorithms,
eventually the approximation algorithm overtook the sampling algorithm.

\vspace*{-0.1cm} 

\section{Discussion}

Numerous issues remain unresolved, briefly discussed below.
Trivial improvements to the $\Trim$ operator are possible, such as the inverse version of the operator used to
generate a lower bound for the empirical results. Other candidate improvements are not performing trimming
(or even stopping a trimming operation) if the current support size is below $1/\varepsilon$, which may increase accuracy but also
the runtime. Another point is that in the combined algorithm, space and time complexity can be reduced by adding some $\Trim$ operations,
especially after processing a parallel node, which is not done in our version. This may reduce accuracy, 
a trade-off yet to be examined.
Another option is, when given a specific threshold, trying for higher accuracy in just the region of the
threshold, but how to do that is non-trivial. For {\em sampling} schemes such methods are known, including adaptive 
sampling \cite{bucher1988adaptive,lipton1990practical},
stratified sampling, and other schemes.
It may be possible to apply such schemes to deterministic algorithms as well - an interesting issue for future work.

Extension to continuous distributions: our algorithm can handle them by
pre-running a version of the $\Trim$ operator on the primitive task distribution. Since one cannot iterate over support values
in a continuous distribution, start with the smallest support value (even if it is $- \infty$), and find the value at which the CDF
increases by $\varepsilon$. This requires access to the inverse of the CDF, which is available, either exactly or approximately,
for many types of distributions. %In fact, this is precisely how Gaussian distributions were handled in our
%empirical evaluation.

We showed that the expectation problem is also NP-hard.
A natural question is on approximation algorithms for the expectation problem, but the answer here is not so obvious.
Sampling algorithms may run into trouble if the target distribution contains major outliers, i.e. values very far from
other values but with extremely low probability. Our approximation algorithm can also be used as-is to estimate
the CDF and then to approximate the expectation, but we do not expect it to perform well because our current $\Trim$
operator only limits the amount of probability mass moved at each location to $\varepsilon$, but does not limit the
``distance'' over which it is moved. The latter may be arbitrarily bad for estimating the expectation.
Possibly adding simple binning schemes to the $\Trim$ operator in addition to limiting the moved probability mass to $\varepsilon$
may work, another issue for future research. 

Related work on computing makespan distributions includes~\cite{hong2013computing}, which examines sum of  Bernoulli distributed r.v.s.
Other work examines both deterministic~\cite{mercier2007discrete} and  Monte-Carlo techniques~\cite{bucher1988adaptive,lipton1990practical}. 
Distribution of maximum of r.v.s was studied in~\cite{devroye1980generating}, with a focus mostly on continuous distributions.

Complexity of finding the probability that the makespan is under a given threshold in task networks was shown to
be NP-hard in \cite{hagstrom1988computational}, even when the completion time of each task has a Bernoulli distribution.
Nevertheless, our results are orthogonal as the source of the complexity in \cite{hagstrom1988computational} is in the graph structure,
whereas in our setting the complexity is due to the size of the support. In fact for linear plans (an NP-hard case in our setting),
the probability of meeting the deadline can be computed in low-order polynomial time for Bernoulli distributions, using straightforward dynamic programming.
Makespan distributions in series parallel networks in the i.i.d. case was examined in~\cite{gutjahr1992average}, without   
considering algorithmic issues. There is also a significant body of work on estimating the makespan of plans and schedules
\cite{herroelen2005project,fu2010towards,beck2007proactive}, within a context of a planner or schedueler. The analysis in these papers is
based on averaging or on limit theorems, and does not provide a guaranteed approximation scheme.

Copmuting the distribution of the makespan in trees is a seemingly trivial problem in probabilistic reasoning \cite{Pearl}.
Given the task network, it is straightforward to represent the distribution using a Bayes network (BN) that has one node
per task, and where the {\em children} of a node $v$ in the task networks are represented by BN nodes that are {\em parents} of
the BN node representing $v$. This results in a tree-shaped BN, where it is well known that probabilistic reasoning can be done in time linear
in the number of nodes, e.g. by belief propagation (message passing) \cite{Pearl,Kim}. The difficulty is in the potentially exponential size
of variable domains, which our algorithm, essentially a limited form
of approximate belief propagation, avoids by trimming.

Looking at makespan distribution computation as probabilistic reasoning leads to interesting issues for future research, such as 
how to handle task completion times that have dependencies, represented as a BN. 
Since reasoning in BNs is NP-hard even for binary-valued variables \cite{Dagum.aij,Cooper.ai}, this is unlikely in general.
But for cases where the BN toplogy is tractable, such as for BNs with bounded treewidth \cite{Bodlaender:2006:TCA:2092758.2092759},
or directed-path singly connected BNs \cite{ShimonyDomshlak.aiRN2003},
a deterministic polynomial-time approximation scheme for the makespan distribution may be achievable.
The research literature contains numerous {\em randomized} approximation schemes that handle depenencies \cite{Pearl,Yuan20061189}, especially for
the case with {\em no evidence}. In fact out original implementation of the sampling scheme in our DARPA robotics challenge entry handled dependent durations.
It is unclear whether such sampling schemes can be adapted to handle dependencies {\em and} arbitrary evidence, such as: ``the completion time of
compound task $X$ in the network is known to be exactly 1 hour from now''.
Finally, one might consider additional commonly used utility functions, such as a ``soft'' deadline: the utility is a constant $U$ before the deadline $T$, decreasing linearly to 0
until $T+G$ for some ``grace'' duration $G$, and 0 thereafter.

\clearpage
\bibliographystyle{plain}
\bibliography{deadline_arxiv}
%\input{deadline_arxiv.bbl}

\end{document}


deadline_arxiv.bbl
\begin{thebibliography}{10}

\bibitem{beck2007proactive}
J~Christopher Beck and Nic Wilson.
\newblock Proactive algorithms for job shop scheduling with probabilistic
  durations.
\newblock {\em J. Artif. Intell. Res.(JAIR)}, 28:183--232, 2007.

\bibitem{Bodlaender:2006:TCA:2092758.2092759}
Hans~L. Bodlaender.
\newblock Treewidth: Characterizations, applications, and computations.
\newblock In {\em Proceedings of the {32Nd} International Conference on
  Graph-{Theoretic} Concepts in Computer Science}, WG'06, pages 1--14, Berlin,
  Heidelberg, 2006. Springer-Verlag.

\bibitem{bonfietti2014disregarding}
Alessio Bonfietti, Michele Lombardi, and Michela Milano.
\newblock Disregarding duration uncertainty in partial order schedules? {Yes},
  we can!
\newblock In {\em Integration of AI and OR Techniques in Constraint
  Programming}, pages 210--225. Springer, 2014.

\bibitem{bucher1988adaptive}
Christian~G Bucher.
\newblock Adaptive sampling{:} an iterative fast {M}onte {C}arlo procedure.
\newblock {\em Structural Safety}, 5(2):119--126, 1988.

\bibitem{buyya2011sla}
Rajkumar Buyya, Saurabh~Kumar Garg, and Rodrigo~N Calheiros.
\newblock {SLA}-oriented resource provisioning for cloud computing:
  {Challenges}, architecture, and solutions.
\newblock In {\em Cloud and Service Computing (CSC), 2011 International
  Conference on}, pages 1--10. IEEE, 2011.

\bibitem{Cooper.ai}
Gregory~F. Cooper.
\newblock The computational complexity of probabilistic inference using
  {B}ayesian belief networks.
\newblock {\em Artificial Intelligence}, 42 (2-3):393--405, 1990.

\bibitem{Dagum.aij}
Paul Dagum and Michael Luby.
\newblock Approximating probabilistic inference in {B}ayesian belief networks
  is {NP}-hard.
\newblock {\em Artificial Intelligence}, 60 (1):141--153, 1993.

\bibitem{devroye1980generating}
Luc Devroye.
\newblock Generating the maximum of independent identically distributed random
  variables.
\newblock {\em Computers \& Mathematics with Applications}, 6(3):305--315,
  1980.

\bibitem{erol1994htn}
Kutluhan Erol, James Hendler, and Dana~S Nau.
\newblock {HTN} planning: Complexity and expressivity.
\newblock In {\em AAAI}, volume~94, pages 1123--1128, 1994.

\bibitem{fu2010towards}
Na~Fu, Pradeep Varakantham, and Hoong~Chuin Lau.
\newblock Towards finding robust execution strategies for {RCPSP/max} with
  durational uncertainty.
\newblock In {\em ICAPS}, pages 73--80, 2010.

\bibitem{gabaldon2002programming}
Alfredo Gabaldon.
\newblock Programming hierarchical task networks in the situation calculus.
\newblock In {\em AIPS’02 Workshop on On-line Planning and Scheduling}, 2002.

\bibitem{Garey:1990:CIG:574848}
Michael~R. Garey and David~S. Johnson.
\newblock {\em Computers and Intractability; {A} Guide to the Theory of
  {NP}-Completeness}.
\newblock W. H. Freeman \& Co., New York, NY, USA, 1990.

\bibitem{gutjahr1992average}
WJ~Gutjahr and G~Ch Pflug.
\newblock Average execution times of series--parallel networks.
\newblock {\em S{\'e}minaire Lotharingien de Combinatoire}, 29:9, 1992.

\bibitem{hagstrom1988computational}
Jane~N Hagstrom.
\newblock Computational complexity of {PERT} problems.
\newblock {\em Networks}, 18(2):139--147, 1988.

\bibitem{herroelen2005project}
Willy Herroelen and Roel Leus.
\newblock Project scheduling under uncertainty: {Survey} and research
  potentials.
\newblock {\em European journal of operational research}, 165(2):289--306,
  2005.

\bibitem{hong2013computing}
Yili Hong.
\newblock On computing the distribution function for the poisson binomial
  distribution.
\newblock {\em Computational Statistics \& Data Analysis}, 59:41--51, 2013.

\bibitem{kelly2008offline}
John~Paul Kelly, Adi Botea, and Sven Koenig.
\newblock Offline planning with {H}ierarchical {T}ask {N}etworks in video
  games.
\newblock In {\em AIIDE}, pages 60--65, 2008.

\bibitem{Kim}
Jin~H. Kim and Judea Pearl.
\newblock A computation model for causal and diagnostic reasoning in inference
  systems.
\newblock In {\em Proceedings of the 6th International Joint Conference on AI},
  1983.

\bibitem{lipton1990practical}
Richard~J Lipton, Jeffrey~F Naughton, and Donovan~A Schneider.
\newblock {\em Practical selectivity estimation through adaptive sampling},
  volume~19.
\newblock ACM, 1990.

\bibitem{mercier2007discrete}
Sophie Mercier.
\newblock Discrete random bounds for general random variables and applications
  to reliability.
\newblock {\em European journal of operational research}, 177(1):378--405,
  2007.

\bibitem{nau2003shop2}
Dana~S Nau, Tsz-Chiu Au, Okhtay Ilghami, Ugur Kuter, J~William Murdock, Dan Wu,
  and Fusun Yaman.
\newblock {SHOP2}: An {HTN} planning system.
\newblock {\em J. Artif. Intell. Res. (JAIR)}, 20:379--404, 2003.

\bibitem{nau1998control}
Dana~S Nau, Stephen~JJ Smith, Kutluhan Erol, et~al.
\newblock Control strategies in {HTN} planning: Theory versus practice.
\newblock In {\em AAAI/IAAI}, pages 1127--1133, 1998.

\bibitem{Pearl}
Judea Pearl.
\newblock {\em Probabilistic Reasoning in Intelligent Systems: {Networks} of
  Plausible Inference}.
\newblock Morgan Kaufmann, San Mateo, CA, 1988.

\bibitem{ShimonyDomshlak.aiRN2003}
S.~E. Shimony and C.~Domshlak.
\newblock Complexity of probabilistic reasoning in directed-path singly
  connected {B}ayes networks.
\newblock {\em Artificial Intelligence}, 151:213--225, 2003.

\bibitem{Yuan20061189}
Changhe Yuan and Marek~J. Druzdzel.
\newblock Importance sampling algorithms for {Bayesian} networks: {Principle}s
  and performance.
\newblock {\em Mathematical and Computer Modelling}, 43(9–10):1189 -- 1207,
  2006.

\end{thebibliography}
