Issues with Horzcat for defining a variable as a row vector

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

[Matthew Tortorella]

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?

[Walter Roberson]

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.

[Matthew Tortorella]

@Walter Roberson in that segment of code I'm trying to redefine what my variables are in terms of all my states x. For instance, in my function handle I described x(1)=x and x(2)=xdot, and so on. I see now that the x output from my ODE is a huge matrix with 15 columns and a bunch of rows, corresponding to the values of each of my 15 states for all the iterations of time. What I would like to do is define the names of each of my states with a variable (for example I would call x(2)=xdot as I did in my function handle) for each iteration in time and then calculate utilde using those defined variables, and then calculate the variables value for the next iteration in time and then recalculate utilde, and then I would like to plot utilde's values for all the ierations of time (hopefully seeing convergence to 0). How could I adjust my variable definition so it is equal to the corresponding state at whatever instance in time?

Answer 1

[Sulaymon Eshkabilov]

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);
...

[Matthew Tortorella]

Is there any advantage to doing it this way rather than the way I did it with the for loop @Sulaymon Eshkabilov?

[Sulaymon Eshkabilov]

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.

Answer 2

[Walter Roberson]

[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:

1. Go back and index all over the place inside the loop; OR
2. 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
3. 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.

[Matthew Tortorella]

Yes i was hoping of doing the inner product and getting a 1x1 scalar. I/m going to try and come up with a vectorized solution I guess because the scalars for each iteration is giving me problems.

[Walter Roberson]

Let's see:

A = reshape(1:20, 4, 5)

A = 4×5

1 5 9 13 17 2 6 10 14 18 3 7 11 15 19 4 8 12 16 20

B = reshape(50:-1:31, 4, 5)

B = 4×5

50 46 42 38 34 49 45 41 37 33 48 44 40 36 32 47 43 39 35 31

C1 = diag(A * B')

C1 = 4×1

1730 1890 2040 2180

C2 = sum(A .* B, 2)

C2 = 4×1

1730 1890 2040 2180

So

udx=Yx*x(8:13)+ex+Bx*rx;

can be coded as

udx = sum(Yx .* thetahat, 2) + some stuff