Produce matrices through for loop

Hello people! I want to produce N number of matrices while k=1:N and I want the result in a way like h(k). I mean h(1),h(2) etc which are matrices that depend on k. I mean creating them through a for loop. I tried the following code but it failed to produce N matrices and it just gave one matrix:

clc;
clear;
close all;

hx = [-1/sqrt(5);-2/sqrt(5)];
hu = [0;0];
N = 25;
Ek1 = zeros(N+1,1);
Ek2 = ones(N+1,1);
Ek3 = ones(N,1);
Ek4 = zeros(N,1);
h = zeros(102,1);
for k=1:N
    if k==1
        h(:,:) = [kron(Ek2,hx);kron(Ek3,hu)];
    else
        h(:,:) = [kron(Ek1,hx);kron(Ek4,hu)];
    end
end

[Stephen23]

" I want to produce N number of matrices..."

Then the best** solution is to use indexing, exactly as your question already shows

** best in the sense simplest, neatest, easiest to write, easiest to debug, and most efficient.

Note that dynamically access variable names is one way that beginners force themselves into writing slow, complex, buggy code that is hard to debug. Read this to know why:

https://www.mathworks.com/matlabcentral/answers/304528-tutorial-why-variables-should-not-be-named-dynamically-eval

Answer 1

[Geoff Hayes]

Ali - if the output matrix of each iteration is of dimension 102x1, then you could store each output as a column in a 102xN matrix. For example,

    h = zeros(102,N);
    for k=1:N
        if k==1
            h(:,k) = [kron(Ek2,hx);kron(Ek3,hu)];
        else
            h(:,k) = [kron(Ek1,hx);kron(Ek4,hu)];
        end
    end

[Ali Esmaeilpour]

Tnx but i dont think it's right. I want to use this code in another code in a for loop and assume N=5 and k=1:N I want to produce five matrices named h(k) to Place them in the optimization FORMULA I mean when k is 1 it Places h(1) in FORMULA etc

[Geoff Hayes]

Well instead of h(1) you would use h(:,1).

I want to produce five matrices named h(k) t

It almost sounds as if you want five distinct matrices. That sort of programming (dynamically creating matrices/arrays) is usually discouraged. In this case, a matrix can be used in place of this.

[Ali Esmaeilpour]

actually I wanted to have h(k) that is in the picture, but the way you edited my code doesn't give h(k) correctly I mean in a modular way that it replaces k in h(k) itself.

[Geoff Hayes]

Ali - I don't understand why having h(:,k) doesn't give you what you want. Is it the colon : that is causing the problem? Would you rather just use a cell array and index as h{k}?

[Ali Esmaeilpour]

when I use your h(:,k) code separatly the result is h 22*5 and when I copy it to my main program and run h is 22*2 I dont know why and I want it 22*1 I mean for k==1 I need 22*1 matrix I wanna be sure that it replaces a 22*1 matrix when k==1 itself but I think it doesn't do that. I give you my full optimization program maybe you can get better picture of my problem. the following program has dimension problem because of h, but I just wanna show you h 22*2 in my main code and the usage of my question in my main code:

clc;
clear;
close all;
%% Time structure
dt=0.01;          %sampling time

%% System structure
A = [1.02  -0.1
     0.1   0.98];
 
B = [0.5 0
    0.05 0.5];

G = [0.3 0
     0 0.3];
 
Qf = [50 0
      0 50];

[A, B] = c2d(A,B,dt);

[nx,nu] = size(B);

%% MPC parameters
Q  = eye(2);

R  = eye(2)*50;

N = 5;

%% Building block matrices
I = eye(2,2);

q = 2;

m = 2;

Qbar = blkdiag(Q,Q,Q,Q,Qf,R,R,R,R);

Fq = blkdiag(sqrt(Q),sqrt(Q),sqrt(Q),sqrt(Q),sqrt(Q),sqrt(Qf),sqrt(R),sqrt(R),sqrt(R),sqrt(R),sqrt(R));

Gxx = zeros(2*N , size(A,1));
for i = 1:N
    Gxx(q*i-q+1:q*i , :) = A^i;
end

Gxu = zeros(2*N , 2*N);
for i = 1:N
   for j = 1:i
      Gxu(q*i-q+1:q*i , m*j-m+1:m*j) = A^(i-j)*(B); 
   end    
