Error when using if-else statement in MATLAB function block. How can I fix this error described below?

2024 5 5
4 个回答
更新时间：2024 5 14
9 次查看（30 天）

0 个投票

function m = fcn(Q, T_sat, T_w, sigma, rho_v, h_fg, M, P_v, R, g, beta, x, nu_novec, alpha_novec, k_m, d_pore, psi, k_f, rporous_outer, height_porous)
q_w = Q/(2*pi*rporous_outer*height_porous); %heat input converted to heat flux
if T_w>=T_sat
m = (((2*sigma)/(2-sigma))*((rho_v*(h_fg)^2)/T_sat)*((M/(2*pi*R*T_sat))^(0.5))*(1- ((P_v)/(2*rho_v*h_fg))));
else
m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x);
%the power of 0.4 is what triggers the complex number
end
When I run this code, I get the following error, which is triggered by the line ' m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x);
I've checked if any of the values result in a negative number, but none of the terms are negative, nor their result. I've run this separately, and this the same line works fine and I get a positive results without an complex number or any error? Why would this occur, does this error have to do anything with the if-else statement used? Any thoughts would be much appreciated.
I think this issue is triggered by the use of ^0.4, which cannot be changed. Why would an error/complex result occur in simulink and not in MATLAB?
Many thanks

1 个评论
显示 -1更早的评论 隐藏 -1更早的评论

Sam Chak
Sam Chak 2024-5-6
Ran in:
Ran in:
Hi @Delika Weragoda
Could you please provide us with the values of the parameters? This will allow us to test the function 'fcn()' and compute 'm'.
m = fcn(1, 10, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
m = Inf
function m = fcn(Q, T_sat, T_w, sigma, rho_v, h_fg, M, P_v, R, g, beta, x, nu_novec, alpha_novec, k_m, d_pore, psi, k_f, rporous_outer, height_porous)
q_w = Q/(2*pi*rporous_outer*height_porous); %heat input converted to heat flux
if T_w>=T_sat
m = (((2*sigma)/(2-sigma))*((rho_v*(h_fg)^2)/T_sat)*((M/(2*pi*R*T_sat))^(0.5))*(1- ((P_v)/(2*rho_v*h_fg))));
else
m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x);
end
end

请先登录，再进行评论。

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

回答（4 个）

Walter Roberson
Walter Roberson 2024-5-6

0 个投票

m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x);
1 2 345 4 5 43 456 5 6 7 654 32 3 210 1 0
Simulink is calculating the expression
(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))
as being negative. When negative is raised to the fractional power 0.4 it results in a complex value.
My guess is that psi is outside the range -1 to +1
If you want a complex result, then code
m = (k_f*complex(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x);
If you do not want a complex result, you will need to protect against it, such as
if abs(psi) >= 1
m = nan;
elseif T_w>=T_sat
m = (((2*sigma)/(2-sigma))*((rho_v*(h_fg)^2)/T_sat)*((M/(2*pi*R*T_sat))^(0.5))*(1- ((P_v)/(2*rho_v*h_fg))));
else
m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x);
end

9 个评论
显示 7更早的评论 隐藏 7更早的评论

