\chapter{Introduction}

Many of the recent advances in modern sciences have been possible due to the improvements made in computational and optimisation theory developed over the course of the last decades. From electrical engineering and mechanics, to economics and molecular modelling, all these sciences have taken advantage of the latest developments in the area of mathematical optimisation. The latter is the process by which the optimal solution to a problem, also known as optimum, is produced.
\\

Since very ancient times, people have been developing models and theories to deal with various optimisations problems in order to achieve the best possible results for specific tasks in everyday life. While interesting in theory, those ideas had very little practical use due to the daunting amount of computational effort required. It was not until the advent of the computer that those early thoughts have been resurrected and resulted in what is now regarded as a growing branch of applied mathematics.
\\

On the other hand, whilst computers have evolved rapidly over the last years and are able to deal with an enormous amount of computation, the problem of precision (i.e. accurate calculation) is still one of great importance even today. The issue arises from the floating-point representation widely available in today's computers which is known to have limited precision. However, certain scientific applications cannot accept such rounding errors as they may well lead to catastrophic events. For example, in computational geometry we may need to deal with infinite amount of precision in order to accurately track the position of objects and points in space. This is known as an emerging trend in exact computation, explained further in this work.

\section{Motivation}

In this context of approximation and computability, differential equations, introduced in the 17th century by Newton and Leibniz, play a central role in modern mathematics and have countless applications in almost all branches of contemporary science. Several numerical approaches for computing the solutions to differential equations have been developed, including the well-known Euler and Runge-Kutta methods. They have been shown, however, to suffer from a great loss of precision as their error estimation is too conservative to be of any practical use \cite{ref1}.
\\

One classical problem is the famous \textit{initial-value problem}, which states that, given an initial condition and an ordinary/partial differential equation that models some evolution of a system, we can determine the unknown function describing the underlying process. For such problems, a novel technique for computing the unique solution up to any desired accuracy has been proposed in \cite{ref1} (such a solution is guaranteed to exist under certain assumptions). 
\\

The basic idea is that, at each stage of the computation, one can obtain two continuous maps which provide upper and lower bounds for the solution, essentially giving the precise error. Furthermore, the two maps have the additional property of being continuous and piecewise linear, making them suitable to be determined by means of a linear programming algorithm, as presented in \cite{ref1}.
\\

A related problem to the one described above formulated again in \cite{ref1}, \cite{ref2}, concerns the idea of \textit{consistency} of two given maps subject to some constraints and will make most of the subject of the present work. A pair of functions $(f, g)$, representing function and derivative information respectively, is called \textit{consistent} if one can find a third map $h$ which is approximated by the first component and whose derivative information is approximated by the second component of the pair. Assuming that $f$ is defined by lower and upper constraints on its value and $g$ is constrained to lie within rectangles (or hyper-rectangles, if referring to higher dimensional spaces) then a \textit{consistency witness} $h$, which is piecewise linear, is guaranteed to exist. In addition, one can determine the least and greatest piecewise linear maps consistent with the information from $f$ and $g$ via a linear programming algorithm as shown in both \cite{ref1} and \cite{ref2}.

\section{Objectives}

This discussion brings us to the aim of this project, which is twofold: 

\begin{itemize}
	\item Firstly, to implement the above optimisation problem in the form of a linear programming algorithm when the derivative constraints are specified as hyper-rectangles and to determine the global bounding surfaces;
	\item Secondly, to study a more general setting in which the derivative constraints are considered to lie within convex polyhedra, rather than hyper-rectangles, which is a more challenging problem.
\end{itemize}

The final outcome of this project should be in the form of a linear programming framework that implements the aforementioned algorithms and allows a user to determine the consistent piecewise linear map given some input functions $f$ and $g$ (the function and derivative information, respectively).

\chapter{Background}

This section provides the necessary background for understanding the theoretical aspects of this project. We begin by introducing the notion of exact computation \cite{ref3}, a new emerging paradigm in the field of modern computation. We also discuss the significance of this work in the context of exact geometric computation \cite{ref4}. Finally, we introduce the main concepts of optimisation and linear programming that the reader should be familiar with.

\section{Exact Computation}

The underlying nature of computation is par excellence numerical: numbers have been at the heart of all calculations since very ancient times, while modern mathematics developed a whole range of theories to formalise many aspects of computable numbers. Early computers had the sole purpose of performing large complex calculations (the so-called number crunchers), which then turned into the original mass-produced computers, in the form of friendly pocket calculators. Although computers have evolved significantly and moved towards more abstracted and higher-level programs, numerical computation remains a major pillar of modern computer technology.
\\

In particular, scientific computation is a rapidly growing inter-disciplinary field that makes use of advanced numerical capabilities. It is considered as adding a new dimension to the classical methods of theory and experimentation, sometimes referred to as the ``third scientific method'' of computation \cite{ref3}.

\subsection{Fixed-Point Paradigm}

At its core, scientific computation is subject to the \textit{fixed-precision paradigm} of computation. Under this representation, numbers are expressed using a fixed number of digits after the decimal/radix point, thus providing fixed computational precision (usually machine-dependant). 
\\

Based on this specification, one can use several approaches to limit the unavoidable rounding errors. For instance, a \textit{mild form} of fixed-precision can be applied such that computations are performed up to a user-specified precision level. Although we can set a very high level of precision that we may think is satisfactory, the build-up of \textit{round-off errors} that accumulate may produce totally unexpected results. We can illustrate this with an example from \cite{ref5}. Consider the sequence $a_n$ defined recursively as:

