Why does 40-variable optimization array exceed maximum array size preference?

I am estimating parameters for a system of ODEs. The goal of my model is to match the results Ftrue of an experiment with input values Fo. The model uses 40 parameters Ak(1) to Ak(40) in its calculations. I used a MATLAB tutorial code to make this optimization problem:

%openExample('optim/FitODEProblemBasedExample')
Ftrue = [0 ... 0]; %experimental end-values
Fo = [0	... 0.014533]; %experimental start-values
Wspan = [0 6];
%soltrue = ode45(@(W,F)diffun(W,F,Ak),Wspan,Fo, Aktrue);

Ak = optimvar('Ak',40,"LowerBound",0,"UpperBound",200);
myfcn = fcn2optimexpr(@RtoODE,Ak,Wspan,Fo,'OutputSize',[1,45]);
obj = sum(sum((myfcn - Ftrue).^2));

prob = optimproblem("Objective",obj);
Ak0.Ak = [50 ... 50];
[Aksol,sumsq] = solve(prob,Ak0);

This code calls on function RtoODE, which calls function diffun.

function solpts = RtoODE(Ak,Wspan,Fo)
sol = ode45(@(W,F)diffun(W,F,Ak),Wspan,Fo);
solpts = deval(sol,Wspan);
end

function dFdW = diffun(~,F,Ak) %simplified sketch of function diffun
%shows relationships between inputs Ak (parameter) and F (material flows)
Flow1 = F(1,:);
...
Flow45 = F(45,:);    
dFdW(1,:) = Ak(1)*Ak(40)*Flow1; %all dFdW outputs are functions of various Ak and F values
...
dFdW(45,:) = Ak(20)*Ak(30)*Flow30 + Ak(1)*Ak(40)*Flow1;
end

When I run the optimization problem, it times out and gives me the following error:

>> FinalOptimization

Error using cat

Requested 45x7x6726720 (15.8GB) array exceeds maximum array size preference. Creation of arrays greater than

this limit may take a long time and cause MATLAB to become unresponsive. See array size limit or preference

panel for more information.

Error in ode45 (line 434)

f3d = cat(3,f3d,zeros(neq,7,chunk,dataType));

Error in RtoODE (line 2)

sol = ode45(@(W,F)diffun(W,F,Ak),Wspan,Fo);

Error in generatedObjective (line 35)

Error in optim.problemdef.OptimizationProblem/compileObjective>@(x)objhandle(x,extraParams)

Error in fmincon (line 552)

initVals.f = feval(funfcn{3},X,varargin{:});

Error in optim.problemdef.OptimizationProblem/solve

Error in FinalOptimization (line 16)

[Aksol,sumsq] = solve(prob,Ak0);

Caused by:

Failure in initial objective function evaluation. FMINCON cannot continue.

The 15.8 GB array seems unrealistic for the size of the optimization problem. I know MATLAB can be used to estimate 30+ parameters in an ODE function, and the function runs without bugs until solve(prob,Ak0). Have I created an infinite loop? Do I need to narrow the Ak bounds further?

[Michael Van de Graaff]

Why does 15.8 GB seem unrealistic? That's just 45x7x6726720*8 bytes.

If that array size seems to large, I would track down where your code is producing those sizes. I see the 45 floating around in your code, but idk where the other two dimensions are coming from.

[Torsten]

What is the total number of measurement data you want to reproduce by your model ?

One more observation:

In Wspan, you must specify the vector of times when measurement data are present.

Just saying Wspan = [0 6] is wrong because it allows the solver to return the solution at times of its own choice.

And you return this arbitrary vector to the optimizer by setting

solpts = deval(sol,Wspan);

[Kelsey Chambers]

Thank you for the responses.

Michael:

