Unable to solve a geometric program using YALMIP

[Rukhsana Ruby]

Hi,

I downloaded YALMIP R20130213 yesterday. I wrote a program to solve a geometric formulation. Program is:

function [] = gp_test2()

sdpvar Es Er y t

Gsd = 0.1;

Gsr = 0.5;

Grd = 0.5;

obj = 1/t;

F = set((t^-1 + Es*Gsd*t^-1 + Es*Gsr*Er*Grd*y^-1*t^-1) <= 1);

F = F + set((Es*Gsr*y^-1 + Er*Grd*y^-1 + y^-1) <= 1);

F = F + set(Er <= 20) + set((Es + Er) <= 35);

s = solvesdp(F, obj, sdpsettings('solver','fmincon-geometric'))

double(Es)

double(Er)

end

I got the following error message:

>> gp_test2

Warning: Trust-region-reflective method does not currently solve this type of

problem,

 using active-set (line search) instead. 

> In fmincon at 437

  In callfmincongp at 132

  In solvesdp at 340

  In gp_test2 at 16

Maximum number of function evaluations exceeded;

 increase OPTIONS.MaxFunEvals.

s = 

    yalmiptime: 0.110963000000002

    solvertime: 0.835329999999999

          info: 'Maximum iterations or time limit exceeded (fmincon-geometric)'

       problem: 3

ans =

     2.227061241501988e-04

ans =

     3.293217186950261e-04

And these are not expected output at all. Please help me in solving this problem as soon as possible.

Regards

Ruby 

[Johan Löfberg]

First, SET is obsolete. Don't use that.

Secondly, using fmincon is a bit annoying, since they seem to change the default algorithm in different version. 'interior-point' seems to work fine here

sdpvar Es Er y t
Gsd = 0.1;
Gsr = 0.5;
Grd = 0.5;
obj = 1/t;

F = t^-1 + Es*Gsd*t^-1 + Es*Gsr*Er*Grd*y^-1*t^-1 <= 1;
F = [F, Es*Gsr*y^-1 + Er*Grd*y^-1 + y^-1 <= 1];
F = [F, Er <= 20, Es + Er <= 35];
s = solvesdp(F, obj, sdpsettings('solver','fmincon-geometric','fmincon.Algorithm','interior-point'))

[Rukhsana Ruby]

Hi Johan,

Thanks for your quick reply. But the output from this program is not expected. Therefore, I have modified 2nd constraint, it should be either "==" or ">=". 

sdpvar Es Er y t

Gsd = 0.1;

Gsr = 0.5;

Grd = 0.5;

obj = 1/t;

F = t^-1 + Es*Gsd*t^-1 + Es*Gsr*Er*Grd*y^-1*t^-1 <= 1;

F = [F, Es*Gsr*y^-1 + Er*Grd*y^-1 + y^-1 == 1];

F = [F, Er <= 20, Es + Er <= 35];

s = solvesdp(F, obj, sdpsettings('solver','fmincon-geometric','fmincon.Algorithm','interior-point'))

However, I got the following error:

s = 

    yalmiptime: 0.107278000000000

    solvertime: 0

          info: 'Solver not applicable (fmincon-geometric)'

       problem: -4

ans =

   NaN

ans =

   NaN

Please suggest me to get rid of this error.

Regards

Ruby

[Johan Löfberg]

It is not a simple posynomial geometric program any longer when you have such an equality constraint.

I would proceed by multiplying out the inverse terms to obtain a polynomial program, and then apply the global solver available in YALMIP (bmibnb)

sdpvar Es Er y t

Gsd = 0.1;

Gsr = 0.5;

Grd = 0.5;

obj = t;

F = y + Es*Gsd*y + Es*Gsr*Er*Grd <= y*t;

F = [F, Es*Gsr + Er*Grd + 1 == y];

F = [F, Er <= 20, Es + Er <= 35];

F = [-100 <= [Es Er y t] <= 100];

s = solvesdp(F, -obj, sdpsettings('solver','bmibnb'))

The solution here seems to indicate that the problem does not have any finite optima, i.e, I can generate arbitrarily good solutions (t tends to infinity, and the rest goes to zero in suitable order)

[Johan Löfberg]

Oops, that should be (adding some bounds for the global solver to be applicable)

F = [F, -100 <= [Es Er y t] <= 100];

