MPC with generic cost and constraints

Generic MPC, sometimes called Economic MPC or simply Nonlinear MPC, refers to Model-Predictive Control with an arbitrary cost function, i.e., not restricted to quadratic cost functions.

The MPC problem type GenericMPCProblem allows much more general cost functions and constraints than the quadratic-programming based MPC problem types, and instances of GenericMPCProblem are solved using Optimization.jl. The examples in this section will all use the IPOPT solver, a good general-purpose solver that supports large, sparse problems with nonlinear constraints.

The GenericMPCProblem also supports integer constraints, or MINLP MPC (mixed-integer nonlinear programming).

Getting started

Below, we design an MPC controller for the nonlinear pendulum-on-a-cart system using the generic interface. Additional examples using this problem type are available in the tutorials

• MPC control of a Continuously Stirred Tank Reactor (CSTR)
• Solving optimal-control problems
• Solving optimal-control problems with MTK models
• Model-Predictive Control for the Research Civil Aircraft system
• MPC with binary or integer variables (MINLP)

The code for this example follows below, the code will be broken down in the following sections.

using DyadControlSystems, Plots
using DyadControlSystems.MPC
using DyadControlSystems.Symbolics
using StaticArrays
using LinearAlgebra

function cartpole(x, u, p, _=0)
    T = promote_type(eltype(x), eltype(u))
    mc, mp, l, g = 1.0, 0.2, 0.5, 9.81

    q  = x[SA[1, 2]]
    qd = x[SA[3, 4]]

    s = sin(q[2])
    c = cos(q[2])

    H = @SMatrix [mc+mp mp*l*c; mp*l*c mp*l^2]
    C = @SMatrix [0 -mp*qd[2]*l*s; 0 0]
    G = @SVector [0, mp * g * l * s]
    B = @SVector [1, 0]
    if T <: Symbolics.Num
        qdd = Matrix(-H) \ Vector(C * qd + G - B * u[1])
        return [qd; qdd]
    else
        qdd = -H \ (C * qd + G - B * u[1])
        return [qd; qdd]
    end
end

nu = 1   # number of controls
nx = 4   # number of states
Ts = 0.1 # sample time
N  = 10  # MPC optimization horizon
x0 = ones(nx)  # Initial state
r  = zeros(nx) # Reference

discrete_dynamics = MPC.rk4(cartpole, Ts)   # Discretize the dynamics
measurement = (x,u,p,t) -> x                # The entire state is available for measurement
dynamics = FunctionSystem(discrete_dynamics, measurement, Ts; x=:x^nx, u=:u^nu, y=:y^nx)

# Create objective
Q1 = Diagonal(@SVector ones(nx))    # state cost matrix
Q2 = 0.1Diagonal(@SVector ones(nu)) # control cost matrix
Q3 = Q2
QN, _ = MPC.calc_QN_AB(Q1, Q2, Q3, dynamics, r) # Compute terminal cost
QN = SMatrix{nx,nx}(QN)

p = (; Q1, Q2, Q3, QN) # Parameter vector

running_cost = StageCost() do si, p, t
    Q1, Q2 = p.Q1, p.Q2 # Access parameters from p
    e = (si.x)
    u = (si.u)
    dot(e, Q1, e) + dot(u, Q2, u)
end

difference_cost = DifferenceCost() do e, p, t
    dot(e, p.Q3, e)
end

terminal_cost = TerminalCost() do ti, p, t
    e = ti.x
    e'p.QN*e
end

objective = Objective(running_cost, terminal_cost, difference_cost)

# Create objective input
x = zeros(nx, N+1)
u = zeros(nu, N)
x, u = MPC.rollout(dynamics, x0, u, p, 0)
oi = ObjectiveInput(x, u, r)

# Create constraints
control_and_state_constraint = StageConstraint([-3, -4], [3, 4]) do si, p, t
    u = (si.u)[]
    x4 = (si.x)[4]
    SA[
        u
        x4
    ]
end

# Create observer, solver and problem
observer = StateFeedback(dynamics, x0)

solver = MPC.IpoptSolver(;
        verbose = false,
        tol = 1e-4,
        acceptable_tol = 1e-1,
        max_iter = 100,
        max_cpu_time = 10.0,
        max_wall_time = 10.0,
        constr_viol_tol = 1e-4,
        acceptable_constr_viol_tol = 1e-1,
        acceptable_iter = 2,
    )

prob = GenericMPCProblem(
    dynamics;
    N,
    observer,
    objective,
    constraints = [control_and_state_constraint],
    p,
    objective_input = oi,
    solver,
    xr = r,
    presolve = true,
);

# Run MPC controller
history = MPC.solve(prob; x0, T = 100, verbose = false)

# Extract matrices
X,E,R,U,Y = reduce(hcat, history)

plot(history)

Specifying cost and constraints

