Lab 5: Open Loop Optimization

Lab 5: Open Loop Optimization#

Learning Objectives#

Compute optimal heater strategy using physics-based model
Verify open loop performance with TC Lab hardware
Estimate unmeasured states (heater temperature) using model

Installation and Setup#

This notebook is using Pyomo and Ipopt. Be sure to follow the installation instructions on the class website.

Before beginning this lab, connect your TCLab device and run the following cell to estimate your ambient temperature.

# experimental parameters
tfinal = 30

# To complete the experiment, change this to True
# After completing the experiment, be sure to change it back to False
# This will prevent you from overwriting your data
run_tamb_tclab = False

from tclab import TCLab, clock, Historian, Plotter, setup
import numpy as np

if run_tamb_tclab:

    TCLab = setup(connected=True)

    # initialize list to store temperature
    T_list = []

    # perform experiment
    with TCLab() as lab:
        lab.U1 = 0
        lab.U2 = 0
        h = Historian(lab.sources)
        p = Plotter(h, tfinal)
        
        # save temperatures to access later
        T_list.append(lab.T1)
        T_list.append(lab.T2)

        for t in clock(tfinal):
            p.update(t)

    # estimate ambient temperature
    Tamb = np.mean(T_list) # deg C
    print(f"The ambient temperature is {Tamb} degrees C")

Exercise 0: Develop the Mathematical Model and Pyomo Simulation#

Recall the two-state model for a single heater/sensor assembly:

\[\begin{align*} C_p^H \frac{dT_H}{dt} & = U_a (T_{amb} - T_H) + U_b (T_S - T_H) + \alpha P u(t) + d(t)\\ C_p^S \frac{dT_S}{dt} & = - U_c (T_S - T_H) \end{align*}\]

Pyomo Model#

The code block below provides a library for performing several analysis tasks with our model. This is a class, which is a more sophisticated way to modularize our code. Briefly, classes store both data and methods to manipualte the data.

Take a few minutes to study this code. Note the missing lines of code to fill in when you get to Exercise 1. A few features:

It supports four modes:
- simulate: solve the model with all of the inputs specified (zero degrees of freedom)
- optimize: determine the best \(u(t)\) to track the desired \(T_{set}(t)\)
- observe: estimate the unmeasured state \(T_{H}(t)\) and disturbance \(d(t)\) from experimental data
- estimate: estimate the model parameters \(U_a\), \(U_b\), \(C_p^H\), \(C_p^S\) from experimental data
The __init__ method is called automatically when you create an instance of the object. The initial method is how you pass in data
The solve method call the numerical optimization algorithm
The plot method plots the Pyomo results
The get_Ts, get_Th, get_U, and get_D methods return data stored in the Pyomo model after it is solved

This notebook walks you through using the library. The take away message is that a class is a convient way to modularize our code and maximize reuse. The four supported modes are very similar.

# uncomment and update value from code cell and run above if returning to notebook with reset kernel
# Tamb = 20.75 # deg C

import matplotlib.pyplot as plt
from scipy import interpolate
import numpy as np

# ensure that IDAES is imported if there is an error calling Ipopt
import idaes

# Protip: only import the functions you need for each library.
# This tells other programs exactly where each function is defined.
# Avoid using from pyomo.environ import * as this imports everything
from pyomo.environ import ConcreteModel, Var, Param, Constraint, TransformationFactory, SolverFactory, Objective, minimize, value, Suffix
from pyomo.dae import DerivativeVar, ContinuousSet, Simulator

##### IMPORTANT #####
# Update the default values for Ua, Ub, CpH, and CpS to match your previous labs
#####################

