I am working on a system identification problem with poor quality input data.
The data describes a heating system, with thermostat temperatures, thermal power inputs, supply/return water temperatures etc.
Some issues with the data are:
frequently missing data
spiking noise caused by activation of heaters
discretization errors
low frequency disturbances (mainly due to sunlight, but other sources exist)
A year ago when I tackled this problem, I kind of rolled my own algorithm based on DSP.resample and tuned it by hand. One of the problems it had was it introduced wiggles/oscillations. So I probably also need antialiasing.
Now I want to setup a proper processing pipeline that I can deploy in production.
I am mainly just looking for some suggestions on algorithms / packages to handle all these problems. Unfortunately, there isnât much I can do about the data source because these are systems used by real people (pesky users right )
If you have any other recommendations I would also be very happy to hear those.
Missing data can be handled by predictor/corrector methods (such as Kalman or particle filters) by not perfoming the measurement correction when there is no data available, only the prediction. You can even handle arbitrary intervals inbetween the obtained data points, hereâs an example
Hereâs another example that has missing data
A common approach to handle this in gaussian filtering is to evaluate the likelihood of the obtained measurement, and if the measurement is several standard deviations away from the expected one, discard it (or inflate the measurement covariance matrix). The last tutorial above computes the Z-score used to evaluate the measurement, but never discards any data.
What does this mean? Quantization?
Something like this?
or any of those?
If you have a state estimator of the system that handles all of those things, system identification can be performed by, e.g., maximizing the likelihood of the data
You may also find it useful to have a look at some methods and tools for outlier detection such as Hampel outlier detection/filter. An implementation in Julia can be found at [ANN] HampelOutliers.jl. Standard filtering/smoothing techniques may unnecessarily project the value of such outlier to the neighbouring samples.
I have indeed studied those, it is a very nice overview. I think state estimation is relatively straightforward in this problem, the tricky part is forecasting states and disturbances. What I couldnât figure out is how to deploy those disturbance models to an actual controller that does MPC. For the sunshine model itâs clear, just extrapolate the timeseries and forecasted solar irradiance (I have that).
But for the low frequency disturbance? The convention in literature seems to be to just treat it like a constant parameter and assume it stays constant in the future. Maybe you have some thoughts on that?
Yes that would be ideal.
Models identified with methods that optimize the next-step ahead prediction diverge too quickly if I simulate longer timeseries. So if I am interested in simulation performance the methods in LowLevelParticleFilters.jl donât necessarily seem applicable here? (if this is a wrong statement please correct me)
The main goal is to get a model that performs well in MPC with a horizon of 24 hours. That behaviour tends to de dominated by what the weather and heating system do, for which I have âgoodâ data and forecasts. I probably should have mentioned itâs still the same problem as my previous post
Its reasonably straightforward if the MPC tools are in the state-space form, as it is the case for most modern MPC software. You just build your MPC with the augmented model. The augmented model is your (deterministic) model of the heating system âconcatenatedâ with all your disturbance models (that is, the state vector will be a concatenation of the heating system states and the disturbance states). Be sure to disable any builtin model augmentation if you DIY this part (for exemple, nint_ym=0 and nint_u=0 options in ModelPredictiveControl.jl).
Yes, a low frequency disturbances is generally approximated as constant value that can be added anywhere in the deterministic model: at the manipulated inputs, at the measured output or a constant parameter used anywhere between the two. The disturbance model of a constant value is Gaussian white-noise integrated only once.
I will let Fredrik answer the other questions related to LowLevelParticleFilters.jl. It seems like a really interesting modeling and control problem!
This is indeed how I implemented it. I think I was trying to understand that when youâre looking forward like in MPC you need good deterministic models to define the âknown/forecasted disturbancesâ. I have good (forecasted) data on these disturbances, but the key question is how do they affect the states.
In the disturbance gallery of LowLevelParticleFilters.jl, I think we talk about âunmeasured disturbancesâ only. These do not seem to be useful in the deterministic MPC case beyond estimating the initial condition, because there is no âcorrectâ step as there is no measured data. I am an EE and not a control engineer so please correct me if thatâs wrong.
This is a great package! Because I need to generate C code for an embedded system, I chose to implement MPC as an optimization model in cvxpy and generate C code with cvxpygen. There does seem to be LinearMPC.jl now, I should try that out
Iâm also an EE, but it morphed into a weird mix of EE+ChemE+BioTech+ControlE+SoftwareE. Iâm no longer sure what I am .
Yes, the stochastic models are for the unmeasured disturbances. But these models will affect the MPC predictions/forecasts beyond the initial state. For exemple, in the 3 plots below, the dots are the predictions at the discrete-time k=90 for the stochastic system:
\frac{y_s(z)}{e(z)} = \frac{1}{1-az^{-1}}
with a=0, 0.5 and 1, respectively from top to bottom.
The a=1 case is integrated Gaussian white noise, and you can see that predictions are kept constant in the future.
If you have access to a known/measured disturbance, you first need to determine how it affects your deterministic model (for exemple, added at the manipulated inputs or a load disturbance). This measured disturbance will then be used for feedforward compensation in the MPC controller, denoted with the vector \mathbf{d} in MPC.jl. At each control period, you will have access to your known disturbance \mathbf{d}(k). But there is still one thing to determine: how you predict \mathbf{d} in the future? The simplest solution is assuming a constant value in the future: \mathbf{\hat{d}}(k+1) = \mathbf{\hat{d}}(k+2) = \dots = \mathbf{d}(k). But the exact same theory as above can also be applied here. A constant \mathbf{d} in the future is just assuming that \mathbf{d} is the output of an integrated white noise. It can be modelled with any system from aforementioned disturbance gallery, but you will need to implement the forecasting yourself.
Yes, the is no code generation for now. I should really work on building an interface with LinearMPC.jl.
Edit: it made me think @darnstrom, am I right in thinking that custom prediction for measured disturbances are not supported right now in LinearMPC.jl (itâs really similar as the implementation of reference preview) ? If yes, it would be a great addition to the package.
Itâs not so much a convention as an explicit model that says itâs constant. If you use the model \dot{x}_d = 0 + w
for your low-frequency disturbance, i.e., an integrator of noise, you are explicitly modeling something that has zero deterministic dynamics (constant), only random fluctuations that are equally likely to go in any direction. There are more elaborate disturbance models that do not have zero dynamics, and in that case, simulating those forward in time (for forecasting) does not in general lead to a constant evolution.
diverge in what sense?
Whenever you have unmeasured disturbances you are typically not interested in pure simulation performance, there must be some mechanism present to account for the unmeasured disturbances, and the only information available about this is the measurements you do have access to. These must feed into the dynamics somehow, which is what the filtering accomplishes. The approach from your previous post on the subject sounds spot on what you want, filtering with 24h forecast at each point, minimizing the forecasting errors along the trajectory.
You can treat forecasted values (from e.g. a weather forecast) as âmeasurementsâ and incorporate those as any other measurement. Typically, those would have a rather large covariance compared to actual sensor measurements.
If the system has future disturbances with uncertainties, a primitive way to handle that in a ârobust MPCâ is as in https://doi.org/10.3384/ecp201709: we tried to control the level in a lake (hydro power reservoir) by manipulating the flood gates, and the disturbance was precipitation. Our partnering power producer provided estimates of flow into the lake based on the dynamics from precipitation through the catchment area and into the lake. This was based on purchased weather prediction 2-3 weeks into the future, with some 50 different realizations of possible predictions. The idea was to operate the flood gates so that the control requirements were satisfied for the worst case scenario.
My colleague in that publication has since done much more work on similar problems, with more efficient algorithms for MPC.
I would guess that it would be relatively straightforward to use your MPC package to do this.
Interesting approach, it fits nicely to weather forecasting! I took a look a the paper quickly and, the deterministic MPC algorithm could be easily implemented in MPC.jl with currently available features. But scenario trees and multi-objective optimization are not supported right now (except maybe indirectly with a custom cost function + a weighted sum of each objective), so the would not be the case for stochastic MPC. Robust MPC through scenario trees are in my long-term goals, but Iâm prioritizing as complete and as efficients-as-possible ânon-robustâ predictive controllers before entering into this rabit hole
The way we did the robust MPC was simply to simulate into the future (say) 50 models of the system with 50 different disturbance sequences, each with the same control input sequence, then compute the cost function for each system, and use the sum of the 50 cost functions as the cost function in the MPC algorithm.
This can easily be done by augmenting the system with (say) 50 times larger state than the deterministic system, 50 times as many disturbance input sequences, and a corresponding output. The number of control inputs is the same as for the deterministic system. Iâm sure this fits into your MPC algorithm.
This is a very intuitive way to do robust MPC, although not the most efficient. But simple to understand, and (I think) simple to implement.
Coincidently I am currently running through pretty much the same⌠well⌠pain.
With a new data-intensive (sysID) project I said to myself that it is perhaps the time to try to learn and adopt the more DataFrame-ish workflow (resembling they Pythonâs Pandas) instead of handling the time series at the very low level as pure vectors (the habit I had developed in Matlab before the advent of their timetables).
First, DataFrames.jl does not support the date/datetime index by itself. I have found the TSFrames.jl package that brings date/datetime index into DataFrames. Mostly the resulting âobjectâ indeed behaves transparently as DataFrame, but there are some glitches.
And I confess that I am struggling with even the elementary steps. Mainly because I am not yet fluent with that dataframe-ish way of working, but also because even very basic functionality is typically provided by some other packages, and it takes quite some time to find them and learn which are usableâŚ
Oh, how I would love to see a time-series extension of the online book https://juliadatascience.io to time-series workflows. How to do imputing (perhaps using Impute.jl, you may find it useful too), how to remove outliers (I am only learning if can use the Hampel filter that I mentioned here earlier, or is there anything reasonably well maintained?) which tool currently best fits into this workflow for smoothing, resampling, the latter being my the step I am currently struggling with. Then how to do cumsum for the columns of the DataFrame/TSFrame (here I got lost in the rabbit holes here in the discussion forum). RollingFunction.jl is in the to-do (to-study) list to learn if and how I could then apply other filters such as Savitzky-Golay.jl.
I am still enjoying it but I admit it is tempting to go back just to âlong vectorsâ.
I found that if I fit the model on 1-step ahead predictions, that fit doesnât perform well if I simulate that model forwards for 24 hours, even if disturbances are mostly absent. I guess thatâs because the dynamics are so slow and the 1-step ahead prediction doesnât capture that. A good observer is not a good predictor in this case
I think this is the right approach, because it also matches how the model would be used in reality. If itâs not too much effort could you maybe give a rough sketch how you would approach this with LowLevelParticleFilters.jl?
function cost(pars)
kf = FilterConstructor(pars)
c = zero(eltype(pars))
for time in all_timepoints[2:end]
update!(kf, u[time-1], y[time-1])
x_sim, u_sim, y_sim = rollout(kf, state(kf), u[time:time+24])
c += sum((y_sim .- y[time:time+24]).^2) / length(y_sim)
end
return c / length(all_timepoints)
end
The sum of squared error cost can be replaced by, e.g., loglikelihood of prediction errors if you step forward with repeated calls to predict! instead of the rollout call above, but this should be a reasonable starting place.
)
.
