First of all I don't see how this estimates the local error defined as the error we make in a single step using correct previous values. What does a lower order method have anything to do with previous correct values? I can accept that this is the error we make using a lower order method instead of a higher order method, but what is the relevance of this?
Both methods are converging to the same true solution, but at slightly different weights. You can think of convergence of rate $p$ as matching the first $p$ terms of the Taylor series. The is because when proving convergence to the true solution, you expand $u(t_{n+1})$ in its Taylor series of $h = t_{n+1} - t_n$ and show that you match the first $p$ terms. So given that, in Taylor series form, $r_{n+1}$ is some polynomial $r_{n+1} = \alpha_p h^p + \mathcal{O}(h^{p+1})$ since both of the methods have exactly the same expansion for the first $p-1$ terms, but only one of the methods matches the $p$th term. Dropping the higher order terms since $h$ is considered small enough, we have $r_{n+1} \approx \alpha_p h^p$. But understand what this term means. Since the $p$th order method has the correct Taylor series expansion, this term is exactly the error of the $p-1$th order method from the true Taylor term, meaning that $\alpha_p$ is the coefficient of the leading (largest) term in the local truncation error of the $p-1$th method and thus a good approximation for the method's local truncation error.
So now we have two things, $r_{n+1}$ is a good approximation to one of the method's local truncation error, and $r_{n+1}\approx \alpha_p \Delta t^p$. The true path to the formula you have up there is a result for the optimal stepsize due to Cechino where you prove that would give you the error/time optimal stepsize in the next step if the error estimate is exact. This is also known as proportional control but doesn't take into account stability concerns so this is not the method actually used in a lot of ODE solver software, but it's a good teaching tool. You can get to this result without the optimal stepsize proof just through a bit of reasoning though. Proportional control makes sense because what it means is that you want to correct the stepsize proportional to to size of the error: the larger the error, the more correction you want. So you could take your $TOL/r_{n+1}$ as the coefficient, that's TOL divided by your error estimate, meaning that if your error estimate is twice the tolerance, you would halve your stepsize. That's reasonable, but recall that $r_{n+1}\approx \alpha_p \Delta t^p$, so dimensionally it makes sense to take the $p$th root so that way this term actually scales like $1/h$. Of course, this is heuristic and the actual proof is more in depth. So there you go, now you have a reasonable error estimate and a reasonable way to make that produce a new stepsize.
But I don't understand why this would work since we update our step size so that $r_{n+1}=TOL$, this would just result in $h_{n+1}=h_n$ after a while.
If the error of the methods were independent of the stepsize or the location in time for the ODE, this would happen. But that's never true. Instead, the local truncation error of the methods very depending on many factors, such as the derivatives of the derivative function $f$, the size of the current state values, etc., and thus $\alpha_p$ is always changing, causing $r_{n+1}$ and the optimal stepsize to always change. But for some simple enough equations, when adaptivity is working correctly, the stepsize will stability around some "optimal" value.