class TCLabPyomo:
    '''
    This class contains the methods used for simulating the TCLab model, optimizing u(t) per a 
    given set point for T, observing a disturbance, and estimating model parameters in the TCLab.
    '''

    def __init__(self, mode, t_data, u_data, d_data, Tset_data, 
                 TS_data, Tamb = Tamb, alpha = 0.00016, P1 = 100,
                 Ua = 0.0261, Ub = 0.0222, CpH = 1.335, CpS = 1.328,
                 integrate_to_initialize=True, obj_weight_observe=0.1,
                 obj_weight_optimize=0.01
    ):
        '''
        This method is automatically called automatically when instantiating class.
        It stores input data in the class and creates the pyomo model.
        
        Arguments:
            mode: specify mode
            t_data: time data
            u_data: input control data
            d_data: disturbance data
            Tset_data: set point data
            TS_data: experimental data
            Tamb: ambient temperature, deg C
            alpha: watts / (units P1 * percent U1)
            P1: max power, P1 units
            Ua: heat transfer coefficient from heater to environment, watts/deg C
            Ub: heat transfer coefficient from heater to sensor, watts/deg C
            CpH: heat capacity of the heater, joules/deg C
            CpS: heat capacity of the sensor, joules/deg C
            integrate_to_initialize: Boolean
            obj_weight_observe: weight for disturbance in objective function, default is 0.1
            obj_weight_optimize: weight for heater temperature in objective function, default is 0.01

        '''
        
        # establish the valid operating modes
        valid_modes = ['simulate','optimize','estimate','observe']

        # raise an error if the user feeds in invalid operating mode
        if mode not in valid_modes:
            raise ValueError("'mode' must be one of the following:"+valid_modes)

        # define mode and data
        self.mode = mode
        self.t_data = t_data
        self.u_data = u_data
        self.d_data = d_data
        self.Tset_data = Tset_data
        self.TS_data = TS_data

        # set parameter values
        self.Tamb = Tamb
        self.Ua = Ua
        self.Ub = Ub
        self.CpH = CpH
        self.CpS = CpS
        self.P1 = P1
        self.alphaP = alpha*P1
        self.integrate_to_initialize = integrate_to_initialize
        self.obj_weight_observe = obj_weight_observe
        self.obj_weight_optimize = obj_weight_optimize

        # create the pyomo model
        self._create_pyomo_model()

        return None


    def _create_pyomo_model(self):
        '''
        Method that creates and defines the pyomo model for each mode.
        Arguments:
            None
        Returns:
            m: the pyomo model
        '''
        
        # create the pyomo model
        m = ConcreteModel()

        # create the time set
        m.t = ContinuousSet(initialize = self.t_data)  # make sure the experimental time grid are discretization points
        # define the heater and sensor temperatures as variables
        m.Th = Var(m.t, bounds=[0, 75], initialize=self.Tamb)
        m.Ts = Var(m.t, bounds=[0, 75], initialize=self.Tamb)

        def helper(my_array):
            '''
            Method that builds a dictionary to help initialization.
            Arguments:
                my_array: an array
            Returns:
                data: a dict {time: array_value}
            '''
            # ensure that the dimensions of array and time data match
            assert len(my_array) == len(self.t_data), "Dimension mismatch."
            data = {}
            for k,t in enumerate(self.t_data):
                data[t] = my_array[k]
            return data

        # for the simulate, observe, estimate modes
        if self.mode in ['simulate', 'observe', 'estimate']:
            # control decision is a parameter initialized with the input control data dict
            m.U = Param(m.t, initialize=helper(self.u_data), default = 0)
        else:
            # otherwise (optimize) control decision is a variable
            m.U = Var(m.t, bounds=(0, 100))

        # for the simulate, optimize, estimate modes
        if self.mode in ['simulate', 'optimize', 'estimate']:
            # if no distrubance data exists, initialize parameter at 0
            if self.d_data is None:
                 m.D = Param(m.t, default = 0)
            # otherwise initialize parameter with disturbance data dict
            else:
                m.D = Param(m.t, initialize=helper(self.d_data))
        # otherwise (observe) the disturbance is a variable
        else:
            m.D = Var(m.t)

        # define parameters that do not depend on mode
        m.Tamb = Param(initialize=self.Tamb)
        m.alphaP = Param(initialize=self.alphaP)

        # for the simulate, optimize, observe modes
        if self.mode in ['simulate', 'optimize', 'observe']:
            # Ua, Ub, CpH, and CpS are parameters
            m.Ua = Param(initialize=self.Ua)
            m.Ub = Param(initialize=self.Ub)
            # 1/CpH and 1/CpS parameters
            m.inv_CpH = Param(initialize=1/self.CpH)
            m.inv_CpS = Param(initialize=1/self.CpS)
        # otherwise (estimate) they are variables
        else:
            m.Ua = Var(initialize=self.Ua, bounds=(1E-5, 1.0))
            m.Ub = Var(initialize=self.Ub, bounds=(1E-5, 1.0))
            m.inv_CpH = Var(initialize=1/self.CpH, bounds=(1E-2, 100))
            m.inv_CpS = Var(initialize=1/self.CpS, bounds=(1E-1, 1000))

        # define variables for change in temperature wrt to time
        m.Thdot = DerivativeVar(m.Th, wrt = m.t)
        m.Tsdot = DerivativeVar(m.Ts, wrt = m.t)

        # define differential equations (model) as contraints
        # moved Cps to the right hand side to diagnose integrator
        m.Th_ode = Constraint(m.t, rule = lambda m, t: 
                            m.Thdot[t] == (m.Ua*(m.Tamb - m.Th[t]) + m.Ub*(m.Ts[t] - m.Th[t]) + m.alphaP*m.U[t] + m.D[t])*m.inv_CpH)
        m.Ts_ode = Constraint(m.t, rule = lambda m, t: 
                            m.Tsdot[t] == (m.Ub*(m.Th[t] - m.Ts[t]) )*m.inv_CpS)

        ##### Simulate to initialize
        # Specify time-varying input data
        # This is a dictionary of the form {time: value}
        # This is used only for the simulate mode,
        # which just initializes the model
        if self.integrate_to_initialize:

            m.var_input = Suffix(direction=Suffix.LOCAL)

            if self.u_data is not None:
                # initialize with data
                m.var_input[m.U] = helper(self.u_data)
            else:
                # otherwise initialize control decision of 0
                m.var_input[m.U] = {0:0}

            if self.d_data is not None:
                # initialize with data
                m.var_input[m.D] = helper(self.d_data)
            else:
                # otherwise initialize disturbance of 0
                m.var_input[m.D] = {0:0}
            

            # Simulate to initiialize
            # Makes the solver more efficient
            sim = Simulator(m, package='scipy') 
            tsim, profiles = sim.simulate(numpoints=100, 
                                        integrator='vode', 
                                        varying_inputs=m.var_input)
            sim.initialize_model()

        # for the optimize mode
        if self.mode == 'optimize':
            # Add requested constraints on U ramping here
            m.Udot = DerivativeVar(m.U, wrt = m.t)

            # Add your solution here

        # for the optimize mode, set point data is a parameter
        if self.mode == 'optimize':
            m.Tset = Param(m.t, initialize=helper(self.Tset_data))
            # otherwise, we are not using it

        # for the estimate and observe modes, experimental data is a parameter
        if self.mode in ['estimate', 'observe']:
            m.Ts_measure = Param(m.t, initialize=helper(self.TS_data))
            # otherwise, we are not using it

        # apply backward finite difference to the model
        TransformationFactory('dae.finite_difference').apply_to(m, scheme='BACKWARD',nfe=len(self.t_data)-1)

        if self.mode == 'optimize':
            # defining the tracking objective function
            m.obj = Objective(expr=sum( (m.Ts[t] - m.Tset[t])**2 + self.obj_weight_optimize*(m.Th[t] - m.Tset[t])**2 for t in m.t), sense=minimize)

        if self.mode == 'observe':
            # define observation (state estimation)
            m.obj = Objective(expr=sum((m.Ts[t] - m.Ts_measure[t])**2 + self.obj_weight_observe*m.D[t]**2 for t in m.t), sense=minimize)

        if self.mode == 'estimate':
            # define parameter estimation objective
            m.obj = Objective(expr=sum((m.Ts[t] - m.Ts_measure[t])**2 for t in m.t), sense=minimize)

        # initial conditions
        if self.TS_data is not None:
            # Initilize with first temperature measurement
            m.Ts[0].fix(self.TS_data[0])
            m.Th[0].fix(self.TS_data[0])
        else:
            # Otherwise, initialize with Tamb
            m.Th[0].fix(m.Tamb)
            m.Ts[0].fix(m.Tamb)

        self.m = m
    

    def solve(self):
        '''
        Solves the pyomo model using ipopt.
        '''
        solver = SolverFactory('ipopt')
        #solver.options['linear_solver'] = 'ma57'
        solver.solve(self.m, tee=True)


    def get_time(self):
        '''
        Returns time data from solved pyomo model.
        '''
        return self.t_data


    def get_Th(self):
        '''
        Returns heater temperature data from solved pyomo model.
        '''
        return np.array([value(self.m.Th[t]) for t in self.t_data])
    

    def get_Ts(self):
        '''
        Returns sensor temperature data from solved pyomo model.
        '''
        return np.array([value(self.m.Ts[t]) for t in self.t_data])
    

    def get_U(self):
        '''
        Returns control decision data from solved pyomo model.
        '''
        return np.array([value(self.m.U[t]) for t in self.t_data])
    

    def get_D(self):
        '''
        Returns disturbance data from solved pyomo model.
        '''
        return np.array([value(self.m.D[t]) for t in self.t_data])
    

    def get_parameters(self):
        '''
        Returns model parameters from solved pyomo model.
        '''
        return value(self.m.Ua), value(self.m.Ub), 1/value(self.m.inv_CpH), 1/value(self.m.inv_CpS)


    def print_parameters(self):
        '''
        Prints out the model parameters from solved pyomo model.
        '''
        Ua, Ub, CpH, CpS = self.get_parameters()

        print("The value of Ua is", round(Ua,4), "Watts/degC.")
        print("The value of Ub is", round(Ub,4), "Watts/degC.")
        print("The value of CpH is", round(CpH,3), "Joules/degC.")
        print("The value of CpS is", round(CpS,3),"Joules/degC.")


    def plot(self):
        '''
        Method to plot the results from the pyomo model.
        '''

        # extract predictions
        Th = self.get_Th()
        Ts = self.get_Ts()
        U = self.get_U()
        D = self.get_D()

        # create figure
        plt.figure(figsize=(10,6))

        # subplot 1: temperatures
        plt.subplot(3, 1, 1)
        if self.TS_data is not None:
            plt.scatter(self.t_data, self.TS_data, marker='.', label="$T_{S}$ measured", alpha=0.5,color='green')
        plt.plot(self.t_data, Th, label='$T_{H}$ predicted')
        plt.plot(self.t_data, Ts, label='$T_{S}$ predicted')
        if self.Tset_data is not None:
            plt.plot(self.t_data, self.Tset_data, label='$T_{set}$')

        plt.title('temperatures')
        plt.ylabel('deg C')
        plt.legend()
        plt.grid(True)

        # subplot 2: control decision
        plt.subplot(3, 1, 2)
        plt.plot(self.t_data, U)
        plt.title('heater power')
        plt.ylabel('percent of max')
        plt.grid(True)

        # subplot 3: disturbance
        plt.subplot(3, 1, 3)
        plt.plot(self.t_data, D)
        plt.title('disturbance')
        plt.ylabel('watts')
        plt.xlabel('time (s)')
        plt.grid(True)

        plt.tight_layout()
        plt.show()

