<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 A Brief Introduction to harold
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.
Step1:  Creating models
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
$$
G_1 = \frac{z-1}{z^3-zs^2+z},\quad 
G_2 = \begin{bmatrix}\frac{1}{s+2}&\frac{s+3}{s^2+s-4}\ 0&\frac{s+2}{s-3}\end{bmatrix},\quad
G_3 = \left[\begin{array}{cc|c}0&1&0\-2&-0.5&1\\hline 0&3.5&1 \end{array}\right],\quad
G_4 = \left[\begin{array}{ccc|cc}-4&-2&2&1&1\0&-5&2&0&1\4&-3&-4&2&-1\\hline -1&5&2&0&0\end{array}\right]
$$
Step2:  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.
Step3:  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}$$
Step4:  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.
Step5:  We could have sliced the array based on the shape of the C array and arrive at the same result
Step6:  Further properties can be accessed via the typical dot syntax, let's check the model data of $G_1$
Step7:  Also the model data of a State representation, for example $G_3$, or just the $A$ matrix of $G_4$
Step8:  Finally we can also create random State models via the random_stable_model
Step9:  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.
Step10:  Conversions of Models
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.
Step11:  Minimal Realization
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.
Step12:  We can see that while on of the $-2$ poles is removed, the others are not cancelled.
Discretization and undiscretization
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.
Step13:  Current known discretization methods are
Step14:  The forward xxx, backward xxx are aliases of &gt;&gt; and &lt;&lt;. 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.
Step15:  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.
Step16:  Model algebra
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.
Step17:  We observe that there is a pole/zero cancellation at 1. which is removed afterwards.
Step18:  Basic Plotting Functionality
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
Step19:  The discrete-time plant plots are automatically drawn as stairs.
Step20:  The plot units can be changed via the keywords available
Step21:  An LQR example
Just to demonstrate the convenience, let's pickup the first section of the LQR example Inverted Pendulum
Step22:  Now, the weights on the position states are increased
<ASSISTANT_TASK:>
Python Code:
import numpy as np
from harold import *
import matplotlib.pyplot as plt
Explanation: A Brief Introduction to harold
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.
End of explanation
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])
Explanation: Creating models
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
$$
G_1 = \frac{z-1}{z^3-zs^2+z},\quad 
G_2 = \begin{bmatrix}\frac{1}{s+2}&\frac{s+3}{s^2+s-4}\ 0&\frac{s+2}{s-3}\end{bmatrix},\quad
G_3 = \left[\begin{array}{cc|c}0&1&0\-2&-0.5&1\\hline 0&3.5&1 \end{array}\right],\quad
G_4 = \left[\begin{array}{ccc|cc}-4&-2&2&1&1\0&-5&2&0&1\4&-3&-4&2&-1\\hline -1&5&2&0&0\end{array}\right]
$$
End of explanation
G1
Explanation: 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.
End of explanation
G5 = Transfer([[[1, 1], 2], [[1, 3], 4]], [1, 5, 1])
print('Numerator :\n', G5.num, '\nDenominator :\n', G5.den)
Explanation: 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}$$
End of explanation
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
Explanation: 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.
End of explanation
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'))
Explanation: We could have sliced the array based on the shape of the C array and arrive at the same result
End of explanation
G1.polynomials
Explanation: Further properties can be accessed via the typical dot syntax, let's check the model data of $G_1$
End of explanation
G3.matrices
G4.a
Explanation: Also the model data of a State representation, for example $G_3$, or just the $A$ matrix of $G_4$
End of explanation
G = random_state_model(5, p=3, m=2)
G
Explanation: Finally we can also create random State models via the random_stable_model
End of explanation
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
Explanation: 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.
End of explanation
H = transfer_to_state(G2)
J = state_to_transfer(G3)
H
J
Explanation: Conversions of Models
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.
End of explanation
Hm = minimal_realization(H)
Hm
Explanation: Minimal Realization
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.
End of explanation
G6 = Transfer([1], [1, 1.4, 1])
G6d = discretize(G6, dt=0.3, method='zoh')
G6d
Explanation: We can see that while on of the $-2$ poles is removed, the others are not cancelled.
Discretization and undiscretization
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.
End of explanation
from harold._global_constants import _KnownDiscretizationMethods as km
print(*km, sep='\n')
Explanation: Current known discretization methods are:
End of explanation
G3.DiscretizedWith is None
G6d.DiscretizedWith
Explanation: The forward xxx, backward xxx are aliases of &gt;&gt; and &lt;&lt;. 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.
End of explanation
G6dcz = undiscretize(G6d)
G6dcz
G6dct = undiscretize(G6d, method='tustin')
G6dct
Explanation: 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.
End of explanation
G1 * G3
G4 @ G2
CL = feedback(G3, G1)
CL_min = minimal_realization(CL)
CL
Explanation: Model algebra
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.
End of explanation
CL_min
Explanation: We observe that there is a pole/zero cancellation at 1. which is removed afterwards.
End of explanation
impulse_response_plot(G4);
Explanation: Basic Plotting Functionality
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
End of explanation
G4_d = discretize(G4, 0.1, method='zoh')
impulse_response_plot(G4_d);
bode_plot(G4);
Explanation: The discrete-time plant plots are automatically drawn as stairs.
End of explanation
bode_plot(G4, use_db=True, use_hz=False, use_degree=False);
nyquist_plot(G2);
Explanation: The plot units can be changed via the keywords available
End of explanation
# 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')
Explanation: An LQR example
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
End of explanation
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')
Explanation: Now, the weights on the position states are increased
End of explanation
<END_TASK>