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

9 views (last 30 days)
Kelsey Chambers on 21 Jan 2022
Commented: Kelsey Chambers on 23 Jan 2022
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
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?
Kelsey Chambers on 23 Jan 2022
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!

Torsten on 23 Jan 2022
Edited: Torsten on 23 Jan 2022
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 on 23 Jan 2022
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!

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by