Process Inputs#

The next cell defines some process inputs that will be used throughout the notebook to demonstrate aspects of process simulation, control, and estimation. These are gathered in one place to make it easier to modify the notebook to test the response under different conditions. These functions are implemented using the interp1d from the scipy library.

%matplotlib inline

tclab_disturbance = interpolate.interp1d(
    [ 0, 300, 400, 9999],   # time
    [ 0, 0, -.5, -.5],      # disturbance value
    fill_value="extrapolate")   # tolerates slight exptrapolation

tclab_input = interpolate.interp1d(
    [ 0, 50, 51, 450, 451, 9999],   # time
    [ 0,  0, 80,  80,   25,   25],  # input value
    fill_value="extrapolate")   # tolerates slight exptrapolation

tclab_setpoint = interpolate.interp1d(
    [0, 50, 150, 450, 550, 9999],   # time 
    [Tamb, Tamb, 60, 60, 35, 35],   # set point value
    fill_value="extrapolate")   # tolerates slight exptrapolation`

t_sim = np.linspace(0, 1000, 201)       # create 201 time points between 0 and 1000 seconds
u_sim = tclab_input(t_sim)              # calculate input signal at time points
d_sim = tclab_disturbance(t_sim)        # calculate disturbance at time points
setpoint_sim = tclab_setpoint(t_sim)    # calculate set point at time points