1. Sorry--I realized I missed a step between drafting this question and publishing. Original problem: I was estimating 40 reaction rates. I simplified to 7 reaction rates to try and reduce the array size. Rather than Ak(1) to Ak(40), I am now using Ak(1) to Ak(7)
2. I thought 15.8 GB was unrealistic because the tutorial code (https://www.mathworks.com/help/optim/ug/fit-ode-problem-based-least-squares.html) found 3 parameters from 6 initial values in a few minutes. I'm now finding 7 parameters from 45 initial values, which isn't much larger. I may be wrong.
3. I have 45 material components in a reactor, shown by F. F changes along the length of the reactor, shown by W, depending on parameter Ak, which affects the rate of reaction. Seven reactions are taking place, so there are 7 values of Ak. I don't know where 6726720 comes from. I think it is the # of iterations/guesses used by the lsq solver, but I don't know how to reduce it.

Torsten:

I designed my function to plot changing F over reactor length W at steady state.

[W,F] = ode45(@(W,F) diffun(W,F), Wspan,Fo);

Once I ran the function, I would find how many length-steps were in the reactor using size() and then find the values of F(1) to F(45) at the last point (the end of the reactor).

>> ReactorModelTestRun8

>> size(F)

ans =

81 45

>> F(81,:)

I thought deval(sol, Wspan) would find the value of F (sol) at W = 0 and W = 6, or the start and finish of the reactor.

Based on your comment, I tried changing solpts to:

function solpts = RtoODE(Ak,Wspan,Fo)
sol = ode45(@(W,F)diffun(W,F,Ak),Wspan,Fo);
lastval = size(sol);
solpts = sol(lastval(1),:);
end

I encountered the same 15.8 GB timeout.

Using show(prob), I made sure that the issue was with solve():

Solve for:

Ak

minimize :

sum(sum((arg2(2) - extraParams{3}).^2))

where:

arg2 = RtoODE(Ak, extraParams{1}, extraParams{2});

extraParams

variable bounds:

0 <= Ak(1) <= 100

0 <= Ak(2) <= 100

0 <= Ak(3) <= 100

0 <= Ak(4) <= 100

0 <= Ak(5) <= 100

0 <= Ak(6) <= 100

0 <= Ak(7) <= 100

I then tried to change the initial guess of Ak(i) = 50 to numbers that would be closer to the true value.

Ak0.Ak = [60 80 60 10 10 80 80];

Running the code produced this error:

Solving problem using lsqnonlin.

Error using optim.problemdef.OptimizationProblem/solve

Index exceeds the number of array elements. Index must not exceed 1.

Error in FinalOptimization (line 20)

[Aksol,sumsq] = solve(prob,Ak0);

Caused by:

Failure in initial objective function evaluation. LSQNONLIN cannot continue.

How does changing the value of initial guesses within the 0 < Ak < 100 range cause this error? It seems to be in solpts, but since I'm now finding the size of F on each iteration, I don't know how that would be out of bounds.

Thank you again!

Answer 1

[Torsten]

It's hard to see what MATLAB does when using the problem-based approach.

I prefer the solver-based approach because it is much more transparent.

Here is the code for your problem - easy to understand and easy to debug:

    function main
      Fo = [0	... 0.014533]; %experimental start-values 
      Ftrue = [0 ... 0]; %experimental end-values 
      Wspan = [0 6];
      Ak0 = [50 ... 50];
      lbAk = [0 ... 0];
      ubAk = [200 ... 200];
      AK = lsqnonlin(@(x)fun_optim(x,Fo,Ftrue,Wspan),Ak0,lbAk,ubAk)
    end  
    
    function res = fun_optim(Ak,Fo,Ftrue,Wspan)
      [t,sol] = ode15s(@(t,y)fun_ode(y,Ak),Wspan,Fo);
      res = sol(end,:) - Ftrue;
    end
    
    function dFdW = fun_ode(F,Ak) 
      Flow1 = F(1);
      ..
      Flow45 = F(45); 
      dFdW = zeros(45,1);   
      dFdW(1) = Ak(1)*Ak(40)*Flow1; 
      ...
      dFdW(45) = Ak(20)*Ak(30)*Flow30 + Ak(1)*Ak(40)*Flow1;
    end

[Kelsey Chambers]

Thank you for putting this together for me, it has solved my problem for now! Using the explicit least squares resulted in this error for an initial guess of all Ak(i) = 50:

Solver stopped prematurely.

lsqnonlin stopped because it exceeded the function evaluation limit,

options.MaxFunctionEvaluations = 7.000000e+02.

Adjusting initial guesses for Ak resulted in:

Local minimum possible.

lsqnonlin stopped because the size of the current step is less than the value of the step size tolerance.

<stopping criteria details>

It looks like function options can be edited to allow a greater range of estimation. I'm having some trouble extracting the Ak(i) values, but that is a problem for another thread. Much appreciated!