1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
|
function [y_out varargout] = user_defined(x_in, varargin)
n = size(x_in, 1);
N = size(x_in, 2);
t_step = 0:10:10000;
m1 = 11;
m2 = length(t_step);
varargout{1} = m1;
varargout{2} = m2;
p = x_in;
y0(1,:) = 0*ones(1, N); % IS IT CORRECT TO PUT IT IN THIS WAY = 0* ?
y0(2,:) = CHP_initial*ones(1, N); % The start value of y2 is the sixth input parameter. HOW To PUT IT?
y0(3,:) = 0*ones(1, N);
y0(4,:) = ((1-(0.14*CHP_initial)*0.24285)*1000)/18*ones(1, N); % The start value of y4 is A FUNCTION OF the sixth input parameter: ((1-(0.14*CHP_initial*0.24285)*1000)/18. Is it correct like this?
y0(5,:) = 0.5*ones(1, N);
y0(7,:) = 0*ones(1, N);
y0(8,:) = 0*ones(1, N);
y0(9,:) = 0*ones(1, N);
y0(10,:) = Tj0*ones(1, N); % The start value of y10 is the fifth input parameter. HOW To PUT IT?
y0(11,:) = 0.26*ones(1, N);
y0 = reshape(y0, m1*N, 1);
[t y] = ode23s(@semibatch, t_step, y0, [], p);
y_1 = reshape(y, m2, m1, N); % What does it mean? should we keep it?
y_out = permute(y_1, [2,1,3]); % What does it mean? should we keep it?
function dydt = semibatch(tin, yin, pin)
dydt = zeros(2*size(pin,2), 1);
if t > tadd % HOW AND WHERE should I put the If test?
F=0;
end
j = 1;
for i = 1:size(pin,2)
dydt(j)=[(1/(1+((1-((y(j+10)-0.12)/y(j+10)))/(((y(j+10)-0.12)/y(j+10))*9))))*((pin(1,i)/(y(j+10)-0.12))*(24-y(j))-((0.15*exp(-(18041.857)*((1/y(j+9))-(0.0029411))))*(0.0017029)*(sqrt(y(j)/y(j+3)))*(y(j)*y(j+1)-(1/(0.96*exp((671.157)*((1/y(j+9))-(0.003298)))))*y(j+2)*y(j+3)))+((0.0009*exp(-(2429.636)*((1/y(j+9))-(0.0029411))))*y(j+2))+(1-((y(j+10)-0.12)/y(j+10)))*(((0.00576*exp(-(7409.189)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+4))+((0.00437*exp(-(1804.1857)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+5))+((0.004*exp(-(1804.1857)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+6))-((0.00339*exp(-(5063.74789)*((1/y(j+9))-(0.0029411))))*(0.0017029)*y(j+7)*y(j)*(sqrt(y(j)/y(j+3)))))/((y(j+10)-0.12)/y(j+10)));
dydt(j+1)= -((0.15*exp(-(18041.857)*((1/y(j+9))-(0.0029411))))*( 0.0017029)*(sqrt(y(j)/y(j+3)))*(y(j)*y(j+1)-(1/(0.96*exp((671.157)*((1/y(j+9))-(0.003298)))))*y(j+2)*y(j+3)))-(pin(1,i)*y(j+1)/(y(j+10)-0.12));
dydt(j+2)= ((-pin(1,i)*y(j+2)/(y(j+10)-0.12))+((0.15*exp(-(18041.857)*((1/y(j+9))-(0.0029411))))*( 0.0017029)*(sqrt(y(j)/y(j+3)))*(y(j)*y(j+1)-(1/(0.96*exp((671.157)*((1/y(j+9))-(0.003298)))))*y(j+2)*y(j+3)))-((0.0009*exp(-(2429.636)*((1/y(j+9))-(0.0029411))))*y(j+2))-((0.001*exp(-(2405.581)*((1/y(j+9))-(0.0029411))))*y(j+2))-(1-((y(j+10)-0.12)/y(j+10)))*(((0.00576*exp(-(7409.189)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+4))+((0.00437*exp(-(1804.1857)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+5))+((0.004*exp(-(1804.1857)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+6))+((0.0592*exp(-(8419.53331)*((1/y(j+9))-(0.0029411))))*( 0.0017029)*y(j+7)*y(j+2)*(sqrt(y(j)/y(j+3)))))/((y(j+10)-0.12)/y(j+10)));
dydt(j+3)= ((0.15*exp(-(18041.857)*((1/y(j+9))-(0.0029411))))*(0.0017029)*(sqrt(y(j)/y(j+3)))*(y(j)*y(j+1)-(1/(0.96*exp((671.157)*((1/y(j+9))-(0.003298)))))*y(j+2)*y(j+3)))+((0.001*exp(-(2405.581)*((1/y(j+9))-(0.0029411))))*y(j+2))-(PIN(1,I)*y(j+3)/(y(j+10)-0.12))-(1-((y(j+10)-0.12)/y(j+10)))*((0.000237*exp(-(18041.857)*((1/y(j+9))-(0.0029411))))*( 0.0017029)*y(j+7)*(sqrt(y(j)*y(j+3))))/((y(j+10)-0.12)/y(j+10));
dydt(j+4)= -((0.00576*exp(-(7409.189)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+4));
dydt(j+5)= -((0.00437*exp(-(1804.1857)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+5));
dydt(j+6)= ((0.00437*exp(-(1804.1857)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+5))-((0.004*exp(-(1804.1857)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+6));
dydt(j+7)= (((0.00576*exp(-(7409.189)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+4))+((0.00437*exp(-(1804.1857)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+5))+((0.004*exp(-(1804.1857)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+6)))-((0.000237*exp(-(18041.857)*((1/y(j+9))-(0.0029411))))*(0.0017029)*y(j+7)*(sqrt(y(j)*y(j+3))))-((0.00339*exp(-(5063.74789)*((1/y(j+9))-(0.0029411))))*( 0.0017029)*y(j+7)*y(j)*(sqrt(y(j)/y(j+3))))-((0.0592*exp(-(8419.53331)*((1/y(j+9))-(0.0029411))))*(0.0017029)*y(j+7)*y(j+2)*(sqrt(y(j)/y(j+3))));
dydt(j+8)= ((0.000237*exp(-(18041.857)*((1/y(j+9))-(0.0029411))))*(0.0017029)*y(j+7)*(sqrt(y(j)*y(j+3))))+((0.00339*exp(-(5063.74789)*((1/y(j+9))-(0.0029411))))*(0.0017029)*y(j+7)*y(j)*(sqrt(y(j)/y(j+3))))+((0.0592*exp(-(8419.53331)*((1/y(j+9))-(0.0029411))))*(0.0017029)*y(j+7)*y(j+2)*(sqrt(y(j)/y(j+3))));
dydt(j+9)=(1/(((y(j+10)-0.12)*1.00+0.12*0.93)*2000))*((-(y(j+10)-0.12)*(((0.15*exp(-(18041.857)*((1/y(j+9))-(0.0029411))))*( 0.0017029)*(sqrt(y(j)/y(j+3)))*(y(j)*y(j+1)-(1/(0.96*exp((671.157)*((1/y(j+9))-(0.003298)))))*y(j+2)*y(j+3)))*-5580+((0.001*exp(-(2405.581)*((1/y(j+9))-(0.0029411))))*y(j+2))*-359000+((0.0009*exp(-(2429.636)*((1/y(j+9))-(0.0029411))))*y(j+2))*-163000)-0.12*((((0.00576*exp(-(7409.189)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+4))+((0.00437*exp(-(1804.1857)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+5))+((0.004*exp(-(1804.1857)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+6)))*-230000+(((0.000237*exp(-(18041.857)*((1/y(j+9))-(0.0029411))))*( 0.0017029)*y(j+7)*(sqrt(y(j)*y(j+3))))+((0.00339*exp(-(5063.74789)*((1/y(j+9))-(0.0029411))))*( 0.0017029)*y(j+7)*y(j)*(sqrt(y(j)/y(j+3))))+((0.0592*exp(-(8419.53331)*((1/y(j+9))-(0.0029411))))*( 0.0017029)*y(j+7)*y(j+2)*(sqrt(y(j)/y(j+3)))))*-90000))+pin(3,i)*(pin(5,i)-y(j+9))+24*pin(1,i)*20*(pin(4,i)-y(j+9)));
dydt(j+10)= pin(1,i);
j=j+11;
end |