plt.figure(figsize=(10,6))
plt.subplot(3, 1, 1)
plt.plot(t_sim, setpoint_sim)
plt.title('setpoint')
plt.ylabel('deg C')
plt.grid(True)

plt.subplot(3, 1, 2)
plt.plot(t_sim, u_sim)
plt.title('heat power input')
plt.ylabel('percent of max')
plt.grid(True)

plt.subplot(3, 1, 3)
plt.plot(t_sim, d_sim)
plt.title('unmeasured disturbance')
plt.ylabel('watts')
plt.xlabel('time (s)')
plt.grid(True)

plt.tight_layout()

Let’s see how well our initial guess at a control strategy will work for us.

\[\begin{align*} C_p^H \frac{dT_H}{dt} & = U_a (T_{amb} - T_H) + U_b (T_S - T_H) + \alpha P u(t) + d(t)\\ C_p^S \frac{dT_S}{dt} & = - U_b (T_S - T_H) \end{align*}\]

subject to initial conditions

\[\begin{align*} T_H(t_0) & = T_{amb} \\ T_S(t_0) & = T_{amb} \end{align*}\]

and prior specification of inputs \(u(t)\) and \(d(t)\).

sim = TCLabPyomo('simulate',
                t_sim,
                u_sim,
                d_sim,
                setpoint_sim,
                None,
)

