Re-write utils.polygcd, utils.polylcm and utils.multi_polylcm and create utils.multi_polygcd

The method currently used to calculate polynomial greatest common divisor and lowest common multiple involves calculating the roots of the polynomials then using logic to determine the gcd and lcm. This method only works well for lower order polynomials (up to approximately 5th order), due to the inherent difficulty in accurately calculating polynomial roots for higher order polynomials. This causes erroneous results in many of the functions in utils when applied to large mimo systems (bigger than 2x2), including utils.poles, utils.zeros and utils.tf2ss, as many of the interim calculations involve large polynomials. A new method should be created that does not rely on calculating roots. Also, Euclid's algorithm (or any method closely related to this) should not be used as it is not applicable to polynomials in floating point arithmetic.

@ilayn has given the following recommendation:

For LCM, I'm using a similar result to Karcanias, Mitrouli, "System theoretic based characterisation and computation of the least common multiple of a set of polynomials", Lin Alg App, 381, 2004. The results is acceptable but never stress tested it to the extremes.

For GCD, Sylvester-matrix based methods are sufficiently good for medium difficulty/size problems.

The following code can be used to test if the new method is working correctly:

import utils
import numpy
s = utils.tf([1,0], 1)
G = utils.mimotf([[1/((s+1)*(s+2)), 1/((s+2)*(s+3)), 1/((s+3)*(s+4)), 1/((s+4)*(s+1))],
                 [1/((s+1)*(s+2)), 1/((s+2)*(s+3)), 1/((s+3)*(s+4)), 1/((s+4)*(s+1))],
                 [1/((s+1)*(s+2)), 1/((s+2)*(s+3)), 1/((s+3)*(s+4)), 1/((s+4)*(s+1))],
                 [1/((s+1)*(s+2)), 1/((s+2)*(s+3)), 1/((s+3)*(s+4)), 1/((s+4)*(s+1))]])
poles = utils.poles(G)
assert poles.shape == (4,)
#There should be exactly 4 poles
Gpoles = [G(i) for i in poles]
AbsoluteMaximums = numpy.array([abs(i).max() for i in Gpoles])
SmallestMaximum = Maximums.min()
assert SmallestMaximum > 10**6
#Evaluating G(poles) should give at least one very large (infinite) element for each pole

alchemyst / Skogestad-Python

Re-write utils.polygcd, utils.polylcm and utils.multi_polylcm and create utils.multi_polygcd #295