Can YALMIP be used to solve a general nonconvex constrained problem?

[zoharl]

Hi,

I have a general nonconvex function with general nonconvex inequality constraints (i.e. it can be anything, not necessarily quadratic - I'm emphasizing this since I only saw in the doc nonconvex QP). I have a feasible starting point, and I' like to minimize an energy. The solver should never leave the feasible domain (i.e. barrier method) and should never increase the energy. So far I used fmincon which failed on both accounts.

I'd like to try other solvers (any recommendations?), but I'm looking to model the problem once at some wrapper (e.g. YALMIP), which would translate it to the other solvers.

Zohar

[Johan Löfberg]

YALMP supports these nonlinear solvers

http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Category.NonlinearProgrammingSolver

However, not only the solver is important, but just as important the way you model things

[zoharl]

I've been through this part of the doc, but I don't seem to find what I want. Off the top of my head, let's consider some general problem:

min_x sin(x)^3 + sqrt(x) + log(x)

s.t.

1 < x

cos^x < x^3

I have a matlab code of the objective and constraint functions - callback which I provide to fmincon.

1. Do you have an example that shows how to write such a general problem in yalmip, let's say using IPOPT? In the past, I formulated SOCP problems with yalmip, but since it has its own variable type (which isn't necessary here, unless I'm looking for auto-diff), I'm not sure if it has the full flexibility of the language (i.e. I can do in the call back whatever I want and return real-double sparse matrices for the objective and constraints functions values and their gradients).

2. Do you happen to know which of the solvers (IPOPT, KNITRO, SNOPT...) support a barrier method (always stays in the feasible region and always decreases the function)?

[Johan Löfberg]

Nothing special here (except you have to use sqrtm to ensure yalmip treats the square root as a general nonlinear operator and doesn't try to do clever socp modelling)

sdpvar x

objective = sin(x)^3 + sqrtm(x) + log(x)

constraints = [1 <= x, cos(x) <= x^3]

% fmincon find a locally optimal solution with objective 2.68

optimize(constraints,objective,sdpsettings('solver','fmincon'))

% Globally optimally search reveals a solution with value 1.59

optimize([constraints,x<=100],objective,sdpsettings('solver','bmibnb'))

Why are you so hell-bent on using an algorithm which always stays in the feasible region? As long as ther final solution is feasible, why would you care?

[zoharl]

1. That's nice, but it still uses sdpvar. I have a long and complex code, and to make it sdpvar friendly won't be that easy (if it's even possible). 

I assume that yalmip would take care of the gradient and hessian?

2. fmincon at some point increases the function and can terminate with a worse point (larger objective), as well as compromise the constraints if terminated prematurely (it doesn't seem to converge). So it seems that I have a tough problem, and just strolling around the space, like fmincon does, won't do. (I thought that interior point supposes to be a barrier method).

[Johan Löfberg]

1. If it doesn't use sdpvar, by definition it isn't YALMIP, so I am not sure why you are discussing on a YALMIP forum. YALMIP computes gradients.

2. OK, might be algorithmic settings that you have to play around with.fmincon implemets a bunch of different algorithm, each with its specific setting etc. Interior-point methods are barrier methods yes, but that doesn't necessarily mean everything works in the interior of the original space. Could be interior in a dual space, or interior of some kind of relaxed space etc.

[zoharl]

Yalmip was the first wrapper that came to mind or otherwise, I thought someone might be able to point me in another direction. I'll post the question in a more general forum (Stackoverflow). Thanks for the tips, I'll consider them.

[Johan Löfberg]

try opti toolbox. it allows you to work on a low-level with callbacks etc, but still switch solver rather easily

[zoharl]

I experimented a little, and summarized my thoughts in

http://stackoverflow.com/questions/36782768/wrappers-for-solvers-of-general-nonconvex-nonlinear-constrained-problems-nlp-f/36845950#36845950

[Johan Löfberg]

Please post a model which I can test. Might be some simple fixes, or at least a good case for future development. Large-scale nonlinear programming is not the main focus of yalmip though

Your comment on interfacing via data as in mosek vs callbacks in ipopt is a bit strange. If you're using mosek, YALMIP (or any other interface) would never supply any derivatives, as all those a done internally in mosek directly from the problem definition (possible as all problems have the same structure), whereas when ipopt is used, the user, i.e. yalmip, has to supply a callback function to compute derivatives. this is always needed regardless of how you interface ipopt, but of course since yalmips auto-diff is a hack and runs in matlab code, it performance suffers on large problems

[zoharl]

I might try to factor my code out to a self-contained test, but we can't expect auto-diff to scale well anyway.

Maybe my phrasing could be better. I meant, as you said, that while yalmip can be perfect for mosek that doesn't require derivatives, this isn't the case in NLP that requires a callback. I'm not sure though what packages such as AMPL do. I assume it's not auto-diff, but rather a script that can be interpreted as an efficient code. Then I'm wondering if auto-diff can be 'translated' to such an efficient implementation (i.e. generate code instead of evaluating an expression tree). (I'm not fully sure about this at the moment and I'll need to think a bit.)

(BTW, I didn't mention that yalmip can write an AMPL model, because it's quite broken - uses norms and puts 'i' and AMPL can't handle them.)

[Johan Löfberg]

Whether an auto-diff from the expression tree can be efficient in yalmip depends on the model.

If you have a simple example where a yalmip model generates a weird model through saveampl, do post it

That holds in general: Development requires bug reports and performance reports from the wild

[zoharl]

Yes, I know, I'm sorry; I need to find the time to do that.

But about the saveampl, from the comment in the doc you didn't sound psyched about it. Anyway, it fails on the very basics; this

clear; clc;
n = 3;
x0 = 5*ones(n,1);
x = sdpvar(n,1);
assign(x, x0);
obj = norm(x); %obj = x'*x;
con = [0 <= x; 1 <= x.^2];
saveampl(con, obj, 'yal.run');

generates

var x {1..4};
minimize obj:  norm(x[1]);
subject to constr1: 0 <= x[1];
subject to constr2: 0 <= x[2];
subject to constr3: 0 <= x[3];
subject to constr4: 0 <= -1+x[1]^2;
subject to constr5: 0 <= -1+x[2]^2;
subject to constr6: 0 <= -1+x[3]^2;
solve;
display x;
display obj;

1. It doesn't initialize x (but I can fix that easily).

2. AMPL doesn't have norm. I can translate it if I look into AMPL - but I kind of trying to avoid that. In my complex model I don't think I'm using norm explicitly, but it's probably inferred somewhere. Also yalip adds complex numbers (that's my interpretation of 123i) - my code might involve complex numbers, I'll need to check, but they aren't supported in AMPL.

All-in-all, if I want to spend time on exporting to AMPL (which I currently just trying the demo), the OptiToolbox export to GAMS (and from there to AMPL) might provide a better starting point.

[Johan Löfberg]

Looks like a consequence of the fact that norm is primarily supposed to be used in convex programs with socp reformulations etc, so it gets compeltely confused here

What do you mean with 123i? I don't see that

[zoharl]

About the complex numbers, I was referring to my real code. But for my mock example, try:

obj = norm(123i*x);

which generates

minimize obj:  norm(0+123i*x[1]);