TCLab: Open-Loop Optimization and Estimation using Pyomo

6.9. TCLab: Open-Loop Optimization and Estimation using Pyomo#

6.9.1. Learning Objectives#

  1. Compute optimal heater strategy using physics-based model

  2. Verify open loop performance with TC Lab simulation

  3. Estimate unmeasured states (heater temperature) using model

import sys
if "google.colab" in sys.modules:
    !wget "https://raw.githubusercontent.com/IDAES/idaes-pse/main/scripts/colab_helper.py"
    import colab_helper
    colab_helper.install_idaes()
    colab_helper.install_ipopt()
    !pip install tclab

6.9.2. 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_b (T_S - T_H) \end{align*}\]

6.9.2.1. 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.

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 #####
# When you do Lab 5, you need to update Ua, Ub, CpH, and CpS to the values
# you previously estimated for your hardware.
#
#####################

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 = 21.0, alpha = 0.00016, P1 = 200,
                 Ua = 0.05, Ub = 0.05, CpH = 5.0, CpS = 1.0,
                 integrate_to_initialize=True, obj_weight=0.1
    ):
        '''
        Method called automatically when instantiating class.
        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, joules/deg C
            integrate_to_initialize: Boolean
            obj_weight: weight for heater temperature or disturbance in objective function, default is 0.1
        Returns:
            None 
        '''
        
        # 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 = obj_weight

        # 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
        # Increased temperature bounds to 85 deg C for digital twin mode
        m.Th = Var(m.t, bounds=[0, 85], initialize=self.Tamb)
        m.Ts = Var(m.t, bounds=[0, 85], 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 and Ub 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, 2.0))
            m.Ub = Var(initialize=self.Ub, bounds=(1E-5, 2.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()

        # In Lab 5, you will add constraints limiting the 
        # rate of change of the heater power (u) 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*(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*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()

6.9.2.2. 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

# Ambient temperature
# In Lab 5, you will update this based on the room's temperature
Tamb = 21.0 # deg C

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()
../_images/ffd55c6dd0c2e5d7c76d451000eb70662c789ee1f382b952c23c1463751f7c07.png

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

Ipopt 3.13.2: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...:     2400
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        0

Total number of variables............................:      802
                     variables with only lower bounds:        0
                variables with lower and upper bounds:      400
                     variables with only upper bounds:        0
Total number of equality constraints.................:      802
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 5.12e-01 0.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  0.0000000e+00 3.21e-15 1.38e+00  -1.0 3.49e+01    -  3.73e-01 1.00e+00h  1

Number of Iterations....: 1

                                   (scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   0.0000000000000000e+00    0.0000000000000000e+00
Constraint violation....:   3.2057689836051395e-15    3.2057689836051395e-15
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   3.2057689836051395e-15    3.2057689836051395e-15


Number of objective function evaluations             = 2
Number of objective gradient evaluations             = 2
Number of equality constraint evaluations            = 2
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 2
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 1
Total CPU secs in IPOPT (w/o function evaluations)   =      0.002
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.
../_images/474f58225f16b477ce9f34ab5319eee0b57f6ce2d3473cc1fadc2599fca50cc5.png

6.9.3. Optimal Control with Knowledge of Disturbances#

An optimal control policy minimizes the differences

\[\begin{align*} \min_{u} \int_{t_0}^{t_f} \|T_S^{SP}(t) - T_S(t)\|^2\ + w \|T_S^{SP}(t) - T_H(t)\|^2\,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)\).

opt = TCLabPyomo('optimize',
                 t_sim,
                 u_sim,
                 d_sim,
                 setpoint_sim,
                 None,
                 obj_weight=0.01
)

opt.solve()
opt.plot()

Ipopt 3.13.2: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...:     2601
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      400

Total number of variables............................:     1003
                     variables with only lower bounds:        0
                variables with lower and upper bounds:      601
                     variables with only upper bounds:        0
Total number of equality constraints.................:      802
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  4.3590553e+04 5.00e-01 7.80e+01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.7609161e+05 2.10e-15 4.92e+01  -1.0 4.52e+01    -  1.94e-02 1.00e+00h  1
   2  4.4312486e+04 2.69e-15 1.82e+01  -1.0 1.35e+02    -  4.73e-03 7.27e-01f  1
   3  4.5351328e+03 2.61e-15 4.43e+00  -1.0 7.82e+01    -  1.66e-01 8.00e-01f  1
   4  4.9045017e+02 2.28e-15 6.94e-01  -1.0 3.28e+01    -  2.55e-01 8.53e-01f  1
   5  2.0157778e+02 2.34e-15 4.42e-01  -1.0 2.40e+01    -  3.55e-01 1.00e+00f  1
   6  1.6575864e+02 2.22e-15 2.83e-01  -1.0 1.93e+01    -  4.28e-01 1.00e+00f  1
   7  1.5100971e+02 2.78e-15 1.57e-01  -1.0 1.88e+01    -  4.68e-01 1.00e+00f  1
   8  1.4418673e+02 2.23e-15 6.00e-02  -1.0 1.50e+01    -  6.40e-01 1.00e+00f  1
   9  1.4018857e+02 2.41e-15 1.19e-02  -1.0 2.89e+01    -  8.20e-01 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  1.3880583e+02 2.21e-15 2.88e-03  -1.0 3.13e+01    -  1.00e+00 8.55e-01f  1
  11  1.3550555e+02 2.19e-15 7.52e-15  -1.7 1.79e+01    -  1.00e+00 1.00e+00f  1
  12  1.3447858e+02 2.64e-15 3.64e-03  -2.5 1.30e+01    -  6.65e-01 1.00e+00f  1
  13  1.3424026e+02 2.16e-15 4.31e-04  -3.8 1.03e+01    -  9.26e-01 1.00e+00f  1
  14  1.3420930e+02 2.55e-15 7.12e-15  -3.8 6.39e+00    -  1.00e+00 1.00e+00f  1
  15  1.3419823e+02 2.71e-15 7.51e-06  -5.7 3.32e+00    -  9.63e-01 1.00e+00f  1
  16  1.3419787e+02 2.48e-15 9.68e-15  -5.7 5.11e+00    -  1.00e+00 1.00e+00f  1
  17  1.3419783e+02 2.28e-15 6.90e-15  -5.7 1.52e+01    -  1.00e+00 1.00e+00f  1
  18  1.3419773e+02 2.61e-15 7.21e-15  -8.6 1.45e-02    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 18

                                   (scaled)                 (unscaled)
Objective...............:   1.3419773060809698e+02    1.3419773060809698e+02
Dual infeasibility......:   7.2051086499059460e-15    7.2051086499059460e-15
Constraint violation....:   2.6088628220262246e-15    2.6088628220262246e-15
Complementarity.........:   8.7906857272874773e-09    8.7906857272874773e-09
Overall NLP error.......:   8.7906857272874773e-09    8.7906857272874773e-09


Number of objective function evaluations             = 19
Number of objective gradient evaluations             = 19
Number of equality constraint evaluations            = 19
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 19
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 18
Total CPU secs in IPOPT (w/o function evaluations)   =      0.013
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.
../_images/268c434f340c87a24a3ca6ce2f06e22078393084ce0620026c36b20475b06089.png

6.9.4. Hardware Implementation#

Implement your calculated solution on the TC Lab hardware (in “digital twin” mode). Save the simulation results to a text file.

def perform_hardware_simulation(opt):
    ''' 
    This function applies the optimized control to the hardware (digital twin).

    Arguments:
        opt: TCLabPyomo object

    Returns:
        h: Historian object
    '''

    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

    # In Lab 5, you will use connected=True
    # And run this experiment on your hardware
    TCLab = setup(connected=False, speedup=20)

    # Set final time to match simulation
    t_final = t_solution[-1]

    # Only implement the controller change every 2 seconds
    t_step = 2.0

    # Run experiments
    with TCLab() as lab:

        # Define sources for plotting
        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)]
        
        # Create a historian to store data
        h = Historian(sources)
        # Define plotting using same names as sources/historian
        p = Plotter(h, t_final, layout=(("T1", "SP1","T1pred"), ("T2",), ("Q1", "Q2")))

        lab.P1 = opt.P1

        # Loop over time steps
        for t in clock(t_final, t_step):
            # Look up the control value via interpolation
            U1 = u_open_loop(t)

            # Notice, we did not apply the distrubance here

            # Set control value to hardware
            lab.Q1(U1)

            # Update the plot
            p.update()

    return h

