Menu
▾
▴
python-control-discuss
|
From: <mur...@us...> - 2012-11-03 17:38:06
|
Revision: 224
http://sourceforge.net/p/python-control/code/224
Author: murrayrm
Date: 2012-11-03 17:38:02 +0000 (Sat, 03 Nov 2012)
Log Message:
-----------
Modified ChangeLog and README file to test out post-commit hooks. If
things are set up correctly, this message should go to the
python-control-discuss mailing list.
Modified Paths:
--------------
README
trunk/ChangeLog
Modified: README
===================================================================
--- README 2012-11-03 05:28:38 UTC (rev 223)
+++ README 2012-11-03 17:38:02 UTC (rev 224)
@@ -2,10 +2,10 @@
RMM, 23 May 09
This directory contains the top level directory for the Python Control
-Systems Library (python-control). This library is still under
-development, but is intended to serve as a wrapper for standard
-control system algorithms in the python programming environment.
+Systems Library (python-control).
- trunk/ latest code; still under development
+ branches/ various branches of the code, as needed
tags/ releases and other stable copies
- branches/ various branches of the code, as needed
+ trunk/ latest code
+ web/ project website pages
+
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2012-11-03 05:28:38 UTC (rev 223)
+++ trunk/ChangeLog 2012-11-03 17:38:02 UTC (rev 224)
@@ -1,3 +1,7 @@
+2012-11-03 Richard Murray <murray@altura.local>
+
+ * ../README: updated readme file as a test of post-commit hook
+
2012-11-02 Richard Murray <murray@altura.local>
* doc/conf.py, setup.py: updated version number to 0.6c
|
|
From: <mur...@us...> - 2012-11-03 21:03:42
|
Revision: 225
http://sourceforge.net/p/python-control/code/225
Author: murrayrm
Date: 2012-11-03 21:03:39 +0000 (Sat, 03 Nov 2012)
Log Message:
-----------
Many small updates to make compatible with python 3.x:
* Updated unit tests to skip certain tests if slycot is not installed
* Fixed some problems with timebase (dt) to avoid comparing None to int
* Fixed exception syntax in a few more spots (E, T to E(T))
* Fixed print syntax in a few more spots (backward compatible with 2.6+)
* Updated FRD to include python 3.x fixes installed in other modules
* tests/test_*.py scripts not working in 3.x, but individual tests OK
* Current version works with nosetests in py2.7, py3.2 ( |
|
From: <mur...@us...> - 2012-11-03 22:51:30
|
Revision: 226
http://sourceforge.net/p/python-control/code/226
Author: murrayrm
Date: 2012-11-03 22:51:28 +0000 (Sat, 03 Nov 2012)
Log Message:
-----------
updated documentation for 0.6c release
Modified Paths:
--------------
trunk/ChangeLog
trunk/doc/analysis.rst
trunk/doc/bdalg_strings.rst
trunk/doc/class_strings.rst
trunk/doc/freqplot.rst
trunk/doc/index.rst
trunk/doc/intro.rst
trunk/doc/modsimp_strings.rst
trunk/doc/modules.rst
trunk/doc/synthesis.rst
trunk/doc/timeresp.rst
trunk/src/bdalg.py
trunk/src/lti.py
trunk/src/statesp.py
Added Paths:
-----------
trunk/doc/utilities.rst
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2012-11-03 21:03:39 UTC (rev 225)
+++ trunk/ChangeLog 2012-11-03 22:51:28 UTC (rev 226)
@@ -1,5 +1,18 @@
2012-11-03 Richard Murray <murray@altura.local>
+ * doc/modules.rst: updated documentation format so that items are
+ broken down by functionality and not internal module name. Still
+ not complete, but most functions should now show up in documentation.
+
+ * doc/index.rst: got rid of todos from documentation (not formatted
+ correctly and shouldn't show up in user documentation)
+
+ * src/lti.py: added missing docstrings
+
+ * doc/intro.rst: small tweaks to text, including adding a link to FBS
+
+2012-11-03 Richard Murray <murray@altura.local>
+
* src/rlocus.py (_RLSortRoots): convert output of range() to
explicit list for python 3 compatability
Modified: trunk/doc/analysis.rst
===================================================================
--- trunk/doc/analysis.rst 2012-11-03 21:03:39 UTC (rev 225)
+++ trunk/doc/analysis.rst 2012-11-03 22:51:28 UTC (rev 226)
@@ -1,2 +1,22 @@
Control System Analysis
-***********************
\ No newline at end of file
+***********************
+
+.. autofunction:: control.acker
+.. autofunction:: control.ctrb
+.. autofunction:: control.care
+.. autofunction:: control.dare
+.. autofunction:: control.dlyap
+.. autofunction:: control.dcgain
+.. autofunction:: control.evalfr
+.. autofunction:: control.gram
+.. autofunction:: control.lyap
+.. autofunction:: control.freqresp
+.. autofunction:: control.margin
+.. autofunction:: control.markov
+.. autofunction:: control.obsv
+.. autofunction:: control.phase_crossover_frequencies
+.. autofunction:: control.pole
+.. autofunction:: control.root_locus
+.. autofunction:: control.stability_margins
+.. autofunction:: control.zero
+
Modified: trunk/doc/bdalg_strings.rst
===================================================================
--- trunk/doc/bdalg_strings.rst 2012-11-03 21:03:39 UTC (rev 225)
+++ trunk/doc/bdalg_strings.rst 2012-11-03 22:51:28 UTC (rev 226)
@@ -3,5 +3,7 @@
.. toctree::
-.. automodule:: bdalg
- :members:
+.. autofunction:: control.feedback
+.. autofunction:: control.negate
+.. autofunction:: control.parallel
+.. autofunction:: control.series
Modified: trunk/doc/class_strings.rst
===================================================================
--- trunk/doc/class_strings.rst 2012-11-03 21:03:39 UTC (rev 225)
+++ trunk/doc/class_strings.rst 2012-11-03 22:51:28 UTC (rev 226)
@@ -1,6 +1,11 @@
Python-Control Classes
**********************
+LTI System Class
+================
+.. automodule:: lti
+ :members:
+
State Space Class
=================
.. automodule:: statesp
Modified: trunk/doc/freqplot.rst
===================================================================
--- trunk/doc/freqplot.rst 2012-11-03 21:03:39 UTC (rev 225)
+++ trunk/doc/freqplot.rst 2012-11-03 22:51:28 UTC (rev 226)
@@ -5,10 +5,10 @@
Plotting routines
=================
-.. autofunction:: freqplot.bode_plot
-.. autofunction:: freqplot.nyquist_plot
-.. autofunction:: freqplot.gangof4_plot
-.. autofunction:: nichols.nichols_plot
+.. autofunction:: control.bode_plot
+.. autofunction:: control.nyquist_plot
+.. autofunction:: control.gangof4_plot
+.. autofunction:: control.nichols_plot
Utility functions
=================
Modified: trunk/doc/index.rst
===================================================================
--- trunk/doc/index.rst 2012-11-03 21:03:39 UTC (rev 225)
+++ trunk/doc/index.rst 2012-11-03 22:51:28 UTC (rev 226)
@@ -17,17 +17,13 @@
:maxdepth: 2
intro
- class_strings
modules
+ class_strings
+ matlab_strings
examples
Indices and tables
==================
* :ref:`genindex`
-* :ref:`modindex`
* :ref:`search`
-
-----------------------------------
-
-.. todolist::
Modified: trunk/doc/intro.rst
===================================================================
--- trunk/doc/intro.rst 2012-11-03 21:03:39 UTC (rev 225)
+++ trunk/doc/intro.rst 2012-11-03 22:51:28 UTC (rev 226)
@@ -10,14 +10,14 @@
Overview of the Toolbox
-----------------------
-The python-control package is a set of python classes and functions
-that implement common operations for the analysis and design of
-feedback control systems. The initial goal is to implement all of the
-functionality required to work through the examples in the textbook
-Feedback Systems by Astrom and Murray. A MATLAB compatibility package
-(control.matlab) is available that provides many of the common
-functions corresponding to commands available in the MATLAB Control
-Systems Toolbox.
+The python-control package is a set of python classes and functions that
+implement common operations for the analysis and design of feedback control
+systems. The initial goal is to implement all of the functionality required
+to work through the examples in the textbook `Feedback Systems
+<http://www.cds.caltech.edu/~murray/FBSwiki>`_ by Astrom and Murray. A
+MATLAB compatibility package (control.matlab) is available that provides
+many of the common functions corresponding to commands available in the
+MATLAB Control Systems Toolbox.
In addition to the documentation here, there is a project wiki that
contains some additional information about how to use the package
@@ -37,10 +37,9 @@
* You must include commas in vectors. So [1 2 3] must be [1, 2, 3].
* Functions that return multiple arguments use tuples
-* Can't use braces for collections; use tuples instead
-* Transfer functions are only implemented for SISO systems (due to
- limitations in the underlying signals.lti class); use state space
- representations for MIMO systems.
+* You cannot use braces for collections; use tuples instead
+* Transfer functions are currently only implemented for SISO systems; use
+ state space representations for MIMO systems.
Getting Started
---------------
@@ -52,15 +51,16 @@
3. To see if things are working correctly, run ipython -pylab and run
the script 'examples/secord-matlab.py'. This should generate a
step response, Bode plot and Nyquist plot for a simple second order
- system.
+ system. (For more detailed tests, run nosetests in the main directory.)
4. To see the commands that are available, run the following commands in
ipython::
>>> import control
- >>> ?control.matlab
+ >>> ?control
5. If you want to have a MATLAB-like environment for running the control
toolbox, use::
>>> from control.matlab import *
+ >>> ?control.matlab
Modified: trunk/doc/modsimp_strings.rst
===================================================================
--- trunk/doc/modsimp_strings.rst 2012-11-03 21:03:39 UTC (rev 225)
+++ trunk/doc/modsimp_strings.rst 2012-11-03 22:51:28 UTC (rev 226)
@@ -1,5 +1,6 @@
Model Simplification Tools
**************************
-.. automodule:: modelsimp
- :members:
+.. autofunction:: control.balred
+.. autofunction:: control.hsvd
+.. autofunction:: control.modred
Modified: trunk/doc/modules.rst
===================================================================
--- trunk/doc/modules.rst 2012-11-03 21:03:39 UTC (rev 225)
+++ trunk/doc/modules.rst 2012-11-03 22:51:28 UTC (rev 226)
@@ -1,12 +1,13 @@
-Python-Control Modules
-**********************
+Python-Control Functions
+************************
.. toctree::
+ creation
bdalg_strings
analysis
freqplot
timeresp
synthesis
modsimp_strings
- matlab_strings
+ utilities
Modified: trunk/doc/synthesis.rst
===================================================================
--- trunk/doc/synthesis.rst 2012-11-03 21:03:39 UTC (rev 225)
+++ trunk/doc/synthesis.rst 2012-11-03 22:51:28 UTC (rev 226)
@@ -1,2 +1,5 @@
Control System Synthesis
-************************
\ No newline at end of file
+************************
+
+.. autofunction:: control.lqr
+.. autofunction:: control.place
Modified: trunk/doc/timeresp.rst
===================================================================
--- trunk/doc/timeresp.rst 2012-11-03 21:03:39 UTC (rev 225)
+++ trunk/doc/timeresp.rst 2012-11-03 22:51:28 UTC (rev 226)
@@ -1,8 +1,13 @@
Time Domain Simulation
**********************
-.. automodule:: timeresp
- :members:
+Time responses
+==============
+.. autofunction:: control.forced_response
+.. autofunction:: control.initial_response
+.. autofunction:: control.step_response
+Phase portraits
+===============
.. autofunction:: phaseplot.phase_plot
.. autofunction:: phaseplot.box_grid
Added: trunk/doc/utilities.rst
===================================================================
--- trunk/doc/utilities.rst (rev 0)
+++ trunk/doc/utilities.rst 2012-11-03 22:51:28 UTC (rev 226)
@@ -0,0 +1,4 @@
+Utility Functions
+*****************
+
+.. autofunction:: control.unwrap
Modified: trunk/src/bdalg.py
===================================================================
--- trunk/src/bdalg.py 2012-11-03 21:03:39 UTC (rev 225)
+++ trunk/src/bdalg.py 2012-11-03 22:51:28 UTC (rev 226)
@@ -204,7 +204,7 @@
parallel
Notes
- ----
+ -----
This function is a wrapper for the feedback function in the StateSpace and
TransferFunction classes. It calls TransferFunction.feedback if `sys1` is a
TransferFunction object, and StateSpace.feedback if `sys1` is a StateSpace
Modified: trunk/src/lti.py
===================================================================
--- trunk/src/lti.py 2012-11-03 21:03:39 UTC (rev 225)
+++ trunk/src/lti.py 2012-11-03 22:51:28 UTC (rev 226)
@@ -1,7 +1,7 @@
"""lti.py
-The lti module contains the Lti parent class to the child classes StateSpace and
-TransferFunction. It is designed for use in the python-control library.
+The Lti module contains the Lti parent class to the child classes StateSpace
+and TransferFunction. It is designed for use in the python-control library.
Routines in this module:
@@ -13,7 +13,6 @@
"""
class Lti:
-
"""Lti is a parent class to linear time invariant control (LTI) objects.
Lti is the parent to the StateSpace and TransferFunction child
@@ -90,7 +89,16 @@
# Check to see if two timebases are equal
def timebaseEqual(sys1, sys2):
- # TODO: add docstring
+ """Check to see if two systems have the same timebase
+
+ timebaseEqual(sys1, sys2)
+
+ returns True if the timebases for the two systems are compatible. By
+ default, systems with timebase 'None' are compatible with either
+ discrete or continuous timebase systems. If two systems have a discrete
+ timebase (dt > 0) then their timebases must be equal.
+ """
+
if (type(sys1.dt) == bool or type(sys2.dt) == bool):
# Make sure both are unspecified discrete timebases
return type(sys1.dt) == type(sys2.dt) and sys1.dt == sys2.dt
@@ -102,7 +110,17 @@
# Check to see if a system is a discrete time system
def isdtime(sys, strict=False):
- # TODO: add docstring
+ """
+ Check to see if a system is a discrete time system
+
+ Parameters
+ ----------
+ sys : LTI system
+ System to be checked
+ strict: bool (default = False)
+ If strict is True, make sure that timebase is not None
+ """
+
# Check to see if this is a constant
if isinstance(sys, (int, float, complex)):
# OK as long as strict checking is off
@@ -122,7 +140,17 @@
# Check to see if a system is a continuous time system
def isctime(sys, strict=False):
- # TODO: add docstring
+ """
+ Check to see if a system is a continuous time system
+
+ Parameters
+ ----------
+ sys : LTI system
+ System to be checked
+ strict: bool (default = False)
+ If strict is True, make sure that timebase is not None
+ """
+
# Check to see if this is a constant
if isinstance(sys, (int, float, complex)):
# OK as long as strict checking is off
Modified: trunk/src/statesp.py
===================================================================
--- trunk/src/statesp.py 2012-11-03 21:03:39 UTC (rev 225)
+++ trunk/src/statesp.py 2012-11-03 22:51:28 UTC (rev 226)
@@ -1,11 +1,10 @@
-# -*- coding: utf-8 -*-
-"""stateSpace.py
+"""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.
+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.
Routines in this module:
@@ -73,7 +72,6 @@
Revised: Kevin K. Chen, Dec 10
$Id$
-
"""
from numpy import all, angle, any, array, asarray, concatenate, cos, delete, \
@@ -88,7 +86,6 @@
from control.lti import Lti, timebaseEqual, isdtime
class StateSpace(Lti):
-
"""The StateSpace class represents state space instances and functions.
The StateSpace class is used throughout the python-control library to
|
|
From: <mur...@us...> - 2012-11-03 23:02:25
|
Revision: 229
http://sourceforge.net/p/python-control/code/229
Author: murrayrm
Date: 2012-11-03 23:02:22 +0000 (Sat, 03 Nov 2012)
Log Message:
-----------
incrementing release count
Modified Paths:
--------------
trunk/doc/conf.py
trunk/setup.py
Modified: trunk/doc/conf.py
===================================================================
--- trunk/doc/conf.py 2012-11-03 23:01:31 UTC (rev 228)
+++ trunk/doc/conf.py 2012-11-03 23:02:22 UTC (rev 229)
@@ -60,7 +60,7 @@
# The short X.Y version.
version = '0.6'
# The full version, including alpha/beta/rc tags.
-release = '0.6c'
+release = '0.6d'
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.
Modified: trunk/setup.py
===================================================================
--- trunk/setup.py 2012-11-03 23:01:31 UTC (rev 228)
+++ trunk/setup.py 2012-11-03 23:02:22 UTC (rev 229)
@@ -4,7 +4,7 @@
from setuptools import setup
setup(name = 'control',
- version = '0.6c',
+ version = '0.6d',
description = 'Python Control Systems Library',
author = 'Richard Murray',
author_email = 'mu...@cd...',
|
|
From: <mur...@us...> - 2012-11-19 05:28:07
|
Revision: 247
http://sourceforge.net/p/python-control/code/247
Author: murrayrm
Date: 2012-11-19 05:28:05 +0000 (Mon, 19 Nov 2012)
Log Message:
-----------
Set up 'config' module that allows user-defined defaults
* Will eventually allow MATLAB defaults if user imports control.matlab first
New module 'canonical' for implementing canconical forms. Needed for
trajectory generation module (coming soon).
Modified Paths:
--------------
trunk/ChangeLog
trunk/src/__init__.py
trunk/src/freqplot.py
trunk/src/matlab.py
trunk/src/statesp.py
trunk/src/xferfcn.py
Added Paths:
-----------
trunk/src/canonical.py
trunk/src/config.py
Property Changed:
----------------
trunk/
trunk/src/
trunk/tests/
Property changes on: trunk
___________________________________________________________________
Modified: svn:ignore
- MANIFEST*
build
dist
+ MANIFEST*
build
dist
__pycache__
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2012-11-11 22:00:29 UTC (rev 246)
+++ trunk/ChangeLog 2012-11-19 05:28:05 UTC (rev 247)
@@ -1,3 +1,31 @@
+2012-11-10 Richard Murray <murray@altura.local>
+
+ * src/canonical.py: new module implementing conversions to selected
+ canonical forms: canonical_form(), reachable_form() implemented
+
+ * src/lti.py (issiso): new function to quickly check if a system is
+ single input, single output or not
+
+2012-11-04 Richard Murray <murray@altura.local>
+
+ * src/__init__.py: added action item to separate out matlab module
+ functionality (eventually)
+
+ * src/matlab.py: added code for setting up MATLAB defaults if matlab
+ submodule is imported instead of main (control) module. See
+ comments in __init.py__ for what needs to be done to implement.
+
+ * src/freqplot.py (bode_plot): moved default values of plotting
+ options to config module to allow easier use override
+
+ * src/config.py: new module to help set package defaults
+
+ * src/xferfcn.py (_convertToTransferFunction): removed extraneous
+ print statements if slycot is not installed
+
+ * src/statesp.py (StateSpace.__str__): fixed error in printing
+ timebase when dt is None (python 3)
+
---- control-0.6c released -----
2012-11-03 Richard Murray <murray@altura.local>
Property changes on: trunk/src
___________________________________________________________________
Modified: svn:ignore
- *.pyc
.*.swp
+ *.pyc
.*.swp
__pycache__
Modified: trunk/src/__init__.py
===================================================================
--- trunk/src/__init__.py 2012-11-11 22:00:29 UTC (rev 246)
+++ trunk/src/__init__.py 2012-11-19 05:28:05 UTC (rev 247)
@@ -63,7 +63,7 @@
from control.dtime import sample_system
from control.freqplot import bode_plot, nyquist_plot, gangof4_plot
from control.freqplot import bode, nyquist, gangof4
-from control.lti import timebase, timebaseEqual, isdtime, isctime
+from control.lti import issiso, timebase, timebaseEqual, isdtime, isctime
from control.margins import stability_margins, phase_crossover_frequencies
from control.mateqn import lyap, dlyap, care, dare
from control.modelsimp import hsvd, modred, balred, era, markov
@@ -77,9 +77,20 @@
from control.xferfcn import TransferFunction
from control.ctrlutil import unwrap, issys
from control.frdata import FRD
+from control.canonical import canonical_form, reachable_form
+# Exceptions
+from exception import *
+
# Import some of the more common (and benign) MATLAB shortcuts
# By default, don't import conflicting commands here
+#! TODO (RMM, 4 Nov 2012): remove MATLAB dependencies from __init__.py
+#!
+#! Eventually, all functionality should be in modules *other* than matlab.
+#! This will allow inclusion of the matlab module to set up a different set
+#! of defaults from the main package. At that point, the matlab module will
+#! allow provide compatibility with MATLAB but no package functionality.
+#!
from control.matlab import ss, tf, ss2tf, tf2ss, drss
from control.matlab import pole, zero, evalfr, freqresp, dcgain
from control.matlab import nichols, rlocus, margin
Added: trunk/src/canonical.py
===================================================================
--- trunk/src/canonical.py (rev 0)
+++ trunk/src/canonical.py 2012-11-19 05:28:05 UTC (rev 247)
@@ -0,0 +1,64 @@
+# canonical.py - functions for converting systems to canonical forms
+# RMM, 10 Nov 2012
+
+import control
+from numpy import zeros, shape, poly
+from numpy.linalg import inv
+
+def canonical_form(sys, form):
+ """Convert a system into canonical form
+
+ Parameters
+ ----------
+ xsys : StateSpace object
+ System to be transformed, with state 'x'
+ form : String
+ Canonical form for transformation. Chosen from:
+ * 'reachable' - reachable canonical form
+ * 'observable' - observable canonical form
+ * 'modal' - modal canonical form [not implemented]
+
+ Outputs
+ -------
+ zsys : StateSpace object
+ System in desired canonical form, with state 'z'
+ T : matrix
+ Coordinate transformation matrix, z = T*x
+ """
+
+ # Call the appropriate tranformation function
+ if form == 'reachable':
+ return reachable_form(xsys)
+ else:
+ raise control.ControlNotImplemented(
+ "Canonical form '%s' not yet implemented" % form)
+
+# Reachable canonical form
+def reachable_form(xsys):
+ # Check to make sure we have a SISO system
+ if not control.issiso(xsys):
+ raise control.ControlNotImplemented(
+ "Canonical forms for MIMO systems not yet supported")
+
+ # Create a new system, starting with a copy of the old one
+ zsys = control.StateSpace(xsys)
+
+ # Generate the system matrices for the desired canonical form
+ zsys.B = zeros(shape(xsys.B)); zsys.B[0, 0] = 1;
+ zsys.A = zeros(shape(xsys.A))
+ Apoly = poly(xsys.A) # characteristic polynomial
+ for i in range(0, xsys.states):
+ zsys.A[0, i] = -Apoly[i+1] / Apoly[0]
+ if (i+1 < xsys.states): zsys.A[i+1, i] = 1
+
+ # Compute the reachability matrices for each set of states
+ Wrx = control.ctrb(xsys.A, xsys.B)
+ Wrz = control.ctrb(zsys.A, zsys.B)
+
+ # Transformation from one form to another
+ Tzx = Wrz * inv(Wrx)
+
+ # Finally, compute the output matrix
+ zsys.C = xsys.C * inv(Tzx)
+
+ return zsys, Tzx
Added: trunk/src/config.py
===================================================================
--- trunk/src/config.py (rev 0)
+++ trunk/src/config.py 2012-11-19 05:28:05 UTC (rev 247)
@@ -0,0 +1,27 @@
+# config.py - package defaults
+# RMM, 4 Nov 2012
+#
+# This file contains default values and utility functions for setting
+# variables that control the behavior of the control package.
+# Eventually it will be possible to read and write configuration
+# files. For now, you can just choose between MATLAB and FBS default
+# values.
+
+# Bode plot defaults
+bode_dB = False # Bode plot magnitude units
+bode_deg = True # Bode Plot phase units
+bode_Hz = False # Bode plot frequency units
+
+# Set defaults to match MATLAB
+def use_matlab_defaults():
+ # Bode plot defaults
+ global bode_dB; bode_dB = True
+ global bode_deg; bode_deg = True
+ global bode_Hz; bode_Hz = True
+
+# Set defaults to match FBS (Astrom and Murray)
+def use_fbs_defaults():
+ # Bode plot defaults
+ global bode_dB; bode_dB = False
+ global bode_deg; bode_deg = True
+ global bode_Hz; bode_Hz = True
Modified: trunk/src/freqplot.py
===================================================================
--- trunk/src/freqplot.py 2012-11-11 22:00:29 UTC (rev 246)
+++ trunk/src/freqplot.py 2012-11-19 05:28:05 UTC (rev 247)
@@ -57,7 +57,7 @@
#
# Bode plot
-def bode_plot(syslist, omega=None, dB=False, Hz=False, deg=True,
+def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
Plot=True, *args, **kwargs):
"""Bode plot for a system
@@ -105,6 +105,12 @@
>>> sys = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.")
>>> mag, phase, omega = bode(sys)
"""
+ # Set default values for options
+ import control.config
+ if (dB is None): dB = control.config.bode_dB
+ if (deg is None): deg = control.config.bode_deg
+ if (Hz is None): Hz = control.config.bode_Hz
+
# If argument was a singleton, turn it into a list
if (not getattr(syslist, '__iter__', False)):
syslist = (syslist,)
Modified: trunk/src/matlab.py
===================================================================
--- trunk/src/matlab.py 2012-11-11 22:00:29 UTC (rev 246)
+++ trunk/src/matlab.py 2012-11-19 05:28:05 UTC (rev 247)
@@ -69,6 +69,18 @@
from scipy.signal import zpk2ss, ss2zpk, tf2zpk, zpk2tf
from numpy import linspace, logspace
+# If configuration is not yet set, import and use MATLAB defaults
+#! NOTE (RMM, 4 Nov 2012): MATLAB default initialization commented out for now
+#!
+#! This code will eventually be used so that import control.matlab will
+#! automatically use MATLAB defaults, while import control will use package
+#! defaults. In order for that to work, we need to make sure that
+#! __init__.py does not include anything in the MATLAB module.
+# import sys
+# if not ('control.config' in sys.modules):
+# import control.config
+# control.config.use_matlab()
+
# Control system library
import control.ctrlutil as ctrlutil
import control.freqplot as freqplot
Modified: trunk/src/statesp.py
===================================================================
--- trunk/src/statesp.py 2012-11-11 22:00:29 UTC (rev 246)
+++ trunk/src/statesp.py 2012-11-19 05:28:05 UTC (rev 247)
@@ -138,6 +138,7 @@
# Here we're going to convert inputs to matrices, if the user gave a
# non-matrix type.
+ #! TODO: [A, B, C, D] = map(matrix, [A, B, C, D])?
matrices = [A, B, C, D]
for i in range(len(matrices)):
# Convert to matrix first, if necessary.
@@ -215,9 +216,10 @@
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 (self.dt > 0):
+ elif (not (self.dt is None) and type(self.dt) != bool and self.dt > 0):
str += "\ndt = " + self.dt.__str__() + "\n"
return str
Modified: trunk/src/xferfcn.py
===================================================================
--- trunk/src/xferfcn.py 2012-11-11 22:00:29 UTC (rev 246)
+++ trunk/src/xferfcn.py 2012-11-19 05:28:05 UTC (rev 247)
@@ -266,6 +266,7 @@
mimo = self.inputs > 1 or self.outputs > 1
if (var == None):
+ #! TODO: replace with standard calls to lti functions
var = 's' if self.dt == None or self.dt == 0 else 'z'
outstr = ""
@@ -294,6 +295,7 @@
# See if this is a discrete time system with specific sampling time
if (not (self.dt is None) and type(self.dt) != bool and self.dt > 0):
+ #! TODO: replace with standard calls to lti functions
outstr += "\ndt = " + self.dt.__str__() + "\n"
return outstr
@@ -945,8 +947,8 @@
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)
+ # print(num)
+ # print(den)
return TransferFunction(num, den, sys.dt)
Property changes on: trunk/tests
___________________________________________________________________
Added: svn:ignore
+ __pycache__
|
|
From: <mur...@us...> - 2013-06-09 14:58:04
|
Revision: 274
http://sourceforge.net/p/python-control/code/274
Author: murrayrm
Date: 2013-06-09 14:57:59 +0000 (Sun, 09 Jun 2013)
Log Message:
-----------
merging changes from Rane van Paassen - see 2013-05-23 ChangeLog entry for details
Modified Paths:
--------------
trunk/ChangeLog
trunk/src/bdalg.py
trunk/src/frdata.py
trunk/src/lti.py
trunk/src/margins.py
trunk/src/matlab.py
trunk/src/modelsimp.py
trunk/src/rlocus.py
trunk/src/statesp.py
trunk/src/timeresp.py
trunk/src/xferfcn.py
trunk/tests/convert_test.py
trunk/tests/matlab_test.py
trunk/tests/slycot_convert_test.py
trunk/tests/statesp_test.py
trunk/tests/timeresp_test.py
trunk/tests/xferfcn_test.py
Removed Paths:
-------------
branches/add-frd/
branches/mix-matrices-and-sys/
Property Changed:
----------------
trunk/
trunk/src/ctrlutil.py
Index: trunk
===================================================================
--- trunk 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk 2013-06-09 14:57:59 UTC (rev 274)
Property changes on: trunk
___________________________________________________________________
Modified: svn:mergeinfo
## -1 +1,2 ##
/branches/add-frd:188-219
+/branches/mix-matrices-and-sys:249-273
\ No newline at end of property
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/ChangeLog 2013-06-09 14:57:59 UTC (rev 274)
@@ -1,3 +1,24 @@
+2013-05-23 Rene van Paassen <ren...@gm...>
+
+ * src/margin.py: re-vamped stability_margins function. Now
+ uses analytical method for TF and state-space, works on FRD
+ objects too, and can optionally return all margins found for
+ the analytical method.
+
+ * src/xferfcn.py and src/frdata.py: Now return array result
+ for fresp when called with an iterable frequency vector
+
+ * src/xferfcn.py and src/statesp.py: added minreal functions
+
+ * src/bdalg.py: implemented append and connect functions (old
+ version, with indexes, not with names)
+
+ * src/matlab.py: added connect, append, minreal
+
+ * src/xferfcn.py and src/statesp.py: added/corrected
+ __rmul__. In most cases a gain matrix (np.matrix) may now be
+ used in combination with ss or tf systems.
+
2012-11-10 Richard Murray <murray@altura.local>
* src/canonical.py: new module implementing conversions to selected
Modified: trunk/src/bdalg.py
===================================================================
--- trunk/src/bdalg.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/bdalg.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -4,10 +4,12 @@
Routines in this module:
+append
series
parallel
negate
feedback
+connect
"""
@@ -236,3 +238,97 @@
sys2 = tf._convertToTransferFunction(sys2)
return sys1.feedback(sys2, sign)
+
+def append(*sys):
+ '''
+ Group models by appending their inputs and outputs
+
+ Forms an augmented system model, and appends the inputs and
+ outputs together. The system type will be the type of the first
+ system given; if you mix state-space systems and gain matrices,
+ make sure the gain matrices are not first.
+
+ Parameters.
+ -----------
+ sys1, sys2, ... sysn: StateSpace or Transferfunction
+ LTI systems to combine
+
+
+ Returns
+ -------
+ sys: LTI system
+ Combined LTI system, with input/output vectors consisting of all
+ input/output vectors appended
+
+ Examples
+ --------
+ >>> sys1 = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.")
+ >>> sys2 = ss("-1.", "1.", "1.", "0.")
+ >>> sys = append(sys1, sys2)
+
+ .. todo::
+ also implement for transfer function, zpk, etc.
+ '''
+ s1 = sys[0]
+ for s in sys[1:]:
+ s1 = s1.append(s)
+ return s1
+
+def connect(sys, Q, inputv, outputv):
+ '''
+ Index-base interconnection of system
+
+ The system sys is a system typically constructed with append, with
+ multiple inputs and outputs. The inputs and outputs are connected
+ according to the interconnection matrix Q, and then the final
+ inputs and outputs are trimmed according to the inputs and outputs
+ listed in inputv and outputv.
+
+ Note: to have this work, inputs and outputs start counting at 1!!!!
+
+ Parameters.
+ -----------
+ sys: StateSpace Transferfunction
+ System to be connected
+ Q: 2d array
+ Interconnection matrix. First column gives the input to be connected
+ second column gives the output to be fed into this input. Negative
+ values for the second column mean the feedback is negative, 0 means
+ no connection is made
+ inputv: 1d array
+ list of final external inputs
+ outputv: 1d array
+ list of final external outputs
+
+ Returns
+ -------
+ sys: LTI system
+ Connected and trimmed LTI system
+
+ Examples
+ --------
+ >>> sys1 = ss("1. -2; 3. -4", "5.; 7", "6, 8", "9.")
+ >>> sys2 = ss("-1.", "1.", "1.", "0.")
+ >>> sys = append(sys1, sys2)
+ >>> Q = sp.mat([ [ 1, 2], [2, -1] ]) # basically feedback, output 2 in 1
+ >>> sysc = connect(sys, Q, [2], [1, 2])
+ '''
+ # first connect
+ K = sp.zeros( (sys.inputs, sys.outputs) )
+ for r in sp.array(Q).astype(int):
+ inp = r[0]-1
+ for outp in r[1:]:
+ if outp > 0 and outp <= sys.outputs:
+ K[inp,outp-1] = 1.
+ elif outp < 0 and -outp >= -sys.outputs:
+ K[inp,-outp-1] = -1.
+ sys = sys.feedback(sp.matrix(K), sign=1)
+
+ # now trim
+ Ytrim = sp.zeros( (len(outputv), sys.outputs) )
+ Utrim = sp.zeros( (sys.inputs, len(inputv)) )
+ for i,u in enumerate(inputv):
+ Utrim[u-1,i] = 1.
+ for i,y in enumerate(outputv):
+ Ytrim[i,y-1] = 1.
+ return sp.matrix(Ytrim)*sys*sp.matrix(Utrim)
Index: trunk/src/ctrlutil.py
===================================================================
--- trunk/src/ctrlutil.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/ctrlutil.py 2013-06-09 14:57:59 UTC (rev 274)
Property changes on: trunk/src/ctrlutil.py
___________________________________________________________________
Modified: svn:mergeinfo
## -1,2 +1,3 ##
/branches/add-frd/src/ctrlutil.py:188-219
+/branches/mix-matrices-and-sys/src/ctrlutil.py:249-273
/trunk/src/ctrlutil.py:140-142
\ No newline at end of property
Modified: trunk/src/frdata.py
===================================================================
--- trunk/src/frdata.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/frdata.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -107,7 +107,7 @@
def __init__(self, *args, **kwargs):
"""Construct a transfer function.
- The default constructor is FRD(w, d), where w is an iterable of
+ The default constructor is FRD(d, w), where w is an iterable of
frequency points, and d is the matching frequency data.
If d is a single list, 1d array, or tuple, a SISO system description
@@ -142,8 +142,8 @@
# The user provided a response and a freq vector
self.fresp = array(args[0], dtype=complex)
if len(self.fresp.shape) == 1:
- self.fresp.reshape(1, 1, len(args[0]))
- self.omega = array(args[1])
+ self.fresp = self.fresp.reshape(1, 1, len(args[0]))
+ self.omega = array(args[1], dtype=float)
if len(self.fresp.shape) != 3 or \
self.fresp.shape[-1] != self.omega.shape[-1] or \
len(self.omega.shape) != 1:
@@ -189,9 +189,9 @@
if mimo:
outstr.append("Input %i to output %i:" % (i + 1, j + 1))
outstr.append('Freq [rad/s] Response ')
- outstr.append('------------ ------------------------')
+ outstr.append('------------ ---------------------')
outstr.extend(
- [ '%12.3f %10.4g + %10.4g' % (w, m, p)
+ [ '%12.3f %10.4g%+10.4gj' % (w, m, p)
for m, p, w in zip(real(self.fresp[j,i,:]), imag(self.fresp[j,i,:]), wt) ])
@@ -340,7 +340,10 @@
"""
# Preallocate the output.
- out = empty((self.outputs, self.inputs), dtype=complex)
+ if getattr(omega, '__iter__', False):
+ out = empty((self.outputs, self.inputs, len(omega)), dtype=complex)
+ else:
+ out = empty((self.outputs, self.inputs), dtype=complex)
if self.ifunc is None:
try:
@@ -350,11 +353,18 @@
"Frequency %f not in frequency list, try an interpolating"
" FRD if you want additional points")
else:
- for i in range(self.outputs):
- for j in range(self.inputs):
- frraw = splev(omega, self.ifunc[i,j], der=0)
- out[i,j] = frraw[0] + 1.0j*frraw[1]
-
+ if getattr(omega, '__iter__', False):
+ for i in range(self.outputs):
+ for j in range(self.inputs):
+ for k,w in enumerate(omega):
+ frraw = splev(w, self.ifunc[i,j], der=0)
+ out[i,j,k] = frraw[0] + 1.0j*frraw[1]
+ else:
+ for i in range(self.outputs):
+ for j in range(self.inputs):
+ frraw = splev(omega, self.ifunc[i,j], der=0)
+ out[i,j] = frraw[0] + 1.0j*frraw[1]
+
return out
# Method for generating the frequency response of the system
Modified: trunk/src/lti.py
===================================================================
--- trunk/src/lti.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/lti.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -12,6 +12,8 @@
timebaseEqual()
"""
+from numpy import absolute, real
+
class Lti:
"""Lti is a parent class to linear time invariant control (LTI) objects.
@@ -64,6 +66,12 @@
self.outputs = outputs
self.dt = dt
+ def damp(self):
+ poles = self.pole()
+ wn = absolute(poles)
+ Z = -real(poles)/wn
+ return wn, Z, poles
+
# Test to see if a system is SISO
def issiso(sys, strict=False):
if isinstance(sys, (int, float, complex)) and not strict:
Modified: trunk/src/margins.py
===================================================================
--- trunk/src/margins.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/margins.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -53,15 +53,37 @@
import numpy as np
import control.xferfcn as xferfcn
from control.freqplot import bode
-from control.lti import isdtime
+from control.lti import isdtime, issiso
+import control.frdata as frdata
+import scipy as sp
-# gain and phase margins
-# contributed by Sawyer B. Fuller <mi...@ca...>
-#! TODO - need to add unit test functions
-def stability_margins(sysdata, deg=True):
+# helper functions for stability_margins
+def _polyimsplit(pol):
+ """split a polynomial with (iw) applied into a real and an
+ imaginary part with w applied"""
+ rpencil = np.zeros_like(pol)
+ ipencil = np.zeros_like(pol)
+ rpencil[-1::-4] = 1.
+ rpencil[-3::-4] = -1.
+ ipencil[-2::-4] = 1.
+ ipencil[-4::-4] = -1.
+ return pol * rpencil, pol*ipencil
+
+def _polysqr(pol):
+ """return a polynomial squared"""
+ return np.polymul(pol, pol)
+
+# Took the framework for the old function by
+# Sawyer B. Fuller <mi...@ca...>, removed a lot of the innards
+# and replaced with analytical polynomial functions for Lti systems.
+#
+# idea for the frequency data solution copied/adapted from
+# https://github.com/alchemyst/Skogestad-Python/blob/master/BODE.py
+# Rene van Paassen <ren...@gm...>
+def stability_margins(sysdata, deg=True, returnall=False):
"""Calculate gain, phase and stability margins and associated
crossover frequencies.
-
+
Usage
-----
gm, pm, sm, wg, wp, ws = stability_margins(sysdata, deg=True)
@@ -76,86 +98,117 @@
bode frequency response data
deg=True: boolean
If true, all input and output phases in degrees, else in radians
-
+ returnall=False: boolean
+ If true, return all margins found. Note that for frequency data or
+ FRD systems, only one margin is found and returned.
+
Returns
-------
- gm, pm, sm, wg, wp, ws: float
+ gm, pm, sm, wg, wp, ws: float or array_like
Gain margin gm, phase margin pm, stability margin sm, and
associated crossover
frequencies wg, wp, and ws of SISO open-loop. If more than
one crossover frequency is detected, returns the lowest corresponding
margin.
- """
- #TODO do this precisely without the effects of discretization of frequencies?
- #TODO assumes SISO
- #TODO unit tests, margin plot
+ When requesting all margins, the return values are array_like,
+ and all margins are returns for linear systems not equal to FRD
+ """
+
+ try:
+ if isinstance(sysdata, frdata.FRD):
+ sys = frdata.FRD(sysdata, smooth=True)
+ elif isinstance(sysdata, xferfcn.TransferFunction):
+ sys = sysdata
+ elif getattr(sysdata, '__iter__', False) and len(sysdata) == 3:
+ mag, phase, omega = sysdata
+ sys = frdata.FRD(mag*np.exp((1j/360.)*phase), omega, smooth=True)
+ else:
+ sys = xferfcn._convertToTransferFunction(sysdata)
+ except Exception, e:
+ print (e)
+ raise ValueError("Margin sysdata must be either a linear system or "
+ "a 3-sequence of mag, phase, omega.")
- if (not getattr(sysdata, '__iter__', False)):
- sys = sysdata
+ # calculate gain of system
+ if isinstance(sys, xferfcn.TransferFunction):
+
+ # check for siso
+ if not issiso(sys):
+ raise ValueError("Can only do margins for SISO system")
+
+ # real and imaginary part polynomials in omega:
+ rnum, inum = _polyimsplit(sys.num[0][0])
+ rden, iden = _polyimsplit(sys.den[0][0])
- # TODO: implement for discrete time systems
- if (isdtime(sys, strict=True)):
- raise NotImplementedError("Function not implemented in discrete time")
+ # test imaginary part of tf == 0, for phase crossover/gain margins
+ test_w_180 = np.polyadd(np.polymul(inum, rden), np.polymul(rnum, -iden))
+ w_180 = np.roots(test_w_180)
+ w_180 = np.real(w_180[(np.imag(w_180) == 0) * (w_180 > 0.)])
+ w_180.sort()
- mag, phase, omega = bode(sys, deg=deg, Plot=False)
- elif len(sysdata) == 3:
- # TODO: replace with FRD object type?
- mag, phase, omega = sysdata
- else:
- raise ValueError("Margin sysdata must be either a linear system or a 3-sequence of mag, phase, omega.")
-
- if deg:
- cycle = 360.
- crossover = 180.
+ # test magnitude is 1 for gain crossover/phase margins
+ test_wc = np.polysub(np.polyadd(_polysqr(rnum), _polysqr(inum)),
+ np.polyadd(_polysqr(rden), _polysqr(iden)))
+ wc = np.roots(test_wc)
+ wc = np.real(wc[(np.imag(wc) == 0) * (wc > 0.)])
+ wc.sort()
+
+ # stability margin was a bitch to elaborate, relies on magnitude to
+ # point -1, then take the derivative. Second derivative needs to be >0
+ # to have a minimum
+ test_wstabn = np.polyadd(_polysqr(rnum), _polysqr(inum))
+ test_wstabd = np.polyadd(_polysqr(np.polyadd(rnum,rden)),
+ _polysqr(np.polyadd(inum,iden)))
+ test_wstab = np.polysub(
+ np.polymul(np.polyder(test_wstabn),test_wstabd),
+ np.polymul(np.polyder(test_wstabd),test_wstabn))
+
+ # find the solutions
+ wstab = np.roots(test_wstab)
+
+ # and find the value of the 2nd derivative there, needs to be positive
+ wstabplus = np.polyval(np.polyder(test_wstab), wstab)
+ wstab = np.real(wstab[(np.imag(wstab) == 0) * (wstab > 0.) *
+ (np.abs(wstabplus) > 0.)])
+ wstab.sort()
+
else:
- cycle = 2 * np.pi
- crossover = np.pi
-
- wrapped_phase = -np.mod(phase, cycle)
-
- # phase margin from minimum phase among all gain crossovers
- neg_mag_crossings_i = np.nonzero(np.diff(mag < 1) > 0)[0]
- mag_crossings_p = wrapped_phase[neg_mag_crossings_i]
- if len(neg_mag_crossings_i) == 0:
- if mag[0] < 1: # gain always less than one
- wp = np.nan
- pm = np.inf
- else: # gain always greater than one
- print("margin: no magnitude crossings found")
- wp = np.nan
- pm = np.nan
- else:
- min_mag_crossing_i = neg_mag_crossings_i[np.argmin(mag_crossings_p)]
- wp = omega[min_mag_crossing_i]
- pm = crossover + phase[min_mag_crossing_i]
- if pm < 0:
- print("warning: system unstable: negative phase margin")
-
- # gain margin from minimum gain margin among all phase crossovers
- neg_phase_crossings_i = np.nonzero(np.diff(wrapped_phase < -crossover) > 0)[0]
- neg_phase_crossings_g = mag[neg_phase_crossings_i]
- if len(neg_phase_crossings_i) == 0:
- wg = np.nan
- gm = np.inf
- else:
- min_phase_crossing_i = neg_phase_crossings_i[
- np.argmax(neg_phase_crossings_g)]
- wg = omega[min_phase_crossing_i]
- gm = abs(1/mag[min_phase_crossing_i])
- if gm < 1:
- print("warning: system unstable: gain margin < 1")
+ # a bit coarse, have the interpolated frd evaluated again
+ def mod(w):
+ """to give the function to calculate |G(jw)| = 1"""
+ return [np.abs(sys.evalfr(w[0])[0][0]) - 1]
- # stability margin from minimum abs distance from -1 point
- if deg:
- phase_rad = phase * np.pi / 180.
+ def arg(w):
+ """function to calculate the phase angle at -180 deg"""
+ return [np.angle(sys.evalfr(w[0])[0][0]) + np.pi]
+
+ def dstab(w):
+ """function to calculate the distance from -1 point"""
+ return np.abs(sys.evalfr(w[0])[0][0] + 1.)
+
+ # how to calculate the frequency at which |G(jw)| = 1
+ wc = np.array([sp.optimize.fsolve(mod, sys.omega[0])])[0]
+ w_180 = np.array([sp.optimize.fsolve(arg, sys.omega[0])])[0]
+ wstab = np.real(
+ np.array([sp.optimize.fmin(dstab, sys.omega[0], disp=0)])[0])
+
+ # margins, as iterables, converted frdata and xferfcn calculations to
+ # vector for this
+ PM = np.angle(sys.evalfr(wc)[0][0], deg=True) + 180
+ GM = 1/(np.abs(sys.evalfr(w_180)[0][0]))
+ SM = np.abs(sys.evalfr(wstab)[0][0]+1)
+
+ if returnall:
+ return GM, PM, SM, wc, w_180, wstab
else:
- phase_rad = phase
- L = mag * np.exp(1j * phase_rad) # complex loop response to -1 pt
- min_Lplus1_i = np.argmin(np.abs(L + 1))
- sm = np.abs(L[min_Lplus1_i] + 1)
- ws = phase[min_Lplus1_i]
+ return (
+ (GM.shape[0] or None) and GM[0],
+ (PM.shape[0] or None) and PM[0],
+ (SM.shape[0] or None) and SM[0],
+ (wc.shape[0] or None) and wc[0],
+ (w_180.shape[0] or None) and w_180[0],
+ (wstab.shape[0] or None) and wstab[0])
- return gm, pm, sm, wg, wp, ws
# Contributed by Steffen Waldherr <wal...@is...>
#! TODO - need to add test functions
Modified: trunk/src/matlab.py
===================================================================
--- trunk/src/matlab.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/matlab.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -89,6 +89,7 @@
from control.statesp import StateSpace, _rss_generate, _convertToStateSpace
from control.xferfcn import TransferFunction, _convertToTransferFunction
from control.lti import Lti #base class of StateSpace, TransferFunction
+from control.frdata import FRD
from control.dtime import sample_system
from control.exception import ControlArgument
@@ -96,11 +97,11 @@
from control.ctrlutil import unwrap
from control.freqplot import nyquist, gangof4
from control.nichols import nichols
-from control.bdalg import series, parallel, negate, feedback
+from control.bdalg import series, parallel, negate, feedback, append, connect
from control.pzmap import pzmap
from control.statefbk import ctrb, obsv, gram, place, lqr
from control.delay import pade
-from control.modelsimp import hsvd, balred, modred
+from control.modelsimp import hsvd, balred, modred, minreal
from control.mateqn import lyap, dlyap, dare, care
__doc__ += r"""
@@ -125,7 +126,7 @@
\* :func:`ss` create state-space (SS) models
\ dss create descriptor state-space models
\ delayss create state-space models with delayed terms
-\ frd create frequency response data (FRD) models
+\* :func:`frd` create frequency response data (FRD) models
\ lti/exp create pure continuous-time delays (TF and
ZPK only)
\ filt specify digital filters
@@ -156,7 +157,7 @@
\* :func:`tf` conversion to transfer function
\ zpk conversion to zero/pole/gain
\* :func:`ss` conversion to state space
-\ frd conversion to frequency data
+\* :func:`frd` conversion to frequency data
\ c2d continuous to discrete conversion
\ d2c discrete to continuous conversion
\ d2d resample discrete-time model
@@ -174,7 +175,7 @@
----------------------------------------------------------------------------
== ========================== ============================================
-\ append group LTI models by appending inputs/outputs
+\* :func:`~bdalg.append` group LTI models by appending inputs/outputs
\* :func:`~bdalg.parallel` connect LTI models in parallel
(see also overloaded ``+``)
\* :func:`~bdalg.series` connect LTI models in series
@@ -200,7 +201,7 @@
\ lti/order model order (number of states)
\* :func:`~pzmap.pzmap` pole-zero map (TF only)
\ lti/iopzmap input/output pole-zero map
-\ damp natural frequency, damping of system poles
+\* :func:`damp` natural frequency, damping of system poles
\ esort sort continuous poles by real part
\ dsort sort discrete poles by magnitude
\ lti/stabsep stable/unstable decomposition
@@ -243,7 +244,7 @@
----------------------------------------------------------------------------
== ========================== ============================================
-\ minreal minimal realization; pole/zero cancellation
+\* :func:`~modelsimp.minreal` minimal realization; pole/zero cancellation
\ ss/sminreal structurally minimal realization
\* :func:`~modelsimp.hsvd` hankel singular values (state contributions)
\* :func:`~modelsimp.balred` reduced-order approximations of LTI models
@@ -569,7 +570,43 @@
else:
raise ValueError("Needs 1 or 2 arguments; received %i." % len(args))
+def frd(*args):
+ '''
+ Construct a Frequency Response Data model, or convert a system
+ frd models store the (measured) frequency response of a system.
+
+ This function can be called in different ways:
+
+ ``frd(response, freqs)``
+ Create an frd model with the given response data, in the form of
+ complex response vector, at matching frequency freqs [in rad/s]
+
+ ``frd(sys, freqs)``
+ Convert an Lti system into an frd model with data at frequencies
+ freqs.
+
+ Parameters
+ ----------
+ response: array_like, or list
+ complex vector with the system response
+ freq: array_lik or lis
+ vector with frequencies
+ sys: Lti (StateSpace or TransferFunction)
+ A linear system
+
+ Returns
+ -------
+ sys: FRD
+ New frequency response system
+
+ See Also
+ --------
+ ss, tf
+ '''
+ return FRD(*args)
+
+
def ss2tf(*args):
"""
Transform a state space system to a transfer function.
@@ -714,7 +751,7 @@
else:
raise ValueError("Needs 1 or 2 arguments; received %i." % len(args))
-def rss(states=1, inputs=1, outputs=1):
+def rss(states=1, outputs=1, inputs=1):
"""
Create a stable **continuous** random state space object.
@@ -751,7 +788,7 @@
return _rss_generate(states, inputs, outputs, 'c')
-def drss(states=1, inputs=1, outputs=1):
+def drss(states=1, outputs=1, inputs=1):
"""
Create a stable **discrete** random state space object.
@@ -852,16 +889,20 @@
return sys.zero()
-def evalfr(sys, omega):
+def evalfr(sys, x):
"""
- Evaluate the transfer function of an LTI system at an angular frequency.
+ Evaluate the transfer function of an LTI system for a single complex
+ number x.
+
+ To evaluate at a frequency, enter x = omega*j, where omega is the
+ frequency in radians
Parameters
----------
sys: StateSpace or TransferFunction
Linear system
- omega: scalar
- Frequency
+ x: scalar
+ Complex number
Returns
-------
@@ -880,15 +921,17 @@
Examples
--------
>>> sys = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.")
- >>> evalfr(sys, 1.)
+ >>> evalfr(sys, 1j)
array([[ 44.8-21.4j]])
>>> # This is the transfer function matrix evaluated at s = i.
.. todo:: Add example with MIMO system
"""
+ if issiso(sys):
+ return sys.horner(x)[0][0]
+ return sys.horner(x)
+
- return sys.evalfr(omega)
-
def freqresp(sys, omega):
"""
Frequency response of an LTI system at multiple angular frequencies.
@@ -1057,10 +1100,11 @@
Returns
-------
- gm, pm, wg, wp : float
- Gain margin gm, phase margin pm (in deg), and associated crossover
- frequencies wg and wp (in rad/sec) of SISO open-loop. If more than
- one crossover frequency is detected, returns the lowest
+ gm, pm, Wcg, Wcp : float
+ Gain margin gm, phase margin pm (in deg), gain crossover frequency
+ (corresponding to phase margin) and phase crossover frequency
+ (corresponding to gain margin), in rad/sec of SISO open-loop.
+ If more than one crossover frequency is detected, returns the lowest
corresponding margin.
Examples
@@ -1083,7 +1127,7 @@
raise ValueError("Margin needs 1 or 3 arguments; received %i."
% len(args))
- return margin[0], margin[1], margin[3], margin[4]
+ return margin[0], margin[1], margin[4], margin[3]
def dcgain(*args):
'''
@@ -1139,10 +1183,48 @@
gain = sys.D - sys.C * sys.A.I * sys.B
return gain
+def damp(sys, doprint=True):
+ '''
+ Compute natural frequency, damping and poles of a system
+
+ The function takes 1 or 2 parameters
+
+ Parameters
+ ----------
+ sys: Lti (StateSpace or TransferFunction)
+ A linear system object
+ doprint:
+ if true, print table with values
+
+ Returns
+ -------
+ wn: array
+ Natural frequencies of the poles
+ damping: array
+ Damping values
+ poles: array
+ Pole locations
+
+ See Also
+ --------
+ pole
+ '''
+ wn, damping, poles = sys.damp()
+ if doprint:
+ print('_____Eigenvalue______ Damping___ Frequency_')
+ for p, d, w in zip(poles, damping, wn) :
+ if abs(p.imag) < 1e-12:
+ print("%10.4g %10.4g %10.4g" %
+ (p.real, 1.0, -p.real))
+ else:
+ print("%10.4g%+10.4gj %10.4g %10.4g" %
+ (p.real, p.imag, d, w))
+ return wn, damping, poles
+
# Simulation routines
# Call corresponding functions in timeresp, with arguments transposed
-def step(sys, T=None, X0=0., input=0, output=0, **keywords):
+def step(sys, T=None, X0=0., input=0, output=None, **keywords):
'''
Step response of a linear system
@@ -1192,7 +1274,7 @@
Examples
--------
- >>> T, yout = step(sys, T, X0)
+ >>> yout, T = step(sys, T, X0)
'''
T, yout = timeresp.step_response(sys, T, X0, input, output,
transpose = True, **keywords)
@@ -1405,3 +1487,4 @@
# TODO: add docstring
# Call the sample_system() function to do the work
return sample_system(sysc, Ts, method)
+
Modified: trunk/src/modelsimp.py
===================================================================
--- trunk/src/modelsimp.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/modelsimp.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -269,6 +269,33 @@
return rsys
+def minreal(sys, tol=None, verbose=True):
+ '''
+ Eliminates uncontrollable or unobservable states in state-space
+ models or cancelling pole-zero pairs in transfer functions. The
+ output sysr has minimal order and the same response
+ characteristics as the original model sys.
+
+ Parameters
+ ----------
+ sys: StateSpace or TransferFunction
+ Original system
+ tol: real
+ Tolerance
+ verbose: bool
+ Print results if True
+
+ Returns
+ -------
+ rsys: StateSpace or TransferFunction
+ Cleaned model
+ '''
+ sysr = sys.minreal(tol)
+ if verbose:
+ print("{nstates} states have been removed from the model".format(
+ nstates=len(sys.pole()) - len(sysr.pole())))
+ return sysr
+
def era(YY, m, n, nin, nout, r):
"""
Calculate an ERA model of order `r` based on the impulse-response data `YY`.
Modified: trunk/src/rlocus.py
===================================================================
--- trunk/src/rlocus.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/rlocus.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -50,9 +50,11 @@
import scipy.signal # signal processing toolbox
import pylab # plotting routines
import control.xferfcn as xferfcn
+from functools import partial
# Main function: compute a root locus diagram
-def root_locus(sys, kvect, xlim=None, ylim=None, plotstr='-', Plot=True):
+def root_locus(sys, kvect, xlim=None, ylim=None, plotstr='-', Plot=True,
+ PrintGain=True):
"""Calculate the root locus by finding the roots of 1+k*TF(s)
where TF is self.num(s)/self.den(s) and each k is an element
of kvect.
@@ -65,7 +67,9 @@
List of gains to use in computing diagram
Plot : boolean (default = True)
If True, plot magnitude and phase
-
+ PrintGain: boolean (default = True)
+ If True, report mouse clicks when close to the root-locus branches,
+ calculate gain, damping and print
Return values
-------------
rlist : list of computed root locations
@@ -80,6 +84,10 @@
# Create the plot
if (Plot):
+ f = pylab.figure()
+ if PrintGain:
+ cid = f.canvas.mpl_connect(
+ 'button_release_event', partial(_RLFeedbackClicks, sys=sys))
ax = pylab.axes();
# plot open loop poles
@@ -165,3 +173,12 @@
sorted[n,ind] = elem
prevrow = sorted[n,:]
return sorted
+
+def _RLFeedbackClicks(event, sys):
+ """Print root-locus gain feedback for clicks on the root-locus plot
+ """
+ s = complex(event.xdata, event.ydata)
+ K = -1./sys.horner(s)
+ if abs(K.real) > 1e-8 and abs(K.imag/K.real) < 0.04:
+ print("Clicked at %10.4g%+10.4gj gain %10.4g damp %10.4g" %
+ (s.real, s.imag, K.real, -1*s.real/abs(s)))
Modified: trunk/src/statesp.py
===================================================================
--- trunk/src/statesp.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/statesp.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -27,6 +27,7 @@
StateSpace.zero
StateSpace.feedback
StateSpace.returnScipySignalLti
+StateSpace.append
_convertToStateSpace
_rss_generate
@@ -338,10 +339,23 @@
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
- else:
- raise TypeError("can't interconnect systems")
+ # 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, 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."""
@@ -370,10 +384,16 @@
else:
s = omega * 1.j
- fresp = self.C * solve(s * eye(self.states) - self.A,
- self.B) + self.D
+ return self.horner(s)
- return array(fresp)
+ 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
# TODO: add discrete time check
@@ -475,6 +495,17 @@
return StateSpace(A, B, C, D, dt)
+ def minreal(self, tol=None):
+ """Calculate a minimal realization, removes unobservable and
+ uncontrollable states"""
+ try:
+ from slycot import tb01pd
+ A, B, C, nr = tb01pd(self.states, self.inputs, self.outputs,
+ self.A, self.B, self.C, tol=tol)
+ return StateSpace(A[:nr,:nr], B[:nr,:], C[:,:nr], self.D)
+ except ImportError:
+ raise TypeError("minreal requires slycot tb01pd")
+
# TODO: add discrete time check
def returnScipySignalLti(self):
"""Return a list of a list of scipy.signal.lti objects.
@@ -497,6 +528,33 @@
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)
+
# TODO: add discrete time check
def _convertToStateSpace(sys, **kw):
"""Convert a system to state space form (if needed).
@@ -757,3 +815,51 @@
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
Modified: trunk/src/timeresp.py
===================================================================
--- trunk/src/timeresp.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/timeresp.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -121,7 +121,7 @@
from copy import deepcopy
import warnings
from control.lti import Lti # base class of StateSpace, TransferFunction
-from control. statesp import StateSpace, _rss_generate, _convertToStateSpace, _mimo2siso
+from control. statesp import StateSpace, _rss_generate, _convertToStateSpace, _mimo2simo, _mimo2siso
from control.lti import isdtime, isctime
# Helper function for checking array-like parameters
@@ -386,15 +386,15 @@
return T, yout, xout
-def step_response(sys, T=None, X0=0., input=0, output=0, \
- transpose = False, **keywords):
+def step_response(sys, T=None, X0=0., input=0, output=None,
+ transpose = False, **keywords):
#pylint: disable=W0622
"""Step response of a linear system
- If the system has multiple inputs or outputs (MIMO), one input and one
- output have to be selected for the simulation. The parameters `input`
- and `output` do this. All other inputs are set to 0, all other outputs
- are ignored.
+ If the system has multiple inputs or outputs (MIMO), one input has
+ to be selected for the simulation. Optionally, one output may be
+ selected. The parameters `input` and `output` do this. All other
+ inputs are set to 0, all other outputs are ignored.
For information on the **shape** of parameters `T`, `X0` and
return values `T`, `yout` see: :ref:`time-series-convention`
@@ -416,7 +416,8 @@
Index of the input that will be used in this simulation.
output: int
- Index of the output that will be used in this simulation.
+ Index of the output that will be used in this simulation. Set to None
+ to not trim outputs
transpose: bool
If True, transpose all input and output arrays (for backward
@@ -447,7 +448,10 @@
>>> T, yout = step_response(sys, T, X0)
"""
sys = _convertToStateSpace(sys)
- sys = _mimo2siso(sys, input, output, warn_conversion=True)
+ if output == None:
+ sys = _mimo2simo(sys, input, warn_conversion=True)
+ else:
+ sys = _mimo2siso(sys, input, output, warn_conversion=True)
if T is None:
if isctime(sys):
T = _default_response_times(sys.A, 100)
@@ -459,13 +463,13 @@
U = np.ones_like(T)
T, yout, _xout = forced_response(sys, T, U, X0,
- transpose=transpose, **keywords)
+ transpose=transpose, **keywords)
return T, yout
-def initial_response(sys, T=None, X0=0., input=0, output=0, transpose=False,
- **keywords):
+def initial_response(sys, T=None, X0=0., input=0, output=None, transpose=False,
+ **keywords):
#pylint: disable=W0622
"""Initial condition response of a linear system
@@ -494,7 +498,8 @@
Index of the input that will be used in this simulation.
output: int
- Index of the output that will be used in this simulation.
+ Index of the output that will be used in this simulation. Set to None
+ to not trim outputs
transpose: bool
If True, transpose all input and output arrays (for backward
@@ -525,7 +530,10 @@
>>> T, yout = initial_response(sys, T, X0)
"""
sys = _convertToStateSpace(sys)
- sys = _mimo2siso(sys, input, output, warn_conversion=True)
+ if output == None:
+ sys = _mimo2simo(sys, input, warn_conversion=False)
+ else:
+ sys = _mimo2siso(sys, input, output, warn_conversion=False)
# Create time and input vectors; checking is done in forced_response(...)
# The initial vector X0 is created in forced_response(...) if necessary
@@ -534,11 +542,11 @@
U = np.zeros_like(T)
T, yout, _xout = forced_response(sys, T, U, X0, transpose=transpose,
- **keywords)
+ **keywords)
return T, yout
-def impulse_response(sys, T=None, X0=0., input=0, output=0,
+def impulse_response(sys, T=None, X0=0., input=0, output=None,
transpose=False, **keywords):
#pylint: disable=W0622
"""Impulse response of a linear system
@@ -568,7 +576,8 @@
Index of the input that will be used in this simulation.
output: int
- Index of the output that will be used in this simulation.
+ Index of the output that will be used in this simulation. Set to None
+ to not trim outputs
transpose: bool
If True, transpose all input and output arrays (for backward
@@ -599,7 +608,10 @@
>>> T, yout = impulse_response(sys, T, X0)
"""
sys = _convertToStateSpace(sys)
- sys = _mimo2siso(sys, input, output, warn_conversion=True)
+ if output == None:
+ sys = _mimo2simo(sys, input, warn_conversion=True)
+ else:
+ sys = _mimo2siso(sys, input, output, warn_conversion=True)
# System has direct feedthrough, can't simulate impulse response numerically
if np.any(sys.D != 0) and isctime(sys):
Modified: trunk/src/xferfcn.py
===================================================================
--- trunk/src/xferfcn.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/xferfcn.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -404,7 +404,8 @@
for k in range(self.inputs): # Multiply & add.
num_summand[k] = polymul(self.num[i][k], other.num[k][j])
den_summand[k] = polymul(self.den[i][k], other.den[k][j])
- num[i][j], den[i][j] = _addSISO(num[i][j], den[i][j],
+ num[i][j], den[i][j] = _addSISO(
+ num[i][j], den[i][j],
num_summand[k], den_summand[k])
return TransferFunction(num, den, dt)
@@ -412,8 +413,48 @@
def __rmul__(self, other):
"""Right multiply two LTI objects (serial connection)."""
- return self * other
+ # Convert the second argument to a transfer function.
+ if isinstance(other, (int, float, complex)):
+ other = _convertToTransferFunction(other, inputs=self.inputs,
+ outputs=self.inputs)
+ else:
+ other = _convertToTransferFunction(other)
+
+ # Check that the input-output sizes are consistent.
+ if other.inputs != self.outputs:
+ raise ValueError("C = A * B: A has %i column(s) (input(s)), but B \
+has %i row(s)\n(output(s))." % (other.inputs, self.outputs))
+ inputs = self.inputs
+ outputs = 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 (self.dt == other.dt):
+ dt = self.dt # use dt from first argument
+ else:
+ raise ValueError("Systems have different sampling times")
+
+ # Preallocate the numerator and denominator of the sum.
+ num = [[[0] for j in range(inputs)] for i in range(outputs)]
+ den = [[[1] for j in range(inputs)] for i in range(outputs)]
+
+ # Temporary storage for the summands needed to find the (i, j)th element
+ # of the product.
+ num_summand = [[] for k in range(other.inputs)]
+ den_summand = [[] for k in range(other.inputs)]
+
+ for i in range(outputs): # Iterate through rows of product.
+ for j in range(inputs): # Iterate through columns of product.
+ for k in range(other.inputs): # Multiply & add.
+ num_summand[k] = polymul(other.num[i][k], self.num[k][j])
+ den_summand[k] = polymul(other.den[i][k], self.den[k][j])
+ num[i][j], den[i][j] = _addSISO(num[i][j], den[i][j],
+ num_summand[k], den_summand[k])
+
+ return TransferFunction(num, den, dt)
+
# TODO: Division of MIMO transfer function objects is not written yet.
def __div__(self, other):
"""Divide two LTI objects."""
@@ -486,9 +527,20 @@
warn("evalfr: frequency evaluation above Nyquist frequency")
else:
s = 1.j * omega
+
+ 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.
+ '''
+
# Preallocate the output.
- out = empty((self.outputs, self.inputs), dtype=complex)
+ if getattr(s, '__iter__', False):
+ out = empty((self.outputs, self.inputs, len(s)), dtype=complex)
+ else:
+ out = empty((self.outputs, self.inputs), dtype=complex)
for i in range(self.outputs):
for j in range(self.inputs):
@@ -620,8 +672,12 @@
newzeros.append(z)
# keep result
- num[i][j] = gain * real(poly(newzeros))
+ if len(newzeros):
+ num[i][j] = gain * real(poly(newzeros))
+ else:
+ num[i][j] = array([gain])
den[i][j] = real(poly(poles))
+
# end result
return TransferFunction(num, den)
@@ -906,6 +962,15 @@
In the latter example, sys's matrix transfer function is [[1., 1., 1.]
[1., 1., 1.]].
+
+ If sys is an array-like type, then it is converted to a constant-gain
+ transfer function.
+
+ >>> sys = _convertToTransferFunction([[1. 0.], [2. 3.]])
+
+ In this example, the numerator matrix will be
+ [[[1.0], [0.0]], [[2.0], [3.0]]]
+ and the denominator matrix [[[1.0], [1.0]], [[1.0], [1.0]]]
"""
from control.statesp import StateSpace
@@ -966,5 +1031,16 @@
den = [[[1] for j in range(inputs)] for i in range(outputs)]
return TransferFunction(num, den)
- else:
- raise TypeError("Can't convert given type to TransferFunction system.")
+
+ # If this is array-like, try to create a constant feedthrough
+ try:
+ D = array(sys)
+ outputs, inputs = D.shape
+ num = [[[D[i,j]] for j in range(inputs)] for i in range(outputs)]
+ den = [[[1] for j in range(inputs)] for i in range(outputs)]
+ return TransferFunction(num, den)
+ except Exception, e:
+ print("Failure to assume argument is matrix-like in"
+ " _convertToTransferFunction, result %s" % e)
+
+ raise TypeError("Can't convert given type to TransferFunction system.")
Modified: trunk/tests/convert_test.py
===================================================================
--- trunk/tests/convert_test.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/tests/convert_test.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -57,7 +57,7 @@
for outputs in range(1, self.maxIO):
# start with a random SS system and transform to TF then
# back to SS, check that the matrices are the same.
- ssOriginal = matlab.rss(states, inputs, outputs)
+ ssOriginal = matlab.rss(states, outputs, inputs)
if (verbose):
self.printSys(ssOriginal, 1)
Modified: trunk/tests/matlab_test.py
===================================================================
--- trunk/tests/matlab_test.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/tests/matlab_test.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -11,9 +11,48 @@
from __future__ import print_function
import unittest
import numpy as np
+from scipy.linalg import eigvals
import scipy as sp
from control.matlab import *
+from control.frdata import FRD
+# for running these through Matlab or Octave
+'''
+siso_ss1 = ss([1. -2.; 3. -4.], [5.; 7.], [6. 8.], [0])
+
+siso_tf1 = tf([1], [1, 2, 1])
+siso_tf2 = tf([1, 1], [1, 2, 3, 1])
+
+siso_tf3 = tf(siso_ss1)
+siso_ss2 = ss(siso_tf2)
+siso_ss3 = ss(siso_tf3)
+siso_tf4 = tf(siso_ss2)
+
+A =[ 1. -2. 0. 0.;
+ 3. -4. 0. 0.;
+ 0. 0. 1. -2.;
+ 0. 0. 3. -4. ]
+B = [ 5. 0.;
+ 7. 0.;
+ 0. 5.;
+ 0. 7. ]
+C = [ 6. 8. 0. 0.;
+ 0. 0. 6. 8. ]
+D = [ 9. 0.;
+ 0. 9. ]
+mimo_ss1 = ss(A, B, C, D)
+
+% all boring, since no cross-over
+margin(siso_tf1)
+margin(siso_tf2)
+margin(siso_ss1)
+margin(siso_ss2)
+
+% make a bit better
+[gm, pm, gmc, pmc] = margin(siso_ss2*siso_ss2*2)
+
+'''
+
class TestMatlab(unittest.TestCase):
def setUp(self):
"""Set up some systems for testing out MATLAB functions"""
@@ -194,6 +233,9 @@
gm, pm, wg, wp = margin(self.siso_tf2);
gm, pm, wg, wp = margin(self.siso_ss1);
gm, pm, wg, wp = margin(self.siso_ss2);
+ gm, pm, wg, wp = margin(self.siso_ss2*self.siso_ss2*2);
+ np.testing.assert_array_almost_equal(
+ [gm, pm, wg, wp], [1.5451, 75.9933, 1.2720, 0.6559], decimal=3)
def testDcgain(self):
#Create different forms of a SISO system
@@ -274,13 +316,16 @@
freqresp(self.siso_tf3, w)
def testEvalfr(self):
- w = 1
- evalfr(self.siso_ss1, w)
+ w = 1j
+ self.assertEqual(evalfr(self.siso_ss1, w), 44.8-21.4j)
evalfr(self.siso_ss2, w)
evalfr(self.siso_ss3, w)
evalfr(self.siso_tf1, w)
evalfr(self.siso_tf2, w)
evalfr(self.siso_tf3, w)
+ np.testing.assert_array_almost_equal(
+ evalfr(self.mimo_ss1, w),
+ np.array( [[44.8-21.4j, 0.], [0., 44.8-21.4j]]))
def testHsvd(self):
hsvd(self.siso_ss1)
@@ -309,12 +354,12 @@
def testRss(self):
rss(1)
rss(2)
- rss(2, 3, 1)
+ rss(2, 1, 3)
def testDrss(self):
drss(1)
drss(2)
- drss(2, 3, 1)
+ drss(2, 1, 3)
def testCtrb(self):
ctrb(self.siso_ss1.A, self.siso_ss1.B)
@@ -372,6 +417,102 @@
for i in range(len(tfdata_1)):
np.testing.assert_array_almost_equal(tfdata_1[i], tfdata_2[i])
+ def testDamp(self):
+ A = np.mat('''-0.2 0.06 0 -1;
+ 0 0 1 0;
+ -17 0 -3.8 1;
+ 9.4 0 -0.4 -0.6''')
+ B = np.mat('''-0.01 0.06;
+ 0 0;
+ -32 5.4;
+ 2.6 -7''')
+ C = np.eye(4)
+ D = np.zeros((4,2))
+ sys = ss(A, B, C, D)
+ wn, Z, p = damp(sys, False)
+ print (wn)
+ np.testing.assert_array_almost_equal(
+ wn, np.array([4.07381994, 3.28874827, 3.28874827,
+ 1.08937685e-03]))
+ np.testing.assert_array_almost_equal(
+ Z, np.array([1.0, 0.07983139, 0.07983139, 1.0]))
+
+ def testConnect(self):
+ sys1 = ss("1. -2; 3. -4", "5.; 7", "6, 8", "9.")
+ sys2 = ss("-1.", "1.", "1.", "0.")
+ sys = append(sys1, sys2)
+ Q= np.mat([ [ 1, 2], [2, -1] ]) # basically feedback, output 2 in 1
+ sysc = connect(sys, Q, [2], [1, 2])
+ print(sysc)
+ np.testing.assert_array_almost_equal(
+ sysc.A, np.mat('1 -2 5; 3 -4 7; -6 -8 -10'))
+ np.testing.assert_array_almost_equal(
+ sysc.B, np.mat('0; 0; 1'))
+ np.testing.assert_array_almost_equal(
+ sysc.C, np.mat('6 8 9; 0 0 1'))
+ np.testing.assert_array_almost_equal(
+ sysc.D, np.mat('0; 0'))
+
+ def testConnect2(self):
+ sys = append(ss([[-5, -2.25], [4, 0]], [[2], [0]],
+ [[0, 1.125]], [[0]]),
+ ss([[-1.6667, 0], [1, 0]], [[2], [0]],
+ [[0, 3.3333]], [[0]]),
+ 1)
+ Q = [ [ 1, 3], [2, 1], [3, -2]]
+ sysc = connect(sys, Q, [3], [3, 1, 2])
+ np.testing.assert_array_almost_equal(
+ sysc.A, np.mat([[-5, -2.25, 0, -6.6666],
+ [4, 0, 0, 0],
+ [0, 2.25, -1.6667, 0],
+ [0, 0, 1, 0]]))
+ np.testing.assert_array_almost_equal(
+ sysc.B, np.mat([[2], [0], [0], [0]]))
+ np.testing.assert_array_almost_equal(
+ sysc.C, np.mat([[0, 0, 0, -3.3333],
+ [0, 1.125, 0, 0],
+ [0, 0, 0, 3.3333]]))
+ np.testing.assert_array_almost_equal(
+ sysc.D, np.mat([[1], [0], [0]]))
+
+
+
+ def testFRD(self):
+ h = tf([1], [1, 2, 2])
+ omega = np.logspace(-1, 2, 10)
+ frd1 = frd(h, omega)
+ assert isinstance(frd1, FRD)
+ frd2 = frd(frd1.fresp[0,0,:], omega)
+ assert isinstance(frd2, FRD)
+
+ def testMinreal(self, verbose=False):
+ """Test a minreal model reduction"""
+ #A = [-2, 0.5, 0; 0.5, -0.3, 0; 0, 0, -0.1]
+ A = [[-2, 0.5, 0], [0.5, -0.3, 0], [0, 0, -0.1]]
+ #B = [0.3, -1.3; 0.1, 0; 1, 0]
+ B = [[0.3, -1.3], [0.1, 0.], [1.0, 0.0]]
+ #C = [0, 0.1, 0; -0.3, -0.2, 0]
+ C = [[0., 0.1, 0.0], [-0.3, -0.2, 0.0]]
+ #D = [0 -0.8; -0.3 0]
+ D = [[0., -0.8], [-0.3, 0.]]
+ # sys = ss(A, B, C, D)
+
+ sys = ss(A, B, C, D)
+ sysr = minreal(sys)
+ self.assertEqual(sysr.states, 2)
+ self.assertEqual(sysr.inputs, sys.inputs)
+ self.assertEqual(sysr.outputs, sys.outputs)
+ np.testing.assert_array_almost_equal(
+ eigvals(sysr.A), [-2.136154, -0.1638459])
+
+ s = tf([1, 0], [1])
+ h = (s+1)*(s+2.00000000001)/(s+2)/(s**2+s+1)
+ hm = minreal(h)
+ 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])
+
+
#! TODO: not yet implemented
# def testMIMOtfdata(self):
# sisotf = ss2tf(self.siso_ss1)
Modified: trunk/tests/slycot_convert_test.py
===================================================================
--- trunk/tests/slycot_convert_test.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/tests/slycot_convert_test.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -37,7 +37,7 @@
for inputs in range(1, self.maxI+1):
for outputs in range(1, self.maxO+1):
for testNum in range(self.numTests):
- ssOriginal = matlab.rss(states, inputs, outputs)
+ ssOriginal = matlab.rss(states, outputs, inputs)
if (verbose):
print('====== Original SS ==========')
print(ssOriginal)
@@ -82,7 +82,7 @@
for testNum in range(self.numTests):
for inputs in range(1,1):
for outputs in range(1,1):
- ssOriginal = matlab.rss(states, inputs, outputs)
+ ssOriginal = matlab.rss(states, outputs, inputs)
tfOriginal_Actrb, tfOriginal_Bctrb, tfOriginal_Cctrb, tfOrigingal_nctrb, tfOriginal_index,\
tfOriginal_dcoeff, tfOriginal_ucoeff = tb04ad(states,inputs,outputs,\
Modified: trunk/tests/statesp_test.py
===================================================================
--- trunk/tests/statesp_test.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/tests/statesp_test.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -5,8 +5,10 @@
import unittest
import numpy as np
+from scipy.linalg import eigvals
import control.matlab as matlab
-from control.statesp import StateSpace
+from control.statesp import StateSpace, _convertToStateSpace
+from control.xferfcn import TransferFunction
class TestStateSpace(unittest.TestCase):
"""Tests for the StateSpace class."""
@@ -134,7 +136,74 @@
np.testing.assert_almost_equal(mag, truemag)
np.testing.assert_almost_equal(phase, truephase)
np.testing.assert_equal(omega, trueomega)
+
+ def testMinreal(self):
+ """Test a minreal model reduction"""
+ #A = [-2, 0.5, 0; 0.5, -0.3, 0; 0, 0, -0.1]
+ A = [[-2, 0.5, 0], [0.5, -0.3, 0], [0, 0, -0.1]]
+ #B = [0.3, -1.3; 0.1, 0; 1, 0]
+ B = [[0.3, -1.3], [0.1, 0.], [1.0, 0.0]]
+ #C = [0, 0.1, 0; -0.3, -0.2, 0]
+ C = [[0., 0.1, 0.0], [-0.3, -0.2, 0.0]]
+ #D = [0 -0.8; -0.3 0]
+ D = [[0., -0.8], [-0.3, 0.]]
+ # sys = ss(A, B, C, D)
+
+ sys = StateSpace(A, B, C, D)
+ sysr = sys.minreal()
+ self.assertEqual(sysr.states, 2)
+ self.assertEqual(sysr.inputs, sys.inputs)
+ self.assertEqual(sysr.outputs, sys.outputs)
+ np.testing.assert_array_almost_equal(
+ eigvals(sysr.A), [-2.136154, -0.1638459])
+
+ def testAppendSS(self):
+ """Test appending two state-space systems"""
+ A1 = [[-2, 0.5, 0], [0.5, -0.3, 0], [0, 0, -0.1]]
+ B1 = [[0.3, -1.3], [0.1, 0.], [1.0, 0.0]]
+ C1 = [[0., 0.1, 0.0], [-0.3, -0.2, 0.0]]
+ D1 = [[0., -0.8], [-0.3, 0.]]
+ A2 = [[-1.]]
+ B2 = [[1.2]]
+ C2 = [[0.5]]
+ D2 = [[0.4]]
+ A3 = [[-2, 0.5, 0, 0], [0.5, -0.3, 0, 0], [0, 0, -0.1, 0],
+ [0, 0, 0., -1.]]
+ B3 = [[0.3, -1.3, 0], [0.1, 0., 0], [1.0, 0.0, 0], [0., 0, 1.2]]
+ C3 = [[0., 0.1, 0.0, 0.0], [-0.3, -0.2, 0.0, 0.0], [0., 0., 0., 0.5]]
+ D3 = [[0., -0.8, 0.], [-0.3, 0., 0.], [0., 0., 0.4]]
+ sys1 = StateSpace(A1, B1, C1, D1)
+ sys2 = StateSpace(A2, B2, C2, D2)
+ sys3 = StateSpace(A3, B3, C3, D3)
+ sys3c = sys1.append(sys2)
+ np.testing.assert_array_almost_equal(sys3.A, sys3c.A)
+ np.testing.assert_array_almost_equal(sys3.B, sys3c.B)
+ np.testing.assert_array_almost_equal(sys3.C, sys3c.C)
+ np.testing.assert_array_almost_equal(sys3.D, sys3c.D)
+
+ def testAppendTF(self):
+ """Test appending a state-space system with a tf"""
+ A1 = [[-2, 0.5, 0], [0.5, -0.3, 0], [0, 0, -0.1]]
+ B1 = [[0.3...
[truncated message content] |
|
From: <mur...@us...> - 2013-06-09 15:50:52
|
Revision: 275
http://sourceforge.net/p/python-control/code/275
Author: murrayrm
Date: 2013-06-09 15:50:50 +0000 (Sun, 09 Jun 2013)
Log Message:
-----------
Added in patches from Ryan Krauss <rk...@si...>:
* now possible to call a transfer function with a complex argument and get
back the value of the transfer function (returns a matrix in MIMO case)
* default argument for second transfer function in feedback() call is now 1.
Note that this will only work for SISO systems (wasn't sure extending to
MIMO was that useful).
Modified Paths:
--------------
trunk/ChangeLog
trunk/src/bdalg.py
trunk/src/frdata.py
trunk/src/statesp.py
trunk/src/xferfcn.py
trunk/tests/bdalg_test.py
trunk/tests/frd_test.py
trunk/tests/xferfcn_test.py
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2013-06-09 14:57:59 UTC (rev 274)
+++ trunk/ChangeLog 2013-06-09 15:50:50 UTC (rev 275)
@@ -1,3 +1,28 @@
+2013-06-09 Richard Murray <murray@altura-2.local>
+
+ * tests/bdalg_test.py (TestFeedback.testScalarSS)
+ (TestFeedback.testScalarTF): added simple test of default argument
+ to feedback() for state space and transfer function systems
+
+ * tests/frd_test.py (TestFRD.testFeedback): added test of default
+ argument to feedback() for FRD type
+
+ * src/bdalg.py (feedback):
+ * src/frdata.py (FRD.feedback):
+ * src/statesp.py (StateSpace.feedback):
+ patched in default value of 1 for other system in feedback().
+ Contributed by Ryan Krauss <rk...@si...>
+
+ * tests/xferfcn_test.py (TestXferFcn.testEvalFrMIMO)
+ (TestXferFcn.testEvalFrSISO): added tests of __call__ version of
+ evalfr.
+
+ * src/xferfcn.py (TransferFunction.__call__): patched in new
+ __call__ functionality from Ryan Krauss <rk...@si...>, but
+ modified to use evalfr (works with MIMO). SISO transfer functions
+ return a single complex number, MIMO transfer functions return a
+ matrix of complex numbers.
+
2013-05-23 Rene van Paassen <ren...@gm...>
* src/margin.py: re-vamped stability_margins function. Now
Modified: trunk/src/bdalg.py
===================================================================
--- trunk/src/bdalg.py 2013-06-09 14:57:59 UTC (rev 274)
+++ trunk/src/bdalg.py 2013-06-09 15:50:50 UTC (rev 275)
@@ -172,7 +172,8 @@
return -sys;
-def feedback(sys1, sys2, sign=-1):
+#! TODO: expand to allow sys2 default to work in MIMO case?
+def feedback(sys1, sys2=1, sign=-1):
"""
Feedback interconnection between two I/O systems.
Modified: trunk/src/frdata.py
===================================================================
--- trunk/src/frdata.py 2013-06-09 14:57:59 UTC (rev 274)
+++ trunk/src/frdata.py 2013-06-09 15:50:50 UTC (rev 275)
@@ -393,7 +393,7 @@
return mag, phase, omega
- def feedback(self, other, sign=-1):
+ def feedback(self, other=1, sign=-1):
"""Feedback interconnection between two FRD objects."""
other = _convertToFRD(other, omega=self.omega)
Modified: trunk/src/statesp.py
===================================================================
--- trunk/src/statesp.py 2013-06-09 14:57:59 UTC (rev 274)
+++ trunk/src/statesp.py 2013-06-09 15:50:50 UTC (rev 275)
@@ -447,7 +447,7 @@
return roots(num)
# Feedback around a state space system
- def feedback(self, other, sign=-1):
+ def feedback(self, other=1, sign=-1):
"""Feedback interconnection between two LTI systems."""
other = _convertToStateSpace(other)
Modified: trunk/src/xferfcn.py
===================================================================
--- trunk/src/xferfcn.py 2013-06-09 14:57:59 UTC (rev 274)
+++ trunk/src/xferfcn.py 2013-06-09 15:50:50 UTC (rev 275)
@@ -81,7 +81,7 @@
# External function declarations
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
+ where, delete, real, poly, poly1d
from scipy.signal import lti
from copy import deepcopy
from warnings import warn
@@ -231,7 +231,21 @@
self.den = den
self._truncatecoeff()
-
+
+ def __call__(self, s):
+ """Evaluate the system's transfer function for a complex vairable
+
+ For a SISO transfer function, returns the value of the
+ transfer function. For a MIMO transfer fuction, returns a
+ matrix of values evaluated at complex variable s."""
+
+ if (self.inputs > 1 or self.outputs > 1):
+ # MIMO transfer function, return a matrix
+ return self.horner(s)
+ else:
+ # SISO transfer function, return a scalar
+ return self.horner(s)[0][0]
+
def _truncatecoeff(self):
"""Remove extraneous zero coefficients from num and den.
@@ -604,7 +618,7 @@
#for now, just give zeros of a SISO tf
return roots(self.num[0][0])
- def feedback(self, other, sign=-1):
+ def feedback(self, other=1, sign=-1):
"""Feedback interconnection between two LTI objects."""
other = _convertToTransferFunction(other)
Modified: trunk/tests/bdalg_test.py
===================================================================
--- trunk/tests/bdalg_test.py 2013-06-09 14:57:59 UTC (rev 274)
+++ trunk/tests/bdalg_test.py 2013-06-09 15:50:50 UTC (rev 275)
@@ -50,6 +50,14 @@
np.testing.assert_array_almost_equal(ans2.C, [[2.5, 0.]])
np.testing.assert_array_almost_equal(ans2.D, [[2.5]])
+ # Make sure default arugments work as well
+ ans3 = feedback(self.sys2, 1)
+ ans4 = feedback(self.sys2)
+ np.testing.assert_array_almost_equal(ans3.A, ans4.A)
+ np.testing.assert_array_almost_equal(ans3.B, ans4.B)
+ np.testing.assert_array_almost_equal(ans3.C, ans4.C)
+ np.testing.assert_array_almost_equal(ans3.D, ans4.D)
+
def testScalarTF(self):
"""Scalar system with transfer function feedback block."""
@@ -61,6 +69,12 @@
np.testing.assert_array_almost_equal(ans2.num, [[[2.5, 5., 7.5]]])
np.testing.assert_array_almost_equal(ans2.den, [[[1., -0.5, -2.]]])
+ # Make sure default arugments work as well
+ ans3 = feedback(self.sys1, 1)
+ ans4 = feedback(self.sys1)
+ np.testing.assert_array_almost_equal(ans3.num, ans4.num)
+ np.testing.assert_array_almost_equal(ans3.den, ans4.den)
+
def testSSScalar(self):
"""State space system with scalar feedback block."""
Modified: trunk/tests/frd_test.py
===================================================================
--- trunk/tests/frd_test.py 2013-06-09 14:57:59 UTC (rev 274)
+++ trunk/tests/frd_test.py 2013-06-09 15:50:50 UTC (rev 275)
@@ -84,6 +84,11 @@
np.testing.assert_array_almost_equal(
f1.feedback(1).freqresp([0.1, 1.0, 10])[0],
h1.feedback(1).freqresp([0.1, 1.0, 10])[0])
+
+ # Make sure default argument also works
+ np.testing.assert_array_almost_equal(
+ f1.feedback().freqresp([0.1, 1.0, 10])[0],
+ h1.feedback().freqresp([0.1, 1.0, 10])[0])
def testFeedback2(self):
h2 = StateSpace([[-1.0, 0], [0, -2.0]], [[0.4], [0.1]],
Modified: trunk/tests/xferfcn_test.py
===================================================================
--- trunk/tests/xferfcn_test.py 2013-06-09 14:57:59 UTC (rev 274)
+++ trunk/tests/xferfcn_test.py 2013-06-09 15:50:50 UTC (rev 275)
@@ -324,6 +324,11 @@
np.testing.assert_array_almost_equal(sys.evalfr(32.),
np.array([[0.00281959302585077 - 0.030628473607392j]]))
+ # Test call version as well
+ np.testing.assert_almost_equal(sys(1.j), -0.5 - 0.5j)
+ np.testing.assert_almost_equal(sys(32.j),
+ 0.00281959302585077 - 0.030628473607392j)
+
def testEvalFrMIMO(self):
"""Evaluate the frequency response of a MIMO system at one frequency."""
@@ -338,6 +343,9 @@
np.testing.assert_array_almost_equal(sys.evalfr(2.), resp)
+ # Test call version as well
+ np.testing.assert_array_almost_equal(sys(2.j), resp)
+
# Tests for TransferFunction.freqresp.
def testFreqRespSISO(self):
|
|
From: <mur...@us...> - 2013-06-10 15:38:54
|
Revision: 281
http://sourceforge.net/p/python-control/code/281
Author: murrayrm
Date: 2013-06-10 15:38:47 +0000 (Mon, 10 Jun 2013)
Log Message:
-----------
docstring updates
Modified Paths:
--------------
trunk/ChangeLog
trunk/src/freqplot.py
trunk/src/xferfcn.py
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2013-06-10 05:05:46 UTC (rev 280)
+++ trunk/ChangeLog 2013-06-10 15:38:47 UTC (rev 281)
@@ -1,3 +1,9 @@
+2013-06-10 Richard Murray <murray@altura-2.local>
+
+ * src/xferfcn.py (TransferFunction.horner): small fix to docstring
+
+ * src/freqplot.py (nyquist_plot): small fix to docstring
+
2013-06-09 Richard Murray <murray@altura-2.local>
* tests/bdalg_test.py (TestFeedback.testScalarSS)
Modified: trunk/src/freqplot.py
===================================================================
--- trunk/src/freqplot.py 2013-06-10 05:05:46 UTC (rev 280)
+++ trunk/src/freqplot.py 2013-06-10 15:38:47 UTC (rev 281)
@@ -204,7 +204,7 @@
Examples
--------
>>> sys = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.")
- >>> real, imag, freq = nyquist(sys)
+ >>> real, imag, freq = nyquist_plot(sys)
"""
# If argument was a singleton, turn it into a list
if (not getattr(syslist, '__iter__', False)):
Modified: trunk/src/xferfcn.py
===================================================================
--- trunk/src/xferfcn.py 2013-06-10 05:05:46 UTC (rev 280)
+++ trunk/src/xferfcn.py 2013-06-10 15:38:47 UTC (rev 281)
@@ -545,10 +545,10 @@
return self.horner(s)
def horner(self, s):
- '''Evaluate the systems's transfer function for a complex variable
+ """Evaluate the systems's transfer function for a complex variable
Returns a matrix of values evaluated at complex variable s.
- '''
+ """
# Preallocate the output.
if getattr(s, '__iter__', False):
|
|
From: <mur...@us...> - 2013-06-25 06:07:53
|
Revision: 288
http://sourceforge.net/p/python-control/code/288
Author: murrayrm
Date: 2013-06-25 06:07:45 +0000 (Tue, 25 Jun 2013)
Log Message:
-----------
updated exception handling to be python2.6+ and python3 compatible
Modified Paths:
--------------
trunk/ChangeLog
trunk/src/margins.py
trunk/src/statesp.py
trunk/src/xferfcn.py
trunk/tests/minreal_test.py
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2013-06-19 01:50:10 UTC (rev 287)
+++ trunk/ChangeLog 2013-06-25 06:07:45 UTC (rev 288)
@@ -1,3 +1,10 @@
+2013-06-24 Richard Murray <murray@altura-2.local>
+
+ * tests/minreal_test.py (TestMinreal.testMinrealBrute),
+ src/xferfcn.py (_convertToTransferFunction), src/statesp.py
+ (StateSpace.__rmul__), src/margins.py (stability_margins): changed
+ exception syntax to be python3 compatibile (works for python2.6+)
+
2013-06-10 Richard Murray <murray@altura-2.local>
* src/xferfcn.py (TransferFunction.horner): small fix to docstring
Modified: trunk/src/margins.py
===================================================================
--- trunk/src/margins.py 2013-06-19 01:50:10 UTC (rev 287)
+++ trunk/src/margins.py 2013-06-25 06:07:45 UTC (rev 288)
@@ -127,7 +127,7 @@
sys = frdata.FRD(mag*np.exp((1j/360.)*phase), omega, smooth=True)
else:
sys = xferfcn._convertToTransferFunction(sysdata)
- except Exception, e:
+ except Exception as e:
print (e)
raise ValueError("Margin sysdata must be either a linear system or "
"a 3-sequence of mag, phase, omega.")
Modified: trunk/src/statesp.py
===================================================================
--- trunk/src/statesp.py 2013-06-19 01:50:10 UTC (rev 287)
+++ trunk/src/statesp.py 2013-06-25 06:07:45 UTC (rev 288)
@@ -351,7 +351,7 @@
D = X * self.D
return StateSpace(self.A, self.B, C, D, self.dt)
- except Exception, e:
+ except Exception as e:
print(e)
pass
raise TypeError("can't interconnect systems")
Modified: trunk/src/xferfcn.py
===================================================================
--- trunk/src/xferfcn.py 2013-06-19 01:50:10 UTC (rev 287)
+++ trunk/src/xferfcn.py 2013-06-25 06:07:45 UTC (rev 288)
@@ -1053,7 +1053,7 @@
num = [[[D[i,j]] for j in range(inputs)] for i in range(outputs)]
den = [[[1] for j in range(inputs)] for i in range(outputs)]
return TransferFunction(num, den)
- except Exception, e:
+ except Exception as e:
print("Failure to assume argument is matrix-like in"
" _convertToTransferFunction, result %s" % e)
Modified: trunk/tests/minreal_test.py
===================================================================
--- trunk/tests/minreal_test.py 2013-06-19 01:50:10 UTC (rev 287)
+++ trunk/tests/minreal_test.py 2013-06-25 06:07:45 UTC (rev 288)
@@ -53,7 +53,7 @@
self.assert_numden_almost_equal(
ht1.num[0][0], ht2.num[0][0],
ht1.den[0][0], ht2.den[0][0])
- except Exception, e:
+ except Exception as e:
# for larger systems, the tf minreal's
# the original rss, but not the balanced one
if n < 6:
|
|
From: <mur...@us...> - 2013-07-15 04:40:20
|
Revision: 289
http://sourceforge.net/p/python-control/code/289
Author: murrayrm
Date: 2013-07-15 04:39:59 +0000 (Mon, 15 Jul 2013)
Log Message:
-----------
Updated sample_system() function to check for proper version of SciPy. If cont2discrete is not avialable (requires scipy 0.10.0+), generate a warning.
Modified Paths:
--------------
trunk/ChangeLog
trunk/src/dtime.py
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2013-06-25 06:07:45 UTC (rev 288)
+++ trunk/ChangeLog 2013-07-15 04:39:59 UTC (rev 289)
@@ -1,3 +1,9 @@
+2013-07-14 Richard Murray <murray@altura-2.local>
+
+ * src/dtime.py (sample_system): added check to make sure
+ cont2discrete is available in scipy.signal
+ (sample_system): added docstring
+
2013-06-24 Richard Murray <murray@altura-2.local>
* tests/minreal_test.py (TestMinreal.testMinrealBrute),
Modified: trunk/src/dtime.py
===================================================================
--- trunk/src/dtime.py 2013-06-25 06:07:45 UTC (rev 288)
+++ trunk/src/dtime.py 2013-07-15 04:39:59 UTC (rev 289)
@@ -47,7 +47,7 @@
"""
-from scipy.signal import zpk2tf, tf2zpk, cont2discrete
+from scipy.signal import zpk2tf, tf2zpk
import numpy as np
from cmath import exp
from warnings import warn
@@ -57,8 +57,39 @@
# Sample a continuous time system
def sample_system(sysc, Ts, method='matched'):
- # TODO: add docstring
+ """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
+ ----------
+ sysc : linsys
+ Continuous time system to be converted
+ Ts : real
+ Sampling period
+ method : string
+ Method to use for conversion: 'matched' (default), 'tustin', 'zoh'
+
+ Returns
+ -------
+ sysd : linsys
+ Discrete time system, with sampling rate Ts
+
+ Notes
+ -----
+ 1. The conversion methods 'tustin' and 'zoh' require the
+ cont2discrete() function, including in SciPy 0.10.0 and above.
+
+ 2. Additional methods 'foh' and 'impulse' are planned for future
+ implementation.
+
+ Examples
+ --------
+ >>> sysc = TransferFunction([1], [1, 2, 1])
+ >>> sysd = sample_system(sysc, 1, method='matched')
+ """
+
# Make sure we have a continuous time system
if not isctime(sysc):
raise ValueError("First argument must be continuous time system")
@@ -77,14 +108,22 @@
sysd = _c2dmatched(sysc, Ts)
elif method == 'tustin':
- sys = [sysc.num[0][0], sysc.den[0][0]]
- scipySysD = cont2discrete(sys, Ts, method='bilinear')
- sysd = TransferFunction(scipySysD[0][0], scipySysD[1], dt)
-
+ try:
+ from scipy.signal import cont2discrete
+ sys = [sysc.num[0][0], sysc.den[0][0]]
+ scipySysD = cont2discrete(sys, Ts, method='bilinear')
+ sysd = TransferFunction(scipySysD[0][0], scipySysD[1], dt)
+ except ImportError:
+ raise TypeError("cont2discrete not found in scipy.signal; upgrade to v0.10.0+")
+
elif method == 'zoh':
- sys = [sysc.num[0][0], sysc.den[0][0]]
- scipySysD = cont2discrete(sys, Ts, method='zoh')
- sysd = TransferFunction(scipySysD[0][0],scipySysD[1], dt)
+ try:
+ from scipy.signal import cont2discrete
+ sys = [sysc.num[0][0], sysc.den[0][0]]
+ scipySysD = cont2discrete(sys, Ts, method='zoh')
+ sysd = TransferFunction(scipySysD[0][0],scipySysD[1], dt)
+ except ImportError:
+ raise TypeError("cont2discrete not found in scipy.signal; upgrade to v0.10.0+")
elif method == 'foh' or method == 'impulse':
raise ValueError("Method not developed yet")
|
|
From: <mur...@us...> - 2013-07-16 05:14:26
|
Revision: 290
http://sourceforge.net/p/python-control/code/290
Author: murrayrm
Date: 2013-07-16 05:14:17 +0000 (Tue, 16 Jul 2013)
Log Message:
-----------
updated division operators for FRD and TransferFunction to be python3 compatible
Modified Paths:
--------------
trunk/ChangeLog
trunk/src/frdata.py
trunk/src/xferfcn.py
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2013-07-15 04:39:59 UTC (rev 289)
+++ trunk/ChangeLog 2013-07-16 05:14:17 UTC (rev 290)
@@ -1,3 +1,8 @@
+2013-07-15 Richard Murray <murray@altura-2.local>
+
+ * src/frdata.py, src/xferfcn.py: updated to use __truediv__ and
+ __rtruediv__ for compatibility with python3
+
2013-07-14 Richard Murray <murray@altura-2.local>
* src/dtime.py (sample_system): added check to make sure
Modified: trunk/src/frdata.py
===================================================================
--- trunk/src/frdata.py 2013-07-15 04:39:59 UTC (rev 289)
+++ trunk/src/frdata.py 2013-07-16 05:14:17 UTC (rev 290)
@@ -1,3 +1,4 @@
+from __future__ import division
"""frdata.py
Frequency response data representation and functions.
@@ -17,8 +18,8 @@
FRD.__rsub__
FRD.__mul__
FRD.__rmul__
-FRD.__div__
-FRD.__rdiv__
+FRD.__truediv__
+FRD.__rtruediv__
FRD.evalfr
FRD.freqresp
FRD.pole
@@ -289,7 +290,7 @@
return FRD(fresp, self.omega)
# TODO: Division of MIMO transfer function objects is not written yet.
- def __div__(self, other):
+ def __truediv__(self, other):
"""Divide two LTI objects."""
if isinstance(other, (int, float, complex)):
@@ -301,12 +302,12 @@
if (self.inputs > 1 or self.outputs > 1 or
other.inputs > 1 or other.outputs > 1):
raise NotImplementedError(
- "FRD.__div__ is currently implemented only for SISO systems.")
+ "FRD.__truediv__ is currently implemented only for SISO systems.")
return FRD(self.fresp/other.fresp, self.omega)
# TODO: Division of MIMO transfer function objects is not written yet.
- def __rdiv__(self, other):
+ def __rtruediv__(self, other):
"""Right divide two LTI objects."""
if isinstance(other, (int, float, complex)):
return FRD(other / self.fresp, self.omega)
@@ -316,7 +317,7 @@
if (self.inputs > 1 or self.outputs > 1 or
other.inputs > 1 or other.outputs > 1):
raise NotImplementedError(
- "FRD.__rdiv__ is currently implemented only for SISO systems.")
+ "FRD.__rtruediv__ is currently implemented only for SISO systems.")
return other / self
Modified: trunk/src/xferfcn.py
===================================================================
--- trunk/src/xferfcn.py 2013-07-15 04:39:59 UTC (rev 289)
+++ trunk/src/xferfcn.py 2013-07-16 05:14:17 UTC (rev 290)
@@ -37,6 +37,7 @@
# Python 3 compatability (needs to go here)
from __future__ import print_function
+from __future__ import division
"""Copyright (c) 2010 by California Institute of Technology
All rights reserved.
@@ -470,7 +471,7 @@
return TransferFunction(num, den, dt)
# TODO: Division of MIMO transfer function objects is not written yet.
- def __div__(self, other):
+ def __truediv__(self, other):
"""Divide two LTI objects."""
if isinstance(other, (int, float, complex)):
@@ -482,7 +483,7 @@
if (self.inputs > 1 or self.outputs > 1 or
other.inputs > 1 or other.outputs > 1):
- raise NotImplementedError("TransferFunction.__div__ is currently \
+ raise NotImplementedError("TransferFunction.__truediv__ is currently \
implemented only for SISO systems.")
# Figure out the sampling time to use
@@ -499,7 +500,7 @@
return TransferFunction(num, den, dt)
# TODO: Division of MIMO transfer function objects is not written yet.
- def __rdiv__(self, other):
+ def __rtruediv__(self, other):
"""Right divide two LTI objects."""
if isinstance(other, (int, float, complex)):
other = _convertToTransferFunction(other, inputs=self.inputs,
@@ -509,7 +510,7 @@
if (self.inputs > 1 or self.outputs > 1 or
other.inputs > 1 or other.outputs > 1):
- raise NotImplementedError("TransferFunction.__rdiv__ is currently \
+ raise NotImplementedError("TransferFunction.__rtruediv__ is currently \
implemented only for SISO systems.")
return other / self
|
|
From: <mur...@us...> - 2014-03-22 18:09:56
|
Revision: 294
http://sourceforge.net/p/python-control/code/294
Author: murrayrm
Date: 2014-03-22 18:09:48 +0000 (Sat, 22 Mar 2014)
Log Message:
-----------
Fixed bug #5 ('dt' instead of 'Ts' in dtime.py) and added some unit tests for dtime and FRD bdalg.
Modified Paths:
--------------
trunk/ChangeLog
trunk/src/bdalg.py
trunk/src/dtime.py
trunk/src/frdata.py
trunk/tests/discrete_test.py
trunk/tests/frd_test.py
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2013-07-16 14:59:20 UTC (rev 293)
+++ trunk/ChangeLog 2014-03-22 18:09:48 UTC (rev 294)
@@ -1,3 +1,17 @@
+2014-03-22 Richard Murray <murray@sidamo.local>
+
+ * tests/discrete_test.py (TestDiscrete.test_sample_system): added
+ conversions using tustin and zoh to catch 'dt' bug #5
+
+ * src/dtime.py (sample_system): Changed 'dt' to 'Ts' to fix bug #4
+
+ * tests/frd_test.py (TestFRD.testbdalg): added unit tests for bdalg
+ functions operating on FRD data
+
+ * src/bdalg.py (series, parallel, negate, feedback): updated
+ documentation to note that FRD ovjects work.
+ (feedback): updated conversions to allow for FRD
+
2013-07-15 Richard Murray <murray@altura-2.local>
* src/frdata.py, src/xferfcn.py: updated to use __truediv__ and
Modified: trunk/src/bdalg.py
===================================================================
--- trunk/src/bdalg.py 2013-07-16 14:59:20 UTC (rev 293)
+++ trunk/src/bdalg.py 2014-03-22 18:09:48 UTC (rev 294)
@@ -56,14 +56,15 @@
import scipy as sp
import control.xferfcn as tf
import control.statesp as ss
+import control.frdata as frd
def series(sys1, sys2):
"""Return the series connection sys2 * sys1 for --> sys1 --> sys2 -->.
Parameters
----------
- sys1: scalar, StateSpace, or TransferFunction
- sys2: scalar, StateSpace, or TransferFunction
+ sys1: scalar, StateSpace, TransferFunction, or FRD
+ sys2: scalar, StateSpace, TransferFunction, or FRD
Returns
-------
@@ -105,8 +106,8 @@
Parameters
----------
- sys1: scalar, StateSpace, or TransferFunction
- sys2: scalar, StateSpace, or TransferFunction
+ sys1: scalar, StateSpace, TransferFunction, or FRD
+ sys2: scalar, StateSpace, TransferFunction, or FRD
Returns
-------
@@ -148,7 +149,7 @@
Parameters
----------
- sys: StateSpace or TransferFunction
+ sys: StateSpace, TransferFunction or FRD
Returns
-------
@@ -179,9 +180,9 @@
Parameters
----------
- sys1: scalar, StateSpace, or TransferFunction
+ sys1: scalar, StateSpace, TransferFunction, FRD
The primary plant.
- sys2: scalar, StateSpace, or TransferFunction
+ sys2: scalar, StateSpace, TransferFunction, FRD
The feedback plant (often a feedback controller).
sign: scalar
The sign of feedback. `sign` = -1 indicates negative feedback, and
@@ -219,13 +220,13 @@
# Check for correct input types.
if not isinstance(sys1, (int, float, complex, tf.TransferFunction,
- ss.StateSpace)):
- raise TypeError("sys1 must be a TransferFunction or StateSpace object, \
-or a scalar.")
+ ss.StateSpace)):
+ raise TypeError("sys1 must be a TransferFunction or StateSpace " +
+ "object, or a scalar.")
if not isinstance(sys2, (int, float, complex, tf.TransferFunction,
- ss.StateSpace)):
- raise TypeError("sys2 must be a TransferFunction or StateSpace object, \
-or a scalar.")
+ ss.StateSpace)):
+ raise TypeError("sys2 must be a TransferFunction or StateSpace " +
+ "object, or a scalar.")
# If sys1 is a scalar, convert it to the appropriate LTI type so that we can
# its feedback member function.
@@ -234,6 +235,8 @@
sys1 = tf._convertToTransferFunction(sys1)
elif isinstance(sys2, ss.StateSpace):
sys1 = ss._convertToStateSpace(sys1)
+ elif isinstance(sys2, frd.FRD):
+ sys1 = ss._convertToFRD(sys1)
else: # sys2 is a scalar.
sys1 = tf._convertToTransferFunction(sys1)
sys2 = tf._convertToTransferFunction(sys2)
Modified: trunk/src/dtime.py
===================================================================
--- trunk/src/dtime.py 2013-07-16 14:59:20 UTC (rev 293)
+++ trunk/src/dtime.py 2014-03-22 18:09:48 UTC (rev 294)
@@ -112,7 +112,7 @@
from scipy.signal import cont2discrete
sys = [sysc.num[0][0], sysc.den[0][0]]
scipySysD = cont2discrete(sys, Ts, method='bilinear')
- sysd = TransferFunction(scipySysD[0][0], scipySysD[1], dt)
+ sysd = TransferFunction(scipySysD[0][0], scipySysD[1], Ts)
except ImportError:
raise TypeError("cont2discrete not found in scipy.signal; upgrade to v0.10.0+")
@@ -121,7 +121,7 @@
from scipy.signal import cont2discrete
sys = [sysc.num[0][0], sysc.den[0][0]]
scipySysD = cont2discrete(sys, Ts, method='zoh')
- sysd = TransferFunction(scipySysD[0][0],scipySysD[1], dt)
+ sysd = TransferFunction(scipySysD[0][0],scipySysD[1], Ts)
except ImportError:
raise TypeError("cont2discrete not found in scipy.signal; upgrade to v0.10.0+")
Modified: trunk/src/frdata.py
===================================================================
--- trunk/src/frdata.py 2013-07-16 14:59:20 UTC (rev 293)
+++ trunk/src/frdata.py 2014-03-22 18:09:48 UTC (rev 294)
@@ -451,7 +451,6 @@
"""
if isinstance(sys, FRD):
-
omega.sort()
if (abs(omega - sys.omega) < FRD.epsw).all():
# frequencies match, and system was already frd; simply use
Modified: trunk/tests/discrete_test.py
===================================================================
--- trunk/tests/discrete_test.py 2013-07-16 14:59:20 UTC (rev 293)
+++ trunk/tests/discrete_test.py 2014-03-22 18:09:48 UTC (rev 294)
@@ -264,6 +264,12 @@
for sysc in (self.siso_ss1, self.siso_ss1c, self.siso_tf1c):
sysd = sample_system(sysc, 1, method='matched')
self.assertEqual(sysd.dt, 1)
+
+ sysd = sample_system(sysc, 1, method='tustin')
+ self.assertEqual(sysd.dt, 1)
+
+ sysd = sample_system(sysc, 1, method='zoh')
+ self.assertEqual(sysd.dt, 1)
# TODO: put in other generic checks
# TODO: check results of converstion
Modified: trunk/tests/frd_test.py
===================================================================
--- trunk/tests/frd_test.py 2013-07-16 14:59:20 UTC (rev 293)
+++ trunk/tests/frd_test.py 2014-03-22 18:09:48 UTC (rev 294)
@@ -10,6 +10,7 @@
from control.xferfcn import TransferFunction
from control.frdata import FRD, _convertToFRD
from control.matlab import bode
+import control.bdalg as bdalg
import control.freqplot
import matplotlib.pyplot as plt
@@ -106,7 +107,43 @@
(f1 / h2).freqresp([0.1, 1.0, 10])[1],
(h1 / h2).freqresp([0.1, 1.0, 10])[1])
# the reverse does not work
-
+
+ def testbdalg(self):
+ # get two SISO transfer functions
+ h1 = TransferFunction([1], [1, 2, 2])
+ h2 = TransferFunction([1], [0.1, 1])
+ omega = np.logspace(-1, 2, 10)
+ f1 = FRD(h1, omega)
+ f2 = FRD(h2, omega)
+
+ np.testing.assert_array_almost_equal(
+ (bdalg.series(f1, f2)).freqresp([0.1, 1.0, 10])[0],
+ (bdalg.series(h1, h2)).freqresp([0.1, 1.0, 10])[0])
+
+ np.testing.assert_array_almost_equal(
+ (bdalg.parallel(f1, f2)).freqresp([0.1, 1.0, 10])[0],
+ (bdalg.parallel(h1, h2)).freqresp([0.1, 1.0, 10])[0])
+
+ np.testing.assert_array_almost_equal(
+ (bdalg.feedback(f1, f2)).freqresp([0.1, 1.0, 10])[0],
+ (bdalg.feedback(h1, h2)).freqresp([0.1, 1.0, 10])[0])
+
+ np.testing.assert_array_almost_equal(
+ (bdalg.negate(f1)).freqresp([0.1, 1.0, 10])[0],
+ (bdalg.negate(h1)).freqresp([0.1, 1.0, 10])[0])
+
+# append() and connect() not implemented for FRD objects
+# np.testing.assert_array_almost_equal(
+# (bdalg.append(f1, f2)).freqresp([0.1, 1.0, 10])[0],
+# (bdalg.append(h1, h2)).freqresp([0.1, 1.0, 10])[0])
+#
+# f3 = bdalg.append(f1, f2, f2)
+# h3 = bdalg.append(h1, h2, h2)
+# Q = np.mat([ [1, 2], [2, -1] ])
+# np.testing.assert_array_almost_equal(
+# (bdalg.connect(f3, Q, [2], [1])).freqresp([0.1, 1.0, 10])[0],
+# (bdalg.connect(h3, Q, [2], [1])).freqresp([0.1, 1.0, 10])[0])
+
def testFeedback(self):
h1 = TransferFunction([1], [1, 2, 2])
omega = np.logspace(-1, 2, 10)
|
|
From: <mur...@us...> - 2014-03-22 18:31:04
|
Revision: 295
http://sourceforge.net/p/python-control/code/295
Author: murrayrm
Date: 2014-03-22 18:30:41 +0000 (Sat, 22 Mar 2014)
Log Message:
-----------
allow FRD for feedback()
Modified Paths:
--------------
trunk/ChangeLog
trunk/src/bdalg.py
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2014-03-22 18:09:48 UTC (rev 294)
+++ trunk/ChangeLog 2014-03-22 18:30:41 UTC (rev 295)
@@ -1,5 +1,9 @@
2014-03-22 Richard Murray <murray@sidamo.local>
+ * src/bdalg.py (feedback): allow FRD for feedback()
+
+2014-03-22 Richard Murray <murray@sidamo.local>
+
* tests/discrete_test.py (TestDiscrete.test_sample_system): added
conversions using tustin and zoh to catch 'dt' bug #5
Modified: trunk/src/bdalg.py
===================================================================
--- trunk/src/bdalg.py 2014-03-22 18:09:48 UTC (rev 294)
+++ trunk/src/bdalg.py 2014-03-22 18:30:41 UTC (rev 295)
@@ -220,13 +220,13 @@
# Check for correct input types.
if not isinstance(sys1, (int, float, complex, tf.TransferFunction,
- ss.StateSpace)):
- raise TypeError("sys1 must be a TransferFunction or StateSpace " +
- "object, or a scalar.")
+ 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)):
- raise TypeError("sys2 must be a TransferFunction or StateSpace " +
- "object, or a scalar.")
+ 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.
|
|
From: <mur...@us...> - 2014-03-22 19:31:47
|
Revision: 296
http://sourceforge.net/p/python-control/code/296
Author: murrayrm
Date: 2014-03-22 19:31:36 +0000 (Sat, 22 Mar 2014)
Log Message:
-----------
allow 5 arguments for matlab.ss() command to create discrete time system
Modified Paths:
--------------
trunk/ChangeLog
trunk/external/yottalab.py
trunk/src/matlab.py
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2014-03-22 18:30:41 UTC (rev 295)
+++ trunk/ChangeLog 2014-03-22 19:31:36 UTC (rev 296)
@@ -1,5 +1,10 @@
2014-03-22 Richard Murray <murray@sidamo.local>
+ * src/matlab.py (ss): allow five arguments to create a discrete time
+ system. From Roberto Bucher, 12 Jan 2014
+
+2014-03-22 Richard Murray <murray@sidamo.local>
+
* src/bdalg.py (feedback): allow FRD for feedback()
2014-03-22 Richard Murray <murray@sidamo.local>
Modified: trunk/external/yottalab.py
===================================================================
--- trunk/external/yottalab.py 2014-03-22 18:30:41 UTC (rev 295)
+++ trunk/external/yottalab.py 2014-03-22 19:31:36 UTC (rev 296)
@@ -1,5 +1,3 @@
-# yottalab.py - Roberto Bucher's yottalab library
-
"""
This is a procedural interface to the yttalab library
@@ -8,104 +6,314 @@
The following commands are provided:
Design and plot commands
- acker - pole placement using Ackermann method
- c2d - contimous to discrete time conversion
+ bb_c2d - contimous to discrete time conversion
+ d2c - discrete to continous time conversion
+ bb_dare - Solve Riccati equation for discrete time systems
+ bb_dlqr - discrete linear quadratic regulator
+ dsimul - simulate discrete time systems
+ dstep - step response (plot) of discrete time systems
+ dimpulse - imoulse response (plot) of discrete time systems
+ bb_step - step response (plot) of continous time systems
+
full_obs - full order observer
red_obs - reduced order observer
comp_form - state feedback controller+observer in compact form
comp_form_i - state feedback controller+observer+integ in compact form
- dsimul - simulate discrete time systems
- dstep - step response (plot) of discrete time systems
- dimpulse - imoulse response (plot) of discrete time systems
- bb_step - step response (plot) of continous time systems
sysctr - system+controller+observer+feedback
- care - Solve Riccati equation for contimous time systems
- dare - Solve Riccati equation for discrete time systems
- dlqr - discrete linear quadratic regulator
- minreal - minimal state space representation
-
+ set_aw - introduce anti-windup into controller
+ bb_dcgain - return the steady state value of the step response
+
"""
from matplotlib.pylab import *
-from control.matlab import *
+from control import *
from numpy import hstack,vstack,pi
from scipy import zeros,ones,eye,mat,shape,size,size, \
arange,real,poly,array,diag
-from scipy.linalg import det,inv,expm,eig,eigvals
+from scipy.linalg import det,inv,expm,eig,eigvals,logm
import numpy as np
import scipy as sp
-from slycot import sb02od, tb03ad
+from slycot import sb02od
+# from scipy.signal import BadCoefficients
+# import warnings
+# warnings.filterwarnings('ignore',category=BadCoefficients)
-def acker(A,B,poles):
- """Pole placemenmt using Ackermann method
+def bb_c2d(sys,Ts,method='zoh'):
+ """Continous to discrete conversion with ZOH method
Call:
- k=acker(A,B,poles)
+ sysd=c2d(sys,Ts,method='zoh')
Parameters
----------
- A, B : State and input matrix of the system
- poles: desired poles
+ sys : System in statespace or Tf form
+ Ts: Sampling Time
+ method: 'zoh', 'bi' or 'matched'
Returns
-------
- k: matrix
- State feedback gains
+ sysd: ss or Tf system
+ Discrete system
"""
- a=mat(A)
- b=mat(B)
- p=real(poly(poles))
- ct=ctrb(A,B)
- if det(ct)==0:
- k=0
- print "Pole placement invalid"
+ flag = 0
+ if isinstance(sys, TransferFunction):
+ sys=tf2ss(sys)
+ flag=1
+
+ a=sys.A
+ b=sys.B
+ c=sys.C
+ d=sys.D
+ n=shape(a)[0]
+ nb=shape(b)[1]
+ nc=shape(c)[0]
+
+ if method=='zoh':
+ ztmp=zeros((nb,n+nb))
+ tmp=hstack((a,b))
+ tmp=vstack((tmp,ztmp))
+ tmp=expm(tmp*Ts)
+ A=tmp[0:n,0:n]
+ B=tmp[0:n,n:n+nb]
+ C=c
+ D=d
+ elif method=='bi':
+ a=mat(a)
+ b=mat(b)
+ c=mat(c)
+ d=mat(d)
+ IT=mat(2/Ts*eye(n,n))
+ A=(IT+a)*inv(IT-a)
+ iab=inv(IT-a)*b
+ tk=2/sqrt(Ts)
+ B=tk*iab
+ C=tk*(c*inv(IT-a))
+ D=d+c*iab
+ elif method=='matched':
+ if nb!=1 and nc!=1:
+ print "System is not SISO"
+ return
+ p=exp(sys.poles*Ts)
+ z=exp(sys.zeros*Ts)
+ infinite_zeros = len(sys.poles) - len(sys.zeros) - 1
+ for i in range(0,infinite_zeros):
+ z=hstack((z,-1))
+ [A,B,C,D]=zpk2ss(z,p,1)
+ sysd=StateSpace(A,B,C,D,Ts)
+ cg = dcgain(sys)
+ dg = dcgain(sysd)
+ [A,B,C,D]=zpk2ss(z,p,cg/dg)
else:
- n=size(p)
- pmat=p[n-1]*a**0
- for i in arange(1,n):
- pmat=pmat+p[n-i-1]*a**i
- k=inv(ct)*pmat
- k=k[-1][:]
- return k
+ print "Method not supported"
+ return
+
+ sysd=StateSpace(A,B,C,D,Ts)
+ if flag==1:
+ sysd=ss2tf(sysd)
+ return sysd
-def c2d(sys,Ts):
+def d2c(sys,method='zoh'):
"""Continous to discrete conversion with ZOH method
Call:
- sysd=c2d(sys,Ts)
+ sysc=c2d(sys,method='log')
Parameters
----------
- sys : System in statespace or Tf form
- Ts: Sampling Time
+ sys : System in statespace or Tf form
+ method: 'zoh' or 'bi'
Returns
-------
- sysd: ss or Tf system
- Discrete system
+ sysc: continous system ss or tf
+
"""
flag = 0
- if type(sys).__name__=='TransferFunction':
+ if isinstance(sys, TransferFunction):
sys=tf2ss(sys)
flag=1
+
a=sys.A
b=sys.B
c=sys.C
d=sys.D
+ Ts=sys.dt
n=shape(a)[0]
nb=shape(b)[1]
- ztmp=zeros((nb,n+nb))
- tmp=hstack((a,b))
- tmp=vstack((tmp,ztmp))
- tmp=expm(tmp*Ts)
- a=tmp[0:n,0:n]
- b=tmp[0:n,n:n+nb]
- sysd=ss(a,b,c,d,Ts)
+ nc=shape(c)[0]
+ tol=1e-12
+
+ if method=='zoh':
+ if n==1:
+ if b[0,0]==1:
+ A=0
+ B=b/sys.dt
+ C=c
+ D=d
+ else:
+ tmp1=hstack((a,b))
+ tmp2=hstack((zeros((nb,n)),eye(nb)))
+ tmp=vstack((tmp1,tmp2))
+ s=logm(tmp)
+ s=s/Ts
+ if norm(imag(s),inf) > sqrt(sp.finfo(float).eps):
+ print "Warning: accuracy may be poor"
+ s=real(s)
+ A=s[0:n,0:n]
+ B=s[0:n,n:n+nb]
+ C=c
+ D=d
+ elif method=='bi':
+ a=mat(a)
+ b=mat(b)
+ c=mat(c)
+ d=mat(d)
+ poles=eigvals(a)
+ if any(abs(poles-1)<200*sp.finfo(float).eps):
+ print "d2c: some poles very close to one. May get bad results."
+
+ I=mat(eye(n,n))
+ tk = 2 / sqrt (Ts)
+ A = (2/Ts)*(a-I)*inv(a+I)
+ iab = inv(I+a)*b
+ B = tk*iab
+ C = tk*(c*inv(I+a))
+ D = d- (c*iab)
+ else:
+ print "Method not supported"
+ return
+
+ sysc=StateSpace(A,B,C,D)
if flag==1:
- sysd=ss2tf(sysd)
- return sysd
+ sysc=ss2tf(sysc)
+ return sysc
+def dsimul(sys,u):
+ """Simulate the discrete system sys
+ Only for discrete systems!!!
+
+ Call:
+ y=dsimul(sys,u)
+
+ Parameters
+ ----------
+ sys : Discrete System in State Space form
+ u : input vector
+ Returns
+ -------
+ y: ndarray
+ Simulation results
+
+ """
+ a=mat(sys.A)
+ b=mat(sys.B)
+ c=mat(sys.C)
+ d=mat(sys.D)
+ nx=shape(a)[0]
+ ns=shape(u)[1]
+ xk=zeros((nx,1))
+ for i in arange(0,ns):
+ uk=u[:,i]
+ xk_1=a*xk+b*uk
+ yk=c*xk+d*uk
+ xk=xk_1
+ if i==0:
+ y=yk
+ else:
+ y=hstack((y,yk))
+ y=array(y).T
+ return y
+
+def dstep(sys,Tf=10.0):
+ """Plot the step response of the discrete system sys
+ Only for discrete systems!!!
+
+ Call:
+ y=dstep(sys, [,Tf=final time]))
+
+ Parameters
+ ----------
+ sys : Discrete System in State Space form
+ Tf : Final simulation time
+
+ Returns
+ -------
+ Nothing
+
+ """
+ Ts=sys.dt
+ if Ts==0.0:
+ "Only discrete systems allowed!"
+ return
+
+ ns=int(Tf/Ts+1)
+ u=ones((1,ns))
+ y=dsimul(sys,u)
+ T=arange(0,Tf+Ts/2,Ts)
+ plot(T,y)
+ grid()
+ show()
+
+def dimpulse(sys,Tf=10.0):
+ """Plot the impulse response of the discrete system sys
+ Only for discrete systems!!!
+
+ Call:
+ y=dimpulse(sys,[,Tf=final time]))
+
+ Parameters
+ ----------
+ sys : Discrete System in State Space form
+ Tf : Final simulation time
+
+ Returns
+ -------
+ Nothing
+
+ """
+ Ts=sys.dt
+ if Ts==0.0:
+ "Only discrete systems allowed!"
+ return
+
+ ns=int(Tf/Ts+1)
+ u=zeros((1,ns))
+ u[0,0]=1/Ts
+ y=dsimul(sys,u)
+ T=arange(0,Tf+Ts/2,Ts)
+ plot(T,y)
+ grid()
+ show()
+
+# Step response (plot)
+def bb_step(sys,X0=None,Tf=None,Ts=0.001):
+ """Plot the step response of the continous system sys
+
+ Call:
+ y=bb_step(sys [,Tf=final time] [,Ts=time step])
+
+ Parameters
+ ----------
+ sys : Continous System in State Space form
+ X0: Initial state vector (not used yet)
+ Ts : sympling time
+ Tf : Final simulation time
+
+ Returns
+ -------
+ Nothing
+
+ """
+ if Tf==None:
+ vals = eigvals(sys.A)
+ r = min(abs(real(vals)))
+ if r < 1e-10:
+ r = 0.1
+ Tf = 7.0 / r
+ sysd=c2d(sys,Ts)
+ dstep(sysd,Tf=Tf)
+
def full_obs(sys,poles):
"""Full order observer of the system sys
@@ -123,14 +331,13 @@
Observer
"""
- if type(sys).__name__=='TransferFunction':
+ if isinstance(sys, TransferFunction):
"System must be in state space form"
return
a=mat(sys.A)
b=mat(sys.B)
c=mat(sys.C)
d=mat(sys.D)
- poles=mat(poles)
L=place(a.T,c.T,poles)
L=mat(L).T
Ao=a-L*c
@@ -139,7 +346,7 @@
m=shape(Bo)
Co=eye(n[0],n[1])
Do=zeros((n[0],m[1]))
- obs=ss(Ao,Bo,Co,Do,sys.Tsamp)
+ obs=StateSpace(Ao,Bo,Co,Do,sys.dt)
return obs
def red_obs(sys,T,poles):
@@ -160,7 +367,7 @@
Reduced order Observer
"""
- if type(sys).__name__=='TransferFunction':
+ if isinstance(sys, TransferFunction):
"System must be in state space form"
return
a=mat(sys.A)
@@ -169,7 +376,7 @@
d=mat(sys.D)
T=mat(T)
P=mat(vstack((c,T)))
- poles=mat(poles)
+ # poles=mat(poles)
invP=inv(P)
AA=P*a*invP
ny=shape(c)[0]
@@ -205,7 +412,7 @@
tmp5=mat(vstack((tmp5,tmp6)))
Dr=invP*tmp5*tmp4
- obs=ss(Ar,Br,Cr,Dr,sys.Tsamp)
+ obs=StateSpace(Ar,Br,Cr,Dr,sys.dt)
return obs
def comp_form(sys,obs,K):
@@ -242,7 +449,7 @@
Bc = hstack((Bu*X,By-Bu*X*K*Dy))
Cc = -X*K*mat(obs.C);
Dc = hstack((X,-X*K*Dy))
- contr = ss(Ac,Bc,Cc,Dc,sys.Tsamp)
+ contr = StateSpace(Ac,Bc,Cc,Dc,sys.dt)
return contr
def comp_form_i(sys,obs,K,Ts,Cy=[[1]]):
@@ -266,7 +473,7 @@
Controller
"""
- if sys.Tsamp==0.0:
+ if sys.dt==0.0:
print "contr_form_i works only with discrete systems!"
return
@@ -274,7 +481,7 @@
nu=shape(sys.B)[1]
nx=shape(sys.A)[0]
no=shape(obs.A)[0]
- ni=shape(Cy)[0]
+ ni=shape(mat(Cy))[0]
B_obsu = mat(obs.B[:,0:nu])
B_obsy = mat(obs.B[:,nu:nu+ny])
@@ -303,134 +510,9 @@
C_ctr=hstack((-X*K*c,-X*Ke))
D_ctr=hstack((zeros((nu,ni)),-X*K*D_obsy))
- contr=ss(A_ctr,B_ctr,C_ctr,D_ctr,sys.Tsamp)
+ contr=StateSpace(A_ctr,B_ctr,C_ctr,D_ctr,sys.dt)
return contr
-def dsimul(sys,u):
- """Simulate the discrete system sys
- Only for discrete systems!!!
-
- Call:
- y=dsimul(sys,u)
-
- Parameters
- ----------
- sys : Discrete System in State Space form
- u : input vector
- Returns
- -------
- y: ndarray
- Simulation results
-
- """
- a=mat(sys.A)
- b=mat(sys.B)
- c=mat(sys.C)
- d=mat(sys.D)
- nx=shape(a)[0]
- ns=shape(u)[1]
- xk=zeros((nx,1))
- for i in arange(0,ns):
- uk=u[:,i]
- xk_1=a*xk+b*uk
- yk=c*xk+d*uk
- xk=xk_1
- if i==0:
- y=yk
- else:
- y=hstack((y,yk))
- y=array(y).T
- return y
-
-def dstep(sys,Tf=10.0):
- """Plot the step response of the discrete system sys
- Only for discrete systems!!!
-
- Call:
- y=dstep(sys, [,Tf=final time]))
-
- Parameters
- ----------
- sys : Discrete System in State Space form
- Tf : Final simulation time
-
- Returns
- -------
- Nothing
-
- """
- Ts=sys.Tsamp
- if Ts==0.0:
- "Only discrete systems allowed!"
- return
-
- ns=int(Tf/Ts+1)
- u=ones((1,ns))
- y=dsimul(sys,u)
- T=arange(0,Tf+Ts/2,Ts)
- plot(T,y)
- grid()
- show()
-
-def dimpulse(sys,Tf=10.0):
- """Plot the impulse response of the discrete system sys
- Only for discrete systems!!!
-
- Call:
- y=dimpulse(sys,[,Tf=final time]))
-
- Parameters
- ----------
- sys : Discrete System in State Space form
- Tf : Final simulation time
-
- Returns
- -------
- Nothing
-
- """
- Ts=sys.Tsamp
- if Ts==0.0:
- "Only discrete systems allowed!"
- return
-
- ns=int(Tf/Ts+1)
- u=zeros((1,ns))
- u[0,0]=1/Ts
- y=dsimul(sys,u)
- T=arange(0,Tf+Ts/2,Ts)
- plot(T,y)
- grid()
- show()
-
-# Step response (plot)
-def bb_step(sys,X0=None,Tf=None,Ts=0.001):
- """Plot the step response of the continous system sys
-
- Call:
- y=bb_step(sys [,Tf=final time] [,Ts=time step])
-
- Parameters
- ----------
- sys : Continous System in State Space form
- X0: Initial state vector (not used yet)
- Ts : sympling time
- Tf : Final simulation time
-
- Returns
- -------
- Nothing
-
- """
- if Tf==None:
- vals = eigvals(sys.A)
- r = min(abs(real(vals)))
- if r == 0.0:
- r = 1.0
- Tf = 7.0 / r
- sysd=c2d(sys,Ts)
- dstep(sysd,Tf=Tf)
-
def sysctr(sys,contr):
"""Build the discrete system controller+plant+output feedback
@@ -449,10 +531,10 @@
as output
"""
- if contr.Tsamp!=sys.Tsamp:
+ if contr.dt!=sys.dt:
print "Systems with different sampling time!!!"
return
- sysf=contr*sys
+ sysf=sys*contr
nu=shape(sysf.B)[1]
b1=mat(sysf.B[:,0])
@@ -470,44 +552,42 @@
Cf=X*mat(sysf.C)
Df=X*d1
- sysc=ss(Af,Bf,Cf,Df,sys.Tsamp)
+ sysc=StateSpace(Af,Bf,Cf,Df,sys.dt)
return sysc
-def care(A,B,Q,R):
- """Solve Riccati equation for discrete time systems
+def set_aw(sys,poles):
+ """Divide in controller in input and feedback part
+ for anti-windup
Usage
=====
- [K, S, E] = care(A, B, Q, R)
+ [sys_in,sys_fbk]=set_aw(sys,poles)
Inputs
------
- A, B: 2-d arrays with dynamics and input matrices
- sys: linear I/O system
- Q, R: 2-d array with state and input weight matrices
+ sys: controller
+ poles : poles for the anti-windup filter
+
Outputs
-------
- X: solution of the Riccati eq.
+ sys_in, sys_fbk: controller in input and feedback part
"""
+# sys=StateSpace(sys);
+ sys=ss(sys);
+ den_old=poly(eigvals(sys.A))
+ den = poly(poles)
+ tmp= tf(den_old,den,sys.dt)
+ tmpss=tf2ss(tmp)
+# sys_in=StateSpace(tmp*sys)
+ sys_in=ss(tmp*sys)
+ sys_in.dt=sys.dt
+# sys_fbk=StateSpace(1-tmp)
+ sys_fbk=ss(1-tmp)
+ sys_fbk.dt=sys.dt
+ return sys_in, sys_fbk
- # Check dimensions for consistency
- nstates = B.shape[0];
- ninputs = B.shape[1];
- if (A.shape[0] != nstates or A.shape[1] != nstates):
- raise ControlDimension("inconsistent system dimensions")
-
- elif (Q.shape[0] != nstates or Q.shape[1] != nstates or
- R.shape[0] != ninputs or R.shape[1] != ninputs) :
- raise ControlDimension("incorrect weighting matrix dimensions")
-
- rcond,X,w,S,T = \
- sb02od(nstates, ninputs, A, B, Q, R, 'C');
-
- return X
-
-
-def dare(A,B,Q,R):
+def bb_dare(A,B,Q,R):
"""Solve Riccati equation for discrete time systems
Usage
@@ -535,13 +615,13 @@
R.shape[0] != ninputs or R.shape[1] != ninputs) :
raise ControlDimension("incorrect weighting matrix dimensions")
- rcond,X,w,S,T = \
+ X,rcond,w,S,T = \
sb02od(nstates, ninputs, A, B, Q, R, 'D');
return X
-def dlqr(*args, **keywords):
+def bb_dlqr(*args, **keywords):
"""Linear quadratic regulator design for discrete systems
Usage
@@ -581,7 +661,7 @@
A = array(args[0].A, ndmin=2, dtype=float);
B = array(args[0].B, ndmin=2, dtype=float);
index = 1;
- if args[0].Tsamp==0.0:
+ if args[0].dt==0.0:
print "dlqr works only for discrete systems!"
return
else:
@@ -595,8 +675,10 @@
R = array(args[index+1], ndmin=2, dtype=float);
if (len(args) > index + 2):
N = array(args[index+2], ndmin=2, dtype=float);
+ Nflag = 1;
else:
N = zeros((Q.shape[0], R.shape[1]));
+ Nflag = 0;
# Check dimensions for consistency
nstates = B.shape[0];
@@ -609,46 +691,154 @@
N.shape[0] != nstates or N.shape[1] != ninputs):
raise ControlDimension("incorrect weighting matrix dimensions")
+ if Nflag==1:
+ Ao=A-B*inv(R)*N.T
+ Qo=Q-N*inv(R)*N.T
+ else:
+ Ao=A
+ Qo=Q
+
#Solve the riccati equation
- X = dare(A,B,Q,R)
+ # (X,L,G) = dare(Ao,B,Qo,R)
+ X = bb_dare(Ao,B,Qo,R)
# Now compute the return value
Phi=mat(A)
H=mat(B)
- K=inv(H.T*X*H+R)*(H.T*X*Phi)
+ K=inv(H.T*X*H+R)*(H.T*X*Phi+N.T)
L=eig(Phi-H*K)
return K,X,L
-def minreal(sys):
- """Minimal representation for state space systems
+def dlqr(*args, **keywords):
+ """Linear quadratic regulator design
+
+ The lqr() function computes the optimal state feedback controller
+ that minimizes the quadratic cost
+
+ .. math:: J = \int_0^\infty x' Q x + u' R u + 2 x' N u
+
+ The function can be called with either 3, 4, or 5 arguments:
+
+ * ``lqr(sys, Q, R)``
+ * ``lqr(sys, Q, R, N)``
+ * ``lqr(A, B, Q, R)``
+ * ``lqr(A, B, Q, R, N)``
+
+ Parameters
+ ----------
+ A, B: 2-d array
+ Dynamics and input matrices
+ sys: Lti (StateSpace or TransferFunction)
+ Linear I/O system
+ Q, R: 2-d array
+ State and input weight matrices
+ N: 2-d array, optional
+ Cross weight matrix
+
+ Returns
+ -------
+ K: 2-d array
+ State feedback gains
+ S: 2-d array
+ Solution to Riccati equation
+ E: 1-d array
+ Eigenvalues of the closed loop system
+
+ Examples
+ --------
+ >>> K, S, E = lqr(sys, Q, R, [N])
+ >>> K, S, E = lqr(A, B, Q, R, [N])
+
+ """
+
+ # Make sure that SLICOT is installed
+ try:
+ from slycot import sb02md
+ from slycot import sb02mt
+ except ImportError:
+ raise ControlSlycot("can't find slycot module 'sb02md' or 'sb02nt'")
+
+ #
+ # Process the arguments and figure out what inputs we received
+ #
+
+ # Get the system description
+ if (len(args) < 4):
+ raise ControlArgument("not enough input arguments")
+
+ elif (ctrlutil.issys(args[0])):
+ # We were passed a system as the first argument; extract A and B
+ #! TODO: really just need to check for A and B attributes
+ A = np.array(args[0].A, ndmin=2, dtype=float);
+ B = np.array(args[0].B, ndmin=2, dtype=float);
+ index = 1;
+ else:
+ # Arguments should be A and B matrices
+ A = np.array(args[0], ndmin=2, dtype=float);
+ B = np.array(args[1], ndmin=2, dtype=float);
+ index = 2;
+
+ # Get the weighting matrices (converting to matrices, if needed)
+ Q = np.array(args[index], ndmin=2, dtype=float);
+ R = np.array(args[index+1], ndmin=2, dtype=float);
+ if (len(args) > index + 2):
+ N = np.array(args[index+2], ndmin=2, dtype=float);
+ else:
+ N = np.zeros((Q.shape[0], R.shape[1]));
+
+ # Check dimensions for consistency
+ nstates = B.shape[0];
+ ninputs = B.shape[1];
+ if (A.shape[0] != nstates or A.shape[1] != nstates):
+ raise ControlDimension("inconsistent system dimensions")
+
+ elif (Q.shape[0] != nstates or Q.shape[1] != nstates or
+ R.shape[0] != ninputs or R.shape[1] != ninputs or
+ N.shape[0] != nstates or N.shape[1] != ninputs):
+ raise ControlDimension("incorrect weighting matrix dimensions")
+
+ # Compute the G matrix required by SB02MD
+ A_b,B_b,Q_b,R_b,L_b,ipiv,oufact,G = \
+ sb02mt(nstates, ninputs, B, R, A, Q, N, jobl='N');
+
+ # Call the SLICOT function
+ X,rcond,w,S,U,A_inv = sb02md(nstates, A_b, G, Q_b, 'D')
+
+ # Now compute the return value
+ K = np.dot(np.linalg.inv(R), (np.dot(B.T, X) + N.T));
+ S = X;
+ E = w[0:nstates];
+
+ return K, S, E
+
+def bb_dcgain(sys):
+ """Return the steady state value of the step response os sys
+
Usage
=====
- [sysmin]=minreal[sys]
+ dcgain=dcgain(sys)
Inputs
------
- sys: system in ss or tf form
+ sys: system
Outputs
-------
- sysfin: system in state space form
+ dcgain : steady state value
"""
+
a=mat(sys.A)
b=mat(sys.B)
c=mat(sys.C)
d=mat(sys.D)
nx=shape(a)[0]
- ni=shape(b)[1]
- no=shape(c)[0]
-
- out=tb03ad(nx,no,ni,a,b,c,d,'R')
-
- nr=out[3]
- A=out[0][:nr,:nr]
- B=out[1][:nr,:ni]
- C=out[2][:no,:nr]
- sysf=ss(A,B,C,sys.D,sys.Tsamp)
- return sysf
-
+ if sys.dt!=0.0:
+ a=a-eye(nx,nx)
+ r=rank(a)
+ if r<nx:
+ gm=[]
+ else:
+ gm=-c*inv(a)*b+d
+ return array(gm)
Modified: trunk/src/matlab.py
===================================================================
--- trunk/src/matlab.py 2014-03-22 18:30:41 UTC (rev 295)
+++ trunk/src/matlab.py 2014-03-22 19:31:36 UTC (rev 296)
@@ -398,7 +398,7 @@
"""
Create a state space system.
- The function accepts either 1 or 4 parameters:
+ The function accepts either 1, 4 or 5 parameters:
``ss(sys)``
Convert a linear system into space system form. Always creates a
@@ -412,7 +412,16 @@
\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.
@@ -429,6 +438,8 @@
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
-------
@@ -457,8 +468,8 @@
"""
- if len(args) == 4:
- return StateSpace(args[0], args[1], args[2], args[3])
+ if len(args) == 4 or len(args) == 5:
+ return StateSpace(*args)
elif len(args) == 1:
sys = args[0]
if isinstance(sys, StateSpace):
@@ -635,7 +646,7 @@
C: array_like or string
Output matrix
D: array_like or string
- Feed forward matrix
+ Feedthrough matrix
Returns
-------
|
|
From: <re...@us...> - 2014-05-23 12:41:56
|
Revision: 299
http://sourceforge.net/p/python-control/code/299
Author: repa
Date: 2014-05-23 12:41:47 +0000 (Fri, 23 May 2014)
Log Message:
-----------
Added c2d functionality for MIMO state-space systems; both in matlab mode
and for python mode; added tests for same
Modified Paths:
--------------
trunk/src/dtime.py
trunk/src/matlab.py
trunk/src/statesp.py
trunk/tests/discrete_test.py
trunk/tests/matlab_test.py
Modified: trunk/src/dtime.py
===================================================================
--- trunk/src/dtime.py 2014-03-23 20:39:48 UTC (rev 298)
+++ trunk/src/dtime.py 2014-05-23 12:41:47 UTC (rev 299)
@@ -94,11 +94,25 @@
if not isctime(sysc):
raise ValueError("First argument must be continuous time system")
- # TODO: impelement MIMO version
+ # If we are passed a state space system, convert to transfer function first
+ if isinstance(sysc, StateSpace) and method == 'zoh':
+
+ try:
+ # try with slycot routine
+ from slycot import mb05nd
+ F, H = mb05nd(sysc.A, Ts)
+ return StateSpace(F, H*sysc.B, sysc.C, sysc.D, Ts)
+ except ImportError:
+ if sysc.inputs != 1 or sysc.outputs != 1:
+ raise TypeError(
+ "mb05nd not found in slycot, or slycot not installed")
+
+ # TODO: implement MIMO version for other than ZOH state-space
if (sysc.inputs != 1 or sysc.outputs != 1):
raise NotImplementedError("MIMO implementation not available")
- # If we are passed a state space system, convert to transfer function first
+ # SISO state-space, with other than ZOH, or failing slycot import,
+ # is handled by conversion to TF
if isinstance(sysc, StateSpace):
warn("sample_system: converting to transfer function")
sysc = _convertToTransferFunction(sysc)
Modified: trunk/src/matlab.py
===================================================================
--- trunk/src/matlab.py 2014-03-23 20:39:48 UTC (rev 298)
+++ trunk/src/matlab.py 2014-05-23 12:41:47 UTC (rev 299)
@@ -1524,8 +1524,29 @@
return (tf.num, tf.den)
# Convert a continuous time system to a discrete time system
-def c2d(sysc, Ts, method):
- # TODO: add docstring
+def c2d(sysc, Ts, method='zoh'):
+ '''
+ Return a discrete-time system
+
+ Parameters
+ ----------
+ sysc: Lti (StateSpace or TransferFunction), continuous
+ System to be converted
+
+ Ts: number
+ Sample time for the conversion
+
+ method: string, optional
+ Method to be applied,
+ 'zoh' Zero-order hold on the inputs (default)
+ 'foh' First-order hold, currently not implemented
+ 'impulse' Impulse-invariant discretization, currently not implemented
+ 'tustin' Bilinear (Tustin) approximation, only SISO
+ 'matched' Matched pole-zero method, only SISO
+ '''
# Call the sample_system() function to do the work
- return sample_system(sysc, Ts, method)
+ sysd = sample_system(sysc, Ts, method)
+ if isinstance(sysc, StateSpace) and not isinstance(sysd, StateSpace):
+ return _convertToStateSpace(sysd)
+ return sysd
Modified: trunk/src/statesp.py
===================================================================
--- trunk/src/statesp.py 2014-03-23 20:39:48 UTC (rev 298)
+++ trunk/src/statesp.py 2014-05-23 12:41:47 UTC (rev 299)
@@ -610,7 +610,7 @@
ssout[3][:sys.outputs, :states],
ssout[4], sys.dt)
except ImportError:
- # TODO: do we want to squeeze first and check dimenations?
+ # TODO: do we want to squeeze first and check dimensions?
# 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))
Modified: trunk/tests/discrete_test.py
===================================================================
--- trunk/tests/discrete_test.py 2014-03-23 20:39:48 UTC (rev 298)
+++ trunk/tests/discrete_test.py 2014-05-23 12:41:47 UTC (rev 299)
@@ -272,6 +272,10 @@
self.assertEqual(sysd.dt, 1)
# TODO: put in other generic checks
+ for sysc in (self.mimo_ss1, self.mimo_ss1c):
+ sysd = sample_system(sysc, 1, method='zoh')
+ self.assertEqual(sysd.dt, 1)
+
# TODO: check results of converstion
# Check errors
Modified: trunk/tests/matlab_test.py
===================================================================
--- trunk/tests/matlab_test.py 2014-03-23 20:39:48 UTC (rev 298)
+++ trunk/tests/matlab_test.py 2014-05-23 12:41:47 UTC (rev 299)
@@ -515,7 +515,26 @@
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])
+ def testSS2cont(self):
+ sys = ss(
+ np.mat("-3 4 2; -1 -3 0; 2 5 3"),
+ np.mat("1 4 ; -3 -3; -2 1"),
+ np.mat("4 2 -3; 1 4 3"),
+ np.mat("-2 4; 0 1"))
+ sysd = c2d(sys, 0.1)
+ np.testing.assert_array_almost_equal(
+ np.mat(
+ """0.742840837331905 0.342242024293711 0.203124211149560;
+ -0.074130792143890 0.724553295044645 -0.009143771143630;
+ 0.180264783290485 0.544385612448419 1.370501013067845"""),
+ sysd.A)
+ np.testing.assert_array_almost_equal(
+ np.mat(""" 0.012362066084719 0.301932197918268;
+ -0.260952977031384 -0.274201791021713;
+ -0.304617775734327 0.075182622718853"""), sysd.B)
+
+
#! TODO: not yet implemented
# def testMIMOtfdata(self):
# sisotf = ss2tf(self.siso_ss1)
|
|
From: <re...@us...> - 2014-07-08 10:57:32
|
Revision: 304
http://sourceforge.net/p/python-control/code/304
Author: repa
Date: 2014-07-08 10:57:18 +0000 (Tue, 08 Jul 2014)
Log Message:
-----------
Correction to the gain margin calculation. Old implementation also returned
gains where arg(H) == 0.
Switched around the w_180 and wc return parameters, to match the order
of gain margin (matching w_180) and phase margin (matching wc) return.
Modified Paths:
--------------
trunk/src/margins.py
trunk/tests/margin_test.py
Modified: trunk/src/margins.py
===================================================================
--- trunk/src/margins.py 2014-07-06 17:06:44 UTC (rev 303)
+++ trunk/src/margins.py 2014-07-08 10:57:18 UTC (rev 304)
@@ -80,6 +80,9 @@
# idea for the frequency data solution copied/adapted from
# https://github.com/alchemyst/Skogestad-Python/blob/master/BODE.py
# Rene van Paassen <ren...@gm...>
+#
+# RvP, July 8, 2014, corrected to exclude phase=0 crossing for the gain
+# margin polynomial
def stability_margins(sysdata, deg=True, returnall=False, epsw=1e-10):
"""Calculate gain, phase and stability margins and associated
crossover frequencies.
@@ -147,7 +150,19 @@
# test imaginary part of tf == 0, for phase crossover/gain margins
test_w_180 = np.polyadd(np.polymul(inum, rden), np.polymul(rnum, -iden))
w_180 = np.roots(test_w_180)
- w_180 = np.real(w_180[(np.imag(w_180) == 0) * (w_180 > epsw)])
+
+ # first remove imaginary and negative frequencies, epsw removes the
+ # "0" frequency for type-2 systems
+ w_180 = np.real(w_180[(np.imag(w_180) == 0) * (w_180 >= epsw)])
+
+ # evaluate response at remaining frequencies, to test for phase 180 vs 0
+ resp_w_180 = np.real(np.polyval(sys.num[0][0], 1.j*w_180) /
+ np.polyval(sys.den[0][0], 1.j*w_180))
+
+ # only keep frequencies where the negative real axis is crossed
+ w_180 = w_180[(resp_w_180 < 0.0)]
+
+ # and sort
w_180.sort()
# test magnitude is 1 for gain crossover/phase margins
@@ -203,14 +218,14 @@
SM = np.abs(sys.evalfr(wstab)[0][0]+1)
if returnall:
- return GM, PM, SM, wc, w_180, wstab
+ return GM, PM, SM, w_180, wc, wstab
else:
return (
(GM.shape[0] or None) and GM[0],
(PM.shape[0] or None) and PM[0],
(SM.shape[0] or None) and SM[0],
+ (w_180.shape[0] or None) and w_180[0],
(wc.shape[0] or None) and wc[0],
- (w_180.shape[0] or None) and w_180[0],
(wstab.shape[0] or None) and wstab[0])
Modified: trunk/tests/margin_test.py
===================================================================
--- trunk/tests/margin_test.py 2014-07-06 17:06:44 UTC (rev 303)
+++ trunk/tests/margin_test.py 2014-07-08 10:57:18 UTC (rev 304)
@@ -17,11 +17,18 @@
self.sys2 = TransferFunction([1], [1, 2, 3, 4])
self.sys3 = StateSpace([[1., 4.], [3., 2.]], [[1.], [-4.]],
[[1., 0.]], [[0.]])
+ s = TransferFunction([1, 0], [1])
+ self.sys4 = (8.75*(4*s**2+0.4*s+1))/((100*s+1)*(s**2+0.22*s+1)) * \
+ 1./(s**2/(10.**2)+2*0.04*s/10.+1)
def test_stability_margins(self):
gm, pm, sm, wg, wp, ws = stability_margins(self.sys1);
gm, pm, sm, wg, wp, ws = stability_margins(self.sys2);
gm, pm, sm, wg, wp, ws = stability_margins(self.sys3);
+ gm, pm, sm, wg, wp, ws = stability_margins(self.sys4);
+ np.testing.assert_array_almost_equal(
+ [gm, pm, sm, wg, wp, ws],
+ [2.2716, 97.5941, 1.0454, 10.0053, 0.0850, 0.4973], 3)
def test_phase_crossover_frequencies(self):
omega, gain = phase_crossover_frequencies(self.sys2)
|
|
From: <re...@us...> - 2014-07-08 10:59:08
|
Revision: 305
http://sourceforge.net/p/python-control/code/305
Author: repa
Date: 2014-07-08 10:59:01 +0000 (Tue, 08 Jul 2014)
Log Message:
-----------
corrected comments on time responses, with the "input" parameter to be
ignored on initial response calculation; initial response does not depend
on input!
Modified Paths:
--------------
trunk/src/matlab.py
trunk/src/timeresp.py
trunk/tests/matlab_test.py
trunk/tests/timeresp_test.py
Modified: trunk/src/matlab.py
===================================================================
--- trunk/src/matlab.py 2014-07-08 10:57:18 UTC (rev 304)
+++ trunk/src/matlab.py 2014-07-08 10:59:01 UTC (rev 305)
@@ -1103,8 +1103,8 @@
"""Root locus plot
The root-locus plot has a callback function that prints pole location,
- gain and damping to the Python consol on mouseclicks on the root-locus
- graph.
+ gain and damping to the Python console on mouseclicks on the root-locus
+ graph.
Parameters
----------
@@ -1112,6 +1112,9 @@
Linear system
klist:
optional list of gains
+
+ Keyword parameters
+ ------------------
xlim : control of x-axis range, normally with tuple, for
other options, see matplotlib.axes
ylim : control of y-axis range
@@ -1165,7 +1168,7 @@
margin: no magnitude crossings found
.. todo::
- better ecample system!
+ better example system!
#>>> gm, pm, wg, wp = margin(mag, phase, w)
"""
@@ -1178,7 +1181,7 @@
raise ValueError("Margin needs 1 or 3 arguments; received %i."
% len(args))
- return margin[0], margin[1], margin[4], margin[3]
+ return margin[0], margin[1], margin[3], margin[4]
def dcgain(*args):
'''
@@ -1279,10 +1282,11 @@
'''
Step response of a linear system
- If the system has multiple inputs or outputs (MIMO), one input and one
- output have to be selected for the simulation. The parameters `input`
- and `output` do this. All other inputs are set to 0, all other outputs
- are ignored.
+ If the system has multiple inputs or outputs (MIMO), one input has
+ to be selected for the simulation. Optionally, one output may be
+ selected. If no selection is made for the output, all outputs are
+ given. The parameters `input` and `output` do this. All other
+ inputs are set to 0, all other outputs are ignored.
Parameters
----------
@@ -1301,7 +1305,7 @@
Index of the input that will be used in this simulation.
output: int
- Index of the output that will be used in this simulation.
+ If given, index of the output that is returned by this simulation.
**keywords:
Additional keyword arguments control the solution algorithm for the
@@ -1326,19 +1330,21 @@
Examples
--------
>>> yout, T = step(sys, T, X0)
+
'''
T, yout = timeresp.step_response(sys, T, X0, input, output,
- transpose = True, **keywords)
+ transpose=True, **keywords)
return yout, T
-def impulse(sys, T=None, input=0, output=0, **keywords):
+def impulse(sys, T=None, input=0, output=None, **keywords):
'''
Impulse response of a linear system
- If the system has multiple inputs or outputs (MIMO), one input and
- one output must be selected for the simulation. The parameters
- `input` and `output` do this. All other inputs are set to 0, all
- other outputs are ignored.
+ If the system has multiple inputs or outputs (MIMO), one input has
+ to be selected for the simulation. Optionally, one output may be
+ selected. If no selection is made for the output, all outputs are
+ given. The parameters `input` and `output` do this. All other
+ inputs are set to 0, all other outputs are ignored.
Parameters
----------
@@ -1381,14 +1387,13 @@
transpose = True, **keywords)
return yout, T
-def initial(sys, T=None, X0=0., input=0, output=0, **keywords):
+def initial(sys, T=None, X0=0., input=None, output=None, **keywords):
'''
Initial condition response of a linear system
- If the system has multiple inputs or outputs (MIMO), one input and one
- output have to be selected for the simulation. The parameters `input`
- and `output` do this. All other inputs are set to 0, all other outputs
- are ignored.
+ If the system has multiple outputs (?IMO), optionally, one output
+ may be selected. If no selection is made for the output, all
+ outputs are given.
Parameters
----------
@@ -1404,10 +1409,11 @@
Numbers are converted to constant arrays with the correct shape.
input: int
- Index of the input that will be used in this simulation.
+ This input is ignored, but present for compatibility with step
+ and impulse.
output: int
- Index of the output that will be used in this simulation.
+ If given, index of the output that is returned by this simulation.
**keywords:
Additional keyword arguments control the solution algorithm for the
@@ -1432,9 +1438,10 @@
Examples
--------
>>> T, yout = initial(sys, T, X0)
+
'''
- T, yout = timeresp.initial_response(sys, T, X0, input, output,
- transpose = True, **keywords)
+ T, yout = timeresp.initial_response(sys, T, X0, output=output,
+ transpose=True, **keywords)
return yout, T
def lsim(sys, U=0., T=None, X0=0., **keywords):
Modified: trunk/src/timeresp.py
===================================================================
--- trunk/src/timeresp.py 2014-07-08 10:57:18 UTC (rev 304)
+++ trunk/src/timeresp.py 2014-07-08 10:59:01 UTC (rev 305)
@@ -389,7 +389,8 @@
def step_response(sys, T=None, X0=0., input=0, output=None,
transpose = False, **keywords):
#pylint: disable=W0622
- """Step response of a linear system
+ """
+ Step response of a linear system
If the system has multiple inputs or outputs (MIMO), one input has
to be selected for the simulation. Optionally, one output may be
@@ -468,15 +469,14 @@
return T, yout
-def initial_response(sys, T=None, X0=0., input=0, output=None, transpose=False,
- **keywords):
+def initial_response(sys, T=None, X0=0., input=None, output=None,
+ transpose=False, **keywords):
#pylint: disable=W0622
"""Initial condition response of a linear system
- If the system has multiple inputs or outputs (MIMO), one input and one
- output have to be selected for the simulation. The parameters `input`
- and `output` do this. All other inputs are set to 0, all other outputs
- are ignored.
+ If the system has multiple outputs (?IMO), optionally, one output
+ may be selected. If no selection is made for the output, all
+ outputs are given.
For information on the **shape** of parameters `T`, `X0` and
return values `T`, `yout` see: :ref:`time-series-convention`
@@ -495,7 +495,8 @@
Numbers are converted to constant arrays with the correct shape.
input: int
- Index of the input that will be used in this simulation.
+ Ignored, has no meaning in initial condition calculation. Parameter
+ ensures compatibility with step_response and impulse_response
output: int
Index of the output that will be used in this simulation. Set to None
@@ -531,9 +532,9 @@
"""
sys = _convertToStateSpace(sys)
if output == None:
- sys = _mimo2simo(sys, input, warn_conversion=False)
+ sys = _mimo2simo(sys, 0, warn_conversion=False)
else:
- sys = _mimo2siso(sys, input, output, warn_conversion=False)
+ sys = _mimo2siso(sys, 0, output, warn_conversion=False)
# Create time and input vectors; checking is done in forced_response(...)
# The initial vector X0 is created in forced_response(...) if necessary
@@ -549,12 +550,13 @@
def impulse_response(sys, T=None, X0=0., input=0, output=None,
transpose=False, **keywords):
#pylint: disable=W0622
- """Impulse response of a linear system
+ """
+ Impulse response of a linear system
- If the system has multiple inputs or outputs (MIMO), one input and one
- output have to be selected for the simulation. The parameters `input`
- and `output` do this. All other inputs are set to 0, all other outputs
- are ignored.
+ If the system has multiple inputs or outputs (MIMO), one input has
+ to be selected for the simulation. Optionally, one output may be
+ selected. The parameters `input` and `output` do this. All other
+ inputs are set to 0, all other outputs are ignored.
For information on the **shape** of parameters `T`, `X0` and
return values `T`, `yout` see: :ref:`time-series-convention`
Modified: trunk/tests/matlab_test.py
===================================================================
--- trunk/tests/matlab_test.py 2014-07-08 10:57:18 UTC (rev 304)
+++ trunk/tests/matlab_test.py 2014-07-08 10:59:01 UTC (rev 305)
@@ -190,8 +190,8 @@
#Test MIMO system, which contains ``siso_ss1`` twice
sys = self.mimo_ss1
x0 = np.matrix(".5; 1.; .5; 1.")
- y_00, _t = initial(sys, T=t, X0=x0, input=0, output=0)
- y_11, _t = initial(sys, T=t, X0=x0, input=1, output=1)
+ y_00, _t = initial(sys, T=t, X0=x0, output=0)
+ y_11, _t = initial(sys, T=t, X0=x0, output=1)
np.testing.assert_array_almost_equal(y_00, youttrue, decimal=4)
np.testing.assert_array_almost_equal(y_11, youttrue, decimal=4)
Modified: trunk/tests/timeresp_test.py
===================================================================
--- trunk/tests/timeresp_test.py 2014-07-08 10:57:18 UTC (rev 304)
+++ trunk/tests/timeresp_test.py 2014-07-08 10:59:01 UTC (rev 305)
@@ -109,8 +109,8 @@
#Test MIMO system, which contains ``siso_ss1`` twice
sys = self.mimo_ss1
x0 = np.matrix(".5; 1.; .5; 1.")
- _t, y_00 = initial_response(sys, T=t, X0=x0, input=0, output=0)
- _t, y_11 = initial_response(sys, T=t, X0=x0, input=1, output=1)
+ _t, y_00 = initial_response(sys, T=t, X0=x0, output=0)
+ _t, y_11 = initial_response(sys, T=t, X0=x0, output=1)
np.testing.assert_array_almost_equal(y_00, youttrue, decimal=4)
np.testing.assert_array_almost_equal(y_11, youttrue, decimal=4)
|
×
Want the latest updates on software, tech news, and AI?
Get latest updates about software, tech news, and AI from SourceForge directly in your inbox once a month.
Thanks for helping keep SourceForge clean.
X