From 323ae3d45e458ef4520e4ac92ac7bf0207431daa Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Thu, 2 Dec 2021 18:39:31 +0100 Subject: [PATCH] vectorize FRD feedback function --- control/frdata.py | 20 +++++++------------- 1 file changed, 7 insertions(+), 13 deletions(-) diff --git a/control/frdata.py b/control/frdata.py index 5e2f3f2e1..e1ce76b9d 100644 --- a/control/frdata.py +++ b/control/frdata.py @@ -536,20 +536,14 @@ def feedback(self, other=1, sign=-1): if (self.noutputs != other.ninputs or self.ninputs != other.noutputs): raise ValueError( "FRD.feedback, inputs/outputs mismatch") - fresp = empty((self.noutputs, self.ninputs, len(other.omega)), - dtype=complex) - # TODO: vectorize this + # TODO: handle omega re-mapping - # TODO: is there a reason to use linalg.solve instead of linalg.inv? - # https://github.com/python-control/python-control/pull/314#discussion_r294075154 - for k, w in enumerate(other.omega): - fresp[:, :, k] = np.dot( - self.fresp[:, :, k], - linalg.solve( - eye(self.ninputs) - + np.dot(other.fresp[:, :, k], self.fresp[:, :, k]), - eye(self.ninputs)) - ) + # reorder array axes in order to leverage numpy broadcasting + myfresp = np.moveaxis(self.fresp, 2, 0) + otherfresp = np.moveaxis(other.fresp, 2, 0) + I_AB = eye(self.ninputs)[np.newaxis, :, :] + otherfresp @ myfresp + resfresp = (myfresp @ linalg.inv(I_AB)) + fresp = np.moveaxis(resfresp, 0, 2) return FRD(fresp, other.omega, smooth=(self.ifunc is not None))