Dynamic range of unknown parameters and form of cost function in scipy.optimize.least_squares
Most solvers are designed for variables in the 1-10 range. A large range can cause numerical problems, but it is not guaranteed to be problematic. Numerical problems sometimes stem from the matrix factorization step of the linear algebra for solving the Newton step, which is more dependent of the magnitude of the derivatives. You may also encounter challenges with termination tolerances for values outside the 1-10 range. Overall, if it looks like it's working, it's probably fine. You could get a slightly better answer by normalizing values.
Division by a degree of freedom can cause difficulties in three ways:

division by zero
discontinuous derivatives around 0
very steep derivatives near 0, or very flat derivatives far from 0

For these reasons, I would recommend \sum^N_{n=1} ( g_n - ( a0' e^(tn/a1) y_n - b0' e^-(tn*b1')) )^2. However, as previously stated, if it's already working it may not be worth the effort to reformulate your problem.
I am using scipy.optimize.least_squares to solve an interval constrained nonlinear least squares optimization problem.  The form of my particular problem is that of finding a0, a1, b0, and b1 such that the cost function:
\sum^N_{n=1} ( g_n - (y_n - b0 e^-(tn/b1)) / a0 e^-(tn/a1) )^2
is minimized where g_n, y_n and t_n are known and there are interval constraints on a0, a1, b0, and b1.
The four unknown parameters span approximately four orders of magnitude (e.g, a0 = 2e-3, a1 = 30, similar for b0 and b1).  I have heard that a high dynamic range of unknown parameters can be numerically problematic for optimization routines.   
My first question is whether four or so orders of magnitude range would be problematic for scipy.optimize.minimize. The routine appears to converge on the data I've applied so far. 
My second question relates to the form of the cost function.  I can equivalently write it as:
\sum^N_{n=1} ( g_n - ( 1/a0 e^(tn/a1) y_n - b0/a0 e^-(tn/b1) +tn/a1) / )^2
=
\sum^N_{n=1} ( g_n - ( a0' e^(tn/a1) y_n - b0' e^-(tn*b1')) )^2
where the new parameters are simple transformations of the original parameters.  Is there any advantage to doing this in terms of numerical stability or the avoidance of local minima?  I haven't proven it, but I wonder whether this new cost function would be convex as opposed to the original cost function.