This forum has been archived. All content is frozen. Please use KDE Discuss instead.

linear least squares with inequality constraint

Tags: None
(comma "," separated)
User avatar
ggael
Moderator
Posts
3447
Karma
19
OS
You are passing evX instead of evB to solve_quadprog.
euans
Registered Member
Posts
3
Karma
0
It looks like I'd made a mistake setting up the constraints. These should be:

CI = Eigen::MatrixXd::Zero(nVar, 2 * 2);
ci0 = Eigen::VectorXd::Zero(2*2);
CI(0, 0) = -1.;
CI(1, 1) = -1.;
CI(0, 2) = 1.;
CI(1, 3) = 1.;
ci0(0) = 1.;
ci0(1) = 1.;
ci0(2) = 1.;
ci0(3) = 1.;

solve_quadprog now returns the correct result when passed this instead:

solve_quadprog(emA, evB, CE, ce0, CI, ci0, evX);

Last edited by euans on Thu Aug 04, 2016 1:45 pm, edited 1 time in total.
euans
Registered Member
Posts
3
Karma
0
Thanks for having a look through this code.

I did also make the mistake of passing solve_quadprog evX instead of evB. Correcting this mistake caused the original code to return the result -1, -1.
kleintheresa
Registered Member
Posts
1
Karma
0
I'm having a similar problem implementing a least squares solution to a problem.
Incidentally, I created a C++ template version of eiquadprog.hpp, which allows you to specify the size of the matrices at compile time, which avoid dynamic memory allocations.
But maybe the bug is in my modified code?

Here is the templated code:
http://s000.tinyupload.com/index.php?fi ... 7920772892

So I have no equality constraints, and inequality between 0 and 440.
The problem I'm trying to find is the thrust firings to produce a desired force and torque of a spacecraft in zero gravity.

So I have X as a Nx6 matrix where the first 3 elements in each row is the position crossed with the unit direction vector (the torque produced by a thrust of unit 1), the second 3 elements are the direction vector (the acceleration), and there are N rows for N thrusters.
For example

X =[ 0 -0.0200 0 1.0000 0 0;
0 0.3361 0 1.0000 0 0;
0 -0.3761 0 1.0000 0 0;
0 -0.0200 -0.3561 1.0000 0 0;
0 -0.0200 0.3561 1.0000 0 0];

y is a desired net torque and thrust I want to achieve, say:
y = [0, 0, 0, 2123.6559802463221, -2.8158062128244268, -249.58517518020722];

So I'm attempting to solve X*F = y .
But the thrusters obviously cannot fire negatively or above their maximum thrust of 440.
so I have F >= 0 and F<= 440 for all N thrusters.

For some reason I could not get my templated code to take 0 equality constraints so I set up one "dummy" equality constraint, basically F*0 = 0.
// equality constraints
CE2 = Eigen::Matrix<double, N, 1>::Zero();
ce02 = Eigen::Matrix<double, 1, 1>::Zero();


// inequality constraints
CI2 << Eigen::Matrix<double, N, N::Identity(), -Eigen::Matrix<double, N, N>::Identity();
ci02 << Eigen::Matrix<double, N, 1>::Zero(), 440*Eigen::Matrix<double, N, 1>::Ones();

And then I have, as before,
G = X.transpose()*X;
g0 = X.transpose()*y;

double cost = Eigen::solve_quadprog<N, 1, N*2>(G, g0 , CE, ce0, CI, ci0, F);

The problem is that this seems to be giving me values outside the constraints. it is giving me F = [2004.3896666303358, 119.26631361598683, -3.4901494507931849e-13, 0, 0]

But that is greater than my constraint of 440N! Also my cost is -2254957.3612179835, which seems like a strange number. Cost should be positive right?
What's going on.

For the record, my G matrix is:
G =[ 1.0004 0.9933 1.0075 1.0004 1.0004
0.9933 1.1130 0.8736 0.9933 0.9933
1.0075 0.8736 1.1415 1.0075 1.0075
1.0004 0.9933 1.0075 1.1272 0.8736
1.0004 0.9933 1.0075 0.8736 1.1272]

and g0 = [-2123.6559802463221, -2123.6559802463221, -2123.6559802463221, -2123.6559802463221, -2123.6559802463221]

Matlab is telling me it is close to singular or badly scaled. How can i fix this problem? Is there some way to scale X so that G is not singular?

Side, note, I did have some code working where I was solving it like this:
G = Eigen::Matrix<double, N,N>::Identity();
g0 = Eigen::Matrix<double, N,1>::Zero();

double cost = Eigen::solve_quadprog<N, 6, N*2>(G, g0, X, y, CI, ci0, F);

basically setting X*F = y as an equality constraint and then minimizing the RSS of F. (y is sign flipped in this version).
That seemed to work ok as long as X*F=y is possible within the constraints, but it breaks down when the constraint is not achievable. I would like a solution which gets me as close as possible to the desired force and torque, hence the least squares problem.
User avatar
ggael
Moderator
Posts
3447
Karma
19
OS
For the record, I have an updated version (almost completely rewrote), there:

https://gitlab.inria.fr/alta/alta/tree/ ... quadprog++

might be worth trying first.
User avatar
ggael
Moderator
Posts
3447
Karma
19
OS
Looking at your numbers, your problem is really ill-posed, it has two zero singular values meaning a 2D subspace of solution. Therefore, QP++ cannot perform the LLT decomposition of X*X'.


Bookmarks



Who is online

Registered users: Bing [Bot], daret, Google [Bot], sandyvee, Sogou [Bot]