Improvement on Bode plot phase

[tallakt]

Just added some code to allow bode plots to start with the (heuristic) correct phase angle. Some systems with double integrators would start with phase +180 which makes the plot kind of nonstandard

[mfalt]

I think this looks reasonable. @baggepinnen @olof3 ?
@JuliaControlBot test-plots

[olof3]

I agree that a poorly unwrapped phase is very annoying in Bode diagrams.

With two integrators the phase currently starts at -180, but you mean that you have some variation of this where the phase starts at +180? A mwe is always good.

Some thoughts,

• It would break for non-rational systems (there should be a test for this!).
• The tolerances for finding zeros and poles at the origin seem very generous. There is a general problem of figuring out the time scales for the plotting functions so would be good with a holistic fix.
• I'm a little bit worried about the numerics of computing zeros and poles of poorly conditioned systems of somewhat higher order. At least with the current naive implementation of zpk. Perhaps tzero would be better.

I would guess that a reasonably robust solution is a bit complicated to get together, but should definitely be doable.

[baggepinnen]

Good points about the tolerances, those should probably be keyword arguments and documented so that one does not get confused if one has very slow dynamics that gets treated as an integrator.

Maybe all the added computations should be opt out so that one can call bode with maximum performance and no surprises.

[coveralls]

Coverage decreased (-0.2%) to 55.463% when pulling dab68bc on tallakt:master into d183cb1 on JuliaControl:master.

[tallakt]

I have fixed the things commented here as far as I can tell. All the suggestions were fine. To me these were details. Having a code that runs faster will always have lower priority than producing the correct plot. Making a bode plot that plots phase -180 as +180 will just not go down well with most users.

The code right now is at the stage I can make it without going really deep into the code. There still are some unresolved questions like why did I have to use zpk_sys.matrix[1] rather than looking at any other values that may be there. Also I have not added any unit tests to this. I'm afraid its out of my time and Julia knowhow budget at this point.

I would not be afraid to recommend including this pull request, or at leat reimplement a code to plot bode plots with the correct phase.

a simple test code to show an errored phase would be:

s=tf("s")
bodeplot(1/s^3)

[olof3]

I'd guess the behavior is inconsistent for MIMO systems? And an error would still thrown for non-rational systems, e.g., DelayLtiSystem? These problems do not seem like details to me.

The naming is inconsistent, for example what count_zeros(zpk_sys.matrix[1].p)) does is counting the number of complex numbers that lie on the imaginary axis (actually poles, which is even more confusing). It be more clear with just count(on_imag_axis, zpk_sys.matrix[1].p)) and then defining on_imag_axis and on_unit_circle. This is also consistent with what is done in hinfnorm. I think it is fine with duplicated definitions of these functions in order to wrap the precision arguments.

init is not a descriptive variable name, phase_init, phase0, \theta_0 would give you an idea what it's about. It would also be good if all the functionality for finding the initial phase was in it's own function. That would make the code more readable and easier to test.

Shouldn't the tolerances also include w[1] (in the case of continuous time)? Or is it perhaps actually this value that should be the tolerance (unless w[1] == 0)?

I'm also current method for finding poles and zeros is probably not robust for systems of moderate order (>10 or so). Not sure if deterministic incorrect phase wrapping or unexpected phase wrapping is to be preferred? Perhaps it should only be the default for systems of low order.

[tallakt]

I'd guess the behavior is inconsistent for MIMO systems? And an error would still thrown for non-rational systems, e.g., DelayLtiSystem? These problems do not seem like details to me.

I agree that all supported kinds should be supported or at least fallback to the original behavior. I don't have an overview over how many types there are. I have omitted a DelayLtiSystem from being processed by the heuristic. I cannot find a way to get the transfer function without the delay to use in the zpk function.

The naming is inconsistent, for example what count_zeros(zpk_sys.matrix[1].p)) does is counting the number of complex

You are overcomplicating. count_zeros count values with absolute value zero. count_unit_circle counts the number of imaginary values that are of absolute value 1. So the naming is very correct and descriptive. The meaning you mention only makes sense in the context that the function is used. I dont mind if these are rewritten to be more clear, but I will not do it in this case.

init is not a descriptive variable name, phase_init, phase0, \theta_0 would give

As we are in the context of the unwrap! function, init is quite sufficient and IMHO better than the other suggestions. reduce uses the init name to the same effect

