How to set model/emulation initial state in MPC control loop?

[Tom Stesco]

Hi, when using the results of a model optimization as the controller input to the system emulation MPC typically uses the first input step then measures the system and reoptimizes the model from this state. I made a simple script (below, warning - it takes over 10 minutes on my laptop) to do this and plot the trajectories with dynamic constraints computed at each time step, but I am having trouble setting the initial states of the model, optimization, and emulation.

I must be missing something. I tried setting the model measurements to the emulation measurements, but still do not know how to correctly specify the initial state for the model and emulation. As result the initial state appears to be fixed for every MPC step which throws off the MPC emulation.

The optimization state output (plot below) for example starts at Tzone=25 degC but drops extremely quickly. The sharp input changes just after initialization in the optimizer also seem problematic, I'm not yet sure if this is due to the initialization problem or the modelica RC model itself.

The model state (simulated after optimization with cntrol input) starts at Tzone = 20 degC. My thinking is this mismatch in initialization of the optimzer causes the Tzone to dip below 20 degC just after initialization.

 

Of course because the state updating is not working the next MPC step has the same problem and the 1-step emulation itself only captures these initial state issues.

Please let me know if there is a good way to set the initial state and then update it after each MPC step with the emulation measurements. Also, if you have any tips on implementation performance they would be most welcome. Thank you!

# -*- coding: utf-8 -*-
#!/usr/bin/env python
# author: Thomas Stesco
# date: 18.05.2018
'''
This is a python 2.7 script for running an MPC control loop on a simplified
RC model of a building using MPCPy and plotting the input/output trajectories. 
'''
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

from mpcpy import variables
from mpcpy import units
from mpcpy import exodata
from mpcpy import models
from mpcpy import systems
from mpcpy import optimization

import pandas as pd
import numpy as np
import datetime
import os

# MPC params config
model_name = "MPC_loop_test"
horizon_steps = 72
step_length_seconds = 3600.

MPC_start_time_str = '4/1/2017 00:00:00' # using local time
MPC_final_time_str = '4/2/2017 00:00:00' # using local time

# output path config
output_dir = os.path.join("/root", "src", "output", model_name)
if not os.path.exists(output_dir):
    os.makedirs(output_dir)

plot_file_type = ".png"
constraint_plot_name = os.path.join(
    output_dir, "constraint_{}_{}"+plot_file_type
)
optimum_trajectory_plot_name = os.path.join(
    output_dir, "optimum_trajectory_{}" + plot_file_type
)
model_trajectory_plot_name = os.path.join(
    output_dir, "model_trajectory_{}" + plot_file_type
)
emulator_trajectory_plot_name = os.path.join(
    output_dir, "emulator_trajectory_{}" + plot_file_type
)

# initialization
script_start_time = datetime.datetime.now() # for runtime stats
horizon_length_seconds = horizon_steps*step_length_seconds
MPC_start_time = datetime.datetime.strptime(MPC_start_time_str, '%m/%d/%Y %H:%M:%S')
MPC_final_time = datetime.datetime.strptime(MPC_final_time_str, '%m/%d/%Y %H:%M:%S')
MPC_duration_seconds = (MPC_final_time-MPC_start_time).total_seconds()

# init step time variables
step_start_time = MPC_start_time
step_start_time_str = MPC_start_time_str
step_final_time = MPC_start_time + datetime.timedelta(seconds=step_length_seconds)
step_final_time_str = step_final_time.strftime('%m/%d/%Y %H:%M:%S')
step_horizon_final_time = MPC_start_time + datetime.timedelta(seconds=horizon_length_seconds)
step_horizon_final_time_str = step_horizon_final_time.strftime('%m/%d/%Y %H:%M:%S')

# MPC config output
print('start_time in local: {}'.format(MPC_start_time_str))
print('final_time in local: {}'.format(MPC_final_time_str))
print('horizon_length_seconds: {}'.format(horizon_length_seconds))
print('step_length_seconds: {}'.format(step_length_seconds))
print('horizon_steps: {}'.format(horizon_steps))
print('MPC duration (seconds): {}'.format(MPC_duration_seconds))
print('MPC total steps: {}'.format(MPC_duration_seconds/step_length_seconds))

weather = exodata.WeatherFromEPW(
    '/opt/MPCPy/doc/userGuide/tutorial/USA_IL_Chicago-OHare.Intl.AP.725300_TMY3.epw'
)
# want to have full year weather DF object for running mean calculations
weather.collect_data('1/1/2017', '31/12/2017')
weather_df = weather.display_data()

# construct control dataframe
control_init_dict = {
    "Time" : pd.date_range(start=step_start_time_str, periods=horizon_steps, freq='H'),
    "input_Q_flow" : np.ones(horizon_steps)*1000,
}

