Hi All,
I am trying to implement the following SDP code in Matlab
clc;
clear all;
close all;
% constants
nsvd = 5; %Number of modes
N_c = 10; %Number of linearly spaced points between 0 and c_max
N_y = 80; %Number of points in y-direction
yP = [0,0.390514938607723,1.46799561745170,3.26332937544719,5.80864526449735,9.03486538711407, ...
12.9515731316936,17.5988700160186,22.6814498827928,28.5869213352291,35.2136148497178, ...
42.6910473000038,50.7978486141225,59.6189626961894,69.2381708298159,79.3724706248449, ...
90.3434973508633,101.808629345078,114.396622785589,126.780118224767,140.158733596590, ...
155.376327116633,170.142777579874,185.191311107031,201.642812698034,218.500954435884, ...
236.105621216221,254.345928171713,273.111292646320,292.887570038162,312.697552616828, ...
332.148735619785,354.602093942407,377.282264867448,398.929900912017,421.336795287591, ...
444.796887896479,469.218039646431,493.509334440579,519.367097025344,544.129893364562, ...
569.381956078757,596.118322292992,623.373679388788,652.190251910176,679.896486024956, ...
708.041489801016,736.529324212737,765.441274927007,796.669335478415,826.371678357109, ...
856.463485961769,886.803585576631,917.363265930359,948.273568432540,982.053780249119, ...
1014.82623722978,1047.92463113639,1081.33604090463,1114.21828509537,1147.37354925624, ...
1180.79557134070,1215.52574717483,1248.74513364622,1282.18717279486,1315.83880940303, ...
1353.17618416708,1388.13333366187,1423.26506828233,1458.55767318165,1493.68520322801, ...
1527.87693735191,1562.18340704703,1596.59121963631,1631.08694288010,1671.06992968758, ...
1708.49317921569,1745.96768177412,1783.47880793829,1821.01191398637]';
log_yP = log(yP);
log_yP(1) = 0;
A = (10^7)*(rand(nsvd,nsvd,N_y,N_c,4) - 0.5); % A(i,j,y,c,:) --> Column 1:uu; 2:vv; 3:ww; 4:uv;
%CVX Code
cvx_begin SDP
variable X_l(nsvd,nsvd,N_c) hermitian %X_l is the weight matrix
variable e; %Variable to be optimized
dual variable Q{4};
expressions num_int(N_y,4) den_int(N_y,4) num(N_y,4) sum_2(N_y,4)
expression E_dns(N_y,4)
E_dns = (rand(N_y,4)-1/2)*(10^7); %E_dns(y,:): Column 1:uu; 2:vv; 3:ww; 4:uv;
% ||g||^2 norm --> integrate {|g log(yP)|^2 d(log(yP))} from yP_min to yP_max
% numm and denn are the norms calculated for the numerator and denominator respectively
% denn is the norm from E_dns
% Variables unaffected by CVX - START
den_int(:,1) = power(abs(E_dns(:,1).*log_yP),2); %For uu
den_int(:,2) = power(abs(E_dns(:,2).*log_yP),2); %For vv
den_int(:,3) = power(abs(E_dns(:,3).*log_yP),2); %For ww
den_int(:,4) = power(abs(E_dns(:,4).*log_yP),2); %For uv
denn(1) = diff(yP,1).' * (den_int(1:N_y-1,1) + den_int(2:N_y,1))/2; %trapz function for uu
denn(2) = diff(yP,1).' * (den_int(1:N_y-1,2) + den_int(2:N_y,2))/2; %trapz function for vv
denn(3) = diff(yP,1).' * (den_int(1:N_y-1,3) + den_int(2:N_y,3))/2; %trapz function for ww
denn(4) = diff(yP,1).' * (den_int(1:N_y-1,4) + den_int(2:N_y,4))/2; %trapz function for uv
% Variables unaffected by CVX - END
minimize(e) %Minimizing error 'e'
subject to
for m = 1,N_y;
for k = 1,N_c;
sum_2(m,1) = sum_2(m,1) + real(trace(A(:,:,m,k,1)*X_l(:,:,k))); %integral term in numerator for uu
sum_2(m,2) = sum_2(m,2) + real(trace(A(:,:,m,k,2)*X_l(:,:,k))); %integral term in numerator for vv
sum_2(m,3) = sum_2(m,3) + real(trace(A(:,:,m,k,3)*X_l(:,:,k))); %integral term in numerator for ww
sum_2(m,4) = sum_2(m,4) + real(trace(A(:,:,m,k,4)*X_l(:,:,k))); %integral term in numerator for uv
end
num(m,1) = E_dns(m,1) - sum_2(m,1); %Numerator before norm for uu
num(m,2) = E_dns(m,2) - sum_2(m,2); %Numerator before norm for vv
num(m,3) = E_dns(m,3) - sum_2(m,3); %Numerator before norm for ww
num(m,4) = E_dns(m,4) - sum_2(m,4); %Numerator before norm for uv
end
num_int(:,1) = power(abs(num(:,1).*log_yP),2); %For uu
num_int(:,2) = power(abs(num(:,2).*log_yP),2); %For vv
num_int(:,3) = power(abs(num(:,3).*log_yP),2); %For ww
num_int(:,4) = power(abs(num(:,4).*log_yP),2); %For uv
numm(1) = diff(yP,1).' * (num_int(1:N_y-1,1) + num_int(2:N_y,1))/2; %trapz function for uu
numm(2) = diff(yP,1).' * (num_int(1:N_y-1,2) + num_int(2:N_y,2))/2; %trapz function for uu
numm(3) = diff(yP,1).' * (num_int(1:N_y-1,3) + num_int(2:N_y,3))/2; %trapz function for uu
numm(4) = diff(yP,1).' * (num_int(1:N_y-1,4) + num_int(2:N_y,4))/2; %trapz function for uu
% X_l is positive semidefinite
for k = 1:N_c
X_l(:,:,k) >= 0; %Constraint numbers 1-N_c
end
% Inequality condition for error
for k = 1:4
Q{k} : numm(k)/denn(k) <= e; %Constraint numbers N_c+1 to N_c+4
end
cvx_end
While the problem is solved, the solver uses only 3 equality constraints whereas, there are more constraints defined in the code. In the code, the constraints display this message when I hover over them “produces a value that might be unused”.
I would greatly appreciate any suggestions and advise to figure out the mistake I am making.
Thanks
Vikram