Describing functions #521

murrayrm · 2021-01-28T20:39:03Z

This PR adds a basic describing function capability, mainly captured in the new describing_function_plot command. Almost all of the new code is in the descfcn.py file, with a few changes elsewhere to support some of the patterns. There is a Jupyter notebook demonstrating the use of the code.

Additional changes:

Added a mirror keyword to nyquist_plot that sets the line style used for plotting the complementary portion of the Nyquist plot. This defaults to --, which matches the FBS convention.
Added an _isstatic method to the NonlinearIOSystem class that allows passing a nonlinear I/O object as the argument for computing describing functions.
Added a __call__ method to the NonlinearIOSystem class that direct calling of a static nonlinear input/output function (returning the value of the function at the given input).

coveralls · 2021-01-28T20:46:30Z

Coverage increased (+0.3%) to 88.024% when pulling 93d97cc on murrayrm:describing_functions into d75c395 on python-control:master.

…be unneccesary and to avoid a merge conbflict with python-control#521

sawyerbfuller · 2021-01-28T23:23:00Z

control/iosys.py

@@ -450,6 +450,10 @@ def issiso(self):
        """Check to see if a system is single input, single output"""
        return self.ninputs == 1 and self.noutputs == 1

+    def isstatic(self):


statespace and xferfcn have methods is_static_gain() I added. should we choose a standard conventiuon? I would lean toward isstatic() because it is shorter and matches issiso(),

I agree we should be uniform. I'm not sure that is_static_gain() would make sense for a nonlinear system, so that is another reason to opt for isstatic().

Also: these are mainly internally used functions, so I think it is OK to stick with isstatic/issiso rather than is_static/is_siso.

Converted to _isstatic in latest push since this really is an internal function, not intended to be seen by the user.

@sawyerbfuller: when we update is_static_gain, we can either change to _isstatic if all usage is internal or change to isstatic (and update descfcn.py for consistency).

renamed in #524

roryyorke

I'll try to review from describing_function_plot on later this weekend.

roryyorke · 2021-01-30T04:09:02Z

control/descfcn.py

+
+# Class for nonlinearities with a built-in describing function
+class DescribingFunctionNonlinearity():
+    """Base class for nonlinear functions with a describing function


Suggested change

"""Base class for nonlinear functions with a describing function

