From 2c632d29f9b84492fbbdf0397b2f4cfffc8e0f8b Mon Sep 17 00:00:00 2001
From: kuanhanl_DH4200 <kuanhanl@andrew.cmu.edu>
Date: Wed, 27 Apr 2022 21:15:49 -0400
Subject: [PATCH 01/26] construct mhe for pipeline

---
 nmpc_examples/mhe/mhe_constructor.py          | 234 ++++++++++
 .../pipeline/control_input_data.json          |   1 +
 nmpc_examples/pipeline/run_pipeline_mhe.py    | 412 ++++++++++++++++++
 3 files changed, 647 insertions(+)
 create mode 100644 nmpc_examples/mhe/mhe_constructor.py
 create mode 100644 nmpc_examples/pipeline/control_input_data.json
 create mode 100644 nmpc_examples/pipeline/run_pipeline_mhe.py

diff --git a/nmpc_examples/mhe/mhe_constructor.py b/nmpc_examples/mhe/mhe_constructor.py
new file mode 100644
index 0000000..27e2129
--- /dev/null
+++ b/nmpc_examples/mhe/mhe_constructor.py
@@ -0,0 +1,234 @@
+#################################################################################
+# The Institute for the Design of Advanced Energy Systems Integrated Platform
+# Framework (IDAES IP) was produced under the DOE Institute for the
+# Design of Advanced Energy Systems (IDAES), and is copyright (c) 2018-2021
+# by the software owners: The Regents of the University of California, through
+# Lawrence Berkeley National Laboratory,  National Technology & Engineering
+# Solutions of Sandia, LLC, Carnegie Mellon University, West Virginia University
+# Research Corporation, et al.  All rights reserved.
+#
+# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and
+# license information.
+#################################################################################
+from pyomo.core.base.set import Set
+from pyomo.core.base.var import Var
+from pyomo.core.base.constraint import Constraint
+from pyomo.core.base.componentuid import ComponentUID
+from pyomo.core.base.expression import Expression
+
+from pyomo.common.collections import ComponentSet
+
+from nmpc_examples.nmpc.dynamic_data.series_data import get_time_indexed_cuid
+
+
+def curr_sample_point(tp, sample_points):
+    for spt in sample_points:
+        if spt >= tp:
+            return spt
+
+
+def construct_measurement_variables_constraints(
+    # block,
+    sample_points,
+    variables,
+):
+
+    meas_set = Set(initialize=range(len(variables)))
+    meas_var = Var(meas_set, sample_points)
+    meas_error_var = Var(meas_set, sample_points, initialize=0.0)
+
+    def meas_con_rule(mod, index, sp):
+        return meas_var[index, sp] == \
+            variables[index][sp] + meas_error_var[index, sp]
+    meas_con = Constraint(
+        meas_set, sample_points, rule=meas_con_rule,
+    )
+
+    return meas_set, meas_var, meas_error_var, meas_con
+
+
+def construct_disturbed_model_constraints(
+    time,
+    sample_points,
+    mod_constraints,
+):
+
+    disturbance_set = Set(initialize=range(len(mod_constraints)))
+    disturbance_var = Var(disturbance_set, sample_points, initialize=0.0)
+
+    def disturbed_con_rule(m, index, i):
+        # try:
+        con_ele = mod_constraints[index][i]
+        # except KeyError:
+        #     if i == 0:
+        #         # Assuem we're using an implicit time discretization.
+        #         return Constraint.Skip
+        #     else:
+        #         raise KeyError(i)
+
+        if not con_ele.equality:
+            raise RuntimeError(con_ele.name, " is not an equality constraint.")
+        if con_ele.lower.value != con_ele.upper.value:
+            raise RuntimeError(con_ele_name, " has different lower and upper "
+                               "bounds.")
+
+        spt = curr_sample_point(i, sample_points)
+
+        return con_ele.body + disturbance_var[index, spt] == \
+            con_ele.upper.value
+    disturbed_con = Constraint(
+        disturbance_set, time, rule=disturbed_con_rule,
+    )
+
+    return disturbance_set, disturbance_var, disturbed_con
+
+
+def activate_disturbed_constraints_based_on_original_constraints(
+    time,
+    sample_points,
+    disturbance_var,
+    mod_constraints,
+    disturbed_con,
+):
+
+    # Deactivate original equalities and activate disturbed equalities
+    for index, con in enumerate(mod_constraints):
+        for tp in time:
+            con_ele = con[tp]
+            if con_ele.active:
+                con_ele.deactivate()
+                disturbed_con[index, tp].activate()
+            else:
+                # If the original equality is not active, deactivate the
+                # disturbed one.
+                disturbed_con[index, tp].deactivate()
+
+    # If all constraints in a specific sample time are not active, fix that
+    # model disturbance at zero.
+    spt_saw_tp = {spt:[] for spt in sample_points}
+    for tp in time:
+        spt = curr_sample_point(tp, sample_points)
+        spt_saw_tp[spt].append(tp)
+
+    for index in range(len(mod_constraints)):
+        for spt in sample_points:
+            con_active_list = [disturbed_con[index, tp].active
+                               for tp in spt_saw_tp[spt]
+                               ]
+            if not any(con_active_list):
+                disturbance_var[index, spt].fix(0.0)
+
+
+def get_error_disturbance_cost(
+    time,
+    sample_points,
+    components,
+    error_dist_var,
+    weight_data=None,
+):
+
+    component_cuids = [
+        get_time_indexed_cuid(comp, sets=(time,))
+        for comp in components
+    ]
+
+    def error_disturbance_rule(m, spt):
+        return sum(
+            weight_data[cuid] * error_dist_var[index, spt]**2
+            for index, cuid in enumerate(component_cuids)
+        )
+    error_disturbance_cost = Expression(
+        sample_points, rule=error_disturbance_rule
+    )
+
+    return error_disturbance_cost
+
+# def get_error_disturbance_cost(
+#     time,
+#     sample_points,
+#     measured_vars,
+#     meas_error_var,
+#     mod_cons,
+#     mod_dist_var,
+#     meas_error_weights,
+#     mod_dist_weights,
+# ):
+
+#     measured_var_cuids = [
+#         get_time_indexed_cuid(var, sets=(time,))
+#         for var in measured_vars
+#     ]
+
+#     mod_con_cuids = [
+#         get_time_indexed_cuid(con, sets=(time,))
+#         for con in mod_cons
+#     ]
+
+#     def error_disturbance_rule(m, spt):
+#         squared_error_sum = sum(
+#             meas_error_weights[cuid] * meas_error_var[index, spt]**2
+#             for index, cuid in enumerate(measured_var_cuids)
+#         )
+
+#         squared_dist_sum = sum(
+#             mod_dist_weights[cuid] * mod_dist_var[index, spt]**2
+#             for index, cuid in enumerate(mod_con_cuids)
+#         )
+
+#         return squared_error_sum + squared_dist_sum
+#     error_disturbance_cost = Expression(sample_points,
+#                                         rule=error_disturbance_rule)
+
+#     return error_disturbance_cost
+
+
+def get_squared_measurement_error(
+    measurement_estimate_map,
+    sample_points,
+    weight_data=None,
+):
+    """
+    This function returns a squared cost expression of measurement error for
+    the given mapping of meausrements to estimates.
+
+    Arguments
+    ---------
+    measurement_estimate_map: ComponentMap
+        Mapping of samplepoint-indexed measurements to estimates
+    sample_points: iterable
+        Set by which to index the squared expression
+    weight_data: dict
+        Optional. Maps measurement names to squared cost weights. If not
+        provided, weights of one are used.
+
+    Returns
+    -------
+    Pyomo Expression, indexed by sample points, containing the sum of weighted
+    squared difference between measurements and estimates, i.e., mearuerment
+    errors.
+
+    """
+    cuids = [
+        get_time_indexed_cuid(var, sets=(sample_points,))
+        for var in measurement_estimate_map.keys()
+    ]
+
+    if weight_data is None:
+        weight_data = {cuid: 1.0 for cuid in cuids}
+    for i, cuid in enumerate(cuids):
+        if cuid not in weight_data:
+            raise KeyError(
+                "Tracking weight dictionary does not contain a key for "
+                "variable\n%s with ComponentUID %s" % (variables[i].name, cuid)
+            )
+
+    def sqr_meas_error_rule(m, t):
+        return sum(
+            weight_data[cuid] * (meas[t] - esti[t])**2
+            for cuid, (meas, esti) in zip(cuids,
+                                          measurement_estimate_map.items())
+        )
+
+    sqr_measurement_error = Expression(sample_points,
+                                        rule=sqr_meas_error_rule)
+    return sqr_measurement_error
diff --git a/nmpc_examples/pipeline/control_input_data.json b/nmpc_examples/pipeline/control_input_data.json
new file mode 100644
index 0000000..3b13b87
--- /dev/null
+++ b/nmpc_examples/pipeline/control_input_data.json
@@ -0,0 +1 @@
+{"fs.pipeline.control_volume.flow_mass[*,1.0]": [553003.1393240632, 535472.4292275416], "fs.pipeline.control_volume.pressure[*,0]": [58.905880452751845, 58.92036532907178]}
\ No newline at end of file
diff --git a/nmpc_examples/pipeline/run_pipeline_mhe.py b/nmpc_examples/pipeline/run_pipeline_mhe.py
new file mode 100644
index 0000000..a109f5b
--- /dev/null
+++ b/nmpc_examples/pipeline/run_pipeline_mhe.py
@@ -0,0 +1,412 @@
+import pyomo.environ as pyo
+from pyomo.dae import ContinuousSet
+from pyomo.dae.flatten import flatten_dae_components
+
+from nmpc_examples.pipeline.pipeline_model import make_model
+
+from nmpc_examples.mhe.mhe_constructor import (
+    construct_measurement_variables_constraints,
+    construct_disturbed_model_constraints,
+    get_error_disturbance_cost,
+    activate_disturbed_constraints_based_on_original_constraints,
+)
+
+from nmpc_examples.nmpc.model_linker import DynamicVarLinker
+from nmpc_examples.nmpc.model_helper import DynamicModelHelper
+from nmpc_examples.nmpc.dynamic_data.series_data import TimeSeriesData
+
+import json
+import matplotlib.pyplot as plt
+
+"""
+Script to run an MHE simulation with a single natural gas pipeline model.
+"""
+
+
+def run_mhe(
+        simulation_horizon=20.0,
+        estimator_horizon=20.0,
+        sample_period=2.0,
+        ):
+    """
+    Run a simple MHE problem on a pipeline model.
+    """
+    #
+    # Make dynamic model
+    #
+    nxfe = 4
+    ntfe_per_sample = 4
+    n_cycles = round(simulation_horizon/sample_period)
+    m_plant = make_model(
+        horizon=sample_period,
+        ntfe=ntfe_per_sample,
+        nxfe=nxfe,
+    )
+    time = m_plant.fs.time
+    t0 = time.first()
+    tf = time.last()
+
+    ipopt = pyo.SolverFactory("ipopt")
+
+    # Fix "inputs" in dynamic model: inlet pressure and outlet flow rate
+    space = m_plant.fs.pipeline.control_volume.length_domain
+    x0 = space.first()
+    xf = space.last()
+    m_plant.fs.pipeline.control_volume.flow_mass[:, xf].fix()
+    m_plant.fs.pipeline.control_volume.pressure[:, x0].fix()
+
+    #
+    # Make steady model for past data
+    #
+    m_steady = make_model(dynamic=False, nxfe=nxfe)
+    m_steady.fs.pipeline.control_volume.flow_mass[:, x0].fix(
+        3.0e5*pyo.units.kg/pyo.units.hr
+    )
+    m_steady.fs.pipeline.control_volume.pressure[:, x0].fix(
+        57.0*pyo.units.bar
+    )
+
+    ipopt.solve(m_steady, tee=True)
+
+    #
+    # Extract data from steady state model
+    #
+    m_steady_helper = DynamicModelHelper(m_steady, m_steady.fs.time)
+    initial_data = m_steady_helper.get_data_at_time(t0)
+    scalar_data = m_steady_helper.get_scalar_variable_data()
+
+    #
+    # Load data into dynamic model
+    #
+    use_linker = False
+    if use_linker:
+        # If I want to use DynamicVarLinker:
+        # (E.g. if we didn't know that names in the steady model would
+        # be valid for the dynamic model.)
+        steady_scalar_vars, steady_dae_vars = flatten_dae_components(
+            m_steady, m_steady.fs.time, pyo.Var
+        )
+        steady_dae_vars_in_plant = [
+            m_plant.find_component(var.referent) for var in steady_dae_vars
+        ]
+        plant_steady_linker = DynamicVarLinker(
+            steady_dae_vars,
+            steady_dae_vars_in_plant,
+            m_steady.fs.time,
+            m_plant.fs.time,
+        )
+        plant_steady_linker.transfer()
+
+    else:
+        # If I want to use DynamicModelHelper:
+        m_plant_helper = DynamicModelHelper(m_plant, m_plant.fs.time)
+        m_plant_helper.load_scalar_data(scalar_data)
+        m_plant_helper.load_data_at_time(initial_data)
+
+    # Solve as a sanity check -- should be square with zero infeasibility
+    res = ipopt.solve(m_plant, tee=True)
+    pyo.assert_optimal_termination(res)
+
+    #
+    # Initialize data structure for simulation data
+    #
+    sim_data = m_plant_helper.get_data_at_time([t0])
+
+    #
+    # Construct dynamic model for estimator
+    #
+    samples_per_estimator = round(estimator_horizon/sample_period)
+    ntfe_per_estimator = ntfe_per_sample * samples_per_estimator
+    m_estimator = make_model(
+        horizon=estimator_horizon,
+        ntfe=ntfe_per_estimator,
+        nxfe=nxfe,
+    )
+    # Fix inputs at all time points
+    m_estimator.fs.pipeline.control_volume.pressure[:, x0].fix()
+    m_estimator.fs.pipeline.control_volume.flow_mass[:, xf].fix()
+
+    # Make helper object w.r.t. time
+    m_estimator_time_helper = DynamicModelHelper(
+        m_estimator, m_estimator.fs.time
+    )
+
+    #
+    # Make sample-point set for measurements and model disturbances
+    #
+    sample_points = [
+        t0 + sample_period*i for i in range(samples_per_estimator+1)
+    ]
+    m_estimator.fs.sample_points = ContinuousSet(initialize=sample_points)
+
+    #
+    # Construct components for measurements and measurement errors
+    #
+    m_estimator.estimation_block = pyo.Block()
+    esti_blo = m_estimator.estimation_block
+
+    cv = m_estimator.fs.pipeline.control_volume
+    measured_variables = [
+        pyo.Reference(cv.flow_mass[:, x]) for x in space.ordered_data()[:-1]
+    ]
+    measurement_info = construct_measurement_variables_constraints(
+        m_estimator.fs.sample_points,
+        measured_variables,
+    )
+    esti_blo.measurement_set = measurement_info[0]
+    esti_blo.measurement_variables = measurement_info[1]
+    # Measurement variables should be fixed all the time
+    esti_blo.measurement_variables.fix()
+    esti_blo.measurement_error_variables = measurement_info[2]
+    esti_blo.measurement_constraints = measurement_info[3]
+
+    #
+    # Construct disturbed model constraints
+    #
+    flatten_momentum_balance = [
+        pyo.Reference(cv.momentum_balance[:,x]) for x in space
+    ]
+    flatten_material_balance = [
+        pyo.Reference(cv.material_balances[:,x,"Vap","natural_gas"])
+        for x in space if x != space.last()
+    ]
+    model_constraints_to_be_disturbed = \
+        flatten_momentum_balance + flatten_material_balance
+
+    model_disturbance_info = construct_disturbed_model_constraints(
+        m_estimator.fs.time,
+        m_estimator.fs.sample_points,
+        model_constraints_to_be_disturbed,
+    )
+    esti_blo.disturbance_set = model_disturbance_info[0]
+    esti_blo.disturbance_variables = model_disturbance_info[1]
+    esti_blo.disturbed_constraints = model_disturbance_info[2]
+
+    activate_disturbed_constraints_based_on_original_constraints(
+        m_estimator.fs.time,
+        m_estimator.fs.sample_points,
+        esti_blo.disturbance_variables,
+        model_constraints_to_be_disturbed,
+        esti_blo.disturbed_constraints,
+    )
+
+    # Make helper object w.r.t. sample points
+    m_estimator_spt_helper = DynamicModelHelper(
+        m_estimator, m_estimator.fs.sample_points
+    )
+
+    #
+    # Construct least square objective to minimize measurement errors
+    # and model disturbances
+    #
+    measurement_error_weights = {
+        pyo.ComponentUID(var.referent): 1.0
+        for var in measured_variables
+    }
+    m_estimator.measurement_error_cost = get_error_disturbance_cost(
+        m_estimator.fs.time,
+        m_estimator.fs.sample_points,
+        measured_variables,
+        esti_blo.measurement_error_variables,
+        measurement_error_weights,
+    )
+
+    model_disturbance_weights = {
+        **{pyo.ComponentUID(con.referent): 10.0
+           for con in flatten_momentum_balance},
+        **{pyo.ComponentUID(con.referent): 20.0
+           for con in flatten_material_balance},
+    }
+    m_estimator.model_disturbance_cost = get_error_disturbance_cost(
+        m_estimator.fs.time,
+        m_estimator.fs.sample_points,
+        model_constraints_to_be_disturbed,
+        esti_blo.disturbance_variables,
+        model_disturbance_weights,
+    )
+
+    m_estimator.squred_error_disturbance_objective = pyo.Objective(
+        expr=(sum(m_estimator.measurement_error_cost.values()) +
+              sum(m_estimator.model_disturbance_cost.values())
+              )
+    )
+
+    #
+    # Initialize dynamic model to initial steady state
+    #
+    # TODO: load a time-series measurements & controls as initialization
+    m_estimator_time_helper.load_scalar_data(scalar_data)
+    m_estimator_time_helper.load_data_at_time(initial_data, m_estimator.fs.time)
+
+    for index, var in enumerate(measured_variables):
+        for spt in m_estimator.fs.sample_points:
+            esti_blo.measurement_variables[index, spt].set_value(var[spt].value)
+
+    #
+    # Initialize data sturcture for estimates
+    #
+    estimate_data = m_estimator_time_helper.get_data_at_time([t0])
+
+    #
+    # Set up a model linker to send measurements to estimator to update
+    # measurement variables
+    #
+    measured_variables_in_plant = [m_plant.find_component(var.referent)
+                                   for var in measured_variables
+    ]
+    flatten_measurements = [
+        pyo.Reference(esti_blo.measurement_variables[index, :])
+        for index in esti_blo.measurement_set
+    ]
+    measurement_linker = DynamicVarLinker(
+        measured_variables_in_plant,
+        flatten_measurements,
+    )
+
+    # set up a model linker to send measurements to estimator to initialize
+    # measured variables
+    #
+    estimate_linker = DynamicVarLinker(
+        measured_variables_in_plant,
+        measured_variables,
+    )
+
+    #
+    # Load control input data for simulation
+    #
+    with open("control_input_data.json", "r") as fr:
+        control_data = json.load(fr)
+
+    control_inputs = TimeSeriesData(
+        control_data, range(len(control_data)), time_set=None
+    )
+
+
+    for i in range(n_cycles):
+        # time.first() in the model corresponds to sim_t0 in "simulation time"
+        # time.last() in the model corresponds to sim_tf in "simulation time"
+        sim_t0 = i*sample_period
+        sim_tf = (i + 1)*sample_period
+
+        #
+        # Load inputs into plant
+        #
+        current_control = control_inputs.get_data_at_time_indices(indices=i)
+        non_initial_plant_time = list(m_plant.fs.time)[1:]
+        m_plant_helper.load_data_at_time(
+            current_control, non_initial_plant_time
+        )
+
+        res = ipopt.solve(m_plant, tee=True)
+        pyo.assert_optimal_termination(res)
+
+        #
+        # Extract time series data from solved model
+        #
+        # Initial conditions have already been accounted for.
+        # Note that this is only correct because we're using an implicit
+        # time discretization.
+        non_initial_time = list(m_plant.fs.time)[1:]
+        model_data = m_plant_helper.get_data_at_time(non_initial_time)
+
+        #
+        # Apply offset to data from model
+        #
+        model_data.shift_time_points(sim_t0)
+
+        #
+        # Extend simulation data with result of new simulation
+        #
+        sim_data.concatenate(model_data)
+
+        #
+        # Load measurements from plant to estimator
+        #
+        plant_tf = m_plant.fs.time.last()
+        estimator_tf = m_estimator.fs.time.last()
+        measurement_linker.transfer(plant_tf, estimator_tf)
+
+        #
+        # Initialize measured variables within the last sample period to
+        # current measurements
+        #
+        last_sample_period = list(m_estimator.fs.time)[-4:]
+        # for index, var in enumerate(measured_variables):
+        #     for tp in last_sampel_period:
+        #         var[tp].set_value(blo.measurement_variables[index, estimator_tf])
+        estimate_linker.transfer(tf, last_sample_period)
+
+        # Load inputs to estimator
+        m_estimator_time_helper.load_data_at_time(
+            current_control, last_sample_period
+        )
+
+        res = ipopt.solve(m_estimator, tee=True)
+        pyo.assert_optimal_termination(res)
+
+        #
+        # Extract estimate data from estimator
+        #
+        estimator_data = m_estimator_time_helper.get_data_at_time([estimator_tf])
+        # Shift time points from "estimator time" to "simulation time"
+        estimator_data.shift_time_points(sim_tf-estimator_tf)
+
+        #
+        # Extend data structure of estimates
+        #
+        estimate_data.concatenate(estimator_data)
+
+        #
+        # Re-initialize estimator model
+        #
+        m_estimator_time_helper.shift_values_by_time(sample_period)
+        m_estimator_spt_helper.shift_values_by_time(sample_period)
+
+        #
+        # Re-initialize model to final values.
+        # This sets new initial conditions, including inputs.
+        #
+        m_plant_helper.copy_values_at_time(source_time=plant_tf)
+
+    return sim_data, estimate_data
+
+
+def plot_states_estimates_from_data(
+        state_data,
+        estimate_data,
+        names,
+        show=False
+        ):
+    state_time = state_data.get_time_points()
+    states = state_data.get_data()
+    estimate_time = estimate_data.get_time_points()
+    estimates = estimate_data.get_data()
+    for i, name in enumerate(names):
+        fig, ax = plt.subplots()
+        state_values = states[name]
+        estimate_values = estimates[name]
+        ax.plot(state_time, state_values, label="Plant states")
+        ax.plot(estimate_time, estimate_values, "o", label="Estimates")
+        ax.set_title(name)
+        ax.set_xlabel("Time (hr)")
+        ax.legend()
+
+        if show:
+            fig.show()
+        else:
+            fname = "state_estimate%s.png" % i
+            fig.savefig(fname, transparent=False)
+
+
+if __name__ == "__main__":
+    simulation_data, estimate_data = run_mhe(
+        simulation_horizon=4.0,
+    )
+    plot_states_estimates_from_data(
+        simulation_data,
+        estimate_data,
+        [
+            pyo.ComponentUID("fs.pipeline.control_volume.flow_mass[*,0.0]"),
+            pyo.ComponentUID("fs.pipeline.control_volume.pressure[*,1.0]"),
+        ],
+    )

