Getting started¶

This guide introduces scaly's modelling workflow through a double-integrator control problem: symbolic expressions, functions, derivatives, optimization, and generated C. It assumes familiarity with Python, NumPy, and basic optimal control, but no experience with symbolic modelling libraries such as CasADi.

The examples use scaly with the IPOPT plugin and a C compiler, as described in Installation.

Symbolic variables and expressions¶

Consider the discrete-time model

\[ z_k = \begin{bmatrix}p_k \\ v_k\end{bmatrix}, \qquad z_{k+1} = f(z_k, u_k) = \begin{bmatrix}p_k + 0.1 v_k \\ v_k + 0.1 u_k\end{bmatrix}. \]

A symbolic variable represents an input whose numerical value is not yet specified. In scaly, it has a name and a fixed shape. For this model, z represents a two-element state vector and u a one-element control vector:

import numpy as np
import scaly as sc

z = sc.sym("z", 2)
u = sc.sym("u", 1)
znext = z + 0.1 * sc.concat([z[1:], u])

All three objects are instances of sc.Expr. z and u are input expressions and znext describes a calculation using them. Unlike an operation on NumPy arrays, this addition does not calculate a numerical result. It creates an expression node that records the addition and its operands.

Together, the inputs and operations form an expression graph. Its nodes represent values, and its edges record which values an operation needs. For znext, the graph records the velocity slice, its concatenation with u, the multiplication by 0.1, and the addition to z. Scaly uses this graph to calculate derivatives and generate code. Assigning a new value to the Python name z later does not change the recorded graph.

Shapes, indexing, broadcasting, and arithmetic follow NumPy conventions. Here, z[1:] has shape (1,), so concatenating it with u produces a vector of shape (2,). * is elementwise multiplication and @ is matrix multiplication. The expression reference lists Expr methods, and array builders include operations such as sc.concat, sc.stack, and sc.sumsqr.

From expressions to a Function¶

An sc.Expr describes a value. An sc.Function gives a calculation named inputs and outputs so that it can be evaluated, composed with other functions, differentiated, or exported as C. Function is the unit of composition in scaly.

The @sc.function decorator constructs one by running a Python body with symbolic inputs:

@sc.function(sc.G(sc.L("z", 2), sc.L("u", 1)), sc.L("znext", ...))
def model(inputs: tuple[sc.Expr, sc.Expr]) -> sc.Expr:
    z, u = inputs
    return z + 0.1 * sc.concat([z[1:], u])

L stands for leaf and declares one named array, G stands for group and combines leaves or other groups into a tree of inputs or outputs. Here, the input group is a tuple containing z and u. The output is one leaf named znext, whose shape is inferred from the returned expression because its declaration uses ....

The decorator creates the symbols, runs the body once, and records the returned expression. After decoration, model is an sc.Function object, rather than the original Python function. It remains callable:

z1 = model((np.array([1.0, 2.0]), np.array([0.5])))
print(z1)
# [1.2  2.05]

A call with numerical arrays evaluates the graph. The first numerical call generates and compiles C. Subsequent calls reuse the compiled code. The Python body does not run again. A call with symbolic expressions instead includes the function in a larger graph, as the trajectory example below demonstrates.

The input tree determines the call structure: sc.G introduces a tuple, whereas a single sc.L takes or returns an array directly. Shape () denotes a scalar, and shape (1,) denotes a one-element vector. These are distinct. The functions guide covers nested groups and other declarations.

Optional type annotations¶

The annotations on model describe the symbolic Python body: a tuple of two Expr inputs and one Expr output. They do not change tracing or numerical evaluation. The example also works without them.

Scaly's typing is designed to give you more checking as you provide more type information. With annotations, an IDE's type checker can check the body's argument and return types against the decorator's declared structure. For example, returning a tuple where the decorator declares one leaf is a type error. The decorated Function also carries symbolic and numerical input and output types, so calls with the wrong tuple structure can be flagged before execution. Array dimensions are checked at runtime as Python annotations cannot encode shapes.

Derivatives are Functions too¶

For the model above,

\[ \frac{\partial f}{\partial z} = \begin{bmatrix} 1 & 0.1 \\ 0 & 1 \end{bmatrix}. \]

sc.jacobian constructs an sc.Function for this derivative:

