[SciPy-Dev] Asking for a code review: pole placement

[Irvin Probst]

Hi all,
I've finished to integrate and test my pole placement code, here is the 
url of the diff. I'd like to know what you think before issuing a pull 
request.

https://github.com/I--P/scipy/compare/pole-placement

As said a few days ago there are two algorithms, KNV method 0 (KNV0) and 
Yang Tits (YT). See code for urls to the reference papers.

Motivations for doing this in pure Python rather than linking some old 
fortran code:

I was in need of a 100% Python dependency free (except Numpy of course) 
pole placement function, integrating it into scipy was not my primary 
goal. Fortran code used for decades will for sure be faster and bug free 
but it was not an option.

Regarding KNV0:

- Licensing issues: I chose this method before knowing Matlab used it 
too. I did read some old matlab source code (circa 2001 release iirc) at 
the beginning of my work I won't lie on that, after Ralf Gommers warned 
me it was not possible to port Matlab code I used my knowledge of their 
code to stay as far as possible from what they did (the algorithm gives 
some degrees of freedom). Given the current status of my code I can't 
see any reason why it would be seen as a rip-off of their code unless 
I'm banned from coding a projection into a vector space because I read 
some Matlab code doing it. See lines 1171, 1172 and 1173 in ltisys.py 
which are the core of this algorithm and are not in any way an idea 
coming from Mathworks folks, these formulas are clearly written in the 
reference publication page 1145 and 1146.

- Complex poles are not supported by the original algorithm. Using the 
ideas from the second algorithm below I managed to reach a success rate 
of 90% with complex poles but this is not reliable enough so the 
relevant code is disabled. Matlab has a 100% success rate in their 2013 
implementation but I don't know how they achieved that (I _*did not*_ 
check their 2013 code).

Regarding Yang Tits algorithm:

- This is an update of the KNV0 algorithm published in the mid 90's, as 
far as I know it is only distributed by Slicot as "contributed MATLAB 
functions" http://slicot.org/contributed-software under a non-free license.

- This implementation is 100% genuine, I did not read the code 
distributed by slicot past the license header until my code was working 
so there are no licensing issues imho. After achieving a working 
implementation of this algorithm I checked the reference implementation 
from Slicot and I'm still puzzled by how they did it. I'm sure there are 
tons of good reasons why their implementation is faster/better than mine 
but I simply can't understand what they did.

- Real poles are supported, sometimes the solution has a better 
conditionning than KNV0's, sometimes not, but it does converge must 
faster all the time. I was considering adding an option to return the 
solution with the best conditionning from the two methods but it might 
confuse users as it would only work on real poles. What do you think ?

- Complex poles are supported, I ran it for hours with random matrices 
and poles and it never failed.

Please let me know if this code match you quality standards and I would 
really be glad if someone could test it and report any bug. The 
reference matlab implementations of these algorithms have been used and 
polished for decades, I'm sure I've missed lots of specific cases.

As a sidenote, if PR #4249 is accepted it should be possible to greatly 
improve the speed of these two algorithms.