unknown problem in solver

[yw]

Hi, Johan, I try to use Yalmip to solve an optimization problem. However, it always reports the error "Unknown problem in solver (try using 'debug'-flag in sdpsettings) (Error using full Out of memory. Type HELP MEMORY for your options.)".

I try many times  by reducing the number of constrains  and many methods to decrease memory occupation, but it didn't work. My decision variables are a matrix with 160*300 dimensions and a vector with 300 dimensions. My computer is Windows 7, 64bit. The memory information from Matlab is as follows:

 

Maximum possible array:     13454 MB (1.411e+10 bytes) *
Memory available for all arrays:     13454 MB (1.411e+10 bytes) *
Memory used by MATLAB:       807 MB (8.465e+08 bytes)
Physical Memory (RAM):      8065 MB (8.457e+09 bytes)

*  Limited by System Memory (physical + swap file) available.

 

Do you know why it happens? Can you provide me some possible solutions? Thank you so much.

[Johan Löfberg]

It happens because the solver runs out of memory. More than that is impossible to say, unless you tell us what kind of problem you are working with, which solver are you using etc

[yw]

Hi, Johan, I have initial value for these two decision variables, and try to get an solution which satisify some constraints and meanwhile are as close as possible to the inital value (every element).  My procedure is as followings:

x=sdpvar(160,300); y=sdpvar(300,1);

F=set(y(1:150,1)==y0(1:150,1))+set(x(160,271:300)==zeros(1,30));

F=set(-ones(150,1)<=sp0*x*sp1 +e-sp2*y<=ones(150,1));

