We can justify the solution to this problem by symmetry. Being i.i.d., the observed $X_i$ are interchangeable, so their coefficients ought to be also. So it is no surprise they are equally weighted. Note that this symmetry argument only works in the case where your observations are of equal variance, and in econometrics you will study situations where this is not the case. It would be useful to learn some machinery that can handle the present case more formally, and can be adapted to handle such variations.
Problems in which you seek to optimize $f(a_1, a_2, \dots, a_n)$ subject to a constraint that $g(a_1, a_2, \dots, a_n) = c$ are well-suited to solution by Lagrange multipliers.
In particular you seek to minimize $f(a_1, a_2, a_3) = \sigma^2 a_1^2 + \sigma^2 a_2^2 + \sigma^2 a_3^2$ subject to $g(a_1, a_2, a_n) = a_1 + a_2 + a_3 = 1$ (so that $c=1$).
Defining $\Lambda = f + \lambda (g - c)$ we see that, since $g(a_1, a_2, a_3)-c = 0$, minimizing $f$ is equivalent to minimizing $\Lambda(a_1, a_2, a_3, \lambda)$. We find the partial derivatives and set them equal to zero:
$$\Lambda(a_1, a_2, a_3, \lambda) = \sigma^2 a_1^2 + \sigma^2 a_2^2 + \sigma^2 a_3^2 + \lambda (a_1 + a_2 + a_3 - 1)$$
$$\frac{\partial \Lambda}{\partial a_1} = 2 \sigma^2 a_1 + \lambda = 0$$
$$\frac{\partial \Lambda}{\partial a_2} = 2 \sigma^2 a_2 + \lambda = 0$$
$$\frac{\partial \Lambda}{\partial a_3} = 2 \sigma^2 a_3 + \lambda = 0$$
$$\frac{\partial \Lambda}{\partial \lambda} = a_1 + a_2 + a_3 - 1 = 0$$
From the partial derivatives with respect to each $a_i$ we can see $a_1 = a_2 = a_3 = -\frac{\lambda}{2 \sigma^2}$, and since each must be equal, the final partial derivative (which is equivalent to the original constraint) ensures that $a_1 = a_2 = a_3 = \frac{1}{3}$. So we must weight the observations $x_1$, $x_2$ and $x_3$ equally when taking their mean.
A little bit of thought shows this generalizes to $n$ observations; for $i=1, 2, \dots, n$ the partial derivative with respect to $a_i$ is simply:
$$\frac{\partial \Lambda}{\partial a_i} = 2 \sigma^2 a_i + \lambda = 0$$
Hence $a_i = -\frac{\lambda}{2 \sigma^2}$ for each $i$, so the weights are equal. Moreover, the final partial derivative would be:
$$\frac{\partial \Lambda}{\partial \lambda} =\sum_{i=1}^{n} a_i - 1 = 0$$
Since the weights total to one, each must be $a_i = \frac{1}{n}$.
Suppose instead that the variances of your observations were unequal, $\text{Var}(X_i)=\sigma^2_i$. Now we minimize $f(a_1, a_2, a_3) = \sigma_1^2 a_1^2 + \sigma_2^2 a_2^2 + \sigma_3^2 a_3^2$ subject to $g(a_1, a_2, a_n) = a_1 + a_2 + a_3 = 1$. This time,
$$\frac{\partial \Lambda}{\partial a_i} = 2 \sigma_i^2 a_i + \lambda = 0$$
$$\frac{\partial \Lambda}{\partial \lambda} =\sum_{i=1}^{n} a_i - 1 = 0$$
Since each $a_i = -\frac{\lambda}{2 \sigma_i^2}$ we can see $a_i \propto \frac{1}{\sigma^2_i}$. So observations with a higher variance are given a lower weighting - this makes sense, as they are in some sense less trustworthy. We give larger weightings to observations with lower variances, since they convey more information about the true value of $\mu$. To ensure weights sum to one, the formula for the weights is:
$$a_i = \frac{1}{\sigma^2_i} \left( \sum_{j=1}^n \frac{1}{\sigma^2_j} \right)^{-1}$$
A confession: at this point I have simply claimed that the weights are "optimal", without checking whether I minimized or maximized the variance! Actually, it is possible to see that the solutions obtained so far are minima in various ways, for example geometrically. But since you are studying econometrics, the method I expect you will encounter is the bordered Hessian, which in the general case, with one constraint, appears as:
$$H^B(f,g) = \begin{vmatrix}
0 & \dfrac{\partial g}{\partial a_1} & \dfrac{\partial g}{\partial a_2} & \cdots & \dfrac{\partial g}{\partial a_n} \\[2.2ex]
\dfrac{\partial g}{\partial a_1} & \dfrac{\partial^2 f}{\partial a_1^2} & \dfrac{\partial^2 f}{\partial a_1\,\partial a_2} & \cdots & \dfrac{\partial^2 f}{\partial a_1\,\partial a_n} \\[2.2ex]
\dfrac{\partial g}{\partial a_2} & \dfrac{\partial^2 f}{\partial a_2\,\partial a_1} & \dfrac{\partial^2 f}{\partial a_2^2} & \cdots & \dfrac{\partial^2 f}{\partial a_2\,\partial a_n} \\[2.2ex]
\vdots & \vdots & \vdots & \ddots & \vdots \\[2.2ex]
\dfrac{\partial g}{\partial a_n} & \dfrac{\partial^2 f}{\partial a_n\,\partial a_1} & \dfrac{\partial^2 f}{\partial a_n\,\partial a_2} & \cdots & \dfrac{\partial^2 f}{\partial a_n^2}
\end{vmatrix}$$
Note that I am following Binmore and Davies, Calculus: Concepts and Methods, in calling the determinant the Hessian, rather than the matrix itself (as Wikipedia does). Observe that the upper left 1-by-1 submatrix of our ($n+1$)-by-($n+1$) matrix is just zero - corresponding to $\frac{\partial^2 \Lambda}{\partial \lambda^2} = 0 $. If in even fuller generality there had been $m$ constraints to apply, then the upper left $m$-by-$m$ submatrix of our ($n+m$)-by-($n+m$) matrix would contain only zeroes. The signs of the leading principal minors (determinants of upper left submatrices) are often sufficient determine whether the solution is a minimum or maximum. The only informative ones consist of the first $2m+1$ rows and columns, or more. If they all have the same sign, and this sign is that of $(-1)^m$, then we have a minimum. If they alternate in signs, with the smallest having sign $(-1)^{m+1}$ (or alternatively, the largest - i.e. the whole bordered Hessian - has sign $(-1)^n$) then we have a maximum. This is covered in Wikipedia but for more detail see e.g. these lecture notes (or many others which can be found by an internet search).
I will return our attention to your case where $m=1$ and $n=3$, and the bordered Hessian is:
$$H^B(f,g) = \begin{vmatrix}
0 & 1 & 1 & 1 \\[2.2ex]
1 & 2\sigma^2 & 0 & 0 \\[2.2ex]
1 & 0 & 2\sigma^2 & 0 \\[2.2ex]
1 & 0 & 0 & 2\sigma^2
\end{vmatrix}$$
In this case $2m+1=3$ so we only need to consider the principal leading minors of sizes 3-by-3 and 4-by-4. The latter is just $H^B_4 = H^B$ itself, and the former is:
$$H^B_3 = \begin{vmatrix}
0 & 1 & 1 \\[2.2ex]
1 & 2\sigma^2 & 0 \\[2.2ex]
1 & 0 & 2\sigma^2
\end{vmatrix}$$
This determinant can easily be reduced to that for a lower-triangular matrix by subtracting $\frac{1}{2 \sigma^2}$ times row 2 from row 1, then subtracting the same multiple of row 3 from row 1:
$$H^B_3 = \begin{vmatrix}
-\frac{1}{\sigma^2} & 0 & 0 \\[2.2ex]
1 & 2\sigma^2 & 0 \\[2.2ex]
1 & 0 & 2\sigma^2
\end{vmatrix} = -4 \sigma^2 < 0$$
By the same procedure, we can show the deteminant $H^B_4$ is negative. The relevant leading principal minors had the same sign as each other and also of $(-1)^m = (-1)^1 = -1$ and so we can see our solution is indeed the minimum variance. There is no difficulty applying the same lower-triangularization technique when there are $n$ observations instead of three, or even to the case with unequal variances: all the appropriate minors are clearly negative after triangularization, since the diagonal entries are positive with the exception of the top-left entry, which will be negative as it is zero subtract various positive multiples of one.
