Category: Engineering

LED Dimmer with Raspberry Pi

In this simple hardware setup, we will explore different coding possibilities using a Raspberry Pi with only a push button and an LED. It will be a good opportunity to compare Python classes vs functions, with code that shows the benefits of the former. More specifically, classes offer modularity and the possibility of reuse through multiple instances.

The diagram shows the very straightforward setup where the button and LED have one of their legs connected to a (distinct) GPIO pin and the other to a ground pin. A 330 Ω resistor must be connected in series with the LED, whose intensity is controlled by a PWM output. That is as simple as it gets.

Note: If the button doesn’t seem to be working, try connecting a different pair of legs, since some of them could be common.

While the dimmer functionality could be achieved using solely electronic circuitry and without a Raspberry Pi, the main point of this post is to explore coding and see what can be done with software!

Let’s start with the requirements of how we want the dimmer to work. By pressing and holding the button, the LED light intensity should increase linearly until reaching its maximum value and then decrease until it reaches its minimum value (LED off). This cycle repeats itself until the button is released, causing the light intensity to stay constant at its last value. By pressing the button again, the cycle resumes from where it stopped.

The plots provide a visual understanding of the requirements. The PWM output controlling the LED intensity ramps up (green) and down (red) while the button is pressed. As it is released, the last PWM output level is then held until the button is pressed again.

In terms of control logic, I use an integer to count ramp transitions. Even values correspond to upward ramps and odd values to downward ramps. Every time the button is released and then pressed again, the counter is reset to zero.

There are two main components in the control logic for the dimmer:

• The `calc_ramp()` function that creates the triangular wave output in three steps. 1) Linear output generation, 2) Converting linear output to a tooth saw wave, thus limiting it between 0 and 1, and 3) Flipping odd count sections to generate a triangular wave.
• Time shifting and counter reset to ensure smooth transitions every time the button is released and then pressed again. When a ramp is interrupted due to a button release event, it it will restart from where it stopped so the effective ramp duration (discounting the elapsed released time) is always the same.

Python Function

The code below shows how the function `calc_ramp()` and the time shifting are integrated inside an execution loop that controls the LED light intensity. As usual, the interface with the Raspberry Pi is done using GPIO Zero classes, more specifically `DigitalInputDevice` for the button pin and `PWMOutputDevice` for the LED pin.

```# Importing modules and classes
import time
import numpy as np
from gpiozero import DigitalInputDevice, PWMOutputDevice
from utils import plot_line

# Assigning dimmer parameter values
pinled = 17  # PWM output (LED input) pin
pinbutton = 5  # button input pin
pwmfreq = 200  # PWM frequency (Hz)
pressedbefore = False  # Previous button pressed state
valueprev = 0  # Previous PWM value
kprev = 0  # Previous PWM ramp segmment counter
tshift = 0  # PWM ramp time shift (to start where it left off)
tramp = 2  # PWM output 0 to 100% ramp time (s)
# Assigning execution loop parameter values
tsample = 0.02  # Sampling period for code execution (s)
tstop = 30  # Total execution time (s)

# Preallocating output arrays for plotting
t = [] # Time (s)
value = [] # PWM output duty cycle value
k = [] # Ramp segment counter

def calc_ramp(t, tramp):
"""
Creates triangular wave output with amplitude 0 to 1 and period 2 x tramp.

"""
# Creating linear output so the value is 1 when t=tramp
valuelin = t/tramp
# Creating time segment counter
k = t//tramp
# Shifting output down by number of segment counts
value = valuelin - k
# Flipping odd count output to create triangular wave
if (k % 2) == 1:
value = 1 - value
return value, k

# Creating button and PWM output (LED input) objects
button = DigitalInputDevice(pinbutton, pull_up=True, bounce_time=0.1)
led = PWMOutputDevice(pinled, frequency=pwmfreq)

# Initializing timers and starting main clock
tpressed = 0
tprev = 0
tcurr = 0
tstart = time.perf_counter()
# Initializing PWM value and ramp time segment counter
valuecurr = 0
kcurr = -1

# Executing loop
print('Running code for', tstop, 'seconds ...')
while tcurr <= tstop:
# Executing code every `tsample` seconds
if (np.floor(tcurr/tsample) - np.floor(tprev/tsample)) == 1:
# Getting button properties only once
pressed = button.is_active
# Checking for button press
if pressed and not pressedbefore:
# Calculating ramp time shift based on last PWM value
if (kprev % 2) == 0:
tshift = valueprev*tramp
else:
tshift = (1-valueprev)*tramp + tramp
# Starting pressed button timer
tpressed = time.perf_counter() - tstart
# Updating previous button state
pressedbefore = True
# Checking for button release
if not pressed and pressedbefore:
# Storing PWM value and ramp segment counter
valueprev = led.value
kprev = kcurr
# Updating previous button state
pressedbefore = False
# Updating PWM output (LED intensity)
if pressed:
valuecurr, kcurr = calc_ramp(tcurr-tpressed+tshift, tramp)
led.value = valuecurr
# Appending current values to output arrays
t.append(tcurr)
value.append(valuecurr)
k.append(kcurr)
# Updating previous time and getting new current time (s)
tprev = tcurr
tcurr = time.perf_counter() - tstart
print('Done.')
# Releasing pins
led.close()
button.close()

# Plotting results
plot_line([t]*2, [value, k],
yname=['PWM Output', 'Segment Counter'], axes='multi')
plot_line([t[1::]], [1000*np.diff(t)], yname=['Sampling Period (ms)'])
```

Python Class

Having the function and the time shifting “mixed” with the code makes for a less clean execution loop section, which can get pretty busy if there are more devices being controlled inside of it. The next piece of code shows how to integrate the two main components of the dimmer control inside a Python class. In addition to removing complexity from the execution loop, it makes it really easy to run multiple dimmers at the same time! This circles back to the modularity that was mentioned at the top of the post. You can create different classes for more complex devices and have them all be part of a module, such as gpiozero_extended.py. In the example below, I create two LED dimmer objects that can be controlled independently.

```# Importing modules and classes
import time
import numpy as np
from gpiozero import DigitalInputDevice, PWMOutputDevice

class Dimmer:
"""
Class that represents an LED dimmer.

"""
def __init__(self, pinbutton=None, pinled=None):
# Assigning GPIO pins
self._pinbutton = pinbutton  # button input pin
self._pinled = pinled  # PWM output (LED input) pin
# Assigning dimmer parameter values
self._pwmfreq = 200  # PWM frequency (Hz)
self._tshift = 0  # PWM ramp time shift (to start where it left off)
self._tramp = 2  # PWM output 0 to 100% ramp time (s)
self._tpressed = 0  # Time button was pressed (s)
self._tstart = 0  # Starting time of dimmer execution
self._pressed = False  # Current button pressed state
self._pressedbefore = False  # Previous button pressed state
self._valueprev = 0  # Previous PWM value
self._kcurr = -1  # Current PWM ramp segment counter
self._kprev = 0  # Previous PWM ramp segmment counter
# Creating button and PWM output (LED input) objects
self._button = DigitalInputDevice(self._pinbutton, pull_up=True, bounce_time=0.1)
self._led = PWMOutputDevice(self._pinled, frequency=self._pwmfreq)

def _calc_ramp_value(self, t):
# Creating linear output so the value is 1 when t=tramp
valuelin = t/self._tramp
# Creating time segment counter
self._kcurr = t//self._tramp
# Shifting output down by number of segment counts
value = valuelin - self._kcurr
# Flipping odd count output to create triangular wave
if (self._kcurr % 2) == 1:
value = 1 - value
return value

def reset_timer(self, tstart):
# Resets dimmer timer to `tstart`
self._tstart = tstart

def update_value(self, t):
# Getting button properties only once
self._pressed = self._button.is_active
# Checking for button press
if self._pressed and not self._pressedbefore:
# Calculating ramp time shift based on last PWM value
if (self._kprev % 2) == 0:
self._tshift = self._valueprev*self._tramp
else:
self._tshift = (1-self._valueprev)*self._tramp + self._tramp
# Starting pressed button timer
self._tpressed = time.perf_counter() - self._tstart
# Updating previous button state
self._pressedbefore = True
# Checking for button release
if not self._pressed and self._pressedbefore:
# Storing PWM value and ramp segment counter
self._valueprev = self._led.value
self._kprev = self._kcurr
# Updating previous button state
self._pressedbefore = False
# Updating PWM output (LED intensity)
if self._pressed:
self._led.value = self._calc_ramp_value(t-self._tpressed+self._tshift)

def __del__(self):
self._button.close()
self._led.close()

# Assigning execution loop parameter values
tsample = 0.02 # Sampling period for code execution (s)
tstop = 30 # Total execution time (s)

# Creating dimmer objects
dimmer1 = Dimmer(pinled=17, pinbutton=5)
dimmer2 = Dimmer(pinled=18, pinbutton=6)

# Initializing timers and starting main clock
tprev = 0
tcurr = 0
tstart = time.perf_counter()
# Resettting dimmer timer
dimmer1.reset_timer(tstart)
dimmer2.reset_timer(tstart)

# Executing loop
print('Running code for', tstop, 'seconds ...')
while tcurr <= tstop:
# Executing code every `tsample` seconds
if (np.floor(tcurr/tsample) - np.floor(tprev/tsample)) == 1:
# Updating dimmer outputs
dimmer1.update_value(tcurr)
dimmer2.update_value(tcurr)
# Updating previous time and getting new current time (s)
tprev = tcurr
tcurr = time.perf_counter() - tstart
print('Done.')

# Releasing pins
del dimmer1
del dimmer2
```