h = perform_hardware_simulation(opt)
../_images/414b8186bededf0c1291b53005f93371a199c75b92300fb94a6540b3fb9516b1.png 
TCLab Model disconnected successfully.
../_images/414b8186bededf0c1291b53005f93371a199c75b92300fb94a6540b3fb9516b1.png 
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 = './data/tclab-open-loop-digital-twin.csv'

data_file = save_tclab_data(h,
                excercise2_tclab_data_file,
                False
)

Successfully saved data to tclab-open-loop-digital-twin.csv

6.9.5. 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*}\]

Using the experimental data generated above, estimate the heater temperature and disturbance.

# local data file
excercise3_tclab_data_file = './data/tclab-open-loop-digital-twin.csv'

# saved data file
# excercise3_tclab_data_file = "https://raw.githubusercontent.com/ndcbe/controls/main/notebooks/data/tclab-open-loop-digital-twin.csv"

import pandas as pd
experiment_data = pd.read_csv(excercise3_tclab_data_file)
experiment_data.head()
Time T1 T2 T1pred SP1 Q1 Q2
0 3.521838 20.9495 20.9495 21.0 21.0 5.000000e+01 0
1 7.544108 20.9495 20.6272 21.0 21.0 9.800000e+00 0
2 9.603698 20.9495 20.9495 21.0 21.0 4.586005e-09 0
3 11.647745 20.9495 20.9495 21.0 21.0 4.746977e-09 0
4 15.576986 21.2718 20.9495 21.0 21.0 5.219476e-09 0 
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

