From 37b06f441c4e821234fcf07d8a778f1707e845bc Mon Sep 17 00:00:00 2001 From: Jed Date: Fri, 17 Feb 2017 13:57:46 -0500 Subject: [PATCH 01/21] Added notebook for doing development work on IO-Cleanup. Added IO cleanup function to xferfcn module. --- control/xferfcn.py | 92 ++++++++------- notebooks/IO_Cleanup.ipynb | 231 +++++++++++++++++++++++++++++++++++++ 2 files changed, 280 insertions(+), 43 deletions(-) create mode 100644 notebooks/IO_Cleanup.ipynb diff --git a/control/xferfcn.py b/control/xferfcn.py index 1559fa047..5fd9d9107 100644 --- a/control/xferfcn.py +++ b/control/xferfcn.py @@ -10,6 +10,7 @@ # Python 3 compatibility (needs to go here) from __future__ import print_function from __future__ import division +from __future__ import absolute_import """Copyright (c) 2010 by California Institute of Technology All rights reserved. @@ -55,6 +56,13 @@ from numpy import angle, any, array, empty, finfo, insert, ndarray, ones, \ polyadd, polymul, polyval, roots, sort, sqrt, zeros, squeeze, exp, pi, \ where, delete, real, poly, poly1d + +from numpy import int, int8, int16, int32, int64 +from numpy import float, float16, float32, float64, float128 +from numpy import complex, complex64, complex128, complex256 + +from copy import deepcopy + import numpy as np from scipy.signal import lti, tf2zpk, zpk2tf, cont2discrete from copy import deepcopy @@ -96,7 +104,7 @@ def __init__(self, *args): where sys is a TransferFunction object (continuous or discrete). """ - + args = deepcopy(args) if len(args) == 2: # The user provided a numerator and a denominator. (num, den) = args @@ -120,48 +128,8 @@ def __init__(self, *args): raise ValueError("Needs 1, 2 or 3 arguments; received %i." % len(args)) - # Make num and den into lists of lists of arrays, if necessary. - # Beware: this is a shallow copy! This should be okay, - # but be careful. - data = [num, den] - for i in range(len(data)): - # Check for a scalar (including 0d ndarray) - if (isinstance(data[i], (int, float, complex)) or - (isinstance(data[i], ndarray) and data[i].ndim == 0)): - # Convert scalar to list of list of array. - if (isinstance(data[i], int)): - # Convert integers to floats at this point - data[i] = [[array([data[i]], dtype=float)]] - else: - data[i] = [[array([data[i]])]] - elif (isinstance(data[i], (list, tuple, ndarray)) and - isinstance(data[i][0], (int, float, complex))): - # Convert array to list of list of array. - if (isinstance(data[i][0], int)): - # Convert integers to floats at this point - #! Not sure this covers all cases correctly - data[i] = [[array(data[i], dtype=float)]] - else: - data[i] = [[array(data[i])]] - elif (isinstance(data[i], list) and - isinstance(data[i][0], list) and - isinstance(data[i][0][0], (list, tuple, ndarray)) and - isinstance(data[i][0][0][0], (int, float, complex))): - # We might already have the right format. Convert the - # coefficient vectors to arrays, if necessary. - for j in range(len(data[i])): - for k in range(len(data[i][j])): - if (isinstance(data[i][j][k], int)): - data[i][j][k] = array(data[i][j][k], dtype=float) - else: - data[i][j][k] = array(data[i][j][k]) - else: - # If the user passed in anything else, then it's unclear what - # the meaning is. - raise TypeError("The numerator and denominator inputs must be \ -scalars or vectors (for\nSISO), or lists of lists of vectors (for SISO or \ -MIMO).") - [num, den] = data + num = tf_clean_parts(num) + den = tf_clean_parts(den) inputs = len(num[0]) outputs = len(num) @@ -1349,3 +1317,41 @@ def tfdata(sys): tf = _convertToTransferFunction(sys) return (tf.num, tf.den) + +def tf_clean_parts(data): + ''' + Return a valid, cleaned up numerator or denominator + for the TransferFunction class. + + Parameters: + data: numerator or denominator of a transfer function. + + Returns: + data: correctly formatted transfer function part. + + ''' + valid_types = (int, int8, int16, int32, int64, + float, float16, float32, float64, float128, + complex, complex64, complex128, complex256) + valid_collection = (list, tuple, ndarray) + + if (isinstance(data, valid_types) or + (isinstance(data, ndarray) and data.ndim == 0)): + return [[array([data], dtype=float)]] + elif (isinstance(data, valid_collection) and + all([isinstance(d, valid_types) for d in data])): + return [[array(data, dtype=float)]] + elif (isinstance(data, list) and + isinstance(data[0], list) and + (isinstance(data[0][0], valid_collection) and + isinstance(data[0][0][0], valid_types))): + for j in range(len(data)): + for k in range(len(data[j])): + data[j][k] = array(data[j][k], dtype=float) + return data + else: + # If the user passed in anything else, then it's unclear what + # the meaning is. + raise TypeError("The numerator and denominator inputs must be \ +scalars or vectors (for\nSISO), or lists of lists of vectors (for SISO or \ +MIMO).") \ No newline at end of file diff --git a/notebooks/IO_Cleanup.ipynb b/notebooks/IO_Cleanup.ipynb new file mode 100644 index 000000000..4aad9b79c --- /dev/null +++ b/notebooks/IO_Cleanup.ipynb @@ -0,0 +1,231 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "from numpy import int, int8, int16, int32, int64\n", + "from numpy import float, float16, float32, float64, float128\n", + "from numpy import complex, complex64, complex128, complex256\n", + "from numpy import all, ndarray, array\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "def tf_clean_parts(data):\n", + " '''\n", + " Return a valid, cleaned up numerator or denominator \n", + " for the TransferFunction class.\n", + " \n", + " Parameters:\n", + " data: numerator or denominator of a transfer function.\n", + " \n", + " Returns:\n", + " data: correctly formatted transfer function part.\n", + " \n", + " '''\n", + " valid_types = (int, int8, int16, int32, int64,\n", + " float, float16, float32, float64, float128,\n", + " complex, complex64, complex128, complex256)\n", + " valid_collection = (list, tuple, ndarray)\n", + "\n", + " if (isinstance(data, valid_types) or\n", + " (isinstance(data, ndarray) and data.ndim == 0)):\n", + " return [[array([data], dtype=float)]]\n", + " elif (isinstance(data, valid_collection) and\n", + " all([isinstance(d, valid_types) for d in data])):\n", + " return [[array(data, dtype=float)]]\n", + " elif (isinstance(data, list) and\n", + " isinstance(data[0], list) and\n", + " (isinstance(data[0][0], valid_collection) and \n", + " isinstance(data[0][0][0], valid_types))):\n", + " for j in range(len(data)):\n", + " for k in range(len(data[j])):\n", + " data[j][k] = array(data[j][k], dtype=float)\n", + " return data\n", + " else:\n", + " # If the user passed in anything else, then it's unclear what\n", + " # the meaning is.\n", + " raise TypeError(\"The numerator and denominator inputs must be \\\n", + "scalars or vectors (for\\nSISO), or lists of lists of vectors (for SISO or \\\n", + "MIMO).\")" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "num = 1\n", + "num_ = tf_clean_parts(num)\n", + "np.testing.assert_array_equal(num_[0][0], array([1.0], dtype=float))" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "num = [1]\n", + "num_ = tf_clean_parts(num)\n", + "np.testing.assert_array_equal(num_[0][0], array([1.0], dtype=float))" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false, + "scrolled": true + }, + "outputs": [], + "source": [ + "num = [1,1]\n", + "num_ = tf_clean_parts(num)\n", + "np.testing.assert_array_equal(num_[0][0], array([1.0, 1.0], dtype=float))" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "num = [[[1,1],[2,2]]]\n", + "num_ = tf_clean_parts(num)\n", + "np.testing.assert_array_equal(num_[0][0], array([1.0, 1.0], dtype=float))\n", + "np.testing.assert_array_equal(num_[0][1], array([2.0, 2.0], dtype=float))" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "num = [[[1.0,1.0],[2.0,2.0]]]\n", + "num_ = tf_clean_parts(num)\n", + "np.testing.assert_array_equal(num_[0][0], array([1.0, 1.0], dtype=float))\n", + "np.testing.assert_array_equal(num_[0][1], array([2.0, 2.0], dtype=float))" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "num = [[array([1,1]),array([2,2])]]\n", + "num_ = tf_clean_parts(num)\n", + "np.testing.assert_array_equal(num_[0][0], array([1.0, 1.0], dtype=float))\n", + "np.testing.assert_array_equal(num_[0][1], array([2.0, 2.0], dtype=float))" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "num = [[array([1.0,1.0]),array([2.0,2.0])]]\n", + "num_ = tf_clean_parts(num)\n", + "np.testing.assert_array_equal(num_[0][0], array([1.0, 1.0], dtype=float))\n", + "np.testing.assert_array_equal(num_[0][1], array([2.0, 2.0], dtype=float))" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[[array([ 1., 2., 3.]), array([ 1., 2., 3.])]]" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "num = [[[1,2,3],[1,2,3]]]\n", + "tf_clean_parts(num)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "ename": "TypeError", + "evalue": "The numerator and denominator inputs must be scalars or vectors (for\nSISO), or lists of lists of vectors (for SISO or MIMO).", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0mnum\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mtf_clean_parts\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnum\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[0;32m\u001b[0m in \u001b[0;36mtf_clean_parts\u001b[0;34m(data)\u001b[0m\n\u001b[1;32m 33\u001b[0m \u001b[0;31m# If the user passed in anything else, then it's unclear what\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 34\u001b[0m \u001b[0;31m# the meaning is.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 35\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"The numerator and denominator inputs must be scalars or vectors (for\\nSISO), or lists of lists of vectors (for SISO or MIMO).\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[0;31mTypeError\u001b[0m: The numerator and denominator inputs must be scalars or vectors (for\nSISO), or lists of lists of vectors (for SISO or MIMO)." + ] + } + ], + "source": [ + "num = [[1,2,3],[1,2,3]]\n", + "tf_clean_parts(num)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} From feb68f5266e046b1a425b65c19780157c42454d7 Mon Sep 17 00:00:00 2001 From: Jed Date: Fri, 17 Feb 2017 13:58:46 -0500 Subject: [PATCH 02/21] Added .eggs directory to gitignore. From .eggs/README.txt: "This directory contains eggs that were downloaded by setuptools to build, test, and run plug-ins. This directory caches those eggs to prevent repeated downloads. However, it is safe to delete this directory." --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index c63a6cf06..ff6a58c0a 100644 --- a/.gitignore +++ b/.gitignore @@ -18,3 +18,4 @@ examples/.ipynb_checkpoints/ .project Untitled*.ipynb *.idea/ +.eggs From 5e2d95d116454338f602d06a955af94abb1029fc Mon Sep 17 00:00:00 2001 From: Jed Date: Sat, 18 Feb 2017 11:23:47 -0500 Subject: [PATCH 03/21] Broke out the numerator/denominator cleanup function into a separate function. Added unit tests to test the cleanup function. Added deep copy to TransferFunction args. --- control/tests/xferfcn_input_test.py | 261 ++++++++++++++++++++++++++++ control/xferfcn.py | 19 +- 2 files changed, 271 insertions(+), 9 deletions(-) create mode 100644 control/tests/xferfcn_input_test.py diff --git a/control/tests/xferfcn_input_test.py b/control/tests/xferfcn_input_test.py new file mode 100644 index 000000000..d66df29a1 --- /dev/null +++ b/control/tests/xferfcn_input_test.py @@ -0,0 +1,261 @@ +#!/usr/bin/env python +# +# xferfcn_test.py - test TransferFunction class +# RMM, 30 Mar 2011 (based on TestXferFcn from v0.4a) + +import unittest +import numpy as np + +from numpy import int, int8, int16, int32, int64 +from numpy import float, float16, float32, float64, float128 +from numpy import all, ndarray, array + +from control.xferfcn import cleanPart + +class TestXferFcnInput(unittest.TestCase): + """These are tests for functionality of cleaning and validating + XferFucnInput.""" + + # Tests for raising exceptions. + def testBadInputType(self): + """Give the part cleaner invalid input type.""" + + self.assertRaises(TypeError, cleanPart, [[0., 1.], [2., 3.]]) + + def testBadInputType2(self): + """Give the part cleaner another invalid input type.""" + self.assertRaises(TypeError, cleanPart, [1,"a"]) + + def testScalar(self): + """Test single scalar value.""" + num = 1 + num_ = cleanPart(num) + + assert isinstance(num_, list) + assert np.all([isinstance(part, list) for part in num_]) + np.testing.assert_array_equal(num_[0][0], array([1.0], dtype=float)) + + def testListScalar(self): + """Test single scalar value in list.""" + num = [1] + num_ = cleanPart(num) + + assert isinstance(num_, list) + assert np.all([isinstance(part, list) for part in num_]) + np.testing.assert_array_equal(num_[0][0], array([1.0], dtype=float)) + + def testTupleScalar(self): + """Test single scalar value in tuple.""" + num = (1) + num_ = cleanPart(num) + + assert isinstance(num_, list) + assert np.all([isinstance(part, list) for part in num_]) + np.testing.assert_array_equal(num_[0][0], array([1.0], dtype=float)) + + def testList(self): + """Test multiple values in a list.""" + num = [1, 2] + num_ = cleanPart(num) + + assert isinstance(num_, list) + assert np.all([isinstance(part, list) for part in num_]) + np.testing.assert_array_equal(num_[0][0], array([1.0, 2.0], dtype=float)) + + def testTuple(self): + """Test multiple values in tuple.""" + num = (1, 2) + num_ = cleanPart(num) + + assert isinstance(num_, list) + assert np.all([isinstance(part, list) for part in num_]) + np.testing.assert_array_equal(num_[0][0], array([1.0, 2.0], dtype=float)) + + def testAllScalarTypes(self): + """Test single scalar value for all valid data types.""" + for dtype in [int, int8, int16, int32, int64, float, float16, float32, float64, float128]: + num = dtype(1) + num_ = cleanPart(num) + + assert isinstance(num_, list) + assert np.all([isinstance(part, list) for part in num_]) + np.testing.assert_array_equal(num_[0][0], array([1.0], dtype=float)) + + def testNpArray(self): + """Test multiple values in numpy array.""" + num = np.array([1, 2]) + num_ = cleanPart(num) + + assert isinstance(num_, list) + assert np.all([isinstance(part, list) for part in num_]) + np.testing.assert_array_equal(num_[0][0], array([1.0, 2.0], dtype=float)) + + def testAllNumpyArrayTypes(self): + """Test scalar value in numpy array of ndim=0 for all data types.""" + for dtype in [int, int8, int16, int32, int64, float, float16, float32, float64, float128]: + num = np.array(1, dtype=dtype) + num_ = cleanPart(num) + + assert isinstance(num_, list) + assert np.all([isinstance(part, list) for part in num_]) + np.testing.assert_array_equal(num_[0][0], array([1.0], dtype=float)) + + def testAllNumpyArrayTypes2(self): + """Test numpy array for all types.""" + for dtype in [int, int8, int16, int32, int64, float, float16, float32, float64, float128]: + num = np.array([1, 2], dtype=dtype) + num_ = cleanPart(num) + + assert isinstance(num_, list) + assert np.all([isinstance(part, list) for part in num_]) + np.testing.assert_array_equal(num_[0][0], array([1.0, 2.0], dtype=float)) + + def testListAllTypes(self): + """Test list of a single value for all data types.""" + for dtype in [int, int8, int16, int32, int64, float, float16, float32, float64, float128]: + num = [dtype(1)] + num_ = cleanPart(num) + assert isinstance(num_, list) + assert np.all([isinstance(part, list) for part in num_]) + np.testing.assert_array_equal(num_[0][0], array([1.0], dtype=float)) + + def testListAllTypes2(self): + """List of list of numbers of all data types.""" + for dtype in [int, int8, int16, int32, int64, float, float16, float32, float64, float128]: + num = [dtype(1), dtype(2)] + num_ = cleanPart(num) + assert isinstance(num_, list) + assert np.all([isinstance(part, list) for part in num_]) + np.testing.assert_array_equal(num_[0][0], array([1.0, 2.0], dtype=float)) + + def testTupleAllTypes(self): + """Test tuple of a single value for all data types.""" + for dtype in [int, int8, int16, int32, int64, float, float16, float32, float64, float128]: + num = (dtype(1),) + num_ = cleanPart(num) + assert isinstance(num_, list) + assert np.all([isinstance(part, list) for part in num_]) + np.testing.assert_array_equal(num_[0][0], array([1.0], dtype=float)) + + def testTupleAllTypes2(self): + """Test tuple of a single value for all data types.""" + for dtype in [int, int8, int16, int32, int64, float, float16, float32, float64, float128]: + num = (dtype(1), dtype(2)) + num_ = cleanPart(num) + assert isinstance(num_, list) + assert np.all([isinstance(part, list) for part in num_]) + np.testing.assert_array_equal(num_[0][0], array([1, 2], dtype=float)) + + def testListListListInt(self): + """ Test an int in a list of a list of a list.""" + num = [[[1]]] + num_ = cleanPart(num) + assert isinstance(num_, list) + assert np.all([isinstance(part, list) for part in num_]) + np.testing.assert_array_equal(num_[0][0], array([1.0], dtype=float)) + + def testListListListFloat(self): + """ Test a float in a list of a list of a list.""" + num = [[[1.0]]] + num_ = cleanPart(num) + assert isinstance(num_, list) + assert np.all([isinstance(part, list) for part in num_]) + np.testing.assert_array_equal(num_[0][0], array([1.0], dtype=float)) + + def testListListListInts(self): + """Test 2 lists of ints in a list in a list.""" + num = [[[1,1],[2,2]]] + num_ = cleanPart(num) + + assert isinstance(num_, list) + assert np.all([isinstance(part, list) for part in num_]) + np.testing.assert_array_equal(num_[0][0], array([1.0, 1.0], dtype=float)) + np.testing.assert_array_equal(num_[0][1], array([2.0, 2.0], dtype=float)) + + def testListListListFloats(self): + """Test 2 lists of ints in a list in a list.""" + num = [[[1.0,1.0],[2.0,2.0]]] + num_ = cleanPart(num) + + assert isinstance(num_, list) + assert np.all([isinstance(part, list) for part in num_]) + np.testing.assert_array_equal(num_[0][0], array([1.0, 1.0], dtype=float)) + np.testing.assert_array_equal(num_[0][1], array([2.0, 2.0], dtype=float)) + + def testListListArray(self): + """List of list of numpy arrays for all valid types.""" + for dtype in int, int8, int16, int32, int64, float, float16, float32, float64, float128: + num = [[array([1,1], dtype=dtype),array([2,2], dtype=dtype)]] + num_ = cleanPart(num) + + assert isinstance(num_, list) + assert np.all([isinstance(part, list) for part in num_]) + np.testing.assert_array_equal(num_[0][0], array([1.0, 1.0], dtype=float)) + np.testing.assert_array_equal(num_[0][1], array([2.0, 2.0], dtype=float)) + + def testTupleListArray(self): + """Tuple of list of numpy arrays for all valid types.""" + for dtype in int, int8, int16, int32, int64, float, float16, float32, float64, float128: + num = ([array([1,1], dtype=dtype),array([2,2], dtype=dtype)],) + num_ = cleanPart(num) + + assert isinstance(num_, list) + assert np.all([isinstance(part, list) for part in num_]) + np.testing.assert_array_equal(num_[0][0], array([1.0, 1.0], dtype=float)) + np.testing.assert_array_equal(num_[0][1], array([2.0, 2.0], dtype=float)) + + def testListTupleArray(self): + """List of tuple of numpy array for all valid types.""" + for dtype in int, int8, int16, int32, int64, float, float16, float32, float64, float128: + num = [(array([1,1], dtype=dtype),array([2,2], dtype=dtype))] + num_ = cleanPart(num) + + assert isinstance(num_, list) + assert np.all([isinstance(part, list) for part in num_]) + np.testing.assert_array_equal(num_[0][0], array([1.0, 1.0], dtype=float)) + np.testing.assert_array_equal(num_[0][1], array([2.0, 2.0], dtype=float)) + + def testTupleTuplesArrays(self): + """Tuple of tuples of numpy arrays for all valid types.""" + for dtype in int, int8, int16, int32, int64, float, float16, float32, float64, float128: + num = ((array([1,1], dtype=dtype),array([2,2], dtype=dtype)), + (array([3,4], dtype=dtype),array([4,4], dtype=dtype))) + num_ = cleanPart(num) + + assert isinstance(num_, list) + assert np.all([isinstance(part, list) for part in num_]) + np.testing.assert_array_equal(num_[0][0], array([1.0, 1.0], dtype=float)) + np.testing.assert_array_equal(num_[0][1], array([2.0, 2.0], dtype=float)) + + def testListTuplesArrays(self): + """List of tuples of numpy arrays for all valid types.""" + for dtype in int, int8, int16, int32, int64, float, float16, float32, float64, float128: + num = [(array([1,1], dtype=dtype),array([2,2], dtype=dtype)), + (array([3,4], dtype=dtype),array([4,4], dtype=dtype))] + num_ = cleanPart(num) + + assert isinstance(num_, list) + assert np.all([isinstance(part, list) for part in num_]) + np.testing.assert_array_equal(num_[0][0], array([1.0, 1.0], dtype=float)) + np.testing.assert_array_equal(num_[0][1], array([2.0, 2.0], dtype=float)) + + def testListListArrays(self): + """List of list of numpy arrays for all valid types.""" + for dtype in int, int8, int16, int32, int64, float, float16, float32, float64, float128: + num = [[array([1,1], dtype=dtype),array([2,2], dtype=dtype)], + [array([3,3], dtype=dtype),array([4,4], dtype=dtype)]] + num_ = cleanPart(num) + + assert len(num_) == 2 + assert np.all([isinstance(part, list) for part in num_]) + assert np.all([len(part) == 2 for part in num_]) + np.testing.assert_array_equal(num_[0][0], array([1.0, 1.0], dtype=float)) + np.testing.assert_array_equal(num_[0][1], array([2.0, 2.0], dtype=float)) + np.testing.assert_array_equal(num_[1][0], array([3.0, 3.0], dtype=float)) + np.testing.assert_array_equal(num_[1][1], array([4.0, 4.0], dtype=float)) + +def suite(): + return unittest.TestLoader().loadTestsFromTestCase(TestXferFcnInput) + +if __name__ == "__main__": + unittest.main() \ No newline at end of file diff --git a/control/xferfcn.py b/control/xferfcn.py index 5fd9d9107..f5c3475a7 100644 --- a/control/xferfcn.py +++ b/control/xferfcn.py @@ -128,8 +128,8 @@ def __init__(self, *args): raise ValueError("Needs 1, 2 or 3 arguments; received %i." % len(args)) - num = tf_clean_parts(num) - den = tf_clean_parts(den) + num = cleanPart(num) + den = cleanPart(den) inputs = len(num[0]) outputs = len(num) @@ -1318,7 +1318,7 @@ def tfdata(sys): return (tf.num, tf.den) -def tf_clean_parts(data): +def cleanPart(data): ''' Return a valid, cleaned up numerator or denominator for the TransferFunction class. @@ -1328,11 +1328,10 @@ def tf_clean_parts(data): Returns: data: correctly formatted transfer function part. - + ; ''' valid_types = (int, int8, int16, int32, int64, - float, float16, float32, float64, float128, - complex, complex64, complex128, complex256) + float, float16, float32, float64, float128) valid_collection = (list, tuple, ndarray) if (isinstance(data, valid_types) or @@ -1341,11 +1340,13 @@ def tf_clean_parts(data): elif (isinstance(data, valid_collection) and all([isinstance(d, valid_types) for d in data])): return [[array(data, dtype=float)]] - elif (isinstance(data, list) and - isinstance(data[0], list) and + elif (isinstance(data, (list, tuple)) and + isinstance(data[0], (list, tuple)) and (isinstance(data[0][0], valid_collection) and - isinstance(data[0][0][0], valid_types))): + all([isinstance(d, valid_types) for d in data[0][0]]))): + data = list(data) for j in range(len(data)): + data[j] = list(data[j]) for k in range(len(data[j])): data[j][k] = array(data[j][k], dtype=float) return data From 509ff8e849f27539016b15e41df6e8e61acd6a29 Mon Sep 17 00:00:00 2001 From: Jed Date: Sat, 18 Feb 2017 11:28:44 -0500 Subject: [PATCH 04/21] Removed test notebook --- notebooks/IO_Cleanup.ipynb | 231 ------------------------------------- 1 file changed, 231 deletions(-) delete mode 100644 notebooks/IO_Cleanup.ipynb diff --git a/notebooks/IO_Cleanup.ipynb b/notebooks/IO_Cleanup.ipynb deleted file mode 100644 index 4aad9b79c..000000000 --- a/notebooks/IO_Cleanup.ipynb +++ /dev/null @@ -1,231 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [ - "from numpy import int, int8, int16, int32, int64\n", - "from numpy import float, float16, float32, float64, float128\n", - "from numpy import complex, complex64, complex128, complex256\n", - "from numpy import all, ndarray, array\n", - "import numpy as np" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "def tf_clean_parts(data):\n", - " '''\n", - " Return a valid, cleaned up numerator or denominator \n", - " for the TransferFunction class.\n", - " \n", - " Parameters:\n", - " data: numerator or denominator of a transfer function.\n", - " \n", - " Returns:\n", - " data: correctly formatted transfer function part.\n", - " \n", - " '''\n", - " valid_types = (int, int8, int16, int32, int64,\n", - " float, float16, float32, float64, float128,\n", - " complex, complex64, complex128, complex256)\n", - " valid_collection = (list, tuple, ndarray)\n", - "\n", - " if (isinstance(data, valid_types) or\n", - " (isinstance(data, ndarray) and data.ndim == 0)):\n", - " return [[array([data], dtype=float)]]\n", - " elif (isinstance(data, valid_collection) and\n", - " all([isinstance(d, valid_types) for d in data])):\n", - " return [[array(data, dtype=float)]]\n", - " elif (isinstance(data, list) and\n", - " isinstance(data[0], list) and\n", - " (isinstance(data[0][0], valid_collection) and \n", - " isinstance(data[0][0][0], valid_types))):\n", - " for j in range(len(data)):\n", - " for k in range(len(data[j])):\n", - " data[j][k] = array(data[j][k], dtype=float)\n", - " return data\n", - " else:\n", - " # If the user passed in anything else, then it's unclear what\n", - " # the meaning is.\n", - " raise TypeError(\"The numerator and denominator inputs must be \\\n", - "scalars or vectors (for\\nSISO), or lists of lists of vectors (for SISO or \\\n", - "MIMO).\")" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "num = 1\n", - "num_ = tf_clean_parts(num)\n", - "np.testing.assert_array_equal(num_[0][0], array([1.0], dtype=float))" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "num = [1]\n", - "num_ = tf_clean_parts(num)\n", - "np.testing.assert_array_equal(num_[0][0], array([1.0], dtype=float))" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": { - "collapsed": false, - "scrolled": true - }, - "outputs": [], - "source": [ - "num = [1,1]\n", - "num_ = tf_clean_parts(num)\n", - "np.testing.assert_array_equal(num_[0][0], array([1.0, 1.0], dtype=float))" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "num = [[[1,1],[2,2]]]\n", - "num_ = tf_clean_parts(num)\n", - "np.testing.assert_array_equal(num_[0][0], array([1.0, 1.0], dtype=float))\n", - "np.testing.assert_array_equal(num_[0][1], array([2.0, 2.0], dtype=float))" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [ - "num = [[[1.0,1.0],[2.0,2.0]]]\n", - "num_ = tf_clean_parts(num)\n", - "np.testing.assert_array_equal(num_[0][0], array([1.0, 1.0], dtype=float))\n", - "np.testing.assert_array_equal(num_[0][1], array([2.0, 2.0], dtype=float))" - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "num = [[array([1,1]),array([2,2])]]\n", - "num_ = tf_clean_parts(num)\n", - "np.testing.assert_array_equal(num_[0][0], array([1.0, 1.0], dtype=float))\n", - "np.testing.assert_array_equal(num_[0][1], array([2.0, 2.0], dtype=float))" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [ - "num = [[array([1.0,1.0]),array([2.0,2.0])]]\n", - "num_ = tf_clean_parts(num)\n", - "np.testing.assert_array_equal(num_[0][0], array([1.0, 1.0], dtype=float))\n", - "np.testing.assert_array_equal(num_[0][1], array([2.0, 2.0], dtype=float))" - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "[[array([ 1., 2., 3.]), array([ 1., 2., 3.])]]" - ] - }, - "execution_count": 10, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "num = [[[1,2,3],[1,2,3]]]\n", - "tf_clean_parts(num)" - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "ename": "TypeError", - "evalue": "The numerator and denominator inputs must be scalars or vectors (for\nSISO), or lists of lists of vectors (for SISO or MIMO).", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", - "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0mnum\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mtf_clean_parts\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnum\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", - "\u001b[0;32m\u001b[0m in \u001b[0;36mtf_clean_parts\u001b[0;34m(data)\u001b[0m\n\u001b[1;32m 33\u001b[0m \u001b[0;31m# If the user passed in anything else, then it's unclear what\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 34\u001b[0m \u001b[0;31m# the meaning is.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 35\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"The numerator and denominator inputs must be scalars or vectors (for\\nSISO), or lists of lists of vectors (for SISO or MIMO).\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", - "\u001b[0;31mTypeError\u001b[0m: The numerator and denominator inputs must be scalars or vectors (for\nSISO), or lists of lists of vectors (for SISO or MIMO)." - ] - } - ], - "source": [ - "num = [[1,2,3],[1,2,3]]\n", - "tf_clean_parts(num)" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.5.2" - } - }, - "nbformat": 4, - "nbformat_minor": 0 -} From d44489f1b0636a1010ad50f4a13a8723254f7b7b Mon Sep 17 00:00:00 2001 From: Jed Date: Sat, 18 Feb 2017 14:11:00 -0500 Subject: [PATCH 05/21] Renamed cleanPart to _cleanPart since it is an internal function, as suggested by @murrayrm --- control/tests/xferfcn_input_test.py | 54 ++++++++++++++--------------- control/xferfcn.py | 8 ++--- 2 files changed, 31 insertions(+), 31 deletions(-) diff --git a/control/tests/xferfcn_input_test.py b/control/tests/xferfcn_input_test.py index d66df29a1..cb9da4427 100644 --- a/control/tests/xferfcn_input_test.py +++ b/control/tests/xferfcn_input_test.py @@ -10,7 +10,7 @@ from numpy import float, float16, float32, float64, float128 from numpy import all, ndarray, array -from control.xferfcn import cleanPart +from control.xferfcn import _cleanPart class TestXferFcnInput(unittest.TestCase): """These are tests for functionality of cleaning and validating @@ -20,16 +20,16 @@ class TestXferFcnInput(unittest.TestCase): def testBadInputType(self): """Give the part cleaner invalid input type.""" - self.assertRaises(TypeError, cleanPart, [[0., 1.], [2., 3.]]) + self.assertRaises(TypeError, _cleanPart, [[0., 1.], [2., 3.]]) def testBadInputType2(self): """Give the part cleaner another invalid input type.""" - self.assertRaises(TypeError, cleanPart, [1,"a"]) + self.assertRaises(TypeError, _cleanPart, [1,"a"]) def testScalar(self): """Test single scalar value.""" num = 1 - num_ = cleanPart(num) + num_ = _cleanPart(num) assert isinstance(num_, list) assert np.all([isinstance(part, list) for part in num_]) @@ -38,7 +38,7 @@ def testScalar(self): def testListScalar(self): """Test single scalar value in list.""" num = [1] - num_ = cleanPart(num) + num_ = _cleanPart(num) assert isinstance(num_, list) assert np.all([isinstance(part, list) for part in num_]) @@ -47,7 +47,7 @@ def testListScalar(self): def testTupleScalar(self): """Test single scalar value in tuple.""" num = (1) - num_ = cleanPart(num) + num_ = _cleanPart(num) assert isinstance(num_, list) assert np.all([isinstance(part, list) for part in num_]) @@ -56,7 +56,7 @@ def testTupleScalar(self): def testList(self): """Test multiple values in a list.""" num = [1, 2] - num_ = cleanPart(num) + num_ = _cleanPart(num) assert isinstance(num_, list) assert np.all([isinstance(part, list) for part in num_]) @@ -65,7 +65,7 @@ def testList(self): def testTuple(self): """Test multiple values in tuple.""" num = (1, 2) - num_ = cleanPart(num) + num_ = _cleanPart(num) assert isinstance(num_, list) assert np.all([isinstance(part, list) for part in num_]) @@ -75,7 +75,7 @@ def testAllScalarTypes(self): """Test single scalar value for all valid data types.""" for dtype in [int, int8, int16, int32, int64, float, float16, float32, float64, float128]: num = dtype(1) - num_ = cleanPart(num) + num_ = _cleanPart(num) assert isinstance(num_, list) assert np.all([isinstance(part, list) for part in num_]) @@ -84,7 +84,7 @@ def testAllScalarTypes(self): def testNpArray(self): """Test multiple values in numpy array.""" num = np.array([1, 2]) - num_ = cleanPart(num) + num_ = _cleanPart(num) assert isinstance(num_, list) assert np.all([isinstance(part, list) for part in num_]) @@ -94,7 +94,7 @@ def testAllNumpyArrayTypes(self): """Test scalar value in numpy array of ndim=0 for all data types.""" for dtype in [int, int8, int16, int32, int64, float, float16, float32, float64, float128]: num = np.array(1, dtype=dtype) - num_ = cleanPart(num) + num_ = _cleanPart(num) assert isinstance(num_, list) assert np.all([isinstance(part, list) for part in num_]) @@ -104,7 +104,7 @@ def testAllNumpyArrayTypes2(self): """Test numpy array for all types.""" for dtype in [int, int8, int16, int32, int64, float, float16, float32, float64, float128]: num = np.array([1, 2], dtype=dtype) - num_ = cleanPart(num) + num_ = _cleanPart(num) assert isinstance(num_, list) assert np.all([isinstance(part, list) for part in num_]) @@ -114,7 +114,7 @@ def testListAllTypes(self): """Test list of a single value for all data types.""" for dtype in [int, int8, int16, int32, int64, float, float16, float32, float64, float128]: num = [dtype(1)] - num_ = cleanPart(num) + num_ = _cleanPart(num) assert isinstance(num_, list) assert np.all([isinstance(part, list) for part in num_]) np.testing.assert_array_equal(num_[0][0], array([1.0], dtype=float)) @@ -123,7 +123,7 @@ def testListAllTypes2(self): """List of list of numbers of all data types.""" for dtype in [int, int8, int16, int32, int64, float, float16, float32, float64, float128]: num = [dtype(1), dtype(2)] - num_ = cleanPart(num) + num_ = _cleanPart(num) assert isinstance(num_, list) assert np.all([isinstance(part, list) for part in num_]) np.testing.assert_array_equal(num_[0][0], array([1.0, 2.0], dtype=float)) @@ -132,7 +132,7 @@ def testTupleAllTypes(self): """Test tuple of a single value for all data types.""" for dtype in [int, int8, int16, int32, int64, float, float16, float32, float64, float128]: num = (dtype(1),) - num_ = cleanPart(num) + num_ = _cleanPart(num) assert isinstance(num_, list) assert np.all([isinstance(part, list) for part in num_]) np.testing.assert_array_equal(num_[0][0], array([1.0], dtype=float)) @@ -141,7 +141,7 @@ def testTupleAllTypes2(self): """Test tuple of a single value for all data types.""" for dtype in [int, int8, int16, int32, int64, float, float16, float32, float64, float128]: num = (dtype(1), dtype(2)) - num_ = cleanPart(num) + num_ = _cleanPart(num) assert isinstance(num_, list) assert np.all([isinstance(part, list) for part in num_]) np.testing.assert_array_equal(num_[0][0], array([1, 2], dtype=float)) @@ -149,7 +149,7 @@ def testTupleAllTypes2(self): def testListListListInt(self): """ Test an int in a list of a list of a list.""" num = [[[1]]] - num_ = cleanPart(num) + num_ = _cleanPart(num) assert isinstance(num_, list) assert np.all([isinstance(part, list) for part in num_]) np.testing.assert_array_equal(num_[0][0], array([1.0], dtype=float)) @@ -157,7 +157,7 @@ def testListListListInt(self): def testListListListFloat(self): """ Test a float in a list of a list of a list.""" num = [[[1.0]]] - num_ = cleanPart(num) + num_ = _cleanPart(num) assert isinstance(num_, list) assert np.all([isinstance(part, list) for part in num_]) np.testing.assert_array_equal(num_[0][0], array([1.0], dtype=float)) @@ -165,7 +165,7 @@ def testListListListFloat(self): def testListListListInts(self): """Test 2 lists of ints in a list in a list.""" num = [[[1,1],[2,2]]] - num_ = cleanPart(num) + num_ = _cleanPart(num) assert isinstance(num_, list) assert np.all([isinstance(part, list) for part in num_]) @@ -175,7 +175,7 @@ def testListListListInts(self): def testListListListFloats(self): """Test 2 lists of ints in a list in a list.""" num = [[[1.0,1.0],[2.0,2.0]]] - num_ = cleanPart(num) + num_ = _cleanPart(num) assert isinstance(num_, list) assert np.all([isinstance(part, list) for part in num_]) @@ -186,7 +186,7 @@ def testListListArray(self): """List of list of numpy arrays for all valid types.""" for dtype in int, int8, int16, int32, int64, float, float16, float32, float64, float128: num = [[array([1,1], dtype=dtype),array([2,2], dtype=dtype)]] - num_ = cleanPart(num) + num_ = _cleanPart(num) assert isinstance(num_, list) assert np.all([isinstance(part, list) for part in num_]) @@ -197,7 +197,7 @@ def testTupleListArray(self): """Tuple of list of numpy arrays for all valid types.""" for dtype in int, int8, int16, int32, int64, float, float16, float32, float64, float128: num = ([array([1,1], dtype=dtype),array([2,2], dtype=dtype)],) - num_ = cleanPart(num) + num_ = _cleanPart(num) assert isinstance(num_, list) assert np.all([isinstance(part, list) for part in num_]) @@ -208,7 +208,7 @@ def testListTupleArray(self): """List of tuple of numpy array for all valid types.""" for dtype in int, int8, int16, int32, int64, float, float16, float32, float64, float128: num = [(array([1,1], dtype=dtype),array([2,2], dtype=dtype))] - num_ = cleanPart(num) + num_ = _cleanPart(num) assert isinstance(num_, list) assert np.all([isinstance(part, list) for part in num_]) @@ -220,7 +220,7 @@ def testTupleTuplesArrays(self): for dtype in int, int8, int16, int32, int64, float, float16, float32, float64, float128: num = ((array([1,1], dtype=dtype),array([2,2], dtype=dtype)), (array([3,4], dtype=dtype),array([4,4], dtype=dtype))) - num_ = cleanPart(num) + num_ = _cleanPart(num) assert isinstance(num_, list) assert np.all([isinstance(part, list) for part in num_]) @@ -232,7 +232,7 @@ def testListTuplesArrays(self): for dtype in int, int8, int16, int32, int64, float, float16, float32, float64, float128: num = [(array([1,1], dtype=dtype),array([2,2], dtype=dtype)), (array([3,4], dtype=dtype),array([4,4], dtype=dtype))] - num_ = cleanPart(num) + num_ = _cleanPart(num) assert isinstance(num_, list) assert np.all([isinstance(part, list) for part in num_]) @@ -244,7 +244,7 @@ def testListListArrays(self): for dtype in int, int8, int16, int32, int64, float, float16, float32, float64, float128: num = [[array([1,1], dtype=dtype),array([2,2], dtype=dtype)], [array([3,3], dtype=dtype),array([4,4], dtype=dtype)]] - num_ = cleanPart(num) + num_ = _cleanPart(num) assert len(num_) == 2 assert np.all([isinstance(part, list) for part in num_]) @@ -258,4 +258,4 @@ def suite(): return unittest.TestLoader().loadTestsFromTestCase(TestXferFcnInput) if __name__ == "__main__": - unittest.main() \ No newline at end of file + unittest.main() diff --git a/control/xferfcn.py b/control/xferfcn.py index f5c3475a7..df384ad74 100644 --- a/control/xferfcn.py +++ b/control/xferfcn.py @@ -128,8 +128,8 @@ def __init__(self, *args): raise ValueError("Needs 1, 2 or 3 arguments; received %i." % len(args)) - num = cleanPart(num) - den = cleanPart(den) + num = _cleanPart(num) + den = _cleanPart(den) inputs = len(num[0]) outputs = len(num) @@ -1318,7 +1318,7 @@ def tfdata(sys): return (tf.num, tf.den) -def cleanPart(data): +def _cleanPart(data): ''' Return a valid, cleaned up numerator or denominator for the TransferFunction class. @@ -1355,4 +1355,4 @@ def cleanPart(data): # the meaning is. raise TypeError("The numerator and denominator inputs must be \ scalars or vectors (for\nSISO), or lists of lists of vectors (for SISO or \ -MIMO).") \ No newline at end of file +MIMO).") From d8d6b02a267e5e4c658353cacf85f05b29ce0c76 Mon Sep 17 00:00:00 2001 From: rssalessio Date: Sun, 26 Mar 2017 17:59:05 +0200 Subject: [PATCH 06/21] Fix minreal not returning a discrete TF --- control/xferfcn.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/control/xferfcn.py b/control/xferfcn.py index 1559fa047..99c8b91d4 100644 --- a/control/xferfcn.py +++ b/control/xferfcn.py @@ -685,7 +685,7 @@ def minreal(self, tol=None): den[i][j] = np.atleast_1d(real(poly(poles))) # end result - return TransferFunction(num, den) + return TransferFunction(num, den, None if self.dt == None else self.dt) def returnScipySignalLTI(self): """Return a list of a list of scipy.signal.lti objects. From 585c0280a8a30c0eefe87e5e5ada3ca2d3724d25 Mon Sep 17 00:00:00 2001 From: rssalessio Date: Wed, 29 Mar 2017 23:24:57 +0200 Subject: [PATCH 07/21] change if to None. add tests --- control/tests/xferfcn_test.py | 13 +++++++++++++ control/xferfcn.py | 2 +- 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/control/tests/xferfcn_test.py b/control/tests/xferfcn_test.py index a6c9ae585..cb33a6190 100644 --- a/control/tests/xferfcn_test.py +++ b/control/tests/xferfcn_test.py @@ -461,6 +461,7 @@ def testMinreal(self): hr = (s+1)/(s**2+s+1) np.testing.assert_array_almost_equal(hm.num[0][0], hr.num[0][0]) np.testing.assert_array_almost_equal(hm.den[0][0], hr.den[0][0]) + np.testing.assert_equal(hm.dt, hr.dt) def testMinreal2(self): """This one gave a problem, due to poly([]) giving simply 1 @@ -474,12 +475,24 @@ def testMinreal2(self): hr = 6205/(s**2+8*s+1241) np.testing.assert_array_almost_equal(H2b.num[0][0], hr.num[0][0]) np.testing.assert_array_almost_equal(H2b.den[0][0], hr.den[0][0]) + np.testing.assert_equal(H2b.dt, hr.dt) def testMinreal3(self): """Regression test for minreal of tf([1,1],[1,1])""" g = TransferFunction([1,1],[1,1]).minreal() np.testing.assert_array_almost_equal(1.0, g.num[0][0]) np.testing.assert_array_almost_equal(1.0, g.den[0][0]) + np.testing.assert_equal(None, g.dt) + + def testMinreal4(self): + """Check minreal on discrete TFs.""" + T = 0.01 + z = TransferFunction([1, 0], [1], T) + h = (z-1.00000000001)*(z+1.0000000001)/((z**2-1)) + hm = h.minreal() + hr = TransferFunction([1], [1], T) + np.testing.assert_array_almost_equal(hm.num[0][0], hr.num[0][0]) + np.testing.assert_equal(hr.dt, hm.dt) @unittest.skipIf(not slycot_check(), "slycot not installed") def testMIMO(self): diff --git a/control/xferfcn.py b/control/xferfcn.py index 99c8b91d4..c22c27c4e 100644 --- a/control/xferfcn.py +++ b/control/xferfcn.py @@ -685,7 +685,7 @@ def minreal(self, tol=None): den[i][j] = np.atleast_1d(real(poly(poles))) # end result - return TransferFunction(num, den, None if self.dt == None else self.dt) + return TransferFunction(num, den, self.dt) def returnScipySignalLTI(self): """Return a list of a list of scipy.signal.lti objects. From d4fc2d0a6f57a841d2fe830f4adab6e92109fb76 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Fri, 22 Dec 2017 19:56:08 -0800 Subject: [PATCH 08/21] BLD: travis update for new/old scipy + improved slycot --- .travis.yml | 27 ++++++++++++++++++--------- 1 file changed, 18 insertions(+), 9 deletions(-) diff --git a/.travis.yml b/.travis.yml index 33aca5703..a3238b2c5 100644 --- a/.travis.yml +++ b/.travis.yml @@ -10,9 +10,19 @@ cache: python: - "2.7" - - "3.3" - - "3.4" - "3.5" + - "3.6" + +# Test against multiple version of SciPy, with and without slycot +# +# Because there were significant changes in SciPy between v0 and v1, we +# test against both of these using the Travis CI environment capability +# +# We also want to test with and without slycot +env: + - SCIPY=0.19.1 SLYCOT= + - SCIPY=0.19.1 SLYCOT=slycot + - SCIPY=1.0.0 SLYCOT= # install required system libraries before_install: @@ -40,16 +50,15 @@ before_install: # Install packages install: - - conda build --python "$TRAVIS_PYTHON_VERSION" conda-recipe - - conda install control --use-local + # Install the desired version of SciPy first, w/ or w/out slycot + - conda install scipy==$SCIPY $SLYCOT + - conda install matplotlib + # Don't use conda for installation of control library [preserves scipy] + - python setup.py install # command to run tests script: - # Before installing Slycot - - python setup.py test - - # Now, get and use Slycot - - conda install slycot + - 'if [ $SLYCOT != "" ]; then python -c "import slycot"; fi' - coverage run setup.py test after_success: From 11c99efdb36fbe7af744d56522c002fba7b3a412 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Sat, 23 Dec 2017 08:40:49 -0800 Subject: [PATCH 09/21] Fixes for SciPy 1.0 compatibility --- control/bdalg.py | 11 ++++++----- control/frdata.py | 11 ++++++----- control/lti.py | 9 +++++---- control/statesp.py | 14 +++++++------- control/xferfcn.py | 33 +++++++++++++++++---------------- 5 files changed, 41 insertions(+), 37 deletions(-) diff --git a/control/bdalg.py b/control/bdalg.py index 9d5aac654..b92496477 100644 --- a/control/bdalg.py +++ b/control/bdalg.py @@ -54,6 +54,7 @@ """ import scipy as sp +import numpy as np from . import xferfcn as tf from . import statesp as ss from . import frdata as frd @@ -221,18 +222,18 @@ def feedback(sys1, sys2=1, sign=-1): """ # Check for correct input types. - if not isinstance(sys1, (int, float, complex, tf.TransferFunction, - ss.StateSpace, frd.FRD)): + if not isinstance(sys1, (int, float, complex, np.number, + tf.TransferFunction, ss.StateSpace, frd.FRD)): raise TypeError("sys1 must be a TransferFunction, StateSpace " + "or FRD object, or a scalar.") - if not isinstance(sys2, (int, float, complex, tf.TransferFunction, - ss.StateSpace, frd.FRD)): + if not isinstance(sys2, (int, float, complex, np.number, + tf.TransferFunction, ss.StateSpace, frd.FRD)): raise TypeError("sys2 must be a TransferFunction, StateSpace " + "or FRD object, or a scalar.") # If sys1 is a scalar, convert it to the appropriate LTI type so that we can # its feedback member function. - if isinstance(sys1, (int, float, complex)): + if isinstance(sys1, (int, float, complex, np.number)): if isinstance(sys2, tf.TransferFunction): sys1 = tf._convertToTransferFunction(sys1) elif isinstance(sys2, ss.StateSpace): diff --git a/control/frdata.py b/control/frdata.py index ff1663b15..313da12b2 100644 --- a/control/frdata.py +++ b/control/frdata.py @@ -49,6 +49,7 @@ """ # External function declarations +import numpy as np from numpy import angle, array, empty, ones, \ real, imag, matrix, absolute, eye, linalg, where, dot from scipy.interpolate import splprep, splev @@ -219,7 +220,7 @@ def __mul__(self, other): """Multiply two LTI objects (serial connection).""" # Convert the second argument to a transfer function. - if isinstance(other, (int, float, complex)): + if isinstance(other, (int, float, complex, np.number)): return FRD(self.fresp * other, self.omega, smooth=(self.ifunc is not None)) else: @@ -245,7 +246,7 @@ def __rmul__(self, other): """Right Multiply two LTI objects (serial connection).""" # Convert the second argument to an frd function. - if isinstance(other, (int, float, complex)): + if isinstance(other, (int, float, complex, np.number)): return FRD(self.fresp * other, self.omega, smooth=(self.ifunc is not None)) else: @@ -272,7 +273,7 @@ def __rmul__(self, other): def __truediv__(self, other): """Divide two LTI objects.""" - if isinstance(other, (int, float, complex)): + if isinstance(other, (int, float, complex, np.number)): return FRD(self.fresp * (1/other), self.omega, smooth=(self.ifunc is not None)) else: @@ -295,7 +296,7 @@ def __div__(self, other): # TODO: Division of MIMO transfer function objects is not written yet. def __rtruediv__(self, other): """Right divide two LTI objects.""" - if isinstance(other, (int, float, complex)): + if isinstance(other, (int, float, complex, np.number)): return FRD(other / self.fresp, self.omega, smooth=(self.ifunc is not None)) else: @@ -449,7 +450,7 @@ def _convertToFRD(sys, omega, inputs=1, outputs=1): return FRD(fresp, omega, smooth=True) - elif isinstance(sys, (int, float, complex)): + elif isinstance(sys, (int, float, complex, np.number)): fresp = ones((outputs, inputs, len(omega)), dtype=float)*sys return FRD(fresp, omega, smooth=True) diff --git a/control/lti.py b/control/lti.py index ab7672d51..aeef5cbb1 100644 --- a/control/lti.py +++ b/control/lti.py @@ -12,6 +12,7 @@ timebaseEqual() """ +import numpy as np from numpy import absolute, real __all__ = ['issiso', 'timebase', 'timebaseEqual', 'isdtime', 'isctime', @@ -96,7 +97,7 @@ def dcgain(self): # Test to see if a system is SISO def issiso(sys, strict=False): - if isinstance(sys, (int, float, complex)) and not strict: + if isinstance(sys, (int, float, complex, np.number)) and not strict: return True elif not isinstance(sys, LTI): raise ValueError("Object is not an LTI system") @@ -114,7 +115,7 @@ def timebase(sys, strict=True): set to False, dt = True will be returned as 1. """ # System needs to be either a constant or an LTI system - if isinstance(sys, (int, float, complex)): + if isinstance(sys, (int, float, complex, np.number)): return None elif not isinstance(sys, LTI): raise ValueError("Timebase not defined") @@ -162,7 +163,7 @@ def isdtime(sys, strict=False): """ # Check to see if this is a constant - if isinstance(sys, (int, float, complex)): + if isinstance(sys, (int, float, complex, np.number)): # OK as long as strict checking is off return True if not strict else False @@ -187,7 +188,7 @@ def isctime(sys, strict=False): """ # Check to see if this is a constant - if isinstance(sys, (int, float, complex)): + if isinstance(sys, (int, float, complex, np.number)): # OK as long as strict checking is off return True if not strict else False diff --git a/control/statesp.py b/control/statesp.py index 13e662f03..816485f77 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -58,6 +58,7 @@ from numpy.random import rand, randn from numpy.linalg import solve, eigvals, matrix_rank from numpy.linalg.linalg import LinAlgError +import scipy as sp from scipy.signal import lti, cont2discrete # from exceptions import Exception import warnings @@ -218,7 +219,7 @@ def __add__(self, other): """Add two LTI systems (parallel connection).""" # Check for a couple of special cases - if (isinstance(other, (int, float, complex))): + if (isinstance(other, (int, float, complex, np.number))): # Just adding a scalar; put it in the D matrix A, B, C = self.A, self.B, self.C; D = self.D + other; @@ -275,7 +276,7 @@ def __mul__(self, other): """Multiply two LTI objects (serial connection).""" # Check for a couple of special cases - if isinstance(other, (int, float, complex)): + if isinstance(other, (int, float, complex, np.number)): # Just multiplying by a scalar; change the output A, B = self.A, self.B C = self.C * other @@ -316,7 +317,7 @@ def __rmul__(self, other): """Right multiply two LTI objects (serial connection).""" # Check for a couple of special cases - if isinstance(other, (int, float, complex)): + if isinstance(other, (int, float, complex, np.number)): # Just multiplying by a scalar; change the input A, C = self.A, self.C; B = self.B * other; @@ -699,11 +700,10 @@ def _convertToStateSpace(sys, **kw): # TODO: do we want to squeeze first and check dimenations? # I think this will fail if num and den aren't 1-D after # the squeeze - lti_sys = lti(squeeze(sys.num), squeeze(sys.den)) - return StateSpace(lti_sys.A, lti_sys.B, lti_sys.C, lti_sys.D, - sys.dt) + A, B, C, D = sp.signal.tf2ss(squeeze(sys.num), squeeze(sys.den)) + return StateSpace(A, B, C, D, sys.dt) - elif isinstance(sys, (int, float, complex)): + elif isinstance(sys, (int, float, complex, np.number)): if "inputs" in kw: inputs = kw["inputs"] else: diff --git a/control/xferfcn.py b/control/xferfcn.py index be21f8509..0d7d2f417 100644 --- a/control/xferfcn.py +++ b/control/xferfcn.py @@ -52,10 +52,11 @@ """ # External function declarations +import numpy as np from numpy import angle, any, array, empty, finfo, insert, ndarray, ones, \ polyadd, polymul, polyval, roots, sort, sqrt, zeros, squeeze, exp, pi, \ where, delete, real, poly, poly1d -import numpy as np +import scipy as sp from scipy.signal import lti, tf2zpk, zpk2tf, cont2discrete from copy import deepcopy from warnings import warn @@ -126,7 +127,7 @@ def __init__(self, *args): data = [num, den] for i in range(len(data)): # Check for a scalar (including 0d ndarray) - if (isinstance(data[i], (int, float, complex)) or + if (isinstance(data[i], (int, float, complex, np.number)) or (isinstance(data[i], ndarray) and data[i].ndim == 0)): # Convert scalar to list of list of array. if (isinstance(data[i], int)): @@ -135,7 +136,7 @@ def __init__(self, *args): else: data[i] = [[array([data[i]])]] elif (isinstance(data[i], (list, tuple, ndarray)) and - isinstance(data[i][0], (int, float, complex))): + isinstance(data[i][0], (int, float, complex, np.number))): # Convert array to list of list of array. if (isinstance(data[i][0], int)): # Convert integers to floats at this point @@ -146,7 +147,8 @@ def __init__(self, *args): elif (isinstance(data[i], list) and isinstance(data[i][0], list) and isinstance(data[i][0][0], (list, tuple, ndarray)) and - isinstance(data[i][0][0][0], (int, float, complex))): + isinstance(data[i][0][0][0], (int, float, complex, + np.number))): # We might already have the right format. Convert the # coefficient vectors to arrays, if necessary. for j in range(len(data[i])): @@ -363,7 +365,7 @@ def __rsub__(self, other): def __mul__(self, other): """Multiply two LTI objects (serial connection).""" # Convert the second argument to a transfer function. - if isinstance(other, (int, float, complex)): + if isinstance(other, (int, float, complex, np.number)): other = _convertToTransferFunction(other, inputs=self.inputs, outputs=self.inputs) else: @@ -410,7 +412,7 @@ def __rmul__(self, other): """Right multiply two LTI objects (serial connection).""" # Convert the second argument to a transfer function. - if isinstance(other, (int, float, complex)): + if isinstance(other, (int, float, complex, np.number)): other = _convertToTransferFunction(other, inputs=self.inputs, outputs=self.inputs) else: @@ -458,7 +460,7 @@ def __rmul__(self, other): def __truediv__(self, other): """Divide two LTI objects.""" - if isinstance(other, (int, float, complex)): + if isinstance(other, (int, float, complex, np.number)): other = _convertToTransferFunction( other, inputs=self.inputs, outputs=self.inputs) @@ -492,7 +494,7 @@ def __div__(self, other): # TODO: Division of MIMO transfer function objects is not written yet. def __rtruediv__(self, other): """Right divide two LTI objects.""" - if isinstance(other, (int, float, complex)): + if isinstance(other, (int, float, complex, np.number)): other = _convertToTransferFunction( other, inputs=self.inputs, outputs=self.inputs) @@ -1128,22 +1130,21 @@ def _convertToTransferFunction(sys, **kw): # Each transfer function matrix row # has a common denominator. den[i][j] = list(tfout[5][i, :]) - # print(num) - # print(den) + except ImportError: # If slycot is not available, use signal.lti (SISO only) if (sys.inputs != 1 or sys.outputs != 1): raise TypeError("No support for MIMO without slycot") - lti_sys = lti(sys.A, sys.B, sys.C, sys.D) - num = squeeze(lti_sys.num) - den = squeeze(lti_sys.den) - # print(num) - # print(den) + # Do the conversion using sp.signal.ss2tf + # Note that this returns a 2D array for the numerator + num, den = sp.signal.ss2tf(sys.A, sys.B, sys.C, sys.D) + num = squeeze(num) # Convert to 1D array + den = squeeze(den) # Probably not needed return TransferFunction(num, den, sys.dt) - elif isinstance(sys, (int, float, complex)): + elif isinstance(sys, (int, float, complex, np.number)): if "inputs" in kw: inputs = kw["inputs"] else: From b6e378ffa7125051a96e81f767c37efdedd78bb8 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Sun, 24 Dec 2017 12:20:38 -0800 Subject: [PATCH 10/21] BLD: updated Travis CI script to build slycot from source --- .travis.yml | 38 +++++++++++++++++++++++++------------- 1 file changed, 25 insertions(+), 13 deletions(-) diff --git a/.travis.yml b/.travis.yml index a3238b2c5..0164195eb 100644 --- a/.travis.yml +++ b/.travis.yml @@ -9,9 +9,9 @@ cache: - $HOME/.local python: - - "2.7" - - "3.5" - "3.6" + - "3.5" + - "2.7" # Test against multiple version of SciPy, with and without slycot # @@ -20,12 +20,20 @@ python: # # We also want to test with and without slycot env: - - SCIPY=0.19.1 SLYCOT= - - SCIPY=0.19.1 SLYCOT=slycot - - SCIPY=1.0.0 SLYCOT= + - SCIPY=scipy SLYCOT=slycot # default, with slycot + - SCIPY=scipy SLYCOT= # default, w/out slycot + - SCIPY="scipy==0.19.1" SLYCOT= # legacy support, w/out slycot # install required system libraries before_install: + # Install gfortran for testing slycot; use apt-get instead of conda in + # order to include the proper CXXABI dependency (updated in GCC 4.9) + # Also need to include liblapack here, to make sure paths are right + - if [[ "$SLYCOT" != "" ]]; then + sudo apt-get update -qq; + sudo apt-get install gfortran liblapack-dev; + fi + # Install display manager to allow testing of plotting functions - export DISPLAY=:99.0 - sh -e /etc/init.d/xvfb start # use miniconda to install numpy/scipy, to avoid lengthy build from source @@ -39,22 +47,26 @@ before_install: - hash -r - conda config --set always_yes yes --set changeps1 no - conda update -q conda - # conda-build must be installed in the conda root environment - - conda install conda-build - conda config --add channels python-control - conda info -a - conda create -q -n test-environment python="$TRAVIS_PYTHON_VERSION" pip coverage - source activate test-environment - # coveralls not in conda repos + # Make sure to look in the right place for python libraries (for slycot) + - export LIBRARY_PATH="$HOME/miniconda/envs/test-environment/lib" + # coveralls not in conda repos => install via pip instead - pip install coveralls # Install packages install: - # Install the desired version of SciPy first, w/ or w/out slycot - - conda install scipy==$SCIPY $SLYCOT - - conda install matplotlib - # Don't use conda for installation of control library [preserves scipy] - - python setup.py install + # Install packages needed by python-control + - conda install $SCIPY matplotlib + # Build slycot from source + # For python 3, need to provide pointer to python library + #! git clone https://github.com/repagh/Slycot.git slycot; + - if [[ "$SLYCOT" != "" ]]; then + git clone https://github.com/python-control/Slycot.git slycot; + cd slycot; python setup.py install; cd ..; + fi # command to run tests script: From 087046fa283e7d1c5ad0156e129933f58007e540 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Thu, 28 Dec 2017 17:13:46 -0800 Subject: [PATCH 11/21] DOC: fix conda-recipe and make_version for conda v3 --- conda-recipe/meta.yaml | 5 +++++ make_version.py | 35 +++++++++++++++++++++-------------- 2 files changed, 26 insertions(+), 14 deletions(-) diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml index 7a578191f..86952af8f 100644 --- a/conda-recipe/meta.yaml +++ b/conda-recipe/meta.yaml @@ -1,7 +1,12 @@ package: name: control + version: {{ GIT_DESCRIBE_TAG }} + +source: + git_url: ../ build: + number: {{ GIT_DESCRIBE_NUMBER }} script: - cd $RECIPE_DIR/.. - $PYTHON make_version.py diff --git a/make_version.py b/make_version.py index bf03dde5c..356f4d747 100644 --- a/make_version.py +++ b/make_version.py @@ -1,3 +1,24 @@ +# make_version.py - generate version information +# +# Author: Clancy Rowley +# Date: 2 Apr 2015 +# Modified: Richard M. Murray, 28 Dec 2017 +# +# This script is used to create the version information for the python- +# control package. The version information is now generated directly from +# tags in the git repository. Now, *before* running setup.py, one runs +# +# python make_version.py +# +# and this generates a file with the version information. This is copied +# from binstar (https://github.com/Binstar/binstar) and seems to work well. +# +# The original version of this script also created version information for +# conda, but this stopped working when conda v3 was released. Instead, we +# now use jinja templates in conda-recipe to create the conda information. +# The current version information is used in setup.py, control/__init__.py, +# and doc/conf.py (for sphinx). + from subprocess import check_output import os @@ -33,19 +54,5 @@ def main(): fd.write('__version__ = "%s.post%s"\n' % (version, build)) fd.write('__commit__ = "%s"\n' % (commit)) - # Write files for conda version number - SRC_DIR = os.environ.get('SRC_DIR', '.') - conda_version_path = os.path.join(SRC_DIR, '__conda_version__.txt') - print("Writing %s" % conda_version_path) - with open(conda_version_path, 'w') as conda_version: - conda_version.write(version) - - conda_buildnum_path = os.path.join(SRC_DIR, '__conda_buildnum__.txt') - print("Writing %s" % conda_buildnum_path) - - with open(conda_buildnum_path, 'w') as conda_buildnum: - conda_buildnum.write(build) - - if __name__ == '__main__': main() From b19ca118dd99420268ee0a29bec78298c41cc0e1 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Thu, 28 Dec 2017 17:33:42 -0800 Subject: [PATCH 12/21] DOC: update docstrings for system creation (see PR #163) --- control/ctrlutil.py | 5 + control/frdata.py | 15 +- control/margins.py | 10 +- control/robust.py | 11 +- control/statesp.py | 25 +- control/statesp.py+ | 1211 +++++++++++++++++++++++++++++++++++++++++++ control/xferfcn.py | 22 +- 7 files changed, 1266 insertions(+), 33 deletions(-) create mode 100644 control/statesp.py+ diff --git a/control/ctrlutil.py b/control/ctrlutil.py index 273257c74..258e5755a 100644 --- a/control/ctrlutil.py +++ b/control/ctrlutil.py @@ -45,6 +45,11 @@ import numpy as np from numpy import pi +# Hack for sphinx.ext.autodoc: if numpy is a mock import, then numpy.pi +# will be assigned to _Mock() and this generates a type error +if not isinstance(pi, float): + pi = 3.14 + __all__ = ['unwrap', 'issys', 'db2mag', 'mag2db'] # Utility function to unwrap an angle measurement diff --git a/control/frdata.py b/control/frdata.py index ff1663b15..9b4c84080 100644 --- a/control/frdata.py +++ b/control/frdata.py @@ -57,7 +57,9 @@ __all__ = ['FRD', 'frd'] class FRD(LTI): - """A class for models defined by Frequency Response Data (FRD) + """FRD(d, w) + + A class for models defined by frequency response data (FRD) The FRD class is used to represent systems in frequency response data form. @@ -80,7 +82,9 @@ class FRD(LTI): epsw = 1e-8 def __init__(self, *args, **kwargs): - """Construct an FRD object + """FRD(d, w) + + Construct an FRD object The default constructor is FRD(d, w), where w is an iterable of frequency points, and d is the matching frequency data. @@ -469,8 +473,9 @@ def _convertToFRD(sys, omega, inputs=1, outputs=1): sys.__class__) def frd(*args): - """ - Construct a Frequency Response Data model, or convert a system + """frd(d, w) + + Construct a frequency response data model frd models store the (measured) frequency response of a system. @@ -500,6 +505,6 @@ def frd(*args): See Also -------- - ss, tf + FRD, ss, tf """ return FRD(*args) diff --git a/control/margins.py b/control/margins.py index 26ec47325..9accd566a 100644 --- a/control/margins.py +++ b/control/margins.py @@ -48,7 +48,6 @@ Date: 14 July 2011 $Id$ - """ import numpy as np @@ -62,7 +61,7 @@ # helper functions for stability_margins def _polyimsplit(pol): """split a polynomial with (iw) applied into a real and an - imaginary part with w applied""" + imaginary part with w applied""" rpencil = np.zeros_like(pol) ipencil = np.zeros_like(pol) rpencil[-1::-4] = 1. @@ -294,8 +293,7 @@ def dstab(w): # Contributed by Steffen Waldherr #! TODO - need to add test functions def phase_crossover_frequencies(sys): - """ - Compute frequencies and gains at intersections with real axis + """Compute frequencies and gains at intersections with real axis in Nyquist plot. Call as: @@ -360,8 +358,8 @@ def margin(*args): Wcp : float Phase crossover frequency (corresponding to gain margin) (in rad/sec) - Margins are of SISO open-loop. If more than one crossover frequency is - detected, returns the lowest corresponding margin. + Margins are of SISO open-loop. If more than one crossover frequency is + detected, returns the lowest corresponding margin. Examples -------- diff --git a/control/robust.py b/control/robust.py index 1b33e93d0..0b9b98c41 100644 --- a/control/robust.py +++ b/control/robust.py @@ -72,7 +72,6 @@ def h2syn(P,nmeas,ncon): >>> K = h2syn(P,nmeas,ncon) """ - #Check for ss system object, need a utility for this? #TODO: Check for continous or discrete, only continuous supported right now @@ -116,11 +115,11 @@ def hinfsyn(P,nmeas,ncon): CL: closed loop system (State-space sys) gam: infinity norm of closed loop system rcond: 4-vector, reciprocal condition estimates of: - 1: control transformation matrix - 2: measurement transformation matrix - 3: X-Ricatti equation - 4: Y-Ricatti equation - TODO: document significance of rcond + 1: control transformation matrix + 2: measurement transformation matrix + 3: X-Ricatti equation + 4: Y-Ricatti equation + TODO: document significance of rcond Raises ------ diff --git a/control/statesp.py b/control/statesp.py index 13e662f03..3a0616c8e 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -68,7 +68,9 @@ __all__ = ['StateSpace', 'ss', 'rss', 'drss', 'tf2ss', 'ssdata'] class StateSpace(LTI): - """A class for representing state-space models + """StateSpace(A, B, C, D[, dt]) + + A class for representing state-space models The StateSpace class is used to represent state-space realizations of linear time-invariant (LTI) systems: @@ -88,15 +90,19 @@ class StateSpace(LTI): means the system timebase is not specified. If 'dt' is set to True, the system will be treated as a discrete time system with unspecified sampling time. - """ def __init__(self, *args): - """Construct a state space object. + """ + StateSpace(A, B, C, D[, dt]) - The default constructor is StateSpace(A, B, C, D), where A, B, C, D are - matrices or equivalent objects. To call the copy constructor, call - StateSpace(sys), where sys is a StateSpace object. + Construct a state space object. + + The default constructor is StateSpace(A, B, C, D), where A, B, C, D + are matrices or equivalent objects. To create a discrete time system, + use StateSpace(A, B, C, D, dt) where 'dt' is the sampling time (or + True for unspecified sampling time). To call the copy constructor, + call StateSpace(sys), where sys is a StateSpace object. """ @@ -110,8 +116,7 @@ def __init__(self, *args): elif len(args) == 1: # Use the copy constructor. if not isinstance(args[0], StateSpace): - raise TypeError("The one-argument constructor can only take in \ -a StateSpace object. Recived %s." % type(args[0])) + raise TypeError("The one-argument constructor can only take in a StateSpace object. Received %s." % type(args[0])) A = args[0].A B = args[0].B C = args[0].C @@ -956,7 +961,8 @@ def _mimo2simo(sys, input, warn_conversion=False): return sys def ss(*args): - """ + """ss(A, B, C, D[, dt]) + Create a state space system. The function accepts either 1, 4 or 5 parameters: @@ -1014,6 +1020,7 @@ def ss(*args): See Also -------- + StateSpace tf ss2tf tf2ss diff --git a/control/statesp.py+ b/control/statesp.py+ new file mode 100644 index 000000000..1007feb0d --- /dev/null +++ b/control/statesp.py+ @@ -0,0 +1,1211 @@ +"""statesp.py + +State space representation and functions. + +This file contains the StateSpace class, which is used to represent linear +systems in state space. This is the primary representation for the +python-control library. + +""" + +# Python 3 compatibility (needs to go here) +from __future__ import print_function +from __future__ import division # for _convertToStateSpace + +"""Copyright (c) 2010 by California Institute of Technology +All rights reserved. + +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. + +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 California Institute of Technology 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: Richard M. Murray +Date: 24 May 09 +Revised: Kevin K. Chen, Dec 10 + +$Id$ +""" + +import numpy as np +from numpy import all, angle, any, array, asarray, concatenate, cos, delete, \ + dot, empty, exp, eye, matrix, ones, pi, poly, poly1d, roots, shape, sin, \ + zeros, squeeze +from numpy.random import rand, randn +from numpy.linalg import solve, eigvals, matrix_rank +from numpy.linalg.linalg import LinAlgError +from scipy.signal import lti, cont2discrete +# from exceptions import Exception +import warnings +from .lti import LTI, timebase, timebaseEqual, isdtime +from .xferfcn import _convertToTransferFunction +from copy import deepcopy + +__all__ = ['StateSpace', 'ss', 'rss', 'drss', 'tf2ss', 'ssdata'] + +class StateSpace(LTI): + """StateSpace(A, B, C, D[, dt]) + + A class for representing state-space models + + The StateSpace class is used to represent state-space realizations of linear + time-invariant (LTI) systems: + + dx/dt = A x + B u + y = C x + D u + + where u is the input, y is the output, and x is the state. + + The main data members are the A, B, C, and D matrices. The class also + keeps track of the number of states (i.e., the size of A). + + Discrete-time state space system are implemented by using the 'dt' instance + variable and setting it to the sampling period. If 'dt' is not None, + then it must match whenever two state space systems are combined. + Setting dt = 0 specifies a continuous system, while leaving dt = None + means the system timebase is not specified. If 'dt' is set to True, the + system will be treated as a discrete time system with unspecified + sampling time. + """ + + def __init__(self, *args): + """ + StateSpace(A, B, C, D[, dt]) + + Construct a state space object. + + The default constructor is StateSpace(A, B, C, D), where A, B, C, D + are matrices or equivalent objects. To create a discrete time system, + use StateSpace(A, B, C, D, dt) where 'dt' is the sampling time (or + True for unspecified sampling time). To call the copy constructor, + call StateSpace(sys), where sys is a StateSpace object. + + """ + + if len(args) == 4: + # The user provided A, B, C, and D matrices. + (A, B, C, D) = args + dt = None; + elif len(args) == 5: + # Discrete time system + (A, B, C, D, dt) = args + elif len(args) == 1: + # Use the copy constructor. + if not isinstance(args[0], StateSpace): + raise TypeError("The one-argument constructor can only take in a StateSpace object. Received %s." % type(args[0])) + A = args[0].A + B = args[0].B + C = args[0].C + D = args[0].D + try: + dt = args[0].dt + except NameError: + dt = None; + else: + raise ValueError("Needs 1 or 4 arguments; received %i." % len(args)) + + A, B, C, D = map(matrix, [A, B, C, D]) + + # TODO: use super here? + LTI.__init__(self, inputs=D.shape[1], outputs=D.shape[0], dt=dt) + self.A = A + self.B = B + self.C = C + self.D = D + + self.states = A.shape[1] + + if 0 == self.states: + # static gain + # matrix's default "empty" shape is 1x0 + A.shape = (0,0) + B.shape = (0,self.inputs) + C.shape = (self.outputs,0) + + # Check that the matrix sizes are consistent. + if self.states != A.shape[0]: + raise ValueError("A must be square.") + if self.states != B.shape[0]: + raise ValueError("A and B must have the same number of rows.") + if self.states != C.shape[1]: + raise ValueError("A and C must have the same number of columns.") + if self.inputs != B.shape[1]: + raise ValueError("B and D must have the same number of columns.") + if self.outputs != C.shape[0]: + raise ValueError("C and D must have the same number of rows.") + + # Check for states that don't do anything, and remove them. + self._remove_useless_states() + + def _remove_useless_states(self): + """Check for states that don't do anything, and remove them. + + Scan the A, B, and C matrices for rows or columns of zeros. If the + zeros are such that a particular state has no effect on the input-output + dynamics, then remove that state from the A, B, and C matrices. + + """ + + # Indices of useless states. + useless = [] + + # Search for useless states. + for i in range(self.states): + if (all(self.A[i, :] == zeros((1, self.states))) and + all(self.B[i, :] == zeros((1, self.inputs)))): + useless.append(i) + # To avoid duplucate indices in useless, jump to the next + # iteration. + continue + if (all(self.A[:, i] == zeros((self.states, 1))) and + all(self.C[:, i] == zeros((self.outputs, 1)))): + useless.append(i) + + # Remove the useless states. + self.A = delete(self.A, useless, 0) + self.A = delete(self.A, useless, 1) + self.B = delete(self.B, useless, 0) + self.C = delete(self.C, useless, 1) + + self.states = self.A.shape[0] + self.inputs = self.B.shape[1] + self.outputs = self.C.shape[0] + + def __str__(self): + """String representation of the state space.""" + + str = "A = " + self.A.__str__() + "\n\n" + str += "B = " + self.B.__str__() + "\n\n" + str += "C = " + self.C.__str__() + "\n\n" + str += "D = " + self.D.__str__() + "\n" + #! TODO: replace with standard calls to lti functions + if (type(self.dt) == bool and self.dt == True): + str += "\ndt unspecified\n" + elif (not (self.dt is None) and type(self.dt) != bool and self.dt > 0): + str += "\ndt = " + self.dt.__str__() + "\n" + return str + + # represent as string, makes display work for IPython + __repr__ = __str__ + + # Negation of a system + def __neg__(self): + """Negate a state space system.""" + + return StateSpace(self.A, self.B, -self.C, -self.D, self.dt) + + # Addition of two state space systems (parallel interconnection) + def __add__(self, other): + """Add two LTI systems (parallel connection).""" + + # Check for a couple of special cases + if (isinstance(other, (int, float, complex))): + # Just adding a scalar; put it in the D matrix + A, B, C = self.A, self.B, self.C; + D = self.D + other; + dt = self.dt + else: + other = _convertToStateSpace(other) + + # Check to make sure the dimensions are OK + if ((self.inputs != other.inputs) or + (self.outputs != other.outputs)): + raise ValueError("Systems have different shapes.") + + # Figure out the sampling time to use + if (self.dt == None and other.dt != None): + dt = other.dt # use dt from second argument + elif (other.dt == None and self.dt != None) or \ + (timebaseEqual(self, other)): + dt = self.dt # use dt from first argument + else: + raise ValueError("Systems have different sampling times") + + # Concatenate the various arrays + A = concatenate(( + concatenate((self.A, zeros((self.A.shape[0], + other.A.shape[-1]))),axis=1), + concatenate((zeros((other.A.shape[0], self.A.shape[-1])), + other.A),axis=1) + ),axis=0) + B = concatenate((self.B, other.B), axis=0) + C = concatenate((self.C, other.C), axis=1) + D = self.D + other.D + + return StateSpace(A, B, C, D, dt) + + # Right addition - just switch the arguments + def __radd__(self, other): + """Right add two LTI systems (parallel connection).""" + + return self + other + + # Subtraction of two state space systems (parallel interconnection) + def __sub__(self, other): + """Subtract two LTI systems.""" + + return self + (-other) + + def __rsub__(self, other): + """Right subtract two LTI systems.""" + + return other + (-self) + + # Multiplication of two state space systems (series interconnection) + def __mul__(self, other): + """Multiply two LTI objects (serial connection).""" + + # Check for a couple of special cases + if isinstance(other, (int, float, complex)): + # Just multiplying by a scalar; change the output + A, B = self.A, self.B + C = self.C * other + D = self.D * other + dt = self.dt + else: + other = _convertToStateSpace(other) + + # Check to make sure the dimensions are OK + if self.inputs != other.outputs: + raise ValueError("C = A * B: A has %i column(s) (input(s)), \ +but B has %i row(s)\n(output(s))." % (self.inputs, other.outputs)) + + # Figure out the sampling time to use + if (self.dt == None and other.dt != None): + dt = other.dt # use dt from second argument + elif (other.dt == None and self.dt != None) or \ + (timebaseEqual(self, other)): + dt = self.dt # use dt from first argument + else: + raise ValueError("Systems have different sampling times") + + # Concatenate the various arrays + A = concatenate( + (concatenate((other.A, zeros((other.A.shape[0], self.A.shape[1]))), + axis=1), + concatenate((self.B * other.C, self.A), axis=1)), axis=0) + B = concatenate((other.B, self.B * other.D), axis=0) + C = concatenate((self.D * other.C, self.C),axis=1) + D = self.D * other.D + + return StateSpace(A, B, C, D, dt) + + # Right multiplication of two state space systems (series interconnection) + # Just need to convert LH argument to a state space object + # TODO: __rmul__ only works for special cases (??) + def __rmul__(self, other): + """Right multiply two LTI objects (serial connection).""" + + # Check for a couple of special cases + if isinstance(other, (int, float, complex)): + # Just multiplying by a scalar; change the input + A, C = self.A, self.C; + B = self.B * other; + D = self.D * other; + return StateSpace(A, B, C, D, self.dt) + + # is lti, and convertible? + if isinstance(other, LTI): + return _convertToStateSpace(other) * self + + # try to treat this as a matrix + try: + X = matrix(other) + C = X * self.C + D = X * self.D + return StateSpace(self.A, self.B, C, D, self.dt) + + except Exception as e: + print(e) + pass + raise TypeError("can't interconnect systems") + + # TODO: __div__ and __rdiv__ are not written yet. + def __div__(self, other): + """Divide two LTI systems.""" + + raise NotImplementedError("StateSpace.__div__ is not implemented yet.") + + def __rdiv__(self, other): + """Right divide two LTI systems.""" + + raise NotImplementedError("StateSpace.__rdiv__ is not implemented yet.") + + # TODO: add discrete time check + def evalfr(self, omega): + """Evaluate a SS system's transfer function at a single frequency. + + self.evalfr(omega) returns the value of the transfer function matrix with + input value s = i * omega. + + """ + # Figure out the point to evaluate the transfer function + if isdtime(self, strict=True): + dt = timebase(self) + s = exp(1.j * omega * dt) + if (omega * dt > pi): + warnings.warn("evalfr: frequency evaluation above Nyquist frequency") + else: + s = omega * 1.j + + return self.horner(s) + + def horner(self, s): + '''Evaluate the systems's transfer function for a complex variable + + Returns a matrix of values evaluated at complex variable s. + ''' + resp = self.C * solve(s * eye(self.states) - self.A, + self.B) + self.D + return array(resp) + + # Method for generating the frequency response of the system + def freqresp(self, omega): + """Evaluate the system's transfer func. at a list of ang. frequencies. + + mag, phase, omega = self.freqresp(omega) + + reports the value of the magnitude, phase, and angular frequency of the + system's transfer function matrix evaluated at s = i * omega, where + omega is a list of angular frequencies, and is a sorted version of the + input omega. + + """ + # when evaluating at many frequencies, much faster to convert to + # transfer function first and then evaluate, than to solve an + # n-dimensional linear system at each frequency + tf = _convertToTransferFunction(self) + return tf.freqresp(omega) + + # Compute poles and zeros + def pole(self): + """Compute the poles of a state space system.""" + + return eigvals(self.A) if self.states else np.array([]) + + def zero(self): + """Compute the zeros of a state space system.""" + + if self.inputs > 1 or self.outputs > 1: + raise NotImplementedError("StateSpace.zeros is currently \ +implemented only for SISO systems.") + + den = poly1d(poly(self.A)) + # Compute the numerator based on zeros + #! TODO: This is currently limited to SISO systems + num = poly1d(poly(self.A - dot(self.B, self.C)) + ((self.D[0, 0] - 1) * + den)) + + return roots(num) + + # Feedback around a state space system + def feedback(self, other=1, sign=-1): + """Feedback interconnection between two LTI systems.""" + + other = _convertToStateSpace(other) + + # Check to make sure the dimensions are OK + if ((self.inputs != other.outputs) or (self.outputs != other.inputs)): + raise ValueError("State space systems don't have compatible \ +inputs/outputs for feedback.") + + # Figure out the sampling time to use + if (self.dt == None and other.dt != None): + dt = other.dt # use dt from second argument + elif (other.dt == None and self.dt != None) or \ + timebaseEqual(self, other): + dt = self.dt # use dt from first argument + else: + raise ValueError("Systems have different sampling times") + + A1 = self.A + B1 = self.B + C1 = self.C + D1 = self.D + A2 = other.A + B2 = other.B + C2 = other.C + D2 = other.D + + F = eye(self.inputs) - sign * D2 * D1 + if matrix_rank(F) != self.inputs: + raise ValueError("I - sign * D2 * D1 is singular to working precision.") + + # Precompute F\D2 and F\C2 (E = inv(F)) + # We can solve two linear systems in one pass, since the + # coefficients matrix F is the same. Thus, we perform the LU + # decomposition (cubic runtime complexity) of F only once! + # The remaining back substitutions are only quadratic in runtime. + E_D2_C2 = solve(F, concatenate((D2, C2), axis=1)) + E_D2 = E_D2_C2[:, :other.inputs] + E_C2 = E_D2_C2[:, other.inputs:] + + T1 = eye(self.outputs) + sign * D1 * E_D2 + T2 = eye(self.inputs) + sign * E_D2 * D1 + + A = concatenate( + (concatenate( + (A1 + sign * B1 * E_D2 * C1, sign * B1 * E_C2), axis=1), + concatenate( + (B2 * T1 * C1, A2 + sign * B2 * D1 * E_C2), axis=1)), + axis=0) + B = concatenate((B1 * T2, B2 * D1 * T2), axis=0) + C = concatenate((T1 * C1, sign * D1 * E_C2), axis=1) + D = D1 * T2 + + return StateSpace(A, B, C, D, dt) + + def minreal(self, tol=0.0): + """Calculate a minimal realization, removes unobservable and + uncontrollable states""" + if self.states: + try: + from slycot import tb01pd + B = empty((self.states, max(self.inputs, self.outputs))) + B[:,:self.inputs] = self.B + C = empty((max(self.outputs, self.inputs), self.states)) + C[:self.outputs,:] = self.C + A, B, C, nr = tb01pd(self.states, self.inputs, self.outputs, + self.A, B, C, tol=tol) + return StateSpace(A[:nr,:nr], B[:nr,:self.inputs], + C[:self.outputs,:nr], self.D) + except ImportError: + raise TypeError("minreal requires slycot tb01pd") + else: + return StateSpace(self) + + + # TODO: add discrete time check + def returnScipySignalLTI(self): + """Return a list of a list of scipy.signal.lti objects. + + For instance, + + >>> out = ssobject.returnScipySignalLTI() + >>> out[3][5] + + is a signal.scipy.lti object corresponding to the transfer function from + the 6th input to the 4th output.""" + + # Preallocate the output. + out = [[[] for j in range(self.inputs)] for i in range(self.outputs)] + + for i in range(self.outputs): + for j in range(self.inputs): + out[i][j] = lti(asarray(self.A), asarray(self.B[:, j]), + asarray(self.C[i, :]), asarray(self.D[i, j])) + + return out + + def append(self, other): + """Append a second model to the present model. The second + model is converted to state-space if necessary, inputs and + outputs are appended and their order is preserved""" + if not isinstance(other, StateSpace): + other = _convertToStateSpace(other) + + if self.dt != other.dt: + raise ValueError("Systems must have the same time step") + + n = self.states + other.states + m = self.inputs + other.inputs + p = self.outputs + other.outputs + A = zeros( (n, n) ) + B = zeros( (n, m) ) + C = zeros( (p, n) ) + D = zeros( (p, m) ) + A[:self.states,:self.states] = self.A + A[self.states:,self.states:] = other.A + B[:self.states,:self.inputs] = self.B + B[self.states:,self.inputs:] = other.B + C[:self.outputs,:self.states] = self.C + C[self.outputs:,self.states:] = other.C + D[:self.outputs,:self.inputs] = self.D + D[self.outputs:,self.inputs:] = other.D + return StateSpace(A, B, C, D, self.dt) + + def __getitem__(self, indices): + """Array style acces""" + if len(indices) != 2: + raise IOError('must provide indices of length 2 for state space') + i = indices[0] + j = indices[1] + return StateSpace(self.A, + self.B[:,j], + self.C[i,:], + self.D[i,j], self.dt) + + def sample(self, Ts, method='zoh', alpha=None): + """Convert a continuous time system to discrete time + + Creates a discrete-time system from a continuous-time system by + sampling. Multiple methods of conversion are supported. + + Parameters + ---------- + Ts : float + Sampling period + method : {"gbt", "bilinear", "euler", "backward_diff", "zoh"} + Which method to use: + + * gbt: generalized bilinear transformation + * bilinear: Tustin's approximation ("gbt" with alpha=0.5) + * euler: Euler (or forward differencing) method ("gbt" with alpha=0) + * backward_diff: Backwards differencing ("gbt" with alpha=1.0) + * zoh: zero-order hold (default) + + alpha : float within [0, 1] + The generalized bilinear transformation weighting parameter, which + should only be specified with method="gbt", and is ignored otherwise + + Returns + ------- + sysd : StateSpace system + Discrete time system, with sampling rate Ts + + Notes + ----- + Uses the command 'cont2discrete' from scipy.signal + + Examples + -------- + >>> sys = StateSpace(0, 1, 1, 0) + >>> sysd = sys.sample(0.5, method='bilinear') + + """ + if not self.isctime(): + raise ValueError("System must be continuous time system") + + sys = (self.A, self.B, self.C, self.D) + Ad, Bd, C, D, dt = cont2discrete(sys, Ts, method, alpha) + return StateSpace(Ad, Bd, C, D, dt) + + def dcgain(self): + """Return the zero-frequency gain + + The zero-frequency gain of a continuous-time state-space + system is given by: + + .. math: G(0) = - C A^{-1} B + D + + and of a discrete-time state-space system by: + + .. math: G(1) = C (I - A)^{-1} B + D + + Returns + ------- + gain : ndarray + An array of shape (outputs,inputs); the array will either + be the zero-frequency (or DC) gain, or, if the frequency + response is singular, the array will be filled with np.nan. + """ + try: + if self.isctime(): + gain = np.asarray(self.D - + self.C.dot(np.linalg.solve(self.A, self.B))) + else: + gain = self.horner(1) + except LinAlgError: + # eigenvalue at DC + gain = np.tile(np.nan,(self.outputs,self.inputs)) + return np.squeeze(gain) + + +# TODO: add discrete time check +def _convertToStateSpace(sys, **kw): + """Convert a system to state space form (if needed). + + 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 = _convertToStateSpace(3.) # Assumes inputs = outputs = 1 + >>> sys = _convertToStateSpace(1., inputs=3, outputs=2) + + In the latter example, A = B = C = 0 and D = [[1., 1., 1.] + [1., 1., 1.]]. + + """ + + from .xferfcn import TransferFunction + import itertools + if isinstance(sys, StateSpace): + if len(kw): + raise TypeError("If sys is a StateSpace, _convertToStateSpace \ +cannot take keywords.") + + # Already a state space system; just return it + return sys + elif isinstance(sys, TransferFunction): + try: + from slycot import td04ad + if len(kw): + raise TypeError("If sys is a TransferFunction, _convertToStateSpace \ + cannot take keywords.") + + # Change the numerator and denominator arrays so that the transfer + # function matrix has a common denominator. + num, den = sys._common_den() + # Make a list of the orders of the denominator polynomials. + index = [len(den) - 1 for i in range(sys.outputs)] + # Repeat the common denominator along the rows. + den = array([den for i in range(sys.outputs)]) + #! TODO: transfer function to state space conversion is still buggy! + #print num + #print shape(num) + ssout = td04ad('R',sys.inputs, sys.outputs, index, den, num,tol=0.0) + + states = ssout[0] + return StateSpace(ssout[1][:states, :states], + ssout[2][:states, :sys.inputs], + ssout[3][:sys.outputs, :states], + ssout[4], sys.dt) + except ImportError: + # No Slycot. Scipy tf->ss can't handle MIMO, but static + # MIMO is an easy special case we can check for here + maxn = max(max(len(n) for n in nrow) + for nrow in sys.num) + maxd = max(max(len(d) for d in drow) + for drow in sys.den) + if 1==maxn and 1==maxd: + D = empty((sys.outputs,sys.inputs),dtype=float) + for i,j in itertools.product(range(sys.outputs),range(sys.inputs)): + D[i,j] = sys.num[i][j][0] / sys.den[i][j][0] + return StateSpace([], [], [], D, sys.dt) + else: + if (sys.inputs != 1 or sys.outputs != 1): + raise TypeError("No support for MIMO without slycot") + + # TODO: do we want to squeeze first and check dimenations? + # I think this will fail if num and den aren't 1-D after + # the squeeze + lti_sys = lti(squeeze(sys.num), squeeze(sys.den)) + return StateSpace(lti_sys.A, lti_sys.B, lti_sys.C, lti_sys.D, + sys.dt) + + elif isinstance(sys, (int, float, complex)): + if "inputs" in kw: + inputs = kw["inputs"] + else: + inputs = 1 + if "outputs" in kw: + outputs = kw["outputs"] + else: + outputs = 1 + + # Generate a simple state space system of the desired dimension + # The following Doesn't work due to inconsistencies in ltisys: + # return StateSpace([[]], [[]], [[]], eye(outputs, inputs)) + return StateSpace(0., zeros((1, inputs)), zeros((outputs, 1)), + sys * ones((outputs, inputs))) + + # If this is a matrix, try to create a constant feedthrough + try: + D = matrix(sys) + outputs, inputs = D.shape + + return StateSpace(0., zeros((1, inputs)), zeros((outputs, 1)), D) + except Exception(e): + print("Failure to assume argument is matrix-like in" \ + " _convertToStateSpace, result %s" % e) + + raise TypeError("Can't convert given type to StateSpace system.") + +# TODO: add discrete time option +def _rss_generate(states, inputs, outputs, type): + """Generate a random state space. + + This does the actual random state space generation expected from rss and + drss. type is 'c' for continuous systems and 'd' for discrete systems. + + """ + + # Probability of repeating a previous root. + pRepeat = 0.05 + # Probability of choosing a real root. Note that when choosing a complex + # root, the conjugate gets chosen as well. So the expected proportion of + # real roots is pReal / (pReal + 2 * (1 - pReal)). + pReal = 0.6 + # Probability that an element in B or C will not be masked out. + pBCmask = 0.8 + # Probability that an element in D will not be masked out. + pDmask = 0.3 + # Probability that D = 0. + pDzero = 0.5 + + # Check for valid input arguments. + if states < 1 or states % 1: + raise ValueError("states must be a positive integer. states = %g." % + states) + if inputs < 1 or inputs % 1: + raise ValueError("inputs must be a positive integer. inputs = %g." % + inputs) + if outputs < 1 or outputs % 1: + raise ValueError("outputs must be a positive integer. outputs = %g." % + outputs) + + # Make some poles for A. Preallocate a complex array. + poles = zeros(states) + zeros(states) * 0.j + i = 0 + + while i < states: + if rand() < pRepeat and i != 0 and i != states - 1: + # Small chance of copying poles, if we're not at the first or last + # element. + if poles[i-1].imag == 0: + # Copy previous real pole. + poles[i] = poles[i-1] + i += 1 + else: + # Copy previous complex conjugate pair of poles. + poles[i:i+2] = poles[i-2:i] + i += 2 + elif rand() < pReal or i == states - 1: + # No-oscillation pole. + if type == 'c': + poles[i] = -exp(randn()) + 0.j + elif type == 'd': + poles[i] = 2. * rand() - 1. + i += 1 + else: + # Complex conjugate pair of oscillating poles. + if type == 'c': + poles[i] = complex(-exp(randn()), 3. * exp(randn())) + elif type == 'd': + mag = rand() + phase = 2. * pi * rand() + poles[i] = complex(mag * cos(phase), + mag * sin(phase)) + poles[i+1] = complex(poles[i].real, -poles[i].imag) + i += 2 + + # Now put the poles in A as real blocks on the diagonal. + A = zeros((states, states)) + i = 0 + while i < states: + if poles[i].imag == 0: + A[i, i] = poles[i].real + i += 1 + else: + A[i, i] = A[i+1, i+1] = poles[i].real + A[i, i+1] = poles[i].imag + A[i+1, i] = -poles[i].imag + i += 2 + # Finally, apply a transformation so that A is not block-diagonal. + while True: + T = randn(states, states) + try: + A = dot(solve(T, A), T) # A = T \ A * T + break + except LinAlgError: + # In the unlikely event that T is rank-deficient, iterate again. + pass + + # Make the remaining matrices. + B = randn(states, inputs) + C = randn(outputs, states) + D = randn(outputs, inputs) + + # Make masks to zero out some of the elements. + while True: + Bmask = rand(states, inputs) < pBCmask + if any(Bmask): # Retry if we get all zeros. + break + while True: + Cmask = rand(outputs, states) < pBCmask + if any(Cmask): # Retry if we get all zeros. + break + if rand() < pDzero: + Dmask = zeros((outputs, inputs)) + else: + Dmask = rand(outputs, inputs) < pDmask + + # Apply masks. + B = B * Bmask + C = C * Cmask + D = D * Dmask + + return StateSpace(A, B, C, D) + +# Convert a MIMO system to a SISO system +# TODO: add discrete time check +def _mimo2siso(sys, input, output, warn_conversion=False): + #pylint: disable=W0622 + """ + Convert a MIMO system to a SISO system. (Convert a system with multiple + inputs and/or outputs, to a system with a single input and output.) + + The input and output that are used in the SISO system can be selected + with the parameters ``input`` and ``output``. All other inputs are set + to 0, all other outputs are ignored. + + If ``sys`` is already a SISO system, it will be returned unaltered. + + Parameters + ---------- + sys: StateSpace + Linear (MIMO) system that should be converted. + input: int + Index of the input that will become the SISO system's only input. + output: int + Index of the output that will become the SISO system's only output. + warn_conversion: bool + If True: print a warning message when sys is a MIMO system. + Warn that a conversion will take place. + + Returns: + + sys: StateSpace + The converted (SISO) system. + """ + if not (isinstance(input, int) and isinstance(output, int)): + raise TypeError("Parameters ``input`` and ``output`` must both " + "be integer numbers.") + if not (0 <= input < sys.inputs): + raise ValueError("Selected input does not exist. " + "Selected input: {sel}, " + "number of system inputs: {ext}." + .format(sel=input, ext=sys.inputs)) + if not (0 <= output < sys.outputs): + raise ValueError("Selected output does not exist. " + "Selected output: {sel}, " + "number of system outputs: {ext}." + .format(sel=output, ext=sys.outputs)) + #Convert sys to SISO if necessary + if sys.inputs > 1 or sys.outputs > 1: + if warn_conversion: + warnings.warn("Converting MIMO system to SISO system. " + "Only input {i} and output {o} are used." + .format(i=input, o=output)) + # $X = A*X + B*U + # Y = C*X + D*U + new_B = sys.B[:, input] + new_C = sys.C[output, :] + new_D = sys.D[output, input] + sys = StateSpace(sys.A, new_B, new_C, new_D, sys.dt) + + return sys + +def _mimo2simo(sys, input, warn_conversion=False): + #pylint: disable=W0622 + """ + Convert a MIMO system to a SIMO system. (Convert a system with multiple + inputs and/or outputs, to a system with a single input but possibly + multiple outputs.) + + The input that is used in the SIMO system can be selected with the + parameter ``input``. All other inputs are set to 0, all other + outputs are ignored. + + If ``sys`` is already a SIMO system, it will be returned unaltered. + + Parameters + ---------- + sys: StateSpace + Linear (MIMO) system that should be converted. + input: int + Index of the input that will become the SIMO system's only input. + warn_conversion: bool + If True: print a warning message when sys is a MIMO system. + Warn that a conversion will take place. + + Returns: + -------- + sys: StateSpace + The converted (SIMO) system. + """ + if not (isinstance(input, int)): + raise TypeError("Parameter ``input`` be an integer number.") + if not (0 <= input < sys.inputs): + raise ValueError("Selected input does not exist. " + "Selected input: {sel}, " + "number of system inputs: {ext}." + .format(sel=input, ext=sys.inputs)) + #Convert sys to SISO if necessary + if sys.inputs > 1: + if warn_conversion: + warnings.warn("Converting MIMO system to SIMO system. " + "Only input {i} is used." + .format(i=input)) + # $X = A*X + B*U + # Y = C*X + D*U + new_B = sys.B[:, input] + new_D = sys.D[:, input] + sys = StateSpace(sys.A, new_B, sys.C, new_D, sys.dt) + + return sys + +def ss(*args): + """ + Create a state space system. + + The function accepts either 1, 4 or 5 parameters: + + ``ss(sys)`` + Convert a linear system into space system form. Always creates a + new system, even if sys is already a StateSpace object. + + ``ss(A, B, C, D)`` + Create a state space system from the matrices of its state and + output equations: + + .. math:: + \dot x = A \cdot x + B \cdot u + + y = C \cdot x + D \cdot u + + ``ss(A, B, C, D, dt)`` + Create a discrete-time state space system from the matrices of + its state and output equations: + + .. math:: + x[k+1] = A \cdot x[k] + B \cdot u[k] + + y[k] = C \cdot x[k] + D \cdot u[ki] + + The matrices can be given as *array like* data types or strings. + Everything that the constructor of :class:`numpy.matrix` accepts is + permissible here too. + + Parameters + ---------- + sys: StateSpace or TransferFunction + A linear system + A: array_like or string + System matrix + B: array_like or string + Control matrix + C: array_like or string + Output matrix + D: array_like or string + Feed forward matrix + dt: If present, specifies the sampling period and a discrete time + system is created + + Returns + ------- + out: :class:`StateSpace` + The new linear system + + Raises + ------ + ValueError + if matrix sizes are not self-consistent + + See Also + -------- + tf + ss2tf + tf2ss + + Examples + -------- + >>> # Create a StateSpace object from four "matrices". + >>> sys1 = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.") + + >>> # Convert a TransferFunction to a StateSpace object. + >>> sys_tf = tf([2.], [1., 3]) + >>> sys2 = ss(sys_tf) + + """ + + if len(args) == 4 or len(args) == 5: + return StateSpace(*args) + elif len(args) == 1: + from .xferfcn import TransferFunction + sys = args[0] + if isinstance(sys, StateSpace): + return deepcopy(sys) + elif isinstance(sys, TransferFunction): + return tf2ss(sys) + else: + raise TypeError("ss(sys): sys must be a StateSpace or \ +TransferFunction object. It is %s." % type(sys)) + else: + raise ValueError("Needs 1 or 4 arguments; received %i." % len(args)) + +def tf2ss(*args): + """ + Transform a transfer function to a state space system. + + The function accepts either 1 or 2 parameters: + + ``tf2ss(sys)`` + Convert a linear system into transfer function form. Always creates + a new system, even if sys is already a TransferFunction object. + + ``tf2ss(num, den)`` + Create a transfer function system from its numerator and denominator + polynomial coefficients. + + For details see: :func:`tf` + + Parameters + ---------- + sys: LTI (StateSpace or TransferFunction) + A linear system + num: array_like, or list of list of array_like + Polynomial coefficients of the numerator + den: array_like, or list of list of array_like + Polynomial coefficients of the denominator + + Returns + ------- + out: StateSpace + New linear system in state space form + + Raises + ------ + ValueError + if `num` and `den` have invalid or unequal dimensions, or if an + invalid number of arguments is passed in + TypeError + if `num` or `den` are of incorrect type, or if sys is not a + TransferFunction object + + See Also + -------- + ss + tf + ss2tf + + Examples + -------- + >>> 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) + + >>> sys_tf = tf(num, den) + >>> sys2 = tf2ss(sys_tf) + + """ + + from .xferfcn import TransferFunction + if len(args) == 2 or len(args) == 3: + # Assume we were given the num, den + return _convertToStateSpace(TransferFunction(*args)) + + elif len(args) == 1: + sys = args[0] + if not isinstance(sys, TransferFunction): + raise TypeError("tf2ss(sys): sys must be a TransferFunction \ +object.") + return _convertToStateSpace(sys) + else: + raise ValueError("Needs 1 or 2 arguments; received %i." % len(args)) + +def rss(states=1, outputs=1, inputs=1): + """ + Create a stable **continuous** random state space object. + + Parameters + ---------- + states: integer + Number of state variables + inputs: integer + Number of system inputs + outputs: integer + Number of system outputs + + Returns + ------- + sys: StateSpace + The randomly created linear system + + Raises + ------ + ValueError + if any input is not a positive integer + + See Also + -------- + drss + + Notes + ----- + If the number of states, inputs, or outputs is not specified, then the + missing numbers are assumed to be 1. The poles of the returned system + will always have a negative real part. + + """ + + return _rss_generate(states, inputs, outputs, 'c') + +def drss(states=1, outputs=1, inputs=1): + """ + Create a stable **discrete** random state space object. + + Parameters + ---------- + states: integer + Number of state variables + inputs: integer + Number of system inputs + outputs: integer + Number of system outputs + + Returns + ------- + sys: StateSpace + The randomly created linear system + + Raises + ------ + ValueError + if any input is not a positive integer + + See Also + -------- + rss + + Notes + ----- + If the number of states, inputs, or outputs is not specified, then the + missing numbers are assumed to be 1. The poles of the returned system + will always have a magnitude less than 1. + + """ + + return _rss_generate(states, inputs, outputs, 'd') + +def ssdata(sys): + ''' + Return state space data objects for a system + + Parameters + ---------- + sys: LTI (StateSpace, or TransferFunction) + LTI system whose data will be returned + + Returns + ------- + (A, B, C, D): list of matrices + State space data for the system + ''' + ss = _convertToStateSpace(sys) + return (ss.A, ss.B, ss.C, ss.D) diff --git a/control/xferfcn.py b/control/xferfcn.py index be21f8509..c89b3254c 100644 --- a/control/xferfcn.py +++ b/control/xferfcn.py @@ -65,7 +65,9 @@ class TransferFunction(LTI): - """A class for representing transfer functions + """TransferFunction(num, den[, dt]) + + A class for representing transfer functions The TransferFunction class is used to represent systems in transfer function form. @@ -87,13 +89,17 @@ class TransferFunction(LTI): """ def __init__(self, *args): - """Construct a transfer function. + """TransferFunction(num, den[, dt]) + + Construct a transfer function. The default constructor is TransferFunction(num, den), where num and den are lists of lists of arrays containing polynomial coefficients. - To crete a discrete time transfer funtion, use TransferFunction(num, - den, dt). To call the copy constructor, call TransferFunction(sys), - where sys is a TransferFunction object (continuous or discrete). + To create a discrete time transfer funtion, use TransferFunction(num, + den, dt) where 'dt' is the sampling time (or True for unspecified + sampling time). To call the copy constructor, call + TransferFunction(sys), where sys is a TransferFunction object + (continuous or discrete). """ @@ -1173,10 +1179,11 @@ def _convertToTransferFunction(sys, **kw): def tf(*args): - """ + """tf(num, den[, dt]) + Create a transfer function system. Can create MIMO systems. - The function accepts either 1 or 2 parameters: + The function accepts either 1, 2, or 3 parameters: ``tf(sys)`` Convert a linear system into transfer function form. Always creates @@ -1221,6 +1228,7 @@ def tf(*args): See Also -------- + TransferFunction ss ss2tf tf2ss From bde823b640480e2bc1d66b4f7c06beb531cad37f Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Thu, 28 Dec 2017 21:00:36 -0800 Subject: [PATCH 13/21] DOC: add more info on class constructors --- doc/control.rst | 16 ++++---- doc/conventions.rst | 91 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 100 insertions(+), 7 deletions(-) diff --git a/doc/control.rst b/doc/control.rst index f3d124ee5..294701f8f 100644 --- a/doc/control.rst +++ b/doc/control.rst @@ -24,14 +24,14 @@ System creation System interconnections ======================= .. autosummary:: - :toctree: generated/ + :toctree: generated/ - append - connect - feedback - negate - parallel - series + append + connect + feedback + negate + parallel + series Frequency domain plotting ========================= @@ -118,6 +118,8 @@ Model simplification tools era markov +.. _utility-and-conversions: + Utility functions and conversions ================================= .. autosummary:: diff --git a/doc/conventions.rst b/doc/conventions.rst index 881a0a6bb..77d6a30dc 100644 --- a/doc/conventions.rst +++ b/doc/conventions.rst @@ -9,6 +9,97 @@ Library conventions The python-control library uses a set of standard conventions for the way that different types of standard information used by the library. +LTI system representation +========================= + +Linear time invariant (LTI) systems are represented in python-control in +state space, transfer function, or frequency response data (FRD) form. Most +functions in the toolbox will operate on any of these data types and +functions for convering between between compatible types is provided. + +State space systems +------------------- +The :class:`StateSpace` class is used to represent state-space realizations +of linear time-invariant (LTI) systems: + +.. math:: + + \frac{dx}{dt} &= A x + B u \\ + y &= C x + D u + +where u is the input, y is the output, and x is the state. + +To create a state space system, use the :class:`StateSpace` constructor: + + sys = StateSpace(A, B, C, D) + +State space systems can be manipulated using standard arithmetic operations +as well as the :func:`feedback`, :func:`parallel`, and :func:`series` +function. A full list of functions can be found in :ref:`function-ref`. + +Transfer functions +------------------ +The :class:`TransferFunction` class is used to represent input/output +transfer functions + +.. math:: + + G(s) = \frac{\text{num}(s)}{\text{den(s)}} + = \frac{a_0 s^n + a_1 s^{n-1} + \cdots + a_n} + {b_0 s^m + b_1 s^{m-1} + \cdots + b_m}, + +where n is generally greater than m (for a proper transfer function). + +To create a transfer function, use the :class:`TransferFunction` +constructor: + + sys = TransferFunction(num, den) + +Transfer functions can be manipulated using standard arithmetic operations +as well as the :func:`feedback`, :func:`parallel`, and :func:`series` +function. A full list of functions can be found in :ref:`function-ref`. + +FRD (frequency response data) systems +------------------------------------- +The :class:`FRD` class is used to represent systems in frequency response +data form. + +The main data members are `omega` and `fresp`, where `omega` is a 1D array +with the frequency points of the response, and `fresp` is a 3D array, with +the first dimension corresponding to the output index of the FRD, the second +dimension corresponding to the input index, and the 3rd dimension +corresponding to the frequency points in omega. + +FRD systems have a somewhat more limited set of functions that are +available, although all of the standard algebraic manipulations can be +performed. + +Discrete time systems +--------------------- +By default, all systems are considered to be continuous time systems. A +discrete time system is created by specifying the 'time base' dt. The time +base argument can be given when a system is constructed: + +* dt = None: continuous time +* dt = number: discrete time system with sampling period 'dt' +* dt = True: discrete time with unspecified sampling period + +Only the :class:`StateSpace` and :class:`TransferFunction` classes allow +explicit representation of discrete time systems. + +Systems must have the same time base in order to be combined. For +continuous time systems, the :func:`sample_system` function or the +:meth:`StateSpace.sample` and :meth:`TransferFunction.sample` methods can be +used to create a discrete time system from a continuous time system. See +:ref:`utility-and-conversions`. + +Conversion between representations +---------------------------------- +LTI systems can be converted between representations either by calling the +constructor for the desired data type using the original system as the sole +argument or using the explict conversion functions :func:`ss2tf` and +:func:`tf2ss`. + .. _time-series-convention: Time series data From b881b77e485be3c7d90a56a0db8d190b56cdf8c1 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Thu, 28 Dec 2017 21:00:58 -0800 Subject: [PATCH 14/21] DOC: improved docstrings for fucntions with variable arguments --- control/bdalg.py | 3 ++- control/margins.py | 4 +++- control/matlab/wrappers.py | 4 +++- control/statefbk.py | 4 +++- control/statesp.py | 3 ++- control/xferfcn.py | 3 ++- 6 files changed, 15 insertions(+), 6 deletions(-) diff --git a/control/bdalg.py b/control/bdalg.py index 9d5aac654..3b2db6051 100644 --- a/control/bdalg.py +++ b/control/bdalg.py @@ -246,7 +246,8 @@ def feedback(sys1, sys2=1, sign=-1): return sys1.feedback(sys2, sign) def append(*sys): - ''' + '''append(sys1, sys2, ... sysn) + Group models by appending their inputs and outputs Forms an augmented system model, and appends the inputs and diff --git a/control/margins.py b/control/margins.py index 9accd566a..5b1ac9131 100644 --- a/control/margins.py +++ b/control/margins.py @@ -336,7 +336,9 @@ def phase_crossover_frequencies(sys): def margin(*args): - """Calculate gain and phase margins and associated crossover frequencies + """margin(sys) + + Calculate gain and phase margins and associated crossover frequencies Parameters ---------- diff --git a/control/matlab/wrappers.py b/control/matlab/wrappers.py index 8c1cfb005..d83890a33 100644 --- a/control/matlab/wrappers.py +++ b/control/matlab/wrappers.py @@ -10,7 +10,9 @@ __all__ = ['bode', 'ngrid', 'dcgain'] def bode(*args, **keywords): - """Bode plot of the frequency response + """bode(syslist[, omega, dB, Hz, deg, ...]) + + Bode plot of the frequency response Plots a bode gain and phase diagram diff --git a/control/statefbk.py b/control/statefbk.py index 977f83ae4..634922131 100644 --- a/control/statefbk.py +++ b/control/statefbk.py @@ -144,7 +144,9 @@ def acker(A, B, poles): return K def lqr(*args, **keywords): - """Linear quadratic regulator design + """lqr(A, B, Q, R[, N]) + + Linear quadratic regulator design The lqr() function computes the optimal state feedback controller that minimizes the quadratic cost diff --git a/control/statesp.py b/control/statesp.py index 3a0616c8e..70d5e4635 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -1052,7 +1052,8 @@ def ss(*args): raise ValueError("Needs 1 or 4 arguments; received %i." % len(args)) def tf2ss(*args): - """ + """tf2ss(sys) + Transform a transfer function to a state space system. The function accepts either 1 or 2 parameters: diff --git a/control/xferfcn.py b/control/xferfcn.py index c89b3254c..79c9036a4 100644 --- a/control/xferfcn.py +++ b/control/xferfcn.py @@ -1273,7 +1273,8 @@ def tf(*args): raise ValueError("Needs 1 or 2 arguments; received %i." % len(args)) def ss2tf(*args): - """ + """ss2tf(sys) + Transform a state space system to a transfer function. The function accepts either 1 or 4 parameters: From b1437f0e025fa732feb9653ff322c1e80b157d5f Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Fri, 29 Dec 2017 08:17:46 -0800 Subject: [PATCH 15/21] unit test from kangwonlee commit facfe9c49d79227dd3bdcd22af9d0aa917ac6008 --- control/tests/tf_input_element_int.py | 46 +++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) create mode 100644 control/tests/tf_input_element_int.py diff --git a/control/tests/tf_input_element_int.py b/control/tests/tf_input_element_int.py new file mode 100644 index 000000000..13d0310e7 --- /dev/null +++ b/control/tests/tf_input_element_int.py @@ -0,0 +1,46 @@ +import unittest + +import numpy as np + +import control as ctl + + +class TestTfInputIntElement(unittest.TestCase): + # currently these do not pass + def test_tf_den_with_numpy_int_element(self): + num = 1 + den = np.convolve([1, 2, 1], [1, 1, 1]) + + sys = ctl.tf(num, den) + + self.assertAlmostEqual(1.0, ctl.dcgain(sys)) + + def test_tf_num_with_numpy_int_element(self): + num = np.convolve([1], [1, 1]) + den = np.convolve([1, 2, 1], [1, 1, 1]) + + sys = ctl.tf(num, den) + + self.assertAlmostEqual(1.0, ctl.dcgain(sys)) + + # currently these pass + def test_tf_input_with_int_element_works(self): + num = 1 + den = np.convolve([1.0, 2, 1], [1, 1, 1]) + + sys = ctl.tf(num, den) + + self.assertAlmostEqual(1.0, ctl.dcgain(sys)) + + def test_ss_input_with_int_element(self): + ident = np.matrix(np.identity(2), dtype=int) + a = np.matrix([[0, 1], + [-1, -2]], dtype=int) * ident + b = np.matrix([[0], + [1]], dtype=int) + c = np.matrix([[0, 1]], dtype=int) + d = 0 + + sys = ctl.ss(a, b, c, d) + sys2 = ctl.ss2tf(sys) + self.assertAlmostEqual(ctl.dcgain(sys), ctl.dcgain(sys2)) From 10cdb7a6bbfa0b9ed9cdff9fac9ea74046003b4d Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Fri, 29 Dec 2017 08:51:00 -0800 Subject: [PATCH 16/21] TST: renamed unit test from kangwonlee --- ...ut_element_int.py => input_element_int_test.py} | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) rename control/tests/{tf_input_element_int.py => input_element_int_test.py} (81%) diff --git a/control/tests/tf_input_element_int.py b/control/tests/input_element_int_test.py similarity index 81% rename from control/tests/tf_input_element_int.py rename to control/tests/input_element_int_test.py index 13d0310e7..c6a6f64a3 100644 --- a/control/tests/tf_input_element_int.py +++ b/control/tests/input_element_int_test.py @@ -1,10 +1,18 @@ -import unittest +# input_element_int_test.py +# +# Author: Kangwon Lee (kangwonlee) +# Date: 22 Oct 2017 +# +# Unit tests contributed as part of PR #158, "SISO tf() may not work +# with numpy arrays with numpy.int elements" +# +# Modified: +# * 29 Dec 2017, RMM - updated file name and added header +import unittest import numpy as np - import control as ctl - class TestTfInputIntElement(unittest.TestCase): # currently these do not pass def test_tf_den_with_numpy_int_element(self): From 68b5dd153d67388b448d96e25edb618aa2864c41 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Fri, 29 Dec 2017 14:56:42 -0800 Subject: [PATCH 17/21] address additional comments from roryyorke in PR #135 --- .gitignore | 6 +++++- control/tests/xferfcn_input_test.py | 4 ++-- control/xferfcn.py | 25 +++++++++++++------------ 3 files changed, 20 insertions(+), 15 deletions(-) diff --git a/.gitignore b/.gitignore index ff6a58c0a..0262ab46f 100644 --- a/.gitignore +++ b/.gitignore @@ -9,6 +9,7 @@ __conda_*.txt record.txt build.log *.egg-info/ +.eggs/ .coverage doc/_build doc/generated @@ -18,4 +19,7 @@ examples/.ipynb_checkpoints/ .project Untitled*.ipynb *.idea/ -.eggs + +# Files created by or for emacs (RMM, 29 Dec 2017) +*~ +TAGS diff --git a/control/tests/xferfcn_input_test.py b/control/tests/xferfcn_input_test.py index cb9da4427..acddc64ae 100644 --- a/control/tests/xferfcn_input_test.py +++ b/control/tests/xferfcn_input_test.py @@ -1,7 +1,7 @@ #!/usr/bin/env python # -# xferfcn_test.py - test TransferFunction class -# RMM, 30 Mar 2011 (based on TestXferFcn from v0.4a) +# xferfcn_input_test.py - test inputs to TransferFunction class +# jed-frey, 18 Feb 2017 (based on xferfcn_test.py) import unittest import numpy as np diff --git a/control/xferfcn.py b/control/xferfcn.py index ebe35621f..982b8c5ba 100644 --- a/control/xferfcn.py +++ b/control/xferfcn.py @@ -10,7 +10,6 @@ # Python 3 compatibility (needs to go here) from __future__ import print_function from __future__ import division -from __future__ import absolute_import """Copyright (c) 2010 by California Institute of Technology All rights reserved. @@ -170,13 +169,6 @@ def __init__(self, *args): if zeronum: den[i][j] = ones(1) - # Check for coefficients that are ints and convert to floats - # TODO - for k in range(den[i][j]): - if (isinstance(data[i], (int, np.int))): - den[i][j][k] = float(den[i][j][k]) - - LTI.__init__(self, inputs, outputs, dt) self.num = num self.den = den @@ -1338,7 +1330,7 @@ def _cleanPart(data): Returns ------- - data: correctly formatted transfer function part. + data: list of lists of ndarrays, with int converted to float ''' valid_types = (int, float, complex, np.number) valid_collection = (list, tuple, ndarray) @@ -1346,10 +1338,10 @@ def _cleanPart(data): if (isinstance(data, valid_types) or (isinstance(data, ndarray) and data.ndim == 0)): # Data is a scalar (including 0d ndarray) - return [[array([data])]] + data = [[array([data])]] elif (isinstance(data, valid_collection) and all([isinstance(d, valid_types) for d in data])): - return [[array(data]] + data = [[array(data)]] elif (isinstance(data, (list, tuple)) and isinstance(data[0], (list, tuple)) and (isinstance(data[0][0], valid_collection) and @@ -1359,10 +1351,19 @@ def _cleanPart(data): data[j] = list(data[j]) for k in range(len(data[j])): data[j][k] = array(data[j][k]) - return data else: # If the user passed in anything else, then it's unclear what # the meaning is. raise TypeError("The numerator and denominator inputs must be \ scalars or vectors (for\nSISO), or lists of lists of vectors (for SISO or \ MIMO).") + + # Check for coefficients that are ints and convert to floats + for i in range(len(data)): + for j in range(len(data[i])): + for k in range(len(data[i][j])): + if (isinstance(data[i][j][k], (int, np.int))): + data[i][j][k] = float(data[i][j][k]) + + return data + From 314c4eb75a7e202a73d7b21feb5ba64e28800002 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Sat, 30 Dec 2017 13:34:25 -0800 Subject: [PATCH 18/21] DOC: fixed typos + address @roryyorke feedback --- doc/conventions.rst | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/doc/conventions.rst b/doc/conventions.rst index 77d6a30dc..9c3b46ce5 100644 --- a/doc/conventions.rst +++ b/doc/conventions.rst @@ -15,7 +15,7 @@ LTI system representation Linear time invariant (LTI) systems are represented in python-control in state space, transfer function, or frequency response data (FRD) form. Most functions in the toolbox will operate on any of these data types and -functions for convering between between compatible types is provided. +functions for converting between between compatible types is provided. State space systems ------------------- @@ -44,11 +44,12 @@ transfer functions .. math:: - G(s) = \frac{\text{num}(s)}{\text{den(s)}} + G(s) = \frac{\text{num}(s)}{\text{den}(s)} = \frac{a_0 s^n + a_1 s^{n-1} + \cdots + a_n} {b_0 s^m + b_1 s^{m-1} + \cdots + b_m}, -where n is generally greater than m (for a proper transfer function). +where n is generally greater than or equal to m (for a proper transfer +function). To create a transfer function, use the :class:`TransferFunction` constructor: @@ -97,7 +98,7 @@ Conversion between representations ---------------------------------- LTI systems can be converted between representations either by calling the constructor for the desired data type using the original system as the sole -argument or using the explict conversion functions :func:`ss2tf` and +argument or using the explicit conversion functions :func:`ss2tf` and :func:`tf2ss`. .. _time-series-convention: From aa6e0dfc05fde5e70d48f1add3f89e01adb1e8d5 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Sat, 30 Dec 2017 13:35:28 -0800 Subject: [PATCH 19/21] replace np.pi with math.pi (avoids Mock() issues) + docstring update in margin --- control/ctrlutil.py | 9 ++------- control/freqplot.py | 23 ++++++++++++----------- control/margins.py | 9 +++++---- control/statesp.py | 7 ++++--- 4 files changed, 23 insertions(+), 25 deletions(-) diff --git a/control/ctrlutil.py b/control/ctrlutil.py index 258e5755a..35035c7e9 100644 --- a/control/ctrlutil.py +++ b/control/ctrlutil.py @@ -43,17 +43,12 @@ # Packages that we need access to from . import lti import numpy as np -from numpy import pi - -# Hack for sphinx.ext.autodoc: if numpy is a mock import, then numpy.pi -# will be assigned to _Mock() and this generates a type error -if not isinstance(pi, float): - pi = 3.14 +import math __all__ = ['unwrap', 'issys', 'db2mag', 'mag2db'] # Utility function to unwrap an angle measurement -def unwrap(angle, period=2*pi): +def unwrap(angle, period=2*math.pi): """Unwrap a phase angle to give a continuous curve Parameters diff --git a/control/freqplot.py b/control/freqplot.py index 315a9635b..c53e83f31 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -44,6 +44,7 @@ import matplotlib.pyplot as plt import scipy as sp import numpy as np +import math from .ctrlutil import unwrap from .bdalg import feedback @@ -128,7 +129,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, else: omega_limits = np.array(omega_limits) if Hz: - omega_limits *= 2.*np.pi + omega_limits *= 2.*math.pi if omega_num: omega = sp.logspace(np.log10(omega_limits[0]), np.log10(omega_limits[1]), num=omega_num, endpoint=True) else: @@ -142,7 +143,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, else: omega_sys = np.array(omega) if sys.isdtime(True): - nyquistfrq = 2. * np.pi * 1. / sys.dt / 2. + nyquistfrq = 2. * math.pi * 1. / sys.dt / 2. omega_sys = omega_sys[omega_sys < nyquistfrq] # TODO: What distance to the Nyquist frequency is appropriate? else: @@ -154,9 +155,9 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, phase = unwrap(phase) nyquistfrq_plot = None if Hz: - omega_plot = omega_sys / (2. * np.pi) + omega_plot = omega_sys / (2. * math.pi) if nyquistfrq: - nyquistfrq_plot = nyquistfrq / (2. * np.pi) + nyquistfrq_plot = nyquistfrq / (2. * math.pi) else: omega_plot = omega_sys if nyquistfrq: @@ -187,7 +188,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, # Phase plot ax_phase = plt.subplot(212, sharex=ax_mag); if deg: - phase_plot = phase * 180. / np.pi + phase_plot = phase * 180. / math.pi else: phase_plot = phase ax_phase.semilogx(omega_plot, phase_plot, *args, **kwargs) @@ -208,8 +209,8 @@ def genZeroCenteredSeries(val_min, val_max, period): ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], 15.), minor=True) else: ylim = ax_phase.get_ylim() - ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], np.pi / 4.)) - ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], np.pi / 12.), minor=True) + ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], math.pi / 4.)) + ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], math.pi / 12.), minor=True) ax_phase.grid(True, which='both') # ax_mag.grid(which='minor', alpha=0.3) # ax_mag.grid(which='major', alpha=0.9) @@ -449,7 +450,7 @@ def default_frequency_range(syslist, Hz=None, number_of_samples=None, feature_pe features_ = features_[features_ != 0.0]; features = np.concatenate((features, features_)) elif sys.isdtime(strict=True): - fn = np.pi * 1. / sys.dt + fn = math.pi * 1. / sys.dt # TODO: What distance to the Nyquist frequency is appropriate? freq_interesting.append(fn * 0.9) @@ -475,12 +476,12 @@ def default_frequency_range(syslist, Hz=None, number_of_samples=None, feature_pe features = np.array([1.]); if Hz: - features /= 2.*np.pi + features /= 2.*math.pi features = np.log10(features) lsp_min = np.floor(np.min(features) - feature_periphery_decade) lsp_max = np.ceil(np.max(features) + feature_periphery_decade) - lsp_min += np.log10(2.*np.pi) - lsp_max += np.log10(2.*np.pi) + lsp_min += np.log10(2.*math.pi) + lsp_max += np.log10(2.*math.pi) else: features = np.log10(features) lsp_min = np.floor(np.min(features) - feature_periphery_decade) diff --git a/control/margins.py b/control/margins.py index 5b1ac9131..7f4f06c23 100644 --- a/control/margins.py +++ b/control/margins.py @@ -50,11 +50,12 @@ $Id$ """ +import math import numpy as np +import scipy as sp from . import xferfcn from .lti import issiso from . import frdata -import scipy as sp __all__ = ['stability_margins', 'phase_crossover_frequencies', 'margin'] @@ -140,7 +141,7 @@ def stability_margins(sysdata, returnall=False, epsw=0.0): sys = sysdata elif getattr(sysdata, '__iter__', False) and len(sysdata) == 3: mag, phase, omega = sysdata - sys = frdata.FRD(mag * np.exp(1j * phase * np.pi/180), + sys = frdata.FRD(mag * np.exp(1j * phase * math.pi/180), omega, smooth=True) else: sys = xferfcn._convertToTransferFunction(sysdata) @@ -336,13 +337,13 @@ def phase_crossover_frequencies(sys): def margin(*args): - """margin(sys) + """margin(sysdata) Calculate gain and phase margins and associated crossover frequencies Parameters ---------- - sysdata: LTI system or (mag, phase, omega) sequence + sysdata : LTI system or (mag, phase, omega) sequence sys : StateSpace or TransferFunction Linear SISO system mag, phase, omega : sequence of array_like diff --git a/control/statesp.py b/control/statesp.py index 70d5e4635..0bf0a37e8 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -51,9 +51,10 @@ $Id$ """ +import math import numpy as np from numpy import all, angle, any, array, asarray, concatenate, cos, delete, \ - dot, empty, exp, eye, matrix, ones, pi, poly, poly1d, roots, shape, sin, \ + dot, empty, exp, eye, matrix, ones, poly, poly1d, roots, shape, sin, \ zeros, squeeze from numpy.random import rand, randn from numpy.linalg import solve, eigvals, matrix_rank @@ -367,7 +368,7 @@ def evalfr(self, omega): if isdtime(self, strict=True): dt = timebase(self) s = exp(1.j * omega * dt) - if (omega * dt > pi): + if (omega * dt > math.pi): warnings.warn("evalfr: frequency evaluation above Nyquist frequency") else: s = omega * 1.j @@ -798,7 +799,7 @@ def _rss_generate(states, inputs, outputs, type): poles[i] = complex(-exp(randn()), 3. * exp(randn())) elif type == 'd': mag = rand() - phase = 2. * pi * rand() + phase = 2. * math.pi * rand() poles[i] = complex(mag * cos(phase), mag * sin(phase)) poles[i+1] = complex(poles[i].real, -poles[i].imag) From fbcdc2123eb12e3a4cd36e58649baf64899dd25f Mon Sep 17 00:00:00 2001 From: arnold Date: Mon, 1 Jan 2018 22:27:53 -0700 Subject: [PATCH 20/21] Use slycot's tb05ad for faster and more accurate frequency response evaluation of state-space systems. If slycot is unavailible, use the built in horner function (instead of converting to a transfer function, as was done before). --- control/statesp.py | 99 ++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 88 insertions(+), 11 deletions(-) diff --git a/control/statesp.py b/control/statesp.py index 13e662f03..4a7700136 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -378,23 +378,100 @@ def horner(self, s): self.B) + self.D return array(resp) + # Method for generating the frequency response of the system - def freqresp(self, omega): - """Evaluate the system's transfer func. at a list of ang. frequencies. + def freqresp(self, omega_s): + """Evaluate the system's transfer func. at a list of freqs, omega. mag, phase, omega = self.freqresp(omega) - reports the value of the magnitude, phase, and angular frequency of the - system's transfer function matrix evaluated at s = i * omega, where - omega is a list of angular frequencies, and is a sorted version of the - input omega. + Reports the frequency response of the system, + + G(j*omega) = mag*exp(j*phase) + + for continuous time. For discrete time systems, the response is + evaluated around the unit circle such that + + G(exp(j*omega*dt)) = mag*exp(j*phase). + + Inputs: + ------ + omega_s: A list of frequencies in radians/sec at which the system + should be evaluated. The list can be either a python list + or a numpy array and will be sorted before evaluation. + + Returns: + ------- + mag: The magnitude (absolute value, not dB or log10) of the system + frequency response. + + phase: The wrapped phase in radians of the system frequency + response. + + omega_s: The list of sorted frequencies at which the response + was evaluated. """ - # when evaluating at many frequencies, much faster to convert to - # transfer function first and then evaluate, than to solve an - # n-dimensional linear system at each frequency - tf = _convertToTransferFunction(self) - return tf.freqresp(omega) + + # In case omega is passed in as a list, rather than a proper array. + omega_s = np.asarray(omega_s) + + numFreqs = len(omega_s) + Gfrf = np.empty((self.outputs, self.inputs, numFreqs), + dtype=np.complex128) + + # Sort frequency and calculate complex frequencies on either imaginary + # axis (continuous time) or unit circle (discrete time). + omega_s.sort() + if isdtime(self, strict=True): + dt = timebase(self) + cmplx_freqs = exp(1.j * omega_s * dt) + if ((omega_s * dt).any() > pi): + warn_message = ("evalfr: frequency evaluation" + " above Nyquist frequency") + warnings.warn(warn_message) + else: + cmplx_freqs = omega_s * 1.j + + # Do the frequency response evaluation. Use TB05AD from Slycot + # if it's available, otherwise use the built-in horners function. + try: + from slycot import tb05ad + + n = np.shape(self.A)[0] + m = self.inputs + p = self.outputs + # The first call both evalates C(sI-A)^-1 B and also returns + # hessenberg transformed matrices at, bt, ct. + result = tb05ad(n, m, p, cmplx_freqs[0], self.A, + self.B, self.C, job='NG') + # When job='NG', result = (at, bt, ct, g_i, hinvb, info) + at = result[0] + bt = result[1] + ct = result[2] + + # TB05AD freqency evaluation does not include direct feedthrough. + Gfrf[:, :, 0] = result[3] + self.D + + # Now, iterate through the remaining frequencies using the + # transformed state matrices, at, bt, ct. + + # Start at the second frequency, already have the first. + for kk, cmplx_freqs_kk in enumerate(cmplx_freqs[1:numFreqs]): + result = tb05ad(n, m, p, cmplx_freqs_kk, at, + bt, ct, job='NH') + # When job='NH', result = (g_i, hinvb, info) + + # kk+1 because enumerate starts at kk = 0. + # but zero-th spot is already filled. + Gfrf[:, :, kk+1] = result[0] + self.D + + except ImportError: # Slycot unavailable. Fall back to horner. + for kk, cmplx_freqs_kk in enumerate(cmplx_freqs): + Gfrf[:, :, kk] = self.horner(cmplx_freqs_kk) + + # mag phase omega_s + return np.abs(Gfrf), np.angle(Gfrf), omega_s # Compute poles and zeros def pole(self): From 66c008822e3d6463d371bcfc63767695d73a556b Mon Sep 17 00:00:00 2001 From: arnold Date: Mon, 1 Jan 2018 22:27:53 -0700 Subject: [PATCH 21/21] Use slycot's tb05ad for faster and more accurate frequency response evaluation of state-space systems. If slycot is unavailable, use the built in horner function (instead of converting to a transfer function, as was done before). --- control/statesp.py | 99 ++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 88 insertions(+), 11 deletions(-) diff --git a/control/statesp.py b/control/statesp.py index 4abc4175f..33beaeea8 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -385,23 +385,100 @@ def horner(self, s): self.B) + self.D return array(resp) + # Method for generating the frequency response of the system - def freqresp(self, omega): - """Evaluate the system's transfer func. at a list of ang. frequencies. + def freqresp(self, omega_s): + """Evaluate the system's transfer func. at a list of freqs, omega. mag, phase, omega = self.freqresp(omega) - reports the value of the magnitude, phase, and angular frequency of the - system's transfer function matrix evaluated at s = i * omega, where - omega is a list of angular frequencies, and is a sorted version of the - input omega. + Reports the frequency response of the system, + + G(j*omega) = mag*exp(j*phase) + + for continuous time. For discrete time systems, the response is + evaluated around the unit circle such that + + G(exp(j*omega*dt)) = mag*exp(j*phase). + + Inputs: + ------ + omega_s: A list of frequencies in radians/sec at which the system + should be evaluated. The list can be either a python list + or a numpy array and will be sorted before evaluation. + + Returns: + ------- + mag: The magnitude (absolute value, not dB or log10) of the system + frequency response. + + phase: The wrapped phase in radians of the system frequency + response. + + omega_s: The list of sorted frequencies at which the response + was evaluated. """ - # when evaluating at many frequencies, much faster to convert to - # transfer function first and then evaluate, than to solve an - # n-dimensional linear system at each frequency - tf = _convertToTransferFunction(self) - return tf.freqresp(omega) + + # In case omega is passed in as a list, rather than a proper array. + omega_s = np.asarray(omega_s) + + numFreqs = len(omega_s) + Gfrf = np.empty((self.outputs, self.inputs, numFreqs), + dtype=np.complex128) + + # Sort frequency and calculate complex frequencies on either imaginary + # axis (continuous time) or unit circle (discrete time). + omega_s.sort() + if isdtime(self, strict=True): + dt = timebase(self) + cmplx_freqs = exp(1.j * omega_s * dt) + if ((omega_s * dt).any() > pi): + warn_message = ("evalfr: frequency evaluation" + " above Nyquist frequency") + warnings.warn(warn_message) + else: + cmplx_freqs = omega_s * 1.j + + # Do the frequency response evaluation. Use TB05AD from Slycot + # if it's available, otherwise use the built-in horners function. + try: + from slycot import tb05ad + + n = np.shape(self.A)[0] + m = self.inputs + p = self.outputs + # The first call both evalates C(sI-A)^-1 B and also returns + # hessenberg transformed matrices at, bt, ct. + result = tb05ad(n, m, p, cmplx_freqs[0], self.A, + self.B, self.C, job='NG') + # When job='NG', result = (at, bt, ct, g_i, hinvb, info) + at = result[0] + bt = result[1] + ct = result[2] + + # TB05AD freqency evaluation does not include direct feedthrough. + Gfrf[:, :, 0] = result[3] + self.D + + # Now, iterate through the remaining frequencies using the + # transformed state matrices, at, bt, ct. + + # Start at the second frequency, already have the first. + for kk, cmplx_freqs_kk in enumerate(cmplx_freqs[1:numFreqs]): + result = tb05ad(n, m, p, cmplx_freqs_kk, at, + bt, ct, job='NH') + # When job='NH', result = (g_i, hinvb, info) + + # kk+1 because enumerate starts at kk = 0. + # but zero-th spot is already filled. + Gfrf[:, :, kk+1] = result[0] + self.D + + except ImportError: # Slycot unavailable. Fall back to horner. + for kk, cmplx_freqs_kk in enumerate(cmplx_freqs): + Gfrf[:, :, kk] = self.horner(cmplx_freqs_kk) + + # mag phase omega_s + return np.abs(Gfrf), np.angle(Gfrf), omega_s # Compute poles and zeros def pole(self):