Optimization&ODE problem // Symptom: Index exceeds matrix dimensions (here we go again!)
1 view (last 30 days)
Show older comments
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
2 Comments
Cris LaPierre
on 8 Feb 2019
Edited: Cris LaPierre
on 8 Feb 2019
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
Accepted Answer
Cris LaPierre
on 8 Feb 2019
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.
2 Comments
Cris LaPierre
on 8 Feb 2019
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, ...
More Answers (0)
See Also
Categories
Find more on Introduction to Installation and Licensing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!