1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156
| clc
clear
funcprot(0)
y_exp(:,1) = [8000 23000 61000 101000 116000 192000 183000 196000]';
y_exp(:,2) = [13000 9000 53000 77000 116000 205000 232000 261000]';
y_exp(:,3) = [6000 26000 46000 67000 136000 192000 236000 303000]';
y_exp(:,4) = [8000 19000 8000 11000 8000 9000 17000 20000]';
y_exp(:,5) = [7000 16000 8000 20000 33000 40000 32000 85000]';
y_exp(:,6) = [9000 16000 15000 16000 36000 74000 57000 79000]';
// data from 12-03-13
y_exp(:,7) = [11000 9000 51000 104000 128000 182000 192000 284000]';
y_exp(:,8) = [9000 7000 14000 26000 46000 100000 114000 162000]';
y_exp(:,9) = [11000 13000 74000 185000 226000 356000 228000 428000]';
y_exp(:,10) = [11000 24000 22000 30000 78000 112000 168000 156000]';
y_exp(:,11) = [13000 27000 111000 201000 222000 272000 316000 332000]';
y_exp(:,12) = [14000 25000 150000 260000 254000 364000 396000 452000]';
// data from 27-03-13
y_exp(:,13) = [10000 10000 16000 33000 46000 84000 75000 104000]';
y_exp(:,14) = [7000 8000 40000 66000 135000 228000 180000 204000]';
y_exp(:,15) = [8000 12000 38000 78000 129000 180000 162000 248000]';
y_exp(:,16) = [6000 10000 11000 14000 11000 28000 23000 27000]';
y_exp(:,17) = [6000 7000 9000 25000 23000 31000 50000 49000]';
y_exp(:,18) = [5000 7000 22000 32000 42000 50000 55000 65000]';
// data from 11-04-13
y_exp(:,19) = [7000 24000 116000 188000 224000 332000 328000 392857]';
y_exp(:,20) = [8000 15000 39000 118000 140000 196000 288000 427429]';
y_exp(:,21) = [8000 12000 43000 103000 144000 192000 300000 480857]';
y_exp(:,22) = [7000 9000 30000 39000 82000 128000 160000 232571]';
y_exp(:,23) = [5000 9000 25000 51000 75000 136000 196000 238857]';
y_exp(:,24) = [8000 11000 27000 45000 99000 148000 184000 257714]';
Fluxi = [148.27 148.71 140.59 51.65 54.36 52.76 148.27 148.71 140.59 233.60 238.9 244.0 148.27 148.71 140.59 51.65 47.12 43.29 241.30 241.80 242.50 ...
141.93 160.74 159.9]';
DOi = [0.425 0.416 0.414 1.286 0.374 0.191 1.325 0.399 0.149 1.358 0.363 0.221 0.384 0.893 0.379 0.194 0.917 0.382 0.855 0.25 0.892 0.885 0.87 0.887]';
//Jb
function [v] = Jb(t, y, k, j, DOi, Fluxi)
g=12.789;
h=0.000196202;
i=54.7014;
v = g + h*y+ i*DOi(j);
endfunction
//Flux
function [z] = Flux(t, y, k, j, DOi, Fluxi)
a=0.902623;
b=0.02971;
c=0.0005421;
d=0.00000557133;
e=0.0000000299771;
f=0.0000000000653467;
v=Jb(t, y, k, j, DOi, Fluxi);
z=Fluxi(j)*(a-b.*v+c.*v.^2-d.*v.^3+e.*v.^4-f.*v.^5);
endfunction
// Model jour
function [ydot] = mymodel1(t, y, k, j, DOi, Fluxi)
z = Flux(t, y, k, j, DOi, Fluxi);
ydot = k(1).*y.*(z./(z+k(2)));
endfunction
//Model nuit
function [ydot] = mymodel2(t, y, k, j, DOi, Fluxi)
ydot=0;
endfunction
y_calc=[];
//Least squares
function [val] = L_Squares(k)
val=[];
u = k(1);
Ks = k(2);
for j=1:24;
y_temp=[]
t = 0;
while t<14;
temp = (t-int(t));
if (temp<0.25) then
if int(t)==0 then
y0 = y_exp(1,j);
//printf('y0-i -- %g\n',y0)
ta=[0:pas:0.25];
y=ode(y0,0,ta,mymodel1)
else
y0 = y_trans3(1,$);
//printf('y0-1 -- %g\n',y0)
ta=[pas:pas:0.25];
y=ode(y0,0,ta,mymodel1)
end
y_trans1=y;
t=t+0.25;
elseif (temp>=0.25)&(temp<0.75) then
tb=[pas:pas:0.5];
y0 = y_trans1(1,$);
//printf('y0-1 -- %g\n',y0)
y=ode(y0,0,tb,mymodel2)
y_trans2=y;
t=t+0.5
elseif (temp>=0.75) then
tc=[pas:pas:0.25];
y0 = y_trans2(1,$);
//printf('y0-1 -- %g\n',y0)
y=ode(y0,0,tc,mymodel1)
y_trans3=y;
y_trans=[y_trans1'; y_trans2'; y_trans3'];
t=t+0.25
y_temp=[y_temp;y_trans];
end
end
y_calc(:,j)=y_temp;
for m=1:8;
n=(m*2/pas+1-2/pas);
y_transb(m,j)=y_calc(n,j);
end
resid = (y_transb(:,j) - y_exp(:,j)).*(y_transb(:,j) - y_exp(:,j)); //matrice of residues
resid = sum(sum(resid));//residue for the experiment j
Residus(1,j)=resid(1,1);
end
val = sum(sum(Residus));//global residue
// printf('val -- %g\n',val)
// val = 1e12
endfunction
// starting guess
u = 0.25;
Ks = 50;
k = [u;Ks];
pas=0.01
//---------------------------------------------------------------------------
tic
options = optimset("TolFun",1D-04,"TolX",1.D-04,"MaxFunEvals",1.0D+06,"MaxIter",1.0D+06);
[param,L_Squares,flag,details] = fminsearch(L_Squares,k,options);
elapsed_time = toc();
parametres = param';
u = parametres(1); Ks = parametres(2);
//clc
parametres
if flag == 1 then
printf('Le système converge\n\n');
else
printf('Le système ne converge pas\n\n');
end;
printf('nombre d itération ---------------- %g\n',details.iterations);
printf('Somme des carrés ----------------- %g \n',L_Squares);
printf('Nombre d appels du modèle ------- %g\n',details.funcCount);
printf('Algorithme de minimisation------------ %s\n',details.algorithm);
printf('µ ');
write(%io(2),u);
printf('Ks :');
write(%io(2),Ks);
if elapsed_time<60 then
printf('Temps écoulé -------------------- %g s \n\n',elapsed_time);
else
printf('Temps écoulé -------------------- %g min \n\n',elapsed_time);
end; |
Partager