Skip to content

Allow np.array or np.matrix for state space matrices, operations #314

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
murrayrm merged 4 commits into from
Jun 30, 2019
Merged

Allow np.array or np.matrix for state space matrices, operations #314

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
0c1e983
allow np.array or np.matrix for state space matrices + unit tests
murrayrm Jun 15, 2019
4c9c2c8
fixed docstring example error + added TODO's in response to review
murrayrm Jun 17, 2019
4f48c5d
allow D=0 to broadcast to correct size for MIMO statespace systems
murrayrm Jun 23, 2019
7c9d9a6
change np.bmat -> np.block in timeresp.py
murrayrm Jun 30, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
91 changes: 46 additions & 45 deletions control/bdalg.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,25 +53,25 @@

"""

import scipy as sp
import numpy as np
from . import xferfcn as tf
from . import statesp as ss
from . import frdata as frd

__all__ = ['series', 'parallel', 'negate', 'feedback', 'append', 'connect']


def series(sys1, *sysn):
"""Return the series connection (... \* sys3 \*) sys2 \* sys1

Parameters
----------
sys1: scalar, StateSpace, TransferFunction, or FRD
sysn: other scalars, StateSpaces, TransferFunctions, or FRDs
sys1 : scalar, StateSpace, TransferFunction, or FRD
sysn : other scalars, StateSpaces, TransferFunctions, or FRDs

Returns
-------
out: scalar, StateSpace, or TransferFunction
out : scalar, StateSpace, or TransferFunction

Raises
------
Expand Down Expand Up @@ -105,18 +105,19 @@ def series(sys1, *sysn):
from functools import reduce
return reduce(lambda x, y:y*x, sysn, sys1)


def parallel(sys1, *sysn):
"""
Return the parallel connection sys1 + sys2 (+ sys3 + ...)

Parameters
----------
sys1: scalar, StateSpace, TransferFunction, or FRD
*sysn: other scalars, StateSpaces, TransferFunctions, or FRDs
sys1 : scalar, StateSpace, TransferFunction, or FRD
*sysn : other scalars, StateSpaces, TransferFunctions, or FRDs

Returns
-------
out: scalar, StateSpace, or TransferFunction
out : scalar, StateSpace, or TransferFunction

Raises
------
Expand Down Expand Up @@ -150,17 +151,18 @@ def parallel(sys1, *sysn):
from functools import reduce
return reduce(lambda x, y:x+y, sysn, sys1)


def negate(sys):
"""
Return the negative of a system.

Parameters
----------
sys: StateSpace, TransferFunction or FRD
sys : StateSpace, TransferFunction or FRD

Returns
-------
out: StateSpace or TransferFunction
out : StateSpace or TransferFunction

Notes
-----
Expand All @@ -177,7 +179,6 @@ def negate(sys):
>>> sys2 = negate(sys1) # Same as sys2 = -sys1.

"""

return -sys;

#! TODO: expand to allow sys2 default to work in MIMO case?
Expand All @@ -187,18 +188,18 @@ def feedback(sys1, sys2=1, sign=-1):

Parameters
----------
sys1: scalar, StateSpace, TransferFunction, FRD
The primary plant.
sys2: scalar, StateSpace, TransferFunction, FRD
The feedback plant (often a feedback controller).
sys1 : scalar, StateSpace, TransferFunction, FRD
The primary process.
sys2 : scalar, StateSpace, TransferFunction, FRD
The feedback process (often a feedback controller).
sign: scalar
The sign of feedback. `sign` = -1 indicates negative feedback, and
`sign` = 1 indicates positive feedback. `sign` is an optional
argument; it assumes a value of -1 if not specified.

Returns
-------
out: StateSpace or TransferFunction
out : StateSpace or TransferFunction

Raises
------
Expand Down Expand Up @@ -256,7 +257,7 @@ def feedback(sys1, sys2=1, sign=-1):
return sys1.feedback(sys2, sign)

def append(*sys):
'''append(sys1, sys2, ..., sysn)
"""append(sys1, sys2, ..., sysn)

Group models by appending their inputs and outputs

Expand All @@ -279,42 +280,40 @@ def append(*sys):

