Since the ODE integration funcitons only use the first output of the ODE functiton, one option is to add them as outputs, then when the function has finished, call it with the final results and get the additional outputs.
That requires slightly re-writing the function declarattion:
function [S,D,kk,k_L_a_o_w,UG] = Sec_model_fun_for_optimization_P_ideal(t,x,Q0,T1,Kc_Total_gases,tauI_Total_gases)
Testing that —
Q0=580;
T1=230;
moleWt=[28;56;84;112;140;168;156]; % mole weight C2,C4,...,C12,C11 [g/mol]
Kc_Total_gases=50;
tauI_Total_gases=50;
y0 = [0.258176232100050 0 0 0 0 0 0.105663461538462 0.0159368044506204 0 0 0 0 0 0.234807692307692 0];
%options = odeset('RelTol',1e-5,'AbsTol',1e-7);
[t,y]= ode23s(@(t,y) Sec_model_fun_for_optimization_P_ideal(t,y,Q0,T1,Kc_Total_gases,tauI_Total_gases),[0 18000],y0); %Sec_model_fun_for_optimization_P_ideal is the function name
figure
plot(t, y)
grid
title('Linear x')
figure
semilogx(t, y)
grid
title('Log x')
endval = numel(t);
[S,D,kk,k_L_a_o_w,UG] = Sec_model_fun_for_optimization_P_ideal(t(endval),y(endval,:),Q0,T1,Kc_Total_gases,tauI_Total_gases)
function [S,D,kk,k_L_a_o_w,UG] = Sec_model_fun_for_optimization_P_ideal(t,x,Q0,T1,Kc_Total_gases,tauI_Total_gases)
%%Dynamic state inputs
n_C2_gas=x(1); %Ethylene mols in gas phase[mol]
n_C4_gas=x(2); %Butene mols in gas phase[mol]
n_C6_gas=x(3); %Hexene mols in gas phase[mol]
n_C8_gas=x(4); %Octene mols in gas phase[mol]
n_C10_gas=x(5); %Decene mols in gas phase[mol]
n_C12_gas=x(6); %Dodecene mols in gas phase[mol]
n_C11_gas=x(7); %Undecene mols in gas phase[mol]
n_C2_liquid=x(8); %Ethylene mols in liquid phase[mol]
n_C4_liquid=x(9); %Butene mols in liquid phase[mol]
n_C6_liquid=x(10); %Hexene mols in liquid phase[mol]
n_C8_liquid=x(11); %Octene mols in liquid phase[mol]
n_C10_liquid=x(12); %Decene mols in liquid phase[mol]
n_C12_liquid=x(13); %Dodecene mols in liquid phase[mol]
n_C11_liquid=x(14); %Undecene mols in liquid phase[mol]
e_integral = x(15); %Error[mol/s]
%% Constants
%Q0=350; % Q_G etylene inflow (ml/min)
Q1=Q0*1e-6/60; % Q_G ethylene inflow (m3/s)
Q2=0; % Q_G butene inflow
Q3=0; % Q_G hexene inflow
Q4=0; % Q_G octene inflow
Q5=0; % Q_G decene inflow
Q6=0; % Q_G dodecene inflow
Q7=0; % Q_G undecane inflow
P1=36e5; % ethylene inflow pressure [Pa]
%T1=230+273.15; % T_Ethylene [K]
T2=230+273.15; % T_ref [K]
R=8.314; % gas constant [J/(mol.K)]
C1=P1/(R*T1); % ethylene inlet gas concentration [mol/m^3]
moleWt = [28;56;84;112;140;168;156];
VR=300e-6; % reactor volume [m^3]
VG=250e-6; % gas volume [m^3]
VL=50e-6; % liquid volume [m^3]
K=[3.24;2.23;1.72;0.2;0.1;0.08;0.09]; % solubility [nondim]
wc=(0.3+0.25); % catalyst weight [g]
kref=[2.224e-4;1.533e-4;3.988e-5;1.914e-7;4.328e-5;...
2.506e-7;4.036e-5;1.062e-6;6.055e-7;]; % rate at Tref=230C [mol/(s.g_cat)]
Eact=[109.5; 15.23; 7.88; 44.45; 9.438; 8.426; 10.91; 12.54; 7.127]; % activation energy [J/mol];
k=kref.*exp(-Eact*(1/T1-1/T2)/R); % rate at T=T2 [mol/(s.g)]
%%Calculation of kL, kH, D_i_M
%%Calculation of diffusivity coefficient
M = [28;56;84;112;140;168;156]; %Species molecular weight[g/mol]
P =P1.*9.86923e-6; %Pressure in atm
v_C=15.9; %Diffusion volume of carbon[cm3]
v_H=2.31; %Diffusion volume of carbon[cm3]
D = zeros(7,1);
for i=1:7
if i==1
a=2;
b=4;
elseif i==2
a=4;
b=8;
elseif i==3
a=6;
b=12;
elseif i==4
a=8;
b=16;
elseif i==5
a=10;
b=20;
elseif i==6
a=12;
b=24;
else
a=11;
b=24;
end
D(i) = 1e-3.*T1.^1.75.*(1./M(i)+1./M(1))./(P.*((a.*v_C).^(1/3)+(b.*v_H).^(1/3)).^2)*1e-4; %Diffusivity coefficient of species i in monomer(C2H4) [m2/s]
end
D_o_w = 0.3173e-4; %Diffusivity of oxygen in water based on ref[m2/s]
%%kki calulations:
d_T = 0.125; %Tank diameter[m]
d_I = 0.0315; %Impeller diameter[m]
UG = Q1./(pi .* (d_T./2)^2); %Superficial gas velocity[m/s]
N_cd = 4.*Q1.^0.5*d_T.^0.25./d_I.^2; %minimum impeller speed [rev/s]
N_I = 14.166667; %impeller speed [rev/s]
k_L_a_o_w = 3.35 .* (N_I./N_cd).^1.464.*UG; %Volumetrice gas-liquid coefficient of water and oxygen[1/s]
lambda_L =5.515e-4; %Parameter of mass-transfer coefficient[-]
lambda_H =2.320e-6; %Parameter of mass-transfer coefficient[-]
sigma_L =4.095; %Parameter of mass-transfer coefficient[-]
sigma_H =2.207; %Parameter of mass-transfer coefficient[-]
kk = zeros(7,1); %Volumetrice gas-liquid coefficient of species i
for i=1:7
if i == 1 || i==2 || i==3
kk(i) = lambda_L.*k_L_a_o_w.*(D(i)./D_o_w).^sigma_L;
else
kk(i) = lambda_H.*k_L_a_o_w.*(D(i)./D_o_w).^sigma_H;
end
end
% Specify initial conditions
xinit=zeros(15,1); % initial state vector
xinit(1)=C1*VR; % initial ethylene in gas (mol)
xinit(14)=36.63/156; % initial undecane in liquid (mol)
xinit(7) = xinit(14)*VG*K(7)/VL; % initial undecane in gas (mol)
xinit(8) = xinit(1)*VL/(K(1)*VG); % initial ethylene in liquid (mol)
xinit(15)=Q1*C1; % initial outflow rate (mol/s)
%%Setpoint
PG_set=36; %Set point pressure [bar]
%Set Point Tracking & Load Rejection
%t1 = 1800; t2 = 3600; t3 = 5400; t4 = 7200; t5 = 9000; t6 = 10800; t7 = 12600; t8 = 14400; t9 = 16200; t10 = 18000;
% Step 1
% +5% change in set point
%if t >= t1 && t <= t2
% nGin_setpoint = 0.616018502797380*1.05;
% Step 2
% -5% change in set point and disturbance
%elseif t >= t3 && t <= t4
% nGin_setpoint = 0.616018502797380*0.95;
% Step 3
% +10% change in setpoint
%elseif t >=t5 && t<= t6
% nGin_setpoint = 0.616018502797380*1.1;
% Step 4
% -10% change in set point
%elseif t >=t7 && t <= t8
% nGin_setpoint = 0.616018502797380*0.9;
% Step 5
% +20% change in set point
%elseif t >= t9 && t <= t10
% nGin_setpoint = 0.616018502797380*1.2;
%end
P_cal = sum(x(1:7))*R*T1/VR*1e-5;
e_P_gases = P_cal - PG_set;
F_G_R = Kc_Total_gases*(e_P_gases+tauI_Total_gases*e_integral);
%Right-hand side evaluation of the dynamic model (DAE set)
S1 = Q1*C1-F_G_R*x(1)/sum(x(1:7))-VR*kk(1)*(x(1)/VG-K(1)*x(8)/VL); % gas phase ethylene (mol/s)
S2 = Q2-F_G_R*x(2)/sum(x(1:7))-VR*kk(2)*(x(2)/VG-K(2)*x(9)/VL); % gas phase butene (mol/s);
S3 = Q3-F_G_R*x(3)/sum(x(1:7))-VR*kk(3)*(x(3)/VG-K(3)*x(10)/VL); % gas phase hexene (mol/s);
S4= Q4-F_G_R*x(4)/sum(x(1:7))-VR*kk(4)*(x(4)/VG-K(4)*x(11)/VL); % gas phase octene (mol/s);
S5= Q5-F_G_R*x(5)/sum(x(1:7))-VR*kk(5)*(x(5)/VG-K(5)*x(12)/VL); % gas phase decene (mol/s);
S6= Q6-F_G_R*x(6)/sum(x(1:7))-VR*kk(6)*(x(6)/VG-K(6)*x(13)/VL); % gas phase dodecene (mol/s);
S7= Q7-F_G_R*x(7)/sum(x(1:7))-VR*kk(7)*(x(7)/VG-K(7)*x(14)/VL); % gas phase undecane (mol/s) ;
S8= VR*kk(1)*(x(1)/VG-K(1)*x(8)/VL)+wc*(-2*k(1)*x(8)^2/VL^2-k(2)*x(8)*x(9)/VL^2-k(3)*x(8)*x(10)/VL^2-k(5)*x(8)*x(11)/VL^2-k(7)*x(8)*x(12)/VL^2);
S9= VR*kk(2)*(x(2)/VG-K(2)*x(9)/VL)+wc*(k(1)*x(8)^2/VL^2-k(2)*x(8)*x(9)/VL^2-2*k(4)*x(9)^2/VL.^2-k(6)*x(9)*x(10)/VL^2-k(8)*x(9)*x(11)/VL^2);
S10= VR*kk(3)*(x(3)/VG-K(3)*x(10)/VL)+wc*(k(2)*x(8)*x(9)/VL^2-k(3)*x(8)*x(10)/VL^2-k(6)*x(9)*x(10)/VL.^2-2*k(9)*x(10)^2/VL^2);
S11= VR*kk(4)*(x(4)/VG-K(4)*x(11)/VL)+wc*(k(3)*x(8)*x(10)/VL^2+k(4)*x(9)^2/VL^2-k(5)*x(8)*x(11)/VL^2-k(8)*x(9)*x(11)/VL^2);
S12= VR*kk(5)*(x(5)/VG-K(5)*x(12)/VL)+wc*(k(5)*x(8)*x(11)/VL^2+k(6)*x(9)*x(10)/VL^2-k(7)*x(8)*x(12)/VL^2);
S13= VR*kk(6)*(x(6)/VG-K(6)*x(13)/VL)+wc*(k(7)*x(8)*x(12)/VL^2+k(8)*x(9)*x(11)/VL^2+k(9)*x(10)^2/VL^2);
S14= VR*kk(7)*(x(7)/VG-K(7)*x(14)/VL);
S15= P_cal-PG_set; %Error
S = ([S1; S2; S3; S4; S5; S6; S7; S8; S9; S10; S11; S12; S13; S14; S15]);
end
See if that does what you wantt.
.