sim.solve()
sim.plot()

Discussion#

How well does the initial control strategy perform? Quantitatively justify your response using plots from the Pyomo model results.

Answer:

Optimal Control with Knowledge of Disturbances#

An optimal control policy minimizes the differences

\[\begin{align*} \min_{u} \int_{t_0}^{t_f} \left( \|T_S^{SP}(t) - T_S(t)\|^2\ + w \|T_S^{SP}(t) - T_H(t)\|^2\, \right) dt \\ \end{align*}\]

subject to constraints

\[\begin{align*} C_p^H \frac{dT_H}{dt} & = U_a (T_{amb} - T_H) + U_c (T_S - T_H) + P u(t) + d(t)\\ C_p^S \frac{dT_S}{dt} & = - U_c (T_S - T_H) \end{align*}\]

initial conditions

\[\begin{align*} T_H(t_0) & = T_{amb} \\ T_S(t_0) & = T_{amb} \end{align*}\]

and prior knowledge of \(d(t)\).

Exercise 1#

The optimal control computed above requires rapid changes in power level. In process systems where control action requires movement of a valve stem position, there are often limits on how fast the manipulated variable can change. Modify the model to include differential inequalities that limit the time rate of change of control to 0.25 \(s^{-1}\).

\[\begin{align*} \frac{du}{dt} & \leq \dot{u}_{max} \\ \frac{du}{dt} & \geq -\dot{u}_{max} \end{align*}\]

where \(\dot{u}_{max}\) is the maximum rate of change.

# Be sure to check and update the ambient temperature!
# Hint: you can do this with the hardware

opt = TCLabPyomo('optimize',
                 t_sim,
                 u_sim,
                 d_sim,
                 setpoint_sim,
                 None,
                 obj_weight_optimize=0.01   # keep default value
)

opt.solve()
opt.plot()

Discussion#

Using what you know about how the simulate and optimize modes operate, desctibe the differences in the controller performance and why those differences are occuring. How does the optimal control with knowledge of disturbances impact the controller performance? Quantitatively justify your response using plots from the Pyomo model results.

Answer:

Exercise 2#

Implement your calculated solution on the TC Lab hardware. Save the simulation results to a text file.

# Recommended: Check that the device is at steady state and ambient temperature
# experimental parameters
tfinal = 30

from tclab import TCLab, clock, Historian, Plotter, setup

# To complete the experiment, change this to True
# After completing the experiment, be sure to change it back to False
# This will prevent you from overwriting your data
run_exercise2_tclab = False

