how to solve for fsolve

how can i degug these code? %optimal gamma_j (SNR of user j)

y = sym('y');
eqn = (1 + y).^L_i - (1 + gamma_i ./ (y + 1)).^L_j;
func = matlabFunction(eqn,'Vars',{y}); 
xo=[0,0];
option = optimset('Display','off');
yj = fsolve(func ,[x0],option);  %#ok<NBRAK>
gamma_j = min(yj, Gamma_j);

Answer 1

[Walter Roberson]

    y = sym('y');

That is a scalar symbol

eqn = (1 + y).^L_i - (1 + gamma_i ./ (y + 1)).^L_j;

We cannot tell how many items that expands to. We can see that it does not distinguish elements of y, such as by trying to use y(2). We cannot prove from that whether y is expected to be scalar or vector inside that context. Nothing in that line contradicts the possibility that y is scalar.

func = matlabFunction(eqn,'Vars',{y});

That is going to produce a function handle in which each time y is unindexed. If any of the values such as gamma_i are non-scalar then eqn would have expanded to a non-scalar before getting to matlabFunction, and y would be used in full at each location. For example

eqn = y*v

where v=[2, 3] would already have expanded to eqn=[y*2,y*3] and when matlabFunction of that is taken, at each point all of y would be used. It would not expand to [y(1)*2,y(2)*3] - not with scalar symbol y.

 xo=[0,0];

but you are passing a vector in for y at the fsolve level. However many elements were implied by gamma_i and so on, you are going to get an output twice as wide as you expect. Except if one of the locations came out to be independent of y... in that case you would probably end up with an error.

[Martin Njagi]

%Let values to K users and F channels
K = 30;
F = 5;

range = 1:1;
maxCMPLTime = zeros(length(range),1);

for i = 1:length(range)
    %random user i location
    mag_min_i = 0;
    mag_max_i = 500; %radius R = 500 m
    mag_i = (mag_max_i - mag_min_i).*sqrt(rand(K/2,1)); %random magnitude
    ang_min_i = 0;
    ang_max_i = 360; %within 360 degree
    ang_i = (ang_max_i - ang_min_i).*rand(K/2,1); %random angle
    d_meter_i = sqrt((mag_i.*cos(ang_i)).^2 + (mag_i.*sin(ang_i)).^2);
    d_i = d_meter_i / 1000;
    %Calculate path loss
    PL_i = -(128.1 + 37.6 .* log10(d_i)); %where d is in km
    
    %random user j location
    mag_min_j = 0;
    mag_max_j = 500; %radius R = 500 m
    mag_j = (mag_max_j - mag_min_j).*sqrt(rand(K/2,1)); %random magnitude
    ang_min_j = 0;
    ang_max_j = 360; %within 360 degree
    ang_j = (ang_max_j - ang_min_j).*rand(K/2,1); %random angle
    d_meter_j = sqrt((mag_j.*cos(ang_j)).^2 + (mag_j.*sin(ang_j)).^2);
    d_j = d_meter_j / 1000;
    %Calculate path loss
    PL_j = -(128.1 + 37.6 .* log10(d_j)); %where d is in km
    
    XY = randn + 1i*randn;
    %Calculate channel gain
    h_i = 10.^(PL_i./10) * (abs(XY/sqrt(2))).^2;
    h_j = 10.^(PL_j./10) * (abs(XY/sqrt(2))).^2;
    
    W = (180 * 10^3) / F; %bandwidth of each channel
    P_dBm = 20; %maximum transmit power of each user
    P_watt = 10 ^ ((P_dBm - 30) / 10);
    N0_dBm = -174; %N0 in dBm/Hz
    N0_WHz = 10 ^ ((N0_dBm - 30) / 10);
    %noise power of each channel
    sigma_sqr = N0_WHz * W; 
    
    %packet size (in bits)
    L_i = randi([8*8, 1024*8], K/2, 1);
    L_j = randi([8*8, 1024*8], K/2, 1);
    
    %optimal gamma_i (SNR of user i)
    gamma_i = h_i * P_watt / sigma_sqr;
    
    Gamma_j = h_j * P_watt / sigma_sqr;
    %optimal gamma_j (SNR of user j)
    y = sym('y');
    eqn = (1 + y).^L_i - (1 + gamma_i ./ (y + 1)).^L_j;
    func = matlabFunction(eqn,'Vars',{y}); 
    
    option = optimset('Display','off');
    yj = fsolve(func ,0,option);  
    gamma_j = min(yj, Gamma_j);
    
    %data rate of user i and j
    R_i = W * log2(1 + gamma_i / (gamma_j + 1));
    R_j = w * log2(1 + gamma_j);
    
    %transmission time
    tau_i = L_i ./ R_i;
    tau_j = L_j ./ R_j;
    
    s_i = sort(tau_i);
    s_j = sort(tau_j);
    
    %Tau_ij = [max(tau_i), max(tau_j)];
    
    %NOMA-SPT scheduling simulation
    
    for  i = 1:K/2F; %#ok<FXSET,SEPEX>
        
        for ii = 1:F;
            CMPLTime(ii) = [max(s_i), max(s_j)]; %#ok<SAGROW>
            disp(ch(ii));
            CMPLTime(ii) = max(ch(ii)); %#ok<SAGROW>
        end
    end
    
    %maximum completion time
    maxCMPLTime(i) = max(CMPLTime);
end

Error using fsolve
Objective function is returning undefined values at initial point. FSOLVE cannot continue.

%list of all max transmission time after looping
allMaxTime = maxCMPLTime(1:length(range));
%average of above list data
average = mean(allMaxTime);

equation = func(y); 
eq = (1 + y).^L_i - (1 + gamma_i ./ (y + 1)).^L_j;

