function [rbar, t, tau, W1] =Optimize_R(T, eta, P1, N, rho1, rho2, rho3, si1, si2)
cvx_begin sdp
cvx_precision high
variable t(4) nonnegative;
variable tau(2) nonnegative;
variable W1(N+1,N+1) complex semidefinite; %W=t1*V1 >= 0
dual variables y1 y2 y3 y4;
expression R12;
expression R02;
expression R13;
expression R1;
expression R2;
%Achievable rate from WD1 to WD2 and The HAP in phase 2 are:
R12=-rel_entr(t(2,1),t(2,1)+real(rho1trace(si1W1)))/log(2); %equation 27
%xlog2(1+y/x)=-rel_entr(x,x+y)/log(2)
R02=-rel_entr(t(2,1),t(2,1)+real(rho2trace(si1*W1)))/log(2); %equation 28
%Achievable rate of WD2 relaying info of WD1 to The HAP in phase 3:
R13=-rel_entr(t(3,1),t(3,1)+real(rho3tau(1,1)))/log(2); %equation 22
%Overall achievable Rates of WD1 and WD2 are:
R1=min(R12,R02+R13); %equation 12
R2=-rel_entr(t(4,1),t(4,1)+real(rho3tau(2,1)))/log(2); %equation 23
rbar = min(R1, R2); %Auxilary Variable
maximize rbar;
subject to
t(1,1)+t(2,1)+t(3,1)+t(4,1)<=T; %(1)
t(1,1) >= 0;
t(2,1)>= 0;
t(3,1)>= 0;
t(4,1)>= 0;
tau(1,1)>= 0;
tau(2,1) >= 0;
diag(W1)==t(1,1); %|[W]n,n|=t1 for n=1,2,...N+1 (26)
W1 >= 0;
y1: -rbar<= R02+R13;
y2: -rbar<=R12;
y3: -rbar<=R2;
rbar>=0;
y4: tau(1,1)+tau(2,1)<=eta*P1*real(trace(si2*W1)); %(29)
W1==semidefinite(N+1);
cvx_end