scipy/scipy — pull request #4295 by I--P (2014-12-19 09:06:22+00:00)
Title: Pole placement

Please find in this PR two pole placement algorithms for linear systems.

I started a code review topic for this on scipy-dev a few days ago and nobody raised any objections so I'm issuing a pull request. Please find below a cut'n paste of the code review request with explanations on what the code does and my concerns regarding licensing issues (it is really important to dig into these imho).

You will see that there are a lot of comments in the code, I hope they are not too verbose but I thought it was important to motivate every choice I made as the original publications are sometimes sloppy on details (see for example the choice of the initial transfer matrix X or what I do when an "improvement" on the transfer matrix yields a null vector, lines 1481 and 1175).

--- CUT'N PASTE STARTS HERE ---

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.

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.