end

Gxw = zeros(2*N , 2*N);
for i = 1:N
    for j = 1:i
    Gxw(q*i-q+1:q*i , m*j-m+1:m*j) = A^(i-j)*(I);    
    end
end

%% Solve using YALMIP

K = sdpvar(repmat(10,1,N),repmat(22,1,N));

alpha = 0.975;

r = 1.96;

Sigmaw = eye(10,10);

Sigma0 = eye(2,2);

xbar0 = ones(2,20);

Cs = [Gxx*xbar0;zeros(12,20)];

hx = [-1/sqrt(5);-2/sqrt(5)];

hu = [0;0];

Ahat = [Gxx*xbar0 eye(10,2);zeros(12,20) zeros(12,2)];

Bhat = [Gxu;eye(12,10)];

Abar = [eye(11,10);zeros(11,10)];

M = Gxx*Sigma0*Gxx' + Gxw*Sigmaw*Gxw';

Mhat = [M zeros(10,12);zeros(12,10) ones(12,12)];

Euhat = ones(22,20);

Ek = [zeros(12,10);eye(10,10)];

g = 3;

Ek1 = zeros(N+1,1);
Ek2 = ones(N+1,1);
Ek3 = ones(N,1);
Ek4 = zeros(N,1);

constraint=[];
 
objective=0;

for k = 1:N   

    if k==1
        h(:,k) = [kron(Ek2,hx);kron(Ek3,hu)];    
    else
        h(:,k) = [kron(Ek1,hx);kron(Ek4,hu)];    
    end
    
    b = sqrt(h'*((Abar+Bhat*K{k}*Ek)*M*(Abar+Bhat*K{k}*Ek)')*h);
    
    objective = objective + 0.5*trace(Fq'*(Ahat+Bhat*K{k})*Mhat*(Ahat+Bhat*K{k})'*Fq);
    
    constraint = [constraint, h'*(Cs+Bhat*K{k}*Euhat)+r*(sqrt(h'*(Abar+Bhat*K{k}*Ek)*M*(Abar+Bhat*K{k}*Ek)'*h))<=g, [(h'*(Abar+Bhat*K{k}*Ek)*M*(Abar+Bhat*K{k}*Ek)'*h) (h'*(Abar+Bhat*K{k}*Ek));((Abar+Bhat*K{k}*Ek)'*h) (Gxx*Sigma0*Gxx' + Gxw*Sigmaw*Gxw')]>=0];

end

options = sdpsettings('solver','sedumi');
optimize (constraint,objective,options);

K = value(K);

[Geoff Hayes]

In these lines of code

b = sqrt(h'*((Abar+Bhat*K{k}*Ek)*M*(Abar+Bhat*K{k}*Ek)')*h); 
objective = objective + 0.5*trace(Fq'*(Ahat+Bhat*K{k})*Mhat*(Ahat+Bhat*K{k})'*Fq);
constraint = [constraint, h'*(Cs+Bhat*K{k}*Euhat)+r*(sqrt(h'*(Abar+Bhat*K{k}*Ek)*M*(Abar+Bhat*K{k}*Ek)'*h))<=g, [(h'*(Abar+Bhat*K{k}*Ek)*M*(Abar+Bhat*K{k}*Ek)'*h) (h'*(Abar+Bhat*K{k}*Ek));((Abar+Bhat*K{k}*Ek)'*h) (Gxx*Sigma0*Gxx' + Gxw*Sigmaw*Gxw')]>=0];

you reference h. Should this be h(:,k) instead (since this code is in the for loop)?

[Ali Esmaeilpour]

I don't know putting this h(:,k) there instead of h could be right because we define an "if else" statement there and h has to be the result of that "if else". when you put h(;,k) there instead of h and run the program you see that h(:,k) is zero and "if else" doesn't work at all.

h(:,k)

ans =

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

[Geoff Hayes]

if all zeros (for k > 1) is the problem then with this line of code

h(:,k) = [kron(Ek1,hx);kron(Ek4,hu)];    

Note that Ek1 and Ek4 are defined to be all zeros...are you sure that this line of code is doing what you expect?

[Ali Esmaeilpour]

I sent the picture of h(k). I wrote this code to have the matrix h(k) which is in the picture. note that i=1 and hx and hu are known and k=1:N

I don't know if the code is wrong. I've asked you guys for that.

[Geoff Hayes]

Maybe the problem is with your definition of [ekn] vector (not sure how best to write that). My thinking is that for any k then

e1n = [1; 0; 0; 0; .... 0];
e2n = [0; 1; 0; 0; .... 0];
e3n = [0; 0; 1; 0; .... 0];
...
enn = [0; 0; 0; 0; .... n];]

