Usage¶

This page shows the intended workflow for both the C++ and Python interfaces. The solver concepts and setup order are the same in both languages; the tabs below only change the interface syntax.

For the shared mathematical background, see Concepts.

Main Types¶

Typical Workflow¶

1. Choose a parameterization¶

Factory methods are provided for

• Move-blocking (moveBlocking())
• Linear interpolation (linearInterp())
• Clamped B-splines (bspline())

When creating a Parameterization object, you can either specify how many control points to use (uniform knots by default, or geometrically spaced knots with knot_ratio) or you can supply a custom knot vector. For unclamped splines, you must use the regular constructor directly. Other bases, such as Laguerre functions, can be supplied as a sampled basis matrix with fromBasis() (see Custom Bases).

If you are unsure what to use, linear interpolation is a good default:

param = affine_mpc.Parameterization.linearInterp(horizon_steps, num_control_points)

auto param = affine_mpc::Parameterization::linearInterp(horizon_steps, num_control_points);

See Input Parameterization for more details.

2. Configure options if needed¶

opts = affine_mpc.Options(
    input_cost=affine_mpc.InputCostMode.ControlPoints,
    input_bounds=affine_mpc.InputBoundMode.ControlPoints,
    slew_initial_input=False,
    slew_bounds=affine_mpc.SlewBoundMode.ControlPoints,
    state_bounds=affine_mpc.StateBoundMode.Off,
)

affine_mpc::Options opts;
opts.input_cost = affine_mpc::InputCostMode::ControlPoints;
opts.input_bounds = affine_mpc::InputBoundMode::ControlPoints;
opts.slew_initial_input = false;
opts.slew_bounds = affine_mpc::SlewBoundMode::ControlPoints;
opts.state_bounds = affine_mpc::StateBoundMode::Off;

The options are grouped into costs (input_cost, use_input_rate_cost), input constraints (input_bounds, refined_bound_interval, slew_initial_input, slew_bounds), state constraints (state_bounds, bounded_states), and general constraints (num_general_constraints, num_terminal_constraints), in this order. In Python they are keyword arguments; in C++20, designated initializers must follow this order. input_bounds chooses how the input limits are enforced: on the control points (ControlPoints, the default), on the sampled inputs (Trajectory), or on the control points after knot insertion (Refined, whose longest knot interval in samples is refined_bound_interval); see Input Bounds. See Concepts for more details on options.

3. Construct the MPC object¶

The formulation chooses how the QP represents the predicted states; the interface is the same for both, and so is the solution (up to the solver tolerance).

• Condensed (Formulation.Condensed in Python, Formulation::Condensed in C++, the default): the states are eliminated with the model, so only the control points are QP variables; usually the best starting point, especially for shorter horizons and moderate dimensions
• Sparse (Formulation.Sparse, Formulation::Sparse): the predicted states are QP variables too, and the model is a set of equality constraints; useful when you want to preserve more explicit sparsity structure or work with larger problems, such as many control points over a long horizon

The best way to know which is faster for your problem is to try them both. In Python, formulation is a keyword-only argument.

With Parameterization:

mpc = affine_mpc.MPC(state_dim=n, input_dim=m, param=param, opts=opts)
# mpc = affine_mpc.MPC(n, m, param, opts, formulation=affine_mpc.Formulation.Sparse)

affine_mpc::MPC mpc(state_dim, input_dim, param, opts);
// affine_mpc::MPC mpc(state_dim, input_dim, param, opts, affine_mpc::Formulation::Sparse);
// affine_mpc::MPC mpc(state_dim, input_dim, param, affine_mpc::Formulation::Sparse);  // default options

No Parameterization:

mpc = affine_mpc.MPC(state_dim, input_dim, horizon_steps, opts)
# mpc = affine_mpc.MPC(state_dim, input_dim, horizon_steps, opts, formulation=affine_mpc.Formulation.Sparse)

affine_mpc::MPC mpc(state_dim, input_dim, horizon_steps, opts);
// affine_mpc::MPC mpc(state_dim, input_dim, horizon_steps, opts, affine_mpc::Formulation::Sparse);

