Error in ode15s: incompatible sizes arrays, while no errors in function

2 次查看（过去 30 天）
Jort Puiman
Jort Puiman 2023-10-2
评论： Jort Puiman 2023-10-3
Ran in:
Hi everyone,
I have script where I try to calculate the position of multiple particles over time. This is not a problem by itself, but the velocity of the particles is dependent on the sum of the particles that are present at each position. After a while, I managed to make a function that at least runs, but now I have an error in ode15s itself. I get similar errors when using other solvers (23s, 45).
The function file is a bit more than 100 lines, but the most important part is at the bottom, where I have my mass balances. Since I want to calculate everything for multiple particle sizes, I give MatLab a 100*4 matrix. MatLab makes it a 400*1 matrix, which I redistribute to the original 100*4 size for every calculation. At the end, I change it back to a 400*1 matrix.
clear
clc
method = 'distribution';
%% parameters
% time
tstep =100;
tend = tstep*100;
t = 0:tstep:tend;
% feed properties
rhow = 1000; % Density water
rhod = 1200; % Preliminary guess density particles
rhop = 1100; % Preliminary guess density protein
etaL = 0.001; % Preliminary guess viscosity continuous phase
g = 9.81; % Gravity
cmax = 0.6; % Maximal volumetric fraction
cini = 0.2; % Initial volumetric fraction
cprot = 0.55; % Preliminary guess protein fraction in algae
rho_in = cini*rhod + (1-cini)*rhow;
% centrifuge properties
r = 0.076; % Centrifugation radius
z = 0.5; % Centrifugation height (m)
Nr = 100; % Parts over radius
Nz = 100; % Parts over height
Ny = 100; % Parts over perpendicular to radius
dr = r./Nr; % Size of radius part (m)
dz = z./Nz; % Size of height part (m)
Vt= pi*r^2*z; % Total volume (m^3)
dy = zeros(length(linspace(0,360,Ny)),1);
Vcell = zeros(Ny,1);
mprot_ini = Vt*(1-cini)*cprot*rhop;
for i = 1:length(dy) % Calculation for size of parts perpendicular to radius
if i>1
dy(i) = pi*0.5*dr*i/Ny;
Vcell(i) = pi*(i*dr)^2*dz/Ny-pi*((i-1)*dr)^2*dz/Ny;
else
Vcell(i) = pi*dr^2*dz/Ny;
end
end
% run parameters
f = 40000/60; % Rotation frequency (1/s)
gc = (2*pi*f)^2*r; % Centrifugal force (m/s^2)
% energetic parameters
Volt = 230;
Ampere = 15;
% economical values
Cp = 5; % Price of protein (€/kg)
Ce = 0.39; % price of energy (€/kWh)
%% Choose method for model
switch method
case 'one size'
dp = 0.5e-6; % Particle size
c = zeros(Nr,1);
c(:) = cini; % Initial condition
case 'distribution'
dp1 = 0.15*10^-6; % Particle size distribution
dp2 = 0.25*10^-6;
dp3 = 0.35*10^-6;
dp4 = 0.45*10^-6;
dp = [dp1 dp2 dp3 dp4];
distribution = [0.1 0.3 0.4 0.2]; % Distribution percentage
c = ones(Nr,1);
cini = cini .* distribution;
c = c.*cini; % Initial conditions
result = zeros(length(t), length(c), length(dp));
end
%% Solver
tic % Measure time of calculation
tolerance = 1e-6;
options = odeset('RelTol',tolerance,'AbsTol',tolerance,'NonNegative',1:100);
% ODE
[t,C] = ode15s((@(t,c)odebatch31(t,c, Nr,dr,cmax,rhod,rhow,dp,g,f,method)),t,c,options);
Re difference too big
Arrays have incompatible sizes for this operation.

