In my model I fixed ShockBtmRF_001 to KuckleRF_001 with two circularEdges - it is not necessary (in previous model two parts was merged), but I like to use it that way to simulate wheel alignment setting using eccentric fasteners between that two parts.
But it not caused serious convergency problems.
When I added green axles, the problems start to rise and right now, just when you load model in it default position solver is unable to solve it.
In some orientations it still works, but really hard.
I also took a look into source of your solver.
I like how it is written - nice clean code, that I can read and understood even with nearly lack of Python knowledge.
But maybe it is good time for me to finally learn Python ?
(as a C/C++ programmer I can't get used to Python indentations)
I was able to do step-by-step simulation by changing solverControlData to do more smaller steps of pos/spinAccuracy, and not requiring accuracy below 0.05, but I looked around whole your code.
I spent only 2h looking into it. so just got some basics of your concept, maybe you would like to use some of my experience with similar problem, that I resolved with similar idea as your.
Main difference was that I was not applying 'move' to rigids to lower an error, but I was applying just a force like a spring - when error was bigger, there was bigger force applied for both rigids to bring then closer to proper position.
Each rigid body had own mass and moment of inertia (I think for such FreeCAD solver it could even be a random - I love randomness in numeric simulations) - so applied force changes a speed of rigid body and rotation around body mass center (when applied not directly into mass center it generates torque).
It would end with endless vibrations, so each constrain also applying a dump force. Each constrain measuring it own stretching/squeezing (could be also rotating etc) speed measured just like: (last_frame_length-this_frame_length)/frame_time and applying dumping force proportional to that speed, but with opposite direction.
To stabilize calculations I added in similar way as you, some 'play' distance that was permitted error for each constrain.
Finally every rigid calculates resultant force and torque that changes speed and rotation of this rigid.
Then each timeframe results in changing position and rotation of each rigid.
With ~100 simple constrains, and few thousands steps per frame I was able to solve it still with 60fps (and it was 10 years ago on one thread) - for sure my constrains was simpler, and permitted errors was higher.
I think it could be good idea to use more physical-like simulation - changing speed and rotation between each convergency step should give result even faster.