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!



Solving integral over the space of two lognormally distributed variables
For $c=0$ there is a lengthy closed form,
$$\int\_a^b\int\_0^{d}(x+e)f(x)f(y)dydx=\frac{1}{4} \left(\text{erf}\left(\frac{\ln d}{\sqrt{2} \sigma}\right)+1\right) $$
$$\qquad\qquad\times\left(e \left(\text{erf}\left(\frac{\ln b}{\sqrt{2} \sigma}\right)-\text{erf}\left(\frac{\ln a}{\sqrt{2} \sigma}\right)\right)+e^{\sigma^2/2} \left(\text{erf}\left(\frac{\sigma^2-\ln a}{\sqrt{2} \sigma}\right)-\text{erf}\left(\frac{\sigma^2-\ln b}{\sqrt{2} \sigma}\right)\right)\right)$$


For $c \neq 0$ it looks pretty hopeless...