if run_exercise2_tclab:

    TCLab = setup(connected=True)

    # 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)

t_solution = opt.get_time()
u_solution = opt.get_U()
TS_solution = opt.get_Ts()
TH_solution = opt.get_Th()

u_open_loop = interpolate.interp1d(
    t_solution,   # time
    u_solution,   # control value
    fill_value="extrapolate")   # tolerates slight exptrapolation

TS_predict = interpolate.interp1d(
    t_solution,   # time
    TS_solution,   # sensor temperature prediction
    fill_value="extrapolate")   # tolerates slight exptrapolation

TH_predict = interpolate.interp1d(
    t_solution,   # time
    TH_solution,   # sensor temperature prediction
    fill_value="extrapolate")   # tolerates slight exptrapolation

from tclab import TCLab, clock, Historian, Plotter, setup

if run_exercise2_tclab:

    TCLab = setup(connected=True)

    t_final = t_solution[-1] # change this to 500 seconds for the actual experiment
    t_step = 0.1
    with TCLab() as lab:
        sources = [("T1", lambda: lab.T1), ("T2", lambda: lab.T2),
                ("T1pred", lambda: TS_predict(t)),
                ("SP1", lambda: tclab_setpoint(t)),
                ("Q1", lab.Q1), ("Q2", lab.Q2)]
        h = Historian(sources)
        p = Plotter(h, t_final, layout=(("T1", "SP1","T1pred"), ("T2",), ("Q1", "Q2")))
        for t in clock(t_final, t_step):
            U1 = u_open_loop(t)
            lab.Q1(U1)
            p.update()

import os.path

def save_tclab_data(h, file_name, overwrite_file=False):
    '''
    Save TCLab data to csv file
    Arguments:
        h: tclab historian objective
        file_name: valid file name as a string
        overwrite_file: bool, if True, overwrite exisiting file
            default is False to safeguard against accidentally rerunning this function
    '''

    if not overwrite_file and os.path.isfile('./'+file_name):
        raise FileExistsError(file_name + ' already exisits. Either choose a new filename or set overwrite_file = True.')
    else:
        h.to_csv(file_name)
        print("Successfully saved data to "+file_name)

    return file_name

excercise2_tclab_data_file = 'tclab-open-loop1.csv'

if run_exercise2_tclab:
    data_file = save_tclab_data(h,
                    excercise2_tclab_data_file, # tip: change this to end in "2" when copying code later in assignment
                    False
    )

Discussion#

How do your optimzized simulation results (with pyomo model) compare to the implementation results (with TCLab hardware)? Why do you think these differences do (or don’t) occur?

Answer:

Estimation/Observation#

“… and now my watch begins …” ―The Night’s Watch oath, Game of Thrones

The trouble with open-loop optimal control is that we cannot anticipate or know the values of unmeasured disturbances, much less the future values of those disturbances. The best we can do is use available data and process models to estimate the process state and disturbances. The state estimation problem is:

\[\begin{align*} \min_{\hat{T}_H, \hat{T}_S, \hat{d}} \int_{t - h}^t (\|\hat{T}_S(t) - T_S(t)\|^2\, + wd(t)^2 ) dt \\ \end{align*}\]

subject to

\[\begin{align*} C_p^H \frac{d\hat{T}_H}{dt} & = U_a (T_{amb} - \hat{T}_H) + U_c (\hat{T}_S - \hat{T}_H) + P u(t) + \hat{d}(t)\\ C_p^S \frac{d\hat{T}_S}{dt} & = - U_c (\hat{T}_S - \hat{T}_H) \end{align*}\]

and initial conditions

\[\begin{align*} T_H(t_0) & = T_{amb} \\ T_S(t_0) & = T_{amb} \end{align*}\]

Exercise 3#

Using the simulation data you saved in Exercise 2, estimate the heater temperature and disturbance when the disturbance has no impact on the objective function (i.e., \(w\) = 0).

# uncomment and run the following line to run exercise 3 if returning to notebook with a restarted kernel
# excercise2_tclab_data_file = 'tclab-open-loop1.csv'

import pandas as pd
experiment_data = pd.read_csv(excercise2_tclab_data_file)
experiment_data.head()

