My try: Split the sum into three parts
$$\sum_{k=1}^n \left(\frac{k}{n}\right)^k = \sum_{1\leq k\leq K} \left(\frac{k}{n}\right)^k + \sum_{K+1 \leq k < n-n^\epsilon} e^{k\log\left(\frac{k}{n}\right)} + \sum_{n-n^\epsilon \leq k \leq n} e^{k\log\left(\frac{k}{n}\right)} $$
for some integer $K$ which defines the order and some small $\epsilon >0$ (say $\epsilon=1/2$). It is easy to see that $k\log\left(\frac{k}{n}\right)$ has a unique minimum at $k=\frac{n}{e}$ somewhere in the range of the middle term for large $n$. Therefore we evaluate the boundary terms of the middle term for some estimate $$k=K+1: \quad \left(\frac{K+1}{n}\right)^{K+1} \\
k=n-n^\epsilon: \quad e^{n(1-n^{\epsilon-1})\log(1-n^{\epsilon-1})} \leq e^{-n^\epsilon + n^{2\epsilon -1}} \, .$$
For fixed $K$ and sufficiently large $n$ the right boundary obviously vanishes exponentially (the optimal $\epsilon$ is $1-\frac{\log 2}{\log n}$ so that $n^\epsilon=n/2$) and the largest value in that range is that for $k=K+1$. Hence the middle term is of order ${\cal O}(n^{-K})$.
For the last term substitute $k\rightarrow k-n$ so that it becomes
$$\sum_{0\leq k \leq n^\epsilon} e^{-k +\left[(n-k)\log\left(1-\frac{k}{n}\right) + k\right]} \, .$$
It remains to estimate the square bracket $$(n-k)\log\left(1-\frac{k}{n}\right) + k = -(n-k) \sum_{m=1}^\infty \frac{1}{m}\left(\frac{k}{n}\right)^m + k \\
= \frac{k^2}{n} - (n-k) \sum_{m=2}^\infty \frac{1}{m}\left(\frac{k}{n}\right)^m \\
= \frac{k^2}{2n} + \sum_{m=2}^\infty \frac{k}{m(m+1)} \left(\frac{k}{n}\right)^m = \sum_{m=1}^\infty \frac{k}{m(m+1)} \left(\frac{k}{n}\right)^m$$
which vanishes vor large $n$. For an order $K$ approximation we can thus write
$$\sum_{0\leq k \leq n^\epsilon} e^{-k +\left[(n-k)\log\left(1-\frac{k}{n}\right) + k\right]} \\
= \sum_{0\leq k \leq n^\epsilon} e^{-k} \left\{ 1 + \sum_{l=1}^\infty \frac{1}{l!} \sum_{m_1=1}^\infty \cdots \sum_{m_l=1}^\infty \frac{k^{l+m_1+\dots+m_l}}{m_1(m_1+1)\cdots m_l(m_l+1)} \frac{1}{n^{m_1+\dots+m_l}} \right\} \\
= \sum_{0\leq k \leq n^\epsilon} e^{-k} \left\{ 1 + \sum_{p=1}^\infty \frac{k^p}{n^p} \sum_{l=1}^p \frac{k^{l}}{l!} \substack{ \sum_{m_1=1}^\infty \cdots \sum_{m_l=1}^\infty \\ m_1+\dots+m_l \, \stackrel{!}{=} \, p }\frac{1}{m_1(m_1+1)\cdots m_l(m_l+1)} \right\} \, .$$
When evaluating the moments
$$ \sum_{0\leq k \leq n^\epsilon} e^{-k} \, k^{p+l} $$
for $p=0,1,2,...,K-1$, the range of summation can be extended up to infinity, because that only introduces an exponentially suppressed error term ${\cal O}\left(n^{(p+l)\epsilon} \, e^{-n^\epsilon}\right)$.
Collecting terms, it is seen that
$$\sum_{k=1}^n \left(\frac{k}{n}\right)^k = a_0 + \sum_{k=1}^{K-1} \frac{k^k + a_k}{n^k} + {\cal O}\left(n^{-K}\right) $$
where
$$a_0 = \frac{e}{e-1} \\
a_k = \sum_{l=1}^k \frac{\sum_{q=0}^\infty q^{k+l} \, e^{-q}}{l!} \substack{ \sum_{m_1=1}^\infty \cdots \sum_{m_l=1}^\infty \\ m_1+\dots+m_l \, \stackrel{!}{=} \, k }\frac{1}{m_1(m_1+1)\cdots m_l(m_l+1)} \, .$$
For $k\geq 2$ the $a_k$ are extremely close to $k^k$, that is less than $0.04\%$ relative error, so that the total coefficient for $k\geq 2$ is in good approximation $2k^k$.
One term beyond leading order we have for $K=2$
$$\sum_{k=1}^n \left(\frac{k}{n}\right)^k = \frac{e}{e-1} + \frac{1+\frac{e(e+1)}{2(e-1)^3}}{n} + {\cal O}(n^{-2}) \, .$$
Increasing the order $K$ also shifts the range of validity to higher $n$, i.e. it is an asymptotic series. The zero, first and fifth order approximations are shown below. The fifth order is visually not distinguishable from the approximation where $a_k=k^k$ has been used for $k\geq 1$. 