These is the whole code which i wanted to use it in NoMA.SPT

[Martin Njagi]

just wanted some steps for debugging the code.Am still an amateur in Matlab.Solution to these will be savior to my academic project

[Walter Roberson]

%Let values to K users and F channels
K = 30;
F = 5;

range = 1:1;
maxCMPLTime = zeros(length(range),1);

for i = 1:1
    %random user i location
    mag_min_i = 0;
    mag_max_i = 500; %radius R = 500 m
    mag_i = (mag_max_i - mag_min_i).*sqrt(rand(K/2,1)); %random magnitude
    ang_min_i = 0;
    ang_max_i = 360; %within 360 degree
    ang_i = (ang_max_i - ang_min_i).*rand(K/2,1); %random angle
    d_meter_i = sqrt((mag_i.*cos(ang_i)).^2 + (mag_i.*sin(ang_i)).^2);
    d_i = d_meter_i / 1000;
    %Calculate path loss
    PL_i = -(128.1 + 37.6 .* log10(d_i)); %where d is in km
    
    %random user j location
    mag_min_j = 0;
    mag_max_j = 500; %radius R = 500 m
    mag_j = (mag_max_j - mag_min_j).*sqrt(rand(K/2,1)); %random magnitude
    ang_min_j = 0;
    ang_max_j = 360; %within 360 degree
    ang_j = (ang_max_j - ang_min_j).*rand(K/2,1); %random angle
    d_meter_j = sqrt((mag_j.*cos(ang_j)).^2 + (mag_j.*sin(ang_j)).^2);
    d_j = d_meter_j / 1000;
    %Calculate path loss
    PL_j = -(128.1 + 37.6 .* log10(d_j)); %where d is in km
    
    XY = randn + 1i*randn;
    %Calculate channel gain
    h_i = 10.^(PL_i./10) * (abs(XY/sqrt(2))).^2;
    h_j = 10.^(PL_j./10) * (abs(XY/sqrt(2))).^2;
    
    W = (180 * 10^3) / F; %bandwidth of each channel
    P_dBm = 20; %maximum transmit power of each user
    P_watt = 10 ^ ((P_dBm - 30) / 10);
    N0_dBm = -174; %N0 in dBm/Hz
    N0_WHz = 10 ^ ((N0_dBm - 30) / 10);
    %noise power of each channel
    sigma_sqr = N0_WHz * W; 
    
    %packet size (in bits)
    L_i = randi([8*8, 1024*8], K/2, 1);
    L_j = randi([8*8, 1024*8], K/2, 1);
    
    %optimal gamma_i (SNR of user i)
    gamma_i = h_i * P_watt / sigma_sqr;
    
    Gamma_j = h_j * P_watt / sigma_sqr;
    %optimal gamma_j (SNR of user j)
    y = sym('y');
    eqn = (1 + y).^L_i - (1 + gamma_i ./ (y + 1)).^L_j;
    func = matlabFunction(eqn,'Vars',{y}); 
    
    option = optimset('Display','off');
  
    vpa((eqn))
        
end

ans = 

Look at those values.

(y+1)^7196 . In order to have a reasonable chance of balancing the other side, that cannot exceed roughly 10^10, which requires that y be less than about 0.002 .

What about the other side? (44168/(y+1) + 1)^6248 and to balance that against (say) 10^10 on the other side, requires y pretty close to 11962791.35361

You also have a series of 15 equations to solve, but only one single y value.

[Martin Njagi]

thanks so much for the idea.Atleast I can now work on it.

[Walter Roberson]

I tried solving an approximation to get an order-of-magnitude estimate for the solution ranges. For the particular one I worked on I was able to tell that the solution was about 504.3 . However, the slope near there was over 10^2100, so there is absolutely no way it could be solved with double precision.

Using a different software package I was able to get a solution fairly quickly when I pushed the number of symbolic digits up to about 3000. MATLAB is unfortunately considerably slower at the operation. Also, finding a starting point it was willing to work with took me a while. For another one from the same batch, 20000 digits of precision was not enough to get a solution in that software package (the slope was above 10^6000)

Finding a solution just might not be practical. But it might be easier to find the minima of the square of the elements.

[Martin Njagi]

%data rate of user i and j
R_i = W * log2(1 + gamma_i / (gamma_j + 1));
R_j = w * log2(1 + gamma_j);

%transmission time
tau_i = L_i ./ R_i;
tau_j = L_j ./ R_j;

s_i = sort(tau_i);
s_j = sort(tau_j);

%Tau_ij = [max(tau_i), max(tau_j)];

%NOMA-SPT scheduling simulation

for  i = 1:K/2F; %#ok<FXSET,SEPEX>
    
    for ii = 1:F;
        CMPLTime(ii) = [max(s_i), max(s_j)]; %#ok<SAGROW>
        disp(ch(ii));
        CMPLTime(ii) = max(ch(ii)); %#ok<SAGROW>
    end
end

%maximum completion time
maxCMPLTime(i) = max(CMPLTime);
end

Error using fsolve

Objective function is returning undefined values at initial point. FSOLVE cannot continue.

%list of all max transmission time after looping
allMaxTime = maxCMPLTime(1:length(range));
%average of above list data
average = mean(allMaxTime);
equation = func(y); 
eq = (1 + y).^L_i - (1 + gamma_i ./ (y + 1)).^L_j; 

these piece gives the details on the transmission time and data rates .How ca I achieve it?

[Walter Roberson]

I do not see anything in that description that relates to the equation you made ??

[Martin Njagi]

the code is just a trial case in Completion_time_Driven Scheduling for the Uplink of NOMA -enabled NB-IoTNetworks

[Walter Roberson]

I do not see any Mathworks example by that name.