[Rukhsana Ruby]

Hi Johan,

Thanks for your reply again. I got the following your error from your given program

* Starting YALMIP global branch & bound.

* Branch-variables : 4

* Upper solver     : fmincon

* Lower solver     : LINPROG

* LP solver        : LINPROG

 Node       Upper      Gap(%)       Lower    Open

s = 

    yalmiptime: NaN

    solvertime: NaN

          info: [1x112 char]

       problem: 9

ans =

   NaN

ans =

   NaN

My original problem was like this:

sdpvar Es Er

obj = 1/(1 + Es*Gsd  +  (Es*Gsr*Er*Grd)/(Es*Gsr + Er*Grd + 1));

F = [ Er <= 20, Es + Er <= 35];

Now i wanted to convert it to gp. Therefore I have replaced denominator of "obj" as "t" and then again the second part of denominator is monomial divided by posynomial which is eventually not posynomial. Therefore, I have replaced (Es*Gsr + Er*Grd + 1) by "y". Thus, the resultant problem looked like that. Furthermore, I want to approximate this problem somehow to GP and this is my goal. Please help me if you have any suggestion on it. I appreciate your quick reply as my deadline for it is very nearby.

Regards

Ruby 

[Johan Löfberg]

Something has crashed hard inside solvesdp (error code 9). Turn on debug to see where it fails

s = solvesdp(F, -obj, sdpsettings('solver','bmibnb','debug',1))

I tried with fmincon as local upper solver and linprog as lp and lower bound solver as you have, and it worked here. What OS and MATLAB version are you running

[Rukhsana Ruby]

Hi Johan,

Thanks again for your reply. I tried the program in a matlab with older version. I will try this tomorrow again in newer version. Actually, I found a way to solve this problem using GP. Essentially, objective function is posynomial divided by posynomial which is not GP to solve. I have to somehow approximate the problem in GP which is usually named as complementary GP, apply single/double condensation in order to solve the approximated problem in iterative way. I will let you know the update soon. 

Regards

Ruby

[Johan Löfberg]

Your conversion is wrong. To make the original model equivalent to your reformulation involving the auxiliary variables, you must reverse the first inequality, or, since I use a global polynomial approach here, just use an equality. The globally optimal solution is found in three iterations in bmibnb

obj = t;

F = y + Es*Gsd*y + Es*Gsr*Er*Grd == y*t;

F = [F, Es*Gsr + Er*Grd + 1 == y];

F = [F, Er <= 20, Es + Er <= 35];

F = [F,0 <= [Es Er y t] <= 1000];

s = solvesdp(F, -obj, sdpsettings('solver','bmibnb'))

>> double(Er)

ans =

   13.8000

>> double(Es)

ans =

   21.2000

[Rukhsana Ruby]

Hi Johan,

Thank you very much for your effort. I am wondering "obj" is 1/t instead of t. That preserves the original setting. Also, could you please tell me what kind of solver is used by bmibnb? I have to wait for the morning to check my program.

Thanks

Ruby

[Rukhsana Ruby]

Hi Johan,

Also, I am wondering, current problem or conversion is not convex. Then result obtained by any kind of solver should stuck to local minima. Although essentially, if we put equality, it is still GP form.

Regards

Ruby

[Johan Löfberg]

I want to avoid division of variables, since many solvers have problems with this (if x at any trial happens to be 0, the solver will have to deal with things like Inf and NaN, and many solvers crash then).

Hence, instead of minimizing 1/t, we can just as well maximize t, which is equivalent to minimizing -t. Hence the objective -t (YALMIP assumes you minimize)

[Johan Löfberg]

The original problem seems to be convex

[Es,Er] = meshgrid(0.01:35, 0.01:20);

Gsd = 0.1;

Gsr = 0.5;

Grd = 0.5;

O = 1./(1 + Es*Gsd  +  (Es*Gsr.*Er*Grd)./(Es*Gsr + Er*Grd + 1));

mesh(Es,Er,O)

