The question about time responses and time series data is a good one. I've been thinking about it too. We don't want to make Python-Control dependent on Pandas but I would argue that the Pandas Series and DataFrame objects are perfectly-designed tools for time series data storage, plotting and analysis (I use them all the time when working with industrial time series data).

If we can find a way to seamlessly use Python-Control and Pandas 'hand-in-glove' so-to-speak, that would be very powerful I think. (This could also save us writing a raft of tools to do all that).

Thanks, I agree, we need to have a good justification to create any new data objects. One thing with two functions, you need to transfer the data from the calculation function to the plotting function. If it's just a couple of things, like 2 arrays, then I think a tuple or dict is fine and it doesn't merit creating an object.

A benefit of an object would be if there is additional 'meta data' that could be automatically passed to the plotting method, such as axis and legend labels. As an example, df.plot() (with no additional arguments), where df is a Pandas dataframe, makes a complete plot with multiple coloured lines, a legend with the names of each series and an x-axis label all from the metadata in the dataframe object.

The case for a new data object gets stronger if it also has additional uses in its own right (i.e. for things other than simply making a plot). I was also wondering, if we make one FreqResp object, can it be used to make a Bode, Nyquist, and Nichols chart? Is that possible? Then, does the FreqResp object have any other uses beyond making plots? In other words, could it be useful in other analysis workflows?

So, in summary I think we need to really challenge the need for any new object and only do it if it has enough benefits over a simpler solution (like two separate functions).

The case for a new data object gets stronger if it also has additional uses in its own right (i.e. for things other than simply making a plot). I was also wondering, if we make one FreqResp object, can it be used to make a Bode, Nyquist, and Nichols chart? Is that possible? Then, does the FreqResp object have any other uses beyond making plots? In other words, could it be useful in other analysis workflows?

I think there already is an object that can store frequency_response data:frd.FrequencyResponseData. It can also be used for system-id applications. (By the way, freqresp is now deprecated and only part of matlab layer).

I think there is a case to be made that a TimeResponseData (Or TimeSeriesData) object, with an optional input field as suggested by @murrayrm could have a similar use.

Then this proposal could include simply adding a plot method to these, and maybe adding one more object for root_locus?

It also says this at the bottom:

I will model the freqresp return structure similarly in a new PR when this has been merged

