Union partially overlapping cubes without boolean modifier
Turns out Python+Blender is not the right tool for the job. I used Python+OpenSCAD instead and so far I have not experienced any issues with the union boolean modifier. Glad to find out there's an easy solution that just works and that there's no need to do any complex stuff in Blender.

From OpenSCAD I can export to STL and import in Blender if I need to do some manual adjustments.