Observability Gramian for discrete-time systems

[billtubbs]

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

python-control/control/statefbk.py

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.

[billtubbs]

Expected output should be this I think:

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

[billtubbs]

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.

[bnavigator]

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.

[billtubbs]

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?

[billtubbs]

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

[billtubbs]

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.

[billtubbs]

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.

[murrayrm]

Fixed in #969.