Optimization problem

[David]

Hello,

I come across the following optimization problem, I need any help to solve it using MATLAB

maximize a*ln(x(1)*(1-x(2)^2)) + b*ln(x(3)*(1-x(4)^2)) + c*ln(x(5))

subject to (a/2)*x(1) + (b/2)*x(3) + (c/2)*x(5) = 1

                -1<=x(2)<=1 ;  -1<=x(4)<=1

                 x(1)>=0; x(3)>=0; x(5>=0)

                 c_1*sqrt(x(1))*x(2) + c_2*sqrt(x(3))*x(4) - c_3*x(1) - c(4)*x(3) + c(5) = 0

where {a,b,c,c_1,c_2,c_3,c_4,c_5} are some positive constants. 

Thank you

[Johan Löfberg]

And where lies your problem?

A basis would be something like this, basically a direct copy of your model

a = rand

b = rand

c = rand

c_1 = rand

c_2 = rand

c_3 = rand

c_4 = rand

c_5 = rand

x = sdpvar(5,1);

Objective = a*log(x(1)*(1-x(2)^2)) + b*log(x(3)*(1-x(4)^2)) + c*log(x(5));

Constraints = [(a/2)*x(1) + (b/2)*x(3) + (c/2)*x(5) == 1,

                -1<=x(2)<=1,

                -1<=x(4)<=1

                 x(1)>=0,

                 x(3)>=0,

                 x(5)>=0,

                 c_1*sqrtm(x(1))*x(2) + c_2*sqrtm(x(3))*x(4) - c_3*x(1) - c_4*x(3) + c_5 == 0]

assign(x,.5)

solvesdp(Constraints,-Objective,sdpsettings('usex0',1))         

Things to note

1. ln is called log in MATLAB

2. You had a variable c(5). I guess you meant c_5

3. sqrt is a convexity aware operator implemented using SOCPs in YALMIP. In this general NLP settings, the sqrtm operator should be used instead.

http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Extra.SQRT

4. Many solvers will struggle with the nonlinear expressions in the log operator, since it easily can happen that the iterate leads to a term which is negative, leading to NaNs in the objective. By initializing the solution to a point far from these critical points. Both fmincon and IPOPT manages to solve the problem easily after this initialization, but fails without it.

[Johan Löfberg]

...and it turns out that the model is small enough for the global solver in YALMIP to be applicable

solvesdp(Constraints,-Objective,sdpsettings('solver','bmibnb','usex0',1,'bmibnb.maxiter',1000))    

With an efficient LP solver installed, YALMIP often declares global optimality in roughly 10 seconds (global solution coincides with the solution fmincon and ipopt finds in the first iteration, indicating this is a pretty simple nonlinear problem (of course this depends on the data, different signs on a b c would probably lead to different behaviour))

[Johan Löfberg]

...and note that the global solver performs much much better if you expand the log terms

Objective = a*(log(x(1))+log(1-x(2)^2)) + b*(log(x(3))+log(1-x(4)^2)) + c*log(x(5));

The convex envelopes are much more tight for this model it looks like.

[David]

Thanks !

I think that the problem was issue number 4 that you have mentioned. I tried to use fmincon and something didn't work. 

So, all I need is to download 'yalmip' software and use your proposed code? BTW, can you recommend for an efficient LP solver?

Thanks a LOT

[Johan Löfberg]

Yes, just download YALMIP are you are almost ready to go. I tried using linprog as a lower bound LP solver here, but it doesn't work. It is not robust enough, so I recommend a better LP solver.

If you are in academia, I recommend you to install Gurobi, Cplex, Mosek or Xpress. They all have academic licenses (full versions free for academics). I used Gurobi here. (Mosek is problematic here since it overloads some fmincon stuff and thus we cannot use fmincon, but if you use ipopt instead as a nonlinear solver it works)

