Introduction
The first section below can be found in standard numerical analysis textbooks. Most current textbooks seem to assume a working environment such as MATLAB or a programming language such C, Python, etc. in which there are at most two choices for working precision: single and double. As a result, WorkingPrecision as a parameter is not much discussed in textbooks and could be considered an issue peculiarly related to Mathematica. I explain its role in establishing reliable convergence criteria, which are set by PrecisionGoal and AccuracyGoal, in the last part of this answer.
###Stability and conditioning
$\def\x{{\bf x}}\def\y{{\bf y}}\def\A{{\bf A}}$
Suppose $\A$ is an algorithm that for an input $\x$ computes an output $\y$ that satisfies a condition
$$F(\x,\y)=0\,.$$
Let us use the same symbol $\|\cdot\|$ to denote "norms" that measure the "size" of $\x$ and $\y$.
We say $\A$ is stable if for a sequence $\x_n\rightarrow\x$, the sequence of corresponding outputs $\y_n\rightarrow\y$.
A measure of stability of a solution $(\x,\y)$ is given by the condition number $K$, which we roughly define as the ratio of the relative change in output to the relative change in input, over all sufficiently small changes in input.
Somewhat more formally,
$$K = \sup_{\Delta\x}\left.{\|\Delta\y\|\over\|\y\|}\right/{\|\Delta\x\|\over\|\x\|} \,,$$
where the supremum is taken over all sufficiently small $\Delta\x$, and $\Delta\y$
corresponds to the solution to $F(\x+\Delta\x,\y+\Delta\y)=0$.
What is considered sufficiently small, just like the norm itself, is context-dependent. The condition number $K$ bounds the factor that multiplies the relative error in the input, the product being equal to relative error in the output.
For a continuous model over the real or complex numbers, the above is probably a sufficient formulation. A floating-point implementation $\A^*$ for a given $\x$ usually computes a solution $\y^*=\y+\Delta\y$ that can be viewed either as a solution corresponding to an input $\x^*=\x+\Delta\x$ satisfying
$$F(\x+\Delta\x,\y+\Delta\y)=0\tag{1}$$
or approximately satisfying the condition, i.e.
$$F(\x,\y+\Delta\y)+\Delta F=0\tag{2}$$
or some other combination of deltas.
It seems to be difficult to give a general definition numerical stability.
The discretization of floating point numbers means there are only finitely many inputs, which complicates the consideration of limits and bounds.
At the least in a numerically stable algorithm, the maximum error $\sup{\|\Delta\y\|}$ would decrease as $\sup\|\Delta\x\|$ decreases.
Even more, one might require the condition number $K$ to be bounded.
One might even want $K$ to be relatively small, which we will discuss below in relation to the working precision. This criterion serves to distinguish well-conditioned problems ($K$ small) from ill-conditioned ones ($K$ large).
WorkingPrecision, PrecisionGoal and AccuracyGoal
WorkingPrecision is fairly straightforward, but
PrecisionGoal
and
AccuracyGoal
are used in different ways in different functions, corresponding roughly to the viewpoints of error in equations (1) and (2) above.
WorkingPrecision selects either MachinePrecision
or sets the number of digits for arbitrary-precision numbers
to be used in the computation.
Currently (2018), MachinePrecision is most commonly what is called (64-bit) double precision; it corresponds to a setting of $MachinePrecision, currently about 15.95, but without the precision-tracking and extra guard digits of arbitrary-precision numbers.
Setting WorkingPrecision has the side effect in most numerical functions of setting the default values for PrecisionGoal -> Automatic and/or AccuracyGoal -> Automatic to half the working precision.
PrecisionGoal and AccuracyGoal are input parameters that are used together to determine a stopping condition in terms of a bound or bounds on the truncation error. As @ilian remarks, they are goals only, since numerical error estimators, which are subject to their own assumptions and errors, can be wrong.
One type of stopping criterion for PrecisionGoal ->pg and AccuracyGoal -> ag, discussed in this answer and Norms in NDSolve, is
$$\|\Delta\y\| < 10^{-ag} + 10^{-pg}\,\|\y\|\tag{3}$$
When $10^{-pg}\,\|\y\| <\mskip-4mu< 10^{-ag}$, the approximation
$$10^{-ag} + 10^{-pg}\,\|\y\| \approx 10^{-ag}$$
implies that if the error indicates that $\y$ is correct to at least $ag$ digits past the decimal point, the algorithm terminates.
And when $10^{-ag} <\mskip-4mu< 10^{-pg}\,\|\y\|$, the approximation
$$10^{-ag} + 10^{-pg}\,\|\y\| \approx 10^{-pg}\,\|\y\|$$
implies that if the error indicates that $\y$ is correct to at least the leading $pg$ digits, the algorithm terminates. Thus depending on how $\|\y\|$ varies,
the stopping criterion can be dominated by precision or accuracy.
This sort of criterion is used in NDSolve and NIntegrate.
The other type of stopping criterion is
$$\|\Delta\y\| < 10^{-pg}\,\|\y\|
\quad\text{and}\quad
\|\Delta F\| < 10^{-ag}\tag{4}$$
This sort is used by FindRoot and related functions. In FindRoot[f[y] == x, {y, y0}], the "accuracy" error $\|\Delta F\|$ is Norm[{f[y]} - {x}].
The relationship between WorkingPrecision and error goals.
Suppose for the condition number and error goals, we have $K = 10^8$ and PrecisionGoal -> 8, AccuracyGoal -> 8. At a working precision (WP) of double precision (DP), or the usual MachinePrecision, we have a machine epsilon, or
$MachineEpsilon, of
$$\epsilon_{WP} = \epsilon_{DP} = 2^{-52} \approx 2 \times 10^{-16}\,.$$
In other words there are nearly 16 digits of precision, or an uncertainty in the input of $\epsilon_{WP}/2$.
The condition number $K$ implies we could lose up to 8 of the digits of precision, leaving around 8 in the computed answer $\y$.
With luck, it will satisfy the precision goal, but we shouldn't be surprised if sometimes it does not. We should not expect that raising PrecisionGoal will improve the result at all. If we're lucky, Mathematica will complain that the error goals cannot be met (e.g. the messages FindRoot::lstol, NIntegrate::slwcon, etc.).
Whether raising AccuracyGoal has any effect depends on how the error is measured and on the size of $\|\y\|$.
For error measured by (3), we have the following.
If $\|\y\| \ge 1$, then AccuracyGoal in effect sets a relative precision goal of at least $ag$, which will not help when raised above the setting of PrecisionGoal. If $\|\y\| <\mskip-4mu<1$, then PrecisionGoal -> 8 implies more than 8 accurate digits after the decimal point,
but AccuracyGoal -> 8 means the algorithm stops at 8;
therefore raising AccuracyGoal could be expected to improve the solution.
For error measured by (4), the effect of raising AccuracyGoal depends on $F(\x,\y)$.
The question of whether for a given solution $\y^*$ there is an update $\delta\y$ so that $F(\x,\y^*+\delta\y)+\delta F=0$ and $\|\delta F\| < \|\Delta F(\x,\y^*)\|$.
For instance, if the norm of the Jacobian of $F$ is too large, it may not be possible at the specified working precision to find an update $\delta\y$ such that $\y + \delta\y \ne \y$ and the error is reduced. In such a case, raising the accuracy goal would not help (see this example).
For both measures, AccuracyGoal can be seen as setting a limit $10^{-ag}$ below which numbers may be assumed to be equivalent to zero. For instance, NIntegrate has as a default AccuracyGoal -> Infinity, which causes it to complain when a numerical integral converges to, but does not exactly equal, zero (NIntegrate::slwcon). Setting AccuracyGoal to a finite value tells how small the approximation needs to be in order to be considered close enough to zero.
Here is another way to look at it:
At a fixed working precision, there are only finitely many possibilities for $\y$.
Among them, there is one with least error.
If that error meets or exceeds the error bounds set by PrecisionGoal and AccuracyGoal, then adjusting the goals will not result in a better solution.
Further, if the conditioning is bad, then the error may vary wildly in the neighborhood of the true solution, making finding the optimal solution like searching for a needle in a haystack.
If, on the other hand, the least possible error is much less than the error bounds,
then increasing PrecisionGoal and AccuracyGoal can result in better solutions.
Increasing working precision increases the number of possibilities for $\y$. It also ameliorates the effects of poor conditioning. For instance, if we set WorkingPrecision -> 24 and $K = 10^{8}$, we might expect around 16 digits of precision in the computed solution. Then raising the goals would have an effect.
In practice, there are other issues. One does not know the optimal answer, sometimes not even approximately. Sometimes one does not have a good idea about the conditioning. (For some problems, you can evaluate the objective function at arbitrary precision numbers and see how many digits are lost, the number lost giving a lower bound on how much to raise WP.) Some of the analysis above assumes the algorithm starts close to the true answer (i.e., the initial error is sufficiently small). Error goals and working precision cannot solve these issues (except maybe for the conditioning).
Here are three reasons for raising WorkingPrecision (the first two were mentioned here):
One cannot achieve the precision and accuracy goals at the current WP. In such a case, raise WP but also manually set PrecisionGoal and AccuracyGoal to their current values (don't raise them).
One wants to raise PrecisionGoal and AccuracyGoal, and a higher WP is needed.
One wants to test numerical stability: If the result is stable under higher and higher settings of WP, then it is evidence of numerical stability at the initial WP.
WorkingPrecisioneffectively sets the maximum accuracy and precision that can be obtained in a computation. However,AccuracyGoalandPrecisionGoalcan be used to obtain a result with a lower accuracy and precision.. This could be useful in a slow computation with, for instance,FindRoot, if you need accuracy and precision less than half theWorkingPrecision. LowAccuracyGoalandPrecisionGoalwould terminate the iterative computation sooner than it otherwise would, saving computer time. $\endgroup$WorkingPrecisionis exactly what it says on the tin: the internal algorithms are effectively fed numbers whose precision is set to that particular setting. This is useful for potentially unstable computations, where it is possible to get a result with, say, only 5 digits of accuracy even when the starting inputs have a precision of 30. $\endgroup$PrecisionGoal/AccuracyGoalshould be thought of as convergence criteria and can have different meaning (usually documented) for different functions. They are indeed just goals, and whether they will be satisfied depends on both the algorithm and the input. Typically a major required ingredient is a sufficiently highWorkingPrecision. In fact, the automatic setting of these is often (heuristically) based on the working precision, e.g. forMachinePrecisioninput the default precision goal forNIntegratewould be around 6, forNDSolvearound 8 etc. But your specific mileage may vary. $\endgroup$