6.9.5.1. Small Objective Weight (\(w=0.01\))#

obs = TCLabPyomo('observe',
                 time_data,
                 u_data,
                 None,
                 None,
                 T1_data,
                obj_weight=0.01
)

obs.solve()
obs.plot()

Ipopt 3.13.2: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...:     4941
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      761

Total number of variables............................:     1903
                     variables with only lower bounds:        0
                variables with lower and upper bounds:      760
                     variables with only upper bounds:        0
Total number of equality constraints.................:     1522
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  2.8926286e+03 5.70e-01 5.72e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.1442845e+01 5.82e-15 7.38e-02  -1.0 5.95e+00    -  7.90e-01 1.00e+00f  1
   2  2.2509680e+00 6.23e-15 2.98e-03  -1.0 1.61e+00    -  9.89e-01 1.00e+00f  1
   3  2.2460375e+00 5.71e-15 7.10e-15  -1.7 1.01e-01    -  1.00e+00 1.00e+00f  1
   4  2.2460008e+00 5.95e-15 7.10e-15  -3.8 4.27e-03    -  1.00e+00 1.00e+00f  1
   5  2.2460008e+00 5.61e-15 7.09e-15  -5.7 3.18e-05    -  1.00e+00 1.00e+00h  1
   6  2.2460008e+00 5.75e-15 7.09e-15  -8.6 3.71e-07    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 6

                                   (scaled)                 (unscaled)
