problem in function inputs

Hi guys. I have this function that gives me the equations for my DAE system.

function S = Sec_model_fun(t,x,Q0,T1)
%%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]
%FGout = x(16);       %Gas outlet[mol/s]



%% Constants
Kc_Total_gases = 1;
tauI_Total_gases =1;
%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]
F0=0.0179;         % gas outflow rate [mmol/s]
F1=F0.*1e-3;       % gas outflow rate [mol/s]
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]
moleWt=[28;56;84;112;140;168;156];              % mole weight C2,C4,...,C12,C11 [g/mol]
wc=(0.3+0.25)*1e-3;                             % catalyst weight [kg]
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)]
kL=1; kH=1;

% 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)
nToti=sum(xinit(1:7));              % initial moles in gas (mol)

%%Setpoint

nGin_setpoint=0.363839693638512;

%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

e_mol_gases = sum(x(1:7)) - nGin_setpoint;
F_G_R = Kc_Total_gases*(e_mol_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*kL*(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*kL*(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*kL*(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*kH*(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*kH*(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*kH*(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*kH*(x(7)/VG-K(7)*x(14)/VL);    % gas phase undecane (mol/s) ;
S8=  VR*kL*(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*kL*(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*kL*(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*kH*(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*kH*(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*kH*(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*kH*(x(7)/VG-K(7)*x(14)/VL);
S15= sum(x(1:7))-(nGin_setpoint);  %Error
    
S = ([S1; S2; S3; S4; S5; S6; S7; S8; S9; S10; S11; S12; S13; S14; S15]);     


end

I want to solve the S1,S2,..,S15 equations for different Q0 and T1. I have tried this code:

clear variables
clc

%Initial Condition
Q0=350; T1=230;
%S = @Sec_model_fun_for_optimization;
%S1 = S(Q0,T1);
y0 = [0.258176232100050 0 0 0 0 0 0.105663461538462 0.0159368044506204 0 0 0 0 0 0.234807692307692 0 Q0 T1];
%options = odeset('RelTol',1e-5,'AbsTol',1e-7);
[t,y]= ode23s(@Sec_model_fun_for_optimization,[0 18000],y0);

but I have got this error:

Not enough input arguments.

Error in Sec_model_fun_for_optimization (line 27)

Q1=Q0*1e-6/60; % Q_G ethylene inflow (m3/s)

As you can see I have done it in a wrong way and MATLAB can't recognize Q0 and T1 as my function inputs. What should I do to solve this problem??

[naiva saeedia]

Guys I guess I have solved my problem. I used this code:

clear variables
clc

%Initial Condition
Q0=100; T1=230;
%S = @Sec_model_fun_for_optimization;
%S1 = S(Q0,T1);
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(t,y,Q0,T1),[0 18000],y0);

and it worked perfectly. I guess It's correct answer??

[Torsten]

I guess It's correct answer??

Yes, that's the way extra parameters can be passed to your function.

[naiva saeedia]

@Torsten

Thanks a lot for your answer🙏🏻

Answer 1

[naiva saeedia]

I have solved the problem. Here is the code:

clear variables
clc
%Initial Condition
Q0=100; T1=230;
%S = @Sec_model_fun_for_optimization;
%S1 = S(Q0,T1);
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(t,y,Q0,T1),[0 18000],y0);%Initial ConditionQ0=100; T1=230;%S = @Sec_model_fun_for_optimization;%S1 = S(Q0,T1);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(t,y,Q0,T1),[0 18000],y0);

[Torsten]

This is the wrong one ...

[naiva saeedia]

@Torsten

Sorry for that..I fixed it. Thanks for pointing it out🙏🏻