Я бы, вероятно, использовал бы только sympy
и их решатель: sympy Library Solver
from sympy import solve, log, exp
from sympy.abc import x,y,z,d
f = (x - y) / log((x - z) / (y - z)) - d
solution = solve(f, x)
И вывод дает мне
[-d*LambertW(-(y - z)*exp(-(y - z)/d)/d) + z]