Lab 4: PI Control#

This lab assignment introduces the use of propotional-integral (PI) control for the temperature control laboratory. In this assignment, you will:

  • testing a simple implementation of PI control

  • identifying issues with the simple implementation

  • test a more sophisticate PI control implementation

  • analyze the results and answer discussion questions

Implementation of a Simple PI Controller#

Given a process variable \(PV\) and setpoint \(SP\), proportional-integral-derivative control determines the value of a manipulated variable MV by the equation

(2)#\[\begin{align} MV(t) & = \bar{MV} + K_P\left(SP(t) - PV(t)\right) + K_I \int_0^t \left(SP(t')-PV(t')\right)dt' \end{align}\]

where \(K_P\) and \(K_I\) are the proportional and integral coefficients, respectively. The value \(\bar{MV}\) is a nominal or initial value of the manipulated variable.

The actual implementation of PI control is normally done by computing how much the \(MV\) should change at each time step. Defining the error at time \(k\) as

(3)#\[\begin{align} e_k & = SP_k - PV_k \end{align}\]

then consecutive values of \(MV\) are given by

(4)#\[\begin{align} MV_{k-1} & = \bar{MV} + K_P e_{k-1} + h K_I \sum_{j=0}^{k-1} e_{j} \\ MV_{k} & = \bar{MV} + K_P e_{k} + h K_I \sum_{j=0}^{k} e_{j} \end{align}\]

Taking differences gives a practical formula (velocity form) for updating the value of \(MV\) in response to measurements

(5)#\[\begin{align} MV_{k} & = MV_{k-1} + K_P(e_{k} - e_{k-1}) + h K_I e_{k} \end{align}\]

The following code defines a Python object that implements this algorithm.

def PI(Kp=1, Ki=0, MV_bar=0):
    ''' Simple proportion-integral controller

    Arguments:
        Kp: proportional gain
        Ki: integral gain
        MV_bar: bias
    
    '''
    # Initialize MV with bias
    MV = MV_bar

    # Initialize previous error as zero
    e_prev = 0

    # Define upper and lower bounds for MV (useful for later in this lab)
    MV_min = 0
    MV_max = 100

    # Infinite loop for the controller
    while True:
        # Yield MV and wait for new t_step, SP, and PV
        t_step, SP, PV = yield MV

        # Calculate error
        e = SP - PV

        # Apply the velocity form of the PI controller
        MV += Kp*(e - e_prev) + t_step*Ki*e

        # Apply upper and lower bounds
        # ** disabled for this simple controller, important later **
        # MV = max(MV_min, min(MV_max, MV)) # Apply upper and lower bounds

        # Store the current error for the next iteration
        e_prev = e

Exercise 1. Tune the PI control for the Temperature Control Lab#

The following cell provides an initial implementation of PI control for heater T1. This is setup for testing with the off-line simulation mode of tclab. Experiment with the simulation to find appropriate values for \(K_P\) and \(K_I\). Your design goal is to achieve the setpoint and stay within a zone of +/- 2 degrees as quickly as possible.

%matplotlib inline
from tclab import setup, clock, Historian, Plotter
# Suggestion: set to "if False:" after you have completed this 
# exercise to prevent yourself from accidently overwriting your results.
if False:

    # Debug this code with connected=False
    # Once you are happy with the code, proceed to the next exercise
    TCLab = setup(connected=False, speedup = 20)

    # Set point
    SP = 50 # deg C

    # Total time
    tfinal = 600 # seconds

    # Time step
    t_step = 1 # seconds

    # Create and intialize the PI controller
    controller1 = PI(Kp=12.0, Ki=0.1) 
    controller1.send(None)

    with TCLab() as lab:
        # Create the historian and plotter
        h = Historian(lab.sources)
        p = Plotter(h, tfinal)

        # Loop over time
        for t in clock(tfinal, t_step):

            # measure the the process variable
            PV = lab.T1

            # send the current time, set point, and process variable to the controller
            # compute the control action
            MV = controller1.send([t_step, SP, PV])

            # set the heat power (acutate)
            lab.Q1(MV)

            # Log the data                                      
            p.update()

        # This saves the results in case you accidently overwrite the experiment
        # After runnng this exercise, please find and save this csv flie
        h.to_csv("lab4-ex1.csv")