If you watch the video, you’ll notice that the light intensity change is not as linear as one would expect, even when the PWM output is. One way to improve this is to add a non-linear transfer function between the ramp function output and the actual PWM set point value.

Also, it would be nice to have the ability to switch the dimmer off by, for instance, pressing and releasing the button for a short period of time. The improved class that does that can be found on my GitHub page.

Curve Fitting with Tangent and Inverse Tangent

The inverse tangent (arctangent) and tangent functions are quite helpful when it comes to fitting data that characterize the response of some types of sensors. Sensors that transition between two saturation states can have their transfer function (relationship between input and output) represented by the arctangent. The tangent function, on the other hand, comes into play when we need to fit the inverse of the sensor transfer function, with the purpose of linearizing its output in the control software.

Let’s use the Python SciPy package for optimization and root finding, more specifically, the `curve_fit` function. which is well suited for the case where nonlinear coefficients need to be determined. Convenient forms for data fitting of arctangent and tangent are shown below:

The coefficientsthat produce the best data fit in a least squares sense can be determined using `curve_fit`. It should be noted that, if you invert the tangent function above, it’s fairly easy to arrive at:

While the fitting problem using the arctangent and tangent functions can be treated completely independently, in some cases, one function does a better job fitting the data than the other one. The coefficient relationships above can be used to go from tangent to arctangent (or vice-versa), when the inverse of the desired function produces a better fit to the (inverted) data than the desired function itself.

Arctangent Function Fit

The following Python code fits data from a photocell sensor as described in Line Tracking Sensor for Raspberry Pi. `xdata` in the code represents the sensor position in mm, while `ydata` represents the photocell normalized voltage output.

```# Importing modules and classes
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# Defining function to generate Matplotlib figure with axes
def make_fig():
#
# Creating figure
fig = plt.figure(
figsize=(5, 4.2),
facecolor='#ffffff',
tight_layout=True)
facecolor='#ffffff',
)
ax.grid(
linestyle=':',
)
# Returning axes handle
return ax

# Defining parametrized arctangent function
def f(x, k, w, x0, y0):
return k * np.arctan(w*(x-x0)) + y0

# Defining data points
xdata = [0, 2, 4, 6, 8, 10, 12, 14, 16]
ydata = [0.432, 0.448, 0.477, 0.547, 0.653, 0.800, 0.908, 0.933, 0.955]

# Estimating initial parameter values
p0 = [1, 1]
p0.append(0.5*(max(xdata)+min(xdata)))
p0.append(0.5*(max(ydata)+min(ydata)))
# Fitting data
output = curve_fit(f, xdata, ydata, p0=p0, full_output=True)
# Extracting relevant information from output
copt = output[0]
res = output[2]['fvec']
numeval = output[2]['nfev']
msg = output[3]

# Plotting data and fitted function
xi = np.linspace(min(xdata), max(xdata), 100)
yi = f(xi, *copt)
ax = make_fig()
ax.set_xlabel('x', fontsize=12)
ax.set_ylabel('y', fontsize=12)
ax.scatter(xdata, ydata)
ax.plot(xi, yi, linestyle=':')

# Calculating residual sum of squares (RSS)
# Calculating root mean squared error (RMSE) from RSS
# Calculating R-square
# Displaying fit info
print('Found solution in {:d} function evaluations.'.format(numeval))
print(msg)
print('R-square = {:1.4f}, RMSE = {:1.4f}'.format(r2, rmse))
print('')
print('Coefficients:')
print('  k  = {:1.3f}'.format(copt[0]))
print('  w  = {:1.3f}'.format(copt[1]))
print('  x0 = {:1.3f}'.format(copt[2]))
print('  y0 = {:1.3f}'.format(copt[3]))
```

Running the program will generate a plot of the data points with the fitted curve, as well as some statistical information about the fit quality and the optimal coefficients.

```Found solution in 43 function evaluations.
Both actual and predicted relative reductions
in the sum of squares are at most 0.000000

R-square = 0.9990, RMSE = 0.0065

Coefficients:
k  = 0.211
w  = 0.377
x0 = 8.565
y0 = 0.701

```

Tangent Function Fit

The same code can be used to fit a tangent function to a data set. All you need to do is replace the arctangent function with the tangent one and have the appropriate `xdata` and `ydata`. In this particular case, the previous data set arrays were just swapped. The program output after the modifications is also shown. The full code can be found on my GitHub repository.

```# Defining parametrized tangent function
def f(x, k, w, x0, y0):
return k * np.tan(w*(x-x0)) + y0

# Defining data points
xdata = [0.432, 0.448, 0.477, 0.547, 0.653, 0.800, 0.908, 0.933, 0.955]
ydata = [0, 2, 4, 6, 8, 10, 12, 14, 16]
```
```Found solution in 54 function evaluations.
Both actual and predicted relative reductions
in the sum of squares are at most 0.000000

R-square = 0.9985, RMSE = 0.1998

Coefficients:
k  = 2.383
w  = 4.896
x0 = 0.696
y0 = 8.377

```

If we compare the coefficients from the arctangent fit with the ones from the tangent fit, using the relationships derived at the top of the post, they don’t quite match. Even when using the exact same data set by swapping `xdata` and `ydata`.

This result shouldn’t be a surprise, since the error between y values and their predictions, used inside the fitting algorithm, are not equivalent even with identical swapped data sets. Additionally, the distribution of the data points among the x coordinates also affects the optimal coefficient values of the fitted function.

Going from Arctangent to Tangent

Let’s take a closer look at the influence of the data point distribution on the resulting fitted function, and when we are better off using the coefficients of the function fitted to the inverse data set.

We want to fit a transfer function to the following data using the tangent function:

```# Defining parametrized tangent function
def f(x, k, w, x0, y0):
return k * np.tan(w*(x-x0)) + y0

# Defining data points
xdata = [0.700, 0.698, 0.691, 0.680, 0.654, 0.587, 0.504, 0.394, 0.341, 0.300, 0.282, 0.275, 0.272]
ydata = [0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24]
```

Using the Python program with the code snippet above substituted into it, we get the fit results below. Even though the R-square is a respectable value, the fitted function doesn’t do a good job capturing the transition between the two end states.

```Found solution in 80 function evaluations.
Both actual and predicted relative reductions
in the sum of squares are at most 0.000000

R-square = 0.9959, RMSE = 0.4766

Coefficients:
k  = -2.143
w  = 6.469
x0 = 0.485
y0 = 12.382

```

Due to the overall shape of the transfer function and the fact that the data was originally collected at constant y increments, there is a high concentration of points along the x axis on both ends of the test data range. That leaves the region of interest (inside the red ellipse) with relatively fewer points, biasing the fitted function towards the low and high x value regions (outside the red ellipse).

One way to tackle this issue is by inverting the data set and using, in this case, the arctangent function:

```# Defining parametrized arctangent function
def f(x, k, w, x0, y0):
return k * np.arctan(w*(x-x0)) + y0

# Defining data points
xdata = [0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24]
ydata = [0.700, 0.698, 0.691, 0.680, 0.654, 0.587, 0.504, 0.394, 0.341, 0.300, 0.282, 0.275, 0.272]
```

The new fit is now done with x values that are evenly distributed, removing the bias that occurred when the tangent function was used.

```Found solution in 60 function evaluations.
Both actual and predicted relative reductions
in the sum of squares are at most 0.000000

R-square = 0.9992, RMSE = 0.0049

Coefficients:
k  = -0.164
w  = 0.354
x0 = 12.253
y0 = 0.486

```

Finally, we can then use the coefficient relationships to go from arctangent, that produces a better overall data fit, to the desired tangent transfer function. The plot below shows the new tangent function in orange (using the arctangent coefficients) in contrast with the original one in blue.

Line Tracking Sensor for Raspberry Pi

