Creating compatible files from MXFunctions that include solvers

[CheeGY]

Hi again!

Thanks for your quick reply for the other issue. I believe using expand() would most likely solve the problem related to some CasADi MXFunctions, but unfortunately, it didn't work out for my case.

I think it is a related, but separate, issue, and hence I'm opening another issue instead of re-opening that.

The issue is: The expand() option is not allowed for some type of MXFunctions, in particular those that include a solver. Here's an (minimal) example, modified from one of the CasADi examples (https://web.casadi.org/blog/opti/):

from casadi import *

opti = casadi.Opti()

x = opti.variable()
y = opti.variable()
r = opti.parameter()

cost = (1 - x)**2 + (y - x**2)**2
con = x**2 + y**2 <= r
opti.minimize(cost)
opti.subject_to(con)
opti.solver('ipopt')

func = opti.to_function('func', [r], [cost])

If I call func.expand(), I encounter the following error:
RuntimeError: Error in Function::operator() for 'solver' [IpoptInterface] at .../casadi/core/function.cpp:1377: .../casadi/core/function_internal.cpp:1832: 'eval_sx' not defined for IpoptInterface

I have also tried to change the solver to 'sqpmethod' with qpsol='qrqp' instead of 'ipopt', but that didn't solve the problem.

Are there other options to get the corresponding .casadi file for func, without using expand()?

---

In the paper, it appears that a few ideas related to this above issue have been brought up:

1. In Section V-A, "MPC Parallelization for the MIT Humanoid", it is stated that "we approximate the solution to an OCP with a strictly fixed number of operations".

2. In the same section, "using the previous solution as an initial guess, a single QP iteration is often sufficient to approximate the solution."

3. "For the MIT Humanoid, we implement the single rigid body model (SRBM) nonlinear MPC controller detailed in [40] entirely in casadi, and demonstrate its subsequent CusADi parallelization across 4,096 environments in IsaacGym, as shown in Fig. 1."

Is the implementation of this "single rigid body model (SRBM) nonlinear MPC controller" available in the repository? I think that might be useful in solving the above issue.

Thanks again!

[se-hwan]

TLDR:

• expand() is necessary
• We use Opti to conveniently form our matrices, we don't create a function around the Opti() instance itself
• I'm working on some clean tutorials and examples for individual parts of the MPC problem, and will add them either to this repo or a separate optimization focused one.

Happy to help!

So in your example, if you want to parallelize the cost function, you might want to try this:

fn_cost = casadi.Function('cost', [x, y], [cost]) fn_cost.save('fn_cost.casadi')

And then the saved output could be codegened with CusADi, and the cost function could be parallelized across any number of environments. Note that your expression for the cost depends only on x and y, not r.

Sorry if it wasn't clear in the paper, we only use the Opti interface to get the expressions for our QP matrices, within an SQP iteration (because it's convenient, we could do it manually to if we wanted).
While CasADi lets you create a function around the whole Opti stack, this includes the solver, which is just an interface to a third party solver (IPOPT, qrqp, etc. which cannot be parallelized with CusADi, because these aren't defined with CasADi symbolics).

As an example, we can use Opti to define all sorts of nonlinear constraints and costs that we want, then start defining the matrices in the QP subproblem. Namely, H = hessian(opti.cost, opti.x), A_ineq = jacobian(opti.g, opti.x) etc.

We end up with six symbolic matrices/vectors - H, c, A_eq, b_eq, A_ineq, b_ineq that define our QP problem.
Then, we actually implement our own QP solver in CasADi. This involves using the penalty method to move the inequality constraints into the cost, then solving the corresponding optimization problem (quadratized cost, linear constraints).

This KKT system is a linear system (Ax=b matrix equation), so we can address it with CasADi's ldl or its solve method (which I believe does some QR decomposition to solve the system). Check out the docs for these functions here: https://web.casadi.org/python-api/#document-api
Of course, if you want to manually implement some more efficient linear solution, please feel free to do so!
We can repeat this process (a fixed number of times) until the QP converges. Note that this is only a single QP for a larger SQP problem, but we mention in the paper that we take "real time iterations", where we can often approximate the solution well with a single QP subproblem.

We don't have a clean version of the SRBM implementation yet - I'm working on making some CasADi tutorials that detail the individual sections (the QP matrices, the linear solver, SRBM dynamics, etc.).
I'm in between whether I should add the examples directly to this repo, or keep the CUDA codegen and optimization code separate.

Hope that clarified some of your questions - I'll leave this issue up so other people can see it.
If you have any other questions or points of confusion, feel free to add more replies or message individually, happy to describe more details if needed.

[CheeGY]

Thanks again for the quick reply!

What I wanted to do is to see if I can solve multiple instances of an optimal control problem on GPU (I guess I was drawn to the Optimal Control part, rather than the Symbolic Expressions part in the title of the paper, "A GPU Parallelization Framework for Symbolic Expressions and Optimal Control".)

Just to be clear(er):

1. The original optimal control problem ((1) or (2), or equivalently (7)) is cast into an SQP framework, which consists of QP subproblems.
2. The paper adopts a RTI approach, which solves only one of these QP subproblems.
3. This QP subproblem is solved by considering the KKT conditions, which works out to be a set of linear equations, which is solved using ldl or solve. And either of those can be saved as a .casadi file if I use expand()?

A follow-up question:
Since CasADi is able to transfer the symbolic OCP, through one of the interfaces, e.g., the ipopt interface shown in the error above, and pass it to the solver, e.g., ipopt or qrqp, I guess the main issue here is (still) with the MXFunction. If the OCP can be written without MXFunction, then expand() is not required. Consequently, the corresponding .casadi file can be obtained and codegen can be done. In other words, we can still use the solvers in CasADi, as long as there are no MXFunction. Is that right, or am I missing something?

---

Looking forward to those tutorials/code! I'm curious to see how the OCP can be constructed (without the MXFunction at the end) such that it can be ported over to CusADi for parallel evaluation on GPU.

In particular, I think it will be super useful if there is an example that illustrates how the OCP in (1) and (2) in the paper can be cast into CasADi symbolics, put together with the (QP) solver, and do multiple parallel evaluations of that OCP on GPU.

[se-hwan]

1. Yes
2. Yup
3. Exactly, the solution from CasADi's LDL can be saved as a CasADi function, which can then be expand()ed, then saved for parallelization.

For your followup, I'm unsure if you mean codegen for GPU parallelization or CasADi's native codegen for C compilation. Yes, CusADi can only parallelize SXFunctions. This isn't a fundamental limitation - we haven't added MX capabilities because almost all MX expressions can be expanded, and SX expressions are typically faster to evaluate.

• For native C compilation, yes, you can codegen the your whole Opti instance and compile it. Under the hood, if you set opti.solver('ipopt'), the .c codegen externally calls ipopt (or whatever solver you chose) libraries.
• For parallelized optimization, there is no point setting opti.solver to anything - we don't use the third party solvers, because we cannot parallelize their computations as part of the function expression graph (essentially a black-box third party node). This is why we make our custom symbolic solver.

Summary:
In general, I think the Opti stack formulates all its expressions as MX, and if you want to solve the single OCP instance with CasADi, there is no problem with using MX expressions with the solvers like ipopt.
CusADi (currently) only parallelizes SX expressions, but you cannot include a third party solver as part of that SX expression.

Hope that clarifies your questions, let me know if anything's still unclear.

[CheeGY]

Thanks!

[se-hwan]

Moving this to discussions since we're talking about implementation details rather than CusADi specific bugs.