Johan de Bruin has attempted to solve this in 2019. His code did not pass the test. It was found that it is possible to rewrite the tf function to work with Polynomial, however numerous other functions in the repo are dependant on the coefficients of numerator and denominator of the fransferfunction. These numerators and denominators are both of type Polynomial and therefore their coefficients are given in the reverse order of that what is expected from the user, namely in increasing order of powers instead of decreasing order, as poly1d gives it. The transfer function can be made to display in decreasing order and takes inputs in decreasing order, however is no obvious way to ensure that the call tf.numerator.coef() would return the values in decreasing order. A possible solution to this would be to replace all instances of tf.numerator.coef with tf.numerator.coef[::-1], however this will cause confusion to anyone who wants code further on it.

His code can be found on the "Polynomial" branch of https://github.com/DeBruinJohan/Skogestad-Python/tree/Polynomial

At the moment our tf object can only handle rational functions with the option one time delay. This causes an issue when calculating transfer functions like 1 + G if G contains a time delay. This issue is for allowing some kind of resolution for this so that we can calculate closed loop transfer functions for systems with delays.

While sympy is a very useful module for control calculations, it is better to do most of our calculations numerically rather than symbolically.

utils.mimotf.mimotf_slice(G, rows, cols) is redundant - it does the same as G[rows, cols].

Hi there! I would like to start by saying that I admire the great work you're doing! I'm new to Python and open source projects in general, and I've been trying to move from Matlab to Python in my control/signal processing projects. I've been looking for a way to work with internal time delays in control systems, and found out about this amazing project!

I'm having a problem when trying to calculate the zeros from a closed-loop transfer function. The zeros that I'm obtaining have significant numerical inaccuracy when compared to the PID-controller's zeros (they should be the same in this case). The following sample code generates the problem that I'm talking about. Also, there is another (minor) problem: when asking the zeros of a MIMO system, it returns the system's poles and zeros (at the end of the code, there is a part where this error occurred for me):

import numpy as np
import robustcontrol as rc


######################### second order system plus dead time (SOPDT) ##########################
K = 0.9
w_n = 2
csi = 0.3
theta = 0.1

# System's transfer function
G = K*rc.utils.tf([w_n**2], [1, 2*csi*w_n, w_n**2], deadtime = theta)


##################################### Ziegler-Nichols PID Tuning ################################

# critical gain Kcr 
K_cr = rc.utils.margins(G)[0][0] 

# critical oscillation period 
P_cr = np.abs(2*np.pi/rc.utils.margins(G)[3][0])

# PID gains 
Kp = 0.6*K_cr
Ki = Kp/(0.5*P_cr)*rc.utils.tf([1], [1, 0])
Kd = Kp*0.125*P_cr* rc.utils.tf([1, 0],[1])

# PID controller transfer function
C_ZN = Kp + Ki + Kd

################################# Calculating system zeros ###########################################

# Delay free model (zero-order padé approx)
G_delay_free = K*rc.utils.tf([w_n**2], [1, 2*csi*w_n, w_n**2])


# closed-loop considering a delay-free model
H = G_delay_free*C_ZN/(1 + G_delay_free*C_ZN)

# Zeros should be equal in both cases, but there is a significant numerical error (the same does not occur in matlab or using python-control package)
H_zeros = H.zeros()
# returns  array([-2.49594595+0.0365365j, -2.49594595-0.0365365j])  

C_ZN_zeros = C_ZN.zeros()
# returns array([-2.49652882+1.43703181e-05j, -2.49652882-1.43703181e-05j])



# When asking the zeros of a mimo transfer function, it return poles and zeros
H_mimo_style = rc.utils.mimotf(H)
H_mimo_style.zeros()

Thanks for your attention and sorry for the bad coding or bad English, I've been trying to exercise both of those things recently. I'm still learning how to use GitHub too! :)

As explained in Skogestad 4.5.2

runall.py was written before pytest was popular. I think we could rewrite it to use pytest instead. I've created a pytest runner for ipynb files in the Dynamics and Control repository.

The code needs to be reworked to be able to upload as a package. This will involve considerable reworking of the directory structures. The process to get a package up is described here. Note that anyone doing this should not upload on my behalf, but rather work on the structure so that it is possible for me (@alchemyst) to upload.

