Guys, here's my question : Yalmip simulation in research project of MIMO beam forming optimization

[yangke Sun]

Hey guys, 

I am recently doing my research project of MIMO stuff. 

Idea: 

To obtain an hermitian complex matrix by solving a convex matrix problem. Unfortunately, i always get "Warning: Solver not found (sedum)". 

The code is draft， sorry about this 

############################## the main function of Yalmip, need to cope the seconde function named test1, where I defined some parameters.######

clear;clc;

warning('off','YALMIP:strict') 

n=64;

W_tire = sdpvar(n,n,'hermitian','complex');

R1= sdpvar(1,1);

R2= sdpvar(1,1);

Re= sdpvar(1,1);

Lambda1 = sdpvar(1,1);

Lambda2 = sdpvar(1,1);

Lambda3 = sdpvar(1,1);

Lambda4 = sdpvar(1,1);

Lambda5 = sdpvar(1,1);

Lambda6 = sdpvar(1,1);

Lambda7 = sdpvar(1,1);

Lambda8 = sdpvar(1,1);

Lambda9 = sdpvar(1,1);

[he,ge,h_wavy,h_chapeau,g_wavy,g_chapeau,e_wavy,e_chapeau,number,power1,power2,sigma] = test1();

g = g_chapeau+g_wavy;

h = h_chapeau + h_wavy;

e = e_chapeau + e_wavy;

sigma1=1;

sigma2=1;

P_vector = zeros(1,number^2);

P = zeros(number^2);

v=1;

for z = 1:number   

    

    for x = 1:number

        

        P_vector(z+number*(x-1)) = 1; 

        

        P(v,:) = P_vector;

        

        P_vector = zeros(1,number^2);

        v = v+1;

    

    end

    

end

% test the kronecker function

% functions with h 

h_tire = kron(h_chapeau, ones(number,1));

h_inversechapeau = kron(ones(number,1),h_chapeau);

delt_h_tire = kron(h_wavy,ones(number,1)) ;

delt_h_inversechapeau = kron(ones(number,1),h_wavy);

% functions with g

g_tire = kron(g_chapeau, ones(number,1));

g_inversechapeau = kron(ones(number,1),g_chapeau);

delt_g_tire = kron(g_wavy,ones(number,1)) ;

delt_g_inversechapeau = kron(ones(number,1),g_wavy);

% functions with e

e_tire = kron(e_chapeau, ones(number,1));

e_inversechapeau = kron(ones(number,1),e_chapeau);

delt_e_tire = kron(e_wavy,ones(number,1)) ;

delt_e_inversechapeau = kron(ones(number,1),e_wavy);

% construct expressions appear in page 7

D_R = kron(eye(number), ones(number,1));

D_L = kron(ones(number,1),eye(number));

G1 = [eye(number) zeros(number) zeros(number)];

G2 = [zeros(number) eye(number) zeros(number)];

G3 = [zeros(number) zeros(number) eye(number)];

c = [h,g,e];

c = c(:);

delt_c = [h_wavy g_wavy e_wavy];

delt_c = delt_c(:);

% test the equation of D_R*G1*delt_c = delt_h_tire

try1 = D_R*G1*delt_c-delt_h_tire;

try2 = D_L*G2*delt_c - delt_g_inversechapeau;

% test the lemma 1 

p = sigma1*randn(number,1);

q = sigma2*randn(number,1);

A = sigma1*sigma2*randn(number);

lemma1_left = transpose(p)*A*q;

lemma1_right = transpose((kron(q,ones(number,1))).*(kron(ones(number,1),p)))*A(:);

S = zeros(number^2,number^2);

SW = zeros(number);

% this matrix is for the use of the lemma 3

SPW = zeros(number);

for m = 1:number

    S(((m-1)*number+1):m*number,:) = [zeros(number,(m-1)*number) eye(number) zeros(number, (number^2-m*number)) ];

    

    SW = SW + S(((m-1)*number+1):m*number,:)*W_tire*(S(((m-1)*number+1):m*number,:))';

   

    SPW = SPW + S(((m-1)*number+1):m*number,:)*P*W_tire*P'*(S(((m-1)*number+1):m*number,:))';

end

lemma2_right = transpose(delt_c)*(transpose(G1)*SW*G1)*conj(delt_c)+...,

2*real(transpose(h_chapeau)*SW*G1*conj(delt_c))+ transpose(h_chapeau)*SW*conj(h_chapeau);

