MPC with Ax update in python, "primal infeasible" problem.

[black...@gmail.com]

Hello!

I have a question which has already been asked in several other threads, but I need clarification. Similar threads here: https://groups.google.com/forum/#!topic/osqp/VlmMAI2Ti1w
and here:
https://groups.google.com/forum/#!topic/osqp/GhR_71HHzVg

I'm using OSQP in python to solve an MPC problem. If I recreate the OSQP instance and call "setup" at every new timestep, the system is solved (fairly quickly!) correctly. However, if I try to update Ax, q, l, and u (my P matrix stays the same), I get a report of primal infeasibility each time.

I'm using sparse matrices. I set up the initial OSQP instance using a "fully occupied" structure matrix A, with 1's in every location where a nonzero value could potentially be, zeros elsewhere. This matrix has on the order of 400-ish nonzeros. I'm using a csc_matrix from scipy.sparse. I call this structure matrix "A_full", and keep track of the row and column indices like this:

self.A_full = ...
self.A_full.eliminate_zeros()
self.A_full.sort_indices()
(self.A_full_row_indices, self.A_full_col_indices) = self.A_full.nonzero()
self.solver = osqp.OSQP()
self.solver.setup(..., A=self.A_full)

When I construct new A matrices at each new timestep, the sparsity is slightly different (more zeros, in varying places) because of the Jacobians of the changing linearized system. I call this new matrix A_current, and I get an array of values that exactly matches the layout of the A_full matrix by using the indices, like this:

self.A_current = ...
Ax_values = self.A_current[self.A_full_row_indices, self.A_full_col_indices].A1
self.solver.update(..., Ax=Ax_values)

When I do this, I get "primal infeasible" errors. The length of Ax_values is exactly the same as the length of self.A_full.data. I read the matlab example code provided in the threads above by Goran, but it doesn't have the same details of sparse layout as python, so it's impossible to compare.

Given the amount of confusion this seems to cause many many people, it would be incredibly helpful to have an update function that could take a sparse matrix A (and P), and just check that it is a proper sparse subset of the original A and P, and then do the data copying correctly.

What am I doing wrong?

Thanks! Chris

[Paul Goulart]

It is a little bit hard to diagnose without seeing a MRE.    Could you recreate the error with a smallish problem (Say a 4x4 A)?    

Some possible things to consider:

1) Do you have a 0/1 indexing issue?   Probably not if this is a python problem.

2) What happens if you solve the problem and then try to overwrite the data with exactly the same values as the first solve?   Does it fail?   What about if you *don't* update the problem data and then call the solve twice in a row"   If the first version fails and the second doesn't, then you are very likely to be updating the problem date incorrectly.

3) Is the column index array the same between calls?

I agree that there is a lack of documentation around this question, and we should add something to the docs showing how to do this for matlab/simulink/python/julia.

Another useful feature might be to provide a function in all interfaces that dumps the P and A matrix that OSQP is using back to the user on request.   If you are willing, you could make a feature request for these on github.

PG

[Christopher Horvath]

Hey Paul, 

Thanks for your quick reply. You outlined the exact steps I'm currently doing, so at least I'm on the right track. If you had an example script of doing this how *you* would do it, that would be immensely helpful. I started off my whole system by evolving the quadcopter MPC control example in your documentation.

Meanwhile, I'll see if I can extract a self-contained example.

Chris

--
You received this message because you are subscribed to a topic in the Google Groups "OSQP" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/osqp/ZFvblAQdUxQ/unsubscribe.
To unsubscribe from this group and all its topics, send an email to osqp+uns...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/osqp/69bc9611-d7b2-478a-bc55-49dea6a0a78e%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Paul Goulart]

Have you had a look at the docs here:

https://osqp.org/docs/examples/update-matrices.html 

On Monday, June 3, 2019 at 8:31:44 PM UTC+1, Christopher Horvath wrote:

Hey Paul, 

Thanks for your quick reply. You outlined the exact steps I'm currently doing, so at least I'm on the right track. If you had an example script of doing this how *you* would do it, that would be immensely helpful. I started off my whole system by evolving the quadcopter MPC control example in your documentation.

Meanwhile, I'll see if I can extract a self-contained example.

