Need help with discrete time SS controllability & phase portrait

I'm very basic in MATLAB, I want to solve the above, I've read some documentation, theory and other Q/A's, but am unable to solve this. I'm attaching the code below till the point I understand and after that it gives an error.

clc;
clear;
A = [1 0.4;0 1];
B = [0.08;0.4];
% I've assumed y=x1
C = [1 0]; 
x0 = [1;1];
x50 = [-5;12];

% I tried using idss(A,B,C,0,x0,1) but gives some error, if I use idss(A,B,C,0,x0,x0,x50,1) it doesn't give error, but introduces some K & idk how to proceed then. I have taken time-step as 1, idk if I should alter that.
% How do I specify both initial (i.e. x(k=0)) & final intended position (i.e. x(k=50))? 
sysss = ss(A,B,C,0,1)

% not sure if this is correct since 't or i' only goes to 50 not infinity.
% Also if I should use lyap or gram?

Wd = lyap (A, B*B') 
i=0;
n=50;
u(i) = -(A^(n-1-i)*B)'*inv(Wd)*(A^(n)*x0-x50);
%            OR        
u(i) = -(B'*(A')^(n-1-i))*inv(Wd)*(A^(n)*x0-x50);

% both give error, and idk how else to solve/proceed
% Once I get the values of u(0) to u(49)-in matrix form I think according to the theory, then how to plot them, again idk
% the syntax for phase portrait (even if I find it, I cant solve the errors)
THANKS!

Answer 1

[Paul]

Do you have a source for the equations you're trying to implement to solve for u?

I thought that the controllability Grammian as the solution of the discrecte Lyapunov equation is the sum over n = 0 to inf, as discussed in

doc gram

But in this problem we only need the "portion" of the Grammian from n = 0 to 50. So the first thing is to find the corrrect equation that solves for u. Then we can figure out the Matlab code. Is this link helpful: General State Transfer and Reachabilty ?

Also, the lyap() function is for continuous systems. The dlyap() function is used for discrete systems. However, neither of these functions need be used for this problem, I don't believe.

[Woohoo PAKFA]

Thanks Paul for the quick reply.

I've attached a slide which I used for getting the u(i) equation.

From what I've been taugh/understood is that the Controllability gramian is defined for 0->n-1 but if we extend it to infinity (assuming it's stable I think, I have to look at the diff b/w BIBO stability/stable i.s.L/asymptotically stable), only then we can use the Lyapunov eqn instead for calculating the gramian.

[Paul]

OK. I didn't realize that Wd defined for finite k is also called the controllability Grammian. And you are correct that if k doesn't go to infinity then you can't use dlyap() or gram() (actually gram() has options that might be applicable here, I didn't check).

Step one is to define the system. The A and B matrices are given. We'll set the C matrix to an identity so that later on we can get the state vector as the output of the system as each state:

A = [1 0.4;0 1];
B = [0.08;0.4];
C = eye(2);
sys = ss(A,B,C,0,-1); % defines sys as a discrete time system with undefined sample time

Step two is to define Wd. The simplest way is to use a loop with Wd initialized to zero

nfinal = 50;
Wd = zeros(2)
for kk = 0:nfinal
% what goes here?
end

Can you fill in the loop?

Step three would then be to use a loop to fill in u(i) for i = 0 - 49.

u = zeros(nfinal+1,1); % actually, add an extra element to u, which will be needed later
for ii = 0:(nfinal-1)
% u(ii+1) =    % Matlab uses 1-based indexing 
end

Can you fill in that loop?

Now you have sys and the vector of inputs u. To verify the solution you can use

doc lsim

and then plot the result.

[Woohoo PAKFA]

I tried, but I'm not too familiar with the syntax, so am getting errors.

nfinal = 50;
Wd = zeros(2)
for kk = 0:nfinal   % Any reason to use 'kk' & not single 'k' ; also nfinal should be 49 (since 0->n-1), right? 
   Wd (kk+1) = ((A^kk)*B*B'*(A')^kk)        % what goes here?
end
Wdc = sum(Wd)


u = zeros(nfinal+1,1); % actually, add an extra element to u, which will be needed later
for ii = 0:(nfinal-1)
% u(ii+1) =    % Matlab uses 1-based indexing 
end
%  similarly I have a hard time understanding the pecularities of the
%  syntax and the syntax itself, I searched other Q/A too, some had done
%  summation without 'for' loop but used '.' before *ing or /ing, so I'm
%  confused

[Paul]

Let's keep the entirety of the code:

A = [1 0.4;0 1];
B = [0.08;0.4];
C = eye(2);
sys = ss(A,B,C,0,-1); % defines sys as a discrete time system with undefined sample time