So on the kth iteration of the for loop, these arrays are all zeros except in the kth position. Your code may look more like

for k=1:N
    Ekn1 = zeros(N,1);
    Ekn1(k) = 1;
    Ekn2 = zeros(N+1,1);
    Ekn2(k) = 1;
    h(:,k) = [kron(Ekn2,hx);kron(Ekn1,hu)];
    % etc.
end

The above may or may not be what you want...it is my interpretation of your above image.

[Ali Esmaeilpour]

eh... excuse me! how come e1n is not ones(N,1) and it is [1;0;0;0;...;n] ??? and other e2n e3n etc are not zeros(N,1)? because i=1 and when k=i, [ekn]=1 else it's zero.

it says if k=i then [ekn]=1 else it is zero and i=1 is constant. another subject is that the picture says ekn is of dimension N*1 so N rows and one column I suppose!

have you consider these? because I can't see any "if else" in your interpretation...

[Geoff Hayes]

where does it say that e1n is all ones and all of the others are all zeros? what does the i in [ekn]i represent if not the ith element of the vector?

[Ali Esmaeilpour]

I think you are right, but if this is the way to construct "h" I should use h(:,k) in my main code or "h"? I mean in those parts h is used like constraint=...

I mean this part:

for k = 1:N   

Ekn1 = zeros(N,1);
Ekn1(k) = 1;
Ekn2 = zeros(N+1,1);
Ekn2(k) = 1;
h(:,k) = [kron(Ekn2,hx);kron(Ekn1,hu)];

b = sqrt(h(:,k)'*((Abar+Bhat*K{k}*Ek)*M*(Abar+Bhat*K{k}*Ek)')*h(:,k));

objective = objective + 0.5*trace(Fq'*(Ahat+Bhat*K{k})*Mhat*(Ahat+Bhat*K{k})'*Fq);

constraint = [constraint, ((h(:,k)'*(Cs+Bhat*K{k}*Euhat))+(r*(sqrt(h(:,k)'*(Abar+Bhat*K{k}*Ek)*M*(Abar+Bhat*K{k}*Ek)'*h(:,k)))))<=g, [(h(:,k)'*(Abar+Bhat*K{k}*Ek)*M*(Abar+Bhat*K{k}*Ek)'*h(:,k)) (h(:,k)'*(Abar+Bhat*K{k}*Ek));((Abar+Bhat*K{k}*Ek)'*h(:,k)) (Gxx*Sigma0*Gxx' + Gxw*Sigmaw*Gxw')]>=0];

end

[Geoff Hayes]

Well this "main" code is inside the for loop so I think that you want to use h(:,k)...but without seeing why you are setting up the constraint in this way (which equation are you replicating?), I can only guess.

[Ali Esmaeilpour]

It's Yalmip code not standard matlab code I mean it's the form of a problem you give to Yalmip it solves it. I used h(:,k) and it says no solver applicable but with h it solves it right or wrong it's not obvious yet.

another weird subject is that when I preallocate h and run it takes forever and no result and without preallocating h's dimension becomes 22*2

I don't know I can't figure out what is this progrram do lol...

[Geoff Hayes]

Not sure why it would run slower when pre-allocating array unless it isn't be pre-allocated correctly...unless your N has changed. I think I was using 102 because your N was 25...maybe that is no longer the case?

[Ali Esmaeilpour]

I pre-allocate with h = zeros((2*N + 1) * size(hx,1), N);

[Geoff Hayes]

and what is N?

[Ali Esmaeilpour]

N can be 5 or 10 or 15 or 20 or 25

[Geoff Hayes]

For the case where N is 5, then

 h = zeros((2*N + 1) * size(hx,1), N);

h is a 22x5 array. Is this correct?

[Ali Esmaeilpour]

yeah tnx my friend i debugged that main code and solve it with sedumi. tnx again for you consideration. cheers!