Optimization&ODE problem // Symptom: Index exceeds matrix dimensions (here we go again!)
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
Lenilein
el 8 de Feb. de 2019
Comentada: Cris LaPierre
el 8 de Feb. de 2019
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 comentarios
Cris LaPierre
el 8 de Feb. de 2019
Editada: Cris LaPierre
el 8 de Feb. de 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
Respuesta aceptada
Cris LaPierre
el 8 de Feb. de 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 comentarios
Cris LaPierre
el 8 de Feb. de 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, ...
Más respuestas (0)
Ver también
Categorías
Más información sobre Introduction to Installation and Licensing en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!