parameters estimation

Question

I have to estimate a parameter (K), but I don't know how I can do it. I think by a regression model (minimum least square?), but I'm not sure. The system is:

dX1/dt = f(t)*X2*X1
dX2/dt = -K*f(t)*X2*X1

where:

X1 and X2 are the state variables, f(t) is a time varying function and K is the parameter that I want to estimate. I know the values of X1 and X2 each sample-time Ts=0.5sec. I'm using Matlab.

Answer 1

So, there are two ways to solve this problem.

The easy, non-rigorous way to solve this problem is to create a function that calls a MATLAB ODE solver using $K$ as a parameter, and returns the solution $(x_{1}(t), x_{2}(t))$ for all times corresponding to measured data points. Then use this function to construct a sum-of-squared error objective function, and minimize using MATLAB's unconstrained optimization solver.

The rigorous way to solve this problem is to look at the literature for solving optimization problems with embedded dynamical systems. As far as I know, there exists no MATLAB solver that will take an optimization problem with an embedded dynamical system as input and return a rigorous global $\varepsilon$-optimum. Such methods normally require a nonconvex global optimization solver (which MATLAB does not have, to my knowledge), interval arithmetic libraries (which do exist for MATLAB; INTLAB is one such library), automatic differentiation tools (which also exist for MATLAB using operator overloading and object-oriented language features), and a lot of mathematical theory.

Since people tend not to want to bother with all of the mathematical theory (there is a lot of it), the first approach tends to be favored by practitioners (it's what I was taught as an undergrad), even though it only yields a local optimum due to the nonconvexity of the problem (in most cases, anyway). Occasionally, the second approach is useful for parameter fitting; both of my advisers used this approach for a chemical kinetics application because the first approach was insufficient. (I'm struggling to remember the name of the paper, though, and being pressed for time, I'll have to come back to this answer later.)

Answer 2

Kalman filter (and its variants) can be used to solve this problem.

Answer 3

You have a so called grey-box system identification problem. If you have the System Identification Toolbox installed, you are basically done.

http://www.mathworks.se/products/sysid/description5.html

http://www.mathworks.se/help/ident/ug/estimating-nonlinear-grey-box-models.html

Answer 4

If the parameter to be found is a single number, then you can do brute-force prediction calculations based on a lot of K's (as Geoff suggests), then plot the difference between computational predictions and measured data.

This gives a "complete picture" of how changing the parameter will influence all the measurements, and is not that costly since solving thousands of ODE's takes only a few minutes on modern computers. You can then visually inspect the results to see what K to choose and how sensitive the results are to changes in K - ie., how certain/uncertain you should be about your result.

At this point, you can do some simple least-squares or statistical analysis on the results if you want, but you'll probably get a better idea of what's going on by just looking at it.

Here's some code to do the brute force calculations and plotting! I release this code to the public domain.

% Some randomish input data (replace with your times tt and function f)
tt=0:.01:1;
Ktrue = 2.31; 
ff =  ifft(fft(exp(randn(1,length(tt)))).* ...
              fft(exp(-(tt-5).^2/(1.5)^2)))';
ff = ff/mean(ff);
f = @(t) interp1(tt,ff,t);

% manufacturing state data X1,X2 (replace with your X1,X2)
X1initial = 0.93; 
X2initial = 1.87;
odefuntrue = @(t,Y) [f(t)*Y(1)*Y(2), -Ktrue*f(t)*Y(1)*Y(2)]';
[T,Ytrue] = ode45(odefuntrue,tt,[X1initial,X2initial]);  
X1 = Ytrue(:,1); 
X2 = Ytrue(:,2);

% Solve the ODE for a lot of different values of K
K=0.01:0.01:5;
Ts = zeros(length(tt),length(K)); 
X1s = Ts; X2s = Ts;
for jj=1:length(K)
    odefun = @(t,Y) [f(t)*Y(1)*Y(2), -K(jj)*f(t)*Y(1)*Y(2)]';
    [T,Y] = ode45(odefun,tt,[X1(1),X2(1)]); 
    X1s(:,jj)=Y(:,1);
    X2s(:,jj)=Y(:,2);
end

% Compute the difference between the measured state, (X1,X2), and
% the predicted state for all different K's, (X1s,X2s)
errors1 = X1*ones(1,length(K)) - X1s;
errors2 = X2*ones(1,length(K)) - X2s;

% Plot the full errors at all times for different values of K
figure;
subplot(1,2,1); 
imagesc(K,tt,errors1);
title('error in X1');
xlabel('K');
ylabel('t');
colorbar;

subplot(1,2,2); 
imagesc(K,tt,errors2);
title('error in X2');
xlabel('K');
ylabel('t');
colorbar;

And here's the results. You can see how the error is zero along the vertical line $K=2.31$ corresponding to the true $K$, and the horizontal line $t=0$, and becomes nonzero as $K$ changes away from the true value and time increases.