2

我写了一个解决颂歌的代码（细胞随着时间的推移而生长，直到达到一定的产品浓度）。因为当 y(3) 的值大于常数 (Pmx) 时我的模型会发生变化，所以我实现了 odeset 事件函数。如果事件发生取决于函数中给出的值。

对于调试，我想设置一个 if 语句来停止函数，并且只使用第一个 ode23 函数结果来绘制（尚未实现）（见代码）

我认为这也是最终排名的问题，不是吗？提前谢谢

function [T,Y,TE,YE,IE, Y1, T1] = SimplePerfusionMode(tf,F)

%% DESCRIPTION
% INPUT:
%   tf: final time [h]
%    F: Feed stream [Lh^-1]
% OUTPUT
%   T:  processing time [h]
%   Y:  states functions (X, S, P)
%   TE: time at which event occurs [h]
%   YE: Solution values corresponding to TE (g/l)
%   IE: Indices into the vector returned by the events function. The values
%   indicate which event the solver detected.
%
% SAMPLE CALL: 
  % [T,Y,TE,YE,IE, Y1, T1] = SimplePerfusionMode(30,15)
% 

%




% Initial values
X0 = 0.06; %kgm-3 
S0 = 110; %kgm-3
P0=0 ;

%%%%%%%%%%%%
%% Solve ODEs until the terminal event P=Pmx
%%%%%%%%%%%%


tspan = [0 tf]; y0 = [X0 S0 P0];



options = odeset('Events',@(t,y) MyEvent(t,y,F,V,Rp,Pmx));
[T,Y,TE,YE,IE] = ode23(@(t,y) process(t,y,mumax,Ksx,Kix,Pix,Pmx,Kis,Pis,Pms,qsmax,Kss,alp,qpmax,Ksp,Kip,Pip,Pmp,F,S,Rx,Rs,Rp,V),tspan,y0,options);


%%%%%%%%%%%%%
% Solving ODE for dx/dt=0
%%%%%%%%%%%%% IMportanT HERE

tspan1=[TE tf]
  y01=[Y(end,1) Y(end,2) Y(end,3)]

if (TE==tf)
    disp('use lower stream feed to guarantee massive growth');
    return
else
  [T1,Y1] = ode23(@(t,y) process2(t,y,mumax,Ksx,Kix,Pix,Pmx,Kis,Pis,Pms,qsmax,Kss,alp,qpmax,Ksp,Kip,Pip,Pmp,F,S,V,Rx,Rp,Rs),tspan1,y01);
end


% --------------------------------------------------------------
% Plotting AREA
not displayed

%% Process Model equations - see 
% --------------------------------------------------------------
function dydt = process(t,y,mumax,Ksx,Kix,Pix,Pmx,Kis,Pis,Pms,qsmax,Kss,alp,qpmax,Ksp,Kip,Pip,Pmp,F,S,Rx,Rp,Rs,V)
    mu = (mumax*y(2)) / (Ksx+y(2));
dydt = [-F/V*(1-Rx)*y(1)+y(1)*mu*(1-(y(3)-Pix)/(Pmx-Pix))*(Kix/(Kix+y(2))); %X
F/V*(S-(1-Rs)*y(2))-(qsmax*(y(2)/(Kss+y(2)))*(1-((y(3)-Pis)/(Pms-Pis)))*(Kis/(Kis+y(2))))*y(1);
     -F/V*(1-Rp)*y(3)+F/V*(1-Rx)*y(1)+y(1)*mu*(1-(y(3)-Pix)/(Pmx-Pix))*(Kix/(Kix+y(2)))*alp+qpmax*(y(2)/(Ksp+y(2)))*(1-(y(3)-Pip)/(Pmp-Pip))*y(1)*(Kip/(Kip+y(2)));
     ];  
end   
% --------------------------------------------------------------
function [value,isterminal,direction] = MyEvent(t,y,F,V,Rp,Pmx)

value = y(3)-Pmx;     % Detect height = 0 STOP the Product inhibitory effect
  isterminal = 1;   % Stop(!!!!) the integration
direction =  1;   % Positive direction only
end
% --------------------------------------------------------------


function dydt = process2(t,y,mumax,Ksx,Kix,Pix,Pmx,Kis,Pis,Pms,qsmax,Kss,alp,qpmax,Ksp,Kip,Pip,Pmp,F,S,V,Rx,Rp,Rs)
    mu = (mumax*y(2)) / (Ksx+y(2));
dydt = [ 0; %X
      F/V*(S-(1-Rs)*y(2))-(qsmax*(y(2)/(Kss+y(2)))*(1-((y(3)-Pis)/(Pms-Pis)))*(Kis/(Kis+y(2))))*y(1);
     -F/V*(1-Rp)*y(3)+qpmax*(y(2)/(Ksp+y(2)))*(1-(y(3)-Pip)/(Pmp-Pip))*y(1)*(Kip/(Kip+y(2)));
    ]; 
    end
% --------------------------------------------------------------

end
4

1 回答 1

0

我不确定您对其中的哪一个感兴趣，但这里有一些可能的解决方案：

如果你的目标是...

  • 在调试时停止程序而不编辑代码：使用断点
  • 通过编辑代码来停止程序：使用keyboard命令
  • 运行代码的不同部分，具体取决于您是自己运行还是在其他地方使用：使用isdeployed命令
  • 根据您是在调试还是运行它来运行不同的代码段：为程序提供一个布尔输入变量，比如说nowDebugging并使用依赖于它的 if 语句。
于 2013-01-13T12:11:03.520 回答