Skip to content

Discrete to continuous #888

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

Open
wants to merge 2 commits into
base: main
Choose a base branch
from

Conversation

kletze
Copy link

@kletze kletze commented Apr 27, 2023

Added a function to convert discrete time LTI systems to continuous time LTI systems.
The function uses the inverse of the zoh-method. A mathematical derivation can be found for example on this page: https://studywolf.wordpress.com/tag/zero-order-hold/

@murrayrm
Copy link
Member

murrayrm commented Apr 30, 2023

I'm not sure this function is really needed. What's the use case where one wants to convert from continuous time to a (zero-order hold) continuous time model?

If we do include this, it needs to use python-control coding conventions and include unit tests (in the control/tests/ directory).

@sawyerbfuller
Copy link
Contributor

thanks for the contribution!

I have not run into a need for this function but it has been suggested before in #294 , and it is available in Matlab https://www.mathworks.com/help/control/ref/dynamicsystem.d2c.html (but not, as far as I can tell, in scipy).

In addition to adding unit tests, I suggest moving this to being a method of the LTI class in lti.py because it applies equally to both ss and tf systems. and adding a function that calls that method. Please also see the example of lti.sample_system to see how to properly pass other system parameters like signal and system names.

@sawyerbfuller
Copy link
Contributor

I just found that Harold has a cont2discrete method — see _undescretize in https://github.com/ilayn/harold/blob/main/harold/_discrete_funcs.py

might have some good ideas.

@ilayn
Copy link

ilayn commented May 4, 2023

Possibly my most useless spent week on a nerdy subject 😃 It is certainly possible but like Richard, I don't know what value it adds to anything. I just generalized this paper. Not that it makes it more useful but I was neck deep in IQCs so seemed simple enough. Didn't publish it but put it in harold instead. You can do arbitrary LFR connections between s and z in case needed.

The trick is to know how you ended up with that particular discrete model. Otherwise if you go --> zoh this way and tustin <-- this way, things get complicated very quickly.

@tobin
Copy link

tobin commented Nov 26, 2024

Hi folks, I actually ended up here precisely because I just started using the controls system toolbox and was looking for an equivalent of the MATLAB d2c function. I have experimental data in the z-domain, but I like to convert it to the s-domain to think about the continuous time poles and zeroes and to integrate with continuous time models. I'm surprised to see that this is not regarded as useful? Or perhaps I misunderstand the discussion.

@murrayrm
Copy link
Member

Seems like a valid use case (which is what I was asking about). A simple implementation is in this PR, but it needs to follow coding conventions, include unit tests, and be documented (in the user manual).

@kletze Are you available to update? If not, perhaps @tobin can take over this PR?

@ilayn
Copy link

ilayn commented Nov 27, 2024

I have experimental data in the z-domain, but I like to convert it to the s-domain to think about the continuous time poles and zeroes and to integrate with continuous time models.

Note that, there is no such unique mapping, so you have to make some decision on how to map it back from z to s domain, hence the discussion above. If you have experimental data it is better to get the rest of the models to z domain.

But I guess fighting this point forever is not useful in particular for the academic use, so probably there is a need for it.

@robertobucher
Copy link
Member

robertobucher commented Nov 27, 2024

I have an (old) implementation that I use, but also in this case it is important to modify the code to follow your conventions.