Delika Weragoda
Delika Weragoda 2024-5-6
移动：Sam Chak 2024-5-6
Hi @Sam Chak/ @Walter Roberson
@Walter Roberson none of the values are negative including psi, I am not sure whether its triggered by Q and T_w, but when I run the same line of code [i.e. m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x)], in the MATLAB function block without an if-else statement and the other condition, I don't get any of these error, works perfectly.
Please see the values of the parameters below,
Q - a signal taken from heat flow sensor block (constant in this case = 10)
T_sat = 307
T_w = taken from temperature sensor block (all positive)
sigma = 0.0002
rho_v = 7.1
h_fg M = 142e3
P_v = 9.15e4
R = 41.55
g = 9.81
beta = 0.002
x = 0.037
nu_novec = 3.2e-7
alpha_novec = 4.12e-8;
k_m = 3160
d_pore = 0.05
psi = 0.46
k_f = 0.075
rporous_outer = 0.0103
height_porous = 0.074
@Sulaymon Eshkabilov thank you for the suggestion, I've tried both does work.
Sam Chak
Sam Chak 2024-5-6
移动：Sam Chak 2024-5-6
Hi @Delika Weragoda
The value for 'h_fg' and sensor data 'T_w' are not provided.
Q = 10;
T_sat = 307;
T_w = 299:305; % <-- this one
sigma = 0.0002;
rho_v = 7.1;
h_fg = % <-- and this one
M = 142e3;
P_v = 9.15e4;
R = 41.55;
g = 9.81;
beta = 0.002;
x = 0.037;
nu_novec = 3.2e-7;
alpha_novec = 4.12e-8;
k_m = 3160;
d_pore = 0.05;
psi = 0.46;
k_f = 0.075;
rporous_outer = 0.0103;
height_porous = 0.074;
for j = numel(T_w)
m(j) = fcn(Q, T_sat, T_w(j), sigma, rho_v, h_fg, M, P_v, R, g, beta, x, nu_novec, alpha_novec, k_m, d_pore, psi, k_f, rporous_outer, height_porous);
end
stem(m), grid on
function m = fcn(Q, T_sat, T_w, sigma, rho_v, h_fg, M, P_v, R, g, beta, x, nu_novec, alpha_novec, k_m, d_pore, psi, k_f, rporous_outer, height_porous)
q_w = Q/(2*pi*rporous_outer*height_porous); %heat input converted to heat flux
if T_w>=T_sat
m = (((2*sigma)/(2-sigma))*((rho_v*(h_fg)^2)/T_sat)*((M/(2*pi*R*T_sat))^(0.5))*(1- ((P_v)/(2*rho_v*h_fg))));
else
m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x);
end
end
Delika Weragoda
Delika Weragoda 2024-5-6
Hi @Sam Chak Sorry, please see the values below.
Q = 10; %This is taken as an input from the model (see below)
T_sat = 307;
T_w = 299:305; % This is also taken as an input from the model (see below)
sigma = 0.0002;
rho_v = 7.1;
h_fg = 142e3
M = 0.2;
P_v = 9.15e4;
R = 41.55;
g = 9.81;
beta = 0.002;
x = 0.037;
nu_novec = 3.2e-7;
alpha_novec = 4.12e-8;
k_m = 3160;
d_pore = 0.05;
psi = 0.46;
k_f = 0.075;
rporous_outer = 0.0103;
height_porous = 0.074;
Please see the model for your reference
Sam Chak
Sam Chak 2024-5-6
Ran in:
Ran in:
@Delika Weragoda
Could you please verify if the m function returns the expected values? So far no errors.
%% Parameters
Q = 10;
T_sat = 307;
sigma = 0.0002;
rho_v = 7.1;
h_fg = 142e3;
M = 142e3;
P_v = 9.15e4;
R = 41.55;
g = 9.81;
beta = 0.002;
x = 0.037;
nu_novec = 3.2e-7;
alpha_novec = 4.12e-8;
k_m = 3160;
d_pore = 0.05;
psi = 0.46;
k_f = 0.075;
rporous_outer = 0.0103;
height_porous = 0.074;
%% Case 1a: 1 deg lesser than T_sat
T_w = 306;
m = fcn(Q, T_sat, T_w, sigma, rho_v, h_fg, M, P_v, R, g, beta, x, nu_novec, alpha_novec, k_m, d_pore, psi, k_f, rporous_outer, height_porous)
m = 22.4468
%% Case 1b: absolute zero
T_w = 0;
m = fcn(Q, T_sat, T_w, sigma, rho_v, h_fg, M, P_v, R, g, beta, x, nu_novec, alpha_novec, k_m, d_pore, psi, k_f, rporous_outer, height_porous)
m = 22.4468
%% Case 2a: 1 deg greater than T_sat
T_w = 308;
m = fcn(Q, T_sat, T_w, sigma, rho_v, h_fg, M, P_v, R, g, beta, x, nu_novec, alpha_novec, k_m, d_pore, psi, k_f, rporous_outer, height_porous)
m = 1.1852e+05
%% Case 2b: % speed of light
T_w = 3e8;
m = fcn(Q, T_sat, T_w, sigma, rho_v, h_fg, M, P_v, R, g, beta, x, nu_novec, alpha_novec, k_m, d_pore, psi, k_f, rporous_outer, height_porous)
m = 1.1852e+05
%% The m-function
function m = fcn(Q, T_sat, T_w, sigma, rho_v, h_fg, M, P_v, R, g, beta, x, nu_novec, alpha_novec, k_m, d_pore, psi, k_f, rporous_outer, height_porous)
q_w = Q/(2*pi*rporous_outer*height_porous); %heat input converted to heat flux
if T_w>=T_sat
m = (((2*sigma)/(2-sigma))*((rho_v*(h_fg)^2)/T_sat)*((M/(2*pi*R*T_sat))^(0.5))*(1- ((P_v)/(2*rho_v*h_fg))));
else
m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x);
end
end
Delika Weragoda
Delika Weragoda 2024-5-6
Hi @Sam Chak
Yes, these are the expected values. I do get them when I run the code in MATLAB, but not in Simulink.
I also tried simply using the the second m function (only) in the function block and it works fine (gets the results as above). However, I wanted to compare this with experimental results so I added an inport block to compare the two results and then the same error pops up, which is very weird.
I am thinking of upgrading MATLAB tomorrow, will let you know how it goes and if the error still persists.
Thanks
Delika
Delika Weragoda
Delika Weragoda 2024-5-7
编辑：Delika Weragoda 2024-5-7
@Sam Chak
I've managed to sort the issue, I've replaced a part of the function (in bold) using another variable (and running the replaced part on workspace) and the model works fine.
m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x);
m = m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(q^(0.4)))))/(x)
Thanks for your help.
Best regards
Delika
Sam Chak
Sam Chak 2024-5-7
Hi @Delika Weragoda,
Thank you for updating us. It's great to hear that you were able to resolve the issue through your own efforts. You might want to consider posting the full solution in the Answer section below and then accepting your self-discovered solution to close the issue.
This could be beneficial for other Simscape users who may encounter a similar issue in the future.
Delika Weragoda
Delika Weragoda 2024-5-9
编辑：Delika Weragoda 2024-5-9
@Sam Chak So, apparently, I've made a mistake in the parenthesis, that is why the model supposedly worked, but I managed to find the issue. I believe this error in complex results in triggered from the heat sensors' input that I take to calculate the value for m, which is dynamically changing.
The variable m, is a function of Q, which is the input taken from the heat flow rate sensor. This is what is causing the error, as this signal will be a power of 0.4. Is there a way to overcome this?
function m = fcn(k_f, g, beta, x, nu_novec, alpha_novec, k_m, d_pore, psi, Q, rporous_outer, height_porous)
q_w = Q/(2*pi*rporous_outer*height_porous); %(heat flux)
m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x);
Please see the model below for your reference
Sam Chak
Sam Chak 2024-5-9
Hi @Delika Weragoda
m = (k_f*(0.5*(((g*beta*q_w*x^4)/(nu_novec*alpha_novec*k_m))*(((d_pore^2)*psi^3/(180*(1-psi^2)))/x^2))^(0.4)))/(x);
I tend to divide long equations and conquer them later. This approach helps 100% prevent human mistakes such as mismatched parentheses. I suggest you rewrite the equation 'm' in multiple terms and then return here for a review.
%% parameters
a = ...;
b = ...;
c = ...;
d = ...;
%% terms
term1 = ...;
term2 = ...;
term3 = ...;
%% equation
m = term1 + (term2/term3)^0.4 ...;

请先登录，再进行评论。

Delika Weragoda
Delika Weragoda 2024-5-14

0 个投票

@Sam Chak
I found the issue. There is a negative heat flow generated due to the thermal mass that I am using in the model. Why would something like this occur? I've tried debugging the model by monitoring the heat flow at different sections of the model. Please see the results below.
Model:
heat transfer subsystem
Heat flow results
The four heat flow sensors shown in the figure above are used to monitor the heat flow at the four locations. The negative values at the first two time steps is what triggers the error.
I would really appreciate if you can give any feedback on how to rectify this please.
Many thanks

类别

在 帮助中心 和 File Exchange 中查找有关 Programmatic Model Editing 的更多信息

产品

版本

R2022b

标签

Community Treasure Hunt

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

Start Hunting!

Translated by