Answers:
@snipperrabbit!!, Blkbird: You are probably correct: No special bonus system for vertex cities, because overall they are already statistically balanced and currently are neither "hyper-valuable" nor "hyper-bad" places to settle down.
@Atrebates: Thanks for the links; interesting.
@hunterai: Thank you.
Bad news: 
I have finally managed to programmatically create the
tile neighbourhood graph for all n. I also implemented the
A* pathfinding algorithm on top of it. I even reduced the tile area deviations to less than +-10%, but then I faced the following:
While studying the shortest tile paths based on the A* algorithm, I saw some artifacts in the form of
deviations from the great circle path (which is the mathematically shortest path) between the start and end tile.
I first thought it was some coding error, but it turned out to be an
intrinsic fact of any tesselation of the sphere (also the pentagon solution): Since your tiles only have a limited number of "preferred directions" (6 in case of hexagonal tiles), the
map requires your unit to go zig-zag, if you want to travel along any other direction. This results in the
direct great circle tile-path having more tiles than another non-direct path. I find this very unattractive and unintuitive. Below is an image, which shows one of the worst cases (look inside the red ellipse and compare the tiles of the zig-zag great circle arc with the outer non-zig-zag shortest tile path):
Of course, as shown previously, you also have this zig-zag problem in 2D, but it is much more obvious on a sphere:
What do you think about this "new" problem? Any ideas how to solve it or at least reduce it?