Visual Studio Code Question posted in
View the official documentation.

Visual Studio Code – Discrepancy in qpOASES Optimal Solutions Between C++ and MATLAB Implementations

I have implemented a quadratic programming problem using qpOASES in MATLAB and C++. The optimal solutions differ despite using the same matrices and vectors in both environments. The MATLAB implementation produces valid solutions (as of my best guess), confirmed by the Webots simulation interface. However, the C++ implementation, when compiled and executed, does not provide a correct solution, leading to performance issues in the Webots environment. Here I have provided a simple yet similar example related to my problem. Upon testing on both platforms, still got different outputs.

MATLAB Code (R2022a)

% Define the matrix H (12 by 12)
H = [1.7201, 0.981765, 0, -1.70402, 0.981765, 0, 1.7161, -1.15451, 0, -1.70402, -1.15451, 0;
     0.981765, 0.573682, 0, -0.981765, 0.569682, 0, 0.981765, -0.656776, 0, -0.981765, -0.656776, 0;
     0, 0, 0.0153304, 0, 0, 0.000762503, 0, 0, 0.00687851, 0, 0, -0.00368936;
     -1.70402, -0.981765, 0, 1.7201, -0.981765, 0, -1.70402, 1.15451, 0, 1.7161, 1.15451, 0;
     0.981765, 0.569682, 0, -0.981765, 0.573682, 0, 0.981765, -0.656776, 0, -0.981765, -0.656776, 0;
     0, 0, 0.000762503, 0, 0, 0.0135145, 0, 0, -0.00171378, 0, 0, 0.00703826;
     1.7161, 0.981765, 0, -1.70402, 0.981765, 0, 1.7201, -1.15451, 0, -1.70402, -1.15451, 0;
     -1.15451, -0.656776, 0, 1.15451, -0.656776, 0, -1.15451, 0.789478, 0, 1.15451, 0.785478, 0;
     0, 0, 0.00687851, 0, 0, -0.00171378, 0, 0, 0.0139643, 0, 0, 0.00137201;
     -1.70402, -0.981765, 0, 1.7161, -0.981765, 0, -1.70402, 1.15451, 0, 1.7201, 1.15451, 0;
     -1.15451, -0.656776, 0, 1.15451, -0.656776, 0, -1.15451, 0.785478, 0, 1.15451, 0.789478, 0;
     0, 0, -0.00368936, 0, 0, 0.00703826, 0, 0, 0.00137201, 0, 0, 0.0160996];

% Define the the vector g (12 by 1)
g = [117.537; 67.4661; -0.393848; -117.52; 67.4661; -0.39092; 117.537; -79.355; -0.393697; -117.52; -79.355; -0.39077];

% Define Matrix A (20 by 12)
A = [0.666667, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0;
    -0.666667, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0;
    0, 0.666667, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0;
    0, -0.666667, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0;
    0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0.666667, 0, 1, 0, 0, 0, 0, 0, 0;
    0, 0, 0, -0.666667, 0, 1, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0.666667, 1, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, -0.666667, 1, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0.666667, 0, 1, 0, 0, 0;
    0, 0, 0, 0, 0, 0, -0.666667, 0, 1, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 0.666667, 1, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, -0.666667, 1, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0;
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0.666667, 0, 1;
    0, 0, 0, 0, 0, 0, 0, 0, 0, -0.666667, 0, 1;
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.666667, 1;
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.666667, 1;
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1];

% Define the constraint vector lower bounds for A matrix lbA (20 by 1)
lbA = zeros(20, 1);

% Define the constraint vector upper bounds for A matrix ubA (20 by 1)
ubA = [150; 150; 150; 150; 120; 150; 150; 150; 150; 120; 150; 150; 150; 150; 120; 150; 150; 150; 150; 120];

% Solve the QP problem using qpOASES
[x,fval,exitflag,iter,lambda] = qpOASES(H, g, A, [], [], lbA, ubA, []);
% Output the result
disp('Optimal solution:');
disp(x);

**C++ Code (Microsoft Visual Studio 2019) **

#include <iostream>
#include <Eigen/Dense>
#include <qpOASES.hpp>

using namespace Eigen;
using namespace qpOASES;

