Example usage for org.apache.commons.math.optimization.linear SimplexSolver optimize

List of usage examples for org.apache.commons.math.optimization.linear SimplexSolver optimize

Introduction

In this page you can find the example usage for org.apache.commons.math.optimization.linear SimplexSolver optimize.

Prototype

public RealPointValuePair optimize(final LinearObjectiveFunction f,
        final Collection<LinearConstraint> constraints, final GoalType goalType,
        final boolean restrictToNonNegative) throws OptimizationException 

Source Link

Usage

From source file:net.sf.tweety.math.opt.solver.ApacheCommonsSimplex.java

@Override
public Map<Variable, Term> solve(ConstraintSatisfactionProblem problem) {
    if (!problem.isLinear())
        throw new IllegalArgumentException("Simplex algorithm is for linear problems only.");
    //this.log.info("Wrapping optimization problem for calling the Apache Commons Simplex algorithm.");            
    // 1.) bring all constraints in linear and normalized form
    Set<Statement> constraints = new HashSet<Statement>();
    for (Statement s : problem)
        constraints.add(s.toNormalizedForm().toLinearForm());
    // 2.) for every constraint we need an extra variable
    int numVariables = problem.getVariables().size();
    // 3.) define mappings from variables to indices
    int index = 0;
    Map<Variable, Integer> origVars2Idx = new HashMap<Variable, Integer>();
    for (Variable v : problem.getVariables())
        origVars2Idx.put(v, index++);//  w  w w  .j  a v a2 s  . c o  m
    // 4.) Check for target function (for constraint satisfaction problems
    //      its empty
    double[] coefficientsTarget = new double[numVariables];
    int i = 0;
    for (; i < numVariables; i++)
        coefficientsTarget[i] = 0;
    double constTerm = 0;
    if (problem instanceof OptimizationProblem) {
        // bring target function in linear form
        Sum t = ((OptimizationProblem) problem).getTargetFunction().toLinearForm();
        for (Term summand : t.getTerms()) {
            // as t is in linear form every summand is a product
            Product p = (Product) summand;
            if (p.getTerms().size() == 1) {
                // p consists of just a constant term
                Constant c = (Constant) p.getTerms().get(0);
                if (c instanceof FloatConstant)
                    constTerm += ((FloatConstant) c).getValue();
                else
                    constTerm += ((IntegerConstant) c).getValue();
            } else {
                // p consists of a variable and a constant
                Variable v = (Variable) ((p.getTerms().get(0) instanceof Variable) ? (p.getTerms().get(0))
                        : (p.getTerms().get(1)));
                Constant c = (Constant) ((p.getTerms().get(0) instanceof Constant) ? (p.getTerms().get(0))
                        : (p.getTerms().get(1)));
                double coefficient = (c instanceof FloatConstant) ? (((FloatConstant) c).getValue())
                        : (((IntegerConstant) c).getValue());
                coefficientsTarget[origVars2Idx.get(v)] += coefficient;
            }
        }
    }
    LinearObjectiveFunction target = new LinearObjectiveFunction(coefficientsTarget, constTerm);
    // 5.) Represent the constraints
    Set<LinearConstraint> finalConstraints = new HashSet<LinearConstraint>();
    for (Statement s : constraints) {
        double[] coefficientsConstraint = new double[numVariables];
        for (i = 0; i < numVariables; i++)
            coefficientsConstraint[i] = 0;
        // as s is in linear form go through the summands
        Sum leftTerm = (Sum) s.getLeftTerm();
        double rest = 0;
        for (Term summand : leftTerm.getTerms()) {
            // as s is in linear form every summand is a product
            Product p = (Product) summand;
            if (p.getTerms().size() == 1) {
                // p consists of just a constant term
                Constant c = (Constant) p.getTerms().get(0);
                if (c instanceof FloatConstant)
                    rest += ((FloatConstant) c).getValue();
                else
                    rest += ((IntegerConstant) c).getValue();
            } else {
                // p consists of a variable and a constant
                Variable v = (Variable) ((p.getTerms().get(0) instanceof Variable) ? (p.getTerms().get(0))
                        : (p.getTerms().get(1)));
                Constant c = (Constant) ((p.getTerms().get(0) instanceof Constant) ? (p.getTerms().get(0))
                        : (p.getTerms().get(1)));
                double coefficient = (c instanceof FloatConstant) ? (((FloatConstant) c).getValue())
                        : (((IntegerConstant) c).getValue());
                coefficientsConstraint[origVars2Idx.get(v)] += coefficient;
            }
        }
        Relationship r = Relationship.EQ;
        if (s instanceof Inequation)
            r = Relationship.GEQ;
        finalConstraints.add(new LinearConstraint(coefficientsConstraint, r, -rest));
    }
    // 6.) Optimize.
    try {
        //this.log.info("Calling the Apache Commons Simplex algorithm.");
        SimplexSolver solver = new SimplexSolver(0.01);
        solver.setMaxIterations(this.MAXITERATIONS);
        RealPointValuePair r = null;
        if (problem instanceof OptimizationProblem) {
            int type = ((OptimizationProblem) problem).getType();
            r = solver.optimize(target, finalConstraints,
                    (type == OptimizationProblem.MINIMIZE) ? (GoalType.MINIMIZE) : (GoalType.MAXIMIZE),
                    this.onlyPositive);
        } else
            r = solver.optimize(target, finalConstraints, GoalType.MINIMIZE, this.onlyPositive);
        //this.log.info("Parsing output from the Apache Commons Simplex algorithm.");
        Map<Variable, Term> result = new HashMap<Variable, Term>();
        for (Variable v : origVars2Idx.keySet())
            result.put(v, new FloatConstant(r.getPoint()[origVars2Idx.get(v)]));
        return result;
    } catch (OptimizationException e) {
        //log.error(e.getMessage());
        throw new ProblemInconsistentException();
    }
}