% define the boundary conditions
x0 = [1 ; 1];
x50 = [-5 ; 12];
nfinal = 50;  % nfinal = 50 because the final state is x50, i.e, n = 50 in the problem statement

% now we have two choices:
% 1. We can initialize Wd as a 2 x 2 matrix and build up the sum incrementally, or
% 2. We can initialize Wd as 2 x 2 x (nfinal + 1) matrix, compute each 2 x 2 pane, and then sum along the third dimension.
% Let's choose the first method
Wd = zeros(2)  % Wd is initialized to a 2x2 matrix

Wd = 2×2

0 0 0 0

for kk = 0:(nfinal-1) % I use double letter variables for counters, just a personal convention
%Wd = Wd + ?  % this should be the summation in the slide
end

Can you fill in the loop?

[Woohoo PAKFA]

A = [1 0.4;0 1];
B = [0.08;0.4];
C = eye(2);
sys = ss(A,B,C,0,-1); % defines sys as a discrete time system with undefined sample time

% define the boundary conditions
x0 = [1 ; 1];
x50 = [-5 ; 12];
nfinal = 50;  % nfinal = 50 because the final state is x50, i.e, n = 50 in the problem statement

% now we have two choices:
% 1. We can initialize Wd as a 2 x 2 matrix and build up the sum incrementally, or
% 2. We can initialize Wd as 2 x 2 x (nfinal + 1) matrix, compute each 2 x 2 pane, and then sum along the third dimension.
% Let's choose the first method
Wd = zeros(2)  % Wd is initialized to a 2x2 matrix
for kk = 0:(nfinal-1) % I use double letter variables for counters, just a personal convention
%Wd = Wd + ?  % this should be the summation in the slide
Wd = Wd + ((A^kk)*B*B'*(A')^kk);
end
Wdc = Wd
inv(Wdc)

u = zeros(nfinal+1,1); % actually, add an extra element to u, which will be needed later
for ii = 0:(nfinal-1)
% u(ii+1) =    % Matlab uses 1-based indexing 
u(ii+1) = -(A^(nfinal-1-i)*B)'*inv(Wdc)*(A^(nfinal)*x0-x50)
end
% after following the synax for "Extract Simulated Response Data" in the doc lsim i'm getting another error.
uM = [u u]
[y,ii+1] = lsim(sys,uM,ii+1)
size(y)
[y,ii+1,x] = lsim(sys,uM,ii+1)
size(x)
plot(t,x)

I got Wd and then u(ii+1) with all values as the same. Now to plot it, it spits out another error

[Paul]

There was an error in the calculation of the u(ii+1). The rhs of the assignment was using i, when it should have been using the counter variable ii. Unless otherwised defined, i = sqrt(-1).

I'm not sure what you saw in the doc page for lsim. Anyway, the correct usage for this case is shown below. Note that it's important to specify the initial condition. Compare that line of code to what's in the doc page. Note that we set up the system so the output y is the state vector x, so we only need the single output form of the call to lsim.

Also, it's almost never needed to use inv(). Better to use the / or \ operators as shown below. Check the doc pages for inv, mrdivide, and mldivide for more discussion.

A = [1 0.4;0 1];
B = [0.08;0.4];
C = eye(2);
sys = ss(A,B,C,0,-1); % defines sys as a discrete time system with undefined sample time
% define the boundary conditions
x0 = [1 ; 1];
x50 = [-5 ; 12];
nfinal = 50;  % nfinal = 50 because the final state is x50, i.e, n = 50 in the problem statement
% now we have two choices:
% 1. We can initialize Wd as a 2 x 2 matrix and build up the sum incrementally, or
% 2. We can initialize Wd as 2 x 2 x (nfinal + 1) matrix, compute each 2 x 2 pane, and then sum along the third dimension.
% Let's choose the first method
Wd = zeros(2)  % Wd is initialized to a 2x2 matrix

Wd = 2×2

0 0 0 0

for kk = 0:(nfinal-1) % I use double letter variables for counters, just a personal convention
%Wd = Wd + ?  % this should be the summation in the slide
Wd = Wd + ((A^kk)*B*B'*(A')^kk);
end
Wdc = Wd

Wdc = 2×2

1.0e+03 * 1.0666 0.0800 0.0800 0.0080

inv(Wdc)

ans = 2×2

0.0038 -0.0375 -0.0375 0.5002

