# Mesh Refinement with Early Termination for Dynamic Feasibility Problems

Eduardo M. G. Vila  Eric C. Kerrigan  Paul Bruce  Department of Electrical and
Electronic Engineering, Imperial College London, SW7 2AZ, UK (e-mail:
<EMAIL_ADDRESS> Department of Electrical and Electronic Engineering,
and Department of Aeronautics, Imperial College London, SW7 2AZ, UK (e-mail:
<EMAIL_ADDRESS> Department of Aeronautics, Imperial College
London, SW7 2AZ, UK (e-mail<EMAIL_ADDRESS>
###### Abstract

We propose a novel early-terminating mesh refinement strategy using an
integrated residual method to solve dynamic feasibility problems. As a
generalization of direct collocation, the integrated residual method is used
to approximate an infinite-dimensional problem into a sequence of finite-
dimensional optimization subproblems. Each subproblem in the sequence is a
finer approximation of the previous. It is shown that these subproblems need
not be solved to a high precision; instead, an early termination procedure can
determine when mesh refinement should be performed. The new refinement
strategy, applied to an inverted pendulum swing-up problem, outperforms a
conventional refinement method by up to a factor of three in function
evaluations.

###### keywords:

constrained control, numerical methods for optimal control, large scale
optimization problems.

††thanks: This work was supported by the Science and Technology Facilities
Council (STFC) under a doctoral training grant (ST/V506722/1)

##  1 Introduction

Dynamical systems are modeled using differential (and algebraic) equations.
The dynamic variables are often constrained by lower and upper bounds, as well
as by initial and/or terminal conditions. Such a formulation constitutes a
dynamic feasibility problem (DFP) which, in general, may have zero, one, or
many feasible solutions. In the latter case, a merit function can be used to
find an _optimal_ solution, resulting in a dynamic optimization problem (DOP).
Many optimal control and system identification problems can be formulated as
DOPs. Moreover, a sequence of DFPs may be used to solve a DOP (Nie and
Kerrigan, 2023; Nita et al., 2022).

In contrast to ordinary optimization problems, where the decision variables
are finite-dimensional, DFPs are infinite-dimensional problems: dynamic
variables are functions of time (i.e., trajectories). In general, these
problems are intractable. Discretization methods are used to _transcribe_ the
DFP into an approximate, finite-dimensional optimization problem, which is
tractable by well-known algorithms.

Within (direct) transcription, the class of _shooting_ methods aims to find
control variables such that the propagated state variables satisfy the problem
constraints. They rely on the sensitivities of a time-marching integrator to
iteratively find a solution. Shooting methods often have highly nonlinear
relations between control variables and state constraints, resulting in poor
convergence rates (Kelly, 2017, Sect. 9.6). The class of _collocation_ methods
takes a different approach: both states and controls are discretized and
become decision variables of the optimization problem. The dynamic equations
are included as constraints, collocated at the discretization points. Direct
collocation methods mainly differ in the numerical integration scheme that is
used. Examples include trapezoidal and Hermite-Simpson collocation (Kelly,
2017, Sect. 3 and 4), as well as certain Runge-Kutta and Gauss-Legendre
families of collocation methods (Betts, 2020, Sect. 4.5).

Recently, the class of _integrated residual_ methods has enjoyed renewed
interest (Neuenhofen and Kerrigan, 2020; Nie and Kerrigan, 2023). Popular
among partial-differential equation solvers, integrated residual methods aim
to minimize the error on the integration of differential equations. This class
of methods can be seen as a generalization of direct collocation, sharing a
favorable convergence rate. Belonging to this class is the _least-squares_
method, which results in optimization problems of slightly increased
dimension; however, the method allows for problems with many algebraic
equations and singular arcs to be solved in cases where direct collocation
fails (Nie and Kerrigan, 2023).

In practice, it is more computationally efficient to solve a _sequence_ of
transcribed optimization problems, rather than a single one. Due to the
presence of discontinuities and regions of high gradients, it is favorable to
start with a coarse discretization (i.e., fewer mesh points), and
progressively refine the mesh towards finer discretizations. The mesh may be
refined by increasing the partitioning of the time domain ($h$-methods),
and/or by increasing the degree of discretization within each interval
($p$-methods).

State-of-the-art mesh refinement strategies solve the transcribed optimization
problems to a high precision. To improve each approximate solution, an
estimate of the local discretization error is used to refine the mesh by $h$
and/or $p$ methods. The general principle is to predict the mesh with the
least number of points necessary to achieve a desired accuracy. The mesh
prediction can use an integer programming technique so as to minimize the
maximum local error (Betts and Huffman, 1998). For problems with smooth
solutions, $p$-methods thrive due to their exponential convergence rate
(Trefethen, 2000). In the presence of non-smoothnesses, an $h$-method
increases convergence by partitioning the time domain into smooth sub-
intervals. Thus, a $p$-then-$h$ approach has shown significant computational
efficiency (Patterson et al., 2015).

Equally important to _how_ the mesh should be refined, is _when_ the
refinement should take place. For shooting methods, strategies have been
proposed where a refinement test is used to terminate each optimization
problem early (Polak, 1997, Sect. 3.3.3). Instead of attempting to predict the
ideal mesh prematurely, this approach gradually refines the mesh, as each
optimization problem progresses. Refinement tests may be of the _filtering_
type, where the early termination criterion is prescribed as a decreasing
tolerance. Alternatively, _convergence_ type tests terminate each optimization
problem once their rates of convergence stagnate, ensuring the convergence
rate of the infinite-dimensional problem is achieved (Polak, 1993).

Our contribution is a progressive mesh refinement strategy for the general
class of integrated residual methods. The proposed algorithm simultaneously
partitions the mesh ($h$-method), as well as adjusts the integration scheme,
ensuring sufficient integration accuracy. We use an early termination test
based on the estimated convergence rate which interrupts the optimization
problem, and triggers the mesh refinement scheme. The novel use of an early
termination test for the class of integrated residual methods results in
significant computational savings, when compared with a predictive refinement
strategy.

The remainder of this paper is organized as follows. Section 2 defines the
DFP, and describes the discretization process using an integrated residual
method. Section 3 outlines two mesh refinement strategies: progressive
refinement (using early termination), and predictive refinement (without early
termination). Section 4 applies both strategies to the inverted-pendulum
swing-up problem, where the performance is compared. Section 5 suggests future
improvements and extensions for the current implementation.

_Notation_ : Considering vectors in $\mathbb{R}^{n}$, $\|\cdot\|$ denotes the
Eucledian norm of a vector, and $\langle\cdot,\cdot\rangle$ denotes the inner
product of two vectors. Given a set of indices
$\mathcal{S}\subseteq\\{1,2,...,n\\}$, $\|\cdot\|_{\mathcal{S}}$ denotes the
Euclidean norm of a vector, only over the components in $\mathcal{S}$, and
$\langle\cdot,\cdot\rangle_{\mathcal{S}}$ denotes the inner product of two
vectors, only over the components in $\mathcal{S}$. $C^{0}(\Omega)$ denotes
the set of continuous functions defined over the domain $\Omega$. Given a
multivariate function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$, $\nabla
f(\cdot)$ denotes the gradient of $f(\cdot)$. The ceiling function
$\lceil\cdot\rceil$ yields the smallest integer greater or equal to the
argument.

##  2 Problem Statement

###  2.1 Dynamic Feasibility Problem

Within a continuous-time domain $\mathcal{T}:=[t_{0},t_{f}]\subset\mathbb{R}$,
we are interested in finding continuous differential dynamic variables
$\boldsymbol{x}:\mathcal{T}\rightarrow\mathbb{R}^{n_{x}}$ (e.g., states of a
system) and algebraic dynamic variables
$\boldsymbol{u}:\mathcal{T}\rightarrow\mathbb{R}^{n_{u}}$ (e.g., control
inputs) which satisfy a _residuals_ function
$r:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\rightarrow\mathbb{R}^{n_{r}}$
(e.g., system dynamics). The dynamic variables are bounded by sets of box
constraints $\mathcal{X}\subseteq\mathbb{R}^{n_{x}}$ and
$\mathcal{U}\subseteq\mathbb{R}^{n_{u}}$, respectively. Moreover,
$\boldsymbol{x}(\cdot)$ may be subject to initial conditions
$\mathcal{X}_{0}\subseteq\mathcal{X}$ as well as terminal conditions
$\mathcal{X}_{f}\subseteq\mathcal{X}$. We formulate the problem as

find $\displaystyle\boldsymbol{x}(\cdot),\,\boldsymbol{u}(\cdot),$ (DFP) s.t.
$\displaystyle
r(\dot{\boldsymbol{x}}(t),\boldsymbol{x}(t),\boldsymbol{u}(t))=0\quad$
$\displaystyle\widetilde{\forall}t\in\mathcal{T},$
$\displaystyle\boldsymbol{x}(t)\in\mathcal{X},\,\boldsymbol{u}(t)\in\mathcal{U}\quad$
$\displaystyle\forall t\in\mathcal{T},$
$\displaystyle\boldsymbol{x}(t_{0})\in\mathcal{X}_{0},\,\boldsymbol{x}(t_{f})\in\mathcal{X}_{f},$
$\displaystyle\boldsymbol{x}(\cdot)\in C^{0}(\mathcal{T}),$

where $\dot{\boldsymbol{x}}:\mathcal{T}\rightarrow\mathbb{R}^{n_{x}}$ is the
time-derivative of $\boldsymbol{x}(\cdot)$. Because $\boldsymbol{x}(\cdot)$
need only be piece-wise differentiable, the residuals function is enforced
‘almost everywhere’ in the Lebesgue sense
($\widetilde{\forall}t\in\mathcal{T}$). The algebraic variables
$\boldsymbol{u}(\cdot)$ are allowed to be discontinous.

We include the dynamic equations implicitly as residuals so as to allow for a
broader class of problems that may contain algebraic equations.

###  2.2 Integrated Residual Method

The transcription process begins by reformulating the constraint satisfaction
problem (DFP) into an optimization problem. Using a least-squares approach,
the normalized integrated squared residual is

$\boldsymbol{f}(\boldsymbol{x},\boldsymbol{u}):=\frac{1}{n_{r}(t_{f}-t_{0})}\int_{\mathcal{T}}\|r(\dot{\boldsymbol{x}}(t),\boldsymbol{x}(t),\boldsymbol{u}(t))\|^{2}\,\mathrm{d}t,$
(1)

which should be minimized. To simplify notation, we define the pair
$\boldsymbol{z}:=(\boldsymbol{x},\boldsymbol{u})$. The integrated residual
$\boldsymbol{f}$ is set as the objective functional of the infinite-
dimensional optimization problem

$\min_{\boldsymbol{z}\in\boldsymbol{\mathcal{Z}}}\quad\boldsymbol{f}(\boldsymbol{z}),$
(P)

where the constraint set is

$\boldsymbol{\mathcal{Z}}:=\\{(\boldsymbol{x},\boldsymbol{u})\,|\,\boldsymbol{x}(t)\in\mathcal{X},\,\boldsymbol{u}(t)\in\mathcal{U}\quad\forall
t\in\mathcal{T},\\\
\boldsymbol{x}(t_{0})\in\mathcal{X}_{0},\,\boldsymbol{x}(t_{f})\in\mathcal{X}_{f},\,\boldsymbol{x}(\cdot)\in
C^{0}(\mathcal{T})\\}.$ (2)

Because (P) is infinite-dimensional, the problem is, in general, intractable.
(P) needs to be discretized by suitable numerical methods so as to achieve an
approximate solution.

###  2.3 Problem Discretization

Firstly, the time domain $\mathcal{T}$ is partitioned by $n_{\mathcal{H}}-1$
points $\\{t_{1},t_{2},...,t_{n_{\mathcal{H}}-1}\\}$, while ensuring

$t_{0}<t_{1}<t_{2}<...<t_{n_{\mathcal{H}}-1}<t_{f}.$ (3)

By change of variables, we define normalized time intervals
$\tau_{i}\in[-1,1]$ such that

$t=\frac{t_{i}+t_{i-1}}{2}+\frac{t_{i}-t_{i-1}}{2}\tau_{i},$ (4)

where $t\in[t_{i-1},t_{i}]$, $\forall i\in\\{1,2,...,n_{\mathcal{H}}\\}$.
Additionally, let the normalized length of each interval be

$h_{i}:=\frac{t_{i+1}-t_{i}}{t_{f}-t_{0}}.$ (5)

The set of interval lengths is
$\mathcal{H}:=\\{h_{1},h_{2},...,h_{n_{\mathcal{H}}}\\}$.

Secondly, the dynamic variables $\boldsymbol{x}(\cdot)$ and
$\boldsymbol{u}(\cdot)$ are approximated by the piece-wise polynomial
functions $\tilde{\boldsymbol{x}}(\cdot,\cdot)$ and
$\tilde{\boldsymbol{u}}(\cdot,\cdot)$. Within each interval,
$\tilde{\boldsymbol{x}}(\cdot,\cdot)$ and
$\tilde{\boldsymbol{u}}(\cdot,\cdot)$ are polynomials in $\tau_{i}$ of degree
$n_{\mathcal{P}_{x}}$ and $n_{\mathcal{P}_{u}}$, respectively. Between
intervals, continuity of the differential variables is enforced by having
$\tilde{\boldsymbol{x}}_{i-1}(1,\cdot)=\tilde{\boldsymbol{x}}_{i}(-1,\cdot)$
$\forall i\in\\{2,3,...,n_{\mathcal{H}}\\}$, as illustrated in Figure 1.

Figure 1: Illustrations of piece-wise Lagrange polynomial approximations of
second and first degree, respectively.

Using Lagrange polynomials in the barycentric form, the approximate dynamic
variables $\tilde{\boldsymbol{x}}(\cdot,\cdot)$ and
$\tilde{\boldsymbol{u}}(\cdot,\cdot)$ are parameterised by the decision
variables $x$ and $u$, respectively. In the $i^{\text{th}}$ sub-interval, the
$j^{\text{th}}$ component of $\tilde{\boldsymbol{x}}(\cdot,\cdot)$ is given by

$\tilde{\boldsymbol{x}}_{j}^{i}(\tau,x)=\frac{\displaystyle\sum_{p=0}^{n_{\mathcal{P}_{x}}}\frac{w_{p}}{\tau-\tau_{p}}x_{p,j}^{i}}{\displaystyle\sum_{p=0}^{n_{\mathcal{P}_{x}}}\frac{w_{p}}{\tau-\tau_{p}}}$
(6)

where $\tau_{p}$ are the interpolation points, and $w_{p}$ are the barycentric
weights. By using Lagrange polynomials in the barycentric form, the
interpolation weights can be pre-calculated and reused, resulting in
computational savings. The interpolation points $\tau_{p}$ are chosen so as to
ensure numerical stability, e.g., Chebyshev or Legendre distributions (Berrut
and Trefethen, 2004). We denote the sets of interpolation points
$\mathcal{P}:=(\mathcal{P}_{x},\mathcal{P}_{u})$, where
$\mathcal{P}_{x}:=\\{\tau_{0},\tau_{1},...,\tau_{n_{\mathcal{P}_{x}}}\\}$ and
$\mathcal{P}_{u}:=\\{\tau_{0},\tau_{1},...,\tau_{n_{\mathcal{P}_{u}}}\\}$.

Thirdly, the mean integrated residual is approximated by a numerical
integration scheme. Any quadrature rule can be used, for the case of Gaussian
quadrature rule of degree $n_{\mathcal{Q}}$ with quadrature nodes
$\tau_{q}\in[-1,1]$, the cost functional (1) at each interval $i$ is
discretised as

$\frac{h_{i}}{2n_{r}}\sum_{q=1}^{n_{\mathcal{Q}}}w_{q}\bigg{\lVert}r\bigg{(}\frac{2\dot{\tilde{\boldsymbol{x}}}_{i}(\tau_{q},x)}{h_{i}(t_{f}-t_{0})},\tilde{\boldsymbol{x}}_{i}(\tau_{q},x),\tilde{\boldsymbol{u}}_{i}(\tau_{q},u)\bigg{)}\bigg{\rVert}^{2},$
(7)

where $\tau_{q}$ are the quadrature points, and $w_{q}$ are the quadrature
weights, which can also be pre-computed and reused. Quadrature rules with high
accuracy include Gauss-Legendre and Clenshaw-Curtis (Trefethen, 2008). Figure
2 illustrates the quadrature sampling of the integrand at quadrature points in
each sub-interval.

Figure 2: Illustration of Gaussian quadrature with $n_{\mathcal{Q}}=5$.

We denote the set of quadrature points
$\mathcal{Q}:=\\{\tau_{0},\tau_{1},...,\tau_{n_{\mathcal{Q}}}\\}$.

Due to change the of variables from $t$ to $\tau_{i}$,
$\dot{\tilde{\boldsymbol{x}}}$ is appropriately scaled in (7). We can use the
sum of (7) for each interval as the discretized objective function for the
transcribed optimization problems.

###  2.4 Sequence of Optimization Problems

In practice, instead of solving a single high-fidelity approximation of (P),
we solve a sequence of approximating problems of increasing accuracy. We
define the _family_ of finite-dimensional optimization problems, parameterized
by the mesh $\mathcal{M}:=(\mathcal{H},\mathcal{P},\mathcal{Q})$,

$\min_{z\in\mathcal{Z}_{\mathcal{M}}}\quad f_{\mathcal{M}}(z),$
($\textbf{P}_{\mathcal{M}}$)

where the pair $z:=(x,u)$, $f_{\mathcal{M}}=\sum_{i=1}^{n_{\mathcal{H}}}$ (7),
and the family of feasible sets

$\displaystyle\mathcal{Z}_{\mathcal{M}}:=\\{(x,u)\,|$ $\displaystyle
x^{i}_{p}\in\mathcal{X}\quad\forall p\in\\{0,1,...,n_{\mathcal{P}_{x}}\\},$
(8) $\displaystyle u^{i}_{p}\in\mathcal{U}\quad\forall
p\in\\{0,1,...,n_{\mathcal{P}_{u}}\\},$ $\displaystyle\forall
i\in\\{1,2,...,n_{\mathcal{H}}\\},$ $\displaystyle
x^{i-1}_{n_{\mathcal{P}_{x}}}=x^{i}_{0}\quad\forall
i\in\\{2,3,...,n_{\mathcal{H}}\\}$ $\displaystyle
x^{1}_{0}\in\mathcal{X}_{0},\,x^{n_{\mathcal{H}}}_{n_{\mathcal{P}_{x}}}\in\mathcal{X}_{f}\\}.$

In practice, we wish to _warm-start_ each problem with the solution of the
previous. Thus we require the constraint sets to be nested, i.e.,
$\mathcal{Z}_{\mathcal{M}_{m}}\subseteq\mathcal{Z}_{\mathcal{M}_{m+1}}\,\forall
m\in\\{1,2,...\\}$, where each $m$ is a mesh refinement iteration.

The structure of the constraints in ($\textbf{P}_{\mathcal{M}}$) can be
simplified (via elimination of variables for the continuity constraints) to
only upper and lower bounds on $z$. Box-constrained optimization algorithms
are tailored for this problem structure. Moreover, algorithms for least-
squares objectives may further exploit the problem structure.

A family of optimality functions $\theta_{\mathcal{M}}(\cdot)$ corresponding
to problems ($\textbf{P}_{\mathcal{M}}$) can be used as measures of
optimality. An optimality function is a nonpositive scalar-valued function
based on optimality conditions of an optimization problem. Minimizers are
zeros of a problem’s optimality functions (Polak, 1997, Sect. 1.1.1). For our
family of box-constrained problems ($\textbf{P}_{\mathcal{M}}$), first-order
quadratic approximations can be used as optimality functions, i.e.,

$\theta_{\mathcal{M}}(z):=\min_{z^{\prime}\in\mathcal{Z}_{\mathcal{M}}}\langle\nabla
f_{\mathcal{M}}(z),z^{\prime}-z\rangle+\frac{1}{2}\|z^{\prime}-z\|^{2},$ (9)

which can easily be evaluated if $\nabla f_{\mathcal{M}}(\cdot)$ is known.
This measure of optimality can be used to quantify the progress of an
iterative algorithm at each step.

##  3 Mesh Refinement

Two strategies to approximate a solution to P are described. Both approaches
contain an outer loop for mesh refinement, which yields a sequence of
transcribed problems
$\\{\textbf{P}_{\mathcal{M}_{0}},\textbf{P}_{\mathcal{M}_{1}},\textbf{P}_{\mathcal{M}_{2}},...\\}$.
Moreover, an inner loop iterates on the transcribed problems using an
optimization algorithm tailored for box constraints.

Firstly, we propose a _progressive_ refinement strategy that uses an early
termination test based on the convergence rate of the transcribed problem. A
quadrature error test is included to ensure accurate integration. Secondly, a
_predictive_ refinement strategy is also described, for comparison. Because no
mesh refinement scheme has yet been developed for integrated residual methods,
we synthesize an $h$-method inspired by the state-of-the-art (Betts, 2020,
Sect. 4.11), extended to include the quadrature error test.

###  3.1 Refinement Strategies

Starting with a coarse mesh $\mathcal{M}_{0}$ and an initial guess
$z_{0}\in\mathcal{Z}_{\mathcal{M}_{0}}$, Algorithm 1 progressively refines the
mesh, terminating each transcribed problem when slow convergence is detected
(as proposed in Polak (1997, Theorem 3.3.24), under convexity assumptions).

Algorithm 1 Refinement with Early Termination

$\mathcal{M}_{0}$, $z_{0}\in\mathcal{Z}_{\mathcal{M}_{0}}$, $\epsilon_{\rho}$,
$\epsilon_{f}$

$\mathcal{M}\leftarrow\mathcal{M}_{0}$

$k\leftarrow 0$

repeat

$z_{k+1}\leftarrow A_{\ref{alg:bci}}(\mathcal{M},z_{k})$ $\triangleright$ Use
Algorithm 5

$\rho_{k}\leftarrow\dfrac{f_{\mathcal{M}}(z_{k+1})}{f_{\mathcal{M}}(z_{k})}$
$\triangleright$ Estimated convergence rate

if $\rho_{k}\geq\epsilon_{\rho}$ then

$k\leftarrow k+1$

else

$(\mathcal{M},z_{k})\leftarrow A_{\ref{alg:progressive}}(\mathcal{M},z_{k})$
$\triangleright$ Use Algorithm 2

end if

until $f_{\mathcal{M}}(z_{k})\leq\epsilon_{f}$

The estimated convergence rate $\rho_{k}$ is the ratio of integrated residuals
between inner loop iterations. Once $\rho_{k}$ exceeds a convergence threshold
$\epsilon_{\rho}$, the optimization is terminated, and Algorithm 2
($A_{\ref{alg:progressive}}$) refines the mesh and interpolates a new $z$
approximation.

Algorithm 2 Progressive Mesh Refinement

$\mathcal{M}$, $z\in\mathcal{Z}_{\mathcal{M}}$, $\epsilon_{q}$, $\epsilon_{f}$

$(\mathcal{H},\mathcal{P},\mathcal{Q})\leftarrow\mathcal{M}$

$e_{q}\leftarrow\dfrac{|f_{\mathcal{M}}(z)-f_{{\mathcal{M}^{+}}}(z)|}{1+f_{{\mathcal{M}^{+}}}(z)}$
$\triangleright$ Estimate quadrature error

if $e_{q}>\epsilon_{q}$ then

increase $n_{\mathcal{Q}}$

else

for $i\in\\{1,2,...,n_{\mathcal{H}}\\}$ do

if $f^{i}_{\mathcal{M}}(z)/h_{i}\geq\epsilon_{f}$ then

Partition $h_{i}\in\mathcal{H}$ into two intervals

end if

end for

end if

$\mathcal{M}\leftarrow(\mathcal{H},\mathcal{P},\mathcal{Q})$

Interpolate $z$ into $\mathcal{Z}_{\mathcal{M}}$

return $(\mathcal{M},z)$

An estimate of the quadrature error can be obtained by using a nested
quadrature rule. To aid notation, $\mathcal{M}^{+}$ denotes a mesh with a
higher degree quadrature than $\mathcal{M}$. Should the local relative
quadrature error $e_{q}$ exceed a tolerance $\epsilon_{q}$, the quadrature
degree is increased (e.g., $n_{\mathcal{Q}}\leftarrow n_{\mathcal{Q}}+1$).

The second mesh refinement approach refines the mesh without early
termination. Betts (2020) models the differential equation error as being
proportional to $h^{p+1}$, where $p$ is the polynomial degree. The synthesized
$h$-method predicts the number of partitions as
$(e_{i}/\epsilon_{f})^{\frac{1}{p+1}}$, where
$e_{i}:=f^{i}_{\mathcal{M}}(z)/h_{i}$ and $\epsilon_{f}$ is the tolerance. The
optimization routine is terminated once the optimality function
$\theta_{\mathcal{M}}(\cdot)$ reaches an optimality threshold
$\epsilon_{\theta}$, as described in Algorithm 3.

Algorithm 3 Refinement without Early Termination

$\mathcal{M}_{0}$, $z_{0}\in\mathcal{Z}_{\mathcal{M}_{0}}$,
$\epsilon_{\theta}$, $\epsilon_{f}$

$\mathcal{M}\leftarrow\mathcal{M}_{0}$

$k\leftarrow 0$

repeat

$z_{k+1}\leftarrow A_{\ref{alg:bci}}(\mathcal{M},z_{k})$ $\triangleright$ Use
Algorithm 5

if $\theta_{\mathcal{M}}(z)\leq\epsilon_{\theta}$ then

$k\leftarrow k+1$

else

$(\mathcal{M},z_{k})\leftarrow A_{\ref{alg:predictive}}(\mathcal{M},z_{k})$
$\triangleright$ Use Algorithm 4

end if

until $f_{\mathcal{M}}(z_{k})\leq\epsilon_{f}$

Once each transcribed problem is solved, Algorithm 4
($A_{\ref{alg:predictive}}$) refines the mesh in a predictive manner. While
the quadrature branch is similar to Algorithm 2, the mesh intervals are
partitioned into a potentially larger number of intervals. The number of
partitions is estimated based on the approximation’s degree of polynomial $p$
(e.g., degree of $\mathcal{P}_{x}$), limited to a maximum of $p_{max}$
partitions.

Algorithm 4 Predictive Mesh Refinement

$\mathcal{M}$, $z\in\mathcal{Z}_{\mathcal{M}}$, $\epsilon_{q}$,
$\epsilon_{f}$, $p$, $p_{max}$

$(\mathcal{H},\mathcal{P},\mathcal{Q})\leftarrow\mathcal{M}$

$e_{q}\leftarrow\dfrac{|f_{\mathcal{M}}(z)-f_{{\mathcal{M}^{+}}}(z)|}{1+f_{{\mathcal{M}^{+}}}(z)}$
$\triangleright$ Estimate quadrature error

if $e_{q}>\epsilon_{q}$ then

increase $n_{\mathcal{Q}}$

else

for $i\in\\{1,2,...,n_{\mathcal{H}}\\}$ do

$e_{i}\leftarrow f^{i}_{\mathcal{M}}(z)/h_{i}$

if $e_{i}\geq\epsilon_{f}$ then

Partition $h_{i}\in\mathcal{H}$ into

$\min\\{\lceil(e_{i}/\epsilon_{f})^{\frac{1}{p+1}}\rceil,p_{max}\\}$ intervals

end if

end for

end if

$\mathcal{M}\leftarrow(\mathcal{H},\mathcal{P},\mathcal{Q})$

Interpolate $z$ into $\mathcal{Z}_{\mathcal{M}}$

return $(\mathcal{M},z)$

###  3.2 Optimization Algorithm

Many optimization methods are compatible with Algorithms 1 and 3. Because the
family of optimization problems $\textbf{P}_{\mathcal{M}}$ can be simplified
into only having box constraints, projected gradient methods (Bertsekas, 1982)
can be used. This class of methods is favorable because the constraints are
enforced at all iterations, can easily be warm-started, and are easy to
implement.

For a set of simple bounds $\mathcal{Z}=\\{z_{\ell}\leq z\leq z_{u}\\}$, the
projection operator $P_{\mathcal{Z}}(\cdot)$ is defined as

$P_{\mathcal{Z}}(z^{j}):=\begin{cases}z_{u}^{j}\quad&\textrm{if}\,z^{j}>z_{u}^{j},\\\
z_{\ell}^{j}&\textrm{if}\,z^{j}<z_{\ell}^{j},\\\
z^{j}&\textrm{otherwise},\end{cases}$ (10)

for each $j^{\text{th}}$ index of $z$. The set of binding indices
$\mathcal{B}$ corresponds to active (or almost active) bound constraints,
defined as

$\mathcal{B}(z):=\big{\\{}j\,|\,z^{j}_{\ell}\leq z^{j}\leq
z^{j}_{\ell}+\epsilon,\,\nabla f(z)^{j}>0\\\ \textrm{or}\quad
z^{j}_{u}-\epsilon\leq z^{j}\leq z^{j}_{u},\,\nabla f(z)^{j}<0\big{\\}},$ (11)

where $\epsilon=\textrm{min}(\epsilon_{z},\|z-P_{\mathcal{Z}}(z-\nabla
f(z))\|)$.

Algorithm 5 describes an iteration of a projected gradient method.

Algorithm 5 Box-constrained iteration

$\mathcal{M}$, $z_{0}$

$m\leftarrow 0$

repeat

$d_{m}$, $\alpha_{m}$ from Bertsekas (1982)

$z_{m+1}\leftarrow P_{\mathcal{Z}_{\mathcal{M}}}(z_{m}+\alpha_{m}d_{m})$

$m\leftarrow m+1$

until $\mathcal{B}(z_{m+1})=\mathcal{B}(z_{m})$

return $z_{m+1}$

Bertsekas (1982) describes how a descent direction $d$ can be obtained from
steepest-descent, conjugate-gradient, and Newton directions. An Armijo-like
line search method is also proposed to compute the step size $\alpha$.

##  4 Example

We assess the efficiency of the mesh refinement strategies by applying them to
a nonlinear problem, whose solution cannot be captured exactly by piece-wise
polynomials.

In this problem (Kelly, 2017, Sect. 6), a simple pendulum is attached to a
cart. The dynamics are modelled by the differential equations

$\ddot{q}_{1}=\frac{u_{1}+m_{2}\ell\sin(q_{2})\dot{q}_{2}^{2}+m_{2}g\cos(q_{2})\sin(q_{2})}{m_{1}+m_{2}\sin^{2}(q_{2})},\,\ddot{q}_{2}=\\\
\frac{u_{1}\cos(q_{2})+m_{2}\ell\cos(q_{2})\sin(q_{2})q_{2}^{2}+(m_{1}+m_{2})g\sin(q_{2})}{m_{1}+m_{2}\sin^{2}(q_{2})},$
(12)

where $q_{1}$ and $\dot{q}_{1}$ are the cart’s position (m) and velocity
(m/s), $q_{2}$ and $\dot{q}_{2}$ are the pendulum’s angle (rad) and angular
rate (rad/s), and $u_{1}$ is the control input (N). The masses for the cart
and pendulum are set to $m_{1}=1$ kg and $m_{2}=0.3$ kg, the pendulum’s length
$\ell=0.5$ m, and gravity’s acceleration $g=9.81$ m/s2. Starting from rest

$\displaystyle q_{1}(0)$ $\displaystyle=0\quad\quad q_{2}(0)=0$ (13)
$\displaystyle\dot{q}_{1}(0)$ $\displaystyle=0\quad\quad\dot{q}_{2}(0)=0,$

the cart must be actuated so as to swing up the pendulum in two seconds, i.e.,
satisfying

$\displaystyle q_{1}(2)$ $\displaystyle=1\quad\quad q_{2}(2)=\pi$ (14)
$\displaystyle\dot{q}_{1}(2)$ $\displaystyle=0\quad\quad\dot{q}_{2}(2)=0.$

The swing-up example is formulated as (DFP) with $q_{1}$, $\dot{q}_{1}$,
$q_{2}$, $\dot{q}_{2}$ as differential dynamic variables, and $u_{1}$ as an
algebraic dynamic variable. The differential equations (12) are included as
residuals in $r(\cdot,\cdot,\cdot)$, and boundary conditions (13), (14) are
enforced as initial and terminal constraints.

The DFP is discretised with piece-wise quadratic polynomials
($n_{\mathcal{P}_{x}}=3$, $n_{\mathcal{P}_{u}}=3$), and is solved using a
projected conjugate gradient calling an Armijo-like line-search (Schwartz and
Polak, 1997). We solve the DFP until an integrated residual tolerance
$\epsilon_{f}=1.15\times 10^{-1}$ is reached. We compare three approaches:

  1. 1.

No mesh refinement, using a sufficiently large number of intervals
$n_{\mathcal{H}}=20$, and quadrature $n_{\mathcal{Q}}=4$;

  2. 2.

The refinement scheme without early termination (Algorithm 3), starting with
$n_{\mathcal{H}}=2$, $n_{\mathcal{Q}}=4$, $\epsilon_{q}=1\times 10^{-1}$,
$\epsilon_{\theta}=2\times 10^{-4}$, and $p_{max}=4$;

  3. 3.

The refinement scheme with early termination (Algorithm 1) starting with
$n_{\mathcal{H}}=2$, $n_{\mathcal{Q}}=4$, $\epsilon_{q}=1\times 10^{-1}$, and
$\epsilon_{\rho}=0.999$.

To ensure a fair comparison, the values of $n_{\mathcal{H}}$,
$\epsilon_{\theta}$, and $\epsilon_{\rho}$ were experimentally chosen to yield
the least number of $\nabla r$ evaluations for each approach.

The three strategies were implemented in the Julia programming language,
solved using an Intel® Core™ i7-11370H laptop running Microsoft Windows
10.0.22621 with 16 GB of RAM. A summary of the results is presented in Table
1, where solve times are the minimum of 10 trials.

Table 1: Summary of performance results | Total $\nabla r$ evaluations | Solve time
---|---|---
No mesh refinement | 44 960 | 4.97 s
MR, no early termination | 40 900 | 4.77 s
MR with early termination | 10 088 | 1.23 s

Figure 3 shows the number of evaluations of the residuals Jacobian function
$\nabla r$ at each iteration of the optimization algorithm.

Figure 3: Performance comparison based on $\nabla r$ evaluations

Evaluations of $\nabla r$ dominate the total computational burden, with mesh
refinement and interpolation procedures becoming insignificant. The colored
area in each plot is therefore representative of the overall performance.

Figure 4 shows the mean integrated residual $f$ against the iteration number
on a log-log plot.

Figure 4: Convergence of the objective function $f$ for each strategy. Mesh
refinements appear as sharp kinks.

The convergence of different approaches can be examined.

The performance of the mesh refinement strategies can be described by a joint
analysis of Figures 3 and 4. Starting with a coarse mesh is advantageous not
only in the computation per iteration (smaller optimization problem), but also
allows for faster convergence on the first iterations.

Each optimization problem need not be solved to high precision, especially
initially. A phenomenon of diminishing returns is observed, where iterating a
small optimization problem no longer provides a significant decrease in
$f_{\mathcal{M}}$, compared to iterating on a larger problem, thus justifying
early termination.

##  5 Conclusions

We have described a novel mesh refinement strategy for solving DFPs using the
general class of integrated residual methods. On the example problem, the
progressive refinement was shown to outperform the conventional predictive
approach by up to a factor of three. Early termination of the optimization
problems was shown to improve the efficiency of the overall numerical method.

The following are suggestions for future research directions:

  * •

A more sophisticated mesh refinement procedure: an $h$-method is modest
compared to state-of-the-art ($hp$ or $ph$) methods. Care must be taken on how
the choice between $h$ and $p$ refinements is made, particularly considering
the sparsity structure of the transcribed optimization problem;

  * •

A selection of matured optimization algorithms: Both first-order and second-
order methods should be explored. Projected search methods may be paired with
recently proposed line-search algorithms that provide stronger, Wolfe-like
conditions (Ferry et al., 2021);

  * •

Optimized _hyper-parameters_ for Algorithms 1, 2, and 5: Zero-order methods
such as Bayesian optimization may be implemented to optimize the performance
and/or robustness of the algorithms;

  * •

A comparison between integrated residuals methods, and state-of-the-art direct
collocation and multiple shooting methods;

  * •

An extension of the refinement method to dynamic optimization problems,
possibly as per Nita et al. (2022).

## References

  * Berrut and Trefethen (2004) Berrut, J.P. and Trefethen, L.N. (2004).  Barycentric Lagrange Interpolation.  _SIAM Review_ , 46(3), 501–517.
  * Bertsekas (1982) Bertsekas, D.P. (1982).  Projected Newton Methods for Optimization Problems with Simple Constraints.  _SIAM Journal on Control and Optimization_ , 20(2), 221–246.  Publisher: Society for Industrial and Applied Mathematics.
  * Betts (2020) Betts, J.T. (2020).  _Practical methods for optimal control using nonlinear programming_.  Number 36 in Advances in design and control. Society for Industrial and Applied Mathematics, Philadelphia, third edition edition.
  * Betts and Huffman (1998) Betts, J.T. and Huffman, W.P. (1998).  Mesh refinement in direct transcription methods for optimal control.  _Optimal Control Applications and Methods_ , 19(1), 1–21.
  * Ferry et al. (2021) Ferry, M.W., Gill, P.E., Wong, E., and Zhang, M. (2021).  Projected-Search Methods for Bound-Constrained Optimization.  _arXiv:2110.08359 [math]_.
  * Kelly (2017) Kelly, M. (2017).  An Introduction to Trajectory Optimization: How to Do Your Own Direct Collocation.  _SIAM Review_ , 59(4), 849–904.
  * Neuenhofen and Kerrigan (2020) Neuenhofen, M.P. and Kerrigan, E.C. (2020).  An Integral Penalty-Barrier Direct Transcription Method for Optimal Control.  In _2020 59th IEEE Conference on Decision and Control (CDC)_ , 456–463.
  * Nie and Kerrigan (2023) Nie, Y. and Kerrigan, E.C. (2023).  Solving dynamic optimization problems to a specified accuracy: An alternating approach using integrated residuals.  _IEEE Transactions on Automatic Control_ , 68(1), 548–555.
  * Nita et al. (2022) Nita, L., Kerrigan, E.C., Vila, E.M.G., and Nie, Y. (2022).  Solving optimal control problems with non-smooth solutions using an integrated residual method and flexible mesh.  In _2022 IEEE 61st Conference on Decision and Control (CDC)_ , 1211–1216.
  * Patterson et al. (2015) Patterson, M.A., Hager, W.W., and Rao, A.V. (2015).  A ph mesh refinement method for optimal control.  _Optimal Control Applications and Methods_ , 36(4), 398–421.
  * Polak (1993) Polak, E. (1993).  On the use of consistent approximations in the solution of semi-infinite optimization and optimal control problems.  _Mathematical Programming_ , 62(1-3), 385–414.
  * Polak (1997) Polak, E. (1997).  _Optimization Algorithms and Consistent Approximations_ , volume 124 of _Applied Mathematical Sciences_.  Springer New York.
  * Schwartz and Polak (1997) Schwartz, A. and Polak, E. (1997).  Family of Projected Descent Methods for Optimization Problems with Simple Bounds.  _Journal of Optimization Theory and Applications_ , 92(1), 1–31.
  * Trefethen (2000) Trefethen, L.N. (2000).  _Spectral Methods in MATLAB_.  Software, Environments, and Tools. Society for Industrial and Applied Mathematics.
  * Trefethen (2008) Trefethen, L.N. (2008).  Is Gauss Quadrature Better than Clenshaw–Curtis?  _SIAM Review_ , 50(1), 67–87.


XXXXXX This is an open access article distributed under the terms of the
Creative Commons Attribution license
(http://creativecommons.org/licenses/by/3.0/), which permits unlimited use,
distribution and reproduction in any medium so long as the original work is
properly cited. XXXX

Patgiri et al.Empirical Study on Airline Delay Analysis and Prediction

# Empirical Study on Airline Delay Analysis and Prediction

Ripon Patgiri1  1  Sajid Hussain1  1  and Aditya Nongmeikapam1  1
11affiliationmark:  National Institute of Technology Silchar, Assam-788010,
India<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
(XXXX; XXXX)

###### Abstract

The Big Data analytics are a logical analysis of very large scale datasets.
The data analysis enhances an organization and improve the decision making
process. In this article, we present Airline Delay Analysis and Prediction to
analyze airline datasets with the combination of weather dataset. In this
research work, we consider various attributes to analyze flight delay, for
example, day-wise, airline-wise, cloud cover, temperature, etc. Moreover, we
present rigorous experiments on various machine learning model to predict
correctly the delay of a flight, namely, logistic regression with L2
regularization, Gaussian Naive Bayes, K-Nearest Neighbors, Decision Tree
classifier and Random forest model. The accuracy of the Random Forest model is
82% with a delay threshold of 15 minutes of flight delay. The analysis is
carried out using dataset from 1987 to 2008, the training is conducted with
dataset from 2000 to 2007 and validated prediction result using 2008 data.
Moreover, we have got recall 99% in the Random Forest model.

###### keywords:

Big Data, Big Data Analytics, Machine Learning, Airline Delay Analysis, Delay
Prediction, Random Forest

††articletype: Research Article EAI.EU

[1]Ripon Patgiri, Corresponding author. Email:

##  1 Introduction

MapReduce [8] is the most famous programming engine for developing parallel
and distributed application. It is the ideal programming language for large-
scale computation with low-cost commodity hardware. It is serving the data-
intensive computing since its inception. MapReduce consists of two functions,
namely, map and reduce; yet it is a very powerful programing engine. Moreover,
MapReduce is also strictly key/value pair programming model. Interestingly,
every problem in data-intensive computation can be converted into key/value
pair easily. MapReduce boomed in the IT industry and change the game of the
industries. It is applied in many fields in Computer Science and it is also
applied in interdisciplinary computing nowadays. For instance, the genome
project. MapReduce has left nothing untouched for large-scale computation. The
MapReduce programming model is easier to develop more customized framework
using map and reduce function. Thus, airline delay analysis is performed using
MapReduce, which takes few seconds to complete. Besides, other Big Data tools
can also be used, for instance, Apache Spark [18].

The airline data is growing day-by-day along with the number of new airlines,
due to the ever growing population of air travelers. There is always a quest
for airline data analytics to support more realistic data analysis with the
combination of weather data [5]. The MapReduce programming model is used to
experiment the behavior of flight delay. In addition, analysis is performed
using the dataset obtained from American Statistical Association [1] and
National Weather Service’s National Digital Forecast Database (NDFD) [13]. The
airline dataset and weather dataset from 1987 to 2008 is huge which cannot be
processed in conventional system. Therefore, MapReduce is the best programming
model to develop the airline delay analytics. Moreover, parallel processing
with low-cost commodity hardware is easier to experiment than HPC, and better
than a conventional systems.

The prediction is performed using various Data Mining techniques which are
used to uncover hidden truth [17]. We deploys conventional Machine learning
models, and the comparison of the results as shown in table 1, namely,
Logistic regression [7] with L2 regularization, K-Nearest Neighbor [6],
Gaussian Naive Bayes [11], decision tree classifier [14], and Random Forest
[10, 4]. The random forest model outperforms the other four models. However,
prediction is a very time-consuming process. As per our experience, more
training data makes the prediction more accurate. However, we reduce the
training data and just used the dataset from 2000 to 2007 to examine the
accuracy of the Machine Learning algorithms. The prediction is evaluated using
the data of 2008.

Table 1: Comparison of Machine Learning models Name | Precision | Recall | F1 score | Accuracy
---|---|---|---|---
Logistic Regression with L2 regularization | 0.83 | 0.67 | 0.74 | 0.62
K-Nearest Neighbors | 0.81 | 0.94 | 0.87 | 0.78
Gaussian Naive Bayes | 0.82 | 0.93 | 0.87 | 0.78
Decision Tree Classifier | 0.81 | 0.92 | 0.87 | 0.77
Random Forest | 0.82 | 0.99 | 0.90 | 0.82

###  1.1 Contribution

The contribution of this work is summarized as follows-

  * •

Development of Airline Analytics, which analyze the large dataset for flight
delays.

  * •

Merging of the weather dataset with the airline dataset to get more realistic
results.

  * •

Analysis report generation helps in the decision making process.

  * •

Prediction of flight delay using various machine learning models.

###  1.2 Organization

The paper is organized as follows- section 2 discusses on the background and
related work on flight delay analysis and prediction. Section 3 exposes the
data sources and various attributes of the dataset in analysis and prediction.
Section 4 provides detailed analysis of airline delay with weather data.
Section 5 predicts the delay of flights using various available machine
learning models. Finally, section 6 concludes the work.

##  2 Background and Related work

The delay in airlines is very frequent, and it occurs due to weather
conditions. The airline delay analysis is an important analysis to know the
behavior of a flight and its future probable delay. There are numerous
opportunities to work for the airline industry. The realistic airline data are
available for free to experiment and develop interesting solutions to existing
dilemmas. However, similar work has been done in the article [15, 9, 5, 3, 2]
for delay analysis. The article [12] analyze the latency of air traffic. There
are various parameters to consider in analyzing the air traffic data. More
precisely, the main motivation is to analyze the airline delay using the
combination of weather data and air traffic data. This creates large sized air
traffic data to analyze. Rebollo and Balakrishnan [15] provides 100 most-
delayed origin-destination pairs from operation dataset. Bloem et al. [3]
works on various cost functions of the airline and analyses the delay.
Ferguson et al. [9] estimate the cost of a delay of flight. Choi et al. [5]
merges weather data and airline data. Choi et al. [5] achieves 80.36 %
accuracy of the Random Forest model using 2005 to 2015 airline and weather
data. Belcastro et al. [2] designs a predictor based on MapReduce. Belcastro
et al. [2] achieves 74.2% accuracy with the delay threshold of 15 minutes,
however, the Random Forest model outperforms this accuracy.

###  2.1 Logistic Regression

The Logistic Regression technique is developed by D. R. Cox [7] in 1958. After
its invention, logistic regression become an integral part of Machine
Learning. The purest form of Logistic Regression by itself is very much
susceptible to the problem of overfitting. To get rid of this problem,
Logistic Regression is often used with L2 regularization. This leads to more
precise and accurate predictions, with the overfitting taken care of.

###  2.2 Gaussian Naive Bayes

It’s a Machine Learning Technique that uses a probabilistic approach for
classification of records. It is one of the two types of Gaussian Classifiers
[11]. It assumes Conditional Independence of the different attributes of the
data used for the classification. The second type of Gaussian Classifier is
called the BBN (Bayesian Belief Network), it performs the classification task
without the help of the assumption of Conditional Independence. In the Naive
Bayes Classifier, we find out the probability of occurrence of each of the
possible classes given the values of the particular attributes. Then, the
class for which this probability would turn out to be the largest will be
declared as the target class for the given attribute values.

For example, let us consider a dataset with two possible target classes,
namely ’yes’ and ’no’ and two attributes ’A’ and ’B’. And, the attribute
values given are $A=a$ and $B=b$. So, the probabilities to be calculated are :

$P(yes/A=a,B=b)$

and

$P(no/A=a,B=b)$

Now,

$P(yes/A=a,B=b)$ $=P(A=a,B=b/yes)P(yes)/P(A=a,B=b)$
$=P(A=a/yes)P(B=b/yes)P(yes)/k$

Here, we have considered the marginal probability

$P(A=a,B=b)$

to be a constant k, since it will be same for all the classes.

###  2.3 K-nearest neighbors

The K-Nearest Neighbors (KNN) is a classifier which clusters K nearest items
[6]. This KNN technique is extremely useful in learning and prediction. The
total nearest neighbors are taken as 10 which provides us the best results. It
locates the nearest neighbors in the data point space, and then, labels the
unknown data point as the same class label, as that of the known class label.

###  2.4 Decision Tree

The decision tree is one of the most basic Machine Learning Techniques [14].
The key feature of the decision tree is the impurities in the dataset. These
impurities suggest how heterogeneous or homogeneous a given set of data is. If
the impurities turn out to be zero, it means that the data set is classified.
These impurities are found by calculating entropy, gini index,
misclassification error, etc for the given data.

Figure 1: Relation among Entropy, Gini Index and Misclassification error.

Entropy-

$E(t)=-\sum_{i}(p(i/t)\log_{2}(p(i/t)))$

Gini Index-

$Gini(t)=1-\sum_{i}(p(i/t))^{2}$

Misclassification error-

$error(t)=1-\max_{i}(p(i/t))$

Where, $p(i/t)$ is the fraction of records belonging to class $i$ at node $t$
and $i$ varies from $0$ to $c-1$. $c$ is the total number of possible classes.

###  2.5 Random Forest

The Random Forest creates a huge set of random trees which is very time
consuming, but more accurate than other models in predicting [10, 4]. The
Random Forest model is a set of decision tree created within a tree for better
results. The features are selected randomly to form a sub-tree. Therefore,
there are many random trees in training periods. These trees are bagged
repeatedly for $T$ times. The Random Forest probability output is

$p(c/v)=\frac{1}{T}\sum_{t}^{T}p_{t}(c/v)$

The total number of trees taken in the experiment is 50.

###  2.6 Confusion Matrix

Table 2: Confusion Matrix |  | Predicted |
---|---|---|---
|  | NO | YES |
Actual | NO | $\sum TN$ | $\sum FP$ |
YES | $\sum FN$ | $\sum TP$ |

Table 2 shows the confusion matrix. The confusion matrix is composed of the
number of actual delayed and predicted delayed flights [16]. The True Negative
(TN) is actually not delayed and predicted to be not delayed. The False
Positive (FP) is actually not delayed and predicted to be delayed. The False
Negative (FN) is actually delayed, but predicted to be not delayed. The True
Positive (TP) is actually delayed and predicted to be delayed. The summation
$\sum$ represents the total number of such events. Based on the Table 2,
precision, recall, F1 score and accuracy are calculated as follows-

$Precision=\frac{\sum TP}{\sum FP+\sum TP}$ $Recall=\frac{\sum TP}{\sum
FN+\sum TP}$ $F1score=2\times\frac{Precision\times Recall}{Precision+Recall}$
$Accuracy=\frac{\sum TN+\sum TP}{\sum TN+\sum FP+\sum FN+\sum TP}$

##  3 Data Source

The data source is the American Statistical Association [1] and National
Weather Service’s National Digital Forecast Database (NDFD) [13]. The data
source contains US flight and weather data. We have taken data from 1987 to
2008 for analysis purpose and from 2000 to 2007 to train the models. Later,
the resultant prediction is validated using 2008 dataset. There are various
attributes in the joined dataset which are enlisted below-

  * •

Year: 1987 to 2008

  * •

Month: 1 (January) to 12 (December)

  * •

Day of month: 1 to 31

  * •

Day of week: 1 (Monday) to 7 (Sunday)

  * •

Unique carrier: Code of the airliner, e.g., AI (Air India)

  * •

Departure delay: The departure delay is measured in minutes.

  * •

Origin: The origin city of flight to the destination.

  * •

Destination: The destination city to reach by flight from origin.

  * •

Distance: The distance is measured in miles.

  * •

Maximum and minimum temperature: The degree of temperature is represented in
Fahrenheit.

  * •

Maximum and minimum visibility: The visibility is represented in miles.

  * •

Mean wind speed: The mean wind speed is represented in m/s.

  * •

Precipitation: It is represented in inches.

  * •

Cloud cover: The cloud cover is fullness of cloud in the sky.

##  4 Delay Analysis

This work uses two different types of datasets, one data set gives the details
regarding the delays in the airline industry and the other dataset gives the
weather information for the different cities in the United States. This work
is a combination of analysis of the delays encountered by the different
airlines in services and the prediction of this delay for the future flights.
The datasets are used for the analysis purpose. The predictor also uses the
same datasets in filtered form. The first dataset is obtained from the
American Statistical Association [1] and it has the attributes are, namely,
year, day of week, day of month, unique carrier name, distance traveled,
scheduled and actual departure time, etc. The second dataset is obtained from
the National Weather Service’s National Digital Forecast Database (NDFD) [13].
It contains information regarding the maximum and minimum temperature,
visibility, wind speed, precipitation, and cloud cover. These two datasets are
joined on the basis of the date using a MapReduce code and served to the
prediction algorithms. The joined data are used for analysis purpose. Now,
let’s discuss the two parts in details next.

Table 3: Experimental environment setup Name | Configuration
---|---
CPU | i7-4770 CPU @ 3.40GHz
RAM | 4GB
HDD | 500 GB each
OS | Ubuntu 16.10 64-bit
Hadoop | V2.7.3
Nodes | 20

The experimental environment is set up using the above configuration as
depicted in Table 3. Hadoop cluster is set up to conduct the experiment. The
total number of nodes is 20 with 4GB RAM of each. This cluster is sufficient
to join the dataset of airline and weather.

Figure 2: Airline-wise delay analysis without weather data. X-axis represents
airlines and Y-axis represents number of delays. Figure 3: Day-wise delay
analysis without weather data. X-axis represents days and Y-axis represents
number of delays. Figure 4: Delay analysis on maximum temperature in
Fahrenheit. X-axis represents maximum temperatures and Y-axis represents
number of delays. Figure 5: Delay analysis on minimum temperature in
Fahrenheit. X-axis represents minimum temperatures and Y-axis represents
number of delays. Figure 6: Delay analysis on maximum visibility. X-axis
represents maximum visibility and Y-axis represents number of delays. Figure
7: Delay analysis on minimum visibility. X-axis represents minimum visibility
and Y-axis represents number of delays. Figure 8: Delay analysis on mean wind
speed. X-axis represents mean wind speed and Y-axis represents number of
delays. Figure 9: Delay analysis on precipitation. X-axis represents
precipitation and Y-axis represents number of delays. Figure 10: Delay
analysis on cloud cover. X-axis represents cloud cover and Y-axis represents
number of delays.

###  4.1 Joining the dataset

MapReduce is used for this task. Our experiment is carried out by joining the
above two described data sets to obtain a unified data set, which is used for
the detailed analysis and prediction purpose. For this, two mappers and a
single reducer are used. The first mapper takes input the airline delay data
as input and the second mapper takes the weather data as input. The outputs of
these mappers are then joined on the basis of the date, i.e. the key, and a
new csv file is generated. This file is the unified one having all the
required fields or attributes. It is this CSV file that is used for the whole
experiment.

###  4.2 Detail analysis

For a manager of an airliner, there are an abundance of queries to be
performed regarding the delay. For instance, in which day, the probability of
a flight delay is the highest? The result can be easily analyzed by the
manager of the airline through the Figure 2. To get this result, we perform
logical analysis on airline data and weather data generated as csv. The
obtained csv file is used to analyze the delays in the airlines on the basis
of different attributes to obtain the results shown in graphical form 4, 5, 6,
7, 8, 9, 10. We spawn many mappers and reducers to analyze the delay. Here,
the mappers are given as input, the previously obtained CSV file in the
KeyValueInputFormat. The mappers find out the delay for a particular key/value
input pair, and output in the key/value format. The intermediate key is the
value of the attribute under consideration and the value is the delay. The
reducer fetches these key/value pairs in shuffled form and sums up all the
delay values for a particular value of the attribute under consideration. The
count of the total number of flights and the number of delayed flights is also
maintained. Finally, the reducer outputs the attribute value, the total number
of flights, the number of delayed flights as the key, the average delay and
the probability of delay as the value. The Record Writer writes the final
result in the mentioned output file in HDFS.

The analysis results are plotted in graphical form as shown in Figure 2, 3, 4,
5, 6, 7, 8, 9, 10. Figures 2, 3 do not consider the weather data. The figure 2
reveals that the airline DL has the highest probability of delay while average
delay airline CO is the highest. Figure 3 shows that most of the flight delays
on Friday. However, other Figures 4, 5, 6, 7, 8, 9, 10 combined airline data
and weather data to analyze various scenarios. Figure 4, 5, 6, 7, 8, 9, 10 are
reported the delay analysis in the case of maximum temperature, minimum
temperature, maximum visibility, minimum visibility, wind speed,
precipitation, and cloud cover respectively. Thus, this analytics enhances the
decision making process of the airline managers. Figure 4 and 5 reports the
delay on maximum and minimum temperature in Fahrenheit. The maximum delay
occurs in below 60 Fahrenheit in Figure 4 and below 40 Fahrenheit in Figure 5.
The figure 6 and 7 report the maximum delay in below 29.67 and below 2
respectively. Figure 8 shows that the average delay become high when the mean
wind speed high. Similarly, more cloud cover and precipitation cause delay as
shown in Figure 9 and 10 respectively.

##  5 Prediction

Figure 11: Prediction analysis and results. X-axis represents precision,
recall, F1 score, and accuracy and Y-axis represents score in an scale of 100.
Figure 12: Comparison among TN, FP, FN, and TN with machine learning model.
X-axis represents TN, FP, FN and TN, and Y-axis represents counts on
prediction.

For this part Python’s NumPy and SciPy Library were used. Python’s excellent
Scikit-learn machine learning package is used to build the predictive models
(Logistic regression, Random Forest, KNN, Gaussian Naive Bayes and Decision
Tree Classifier) and compare their performance. First the previously mentioned
CSV file (combination of weather and delay data) is loaded as feature matrix
by filtering out all the non-numeric attributes. The matrix consists of all
the data from the airline delay joined with the weather data. This feature
matrix is of $(1566226)\times(8)$ size. The training dataset contained all the
data from 2001 to 2007 and the test dataset contained the 2008 data. In the
first Iteration Logistic Regression with L2 regularization was used. The delay
is classified into binary classes:

  * •

Delayed (more then 15 mins delay)

  * •

Not Delayed (less then 15 mins delay).

After that, the Random forest model (50 trees) is used, it gives a better
accuracy, and increase in true positives which signifies improvement in our
model as shown in the table 1.

###  5.1 Resultant confusion matrix

|  | Predicted |
---|---|---|---
|  | NO | YES |
Actual | NO | $17981$ | $23129$ |
YES | $57915$ | $116583$ |
(a) Confusion Matrix of logistic regression with L2 regularization

|  | Predicted |
---|---|---|---
|  | NO | YES |
Actual | NO | $3437$ | $37673$ |
YES | $9810$ | $164688$ |
(b) Confusion Matrix of K-NN

|  | Predicted |
---|---|---|---
|  | NO | YES |
Actual | NO | $5191$ | $35919$ |
YES | $11394$ | $163104$ |
(c) Confusion Matrix of Gaussian Naive Bayes

|  | Predicted |
---|---|---|---
|  | NO | YES |
Actual | NO | $4176$ | $36934$ |
YES | $13173$ | $161325$ |
(d) Confusion Matrix of Decision tree

|  | Predicted |
---|---|---|---
|  | NO | YES |
Actual | NO | $457$ | $38653$ |
YES | $532$ | $175966$ |
(e) Confusion Matrix of Random Forest

Figure 13: Confusion matrix

The table 13a, 13b 13c, 13d and 13e show the confusion matrix after
experimentation using logistic regression with L2 regularization, K-NN,
Gaussian Naive Bayes, Decision tree and Random Forest model. The precision,
recall, F1 score, and accuracy of logistic regression with L2 regularization
are 0.83, 0.67, 0.74, and 0.62 respectively. The K-NN gives precision=0.81,
recall =0.94, F1 score = 0.87, accuracy = 0.78. The Gaussian Naive Bayes gives
precision = 0.82, recall = 0.93, F1 = 0.87, and accuracy = 0.78. The Decision
tree gives precision = 0.81, recall = 0.92, F1 = 0.87, and accuracy = 0.77.
The random forest achieves the highest and these are precision = 0.82, recall
= 0.99, F1 = 0.90, and accuracy = 0.82.

Similarly, to these two models the other Machine Learning Techniques described
above are also used for the same prediction purpose, and their results and
performances are mentioned using tables and graphs. All these Machine Learning
Techniques were used directly from the previously mentioned Python libraries.
The total number of True Negative, False Positive, False Negative and True
Positive is depicted in the figure 12.

###  5.2 Prediction results and analysis

We train using air traffic dataset from 2000 to 2007 and analyze the
prediction using 2008 dataset. In this case, we omit all data from 1987 to
2000 to see how accurately can we predict airline delay with a very small set
of data. We have split the training and testing dataset into 80:20 ratio.
Surprisingly, the Random Forest model performs outstanding with an accuracy of
81% using a very small set of data. The K-Nearest Neighbor and Gaussian Naive
Bayes perform with an accuracy of 78% while decision tree gives 77%. We
demonstrate the prediction results in figure 11. As depicted in the figure 11,
we compare the performance of the Random Forest model, Logistic Regression
With L2 Regularization, K-Nearest neighbor, Gaussian Naive Bayes model and
Decision Tree model. In this comparison, the Random Forest model outperforms
the other four models.

##  6 Conclusion

In this airline delay analysis and prediction, we have found that an airline
can be scheduled in a better way if we uses analytics. As we have shown in the
logical analysis reports, airline manager can perform analysis on various
queries on delays. Moreover, the flight delay is predicted by training a very
small set of data from 2000 to 2007 which is validated using 2008 dataset.
Thus, a flight can be scheduled, organized and analyzed in a much better way.
The Random Forest outperforms other four machine learning models, however,
K-NN and Gaussian Naive Bayes perform similar. The accuracy of the Random
Forest model is 82% with a threshold of 15 minutes. However, the work can also
be carried out by Apache Mahout for machine learning library with Spark which
is our future work. Finally, we believe that this work creates understanding
in experimentation using weather and airline data for the development of new
Airline Delay Analytics.

## References

  * [1] A. S. Association.  Airline on-time performance.  Retrieved on 20/03/2017 from http://stat-computing.org/dataexpo/2009/, 2009.
  * [2] L. Belcastro, F. Marozzo, D. Talia, and P. Trunfio.  Using scalable data mining for predicting flight delays.  ACM Trans. Intell. Syst. Technol., 8(1):5:1–5:20, 2016.
  * [3] M. Bloem, H. Huang, and N. Bambos.  Approximating the likelihood of historical airline actions to evaluate airline delay cost functions.  In 2012 IEEE 51st IEEE Conference on Decision and Control (CDC), pages 508–513, Maui, HI, USA, 2012. IEEE.
  * [4] L. Breiman.  Random forests.  Machine Learning, 45(1):5–32, 2001.
  * [5] S. Choi, Y. J. Kim, S. Briceno, and D. Mavris.  Prediction of weather-induced airline delays based on machine learning algorithms.  In 2016 IEEE/AIAA 35th Digital Avionics Systems Conference (DASC), pages 1–6, Sacramento, CA, USA, 2016. IEEE.
  * [6] T. Cover and P. Hart.  Nearest neighbor pattern classification.  IEEE Trans. Inf. Theor., 13(1):21–27, 2006.
  * [7] D. R. Cox.  The regression analysis of binary sequences.  Journal of the Royal Statistical Society. Series B (Methodological), 20(2):215–242, 1958.
  * [8] J. Dean and S. Ghemawat.  Mapreduce: Simplified data processing on large clusters.  In Proceedings of the 6th Conference on Symposium on Opearting Systems Design & Implementation - Volume 6, pages 10–10, Berkeley, CA, USA, 2004. USENIX Association.
  * [9] J. Ferguson, A. Q. Kara, K. Hoffman, and L. Sherry.  Estimating domestic us airline cost of delay based on european model.  Transportation Research Part C: Emerging Technologies, 33:311 – 323, 2013.
  * [10] T. K. Ho.  Random decision forests.  In Proceedings of the Third International Conference on Document Analysis and Recognition (Volume 1) - Volume 1, ICDAR ’95, pages 278–, Washington, DC, USA, 1995. IEEE Computer Society.
  * [11] G. H. John and P. Langley.  Estimating continuous distributions in bayesian classifiers.  In Proceedings of the Eleventh Conference on Uncertainty in Artificial Intelligence, pages 338–345, San Francisco, CA, USA, 1995. Morgan Kaufmann Publishers Inc.
  * [12] Z. Liu and J. Heer.  The effects of interactive latency on exploratory visual analysis.  IEEE Transactions on Visualization and Computer Graphics, 20(12):2122–2131, 2014.
  * [13] N. W. S. (NWS).  National digital forecast database (ndfd).  Retrieved on 20/03/2017 from http://www.nws.noaa.gov/ndfd/index.htm, 2017\.
  * [14] J. R. Quinlan.  Induction of decision trees.  Mach. Learn., 1(1):81–106, 1986.
  * [15] J. J. Rebollo and H. Balakrishnan.  Characterization and prediction of air traffic delays.  Transportation Research Part C: Emerging Technologies, 44:231 – 241, 2014.
  * [16] M. Sharma, S. Sharma, and G. Singh.  Performance analysis of statistical and supervised learning techniques in stock data mining.  Data, 3(4), 2018.
  * [17] M. Sharma, G. Singh, and R. Singh.  Stark assessment of lifestyle based human disorders using data mining based learning techniques.  IRBM, 38(6):305 – 324, 2017.
  * [18] M. Zaharia, M. Chowdhury, M. J. Franklin, S. Shenker, and I. Stoica.  Spark: Cluster computing with working sets.  In Proceedings of the 2Nd USENIX Conference on Hot Topics in Cloud Computing, pages 10–10, Berkeley, CA, USA, 2010. USENIX Association.