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.