4. Configure the MPC object¶

The order of these steps does not matter, but they all must happen before initializing the solver.

Set the model¶

Discrete:

mpc.setModelDiscrete(Ad, Bd, wd)

mpc.setModelDiscrete(Ad, Bd, wd);

Continuous:

mpc.setModelContinuous2Discrete(Ac, Bc, wc, dt)

mpc.setModelContinuous2Discrete(Ac, Bc, wc, dt);

The continuous model is discretized exactly under a zero-order hold on the input (the input is constant over each step of length dt). The matrix exponential is computed with scaling and squaring, so the result stays accurate even when the dynamics are fast relative to dt. An optional fifth argument, tol (default 1e-10), sets the truncation tolerance of the underlying Taylor series. The default gives results close to machine precision. If you update the affine approximation and rediscretize the model every time step, a looser value such as 1e-6 saves a few matrix products while keeping errors far below typical model error. If only the model offset changes between solves (with \(A\) and \(B\) fixed), set it with setModelOffsetDiscrete(wd) or setModelOffsetContinuous2Discrete(wc) instead, which avoids rebuilding and refactoring the QP (see Runtime Updates). In Python, model matrices written in place into Fortran-order arrays avoid a copy on every call (see Interface Notes).

The same discretization is available on its own as Discretizer, e.g. to simulate the plant or to discretize a model before passing it to several MPC objects:

disc = affine_mpc.Discretizer(state_dim, input_dim)
Ad, Bd, wd = disc.discretize(Ac, Bc, wc, dt)

affine_mpc::Discretizer disc{state_dim, input_dim};
disc.discretize(Ac, Bc, wc, dt, Ad, Bd, wd); // writes into Ad, Bd, wd

Set input limits¶

mpc.setInputLimits(u_min, u_max)

mpc.setInputLimits(u_min, u_max);

Set enabled optional parameters¶

if opts.slew_initial_input:
    mpc.setSlewRateInitial(u0_slew)
if opts.slew_initial_input or opts.use_input_rate_cost:
    mpc.setPreviousInput(u_prev)  # defaults to zeros
if opts.slew_bounds != affine_mpc.SlewBoundMode.Off:
    mpc.setSlewRate(u_slew)  # largest change per sample step
if opts.state_bounds != affine_mpc.StateBoundMode.Off:
    mpc.setStateLimits(x_min, x_max)  # one entry per bounded state
if opts.state_bounds == affine_mpc.StateBoundMode.Soft:
    mpc.setStateBoundPenalties(l1, l2)  # one entry per bounded state
if opts.num_general_constraints > 0:
    mpc.setGeneralConstraints(Hx, Hu, h_min, h_max)  # h_min <= Hx x_k + Hu u_k <= h_max
if opts.num_terminal_constraints > 0:
    mpc.setTerminalConstraints(Hf, hf_min, hf_max)  # hf_min <= Hf x_T <= hf_max

if (opts.slew_initial_input) {
  mpc.setSlewRateInitial(u0_slew);
}
if (opts.slew_initial_input || opts.use_input_rate_cost) {
  mpc.setPreviousInput(u_prev); // defaults to zeros
}
if (opts.slew_bounds != affine_mpc::SlewBoundMode::Off) {
  mpc.setSlewRate(u_slew); // largest change per sample step
}
if (opts.state_bounds != affine_mpc::StateBoundMode::Off) {
  mpc.setStateLimits(x_min, x_max); // one entry per bounded state
}
if (opts.state_bounds == affine_mpc::StateBoundMode::Soft) {
  mpc.setStateBoundPenalties(l1, l2); // one entry per bounded state
}
if (opts.num_general_constraints > 0) {
  // h_min <= Hx x_k + Hu u_k <= h_max
  mpc.setGeneralConstraints(Hx, Hu, h_min, h_max);
}
if (opts.num_terminal_constraints > 0) {
  mpc.setTerminalConstraints(Hf, hf_min, hf_max); // hf_min <= Hf x_T <= hf_max
}

