Need a better guess y0 for consistent initial conditions.

12 vues (au cours des 30 derniers jours)
Dursman Mchabe
Dursman Mchabe le 14 Jan 2019
Commenté : Dursman Mchabe le 16 Jan 2019
Hi all,
I get the the error message:
Error using daeic12 (line 166)
Need a better guess y0 for consistent initial conditions.
Error in ode15s (line 310)
[y,yp,f0,dfdy,nFE,nPD,Jfac] = daeic12(odeFcn,odeArgs,t,ICtype,Mt,y,yp0,f0,...
Error in IgorWaterODE15s/kinetics (line 58)
[t,Cv] = ode15s(@DifEq,tspan,c0,options);
Error in lsqcurvefit (line 213)
initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});
Error in IgorWaterODE15s (line 166)
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
Caused by:
Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.
When I run the code below:
function IgorWaterODE15s
function C =kinetics(theta,t)
% Initial conditions
%c0 = [1; 1; 1; 1];
c0 = [0;0;0;0];
%c0 = [3.16E-02;0.33;1.62E-04;1.35E-04];
% Time
t=[0
960
1920
2880
3840
4800
5760
6720
7680
8640
9600
10560
11520
12480
13440
14400
15360
16320
17280
18240
19200
20160
21120
22080
23040
24000
24960
25920
26880
27840
28800
29760
30720];
tspan = linspace(min(t), max(t),33);
tspan = tspan';
% A constant, singular mass matrix
M = [1 0 0 0
0 1 0 0
0 0 0 0
0 0 0 0];
% Options
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',[1e-6 1e-6 1e-6 1e-6],'Vectorized','on');
% Solving
[t,Cv] = ode15s(@DifEq,tspan,c0,options);
%
function dC = DifEq(t,c)
% theta(1) = 5.97E-04;
% theta(2) =1.95E-03;
% theta(3) =8.314;
% theta(4) =323.15;
% theta(5) =5.3E-8;
% theta(6) =143.40;
% theta(7) =6.24;
% theta(8) =5.68E-05;
% theta(9) =2.04E-09;
% theta(10) =2.85E-09;
% theta(11) =1.70E-09;
% theta(12) =5.00E-06;
% theta(13) =4.72E-03;
dC = [(theta(1)./theta(2)).*((1930.65./1000000)- (c(1,:) .*theta(3) .*theta(4) ./ 875000 )) - (c(1,:).* theta(3).* theta(4) - theta(6).* ((c(2,:).* c(4,:).^2)./(c(4,:).^2 + theta(7).* c(4,:) + theta(7).* theta(8))))./((1./theta(12)) + theta(6)./((1 + theta(9).* ((theta(7).* theta(6).* c(1,:).* theta(3).* theta(4) ./ c(3,:)) - ((c(2,:).* theta(7).* c(4,:))./(c(4,:).^2 + theta(7) .* c(4,:) + theta(7) .*theta(8))))./2.85E-09.* ((theta(6).* c(1,:).* theta(3).* theta(4)) - ((c(2,:).* c(4,:).^2)./(c(4,:).^2 + theta(7).* c(4,:) + theta(7).* theta(8)))) + theta(11).* ((theta(8).* theta(7).* theta(6).* c(1,:).* theta(3).* theta(4) ./ (c(3,:)).^2) - ((c(2,:).* theta(7).* theta(8))./(c(4,:).^2 + theta(7) .* c(4,:) + theta(7) .*theta(8))))./ 2.85E-09.* ((theta(6).* c(1,:).* theta(3).* theta(4)) - ((c(2,:).* c(4,:).^2)./(c(4,:).^2 + theta(7).* c(4,:) + theta(7).* theta(8))))).* theta(13)))
(c(1,:).* theta(3).* theta(4) - theta(6).* ((c(2,:).* c(4,:).^2)./(c(4,:).^2 + theta(7).* c(4,:) + theta(7).* theta(8))))./((1./theta(12)) + theta(6)./((1 + theta(9).* ((theta(7).* theta(6).* c(1,:).* theta(3).* theta(4) ./ c(3,:)) - ((c(2,:).* theta(7).* c(4,:))./(c(4,:).^2 + theta(7) .* c(4,:) + theta(7) .*theta(8))))./2.85E-09.* ((theta(6).* c(1,:).* theta(3).* theta(4)) - ((c(2,:).* c(4,:).^2)./(c(4,:).^2 + theta(7).* c(4,:) + theta(7).* theta(8)))) + theta(11).* ((theta(8).* theta(7).* theta(6).* c(1,:).* theta(3).* theta(4) ./ (c(3,:)).^2) - ((c(2,:).* theta(7).* theta(8))./(c(4,:).^2 + theta(7) .* c(4,:) + theta(7) .*theta(8))))./ 2.85E-09.* ((theta(6).* c(1,:).* theta(3).* theta(4)) - ((c(2,:).* c(4,:).^2)./(c(4,:).^2 + theta(7).* c(4,:) + theta(7).* theta(8))))).* theta(13)))
- c(3,:) + (theta(7).* theta(6).* c(1,:).* theta(3).* theta(4) ./ c(3,:)) + 2.* (theta(8).* theta(7).* theta(6).* c(1,:).* theta(3).* theta(4) ./ (c(3,:)).^2)+ (theta(5)./c(3,:))
- c(4,:) + ((c(2,:).* theta(7).* c(4,:))./(c(4,:).^2 + theta(7) .* c(4,:) + theta(7) .*theta(8))) + 2.*((c(2,:).* theta(7).* theta(8))./(c(4,:).^2 + theta(7) .* c(4,:) + theta(7) .*theta(8))) + (theta(5)./c(4,:))];
end
C=Cv
end
t=[0
960
1920
2880
3840
4800
5760
6720
7680
8640
9600
10560
11520
12480
13440
14400
15360
16320
17280
18240
19200
20160
21120
22080
23040
24000
24960
25920
26880
27840
28800
29760
30720];
c= [3.16E-02 0.33 1.62E-04 1.35E-04
3.56E-02 0.64 2.69E-04 2.24E-04
3.83E-02 0.93 4.46E-04 3.72E-04
4.12E-02 1.18 2.51E-03 2.09E-03
4.34E-02 1.40 0.43 0.35
4.58E-02 1.60 1.20 1.00
4.74E-02 1.78 1.70 1.41
4.92E-02 1.94 1.99 1.66
5.08E-02 2.08 2.81 2.34
5.21E-02 2.20 2.95 2.45
5.34E-02 2.31 3.23 2.69
5.49E-02 2.41 3.31 2.75
5.56E-02 2.49 3.38 2.82
5.63E-02 2.57 3.16 2.63
5.72E-02 2.63 3.08 2.57
5.80E-02 2.69 3.38 2.82
5.88E-02 2.74 3.23 2.69
5.92E-02 2.78 3.23 2.69
5.98E-02 2.82 3.62 3.02
6.03E-02 2.85 3.71 3.09
6.06E-02 2.87 3.71 3.09
6.10E-02 2.89 3.54 2.95
6.14E-02 2.91 3.54 2.95
6.15E-02 2.93 3.71 3.09
6.18E-02 2.94 3.79 3.16
6.20E-02 2.95 3.71 3.09
6.22E-02 2.96 3.88 3.24
6.24E-02 2.96 3.88 3.24
6.25E-02 2.97 3.88 3.24
6.26E-02 2.97 4.07 3.39
6.27E-02 2.97 4.07 3.39
6.28E-02 2.98 4.16 3.47
6.29E-02 2.98 4.16 3.47];
%theta0=[1;1;1;1;1;1;1;1;1;1;1;1;1];
%theta0 = [0;0;0;0;0;0;0;0;0;0;0;0;0];
theta0 = [0;1;0;1;0;1;0;1;0;1;0;1;0];
%theta0 = [5.97E-04;1.95E-03;8.314;323.15;5.3E-8;143.40;6.24;5.68E-05;2.04E-09;2.85E-09;1.70E-09;5.00E-06;4.72E-03];
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
%[theta]=lsqcurvefit(@kinetics,theta0,t,c);
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t),33);
tv = tv';
Cv = kinetics(theta, tv)
figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cv);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend('Exp-SO_{2,gas}','Exp-S_{total}','Exp-H^+_{liquid}','Exp-H^+_{Interface}','Mod-SO_{2,gas}','Mod-S_{total}','Mod-H^+_{liquid}','Mod-H^+_{Interface}')
end
What better guess can I try?
  36 commentaires