\[
a_n = \left\{\begin{array}{cll}
\displaystyle\frac {11}2    & , & n = 0 \\\\
\displaystyle\frac {61}{11} & , & n = 1 \\\\
111 - \displaystyle\frac{1130 - \frac {3000}{a_{n-1}}}{a_n} & , & n \ge 2
\end{array}\right. .
\]
\\

Using the Unix utility \lit*{bc}, we can compute the terms of the sequence $a_n$ up to some fixed precision of $k$ decimal places, which we shall denote by $a_n^{(k)}$. Performing calculations with 5 decimal places, gives (note that the results have been rounded for presentation purposes, as in \cite{ref5}):
\\

\[\begin{array}{c}
\begin{tabular}{|c|r|}
\hline
$a_0^{(5)}$ & $5.500$  \\ \hline
$a_1^{(5)}$ & $5.500$  \\ \hline
$a_2^{(5)}$ & $5.500$  \\ \hline
$a_3^{(5)}$ & $5.500$  \\ \hline
$a_4^{(5)}$ & $ 5.648$ \\ \hline
$a_5^{(5)}$ & $ 5.242$ \\ \hline
\end{tabular}

\hspace{2cm}

\begin{tabular}{|c|r|}
\hline
$a_6^{(5)}$    & $ -3.241$ \\ \hline
$a_7^{(5)}$    & $  283.1$ \\ \hline
$a_8^{(5)}$    & $103.738$ \\ \hline
$a_9^{(5)}$    & $100.209$ \\ \hline
$a_{10}^{(5)}$ & $100.012$ \\ \hline
$a_{11}^{(5)}$ & $100.001$ \\ \hline
\end{tabular}
\end{array}\]
\\

One would therefore believe that the sequence converges to $100$. However, computing the number $a_{100}$ with higher and higher precisions will contradict our expectations (the ``exponents'' below indicate repeating digits, e.g. $1.2^34=1.2224$):

\[\begin{array}{c}
\begin{tabular}{|c|l|}
\hline
$a_{100}^{(5)}$   & $100.0^41$           \\ \hline
$a_{100}^{(30)}$  & $100.0^{29}1$        \\ \hline
$a_{100}^{(60)}$  & $100.0^{57}997$      \\ \hline
$a_{100}^{(100)}$ & $100.0^{17}98\ldots$ \\ \hline
\end{tabular}

\hspace{2cm}

\begin{tabular}{|c|l|}
\hline
$a_{100}^{(110)}$ & $100.0^{7}92$        \\ \hline
$a_{100}^{(120)}$ & $-3.790\ldots$       \\ \hline
$a_{100}^{(130)}$ & $5.9^{7}8697\ldots$  \\ \hline
$a_{100}^{(140)}$ & $5.9^{7}87925\ldots$ \\ \hline
\end{tabular}
\end{array}\]
\\

Here we can notice that if the precision is either $5$, $30$, $60$, $100$ or even $110$ decimal places, then our expectation from above holds ($a_n\to 100$ as $n$ grows larger). However, when increasing the precision even further, we obtain rather spurious results. 

The actual limit of the sequence is equal to $6$, since we can evaluate the general term $a_n$ as:

\[a_n = \frac {6^{n+1} + 5^{n+1}}{6^n + 5^n}.\]
\\

This example therefore shows how even the \textit{mild form} of fixed-precision can lead to flawed results. It is therefore unclear what level of precision needs to be set in advance to a program such as \lit*{bc} in order to obtain the correct answer. As we could see, up until 110 decimal places we got roughly the same (wrong) result and we actually had to consider precision above 130 decimal places to arrive at the right answer.

\subsection{Floating-Point Representation}

Similar to the fixed precision paradigm, modern computers use floating-point arithmetic to perform real number calculations. This representation can be seen as a trade-off between range and precision when approximating real number: with a given precision, the floating-point model is able to represent both numbers of small magnitude with many bits of significance or conversely, large magnitude and few bits of significance. While many different representations have been developed in the past, the one defined by IEEE 754 has emerged as the industry standard. This standard is a step forward towards addressing the issue of portability and therefore makes errors in floating point computation \textit{machine-independent}.
\\

However, this representation is essentially just as inaccurate as the fixed-point model, since we must approximate real numbers by their nearest representable one. Thus, the same rounding errors will occur in practice, e.g. when small inaccuracies propagate in successive iterations like the one that we have just seen when computing a simple mathematical limit.

\subsection{Towards an Alternative to the f.p. Paradigm} 

Even with industry-leading standards, rather intractable problems arise from the presence of round-off errors and compromise the \textit{robustness} of the f.p. paradigm. We can, at the arithmetic level, increase precision via techniques such as double extended precision, guard-bits, gradual underflow etc, which can usually be implemented in hardware. Interval arithmetic is another well-known method. Looking from a geometric perspective, an idea would be to divide the input and computed data into combinatorial and numerical, and to give precedence to the former when making decisions. An argument in favour of this approach is that we can allow the numerical data to be perturbed in order to maintain the combinatorial data. This avoids ``topological inconsistencies'' and can be implemented for simple cases, but the intractable nature of combinatorial problems makes it rather difficult to deal with more general cases \cite{ref3}.
\\

Therefore, in the light of these non-robustness issues, one needs a completely different approach to handle cases in which the correct and exact answer is required. A new direction in the literature of computation is called, unsurprisingly, the \textit{exact computation paradigm}. According, to \cite{ref3} or \cite{ref4}, this paradigm assumes a computational process that:

\begin{enumerate}
	\item represents all the underlying mathematical objects \textit{exactly};
	\item all branching decisions are error-free.
\end{enumerate}

As a result, multi-precision arithmetic is a necessary condition (but not sufficient) for exact computation. A different issue is that exact computation will naturally come at the expense of performance. It therefore makes sense to target only those applications which are not \textit{cycle-critical} (i.e. we afford to incur some sort of computational slow-down). For example, exact computation cannot be avoided in computational number theory and in many aspects of algebra (e.g. testing the irreducibility of a polynomial).
\\

Finally, we conclude this section with the idea of \textit{weak exact computation}, similar in intent with the mild form of the f.p. paradigm. As explained in the numerical example above, one starts the computations using some fixed bound $k$ on the precision. However, one would need to perform a thorough analysis as to what minimal values of $k$ will indeed give the exact results. This is clearly a non-trivial problem and several suggestions have been proposed: in the same \cite{ref3}, the theory of root bounds has been developed; in \cite{ref5}, the exact real arithmetic offers lower and upper bounds guarantees that are trustworthy.

\section{Linear Programming}

A great deal of problems encountered in the real world involve maximisation or minimisation of certain quantities. Most often, these quantities we seek to optimise are profits (in the case of maximisation) as well as costs (in the case of minimisation). In linear programming we therefore aim to minimise/maximise a certain linear function (which we usually call objective function) subject to some linear constraints that describe the restrictions of our problem. 
\\

One would immediately think that calculus, developed by Leibniz and Newton in the 17th century, can deal with this types of problems very elegantly given the arsenal of techniques readily available. However, calculus is inadequate since it can only imply that the maxima and/or minima of some objective function lie on the boundaries of the sets determined by the constraints. Thus, a special set of techniques and algorithms is needed to deal with such linear programming problems (this is what the term ``programming'' really refers to in this case).
\\

As such, linear programming plays a very important role in the fields of mathematics and business, finding many applications in management science and operations research.

\subsection{Canonical Forms for LP Problems}

As already discussed, there are two ways to define a linear programming problem: either as a minimisation or maximisation problem. The canonical forms for each of those is presented below along with some standard terminology \cite{ref6}. 
\\

\begin{definition} A problem is said to be a \textit{canonical maximisation linear programming problem} if it has the following form:

\[
\begin{array}{ll}
\text{Maximise}   & f(x_1, x_2, \ldots, x_n) = c_1x_1 + c_2x_2 + \ldots + c_nx_n - d \\
\text{subject to} & a_{11}x_1 + a_{12}x_2 + \ldots + a_{1n}x_n \le b_1 \\
	              & a_{21}x_1 + a_{22}x_2 + \ldots + a_{2n}x_n \le b_2 \\
	              & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \vdots \\
	              & a_{m1}x_1 + a_{m2}x_2 + \ldots + a_{mn}x_n \le b_m \\
	              & x_1, x_2, \ldots x_n \ge 0.
\end{array}
\]

\end{definition}

\begin{definition} A problem is said to be a \textit{canonical minimisation linear programming problem} if it has the following form:

\[
\begin{array}{ll}
\text{Minimise}   & g(x_1, x_2, \ldots, x_n) = c_1x_1 + c_2x_2 + \ldots + c_nx_n - d \\
\text{subject to} & a_{11}x_1 + a_{12}x_2 + \ldots + a_{1n}x_n \ge b_1 \\
	              & a_{21}x_1 + a_{22}x_2 + \ldots + a_{2n}x_n \ge b_2 \\
	              & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \vdots \\
	              & a_{m1}x_1 + a_{m2}x_2 + \ldots + a_{mn}x_n \ge b_m \\
	              & x_1, x_2, \ldots x_n \ge 0.
\end{array}
\]

\end{definition}

In both canonical forms, the first $m$ constraints are said to be \textit{main constraints}, while the second $n$ constraints are called \textit{non-negativity constraints}.

\begin{definition} The linear function $f$ and $g$ in Definition 1.1 and 1.2 are said to be \textit{objective functions}.
\end{definition}

\begin{definition} The set of all points $(x_1, x_2, \ldots, x_n)$ satisfying the $m + n$ constraints of a canonical maximisation/minimisation linear programming problem is said to be the \textit{constraint} set of the problem. Any element of the constraint set is said to be a \textit{feasible point} or \textit{feasible solution}.
\end{definition}

\begin{definition} Any feasible solution of a canonical maximisation (respectively minimisation) linear programming problem which maximises (respectively minimises) the objective function is said to be an \textit{optimal solution}.
\end{definition}

\subsection{Polyhedral Convex Sets}

Let us state some additional results that are essential towards assessing the optimality of feasible solutions in linear programming (see also \cite{ref6}).

\begin{changemargin}{0cm}{0cm}

\begin{definition} Let $\mathbf{x} = (x_1, x_2, \ldots, x_n)$, $\mathbf{y} = (y_1, y_2, \ldots, y_n)\in\mathbb{R}^n$. Then the \textit{line segment} between $\mathbf{x}$ and $\mathbf{y}$ (including the endpoints) is given by:

\[\alpha\mathbf{x} + (1-\alpha)\mathbf{y},\ 0\le \alpha\le 1.\] 

\end{definition}

\begin{definition} Let $S\subset \mathbb{R}^n$. The set $S$ is said to be \textit{convex} if for all $\mathbf{x} = (x_1, x_2, \ldots, x_n)$, $\mathbf{y} = (y_1, y_2, \ldots, y_n)\in S$ it holds that:

\[\alpha\mathbf{x} + (1-\alpha)\mathbf{y},\ 0\le \alpha\le 1.\]

\end{definition}

\begin{definition} The set of points $(x_1, x_2, \ldots, x_n)\in\mathbb{R}^n$ satisfying an equation of the form

\[a_1x_1 + a_2x_2 + \ldots + a_nx_n = b\]

is said to be a \textit{hyperplane} of $\mathbb{R}^n$. The set of points $(x_1, x_2, \ldots, x_n)\in\mathbb{R}^n$ satisfying an inequality of the form

\[a_1x_1 + a_2x_2 + \ldots + a_nx_n \lessgtr b\]

is said to be a \textit{closed half-space} of $\mathbb{R}^n$.
\end{definition}

\begin{theorem} The constraint set of a canonical maximisation/minimisation linear programming problem is convex. Such a set is said to be a polyhedral convex set.	
\end{theorem} 

\begin{definition} Let $\mathbf{x} = (x_1, x_2, \ldots, x_n)\in\mathbb{R}^n$. The norm of $\mathbf{x}$, denoted $\Vert \mathbf{x} \Vert$, is given by:

\[\Vert \mathbf{x} \Vert = \sqrt{x_1^2 + x_2^2 + \ldots + x_n^2}.\]

\end{definition}

\begin{definition} Let $r\ge 0$. The set of points $\mathbf{x} = (x_1, x_2, \ldots, x_n)\in\mathbb{R}^n$ such that

\[\Vert \mathbf{x} \Vert \le r\]

is said to be the \textit{closed ball of radius} $r$ \textit{centred at the origin}.
\end{definition}

\begin{definition} A set $S\subset\mathbb{R}^n$ is said to be \textit{bounded} if there exists $r\ge 0$ such that every element of $S$ is contained in the closed ball or radius $r$ centred at the origin. A set $S\subset\mathbb{R}^n$ is said to be \textit{unbounded} if it is not bounded.
\end{definition}

\begin{definition} Let $S$ be a convex set in $\mathbb{R}^n$. A point $e\in S$ is said to be an extreme point of $S$ if there do not exist $\mathbf{x}, \mathbf{y}\in S$ and $\alpha\in (0, 1)$ such that

\[\textbf{e} = \alpha \textbf{x} + (1-\alpha) \textbf{y}.\]

\end{definition}

\begin{theorem}\label{T1} If the constraint set $S$ of a canonical maximisation/minimisation linear programming program is bounded, then the maximum/minimum value of the objective function is attained at an extreme point of $S$.	
\end{theorem}

\begin{theorem}\label{T2} If the constraint set $S$ of a canonical maximisation/minimisation linear programming problem is unbounded, then there exists some $M\in\mathbb{R}$ such that the objective function $f$ satisfies $f(x_1, x_2, \ldots, x_n)\lessgtr M$ for all $(x_1, x_2, \ldots, x_n)\in S$, i.e. $f$ is bounded above/below (by $M$), them the maximum/minimum value of the objective function is attained at an extreme point of $S$. 	
\end{theorem}

\end{changemargin}

\subsection{A Geometric Approach}

Let us illustrate the above concepts with a typical example of resource allocation problem.

\begin{example} A furniture company manufactures sofas and armchairs. The production of one sofa requires 2 hours in the parts division and 1 hour on the assembly line of the company. The production of one armchair requires 1 hour in the parts division and 2 hours on the assembly line. The parts department of the company operates at most 8 hours per day, while the assembly division is active at most 10 hours per day. Knowing that the profit of selling one sofa is \pounds 30 and the profit of selling one armchair is \pounds 50, what are the optimal quantities of sofas and armchairs that the company should produce in order to maximise profits?	
\end{example}

We start by translating the problem in mathematical terms. Note that we are interested in the number of sofas and armchairs to be produced, so let us denote:

\[\begin{array}{lll}
x_1 & = & \text{\# of sofas per day;} \\ 
x_2 & = & \text{\# of armchairs per day.}
\end{array}
\]

According to the above definitions, $x_1$ and $x_2$ are the \textbf{decision variables} of our problem. Next, we wish to maximise the profits of the company, namely: 

\[f(x_1, x_2) = 30 x_1 + 50 x_2.\]

This is the \textbf{objective function} that we want to optimise. But we cannot have unlimited quantities and hence infinite profits, since the company is constrained by the availability of the parts and assembly divisions. We observe that 2 hours are spent in the parts division for one sofa and 1 hour is spent in the same division for one armchair. Together with the constraint of 8 hours per day in this department, we derive the constraint:

\[2x_1 + x_2 \le 8.\]

Similarly, we can derive the other constraint for the assembly line, which gives a second inequality:

\[x_1 + 2x_2 \le 10.\]

Finally, we also have the obvious constrains:

\[
\begin{array}{c}
x_1 \ge 0, \\x_2 \ge 0. 	
\end{array}
\]

since the company cannot produce negative quantities of sofas or armchairs. The above 4 constraints represent the \textbf{feasible set} that enforce the admissible production plans $x = (x_1, x_2)$. Note that, for example, $x = (20, 30)$ cannot belong to this feasible set, as there is not enough availability in both the parts and assembly divisions to satisfy that level of production. The final optimisation problem is a maximisation problem that can be formulated mathematically as follows:

\begin{figure}[H]
	\centerline{\includegraphics[scale=0.5]{pic1}}
\label{figure1}
\caption{Feasible set of points satisfying all constraints}
\label{figure1}
\end{figure}

\[
\begin{array}{ll}
\text{Maximise}   & f(x_1, x_2) = 30x_1 + 50x_2 \\
\text{subject to} & 2x_1 + x_2 \le 8 \\
	              & x_1 + 2x_2 \le 10 \\
	              & x_1 \ge 0 \\
	              & x_2 \ge 0.
\end{array}
\]

The set of points $(x_1, x_2)$ satisfying the above constraints is given by the shaded region from Figure \ref{figure1}.
\\

\begin{figure}[!htb]
\minipage{0.5\textwidth}
	\centerline{\includegraphics[scale=0.35]{pic2}}
\endminipage\hfill
\minipage{0.5\textwidth}
	\centerline{\includegraphics[scale=0.35]{pic3}}
\endminipage
\end{figure}

\vspace{0.3cm}

\begin{figure}[!htb]
\minipage{0.5\textwidth}
	\centerline{\includegraphics[scale=0.35]{pic4}}
\endminipage\hfill
\minipage{0.5\textwidth}
	\centerline{\includegraphics[scale=0.35]{pic5}}
\endminipage
\caption{All individual constraints}
\label{figure2}
\end{figure}

The region in Figure 1 above was obtained by intersecting the four 2-dimensional regions determined by the constraints $2x_1 + x_2 \le 8$, $x_1 + 2x_2 \le 10$, $x_1 \ge 0$, and $x_2 \ge 0$ (Figure \ref{figure2}). Therefore, the question boils down to determining which point in the shaded region in Figure 1 maximises $f(x_1, x_2)$.
\\

Since the constraint set is bounded, we can use Theorem \ref{T1} to deduce that the maximum value of $f(x, y)$ is attained at an extreme point in Figure 1, namely at either $(0, 0)$, $(4, 0)$, $(0,5)$, or $(2,4)$. A simple analysis shows that the maximum is attained at $(2, 4)$ where $f(2, 4) = 260$, and as a result the company should produce $2$ sofas and $4$ armchairs in order to maximise their profits (and get $\pounds 260$ in return).

\subsection{Issues and Alternative Approaches}

Having seen the geometric method which allows us to reason about LP problems, we can already notice some of its limitations. In real-life LP problems, one would typically deal with potentially tens if not hundreds of variables. 
\\

For example, assuming $m$ constraints and $n$ decision variables, we will have at most

\[{m + n \choose n} = \frac {(m + n)!}{m!\cdot n!}\]

\noindent candidate extrema points to be tested. This clearly poses a problem and makes the visualisation of the constraint set impossible.
\\

Another more relevant problem concerns the cases when the constraint set is unbounded. In such scenarios, according to Theorem \ref{T2} we would also need to ensure that the objective function is bounded by above or below, depending on whether the LP problem at hand is a maximisation or a minimisation problem, respectively. This is largely problem dependant and can complicate the analysis significantly. One such example is given by the following LP problem:

\[
\begin{array}{ll}
\text{Maximise}   & f(x, y, z) = 2x + 3y + 4z \\
\text{subject to} &  y + 5x \le 10 \\
	              & 2y + 3z \le 15 \\
	              & x, y, z \ge 0
\end{array}
\]

Here there is no restriction on $x$ other than being positive, and as such for $x\to\infty$ we will have $f(x, y, z)\to\infty$, and hence the objective function is unbounded.
\\

Finally, we end this section on linear programming noting that several algorithms have been developed to deal with the issues that we have just noted. The famous \textit{simplex algorithm} is able to find optimal solution to LP problems without testing a large number of candidate extrema points. It is also able to detect edge cases where the constraint sets are empty and the objective functions are unbounded (see \cite{ref6}). Several other optimisations such as branch and bound and cutting plane algorithms can further be used to further enhance the solution finding to LP problems \cite{ref7}.

\section{CVXOPT Framework}

For the purposes of implementing the consistency properties as described in the introduction section, we decided to use an open-source framework named CVXOPT \cite{ref8}. The main reason for choosing this software package was that it allows straightforward development of linear/convex optimisation problems using the Python programming language. This framework can be used either via an interactive Python interpreter or can be integrated within existing Python modules.
\\

Previous experience with GNU's GLPK package \cite{ref9}, which uses a rather obfuscated syntax with many intricate constructs, made us choose an alternative package which can leverage the strengths of the Python programming language.

\subsection{Simple Example}

We can illustrate the strengths of the CVXOPT framework with a simple example. Given the following minimisation problem:

\[
\begin{array}{ll}
\text{Minimise}   & f(x_1, x_2) = 2x_1 + x_2 \\
\text{subject to} & -x_1 + x_2 \le 1 \\
	              & x_1 + x_2 \ge 2 \\
	              & x_2 \ge 0 \\
	              & x_1 - 2x_2 \le 4,
\end{array}
\]

\noindent we can retrieve the optimal solution vector $\mathbf{x}$ using the following set of commands:

\begin{lstlisting}
>>> from cvxopt import matrix, solvers
>>> A = matrix([[-1.0, -1.0, 0.0, 1.0], [1.0, -1.0, -1.0, -2.0]])
>>> b = matrix([1.0, -2.0, 0.0, 4.0])
>>> c = matrix([2.0, 1.0])
>>> sol=solvers.lp(c,A,b)

     pcost       dcost       gap    pres   dres   k/t
 0:  2.6471e+00 -7.0588e-01  2e+01  8e-01  2e+00  1e+00
 1:  3.0726e+00  2.8437e+00  1e+00  1e-01  2e-01  3e-01
 2:  2.4891e+00  2.4808e+00  1e-01  1e-02  2e-02  5e-02
 3:  2.4999e+00  2.4998e+00  1e-03  1e-04  2e-04  5e-04
 4:  2.5000e+00  2.5000e+00  1e-05  1e-06  2e-06  5e-06
 5:  2.5000e+00  2.5000e+00  1e-07  1e-08  2e-08  5e-08

>>> print(sol[`x'])
[ 5.00e-01]
[ 1.50e+00]

\end{lstlisting}

\section{Domain Theory}

We introduce basic domain theory knowledge that will become useful when dealing with the problem of consistency in a later section. We use \cite{ref5} and \cite{ref10} when referencing these results, unless otherwise stated.

\subsection{Introduction}

Domain theory was introduced in 1970 by Dana Scott as a mathematical theory of programming languages.  The basic idea of domain theory is to provide better and better approximations to an object by means of simple recursion. It has applications in the field of computer science, e.g. solving canonically fixed point equations or recursive equations of procedures and data structures.

\subsection{Main Results}

\begin{definition} A \textit{partial order} (or a \textit{partially ordered set} or \textit{poset}) $(D, \le)$ is a set $D$ together with a binary relation $\le$ which is:

\begin{enumerate}
	\item reflexive: $a\le a$,
	\item anti-symmetric: $a\le b \wedge b\le a \implies a = b$, and
	\item transitive: $a\le b \wedge b\le c \implies a\le c$.
\end{enumerate}

Most often, the binary relation of a partial order is written as $\sqsubseteq$. Then $a\sqsubseteq b$ can be interpreted as $a$ having less information than $b$.
\end{definition} 

\begin{definition} A subset $A$ of an ordered set $(P, \sqsubseteq)$ is an \textit{upper set} if $x\in A\implies y\in A$, for all $y\sqsupseteq x$. We denote by $\uparrow A$ the set of all elements above some element of $A$. For convenience, we abbreviate $\uparrow \{x\}$ as $\uparrow x$. The dual notions are lower set and $\downarrow A$ \cite{ref11}.
\end{definition} 

\begin{definition} A \textit{cpo} (\textit{complete partial order}) is a poset with a least element denoted $\bot$ such that every chain $d_0 \sqsubseteq d_1 \sqsubseteq d_2 \sqsubseteq \ldots$ has a least upper bound, denoted by $\bigsqcup_{i\ge 0} d_i$ or, more conveniently, $\bigsqcup_i d_i$. Hence, $\bigsqcup_i d_i$ gives precisely the total information contained in the chain (and no more).
\end{definition} 

\begin{definition} A non-empty subset $A\subset P$ is said to be \textit{directed} if for any pair of elements $a, b\in A$ there exists $c\in A$ such that $a\sqsubseteq c$ and $b\sqsubseteq c$.
\end{definition} 

\begin{definition} A \textit{directed complete partial order} (\textit{dcpo}) or a domain is a partial order in which every directed subset has a least upper bound (lub).
\end{definition} 

\begin{definition} A dcpo is said to be \textit{pointed} if it has a least element which is denoted by $\bot$ and is called bottom.
\end{definition}

\begin{definition} Given two elements $a$ and $b$ of a dcpo we say that $a$ is \textit{way-below} or \textit{approximates} $b$, denoted by $a\ll b$, if for every directed subset $A$ with $b\in A$ there exists $c\in A$ with $a\sqsubseteq c$.
\end{definition} 

\begin{definition} A basis of a domain $D$ is a subset $B\subset D$ such that for every element $x\in D$ of the domain the set $B_x = \{y\in B\ |\ y \ll x\}$ of elements in the basis way-below $x$ is directed with $x = \bigsqcup B_x$. 
\end{definition} 

\begin{definition} A dcpo with a (countable) basis is said to be an \textit{($\omega$)-continuous} domain. 
\end{definition} 

\begin{definition} A function $f : D \to E$ between dcpo's is said to be \textit{Scott-continuous} if and only if it is \textit{monotone} (i.e. $a\sqsubseteq b \implies f(a) \sqsubseteq f(b)$) and preserves' lub's of directed sets i.e. for any directed $A\subseteq D$, we have $f\left(\bigsqcup_{a\in A} a\right) = \bigsqcup_{a\in A} f(a)$. Moreover, if $D$ is an $\omega$-continuous dcpo, then $f$ is continuous if and only if it is monotone and preserves the lub's of increasing sequences (i.e. $f\left(\bigsqcup_{i\in\omega} x_i\right) = \bigsqcup_{i\in\omega} f(x_i)$, for any increasing $(x_i)_{i\in\omega}$).
\end{definition} 

\begin{definition} Let $D$ be a dcpo. A subset $A$ is called \textit{(Scott-)closed} if it is a lower set and is closed under suprema of directed subsets. Complements of closed sets are called \textit{(Scott-)open}; they are the elements of $\omega_D$, the \textit{Scott-topology} on $D$ \cite{ref11}.
\end{definition} 

\section{Consistency of two functions}

In this section we focus on the main idea of the proposed project, which involves the idea of \textit{consistency} of two given maps subject to some constraints. We begin by stating the main notations and terminologies used throughout this section. We then present the general results for the $n$-dimensional case and illustrate with examples for the particular 1D and 2D cases \cite{ref1}, \cite{ref2}.

\subsection{Notations and terminology}

We denote by $\mathbb{R}$ the set of real numbers and by $\mathbf{I}\mathbb{R} = \{[a, b]\ |\ a\le b\in\mathbb{R}\}\cup\{\mathbb{R}\}$ the interval domain, i.e. the set of compact, nonempty intervals, equipped with a least element $\bot = \mathbb{R}$, ordered by reverse inclusion. It has a canonical basis consisting of all compact intervals with rational end points augmented with $\bot$. We write a non-bottom element $v\in\mathbf{I}\mathbb{R}$ as $v = [v^-, v^+]$ and we identify any real number $x\in\mathbb{R}$ with the singleton $\{x\}\subset\mathbb{R}$.
\\

Also, we denote by $\mathbf{I}\mathbb{R}^n$ the product domain consisting of all non-empty compact hyper- rectangles with faces parallel to the standard coordinate planes ordered with reverse inclusion and augmented with the whole space $\mathbb{R}^n$ as the bottom element. It has a canonical basis consisting of all its rational (compact) hyper-rectangles and the bottom element. We denote the continuous Scott domain of the nonempty, compact and convex subsets of $\mathbb{R}^n$, taken together with $\mathbb{R}^n$ as the bottom element and ordered by reverse inclusion, by $\mathbf{C}\mathbb{R}^n$. We will use a canonical basis of $\mathbf{C}\mathbb{R}^n$, consisting of rational convex compact polyhedra together with the set $\mathbb{R}^n$ as the bottom element.
\\

For an open subset $U\subset\mathbb{R}^n$, let $C^0(U)$ be the function space of all continuous functions of type $U\to\mathbb{R}$. We will also use domains of function spaces of the form $(U\to D)$ where $D$ is a countably based continuous dcpo, which is either $\mathbf{I}\mathbb{R}$, $\mathbf{I}\mathbb{R}^n$ or $\mathbf{C}\mathbb{R}^n$ in our case. For the sake of convenience, denote $D^0(U) = U\to\mathbf{I}\mathbb{R}$. A function $f\in D^0(U)$ is given by a pair of respectively lower and upper semi-continuous functions $f^-, f^+:U\to\mathbb{R}$ with $f(x) = [f^-(x), f^+(x)]$ when $f(x)\ne\bot$ for all $x\in U$. Recall that given an open subset $a\subset U$ and an element $b\in D$, the single step function $b_{\chi_a} : X\to D$ is deﬁned as $(b_{\chi_a})(x) = b$ if $x\in a$ and $\bot$ otherwise. Single-step functions are continuous with respect to the Scott topology. Any finite set of single-step functions that are bounded in the function space $U\to D$ has a least upper bound, called a step function; the set of step functions provides a basis for the continuous Scott domain $U\to D$. This basis in turn gives a countable and canonical basis of rational step functions for $U\to D$, where $D = \mathbf{I}\mathbb{R}$, $\mathbf{I}\mathbb{R}^n$ or $\mathbf{C}\mathbb{R}^n$, generated by single-step functions of the form $b_{\chi_a}$ where $a$ is a rational open hyper-rectangle with faces parallel to the coordinate hyper-planes of $\mathbb{R}^n$ and $b$ is a rational interval for $D = \mathbf{I}\mathbb{R}$, a rational hyper-rectangle for $D = \mathbf{I}\mathbb{R}^n$ and a rational compact convex polyhedron in $\mathbb{R}^n$ for $D = \mathbf{C}\mathbb{R}^n$.
\\

Finally, we introduce the following definition of \textit{interval Lipschitz constant} functions which generalises the well-known \textit{Lipschitz functions}:

\begin{definition}
	The continuous function $f:U\to\mathbb{R}$ has a \textit{non-empty, convex and compact set-values Lipschitz constant} $b\in\mathbf{C}\mathbb{R}^n$ in an open subset $a\subset U$ if for all $x, y\in a$ we have: $b\cdot (x-y) \sqsubseteq f(x) - f(y)$. The single step tie $\delta (a, b)\subseteq C^0(U)$ of $a$ with $b$ is the collection of all partial functions $f$ on $U$ with $a\subset \text{dom}(f)\subset U$ in $C^0(U)$ which have $b$ as non-empty convex compact set-values Lipschitz constant in $a$. 	
\end{definition}


\subsection{The property of consistency}

In simple terms, a pair of functions $(f, g)$, representing function and derivative information respectively, is called \textit{consistent} if one can find a third map $h$ which is approximated by the first component and whose derivative information is approximated by the second component of the pair. 
\\

The function information is given by a number of \textit{step functions} whose values are given by a finite set of pairs $(a_i, b_i)_{i\in I}$, where $a_i\subseteq\mathbb{R}^n$ is a rational hyper-rectangle and $b_i\subseteq\mathbb{R}$ is a compact interval such that $b_i$ and $b_j$ have non empty intersection whenever this is the case for the interiors of $a_i$ and $a_j$. The derivative information for the $n$ partial derivatives is given as a finite set of pairs $(a_i, b_i)_{i\in I}$, where $a_i$ are as above, but $b_j$ represent rational compact polyhedra. 
\\

In the particular case when the derivative constraints lie within hyper-rectangles with faces parallel to the coordinate planes, there is always a \textit{piecewise linear witness} for consistency. In addition, one can find the least and greatest piecewise linear maps consistent with the information from $f$ and $g$ via a linear programming algorithm as described in both \cite{ref1} and \cite{ref2}.
\\

We now formalise the \textit{consistency} relation under the definitions and notations used thus far:

\begin{definition}
	The consistency relation $\mathsf{Cons} \subset D^0(U)\times (U\to\mathbf{C}\mathbb{R}^n)$ is defined by:
	
	\[
		(f, g)\in\mathsf{Cons} \text{ if } \uparrow f\cap \textstyle\int g \ne \emptyset.
	\]
\end{definition}

For a consistent $(f, g)$, we think of $f$ as the \textit{function part} or the \textit{function approximation} and $g$ as the \textit{derivative part} or the \textit{derivative approximation}. It can be shown that the $\mathsf{Cons}$ is Scott closed \cite{ref2}. Also, in the one-dimensional case it holds that (Section 2, \cite{ref1}):

\begin{theorem} The following duality property holds:
\[h\in \int g\iff g\sqsubseteq \frac{dh}{dx}.\]
\end{theorem}

Furthermore, for a $n$-dimensional space, $n\ge 2$, we have that (Corrolary 2.8 from \cite{ref2}):

\begin{theorem} The following duality property holds: 
\[h\in \int g\iff g\sqsubseteq \mathcal{L}h.\]
\end{theorem}

Now, according to \cite{ref1} and \cite{ref2}, we have the least and upper consistency witnesses defined by the following:

\begin{proposition} Let $O$ be a connected component $\text{dom}(g)$ and let $R(U)$ be the set of partial maps of $U$ into the extended real line $\mathbb{R}\cup \{-\infty, \infty\}$. Consider the two dcpos $(R(U), \le)$ and $(R(U), \ge)$ having pointwise ordering inherited from the extended real line. Then the maps $s : D^0(O) \times (U\to \mathbf{C}\mathbb{R}^n)\to (R(U), \le)$ and $t : D^0(O) \times (U\to \mathbf{C}\mathbb{R}^n)\to (R(U), \ge)$ defined by:
	\[
	\begin{array}{c}
		\displaystyle s(f, g) = \inf \left\{h:\text{dom}(g) \to\mathbb{R}\ \left|\ h\in\int g, h\ge f^-\right.\right\} \\\\
		\displaystyle t(f, g) = \sup \left\{h:\text{dom}(g) \to\mathbb{R}\ \left|\ h\in\int g, h\le f^+\right.\right\}
	\end{array}
	\]
	
\noindent represent the least primitive map of $g$ that is greater than the lower part of $f$ and the greatest primitive map of $g$ that is less than the upper part of $f$, respectively.
	
\end{proposition}

\begin{corollary}
	The following 3 conditions are equivalent:
	\begin{enumerate}
		\item $(f, g)\in \mathsf{Cons}$;
		\item $s(f, g) \le t(f, g)$;
		\item There exists a locally Lipschitz function $h:\text{dom}(g)\to\mathbb{R}$ with $g\sqsubseteq\mathcal{L}h$ and $f\sqsubseteq h$ on $\text{dom}(g)$.  
	\end{enumerate}
\end{corollary}

Furthermore, the maps $s$ and $t$ are Scott continuous and the relation $\mathsf{Cons
}$ is Scott closed. 

\subsection{Consistency for the one-dimensional case}

Figure \ref{1d} shows the piecewise linear maps $s$ and $t$ which provide the least and greatest bounds for the consistency witness, with the given function and derivative approximations (the latter being constrained to rectangles in the Euclidean plane).   

\begin{wrapfigure}{l}{0.5\textwidth}
	\vspace{-30pt}
	\begin{center}
    	\includegraphics[scale=0.2]{pic7}
    	\caption{Consistency for the 1D case.}
    	\label{1d}
  	\end{center}
  	%\vspace{-30pt}
\end{wrapfigure}

A linear algorithm in the number of partitions induced by $(f, g)$ (i.e. the partition of points derived from intersecting all rectangle constraints, within the domain $[0, 1]$) is given in \cite{ref1}, Algorithm 3.3, to obtain the map $s(f, g)$ (a similar one computed the upper bound $t(f, g)$). 

\begin{algorithm} 
	The function updating algorithm consists of an initialisation step and two other main steps (see Figure \ref{1dAlgo}). The initialisation process determines the common partition points $\{y_0, \ldots y_n\}$ of $(f, g)$. On each interval $(y_{k-1}, y_k)$, the functions $g^-$ and $g^+$ are constant, with $g^-\vert_{(y_{k-1}, y_k)} = \lambda t.e_k^-$ and $g^+\vert_{(y_{k-1}, y_k)} = \lambda t.e_k^+$, where $e_k^-, e_k^+\in\mathbb{R}$. Furthermore, on each interval $(y_{k-1}, y_k)$, the map $f^-$ has a constant slope, $a_k$ say, i.e. $f^-\vert_{(y_{k-1}, y_k)} = f_k^-$, with $f_k^-(x) = a_kx + b_k$.
\end{algorithm}

\noindent \textbf{Input:} $f, g :[0, 1]\to\mathbf{I}\mathbb{R}$ where $f$ is a linear step function and $g$ is a step function.

\noindent \textbf{Output:} Continuous function $s(f, g) : [0, 1]\to\mathbf{I}\mathbb{R}$ which represents the least function consistent with the information from $f$ and $g$.
\\

\noindent \textbf{Initialisation:} 

\noindent $\{y_0,\ldots y_n\}$ \lit*{:  induced-partition-of} $(f, g)$ 

\noindent \textbf{Part 1:}

\noindent $u(y_0) := f^-(y_0^+)$

\noindent \lit*{for } $k = 1\ldots n$ and $\forall x\in [y_{k-1}, y_k)$

\noindent \hspace{0.5cm} $u(x) := \max\{f^-(x), u(y_{k-1}) + (x - y_{k-1}) e_k^-\}$

\noindent \hspace{0.5cm} $u(y_k) := \max\{\overline{\lim} f^-(y_k), u(y_{k-1}) + (y_k - y_{k-1}) e_k^-\}$

\noindent \textbf{Part 2:}

\noindent $s(y_n) := u(y_n)$

\noindent \lit*{for } $k = n\ldots 1$ and $\forall x\in [y_{k-1}, y_k)$

\noindent \hspace{0.5cm} $s(f, g)(x) := \max\{u(x), s(y_k) + (x - y_{k-1}) e_k^+\}$

\begin{figure}[H]
	\centerline{\includegraphics[scale=0.5]{pic8}}
	\caption{The function updating algorithm} 
	\label{1dAlgo}
\end{figure}

\subsection{Consistency for the $n$-dimensional case, $n\ge 2$}

The framework for studying consistency for an $n$-dimensional setup ($n\ge 2$) is much more challenging than what we have just seen already. In this section we describe the algorithm from \cite{ref2} which works for the general case $n\ge 2$. However, for the purposes of presentation, we will only deal with the case when $n = 2$, which is easier to reason about (from a geometric perspective at least). Notice that, as with the 1D case, the derivative constraints are given as hyper-rectangles, to simplify the computational framework.
\\

Figure \ref{2d} shows two examples of consistent tuples, namely the least and greatest witnesses that enforce the consistency property described so far. In the example from the left hand side, there is one hyper-rectangle for the function approximation, while the derivative approximation in the $x$ and $y$ directions is given my the constant intervals $[n, N]$ and $[m, M]$, where $n, m > 0$. In the other instance, on the right, there are two intersecting hyper-rectangles for the function approximation and the derivative approximations are the constant intervals $[0, 0]$ and $[m, M]$, where $m > 0$.

\begin{figure}[H]
	\centerline{\includegraphics[scale=0.27]{pic6}}
	\caption{Consistency for the 2D case}
	\label{2d}
\end{figure}

We now proceed to explaining the framework, as presented in \cite{ref2}, Section 3. Let us introduce the following linear programming algorithm:

\begin{algorithm} In a similar fashion to the induced partitions for the 1D case, we impose a grid $(p_0, \ldots p_k)\times (q_0, \ldots q_l)$ on the unit square $U$ such that the function approximation given by the step function $f : U\to\mathbf{I}\mathbb{R}$ and the derivative approximation given by the step function $g : U\to\mathbf{C}\mathbb{R}^n$ are constant respectively with values $c_{ij}\in\mathbf{I}\mathbb{R}$ and $b_{ij}\in\mathbf{C}\mathbb{R}^n$ inside every sub-rectangle $(p_i, p_{i+1})\times (q_j, q_{j+1})$, for $i=\overline{0, k-1}$ and $j = \overline{0, l-1}$, defined by adjacent grid points. Note that if $c_{ij} = \bot$ or $b_{ij} = \bot$ then $f$ or $g$ are undefined in the sub-rectangle and their contribution can be ignored in the following analysis.

If the pair of step functions $(f, g)\in (U\to\mathbf{I}\mathbb{R})\times (U\to\mathbf{C}\mathbb{R}^n)$ is consistent then we can find values $h_{i, j}\in\mathbb{R}$ at all the grid points $(p_i, q_j)$ for $i = \overline{0, k}$ and $j = \overline{0, l}$, such that for $i = \overline{0, k-1}$ and $j = \overline{0, l-1}$ they satisfy the following two conditions:

\begin{enumerate}
	\item $c_{ij}^- \le h_{s, t} \le c_{ij}^+$, for $s = i, i + 1$ and $t = j, j + 1$;
	\item $c := \displaystyle\left(\frac {h_{i+1, j}-h_{i, j}}{p_{i+1} - p_i}, \frac {h_{i, j+1}-h_{i, j}}{q_{j+1} - q_j}\right)\in b_{ij}^1 \times b_{ij}^2 = b_{ij}$.
\end{enumerate}

\end{algorithm}

\chapter{Project Plan}

This section provides information on the next steps of the project, and the order in which the implementation tasks will be handled. Given the reading I have carried out so far, as detailed in the background section above, I expect to complete the project in the following stages:

\begin{enumerate}
	\item Start with the 1D case for consistency, when derivative constraints lie within rectangles and implement the linear programming algorithm as presented in section 3.3 from \cite{ref1};
	\item Move on towards the $n$D case for consistency, $n\ge 2$, when the restrictions for the derivative information are given by compact hyper-rectangles with faces parallel to the coordinate planes. The algorithm for this more generic case is presented in \cite{ref2};
	\item Implement the 2D case for consistency, when the constraints lie in convex (hyper)-quadrilaterals, for which a known linear programming algorithm exists;
	\item Finally, investigate if the latter case can be generalised to convex polygons/polyhedra constraints and whether the result still holds in higher dimensional spaces.
\end{enumerate}

Below is a table giving the timeline I expect to follow for each of the above-mentioned tasks, until the project deadline. Any incomplete work by the Easter break will need to be completed during the holiday, to ensure the completion of all milestones of the project.

\begin{center}
\begin{tabular}{|l|l|p{5.3cm}|}
\hline
\multicolumn{1}{|c|}{\textbf{Term \& Weeks}} & \multicolumn{1}{|c|}{\textbf{Dates}} & \multicolumn{1}{|c|}{\textbf{Description}} \\ \hline \hline
Term 2, Weeks 4-5 & 1 FEB - 14 FEB & Starting first task (consistency for the 1D case)                                                                                                                                                               \\ \hline
Term 2, Weeks 6-7 & 15 FEB - 28 FEB & Proceeding to the second task (consistency for the n-dimensional case, $n\ge 2$)                                                                                                                                                     \\ \hline
Term 2, Weeks 8-9 & 29 FEB - 13 MAR & Moving on to a more general case as part of task 3 (consistency for the 2D case with convex quadrilateral constraints for derivatives)                                                                                                        \\ \hline
Term 2, Weeks 9-11 & 14 MAR - 27 MAR & Exams Preparation                                                                                                                                                                                                                                                  \\ \hline
Easter Break & 28 MAR - 24 APR & Code tweaking, refining and catching up on any implementation that is left incomplete                                                                                                                           \\ \hline
Term 3, Weeks 1-4 & 25 APR - 22 MAY & Investigating the most challenging generalisation as described in task 4 (consistency in higher dimensions, where derivative constraints are specified as convex polyhedra, rather than hyper-rectangles/quadrilaterals) 
    \\ \hline
Term 3, Weeks 5-7 & 23 MAY - 12 JUN & Consolidating final report, findings and evaluation                                                                                                                                                                      \\ \hline
Term 3, Weeks 8-9 & 13 JUN - 24 JUN & Final report submission and presentation                                                                                                                                                                                                                       \\ \hline
\end{tabular}
\end{center}

\chapter{Evaluation Plan}

In terms of evaluation, the implementation of the linear programming algorithms must be verified against a relevant and thorough test-suite. After each stage of the project, a number of test-suites will be written so as to validate the work and make sure I can proceeded confidently to the next challenges.
\\

The plan is to also validate some of the results presented in \cite{ref1}, \cite{ref2} in order to have a satisfying product by the end of the project. In addition, it would probably be worthwhile to have some sort of GUI to demonstrate with relevant examples the implementation of the linear programming algorithm. This will be useful, for example, at the end of task 2, when we want to find the minimal and maximal bounding surfaces, once the consistent witness is determined.
\\

In addition, when the linear programming framework will be fully implemented at the end of stage 3, it should become relatively straightforward to extend it in order to assess whether the consistency property holds in more general settings. Assuming the results hold for convex shaped objects and a consistent witness can be determined, we will be able to conjecture or even prove our findings. Should counterexamples be found, we can be certain that the problem of consistency will not hold for those specific cases.
\\

Finally, the framework to be developed should be convenient to use to allow further extensions to be carried out easily in the future. This could include, as mentioned before, some graphical additions to allow the visualisation of the piecewise linear maps. In order to enhance the overall use of our framework, we may also consider ways to add more constraints at runtime or modify input parameters on-the-fly and have the results displayed in real-time.








