chc.control

Stability tier: stable.

Offline optimal control over a box: projected gradient, and a bound-constrained quasi-Newton.

The projected gradient also takes linear constraints on the whole action sequence (LinearConstraint -- a budget, a rate limit), projected onto together with the box so that every iterate stays feasible; see docs/adr/0001-linear-action-constraints.md.

Both minimise the same Bolza objective from the same discrete adjoint (chc.adjoint.control_gradient_adjoint); they differ only in how they use it. projected_gradient_control takes backtracked steepest-descent steps in variables scaled action by action, so it plans a problem the same in any units of its levers and of its cost (docs/adr/0059-the-planner-reads-the-problem-in-any-units.md). It is monotone and dependency-free but first-order: the scale is a diagonal, one number per action, so an instance whose actions act on the cost together still needs many steps -- thousands, not hundreds. Its steps is therefore a cap rather than a bill: the descent runs inside one compiled program that stops the moment the line search fails, so an unused step costs nothing and the default is loose enough for the stopping rule to decide. lbfgs_box_control hands the same gradient, in the same scaled variables, to SciPy's L-BFGS-B, which curves the step with a limited-memory secant approximation and reaches stationarity in about a hundred iterations or fewer rather than thousands -- but crosses the Python boundary on each one, so it cannot be compiled and is the reference rather than the workhorse. nlp_solver_certificate measures the gap rather than asserting it in prose.

Bound = float | Array module-attribute

One side of the action box: a scalar shared by every lever, or a per-lever array.

A real actuator set is rarely a cube -- a budget and a discount move in different units and over different ranges -- so a single scalar pair forces the caller to widen every lever to the loosest one, which is a larger feasible set than the plant has.

Blocks = tuple[tuple[Array, Array, Array, Array], ...] module-attribute

Constraint rows coloured into mutually orthogonal classes: (rows, lower, upper, |row|^2).

Duals = tuple[Array, tuple[Array, ...]] module-attribute

Dykstra's increments: one full vector for the box, one multiplier per constraint row.

SolverStatus = Literal['converged', 'max_iterations', 'no_progress'] module-attribute

Why the descent stopped -- reported, because "it stopped" and "it arrived" are different claims.

converged is the method's own stopping rule: the backtracking line search could not find a step that lowered the objective by more than tol times its value at the guess clipped to the box -- the cost in projected_gradient_solve, the cost and its penalties in chc.support.pessimistic_solve, which stops only where a step in the planner's own scale fails as well -- and the plan's first-order residual, in the descent's variables at the guess of zeros clipped to the box, is within 10 sqrt(2 max(tol, 4 eps)), the chc.plan.CausalPlan.stationarity_tolerance of the default tol, or no trial falls by the gradients at its ends either (_descend). Under linear constraints both also stop at a trial that moves the plan by the projection's error alone. max_iterations means the budget ran out first, so the answer is wherever the descent happened to be. no_progress means not one step was accepted, so the result is the caller's initial guess -- which is either an already-optimal guess or a badly scaled problem, and SolverResult.stationarity is what tells the two apart.

The status deliberately does not claim optimality. A stalled line search is a statement about steps, not about gradients, and on an ill-scaled instance the two come apart; the stationarity residual is returned beside it so the caller judges rather than trusts a label.

LinearConstraint dataclass

lower <= matrix @ us.ravel() <= upper: linear inequalities over the action sequence.

The actions are flattened row-major, so column k * m + j of matrix is lever j at step k. An absent side is -inf or inf, which lets one type state a half-space, a band or an equality: a total budget is a single row of ones with upper the budget, and a limit on how fast a lever may move is a band on consecutive differences (rate_limit). lower and upper may be scalars shared by every row.

Every iterate of the solve satisfies the constraints, not only the answer: the step is projected onto box ∩ constraints by Dykstra's algorithm, so a solve stopped by its step budget still hands back an executable plan. Why that, and not an augmented Lagrangian, is measured in docs/adr/0001-linear-action-constraints.md.

Raises:

Type	Description
ValueError	on a row of zeros (it constrains nothing, or nothing satisfies it, and either way it is a mistake), a non-finite coefficient, a NaN bound, or lower > upper.

rate_limit(horizon, caps) classmethod

|u[k+1, j] - u[k, j]| <= caps[j] between consecutive planned steps, per lever.

caps holds one entry per lever and inf leaves that lever free; only finite caps produce rows, so an all-inf call returns a constraint with none. The band is between planned steps only: the move from the action already applied to the first planned one is not constrained, because the solver is not told what that action was.

SolverResult dataclass

A finished box-constrained solve, with the evidence for how far to trust it.

SolverComparison dataclass

One instance of the box-constrained OC problem solved both ways.

NLPSolverCertificate dataclass

Evidence for where the first-order solver suffices and where it stops short.

broadcast_box(bound, shape, name, dtype)

Expand a scalar or per-lever bound to the full (horizon, m) action shape.

Accepts a scalar, a (m,) per-lever vector, or an already-full (horizon, m) schedule. A 1-D array is read as per-lever, never as per-step: for m == 1 a (horizon,) array would broadcast along the lever axis and silently constrain the wrong axis, so time-varying bounds have to be spelled out in two dimensions.

dtype is the actions' dtype, and it is required rather than defaulted: the box is validated in float64 but returned in the caller's precision, because a float64 box clipped against float32 actions promotes the answer and silently changes what the caller asked for.

Raises:

Type	Description
ValueError	on a shape that is neither, or on an inverted box. jnp.clip with lo > hi returns hi everywhere without complaint, which is a wrong answer rather than a failure, so the ordering is checked here where it can still be reported.

check_box(lo, hi)

Reject an inverted or empty box, or a nan bound, before jnp.clip turns it into a silent answer.

project_box(us, lo, hi)

Euclidean projection onto the box [lo, hi] (elementwise clip).

lo and hi may be scalars or per-lever arrays; both broadcast against us.

projected_gradient_solve(dyn, x0, us0, dt, cost, u_lo, u_hi, steps=10000, lr0=1.0, tol=1e-14, *, constraints=())

projected_gradient_control, returning why it stopped as well as where.

Same solve, same numbers; the tuple-returning function below is a thin wrapper on this one. It exists because the compact form cannot distinguish a descent that reached its stopping rule from one that ran out of budget, and a planner that acts on the second without knowing is acting on an unfinished solve. Costs one extra gradient evaluation, for the stationarity residual -- negligible against the ones the descent already spent.

lr0 is the line search's first step in the descent's scaled variables, and tol a fall of the cost relative to its value at the guess; projected_gradient_control says how.

constraints add linear rows over the whole sequence (LinearConstraint); the stationarity residual then projects onto the box and the rows, and SolverResult.constraint_violation reports how closely the rows hold. The residual is in the caller's units, so it is zero at a KKT point in any units but its size is not.

projected_gradient_control(dyn, x0, us0, dt, cost, u_lo, u_hi, steps=10000, lr0=1.0, tol=1e-14, *, constraints=())

Minimise J over the control sequence subject to box constraints.

Uses the discrete adjoint (chc.adjoint.control_gradient_adjoint) for the gradient and backtracking line search, so no accepted step raises the cost. Returns the optimised controls and the cost history (length 1 + accepted steps).

The step is scaled action by action, so a problem is planned the same in any units of each lever and of the cost. The descent runs in v_kj = sigma_kj u_kj: sigma_kj^2 is the cost's Gauss-Newton curvature along lever j at step k alone, over the cost at the guess clipped to the box. The box and the rows are carried into v and the projection is Euclidean there, so the descent stops at the problem's own KKT points.

• lr0 is the first step the line search tries, in v: 1 is the Newton step along each action alone. It is halved until the cost falls by a third of what the gradient predicts for the step, as a Newton step on a quadratic does; where no halving does, the longest step that lowers the cost is taken.
• tol is relative: a step counts where it lowers the cost by more than tol times the cost at the clipped guess, and from a guess that costs nothing any decrease counts. A cost mostly out of a plan's reach -- a target no plan comes near -- stops further from its optimum in what a plan can change, and may want a smaller tol. Where no step counts but the plan's first-order residual is beyond 10 sqrt(2 max(tol, 4 eps)), a step counts where the gradients at its ends read a fall and the cost does not rise (docs/adr/0097-a-plan-reads-its-first-order-check-on-its-face.md).