Had a quick look at the Julia proposal. Looks like functions like lsim, step, etc. will in future return a SimResult(y, t, x, uout, sys) struct (Julia doesn't have object-oriented programming, i.e. classes) and this struct can then be used as a single argument in a new method of the plot function, so you can do:

sim_result = lsim(my_sys, u, t)
plot(sim_result)

or simply:

plot(lsim(my_sys, u, t))

For comparison, here is what this could look like in Python Control:

sim_result = input_output_response(my_sys, u, t)
sim_result.plot()

or

input_output_response(my_sys, u, t).plot()

Perhaps, we could preserve the old usage of returning a tuple for the matlab-like functions lsim, step, etc:

yout, t, xout = lsim(sys, u, t, x0)

I've been playing around with some ideas for time responses and I have an implementation that allows the following to work:

response = ct.input_output_response(sys, t, u)
response.plot()

or you can say

t, y = ct.input_output_response(sys, t, u)
plt.plot(t, y)

This is all done by defining a new class InputOutputResponse and then changing the time response functions to return an element of that class. By implementing an __iter__ method in the class, you can allow the object to be "assigned" to a tuple.

So far, everything is backward compatible so that all unit tests run without change. But we get the advantage of having an object returned to which we can add additional methods (like plot).

I'll do a bit more work on it and post a draft PR so that people can take a look. I think it fits the intent of the paradigm here.

For the input-output response, it would be nice (IMO) if response.plot() produced a figure with two subplots: upper plot of (t, y), lower plot of (t, u), with appropriate titles ('Outputs', 'Inputs') and legends ('y_1', 'y_2', ...).

Thanks for the feedback @sawyerbfuller.

an alternative to freqresp(sys, w) is the more pythonic mag, phase, w = sys.frequency_response(w)

I agree that providing the most commonly-used functions as methods of sys is a good idea. The question is where to draw the line. Here are some other things that could be done (which doesn't mean they should be!):

• sys.frequency_response(w)
• sys.margin()
• sys.pzmap()
• sys.root_locus()
• sys.bode_plot()
• sys.nyquist_plot()
• ...etc.

Or we could have some kind of catch-all plot method (similar to Pandas):

• sys.plot(kind="nyquist")

or

• sys.plot.nyquist()

I don't plan to do anything with the sys object at the moment unless others think it's important. We can always add these later once the basic plotting functions are sorted out.

I think there might now be a singular_values_plot Function.

I see that now. I'll add this to the table above. Thanks.

@billtubbs i am also in favor of everything suggested by @murrayrm above.

In regards to the frd Class, it would be nice if we could re-use it to hold data coming out of eg frequency_response, but it does a lot of other things eg can be multiplied by a transfer function. Maybe it would become ungainly?

@billtubbs in regards to frequency_response, my recollection was that there was consensus that it was an accident that it outputs mag, phase, w. It should change to output response as complex number. If you want mag, phase, w, maybe you could use the new plot-free bode command.

In regards to bode vs. nyquist plots, it turns out they both kind of want similar frequency points( eg both need lots of frequency points around resonant peaks) so it might be possible to do as you suggest and use the same data for each. This is how it works now.

For the return value for sys.frequency_response, you can set things up so that this is a FrequencyResponseData object and then set up an __iter__ method that returns either mag, phase, omega (default) or s, omega. Similarly, you can use the property decorator to allow you to have attributes that convert data types. So something like the following should be possible:

fr = sys.frequency_response()
d, w = fr.d, fr.w                         # d is complex, w is real
mag, phase, w = fr                        # implement via __iter__
mag, phase, w = fr.mag, fr.phase, fr.w    # implement via @property
d, w = fr(return_complex=True)            # implement via __iter__
fr.plot(kind='bode')
ct.bode_plot(sys)                         # same plot
ct.bode_plot(fr)                          # same plot
fr.plot(kind='nyquist')
ct.nyquist_plot(sys)                      # same plot
ct.nyquist_plot(fr)                       # same plot

Doing things this way will make everything backward compatible with the current code but also alow us to be consistent with the new plotting paradigm. (This is also compatible with the way that TimeResponseData works and you can look there for an example of how to implement.)

Regarding Bode vs Nyquist: there are two differences that I can think of that will need some thought:

• For Nyquist plots, we want to capture the value of the transfer function at (and near) s = 0, but for Bode plots this point is not shown, so should be omitted when plotting.
• For Nyquist plots, we add indentations around poles on the imaginary axis. These aren't really part of the frequency response, and so we would need to think through how they would be represented in the FrequencyResponseData object (perhaps a separate attribute?).

Finally, the current code for nyquist_plot returns the number of encirclements. If we change this to return handles to the figure and axis, we should think about how/where to return the number of encirclements.

Thanks for the feedback. This looks like a good path forward. Quick question:

d, w = fr(return_complex=True)            # implement via __iter__

Why not

d, w = sys.frequency_response(return_complex=True)
d, w = frequency_response(syslist, return_complex=True)

Oh, I see that TimeResponseData uses this approach with return_x for backwards compatibility. Okay, not a bad idea.

I agree that those variants should also work. The main thing is that the __iter__ method needs to know how many (and what type) of values to return. So we can keep the default as mag, phase, w for backward compatibility, but implement other options as well (with keywords).

They allow bode, nyquist and nichols plots on FRD sys objects. I tested them and they look okay with the example data above.

For a system with imaginary poles I'm not sure. I tried this:

>> sys = tf(1, [1 2*0.2 1]);
>> [re,im,wout] = nyquist(sys);
>> fr_data = squeeze(re) + squeeze(im)*i;
>> fr_sys = frd(fr_data, wout);
>> subplot(1, 2, 1)
>> nyquist(sys)
>> subplot(1, 2, 2)
>> nyquist(fr_sys)

and the two plots are identical. But when I generated some random FRD data from this system and tried to make the nyquist plot it had to think for a very long time and the resulting plot looked odd. Is this what you meant?

>> fr_data2 = squeeze(re) + 0.01*randn(size(re)) + (squeeze(im) + 0.01*randn(size(im)))*i;
>> fr_sys2 = frd(fr_data2, wout);

Oh wait, you said pure imaginary. Sorry.
I tested this system:

>> sys = tf(1, conv([-0.5i 1],[0.5i 1]));
>> figure(1)
>> subplot(1,2,1)
>> nyquist(sys)
>> [re,im,wout] = nyquist(sys);
>> fr_data = squeeze(re) + squeeze(im)*i;
>> fr_sys = frd(fr_data, wout);
>> subplot(1,2,2)
>> nyquist(fr_sys)

and get the following:

Ah, note the scaling on the axes above (10^-8, 10^7). When you re-scale the MATLAB nyquist plot it looks like this:

Here is Python control Nyquist plot with same axis scaling:

What is going on here?

That's the infinite circle feature (#534), see also #671

Ah right. Neat. To represent "to infinity and back" is it?

How about this solution:

bode_mag_plot(syslist)  # quick plot
bode_mag_plot(syslist, ax=ax)  # adding to existing figure
ax = bode_mag_plot(syslist)  # returning axis handle

and same for bode_phase_plot()

And then in the new paradigm:

fr = frequency_response(syslist)
# quick plot
fr.plot(kind='bode')
fr.plot(kind='bode_mag')
fr.plot(kind='bode_phase')
# or
fig, axes = fr.plot(kind='bode')
fr.plot(kind='bode_mag', ax=ax)
ax = fr.plot(kind='bode_phase')
... etc

Some ideas:

• count[, contour] = nyquist_analysis(sys, return_contour=True/False)
• count[, contour] = nyquist_criterion(sys, return_contour=True/False)
• contour[, count] = nyquist_curve(sys, return_count=True/False)

Also, we should keep in mind that for a MIMO system, the Nyquist criterion is different (one looks at net encirclements of the origin for $\det(I + L(s))$. This isn't currently implemented, but we should make sure that we don't get in the way of that.

How do I post as a draft PR? Just give you a link to my fork? All the unit tests have been updated so it's a working branch.

Yes I need to write the documentation. And I am working on a Jupyter notebook as a tutorial to all the plotting commands.

Like a normal PR, but there will be a point in the process where you can mark it as draft. Once it is there, everyone can see and comment, but can't merge. As you update your branch on GitHub, the PR gets updated.

Sounds good. I should probably merge my branch with my master first I guess, and then submit my master for the draft pull request.

It's usually better to work in a branch and then submit the PR from the branch (in your fork). Just make sure to build (or rebase) on top of the latest (upstream) master.

For this case (singular_values_plot) I think it can probably be removed. But for something like a Bode plot, there are two axes => you need labels (I think) in order to get the magnitude and phase axes (subplots).

More generally, think the question is what we want the behavior to be if you do something like this:

ct.bode_plot(sys1)
ct.bode_plot(sys2)

versus something like this:

ct.bode_plot(sys1)
ct.nyquist_plot(sys1)

I think the behavior in the first case should be to plot the second system on the same axes as the first. But for the second case it seems like we should get a new set of axes (probably on the same figure, clearing the old one?).

I recommend we follow the conventions of Matplotlib and Pandas and not try to make our plot functions too 'clever'.

The convention when plotting in pyplot mode is that it will always try to add the new plot to whatever is in the current axis, regardless of what that is or if it makes sense. The only time a plot function creates a new figure is when none exists.

I checked in Pandas and it does indeed try to add anything to the current plot, even when this makes no sense:

This is in fact what happens in my current branch:

When something like this happens, the user should think "oh, I should have opened a new figure before making the second plot" which should lead them to the solution (which is more explicit):

sys = tf(1, [1, 2, 1])
ct.bode_plot(sys)
plt.figure()
ct.nyquist_plot(sys)

I get why it was being handled differently before (with a hidden labelling system) to replicate MATLAB behaviour but to me that's not a fool-proof method and it's a bit confusing for the user who may not know what determines whether a plot is added to an existing axes or when it is erased first or when a new figure is created.

The nice thing about the changes we're introducing now is that if they want, the user can take full control of where the plot goes using the new arguments and returned objects:

fig, axes = plt.subplots(2, 2)
ct.bode_plot(sys, axes=axes[:, 0])
ct.nyquist_plot(sys, ax=axes[0, 1])
ct.step_response(sys, ax=axes[1, 1])
plt.show()

But I'm happy to be out-voted on this one and re-introduce the 'smart' replacement of any existing plots/lines.

B.t.w. I can't see any reason why you would want to clear any existing content in an axis or figure and replace it. In matplotlib, plots aren't shown until the end (unlike Matlab), so the content that was wiped would never be seen. After a plot is rendered, I don't think you can change it can you?

Seems fine to take the matplotlib/pandas approach. But I do think that most people would expect that calling bode_plot twice in a row would do the usual thing of overlapping the plots without extra arguments. So

bode_plot(sys1)
bode_plot(sys2)

rather than

fig, ax = bode_plot(sys1)
bode_plot(sys2, ax=ax)

(though the latter should also work, of course).

Happy to be overridden and you are writing the code, so you get to choose -:).

Yes, no problem with that. That is what it does (if no axes argument is given, it gets the current figure and adds the bode plots to the first two axes).

The only debate was what to do if someone tries to add a nyquist plot to a bode plot, or some other strange combination.

The link points to the MPL 2.0.2 documentation which is a bit outdated. The closest I could find in the current doc is Usage Patterns.

I totally agree with the conclusion: We should do it the same way as most of the matplotlib.pyplot functions and matplotlib.axes.Axes methods do it: Return Artists collections.

Hi @roryyorke, I'm currently working on this. If you can go through the table above which is a list of the proposed changes we're making (also compares to the latest Julia ControlSystems equivalents) and see if you think these make sense. As well as having every plot function return a fig and/or axes, we have to come up with alternative ways to calculate the data that these functions previously returned. This is perhaps the trickier part of this proposal. (Although we will always have the original MATLAB equivalents which do plotting and return data in one function).

HI @billtubbs, thanks for working on this.

The changes on the whole make sense; splitting computation from presentation is a good idea.

I'm not sure about using custom classes to represent results; it does look tidy, and as you mentioned it could be handy if one needs metadata.

I love using pandas, but it's a fairly heavy dependency. A pandas-compatibility module, separate from plotting, could be useful, to generate DataFrames from time series and frequency response data.

FWIW, one of my favourite Matlab plotting functions, bodemag, is missing; however, we could add it later (it's hardly complicated if one has bode).

To expand on my suggestion of returning a list of Artists instead of the figure or axes, using this would look something like:

h1 = nyquist_plot(sys1)
ax = h1[0].axes
h2 = nyquist_plot(sys2, ax=ax)
h1[0].set_label('System 1')
h2[0].set(label='System 2', linewidth=5, alpha=0.5, linestyle='--')
ax.legend()

The equivalent using a fig&axes return isn't any more complicated (see below), but I don't know if there's any guarantee on the ordering of Artists in the list returned by Axes.get_children; I googled but could find no discussion of this.

plt.close('all')
fig, ax = nyquist_plot2(sys1)
nyquist_plot2(sys2, ax=ax)
h1, h2 = ax.get_children()[:2]
h1.set_label('System 1')
h2.set(label='System 2', linewidth=5, alpha=0.5, linestyle='--')
ax.legend()

@roryyorke I agree with your point about bode_mag. In this new paradigm we should not loose any functionality and I think we should keep the underlying functional approach. Here's one way we could achieve the bode_mag plot functionality in the new paradigm:

fr = freqresp(sys, w)

plt.figure(1)
fr.plot()  # make a Bode plot (by default, kind='bode'), returns fig

plt.figure(2)
fr.plot(kind='mag')  # make the mag plot, returns ax

plt.figure(3)
fr.plot(kind='phase')  # make the phase plot, returns ax

However, we could retain:

plt.figure(1)
ct.bode_plot(sys)  # returns (fig, axes) not (mag, phase, w) as currently

plt.figure(2)
ct.mag_plot(sys)  # returns ax

plt.figure(3)
ct.phase_plot(sys)  # returns ax

Actually, I found more work I did dating to Apr 30, 2022 titled 'work in progress' which I have now merged into the master of this fork:
https://github.com/basicmachines/python-control

I added some temporary files that I was using to guide the code development that might be useful to understand what I was doing:

• plotting_readme.md
• Tests-of-Plotting-Tool-Modifications.ipynb

B.t.w. I tried to contact someone on the Julia ControlSystems project (@murrayrm you gave me the name) but they never replied. I was hoping to collaborate and co-ordinate the changes a bit—see table above. I think this would still be worthwhile as some consistency between the projects would be good (but not essential).

I think it's only the last 6 rows of that table (in plotting_readme.md) that I didn't complete. Although there will be more work to do to integrate the stale fork now.

I recall one of the biggest challenges I encountered but I can't remember if I overcame it or not:

• In the Pandas plotting paradigm, the user can create and prepare a figure and then individual subplots may be added by the user using a keyword (ax=..). However, all Pandas plots use a cartesian grid system (i.e. axis projection). The challenge we have is that in some cases we need a polar grid. Changing an existing axis that the user has created from cartesian to polar is difficult. But I think @bnavigator found a solution and I can't remember if I got it working, sorry!

In matplotlib, you're supposed to specify the projection type when you create the figure. E.g.:

fig, ax = plt.subplots(subplot_kw=dict(projection='polar'))

But if you look in the notebook linked above, I have some code as follows, which suggests I did get it working (almost, look at the plots, there are still some bugs) using a position keyword argument or in some cases the desired ax keyword:

# Create figure with subplots
fig = plt.figure(figsize=(9,8))
sgrid(position=(2, 2, 1))
ax2 = fig.add_subplot(2, 2, 2)
zgrid(ax=ax2)  # alternative for zgrids
#zgrid(position=(2, 2, 2))
sgrid(position=(2, 2, 3))
ax4 = fig.add_subplot(2, 2, 4)
zgrid(ax=ax4)  # alternative for zgrids
plt.tight_layout()
savefig('grids_subplots')

This is probably where I got bogged down and should have called for help!

Another question: If you define a system with named input and output signals, do the names appear on the plot axes?

Yup.

sys = ct.rss(4, inputs=['velocity', 'steering'], outputs=['xpos', 'ypos'], name='car'])
ct.step_response(sys).plot()

Final thing for me, would be to pass unrecognized keyword arguments on to the plot function.
E.g.

sys = ct.rss(4, 2, 2)
ct.step_response(sys).plot(linewidth=2, linestyle='--')

Good point. I've got that on my list to implement prior to the non-draft PR, but haven't done anything on it yet. One question is going to be how to handle multi-signal and/or multi-trace plots. So in your example above, I would assume that all lines are width 2 and dashed. But what if I want to combine outputs onto a single axes (using combine_outputs=True) and then specify different colors/width/styles for each? Perhaps

ct.step_response(sys).plot(
    output_color=['k', 'b'], trace_linestyle=['-', --'], combine_signals=True, combine_traces=True)

which would produce something like this (created by updating the linestyles manually, using the returned lines):

(the trailing 'sys[65]' in the legend will go away when implemented).

See PR #920 for the changes that implement the time response functions above.

Moving on to frequency response functions (Bode, Nyquist) next.