Closed form of the integral of the difference of two Gaussian CDFs? Problem
I'm trying to find the simplest form of the difference of two Guassian CDFs, i.e.

$$
\int_{-\infty}^\infty \left(  \Phi\left(ax+b \right) - \Phi\left(cx+d \right)  \right) dx
$$

for $\Phi(\cdot)$ : standard norm CDF. Unfortunately there is no spoon-feeding solution for this problem in the Wiki page for the List of Integrals of Gaussian functions. For those who may flag this post, this post does not seem to address my question due to some misunderstandings on OP's part. 

Try
But I tried with the indefinite integral of $\int \Phi(ax+b)dx$, which is given by
$$
\int \Phi(ax+b)dx = \frac{1}{b} \left ((a+bx)\Phi(a+bx) + \phi(a+bx)\right) + C
$$
whose source is the above Wiki page. Thus 
$$
\begin{aligned}
\int \left( \Phi(ax+b) - \Phi(cx+d)\right) dx &= \frac{1}{b} \left ((a+bx)\Phi(a+bx) + \phi(a+bx)\right) \\ &- \frac{1}{d} \left ((c+dx)\Phi(c+dx) + \phi(c+dx)\right) + \tilde{C}
\end{aligned}
$$
but in view of definite integral, the tricky part is
$$
 x \left[ \Phi(a+bx) - \Phi(c+dx) \right]  \Large|\normalsize_{-\infty}^\infty
$$
so my struggling does not seem to end up with something.

Speculation
In search of an answer, I kind of found out that if $a=c$, 
$$
\int_{-\infty}^\infty \left(  \Phi\left(ax+b \right) - \Phi\left(ax+d \right)  \right) dx = (b-d)/a
$$
without any proof. The way I found out this is via R. The following code illustrates how I speculated it(for some arbitrary mu, sd). 
# mu, sd can be arbitrary numbers

a = 0.5; c = 0.5
b = 2; d = -4

lb = -1000; ub = 1000 # suff. large lower/upper bounds for test integration

f = function(theta){
  CDF1 = pnorm(a*theta + b, mean=mu, sd=sd)
  CDF2 = pnorm(c*theta + d, mean=mu, sd=sd)
  return(CDF1-CDF2)
}

f = Vectorize(f)
(val = integrate(f, lower=lb, upper=ub)$value)
# the result is (b+d)/a, for any b,d if a=c

But for the cases in which $a$, $c$ are different, I currently have no idea. 

Question
I would like to find the general form for the cases in which $a \neq c$ with a proof. I may be missing a trivial part, but I would like to ask for a help. Any helping hands will be greeted. Thanks.
 A: A geometrical intuition
A geometrical intuition to accompany Whuber's answer is the following:

It relates to answers to other questions here and here.
The mean of the variable relates to the area of the gray striped surfaces, which can be computed in two directions

*

*The vertical stripes: as an integral of the quantile function with the quantiles $dp$ as integrand
$$E[X] = \int_0^1 Q_X(p) dp$$


*The horizontal stripes: as an integral with the CDF and the variable $dz$ as integrand
$$E[X] = -\int_{-\infty}^0 F_X(z) dz + \int_{0}^{\infty} 1-F_X(z) dz$$
Then the difference $E[X]-E[Y]$ is
$$\begin{array}{}
E[X]-E[Y] &=& \int_{0}^1 Q_X(p)-Q_Y(p)  dp \\
\end{array}$$
which relates to the difference of the quantile functions, and the blue and red areas marked by the vertical stripes.
But we can just as well compute this in the horizontal direction by taking the difference between the CDF functions.
$$\begin{array}{}
E[X]-E[Y] &=& -\int_{-\infty}^0 F_X(z)-F_Y(z) dz + \int_{0}^{\infty} -F_X(z)+F_Y(z) dz \\
&=& - \int_{-\infty}^\infty F_X(z)-F_Y(z) dz 
\end{array}$$
Applied to Gaussian CDF
With the general equation
$$\begin{array}{}
E[X]-E[Y] &=& - \int_{-\infty}^\infty F_X(z)-F_Y(z) dz 
\end{array}$$
you will get for the CDF's of the Gaussian distribution $\Phi\left(ax+b \right)$ which has the mean $\mu(a,b) = -b/a$
$$\int_{-\infty}^\infty \left(  \Phi\left(ax+b \right) - \Phi\left(cx+d \right)  \right) dx = -[\mu(a,b) - \mu(c,d)] = -b/a - c/d$$
