How to use one fmincon optimizer in a loop and one optimizer out of that loop?

I am trying to optimize k and Ea values for a set of reactions (ode equations) for two different data sets (t,T). We know that k and Ea values can be different for each data set based on the temperature. Also, there is a stoichiometric factor (alpha) that should be optimized which logically, should be the same for both data sets.

So, I need to optimize the k and Ea values first inside a loop and then take the optimized k and Ea values to be utilized for an outer loop optimization part optimizing the stoichiometric factor (alpha).

My problem is that I do ot know haw to take the optimized k and Ea values and substitute them in ode equations to be utilized for the outer loop optimization part. Also, I do not know what the lenght of "i" should be exactly. I would be thankful if anyone could help me with that.

Here is that part of my code:

function Res = Alpha_objfunc(alpha)
for i = 1:1
    k0_vec0 = [0.63, 0.03, 0.951; 17000, 10000, 100000]; % combination of k and Ea (the second raw)
    A=[];
    b=[];
    Aeq=[];
    beq=[];
    lb=[];
    ub=[];
    options=optimoptions('fmincon','Display','iter','Algorithm','sqp');
    [k0_vec,fval]=fmincon(@Obj_func60,k0_vec0,A,b,Aeq,beq,lb,ub,@Con_func60,options);
    % k0_vec0=k0_vec;
    while fval > 0.000001  % I am not sure if it is useful or not
        
        Ode_Cell_deg60;
        k0_vec = k0_vec(i);
    end
    
    fvalnew = fval(end);
end

Res = fvalnew;  % Taking Rseiduals for alpha optimization part
end

Answer 1

[Torsten]

Why don't you optimize all parameters in one call to fmincon ?

It's not a problem that Ea and k are different for each dataset, but that they share a common value for alpha.

Further I suggest using "lsqcurvefit" instead of "fmincon" which is especially designed for such fitting problems.

Take a look at

https://de.mathworks.com/matlabcentral/answers/43439-monod-kinetics-and-curve-fitting

to see how to proceed.

[Farangis]

Because I think this way it will give different alpha for each dataset.

I have not used "lsqcurvefit" before but if you say that it is designed for such problems I will definitely read its instructions to use it.

Thank you very much for your comment and suggestion.

[Matt J]

Because I think this way it will give different alpha for each dataset.

It can't because the model you give to the solver will have only one alpha parameter in it.

[Torsten]

You cannot get different alpha values if you only define one alpha value to be fitted and use that single value for both data sets.

[Farangis]

Thank you for the comment, so I need to review my code again. Maybe I did not define the data set in a correct way.

[Torsten]

Maybe I did not define the data set in a correct way.

You mean the parameter set here, don't you ?

[Farangis]

No, actually I have an ode function that gives t (time) and C (Consentration) and there are k values that changes versus T (Temperature) and also temperature changes versus t. I have two different data set of (t , T) and I have to solve the ode equations to obtain C for both (t , T) data set in a one code.

Finally, I need to optimize k and Ea values for the two data set (which will be different from each other) and "alpha" (which should be the same for each data set).

So, I think I do not know how to define these two (t , T) data set for ode solver.

Here is my ode:

function dcdt = Cell_deg_60C(t,c,k0_vec,Ea_vec,alpha)

R = 8.314; % J/K.mol

%% For T=60 C

% I need to define T data in a vector like it T = [320.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15]

% and For T=80 C ---> T= [321.15 334.15 348.15 353.15 353.15 353.15 353.15 353.15 353.15 353.15]

%% here I just defined one data set of T, because I have no idea how to define the second one

if t <= 720

T = 297.15 + ((320.15-297.15)*(t/540)); % K , s

else

T = 333.15;

end

%% ------------------------------------------------------------------------------------------------

k01 = k0_vec(1);

k02 = k0_vec(2);

k03 = k0_vec(3);

Ea1 = Ea_vec(1);

Ea2 = Ea_vec(2);

Ea3 = Ea_vec(3);

k1 = (k01.*exp(-Ea1./(R.*T)));

k2 = (k02.*exp(-Ea2./(R.*T)));

k3 = (k03.*exp(-Ea3./(R.*T)));

%% D=1 G=2 M=3 O=4

dcdt = zeros(4,1);

dcdt(1) = - k1.*c(1) - k3.*c(1);

dcdt(2) = k1.*c(1);

dcdt(3) = k1.*c(1) + 2*k3.*c(1) - k2.*c(3);

dcdt(4) = alpha*k2*c(3);

dcdt = [dcdt(1);dcdt(2);dcdt(3);dcdt(4)];

% here is odesolver mfile

c0 = zeros(4,1);

c0(1) = 0.14; % mol/l

c0(2) = 0.00045;

c0(3) = 0.01529;

c0(4) = 0.00101;

k0_vec = [0.63, 0.03, 0.951];

Ea_vec = [17000, 10000, 100000];

alpha = 1.44;

tspan = 60*[

9

12

27

42

57

72

87

102

117

132

]; % second

[t,c] = ode45(@(t,c) Cell_deg_60C(t,c,k0_vec,Ea_vec,alpha),tspan,c0);

plot (t/60,c);

legend('Disac','GISA','Monosac','Otheracids')

xlabel('time (min)');

ylabel('C at T=60C (mol/L)');