I'm trying to algorithmically solve the Time Delay of Arrival problem as part of some mathematics research. The problem is as follows:

Given the location of three receivers in a plane (A, B, and C), the velocity of some signal, and the time at which each receiver "saw" the signal, determine the location of the source. We can assume the source is also in the plane, and we know nothing about it's location or distance, nor the time the signal was emitted.

Here's an excerpt of what I've found on the topic:

These linear system of equation [referring to the closed-form linear solution of the problem] are easily written in matrix form as below. (i = 1,2,3 .. N-1) \begin{gather} \begin{bmatrix} x_0-x_1 & y_0-y_1 & z_0-z_1 & d_{01}\\ x_0-x_2 & y_0-y_2 & z_0-z_2 & d_{02}\\ x_0-x_3 & y_0-y_3 & z_0-z_3 & d_{03}\\ ... & ... & ... & ...\\ x_0-x_n & y_0-y_n & z_0-z_n & d_{0n}\\ \end{bmatrix} \begin{bmatrix} x\\ y\\ z\\ d_0\\ \end{bmatrix} = \begin{bmatrix} \frac{1}{2}(x_0^{2} - x_1^{2} + y_0^{2} - y_1^{2} + z_0^{2} - z_1^{2} + d_{01}^{2}) \\ \frac{1}{2}(x_0^{2} - x_2^{2} + y_0^{2} - y_2^{2} + z_0^{2} - z_2^{2} + d_{02}^{2}) \\ \frac{1}{2}(x_0^{2} - x_3^{2} + y_0^{2} - y_3^{2} + z_0^{2} - z_3^{2} + d_{03}^{2}) \\ ... \\ \frac{1}{2}(x_0^{2} - x_n^{2} + y_0^{2} - y_n^{2} + z_0^{2} - z_n^{2} + d_{0n}^{2}) \\ \end{bmatrix} \end{gather} It's an ordinary $A\vec{x}=\vec{b}$ equation. Matrice A and vector $\vec{b}$ are known. The problem requires determination of $\vec{x}$ such that $A\vec{x} \approx \vec{b}$ minimizing sum of squares of residuals. It's called least square regression.

Assuming that error is $\varepsilon$ $$ \varepsilon(\vec{x}) = A\vec{x} - \vec{b} $$ Minimize square of error function $\parallel A\vec{x} - \vec{b}\parallel^2$ by deriving it with respect to $\vec{x}$ and then equate zero. \begin{equation} \begin{split} \frac{\partial}{\partial \vec{x}} \parallel A\vec{x}-\vec{b}\parallel^2 &=0\\ \frac{\partial}{\partial \vec{x}} [(A\vec{x}-\vec{b})^T(A\vec{x}-\vec{b})] &= 0\\ \frac{\partial}{\partial \vec{x}} [\vec{x}^T A^T A \vec{x} - 2\vec{x}^TA^T \vec{b} + \vec{b}^T \vec{b} ] &= 0\\ 2A^T A \vec{x} - 2A^T \vec{b} &= 0 \\ \Rightarrow \vec{x} &= (A^TA)^{-1} A^T\vec{b} \end{split} \end{equation} The vector $\vec{x}$ is solution.

If anchors are not placed uniformly then matrice $A$ is not singular so it can be invertible and solution exists, otherwise QR decomposition can be used. It's another topic.

Admittedly, this is *way* outside my area of expertise and it's been some time since I've had any linear algebra, etc.

I'm looking for a "more detailed" explanation, in simpler terms, of this information. What is this "error" ($\epsilon$), and how is it related to the solution? What's going on with that bit at the end ("minimize the square of the error function")?

EDIT: My source for the above quote is an answer to a similar question at math SE. Effectively the same information can be found in related research works, Wikipedia, etc.