The GenericMPCProblem requires the specification of an Objective, which internally contains one or many cost functions, such as

• StageCost
• DifferenceCost
• TerminalCost
• TrajectoryCost

Each cost object takes a function as its first argument that computes the cost based on the relevant optimization variables. We illustrate with an example where we create a stage cost that computes $x(t)^T Q_1 x(t) + u(t)^T Q_2 u(t)$:

running_cost = StageCost() do si, p, t
    x = si.x
    u = si.u
    dot(x, Q1, x) + dot(u, Q2, u)
end

This uses the Julia do-syntax to create an anonymous function that takes the tuple (si, p, t). si is of type MPC.StageInput, a structure containing vectors x, u, r, all at the stage time t. A stage refers to a single instant in time in the optimization horizon. While most cost and constraint types passes a MPC.StageInput as the first argument to the cost/constraint function, the TrajectoryCost and TrajectoryConstraint are passed an MPC.ObjectiveInput which contains the entire trajectories $x \in \mathbb{R}^{n_x \times N+1}$ and $u \in \mathbb{R}^{n_u \times N}$.

One or many cost functions are finally packaged into an Objective, illustrated in the comprehensive example below.

Constraints are similarly defined to take a function that computes the constrained output as first argument. We illustrate with an example that creates control and state constraints corresponding to

\[\begin{aligned} -3 &\leq u_1(t) \leq 3 \\ -4 &\leq x_2(t) + x_4(t) \leq 4 \end{aligned}\]

using StaticArrays
control_and_state_constraint = StageConstraint([-3, -4], [3, 4]) do si, p, t
    u = si.u[1]
    x2 = si.x[2]
    x4 = si.x[4]
    SA[
        u
        x2 + x4
    ]
end

The full signature of StageConstraint is

StageConstraint(fun, lb, ub)

where fun is a function from (stage_input, parameters, time) to constrained output and lb, ub are the lower and upper bounds of the constrained output. In this case, our constrained output is a static array (for high performance) containing [u[1], x[2]+x[4]], which are the expressions we wanted to constrain, i.e., the constrained output.

If simple bounds on states and control inputs are desired, the function BoundsConstraint can be used instead.

The available constraint types are

• BoundsConstraint
• StageConstraint
• TerminalStateConstraint
• TrajectoryConstraint
• Integer and binary constraints are handled using the int_x and int_u keywords to GenericMPCProblem, see MPC with binary or integer variables (MINLP) for an example.

Specifying the discretization

The GenericMPCProblem allows the user to select the discretization method used in the transcription from continuous to discrete time. The available choices are

• MultipleShooting (the default if none is chosen)
• CollocationFinE
• Trapezoidal

More details on these choices are available under Discretization.

ModelingToolkit Interface

For users modeling systems with ModelingToolkit.jl, DyadControlSystems provides a streamlined interface to create MPC problems directly from symbolic MTK models. This approach offers several advantages:

• Cleaner code: Express costs and constraints symbolically without manual function coding
• Automatic differentiation: MTK handles differentiation through the symbolic expressions
• Symbolic simplification: MTK can optimize expressions before code generation
• Less error-prone: Symbolic expressions are easier to read and verify

Creating an MPC Problem from an MTK Model

The workflow consists of two simple steps:

1. Define costs and constraints symbolically using symbolic expressions
2. Construct the GenericMPCProblem directly from the MTK model

Example: Double Integrator with MTK

using DyadControlSystems, DyadControlSystems.MPC, Plots
using ModelingToolkit
using LinearAlgebra
t = ModelingToolkit.t_nounits
D = ModelingToolkit.D_nounits

# Step 1: Define the MTK model
vars = @variables x(t) [description="Position"] v(t) [description="Velocity"] u(t) [description="Control input"] rx(t) [description="Reference position"] rv(t) [description="Reference velocity"] e(t)[1:2] [description="Tracking error"]
e = collect(e)

@named model = ODESystem([
    D(x) ~ v
    D(v) ~ u
    e ~ [rx-x, rv-v]
], t, vars, [])

# Step 2: Define costs symbolically
Q1 = 1.0I(2)   # State weight
Q2 = 0.1I(1)   # Control weight
QN = 10.0I(2)  # Terminal weight

stage_cost = SymbolicStageCost(dot(e, Q1, e) + dot(u, Q2, u))
terminal_cost = SymbolicTerminalCost(dot(e, QN, e))
objective = Objective(stage_cost, terminal_cost)

# Step 3: Define constraints symbolically
N = 10
u_constraint = SymbolicStageConstraint(u, [-1.0], [1.0], N)  # -1 ≤ u ≤ 1

# Step 4: Create the MPC problem directly from MTK model
xr = [0.1, 0]  # Reference: appearing in e ~ [rx-x, rv-v]
solver = MPC.UnoSolver(verbose=false)