However, if you simply plug the fractional objective into YALMIP, YALMIP will make some internal reformulation (essentially to the model I've described above) in order to remove expressions inside fractions etc, in order to put the model in a normalized format for which YALMIP easily can compute Jacobians etc. This reformulation can destroy convexity (as it introduces the nonlinear equalities we use above)

The solver bmibnb I use is a global solver, hence it returns the globally optimal solution, despite nonconvexity. The price you pay is the possibility that it is very slow. For your problem though, it is solved without effort (probably since the underlying problem actually is convex)

BTW, With GP form, it is almost universally understood that it is a GP which can be formulated as a convex problem using the standard exponential/logarithmic reformulation, hence a model with equality constraints involving posynomials is not what most people would call standard GP form (standard GP form allows equalities with single monomial terms)

[Rukhsana Ruby]

Hi Rohan,

Thank you very much again for your kind effort. This is simplified problem i wanted to solve, therefore it has only 2 variables. But the original problem is very big, it has at least 10 variables. I have to wait for the morning to write code. Also, I am wondering, according  to GP specification, posynomial constraint can be both equality or less than equality, that means equality is permitted on the posynomial constraint. I see your spec that, bmibnb uses branch and bound technique to find the global solution.

Regards

Ruby

[Johan Löfberg]

A posynomial program with general posynomial equalities cannot be reformulated to a convex program as far as I know. The equalities may only be monomial equalities. (There are some special cases with equalities, but I think that is the case when the equalitiy is linear and has some structure in relation to the overall problem)

See here for instance for a typical definition

http://en.wikipedia.org/wiki/Geometric_programming

[Rukhsana Ruby]

Hi Johan,

Could you please tell me how you have concluded that original problem is convex? I checked myself and found it is non-convex.

Thanks

Ruby

[Johan Löfberg]

I looked at the plot generated above (I assume Es and Er are non-negative)

[Johan Löfberg]

Scratch that, zooming in around zero shows non-convexity. Sorry for confusing you. Summary: It is close to convex, and that is probably why the global solver easily finds and proves global optimality.

[Es,Er] = meshgrid(0:0.1:2, 0:0.1:2);

[Rukhsana Ruby]

Hi Johan,

Thanks. Yes, Es and Er must be greater than or equal to zero. I calculated the hessian of this objective function. Moreover, according to the geometric definition, objective function has a term which monomial divided by posynomial results in having non-convex property.

Thanks

Ruby

[Johan Löfberg]

Monomial divided by posynomial is not necessarily non-convex as I can see (but neither is it guaranteed to be convex)

More importantly though, convex and non-convexity is not really the important feature when dealing with GPs. The important feature is if they can be transformed to a problem which is convex, through the use of logarithmic variables changes.

[Rukhsana Ruby]

Ok, Thanks!

[Rukhsana Ruby]

Hi Rohan,

I solved the problem using GP. I just jave used one tricks for solving generalized GP mentioned in Boyed's GP tutorial. Code is below. Also, the result obtained by bmibnb is not globally optimal.

function [] = gp_test2()

sdpvar Es Er y t

Gsd = 0.5;

Gsr = 0.5;

Grd = 0.1;

obj = 1/t;

F = (1 + Es*Gsd + Es*Gsr*Er*Grd*y^-1) <= t;

F = [F, (Es*Gsr + Er*Grd + 1) <= y];

F = [F, Er <= 20, (Es + Er) <= 35];

F = [F,0 <= [Es Er y t] <= 100000];

s = solvesdp(F, obj, sdpsettings('solver','fmincon-geometric', 'fmincon.Algorithm', 'interior-point'))

double(Es)

double(Er)

double(t)

double(y)

obj = 1/Es + 1/Er;

F = (1 + Es*Gsd + Es*Gsr*Er*Grd*y^-1) <= t;

F = [F, (Es*Gsr + Er*Grd + 1) <= y];

F = [F, Er <= 20, (Es + Er) <= 35];

F = [F,0 <= [Es Er y t] <= 100000];

F = [F, 1/t <= 1/100000];

s = solvesdp(F, obj, sdpsettings('solver','fmincon-geometric', 'fmincon.Algorithm', 'interior-point'))

double(Es)

double(Er)

res = 1 + double(Es)*Gsd + (double(Es)*Gsr*double(Er)*Grd)/(double(Es)*Gsr + double(Er)*Grd + 1)

end

Thanks for your all helps. I appreaciate it.

Regards

Ruby

[Johan Löfberg]

To begin with, I don't know where the number 13.8 and 21.2 came from. I must have done something different and reported the wrong numbers. The solution I get now is Es=35 and Er = 0. It yields the objective 0.054054

Your first model is obviously flawed

sdpvar Es Er y t

Gsd = 0.5;

Gsr = 0.5;

Grd = 0.1;

obj = 1/t;

F = (1 + Es*Gsd + Es*Gsr*Er*Grd*y^-1) <= t;

F = [F, (Es*Gsr + Er*Grd + 1) <= y];

F = [F, Er <= 20, (Es + Er) <= 35];

F = [F,0 <= [Es Er y t] <= 100000];

Since you effectively maximize t, the optimal solution is trivially infinity (or 100000 which you have added as an upper bound for some reason). Hence, the first constraint is theoretically redundant since t goes to infinity. Hence, you simply compute any solution which satisfies

F = [(Es*Gsr + Er*Grd + 1) <= 100000];

F = [F, Er <= 20, (Es + Er) <= 35];

F = [F,0 <= [Es Er] <= 100000];

And you get some arbitrary solution satisfying this

The second model is just as weird. First, I don't see how the objective in any sense relates to the original model. Secondly, you say t<=1000000 and 1/t <= 1/100000, which means you are saying t == 100000, hence the model is

obj = 1/Es + 1/Er;

F = (1 + Es*Gsd + Es*Gsr*Er*Grd*y^-1) <= 100000;

F = [F, (Es*Gsr + Er*Grd + 1) <= y];

F = [F, Er <= 20, (Es + Er) <= 35];

F = [F,0 <= [Es Er y] <= 100000];

s = solvesdp(F, obj, sdpsettings('solver','fmincon-geometric', 'fmincon.Algorithm', 'interior-point'))

The solution I get from your model here is Er=Es=17.5, which is sub-optimal since the objective (the one you talked about initially) evaluates to 0.09024

[Johan Löfberg]

The 13.8/21.2 solution was for the other other data, you flipped Gsd and Grd compared to your initial model.

[Rukhsana Ruby]

Hi Rohan,

Thank you very much for your explanation. I got it. Since the problem is non-convex, I have to find some analytical way to get closer optimal solution. i understand "bmibnb" gives the global optomal solution, however it searches exhaustively to find the global solution which is not my goal.

Regards

Ruby

[Johan Löfberg]

Why is not a global solver a good choice, if it is fast anyway? For this model, bmibnb only requires three calls to fmincon + a couple of linear programs to prove optimality.

[Rukhsana Ruby]

Hi Rohan,

Thanks for your reply. This is simplified problem with only 2 variables. Original problem has at least 10 variables. What I am thinking now is below:

My goal (with 2 variables) is to maximize:

1 + Es*Gsd + (Es*Gsr*Er*Grd)/(Es*Gsr + Er*Grd + 1)

= [(1 + Es*Gsd)(Es*Gsr + Er*Grd + 1) + Es*Gsr*Er*Grd]/[Es*Gsr + Er*Grd + 1]

which means to minimize:

[Es*Gsr + Er*Grd + 1]/[(1 + Es*Gsd)(Es*Gsr + Er*Grd + 1) + Es*Gsr*Er*Grd]

denominator and numerator both are posynomial. If I somehow can approximate denominator by monomial, I will be able to solve this problem using GP. Please see boyd tutorial (page 105):

http://www.stanford.edu/~boyd/papers/gp_tutorial.html 

Regards

Ruby

[Johan Löfberg]

To be correct: If you are able to approximate the denominator with a monomial, you will be able to solve an approximation of this problem using GP. You are thus not necessarily computing the optimal solution to the original problem, but may end up with a worse solution than you would have obtained by simply applying a local nonlinear solver on the original problem.

I would be interested to see the 10D problem.

[Rukhsana Ruby]

Hi Rohan,

Original problem is to maximize:

\mult_{i}^N[(1 + Esi*G_{si,d})(Esi*G_{si,r} + Eri*Grd + 1) + Esi*G_{si,r}*Eri*Grd]/[Esi*G_{si,r} + Eri*Grd + 1]

which means to minimize:

\mult_{i}^N[Esi*G_{si,r} + Eri*Grd + 1]/[(1 + Esi*G_{si,d})(Esi*G_{si,r} + Eri*Grd + 1) + Esi*G_{si,r}*Eri*Grd]

All G's will be given, variables are E's. If N=5, number of variables is 10. I hope, you get my latex syntax. If not, I will compile it and will send it you the pdf file.

Regards

Ruby

[Rukhsana Ruby]

Hi Rohan,

Here is the code for solving 10D problem.

function [] = gp_srcs_test()

    rx = 0;

    ry = 3;

    dx = 0;

    dy = 6;

    

    relay_dst_dist = sqrt((rx - dx)^2 + (ry-dy)^2);

    

    sx = [-2;-1;0;1;2];

    sy = [1;1;1;1;1;1];

    

    num_srcs = length(sx);

    srcs_relay_rayleigh = zeros(num_srcs, 1);

    srcs_relay_rayleigh(:, 1) = 1.4;

    relay_dst_rayleigh = 1.4;

    srcs_dst_rayleigh = zeros(num_srcs, 1);

    srcs_dst_rayleigh(:, 1) = 1.4;

    

    srcs_relay_dist = zeros(num_srcs, 1);

    srcs_dst_dist = zeros(num_srcs, 1);

    srcs_relay_snr = zeros(num_srcs, 1);

    srcs_dst_snr = zeros(num_srcs, 1);

    

    for i=1:num_srcs

        srcs_relay_dist(i, 1) = sqrt((sx(i, 1) - rx)^2 + (sy(i, 1) - ry)^2);

        srcs_dst_dist(i, 1) = sqrt((sx(i, 1) - dx)^2 + (sy(i, 1) - dy)^2);

        srcs_relay_snr(i, 1) = calculateSNR(srcs_relay_dist(i, 1)/10, srcs_relay_rayleigh(i, 1));

        srcs_dst_snr(i, 1) = calculateSNR(srcs_dst_dist(i, 1)/10, srcs_dst_rayleigh(i, 1));

    end

    relay_dst_snr = calculateSNR(relay_dst_dist/10, relay_dst_rayleigh);

    

    sdpvar Es1 Es2 Es3 Es4 Es5 Er1 Er2 Er3 Er4 Er5 t1 t2 t3 t4 t5 y1 y2 y3 y4 y5

    

    %sdpvar Es1 Es2 Er1 Er2 t1 t2 y1 y2

    

    %obj = t1*t2*t3*t4*t5;

    F = y1 + Es1*srcs_dst_snr(1)*y1 + Es1*srcs_relay_snr(1)*Er1*relay_dst_snr == y1*t1;

    F = [F, Es1*srcs_relay_snr(1) + Er1*relay_dst_snr + 1 == y1];

    

    F = [F, y2 + Es2*srcs_dst_snr(2)*y2 + Es2*srcs_relay_snr(2)*Er2*relay_dst_snr == y2*t2];

    F = [F, Es2*srcs_relay_snr(2) + Er2*relay_dst_snr + 1 == y2];

    

    F = [F, y3 + Es3*srcs_dst_snr(3)*y3 + Es3*srcs_relay_snr(3)*Er3*relay_dst_snr == y3*t3];

    F = [F, Es3*srcs_relay_snr(3) + Er3*relay_dst_snr + 1 == y3];

    

    F = [F, y4 + Es4*srcs_dst_snr(4)*y4 + Es4*srcs_relay_snr(4)*Er4*relay_dst_snr == y4*t4];

    F = [F, Es4*srcs_relay_snr(4) + Er4*relay_dst_snr + 1 == y4];

    

    F = [F, y5 + Es5*srcs_dst_snr(5)*y5 + Es5*srcs_relay_snr(5)*Er5*relay_dst_snr == y5*t5];

    F = [F, Es5*srcs_relay_snr(5) + Er5*relay_dst_snr + 1 == y5];

%     F = [F, Er1+Er2 <= 20, Es1 <= 20, Es2 <= 20, Es1+Es2+Er1+Er2 <= 50];

%     F = [F,0 <= [Es1 Es2 Er1 Er2 t1 t2 y1 y2] <= 1000];

    

    F = [F, Er1+Er2+Er3+Er4+Er5 <= 50, Es1 <= 20, Es2 <= 20, Es3 <= 20, Es4 <= 20, Es5 <= 20, Es1+Es2+Es3+Es4+Es5+Er1+Er2+Er3+Er4+Er5 <= 125];

    F = [F,0 <= [Es1 Es2 Es3 Es4 Es5 Er1 Er2 Er3 Er4 Er5 t1 t2 t3 t4 t5 y1 y2 y3 y4 y5] <= 1000];

    

    s = solvesdp(F, -t1*t2*t3*t4*t5, sdpsettings('solver','bmibnb'))

    

    double(Es1)

    double(Es2)

    double(Es3)

    double(Es4)

    double(Es5)

    double(Er1)

    double(Er2)

    double(Er3)

    double(Er4)

    double(Er5)

    

    

end

function [snr] = calculateSNR(dist, rayleighFading)

   lambda = 3.76;

   constK = 20.1; %28.1

   logNShadPower = getLogNormalShadowingPower();

   Hki = (-constK - 10 * lambda * log10(dist)) - logNShadPower + 10*log10(rayleighFading);

   snr = 10^(Hki/10);

   

end

function logNShadPower = getLogNormalShadowingPower()

    logNShadPower = 10.5;

end

function [tot_thrput] = calc_thrput(Es, srcs_dst_snr, srcs_relay_snr, Er, relay_dst_snr, sigma_2)

    num_srcs = length(Es);

    tot_thrput = 0;

for i=1:num_srcs 

        tmp1 =  Es(i, 1)*srcs_dst_snr(i, 1)/sigma_2;

        tmp2 = (Es(i)*Er*srcs_relay_snr(i)*relay_dst_snr)/(sigma_2*(sigma_2 + Es(i)*srcs_relay_snr(i) + Er*relay_dst_snr));

        tmp = log2(1 + tmp1 + tmp2);

        tot_thrput = tot_thrput + tmp;

    end

   

end

While running the code, i get the following error:

??? Error using ==> t1

Too many output arguments.

Error in ==> gp_srcs_test at 38

    F = y1 + Es1*srcs_dst_snr(1)*y1 + Es1*srcs_relay_snr(1)*Er1*relay_dst_snr ==

    y1*t1;

 

Please help me in solving this out.

Regards

Ruby

[Rukhsana Ruby]

Hi Rohan,

I think, its weekend for you. Once you are back, if could reply back to me on this problem, that is really helpful for me.

Regards

Ruby

[Johan Löfberg]

The error message means you have a function in your path called t1 (see what "which t1 -all" returns)

This problem is significantly tougher to solve. It does not close the gap in 100 iterations so you must crank up the number of iterations

s = solvesdp(F, -x1*x2*x3*x4*x5, sdpsettings('solver','bmibnb','bmibnb.maxiter',10000))

BTW, it might lead to an easier problem if you perform a logarithmic transform on the objective (leading to somewhat stronger LP relaxations since they are univariate and probably a bit tighter than the model derived for the product of several variables (it is bilinearized and cuts for every bilinear term is introduced))

s = solvesdp(F, -sum(log([t1 t2 t3 t4 t5])), sdpsettings('solver','bmibnb','bmibnb.maxiter',10000))    

Required ~5000 iterations to close the gap (optimal solution found in first iteration...)

[Rukhsana Ruby]

Hi Rohan,

Thanks for your quick reply. Alright, I will do so. Actually, the original problem has logarithm you have in the obj variable. That is even nicer for me. I will get back to you soon once i go to school.

Regards

Ruby

[Rukhsana Ruby]

Hi Johan,

Yes, I see the results. Thanks for your help. I am wondering one thing, how you have concluded that results are globally optimal. 

Regards

Ruby

[Rukhsana Ruby]

Hi Johan,

If I change the parameters, I get the results as NaN. Also, I think, for all kinds of parameters, it does not give optimal solution. For example, for the following code, I get NaN.

function [] = gp_srcs_test()

    rx = 0;

    ry = 3;

    dx = 0;

    dy = 6;

    

    relay_dst_dist = sqrt((rx - dx)^2 + (ry-dy)^2);

    

    sx = [-3;-1;0;1;3];

    sy = [2;2;2;2;2];

    

    num_srcs = length(sx);

    srcs_relay_rayleigh = zeros(num_srcs, 1);

    srcs_relay_rayleigh(:, 1) = 1.4;

    relay_dst_rayleigh = 1.4;

    srcs_dst_rayleigh = zeros(num_srcs, 1);

    srcs_dst_rayleigh(:, 1) = 1.4;

    

    srcs_relay_dist = zeros(num_srcs, 1);

    srcs_dst_dist = zeros(num_srcs, 1);

    srcs_relay_snr = zeros(num_srcs, 1);

    srcs_dst_snr = zeros(num_srcs, 1);

    

    for i=1:num_srcs

        srcs_relay_dist(i, 1) = sqrt((sx(i, 1) - rx)^2 + (sy(i, 1) - ry)^2);

        srcs_dst_dist(i, 1) = sqrt((sx(i, 1) - dx)^2 + (sy(i, 1) - dy)^2);

        srcs_relay_snr(i, 1) = calculateSNR(srcs_relay_dist(i, 1)/10, srcs_relay_rayleigh(i, 1));

        srcs_dst_snr(i, 1) = calculateSNR(srcs_dst_dist(i, 1)/10, srcs_dst_rayleigh(i, 1));

    end

    relay_dst_snr = calculateSNR(relay_dst_dist/10, relay_dst_rayleigh);

    

    sdpvar Es1 Es2 Es3 Es4 Es5 Er1 Er2 Er3 Er4 Er5 x1 x2 x3 x4 x5 y1 y2 y3 y4 y5

    

    %sdpvar Es1 Es2 Er1 Er2 t1 t2 y1 y2

    

    %obj = t1*t2*t3*t4*t5;

    F = y1 + Es1*srcs_dst_snr(1)*y1 + Es1*srcs_relay_snr(1)*Er1*relay_dst_snr == y1*x1;

    F = [F, Es1*srcs_relay_snr(1) + Er1*relay_dst_snr + 1 == y1];

    

    F = [F, y2 + Es2*srcs_dst_snr(2)*y2 + Es2*srcs_relay_snr(2)*Er2*relay_dst_snr == y2*x2];

    F = [F, Es2*srcs_relay_snr(2) + Er2*relay_dst_snr + 1 == y2];

    

    F = [F, y3 + Es3*srcs_dst_snr(3)*y3 + Es3*srcs_relay_snr(3)*Er3*relay_dst_snr == y3*x3];

    F = [F, Es3*srcs_relay_snr(3) + Er3*relay_dst_snr + 1 == y3];

    

    F = [F, y4 + Es4*srcs_dst_snr(4)*y4 + Es4*srcs_relay_snr(4)*Er4*relay_dst_snr == y4*x4];

    F = [F, Es4*srcs_relay_snr(4) + Er4*relay_dst_snr + 1 == y4];

    

    F = [F, y5 + Es5*srcs_dst_snr(5)*y5 + Es5*srcs_relay_snr(5)*Er5*relay_dst_snr == y5*x5];

    F = [F, Es5*srcs_relay_snr(5) + Er5*relay_dst_snr + 1 == y5];

%     F = [F, Er1+Er2 <= 20, Es1 <= 20, Es2 <= 20, Es1+Es2+Er1+Er2 <= 50];

%     F = [F,0 <= [Es1 Es2 Er1 Er2 t1 t2 y1 y2] <= 1000];

    

    F = [F, Er1+Er2+Er3+Er4+Er5 <= 50, Es1 <= 20, Es2 <= 20, Es3 <= 20, Es4 <= 20, Es5 <= 20, Es1+Es2+Es3+Es4+Es5+Er1+Er2+Er3+Er4+Er5 <= 125];

    F = [F,0 <= [Es1 Es2 Es3 Es4 Es5 Er1 Er2 Er3 Er4 Er5 x1 x2 x3 x4 x5 y1 y2 y3 y4 y5] <= 1000];

    

    % sum(log([x1 x2 x3 x4 x5]))

    

    s = solvesdp(F, -x1*x2*x3*x4*x5, sdpsettings('solver','bmibnb','bmibnb.maxiter',1000000000)) 

Regards

Ruby

[Johan Löfberg]

Because the solver terminated after having closed the gap completely, hence it is a globally optimal solution (since bmibnb is a global solver).

(Assuming of course that the solver/developer didn't do any mistakes, but that issue is the same for all solvers and situations. How can you trust an LP solver to return the global solution: Because the solver was designed to do so and claims it has done so.)

[Johan Löfberg]

If you get NaN, then something has crashed, since the solver cannot terminate normally within reasonable time with 1000000000 iterations done. If this reproduces, turn on the 'debug' flag in sdpsettings and see something crashes

[Rukhsana Ruby]

Hi Johan,

Thanks for your reply. Sure, I will turn the debug flag on and will let you know what exactly is crashing.

Regards

Ruby