Trouble solving "simple" problem

aknight0 · April 1, 2020, 12:16am

I’m having trouble solving what seems to be a fairly straightforward problem. I’m using a 2D single integrator model, trying to move from (0,0) towards (10,5), with a two norm input constraint, something like:
opt_problem

What I would expect is a straight line from the origin pointing towards (10,5). With a one second horizon, the answer should be (0.89443,0.44721). For some reason though, Acados is unable to converge the nonlinear inequality constraint. The solver doesn’t seem to make progress towards the solution, and doesn’t seem to respect the constraint bounds, though without a converged solution that may not mean much.

What is strange is that IPOPT is able to solve the same problem in just a handful of iterations, producing the expected result. Even when initialized with the solution from IPOPT, Acados does not converge, which makes me wonder if I am missing something. I tried with HPIPM, OSQP, and QPOASES, with similar results from all three.

Is there a reason this specific type of problem might be difficult for Acados to solve?

OCP script:

gist.github.com

https://gist.github.com/abe-mart/ca631fdbd225a93791b8e3e151edae3e

run_ocp.m

%% Control 2D single integrator

clear all

%% discretization
N = 10;
T = 1; % time horizon length
x0 = [0; 0];

nlp_solver = 'sqp'; % sqp, sqp_rti

This file has been truncated. show original

Model:

gist.github.com

https://gist.github.com/abe-mart/acc983005a3baf8cfd3576836b4227c5

model_SI_2D.m

function model = model_SI_2D()
%2D_SI_model Super simple 2D single integrator model

import casadi.*

%% Differential states
x = SX.sym('x',2);

sym_x = x;
sym_xdot = SX.sym('xdot',numel(sym_x),1);

This file has been truncated. show original

FreyJo · April 2, 2020, 10:29am

Hi @aknight0,

I just had a look at your example.
The problem is that the Hessian with respect to u is basically indefinite.

You can recover convergence by:
a) adding a penalty on u.
In your example:

expr_ext_cost = (x(1)-10)^2 + (x(2)-5)^2 + 1 * dot(u,u);

b) I recently implemented the option to add Levenberg-Marquadt regularization term to the Hessian, basically adding a multiple of the unit matrix to the Hessian.
It is not merged yet, but you can find it here:

With this feature, you can also do:

ocp_opts.set('levenberg_marquadt', 1.0);

Cheers!

aknight0 · April 1, 2020, 9:21pm

Brilliant! I pulled your Levenberg-Marquadt branch down, and using that option solved my problem. If I recall correctly, I think IPOPT does something similar to try to regularize the Hessian? Now I’m getting equivalent answers from IPOPT and Acados, which is great. Is the 1.0 you set a binary switch, or can this be applied in varying magnitudes?

As a side note, adding a penalty on u also helped, but took more iterations and also required that I reduce the solver step size. The regularization method solved in fewer iterations with the default step size.

Thanks for putting up with my endless stream of questions! This is a great project, and you’ve been extremely helpful.

FreyJo · April 2, 2020, 8:49am

Great!
The Levenberg-Marquadt value is the factor for the unit matrix that is added to the Hessian.
Both options to fix the convergence issue have some tuning parameter and I would not say one is strictly better than the other.

Happy to hear that!