Robots that can follow a line on the floor (more specifically AGVs – Automated Guided Vehicles) have been around for quite some time now, moving around in factories and warehouses. The most basic ones follow a magnetic embedded tape or an actual line painted on the floor. In this post, we will take a closer look at a sensor whose output is a continuous value that can be used with the Raspberry Pi. There are quite a few so-called line tracking sensors available for the Pi. However, most of them provide a digital output, i.e., they detect whether you are either on (one) or off (zero) the line, putting the burden on the control system that has to deal with only two possible sensor states. A sensor with a continuous output is much more suitable for what is primarily a set-point tracking control problem.

The pictorial representation on the left shows the two electronic components that were used to make the sensor: a white LED (Keyestudio KS0016) and a photocell (KS0028). The former provides a consistent light source and the latter senses the reflected light intensity off the surface the sensor is pointing at.

In you are mounting the two components to a metal bracket, it is VERY IMPORTANT to use an insulator as shown in the picture! Otherwise, you will short the soldered tips on the circuit boards and permanently damage your Raspberry Pi. I used the paper from a cereal box, folded to the desired thickness.

The wiring schematic below shows how the LED and the photocell are integrated with the Raspberry Pi. Since the photocell output is continuous, an MCP3008 analog-to-digital converter chip is also required.

The connection of the LED to the Pi is pretty straightforward, with the signal pin connected to a GPIO pin of your choosing. In my case, I went with pin 18.

Taking a closer look at the MCP3008 connections, I used 3.3V as the reference voltage for the chip (instead of 5V) for a better absolute resolution.

Since both components require a 5V supply voltage, for the KS0028, a voltage divider must be used to bring the maximum photocell output down to 3.3V (or less) before connecting it to the MCP3008. After doing some testing, a 1 kΩ / 4 kΩ split gave me the ratio that best utilized the 3.3V range.

Sensor Characterization

The actual sensor setup is shown next, where I used black electrical tape (approx. 18 mm wide) as the line to be followed. The test procedure for the sensor characterization is illustrated on the right at its initial position (marker 0 on the actual sensor picture) and final position (after moving the bracket in 2 mm steps until reaching 18 mm). As a side note, the LED tip sits at approximately 5 mm above the surface.

The table below shows the normalized sensor output (where 1 would be the reference voltage of 3.3V) for three sets of measurements at 2 mm position increments. Observe that the row at 18 mm was omitted, since its values were basically the same as the row at 16 mm. That means that the sensor was already being exposed to the maximum amount of reflected light off the white surface at 16 mm. Additionally, the sensor placed fully above the black line (position 0) gives the lowest output for the reflected light, which is about 0.43.

The plot on the left shows the bracket position as a function of the sensor output and the fitted cubic polynomial. The plot on the right shows the transfer function that characterizes the sensor, which is basically the fitted polynomial with the appropriate shift and offset.

A few notes on the transfer function:

• Technically speaking, the sensor transfer function provides the relationship between the physical input (position, in this case) and the sensor output (normalized voltage). In this particular post, I will be calling that one “sensor output function“. The relationship between sensor output and physical input that is to be implemented in code will then be called transfer function.
• The variability on the sensor output is more likely due to the variability of the position, since eyeballing the exact 2 mm steps when moving the bracket has plenty of uncertainty to it.
• The transfer function is basically the inverse of the physical sensor output function and would reside inside the software bridging the real world and the control system. Moreover, it linearizes the sensor behavior! For example, if the bracket is at the position 2 mm, the sensor output would be 0.85. However, once that measurement is brought into the control software and input into the transfer function, the output is an estimation of the original position of 2 mm.

The Python code below was used to generate the previous plots and determine the coefficients for the transfer function. The next section shows some polynomial trickery used to calculate the shift and offset for the fitted cubic function, centered around the position at which the LED is right on the transition between the white and black surfaces (the inflection at 0.7).

```# Importing modules and classes
import numpy as np
import matplotlib.pyplot as plt
from numpy.polynomial import Polynomial as P

# Defining function to generate Matplotlib figure with axes
def make_fig():
#
# Creating figure
fig = plt.figure(
figsize=(5, 4.2),
facecolor='#ffffff',
tight_layout=True)
facecolor='#ffffff',
)
ax.grid(
linestyle=':',
)
# Returning axes handle
return ax

# Defining data values u (sensor output) and s (bracket position)
# from sensor characterization test
u = [0.432, 0.448, 0.477, 0.547, 0.653, 0.800, 0.908, 0.933, 0.955]
s = [0, 2, 4, 6, 8, 10, 12, 14, 16]

# Fitting 3rd order polynomial to the data
p, stats = P.fit(u, s, 3, full=True)
# Evaluating polynomial for plotting
ui = np.linspace(0.45, 0.95, 100)
si = p(ui)
# Plotting data and fitted polynomial
ax = make_fig()
ax.set_xlabel('Sensor Output (u)', fontsize=12)
ax.set_ylabel('Position (mm)', fontsize=12)
ax.scatter(u, s)
ax.plot(ui, si, linestyle=':')

# Getting polynomial coefficients
c = p.convert().coef
# Calculating coefficients for transfer function
# (shift and offset representation)
k = c[3]
u0 = -c[2]/(3*k)
w = c[1] - 3*k*u0**2
s0 = c[0] + k*u0**3 + w*u0
# Defining transfer function for plotting
stf = k*(ui-u0)**3 + w*(ui-u0)
# Plotting transfer function
ax = make_fig()
ax.set_xlabel('Sensor Output (u)', fontsize=12)
ax.set_ylabel('Transfer Function Output (mm)', fontsize=12)
ax.plot(ui, stf)

# Calculating root mean squared error from residual sum of squares
# Calculating R-square
# Displaying fit info
print('R-square = {:1.3f}'.format(r2))
print('RMSE (mm) = {:1.3f}'.format(rmse))
print('')
print('Shift (-) = {:1.3f}, Offset (mm) = {:1.3f}'.format(u0, s0))
```

Polynomial Trickery

Let’s see how we can determine the shift and offset of a third order polynomial that is expressed in its most general form:

Due to the nature of the sensor response, it is reasonable to assume that the corresponding polynomial in its shift and offset form is odd (other than the offset term) and can be expressed as:

Both representations have 4 coefficients that need to be determined (or). The polynomial fitting gave us the first set, so we need to determineas a function of. To do so, let’s expand the second polynomial above and regroup its terms in descending powers of x:

Finally, all we have to do is equate (solving for the variables of interest) the corresponding coefficients of the polynomial above with the one at the top of this section:

Which is exactly what’s in the Python code that produced the transfer function for the line tracking sensor.

Sensor Python Class

From a control software standpoint, the best way to use the line tracking sensor with a robot is to package its transfer function inside a class. That way, it’s much easier to get the sensor position reading within a control loop. The Python code below contains the `LineSensor` class, in which the transfer function can be easily identified. The motor speed control post can be a starting point for a closed-loop control using the `PID` class. That one and a more complete version of `LineSensor` can be found in `gpiozero_extended.py`.

```# Importing modules and classes
import time
import numpy as np
from gpiozero import MCP3008, DigitalOutputDevice

class LineSensor:
"""
Class that implements a line tracking sensor.

"""
def __init__(self, lightsource, photosensor):
# Defining transfer function parameters
self._k = 335.5
self._w = 5.821
self._u0 = 0.700
# Creating GPIOZero objects
self._lightsource = DigitalOutputDevice(lightsource)
self._photosensor = MCP3008(
channel=photosensor,
clock_pin=11, mosi_pin=10, miso_pin=9, select_pin=8)
# Turning light source on
self._lightsource.value = 1

@property
def position(self):
# Getting sensor output
u = self._photosensor.value
# Calculating position using transfer function
return self._k*(u-self._u0)**3 + self._w*(u-self._u0)

@position.setter
def position(self, _):
print('"position" is a read only attribute.')

def __del__(self):
self._lightsource.close()
self._photosensor.close()

# Assigning some parameters
tsample = 0.1  # Sampling period for code execution (s)
tdisp = 0.5  # Output display period (s)
tstop = 60  # Total execution time (s)
# Creating line tracking sensor object on GPIO pin 18 for
# light source and MCP3008 channel 0 for photocell output
linesensor = LineSensor(lightsource=18, photosensor=0)

# Initializing variables and starting main clock
tprev = 0
tcurr = 0
tstart = time.perf_counter()
# Execution loop
print('Running code for', tstop, 'seconds ...')
while tcurr <= tstop:
# Getting current time (s)
tcurr = time.perf_counter() - tstart
# Doing I/O and computations every `tsample` seconds
if (np.floor(tcurr/tsample) - np.floor(tprev/tsample)) == 1:
# Getting sensor position
poscurr = linesensor.position
#
# Insert control action here using gpiozero_extended.PID()
#
# Displaying sensor position every `tdisp` seconds
if (np.floor(tcurr/tdisp) - np.floor(tprev/tdisp)) == 1:
print("Position = {:0.1f} mm".format(poscurr))
# Updating previous time value
tprev = tcurr

print('Done.')
# Deleting sensor
del linesensor
```

