Skip to content

Commit af08a30

Browse files
murrayrmmurrayrm
authored
Merge branch 'master' into fix-latex
2 parents c268c7c + 4785145 commit af08a30

File tree

13 files changed

+466
-182
lines changed

13 files changed

+466
-182
lines changed

‎.travis.yml

Lines changed: 45 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -24,23 +24,56 @@ python:
2424
#
2525
# We also want to test with and without slycot
2626
env:
27-
- SCIPY=scipy SLYCOT=slycot # default, with slycot
27+
- SCIPY=scipy SLYCOT=conda # default, with slycot via conda
2828
- SCIPY=scipy SLYCOT= # default, w/out slycot
2929
- SCIPY="scipy==0.19.1" SLYCOT= # legacy support, w/out slycot
3030

31-
# Exclude combinations that are very unlikely (and don't work)
31+
# Add optional builds that test against latest version of slycot
32+
jobs:
33+
include:
34+
- name: "linux, Python 2.7, slycot=source"
35+
os: linux
36+
dist: xenial
37+
services: xvfb
38+
python: "2.7"
39+
env: SCIPY=scipy SLYCOT=source
40+
- name: "linux, Python 3.7, slycot=source"
41+
os: linux
42+
dist: xenial
43+
services: xvfb
44+
python: "3.7"
45+
env: SCIPY=scipy SLYCOT=source
46+
3247
matrix:
48+
# Exclude combinations that are very unlikely (and don't work)
3349
exclude:
3450
- python: "3.7" # python3.7 should use latest scipy
3551
env: SCIPY="scipy==0.19.1" SLYCOT=
3652

53+
allow_failures:
54+
- name: "linux, Python 2.7, slycot=source"
55+
os: linux
56+
dist: xenial
57+
services: xvfb
58+
python: "2.7"
59+
env: SCIPY=scipy SLYCOT=source
60+
- name: "linux, Python 3.7, slycot=source"
61+
os: linux
62+
dist: xenial
63+
services: xvfb
64+
python: "3.7"
65+
env: SCIPY=scipy SLYCOT=source
66+
3767
# install required system libraries
3868
before_install:
3969
# Install gfortran for testing slycot; use apt-get instead of conda in
4070
# order to include the proper CXXABI dependency (updated in GCC 4.9)
41-
- if [[ "$SLYCOT" != "" ]]; then
71+
# Note: these commands should match the slycot .travis.yml configuration
72+
- if [[ "$SLYCOT" = "source" ]]; then
4273
sudo apt-get update -qq;
74+
sudo apt-get install liblapack-dev libblas-dev;
4375
sudo apt-get install gfortran;
76+
sudo apt-get install cmake;
4477
fi
4578
# use miniconda to install numpy/scipy, to avoid lengthy build from source
4679
- if [[ "$TRAVIS_PYTHON_VERSION" == "2.7" ]]; then
@@ -57,9 +90,8 @@ before_install:
5790
- conda info -a
5891
- conda create -q -n test-environment python="$TRAVIS_PYTHON_VERSION" pip coverage
5992
- source activate test-environment
60-
# Install openblas if slycot is being used
61-
# also install scikit-build for the build process
62-
- if [[ "$SLYCOT" != "" ]]; then
93+
# Install scikit-build for the build process if slycot is being used
94+
- if [[ "$SLYCOT" = "source" ]]; then
6395
conda install openblas;
6496
conda install -c conda-forge scikit-build;
6597
fi
@@ -72,13 +104,15 @@ before_install:
72104
install:
73105
# Install packages needed by python-control
74106
- conda install $SCIPY matplotlib
75-
# Build slycot from source
76-
# For python 3, need to provide pointer to python library
77-
# Use "Unix Makefiles" as generator, because Ninja cannot handle Fortran
78-
#! git clone https://github.com/repagh/Slycot.git slycot;
79-
- if [[ "$SLYCOT" != "" ]]; then
107+
108+
# Figure out how to build slycot
109+
# source: use "Unix Makefiles" as generator; Ninja cannot handle Fortran
110+
# conda: use pre-compiled version of slycot on conda-forge
111+
- if [[ "$SLYCOT" = "source" ]]; then
80112
git clone https://github.com/python-control/Slycot.git slycot;
81113
cd slycot; python setup.py install -G "Unix Makefiles"; cd ..;
114+
elif [[ "$SLYCOT" = "conda" ]]; then
115+
conda install -c conda-forge slycot;
82116
fi
83117

