I want to estimate the value of this sequence for large $n$ (with a reasonable lower bound and upper bound).
That is, can we find a function $f(n)$ such that $$ \frac{1^n + 2^{n-1} + 3^{n-2} + \cdots + n^1}{f(n)} \rightarrow 1?$$
I tried to approximate $ (n - x)^x $ as a Gaussian, look at how close it is visually, but the manipulation escapes me.
Here is the case where $n=20$:


Addendum: This post might be relevant.
If you let $m_n>1$ be the unique solution to $m_n = \frac{n}{1+\ln(m_n)}$ and $S_n=\sum_{k=1}^n k^{n-k}$ (I have shifted n by one as in your plots), then $$ \lim_{n\rightarrow +\infty} \frac{1}{S_n} m_n^{n-m_n+1} \sqrt{\frac{2\pi}{m_n+n}} =1.$$ This is what you would suspect using the method of Laplace, writing $x^{n-x}=e^{f_n(x)}$ approximating $f_n$ by a parabola at the maximum (which is for $x=m_n$) and calculating the resulting gauss integral. I didn't check rigorously if the third order correction to the Laplace formula is negligible but numerically the above seems to hold.
EDIT: Let me outline a sketch of a proof (filling in details would probably amount to writing about 3 pages, dense with formulae). It follows more or less known lines in the case of Stirlings formula.
Step 1: Show that $\sum_{k=1}^n k^{n-k}$ is equivalent to the integral $\int_1^n e^{f_n(x)}dx$ with $f_n(x)=(n-x)\ln(x)$ as $n\rightarrow +\infty$. This is not so difficult.
Step 2: We write $m=m_n$. Substituting $x=m+ t \frac{m}{\sqrt{m+n}}$ and using $\ln m = \frac{n-m}{m}$ one obtains after some algebra: \begin{align} f_n(x) & = m \left( \frac{n-m}{m} - \frac{t}{\sqrt{n+m}}\right) \left( \frac{n-m}{m} + \ln \left(1+ \frac{t}{\sqrt{n+m}} \right)\right)\\ & = m\left(\frac{n-m}{m}\right)^2 -\frac{t^2}{2} + O\left(\frac{t^3}{\sqrt{n}}\right)\end{align} and then $ \int_1^n e^{f_n(x)} dx = \frac{m_n^{n-m_n+1}}{\sqrt{n+m}} \int_{-l_n}^{u_n} \exp \left( -\frac{t^2}{2} + O \left( \frac{t^3}{\sqrt{n}}\right)\right)\; dt.$ Here, $l_n,u_n\rightarrow +\infty$ as $n$ goes to infinity and pointwise the integrand goes to $e^{-t^2/2}$.
Step 3: Show that the integrand is bounded uniformly in $n$ by an integrable function and apply Lebesgue dominated convergence to conclude the proof.
For Step 3, note that a Gaussian function would not work as a dominating function due to the logarithmic tail in $f_n$. You need to construct the dominating function in a more clever way. This is where calculations become more ugly and I leave it aside.