在ojAlgo中求解二次程序后如何获得乘法器

问题描述

我实现了一个顺序二次编程（SQP）优化器，并使用ojAlgo来解决二次编程（QP）子问题。

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

在随附的示例代码中，它解决了QP结果。getMultipliers（）仅返回一个空的Optional。

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,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,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,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);
    }
   
}

更新1： 这是我使用org.ojalgo.optimisation.convex.ConvexSolver.getBuilder（）;重做的示例。

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',x2 + x3 = 0
//  A*x = b

//objectiveGradient
      Primitive64Store gStore = Primitive64Store.FACTORY.rows(new double[][]{
         {-8},{-3}
      });
//objectiveHessian
      Primitive64Store HStore = Primitive64Store.FACTORY.rows(new double[][]{
         {6,4}
      });
// constraint equations
      Primitive64Store AStore = Primitive64Store.FACTORY.rows(new double[][]{
         {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);

   }

}

解决方法

我认为这是一个Optional，因为（有时）太乱了，无法将拉格朗日乘子从求解器映射到模型的约束条件。

如果您要实现SQP求解器，我建议您不要根据ExpressionsBasedModel来实现它，而应直接委托给凸求解器。构建实现org.ojalgo.optimisation.Optimisation.Solver的东西，并将其委托给org.ojalgo.optimisation.convex包中的各个类。然后，您可以直接使用矩阵，向量和乘数进行编码。

要使该求解器在ExpressionsBasedModel中可用，您还可以实现org.ojalgo.optimisation.Optimisation.Integration并通过调用ExpressionsBasedModel.addPreferredSolver(myIntegeration)或ExpressionsBasedModel.addFallbackSolver(myIntegeration)进行注册。

实现求解器并使其在建模工具中可用是两件分开的事情。