To limit only some states, list them in the options, e.g. bounded_states=[3, 4] in Python or opts.bounded_states = {3, 4}; in C++ (see State Saturation).

For linear constraints on combinations of states and inputs, set the number of rows in the options, e.g. num_general_constraints=2 in Python or opts.num_general_constraints = 2; in C++, and pass \(H_x\) (\(2 \times n\)) and \(H_u\) (\(2 \times m\)) to setGeneralConstraints() (see General Linear Constraints).

Set weights¶

Without input cost:

mpc.setStateWeights(Q_diag, Qf_diag)
mpc.setStateWeights(Q_diag)  # Qf = Q

mpc.setStateWeights(Q_diag, Qf_diag);
mpc.setStateWeights(Q_diag); // Qf = Q

With input cost:

# set together
mpc.setWeights(Q_diag, Qf_diag, R_diag)
mpc.setWeights(Q_diag, R_diag)  # Qf = Q

# OR set individually
mpc.setStateWeights(Q_diag, Qf_diag)
mpc.setInputWeights(R_diag)

// set together
mpc.setWeights(Q_diag, Qf_diag, R_diag);
mpc.setWeights(Q_diag, R_diag); // Qf = Q

// OR set individually
mpc.setStateWeights(Q_diag, Qf_diag);
mpc.setInputWeights(R_diag);

With InputCostMode::ControlPoints, the input weight R weights the control points, not the sampled inputs (InputCostMode::Trajectory weights the sampled inputs). With non-uniform knots, a control point on a short early interval is weighted the same as one on a long late interval, even though it affects fewer samples. Keep this in mind when comparing input penalties across parameterizations.

With input rate cost (use_input_rate_cost), which penalizes changes of the sampled inputs instead:

mpc.setInputRateWeights(R_rate_diag)  # defaults to ones

mpc.setInputRateWeights(R_rate_diag); // defaults to ones

Set references¶

State step reference:

mpc.setReferenceState(x_step)

mpc.setReferenceState(x_step);

State trajectory reference: In Python, trajectories are arrays with one row per step, e.g. x_traj[k] is the reference for \(\bm{x}_{k+1}\) in a (horizon_steps, state_dim) array. In C++, they are matrices with one column per step, e.g. x_traj.col(k) in a state_dim x horizon_steps matrix. Both also accept the stacked layout (horizon_steps * state_dim values, step after step), which has the same memory as either one, so a whole matrix or array is copied in one pass.

mpc.setReferenceStateTrajectory(x_traj)  # (horizon_steps, state_dim) or flat

mpc.setReferenceStateTrajectory(x_traj); // state_dim x horizon_steps, or stacked

If input cost is enabled, input references can also be configured. setReferenceInput() sets a constant reference in both modes (with InputCostMode::ControlPoints and a custom basis it throws; use setReferenceInputControlPoints() there). setReferenceInputControlPoints() works only with InputCostMode::ControlPoints, and setReferenceInputTrajectory() only with InputCostMode::Trajectory (see Input Cost):

mpc.setReferenceInput(u_step)
# OR, with InputCostMode.ControlPoints, (num_control_points, input_dim)
mpc.setReferenceInputControlPoints(control_points)
# OR, with InputCostMode.Trajectory, (horizon_steps, input_dim)
mpc.setReferenceInputTrajectory(u_traj)

mpc.setReferenceInput(u_step);
// OR, with InputCostMode::ControlPoints, input_dim x num_control_points
mpc.setReferenceInputControlPoints(control_points);
// OR, with InputCostMode::Trajectory, input_dim x horizon_steps
mpc.setReferenceInputTrajectory(u_traj);

5. Initialize the solver¶

if not mpc.initializeSolver():
    raise RuntimeError("Failed to initialize solver")

if (!mpc.initializeSolver()) {
  throw std::runtime_error("Failed to initialize MPC solver");
}