F=F+ set(-ones(1,300)<=sp3*x*sp4 -y'<=ones(1,300));

F=F+set(sp5*x*sp6==be0);

F=F+set(sp7*x*sp8==be1);

F=F+set(sp9*x*sp10==be2);

g=sum(sum(((x-x0).^2)./((x0+ones(160,300)).^2)))+sum(((y-y0).^2)./((y0+ones(300,1)).^2));

sol=solvesdp(F,g);

Except x and y, others are all known matrices or vector (only e is vector).  variables started with sp are sparse matrices. 

About the solver, I try the 'quadprog', and also try not assigning any specific solver.

The total solvers in my matlab include:

 

 Searching for installed solvers       |
+++++++++++++++++++++++++++++++++++++++++++++++
|        Solver|   Version/module|      Status|
+++++++++++++++++++++++++++++++++++++++++++++++
|       LINPROG|                 |       found|
|      QUADPROG|                 |       found|
|        LMILAB|                 |       found|
|       FMINCON|        geometric|       found|
|       FMINCON|         standard|       found|
|    FMINSEARCH|                 |       found|
|           BNB|                 |       found|
|      BINTPROG|                 |       found|
|        CUTSDP|                 |       found|
|        BMIBNB|                 |       found|
|         KKTQP|                 |       found|
|          NONE|                 |       found|
|     LSQNONNEG|                 |       found|
|        LSQLIN|                 |       found|

[Johan Löfberg]

As it says on the front page, please supply reproducible code. The code above does not run (create suitable code to generate (possibly) random x0 etc)

BTW, note that

sp3*x*sp4 -y'

is different from

sp3*x*sp4-y'

I.e, there is a dangerous space in there. Compare with

[1 -1]
[1-1]

[yw]

Hi, Johan, thank you very much for your reply and correction. Sorry for only supplying parts of the code before. Now, I added the reproduible code by using random generator.  To avoid the conflicts probably caused by randoms, I decrease both constraints and dimensitions of decision variables.  The full code is as follows:

 

x0=500*rand(154, 290);  

y0=1000*rand(272,1);

e=0.1*y0(137:272,1);

be0=800*rand(35, 38); be0(35,35:38)=zeros(1,4);

 

q0(1:272,1:34)=repmat(eye(34),8,1); q0(273:288, 35:36)=repmat(eye(2),8,1);

q0(289:290,37:38)=eye(2);

q1(1:17,1:136)=repmat(eye(17),1,8);  q1(18:35,137:154)=eye(18);

q2=zeros(136,272); q2(1:136, 137:272)=eye(136);

q3=ones(290,1); q4=zeros(136,154); q4(1:136,1:136)=eye(136);

q5=ones(1,154); q6=zeros(290,272); q6(1:272,1:272)=eye(272);

sp0=sparse(q0); sp1=sparse(q1); sp2=sparse(q2);

sp3=sparse(q3); sp4=sparse(q4); sp5=sparse(q5); sp6=sparse(q6);

clearvars q0 q1 q2 q3 q4 q5 q6;

 

x=sdpvar(154,290); y=sdpvar(272,1);

F=set(y(1:136,1)==y0(1:136,1))+set(x(154,273:290)==zeros(1,18));

F=F+set(sp1*x*sp0==be0); 

F=set(-ones(136,1)<=sp4*x*sp3+e-sp2*y<=ones(136,1));

F=F+set(-ones(1,272)<=sp5*x*sp6-y'<=ones(1,272));

g=sum(sum(((x-x0)./(x0+ones(154,290))).^2))+sum(((y-y0)./(y0+ones(272,1))).^2);

[yw]

Besides, be careful, because it will take a long time before the solvers give any result.

[Johan Löfberg]

Several issues with your code (not all related to the problem though)

1. Don't use SET

http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Commands.Set

2. Why are you working with dense matrices, and then convert to sparse. Create sparse to begin with

q0(1:272,1:34)=repmat(speye(34),8,1); q0(273:288, 35:36)=repmat(speye(2),8,1);

q0(289:290,37:38)=speye(2);
q2=spalloc(136,272,0); q2(1:136, 137:272)=speye(136);
etc...

3. It looks to me as if you have variables naturally divided into different logical parts, but you have concatenated them into two large things x and y, and then you model the problem by extracting suitable parts from x and y. It looks like the model can be described much more naturally

4. YALMIP vectorizes matrix/scalar things like MATLAB. Hence this

x(154,273:290)==zeros(1,18)
-ones(1,272)<=sp5*x*sp6-y'<=ones(1,272)

can just as well be written as

x(154,273:290)==0
-1<=sp5*x*sp6-y'<=1

5. x has no structure (it is a fully parameterized matrix), so why waste decision variables on things you know are 0?

x=sdpvar(154,290);
x(154,273:290)==zeros(1,18)

Simply place zeros in x instead

x=sdpvar(154,290);
x(154,273:290)=0

6. Your actual problem comes from the fact that quadprog uses a dense method here. You have to tweak the options to force it to work with a sparse model (although I would recommend you to install a better solver such as gurobi)

solvesdp(F,g,sdpsettings('solver','quadprog','quadprog.largescale','on','quadprog.algorithm','interior-point-convex'))

7. However, the main problem with the code, once using the correct solver, is due to a performance issue in YALMIP. When creating the symbolic representation of your quadratic objective involving ~45000 monomials, it takes a lot of overhead. By exploiting the fact that YALMIP is much much faster in working with pure quadratic forms, I would write the last stage as

z1 = sdpvar(154,290);

z2 = sdpvar(272,1);

g = sum(sum(z1.^2)) + sum(z2.^2);

F = [F, z1 == (x-x0)./(x0+1)];

F = [F, z2 == (y-y0)./(y0+1)];

solvesdp(F,g,sdpsettings('solver','quadprog','quadprog.largescale','on','quadprog.algorithm','interior-point-convex'))

The whole code, involving the solution, now runs in a couple of seconds.

[yw]

Thank you so much. Yes, it works. It is my first time to use Yalmip. I have to say that it is really intelligent.

 Two more questions: Besides gurobi, which solvers are also good choices for solving optimization with large scale of variables? 

 If I change the target function to g=sum(sum(abs(z1))+sum(abs(z2); which solvers are better?

[Johan Löfberg]

cplex, gurobi and mosek are all strong lp/qp solvers with free licenses for academia

[Edi]

Hello Johan,

I have the same problem as yw. I am new user of yalmip and it seems that I have got an error saying  "Unknown problem in solver (try using 'debug'-flag in sdpsettings) (Error using full Out of memory. Type HELP MEMORY for your options.)".

The code is as following which is a simple optimization problem.

Your suggestions are very valuable to me.

  

N=10;

hk_rand=abs(randn(1,N));

hk=hk_rand;

No=0.1;

A=0.11;

Pt=10;

sigma=0.3;

   L=13.5                                         %length of bridge

    sk=[];

 

     for j=1:N

    

       

       zk=(j-1)*L/N;

%----------------------Adjust the SNR for each sensor----------------------------------

sk1=0.15*cos(pi*zk/L);                                %adjusting SNR 

     sk=[sk sk1];

     end

     

     

     

SNR_i=(sk.^2)./sigma.^2;

pk = sdpvar(1,N);                         %transmit power

     qk=(1/2)*log((1+(pk.*hk/No)));

     

    sigmaq1=(A^2)./(3*(2^(2*qk)));    %quantization variance

 

   h3=((2*N*sigma.^4.*(1+2*SNR_i)+sigmaq1).*(2*N*sigma.^2)-(2*N.*sigma.^4+sigmaq1).*(2*N*sigma.^2.*(1+2*SNR_i)))./(8*N*sigma.^4.*SNR_i); %with quantization

   a3=-((2*N*sigma.^4+sigmaq1)-2*N*sigma.^4.*(1+2*SNR_i)-sigmaq1)./((2*N.*sigma.^4.*(1+2*SNR_i)+sigmaq1).*(2*N+sigmaq1));

   

    E_T_H0=N^2.*sigma.^4+2.*N.*sigma.^4+sigmaq1-2.*h3.*N.*sigma.^2+h3.^2;

    E_T_H1=N^2.*sk.^4+2.*N^2.*sk.^2.*sigma.^2+4.*N.*sigma.^2.*sk.^2+N^2.*sigma.^4+2.*N.*sigma.^4+sigmaq1-2.*h3.*(N.*sk.^2+N.*sigma.^2)+h3.^2;              %E{T¦H1}

    eta_i=a3.*(E_T_H1-E_T_H0);

Var_T_H1=4*(N^2.*sk.^4+2*N^2.*sk.^2.*sigma.^2+N^2.*sigma.^4).*(2*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2)+2.*(2.*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2).^2 ...

 +4.*h3.^2.*(2.*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2)-8.*h3.*(N.*sk.^2+N.*sigma.^2).*(2.*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2);

Var_T_H11=(a3.^2).*Var_T_H1;

  

  

    g=(a3.*eta_i);

    f=((a3.^2).*Var_T_H1);

   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  

    Object_ive=sum(g)/sum(f);

    Constraints = [sum(pk)<= Pt, pk>=0];

    

    %solvesdp( Constraints,Object_ive,sdpsettings('solver','quadprog','quadprog.largescale','on','quadprog.algorithm','interior-point-convex'))

    solvesdp( Constraints,Object_ive,sdpsettings('debug','flag'))

Thanks

Edi

[Johan Löfberg]

The fact that you have been trying to use quadprog to solve the problem indicates that you have failed completely to understand just how hard this problem is. What makes you think it is solvable using a QP solver? You have loads of nonconvex operations involving divisions, powers, multiplications etc on your decision variables. In the current form, your model is probably not solvable by any nonlinear solver.

Assigning an initial guess, and the solving with fmincon makes fmincon terminate directly and simply return the initial guess

pk = sdpvar(1,N);                        
assign(pk,.1);

...

solvesdp(Constraints,Object_ive,sdpsettings('solver','fmincon','usex0',1))

[Edi]

Dear Johan,

Thank you for your very fast reply. Sorry, the "quadprog" is commented and not used. I didn't specify the solver as it is shown "solvesdp( Constraints,Object_ive,sdpsettings('debug','flag'))."

Actually yes, I have a lot of operations that does not preserve the convexity but the function f and g can be shown to be convex with respect to pk.

The problem than reduces to a convex-convex fractional problem.

I have been trying to solve the problem using cvx combined with Dinkelbach's algorithm but I failed to do so.

I am not sure  if this can be solved by yumli though.

Sorry for the inconvenience caused.

Your valuable comments are welcome.

Regards

Edmond

[Mark L. Stone]

Do you know a feasible starting value?  If so, you can give fmincon (or perhaps knitro, if you have it) a fighting chance.

With regard to CVX, did you ever get to the point where CVX accepted your model and invoked a solver, but then the solver failed, or did you not succeed in getting CVX to accept your formulation?  Of course, for CVX, it is necessary, but not sufficient, that apart form integer constraints, the model must be convex; it must also be expressible as a Disciplined Convex Program using CVX's ruleset and available functions, in which case, apart from the integer constraints, it will necessarily be convex.

[Johan Löfberg]

CVX would never get started on this model. It is highly undisciplined...

The code I gave above sets an initial guess and calls fmincon, but the problem seems very ill-posed so essentially all random initial guess yield objective values around 0.36 with first-order optimality within the tolerance of fmincon, so fmincon terminates directly

[Johan Löfberg]

knitro manages better to move away from the initial guess

p0 = rand(1,10);p0 = 9.99*p0/sum(p0)
assign(pk,p0);
double(Object_ive)
solvesdp(Constraints,Object_ive,sdpsettings('solver','knitro','usex0',1))

ans =

   0.359780860450551

>> double(Object_ive)

ans =

   0.359331956065639

 [double(pk)' double(p0)']

ans =

   1.233239664201927   1.268054055747984
   1.530276044100641   1.532929627773902
   0.000048320845871   0.060997533699532
   1.570859697381596   1.571505652792795
   1.264442631845528   1.266252936698882
   1.376820312578258   1.579469386563020
   0.903066155919263   0.904380888201387
   0.358447969130552   0.416595368954150
   1.148274914472070   1.211664745351628
   0.000017362436041   0.178149804216722

ipopt stalls, and fmincon exits immediately as I said earlier

[Edi]

Dear Johan

Thank you for your valuable comments. I think I managed to work somehow. However I have a problem in that: I want to max f(p)/g(p); If that was convex that would be easily converted. When I minimize f(p)/g(p) it seems that it works as expected.

I am not sure how to convert a maximization problem into a minimization problem that Yalmip understands.

Here is the code:

N=10;

 hk_rand=abs(randn(1,N));

 hk=hk_rand;

%hk=[ 0.4361    0.4017    0.8862    0.3827    1.2176    0.5227    1.6672    0.2297    0.5548    1.1178];

No=0.1;

A=0.5;

Pt=10;

sigma=0.3;

   L=13.5                                         %length of bridge

    sk=[];

%  

%      for j=1:N

%     

%        

%        zk=(j-1)*L/N;

% %----------------------Adjust the SNR for each sensor----------------------------------

% sk1=0.15*cos(pi*zk/L);                                %adjusting SNR 

%      sk=[sk sk1];

%      end

     %sk=randn(1,N);

     sk=1;

SNR_i=(sk.^2)./sigma.^2;

pk = sdpvar(1,N);                       

p0 = rand(1,10);p0 = 9.99*p0/sum(p0)

assign(pk,p0);

%assign(pk,.5);

     qk=(1/2)*log2((1+(pk.*hk/No)));

     

    sigmaq1=(A^2)./(3*(2^(2*qk)));

 

   h3=((2*N*sigma.^4.*(1+2*SNR_i)+sigmaq1).*(2*N*sigma.^2)-(2*N.*sigma.^4+sigmaq1).*(2*N*sigma.^2.*(1+2*SNR_i)))./(8*N*sigma.^4.*SNR_i); %with quantization

   %a3=-((2*N*sigma.^4+sigmaq1)-2*N*sigma.^4.*(1+2*SNR_i)-sigmaq1)./((2*N.*sigma.^4.*(1+2*SNR_i)+sigmaq1).*(2*N+sigmaq1));

   a3=1;

   

    E_T_H0=N^2.*sigma.^4+2.*N.*sigma.^4+sigmaq1-2.*h3.*N.*sigma.^2+h3.^2;

    E_T_H1=N^2.*sk.^4+2.*N^2.*sk.^2.*sigma.^2+4.*N.*sigma.^2.*sk.^2+N^2.*sigma.^4+2.*N.*sigma.^4+sigmaq1-2.*h3.*(N.*sk.^2+N.*sigma.^2)+h3.^2;              %E{T¦H1}

eta_i=(E_T_H1-E_T_H0).^2;

Var_T_H1=4*(N^2.*sk.^4+2*N^2.*sk.^2.*sigma.^2+N^2.*sigma.^4).*(2*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2)+2.*(2.*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2).^2 ...

 +4.*h3.^2.*(2.*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2)-8.*h3.*(N.*sk.^2+N.*sigma.^2).*(2.*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2);

  

  

    f=((a3.^2).*eta_i);

    g=((a3.^2).*Var_T_H1);

   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  

    Object_ive=(sum(f)/sum(g));

    Constraints = [sum(pk)<= Pt, pk>=0]

    

   double(Object_ive)

    

    

    

%sol4=solvesdp(Constraints,Object_ive,sdpsettings('solver','fmincon','usex0',1))

ops = sdpsettings('solver','bmibnb')

sol = solvesdp(Constraints,Object_ive,ops)

ans =

    0.5895    0.4044

    0.5826    0.5130

    0.7912    0.8214

    1.3772    1.1766

    0.6457    1.4474

    2.0255    0.2997

    0.6502    0.9557

    1.5144    1.8107

    1.2317    2.0166

    0.5921    0.5444

I agree as well that it is impossible to be solved with cvx this problem. There are a lot of issues....

I would be grateful if any help regrading the max to min conversion of optimization

Thanks

[Johan Löfberg]

minimize -f/g, or minimize g/f

[Johan Löfberg]

This...

    qk=(1/2)*log2((1+(pk.*hk/No)));
    sigmaq1=(A^2)./(3*(2^(2*qk)));

Hint: What is 2^(2*0.5*log2(1+x))

[Edi]

sigmaq1 is the quantization variance, where qk is the number of quantized bits.

[Johan Löfberg]

What I was hinting at was the fact that the 2^ operator and the log2 operator cancel each other, so your model can be simplified significantly (this is not eliminated by YALMIP, so YALMIP has to do automatic differentiation and all that stuff with expressions involving log2 and powers and work with a very complex expression which really is linear)

[Edi]

Yeah that is correct.... Now it is simplier and  faster. But still I tried even before to min g/f but the solution it gives is pk=0 which seems to be a failure. I know in reality what a success should look like. I am not sure if that is the problem with the solver.

[Johan Löfberg]

Which solver?

[Edi]

The "bmibnb" which is called by Yalmi.

What do you suggest?

Thanks

[Johan Löfberg]

Does it really run? What is displayed? It crashes here in call to lower bound solver

[Johan Löfberg]

Seem to be a problem when I use QP-capable solvers to solve the LP relaxation. When using an LP-solver, the problem is solved to global optimality in no time

The optimal solution is pk=0. You say it shouldn't be, but then the model is incorrectly specified. Simply search randomly and you will see that you never beat pk=0

best = inf;
for i = 1:1000000
pk = rand(1,10);pk = 9.99*p0/sum(p0);

    qk=(1/2)*log2((1+(pk.*hk/No)));        
    sigmaq1=(A^2)./(3*(1+(pk.*hk/No)));

 
   h3=((2*N*sigma.^4.*(1+2*SNR_i)+sigmaq1).*(2*N*sigma.^2)-(2*N.*sigma.^4+sigmaq1).*(2*N*sigma.^2.*(1+2*SNR_i)))./(8*N*sigma.^4.*SNR_i); %with quantization
   %a3=-((2*N*sigma.^4+sigmaq1)-2*N*sigma.^4.*(1+2*SNR_i)-sigmaq1)./((2*N.*sigma.^4.*(1+2*SNR_i)+sigmaq1).*(2*N+sigmaq1));
   a3=1;
   
    E_T_H0=N^2.*sigma.^4+2.*N.*sigma.^4+sigmaq1-2.*h3.*N.*sigma.^2+h3.^2;
    E_T_H1=N^2.*sk.^4+2.*N^2.*sk.^2.*sigma.^2+4.*N.*sigma.^2.*sk.^2+N^2.*sigma.^4+2.*N.*sigma.^4+sigmaq1-2.*h3.*(N.*sk.^2+N.*sigma.^2)+h3.^2;              %E{T¦H1}

eta_i=(E_T_H1-E_T_H0).^2;

Var_T_H1=4*(N^2.*sk.^4+2*N^2.*sk.^2.*sigma.^2+N^2.*sigma.^4).*(2*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2)+2.*(2.*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2).^2 ...
 +4.*h3.^2.*(2.*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2)-8.*h3.*(N.*sk.^2+N.*sigma.^2).*(2.*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2);

    f=((a3.^2).*eta_i);
    g=((a3.^2).*Var_T_H1);

   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 
    Object_ive=(sum(g)/sum(f));
   
    if Object_ive < best
        best = Object_ive
        bestp = pk;
    end
    best
end

[Edi]

when I minimize f/g it works and it assigned the values . it runs with no error at all. Here is the displayed:

p0 =

    1.7811    1.4237    0.1572    1.4204    0.8276    1.3390    0.9394    0.2220    1.3443    0.5351

Sigmonial matrix variable 1x10 (full, real, 10 variables, values in range [0.034798,0.080769])

++++++++++++++++++++++++++++++++++++++++++++++++++++++++

|   ID|      Constraint|                           Type|

++++++++++++++++++++++++++++++++++++++++++++++++++++++++

|   #1|   Numeric value|    Element-wise inequality 1x1|

|   #2|   Numeric value|   Element-wise inequality 10x1|

++++++++++++++++++++++++++++++++++++++++++++++++++++++++

ans =

    8.5195

* Starting YALMIP global branch & bound.

* Branch-variables : 21

* Upper solver     : fmincon

* Lower solver     : LINPROG

* LP solver        : LINPROG

 Node       Upper      Gap(%)       Lower    Open

    1 :   8.331E+000     0.00     8.331E+000   0  Infeasible  

* Finished.  Cost: 8.3308 Gap: 0

* Timing: 16% spent in upper solver

*         0% spent in lower solver

*         61% spent in LP-based domain reduction

sol3 = 

    yalmiptime: 0.7260

    solvertime: 3.1940

          info: 'Successfully solved (BMIBNB)'

       problem: 0

this is the result:

 [double(pk)' double(p0)']

ans =

    1.0008    1.7811

    0.5354    1.4237

    0.9038    0.1572

    0.1037    1.4204

    0.7202    0.8276

    3.5041    1.3390

    1.1655    0.9394

    0.1205    0.2220

    1.5976    1.3443

    0.2661    0.5351

[Edi]

I think I got what was the problem. I will implement it tomorrow and I will contact you again. It seems that it works now....

Thank you for your help. And congratulations on developing this tool. I bet it is amazing. Thank you for your valuable contribution in research.

Have a good night.

[Edi]

Dear Johan

I am trying now to improve my skills using yalmip. I have the following code now and I get the error:

Error using  * 

Inner matrix dimensions must agree.

Error in bmibnb>diagonalize_quadratic_program (line 480)

    pnew.x0 = V*p.x0(1:2);

Error in bmibnb (line 116)

p = diagonalize_quadratic_program(p);

Error in solvesdp (line 355)

    eval(['output = ' solver.call '(interfacedata);']);

Error in yalmi_edmond_test (line 138)

  sol = solvesdp(Constraints,Object_ive,ops)

%*****************====================******************************

The code is:

N=20;

 hk_rand=2*abs(randn(1,N));

 hk=hk_rand;

%hk=[ 0.4361    0.4017    0.8862    0.3827    1.2176    0.5227    1.6672    0.2297    0.5548    1.1178];

No=0.1;

A=1;

Pt=10;

sigma=0.36;

   L=13.5                                         %length of bridge

    sk=[];

 

     for j=1:N

    

       

       zk=(j-1)*L/N;

%----------------------Adjust the SNR for each sensor----------------------------------

sk1=0.1*cos(pi*zk/L);                                %adjusting SNR 

     sk=[sk sk1];

     end

     %sk=randn(1,N);

     %sk=1;

SNR_i=(sk.^2)./sigma.^2;

pk = sdpvar(1,N);                       

p0 = rand(1,N);p0 = 9.99*p0/sum(p0)

assign(pk,p0);

%assign(pk,.5);

     qk=(1/2)*log2((1+(pk.*hk/No)));

     

    %sigmaq1=(A^2)./(3*(2^(2*qk)));

    

       sigmaq1=(A^2)./(3*(1+pk.*hk/No))

 

   h3=((2*N*sigma.^4.*(1+2*SNR_i)+sigmaq1).*(2*N*sigma.^2)-(2*N.*sigma.^4+sigmaq1).*(2*N*sigma.^2.*(1+2*SNR_i)))./(8*N*sigma.^4.*SNR_i); %with quantization

   a3=-((2*N*sigma.^4+sigmaq1)-2*N*sigma.^4.*(1+2*SNR_i)-sigmaq1)./((2*N.*sigma.^4.*(1+2*SNR_i)+sigmaq1).*(2*N+sigmaq1));

   %a3=a3./(norm(a3));

   %a3=1;

   

    E_T_H0=N^2.*sigma.^4+2.*N.*sigma.^4+sigmaq1-2.*h3.*N.*sigma.^2+h3.^2;

    E_T_H1=N^2.*sk.^4+2.*N^2.*sk.^2.*sigma.^2+4.*N.*sigma.^2.*sk.^2+N^2.*sigma.^4+2.*N.*sigma.^4+sigmaq1-2.*h3.*(N.*sk.^2+N.*sigma.^2)+h3.^2;              %E{T¦H1}

eta_i=(E_T_H1-E_T_H0).^2;

Var_T_H0=4*((N*sigma.^2).^2).*(2*N.*sigma.^4+sigmaq1)+2.*(2.*N.*sigma.^4+sigmaq1).^2+4.*h3.^2.*(2.*N.*sigma.^4+sigmaq1)-8.*h3.*N.*sigma.^2.*(2.*N.*sigma.^4+sigmaq1);

Var_T_H1=4*(N^2.*sk.^4+2*N^2.*sk.^2.*sigma.^2+N^2.*sigma.^4).*(2*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2)+2.*(2.*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2).^2 ...

 +4.*h3.^2.*(2.*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2)-8.*h3.*(N.*sk.^2+N.*sigma.^2).*(2.*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2);

  

  

    f=((a3.^2).*eta_i);

    g=((a3.^2).*Var_T_H1);

%===================================================================

kkk_1=(sum(double(E_T_H1)).^2)./sum(double(Var_T_H1));             %shape parameter Under H1

ooo_1=sum(double(Var_T_H1))./sum(double(E_T_H1));              %scale parameter under H1

kkk_0=(sum(double(E_T_H0)).^2)./sum(double(Var_T_H0));        %shape parameter Under H0

ooo_0=sum(double(Var_T_H0))./sum(double(E_T_H0));        %scale parameter under H0

Pfa=0.1;

Thresh= gaminv(Pfa,kkk_0,ooo_0);

Pd=1-gamcdf(Thresh,kkk_1,ooo_1);

Object_ive=1/Pd

    Constraints = [sum(pk)== Pt, pk>=0]

    

   double(Object_ive)

    

    

    

%sol4=solvesdp(Constraints,Object_ive,sdpsettings('solver','fmincon','usex0',1))

  ops = sdpsettings('solver','bmibnb','bmibnb.maxiter',50,'usex0',1,'debug','flag')

  sol = solvesdp(Constraints,Object_ive,ops)

I cant see what is causing the problem. Your help is very appreciated.

Kind Regards

Edmond

[Johan Löfberg]

It is due to a bug in bmibnb

Simply turn off usex0 and it goes through

[Johan Löfberg]

The reason is partly due to a trivial bug from some remaining debug experiment (remove the (1:2) on line 480 in the crashing file). However, it still crashes elsewhere. The reason is that BMIBNB gets confused (i.e., another bug) as the problem is too simple. It is a convex QP. No reason to use BMIBNB

 sol = solvesdp(Constraints,Object_ive)
Optimize a model with 21 rows, 20 columns and 40 nonzeros
Presolve removed 21 rows and 20 columns
Presolve time: 0.00s
Presolve: All rows and columns removed
Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    0.0000000e+00   0.000000e+00   0.000000e+00      0s

Solved in 0 iterations and 0.00 seconds
Optimal objective  0.000000000e+00

sol =

    yalmiptime: 2.000000000000028e-01
    solvertime: 6.000000000000227e-03
          info: 'Successfully solved (GUROBI-GUROBI)'
       problem: 0

[Edi]

Hello Johan

I am trying to solve the similar problem but now I slightly changed my objective function and I get the following error:  solvertime: 0

       problem: 14

          info: 'Convexity check failed (Model not available in objective at level 1)'

    yalmiptime: 0.2150

From previous post you have stated that BMIBNB gets confused. Can you please recommend any other solver to solve this problem? I tried the kktqp but same error.

N=20;

 hk_rand=2*abs(randn(1,N));

 hk=hk_rand;

%hk=[ 0.4361    0.4017    0.8862    0.3827    1.2176    0.5227    1.6672    0.2297    0.5548    1.1178];

No=0.1;

A=1;

Pt=20;

sigma=0.36;

   L=13.5                                         %length of bridge

    sk=[];

%  

     for j=1:N

    

       

       zk=(j-1)*L/N;

%----------------------Adjust the SNR for each sensor----------------------------------

sk1=0.1*cos(pi*zk/L);                                %adjusting SNR 

     sk=[sk sk1];

     end

     %sk=randn(1,N);

     %sk=1;

SNR_i=(sk.^2)./sigma.^2;

pk = sdpvar(1,N);                       

p0 = rand(1,N);p0 = 9.99*p0/sum(p0)

assign(pk,p0);

%assign(pk,.5);

     qk=(1/2)*log2((1+(pk.*hk/No)));

     

    %sigmaq1=(A^2)./(3*(2^(2*qk)));

    

       sigmaq1=(A^2)./(3*(1+pk.*hk/No))

 

   h3=((2*N*sigma.^4.*(1+2*SNR_i)+sigmaq1).*(2*N*sigma.^2)-(2*N.*sigma.^4+sigmaq1).*(2*N*sigma.^2.*(1+2*SNR_i)))./(8*N*sigma.^4.*SNR_i); %with quantization

   a3=-((2*N*sigma.^4+sigmaq1)-2*N*sigma.^4.*(1+2*SNR_i)-sigmaq1)./((2*N.*sigma.^4.*(1+2*SNR_i)+sigmaq1).*(2*N+sigmaq1));

   %a3=a3./(norm(a3));

   %a3=1;

   

    E_T_H0=N^2.*sigma.^4+2.*N.*sigma.^4+sigmaq1-2.*h3.*N.*sigma.^2+h3.^2;

    E_T_H1=N^2.*sk.^4+2.*N^2.*sk.^2.*sigma.^2+4.*N.*sigma.^2.*sk.^2+N^2.*sigma.^4+2.*N.*sigma.^4+sigmaq1-2.*h3.*(N.*sk.^2+N.*sigma.^2)+h3.^2;              %E{T¦H1}

eta_i=(E_T_H1-E_T_H0);

Var_T_H0=4*((N*sigma.^2).^2).*(2*N.*sigma.^4+sigmaq1)+2.*(2.*N.*sigma.^4+sigmaq1).^2+4.*h3.^2.*(2.*N.*sigma.^4+sigmaq1)-8.*h3.*N.*sigma.^2.*(2.*N.*sigma.^4+sigmaq1);

Var_T_H1=4*(N^2.*sk.^4+2*N^2.*sk.^2.*sigma.^2+N^2.*sigma.^4).*(2*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2)+2.*(2.*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2).^2 ...

 +4.*h3.^2.*(2.*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2)-8.*h3.*(N.*sk.^2+N.*sigma.^2).*(2.*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2);

  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    

 Pfa=0.1;

    

 Object_ive=(qfuncinv(Pfa)*sqrt(sum(((a3.^2).*(Var_T_H0))))-sum((a3.*eta_i)))/(sqrt(sum((a3.^2).*(Var_T_H1))));

   

  

    Constraints = [sum(pk)== Pt, pk>=0]

    

   double(Object_ive)

    

    

    

%sol4=solvesdp(Constraints,Object_ive,sdpsettings('solver','fmincon','usex0',1))

% 

%   ops = sdpsettings('solver','bmibnb','bmibnb.maxiter',50,'usex0',1)

%   sol = solvesdp(Constraints,Object_ive,ops)

ops3= sdpsettings('solver','kktqp');

sol3 = solvesdp(Constraints,Object_ive,ops3)

% ops4= sdpsettings('solver','qpip');

% solvesdp(Constraints,Object_ive,ops4);

Power_test=double(pk);

figure

subplot(4,1,1)

bar(hk),colormap(cool);

 xlabel('sensors'); ylabel('hk^2');

grid on

colormap([1 1 0])

subplot(4,1,2)

bar(Power_test)

xlabel('sensors'); ylabel('pk*');

grid on

subplot(4,1,3)

bar(a3)

xlabel('sensors'); ylabel('combiner (a3)');

grid on

subplot(4,1,4)

bar(qk)

xlabel('sensors'); ylabel('quantization bits');

%=====================================================================

%============================================================

 %ops1 = sdpsettings('solver','bmibnb','bmibnb.maxiter',20);

 %ops1 = sdpsettings(ops1,'bmibnb.uppersolver','fmincon');

%  ops2 = sdpsettings('solver','moment','moment.order',2)

% % 

  %sol_test1 = solvesdp(Constraints,Object_ive,ops1)

% sol_test2 = solvesdp(Constraints,Object_ive,ops2)

% 

% 

% 

% 

% 

 %Power_test=double(pk);

% 

% figure(1)

% subplot(3,1,1)

% bar(hk),colormap(cool);

%  xlabel('sensors'); ylabel('hk^2');

% grid on

% colormap([1 1 0])

% subplot(3,1,2)

% bar(Power_test)

% xlabel('sensors'); ylabel('pk*');

% grid on

% subplot(3,1,3)

% bar(a3)

% xlabel('sensors'); ylabel('combiner (a3)');

% % 

% 

%==========================New Code=============================

% N=10;

%  hk_rand=abs(randn(1,N));

%  hk=hk_rand;

% %hk=[ 0.4361    0.4017    0.8862    0.3827    1.2176    0.5227    1.6672    0.2297    0.5548    1.1178];

% No=0.1;

% A=0.5;

% Pt=10;

% sigma=0.3;

% 

%    L=13.5                                         %length of bridge

% 

%     sk=[];

% %  

% %      for j=1:N

% %     

% %        

% %        zk=(j-1)*L/N;

% % %----------------------Adjust the SNR for each sensor----------------------------------

% % sk1=0.15*cos(pi*zk/L);                                %adjusting SNR 

% %      sk=[sk sk1];

% %      end

%      %sk=randn(1,N);

%      sk=1;

% SNR_i=(sk.^2)./sigma.^2;

% 

% pk = sdpvar(1,N);                       

% p0 = rand(1,10);p0 = 9.99*p0/sum(p0)

% assign(pk,0.1);

% for i = 1:100

% pk = rand(1,10);pk = 9.99*p0/sum(p0);

% 

%     qk=(1/2)*log2((1+(pk.*hk/No)));        

%     sigmaq1=(A^2)./(3*(1+(pk.*hk/No)));

% 

%  

%    h3=((2*N*sigma.^4.*(1+2*SNR_i)+sigmaq1).*(2*N*sigma.^2)-(2*N.*sigma.^4+sigmaq1).*(2*N*sigma.^2.*(1+2*SNR_i)))./(8*N*sigma.^4.*SNR_i); %with quantization

%    %a3=-((2*N*sigma.^4+sigmaq1)-2*N*sigma.^4.*(1+2*SNR_i)-sigmaq1)./((2*N.*sigma.^4.*(1+2*SNR_i)+sigmaq1).*(2*N+sigmaq1));

%    a3=1;

%    

%     E_T_H0=N^2.*sigma.^4+2.*N.*sigma.^4+sigmaq1-2.*h3.*N.*sigma.^2+h3.^2;

%     E_T_H1=N^2.*sk.^4+2.*N^2.*sk.^2.*sigma.^2+4.*N.*sigma.^2.*sk.^2+N^2.*sigma.^4+2.*N.*sigma.^4+sigmaq1-2.*h3.*(N.*sk.^2+N.*sigma.^2)+h3.^2;              %E{T¦H1}

% 

% eta_i=(E_T_H1-E_T_H0).^2;

% 

% Var_T_H1=4*(N^2.*sk.^4+2*N^2.*sk.^2.*sigma.^2+N^2.*sigma.^4).*(2*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2)+2.*(2.*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2).^2 ...

%  +4.*h3.^2.*(2.*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2)-8.*h3.*(N.*sk.^2+N.*sigma.^2).*(2.*N.*sigma.^4+4.*N.*sigma.^2.*sk.^2);

% 

%     f=((a3.^2).*eta_i);

%     g=((a3.^2).*Var_T_H1);

% 

%    

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  

%     Object_ive=(sum(f)/sum(g));

%     

%     if Object_ive < best

%         best = Object_ive

%         bestp = pk;

%     end

%     best

% end

% 

% 

[Johan Löfberg]

Define "slightly changed"

[Johan Löfberg]

OK, now I see in your code.
http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Extra.SQRT

[Edi]

Hi Johan,

Can you please help me how to define an integer constraint using yalmip?

I am trying to understand how yalmip does the reprocessing and bound-propagation for improving the bound. What are the criterion  used for branching (i.e., which rectangle to branch, where to branch it)?

Secondly, How the convex relaxation is done?

I want to know more about this interesting solver "bmibnp".

Thanks a lot

Edmond

[Johan Löfberg]

You define integrality constraints by defining variables using intvar or using the command integer.

bmibnb implements a standard spatial branch&bound strategy
http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Solvers.BMIBNBTheory

[Edi]

Dear Johan

Is the command ''prod'' available in yalmip? I am trying to run this part of the code:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Lk=0.01;                  

Emax=1*ones(1,M);                  

t_N=2;       

             

pk = sdpvar(M,t_N); 

assign(pk,.5);

for i=1:M

    for j=1:t_N

        

        y(i,j)=((((1+(pk(i,j)*(hk(i)^2)/SNR_i(i)))^-Lk)));

        

   

    end

end

y_true=prod((y'))                                  % (1)

for ii=1:M

    for jj=1:t_N

        

   

    Object_ive(ii)=(((3*N^2)*(sigma(ii)^4)*(SNR_i(ii)^2))./((6*N*(sigma(ii)^4).*(1+2*SNR_i(ii))+((A^2)/y_true(ii))))); 

    

       end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

It gives an error in the Object_ive function. However one I replace (1) (i.e., y_true=prod((y')) with y_true=prod(double(y')) it works but I am not sure if the optimization problem is solved successfully. (I have a doubt)

Your help and comment is appreciated.

Kind Regards

Edmond

[Johan Löfberg]

Yes, product is supported, however there appears to be a bug since it always returns a scalar despite the argument being a matrix

BTW, how do you intend to solve the problem? You are creating a pretty nasty problem

[Johan Löfberg]

as a work-around, you would have to do

yprod=[];
for i=1:size(y,2)
 yprod=[yprod prod(y(:,i))];
end

[Edi]

I got it, thank you very much for your help.

Regards

[Edi]

Hi Johan

I am trying to solve a convex problem using yalmip. If I use the "bmibnb" global solver it returns to me the solution with a reltively very small gap. However, I was trying instead of using branch and bound method, trying to find a possible global optimum solver.

When I use"fmincon" solver it gives an error:

++++++++++++++++++++++++++++++++++++++++++++++++++++++++

|   ID|      Constraint|                           Type|

++++++++++++++++++++++++++++++++++++++++++++++++++++++++

|   #1|   Numeric value|   Element-wise inequality 40x1|

|   #2|   Numeric value|   Element-wise inequality 40x1|

|   #3|   Numeric value|   Element-wise inequality 40x1|

++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Your initial point x0 is not strictly between bounds lb and ub; FMINCON

 shifted x0 to strictly satisfy the bounds.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Can you please suggest what might be the problem?

Thanks

Edi

[Johan Löfberg]

Some of the algorithms in fmincon require a feasible initial point. Hence, you either have to assign a feasible solution to all variables (using assign) or change the algorithm used in fmincon

Your question is a bit strange though. bmibnb is a global solver (using, e.g., fmincon as a local solver), fmincon is a local solver.

[Edi]

Yeah I assigned an feasible solution at the beginning. Here is the code: 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N=20;         

M=10;        

  hk=[0.05 0.07 0.6 0.45 1.4 1.1 1 0.75 0.73  0.7 ];

  

No=0.1;

A=3;

sigma=0.36*ones(1,M);

   L=13.5;                                        

    

     sk=0.2*ones(1,M);

     

SNR_i=(sk.^2)./sigma.^2;

%==================================================================

t_N=4;

Lk=2;                           

Emax=1000*ones(t_N,M);                  

                             

Ek=10*rand(M,t_N);                    

%==========================================================================

pk = sdpvar(M,t_N);                            

assign(pk,.5);

for i=1:M

    for j=1:t_N

        

        rr(i,j)=((((1+(pk(i,j)*(hk(i)^2)/SNR_i(i)))^-Lk)));

        

   

    end

end

yprod=[];

for i=1:M

    yprod=[yprod prod(rr(i,:))];

end

ws=yprod;

for ii=1:M

    for jj=1:t_N

        

   

    Object_ive(ii)=(((3*N^2)*(sigma(ii)^4)*(SNR_i(ii)^2))./((6*N*(sigma(ii)^4).*(1+2*SNR_i(ii))+((A^2)*ws(ii))))); 

    

       end

end

 Objective_final=-sum(Object_ive);

 double(Objective_final)

Constraints = [(cumsum(Lk*pk'))<=(cumsum(Ek')), pk>=0,(cumsum(Ek'))-(cumsum(Lk*pk'))<=(Emax)]

 

 ops4= sdpsettings('solver','fmincon');

 solvesdp(Constraints,Objective_final,ops4);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[Johan Löfberg]

You have to tell YALMIP to use the initial solution

ops4= sdpsettings('solver','fmincon','usex0');

[Johan Löfberg]

should be

ops4= sdpsettings('solver','fmincon','usex0',1);