I must say I am taken aback by the hostile reactions I am getting here, nitpicking at details while there is a glaring hole in the functionality of ControlSystems. To say it like it is, I will probably not try to create further pull requests for this project. But on the other hand, I do use the project, and I like it, and getting the bode plot to make sense would allow me to share them with my colleagues. Right now it's just too embarrasing to present a bode plot with positive phase at w=0. So I hope you will at least accept the work thus far, or perhaps reimplement in a fashion to suit your tastes.

And honestly thanks for the constructive feedback I have had.

[mfalt]

I must say I am taken aback by the hostile reactions I am getting here, nitpicking at details while there is a glaring hole in the functionality of ControlSystems. To say it like it is, I will probably not try to create further pull requests for this project. But on the other hand, I do use the project, and I like it, and getting the bode plot to make sense would allow me to share them with my colleagues. Right now it's just too embarrasing to present a bode plot with positive phase at w=0. So I hope you will at least accept the work thus far, or perhaps reimplement in a fashion to suit your tastes.

I am sorry that you feel this way. We very much appriciate input and help for anyone that is interested.
I think part of the problem here is that we (or at least I) don't see the current behavior as a critical problem, which makes the drawbacks problematic. It is possible that I am wrong here though.
I talked a bit with @olof3 about this over lunch. We agree that it would be good if 1/s^3 started at -270 instead of +90. However, the current solution comes with quite a few drawbacks. I do think however, that it is a bit early/unnessesary to start discussing whether the naming of the variables and functions are perfect. The following are the main drawbacks that I can see, that have to be weighed against the benefit, or be fixed.

1. It requires calculating the poles and zeros of the system (which can be comparatively expensive, solving eigenvalues)
2. It doesn't work on delay systems (or maybe as importantly, it is not consisten between system types)
3. It doesn't work on MIMO systems (you mentioned that you had to use zpk_sys.matrix[1], that is the transfer function to output 1 from input1 (matrix[1,1]) )
4. It doesn't respect the type-parameters of sys or w (for example, promoting Float32 to Float64) This seems to be fixed in the unwrap function, although it is wrong for MIMO systems, it currently only useds the phase of one of the systems, and even for SISO systems, prev seems to change between a scalar and an array, causing type instability and possibly performance issues.
5. It is numerically unstable

That doesn't mean that we should't implement this.

1. As mentioned above, could be solved by a kwarg, possibly that defaults to the new suggestion. I have no problem with defaulting bode to something slower, if it can be disabled, since it is often used (but not nessesarily, for plotting), and since freqresp is a fast alternative.
2. Might not be critical, as long as we document the feature properly.
3. Can probably(?) be fixed quite easily, either by you, or by some else.
4. Should definetely be fixed, and with a keyword argument for enabling it, we need to make sure that it is type-stable
5. Is a bit trickier. It will produce the "wrong" result when the poles are close to 0. This is a drawback that might be ok.
   It will also create a discontinuity in the function that wasn't present before, when a pole or zero moves towards zero. This could be a problem for example when optimizing over the bode function. See for example https://github.com/JuliaControl/ControlExamples.jl/blob/master/autotuning.ipynb (which might be outdated) on how we sometimes do automatic differentiation and optimization over outputs that you might not expect. How problematic this last point is could be discussed, and if we allow to disable the suggested feautue, it doesn't have to ba a problem.

I want to thank you again for finding this, and helping in creating a solution to the problem. For this package to grow and continue to be useful, we need people. Although far from everything in this package is perfect, we try our best to keep it both efficient and generic enough to be used with other packages in the Julia ecosystem, this is probably why we have so many comments.

[baggepinnen]

I think Mattias summarizes quite well in his last post. We have developed and maintained this toolbox over quite some years now, and it has lately been driven by what we need (at least for me). We also have a lot of workflows using this toolbox, both in research and in teaching. The reason I gave feedback on potential performance regressions in the implementation of bode is that we use the function to produce interactive graphical user interfaces, as well as as optimization targets in both sys.id. and controller optimization. Solving what I consider mostly an aesthetic issue does not warrant significant regressions in performance for a function which is used so extensively in performance sensitive applications.

We have also worked quite hard trying to make the toolbox type stable and not performing unintended type promotions. Indeed, one of the largest PRs pulled the last few years handled this specifically
#118
The conversation around this PR (and most other on this project as well) was also quite extensive, and served exactly the intended purpose: many eyes with different perspectives looking at and evaluating the contributions in order to end up with something more robust.