Objective...............:   2.2460008023069951e+00    2.2460008023069951e+00
Dual infeasibility......:   7.0860308655916949e-15    7.0860308655916949e-15
Constraint violation....:   5.7547825012271847e-15    5.7547825012271847e-15
Complementarity.........:   2.5059270518034706e-09    2.5059270518034706e-09
Overall NLP error.......:   2.5059270518034706e-09    2.5059270518034706e-09


Number of objective function evaluations             = 7
Number of objective gradient evaluations             = 7
Number of equality constraint evaluations            = 7
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 7
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 6
Total CPU secs in IPOPT (w/o function evaluations)   =      0.010
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.
../_images/f777fb9c012370d399828059ade0ab962d9bc5a12c007d00211005515cef91c2.png

6.9.5.2. Medium Objective Weight (\(w=0.1\))#

obs = TCLabPyomo('observe',
                 time_data,
                 u_data,
                 None,
                 None,
                 T1_data,
                obj_weight=0.1
)

obs.solve()
obs.plot()

Ipopt 3.13.2: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...:     4941
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      761

Total number of variables............................:     1903
                     variables with only lower bounds:        0
                variables with lower and upper bounds:      760
                     variables with only upper bounds:        0
Total number of equality constraints.................:     1522
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  2.8926286e+03 5.70e-01 5.72e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.3103007e+01 7.83e-15 7.69e-02  -1.0 5.68e+00    -  7.83e-01 1.00e+00f  1
   2  4.7390881e+00 5.44e-15 1.38e-03  -1.0 4.96e-01    -  9.95e-01 1.00e+00f  1
   3  4.7384162e+00 5.49e-15 6.98e-15  -1.7 1.64e-02    -  1.00e+00 1.00e+00f  1
   4  4.7383858e+00 5.81e-15 7.08e-15  -3.8 1.51e-03    -  1.00e+00 1.00e+00f  1
   5  4.7383858e+00 5.22e-15 7.09e-15  -5.7 1.14e-05    -  1.00e+00 1.00e+00h  1
   6  4.7383858e+00 6.19e-15 7.13e-15  -8.6 1.40e-07    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 6

                                   (scaled)                 (unscaled)
Objective...............:   4.7383857509470761e+00    4.7383857509470761e+00
Dual infeasibility......:   7.1264405509564865e-15    7.1264405509564865e-15
Constraint violation....:   6.1884975170595391e-15    6.1884975170595391e-15
Complementarity.........:   2.5059147960052462e-09    2.5059147960052462e-09
Overall NLP error.......:   2.5059147960052462e-09    2.5059147960052462e-09


Number of objective function evaluations             = 7
Number of objective gradient evaluations             = 7
Number of equality constraint evaluations            = 7
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 7
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 6
Total CPU secs in IPOPT (w/o function evaluations)   =      0.009
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.
../_images/32edc29848bcf08926061b23e246a734e475221c8ac7b7756dc4d69a56b58be0.png

6.9.5.3. Medium Large Objective Weight (\(w=1\))#

obs = TCLabPyomo('observe',
                 time_data,
                 u_data,
                 None,
                 None,
                 T1_data,
                obj_weight=1.0
)

obs.solve()
obs.plot()

Ipopt 3.13.2: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...:     4941
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      761

Total number of variables............................:     1903
                     variables with only lower bounds:        0
                variables with lower and upper bounds:      760
                     variables with only upper bounds:        0
Total number of equality constraints.................:     1522
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  2.8926286e+03 5.70e-01 5.72e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  2.7708675e+01 5.85e-15 8.95e-02  -1.0 6.25e+00    -  7.68e-01 1.00e+00f  1
   2  1.9664761e+01 6.13e-15 1.10e-03  -1.0 3.55e-01    -  9.96e-01 1.00e+00f  1
   3  1.9664286e+01 5.47e-15 7.08e-15  -1.7 4.16e-03    -  1.00e+00 1.00e+00f  1
   4  1.9664258e+01 5.92e-15 7.07e-15  -3.8 8.71e-04    -  1.00e+00 1.00e+00f  1
   5  1.9664258e+01 5.92e-15 7.14e-15  -5.7 6.56e-06    -  1.00e+00 1.00e+00h  1
   6  1.9664258e+01 5.50e-15 7.10e-15  -8.6 8.07e-08    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 6

                                   (scaled)                 (unscaled)