# make control df
control_init_df = pd.DataFrame.from_dict(data=control_init_dict)
control_init_df = control_init_df.set_index('Time')

# init control
control_variable_map = {'input_Q_flow' : ('input_Q_flow', units.W)}
control_init = exodata.ControlFromDF(
    df=control_init_df,
    variable_map=control_variable_map,
    tz_name=weather.tz_name
)
control_init.collect_data(step_start_time_str, step_final_time_str)

# set measurements
measurements = {
    'Tzone' : {},
    'input_Q_flow' : {},
    'weaTDryBul' : {},
    'heatCapacitor.T' : {},
}
measurements['Tzone']['Sample'] = variables.Static('sample_rate', 
                                                   3600,
                                                   units.s)
measurements['input_Q_flow']['Sample'] = variables.Static('sample_rate',
                                                   3600,
                                                   units.s)
measurements['weaTDryBul']['Sample'] = variables.Static('sample_rate',
                                                   3600,
                                                   units.s)
measurements['heatCapacitor.T']['Sample'] = variables.Static('sample_rate',
                                                   3600,
                                                   units.s)

# model setup
moinfo = ('/root/modelica/SimpleBuilding.mo', 
          'SimpleBuilding.RC', 
          {})

# init emulator
emulation = systems.EmulationFromFMU(
    measurements = measurements,
    moinfo = moinfo,
    weather_data = weather.data,
    control_data = control_init.data,
    tz_name = weather.tz_name
)
emulation_meas = pd.DataFrame()

# init model
model = models.Modelica(
    models.JModelica,
    models.RMSE,
    measurements = measurements,
    moinfo = moinfo,
    parameter_data = {},
    weather_data = weather.data,
    control_data = control_init.data,
    tz_name = weather.tz_name
)

# run MPC loop 
while (step_start_time < MPC_final_time):
    print('='*70)
    print("step start time: {}".format(step_start_time))
    print("step final time: {}".format(step_final_time))
    print("step horizon final time: {}".format(step_horizon_final_time))

    # get_updated_disturbances
    weather.collect_data(step_start_time_str, step_horizon_final_time_str)
    # get updated constraints
    constraint_variable_map = {
        'Qflow_min' : ('input_Q_flow', 'GTE', units.W),
        'Qflow_max' : ('input_Q_flow', 'LTE', units.W),
        'T_min' : ('Tzone', 'GTE', units.degC),
        'T_max' : ('Tzone', 'LTE', units.degC)
    }
    # construct constraint dataframe
    start_time_UTC = pd.to_datetime(step_start_time_str).tz_localize(weather.tz_name)
    final_time_UTC = pd.to_datetime(step_horizon_final_time_str).tz_localize(weather.tz_name)
    constraint_base = {
        "Time" : pd.date_range(start=step_start_time_str, periods=horizon_steps+1, freq='H'),
        "Qflow_min" : np.ones(horizon_steps+1)*0,
        "Qflow_max" : np.ones(horizon_steps+1)*4000
    }
    # make constraints object
    T_min = np.ones(horizon_steps+1)*20
    T_max = np.ones(horizon_steps+1)*25
    constraint_T = {
        'T_min': T_min,
        'T_max': T_max
    }

    # concatenate constraint dicts
    constraint_data = constraint_T.copy()  
    constraint_data.update(constraint_base)
    # make constraint_df
    constraint_df = pd.DataFrame.from_dict(data=constraint_data)
    constraint_df = constraint_df.set_index('Time')
    constraints = exodata.ConstraintFromDF(
        df=constraint_df, 
        variable_map=constraint_variable_map, 
        tz_name=weather.tz_name
    )
    constraints.collect_data(step_start_time_str, step_horizon_final_time_str)

    # T constraint plot
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(T_max, label='T_max', color="red")
    ax.plot(T_min, label='T_min', color="blue")
    ax.set_title('T')
    ax.set_ylabel('Temperature (degC)')
    ax.set_xlabel('Horizon steps')
    ax.legend()
    plt.savefig(constraint_plot_name.format("T", 
        (step_final_time-MPC_start_time).total_seconds()
    ))
    # optimization problem formulation
    opt_problem = optimization.Optimization(
        model, 
        optimization.EnergyMin,
        optimization.JModelica,
        'input_Q_flow',
        constraint_data = constraints.data,
        tz_name=weather.tz_name
    )
    opt_problem.optimize(
        step_start_time_str, 
        step_horizon_final_time_str, 
        res_control_step=1.0
    )

    # trajectory plot
    opt_meas = opt_problem.display_measurements('Simulated')
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(opt_meas['Tzone']-273.15, label='Tzone', color="green")
    ax.plot(weather.display_data()['weaTDryBul'], label='weaTDryBul', color="blue")
    ax2 = ax.twinx()
    ax2.plot(opt_meas['input_Q_flow'], label='input_Q_flow', color="red")
    ax.set_title('Optimizer Trajectory')
    ax.set_ylabel('Temperature (degC)')
    ax2.set_ylabel('Input Heat (W)')
    ax.set_xlabel('Time')
    fig.legend()
    plt.savefig(optimum_trajectory_plot_name.format(
        (step_final_time-MPC_start_time).total_seconds()
    ))

    # NOTE: model/emulation.control_data is updated by the optimization 
    # results to be the new optimum trajectory 
    model.simulate(step_start_time_str, step_horizon_final_time_str)

    # model plot
    model_meas = model.display_measurements('Simulated')
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(model_meas['Tzone']-273.15, label='Tzone', color="green")
    ax.plot(weather.display_data()['weaTDryBul'], label='weaTDryBul', color="blue")
    ax2 = ax.twinx()
    ax2.plot(model_meas['input_Q_flow'], label='input_Q_flow', color="red")
    ax.set_title('Simulated model Trajectory')
    ax.set_ylabel('Temperature (degC)')
    ax2.set_ylabel('Input Heat (W)')
    ax.set_xlabel('Time')
    fig.legend()
    plt.savefig(model_trajectory_plot_name.format(
        (step_final_time-MPC_start_time).total_seconds()
    ))
    emulation.collect_measurements(step_start_time_str, step_final_time_str)
    model.measurements = emulation.measurements
    # accumulate emulator state measurements for single "actual" step
    emulation_meas = emulation_meas.append(emulation.display_measurements('Simulated'))

    # increment step times
    current_iteration = 0
    step_start_time = step_start_time + datetime.timedelta(seconds=step_length_seconds)
    step_final_time = step_final_time + datetime.timedelta(seconds=step_length_seconds)
    step_horizon_final_time = step_start_time \
        + datetime.timedelta(seconds=horizon_length_seconds)
    step_start_time_str = step_start_time.strftime('%m/%d/%Y %H:%M:%S')
    step_final_time_str = step_final_time.strftime('%m/%d/%Y %H:%M:%S')
    step_horizon_final_time = step_horizon_final_time.strftime('%m/%d/%Y %H:%M:%S')

# plot emulation results
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(emulation_meas['Tzone']-273.15, label='Tzone', color="green")
ax.plot(weather.display_data()['weaTDryBul'], label='weaTDryBul', color="blue")
ax2 = ax.twinx()
ax2.plot(emulation_meas['input_Q_flow'], label='input_Q_flow', color="red")
ax.set_title('Emulator Trajectory')
ax.set_ylabel('Temperature (degC)')
ax2.set_ylabel('Input Heat (W)')
ax.set_xlabel('Time')
fig.legend()
plt.savefig(emulator_trajectory_plot_name.format(
    (step_final_time-MPC_start_time).total_seconds()
))

print("Completed")
print("output in {}".format(output_dir))
print("Run time: {} seconds".format(datetime.datetime.now() - script_start_time))

​

[dhb...@gmail.com]

Hey Tom, 

Right now there is no formal construct to initialize the model for simulation and optimization.  What I recommend is to create parameters in the Modelica model to set initialization of certain variables.  Then, in MPCPy, you can use parameter_data to update the parameters before each optimization or model simulation.  It is not in the tutorial example, but you can see an example usage in the unittest "unittests/test_optimization.OptimizeSimpleFromJModelica.test_set_parameters".  So, for example, to initialize the temperature of the capacitor of an RC, try in the Modelica model:

// Define a parameter
parameter Modelica.SIunits.Temperature T0 = 20+273.15 "Initial temperature of C";

// Add modifier to heat capacitor for initialization

Modelica.Thermal.Heat Transfer.Components.HeatCapacitor cap(T(start=T0, fixed=true)) "Heat capacitor with initial temperature T0";

For emulation, I recommend using the 'continue' keyword as the start time after the initial emulation step has been taken.  Using the 'continue' keyword will continue the emulation from where it last left off (e.g. does not reset the fmu) until the final time that specify.  You should still set a start time for the first iteration, however.  It is not in the tutorial example, but you can see an example usage in the unittest "unittests/test_systems.EmulationFromFMU.test_collect_measurements_continue".  Depending on the last time you pulled the latest master branch, this will also speed up your loop a bit, because you will not be reloading the emulation fmu each iteration, which can take a little time.  Although a recent commit to the master fixed this bug so that even if you weren't using the 'continue' keyword, the fmu is not reloaded each time, but instead reset.  Either way, I recommend pulling the latest master.

Let me know if this helps you move forward with your test.

Best Regards, 

Dave

[Tom Stesco]

Hi Dave,

Thanks again for the help and pointing me to the specific unittests. I'll use the parameter_data object for now and try using the 'continue' feature with the emulation. I'll post the results here afterwards.

Best,

Tom

[Tom Stesco]

Hi, to follow up:

Using the parameter_data and ‘continue’ argument with the emulation gives the expected results. There were however issues with feasibility of the optimal control input, which tended to dip below the T_min constraint. I fixed this by changing the OCP option number of communication points ('ncp' kwarg) from the default of 3 to 6 (still trying to figure out exactly why this helps). The model trajectory now matches the optimal trajectory, similarly for the emulation trajectory (though this appears flat due to the shorter emulation duration).

The optimizer trajectory has a strange response on the initial iteration from the boundary state condition, but this does not appear in the later MPC iterations.

Also, I found that if I didn’t re-initialize the OCP for each iteration it would not use the updated constraint data (which I made to be time variant). In case anyone is interested here are excerpts of the code I used based on the test code in “unittests/test_optimization.OptimizeSimpleFromJModelica.test_set_parameters” and “unittests/test_systems.EmulationFromFMU.test_collect_measurements_continue” for the OCP reinitialization and parameter updates within the MPC loop.

    <code cropped>
    # reinitialize OCP with updated constrains

    opt_problem = optimization.Optimization(
        model, 
        optimization.EnergyMin,
        optimization.JModelica,
        'input_Q_flow'

,
        constraint_data = constraints.data,
        tz_name=weather.tz_name
    )
    # change default ncp
    opt_options = opt_problem.get_optimization_options()
    opt_options['n_cp'] = 6
    opt_problem.set_optimization_options(opt_options)
    # solve OCP
    opt_problem.optimize(
        step_start_time_str, 
        step_horizon_final_time_str, 
        res_control_step=1.0
    )

    <code cropped>
    if first_step:
        emulation.collect_measurements(step_start_time_str, step_final_time_str);
    else:
        emulation.collect_measurements('continue', step_final_time_str);
    first_step = False

    model.measurements = emulation.measurements
    # update T0 parameter with final emulation state
    T_final = emulation.display_measurements('Simulated')['Tzone'][-1]
    parameter_data['T0']['Value'] = variables.Static('T0', T_final, units.K);
    opt_problem.Model.parameter_data = parameter_data
    model.parameter_data = parameter_data

[David Blum]

Hey Tom, 

Thanks for posting the results.  I'm glad things are working better for you.  A couple of things:

I've found that there can be those control input oscillations the moment a state hits a constraint (notice the oscillations are when the temperature reaches 20 C).  You may try adding to your objective a penalty on oscillations of the control input (as is found in control literature).  For a scaler input, this would look like k*u^2 (k is a gain and u is the input signal - for a vector of inputs, it would look something like u*K*u_transpose where K is a gain matrix).  You can do this by adding this term in your Modelica model.  It may help to define a new output, maybe "J", and switch this with 'input_q_flow' in your instantiation of the optimization problem.  So in the end, J = input_q_flow + k*input_q_flow^2.  You'll have to tune k a bit.  The second term should be relatively small compared to the first (so that you are still minimizing your energy input), but not insignificant enough that it doesn't penalize larger oscillations.  Let me know if you find this works.

To clarify, if you had ncp = 3 you were getting an infeasible problem result from IPOPT?  ncp increases the number of support points for the collocation method.  But I wouldn't imagine changing from 6 to 3 would make the problem infeasible.

Best Regards, 

Dave

[Tom Stesco]

Hi Dave,

Thanks for adding some ideas about the oscillation. I'm working on adding constraint relaxation for the temperature bounds now via the objective function and perhaps can address this also.

I should clarify that changing n_cp collocation points parameter does not immediately cause IPOPT optimization to become infeasible, but rather in the first MPC iteration when the optimal trajectory for the horizon (72 hour in this case) is input to the emulation for a single time step (1 hour in this case), then the Tzone is extracted and set as the new T0 for the next MPC iteration, Tzone would not always be within the T constraints (i.e. just below T_min) which caused the this next MPC iteration to be infeasible with IPOPT.

It could be that without having the NLP evaluated at (or sufficiently near) t=1hr via collocation, the interpolated polynomial gives a result outside the hard constraints because the T_min constraint at t=2 is slightly higher. I'm currently looking at some references on the LG/LGR collocation discretization in JModelica.org and how the results are extracted. In any case relaxing the constraint should solve this issue and provide better behaviour.

[dhb...@gmail.com]

Hey Tom, 

Yes, I've had issues in the past with the initial state of the optimization being infeasible.  Certainly implementing soft constraints is an option, as is dynamically expanding the hard constraints in a cone around the initial point.

Best Regards, 

Dave