Skip to content

This issue was moved to a discussion.

You can continue the conversation there. Go to discussion →

Scaling for control.impulse_response discrete time is not correct #974

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
cdboschen opened this issue Mar 2, 2024 · 6 comments
Closed

Comments

@cdboschen
Copy link

cdboschen commented Mar 2, 2024

Tested on control 0.9.2 and again on 0.9.4

The resulting impulse response (unit sample response) for the following transfer function (accumulator scaled by 0.1) is not properly scaled:

$$H(z) = \frac{Tz}{z-1}$$

With the sample duration set to $T=0.1$ is being scaled again by $1/T$ such that the result is a step with every output sample set to 1. This is not correct as the result should be 0.1 for every output sample.

To further clarify below is the implementation diagram:

image

If we pass a unit sample into this, it would scale by $T$ to be $0.1$ per the transfer function and block diagram shown and then the output would go to $0.1$ and then repeat for every subsequent sample. I confirmed with Octave that the result is as I expect. Further if I use con.step_response on the same transfer function, the result is also what I expect for that case (the output increments by 0.1).

Code to repeat problem:

import control as con
dt=0.1
accum = con.tf([dt, 0], [1, -1], dt)
t1, y1 = con.impulse_response(accum)     # doesn't scale properly
t2, y2 = con.step_response(accum)           # scales properly