Q1 = power2*transpose(G2)*transpose(D_R)*diag(h_inversechapeau)*W_tire*...,

diag(conj(g_tire))*D_L*G1+power2*transpose(G1)*transpose(D_L)*...,

diag(g_tire)*W_tire*diag(conj(h_inversechapeau))*D_R*G2+...,

power2*transpose(G2)*transpose(D_R)*diag(h_inversechapeau)*W_tire*...,

diag(conj(h_inversechapeau))*D_R*G2+power2*transpose(G1)*...,

transpose(D_L)*diag(g_tire)*W_tire*diag(conj(g_tire))*D_L*G1;

q1 = power2*transpose(g_tire.*h_inversechapeau)*W_tire*...,

diag(conj(h_inversechapeau))*D_R*G2+power2*transpose(g_tire.*...,

h_inversechapeau)*W_tire*diag(conj(g_tire))*D_L*G1;

k1 = power2*transpose(g_tire.*h_inversechapeau)*W_tire*conj(g_tire.*h_inversechapeau);

a1 = transpose(delt_c)*Q1*conj(delt_c)+2*real(q1*conj(delt_c))+k1;

Q2=power1*transpose(G1)*transpose(D_R)*diag(h_inversechapeau)*W_tire*...,

diag(conj(h_tire))*D_L*G1+power1*transpose(G1)*transpose(D_L)*...,

diag(h_tire)*W_tire*diag(conj(h_inversechapeau))*D_R*G1+...,

power1*transpose(G1)*transpose(D_R)*diag(h_inversechapeau)*W_tire*...,

diag(conj(h_inversechapeau))*D_R*G1+power1*transpose(G1)*...,

transpose(D_L)*diag(h_tire)*W_tire*diag(conj(h_tire))*D_L*G1+...,

sigma^2*transpose(G1)*SW*G1;

q2=sigma^2*transpose(h_chapeau)*SW*G1;

k2=sigma^2*(transpose(h_chapeau)*SW*conj(h_chapeau)+1);

b1 = transpose(delt_c)*Q2*conj(delt_c)+2*real(q2*conj(delt_c))+k2;

Q3=power1*transpose(G1)*transpose(D_R)*diag(g_inversechapeau)*W_tire*...,

diag(conj(h_tire))*D_L*G2+power1*transpose(G2)*transpose(D_L)*...,

diag(h_tire)*W_tire*diag(conj(g_inversechapeau))*D_R*G1+...,

power1*transpose(G1)*transpose(D_R)*diag(g_inversechapeau)*W_tire*...,

diag(conj(g_inversechapeau))*D_R*G1+power1*transpose(G2)*...,

transpose(D_L)*diag(h_tire)*W_tire*diag(conj(h_tire))*D_L*G2;

q3=power1*transpose(h_tire.*g_inversechapeau)*W_tire*...,

diag(conj(g_inversechapeau))*D_R*G1+power1*transpose(h_tire.*...,

g_inversechapeau)*W_tire*diag(conj(h_tire))*D_L*G2;

k3=power1*transpose(h_tire.*g_inversechapeau)*W_tire*conj(h_tire.*g_inversechapeau);

a2=transpose(delt_c)*Q3*conj(delt_c)+2*real(q3*conj(delt_c))+k3;

Q4 = power2*transpose(G2)*transpose(D_R)*diag(g_inversechapeau)*W_tire*...,

diag(conj(g_tire))*D_L*G2+power2*transpose(G2)*transpose(D_L)*...,

diag(g_tire)*W_tire*diag(conj(g_inversechapeau))*D_R*G2+...,

power2*transpose(G2)*transpose(D_R)*diag(g_inversechapeau)*W_tire*...,

diag(conj(g_inversechapeau))*D_R*G2+power2*transpose(G2)*...,

transpose(D_L)*diag(g_tire)*W_tire*diag(conj(g_tire))*D_L*G2+...,

sigma^2*transpose(G2)*SW*G2;

q4=sigma^2*transpose(g_chapeau)*SW*G2;

k4=sigma^2*(transpose(g_chapeau)*SW*conj(g_chapeau)+1);

b2=transpose(delt_c)*Q4*conj(delt_c)+2*real(q4*conj(delt_c))+k4;

Q5 = power1*transpose(G3)*transpose(D_L)*diag(h_tire)*W_tire*...,

diag(conj(h_tire))*D_L*G3+power1*transpose(G1)*transpose(D_R)*...,

