Why i getting complex double in fitting ode with experimental data

I need to optimize kinetic parameter in following equations:

dca/dw=rA
dcb/dw=rB
dcc/dw=rC
dcd/dw=rD
dce/dw=rE
dcg/dw=rG
dP/dw=(-alfa / 2) * (var(8) / T0) * (var(7) / (var(7) / P0)) * (Ft / Ft0);
dT/dW=(((U*a*(Ta-var(8))/roB*10)+((-dHr1)*(-rA1)+(-dHr2)*(-rA2)+(-dHr3)*(-rC3)+(-dHr4)*(-rA4))))/(var(1)*cpA+var(2)*cpB+var(3)*cpC+var(4)*cpD+var(5)*cpE+var(6)*cpG);
where:
rA = -(rA1 + rA2 + rA4);
rB1 = 3.5 * rA1;
rB2 = 6.5 * rA2;
rB3 = 3 * rC3;
rB4 = 15 / 6 * rA4;
rB = -(rB1 + rB2 + rB3 + rB4);
rC1 =rA1;
rC = rC1+rC3;
rD1 = 4 * rA1;
rD2 = 5 * rA2;
rD3 = rC3;
rD4 = 14 / 6 * rA4;
rD = rD1 + rD2 + rD3 + rD4;
rE2 =4 * rA2;
rE3 =4 * rC3;
rE = rE2+rE3;
rG =8 / 6 * rA4;

I used this:

https://www.mathworks.com/matlabcentral/answers/uploaded_files/64693/Igor_Moura.m

and this:

https://www.mathworks.com/matlabcentral/answers/314965-parameter-estimation-for-a-system-of-differential-equations?s_tid=prof_contriblnk

to make function:

