Translating the center of a group of objects to the world origin
Use the matrix world.


  is there an alternative to my approach to move the center of a group
  of objects to the origin?


An objects matrix world gives us access to global coordinates.  Hence could look at all the global coordinates of objects, get the mean and subtract it from each.

import bpy
from bpy import context
from mathutils import Vector

obs = context.selected_objects

o = sum((o.matrix_world.translation for o in obs), Vector()) / len(obs) 

for ob in obs:
    ob.matrix_world.translation -= o


On parenting. In the UI we are used to seeing child coordinates such that the parent inverse matrix is identity.  To do this in code would set the matrix local of each child to be in the object space of the parent.  To do this we multiply its world matrix by the inverse of the parents.

Taking up script above, from mean calculation

bpy.ops.object.empty_add(location=o)
mt = context.object
mwi = mt.matrix_world.inverted()

for ob in obs:
    ob.parent = mt
    # alter ob.matrix_local
    ob.matrix_local = mwi @ ob.matrix_world
    # instead of
    # ob.matrix_parent_inverse = mwi
    # which retains matrix_local

mt.location = (0, 0, 0)
bpy.data.objects.remove(mt)


Note the difference here is that the local coordinates are being altered directly to what is the result of setting parent inverse (and keeping local coords same).  

The latter is the reason the objects would return to original location when parent  removed,  parent inverse isn't used in calculations.