Optimization&ODE problem // Symptom: Index exceeds matrix dimensions (here we go again!)

Dear Community,

Since your help brought me so much forward, I have become more and more ambitious! Right now I am trying my hand at the optimization of a parameter used within an ODE system. For that I read following: https://de.mathworks.com/help/gads/optimize-an-ode-in-parallel.html

However my problem is a little different I believe because the parameter that should be optimized is not a variable output of the ODE system (y in my code, see below) but a parameter (hap) that is used as constant for solving the ODE.

Because of that difference I am not sure how and where to declare this parameter as I want to feed the initial values for this parameter only once I define the optimization expression. And the symptom of this problem is that I get the error message: Index exceeds matrix dimensions

Also, the objective function (defined at bottom of my code) is defined according to a variable (ms) that is different from the variable to be optimized (however its value depends on it) which will at least cause a problem when defining the initial condition of the optimization...

Just to make sure it is clear: the result I expect from this code is the matrix hapA.

Please let me know if you would like me to specify something more to make it easier to understand my problem...

Thanks!!

s=127;
LPool1=zeros(s,1);

% Defining the matrix LPool
for k=1:numel(LPool1)           
        if rem(k,2)==1 && k<98; LPool1(k)= 5*(k-1)/2; elseif rem(k,2)==1 && k>97; LPool1(k)= 5*(97-1)/2+6*(k-97)/2; elseif k<98; LPool1(k)=(k/2)*2+(k/2-1)*3; elseif k>97; LPool1(k)=(96/2)*2+(96/2)*3 +((k-96)/2)*3+((k-96)/2-1)*3;  end             
end
LPool=LPool1.';

y0=[1.1,40];                    
LA=zeros([],1);yA=zeros([],2); hapA=[];

% Solving ODE system with Loop
for k=1:numel(LPool)-1      
hap=hapA(k);
if k<97;[L,y]=ode45(@(L,y) myODE(L,y),LPool(k:k+1),y0); y0=y(end,:); y1(1)=0.47; y1(2)=y0(2); else; [L,y]=ode45(@(L,y) myODE(L,y,hap),LPool(k:k+1),y1); y1=y(end,:); end
LA=[LA;L];yA=[yA;y];hapA=[hapA;hap];
end

mscA=yA(:,1)*10;
ms=sum(mscA);

%%%%%%%%%% OPTIMIZATION %%%%%%%%%%%%%%%%%%%%%%
hap1=[0.017;0.0197]; hap2=[0.0163;0.0192];
hapA0=[repmat(hap1,n1(1),1);repmat(hap2,n1(2),1)];  % Confirmed heat transfer parameter between air and sheet

[hapA] = patternsearch(@objfun,hapA0);

%%%%%%%%%% FUNCTIONS %%%%%%%%%%%%%%%%%%%%%
function dy = myODE(L,y,hap)
u=y(1); Tp=y(2);
dudL = myODE1(L,u,Tp,hap); dTpdL = myODE2(L,u,Tp,hap);
dy = [dudL;dTpdL];
end

function dudL = myODE1(~,~,Tp,~)    
dudL=-hap/(Tp+273.15);  % u is in kg water/kg fiber
end

function dTpdL=myODE2(~,u,Tp,hap)
dTpdL=((1+0.955*u)*(100-Tp)+hap*(100-Tp))/(hap*u);                                                                                           
end

function f=objfun(ms)
msC=20.97;
f=msC-ms;
end    

[Cris LaPierre]

Code is much more readable if you put each expression on its own line. It also makes it easier to debug.

s=127;
LPool1=zeros(1,s);
% Defining the matrix LPool
for k=1:numel(LPool1) 
    if rem(k,2)==1 && k<98
        LPool1(k)= 5*(k-1)/2; 
    elseif rem(k,2)==1 && k>97
        LPool1(k)= 5*(97-1)/2+6*(k-97)/2; 
    elseif k<98
        LPool1(k)=(k/2)*2+(k/2-1)*3; 
    elseif k>97
        LPool1(k)=(96/2)*2+(96/2)*3 +((k-96)/2)*3+((k-96)/2-1)*3;  
    end             
end

y0=[1.1,40];
LA=zeros([],1);
yA=zeros([],2);
% Solving ODE system with Loop
for k=1:numel(LPool)-1
    hap=hapA(k);
    if k<97
        [L,y]=ode45(@myODE,LPool(k:k+1),y0);
        y0=y(end,:);
        y1(1)=0.47;
        y1(2)=y0(2); 
    else
        [L,y]=ode45(@myODE,LPool(k:k+1),y1,[],hap); 
        y1=y(end,:); 
    end
    LA=[LA;L];
    yA=[yA;y];
    hapA=[hapA;hap];
end
mscA=yA(:,1)*10;
ms=sum(mscA);
%%%%%%%%%% OPTIMIZATION %%%%%%%%%%%%%%%%%%%%%%
hap1=[0.017;0.0197]; 
hap2=[0.0163;0.0192];
hapA0=[repmat(hap1,n1(1),1);
repmat(hap2,n1(2),1)];  % Confirmed heat transfer parameter between air and sheet
[hapA] = patternsearch(@objfun,hapA0);
%%%%%%%%%% FUNCTIONS %%%%%%%%%%%%%%%%%%%%%