time_data = experiment_data['Time'].array # grab column from pandas, convert to array
time_data -= time_data[0] # make initial time 0

u_data = experiment_data['Q1'].array
T1_data = experiment_data['T1'].array

obs = TCLabPyomo('observe',
                 time_data,
                 u_data,
                 None,
                 None,
                 T1_data,
                 obj_weight_observe=0.1   # keep default value
)

obs.solve()
obs.plot()

Discussion#

Interpret your estimates for the heater temperature and disturbance.

Answer:

Parameter Estimation#

It is also possible our model does not match the physics of the system (e.g., wrong assumptions) or our estimated parameters are no longer valid or both. To investigate the later, we can repeat the parameter estimation problem for earlier in the semester:

\[\begin{align*} \min_{\hat{T}_H, \hat{T}_S, \hat{U_a}, \hat{U_b}, \hat{C_p^H}, \hat{C_p^S}} \int_{t - h}^t \|\hat{T}_S(t) - T_S(t)\|^2\,dt \\ \end{align*}\]

subject to

\[\begin{align*} \hat{C_p^H} \frac{d\hat{T}_H}{dt} & = \hat{U_a} (T_{amb} - \hat{T}_H) + \hat{U_b} (\hat{T}_S - \hat{T}_H) + \alpha P u(t) + d(t)\\ \hat{C_p^S} \frac{d\hat{T}_S}{dt} & = - \hat{U_b} (\hat{T}_S - \hat{T}_H) \end{align*}\]

and initial conditions

\[\begin{align*} T_H(t_0) & = T_{amb} \\ T_S(t_0) & = T_{amb} \end{align*}\]

and assumption disturbance is zero

(8)#\[\begin{equation} d(t) = 0 \end{equation}\]

Exercise 4#

Using data from your experiment, restimate the model parameters \(U_a\), \(U_b\), \(C_p^H\), \(C_p^S\).

NOTE: Depending on your operating system, the following code cell may take up to 2 minutes to solve.

est = TCLabPyomo('estimate',
                 time_data,
                 u_data,
                 None,
                 None,
                 T1_data)

est.solve()
est.plot()

Ua, Ub, CpH, CpS = est.get_parameters()

est.print_parameters()

Discussion#

Are these parameters close to what you expected (i.e., your Lab 3 results)? Why or why not?

Answer:

Exercise 5#

Compare the profiles of \(T_H\), \(T_S\), \(Q_1\), and \(d\) (disturbance) across:

The open loop optimization (prediction)
The hardware results
The state estimation

Write a short paragraph or a bulleted list of observations.

Your Observations:

Answer:

Compare the mean absolute and mean squared tracking error \(T_S - T_{set}\) and \(T_H - T_{set}\) across:

The open loop optimization
The state estimation (using hardware results)

# Hint: pay attention to how the set point was calculated and the size of your temperature data sets.
# Add your solution here

Write a short paragraph or a bulleted list of observations that interpret the above data.

Your Observations:

Answer:

Excerise 6. Chocolate Tempering#

Repeat the procedure from this lab for the dark chocolate tempering profile. You may optionally use the new model parameters you estimated in Exercise 4.

Open Loop Optimization#

# experimental parameters
tfinal = 30

# To complete the experiment, change this to True
# After completing the experiment, be sure to change it back to False
# This will prevent you from overwriting your data
run_tamb2_tclab = False

from tclab import TCLab, clock, Historian, Plotter, setup
import numpy as np

if run_tamb2_tclab:

    TCLab = setup(connected=True)

    # initialize list to store temperature
    T_list = []

    # perform experiment
    with TCLab() as lab:
        lab.U1 = 0
        lab.U2 = 0
        h = Historian(lab.sources)
        p = Plotter(h, tfinal)
        
        # save temperatures to access later
        T_list.append(lab.T1)
        T_list.append(lab.T2)

        for t in clock(tfinal):
            p.update(t)

    # estimate ambient temperature
    Tamb = np.mean(T_list) # deg C
    print(f"The ambient temperature is {Tamb} degrees C")

# uncomment and update value from code cell and run above if returning to notebook with reset kernel
# Tamb = 22.36 # deg C

