.. _getting started:
Getting started
===============
Most solvers in this library solve minimization problems of the following form:
.. math::
\newcommand\mymathbb[1]
{ {\rm I\mathchoice{\hspace{-2pt}}{\hspace{-2pt}}
{\hspace{-1.75pt}}{\hspace{-1.7pt}}#1} }
\newcommand{\Re}{\mymathbb R}
\begin{aligned}
& \underset{x}{\text{minimize}}
& & f(x) &&&& f : \Re^n \rightarrow \Re \\
& \text{subject to}
& & \underline{x} \le \phantom{g(}x\phantom{)} \le \overline{x} \\
&&& \underline{z} \le g(x) \le \overline{z} &&&& g : \Re^n \rightarrow \Re^m
\end{aligned}
The objective function :math:`f(x)` and the constraints function :math:`g(x)`
should have a Lipschitz-continuous gradient.
The two types of constraints are handled differently: the box constraints on
:math:`x` are taken into account directly, while the general constraints
:math:`g(x)` are relaxed using an Augmented Lagrangian Method (ALM).
Equality constraints can be expressed by setting
:math:`\underline{z}_i = \overline{z}_i`.
For ease of notation, the box constraints on :math:`x` and :math:`g(x)`
are represented as the inclusion in rectangular sets :math:`C` and :math:`D`
respectively:
.. math::
\begin{aligned}
& \underset{x}{\text{minimize}}
& & f(x) \\
& \text{subject to}
& & \phantom{g(}x\phantom{)} \in C \\
&&& g(x) \in D
\end{aligned}
Simple example
--------------
Consider the following simple two-dimensional problem, with the
Rosenbrock function (parametrized by :math:`p`) as the cost function:
.. math::
\begin{aligned}
& \underset{x_1,x_2}{\text{minimize}}
& & (1 - x_1)^2 + p\,(x_2 - x_1^2)^2 \\
&\text{subject to}
& & \begin{aligned}[t]
-0.25 &\le x_1 \le 1.5 \\
-0.5 &\le x_2 \le 2.5 \\
\end{aligned} \\
&&& \begin{aligned}[t]
(x_1 - 0.5)^3 - x_2 + 1 &\le 0 \\
x_1 + x_2 - 1.5 &\le 0 \\
\end{aligned}
\end{aligned}
In other words,
.. math::
\begin{aligned}
f(x) &= (1 - x_1)^2 + p\,(x_2 - x_1^2)^2 \\
C &= [-0.25, 1.5] \times [-0.5, 2.5] \\
g(x) &= \begin{pmatrix}
(x_1 - 0.5)^3 - x_2 + 1 \\
x_1 + x_2 - 1.5
\end{pmatrix} \\
D &= [-\infty, 0] \times [-\infty, 0]
\end{aligned}
Problem description
^^^^^^^^^^^^^^^^^^^
The objective function and the constraints are defined as
`CasADi `_ expressions with symbolic variables, and
then converted into CasADi functions. The arguments of these functions are the
decision variables and an optional parameter vector.
.. testcode::
# %% Build the problem (CasADi code, independent of alpaqa)
import casadi as cs
# Make symbolic decision variables
x1, x2 = cs.SX.sym("x1"), cs.SX.sym("x2")
x = cs.vertcat(x1, x2) # Collect decision variables into one vector
# Make a parameter symbol
p = cs.SX.sym("p")
# Objective function f and the constraints function g
f = (1 - x1) ** 2 + p * (x2 - x1**2) ** 2
g = cs.vertcat(
(x1 - 0.5) ** 3 - x2 + 1,
x1 + x2 - 1.5,
)
# Define the bounds
C = [-0.25, -0.5], [1.5, 2.5] # -0.25 <= x1 <= 1.5, -0.5 <= x2 <= 2.5
D = [-cs.inf, -cs.inf], [0, 0] # g1 <= 0, g2 <= 0
Next, we construct the alpaqa-specific minimization problem, using the
:py:func:`alpaqa.minimize` function.
The gradients of the problem functions are computed using CasADi, and they are
compiled as efficient C functions. All of this happens inside of the
:py:meth:`alpaqa.MinimizationProblemDescription.compile`
function, which returns an instance of
:py:class:`alpaqa.CasADiProblem` that can later be passed to a
solver.
.. testcode::
# %% Generate and compile C code for the objective and constraints using alpaqa
from alpaqa import minimize
problem = (
minimize(f, x) # Objective function f(x)
.subject_to_box(C) # Box constraints x ∊ C
.subject_to(g, D) # General ALM constraints g(x) ∊ D
.with_param(p, [1]) # Parameter with default value (can be changed later)
).compile()
.. testoutput::
:options: +ELLIPSIS
:hide:
...
Numerical values of the problem (like the bounds and the parameters) can be
specified when generating the problem, or can be modified after loading it:
.. testcode::
# You can change the bounds and parameters after loading the problem
problem.param = [10.0]
problem.general_bounds.lower[1] = -1e20
Selecting a solver
^^^^^^^^^^^^^^^^^^
The solvers in this package consist of an inner solver that can handle box
constraints, such as `PANOC `_,
and an outer ALM solver that relaxes the general constraints :math:`g(x) \in D`.
Solvers can be composed easily, for instance:
.. testcode::
# %% Build a solver with the default parameters
import alpaqa as pa
inner_solver = pa.PANOCSolver()
solver = pa.ALMSolver(inner_solver)
Each solver has its own set of optional parameters that can be specified using
a dictionary, for example:
.. testcode::
# %% Build a solver with custom parameters
inner_solver = pa.PANOCSolver(
panoc_params={
'max_iter': 1000,
'stop_crit': pa.PANOCStopCrit.FPRNorm,
},
lbfgs_params={
'memory': 10,
},
)
solver = pa.ALMSolver(
alm_params={
'tolerance': 1e-10,
'dual_tolerance': 1e-10,
'initial_penalty': 50,
'penalty_update_factor': 20,
},
inner_solver=inner_solver,
)
For a full overview and description of all parameters, see the
:ref:`Parameters cpp` page. For example, you can find the documentation for
:cpp:class:`alpaqa::PANOCParams`, :cpp:class:`alpaqa::LBFGSParams`, and
:cpp:class:`alpaqa::ALMParams`.
Some inner solvers can be configured with alternative fast directions. For
example, the PANOC solver uses the :cpp:class:`alpaqa::StructuredLBFGSDirection`
by default, but can also make use of e.g. :cpp:class:`alpaqa::LBFGSDirection`
or :cpp:class:`alpaqa::AndersonDirection`.
.. testcode::
# %% Build a solver with alternative fast directions
direction = pa.LBFGSDirection({'memory': 10})
inner_solver = pa.PANOCSolver({"stop_crit": pa.FPRNorm}, direction)
solver = pa.ALMSolver(
{
'tolerance': 1e-10,
'dual_tolerance': 1e-10,
'initial_penalty': 50,
'penalty_update_factor': 20,
},
inner_solver,
)
.. image:: ../img/classes-light.svg
:width: 100%
:alt: Different solver classes
:class: only-light
.. image:: ../img/classes-dark.svg
:width: 100%
:alt: Different solver classes
:class: only-dark
Solving the problem
^^^^^^^^^^^^^^^^^^^
Finally, you can obtain a solution by passing the problem specification to the
solver. This returns the local minimizer :math:`x_\mathrm{sol}`, the
corresponding Lagrange multipliers :math:`y_\mathrm{sol}` of the general
constraints :math:`g(x) \in D`, and a dictionary containing solver statistics.
.. testcode::
# %% Compute a solution
x_sol, y_sol, stats = solver(problem)
Optionally, you can supply an initial guess for both the decision variables
:math:`x` and the Lagrange multipliers :math:`y`. If no initial guess is
specified, the default initial values for :code:`x0` and :code:`y0` are zero.
.. testcode::
:hide:
import numpy as np
np.set_printoptions(precision=5) # make doctest predictable
.. testcode::
# %% Compute a solution
# Set initial guesses at arbitrary values
x0 = [0.1, 1.8] # decision variables
y0 = [0.0, 0.0] # Lagrange multipliers for g(x)
# Solve the problem
x_sol, y_sol, stats = solver(problem, x0, y0)
# Print the results
print(stats["status"])
print(f"Solution: {x_sol}")
print(f"Multipliers: {y_sol}")
print(f"Cost: {problem.eval_objective(x_sol):.5f}")
This will print something similar to:
.. testoutput::
SolverStatus.Converged
Solution: [-0.25 0.57813]
Multipliers: [10.3125 0. ]
Cost: 4.22119
It is highly recommended to always check the solver status in
:code:`stats["status"]` to make sure that it actually converged.
The :code:`stats` variable contains some other solver statistics as well, for
both the outer and the inner solver. You can find a full overview in the
documentation of :cpp:class:`alpaqa::ALMSolver::Stats`
and :cpp:class:`alpaqa::InnerStatsAccumulator\`.
.. figure:: ../img/example_minimal.svg
:width: 100%
:alt: Contour plot of the result
A contour plot of the objective function, with the constraints shown in
black, and the iterates produced by the solver in red. The red circle
indicates the optimal solution.