Given a series of points on the Earth's surface: which is the geodesic (great-circle) that minimizes the total cross distance between it and the reference points?

There may be an analytical solution but it seems complicated enough, so I have taken a numerical approach. The geodesic formulas are taken from movable-type.

Attached is a Jupiter Notebook that provides a solution based on a double iteration:

• it starts navigating through a random selection of great-circles, looking for those which minimize the total sum,
• then refines the result with a gradient search with the same objective.

(In my trials the core of the job is done by the random walk, but this clearly depends on the number and conditions of the iterations.)

The script provides visual checks on the convergence of the random approach:

Finally, it also generates a representation of the added cross distance from the set of points to "all" great-circles:

It allows matplotlib visualizations:

And (manual) zooming_in: