In [1]:
from controlboros import StateSpaceBuilder
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
In this notebook, we shall simulate the step response of a simple control loop using scipy.signal and controlboros.
The control loop consists of a controller $C(s)$ and a plant $H(s)$ defined by the following transfer functions:
$$ \begin{aligned} C(s) &= \frac{0.1 (5 s + 1)}{s + 1} & H(s) &= \frac{1}{(5 s + 1)(s^2 + 0.2 s + 1)} = \frac{1}{5 s^3 + 2 s^2 + 5.2 s + 1} \end{aligned} $$The closed-loop transfer function is given by:
$$T(s) = \frac{C(s) H(s)}{1 + C(s) H(s)} = \frac{0.1}{s^3 + 1.2 s^2 + 1.2 s + 1.1}$$We want to simulate from 0 to 10 seconds and we use the same coarse sample time (100 ms) for both scipy.signal and controlboros:
In [2]:
t_begin, t_end = 0.0, 10.0
dt = 0.1
Ok. Now we compute the reference step response using scipy.signal:
In [3]:
t_ref, y_ref = signal.step(
([0.1], [1.0, 1.2, 1.2, 1.1]),
T=np.arange(t_begin, t_end, dt),
)
Now we create two controlboros.StateSpace models for the plant and controller using the builder pattern.
Notice that discretise() is just a wrapper around scipy.signal.cont2discrete(). It uses the zero-order hold method per default. You can play around and see the difference e.g. when using the Tustin's approximation (method="bilinear").
In [4]:
ctrl = StateSpaceBuilder().from_tf([0.5, 0.1], [1.0, 1.0])\
.discretise(dt)\
.build()
plant = StateSpaceBuilder().from_tf([1.0], [5.0, 2.0, 5.2, 1.0])\
.discretise(dt)\
.build()
And we initalise variables and arrays for our signals:
In [5]:
t_cb = np.arange(t_begin, t_end, dt)
setpoint = np.array([1.0]) # Setpoint, constantly 1.0 because step response
command = np.array([0.0]) # Controller output
feedback = np.array([0.0]) # Feedback value
y_cb = np.zeros((len(t_cb),)) # Array for the step response
We're ready to run the main loop. Notice how we resolve the control loop by using a unit delay:
In [6]:
# Reset the inital state of the systems,
# helpful if you run this cell multiple times!
ctrl.set_state_to_zero()
plant.set_state_to_zero()
for i in range(len(t_cb)):
command = ctrl.push_stateful(setpoint - feedback)
y_cb[i] = plant.push_stateful(command)
feedback = y_cb[i] # Unit delay!
Plot and compare results:
In [7]:
plt.figure(figsize=(7, 5))
plt.plot(t_ref, y_ref)
plt.step(t_cb, y_cb, where="post")
plt.xlabel("Time (s)")
plt.ylabel("Output")
plt.legend(["scipy.signal", "controlboros"])
plt.grid()
plt.show()
Obviously, the step response simulated with controlboros is very inaccurate due to 100 ms time delay in the feedback loop. We can alleviate (but never solve!) this problem by using a finer sample time, e.g. 1 ms:
In [8]:
dt_fine = 1.0e-3
ctrl_fine = StateSpaceBuilder().from_tf([0.5, 0.1], [1.0, 1.0])\
.discretise(dt_fine)\
.build()
plant_fine = StateSpaceBuilder().from_tf([1.0], [5.0, 2.0, 5.2, 1.0])\
.discretise(dt_fine)\
.build()
t_cb_fine = np.arange(t_begin, t_end, dt_fine)
setpoint = np.array([1.0]) # Setpoint, constantly 1.0 because step response
command = np.array([0.0]) # Controller output
feedback = np.array([0.0]) # Feedback value
y_cb_fine = np.zeros((len(t_cb_fine),)) # Array for the step response
# Reset the inital state of the systems
ctrl_fine.set_state_to_zero()
plant_fine.set_state_to_zero()
for i in range(len(t_cb_fine)):
command = ctrl_fine.push_stateful(setpoint - feedback)
y_cb_fine[i] = plant_fine.push_stateful(command)
feedback = y_cb_fine[i]
We plot the results and see that the controlboros solution with 1 ms step size is very close to the one computed with scipy.signal:
In [9]:
plt.figure(figsize=(7, 5))
plt.plot(t_ref, y_ref)
plt.step(t_cb, y_cb, where="post")
plt.step(t_cb_fine, y_cb_fine, "r", where="post")
plt.xlabel("Time (s)")
plt.ylabel("Output")
plt.legend(["scipy.signal", "controlboros, 100 ms", "controlboros, 1 ms"])
plt.grid()
plt.show()