From: Oddvar G. <od...@pr...> - 2008-12-19 14:48:35
|
Changing the sign of the c vector did fix the odd results of the new solver. I had not noticed the change. Further, I found a bug in my own code (sorry!). I change the sign of the elements in c according to the sign of the corresponding vx2 variable, this is to force the variables towards zero. This causes oscillations around zero if the variable is (near) zero, vx2_0 in this case. If I detect oscillations like this I set the variable to zero and change the constraints 0.0 < variable < 0.0. The faulty detection of the oscillation was the reason the version with AS did not converge. After I fixed it, the AS results had the same accuracy as Matlabs, and about the same performance. (I got latest version from CVS today). I have done a few more test (changing tau) where the AS fails [60 10 0] and matlab does not. This can probably be fixed by the new Builder system, but I have to look into it some more later. Any hints on how to use the enforce method? - Oddvar -----Original Message----- From: Anders Peterson [mailto:an...@op...] Sent: 19. desember 2008 12:12 To: Anders Peterson Subject: Re: [ojAlgo-user] v25.1 On 18 dec 2008, at 20.15, Oddvar Grønning wrote: > Hi, > > I've been comparing the 25.1 ActiveSetSolver version with the one > currently in the CVS. Both of them are again compared to Matlabs > quadprog QP solver. > > My system is quite complex, I've reduced the example code as much > as I possibly can, but it is still quite large. You asked for test > cases, so here it goes.:-) It'll take me some time to study this. Don't think I have that time before Xmas. > I've made a SQP for nonlinear problems using the ActiveSetSolver as > QP engine. This means that I am linearizing and iterating till a > good enough solution is found. > > What I have found is that Matlab produces results very close to the > expected result, within 2 decimals. The AS 25.1 version produces > results within 1 decimal of expected results, but uses way to much > cpu-time. The latest AS version is fast, but the results are > different from Matlabs and the 25.1 version. See simulations > results below. I had set the iterations limit too tight - perhaps that was part of the problem. Also changed a couple of other things yesterday. Get the latest source from CVS. Have you tried the new builder/presolver? http://ojalgo.org/generated/org/ojalgo/optimisation/quadratic/ QuadraticSolver.Builder.html Calling "balance" and maybe "enforce" will change the quality of the solution. When you move from using the AS constructor with your BasicMatrix[] to using the QuadraticSolver.Builder you should change the sign of "C". I was able to get results, on your previous test cases, much closer to the MatLab-solution by simply "balancing" and "enforcing". But you should really try to do this earlier in your program. Once you end up with a number like -3.5E-15 you have already lost some precision. You can't just scale it back up to a larger number. > I think I need to explain the system a bit. There are 10 variables > to find optimum values for, the first 5 are low cost, they can have > any value within the constraints. The next 5 have high cost and > should be as close to zero as possible. There are nonlinear > relationships between the variables, defined as equality > constraints. 3 slack variables are used for the equality > constraints to ensure a feasible solution. This makes a total of 13 > variables.. > > The slack variables s are defined: > > s = tau - Beta(vx1)*vx2, where tau is a vector and Beta is a > matrix, we want s->0. vx1 are the first 5 low cost states and vx2 > are the high cost states. > > In the example I set tau to be [60 0 0]'. The expected result is: > > vx1: 1.5708; 0; 0; 0; 0; vx2: 0; 5; 5; 25; 25 > > In the first call to the SQP all states are zero. In the second > call the result variables from the first call are used as start > values. > > Matlab produces: > > i = 1 > > Solution found after 11 iterations > Result vx1: 1.5708; -1.0601e-005; -0.002093; 0.00073549; > -0.00031633; vx2: 0; 4.9991; 4.9947; 24.9996; 25.0065 > Elapsed time is 0.641000 seconds. > > i = 2 > > Solution found after 1 iterations > Result vx1: 1.5708; -0.0020056; -9.2987e-005; -0.0012645; > 0.0016837; vx2: -5.0796e-017; 4.9991; 4.9947; 24.9997; 25.0064 > Elapsed time is 0.015000 seconds. > > Version 25.1 of AS produces: > > Simstep 0 > Aborting, too long time consumed ( > 5 sec). numIt: 41 > Result: vx1_0: 1.5707672233183787 vx1_1: 1.6823620274180337E-4 > vx1_2: -0.00196511057420812 vx1_3: -0.002354032968486829 vx1_4: > 0.0015265945920139418 vx2_0: 0.029590231612839257 vx2_1: > 4.999569087984196 vx2_2: 4.990762180635852 vx2_3: 25.10133831261425 > vx2_4: 24.908290373430617 > Elapsed time is: 5.172 seconds > > Simstep 1 > Aborting, too long time consumed ( > 5 sec). numIt: 17 > Result: vx1_0: 1.5707672230998035 vx1_1: 1.6823670234216445E-4 > vx1_2: -0.001965110213385637 vx1_3: -0.0023328019515429285 vx1_4: > 0.0015265967291932642 vx2_0: 0.0290511558573982 vx2_1: > 4.999360102203392 vx2_2: 4.990900599987205 vx2_3: > 25.099688807764156 vx2_4: 24.910009190020308 > Elapsed time is: 5.015 seconds > > Latest AS produces: > > Simstep 0 > Solution found after 50 iteration(s) > Result: vx1_0: 1.5708199749750982 vx1_1: 0.0920372744518574 vx1_2: > 0.09394993946428473 vx1_3: -0.09067631072129541 vx1_4: > 0.0808235035333966 vx2_0: -1.427792446101636 vx2_1: > 4.557265828873959 vx2_2: 5.000000000931323 vx2_3: > 20.727837553658237 vx2_4: 29.939015120581345 > Elapsed time is: 1.578 seconds > > Simstep 1 > Solution found after 1 iteration(s) > Result: vx1_0: 1.5708199751971428 vx1_1: 0.09403727461728063 > vx1_2: 0.09594993965746354 vx1_3: -0.09267630822218338 vx1_4: > 0.08184556774440976 vx2_0: -1.4450338119993977 vx2_1: > 4.551633923447908 vx2_2: 5.000000000931323 vx2_3: > 20.680542518842206 vx2_4: 30.0 > Elapsed time is: 0.078 seconds > > The odd thing is that the results from the latest ActiveSetSolver, > is not all wrong, it meets the iteration stop requirement for s > (which the 25.1 version does not), but the solution is not optimal. > The c vector is designed to drive vx2 to a minimum (towards zero), > but this does not seem happen in this case. This is the part I'll have to look into a bit later... > I'll be happy to supply the Matlab code if you want. Don't have MatLab. I prefer junit test cases where you also supply the expected/MatLab solution and/or describes what goes wrong. > Java program was run on JRE1.6.0_01. I'm on Java 1.5 on a Mac. /Anders > - Oddvar > > -----Original Message----- >> From: Anders Peterson [mailto:an...@op...] > Sent: 16. desember 2008 14:47 > To: Anders Peterson > Subject: Re: [ojAlgo-user] v25.1 > > On 26 nov 2008, at 09.06, Oddvar Grønning wrote: > >> >>>> Do you have any hints/links that can put me on the right track to >>>> pre-solve the system myself? >> >>> Don't feel I know the best way to fix your problem. >> >>> The AE matrix contains some very small elements like for instance >>> -3.4914813388431334E-15. This is a constraint parameter?! Having >>> parameters with very different magnitude could be a problem, but in >>> this case the number itself is just too small. >> >> Ok, I see. Just so you don't think I take the numbers out of the >> air :-): That small number is the result of a cosine call with an >> angle close to pi/2.. The angle is actually a direction of a moving >> object, and can take the value of any number between - pi and pi. >> Such numbers are also common in control systems (i.e. Model >> Predictive Control) due to rounding errors. >> I'll try to remove such numbers before calling the QP solver. >> >>>> Or do you have any plans for developing a pre-solver? >> >>> Yes... > > > I'm interested in getting feedback on what's new / in CVS. > > /Anders > ---------------------------------------------------------------------- > -------- > SF.Net email is Sponsored by MIX09, March 18-20, 2009 in Las Vegas, > Nevada. > The future of the web can't happen without you. Join us at MIX09 > to help > pave the way to the Next Web now. Learn more and register at > http://ad.doubleclick.net/clk;208669438;13503038;i?http:// > 2009.visitmix.com/ > _______________________________________________ > ojAlgo-user mailing list > ojA...@li... > https://lists.sourceforge.net/lists/listinfo/ojalgo-user > > > <AllocationTest.java>------------------------------------------------- > ----------------------------- > SF.Net email is Sponsored by MIX09, March 18-20, 2009 in Las Vegas, > Nevada. > The future of the web can't happen without you. Join us at MIX09 > to help > pave the way to the Next Web now. Learn more and register at > http://ad.doubleclick.net/clk;208669438;13503038;i?http:// > 2009.visitmix.com/_______________________________________________ > ojAlgo-user mailing list > ojA...@li... > https://lists.sourceforge.net/lists/listinfo/ojalgo-user ------------------------------------------------------------------------------ _______________________________________________ ojAlgo-user mailing list ojA...@li... https://lists.sourceforge.net/lists/listinfo/ojalgo-user |