commons-user mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Luc Maisonobe <...@apache.org>
Subject [math] Re: solving a linear optimization problem with commons math ?
Date Fri, 04 Jan 2008 20:13:42 GMT
Andreas Andreakis a écrit :
> hi,
> 
> sorry to send you an email regarding commons math directly. I tried to 
> subscribe to the mailing list, but for some reason, I was considered as 
> spam.
> Thus I chose two members randomly.
> 
> Well, anyway :) I want to find a solution to the following problem:
> 
> (explained by example)
> 1 *c1 + 4 * c2 + 4* c3 + 4 * c4 = 2.5
> 1 *c1 + 4 * c2 + 1* c3 + 4 * c4 = 3.0
> 
> here the constraints:
> 1) c1+c2+c3+c4 = 1 (of course this could be expressed as an additional 
> equation: c1 + c2 + c3 + c4 = 1.0)
> 2) c1,c2,c3 and c4 are all positive
> 
> Important is that I need only one and feasible solution, thus the 
> computation should be as fast as possible. So there is no need to search 
> for every possible result.
> And it needs to work with doubles (not only integers, because some 
> libraries work only with integers)
> 
> Actually the solution needs to scale for several dozens c´s ( there will 
> not be only 4 like c1,c2,c3,c4), but more like 50 c1, c2,...,c49, c50
> And there will not be only two equations (like in the example) but 
> several hundred.
> 
> 
> So the questions are:
> 1) can I find a solution using commons math. There is an optimization 
> package, but I could not figure out how to use it?

In this case, the estimation package is more suited to your needs. There 
are much better methods for such problems, but they are not available in 
commons-math.

The Levenberg-Marquardt solver implemented in thie estimation package 
should support both under- and over-determined problems (i.e. either 
more equations than variables or more variables than equations). The 
case with more variable than equations is still considered experimental, 
though. I did use it for such non homogeneous linear problems without 
problems myself.

This solver does not completely fit your needs as is because it 
currently do not handle constraints. You will have to set up an 
intermediate unconstrained problem that can be solved and from which the 
solution of your constrained problem can be derived. Here are a few 
tricks to do that, take care, they may induce numerical problems:

For equality constraints, either add the constraint to the equation set 
or remove one parameter (say c1) from the parameter set and compute it 
as needed by solving the constraint each time from the values of the 
other ones, i.e. compute c1 = 1 - f(c2, c3, ..., cn).

For simple bounds constraints on parameters (such as ck > ak), replace 
the real parameter ck by an unconstrained parameter pk using
   ck = ak + pk² (or ak - pk² if the constraint was ck < ak)

For interval constraints on parameters (such as lowk < ck < highk), 
replace the real parameter ck by an unconstrained parameter pk using
   ck = [ lowkk + highk ] / 2 + atan(pk) * [ highk - lowk ] / pk

Once you have changed your problem using these pk, it is an 
unconstrained problem that can be solved by our implementation of 
Levenberg-Marquardt. Once you have the pk values, you compute the ck values.


> 2) how will the performance be ?

Since the problem is linear, the performance should be fair.

> 3) could you give me a code snipped which is able to compute the given 
> example ? (I will generalize it then myself)

Try something along these lines. Take care, I wrote this exemple 
directly in the mail, I did not attempt to even compile it ;-)

Hope this helps,
Luc

import org.apache.commons.math.estimation.EstimatedParameter;
import org.apache.commons.math.estimation.EstimationException;
import org.apache.commons.math.estimation.LevenbergMarquardtEstimator;
import org.apache.commons.math.estimation.SimpleEstimationProblem;
import org.apache.commons.math.estimation.WeightedMeasurement;


public class LinearConstrainedProblem extends SimpleEstimationProblem {

     public static void main(String[] args) {

         try {
             // create the uninitialized problem
             LinearConstrainedProblem pb =
                new LinearConstrainedProblem(0.25, 0.25, 0.25, 0.25)

             // 1 *c1 + 4 * c2 + 4* c3 + 4 * c4 = 2.5
             pb.addEquation(1.0, 4.0, 4.0, 4.0, 2.5);

             // 1 *c1 + 4 * c2 + 1* c3 + 4 * c4 = 3.0
             pb.addEquation(1.0, 4.0, 1.0, 4.0, 3.0);

             // solve the problem, using a Levenberg-Marquardt algorithm
             // with default settings
             LevenbergMarquardtEstimator estimator =
               new LevenbergMarquardtEstimator();
             estimator.estimate(pb);
             System.out.println("RMS after solve = "
                              + estimator.getRMS(pb));

             // display the results
             System.out.println("c1 = " + pb.getC1()
                     + ", c2 = " + pb.getC2()
                     + ", c3 = " + pb.getC3()
                     + ", c4 = " + pb.getC4());

         } catch (EstimationException ee) {
             System.err.println(ee.getMessage());
         }

     }

     /** Fitter for a linear constrained model.
      * @param c1Guess start value for c1 (in fact, ignored)
      * @param c2Guess start value for c2
      * @param c3Guess start value for c3
      * @param c4Guess start value for c4
      */
     public LinearConstrainedProblem() {

         // one parameter (c1) is computed from the other ones
         // and the other ck are derived from pk using ck = pk*pk
         p2 = new EstimatedParameter("p2", Math.sqrt(c2Guess));
         p3 = new EstimatedParameter("p3", Math.sqrt(c3Guess));
         p4 = new EstimatedParameter("p4", Math.sqrt(c4Guess));

         // provide the parameters to the base class which
         // implements the getAllParameters and
         // getUnboundParameters methods
         addParameter(c2);
         addParameter(c3);
         addParameter(c4);

     }

     /** Add an equation.
      * @param a1 coeff for c1
      * @param a2 coeff for c2
      * @param a3 coeff for c3
      * @param a4 coeff for c4
      * @param rhs right hand side of the equation
      */
     public void addEquation(double a1, double a2, double a3,
                             double a4, double rhs) {
         addMeasurement(new Equation(a1, a2, a3, a4, rhs));
     }

     public double getC1() {
         // this solves the equality constraint
         return 1.0 - getC2() - getC3() - getC4();
     }

     public double getC2() {
       double p = p2.getEstimate();
       return p * p;
     }

     public double getC3() {
       double p = p3.getEstimate();
       return p * p;
     }

     public double getC4() {
       double p = p4.getEstimate();
       return p * p;
     }


     /** Internal measurements class, representing one equation. */
     private class Equation extends WeightedMeasurement {

         public Equation(double a1, double a2, double a3,
                         double a4, double rhs) {
             super(1.0, rhs);
             this.a1 = a1;
             this.a2 = a2;
             this.a3 = a3;
             this.a4 = a4;
         }

         public double getTheoreticalValue() {
             return a1 * getC1() + a2 * getC2()
                  + a3 * getC3() + a4 * getC4();
         }

         public double getPartial(EstimatedParameter parameter) {
           if (parameter == p2) {
             return 2.0 * p2.getEstimate() * (a2 - a1);
           } else if (parameter == p3) {
             return 2.0 * p3.getEstimate() * (a3 - a1);
           } else {
             return 2.0 * p4.getEstimate() * (a4 - a1);
           }
         }

         private final double a1;
         private final double a2;
         private final double a3;
         private final double a4;

     }

     private EstimatedParameter p2;
     private EstimatedParameter p3;
     private EstimatedParameter p4;

}

> 
> kind regards,
> Andreas



---------------------------------------------------------------------
To unsubscribe, e-mail: user-unsubscribe@commons.apache.org
For additional commands, e-mail: user-help@commons.apache.org


Mime
View raw message