You can pass an OSQPSettings object into initializeSolver() if you want to customize the default solver settings. It has the fields of OSQP's settings (see OSQP's documentation). In Python, OSQPSettings() gives the recommended settings and OSQPSettings.fromOSQPDefaults() OSQP's defaults; in C++, OSQPSolver::getRecommendedSettings() gives the recommended settings and OSQPSettings{} (or OSQPSolver::getDefaultSettings()) OSQP's defaults. getRecommendedSettings() leaves max_iter at OSQP's default (4000). Long horizons solved from a cold start (e.g. \(T = 500\) with degree 2) can need more, so raise it if you see MaxIterReached:

settings = affine_mpc.OSQPSettings()  # recommended settings
settings.max_iter = 20000
mpc.initializeSolver(settings)

affine_mpc::OSQPSettings settings =
    affine_mpc::OSQPSolver::getRecommendedSettings();
settings.max_iter = 20000;
if (!mpc.initializeSolver(settings)) { /* handle the error */ }

With scaling on (the default), OSQP computes its problem scaling in initializeSolver() and again whenever \(P\) or \(A_{\mathrm{qp}}\) change (a model, weight, or constraint matrix update). At those times the library gives OSQP a cost vector \(\bm{q}\) too small to affect the cost scaling, so that scaling depends only on \(P\), not on any particular \(x_0\) or reference. (OSQP treats \(\bm{q} = 0\) as \(\|\bm{q}\|_\infty = 1\), which would cap the cost scaling at 1 and slow OSQP down when the weights are small.)

Choosing the QP solver¶

OSQP is the default solver for both formulations. For the condensed formulation, you can use DAQP instead, a dense dual active-set solver, by passing DAQPSettings to initializeSolver():

settings = affine_mpc.DAQPSettings()  # DAQP's defaults
if not mpc.initializeSolver(settings):
    raise RuntimeError("Failed to initialize solver")

affine_mpc::DAQPSettings settings; // DAQP's defaults
if (!mpc.initializeSolver(settings)) { /* handle the error */ }

The Python wheels include DAQP. In C++ it is optional: build the library with -DAFFINE_MPC_WITH_DAQP=ON (an installed DAQP is used if CMake finds one, otherwise it is downloaded and built with the library). DAQPSettings and the initializeSolver() overload exist in every build, but without DAQP the overload throws std::logic_error (RuntimeError in Python); hasDAQP() tells whether the library has it. With the sparse formulation it throws std::invalid_argument (ValueError).

In C++, initializeSolver(std::unique_ptr<QPSolver>) also takes a solver object, such as std::make_unique<DAQPSolver>(settings) (from affine_mpc/daqp_solver.hpp) or any other implementation of the QPSolver interface.

Prefer DAQP for condensed problems of moderate size, where its active-set iterations return the solution to its tolerance (1e-6 on the constraints by default) instead of the approximate solution of OSQP's first-order iterations. In our measurements on this library's mass-spring-damper example (examples/sim.cpp), a solve() call with DAQP was more than ten times faster than with OSQP. DAQP's work grows with the cube of the number of QP variables (control points and slacks) and with the number of active constraints, so OSQP, or the sparse formulation, remains the choice for long horizons with many control points.

DAQP needs a positive definite cost matrix P for its standard iterations. When P is only positive semidefinite (for example with InputCostMode::Off and some states unweighted, or soft state bounds with only an L1 penalty), DAQP detects this while factorizing P and switches to proximal-point iterations, which still find a minimizer. DAQPSettings::eps_prox controls this: negative (the default) uses them only when needed, positive always, and zero never, in which case a singular P makes initializeSolver() return false (or solve() return UpdateFailed after a weight or model change).

resetWarmStart() empties DAQP's working set (the constraints it starts from), and SolveInfo::qp_solve_time is DAQP's own solve time, which excludes factorizing P after a weight or model change.

6. Solve in the control loop¶

Configuration misuse tends to throw exceptions, but solve() returns a SolveStatus type rather than throwing since it is in the main runtime loop. Common cases include:

• Success
• NotInitialized (solve() was called prior to initializeSolver())
• UpdateFailed (the solver rejected an update to the QP data, often from non-finite model or weight values; the QP was not solved, and the update is retried at the next solve, which usually recovers once the values are valid again)
• OSQP-derived failure conditions (see OSQP documentation)