If you want to try another NLP solver than fmincon (which I used here), I recommend you to download the OPTI Toolbox. It comes with precompiled binaries for a large number of solver, such as the nonlinear solver IPOPT (which isn't any faster than fmincon on this small problem). OPTI Toolbox also comes with a range of LP solver. I tried to use them here and they do work, but some of them seem to be a bit shaky, and some work but run almost an order of magnitude slower than Gurobi, which makes the global search slow, as 50% of the time is spent in computing lower bounds using LPs.

Note though, this whole discussion is mainly relevant if you want to use the global solver. It looks like a local solver almost always returns the global solution anyway, so unless you want a guaranteed global solution, you can just use fmincon or ipopt directly. The reason is that the objective actually is convex. The last constraint is nonconvex, but maybe it is easy to show that you can relax the equality to an >=, which leads to a convex set. When I did this, the solution was still tight at the optima.

[David]

Johan,

I tried to use the following code (as you suggested)

x = sdpvar(5,1);

          

            Objective = 0.25*(log(x(1))+log(1-x(2)^2)) + ((frq-0.5*pi)/(2*pi))*(log(x(3))+log(1-x(4)^2)) + ((pi-frq)/(2*pi))*log(x(5));

            Constraints = [0.5*x(1) + (frq-0.5*pi)/(pi)*x(3) + (pi-frq)/(pi)*x(5) == 1,

                -1<=x(2)<=1,

                -1<=x(4)<=1

                 x(1)>=0,

                 x(3)>=0,

                 x(5)>=0,

                 0.5*sqrtm(x(1))*x(2)*sqrt(1+1/var) + ((frq-0.5*pi)/(pi))*sqrtm(x(3))*x(4)*sqrt(1/var) - 0.25*(x(1)-1) -((frq-0.5*pi)/(2*pi))*(x(3)-1) -0.5 == 0]

            assign(x,.5)

and first I use

            solvesdp(Constraints,-Objective,sdpsettings('usex0',1))

and then

            solvesdp(Constraints,-Objective,sdpsettings('solver','bmibnb','usex0',1,'bmibnb.maxiter',1000)) 

For some reason I obtained two different maximizers. In particaular, the solution of the second option is [0 0 0 0 0] which is not reasonable due to the first constraint. What i'm doing wrong?

Thanks

[Johan Löfberg]

Show the complete output displayed by bmibnb, i.e., something like

* Starting YALMIP global branch & bound.

* Branch-variables : 7

* Upper solver     : fmincon

* Lower solver     : LINPROG

* LP solver        : LINPROG

 Node       Upper      Gap(%)       Lower    Open

    1 :    3.148E-07     0.00      3.148E-07   0  Infeasible  

* Finished.  Cost: 3.1483e-07 Gap: 0

* Timing: 0% spent in upper solver

*         0% spent in lower solver

*         11% spent in LP-based domain reduction

I get the same solution using both local and global approach

>> double(x)

ans =

    1.0002

    0.0008

    0.9987

    0.0005

    1.0008

[David]

In the first approach I get the same results you do.

In the second approach

* Starting YALMIP global branch & bound.

* Branch-variables : 7

* Upper solver     : fmincon

* Lower solver     : LINPROG

* LP solver        : LINPROG

 Node       Upper      Gap(%)       Lower    Open

    1 :          Inf      Inf            Inf   0  Infeasible  

* Finished.  Cost: Inf Gap: Inf

* Timing: 0% spent in upper solver

*         0% spent in lower solver

*         15% spent in LP-based domain reduction

ans = 

    yalmiptime: 0.5620

    solvertime: 0.4460

          info: 'Infeasible problem (BMIBNB)'

       problem: 1

K>> double(x)

ans =

     0

     0

     0

     0

     0

BTW, I know that it is a stupid question, but regard you last comment 

"Note though, this whole discussion is mainly relevant if you want to use the global solver. It looks like a local solver almost always returns the global solution anyway, so unless you want a guaranteed global solution, you can just use fmincon or ipopt directly. The reason is that the objective actually is convex. The last constraint is nonconvex, but maybe it is easy to show that you can relax the equality to an >=, which leads to a convex set. When I did this, the solution was still tight at the optima."

When you say "just use fmincon", you actually meant that I need to use your first approach - which performed by yalmip software using fmincon, right? Because, any previous experience with fmincon (as is) was unsuccessful for me :(

[Johan Löfberg]

Looks like LINPROG fails miserably on your machine claiming the root relaxation is infeasible, which it isn't of course (which MATLAB/OS are you using?). As I said, install gurobi or equivalent solver if you intend to solve LPs/QPs etc in the future

Yes, I mean use fmincon using YALMIP, i.e., do not use bmibnb to explicitly search for a global solution. You will get a global solution anyway most likely, since the problem is convex (almost, disregarding the last equality, which it looks like you can relax). The fact that you haven't managed to compute solutions using fmincon directly is probably due to poor initialization, maybe you supplied incorrect derivatives, etc. In addition to that, YALMIP performs some clean up of the model, and tries to put it in a format that fmincon likes. Many users claim fmincon works terribly, but with a well formulated model, I cannot agree on that. It typically performs pretty well when used through YALMIP.

[David]

I use MATLAB R2012a and Windows 7

OK great! you helped me very much!

Best

[Johan Löfberg]

That is odd, as I use the same here. Anyway, LINPROG is clearly not robust enough.

[David]

Johan,

Now I noticed to the following problem. Using the previous code:

x = sdpvar(5,1);

          

            Objective = 0.25*(log(x(1))+log(1-x(2)^2)) + ((frq-0.5*pi)/(2*pi))*(log(x(3))+log(1-x(4)^2)) + ((pi-frq)/(2*pi))*log(x(5));

            Constraints = [0.5*x(1) + (frq-0.5*pi)/(pi)*x(3) + (pi-frq)/(pi)*x(5) == 1,

                -1<=x(2)<=1,

                -1<=x(4)<=1

                 x(1)>=0,

                 x(3)>=0,

                 x(5)>=0,

                 0.5*sqrtm(x(1))*x(2)*sqrt(1+1/var) + ((frq-0.5*pi)/(pi))*sqrtm(x(3))*x(4)*sqrt(1/var) - 0.25*(x(1)-1) -((frq-0.5*pi)/(2*pi))*(x(3)-1) -0.5 == 0]

            assign(x,.5)

            solvesdp(Constraints,-Objective,sdpsettings('usex0',1))

I obtain the solution

1.0002

0.0008

0.9987

0.0005

1.0008

as you did. However, this solution is not satisfying the nonlinear constraint:

0.5*sqrtm(x(1))*x(2)*sqrt(1+1/var) + ((frq-0.5*pi)/(pi))*sqrtm(x(3))*x(4)*sqrt(1/var) - 0.25*(x(1)-1) -((frq-0.5*pi)/(2*pi))*(x(3)-1) -0.5 == 0

(the left term is equal to minus half using the above solution).

Where is my mistake?

[David]

Weird !

I changed 

Constraints = [0.5*x(1) + (frq-0.5*pi)/(pi)*x(3) + (pi-frq)/(pi)*x(5) == 1,

                 -1<=x(2)<=1,

                 -1<=x(4)<=1,

                  x(1)>=0,

                  x(3)>=0,

                  x(5)>=0,

                  0.5*sqrtm(x(1))*x(2)*sqrt(1+1/var) + ((frq-0.5*pi)/(pi))*sqrtm(x(3))*x(4)*sqrt(1/var) - 0.25*(x(1)-1) -((frq-0.5*pi)/(2*pi))*(x(3)-1) -0.5 == 0]

to the following equivalent form

Constraints = set(0.5*sqrtm(x(1))*x(2)*sqrt(1+1/var) + ((frq-0.5*pi)/(pi))*sqrtm(x(3))*x(4)*sqrt(1/var) - 0.25*(x(1)-1) -((frq-0.5*pi)/(2*pi))*(x(3)-1) -0.5 == 0)

Constraints = Constraints + set(0.5*x(1) + (frq-0.5*pi)/(pi)*x(3) + (pi-frq)/(pi)*x(5) == 1);

and now I obtain another solution, which satisfy both constraints. 

[Johan Löfberg]

It's due to misplaced spaces in your code. Note the difference in matlab between these three expression

[1 +2]
[1+2]
[1 + 2]

Your last equality evaluates to a vector with three elements