From source file:emlab.role.AbstractEnergyProducerRole.java

/**
 * The fuel mix is calculated with a linear optimization model of the possible fuels and the requirements.
 * //from   w  w  w  .  java  2 s.  c om
 * @param substancePriceMap
 *            contains the possible fuels and their market prices
 * @param minimumFuelMixQuality
 *            is the minimum fuel quality needed for the power plant to work
 * @param efficiency
 *            of the plant determines the need for fuel per MWhe
 * @param co2TaxLevel
 *            is part of the cost for CO2
 * @param co2AuctionPrice
 *            is part of the cost for CO2
 * @return the fuel mix
 */
public Set<SubstanceShareInFuelMix> calculateFuelMix(PowerPlant plant, Map<Substance, Double> substancePriceMap,
        double co2Price) {

    double efficiency = plant.getActualEfficiency();

    Set<SubstanceShareInFuelMix> fuelMix = (plant.getFuelMix() == null) ? new HashSet<SubstanceShareInFuelMix>()
            : plant.getFuelMix();

    int numberOfFuels = substancePriceMap.size();
    if (numberOfFuels == 0) {
        logger.info("No fuels, so no operation mode is set. Empty fuel mix is returned");
        return new HashSet<SubstanceShareInFuelMix>();
    } else if (numberOfFuels == 1) {
        SubstanceShareInFuelMix ssifm = null;
        if (!fuelMix.isEmpty()) {
            ssifm = fuelMix.iterator().next();
        } else {
            ssifm = new SubstanceShareInFuelMix().persist();
            fuelMix.add(ssifm);
        }

        Substance substance = substancePriceMap.keySet().iterator().next();

        ssifm.setShare(calculateFuelConsumptionWhenOnlyOneFuelIsUsed(substance, efficiency));
        ssifm.setSubstance(substance);
        logger.info("Setting fuel consumption for {} to {}", ssifm.getSubstance().getName(), ssifm.getShare());

        return fuelMix;
    } else {

        double minimumFuelMixQuality = plant.getTechnology().getMinimumFuelQuality();

        double[] fuelAndCO2Costs = new double[numberOfFuels];
        double[] fuelDensities = new double[numberOfFuels];
        double[] fuelQuality = new double[numberOfFuels];

        int i = 0;
        for (Substance substance : substancePriceMap.keySet()) {
            fuelAndCO2Costs[i] = substancePriceMap.get(substance) + substance.getCo2Density() * (co2Price);
            fuelDensities[i] = substance.getEnergyDensity();
            fuelQuality[i] = substance.getQuality() - minimumFuelMixQuality;
            i++;
        }

        logger.info("Fuel prices: {}", fuelAndCO2Costs);
        logger.info("Fuel densities: {}", fuelDensities);
        logger.info("Fuel purities: {}", fuelQuality);

        // Objective function = minimize fuel cost (fuel
        // consumption*fuelprices
        // + CO2 intensity*co2 price/tax)
        LinearObjectiveFunction function = new LinearObjectiveFunction(fuelAndCO2Costs, 0d);

        List<LinearConstraint> constraints = new ArrayList<LinearConstraint>();

        // Constraint 1: total fuel density * fuel consumption should match
        // required energy input
        constraints.add(new LinearConstraint(fuelDensities, Relationship.EQ, (1 / efficiency)));

        // Constraint 2&3: minimum fuel quality (times fuel consumption)
        // required
        // The equation is derived from (example for 2 fuels): q1 * x1 / (x1+x2) + q2 * x2 / (x1+x2) >= qmin
        // so that the fuelquality weighted by the mass percentages is greater than the minimum fuel quality.
        constraints.add(new LinearConstraint(fuelQuality, Relationship.GEQ, 0));

        try {
            SimplexSolver solver = new SimplexSolver();
            RealPointValuePair solution = solver.optimize(function, constraints, GoalType.MINIMIZE, true);

            logger.info("Succesfully solved a linear optimization for fuel mix");

            int f = 0;
            Iterator<SubstanceShareInFuelMix> iterator = plant.getFuelMix().iterator();
            for (Substance substance : substancePriceMap.keySet()) {
                double share = solution.getPoint()[f];

                SubstanceShareInFuelMix ssifm;
                if (iterator.hasNext()) {
                    ssifm = iterator.next();
                } else {
                    ssifm = new SubstanceShareInFuelMix().persist();
                    fuelMix.add(ssifm);
                }

                double fuelConsumptionPerMWhElectricityProduced = convertFuelShareToMassVolume(share);
                logger.info("Setting fuel consumption for {} to {}", substance.getName(),
                        fuelConsumptionPerMWhElectricityProduced);
                ssifm.setShare(fuelConsumptionPerMWhElectricityProduced);
                ssifm.setSubstance(substance);
                f++;
            }

            logger.info("If single fired, it would have been: {}",
                    calculateFuelConsumptionWhenOnlyOneFuelIsUsed(substancePriceMap.keySet().iterator().next(),
                            efficiency));
            return fuelMix;
        } catch (OptimizationException e) {
            logger.warn(
                    "Failed to determine the correct fuel mix. Adding only fuel number 1 in fuel mix out of {} substances and minimum quality of {}",
                    substancePriceMap.size(), minimumFuelMixQuality);
            logger.info("The fuel added is: {}", substancePriceMap.keySet().iterator().next().getName());

            // Override the old one
            fuelMix = new HashSet<SubstanceShareInFuelMix>();
            SubstanceShareInFuelMix ssifm = new SubstanceShareInFuelMix().persist();
            Substance substance = substancePriceMap.keySet().iterator().next();

            ssifm.setShare(calculateFuelConsumptionWhenOnlyOneFuelIsUsed(substance, efficiency));
            ssifm.setSubstance(substance);
            logger.info("Setting fuel consumption for {} to {}", ssifm.getSubstance().getName(),
                    ssifm.getShare());
            fuelMix.add(ssifm);
            return fuelMix;
        }
    }
}