Objective...............:   1.9664257877518683e+01    1.9664257877518683e+01
Dual infeasibility......:   7.1010527434483437e-15    7.1010527434483437e-15
Constraint violation....:   5.4950690716677132e-15    5.4950690716677132e-15
Complementarity.........:   2.5059104850534475e-09    2.5059104850534475e-09
Overall NLP error.......:   2.5059104850534475e-09    2.5059104850534475e-09


Number of objective function evaluations             = 7
Number of objective gradient evaluations             = 7
Number of equality constraint evaluations            = 7
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 7
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 6
Total CPU secs in IPOPT (w/o function evaluations)   =      0.010
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.
../_images/0c744bf0be5a09c5141f6488525d08b3d930440921443b9176c8e1101cb57e50.png

6.9.5.4. Extra Large Objective Weight (\(w=10\))#

obs = TCLabPyomo('observe',
                 time_data,
                 u_data,
                 None,
                 None,
                 T1_data,
                obj_weight=10.0
)

obs.solve()
obs.plot()

Ipopt 3.13.2: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...:     4941
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      761

Total number of variables............................:     1903
                     variables with only lower bounds:        0
                variables with lower and upper bounds:      760
                     variables with only upper bounds:        0
Total number of equality constraints.................:     1522
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  2.8926286e+03 5.70e-01 5.72e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.3478367e+02 5.77e-15 9.18e-02  -1.0 6.37e+00    -  7.66e-01 1.00e+00f  1
   2  1.2757689e+02 6.00e-15 1.24e-03  -1.0 2.83e-01    -  9.96e-01 1.00e+00f  1
   3  1.2757652e+02 5.44e-15 7.31e-15  -1.7 2.98e-03    -  1.00e+00 1.00e+00f  1
   4  1.2757649e+02 5.57e-15 8.45e-15  -3.8 5.30e-04    -  1.00e+00 1.00e+00f  1
   5  1.2757649e+02 5.77e-15 8.81e-15  -5.7 3.98e-06    -  1.00e+00 1.00e+00h  1
   6  1.2757649e+02 6.33e-15 7.25e-15  -8.6 4.92e-08    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 6

                                   (scaled)                 (unscaled)
Objective...............:   1.2757649007576410e+02    1.2757649007576410e+02
Dual infeasibility......:   7.2463869543232571e-15    7.2463869543232571e-15
Constraint violation....:   6.3313131075306715e-15    6.3313131075306715e-15
Complementarity.........:   2.5059079948903840e-09    2.5059079948903840e-09
Overall NLP error.......:   2.5059079948903840e-09    2.5059079948903840e-09


Number of objective function evaluations             = 7
Number of objective gradient evaluations             = 7
Number of equality constraint evaluations            = 7
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 7
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 6
Total CPU secs in IPOPT (w/o function evaluations)   =      0.010
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.
../_images/a8176ae8588e54db90e4e9b5a4f9cd33a168a2cd05f8122ba44ca2d42fc88703.png

6.9.6. 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 (optionally) the assumption disturbance is zero

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

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

NOTE: Depending on your computer’s hardward and the version of Ipopt installed, the following code cell may take up to 2 minutes to solve.

6.9.6.1. Assume no disturbance#

est = TCLabPyomo('estimate',
                 time_data,
                 u_data,
                 None, # assume d(t) = 0
                 None,
                 T1_data)

est.solve()
est.plot()

Ipopt 3.13.2: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...:     6462
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:     3043

Total number of variables............................:     1526
                     variables with only lower bounds:        0
                variables with lower and upper bounds:      764
                     variables with only upper bounds:        0