Chris

To unsubscribe from this group and all its topics, send an email to os...@googlegroups.com.

[Christopher Horvath]

The problem with the example is that the original matrix and the new matrix have exactly the same sparse structure, so the data is aligned already. When you have a dynamic system, sometimes zeros appear in the matrix that weren't there before - I'm doing a steering model of a vehicle, and when the orientation goes to zero, some of the fields that were not zero before become zero.

I've got a solution, I'm trying to make it efficient and then I'll share it with this thread.

To unsubscribe from this group and all its topics, send an email to osqp+uns...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/osqp/ce08720f-b486-4281-ac03-57ab0cbf2245%40googlegroups.com.

[Christopher Horvath]

I have a solution to this problem. The problem is that the "nonzero" method of the csc_matrix in the sparse library returns coordinates sorted in a row-major way, but because we are using csc matrixes, the data in the vector is stored in a column-major way.  So you see the problem in this simple example:

B = sparse.csc_matrix([[1, 0, 0, 2], [3, 4, 0, 5], [0, 6, 7, 0], [0, 0, 8, 9]])

print(B.data)

    # [1 3 4 6 7 8 2 5 9]

(row_indices, col_indices) = B.nonzero()

print(B[row_indices, col_indices].A1)

    # [1 2 3 4 5 6 7 8 9]

Unfortunately you can't get nonzero to sort the indices the other way, so I resorted to this goofy hack.

Bt = B.transpose(copy=True) # Have to copy or it transposes in place

Bt.sort_indices()

# note that col,row are transposed

(col_indices, row_indices) = Bt.nonzero()

print(B[row_indices col_indices].A1)

    # [1 3 4 6 7 8 2 5 9]

    # Same as B.data!!!

So in my code now, I compute my "full structure matrix" that has 1's anywhere data *could ever be nonzero, and then get the indices from the transposed version of that:

self.A_full = compute_full_structure_matrix()

A_full_transposed = self.A_full.transpose(copy=True)

A_full_transposed.sort_indices()

(self.A_full_col_indices, self.A_full_row_indices) = A_full_transposed.nonzero()

# Then later, I can update from A_current

self.A_current = get_state_matrix_A()

Ax_values = self.A_current[self.A_full_row_indices, self.A_full_col_indices].A1

self.osqp_solver.update(Ax=Ax_values)

I just tried this and it works. It's hard to work this into the existing python API without storing the row and column indices from setup, but that could be done optionally. I could submit a PR with something like that, if there's any interest.

[Goran Banjac]

Hi Chris,

Thanks for using OSQP.

You need the nonzero elements sorted in a column-major way because the CSC format stores the data that way. This is why A.data in this example returns the nonzero elements sorted correctly.

Thanks for being willing to submit a PR, but I am not sure we want to include the procedure above in the OSQP code. There are several reasons for this.

1. Performing all these matrix operations would reduce the efficiency of the solver. For instance, if the matrix has a huge number of nonzero elements and we want to update only a few of them, the procedure you suggested above would slow down the update method considerably.
2. We would like to keep all the OSQP interfaces as similar as possible. The fact that scipy nonzero method returns the indices in a different order does not justify changing the entire interface. On the other hand, Matlab find function returns row and column indices of nonzero elements in a column-major way.
3. The best way of extracting the nonzero elements (in the correct order) probably depends on the application. For instance, you could have used the fully occupied structure matrix MASK (with ones at possible nonzero locations) and then add eps*MASK to A_new, where eps is say 1e-12. Then you can extract the nonzero values (in the correct order) by calling A_new.data and then subtract eps from your data vector.

I know that the interface for matrix updates is a bit confusing, but it is definitely more efficient than requiring the user to send the whole matrix that needs to be updated.

Best,

Goran

[Christopher Horvath]

Thanks for your detailed response! The mask & eps trick is cool, I'll give that a try and compare performance, though the approach I already have is doing great so I'm not too worried.

This is a great library, thanks for all the work you've put into it.

Chris

To unsubscribe from this group and all its topics, send an email to osqp+uns...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/osqp/fa833427-01d5-4411-8fcf-32e8a8dad17e%40googlegroups.com.