Final Remarks

First, it’s important to check how well the transfer function estimates the actual position of the sensor. That can be done by running the code from the previous section placing the bracket at the original test positions.

The graph on the left shows measured values (from the code output) vs. desired values (by moving the bracket onto each test marker). The dashed line represents the 45-degree line on which measured values match actual ones. You can find the code for it here.

The error is more accentuated away from the center, which seems to be where the cubic fit deviates more from the test points.

Second, the ambient light intensity does affect the sensor reading and therefore will compromise the validity of the transfer function. It is important then to determine the transfer function with lighting conditions that represent where the robot will be running. Preferably, with the sensor mounted on the robot!

Finally, for this specific sensor, the width of the line should be about 20 mm. If the line is too thin, or if the robot is moving too fast for the turn it’s supposed to make, the position deviation may exceed the line thickness. In that case, the sensor would again be “seeing” a white surface on the wrong side of the line(!), causing the negative feedback loop of the control system to become a positive one, moving the robot further away from the line.

DC Motor Characterization (2 of 2)

Picking it up from where we stopped in the previous post, let’s see how we can determine the effective moment of inertia, the effective damping and the torque constant for the reduced order model shown in the diagram.

Before we start, it’s a good idea to revisit a typical set of torque/speed curves for a DC motor. They will come in handy to better understand the parameter determination (identification) tests that are explained further down the post.

The graph shows torque/speed curves at constant input voltages for a DC motor, for arbitrarily chosen levels of 40, 60, 80 and 100%.

The torque values on the zero-speed y-axis represent the motor stall torque at different input voltages. Conversely, the speed values on the zero-torque x-axis represent the unloaded maximum motor speed, also at different input voltages.

These two special sets of data points will be used to determineandrespectively.

This time around, I’ll use a LEGO Mindstorms EV3 large motor. The basic reason behind such hardware choice is the easiness to build different test fixtures based on the parameter that’s being identified.

Determination of DC Motor Torque Constant

The torque constantcan be determined through a series of torque stall tests with varying motor input levels. The images below show the LEGO fixture that was built to accomplish the task.

The fixture design allows for the motor angular motion to be converted to a linear motion through the pulley system and the Nylon line. Therefore, the torque can be calculated from the measured linear force using the spring gauge.

Since the radius of the large pulley is known, the measured force in N can be converted to a measured torque in N.m.

The low budget gauge accuracy was verified with known masses and proved to be meet the standards of this blog. As a side note, calibration/verification of your measurement gauges should always be a part of any experimental testing!

Because of the nature of the force gauge, where the measurement readings had to be done visually by myself, it wasn’t possible to automate this test. Nonetheless, I used the pyev3 package to create a motor object and interactively set output levels.

It should be noted that the motor input can range between -100 and 100 %. The values are likely PWM duty cycle and seem to have a linear relationship with armature voltage (to which I don’t have access through the EV3 brick).

Finally, the mass of the spring gauge was taking into account when calculating the torque values.

The stall torque was measured for 4 large EV3 motors, where each point on the plot is the average of three torque values calculated from the corresponding force reading.

Four input levels were used, i.e., 40, 60, 80, and 100%. The dashed lines were fitted to the data points for each motor. Now, if we remember from the first test of the previous post that:

Then the slope of each line represents the torque constant in N.m/%. The average value is:

Determination of DC Motor Damping

Similarly to the torque constant determination, four large EV3 motors were used to generate the damping data. This time around however, because no torque load is applied to the motors, no fixture was required for the testing other than my hand holding the motor.

The plot shows the results for the zero-load tests, where each point is the average of three measurements for a given motor. The same input levels of the previous section were used.

If we reference back to the second test of the DC Motor Characterization (1 of 2) post, where:

The slope of the fitted dashed lines represents the speed constant in rad/s/%. The average value is:

Since we already know the value of the torque constant, the damping constant can then be calculated, leading to:

The Python code used to generate the plots and data fits above, as well as the automated test for the determination of the motor speed constant, can be found here.

Determination of DC Motor Inertia

The final parameter to be identified is the effective moment of inertia. As the name suggests, it accounts for the rotor inertia and any rotating masses attached to the output shaft. In the particular case of the EV3 motor, there’s a gearbox between the rotor and the output shaft. Therefore the moment of inertia of all the gears and the effect of the gear ratios will be accounted for, once we determine the effective moment of inertiareduced to the output shaft.

The moment of inertia will be estimated using the third test from the previous post as illustrated on the plot to the left. The idea is to determine the time constantof the approximated first order system and use

to determine since is already known.is the time to reach 63.2% of the steady-state speed value.

Some signal processing was used, as suggested in the graph above. A zero-shift low-pass filter was applied to the motor angular speed signal and thenwas calculated as the median value of the last 0.5 seconds of the signal. With the steady-state speed value in hands, can be easily determined and therefore .

Four different motors were tested at four input levels (40, 60, 80, and 100%), where five measurements were taken at each level. With a median time constantof 0.076 s, and using the the formula above with the already determined , The effective moment of inertia at the motor shaft output is:

Note that the effective moment of inertia is a fairly large number for such a small motor. Remember however that the moment of inertia at the shaft output is proportional to the squared value of the gear ratio. Which means, for our particular LEGO motor, if the gear ratio is for example 1:30, the rotor moment of inertia is 900 times smaller that the value at the gearbox output!

More on Effective Moment of Inertia

Let’s explore what happens if we attach an additional rotating mass to the shaft output, thus modifying the effective moment of inertia. One would expect the response time constant to increase as the rotating inertia increases.

To that end, the original test fixture was modified to accommodate two balanced counter weights. Such setup creates a zero load torque while allowing for different load inertias (by changing the mass of the counter weights).

Just for the sake of this exercise, let’s call the original effective moment of inertia of the LEGO large motor at the output shaft asand the additional inertia due to the counter-weights as . Note that

whereis the total mass of the counter-weights andis the radius of the pulley. The moment of inertia of the pulley is negligible when compared to that of the counter-weights.

Two sets of weights were used, i.e., 0.149 and 0.249 kg with a corresponding moment of inertia of 0.00023 and 0.00045 kg.m2 (the pulley radius is 39.25 mm). Using the formula below with all the known quantities allows us to calculate the time constant of the system with motor plus counter-weights.

In our case, the calculated values ofare 0.086 and 0.098 s. Running the test shown on the left, with four different motors and four different output levels, in fact produced measured values within +/- 3% of the corresponding calculated ones!

DC Motor Model Confirmation

It’s always a good idea to check whether the motor characterization parameters can be used to build a model that represents the actual system behavior. To do that, a set of counter-weight masses that haven’t been tested during the characterization was chosen, resulting in the following motor model.

Where the load torque due to the uneven counter-weight masses is 0.055 N.m and the effective inertia is 0.0014 + 0.00034 = 0.00174 kg.m2. The plot on the left shows good agreement between the actual LEGO fixture test results and the Simulink simulation of the reduced-order model above, for a 70% duty-cycle step input.

DC Motor Characterization (1 of 2)

In this first of two posts, we’ll go over the simplifications of the equations that govern the dynamic behavior of a DC motor. By doing so, we should be able to use a non-intrusive and fairly straightforward approach to determine the parameters that characterize the motor from a control systems stand point.

While this might get a bit arid at times, I’ll try to keep it interesting by using some real motor parameters and Simulink to show the impact of the simplifications as we move towards our final dynamic system equation. Let’s start with the well known coupled differential equations (which are essentially a mathematical representation) for the electric motor:

Where (in SI units, for good measure):

By applying the Laplace transformation to the equations above and rearranging them in a way that they are solved for the angular velocity and current, we get:

and therefore

The rearranged equations can then be represented in the form of a block diagram, which will come in handy when we talk about the model simplification. Also, note the negative “feedback” nature of the back EMF term(which makes the rotor speed inherently stable to small disturbances in the load torque).

The DC motor has two time constants associated with its dynamic behavior. The electrical time constant(directly obtained from one of the first-order blocks in the diagram) and the mechanical time constant, which is not so intuitive and we’ll see later where it comes from.

Especially for smaller motors, the electrical time constant is substantially faster than the mechanical time constant. That is generally accompanied by, which means that the corresponding inductance term in the first-order block can be neglected. The block diagram of the resulting system is shown next.

It’s time to put some numbers into these two models and observe how they behave when simulated using Simulink. In order to do that, I found myself a motor data sheet for permanent magnet DC motors from Moog. These tend to be very high end components when compared to the DC motors I’ve been using in my previous posts. I.e. they are designed for a more professional application. Not for the “DIY-finite-budget-hobbyist”, like myself!

With that said, I chose the smallest and the largest motors from the catalog, whose parameters are summarized below:

Taking a closer look at the last two rows, the electrical time constant of the smaller motor is significantly faster than its mechanical time constant, especially when compared to the larger motor. We will see how that plays out in a bit.

Also, for the simulation, the two inputs are the armature voltage and the load (or disturbance) torque, the latter accounting for the friction torque as well. The outputs are shown in the Simulink model picture below. The armature current and angular speed are the ones of most interest for checking the validity of the model against the data sheet values. In other words, when running the model with only voltage and torque as inputs, everything else should “fall in place”, matching the table data.

Small DC Motor Simulation Results

The plots show the simulation results for a step response with the armature voltage (12 V) and rated torque (0.05 N.m) for the small C23-L33-W10 motor. Two Simulink models were used: the complete model (shown above) and the simplified one (not shown here). Both can be found on my GitHub page.

The first takeaway is that the final rotor speed and armature current agree well with the corresponding rated values in the table. Also, the response time of 25.6 ms correlates fairly well with what is, in this particular case, the mechanical response time in the table.

While the simulated steady-state rotor speed is very close to the table value (502 vs 492 rpm), the armature current correlation is not so stellar (4.0 vs. 4.75 A) yet still satisfactory. It should be noted that the table values are probably obtained from actual motor tests or from a more sophisticated model.

The second takeaway is that the reduced order model produces virtually the same results as the complete model.

Large DC Motor Simulation Results

The large motor C42-L90-W30 simulation results show that, when the electrical and mechanical time constants are comparable, the simplified model cannot capture the dynamic behavior of the motor, as it was the case with the small motor.

However, the steady-state values still correlate fairly well with the published table values. A rotor speed of 146 rad/s vs. 159 rad/s, and an armature current of 4.37 A vs. 5 A.

As alluded to in the previous case, even the complete model used in the simulation is still an approximation when it comes to accurately describing the behavior of the DC motor.

Simplified DC Motor Model

With the simulation results in mind, we can use the simplified (or reduced order) model to describe the dynamic response of the types of DC motors that I’ll be dealing with throughout my blog.

From the block diagram it’s possible to write:

And rearranging:

Doing the inverse Laplace Transform to put in back in differential equation form:

or

Where:

As a quick note, if we take a closer look at the reduced model, we can identify the mechanical time constant, in the particular case where there are no additional rotating inertias reduced to the rotor shaft.

At this point, the motor can be characterized if we determine the values for,and. There are three non-intrusive tests that can be done with an actual motor, from which the three aforementioned parameters can be extracted.

The first test is a torque stall test () which reduces our last differential equation to:

From whichcan be determined for a given input voltage and the corresponding load torque that causes the DC motor to stall.

The second test is a steady-state test () with no load torque, which gives us:

Sinceis known from the previous test and the angular velocity can be measured for a given input voltage, the parametercan therefore be determined.

The third test is a transient step-input test with no load, which reduces our simplified equation to:

By measuring the response time of the first-order system,can be finally determined (sinceis already known at this point).

Final Remarks

We have seen how it’s possible to determine, or better said estimate, the main parameters , andfor the reduced order model of a DC motor. This can be accomplished without having to disassemble the motor to measure the rotor moment of inertia, or having a more sophisticated electrical measurement setup to determine the armature inductance.

The three tests described above can be performed on a simple test bench. Part 2 of this post will go over an experimental setup with the LEGO EV3 large motor and some automated testing using the pyev3 python package. Stay tuned!

Event Detection in Signal Processing

One thing I always come across when analyzing time-based data is the need to identify events such as changes in signal level, triggers, among others. Many years ago I created a MATLAB function that can find clusters of data values inside a time series. It returns the indexes of the starting point of each cluster, with the corresponding number of samples within each cluster. I have been using the Python version of the function quite a bit recently. Apparently it never gets old!

By doing the proper manipulation of the data, virtually any type of event can be detected and extracted from the time series.

The example on the left shows a generic signal where we search for clusters of data (or events) with values greater than zero. In most cases, it’s good practice to apply a post-processing filter to the signal as the very first step of the data analysis.

The code for the python function used to find the clusters is shown below. Let’s get started by showing some basic usage and then some more interesting cases.

```def find_cluster(x, xval):
"""
Find clusters of data in an ndarray that satisfy a certain condition.

:param x: The array containing the data for the cluster search.
:type x: ndarray

:param xval: The value of x that has to be satisfied for clustering.
:type xval: integer, float

:returns: 2-tuple

* i0:
The index of each cluster starting point.

* clustersize:
The corresponding lengths of each cluster.

:rtype: (list, list)

Example
-------
>>> x = np.int32(np.round(np.random.rand(20)+0.1))
>>> i0, clustersize = find_cluster(x, 1)

"""
# Cluster information list
a = []
# Initial (place holder) values for cluster start and end points
kstart = -1
kend = -1
# Going through each value of x
for i, xi in enumerate(x):
if xi == xval:
# Assigning cluster starting point
if kstart == -1:
kstart = i
# Assigning cluster end point for particular case
# when there is an xval in the last position of x
if i == len(x)-1:
kend = i
else:
# Assigning cluster end point
if kstart != -1 and kend == -1:
kend = i-1
# Updating cluster information list
# and resetting kstart and kend
if kstart != -1 and kend != -1:
a.append(kstart)
a.append(kend)
kstart = -1
kend = -1
# Assigning indeces of cluster starting points
# (Every other list element starting from position 0)
i0 = a[0:-1:2]
# Assigning cluster sizes
# (Every other list element starting from position 1)
clustersize = list(np.array(a[1::2]) - np.array(i0) + 1)
# Case where cluster size is ZERO
if len(i0) == 0:
i0 = []
clustersize = []
return i0, clustersize
```

The example in the doc string above generates a random array `x` of zeros and ones, in which we use the function to find the clusters of values equal to 1. By running it, one should get something like the results shown below. `i0` contains the starting indexes of the clusters of 1s and `clustersize `contains the number of elements in each cluster.

```>>> x = np.int32(np.round(np.random.rand(20)+0.1))
>>> i0, clustersize = find_cluster(x, 1)

>>> x
array([1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0])

>>> i0
[0, 4, 6, 13, 17]

>>> clustersize
[3, 1, 6, 2, 2]
```

Finding Signal Peaks

Let’s use the signal from the top of the post and find the signal peaks above 1 and below -1, as illustrated on the left.

The code snippet used to find the peaks as well as the peak values can be seen below. The complete code for it and for all the examples in this post can be found here.

```# Finding positive peaks greater than 1
imax = []
icl, ncl = find_cluster(x>1, 1)
for ik, nk in zip(icl, ncl):
imax.append(ik+np.argmax(x[ik:ik+nk])-1)
# Finding negative peaks smaller than 1
imin = []
icl, ncl = find_cluster(x<-1, 1)
for ik, nk in zip(icl, ncl):
imin.append(ik+np.argmin(x[ik:ik+nk])-1)
```
```>>> np.round(x[imax], 2)
array([1.3 , 1.03, 1.47, 1.06, 1.18, 1.53, 1.21, 1.19])

>>> np.round(x[imin], 2)
array([-1.25, -1.19, -1.53, -1.24, -1.07, -1.47, -1.08, -1.3 ])
```

Finding Trigger Edges

`find_cluster` can be very useful for determining trigger events in a recorded time series. Unlike most other signals, I usually do not filter a trigger signal, since having a sharp rising edge (at the expense of noise) is critical for a more exact event determination.

The derivative of the trigger signal is then used to determine the rising edges, as shown in the second plot on the left. A threshold value of 0.5 for the derivative (normalized by its maximum value) can be used to locate the trigger edges.

The code snippet below shows how to do it for a fictitious recorded data set.

```# Calculating and normalizing trigger signal derivative
dxdt = dxdt/np.max(dxdt)

# Finding trigger event times using derivative theshold
icl, ncl = find_cluster(dxdt>0.5, 1)
ttrigger = t[icl]
```
```>>> ttrigger
array([0.0845, 0.1975, 0.3065, 0.402 , 0.485 , 0.559 , 0.6305, 0.7055, 0.796 , 0.9495])
```

DC Motor Position Tracking

At the end of the post Motor Position Control with Raspberry Pi, we saw that there are different ways to move a motor shaft between two angular positions. In this post, we will take a deeper look into three types of paths between an initial and a final position. While I’ll be exploring the concept using a system with one degree of freedom, with proper normalization and parametrization, the idea can be applied to any trajectory between two points in space.

The first observation to be made is that the path that our closed-loop system is going to follow is calculated before hand, i.e., the initial and final positionsand(shown in the diagram on the left), as well as other parameters that define the trajectory, are known before the position tracking is executed.

In our case, we’ll use the acceleration time(at the beginning and end of the motion) and the maximum allowed velocityas the additional parameters that are sufficient to fully determine the path.

Maximum velocity and acceleration time (or conversely the maximum acceleration) are a reasonable choice, since they are related to the physical limitations of the system.

The total motion time is therefore automatically determined once the four parameters above are specified. Even though there are different ways to go about calculating the path, I’ll use the velocity profile as the primary trace, with position and acceleration being calculated respectively as the integral and derivative of the velocity. Let’s then go over the three velocity profiles for this post: linear, quadratic, and trigonometric.

Linear Speed Profile

Using the illustration below it’s not too difficult to arrive at the following three equations that describe the velocity as a function of time.

The timeat maximum speed is given by:

Where it can be calculated by using the fact that the area under the speed trace is equal to .

As mentioned before, and can be used interchangeably, where, in this case, . The Python function below shows one possible implementation of the path calculation, where the position and acceleration traces are calculated through numerical integration and derivative, respectively. I found this approach to be more robust to the path discretization, specially at the speed profile transitions.

```def path_linear(xstart, xstop, vmax, ta, nstep=1000):
#
# Adjusting velocity based on end points
if xstop < xstart:
vmax = -vmax
# Assigning time at max velocity
if (xstop-xstart)/vmax > ta:
# There's enough time for constant velocity section
tmax = (xstop-xstart)/vmax - ta
else:
# There isn't (triangular velocity profile)
tmax = 0
vmax = (xstop-xstart)/ta
# Assigning important time stamps
t1 = ta  # End of acceleration section
t2 = ta + tmax  # End of constant velocity section
t3 = 2*ta + tmax  # End of motion (deceleration) section
# Creating time array for path discretization
t = np.linspace(0, t3, nstep+1)
# Finding transition indices
i1 = np.nonzero(t<=t1)
i2 = np.nonzero((t>t1)&(t<=t2))
i3 = np.nonzero(t>t2)
# Calculating accelereration section array
v1 = vmax/ta*t[i1]
# Calculating constant velocity section array
if len(i2[0]) > 0:
v2 = np.array([vmax]*len(i2[0]))
else:
v2 = np.array([])
# Calculating deceleration section array
v3 = vmax - vmax/ta*(t[i3]-t2)
# Concatenating arrays
v = np.concatenate((v1, v2, v3))
# Calculating numeric integral (position)
x = xstart + t[-1]/nstep * np.cumsum(v)
# Calculating numeric derivative (acceleration)
# Returning time, position, velocity, and acceleration arrays
return t, x, v, a
```

Unlike the linear speed profile, which uses constant acceleration, the quadratic (or parabolic) speed profile uses linearly changing acceleration at the beginning and end of the motion. The equations for this type of trajectory are shown below.

An easy way to obtain the first and third velocity equations above is by writing the corresponding acceleration functions and integrating them.

Similarly to the linear case, the time at maximum speed can be obtained from the area under the velocity curve. (Tip: integrate the first velocity equation and evaluate it at, multiply the result by 2 to account for the deceleration portion, and add the area of the rectangle under)

The Python function for this case follows the same structure as the previous one. You can find it on my GitHub page.

Trigonometric Speed Profile

The third profile type uses a cosine function for the velocity during the acceleration and deceleration phases of the trajectory. By choosing the appropriate function, it is possible to add an interesting feature to this last profile type: zero initial and final accelerations.

Also through integration of the velocity function and knowing that the area under the curve is, it is possible to solve for the time at maximum velocity:

Like the two previous cases, the Python function that implements this type of path calculation is located on the GitHub page.

Practical Application for Position Tracking

When it comes to choosing the path calculation method, two factors have to be kept in mind: how fast do we want to (or can) go from the initial to the final position and how much acceleration change (jerk) do we want to apply to the system. One good way to go about it is to test all three and quantify how the overall closed-loop system behaves.

Just to explore a little bit, the plot on the left shows an overlay of the three types of path for the examples shown above. In this scenario, the maximum velocity is limited to 1 (m/s, for example) and the total time set to 1.4s, by proper choice of. Even though the acceleration profiles are very distinct, the difference between the three types of trajectory is visually very subtle.

In contrast, if we now limit both the maximum velocity to 1 m/s and the maximum acceleration to 4 m/s2, the total motion duration is fairly different. The linear velocity profile is the fastest one (at the expense of having the largest jerk). Conversely, the trigonometric method takes the longest time with the least amount of jerk.

Let’s use the same setup with a Raspberry Pi and a DC motor that we had in the post Motor Position Control with Raspberry Pi to test the different types of paths for the positioning of the DC motor shaft. The Python program below shows how the quadratic path is generated prior to the execution loop and then sampled during the execution loop.

The full code with the plotting features and the module `path.py`, containing the three types of path calculation functions, can be found here.

```# Importing modules and classes
import time
import numpy as np
from scipy.interpolate import interp1d
from gpiozero_extended import Motor, PID

# Setting path parameters
thetastart = 0  # Start shaft angle (rad)
thetaend = 3*np.pi  # End shaft angle (rad)
wmax = 4*np.pi  # Max angular speed (rad/s)
ta = 0.45  # Acceleration time (s)
tsample = 0.01  # Sampling period (s)

# Calculating path
t, x, _, _ = path_quad(thetastart, thetaend, wmax, ta)
f = interp1d(t, 180/np.pi*x)  # Interpolation function
tstop = t[-1]  # Execution duration (s)

# Creating PID controller object
kp = 0.036
ki = 0.260
kd = 0.0011
taupid = 0.01
pid = PID(tsample, kp, ki, kd, tau=taupid)

# Creating motor object using GPIO pins 16, 17, and 18
# (using SN754410 quadruple half-H driver chip)
# Integrated encoder on GPIO pins 24 and 25.
mymotor = Motor(
enable1=16, pwm1=17, pwm2=18,
encoder1=24, encoder2=25, encoderppr=300.8)
mymotor.reset_angle()

# Initializing variables and starting clock
thetaprev = 0
tprev = 0
tcurr = 0
tstart = time.perf_counter()

# Running execution loop
print('Running code for', tstop, 'seconds ...')
while tcurr <= tstop:
# Pausing for `tsample` to give CPU time to process encoder signal
time.sleep(tsample)
# Getting current time (s)
tcurr = time.perf_counter() - tstart
# Getting motor shaft angular position
thetacurr = mymotor.get_angle()
# Interpolating set point angle at current time step
if tcurr <= tstop:
thetaspcurr = f(tcurr)
# Calculating closed-loop output
ucurr = pid.control(thetaspcurr, thetacurr)
# Assigning motor output
mymotor.set_output(ucurr)
# Updating previous values
thetaprev = thetacurr
tprev = tcurr

print('Done.')
# Stopping motor and releasing GPIO pins
mymotor.set_output(0, brake=True)
del mymotor
```

The plots show the set point vs. actual shaft position for a 450 degrees rotation, using the quadratic velocity path.

The derivative of the signal, i.e., the shaft speed, is a useful gage to see how well the closed-loop system can follow the trajectory. The position error between set point and measured value is also a good metric. In this case, it is within +/- 5 degrees with an RMSE (Root Mean Squared Error) of about 2.5 degrees. Although not a stellar result, it’s not too bad, given the budget hardware and simple control (PID) that were used.

Finally, the graph at the bottom shows a comparison of the position tracking error for the three types of path explored in this post for the same 450 degree rotation. The average RMSE values for the linear, quadratic, and trigonometric speed profiles were respectively 3.3, 2.4, and 3.6 degrees.

Even though the trigonometric profile has the smoothest acceleration, it was outperformed by the other two for this particular system.

The trigonometric one is usually beneficial for dynamic systems that are more susceptible to vibration when subjected to a more “jerky” path profile.

Motor Position Control with Raspberry Pi

Along the lines of the Motor Speed Control post, let’s reuse some of our Python classes to control the angular position of a DC motor. We’ll also look into how to tune the PID using the Ziegler-Nichols method, as well as different ways to apply a position set point input.

Like before, the starting point is the Raspberry Pi setup that can handle both a PWM output and an encoder input. Unlike the speed control scenario, there’s no need to calculate and filter the speed from the measured angular position.

Most of the work presented in this post was done using the Python code below. It uses the Motor class and the Digital PID class as its building blocks. On my GitHub page, you can find the complete version of this program, a slight variation of it (used for the PID tuning), and all the classes that it uses.

```# Importing modules and classes
import time
import numpy as np
from gpiozero_extended import Motor, PID

# Setting general parameters
tstop = 1  # Execution duration (s)
tsample = 0.01  # Sampling period (s)
thetamax = 180  # Motor position amplitude (deg)

# Setting motion parameters
# (Valid options: 'sin', 'cos')
option = 'cos'
if option == 'sin':
T = 2*tstop  # Period of sine wave (s)
theta0 = thetamax  # Reference angle
elif option == 'cos':
T = tstop  # Period of cosine wave (s)
theta0 = 0.5*thetamax  # Reference angle

# Creating PID controller object
kp = 0.036
ki = 0.379
kd = 0.0009
taupid = 0.01
pid = PID(tsample, kp, ki, kd, tau=taupid)

# Creating motor object using GPIO pins 16, 17, and 18
# (using SN754410 quadruple half-H driver chip)
# Integrated encoder on GPIO pins 24 and 25.
mymotor = Motor(
enable1=16, pwm1=17, pwm2=18,
encoder1=24, encoder2=25, encoderppr=300.8)
mymotor.reset_angle()

# Initializing variables and starting clock
thetaprev = 0
tprev = 0
tcurr = 0
tstart = time.perf_counter()

# Running execution loop
print('Running code for', tstop, 'seconds ...')
while tcurr <= tstop:
# Pausing for `tsample` to give CPU time to process encoder signal
time.sleep(tsample)
# Getting current time (s)
tcurr = time.perf_counter() - tstart
# Getting motor shaft angular position
thetacurr = mymotor.get_angle()
# Calculating current set point angle
if option == 'sin':
thetaspcurr = theta0 * np.sin((2*np.pi/T) * tcurr)
elif option == 'cos':
thetaspcurr = theta0 * (1-np.cos((2*np.pi/T) * tcurr))
# Calculating closed-loop output
ucurr = pid.control(thetaspcurr, thetacurr)
# Assigning motor output
mymotor.set_output(ucurr)
# Updating previous values
thetaprev = thetacurr
tprev = tcurr

print('Done.')
# Stopping motor and releasing GPIO pins
mymotor.set_output(0, brake=True)
del mymotor
```

PID Tuning using Ziegler-Nichols

Unfortunately, the PID tuning method presented in the previous post cannot be used to obtain an initial set of controller gains. The first-order approximation used for the speed response doesn’t work for the angular position, since it increases linearly once a steady-state speed is reached. In other words, the DC motor system with position as an output is an integrator of the corresponding system with speed as an output.

That brings us to the Ziegler-Nichols method, which can be used with a closed-loop system. In our case, the feedback loop is closed on the motor shaft position. The method is based on starting the tuning process only with a proportional gain and then increasing its value, until the closed-loop system becomes marginally stable. I.e., oscillates with a constant amplitude. The four plots below show the response and corresponding controller output for increasing values of.

The first question that comes to mind is what should the starting value ofbe. A good guess is based on the size of the step input and the saturation limit of the controller. The very instant the step input is applied, the position error going into the controller is equal to the step size. Therefore, at that moment, the controller output is given by . As shown in the first plot, for , the initial controller output is 0.6 for a step input of 20 degrees. The choice of the step size and the initial proportional gain are such that the controller doesn’t saturate (in our case, 1 is the output limit) and there’s still room to increase the gain without running “too much” into the saturation zone.

The ultimate proportional gainis achieved when the system oscillates with a constant amplitude. In our example, . Under this condition, the ultimate oscillation period can be extracted from the last plot above and is calculated as.

The Ziegler-Nichols method states that for a 1/4 decay ratio of the closed-loop system response, the PID gains are given by:

The graph shows the closed-loop response for a shaft position step input of 60 degrees. The PID gains were obtained after substitutingandinto the formulas above.

Observe the 1/4 decay ratio of the amplitude. Also, the fact that the motor position cannot meet exactly 60 degrees causes the integral error to build up, resulting in small position “jumps” around the 1.4 and 1.7 s marks.

Starting from the gains defined by using the Ziegler-Nichols method, let’s do some manual tuning aiming to decrease the oscillation of the closed-loop response. That can be achieved by decreasing the integral gain and increasing the derivative gain.

The plot on the left shows the system response for values that were obtained after a few attempts. The idea is to later compare the effect of the two sets of gains on the response (when applying different types of position set point inputs).

Sine and Cosine Path Functions

While using step changes to move from one position to another is fine, there are situations where we want the position to follow a pre-determined path. By choosing the right combination of PID gains and characteristics of the path function, it’s possible to greatly improve the tracking ability of the position controller.

In this post let’s take a look at two trigonometric functions (sine and cosine), that can be used to move the motor to a given set point angle and then back to its initial position. By inspecting the code above, we can identify the sine and cosine path functions shown below.

and

By choosing suitable values of the amplitudeand period, the two functions can produce motion paths that have the same set point (maximum) angle and the same motion duration to return to the original angular position.

The next two plots show the general closed-loop tracking performance for the sine-based path function, using the two sets of PID gains obtained in the previous section. Notice the more pronounced oscillatory behavior for the more aggressive set of gains. Additionally, due to the nature of the sine-based path, the initial and final speeds are not zero.

The cosine-based path function on the other hand (shown below) has zero initial and final speeds. That characteristic can greatly reduce the system oscillation, especially at the very beginning of the motion.

Final Remarks

When it comes to position control, the DC motor is usually part of a machine or even a robot, where it is important to understand the dynamic behavior of the complete system. Speed and accuracy are usually conflicting goals that have to be compromised to obtain the best performance of the closed-loop control system. However, having a suitable path function can minimize undesired system oscillation caused by the input excitation (or path).

In future posts we will take a deeper look into path functions, particularly the ones that in addition to having zero end-point speeds also have zero end-point accelerations. Stay tuned!

PID Tuning

Tuning a PID controller consists of choosing the best values for the proportional, integral, and derivative gains in order to satisfy some performance criteria of the closed-loop system. If we’re able to characterize or model the plant that we want to control, then it’s possible do find the best gains either through direct calculation or through simulation.

If we go back to the day when computers weren’t exactly around and there was no such thing as a digital PID, folks had to come up with methods that could be applied directly to the plant or process. Developed respectively in the 1940s and 50s, the well known Ziegler-Nichols and Cohen-Coon methods were the way to go about PID tuning back then. One problem with those two methods is that they result in a fairly aggressive set of gains that can potentially lead to an unstable (or a reduced stability margin) closed-loop system. That undesired characteristic becomes even more critical when a digital PID is employed.

Because we will be dealing mostly with DC motors, I’ll go over a method where the motor speed response is approximated by a first order system. The identified open-loop system is then used to determine the PID gains for the corresponding closed-loop system. In a future post we will go over characterizing the motor individual plant parameters, thus obtaining a better model.

So, let’s try to keep the math under control and start with characterization of the first-order system.

First-Order System

If we look back at the digital filter post, we can extend the low-pass filter equation to a more general form of the first-order system, given by:

where the steady-state gainis the ratio between system output and input. As seen in that same post,is the response time of the system. By determining those two parameters, the first-order system can be fully characterized.

Closed-Loop System

The block diagram below represents the closed-loop system with a PID controller and our first-order system, the latter now expressed as a differential equation. While technically correct, no one uses a time-based representation of the system as shown in the diagram. Mostly because it is really difficult to work with differential equations whose inputs and outputs are all interconnected in some way.

Thankfully, a pretty smart guy named Pierre-Simon Laplace (1749 – 1827) came up with the Laplace transform. It allows for mere mortals to treat complicated convolution problems involving differential equations as simple algebraic ones. We can then redraw the diagram using the corresponding Laplace transforms, as shown in the more familiar form below.

Because we will be doing some algebraic manipulation let’s use the following transfer function notation:

Where

and

Doing the usual “with a few steps, it can be shown that” the closed-loop system transfer function in our case is:

It has the characteristic equation of a second-order system, with damping ratio(zeta) and natural frequency

or

Where

and

After all this mumbo jumbo, we’re finally getting to a point where, for a given first-order plant (with knownand), we can determine the PID gains,andthat would give us the desired closed-loop response, characterized byand. In other words, by selecting a damping ratio and a natural frequency, we can determine a set of PID gains that make our closed-loop system behave as desired.

PID Gains Selection

If we take a look at the last two equations, there are three unknowns (the PID gains) for two known parameters (and ). A good starting point is to make . That way, we can solve the two following equations, for and.

and

Where:

• in rad/s is the natural frequency of the closed-loop system. The larger its value the faster the system’s response time. In practical terms, you cannot make this number arbitrarily high since you may run into saturation of the controller output. While, for most applications, short saturation periods are fine, it’s better to keep them at a minimum.
• is the damping ratio of the closed-loop system. Values close to 1 (critical damping ratio) will keep the oscillations and overshoot to a step input under control. In some closed-loop positioning applications, overshoot isn’t an option. Hence, a damping ratio greater than 1 is required.

