Fitting a fixed nonlinear ODE model to closed-loop greenhouse data: which estimator avoids the controller-induced bias?

Question

As an Amazon Associate, I earn from qualifying purchases. For more practice problems like this, see Schaum’s Outline of Statistics, 6th Edition.

1. What is being asked?

A greenhouse climate model is already fixed:

28 ordinary differential equations (ODEs) describing temperature, humidity, CO₂, etc.
≈200 physical parameters – the student is only allowed to scale a handful of them (multiplicative factors).

The data come from a closed‑loop experiment: a controller (the climate computer) continuously reads the same variables that we want to fit and decides how much heating, ventilation, shading, etc., to apply. Consequently the input signals are correlated with the disturbances (measurement noise, unmodelled dynamics, weather fluctuations).

The student tried a deterministic “output‑error” (OE) fit (simulation error minimisation) and obtained a biased, poorly predictive model.

Question:
Given a known nonlinear ODE structure and closed‑loop data, which estimation strategy should be used instead of plain OE least‑squares?

Possible candidates mentioned are

prediction‑error (PE) with a stochastic noise model,
multiple‑shooting,
Bayesian parameter estimation,
an estimator that explicitly uses the known controller,

and the student would also like a concrete illustration on a physical model.

2. Step‑by‑step answer

Below we develop a complete, practical recipe for this class of problems.
We start from the statistical description of the data‑generating system, identify why OE is biased, and then construct an estimator that eliminates—or at least greatly reduces—that bias.
Every step is written out explicitly so you can reproduce it on the greenhouse model (or any other nonlinear ODE system).

2.1 Write down the closed‑loop system in statistical form

Symbol	Meaning
(x(t)\in\mathbb R^{n})	State vector (solution of the 28 ODEs).
(\theta\in\mathbb R^{p})	Vector of unknown parameters (the 8 scaling factors).
(u(t)\in\mathbb R^{m})	Vector of manipulated inputs (heater power, vent opening, screen position).
(w(t))	Process disturbance (unmodelled heat gains, sensor drift, etc.).
(v(t))	Measurement noise on the observed outputs.
(y(t))	Measured outputs (temperature, humidity, CO₂).
(C)	Output matrix (selects the observable components).
(\mathcal{K}(\cdot))	The known controller (digital, possibly PI/PID).

The continuous‑time ODE model (discretised at the 5‑minute sampling period, (\Delta=5) min) is

[ \begin{aligned} \dot x(t) &= f\bigl(x(t),u(t),\theta\bigr) + w(t) ,
y(t) &= C\,x(t) + v(t) . \end{aligned} \tag{1} ]

Because the system runs in closed loop, the inputs are generated by the controller

[ u(t) = \mathcal K!\bigl(y(t),\; r(t)\bigr), \tag{2} ]

where (r(t)) are the reference signals (the set‑points for temperature, humidity, CO₂).
Equation (2) makes the inputs statistically dependent on the disturbances (w(t),v(t)); consequently a naïve OE estimator that treats (u(t)) as exogenous will be biased (Forssell & Ljung, 1999).

2.2 Why the ordinary output‑error (OE) estimator fails

The OE estimator solves

[ \hat\theta_{\text{OE}} = \arg\min_\theta \sum_{k=1}^{N}| y(k) - C\hat x(k;\theta) |^2 , \tag{3} ]

where (\hat x(k;\theta)) is the simulated state trajectory obtained by integrating (1) with the measured inputs (u(k)) and ignoring any noise model.

Because (u(k)) depends on (y(k)), the term (C\hat x(k;\theta)) is correlated with the very disturbances we are trying to “explain”. The minimisation therefore drives the parameters toward values that absorb part of the disturbance (i.e., a biased solution).

The bias does not disappear even with infinitely long data, unless the disturbance never reaches the inputs (the special case where the controller is feed‑forward only).

2.3 General remedy: Closed‑loop prediction‑error (PE) identification

The classic solution for linear systems is the prediction‑error method (PEM) with an explicit stochastic noise model and, if needed, instrumental variables.
The same idea extends to nonlinear ODEs.

2.3.1 Augment the model with a stochastic noise description

Add a dynamic noise model (usually a linear filter) that captures the correlation structure of (w(t)) and (v(t)). A minimal but effective choice is a first‑order ARMA (or a low‑order AR) process:

[ \begin{aligned} e(t) &\triangleq y(t) - C\,x(t) ,
e(t) &= H(q,\psi) \, e(t) + v(t), \qquad H(q,\psi) = \frac{B(q,\psi)}{A(q,\psi)}, \end{aligned} \tag{4} ]

where (q^{-1}) is the forward shift operator, (\psi) are the noise‑model parameters, and (v(t)) is assumed white Gaussian.

The full stochastic model (sometimes called a state‑space‑with‑noise model) becomes

[ \boxed{ \begin{aligned} \dot x(t) &= f\bigl(x(t),u(t),\theta\bigr) + w(t),
e(t) &= H(q,\psi)\, v(t),
y(t) &= C\,x(t) + e(t),
u(t) &= \mathcal K!\bigl(y(t), r(t)\bigr) . \end{aligned}} \tag{5} ]