84118
# command to run tests

‎control/bdalg.py

Lines changed: 0 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -169,11 +169,6 @@ def negate(sys):
169169
This function is a wrapper for the __neg__ function in the StateSpace and
170170
TransferFunction classes. The output type is the same as the input type.
171171
172-
If both systems have a defined timebase (dt = 0 for continuous time,
173-
dt > 0 for discrete time), then the timebase for both systems must
174-
match. If only one of the system has a timebase, the return
175-
timebase will be set to match it.
176-
177172
Examples
178173
--------
179174
>>> sys2 = negate(sys1) # Same as sys2 = -sys1.

‎control/iosys.py

Lines changed: 7 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -212,12 +212,12 @@ def __repr__(self):
212212

213213
def __str__(self):
214214
"""String representation of an input/output system"""
215-
str = "System: " + (self.name if self.name else "(none)") + "\n"
216-
str += "Inputs (%d): " % self.ninputs
215+
str = "System: " + (self.name if self.name else "(None)") + "\n"
216+
str += "Inputs (%s): " % self.ninputs
217217
for key in self.input_index: str += key + ", "
218-
str += "\nOutputs (%d): " % self.noutputs
218+
str += "\nOutputs (%s): " % self.noutputs
219219
for key in self.output_index: str += key + ", "
220-
str += "\nStates (%d): " % self.nstates
220+
str += "\nStates (%s): " % self.nstates
221221
for key in self.state_index: str += key + ", "
222222
return str
223223

@@ -317,13 +317,8 @@ def __add__(sys1, sys2):
317317
ninputs = sys1.ninputs
318318
noutputs = sys1.noutputs
319319

320-
# Make sure timebase are compatible
321-
dt = _find_timebase(sys1, sys2)
322-
if dt is False:
323-
raise ValueError("System timebases are not compabile")
324-
325320
# Create a new system to handle the composition
326-
newsys = InterconnectedSystem((sys1, sys2), dt=dt)
321+
newsys = InterconnectedSystem((sys1, sys2))
327322