Examples
--------
>>> sys1 = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.")
>>> sys2 = ss("-1.", "1.", "1.", "0.")
>>> 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
"""Index-based interconnection of an LTI 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.
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!!!!
NOTE: Inputs and outputs are indexed starting at 1 and negative values
correspond to a negative feedback interconnection.

Parameters
----------
sys: StateSpace Transferfunction
sys : StateSpace Transferfunction
System to be connected
Q: 2d array
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
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
no connection is made. Inputs and outputs are indexed starting at 1.
inputv : 1D array
list of final external inputs
outputv: 1d array
outputv : 1D array
list of final external outputs

Returns
Expand All @@ -324,28 +323,30 @@ def connect(sys, Q, inputv, outputv):

Examples
--------
>>> sys1 = ss("1. -2; 3. -4", "5.; 7", "6, 8", "9.")
>>> sys2 = ss("-1.", "1.", "1.", "0.")
>>> 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
>>> Q = [[1, 2], [2, -1]]) # negative feedback interconnection

@roryyorke roryyorke Jun 16, 2019

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

trailing parenthesis here looks wrong. Do we run these with doctest?

@murrayrm murrayrm Jun 17, 2019

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed.

We are not currently using doctest, but I ran the revised version manually to make sure it was OK.

>>> sysc = connect(sys, Q, [2], [1, 2])
'''

"""
# first connect
K = sp.zeros( (sys.inputs, sys.outputs) )
for r in sp.array(Q).astype(int):
K = np.zeros((sys.inputs, sys.outputs))
for r in np.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)
sys = sys.feedback(np.array(K), sign=1)

# now trim
Ytrim = sp.zeros( (len(outputv), sys.outputs) )
Utrim = sp.zeros( (sys.inputs, len(inputv)) )
Ytrim = np.zeros((len(outputv), sys.outputs))
Utrim = np.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)

return Ytrim * sys * Utrim
13 changes: 13 additions & 0 deletions control/config.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,13 +7,17 @@
# files. For now, you can just choose between MATLAB and FBS default
# values.

import warnings

# Bode plot defaults
bode_dB = False # Bode plot magnitude units
bode_deg = True # Bode Plot phase units
bode_Hz = False # Bode plot frequency units
bode_number_of_samples = None # Bode plot number of samples
bode_feature_periphery_decade = 1.0 # Bode plot feature periphery in decades

# State space module variables
_use_numpy_matrix = True # Decide whether to use numpy.marix

def reset_defaults():
"""Reset configuration values to their default values."""
Expand All @@ -22,6 +26,7 @@ def reset_defaults():
global bode_Hz; bode_Hz = False
global bode_number_of_samples; bode_number_of_samples = None
global bode_feature_periphery_decade; bode_feature_periphery_decade = 1.0
global _use_numpy_matrix; _use_numpy_matrix = True


# Set defaults to match MATLAB
Expand All @@ -36,6 +41,7 @@ def use_matlab_defaults():
global bode_dB; bode_dB = True
global bode_deg; bode_deg = True
global bode_Hz; bode_Hz = True
global _use_numpy_matrix; _use_numpy_matrix = True


# Set defaults to match FBS (Astrom and Murray)
Expand All @@ -52,3 +58,10 @@ def use_fbs_defaults():
global bode_deg; bode_deg = True
global bode_Hz; bode_Hz = False


# Decide whether to use numpy.matrix for state space operations
def use_numpy_matrix(flag=True, warn=True):
if flag and warn:
warnings.warn("Return type numpy.matrix is soon to be deprecated.",
stacklevel=2)
global _use_numpy_matrix; _use_numpy_matrix = flag
17 changes: 11 additions & 6 deletions control/frdata.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@
from warnings import warn
import numpy as np
from numpy import angle, array, empty, ones, \
real, imag, matrix, absolute, eye, linalg, where, dot
real, imag, absolute, eye, linalg, where, dot
from scipy.interpolate import splprep, splev
from .lti import LTI

Expand Down Expand Up @@ -81,6 +81,10 @@ class FRD(LTI):

"""

# Allow NDarray * StateSpace to give StateSpace._rmul_() priority
# https://docs.scipy.org/doc/numpy/reference/arrays.classes.html
__array_priority__ = 11 # override ndarray and matrix types

epsw = 1e-8

def __init__(self, *args, **kwargs):
Expand Down Expand Up @@ -436,12 +440,13 @@ def feedback(self, other=1, sign=-1):
# TODO: vectorize this
# TODO: handle omega re-mapping
for k, w in enumerate(other.omega):
fresp[:, :, k] = self.fresp[:, :, k].view(type=matrix)* \
fresp[:, :, k] = np.dot(
self.fresp[:, :, k],
linalg.solve(
eye(self.inputs) +
other.fresp[:, :, k].view(type=matrix) *
self.fresp[:, :, k].view(type=matrix),
eye(self.inputs))
eye(self.inputs)
+ np.dot(other.fresp[:, :, k], self.fresp[:, :, k]),
eye(self.inputs))

@roryyorke roryyorke Jun 16, 2019 •
edited
Loading

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I know this is not the point of this change, but is linalg(X, eye(shape)) any better than np.linalg.inv? I'd have thought a "right divide" (linalg.solve with appropriate transposes) would have been a more standard approach.

@murrayrm murrayrm Jun 17, 2019

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree that this could use sp.linalg.solve instead, but since I wasn't the original author I just did the conversion to remove the np.matrix dependence. I've added a TODO for someone to fix this, so at least we don't forget about it (in principle).

)

return FRD(fresp, other.omega, smooth=(self.ifunc is not None))

Expand Down
Loading