function C=kinetics(par,W)
var0=[0.001808;0.02258704;0;0.00057911;0.00003447;0;134;431.15];
[W,Cv]=ode45(@DifEq,W,var0);
function dC=DifEq(W,var)
dcdt=zeros(8,1);
T0=431.15;
a=190.48;
Ta=673.15;
du=0.021;
U=1.22*((var(8)-Ta)/du)^(1/4);
P0=134;
pi = 3.14;
dv = 25 / 1000;
delta = 2 / 1000;
Hkat = 3.25;
roB = 0.965517241/ ((pi / 4) * ((dv - 2 * delta) ^ 2 * Hkat));
Ft=var(1)+var(2)+var(3)+var(4)+var(5)+var(6);
Ct0=P0/(8.314*T0);
Ca= Ct0 * (var(1) / Ft) * (T0 / var(8)) * (var(7) / P0);
Cb= Ct0 * (var(2) / Ft) * (T0 / var(8)) * (var(7) / P0);
Cc= Ct0 * (var(3) / Ft) * (T0 / var(8)) * (var(7) / P0);
Cd= Ct0 * (var(4) / Ft) * (T0 / var(8)) * (var(7) / P0);
Ce= Ct0 * (var(5) / Ft) * (T0 / var(8)) * (var(7) / P0);
Cg= Ct0 * (var(6) / Ft) * (T0 / var(8)) * (var(7) / P0);
FA0= 0.001808;
FB0= 0.02258704;
FC0= 0;
FD0= 0.00057911;
FE0= 0.00003447;
FG0 = 0;
Ft0=FA0+FB0+FC0+FD0+FE0+FG0;
G = 0.0007189;
epsilon = 0.45;
roO = 1.160147;
Dp = 0.003461538;
Ac = 0.00035;
viskA = ((261.5278 + 12.55191 * (var(7) / 100) - 0.55997 * (var(8) - 273.15)) / 1000000) * 3600;
yA = 0.016405;
Ma = 58;
viskB = ((17.87683461 + 4.238289936 * (var(7) / 100) + 0.016904644 * (var(8) - 273.15)) / 1000000) * 3600;
yB = 0.773086;
Mb = 28;
viskC = ((28.30049 + 5.541702 * (var(7) / 100) - 0.00358 * (var(8) - 273.15)) / 1000000) * 3600;
yC = 0.000313;
Mc = 44;
viskE = ((730.3278 + 11.11223 * (var(7) / 100) - 1.11653 * (var(8) - 273.15)) / 1000000) * 3600;
yE = 0.005255;
Me = 18;
viskD = ((-5.38734 + 2.683369 * (var(7) / 100) + 0.877204 * (var(8) - 273.15)) / 1000000) * 3600;
yD = 0.204942;
Md = 32;
viskosmj = (((viskA * yA * (Ma ^ 0.5)) + (viskB * yB * (Mb ^ 0.5)) + (viskC * yC * (Mc ^ 0.5)) + (viskD * yD * (Md ^ 0.5)) + (viskE * yE * (Me ^ 0.5))) / ((yA * (Ma ^ 0.5)) + (yB * (Mb ^ 0.5)) + (yC * (Mc ^ 0.5)) + (yD * (Md ^ 0.5)) + (yE * (Me ^ 0.5)))) / 1000;
alfa = ((2 * G * (1 - epsilon)) / (roO * Dp * epsilon ^ 3 * roB * Ac * P0)) * ((150 * (1 - epsilon) * viskosmj) / Dp + 1.75 * G) * (760 / (60 * 10));
cpA1 = 9.487;
cpA2 = 3.313e-1;
cpA3 = -0.0001108;
cpA4 = -0.000000002821;
cpB1 = 28.106;
cpB2 = -0.00000368;
cpB3 = 1.475e-5;
cpB4 = -0.00000001065;
cpC1 = -13.075;
cpC2 = 0.3484;
cpC3 = -0.0002184;
cpC4 = 0.00000004839;
cpD1 = 32.243;
cpD2 = 0.001923;
cpD3 = 0.00001055;
cpD4 = -0.000000003596;
cpE1 = 19.975;
cpE2 = 0.07343;
cpE3 = -0.00005601;
cpE4 = 0.00000001715;
cpG1 = 1.742;
cpG2 = 3.19e-1;
cpG3 = -0.0002352;
cpG4 = 6.975e-8;
cpA = cpA1 + cpA2 * var(8) + cpA3 * var(8) ^ 2 + cpA4 * (var(8) ^ 3);
cpB = cpB1 + cpB2 * var(8) + cpB3 * var(8) ^ 2 + cpB4 * (var(8) ^ 3);
cpC = cpC1 + cpC2 * var(8) + cpC3 * var(8) ^ 2 + cpC4 * (var(8) ^ 3);
cpD = cpD1 + cpD2 * var(8) + cpD3 * var(8) ^ 2 + cpD4 * (var(8) ^ 3);
cpE = cpE1 + cpE2 * var(8) + cpE3 * var(8) ^ 2 + cpE4 * (var(8) ^ 3);
cpG = cpG1 + cpG2 * var(8) + cpG3 * var(8) ^ 2 + cpG4 * (var(8) ^ 3);
rA1=((par(1) * par(2) * Ca * (Cb^par(3)))) / (1 + par(2) * Ca);
rA2=(par(4) * (Cb ^ par(3)));
rC3=(par(5) * Cc * ((Cb ^ par(6)))/ (Ca ^ par(7)));
rA4 =((par(8) * par(9)*Ca* (Cb ^ par(10)))) / (1 + par(9)*Ca);
dHr1 = -1240070 + (8.039 * (var(8) - 298.15) + (0.024805 / 2) * (var(8) ^ 2 - 298.15 ^ 2) - (0.00012 / 3) * (var(8) ^ 3 - 298.15 ^ 3) + (7.4102e-8 / 4) * (var(8) ^ 4 - 298.15 ^ 4));
dHr2 = -2646628.37 + (48.939 * (var(8) - 298.15) - (0.0279411 / 2) * (var(8) ^ 2 - 298.15 ^ 2) - (0.0001551 / 3) * (var(8) ^ 3 - 298.15 ^ 3) + (1.22801e-7 / 4) * (var(8) ^ 4 - 298.15 ^ 4));
dHr3 = -1428408.55 + (40.9 * (var(8) - 298.15) - (0.026385 / 2) * (var(8) ^ 2 - 298.15 ^ 2) - (1.19333 / 3) * (var(8) ^ 3 - 298.15 ^ 3) + (1.2141E-8 / 4) * (var(8) ^ 4 - 298.15 ^ 4));
dHr4 = -5322.22+(1005.905*(var(8)-298.15)+4.595567/2*(var(8)^2-298.15^2)-0.00224/3*(var(8)^3-298.15^3)+4.37E-7/4*(var(8)^4-298.15^4));
rA = -(rA1 + rA2 + rA4);
rB1 = 3.5 * rA1;
rB2 = 6.5 * rA2;
rB3 = 3 * rC3;
rB4 = 15 / 6 * rA4;
rB = -(rB1 + rB2 + rB3 + rB4);
rC1 =rA1;
rC = rC1-rC3;
rD1 = 4 * rA1;
rD2 = 5 * rA2;
rD3 = rC3;
rD4 = 14 / 6 * rA4;
rD = rD1 + rD2 + rD3 + rD4;
rE2 =4 * rA2;
rE3 =4 * rC3;
rE = rE2+rE3;
rG =8 / 6 * rA4;
dcdt(1)=rA;
dcdt(2)=rB;
dcdt(3)=rC;
dcdt(4)=rD;
dcdt(5)=rE;
dcdt(6)=rG;
dcdt(7)=(-alfa / 2) * (var(8) / T0) * (var(7) / (var(7) / P0)) * (Ft / Ft0);
dcdt(8)=(((U*a*(Ta-var(8))/roB*10)+((-dHr1)*(-rA1)+(-dHr2)*(-rA2)+(-dHr3)*(-rC3)+(-dHr4)*(-rA4))))/(var(1)*cpA+var(2)*cpB+var(3)*cpC+var(4)*cpD+var(5)*cpE+var(6)*cpG);
dC=dcdt;
    end
