Techniques that match piecewise physical model with empirical data

Question

I have a physical model that describes some experimental data. The problem that I have is that the model equations are defined piecewise and are valid at different stages. I would like to be able to estimate the parameters in the model (in a min least squares sense) but also would need to know when to use which equations (i.e. be able to also estimate the stages). Important to note that the number of stages are not necessarily fixed but the available stages are. So one could repeat stages to obtain a better fit. The model consists in a bunch of straight lines.

$\displaystyle f(x)=\left\{ \begin{array}{ll} c_1x+c_2 & x \in [a,b) \\ c_3x+c_4 & x \in [b,c) \\ ... \end{array} \right. $

Is there any accepted technique to tackle such a task?

I thought about running some sort of recursive program that guesses the stages and then runs Least Squares to optimize the parameters. But maybe there's already something out there...

Thank you!

Answer 1

We can focus this as an optimization problem:

Given a table of $N$ points in $d$ from $x_{min}$ to $x_{max}$ and a generic piecewise continuous function

$$ f(x)=\sum_{k=1}^{n-1}(m_k(x-\mu_k)+c_k)(\theta(x-\mu_k)-\theta(x-\mu_{k+1})) $$

here $m_k, \mu_k, c_k$ are parameters to determine and $\theta(x)$ is the Heaviside step function.

now defining a set of restrictions as

$$ \mathcal{R}=\cases{ m_k(\mu_{k+1}-\mu_k)+c_k = c_{k+1}\\ \mu_1\le \mu_k \lt \mu_{k+1}\lt \mu_n } $$

with the error function

$$ \mathcal{E}=\sum_{k=1}^{N}\left(d_k(2)-f(d_k(1))\right)^2 $$

so the problem can be stated as

$$ \min_{\phi}\mathcal{E}(\phi),\ \ \ \text{s. t.}\ \ \ \phi\in \mathcal{R}(\phi) $$

with $\phi = \{m_k,\mu_k,c_k\}$

Follows a basic MATHEMATICA script showing the procedure.

Clear[f, F]
n = 5;
xmax = 2 Pi;
dx = 0.2;
np = Round[xmax/dx];
d = Table[{k dx, Sin[k dx]}, {k, 0, np}];
f[x_?NumericQ] := Sum[(m[k] (x - mu[k]) + c[k]) (UnitStep[x - mu[k]] - UnitStep[x - mu[k + 1]]), {k, 1, n - 1}]
resteq = Table[m[k] (mu[k + 1] - mu[k]) + c[k] - c[k + 1] == 0, {k, 1, n - 1}];
restrmu = Join[{mu[1] == 0, mu[n] == xmax}, Table[mu[k] < mu[k + 1], {k, 1, n - 1}]];
vars = Join[Join[Table[m[k], {k, 1, n - 1}], Table[c[k], {k, 1, n}]], Table[mu[k], {k, 1, n}]];
restrs = Join[resteq, restrmu];

obj = Sum[(d[[k, 2]] - f[d[[k, 1]]])^2, {k, 1, np}];
sol = NMinimize[Join[{obj}, restrs], vars]

F[x_] := Sum[(m[k] (x - mu[k]) + c[k]) (UnitStep[x - mu[k]] - UnitStep[x - mu[k + 1]]), {k, 1, n - 1}]
f0 = F[x] /. sol[[2]]
gr1 = Plot[f0, {x, 0, xmax}];
gr2 = ListPlot[d, PlotStyle -> Red];
Show[gr2, gr1, PlotRange -> All]