试用软件

Why my solver still simulating a model, even there is no input value given to the S-function?

1 次查看（过去 30 天）

显示更早的评论

fiza fiza 2017-1-9

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/319564-why-my-solver-still-simulating-a-model-even-there-is-no-input-value-given-to-the-s-function

编辑： fiza fiza 2017-1-9

在 MATLAB Online 中打开

Hi..

I try to simulate a dynamic system with 9 state derivatives (4 input and 9 output) using level-2 S-function in SIMULINK. This simulation works well.

But, when I delete one/more input value(s) that been given to this model, my solver (ode23s) still solve the model and the simulation run. Why could this happen? Could anybody tell me, since I really confuse right now.. :-(

Here is my S-function:

%%%%%%%%%%%%%%%%%%%%%%%%%%%

function level2_blood(block)

setup(block);

function setup(block)

% Register number of ports

block.NumInputPorts = 4;

block.NumOutputPorts = 9;

% Setup port properties to be inherited or dynamic

block.SetPreCompInpPortInfoToDynamic;

block.SetPreCompOutPortInfoToDynamic;

% Override input port properties

block.InputPort(1).Dimensions = 1;

block.InputPort(1).DatatypeID = 0; % double

block.InputPort(1).Complexity = 'Real';

block.InputPort(1).DirectFeedthrough = true;

block.InputPort(2).Dimensions = 1;

block.InputPort(2).DatatypeID = 0; % double

block.InputPort(2).Complexity = 'Real';

block.InputPort(2).DirectFeedthrough = true;

block.InputPort(3).Dimensions = 1;

block.InputPort(3).DatatypeID = 0; % double

block.InputPort(3).Complexity = 'Real';

block.InputPort(3).DirectFeedthrough = true;

block.InputPort(4).Dimensions = 1;

block.InputPort(4).DatatypeID = 0; % double

block.InputPort(4).Complexity = 'Real';

block.InputPort(4).DirectFeedthrough = true;

% Override output port properties

block.OutputPort(1).Dimensions = 1;

block.OutputPort(1).DatatypeID = 0; % double

block.OutputPort(1).Complexity = 'Real';

% Override output port properties

block.OutputPort(2).Dimensions = 1;

block.OutputPort(2).DatatypeID = 0; % double

block.OutputPort(2).Complexity = 'Real';

% Override output port properties

block.OutputPort(3).Dimensions = 1;

block.OutputPort(3).DatatypeID = 0; % double

block.OutputPort(3).Complexity = 'Real';

% Override output port properties

block.OutputPort(4).Dimensions = 1;

block.OutputPort(4).DatatypeID = 0; % double

block.OutputPort(4).Complexity = 'Real';

% Override output port properties

block.OutputPort(5).Dimensions = 1;

block.OutputPort(5).DatatypeID = 0; % double

block.OutputPort(5).Complexity = 'Real';

% Override output port properties

block.OutputPort(6).Dimensions = 1;

block.OutputPort(6).DatatypeID = 0; % double

block.OutputPort(6).Complexity = 'Real';

% Override output port properties

block.OutputPort(7).Dimensions = 1;

block.OutputPort(7).DatatypeID = 0; % double

block.OutputPort(7).Complexity = 'Real';

% Override output port properties

block.OutputPort(8).Dimensions = 1;

block.OutputPort(8).DatatypeID = 0; % double

block.OutputPort(8).Complexity = 'Real';

% Override output port properties

block.OutputPort(9).Dimensions = 1;

block.OutputPort(9).DatatypeID = 0; % double

block.OutputPort(9).Complexity = 'Real';

% Register parameters

block.NumDialogPrms = 9;

% Set up the continuous states.

block.NumContStates = 9;

block.SampleTimes = [0 0];

block.SimStateCompliance = 'DefaultSimState';

block.RegBlockMethod('SetInputPortSamplingMode',@SetInputPortSamplingMode);

block.RegBlockMethod('InitializeConditions', @InitializeConditions);

block.RegBlockMethod('Outputs', @Outputs); % Required

block.RegBlockMethod('Derivatives', @Derivatives);

block.RegBlockMethod('Terminate', @Terminate); % Required

function SetInputPortSamplingMode(block, idx, fd)

 block.InputPort(idx).SamplingMode = fd;
 block.InputPort(idx).SamplingMode = fd;
 block.InputPort(idx).SamplingMode = fd;
 block.InputPort(idx).SamplingMode = fd;

block.OutputPort(1).SamplingMode = fd;

block.OutputPort(2).SamplingMode = fd;

block.OutputPort(3).SamplingMode = fd;

block.OutputPort(4).SamplingMode = fd;

block.OutputPort(5).SamplingMode = fd;

block.OutputPort(6).SamplingMode = fd;

block.OutputPort(7).SamplingMode = fd;

block.OutputPort(8).SamplingMode = fd;

block.OutputPort(9).SamplingMode = fd;

function InitializeConditions(block)

%initialize continuous states

carb_init= block.DialogPrm(1).Data;

pco2pl_init= block.DialogPrm(2).Data;

pco2rbc_init= block.DialogPrm(3).Data;

hco3pl_init= block.DialogPrm(4).Data;

hco3rbc_init= block.DialogPrm(5).Data;

hpl_init= block.DialogPrm(6).Data;

hrbc_init= block.DialogPrm(7).Data;

po2_init= block.DialogPrm(8).Data;

phvirt_init= block.DialogPrm(9).Data;

block.ContStates.Data(1) = carb_init;

block.ContStates.Data(2) = pco2pl_init;

block.ContStates.Data(3) = pco2rbc_init;

block.ContStates.Data(4) = hco3pl_init;

block.ContStates.Data(5) = hco3rbc_init;

block.ContStates.Data(6)= hpl_init;

block.ContStates.Data(7) = hrbc_init;

block.ContStates.Data(8) = po2_init;

block.ContStates.Data(9) = phvirt_init;

%end InitializeConditions

%% Outputs:

function Outputs(block)

block.OutputPort(1).Data = block.ContStates.Data(1);

block.OutputPort(2).Data = block.ContStates.Data(2);

block.OutputPort(3).Data = block.ContStates.Data(3);

block.OutputPort(4).Data = block.ContStates.Data(4);

block.OutputPort(5).Data = block.ContStates.Data(5);

block.OutputPort(6).Data = block.ContStates.Data(6);

block.OutputPort(7).Data = block.ContStates.Data(7);

block.OutputPort(8).Data = block.ContStates.Data(8);

block.OutputPort(9).Data = block.ContStates.Data(9);

%end Outputs

%% Derivatives:

function Derivatives(block) Xa(1)= block.ContStates.Data(1);

Xa(2)= block.ContStates.Data(2);

Xa(3)= block.ContStates.Data(3);

Xa(4)= block.ContStates.Data(4);

Xa(5)= block.ContStates.Data(5);

Xa(6)= block.ContStates.Data(6);

Xa(7)= block.ContStates.Data(7);

Xa(8)= block.ContStates.Data(8);

Xa(9)= block.ContStates.Data(9);

%model parameter values

Qb= block.InputPort(1).Data;

Vb= block.InputPort(2).Data;

PO2_g= block.InputPort(3).Data;

pco2_g= block.InputPort(4).Data;

alpha_co2= 3e-3;%co2 solubility

co2_pl= alpha_co2 * Xa(2);

co2_rbc= alpha_co2 * Xa(3);

beta_pl= 6e-3; %plasma buffer capacity

beta_rbc= 57.7e-3; % RBC buffer capacity

cat= 13000; %CA catalysis factor

K= 5.5e-4; %CA dissociation equilibrium constant

Ku= 0.12; %CO2 hydratation reaction forward constant rate, per second

Kv= 89; %co2 hydratation reaction reverse rate constant

Ka= 5000; %CO2-Hb forward constant

Kc= 2.4e-5; %carbamate ionisation constant

Kzo= 8.4e-9; %oxygenated Hb amino group ionisation constant

Kzr= 7.2e-8; %reduced Hb amino groupionisation constant

T_rbc= 0.001; %tau_rbc, half time of red cell membrane diffusion

T_hco3=0.2; %tau_hco3, half time of rbc membrane CL shift

hct= 0.3; % hematocrit

cap= 0.201; %O2 capacity of hemoglobin (per ml per blood)

vrbc= Vb*hct;

vpl= Vb*(1-hct);

Q_rbc= Qb*hct;

Q_pl= Qb* (1-hct);

a1= -8532;

a2= 2121;

a3= -67.07;

a4= 936000;

a5= -31350;

a6= 2396;

a7= -67.10;

T= 28;

ph_virt= Xa(9)

syms x

S= ((a1*x) + (a2*(x^2)) + (a3*(x^3)) + (x^4))/ (a4 + a5*x + (a6*(x^2)) + (a7*(x^3)) + (x^4));

diff_S= jacobian(S,x);

diff_simple= simplify(diff_S);

eqn= Xa(8) * (10^(0.024*(37-T) + 0.4*(ph_virt - 7.4)+ 0.06*(log(40/Xa(3)))));

int_subs= subs(diff_simple,x,eqn);

h= double(int_subs);

syms P

X= P * (10^(0.024*(37-T) + 0.4*(ph_virt - 7.4)+ 0.06*(log(40/Xa(3)))));

diff_X= jacobian(X,P);

diff_simpleP= simplify (diff_X);

int_subsP= subs(diff_simpleP,P,Xa(8));

g= double(int_subsP);

m= h*g;

x= Xa(8) * (10^(0.024*(37-T) + 0.4*(ph_virt - 7.4)+ 0.06*(log(40/Xa(3)))));

S= ((a1*x) + (a2*(x^2)) + (a3*(x^3)) + (x^4))/ (a4 + a5*x + (a6*(x^2)) + (a7*(x^3)) + (x^4));

dSdt= gradient (S);

r= (0.058*Xa(9)-0.437)*S - (0.529*Xa(9)) +4.6;

Rhco3_pl= (-Ku*alpha_co2* Xa(2))+ ((Kv/K)*Xa(6)*Xa(4));

Rhco3_rbc= cat*((-Ku *alpha_co2 *Xa(3))+ ((Kv/K)*Xa(7) *Xa(5)));

Dco2_rbc=((0.693*alpha_co2)/T_rbc)*((vrbc*vpl)/(vrbc + vpl));

Dhco3_rbc= (0.693/T_hco3)*((vrbc*vpl)/(vrbc + vpl));

carb_in= 0.00235; %from misgeld thesis

co2_plin= 0.00138; %from misgeld thesis

co2_rbcin= 0.00138; %from misgeld thesis

hco3_plin= 0.0263;%from misgeld thesis

hco3_rbcin= 0.0182; %from misgeld thesis

H_plin= 42.3e-9; %from misgeld thesis

H_rbcin= 61e-9;%from misgeld thesis

Hb= 20.58e-3;

Hb_rbc= 20.7e-3;

DO2_m= 11.291e-6; % total membrane O2 diffusing capacity of the oxygenator, Hill et al. 1973

Dco2_m= 414.64e-6; % total membrane co2 diffusing capacity of the oxygenator, Hill et al. 1973

alpha_O2= 1.35e-6;

cap_b= hct*Hb_rbc;

O2_bloodin= 0.0068;

O2_blood= (1.35e-6 * Xa(8))+ (cap_b * S);

TpH= 100000;

dcarbdt= (Q_rbc*(carb_in - Xa(1))+ Ka*co2_rbc vrbc(Hb - Xa(1))* (( (Kzo* S /(Kzo + Xa(7))) + ((Kzr*(1-S)))/ (Kzr + Xa(7))))- (vrbc*(Ka* Xa(1) * Xa(7))/Kc))/vrbc;

dpco2pldt= (Q_pl*(co2_plin-co2_pl)+ Dco2_m*(pco2_g - Xa(2)) + Dco2_rbc*(Xa(3) - Xa(2))+ vpl*Rhco3_pl)/ vpl*alpha_co2;

dpco2rbcdt= (Q_rbc*(co2_rbcin - co2_rbc) + Dco2_rbc*(Xa(2)-Xa(3))+ vrbc*Rhco3_rbc - vrbc*dcarbdt)/ vrbc*alpha_co2;

dhco3_pldt= (Q_pl*(hco3_plin -Xa(4)) - Dhco3_rbc*(Xa(4) - ((Xa(5))/r))- vpl*Rhco3_pl)/ vpl;

dhco3_rbcdt= (Q_rbc*(hco3_rbcin - Xa(5))+ Dhco3_rbc*(Xa(4) - ((Xa(5))/r))- (vrbc*Rhco3_rbc)) / vrbc;

dH_pldt= (Q_pl*(H_plin - Xa(6))- (vpl*(2.303/beta_pl)*(Xa(6))*Rhco3_pl))/ vpl; dH_rbcdt= (Q_rbc*(H_rbcin - Xa(7))- (vrbc*(2.303/beta_rbc)*Xa(7))*(-Rhco3_rbc + 1.5* dcarbdt - 0.6*cap*dSdt))/vrbc;

dPO2dt= (Qb*(O2_bloodin - O2_blood)+ (DO2_m * (PO2_g - Xa(8))))./ (Vb * (alpha_O2 + (cap_b * m)));

dphvirtdt= (-Xa(9)-log(Xa(7)*((0.058*Xa(9)-0.437)*S - (0.529*Xa(9)) +4.6)))/TpH;

% matrix composite

block.Derivatives.Data(1) = dcarbdt;

block.Derivatives.Data(2)= dpco2pldt;

block.Derivatives.Data(3)= dpco2rbcdt;

block.Derivatives.Data(4)= dhco3_pldt;

block.Derivatives.Data(5)= dhco3_rbcdt;

block.Derivatives.Data(6)= dH_pldt;

block.Derivatives.Data(7)= dH_rbcdt;

block.Derivatives.Data(8)= dPO2dt;

block.Derivatives.Data(9)= dphvirtdt;

%end Derivatives

%% Terminate:

function Terminate(block)

%end Terminate

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Thank you..

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

回答（0 个）

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

类别

Physical Modeling Simscape Multibody Multibody Modeling Assembly

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

标签

Community Treasure Hunt

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

Translated by

试用软件