function dy = myODE(L,y,hap)
u=y(1); 
Tp=y(2);
dudL = myODE1(L,u,Tp,hap); 
dTpdL = myODE2(L,u,Tp,hap);
dy = [dudL;dTpdL];
end

function dudL = myODE1(~,~,Tp,~) 
dudL=-hap/(Tp+273.15);  % u is in kg water/kg fiber
end

function dTpdL=myODE2(~,u,Tp,hap)
dTpdL=((1+0.955*u)*(100-Tp)+hap*(100-Tp))/(hap*u);                                                                                           
end

function f=objfun(ms)
msC=20.97;
f=msC-ms;
end    

[Lenilein]

Yes you are right, sorry for that, I will remember next time!

Answer 1

[Cris LaPierre]

You define hapA in line 9:

hapA=[];

But then immediately try to index into it. But you just defined it as empty, so you get the error.

hap=hapA(k);

This page shows you how to determe why you are getting error messages.

[Lenilein]

Thanks Cris, it is clear now why I got this error.

Ideally I would like to not define values for hapA initially as this the output I expect from the optimization.

I also tried removing the declaration of hap from the for loop and passing hap as input variable for the function myODE (see below) but that leads (quite logically) to a problem when solving the ODE:

>> Simplified

Not enough input arguments.

Error in Simplified>myODE (line 33)

dudL = myODE1(L,u,Tp,hap); dTpdL = myODE2(L,u,Tp,hap);

Error in Simplified>@(L,y)myODE(L,y)

Error in odearguments (line 90)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in Simplified (line 16)

if k<97;[L,y]=ode45(@(L,y) myODE(L,y),LPool(k:k+1),y0); y0=y(end,:); y1(1)=0.47; y1(2)=y0(2); else; [L,y]=ode45(@(L,y)

myODE(L,y),LPool(k:k+1),y1); y1=y(end,:); end

s=127;
LPool1=zeros(s,1);
n1=[48,15];                     % Number of ylinders in pre- and postdrying sections

% Defining the matrix LPool
for k=1:numel(LPool1)           
        if rem(k,2)==1 && k<98; LPool1(k)= 5*(k-1)/2; elseif rem(k,2)==1 && k>97; LPool1(k)= 5*(97-1)/2+6*(k-97)/2; elseif k<98; LPool1(k)=(k/2)*2+(k/2-1)*3; elseif k>97; LPool1(k)=(96/2)*2+(96/2)*3 +((k-96)/2)*3+((k-96)/2-1)*3;  end             
end
LPool=LPool1.';

y0=[1.1,40];                    
LA=zeros([],1); yA=zeros([],2);

% Solving ODE system with Loop
for k=1:numel(LPool)-1  
if k<97;[L,y]=ode45(@(L,y) myODE(L,y),LPool(k:k+1),y0); y0=y(end,:); y1(1)=0.47; y1(2)=y0(2); else; [L,y]=ode45(@(L,y) myODE(L,y),LPool(k:k+1),y1); y1=y(end,:); end

LA=[LA;L];yA=[yA;y];hapA=[hapA;hap];
end

mscA=yA(:,1)*10;
ms=sum(mscA);

%%%%%%%%%% OPTIMIZATION %%%%%%%%%%%%%%%%%%%%%%

hap1=[0.017;0.0197]; hap2=[0.0163;0.0192];
hapA0=[repmat(hap1,n1(1),1);repmat(hap2,n1(2),1)];  % Confirmed heat transfer parameter between air and sheet
[hapsol] = patternsearch(@objfun,hapA0);

%%%%%%%%%% FUNCTIONS %%%%%%%%%%%%%%%%%%%%%
function dy = myODE(L,y,hap)
u=y(1); Tp=y(2);
dudL = myODE1(L,u,Tp,hap); dTpdL = myODE2(L,u,Tp,hap);
dy = [dudL;dTpdL];
end

function dudL = myODE1(~,~,Tp,hap)    
dudL=-hap/(Tp+273.15);  % u is in kg water/kg fiber
end

function dTpdL=myODE2(~,u,Tp,hap)
dTpdL=((1+0.955*u)*(100-Tp)+hap*(100-Tp))/(hap*u);                                                                                           
end

function f=objfun(ms)
msC=20.97;
f=msC-ms;
end    

So my question is rather regarding the way to general procedure with the optimization of a parameter fixed for each iteration loop of the ODE solving procedure than regarding the error message that do make sense to me...

[Cris LaPierre]

Let's start small. Try to get the following subset of your code working as you want first.

n1=[48,15];       % Number of ylinders in pre- and postdrying sections

%%%%%%%%%% OPTIMIZATION %%%%%%%%%%%%%%%%%%%%%%
hap1=[0.017;0.0197]; 
hap2=[0.0163;0.0192];
hapA0=[repmat(hap1,n1(1),1);
repmat(hap2,n1(2),1)];  % Confirmed heat transfer parameter between air and sheet
[hapA] = patternsearch(@objfun,hapA0);
%%%%%%%%%% FUNCTIONS %%%%%%%%%%%%%%%%%%%%%

function f=objfun(ms)
msC=20.97;
f=msC-ms;
end  

It currently gives the following error at [hapA] = patternsearch(@objfun,hapA0);

Error using poptimfcnchk (line 32)
Your objective function must return a scalar value.


Error in patternsearch (line 262)
[Iterate,OUTPUT.funccount] = poptimfcnchk(FUN,nonlcon,X0,Iterate, ...