solve bandwidth by bisection for zero-crossing

SCLiao47 · SCLiao47 · commit bc2ddff63716 · 2023-05-01T00:32:56.000-04:00
diff --git a/control/lti.py b/control/lti.py
@@ -20,7 +20,7 @@
 from .namedio import NamedIOSystem, isdtime
 
 __all__ = ['poles', 'zeros', 'damp', 'evalfr', 'frequency_response',
-           'freqresp', 'dcgain', 'pole', 'zero']
+           'freqresp', 'dcgain', 'bandwidth', 'pole', 'zero']
 
 
 class LTI(NamedIOSystem):
@@ -202,6 +202,39 @@ def _dcgain(self, warn_infinite):
         else:
             return zeroresp
 
+    def bandwidth(self, dbdrop=-3):
+        """Return the bandwidth"""
+        raise NotImplementedError("bandwidth not implemented for %s objects" %
+                                  str(self.__class__))
+
+    def _bandwidth(self, dbdrop=-3):
+        # check if system is SISO and dbdrop is a negative scalar
+        if (not self.issiso()) and (dbdrop >= 0):
+            raise ValueError("NOT sure what to raise #TODO ")
+
+        # result = scipy.optimize.root(lambda w: np.abs(self(w*1j)) - np.abs(self.dcgain())*10**(dbdrop/20), x0=1)
+        # # this will probabily fail if there is a resonant frequency larger than the bandwidth, the initial guess can be around that peak
+
+        # use bodeplot to identify the 0-crossing bracket
+        from control.freqplot import _default_frequency_range
+        omega = _default_frequency_range(self)
+        mag, phase, omega = self.frequency_response(omega)
+
+        dcgain = self.dcgain()
+        idx_out = np.nonzero(mag - dcgain*10**(dbdrop/20) < 0)[0][0]
+
+        # solve for the bandwidth, use scipy.optimize.root_scalar() to solve using bisection
+        import scipy
+        result = scipy.optimize.root_scalar(lambda w: np.abs(self(w*1j)) - np.abs(dcgain)*10**(dbdrop/20), 
+                                            bracket=[omega[idx_out-1], omega[idx_out]],
+                                            method='bisect')
+
+        # check solution
+        if result.converged:
+            return np.abs(result.root)
+        else:
+            raise Exception(result.message)
+
     def ispassive(self):
         # importing here prevents circular dependancy
         from control.passivity import ispassive
@@ -499,6 +532,33 @@ def dcgain(sys):
     return sys.dcgain()
 
 
+def bandwidth(sys, dbdrop=-3):
+    """Return the first freqency where the gain drop by dbdrop of the system.
+
+    Parameters
+    ----------
+    sys: StateSpace or TransferFunction
+        Linear system
+    dbdrop : float, optional
+        By how much the gain drop in dB (default = -3) that defines the
+        bandwidth. Should be a negative scalar
+
+    Returns
+    -------
+    bandwidth : #TODO data-type
+        The first frequency where the gain drops below dbdrop of the dc gain
+        of the system.
+
+    Example
+    -------
+    >>> G = ct.tf([1], [1, 2])
+    >>> ct.bandwidth(G)
+    0.9976
+
+    """
+    return sys.bandwidth(dbdrop)
+
+
 # Process frequency responses in a uniform way
 def _process_frequency_response(sys, omega, out, squeeze=None):
     # Set value of squeeze argument if not set
diff --git a/control/statesp.py b/control/statesp.py
@@ -1424,6 +1424,10 @@ def dcgain(self, warn_infinite=False):
         """
         return self._dcgain(warn_infinite)
 
+    def bandwidth(self, dbdrop=-3):
+        """Return the bandwith"""
+        return self._bandwidth(dbdrop)
+
     def dynamics(self, t, x, u=None, params=None):
         """Compute the dynamics of the system
 
diff --git a/control/xferfcn.py b/control/xferfcn.py
@@ -1246,6 +1246,10 @@ def dcgain(self, warn_infinite=False):
 
         """
         return self._dcgain(warn_infinite)
+    
+    def bandwidth(self, dbdrop=-3):
+        """Return the bandwith"""
+        return self._bandwidth(dbdrop)
 
     def _isstatic(self):
         """returns True if and only if all of the numerator and denominator
