ODE Function Parameter Fit - Index exceeds bounds/Array arguments do not agree

5 views (last 30 days)
Fitting parameters for an ODE function using this tutorial: https://www.mathworks.com/help/optim/ug/fit-ode-problem-based-least-squares.html
I have an ODE with manually-tuned rate constants Ak(1) to Ak(40) that finds material flows F(1) to F(45) over space Wspan. I need to automatically tune these parameters so I can test the effects of different initial values Fo(1) to Fo(45).
I used the linked tutorial to code parameter tuning: https://www.mathworks.com/help/optim/ug/fit-ode-problem-based-least-squares.html. I copied the tutorial closely to minimize bugs. The simplified code is shown below. solpts and RtoODE are in separate .m files.
%list experimenal end values
Ftrue(1) = 0;
...
Ftrue(45) =0;
%list experimental start values
Fo(1) = 0;
...
Fo(45) = 0.014533;
Wspan = [0 6]; %start and end points of ODE solution -> solve from W = 0 to W = 6
Ak = optimvar('Ak',40,"LowerBound",0,"UpperBound",200); %create vector for 40 parameters
myfcn = fcn2optimexpr(@NTPReactor10, Ak, Wspan, Fo,'OutputSize',[1,45]); %plug into ODE function RtoODE
% RtoODE takes parameters Ak and initial material flows Fo to find flow F
% F calculated over W = 0 to W = 6
obj = sum(sum((myfcn - Ftrue).^2)); %least-squares objective function
prob = optimproblem("Objective",obj); %set optimization problem to minimize obj
%set initial guess to 50 for all Ak
Ak0.Ak = [50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50];
[Aksol,sumsq] = solve(prob,Ak0); %solve optimization problem
function solpts = RtoODE(Ak,Wspan,Fo) %runs RtoODE with optimization variable Ak
soltn = ode45(@(W,F)NTPReactor10(W,F,Ak),Wspan, Fo); %runs original ODE function NTPReactor10 with parameter Ak
finalval = deval(soltn,Wspan); %finds value of RtoODE at Wspan = 6
end
function dFdW = RtoODE(Ak,W,F) %extremely simplified sketch of function
%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
With this code, I get the error:
Index in position 1 exceeds array bounds (must not exceed 1).
Error in NTPReactor10 (line 14)
Flow2 = F(2,:);
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 ReactorModelTestRun9 (line 144)
[Aksol,sumsq] = solve(prob,Ak0);
Caused by:
Failure in initial objective function evaluation. FMINCON cannot continue.
I assumed my OutputSize had been entered incorrectly, and should be [45,1] instead. When I switched it to [45,1], I got the error:
Error using optim.internal.problemdef.ElementwiseOperator/checkIsValid
Argument dimensions 45-by-1 and 1-by-45 must agree.
Error in optim.internal.problemdef.Minus.getMinusOperator
Error in -
Error in ReactorModelTestRun9 (line 141)
obj = sum(sum((myfcn - Ftrue).^2)); %least-squares objective function
I assumed [45,1] was correct, but it didn't align with the Fo array that I used in my original (manually-tuned) function. Substituting transpose(myfcn) or transpose(Ftrue) into obj gave the original error (index cannot exceed 1). Substituting transpose(Fo) in fcn2optimexpr gave the second error (1x45 and 45x1).
Given that the original ODE function NTPReactor10 functions with Wspan and Fo as inputs, how can I get fcn2optimexpr to run the function and have the arrays align?
I am an inexperienced user, please feel free to link existing Answers that may be applicable.

Accepted Answer

Kelsey Chambers
Kelsey Chambers on 21 Jan 2022
Edit: fixed initial/final value matrix and ode45 syntax. The following code now works
%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);
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

More Answers (0)

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!