The problem with xddotdot and exdot is because they have a time term in them which is what is making them huge column vectors. I am evaluating those terms at every instance in time from the ODE but my code right now is multipllying all of t for the whole ODE. So how can I multiply by an indice of time not all time?
Issues with Horzcat for defining a variable as a row vector
2 views (last 30 days)
Show older comments
Matthew Tortorella
on 5 Nov 2021
Commented: Walter Roberson
on 6 Nov 2021
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 an error stating "Dimensions of arrays being concatenated are not consistent." for Yx. When I look through that line I believe its because of a dimensions issue with my xddotdot and exdot terms. They are huge column vectors when they should just be a scalar because everything going into them is a scalar. My code is pasted below. Any suggestions would be much appreciated. Thanks.
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
Agammat=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 1 for constant gamma theta
Agammaa=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 1 for constant gamma a
%CL Constants
lowCL=0;
highCL=1;
AKtheta=(highCL-lowCL).*rand(100,1) + lowCL;%column vector of 100 uniformly distributed decimals between 0 and 1 for constant Ktheta
AKCLtheta=(highCL-lowCL).*rand(100,1) + lowCL;%column vector of 100 uniformly distributed decimals between 0 and 1 for constant KCLtheta
AKa=(highCL-lowCL).*rand(100,1) + lowCL;%column vector of 100 uniformly distributed decimals between 0 and 1 for constant Ka
AKCLa=(highCL-lowCL).*rand(100,1) + lowCL;%column vector of 100 uniformly distributed decimals between 0 and 1 for constant KCLa
%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);
gammat=Agammat(ii);
gammaa=Agammaa(ii);
%CL Constants
Ktheta=AKtheta(ii);
KCLtheta=AKCLtheta(ii);
Ka=AKa(ii);
KCLa=AKCLa(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);
tspan= [0 100];%time span for elapsed simulation
options = odeset('reltol', 1e-5, 'abstol', 1e-8 );%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), tspan, x0, options);
xdot=x(2);
xdotdot=x(3);
phidot=x(5);
ux=x(6);
uphi=x(7);
thetahat=x(8:13);
ahatx=x(14);
ahatphi=x(15);
ahat=[ahatx;ahatphi];
%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-thetahat;
atilde=a-ahat;
%Tracking Errors
ex=xd-x;
exdot=xddot-xdot;
exdotdot=xddotdot-xdotdot;
ephi=phid-phi;
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
yx=[xdotdot (sin(phi))^(2)*xdotdot-g*cos(phi)*sin(phi) phidot^(2)*sin(phi) 0 xdot 0];
yphi=[0 0 g*sin(phi)+cos(phi)*xdotdot (-1*cphi*phidot-mphi*l*g*sin(phi)-mphi*l*cos(phi)*xdotdot+uphi)/(mphi*l^(2)) 0 phidot];
Yx=[xddotdot+alphax*exdot sin(phi)^(2)*xdotdot-g*cos(phi)*sin(phi)+sin(phi)^(2)*alphax*exdot+sin(phi)*cos(phi)*phidot*rx phidot^(2)*sin(phi) 0 xdot 0];
Yphi=[0 0 g*sin(phi)+cos(phi)*xdotdot phiddotdot+alphaphi*ephidot 0 phidot];
udx=Yx*thetahat+ex+Bx*rx;
udphi=Yphi*thetahat+ephi+Bphi*rphi;
% udxdot=Yx*x(8:13)+exdot+Bx*rxdot;
% udphidot=Yphi*x(8:13)+ephidot+Bphi*rphidot;
%
utildex=udx-ux;
utildephi=udphi-uphi;
% mux=udxdot+ahatx*ux+rx+Jx*utildex;
% muphi=udphidot+ahatphi*uphi+rphi+Jphi*utildephi;
% xdotdotdot=(-1*cx*xdotdot+mphi*g*phidot*(cos(phi)^(2)-sin(phi)^(2))+mphi*l*phidot*(2*((-1*cphi*phidot-mphi*l*g*sin(phi)-mphi*l*cos(phi)*xdotdot+uphi)/(mphi*l^(2)))*sin(phi)+phidot^(2)*cos(phi))+(-1*ax*ux+(mux))-2*mphi*cos(phi)*phidot*xdotdot)/(mx+mphi*sin(phi)^(2));
%
%Tracking Error Plot
norme=zeros(1,length(t));
for jj=1:length(t)
norme(jj)= norm([ex(jj) ephi(jj)]);
end
%Filtered Tracking Error Plot
normr=zeros(1,length(t));
for kk=1:length(t)
normr(kk)= norm([rx(kk) rphi(kk)]);
end
%Estimation Errors
normutilde=zeros(1,length(t));
for aa=1:length(t)
normutilde(aa)= norm([utildex(aa) utildephi(aa)]);
end
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;
udxdot=Yx*x(8:13)+exdot+Bx*rxdot;
udphidot=Yphi*x(8:13)+ephidot+Bphi*rphidot;
utildex=udx-x(6);
utildephi=udphi-x(7);
mux=udxdot+x(14)*x(6)+rx+Jx*utildex;
muphi=udphidot+x(15)*x(7)+rphi+Jphi*utildephi;
% yxa=[x(6) 0];
% yphia=[0 x(7)];
%
% thetahatdot=gammat.*transpose(Yx).*rx+gammat.*transpose(Yphi).*rphi+gammat*Ktheta.*(transpose(yx)*x(6)-transpose(yx)*yx*x(8:13))+gammat*Ktheta.*(transpose(yphi)*x(7)-transpose(yphi)*yphi*x(8:13));
% ahatdot=gammaa*utildex*x(6)+gamma*utildephi*x(7)+gammaa*Ka.*(transpose(yxa)*(mux-(-1*ax*x(6)+mux))-transpose(yxa)*yxa*x(14:15))+gammaa*Ka.*(transpose(yphia)*(muphi-(-1*aphi*x(7)+muphi))-transpose(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
3 Comments
Walter Roberson
on 5 Nov 2021
No, the bug is in
xdot=x(2);
xdotdot=x(3);
phidot=x(5);
ux=x(6);
uphi=x(7);
thetahat=x(8:13);
ahatx=x(14);
ahatphi=x(15);
You should be extracting all columns of x for each of those, to get columns that are the same size as the t -- they should all be long column vectors, the same number of rows as t has. (Be careful with thetahat that you get the right number of columns.)
What you have now is definitely wrong. Remember that x(2) for a 2D array, is the same as x(2,1), which is the first x output corresponding to the second time output.
If you are trying to determine some kind of "settled" value, then you should be extracting x(end,2), x(end,3) and so on, and probably t(end) to go along with that.
Accepted Answer
Sulaymon Eshkabilov
on 5 Nov 2021
You had better use this vectorized approach instead of a loop based calculation of norm values, i.e. norme, normr, normutilde:
...
%Tracking Error Plot
% norme=zeros(1,length(t));
% for jj=1:length(t)
% norme(jj)= norm([ex(jj) ephi(jj)]);
% end
norme = sqrt(ex.^2+ephi.^2);
%Filtered Tracking Error Plot
% normr=zeros(1,length(t));
% for kk=1:length(t)
% normr(kk)= norm([rx(kk) rphi(kk)]);
% end
normr= sqrt(rx.^2+rphi.^2);
%Estimation Errors
% normutilde=zeros(1,length(t));
% for aa=1:length(t)
% normutilde(aa)= norm([utildex(aa) utildephi(aa)]);
% end
normutilde= sqrt(utildex.^2+utildephi.^2);
...
2 Comments
Sulaymon Eshkabilov
on 6 Nov 2021
This way (vectorization) speeds up the calculation significantly if there is a large number of simulations (iterations). In your exercise, the size of t is large that makes the whole simulation for a number of times slower in your [for end] loop based calcs than this proposed vectorization based calcs.
More Answers (1)
Walter Roberson
on 6 Nov 2021
[t, xout] = 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), tspan, x0, options);
num_t = length(t);
for jj = 1 : num_t
x = xout(jj,1);
xdot = xout(jj,2);
xdotdot = xout(jj,3);
phidot = xout(jj,5);
ux = xout(jj,6);
uphi = xout(jj,7);
thetahat = xout(jj,8:13);
ahatx = xout(jj,14);
ahatphi = xout(jj,14);
When you got to xddotdot, use t(K) .
Eventually you will get to
%Tracking Error Plot
and the for jj loop above should likely end just before that.
Then in the lines below that, use num_t instead of length(t)
norme=zeros(1,num_t);
for jj=1:num_t
norme(jj)= norm([ex(jj) ephi(jj)]);
end
and so on.
You will promptly have the problem that your ex and ephi are scalars. You rejected the suggestion to vectorize everything, so you are dealing with scalars; each iteration of the first for jj loop, you are creating scalars only, and not recording any of the values, so the vector of values will not exist when you get to the plotting.
So you now have three choices:
- Go back and index all over the place inside the loop; OR
- Do not index all over the place inside the loop, but at the end of the loop, take copies of the scalars you want to hold on to, storing them in indexed variables; OR
- Mostly vectorize instead of looping
Example of vectorization:
%the below assumes that phi and so on are colum vectors
Z = zeros(size(phi,1));
yx=[xdotdot, (sin(phi)).^(2).*xdotdot-g.*cos(phi).*sin(phi), phidot.^(2).*sin(phi), Z, xdot, Z];
This would be a num_t x 6 array.
Then a few lines later you currently have
udx=Yx*x(8:13)+ex+Bx*rx;
is that x(8:13) intended to be the same as thetahat, which is num_t x 6 ? When you wrote the code, were you picturing doing element-by-element multiplication of a 1 x 6 and a 1 x 6, getting out a 1 x 6? Or were you picturing doing algebraic matrix multiplication between a 1 x 6 and a 6 x 1, getting out a 1 x 1 scalar ? If you were thinking in terms of algebraic matrix multiplication, then you should recode.
Unfortunately, there is noise at my place right now and I cannot concentrate enough to come up with vectorized replacment code for the case where you were thinking in terms of inner product... the ones that are coming to mind would get you a num_t by num_t output instead of a num_t by 1 output.
2 Comments
Walter Roberson
on 6 Nov 2021
Let's see:
A = reshape(1:20, 4, 5)
B = reshape(50:-1:31, 4, 5)
C1 = diag(A * B')
C2 = sum(A .* B, 2)
So
udx=Yx*x(8:13)+ex+Bx*rx;
can be coded as
udx = sum(Yx .* thetahat, 2) + some stuff
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!