The tf class is perhaps reaching the end of its usefulness. We don't really need to maintain a whole separate class from the ones found in control and harold or sympy.

Implement code which will obtain a minimal realisation for state space models.

Clone code refers to code which has been copied and pasted with perhaps only minor modifications. The current code base has many clone code snippets. You can find the clones by running Clone Digger. I would like to have clone digger return an empty report.

It should be simple to calculate MIMO feedback using utils.feedback

The method currently used to calculate polynomial greatest common divisor and lowest common multiple involves calculating the roots of the polynomials then using logic to determine the gcd and lcm. This method only works well for lower order polynomials (up to approximately 5th order), due to the inherent difficulty in accurately calculating polynomial roots for higher order polynomials. This causes erroneous results in many of the functions in utils when applied to large mimo systems (bigger than 2x2), including utils.poles, utils.zeros and utils.tf2ss, as many of the interim calculations involve large polynomials. A new method should be created that does not rely on calculating roots. Also, Euclid's algorithm (or any method closely related to this) should not be used as it is not applicable to polynomials in floating point arithmetic.

@ilayn has given the following recommendation:

For LCM, I'm using a similar result to Karcanias, Mitrouli, "System theoretic based characterisation and computation of the least common multiple of a set of polynomials", Lin Alg App, 381, 2004. The results is acceptable but never stress tested it to the extremes.

For GCD, Sylvester-matrix based methods are sufficiently good for medium difficulty/size problems.

The following code can be used to test if the new method is working correctly:

import utils
import numpy
s = utils.tf([1,0], 1)
G = utils.mimotf([[1/((s+1)*(s+2)), 1/((s+2)*(s+3)), 1/((s+3)*(s+4)), 1/((s+4)*(s+1))],
                 [1/((s+1)*(s+2)), 1/((s+2)*(s+3)), 1/((s+3)*(s+4)), 1/((s+4)*(s+1))],
                 [1/((s+1)*(s+2)), 1/((s+2)*(s+3)), 1/((s+3)*(s+4)), 1/((s+4)*(s+1))],
                 [1/((s+1)*(s+2)), 1/((s+2)*(s+3)), 1/((s+3)*(s+4)), 1/((s+4)*(s+1))]])
poles = utils.poles(G)
assert poles.shape == (4,)
#There should be exactly 4 poles
Gpoles = [G(i) for i in poles]
AbsoluteMaximums = numpy.array([abs(i).max() for i in Gpoles])
SmallestMaximum = Maximums.min()
assert SmallestMaximum > 10**6
#Evaluating G(poles) should give at least one very large (infinite) element for each pole

We aim for high-quality readable code. We don't want Fs in Code Climate. Click here and keep fixing things until there's nothing in that list.

At least utils.py and utilsplot.py should be moved to a subdirectory in preparation for packaging.

Note that this involves more work than simply converting the .py files to .ipynb. Currently every one of the .py files are part of the testing process, so to complete this task you also have to rework the testing machinery.

The goal is also specifically to replace the .py files with the .ipynb files and not to duplicate them as this increases the maintenance load.

Currently tests fail for newer matplotilb, this needs to be fixed and the limit removed from the environment.yml file.

As explained in Skogestad 4.4.3

There is a procedure highlighted in this paper

Currently opening a file and pressing the play button doesn't work as expected. We need a doc for this.

The PyPI compliance milestone needs to be fleshed out such that we can do each task relatively easily.

Once #317 is closed, we should change utilsplot.py to run the doctests (see utils.py for how this works) and add to .travis.yml so that it can run automatically.

hypothesis is a testing framework which generates cases automatically.

With the new mimotf object, an error occurs when a plant functions is called with a tf element only containing constants. Consider G22 in the the code below.

from utils import tf, mimotf

s = tf([1, 0])

G11 = (s - 2.5) / (s - 2)
G12 = -(0.1 * s + 1) / (s + 2)
G21 = (s - 2.5) / (0.1 * s + 1)
G22 = 1

G = mimotf([[G11, G12], [G21, G22]])
print G(0)

We need to make the current code compatable with Python3 and Python 2. This procedure illustrates how.