int main()
{
    USING_NAMESPACE_QPOASES;

    // Initialize the qpOASES library
    QProblem qp(12, 20); // 12 variables, 20 constraints

    // Set up the matrices and vectors as shown in the images
    MatrixXd H(12, 12);
    H << 1.7201, 0.981765, 0, -1.70402, 0.981765, 0, 1.7161, -1.15451, 0, -1.70402, -1.15451, 0,
        0.981765, 0.573682, 0, -0.981765, 0.569682, 0, 0.981765, -0.656776, 0, -0.981765, -0.656776, 0,
        0, 0, 0.0153304, 0, 0, 0.000762503, 0, 0, 0.00687851, 0, 0, -0.00368936,
        -1.70402, -0.981765, 0, 1.7201, -0.981765, 0, -1.70402, 1.15451, 0, 1.7161, 1.15451, 0,
        0.981765, 0.569682, 0, -0.981765, 0.573682, 0, 0.981765, -0.656776, 0, -0.981765, -0.656776, 0,
        0, 0, 0.000762503, 0, 0, 0.0135145, 0, 0, -0.00171378, 0, 0, 0.00703826,
        1.7161, 0.981765, 0, -1.70402, 0.981765, 0, 1.7201, -1.15451, 0, -1.70402, -1.15451, 0,
        -1.15451, -0.656776, 0, 1.15451, -0.656776, 0, -1.15451, 0.789478, 0, 1.15451, 0.785478, 0,
        0, 0, 0.00687851, 0, 0, -0.00171378, 0, 0, 0.0139643, 0, 0, 0.00137201,
        -1.70402, -0.981765, 0, 1.7161, -0.981765, 0, -1.70402, 1.15451, 0, 1.7201, 1.15451, 0,
        -1.15451, -0.656776, 0, 1.15451, -0.656776, 0, -1.15451, 0.785478, 0, 1.15451, 0.789478, 0,
        0, 0, -0.00368936, 0, 0, 0.00703826, 0, 0, 0.00137201, 0, 0, 0.0160996;

    VectorXd g(12);
    g << 117.537, 67.4661, -0.393848, -117.52, 67.4661, -0.39092, 117.537, -79.355, -0.393697, -117.52, -79.355, -0.39077;

    MatrixXd A(20, 12);
    A << 0.666667, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        -0.666667, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0.666667, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, -0.666667, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0.666667, 0, 1, 0, 0, 0, 0, 0, 0,
        0, 0, 0, -0.666667, 0, 1, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0.666667, 1, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, -0.666667, 1, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0.666667, 0, 1, 0, 0, 0,
        0, 0, 0, 0, 0, 0, -0.666667, 0, 1, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0.666667, 1, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, -0.666667, 1, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0.666667, 0, 1,
        0, 0, 0, 0, 0, 0, 0, 0, 0, -0.666667, 0, 1,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.666667, 1,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.666667, 1,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1;

    VectorXd lbA = VectorXd::Zero(20);
    VectorXd ubA(20);
    ubA << 150, 150, 150, 150, 120, 150, 150, 150, 150, 120, 150, 150, 150, 150, 120, 150, 150, 150, 150, 120;

    // Set the number of working set recalculations (iterations)
    int nWSR = 100;

    // Solve the QP problem
    returnValue status = qp.init(H.data(), g.data(), A.data(), nullptr, nullptr, lbA.data(), ubA.data(), nWSR);

    // Retrieve the solution
    VectorXd x(12);
    qp.getPrimalSolution(x.data());

    

    // Check if the problem was solved successfully
    if (status == SUCCESSFUL_RETURN)
    {
        qp.getPrimalSolution(x.data());
        std::cout << "Successful: Optimal solution: " << std::endl;
    }
    else
    {
        std::cout << "Failed to solve the QP problem. Status: " << status << std::endl;
    }

    // Output the result
    std::cout << "Optimal solution:" << std::endl << x << std::endl;

    return 0;
}

Output in MATLAB and C++ after qpOASES optimization

The optimal solution vector obtained after the qpOASES operation in both environments is:
MATLAB, C++ output

Given that both implementations use the same matrices and vectors (H, g, A, lbA ubA), I am confused by the discrepancy in the solutions. Could there be an issue with the configuration or usage of qpOASES in the C++ code? Any insights or suggestions to resolve this issue would be greatly appreciated.

Tags:

2

Answers


  1. Chosen as BEST ANSWER
    0

    Thank you for your answer. I understand that by default, Eigen stores matrices in column-major order, similar to MATLAB. However, I noticed that for qpOASES, I need to explicitly define the constraint matrix A in row-major order using Eigen's RowMajor specifier to get the correct results. But I have 2 questions:

    1. Why do I need to define matrix A in row-major order in C++ when using qpOASES, but I don't need to do this in MATLAB? If MATLAB also stores matrices in column-major order, why does qpOASES work correctly without explicitly specifying the row-major format?
    2. Why does the Hessian matrix H produce identical results regardless of whether it's stored in a row-major or column-major format in Eigen, while matrix A does not?

  2. 0

    This is a matrix data conversion problem with Eigen. If you store your constraint matrix data like

    double A_data[12*20] = {
    0.666667, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    -0.666667, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0.666667, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, -0.666667, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0.666667, 0, 1, 0, 0, 0, 0, 0, 0,
    0, 0, 0, -0.666667, 0, 1, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0.666667, 1, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, -0.666667, 1, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0.666667, 0, 1, 0, 0, 0,
    0, 0, 0, 0, 0, 0, -0.666667, 0, 1, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0.666667, 1, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, -0.666667, 1, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0.666667, 0, 1,
    0, 0, 0, 0, 0, 0, 0, 0, 0, -0.666667, 0, 1,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.666667, 1,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.666667, 1,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1};

    and pass A_data instead of A.data() to qp.init, you get what you expect. The A.data() you are using is in column major format, i.e., the array starts with 0.666667, -0.666667, 0, …

    The problem does not appear for the Hessian because H is symmetric.

    Alternatively, you can define your constraint matrix as

    Eigen::Matrix<double, Dynamic, Dynamic, RowMajor> A(20, 12);

    (see https://eigen.tuxfamily.org/dox/group__TopicStorageOrders.html).

    Login or Signup to reply.
Please signup or login to give your own answer.

Back To Top
Search