如何在ojAlgo中求解二次规划后求乘子

Question

我实现了一个序列二次规划(SQP)优化器，并使用ojAlgo解决二次规划(QP)子问题。

我的问题是:如何获得QP解的“拉格朗日乘数”？

在附带的示例代码中，解决QP的result.getMultipliers()只返回一个空的可选项。

package com.mycompany.testojalgo;

import java.math.BigDecimal;
import java.util.ArrayList;
import java.util.List;
import java.util.Optional;
import org.ojalgo.matrix.Primitive64Matrix;
import org.ojalgo.optimisation.Expression;
import org.ojalgo.optimisation.ExpressionsBasedModel;
import org.ojalgo.optimisation.Optimisation;
import org.ojalgo.optimisation.Variable;
import org.ojalgo.structure.Access1D;
import org.ojalgo.type.StandardType;
import org.ojalgo.type.context.NumberContext;
public class ojAlgoQP {
   
   public static void main(String[] args) {
      testOjAlgoQuadraticProgramming();
   }

    public static void testOjAlgoQuadraticProgramming() {
//  QP Example 16.2 p453 in 'Numerical Optimization', 2ed, (2006), Jorge Nocedal and Stephen J. Wright.
//  minimize function F(x1,x2,x3) = 3*x1*x1 + 2*x1*x2 + x1*x3 + 2.5*x2*x2 + 2*x2*x3 + 2*x3*x3 - 8*x1 - 3*x2 - 3*x3
//  x = [x1, x2, x3]'
//  F(x) = 1/2*x'*H*x + x'*g
//  constraints x1 + x3 = 3, x2 + x3 = 0
//  A*x = b

//objectiveGradient
        Primitive64Matrix g = Primitive64Matrix.FACTORY.rows(new double[][]{
            {-8}, {-3}, {-3}
        });
//objectiveHessian
        Primitive64Matrix H = Primitive64Matrix.FACTORY.rows(new double[][]{
            {6, 2, 1},
            {2, 5, 2},
            {1, 2, 4}
        });

        Variable x1 = new Variable("x1");
        Variable x2 = new Variable("x2");
        Variable x3 = new Variable("x3");
// constraint equations
        Primitive64Matrix A = Primitive64Matrix.FACTORY.rows(new double[][]{
            {1, 0, 1},
            {0, 1, 1}
        });
// required constraint values
        Primitive64Matrix b = Primitive64Matrix.FACTORY.rows(new double[][]{
            {3}, {0}
        });
        
        List<Variable> variables = new ArrayList<>();
        variables.add(x1);
        variables.add(x2);
        variables.add(x3);
        
        ExpressionsBasedModel model = new ExpressionsBasedModel(variables);  
        
        Expression energy = model.addExpression("Energy");
        energy.setLinearFactors(variables, g);
//divide by two to express function using hessian.        
        energy.setQuadraticFactors(variables, H.divide(2));
        energy.weight(BigDecimal.ONE);
        
//create constraint equations
        for (int i = 0; i < A.countRows(); i++) {
            Expression expression = model.addExpression("Constraint#"+i);
            for (int j = 0; j < A.countColumns(); j++) {
                expression.set(variables.get(j), A.get(i, j));
            }
            expression.level(b.get(i));
        }
        
        Optimisation.Result result = model.minimise();
        
        NumberContext accuracy = StandardType.PERCENT.withPrecision(1);
        boolean ok = model.validate(result, accuracy);        
        Optimisation.State v = result.getState();
        
// How do I get the multipliers
        Optional<Access1D<?>> multipliers = result.getMultipliers();
        double value1 = result.getValue();         

// Get result and check value and constraint
        Primitive64Matrix x = Primitive64Matrix.FACTORY.rows(new double[][]{
            {x1.getValue().doubleValue()}, {x2.getValue().doubleValue()}, {x3.getValue().doubleValue()}
        });
//divide by two to express function using hessian, again.  
        Primitive64Matrix value = x.transpose().multiply(H.divide(2)).multiply(x).add(x.transpose().multiply(g));
        Primitive64Matrix residual= A.multiply(x).subtract(b);
    }
   
}