How to respond to a failed solve depends on the application. For example, you could reuse the previously applied input or have a fallback control law like LQR. After a failed solve (any status other than Success or SolvedInaccurate), the result getters such as getNextInput() return the solver's last iterate, which is not a solution and may violate the constraints (after UpdateFailed they still return the previous solution). The previous input used by slew_initial_input and use_input_rate_cost is only updated after Success or SolvedInaccurate, so after a failure call setPreviousInput() with the input you actually applied. The failed solve's iterates may also have diverged, which would slow down the next solve since the solver warm starts from them, so call resetWarmStart() to start the next solve from zero.

MaxIterReached is a partial exception. With B-spline parameterizations its iterate is often close to the solution, but it is not guaranteed to satisfy the constraints. Applying getNextInput() (clamped to the input limits) can then be a reasonable fallback, but check it for your problem, and still call setPreviousInput() with the applied input.

With OSQP, SolvedInaccurate means it stopped at max_iter or its time limit with residuals within 10 times eps_abs and eps_rel, so treat it as a solve with looser tolerances; SolveInfo's residuals say how loose.

OSQP's first-order iterations are also less accurate when \(P\) is badly conditioned (its largest eigenvalue many orders of magnitude above its smallest), for example with nearly dependent custom basis functions, no input weight, or slow dynamics over a long horizon. Directions in which the cost barely changes are then weakly determined, so OSQP can stop with a first input that is noticeably off even when it reports Success, and warm starting can carry that error from step to step, which shows up in closed loop as input chatter or a small limit cycle. Tightening eps_abs and eps_rel, adding an input or input-rate weight, or using DAQP for a condensed problem fixes this; polishing does not help when no constraints are active. To check, compute the condition number of \(P\) from getQP() (e.g. np.linalg.cond(mpc.getQP().P.toarray()), outside the control loop).

status = mpc.solve(xk)
if status == affine_mpc.SolveStatus.Success:
    uk = mpc.getNextInput()
else:
    uk = u_fallback
    mpc.setPreviousInput(uk)  # if slew_initial_input or use_input_rate_cost is enabled
    mpc.resetWarmStart()

affine_mpc::SolveStatus status = mpc.solve(xk);
if (status == affine_mpc::SolveStatus::Success) {
  mpc.getNextInput(uk);
} else {
  uk = u_fallback;
  mpc.setPreviousInput(uk); // if slew_initial_input or use_input_rate_cost is enabled
  mpc.resetWarmStart();
}

7. Retrieve results in the control loop¶

To get the first optimized input (the next to apply):

uk = mpc.getNextInput()

Eigen::VectorXd uk(input_dim);
mpc.getNextInput(uk);

You can also access the full input and predicted state trajectories:

u_traj = mpc.getInputTrajectory()  # (horizon_steps, input_dim)
u_control_points = mpc.getInputControlPoints()  # (num_control_points, input_dim)
x_pred = mpc.getPredictedStateTrajectory()  # (horizon_steps, state_dim), x_1 first

Eigen::MatrixXd u_traj(input_dim, horizon_steps); // column k is u_k
mpc.getInputTrajectory(u_traj);

Eigen::MatrixXd u_control_points(input_dim, num_control_points);
mpc.getInputControlPoints(u_control_points);

Eigen::MatrixXd x_pred(state_dim, horizon_steps); // column k is x_{k+1}
mpc.getPredictedStateTrajectory(x_pred);

Solver diagnostics from the last solve (status, iterations, objective, residuals, and timing) are available as a SolveInfo:

info = mpc.getSolveInfo()
print(info.iterations, info.qp_solve_time, info.mpc_solve_time)

const affine_mpc::SolveInfo info = mpc.getSolveInfo();
std::cout << info.iterations << ' ' << info.qp_solve_time << ' '
          << info.mpc_solve_time << '\n';

