Error: Failure at t=1.776273e-02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (6.310588e-17) at time t.

1 vue (au cours des 30 derniers jours)
I get this error while running my program in ode23s
clc
clear all
close all
t1 = 10;
x0 = [10 0 0 0 0 0 0 0 0 0 0 0 0]';
% for i = 1:length(t1)
tspan = (0:0.1:t1);
[t,y] = ode23s(@(t,y) mysystem(t,y),tspan,x0);
u = y(:,1);
v = y(:,2);
w = y(:,3);
p = y(:,4);
q = y(:,5);
r = y(:,6);
e0 = y(:,7);
e1 = y(:,8);
e2 = y(:,9);
e3 = y(:,10);
phi = y(:,11);
theta = y(:,12);
si = y(:,13);
x0 = y(end,:);
% end
figure,plot(t,u,'r','linewidth',1.5),hold on,plot(t,v,'--','linewidth',1.5),hold on,plot(t,w,'b','linewidth',1.5)
legend('u','v','w'),grid on;
figure,plot(t,p,'r','linewidth',1.5),hold on,plot(t,q,'--','linewidth',1.5),hold on,plot(t,r,'b','linewidth',1.5)
legend('p','q','r'),grid on;
%%
function dydt =mysystem(t,y)
m = 86816;
l = 129.5;
d = 32;
b = 16;
Ix = 12470787;
Iy = 55472495;
Iz = 49248564;
Ixz = 1648250;
ax = 60.39;
ay = 0.01;
az = -10.43;
k1 = 0.085;
k2 = 0.865;
k_dash = 0.61;
Xudot = -k1*m;
Yvdot = -k2*m;
Zwdot = Yvdot;
Iy_bar = m*(l^2+d^2)/20;
Mqdot = -k_dash*Iy_bar;
Nrdot = Mqdot;
mx = m-Xudot;
my = m-Yvdot;
mz = m-Zwdot;
my = mz;
Jx = Ix;
Jy = Iy-Mqdot;
Jz = Iz-Nrdot;
Jxz = Ixz;
M = [mx 0 0 0 m*az 0;0 my 0 -m*az 0 m*ax;0 0 mz 0 -m*ax 0;0 -m*az 0 Jx 0 -Jxz;m*az 0 -m*ax 0 Jy 0;0 m*ax 0 -Jxz 0 Jz];
mu = 0;
V = 70870.2;
% S = V^(2/3);
S = 1;
row = 1.225;
U0 = 25;
a1 = 60.23;
a2 = 69.31;
lf1 = 109.1;
lf2 = 113.7;
lf3 = 13.8;
lh = 103.6;
lgx = 45.3;
lgz = 17.6;
xcv = 63.47;
Sh = 1689.3;
Sf = 709.7;
Sg = 24.7;
etaf = 0.3789;
etak = 1.0518;
CDh0 = 0.0250;
CDf0 = 0.006;
CDg0 = 0.01;
CDch = 0.50;
CDcf = 1;
CDcg = 1;
Ctf = 0;
doalpha = 5.73;
dodelta = 1.24;
f = (lh-a1)/a2;
I1 = pi*((b^2)/(V^(2/3)))*(1-f^2);
I3 = pi*((b^2)/(3*l*V^(2/3)))*(a1-2*a2*f^3-3*a1*f^2)-xcv*I1/l;
J1 = (b/(2*V^(2/3)))*(pi*a1+a2*(sqrt(1-f^2))*f+2*a2*asin(f));
J2 = (J1/l)*(a1-xcv)+((b*2)/(3*l*V^(2/3)))*(a2^2-a1^2-a2^2*((1-f^2)^(2/3)));
Cx1 = -(CDh0*Sh+CDf0*Sf+CDg0*Sg);
Cx2 = (k2-k1)*etak*I1*Sh;
Cx3 = Ctf*Sf;
Cy1 = -Cx2;
Cy2 = -0.5*doalpha*Sf*etaf;
Cy3 = -(CDch*J1*Sh+CDcf*Sf+CDcg*Sg);
Cy4 = -0.5*dodelta*Sf*etaf;
Cz1 = -Cx2;
Cz2 = Cy2;
Cz3 = -(CDch*J1*Sh+CDcf*Sf);
Cz4 = Cy4;
Cl1 = dodelta*Sf*etaf*lf3;
Cl2 = CDcg*Sg*lgz;
Cl3 = CDcg*Sg*lgz*l^2;
Cl4 = -(CDcg*Sg*lgz+CDcf*Sf*lf3)*d^2;
Cm1 = -(k2-k1)*etak*I3*Sh*l;
Cm2 = -0.5*doalpha*Sf*etaf*lf1;
Cm3 = -(CDch*J2*Sh*l+CDcf*Sf*lf2);
Cm4 = -0.5*dodelta*Sf*etaf*lf1;
Cm5 = -CDcf*Sf*lf2*l^2;
Cn1 = -Cm1;
Cn2 = -Cm2;
Cn3 = CDch*J2*Sh*l+CDcf*Sf*lf2+CDcg*Sg*lgx;
Cn4 = -Cm4;
Cn5 = -(CDcf*Sf*lf2+CDcg*Sg*lgx)*l^2;
g = 9.81;
B = V*row*g;
H = 1;
W = B+H*g;
dx = 3.5;
dy = 7.87;
dz = 21;
%% inputs
Tds = 1;
Tdp = 0;
% Tds = [ones(30,1)*1000;ones(30,1)*0;ones(30,1)*0];
% Tdp = [ones(30,1)*1000;ones(30,1)*0;ones(30,1)*0];
deltart = 0;%[ones(Tend,1)*10*pi/180]*0;
deltarb = deltart;
deltael = 0;%zeros(Tend,1);
deltaer = deltael;
% Vtotal = sqrt(y(1)^2+y(2)^2+y(3)^2);
% alpha = atan(y(3)/y(1));
% beta = asin(y(2)/Vtotal);
%
% Vtotal = sqrt(y(1)^2+y(2)^2+y(3)^2);
% alpha = atan(y(3)/y(1));
% beta = asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2)));
dydt = zeros(13,1);
% for i = 1:90
% fval(1) = -mz*y(3)*y(5)+my*y(6)*y(2)+m*(ax*(y(5)^2+y(6)^2)-az*y(6)*y(4))+0.5*row*U0^2*S*(Cx1*(cos(alpha))^2*(cos(beta))^2+Cx2*(sin(2*alpha)*sin(alpha/2)+sin(2*beta)*sin(beta/2)+Cx3))+(Tds(i)+Tdp(i))*cos(mu)+DCM(3,1)*(W-B);
% fval(2) = -mx*y(1)*y(6)+mz*y(4)*y(3)+m*(-ax*y(4)*y(5)-az*y(6)*y(5))+0.5*row*U0^2*S*(Cy1*cos(beta/2)*sin(2*beta)+Cy2*sin(2*beta)+Cy3*sin(beta)*sin(abs(beta))+Cy4*(deltart+deltarb))+DCM(3,2)*(W-B);
% fval(3) = -my*y(2)*y(4)+mx*y(5)*y(1)+m*(-ax*y(6)*y(4)+az*(y(5)^2+y(4)^2))+0.5*row*U0^2*S*(Cz1*cos(alpha/2)*sin(2*alpha)+Cz2*sin(2*alpha)+Cz3*sin(alpha)*sin(abs(alpha))+Cz4*(deltael+deltaer))-(Tdp(i)+Tds(i))*sin(mu)+DCM(3,3)*(W-B);
% fval(4) = -y(6)*y(5)*(Jz-Jy)+Jxz*y(4)*y(5)+m*az*(y(1)*y(6)-y(4)*y(3))+0.5*row*U0^2*S*l*(Cl1*(deltael-deltaer+deltarb-deltart)+Cl2*(sin(beta)*sin(abs(beta)))+0.5*row*S*(Cl3*y(6)*abs(y(6))+Cl4*y(4)*abs(y(4))))+(Tdp(i)-Tds(i))*sin(mu)*dy-DCM(3,2)*az*W;
% fval(5) = -y(4)*y(6)*(Jx-Jz)+Jxz*(y(6)^2-y(4)^2)+m*(ax*(y(2)*y(4)-y(5)*y(1))-az*(y(3)*y(5)-y(6)*y(2)))+0.5*row*U0^2*S*l*(Cm1*cos(alpha/2)*sin(2*alpha)+Cm2*sin(2*alpha)+Cm3*sin(alpha)*sin(abs(alpha))+Cm4*(deltael+deltaer)+0.5*row*S*Cm5*y(5)*abs(y(5)))+(Tds(i)+Tdp(i))*(dz*cos(mu)-dx*sin(mu))+(DCM(3,1)*az-DCM(3,3)*ax)*W;
% fval(6) = -y(5)*y(4)*(Jy-Jx)-Jxz*y(5)*y(6)+m*(-ax*(y(1)*y(6)-y(4)*y(3)))+0.5*row*U0^2*S*l*(Cn1*cos(beta/2)*sin(2*beta)+Cn2*sin(2*beta)+Cn3*sin(beta)*sin(abs(beta))+Cn4*(deltart+deltarb)+0.5*row*S*Cm5*y(6)*abs(y(6)))+(Tdp(i)-Tds(i))*cos(mu)*dy+DCM(3,2)*ax*W;
% end
dydt(1) = -mz*y(3)*y(5)+my*y(6)*y(2)+m*(ax*(y(5)^2+y(6)^2)-az*y(6)*y(4))+0.5*row*U0^2*S*(Cx1*(cos(atan(y(3)/y(1))))^2*(cos(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2)))))^2+Cx2*(sin(2*(atan(y(3)/y(1))))*sin((atan(y(3)/y(1)))/2)+sin(2*(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2)))))*sin((asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2))))/2)+Cx3))+(Tds+Tdp)*cos(mu)+2*(y(8)*y(10)-y(7)*y(9))*(W-B);
dydt(2) = -mx*y(1)*y(6)+mz*y(4)*y(3)+m*(-ax*y(4)*y(5)-az*y(6)*y(5))+0.5*row*U0^2*S*(Cy1*cos((asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2))))/2)*sin(2*(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2)))))+Cy2*sin(2*(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2)))))+Cy3*sin(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2))))*sin(abs(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2)))))+Cy4*(deltart+deltarb))+2*(y(7)*y(8)-y(9)*y(10))*(W-B);
dydt(3) = -my*y(2)*y(4)+mx*y(5)*y(1)+m*(-ax*y(6)*y(4)+az*(y(5)^2+y(4)^2))+0.5*row*U0^2*S*(Cz1*cos((atan(y(3)/y(1)))/2)*sin(2*(atan(y(3)/y(1))))+Cz2*sin(2*(atan(y(3)/y(1))))+Cz3*sin(atan(y(3)/y(1)))*sin(abs(atan(y(3)/y(1))))+Cz4*(deltael+deltaer))-(Tdp+Tds)*sin(mu)+((y(7)^2)-(y(8)^2)-(y(9)^2)+(y(10)^2))*(W-B);
dydt(4) = -y(6)*y(5)*(Jz-Jy)+Jxz*y(4)*y(5)+m*az*(y(1)*y(6)-y(4)*y(3))+0.5*row*U0^2*S*l*(Cl1*(deltael-deltaer+deltarb-deltart)+Cl2*(sin(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2))))*sin(abs(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2))))))+0.5*row*S*(Cl3*y(6)*abs(y(6))+Cl4*y(4)*abs(y(4))))+(Tdp-Tds)*sin(mu)*dy-2*(y(7)*y(8)-y(9)*y(10))*az*W;
dydt(5) = -y(4)*y(6)*(Jx-Jz)+Jxz*(y(6)^2-y(4)^2)+m*(ax*(y(2)*y(4)-y(5)*y(1))-az*(y(3)*y(5)-y(6)*y(2)))+0.5*row*U0^2*S*l*(Cm1*cos((atan(y(3)/y(1)))/2)*sin(2*(atan(y(3)/y(1))))+Cm2*sin(2*(atan(y(3)/y(1))))+Cm3*sin(atan(y(3)/y(1)))*sin(abs(atan(y(3)/y(1))))+Cm4*(deltael+deltaer)+0.5*row*S*Cm5*y(5)*abs(y(5)))+(Tds+Tdp)*(dz*cos(mu)-dx*sin(mu))+2*(y(8)*y(10)-y(7)*y(9))*az-((y(7)^2)-(y(8)^2)-(y(9)^2)+(y(10)^2))*ax*W;
dydt(6) = -y(5)*y(4)*(Jy-Jx)-Jxz*y(5)*y(6)+m*(-ax*(y(1)*y(6)-y(4)*y(3)))+0.5*row*U0^2*S*l*(Cn1*cos((asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2))))/2)*sin(2*(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2)))))+Cn2*sin(2*(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2)))))+Cn3*sin(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2))))*sin(abs(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2)))))+Cn4*(deltart+deltarb)+0.5*row*S*Cm5*y(6)*abs(y(6)))+(Tdp-Tds)*cos(mu)*dy+2*(y(7)*y(8)-y(9)*y(10))*ax*W;
dydt(7) = 0.5*(-y(4)*y(8)-y(5)*y(9)-y(6)*y(10));
dydt(8) = 0.5*(y(4)*y(7)-y(5)*y(10)+y(6)*y(9));
dydt(9) = 0.5*(y(4)*y(10)+y(5)*y(7)-y(6)*y(8));
dydt(10) = 0.5*(-y(4)*y(9)+y(5)*y(8)+y(6)*y(7));
dydt(11) = y(4)+y(5)*tan(y(12))*sin(y(11))+y(6)*tan(y(12))*cos(y(12));
dydt(12) = y(5)*cos(y(11))-y(6)*sin(y(11));
dydt(13) = y(5)*sin(y(11))*sec(y(12))+y(6)*cos(y(11))*sec(y(12));
% fval(7) = 0.5*(-y(4)*y(8)-y(5)*y(9)-y(6)*y(10));
% fval(8) = 0.5*(y(4)*y(7)-y(5)*y(10)+y(6)*y(9));
% fval(9) = 0.5*(y(4)*y(10)+y(5)*y(7)-y(6)*y(8));
% fval(10) = 0.5*(-y(4)*y(9)+y(5)*y(8)+y(6)*y(7));
dydt(1:6) = M\dydt(1:6);
end

Réponses (1)

Shadaab Siddiqie
Shadaab Siddiqie le 15 Avr 2021
From my understanding you are getting an error while running above code.
  • This is because MATLAB is trying to reduce the time step to a really small amount in order to counter the abrupt change due to the discontinuity in the reference signal.
  • For sharp discontinuities,it might not be possible to avoid this warning. However for non discontinous inputs we can set relative and absolute tolerance to a higher number than the default setting.
We can set the tolerances to a higher value for example by using the following commands:
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);
[t,y] = ode23s(@objectivefunction,tspan,y0,options);

Catégories

En savoir plus sur Programming dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by