C=Cv;
end

and i made script:

W=[0
0.118832891
0.148541114
0.178249337
0.20795756
0.237665782
0.267374005
0.297082228
0.326790451
0.356498674
0.356498675
0.386206896
0.386206897
0.415915119
0.415915120
0.445623342
0.475331565
0.505039788
0.53474801
0.564456233
0.594164456
0.623872679
0.653580902
0.683289124
0.712997347
0.74270557
0.772413793
0.831830238
0.89124668
0.95066313
0.965517241];
var=[0.001808	0.022587	0	        0.000579	0.00003447	0	        134	    431.15
0.0017568	0.0225492	0.000048	0.00162	    0.0003141	0.000979	131.75	685.15
0.0017056	0.0225114	0.000096	0.002661	0.00059373	0.001958	129.5	685.15
0.0016544	0.0224736	0.000144	0.003702	0.00087336	0.002937	127.25	686.15
0.0016032	0.0224358	0.000192	0.004743	0.00115299	0.003916	125	    683.15
0.001552	0.022398	0.00024	    0.005784	0.00143262	0.004895	122.75	685.15
0.0015008	0.0223602	0.000288	0.006825	0.00171225	0.005874	120.5	685.15
0.0014496	0.0223224	0.000336	0.007866	0.00199188	0.006853	118.25	686.15
0.0013984	0.0222846	0.000384	0.008907	0.00227151	0.007832	116	    685.15
0.0013472	0.0222468	0.000432	0.009948	0.00255114	0.008811	113.75	685.15
0.001296	0.022209	0.00048	    0.010989	0.00283077	0.00979	    111.5	684.15
0.0012448	0.0221712	0.000528	0.01203	    0.0031104	0.010769	109.25	686.15
0.0011936	0.0221334	0.000576	0.013071	0.00339003	0.011748	107	    684.15
0.0011424	0.0220956	0.000624	0.014112	0.00366966	0.012727	104.75	684.15
0.0010912	0.0220578	0.000672	0.015153	0.00394929	0.013706	102.5	685.15
0.00104	    0.02202	    0.00072	    0.016194	0.00422892	0.014685	100.25	686.15
0.0009888	0.0219822	0.000768	0.017235	0.00450855	0.015664	98	    685.15
0.0009376	0.0219444	0.000816	0.018276	0.00478818	0.016643	95.75	684.15
0.0008864	0.0219066	0.000864	0.019317	0.00506781	0.017622	93.5	685.15
0.0008352	0.0218688	0.000912	0.020358	0.00534744	0.018601	91.25	683.15
0.000784	0.021831	0.00096	    0.021399	0.00562707	0.01958	    89	    686.15
0.0007328	0.0217932	0.001008	0.02244	    0.0059067	0.020559	86.75	685.15
0.0006816	0.0217554	0.001056	0.023481	0.00618633	0.021538	84.5	685.15
0.0006304	0.0217176	0.001104	0.024522	0.00646596	0.022517	82.25	685.15
0.0005792	0.0216798	0.001152	0.025563	0.00674559	0.023496	80	    686.15
0.000528	0.021642	0.0012	    0.026604	0.00702522	0.024475	77.75	686.15
0.0004768	0.0216042	0.001248	0.027645	0.00730485	0.025454	75.5	683.15
0.0004256	0.0215664	0.001296	0.028686	0.00758448	0.026433	73.25	686.15
0.0003744	0.0215286	0.001344	0.029727	0.00786411	0.027412	71	    686.15
0.0003232	0.0214908	0.001392	0.030768	0.00814374	0.028391	68.75	683.15
0.000272	0.021453	0.00144	    0.031809	0.00842337	0.02937	    66.5	685.15
]; 
par0=[0.0001;2250;0.2;0.00001;0.00001;0.1;1;0.001;1;1];
[par,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,par0,W,var);
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(par)
    fprintf(1, '\t\tpar(%d) = %8.5f\n', k1, par(k1))
