problematic step() output: unstable instead of sinusoidal #384

Activity

1. murrayrm commented on Mar 24, 2020

   

   murrayrm

   on Mar 24, 2020

   Member

   More actions

   The system is ill-conditioned since the closed loop system has a pure imaginary pole => small numerical errors can generate problems.

   The main problem in the python-control version is that it is computing the time range for the step response differently. If you look closely at the units on the x axis for your first plot, you will see they go to 3e15 (!). If you give the step response an explicit time range, you get something that looks like Octave:

   from control.matlab import *
   import numpy as np
   import matplotlib.pyplot as plt
   
   Gp = tf(6, [1, 6, 11, 6])
   Gc = 10
   Gor = Gp*Gc/(1 + Gp*Gc)
   y, t = step(Gor, np.linspace(0, 10, 1000))
   plt.plot(t, y)
   

   

   There is still a problem here to be fixed, which is the computation of the default time for the simulation in cases where there is a pure imaginary pole. Currently this computation is coming from scipy.signal.lti._default_response_times. Related to PR #383 and issue #287.
2. murrayrm

   self-assigned this

   on Mar 24, 2020
3. prashant2 commented on Mar 26, 2020

   

   prashant2

   on Mar 26, 2020

   Author

   More actions

   Thanks. Missed the e15 part.

   As you have mentioned, a problematic default sampling rate gave me the following aliased output even for the simple damped oscillator, G = s/[(s + 0.01)^2 + 1] (this's similar to the 'goto example' mentioned by ilayn in #287, but a bit simpler.)

   a = -0.01
   G1 = tf([1,0],[1,-2*a, 1+a**2])
   y,t = step(G1)
   
   ytrue_func = lambda t : np.exp(a*t)*np.sin(t)
   t_true = linspace(0,700,2000)
   ytrue = ytrue_func(t_true)

   

   Somehow all cases like this seem to be using a fixed 100 time steps, instead of using more to avoid aliasing.
   Scipy output is also identical - does the step function use scipy lti.signal.step?

   The code at https://github.com/ilayn/harold/blob/master/harold/_time_domain.py#L280-L454 mentioned in #287 works very well for this case with a good dt of 0.25 and ~1800 steps.
4. murrayrm commented on Mar 26, 2020

   

   murrayrm

   on Mar 26, 2020

   Member

   More actions

   Thanks for pointing out the Harold code, @prashant2. I'm going to spend some time trying to fix up this problem as part of PR #383.
5. ilayn commented on Mar 26, 2020

   

   ilayn

   on Mar 26, 2020

   Last edited by ilayn

   More actions

   This is what I get from the current master of harold.

   Gp = harold.Transfer(6,[1,6,11,6])
   Gc = 10
   Gor = Gp*Gc/(1 + Gp*Gc)  # should have used feedback() command here anyways
   harold.step_response_plot(Gor)

   

   and for the other example

    a = -0.01
   G1 = harold.Transfer([1,0],[1,-2*a, 1+a**2])
   step_response_plot(G1)

   Note that the dynamics are still decaying until 400s so I think Octave doesn't take that into account (in case you haven't specifically zoomed in).

   

   If this corona WFH thing lets me do some open source work, I'll try to release 1.0.2 soon. Please let me know if you encounter some nonsensical results.

   Richard it would be great if you can look into the repeated pole sensitivity issue. The rest you can copy/paste at will if the code is useful (I can delicense it if needed no attribution etc. required).
6. sawyerbfuller commented on Aug 19, 2020

   

   sawyerbfuller

   on Aug 19, 2020

   Contributor

   More actions

   New code that incorporates @ilayn's code, in PR #454 seems to fix this issue:

   

   as well as the second example.

   

   Should be resolved when that PR is merged.
7. bnavigator

   linked a pull request that will close this issueimproved default time vector for time response functions that takes zeros into account. #454on Aug 19, 2020
8. murrayrm

   closed this as completedin #454on Oct 17, 2020