Error in ode15s (line 553)
rhs = hinvGak*ode(tnew,ynew) - Mtnew*(psi+difkp1);
for i = 0:length(dp)-1
result(:,:,i+1) = C(:,Nr*i+1:Nr*i+Nr);
end
toc % Measure time of calculation
sum_result = sum(result, 3); % Add distributions together to get local concentration
% Mass distribution loop
for j = size(sum_result,2):-1:1
for i = 1:size(sum_result,1)
if sum_result(i,j)>cmax
excess(i,j) = sum_result(i,j)-cmax;
sum_result(i,j) = cmax;
sum_result(i,j-1) = sum_result(i,j-1)+excess(i,j);
end
end
end
cmax_reached = zeros(length(t),1);
for i = 1:length(t)
[cmax, max_col] = max(sum_result(i,:));
cmax_reached(i) = max_col;
if cmax_reached(i) == 1
cmax_reached(i) = 100;
end
end
%% Control on c
Mbegin = sum(sum_result(1,:))*rhod + (1-sum(sum_result(1,:)))*rhow;
Mend = sum(sum_result(end,:))*rhod + (1-sum(sum_result(end,:)))*rhow;
Mbal = sum( Mend) - sum( Mbegin );
if Mbegin-Mend>1e-3 || Mbegin-Mend < -1e-3
fprintf('Unequal mass at beginning and end \n')
fprintf('Mass balance = %f kg/m3 \n',Mbal)
end
%% figure for check
figure(1)
xaxis = linspace(0,r*100,Nr);
yaxis = t/60;
surf(xaxis,yaxis,sum_result)
xlabel('r [cm]')
ylabel('time [min]')
zlabel('volumetric fraction [-]')
colorbar
%% Protein concentration and figures for checking
c_water = 1-sum_result;
c_prot = c_water*cprot;
pellet_thickness = (Nr-cmax_reached)*dr;
polynomal = polyfit(t, pellet_thickness, 3);
polyval = polyval(polynomal,t);
pos = 0.2*Nr; % choose percentage of time
tpos = pos*tstep;
figure(2)
plot(xaxis,c_prot(pos,:), xaxis,sum_result(pos,:))
title(['concentration through centrifuge at ', num2str(tpos/60,'%.1f'),' min'])
legend('protein','cell debris')
xlabel('r [cm]')
ylabel('fraction[-]')
% figure(3)
% plot(t, pellet_thickness,t,polyval, 'or')
% legend('pellet size', 'polynomal')
% xlabel('time [s])')
% ylabel('pellet size [m]')
%% calculation for harvestable proteins
% Let's say that proteins can be harvested when the concentration of cell
% cebris is lower than 0.1. Now, we can calculate the protein in
% each cell, translate it to volumetric cells and therefore to a mass, and
% lastly calculate the total mass that can be harvested.
c_prot_harvest = c_water*cprot;
c_prot_harvest(sum_result>0.1)=0;
m_prot_harvest = c_prot_harvest.*Vcell'*rhop*Ny*Nz;
figure(4)
yyaxis left
plot(xaxis,m_prot_harvest(pos,:))
ylabel('mass [kg]')
yyaxis right
plot(xaxis,c_prot_harvest(pos,:))
title(['concentration through centrifuge at ', num2str(tpos/60,'%.1f'),' min'])
legend('protein mass','protein fraction')
xlabel('r [cm]')
ylabel('fraction [-]')
harvest = zeros(length(t),1); slope = zeros(length(t)-1,1);
for i = 1:length(t)
harvest(i) = sum(m_prot_harvest(i,:));
end
flag = false;
for i = 20:length(t)
slope(i) = (harvest(i)-harvest(i-1))/tstep;
if all(harvest(i)-harvest(i-1)<0.00005) && ~flag
disp(['Optimal centrifugal time is ', num2str((i*tstep)/60, '%.1f'), ' min'])
disp(['Protein harvest at this time is ', num2str(harvest(i), '%.2f'),' kg'])
disp(['The yield on protein is ', num2str(((harvest(i)/mprot_ini)*100), '%.2f'),'%'])
flag = true;
revenue(i) = harvest(i)*Cp;
energy = (rho_in*(2*pi*f*r)^2*Vt./t); % J
costs = energy/1000*Ce;
profit(i) = revenue(i)-costs;
disp(['The profit is €',num2str(profit(i), '%.2f')])
end
end
revenue = harvest*Cp;
energy = rho_in*(2*pi*f*r)^2*Vt./t; % J
costs = energy/1000*Ce;
profit = revenue-costs;
figure(5)
plot(t/60,(harvest/mprot_ini))
xlabel('time [min]')
ylabel('mass [kg]')
title('Protein harvest over time')
yyaxis right
plot(t/60, profit)
ylabel('profit [€]')
function [dcdt] = odebatch31(~,c, Nr,dr,cmax,rhod,rhow,dp,g,f,method)
%% zeros matrices
v = zeros(Nr,length(dp));
excess = zeros(Nr,1);
dcdt = zeros(Nr,length(dp));
Re = zeros(Nr,length(dp));
eta = zeros(Nr,1);
vgs = zeros(Nr,length(dp));
Ar = zeros(Nr,length(dp));
n_check = zeros(Nr,length(dp));
n = zeros(Nr,length(dp));
v_check = zeros(Nr,length(dp));
Re_check = zeros(Nr,length(dp));
csum = zeros(Nr,1);
c = reshape(c, Nr, length(dp));
r = zeros(Nr,1);
gc = zeros(Nr,1);
%% Mass distribution loop
switch method
case 'one size'
for i=Nr:-1:2
if c(i) >= cmax
excess(i) = c(i)-cmax;
c(i) = cmax;
c(i-1) = c(i-1)+excess(i);
end
end
% Move location of wall condition
saturated = zeros(Nr, 1);
for i = 1:Nr
if c(i) == cmax
saturated(i) = 1;
end
end
wall_coords = find(saturated == 1);
if isempty(wall_coords)
wall_condition = Nr;
else
wall_condition = wall_coords(1)-1;
end
case 'distribution'
for i=Nr:-1:2
clocal(i,:) = c(i,:);
clocal_sum(i) = sum(clocal(i));
if clocal_sum(i) >= cmax
distribution(i,:) = clocal(i,:)./clocal(i)*cmax;
c(i,:) = distribution(i,:);
excess(i,:) = clocal(i,:)-distribution(i,:);
c(i-1,:) = c(i-1,:)+excess(i,:);
end
end
% Move location of wall condition
saturated = zeros(Nr, 1);
for i = 1:Nr
if c(i) == cmax
saturated(i) = 1;
end
end
wall_coords = find(saturated == 1);
if isempty(wall_coords)
wall_condition = Nr;
else
wall_condition = wall_coords(1)-1;
end
end
%% Main loop
for i=1:wall_condition
%% find values for velocity and Reynolds check
for j = 1:length(dp)
dparticle(j) = dp(1,j);
eta(i) = exp(2.5*(sum(c(i,j)))./(1-sum(c(i,j))/cmax));
r(i) = r(i)+dr;
gc(i) = (2*pi*f)^2*r(i);
eta(i) = 0.001;
vgs(i,j) = ((rhod-rhow).*dparticle(j).^2.*gc(i))./(18*eta(i));
Ar(i,j) = (dparticle(j).^3*rhow*(rhod-rhow)*g)./eta(i).^3;
n_check(i,j) = (4.8+2.4*0.015.*Ar(i,j).^0.5)/(0.015.*Ar(i,j).^0.5+1);
v_check(i,j) = (1-c(i,j)).^n_check(i,j).*vgs(i,j);
Re_check(i,j) = (rhod*v_check(i,j).*dparticle(j))./eta(i);
n(i,j) = (4.8-2.4*Re_check(i,j).^(3/4)*0.015)/(1+Re_check(i,j).^(3/4)*0.015);
v(i,j) = (1-sum(c(i,j))).^n(i,j).*vgs(i,j);
Re(i,j) = (rhod.*v(i,j).*dparticle(j))./eta(i);
if Re(i,j)/Re_check(i,j)<0.95 || Re(i,j)/Re_check(i,j)>1.05
fprintf('Re difference too big \n')
return
end
if Re(i,j)>1
fprintf('Re too big for stokes \n')
return
end
%% ODE calculation
if i > 1
dcdt(i,j) = (1/dr).*((c(i-1,j).*v(i-1,j))-c(i,j).*v(i,j)); % general ode
else
dcdt(i,j)=(1/dr).*-(c(i,j).*v(i,j)); % initial condition
end
if i == wall_condition % wall condition
dcdt(i,j) = (1/dr).*(c(i-1,j).*v(i-1,j));
end
csum(i) = sum(c(i,:));
end
end
dcdt = dcdt(:);
end
When running the code, I do not get any errors, only after everything is run through the function and the ode is calculated, I get the following error:
Arrays have incompatible sizes for this operation.
Error in ode15s (line 553)
rhs = hinvGak*ode(tnew,ynew) - Mtnew*(psi+difkp1);
Error in centrifugation31 (line 83)
[t,C] = ode15s((@(t,c)odebatch31(t,c, Nr,dr,cmax,rhod,rhow,dp,g,f,method)),t,c,options);
I attached my script and function, so that you can maybe understand it a bit better.
I am looking forward to a reply, and a big thanks to anyone that wants to take a look!

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2023-10-2
编辑：Torsten 2023-10-2
if Re(i,j)/Re_check(i,j)<0.95 || Re(i,j)/Re_check(i,j)>1.05
fprintf('Re difference too big \n')
return
end
if Re(i,j)>1
fprintf('Re too big for stokes \n')
return
end
You cannot simply return from odebatch31 ; the solver will crash in this case because it expects at least any output for dcdt of size 400x1.
When I comment these lines of your code, I get an error in this if-statement:
if clocal_sum(i) >= cmax
distribution(i,:) = clocal(i,:)./clocal(i)*cmax;
c(i,:) = distribution(i,:);
excess(i,:) = clocal(i,:)-distribution(i,:);
c(i-1,:) = c(i-1,:)+excess(i,:);
end
The line
distribution(i,:) = clocal(i,:)./clocal(i)*cmax;
looks incorrect.
  1 个评论
Jort Puiman
Jort Puiman 2023-10-3
Thank you so much! I completely overlooked that return would have this effect.
For the other comment, I am aware that my 'distribution' method is not working properly. I have a post-ode fix, but this is not how it should be. I will look at how to fix this.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2023b

Community Treasure Hunt

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

Start Hunting!

Translated by