Generate random ellipsoids
Cobbled together a test script, that has parts that answer both of your q's.  Scale a sphere on x, y, z to make ellipsoids. (set sphere_only=False in script)  

The 'within_touch' method looks at two spheres and if they are closer than the sum of their two radii, they are overlapping.  The distance between the spheres is the length of the vector created by subtracting ones location from the other. 

If a pair of spheres are found to be overlapping the method returns the required vector to move to make them touch (with a bit of tol)

While testing added an object color, the bluer the harder to place.

import bpy
context = bpy.context
from random import uniform
from mathutils import Vector

count = 200

# random scale 
min_scale = 0.001
max_scale = 1
# scale_factor .. scale down to help fit
scale_factor = 0.1
domain = 1
only_spheres = True

tries_per_sphere = 256

use_context_object = False
obj = context.object

def random_vector(a, b):
    return Vector([uniform(a, b) for c in "xyz"])

def checkbounds(sphere, fix):
    #return True
    loc = sphere.location + fix
    #dm = max(sphere.dimensions) / 2
    dm = max(sphere.scale)
    _min = -domain + dm
    _max =  domain - dm
    for i, v in enumerate(loc):
        if v &lt; _min:
            loc[i] = _min
        if v &gt; _max:
            loc[i] = _max
    return loc

def inbounds(sphere, loc):
    #return True
    dm = max(sphere.scale)
    return -domain + dm &lt; min(loc) and  max(loc) &lt; domain - dm

def random_sphere(sphere):
    # random scale and rotation for sphere based on settings
    if only_spheres:
        scale = uniform(min_scale , max_scale) * Vector((1, 1, 1))
    else:
        scale = random_vector(min_scale, max_scale)
    scale = Vector([min(s, domain) for s in scale])
    dom = domain - max(scale)    
    sphere.scale = scale

    sphere.location = random_vector(-dom, dom) 

scene = context.scene

mat = bpy.data.materials.get("randobjcol")
if not mat:
    mat = bpy.data.materials.new("randobjcol")
    mat.use_object_color = True

boundbox = bpy.data.objects.get("BoundBOX")
if not boundbox:
    bpy.ops.mesh.primitive_cube_add()
    boundbox = context.scene.objects.active
    boundbox.name = "BoundBOX"
    boundbox.draw_type = 'WIRE'
    boundbox.hide_select = True
    boundbox.location = (0, 0, 0)

if not scene.objects.get("BoundBOX"):
    scene.objects.link(boundbox)

boundbox.scale = domain * Vector((1, 1, 1))

spheres = []
if use_context_object and obj:
    sphere = obj
else:
    bpy.ops.mesh.primitive_uv_sphere_add()
    sphere = context.scene.objects.active
    bpy.ops.object.shade_smooth()
    scene.objects.unlink(sphere)

sphere.active_material = mat

for i in range(count):
    spheres.append(sphere)
    random_sphere(sphere)    
    sphere = sphere.copy()
    #context.scene.objects.link(sphere)

# update the spheres dimensions    
#context.scene.update()        

tries = 0
t_count = 0
r_count = 0
moves = 0
resets = 0
pfix = Vector()

def within_touch(s1, s2):
    def r(s):
        return max(s.scale)
        #return max(s.dimensions) / 2
        #return sum(s.dimensions) / 3 #  average.

    d = (s2.location - s1.location)

    r1, r2 = r(s1), r(s2)

    radsum = r(s1) + r(s2)
    if d.length &lt;= 0.0001:
        print("SAME")
        return Vector((r1, r1, r1))
        # same position    
    elif d.length &lt; radsum:
        # return a vector to move away
        d.length = radsum - d.length
        return  -d
    return Vector()

sphere = spheres.pop()
arranged_spheres = [sphere]
sphere = spheres.pop()

while sphere and tries &lt; tries_per_sphere * count:
    if t_count &gt; 5 and abs(min_scale - max_scale) &gt; 0.00001:
        v = scale_factor * sphere.scale
        sphere.scale = Vector([max(s, min_scale) for s in v])        
    #scene.update()
    touchers = [s for s in [within_touch(sphere, s) for s in arranged_spheres] if s.length &gt; 0]
    if len(touchers):
        fix = Vector()
        for v in touchers:
            if sphere.color[0]:
                sphere.color[0] -= 0.1
            elif sphere.color[1]:
                sphere.color[1] -= 0.01
            else:
                sphere.color[2] -= 0.001
            fix += v

        #if fix.length &lt; 0.0001 or not inbounds(sphere, loc) or r_count &gt; tries_per_sphere:
        if (fix &lt; 0.0001) or r_count &gt; tries_per_sphere / 2:
            r_count = 0
            random_sphere(sphere)
            resets += 1            
            #sphere.color[1] = 0
        else:
            loc = checkbounds(sphere, fix)
            sphere.location = loc
            moves += 1
        t_count += 1
        r_count += 1
        tries += 1
    else:
        print("Arranged ", sphere.name, "moves:", moves, "rand moves", resets, "fixes:", t_count)
        moves = 0
        t_count = 0
        r_count = 0
        resets = 0
        arranged_spheres.append(sphere)
        sphere =  spheres.pop() if len(spheres) else None

