2

私は頌歌（特定の製品濃度までの時間の経過に伴う細胞の成長）を解くためのコードを書きました。y(3) の値が定数 (Pmx) より大きい場合にモデルが変更されるため、odeset イベント関数を実装しました。イベントが発生する IF は、関数で指定された値によって異なります。

デバッグのために、if ステートメントを設定して関数を停止し、最初の ode23 関数の結果のみを使用してプロットしたい (まだ実装されていない) (コードを参照)

最終順位の問題もあると思いますよね？thx事前に

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 に答える