Exercise 2. Hardware testing the PI controller#

  • Copy and paste the above code into the cell below. Connect the code to the tclab hardware by changing connected=True. Adjust the experiment horizon to 1200 seconds to provide plenty of time for testing.

  • Test your controller.

  • After the controller has achieved the setpoint, introduce a disturbance. An example of a disturbance would be to increase air flow around the device (e.g., fan with a piece of paper), or to touch the heater with something thermally conductive (be careful, don’t use your finger. 50 deg C is hot enough to burn your skin.)

# Suggestion: set to False after you have completed this 
# exercise to prevent yourself from accidently overwriting your results.
if False: 

    TCLab = setup(connected=True)

    # Add your solution here

    with TCLab() as lab:
        # Create the historian and plotter
        h = Historian(lab.sources)
        p = Plotter(h, tfinal)
        # Add your solution here

        # This saves the results in case you accidently overwrite the experiment
        # After runnng this exercise, please find and save this csv flie
        h.to_csv("lab4-ex2.csv")

Discribe the distrubance you introduced. Note the time you introduced the distrubance. Is the distrubance visible in the data?

Does the performance match the simulation (Exercise 1)?

Answer:

Do you see any shortcomings in this control implementation with the naive “textbook” PI control? Please refer to specific aspects of your data to support your discussion.

Answer:

Exercise 3. Improved PI Control#

Now implmement the most practical version of the PI controller. You may optionally use the trick above to simulate the system to debug and then try it on your hardware.

# Optional: test your implementation here using connected=False and speedup=20 (or similar)

Test your controller with the TCLab hardware. Be sure to wait until your TCLab cools to room temperature before starting.

# Check that the device is at steady state 
# length of check
tfinal = 30 # seconds

# perform experiment
with TCLab() as lab:
    # No power to either heater
    lab.U1 = 0
    lab.U2 = 0

    # Initialize plotter and historian
    h = Historian(lab.sources)
    p = Plotter(h, tfinal)

    # Loop over time
    for t in clock(tfinal):
        
        # Take measurements and update plots
        p.update(t)
def PI_Bumpless(Kp=1, Ki=0, MV_bar=0):
    ''' Improved PI controller with:
        - Two anti-integral windup protections
        - Bumpless transition from manual to automatic control
    
    Arguments:
        Kp: proportional gain
        Ki: integral gain
        MV_bar: bias

    Usage:
        controller = PI_Bumpless(Kp, Ki, MV_bar)
        controller.send(None)
        MV = controller.send([t_step, SP, PV, MV])

    '''
    MV = MV_bar
    PV_prev = None
    MV_min = 0
    MV_max = 100
    while True:
        t_step, SP, PV, MV = yield MV
        # Add your solution here
        PV_prev = PV
# Suggestion: set to False after you've completed this exercise to prevent yourself from accidently overwriting your results.
if False: 

    # Add your solution here
    
        # This saves the results in case you accidently overwrite the experiment
        # After runnng this exercise, please find and save this csv flie
        h.to_csv("lab4-ex3.csv")

How do your results compare to the simple PI controller? Does this make sense? Use the mathematical formula for both controllers in your explanation.

Answer:

Excercise 4. Chocolate Tempering Revisited#

