Design of FIR Filters with Arbitrary Magnitude and Phase Responses
For the sake of clarity, let me point out that you do not approximate $(1)$ by $(2)$, but the other way around: you approximate a desired frequency response, i.e., your specification, by the filter's frequency response, and this is done by choosing the coefficients in the vector $\mathbf{h}$ such that some error measure (e.g., least squares) is minimized.

Since the frequency response $(1)$ is linear in $\mathbf{h}$, the problem is a linear least squares problem, which can be solved by solving a system of linear equations. This is true for the continuous error measure $(3)$ as well as for the discrete error measure

$$E(\mathbf{h})=\sum_{k=1}^K|H(e^{j\omega_k})-D(e^{j\omega_k})|^2\tag{1}$$

where $K$ is the number of discrete frequency points $\omega_k$. In general we require $K\gg N$, where $N$ is the number of FIR filter coefficients in the vector $\mathbf{h}$.

In the general case it is more convenient to use the discrete error measure $(1)$. There are only a few special cases in which the integral of the continuous error measure (Eq. $(3)$ in your question) can be solved analytically, otherwise you have to resort to numerical approximations anyway. And even if an analytical solution is available, the difference when minimizing the discrete error measure is negligible if the number of frequency points $K$ is sufficiently large.

Minimizing the discrete squared error $(1)$ is equivalent to solving the following overdetermined linear system of equations in a least squares sense:

$$H(e^{j\omega_k})\stackrel{!}{=}D(e^{j\omega_k}),\qquad k=1,2,\ldots,K\tag{2}$$

Formulation $(2)$ is convenient because many software packages such as Matlab have syntactically elegant and computationally efficient methods to solve such systems.

I've implemented a complex FIR filter design based on solving $(2)$ in the (extremely simple) Matlab/Octave function cfirls.m.