Skip to main content
This is a pre-production deployment of Warehouse. Changes made here affect the production instance of PyPI (
Help us improve Python packaging - Donate today!

A control systems library for Python3

Project Description

|Travis-CI| |License| |Join the chat at Gitter| |Coveralls|


A (yet another) systems and controls module for Python3. See ````
file for details (spoiler alert it's MIT Licensed).


This is a work in progress. It's just almost about to have the first release,
version 0.1. So, please adjust the expectations accordingly.


harold library is written with the ambition of becoming a daily work-horse of
a control engineer/student/researcher with complete access to the source
code and also with full rights in terms of open source.

The common licensing problems where there are more users than the available
number of licenses or, in a commercial setting, exporting commercial code
with the product, or building a product using commercial software and so on
are naturally avoided.

In short, the purpose of this library is to disclose the development of
control algorithms. harold fully supports the mantra of reproducible research
and thus aims to provide the means of accessible and transparent computational
development tools.

harold also follows the Pythonic way of naming things and hence deviates from
the matlab way of working. This might add some learning effort on the user
but hopefully the function names are self-explanatory enough to ease the pain.

harold currently comes with two representation models: State and Transfer.
All is based around such object classes and the rest is just (lots of)

To demonstrate some examples, you create a state space or transfer function
representation via the following syntax::

G = State([[0,1],[-2,-1]],[[0],[1]],[[1,0]]) # you can skip D if zero
H = Transfer([1,2,1],[1,4,5,6])
J = State(5) # Gain models are supported
K = Transfer(3, dt=0.1) # in discrete time too.

If we type `H` we get::

Continous-Time Transfer function with:
1 input(s) and 1 output(s)

Poles(real) Poles(imag) Zeros(real) Zeros(imag)
------------- ------------- ------------- -------------
-0.5 1.32288 -1 0
-0.5 -1.32288 -1 0
-3 0

End of Transfer object description

This might be better represented in the future versions. Or you can implement
your own infostring because you might be computing the damping of the
zeros/poles all the time. Important to note is the matrix syntax of NumPy
is a little awkward and we don't have much to change that now. However,
in turn, `harold` tries to be as forgiving as possible expecting the least
formal representation of arrays.

For example, for MIMO Transfer representations you can simply give a common
numerator or denominator and it will be completed internally to a matching
size, e.g.,::

G = Transfer(1,[ [ [1,1],[1,2,1] ],[ [2,3],[1,3,3,1] ] ])

Notice how we had to increase the nesting of the bracket level to make it
contain arrays for a 2x2 system. This is, say in matlab, is done via curly
braces but for Python, brackets are natural for a list syntax.

While we are at it, here is a Bode Plot of a 348-state Clamped Beam model
from the SLICOT collections::

import as sio # needed to read .mat files
M = sio.loadmat('beam.mat')
A = M['A'].todense()
B = M['B']
C = M['C']
G = State(A,B,C)

.. image::

By the way, you can't convert this model to a Transfer model (including in
matlab etc.) because of the numerical error build-up. Thus, harold tries to
stay with whatever is given to it.

What can I do with `harold`?

Pretty much everything that is related to control systems
number-crunching as long as I can keep up with the implementation. The
testing phase is coming to an end, however, certain functionality is
(hopefully) in place:

- State, Transfer representations ✓
- Pole/zero computations ✓
- Minimal realizations ✓
- Controllabe/Observable Hessenberg Forms ✓
- Pole/Zero cancellation distance computations ✓
- Discretizations/Undiscretizations ✓
- Frequency Response calculations ✓
- LTI Algebra (add,mul,sub,feedback) ✓ (except feedback and series)
- Polynomial Operations ✓
- Riccati and Lyapunov equation solvers ✓ (Riccati solver is moved to SciPy)
- H2- and HInfinity norms of systems ✓
- Bode, Nyquist static plots ✓

The coded-but-not-tested parts are the essential synthesis techniques

- Pole placement (uses the scipy versions)
- LQR/LQG (the functionality is there but needs to be wrapped in a function)

Still waiting to be coded

- Step, Impulse responses
- Hinf, H2 controllers
- Lead, Lag and notch filters with semi-recommendations about where to place them

and so on.

For example, you would like to have a row compressed, observable Hessenberg
form such that the rows of C matrix is compressed on the right end. You might
wonder why you would need such a representation. Welcome to the weird world
of control theory!::

M = np.array([[-6.5, 0.5, 6.5, -6.5, 0. , 1. , 0. ],
[-0.5, -5.5, -5.5, 5.5, 2. , 1. , 2. ],
[-0.5, 0.5, 0.5, -6.5, 3. , 4. , 3. ],
[-0.5, 0.5, -5.5, -0.5, 3. , 2. , 3. ],
[ 1. , 1. , 0. , 0. , 0. , 0. , 0. ]])
# Slice the matrix into 4 pieces with south west piece has (1,4) shape
A , B , C , D = matrix_slice(M, (1, 4), corner='sw')
a,b,c,k = staircase(A, B, C, form='o', invert=True)

array([[ 0. , -6. , -7.48528137, 0. ],
[-6. , 0. , 9.48528137, 0. ],
[ 0. , 0. , -6. , 0. ],
[ 0. , 0. , 1.41421356, -6. ]])

array([[ 3. , 2. , 3. ],
[ 3. , 4. , 3. ],
[ 1.41421356, 0. , 1.41421356],
[ 1.41421356, 1.41421356, 1.41421356]])

array([[ 0. , 0. , 0. , 1.41421356]])

array([ 1., 1.])

Here `k` is the size of the block that is identified as observable at each
step of the staircase. We can deduce that two of the modes are already
unobservable since the upper left 2x2 block does not interact with the lower
right two modes (A21 block is identically zero). Let's check the minimality

a,b,c = minimalrealization(A,B,C)

array([[-6. , 0. ],
[ 1.41421356, -6. ]])

array([[ 1.41421356, 0. , 1.41421356],
[ 1.41421356, 1.41421356, 1.41421356]])

array([[ 0. , 1.41421356]])

which gives a 2x2 system as we have expected.

Suppose you have some polynomials and would like to compute the LCM/GCD. Then
you can go about it via::

a , b = haroldlcm([1,3,0,-4], [1,-4,-3,18], [1,-4,3], [1,-2,-8])

which returns::

(array([ 1., -7., 3., 59., -68., -132., 144.])

[array([ 1., -10., 33., -36.]),
array([ 1., -3., -6., 8.]),
array([ 1., -3., -12., 20., 48.]),
array([ 1., -5., 1., 21., -18.])]

Here `a` is the least common multiple and `b` is the array of polynomials
that are needed to be multiplied by the original polynomials (in the order
of appearance) to obtain the LCM.

Another point-of-interest is the interactive plots that are promising. That
would hopefully minimize the right-click mania that follows almost every
plotting command in every commercial software for Bode, Nyquist, Sensitivity,
Coherence and others.

What about ...?

Yes, yes, LMIs are coming. I have to learn cvxpy a bit faster. Other stuff
you need to let me know what the need is.


See the `Sphinx documentation`_ .

Useful Links

1- There is already an almost-matured control toolbox which is led by
Richard Murray et al. (`click for the Github page`_ ) and it can perform
already most of the essential tasks. Hence, if you want to have
something that resembles the basics of matlab control toolbox give it a
try. However, it is mostly limited to SISO tools and also relies on
SLICOT library which can lead to licensing problems for nontrivial tasks.

2- By the way, if you are interested in robust control you would
probably appreciate the `Skogestad-Python`_ project.

Help Wanted!

harold is built on rainy days and boring evenings. If you are desperately
missing out a feature, don't be shy and contact. User-feedback has higher
priority over the general development. Or shout out in the Gitter chatroom.

Or if you want to jump into development, PR submissions are more than welcome.


If you have questions/comments feel free to shoot one to

.. _click for the Github page:
.. _Sphinx documentation:
.. _Skogestad-Python:

.. |License| image::
.. |Join the chat at Gitter| image::
.. |Travis-CI| image::
.. |Coveralls| image::


+ Sanitized the circular dependencies a bit more
+ minimal_realization is changed to accept models instead of A,B,C triplet
+ minimal_realization for Transfer uses the pole zero cancellation check
+ more housekeeping and bug fixes
+ added damping, natural frequency properties of poles
- state_to_transfer does not return minimal realizations (per request)

+ More tests
+ bode, nyquist plots with matplotlib
+ Rewritten the transmission_zeros to improve accuracy
- Removed the single file and replaced it with modular files.
- Refactored Riccati solvers to SciPy official repo
- Lyapunov solver safety net is moot. Created PR #6775 in SciPy

+ Added Riccati solvers
+ More documentation
+ Added safety net for lyapunov solvers in case there is no solution

+ Added Lyapunov solvers
+ Fixed many bugs
- Removed block diag and switched to scipy version

+ Initial versioning and packaging.
+ Adding documentation and Sphinx integration.
+ Basically everything there is.

Release History

This version
History Node


History Node


History Node


History Node


History Node


History Node


History Node


Download Files

Download the file for your platform. If you're not sure which to choose, learn more about installing packages.

Filename, Size & Hash SHA256 Hash Help File Type Python Version Upload Date
(80.4 kB) Copy SHA256 Hash SHA256
Wheel py3 May 1, 2017
(959.2 kB) Copy SHA256 Hash SHA256
Source None May 1, 2017

Supported By

Elastic Elastic Search Pingdom Pingdom Monitoring Dyn Dyn DNS Sentry Sentry Error Logging CloudAMQP CloudAMQP RabbitMQ Heroku Heroku PaaS Kabu Creative Kabu Creative UX & Design Fastly Fastly CDN DigiCert DigiCert EV Certificate Google Google Cloud Servers DreamHost DreamHost Log Hosting