I have worked out (with manual debugging with print statements) that the problem was with the following code part:
name_cluster_head = '{0}_{1}'.format(i, output_name)
out_cluster_head = os.path.join(output_fds, name_cluster_head)
arcpy.CopyFeatures_management(facilities_sublayer, out_cluster_head)
More specifically with arcpy.CopyFeatures_management. I have not really found any related issues online, so I have tried to avoid the function with:
arcpy.FeatureClassToFeatureClass_conversion(facilities_sublayer, output_fds, name_cluster_head)
This resulted in the following error message:
So it became clear that the problem was in the naming. The name starts with the number and it somehow creates the problem with arcpy. After changing to:
name_cluster_head = 'Cluster_head_{0}_{1}'.format(i, output_name)
The issue was solved.