model_jac = sc.jacobian(model, "znext", "z")
print(model_jac((np.array([1.0, 2.0]), np.array([0.5]))))
# [[1.  0.1]
#  [0.  1. ]]

The strings select the output and input by their declared names. model_jac keeps the input structure of model, even though this particular Jacobian is constant. For a nonlinear model, the supplied values determine where the Jacobian is evaluated.

Scaly applies differentiation rules to the expression graph, producing another graph. It does not approximate derivatives by perturbing numerical inputs. The resulting function can be composed and exported like any other function. The derivatives guide covers gradients, Hessians, and products with derivative matrices. The sparsity guide covers calculating only entries that may be nonzero.

Composing a trajectory model¶

A function can reuse model to compute a trajectory and its cost:

\[ z_{k+1} = f(z_k,u_k), \qquad J(z_0, U) = \sum_{k=0}^{N-1}\left(\lVert z_k\rVert^2 + 0.1 u_k^2\right) + 10\lVert z_N\rVert^2, \qquad U=(u_0,\ldots,u_{N-1}), \quad N=20. \]

N = 20


@sc.function(
    sc.G(sc.L("z0", 2), sc.L("us", N)),
    sc.G(sc.L("zN", ...), sc.L("cost", ...)),
)
def rollout(inputs: tuple[sc.Expr, sc.Expr]) -> tuple[sc.Expr, sc.Expr]:
    z, us = inputs
    cost = sc.const(0.0)
    for k in range(N):
        u = us[k : k + 1]
        cost = cost + sc.sumsqr(z) + 0.1 * sc.sumsqr(u)
        z = model((z, u))
    return z, cost + 10.0 * sc.sumsqr(z)


z0 = np.array([1.0, 0.0])
us0 = np.zeros(N)
zN, cost = rollout((z0, us0))
print(zN, float(cost))
# [1. 0.] 30.0

sc.const(0.0) creates a constant expression, and sc.sumsqr(z) expresses \(\lVert z\rVert^2\). The symbolic call model((z, u)) records a call to the existing Function inside rollout. The two output leaves give the numerical result its tuple structure.

Python runs while the decorator builds the graph. Consequently, this for loop creates N successive calls to model. It does not record a loop. The horizon is fixed when rollout is defined, while z0 and us can change on every evaluation. This distinction between Python execution and recorded operations also matters for repeated independent calculations, covered in the bonus section on vmap.

Optimization problems and solvers¶

The trajectory model gives a single-shooting formulation with controls as decision variables and the measured initial state \(\bar z\) as a parameter:

\[ \begin{aligned} \min_U \quad & J(\bar z,U) \\ \text{subject to}\quad & z_0 = \bar z, \\ & z_{k+1}=f(z_k,u_k), && k=0,\ldots,N-1, \\ & z_N=0, \\ & -2 \le u_k \le 2, && k=0,\ldots,N-1. \end{aligned} \]

The recurrence is already built into rollout. sc.problem adds the objective, terminal equality, and control bounds:

@sc.problem(vars=sc.L("us", N), params=sc.L("z0", 2))
def control_problem(us: sc.Expr, z0: sc.Expr) -> sc.ProblemSpec[sc.Expr]:
    zN, cost = rollout((z0, us))
    return sc.ProblemSpec(
        minimize=cost,
        eq=(zN,),
        lb=sc.const(-2.0),
        ub=sc.const(2.0),
    )


solve = sc.solver(control_problem, "ipopt", options={"print_level": 0})

Inside control_problem, us and z0 are symbolic sc.Expr inputs. The cost, constraint expressions, and bounds returned in sc.ProblemSpec are also sc.Expr objects. The builder records their dependence on variables and parameters without evaluating them numerically.

The decorator creates an sc.Problem, which describes the optimization problem independently of the solver. vars declares what the solver may change, and params declares what stays fixed during each solve. Expressions in eq are constrained to zero. The scalar lb and ub bounds apply to every control.

sc.solver selects a backend and returns an sc.Solver. Scaly constructs the objective, constraint, and derivative functions that backend needs. The numerical call takes the problem parameters:

us_opt, lam_box, lam_eq, lam_ineq = solve(z0)
status = solve.stats().to_solver_status()
print(status.name)
assert status.ok

zN_opt, cost_opt = rollout((z0, us_opt))
print(np.round(zN_opt, 6))
# Approximately [0. 0.]

