I agree with @Kermittfrog's comment, that this only works if you do not impose any budget constraint (in the sense that your weights sum up to one). Other than that, I am sorry that I can not precisely answer your question where it was first derived (tbh: I am not even sure if it was explicity derived at all somewhere because it simply follows from the very definition of a target-return constraint optimziation problem). In brief:
Let $\mu$ be a $n$-dimensional vector that contains the single-asset returns, and $\Sigma$ the corresponding $n \times n$ covariance matrix. Moreover, let $\mathbf{w}$ be a $n$-dimensional vector that assigns the portfolio weights to the $n$ individual assets. (Here we do not impose summation to unity, see comment above!) Furthermore, express by $r$ the target return an investor expects to achieve on her portfolio (thus, $r$ is a scalar). In other words, the optimal solution must in addition satisfy $r = \mu^{T} \mathbf{w}$, where superscript $T$ denotes a vector/matrix transpose.
Now, we seek to minimize (for mathematical convenience: half of) the portfolio variance, which is given by $\mathbf{w}^{T}\Sigma\mathbf{w}$. And we will minimize it under the target return restriction presented above. As we have formulated it as an equality constraint, we can use the Lagrangian method to solve this problem. Defining the Lagrangian multiplier $\lambda$, we can write the first-order condition of the optimization problem (with respect to the weight vector) as follows:
$$ \Sigma \mathbf{w} - \lambda \mu = \mathbf{0} $$
Note that $\mathbf{0}$ describes the all-zero vector of dimension $n$. Solving this for $\mathbf{w}$ is easy, and defining by $\Sigma^{-1}$ the inverse of our covariance matrix, we obtain that:
$$ \mathbf{w} = \Sigma^{-1} \lambda \mu$$
Moreover, recall our optimization constraint, which will help us to find $\lambda$; it reads $ r = \mu^{T} \mathbf{w} $. If you plug in the solution for $\mathbf{w}$ and rearrange, you will get that $\lambda = \dfrac{r}{\mu^{T} \Sigma^{-1} \mu}$.
Clearly, you see that this Lagrangian "slack" parameter depends on your choice of target return, $r$. Plugging this Lambda back into our weight solution, you obtain the proposed result:
$$ \mathbf{w}(r) = \dfrac{r \Sigma^{-1} \mu}{\mu^{T} \Sigma^{-1} \mu} $$.
Finally, in a similar fashion proposed in the answer to your other post (Closed-form analytical solution for the variance of the minimum-variance portfolio?), once you have your ($r$-dependent) weights, the ($r$-dependent) mean and variance of the portfolio follow therefrom.
