Error when function returns a vector of length 51, but the length of initial conditions vector is 3.

The error results look like this:

    Error using odearguments (line 95)
COBASIMULASIKP_4/FUNGSIDIFF returns a vector of length 51, but the length of initial conditions vector is 3. The vector returned by COBASIMULASIKP_4/FUNGSIDIFF
and the initial conditions vector must have the same number of elements.
-
Error in ode15s (line 150)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
-
Error in cobasimulasiKP_4/fungsi (line 71)
        [t S]=ode15s(@fungsidiff,tdata,[FEB0 FST0 T0]',[],A1,Ea1,A2,Ea2,A3,Ea3);
-
Error in fminsearch (line 200)
fv(:,1) = funfcn(x,varargin{:});
-
Error in cobasimulasiKP_4 (line 51)
[kons fval]=fminsearch(@fungsi,[A10 Ea10 A20 Ea20 A30 Ea30]',options)

Using fminseacrh, The script of the Matlab looks like this:

=====

%Searching for constant A1, Ea1, A2, Ea2, A3, Ea3

%Initial Guess

A10=100;

Ea10=100;

A20=100;

Ea20=100;

A30=100;

Ea30=100;

---

FEB0=FEBdata(1); FST0=FSTdata(1); T0=Tdata(1);

---

options=optimset('display','iter')

[kons fval]=fminsearch(@fungsi,[A10 Ea10 A20 Ea20 A30 Ea30]',options)

fprintf('Value of A1 is %2.4f\n' ,kons(1))

fprintf('Value of Ea1 is %2.4f\n' ,kons(2))

fprintf('Value of A2 is %2.4f\n' ,kons(3))

fprintf('Value of Ea2 is %2.4f\n' ,kons(4))

fprintf('Value of A3 is %2.4f\n' ,kons(5))

fprintf('Value of Ea3 is %2.4f\n' ,kons(6))

fprintf('Value of SSE is %2.4f\n',fval)

    function minimasi=fungsi(kons)
        A1=kons(1);
        Ea1=kons(2);
        A2=kons(3);
        Ea2=kons(4);
        A3=kons(5);
        Ea3=kons(6);

        [t S]=ode15s(@fungsidiff,tdata,[FEB0 FST0 T0]',[],A1,Ea1,A2,Ea2,A3,Ea3);

        SSE1=(FEBdata-S(:,1)).^2.*mean(FEBdata);
        SSE2=(FSTdata-S(:,2)).^2.*mean(FSTdata);
        SSE3=(Tdata-S(:,3)).^2.*mean(Tdata);
        minimasi=sum(SSE1+SSE2+SSE3);
    end

    function dSdr=fungsidiff(t,S,A1,Ea1,A2,Ea2,A3,Ea3)
        FEB=S(1);
        FST=S(2);
        T=S(3);
        PEB=FEB./FTdata;
        PST=FST./FTdata;
        KEQ=PST./PEB;
        k1 = A1.*exp(-Ea1/8.314./T);
        k2 = A2.*exp(-Ea2/8.314./T);
        k3 = A3.*exp(-Ea3/8.314./T);
        rA1 = -k1.*((PEB)-(PST.*KEQ));
        rA2 = -k2.*PEB;
        rA3 = -k3.*PEB;
        dFEBdr = -((rA1)+(rA2)+(rA3)).*(Ab*z*rhoKat);
        dFSTdr = (rA1).*(Ab*z*rhoKat);
        dTdr = (((rA1).*(delHr1))+((rA2).*(delHr2))+((rA3).*(delHr3))).*(Ab*z*rhoKat)./((FEB.*CpEB)+(FST.*CpST));
(line 95) dSdr = [dFEBdr dFSTdr dTdr];
        dSdr = dSdr(:);
    end

For the data: (there are 17 data each variables)

%data for initial flow rate input reactor 1

tdata=[3; 35; 68; 98; 126; 154; 189; 216; 245; 279; 308; 343; 371; 405; 434; 462; 491]'; %DOS

FTdata = [29.2999; 45.7000; 46.5011; 46.2644; 46.8000; 46.4992; 46.9000; 46.8998; 46.4901; 46.2902; 46.2912; 46.2920; 46.1917; 46.1001; 46.0879; 46.0998; 46.1001]'; %total flow rate, ton/h

FEBdata = [29.1080; 45.0561; 45.7191; 45.4445; 46.0436; 45.8337; 46.1378; 46.0039; 45.9417; 45.7248; 45.6731; 45.4102; 45.5947; 45.4548; 45.4932; 45.4213; 45.4751]'; %EB (reactant) component flow rate input, ton/h

FSTdata = [0.0438; 0.2624; 0.3380; 0.3055; 0.3031; 0.3073; 0.4014; 0.5309; 0.2379; 0.2613; 0.2662; 0.5721; 0.3088; 0.3231; 0.2852; 0.3549; 0.2701]'; %ST (product) component flow rate input, ton/h

Tdata = [873.2854; 893.1460; 887.5787; 892.0498; 895.0469; 895.2626; 896.1332; 896.9493; 899.1358; 900.0131; 900.6928; 901.8408; 902.0395; 902.2294; 902.2361; 902.2125; 904.1859]'; %reactor 1 temperature inlet, Kelvin

---

D = 4.7124; %cone diameter (catalyst bed inner diameter), meter

Ab = (22/7)*D; %circumference inside catalyst bed, meter

z = 11.127; %height of catalyst bed, meter

rhoKat = 1422; %catalyst density (bulk), kg/m3

---

delHr1_298 = 117.5704 ; %heat reaction reference 298K reaction 1, kJ/mol = J/gmol

delHr2_298 = 101.5 ; %heat reaction reference 298K reaction 2, kJ/mol = J/gmol

delHr3_298 = -65.06 ; %heat reaction reference 298K reaction 3, kJ/mol = J/gmol

---

CpEB = (-43.1+(7.072e-01.*Tdata+(-4.811e-04.*Tdata.^(2)+(1.301e-07.*Tdata.^(3)))));

CpST = (-28.25+(6.159e-02.*Tdata+(-4.023e-04.*Tdata.^(2)+(9.935e-08.*Tdata.^(3)))));

---

delHr1 = (delHr1_298)+((-1.6784e+04)+(72.828.*Tdata+(-6.3452e-02.*Tdata.^(2)+(2.8552e-05.*Tdata.^(3)... +(-5.7763e-09.*Tdata.^(4)))))); %J/gmol

delHr2 = (delHr2_298)+(((-8.2097e+04)+(45.474.*Tdata+(-7.2449e-02.*Tdata.^(2)+(4.4272e-05.*Tdata.^(3)... +(-1.0313e-08.*Tdata.^(4))))))); %J/gmol

delHr3 = (delHr3_298)+(((-1.2219e+04)+(79.247.*Tdata+(-1.5784e-01.*Tdata.^(2)+(1.0636e-04.*Tdata.^(3)... +(-2.4789e-08.*Tdata.^(4))))))); %J/gmol

---

My guess for the error is for rA1 value has the value of PEB, PST, and KEQ which they have 17 data each. 17 data times 3 (variables) equals to 51 data. However, I don't know how to change length vector of 51 into 3 (like initial condition).

Answer 1

[Torsten]

Take a look at the example

"ODE with Time-Dependent Terms"

under

https://de.mathworks.com/help/matlab/ref/ode45.html

It will show you how to proceed.

Best wishes

Torsten.