problems with integration tolerances when using ode15s
显示 更早的评论
I am working on simulating a dynamic system using a Concurrent Learning controller, but for now all I am doing is plotting the norm of my errors (e and r) and the norm of one of my estimates (utilde). When I run my code I get a yellow warning stating "Warning: Failure at t=1.879973e+01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (5.684342e-14) at time t." I'm not sure if this is due to my reltol and abstol values or a problem within my code but it is effecting how my plots should look. The three plots should show convergence to zero, but norm utilde(t) plots strangely while norm e and norm r seem to plot normally. My code is posted below. Any help of suggestions are greatly appreciated.
clc;
clear;
close all;
low=0;%stated that all constants mx,mphi,cx,cphi,1,ax,and aphi were greater than zero
high=1;
Amx=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 1 for constant mx
Amphi=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 1 for constant mphi
Acx=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 1 for constant cx
Acphi=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 1 for constant cphi
Al=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 1 for constant l
Aax=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 1 for constant ax
Aaphi=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 1 for constant aphi
%Feedback terms
lowc=0;
highc=5;
ABx=(highc-lowc).*rand(100,1) + lowc;%column vector of 100 uniformly distributed decimals between 0 and 5 for constant Bx
ABphi=(highc-lowc).*rand(100,1) + lowc;%column vector of 100 uniformly distributed decimals between 0 and 5 for constant Bphi
Aalphax=(highc-lowc).*rand(100,1) + lowc;%column vector of 100 uniformly distributed decimals between 0 and 5 for constant alphax
Aalphaphi=(highc-lowc).*rand(100,1) + lowc;%column vector of 100 uniformly distributed decimals between 0 and 5 for constant alphaphi
AJx=(highc-lowc).*rand(100,1) + lowc;%column vector of 100 uniformly distributed decimals between 0 and 5 for constant Jx
AJphi=(highc-lowc).*rand(100,1) + lowc;%column vector of 100 uniformly distributed decimals between 0 and 5 for constant Jphi
%Creating xd and phid equations and their derivatives
lowxd=-1;%left limit of xd set to 1m
highxd=1;%right limit of xd set to 1m
Axdbar=(highxd-lowxd).*rand(100,1)+lowxd;%column vector of 100 uniformly distributed decimals between -1 and 1 for variable xdbar
Ax=(highxd-lowxd).*rand(100,1)+lowxd;%column vector of 100 uniformly distributed decimals between -1 and 1 for variable x
lowfreq=0;%lowest frequency set to 0 Hz
highfreq=0.5;%highest frequency set to 0.5 Hz
Afxd=(highfreq-lowfreq).*rand(100,1)+lowfreq;%column vector of 100 uniformlay distributed decimals between 0 and 0.5 Hz for variable fxd
lowphi=-1*(pi/2);
highphi=pi/2;
Aphidbar=(highphi-lowphi).*rand(100,1)+lowphi;%column vector of 100 uniformly distributed decimals between -pi/2 and pi/2 for variable phidbar
Aphi=(highphi-lowphi).*rand(100,1)+lowphi;%column vector of 100 uniformly distributed decimals between -pi/2 and pi/2 for variable phi
Afphid=(highfreq-lowfreq).*rand(100,1)+lowfreq;%column vector of 100 uniformlay distributed decimals between 0 and 0.5 Hz for variable fphid
for ii=1:100 %signifies 100 iterations of each constant and variable
%Constants
mx=Amx(ii);
mphi=Amphi(ii);
cx=Acx(ii);
cphi=Acphi(ii);
l=Al(ii);
ax=Aax(ii);
aphi=Aaphi(ii);
%Trajectory Variables
xdbar=Axdbar(ii);
x=Ax(ii);
fxd=Afxd(ii);
phidbar=Aphidbar(ii);
phi=Aphi(ii);
fphid=Afphid(ii);
g=9.8;%gravitational constant
%Feedback Terms
Bx=ABx(ii);
Bphi=ABphi(ii);
alphax=Aalphax(ii);
alphaphi=Aalphaphi(ii);
Jx=AJx(ii);
Jphi=AJphi(ii);
gammat=eye(6);
gammaa=eye(2);
Ktheta=eye(6);
Ka=eye(2);
tspan= [0 50];%time span for elapsed simulation
options = odeset('reltol', 1e-3, 'abstol', 1e-5 );%default options for relative and absolute tolerance. The tolerances are used to limit the local discretization error during integration
x0=[x 0 0 phi 0 0 0 0 0 0 0 0 0 0 0];
[t, x] = ode15s(@(t,x) dynamics(t,x,phi,mx,mphi,cx,cphi,l,ax,aphi,Bx,Bphi,alphax,alphaphi,Jx,Jphi,xdbar,fxd,phidbar,fphid,g,gammat,gammaa,Ktheta,Ka), tspan, x0, options);
t=transpose(t);
x=transpose(x);
x1=x(1,:);
xdot=x(2,:);
xdotdot=x(3,:);
phi1=x(4,:);
phidot=x(5,:);
ux=x(6,:);
uphi=x(7,:);
thetahat=x(8:13,:);
ahat=x(14:15,:);
%Trajectory equations
xd=xdbar*sin(2*pi*fxd.*t);%xd(t)
xddot=2*pi*fxd*xdbar*cos(2*pi*fxd.*t);%xd(t)dot
xddotdot=-1*(2*pi*fxd)^(2)*xdbar*sin(2*pi*fxd.*t);%xd(t)dotdot
phid=phidbar*sin(2*pi*fphid.*t);%phid(t)
phiddot=2*pi*fphid*phidbar*cos(2*pi*fphid.*t);%phid(t)dot
phiddotdot=-1*(2*pi*fphid)^(2)*phidbar*sin(2*pi*fphid.*t);%phid(t)dotdot
%Estimation Errors
sizet=size(t,2);
theta=[mx*ones(1,sizet); mphi*ones(1,sizet); mphi*l*ones(1,sizet); mphi*(l)^(2)*ones(1,sizet); cx*ones(1,sizet); cphi*ones(1,sizet)];
a=[ax*ones(1,sizet); aphi*ones(1,sizet)];
thetatilde=theta-thetahat;
atilde=a-ahat;
%Tracking Errors
ex=xd-x1;
exdot=xddot-xdot;
exdotdot=xddotdot-xdotdot;
ephi=phid-phi1;
ephidot=phiddot-phidot;
ephidotdot=phiddotdot-((-1*cphi.*phidot-mphi*l*g*sin(phi)-mphi*l*cos(phi).*xdotdot+uphi)/(mphi*l^(2)));
%Filtered Tracking Errors
rx=exdot+alphax.*ex;
rxdot=exdotdot+alphax.*exdot;
rphi=ephidot+alphaphi.*ephi;
rphidot=ephidotdot+alphaphi.*ephidot;
%Design Inputs
yxthetahat=zeros(1,length(t));
yphithetahat=zeros(1,length(t));
Yxthetahat=zeros(1,length(t));
Yphithetahat=zeros(1,length(t));
yxathetahatdot=zeros(1,length(t));
yphiathetahatdot=zeros(1,length(t));
Yxthetahatdot=zeros(1,length(t));
Yphithetahatdot=zeros(1,length(t));
for jj=1:length(t)
yx=[xdotdot(1,jj) (sin(phi))^(2)*xdotdot(1,jj)-g*cos(phi)*sin(phi) phidot(1,jj)^(2)*sin(phi) 0 xdot(1,jj) 0];
yphi=[0 0 g*sin(phi)+cos(phi)*xdotdot(1,jj) (-1*cphi*phidot(1,jj)-mphi*l*g*sin(phi)-mphi*l*cos(phi)*xdotdot(1,jj)+uphi(1,jj))/(mphi*l^(2)) 0 phidot(1,jj)];
Yx=[xddotdot(1,jj)+alphax*exdot(1,jj) sin(phi)^(2)*xdotdot(1,jj)-g*cos(phi)*sin(phi)+sin(phi)^(2)*alphax*exdot(1,jj)+sin(phi)*cos(phi)*phidot(1,jj)*rx(1,jj) phidot(1,jj)^(2)*sin(phi) 0 xdot(1,jj) 0];
Yphi=[0 0 g*sin(phi)+cos(phi)*xdotdot(1,jj) phiddotdot(1,jj)+alphaphi*ephidot(1,jj) 0 phidot(1,jj)];
yxa=[ux(1,jj) 0];
yphia=[0 uphi(1,jj)];
% thetahatdot=gammat*Yx.'*rx+gammat*Yphi.'*rphi+gammat*Ktheta*(yx.'*x(6)-yx.'*yx*x(8:13))+gammat*Ktheta*(yphi.'*x(7)-yphi.'*yphi*x(8:13));
yxthetahat(jj)=yx*thetahat(:,jj);
yphithetahat(jj)=yphi*thetahat(:,jj);
Yxthetahat(jj)=Yx*thetahat(:,jj);
Yphithetahat(jj)=Yphi*thetahat(:,jj);
% Yxthetahatdot(jj)=Yx*thetahatdot
end
udx=Yxthetahat+ex+Bx.*rx;
udphi=Yphithetahat+ephi+Bphi.*rphi;
utildex=udx-ux;
utildephi=udphi-uphi;
% udxdot=Yx*thetahatdot+exdot+Bx*rxdot;
% udphidot=Yphi*thetahatdot+ephidot+Bphi*rphidot;
% mux=udxdot+x(14)*x(6)+rx+Jx*utildex;
% muphi=udphidot+x(15)*x(7)+rphi+Jphi*utildephi;
% ahatdot=gammaa*utildex*x(6)+gammaa*utildephi*x(7)+gammaa*Ka*(yxa.'*(mux-(-1*ax*x(6)+mux))-yxa.'*yxa*x(14:15))+gammaa*Ka*(yphia.'*(muphi-(-1*aphi*x(7)+muphi))-yphia.'*yphia*x(14:15));
%
%Tracking Error Plot
norme=sqrt(ex.^2+ephi.^2);
%Filtered Tracking Error Plot
normr=sqrt(rx.^2+rphi.^2);
%Estimation Errors
normutilde=sqrt(utildex.^2+utildephi.^2);
figure(1)
plot(t,norme)
xlabel('Time')
ylabel('norm e(t)')
hold on
figure(2)
plot(t,normr)
xlabel('Time')
ylabel('norm r(t)')
hold on
figure(3)
plot(t, normutilde)
xlabel('Time')
ylabel('norm utilde(t)')
hold on
% figure(4)
% plot(t,normthetatilde)
% xlabel('Time')
% ylabel('norm thetatilde(t)')
% hold on
% figure(5)
% plot(t,normatilde)
% xlabel('Time')
% ylabel('norm atilde(t)')
% hold on
end
hold off
function dstate = dynamics(t,x,phi,mx,mphi,cx,cphi,l,ax,aphi,Bx,Bphi,alphax,alphaphi,Jx,Jphi,xdbar,fxd,phidbar,fphid,g,gammat,gammaa,Ktheta,Ka)
%Trajectory equations
xd=xdbar*sin(2*pi*fxd*t);%xd(t)
xddot=2*pi*fxd*xdbar*cos(2*pi*fxd*t);%xd(t)dot
xddotdot=-1*(2*pi*fxd)^(2)*xdbar*sin(2*pi*fxd*t);%xd(t)dotdot
phid=phidbar*sin(2*pi*fphid*t);%phid(t)
phiddot=2*pi*fphid*phidbar*cos(2*pi*fphid*t);%phid(t)dot
phiddotdot=-1*(2*pi*fphid)^(2)*phidbar*sin(2*pi*fphid*t);%phid(t)dotdot
%Estimation Errors
theta=[mx;mphi;mphi*l;mphi*(l)^(2);cx;cphi];
a=[ax;aphi];
thetatilde=theta-x(8:13);
atilde=a-x(14:15);
%Tracking Errors
ex=xd-x(1);
exdot=xddot-x(2);
exdotdot=xddotdot-x(3);
ephi=phid-x(4);
ephidot=phiddot-x(5);
ephidotdot=phiddotdot-((-1*cphi*x(5)-mphi*l*g*sin(phi)-mphi*l*cos(phi)*x(3)+x(7))/(mphi*l^(2)));
%Filtered Tracking Errors
rx=exdot+alphax*ex;
rxdot=exdotdot+alphax*exdot;
rphi=ephidot+alphaphi*ephi;
rphidot=ephidotdot+alphaphi*ephidot;
%Design Inputs
yx=[x(3) (sin(phi))^(2)*x(3)-g*cos(phi)*sin(phi) x(5)^(2)*sin(phi) 0 x(2) 0];
yphi=[0 0 g*sin(phi)+cos(phi)*x(3) (-1*cphi*x(5)-mphi*l*g*sin(phi)-mphi*l*cos(phi)*x(3)+x(7))/(mphi*l^(2)) 0 x(5)];
Yx=[xddotdot+alphax*exdot sin(phi)^(2)*x(3)-g*cos(phi)*sin(phi)+sin(phi)^(2)*alphax*exdot+sin(phi)*cos(phi)*x(5)*rx x(5)^(2)*sin(phi) 0 x(2) 0];
Yphi=[0 0 g*sin(phi)+cos(phi)*x(3) phiddotdot+alphaphi*ephidot 0 x(5)];
udx=Yx*x(8:13)+ex+Bx*rx;
udphi=Yphi*x(8:13)+ephi+Bphi*rphi;
yxa=[x(6) 0];
yphia=[0 x(7)];
utildex=udx-x(6);
utildephi=udphi-x(7);
thetahatdot=gammat*Yx.'*rx+gammat*Yphi.'*rphi+gammat*Ktheta*(yx.'*x(6)-yx.'*yx*x(8:13))+gammat*Ktheta*(yphi.'*x(7)-yphi.'*yphi*x(8:13));
udxdot=Yx*thetahatdot+exdot+Bx*rxdot;
udphidot=Yphi*thetahatdot+ephidot+Bphi*rphidot;
mux=udxdot+x(14)*x(6)+rx+Jx*utildex;
muphi=udphidot+x(15)*x(7)+rphi+Jphi*utildephi;
ahatdot=gammaa*utildex*yxa.'+gammaa*utildephi*yphia.'+gammaa*Ka*(yxa.'*(mux-(-1*ax*x(6)+mux))-yxa.'*yxa*x(14:15))+gammaa*Ka*(yphia.'*(muphi-(-1*aphi*x(7)+muphi))-yphia.'*yphia*x(14:15));
%Define State Variables
dstate = zeros(15,1);%initialize system matrix
dstate(1)=x(2);%xdot
dstate(2)=x(3);%xdotdot
dstate(3)=(-1*cx*x(3)+mphi*g*x(5)*(cos(phi)^(2)-sin(phi)^(2))+mphi*l*x(5)*(2*((-1*cphi*x(5)-mphi*l*g*sin(phi)-mphi*l*cos(phi)*x(3)+x(7))/(mphi*l^(2)))*sin(phi)+x(5)^(2)*cos(phi))+(-1*ax*x(6)+(mux))-2*mphi*cos(phi)*x(5)*x(3))/(mx+mphi*sin(phi)^(2));%xdotdotdot
dstate(4)=x(5);%phidot
dstate(5)=(-1*cphi*x(5)-mphi*l*g*sin(phi)-mphi*l*cos(phi)*x(3)+x(7))/(mphi*l^(2));%phidotdot
dstate(6)=-1*ax*x(6)+(mux);%uxdot
dstate(7)=-1*aphi*x(7)+(muphi);%uphidot
dstate(8:13)=thetahatdot;%thetahatdot
dstate(14:15)=ahatdot;%ahatdot
end
采纳的回答
Sulaymon Eshkabilov
2021-11-9
0 个投票
The answer lies in magnitudes of the calculated (integration) values of the variables that don't allow to have the set tolerances without reducing the time step. Therefore, it is better to loosen the tolerances, i.e. reltol and abstol. As it is now, the tolerances are quite significantly tight.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Mathematics 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
