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.