Hi Allen, below is the script I’m using, I need to optimise Beta, Kco, BetaAC and ACKco, at this stage I’m fitting them manually.

Any suggestions?

clear all

#--------Slurry properties------------------------#

global c = 8.7;

global Ore = 1/1000;

vol = Ore*3;

#---------------------Integration function----------------------------------#

function Ddt = myode(t,C,tyd,temp)

AcOn = 1;

#---------------------Experimental parameters-------------------------------#

global Ore = 1/1000;

global dc = 75/1000000;

global ro = 0.450;

vol = Ore*3;

#--------------First Thermodynamic & Kinetic parameters----------------------#

global Beta = 0.45; #This parameter needs to be optimised

global n = 1;

global Kco = 80.8*10^-8#This parameter needs to be optimised

global EaKc = 0;

global Hb = -14000;

#---------Thermodynamic & Kinetic paratmeters-----------#

global BetaAC = 50; #This parameter needs to be optimised

global nAC = 1

global ACKco = 5.5*10^-5; #This parameter needs to be optimised

global AcEaKc = 0;

global AcHb = -14000;

ACConl = 30;

MassAC1 = Ore*ACConl/1000000;

dcAC = 2.265/1000;

#-------------Dynamic modelling equations---------------------------------#

T = interp1(tyd,temp,t);

B = @(T) Beta*exp(-Hb/(8.314*(T+273)));

Kc = @(T) Kco*exp(-EaKc/(8.314*(T+273)));

BAC = @(T) BetaAC*exp(-AcHb/(8.314*(T+273)));

KcAC = @(T) ACKco*exp(-AcEaKc/(8.314*(T+273)));

Cs = @(C) ((C(2))/B(T))^(1/n);

CsAC = @(C) (C(3)/BAC(T))^(1/nAC);

dcbdt = -6*Kc(T)*Ore*(C(1)-Cs(C))/(ro*dc*vol) - AcOn*((6*KcAC(T)*MassAC1*(C(1)-CsAC(C))/(ro*dcAC*vol)));

dqdt = -(6*Kc(T)*Ore*(C(1)-Cs(C))/(ro*dc*vol))*vol/(Ore);

dACdt = AcOn*(6*KcAC(T)*MassAC1*(C(1)-CsAC(C))/(ro*dcAC*vol))*vol/(MassAC1);

Ddt = [dcbdt,dqdt,dACdt];

endfunction

#-------------------Experimental data--------------------#

graf = [0 0.18 14

43 0.13 50

120 0.13 70

450 0.15 70

810 0.12 70

840 0.07 50

870 0.04 14

880 0.03 14

1440 0.03 14

];

#---------Used to determine the size of integration------#

plek = size(graf);

plekhou = plek(1,1);

range = linspace(graf(1,1),graf(plekhou,1));

loading = (graf(1,2) - graf(plek(1,1),2))*vol/Ore +c ;

#-------------Carbon Loading trend ------------------#

koolstofgraf(1,1) = c;

for a = 1:plek-1

koolstofgraf(a+1,1) = koolstofgraf(a,1)+(graf(a,2)-graf(a+1,2))*vol/Ore;

endfor

#-------------Temperature trend------------------#

global temp = graf(:,3);

global tyd = graf(:,1);

#-------------------Initialvalue & Integration-------------------#

init = [graf(1,2),c,0];

opts = odeset('RelTol',1e-2, 'AbsTol',1e-4);

[t,ant] = ode45(@(t,C) myode(t, C, tyd , temp), range, init, opts);

Sent from Mail for Windows 10

Ruan,

From: Help-octave <help-octave-bounces+allen.windhorn=[hidden email]> On Behalf Of Ruan Pieterse

> I’m in desperate need of help, I’m working on my final masters project,

> as part of this I have a non linear ODE which I need to integrate, at the

> same rate I need to fit this data to experimental data I’ve generated.

> At this stage I have a function which I’m passing to ODE45, with 4

> parameters which I need to optimise. I want to use lsqcurvefit but I have

> no idea how to incorporate this.

Probably an optimization algorithm will be easier to work into the method?

Start with some reasonable values for your parameters and run the ODE --

then measure the LSQ error with the experimental data. This is the number

that you are trying to minimize. Feed the ODE error to the optimization

algorithm as the fitting function, and let it adjust the parameters for best

fit. Could be gradient method or genetic algorithm, not sure what will be

the best for this problem.

Can you give us your ODE and some of the data? This sounds like a very

interesting problem.

Regards,

Allen