From 802b47b1103915440fc469c16fed9c230ca8d246 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Sat, 4 Jan 2020 22:39:21 -0800 Subject: [PATCH 1/7] DOC: fix typos, equation numbering in doc/flatsys.rst --- doc/flatsys.rst | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/doc/flatsys.rst b/doc/flatsys.rst index ed65cfd01..f085347a6 100644 --- a/doc/flatsys.rst +++ b/doc/flatsys.rst @@ -27,6 +27,7 @@ and we can write the solutions of the nonlinear system as functions of .. math:: x &= \beta(z, \dot z, \dots, z^{(q)}) \\ u &= \gamma(z, \dot z, \dots, z^{(q)}). + :label: flat2state For a differentially flat system, all of the feasible trajectories for the system can be written as functions of a flat output :math:`z(\cdot)` and @@ -52,7 +53,7 @@ and we see that the initial and final condition in the full state space depends on just the output :math:`z` and its derivatives at the initial and final times. Thus any trajectory for :math:`z` that satisfies these boundary conditions will be a feasible trajectory for the -system, using equation~\eqref{eq:trajgen:flat2state} to determine the +system, using equation :eq:`flat2state` to determine the full state space and input trajectories. In particular, given initial and final conditions on :math:`z` and its @@ -142,7 +143,7 @@ For more general systems, the `FlatSystem` object must be created manually In addition to the flat system descriptionn, a set of basis functions :math:`\phi_i(t)` must be chosen. The `FlatBasis` class is used to represent the basis functions. A polynomial basis function of the form 1, :math:`t`, -:math:`t^2:, ... can be computed using the `PolyBasis` class, which is +:math:`t^2`, ... can be computed using the `PolyBasis` class, which is initialized by passing the desired order of the polynomial basis set: polybasis = control.flatsys.PolyBasis(N) @@ -225,9 +226,9 @@ derived *Feedback Systems* by Astrom and Murray, Example 3.11. To find a trajectory from an initial state :math:`x_0` to a final state :math:`x_\text{f}` in time :math:`T_\text{f}` we solve a point-to-point -trajectory generation problem. We also set the initial and final inputs, whi -ch sets the vehicle velocity :math:`v` and steering wheel angle :math:`\delta` -at the endpoints. +trajectory generation problem. We also set the initial and final inputs, which +sets the vehicle velocity :math:`v` and steering wheel angle :math:`\delta` at +the endpoints. .. code-block:: python From 7b9d1a9d080559abdcca90977fec8d1e753ad5a5 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Fri, 17 Jan 2020 22:06:56 -0800 Subject: [PATCH 2/7] small updates to statefbk.py docstrings --- control/statefbk.py | 35 ++++++++++++++++++----------------- 1 file changed, 18 insertions(+), 17 deletions(-) diff --git a/control/statefbk.py b/control/statefbk.py index c9e01a52f..91408af25 100644 --- a/control/statefbk.py +++ b/control/statefbk.py @@ -134,16 +134,17 @@ def place_varga(A, B, p, dtime=False, alpha=None): Optional Parameters --------------- - dtime: False for continuous time pole placement or True for discrete time. - The default is dtime=False. - alpha: double scalar - If DICO='C', then place_varga will leave the eigenvalues with real - real part less than alpha untouched. - If DICO='D', the place_varga will leave eigenvalues with modulus - less than alpha untouched. + dtime : bool + False for continuous time pole placement or True for discrete time. + The default is dtime=False. - By default (alpha=None), place_varga computes alpha such that all - poles will be placed. + alpha : double scalar + If DICO='C', then place_varga will leave the eigenvalues with real real + part less than alpha untouched. If DICO='D', the place_varga will + leave eigenvalues with modulus less than alpha untouched. + + By default (alpha=None), place_varga computes alpha such that all poles + will be placed. Returns ------- @@ -153,14 +154,13 @@ def place_varga(A, B, p, dtime=False, alpha=None): Algorithm --------- - This function is a wrapper for the slycot function sb01bd, which - implements the pole placement algorithm of Varga [1]. In contrast to - the algorithm used by place(), the Varga algorithm can place - multiple poles at the same location. The placement, however, may not - be as robust. + This function is a wrapper for the slycot function sb01bd, which + implements the pole placement algorithm of Varga [1]. In contrast to the + algorithm used by place(), the Varga algorithm can place multiple poles at + the same location. The placement, however, may not be as robust. - [1] Varga A. "A Schur method for pole assignment." - IEEE Trans. Automatic Control, Vol. AC-26, pp. 517-519, 1981. + [1] Varga A. "A Schur method for pole assignment." IEEE Trans. Automatic + Control, Vol. AC-26, pp. 517-519, 1981. Examples -------- @@ -171,6 +171,7 @@ def place_varga(A, B, p, dtime=False, alpha=None): See Also: -------- place, acker + """ # Make sure that SLICOT is installed @@ -213,7 +214,7 @@ def place_varga(A, B, p, dtime=False, alpha=None): # but does the trick alpha = -2*abs(min(system_eigs.real)) elif dtime and alpha < 0.0: - raise ValueError("Need alpha > 0 when DICO='D'") + raise ValueError("Discrete time systems require alpha > 0") # Call SLICOT routine to place the eigenvalues From 4f49094e4fa8bba9455c243abbd75463d59bfec8 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Sun, 19 Jan 2020 22:59:56 -0800 Subject: [PATCH 3/7] update setup.py contract info --- setup.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/setup.py b/setup.py index ea825d471..ec16d7135 100644 --- a/setup.py +++ b/setup.py @@ -34,10 +34,10 @@ setup( name='control', version=version, - author='Richard Murray', - author_email='murray@cds.caltech.edu', - url='http://python-control.sourceforge.net', - description='Python control systems library', + author='Python Control Developers', + author_email='python-control-developers@lists.sourceforge.net', + url='http://python-control.org', + description='Python Control Systems Library', long_description=long_description, packages=find_packages(), classifiers=[f for f in CLASSIFIERS.split('\n') if f], From 17c1840509cf318225bbe48a576df37517f4a70d Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Sun, 17 May 2020 13:14:11 -0700 Subject: [PATCH 4/7] small fixes to iosys docstrings (from samlaf) --- control/iosys.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/control/iosys.py b/control/iosys.py index 9b5b89bc8..a90b5193c 100644 --- a/control/iosys.py +++ b/control/iosys.py @@ -780,7 +780,7 @@ def __init__(self, syslist, connections=[], inplist=[], outlist=[], The InterconnectedSystem class is used to represent an input/output system that consists of an interconnection between a set of subystems. - The outputs of each subsystem can be summed together to to provide + The outputs of each subsystem can be summed together to provide inputs to other subsystems. The overall system inputs and outputs can be any subset of subsystem inputs and outputs. From 97c4a65252ade16eebe3b87331789813fe9f9831 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Mon, 21 Dec 2020 21:58:33 -0800 Subject: [PATCH 5/7] TRV: fix indentation in docstring --- control/statefbk.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/control/statefbk.py b/control/statefbk.py index 91408af25..dadbe560e 100644 --- a/control/statefbk.py +++ b/control/statefbk.py @@ -139,12 +139,12 @@ def place_varga(A, B, p, dtime=False, alpha=None): The default is dtime=False. alpha : double scalar - If DICO='C', then place_varga will leave the eigenvalues with real real - part less than alpha untouched. If DICO='D', the place_varga will - leave eigenvalues with modulus less than alpha untouched. + If DICO='C', then place_varga will leave the eigenvalues with real + part less than alpha untouched. If DICO='D', the place_varga will + leave eigenvalues with modulus less than alpha untouched. - By default (alpha=None), place_varga computes alpha such that all poles - will be placed. + By default (alpha=None), place_varga computes alpha such that all + poles will be placed. Returns ------- From 3e6347f2f4f18817a822546af13d49c2d63d1d93 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Mon, 21 Dec 2020 22:18:11 -0800 Subject: [PATCH 6/7] update lqr() return type documentation (addresses #418) --- control/statefbk.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/control/statefbk.py b/control/statefbk.py index dadbe560e..b9a06f502 100644 --- a/control/statefbk.py +++ b/control/statefbk.py @@ -374,9 +374,9 @@ def lqr(*args, **keywords): Returns ------- - K: 2D array + K: 2D array (or matrix) State feedback gains - S: 2D array + S: 2D array (or matrix) Solution to Riccati equation E: 1D array Eigenvalues of the closed loop system @@ -390,6 +390,12 @@ def lqr(*args, **keywords): -------- lqe + Notes + ----- + The return type for `K` and `S` depends on the default class set for + state space operations. By default, this is the Numpy `matrix` + class in this release, but this can be reconfigured using the + :func:`~control.use_numpy_matrix` function. """ # Make sure that SLICOT is installed From 6ff59f091421a26f9816aecf82c89c4aae0dc3f7 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Thu, 24 Dec 2020 09:19:55 -0800 Subject: [PATCH 7/7] Added 2D array/matrix note to all functions using _ssmatrix + numpydoc style fixes + PEP8 --- control/config.py | 7 +- control/mateqn.py | 348 +++++++++++++++++++++++++++----------------- control/statefbk.py | 276 ++++++++++++++++++++--------------- 3 files changed, 379 insertions(+), 252 deletions(-) diff --git a/control/config.py b/control/config.py index c06d38b6e..21840231b 100644 --- a/control/config.py +++ b/control/config.py @@ -154,6 +154,11 @@ class and functions. If flat is `False`, then matrices are of the Numpy `matrix` class. Set `warn` to false to omit display of the warning message. + Notes + ----- + 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. """ if flag and warn: warnings.warn("Return type numpy.matrix is soon to be deprecated.", @@ -179,4 +184,4 @@ def use_legacy_defaults(version): third_digit = int(version[4]) use_numpy_matrix(True) # alternatively: set_defaults('statesp', use_numpy_matrix=True) else: - raise ValueError('''version number not recognized. Possible values range from '0.1' to '0.8.4'.''') \ No newline at end of file + raise ValueError('''version number not recognized. Possible values range from '0.1' to '0.8.4'.''') diff --git a/control/mateqn.py b/control/mateqn.py index 87dd00dab..0b129fd9e 100644 --- a/control/mateqn.py +++ b/control/mateqn.py @@ -1,45 +1,42 @@ -""" mateqn.py - -Matrix equation solvers (Lyapunov, Riccati) - -Implementation of the functions lyap, dlyap, care and dare -for solution of Lyapunov and Riccati equations. """ +# mateqn.py - Matrix equation solvers (Lyapunov, Riccati) +# +# Implementation of the functions lyap, dlyap, care and dare +# for solution of Lyapunov and Riccati equations. +# +# Author: Bjorn Olofsson # Python 3 compatibility (needs to go here) from __future__ import print_function -"""Copyright (c) 2011, All rights reserved. +# Copyright (c) 2011, All rights reserved. -Redistribution and use in source and binary forms, with or without -modification, are permitted provided that the following conditions -are met: +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions +# are met: -1. Redistributions of source code must retain the above copyright - notice, this list of conditions and the following disclaimer. +# 1. Redistributions of source code must retain the above copyright +# notice, this list of conditions and the following disclaimer. -2. Redistributions in binary form must reproduce the above copyright - notice, this list of conditions and the following disclaimer in the - documentation and/or other materials provided with the distribution. +# 2. Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in the +# documentation and/or other materials provided with the distribution. -3. Neither the name of the project author nor the names of its - contributors may be used to endorse or promote products derived - from this software without specific prior written permission. +# 3. Neither the name of the project author nor the names of its +# contributors may be used to endorse or promote products derived +# from this software without specific prior written permission. -THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL CALTECH -OR THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, -SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT -LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF -USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND -ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT -OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF -SUCH DAMAGE. - -Author: Bjorn Olofsson -""" +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS +# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL CALTECH +# OR THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF +# USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND +# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +# OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF +# SUCH DAMAGE. from numpy import shape, size, asarray, copy, zeros, eye, dot, \ finfo, inexact, atleast_2d @@ -49,7 +46,10 @@ __all__ = ['lyap', 'dlyap', 'dare', 'care'] -#### Lyapunov equation solvers lyap and dlyap +# +# Lyapunov equation solvers lyap and dlyap +# + def lyap(A, Q, C=None, E=None): """X = lyap(A, Q) solves the continuous-time Lyapunov equation @@ -59,13 +59,13 @@ def lyap(A, Q, C=None, E=None): where A and Q are square matrices of the same dimension. Further, Q must be symmetric. - X = lyap(A,Q,C) solves the Sylvester equation + X = lyap(A, Q, C) solves the Sylvester equation :math:`A X + X Q + C = 0` where A and Q are square matrices. - X = lyap(A,Q,None,E) solves the generalized continuous-time + X = lyap(A, Q, None, E) solves the generalized continuous-time Lyapunov equation :math:`A X E^T + E X A^T + Q = 0` @@ -73,6 +73,24 @@ def lyap(A, Q, C=None, E=None): where Q is a symmetric matrix and A, Q and E are square matrices of the same dimension. + Parameters + ---------- + A : 2D array + Dynamics matrix + C : 2D array, optional + If present, solve the Slyvester equation + E : 2D array, optional + If present, solve the generalized Laypunov equation + + Returns + ------- + Q : 2D array (or matrix) + Solution to the Lyapunov or Sylvester equation + + Notes + ----- + The return type for 2D arrays depends on the default class set for + state space operations. See :func:`~control.use_numpy_matrix`. """ # Make sure we have access to the right slycot routines @@ -128,7 +146,8 @@ def lyap(A, Q, C=None, E=None): # Solve the Lyapunov equation by calling Slycot function sb03md try: - X,scale,sep,ferr,w = sb03md(n,-Q,A,eye(n,n),'C',trana='T') + X, scale, sep, ferr, w = \ + sb03md(n, -Q, A, eye(n, n), 'C', trana='T') except ValueError as ve: if ve.info < 0: e = ValueError(ve.message) @@ -153,13 +172,14 @@ def lyap(A, Q, C=None, E=None): raise ControlArgument("Q must be a quadratic matrix.") if (size(C) > 1 and shape(C)[0] != n) or \ - (size(C) > 1 and shape(C)[1] != m) or \ - (size(C) == 1 and size(A) != 1) or (size(C) == 1 and size(Q) != 1): + (size(C) > 1 and shape(C)[1] != m) or \ + (size(C) == 1 and size(A) != 1) or \ + (size(C) == 1 and size(Q) != 1): raise ControlArgument("C matrix has incompatible dimensions.") # Solve the Sylvester equation by calling the Slycot function sb04md try: - X = sb04md(n,m,A,Q,-C) + X = sb04md(n, m, A, Q, -C) except ValueError as ve: if ve.info < 0: e = ValueError(ve.message) @@ -178,14 +198,14 @@ def lyap(A, Q, C=None, E=None): elif C is None and E is not None: # Check input data for consistency if (size(Q) > 1 and shape(Q)[0] != shape(Q)[1]) or \ - (size(Q) > 1 and shape(Q)[0] != n) or \ - (size(Q) == 1 and n > 1): + (size(Q) > 1 and shape(Q)[0] != n) or \ + (size(Q) == 1 and n > 1): raise ControlArgument("Q must be a square matrix with the same \ dimension as A.") if (size(E) > 1 and shape(E)[0] != shape(E)[1]) or \ - (size(E) > 1 and shape(E)[0] != n) or \ - (size(E) == 1 and n > 1): + (size(E) > 1 and shape(E)[0] != n) or \ + (size(E) == 1 and n > 1): raise ControlArgument("E must be a square matrix with the same \ dimension as A.") @@ -201,8 +221,9 @@ def lyap(A, Q, C=None, E=None): # Solve the generalized Lyapunov equation by calling Slycot # function sg03ad try: - A,E,Q,Z,X,scale,sep,ferr,alphar,alphai,beta = \ - sg03ad('C','B','N','T','L',n,A,E,eye(n,n),eye(n,n),-Q) + A, E, Q, Z, X, scale, sep, ferr, alphar, alphai, beta = \ + sg03ad('C', 'B', 'N', 'T', 'L', n, + A, E, eye(n, n), eye(n, n), -Q) except ValueError as ve: if ve.info < 0 or ve.info > 4: e = ValueError(ve.message) @@ -235,7 +256,7 @@ def lyap(A, Q, C=None, E=None): return _ssmatrix(X) -def dlyap(A,Q,C=None,E=None): +def dlyap(A, Q, C=None, E=None): """ dlyap(A,Q) solves the discrete-time Lyapunov equation :math:`A X A^T - X + Q = 0` @@ -275,27 +296,27 @@ def dlyap(A,Q,C=None,E=None): # Reshape 1-d arrays if len(shape(A)) == 1: - A = A.reshape(1,A.size) + A = A.reshape(1, A.size) if len(shape(Q)) == 1: - Q = Q.reshape(1,Q.size) + Q = Q.reshape(1, Q.size) if C is not None and len(shape(C)) == 1: - C = C.reshape(1,C.size) + C = C.reshape(1, C.size) if E is not None and len(shape(E)) == 1: - E = E.reshape(1,E.size) + E = E.reshape(1, E.size) # Determine main dimensions if size(A) == 1: n = 1 else: - n = size(A,0) + n = size(A, 0) if size(Q) == 1: m = 1 else: - m = size(Q,0) + m = size(Q, 0) # Solve standard Lyapunov equation if C is None and E is None: @@ -315,7 +336,8 @@ def dlyap(A,Q,C=None,E=None): # Solve the Lyapunov equation by calling the Slycot function sb03md try: - X,scale,sep,ferr,w = sb03md(n,-Q,A,eye(n,n),'D',trana='T') + X, scale, sep, ferr, w = \ + sb03md(n, -Q, A, eye(n, n), 'D', trana='T') except ValueError as ve: if ve.info < 0: e = ValueError(ve.message) @@ -336,13 +358,13 @@ def dlyap(A,Q,C=None,E=None): raise ControlArgument("Q must be a quadratic matrix") if (size(C) > 1 and shape(C)[0] != n) or \ - (size(C) > 1 and shape(C)[1] != m) or \ - (size(C) == 1 and size(A) != 1) or (size(C) == 1 and size(Q) != 1): + (size(C) > 1 and shape(C)[1] != m) or \ + (size(C) == 1 and size(A) != 1) or (size(C) == 1 and size(Q) != 1): raise ControlArgument("C matrix has incompatible dimensions") # Solve the Sylvester equation by calling Slycot function sb04qd try: - X = sb04qd(n,m,-A,asarray(Q).T,C) + X = sb04qd(n, m, -A, asarray(Q).T, C) except ValueError as ve: if ve.info < 0: e = ValueError(ve.message) @@ -361,14 +383,14 @@ def dlyap(A,Q,C=None,E=None): elif C is None and E is not None: # Check input data for consistency if (size(Q) > 1 and shape(Q)[0] != shape(Q)[1]) or \ - (size(Q) > 1 and shape(Q)[0] != n) or \ - (size(Q) == 1 and n > 1): + (size(Q) > 1 and shape(Q)[0] != n) or \ + (size(Q) == 1 and n > 1): raise ControlArgument("Q must be a square matrix with the same \ dimension as A.") if (size(E) > 1 and shape(E)[0] != shape(E)[1]) or \ - (size(E) > 1 and shape(E)[0] != n) or \ - (size(E) == 1 and n > 1): + (size(E) > 1 and shape(E)[0] != n) or \ + (size(E) == 1 and n > 1): raise ControlArgument("E must be a square matrix with the same \ dimension as A.") @@ -378,8 +400,9 @@ def dlyap(A,Q,C=None,E=None): # Solve the generalized Lyapunov equation by calling Slycot # function sg03ad try: - A,E,Q,Z,X,scale,sep,ferr,alphar,alphai,beta = \ - sg03ad('D','B','N','T','L',n,A,E,eye(n,n),eye(n,n),-Q) + A, E, Q, Z, X, scale, sep, ferr, alphar, alphai, beta = \ + sg03ad('D', 'B', 'N', 'T', 'L', n, + A, E, eye(n, n), eye(n, n), -Q) except ValueError as ve: if ve.info < 0 or ve.info > 4: e = ValueError(ve.message) @@ -412,10 +435,14 @@ def dlyap(A,Q,C=None,E=None): return _ssmatrix(X) -#### Riccati equation solvers care and dare +# +# Riccati equation solvers care and dare +# + + def care(A, B, Q, R=None, S=None, E=None, stabilizing=True): - """ (X,L,G) = care(A,B,Q,R=None) solves the continuous-time algebraic Riccati - equation + """(X, L, G) = care(A, B, Q, R=None) solves the continuous-time + algebraic Riccati equation :math:`A^T X + X A - X B R^{-1} B^T X + Q = 0` @@ -425,16 +452,39 @@ def care(A, B, Q, R=None, S=None, E=None, stabilizing=True): matrix G = B^T X and the closed loop eigenvalues L, i.e., the eigenvalues of A - B G. - (X,L,G) = care(A,B,Q,R,S,E) solves the generalized continuous-time - algebraic Riccati equation + (X, L, G) = care(A, B, Q, R, S, E) solves the generalized + continuous-time algebraic Riccati equation :math:`A^T X E + E^T X A - (E^T X B + S) R^{-1} (B^T X E + S^T) + Q = 0` - where A, Q and E are square matrices of the same - dimension. Further, Q and R are symmetric matrices. If R is None, - it is set to the identity matrix. The function returns the - solution X, the gain matrix G = R^-1 (B^T X E + S^T) and the - closed loop eigenvalues L, i.e., the eigenvalues of A - B G , E.""" + where A, Q and E are square matrices of the same dimension. Further, Q + and R are symmetric matrices. If R is None, it is set to the identity + matrix. The function returns the solution X, the gain matrix G = R^-1 + (B^T X E + S^T) and the closed loop eigenvalues L, i.e., the eigenvalues + of A - B G , E. + + Parameters + ---------- + A, B, Q : 2D arrays + Input matrices for the Riccati equation + R, S, E : 2D arrays, optional + Input matrices for generalized Riccati equation + + Returns + ------- + X : 2D array (or matrix) + Solution to the Ricatti equation + L : 1D array + Closed loop eigenvalues + G : 2D array (or matrix) + Gain matrix + + Notes + ----- + The return type for 2D arrays depends on the default class set for + state space operations. See :func:`~control.use_numpy_matrix`. + + """ # Make sure we can import required slycot routine try: @@ -455,35 +505,35 @@ def care(A, B, Q, R=None, S=None, E=None, stabilizing=True): # Reshape 1-d arrays if len(shape(A)) == 1: - A = A.reshape(1,A.size) + A = A.reshape(1, A.size) if len(shape(B)) == 1: - B = B.reshape(1,B.size) + B = B.reshape(1, B.size) if len(shape(Q)) == 1: - Q = Q.reshape(1,Q.size) + Q = Q.reshape(1, Q.size) if R is not None and len(shape(R)) == 1: - R = R.reshape(1,R.size) + R = R.reshape(1, R.size) if S is not None and len(shape(S)) == 1: - S = S.reshape(1,S.size) + S = S.reshape(1, S.size) if E is not None and len(shape(E)) == 1: - E = E.reshape(1,E.size) + E = E.reshape(1, E.size) # Determine main dimensions if size(A) == 1: n = 1 else: - n = size(A,0) + n = size(A, 0) if size(B) == 1: m = 1 else: - m = size(B,1) + m = size(B, 1) if R is None: - R = eye(m,m) + R = eye(m, m) # Solve the standard algebraic Riccati equation if S is None and E is None: @@ -492,13 +542,13 @@ def care(A, B, Q, R=None, S=None, E=None, stabilizing=True): raise ControlArgument("A must be a quadratic matrix.") if (size(Q) > 1 and shape(Q)[0] != shape(Q)[1]) or \ - (size(Q) > 1 and shape(Q)[0] != n) or \ - size(Q) == 1 and n > 1: + (size(Q) > 1 and shape(Q)[0] != n) or \ + size(Q) == 1 and n > 1: raise ControlArgument("Q must be a quadratic matrix of the same \ dimension as A.") if (size(B) > 1 and shape(B)[0] != n) or \ - size(B) == 1 and n > 1: + size(B) == 1 and n > 1: raise ControlArgument("Incompatible dimensions of B matrix.") if not _is_symmetric(Q): @@ -514,7 +564,7 @@ def care(A, B, Q, R=None, S=None, E=None, stabilizing=True): # Solve the standard algebraic Riccati equation by calling Slycot # functions sb02mt and sb02md try: - A_b,B_b,Q_b,R_b,L_b,ipiv,oufact,G = sb02mt(n,m,B,R) + A_b, B_b, Q_b, R_b, L_b, ipiv, oufact, G = sb02mt(n, m, B, R) except ValueError as ve: if ve.info < 0: e = ValueError(ve.message) @@ -568,7 +618,7 @@ def care(A, B, Q, R=None, S=None, E=None, stabilizing=True): # Return the solution X, the closed-loop eigenvalues L and # the gain matrix G - return (_ssmatrix(X) , w[:n] , _ssmatrix(G)) + return (_ssmatrix(X), w[:n], _ssmatrix(G)) # Solve the generalized algebraic Riccati equation elif S is not None and E is not None: @@ -577,31 +627,31 @@ def care(A, B, Q, R=None, S=None, E=None, stabilizing=True): raise ControlArgument("A must be a quadratic matrix.") if (size(Q) > 1 and shape(Q)[0] != shape(Q)[1]) or \ - (size(Q) > 1 and shape(Q)[0] != n) or \ - size(Q) == 1 and n > 1: + (size(Q) > 1 and shape(Q)[0] != n) or \ + size(Q) == 1 and n > 1: raise ControlArgument("Q must be a quadratic matrix of the same \ dimension as A.") if (size(B) > 1 and shape(B)[0] != n) or \ - size(B) == 1 and n > 1: + size(B) == 1 and n > 1: raise ControlArgument("Incompatible dimensions of B matrix.") if (size(E) > 1 and shape(E)[0] != shape(E)[1]) or \ - (size(E) > 1 and shape(E)[0] != n) or \ - size(E) == 1 and n > 1: + (size(E) > 1 and shape(E)[0] != n) or \ + size(E) == 1 and n > 1: raise ControlArgument("E must be a quadratic matrix of the same \ dimension as A.") if (size(R) > 1 and shape(R)[0] != shape(R)[1]) or \ - (size(R) > 1 and shape(R)[0] != m) or \ - size(R) == 1 and m > 1: + (size(R) > 1 and shape(R)[0] != m) or \ + size(R) == 1 and m > 1: raise ControlArgument("R must be a quadratic matrix of the same \ dimension as the number of columns in the B matrix.") if (size(S) > 1 and shape(S)[0] != n) or \ - (size(S) > 1 and shape(S)[1] != m) or \ - size(S) == 1 and n > 1 or \ - size(S) == 1 and m > 1: + (size(S) > 1 and shape(S)[1] != m) or \ + size(S) == 1 and n > 1 or \ + size(S) == 1 and m > 1: raise ControlArgument("Incompatible dimensions of S matrix.") if not _is_symmetric(Q): @@ -624,7 +674,8 @@ def care(A, B, Q, R=None, S=None, E=None, stabilizing=True): else: sort = 'U' rcondu, X, alfar, alfai, beta, S_o, T, U, iwarn = \ - sg02ad('C', 'B', 'N', 'U', 'N', 'N', sort, 'R', n, m, 0, A, E, B, Q, R, S) + sg02ad('C', 'B', 'N', 'U', 'N', 'N', sort, + 'R', n, m, 0, A, E, B, Q, R, S) except ValueError as ve: if ve.info < 0 or ve.info > 7: e = ValueError(ve.message) @@ -662,14 +713,14 @@ def care(A, B, Q, R=None, S=None, E=None, stabilizing=True): raise e # Calculate the closed-loop eigenvalues L - L = zeros((n,1)) + L = zeros((n, 1)) L.dtype = 'complex64' for i in range(n): L[i] = (alfar[i] + alfai[i]*1j)/beta[i] # Calculate the gain matrix G if size(R_b) == 1: - G = dot(1/(R_b), dot(asarray(B_b).T, dot(X,E_b)) + asarray(S_b).T) + G = dot(1/(R_b), dot(asarray(B_b).T, dot(X, E_b)) + asarray(S_b).T) else: G = solve(R_b, dot(asarray(B_b).T, dot(X, E_b)) + asarray(S_b).T) @@ -681,8 +732,9 @@ def care(A, B, Q, R=None, S=None, E=None, stabilizing=True): else: raise ControlArgument("Invalid set of input parameters.") + def dare(A, B, Q, R, S=None, E=None, stabilizing=True): - """ (X,L,G) = dare(A,B,Q,R) solves the discrete-time algebraic Riccati + """(X, L, G) = dare(A, B, Q, R) solves the discrete-time algebraic Riccati equation :math:`A^T X A - X - A^T X B (B^T X B + R)^{-1} B^T X A + Q = 0` @@ -692,8 +744,8 @@ def dare(A, B, Q, R, S=None, E=None, stabilizing=True): matrix G = (B^T X B + R)^-1 B^T X A and the closed loop eigenvalues L, i.e., the eigenvalues of A - B G. - (X,L,G) = dare(A,B,Q,R,S,E) solves the generalized discrete-time algebraic - Riccati equation + (X, L, G) = dare(A, B, Q, R, S, E) solves the generalized discrete-time + algebraic Riccati equation :math:`A^T X A - E^T X E - (A^T X B + S) (B^T X B + R)^{-1} (B^T X A + S^T) + Q = 0` @@ -701,6 +753,28 @@ def dare(A, B, Q, R, S=None, E=None, stabilizing=True): R are symmetric matrices. The function returns the solution X, the gain matrix :math:`G = (B^T X B + R)^{-1} (B^T X A + S^T)` and the closed loop eigenvalues L, i.e., the eigenvalues of A - B G , E. + + Parameters + ---------- + A, B, Q : 2D arrays + Input matrices for the Riccati equation + R, S, E : 2D arrays, optional + Input matrices for generalized Riccati equation + + Returns + ------- + X : 2D array (or matrix) + Solution to the Ricatti equation + L : 1D array + Closed loop eigenvalues + G : 2D array (or matrix) + Gain matrix + + Notes + ----- + The return type for 2D arrays depends on the default class set for + state space operations. See :func:`~control.use_numpy_matrix`. + """ if S is not None or E is not None or not stabilizing: return dare_old(A, B, Q, R, S, E, stabilizing) @@ -712,6 +786,7 @@ def dare(A, B, Q, R, S=None, E=None, stabilizing=True): L = eigvals(A - B.dot(G)) return _ssmatrix(X), L, _ssmatrix(G) + def dare_old(A, B, Q, R, S=None, E=None, stabilizing=True): # Make sure we can import required slycot routine try: @@ -732,33 +807,33 @@ def dare_old(A, B, Q, R, S=None, E=None, stabilizing=True): # Reshape 1-d arrays if len(shape(A)) == 1: - A = A.reshape(1,A.size) + A = A.reshape(1, A.size) if len(shape(B)) == 1: - B = B.reshape(1,B.size) + B = B.reshape(1, B.size) if len(shape(Q)) == 1: - Q = Q.reshape(1,Q.size) + Q = Q.reshape(1, Q.size) if R is not None and len(shape(R)) == 1: - R = R.reshape(1,R.size) + R = R.reshape(1, R.size) if S is not None and len(shape(S)) == 1: - S = S.reshape(1,S.size) + S = S.reshape(1, S.size) if E is not None and len(shape(E)) == 1: - E = E.reshape(1,E.size) + E = E.reshape(1, E.size) # Determine main dimensions if size(A) == 1: n = 1 else: - n = size(A,0) + n = size(A, 0) if size(B) == 1: m = 1 else: - m = size(B,1) + m = size(B, 1) # Solve the standard algebraic Riccati equation if S is None and E is None: @@ -767,13 +842,13 @@ def dare_old(A, B, Q, R, S=None, E=None, stabilizing=True): raise ControlArgument("A must be a quadratic matrix.") if (size(Q) > 1 and shape(Q)[0] != shape(Q)[1]) or \ - (size(Q) > 1 and shape(Q)[0] != n) or \ - size(Q) == 1 and n > 1: + (size(Q) > 1 and shape(Q)[0] != n) or \ + size(Q) == 1 and n > 1: raise ControlArgument("Q must be a quadratic matrix of the same \ dimension as A.") if (size(B) > 1 and shape(B)[0] != n) or \ - size(B) == 1 and n > 1: + size(B) == 1 and n > 1: raise ControlArgument("Incompatible dimensions of B matrix.") if not _is_symmetric(Q): @@ -790,7 +865,7 @@ def dare_old(A, B, Q, R, S=None, E=None, stabilizing=True): # Solve the standard algebraic Riccati equation by calling Slycot # functions sb02mt and sb02md try: - A_b,B_b,Q_b,R_b,L_b,ipiv,oufact,G = sb02mt(n,m,B,R) + A_b, B_b, Q_b, R_b, L_b, ipiv, oufact, G = sb02mt(n, m, B, R) except ValueError as ve: if ve.info < 0: e = ValueError(ve.message) @@ -839,15 +914,15 @@ def dare_old(A, B, Q, R, S=None, E=None, stabilizing=True): # Calculate the gain matrix G if size(R_b) == 1: - G = dot(1/(dot(asarray(B_ba).T, dot(X, B_ba)) + R_ba), \ - dot(asarray(B_ba).T, dot(X, A_ba))) + G = dot(1/(dot(asarray(B_ba).T, dot(X, B_ba)) + R_ba), + dot(asarray(B_ba).T, dot(X, A_ba))) else: - G = solve(dot(asarray(B_ba).T, dot(X, B_ba)) + R_ba, \ - dot(asarray(B_ba).T, dot(X, A_ba))) + G = solve(dot(asarray(B_ba).T, dot(X, B_ba)) + R_ba, + dot(asarray(B_ba).T, dot(X, A_ba))) # Return the solution X, the closed-loop eigenvalues L and # the gain matrix G - return (_ssmatrix(X) , w[:n], _ssmatrix(G)) + return (_ssmatrix(X), w[:n], _ssmatrix(G)) # Solve the generalized algebraic Riccati equation elif S is not None and E is not None: @@ -856,31 +931,31 @@ def dare_old(A, B, Q, R, S=None, E=None, stabilizing=True): raise ControlArgument("A must be a quadratic matrix.") if (size(Q) > 1 and shape(Q)[0] != shape(Q)[1]) or \ - (size(Q) > 1 and shape(Q)[0] != n) or \ - size(Q) == 1 and n > 1: + (size(Q) > 1 and shape(Q)[0] != n) or \ + size(Q) == 1 and n > 1: raise ControlArgument("Q must be a quadratic matrix of the same \ dimension as A.") if (size(B) > 1 and shape(B)[0] != n) or \ - size(B) == 1 and n > 1: + size(B) == 1 and n > 1: raise ControlArgument("Incompatible dimensions of B matrix.") if (size(E) > 1 and shape(E)[0] != shape(E)[1]) or \ - (size(E) > 1 and shape(E)[0] != n) or \ - size(E) == 1 and n > 1: + (size(E) > 1 and shape(E)[0] != n) or \ + size(E) == 1 and n > 1: raise ControlArgument("E must be a quadratic matrix of the same \ dimension as A.") if (size(R) > 1 and shape(R)[0] != shape(R)[1]) or \ - (size(R) > 1 and shape(R)[0] != m) or \ - size(R) == 1 and m > 1: + (size(R) > 1 and shape(R)[0] != m) or \ + size(R) == 1 and m > 1: raise ControlArgument("R must be a quadratic matrix of the same \ dimension as the number of columns in the B matrix.") if (size(S) > 1 and shape(S)[0] != n) or \ - (size(S) > 1 and shape(S)[1] != m) or \ - size(S) == 1 and n > 1 or \ - size(S) == 1 and m > 1: + (size(S) > 1 and shape(S)[1] != m) or \ + size(S) == 1 and n > 1 or \ + size(S) == 1 and m > 1: raise ControlArgument("Incompatible dimensions of S matrix.") if not _is_symmetric(Q): @@ -904,7 +979,8 @@ def dare_old(A, B, Q, R, S=None, E=None, stabilizing=True): else: sort = 'U' rcondu, X, alfar, alfai, beta, S_o, T, U, iwarn = \ - sg02ad('D', 'B', 'N', 'U', 'N', 'N', sort, 'R', n, m, 0, A, E, B, Q, R, S) + sg02ad('D', 'B', 'N', 'U', 'N', 'N', sort, + 'R', n, m, 0, A, E, B, Q, R, S) except ValueError as ve: if ve.info < 0 or ve.info > 7: e = ValueError(ve.message) @@ -941,18 +1017,18 @@ def dare_old(A, B, Q, R, S=None, E=None, stabilizing=True): e.info = ve.info raise e - L = zeros((n,1)) + L = zeros((n, 1)) L.dtype = 'complex64' for i in range(n): L[i] = (alfar[i] + alfai[i]*1j)/beta[i] # Calculate the gain matrix G if size(R_b) == 1: - G = dot(1/(dot(asarray(B_b).T, dot(X,B_b)) + R_b), \ - dot(asarray(B_b).T, dot(X,A_b)) + asarray(S_b).T) + G = dot(1/(dot(asarray(B_b).T, dot(X, B_b)) + R_b), + dot(asarray(B_b).T, dot(X, A_b)) + asarray(S_b).T) else: - G = solve(dot(asarray(B_b).T, dot(X,B_b)) + R_b, \ - dot(asarray(B_b).T, dot(X,A_b)) + asarray(S_b).T) + G = solve(dot(asarray(B_b).T, dot(X, B_b)) + R_b, + dot(asarray(B_b).T, dot(X, A_b)) + asarray(S_b).T) # Return the solution X, the closed-loop eigenvalues L and # the gain matrix G diff --git a/control/statefbk.py b/control/statefbk.py index b9a06f502..d07410bfa 100644 --- a/control/statefbk.py +++ b/control/statefbk.py @@ -47,7 +47,8 @@ from .statesp import _ssmatrix from .exception import ControlSlycot, ControlArgument, ControlDimension -__all__ = ['ctrb', 'obsv', 'gram', 'place', 'place_varga', 'lqr', 'lqe', 'acker'] +__all__ = ['ctrb', 'obsv', 'gram', 'place', 'place_varga', 'lqr', 'lqe', + 'acker'] # Pole placement @@ -58,19 +59,18 @@ def place(A, B, p): Parameters ---------- - A : 2-d array + A : 2D array Dynamics matrix - B : 2-d array + B : 2D array Input matrix - p : 1-d list + p : 1D list Desired eigenvalue locations Returns ------- - K : 2-d array + K : 2D array (or matrix) Gain such that A - B K has eigenvalues given in p - Notes ----- Algorithm @@ -83,6 +83,9 @@ def place(A, B, p): The algorithm will not place poles at the same location more than rank(B) times. + The return type for 2D arrays depends on the default class set for + state space operations. See :func:`~control.use_numpy_matrix`. + References ---------- .. [1] A.L. Tits and Y. Yang, "Globally convergent algorithms for robust @@ -98,6 +101,11 @@ def place(A, B, p): See Also -------- place_varga, acker + + Notes + ----- + The return type for 2D arrays depends on the default class set for + state space operations. See :func:`~control.use_numpy_matrix`. """ from scipy.signal import place_poles @@ -125,11 +133,11 @@ def place_varga(A, B, p, dtime=False, alpha=None): Required Parameters ---------- - A : 2-d array + A : 2D array Dynamics matrix - B : 2-d array + B : 2D array Input matrix - p : 1-d list + p : 1D list Desired eigenvalue locations Optional Parameters @@ -139,19 +147,19 @@ def place_varga(A, B, p, dtime=False, alpha=None): The default is dtime=False. alpha : double scalar - If DICO='C', then place_varga will leave the eigenvalues with real - part less than alpha untouched. If DICO='D', the place_varga will - leave eigenvalues with modulus less than alpha untouched. + If `dtime` is false then place_varga will leave the eigenvalues with + real part less than alpha untouched. If `dtime` is true then + place_varga will leave eigenvalues with modulus less than alpha + untouched. - By default (alpha=None), place_varga computes alpha such that all + By default (alpha=None), place_varga computes alpha such that all poles will be placed. Returns ------- - K : 2D array + K : 2D array (or matrix) Gain such that A - B K has eigenvalues given in p. - Algorithm --------- This function is a wrapper for the slycot function sb01bd, which @@ -162,6 +170,11 @@ def place_varga(A, B, p, dtime=False, alpha=None): [1] Varga A. "A Schur method for pole assignment." IEEE Trans. Automatic Control, Vol. AC-26, pp. 517-519, 1981. + Notes + ----- + The return type for 2D arrays depends on the default class set for + state space operations. See :func:`~control.use_numpy_matrix`. + Examples -------- >>> A = [[-1, -1], [0, 1]] @@ -183,8 +196,7 @@ def place_varga(A, B, p, dtime=False, alpha=None): # Convert the system inputs to NumPy arrays A_mat = np.array(A) B_mat = np.array(B) - if (A_mat.shape[0] != A_mat.shape[1] or - A_mat.shape[0] != B_mat.shape[0]): + if (A_mat.shape[0] != A_mat.shape[1] or A_mat.shape[0] != B_mat.shape[0]): raise ControlDimension("matrix dimensions are incorrect") # Compute the system eigenvalues and convert poles to numpy array @@ -216,15 +228,15 @@ def place_varga(A, B, p, dtime=False, alpha=None): elif dtime and alpha < 0.0: raise ValueError("Discrete time systems require alpha > 0") - # Call SLICOT routine to place the eigenvalues - A_z,w,nfp,nap,nup,F,Z = \ + A_z, w, nfp, nap, nup, F, Z = \ sb01bd(B_mat.shape[0], B_mat.shape[1], len(placed_eigs), alpha, A_mat, B_mat, placed_eigs, DICO) # Return the gain matrix, with MATLAB gain convention return _ssmatrix(-F) + # contributed by Sawyer B. Fuller def lqe(A, G, C, QN, RN, NN=None): """lqe(A, G, C, QN, RN, [, N]) @@ -246,33 +258,42 @@ def lqe(A, G, C, QN, RN, NN=None): .. math:: x_e = A x_e + B u + L(y - C x_e - D u) - produces a state estimate that x_e that minimizes the expected squared error - using the sensor measurements y. The noise cross-correlation `NN` is set to - zero when omitted. + produces a state estimate that x_e that minimizes the expected squared + error using the sensor measurements y. The noise cross-correlation `NN` + is set to zero when omitted. Parameters ---------- - A, G: 2-d array + A, G : 2D array Dynamics and noise input matrices - QN, RN: 2-d array + QN, RN : 2D array Process and sensor noise covariance matrices - NN: 2-d array, optional + NN : 2D array, optional Cross covariance matrix Returns ------- - L: 2D array + L : 2D array (or matrix) Kalman estimator gain - P: 2D array + P : 2D array (or matrix) Solution to Riccati equation .. math:: A P + P A^T - (P C^T + G N) R^{-1} (C P + N^T G^T) + G Q G^T = 0 - E: 1D array + E : 2D array (or matrix) Eigenvalues of estimator poles eig(A - L C) + Notes + ----- + The return type for 2D arrays depends on the default class set for + state space operations. See :func:`~control.use_numpy_matrix`. + + The return type for `E` differs from the equivalent return values in the + :func:`~control.lqr`, :func:`~control.care`, and other similar + functions. The return type will be changed to a 1D array in a future + release. Examples -------- @@ -282,18 +303,19 @@ def lqe(A, G, C, QN, RN, NN=None): See Also -------- lqr + """ # TODO: incorporate cross-covariance NN, something like this, # which doesn't work for some reason - #if NN is None: + # if NN is None: # NN = np.zeros(QN.size(0),RN.size(1)) - #NG = G @ NN + # NG = G @ NN - #LT, P, E = lqr(A.T, C.T, G @ QN @ G.T, RN) - #P, E, LT = care(A.T, C.T, G @ QN @ G.T, RN) + # LT, P, E = lqr(A.T, C.T, G @ QN @ G.T, RN) + # P, E, LT = care(A.T, C.T, G @ QN @ G.T, RN) A, G, C = np.array(A, ndmin=2), np.array(G, ndmin=2), np.array(C, ndmin=2) - QN, RN = np.array(QN, ndmin=2), np.array(RN, ndmin=2) + QN, RN = np.array(QN, ndmin=2), np.array(RN, ndmin=2) P, E, LT = care(A.T, C.T, np.dot(np.dot(G, QN), G.T), RN) return _ssmatrix(LT.T), _ssmatrix(P), _ssmatrix(E) @@ -307,16 +329,20 @@ def acker(A, B, poles): Parameters ---------- - A, B : 2-d arrays + A, B : 2D arrays State and input matrix of the system - poles: 1-d list + poles : 1D list Desired eigenvalue locations Returns ------- - K: matrix + K : 2D array (or matrix) Gains such that A - B K has given eigenvalues + Notes + ----- + The return type for 2D arrays depends on the default class set for + state space operations. See :func:`~control.use_numpy_matrix`. """ # Convert the inputs to matrices a = _ssmatrix(A) @@ -334,13 +360,14 @@ def acker(A, B, poles): # TODO: compute pmat using Horner's method (O(n) instead of O(n^2)) n = np.size(p) pmat = p[n-1] * np.linalg.matrix_power(a, 0) - for i in np.arange(1,n): + for i in np.arange(1, n): pmat = pmat + np.dot(p[n-i-1], np.linalg.matrix_power(a, i)) K = np.linalg.solve(ct, pmat) K = K[-1][:] # Extract the last row return _ssmatrix(K) + def lqr(*args, **keywords): """lqr(A, B, Q, R[, N]) @@ -363,39 +390,37 @@ def lqr(*args, **keywords): Parameters ---------- - A, B: 2-d array + A, B : 2D array Dynamics and input matrices - sys: LTI (StateSpace or TransferFunction) + sys : LTI (StateSpace or TransferFunction) Linear I/O system - Q, R: 2-d array + Q, R : 2D array State and input weight matrices - N: 2-d array, optional + N : 2D array, optional Cross weight matrix Returns ------- - K: 2D array (or matrix) + K : 2D array (or matrix) State feedback gains - S: 2D array (or matrix) + S : 2D array (or matrix) Solution to Riccati equation - E: 1D array + E : 1D array Eigenvalues of the closed loop system - Examples - -------- - >>> K, S, E = lqr(sys, Q, R, [N]) - >>> K, S, E = lqr(A, B, Q, R, [N]) - See Also -------- lqe Notes ----- - The return type for `K` and `S` depends on the default class set for - state space operations. By default, this is the Numpy `matrix` - class in this release, but this can be reconfigured using the - :func:`~control.use_numpy_matrix` function. + The return type for 2D arrays depends on the default class set for + state space operations. See :func:`~control.use_numpy_matrix`. + + Examples + -------- + >>> K, S, E = lqr(sys, Q, R, [N]) + >>> K, S, E = lqr(A, B, Q, R, [N]) """ # Make sure that SLICOT is installed @@ -416,26 +441,26 @@ class in this release, but this can be reconfigured using the try: # If this works, we were (probably) passed a system as the # first argument; extract A and B - A = np.array(args[0].A, ndmin=2, dtype=float); - B = np.array(args[0].B, ndmin=2, dtype=float); - index = 1; + A = np.array(args[0].A, ndmin=2, dtype=float) + B = np.array(args[0].B, ndmin=2, dtype=float) + index = 1 except AttributeError: # Arguments should be A and B matrices - A = np.array(args[0], ndmin=2, dtype=float); - B = np.array(args[1], ndmin=2, dtype=float); - index = 2; + A = np.array(args[0], ndmin=2, dtype=float) + B = np.array(args[1], ndmin=2, dtype=float) + index = 2 # Get the weighting matrices (converting to matrices, if needed) - Q = np.array(args[index], ndmin=2, dtype=float); - R = np.array(args[index+1], ndmin=2, dtype=float); + Q = np.array(args[index], ndmin=2, dtype=float) + R = np.array(args[index+1], ndmin=2, dtype=float) if (len(args) > index + 2): - N = np.array(args[index+2], ndmin=2, dtype=float); + N = np.array(args[index+2], ndmin=2, dtype=float) else: - N = np.zeros((Q.shape[0], R.shape[1])); + N = np.zeros((Q.shape[0], R.shape[1])) # Check dimensions for consistency - nstates = B.shape[0]; - ninputs = B.shape[1]; + nstates = B.shape[0] + ninputs = B.shape[1] if (A.shape[0] != nstates or A.shape[1] != nstates): raise ControlDimension("inconsistent system dimensions") @@ -445,33 +470,39 @@ class in this release, but this can be reconfigured using the raise ControlDimension("incorrect weighting matrix dimensions") # Compute the G matrix required by SB02MD - A_b,B_b,Q_b,R_b,L_b,ipiv,oufact,G = \ - sb02mt(nstates, ninputs, B, R, A, Q, N, jobl='N'); + A_b, B_b, Q_b, R_b, L_b, ipiv, oufact, G = \ + sb02mt(nstates, ninputs, B, R, A, Q, N, jobl='N') # Call the SLICOT function - X,rcond,w,S,U,A_inv = sb02md(nstates, A_b, G, Q_b, 'C') + X, rcond, w, S, U, A_inv = sb02md(nstates, A_b, G, Q_b, 'C') # Now compute the return value # We assume that R is positive definite and, hence, invertible - K = np.linalg.solve(R, np.dot(B.T, X) + N.T); - S = X; - E = w[0:nstates]; + K = np.linalg.solve(R, np.dot(B.T, X) + N.T) + S = X + E = w[0:nstates] return _ssmatrix(K), _ssmatrix(S), E + def ctrb(A, B): """Controllabilty matrix Parameters ---------- - A, B: array_like or string + A, B : array_like or string Dynamics and input matrix of the system Returns ------- - C: matrix + C : 2D array (or matrix) Controllability matrix + Notes + ----- + The return type for 2D arrays depends on the default class set for + state space operations. See :func:`~control.use_numpy_matrix`. + Examples -------- >>> C = ctrb(A, B) @@ -484,28 +515,34 @@ def ctrb(A, B): n = np.shape(amat)[0] # Construct the controllability matrix - ctrb = np.hstack([bmat] + [np.dot(np.linalg.matrix_power(amat, i), bmat) - for i in range(1, n)]) + ctrb = np.hstack( + [bmat] + [np.dot(np.linalg.matrix_power(amat, i), bmat) + for i in range(1, n)]) return _ssmatrix(ctrb) + def obsv(A, C): """Observability matrix Parameters ---------- - A, C: array_like or string + A, C : array_like or string Dynamics and output matrix of the system Returns ------- - O: matrix + O : 2D array (or matrix) Observability matrix + Notes + ----- + The return type for 2D arrays depends on the default class set for + state space operations. See :func:`~control.use_numpy_matrix`. + Examples -------- >>> O = obsv(A, C) - - """ + """ # Convert input parameters to matrices (if they aren't already) amat = _ssmatrix(A) @@ -517,21 +554,22 @@ def obsv(A, C): for i in range(1, n)]) return _ssmatrix(obsv) -def gram(sys,type): + +def gram(sys, type): """Gramian (controllability or observability) Parameters ---------- - sys: StateSpace - State-space system to compute Gramian for - type: String - Type of desired computation. - `type` is either 'c' (controllability) or 'o' (observability). To compute the - Cholesky factors of gramians use 'cf' (controllability) or 'of' (observability) + sys : StateSpace + System description + type : String + Type of desired computation. `type` is either 'c' (controllability) + or 'o' (observability). To compute the Cholesky factors of Gramians + use 'cf' (controllability) or 'of' (observability) Returns ------- - gram: array + gram : 2D array (or matrix) Gramian of system Raises @@ -545,22 +583,27 @@ def gram(sys,type): if slycot routine sb03md cannot be found if slycot routine sb03od cannot be found + Notes + ----- + The return type for 2D arrays depends on the default class set for + state space operations. See :func:`~control.use_numpy_matrix`. + 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 + >>> Wc = gram(sys, 'c') + >>> Wo = gram(sys, 'o') + >>> Rc = gram(sys, 'cf'), where Wc = Rc' * Rc + >>> Ro = gram(sys, 'of'), where Wo = Ro' * Ro """ - #Check for ss system object - if not isinstance(sys,statesp.StateSpace): + # Check for ss system object + if not isinstance(sys, statesp.StateSpace): raise ValueError("System must be StateSpace!") if type not in ['c', 'o', 'cf', 'of']: raise ValueError("That type is not supported!") - #TODO: Check for continous or discrete, only continuous supported right now + # TODO: Check for continous or discrete, only continuous supported for now # if isCont(): # dico = 'C' # elif isDisc(): @@ -568,50 +611,53 @@ def gram(sys,type): # else: dico = 'C' - #TODO: Check system is stable, perhaps a utility in ctrlutil.py - # or a method of the StateSpace class? + # TODO: Check system is stable, perhaps a utility in ctrlutil.py + # or a method of the StateSpace class? if np.any(np.linalg.eigvals(sys.A).real >= 0.0): raise ValueError("Oops, the system is unstable!") - if type=='c' or type=='o': - #Compute Gramian by the Slycot routine sb03md - #make sure Slycot is installed + if type == 'c' or type == 'o': + # Compute Gramian by the Slycot routine sb03md + # make sure Slycot is installed try: from slycot import sb03md except ImportError: raise ControlSlycot("can't find slycot module 'sb03md'") - if type=='c': + if type == 'c': tra = 'T' - C = -np.dot(sys.B,sys.B.transpose()) - elif type=='o': + C = -np.dot(sys.B, sys.B.transpose()) + elif type == 'o': tra = 'N' - C = -np.dot(sys.C.transpose(),sys.C) + C = -np.dot(sys.C.transpose(), sys.C) n = sys.states - U = np.zeros((n,n)) + U = np.zeros((n, n)) A = np.array(sys.A) # convert to NumPy array for slycot - X,scale,sep,ferr,w = sb03md(n, C, A, U, dico, job='X', fact='N', trana=tra) + X, scale, sep, ferr, w = sb03md( + n, C, A, U, dico, job='X', fact='N', trana=tra) gram = X return _ssmatrix(gram) - elif type=='cf' or type=='of': - #Compute cholesky factored gramian from slycot routine sb03od + elif type == 'cf' or type == 'of': + # Compute cholesky factored gramian from slycot routine sb03od try: from slycot import sb03od except ImportError: raise ControlSlycot("can't find slycot module 'sb03od'") - tra='N' + tra = 'N' n = sys.states - Q = np.zeros((n,n)) + Q = np.zeros((n, n)) A = np.array(sys.A) # convert to NumPy array for slycot - if type=='cf': + if type == 'cf': m = sys.B.shape[1] B = np.zeros_like(A) - B[0:m,0:n] = sys.B.transpose() - X,scale,w = sb03od(n, m, A.transpose(), Q, B, dico, fact='N', trans=tra) - elif type=='of': + B[0:m, 0:n] = sys.B.transpose() + X, scale, w = sb03od( + n, m, A.transpose(), Q, B, dico, fact='N', trans=tra) + elif type == 'of': m = sys.C.shape[0] C = np.zeros_like(A) - C[0:n,0:m] = sys.C.transpose() - X,scale,w = sb03od(n, m, A, Q, C.transpose(), dico, fact='N', trans=tra) + C[0:n, 0:m] = sys.C.transpose() + X, scale, w = sb03od( + n, m, A, Q, C.transpose(), dico, fact='N', trans=tra) gram = X return _ssmatrix(gram)