Adapt the novel step eq. form also to other DCS types.

Either remove or make robust options currently listed as experimental.

We have this for the plain numerical time, to make an ell_1 relaxation of the terminal time constraint.
It would be nice to have the same also for physical time constraints (maybe even also for the stagewise constraints - but to be discussed).

I want also to test with this the ease of use of vdx, once we migrated to it.

Hi,
Thanks for the nice tool and quite interesting research papers.

I am not sure how to formulate one hybrid dynamics problem in order to solve it with Nosnoc. In simple terms, it is as follows:

master; % integer, influences the speed capabilities by switching the profiles 
u; % continuous control

if master==0
    f = profile_0; % define dynamics profile 0
elseif master==1
    f = profile_1; % define dynamics profile 1 
end

% goal: achieve *maximum speed*. 
% Determine best `u` (continuous control) and `master` (profile sequence e.g. profile_0, profile_1, profile_0)

I have checked examples, but I start to think you might not be dealing with these cases. I saw that you are defining the gap function. However, in this case the gap function should be kinda determined by the algorithm and it might be changing in time as well.

Would really appreciate feedback.

Best,

If a problem gets infeasible, we can increase the sigma value once or twice to see will the homotopy loop recover.
This should be an option called, e.g., nonmonotone_homotopy_parameter, and we need a "kappa_2" for how much we update it.

If problem fails: sigma_k = kappa_2*sigma_k, % kappa2 >1

Currently all constraints in nosnoc.model.base.g_path are either evaluated on control stage boundaries, finite element boundaries, and/or rk stages uniformly. This can in some cases cause duplicate constraints. As a result we should a priori split the path constraints based which of the three levels they make sense to apply.

Currently it is looped in ControlStage and it might not even work. Moreover the current for loop will be broken in newer matlab versions. R2024a already throws a warning.

In CLS, we have the option to lift the velocity state for numerical inversion of M. Add analogous functionality to the cls time freezing reformulations.

I believe this is only easy to do for the multicontact reformulation, but this needs more thought.

Most examples use the old nosnoc. As new design is coming in, we should translate all examples and update the comments such that they match the new features.
For example, tutorial examples should be more timely (and maybe have only two files - one ocp and one simulation, or one per major feature: stewart/heaviseide, pds, time freezing, cls)

Based on discussion in #62

We want to make the front end of mpccsol more clear and as such need to consolidate the relaxation types we provide.

Based on #62

Re-implement functionality to calculate the type of stationarity point in RelaxationSolver.

Add the three examples from the video.

The path complementarity constraints should have optional upper bounds, as this can improve convergence.

Based on #62

NosnocIntegrator needs its own separate options.

It would be useful to decouple the nosnoc specific sections of the homotopy solver from the core loop and provide it as an interface to solve generic MPCCS.

As we will have a new integration and OCP class after the refactoring, it would be nice to add afterwards an mpc loop that solves a sequence of OCPs.

This will also have some options, e.g. how to initialize? To use an early homotopy iterate or the solution of the previous step?
For having the new x0, use the same integrator (e..g., just shift ocp solution) or simulate with a different nosnoc integrator (e.g. higher order, more intermediate FEs).

Currently time-freezing only supports the Heaviside-step reformulation. Add functionality to be able to use the Stewart reformulation as well.

First implementation of new design for nosnoc. For details see https://github.com/nurkanovic/nosnoc/tree/class-org-plan/organization_proposal .

Verify correctness of reformulations via more granular unit tests.

This option was experimental in the mpcc homotopy and proved not to bring much.
I suggest removing it completely as it adds no value but requires maintenance.

Currently Stewart and Heaviside problems use Options.time_freezing to identify a physical time state. This makes the resulting implementation confusing. A possible solution to this is just adding a new option that specifies the existence of a physical time state.

Based on conversation in #62.

Re-implement the fixed active set NLP extraction for finding a polished solution.

To reduce the size of the repository we should use syscop.de links for the papers.

Probably with the syntax:

ocp_solver.set(var, field, varargin)

where varargin contains indices and the value to set.

New README should contain a clearer explanation of the steps that NOSNOC uses to solve an OCP/simulate a system.

Currently the results struct has lots of noise in it with possibly empty entries. Remove these and document what we store in the results.

This may be irrelevant if we use vdx.

Currently the reference values for quadratic costs are fixed at model time. Make them time varying parameters to allow for more flexibility.

This includes duplicate names for irk schemes and other things that have been moved to Enums.

By introducing a slack variable, MPVCs are equivalent to MPCCs. This can always be done manually, but it would be maybe useful to automate this step if the user indicates that the given problem is a mpvc (or we automatically recognize it by looking at the bounds of G and H) and transform into an mpcc.
maybe it makes sense to rename everything that is mpcc into mpec.

Currently, we support ell_1 relaxation of some terminal constraints

1. the terminal constraints of the ocp
2. some of the FESD/time-freezing terminal constraints to meet a certain final numerical/physical time.

I would still steer the relaxation of these constraints with appropriate separate options, as we have them now.

In addition to this, it would be nice to be able to relax other constraints (e.g. using vdx internally to define according variables and constraints). In particular it would be nice to relax:

1. Path constraints
2. Dynamics constraints (maybe to separate even between the differential and algebraic part of the dynamics)

Choices for relaxing:

1. ell 2 (every constraint gets a separate slack), maybe we can have upper bounds on slacks to avoid unbounded problems
2. ell1 (every constraint gets a separate slack)
3. ell_inf (every constraint group gets one slack)

Choice for penalty parameter:

1. fixed value (separate parameter for every constraint group, by groups i mean: path, dynamics, terminal, terminal_algoritmic (nice and instructive names are a plus))
2. use rho = 1/sigma - we should decide if we can mix fixed values and homotopy penalities. one one hand it gives a lot of flexibility, on the other hand too many choices might be confusing.

Based on #62

Currently indexes are tracked independently at different levels and must be handled somewhat manually.
Several things need to be done to make maintenance of NOSNOC easier:

• Remove FiniteElement and ControlStage classes
• Use a smarter index tracking utility to avoid pre-definition, etc. Something similar to CasADi "structures" in casadi.tools

Currently, we have for every FE an h.
Alternatively, we can define a clock state t, and have a "control variable" h_n to determine the length of a finite element.
We simply integrate the clock state as t_{n+1} = t_n + h_n; and enforce a "terminal" constraint
t_N = T_final instead of sum_n h_n = T_final
The intermediate variables t_n help eliminate sum_n h_n = T_final, which kills the block banded structure.

Combine use_speed_of_time_variables and local_speed_of_time_variable into SpeedOfTimeMode with entries NONE, LOCAL, GLOBAL

The current interface for mpccsol is designed to be similar to nlpsol i.e. mpccsol(solver_name, plugin, problem, options). In its current state we only support nosnoc.solver.MpccSolver as a backend for mpccsol which is a relaxation solver for MPCCs. This means plugin is used to pass the relaxation method to mpccsol. This is in my opinion problematic as it both makes changing the relaxation mode a bit unintuitive and will become more confusing as other solver backends are added.

We have as I see it several options:

1. Continue as is but then add other plugins.
2. Drop the plugin as a part of the interface and simply pass it as part of the options.
3. Change the plugin name for the current MpccSolver backend to something like relaxation (name open to change) and have the relaxation mode as part of the solver options.

I view 3 as the best option as it is relatively ergonomic and makes it clear which solver backend is being used. 2 is also passable for me though it moves away from the nlpsol interface which may be undesirable in terms of user ergonomics.

clear all
clc
close all
import casadi.*
%% model parameters
e = 0.2; u_max = 9; beta = 0.0; mu=0;
%% NOSNOC settings
problem_options = NosnocProblemOptions();
solver_options = NosnocSolverOptions(); %% Optionally call this function to have an overview of all options.

problem_options.n_s = 3;
solver_options.mpcc_mode = 'Scholtes_ineq';
solver_options.homotopy_update_rule = 'superlinear';
% problem_options.step_equilibration = 'direct_homotopy';
problem_options.use_fesd = 1;
problem_options.time_freezing = 1;

problem_options.local_speed_of_time_variable = 1;
problem_options.stagewise_clock_constraint = 0;
problem_options.nonsmooth_switching_fun = 0;
problem_options.pss_lift_step_functions = 0;

problem_options.no_initial_impacts = 0;
% solver
solver_options.opts_casadi_nlp.ipopt.max_iter = 1e5;
solver_options.opts_casadi_nlp.ipopt.tol = 1e-7;
solver_options.opts_casadi_nlp.ipopt.acceptable_tol = 1e-5;
solver_options.use_previous_solution_as_initial_guess = 0;
solver_options.N_homotopy = 6;
solver_options.print_level = 3;
solver_options.comp_tol = 1e-6;
solver_options.elastic_scholtes = 1;
%% model parameters
g=9.81;
q_target = [4;0.5];

%% solver
model = NosnocModel();
model.mu_f = mu;

problem_options.T = 4;
problem_options.N_stages = 20;
problem_options.N_finite_elements = 3;
% model equations
q = SX.sym('q',2);
v = SX.sym('v',2);
u = SX.sym('u',2);

model.lbx = [0;-0.001;0;-inf]; model.ubx = [4;inf;inf;inf];
model.x0 = [0;0.5;0;0];
model.x = [q;v]; model.u = u; model.e = e ;
model.f_c = q(2)-0;
model.f_v = [u-[0;g]-betavsqrt(v(1)^2^2+v(2)^2+1e-3)];
model.J_tangent = [1; 0];
model.dims.n_dim_contact = 2;

% Objective and constraints
model.f_q = u'u; model.f_q_T = 100v'*v;
model.g_path = u'*u-u_max^2;
model.g_terminal = q-[q_target];

mpcc = NosnocMPCC(problem_options, model);
solver = NosnocSolver(mpcc, solver_options);
[results,stats] = solver.solve();

Currently with the new solver syntax TestMpccCrossComp is a test that really tests multiple modules. We should re-architect this.

Rework how we slice the step equilibration options. Particularly with the introduction of relaxation mode support in vdx we can perhaps build the step equilibration out of several building blocks.