objective is the QP objective, which omits constant terms of the MPC cost, so use it to compare solves with the same initial state and references rather than as the absolute MPC cost. qp_solve_time is the time the QP solver spent solving (for OSQP zero if it is built without profiling), and mpc_solve_time is the whole solve() call, so mpc_solve_time - qp_solve_time is the library's own overhead (QP assembly, data updates, and refactorization). Setters only store their values and the QP is updated in solve(), so mpc_solve_time covers essentially all of the MPC's work. The exception is setModelContinuous2Discrete(), which computes a matrix exponential inside the setter, so time it yourself if the model changes every step. If the model or weights change every step, the optional prepare() can do the matrix work before the next measurement (see Runtime Updates). prepare_time is the total time of the prepare() calls since the previous solve (zero if there were none), so prepare_time + mpc_solve_time is the MPC's work for the step while mpc_solve_time alone is its latency. primal_residual and dual_residual are the solver's absolute residuals (NaN when it does not report them, as with DAQP); for example, they tell how far a SolvedInaccurate solution is from optimal (see Accessing the QP).

8. (Optional) Logging¶

Create a logger bound to one MPC object:

logger = affine_mpc.MPCLogger(
    mpc,
    save_dir,
    ts,
    prediction_stride=1,
    mode=affine_mpc.MPCLogger.Mode.NpzCompressed,
)

affine_mpc::MPCLogger logger{&mpc, "/tmp/affine_mpc_example", dt, 1, false,
                             "log",
                             affine_mpc::MPCLogger::Mode::NpzCompressed};

Inside the control loop:

logger.logStep(t, xk)  # or logStep(t, xk, loop_time) to log your own timing

logger.logStep(t, xk); // or logStep(t, xk, loop_time) for your own timing

See Logging for output format and workflow details.

9. (Optional) Inspect the QP¶

getQP() returns the QP the last solve used, and assembleQP(x0, qp) assembles it for a given initial state without solving, for analysis or for handing to another solver. getQPVariableBlocks() and getQPConstraintBlocks() say which variables and rows are which. These are analysis and debugging tools: they copy or allocate, so keep them out of the control loop. getDualSolution() gives the multipliers of the rows from the last solve; it only fills the vector you pass in C++. See Accessing the QP.

Interface Notes¶

uk = mpc.getNextInput()  # most common
u_traj = mpc.getInputTrajectory()
u_control_points = mpc.getInputControlPoints()
x_pred = mpc.getPredictedStateTrajectory()

uk = np.empty(m)
mpc.getNextInput(uk)  # uk is also returned
x_pred = np.empty((T, n))  # or flat np.empty(T * n)
mpc.getPredictedStateTrajectory(x_pred)

Ad = np.zeros((n, n), order="F")  # allocate once, outside the loop
Bd = np.zeros((n, m), order="F")
wd = np.zeros(n)
while running:
    affine_model(x, u, Ad, Bd, wd)  # your code writes the entries in place, e.g. Ad[1, 0] = -k / mass
    mpc.setModelDiscrete(Ad, Bd, wd)

Eigen::VectorXd uk(input_dim);
mpc.getNextInput(uk);

Model Propagation Helper¶

You can propagate the internal model one step with:

x_next = mpc.propagateModel(xk, uk, x_next)

mpc.propagateModel(xk, uk, x_next);

Practical Guidance¶

• Do not call solve() before initializeSolver().
• If an option enables a constraint, call the corresponding setter before initialization.
• Fully configure the model, limits, weights, and references before calling initializeSolver().
• QP matrix sparsity is fixed after initialization, so later updates must not introduce new nonzero structure.
• Runtime updates to model terms and weights must preserve the initialized QP sparsity pattern.
• If a model coefficient may become nonzero later, initialize with that structure already present, e.g. with a representative model set before initializeSolver(); with the sparse formulation you can instead pass ModelSparsity::Full (see Model Entries).
• With the condensed formulation, the same holds for cost weights and constraint matrix entries (the sparse formulation tracks every one of them).
• To debug a suspected sparsity problem, call isWithinSparsityPattern() after solve() or prepare().
• If you enable slew_initial_input or use_input_rate_cost, provide the previous input before solving; after each successful solve, it is automatically set to the next input.