The readthedocs page builds correctly but is empty. We need at least a proper landing page so that users can get their bearings.

The project has been registered on Read the Docs. The documentation build is succeeding on a local machine, but failing due to the dependencies on Numpy.

The current testing code uses hard-coded dependencies. I think we should move this to something like conda env

Hi Carl, I followed this link on how to submit wiki changes https://gist.github.com/larrybotha/10650410 please advise if there is a different way of doing this that you prefer.

Link to new wiki: https://github.com/Bern-Alberts/Skogestad-Python.wiki.git

The doctests in utilsplot.py are not passing

from utils import InternalDelay, tf, mimotf

NoModel = tf([0], [1])

G11 = tf([0.0039 ], [85, 1], deadtime=85)
G12 = NoModel
G13 = NoModel
G14 = tf([279.067728070848, 0.862100005149841], [7122.56780575453, 186.547394918102, 1],deadtime= 45)
G15 = tf([0.009], [30, 1], deadtime=180)
G16 = NoModel
G21 = tf([0.0032], [160, 1], deadtime=85)
G22 = tf([-0.1],[ 85, 1], deadtime=55)
G23 = tf([-0.055], [160, 1], deadtime=100)
G24 = NoModel
G25 = NoModel
G26 = NoModel
G31 = tf([0.0033],[120, 1], deadtime=80)
G32 = NoModel
G33 = tf([-0.065], [120, 1], deadtime=95)
G34 = NoModel
G35 = NoModel
G36 = NoModel

G = mimotf([
[G11, G12, G13, G14, G15, G16],
[G21, G22, G23, G24, G25, G26],
[G31, G32, G33, G34, G35, G36]
])
Gd = InternalDelay(G)

Running the above code will result in the following error.

Skogestad-Python\robustcontrol\InternalDelay.py:185: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray
num, den, delays = [numpy.array(i) for i in [num, den, delays]]
Traceback (most recent call last):
File "c:\Users\user1\Documents\James\Skogestad-Python\robustcontrol\trasfer_func.py", line 33, in
Gd = InternalDelay(G)
File "C:\Users\user1\Documents\James\Skogestad-Python\robustcontrol\InternalDelay.py", line 81, in init
matrices = InternalDelay.from_mimotf(sys).get_matrices()
File "C:\Users\user1\Documents\James\Skogestad-Python\robustcontrol\InternalDelay.py", line 136, in from_mimotf
return InternalDelay.from_tf_coefficients(num, den, delays)
File "C:\Users\user1\Documents\James\Skogestad-Python\robustcontrol\InternalDelay.py", line 214, in from_tf_coefficients
Gss_i = harold.transfer_to_state(harold.Transfer(num_i, den_i))
File "C:\Miniconda3\envs\control\lib\site-packages\harold_classes.py", line 2729, in transfer_to_state
D[rows, cols] = num[rows][cols]/den[rows][cols]
IndexError: index 3 is out of bounds for axis 1 with size 3

Observation:
If I remove any one of the delay from first output (G11, G14 or G15) this error will not occur.

With the previous version of utils it was not possible to generate the state-space representation for the scaled transfer function in the attachment. The updated version permits this, however the shape of C is incorrect. The following code reproduces the problem:

import numpy as np
import utils

s = utils.tf([1, 0], 1)

G = utils.mimotf([
    [0.66/(6.7*s + 1), -0.61/(8.64*s + 1), -0.0049/(9.06*s + 1)],
    [1.11/(3.25*s + 1), -2.36/(5.0*s + 1), -0.012/(7.09*s + 1)],
    [-34.68/(8.15*s + 1), 46.2/(10.9*s + 1), 0.87*(11.61*s + 1)/((3.89*s + 1)*(18.8*s + 1))]
    ])

A, B, C, D = utils.tf2ss(G)

Nstates = A.shape[0]
Ninputs = B.shape[1]
Noutputs = C.shape[0]

assert A.shape[0] == A.shape[1], "A must be square"
assert B.shape[0] == Nstates, "B must have the same number of rows as A"
assert C.shape[1] == Nstates, "C must have the same number or columns as A"
assert D.shape[0] == Ninputs, "D must have the same number of rows as C"

There is currently code in tf.simplify which attempts to round transfer function coefficients and keep them as integers. This causes significant pain to people who have real values like 1.231e-2 in their tranfer function description.

