The energy of an array of 3D coherent strained islands on a lattice-mismatched substrate equals: E = Delta E(EL)(V) + Delta E(?)(FACETS)((RENORM)) + Delta(EL)(EDGES) + E(EDGES) + E(INTER), where Delta E(EL)(V), is the volume elastic relaxation energy, Delta E(FACETS)(RENORM) is the change of the surface energy of the system due to the formation of islands, which includes the strain-induced renormalization of the surface energy of the island facets and of the planar surface, Delta E(EL)(EDGES) is the contribution of the island edges to the elastic relaxation energy, E(EDGES) is the short-range energy of the edges, and E(INTER) is the energy of the elastic interaction between islands via the substrate. The energy Delta E(EL)(EDGES) approximate to -L(-2) . ln L always has a minimum as a function of the size of the islands L, and the total energy E = E(L) may have a minimum at an optimum size L(opr). E(INTER) is the driving force for the lateral ordering of 3D islands. Among different arrays of islands on the (001)- surface of a cubic crystal, the total energy is minimum for the periodic square lattice with primitive lattice vectors along the ''soft'' directions [100] and [010]. Thus, a periodic square lattice of equal-shaped and equal-sized 3D islands is, under certain conditions, the stable array of islands which do not undergo ripening. The theory explains the spontaneous formation of ordered arrays of 3D islands in the InAs/GaAs(001) system.