end
Wv = linspace(min(W), max(W));
Cfit = kinetics(par, Wv);
figure(1)
plot(W, var, 'p')
hold on
plot(Wv,Cfit);
hold off
grid
xlabel('Mass, kg')
ylabel('All, K')

I write t instead of w.

I got complex double, and i got only good agreement for 8th equation.

Please, try my MATLAB code and if you know what is problem, i will appreciate help.

Thanks in advance.

[Walter Roberson]

If var(8) becomes less than Ta then after you subtract that would be negative and raising to 1/4 would be complex.

ode45 does not have any bounds for variables except what you enforce with event functions.

[Bosnian Kingdom]

I raised Ta to 693, because var(8) has values form 431 to 685. I eliminated that problem. I don't understand: ode45 does not have any bounds for variables except what you enforce with event functions. What is your suggest?

[Walter Roberson]

Inside DifEq, var(8) is your 8th boundary condition -- the current value of the 8th parameter, what-ever role that might happen to have in your equations.

You give ode45 initial values for the boundary conditions, but after that they evolve according to the differential equations, and ode45 makes no attempt to control them. If the differential equations happen to have var(8) become 1E10 in the 17th time step, and -1E20 in the 29th time step, then ode45 might insert extra time steps to ensure that it meets integration tolerances, but it will make no attempt to say "Hold on, var(8) has to be strictly between 0 and 634!"

