diff --git a/.github/conda-env/doctest-env.yml b/.github/conda-env/doctest-env.yml new file mode 100644 index 000000000..f46b239cd --- /dev/null +++ b/.github/conda-env/doctest-env.yml @@ -0,0 +1,19 @@ +name: test-env +dependencies: + - conda-build # for conda index + - pip + - coverage + - coveralls + - pytest + - pytest-cov + - pytest-timeout + - pytest-xvfb + - numpy + - matplotlib + - scipy + - sphinx + - sphinx_rtd_theme + - ipykernel + - nbsphinx + - docutils + - numpydoc diff --git a/.github/workflows/doctest.yml b/.github/workflows/doctest.yml new file mode 100644 index 000000000..62638d104 --- /dev/null +++ b/.github/workflows/doctest.yml @@ -0,0 +1,47 @@ +name: Doctest + +on: [push, pull_request] + +jobs: + doctest-linux: + # doctest needs to run only on + # latest-greatest platform with full options + runs-on: ubuntu-latest + + steps: + - name: Checkout python-control + uses: actions/checkout@v3 + + - name: Setup Conda + uses: conda-incubator/setup-miniconda@v2 + with: + python-version: 3.11 + activate-environment: test-env + environment-file: .github/conda-env/doctest-env.yml + miniforge-version: latest + miniforge-variant: Mambaforge + channels: conda-forge + channel-priority: strict + auto-update-conda: false + auto-activate-base: false + + - name: Install full dependencies + shell: bash -l {0} + run: | + mamba install cvxopt pandas slycot + + - name: Run doctest + shell: bash -l {0} + env: + PYTHON_CONTROL_ARRAY_AND_MATRIX: ${{ matrix.array-and-matrix }} + MPLBACKEND: ${{ matrix.mplbackend }} + working-directory: doc + run: | + make html + make doctest + + - name: Archive results + uses: actions/upload-artifact@v3 + with: + name: doctest-output + path: doc/_build/doctest/output.txt diff --git a/control/__init__.py b/control/__init__.py index 649a90861..340370c8b 100644 --- a/control/__init__.py +++ b/control/__init__.py @@ -42,6 +42,23 @@ """ The Python Control Systems Library :mod:`control` provides common functions for analyzing and designing feedback control systems. + +Documentation is available in two forms: docstrings provided with the code, +and the python-control users guide, available from `the python-control +homepage `_. + +The docstring examples assume that the following import commands:: + + >>> import numpy as np + >>> import control as ct + +Available subpackages +--------------------- +flatsys + Differentially flat systems +optimal + Optimization-based control + """ # Import functions from within the control system library diff --git a/control/bdalg.py b/control/bdalg.py index d1baaa410..0b1d481c8 100644 --- a/control/bdalg.py +++ b/control/bdalg.py @@ -99,9 +99,17 @@ def series(sys1, *sysn): Examples -------- - >>> sys3 = series(sys1, sys2) # Same as sys3 = sys2 * sys1 - - >>> sys5 = series(sys1, sys2, sys3, sys4) # More systems + >>> G1 = ct.rss(3) + >>> G2 = ct.rss(4) + >>> G = ct.series(G1, G2) # Same as sys3 = sys2 * sys1 + >>> G.ninputs, G.noutputs, G.nstates + (1, 1, 7) + + >>> G1 = ct.rss(2, inputs=2, outputs=3) + >>> G2 = ct.rss(3, inputs=3, outputs=1) + >>> G = ct.series(G1, G2) # Same as sys3 = sys2 * sys1 + >>> G.ninputs, G.noutputs, G.nstates + (2, 1, 5) """ from functools import reduce @@ -146,9 +154,17 @@ def parallel(sys1, *sysn): Examples -------- - >>> sys3 = parallel(sys1, sys2) # Same as sys3 = sys1 + sys2 - - >>> sys5 = parallel(sys1, sys2, sys3, sys4) # More systems + >>> G1 = ct.rss(3) + >>> G2 = ct.rss(4) + >>> G = ct.parallel(G1, G2) # Same as sys3 = sys1 + sys2 + >>> G.ninputs, G.noutputs, G.nstates + (1, 1, 7) + + >>> G1 = ct.rss(3, inputs=3, outputs=4) + >>> G2 = ct.rss(4, inputs=3, outputs=4) + >>> G = ct.parallel(G1, G2) # Add another system + >>> G.ninputs, G.noutputs, G.nstates + (3, 4, 7) """ from functools import reduce @@ -174,7 +190,13 @@ def negate(sys): Examples -------- - >>> sys2 = negate(sys1) # Same as sys2 = -sys1. + >>> G = ct.tf([2], [1, 1]) + >>> G.dcgain() + 2.0 + + >>> Gn = ct.negate(G) # Same as sys2 = -sys1. + >>> Gn.dcgain() + -2.0 """ return -sys @@ -222,6 +244,14 @@ def feedback(sys1, sys2=1, sign=-1): the corresponding feedback function is used. If `sys1` and `sys2` are both scalars, then TransferFunction.feedback is used. + Examples + -------- + >>> G = ct.rss(3, inputs=2, outputs=5) + >>> C = ct.rss(4, inputs=5, outputs=2) + >>> T = ct.feedback(G, C, sign=1) + >>> T.ninputs, T.noutputs, T.nstates + (2, 5, 7) + """ # Allow anything with a feedback function to call that function try: @@ -278,9 +308,17 @@ def append(*sys): Examples -------- - >>> sys1 = ss([[1., -2], [3., -4]], [[5.], [7]], [[6., 8]], [[9.]]) - >>> sys2 = ss([[-1.]], [[1.]], [[1.]], [[0.]]) - >>> sys = append(sys1, sys2) + >>> G1 = ct.rss(3) + >>> G2 = ct.rss(4) + >>> G = ct.append(G1, G2) + >>> G.ninputs, G.noutputs, G.nstates + (2, 2, 7) + + >>> G1 = ct.rss(3, inputs=2, outputs=4) + >>> G2 = ct.rss(4, inputs=1, outputs=4) + >>> G = ct.append(G1, G2) + >>> G.ninputs, G.noutputs, G.nstates + (3, 8, 7) """ s1 = ss._convert_to_statespace(sys[0]) @@ -323,11 +361,11 @@ def connect(sys, Q, inputv, outputv): Examples -------- - >>> sys1 = ss([[1., -2], [3., -4]], [[5.], [7]], [[6, 8]], [[9.]]) - >>> sys2 = ss([[-1.]], [[1.]], [[1.]], [[0.]]) - >>> sys = append(sys1, sys2) - >>> Q = [[1, 2], [2, -1]] # negative feedback interconnection - >>> sysc = connect(sys, Q, [2], [1, 2]) + >>> G = ct.rss(7, inputs=2, outputs=2) + >>> K = [[1, 2], [2, -1]] # negative feedback interconnection + >>> T = ct.connect(G, K, [2], [1, 2]) + >>> T.ninputs, T.noutputs, T.nstates + (1, 2, 7) Notes ----- diff --git a/control/bench/time_freqresp.py b/control/bench/time_freqresp.py index 3ae837082..4da2bcdc4 100644 --- a/control/bench/time_freqresp.py +++ b/control/bench/time_freqresp.py @@ -3,12 +3,14 @@ from numpy import logspace from timeit import timeit -nstates = 10 -sys = rss(nstates) -sys_tf = tf(sys) -w = logspace(-1,1,50) -ntimes = 1000 -time_ss = timeit("sys.freqquency_response(w)", setup="from __main__ import sys, w", number=ntimes) -time_tf = timeit("sys_tf.frequency_response(w)", setup="from __main__ import sys_tf, w", number=ntimes) -print("State-space model on %d states: %f" % (nstates, time_ss)) -print("Transfer-function model on %d states: %f" % (nstates, time_tf)) + +if __name__ == '__main__': + nstates = 10 + sys = rss(nstates) + sys_tf = tf(sys) + w = logspace(-1,1,50) + ntimes = 1000 + time_ss = timeit("sys.freqquency_response(w)", setup="from __main__ import sys, w", number=ntimes) + time_tf = timeit("sys_tf.frequency_response(w)", setup="from __main__ import sys_tf, w", number=ntimes) + print("State-space model on %d states: %f" % (nstates, time_ss)) + print("Transfer-function model on %d states: %f" % (nstates, time_tf)) diff --git a/control/canonical.py b/control/canonical.py index e714e5b8d..c84ac1f8f 100644 --- a/control/canonical.py +++ b/control/canonical.py @@ -17,6 +17,7 @@ __all__ = ['canonical_form', 'reachable_form', 'observable_form', 'modal_form', 'similarity_transform', 'bdschur'] + def canonical_form(xsys, form='reachable'): """Convert a system into canonical form @@ -36,6 +37,24 @@ def canonical_form(xsys, form='reachable'): System in desired canonical form, with state 'z' T : (M, M) real ndarray Coordinate transformation matrix, z = T * x + + Examples + -------- + >>> Gs = ct.tf2ss([1], [1, 3, 2]) + >>> Gc, T = ct.canonical_form(Gs) # default reachable + >>> Gc.B + array([[1.], + [0.]]) + + >>> Gc, T = ct.canonical_form(Gs, 'observable') + >>> Gc.C + array([[1., 0.]]) + + >>> Gc, T = ct.canonical_form(Gs, 'modal') + >>> Gc.A # doctest: +SKIP + array([[-2., 0.], + [ 0., -1.]]) + """ # Call the appropriate tranformation function @@ -65,6 +84,15 @@ def reachable_form(xsys): System in reachable canonical form, with state `z` T : (M, M) real ndarray Coordinate transformation: z = T * x + + Examples + -------- + >>> Gs = ct.tf2ss([1], [1, 3, 2]) + >>> Gc, T = ct.reachable_form(Gs) # default reachable + >>> Gc.B + array([[1.], + [0.]]) + """ # Check to make sure we have a SISO system if not issiso(xsys): @@ -119,6 +147,14 @@ def observable_form(xsys): System in observable canonical form, with state `z` T : (M, M) real ndarray Coordinate transformation: z = T * x + + Examples + -------- + >>> Gs = ct.tf2ss([1], [1, 3, 2]) + >>> Gc, T = ct.observable_form(Gs) + >>> Gc.C + array([[1., 0.]]) + """ # Check to make sure we have a SISO system if not issiso(xsys): @@ -177,6 +213,20 @@ def similarity_transform(xsys, T, timescale=1, inverse=False): zsys : StateSpace object System in transformed coordinates, with state 'z' + + Examples + -------- + >>> Gs = ct.tf2ss([1], [1, 3, 2]) + >>> Gs.A + array([[-3., -2.], + [ 1., 0.]]) + + >>> T = np.array([[0, 1], [1, 0]]) + >>> Gt = ct.similarity_transform(Gs, T) + >>> Gt.A + array([[ 0., 1.], + [-2., -3.]]) + """ # Create a new system, starting with a copy of the old one zsys = StateSpace(xsys) @@ -244,7 +294,8 @@ def _bdschur_condmax_search(aschur, tschur, condmax): Iterates mb03rd with different pmax values until: - result is non-defective; - - or condition number of similarity transform is unchanging despite large pmax; + - or condition number of similarity transform is unchanging + despite large pmax; - or condition number of similarity transform is close to condmax. Parameters @@ -283,22 +334,25 @@ def _bdschur_condmax_search(aschur, tschur, condmax): # get lower bound; try condmax ** 0.5 first pmaxlower = condmax ** 0.5 - amodal, tmodal, blksizes, eigvals = mb03rd(aschur.shape[0], aschur, tschur, pmax=pmaxlower) + amodal, tmodal, blksizes, eigvals = mb03rd( + aschur.shape[0], aschur, tschur, pmax=pmaxlower) if np.linalg.cond(tmodal) <= condmax: reslower = amodal, tmodal, blksizes, eigvals else: pmaxlower = 1.0 - amodal, tmodal, blksizes, eigvals = mb03rd(aschur.shape[0], aschur, tschur, pmax=pmaxlower) + amodal, tmodal, blksizes, eigvals = mb03rd( + aschur.shape[0], aschur, tschur, pmax=pmaxlower) cond = np.linalg.cond(tmodal) if cond > condmax: - msg = 'minimum cond={} > condmax={}; try increasing condmax'.format(cond, condmax) + msg = f"minimum {cond=} > {condmax=}; try increasing condmax" raise RuntimeError(msg) pmax = pmaxlower # phase 1: search for upper bound on pmax for i in range(50): - amodal, tmodal, blksizes, eigvals = mb03rd(aschur.shape[0], aschur, tschur, pmax=pmax) + amodal, tmodal, blksizes, eigvals = mb03rd( + aschur.shape[0], aschur, tschur, pmax=pmax) cond = np.linalg.cond(tmodal) if cond < condmax: pmaxlower = pmax @@ -319,7 +373,8 @@ def _bdschur_condmax_search(aschur, tschur, condmax): # phase 2: bisection search for i in range(50): pmax = (pmaxlower * pmaxupper) ** 0.5 - amodal, tmodal, blksizes, eigvals = mb03rd(aschur.shape[0], aschur, tschur, pmax=pmax) + amodal, tmodal, blksizes, eigvals = mb03rd( + aschur.shape[0], aschur, tschur, pmax=pmax) cond = np.linalg.cond(tmodal) if cond < condmax: @@ -370,6 +425,15 @@ def bdschur(a, condmax=None, sort=None): If `sort` is 'discrete', the blocks are sorted as for 'continuous', but applied to log of eigenvalues (i.e., continuous-equivalent eigenvalues). + + Examples + -------- + >>> Gs = ct.tf2ss([1], [1, 3, 2]) + >>> amodal, tmodal, blksizes = ct.bdschur(Gs.A) + >>> amodal #doctest: +SKIP + array([[-2., 0.], + [ 0., -1.]]) + """ if condmax is None: condmax = np.finfo(np.float64).eps ** -0.5 @@ -382,7 +446,8 @@ def bdschur(a, condmax=None, sort=None): return a.copy(), np.eye(a.shape[1], a.shape[0]), np.array([]) aschur, tschur = schur(a) - amodal, tmodal, blksizes, eigvals = _bdschur_condmax_search(aschur, tschur, condmax) + amodal, tmodal, blksizes, eigvals = _bdschur_condmax_search( + aschur, tschur, condmax) if sort in ('continuous', 'discrete'): @@ -426,9 +491,11 @@ def modal_form(xsys, condmax=None, sort=False): xsys : StateSpace object System to be transformed, with state `x` condmax : None or float, optional - An upper bound on individual transformations. If None, use `bdschur` default. + An upper bound on individual transformations. If None, use + `bdschur` default. sort : bool, optional - If False (default), Schur blocks will not be sorted. See `bdschur` for sort order. + If False (default), Schur blocks will not be sorted. See `bdschur` + for sort order. Returns ------- @@ -436,6 +503,15 @@ def modal_form(xsys, condmax=None, sort=False): System in modal canonical form, with state `z` T : (M, M) ndarray Coordinate transformation: z = T * x + + Examples + -------- + >>> Gs = ct.tf2ss([1], [1, 3, 2]) + >>> Gc, T = ct.modal_form(Gs) # default reachable + >>> Gc.A # doctest: +SKIP + array([[-2., 0.], + [ 0., -1.]]) + """ if sort: diff --git a/control/config.py b/control/config.py index 37763a6b8..abf2077af 100644 --- a/control/config.py +++ b/control/config.py @@ -65,9 +65,17 @@ def set_defaults(module, **keywords): """Set default values of parameters for a module. The set_defaults() function can be used to modify multiple parameter - values for a module at the same time, using keyword arguments: - - control.set_defaults('module', param1=val, param2=val) + values for a module at the same time, using keyword arguments. + + Examples + -------- + >>> ct.defaults['freqplot.number_of_samples'] + 1000 + >>> ct.set_defaults('freqplot', number_of_samples=100) + >>> ct.defaults['freqplot.number_of_samples'] + 100 + >>> # do some customized freqplotting + >>> ct.reset_defaults() """ if not isinstance(module, str): @@ -80,7 +88,22 @@ def set_defaults(module, **keywords): def reset_defaults(): - """Reset configuration values to their default (initial) values.""" + """Reset configuration values to their default (initial) values. + + Examples + -------- + >>> ct.defaults['freqplot.number_of_samples'] + 1000 + >>> ct.set_defaults('freqplot', number_of_samples=100) + >>> ct.defaults['freqplot.number_of_samples'] + 100 + + >>> # do some customized freqplotting + >>> ct.reset_defaults() + >>> ct.defaults['freqplot.number_of_samples'] + 1000 + + """ # System level defaults defaults.update(_control_defaults) @@ -181,6 +204,12 @@ def use_matlab_defaults(): rad/sec, with grids * State space class and functions use Numpy matrix objects + Examples + -------- + >>> ct.use_matlab_defaults() + >>> # do some matlab style plotting + >>> ct.reset_defaults() + """ set_defaults('freqplot', dB=True, deg=True, Hz=False, grid=True) set_defaults('statesp', use_numpy_matrix=True) @@ -195,6 +224,12 @@ def use_fbs_defaults(): frequency in rad/sec, no grid * Nyquist plots use dashed lines for mirror image of Nyquist curve + Examples + -------- + >>> ct.use_fbs_defaults() + >>> # do some FBS style plotting + >>> ct.reset_defaults() + """ set_defaults('freqplot', dB=False, deg=True, Hz=False, grid=False) set_defaults('nyquist', mirror_style='--') @@ -222,6 +257,13 @@ class and functions. If flat is `False`, then matrices are Prior to release 0.9.x, the default type for 2D arrays is the Numpy `matrix` class. Starting in release 0.9.0, the default type for state space operations is a 2D array. + + Examples + -------- + >>> ct.use_numpy_matrix(True, False) + >>> # do some legacy calculations using np.matrix + >>> ct.reset_defaults() + """ if flag and warn: warnings.warn("Return type numpy.matrix is deprecated.", @@ -236,6 +278,14 @@ def use_legacy_defaults(version): ---------- version : string Version number of the defaults desired. Ranges from '0.1' to '0.8.4'. + + Examples + -------- + >>> ct.use_legacy_defaults("0.9.0") + (0, 9, 0) + >>> # do some legacy style plotting + >>> ct.reset_defaults() + """ import re (major, minor, patch) = (None, None, None) # default values diff --git a/control/ctrlutil.py b/control/ctrlutil.py index 35035c7e9..fa7b91ee5 100644 --- a/control/ctrlutil.py +++ b/control/ctrlutil.py @@ -65,10 +65,17 @@ def unwrap(angle, period=2*math.pi): Examples -------- - >>> import numpy as np - >>> theta = [5.74, 5.97, 6.19, 0.13, 0.35, 0.57] - >>> unwrap(theta, period=2 * np.pi) - [5.74, 5.97, 6.19, 6.413185307179586, 6.633185307179586, 6.8531853071795865] + >>> # Already continuous + >>> theta1 = np.array([1.0, 1.5, 2.0, 2.5, 3.0]) * np.pi + >>> theta2 = ct.unwrap(theta1) + >>> theta2/np.pi # doctest: +SKIP + array([1. , 1.5, 2. , 2.5, 3. ]) + + >>> # Wrapped, discontinuous + >>> theta1 = np.array([1.0, 1.5, 0.0, 0.5, 1.0]) * np.pi + >>> theta2 = ct.unwrap(theta1) + >>> theta2/np.pi # doctest: +SKIP + array([1. , 1.5, 2. , 2.5, 3. ]) """ dangle = np.diff(angle) @@ -78,7 +85,20 @@ def unwrap(angle, period=2*math.pi): return angle def issys(obj): - """Return True if an object is a system, otherwise False""" + """Return True if an object is a Linear Time Invariant (LTI) system, + otherwise False + + Examples + -------- + >>> G = ct.tf([1], [1, 1]) + >>> ct.issys(G) + True + + >>> K = np.array([[1, 1]]) + >>> ct.issys(K) + False + + """ return isinstance(obj, lti.LTI) def db2mag(db): @@ -98,6 +118,14 @@ def db2mag(db): mag : float or ndarray corresponding magnitudes + Examples + -------- + >>> ct.db2mag(-40.0) # doctest: +SKIP + 0.01 + + >>> ct.db2mag(np.array([0, -20])) # doctest: +SKIP + array([1. , 0.1]) + """ return 10. ** (db / 20.) @@ -118,5 +146,13 @@ def mag2db(mag): db : float or ndarray corresponding values in decibels + Examples + -------- + >>> ct.mag2db(10.0) # doctest: +SKIP + 20.0 + + >>> ct.mag2db(np.array([1, 0.01])) # doctest: +SKIP + array([ 0., -40.]) + """ return 20. * np.log10(mag) diff --git a/control/delay.py b/control/delay.py index b5867ada8..d22e44107 100644 --- a/control/delay.py +++ b/control/delay.py @@ -74,6 +74,18 @@ def pade(T, n=1, numdeg=None): Ed. pp. 572-574 2. M. Vajta, "Some remarks on Padé-approximations", 3rd TEMPUS-INTCOM Symposium + + Examples + -------- + >>> delay = 1 + >>> num, den = ct.pade(delay, 3) + >>> num, den + ([-1.0, 12.0, -60.0, 120.0], [1.0, 12.0, 60.0, 120.0]) + + >>> num, den = ct.pade(delay, 3, -2) + >>> num, den + ([-6.0, 24.0], [1.0, 6.0, 18.0, 24.0]) + """ if numdeg is None: numdeg = n diff --git a/control/descfcn.py b/control/descfcn.py index 41d273f38..d0f48618c 100644 --- a/control/descfcn.py +++ b/control/descfcn.py @@ -120,6 +120,14 @@ def describing_function( TypeError If A[i] < 0 or if A[i] = 0 and the function F(0) is non-zero. + Examples + -------- + >>> F = lambda x: np.exp(-x) # Basic diode description + >>> A = np.logspace(-1, 1, 20) # Amplitudes from 0.1 to 10.0 + >>> df_values = ct.describing_function(F, A) + >>> len(df_values) + 20 + """ # If there is an analytical solution, trying using that first if try_method and hasattr(F, 'describing_function'): @@ -239,10 +247,10 @@ def describing_function_plot( Examples -------- >>> H_simple = ct.tf([8], [1, 2, 2, 1]) - >>> F_saturation = ct.descfcn.saturation_nonlinearity(1) + >>> F_saturation = ct.saturation_nonlinearity(1) >>> amp = np.linspace(1, 4, 10) - >>> ct.describing_function_plot(H_simple, F_saturation, amp) - [(3.344008947853124, 1.414213099755523)] + >>> ct.describing_function_plot(H_simple, F_saturation, amp) # doctest: +SKIP + [(3.343844998258643, 1.4142293090899216)] """ # Decide whether to turn on warnings or not @@ -354,6 +362,16 @@ class saturation_nonlinearity(DescribingFunctionNonlinearity): functions will not have zero bias and hence care must be taken in using the nonlinearity for analysis. + Examples + -------- + >>> nl = ct.saturation_nonlinearity(5) + >>> nl(1) + 1 + >>> nl(10) + 5 + >>> nl(-10) + -5 + """ def __init__(self, ub=1, lb=None): # Create the describing function nonlinearity object @@ -407,6 +425,20 @@ class relay_hysteresis_nonlinearity(DescribingFunctionNonlinearity): `-c <= x <= c`, the value depends on the branch of the hysteresis loop (as illustrated in Figure 10.20 of FBS2e). + Examples + -------- + >>> nl = ct.relay_hysteresis_nonlinearity(1, 2) + >>> nl(0) + -1 + >>> nl(1) # not enough for switching on + -1 + >>> nl(5) + 1 + >>> nl(-1) # not enough for switching off + 1 + >>> nl(-5) + -1 + """ def __init__(self, b, c): # Create the describing function nonlinearity object @@ -462,6 +494,22 @@ class friction_backlash_nonlinearity(DescribingFunctionNonlinearity): center, the output is unchanged. Otherwise, the output is given by the input shifted by `b/2`. + Examples + -------- + >>> nl = ct.friction_backlash_nonlinearity(2) # backlash of +/- 1 + >>> nl(0) + 0 + >>> nl(1) # not enough to overcome backlash + 0 + >>> nl(2) + 1.0 + >>> nl(1) + 1.0 + >>> nl(0) # not enough to overcome backlash + 1.0 + >>> nl(-1) + 0.0 + """ def __init__(self, b): diff --git a/control/dtime.py b/control/dtime.py index 724eafb76..6197ae8af 100644 --- a/control/dtime.py +++ b/control/dtime.py @@ -109,8 +109,13 @@ def sample_system(sysc, Ts, method='zoh', alpha=None, prewarp_frequency=None, Examples -------- - >>> sysc = TransferFunction([1], [1, 2, 1]) - >>> sysd = sample_system(sysc, 1, method='bilinear') + >>> Gc = ct.tf([1], [1, 2, 1]) + >>> Gc.isdtime() + False + >>> Gd = ct.sample_system(Gc, 1, method='bilinear') + >>> Gd.isdtime() + True + """ # Make sure we have a continuous time system @@ -150,8 +155,13 @@ def c2d(sysc, Ts, method='zoh', prewarp_frequency=None): Examples -------- - >>> sysc = TransferFunction([1], [1, 2, 1]) - >>> sysd = sample_system(sysc, 1, method='bilinear') + >>> Gc = ct.tf([1], [1, 2, 1]) + >>> Gc.isdtime() + False + >>> Gd = ct.sample_system(Gc, 1, method='bilinear') + >>> Gd.isdtime() + True + """ # Call the sample_system() function to do the work diff --git a/control/exception.py b/control/exception.py index 575c78c0a..e4758cc49 100644 --- a/control/exception.py +++ b/control/exception.py @@ -63,6 +63,7 @@ class ControlNotImplemented(NotImplementedError): # Utility function to see if slycot is installed slycot_installed = None def slycot_check(): + """Return True if slycot is installed, otherwise False.""" global slycot_installed if slycot_installed is None: try: @@ -76,6 +77,7 @@ def slycot_check(): # Utility function to see if pandas is installed pandas_installed = None def pandas_check(): + """Return True if pandas is installed, otherwise False.""" global pandas_installed if pandas_installed is None: try: @@ -88,6 +90,7 @@ def pandas_check(): # Utility function to see if cvxopt is installed cvxopt_installed = None def cvxopt_check(): + """Return True if cvxopt is installed, otherwise False.""" global cvxopt_installed if cvxopt_installed is None: try: diff --git a/control/flatsys/__init__.py b/control/flatsys/__init__.py index 157800073..6345ee2b9 100644 --- a/control/flatsys/__init__.py +++ b/control/flatsys/__init__.py @@ -50,6 +50,12 @@ function can be used to solve an optimal control problem with trajectory and final costs or constraints. +The docstring examples assume that the following import commands:: + + >>> import numpy as np + >>> import control as ct + >>> import control.flatsys as fs + """ # Basis function families diff --git a/control/frdata.py b/control/frdata.py index c78607a07..83873a120 100644 --- a/control/frdata.py +++ b/control/frdata.py @@ -102,7 +102,7 @@ class FrequencyResponseData(LTI): second dimension corresponding to the input index, and the 3rd dimension corresponding to the frequency points in omega. For example, - >>> frdata[2,5,:] = numpy.array([1., 0.8-0.2j, 0.2-0.8j]) + >>> frdata[2,5,:] = numpy.array([1., 0.8-0.2j, 0.2-0.8j]) # doctest: +SKIP means that the frequency response from the 6th input to the 3rd output at the frequencies defined in omega is set to the array above, i.e. the rows @@ -673,8 +673,17 @@ def _convert_to_FRD(sys, omega, inputs=1, outputs=1): scalar, then the number of inputs and outputs can be specified manually, as in: + >>> import numpy as np + >>> from control.frdata import _convert_to_FRD + + >>> omega = np.logspace(-1, 1) >>> frd = _convert_to_FRD(3., omega) # Assumes inputs = outputs = 1 - >>> frd = _convert_to_FRD(1., omegs, inputs=3, outputs=2) + >>> frd.ninputs, frd.noutputs + (1, 1) + + >>> frd = _convert_to_FRD(1., omega, inputs=3, outputs=2) + >>> frd.ninputs, frd.noutputs + (3, 2) In the latter example, sys's matrix transfer function is [[1., 1., 1.] [1., 1., 1.]]. @@ -755,5 +764,17 @@ def frd(*args): See Also -------- FRD, ss, tf + + Examples + -------- + >>> # Create from measurements + >>> response = [1.0, 1.0, 0.5] + >>> freqs = [1, 10, 100] + >>> F = ct.frd(response, freqs) + + >>> G = ct.tf([1], [1, 1]) + >>> freqs = [1, 10, 100] + >>> F = ct.frd(G, freqs) + """ return FRD(*args) diff --git a/control/freqplot.py b/control/freqplot.py index a749c3f6d..c20fb189e 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -174,8 +174,8 @@ def bode_plot(syslist, omega=None, Examples -------- - >>> sys = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.") - >>> mag, phase, omega = bode(sys) + >>> G = ct.ss([[-1, -2], [3, -4]], [[5], [7]], [[6, 8]], [[9]]) + >>> Gmag, Gphase, Gomega = ct.bode_plot(G) """ # Make a copy of the kwargs dictionary since we will modify it @@ -692,8 +692,9 @@ def nyquist_plot( Examples -------- - >>> sys = ss([[1, -2], [3, -4]], [[5], [7]], [[6, 8]], [[9]]) - >>> count = nyquist_plot(sys) + >>> G = ct.zpk([], [-1, -2, -3], gain=100) + >>> ct.nyquist_plot(G) + 2 """ # Check to see if legacy 'Plot' keyword was used @@ -1258,6 +1259,13 @@ def gangof4_plot(P, C, omega=None, **kwargs): Returns ------- None + + Examples + -------- + >>> P = ct.tf([1], [1, 1]) + >>> C = ct.tf([2], [1]) + >>> ct.gangof4_plot(P, C) + """ if not P.issiso() or not C.issiso(): # TODO: Add MIMO go4 plots. @@ -1401,15 +1409,13 @@ def singular_values_plot(syslist, omega=None, Examples -------- - >>> import numpy as np + >>> omegas = np.logspace(-4, 1, 1000) >>> den = [75, 1] - >>> sys = TransferFunction( - [[[87.8], [-86.4]], [[108.2], [-109.6]]], [[den, den], [den, den]]) - >>> omega = np.logspace(-4, 1, 1000) - >>> sigma, omega = singular_values_plot(sys, plot=True) - >>> singular_values_plot(sys, 0.0, plot=False) - (array([[197.20868123], - [ 1.39141948]]), array([0.])) + >>> G = ct.tf([[[87.8], [-86.4]], [[108.2], [-109.6]]], + ... [[den, den], [den, den]]) + >>> sigmas, omegas = ct.singular_values_plot(G, omega=omegas, plot=False) + + >>> sigmas, omegas = ct.singular_values_plot(G, 0.0, plot=False) """ @@ -1625,9 +1631,10 @@ def _default_frequency_range(syslist, Hz=None, number_of_samples=None, Examples -------- - >>> from matlab import ss - >>> sys = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.") - >>> omega = _default_frequency_range(sys) + >>> G = ct.ss([[-1, -2], [3, -4]], [[5], [7]], [[6, 8]], [[9]]) + >>> omega = ct._default_frequency_range(G) + >>> omega.min(), omega.max() + (0.1, 100.0) """ # Set default values for options @@ -1755,11 +1762,6 @@ def gen_prefix(pow1000): 'y'][8 - pow1000] # yocto (10^-24) -def find_nearest_omega(omega_list, omega): - omega_list = np.asarray(omega_list) - return omega_list[(np.abs(omega_list - omega)).argmin()] - - # Function aliases bode = bode_plot nyquist = nyquist_plot diff --git a/control/grid.py b/control/grid.py index 07ca4a59d..785ec2743 100644 --- a/control/grid.py +++ b/control/grid.py @@ -156,7 +156,7 @@ def nogrid(): def zgrid(zetas=None, wns=None, ax=None): - '''Draws discrete damping and frequency grid''' + """Draws discrete damping and frequency grid""" fig = plt.gcf() if ax is None: diff --git a/control/iosys.py b/control/iosys.py index 6b0f6cfaa..97e1fb63d 100644 --- a/control/iosys.py +++ b/control/iosys.py @@ -2371,12 +2371,14 @@ def ss(*args, **kwargs): Examples -------- - >>> # Create a Linear I/O system object from from for matrices - >>> sys1 = ss([[1, -2], [3 -4]], [[5], [7]], [[6, 8]], [[9]]) + Create a Linear I/O system object from matrices. - >>> # Convert a TransferFunction to a StateSpace object. - >>> sys_tf = tf([2.], [1., 3]) - >>> sys2 = ss(sys_tf) + >>> G = ct.ss([[-1, -2], [3, -4]], [[5], [7]], [[6, 8]], [[9]]) + + Convert a TransferFunction to a StateSpace object. + + >>> sys_tf = ct.tf([2.], [1., 3]) + >>> sys2 = ct.ss(sys_tf) """ # See if this is a nonlinear I/O system @@ -2495,6 +2497,15 @@ def drss(*args, **kwargs): function calls :func:`rss` using either the `dt` keyword provided by the user or `dt=True` if not specified. + Examples + -------- + >>> G = ct.drss(states=4, outputs=2, inputs=1) + >>> G.ninputs, G.noutputs, G.nstates + (1, 2, 4) + >>> G.isdtime() + True + + """ # Make sure the timebase makes sense if 'dt' in kwargs: @@ -2585,10 +2596,12 @@ def tf2io(*args, **kwargs): -------- >>> num = [[[1., 2.], [3., 4.]], [[5., 6.], [7., 8.]]] >>> den = [[[9., 8., 7.], [6., 5., 4.]], [[3., 2., 1.], [-1., -2., -3.]]] - >>> sys1 = tf2ss(num, den) + >>> sys1 = ct.tf2ss(num, den) - >>> sys_tf = tf(num, den) - >>> sys2 = tf2ss(sys_tf) + >>> sys_tf = ct.tf(num, den) + >>> G = ct.tf2ss(sys_tf) + >>> G.ninputs, G.noutputs, G.nstates + (2, 2, 8) """ # Convert the system to a state space system @@ -2765,26 +2778,25 @@ def interconnect( Examples -------- - >>> P = control.LinearIOSystem( - >>> control.rss(2, 2, 2, strictly_proper=True), name='P') - >>> C = control.LinearIOSystem(control.rss(2, 2, 2), name='C') - >>> T = control.interconnect( - >>> [P, C], - >>> connections = [ - >>> ['P.u[0]', 'C.y[0]'], ['P.u[1]', 'C.y[1]'], - >>> ['C.u[0]', '-P.y[0]'], ['C.u[1]', '-P.y[1]']], - >>> inplist = ['C.u[0]', 'C.u[1]'], - >>> outlist = ['P.y[0]', 'P.y[1]'], - >>> ) + >>> P = ct.rss(2, 2, 2, strictly_proper=True, name='P') + >>> C = ct.rss(2, 2, 2, name='C') + >>> T = ct.interconnect( + ... [P, C], + ... connections = [ + ... ['P.u[0]', 'C.y[0]'], ['P.u[1]', 'C.y[1]'], + ... ['C.u[0]', '-P.y[0]'], ['C.u[1]', '-P.y[1]']], + ... inplist = ['C.u[0]', 'C.u[1]'], + ... outlist = ['P.y[0]', 'P.y[1]'], + ... ) For a SISO system, this example can be simplified by using the :func:`~control.summing_block` function and the ability to automatically interconnect signals with the same names: - >>> P = control.tf(1, [1, 0], inputs='u', outputs='y') - >>> C = control.tf(10, [1, 1], inputs='e', outputs='u') - >>> sumblk = control.summing_junction(inputs=['r', '-y'], output='e') - >>> T = control.interconnect([P, C, sumblk], inputs='r', outputs='y') + >>> P = ct.tf(1, [1, 0], inputs='u', outputs='y') + >>> C = ct.tf(10, [1, 1], inputs='e', outputs='u') + >>> sumblk = ct.summing_junction(inputs=['r', '-y'], output='e') + >>> T = ct.interconnect([P, C, sumblk], inputs='r', outputs='y') Notes ----- @@ -2986,10 +2998,12 @@ def summing_junction( Examples -------- - >>> P = control.tf2io(ct.tf(1, [1, 0]), inputs='u', outputs='y') - >>> C = control.tf2io(ct.tf(10, [1, 1]), inputs='e', outputs='u') - >>> sumblk = control.summing_junction(inputs=['r', '-y'], output='e') - >>> T = control.interconnect((P, C, sumblk), inputs='r', outputs='y') + >>> P = ct.tf2io(1, [1, 0], inputs='u', outputs='y') + >>> C = ct.tf2io(10, [1, 1], inputs='e', outputs='u') + >>> sumblk = ct.summing_junction(inputs=['r', '-y'], output='e') + >>> T = ct.interconnect([P, C, sumblk], inputs='r', outputs='y') + >>> T.ninputs, T.noutputs, T.nstates + (1, 1, 2) """ # Utility function to parse input and output signal lists diff --git a/control/lti.py b/control/lti.py index 1bc08229d..6dc3bc62c 100644 --- a/control/lti.py +++ b/control/lti.py @@ -324,6 +324,13 @@ def damp(sys, doprint=True): wn = abs(s) Z = -real(s)/wn. + Examples + -------- + >>> G = ct.tf([1], [1, 4]) + >>> wn, damping, poles = ct.damp(G) + _____Eigenvalue______ Damping___ Frequency_ + -4 1 4 + """ wn, damping, poles = sys.damp() if doprint: @@ -386,10 +393,8 @@ def evalfr(sys, x, squeeze=None): Examples -------- - >>> sys = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.") - >>> evalfr(sys, 1j) - array([[ 44.8-21.4j]]) - >>> # This is the transfer function matrix evaluated at s = i. + >>> G = ct.ss([[-1, -2], [3, -4]], [[5], [7]], [[6, 8]], [[9]]) + >>> fresp = ct.evalfr(G, 1j) # evaluate at s = 1j .. todo:: Add example with MIMO system @@ -449,12 +454,8 @@ def frequency_response(sys, omega, squeeze=None): Examples -------- - >>> sys = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.") - >>> mag, phase, omega = freqresp(sys, [0.1, 1., 10.]) - >>> mag - array([[[ 58.8576682 , 49.64876635, 13.40825927]]]) - >>> phase - array([[[-0.05408304, -0.44563154, -0.66837155]]]) + >>> G = ct.ss([[-1, -2], [3, -4]], [[5], [7]], [[6, 8]], [[9]]) + >>> mag, phase, omega = ct.freqresp(G, [0.1, 1., 10.]) .. todo:: Add example with MIMO system @@ -488,6 +489,12 @@ def dcgain(sys): the origin, (nan + nanj) if there is a pole/zero cancellation at the origin. + Examples + -------- + >>> G = ct.tf([1], [1, 2]) + >>> ct.dcgain(G) # doctest: +SKIP + 0.5 + """ return sys.dcgain() diff --git a/control/margins.py b/control/margins.py index 662634086..28daaf358 100644 --- a/control/margins.py +++ b/control/margins.py @@ -476,9 +476,9 @@ def phase_crossover_frequencies(sys): Examples -------- - >>> tf = TransferFunction([1], [1, 2, 3, 4]) - >>> phase_crossover_frequencies(tf) - (array([ 1.73205081, 0. ]), array([-0.5 , 0.25])) + >>> G = ct.tf([1], [1, 2, 3, 4]) + >>> x_omega, x_gain = ct.phase_crossover_frequencies(G) + """ # Convert to a transfer function tf = xferfcn._convert_to_transfer_function(sys) @@ -537,8 +537,8 @@ def margin(*args): Examples -------- - >>> sys = tf(1, [1, 2, 1, 0]) - >>> gm, pm, wcg, wcp = margin(sys) + >>> G = ct.tf(1, [1, 2, 1, 0]) + >>> gm, pm, wcg, wcp = ct.margin(G) """ if len(args) == 1: diff --git a/control/matlab/timeresp.py b/control/matlab/timeresp.py index 58b5e589d..5420bfdf4 100644 --- a/control/matlab/timeresp.py +++ b/control/matlab/timeresp.py @@ -46,7 +46,11 @@ def step(sys, T=None, X0=0., input=0, output=None, return_x=False): Examples -------- - >>> yout, T = step(sys, T, X0) + >>> from control.matlab import step, rss + + >>> G = rss(4) + >>> yout, T = step(G) + ''' from ..timeresp import step_response @@ -115,7 +119,11 @@ def stepinfo(sysdata, T=None, yfinal=None, SettlingTimeThreshold=0.02, Examples -------- - >>> S = stepinfo(sys, T) + >>> from control.matlab import stepinfo, rss + + >>> G = rss(4) + >>> S = stepinfo(G) + """ from ..timeresp import step_info @@ -166,7 +174,11 @@ def impulse(sys, T=None, X0=0., input=0, output=None, return_x=False): Examples -------- - >>> yout, T = impulse(sys, T) + >>> from control.matlab import rss, impulse + + >>> G = rss() + >>> yout, T = impulse(G) + ''' from ..timeresp import impulse_response @@ -214,7 +226,10 @@ def initial(sys, T=None, X0=0., input=None, output=None, return_x=False): Examples -------- - >>> yout, T = initial(sys, T, X0) + >>> from control.matlab import initial, rss + + >>> G = rss(4) + >>> yout, T = initial(G) ''' from ..timeresp import initial_response @@ -261,7 +276,12 @@ def lsim(sys, U=0., T=None, X0=0.): Examples -------- - >>> yout, T, xout = lsim(sys, U, T, X0) + >>> from control.matlab import rss, lsim + + >>> G = rss(4) + >>> T = np.linspace(0,10) + >>> yout, T, xout = lsim(G, T=T) + ''' from ..timeresp import forced_response diff --git a/control/matlab/wrappers.py b/control/matlab/wrappers.py index 22ec95f39..1c4ee053a 100644 --- a/control/matlab/wrappers.py +++ b/control/matlab/wrappers.py @@ -26,11 +26,13 @@ def bode(*args, **kwargs): a list of systems can be entered, or several systems can be specified (i.e. several parameters). The sys arguments may also be interspersed with format strings. A frequency argument (array_like) - may also be added, some examples: - * >>> bode(sys, w) # one system, freq vector - * >>> bode(sys1, sys2, ..., sysN) # several systems - * >>> bode(sys1, sys2, ..., sysN, w) - * >>> bode(sys1, 'plotstyle1', ..., sysN, 'plotstyleN') # + plot formats + may also be added, some examples:: + + >>> bode(sys, w) # one system, freq vector # doctest: +SKIP + >>> bode(sys1, sys2, ..., sysN) # several systems # doctest: +SKIP + >>> bode(sys1, sys2, ..., sysN, w) # doctest: +SKIP + >>> bode(sys1, 'plotstyle1', ..., sysN, 'plotstyleN') # + plot formats # doctest: +SKIP + omega: freq_range Range of frequencies in rad/s dB : boolean @@ -44,6 +46,8 @@ def bode(*args, **kwargs): Examples -------- + >>> from control.matlab import ss, bode + >>> sys = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.") >>> mag, phase, omega = bode(sys) @@ -51,10 +55,10 @@ def bode(*args, **kwargs): Document these use cases - * >>> bode(sys, w) - * >>> bode(sys1, sys2, ..., sysN) - * >>> bode(sys1, sys2, ..., sysN, w) - * >>> bode(sys1, 'plotstyle1', ..., sysN, 'plotstyleN') + * >>> bode(sys, w) # doctest: +SKIP + * >>> bode(sys1, sys2, ..., sysN) # doctest: +SKIP + * >>> bode(sys1, sys2, ..., sysN, w) # doctest: +SKIP + * >>> bode(sys1, 'plotstyle1', ..., sysN, 'plotstyleN') # doctest: +SKIP """ from ..freqplot import bode_plot diff --git a/control/modelsimp.py b/control/modelsimp.py index b1c1ae31c..f7b15093d 100644 --- a/control/modelsimp.py +++ b/control/modelsimp.py @@ -84,7 +84,10 @@ def hsvd(sys): Examples -------- - >>> H = hsvd(sys) + >>> G = ct.tf2ss([1], [1, 2]) + >>> H = ct.hsvd(G) + >>> H[0] + 0.25 """ # TODO: implement for discrete time systems @@ -135,7 +138,11 @@ def modred(sys, ELIM, method='matchdc'): Examples -------- - >>> rsys = modred(sys, ELIM, method='truncate') + >>> G = ct.rss(4) + >>> Gr = ct.modred(G, [0, 2], method='matchdc') + >>> Gr.nstates + 2 + """ # Check for ss system object, need a utility for this? @@ -255,7 +262,10 @@ def balred(sys, orders, method='truncate', alpha=None): Examples -------- - >>> rsys = balred(sys, orders, method='truncate') + >>> G = ct.rss(4) + >>> Gr = ct.balred(G, orders=2, method='matchdc') + >>> Gr.nstates + 2 """ if method != 'truncate' and method != 'matchdc': @@ -386,7 +396,7 @@ def era(YY, m, n, nin, nout, r): Examples -------- - >>> rsys = era(YY, m, n, nin, nout, r) + >>> rsys = era(YY, m, n, nin, nout, r) # doctest: +SKIP """ raise NotImplementedError('This function is not implemented yet.') @@ -446,10 +456,10 @@ def markov(Y, U, m=None, transpose=False): Examples -------- - >>> T = numpy.linspace(0, 10, 100) - >>> U = numpy.ones((1, 100)) - >>> T, Y, _ = forced_response(tf([1], [1, 0.5], True), T, U) - >>> H = markov(Y, U, 3, transpose=False) + >>> T = np.linspace(0, 10, 100) + >>> U = np.ones((1, 100)) + >>> T, Y = ct.forced_response(ct.tf([1], [1, 0.5], True), T, U) + >>> H = ct.markov(Y, U, 3, transpose=False) """ # Convert input parameters to 2D arrays (if they aren't already) diff --git a/control/optimal.py b/control/optimal.py index 9a8218a91..77be02d2a 100644 --- a/control/optimal.py +++ b/control/optimal.py @@ -6,6 +6,12 @@ """The :mod:`~control.optimal` module provides support for optimization-based controllers for nonlinear systems with state and input constraints. +The docstring examples assume that the following import commands:: + + >>> import numpy as np + >>> import control as ct + >>> import control.optimal as obc + """ import numpy as np diff --git a/control/phaseplot.py b/control/phaseplot.py index 6a4be5ca6..91d7b79b0 100644 --- a/control/phaseplot.py +++ b/control/phaseplot.py @@ -115,9 +115,6 @@ def phase_plot(odefun, X=None, Y=None, scale=1, X0=None, T=None, -------- box_grid : construct box-shaped grid of initial conditions - Examples - -------- - """ # diff --git a/control/robust.py b/control/robust.py index aa5c922dc..a0e53d199 100644 --- a/control/robust.py +++ b/control/robust.py @@ -70,7 +70,24 @@ def h2syn(P, nmeas, ncon): Examples -------- - >>> K = h2syn(P,nmeas,ncon) + >>> # Unstable first order SISI system + >>> G = ct.tf([1], [1, -1], inputs=['u'], outputs=['y']) + >>> max(G.poles()) < 0 # Is G stable? + False + + >>> # Create partitioned system with trivial unity systems + >>> P11 = ct.tf([0], [1], inputs=['w'], outputs=['z']) + >>> P12 = ct.tf([1], [1], inputs=['u'], outputs=['z']) + >>> P21 = ct.tf([1], [1], inputs=['w'], outputs=['y']) + >>> P22 = G + >>> P = ct.interconnect([P11, P12, P21, P22], + ... inplist=['w', 'u'], outlist=['z', 'y']) + + >>> # Synthesize H2 optimal stabilizing controller + >>> K = ct.h2syn(P, nmeas=1, ncon=1) + >>> T = ct.feedback(G, K, sign=1) + >>> max(T.poles()) < 0 # Is T stable? + True """ # Check for ss system object, need a utility for this? @@ -134,7 +151,23 @@ def hinfsyn(P, nmeas, ncon): Examples -------- - >>> K, CL, gam, rcond = hinfsyn(P,nmeas,ncon) + >>> # Unstable first order SISI system + >>> G = ct.tf([1], [1,-1], inputs=['u'], outputs=['y']) + >>> max(G.poles()) < 0 + False + + >>> # Create partitioned system with trivial unity systems + >>> P11 = ct.tf([0], [1], inputs=['w'], outputs=['z']) + >>> P12 = ct.tf([1], [1], inputs=['u'], outputs=['z']) + >>> P21 = ct.tf([1], [1], inputs=['w'], outputs=['y']) + >>> P22 = G + >>> P = ct.interconnect([P11, P12, P21, P22], inplist=['w', 'u'], outlist=['z', 'y']) + + >>> # Synthesize Hinf optimal stabilizing controller + >>> K, CL, gam, rcond = ct.hinfsyn(P, nmeas=1, ncon=1) + >>> T = ct.feedback(G, K, sign=1) + >>> max(T.poles()) < 0 + True """ @@ -221,28 +254,33 @@ def _size_as_needed(w, wname, n): def augw(g, w1=None, w2=None, w3=None): """Augment plant for mixed sensitivity problem. - Parameters - ---------- - g: LTI object, ny-by-nu - w1: weighting on S; None, scalar, or k1-by-ny LTI object - w2: weighting on KS; None, scalar, or k2-by-nu LTI object - w3: weighting on T; None, scalar, or k3-by-ny LTI object - p: augmented plant; StateSpace object - If a weighting is None, no augmentation is done for it. At least one weighting must not be None. If a weighting w is scalar, it will be replaced by I*w, where I is ny-by-ny for w1 and w3, and nu-by-nu for w2. + Parameters + ---------- + g: LTI object, ny-by-nu + Plant + w1: None, scalar, or k1-by-ny LTI object + Weighting on S + w2: None, scalar, or k2-by-nu LTI object + Weighting on KS + w3: None, scalar, or k3-by-ny LTI object + Weighting on T + Returns ------- - p: plant augmented with weightings, suitable for submission to hinfsyn or h2syn. + p: StateSpace + Plant augmented with weightings, suitable for submission to hinfsyn or + h2syn. Raises ------ ValueError - - if all weightings are None + If all weightings are None See Also -------- @@ -331,25 +369,35 @@ def mixsyn(g, w1=None, w2=None, w3=None): Parameters ---------- - g: LTI; the plant for which controller must be synthesized - w1: weighting on s = (1+g*k)**-1; None, or scalar or k1-by-ny LTI - w2: weighting on k*s; None, or scalar or k2-by-nu LTI - w3: weighting on t = g*k*(1+g*k)**-1; None, or scalar or k3-by-ny LTI - At least one of w1, w2, and w3 must not be None. + g: LTI + The plant for which controller must be synthesized + w1: None, or scalar or k1-by-ny LTI + Weighting on S = (1+G*K)**-1 + w2: None, or scalar or k2-by-nu LTI + Weighting on K*S + w3: None, or scalar or k3-by-ny LTI + Weighting on T = G*K*(1+G*K)**-1; Returns ------- - k: synthesized controller; StateSpace object - cl: closed system mapping evaluation inputs to evaluation outputs; if - p is the augmented plant, with - [z] = [p11 p12] [w], - [y] [p21 g] [u] - then cl is the system from w->z with u=-k*y. StateSpace object. - - info: tuple with entries, in order, - - gamma: scalar; H-infinity norm of cl - - rcond: array; estimates of reciprocal condition numbers - computed during synthesis. See hinfsyn for details + k: StateSpace + Synthesized controller; + cl: StateSpace + Closed system mapping evaluation inputs to evaluation outputs. + + Let p be the augmented plant, with:: + + [z] = [p11 p12] [w] + [y] [p21 g] [u] + + then cl is the system from w->z with `u = -k*y`. + + info: tuple + gamma: scalar + H-infinity norm of cl + rcond: array + Estimates of reciprocal condition numbers + computed during synthesis. See hinfsyn for details If a weighting w is scalar, it will be replaced by I*w, where I is ny-by-ny for w1 and w3, and nu-by-nu for w2. diff --git a/control/sisotool.py b/control/sisotool.py index 8a3b3d9f7..c1b260d08 100644 --- a/control/sisotool.py +++ b/control/sisotool.py @@ -81,8 +81,8 @@ def sisotool(sys, initial_gain=None, xlim_rlocus=None, ylim_rlocus=None, Examples -------- - >>> sys = tf([1000], [1,25,100,0]) - >>> sisotool(sys) + >>> G = ct.tf([1000], [1, 25, 100, 0]) + >>> ct.sisotool(G) # doctest: +SKIP """ from .rlocus import root_locus diff --git a/control/statefbk.py b/control/statefbk.py index c76a4e31a..20d99dff6 100644 --- a/control/statefbk.py +++ b/control/statefbk.py @@ -122,7 +122,7 @@ def place(A, B, p): -------- >>> A = [[-1, -1], [0, 1]] >>> B = [[0], [1]] - >>> K = place(A, B, [-2, -5]) + >>> K = ct.place(A, B, [-2, -5]) See Also -------- @@ -375,8 +375,8 @@ def lqr(*args, **kwargs): Examples -------- - >>> K, S, E = lqr(sys, Q, R, [N]) - >>> K, S, E = lqr(A, B, Q, R, [N]) + >>> K, S, E = lqr(sys, Q, R, [N]) # doctest: +SKIP + >>> K, S, E = lqr(A, B, Q, R, [N]) # doctest: +SKIP """ # @@ -520,9 +520,8 @@ def dlqr(*args, **kwargs): Examples -------- - >>> K, S, E = dlqr(dsys, Q, R, [N]) - >>> K, S, E = dlqr(A, B, Q, R, [N]) - + >>> K, S, E = dlqr(dsys, Q, R, [N]) # doctest: +SKIP + >>> K, S, E = dlqr(A, B, Q, R, [N]) # doctest: +SKIP """ # @@ -993,7 +992,10 @@ def ctrb(A, B): Examples -------- - >>> C = ctrb(A, B) + >>> G = ct.tf2ss([1], [1, 2, 3]) + >>> C = ct.ctrb(G.A, G.B) + >>> np.linalg.matrix_rank(C) + 2 """ @@ -1029,7 +1031,11 @@ def obsv(A, C): Examples -------- - >>> O = obsv(A, C) + >>> G = ct.tf2ss([1], [1, 2, 3]) + >>> C = ct.obsv(G.A, G.C) + >>> np.linalg.matrix_rank(C) + 2 + """ # Convert input parameters to matrices (if they aren't already) @@ -1078,10 +1084,11 @@ def gram(sys, type): Examples -------- - >>> Wc = gram(sys, 'c') - >>> Wo = gram(sys, 'o') - >>> Rc = gram(sys, 'cf'), where Wc = Rc' * Rc - >>> Ro = gram(sys, 'of'), where Wo = Ro' * Ro + >>> G = ct.rss(4) + >>> Wc = ct.gram(G, 'c') + >>> Wo = ct.gram(G, 'o') + >>> Rc = ct.gram(G, 'cf') # where Wc = Rc' * Rc + >>> Ro = ct.gram(G, 'of') # where Wo = Ro' * Ro """ diff --git a/control/statesp.py b/control/statesp.py index 9fff28d27..41f92ae21 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -1217,8 +1217,8 @@ def returnScipySignalLTI(self, strict=True): For instance, - >>> out = ssobject.returnScipySignalLTI() - >>> out[3][5] + >>> out = ssobject.returnScipySignalLTI() # doctest: +SKIP + >>> out[3][5] # doctest: +SKIP is a :class:`scipy.signal.lti` object corresponding to the transfer function from the 6th input to the 4th output. @@ -1362,8 +1362,8 @@ def sample(self, Ts, method='zoh', alpha=None, prewarp_frequency=None, Examples -------- - >>> sys = StateSpace(0, 1, 1, 0) - >>> sysd = sys.sample(0.5, method='bilinear') + >>> G = ct.ss(0, 1, 1, 0) + >>> sysd = G.sample(0.5, method='bilinear') """ if not self.isctime(): @@ -1526,14 +1526,7 @@ def _convert_to_statespace(sys): If sys is already a state space, then it is returned. If sys is a transfer function object, then it is converted to a state space and - returned. If sys is a scalar, then the number of inputs and outputs can - be specified manually, as in: - - >>> sys = _convert_to_statespace(3.) # Assumes inputs = outputs = 1 - >>> sys = _convert_to_statespace(1., inputs=3, outputs=2) - - In the latter example, A = B = C = 0 and D = [[1., 1., 1.] - [1., 1., 1.]]. + returned. Note: no renaming of inputs and outputs is performed; this should be done by the calling function. @@ -1898,10 +1891,10 @@ def tf2ss(*args, **kwargs): -------- >>> num = [[[1., 2.], [3., 4.]], [[5., 6.], [7., 8.]]] >>> den = [[[9., 8., 7.], [6., 5., 4.]], [[3., 2., 1.], [-1., -2., -3.]]] - >>> sys1 = tf2ss(num, den) + >>> sys1 = ct.tf2ss(num, den) - >>> sys_tf = tf(num, den) - >>> sys2 = tf2ss(sys_tf) + >>> sys_tf = ct.tf(num, den) + >>> sys2 = ct.tf2ss(sys_tf) """ diff --git a/control/stochsys.py b/control/stochsys.py index 000c1b6ab..65b4ebd95 100644 --- a/control/stochsys.py +++ b/control/stochsys.py @@ -108,8 +108,8 @@ def lqe(*args, **kwargs): Examples -------- - >>> L, P, E = lqe(A, G, C, QN, RN) - >>> L, P, E = lqe(A, G, C, Q, RN, NN) + >>> L, P, E = lqe(A, G, C, QN, RN) # doctest: +SKIP + >>> L, P, E = lqe(A, G, C, Q, RN, NN) # doctest: +SKIP See Also -------- @@ -240,8 +240,8 @@ def dlqe(*args, **kwargs): Examples -------- - >>> L, P, E = dlqe(A, G, C, QN, RN) - >>> L, P, E = dlqe(A, G, C, QN, RN, NN) + >>> L, P, E = dlqe(A, G, C, QN, RN) # doctest: +SKIP + >>> L, P, E = dlqe(A, G, C, QN, RN, NN) # doctest: +SKIP See Also -------- diff --git a/control/timeresp.py b/control/timeresp.py index 24f553a32..638a07329 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -916,7 +916,9 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False, Examples -------- - >>> T, yout, xout = forced_response(sys, T, u, X0) + >>> G = ct.rss(4) + >>> T = np.linspace(0, 10) + >>> T, yout = ct.forced_response(G, T=T) See :ref:`time-series-convention` and :ref:`package-configuration-parameters`. @@ -1328,7 +1330,8 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None, Examples -------- - >>> T, yout = step_response(sys, T, X0) + >>> G = ct.rss(4) + >>> T, yout = ct.step_response(G) """ # Create the time and input vectors @@ -1440,9 +1443,8 @@ def step_info(sysdata, T=None, T_num=None, yfinal=None, Examples -------- - >>> from control import step_info, TransferFunction - >>> sys = TransferFunction([-1, 1], [1, 1, 1]) - >>> S = step_info(sys) + >>> sys = ct.TransferFunction([-1, 1], [1, 1, 1]) + >>> S = ct.step_info(sys) >>> for k in S: ... print(f"{k}: {S[k]:3.4}") ... @@ -1460,15 +1462,14 @@ def step_info(sysdata, T=None, T_num=None, yfinal=None, characteristics for the second input and specify a 5% error until the signal is considered settled. - >>> from numpy import sqrt - >>> from control import step_info, StateSpace - >>> sys = StateSpace([[-1., -1.], + >>> from math import sqrt + >>> sys = ct.StateSpace([[-1., -1.], ... [1., 0.]], ... [[-1./sqrt(2.), 1./sqrt(2.)], ... [0, 0]], ... [[sqrt(2.), -sqrt(2.)]], ... [[0, 0]]) - >>> S = step_info(sys, T=10., SettlingTimeThreshold=0.05) + >>> S = ct.step_info(sys, T=10., SettlingTimeThreshold=0.05) >>> for k, v in S[0][1].items(): ... print(f"{k}: {float(v):3.4}") RiseTime: 1.212 @@ -1686,7 +1687,8 @@ def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None, Examples -------- - >>> T, yout = initial_response(sys, T, X0) + >>> G = ct.rss(4) + >>> T, yout = ct.initial_response(G) """ squeeze, sys = _get_ss_simo(sys, input, output, squeeze=squeeze) @@ -1801,7 +1803,8 @@ def impulse_response(sys, T=None, X0=0., input=None, output=None, T_num=None, Examples -------- - >>> T, yout = impulse_response(sys, T, X0) + >>> G = ct.rss(4) + >>> T, yout = ct.impulse_response(G) """ # Convert to state space so that we can simulate diff --git a/control/xferfcn.py b/control/xferfcn.py index 6edb26858..89e2546f8 100644 --- a/control/xferfcn.py +++ b/control/xferfcn.py @@ -121,7 +121,7 @@ class TransferFunction(LTI): The attribues 'num' and 'den' are 2-D lists of arrays containing MIMO numerator and denominator coefficients. For example, - >>> num[2][5] = numpy.array([1., 4., 8.]) + >>> num[2][5] = numpy.array([1., 4., 8.]) # doctest: +SKIP means that the numerator of the transfer function from the 6th input to the 3rd output is set to s^2 + 4s + 8. @@ -152,7 +152,7 @@ class TransferFunction(LTI): discrete time. These can be used to create variables that allow algebraic creation of transfer functions. For example, - >>> s = TransferFunction.s + >>> s = ct.TransferFunction.s >>> G = (s + 1)/(s**2 + 2*s + 1) """ @@ -475,7 +475,8 @@ def __str__(self, var=None): denstr = _tf_polynomial_to_string(self.den[no][ni], var=var) elif self.display_format == 'zpk': z, p, k = tf2zpk(self.num[no][ni], self.den[no][ni]) - numstr = _tf_factorized_polynomial_to_string(z, gain=k, var=var) + numstr = _tf_factorized_polynomial_to_string( + z, gain=k, var=var) denstr = _tf_factorized_polynomial_to_string(p, var=var) # Figure out the length of the separating line @@ -532,7 +533,8 @@ def _repr_latex_(self, var=None): denstr = _tf_polynomial_to_string(self.den[no][ni], var=var) elif self.display_format == 'zpk': z, p, k = tf2zpk(self.num[no][ni], self.den[no][ni]) - numstr = _tf_factorized_polynomial_to_string(z, gain=k, var=var) + numstr = _tf_factorized_polynomial_to_string( + z, gain=k, var=var) denstr = _tf_factorized_polynomial_to_string(p, var=var) numstr = _tf_string_to_latex(numstr, var=var) @@ -907,8 +909,8 @@ def returnScipySignalLTI(self, strict=True): For instance, - >>> out = tfobject.returnScipySignalLTI() - >>> out[3][5] + >>> out = tfobject.returnScipySignalLTI() # doctest: +SKIP + >>> out[3][5] # doctest: +SKIP is a :class:`scipy.signal.lti` object corresponding to the transfer function from the 6th input to the 4th output. @@ -996,7 +998,7 @@ def _common_den(self, imag_tol=None, allow_nonproper=False): Examples -------- - >>> num, den, denorder = sys._common_den() + >>> num, den, denorder = sys._common_den() # doctest: +SKIP """ @@ -1178,7 +1180,7 @@ def sample(self, Ts, method='zoh', alpha=None, prewarp_frequency=None, Examples -------- - >>> sys = TransferFunction(1, [1,1]) + >>> sys = ct.tf(1, [1, 1]) >>> sysd = sys.sample(0.5, method='bilinear') """ @@ -1235,6 +1237,13 @@ def dcgain(self, warn_infinite=False): For real valued systems, the empty imaginary part of the complex zero-frequency response is discarded and a real array or scalar is returned. + + Examples + -------- + >>> G = ct.tf([1], [1, 4]) + >>> G.dcgain() + 0.25 + """ return self._dcgain(warn_infinite) @@ -1263,8 +1272,8 @@ def _isstatic(self): #: #: Example #: ------- - #: >>> s = TransferFunction.s - #: >>> G = (s + 1)/(s**2 + 2*s + 1) + #: >>> s = TransferFunction.s # doctest: +SKIP + #: >>> G = (s + 1)/(s**2 + 2*s + 1) # doctest: +SKIP #: #: :meta hide-value: s = None @@ -1276,8 +1285,8 @@ def _isstatic(self): #: #: Example #: ------- - #: >>> z = TransferFunction.z - #: >>> G = 2 * z / (4 * z**3 + 3*z - 1) + #: >>> z = TransferFunction.z # doctest: +SKIP + #: >>> G = 2 * z / (4 * z**3 + 3*z - 1) # doctest: +SKIP #: #: :meta hide-value: z = None @@ -1613,15 +1622,15 @@ def tf(*args, **kwargs): >>> # (3s + 4) / (6s^2 + 5s + 4). >>> num = [[[1., 2.], [3., 4.]], [[5., 6.], [7., 8.]]] >>> den = [[[9., 8., 7.], [6., 5., 4.]], [[3., 2., 1.], [-1., -2., -3.]]] - >>> sys1 = tf(num, den) + >>> sys1 = ct.tf(num, den) >>> # Create a variable 's' to allow algebra operations for SISO systems - >>> s = tf('s') - >>> G = (s + 1) / (s**2 + 2*s + 1) + >>> s = ct.tf('s') + >>> G = (s + 1)/(s**2 + 2*s + 1) >>> # Convert a StateSpace to a TransferFunction object. - >>> sys_ss = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.") - >>> sys2 = tf(sys1) + >>> sys_ss = ct.ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.") + >>> sys2 = ct.tf(sys1) """ @@ -1700,12 +1709,13 @@ def zpk(zeros, poles, gain, *args, **kwargs): Examples -------- - >>> from control import tf - >>> G = zpk([1],[2, 3], gain=1, display_format='zpk') - >>> G - s - 1 - --------------- - (s - 2) (s - 3) + >>> G = ct.zpk([1], [2, 3], gain=1, display_format='zpk') + >>> print(G) # doctest: +SKIP + + s - 1 + --------------- + (s - 2) (s - 3) + """ num, den = zpk2tf(zeros, poles, gain) return TransferFunction(num, den, *args, **kwargs) @@ -1773,14 +1783,14 @@ def ss2tf(*args, **kwargs): Examples -------- - >>> A = [[1., -2], [3, -4]] - >>> B = [[5.], [7]] - >>> C = [[6., 8]] - >>> D = [[9.]] - >>> sys1 = ss2tf(A, B, C, D) - - >>> sys_ss = ss(A, B, C, D) - >>> sys2 = ss2tf(sys_ss) + >>> A = [[-1, -2], [3, -4]] + >>> B = [[5], [6]] + >>> C = [[7, 8]] + >>> D = [[9]] + >>> sys1 = ct.ss2tf(A, B, C, D) + + >>> sys_ss = ct.ss(A, B, C, D) + >>> sys2 = ct.ss2tf(sys_ss) """ diff --git a/doc/Makefile b/doc/Makefile index b2f9eaeed..6e1012343 100644 --- a/doc/Makefile +++ b/doc/Makefile @@ -16,9 +16,10 @@ help: # Rules to create figures FIGS = classes.pdf -classes.pdf: classes.fig; fig2dev -Lpdf $< $@ +classes.pdf: classes.fig + fig2dev -Lpdf $< $@ # Catch-all target: route all unknown targets to Sphinx using the new # "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). -html pdf clean: Makefile $(FIGS) +html pdf clean doctest: Makefile $(FIGS) @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/doc/classes.png b/doc/classes.png new file mode 100644 index 000000000..25724b43f Binary files /dev/null and b/doc/classes.png differ diff --git a/doc/conf.py b/doc/conf.py index e2e420104..7cfa4b4f6 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -56,7 +56,8 @@ extensions = [ 'sphinx.ext.autodoc', 'sphinx.ext.todo', 'sphinx.ext.napoleon', 'sphinx.ext.intersphinx', 'sphinx.ext.imgmath', - 'sphinx.ext.autosummary', 'nbsphinx', 'numpydoc', 'sphinx.ext.linkcode' + 'sphinx.ext.autosummary', 'nbsphinx', 'numpydoc', + 'sphinx.ext.linkcode', 'sphinx.ext.doctest' ] # scan documents for autosummary directives and generate stub pages for each. @@ -269,3 +270,13 @@ def linkcode_resolve(domain, info): author, 'PythonControlLibrary', 'One line description of project.', 'Miscellaneous'), ] + +# -- Options for doctest ---------------------------------------------- + +# Import control as ct +doctest_global_setup = """ +import numpy as np +import control as ct +import control.optimal as obc +import control.flatsys as fs +""" diff --git a/doc/control.rst b/doc/control.rst index 1ddf7f80a..8dc8a09a4 100644 --- a/doc/control.rst +++ b/doc/control.rst @@ -186,6 +186,8 @@ Utility functions and conversions reachable_form reset_defaults sample_system + set_defaults + similarity_transform ss2tf ssdata tf2ss diff --git a/doc/optimal.rst b/doc/optimal.rst index 807b9b9c6..e84cd2d98 100644 --- a/doc/optimal.rst +++ b/doc/optimal.rst @@ -225,18 +225,18 @@ penalizes the state and input using quadratic cost functions:: Q = np.diag([0, 0, 0.1]) # don't turn too sharply R = np.diag([1, 1]) # keep inputs small P = np.diag([1000, 1000, 1000]) # get close to final point - traj_cost = opt.quadratic_cost(vehicle, Q, R, x0=xf, u0=uf) - term_cost = opt.quadratic_cost(vehicle, P, 0, x0=xf) + traj_cost = obc.quadratic_cost(vehicle, Q, R, x0=xf, u0=uf) + term_cost = obc.quadratic_cost(vehicle, P, 0, x0=xf) We also constraint the maximum turning rate to 0.1 radians (about 6 degees) and constrain the velocity to be in the range of 9 m/s to 11 m/s:: - constraints = [ opt.input_range_constraint(vehicle, [8, -0.1], [12, 0.1]) ] + constraints = [ obc.input_range_constraint(vehicle, [8, -0.1], [12, 0.1]) ] Finally, we solve for the optimal inputs:: timepts = np.linspace(0, Tf, 10, endpoint=True) - result = opt.solve_ocp( + result = obc.solve_ocp( vehicle, timepts, x0, traj_cost, constraints, terminal_cost=term_cost, initial_guess=u0)