Creating large numbers of unique objects

I am working on a script to generate several thousand rocks. For this I am using a low_poly_rockgenerator add-on. Once this section of the code is complete, I need to apply rigid body physics to each rock for a physics simulation, where the rocks fall into a container and settle (including colliding with each other and the container).

As I add more rocks to the scene, the generation of each rock slows, leading to my desired outcome taking several weeks. I have tried some instancing and dupliverts, but this leaves me with one mesh encompassing all rocks, making the physics simulation impossible.

Ideally I would like to have unique objects in each case but from what Ive read that seems unlikely, I am looking for any advice to help speed up the process of particle generation.

Below is the function which is responsible for the rock generation, with the only variables being the number of rocks, a size distribution for the rocks and a list of positions (previously determined). Various other parameters are determined in the function which are then used in the rock generation

Any advice is greatly appreciated

def rockgenerator(particles,minsize,maxsize,cords):
    for i in range(particles):
        subdivisions=random.randint(4,6)
        radius=random.uniform(minsize,maxsize)
        xskew=random.uniform(0.8,1)
        yskew=random.uniform(0.8,1)
        zskew= 1
        sharpness=random.uniform(0.2,0.5)
        maxangle=random.uniform(90,180)
        Dispstrength=random.uniform(0.25,1.5)
        bpy.ops.mesh.add_low_poly_rock(subdivisions_prop=subdivisions,radius_prop=radius,x_scale_prop=xskew,y_scale_prop=yskew,z_scale_prop=zskew,decimate_sharpness_prop=sharpness,decimate_planar_angle_prop=maxangle,displace_strength_prop=Dispstrength)
        bpy.ops.rigidbody.object_add()
        bpy.context.object.rigid_body.collision_shape = 'CONVEX_HULL'
        bpy.context.object.rigid_body.mesh_source = 'FINAL'
        bpy.context.object.rigid_body.use_margin = True
        bpy.context.object.rigid_body.collision_margin=0
        bpy.context.object.rigid_body.linear_damping = 0.98
        bpy.context.object.rotation_euler[0] = random.uniform(0,6.265732)
        bpy.context.object.rotation_euler[1] = random.uniform(0,6.265732)
        bpy.context.object.rotation_euler[2] = random.uniform(0,6.265732)
        loc=random.choice(cords)
        bpy.context.object.location[0] = float(loc[0])
        bpy.context.object.location[1] = float(loc[1])
        bpy.context.object.location[2] = float(loc[2])
        cords.remove(loc) 
    return cords