def d2c(sys,method='zoh'):
    """
    Continous to discrete conversion with ZOH method

    Call:
    sysc=d2c(sys,method='zoh')

    Parameters
    ----------
    sys :   System in statespace or Tf form 
    method: 'zoh' or 'tustin'

    Returns
    -------
    sysc: continous system ss or tf
    

    """
    flag = 0
    if isinstance(sys, ct.TransferFunction):
        sys = ct.tf2ss(sys)
        flag = 1

    a = sys.A
    b = sys.B
    c = sys.C
    d = sys.D
    Ts = sys.dt
    n = np.shape(a)[0]
    nb = np.shape(b)[1]
    nc = np.shape(c)[0]
    tol=1e-12
    
    if method=='zoh':
        if n==1 and a[0,0]==1:
            A = 0
            B = b/sys.dt
            C = c
            D = d
        else:
            tmp1 = np.hstack((a,b))
            tmp2 = np.hstack((np.zeros((nb, n)), np.eye(nb)))
            tmp = np.vstack((tmp1, tmp2))
            s = la.logm(tmp)
            s = s/Ts
            if la.norm(np. imag(s), np.inf) > np.sqrt(sp.finfo(float).eps):
                print('Warning: accuracy may be poor')
            s = np.real(s)
            A = s[0:n,0:n]
            B = s[0:n,n:n+nb]
            C = c
            D = d
    elif method=='foh':
        a = np.asmatrix(a)
        b = np.asmatrix(b)
        c = np.asmatrix(c)
        d = np.asmatrix(d)
        Id = np.asmatrix(np.eye(n))
        A = la.logm(a)/Ts
        A = np.real(np.around(A,12))
        Amat = np.asmatrix(A)
        B = (a-Id)**(-2)*Amat**2*b*Ts
        B = np.real(np.around(B,12))
        Bmat = np.asmatrix(B)
        C = c
        D = d - C*(Amat**(-2)/Ts*(a-Id)-Amat**(-1))*Bmat
        D = np.real(np.around(D,12))
    elif method=='tustin':
        a = np.asmatrix(a)
        b = np.asmatrix(b)
        c = np.asmatrix(c)
        d = np.asmatrix(d)
        poles = la.eigvals(a)
        if any(abs(poles-1)<200*sp.finfo(float).eps):
            print('d2c: some poles very close to one. May get bad results.')
        
        I = np.asmatrix(np.eye(n,n))
        tk = 2 / np.sqrt (Ts)
        A = (2/Ts)*(a-I)*la.inv(a+I)
        iab = la.inv(I+a)*b
        B = tk*iab
        C = tk*(c*la.inv(I+a))
        D = d- (c*iab)
    else:
        print('Method not supported')
        return
    
    sysc=ct.StateSpace(A,B,C,D)
    if flag==1:
        sysc = ct.ss2tf(sysc)
    return sysc

@robertobucher
Copy link
Member

The function discrete_to_continuous is important because of the design of discrete controllers, but it is important to have other methods as only zoh (in particular in education).

Usually we proposed to our students the following method in order to design a discrete controller:

  1. Convert the continuous plant to a discrete plant
  2. Convert the discrete plant to continuous plant using "tustin"
  3. Design the controller in continuous time domain with the usual methods (Bode, Nyquist etc.) with the new continuous transfer function
  4. Convert the continuous controller to discrete using "tustin"

@murrayrm
Copy link
Member

murrayrm commented Nov 27, 2024

If someone takes a crack at updating this PR to include more methods, in addition to utilizing @robertobucher's code, you might take a look at @ilayn's implementation of undiscretize in harold.

@tobin
Copy link

tobin commented Nov 27, 2024 via email

@kletze
Copy link
Author

kletze commented Nov 28, 2024

My original use case for this PR was to change the sampling time of models estimated in the discrete-time domain from measured values.
Since Matlab provides the d2d method, I figured it would be kind of handy, to implement this functionality here as well.
A look at the GNU Octave source code shows that within the d2d function, it simply first calls the d2c function and then the c2d function with the corresponding sampling time: https://github.com/gnu-octave/pkg-control/blob/main/inst/%40lti/d2d.m.
However, since the pull request for the “d2c” function didn’t really gain broad acceptance and @ilyan had already implemented the function in
https://github.com/ilayn/harold/blob/main/harold/_discrete_funcs.py, and I was deeply involved in other projects, the PR somewhat fell off my radar.
At the time, it was sufficient for my purposes to implement only the ZOH method.
If the PR is updated here, it definately makes sense to implement the other methods as well, as @robertobucher has shown in his code.
I won’t be available for the next two weeks. After that, I could take care of the matter, and if @tobin joins in before then, all the better.

@KybernetikJo
Copy link
Contributor

The bilinear transformation ('tustin') is implemented in slicot, it is part of slycot 0.6.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

7 participants