"""Base class for nonlinear systems with a describing function

roryyorke · 2021-01-30T04:09:53Z

control/descfcn.py

+    """Base class for nonlinear functions with a describing function
+
+    This class is intended to be used as a base class for nonlinear functions
+    that have a analytically defined describing function (accessed via the


Suggested change

that have a analytically defined describing function (accessed via the

that have an analytically defined describing function (accessed via the

roryyorke · 2021-01-30T04:40:44Z

control/descfcn.py

+#
+
+"""The :mod:~control.descfcn` module contains function for performing
+closed loop analysis of systems with static nonlinearities using


I find the static terminology a little confusing, especially given that some of the describing functions declared here, e.g. hysteresis, return False in their _isstatic methods. Maybe "non-linearities with frequency-independent describing functions", though that is rather long-winded.

roryyorke · 2021-01-30T04:45:38Z

control/descfcn.py

+    implement a `call` method that evaluates the nonlinearity at a given point
+    and an `_isstatic` method that is `True` if the nonlinearity has no
+    internal state.
+


How about: "Subclasses should override the describing_function, __call__, and _isstatic methods."

It looks like __call__ 's inputs and outputs must be scalar, since __call__ is allowed to change instance memory, so the sequence of calls matters.

Can describing_function be called with a numpy.ndarray of amplitude values? This seems like a reasonable thing to want to do.

If I understand correctly, if __call__ updates any instance state, then _isstatic must return False -- it would be good to mention that too.

I've updated the documentation a bit to address the comments above.

In terms of calling DescribingFunctionNonlinearity.description_function with an array: the intent is that this is a low level function that just evaluates the describing function at a scalar value. The describing_function function has the code to evaluate on an array.

roryyorke · 2021-01-30T04:49:19Z

control/descfcn.py

+
+    """
+    # If there is an analytical solution, trying using that first
+    if try_method and hasattr(F, 'describing_function'):


the hasattr check will succeed for any base class of DescribingFunctionNonlinearity - shouldn't you also catch the NotImplementedError ?

Good catch. Fixed.

roryyorke · 2021-01-30T05:15:46Z

control/descfcn.py

+    #
+
+    # Evaluate over a full range of angles
+    theta = np.linspace(0, 2*np.pi, num_points)


this is not something I'm not at all sure about: should this omit the final point, so it's more like a discrete Fourier transform?

Good catch. We are approximating the integral by a discrete sum with the assumption that the function is constant on each small interval. But this means that we don't want to include the value at 2pi (since we already captured the segment just before 2pi in the previous evaluation).

roryyorke · 2021-01-30T05:17:21Z

control/descfcn.py

+        # Save the scaling factor for to make the formulas simpler
+        scale = dtheta / np.pi / a
+
+        # Evaluate the function (twice) along a sinusoid (for internal state)


I don't see how this is evaluating the function twice over the sinusoid (assume that means two periods of sinusoid).

Left over from an earlier iteration (before I implemented _isstatic). Fixed.

roryyorke · 2021-01-30T05:25:34Z

control/descfcn.py

+    #
+    #    N(A) = M_1(A) e^{j \phi_1(A)} / A
+    #
+    # To compute this, we compute F(\theta) for \theta between 0 and 2 \pi,


Suggested change

# To compute this, we compute F(\theta) for \theta between 0 and 2 \pi,

# To compute this, we compute F(A\sin\theta) for \theta between 0 and 2 \pi,

roryyorke · 2021-01-30T05:26:06Z

control/descfcn.py

+    #
+    # and then integrate the product against \sin\theta and \cos\theta to obtain
+    #
+    #   \int_0^{2\pi} F(a\sin\theta) \sin\theta d\theta = M_1 \pi \cos\phi


Suggested change

# \int_0^{2\pi} F(a\sin\theta) \sin\theta d\theta = M_1 \pi \cos\phi

# \int_0^{2\pi} F(A\sin\theta) \sin\theta d\theta = M_1 \pi \cos\phi

roryyorke · 2021-01-30T05:26:21Z

control/descfcn.py

+    # and then integrate the product against \sin\theta and \cos\theta to obtain
+    #
+    #   \int_0^{2\pi} F(a\sin\theta) \sin\theta d\theta = M_1 \pi \cos\phi
+    #   \int_0^{2\pi} F(a\sin\theta) \cos\theta d\theta = M_1 \pi \sin\phi


Suggested change

# \int_0^{2\pi} F(a\sin\theta) \cos\theta d\theta = M_1 \pi \sin\phi

# \int_0^{2\pi} F(A\sin\theta) \cos\theta d\theta = M_1 \pi \sin\phi

roryyorke

Searching for the intersections will be O(m * n).

We could look libraries we can pull in to find the intersection of H(jw) and N(A). I googled "numpy polyline intersection" and found pointers to https://pypi.org/project/Shapely/, but that may be over the top.

We could also look at using scipy.optimize.root, though I don't see an obvious way to bound the domain of search. And we'd have to handle the multiple roots ourselves.

When N(A) is strictly real (saturation, deadzone, relay), we could probably specialize the solver.

roryyorke · 2021-01-30T08:01:17Z

control/descfcn.py

+        self.ub = ub
+
+    def __call__(self, x):
+        return np.maximum(self.lb, np.minimum(x, self.ub))


Suggested change

return np.maximum(self.lb, np.minimum(x, self.ub))

return np.clip(x, self.lb, self.ub)

roryyorke · 2021-01-30T08:14:13Z

control/freqplot.py

@@ -508,7 +508,7 @@ def gen_zero_centered_series(val_min, val_max, period):

 def nyquist_plot(syslist, omega=None, plot=True, label_freq=0,
                 arrowhead_length=0.1, arrowhead_width=0.1,
-                 color=None, *args, **kwargs):
+                 mirror='--', color=None, *args, **kwargs):


This new argument should be documented. Are we going to allow None for no mirror? Below you checked for False.

Updated documentation. Should either be a linestyle or False if you don't want to plot the mirror image at all. Setting it to (explicitly) to None will generate an error (from matplotlib).

roryyorke · 2021-01-30T08:23:23Z

control/descfcn.py

+        if self.lb <= A and A <= self.ub:
+            return 1.
+        else:
+            alpha, beta = math.asin(self.ub/A), math.asin(-self.lb/A)


FWIW, the expression in Slotine & Li, and also https://ocw.mit.edu/courses/aeronautics-and-astronautics/16-30-estimation-and-control-of-aerospace-systems-spring-2004/readings/gelb_ch2_ocr.pdf and https://ocw.mit.edu/resources/res-6-010-electronic-feedback-systems-spring-2013/textbook/MITRES_6-010S13_chap06.pdf are different to what you have here, but they may be equivalent. All these references restrict themselves to saturation with odd symmetry, though.

It's equivalent (I checked, but also the unit tests check this).

roryyorke · 2021-01-30T08:30:48Z

control/descfcn.py

+
+# Backlash nonlinearity (#48 in Gelb and Vander Velde, 1968)
+class backlash_nonlinearity(DescribingFunctionNonlinearity):
+    """Backlash nonlinearity for use in describing function analysis


Friction-controlled backlash? (https://ocw.mit.edu/courses/aeronautics-and-astronautics/16-30-estimation-and-control-of-aerospace-systems-spring-2004/readings/gelb_ch2_ocr.pdf p. 68?)

I saw that, but couldn't figure out what the "friction-controlled" had to do with it exactly. It seemed to me that you just had backlash as in a gear and so the output just followed the input (modulo the backlash action).

Maybe it's not important, but I remember a colleague always talking about the two backlash types (friction-controlled and inertia-controlled), and I thought it was standard nomenclature. The inertia-controlled case is discussed in the reference above immediately after friction-controlled.

Sorry, I was looking at something different (probably just the table in Appendix B, which didn't have the explanation). I'll update.

roryyorke · 2021-01-30T08:33:57Z

control/descfcn.py

+        # Convert input to an array
+        x_array = np.array(x)
+
+        y = []


there's an inconsistency here: you've allowed a vector input for this __call__, but not the others.

Good catch.

roryyorke · 2021-01-30T08:41:37Z

control/tests/descfcn_test.py

+
+
+# Static function via a class
+class saturation_class():


Suggested change

class saturation_class():

class saturation_class:

roryyorke · 2021-01-30T08:44:37Z

control/tests/descfcn_test.py

+    np.testing.assert_array_equal(satsys([0.]), [0.])
+
+    # Test SIMO nonlinearity
+    def _simofcn(t, x, u, params={}):


Suggested change

def _simofcn(t, x, u, params={}):

def _simofcn(t, x, u, params):

roryyorke · 2021-01-30T08:45:06Z

control/tests/descfcn_test.py

+        return np.array([np.sin(u[0]) * np.cos(u[1])])
+    miso_sys = ct.NonlinearIOSystem(None, outfcn=_misofcn, input=2, output=1)
+    np.testing.assert_array_equal(miso_sys([0, 0]), [0])
+    np.testing.assert_array_equal(miso_sys([0, 0]), [0])


Suggested change

np.testing.assert_array_equal(miso_sys([0, 0]), [0])

roryyorke · 2021-01-30T08:50:07Z

control/tests/descfcn_test.py

+
+
+# Test saturation describing function in multiple ways
+def test_saturation_describing_function(satsys):


I think you mean to test with the function saturation here, but this tests saturation_class twice.

roryyorke · 2021-01-30T08:52:50Z

doc/descfcn.rst

+linear system and a static nonlinearity, it is possible to obtain a
+generalization of Nyquist's stability criterion based on the idea of
+describing functions.  The basic concept involves approximating the
+response of a static nonlinearity to an input :math:`u = A e^{\omega


e^{j\omega t} ?

murrayrm force-pushed the describing_functions branch from e180716 to e82e2af Compare January 28, 2021 20:42

sawyerbfuller added a commit to sawyerbfuller/python-control that referenced this pull request Jan 28, 2021

revert some freqplot.nyquist_plot changes because they turned out to …

d195e06

…be unneccesary and to avoid a merge conbflict with python-control#521

sawyerbfuller mentioned this pull request Jan 28, 2021

Freqplot improvements #522

Merged

sawyerbfuller reviewed Jan 28, 2021

View reviewed changes

roryyorke reviewed Jan 30, 2021

View reviewed changes

murrayrm force-pushed the describing_functions branch 3 times, most recently from aba5b43 to f5bca43 Compare January 31, 2021 18:50

murrayrm mentioned this pull request Feb 6, 2021

Nyquist plot improvements #534

Merged

murrayrm added 8 commits February 6, 2021 13:09

initial implementation of describing function computation

cda5640

add intersection computation + Jupyter notebook, documentation

c85491e

change nltools to descfcn, update docs/ and docstrings

df43e69

add unit tests for exceptions/warnings + cleanup

776a814

make _isstatic internal; updated doc/ for descfcn

2ac5d9d

address @roryyorke review comments

17ed781

additional updates based on @roryyorke review feedback

b29cedf

small fixes rebasing against master

93d97cc

murrayrm force-pushed the describing_functions branch from f5bca43 to 93d97cc Compare February 6, 2021 21:23

murrayrm merged commit c04431f into python-control:master Feb 7, 2021

murrayrm deleted the describing_functions branch February 8, 2021 00:25

murrayrm mentioned this pull request Feb 11, 2021

fix isssues arising in merge of descfcn/iosys changes #540

Merged

murrayrm added this to the 0.9.0 milestone Mar 20, 2021

	"""Base class for nonlinear functions with a describing function
	"""Base class for nonlinear systems with a describing function

	that have a analytically defined describing function (accessed via the
	that have an analytically defined describing function (accessed via the

	# To compute this, we compute F(\theta) for \theta between 0 and 2 \pi,
	# To compute this, we compute F(A\sin\theta) for \theta between 0 and 2 \pi,

	# \int_0^{2\pi} F(a\sin\theta) \sin\theta d\theta = M_1 \pi \cos\phi
	# \int_0^{2\pi} F(A\sin\theta) \sin\theta d\theta = M_1 \pi \cos\phi

	# \int_0^{2\pi} F(a\sin\theta) \cos\theta d\theta = M_1 \pi \sin\phi
	# \int_0^{2\pi} F(A\sin\theta) \cos\theta d\theta = M_1 \pi \sin\phi

	return np.maximum(self.lb, np.minimum(x, self.ub))
	return np.clip(x, self.lb, self.ub)

	def _simofcn(t, x, u, params={}):
	def _simofcn(t, x, u, params):



		# Test saturation describing function in multiple ways
		def test_saturation_describing_function(satsys):

Describing functions #521

Describing functions #521

Conversation

murrayrm commented Jan 28, 2021 • edited Loading

coveralls commented Jan 28, 2021 • edited Loading

sawyerbfuller Jan 28, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

roryyorke left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

roryyorke left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

murrayrm commented Jan 28, 2021 •

edited

Loading

coveralls commented Jan 28, 2021 •

edited

Loading

sawyerbfuller Jan 28, 2021 •

edited

Loading