Observability Gramian for discrete-time systems

I noticed in the source code for the gram function that it says only continuous-time systems are supported for now:

Lines 1122 to 1127 in 85328d9

    
           # TODO: Check for continuous or discrete, only continuous supported for now 
        
               # if isCont(): 
        
               #    dico = 'C' 
        
               # elif isDisc(): 
        
               #    dico = 'D' 
        
               # else:

I checked and it doesn't raise an exception when you pass a discrete time system but I don't think this means it is calculating a gramian for a discrete-time system:

A = [[-0.31,  0.21], [-0.68, -0.57]]
B = [[1.23], [1.42]]
C = [[ 1.32, -0.55]]
D = [[0.]]
sysd = ct.ss(A, B, C, D, dt=1)
assert sysd.isdtime()
Wo = ct.gram(sysd, 'o')
print(Wo)

array([[ 3.24632551, -0.19876604],
       [-0.19876604,  0.19212128]])

What would it take to get it working for discrete time systems? Is it simply a different call to the same slycot module, sb03md, or is it more involved?

Since I need it I'd be happy to try implementing it.

Expected output should be this I think:

array([[ 2.06193690, -0.57482614],
       [-0.57482614,  0.78661650]])

Where is the documentation on slycot? I tried googling and found the API but couldn't find anything specific on sb03md.

From wikipedia, the two equations that need to be solved are quite similar:

Observability Gramian for a continuous-time system:
$$A^T W_o + W_o A = -C^T C$$

Observability Gramian for a discrete-time system:
$$A^T W_{o} A - W_{o} = -C^T C$$

So I guess I just need to figure out how to do the second one.

Any advice/pointers would be appreciated.

There is no hosted doc of Slycot (yet) (python-control/Slycot#202).

You can look at the documentation of SLICOT's SB03MD
There are the docstrings, which are used by IDEs like spyder to render documentation on-the-fly.

You are probably looking for the dico='D' parameter of slycot.sb03md57.

Okay well that was easy. If I simply allow dico='D' option it calculates the correct Gramians.

The code is all there, just commented-out:

    # TODO: Check for continuous or discrete, only continuous supported for now
        # if isCont():
        #    dico = 'C'
        # elif isDisc():
        #    dico = 'D'
        # else:

I simply replaced above with:

    if sys.isctime():
        dico = 'C'
    elif sys.isdtime():
        dico = 'D'
    else:
        raise ValueError("sys")

Is there a reason this wasn't implemented? Maybe at the time of writing isCont and isDisc didn't exist and the author was planning to come back and finish it but didn't perhaps?

Unless anyone knows of a problem with this, I'll just go ahead and implement it, including updated unit tests.

I noticed this check for stability in the continuous-time case:

    # Check for continuous or discrete
    if sys.isctime():
        dico = 'C'

        # TODO: Check system is stable, perhaps a utility in ctrlutil.py
        # or a method of the StateSpace class?
        if np.any(np.linalg.eigvals(sys.A).real >= 0.0):
            raise ValueError("Oops, the system is unstable!")

Is this the right way to do this for discrete-time systems?

    elif sys.isdtime():
        dico = 'D'

        if np.any(np.abs(sys.poles()) >= 1.):
            raise ValueError("Oops, the system is unstable!")

Correct me if I'm wrong but we don't have an sys.isstable method yet.

I've submitted the pull request.

One other thing I noticed, the gram function has an argument type which is a reserved Python keyword.

def gram(sys, type):

Should we change this to something else?

Or we could replace gram with obsv_gram and ctrl_gram.

Fixed in #969.

Navigation Menu

Repository navigation

Observability Gramian for discrete-time systems #967

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

	# TODO: Check for continuous or discrete, only continuous supported for now
	# if isCont():
	# dico = 'C'
	# elif isDisc():
	# dico = 'D'
	# else:

Observability Gramian for discrete-time systems #967

Description

Activity

billtubbs commented on Feb 8, 2024

billtubbs commented on Feb 10, 2024

bnavigator commented on Feb 11, 2024

billtubbs commented on Feb 11, 2024

billtubbs commented on Feb 11, 2024

billtubbs commented on Feb 11, 2024

billtubbs commented on Feb 11, 2024

murrayrm commented on Mar 30, 2024

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions