This notebook contains material from cbe30338-2021; content is available on Github.
Setpoints are functions of time that establish target values for key control variables. This notebook describes typical nomenclature used in describing setpoint functions, and shows how to creatd setpoint functions in Python.
Example descriptions from commercial vendors:
Common descriptions for setpoint functions include so-called step changes, and ramp and soak periods.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
# specify a setpoint profile
profile = [
(0, 25),
(50, 35),
(120, 35),
(120, 40),
(170, 50),
(170, 50),
(270, 50),
(370, 40),
(450, 40),
(600, 25),
]
# convert to numpy array to provide simple access to columns
profile = np.array(profile)
# create an annotated plot
fig, ax = plt.subplots(1, 1, figsize=(10, 5))
ax.plot(profile[:, 0], profile[:, 1], lw=4)
ax.annotate("Step", xy=(120, 37.5), xytext=(180, 37.5), fontsize=20,
va="center", arrowprops=dict(facecolor='black', shrink=0.05))
ax.annotate("Ramp", xy=(320, 45), xytext=(380, 48), fontsize=20,
va="center", arrowprops=dict(facecolor='black', shrink=0.05))
ax.annotate("Soak/Dwell", xy=(210, 50), xytext=(210, 55), fontsize=20,
ha="center", arrowprops=dict(facecolor='black', shrink=0.05))
ax.set_ylim(25, 60)
ax.set_title("Sample Setpoint Profile")
ax.set_xlabel("Time")
ax.grid(True)
Study Question: Classify all of the segments in the sample setpoint profile.
Study Question: What is the ramp rate of the first ramp in the example above.
Study Question: Modify the data in the above example to remove the step. Replace it with a single raamp from the initial condition to the soak period that begins at t=170 at a temperature of 170C.
For feedback control we would like functions that return the value of a setpoint for any point in time. Functions are in the form $SP_1(t)$ and $SP_2(t)$, for example, are straightfoward to use inside in control applications.
In the section we show how to write a function that accepts points defining a piecewise linear setpoint profile, then produce a function to compute the setpoint for any point in time.
Describiing the setpoint as a series of a step/ramp/soak periods naturally leads a piecewise linear function. The start and end of each line segment are spceified by (time, value) pairs. Ordering these points into a list provides a straightforward specification of the setpoint,
Here we show the points for a typical setpoint.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
profile = [
(0, 25),
(50, 35),
(120, 35),
(120, 40),
(170, 50),
(170, 50),
(270, 50),
(370, 40),
(450, 40),
(600, 25),
]
t = [t for t, _ in profile]
y = [y for _, y in profile]
fig, ax = plt.subplots(1, 1)
ax.plot(t, y, '.', ms=10)
[<matplotlib.lines.Line2D at 0x7fe217098fd0>]
Python functions are frequently written to accept data, perform calculations, and return values. What may be less familiar is that functions can also return function. This is just what we neeed - a function that accepts a series of (time, value) pairs describing a setpoint profile, then returns a function that can be used to find values of the setpoint at any point in time.
def create_setpoint_function(profile):
profile = np.array(profile)
ti = profile[:, 0]
yi = profile[:, 1]
# define a function to interpolate time and values
def setpoint_function(t):
return np.interp(t, ti, yi)
# return that function
return setpoint_function
The following example creates a setpoint function that produces setpoints corresponding the profile described in the introduction to this notebook.
sp = create_setpoint_function(profile)
# print select values
t = 100
print(f"at time = {t:3d} setpoint = {sp(t)}")
at time = 100 setpoint = 35.0
We can use the setpoint function to create plots.
# compute setpoint values
t = np.linspace(0, 600, 600)
y = sp(t)
# create a plot
fix, ax = plt.subplots(1, 1, figsize=(10, 5))
ax.plot(t, y)
ax.set_xlabel("time / seconds")
ax.set_title("setpoint function")
ax.grid(True)
The next cell demonstrates the use of create_setpoint_function to create multiple independent setpoint functions.
T_amb = 21.0
sp_1 = create_setpoint_function([[0, T_amb], [20, T_amb], [60, 50], [100, 50], [140, T_amb]])
sp_2 = create_setpoint_function([[0, T_amb], [0, 45], [120, 35], [200, T_amb]])
# create plot axes
fig, ax = plt.subplots(2, 1)
# plot setpoint functions
t = np.linspace(-1, 250, 250)
ax[0].plot(t, sp_1(t))
ax[1].plot(t, sp_2(t))
[<matplotlib.lines.Line2D at 0x7fe21721f1c0>]
The goal of this next section is to create a function that returns the value of a setpoint profile for a PCR thermal cycler. We'll break this into a series of steps:
Here's an example of a PCR thermal cycler protocol:
The details of these protocols vary depending on the nature of the test, the reagents used, and the type of detection that will be employed. In real-time PCR, the number of cycling steps will end once a positive result is obtained.
The PCR protocol is a series of (time, temperature) pairs. The code in the next cell represents a protocol as a sequence of (time, temperature) pairs in a Python list. The list is constructed by concatenating subprotocols denoting the activation, cycling and extension steps in a PCR protocol
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
# number of cycles
n_cycles = 5
# enter steps in the protocol as a list of (time, temperature) pairs
activation = [(900, 95)]
cycling = [(20, 94), (20, 60), (30, 72)]*n_cycles
extension = [(600, 72)]
finish = [[0, 30]]
# concatenate into a list of (time, temperature) intervals
protocol = np.concatenate([activation, cycling, extension, finish])
protocol
array([[900, 95],
[ 20, 94],
[ 20, 60],
[ 30, 72],
[ 20, 94],
[ 20, 60],
[ 30, 72],
[ 20, 94],
[ 20, 60],
[ 30, 72],
[ 20, 94],
[ 20, 60],
[ 30, 72],
[ 20, 94],
[ 20, 60],
[ 30, 72],
[600, 72],
[ 0, 30]])
Each step in a PCR protocol consists of an initial ramp to the specified temperature, followed by a soak for the specified time and temperature. We begin the coding by demonstrating how to interpret the protocol specification as a series of ramp and soak periods.
# each soak period is preceeded by a ramp
for time, temp in protocol:
print(f"Ramp to {temp}C")
print(f"Soak at {temp}C for {time:3d} seconds." )
Ramp to 95C Soak at 95C for 900 seconds. Ramp to 94C Soak at 94C for 20 seconds. Ramp to 60C Soak at 60C for 20 seconds. Ramp to 72C Soak at 72C for 30 seconds. Ramp to 94C Soak at 94C for 20 seconds. Ramp to 60C Soak at 60C for 20 seconds. Ramp to 72C Soak at 72C for 30 seconds. Ramp to 94C Soak at 94C for 20 seconds. Ramp to 60C Soak at 60C for 20 seconds. Ramp to 72C Soak at 72C for 30 seconds. Ramp to 94C Soak at 94C for 20 seconds. Ramp to 60C Soak at 60C for 20 seconds. Ramp to 72C Soak at 72C for 30 seconds. Ramp to 94C Soak at 94C for 20 seconds. Ramp to 60C Soak at 60C for 20 seconds. Ramp to 72C Soak at 72C for 30 seconds. Ramp to 72C Soak at 72C for 600 seconds. Ramp to 30C Soak at 30C for 0 seconds.
The next step is to introduce a ramp rate and add variables to track the start time and temperature for each segment. We will assume all ramp rates have the same absolute value.
# add varibles to track current time and temperature
# ramp period is determined by a "ramp_rate"
ramp_rate = 2.5
time_now = 0.0
temp_now = 21.0
for time, temp in protocol:
print(f"Time = {time_now:6.1f}: Ramp to {temp}C")
time_now += np.abs((temp - temp_now)/ramp_rate)
temp_now = temp
print(f"Time - {time_now:6.1f}: Soak at {temp}C for {time} seconds")
time_now += time
Time = 0.0: Ramp to 95C Time - 29.6: Soak at 95C for 900 seconds Time = 929.6: Ramp to 94C Time - 930.0: Soak at 94C for 20 seconds Time = 950.0: Ramp to 60C Time - 963.6: Soak at 60C for 20 seconds Time = 983.6: Ramp to 72C Time - 988.4: Soak at 72C for 30 seconds Time = 1018.4: Ramp to 94C Time - 1027.2: Soak at 94C for 20 seconds Time = 1047.2: Ramp to 60C Time - 1060.8: Soak at 60C for 20 seconds Time = 1080.8: Ramp to 72C Time - 1085.6: Soak at 72C for 30 seconds Time = 1115.6: Ramp to 94C Time - 1124.4: Soak at 94C for 20 seconds Time = 1144.4: Ramp to 60C Time - 1158.0: Soak at 60C for 20 seconds Time = 1178.0: Ramp to 72C Time - 1182.8: Soak at 72C for 30 seconds Time = 1212.8: Ramp to 94C Time - 1221.6: Soak at 94C for 20 seconds Time = 1241.6: Ramp to 60C Time - 1255.2: Soak at 60C for 20 seconds Time = 1275.2: Ramp to 72C Time - 1280.0: Soak at 72C for 30 seconds Time = 1310.0: Ramp to 94C Time - 1318.8: Soak at 94C for 20 seconds Time = 1338.8: Ramp to 60C Time - 1352.4: Soak at 60C for 20 seconds Time = 1372.4: Ramp to 72C Time - 1377.2: Soak at 72C for 30 seconds Time = 1407.2: Ramp to 72C Time - 1407.2: Soak at 72C for 600 seconds Time = 2007.2: Ramp to 30C Time - 2024.0: Soak at 30C for 0 seconds
Finally, we create a two column numpy array representing the setpoint profile. The first row in the array is the starting time and temperature. Each subsequent row constains the ending time and temperature of a segment.
# store the data in a list of time, temperature pairs marking the end of each period
ramp_rate = 0.5 # deg/sec
time_now = 0.0
temp_now = 21.0
# intialze a list with the starting time and temperature
SP_list = [[time_now, temp_now]]
# append the ending time and temperature to the list
for time, temp in protocol:
# ramp
time_now += np.abs((temp - temp_now)/ramp_rate)
temp_now = temp
SP_list.append([time_now, temp_now])
# soak
time_now += time
SP_list.append([time_now, temp_now])
# convert list to numpy array to access columns
SP_array = np.array(SP_list)
SP_array
array([[ 0., 21.],
[ 148., 95.],
[1048., 95.],
[1050., 94.],
[1070., 94.],
[1138., 60.],
[1158., 60.],
[1182., 72.],
[1212., 72.],
[1256., 94.],
[1276., 94.],
[1344., 60.],
[1364., 60.],
[1388., 72.],
[1418., 72.],
[1462., 94.],
[1482., 94.],
[1550., 60.],
[1570., 60.],
[1594., 72.],
[1624., 72.],
[1668., 94.],
[1688., 94.],
[1756., 60.],
[1776., 60.],
[1800., 72.],
[1830., 72.],
[1874., 94.],
[1894., 94.],
[1962., 60.],
[1982., 60.],
[2006., 72.],
[2036., 72.],
[2036., 72.],
[2636., 72.],
[2720., 30.],
[2720., 30.]])
fig, ax = plt.subplots(1, 1, figsize=(12, 5))
ax.plot(SP_array[:,0], SP_array[:,1], '.')
[<matplotlib.lines.Line2D at 0x7fe2158d2160>]
To implement control algorithms we will need to find values of the setpoint function at arbitrary points in time, not just at the start and finish of ramp or soak periods. The standard numpy library includes a linear interpolation function numpy.interp well suited to this purpose.
# create function to interpolate
def SP(t):
return np.interp(t, SP_array[:,0], SP_array[:, 1])
# plotting the setpoint function
t = np.linspace(-100, max(SP_array[:,0]), 2000)
fig, ax = plt.subplots(1, 1, figsize=(12,5))
ax.plot(t, SP(t))
ax.set_title("Setpoint function")
Text(0.5, 1.0, 'Setpoint function')
The last step is to put all of these steps together to create a function that generates a setpoint function. This is an example of a very powerful coding technique of nested functions.
# create the setpoint function from a specified PCR
def PCR_setpoint(protocol):
ramp_rate = 0.5 # deg/sec
time_now = 0.0
temp_now = 21.0
# intialze a list with the starting time and temperature
SP_list = [[time_now, temp_now]]
# append the ending time and temperature to the list
for time, temp in protocol:
# ramp
time_now += np.abs((temp - temp_now)/ramp_rate)
temp_now = temp
SP_list.append([time_now, temp_now])
# soak
time_now += time
SP_list.append([time_now, temp_now])
# convert list to numpy array to access columns
SP_array = np.array(SP_list)
def SP(t):
return np.interp(t, SP_array[:,0], SP_array[:, 1])
return SP
# create a setpoint function
setpoint = PCR_setpoint(protocol)
# plotting the setpoint function
t = np.linspace(0, 3000, 3000)
fig, ax = plt.subplots(1, 1, figsize=(12,5))
ax.plot(t, setpoint(t))
ax.set_title("Setpoint function")
Text(0.5, 1.0, 'Setpoint function')
Study Question: Change the protocol to include 30 thermal cycles, then create the setpoint function with PCR_setpoint() and plot the results.
Study Question: To better reflect the unequal heating and cooling rates available in most PCR devices, modify PCR_setpoint() to provide differing ramp rates for positive going and negative going ramps. Demonstrate the result using a postive ramp_rate of 2.5 degC/sec and a negative ramp_rate of -0.5 degC/sec.