I have added a flag to the tf object which allows one to bypass this part of the simplification. However, this will still cause issues when doing "normal" operations with these objects, as they are simplfied on creation.

If we don't do this, we need to figure out how else to normalise the transfer functions so that we can get repeatable output. Is tf([1], [3, 3]) better than for instance tf([0.3333], [1, 1])?

Coveralls hasn't updated since 2015. Something is wrong - it may be something to do with the build script. Coveralls is configured to run after success in .travis.yml.

This function seems to work for some inputs but not for others. The issue I have picked up specifically is that recalculating the minimal realisation can lead to a different eigenvalues of the A matrix, and therefore different poles. Test code is as follows:

import utils
import numpy

s = utils.tf([1, 0], 1)
U2Y1nd = -0.0499*(9.41*s+1)/(48.2*s**2+23.1*s+1)
U1Y2nd = 0.203/(10*s+1)
U2Y2nd = 0.0463*(27.5*s+1)/(27.3*s**2+10.4*s+1)
U3Y2nd = 0.121/(7.71*s+1)
U4Y2nd = 0.156/(2.25*s+1)
U3Y3nd = -1/(0.04*s**2+0.4*s+1)
U1Y4nd = 0.549*(-4.65*s+1)/(1.06*s**2+2.05*s+1)
U2Y4nd = -0.296*(18*s+1)/(0.58*s**2+2.76*s+1)
U4Y4nd = -1.24*(-16.1*s+1)/(21.8*s**2+9.33*s+1)

Gnd = utils.mimotf([[0, U2Y1nd, 0, 0],[U1Y2nd, U2Y2nd, U3Y2nd, U4Y2nd], [0, 0, U3Y3nd, 0], [U1Y4nd, U2Y4nd, 0, U4Y4nd]])
A, B, C, _ = utils.tf2ss(Gnd)

Amin1, _, _ = utils.minimal_realisation(A, B, C)
Amin2, _, _ = utils.minimal_realisation(A, B, C)

Eigs1, _ = numpy.linalg.eig(Amin1)
Eigs2, _ = numpy.linalg.eig(Amin2)

assert numpy.array_equal(sorted(Eigs1, key=lambda x: x.real), sorted(Eigs2, key=lambda x: x.real)), "Eigenvalues should be equal" 

I believe the problem is in the remove_uncontrollable_or_unobservable_states function, specifically in the following lines:

    elif uncontrollable is False and unobservable is True:
        P[0:rank, :] = con_or_obs_matrix[0:rank, :]

        # matrix will replace the dependent columns in P to make P invertible
        replace_matrix = numpy.matrix(numpy.random.random((m, n_states)))

        # make P invertible
        P[rank:n_states, :] = replace_matrix

        P_inv = numpy.linalg.inv(P)

    A_new = P*a*P_inv
    A_new = numpy.delete(A_new, numpy.s_[rank:n_states], 1)
    A_new = numpy.delete(A_new, numpy.s_[rank:n_states], 0)

    B_new = P*b
    B_new = numpy.delete(B_new, numpy.s_[rank:n_states], 0)

    C_new = c*P_inv
    C_new = numpy.delete(C_new, numpy.s_[rank:n_states], 1)

given that we entered the elif part of the if statement. I suspect that the issue is in the deletes, are these supposed to be different depending on which path we went through in the if statement? Unfortunately I don't know what the algorithm should look like.

We currently have a tf.step method to handle step responses, but we should probably also add tf.lsim which would work similarly to scipy.signal.lsim.

The polygcd and polylcm methods currently rely on the Euclidean algorithm, however this only works for numbers that can be represented exactly in binary, so most float values will cause the algorithm to break. Forcing the coefficient of the polynomials to integers won't work, because any rounding will also cause the Euclidean algorithm to break. polygcd is used to get the zeros of mimotf objects and polylcm is used in tf2ss, so both of these methods will probably have issues if using floating point coefficients. A new approach for the gcd and lcm that doesn't rely on the Euclidean algorithm is needed.

It is more robust to use isinstance(variable, typetocheck) than type(variable) == typetocheck.

Not working at the moment due to requiring a virtualenv install.