Skip to content

Improved speed of ctrb and obsv functions #941

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 8 commits into from
Dec 1, 2023
34 changes: 26 additions & 8 deletions control/statefbk.py
Original file line number Diff line number Diff line change
Expand Up @@ -990,13 +990,15 @@ def _control_output(t, states, inputs, params):
return ctrl, closed


def ctrb(A, B):
def ctrb(A, B, t=None):
"""Controllabilty matrix.

Parameters
----------
A, B : array_like or string
Dynamics and input matrix of the system
t : None or integer
maximum time horizon of the controllability matrix, max = A.shape[0]

Returns
-------
Expand All @@ -1016,22 +1018,30 @@ def ctrb(A, B):
amat = _ssmatrix(A)
bmat = _ssmatrix(B)
n = np.shape(amat)[0]
m = np.shape(bmat)[1]

if t is None or t > n:
t = n

# Construct the controllability matrix
ctrb = np.hstack(
[bmat] + [np.linalg.matrix_power(amat, i) @ bmat
for i in range(1, n)])
ctrb = np.zeros((n, t * m))
ctrb[:, :m] = bmat
for k in range(1, t):
ctrb[:, k * m:(k + 1) * m] = np.dot(amat, ctrb[:, (k - 1) * m:k * m])

return _ssmatrix(ctrb)


def obsv(A, C):
def obsv(A, C, t=None):
"""Observability matrix.

Parameters
----------
A, C : array_like or string
Dynamics and output matrix of the system

t : None or integer
maximum time horizon of the controllability matrix, max = A.shape[0]

Returns
-------
O : 2D array (or matrix)
Expand All @@ -1050,10 +1060,18 @@ def obsv(A, C):
amat = _ssmatrix(A)
cmat = _ssmatrix(C)
n = np.shape(amat)[0]
p = np.shape(cmat)[0]

if t is None or t > n:
t = n

# Construct the observability matrix
obsv = np.vstack([cmat] + [cmat @ np.linalg.matrix_power(amat, i)
for i in range(1, n)])
obsv = np.zeros((t * p, n))
obsv[:p, :] = cmat

for k in range(1, t):
obsv[k * p:(k + 1) * p, :] = np.dot(obsv[(k - 1) * p:k * p, :], amat)

return _ssmatrix(obsv)


Expand Down
16 changes: 16 additions & 0 deletions control/tests/statefbk_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,14 @@ def testCtrbMIMO(self):
Wc = ctrb(A, B)
np.testing.assert_array_almost_equal(Wc, Wctrue)

def testCtrbT(self):
A = np.array([[1., 2.], [3., 4.]])
B = np.array([[5., 6.], [7., 8.]])
t = 1
Wctrue = np.array([[5., 6.], [7., 8.]])
Wc = ctrb(A, B, t=t)
np.testing.assert_array_almost_equal(Wc, Wctrue)

def testObsvSISO(self):
A = np.array([[1., 2.], [3., 4.]])
C = np.array([[5., 7.]])
Expand All @@ -62,6 +70,14 @@ def testObsvMIMO(self):
Wotrue = np.array([[5., 6.], [7., 8.], [23., 34.], [31., 46.]])
Wo = obsv(A, C)
np.testing.assert_array_almost_equal(Wo, Wotrue)

def testObsvT(self):
A = np.array([[1., 2.], [3., 4.]])
C = np.array([[5., 6.], [7., 8.]])
t = 1
Wotrue = np.array([[5., 6.], [7., 8.]])
Wo = obsv(A, C, t=t)
np.testing.assert_array_almost_equal(Wo, Wotrue)

def testCtrbObsvDuality(self):
A = np.array([[1.2, -2.3], [3.4, -4.5]])
Expand Down