@fujikit, thank you for sharing. it is quite interesting work.
I have looked into your code and I see you are using an analytic solution for IK based on the fact that the shoulder, elbow, and wrist-1 joints are parallel.
The robot has calibration data that slightly changes the DH parameters. Perhaps, you should call get_inverse_kin with the qnear from your solutions to get a more accurate IK solution. or alternatively, you can implement a Jacobian based IK with initial conditions from the analytic solver to consider the calibration data.
check out this post UR10 Calibration.conf DH model generation