Both used to be absolute: lr0 a step of 0.2 action units per unit of gradient and tol a fall of 1e-9 in the cost's own units, so in units a thousand times a problem's own the descent ran out of steps (docs/adr/0059-the-planner-reads-the-problem-in-any-units.md). A fall eps in the cost's own units is tol = eps / |J|, with J the cost of the guess clipped to the box.

The descent runs inside a single compiled program (_projected_gradient_scan); only the trim to the accepted prefix happens on the host, so the return shape stays what a Python loop would have produced.

u_lo and u_hi are a scalar shared by every lever, a per-lever (m,) array, or a full (horizon, m) schedule. They are normalised to the action shape here and enter the compiled program as arrays, so a caller sweeping boxes compiles once rather than once per box value. dt, steps and the line-search scalars stay static to the compilation, the same convention chc.cost.total_cost already uses for dt.

constraints intersect the box with linear rows over the whole sequence (LinearConstraint), and every iterate stays inside both. The projection holds the rows to its tolerance, and a trial it moves off the plan by that alone is not a step: from a plan the rows and the box pin to one point, the descent stays on it.

Returns only where the descent landed. projected_gradient_solve returns why it stopped as well, which is what a caller needs before acting on the answer.

lbfgs_box_control(dyn, x0, us0, dt, cost, u_lo, u_hi, steps=300)

Minimise J over the control sequence subject to box constraints, by L-BFGS-B.

Same objective, same discrete-adjoint gradient and the same (controls, cost history) return as projected_gradient_control -- a drop-in whose step is curved by a limited-memory secant approximation instead of scaled by a backtracked constant. steps caps L-BFGS-B iterations, not gradient evaluations, and it converges well inside the default on the instances nlp_solver_certificate sweeps.

It minimises J over its value at the guess clipped to the box, in the planner's variables v = sigma * u read there (_action_units), and reads the answer back as the planner does (_actions_of), so it reads a problem the same in any units of its levers and of its cost. L-BFGS-B's stopping rules, a projected gradient under pgtol and a fall under ftol of the larger of the objective and 1, are pure numbers there; in the caller's units they are a gradient and a cost.

SciPy is the trust boundary: the objective and gradient cross it as float64 NumPy and the answer is cast back to us0's dtype, so a float32 caller is not silently promoted.

box_stationarity(dyn, x0, us, dt, cost, u_lo, u_hi)

First-order optimality residual ||u - P_box(u - grad J)|| -- zero exactly at a KKT point.

The convergence measure that does not depend on knowing the optimum: it is the natural stationarity map of a box-constrained problem, so it separates "the solver stopped" from "the solver arrived" without a reference solution and without a wall clock.

nlp_solver_certificate(horizon=40, dt=0.1, pg_steps=50, seed=3)

Sweep the control weight and confirm both halves of the claim about the first-order solver.

The instances differ only in R, which sets the conditioning of the reduced Hessian: a large control weight makes the objective well conditioned and steepest descent adequate, a small one makes it ill conditioned and a short step budget nowhere near enough. Asserting both is the point -- the certificate fails if the well-conditioned instance develops a gap (the comparison would be measuring something other than conditioning) and equally if the ill-conditioned one stops showing one (the reframing would no longer be justified).

pg_steps is deliberately far below the shipped default: the effect being exhibited is what a short first-order budget costs on an ill-conditioned instance. At the shipped cap the projected gradient ends below L-BFGS-B's cost on every instance, each on its own rule, the ill-conditioned one after 8 083 steps: by a relative 1.1e-9 to 1.0e-7, the certificate's relative_gap, or 4.4e-9 to 2.0e-7 of the cost. 50 steps is short of the 180 where the well-conditioned instance stops on its own rule, so that every instance shows a truncated descent. This measures the conditioning, not the library's behaviour.