Skip to content
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 13 additions & 13 deletions control/modelsimp.py
Original file line number Diff line number Diff line change
Expand Up @@ -377,14 +377,14 @@ def minimal_realization(sys, tol=None, verbose=True):
def _block_hankel(Y, m, n):
"""Create a block Hankel matrix from impulse response."""
q, p, _ = Y.shape
YY = Y.transpose(0,2,1) # transpose for reshape
YY = Y.transpose(0, 2, 1) # transpose for reshape

H = np.zeros((q*m,p*n))
H = np.zeros((q*m, p*n))

for r in range(m):
# shift and add row to Hankel matrix
new_row = YY[:,r:r+n,:]
H[q*r:q*(r+1),:] = new_row.reshape((q,p*n))
new_row = YY[:, r:r+n, :]
H[q*r:q*(r+1), :] = new_row.reshape((q, p*n))

return H

Expand Down Expand Up @@ -480,22 +480,22 @@ def eigensys_realization(arg, r, m=None, n=None, dt=True, transpose=False):
if (l-1) < m+n:
raise ValueError("not enough data for requested number of parameters")

H = _block_hankel(YY[:,:,1:], m, n+1) # Hankel matrix (q*m, p*(n+1))
Hf = H[:,:-p] # first p*n columns of H
Hl = H[:,p:] # last p*n columns of H
H = _block_hankel(YY[:, :, 1:], m, n+1) # Hankel matrix (q*m, p*(n+1))
Hf = H[:, :-p] # first p*n columns of H
Hl = H[:, p:] # last p*n columns of H

U,S,Vh = np.linalg.svd(Hf, True)
Ur =U[:,0:r]
Vhr =Vh[0:r,:]
Ur =U[:, 0:r]
Vhr =Vh[0:r, :]

# balanced realizations
Sigma_inv = np.diag(1./np.sqrt(S[0:r]))
Ar = Sigma_inv @ Ur.T @ Hl @ Vhr.T @ Sigma_inv
Br = Sigma_inv @ Ur.T @ Hf[:,0:p]*dt # dt scaling for unit-area impulse
Cr = Hf[0:q,:] @ Vhr.T @ Sigma_inv
Dr = YY[:,:,0]
Br = Sigma_inv @ Ur.T @ Hf[:, 0:p]*dt # dt scaling for unit-area impulse
Cr = Hf[0:q, :] @ Vhr.T @ Sigma_inv
Dr = YY[:, :, 0]

return StateSpace(Ar,Br,Cr,Dr,dt), S
return StateSpace(Ar, Br, Cr, Dr, dt), S


def markov(*args, m=None, transpose=False, dt=None, truncate=False):
Expand Down
Loading