Why my solver still simulating a model, even there is no input value given to the S-function?
1 回表示 (過去 30 日間)
古いコメントを表示
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 件のコメント
回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Assembly についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!