Finally, note that the PID gains should always be positive. Therefore,.

DC Motor PID Tuning

It’s time to use the hardware setup and code from the Motor Speed Control post to do some exploration of this PID tuning method. First, a step input is applied to the open-loop motor in order to determine the first order response parametersand.

For our little DC motor, with the Raspberry Pi setup, the steady state gain is 20.9/0.5 = 41.8 rad/s/u. Remember that, in our case, u is the normalized PWM output between -1 and 1. The response time of the open-loop system is 0.184 s.

To improve the parameter extraction, a zero-phase filter was used. The code that automatically determines the parameters from the system response can be found on my GitHub page.

Next step consists of choosing values for the closed-loop system’s natural frequencyand damping ratio. The graphs on the left show the system response to a step input for a natural frequency of 2 Hz () and. Using the last two formulas for the integral and proportional gains

,

with all the values we gathered so far, would result inand. We can observe that the motor now takes 0.11 s to reach the target speed (instead of 0.184 to reach 63% of it). That comes of course at expense of using the full capacity of the PWM with some brief saturation.

The reasoning behind my choice ofandis that I wanted a faster response with limited saturation. 2 Hz was the right amount for that. I also wanted a damping ratio slightly greater than one to avoid any oscillations. At this point you might be asking yourself: why is the closed-loop response showing an oscillatory behavior then! The answer is a combination of:

• The motor is in fact a second order system, not a first order one as our approximation suggests. So, higher order modes may start showing up in the closed-loop response
• The PID is discrete, not continuous (as that was the basis for all the deductions above). And we’ve seen that the digitalization of the PID tends to promote a more oscillatory behavior
• The actual closed-loop system includes the transfer function of a low-pass filter for the motor speed. So, the overall system is not composed only by a PID and a plant, as we initially assumed.

Acknowledging that the choice of gains was fairly straightforward, the resulting closed-loop response is quite good! And we haven’t incorporated the derivative gain yet. While it’s usage is not always necessary, it tends to be quite helpful in the case of closed-loop DC motor control. A good starting point foris a couple of orders of magnitude smaller than the proportional gain. So, after starting with 0.001, I ended up with a value of 0.005. The addition of the derivative term really decreased the overshoot and virtually removed the oscillations.

Finally, the next set of graphs shows the results when choosing different natural frequencies and damping ratios. Notice the effect on the response time, the overshoot and the saturation of the controller output.

Motor Speed Control with Raspberry Pi

Let’s use the Motor class with encoder and the Digital PID class to create a closed-loop speed controller for a DC motor. Those two classes have all the required attributes and methods that we need to program an execution loop for a PID controller with very few lines of code.

The starting point is the Raspberry Pi setup that was used to implement encoder capability into the Motor class. I took the artistic liberty to create the diagram below, where I combine visual elements of the PID loop, the hardware connections, and the class methods used for the I/O between the Pi computer and the physical system.

Going back to the concept of execution loops, the I/O tasks are taken care by the Motor class methods `get_angle()` and `set_output()`. The computation tasks happening in between are essentially the shaft speed calculation, low-pass filtering, and PID control action. The Python code below implements the execution loop, where the different elements in the Pi box above can be identified. The complete program, which includes the plotting features, along with the Python module with all necessary classes, can be found on my GitHub page.

```# Importing modules and classes
import time
import numpy as np
from gpiozero_extended import Motor, PID

# Setting general parameters
tstop = 2  # Execution duration (s)
tsample = 0.01  # Sampling period (s)
wsp = 20  # Motor speed set point (rad/s)
tau = 0.1  # Speed low-pass filter response time (s)

# Creating PID controller object
kp = 0.15
ki = 0.35
kd = 0.01
taupid=0.01
pid = PID(tsample, kp, ki, kd, umin=0, tau=taupid)

# Creating motor object using GPIO pins 16, 17, and 18
# (using SN754410 quadruple half-H driver chip)
# Integrated encoder on GPIO pins 24 and 25.
mymotor = Motor(
enable1=16, pwm1=17, pwm2=18,
encoder1=24, encoder2=25, encoderppr=300.8)
mymotor.reset_angle()

# Initializing previous and current values
ucurr = 0  # x[n] (step input)
wfprev = 0  # y[n-1]
wfcurr = 0  # y[n]

# Initializing variables and starting clock
thetaprev = 0
tprev = 0
tcurr = 0
tstart = time.perf_counter()

# Running execution loop
print('Running code for', tstop, 'seconds ...')
while tcurr <= tstop:
# Pausing for `tsample` to give CPU time to process encoder signal
time.sleep(tsample)
# Getting current time (s)
tcurr = time.perf_counter() - tstart
# Getting motor shaft angular position: I/O (data in)
thetacurr = mymotor.get_angle()
wcurr = np.pi/180 * (thetacurr-thetaprev)/(tcurr-tprev)
# Filtering motor speed signal
wfcurr = tau/(tau+tsample)*wfprev + tsample/(tau+tsample)*wcurr
wfprev = wfcurr
# Calculating closed-loop output
ucurr = pid.control(wsp, wfcurr)
# Assigning motor output: I/O (data out)
mymotor.set_output(ucurr)
# Updating previous values
thetaprev = thetacurr
tprev = tcurr

print('Done.')
# Stopping motor and releasing GPIO pins
mymotor.set_output(0, brake=True)
del mymotor
```

Motor Speed Filtering

Now that the Python code is laid out, let’s go over some of its components. The motor shaft speed is calculated as the derivative of the angular position of the shaft using a backward difference formula `wcurr = np.pi/180 * (thetacurr-thetaprev)/(tcurr-tprev)`.

The encoder that’s being used has a resolution of `encoderppr=300.8`, which corresponds to approximately 360/300.8 = 1.2 degrees. The difference between two consecutive time steps `tcurr` and `tprev `is about `tsample = 0.01 `seconds. That gives us a not so great speed resolution of `np.pi/180 * (1.2)/(0.01) = 2` rad/s (or 19 rpm!). As a side note, the speed resolution could be improved by increasing either the encoder PPR or the sampling rate. However the first option might reach the limit of the GPIO zero RotaryEncoder class, while the latter will eventually make the closed-loop system unstable. Therefore a compromise has to be reached given the constraints at hand.

One way reduce the impact of this limitation is to apply a low-pass digital filter to the motor speed signal: `wfcurr = tau/(tau+tsample)*wfprev + tsample/(tau+tsample)*wcurr`. The graphs below show the raw and the filtered signal as well as the controller output.

On the left, `tau = 0` and therefore no filtering is being done. Observe how the controller output bounces between saturation limits. On the right, `tau = 0.1`, causing the control output to be less erratic, which in turn makes the raw speed more stable. One can also see the “discretization effect” of the 2 rad/s speed resolution on the top row plots.

As discussed in the post Encoder with Raspberry Pi, we cannot use a more accurate sampling period during the execution loop due to the intensive computing resources of the RotaryEncoder class.

However, as shown on the left, having to rely on the command `time.sleep(tsample)` inside the execution loop, still produces repeatability of the sampling period within 1 ms.

PID Derivative Term Filtering

If you notice in the Python code above, I am also filtering the PID derivative when I created the PID class instance: `pid = PID(tsample, kp, ki, kd, umin=0, tau=taupid)`, where `taupid = 0.01` seconds.

If we turn it off by making `taupid = 0`, we obtain the results on the left, where the controller output shows a bigger oscillation amplitude. However, the effect of turning the derivative term filter off is not as pronounced as the effect of not filtering the motor speed signal.

By comparing all the system responses so far, even though the PID gains were never changed from case to case, the closed-loop system response is changing. Ultimately, by adding low-pass filters, the closed-loop dynamic system is different from its original form without the filters.

PID Saturation and Anti-Windup

As it can be seen in the previous graphs, the PID output is saturating right after the step input. We can observe the effect of PID saturation and the integral term windup (which has been avoided so far) by creating the PID class instance: `pid = PID(tsample, kp, ki, kd, umin=-10, umax=10, tau=taupid)`.

Because the controller output (PWM in the Motor class) is in fact limited between -1 and 1, applying -10 and 10 as the new PID controller limits will effectively turn off the anti-windup feature and cause the integral term to wind up.

The plots on the left show the resulting speed overshoot due to that modification as the actual PWM output is saturated at 1.

Final Remarks

We have been building up our knowledge with the previous posts (namely: Encoder with Raspberry Pi, Digital Filtering, and Digital PID Controller) so we could finally run a DC motor in closed-loop speed control. While this could have been achieved from the very beginning, I believe it was more educational to do it in smaller parts through the use of Python classes.

When building more complex automation systems, having a modularized approach with classes that can be tested and used independently, as well as combined in different ways, will make things much easier to tackle in the long run.