I do see in the docs "For continuous time systems, the initial condition is altered to account for the initial impulse. For discrete-time systems, the impulse is sized so that it has unit area", but would argue that it should not be rescaled especially if Octave and MATLAB are not doing that (I haven't confirmed MATLAB) and should represent the output of the discrete time system described by $H(z)$ alone without modification when a unit sample is applied to the input.

@murrayrm
Copy link
Member

murrayrm commented Mar 2, 2024

There is some discussion of the current approach in #812.

Something to try in python-control, MATLAB, and Octave is the following:

import control as ct
import matplotlib.pyplot as plt
sys = ct.tf(4, [1, 2, 10])  # example from MATLAB documentation
for dt in [0, 0.5, 0.2, 0.05]:
    sys_dt = sys if dt == 0 else ct.c2d(sys, dt)
    impulse_resp = ct.impulse_response(sys_dt)
    plt.plot(impulse_resp.time, impulse_resp.outputs, label=f"{dt=}")
plt.legend()
plt.show()

The output in python-control is
impulse

This shows that the scaling makes sense when you think about the discrete time system being a sampled version of the continuous system (as described in #812).

@cdboschen
Copy link
Author

cdboschen commented Mar 2, 2024

Ah thank you! I was trying to find if this was already posted. I see now from a friend that MATLAB also scales it this way.

Regarding your comment on scaling, yes it totally makes sense in the mapping of an ODE from continuous time IF we were doing that. If I have a discrete time system then there is no further mapping and any connection to some continuous time system should already be in the scaling for that function. My example of an accumulator was exactly that: I included the $T=0.1$ in the transfer function as the "dt" you would see in an integrator (the accumulator is a discrete time integrator or mapping of the continuous time integral to discrete time). Similarly to do the dx/dt for differentiation, you scale by $1/T$. So if I had a continuous time function and wanted to differentiate it in discrete time, then yes I would multiply by $1/T$ as is done. And we note that an impulse response is the integration of the step response. However if I have an arbitrary discrete time system and want it to match some continuous time system at a specific update rate, then the scaling mentioned should already be included in the transfer function. Case in point is the step response: I already included the $T=0.1$ and the step increased by $0.1$ in each step. The function didn't multiply by $0.1$ as would be required to properly represent an accumulator increasing at period $T=0.1$.

If it were to do the scaling as part of the c2d, that would make sense to me. But when I already have a discrete time system that I am modelling with its specific gain constants whatever they may be, it seems better (to me) that it represent what we would expect given the definition of the "unit sample response" as a discrete time impulse response.

(Should I be writing this under #812 and not here?)

@murrayrm
Copy link
Member

murrayrm commented Mar 3, 2024

(Should I be writing this under #812 and not here?)

#812 is merged, so better to keep the comments here.

@ilayn
Copy link

ilayn commented Mar 3, 2024

Hi @cdboschen

This is typically tied to the common problem of defining the dirac function as a standalone unit sample function however that is not the definition of a unit impulse function. Hence, there is no right or wrong, it's just the definition. The discrete version is the Kronecker delta and both have only mathematical meaning under an integral or summation sign because they are technically generalized functions (or tempered distributions and so on depending on who you are hanging out with). They are not functions of time. And they should have the property

$$ \sum_{j}\delta_{ij} a_{j} = a_i $$

Because an impulse does not carry the value of 1 throughout the whole sample but when summed up (again only under summation sign) equals to unity mimicking the dirac delta

$$ \int _{-\infty }^{\infty }f(x),\delta (dx)=f(0) $$

That's why it is not a proper function of $\mathbb{z}$, similar to the dirac delta. Thus you have to scale it with the sampling period to keep that unity behavior. What you want is the unit sample function response and can be done by supplying [1, 0, 0, 0...] which is exactly what it is.

Regarding your comment on scaling, yes it totally makes sense in the mapping of an ODE from continuous time IF we were doing that. If I have a discrete time system then there is no further mapping and any connection to some continuous time system should already be in the scaling for that function.

By mentioning that there is a continuous time system, we don't mean that we always start from a continuous one and derive all discrete-time models from. What I meant over #812, is that there is an "equivalence-class" between such systems that you can argue about an impulse response of all such systems which should result in the same behavior.

And we note that an impulse response is the integration of the step response.

I think you meant the other way around.

In control software generally you don't need to include $T$ in discrete transfer functions. The software takes care of it by the dt keyword. So if you want a discrete backwards euler integrator tf(1, [1, 0], dt=dt) is sufficient and it will arrange things for you (I haven't checked it but should). Hence you don't need to add it yourself or if you want to model everything then you should keep track of $T$ everywhere.

Note that if this is an integrator, then either it is your interpretation since you already have the information that it is linked to continuous time integrator or it is an information that you have beforehand. However suppose someone hands you a Tustin approximation of it $\frac{T}{2} \frac{z+1}{z-1}$ and did not tell you that, then you would not know if it is an integrator or something else. Hence we have to be a bit careful about talking about continuous equivalences and behaviors and etc.

I am not suggesting that your argument is wrong or inaccurate but I am trying to portray how the convention is typically constructed.

@cdboschen
Copy link
Author

cdboschen commented Mar 3, 2024

Thank you @ilayn and I understand your points (and yes, you are correct, I meant the other way around). I really appreciate you taking the time to respond to my previous comments and helping me understand the broader perspective with the detail that you wrote, thank you!

I'm thinking of applications where I am simulating a discrete time system, and for a step and "impulse" response, I really want the result to be for the unit sample function (so the unit sample response specifically), and for my model to match the specific results I would get without "hidden" gain modifications (hence my examples). An "impulse response" does not really exist in discrete time so I would have expected the unit sample response from a con.impulse_response function.

To your point about "interpretation"; I am using the discrete time models as they are without concern of a mapping or equivalence to a continuous time model: my interpretation of the Tustin approximation of an integrator if that was handed to me as $H(z)$ is simply an accumulator with a gain factor followed by a two sample moving average, so to me the gain modifications beyond what I write in the transfer function are a nuisance. (Going back to if I want the true unit sample response for the block diagram as it is, what would the result be?). So my concern isn't how it was mapped but more to say that I'm not concerned with any equivalence to continuous time-- I want the results to be exactly what my system produces exactly at each dt. (These are models of discrete time hardware or software implementations where the gain values are obvious and precise - so would be easy for me to write $H(z)$ with those gains rather than have something derive it for me from dt. I see the value of having that happen in functions that do map from s to z (such as c2d), but less so in functions providing the results from transfer functions that are in discrete time alone. I just wanted that part of my concern to be clear.

So, that's my two cents, but I also get that there must be value in having the dt automatically scale when using the discrete time transfer functions to model continuous time systems. I can get to what I need without too much grief by simply always using dt=1 for my applications and scaling the time axis myself, and perhaps that is the best solution. (And I see that if it was changed to what I really want, then the inclusion of dt is trivial since it would only serve to scale the t that the function returns and have no other effect).

@cdboschen
Copy link
Author

(Should I be writing this under #812 and not here?)

#812 is merged, so better to keep the comments here.

Ok understood I deleted my comment there to clean up.

@python-control python-control locked and limited conversation to collaborators Mar 3, 2024
@murrayrm murrayrm converted this issue into discussion #977 Mar 3, 2024

This issue was moved to a discussion.

You can continue the conversation there. Go to discussion →

Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants