Problem 641360 · medium · Level 06 Heuristics & Optimization

A Model of the Hot-Water Cylinder, Tested on a New Day

system identification · least squares · ARX model · validation · numpy · simulation error

A heat pump's controller needs a model of its hot-water cylinder. Heat flows from the heating coil into the water and from the water into the sensor pocket, two stages of lag, so a second-order model is a sensible guess. Once a minute the log records the heater power u[k] (kW) and the temperature rise at the sensor y[k] (°C above the start). The model is

y[k] = a1 * y[k-1] + a2 * y[k-2] + b1 * u[k-1] + b2 * u[k-2]

with four unknown numbers. Every k >= 2 of a log gives one equation: the row [y[k-1], y[k-2], u[k-1], u[k-2]] times [a1, a2, b1, b2] should equal y[k]. Stack the rows into a matrix A and the targets into a vector b; the least-squares fit is the parameter vector that makes the sum of squared misses smallest, numpy.linalg.lstsq(A, b, rcond=None)[0].

A model must be judged on data it has not seen. Fit on the training log (u_train, y_train), then on the test log (u_test, y_test, recorded on another day) compute two errors over k = 2, ..., len(y_test) - 1:

  • the one-step error: predict y[k] from the measured y_test[k-1], y_test[k-2] and the inputs;
  • the simulation error: run the model on its own from s[0] = y_test[0], s[1] = y_test[1], each step using its own previous outputs s[k-1], s[k-2] and the measured inputs, and compare s[k] with y_test[k].

Each error is an RMS: sqrt(mean((prediction - y_test[k]) ** 2)). Write identify(u_train, y_train, u_test, y_test) that returns (a1, a2, b1, b2, rms_one_step, rms_simulation).

Examples

Input:  u_train = [1, 1, 1, 0, 0, 0, 2, 2, 2, 0, 0, 0]
        y_train = [-0.05, 0.0, 0.45, 1.23, 1.97, 2.33, 2.51, 3.06, 4.41, 6.09, 7.3, 7.99]
        u_test  = [1, 1, 1, 2, 2, 2, 3, 3, 3, 0]
        y_test  = [-0.02, -0.03, 0.51, 1.34, 2.55, 4.25, 6.26, 8.58, 11.31, 14.2]
Output: (1.5270838093646282, -0.5309992338826182, 0.26473230769363126, 0.2289205076808326, 0.0650916888407576, 0.49195478522803704)
Explanation: one minute ahead the model is within 0.07 °C of the new day's readings, about the sensor's
own noise. Left to run on its own for eight minutes, its small errors add up to 0.49 °C.

Input:  u_train = [0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2]
        y_train = [0.04, -0.07, -0.02, -0.04, -0.05, -0.04, -0.03, -0.01, -0.05, 0.02, 0.57, 1.8]
        u_test  = [1, 1, 1, 2, 2, 2, 1, 1, 1, 1]
        y_test  = [-0.05, 0.05, 0.49, 1.29, 2.5, 4.27, 6.14, 7.97, 9.34, 10.37]
Output: (0.33206379433468125, 0.44168055224946673, 0.2927213758628896, 0.5082236372292311, 2.022597057482767, 3.7421794295128215)
Explanation: the heater was off for most of this training log, so it holds almost nothing about how the
cylinder responds. The fit still finds four numbers, but they fail badly on the new day.

Some tests use sim_cylinder(n, seed, noise), defined for you, which returns a simulated log (u, y) of n minutes with readings rounded to 0.01 °C after adding noise of standard deviation noise. Try identify(*sim_cylinder(300, 1, 0.05), *sim_cylinder(100, 2, 0.05)) with Run, and plot a test log's y next to your simulation.

Constraints

  • 4 <= len(u_train) == len(y_train) and 3 <= len(u_test) == len(y_test); every training log in the tests has a unique least-squares solution
  • answers are compared with a tolerance of 1e-6

Goals

  • Set up a second-order input-output model as a least-squares problem and solve it with numpy.linalg.lstsq
  • Judge the model on fresh data it was not fitted to
  • Tell apart the one-step prediction error and the far larger error of a free-running simulation
Starting Python…