# Be sure to check and update the ambient temperature!

# Determine the optimal control scheme
# Hint: the correct implementation of the set point can be found in Lab 2 solutions
# Add your solution here

Hardware Implementation#

# Recommended: Check that the device is at steady state and ambient temperature
# experimental parameters
tfinal = 30

# To complete the experiment, change this to True
# After completing the experiment, be sure to change it back to False
# This will prevent you from overwriting your data
run_exercise6_tclab = False

from tclab import TCLab, clock, Historian, Plotter, setup

if run_exercise6_tclab:

    TCLab = setup(connected=True)

    # 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)

# Implement optimal solution on your TCLab
## BEGIN SOLUTION

t_solution = choc_opt.get_time()
u_solution = choc_opt.get_U()
TS_solution = choc_opt.get_Ts()
TH_solution = choc_opt.get_Th()

u_open_loop = interpolate.interp1d(
    t_solution,   # time
    u_solution,   # control value
    fill_value="extrapolate")   # tolerates slight exptrapolation

TS_predict = interpolate.interp1d(
    t_solution,   # time
    TS_solution,   # sensor temperature prediction
    fill_value="extrapolate")   # tolerates slight exptrapolation

TH_predict = interpolate.interp1d(
    t_solution,   # time
    TH_solution,   # sensor temperature prediction
    fill_value="extrapolate")   # tolerates slight exptrapolation

from tclab import TCLab, clock, Historian, Plotter, setup

if run_exercise6_tclab:

    TCLab = setup(connected=True)

    t_final = t_chocolate[-1]
    t_step = 0.1
    with TCLab() as lab:
        sources = [("T1", lambda: lab.T1), ("T2", lambda: lab.T2),
                ("T1pred", lambda: TS_predict(t)),
                ("SP1", lambda: chocolate_setpoint(t)),
                ("Q1", lab.Q1), ("Q2", lab.Q2)]
        h = Historian(sources)
        p = Plotter(h, t_final, layout=(("T1", "SP1", "T1pred"), ("T2",), ("Q1", "Q2")))
        for t in clock(t_final, t_step):
            U1 = u_open_loop(t)
            lab.Q1(U1)
            p.update()

### END SOLUTION

# Save your experimental data as data_file2.

excercise6_tclab_data_file = 'tclab-open-loop2.csv'

# Add your solution here

Discussion#

How does the open loop optimization control scheme compare to the other implememntations we have seen this semester (e.g., relay, PI, etc.)? Justify your answer using knowlege of the mathematical models.

Answer:

Does your hardware implementation make sense based on your optimization results from the Pyomo model? Justify your answer.

Answer:

State Estimation#

# uncomment and run the following line to run exercise 6 if returning to notebook with a restarted kernel
# excercise6_tclab_data_file = 'tclab-open-loop2.csv'

experiment_data2 = pd.read_csv(excercise6_tclab_data_file)

time_data2 = experiment_data2['Time'].array # grab column from pandas, convert to array
time_data2 -= time_data2[0] # make initial time 0

# Estimate TH and d from data
# Add your solution here

Discussion#

Interpret your estimates for the heater temperature and disturbance. Do these makes sense? How do you know?

Answer:

Parameter Estimation#

NOTE: Depending on your operating system, the following code cell may take up to 2 minnutes to solve.

# Estimate model parameters
# Add your solution here

est2.print_parameters()

Discussion#

How do your parameters compare to those found in previous labs? Does this make sense?

Answer:

Lab 5: Open Loop Optimization

Contents

Lab 5: Open Loop Optimization#

Learning Objectives#

Installation and Setup#

Exercise 0: Develop the Mathematical Model and Pyomo Simulation#

Pyomo Model#

Process Inputs#

Discussion#

Optimal Control with Knowledge of Disturbances#

Exercise 1#

Discussion#

Exercise 2#

Discussion#

Estimation/Observation#

Exercise 3#

Discussion#

Parameter Estimation#

Exercise 4#

Discussion#

Exercise 5#

Excerise 6. Chocolate Tempering#

Open Loop Optimization#

Hardware Implementation#

Discussion#

State Estimation#

Discussion#

Parameter Estimation#

Discussion#