The only way to get ode45 to enforce "var(8) has to be strictly between 0 and 634" is to use the ode options to add an "event function" to test that the values are inside expected bounds, and to signal termination of the ode if not. See the ballode example, which uses this to bounce a ball off of the floor (velocity and accerlation say that the position of the ball should go negative, but the event function catches negative position and the controlling routine responds by changing changing the initial conditions.

[Bosnian Kingdom]

Inside DifEq, var(8) is your 8th boundary condition -- the current value of the 8th parameter, what-ever role that might happen to have in your equations. - var(8) is temperature and in the equation for the tubular reactor dT/dW, T is included in equation.

Var(8) is show good agreement with experimnetal data (see Cfit), but var(1) to var(7) are problem, there are not even close to experimental data (they are even negative, although they must be positive).

[Walter Roberson]

Your fitting parameters, par, are all complex valued. Your program misleads you about that because it uses fprintf(), and fprintf() only displays the real portion of a value.

The problem starts with the very first call to kinectics, invoked by lsqcurvefit with the initial parameters. You invoke ode45 passing in boundary conditions var0, which starts var(8) as 431.15 . When you go to calculate U, you calculate ((var(8)-Ta)/du)^(1/4) . du is a postive constant from a few lines above, so the base value is negative of var(8)-Ta is negative. Ta is 673.15 which is less than 431.15 so indeed ((var(8)-Ta)/du) is negative, giving you a complex value when raised to 1/4, leading to complex U, which polutes everything else.

Raising Ta does not help, as that makes (var(8)-Ta) more negative.

In my tests, no matter which Ta value I choose, at some point var(8) become less than it. Even if I made Ta 0, then eventually var(8) would become negative. If I made du negative against the off-chance that you had reversed the slope, then var(8) would increase until there was a problem.

[Bosnian Kingdom]

First of all, thank you very much. I appreciate your help. I tried to fix the value for U (I can do it, for example I wrote U = 400) and you are right there is a problem. Now I do not get a double complex solution, but I get double solutions (see Cfit), but again I have poor agreement with experimental values (compare Cfit and var) for the first seven equations, that is, var (1) to var (7). Is it possible that lsqurvefit can not do this?

I also got this:

Solver stopped prematurely.

lsqcurvefit stopped because it exceeded the function evaluation limit,

options.MaxFunctionEvaluations = 1000 (the default value).

[Walter Roberson]

You can add options to lsqcurvefit to permit to try longer. That might help. But so far 15000 function evaluations hasn't been enough.

[Bosnian Kingdom]

Can you write me that option to try? Thanks.

[Walter Roberson]

Modified version attached. But one million function evaluations was not enough to get the curve fitting to converge.

[Bosnian Kingdom]

Thank you very much. I copied and run, but i got this:

Error using lsqcurvefit (line 251)

Function value and YDATA sizes are not equal.

Error in solve_centi (line 72)

[par,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics, X0, XDATA,

YDATA, LB, UB, OPTIONS);

What is problem?

[Walter Roberson]

Please put a breakpoint in at line 72 and when it is encountered, ask

size(XDATA)
size(YDATA)

I get 31 x 1 and 31 x 8 and it is running for me.

Which MATLAB release are you using?

[Bosnian Kingdom]

Now is ok. I use MATLAB 2016b.

I got negative parameters and samotehnig +0 I:

par(1) = -2.5801 +0 I

par(2) = 2517.5 +0 I

par(3) = 0.16724 +0 I

par(4) = 0.99308 +0 I

par(5) = 0.21907 +0 I

par(6) = -0.056637 +0 I

par(7) = 0.87039 +0 I

par(8) = -0.061563 +0 I

par(9) = 1.4165 +0 I

par(10) = 0.80325 +0 I

Is this to complicate for MATLAB or?

[Walter Roberson]

If the curve fitting ran out of iterations then it might return the initial parameters or it might return the best it had found.

The + 0I is good and is telling you that the imaginary components of the parameters are exactly 0. I put in the code to display the imaginary components of the parameters because earlier we had not realized that we were getting complex valued parameters. It is too much bother to create code that displays the imaginary component only if it is nonzero, easier to always display it.

[Bosnian Kingdom]

Is it possible that there are no values for the par(1) to par(10) for which would be good fitting of the experimental data? All par (par(1) to par(10)) must be positiv.

[Walter Roberson]

Earlier I created the symbolic form of the equations with those var0 boundary conditions. I had hoped that it might be practical to dsolve to find theoretical solutions that could then be optimized. Unfortunately most of the variables involved a ratio of sums of some of the terms, raised to 1/5 or 3/5. dsolve for maple was taking too long so I stopped that. I do not know for sure that closed form solutions are not available but it seems likely. Which leaves simulation.

I have a suspicion that this all can be rephrased as a boundary value problem, but I do not know if doing that would help.

I could perhaps rip it apart and feed it to my pet optimization program, but with 10 parameters the search would not be fast.

You mentioned that the parameters must be positive. Are there any other known bounds? Even order if magnitude would help: sending some of the parameters large enough would drive a number of terms to effectively zero, resulting in a fitting around the remaining terms, which is not necessarily invalid but you can waste a lot of time searching towards infinity looking for a combination that gets you a one bit improvement in results.

[Bosnian Kingdom]

Thank you for trying.

Par1 should be around 2.191e-4

Par 2 around 2616

Par3 around 0.2298

Par4 around 7.028e-5

Par5 around 4.989e-5

Par6 around 0.6345

Par7 around 1.151

Par8, Par9 and Par10 i have no data because I used one kinetic model (with Par1-Par7) and i aded reaction kinetic with Par8, Par9 and Par10.

I hope so this will hep.

Answer 1

[Stephane Dauvillier]

Sorry, I may have been confusing. The LB and UB will be the lower and upper bound of your parameter. Youtr parameter is a 10 elements vector so LB and UB should have 1à elements (that's why you should have a warning message saying

Warning: Length of lower bounds is < length(x); filling in missing lower bounds with -Inf.

I was confused by the 2 variables named both var.

Knowing that, can you set boundaries or do you need to specify constraint (see the docuement page here) ?

For your problem, the column is the "Least Squares" but depending on your constraint you lay have to change the solver.

[Bosnian Kingdom]

Thank you for your efforts. I did as you said:

LB (1) = [0]; UB (1) = [0.0001];

LB (2) = [0]; UB (2) = [2600]

LB (3) = [0]; UB (3) = [0.25];

LB (4) = [0]; UB (4) = [0.00008];

LB (5) = [0]; UB (5) = [0.00005];

LB (6) = [0]; UB (6) = [0.7];

LB (7) = [0]; UB (7) = [2];

LB (8) = [0]; UB (8) = [1];

LB (9) = [0]; UB (9) = [1];

LB (10) = [0]; UB (10) = [1];

And I got this:

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to

its initial value is less than the selected value of the function tolerance.

<stopping criteria details>

Rate Constants:

pair (1) = 1.2517e-12 +0

pair (2) = 127 +0

pair (3) = 0.096206 +0 I

pair (4) = 8.4638e-11 +0

pair (5) = 5e-05 +0 I

pair (6) = 0.028507 +0 I

pair (7) = 0.46484 +0

pair (8) = 0.99999 +0

pair (9) = 1 +0 I

pair (10) = 1.6442e-05 +0

So the parameters pair (1), pair (3), and pair (10) come out of the upper boundaries.

However, there is better agreement with the experimental values ​​for var (1), var (2) var (8).

Sorry, I do not understand this part, how can I do it?

For your problem, the column is the "Least Squares" but depending on your constraint you lay to change the solver.

It is important for me to get par1-par10 for which there is good agreement with experimental values ​​(var and Cfit), it's not a problem that they come out of the upper border.

Answer 2

[Stephane Dauvillier]

Hi it's seems that the variable U is the first to become complex (just put a breakpoint on the following line):

U=1.22*((var(8)-Ta)/du)^(1/4);

In fact on first iteration (var(8)-Ta)/du is negative, which explains why you have a complex number for the variable U

You may want to specify boundaries on var to avoid this

[Bosnian Kingdom]

Hi, you are ight, but i can't, because var(8) have initilal values of 431.15, so I fixed values for U, and then there is no complex results, but I still do not get good agreement with experimental values (see var and Cfit).

Answer 3

[Stephane Dauvillier]

Hi, You can specify lower and upper bounds for your var.

By the way is it normal that your function DifEq compute the following variable but don't use it :

        Cd= Ct0 * (var(4) / Ft) * (T0 / var(8)) * (var(7) / P0);
        Ce= Ct0 * (var(5) / Ft) * (T0 / var(8)) * (var(7) / P0);
        Cg= Ct0 * (var(6) / Ft) * (T0 / var(8)) * (var(7) / P0);

[Bosnian Kingdom]

Hi, in link bellow (created by Mrs. Roberson), there are lower and upper bound are set:

https://www.mathworks.com/matlabcentral/answers/uploaded_files/218595/ode_to_W.m

Is this ok:

X0 = par0;

XDATA = W;

YDATA = var;

LB = [];

UB = [];

OPTIONS = optimset('lsqcurvefit');

OPTIONS = optimset(OPTIONS,'MaxFunEvals', 1e6, 'MaxIter', 5e5); %cannot combine this with other optimset call [par,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics, X0, XDATA, YDATA, LB, UB, OPTIONS);

You're right, I do not use it anywhere, I put it in to see the values for the output concentrations (Cd, Ce and Cg). I can throw them out.

Answer 4

[Stephane Dauvillier]

Yes, except you have to modify LB and UB to specify these boundaries. These should be vectors (for instance LB(1) and UB(1) will be the lower and upper boundaries for var(1) and so on).

If one of the element has no lower or upper boundaries, just set it to -inf or +inf

[Bosnian Kingdom]

I did it:

LB(1)=[0.000272]; UB(1)=[0.001808];

LB(2)=[0.022587];UB(2)=[0.022587]

LB(3)=[0];UB(3)=[0.00144];

LB(4)=[0.000579];UB(4)=[0.031809];

LB(5)=[0.00003447];UB(5)=[0.00842337];

LB(6)=[0];UB(6)=[0.02937];

LB(7)=[66.5];UB(5)=[134];

LB(8)=[431.15];UB(8)=[685.15];

I put lower bond from the:

var=[0.001808 0.022587 0 0.000579 0.00003447 0 134 431.15

0.0017568 0.0225492 0.000048 0.00162 0.0003141 0.000979 131.75 685.15

0.0017056 0.0225114 0.000096 0.002661 0.00059373 0.001958 129.5 685.15

0.0016544 0.0224736 0.000144 0.003702 0.00087336 0.002937 127.25 686.15

0.0016032 0.0224358 0.000192 0.004743 0.00115299 0.003916 125 683.15

0.001552 0.022398 0.00024 0.005784 0.00143262 0.004895 122.75 685.15

0.0015008 0.0223602 0.000288 0.006825 0.00171225 0.005874 120.5 685.15

0.0014496 0.0223224 0.000336 0.007866 0.00199188 0.006853 118.25 686.15

0.0013984 0.0222846 0.000384 0.008907 0.00227151 0.007832 116 685.15

0.0013472 0.0222468 0.000432 0.009948 0.00255114 0.008811 113.75 685.15

0.001296 0.022209 0.00048 0.010989 0.00283077 0.00979 111.5 684.15

0.0012448 0.0221712 0.000528 0.01203 0.0031104 0.010769 109.25 686.15

0.0011936 0.0221334 0.000576 0.013071 0.00339003 0.011748 107 684.15

0.0011424 0.0220956 0.000624 0.014112 0.00366966 0.012727 104.75 684.15

0.0010912 0.0220578 0.000672 0.015153 0.00394929 0.013706 102.5 685.15

0.00104 0.02202 0.00072 0.016194 0.00422892 0.014685 100.25 686.15

0.0009888 0.0219822 0.000768 0.017235 0.00450855 0.015664 98 685.15

0.0009376 0.0219444 0.000816 0.018276 0.00478818 0.016643 95.75 684.15

0.0008864 0.0219066 0.000864 0.019317 0.00506781 0.017622 93.5 685.15

0.0008352 0.0218688 0.000912 0.020358 0.00534744 0.018601 91.25 683.15

0.000784 0.021831 0.00096 0.021399 0.00562707 0.01958 89 686.15

0.0007328 0.0217932 0.001008 0.02244 0.0059067 0.020559 86.75 685.15

0.0006816 0.0217554 0.001056 0.023481 0.00618633 0.021538 84.5 685.15

0.0006304 0.0217176 0.001104 0.024522 0.00646596 0.022517 82.25 685.15

0.0005792 0.0216798 0.001152 0.025563 0.00674559 0.023496 80 686.15

0.000528 0.021642 0.0012 0.026604 0.00702522 0.024475 77.75 686.15

0.0004768 0.0216042 0.001248 0.027645 0.00730485 0.025454 75.5 683.15

0.0004256 0.0215664 0.001296 0.028686 0.00758448 0.026433 73.25 686.15

0.0003744 0.0215286 0.001344 0.029727 0.00786411 0.027412 71 686.15

0.0003232 0.0214908 0.001392 0.030768 0.00814374 0.028391 68.75 683.15

0.000272 0.021453 0.00144 0.031809 0.00842337 0.02937 66.5 685.15

];

But i get this:

par(1) = 0.0001 +0 I

par(2) = 2250 +0 I

par(3) = 0.2 +0 I

par(4) = 1e-05 +0 I

par(5) = 1e-05 +0 I

par(6) = 0.1 +0 I

par(7) = 1 +0 I

par(8) = 0.001 +0 I

par(9) = 1 +0 I

par(10) = 1 +0 I

As a result for the par gives me initial values:

par0=[0.0001;2250;0.2;0.00001;0.00001;0.1;1;0.001;1;1];

And the is no good agreement with experimnetal data except for var(8).

Answer 5

[Stephane Dauvillier]

Have follow the hyperlink:

Knowing that, can you set boundaries or do you need to specify constraint (see the docuement page here) ?

Your objective function is least square since your want to minimize the square root norm.

[Bosnian Kingdom]

since I use the reaction kinetics in which the known parameters are:

Par1 should be around 2.191e-4

Couple 2 around 2616

Par3 around 0.2298

Par4 around 7,028e-5

Par5 around 4.989e-5

Par6 around 0.6345

Par7 around 1.151

Par8, Par9 and Par10 i have no data because I used one kinetic model (with Par1-Par7) and i aded reaction kinetic with Par8, Par9 and Par10.

I do not have to get the approximate values listed above, but it might just be easier to fix it if I set limits on the above values.