Why does the redshift formula hold true when measuring the expansion of the universe? As far as I understand, photons get 'stretched' as the universe expands, and the expansion rate is not constant. This means that the redshift observed on Earth is the sum of individual redshifts arising from the expansion of space at each instant. Therefore, I'd think the 'velocity' of the object we are looking at at the moment the observed light was emitted should be proportional to $$λ_0\int a(t)λ(t)dt$$ (where $λ_0$ is the original wavelength of the emitted light, $λ(t)$ is the wavelength at time $t$, and $a(t)$ is the scale factor at time $t$), with $a(t)$ behaving in a manner that we are trying to find. If that were the case, it would render drawing any reasonable conclusions based on redshift data practically impossible, as it would create a logical loop (in order to find $a(t)$, we need to know how it changes with time, which we can only determine by knowing what $a(t)$ is). Where am I wrong?
 A: In the following, I will adopt the notation of $a(t)$ for the scale factor, and take $a(t_0) = 1$ as an arbitrary unit for it. Then the expansion causes a redshift given simply by
$$ z = \frac{R(t_0)}{R(t_e)} - 1 = \frac{1}{a(t_e)} - 1, $$
where $t_e$ is the time of emission of the observed photon.
On the other hand, computing the present-day distance to the object that emitted the photon requires an integral over the whole time of flight:
$$
d = \int_{t_e}^{t_0} \frac{cdt}{a(t)} $$
because the distance increments $cdt$ travelled at some time in the past have subsequently "stretched" by $1/a(t)$; equivalently, we can think of this as calculating the distance travelled in co-moving coordinates, which with the convention I adopted ($a_0 = 1$) are the same as current coordinates. All this you seem to already know, but I include it for completeness.
As Cristoph's answer already suggested (I'll try not to repeat too much of it), if we assume a particular expansion model, we can now find redshift and distance, both as functions of $t_e$; if we then eliminate $t_e$ from these two equations, we get distance as function of redshift. No-one is claiming that it has to be linear in general. But for any reasonable expansion model, it will be linear for small enough distances/recent enough emission times. So, first question: what is, after all, "not too far"/"not too long ago?"
Well, that depends on the model. You've already calculated that for linear expansion we get a logarithmic curve: if $a(t) = 1 + H_0(t - t_0)$, then $d = (c/H_0)\ln(1 + z)$. The approximation $\ln(1+z) \approx z$ holds for $z \ll 1$, which in turn implies $|a-1| \ll 1$, and then $|t - t_0| \ll 1/H_0$. With an exponential model, $a(t) = e^{H_0 t}$, one finds $d = (c/H_0)z$: that one turns out to be linear for all distances/times. Any model I've come across, linearity inevitably ends up being a valid approximation for small redshifts $z \ll 1$, whatever that means in terms of time/distance, and this is why cosmologists often treat redhsift itself as a (model-independent) measure of "nearness" (in an abstract sense). But for any particular model, we do have a theoretically predicted redshift-distance relation.
Which brings us to the second part of your question, how do we know which model is the right one, so that we can use it to deduce distance from redshift alone. If we only had redshift as observable data, it would be - as you say - a logical loop. But we do have other things to work with. Cristoph touched on that, too: first we need to judge our models based on other input. Theoretically, the expansion is driven by the distribution of energy in the universe. While there are still a lot of open questions in that respect, we do know a little bit. And, very importantly, we have other ways to infer something about distance. At such large scales, usually they involve arguments of homogeneity and isotropy of the universe.
For example, the mean density of "stuff" in the universe is assumed to be mostly constant. That means that if we constructed a cumulative frequency function for number of observed galaxies up to a given redshift, we would expect that number to be proportional to the volume of a sphere whose radius is the distance corresponding to that redshift. So we can get distance through volume. If we assume the mean brightness of a galaxy is also mostly the same everywhere, we can use apparent brightness to deduce distance, this time through area (assuming all directions equivalent, light from a point source spreads out evenly over the surface of that same sphere, so intensity is inversely proportional to its surface area). Finally, as Cristoph mentioned, a direct radial estimate might come through time itself, if we expect mean age of stars or galaxies to be about the same everywhere (information about age at time of emission is present in the spectral properties of the light). These are just the ones that come to my mind, there may be other methods. Interestingly, since they all give a different "type" of distance (geometrically speaking), they allow us to infer something about curvature as well, which plays back into the metric equations and thus links back to the theoretical foundation of the models. 
So all these considerations are interconnected, and the thought process is not just linear from observation to conclusion. We come up with models, explore their implications, compare to observations, if necessary adjust the models, and do this with many models in parallel. The scientific process is often one of guessing and fitting; much like a young kid, who knows nothing about equations, can still solve a riddle like "what number plus 5 is equal to twice that number?"
A: In the Friedmann model, proper distances at constant cosmological time evolve according to
$$ r = \frac a{a_0} r_0 $$
The same goes for wavelengths, ie
$$ 1 + z = \frac \lambda{\lambda_0} = \frac a{a_0} $$
which you can also derive kinematically by parallel-transporting momentum vectors along null geodesics.
This means cosmological redshifts directly correspond to scale factors. To experimentally measure the time evolution of the scale factor, you'd then have to determine the age of the object in some other fashion (perhaps the age of the oldest stars in a given galaxy?). On the flip side, if you have fixed the parameters of your cosmological model and computed the evolution of the scale factor by solving Einstein's field equations, you can then determine the age of an object from its redshift.
Hubble's law relating velocities and distances follows directly from the first equation:
$$ \dot r = \frac {\dot a}{a_0} r_0 = \frac {\dot a}a r \equiv H r $$
However, the observed linear relation between redshift and distance is only an approximation:
$$ 1 + z \approx \frac{a + \dot a \Delta t}a = 1 + H \Delta t$$
and thus
$$ z \approx H \Delta t \approx H \Delta r c^{-1} $$
under the additional assumption $\Delta r \approx c \Delta t$.
Note that historically, this has been interpreted in terms of Doppler shifts (eg in Hubble's famous paper, even though Lemaître had already published his interpretation). Such an interpretation also yields recession velocities
$$ v \approx H \Delta r $$
as special-relativistic Doppler shifts reduce to $z \approx \frac vc$ if $v \ll c$.
