- You will see updates in your activity feed.
- You may receive emails, depending on your notification preferences.

- 17 questions asked
- 0 answers
- 0 accepted answers
- reputation: 0

Asked by Ali Esmaeilpour
on 27 May 2019

- 17 questions asked
- 0 answers
- 0 accepted answers
- reputation: 0

- 17 questions asked
- 0 answers
- 0 accepted answers
- reputation: 0

Accepted Answer by Geoff Hayes

- 0 questions asked
- 2,740 answers
- 1,468 accepted answers
- reputation: 8,154

2 views (last 30 days)

2 views (last 30 days)

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

- 0 questions asked
- 2,740 answers
- 1,468 accepted answers
- reputation: 8,154

Answer by Geoff Hayes
on 27 May 2019

- 0 questions asked
- 2 740 answers
- 1 468 accepted answers
- reputation: 8 154

Accepted Answer

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

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.

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

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)?

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

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?

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.

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.

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

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

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

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

and what is N?

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

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?

Sign in to comment.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply TodayUnable to complete the action because of changes made to the page. Reload the page to see its updated state.

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

Select web siteYou can also select a web site from the following list:

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

- América Latina (Español)
- Canada (English)
- United States (English)

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

## 1 Comment

## Stephen Cobeldick (view profile)

## Direct link to this comment

https://fr.mathworks.com/matlabcentral/answers/464248-produce-matrices-through-for-loop#comment_709326

Sign in to comment.