[![enter image description here][1]][1]
_Test run, settings as posted_

**https://blender.stackexchange.com/questions/89643/algorithm-for-random-distribution-of-circles-in-a-square**

> This simple script to place a number of `num_circles` circles with
> radius `radius` in a refined rectangle (aligned to global axes) speficied
> by `xmin, xmax, ymin` and `ymax` by trying random locations until it not
> colliding with an existing circle. This is certainly not the most
> efficient way of doing it, but it should be fine for smaller problems
> and if the relative area covered is not too high.


Test code after making a couple of edits to https://blender.stackexchange.com/a/89658/15543 to make it add cylinders of random height to a rectangle.

    from mathutils import Vector, Matrix
    from mathutils.noise import random, seed_set
    import bpy

    # Specify parameters
    seed = 0
    xmin = -10
    xmax = 10
    ymin = -10
    ymax = 10
    radius = 1
    num_circles = 68
    max_tries = 10000
    max_height = 4
    # Init
    seed_set(seed)
    sx = xmax - xmin - 2 * radius
    sy = ymax - ymin - 2 * radius
    xminm = xmin + radius
    yminm = ymin + radius
    existing_locations = []
    sce = bpy.context.scene
    bpy.ops.mesh.primitive_cylinder_add(
        location=(0, 0, 0),
        depth=1,
        radius=radius)
    ref_circle = bpy.context.object
    ref_circle.data.transform(Matrix.Translation((0, 0, 0.5)))

    # Loop
    for i in range(num_circles):
        j = 0
        searchOn = True
        while searchOn:
            if j > max_tries:
                bpy.ops.object.select_all(action='DESELECT')
                ref_circle.select_set(True)
                bpy.ops.object.delete()
                raise ValueError('Found no more room for another circle')
                break
            j += 1
            new_location = (xminm + random() * sx,
                            yminm + random() * sy,
                            0)
            for existing_location in existing_locations:
                if (Vector(existing_location) - Vector(new_location)).length < 2 * radius:
                    break
            else:
                new_circle = ref_circle.copy()
                new_circle.dimensions.z = max_height * random()
                new_circle.location = new_location
                sce.collection.objects.link(new_circle)
                existing_locations.append(new_location)
                searchOn = False

    ref_circle.select_set(True)
    bpy.ops.object.delete()
    sce.update()


  [1]: https://i.sstatic.net/98m94.png