Use casadi.interpolant in Model Formulation

Jannis · December 3, 2024, 3:40pm

Hi

as part of my model formulation I want to use a look-up table.

I first define my look up talbe:

PH2Fcn = interpolant(‘LUT’,‘bspline’,{FCS.Pnet},FCS.PH2);

then I create a symbolic variable:

PFCSBus = SX.sym(‘PFCSBus’)

when I call:

PH2 = PH2Fcn(PFCSBus );

I get the following error:

Error in Function::call for ‘LUT’ [BSplineInterpolant] at …/casadi/core/function.cpp:370:
…/casadi/core/function_internal.cpp:2024: ‘eval_sx’ not defined for BSplineInterpolant

After searching for a bit I realized that casadi.interpolation only support MX variables. So I tried to rewrite my problem formulation using MX variables, but this seems to be a problem for ACADOS:

constraint detection only works for casadi.SX!

The two solution that I come up with are, first, to approximate the look-up talbe using a polynomial, but this leads to a non ideal solution and second, to manuale code a linear interpolation routine. This works, but then I always get a solver status == 2, when simulation the MPC.

Does anyone know of a better solution for the descriped problem?

Thanks Jannis

kaethe · December 3, 2024, 6:25pm

Hi Jannis,

no need to approximate the look-up table. Just use MX symbolics and set the constraint type to 'BGH' explicitly (The defualt type is 'AUTO' which results in acados trying to detect the constraint type. This only works for SX atm)

Best, Katrin

Jannis · December 6, 2024, 1:02pm

Hi Katrin,

thank you for the response. That solved the Problem. But now the solver status is most of the times equal to two (= ACADOS_MAXITER).

I’ve tried to solve this, but I have no idea what the possible reasons for this could be. I there a proper way of debugging this solver status?

Thanks Jannis

kaethe · December 6, 2024, 1:07pm

Maybe simply increasing nlp_solver_max_iter might fix it

As a first step you should check the residuals with solver.print_statistics() to check whether it’s the dynamics or inequality constraints that cause the problem.

https://docs.acados.org/troubleshooting/index.html

Jannis · December 6, 2024, 1:48pm

Sorry, I should have mentioned in my previous post that increasing ocp.solver_options.nlp_solver_max_iter does not help, only the calculation time increases but solver status = 2 remains.

Currently I just use the generated s-function for simulation, but I will try to build a script-based simulation and check the solver statistics.

Thanks for the hint

Jannis

david0oo · December 10, 2024, 10:49am

Hi Jannis,

it could also be that full step SQP does not converge for your problem. In that case, you could try to use the new Funnel globalization which is available in acados. For an example on how to use it, you can consult:

github.com

acados/acados/blob/master/examples/acados_matlab_octave/convex_problem_globalization_needed/convex_problem_globalization_necessary.m

%
% Copyright (c) The acados authors.
%
% This file is part of acados.
%
% The 2-Clause BSD License
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE

This file has been truncated. show original

Additionally, I would set the option
ocp.solver_options.globalization_full_step_dual = true;
Looking at the statistics table always helps, you can also show the statistics, if you set
ocp.solver_options.print_level = 1;

I hope this helps.

David

Jannis · December 16, 2024, 11:50am

Hi Katrin, Hi David,

this is an example of printed solver stats for solver status = 2. I am not quite sure how to interpret these information and how it helps me to debug my solver.

iter	res_stat	res_eq		res_ineq	res_comp	qp_stat	qp_iter	alpha
0	2.404421e-02	9.661000e-03	0.000000e+00	1.940592e-01	0	0	0.000000e+00
1	3.098387e-02	1.167844e-02	0.000000e+00	7.568737e-02	0	5	8.235430e-02
2	2.441422e-02	1.046383e-02	0.000000e+00	2.627153e-02	0	5	1.680700e-01
3	4.967150e-02	9.860611e-03	0.000000e+00	1.493769e-01	0	5	5.764801e-02
4	2.441411e-02	9.673542e-03	0.000000e+00	1.494973e-01	0	5	5.764801e-02
5	4.439071e-02	1.167613e-02	0.000000e+00	7.333520e-02	0	5	1.176490e-01
6	2.294103e-02	1.030245e-02	0.000000e+00	8.656656e-02	0	5	1.176490e-01
7	2.404419e-02	9.708531e-03	0.000000e+00	1.133045e-01	0	5	5.764801e-02
8	4.064227e-02	1.263835e-02	0.000000e+00	6.634963e-02	0	5	1.680700e-01
9	2.876244e-02	1.190977e-02	0.000000e+00	1.481476e-01	0	5	5.764801e-02
10	2.876244e-02	1.190977e-02	0.000000e+00	1.481476e-01	0	6	5.764801e-02

I also tried to set ocp.solver_options.globalization_full_step_dual = true;, but this does not seem to change anything.

Thanks for the help

Jannis

david0oo · December 16, 2024, 1:23pm

Hi Jannis,

thanks for sending the a log. Unfortunately, it is not very insightful for me at the moment. Could you share your code, then I can try to figure out a solution.
I have two ideas that you could try out and a question:

What Hessian matrix are you using? If you use the exact Hessian, you should use some kind of regularization. Here I recommend to use solver_options.regularize_method = 'MIRROR'

Use the funnel globalization: This means use the following options
solver_options.nlp_solver_type = 'SQP'
solver_options.globalization = 'FUNNEL_L1PEN_LINESEARCH' solver_options.globalization_full_step_dual = True solver_options.print_level = 1

Can you post the output here for the funnel method?

The res_ineq is always zero. Does this mean that you do not have constraints? If the only constraints are your dynamics, you could use DDP as solver. This can be done using the following options:
ocp.solver_options.nlp_solver_type = 'DDP'
ocp.solver_options.globalization = 'MERIT_BACKTRACKING'

I hope this helps. If it does not work, I would appreciate sending a minimal example such that I could try by myself to resolve your issue.

Best,
David

kaethe · December 16, 2024, 3:41pm

Adding to David’s suggestions, I recommend to solve the QP subproblems more accurately than the NLP, i.e. set the values of qp_solver_tol_stat to nlp_solver_tol_stat/10. and similarly for the equality, inequality and complementarity tolerances.

The inequality and complementarity residuals look indeed a bit weird. While the inequality constraints seem to be satisfied from the beginning, the complementarity residual is not really converging. What kind of inequalities did you impose? Are they active at the solution that you obtain?

Jannis · December 17, 2024, 2:27pm

Hi David,

to your two questions

I tried to use the funnle globalization, but It does not seem to work. I get the following error:

FUNNEL_L1PEN_LINESEARCH only supports SQP.

although I am using:

ocp.solver_options.nlp_solver_type = ‘SQP’
I have upper and lower bound on inputs and states and a soften equality constraint.

Here you can find a (not so minimal) minimal Example: minimalExample.zip - Google Drive

When Acados is initialized you should be able to run minimalExample.m

THanks for the Help
Jannis

kaethe · December 18, 2024, 1:36pm

I think you need to update acados. On current master, your example works also with funnel globalization (still running into maximum iterations though…)

Is there anything special about your model? The equality constraints are not converging. Did you try to simulate the system using an acados integrator separately?