The code I wrote for SQP method optimization is not working.
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
clc;clear all;
syms Rd Rb bd u1 u2 u3 u4 u5 u6 u7 u8 u9 u10 s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 real
R0=1;
gama=0.1;
eps1=1e-5;
eps2=1e-5;
dd=1;
V=1;
% x=[20,40,3,5];
% x0=[0.35,0.42,0.7,0.8];
x0=[0.050,0.020,0.005];
hisx(1,:)=x0;
n=1;
H(:,:,1)=eye(3);
while (dd>=eps2)&&(V>=eps1)&&(n<=50)
if n==1
Rd=x0(1);Rb=x0(2);bd=x0(3);
dR=(Rd-Rb);
Td=4/3*pi*40000*dR^3+pi*3.8*45/0.001*dR^4+2*pi*Rd^2*bd*[15+0.1*(45*Rd/0.001)];
Td0=4/3*pi*15*dR^3+pi*0.1*45/0.001*dR^4+2*pi*Rd^2*bd*[15+0.1*(45*Rd/0.001)];
muJ=Rd^4*bd/[(0.069)^4*0.022-(0.065)^4*(bd+0.002)];
f=-Td;
g1=1-Td/(50*Td0);
g2=1-muJ/0.25;
g3=muJ/0.5-1;
g4=Rb/Rd-1;
g5=1-Rd/0.030;
g6=Rd/0.068-1;
g7=1-Rb/0.010;
g8=Rb/0.020-1;
g9=1-bd/0.003;
g10=bd/0.012-1;
gf=[(2307309447473265*(Rb - Rd)^3)/1073741824 - (17271088261788183*(Rb - Rd)^2)/34359738368 - 9000*pi*Rd^2*bd - 4*pi*Rd*bd*(4500*Rd + 15);
(17271088261788183*(Rb - Rd)^2)/34359738368 - (2307309447473265*(Rb - Rd)^3)/1073741824;
-2*pi*Rd^2*(4500*Rd + 15)]';
gg1=[(((2307309447473265*(Rb - Rd)^4)/4294967296 - (5757029420596061*(Rb - Rd)^3)/34359738368 + 2*pi*Rd^2*bd*(4500*Rd + 15))*(3000*pi*(Rb - Rd)^2 - 900000*pi*(Rb - Rd)^3 + 450000*pi*Rd^2*bd + 200*pi*Rd*bd*(4500*Rd + 15)))/(- 1000*pi*(Rb - Rd)^3 + 225000*pi*(Rb - Rd)^4 + 100*pi*Rd^2*bd*(4500*Rd + 15))^2 - ((17271088261788183*(Rb - Rd)^2)/34359738368 - (2307309447473265*(Rb - Rd)^3)/1073741824 + 9000*pi*Rd^2*bd + 4*pi*Rd*bd*(4500*Rd + 15))/(225000*pi*(Rb - Rd)^4 - 1000*pi*(Rb - Rd)^3 + 100*pi*Rd^2*bd*(4500*Rd + 15));
((17271088261788183*(Rb - Rd)^2)/34359738368 - (2307309447473265*(Rb - Rd)^3)/1073741824)/(225000*pi*(Rb - Rd)^4 - 1000*pi*(Rb - Rd)^3 + 100*pi*Rd^2*bd*(4500*Rd + 15)) - ((3000*pi*(Rb - Rd)^2 - 900000*pi*(Rb - Rd)^3)*((2307309447473265*(Rb - Rd)^4)/4294967296 - (5757029420596061*(Rb - Rd)^3)/34359738368 + 2*pi*Rd^2*bd*(4500*Rd + 15)))/(- 1000*pi*(Rb - Rd)^3 + 225000*pi*(Rb - Rd)^4 + 100*pi*Rd^2*bd*(4500*Rd + 15))^2;
(100*pi*Rd^2*(4500*Rd + 15)*((2307309447473265*(Rb - Rd)^4)/4294967296 - (5757029420596061*(Rb - Rd)^3)/34359738368 + 2*pi*Rd^2*bd*(4500*Rd + 15)))/(- 1000*pi*(Rb - Rd)^3 + 225000*pi*(Rb - Rd)^4 + 100*pi*Rd^2*bd*(4500*Rd + 15))^2 - (2*pi*Rd^2*(4500*Rd + 15))/(225000*pi*(Rb - Rd)^4 - 1000*pi*(Rb - Rd)^3 + 100*pi*Rd^2*bd*(4500*Rd + 15))]';
gg2=[(16*Rd^3*bd)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000);
0;
(4*Rd^4)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000) - (2634287287446093*Rd^4*bd)/(36893488147419103232*((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000)^2)]';
gg3=[-(8*Rd^3*bd)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000);
0;
(2634287287446093*Rd^4*bd)/(73786976294838206464*((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000)^2) - (2*Rd^4)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000)]';
gg4=[-Rb/Rd^2;
1/Rd;
0]';
gg5=[-100/3;
0;
0]';
gg6=[250/17;
0;
0]';
gg7=[0;
-100;
0]';
gg8=[0;
50;
0]';
gg9=[0;
0;
-1000/3]';
gg10=[0;
0;
250/3]';
a=[g1,g2,g3,g4,g5,g6,g7,g8,g9,g10];
V0=max([0 a]);
A=[gg1;gg2;gg3;gg4;gg5;gg6;gg7;gg8;gg9;gg10];
options = optimoptions('quadprog','Algorithm','interior-point-convex','Display','off');
[d,fval,exitflag,output,lambda] = quadprog(eye(3),gf,A,-a',[],[],[],[],[],options);
dd=norm(d);
r0=sum(lambda.ineqlin);
b=[R0 r0];
R=max(b);
beta0=gama*dd^2;
phi0=f+R*V0;
phi=1;
phiii=0;
k=1;
while (phi>phiii)
t0=(1/2)^(k-1);
Rd=x0(1)+t0*d(1);
Rb=x0(2)+t0*d(2);
bd=x0(3)+t0*d(3);
dR=(Rd-Rb);
Td=4/3*pi*40000*dR^3+pi*3.8*45/0.001*dR^4+2*pi*Rd^2*bd*[15+0.1*(45*Rd/0.001)];
Td0=4/3*pi*15*dR^3+pi*0.1*45/0.001*dR^4+2*pi*Rd^2*bd*[15+0.1*(45*Rd/0.001)];
muJ=Rd^4*bd/[(0.069)^4*0.022-(0.065)^4*(bd+0.002)];
f=-Td;
g1=1-Td/(50*Td0);
g2=1-muJ/0.25;
g3=muJ/0.5-1;
g4=Rb/Rd-1;
g5=1-Rd/0.030;
g6=Rd/0.068-1;
g7=1-Rb/0.010;
g8=Rb/0.020-1;
g9=1-bd/0.003;
g10=bd/0.012-1;
a=[g1,g2,g3,g4,g5,g6,g7,g8,g9,g10];
V=max([0 a]);
phi=f+R*V;
phiii=phi0-t0*beta0;
k=k+1;
end
alpha=t0;
x=[Rb Rd bd]
R0=R;
n=n+1
hisx(n,:)=x;
dd
V
f
a
pause(1)
s=alpha*d
z=H*s
Rd=x(1);Rb=x(2);bd=x(3);
u1=lambda.ineqlin(1);u2=lambda.ineqlin(2);u3=lambda.ineqlin(3);u4=lambda.ineqlin(4);
u5=lambda.ineqlin(5);u6=lambda.ineqlin(6);u7=lambda.ineqlin(7);u8=lambda.ineqlin(8);
u9=lambda.ineqlin(9);u10=lambda.ineqlin(10);
yy1=[(250*u6)/17 - (100*u5)/3 - u1*(((17271088261788183*(Rb - Rd)^2)/34359738368 - (2307309447473265*(Rb - Rd)^3)/1073741824 + 9000*pi*Rd^2*bd + 4*pi*Rd*bd*(4500*Rd + 15))/(225000*pi*(Rb - Rd)^4 - 1000*pi*(Rb - Rd)^3 + 100*pi*Rd^2*bd*(4500*Rd + 15)) - (((2307309447473265*(Rb - Rd)^4)/4294967296 - (5757029420596061*(Rb - Rd)^3)/34359738368 + 2*pi*Rd^2*bd*(4500*Rd + 15))*(3000*pi*(Rb - Rd)^2 - 900000*pi*(Rb - Rd)^3 + 450000*pi*Rd^2*bd + 200*pi*Rd*bd*(4500*Rd + 15)))/(- 1000*pi*(Rb - Rd)^3 + 225000*pi*(Rb - Rd)^4 + 100*pi*Rd^2*bd*(4500*Rd + 15))^2) - (17271088261788183*(Rb - Rd)^2)/34359738368 + (2307309447473265*(Rb - Rd)^3)/1073741824 - 9000*pi*Rd^2*bd - (Rb*u4)/Rd^2 + (16*Rd^3*bd*u2)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000) - (8*Rd^3*bd*u3)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000) - 4*pi*Rd*bd*(4500*Rd + 15);
50*u8 - 100*u7 + u4/Rd + u1*(((17271088261788183*(Rb - Rd)^2)/34359738368 - (2307309447473265*(Rb - Rd)^3)/1073741824)/(225000*pi*(Rb - Rd)^4 - 1000*pi*(Rb - Rd)^3 + 100*pi*Rd^2*bd*(4500*Rd + 15)) - ((3000*pi*(Rb - Rd)^2 - 900000*pi*(Rb - Rd)^3)*((2307309447473265*(Rb - Rd)^4)/4294967296 - (5757029420596061*(Rb - Rd)^3)/34359738368 + 2*pi*Rd^2*bd*(4500*Rd + 15)))/(- 1000*pi*(Rb - Rd)^3 + 225000*pi*(Rb - Rd)^4 + 100*pi*Rd^2*bd*(4500*Rd + 15))^2) + (17271088261788183*(Rb - Rd)^2)/34359738368 - (2307309447473265*(Rb - Rd)^3)/1073741824;
(250*u10)/3 - (1000*u9)/3 - u3*((2*Rd^4)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000) - (2634287287446093*Rd^4*bd)/(73786976294838206464*((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000)^2)) + u2*((4*Rd^4)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000) - (2634287287446093*Rd^4*bd)/(36893488147419103232*((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000)^2)) - u1*((2*pi*Rd^2*(4500*Rd + 15))/(225000*pi*(Rb - Rd)^4 - 1000*pi*(Rb - Rd)^3 + 100*pi*Rd^2*bd*(4500*Rd + 15)) - (100*pi*Rd^2*(4500*Rd + 15)*((2307309447473265*(Rb - Rd)^4)/4294967296 - (5757029420596061*(Rb - Rd)^3)/34359738368 + 2*pi*Rd^2*bd*(4500*Rd + 15)))/(- 1000*pi*(Rb - Rd)^3 + 225000*pi*(Rb - Rd)^4 + 100*pi*Rd^2*bd*(4500*Rd + 15))^2) - 2*pi*Rd^2*(4500*Rd + 15)];
Rd=x0(1);Rb=x0(2);bd=x0(3);
yy0=[(250*u6)/17 - (100*u5)/3 - u1*(((17271088261788183*(Rb - Rd)^2)/34359738368 - (2307309447473265*(Rb - Rd)^3)/1073741824 + 9000*pi*Rd^2*bd + 4*pi*Rd*bd*(4500*Rd + 15))/(225000*pi*(Rb - Rd)^4 - 1000*pi*(Rb - Rd)^3 + 100*pi*Rd^2*bd*(4500*Rd + 15)) - (((2307309447473265*(Rb - Rd)^4)/4294967296 - (5757029420596061*(Rb - Rd)^3)/34359738368 + 2*pi*Rd^2*bd*(4500*Rd + 15))*(3000*pi*(Rb - Rd)^2 - 900000*pi*(Rb - Rd)^3 + 450000*pi*Rd^2*bd + 200*pi*Rd*bd*(4500*Rd + 15)))/(- 1000*pi*(Rb - Rd)^3 + 225000*pi*(Rb - Rd)^4 + 100*pi*Rd^2*bd*(4500*Rd + 15))^2) - (17271088261788183*(Rb - Rd)^2)/34359738368 + (2307309447473265*(Rb - Rd)^3)/1073741824 - 9000*pi*Rd^2*bd - (Rb*u4)/Rd^2 + (16*Rd^3*bd*u2)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000) - (8*Rd^3*bd*u3)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000) - 4*pi*Rd*bd*(4500*Rd + 15);
50*u8 - 100*u7 + u4/Rd + u1*(((17271088261788183*(Rb - Rd)^2)/34359738368 - (2307309447473265*(Rb - Rd)^3)/1073741824)/(225000*pi*(Rb - Rd)^4 - 1000*pi*(Rb - Rd)^3 + 100*pi*Rd^2*bd*(4500*Rd + 15)) - ((3000*pi*(Rb - Rd)^2 - 900000*pi*(Rb - Rd)^3)*((2307309447473265*(Rb - Rd)^4)/4294967296 - (5757029420596061*(Rb - Rd)^3)/34359738368 + 2*pi*Rd^2*bd*(4500*Rd + 15)))/(- 1000*pi*(Rb - Rd)^3 + 225000*pi*(Rb - Rd)^4 + 100*pi*Rd^2*bd*(4500*Rd + 15))^2) + (17271088261788183*(Rb - Rd)^2)/34359738368 - (2307309447473265*(Rb - Rd)^3)/1073741824;
(250*u10)/3 - (1000*u9)/3 - u3*((2*Rd^4)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000) - (2634287287446093*Rd^4*bd)/(73786976294838206464*((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000)^2)) + u2*((4*Rd^4)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000) - (2634287287446093*Rd^4*bd)/(36893488147419103232*((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000)^2)) - u1*((2*pi*Rd^2*(4500*Rd + 15))/(225000*pi*(Rb - Rd)^4 - 1000*pi*(Rb - Rd)^3 + 100*pi*Rd^2*bd*(4500*Rd + 15)) - (100*pi*Rd^2*(4500*Rd + 15)*((2307309447473265*(Rb - Rd)^4)/4294967296 - (5757029420596061*(Rb - Rd)^3)/34359738368 + 2*pi*Rd^2*bd*(4500*Rd + 15)))/(- 1000*pi*(Rb - Rd)^3 + 225000*pi*(Rb - Rd)^4 + 100*pi*Rd^2*bd*(4500*Rd + 15))^2) - 2*pi*Rd^2*(4500*Rd + 15)];
y=yy1-yy0
kesi1=dot(s,y)
kesi2=dot(s,z)
if kesi1>=0.2*kesi2
teta=1
else
teta=0.8*kesi2/(kesi2-kesi1)
end
w=teta*y+(1-teta)*z
kesi3=dot(s,w)
D=(1/kesi3)*((w)*(w'))
E=(1/kesi2)*((z)*(z'))
H(:,:,n)=H(:,:,n-1)+D-E
else
Rd=x(1);Rb=x(2);bd=x(3);
dR=(Rd-Rb);
Td=4/3*pi*40000*dR^3+pi*3.8*45/0.001*dR^4+2*pi*Rd^2*bd*[15+0.1*(45*Rd/0.001)];
Td0=4/3*pi*15*dR^3+pi*0.1*45/0.001*dR^4+2*pi*Rd^2*bd*[15+0.1*(45*Rd/0.001)];
muJ=Rd^4*bd/[(0.069)^4*0.022-(0.065)^4*(bd+0.002)];
f=-Td;
g1=1-Td/(50*Td0);
g2=1-muJ/0.25;
g3=muJ/0.5-1;
g4=Rb/Rd-1;
g5=1-Rd/0.030;
g6=Rd/0.068-1;
g7=1-Rb/0.010;
g8=Rb/0.020-1;
g9=1-bd/0.003;
g10=bd/0.012-1;
gf=[(2307309447473265*(Rb - Rd)^3)/1073741824 - (17271088261788183*(Rb - Rd)^2)/34359738368 - 9000*pi*Rd^2*bd - 4*pi*Rd*bd*(4500*Rd + 15);
(17271088261788183*(Rb - Rd)^2)/34359738368 - (2307309447473265*(Rb - Rd)^3)/1073741824;
-2*pi*Rd^2*(4500*Rd + 15)]';
gg1=[(((2307309447473265*(Rb - Rd)^4)/4294967296 - (5757029420596061*(Rb - Rd)^3)/34359738368 + 2*pi*Rd^2*bd*(4500*Rd + 15))*(3000*pi*(Rb - Rd)^2 - 900000*pi*(Rb - Rd)^3 + 450000*pi*Rd^2*bd + 200*pi*Rd*bd*(4500*Rd + 15)))/(- 1000*pi*(Rb - Rd)^3 + 225000*pi*(Rb - Rd)^4 + 100*pi*Rd^2*bd*(4500*Rd + 15))^2 - ((17271088261788183*(Rb - Rd)^2)/34359738368 - (2307309447473265*(Rb - Rd)^3)/1073741824 + 9000*pi*Rd^2*bd + 4*pi*Rd*bd*(4500*Rd + 15))/(225000*pi*(Rb - Rd)^4 - 1000*pi*(Rb - Rd)^3 + 100*pi*Rd^2*bd*(4500*Rd + 15));
((17271088261788183*(Rb - Rd)^2)/34359738368 - (2307309447473265*(Rb - Rd)^3)/1073741824)/(225000*pi*(Rb - Rd)^4 - 1000*pi*(Rb - Rd)^3 + 100*pi*Rd^2*bd*(4500*Rd + 15)) - ((3000*pi*(Rb - Rd)^2 - 900000*pi*(Rb - Rd)^3)*((2307309447473265*(Rb - Rd)^4)/4294967296 - (5757029420596061*(Rb - Rd)^3)/34359738368 + 2*pi*Rd^2*bd*(4500*Rd + 15)))/(- 1000*pi*(Rb - Rd)^3 + 225000*pi*(Rb - Rd)^4 + 100*pi*Rd^2*bd*(4500*Rd + 15))^2;
(100*pi*Rd^2*(4500*Rd + 15)*((2307309447473265*(Rb - Rd)^4)/4294967296 - (5757029420596061*(Rb - Rd)^3)/34359738368 + 2*pi*Rd^2*bd*(4500*Rd + 15)))/(- 1000*pi*(Rb - Rd)^3 + 225000*pi*(Rb - Rd)^4 + 100*pi*Rd^2*bd*(4500*Rd + 15))^2 - (2*pi*Rd^2*(4500*Rd + 15))/(225000*pi*(Rb - Rd)^4 - 1000*pi*(Rb - Rd)^3 + 100*pi*Rd^2*bd*(4500*Rd + 15))]';
gg2=[(16*Rd^3*bd)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000);
0;
(4*Rd^4)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000) - (2634287287446093*Rd^4*bd)/(36893488147419103232*((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000)^2)]';
gg3=[-(8*Rd^3*bd)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000);
0;
(2634287287446093*Rd^4*bd)/(73786976294838206464*((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000)^2) - (2*Rd^4)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000)]';
gg4=[-Rb/Rd^2;
1/Rd;
0]';
gg5=[-100/3;
0;
0]';
gg6=[250/17;
0;
0]';
gg7=[0;
-100;
0]';
gg8=[0;
50;
0]';
gg9=[0;
0;
-1000/3]';
gg10=[0;
0;
250/3]';
a=[g1,g2,g3,g4,g5,g6,g7,g8,g9,g10];
V0=max([0 a]);
A=[gg1;gg2;gg3;gg4;gg5;gg6;gg7;gg8;gg9;gg10];
options = optimoptions('quadprog','Algorithm','interior-point-convex','Display','off');
[d,fval,exitflag,output,lambda] = quadprog(H(:,:,n),gf,A,-a',[],[],[],[],[],options);
dd=norm(d);
r0=sum(lambda.ineqlin);
b=[R0 r0];
R=max(b);
beta0=gama*dd^2;
phi0=f+R*V0;
phi=1;
phiii=0;
k=1;
xold=x;
while (phi>phiii)
t0=(1/2)^(k-1);
% x1=x0(1)+t0*d(1);
% x2=x0(2)+t0*d(2);
% x3=x0(3)+t0*d(3);
% x4=x0(4)+t0*d(4);
Rd=x(1)+t0*d(1);
Rb=x(2)+t0*d(2);
bd=x(3)+t0*d(3);
% x1=x(1);x2=x(2);x3=x(3);x4=x(4);
dR=(Rd-Rb);
Td=4/3*pi*40000*dR^3+pi*3.8*45/0.001*dR^4+2*pi*Rd^2*bd*[15+0.1*(45*Rd/0.001)];
Td0=4/3*pi*15*dR^3+pi*0.1*45/0.001*dR^4+2*pi*Rd^2*bd*[15+0.1*(45*Rd/0.001)];
muJ=Rd^4*bd/[(0.069)^4*0.022-(0.065)^4*(bd+0.002)];
f=-Td;
g1=1-Td/(50*Td0);
g2=1-muJ/0.25;
g3=muJ/0.5-1;
g4=Rb/Rd-1;
g5=1-Rd/0.030;
g6=Rd/0.068-1;
g7=1-Rb/0.010;
g8=Rb/0.020-1;
g9=1-bd/0.003;
g10=bd/0.012-1;
a=[g1,g2,g3,g4,g5,g6,g7,g8,g9,g10];
V=max([0 a]);
phi=f+R*V;
phiii=phi0-t0*beta0;
k=k+1;
end
alpha=t0;
x=[Rb Rd bd];
xnew=x;
R0=R;
n=n+1
hisx(n,:)=x;
dd
V
f
a
pause(1)
s=alpha*d;
z=H(:,:,n-1)*s;
Rd=xnew(1);Rb=xnew(2);bd=xnew(3);
u1=lambda.ineqlin(1);u2=lambda.ineqlin(2);u3=lambda.ineqlin(3);u4=lambda.ineqlin(4);
u5=lambda.ineqlin(5);u6=lambda.ineqlin(6);u7=lambda.ineqlin(7);u8=lambda.ineqlin(8);
u9=lambda.ineqlin(9);u10=lambda.ineqlin(10);
yy1=[(250*u6)/17 - (100*u5)/3 - u1*(((17271088261788183*(Rb - Rd)^2)/34359738368 - (2307309447473265*(Rb - Rd)^3)/1073741824 + 9000*pi*Rd^2*bd + 4*pi*Rd*bd*(4500*Rd + 15))/(225000*pi*(Rb - Rd)^4 - 1000*pi*(Rb - Rd)^3 + 100*pi*Rd^2*bd*(4500*Rd + 15)) - (((2307309447473265*(Rb - Rd)^4)/4294967296 - (5757029420596061*(Rb - Rd)^3)/34359738368 + 2*pi*Rd^2*bd*(4500*Rd + 15))*(3000*pi*(Rb - Rd)^2 - 900000*pi*(Rb - Rd)^3 + 450000*pi*Rd^2*bd + 200*pi*Rd*bd*(4500*Rd + 15)))/(- 1000*pi*(Rb - Rd)^3 + 225000*pi*(Rb - Rd)^4 + 100*pi*Rd^2*bd*(4500*Rd + 15))^2) - (17271088261788183*(Rb - Rd)^2)/34359738368 + (2307309447473265*(Rb - Rd)^3)/1073741824 - 9000*pi*Rd^2*bd - (Rb*u4)/Rd^2 + (16*Rd^3*bd*u2)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000) - (8*Rd^3*bd*u3)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000) - 4*pi*Rd*bd*(4500*Rd + 15);
50*u8 - 100*u7 + u4/Rd + u1*(((17271088261788183*(Rb - Rd)^2)/34359738368 - (2307309447473265*(Rb - Rd)^3)/1073741824)/(225000*pi*(Rb - Rd)^4 - 1000*pi*(Rb - Rd)^3 + 100*pi*Rd^2*bd*(4500*Rd + 15)) - ((3000*pi*(Rb - Rd)^2 - 900000*pi*(Rb - Rd)^3)*((2307309447473265*(Rb - Rd)^4)/4294967296 - (5757029420596061*(Rb - Rd)^3)/34359738368 + 2*pi*Rd^2*bd*(4500*Rd + 15)))/(- 1000*pi*(Rb - Rd)^3 + 225000*pi*(Rb - Rd)^4 + 100*pi*Rd^2*bd*(4500*Rd + 15))^2) + (17271088261788183*(Rb - Rd)^2)/34359738368 - (2307309447473265*(Rb - Rd)^3)/1073741824;
(250*u10)/3 - (1000*u9)/3 - u3*((2*Rd^4)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000) - (2634287287446093*Rd^4*bd)/(73786976294838206464*((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000)^2)) + u2*((4*Rd^4)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000) - (2634287287446093*Rd^4*bd)/(36893488147419103232*((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000)^2)) - u1*((2*pi*Rd^2*(4500*Rd + 15))/(225000*pi*(Rb - Rd)^4 - 1000*pi*(Rb - Rd)^3 + 100*pi*Rd^2*bd*(4500*Rd + 15)) - (100*pi*Rd^2*(4500*Rd + 15)*((2307309447473265*(Rb - Rd)^4)/4294967296 - (5757029420596061*(Rb - Rd)^3)/34359738368 + 2*pi*Rd^2*bd*(4500*Rd + 15)))/(- 1000*pi*(Rb - Rd)^3 + 225000*pi*(Rb - Rd)^4 + 100*pi*Rd^2*bd*(4500*Rd + 15))^2) - 2*pi*Rd^2*(4500*Rd + 15)];
Rd=xold(1);Rb=xold(2);bd=xold(3);
yy0=[(250*u6)/17 - (100*u5)/3 - u1*(((17271088261788183*(Rb - Rd)^2)/34359738368 - (2307309447473265*(Rb - Rd)^3)/1073741824 + 9000*pi*Rd^2*bd + 4*pi*Rd*bd*(4500*Rd + 15))/(225000*pi*(Rb - Rd)^4 - 1000*pi*(Rb - Rd)^3 + 100*pi*Rd^2*bd*(4500*Rd + 15)) - (((2307309447473265*(Rb - Rd)^4)/4294967296 - (5757029420596061*(Rb - Rd)^3)/34359738368 + 2*pi*Rd^2*bd*(4500*Rd + 15))*(3000*pi*(Rb - Rd)^2 - 900000*pi*(Rb - Rd)^3 + 450000*pi*Rd^2*bd + 200*pi*Rd*bd*(4500*Rd + 15)))/(- 1000*pi*(Rb - Rd)^3 + 225000*pi*(Rb - Rd)^4 + 100*pi*Rd^2*bd*(4500*Rd + 15))^2) - (17271088261788183*(Rb - Rd)^2)/34359738368 + (2307309447473265*(Rb - Rd)^3)/1073741824 - 9000*pi*Rd^2*bd - (Rb*u4)/Rd^2 + (16*Rd^3*bd*u2)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000) - (8*Rd^3*bd*u3)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000) - 4*pi*Rd*bd*(4500*Rd + 15);
50*u8 - 100*u7 + u4/Rd + u1*(((17271088261788183*(Rb - Rd)^2)/34359738368 - (2307309447473265*(Rb - Rd)^3)/1073741824)/(225000*pi*(Rb - Rd)^4 - 1000*pi*(Rb - Rd)^3 + 100*pi*Rd^2*bd*(4500*Rd + 15)) - ((3000*pi*(Rb - Rd)^2 - 900000*pi*(Rb - Rd)^3)*((2307309447473265*(Rb - Rd)^4)/4294967296 - (5757029420596061*(Rb - Rd)^3)/34359738368 + 2*pi*Rd^2*bd*(4500*Rd + 15)))/(- 1000*pi*(Rb - Rd)^3 + 225000*pi*(Rb - Rd)^4 + 100*pi*Rd^2*bd*(4500*Rd + 15))^2) + (17271088261788183*(Rb - Rd)^2)/34359738368 - (2307309447473265*(Rb - Rd)^3)/1073741824;
(250*u10)/3 - (1000*u9)/3 - u3*((2*Rd^4)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000) - (2634287287446093*Rd^4*bd)/(73786976294838206464*((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000)^2)) + u2*((4*Rd^4)/((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000) - (2634287287446093*Rd^4*bd)/(36893488147419103232*((2634287287446093*bd)/147573952589676412928 - 546584892005391387/1180591620717411303424000)^2)) - u1*((2*pi*Rd^2*(4500*Rd + 15))/(225000*pi*(Rb - Rd)^4 - 1000*pi*(Rb - Rd)^3 + 100*pi*Rd^2*bd*(4500*Rd + 15)) - (100*pi*Rd^2*(4500*Rd + 15)*((2307309447473265*(Rb - Rd)^4)/4294967296 - (5757029420596061*(Rb - Rd)^3)/34359738368 + 2*pi*Rd^2*bd*(4500*Rd + 15)))/(- 1000*pi*(Rb - Rd)^3 + 225000*pi*(Rb - Rd)^4 + 100*pi*Rd^2*bd*(4500*Rd + 15))^2) - 2*pi*Rd^2*(4500*Rd + 15)];
y=yy1-yy0;
kesi1=dot(s,y);
kesi2=dot(s,z);
if kesi1>=0.2*kesi2
teta=1;
else
teta=0.8*kesi2/(kesi2-kesi1);
end
w=teta*y+(1-teta)*z;
kesi3=dot(s,w);
D=(1/kesi3)*((w)*(w'));
E=(1/kesi2)*((z)*(z'));
H(:,:,n)=H(:,:,n-1)+D-E;
end
hisf(n,:)=f
end
disp('end')
Antworten (0)
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!