Menu
▾
▴
[python-control-discuss] SF.net SVN: python-control:[161] trunk
|
From: <mur...@us...> - 2011-06-18 19:13:17
|
Revision: 161
http://python-control.svn.sourceforge.net/python-control/?rev=161&view=rev
Author: murrayrm
Date: 2011-06-18 19:13:09 +0000 (Sat, 18 Jun 2011)
Log Message:
-----------
Added Eike's Welk's time response functions, including MIMO
capability. To maintain MATLAB compatibility, the underlying
functions have been moved to a new module, timeresp, and renamed
ForcedResponse (lsim), InitialResponse (initial), ImpulseResponse
(impulse), StepResponse (step). MATLAB version of the functions (with
the usual names: lsim, impulse, initial, step) use MATLAB conventions
for time vectors.
See ChangeLog for detailed list of changes.
Modified Paths:
--------------
trunk/ChangeLog
trunk/doc/conf.py
trunk/examples/pvtol-lqr.py
trunk/examples/pvtol-nested-ss.py
trunk/examples/pvtol-nested.py
trunk/src/__init__.py
trunk/src/matlab.py
trunk/src/statesp.py
trunk/tests/matlab_test.py
Added Paths:
-----------
trunk/src/timeresp.py
trunk/tests/test_control_matlab.py
trunk/tests/timeresp_test.py
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2011-06-18 00:20:49 UTC (rev 160)
+++ trunk/ChangeLog 2011-06-18 19:13:09 UTC (rev 161)
@@ -1,3 +1,66 @@
+2011-06-18 Richard Murray <mu...@ma...>
+
+ * src/timeresp.py, src/matlab.py: moved documentation about time
+ series convention from matlab.py to timeresp.py
+
+ * examples/pvtol-nested-ss.py: Fixed bug in call to step (wrong
+ second argument)
+
+ * tests/matlab_test.py: Updated tests to use MATLAB time response
+ conventions.
+
+ * tests/timeresp_test.py: Created unit tests for timeresp module,
+ based on matlab_test.py
+
+2011-06-17 Richard Murray <mu...@ma...>
+
+ * src/timeresp.py (ForcedResponse): swapped order of input and time
+ arguments for linear response, following Eike's comment "T must
+ always be supplied by the user, but U has a useful default value of
+ 0."
+
+ * src/matlab.py: moved code for lsim, initial, step, and impulse to
+ timeresp.py and put in new routes that call timeresp.* versions of
+ the functions with transposeData set to True.
+
+ * src/timesim.py (_check_convert_array): added transpose argument
+ that will transpose input data before processing it.
+
+ * src/timesim.py: renamed lsim, initial, step, and impulse functions
+ to ForcedResponse, InitialResponse, StepResponse and
+ ImpulseResponse. These versions use Eike Welk's input ordering.
+
+ * examples/pvtol-nested.py: calls to step() had screwed up inputs.
+ Fixed.
+
+2011-06-17 Richard Murray <mu...@ma...>
+
+ * src/matlab.py: added MIMO extensions from Eike Welk on 12 Jun
+ 2011: adds MIMO capabilities for ``lsim``, ``step``, ``impulse``,
+ ``initial``
+
+ * src/matlab.py: added changes from Eike Welk on 12 May 2011:
+
+ - An implementation of the four simulation functions ``lsim``,
+ ``step``, ``initial``, and ``impulse`` of the module ``matlab``.
+
+ - Adds a function ``dcgain`` to the ``matlab`` module, which
+ computes the gain of a linear system for steady state and
+ constant input.
+
+ - The patch contains a bug fix for class ``StateSpace``, which
+ enables it to work properly together with Scipy's ``signal``
+ module.
+
+ - The simulation functions' return values are changed (back?) to
+ arrays, because matrices confuse Matplotlib.
+
+ - New times series convention: see _time-series-convention section
+ of matlab documentation
+
+ - SISO simulation data are squeezed on output. To turn this off,
+ pass the option squeeze=False
+
---- control-0.4c released -----
2011-06-17 Richard Murray <mu...@dh...>
Modified: trunk/doc/conf.py
===================================================================
--- trunk/doc/conf.py 2011-06-18 00:20:49 UTC (rev 160)
+++ trunk/doc/conf.py 2011-06-18 19:13:09 UTC (rev 161)
@@ -26,7 +26,8 @@
# Add any Sphinx extension module names here, as strings. They can be extensions
# coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
-extensions = ['sphinx.ext.autodoc','numpydoc']
+extensions = ['sphinx.ext.autodoc', 'numpydoc', 'sphinx.ext.pngmath',
+ 'sphinx.ext.intersphinx']
# Add any paths that contain templates here, relative to this directory.
templates_path = ['_templates']
@@ -87,7 +88,11 @@
# A list of ignored prefixes for module index sorting.
#modindex_common_prefix = []
+#This config value contains the locations and names of other projects that
+#should be linked to in this documentation.
+intersphinx_mapping = {'scipy':('http://docs.scipy.org/doc/scipy/reference/', None)}
+
# -- Options for HTML output ---------------------------------------------------
# The theme to use for HTML and HTML Help pages. See the documentation for
Modified: trunk/examples/pvtol-lqr.py
===================================================================
--- trunk/examples/pvtol-lqr.py 2011-06-18 00:20:49 UTC (rev 160)
+++ trunk/examples/pvtol-lqr.py 2011-06-18 19:13:09 UTC (rev 161)
@@ -119,7 +119,7 @@
subplot(221); title("Identity weights")
# plot(T, Y[:,1, 1], '-', T, Y[:,2, 2], '--'); hold(True);
-plot(Tx.T, Yx[0,:].T, '-', Ty.T, Yy[0,:].T, '--'); hold(True);
+plot(Tx.T, Yx.T, '-', Ty.T, Yy.T, '--'); hold(True);
plot([0, 10], [1, 1], 'k-'); hold(True);
axis([0, 10, -0.1, 1.4]);
@@ -141,9 +141,9 @@
[T3, Y3] = step(H1cx, T=linspace(0,10,100));
subplot(222); title("Effect of input weights")
-plot(T1.T, Y1[0,:].T, 'b-'); hold(True);
-plot(T2.T, Y2[0,:].T, 'b-'); hold(True);
-plot(T3.T, Y3[0,:].T, 'b-'); hold(True);
+plot(T1.T, Y1.T, 'b-'); hold(True);
+plot(T2.T, Y2.T, 'b-'); hold(True);
+plot(T3.T, Y3.T, 'b-'); hold(True);
plot([0 ,10], [1, 1], 'k-'); hold(True);
axis([0, 10, -0.1, 1.4]);
@@ -162,7 +162,7 @@
subplot(223); title("Output weighting")
[T2x, Y2x] = step(H2x, T=linspace(0,10,100));
[T2y, Y2y] = step(H2y, T=linspace(0,10,100));
-plot(T2x.T, Y2x[0,:].T, T2y.T, Y2y[0,:].T)
+plot(T2x.T, Y2x.T, T2y.T, Y2y.T)
ylabel('position');
xlabel('time'); ylabel('position');
legend(('x', 'y'), loc='lower right');
@@ -185,7 +185,7 @@
# step(H3x, H3y, 10);
[T3x, Y3x] = step(H3x, T=linspace(0,10,100));
[T3y, Y3y] = step(H3y, T=linspace(0,10,100));
-plot(T3x.T, Y3x[0,:].T, T3y.T, Y3y[0,:].T)
+plot(T3x.T, Y3x.T, T3y.T, Y3y.T)
title("Physically motivated weights")
xlabel('time');
legend(('x', 'y'), loc='lower right');
Modified: trunk/examples/pvtol-nested-ss.py
===================================================================
--- trunk/examples/pvtol-nested-ss.py 2011-06-18 00:20:49 UTC (rev 160)
+++ trunk/examples/pvtol-nested-ss.py 2011-06-18 19:13:09 UTC (rev 161)
@@ -144,10 +144,10 @@
# 'EdgeColor', color, 'FaceColor', color);
figure(9);
-(Tvec, Yvec) = step(T, None, linspace(1, 20));
+(Tvec, Yvec) = step(T, linspace(1, 20));
plot(Tvec.T, Yvec.T); hold(True);
-(Tvec, Yvec) = step(Co*S, None, linspace(1, 20));
+(Tvec, Yvec) = step(Co*S, linspace(1, 20));
plot(Tvec.T, Yvec.T);
#TODO: PZmap for statespace systems has not yet been implemented.
Modified: trunk/examples/pvtol-nested.py
===================================================================
--- trunk/examples/pvtol-nested.py 2011-06-18 00:20:49 UTC (rev 160)
+++ trunk/examples/pvtol-nested.py 2011-06-18 19:13:09 UTC (rev 161)
@@ -134,10 +134,10 @@
# 'EdgeColor', color, 'FaceColor', color);
figure(9);
-(Tvec, Yvec) = step(T, None, linspace(1, 20));
+(Tvec, Yvec) = step(T, linspace(1, 20));
plot(Tvec.T, Yvec.T); hold(True);
-(Tvec, Yvec) = step(Co*S, None, linspace(1, 20));
+(Tvec, Yvec) = step(Co*S, linspace(1, 20));
plot(Tvec.T, Yvec.T);
figure(10); clf();
Modified: trunk/src/__init__.py
===================================================================
--- trunk/src/__init__.py 2011-06-18 00:20:49 UTC (rev 160)
+++ trunk/src/__init__.py 2011-06-18 19:13:09 UTC (rev 161)
@@ -66,3 +66,5 @@
from modelsimp import *
from rlocus import *
from mateqn import *
+from timeresp import ForcedResponse, InitialResponse, StepResponse, \
+ ImpulseResponse
Modified: trunk/src/matlab.py
===================================================================
--- trunk/src/matlab.py 2011-06-18 00:20:49 UTC (rev 160)
+++ trunk/src/matlab.py 2011-06-18 19:13:09 UTC (rev 161)
@@ -1,3 +1,4 @@
+# -*- coding: utf-8 -*-
"""matlab.py
MATLAB emulation functions.
@@ -18,6 +19,8 @@
"""Copyright (c) 2009 by California Institute of Technology
All rights reserved.
+Copyright (c) 2011 by Eike Welk
+
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:
@@ -58,8 +61,9 @@
# Libraries that we make use of
import scipy as sp # SciPy library (used all over)
import numpy as np # NumPy library
-import scipy.signal as signal # Signal processing library
+from scipy.signal.ltisys import _default_response_times
from copy import deepcopy
+import warnings
# Import MATLAB-like functions that are defined in other packages
from scipy.signal import zpk2ss, ss2zpk, tf2zpk, zpk2tf
@@ -68,8 +72,10 @@
# Control system library
import ctrlutil
import freqplot
+import timeresp
from statesp import StateSpace, _rss_generate, _convertToStateSpace
from xferfcn import TransferFunction, _convertToTransferFunction
+from lti import Lti #base class of StateSpace, TransferFunction
from exception import ControlArgument
# Import MATLAB-like functions that can be used as-is
@@ -86,188 +92,216 @@
__doc__ += """
The control.matlab module defines functions that are roughly the
equivalents of those in the MATLAB Control Toolbox. Items marked by a
-\* are currently implemented; those marked with a ``-`` are not planned
+``*`` are currently implemented; those marked with a ``-`` are not planned
for implementation.
-Creating linear models.
-\* tf - create transfer function (TF) models
- zpk - create zero/pole/gain (ZPK) models.
-\* 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
- lti/exp - create pure continuous-time delays (TF and ZPK only)
- filt - specify digital filters
-- lti/set - set/modify properties of LTI models
-- setdelaymodel - specify internal delay model (state space only)
-
-Data extraction
- lti/tfdata - extract numerators and denominators
- lti/zpkdata - extract zero/pole/gain data
- lti/ssdata - extract state-space matrices
- lti/dssdata - descriptor version of SSDATA
- frd/frdata - extract frequency response data
- lti/get - access values of LTI model properties
- ss/getDelayModel - access internal delay model (state space only)
-
-Conversions
-\* tf - conversion to transfer function
- zpk - conversion to zero/pole/gain
-\* ss - conversion to state space
- frd - conversion to frequency data
- c2d - continuous to discrete conversion
- d2c - discrete to continuous conversion
- d2d - resample discrete-time model
- upsample - upsample discrete-time LTI systems
-\* ss2tf - state space to transfer function
- ss2zpk - transfer function to zero-pole-gain
-\* tf2ss - transfer function to state space
- tf2zpk - transfer function to zero-pole-gain
- zpk2ss - zero-pole-gain to state space
- zpk2tf - zero-pole-gain to transfer function
-
-System interconnections
- append - group LTI models by appending inputs and outputs
-\* parallel - connect LTI models in parallel (see also overloaded +)
-\* series - connect LTI models in series (see also overloaded \*)
-\* feedback - connect lti models with a feedback loop
- lti/lft - generalized feedback interconnection
- lti/connect - arbitrary interconnection of lti models
- sumblk - specify summing junction (for use with connect)
- strseq - builds sequence of indexed strings (for I/O naming)
-
-System gain and dynamics
- dcgain - steady-state (D.C.) gain
- lti/bandwidth - system bandwidth
- lti/norm - h2 and Hinfinity norms of LTI models
-\* lti/pole - system poles
-\* lti/zero - system (transmission) zeros
- lti/order - model order (number of states)
-\* pzmap - pole-zero map (TF only)
- lti/iopzmap - input/output pole-zero map
- damp - natural frequency and damping of system poles
- esort - sort continuous poles by real part
- dsort - sort discrete poles by magnitude
- lti/stabsep - stable/unstable decomposition
- lti/modsep - region-based modal decomposition
-
-Time-domain analysis
-\* step - step response
- stepinfo - step response characteristics (rise time, ...)
-\* impulse - impulse response
- initial - free response with initial conditions
-\* lsim - response to user-defined input signal
- lsiminfo - linear response characteristics
- gensig - generate input signal for LSIM
- covar - covariance of response to white noise
-
-Frequency-domain analysis
-\* bode - Bode plot of the frequency response
- lti/bodemag - Bode magnitude diagram only
- sigma - singular value frequency plot
-\* nyquist - Nyquist plot
-\* nichols - Nichols plot
- margin - gain and phase margins
- lti/allmargin - all crossover frequencies and related gain/phase margins
-\* lti/freqresp - frequency response over a frequency grid
-\* lti/evalfr - evaluate frequency response at given frequency
-
-Model simplification
- minreal - minimal realization and pole/zero cancellation
- ss/sminreal - structurally minimal realization (state space)
-\* lti/hsvd - hankel singular values (state contributions)
-\* lti/balred - reduced-order approximations of LTI models
-\* ss/modred - model order reduction
-
-Compensator design
-\* rlocus - evans root locus
-\* place - pole placement
- estim - form estimator given estimator gain
- reg - form regulator given state-feedback and estimator gains
-
-LQR/LQG design
- ss/lqg - single-step LQG design
-\* lqr - linear-Quadratic (LQ) state-feedback regulator
- dlqr - discrete-time LQ state-feedback regulator
- lqry - lq regulator with output weighting
- lqrd - discrete LQ regulator for continuous plant
- ss/lqi - linear-Quadratic-Integral (LQI) controller
- ss/kalman - Kalman state estimator
- ss/kalmd - discrete Kalman estimator for continuous plant
- ss/lqgreg - build LQG regulator from LQ gain and Kalman estimator
- ss/lqgtrack - build LQG servo-controller
- augstate - augment output by appending states
-
-State-space (SS) models
-\* rss - random stable continuous-time state-space models
-\* drss - random stable discrete-time state-space models
- ss2ss - state coordinate transformation
- canon - canonical forms of state-space models
-\* ctrb - controllability matrix
-\* obsv - observability matrix
-\* gram - controllability and observability gramians
- ss/prescale - optimal scaling of state-space models.
- balreal - gramian-based input/output balancing
- ss/xperm - reorder states.
-
-Frequency response data (FRD) models
- frd/chgunits - change frequency vector units
- frd/fcat - merge frequency responses
- frd/fselect - select frequency range or subgrid
- frd/fnorm - peak gain as a function of frequency
- frd/abs - entrywise magnitude of the frequency response
- frd/real - real part of the frequency response
- frd/imag - imaginary part of the frequency response
- frd/interp - interpolate frequency response data
- mag2db - convert magnitude to decibels (dB)
- db2mag - convert decibels (dB) to magnitude
-
-Time delays
- lti/hasdelay - true for models with time delays
- lti/totaldelay - total delay between each input/output pair
- lti/delay2z - replace delays by poles at z=0 or FRD phase shift
-\* pade - pade approximation of time delays
-
-Model dimensions and characteristics
- class - model type ('tf', 'zpk', 'ss', or 'frd')
- isa - test if model is of given type
- tf/size - model sizes
- lti/ndims - number of dimensions
- lti/isempty - true for empty models
- lti/isct - true for continuous-time models
- lti/isdt - true for discrete-time models
- lti/isproper - true for proper models
- lti/issiso - true for single-input/single-output models
- lti/isstable - true for models with stable dynamics
- lti/reshape - reshape array of linear models
-
-Overloaded arithmetic operations
-\* + and - - add and subtract systems (parallel connection)
-\* \* - multiply systems (series connection)
- / - left divide -- sys1\sys2 means inv(sys1)\*sys2
-- \ - right divide -- sys1/sys2 means sys1\*inv(sys2)
- ^ - powers of a given system
- ' - pertransposition
- .' - transposition of input/output map
- .\* - element-by-element multiplication
- [..] - concatenate models along inputs or outputs
- lti/stack - stack models/arrays along some array dimension
- lti/inv - inverse of an LTI system
- lti/conj - complex conjugation of model coefficients
-
-Matrix equation solvers and linear algebra
-\* lyap, dlyap - solve Lyapunov equations
- lyapchol, dlyapchol - square-root Lyapunov solvers
-\* care, dare - solve algebraic Riccati equations
- gcare, gdare - generalized Riccati solvers
- bdschur - block diagonalization of a square matrix
-
-Additional functions
-\* gangof4 - generate the Gang of 4 sensitivity plots
-\* linspace - generate a set of numbers that are linearly spaced
-\* logspace - generate a set of numbers that are logarithmically spaced
-\* unwrap - unwrap a phase angle to give a continuous curve
-
+== ========================== ================================================
+**Creating linear models.**
+--------------------------------------------------------------------------------
+\* :func:`tf` create transfer function (TF) models
+\ zpk create zero/pole/gain (ZPK) models.
+\* :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
+\ lti/exp create pure continuous-time delays (TF and ZPK
+ only)
+\ filt specify digital filters
+\- lti/set set/modify properties of LTI models
+\- setdelaymodel specify internal delay model (state space only)
+\
+**Data extraction**
+--------------------------------------------------------------------------------
+\ lti/tfdata extract numerators and denominators
+\ lti/zpkdata extract zero/pole/gain data
+\ lti/ssdata extract state-space matrices
+\ lti/dssdata descriptor version of SSDATA
+\ frd/frdata extract frequency response data
+\ lti/get access values of LTI model properties
+\ ss/getDelayModel access internal delay model (state space only)
+\
+**Conversions**
+--------------------------------------------------------------------------------
+\* :func:`tf` conversion to transfer function
+\ zpk conversion to zero/pole/gain
+\* :func:`ss` conversion to state space
+\ frd conversion to frequency data
+\ c2d continuous to discrete conversion
+\ d2c discrete to continuous conversion
+\ d2d resample discrete-time model
+\ upsample upsample discrete-time LTI systems
+\* :func:`ss2tf` state space to transfer function
+\ ss2zpk transfer function to zero-pole-gain
+\* :func:`tf2ss` transfer function to state space
+\ tf2zpk transfer function to zero-pole-gain
+\ zpk2ss zero-pole-gain to state space
+\ zpk2tf zero-pole-gain to transfer function
+\
+**System interconnections**
+--------------------------------------------------------------------------------
+\ append group LTI models by appending inputs and outputs
+\* :func:`parallel` connect LTI models in parallel
+ (see also overloaded +)
+\* :func:`series` connect LTI models in series
+ (see also overloaded \*)
+\* :func:`feedback` connect lti models with a feedback loop
+\ lti/lft generalized feedback interconnection
+\ lti/connect arbitrary interconnection of lti models
+\ sumblk specify summing junction (for use with connect)
+\ strseq builds sequence of indexed strings
+ (for I/O naming)
+\
+**System gain and dynamics**
+--------------------------------------------------------------------------------
+\* :func:`dcgain` steady-state (D.C.) gain
+\ lti/bandwidth system bandwidth
+\ lti/norm h2 and Hinfinity norms of LTI models
+\* :func:`pole` system poles
+\* :func:`zero` system (transmission) zeros
+\ lti/order model order (number of states)
+\* :func:`pzmap` pole-zero map (TF only)
+\ lti/iopzmap input/output pole-zero map
+\ damp natural frequency and damping of system poles
+\ esort sort continuous poles by real part
+\ dsort sort discrete poles by magnitude
+\ lti/stabsep stable/unstable decomposition
+\ lti/modsep region-based modal decomposition
+\
+**Time-domain analysis**
+--------------------------------------------------------------------------------
+\* :func:`step` step response
+\ stepinfo step response characteristics (rise time, ...)
+\* :func:`impulse` impulse response
+\* :func:`initial` free response with initial conditions
+\* :func:`lsim` response to user-defined input signal
+\ lsiminfo linear response characteristics
+\ gensig generate input signal for LSIM
+\ covar covariance of response to white noise
+\
+**Frequency-domain analysis**
+--------------------------------------------------------------------------------
+\* :func:`bode` Bode plot of the frequency response
+\ lti/bodemag Bode magnitude diagram only
+\ sigma singular value frequency plot
+\* :func:`nyquist` Nyquist plot
+\* :func:`nichols` Nichols plot
+\* :func:`margin` gain and phase margins
+\ lti/allmargin all crossover frequencies and related gain/phase
+ margins
+\* :func:`freqresp` frequency response over a frequency grid
+\* :func:`evalfr` evaluate frequency response at given frequency
+\
+**Model simplification**
+--------------------------------------------------------------------------------
+\ minreal minimal realization and pole/zero cancellation
+\ ss/sminreal structurally minimal realization (state space)
+\* :func:`lti/hsvd` hankel singular values (state contributions)
+\* :func:`lti/balred` reduced-order approximations of LTI models
+\* :func:`ss/modred` model order reduction
+\
+**Compensator design**
+--------------------------------------------------------------------------------
+\* :func:`rlocus` evans root locus
+\* :func:`place` pole placement
+\ estim form estimator given estimator gain
+\ reg form regulator given state-feedback and
+ estimator gains
+\
+**LQR/LQG design**
+--------------------------------------------------------------------------------
+\ ss/lqg single-step LQG design
+\* :func:`lqr` linear-Quadratic (LQ) state-feedback regulator
+\ dlqr discrete-time LQ state-feedback regulator
+\ lqry lq regulator with output weighting
+\ lqrd discrete LQ regulator for continuous plant
+\ ss/lqi linear-Quadratic-Integral (LQI) controller
+\ ss/kalman Kalman state estimator
+\ ss/kalmd discrete Kalman estimator for continuous plant
+\ ss/lqgreg build LQG regulator from LQ gain and Kalman
+ estimator
+\ ss/lqgtrack build LQG servo-controller
+\ augstate augment output by appending states
+\
+**State-space (SS) models**
+--------------------------------------------------------------------------------
+\* :func:`rss` random stable continuous-time state-space models
+\* :func:`drss` random stable discrete-time state-space models
+\ ss2ss state coordinate transformation
+\ canon canonical forms of state-space models
+\* :func:`ctrb` controllability matrix
+\* :func:`obsv` observability matrix
+\* :func:`gram` controllability and observability gramians
+\ ss/prescale optimal scaling of state-space models.
+\ balreal gramian-based input/output balancing
+\ ss/xperm reorder states.
+\
+**Frequency response data (FRD) models**
+--------------------------------------------------------------------------------
+\ frd/chgunits change frequency vector units
+\ frd/fcat merge frequency responses
+\ frd/fselect select frequency range or subgrid
+\ frd/fnorm peak gain as a function of frequency
+\ frd/abs entrywise magnitude of the frequency response
+\ frd/real real part of the frequency response
+\ frd/imag imaginary part of the frequency response
+\ frd/interp interpolate frequency response data
+\ mag2db convert magnitude to decibels (dB)
+\ db2mag convert decibels (dB) to magnitude
+\
+**Time delays**
+--------------------------------------------------------------------------------
+\ lti/hasdelay true for models with time delays
+\ lti/totaldelay total delay between each input/output pair
+\ lti/delay2z replace delays by poles at z=0 or FRD phase
+ shift
+\* :func:`pade` pade approximation of time delays
+\
+**Model dimensions and characteristics**
+--------------------------------------------------------------------------------
+\ class model type ('tf', 'zpk', 'ss', or 'frd')
+\ isa test if model is of given type
+\ tf/size model sizes
+\ lti/ndims number of dimensions
+\ lti/isempty true for empty models
+\ lti/isct true for continuous-time models
+\ lti/isdt true for discrete-time models
+\ lti/isproper true for proper models
+\ lti/issiso true for single-input/single-output models
+\ lti/isstable true for models with stable dynamics
+\ lti/reshape reshape array of linear models
+\
+**Overloaded arithmetic operations**
+--------------------------------------------------------------------------------
+\* \+ and - add and subtract systems (parallel connection)
+\* \* multiply systems (series connection)
+\ / right divide -- sys1/sys2 means sys1\*inv(sys2)
+\- \\ left divide -- sys1\\sys2 means inv(sys1)\*sys2
+\ ^ powers of a given system
+\ ' pertransposition
+\ .' transposition of input/output map
+\ .\* element-by-element multiplication
+\ [..] concatenate models along inputs or outputs
+\ lti/stack stack models/arrays along some array dimension
+\ lti/inv inverse of an LTI system
+\ lti/conj complex conjugation of model coefficients
+\
+**Matrix equation solvers and linear algebra**
+--------------------------------------------------------------------------------
+\* lyap, dlyap solve Lyapunov equations
+\ lyapchol, dlyapchol square-root Lyapunov solvers
+\* care, dare solve algebraic Riccati equations
+\ gcare, gdare generalized Riccati solvers
+\ bdschur block diagonalization of a square matrix
+\
+**Additional functions**
+--------------------------------------------------------------------------------
+\* :func:`gangof4` generate the Gang of 4 sensitivity plots
+\* :func:`linspace` generate a set of numbers that are linearly
+ spaced
+\* :func:`logspace` generate a set of numbers that are
+ logarithmically spaced
+\* :func:`unwrap` unwrap a phase angle to give a continuous curve
+== ========================== ================================================
"""
def ss(*args):
@@ -752,9 +786,9 @@
def ngrid():
"""Nichols chart grid.
- Usage
- =====
- ngrid()
+ Examples
+ --------
+ >>> ngrid()
"""
from nichols import nichols_grid
nichols_grid()
@@ -815,134 +849,81 @@
% len(args))
return margins[0], margins[1], margins[3], margins[4]
-#
-# Modifications to scipy.signal functions
-#
-# Redefine lsim to use lsim2
-def lsim(sys, U=None, T=None, X0=None, **keywords):
- """Simulate the output of a linear system
- Examples
- --------
- >>> T, yout, xout = lsim(sys, u, T, X0)
+def dcgain(*args):
+ '''
+ Compute the gain of the system in steady state.
- Parameters
- ----------
- sys: StateSpace, or TransferFunction
- LTI system to simulate
- u: input array giving input at each time T
- T: time steps at which the input is defined
- X0: initial condition (optional, default = 0)
+ The function takes either 1, 2, 3, or 4 parameters:
- Returns
- -------
- T: time values of the output
- yout: response of the system
- xout: time evolution of the state vector
- """
- # Convert the system to an signal.lti for simulation
- #! This should send a warning for MIMO systems
- ltiobjs = sys.returnScipySignalLti()
- ltiobj = ltiobjs[0][0]
-
- return sp.signal.lsim2(ltiobj, U, T, X0, **keywords)
-
-#! Redefine step to use lsim2
-#! Not yet implemented
-def step(*args, **keywords):
- """Step response of a linear system
-
- Examples
- --------
- >>> T, yout = step(sys, T, X0)
-
Parameters
----------
- sys: StateSpace, or TransferFunction
- T: array
- T is the time vector (optional; autocomputed if not given)
- X0: array
- X0 is the initial condition (optional; default = 0)
+ A, B, C, D: array-like
+ A linear system in state space form.
+ Z, P, k: array-like, array-like, number
+ A linear system in zero, pole, gain form.
+ num, den: array-like
+ A linear system in transfer function form.
+ sys: Lti (StateSpace or TransferFunction)
+ A linear system object.
Returns
-------
- T: array
- Time values of the output
- yout: array
- response of the system
- """
- sys = args[0]
- ltiobjs = sys.returnScipySignalLti()
- ltiobj = ltiobjs[0][0]
+ gain: matrix
+ The gain of each output versus each input:
+ :math:`y = gain \cdot u`
+
+ Notes
+ -----
+ This function is only useful for systems with invertible system
+ matrix ``A``.
+
+ All systems are first converted to state space form. The function then
+ computes:
+
+ .. math:: gain = - C \cdot A^{-1} \cdot B + D
+ '''
+ #Convert the parameters to state space form
+ if len(args) == 4:
+ A, B, C, D = args
+ sys = ss(A, B, C, D)
+ elif len(args) == 3:
+ Z, P, k = args
+ A, B, C, D = zpk2ss(Z, P, k)
+ sys = ss(A, B, C, D)
+ elif len(args) == 2:
+ num, den = args
+ sys = tf2ss(num, den)
+ elif len(args) == 1:
+ sys, = args
+ sys = ss(sys)
+ else:
+ raise ValueError("Function ``dcgain`` needs either 1, 2, 3 or 4 "
+ "arguments.")
+ #gain = - C * A**-1 * B + D
+ gain = sys.D - sys.C * sys.A.I * sys.B
+ return gain
- out = sp.signal.step(ltiobj, **keywords)
- yout = []
- yout.append(np.mat(out[0]))
- yout.append(out[1])
- yout = tuple(yout)
- return yout
+# Simulation routines
+# Call corresponding functions in timeresp, with arguments transposed
-# Redefine initial to use lsim2
-#! Not yet implemented (uses step for now)
-def initial(*args, **keywords):
- """Initial condition response of a linear system
+def step(sys, T=None, X0=0., input=0, output=0, **keywords):
+ T, yout = timeresp.StepResponse(sys, T, X0, input, output,
+ transpose = True, **keywords)
+ return T, yout
- Examples
- --------
- >>> T, yout = initial(sys, T, X0)
+def impulse(sys, T=None, X0=0., input=0, output=0, **keywords):
+ T, yout = timeresp.ImpulseResponse(sys, T, X0, input, output,
+ transpose = True, **keywords)
+ return T, yout
- Parameters
- ----------
- sys: StateSpace, or TransferFunction
- T: array
- T is the time vector (optional; autocomputed if not given)
- X0: array
- X0 is the initial condition (optional; default = 0)
+def initial(sys, T=None, X0=0., input=0, output=0, **keywords):
+ T, yout = timeresp.InitialResponse(sys, T, X0, input, output,
+ transpose = True, **keywords)
+ return T, yout
- Returns
- -------
- T: array
- Time values of the output
- yout: array
- response of the system
-
- """
- sys = args[0]
- ltiobjs = sys.returnScipySignalLti()
- ltiobj = ltiobjs[0][0]
-
- yout = sp.signal.initial(ltiobj, **keywords)
- return np.mat(yout)
-
-# Redefine impulse to use initial()
-#! Not yet implemented (uses impulse for now)
-def impulse(*args, **keywords):
- """Impulse response of a linear system
-
- Examples
- --------
- >>> T, yout = impulse(sys, T, X0)
-
- Parameters
- ----------
- sys: StateSpace, or TransferFunction
- T: array
- T is the time vector (optional; autocomputed if not given)
- X0: array
- X0 is the initial condition (optional; default = 0)
-
- Returns
- -------
- T: array
- Time values of the output
- yout: array
- response of the system
-
- """
- sys = args[0]
- ltiobjs = sys.returnScipySignalLti()
- ltiobj = ltiobjs[0][0]
-
- yout = sp.signal.impulse(ltiobj, **keywords)
- return np.mat(yout)
+def lsim(sys, U=0., T=None, X0=0., **keywords):
+ T, yout, xout = timeresp.ForcedResponse(sys, T, U, X0,
+ transpose = True, **keywords)
+ return T, yout, xout
Modified: trunk/src/statesp.py
===================================================================
--- trunk/src/statesp.py 2011-06-18 00:20:49 UTC (rev 160)
+++ trunk/src/statesp.py 2011-06-18 19:13:09 UTC (rev 161)
@@ -1,3 +1,4 @@
+# -*- coding: utf-8 -*-
"""stateSpace.py
State space representation and functions.
@@ -72,8 +73,9 @@
"""
-from numpy import all, angle, any, array, concatenate, cos, delete, dot, \
- empty, exp, eye, matrix, ones, pi, poly, poly1d, roots, shape, sin, zeros
+from numpy import all, angle, any, array, asarray, concatenate, cos, delete, \
+ dot, empty, exp, eye, matrix, ones, pi, poly, poly1d, roots, shape, sin, \
+ zeros
from numpy.random import rand, randn
from numpy.linalg import inv, det, solve
from numpy.linalg.linalg import LinAlgError
@@ -427,8 +429,8 @@
for i in range(self.outputs):
for j in range(self.inputs):
- out[i][j] = lti(self.A, self.B[:, j], self.C[i, :],
- self.D[i, j])
+ out[i][j] = lti(asarray(self.A), asarray(self.B[:, j]),
+ asarray(self.C[i, :]), asarray(self.D[i, j]))
return out
Added: trunk/src/timeresp.py
===================================================================
--- trunk/src/timeresp.py (rev 0)
+++ trunk/src/timeresp.py 2011-06-18 19:13:09 UTC (rev 161)
@@ -0,0 +1,652 @@
+# timesim.py - time-domain simulation routes
+"""timesim.py
+
+Time domain simulation.
+
+Copyright (c) 2010 by SciPy Developers
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+1. Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in the
+ documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the California Institute of Technology nor
+ the names of its contributors may be used to endorse or promote
+ products derived from this software without specific prior
+ written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL CALTECH
+OR THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
+USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
+OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGE.
+
+Author: Eike Welk
+Date: 12 May 2011
+$Id: matlab.py 157 2011-06-17 23:51:46Z murrayrm $
+
+.. _time-series-convention:
+
+Convention for Time Series
+--------------------------
+
+This is a convention for function arguments and return values that
+represent time series: sequences of values that change over time. It is used
+throughout the library, for example in the functions :func:`lsim`, :func:`step`,
+:func:`impulse`, and :func:`initial`.
+
+.. note::
+ This convention is different from the convention used in the library
+ :mod:`scipy.signal`. In Scipy's convention the meaning of rows and columns
+ is interchanged. All 2D values must be transposed when they are used with
+ functions from :mod:`scipy.signal`.
+
+Types:
+
+ * **Arguments** can be **arrays**, **matrices**, or **nested lists**.
+ * **Return values** are **arrays** (not matrices).
+
+The time vector is either 1D, or 2D with shape (1, n)::
+
+ T = [[t1, t2, t3, ..., tn ]]
+
+Input, state, and output all follow the same convention. Columns are different
+points in time, rows are different components. When there is only one row, a
+1D object is accepted or returned, which adds convenience for SISO systems::
+
+ U = [[u1(t1), u1(t2), u1(t3), ..., u1(tn)]
+ [u2(t1), u2(t2), u2(t3), ..., u2(tn)]
+ ...
+ ...
+ [ui(t1), ui(t2), ui(t3), ..., ui(tn)]]
+
+ Same for X, Y
+
+The initial conditions are either 1D, or 2D with shape (j, 1)::
+
+ X0 = [[x1]
+ [x2]
+ ...
+ ...
+ [xj]]
+
+So, U[:,2] is the system's input at the third point in time; and U[1] or U[1,:]
+is the sequence of values for the system's second input.
+
+As all simulation functions return *arrays*, plotting is convenient::
+
+ t, y = step(sys)
+ plot(t, y)
+
+The output of a MIMO system can be plotted like this::
+
+ t, y, x = lsim(sys, u, t)
+ plot(t, y[0], label='y_0')
+ plot(t, y[1], label='y_1')
+
+The convention also works well with the state space form of linear systems. If
+``D`` is the feedthrough *matrix* of a linear system, and ``U`` is its input
+(*matrix* or *array*), then the feedthrough part of the system's response,
+can be computed like this::
+
+ ft = D * U
+
+----------------------------------------------------------------
+
+"""
+
+# Libraries that we make use of
+import scipy as sp # SciPy library (used all over)
+import numpy as np # NumPy library
+from scipy.signal.ltisys import _default_response_times
+from copy import deepcopy
+import warnings
+from statesp import StateSpace, _rss_generate, _convertToStateSpace
+from lti import Lti # base class of StateSpace, TransferFunction
+
+#
+# Solve the systems's differential equations
+#
+
+def _check_convert_array(in_obj, legal_shapes, err_msg_start, squeeze=False,
+ transpose=False):
+ """
+ Helper function for checking array-like parameters.
+
+ * Check type and shape of ``in_obj``.
+ * Convert ``in_obj`` to an array if necessary.
+ * Change shape of ``in_obj`` according to parameter ``squeeze``.
+ * If ``in_obj`` is a scalar (number) it is converted to an array with
+ a legal shape, that is filled with the scalar value.
+
+ The function raises an exception when it detects an error.
+
+ Parameters:
+
+ in_obj: array like object
+ The array or matrix which is checked.
+
+ legal_shapes: list of tuple
+ A list of shapes that in_obj can legally have.
+ The special value "any" means that there can be any
+ number of elements in a certain dimension.
+
+ * ``(2, 3)`` describes an array with 2 rows and 3 columns
+ * ``(2, "any")`` describes an array with 2 rows and any number of
+ columns
+
+ err_msg_start: str
+ String that is prepended to the error messages, when this function
+ raises an exception. It should be used to identify the argument which
+ is currently checked.
+
+ squeeze: bool
+ If True, all dimensions with only one element are removed from the
+ array. If False the array's shape is unmodified.
+
+ For example:
+ ``array([[1,2,3]])`` is converted to ``array([1, 2, 3])``
+
+ transpose: bool
+ If True, assume that input arrays are transposed for the standard
+ format. Used to convert MATLAB-style inputs to our format.
+
+ Returns:
+
+ out_array: array
+ The checked and converted contents of ``in_obj``.
+ """
+ #convert nearly everything to an array.
+ out_array = np.asarray(in_obj)
+ if (transpose):
+ out_array = np.transpose(out_array)
+
+ #Test element data type, elements must be numbers
+ legal_kinds = set(("i", "f", "c")) #integer, float, complex
+ if out_array.dtype.kind not in legal_kinds:
+ err_msg = "Wrong element data type: '{d}'. Array elements " \
+ "must be numbers.".format(d=str(out_array.dtype))
+ raise TypeError(err_msg_start + err_msg)
+
+ #If array is zero dimensional (in_obj is scalar):
+ #create array with legal shape filled with the original value.
+ if out_array.ndim == 0:
+ for s_legal in legal_shapes:
+ #search for shape that does not contain the special symbol any.
+ if "any" in s_legal:
+ continue
+ the_val = out_array[()]
+ out_array = np.empty(s_legal, 'd')
+ out_array.fill(the_val)
+ break
+
+ #Test shape
+ def shape_matches(s_legal, s_actual):
+ """Test if two shape tuples match"""
+ #Array must have required number of dimensions
+ if len(s_legal) != len(s_actual):
+ return False
+ #All dimensions must contain required number of elements. Joker: "all"
+ for n_legal, n_actual in zip(s_legal, s_actual):
+ if n_legal == "any":
+ continue
+ if n_legal != n_actual:
+ return False
+ return True
+
+ #Iterate over legal shapes, and see if any matches out_array's shape.
+ for s_legal in legal_shapes:
+ if shape_matches(s_legal, out_array.shape):
+ break
+ else:
+ legal_shape_str = " or ".join([str(s) for s in legal_shapes])
+ err_msg = "Wrong shape (rows, columns): {a}. Expected: {e}." \
+ .format(e=legal_shape_str, a=str(out_array.shape))
+ raise ValueError(err_msg_start + err_msg)
+
+ #Convert shape
+ if squeeze:
+ out_array = np.squeeze(out_array)
+ #We don't want zero dimensional arrays
+ if out_array.shape == tuple():
+ out_array = out_array.reshape((1,))
+
+ return out_array
+
+# Forced response of a linear system
+def ForcedResponse(sys, T=None, U=0., X0=0., transpose=False, **keywords):
+ """Simulate the output of a linear system.
+
+ As a convenience for parameters `U`, `X0`:
+ Numbers (scalars) are converted to constant arrays with the correct shape.
+ The correct shape is inferred from arguments `sys` and `T`.
+
+ For information on the **shape** of parameters `U`, `T`, `X0` and
+ return values `T`, `yout`, `xout` see: :ref:`time-series-convention`
+
+ Parameters
+ ----------
+ sys: Lti (StateSpace, or TransferFunction)
+ LTI system to simulate
+
+ T: array-like
+ Time steps at which the input is defined, numbers must be (strictly
+ monotonic) increasing.
+
+ U: array-like or number, optional
+ Input array giving input at each time `T` (default = 0).
+
+ If `U` is ``None`` or ``0``, a special algorithm is used. This special
+ algorithm is faster than the general algorithm, which is used otherwise.
+
+ X0: array-like or number, optional
+ Initial condition (default = 0).
+
+ transpose: bool
+ If True, transpose all input and output arrays (for backward
+ compatibility with MATLAB and scipy.signal.lsim)
+
+ **keywords:
+ Additional keyword arguments control the solution algorithm for the
+ differential equations. These arguments are passed on to the function
+ :func:`scipy.integrate.odeint`. See the documentation for
+ :func:`scipy.integrate.odeint` for information about these
+ arguments.
+
+ Returns
+ -------
+ T: array
+ Time values of the output.
+ yout: array
+ Response of the system.
+ xout: array
+ Time evolution of the state vector.
+
+ See Also
+ --------
+ StepResponse, InitialResponse, ImpulseResponse
+
+ Examples
+ --------
+ >>> T, yout, xout = ForcedResponse(sys, T, u, X0)
+ """
+ if not isinstance(sys, Lti):
+ raise TypeError('Parameter ``sys``: must be a ``Lti`` object. '
+ '(For example ``StateSpace`` or ``TransferFunction``)')
+ sys = _convertToStateSpace(sys)
+ A, B, C, D = np.asarray(sys.A), np.asarray(sys.B), np.asarray(sys.C), \
+ np.asarray(sys.D)
+# d_type = A.dtype
+ n_states = A.shape[0]
+ n_inputs = B.shape[1]
+
+ #Test if T has shape (n,) or (1, n);
+ #T must be array-like and values must be increasing.
+ #The length of T determines the length of the input vector.
+ if T is None:
+ raise ValueError('Parameter ``T``: must be array-like, and contain '
+ '(strictly monotonic) increasing numbers.')
+ T = _check_convert_array(T, [('any',), (1,'any')],
+ 'Parameter ``T``: ', squeeze=True,
+ transpose = transpose)
+ if not all(T[1:] - T[:-1] > 0):
+ raise ValueError('Parameter ``T``: time values must be '
+ '(strictly monotonic) increasing numbers.')
+ n_steps = len(T) #number of simulation steps
+
+ #create X0 if not given, test if X0 has correct shape
+ X0 = _check_convert_array(X0, [(n_states,), (n_states,1)],
+ 'Parameter ``X0``: ', squeeze=True)
+
+ #Solve the differential equation, copied from scipy.signal.ltisys.
+ dot, squeeze, = np.dot, np.squeeze #Faster and shorter code
+ #Faster algorithm if U is zero
+ if U is None or (isinstance(U, (int, float)) and U == 0):
+ #Function that computes the time derivative of the linear system
+ def f_dot(x, _t):
+ return dot(A,x)
+
+ xout = sp.integrate.odeint(f_dot, X0, T, **keywords)
+ yout = dot(C, xout.T)
+ #General algorithm that interpolates U in between output points
+ else:
+ #Test if U has correct shape and type
+ legal_shapes = [(n_steps,), (1,n_steps)] if n_inputs == 1 else \
+ [(n_inputs, n_steps)]
+ U = _check_convert_array(U, legal_shapes,
+ 'Parameter ``U``: ', squeeze=False,
+ transpose=transpose)
+ #convert 1D array to D2 array with only one row
+ if len(U.shape) == 1:
+ U = U.reshape(1,-1) #pylint: disable=E1103
+
+ #Create a callable that uses linear interpolation to
+ #calculate the input at any time.
+ compute_u = sp.interpolate.interp1d(T, U, kind='linear', copy=False,
+ axis=-1, bounds_error=False,
+ fill_value=0)
+
+ #Function that computes the time derivative of the linear system
+ def f_dot(x, t):
+ return dot(A,x) + squeeze(dot(B,compute_u([t])))
+
+ xout = sp.integrate.odeint(f_dot, X0, T, **keywords)
+ yout = dot(C, xout.T) + dot(D, U)
+
+ yout = squeeze(yout)
+ xout = xout.T
+
+ # See if we need to transpose the data back into MATLAB form
+ if (transpose):
+ T = np.transpose(T)
+ yout = np.transpose(yout)
+ xout = np.transpose(xout)
+
+ return T, yout, xout
+
+def StepResponse(sys, T=None, X0=0., input=0, output=0, \
+ 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.
+
+ For information on the **shape** of parameters `T`, `X0` and
+ return values `T`, `yout` see: :ref:`time-series-convention`
+
+ Parameters
+ ----------
+ sys: StateSpace, or TransferFunction
+ LTI system to simulate
+
+ T: array-like object, optional
+ Time vector (argument is autocomputed if not given)
+
+ X0: array-like or number, optional
+ Initial condition (default = 0)
+
+ Numbers are converted to constant arrays with the correct shape.
+
+ input: int
+ Index of the input that will be used in this simulation.
+
+ output: int
+ Index of the output that will be used in this simulation.
+
+ transpose: bool
+ If True, transpose all input and output arrays (for backward
+ compatibility with MATLAB and scipy.signal.lsim)
+
+ **keywords:
+ Additional keyword arguments control the solution algorithm for the
+ differential equations. These arguments are passed on to the function
+ :func:`lsim`, which in turn passes them on to
+ :func:`scipy.integrate.odeint`. See the documentation for
+ :func:`scipy.integrate.odeint` for information about these
+ arguments.
+
+ Returns
+ -------
+ T: array
+ Time values of the output
+
+ yout: array
+ Response of the system
+
+ See Also
+ --------
+ lsim, initial, impulse
+
+ Examples
+ --------
+ >>> T, yout = step(sys, T, X0)
+ """
+ sys = _convertToStateSpace(sys)
+ sys = _mimo2siso(sys, input, output, warn_conversion=True)
+ if T is None:
+ T = _default_response_times(sys.A, 100)
+ U = np.ones_like(T)
+
+ T, yout, _xout = ForcedResponse(sys, T, U, X0,
+ transpose = transpose, **keywords)
+
+ return T, yout
+
+
+def InitialResponse(sys, T=None, X0=0., input=0, output=0, 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.
+
+ For information on the **shape** of parameters `T`, `X0` and
+ return values `T`, `yout` see: :ref:`time-series-convention`
+
+ Parameters
+ ----------
+ sys: StateSpace, or TransferFunction
+ LTI system to simulate
+
+ T: array-like object, optional
+ Time vector (argument is autocomputed if not given)
+
+ X0: array-like object or number, optional
+ Initial condition (default = 0)
+
+ Numbers are converted to constant arrays with the correct shape.
+
+ input: int
+ Index of the input that will be used in this simulation.
+
+ output: int
+ Index of the output that will be used in this simulation.
+
+ transpose: bool
+ If True, transpose all input and output arrays (for backward
+ compatibility with MATLAB and scipy.signal.lsim)
+
+ **keywords:
+ Additional keyword arguments control the solution algorithm for the
+ differential equations. These arguments are passed on to the function
+ :func:`lsim`, which in turn passes them on to
+ :func:`scipy.integrate.odeint`. See the documentation for
+ :func:`scipy.integrate.odeint` for information about these
+ arguments.
+
+
+ Returns
+ -------
+ T: array
+ Time values of the output
+ yout: array
+ Response of the system
+
+ See Also
+ --------
+ lsim, step, impulse
+
+ Examples
+ --------
+ >>> T, yout = InitialResponsesys, T, X0)
+ """
+ sys = _convertToStateSpace(sys)
+ sys = _mimo2siso(sys, input, output, warn_conversion=True)
+ #Create time and input vectors; checking is done in ForcedResponse(...)
+ #The initial vector X0 is created in ForcedResponse(...) if necessary
+ if T is None:
+ T = _default_response_times(sys.A, 100)
+ U = np.zeros_like(T)
+
+ T, yout, _xout = ForcedResponse(sys, T, U, X0, **keywords)
+ return T, yout
+
+
+def ImpulseResponse(sys, T=None, X0=0., input=0, output=0,
+ transpose=False, **keywords):
+ #pylint: disable=W0622
+ """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.
+
+ For information on the **shape** of parameters `T`, `X0` and
+ return values `T`, `yout` see: :ref:`time-series-convention`
+
+ Parameters
+ ----------
+ sys: StateSpace, TransferFunction
+ LTI system to simulate
+
+ T: array-like object, optional
+ Time vector (argument is autocomputed if not given)
+
+ X0: array-like object or number, optional
+ Initial condition (default = 0)
+
+ Numbers are converted to constant arrays with the correct shape.
+
+ input: int
+ Index of the input that will be used in this simulation.
+
+ output: int
+ Index of the output that will be used in this simulation.
+
+ transpose: bool
+ If True, transpose all input and output arrays (for backward
+ compatibility with MATLAB and scipy.signal.lsim)
+
+ **keywords:
+ Additional keyword arguments control the solution algorithm for the
+ differential equations. These arguments are passed on to the function
+ :func:`lsim`, which in turn passes them on to
+ :func:`scipy.integrate.odeint`. See the documentation for
+ :func:`scipy.integrate.odeint` for information about these
+ arguments.
+
+
+ Returns
+ -------
+ T: array
+ Time values of the output
+ yout: array
+ Response of the system
+
+ See Also
+ --------
+ lsim, step, initial
+
+ Examples
+ --------
+ >>> T, yout = ImpulseResponse(sys, T, X0)
+ """
+ sys = _convertToStateSpace(sys)
+ sys = _mimo2siso(sys, input, output, warn_conversion=True)
+
+ #System has direct feedthrough, can't simulate impulse response numerically.
+ if np.any(sys.D != 0):
+ warnings.warn('System has direct feedthrough: ``D != 0``. The infinite '
+ 'impulse at ``t=0`` does not appear in the output. \n'
+ 'Results may be meaningless!')
+
+ #create X0 if not given, test if X0 has correct shape.
+ #Must be done here because it is used for computations here.
+ n_states = sys.A.shape[0]
+ X0 = _check_convert_array(X0, [(n_states,), (n_states,1)],
+ 'Parameter ``X0``: \n', squeeze=True)
+
+ #Compute new X0 that contains the impulse
+ #We can't put the impulse into U because there is no numerical
+ #representation for it (infinitesimally short, infinitely high).
+ #See also: http://www.mathworks.com/support/tech-notes/1900/1901.html
+ B = np.asarray(sys.B).squeeze()
+ new_X0 = B + X0
+
+ #Compute T and U, no checks necessary, they will be checked in lsim
+ if T is None:
+ T = _default_response_times(sys.A, 100)
+ U = np.zeros_like(T)
+
+ T, yout, _xout = ForcedResponse(sys, T, U, new_X0, \
+ transpose=transpose, **keywords)
+ return T, yout
+
+#! TODO: this function probably belongs in a different file
+def _mimo2siso(sys, input, output, warn_conversion):
+ #pylint: disable=W0622
+ """
+ Convert a MIMO system to a SISO system. (Convert a system with multiple
+ inputs and/or outputs, to a system with a single input and output.)
+
+ The input and output that are used in the SISO system can be selected
+ with the parameters ``input`` and ``output``. All other inputs are set
+ to 0, all other outputs are ignored.
+
+ If ``sys`` is already a SISO system, it will be returned unaltered.
+
+ Parameters:
+
+ sys...
[truncated message content] |