328323
# Set up the input map
329324
newsys.set_input_map(np.concatenate(
@@ -937,6 +932,7 @@ def __init__(self, syslist, connections=[], inplist=[], outlist=[],
937932
system_count = 0
938933
for sys in syslist:
939934
# Make sure time bases are consistent
935+
# TODO: Use lti._find_timebase() instead?
940936
if dt is None and sys.dt is not None:
941937
# Timebase was not specified; set to match this system
942938
dt = sys.dt
@@ -948,7 +944,7 @@ def __init__(self, syslist, connections=[], inplist=[], outlist=[],
948944
sys.nstates is None:
949945
raise TypeError("System '%s' must define number of inputs, "
950946
"outputs, states in order to be connected" %
951-
sys)
947+
sys.name)
952948

953949
# Keep track of the offsets into the states, inputs, outputs
954950
self.input_offset.append(ninputs)

‎control/mateqn.py

Lines changed: 23 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -41,7 +41,8 @@
4141
Author: Bjorn Olofsson
4242
"""
4343

44-
from numpy import shape, size, array, asarray, copy, zeros, eye, dot
44+
from numpy import shape, size, asarray, copy, zeros, eye, dot, \
45+
finfo, inexact, atleast_2d
4546
from scipy.linalg import eigvals, solve_discrete_are, solve
4647
from .exception import ControlSlycot, ControlArgument
4748
from .statesp import _ssmatrix
@@ -122,7 +123,7 @@ def lyap(A, Q, C=None, E=None):
122123
if size(Q) > 1 and shape(Q)[0] != shape(Q)[1]:
123124
raise ControlArgument("Q must be a quadratic matrix.")
124125

125-
if not (asarray(Q) == asarray(Q).T).all():
126+
if not _is_symmetric(Q):
126127
raise ControlArgument("Q must be a symmetric matrix.")
127128

128129
# Solve the Lyapunov equation by calling Slycot function sb03md
@@ -188,7 +189,7 @@ def lyap(A, Q, C=None, E=None):
188189
raise ControlArgument("E must be a square matrix with the same \
189190
dimension as A.")
190191

191-
if not (asarray(Q) == asarray(Q).T).all():
192+
if not _is_symmetric(Q):
192193
raise ControlArgument("Q must be a symmetric matrix.")
193194

194195
# Make sure we have access to the write slicot routine
@@ -309,7 +310,7 @@ def dlyap(A,Q,C=None,E=None):
309310
if size(Q) > 1 and shape(Q)[0] != shape(Q)[1]:
310311
raise ControlArgument("Q must be a quadratic matrix.")
311312

312-
if not (asarray(Q) == asarray(Q).T).all():
313+
if not _is_symmetric(Q):
313314
raise ControlArgument("Q must be a symmetric matrix.")
314315

315316
# Solve the Lyapunov equation by calling the Slycot function sb03md
@@ -371,7 +372,7 @@ def dlyap(A,Q,C=None,E=None):
371372
raise ControlArgument("E must be a square matrix with the same \
372373
dimension as A.")
373374

374-
if not (asarray(Q) == asarray(Q).T).all():
375+
if not _is_symmetric(Q):
375376
raise ControlArgument("Q must be a symmetric matrix.")
376377

377378
# Solve the generalized Lyapunov equation by calling Slycot
@@ -500,10 +501,10 @@ def care(A, B, Q, R=None, S=None, E=None, stabilizing=True):
500501
size(B) == 1 and n > 1:
501502
raise ControlArgument("Incompatible dimensions of B matrix.")
502503

503-
if not (asarray(Q) == asarray(Q).T).all():
504+
if not _is_symmetric(Q):
504505
raise ControlArgument("Q must be a symmetric matrix.")
505506

506-
if not (asarray(R) == asarray(R).T).all():
507+
if not _is_symmetric(R):
507508
raise ControlArgument("R must be a symmetric matrix.")
508509

509510
# Create back-up of arrays needed for later computations
@@ -603,10 +604,10 @@ def care(A, B, Q, R=None, S=None, E=None, stabilizing=True):
603604
size(S) == 1 and m > 1:
604605
raise ControlArgument("Incompatible dimensions of S matrix.")
605606

606-
if not (asarray(Q) == asarray(Q).T).all():
607+
if not _is_symmetric(Q):
607608
raise ControlArgument("Q must be a symmetric matrix.")
608609

609-
if not (asarray(R) == asarray(R).T).all():
610+
if not _is_symmetric(R):
610611
raise ControlArgument("R must be a symmetric matrix.")
611612

612613
# Create back-up of arrays needed for later computations
@@ -775,10 +776,10 @@ def dare_old(A, B, Q, R, S=None, E=None, stabilizing=True):
775776
size(B) == 1 and n > 1:
776777
raise ControlArgument("Incompatible dimensions of B matrix.")
777778

778-
if not (asarray(Q) == asarray(Q).T).all():
779+
if not _is_symmetric(Q):
779780
raise ControlArgument("Q must be a symmetric matrix.")
780781

781-
if not (asarray(R) == asarray(R).T).all():
782+
if not _is_symmetric(R):
782783
raise ControlArgument("R must be a symmetric matrix.")
783784

784785
# Create back-up of arrays needed for later computations
@@ -882,10 +883,10 @@ def dare_old(A, B, Q, R, S=None, E=None, stabilizing=True):
882883
size(S) == 1 and m > 1:
883884
raise ControlArgument("Incompatible dimensions of S matrix.")
884885

885-
if not (asarray(Q) == asarray(Q).T).all():
886+
if not _is_symmetric(Q):
886887
raise ControlArgument("Q must be a symmetric matrix.")
887888

888-
if not (asarray(R) == asarray(R).T).all():
889+
if not _is_symmetric(R):
889890
raise ControlArgument("R must be a symmetric matrix.")
890891

891892
# Create back-up of arrays needed for later computations
@@ -960,3 +961,12 @@ def dare_old(A, B, Q, R, S=None, E=None, stabilizing=True):
960961
# Invalid set of input parameters
961962
else:
962963
raise ControlArgument("Invalid set of input parameters.")
964+
965+
966+
def _is_symmetric(M):
967+
M = atleast_2d(M)
968+
if isinstance(M[0, 0], inexact):
969+
eps = finfo(M.dtype).eps
970+
return ((M - M.T) < eps).all()
971+
else:
972+
return (M == M.T).all()

‎control/pzmap.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -81,7 +81,7 @@ def pzmap(sys, Plot=True, grid=False, title='Pole Zero Map'):
8181
The system's zeros.
8282
"""
8383
# Get parameter values
84-
Plot = config._get_param('rlocus', 'Plot', grid, True)
84+
Plot = config._get_param('rlocus', 'Plot', Plot, True)
8585
grid = config._get_param('rlocus', 'grid', grid, False)
8686

8787
if not isinstance(sys, LTI):

‎control/statefbk.py

Lines changed: 76 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -43,10 +43,11 @@
4343
import numpy as np
4444
import scipy as sp
4545
from . import statesp
46+
from .mateqn import care
4647
from .statesp import _ssmatrix
4748
from .exception import ControlSlycot, ControlArgument, ControlDimension
4849

49-
__all__ = ['ctrb', 'obsv', 'gram', 'place', 'place_varga', 'lqr', 'acker']
50+
__all__ = ['ctrb', 'obsv', 'gram', 'place', 'place_varga', 'lqr', 'lqe', 'acker']
5051

5152

5253
# Pole placement
@@ -219,6 +220,76 @@ def place_varga(A, B, p, dtime=False, alpha=None):
219220
# Return the gain matrix, with MATLAB gain convention
220221
return _ssmatrix(-F)
221222

223+
# contributed by Sawyer B. Fuller <minster@uw.edu>
224+
def lqe(A, G, C, QN, RN, NN=None):
225+
"""lqe(A, G, C, QN, RN, [, N])
226+
227+
Linear quadratic estimator design (Kalman filter) for continuous-time
228+
systems. Given the system
229+
230+
Given the system
231+
.. math::
232+
x = Ax + Bu + Gw
233+
y = Cx + Du + v
234+
235+
with unbiased process noise w and measurement noise v with covariances
236+
237+
.. math:: E{ww'} = QN, E{vv'} = RN, E{wv'} = NN
238+
239+
The lqe() function computes the observer gain matrix L such that the
240+
stationary (non-time-varying) Kalman filter
241+
242+
.. math:: x_e = A x_e + B u + L(y - C x_e - D u)
243+
244+
produces a state estimate that x_e that minimizes the expected squared error
245+
using the sensor measurements y. The noise cross-correlation `NN` is set to
246+
zero when omitted.
247+
248+
Parameters
249+
----------
250+
A, G: 2-d array
251+
Dynamics and noise input matrices
252+
QN, RN: 2-d array
253+
Process and sensor noise covariance matrices
254+
NN: 2-d array, optional
255+
Cross covariance matrix
256+
257+
Returns
258+
-------
259+
L: 2D array
260+
Kalman estimator gain
261+
P: 2D array
262+
Solution to Riccati equation
263+
.. math::
264+
A P + P A^T - (P C^T + G N) R^-1 (C P + N^T G^T) + G Q G^T = 0
265+
E: 1D array
266+
Eigenvalues of estimator poles eig(A - L C)
267+
268+
269+
Examples
270+
--------
271+
>>> K, P, E = lqe(A, G, C, QN, RN)
272+
>>> K, P, E = lqe(A, G, C, QN, RN, NN)
273+
274+
See Also
275+
--------
276+
lqr
277+
"""
278+
279+
# TODO: incorporate cross-covariance NN, something like this,
280+
# which doesn't work for some reason
281+
#if NN is None:
282+
# NN = np.zeros(QN.size(0),RN.size(1))
283+
#NG = G @ NN
284+
285+
#LT, P, E = lqr(A.T, C.T, G @ QN @ G.T, RN)
286+
#P, E, LT = care(A.T, C.T, G @ QN @ G.T, RN)
287+
A, G, C = np.array(A, ndmin=2), np.array(G, ndmin=2), np.array(C, ndmin=2)
288+
QN, RN = np.array(QN, ndmin=2), np.array(RN, ndmin=2)
289+
P, E, LT = care(A.T, C.T, np.dot(np.dot(G, QN), G.T), RN)
290+
return _ssmatrix(LT.T), _ssmatrix(P), _ssmatrix(E)
291+
292+
222293
# Contributed by Roberto Bucher <roberto.bucher@supsi.ch>
223294
def acker(A, B, poles):
224295
"""Pole placement using Ackermann method
@@ -307,6 +378,10 @@ def lqr(*args, **keywords):
307378
>>> K, S, E = lqr(sys, Q, R, [N])
308379
>>> K, S, E = lqr(A, B, Q, R, [N])
309380
381+
See Also
382+
--------
383+
lqe
384+
310385
"""
311386

312387
# Make sure that SLICOT is installed

0 commit comments

Comments
 (0)