Constrained transcription (pybnf.transcription)

A transcription restates the fit PyBNF was asked to run as a larger problem with internal auxiliary variables and equality constraints that tie them back together. At the solution the two problems coincide; on the way there the enlarged one can be far better conditioned. Multiple shooting is the first consumer (issue #563): split an experiment at knots, introduce segment-start states z_j, simulate each segment independently, and impose continuity c_j = Phi_j(z_j, theta) - z_{j+1} = 0. Direct collocation, latent-state estimation, and path constraints are the same shape.

The pybnf.transcription package is that shape and nothing else. It holds no dynamics, makes no simulator call, defines no configuration key and no job_type; its only dependency inside PyBNF is pybnf.printing, for the exception base class. A consumer supplies two methods — score-and-differentiate the objective at an augmented point, and linearise the equality constraints there — plus a certification hook, and gets the augmented-Lagrangian outer loop, the penalty schedule, the transcription homotopy, best-iterate certification, and the reporting.

The design and the measurements behind its defaults are recorded in ADR-0109.

Errors

Constrained-transcription exceptions (#563).

A single dependency-free home for the layer’s exception, mirroring pybnf.gradient.errors: the layout, the equality interface, and the outer loop all raise it, and a consumer catches it without importing any of them.

exception pybnf.transcription.errors.TranscriptionError(log_message, user_message=None, hint=None)

A constrained transcription could not be built, stepped, or reconciled.

Raised for the structural faults of the layer – a variable block whose name collides with a reported free parameter, a Jacobian block outside its declared shape, a non-positive constraint scale, a penalty schedule that cannot start. These are consumer wiring errors, not bad points to be stepped over: an augmented layout that does not describe what the residual assembly is about to write into is silently wrong in a way no fit result would reveal.

It is a PybnfError so the message reaches the user intact when a fit is driven from a configuration rather than from a test.

Augmented variable layout

The augmented variable layout: reported free parameters plus internal auxiliary blocks (#563).

A constrained transcription solves a different problem from the one the user asked about. Multiple shooting (#563’s first consumer) splits an experiment at knots and adds one segment-start state z_j per knot, so the decision vector grows from the fit’s k free parameters to k + sum_j dim(z_j). Direct collocation would add a state per collocation node; latent-state estimation adds a state per unobserved species. In every case the added coordinates are internal to the transcription: they are searched, they are bounded, they carry gradient columns – and they are not biological fit parameters. Reporting them in sorted_params_*.txt would claim the fit estimated 3x as many quantities as it did, and would put a quantity with no scientific meaning next to ones that have it.

AugmentedLayout is the bookkeeping that keeps those two populations apart while letting one flat vector carry both. It owns exactly one thing: the map between

u_aug = [ u_reported | z_1 | z_2 | ... | z_K ]

and its named parts. The reported block is always first and always contiguous, so u_aug[:n_reported] is the vector every existing PyBNF seam already understands (a PSet’s coordinates in sampling space, ADR-0029) with no slicing ceremony – and a consumer that forgets to unpack gets the reported parameters, not a silently misaligned mixture.

Space

Every coordinate is in the space the optimizer walks. For the reported block that is the free parameters’ sampling space u (log10(theta) for a logvar), exactly as trf / lbfgs / gntr step in; for an internal block it is whatever space that block declares its bounds in. The layout does not transform anything – the d theta/d u chain rule stays where it already lives, in pybnf.gradient.assembly, applied once when the Jacobian is built. A block that wants to be searched in log space says so by being built in log space.

The homotopy seam

The #563 prototype’s central finding is that the segment-count homotopy is the mechanism, not a refinement to add later (issue #563, finding 5.2): coarsening 4 -> 2 -> 1 is what converts a segmented stage that scores worse than a flat line into a solve. A homotopy is a sequence of transcriptions of the same fit, so the layer needs a way to carry a point from one layout to the next. AugmentedLayout.carry_over() is that: the reported block always survives (it is the same fit), an internal block survives iff the next layout still declares a block of that name and size, and a block the next layout adds is seeded from its own VariableBlock.initial. Matching is by name, which is what makes the rule generic – the layout never learns what a knot is.

class pybnf.transcription.layout.AugmentedLayout(reported_names, lower, upper, blocks=())

The map between a flat augmented vector and its reported / internal parts.

Parameters:

• reported_names – The fit’s free parameters, in the order every existing PyBNF seam already uses (Configuration.variables). These occupy the leading, contiguous slice of every augmented vector.

• lower – Reported-block lower bounds, in sampling space.

• upper – Reported-block upper bounds, in sampling space.

• blocks – The internal VariableBlocks, in the order they are laid out after the reported block.

The layout is immutable and cheap to build, so a homotopy builds one per stage.

block(block_name)

The VariableBlock named block_name.

carry_over(u, target)

Move a point from this layout into target’s – one step of the homotopy.

The reported block always survives: it is the same fit, and its value is the whole reason the previous stage ran. An internal block survives iff target declares a block of the same name and size; a block target adds is seeded from its own VariableBlock.initial; a block target dropped is discarded (that is what coarsening is). Carried values are clipped into the target block’s box, since two stages need not bound an auxiliary variable identically.

A name that matches with a different size is a consumer bug – two stages disagree about what that block means – and raises rather than being silently reseeded.

describe()

One line for the run log: how many coordinates the transcription added, and where.

embed_gradient(gradient)

Zero-pad a reported-space gradient into augmented space.

For the corner where a term genuinely has no dependence on the auxiliary variables – a prior, a parameter-only penalty. A term that does depend on them (the data fit of a multiple-shooting segment, which reads z_j through the IC route) must be assembled in augmented space directly, not embedded.

embed_jacobian(jacobian)

Zero-pad a reported-space (m, k) Jacobian’s columns into augmented space.

initial_point(reported)

The augmented start point: reported as given, every internal block at its declared VariableBlock.initial.

internal_of(u, block_name)

Just block block_name’s values.

is_internal(index)

Whether coordinate index is an internal auxiliary variable rather than a reported free parameter – the predicate any reporting path filters on.

property lower

Stacked lower bounds over the whole augmented vector.

property n_internal

Number of internal auxiliary coordinates the transcription added.

property n_reported

Number of reported free parameters – the fit’s own k.

property names

Every coordinate’s name: the reported free parameters, then each block’s qualified component names. Guaranteed unique, and guaranteed to mark which coordinates are internal (they alone contain QUALIFIER).

pack(reported, internals=None)

Build an augmented vector from the reported parameters and a {block name: values} mapping. Every declared block must be supplied.

reported_of(u)

Just the reported free parameters – the only part of the vector that is a fit result. Every reporting, certification, and PSet path goes through this.

property reported_slice

The leading slice every augmented vector carries the reported parameters in.

property size

Length of an augmented vector.

slice_of(block_name)

The slice block_name occupies, raising rather than returning a wrong one.

unpack(u)

Split an augmented vector into (reported, {block name: values}).

property upper

Stacked upper bounds over the whole augmented vector.

pybnf.transcription.layout.QUALIFIER = '::'

Separator between an internal block’s name and one component’s label in a qualified name ('seg2::A_state'). Chosen because no PyBNF free-parameter name can contain it, so a qualified internal name can never be mistaken for – or collide with – a reported one.

class pybnf.transcription.layout.VariableBlock(name, labels, lower, upper, initial)

One named group of internal auxiliary variables.

Parameters:

• name – The block’s identity. Unique within a layout, and the key AugmentedLayout.carry_over() matches on across a homotopy stage change.

• labels – One label per component, for diagnostics and defect reports (for multiple shooting: the state names, e.g. ('A_state', 'Y_state', 'Z_state')).

• lower – Per-component lower bounds, in the block’s own space.

• upper – Per-component upper bounds.

• initial – The value a layout that newly introduces this block starts it at – the consumer’s best guess for the auxiliary variable (for multiple shooting, the state read off a nominal trajectory at the knot). Also what AugmentedLayout.initial_point() seeds.

Bounds are part of the block because the inner optimizers this layer feeds are bound-constrained (the Coleman-Li reflective step in trf / gntr), and a segment-start concentration that is allowed to go negative is not a state the simulator can restart from.

clipped(values)

values projected into this block’s box – what a consumer applies after a AugmentedLayout.carry_over() whose source stage had looser bounds.

property qualified_names

This block’s components as '<block>::<label>' – the names that appear in a defect report or a diagnostic, and that are guaranteed disjoint from every reported free-parameter name.

Equality residuals and Jacobians

The equality residual / Jacobian interface (#563).

What a constrained transcription adds to a fit is a vector of equality constraints c(u_aug) = 0 on the augmented variables. For multiple shooting they are the continuity defects c_j = Phi_j(z_j, theta) - z_{j+1}; for direct collocation they would be the per-node collocation equations; for latent-state estimation, the state-transition residuals. This module fixes the interface all of them present to the outer loop, and nothing else – it contains no dynamics, no simulator, and no optimizer.

Three decisions are worth naming.

The Jacobian is block-sparse, not dense. A transcription’s constraint Jacobian is mostly zeros with a strong, known structure: a continuity row for segment j reads theta, z_j, and z_{j+1} and nothing else, so dc/du is (constraint group x variable block) blocks against a background of exact structural zeros. BlockJacobian stores exactly those blocks – one dense (rows x cols) array per non-zero region, both ranges contiguous because AugmentedLayout lays every block out contiguously – and implements matvec(), rmatvec(), and gram() block-wise. Dense assembly (to_dense()) exists because today’s inner optimizers consume dense linear algebra (gntr eigen-decomposes its Hessian), but the structure is preserved rather than discarded on the way in, which is what leaves room for the condensing seam: eliminating the z block-by-block to recover a dense system of the fit’s own dimension k instead of k + sum_j dim(z_j). That is the standard multiple-shooting condensation and the reason the representation is block-structured now, before there is a model big enough to need it. Nothing in this module assumes condensing exists; nothing in it prevents adding it.

Blocks accumulate. Two blocks covering the same region add, exactly as pybnf.gradient.routing.route_experiment() folds two chain-rule paths reaching one sensitivity column. Every operation here is linear in the block list – to_dense, matvec, rmatvec, gram – so additive is the only semantics that makes them agree with each other, and a consumer that reaches one region by two paths writes two blocks rather than pre-summing them.

The constraints are scaled, and the outer loop sees only the scaled ones. A continuity defect is a difference of states, so its units are the state’s, and a model whose species span six orders of magnitude would otherwise hand the penalty term a condition number for free. Each constraint carries a strictly positive scale (s_i, typically the state’s own magnitude), and EqualityModel exposes c_i / s_i. One penalty parameter then means the same thing for every constraint, the feasibility tolerance is dimensionless, and the defect report the issue asks for is comparable across states. Scaling a constraint is an exact reparameterisation – lambda absorbs s – so nothing downstream has to know it happened. It matters most in the corner the #563 thread flags as the hard part of the motivating problem: with one observed state of three, the unobserved segment-start states are determined by continuity alone, so the conditioning of the constraint block is the conditioning of the whole inner problem.

class pybnf.transcription.equality.BlockJacobian(shape, blocks)

A constraint Jacobian dc/du_aug held as its non-zero blocks.

Parameters:

• shape – (n_constraints, layout.size).

• blocks – The JacobianBlocks. Regions that overlap add.

Everything outside a block is a structural zero – not a small number, a zero – which is what makes the block-wise operations exact rather than approximate.

property density

Stored fraction of the dense shape. The number a condensing seam would act on.

gram()

J.T @ J – the Gauss-Newton curvature of the penalty term rho/2 ||c||^2.

Accumulated over pairs of blocks whose row ranges intersect, since only those contribute: two blocks on disjoint rows are orthogonal by construction. For the block-diagonal-plus-arrowhead structure a transcription actually produces, almost every pair is disjoint.

matvec(v)

J @ v, touching only the stored blocks.

property nnz

Stored entries – the structural zeros are not among them.

rmatvec(y)

J.T @ y – how the constraint gradient J^T (lambda + rho c) is formed.

scaled(row_scales)

This Jacobian with row i divided by row_scales[i] – the Jacobian of the scaled constraints c_i / s_i.

to_dense()

The dense (m, n) array. For the inner optimizers that consume dense linear algebra; the block structure is kept on the way in so a future condensation does not have to rediscover it.

class pybnf.transcription.equality.EqualityModel(residual, jacobian, scales=None, names=None)

The equality constraints linearised at one augmented point.

Parameters:

• residual – The raw defect c(u), in the constraints’ own units.

• jacobian – dc/du_aug as a BlockJacobian.

• scales – Strictly positive per-constraint scales; None means all ones.

• names – Per-constraint labels for the defect report (for multiple shooting, '<knot block>::<state>').

The outer loop reads scaled_residual and scaled_jacobian, never the raw pair – see the module docstring. The raw pair is kept because it is what the consumer reports in the model’s own units when a user asks how far from continuous the fit was.

property defect_norm

max_i |c_i / s_i| – the scaled infinity norm. The infinity norm rather than the 2-norm so the feasibility tolerance is a statement about the worst constraint and does not loosen as segments are added.

property defect_rms

Root-mean-square scaled defect – the aggregate companion to defect_norm.

is_finite()

Whether this linearisation can be stepped from at all.

property scaled_jacobian

The Jacobian of scaled_residual.

property scaled_residual

c_i / s_i – the dimensionless defect the outer loop and its feasibility test are defined on.

worst(count=5)

The count largest scaled defects as (name, value), worst first – the “report scaled continuity defects” the issue asks a converged run to print.

class pybnf.transcription.equality.EqualitySystem

What a transcription implements to declare its equality constraints.

Two obligations: say how many constraints there are and what they are called (once, at construction – the layout and the constraint list are both static within a homotopy stage), and linearise them at a point. Everything else – the penalty, the multipliers, the inner optimizer, the certification – is somebody else’s.

equality_at() is the only place a consumer touches a simulator, which is exactly why this layer is testable without one: an offline system implements it in closed form. The method is named to match equality_at(), so a class that implements this ABC is the constraint half of a transcription problem and the two compose by inheritance rather than by an adapter.

abstract property constraint_names

One label per constraint, in residual order.

empty_model()

The zero-constraint model – what a system with nothing to enforce returns.

The final stage of a segment-count homotopy has exactly this shape: coarsened to one segment there are no knots, hence no continuity constraints, and the augmented problem is the ordinary single-shoot problem. That stage is not a special case to be branched around; it is this model, and the outer loop reduces to one inner solve on it.

abstractmethod equality_at(u)

The EqualityModel at augmented point u.

abstract property layout

The AugmentedLayout these constraints are written against.

class pybnf.transcription.equality.JacobianBlock(rows, cols, values)

One dense, non-zero region of a constraint Jacobian.

Parameters:

• rows – Contiguous constraint rows this region covers (a slice).

• cols – Contiguous augmented-variable columns it covers (a slice).

• values – The (len(rows), len(cols)) derivatives.

For multiple shooting a segment contributes three blocks per knot: dPhi_j/dtheta over the reported columns (from the PARAM sensitivity route), dPhi_j/dz_j over that segment’s own state block (the IC route – which the #563 prototype found is an IC contribution with chain-rule factor 1, so it needs no new residual math), and the constant -I over the next knot’s state block.

The augmented Lagrangian at a point

The augmented Lagrangian at a point: the smooth problem an inner optimizer walks (#563).

The outer loop (pybnf.transcription.outer) fixes the multipliers lambda and the penalty rho; what is left is an ordinary bound-constrained smooth minimisation,

L_A(u) = f(u) + lambda^T c(u) + rho/2 ||c(u)||^2

over the augmented vector, with c the scaled equality defects (pybnf.transcription.equality). This module builds that problem’s local model, in every form PyBNF’s existing optimizers consume, so that “optimizer-agnostic” is a property of the interface rather than an aspiration.

The three forms, and why the least-squares one is exact

• scalar – value and gradient, for a quasi-Newton method (lbfgs);

• least-squares – a stacked residual and Jacobian, for a trust-region least-squares method (trf);

• Gauss-Newton – gradient and a PSD hessian, for the EFIM trust-region method (gntr).

The least-squares form is the one worth spelling out. Completing the square on the two constraint terms,

lambda^T c + rho/2 ||c||^2  ==  rho/2 ||c + lambda/rho||^2  -  ||lambda||^2/(2 rho)

so with an objective that carries an exact least-squares residual (0.5||r_f||^2 == f, which is what least_squares_exact certifies) the whole augmented Lagrangian is a sum of squares up to a constant:

r_aug = [ r_f ; sqrt(rho) (c + lambda/rho) ], J_aug = [ J_f ; sqrt(rho) J_c ], 0.5 ||r_aug||^2 == L_A + ||lambda||^2/(2 rho).

The offset does not depend on u, so it changes no step, no gradient, and no accept test – but it is reported (AugmentedModel.residual_offset) rather than left for a caller to rediscover when 0.5||r||^2 does not equal the value it was told. The shifted form is also the numerically better one: it keeps the multiplier inside the square instead of adding a large linear term to a large quadratic one.

The Gauss-Newton curvature is H_f + rho J_c^T J_c. It drops the exact augmented Lagrangian’s sum_i (lambda_i + rho c_i) grad^2 c_i term, which needs constraint second derivatives nobody assembles – the same Gauss-Newton omission trf and gntr already make on the data term, and the standard one for augmented-Lagrangian SQP.

The invariant the noise scale depends on

The constraint terms never enter the likelihood. f is the fit’s own objective and nothing this module does is folded back into it: AugmentedModel.objective_value is f alone, and it is the only quantity a certification or a reported score is allowed to read. This is not tidiness. 13 of the 23 slugs in the motivating benchmark corpus estimate at least one noise scale, and an estimated sigma is fitted to the residuals it is given – so a sigma that could see continuity defects would absorb constraint violation as measurement noise, and the reported objective would stop being comparable to a single-shoot one. Keeping the penalty strictly outside f is what makes the certification step meaningful; with noise_profiling = 1 (ADR-0108) the profiled scale is defined by the data residuals alone and the separation is structural.

class pybnf.transcription.augmented.AugmentedModel(objective, equality, multipliers)

The augmented Lagrangian and its derivatives at one point, in all three forms.

Built by AugmentedSubproblem.at(); an inner optimizer asks for one per point it visits and reads whichever form it steps from.

property defect

The scaled equality defects c(u).

property gradient

grad f + J_c^T (lambda + rho c).

hessian()

The Gauss-Newton curvature H_f + rho J_c^T J_c, or None when the objective supplies neither a Hessian nor a residual to build one from.

property least_squares_exact

Whether the stacked residual models the whole augmented Lagrangian (up to residual_offset). Inherited from the objective: the constraint rows are exact squares by construction, so the only question is whether f is.

property objective_value

f(u) alone – the fit’s own objective, with no constraint term in it. The only value a certification or a reported score may read.

residual_model()

The stacked (r_aug, J_aug) a least-squares inner optimizer consumes, or None when the objective carries no residual of its own.

0.5||r_aug||^2 == value + residual_offset when least_squares_exact; when it is False the constraint rows are still exact and the whole model is not, which is the same signal GradientResult already carries – a caller that needs an exact model must step from the scalar or Gauss-Newton form instead.

property residual_offset

||lambda||^2/(2 rho), the constant by which 0.5||r_aug||^2 exceeds value. Zero on the first outer iteration and whenever there are no constraints.

property value

L_A = f + lambda^T c + rho/2 ||c||^2.

class pybnf.transcription.augmented.AugmentedSubproblem(problem, multipliers)

The augmented Lagrangian at fixed multipliers: a plain bound-constrained smooth minimisation, and the object handed to an inner optimizer.

This is the whole of the optimizer-agnostic contract. An inner solver is any callable

solve(subproblem, u0, tolerance) -> InnerOutcome

that reads lower / upper / size, calls at() at the points it visits, and returns where it stopped. It may step from the scalar form, the stacked least-squares form, or the Gauss-Newton form; the subproblem does not know or care which, and never calls back into the outer loop.

at(u)

The AugmentedModel at u. One call per visited point: the objective and the constraints are linearised together, which for a simulator-backed consumer is one pass of segment simulations rather than two.

value_and_gradient(u)

Convenience for a scalar inner optimizer (lbfgs-shaped).

class pybnf.transcription.augmented.Multipliers(values, penalty)

The outer loop’s state: one multiplier per (scaled) constraint, and the penalty.

Immutable – an update produces a new instance – so an outer iterate can record the multipliers it was solved under without defensive copying.

updated(defect, penalty=None, clamp=None)

The first-order multiplier update lambda <- lambda + rho c.

clamp bounds the result componentwise; a multiplier that runs away is the classic symptom of an infeasible or badly scaled constraint, and letting it do so turns the inner problem into an unrecoverable one rather than a slow one.

with_penalty(penalty)

The same multipliers under a raised penalty.

classmethod zeros(n_constraints, penalty)

The start of an outer loop: no multiplier information yet, so the first inner solve is a plain quadratic-penalty solve.

class pybnf.transcription.augmented.ObjectiveModel(value, gradient, residual=None, jacobian=None, hessian=None, least_squares_exact=False)

The fit’s own objective, linearised at one augmented point.

This is deliberately the same shape as GradientResult, because for #563’s first consumer it is one: the #563 prototype established that a segment-start state enters the data fit as an IC route with chain-rule factor 1, so assemble_gradient_and_fisher_hessian() builds a segment’s gradient column and Fisher block for an auxiliary variable with no new residual math. Use from_gradient_result() for that path.

Parameters:

• value – f(u) – the fit’s objective. The certified quantity; see the module docstring on why the penalty never enters it.

• gradient – df/du_aug, length layout.size.

• residual – The objective’s own least-squares residual, or None if it has none.

• jacobian – The matching (n_obs, layout.size) residual Jacobian.

• hessian – A PSD curvature model (the EFIM), or None.

• least_squares_exact – Whether 0.5||residual||^2 == value. Only then does the stacked least-squares form model the whole augmented Lagrangian.

classmethod from_gradient_result(value, grad, layout=None)

Adapt an assembled GradientResult.

Duck-typed, so this module stays free of the gradient package (and of everything it pulls in) – which is what makes the whole transcription layer importable and testable with no simulation backend present. value comes from the objective’s own evaluate; a GradientResult carries derivatives, not a score.

If layout is given and the assembled columns are the reported parameters only, the gradient and Jacobian are zero-padded into augmented space. A gradient already assembled over the augmented free-parameter list is used as-is.

class pybnf.transcription.augmented.TranscriptionProblem

What a constrained transcription implements for the outer loop.

One homotopy stage is one of these. It owns the augmented layout, can score and differentiate the fit’s objective at an augmented point, can linearise its equality constraints there, and – the step that keeps the whole method honest – can certify a reported parameter vector by reconstructing it through the fit’s ordinary unsegmented path.

Everything the layer does is expressed through these four methods, so an offline implementation with closed-form dynamics exercises the same code path a simulator-backed one does.

augmented_at(u, multipliers)

The AugmentedModel at u – one objective linearisation, one constraint linearisation, combined. The single entry point an inner optimizer’s evaluation goes through.

certify(reported)

Reconstruct reported through the fit’s ordinary single-shoot path and return a Certificate.

Returning None (the default) means this problem cannot certify, and the outer loop marks its whole result uncertified. That is a legitimate state for the last stage of a segment homotopy – coarsened to one segment the transcription already is the single-shoot problem, so its objective is its own certificate – but for any stage with constraints an uncertified score is not a fit result: the objective at an infeasible augmented point is computed on trajectories that do not join up.

abstractmethod equality_at(u)

The EqualityModel at u.

abstract property layout

This stage’s AugmentedLayout.

property name

A short label for the stage trace (multiple shooting uses 'm=4').

abstractmethod objective_at(u)

The ObjectiveModel at augmented point u.

The outer loop

The optimizer-agnostic augmented-Lagrangian outer loop (#563).

The outer loop owns the multipliers and the penalty; an inner solver – any callable matching the contract below – owns the search. Separating them is the point: the transcription layer has to work with gntr (the #563 MVP’s inner optimizer), with trf and lbfgs, and offline with neither, and none of those may need to know that a multiplier exists.

The inner-solver contract

outcome = inner_solver(subproblem, u0, tolerance)

subproblem is an AugmentedSubproblem – fixed multipliers, a box, and at(u) giving the scalar / least-squares / Gauss-Newton forms of the augmented Lagrangian. tolerance is the outer loop’s current inner-optimality target omega_k, which starts loose and tightens; a solver that ignores it is correct but slower. The return is an InnerOutcome: where it stopped, whether it converged, how much it spent. The solver never calls back into the outer loop, and the outer loop never inspects the solver.

The schedule, and why it starts tight

The update is the classical Hestenes-Powell first-order rule inside the Conn-Gould-Toint / LANCELOT test-and-tighten frame (Nocedal & Wright, Algorithm 17.4): solve the subproblem to omega_k; if the scaled defect met its target eta_k, accept the step and update lambda <- lambda + rho c, then tighten both targets; otherwise keep lambda, raise rho, and reset the targets from the new rho.

The defaults depart from the obvious reading of the multiple-shooting literature, on measurement. Balsa-Canto et al. argue that the method’s benefit comes from allowing discontinuity, which suggests starting the penalty loose. On the motivating problem the #563 prototype measured the opposite (issue #563, finding 5.1): from one start, rho_0 = 0.1, gamma = 3 reached -178.38 in 124 s while rho_0 = 10, gamma = 5 reached -200.70 in 62 s – better and at half the cost. Too loose is not merely ineffective; it is expensive, because the inner solve on a nearly-unconstrained subproblem never converges and burns its whole budget every outer iteration. So PenaltySchedule starts at rho_0 = 10 and grows by 5.

Stalling

A run has two ways of going nowhere, and the guard against them has to span both branches of the schedule.

On the accepting branch: the convergence test needs the defect and the first-order optimality below tolerance in one iterate, and once a run is feasible the schedule’s inner tolerance is already floored – so an inner solver that cannot drive the augmented Lagrangian’s optimality lower re-solves a near-identical subproblem every remaining outer iteration.

On the penalty-raising branch the failure is worse, because raising rho is only justified if the previous inner solve did something. An inner solver that fails on an ill-conditioned subproblem and returns its own start point leaves the defect exactly where it was – which reads as “not feasible enough”, raises rho by gamma, and hands the same solver a strictly harder problem. Measured on the offline shooting problem, that death spiral runs the penalty from 1.25e3 to the 1e8 ceiling over ~15 outer iterations during which the point never moves at all and the augmented gradient grows to 2e6. So a penalty raise does not reset the stall counter.

Progress is the scaled defect improving or the point moving – deliberately not the optimality improving, which is not comparable across a change of rho (a raise scales the augmented gradient by gamma).

Stalling is a state to report, not a failure: a feasible, stalled run has found whatever it found and could not certify a KKT residual for it, which is different from failing and different again from having converged.

The defect report

The aggregate scaled defect norm says how far an iterate is from feasible; it does not say where. Every iterate therefore also carries the largest individual scaled defects by name (WORST_DEFECTS of them) and their RMS, which for multiple shooting reads experiment1@1/2::Z_state – the knot that did not close and the state it did not close in. That is the “report scaled continuity defects” (plural) the issue asks for, and it is stated here rather than in the consumer because the scaling is what makes the numbers comparable across states of different magnitude, and the scaling lives in this layer.

Best-iterate certification

The loop certifies every outer iterate, not just the last, and reports the best (issue #563, finding 5.3: on one prototype start the final stage held -147.0 while an earlier iterate certified at -196.3). Certification is the transcription’s honesty mechanism – reconstruct the reported parameters through the fit’s ordinary unsegmented path and score them there – so it is also the only ranking key that is comparable with an ordinary fit’s. The augmented objective at an infeasible point is computed on trajectories that do not join up and is not a fit result; OuterResult says so explicitly when the problem could not certify (OuterResult.certified).

class pybnf.transcription.outer.AugmentedLagrangian(problem, inner_solver, schedule=None, max_outer=25, shared_best=None, stop_check=None, on_iterate=None, max_stall=3)

The outer loop.

Parameters:

• problem – The TranscriptionProblem – one homotopy stage.

• inner_solver – The inner-solver callable (see the module docstring).

• schedule – The PenaltySchedule; the defaults carry finding 5.1.

• max_outer – Cap on outer iterations.

• max_stall – Consecutive outer iterations that neither improve the scaled defect nor move the point before the run stops with stop_reason = 'stalled'. See Stalling in the module docstring for what this guards and why it spans both branches of the schedule.

• shared_best – An optional CertifiedBest that every iterate is also offered to, so a homotopy tracks one best across all its stages. Each run() keeps its own local best regardless, which is what OuterResult.best reports.

• stop_check – Optional zero-argument callable returning True when the run must stop – the seam a wall-clock budget (wall_time_fit, ADR-0093) plugs into without this module importing it.

• on_iterate – Optional callback given each CertifiedIterate as it is produced, for progress logging.

STALL_FACTOR = 0.9

Relative improvement in the scaled defect that counts as progress for the stall detector. Deliberately generous: the point is to catch a loop that has stopped moving, not to police the rate at which it moves.

STALL_STEP = 1e-10

Relative step below which an outer iterate counts as not having moved.

run(u0, multipliers=None)

Solve from augmented start point u0; return an OuterResult.

class pybnf.transcription.outer.Certificate(objective, accepted=True, certified=True, detail='')

The verdict of reconstructing a parameter vector through the fit’s ordinary path.

Parameters:

• objective – The score the reconstruction produced – comparable with any ordinary PyBNF fit’s, which is the whole point.

• accepted – Whether that score is usable. A reconstruction that fails to simulate, or whose objective is not finite, is rejected: it is not a fit result, and ranking it would let a transcription report a number no single-shoot run can reproduce.

• certified – Whether the score came from the ordinary unsegmented path. False marks a score taken from the augmented problem itself – legitimate only where the transcription is the ordinary problem (the one-segment stage of a homotopy).

• detail – Free text for the log (a rejection reason, a defect norm, a simulation error).

classmethod accept(objective, detail='')

A reconstruction that reproduced a finite objective.

classmethod reject(detail)

A reconstruction that did not.

classmethod uncertified(objective, detail='')

The augmented problem’s own objective, standing in for a reconstruction that was not performed.

class pybnf.transcription.outer.CertifiedBest

The best certified iterate seen so far – across outer iterations, and across homotopy stages (pybnf.transcription.homotopy shares one of these over the whole ladder).

Only Certificate.accepted records compete. Ties keep the earlier one, so a later iterate has to be strictly better to displace an established result.

offer(record)

Consider record; return True if it became the best.

class pybnf.transcription.outer.CertifiedIterate(stage, iteration, reported, point, certificate, defect_norm, objective_value, augmented_value, penalty, optimality, defect_rms=0.0, worst_defects=(), n_constraints=0)

One outer iterate, with the certificate earned by its reported parameters.

defect_report(count=None)

The per-constraint defect breakdown, as one line, or '' when there are no constraints (the unsegmented rung of a ladder).

defect_rms

Root-mean-square scaled defect, the aggregate companion to the infinity norm.

n_constraints

How many constraints there were in total, so a report that lists only the worst few can say so (“the 8 largest of 1068”) rather than reading as the whole set.

property score

The ranking key: the certified objective.

worst_defects

The largest individual scaled defects as (name, value), worst first – worst(). For multiple shooting a name is '<experiment>@<fraction>::<state>', so this says which knot did not close and in which state rather than only how far off the worst one was.

class pybnf.transcription.outer.InnerOutcome(point, converged=True, n_evaluations=0, message='')

Where an inner solver stopped.

Parameters:

• point – The augmented vector it stopped at.

• converged – Whether it met the tolerance it was given (as opposed to running out of iterations or budget). The outer loop will not declare convergence off an inner solve that did not converge, but it will keep going – an approximate inner minimisation is what the method is designed around.

• n_evaluations – Model evaluations spent, for the cost accounting the #563 acceptance benchmark reports.

• message – Free text for the run log.

class pybnf.transcription.outer.OuterResult(stage, iterates, best, final_point, multipliers, converged, stop_reason, certified, n_inner_evaluations, n_outer_evaluations, defect_norm, optimality, defect_rms=0.0, worst_defects=())

What one augmented-Lagrangian run produced.

defect_report(count=None)

The per-constraint defect breakdown at final_point, as one line.

The aggregate defect_norm says how far from continuous the transcription ended; this says where. Empty when there are no constraints, which is the unsegmented rung of a ladder rather than a degenerate case.

defect_rms

RMS scaled defect at final_point.

optimality

Projected-gradient first-order optimality at final_point, measured under the last multipliers – the other half of the KKT test the defect norm starts.

summary()

One line for the run log.

worst_defects

The largest individual scaled defects at final_point, worst first.

class pybnf.transcription.outer.PenaltySchedule(initial_penalty=10.0, growth=5.0, max_penalty=100000000.0, optimality_tol=1e-06, feasibility_tol=1e-06, multiplier_clamp=10000000000.0)

The penalty / tolerance schedule of the outer loop.

Parameters:

• initial_penalty – rho_0. Default 10.0 – tight, per the finding in the module docstring, not the loose start the literature’s motivation suggests.

• growth – gamma, the factor rho is multiplied by when an outer iteration fails its feasibility target. Default 5.0.

• max_penalty – The ceiling. Reaching it without meeting the feasibility target stops the loop (stop_reason = 'penalty_ceiling') rather than grinding on an ill-conditioned subproblem: a penalty that large means the constraints are infeasible at this transcription, which is information, not a reason to keep going.

• optimality_tol – The projected-gradient first-order tolerance the run must reach to declare convergence. Default 1e-6, deliberately looser than trf / gntr’s 1e-8: those measure the fit’s gradient, while this measures the augmented Lagrangian’s, whose penalty term carries a factor of rho. Grinding that to 1e-8 needs penalty raises that make the subproblem worse conditioned than the answer requires – and the answer is certified by reconstruction, not by a KKT residual, so the extra digits buy nothing a certificate does not already establish.

• feasibility_tol – The scaled-defect target for convergence. Dimensionless, because the constraints are scaled (pybnf.transcription.equality). It also floors the schedule’s own feasibility target: a point already this feasible is never a reason to raise the penalty.

• multiplier_clamp – Componentwise bound on lambda. A multiplier that runs away is the signature of an infeasible or badly scaled constraint; clamping keeps the subproblem solvable so the loop can reach its ceiling and report that, rather than failing inside an inner solve.

eta (the feasibility target a step must meet to earn a multiplier update) and omega (the inner-optimality target) follow Algorithm 17.4 with mu = 1/rho: omega = mu and eta = mu**0.1 on a reset, eta *= mu**0.9 and omega *= mu on a success.

TARGET_FLOOR = 1e-14

Floors under the two targets, so a large penalty cannot drive either to a value no inner solver can meet and stall the loop at “not converged” forever.

raised(penalty)

The next penalty, capped at max_penalty.

reset_targets(penalty)

(eta, omega) for a freshly raised penalty.

tighten(penalty, eta, omega)

(eta, omega) after an outer iteration that met its feasibility target.

pybnf.transcription.outer.WORST_DEFECTS = 8

How many individual scaled defects an iterate carries for the report. A defect list is one entry per (knot, state), so a big model at a fine rung has thousands of them; the worst few are what answers “which knot did not close, and in which state”, and keeping every iterate’s whole vector would grow a run’s memory with the model’s state count for numbers nobody reads. The aggregate defect_norm and defect_rms are kept in full alongside.

pybnf.transcription.outer.projected_gradient_norm(point, gradient, lower, upper)

First-order optimality of a bound-constrained problem at point: ||P_[l,u](x - g) - x||_inf.

The standard projected-gradient stationarity measure – zero exactly at a KKT point of the box-constrained subproblem, and (unlike a raw gradient norm) correctly reading as zero for a coordinate pinned at a bound by a gradient pushing outward, which is what trf / gntr’s Coleman-Li optimality test measures too.

The outer loop measures this rather than trusting the inner solver’s own converged flag. An inner solver that stopped on its iteration cap at a point that happens to be stationary should end the run; one that reports success against a loose internal tolerance should not.

The transcription homotopy

The transcription homotopy: solve a ladder of transcriptions, not one (#563).

A constrained transcription has a knob – how many segments, how many collocation nodes, how much of the state is made auxiliary – and the whole method’s difficulty runs along it. The obvious MVP fixes the knob and solves once. The #563 prototype measured that that is the wrong shape, twice over, and this module is the consequence.

The ladder is the mechanism, not a refinement. On the motivating problem the stage trace is the result (issue #563, finding 5.2):

m=8: -132.32   m=4: -166.52   m=2: -221.91   m=1: -248.07

Every segmented stage scores worse than a flat line (-164.68); the coarsening is what converts them, and the run that produced this trace is the first solve of a problem fifteen independent global searches did not solve. Fixing the segment count and solving once reaches m=4 and stops.

And it starts in the middle, not at the far end. The original formulation proposed starting with many short segments – the easiest landscape – and coarsening toward one. Measured, m=8 routinely certifies worse than its own start: with one observed state of three and ~14 points per segment, the segmented problem is under-determined, so the data term is satisfiable without correct dynamics and the auxiliary states absorb the rest. Over eight paired starts, 4-2-1 had the best tail at moderate cost and 8-4-2-1 did not. coarsening_ladder() therefore defaults to starting at 4.

Two rules the driver follows, and the asymmetry between them

Continue from the last point; report the best certified one. A stage seeds the next from where it finished – that is what continuation means, and the trace above is what it buys. But the answer the run reports is the best certified iterate across every stage (CertifiedBest), because a stage’s final point is not reliably its best one: on one prototype start the final stage held -147.0 while an earlier iterate certified at -196.3 (finding 5.3). The two rules disagree on purpose.

Multipliers are not carried between stages. Coarsening changes the constraint set – different knots, different count, different meaning – so a multiplier estimated for a constraint that no longer exists is not an estimate of anything. Each stage restarts from lambda = 0 at rho_0. The auxiliary variables do carry over, by name (carry_over()), which is where the continuation information actually lives.

class pybnf.transcription.homotopy.HomotopyResult(stages, best, stop_reason, certified)

What a whole ladder produced.

property reported

The reported free parameters of the best certified iterate – the fit result. None if nothing certified.

trace()

The ladder in one line – 'm=4: -166.52   m=2: -221.91   m=1: -248.07'.

The single most informative artifact a homotopy produces: it shows whether the coarsening is converting the segmented stages, which is the mechanism the method rests on, and it shows it before the run finishes.

class pybnf.transcription.homotopy.StageResult(name, outer, final_point, layout)

One rung of the ladder.

property best_score

The best objective certified within this stage – the number the trace prints.

pybnf.transcription.homotopy.coarsening_ladder(start=4, factor=2, stop=1)

The default homotopy ladder: (4, 2, 1).

Parameters:

• start – The finest transcription’s knob value. Defaults to 4 rather than the largest affordable number – see the module docstring on why the fine end is the wrong end.

• factor – The coarsening factor between stages.

• stop – The coarsest stage, reached exactly. 1 is the ordinary unsegmented problem, so a ladder always ends by solving the fit the user actually asked about.

Values are integers, strictly decreasing, ending at stop.

pybnf.transcription.homotopy.run_homotopy(stages, inner_solver, reported_start, schedule=None, max_outer=25, stop_check=None, on_iterate=None, on_stage=None)

Run an augmented-Lagrangian solve on each stage in turn, warm-starting down the ladder.

Parameters:

• stages – The transcriptions, finest first. Each item is either a TranscriptionProblem or a callable (reported) -> TranscriptionProblem – the callable form exists because a stage’s auxiliary variables usually have to be seeded from the incoming parameters (for multiple shooting, by reading a nominal trajectory at the current theta), which is not knowable when the ladder is written.

• inner_solver – The inner-solver callable, passed through to AugmentedLagrangian.

• reported_start – The fit’s start point, in sampling space, over the reported free parameters.

• schedule – The shared PenaltySchedule.

• max_outer – Outer-iteration cap per stage.

• stop_check – Zero-argument callable; True ends the ladder between or within stages (the wall-clock-budget seam).

• on_iterate – Per-iterate callback, passed through.

• on_stage – Optional callback given each StageResult as it completes.

Returns a HomotopyResult whose best is the best certified iterate over the whole ladder.