diag(e_inversechapeau)*W_tire*diag(conj(e_inversechapeau))*D_R*G1+...,

power1*transpose(G3)*transpose(D_L)*diag(h_tire)*W_tire*...,

diag(conj(e_inversechapeau))*D_R*G1+power1*transpose(G1)*...,

transpose(D_R)*diag(e_inversechapeau)*W_tire*diag(conj(h_tire))*D_L*G3+...,

power2*transpose(G3)*transpose(D_L)*diag(g_tire)*W_tire*...,

diag(conj(g_tire))*D_L*G3+power2*transpose(G2)*transpose(D_R)*...,

diag(e_inversechapeau)*W_tire*diag(conj(e_inversechapeau))*D_R*G2+...,

power2*transpose(G3)*transpose(D_L)*diag(g_tire)*W_tire*...,

diag(conj(e_inversechapeau))*D_R*G2+power2*transpose(G2)*...,

transpose(D_R)*diag(e_inversechapeau)*W_tire*diag(conj(g_tire))*D_L*G3;

q5 = power1*transpose(h_tire.*e_inversechapeau)*W_tire*diag(conj(h_tire))*D_L*G3+...,

    power1*transpose(h_tire.*e_inversechapeau)*W_tire*diag(conj(e_inversechapeau))*D_R*G1+...,

    power2*transpose(g_tire.*e_inversechapeau)*W_tire*diag(conj(g_tire))*D_L*G3+...,

    power2*transpose(g_tire.*e_inversechapeau)*W_tire*diag(conj(e_inversechapeau))*D_R*G2;

k5= power1*(abs(he)^2+transpose(h_tire.*e_inversechapeau)*W_tire*conj(h_tire.*e_inversechapeau))+...,

    power2*(abs(ge)^2+transpose(g_tire.*e_inversechapeau)*W_tire*conj(g_tire.*e_inversechapeau));

ae = transpose(delt_c)*Q5*conj(delt_c)+2*real(q5*conj(delt_c))+k5;

% ae ae_expression have been verified

