Problem applying fminsearch to piecewise numerical solution genereated by ode45

I'm trying use fminsearch to optimize chemical reaction process in terms of reactant and product components.

A + 3B ---> C

The overall process is piecewise, composed of 5 small processes that alternate between semi-batch and

batch as follows:

1) Semi-batch 2) Batch 3) Semi-batch 4) Batch 5) Semi-batch

During semi-batch, component B is admixed into the reactor at constant volumetric flow rate, so the total volume changes (increasing). During batch, component B is no longer being added, so the total volume in the reactor is constant. Hence, I need to use Ode45 five times, since each batch or semi-batch process is modeled by a separate set of coupled differential equation.

Here's the code I have so far. Note that the results of the first use of Ode45 are used as the initial conditions for the next use of Ode45, and so on.

%----------Declare constants

A = 10^15; % Frequency factor % Frequency factor [=] 1/min
Ea = 100; % Activation energy %  Activation energy [=] kJ/mol

%----------Solve coupled system of differential equations

Vo = 3.5;
Cao = 30/Vo;
Cbo = 0;
Cco = 0;

[t0 f0] = ode45('sim3a',[0 70],[Vo Cao Cbo Cco],[],A,Ea);

V1 = f0(:,1);
Ca1 = f0(:,2);
Cb1 = f0(:,3);
Cc1 = f0(:,4);

[t1 f1] = ode45('sim3b',[70 90],[V1(end) Ca1(end) Cb1(end) Cc1(end)],[],A,Ea);

V2 = f1(:,1);
Ca2 = f1(:,2);
Cb2 = f1(:,3);
Cc2 = f1(:,4);

[t2 f2] = ode45('sim3c',[90 160],[V2(end) Ca2(end) Cb2(end) Cc2(end)],[],A,Ea);

V3 = f2(:,1);
Ca3 = f2(:,2);
Cb3 = f2(:,3);
Cc3 = f2(:,4);

[t3 f3] = ode45('sim3d',[160 180],[V3(end) Ca3(end) Cb3(end) Cc3(end)],[],A,Ea);

V4 = f3(:,1);
Ca4 = f3(:,2);
Cb4 = f3(:,3);
Cc4 = f3(:,4);

[t4 f4] = ode45('sim3e',[180 280],[V4(end) Ca4(end) Cb4(end) Cc4(end)],[],A,Ea);

V5 = f4(:,1);
Ca5 = f4(:,2);
Cb5 = f4(:,3);
Cc5 = f4(:,4);

Okay, so the above code works fine by itself. It executes and outputs beautiful graphs that make sense. The problem arises when I try to optimize the model: I attempt to fit it to experimental data by wrapping the whole thing into an fminsearch function. The two parameters that are optimized are Frequency factor A and Activation energy Ea.

I apply the principle of minimizing the sum of residuals squared (SSR). For this model, I only know the initial and final concentrations as data points (nothing in between). I'm attempting to use fminsearch to tweak with A and Ea such that my model's final concentrations match my experimental final concentrations.

function fminsearch = sim3Cmdb(x)

%----------Declare constants

A = x(1); % Frequency factor [=] 1/min
Ea = x(2); % Activation energy [=] kJ/mol

%----------Solve coupled system of differential equations

... Same exact stuff as above
...
...
...
...

%----------fminsearch

Ca_exp  = 1.32; % Experimental final value component A [=] mol/L
Cb_exp = 0.5; % Experimental final value component B [=] mol/L
Cc_exp = 3.26; % Experimental final value component C [=] mol/L

SSR = (Ca5(end) - Ca_exp)^2 + (Cb5(end) - Cb_exp)^2 + (Cc5(end) - Cc_exp)^2;

fminsearch = SSR

end

My problem: Once I execute the simulation, fminsearch begins iterating; it then gets slower and slower and slower. I left the simulation running all night, only to see in the morning that my PC ran out of memory.

Any ideas on how to fix this? I'm pretty new to Matlab.

[Star Strider]

If you have the Statistics or Optimization Toolboxes, I suggest instead using nlinfit or lsqcurvefit respectively. Then create an objective function that takes A and Ea as parameters (so do not define them as constants in your objective function), then fit the final output of your objective function to your data.

That said, you are fitting two parameters to three data points, so there are no guarantees that you will either be able to determine A and Ea with either accuracy or precision, or that your estimates for them will be unique in different simulation runs.

Answer 1

[Matt Tearle]

It looks like you're on the right track.

I'm a little suspicious of the use of the variable name fminsearch in sim3Cmdb and the out of memory issue. Variables in functions are local, but if fminsearch is calling sim3Cmdb which then references something called fminsearch... Even if it isn't causing problems, it's definitely confusing.

So first thing: remove the line

fminsearch = SSR

and change the function declaration to

function SSR = sim3Cmdb(x)

(Also, I don't know if the missing semicolon on the line fminsearch = SSR is actually in your code, or just a cut/paste artifact. But if you were missing the semicolon, that would cause a huge slowdown, due to the gajillion printouts to the Command Window.)

Before trying to do the minimization, you should test and play a bit with sim3Cmdb. First just call it

sim3Cmdb([1e15,100])

and make sure you get back a single positive number.

Then, given that you have two parameters, you could perhaps try evaluating sim3Cmdb on a grid and plotting the result as a surface or contour plot. This should give you an idea of where the minimum occurs. Ugly but effective approach:

n = 50;  % number of points in grid -- adjust as desired
A = linspace(Amin,Amax,n);    % fill in appropriate values for Amin/Amax
Ea = linspace(Eamin,Eamax,n); % ditto Eamin/Eamax
err = zeros(n);
for k = 1:n
  for j = 1:n
    err(j,k) = sim3Cmdb([A(k),Ea(j)]);
  end
end
surf(A,Ea,err)

The functions tic and toc may be of use for determining how long this should take to run. Do

tic, sim3Cmdb([1e15,100]); toc

to see how long it takes to run your function once. (If it's quick, do it a few times to get an idea of the variability, and ignore the first run as there's often a bit of overhead. If it's at least several seconds, just take that value, because the overhead and variability will be negligible in comparison.) The code above will take n^2 times as long. So choose n accordingly.

Once you have a decent idea of where the minimum should be, you can use fminsearch to locate it precisely.

One last comment: seeing parameters like 1e15 and 100 together makes me squirm. I don't know what your ODEs look like, and you said your function is working OK and giving a good-looking plot, so maybe there's no problem, but I'd be wary of the possibility of numerical precision issues. Is there any way to rescale units? If you have a situation where, for example, you were adding terms that were O(A) to terms that were O(Ea), you're asking for roundoff error. Just something to consider.