2.3.2 One‑step‑ahead predictor (the quantity that the PE cost uses)

From (5) we can write the one‑step‑ahead predictor of the output, conditional on all past data (\mathcal I_{t-1}):

[ \hat y(t|t-1;\theta,\psi) = C\,\hat x(t|t-1) + \hat e(t|t-1), \tag{6} ]

where

(\hat x(t t-1)) is obtained by integrating the ODE with the predicted input

[ \hat u(t-1) = \mathcal K!\bigl(\hat y(t-1|t-2), r(t-1)\bigr), ]

i.e. we propagate the controller inside the predictor, not with the measured (u).

(\hat e(t

t-1)) is the output of the linear filter (H(q,\psi)) applied to past prediction errors (the usual Kalman‑type predictor for an ARMA filter).

The predictor is unbiased because the input used inside the predictor is generated by the same controller that generated the real inputs, but it depends only on past data – the future disturbance cannot influence it.

2.3.3 Prediction‑error cost

Define the prediction error

[ \varepsilon(t;\theta,\psi) \;=\; y(t) - \hat y(t|t-1;\theta,\psi) . \tag{7} ]

The PEM estimator solves

[ (\hat\theta_{\text{PE}},\hat\psi_{\text{PE}}) \;=\; \arg\min_{\theta,\psi} \; \frac{1}{N}\sum_{t=1}^{N}|\varepsilon(t;\theta,\psi)|^{2}. \tag{8} ]

Because the predictor respects the closed‑loop relationship (2), the estimator is consistent under the usual regularity conditions (identifiability, persistent excitation, correct noise model order).

2.4 Practical implementation steps for the greenhouse case

Below we list the concrete actions you have to perform, with enough detail to code them in MATLAB, Python (SciPy + CasADi), or Julia (DifferentialEquations + Optim).

Step	Action	Details & Tips
1	Gather the data	Build a single time‑indexed table containing: - Measured outputs (y(t)) (T, RH, CO₂). - Reference set‑points (r(t)). - All actuator commands (u(t)) as recorded (they will be used only for validation). - Optional exogenous weather variables (outside temperature, solar radiation) – treat them as known disturbances that enter the ODE as additional inputs (they are not controlled).
2	Discretise the ODE	Use a high‑order integrator (e.g., 4th‑order Runge‑Kutta) with step size equal to the sampling interval (5 min). For a parameter vector (\theta), define a function `simulate_one_step(x,u,θ,Δ)` that returns (x^{+}).
3	Implement the known controller	Translate the greenhouse controller logic (often a PID on temperature, humidity, CO₂) into a deterministic function `controller(y_prev, r, θc)` that outputs the next input (u). Keep the controller parameters (\theta_c) fixed (they are known).
4	Choose a noise model	Start with a simple AR(1): (e(t)=a\,e(t-1)+v(t)). If residuals after identification are still autocorrelated, increase to AR(2) or ARMA(1,1).
5	Build the one‑step‑ahead predictor (Algorithm)	1. Initialise (x(0)) (e.g., set to the measured state at the first sample). 2. For each time step (t=1\ldots N): a. Compute the predicted input (\hat u(t-1)=\mathcal K(\hat y(t-1	t-2),r(t-1))). b. Propagate the state: (\hat x(t	t-1)=\text{integrate}( \hat x(t-1	t-2),\hat u(t-1),\theta )). c. Compute the filter prediction (\hat e(t	t-1)=a\,\hat e(t-1	t-2)). d. Form (\hat y(t	t-1)=C\hat x(t	t-1)+\hat e(t	t-1)).
6	Define the objective	Implement a function `cost(θ,ψ)` that runs the predictor over the whole dataset, computes (\varepsilon(t)=y(t)-\hat y(t	t-1)), and returns the sum of squares. Add optional regularisation (e.g., ridge) if parameters are near bounds.
7	Choose an optimiser	Because the cost is non‑convex, use a global‑local scheme: - Global stage: differential evolution / particle swarm to locate a promising region. - Local stage: Levenberg‑Marquardt (LM) or BFGS on the same cost (LM works well when you can provide Jacobians). - For the ODE part you can obtain sensitivities (\partial \hat x / \partial \theta) by forward sensitivity equations or by automatic differentiation (AD).
8	Check identifiability	Compute the Fisher information matrix (FIM) at the final estimate: (F = \frac{1}{N}\sum_t \frac{\partial \varepsilon(t)}{\partial \theta}\frac{\partial \varepsilon(t)}{\partial \theta}^\top). Small eigenvalues reveal unidentifiable combinations – drop or re‑parameterise those directions (this explains the “one parameter not identifiable” the student observed).
9	Validate	Closed‑loop simulation (using the measured inputs) – compute MAE on the hold‑out set as a sanity check (should be better than 0.7

Original question: Fitting a fixed nonlinear ODE model to closed-loop greenhouse data: which estimator avoids the controller-induced bias? on Cross Validated (Stats Stack Exchange), licensed CC BY-SA.