Revisit the chocolate tempering example described in Lab 2 (https://ndcbe.github.io/controls/assignments/Lab-2-Relay-Control.html). Implement this heating curve using PI control. Be sure to wait until your TCLab cools to room temperature before starting.

# Optional: test your implementation here using connected=False and speedup=20 (or similar)
# Check that the device is at steady state 
# experimental parameters
tfinal = 30

# perform experiment
with TCLab() as lab:
    lab.U1 = 0
    lab.U2 = 0
    h = Historian(lab.sources)
    p = Plotter(h, tfinal)
    for t in clock(tfinal):
        p.update(t)
# Suggestion: set to False after you've completed this exercise to prevent yourself from accidently overwriting your results.
if True: 
    # Add your solution here
            
        # This saves the results in case you accidently overwrite the experiment
        # After runnng this exercise, please find and save this csv flie
        h.to_csv("lab4-ex4.csv")

How do your results compare to the chocolate tempering curve using relay control? Does this make sense?

Answer:

Exercise 5. Stability Analysis Revisited#

Repeat the stability analysis for the PI controller that was performed in the notes (https://ndcbe.github.io/controls/notebooks/03.09-PI-Controller-Analysis.html). Identify values of Kp and Ki for (a) no oscillations and (b) with oscillations. Include a justification for these values in the Markdown cells. Make and reference plots in your justification.

import matplotlib.pyplot as plt
import numpy as np
from scipy.linalg import eig

# instructions:
# 1. update the parameters to match your TCLab device
# 2. moodify with code from our lecture

# parameters
T_amb = 21 # deg C <-- update based on your experimental data above
alpha = 0.00016 # watts / (units P1 * percent U1)
P1 = 100 # P1 units
U1 = 50 # steady state value of u1 (percent)

# fitted model parameters
# adjust using your values from Lab 3
Ua = 0.0261 # watts/deg C
Ub = 0.0222 # watts/deg C
CpH = 1.335 # joules/deg C
CpS = 1.328 # joules/deg C


def calc_eig(Kp,Ki,verbose=True):
    ''' Calculate the eigenvalues of the PI controller

    Arguments:
        Kp: proportional gain
        Ki: integral gain
        verbose: print the eigenvalues and eigenvectors

    Returns:
        w: vector of eigenvalues (complex numbers)
    
    '''

    # Add your solution here
# Create plots to analyze the stability using eigenvalues and eigenvectors

# Add your solution here

Analysis: Using the plot, identify values for \(K_P\) and \(K_I\) that (1) overshoot and oscillate (underdamped) and (b) do not oscillate (overdamped).

Overdamped: \(K_P = \), \(K_I= \) (from plot)

## Overdamped
#  Compute and print the eigenvalues and eigenvectors here

# Add your solution here

Justify why these Kp and Ki values are overdamped. Refer to the plot and the eigenvalues/eigenvectors printed above.

Answer:

Underdamped: \(K_P = \), \(K_I= \) (from plot).

# Underdamped
# Compute and print the eigenvalues and eigenvectors here

# Add your solution here

Justify why these Kp and Ki values are underdamped. Refer to the plot and the eigenvalues/eigenvectors printed above.

Answer:

Exercise 6. Simulate a Disturbance#

Using your both sets of controller gains from Exercise 5, simulate the system using a test case similar to class (https://ndcbe.github.io/controls/notebooks/03.06-Proportional-Integral-Control.html).

def SP(t):
    '''
    This function defines the setpoint for T1 as: 
        SP is 25 degC for t <= 20 seconds
        SP is 40 degC for t > 20 seconds 
    '''

    # Add your solution here

def DV(t):
    '''
    This function uses heat 2 as a disurbance:
        DV is 0% for t <= 200 seconds
        DV is 100% for t > 200 seconds
    '''
    # Add your solution here
def simulate_system(Kp,Ki):
    ''' Define function to simulate Exercise 6

    Arguments:
        Kp: proportional gain
        Ki: integral gain

    Action:
        Simulate the TCLab and creates plots
    
    '''

    # Create the TCLab object
    # you may use connected=False for debugging
    # but you need to conduct a real experiment for the final submission
    # TCLab = setup(connected=False, speedup=60)
    TCLab = setup(connected=True)

    # Add your solution here
# Underdamped system

# Add your solution here
# Overdamped system

# Add your solution here

How does the controller perform? Compare the eigenvalue analysis.

Answer:

Declarations#

Collaboration: If you worked with any classmates, please give their names here. Describe the nature of the collaboration.

Generative AI: If you used any Generative AI tools, please elaborate here.