Produce matrices through for loop
2 次查看(过去 30 天)
显示 更早的评论
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
1 个评论
Stephen23
2019-5-28
编辑:Stephen23
2019-5-28
" 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:
采纳的回答
Geoff Hayes
2019-5-27
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
21 个评论
Ali Esmaeilpour
2019-5-27
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
2019-5-27
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
2019-5-28
编辑:Ali Esmaeilpour
2019-5-28
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.
![1.jpg](https://www.mathworks.com/matlabcentral/answers/uploaded_files/221413/1.jpeg)
Geoff Hayes
2019-5-28
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
2019-5-28
编辑:Geoff Hayes
2019-5-29
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
2019-5-29
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
2019-5-29
编辑:Ali Esmaeilpour
2019-5-29
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
2019-5-29
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
2019-5-29
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.
![1.jpg](https://www.mathworks.com/matlabcentral/answers/uploaded_files/221914/1.jpeg)
Geoff Hayes
2019-5-30
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
2019-5-30
编辑:Ali Esmaeilpour
2019-5-30
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
2019-5-30
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
2019-5-30
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
2019-5-30
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
2019-5-30
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
2019-5-30
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?
Geoff Hayes
2019-5-31
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
2019-5-31
yeah tnx my friend i debugged that main code and solve it with sedumi. tnx again for you consideration. cheers!
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- 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)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)