cureos / jcobyla

Java port of COBYLA2 nonlinear constrained optimization method
35 stars 17 forks source link

Successful termination when constraints violated #1

Closed stkiese closed 8 years ago

stkiese commented 9 years ago

I have a similar issue as described in https://github.com/scipy/scipy/issues/2891 for the python port of cobyla2: the minimization can terminate without actually satisfying the constraints.

I've uploaded an example to pastebin: http://pastebin.com/B7Z0Y95h

In my case, reducing rhobeg (solution in https://github.com/scipy/scipy/issues/2891 ) won't help. With rhobeg=1 the constraints are violated by 0.091. With rhobeg=0.001 the constraints are violated by 0.254. With rhobeg=0.00001 the constraints are violated by 0.267 ... (with simultaneous increase from 710 to 27153 iterations).

In my opinion this is a bug.

Thanks Stefan

anders9ustafsson commented 9 years ago

Which constraint(s) is/are violated? RHOBEG 1 sounds large if the sum of variables should be unity, but on the other hand RHOBEG 0.001 or smaller sounds relatively small to start with. How large violations do you see if you start with RHOBEG 0.1 for example?

stkiese commented 9 years ago

Thanks for your fast reply :-). It seems that the

6th (the 6th input var has to be greater or equal to zero), (9th or 10th) (the sum of all input vars has to be approximately 1) and 11st (to complex to explain in one sentence ;-) )

are violated. It depends on rhobeg if the 9th or the 10th constraint is violated.

rhobeg=0.75 => MAXCV=0.0543 rhobeg=0.5 => MAXCV=0.0287 rhobeg=0.25 => MAXCV=0.0685 rhobeg=0.1 => MAXCV=0.0925 rhobeg=0.05 => MAXCV=0.0735

Is there any rule or recommendation how to choose rhobeg?

Thanks!

anders9ustafsson commented 9 years ago

If I am not mistaken, RHOBEG should be a suitable step size for variable changes in the first round of iterations, so if the sum of your 8 variables should be equal to unity, a first estimate for RHOBEG would probably be around 0.1. Apparently RHOBEG 0.5 seem to give the least violation in your case, though.

COBYLA is probably not really suited for equality constraints, so I guess this is why you see these violations. If you would change the equality tolerance to for example 0.02 instead, i.e. 0.99 <= sum <= 1.01, would this improve the level of violation?

stkiese commented 9 years ago

Thanks for the explanation.

Unfortunately this is not improving the level of violation. Even if I check for 0.9 <= sum <= 1.1 or set the 9th and/or 10th constraint to 0, the violation level is 0.0x for different rhobeg.

By the way, there is a solution that satisfies all constraints:

{0.705, 0.0075, 0.075, 0, 0.0375, 0, 0.05, 0.125}

But even if I start with this vector, I have to decrease rhobeg to 0.0002 to get a solution with a violation level lower than 0.001.

I'm using JCobyla in a program that generates (based on user input) a lot of different optimization problems like my example. The problems differ in the matrix, the maximum and the used vectornorm in constraint 11 and in the number of variables. Therefore it is very unsatisfying if I have to try different rhobeg and input vectors to get a suitable result.

Do you have any idea how to deal with this? Is (J)Cobyla just unsuitable for that kind of optimization problem? Or does JCobyla calculate solutions only approximated and does not guarantee to satisfy the constraints? I'm also able to provide more examples.

Thank you!

anders9ustafsson commented 9 years ago

I cannot say for certain, but I would think that COBYLA would have difficulties with equality constraints in general. In principle COBYLA will satisfy all constraints after a sufficient number of iterations, but it might get trapped in an infeasible solution that it cannot escape from, depending for example on the RHO update procedure.

You might want to try another optimizer, although I cannot name another Java implementation that supports nonlinear constraints out of my head. Maybe some global optimizer is available elsewhere?

stkiese commented 9 years ago

ok, thanks for your explanation and your time! I've tried another java cobyla implementation https://github.com/xypron/jcobyla but this implementation exhibits the same behavior. Therefore I think that this really could be a general issue of cobyla. Or do I use JCobyla in a wrong way?

I've also created an issue for the xypron implementation. Maybe xypron has an additional idea.

anders9ustafsson commented 9 years ago

The xypron implementation appears to be a direct copy of this implementation, so you are expected to exhibit the same behavior.

I do not see anything immediately wrong in the way you are using JCobyla. Since it is a fairly small problem, would you be able to test it using other language implementations of COBYLA, for example the C implementation or C#? (If you try out C# though, please be aware that the Java implementation originates from the C# implementation.)

In the C# library CSNumerics you'll also find alternative solvers, for example the LINCOA algorithm that supports nonlinear objective and linear constraints. I haven't examined it too closely, but if also your last constraint is linear, LINCOA might be an alternative solver. Unfortunately, there is not yet any Java implementation of LINCOA as far as I know.

Best regards, Anders @ Cureos

stkiese commented 9 years ago

I've tried your C# implementation with different rhobeg

http://pastebin.com/SaSgGfyU

The results are also violating the constraint but differ from JCobyla in some cases:

rhobeg=1 => MAXCV=0.108 rhobeg=0.75 => MAXCV=0.054 rhobeg=0.5 => MAXCV=0.034 rhobeg=0.25 => MAXCV=0.039 rhobeg=0.1 => MAXCV=0.092 rhobeg=0.05 => MAXCV=0.073 rhobeg=0.001 => MAXCV=0.255

Furthermore I've tried the C# implementation NLoptNet https://github.com/BrannonKing/NLoptNet . I haven't figured out how to set rhobeg in this implementation, but after testing with some different input vectors, it seems that rhobeg is set to 1. My implementation

http://pastebin.com/hVugXh7s

calculates for every tested input vector the same maximal violation as the csnumerics implementation with rhobeg=1 does. In difference to your implementation, NLoptNet returns as result infinity. I think NLoptNet recognizes that the optimization has failed.

Unfortunately I have no C compiler installed. Also the C implementation looks quite different than the Java or C# implementation. Therefore and due to I need a Java optimizer, I haven't tried the C implementation.

In two special cases, the last constraint is actually linear. These are also the two special cases that are causing the worst results. Therefore it could actually be an improvement to use for these cases LINCOA. Unfortunately I have to use Java. Do you know a Java implementation of LINCOA?

Thanks Stefan

anders9ustafsson commented 9 years ago

As far as I know, no-one else has ported LINCOA to Java. I have planned to do it myself eventually using the C# code as a base, but I haven't gotten that far yet.

Maybe you can test your specific problem on C# LINCOA to see if this algorithm would solve your problems. If it does, maybe you could take on the effort yourself to port LINCOA to Java?

Best regards, Anders @ Cureos

stkiese commented 9 years ago

I had to redefine the optimization problem to ensure compatibility with lincoa:

http://pastebin.com/n921tAcQ

If the start vector isn't violating the constraints, lincoa solves the problem successfully. If the start vector violates the constraints, lincoa is also not able to solve the problem without violating one or more constraints. Therefore I don't think that porting lincoa to java would solve my issue :-(.

Kind regards Stefan

anders9ustafsson commented 9 years ago

I am sorry to hear that LINCOA did not work out well for you, but I think it is great that you have been able to set up this test case, many thanks for this excellent example. Would you mind if I included it in the unit tests for LINCOA (in case the code eventually is improved to solve your type of problem)? You will of course be sufficiently acknowledged in the source code.

Thanks in advance, and best regards, Anders @ Cureos

stkiese commented 9 years ago

No, I wouldn't mind if you include my example as unit test for LINCOA. Do you need any additional information like the formal problem description?