for s in arranged_spheres:
    if not scene.objects.get(s.name):
        scene.objects.link(s)

print("Arranged % d of %d" % (len(arranged_spheres), count))
print("TRIES", tries)






Edit

Following shows size 0.142 spheres placed in a domain 1 cube.  Placed 197 of 200 on this run, with settings shown.



Update for 2.8

import bpy
context = bpy.context
from random import uniform
from mathutils import Vector

count = 200

# random scale 
min_scale = 0.1
max_scale = 0.1
# scale_factor .. scale down to help fit
scale_factor = 0.1
domain = 1
only_spheres = True

tries_per_sphere = 256

use_context_object = False
obj = context.object

def random_vector(a, b):
    return Vector([uniform(a, b) for c in "xyz"])

def checkbounds(sphere, fix):
    #return True
    loc = sphere.location + fix
    #dm = max(sphere.dimensions) / 2
    dm = max(sphere.scale)
    _min = -domain + dm
    _max =  domain - dm
    for i, v in enumerate(loc):
        if v &lt; _min:
            loc[i] = _min
        if v &gt; _max:
            loc[i] = _max
    return loc

def inbounds(sphere, loc):
    #return True
    dm = max(sphere.scale)
    return -domain + dm &lt; min(loc) and  max(loc) &lt; domain - dm

def random_sphere(sphere):
    # random scale and rotation for sphere based on settings
    if only_spheres:
        scale = uniform(min_scale , max_scale) * Vector((1, 1, 1))
    else:
        scale = random_vector(min_scale, max_scale)
    scale = Vector([min(s, domain) for s in scale])
    dom = domain - max(scale)    
    sphere.scale = scale

    sphere.location = random_vector(-dom, dom) 

scene = context.scene
view_layer = context.view_layer
mat = bpy.data.materials.get("randobjcol")
if not mat:
    mat = bpy.data.materials.new("randobjcol")
#    mat.use_object_color = True

boundbox = bpy.data.objects.get("BoundBOX")
if not boundbox:
    bpy.ops.mesh.primitive_cube_add()
    boundbox = context.object
    boundbox.name = "BoundBOX"
    boundbox.display_type = 'WIRE'
    boundbox.hide_select = True
    boundbox.location = (0, 0, 0)

if not scene.collection.objects.get("BoundBOX"):
    scene.collection.objects.link(boundbox)

boundbox.scale = domain * Vector((1, 1, 1))

spheres = []
if use_context_object and obj:
    sphere = obj
else:
    bpy.ops.mesh.primitive_uv_sphere_add()
    sphere = context.object
    #bpy.ops.object.shade_smooth()
    context.collection.objects.unlink(sphere)

sphere.active_material = mat

for i in range(count):
    spheres.append(sphere)
    random_sphere(sphere)    
    sphere = sphere.copy()
    #context.scene.objects.link(sphere)

# update the spheres dimensions    
#context.scene.update()        

tries = 0
t_count = 0
r_count = 0
moves = 0
resets = 0
pfix = Vector()

def within_touch(s1, s2):
    def r(s):
        return max(s.scale)
        #return max(s.dimensions) / 2
        #return sum(s.dimensions) / 3 #  average.

    d = (s2.location - s1.location)

    r1, r2 = r(s1), r(s2)

    radsum = r(s1) + r(s2)
    if d.length &lt;= 0.0001:
        print("SAME")
        return Vector((r1, r1, r1))
        # same position    
    elif d.length &lt; radsum:
        # return a vector to move away
        d.length = radsum - d.length
        return  -d
    return Vector()

sphere = spheres.pop()
arranged_spheres = [sphere]
sphere = spheres.pop()

while sphere and tries &lt; tries_per_sphere * count:
    if t_count &gt; 5 and abs(min_scale - max_scale) &gt; 0.00001:
        v = scale_factor * sphere.scale
        sphere.scale = Vector([max(s, min_scale) for s in v])        
    #scene.update()
    touchers = [s for s in [within_touch(sphere, s) for s in arranged_spheres] if s.length &gt; 0]
    if len(touchers):
        fix = Vector()
        for v in touchers:
            if sphere.color[0]:
                sphere.color[0] -= 0.1
            elif sphere.color[1]:
                sphere.color[1] -= 0.01
            else:
                sphere.color[2] -= 0.001
            fix += v

        #if fix.length &lt; 0.0001 or not inbounds(sphere, loc) or r_count &gt; tries_per_sphere:
        if (fix &lt; 0.0001) or r_count &gt; tries_per_sphere / 2:
            r_count = 0
            random_sphere(sphere)
            resets += 1            
            #sphere.color[1] = 0
        else:
            loc = checkbounds(sphere, fix)
            sphere.location = loc
            moves += 1
        t_count += 1
        r_count += 1
        tries += 1
    else:
        print("Arranged ", sphere.name, "moves:", moves, "rand moves", resets, "fixes:", t_count)
        moves = 0
        t_count = 0
        r_count = 0
        resets = 0
        arranged_spheres.append(sphere)
        sphere =  spheres.pop() if len(spheres) else None

for s in arranged_spheres:
    if not context.collection.objects.get(s.name):
        context.collection.objects.link(s)

print("Arranged % d of %d" % (len(arranged_spheres), count))
print("TRIES", tries)