%be = sigma^2*((transpose(e_chapeau+e_wavy)*W)*(transpose(e_chapeau+e_wavy)*W)'+2);

Q6 = sigma^2*transpose(G3)*SW*G3;

q6 = sigma^2*transpose(e_chapeau)*SW*G3;

k6 = sigma^2*(transpose(e_chapeau)*SW*conj(e_chapeau)+2);

be = transpose(delt_c)*Q6*conj(delt_c)+2*real(q6*conj(delt_c))+k6;

QR= power1*transpose(G1)*SPW*G1+power2*transpose(G2)*SPW*G2;

qR=power1*transpose(h_chapeau)*SPW*G1+power2*transpose(g_chapeau)*SPW*G2;

kR = power1*transpose(h_chapeau)*SPW*conj(h_chapeau)+power2*transpose(g_chapeau)*SPW*conj(g_chapeau)+sigma^2*trace(W_tire);

PR = 100;

gamma1_wavy = 2^(2*R1)-1;

gamma2_wavy = 2^(2*R2)-1;

gammae_wavy = 2^(2*Re)-1;

Epsilon_h = 0.001*rand(1);

Epsilon_g = 0.001*rand(1);

Epsilon_e= 0.001*rand(1);

F = set(W_tire>0);

F = F+set(Lambda1>0);

F = F+set(Lambda2>0);

F = F+set(Lambda3>0);

F = F+set(Lambda4>0);

F = F+set(Lambda5>0);

F = F+set(Lambda6>0);

F = F+set(Lambda7>0);

F = F+set(Lambda8>0);

F = F+set(Lambda9>0);

F = F+set(R1>0);

F = F+set(R2>0);

F = F+set(Re>0);

F = F + set([Q1-gamma1_wavy*Q2+Lambda1*transpose(G1)*G1+Lambda2*transpose(G2)*G2 q1'-gamma1_wavy*q2';q1-gamma1_wavy*q2 k1-gamma1_wavy*k2-Lambda1* Epsilon_h^2-Lambda2*Epsilon_g^2]>0);

F = F + set([Q3-gamma2_wavy*Q4+Lambda3*transpose(G1)*G1+Lambda4*transpose(G2)*G2  q3'-gamma2_wavy*q4';q3-gamma2_wavy*q4  k3-gamma2_wavy*k4-Lambda3*Epsilon_h^2-Lambda4*Epsilon_g^2] >0);

F = F + set([Q5-gammae_wavy*Q4+Lambda5*transpose(G1)*G1+Lambda6*transpose(G2)*G2+Lambda7*transpose(G3)*G3  q5'-gammae_wavy*q6';q5-gammae_wavy*q6  k5-gammae_wavy*k4-Lambda5*Epsilon_h^2-Lambda6*Epsilon_g^2-Lambda7*Epsilon_e^2]>0);

F = F + set([-QR+Lambda8*transpose(G1)*G1+Lambda9*transpose(G2)*G2  -qR';qR PR-kR-Lambda8*Epsilon_h^2-Lambda8*Epsilon_g^2] > 0);

ops = sdpsettings('solver','sedumi','verbose',0);

solvesdp(F,-(R1+R2-Re),ops);

###################### test1 function###########

function [he,ge,h_wavy,h_chapeau,g_wavy,g_chapeau,e_wavy,e_chapeau,number,power1,power2,Sigma] = test1()

clc;

clear;

% generation of two single-antenna source nodes 

number = 8;

s1 = rand(1);

s2 = rand(1);

Sigma = 1;

 

hr = -1+2*rand(number,1);

hi = -1+2*rand(number,1);

h = complex(hr,hi);

gr = -1+2*rand(number,1);

gi = -1+2*rand(number,1);

g = complex(gr,gi);

her= rand(1);

hei= rand(1);

he = complex(her,hei);

ger= rand(1);

gei= rand(1);

ge = complex(ger,gei);

e = rand(number,1);

power1 = 10;

power2 = 10;

%  Wr = randn(number) ;

% % 

%  Wi = randn(number) ;

% % 

%  W = complex(Wr, Wi);

Eta = Sigma*rand(number,1);

Etae = Sigma*rand(1);

n1 = Sigma*rand(1);

n2 = Sigma*rand(1);

ne = Sigma*rand(1);

% the signals received at R 

x= sqrt(power1)*h*s1 + sqrt(power2)*g*s2 + Eta;

% the signals received at E

% xe = sqrt(power1)*he*s1 + sqrt(power2)*ge*s2 + Etae;

% 

% y1 = transpose(h)*W*x + n1;

% y2 = transpose(g)*W*x + n2;

% y3 = transpose(e)*W*x + ne;

% firstly represent the errors in the CSI which have bounded norms

% in the beginning, we assumed that the generated matrix fit the bounded

% norms, otherwise, we need to re-generate the martix until it reachs the

% bounds norm.

sigma_error = 0.0002;

h_wavyr =  sigma_error*randn(number,1);

h_wavyi =  sigma_error*randn(number,1);

h_wavy = complex(h_wavyr,h_wavyi);

h_chapeau = h+h_wavy;

g_wavyr =  sigma_error*rand(number,1);

g_wavyi =  sigma_error*rand(number,1);

g_wavy = complex(g_wavyr,g_wavyi);

g_chapeau = g+g_wavy;

e_wavyr =  sigma_error*rand(number,1);

e_wavyi =  sigma_error*rand(number,1);

e_wavy = complex(e_wavyr,e_wavyi);

e_chapeau = e+e_wavy;

 end

Thank you very much for your help and your time,

[yangke Sun]

Hi, Johan, thank you in advance for your help, It is a robust Two way relay beam forming 

[Johan Löfberg]

If it says "Warning: Solver not found (sedum)", well then obviously you haven't installed SeDuMi properly (i.e installed the correct binaries, and made those visible by adding suitable directories to your MATLAB path.

BTW, don't use warning('off','YALMIP:strict'). Fix your code instead.
http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Blog.Prepare-your-code

...and SET is obsolete
http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Commands.set

...and as a basic rule, never turn off display until you know everything works

...and vectorize. Instead of lambd1, lambda2, lambda3 etc, why not define a vector and update code accordingly

lambda=sdpvar(9,1),
F = [F, lambda >=0]

However, the problem you have cannot be solved using SeDuMi as it is not linear in your decision variables. E.g, this is bilinear

[yangke Sun]

Hi johan, 

Thank you for your time and your answer. It is really helpful.

However, when you said that the problem cannot be solved, because it is bilinear. Do you have any suggestion for doing this? 

Thanks!

Yangke

[Johan Löfberg]

Not really. It is a much harder (close to intractable) problem. Some iterative scheme, fixing some variables, path-following, linearization etc