I am interested in the closed form solution to the following problem:

$\int_a^b\int_0^{cx+d}(x+e)f(x)f(y)dydx$, where $f(.)$ is the pdf of the lognormal distribution with mean $0$ and variance $\sigma^2$. $a,b,c,d,e$ are constants. Say $F(.)$ is the corresponding cdf.

I got as far as $\int_a^b\int_0^{cx+d}(x+e)f(x)f(y)dydx=\int_a^bxf(x)F(cx+d)dx+e\int_a^bf(x)F(cx+d)dx$.

But I am not sure if generating a closed form solution to the above integral is possible. Given this, I am happy to approximate the lognormal cdf $F(x)$ to $G(x)=\frac{x^\alpha}{1+x^\alpha}$, with pdf $g(x)=\frac{\alpha x^{\alpha-1}}{(1+x^\alpha)^2}$, where $\alpha>1$ and $\alpha$ is smaller the larger the $\sigma$. The idea for this approximation came from https://www.jstor.org/stable/3621676?seq=1#metadata_info_tab_contents (Mihálykó, Csaba, and Tibor Blickle. “84.48 On the Approximation of the Lognormal Distribution.” The Mathematical Gazette, vol. 84, no. 500, 2000, pp. 308–311.).

Then my problem becomes one of finding closed form solutions to:

$\int_a^bxg(x)G(cx+d)dx=\int_a^bx\frac{\alpha x^{\alpha-1}}{(1+x^\alpha)^2}\frac{(cx+d)^\alpha}{1+(cx+d)^\alpha}dx$

and

$\int_a^bg(x)G(cx+d)dx=\int_a^b\frac{\alpha x^{\alpha-1}}{(1+x^\alpha)^2}\frac{(cx+d)^\alpha}{1+(cx+d)^\alpha}dx$.

But I am still lost as to how to exactly do this. Detailed steps would be appreciated, thanks!