u = zeros(nfinal+1,1); % actually, add an extra element to u, which will be needed later
for ii = 0:(nfinal-1)
% u(ii+1) =    % Matlab uses 1-based indexing 
%u(ii+1) = -(A^(nfinal-1-i)*B)'*inv(Wdc)*(A^(nfinal)*x0-x50);
%u(ii+1) = -(A^(nfinal-1-ii)*B)'*inv(Wdc)*(A^(nfinal)*x0-x50);  % replace i with the counter ii!
 u(ii+1) = -(A^(nfinal-1-ii)*B)'/(Wdc)*(A^(nfinal)*x0-x50);     % use / instead of inv
end
% after following the synax for "Extract Simulated Response Data" in the doc lsim i'm getting another error.
%uM = [u u]
%[y,ii+1] = lsim(sys,uM,ii+1)
% size(y)
% [y,ii+1,x] = lsim(sys,uM,ii+1)
% size(x)
% plot(t,x)
y = lsim(sys,u,[],x0);
plot(y(:,1),y(:,2),y(1,1),y(1,2),'x',y(end,1),y(end,2),'o')
xlabel('x1')
ylabel('x2')
legend('trajectory','x0','xfinal')

plot(0:49,u(1:50)),xlabel('n'),ylabel('u')

[Woohoo PAKFA]

Hey Paul,

Finally got some output without errors thanks to your time and help, I'm really grateful for your step-by-step instructions!!

I have some doubts regarding the syntax though, if you wouldn't mind clarifying those:

1. I missed the 'ii' in the rhs, so no wonder I was getting all the same values, my bad.

   2.for kk = 0:(nfinal-1) % I use double letter variables for counters, just a personal convention
     Wd = Wd + ((A^kk)*B*B'*(A')^kk); 
     end
% When writing this for loop, does MATLAB consider different values of 'Wd' for each subsequent summation,
% as in when kk=1 . . ., the 1st Wd (lhs) is just assigning a value to Wd, which
% is equal to the 2nd Wd (rhs-this it recalls from the previous calculation?) + second term. . . ., right?

  3.u = zeros(nfinal+1,1); % actually, add an extra element to u, which will be needed later
   % why the extra element?
    for ii = 0:(nfinal-1)
   % u(ii+1) =    % Matlab uses 1-based indexing-> Could you explain this a bit more?
     u(ii+1) = -(A^(nfinal-1-ii)*B)'/(Wdc)*(A^(nfinal)*x0-x50);     % use / instead of inv
    end
 % Could you explain this formula and how it works, I understood that u is
 % initiated as a size-51x1 matrix, after that i'm a bit confused.

4. y = lsim(sys,u,[],x0); 
   % so this equation basically generates the response data (x1,x2) in 'y' using the time solution equation & [] means that the time interval isn't specified?
   plot(y(:,1),y(:,2),y(1,1),y(1,2),'x',y(end,1),y(end,2),'o') 
   % have some trouble understanding the y(1,1),y(1,2) & the y(end,1),y(end,2) parts, we just have to plot x1 vs x2, so what are these terms?
   xlabel('x1')
   ylabel('x2')
   legend('trajectory','x0','xfinal')
   plot(0:49,u(1:50)),xlabel('n'),ylabel('u')

Sorry for the delayed response, I was also doing another questions, it was simpler than this and I did it myself after a few errors and reading other Q/A's online.

Again a HUGE thanks PAUL for your help!! I've realized that I can't run away from coding no matter how much I hated it (Mech background), and since I'm also learning ML & NN now, I think it's better to learn it.

[Paul]

2. Yes, that line basically says that each time through the loop assign as the new value of Wd the old value of Wd + that other stuff. It is implementing the summation in the middle of the slide that defines Wd.

3. Matlab arrays start with an index of 1. So we can't use the u(0) for the variable in the problem of u0. Instead, u(1) in the Matlab array is u0 in the problem. I thought that you basically wrote that equation for u(ii+1). All it is doing is filling in the values of u0, u1, u2, ... u49 into the Matlab array u using the equation in the middle of the slide, keeping in mind that u0 maps to u(1), u1 to u(2), etc. From the problem statement, you need to find u0 - 49, which are fifty elements of u. But x50 is the is the 51st element of the state vector (counting from 0), so we need to make sure that the Matlab array u has 51 elements so that lsim will compute 51 iterations of the state vector, x0 - x50. The 51st element of the Matlab array u isn't actually used for any computation, which is why it can stay zero.

4. The lsim comand using [] as the time input just means that "time vector" elements are just the discrete integers 0,1, ...50. Each row of y is an iteration of the state vector. y(1,1) and y(1,2) correspond to the first and second elements of the initial state vector (x0), and y(end,1) and y(end,2) are the elements of the final (or "end" ) state vector (x50). I added those single points on the plot with the markers (the x and the o) to show the initial and ending position on the state trajectory. Just adds a bit of clarity, IMO.