Derivation of Yule-Walker Equation for LPC
Yes, your derivation is correct.
There is just a minor inconsistency in your definition of the auto-correlation function $r(k)$. If you choose to define
$$r(k) = E[s(n-k)s^*(n)]$$
then $E[s(n-k)s^*(n-i)]$ equals $r(k-i)$ or, equivalently, $r^*(i-k)$, but not $r(i-k)$. So the last equation should read
$$\sum_{i=1}^pa^*_i r^*(i-k) = -r(k)$$
You take the derivative w.r.t. to $a_k$ because these are the coefficients that you can choose as to minimize the error. They are the unknowns in your problem.
The fact that $e^*(n)$ can be viewed as a constant when taking the derivative w.r.t. $a_k$ has to do with how computing the gradient w.r.t. complex variables works. It is equivalent to taking the derivatives w.r.t. to the real and imaginary part of the coefficients separately, but the complex technique is more elegant and easier (at least as soon as you understand it). I've explained this technique in this answer.
The complex conjugates just pop up because you consider $|e(n)|^2$ in your optimization objective: $|e(n)|^2=e(n)e^*(n)$. You wouldn't need them at all if you only considered real-valued coefficients $a_k$ and real-valued signals $s(n)$.