Accepting a PR which messes with types and performance just to introduce more maintenance burden to fix it when it breaks some existing workflow is not sustainable, hence the extensive feedback on what seems like "details" to some.

To finish off, there definitely should be a test that exposes the previous, unwanted behavior.

If our comments came off as hostile, that's very unfortunate, and I don't think anyone intended to alienate anyone. This project is important to many, and the fact that three maintainers gave detailed reviews on a PR within so short time hopefully reflects how keen we are on making this toolbox better.

[tallakt]

The difficult thing about this discussion is that I received both excellect feedback and what I consider nitpicking. It can’t be generalized.

I dont consider a bode plot starting at phase +90 degrees for ‘1/s^3’ an aestetic issue. If you believe this, then perhaps ask some people who use ControlSystems to make bode plots. ‘1/s^3’ has 270 degree phase delay. Not 90 degree advance. I can tell you my two colleagues who saw these plots would not consider them seriously for this reason.

I guess this ends up a discussion: is there a difference between a double derivation and a double integration for phase. Yes - and thats why you have unwrap! there in the first place

Having this problem is also quite a big deal if you plot an uncompensated system along with a compensated system in the same plot. With slightly bad luck (that is most times for the system Im working on) there is almost 360 degree offset between the two phase plots.

I appreciate a lot of the feedback I got, and you might see that I already went through many hoops to satisfy your requests. But at one point enough is enough. If you are so nitpicky just say you will do it yourselves. Or decide not to do it, thats ok too.

If you want me to continue working on this, perhaps provide constructive feedback such as: «We have several types of systems as defined in files ..., we need the bode function to run for all of these»

[tallakt]

The difficult thing about this discussion is that I received both excellect feedback and what I consider nitpicking. It can’t be generalized.

I dont consider a bode plot starting at phase +90 degrees for ‘1/s^3’ an aestetic issue. If you believe this, then perhaps ask some people who use ControlSystems to make bode plots. ‘1/s^3’ has 270 degree phase delay. Not 90 degree advance. I can tell you my two colleagues who saw these plots would not consider them seriously for this reason.

I guess this ends up a discussion: is there a difference between a double derivation and a double integration for phase. Yes - and thats why you have unwrap! there in the first place

Having this problem is also quite a big deal if you plot an uncompensated system along with a compensated system in the same plot. With slightly bad luck (that is most times for the system Im working on) there is almost 360 degree offset between the two phase plots.

I appreciate a lot of the feedback I got, and you might see that I already went through many hoops to satisfy your requests. But at one point enough is enough. If you are so nitpicky just say you will do it yourselves. Or decide not to do it, thats ok too.

If you want me to continue working on this, perhaps provide constructive feedback such as: «We have several types of systems as defined in files ..., we need the bode function to run for all of these»

[olof3]

We are all very thankful that you take you have taken the time to contribute on a problem that annoys you, in that case there are surely many other people that are annoyed by it. As it turns out there was no trivial fix for this problem, so it is natural that we discuss potential problems and what to do about them.

Back to the technical things,

You are overcomplicating. count_zeros count values with absolute value zero. count_unit_circle counts the number of imaginary values that are of absolute value 1.

My bad, the naming is consistent, but is the behavior so? Shouldn't s = 0 in continuous time correspond to z = 1 in discrete time (and not abs(z) = 1)?

As we are in the context of the unwrap! function, init is quite sufficient and IMHO better than the other suggestions. reduce uses the init name to the same effect

To me it is not clear what init is, in many instances I think it would be convenient to unwrap from the end or from the middle of the frequency response, so it would be reasonable to assume that init corresponds to starting index for the unwrapping.

[olof3]

@mfalt

2. It doesn't work on delay systems (or maybe as importantly, it is not consisten between system types)

Couldn't one replace the delays with 1 and then analyze that system? That shouldn't be more than a line of code if the LFT functionality is used. I think the functionality should be consistent for delay systems, otherwise it's quite weird.

[tallakt]

I was not able to figure out how to convert a DelayLtiSystem to a non-delayed system. I could also implement this if I was helped a bit or perhaps figure it out myself later.

As the phase is only important for the actual bode plot, perhaps a good option to preserve computing speed for other used of the bode function would be to make the initial phase heuristic optional, then use it as default only in the bodeplot function