From 61f3ff581a91ee10f7d4ecd5fa0268afe62966d6 Mon Sep 17 00:00:00 2001
From: kuanhanl_DH4200 <kuanhanl@andrew.cmu.edu>
Date: Thu, 28 Apr 2022 14:14:57 -0400
Subject: [PATCH 02/26] add function descriptions

---
 nmpc_examples/mhe/mhe_constructor.py | 101 +++++++++++++++++++++++++--
 1 file changed, 95 insertions(+), 6 deletions(-)

diff --git a/nmpc_examples/mhe/mhe_constructor.py b/nmpc_examples/mhe/mhe_constructor.py
index 27e2129..6d7f4a9 100644
--- a/nmpc_examples/mhe/mhe_constructor.py
+++ b/nmpc_examples/mhe/mhe_constructor.py
@@ -16,23 +16,62 @@
 from pyomo.core.base.componentuid import ComponentUID
 from pyomo.core.base.expression import Expression
 
-from pyomo.common.collections import ComponentSet
-
 from nmpc_examples.nmpc.dynamic_data.series_data import get_time_indexed_cuid
 
 
 def curr_sample_point(tp, sample_points):
+    """
+    This function returns the sample point to which the given time element
+    belongs.
+
+    Arguments
+    ---------
+    tp: float
+        Time element
+    sample_points: iterable
+        Set of sample points
+
+    Returns
+    -------
+    A sample point to which the given time element belongs.
+
+    """
     for spt in sample_points:
         if spt >= tp:
             return spt
 
 
 def construct_measurement_variables_constraints(
-    # block,
     sample_points,
     variables,
 ):
