英文:
How to get multipliers after solving a quadratic program in ojAlgo
问题
我实现了一个顺序二次规划（Sequential quadratic programming，SQP）优化器，并使用ojAlgo来解决二次规划（Quadratic programming，QP）子问题。
我的问题是：
如何获得QP解的"拉格朗日乘数"（Lagrange multipliers）？
在附加的示例代码中，对解决QP问题的代码 result.getMultipliers() 仅返回一个空的Optional。
// ...（之前的代码部分）
// How do I get the multipliers
Optional<Access1D<?>> multipliers = result.getMultipliers();
double value1 = result.getValue();
// ...（之后的代码部分）
更新 1:
这是我使用 org.ojalgo.optimisation.convex.ConvexSolver.getBuilder() 进行重新编写的示例：
// ...（之前的代码部分）
// 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);
// ...（之后的代码部分）
注意：我只返回翻译好的部分，代码的其余部分被省略。
英文:
I implement a Sequential quadratic programming (SQP) optimizer and use ojAlgo for the quadratic programming (QP) subproblem.
My question is:
How do I get hold of the "Lagrange multipliers" for the QP solution?
In the attached example code that solve an QP result.getMultipliers() only return an empty 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, 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);
}
}
Update 1:
Here is my reworked example using 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', 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);
}
}
答案1
得分: 0
我认为这是一个Optional，因为有时将求解器中的Lagrange乘数映射到模型的约束条件中会变得混乱不堪。
如果您正在实现一个SQP求解器，我建议您不要以ExpressionsBasedModel为基础进行实现，而是直接委托给凸优化求解器。构建一个实现org.ojalgo.optimisation.Optimisation.Solver接口的内容，并委托给org.ojalgo.optimisation.convex包中的各种类。然后，您可以直接使用矩阵、向量和乘数进行编码。
为了使该求解器可以被ExpressionsBasedModel使用，您还需要实现一个org.ojalgo.optimisation.Optimisation.Integration接口，并通过调用ExpressionsBasedModel.addPreferredSolver(myIntegration)或ExpressionsBasedModel.addFallbackSolver(myIntegration)将其注册到模型中。
实现求解器和使其能够从建模工具中使用是两个不同的事情。
英文:
I believe that is an Optional because it was (sometimes) too messy to map the Lagrange multipliers from the solver to the constraints of the model.
If you're implementing an SQP solver may I suggest that you don't implement it in terms of ExpressionsBasedModel, but delegate to the convex solvers directly. Build something that implements org.ojalgo.optimisation.Optimisation.Solver and delegate to the various classes in the org.ojalgo.optimisation.convex package. Then you code more directly with the matrices, vectors and multipliers.
To make that solver usable by ExpressionsBasedModel you also implement an org.ojalgo.optimisation.Optimisation.Integration and register that by calling ExpressionsBasedModel.addPreferredSolver(myIntegeration) or ExpressionsBasedModel.addFallbackSolver(myIntegeration).
Implementing a solver and making it usable from the modelling tool are two separate things.
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论