Re: [ojAlgo-user] v25.1

[Anders P. <an...@op...>]

On 19 dec 2008, at 15.47, Oddvar Grønning wrote:

> 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).


Having roughly the same performance and accuracy as MatLab seems ok  
to me.


> 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?


I'm testing a new/simpler/different implementation of  
NumberContext.enforce(double), and changed Builder.enforce 
(NumberContext) to Builder.round(int) - they do the same thing with  
this new implementation. They simply round to a specified number of  
decimals.

Note that if you call more than one of the methods balance, eliminate  
and round it matters in which order you do it.  There is no correct  
or wrong order; just different results.


/Anders


> -----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
>
>
> <AllocationTest.java>------------------------------------------------- 
> -----------------------------
> _______________________________________________
> ojAlgo-user mailing list
> ojA...@li...
> https://lists.sourceforge.net/lists/listinfo/ojalgo-user