Using zpk with complex poles requires specific order for conjugates

[albheim]

Is this documented somewhere?

Cause I wanted to use zpk to create a system and had some complex values that I threw in in some random order and got this error about not having matching complex conjugates. Looking at the code it seems that the roots2real_poly_factors method requires the pole with positive imaginary part to come first among the two and the other right after. I don't really see the reason why we require an order at all? It should not be very tricky to just keep a set of non-matched complex poles and check through them for matches for every new complex pole we hit, right? Or at least we should document it in the docstring for zpk so people understand why they get this error.

[mfalt]

The reason behind this is that it is not trivial to do the matching when there are numerical errors. Say for example poles at (1+x*im),(1-x*im),1+0*im for some very small x. If you are not very careful, under numerical errors, you might match (1+x*im)and 1+0*im. But now you are left with (1-x*im) which is not close enough to the imaginary axis to be considered "real", and you therefore end up with a complex valued system, or an error, even when the system you were manipulating was real. The good thing is that for real matrices, eigenvalue solvers like LAPACK will reduce their problems to 2x2 matrices, and therefore always know which eigenvalues belong together. The result is that we can rely on this specific structure internally.

This became a problem when julia changed the behaviour of the eigenvalue solver in v1.2, to do some random ordering by default. You can therefore find calls to eigvalsnosort (found in utilities.jl), which asks julia not to sort them, and we therefore get the order as supplied by LAPACK.

I do agree however, that we should either document this requirement, or sort the poles that the user supplies when creating a system. Especially if the user tries to create a TransferFunction{Continous, SisoZpk{<:Real,<:Complex}}, for example in https://github.com/JuliaControl/ControlSystems.jl/blob/4d1efa16ba8a35263203c2678846f06af7803023/src/types/SisoTfTypes/SisoZpk.jl#L117

[olof3]

I think that the error message from zpk is already reasonable clear about what the cause and the fix is, but it could definitely be made more clear and documented further.

julia> zpk([], [-1 + im, 1, -1 - im], 1)
ERROR: AssertionError: zpk model should be real-valued, but poles do not come in conjugate pairs.

I don't think this is a big nuisance, but if many others do, then I'm not against fixing this. We could have it so that that the zpk method where the user specifies poles and zeros does some sorting, but those that convert StateSpace and TransferFunction models do not. I guess it comes down to how forgiving the tool box should be vs. amount of code + efficiency + maintenance. In this instance the penalty is certainly quite small.

If we don't make this change, we should at least change so that the order of e.g. -1 + im and -1 - im doesn't matter, right now zpk([], [-1 - im, -1 + im], 1) gives the terrible error message ERROR: AssertionError: Found pole without matching conjugate.

Perhaps there should also be better help for creating complex-coefficient system with zpk.

[albheim]

The reason behind this is that it is not trivial to do the matching when there are numerical errors. Say for example poles at (1+x*im),(1-x*im),1+0*im for some very small x. If you are not very careful, under numerical errors, you might match (1+x*im)and 1+0*im. But now you are left with (1-x*im) which is not close enough to the imaginary axis to be considered "real", and you therefore end up with a complex valued system, or an error, even when the system you were manipulating was real.

Ahh, makes sense. Didn't think too far there...

I do agree however, that we should either document this requirement, or sort the poles that the user supplies when creating a system. Especially if the user tries to create a TransferFunction{Continous, SisoZpk{<:Real,<:Complex}}, for example in...

I don't this this is a big nuisance, but if many others do, then I'm not against fixing this. We could have it so that that the zpk method where the user specifies poles and zeros does some sorting, but those that convert StateSpace and TransferFunction models do not. I guess it comes down to how forgiving the tool box should be vs. amount of code + efficiency + maintenance. In this instance the penalty is certainly quite small.

Sounds reasonable to me, if we assume the user supplies numbers without any numercial funniness we could just ad something like this sort!(poles, by = z -> (abs(imag(z)), real(z), imag(z))) which should first clump up things by imaginary an real part and then order them in negative positive imaginary part order within the pairs.

If we don't make this change, we should at least change so that the order of e.g. -1 + im and -1 - im doesn't matter, right now zpk([], [-1 - im, -1 + im], 1) gives the terrible error message ERROR: AssertionError: Found pole without matching conjugate.

Yeah, this is the error I ran into and if someone who is not keen on checking the source encounters this it is a bit confusing since there clearly is a conjugate pair.

[olof3]

Yeah, this is the error I ran into and if someone who is not keen on checking the source encounters this it is a bit confusing since there clearly is a conjugate

Let's fix this asap.

if we assume the user supplies numbers without any numercial funniness we could just ad something like this sort!(poles, by = z -> (abs(imag(z)), real(z), imag(z)) which should first clump up things by imaginary an real part and then order them in negative positive imaginary part order within the pairs.

Doesn't this run into the problems for double complex poles?

Perhaps it is better after all to try to pair the roots up, and in the case that they don't come after each other, to just search through the reminding roots to see if there is a match somewhere. Using a vectors of Bools to keep track of the roots that have been used, this approch would have a no impact if they do come in pairs, and be reasonably fast if they come in an arbitrary order.

[olof3]

But now you are left with (1-x*im) which is not close enough to the imaginary axis to be considered "real", and you therefore end up with a complex valued system

Unless the gain k is complex you will never end up with a complex-coefficient system, but you'll get an error instead. (This is how it should work at least)

[albheim]

if we assume the user supplies numbers without any numercial funniness we could just ad something like this sort!(poles, by = z -> (abs(imag(z)), real(z), imag(z)) which should first clump up things by imaginary an real part and then order them in negative positive imaginary part order within the pairs.

Doesn't this run into the problems for double complex poles?

Yeah, you are right.

Perhaps it is better after all to try to pair the roots up, and in the case that they don't come after each other, to just search through the reminding roots to see if there is a match somewhere. Using a vectors of Bools to keep track of the roots that have been used, this approch would have a no impact if they do come in pairs, and be reasonably fast if they come in an arbitrary order.

How do we do the matching in case we have numerical problems though? Don't we need to check all numbers and make sure all are matched in a was that minimizes differences between matches? Otherwise we could imagine having a chain where 1+im matches to 1+eps-im, and then 1+eps+im is matched to 1+2eps-im and so on. Then we might be left with 1-im and 1+N*eps+im. Maybe not a very likely scenario, and maybe not even something that is bad enough to care about?

Do we want users to be able to supply poles/zeros in any order? And do we assume user supplied poles/zeros are nice in that they can easily be sorted? Or should we just add in the docstring that complex numbers should be supplied in positive imaginary part first and negative directly after.

Unless the gain k is complex you will never end up with a complex-coefficient system, but you'll get an error instead. (This is how it should work at least)

Isn't that exactly what @mfalt wanted to show though? If you don't do the matching correct you might have a system that should work but will error because a pole is left over so it looks like it is complex?