Torsten
Torsten le 16 Jan 2019
This is the code I tested, and it worked. Since I don't have the optimization toolbox, I can't use "fsolve" or "lsqcurvefit":
function main
t=[0;960;1920;2880;3840;4800;5760;6720;7680;8640;9600;10560;11520;12480;13440;14400;15360;16320;17280;18240;19200;20160;21120;22080;23040;24000;24960;25920;26880;27840;28800;29760;30720];
c= [3.16E-02 0.33 1.62E-04 1.35E-04; ...
3.56E-02 0.64 2.69E-04 2.24E-04; ...
3.83E-02 0.93 4.46E-04 3.72E-04; ...
4.12E-02 1.18 2.51E-03 2.09E-03; ...
4.34E-02 1.40 0.43 0.35; ...
4.58E-02 1.60 1.20 1.00; ...
4.74E-02 1.78 1.70 1.41; ...
4.92E-02 1.94 1.99 1.66; ...
5.08E-02 2.08 2.81 2.34; ...
5.21E-02 2.20 2.95 2.45; ...
5.34E-02 2.31 3.23 2.69; ...
5.49E-02 2.41 3.31 2.75; ...
5.56E-02 2.49 3.38 2.82; ...
5.63E-02 2.57 3.16 2.63; ...
5.72E-02 2.63 3.08 2.57; ...
5.80E-02 2.69 3.38 2.82; ...
5.88E-02 2.74 3.23 2.69; ...
5.92E-02 2.78 3.23 2.69; ...
5.98E-02 2.82 3.62 3.02; ...
6.03E-02 2.85 3.71 3.09; ...
6.06E-02 2.87 3.71 3.09; ...
6.10E-02 2.89 3.54 2.95; ...
6.14E-02 2.91 3.54 2.95; ...
6.15E-02 2.93 3.71 3.09; ...
6.18E-02 2.94 3.79 3.16; ...
6.20E-02 2.95 3.71 3.09; ...
6.22E-02 2.96 3.88 3.24; ...
6.24E-02 2.96 3.88 3.24; ...
6.25E-02 2.97 3.88 3.24; ...
6.26E-02 2.97 4.07 3.39; ...
6.27E-02 2.97 4.07 3.39; ...
6.28E-02 2.98 4.16 3.47; ...
6.29E-02 2.98 4.16 3.47];
theta0 = [5.97E-04;1.95E-03;8.314;323.15;5.3E-8;143.40;6.24;5.68E-05;2.04E-09;2.85E-09;1.70E-09;5.00E-06;4.72E-03];
C = kinetics(theta0,t)
end
function C = kinetics(theta,t)
c0 = [3.16E-02;0.33];
fun = @(x) sum([-x(1)+(theta(7)*theta(6)*c0(1)*theta(3)*theta(4)/x(1))+2.0*(theta(8)*theta(7)*theta(6)*c0(1)*theta(3)*theta(4)/(x(1))^2)+(theta(5)/x(1)); -x(2)+((c0(2)*theta(7)*x(2))/(x(2)^2+theta(7)*x(2)+theta(7)*theta(8)))+2.0*((c0(2)*theta(7)*theta(8))/(x(2)^2+theta(7)*x(2)+theta(7)*theta(8)))+(theta(5)/x(2))].^2);
sol = fminsearch(fun,[1;1]);
fun(sol)
c0 = [c0;sol]
tspan = t;
% A constant, singular mass matrix
M = [1 0 0 0;0 1 0 0;0 0 0 0;0 0 0 0];
% Options
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',[1e-6 1e-6 1e-6 1e-6]);
% Solving
[t,C] = ode15s(@(t,y)DifEq(t,y,theta),tspan,c0,options);
end
function dC = DifEq(t,c,theta)
dC = zeros(4,1);
dC(1) = (theta(1)/theta(2))*((1930.65/1000000)- (c(1)*theta(3)*theta(4) / 875000 )) - (c(1)* theta(3)* theta(4) - theta(6)* ((c(2)* c(4)^2)/(c(4)^2 + theta(7)* c(4) + theta(7)* theta(8))))/((1.0/theta(12)) + theta(6)/((1 + theta(9)* ((theta(7)* theta(6)* c(1)* theta(3)* theta(4) / c(3)) - ((c(2)* theta(7)* c(4))/(c(4)^2 + theta(7)* c(4) + theta(7) *theta(8))))/2.85E-09* ((theta(6)* c(1)* theta(3)* theta(4)) - ((c(2)* c(4)^2)/(c(4)^2 + theta(7)* c(4) + theta(7)* theta(8)))) + theta(11)* ((theta(8)* theta(7)* theta(6)* c(1)* theta(3)* theta(4) / (c(3))^2) - ((c(2)* theta(7)* theta(8))/(c(4)^2 + theta(7) * c(4) + theta(7) *theta(8))))/ 2.85E-09* ((theta(6)* c(1)* theta(3)* theta(4)) - ((c(2)* c(4)^2)/(c(4)^2 + theta(7)* c(4) + theta(7)* theta(8)))))* theta(13)));
dC(2) = (c(1)* theta(3)* theta(4) - theta(6)* ((c(2)* c(4)^2)/(c(4)^2 + theta(7)* c(4) + theta(7)* theta(8))))/((1.0/theta(12)) + theta(6)/((1 + theta(9)* ((theta(7)* theta(6)* c(1)* theta(3)* theta(4) / c(3)) - ((c(2)* theta(7)* c(4))/(c(4).^2 + theta(7) * c(4) + theta(7) *theta(8))))/2.85E-09* ((theta(6)* c(1)* theta(3)* theta(4)) - ((c(2)* c(4)^2)/(c(4)^2 + theta(7)* c(4) + theta(7)* theta(8)))) + theta(11)* ((theta(8)* theta(7)* theta(6)* c(1)* theta(3)* theta(4) / (c(3))^2) - ((c(2)* theta(7)* theta(8))/(c(4)^2 + theta(7) * c(4) + theta(7) *theta(8))))/ 2.85E-09* ((theta(6)* c(1)* theta(3)* theta(4)) - ((c(2)* c(4)^2)/(c(4)^2 + theta(7)* c(4) + theta(7)* theta(8)))))* theta(13)));
dC(3) = - c(3) + (theta(7)* theta(6)* c(1)* theta(3)* theta(4) / c(3)) + 2.0* (theta(8)* theta(7)* theta(6)* c(1)* theta(3)* theta(4) / (c(3))^2)+ (theta(5)/c(3));
dC(4) = - c(4) + ((c(2)* theta(7)* c(4))/(c(4)^2 + theta(7) * c(4) + theta(7) *theta(8))) + 2.0*((c(2)* theta(7)* theta(8))/(c(4)^2 + theta(7) * c(4) + theta(7) *theta(8))) + (theta(5)/c(4));
end
Dursman Mchabe
Dursman Mchabe le 16 Jan 2019
It works!!!! Thank you! Thank You! Thank you!

Connectez-vous pour commenter.

Réponses (0)

Catégories

En savoir plus sur Model Predictive Control Toolbox dans Help Center et File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by