clear all
close all
clc
for i=1:1
Ss(1)=0.000202028;
li=0.0;
Ma=log(0.23178)+8.2+2*0.00132-li;
Mb=log(0.04)+8.2+0.00132-li;
Mc=log(0.5087)+2*15.43+0.00132-li;
Md=log(0.0006)+2*0.12-li;
Me=log(4e-5)+15.43+0.12-li;
Mf=log(7e-12)+012-li;
Mg=log(3.082e-9)+2*0.00132-li;
Mh=log(1.203e-6)+2*15.43-li;
Mi=log(3.004e-8)+15.43+0.00132-li;
Mq=log(3.03e-10)+0.00132-li;
Mk=log(2e-4)+15.43-li;
Ml=log(2e-4)+3.92+2*0.00132-li;
Maa=log(2.028e-6)+3.92+3*0.00132-li;
Mbb=log(0.0372)+60.7+2*0.12-li;
Mcc=log(5.3e-3)+60.7-li;
Mdd=log(7.2e-8)+2*15.43+2*0.00132-li;
Mee=log(7.2e-2)+8.2+3*15.43+0.12-li;
Mff=log(1.113e-12)+6*8.2+6*15.43+0.00132-li;
Mgg=log(5.42e-10)+7*8.2+8*15.43-li;
Mhh=log(9.93e-9)+8.2+4*0.12-li;
Mii=log(0.103976653094887)+4*8.2+8*15.43+2*0.00132-li;
functiun=@(x)[Ma+log(x(1)/(sum(x(1:21))))+x(22)+2*x(27)-1.348134673384821e+01-li;
Mb+log(x(2)/(sum(x(1:21))))+x(22)+x(27)-9.964888350263596e+00-li;
Mc+log(x(3)/(sum(x(1:21))))+2*x(23)+x(27)-6.037084634569317e+01-li;
Md+log(x(4)/(sum(x(1:21))))+2*x(24)+1.435716180549626e+01-li;
Me+log(x(5)/(sum(x(1:21))))+x(23)+x(24)-1.084673779229932e+01-li;
Mf+log(x(6)/(sum(x(1:21))))+x(24)+3.925022193374647e+01-li;
Mg+log(x(7)/(sum(x(1:21))))+2*x(27)+3.919009420009697e+01-li;
Mh+log(x(8)/(sum(x(1:21))))+2*x(23)-3.445861575805654e+01-li;
Mi+log(x(9)/(sum(x(1:21))))+x(23)+x(27)+3.778832020100956e+00-li;
Mq+log(x(10)/(sum(x(1:21))))+x(27)+4.383193662083836e+01-li;
Mk+log(x(11)/(sum(x(1:21))))+x(23)-1.382561361716752e+01-li;
Ml+log(x(12)/(sum(x(1:21))))+x(26)+2*x(27)+9.189106382832476e+00-li;
Maa+log(x(13)/(sum(x(1:21))))+x(26)+3*x(27)+1.836900094447068e+01-li;
Mbb+log(x(14)/(sum(x(1:21))))+x(25)+2*x(24)-1.152971069645939e+02-li;
Mcc+log(x(15)/(sum(x(1:21))))+x(25)-1.109199030831519e+02-li;
Mdd+log(x(16)/(sum(x(1:21))))+2*x(23)+2*x(27)-2.883208056413929e+01-li;
Mee+log(x(17)/(sum(x(1:21))))+x(22)+3*x(23)+x(24)-1.039578216800678e+02-li;
Mff+log(x(18)/(sum(x(1:21))))+6*x(22)+6*x(23)+x(27)-2.285147159127297e+02-li;
Mgg+log(x(19)/(sum(x(1:21))))+7*x(22)+8*x(23)-3.190084897710221e+02-li;
Mhh+log(x(20)/(sum(x(1:21))))+x(22)+4*x(24)+1.949541071777866e+01-li;
Mii+log(x(21)/(sum(x(1:21))))+4*x(22)+8*x(23)+2*x(27)-3.079581022109709e+02-li;
x(1)+x(2)+x(17)+6*x(18)+7*x(19)+x(20)+4*x(21)-.7596866261102260-li;
2*x(3)+x(5)+2*x(8)+x(9)+x(11)+2*x(16)+3*x(17)+6*x(18)+8*x(19)+8*x(21)-2.0654558091417741-li;
2*x(4)+x(5)+x(6)+2*x(14)+x(17)+4*x(20)-0.1476400397270000-li;
x(14)+x(15)-0.0425-li;
x(12)+x(13)-Ss(i)-li;
2*x(1)+x(2)+x(3)+2*x(7)+x(9)+x(10)+2*x(12)+3*x(13)+2*x(16)+x(18)+2*x(21)-1.2206195706978869-li];
format long
x0=11*[0.1,0.003,0.46,0.0001,3e-13*ones(1,6),0.0001,0.0001,2e-8,0.03,0.005,7e-8,7e-2,1e-12*ones(1,3),0.1,2e-3*ones(1,6)];
lb=10*[0.2,0.02,0.48,0.0001,3e-13*ones(1,6),0.0001,0.0001,2e-8,0.03,0.005,7e-8,7e-2,1e-12*ones(1,3),0.1,2e-3*ones(1,6)];
ub=12*[0.2,0.04,0.43,0.08*ones(1,17),0.092,80*ones(1,6)];
options = optimoptions('fsolve','display','final',...
'SubproblemAlgorithm','cg','TolPCG',1e-60,'PlotFcn',@optimplotfirstorderopt);
options.TolFun=1e-60;
options.TolX=1e-70;
options.FinDiffRelStep=1e-35;
options.StepTolerance = 9.000000e-30;
options.MaxFunctionEvaluations = 40000;
options.MaxIterations=4000;
options.FunctionTolerance = 1.000000e-60;
options.DiffMinChange=1e-20;
options.FiniteDifferenceStepSize=eps^(1/100);
[x]=fsolve(functiun,x0,options)
X(1,1:21)=[x(1:21)./(sum(x(1:21)))];
gama=x(22:27);
Ss(i+1)=Ss(i)+0.01;
Yi(i,:)=[Ss(i),X];
end
Yi
