Least squares fitting to a differential equation

Posted in hacks and kludges on Saturday, April 27 2013

Recently, in a class on process modeling, I had the task to fit some model parameters using some experimental data for a CSTR. The model for the CSTR being used was a set of nonlinear differential equations. I looked at the answer given and I don't like it, it uses a bunch of global variables and is kind of a mess. What follows is how I used ode45 in matlab and nonlinear fitting to solve the same problem.

First off there are a couple ways I could have done this, for one I could have linearized the differential equations and come up with a function of the sort $ y = f(x;p)$ where the parameters are p. In this case, however, I chose to use the fmincon function. fmincon minimizes some function f using a gradient method, I want to minimize the sum of squares error with the input being the model parameters.

Before I do all that I should build my model for the CSTR. This is a simple model where $ A \rightarrow products$ in a first order reaction, we assume perfect mixing, and the CSTR is jacketed with some heat transfer coefficient UA. What we are fitting for is the Arrhenius preexponential and the heat of reaction, for some reason we already know the activation energy (though I suppose that could have been fitted for).

Anyways the model is:

$$ \dot{C_{A}} = q (C_{A0}-C_{A}) - k_{0} exp({{-E_{a}} \over {R T}}) C_{A}$$

$$ \dot{T} = q (T_{0}-T) + \frac{H}{\rho C_{p}} k_{0} exp({{-E_{a}} \over {R T}}) C_{A} - \frac{UA}{\rho C_{p} V} (T-T_{c})$$

Where q is the volumetric flowrate, V the volume of the CSTR, T the temperature of the CSTR, T0 the temperature of the inlet, and Tc the temperature of the cooling water in the jacket.

Suppose my outputs are the concentration of A and the temperature, this is a vector $ x$, I need a function that gives me $ x = f(x,t;p)$, which I don't know how to write directly because my model is nonlinear. However I can do this using nested functions and ode45 in matlab. I call this fitfun:

function [ ts, xs ] = fitfun( params, xdata, tdata )
%System Parameters
H = params(2)*5960; %kcal/kmol
k0 = params(1)*9703*3600; %h^-1
UA = 150; %kcal/(m3*C*h)
CA0 = 10;
T0 = 298.15;% K
q = 1;
Tc = 298.15;% K
Ea = 11843; % kcal/mol
R = 1.9858775; % kcal/mol*K
rCp = 500;

% This is the actual model giving \dot{x} = f(x,t)
    function eq = cstr(t, x)
        CA = x(1);
        T = x(2);
        eq1 = q*(CA0-CA) - k0*exp(-Ea/(R*T))*CA;
        eq2 = q*(T0-T) + (H/rCp)*k0*exp(-Ea/(R*T))*CA - (UA/rCp)*(T-Tc);
        eq = [eq1; eq2];

% Using the data given to find the time interval and the initial conditions
t0 = min(tdata);
tf = max(tdata);
space = linspace(t0,tf,length(tdata));
x0 = xdata(1,:);

% Finally ode45 spits out the matching data points
[ts, xs] = ode45(@cstr,space,x0,[]);


Since the differential equations need initial conditions I just pull that from the experimental data, along with the time domain to integrate over.

This is a start but to use fmincon I need to create an objective function. In this case lets suppose I a more concerned with fitting the second equation, the temperature dependence, I can create an objective function that like so:

function [ SSE ] = objective( params, xdata, tdata )
    [ts, xs] = fitfun(params, xdata, tdata);
    SSE = (xdata(:,2)-xs(:,2))'*(xdata(:,2)-xs(:,2));

Now I can actually do the least squares fit, by running a script like the following:

load cstr_data.mat;

% perform least squares fit
objfun = @(params) objective(params, xdata, tdata);
[p fval] = fmincon(objfun, [1 1], [], [], [], [], [0.25 0.5], [1.75 1.5]);

If you are wondering why I didn't roll the objective function and the fit function into one, this was just so I could call the fit function with the final parameters and create some nice fancy plots. This isn't the perfect solution, and it is limited in reusability, since there are some hard coded constants in there, but it works pretty well.