From source file:org.rascalmpl.library.analysis.linearprogramming.LinearProgramming.java

public IValue llOptimize(IBool minimize, IBool nonNegative, ISet constraints, IConstructor f) {
    SimplexSolver solver = new SimplexSolver();
    ArrayList<LinearConstraint> constraintsJ = new ArrayList<LinearConstraint>(constraints.size());
    for (IValue v : constraints) {
        constraintsJ.add(convertConstraint((IConstructor) v));
    }//from w ww  .j  a  va2s.co m
    LinearObjectiveFunction fJ = convertLinObjFun(f);
    GoalType goal = minimize.getValue() ? GoalType.MINIMIZE : GoalType.MAXIMIZE;
    IValueFactory vf = values;
    boolean nonNegativeJ = nonNegative.getValue();
    try {
        RealPointValuePair res = solver.optimize(fJ, constraintsJ, goal, nonNegativeJ);
        return vf.constructor(Maybe.Maybe_just, vf.constructor(LLSolution_llSolution,
                convertToRealList(res.getPoint(), vf), vf.real(res.getValue())));
    } catch (Exception e) {
        return vf.constructor(Maybe.Maybe_nothing);
    }

}

From source file:rb.app.RBnSCMCSystem.java

public double get_SCM_LB() {
    //return 0.01;

    double min_J_obj = 0.;
    double[] min_Jlocal_obj = new double[n_mubar];

    // Sort the indices of mu_bar based on distance from current_parameters
    List<Integer> sortedmubarIndices = getSorted_CJ_Indices(mu_bar);
    int icount = 0;
    //while ((min_J_obj<=0) && (icount < sortedmubarIndices.size())){
    while ((min_J_obj <= 0) && (icount < sortedmubarIndices.size())) {
        int imubar = sortedmubarIndices.get(icount);

        // First, declare the constraints
        Collection constraints = new ArrayList();

        // Add bounding box constraints for the get_Q_a() variables
        for (int q = 0; q < get_Q_a(); q++) {
            double[] index = new double[get_Q_a() * 2];
            index[q] = 1.;/*from  w  w  w.  j  ava 2  s . c o  m*/

            constraints.add(new LinearConstraint(index, Relationship.GEQ, B_min[q] / beta_bar[imubar]));
            constraints.add(new LinearConstraint(index, Relationship.LEQ, B_max[q] / beta_bar[imubar]));

            index[q] = 0.;
            index[q + get_Q_a()] = 1.;

            constraints.add(new LinearConstraint(index, Relationship.GEQ, B_min[q] / beta_bar[imubar]));
            constraints.add(new LinearConstraint(index, Relationship.LEQ, B_max[q] / beta_bar[imubar]));
        }

        // Save the current_parameters since we'll change them in the loop below
        save_current_parameters();

        // Add the constraint rows
        if (n_muhat[imubar] > 0) {
            for (int imuhat = 0; imuhat < n_muhat[imubar]; imuhat++) {
                current_parameters = mu_hat[imubar].get(imuhat);

                double[] constraint_row = new double[get_Q_a() * 2];
                for (int q = 0; q < get_Q_a(); q++) {
                    Complex theta_q_a = complex_eval_theta_q_a(q);
                    constraint_row[q] = theta_q_a.getReal() * beta_bar[imubar];
                    constraint_row[q + get_Q_a()] = theta_q_a.getImaginary() * beta_bar[imubar];
                }

                constraints
                        .add(new LinearConstraint(constraint_row, Relationship.GEQ, beta_hat[imubar][imuhat]));
            }
        }

        // Now load the original parameters back into current_parameters
        // in order to set the coefficients of the objective function
        reload_current_parameters();

        // Create objective function object
        double[] objectiveFn = new double[get_Q_a() * 2];
        for (int q = 0; q < get_Q_a(); q++) {
            Complex theta_q_a = complex_eval_theta_q_a(q);
            objectiveFn[q] = theta_q_a.getReal() * beta_bar[imubar];
            objectiveFn[q + get_Q_a()] = theta_q_a.getImaginary() * beta_bar[imubar];
        }
        LinearObjectiveFunction f = new LinearObjectiveFunction(objectiveFn, 0.);

        try {
            SimplexSolver solver = new SimplexSolver(1e-6);
            RealPointValuePair opt_pair = solver.optimize(f, constraints, GoalType.MINIMIZE, false);
            min_Jlocal_obj[icount] = opt_pair.getValue();
        } catch (OptimizationException e) {
            Log.e("RBSCMSYSTEM_TAG", "Optimal solution not found");
            e.printStackTrace();
        } catch (Exception e) {
            Log.e("RBSCMSYSTEM_TAG", "Exception occurred during SCM_LB calculation");
            e.printStackTrace();
        }

        min_J_obj = min_J_obj > min_Jlocal_obj[icount] ? min_J_obj : min_Jlocal_obj[icount];
        icount++;
    }
    return min_J_obj;
}