org.ojalgo.optimisation.convex.ConvexSolver.getBuilder()；更新1：更新1:这里是我使用重新编写的示例

package com.mycompany.testojalgo;

import java.util.Optional;
import org.ojalgo.matrix.store.MatrixStore;
import org.ojalgo.matrix.store.Primitive64Store;
import org.ojalgo.optimisation.Optimisation;
import org.ojalgo.optimisation.convex.ConvexSolver;
import org.ojalgo.structure.Access1D;

public class ojAlgoQP {

   public static void main(String[] args) {
      testOjAlgoQuadraticProgramming2();
   }

   public static void testOjAlgoQuadraticProgramming2() {
//  QP Example 16.2 p453 in 'Numerical Optimization', 2ed, (2006), Jorge Nocedal and Stephen J. Wright.
//  minimize function F(x1,x2,x3) = 3*x1*x1 + 2*x1*x2 + x1*x3 + 2.5*x2*x2 + 2*x2*x3 + 2*x3*x3 - 8*x1 - 3*x2 - 3*x3
//  x = [x1, x2, x3]'
//  F(x) = 1/2*x'*H*x + x'*g
//  constraints x1 + x3 = 3, x2 + x3 = 0
//  A*x = b

//objectiveGradient
      Primitive64Store gStore = Primitive64Store.FACTORY.rows(new double[][]{
         {-8}, {-3}, {-3}
      });
//objectiveHessian
      Primitive64Store HStore = Primitive64Store.FACTORY.rows(new double[][]{
         {6, 2, 1},
         {2, 5, 2},
         {1, 2, 4}
      });
// constraint equations
      Primitive64Store AStore = Primitive64Store.FACTORY.rows(new double[][]{
         {1, 0, 1},
         {0, 1, 1}
      });
// required constraint values
      Primitive64Store bStore = Primitive64Store.FACTORY.rows(new double[][]{
         {3}, {0}
      });
      ConvexSolver.Builder builder = ConvexSolver.getBuilder();
      builder.equalities(AStore, bStore);
      builder.objective(HStore, gStore.negate());
      ConvexSolver solver = builder.build();
      Optimisation.Result result = solver.solve();

// How do I get the multipliers?  multipliers = Optional.empty
      Optional<Access1D<?>> multipliers = result.getMultipliers();
// value1 = -3.5
      double value1 = result.getValue();

// Verify result:
// x= [2.0, -0.9999999999999996, 0.9999999999999997]';
// value = -3.5
// residual =[-4.440892098500626E-16, 1.1102230246251565E-16]'
      Primitive64Store x = Primitive64Store.FACTORY.column(result.toRawCopy1D());
      MatrixStore<Double> value = x.transpose().multiply(HStore.multiply(0.5)).multiply(x).add(x.transpose().multiply(gStore));
      MatrixStore<Double> residual = AStore.multiply(x).subtract(bStore);

   }

}

Answer 1

我认为这是一个Optional，因为它(有时)过于混乱，无法将拉格朗日乘数从求解器映射到模型的约束。

如果您正在实现一个SQP求解器，我建议您不要用ExpressionsBasedModel来实现它，而是直接委托凸求解器。构建实现org.ojalgo.optimisation.Optimisation.Solver并将其委托给org.ojalgo.optimisation.convex包中的各个类的东西。然后，用矩阵、向量和乘子更直接地编码。

为了使ExpressionsBasedModel可以使用该求解器，您还可以实现一个org.ojalgo.optimisation.Optimisation.Integration，并通过调用ExpressionsBasedModel.addPreferredSolver(myIntegeration)或ExpressionsBasedModel.addFallbackSolver(myIntegeration)注册它。

实现一个求解器并使它从建模工具中可用是两件不同的事情。