| id | openrust-basic |
|---|---|
| title | Introduction |
| description | How to use OpEn directly in Rust |
OpEn can solve problems of the form:
where
The definition of an optimization problem consists in specifying the following three componenets:
- the cost function
$f$ as a Rust function - the gradient of
$f$ ,$\nabla f$ , as a Rust function - the set of constraints,
$U$ , as an implementation of a trait
The cost function f is a Rust function of type |u: &[f64], cost: &mut f64| -> Result<(), SolverError>. The first argument, u, is the argument of the function. The second argument, is a mutable reference to the result (cost). The function returns a status code of the type Result<(), SolverError> and the status code Ok(()) means that the computation was successful. Other status codes can be used to encode errors/exceptions as defined in the SolverError enum.
As an example, consider the cost function
let f = |u: &[f64], c: &mut f64| -> Result<(), SolverError> {
*c = 5.0 * u[0] - u[1].powi(2);
Ok(())
};The gradient of the cost is a function df with signature |u: &[f64], grad: &mut [f64]| -> Result<(), SolverError>. The first argument, u, is again the argument of the function. The second argument, is a mutable reference to the result (gradient). The function returns again a status code (same as above).
For the cost function
This function can be implemented as follows:
let df = |u: &[f64], grad: &mut [f64]| -> Result<(), SolverError> {
grad[0] = 5.0;
grad[1] = -2.0*u[1];
Ok(())
};Constraints implement the namesake trait, Constraint. Implementations of Constraint implement the method project which computes projections on the set of constraints. This way, users can implement their own constraints. OpEn comes with the following implementations of Constraint:
| Constraint | Explanation |
|---|---|
AffineSpace |
|
Ball1 |
|
Ball2 |
|
BallP |
|
Sphere2 |
|
BallInf |
|
Halfspace |
|
Hyperplane |
|
Rectangle |
|
Simplex |
|
NoConstraints |
|
FiniteSet |
|
SecondOrderCone |
|
Zero |
|
CartesianProduct |
Cartesian products of any of the above |
These are the most common constraints in practice.
The construction of a constraint is very easy. Here is an example of a Euclidean ball centered at the origin with given radius:
let radius = 0.5;
let bounds = constraints::Ball2::new(None, radius);Having defined a cost function, its gradient and the constraints, we may define an optimization problem. Here is an example:
let problem = Problem::new(&bounds, df, f);Note that problem now owns the gradient and the cost function
and borrows the constraints.
OpEn uses two essential structures: (i) a cache and (ii) an optimizer.
Rarely will one need to solve a single optimization problem in an embedded engineering application.
The solution of an optimization problem requires the allocation of memory.
A cache (see PANOCCache) allows for multiple instances of a problem to have a common workspace, so that we will
not need to free and reallocate memory unnecessarily.
A cache object will be reused once we need to solve a similar problem; this is the case with model
predictive control, where an optimal control problem needs to be solved at every time instant.
let n = 50;
let lbfgs_memory = 10;
let tolerance = 1e-6;
let mut panoc_cache = PANOCCache::new(n, tolerance, lbfgs_memory);The last necessary step is the construction of an Optimizer. An Optimizer uses an instance of the problem adn the cache to run the algorithm and solve the optimization problem. An optimizer may have additional parameters such as the maximum number of iterations, which can be configured. Here is an example:
let mut panoc = PANOCOptimizer::new(problem, &mut panoc_cache)
.with_max_iter(max_iters);We may then call the solver using the method solve providing an initial guess:
let status = panoc.solve(&mut u).unwrap();This will return the solver status (if it is successful) and will update u with the solution.
Note: The algorithm may, in general, return an error (instance of SolverError), this is
why the caller function should not use unwrap directly, but should rather check whether the
solver returned an error. The error can be propagated upstream using let status = panoc.solve(&mut u)?;.
Method solve, if successful, returns an object of type SolverStatus which stores information about the solver, namely, whether the required tolerance was attained (has_converged()), the exit status (exit_status()), the number of iterations (iterations()), the norm of the residual (a measure of the accuracy of the solution) (norm_fpr()) and the cost value at the approximate solution (cost_value()).
The minimization of the Rosenbrock function
is a challenging problem in optimization.
The Rosenbrock function in two dimensions with parameters
with gradient
that is,
pub fn rosenbrock_cost(a: f64, b: f64, u: &[f64]) -> f64 {
(a - u[0]).powi(2) + b * (u[1] - u[0].powi(2)).powi(2)
}
pub fn rosenbrock_grad(a: f64, b: f64, u: &[f64], grad: &mut [f64]) {
grad[0] = 2.0 * (u[0] - a) - 4.0 * b * u[0] * (-u[0].powi(2) + u[1]);
grad[1] = 2.0 * b * (u[1] - u[0].powi(2));
}Here, we minimize the above two-dimensional Rosenbrock function with parameters
The required tolerance is
fn main() {
/* USER PARAMETERS */
let tolerance = 1e-14;
let a = 1.0;
let b = 200.0;
let n = 2;
let lbfgs_memory = 10;
let max_iters = 80;
let mut u = [-1.5, 0.9];
let radius = 1.0;
// define the cost function and its gradient
let df = |u: &[f64], grad: &mut [f64]| -> Result<(), SolverError> {
rosenbrock_grad(a, b, u, grad);
Ok(())
};
let f = |u: &[f64], c: &mut f64| -> Result<(), SolverError> {
*c = rosenbrock_cost(a, b, u);
Ok(())
};
// define the constraints
let bounds = constraints::Ball2::new(None, radius);
/* PROBLEM STATEMENT */
let problem = Problem::new(&bounds, df, f);
let mut panoc_cache = PANOCCache::new(n, tolerance, lbfgs_memory);
let mut panoc = PANOCOptimizer::new(problem, &mut panoc_cache)
.with_max_iter(max_iters);
// Invoke the solver
let status = panoc.solve(&mut u);
}This example can be found in examples/panoc_ex1.rs.
In embedded applications, we typically need to solve parametric problems, that is, problems of the form
where
As mentioned above, a cache needs to be constructed once.
Then, every time the value of
Let us give a complete example. We consider the problem of minimizing the two-dimensional
Rosenbrock function,
In the following example, we solve 100 problems while modifying the parameters,
Note also that the initial guess of each problem is the solution of the previous problem.
use optimization_engine::{constraints::*, panoc::*, *};
fn rosenbrock_cost(a: f64, b: f64, u: &[f64]) -> f64 {
(a - u[0]).powi(2) + b * (u[1] - u[0].powi(2)).powi(2)
}
fn rosenbrock_grad(a: f64, b: f64, u: &[f64], grad: &mut [f64]) {
grad[0] = 2.0 * u[0] - 2.0 * a - 4.0 * b * u[0] * (-u[0].powi(2) + u[1]);
grad[1] = b * (-2.0 * u[0].powi(2) + 2.0 * u[1]);
}
fn main() {
let tolerance = 1e-6;
let mut a = 1.0;
let mut b = 100.0;
let n = 2;
let lbfgs_memory = 10;
let max_iters = 100;
let mut u = [-1.5, 0.9];
let mut radius = 1.0;
// the cache is created only ONCE
let mut panoc_cache = PANOCCache::new(n, tolerance, lbfgs_memory);
let mut i = 0;
while i < 100 {
// update the values of `a`, `b` and `radius`
b *= 1.01;
a -= 1e-3;
radius += 0.001;
// define the cost function and its gradient
let df = |u: &[f64], grad: &mut [f64]| -> Result<(), SolverError> {
if a < 0.0 || b < 0.0 {
Err(SolverError::Cost)
} else {
rosenbrock_grad(a, b, u, grad);
Ok(())
}
};
let f = |u: &[f64], c: &mut f64| -> Result<(), SolverError> {
if a < 0.0 || b < 0.0 {
Err(SolverError::Cost)
} else {
*c = rosenbrock_cost(a, b, u);
Ok(())
}
};
// define the bounds at every iteration
let bounds = constraints::Ball2::new(None, radius);
// the problem definition is updated at every iteration
let problem = Problem::new(&bounds, df, f);
// updated instance of the solver
let mut panoc = PANOCOptimizer::new(problem, &mut panoc_cache)
.with_max_iter(max_iters);
let status = panoc.solve(&mut u).unwrap();
i += 1;
// print useful information
println!(
"parameters: (a={:.4}, b={:.4}, r={:.4}), iters = {}",
a, b, radius, status.iterations()
);
println!("u = {:#.6?}", u);
}
}The above function will print
parameters: (a=0.9990, b=101.0000, r=1.0010), iters = 43
u = [
0.786980,
0.618599
]
parameters: (a=0.9980, b=102.0100, r=1.0020), iters = 5
u = [
0.787543,
0.619500
]
parameters: (a=0.9970, b=103.0301, r=1.0030), iters = 5
u = [
0.788107,
0.620400
]
parameters: (a=0.9960, b=104.0604, r=1.0040), iters = 5
u = [
0.788670,
0.621301
]
...
In a real time implementation, it is important for our parametric optimizer to meet the maximum runtime requirements. For example, in digital optimization-based control or estimation applications, the computation needs to complete well within the sampling period of the system.
OpEn provides the method with_max_duration in PANOCOPtimizer that allows us to set
a maximum time duration, after which the solver stops. If it fails to converge because of
the imposition of a maximum allowed duration, the exit status will be
ExitStatus::NotConvergedOutOfTime.