TL;DR: the average position is the midpoint between centre of ellipse and empty focus.
Prelude: I'm no mathematician. Expect errors. Get ready for many hand-wavy arguments. Some may be true, but some may be not.
I'm trying to find the average position of a satellite in an elliptical orbit, first analytically, then numerically.
Nomenclature:
y - vertical axis
x - horizontal axis
T - orbital period
r - orbit radius
TrA - true anomaly
MnA - mean anomaly
EcA - eccentric anomaly
n - mean motion
SMa - semi-major axis
SMi - semi-minor axis
ecc - eccentricity
μecc - an auxiliary number, defined as \( \mu_{ecc} = \sqrt{\frac{1-ecc}{1+ecc}} \)
Animation from:
https://commons.wikimedia.org/wiki/File:Elliptic_orbit.gif
Mathematical argument:
First, we agree to that the average position is somewhere on the axis through the two ellipse foci, and thus on the axis through the periapsis and apoapsis.
In mathematical language,
\(<y> = 0\),
where the brackets denote the average/mean.
Then, the average x-position is the mean of all x positions in the orbit:
\(<x> = \frac{1}{T} \int^T_0 x(t) dt\).
As the equation for an elliptical orbit is given in polar coordinates as
\(r(TrA) = \frac{SMa \left(1 - ecc^2\right)}{1+ecc \cos TrA}\),
the x-position is
\(x(TrA) = r(TrA) \cdot \cos TrA\).
Now, we must link the angle TrA to time t. The link will go through the mean anomaly MnA.
The mean anomaly is equal to the time from periapsis times the mean motion \(n = \frac{2 \pi}{T}\), and the mean anomaly is from the Kepler Equation related to the eccentric anomaly, which is itself related to the true anomaly!
\(MnA = n \cdot t = EcA - ecc \cdot \sin EcA \Rightarrow t = \frac{EcA - ecc \cdot \sin EcA}{n}\)
We have that \(\tan \frac{EcA}{2} = \sqrt{ \frac{1-ecc}{1+ecc}} \tan \frac{TrA}{2} \Rightarrow EcA = 2 \arctan\left( \sqrt { \frac{1-ecc}{1+ecc}} \tan \frac{TrA}{2}\right)\), so we then finally have a big equation to relate time t to the true anomaly TrA:
\(t = \frac{2 \arctan\left( \sqrt{\frac{1-ecc}{1+ecc}} \tan \frac{TrA}{2}\right) -ecc \sin \left(2 \arctan\left( \sqrt{\frac{1-ecc}{1+ecc}} \tan \frac{TrA}{2}\right)\right)}{n}\).
Due to the LaTeX on Orbiter-Forum being limited to 400 characters (a limit I'll hit later), I'll have to define the number \(\mu_{ecc} = \sqrt{\frac{1-ecc}{1+ecc}}\).
We must also convert the integral from dt to dTrA, which we can do from Kepler's second equation, which gives (look on Wikipedia if you don't believe me):
\(dt = r^2 \cdot \frac{T}{2\pi SMa SMi} dTrA\).
We now get the original integral to become
\(<x> = \frac{1}{T} \int^T_0 x(t) dt = \frac{1}{T} \int^T_0 r(t) \cdot \cos(TrA(t)) dt = \frac{1}{2 \pi SMa SMi} \int^{2 \pi}_0 r(t)^3 \cdot \cos(TrA(t)) dTrA\)

(uploaded as image, as forum doesn't support such a long equation).
This is where I stopped. I tried to let WolframAlpha and GeoGebra give it a go, but neither could. For all I know, I have done a mistake by now, so the equation may be wrong. But as you understand: finding an analytical answer is not the way for me to go.
Let's go to numerical "solving"!
Numerical argument:
Much simpler for me, especially considering my recent programming efforts with orbital calculations.
Using Python, I do much the same as what I attempted to generalise above: put an object in orbit with some semi-major axis and eccentricity, and then calculate the position every second (or other interval), and find the average.
Fairly easy, and we get the figure below:
Here we see the position of the satellite every 500 seconds, from the perigee at 160 km, with eccentricity 0.800. The three crosses from the left to right are the empty focus (x=PeR-ApR in graph coordinates), centre of the ellipse (x=PeR-SMa in graph coordinates), and the focus with the Earth's centre (x=0 in graph coordinates).
The average position is the green dot.
Wait a second, doesn't that seem to be in the middle between the left focus and the centre of the ellipse?
Calculating the ratio between (average position - centre of ellipse) and (left focus - centre of ellipse), we find that value to be 0.5067.
Reducing the step size to 10 seconds, and checking elliptical orbits with eccentricity 0.1 -- 0.9 and PeA 100 km, 1 000 km, 10 000 km and 100 000 km, we get the following ratios:
As you can see, when increasing steps (i.e. longer period, either by higher perigee or eccentricity), we get values closer and closer to 1/2.
So using this "evidence", we can quite safely say that the average position of an object in a Keplerian orbit is the midpoint between the empty focus and ellipse centre, or precisely \(<x> = PeR - 0.5 \cdot SMa - 0.5 \cdot ApR = SMa \cdot \left(1.5 ecc\right)\) in our coordinate system, which is distance \(1.5 \cdot SMa \cdot ecc\) from the planet centre in the direction of the apoapsis for any elliptical orbit.
Coda:
I assume the average distance between two independent satellites is the distance between their average positions (the point I've found above). As is mentioned by kuddel, this is not true for objects with equal periods, and maybe also not for objects with rational period fractions. But for every irrational period fraction, I assume this mean-point-to-mean-point distance to be the average distance when we let the time pass.
EDIT 25.09.2020: changed from old forum
[math] - [/math]
to
\( - \)
formatting.