Penalty method in Rust

[alphaville]

Is your feature request related to a problem? Please describe. Implement the penalty method in Rust.

Describe the solution you'd like

-We first need to define a type (a closure) for penalty functions. These should not depend on the weights.

• We need something like icasadi to be able to interface such penalty functions written in C (generated by CasADi).
• We need a new trait that will define the penalty method and will allow the use of any numerical algorithm for solving the individual problems.

Describe alternatives you've considered We can either extend icasadi into a generic code generation tool, or create a new module for penalty functions. I'd prefer the first option.

Additional context We should keep in mind that the next step should be the implementation of the Augmented Lagrangian method. In fact, the user may choose to use the penalty method on some constraints and AL on other ones.

[korken89]

This is high on my list as well, collision avoidance in drones work a lot faster (and more stable) with the penalty method

[alphaville]

Notes:

The penalty method and the augmented Lagrangian method can be used to solve (parametric, in our case) problems of the form

Minimize_{u in U(p)} phi(u; p)
subject to: c_i(u; p) = 0, 
            for i=1,...,N

The penalty method requires the solution of problems of the form

P_k: Minimize_{u in U(p)} phi(u, p) +  SUM_i μ_{i,k} g(c_i(u)) =: F(u, q),

where q is the augmented parameter vector q = (p, mu) and g can be the identity function, g(z)=z, or g(z) = z^2 and μ_{i,k}>0 are positive scalars.

Problem P_0 is first solved with some arbitrary μ_{i,0}>0 from which we obtain a solution u^*_0. Then, if g(c_i(u^*_0)) > epsilon, we increase the value of mu_{i, 1} by multiplying it by a>1 (e.g., a=10). The solution at iteration k is used as the starting point at iteration k+1.

In order to encode inequality constraints of the form c_i(u) <=0, we can use c_i'(u) = max{0, c_i(u)} and require that g(c_i'(u)) = 0. The typical choice of g will be g(z) = z^, but it's good to allow the user to provide some other function if necessary.

The user will have to specify:

• A function c(u) = (c_1(u), ...., c_N(u)), i.e., a mapping from R^n to R^N
• Function g [default: g(z) = z^2]
• Function F(u, q) and its gradient (which will be computed by the code generator using CasADi)
• The initial weights μ_{i,0} [default: vector of ones]
• The multiplicative factor a>1 [default: a=10]
• The maximum number of outer iterations (of the penalty method) [default: 6]

Associated branches:

• feature/homotopy

See also:

• This whitepaper on overleaf

[alphaville]

Update: currently under development (see #21)