Expected value of applying the sigmoid function to a normal distribution

Apart from the the MacLaurin approximation, the usual way to compute that integral in Statistics is to approximate the sigmoid with a probit function. More specifically $\mathrm{sigm}(a) \approx \Phi(\lambda a)$ with $\lambda^2=\pi/8$. Then the result would be: $$\int \mathrm{sigm}(x) \, N(x \mid \mu,\sigma^2) \, dx \approx \int \Phi(\lambda x) \, N(x \mid \mu,\sigma^2) \, dx = \Phi\left(\frac{\mu}{\sqrt{\lambda^{-2} + \sigma^2}}\right).$$

Since I do not have enough reputation to comment, I'll instead add a new answer. @korkinof's answer is almost correct. The final integral evaluates to the following: \begin{equation} \int_x sigmoid(x) \mathcal{N}(x; \mu, \sigma^2) \approx \int_x \Phi(\lambda x) \mathcal{N}(x; \mu, \sigma^2) = \Phi\left(\frac{\lambda \mu}{\sqrt{1 + \lambda^2 \sigma^2}}\right). \end{equation} I verified my answer through simulation.

I doubt that there's a closed-form solution. However, here's a series in powers of $\sigma$:

$$ \left( {{\rm e}^{-{\mu}}}+1 \right) ^{-1}+{\frac { \left( { {\rm e}^{-{\mu}}}-1 \right) {{\rm e}^{-{\mu}}}}{2\, \left( {{\rm e} ^{-{\mu}}}+1 \right) ^{3}}}{{\sigma}}^{2}+{\frac { \left( { {\rm e}^{-3\,{\mu}}}-11\,{{\rm e}^{-2\,{\mu}}}+11\,{{\rm e}^{-{ \mu}}}-1 \right) {{\rm e}^{-{\mu}}}}{8\, \left( {{\rm e}^{-{\mu} }}+1 \right) ^{5}}}{{\sigma}}^{4}+{\frac {{{\rm e}^{-{\mu} }} \left( {{\rm e}^{-5\,{\mu}}}-57\,{{\rm e}^{-4\,{\mu}}}+302\,{ {\rm e}^{-3\,{\mu}}}-302\,{{\rm e}^{-2\,{\mu}}}+57\,{{\rm e}^{-{ \mu}}}-1 \right) }{48\, \left( {{\rm e}^{-{\mu}}}+1 \right) ^{7}}}{{ \sigma}}^{6}+{\frac {{{\rm e}^{-{\mu}}} \left( {{\rm e}^{-7\,{\mu}}}-247\,{{\rm e}^{-6\,{\mu}}}+4293\,{ {\rm e}^{-5\,{\mu}}}-15619\,{{\rm e}^{-4\,{\mu}}}+15619\,{ {\rm e}^{-3\,{\mu}}}-4293\,{{\rm e}^{-2\,{\mu}}}+247\,{{\rm e}^{ -{\mu}}}-1 \right) }{384\, \left( {{\rm e}^{-{\mu}}}+1 \right) ^{9}}} {{\sigma}}^{8}+O \left( {{\sigma}}^{10} \right) $$

EDIT: To obtain this, first do the change of variables $x = \mu + \sigma t$. The integral becomes $$ \frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty \dfrac{e^{-t^2/2}}{1 + e^{-\mu - \sigma t}}\ dt $$ Now take the Maclaurin series $$\frac{1}{1+e^{-\mu - \sigma t}} = \frac{1}{1+e^{-\mu}} + \frac{e^{-\mu} \sigma t}{(1+e^{-\mu})^2} + \frac{e^{-\mu} ( e^{-\mu} - 1) \sigma^2 t^2}{(1+e^{-\mu})^3} + \ldots$$ and integrate term by term.