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

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.