harold is a control systems toolbox that provides essential control engineering tools. In this notebook, we cover the basics of working with models and standard operations using them.
Let us first import the common libraries that will be used throughout this document. Star importing is generally frowned upon due to possible name clashes however since we are exclusively working with harold this simplifies the coding.
In [1]:
import numpy as np
from harold import *
import matplotlib.pyplot as plt
harold defines two classes to represent the dynamic models, namely, Transfer and State classes. The usage is straightforward. Since NumPy arrays are a bit laborious to enter, the inputs are designed to be as flexible as possible. We are going to define a few models for sake of demonstration
In [2]:
G1 = Transfer([1, -1],[1, -2, 1, 0], dt=0.1) # discrete
G2 = Transfer([[1, [1, 3]],[0, [1, 2]]], [[[1, 2], [1, 0, -4]],[1, [1, -3]]])
G3 = State([[0, 1], [-0.1, -0.5]], [[0], [1]], [0, 3.5], 1, dt=0.1) # discrete
G4 = State([[-4, -2, 2], [0, -5, 2], [4, -3, -4]], [[1, 1], [0, 1], [2, -1]], [-1, 5, 2])
Internally, the list objects are converted to NumPy arrays. Notice that when we are defining $G_4$, we skipped the $D$ matrix since it is zero. Discrete-time models can be defined with the same syntax but also adding dt keyword with a value in seconds, for example, State(a, b, c, d, dt=1). Getting used to using keywords in Python, both improves the readability and also avoids problems such as State(a, b, c, 1) which would take 1 as the $D$ element instead of the sampling period.
The default __repr__ method is to provide basic information about the objects however can be converted to any other custom string.
In [3]:
G1
Out[3]:
Now, we can directly identify the pole zero structure right away. Though MIMO state representations should be preferred exclusively, MIMO transfer models explicitly require a list of lists structure. Additionally, common denominators can also be given and it will be expanded to every entry in the denominator list of lists, say, $$G_5 = \frac{1}{s^2+5s+1}\begin{bmatrix}s+1&2\\s+3&4\end{bmatrix}$$
In [4]:
G5 = Transfer([[[1, 1], 2], [[1, 3], 4]], [1, 5, 1])
print('Numerator :\n', G5.num, '\nDenominator :\n', G5.den)
One typical inconvenience entering State models is that, it is much easier to type the whole array $\begin{bmatrix}A&B\\C&D\end{bmatrix}$ and slice the individual arrays. However, this needs the extra slice indices which might be also laborious. The matrix_slice() function is provided for this purpose. Hence, suppose that the model data for $G_4$ above is given as an array already. Then creating $G_4$ is basically providing the slices and using the *-argument expansion of Python.
In [5]:
M = np.array([[-4, -2, 2, 1, 1], [0, -5, 2, 0, 1], [4, -3, -4, 2, -1], [-1, 5, 2, 0, 0]])
G4 = State(*matrix_slice(M, corner_shape=[3, 3], corner='nw'))
G4.matrices
Out[5]:
We could have sliced the array based on the shape of the C array and arrive at the same result
In [6]:
M = np.array([[-4, -2, 2, 1, 1],
[0, -5, 2, 0, 1],
[4, -3, -4, 2, -1],
[-1, 5, 2, 0, 0]])
G4 = State(*matrix_slice(M, corner_shape=[1, 3], corner='sw'))
Further properties can be accessed via the typical dot syntax, let's check the model data of $G_1$
In [7]:
G1.polynomials
Out[7]:
Also the model data of a State representation, for example $G_3$, or just the $A$ matrix of $G_4$
In [8]:
G3.matrices
Out[8]:
In [9]:
G4.a
Out[9]:
Finally we can also create random State models via the random_stable_model
In [10]:
G = random_state_model(5, p=3, m=2)
G
Out[10]:
We can also force the random model to have more chance to have oscillatory modes by changing the probability distribution of the selected poles. Assume that we want a discrete-time model with mostly poles on the imaginary axis and occasional integrators.
In [11]:
G = random_state_model(20, p=3, m=2, dt=0.01, prob_dist=[0, 0, 0.1, 0.9],stable=False) # 90% osc. modes, 10% integrators
G
Out[11]:
In order to convert one model representation to another, the relevant functions are transfer_to_state() and state_to_transfer(). Notice that, typically the conversion from a Transfer model to a State model leads to a nonminimal representations.
In [12]:
H = transfer_to_state(G2)
J = state_to_transfer(G3)
In [13]:
H
Out[13]:
In [14]:
J
Out[14]:
The minimality is an essential property of any representations for reliable computations and hence we can use minimal_realization() function. This function uses a distance metric to the closest rank-deficient matrix pencil for State models and a straightforward walk-over the poles and zeros for Transfer representations. As usual, the tolerance for decision can be adjusted.
Note that, MIMO poles and zeros don't necessarily cancel each other even though the values are identical. This is because for MIMO representations pole and zero directionalities can be different. Let's see this for the representation $H$ we have obtained above.
In [15]:
Hm = minimal_realization(H)
Hm
Out[15]:
We can see that while on of the $-2$ poles is removed, the others are not cancelled.
Discretization of models is more of an art than science. Though the methods are quite sound, selecting the right method and the relevant sampling period is up to designer's choice.
In [16]:
G6 = Transfer([1], [1, 1.4, 1])
G6d = discretize(G6, dt=0.3, method='zoh')
G6d
Out[16]:
Current known discretization methods are:
In [17]:
from harold._global_constants import _KnownDiscretizationMethods as km
print(*km, sep='\n')
The forward xxx, backward xxx are aliases of >> and <<. Similarly, bilinear is an alias of tustin.
Moreover, converting the models back and forth after each manipulation necessitates some sort of a memory of what has happened before. For this purpose, harold encodes the method of discretization into the model such that when converted back the discretization information is searched first and then if not found, the defaults are used. Let's start by using $G_3$ above which we have not provided any method but for $G_{6d}$ the default method is already implemented.
In [18]:
G3.DiscretizedWith is None
Out[18]:
In [19]:
G6d.DiscretizedWith
Out[19]:
Now, if we convert back these models, $G_3$ will be converted using the default method tustin however $G_{6d}$ will be converted via zero-order hold method. Had this information not present, the resulting continuous model would be slightly different than what we started with.
In [20]:
G6dcz = undiscretize(G6d)
G6dcz
Out[20]:
In [21]:
G6dct = undiscretize(G6d, method='tustin')
G6dct
Out[21]:
Basic algebraic operations and feedback operations are implemented via typical *,+,-,@ and feedback() function. The shape and sampling time compatibility is checked before the operations and errors are raised in case of mismatches.
One exception is the operation with a SISO model over a MIMO model. For example, $G_1$ is a SISO model and $G_3$ is a sampling time matching MIMO model, hence following is allowed which multiplies each entry of $G_3$ with $G_1$. For multiplication of MIMO models matrix multiplication rules and hence matrix multiplication operator @ is followed.
In [22]:
G1 * G3
Out[22]:
In [23]:
G4 @ G2
Out[23]:
In [24]:
CL = feedback(G3, G1)
CL_min = minimal_realization(CL)
CL
Out[24]:
We observe that there is a pole/zero cancellation at 1. which is removed afterwards.
In [25]:
CL_min
Out[25]:
The typical plots that are needed are already available. For frequency domain plotting, the default units are Hz and powers of ten for amplitudes. Motivation for these choices can be found in the book Feedback Systems
In [26]:
impulse_response_plot(G4);
The discrete-time plant plots are automatically drawn as stairs.
In [27]:
G4_d = discretize(G4, 0.1, method='zoh')
impulse_response_plot(G4_d);
In [28]:
bode_plot(G4);
The plot units can be changed via the keywords available
In [29]:
bode_plot(G4, use_db=True, use_hz=False, use_degree=False);
In [30]:
nyquist_plot(G2);
Just to demonstrate the convenience, let's pickup the first section of the LQR example Inverted Pendulum: State-Space Methods for Controller Design from the classic Control Tutorials for MATLAB
In [31]:
# Define some parameters
M, m, b, I, g, l = 0.5, 0.2, 0.1, 0.006, 9.8, 0.3
p = I*(M+m) + M*m*l**2 #denominator for the A and B matrices
A = np.array([[0, 1, 0, 0], [0, -(I+m*l**2)*b/p, (m**2*g*l**2)/p, 0],[0, 0, 0, 1], [0, -(m*l*b)/p, m*g*l*(M+m)/p, 0]])
B = np.array([[0], [(I+m*l**2)/p], [0], [m*l/p]])
C = np.array([[1, 0, 0, 0], [0, 0, 1, 0]])
sys_ss = State(A,B,C)
print('The system is Kalman controllable:', is_kalman_controllable(sys_ss))
Q = C.T @ C
K, X, eigs = lqr(sys_ss, Q) # R = 1 if omitted
print('Controller K gains : ', K)
sys_cl = State(A-B@K, B, C)
t = np.arange(0, 5, 0.01)
r =0.2*np.ones(len(t))
y, t = simulate_linear_system(sys_cl, u=r, t=t)
fig, ax = plt.subplots(1, 1)
ax.plot(t,y);
ax.grid(which='both')
Now, the weights on the position states are increased
In [32]:
Q[[0, 2], [0, 2]] = [5000, 100]
K, X, eigs = lqr(sys_ss, Q)
print('New Controller gains : ', K)
sys_cl = State(A-B@K, B, C)
y, t = simulate_linear_system(sys_cl, u=r, t=t)
fig, ax = plt.subplots(1, 1)
ax.plot(t,y);
ax.grid(which='both')
There are many other functions already available such as $H_\infty,H_2$ norms, LQR design with output weights, staircase/hessenberg forms, kalman tests and so on. Please don't hesitate to provide feedback by opening issues or submitting PRs to the code repository on Github. That would actually shape where the development time is going to be spent on.
I am especially interested in the pain points of the current software tools for control such that we can redesign the API and make the software more useful.