I would not use Poly for this implementation; It doesn't work well with time-delay terms.

>>> G = TransferFunction(s*exp(-s), s**2 + 2*s + 1, s)
>>> pole_zero(G)
Traceback (most recent call last):
  File "/home/namannimmo/oss/sympy/sympy/polys/polyutils.py", line 211, in _parallel_dict_from_expr_if_gens
    monom[indices[base]] = exp
KeyError: exp(-s)

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/namannimmo/oss/sympy/sympy/physics/control/control_plots.py", line 28, in pole_zero
    num_poly = Poly(system.num, system.var).all_coeffs()
  File "/home/namannimmo/oss/sympy/sympy/polys/polytools.py", line 164, in __new__
    return cls._from_expr(rep, opt)
  File "/home/namannimmo/oss/sympy/sympy/polys/polytools.py", line 293, in _from_expr
    rep, opt = _dict_from_expr(rep, opt)
  File "/home/namannimmo/oss/sympy/sympy/polys/polyutils.py", line 368, in _dict_from_expr
    rep, gens = _dict_from_expr_if_gens(expr, opt)
  File "/home/namannimmo/oss/sympy/sympy/polys/polyutils.py", line 307, in _dict_from_expr_if_gens
    (poly,), gens = _parallel_dict_from_expr_if_gens((expr,), opt)
  File "/home/namannimmo/oss/sympy/sympy/polys/polyutils.py", line 216, in _parallel_dict_from_expr_if_gens
    raise PolynomialError("%s contains an element of "
sympy.polys.polyerrors.PolynomialError: exp(-s) contains an element of the set of generators.

The reason behind using Poly() was speed. I was using NumPy for root-locus and getting a huge performance advantage because- (a) SymPy is a symbolic computation library and NumPy is a numeric computation library. (b) Matplotlib requires NumPy for internal numeric computation and hence NumPy is also a listed dependency for matplotlib. It is installed with matplotlib.

Also just to be clear, these plots will not support transfer functions with time delay. Maybe I'll allow time-delay in Pole-Zero because time delay doesn't affect the poles and zeroes of the system.

It is possible to use sympy with other plotting backends that do not necessarily use numpy:
https://pypi.org/project/sympy-plot-backends/

Generally when sympy has numeric routines they are more accurate than numpy's:

In [8]: p = expand(prod(x - i for i in range(1, 21)))

In [9]: nroots(p)
Out[9]: 
[1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15
.0, 16.0, 17.0, 18.0, 19.0, 20.0]

In [10]: numpy.roots(Poly(p).all_coeffs())
Out[10]: 
array([19.99980929, 19.00190982, 17.99092135, 17.02542715, 15.94628672,
       15.0754938 , 13.91475559, 13.07431403, 11.95328325, 11.02502293,
        9.99041304,  9.00291529,  7.99935583,  7.000102  ,  5.99998925,
        5.00000067,  3.99999998,  3.        ,  2.        ,  1.        ])

Plotting backend is not an issue here. Whichever backend we may use, it requires efficient sampling to generate plots. Some form of numerical computation will always occur under the hood. SymPy undoubtedly is more accurate in numerical operations. Still, during sampling with a large number of points (10000), accuracy hardly matters, especially when a much faster approach (more than 1000 times) is accurate up to 3-4 decimal places.

is accurate up to 3-4 decimal places.

The accuracy can degrade much more than that.

The idea with using a symbolic library is that you should have greater confidence in things like the accuracy of any computed result even if that means waiting longer for things to be computed. It also means that different approaches are often better to make something efficient.

As an example your root_locus function computes the roots of 10000 different polynomials. Why 10000? No one needs that many points in the plot. Many of the lines are just joining points on the real axis. To plot a straight line you only need to know its start and end and those can be computed exactly. Then you don't need to evaluate the polynomial at 10000 points.

I'm not sure I've understood the calculation that you want precisely but I think it's something like this:

In [35]: eq
Out[35]: k⋅(s + 7) + s⋅(s + 5)⋅(s + 15)⋅(s + 20)

In [36]: discriminant(eq)
Out[36]: 
      4          3             2                               
- 27⋅k  + 47088⋅k  - 37907500⋅k  + 3881250000⋅k + 1265625000000

In [37]: nroots(discriminant(eq))
Out[37]: [-130.564918275438, 307.765516989741, 783.399700642849 - 743.513266939402⋅ⅈ, 783.399700642849 + 743.513266939402⋅ⅈ]

The positive real roots for k are the points on the real axis where lines start and end. This now shows you what values of k are interesting for the plot.

For TransferFunction(2*s**2 + 5*s + 1, s**2 + 3*s + 5, s)

• If num=10000 (sampled 10000 times, default)

  

  The scatter plot looks like this.
  And the actual root-locus plot (corresponding to the regression of these sampled points) looks like this-

  

Similarly,

• If num=1000

  

  

• If num=300

  

  

To plot a straight line you only need to know its start and end and those can be computed exactly. Then you don't need to evaluate the polynomial at 10000 points.

Straight Lines in root locus plots are not very common. In root-locus plots, we plot the movement of the roots of the polynomial den+ k*num when k varies from 0 to inf. As you can tell when k is 0 we get the roots of den (poles) and as k approaches inf, the roots of den + k*num are roots of num only (zeros). Therefore, when num and den have the same degree, the root locus plot is a path from poles to zeros of the transfer function. This is not always a straight line.

Only one is a straight line path here and knowing this fact can't help us in any way. I don't understand how knowing the discriminant of the polynomial will help us in knowing the points on the rl plot. Maybe I'm missing something.

Straight Lines in root locus plots are not very common.

Your example with num=1000 illustrates the point I wanted to make perfectly. To generate the plot you compute the roots of the polynomial for 1000 different values of k ranging from 0 to kmax where kmax is arbitrarily chosen to be 400. Out of those 1000 values of k only 4 have non-real roots so the other 996 polynomials whose roots are computed are only used to plot two straight lines. The straight line could be computed more efficiently and more accurately by simply computing its start and end points which you can do like this:

In [93]: k, s = symbols('k, s')

In [94]: n, d = 2*s**2 + 5*s + 1, s**2 + 3*s + 5

In [95]: eq = d + k*n

In [96]: eq
Out[96]: 
  ⎛   2          ⎞    2          
k⋅⎝2⋅s  + 5⋅s + 1⎠ + s  + 3⋅s + 5

In [97]: nroots(n) # k = oo
Out[97]: [-2.28077640640441, -0.219223593595585]

In [98]: discriminant(eq, s)
Out[98]: 
    2            
17⋅k  - 14⋅k - 11

In [99]: nroots(discriminant(eq, s))
Out[99]: [-0.491899499749248, 1.31542891151395]

In [100]: roots(discriminant(eq, s), multiple=True)
Out[100]: 
⎡7    2⋅√59  7    2⋅√59⎤
⎢── - ─────, ── + ─────⎥
⎣17     17   17     17 ⎦

In [101]: eq.subs(k, roots(discriminant(eq, s), multiple=True)[1])
Out[101]: 
 2         ⎛7    2⋅√59⎞ ⎛   2          ⎞    
s  + 3⋅s + ⎜── + ─────⎟⋅⎝2⋅s  + 5⋅s + 1⎠ + 5
           ⎝17     17 ⎠                     

In [102]: roots(eq.subs(k, roots(discriminant(eq, s), multiple=True)[1]), s)
Out[102]: {-9 + √59: 2}

In [103]: list(_)[0].n()
Out[103]: -1.31885425213139

Your suggestion to use num=10000 and compute the roots of 10000 polynomials is a workaround for having the wrong value of kmax. In fact by computing the roots of the discriminant it is possible to compute what the relevant values of k are.

Fine. I get your point now. But there is still an issue. You see, with your calculations, we get these points:

Which is basically poles and zeroes of the Transfer Function + The point where the discriminant is zero (equal roots). I agree that for the two straight lines on the axis, this will work fine but what about the two symmetric arms of the curve?
We cannot plot them by simply knowing the start and endpoints. We still need to compute the roots of the polynomials from 0 to k_c (value of k for which the 0 discriminant point is reached) in uniform intervals.

We still need to compute the roots of the polynomials from 0 to k_c

Yes, but we don't need 10000 points to plot the curve. The num=10000 plot above shows that an acceptable plot is obtained from only 40 points (since the other 9960 points are on the real axis).

I expect that uniform spacing is not optimal. Naturally I would guess that something like logarithmic or harmonic spacing would be better. You can see in the num=10000 plot that using uniform intervals for k does not lead to a uniform distance between the points in the plot.

Ideally the curve would be generated using pseudo-arc-length continuation so you can generate them at a fixed 2D distance within the plot but that's probably more complicated than it is worth.