Approximation and IIR Filters
I think it's instructive to try to understand how this error function is derived and why it makes sense. First, the desired impulse response $h_d[n]$ is given implicitly by the input-output relation of the desired system:
$$y[n]=\sum_{k=1}^{5}x[n-k]=\sum_{k=-\infty}^{\infty}h_d[k]x[n-k]\tag{1}$$
The given error function is the error function minimized by Prony's method for the design of IIR filters. We try to approximate a given transfer function $H_d(z)$ by an IIR filter $H(z)=B(z)/A(z)$:
$$B(z)\stackrel{!}{=}H_d(z)A(z)\tag{2}$$
For the given example we get (with $B(z)=1$)
$$1\stackrel{!}{=}\left(h_d[0]+h_d[1]z^{-1}+\ldots+h_d[N-1]z^{-(N-1)}\right)\left(1+a_1z^{-1}+z_2z^{-2}\right)\tag{3}$$
For equation $(3)$ to be satisfied, we need all coefficients associated with negative powers of $z$ to vanish:
$$\begin{align}h_d[1]+h_d[0]a_1&amp;\stackrel{!}{=}0\\
h_d[2]+h_d[1]a_1+h_d[0]a_2&amp;\stackrel{!}{=}0\\
h_d[3]+h_d[2]a_1+h_d[1]a_2&amp;\stackrel{!}{=}0\\\vdots\end{align}\tag{4}$$
In practice, we can solve $(4)$ in an approximate way by minimizing the sum of the squares of the left-hand side, leading to the given error function
$$E=\sum_{n=-\infty}^{\infty}\big(h_d[n]+h_d[n-1]a_1+h_d[n-2]a_2\big)^2\tag{5}$$
Finally, the optimal coefficients $a_1$ and $a_2$ are obtained by taking the derivative of $(5)$ w.r.t. $a_1$ and $a_2$ and equating them to zero. This results in two linear equations with two unknowns.