prob = GenericMPCProblem(model;
    N=N,
    Ts=0.1,
    inputs=[u],                   # Control inputs
    disturbance_inputs=[rx, rv],  # Reference/disturbance signals. This array determines the order in which the references appear in `xr`
    objective=objective,
    constraints=[u_constraint],
    xr=xr,
    solver=solver,
)

# Step 5: Simulate
history = MPC.solve(prob; x0=[1.0, 0.5], T=60, verbose=false)
plot(history, idxs=[model.x, model.v], layout=2, sp=1)
plot!(history, idxs=[model.rx, model.rv], sp=1, l=:dash, c=[1 2])
plot!(history, idxs=[model.u], sp=2)

using Test
@test history[model.e][end] ≈ zeros(2) atol=0.03

Test Passed

Symbolic Cost Types

All standard cost types support symbolic construction:

StageCost - Cost at each time step:

stage_cost = SymbolicStageCost(x^2 + v^2 + 0.1*u^2)

TerminalCost - Cost at the end of the horizon:

terminal_cost = SymbolicTerminalCost(10.0 * (x^2 + v^2))

DifferenceCost - Penalize changes over time:

# Penalize Δu = u(k) - u(k-1)
Q3 = 0.01
difference_cost = SymbolicDifferenceCost(u) do Δu, p, t
    Q3 * sum(abs2, Δu)
end

Symbolic Constraint Types

StageConstraint - Constraints at each time step:

# Single constraint
u_constraint = SymbolicStageConstraint(u, [-1], [1], N)

# Multiple constraints
multi_constraint = SymbolicStageConstraint([u, x], [-1, -5], [1, 5], N)

# Nonlinear constraints
nonlinear_constraint = SymbolicStageConstraint(x^2 + v^2, [0], [4], N)

TerminalStateConstraint - Terminal state constraints:

# Force terminal state to origin
terminal_constraint = SymbolicTerminalStateConstraint([x, v], [0, 0], [0, 0])

DifferenceConstraint - Bound rate of change:

# Limit control input rate: -0.5 ≤ Δu ≤ 0.5
rate_constraint = SymbolicDifferenceConstraint(u, [-0.5], [0.5], N) do Δu, p, t
    Δu  # Identity function for simple box constraints on the difference
end

Comparison: Manual vs. MTK Interface

Manual approach (traditional):

running_cost = StageCost() do si, p, t
    x, v = si.x
    u = si.u[1]
    x^2 + v^2 + 0.1*u^2
end

MTK approach (symbolic):

stage_cost = StageCost(x^2 + v^2 + 0.1*u^2)

The MTK approach is more concise and the symbolic expression is easier to verify for correctness. Both approaches produce identical results. The key difference is that with the MTK interface, you don't need to manually extract variables from the StageInput object - you work directly with symbolic variables.

How It Works

The MTK interface uses a lazy evaluation pattern:

1. When you call SymbolicStageCost(expr), it returns a SymbolicStageCost that stores the expression
2. When you create the GenericMPCProblem, it converts symbolic types to actual function-based types
3. This conversion uses build_explicit_observed_function from MTK to generate optimized code
4. The generated functions are then wrapped to match the expected MPC signatures

This design provides:

• Clean API: No configuration objects to pass around
• Type safety: Multiple dispatch ensures correct conversion
• Extensibility: Easy to add new cost/constraint types
• Performance: Generated code is as efficient as hand-written functions

Discretization Methods

The MTK interface supports all three discretization methods through the discretization parameter:

MultipleShooting (default):

prob = GenericMPCProblem(model; ..., discretization=MultipleShooting)

• Uses explicit RK4 integration
• Control inputs are piecewise constant (zero-order hold)
• Best for non-stiff ODE systems
• Fewest optimization variables

Trapezoidal:

prob = GenericMPCProblem(model; ..., discretization=Trapezoidal)

• Uses implicit trapezoidal integration
• Control inputs are piecewise affine (first-order hold)
• Supports both ODEs and DAEs
• More accurate for stiff systems

CollocationFinE:

prob = GenericMPCProblem(model; ...,
    discretization=CollocationFinE,
    discretization_kwargs=(n_colloc=5, roots_c="Legendre")
)

• Uses orthogonal collocation on finite elements
• Supports both ODEs and DAEs
• Most accurate for smooth trajectories
• More optimization variables but higher-order accuracy
• Use discretization_kwargs to specify collocation points and method

Notes and Limitations

• State variables are automatically determined by MTK; you cannot manually specify them
• Observer: The MTK interface uses state feedback with RK4 discretization for state prediction
• Scaling: The scale_x and scale_u parameters are not yet supported with the MTK interface
• Parameters: System parameters defined in the MTK model are automatically passed through to the generated functions
• Mixed usage: You can mix symbolic and manually-defined costs/constraints in the same problem
• Discretization: The observer always uses RK4 discretization regardless of the main discretization method