+    """
+    This function constructs components for measurments, including:
+    1. Index set for measurements
+    2. Variables for measurements and measurement errors
+    2. Constraints for measured variables, measurements, and measurement
+    errors.
+
+    Arguments
+    ---------
+    sample_points: iterable
+        Set of sample points
+    variables: list
+        List of time-indexed variables that are measured
 
+    Returns
+    -------
+    meas_set: Pyomo Set
+        Index set for measurements
+    meas_var: Pyomo Var
+        Measurement variable, indexed by meas_set and sample_points
+    meas_error_var: Pyomo Var
+        Measurement error variable, indexed by meas_set and sample_points
+    meas_con: Pyomo Constraint
+        measurement constraint, indexed by meas_set and sample_points,
+            **measurement == measured_var + measurement_error**
+
+    """
     meas_set = Set(initialize=range(len(variables)))
     meas_var = Var(meas_set, sample_points)
     meas_error_var = Var(meas_set, sample_points, initialize=0.0)
@@ -52,7 +91,34 @@ def construct_disturbed_model_constraints(
     sample_points,
     mod_constraints,
 ):
+    """
+    This function constructs components for model disturbances, including:
+    1. Index set for model disturbances
+    2. Variables for model disturbances
+    3. Disturbed model constraints, consisting of original model constraints
+    and model disturbances.
+
+    Arguments
+    ----------
+    time: iterable
+        Set by which to index model constraints
+    sample_points: iterable
+        Set of sample points
+    mod_constraints : list
+        List of model constraints to add model disturbances
 
+    Returns
+    -------
+    disturbance_set: Pyomo Set
+        Index set of model disturbances
+    disturbance_var: Pyomo Var
+        Model disturbance variable, indexed by disturbance_set and time
+    disturbed_con: Pyomo Constraint
+        Model constraints with model disturbances, indexed by disturbance_set
+        and time,
+            ** original_constraint + disturbance == 0.0 **
+
+    """
     disturbance_set = Set(initialize=range(len(mod_constraints)))
     disturbance_var = Var(disturbance_set, sample_points, initialize=0.0)
 
@@ -71,6 +137,9 @@ def disturbed_con_rule(m, index, i):
         if con_ele.lower.value != con_ele.upper.value:
             raise RuntimeError(con_ele_name, " has different lower and upper "
                                "bounds.")
+        if con_ele.upper.value != 0.0:
+            raise RuntimeError(con_ele_name, " is an equality but its bound "
+                               "is not zero.")
 
         spt = curr_sample_point(i, sample_points)
 
@@ -90,7 +159,26 @@ def activate_disturbed_constraints_based_on_original_constraints(
     mod_constraints,
     disturbed_con,
 ):
+    """
+    This function activate the model constraint and deactivate the original
+    constraint, if the original constarint is active. Also, if all model
+    constraints within a specific sample period are not active, fix the
+    disturbance variable at zero.
+
+    Parameters
+    ----------
+    time: iterable
+        Set which indexes model constraints
+    sample_points: iterable
+        Set of sample points
+    disturbance_var: Pyomo Var
+        Model disturbances
+    mod_constraints: list
+        List of original model constraints
+    disturbed_con: Pyomo Constraint
+        Model constraints with model disturbances
 
+    """
     # Deactivate original equalities and activate disturbed equalities
     for index, con in enumerate(mod_constraints):
         for tp in time:
@@ -112,9 +200,9 @@ def activate_disturbed_constraints_based_on_original_constraints(
 
     for index in range(len(mod_constraints)):
         for spt in sample_points:
-            con_active_list = [disturbed_con[index, tp].active
-                               for tp in spt_saw_tp[spt]
-                               ]
+            con_active_list = [
+                disturbed_con[index, tp].active for tp in spt_saw_tp[spt]
+            ]
             if not any(con_active_list):
                 disturbance_var[index, spt].fix(0.0)
 
@@ -143,6 +231,7 @@ def error_disturbance_rule(m, spt):
 
     return error_disturbance_cost
 
+
 # def get_error_disturbance_cost(
 #     time,
 #     sample_points,

From 1edaf7e11e81b43bf354d1153a6b3bec88a71278 Mon Sep 17 00:00:00 2001
From: kuanhanl_DH4200 <kuanhanl@andrew.cmu.edu>
Date: Fri, 29 Apr 2022 11:21:38 -0400
Subject: [PATCH 03/26] add description for the function of adding cost
 expression

---
 nmpc_examples/mhe/mhe_constructor.py | 41 ++++++++++++++++++++++++++--
 1 file changed, 39 insertions(+), 2 deletions(-)

diff --git a/nmpc_examples/mhe/mhe_constructor.py b/nmpc_examples/mhe/mhe_constructor.py
index 6d7f4a9..619872c 100644
--- a/nmpc_examples/mhe/mhe_constructor.py
+++ b/nmpc_examples/mhe/mhe_constructor.py
@@ -168,7 +168,7 @@ def activate_disturbed_constraints_based_on_original_constraints(
     Parameters
     ----------
     time: iterable
-        Set which indexes model constraints
+        Time set which indexes model constraints
     sample_points: iterable
         Set of sample points
     disturbance_var: Pyomo Var
@@ -214,12 +214,50 @@ def get_error_disturbance_cost(
     error_dist_var,
     weight_data=None,
 ):
+    """
+    #TODO: This function is similar to "get_tracking_cost_from_constant_setpoint",
+    #Should I merge this one to that one?
+    This function return a squared cost expression for measurement errors or
+    model disturbances.
+
+    Arguments
+    ---------
+    time: iterable
+        Time set which indexes components (measured variables or
+        model constraints)
+    sample_points: iterable
+        Set of sample points
+    components: List
+        List of components (measured variables or model constraints)
+    error_dist_var: Pyomo Var
+        Variable of measurement error
+        from "construct_measurement_variables_constraints"
+        or model disturbance
+        from "construct_disturbed_model_constraints"
+    weight_data: dict
+        Optional. Maps variable names to squared cost weights. If not provided,
+        weights of one are used.
+
+    Returns
+    -------
+    Pyomo Expression, indexed by sample_points, containing the sum of squared
+    errors or disturbances.
 
+    """
     component_cuids = [
         get_time_indexed_cuid(comp, sets=(time,))
         for comp in components
     ]
 
+    if weight_data is None:
+        weight_data = {cuid: 1.0 for cuid in cuids}
+    for i, cuid in enumerate(cuids):
+        if cuid not in weight_data:
+            raise KeyError(
+                "Error/disturbance weight dictionary does not contain a key for "
+                "variable\n%s with ComponentUID %s" % (variables[i].name, cuid)
+            )
+
     def error_disturbance_rule(m, spt):
         return sum(
             weight_data[cuid] * error_dist_var[index, spt]**2
@@ -228,7 +266,6 @@ def error_disturbance_rule(m, spt):
     error_disturbance_cost = Expression(
         sample_points, rule=error_disturbance_rule
     )
-
     return error_disturbance_cost
 
 

From f11bd449ef1c2b2a0e2915e70ba387ac831fd4d7 Mon Sep 17 00:00:00 2001
From: kuanhanl_DH4200 <kuanhanl@andrew.cmu.edu>
Date: Fri, 29 Apr 2022 11:23:28 -0400
Subject: [PATCH 04/26] delete unused codes

---
 nmpc_examples/mhe/mhe_constructor.py | 91 ----------------------------
 1 file changed, 91 deletions(-)

diff --git a/nmpc_examples/mhe/mhe_constructor.py b/nmpc_examples/mhe/mhe_constructor.py
index 619872c..77d6cce 100644
--- a/nmpc_examples/mhe/mhe_constructor.py
+++ b/nmpc_examples/mhe/mhe_constructor.py
@@ -267,94 +267,3 @@ def error_disturbance_rule(m, spt):
         sample_points, rule=error_disturbance_rule
     )
     return error_disturbance_cost
-
-
-# def get_error_disturbance_cost(
-#     time,
-#     sample_points,
-#     measured_vars,
-#     meas_error_var,
-#     mod_cons,
-#     mod_dist_var,
-#     meas_error_weights,
-#     mod_dist_weights,
-# ):
-
-#     measured_var_cuids = [
-#         get_time_indexed_cuid(var, sets=(time,))
-#         for var in measured_vars
-#     ]
-
-#     mod_con_cuids = [
-#         get_time_indexed_cuid(con, sets=(time,))
-#         for con in mod_cons
-#     ]
-
-#     def error_disturbance_rule(m, spt):
-#         squared_error_sum = sum(
-#             meas_error_weights[cuid] * meas_error_var[index, spt]**2
-#             for index, cuid in enumerate(measured_var_cuids)
-#         )
-
-#         squared_dist_sum = sum(
-#             mod_dist_weights[cuid] * mod_dist_var[index, spt]**2
-#             for index, cuid in enumerate(mod_con_cuids)
-#         )
-
-#         return squared_error_sum + squared_dist_sum
-#     error_disturbance_cost = Expression(sample_points,
-#                                         rule=error_disturbance_rule)
-
-#     return error_disturbance_cost
-
-
-def get_squared_measurement_error(
-    measurement_estimate_map,
-    sample_points,
-    weight_data=None,
-):
-    """
-    This function returns a squared cost expression of measurement error for
-    the given mapping of meausrements to estimates.
-
-    Arguments
-    ---------
-    measurement_estimate_map: ComponentMap
-        Mapping of samplepoint-indexed measurements to estimates
-    sample_points: iterable
-        Set by which to index the squared expression
-    weight_data: dict
-        Optional. Maps measurement names to squared cost weights. If not
-        provided, weights of one are used.
-
-    Returns
-    -------
-    Pyomo Expression, indexed by sample points, containing the sum of weighted
-    squared difference between measurements and estimates, i.e., mearuerment
-    errors.
-
-    """
-    cuids = [
-        get_time_indexed_cuid(var, sets=(sample_points,))
-        for var in measurement_estimate_map.keys()
-    ]
-
-    if weight_data is None:
-        weight_data = {cuid: 1.0 for cuid in cuids}
-    for i, cuid in enumerate(cuids):
-        if cuid not in weight_data:
-            raise KeyError(
-                "Tracking weight dictionary does not contain a key for "
-                "variable\n%s with ComponentUID %s" % (variables[i].name, cuid)
-            )
-
-    def sqr_meas_error_rule(m, t):
-        return sum(
-            weight_data[cuid] * (meas[t] - esti[t])**2
-            for cuid, (meas, esti) in zip(cuids,
-                                          measurement_estimate_map.items())
-        )
-
-    sqr_measurement_error = Expression(sample_points,
-                                        rule=sqr_meas_error_rule)
-    return sqr_measurement_error

From 79ce495e2513c9cccc08300b7a4090168f806aa1 Mon Sep 17 00:00:00 2001
From: kuanhanl_DH4200 <kuanhanl@andrew.cmu.edu>
Date: Sun, 1 May 2022 14:44:33 -0400
Subject: [PATCH 05/26] fix a bug

---
 nmpc_examples/mhe/mhe_constructor.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/nmpc_examples/mhe/mhe_constructor.py b/nmpc_examples/mhe/mhe_constructor.py
index 77d6cce..ad06e7e 100644
--- a/nmpc_examples/mhe/mhe_constructor.py
+++ b/nmpc_examples/mhe/mhe_constructor.py
@@ -251,7 +251,7 @@ def get_error_disturbance_cost(
 
     if weight_data is None:
         weight_data = {cuid: 1.0 for cuid in cuids}
-    for i, cuid in enumerate(cuids):
+    for i, cuid in enumerate(component_cuids):
         if cuid not in weight_data:
             raise KeyError(
                 "Error/disturbance weight dictionary does not contain a key for "

From 5c290576a189619dbecb5606969a6e3a62d8a24c Mon Sep 17 00:00:00 2001
From: kuanhanl_DH4200 <kuanhanl@andrew.cmu.edu>
Date: Sun, 1 May 2022 14:45:09 -0400
Subject: [PATCH 06/26] redefine time points for control inputs

---
 nmpc_examples/pipeline/run_pipeline_mhe.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/nmpc_examples/pipeline/run_pipeline_mhe.py b/nmpc_examples/pipeline/run_pipeline_mhe.py
index a109f5b..0ffb890 100644
--- a/nmpc_examples/pipeline/run_pipeline_mhe.py
+++ b/nmpc_examples/pipeline/run_pipeline_mhe.py
@@ -277,11 +277,11 @@ def run_mhe(
     with open("control_input_data.json", "r") as fr:
         control_data = json.load(fr)
 
+    control_input_time = [i*sample_period for i in range(len(control_data))]
     control_inputs = TimeSeriesData(
-        control_data, range(len(control_data)), time_set=None
+        control_data, control_input_time, time_set=None
     )
 
-
     for i in range(n_cycles):
         # time.first() in the model corresponds to sim_t0 in "simulation time"
         # time.last() in the model corresponds to sim_tf in "simulation time"
@@ -291,7 +291,7 @@ def run_mhe(
         #
         # Load inputs into plant
         #
-        current_control = control_inputs.get_data_at_time_indices(indices=i)
+        current_control = control_inputs.get_data_at_time(time=sim_t0)
         non_initial_plant_time = list(m_plant.fs.time)[1:]
         m_plant_helper.load_data_at_time(
             current_control, non_initial_plant_time

From 4c4519b8ccd3273d6453b4418051dbebb325a6b3 Mon Sep 17 00:00:00 2001
From: kuanhanl_DH4200 <kuanhanl@andrew.cmu.edu>
Date: Sun, 1 May 2022 15:48:41 -0400
Subject: [PATCH 07/26] import control data with os module

---
 nmpc_examples/pipeline/run_pipeline_mhe.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/nmpc_examples/pipeline/run_pipeline_mhe.py b/nmpc_examples/pipeline/run_pipeline_mhe.py
index 0ffb890..e34252e 100644
--- a/nmpc_examples/pipeline/run_pipeline_mhe.py
+++ b/nmpc_examples/pipeline/run_pipeline_mhe.py
@@ -16,6 +16,7 @@
 from nmpc_examples.nmpc.dynamic_data.series_data import TimeSeriesData
 
 import json
+import os
 import matplotlib.pyplot as plt
 
 """
@@ -274,7 +275,9 @@ def run_mhe(
     #
     # Load control input data for simulation
     #
-    with open("control_input_data.json", "r") as fr:
+    directory = os.path.abspath(os.path.dirname(__file__))
+    filepath = os.path.join(directory, "control_input_data.json")
+    with open(filepath, "r") as fr:
         control_data = json.load(fr)
 
     control_input_time = [i*sample_period for i in range(len(control_data))]

From 16ba7452eeb2a9e8057e03379a1c58a12c9f02cf Mon Sep 17 00:00:00 2001
From: kuanhanl_DH4200 <kuanhanl@andrew.cmu.edu>
Date: Sun, 1 May 2022 15:49:22 -0400
Subject: [PATCH 08/26] add test for pipeline mhe

---
 .../pipeline/tests/test_pipeline_mhe.py       | 57 +++++++++++++++++++
 1 file changed, 57 insertions(+)
 create mode 100644 nmpc_examples/pipeline/tests/test_pipeline_mhe.py

diff --git a/nmpc_examples/pipeline/tests/test_pipeline_mhe.py b/nmpc_examples/pipeline/tests/test_pipeline_mhe.py
new file mode 100644
index 0000000..5366850
--- /dev/null
+++ b/nmpc_examples/pipeline/tests/test_pipeline_mhe.py
@@ -0,0 +1,57 @@
+import pyomo.common.unittest as unittest
+import pyomo.environ as pyo
+
+from nmpc_examples.pipeline.run_pipeline_mhe import run_mhe
+
+class TestPipelineMHE(unittest.TestCase):
+
+    def test_small_mhe_simulation(self):
+        simulation_data, estimate_data = run_mhe(
+            simulation_horizon=4.0,
+            sample_period=2.0,
+        )
+
+        # These values were not reproduced with another model, just taken
+        # from this MHE simulation. They are here to make sure the result
+        # doesn't change.
+        pred_inlet_flow = [
+            3.000e5, 4.927e5, 4.446e5, 4.460e5, 4.557e5, 4.653e5,
+            4.722e5, 4.785e5, 4.842e5,
+        ]
+        pred_outlet_pressure = [
+            50.91, 47.82, 46.19, 45.06, 43.58, 43.08, 42.62, 42.21, 41.83,
+        ]
+
+        pred_estimate_inlet_flow = [3.00000e5, 4.557e5, 4.842e5]
+        pred_estimate_outlet_pressure = [50.91, 43.58, 41.83]
+
+        actual_inlet_flow = simulation_data.get_data()[
+            pyo.ComponentUID("fs.pipeline.control_volume.flow_mass[*,0.0]")
+        ]
+        actual_outlet_pressure = simulation_data.get_data()[
+            pyo.ComponentUID("fs.pipeline.control_volume.pressure[*,1.0]")
+        ]
+
+        estimate_inlet_flow = estimate_data.get_data()[
+            pyo.ComponentUID("fs.pipeline.control_volume.flow_mass[*,0.0]")
+        ]
+        estimate_outlet_pressure = estimate_data.get_data()[
+            pyo.ComponentUID("fs.pipeline.control_volume.pressure[*,1.0]")
+        ]
+
+        self.assertStructuredAlmostEqual(
+            pred_inlet_flow, actual_inlet_flow, reltol=1e-3
+        )
+        self.assertStructuredAlmostEqual(
+            pred_outlet_pressure, actual_outlet_pressure, reltol=1e-3
+        )
+        self.assertStructuredAlmostEqual(
+            pred_estimate_inlet_flow, estimate_inlet_flow, reltol=1e-3
+        )
+        self.assertStructuredAlmostEqual(
+            pred_estimate_outlet_pressure, estimate_outlet_pressure, reltol=1e-3
+        )
+
+
+if __name__ == "__main__":
+    unittest.main()

From 3c4c01c16922c7410dd472517c9c84d77b22031a Mon Sep 17 00:00:00 2001
From: kuanhanl_DH4200 <kuanhanl@andrew.cmu.edu>
Date: Mon, 2 May 2022 11:05:03 -0400
Subject: [PATCH 09/26] minor changes

---
 nmpc_examples/pipeline/run_pipeline_mhe.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/nmpc_examples/pipeline/run_pipeline_mhe.py b/nmpc_examples/pipeline/run_pipeline_mhe.py
index e34252e..ce0d6dd 100644
--- a/nmpc_examples/pipeline/run_pipeline_mhe.py
+++ b/nmpc_examples/pipeline/run_pipeline_mhe.py
@@ -57,7 +57,7 @@ def run_mhe(
     m_plant.fs.pipeline.control_volume.pressure[:, x0].fix()
 
     #
-    # Make steady model for past data
+    # Make steady model for plant's and estimator's initializations
     #
     m_steady = make_model(dynamic=False, nxfe=nxfe)
     m_steady.fs.pipeline.control_volume.flow_mass[:, x0].fix(
@@ -148,7 +148,7 @@ def run_mhe(
 
     cv = m_estimator.fs.pipeline.control_volume
     measured_variables = [
-        pyo.Reference(cv.flow_mass[:, x]) for x in space.ordered_data()[:-1]
+        pyo.Reference(cv.flow_mass[:, x]) for x in list(space)[:-1]
     ]
     measurement_info = construct_measurement_variables_constraints(
         m_estimator.fs.sample_points,
@@ -169,7 +169,7 @@ def run_mhe(
     ]
     flatten_material_balance = [
         pyo.Reference(cv.material_balances[:,x,"Vap","natural_gas"])
-        for x in space if x != space.last()
+        for x in list(space)[:-1]
     ]
     model_constraints_to_be_disturbed = \
         flatten_momentum_balance + flatten_material_balance
@@ -264,7 +264,7 @@ def run_mhe(
         flatten_measurements,
     )
 
-    # set up a model linker to send measurements to estimator to initialize
+    # Set up a model linker to send measurements to estimator to initialize
     # measured variables
     #
     estimate_linker = DynamicVarLinker(
@@ -333,7 +333,7 @@ def run_mhe(
         # Initialize measured variables within the last sample period to
         # current measurements
         #
-        last_sample_period = list(m_estimator.fs.time)[-4:]
+        last_sample_period = list(m_estimator.fs.time)[-ntfe_per_sample:]
         # for index, var in enumerate(measured_variables):
         #     for tp in last_sampel_period:
         #         var[tp].set_value(blo.measurement_variables[index, estimator_tf])

From 6894c668958d2fa4c7c65af593094c79d9bb2873 Mon Sep 17 00:00:00 2001
From: kuanhanl_DH4200 <kuanhanl@andrew.cmu.edu>
Date: Mon, 2 May 2022 11:05:41 -0400
Subject: [PATCH 10/26] fix variable bugs

---
 nmpc_examples/mhe/mhe_constructor.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/nmpc_examples/mhe/mhe_constructor.py b/nmpc_examples/mhe/mhe_constructor.py
index ad06e7e..2289fae 100644
--- a/nmpc_examples/mhe/mhe_constructor.py
+++ b/nmpc_examples/mhe/mhe_constructor.py
@@ -250,12 +250,12 @@ def get_error_disturbance_cost(
     ]
 
     if weight_data is None:
-        weight_data = {cuid: 1.0 for cuid in cuids}
+        weight_data = {cuid: 1.0 for cuid in component_cuids}
     for i, cuid in enumerate(component_cuids):
         if cuid not in weight_data:
             raise KeyError(
                 "Error/disturbance weight dictionary does not contain a key for "
-                "variable\n%s with ComponentUID %s" % (variables[i].name, cuid)
+                "variable\n%s with ComponentUID %s" % (components[i].name, cuid)
             )
 
     def error_disturbance_rule(m, spt):

From 465ae923676078ca6a5cad1a7392af3cfcf91556 Mon Sep 17 00:00:00 2001
From: Robert Parker <robbybparker@gmail.com>
Date: Tue, 31 May 2022 16:12:26 -0400
Subject: [PATCH 11/26] adjust test to use a smaller estimator model and
 therefore run faster

---
 nmpc_examples/pipeline/tests/test_pipeline_mhe.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/nmpc_examples/pipeline/tests/test_pipeline_mhe.py b/nmpc_examples/pipeline/tests/test_pipeline_mhe.py
index 5366850..949c5bd 100644
--- a/nmpc_examples/pipeline/tests/test_pipeline_mhe.py
+++ b/nmpc_examples/pipeline/tests/test_pipeline_mhe.py
@@ -9,6 +9,12 @@ def test_small_mhe_simulation(self):
         simulation_data, estimate_data = run_mhe(
             simulation_horizon=4.0,
             sample_period=2.0,
+            # NOTE: When I reduce the estimator horizon to 4.0, the test
+            # runs much faster, but the results here don't change.
+            # Is this what I expect?
+            # It seems plausible, as the estimated values are the same
+            # as those simulated...
+            estimator_horizon=4.0,
         )
 
         # These values were not reproduced with another model, just taken

From f9b05294e64c358d56481dcf9c57bf0bf665b0d3 Mon Sep 17 00:00:00 2001
From: Robert Parker <robbybparker@gmail.com>
Date: Tue, 31 May 2022 16:16:17 -0400
Subject: [PATCH 12/26] construct control_input_data in a python function
 rather than leaving as json

---
 .../pipeline/control_input_data.json          |  1 -
 nmpc_examples/pipeline/run_pipeline_mhe.py    | 25 ++++++++++++-------
 2 files changed, 16 insertions(+), 10 deletions(-)
 delete mode 100644 nmpc_examples/pipeline/control_input_data.json

diff --git a/nmpc_examples/pipeline/control_input_data.json b/nmpc_examples/pipeline/control_input_data.json
deleted file mode 100644
index 3b13b87..0000000
--- a/nmpc_examples/pipeline/control_input_data.json
+++ /dev/null
@@ -1 +0,0 @@
-{"fs.pipeline.control_volume.flow_mass[*,1.0]": [553003.1393240632, 535472.4292275416], "fs.pipeline.control_volume.pressure[*,0]": [58.905880452751845, 58.92036532907178]}
\ No newline at end of file
diff --git a/nmpc_examples/pipeline/run_pipeline_mhe.py b/nmpc_examples/pipeline/run_pipeline_mhe.py
index ce0d6dd..22b4be3 100644
--- a/nmpc_examples/pipeline/run_pipeline_mhe.py
+++ b/nmpc_examples/pipeline/run_pipeline_mhe.py
@@ -24,6 +24,21 @@
 """
 
 
+def get_control_inputs(sample_period=2.0):
+    n_samples = 2
+    control_input_time = [i*sample_period for i in range(n_samples)]
+    control_input_data = {
+        'fs.pipeline.control_volume.flow_mass[*,1.0]': [
+            553003.1393240632, 535472.4292275416
+        ],
+        'fs.pipeline.control_volume.pressure[*,0]': [
+            58.905880452751845, 58.92036532907178
+        ],
+    }
+    series = TimeSeriesData(control_input_data, control_input_time)
+    return series
+
+
 def run_mhe(
         simulation_horizon=20.0,
         estimator_horizon=20.0,
@@ -275,15 +290,7 @@ def run_mhe(
     #
     # Load control input data for simulation
     #
-    directory = os.path.abspath(os.path.dirname(__file__))
-    filepath = os.path.join(directory, "control_input_data.json")
-    with open(filepath, "r") as fr:
-        control_data = json.load(fr)
-
-    control_input_time = [i*sample_period for i in range(len(control_data))]
-    control_inputs = TimeSeriesData(
-        control_data, control_input_time, time_set=None
-    )
+    control_inputs = get_control_inputs()
 
     for i in range(n_cycles):
         # time.first() in the model corresponds to sim_t0 in "simulation time"

From c379a638ad38b3fe7aca9d03a3459803cfd7adf2 Mon Sep 17 00:00:00 2001
From: Robert Parker <robbybparker@gmail.com>
Date: Tue, 31 May 2022 16:17:53 -0400
Subject: [PATCH 13/26] comment about getting json file if necessary

---
 nmpc_examples/pipeline/run_pipeline_mhe.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/nmpc_examples/pipeline/run_pipeline_mhe.py b/nmpc_examples/pipeline/run_pipeline_mhe.py
index 22b4be3..bf93a99 100644
--- a/nmpc_examples/pipeline/run_pipeline_mhe.py
+++ b/nmpc_examples/pipeline/run_pipeline_mhe.py
@@ -36,6 +36,8 @@ def get_control_inputs(sample_period=2.0):
         ],
     }
     series = TimeSeriesData(control_input_data, control_input_time)
+    # Note that if we want a json representation of this data, we
+    # can always call json.dump(fp, series.to_serializable()).
     return series
 
 

From c80f12212ebc93395d6c408918f6a6a2b310df1b Mon Sep 17 00:00:00 2001
From: Robert Parker <robbybparker@gmail.com>
Date: Tue, 31 May 2022 18:05:50 -0400
Subject: [PATCH 14/26] comment on disturbance cost function

---
 nmpc_examples/mhe/mhe_constructor.py | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/nmpc_examples/mhe/mhe_constructor.py b/nmpc_examples/mhe/mhe_constructor.py
index 2289fae..5990672 100644
--- a/nmpc_examples/mhe/mhe_constructor.py
+++ b/nmpc_examples/mhe/mhe_constructor.py
@@ -244,6 +244,14 @@ def get_error_disturbance_cost(
     errors or disturbances.
 
     """
+    # In the current formulation, we would like a constant setpoint of zero for
+    # the error variables. We could imagine this as a time-varying setpoint
+    # for the measurable variables. The former is easier, however.
+    # How would I like to parameterize this function?
+    # (variables, time (sample points), setpoint, weights)
+    # Literally should just call the existing function.
+    # TODO: Refactor this function so it does that.
+    # Need time set to get the time-indexed CUID. Will this complicate things?
     component_cuids = [
         get_time_indexed_cuid(comp, sets=(time,))
         for comp in components

From 12b399acfb825635158499d6935e3b44c0cbd328 Mon Sep 17 00:00:00 2001
From: Robert Parker <robbybparker@gmail.com>
Date: Tue, 31 May 2022 18:06:39 -0400
Subject: [PATCH 15/26] make sure code work swith use_linker=True

---
 nmpc_examples/pipeline/run_pipeline_mhe.py | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/nmpc_examples/pipeline/run_pipeline_mhe.py b/nmpc_examples/pipeline/run_pipeline_mhe.py
index bf93a99..ac6cf83 100644
--- a/nmpc_examples/pipeline/run_pipeline_mhe.py
+++ b/nmpc_examples/pipeline/run_pipeline_mhe.py
@@ -97,6 +97,11 @@ def run_mhe(
     # Load data into dynamic model
     #
     use_linker = False
+    # Need to create the plant helper regardless of whether I use it for
+    # this particular data transfer, as it is used to initialize and extend
+    # the simulation data structures.
+    m_plant_helper = DynamicModelHelper(m_plant, m_plant.fs.time)
+    m_plant_helper.load_scalar_data(scalar_data)
     if use_linker:
         # If I want to use DynamicVarLinker:
         # (E.g. if we didn't know that names in the steady model would
@@ -117,8 +122,8 @@ def run_mhe(
 
     else:
         # If I want to use DynamicModelHelper:
-        m_plant_helper = DynamicModelHelper(m_plant, m_plant.fs.time)
-        m_plant_helper.load_scalar_data(scalar_data)
+        #m_plant_helper = DynamicModelHelper(m_plant, m_plant.fs.time)
+        #m_plant_helper.load_scalar_data(scalar_data)
         m_plant_helper.load_data_at_time(initial_data)
 
     # Solve as a sanity check -- should be square with zero infeasibility

From f019410c17bdc8083816a14860248d614f31b166 Mon Sep 17 00:00:00 2001
From: Robert Parker <robbybparker@gmail.com>
Date: Tue, 31 May 2022 18:06:53 -0400
Subject: [PATCH 16/26] update log

---
 log/log.tex | 25 +++++++++++++++++++++++++
 1 file changed, 25 insertions(+)

diff --git a/log/log.tex b/log/log.tex
index a5d40f3..e058066 100644
--- a/log/log.tex
+++ b/log/log.tex
@@ -10,6 +10,31 @@
 \begin{document}
 \maketitle
 
+\section{May 31, 2022}
+
+Modifying KH's branch to make sure it makes use of my existing tools whenever
+possible.
+\begin{itemize}
+  \item The test runs quickly. I can now make changes and run it to make sure I
+    didn't horribly break anything.
+  \item Remove json file for control inputs. -- {\bf Done.}
+  \item Make sure code works with both model linker and model helper for data
+    transfer between models -- {\bf Done.} The problem was that I wasn't
+    initializing the scalar variables.
+  \item Now I need to replace \texttt{get\_error\_disturbance\_cost} with
+    the time-varying setpoint tracking cost.
+\end{itemize}
+
+Recall that this all falls into the category of making sure my data structures
+hold up and can be used in KH's example. I need to make sure they can be used
+everywhere I think they can be used.
+These modifications would be sufficient for a PR into KH's branch. Then another
+PR could incorporate whatever new functionality I decide to add.
+Functionality to add:
+\begin{itemize}
+  \item Data structure for scalar data for time-indexed variable.
+\end{itemize}
+
 \section{Mar 17, 2022}
 
 Where did the code in the \texttt{simple\_pipeline} directory go?

From 34d3544fb02e0e127718c6f18178ed6f310dacd8 Mon Sep 17 00:00:00 2001
From: Robert Parker <robbybparker@gmail.com>
Date: Wed, 1 Jun 2022 13:28:38 -0400
Subject: [PATCH 17/26] add small get_error_cost function

---
 nmpc_examples/mhe/mhe_constructor.py | 56 +++++++++++++++++++++++++++-
 1 file changed, 54 insertions(+), 2 deletions(-)

diff --git a/nmpc_examples/mhe/mhe_constructor.py b/nmpc_examples/mhe/mhe_constructor.py
index 5990672..7dbcaa4 100644
--- a/nmpc_examples/mhe/mhe_constructor.py
+++ b/nmpc_examples/mhe/mhe_constructor.py
@@ -252,6 +252,8 @@ def get_error_disturbance_cost(
     # Literally should just call the existing function.
     # TODO: Refactor this function so it does that.
     # Need time set to get the time-indexed CUID. Will this complicate things?
+
+    # time should not be strictly necessary here (we never expect vardata, do we?)
     component_cuids = [
         get_time_indexed_cuid(comp, sets=(time,))
         for comp in components
@@ -266,12 +268,62 @@ def get_error_disturbance_cost(
                 "variable\n%s with ComponentUID %s" % (components[i].name, cuid)
             )
 
+    # We are assuming that model variables and error variables are in the
+    # same order.
     def error_disturbance_rule(m, spt):
         return sum(
             weight_data[cuid] * error_dist_var[index, spt]**2
             for index, cuid in enumerate(component_cuids)
         )
-    error_disturbance_cost = Expression(
-        sample_points, rule=error_disturbance_rule
+    #error_disturbance_cost = Expression(
+    #    sample_points, rule=error_disturbance_rule
+    #)
+    from nmpc_examples.nmpc.cost_expressions import (
+        get_tracking_cost_from_constant_setpoint,
+    )
+    from pyomo.core.base.reference import Reference
+    error_var_refs = [
+        # This is a downside of having error_dist_var be two-dimensional.
+        # It is not easy to get the number of vars it refers to.
+        # Note that we are using len(components) here...
+        Reference(error_dist_var[idx, :]) for idx in range(len(components))
+    ]
+    error_var_cuids = [get_time_indexed_cuid(var) for var in error_var_refs]
+    setpoint_data = {cuid: 0.0 for cuid in error_var_cuids}
+    #weight_data = {cuid: 1.0 for cuid in error_var_cuids}
+    # Need to get weight_data from the function argument.
+    # Map component cuid -> error var cuid somehow...
+    # Component cuid -> component index (= error var index) -> error var cuid
+    cuid_index_map = {cuid: idx for idx, cuid in enumerate(component_cuids)}
+    weight_data = {
+        error_var_cuids[cuid_index_map[cuid]]: val
+        for cuid, val in weight_data.items()
+    }
+    error_disturbance_cost = get_tracking_cost_from_constant_setpoint(
+        error_var_refs,
+        sample_points,
+        setpoint_data,
+        weight_data=weight_data,
     )
     return error_disturbance_cost
+
+
+def get_error_cost(variables, time, weight_data=None):
+    """
+    """
+    # I know that these variables are indexed by time (aren't VarData), so I
+    # don't need to send the time set to get_time_indexed_cuid
+    cuids = [get_time_indexed_cuid(var) for var in variables]
+    if weight_data is None:
+        weight_data = {cuid: 1.0 for cuid in cuids}
+    setpoint_data = {cuid: 0.0 for cuid in cuids}
+    from nmpc_examples.nmpc.cost_expressions import (
+        get_tracking_cost_from_constant_setpoint,
+    )
+    error_cost = get_tracking_cost_from_constant_setpoint(
+        variables,
+        time,
+        setpoint_data,
+        weight_data=weight_data,
+    )
+    return error_cost

From 6ae487218a53a397b234182d89ccd746ccfe17d8 Mon Sep 17 00:00:00 2001
From: Robert Parker <robbybparker@gmail.com>
Date: Wed, 1 Jun 2022 13:29:03 -0400
Subject: [PATCH 18/26] remove get_error_disturbance_cost

---
 nmpc_examples/mhe/mhe_constructor.py | 101 ---------------------------
 1 file changed, 101 deletions(-)

diff --git a/nmpc_examples/mhe/mhe_constructor.py b/nmpc_examples/mhe/mhe_constructor.py
index 7dbcaa4..11856cd 100644
--- a/nmpc_examples/mhe/mhe_constructor.py
+++ b/nmpc_examples/mhe/mhe_constructor.py
@@ -207,107 +207,6 @@ def activate_disturbed_constraints_based_on_original_constraints(
                 disturbance_var[index, spt].fix(0.0)
 
 
-def get_error_disturbance_cost(
-    time,
-    sample_points,
-    components,
-    error_dist_var,
-    weight_data=None,
-):
-    """
-    #TODO: This function is similar to "get_tracking_cost_from_constant_setpoint",
-    #Should I merge this one to that one?
-    This function return a squared cost expression for measurement errors or
-    model disturbances.
-
-    Arguments
-    ---------
-    time: iterable
-        Time set which indexes components (measured variables or
-        model constraints)
-    sample_points: iterable
-        Set of sample points
-    components: List
-        List of components (measured variables or model constraints)
-    error_dist_var: Pyomo Var
-        Variable of measurement error
-        from "construct_measurement_variables_constraints"
-        or model disturbance
-        from "construct_disturbed_model_constraints"
-    weight_data: dict
-        Optional. Maps variable names to squared cost weights. If not provided,
-        weights of one are used.
-
-    Returns
-    -------
-    Pyomo Expression, indexed by sample_points, containing the sum of squared
-    errors or disturbances.
-
-    """
-    # In the current formulation, we would like a constant setpoint of zero for
-    # the error variables. We could imagine this as a time-varying setpoint
-    # for the measurable variables. The former is easier, however.
-    # How would I like to parameterize this function?
-    # (variables, time (sample points), setpoint, weights)
-    # Literally should just call the existing function.
-    # TODO: Refactor this function so it does that.
-    # Need time set to get the time-indexed CUID. Will this complicate things?
-
-    # time should not be strictly necessary here (we never expect vardata, do we?)
-    component_cuids = [
-        get_time_indexed_cuid(comp, sets=(time,))
-        for comp in components
-    ]
-
-    if weight_data is None:
-        weight_data = {cuid: 1.0 for cuid in component_cuids}
-    for i, cuid in enumerate(component_cuids):
-        if cuid not in weight_data:
-            raise KeyError(
-                "Error/disturbance weight dictionary does not contain a key for "
-                "variable\n%s with ComponentUID %s" % (components[i].name, cuid)
-            )
-
-    # We are assuming that model variables and error variables are in the
-    # same order.
-    def error_disturbance_rule(m, spt):
-        return sum(
-            weight_data[cuid] * error_dist_var[index, spt]**2
-            for index, cuid in enumerate(component_cuids)
-        )
-    #error_disturbance_cost = Expression(
-    #    sample_points, rule=error_disturbance_rule
-    #)
-    from nmpc_examples.nmpc.cost_expressions import (
-        get_tracking_cost_from_constant_setpoint,
-    )
-    from pyomo.core.base.reference import Reference
-    error_var_refs = [
-        # This is a downside of having error_dist_var be two-dimensional.
-        # It is not easy to get the number of vars it refers to.
-        # Note that we are using len(components) here...
-        Reference(error_dist_var[idx, :]) for idx in range(len(components))
-    ]
-    error_var_cuids = [get_time_indexed_cuid(var) for var in error_var_refs]
-    setpoint_data = {cuid: 0.0 for cuid in error_var_cuids}
-    #weight_data = {cuid: 1.0 for cuid in error_var_cuids}
-    # Need to get weight_data from the function argument.
-    # Map component cuid -> error var cuid somehow...
-    # Component cuid -> component index (= error var index) -> error var cuid
-    cuid_index_map = {cuid: idx for idx, cuid in enumerate(component_cuids)}
-    weight_data = {
-        error_var_cuids[cuid_index_map[cuid]]: val
-        for cuid, val in weight_data.items()
-    }
-    error_disturbance_cost = get_tracking_cost_from_constant_setpoint(
-        error_var_refs,
-        sample_points,
-        setpoint_data,
-        weight_data=weight_data,
-    )
-    return error_disturbance_cost
-
-
 def get_error_cost(variables, time, weight_data=None):
     """
     """

From f993c77b89bdff0815eb81cf073368dfee4746aa Mon Sep 17 00:00:00 2001
From: Robert Parker <robbybparker@gmail.com>
Date: Wed, 1 Jun 2022 14:35:19 -0400
Subject: [PATCH 19/26] update to use new get_error_cost function

---
 nmpc_examples/pipeline/run_pipeline_mhe.py | 86 ++++++++++++++++------
 1 file changed, 62 insertions(+), 24 deletions(-)

diff --git a/nmpc_examples/pipeline/run_pipeline_mhe.py b/nmpc_examples/pipeline/run_pipeline_mhe.py
index ac6cf83..e17545e 100644
--- a/nmpc_examples/pipeline/run_pipeline_mhe.py
+++ b/nmpc_examples/pipeline/run_pipeline_mhe.py
@@ -1,15 +1,19 @@
 import pyomo.environ as pyo
 from pyomo.dae import ContinuousSet
 from pyomo.dae.flatten import flatten_dae_components
+from pyomo.common.collections import ComponentMap
 
 from nmpc_examples.pipeline.pipeline_model import make_model
 
 from nmpc_examples.mhe.mhe_constructor import (
     construct_measurement_variables_constraints,
     construct_disturbed_model_constraints,
-    get_error_disturbance_cost,
+    get_error_cost,
     activate_disturbed_constraints_based_on_original_constraints,
 )
+from nmpc_examples.nmpc.cost_expressions import (
+    get_tracking_cost_from_time_varying_setpoint,
+)
 
 from nmpc_examples.nmpc.model_linker import DynamicVarLinker
 from nmpc_examples.nmpc.model_helper import DynamicModelHelper
@@ -222,31 +226,65 @@ def run_mhe(
     # Construct least square objective to minimize measurement errors
     # and model disturbances
     #
-    measurement_error_weights = {
-        pyo.ComponentUID(var.referent): 1.0
-        for var in measured_variables
-    }
-    m_estimator.measurement_error_cost = get_error_disturbance_cost(
-        m_estimator.fs.time,
-        m_estimator.fs.sample_points,
-        measured_variables,
-        esti_blo.measurement_error_variables,
-        measurement_error_weights,
-    )
+    # This flag toggles between two different objective formulations.
+    # I included it just to test that we can support both.
+    error_var_objective = True
+    if measurement_var_objective:
+        error_vars = [
+            pyo.Reference(esti_blo.measurement_error_variables[idx, :])
+            for idx in esti_blo.measurement_set
+        ]
+        # This cost function penalizes the square of the "error variables"
+        m_estimator.measurement_error_cost = get_error_cost(
+            error_vars, m_estimator.fs.sample_points
+        )
+    else:
+        measurement_map = ComponentMap(
+            (var, [
+                esti_blo.measurement_variables[i, t]
+                for t in m_estimator.fs.sample_points
+            ])
+            for i, var in enumerate(measured_variables)
+        )
+        setpoint_data = TimeSeriesData(
+            measurement_map, m_estimator.fs.sample_points
+        )
+        # This cost function penalizes the difference between measurable
+        # estimates and their corresponding measurements.
+        error_cost = get_tracking_cost_from_time_varying_setpoint(
+            measured_variables, m_estimator.fs.sample_points, setpoint_data
+        )
+        m_estimator.measurement_error_cost = error_cost
 
-    model_disturbance_weights = {
-        **{pyo.ComponentUID(con.referent): 10.0
-           for con in flatten_momentum_balance},
-        **{pyo.ComponentUID(con.referent): 20.0
-           for con in flatten_material_balance},
-    }
-    m_estimator.model_disturbance_cost = get_error_disturbance_cost(
-        m_estimator.fs.time,
-        m_estimator.fs.sample_points,
-        model_constraints_to_be_disturbed,
-        esti_blo.disturbance_variables,
-        model_disturbance_weights,
+    #
+    # Construct disturbance cost expression
+    #
+    disturbance_vars = [
+        pyo.Reference(esti_blo.disturbance_variables[idx, :])
+        for idx in esti_blo.disturbance_set
+    ]
+
+    # We know what order we sent constraints to the disturbance constraint
+    # function, so we know which indices correspond to which balance equations
+    momentum_bal_indices = set(range(len(flatten_momentum_balance)))
+    material_bal_indices = set(range(
+        len(flatten_momentum_balance),
+        len(flatten_momentum_balance) + len(flatten_material_balance),
+    ))
+    weights = {}
+    for idx in esti_blo.disturbance_set:
+        slc = disturbance_vars[idx].referent
+        if idx in momentum_bal_indices:
+            weight = 10.0
+        elif idx in material_bal_indices:
+            weight = 20.0
+        else:
+            weight = 1.0
+        weights[pyo.ComponentUID(slc)] = weight
+    m_estimator.model_disturbance_cost = get_error_cost(
+        disturbance_vars, m_estimator.fs.sample_points, weight_data=weights
     )
+    ###
 
     m_estimator.squred_error_disturbance_objective = pyo.Objective(
         expr=(sum(m_estimator.measurement_error_cost.values()) +

From 69c0e27b6df39fe06cdc140f07d6fc1af64d16c8 Mon Sep 17 00:00:00 2001
From: Robert Parker <robbybparker@gmail.com>
Date: Wed, 1 Jun 2022 14:36:24 -0400
Subject: [PATCH 20/26] update log

---
 log/log.tex | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/log/log.tex b/log/log.tex
index e058066..000c6f0 100644
--- a/log/log.tex
+++ b/log/log.tex
@@ -10,6 +10,20 @@
 \begin{document}
 \maketitle
 
+\section{June 1, 2022}
+\begin{itemize}
+  \item I can use a small, simple \texttt{get\_error\_cost} function in place
+    of KH's \texttt{get\_error\_disturbance\_cost}. Now I should delete that
+    old function.
+  \item It seems like my setpoint tracking cost function holds up to this
+    challenge.
+  \item What if we wanted a time-varying setpoint between measurements
+    and estimated measurable variables? I should probably test this?
+    -- {\bf Done.} Using a time-varying tracking cost between measurements
+    and estimated measurable variables seems to work.
+  \item Regardless, \texttt{get\_error\_disturbance\_cost} can go.
+\end{itemize}
+
 \section{May 31, 2022}
 
 Modifying KH's branch to make sure it makes use of my existing tools whenever

From 4463b314e28d44b1faeeb124e8895ecabdf768a1 Mon Sep 17 00:00:00 2001
From: Robert Parker <robbybparker@gmail.com>
Date: Mon, 25 Jul 2022 17:21:36 -0400
Subject: [PATCH 21/26] fix name of flag

---
 nmpc_examples/pipeline/run_pipeline_mhe.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/nmpc_examples/pipeline/run_pipeline_mhe.py b/nmpc_examples/pipeline/run_pipeline_mhe.py
index e17545e..778110e 100644
--- a/nmpc_examples/pipeline/run_pipeline_mhe.py
+++ b/nmpc_examples/pipeline/run_pipeline_mhe.py
@@ -229,7 +229,7 @@ def run_mhe(
     # This flag toggles between two different objective formulations.
     # I included it just to test that we can support both.
     error_var_objective = True
-    if measurement_var_objective:
+    if error_var_objective:
         error_vars = [
             pyo.Reference(esti_blo.measurement_error_variables[idx, :])
             for idx in esti_blo.measurement_set

From 2f7c36b026265222e22572ae08e074f314f6a29d Mon Sep 17 00:00:00 2001
From: Robert Parker <robbybparker@gmail.com>
Date: Mon, 25 Jul 2022 17:53:44 -0400
Subject: [PATCH 22/26] remove comments and simplify slightly

---
 nmpc_examples/pipeline/run_pipeline_mhe.py | 8 +-------
 1 file changed, 1 insertion(+), 7 deletions(-)

diff --git a/nmpc_examples/pipeline/run_pipeline_mhe.py b/nmpc_examples/pipeline/run_pipeline_mhe.py
index 778110e..fc47ffc 100644
--- a/nmpc_examples/pipeline/run_pipeline_mhe.py
+++ b/nmpc_examples/pipeline/run_pipeline_mhe.py
@@ -101,9 +101,6 @@ def run_mhe(
     # Load data into dynamic model
     #
     use_linker = False
-    # Need to create the plant helper regardless of whether I use it for
-    # this particular data transfer, as it is used to initialize and extend
-    # the simulation data structures.
     m_plant_helper = DynamicModelHelper(m_plant, m_plant.fs.time)
     m_plant_helper.load_scalar_data(scalar_data)
     if use_linker:
@@ -126,8 +123,6 @@ def run_mhe(
 
     else:
         # If I want to use DynamicModelHelper:
-        #m_plant_helper = DynamicModelHelper(m_plant, m_plant.fs.time)
-        #m_plant_helper.load_scalar_data(scalar_data)
         m_plant_helper.load_data_at_time(initial_data)
 
     # Solve as a sanity check -- should be square with zero infeasibility
@@ -273,14 +268,13 @@ def run_mhe(
     ))
     weights = {}
     for idx in esti_blo.disturbance_set:
-        slc = disturbance_vars[idx].referent
         if idx in momentum_bal_indices:
             weight = 10.0
         elif idx in material_bal_indices:
             weight = 20.0
         else:
             weight = 1.0
-        weights[pyo.ComponentUID(slc)] = weight
+        weights[esti_blo.disturbance_variables[idx, :]] = weight
     m_estimator.model_disturbance_cost = get_error_cost(
         disturbance_vars, m_estimator.fs.sample_points, weight_data=weights
     )

From a689f664f705350d35f378d3eb6ed6a97eb81b9c Mon Sep 17 00:00:00 2001
From: Robert Parker <robbybparker@gmail.com>
Date: Mon, 25 Jul 2022 17:54:14 -0400
Subject: [PATCH 23/26] remove old commented code

---
 nmpc_examples/pipeline/run_pipeline_mhe.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/nmpc_examples/pipeline/run_pipeline_mhe.py b/nmpc_examples/pipeline/run_pipeline_mhe.py
index fc47ffc..dc74416 100644
--- a/nmpc_examples/pipeline/run_pipeline_mhe.py
+++ b/nmpc_examples/pipeline/run_pipeline_mhe.py
@@ -380,9 +380,6 @@ def run_mhe(
         # current measurements
         #
         last_sample_period = list(m_estimator.fs.time)[-ntfe_per_sample:]
-        # for index, var in enumerate(measured_variables):
-        #     for tp in last_sampel_period:
-        #         var[tp].set_value(blo.measurement_variables[index, estimator_tf])
         estimate_linker.transfer(tf, last_sample_period)
 
         # Load inputs to estimator

From 002491478bd9f09de1ced9ebf21e61c521b6e89b Mon Sep 17 00:00:00 2001
From: Robert Parker <robbybparker@gmail.com>
Date: Mon, 25 Jul 2022 17:54:41 -0400
Subject: [PATCH 24/26] comment on confusing name

---
 nmpc_examples/mhe/mhe_constructor.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/nmpc_examples/mhe/mhe_constructor.py b/nmpc_examples/mhe/mhe_constructor.py
index 11856cd..c345c19 100644
--- a/nmpc_examples/mhe/mhe_constructor.py
+++ b/nmpc_examples/mhe/mhe_constructor.py
@@ -212,6 +212,10 @@ def get_error_cost(variables, time, weight_data=None):
     """
     # I know that these variables are indexed by time (aren't VarData), so I
     # don't need to send the time set to get_time_indexed_cuid
+    #
+    # NOTE: This name is slightly confusing. Something better might be
+    # get_cost_from_error_variables
+    # error_cost seems to imply an "error" between two quantities...
     cuids = [get_time_indexed_cuid(var) for var in variables]
     if weight_data is None:
         weight_data = {cuid: 1.0 for cuid in cuids}

From 4ae5ab77e058d64a3174d0860517af16ef6e7c3d Mon Sep 17 00:00:00 2001
From: Robert Parker <robbybparker@gmail.com>
Date: Mon, 25 Jul 2022 17:55:51 -0400
Subject: [PATCH 25/26] update log

---
 log/log.tex | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/log/log.tex b/log/log.tex
index eb21970..9e96779 100644
--- a/log/log.tex
+++ b/log/log.tex
@@ -53,6 +53,13 @@ \subsection{Testing this functionality}
 
 \begin{itemize}
   \item NMPC and open-loop have been tested.
+  \item KH's MHE code.
+    \begin{itemize}
+      \item Tests pass.
+      \item Do the new data structure buy me anything here?
+      \item They make the code for sending weights to the disturbance objective
+        slightly easier to understand.
+    \end{itemize}
 \end{itemize}
 
 \section{June 9, 2022}

From 76825face9e8d205c9de0c3b9c4c2037616b3a5d Mon Sep 17 00:00:00 2001
From: Robert Parker <robbybparker@gmail.com>
Date: Wed, 31 Aug 2022 16:13:58 -0400
Subject: [PATCH 26/26] remove additions from git merge

---
 log/log.tex | 6 ------
 1 file changed, 6 deletions(-)

diff --git a/log/log.tex b/log/log.tex
index 487ed07..759a3de 100644
--- a/log/log.tex
+++ b/log/log.tex
@@ -110,7 +110,6 @@ \subsection{Testing this functionality}
 
 \begin{itemize}
   \item NMPC and open-loop have been tested.
-<<<<<<< HEAD
   \item KH's MHE code.
     \begin{itemize}
       \item Tests pass.
@@ -118,10 +117,8 @@ \subsection{Testing this functionality}
       \item They make the code for sending weights to the disturbance objective
         slightly easier to understand.
     \end{itemize}
-=======
   \item KH's MHE has been tested.
   \item Now I need to test the code in \texttt{gas\_distribution}
->>>>>>> mhe-patch
 \end{itemize}
 
 \section{June 9, 2022}
@@ -160,8 +157,6 @@ \section{June 9, 2022}
   \item In scripts for pipeline NMPC
 \end{itemize}
 
-<<<<<<< HEAD
-=======
 \section{June 1, 2022}
 
 Kuan-Han's MHE example appears to work fine using my data structures
@@ -175,7 +170,6 @@ \section{June 1, 2022}
   \item Maybe some functions for adding noise to data.
 \end{itemize}
 
->>>>>>> mhe-patch
 \section{June 1, 2022}
 \begin{itemize}
   \item I can use a small, simple \texttt{get\_error\_cost} function in place