From source file:rb.app.RBSCMSystem.java

/**
* @return the SCM lower bound for the current parameters.
*///from www.ja  va  2  s  . c o m
public double get_SCM_LB() {
    double min_J_obj = 0.;

    try {

        // First, declare the constraints
        Collection constraints = new ArrayList();

        // Add bounding box constraints for the get_Q_a() variables
        for (int q = 0; q < get_Q_a(); q++) {
            double[] index = new double[get_Q_a()];
            index[q] = 1.;

            constraints.add(new LinearConstraint(index, Relationship.GEQ, B_min[q]));
            constraints.add(new LinearConstraint(index, Relationship.LEQ, B_max[q]));
        }

        // Sort the indices of C_J based on distance from current_parameters
        List<Integer> sortedIndices = getSorted_CJ_Indices();

        // Save the current_parameters since we'll change them in the loop below
        save_current_parameters();

        // Add the constraint rows
        int n_rows = Math.min(SCM_M, C_J.size());
        int count = 1;

        if (n_rows > 0) {
            for (Iterator it = sortedIndices.iterator(); it.hasNext();) {
                Integer mu_index = (Integer) it.next();

                get_current_parameters_from_C_J(mu_index);
                // current_parameters = C_J.get(mu_index);

                double[] constraint_row = new double[get_Q_a()];
                for (int q = 0; q < get_Q_a(); q++) {
                    constraint_row[q] = eval_theta_q_a(q);
                }

                constraints.add(
                        new LinearConstraint(constraint_row, Relationship.GEQ, C_J_stability_vector[mu_index]));

                if (count >= n_rows)
                    break;

                count++;

            }
        }

        // Now load the original parameters back into current_parameters
        // in order to set the coefficients of the objective function
        reload_current_parameters();

        // Create objective function object
        double[] objectiveFn = new double[get_Q_a()];
        for (int q = 0; q < get_Q_a(); q++) {
            objectiveFn[q] = eval_theta_q_a(q);
        }
        LinearObjectiveFunction f = new LinearObjectiveFunction(objectiveFn, 0.);

        SimplexSolver solver = new SimplexSolver();
        RealPointValuePair opt_pair = solver.optimize(f, constraints, GoalType.MINIMIZE, false);
        min_J_obj = opt_pair.getValue();
    } catch (OptimizationException e) {
        Log.e("DEBUG_TAG", "Optimal solution not found");
        e.printStackTrace();
    } catch (Exception e) {
        Log.e("DEBUG_TAG", "Exception occurred during SCM_LB calculation");
        e.printStackTrace();
    }

    Log.d(DEBUG_TAG, "SCM val = " + min_J_obj);
    return min_J_obj;
}