Total number of equality constraints.................:     1522
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  2.8926286e+03 5.70e-01 7.75e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  2.0120167e+03 4.65e-01 1.02e+01  -1.0 1.05e+01    -  4.05e-01 1.67e-01f  1
   2  1.0717968e+02 1.32e-01 7.78e+03  -1.0 6.65e+00    -  3.88e-02 1.00e+00f  1
   3  2.2025371e+01 1.71e-02 1.30e+03  -1.0 1.37e+00    -  9.48e-01 1.00e+00f  1
   4  1.3705448e+01 7.71e-03 2.86e+02  -1.0 9.40e-01    -  9.93e-01 1.00e+00f  1
   5  1.3645057e+01 5.79e-04 1.21e+01  -1.0 1.56e-01    -  1.00e+00 1.00e+00h  1
   6  1.3645108e+01 1.86e-01 3.13e+00  -1.0 1.10e+01    -  1.00e+00 1.00e+00f  1
   7  1.3645719e+01 3.44e-02 8.27e+00  -1.0 4.37e+00    -  1.00e+00 1.00e+00h  1
   8  1.3648763e+01 2.46e-02 7.50e-01  -1.0 4.68e-01    -  8.49e-01 1.00e+00H  1
   9  1.3648007e+01 3.03e-02 5.69e-01  -1.7 6.02e-01    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  1.3658732e+01 4.77e-03 2.40e-01  -1.7 3.44e-01  -2.0 1.00e+00 1.00e+00h  1
  11  1.3657205e+01 7.18e-03 5.00e-02  -1.7 3.54e-01    -  1.00e+00 1.00e+00h  1
  12  1.3657195e+01 1.99e-03 8.80e-02  -2.5 2.36e-01    -  1.00e+00 1.00e+00h  1
  13  1.3657194e+01 3.19e-03 4.04e-01  -2.5 4.98e+00    -  1.00e+00 3.60e-02h  5
  14  1.3657153e+01 6.82e-04 3.23e-03  -2.5 1.55e-01    -  1.00e+00 1.00e+00h  1
  15  1.3657153e+01 3.89e-04 6.66e-03  -3.8 1.14e-01    -  1.00e+00 1.00e+00h  1
  16  1.3657152e+01 1.26e-07 7.73e-06  -3.8 2.32e-03  -2.5 1.00e+00 1.00e+00h  1
  17  1.3657152e+01 1.30e-06 1.73e-05  -5.7 6.61e-03    -  1.00e+00 1.00e+00h  1
  18  1.3657152e+01 8.05e-07 3.77e-06  -8.6 5.23e-03    -  1.00e+00 1.00e+00h  1
  19  1.3657152e+01 5.63e-13 4.11e-09  -8.6 3.70e-06  -3.0 1.00e+00 1.00e+00h  1

Number of Iterations....: 19

                                   (scaled)                 (unscaled)
Objective...............:   1.3657151862902243e+01    1.3657151862902243e+01
Dual infeasibility......:   4.1069780432464661e-09    4.1069780432464661e-09
Constraint violation....:   5.6343818499726694e-13    5.6343818499726694e-13
Complementarity.........:   2.5059066030777639e-09    2.5059066030777639e-09
Overall NLP error.......:   4.1069780432464661e-09    4.1069780432464661e-09


Number of objective function evaluations             = 26
Number of objective gradient evaluations             = 20
Number of equality constraint evaluations            = 26
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 20
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 19
Total CPU secs in IPOPT (w/o function evaluations)   =      0.059
Total CPU secs in NLP function evaluations           =      0.005

EXIT: Optimal Solution Found.
../_images/0ca50f932fa0dad0aba5968f520226f8601e04cabdc40171fd48bb893e845d0e.png 
est.print_parameters()

The value of Ua is 0.0535 Watts/degC.
The value of Ub is 0.0148 Watts/degC.
The value of CpH is 6.911 Joules/degC.
The value of CpS is 0.318 Joules/degC.
Contents