Approximate $\int_a^b \frac{1}{\sqrt{2 \pi \sigma^2}}e^{-(x-\mu)^2/2 \sigma^2}\log(1+e^{-x}) \ \ dx $ I am trying to find an approximation to
$$
I = \int_a^b \frac{1}{\sqrt{2 \pi \sigma^2}}e^{-(x-\mu)^2/2 \sigma^2}\log(1+e^{-x}) \ \ dx.
$$
My attempt is as follows:
$$
\begin{align}
I &= \int_a^b \frac{1}{\sqrt{2 \pi \sigma^2}}e^{-(x-\mu)^2/2 \sigma^2} \left( \sum_{i=1}^\infty \frac{e^{-ix}}{i} (-1)^{(i+1)} \right)\ dx\\
&= \sum_{i=1}^\infty \frac{(-1)^{(i+1)}}{i}\int_a^b \frac{1}{\sqrt{2 \pi \sigma^2}}e^{-(x-\mu)^2/2 \sigma^2}  e^{-ix} \ \ dx\\
&= \sum_{i=1}^\infty \frac{(-1)^{(i+1)}}{i} k_i \int_a^b \frac{1}{\sqrt{2 \pi \sigma^2}}e^{-(x-(\mu-i \sigma^2))^2/2 \sigma^2} \ \ dx,\\
\end{align}
$$
where 
$$
k_i = e^{(\mu -i \sigma^2)^2-\mu^2}.
$$
The $k_i$ increases exponentially with increasing $i$ and thus makes the sum divergent. I don't understand why this is happening although this sum should be finite (because I don't see any problem with the original integral).
P.S. I used MacLauren series in approximating natural logarithm.

For the expansion to make sense we need $e^{-x}<1$, so let's assume
$0<a<b$. 
If you do the remaining integral you'll find 
$$\begin{equation*}
I = \sum_{n=1}^\infty \frac{(-1)^{n+1}}{2n} 
\exp\left({-n \mu + \frac{1}{2}n^2\sigma^2}\right)
\left(
\mathrm{erf}\left(\frac{b-\mu+n\sigma^2}{\sqrt{2}\sigma}\right)
- \mathrm{erf}\left(\frac{a-\mu+n\sigma^2}{\sqrt{2}\sigma}\right)
\right),\tag{1}
\end{equation*}$$
where $\mathrm{erf}$ is the error function.
The summand may look badly behaved for large $n$ but the bad behavior is tamed by the asymptotic behavior of $\mathrm{erf}$. 
In fact, for large $n$ the summand goes like 
$$\begin{equation*}
\frac{1}{\sqrt{2\pi}\sigma}
\left[
\exp\left({-\frac{(b-\mu)^2}{2\sigma^2}}\right) \frac{(-1)^n e^{-n b}}{n^2}
- \exp\left({-\frac{(a-\mu)^2}{2\sigma^2}}\right) \frac{(-1)^n e^{-n a}}{n^2} 
\right].\tag{2}
\end{equation*}$$
Note that for $a>0$ the sum $\sum_n (-1)^n e^{-na}/{n^2}$ is absolutely convergent. 
This sum is related to the dilogarithm.
For $(a-\mu+\sigma^2)/\sigma \gg 1$, the integral is well approximated by 
$$\begin{equation*}
I\approx \frac{1}{\sqrt{2\pi}\sigma}
\left[
\exp\left({-\frac{(b-\mu)^2}{2\sigma^2}}\right) \mathrm{Li}_2(-e^{-b})
- \exp\left({-\frac{(a-\mu)^2}{2\sigma^2}}\right) \mathrm{Li}_2(-e^{-a})
\right].\tag{3}
\end{equation*}$$
In general you can cut the sum in (1) off at some appropriate $n$ dependent on your choice of the various parameters.
The higher order terms are exponentially suppressed, so this should work quite well for a good choice of $n$.

Figure 1. Plot of $I(a)$ (solid) and the fit using $(3)$ (dashed) for $\mu=2$, $\sigma=4$, and $b=4$.