The result contains the controls and multipliers for variable bounds, equality constraints, and inequality constraints. lam_ineq is empty here because the problem has only variable bounds and equalities. Solver status is separate from these arrays and indicates whether the solve succeeded.

Initial variables and multipliers default to zero. solve(z0, x0=us0) supplies an initial control sequence, and warm= accepts a previous result for a warm start. The solver guide describes backend-specific warm-start behavior and the underlying solve.function, an sc.Function with explicit inputs for the initial variables and multipliers.

Generated C for deployment¶

The same model can be used from Python during development and exported for integration into a C or C++ application. For example, a ROS node can evaluate the dynamics or run an optimization solver without calling Python. Standalone model code can also be compiled with the toolchain used by an embedded target.

from pathlib import Path
from scaly.codegen import write_module

write_module(model, Path("generated"))

This writes generated/model.c and generated/model.h. Passing solve instead exports the solver and its model calculations. That code also needs the native IPOPT libraries. Exporting files ahead of the application build is the ahead-of-time (AOT) path, in contrast to just-in-time (JIT) compilation on the first Python call.

The code generation guide describes export options, compilation, and the generated C and C++ APIs, including argument buffers, working memory, and sparse output layouts.

Bonus: repeated stages with `vmap`¶

sc.vmap records repeated independent calls to a function as one mapped operation. The compiler can then emit a loop instead of a separate call site for every stage. A multiple-shooting formulation makes this useful for the same control problem by including states among the decision variables:

\[ \begin{aligned} \min_{Z,U}\quad & \sum_{k=0}^{N-1}\left(\lVert z_k\rVert^2+0.1u_k^2\right) + 10\lVert z_N\rVert^2 \\ \text{subject to}\quad & z_0=\bar z, \qquad z_N=0, \\ & f(z_k,u_k)-z_{k+1}=0, && k=0,\ldots,N-1, \\ & -2\le u_k\le2, && k=0,\ldots,N-1. \end{aligned} \]

Each dynamics residual, or defect, now takes candidate states as inputs. Its evaluation does not depend on the result of another defect evaluation:

@sc.function(
    sc.G(sc.L("z", 2), sc.L("u", 1), sc.L("znext", 2)),
    sc.L("defect", ...),
)
def defect(inputs: tuple[sc.Expr, sc.Expr, sc.Expr]) -> sc.Expr:
    z, u, znext = inputs
    return model((z, u)) - znext


@sc.problem(vars=sc.L("w", 3 * N + 2), params=sc.L("z0", 2))
def multiple_shooting(w: sc.Expr, z0: sc.Expr) -> sc.ProblemSpec[sc.Expr]:
    states = w[: 2 * (N + 1)]
    controls = w[2 * (N + 1) :]
    defects = sc.vmap(defect, N, {
        "z": states[:-2],
        "u": controls,
        "znext": states[2:],
    })
    return sc.ProblemSpec(
        minimize=sc.sumsqr(states[:-2]) + 0.1 * sc.sumsqr(controls)
                 + 10.0 * sc.sumsqr(states[-2:]),
        eq=(states[:2] - z0, defects, states[-2:]),
        ineq=(sc.bounded(controls, lo=-2.0, hi=2.0),),
    )

w stacks N + 1 two-element states followed by N controls. The mapping keys are the input names of defect. Scaly splits each supplied expression into N chunks of the corresponding input size: two elements for each state and one for each control. states[:-2] supplies stages 0 through N - 1, and states[2:] supplies stages 1 through N. sc.bounded expresses the control limits as inequalities within the larger decision vector.

The difference in generated structure is roughly as follows. These are sketches with simplified signatures, before any inlining or other compiler optimization:

// Single shooting: the Python loop creates N call sites.
model(z0, u0, z1);
model(z1, u1, z2);
/* ... */
model(z19, u19, z20);

// Multiple shooting: vmap keeps one loop body.
for (int k = 0; k < N; ++k) {
    defect(states[k], controls[k], states[k + 1], defects[k]);
}

The mapped stage code stays one loop body as the horizon grows, including in its derivatives. Numerical work and storage still grow with N, as can sparsity tables in the generated header. vmap cannot replace the sequential recurrence in rollout: it applies when the calls can be evaluated independently. The functions guide covers shared inputs and explicit slice mappings.