Paper deep dive

ConicFlow: constrained Inference-time optimal control for flow matching models

Shervin Khalafi, Kasra Fallah, Sung Cho, Dongsheng Ding, James Anderson, Alejandro Ribeiro

Year: 2026Venue: arXiv preprintArea: cs.AIType: PreprintEmbeddings: 162

Intelligence

Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 91%

Last extracted: 10/14/2026, 9:46:35 AM

Summary

The paper introduces ConicFlow, an inference-time method for constrained sampling in flow matching models. It formulates the problem as an optimal control task balancing reward maximization against deviation from pretrained dynamics, subject to hard terminal constraints. To ensure computational tractability, it employs receding-horizon control with a novel terminal-estimate variant that uses pretrained flow maps to predict terminal states, avoiding local truncation errors. The method uses a receding-horizon primal-dual algorithm with theoretical guarantees on suboptimality and convergence.

Entities (6)

ConicFlow · method · 100%Flow Matching Models · model-class · 95%Primal-Dual Algorithm · algorithm · 90%Receding-Horizon Control · algorithm-technique · 90%Terminal-Estimate Variant · method-variant · 90%Flow Maps · mathematical-concept · 85%

Relation Signals (5)

ConicFlow → solves → Flow Matching Models

confidence 95% · We propose ConicFlow, which frames constrained sampling from a flow matching model at inference time

ConicFlow → uses → Receding-Horizon Control

confidence 95% · We discretize this problem and adapt receding-horizon control, optimizing only a short window of the trajectory at each step

ConicFlow → employs → Primal-Dual Algorithm

confidence 90% · address the resulting sequence of receding-horizon subproblems with a receding-horizon primal-dual algorithm

Terminal-Estimate Variant → ispartof → ConicFlow

confidence 90% · We further propose a terminal-estimate variant... Experiments... show that ConicFlow achieves high reward

Terminal-Estimate Variant → utilizes → Flow Maps

confidence 90% · proposes a terminal-estimate variant that avoids accumulating local truncation errors by evaluating the reward and constraints at the predicted terminal point... uses the pretrained flow map

Cypher Suggestions (0)

No Cypher suggestions yet.

Abstract

Abstract:We propose ConicFlow, which frames constrained sampling from a flow matching model at inference time as an optimal control problem that balances deviation from the pretrained trajectory against reward, subject to hard terminal constraints. We discretize this problem and adapt receding-horizon control, optimizing only a short window of the trajectory at each step for computational tractability. We address the resulting sequence of receding-horizon subproblems with a receding-horizon primal-dual algorithm. Its dual multipliers capture the difficulty of satisfying each constraint and are refined throughout sampling, steering generation toward feasible final samples. We further propose a terminal-estimate variant that avoids accumulating local truncation errors by evaluating the reward and constraints at the predicted terminal point rather than at the end of the current horizon. Under standard regularity assumptions, we establish finite-horizon guarantees on suboptimality. Under convex-affine structure, we derive a sharper bound that scales linearly with omitted horizon length, and we show that the terminal-estimate variant is provably tighter when the flow dynamics change slowly. Experiments in both low- and high-dimensional settings with a recent large-scale flow matching model show that ConicFlow achieves high reward and strong constraint satisfaction while deviating less from the pretrained flow than recent baselines.

PDF

Open source PDF →Open local PDF →

Trouble viewing inline? Open PDF directly →

Full Text

158,788 characters extracted from source content.

Expand or collapse full text

mnlargesymbols’164 mnlargesymbols’171 ConicFlow: CONstrained Inference-time optimal Control for FLOW matching Models Shervin Khalafi Kasra Fallah11footnotemark: 1 University of Pennsylvania Columbia University shervink@engineering.upenn.edu kasra.fallah@columbia.edu Sung Cho Dongsheng Ding University of Pennsylvania University of Tennessee, Knoxville sungcho9@engineering.upenn.edu dongshed@utk.edu James Anderson Alejandro Ribeiro Columbia University University of Pennsylvania ja3451@columbia.edu aribeiro@engineering.upenn.edu †thanks: Equal contribution. Abstract We propose ConicFlow, which frames constrained sampling from a flow matching model at inference time as an optimal control problem that balances deviation from the pretrained trajectory against reward, subject to hard terminal constraints. We discretize this problem and adapt receding-horizon control, optimizing only a short window of the trajectory at each step for computational tractability. We address the resulting sequence of receding-horizon subproblems with a receding-horizon primal–dual algorithm. Its dual multipliers capture the difficulty of satisfying each constraint and are refined throughout sampling, steering generation toward feasible final samples. We further propose a terminal-estimate variant that avoids accumulating local truncation errors by evaluating the reward and constraints at the predicted terminal point rather than at the end of the current horizon. Under standard regularity assumptions, we establish finite-horizon guarantees on suboptimality. Under convex-affine structure, we derive a sharper bound that scales linearly with omitted horizon length, and we show that the terminal-estimate variant is provably tighter when the flow dynamics change slowly. Experiments in both low- and high-dimensional settings with a recent large-scale flow matching model show that ConicFlow achieves high reward and strong constraint satisfaction while deviating less from the pretrained flow than recent baselines. 11 1 Our implementation code can be found at: https://github.com/shervinkhalafi/conicflow 1 Introduction Diffusion models and flow matching models generate high-quality samples across a variety of modalities (Ho et al., 2020; Lipman et al., 2022; Rombach et al., 2022). However, sample quality alone is often insufficient: many applications require samples that both achieve high reward and satisfy hard constraints. Rewards can encode user preferences, such as image aesthetics evaluated by models trained on human feedback (Wu et al., 2023), consistency with observations in inverse problems (Ben-Hamu et al., 2024), or safety-related criteria (Schramowski et al., 2023). Optimizing these rewards alone, however, does not guarantee satisfaction of hard constraints. The challenge is to achieve high reward and terminal feasibility while limiting deviations from the pretrained dynamics. In both theory and practice, ensuring that samples generated by these models satisfy hard constraints is underexplored. Existing approaches fall broadly into two categories. One approach is to adapt model parameters through fine-tuning (Black et al., 2024; Domingo i Enrich et al., 2025), but it can be prohibitively expensive. Furthermore, such an approach would require finetuning again for each new reward/constraint. Guidance, inference-time steering, and projection based methods for constraint satisfaction (Ben-Hamu et al., 2024; Christopher et al., 2026), avoid the computational overhead of fine-tuning, usually at the cost of additional overhead during sampling. Our approach is an inference-time approach. We provide extended discussion on related approaches in Appendix A. A pre-trained flow matching model describes the dynamics of the sample evolving across time. We need to deviate from these dynamics in order to satisfy constraints and high reward at the endpoint, meanwhile ideally we would want these deviations to be small, because we still want to end up reasonably close to the data manifold. Given this description, the problem naturally lends itself to an optimal control formulation, where the deviation from the flow dynamics can be viewed as the control input that we are optimizing over. This leads to the continuous-time constrained inference problem (P-C) where r: ℝd→ℝR^d encodes a soft preference over the generated sample, while h: ℝd→ℝmR^d ^m represents m hard terminal requirements. minu⁡(t) _u(t) JC(u):=∫01‖u⁡(t)‖22λdt−r⁡(x1u) J_C(u)\; :=\; _0^1 \|u(t)\|^22λ\,dt-r(x_1^u) (P-C) s.t. .t. x˙tu=bt(xtu)+u⁡(t),given x0, x_t^u\;=\;b_t(x_t^u)+u(t), x_0, h⁡(x1u)≤ 0. h(x_1^u)\;≤\;0. Full-horizon constrained control can be computationally demanding, particularly for high-dimensional flow models. We therefore adopt receding-horizon control (RHC), which replaces the full problem by a sequence of shorter optimizations and replans after each intervention. (García et al., 1989; Mayne et al., 2000; Mattingley et al., 2011; Bemporad et al., 2000). This improves tractability at the cost of horizon-truncation error, motivating our formulations (RH-H) and (RH-T) in Section 3. Specifically, in RHC, rather than optimizing the full trajectory at once, at each step of the dynamics (in the case of flow matching, at each sampling step) we optimize over a short window and then apply the first optimized control action before shifting the window forward by one step and repeating this procedure until we reach terminal time. Note that in (P-C), both the reward and the constraints are evaluated only at the terminal time. Classical receding-horizon control introduced in (RH-H) instead evaluates them at the end of each truncated window, which creates a mismatch with the original objective. We therefore propose (RH-T), which uses the pretrained flow map (Boffi et al., 2025) to predict the terminal state reached from the end of each window and evaluates the reward and constraints at this prediction. To the best of our knowledge, this is the first work to use flow maps within receding-horizon control, allowing each subproblem to account for the full-horizon objective without extending the control horizon. Our contributions are fourfold. (i) Constrained-control formulation and theory. We formulate reward-guided generation with hard terminal constraints as an optimal control problem that penalizes deviations from pretrained dynamics, and bound the errors introduced by discretization and horizon truncation. (i) Terminal prediction through flow maps. We propose a variant of receding-horizon control (RH-T) that uses the pretrained flow map to predict the uncontrolled continuation beyond each optimization window. This enables evaluation of rewards and constraints at the final sampling time without extending the optimized control horizon. (i) Primal-dual algorithm and convergence guarantees. We develop a receding-horizon primal-dual algorithm and, under affine dynamics, establish sublinear convergence rates. We further quantify how local optimization errors propagate and establish convergence to the exact receding-horizon trajectory under uniform bounds. (iv) Experiments. On synthetic data and image editing, ConicFlow improves the reward-control trade-off and achieves near-100%100\% empirical feasibility with up to a 2020-fold reduction in control energy relative to the baseline. Preliminaries on flow matching and flow maps. A flow matching model generates samples by evolving an initial state x0∼p0x_0 p_0 according to x˙t=bt(xt) x_t=b_t(x_t) for t∈[0,1]t∈[0,1], where b: [0,1]×ℝd→ℝd[0,1]×R^d ^d is a learned velocity field. Training aims to transport the source distribution p0p_0 to the target distribution pdatap_data at t=1t=1, and can be formulated using stochastic interpolants (Albergo et al., 2025). Assuming the ODE has unique solutions over the time interval, its flow map Xs,t:ℝd→ℝdX_s,t:R^d ^d sends a state at time s to the corresponding state at time t: Xs,t(xs)=xtX_s,t(x_s)=x_t, with Xs,s(x)=xX_s,s(x)=x for 0≤s≤t≤10≤ s≤ t≤ 1. Learned approximations of these solution operators enable direct prediction across time intervals, accelerating sampling without explicitly integrating every intermediate step (Boffi et al., 2025). 2 Problem Formulation Problem statement. We consider the problem of steering a pretrained flow model toward samples that achieve a desired terminal reward while simultaneously satisfying hard output constraints. Rather than modifying or retraining the underlying flow model, we introduce an inference-time control utu_t that perturbs the nominal vector field. This leads to the continuous-time constrained inference problem (P-C). The reward r and the constraint h play distinct roles: The constraint specifies strict properties that the generated output must satisfy, while the reward ranks feasible terminal samples. The quadratic energy term penalizes deviation from the pretrained generative dynamics. In particular, smaller λ places a larger penalty on intervention and keeps the controlled trajectory closer to the nominal flow, while larger λ permits stronger inference-time steering. Discrete-time constrained inference. In practice, flow-based inference is carried out using finitely many vector-field evaluations. We introduce a uniform forward-Euler discretization of the flow dynamics in (P-C) over t0<t1<⋯<tNt_0<t_1<·s<t_N: xk+1=xk+Δt(btk(xk)+uk)x_k+1=x_k+ t (b_t_k(x_k)+u_k ), where t0=0t_0=0, tN=1t_N=1, and Δt=1/N t=1/N. By applying Euler discretization, we approximate (P-C) as a discrete-time constrained inference problem, minu0:N−1 _u_0:N-1 JD(u):=∑k= 0N−1Δt‖uk‖22λ−r⁡(xNu) J_D(u)\; :=\; _k\,=\,0^N-1 t \|u_k\|^22λ-r(x_N^u) (P-D) s.t. .t. xk+1=xk+Δt(btk(xk)+uk),given x0, x_k+1\;=\;x_k+ t (b_t_k(x_k)+u_k ), x_0, h⁡(xNu)≤ 0. h(x_N^u)\;≤\;0. The constrained inference problem (P-D) is finite-dimensional. To capture the approximation error, we prove that the optimal value of the discrete-time problem (P-D) approximates that of the continuous-time problem (P-C) with error O⁡(1/N)O(1/N) in Lemma 1. Let uC⋆u_C be an optimal control for (P-C), and uD⋆(N)u_D (N) be an optimal control for the N-step discretized problem (P-D). Assumption 1 (Regularity, feasibility, and repairability). (i) On the relevant energy-bounded reachable set, bt(x)b_t(x) is LbL_b-Lipschitz in state and LbtL_b^t-Lipschitz in time, while r and h are LrL_r- and LhL_h-Lipschitz, respectively; (i) there exist u0,uC⋆∈L2u^0,u_C ∈ L^2 such that h⁡(x1u0)≤0h(x_1^u^0)≤ 0, JC(u0)<∞J_C(u^0)<∞, and uC⋆∈argminu:h⁡(x1u)≤0JC(u)u_C ∈ _u:\,h(x_1^u)≤ 0J_C(u); (i) there exist κ,ρ0>0κ, _0>0, independent of N, such that any relevant continuous- or discrete-time control with terminal violation at most ρ0 _0 admits a feasible perturbation whose control norm is at most κ times the violation. We state Assumption 1 in a compact form with its precise statement deferred to Appendix B.1. Lemma 1 (Discretization error). Let Assumption 1 hold. Then, there exist constants N1∈ℕN_1 , and Cdisc>0C_disc>0 that is independent of N, such that for all N≥N1N≥ N_1, |JD(uD⋆(N))−JC(uC⋆)|≤CdiscN. |J_D (u_D (N) )-J_C(u_C ) |≤ C_discN. (1) Lemma 1’s proof is provided in Appendix B.1. In Lemma 1, we show that the objective gap between the discrete- and continuous-time problems, (P-D) and (P-C), decays to zero at a rate O⁡(1/N)O(1/N). Path-space constrained inference. Although (P-D) is finite-dimensional, solving it directly requires a full-horizon constrained nonlinear optimization over all N control inputs, with repeated evaluations of the pretrained flow map. For high-dimensional flow maps, this can be computationally prohibitive at inference time. We therefore rewrite Problem (P-D) in terms of the generated trajectory rather than the control sequence. The control input uku_k in (P-D) is directly expressed by uk=(xk+1−xk)/Δt−btk(xk)u_k=(x_k+1-x_k)/ t-b_t_k(x_k). Thus, we can remove the explicit flow dynamics in (P-D) and introduce a discrete-time path-space problem as follows. minx1:N _x_1:N JD(u):=∑k= 0N−112λΔt‖xk+1−xk−Δtbtk(xk)‖2−r⁡(xN) J_D(u)\; :=\; _k\,=\,0^N-1 12λ t \|x_k+1-x_k- t\,b_t_k(x_k) \|^2-r(x_N) (P-DP) s.t. .t. h⁡(xN)≤0,given x0. h(x_N)≤ 0, x_0. Like Problem (P-D), (P-DP) remains finite-dimensional, but it does not include constraints from the flow dynamics. Each stage cost of (P-DP) measures the discrepancy between the optimized displacement xk+1−xkx_k+1-x_k and the pretrained flow drift Δtbtk(xk) t\,b_t_k(x_k). Thus, path-space optimization can be interpreted as deforming the nominal generative trajectory while penalizing the magnitude of that deformation. However, solving (P-DP) still requires optimizing the entire N-step generation trajectory. For flow matching models where the vector field is typically given by a deep neural network, repeatedly optimizing such a full path requires a large number of evaluations of the high-dimensional vector field. This becomes more costly as the number of sampling steps grows. To address this challenge, we will establish a receding-horizon primal-dual approach that allows us to retain the reward and constraint structure of the full inference problem (P-C) while optimizing a short portion of the future trajectory at each time step during inference. 3 Receding-horizon primal-dual To make (P-DP) more tractable, we first propose to replace the whole trajectory with a local subset of it. At generation step j, starting from current state xj=x¯jx_j= x_j, we select next HjH_j states and optimize over them: (xj+1,…,xj+Hj)(x_j+1,…,x_j+H_j), where Hj≪NH_j N. We apply the first induced control, and then repeat the optimization from the resulting state. This produces a receding-horizon procedure analogous to model predictive control (García et al., 1989; Mayne et al., 2000; Alonso et al., 2022), but expressed directly in the path space of the pretrained flow. Where should the terminal objective and constraint be evaluated? Denote a local path xj+1:j+Hj=(xj+1,…,xj+Hj)x_j+1:j+H_j=(x_j+1,…,x_j+H_j) and abbreviate it as jx_j. Restricting the optimization horizon introduces an important design choice. A direct truncation evaluates the reward and constraint at the end of the local optimization window: minj _x_j ∑ℓ=j+Hj−112λΔt‖xℓ+1−xℓ−Δtbtℓ(xℓ)‖2−r⁡(xj+Hj) _ \,=\,j^j+H_j-1 12λ t \|x_ +1-x_ - t\,b_t_ [rgb]0.66,0,0.01 (x_ ) \|^2-r(x_j+H_j) (RH-H) s.t. .t. h⁡(xj+Hj)≤0. h(x_j+H_j)≤ 0. Formulation (RH-H) is computationally attractive, since it involves shorter paths to optimize. Meanwhile, (RH-H) relaxes the original inference problem (P-C): the terminal reward and feasibility requirement are now imposed at time tj+Hjt_j+H_j rather than at the final generation time. Once the local cost has been optimized at j⋆x_j , the first residual intervention is recovered as uj⋆=(xj+1⋆−x¯j)/Δt−btj(x¯j)u_j =(x_j+1 - x_j)/ t-b_t_j( x_j). Only this first intervention is executed before shifting the window forward by one step and replanning (Mayne et al., 2000; Alonso et al., 2022). A distinctive feature of our setting is that the pretrained flow itself provides a natural prediction of how a locally optimized state would evolve up to the true terminal time. Rather than evaluating reward and feasibility prematurely at xj+Hjx_j+H_j, we can propagate this state under the uncontrolled flow and evaluate the original terminal objective at the resulting prediction. This brings each truncated subproblem closer to the full-horizon problem (P-DP) without extending the control horizon. To make this precise, we introduce the discrete flow map induced by the pretrained flow, Xi,j:ℝd→ℝdX_i,j:R^d ^d, defined recursively by Xi,i(x)=xX_i,i(x)=x and Xi,k+1(x)=Xi,k(x)+Δtbtk(Xi,k(x))X_i,k+1(x)=X_i,k(x)+ t\,b_t_k (X_i,k(x) ) for k≥ik≥ i, so that Xi,j(x)X_i,j(x) is the state obtained by propagating x from tit_i to tjt_j under the pretrained flow. Given the endpoint xj+Hjx_j+H_j of the local window, the terminal prediction is then x^N=Xj+Hj,N(xj+Hj) x_N=X_j+H_j,N(x_j+H_j). Evaluating the terminal objective at x^N x_N in place of xj+Hjx_j+H_j yields our terminal-rollout formulation, RH-T, described below. This gives the terminal-rollout formulation minj _x_j ∑ℓ=j+Hj−112λΔt‖xℓ+1−xℓ−Δtbtℓ(xℓ)‖2−r⁡(x^N) _ \,=\,j^j+H_j-1 12λ t \|x_ +1-x_ - t\,b_t_ (x_ ) \|^2-r ( x_N ) (RH-T) s.t. .t. h⁡(x^N)≤0. h ( x_N )≤ 0. Compared with (RH-H), (RH-T) truncates only the control horizon while retaining the terminal objective and constraint. Our analysis in Section 4.1 identifies when this yields a tighter guarantee; the comparison is detailed in Remark 1. Receding-Horizon Primal-Dual Approach. Directly solving a new hard-constrained nonlinear program at every generation step itself is expensive. We therefore introduce a primal–dual approach in which terminal constraints are represented by dual variables μ∈ℝ+mμ _+^m. To unify our analysis, we represent an approximate terminal state by Φj(x) _j(x), which is x or Xj+Hj,N(x)X_j+H_j,N(x) for (RH-H) and (RH-T), respectively. Denote ℰj(j):=∑ℓ=j+Hj−112λΔt‖xℓ+1−xℓ−Δtbtℓ(xℓ)‖2E_j(x_j) := _ \,=\,j^j+H_j-1 12λ t \|x_ +1-x_ - t\,b_t_ (x_ ) \|^2 and gj(j):=h⁡(Φj(xj+Hj))g_j(x_j) :=h ( _j(x_j+H_j) ). For each j, we have a frozen local problem, pj⋆(x¯j):=minjφj(j;x¯j):=ℰj(xj:j+Hj)−r(Φj(xj+Hj))s.t.gj(j)≤0.p_j ( x_j)\;\, :=\;\, _x_j\;\; \ _j(x_j; x_j) :=E_j(x_j:j+H_j)-r ( _j(x_j+H_j) ) \ .t. g_j(x_j)≤ 0. (RHjRH_j) The Lagrangian for ( RH j ) is given by ℒj(j,μ,x¯j)=φj(j,x¯j)+μ⊤gj(j)L_j(x_j,μ; x_j)= _j(x_j; x_j)+μ g_j(x_j). Define the dual function by dj(μ,x¯j):=minj⁡ℒj(j,μ,x¯j)d_j(μ; x_j) := _x_jL_j(x_j,μ; x_j). The associated dual problem is maxμ≥ 0⁡dj(μ,x¯j) _μ\,≥\,0d_j(μ; x_j), whose solution set is denoted by ℳj⋆(x¯j)M_j ( x_j). At each receding-horizon window, we hold the current state x¯j x_j fixed and alternate primal minimization with projected dual ascent: j(s)∈arg⁡minjℒj(j,μj(s),x¯j)x_j^(s)∈ _x_jL_j(x_j, _j^(s); x_j) and μj(s+1)=[μj(s)+ημgj(j(s))]+ _j^(s+1)= [ _j^(s)+ _μg_j(x_j^(s)) ]_+, where ημ>0 _μ>0 is the dual step size. The primal step optimizes the planned path for the current multipliers, while the dual step increases the weight of violated constraints and decreases that of strictly satisfied constraints. After K dual updates, we recompute the primal response at μj(K) _j^(K) and execute only its first transition. We then advance the window and warm-start its multiplier using μj+1(0)=μj(K) _j+1^(0)= _j^(K). This transfers the learned constraint weights between successive windows without assuming that their optimal multipliers coincide. We summarize the receding-horizon inference procedure in Algorithm 1 and establish its theoretical guarantees for solving the constrained inference problem (P-C) in Section 4. Algorithm 1 Receding-Horizon Primal-Dual Algorithm 1: Pretrained flow btb_t, reward r, constraints h, horizon H, local optimization steps K≥1K≥ 1, dual step size ημ _μ, a choice of (RH-H) and (RH-T) through the terminal map Φj _j. 2: Initialize x0x_0 and μ0(0)≥0 _0^(0)≥ 0 3: for j=0,…,N−1j=0,…,N-1 do 4: Set Hj←min⁡H,N−jH_j← \H,N-j\, x¯j←xj x_j← x_j and initial μj(0)≥0 _j^(0)≥ 0. 5: for s=0,1,…,K−1s=0,1,…,K-1 do 6: Compute the primal response: j(s)←argminjℒj(j,μj(s),x¯j)x_j^(s)\;←\; _x_jL_j(x_j, _j^(s); x_j). 7: Perform the projected dual update: μj(s+1)←[μj(s)+ημgj(j(s))]+ _j^(s+1)\;←\; [ _j^(s)+ _μg_j(x_j^(s)) ]_+. 8: end for 9: Compute the primal response at μj(K) _j^(K): j(K)←argminjℒj(j,μj(K),x¯j)x_j^(K)\;←\; _x_jL_j(x_j, _j^(K); x_j). 10: Execute the first state of j(K)x_j^(K): xj+1←xj+1(K)x_j+1← x_j+1^(K), and recover the residual control: uj←(xj+1(K)−x¯j)/Δt−btj(x¯j)u_j←(x_j+1^(K)- x_j)/ t-b_t_j( x_j). 11: If j<N−1j<N-1, warm start the next physical window with μj+1(0)=μj(K) _j+1^(0)= _j^(K). 12: end forreturn The executed path (x0,…,xN)(x_0,…,x_N) and residual controls (u0,…,uN−1)(u_0,…,u_N-1). 4 Theoretical Analysis To bridge the original continuous-time problem and the implementable Algorithm 1, we proceed in two steps. First, we bound the approximation error of (RH-H) and (RH-T) relative to the full discretized problem (P-D). We then bound the error introduced by finitely solving the local primal–dual problem ( RH j ). Combined with Lemma 1, these results connect Algorithm 1 to the original continuous formulation (P-C). 4.1 Bounding the Suboptimality We first quantify the error introduced by restricting the control horizon. Specifically, we compare the idealized constrained trajectories of (RH-H) and (RH-T) with the optimal solution ⋆ x of the full discretized problem (P-D). Both approximations optimize only HjH_j future steps, but differ in where the terminal reward and constraint are evaluated. Throughout this section, Assumption 1 remains in force. For the receding-horizon comparison, we additionally impose the following structure. Unlike the continuous-to-discrete analysis, terminal feasibility here is handled through a uniform Slater condition rather than local repairability. Assumption 2 (Affine flow and bounded amplification). The discretized flow is affine, xk+1=xk+Δt(Akxk+ck+uk)x_k+1=x_k+ t(A_kx_k+c_k+u_k). On the relevant reachable set, ‖Akx+ck‖≤M\|A_kx+c_k\|≤ M, and further let κA:=sup0≤k≤ℓ≤N‖∏s=kℓ−1(I+ΔtAs)‖<∞, _A\; :=\; _0\,≤\,k\,≤\, \,≤\,N\; \| _s=k -1(I+ tA_s) \|\;<\;∞, with M and κA _A independent of the horizon. Assumption 3 (Convexity and uniform feasibility). The reward r is concave and LrL_r-Lipschitz, and h is componentwise convex and LhL_h-Lipschitz. Either subproblem (RH-H) or (RH-T) admits a strictly feasible trajectory with uniform margin γ>0γ>0 and objective at most J¯D J_D. Let J¯D J_D be a uniform lower bound on the objectives of (RH-H) or (RH-T) and define χh:=Lh(JD¯−JD¯)/γ _h :=L_h( J_D- J_D)/γ and Gh:=Lr+χh.G_h :=L_r+ _h. Here, GhG_h represents the effective terminal sensitivity, with LrL_r capturing reward variation and χh _h the additional sensitivity required to maintain feasibility. We first consider (RH-H), which evaluates reward and feasibility at the end of the truncated horizon. Consequently, its approximation error depends on how far the nominal flow can move over the omitted portion of the trajectory, captured by the drift bound M. Lemma 2 (Constrained RH-H bound). Let Assumptions 1–3 hold, and M<∞M<∞, Gh<∞G_h<∞, and γ>0γ>0. Then, the optimal trajectory H,⋆ x^H, of (RH-H) satisfies 0≤J⁡(H,⋆)−J⁡(⋆)≤CHCN−HN, where CHC:=2max⁡GhM+12λGh2,M22λ.0\;≤\;J( x^H, )-J( x )\;≤\;C_H C N-HN, where \;C_H C :=2 \G_hM+ 12λ G_h^2, M^22λ \. Thus, (RH-H) incurs a horizon error linear in the omitted fraction (N−H)/N(N-H)/N. Its constant depends on the absolute motion of the pretrained flow through M, reflecting that reward and feasibility are evaluated before the true terminal time. (RH-T) instead propagates the local endpoint to the original terminal time N using the pretrained flow map. This removes dependence on the absolute nominal motion after the control horizon; the relevant quantity is instead how strongly the remaining flow amplifies perturbations, measured by κA _A. Lemma 3 (Constrained RH-T bound). Let Assumptions 1–3 hold, and assume κA<∞ _A<∞, Gh<∞G_h<∞, γ>0γ>0. Then, the optimal trajectory ^H,⋆ x^H, of (RH-T) satisfies 0≤J⁡(^H,⋆)−J⁡(⋆)≤CTCN−HN, where CTC:=12κA2Gh2λ.0\;≤\;J( x^H, )-J( x )\;≤\;C_T C N-HN, where \;C_T C := 12 _A^2G_h^2λ. -5.69046pt Both formulations thus achieve an O⁡((N−H)/N)O((N-H)/N) approximation to the full discretized problem, but their constants capture different properties of the pretrained dynamics: (RH-H) depends on nominal flow motion, whereas (RH-T) depends on perturbation amplification. The resulting comparison, including conditions under which (RH-T) yields the tighter guarantee, is detailed in Remark 1. Remark 1 (When is terminal rollout beneficial?). For our constrained problem, κA≤2 _A≤ 2 implies CTC≤CHCC_T C≤ C_H C. Thus, when the remaining pretrained flow does not strongly amplify intermediate-state perturbations, terminal rollout gives the tighter guarantee. The analogous unconstrained comparison is given in Appendix B.2. Intuitively, (RH-T) is most useful when the state continues to evolve after the control horizon, but this evolution remains predictable. 4.2 Convergence of Dynamic Primal–Dual In Lemmas 2–3, we assume that each subproblem of ( RH j ) is solved, exactly. We next quantify the error introduced by a finite number of primal–dual updates in Algorithm 1. For the jjth subproblem ( RH j ), we denote by BjB_j a lower block diagonal matrix with off-diagonals −(Id+ΔtAj+1)-(I_d+ tA_j+1), SjS_j a row vector with first block Id+ΔtAjI_d+ tA_j, and QjQ_j the linear part of the terminal map Φj _j (Qj=IQ_j=I for (RH-H)). Using the strong convexity and dual smoothness, we first show that the inner-loop primal-dual iterates (Lines 4–7 of Algorithm 1) converge to the optimal, sublinearly. Lemma 4 (Local primal–dual convergence). Let Assumptions 2–3 hold for the jjth frozen subproblem problem ( RH j ). Then, if 0<ημ<2/Ld,j0< _μ<2/L_d,j, then the inner-loop primal-dual iterates (Lines 4–7 of Algorithm 1) converge their optima, i.e., μj(s)→μ¯j∈ℳj⋆ _j^(s)→ μ_j _j and j(s)→j⋆x_j^(s) _j . Moreover, if 0<ημ≤1/Ld,j0< _μ≤ 1/L_d,j, then the gaps between the iterates (j(s),μj(s))(x_j^(s), _j^(s)) and the optima are bounded by pj⋆−dj(μj(s))≤Rj22ημs,‖j(s)−j⋆‖≤Rjσjημs,‖[gj(j(s))]+‖≤Lg,jRjσjημsp_j -d_j( _j^(s))\;≤\; R_j^22 _μs,\; \;\|x_j^(s)-x_j \|\;≤\; R_j _j _μs,\; \;\|[g_j(x_j^(s))]_+\|\;≤\; L_g,jR_j _j _μs where Rj:=dist⁡(μj(0),ℳj⋆)R_j :=dist( _j^(0),M_j ), σj:=σmin2(Bj)/(λΔt) _j := _ ^2(B_j)/(λ t) is the strong convexity modulus, Lg,j:=Lh‖Qj‖L_g,j :=L_h\|Q_j\| and Ld,j:=Lg,j2/σjL_d,j :=L_g,j^2/ _j are the Lipschitz constants of the constraint and dual functions. We defer the proof of Lemma 4 to Appendix B.5. In Lemma 4, we show that the primal-dual gap pj⋆−dj(μj(s))p_j -d_j( _j^(s)), and the primal gap ‖j(s)−j⋆‖\|x_j^(s)-x_j \| and constraint violation ‖[gj(j(s))]+‖\|[g_j(x_j^(s))]_+\| decay to zero, sublinearly. This result shows that the primal-dual update (Lines 4–7 of Algorithm 1) solves the subproblem ( RH j ) up to an error ⁡(1/s)O(1/ s). We next quantify how these errors propagate through the receding-horizon trajectory, and establish the global convergence of Algorithm 1 in Theorem 1; see Appendix B.6 for proofs. Denote the first state of the solution j⋆(ξ)x_j (ξ) to ( RH j ) initialized at ξ as πj(ξ) _j(ξ). We let the exact trajectory generated by πj _j be xj+1RH=πj(xjRH)x_j+1 RH= _j(x_j RH) with x0RH=x0x_0 RH=x_0. Theorem 1 (Global convergence). Let Assumptions 2–3 hold for the jjth frozen subproblem problem ( RH j ). Then, if 0<ημ≤1/L¯d0< _μ≤ 1/ L_d and each inner-loop primal-dual update of Algorithm 1 take K iterations, then the generated trajectory xjx_j for j=1,…,Nj=1,…,N satisfies ‖xj−xjRH‖≤R¯Γj(κ¯)σ¯ημK,‖[h⁡(xN)]+‖≤LhR¯σ¯ημK\|x_j-x_j RH\|≤ R\, _j( κ) σ\, _μK,\; \;\|[h(x_N)]_+\|≤ L_h R σ\, _μK where L¯d:=maxj⁡Ld,j L_d := _jL_d,j, R¯:=maxj⁡dist⁡(μj(0),ℳj⋆(xj)) R := _jdist( _j^(0),M_j (x_j)), σ¯:=minj⁡σj σ := _j _j, κj:=‖Bj⊤Sj‖λΔtσj _j := \|B_j S_j\|λ t\, _j κ¯:=maxj⁡κj κ := _j _j, and Γj(κ¯):=∑q= 0j−1κ¯q. _j( κ) := _q\,=\,0^j-1 κ^q. In Theorem 1, we prove that the trajectories generated by Algorithm 1 converge to the exact trajectory induced by the policy πj _j, in terms of the primal gap ‖xj−xjRH‖\|x_j-x_j RH\| and the constraint violation ‖[h⁡(xN)]+‖\|[h(x_N)]_+\|. This implies that a finite iterations of primal-dual optimization contribute an O⁡(1/K)O(1/ K) error, modulated by the sensitivity Γj(κ¯) _j( κ) of the exact policy. Corollary 1 (Error decomposition). Let the assumptions of Theorem 1 hold. Suppose J is LJL_J-Lipschitz on the relevant set. Then, the generated trajectory x(K)x^(K) by Algorithm 1 satisfies J⁡((K))−JC(uC⋆)≤CdiscN+ChorN−HN+LJNR¯ΓN(κ¯)σ¯ημK,J( x^(K))-J_C(u_C )\;≤\; C_discN+C_ hor N-HN+ L_J N\, R\, _N( κ) σ\, _μK, where Chor:=CHCC_ hor :=C_H C for (RH-H) and Chor:=CTCC_ hor :=C_T C for (RH-T). The proof of Corollary 1 is straightforward by adding and subtracting JD(uD⋆(N))J_D (u_D (N) ) and J⁡(RH)J( x RH), and applying Lemmas 1–3 and Theorem 1. The bound separates the discretization error O⁡(1/N)O(1/N) and horizon bias O⁡((N−H)/N)O((N-H)/N) from the O⁡(1/K)O(1/ K) finite-solver error. Taking H=NH=N eliminates the horizon bias O⁡((N−H)/N)O((N-H)/N), and letting K→∞K→∞ for fixed N eliminates the finite-solver error, leaving only the discretization error bounded by O⁡(1/N)O(1/N), which vanishes as N→∞N→∞, yielding global convergence to the continuous-time optimum. 5 Experiments We evaluate ConicFlow in two regimes: unconstrained reward-guided generation, where the objective trades terminal reward against control energy, and constrained generation, where the sample must additionally satisfy hard terminal constraints. We compare against the state-of-the-art approach in each regime: (i) FMRG (Huang et al., 2026) and (i) HardFlow (Li et al., 2026). FMRG is a training-free, single-trajectory method that uses soft reward guidance, and it is theoretically near-optimal only when small deviations from the pretrained flow suffice. HardFlow enforces constraints at terminal time, but with a different methodology compared to our approach, and it has been shown to outperform projection-based alternatives. See Appendix A for more discussion on related approaches. Unconstrained low-dimensional setting. We first use a flow map model trained on 2D samples from the synthetic two moons dataset and consider two rewards, both taking values in [0,1][0,1]. The first reward increases closer to the upper moon. As shown in Figure 1 (left), all methods saturate the reward at small λ (i.e., large control-energy penalties). On average, small deviations from the pretrained flow suffice to land trajectories on the upper moon. Both ConicFlow variants: (RH-H) and (RH-T) achieve the best trade-off between terminal reward and control energy, closely matched by FMRG, which is theoretically near-optimal in the small-λ regime (Huang et al., 2026). Figure 1: Trade-off between terminal reward and total control cost, traced by varying the control-cost weight λ (upper left is better). λ value displayed on each point. Left: easy reward (upper moon). Middle: difficult reward (off-manifold). Right: uncontrolled sample distribution and both rewards. To study settings that need large deviations (thus large λ), the second reward peaks off the pretrained data manifold (middle in Figure 1). FMRG fails outside the small-λ regime as its control energy increases without noticeable reward gains. ConicFlow (RH-H) with H=4H=4 matches HardFlow (and improves on it for small deviations), while increasing the horizon to H=8H=8 improves the trade-off across the board, outperforming all baselines and saturating the reward with the least control energy. Neither (RH-H) nor (RH-T) consistently outperforms the other. In short, FMRG does well for small λ and HardFlow for larger λ. ConicFlow with a short horizon combines the strengths of FMRG and HardFlow across all λ, and improves further as H increases. Unconstrained image setting. In high-dimensional, prompt-conditioned image generation, we use FLUX.1-dev with the flow-map LoRA released by (Huang et al., 2026), distilled from the base model as our flow matching model (Boffi et al., 2025). We consider a well-established task: reward-guided image editing (Li et al., 2026; Kim et al., 2022; Patashnik et al., 2021). A base prompt c generates the sample, while the reward r⁡(x,c′)r(x;c ) is evaluated against a different edit prompt c′c . The control must introduce the attribute named by c′c while deviating as little as possible from the original image. We use three commonly used rewards CLIP (Radford et al., 2021), ImageReward (Xu et al., 2023), and PickScore (Kirstain et al., 2023). The RH-T variant of ConicFlow was used across all image experiments as it led to better performance. Figure 2: Trade-off between terminal reward and total control cost. Left: Reward vs. control energy on the editing task, traced by sweeping λ (upper left is better). Right: Edits from each method at approximately matched reward; ConicFlow reaches the same reward with much lower control energy. ConicFlow achieves the best trade-off between reward and control energy (Figure 2, left). At matched reward, its edits stay closer to the original image (Figure 2, right). Table 1 quantifies this across multiple rewards, seeds and prompts, each describing a new optimization problem. Table 1: Mean reward gain over a shared energy budget. Each method’s reward–energy frontier, averaged over the control-energy range attained by all three methods and measured relative to the uncontrolled sample (units are those of rewards). Mean over n=4n=4 (prompt, seed) pairs per reward. Method CLIP (↑)( ) ImageReward (↑)( ) PickScore (↑)( ) ConicFlow 0.130±0.0180.130± 0.018 2.048±0.4552.048± 0.455 6.158±0.7286.158± 0.728 HardFlow 0.118±0.0210.118± 0.021 1.279±0.3741.279± 0.374 1.736±0.4721.736± 0.472 FMRG 0.065±0.0080.065± 0.008 1.278±0.5041.278± 0.504 2.670±0.0732.670± 0.073 Constrained image setting. We consider the same image editing task, but now enforce the desired attribute through a hard constraint r⁡(x,c′)≥τr(x;c )≥τ, i.e., h⁡(x)=τ−r⁡(x,c′)≤0h(x)=τ-r(x;c )≤ 0, using PickScore as r. In the multi-constraint setting, we impose one such constraint per edit prompt, hk(x)=τk−r⁡(x,ck′)≤0h_k(x)= _k-r(x;c _k)≤ 0 for k=1,…,Kk=1,…,K, requiring the sample to exhibit all K attributes simultaneously. For each problem, the threshold is set to τ=m⋅r⁡(xunc,c′)τ=m· r(x^unc;c ), a multiple of the score achieved by the uncontrolled sample, with m ranging from 1.11.1 (easier) to 1.51.5 (harder). We report the control energy and the feasibility rate, i.e., the fraction of (prompt, seed) problems whose final sample satisfies the constraint. Since FMRG has no mechanism for enforcing hard constraints, we do not include it in this comparison. Figure 3: Constrained image editing comparison. Top left: control energy vs. threshold multiplier m, median over (prompt, seed) pairs. Bottom left: feasibility rate vs. m. Middle: single-constraint edits with increasing m. Right: multi-constraint edits with increasing m. As m increases, both methods expend more control energy (Figure 3, top left), but ConicFlow satisfies the constraints with substantially less energy (up to 20×20×). ConicFlow also maintains a feasibility rate close to 100%100\%, whereas HardFlow’s drops as the threshold grows (bottom left). The examples demonstrate the same message: as m increases, both methods express the edit concept more strongly, but HardFlow deviates increasingly from the uncontrolled image, eventually producing visible noise and artifacts. This is most pronounced in the multi-constraint setting (right). See Appendix C for additional implementation details and results. 6 Conclusion We introduced ConicFlow, a receding-horizon primal-dual framework for constrained inference with pretraned flow matching models. This framework combines reward optimization, control-energy regularization, and explicit terminal constraints. We adapt receding-horizon control and propose a terminal-estimate variant that uses the pretrained flow map to evaluate the reward and constraints at the predicted terminal state rather than at the end of the truncated window. Our analysis quantifies discretization and horizon truncation errors and establishes primal-dual convergence guarantees under convex-affine assumptions, including the propagation of local optimization errors across successive windows. Experiments on synthetic data and image editing demonstrate that this approach can improve the reward–control trade-off while maintaining high empirical feasibility. Acknowledgment James Anderson is partially funded by NSF grants EECS 2144634 and CNS 2535097, and the Center of AI Technology (CAIT) in collaboration with Amazon. References Albergo et al. (2025) M. Albergo, N. M. Boffi, and E. Vanden-Eijnden Stochastic interpolants: a unifying framework for flows and diffusions. Journal of Machine Learning Research 26 (209), p. 1–80. Cited by: Appendix A, §1. Alonso et al. (2022) C. A. Alonso, J. S. Li, J. Anderson, and N. Matni Distributed and localized model-predictive control—part i: synthesis and implementation. IEEE Transactions on Control of Network Systems 10 (2), p. 1058–1068. External Links: Document Cited by: Appendix A, §3, §3. Bemporad et al. (2000) A. Bemporad, M. Morari, V. Dua, and E. N. Pistikopoulos The explicit solution of model predictive control via multiparametric quadratic programming. In Proceedings of the 2000 American control conference. ACC (IEEE Cat. No. 00CH36334), Vol. 2, p. 872–876. Cited by: §1. Ben-Hamu et al. (2024) H. Ben-Hamu, O. Puny, I. Gat, B. Karrer, U. Singer, and Y. Lipman D-flow: differentiating through flows for controlled generation. arXiv preprint arXiv:2402.14017. Cited by: Appendix A, §1, §1. Black et al. (2024) K. Black, M. Janner, Y. Du, I. Kostrikov, and S. Levine Training diffusion models with reinforcement learning. In International Conference on Learning Representations, Vol. 2024, p. 4965–4987. Cited by: Appendix A, §1. Boffi et al. (2024) N. M. Boffi, M. S. Albergo, and E. Vanden-Eijnden Flow map matching with stochastic interpolants: a mathematical framework for consistency models. arXiv preprint arXiv:2406.07507. Cited by: Appendix A. Boffi et al. (2025) N. M. Boffi, M. S. Albergo, and E. Vanden-Eijnden How to build a consistency model: learning flow maps via self-distillation. arXiv preprint arXiv:2505.18825. Cited by: Appendix A, §1, §1, §5. Bouvier et al. (2025) J. Bouvier, K. Ryu, K. Nagpal, Q. Liao, K. Sreenath, and N. Mehr DDAT: diffusion policies enforcing dynamically admissible robot trajectories. In Robotics: Science and Systems (RSS), Cited by: Appendix A. Cagienard et al. (2007) R. Cagienard, P. Grieder, E. C. Kerrigan, and M. Morari Move blocking strategies in receding horizon control. Journal of Process Control 17 (6), p. 563–570. External Links: Document Cited by: Appendix A. Christopher et al. (2026) J. K. Christopher, J. E. Warner, and F. Fioretto Constraint-aware flow matching: decision aligned end-to-end training for constrained sampling. arXiv preprint arXiv:2605.12754. Cited by: §1. Christopher et al. (2024) J. K. Christopher, S. Baek, and F. Fioretto Constrained synthesis with projected diffusion models. In Advances in Neural Information Processing Systems (NeurIPS), Vol. 37, p. 89307–89333. Cited by: Appendix A. Chung et al. (2023) H. Chung, J. Kim, M. T. Mccann, M. L. Klasky, and J. C. Ye Diffusion posterior sampling for general noisy inverse problems. In International Conference on Learning Representations (ICLR), Cited by: Appendix A. Dhariwal and Nichol (2021) P. Dhariwal and A. Nichol Diffusion models beat GANs on image synthesis. In Advances in Neural Information Processing Systems (NeurIPS), Vol. 34. Cited by: Appendix A. Domingo i Enrich et al. (2025) C. Domingo i Enrich, M. Drozdzal, B. Karrer, and R. T. Chen Adjoint matching: fine-tuning flow and diffusion generative models with memoryless stochastic optimal control. In International Conference on Learning Representations, Vol. 2025, p. 53791–53846. Cited by: Appendix A, §1. Fan et al. (2023) Y. Fan, O. Watkins, Y. Du, H. Liu, M. Ryu, C. Boutilier, P. Abbeel, M. Ghavamzadeh, K. Lee, and K. Lee DPOK: reinforcement learning for fine-tuning text-to-image diffusion models. In Advances in Neural Information Processing Systems (NeurIPS), Vol. 36. Cited by: Appendix A. Gao et al. (2026) W. Gao, M. Li, and Q. Li Terminally constrained flow-based generative models from an optimal control perspective. arXiv preprint arXiv:2601.09474. Cited by: Appendix A. García et al. (1989) C. E. García, D. M. Prett, and M. Morari Model predictive control: theory and practice—a survey. Automatica 25 (3), p. 335–348. External Links: Document Cited by: Appendix A, §1, §3. Goulart et al. (2006) P. J. Goulart, E. C. Kerrigan, and J. M. Maciejowski Optimization over state feedback policies for robust control with constraints. Automatica 42 (4), p. 523–533. External Links: Document Cited by: Appendix A. Goulart and Kerrigan (2007) P. J. Goulart and E. C. Kerrigan Output feedback receding horizon control of constrained systems. International Journal of Control 80 (1), p. 8–20. External Links: Document Cited by: Appendix A. Ho et al. (2020) J. Ho, A. Jain, and P. Abbeel Denoising diffusion probabilistic models. Advances in neural information processing systems 33, p. 6840–6851. Cited by: §1. Ho and Salimans (2022) J. Ho and T. Salimans Classifier-free diffusion guidance. arXiv preprint arXiv:2207.12598. Cited by: Appendix A. Huang et al. (2026) J. Y. Huang, J. Lin, S. Shah, K. Nair, and N. M. Boffi How to guide your flow: few-step alignment via flow map reward guidance. arXiv preprint arXiv:2604.27147. Cited by: Appendix A, Appendix A, §5, §5, §5. Kim et al. (2022) G. Kim, T. Kwon, and J. C. Ye Diffusionclip: text-guided diffusion models for robust image manipulation. In 2022 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), p. 2416–2425. Cited by: §5. Kirstain et al. (2023) Y. Kirstain, A. Polyak, U. Singer, S. Matiana, J. Penna, and O. Levy Pick-a-pic: an open dataset of user preferences for text-to-image generation. Advances in neural information processing systems 36, p. 36652–36663. Cited by: §C.1.2, §5. Kurtz and Davydov (2026) V. Kurtz and A. Davydov Constrained flow matching via lagrangian dual flows. arXiv preprint arXiv:2607.04513. Cited by: Appendix A, Appendix A. Li et al. (2026) Z. Li, K. Alim, and N. Azizan HardFlow: hard-constrained sampling for flow-matching models via trajectory optimization. IEEE Transactions on Pattern Analysis and Machine Intelligence. Cited by: Appendix A, Appendix A, Appendix A, §5, §5. Lipman et al. (2022) Y. Lipman, R. T. Chen, H. Ben-Hamu, M. Nickel, and M. Le Flow matching for generative modeling. arXiv preprint arXiv:2210.02747. Cited by: Appendix A, §1. Liu et al. (2022) X. Liu, C. Gong, and Q. Liu Flow straight and fast: learning to generate and transfer data with rectified flow. arXiv preprint arXiv:2209.03003. Cited by: Appendix A. Liu et al. (2023) X. Liu, L. Wu, S. Zhang, C. Gong, W. Ping, and Q. Liu FlowGrad: controlling the output of generative ODEs with gradients. In IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), p. 24335–24344. Cited by: Appendix A. Luan et al. (2026) H. Luan, Y. X. Goh, S. Ng, and C. K. Ling Projected coupled diffusion for test-time constrained joint generation. In International Conference on Learning Representations (ICLR), Cited by: Appendix A. Mattingley et al. (2011) J. Mattingley, Y. Wang, and S. Boyd Receding horizon control: automatic generation of high-speed solvers. IEEE Control Systems Magazine 31 (3), p. 52–65. External Links: Document Cited by: Appendix A, §1. Mayne et al. (2000) D. Q. Mayne, J. B. Rawlings, C. V. Rao, and P. O. M. Scokaert Constrained model predictive control: stability and optimality. Automatica 36 (6), p. 789–814. External Links: Document Cited by: Appendix A, §1, §3, §3. Patashnik et al. (2021) O. Patashnik, Z. Wu, E. Shechtman, D. Cohen-Or, and D. Lischinski Styleclip: text-driven manipulation of stylegan imagery. In 2021 IEEE/CVF international conference on computer vision (ICCV), p. 2065–2074. Cited by: §5. Radford et al. (2021) A. Radford, J. W. Kim, C. Hallacy, A. Ramesh, G. Goh, S. Agarwal, G. Sastry, A. Askell, P. Mishkin, J. Clark, et al. Learning transferable visual models from natural language supervision. In International conference on machine learning, p. 8748–8763. Cited by: §C.1.2, §5. Rombach et al. (2022) R. Rombach, A. Blattmann, D. Lorenz, P. Esser, and B. Ommer High-resolution image synthesis with latent diffusion models. In 2022 IEEE/CVF conference on computer vision and pattern recognition (CVPR), p. 10674–10685. Cited by: §1. Römer et al. (2025) R. Römer, A. von Rohr, and A. P. Schoellig Diffusion predictive control with constraints. In Proceedings of the 7th Annual Learning for Dynamics & Control Conference (L4DC), Proceedings of Machine Learning Research, Vol. 283, p. 791–803. Cited by: Appendix A. Schramowski et al. (2023) P. Schramowski, M. Brack, B. Deiseroth, and K. Kersting Safe latent diffusion: mitigating inappropriate degeneration in diffusion models. In 2023 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), p. 22522–22531. Cited by: §1. Uehara et al. (2024) M. Uehara, Y. Zhao, K. Black, E. Hajiramezanali, G. Scalia, N. L. Diamant, A. M. Tseng, T. Biancalani, and S. Levine Fine-tuning of continuous-time diffusion models as entropy-regularized control. arXiv preprint arXiv:2402.15194. Cited by: Appendix A. Wang et al. (2025) L. Wang, C. Cheng, Y. Liao, Y. Qu, and G. Liu Training free guided flow matching with optimal control. In International Conference on Learning Representations (ICLR), Cited by: Appendix A. Wu et al. (2023) X. Wu, Y. Hao, K. Sun, Y. Chen, F. Zhu, R. Zhao, and H. Li Human preference score v2: a solid benchmark for evaluating human preferences of text-to-image synthesis. arXiv preprint arXiv:2306.09341. Cited by: §1. Xu et al. (2023) J. Xu, X. Liu, Y. Wu, Y. Tong, Q. Li, M. Ding, J. Tang, and Y. Dong Imagereward: learning and evaluating human preferences for text-to-image generation. Advances in Neural Information Processing Systems 36, p. 15903–15935. Cited by: §C.1.2, §5. Zampini et al. (2025) S. Zampini, J. K. Christopher, L. Oneto, D. Anguita, and F. Fioretto Training-free constrained generation with stable diffusion models. In Advances in Neural Information Processing Systems (NeurIPS), Cited by: Appendix A. Zhang et al. (2025) J. Zhang, L. Zhao, A. Papachristodoulou, and J. Umenberger Constrained diffusers for safe planning and control. In Advances in Neural Information Processing Systems (NeurIPS), Cited by: Appendix A. Appendix Table of Contents Appendix A Related Work Inference-time control of pretrained diffusion and flow matching models has developed along two largely separate lines: soft guidance methods, which trade reward against fidelity to the pretrained dynamics but offer no feasibility guarantee, and hard-constraint methods, which enforce feasibility, most often by projecting onto the constraint set during sampling. ConicFlow belongs to the second line. It is a hard-constrained sampling method whose terminal constraints are enforced through dual multipliers, and whose objective additionally accommodates a reward and a control-energy penalty so that feasible samples are ranked by preference and stay close to the pretrained dynamics. To evaluate it we compare against a recent state-of-the-art representative from each line: Flow Map Reward Guidance (FMRG) (Huang et al., 2026) for soft guidance, and HardFlow (Li et al., 2026) for hard-constrained sampling in the tradition of projection-based methods. We review both families below and then discuss the control-theoretic tools we build on. Flow matching and flow maps. The pretrained models we steer are trained with flow matching (Lipman et al., 2022; Liu et al., 2022) or, equivalently, stochastic interpolants (Albergo et al., 2025), which learn a velocity field btb_t whose ODE transports a simple base distribution to the data distribution. Flow maps (Boffi et al., 2024; Boffi et al., 2025), which directly predict the ODE solution between two arbitrary times, were introduced to accelerate sampling; we use a pretrained flow map both to integrate the dynamics and, in the RH-T variant, to predict the terminal sample from an intermediate state. Soft guidance and reward-tilted sampling. Classifier guidance (Dhariwal and Nichol, 2021) and classifier-free guidance (Ho and Salimans, 2022) steer diffusion models by adding a (learned or amortized) gradient of a conditional likelihood to the score. Training-free variants replace the classifier with the gradient of a reward evaluated at a one-step estimate of the clean sample, as in diffusion posterior sampling (Chung et al., 2023) and its flow-based counterparts. A second line differentiates through the sampling ODE itself: FlowGrad (Liu et al., 2023) and D-Flow (Ben-Hamu et al., 2024) optimize the initial noise or the whole trajectory by backpropagating a terminal loss, and OC-Flow (Wang et al., 2025) recasts this as an optimal control problem solved via the Pontryagin maximum principle, with convergence guarantees in Euclidean space and on SO⁡(3)SO(3). OC-Flow is the closest precursor to our unconstrained formulation, and it shows that the KL divergence between the controlled and uncontrolled terminal distributions is bounded by the expected control cost, which justifies the quadratic energy term in (P-C); it does not, however, handle constraints, and solving the full-horizon adjoint problem is considerably more expensive than a receding-horizon scheme. Most recently, FMRG (Huang et al., 2026) reformulates guidance as deterministic optimal control and shows that the flow map arises naturally in the optimal feedback, yielding a single-trajectory, training-free method that matches or surpasses multi-particle guidance schemes with as few as three function evaluations at text-to-image scale. FMRG is provably near-optimal in the small-λ regime (small deviations from the pretrained flow), which is exactly the regime in which we observe it performing best (Section 5); it has no mechanism for enforcing hard terminal constraints, and its greedy, myopic control degrades when the reward lies far from the data manifold. We therefore use FMRG as the representative soft-guidance baseline, and, since it subsumes DPS-style methods as special cases, we do not compare against those separately. Fine-tuning for reward alignment. An orthogonal approach adapts the model weights to a reward, via reinforcement learning (Black et al., 2024; Fan et al., 2023) or stochastic optimal control (Domingo i Enrich et al., 2025; Uehara et al., 2024). Fine-tuning amortizes the cost of steering over many samples but must be repeated for every new reward or constraint, can be prohibitively expensive for large models, and offers no per-sample feasibility guarantee. Our method is inference-time only and complementary to these approaches. Hard-constrained sampling via projection. The dominant approach to hard constraints projects the iterate onto the feasible set at every (or every late) sampling step. Projected diffusion models (Christopher et al., 2024) interpret reverse diffusion as projected Langevin dynamics and require a projection operator for the constraint set at each step; Zampini et al. (2025) extend this to latent (Stable Diffusion) models and, notably, to copyright-style constraints on images, but still require a function measuring distance to the feasible set and a projection in latent space. Projected coupled diffusion (Luan et al., 2026) applies the same per-step projection to the joint generation of correlated samples from several pretrained models. In robotics, diffusion predictive control with constraints (Römer et al., 2025) and DDAT (Bouvier et al., 2025) project generated trajectories onto sets that are dynamically admissible for the robot, either at every step or only late in sampling; these methods are tailored to state–action constraints and system dynamics and are less directly comparable to the general terminal-constraint problem we consider. As several authors have noted (Li et al., 2026; Kurtz and Davydov, 2026), projecting the entire sampling path is overly restrictive: intermediate iterates need not be feasible, exact projection onto a general nonlinear set is itself an optimization problem, and projecting inside an ODE produces a discontinuous, projected dynamical system that degrades sample quality. Hard constraints via penalties and dual dynamics. Several recent works instead enforce terminal constraints through the sampling dynamics themselves. TOCFlow (Gao et al., 2026) formulates terminally constrained sampling as an optimal control problem, characterizes the value function through a Hamilton–Jacobi–Bellman equation, and shows that, as the control penalty vanishes, the terminal law converges to a generalized Wasserstein projection onto the constraint manifold x:h⁡(x)=0\x:h(x)=0\. Algorithmically, however, constraints enter through a quadratic terminal penalty ‖h⁡(x)‖2\|h(x)\|^2, so exact feasibility is recovered only in the vanishing-penalty limit; the theoretical analysis is formulated around a smooth zero-level-set constraint manifold. Lagrangian dual flows (Kurtz and Davydov, 2026) augment the generative ODE with a dual co-state (and a slack variable for inequalities) so that constraint violation is driven to zero along the flow without projections or inner optimization; the authors also show why quadratic penalties of the TOCFlow type are a poor substitute for genuine dual variables. This is conceptually close to our dynamic primal–dual scheme, but LDF flows the multiplier continuously alongside the sample with no notion of look-ahead, and has so far been demonstrated only on low-dimensional and MNIST-scale problems. Constrained diffusers (Zhang et al., 2025) bring projected, primal–dual, and augmented-Lagrangian updates to reverse Langevin sampling for safe planning; they target a different, distribution-level problem, maxλ≥0minqKL(q∥p)+λ⊤x∼q[g(x)] _λ≥ 0 _qKL(q\,\|\,p)+λ E_x q[g(x)], whose dual update requires an expectation over samples at every step, and their setting is diffusion-based trajectory planning rather than deterministic flows. Hard constraints via trajectory optimization. HardFlow (Li et al., 2026) is the work most closely related to ours. It reformulates hard-constrained sampling from a flow matching model as a trajectory optimization problem in which constraints are imposed only at the terminal time, uses model predictive control ideas to decompose the long-horizon problem into a sequence of single-step subproblems, allows integral costs that limit distribution shift, and provides a control-theoretic bound on the gap between its surrogate and the ideal problem. It substantially outperforms projection-based methods on robotics, PDE boundary-control, and text-guided image-editing tasks, and it is therefore our representative of the hard-constraint family. The differences from ConicFlow are, first, that HardFlow optimizes a single step at a time and solves each constrained subproblem with an off-the-shelf solver, whereas we optimize over a horizon H≥1H≥ 1 with a dynamic primal–dual scheme whose multipliers are warm-started across windows; second, that HardFlow evaluates the constraint at the end of the current step, whereas our RH-T variant evaluates it at the predicted terminal point; and third, that we provide truncation and finite-iteration guarantees for both variants. Empirically, HardFlow is competitive precisely when large deviations from the pretrained flow are required, complementing FMRG’s strength in the small-deviation regime, and ConicFlow combines the strengths of both (Section 5). Receding-horizon and constrained control. Receding-horizon control provides a classical framework for handling finite-horizon objectives and hard constraints through repeated online optimization (García et al., 1989; Mayne et al., 2000). Several works have addressed its computational and robustness limitations. In particular, Kerrigan and collaborators studied move blocking to reduce the online decision dimension while retaining feasibility and stability guarantees Cagienard et al. (2007), convex parametrizations of feedback policies for robust constrained finite-horizon control (Goulart et al., 2006), and constrained output-feedback receding-horizon control under process and measurement uncertainty (Goulart and Kerrigan, 2007). Related work has developed high-speed RHC solvers Mattingley et al. (2011) and scalable localized MPC architectures Alonso et al. (2022). Our setting differs in that the controlled dynamics are given by a pretrained generative flow and the objective and hard constraints are imposed on the generated terminal sample. Appendix B Proofs B.1 Proof of Lemma 1 Proof. Proof strategy. We prove the result by comparing the continuous and discretized optimal values in both directions. The argument proceeds in four steps. First, the quadratic control penalty together with the Lipschitz terminal reward yields a coercive lower bound on the objective. Using the finite-cost feasible trajectory from Assumption 1 i, this gives a uniform energy bound on all controls relevant to the comparison and, consequently, a common bounded region containing their trajectories. Second, we show that replacing a continuous control by its block averages and propagating the dynamics with Euler discretization changes the terminal state by O⁡(Δt)O( t). Conversely, the piecewise-constant interpolation of a discrete control incurs the same order of terminal-state error. Third, the resulting O⁡(Δt)O( t) terminal constraint violation is corrected using Assumption 1 i, requiring only an O⁡(Δt)O( t) perturbation of the control. This construction also establishes feasibility and existence of the discrete optimum for sufficiently large N. Finally, we apply these constructions in both directions. Discretizing and repairing the continuous optimum gives an upper bound on the discrete optimal value, while interpolating and repairing the discrete optimum gives the reverse inequality. Combining the two yields |JD(UD⋆(N))−JC(uC⋆)|=O⁡(Δt)=O⁡(1N). |J_D (U_D (N) )-J_C(u_C ) |=O( t)=O\! ( 1N ). We next make the quantitative content of Assumptions 1 from i to i explicit and introduce the constants used in the proof. Quantitative form of the assumptions. By Assumption 1 i, there exist constants Lb,Lbt,Lr,Lh>0L_b,L_b^t,L_r,L_h>0 such that ‖bt(x)−bs(y)‖≤Lb‖x−y‖+Lbt|t−s|,\|b_t(x)-b_s(y)\|≤ L_b\|x-y\|+L_b^t|t-s|, and |r⁡(x)−r⁡(y)|≤Lr‖x−y‖,‖h⁡(x)−h⁡(y)‖≤Lh‖x−y‖.|r(x)-r(y)|≤ L_r\|x-y\|, \|h(x)-h(y)\|≤ L_h\|x-y\|. We also define B0:=supt∈[0,1]‖bt(x0)‖<∞.B_0:= _t∈[0,1]\|b_t(x_0)\|<∞. By Assumption 1 i, there exists a feasible continuous-time control u0u^0 satisfying h⁡(x1u0)≤0,JC(u0)<∞,h(x_1^u^0)≤ 0, J_C(u^0)<∞, and the continuous problem admits an optimal control uC⋆u_C . Assumption 1 i is used in the following quantitative form. On the energy ball introduced below, there exist constants κ,ρ0>0κ, _0>0, independent of N, such that any continuous or discrete control whose terminal violation is at most ρ0 _0 can be made feasible by a control perturbation whose norm is at most κ times the violation. For continuous controls we use the L2L^2 norm, while for discrete controls we use ‖U‖ℓΔt22:=∑k=0N−1Δt‖uk‖2.\|U\|_ ^2_ t^2:= _k=0^N-1 t\,\|u_k\|^2. A priori energy bound. We first derive the energy radius on which the remainder of the proof is carried out. For any continuous control u, the dynamics and Assumption 1 i imply ‖xtu−x0‖≤B0t+Lb∫0t‖xsu−x0‖s+∫0t‖us‖s.\|x_t^u-x_0\|≤ B_0t+L_b _0^t\|x_s^u-x_0\|\,ds+ _0^t\|u_s\|\,ds. By Cauchy–Schwarz and Grönwall’s inequality, sup0≤t≤1‖xtu−x0‖≤eLb(B0+‖u‖L2). _0≤ t≤ 1\|x_t^u-x_0\|≤ e^L_b (B_0+\|u\|_L^2 ). Hence, by Lipschitz continuity of r, r⁡(x1u)≤|r⁡(x0)|+LreLb(B0+‖u‖L2).r(x_1^u)≤|r(x_0)|+L_re^L_b (B_0+\|u\|_L^2 ). Define ar:=LreLb,βr:=|r⁡(x0)|+arB0+λar2.a_r:=L_re^L_b, _r:=|r(x_0)|+a_rB_0+λ a_r^2. Young’s inequality gives ar‖u‖L2≤14λ‖u‖L22+λar2.a_r\|u\|_L^2≤ 14λ\|u\|_L^2^2+λ a_r^2. Therefore, JC(u)≥14λ‖u‖L22−βr.J_C(u)≥ 14λ\|u\|_L^2^2- _r. The same argument applied to the Euler recursion gives JD(U)≥14λ‖U‖ℓΔt22−βr.J_D(U)≥ 14λ\|U\|_ ^2_ t^2- _r. Define R⋆:=4λ(JC(u0)+βr),R¯:=R⋆+1.R_ :=4λ (J_C(u^0)+ _r ), R:=R_ +1. Since uC⋆u_C is optimal and u0u^0 is feasible, JC(uC⋆)≤JC(u0).J_C(u_C )≤ J_C(u^0). Using the coercivity estimate above, ‖uC⋆‖L22≤R⋆<R¯.\|u_C \|_L^2^2≤ R_ < R. Step 1: Uniform confinement of relevant trajectories. Consider any continuous control satisfying ‖u‖L22≤R¯.\|u\|_L^2^2≤ R. The preceding trajectory estimate gives sup0≤t≤1‖xtu−x0‖≤(B0+R¯)eLb. _0≤ t≤ 1\|x_t^u-x_0\|≤ (B_0+ R )e^L_b. Define ρR¯:=(B0+R¯)eLb, _ R:= (B_0+ R )e^L_b, and MR¯:=B0+LbρR¯.M_ R:=B_0+L_b _ R. Then ‖bt(xtu)‖≤MR¯.\|b_t(x_t^u)\|≤ M_ R. The same conclusion holds for Euler trajectories. Let U=(u0,…,uN−1)U=(u_0,…,u_N-1) satisfy ‖U‖ℓΔt22≤R¯.\|U\|_ ^2_ t^2≤ R. The Euler recursion gives ‖xk+1U−x0‖≤(1+LbΔt)‖xkU−x0‖+B0Δt+Δt‖uk‖.\|x_k+1^U-x_0\|≤(1+L_b t)\|x_k^U-x_0\|+B_0 t+ t\,\|u_k\|. Moreover, ∑k=0N−1Δt‖uk‖≤(∑k=0N−1Δt‖uk‖2)1/2≤R¯. _k=0^N-1 t\,\|u_k\|≤ ( _k=0^N-1 t\,\|u_k\|^2 )^1/2≤ R. Iterating the recursion and using (1+LbΔt)N≤eLb(1+L_b t)^N≤ e^L_b gives max0≤k≤N⁡‖xkU−x0‖≤ρR¯. _0≤ k≤ N\|x_k^U-x_0\|≤ _ R. Thus, both continuous and Euler trajectories generated by controls in the energy ball R¯ R remain uniformly bounded, and the drift along all such trajectories is bounded by MR¯M_ R. Step 2: Continuous–discrete approximation. We next establish the approximation estimates required in the remainder of the proof. Continuous-to-discrete approximation. Let u∈L2([0,1],ℝd)u∈ L^2([0,1];R^d) satisfy ‖u‖L22≤R¯,\|u\|_L^2^2≤ R, and define its block averages by u¯k:=1Δt∫tktk+1uss. u_k:= 1 t _t_k^t_k+1u_s\,ds. By Jensen’s inequality, ∑k=0N−1Δt‖u¯k‖2≤∫01‖ut‖2t≤R¯. _k=0^N-1 t\,\| u_k\|^2≤ _0^1\|u_t\|^2\,dt≤ R. Let xu(t)x^u(t) denote the exact continuous trajectory and let xkΔx_k denote the Euler trajectory driven by U¯=(u¯0,…,u¯N−1) U=( u_0,…, u_N-1): xk+1Δ=xkΔ+Δt(btk(xkΔ)+u¯k).x_k+1 =x_k + t (b_t_k(x_k )+ u_k ). Set ek:=xu(tk)−xkΔ.e_k:=x^u(t_k)-x_k . By construction, ∫tktk+1(us−u¯k)s=0. _t_k^t_k+1 (u_s- u_k )\,ds=0. Thus the integrated control contribution is reproduced exactly over each grid interval. Subtracting the Euler recursion from the exact dynamics gives ek+1=ek+Δt(btk(xu(tk))−btk(xkΔ))+τk,e_k+1=e_k+ t (b_t_k(x^u(t_k))-b_t_k(x_k ) )+ _k, where τk:=∫tktk+1[bs(xu(s))−btk(xu(tk))]s. _k:= _t_k^t_k+1 [b_s(x^u(s))-b_t_k(x^u(t_k)) ]ds. Define qk:=∫tktk+1‖us‖2s.q_k:= _t_k^t_k+1\|u_s\|^2\,ds. For s∈[tk,tk+1]s∈[t_k,t_k+1], Step 1 and Cauchy–Schwarz imply ‖xu(s)−xu(tk)‖≤MR¯(s−tk)+s−tkqk.\|x^u(s)-x^u(t_k)\|≤ M_ R(s-t_k)+ s-t_k\, q_k. Therefore, ‖τk‖≤LbMR¯Δt2+12LbtΔt2+LbΔt3/2qk.\| _k\|≤ L_bM_ R t^2+ 12L_b^t t^2+L_b t^3/2 q_k. Moreover, ∑k=0N−1qk≤N(∑k=0N−1qk)1/2≤R¯Δt. _k=0^N-1 q_k≤ N ( _k=0^N-1q_k )^1/2≤ R t. Consequently, ∑k=0N−1‖τk‖≤(LbMR¯+12Lbt+LbR¯)Δt. _k=0^N-1\| _k\|≤ (L_bM_ R+ 12L_b^t+L_b R ) t. Since ‖ek+1‖≤(1+LbΔt)‖ek‖+‖τk‖,\|e_k+1\|≤(1+L_b t)\|e_k\|+\| _k\|, discrete Grönwall yields max0≤k≤N⁡‖xu(tk)−xkΔ‖≤C0Δt, _0≤ k≤ N\|x^u(t_k)-x_k \|≤ C_0 t, where C0:=eLb(LbMR¯+12Lbt+LbR¯).C_0:=e^L_b (L_bM_ R+ 12L_b^t+L_b R ). Discrete-to-continuous interpolation. Conversely, let U=(u0,…,uN−1)U=(u_0,…,u_N-1) satisfy ‖U‖ℓΔt22≤R¯,\|U\|_ ^2_ t^2≤ R, and define its piecewise-constant interpolation u^(t):=uk,t∈[tk,tk+1). u(t):=u_k, t∈[t_k,t_k+1). Then ‖u^‖L22=‖U‖ℓΔt22.\| u\|_L^2^2=\|U\|_ ^2_ t^2. Since the block average of u u on each interval is exactly uku_k, the preceding argument applies directly and gives max0≤k≤N⁡‖x^(tk)−xkU‖≤C0Δt, _0≤ k≤ N\| x(t_k)-x_k^U\|≤ C_0 t, where x x is the exact continuous trajectory generated by u u. Sensitivity to control perturbations. Let u,v∈L2([0,1],ℝd)u,v∈ L^2([0,1];R^d). Then ‖xtu−xtv‖≤Lb∫0t‖xsu−xsv‖s+∫0t‖us−vs‖s.\|x_t^u-x_t^v\|≤ L_b _0^t\|x_s^u-x_s^v\|\,ds+ _0^t\|u_s-v_s\|\,ds. Grönwall and Cauchy–Schwarz imply sup0≤t≤1‖xtu−xtv‖≤eLb‖u−v‖L2. _0≤ t≤ 1\|x_t^u-x_t^v\|≤ e^L_b\|u-v\|_L^2. Similarly, for two discrete controls U and V, max0≤k≤N⁡‖xkU−xkV‖≤eLb‖U−V‖ℓΔt2. _0≤ k≤ N\|x_k^U-x_k^V\|≤ e^L_b\|U-V\|_ ^2_ t. Thus Step 2 provides the two estimates used below: Euler endpoint error=O⁡(Δt),control perturbation δ⟹state perturbation O(δ).Euler endpoint error=O( t), perturbation δ perturbation O(δ). Step 3: Discrete feasibility, existence, and energy bound. Let U¯0 U^0 denote the block average of the feasible reference control u0u^0. From the coercivity estimate, ‖u0‖L22≤R⋆,\|u^0\|_L^2^2≤ R_ , and Jensen’s inequality gives ‖U¯0‖ℓΔt22≤R⋆<R¯.\| U^0\|_ ^2_ t^2≤ R_ < R. By Step 2, ‖xNU¯0−x1u0‖≤C0Δt.\|x_N U^0-x_1^u^0\|≤ C_0 t. Since u0u^0 is feasible and h is LhL_h-Lipschitz, ‖[h⁡(xNU¯0)]+‖≤LhC0Δt.\|[h(x_N U^0)]_+\|≤ L_hC_0 t. Define N0:=⌈max⁡1,LhC0ρ0⌉.N_0:= \1, L_hC_0 _0 \ . For N≥N0N≥ N_0, Assumption 1 i provides a feasible discrete control U~0 U^0 satisfying ‖U~0−U¯0‖ℓΔt2≤δN,\| U^0- U^0\|_ ^2_ t≤ _N, where δN:=κLhC0Δt. _N:=κ L_hC_0 t. The change in quadratic control cost satisfies 12λ(‖U~0‖ℓΔt22−‖U¯0‖ℓΔt22)≤R¯λδN+12λδN2. 12λ (\| U^0\|_ ^2_ t^2-\| U^0\|_ ^2_ t^2 )≤ Rλ _N+ 12λ _N^2. Moreover, the control-sensitivity estimate and Lipschitz continuity of r give |r⁡(xNU~0)−r⁡(xNU¯0)|≤LreLbδN. |r(x_N U^0)-r(x_N U^0) |≤ L_re^L_b _N. Using also the Euler endpoint estimate and the fact that block averaging cannot increase the quadratic control cost, JD(U~0)≤JC(u0)+LrC0Δt+(R¯λ+LreLb)δN+12λδN2.J_D( U^0)≤ J_C(u^0)+L_rC_0 t+ ( Rλ+L_re^L_b ) _N+ 12λ _N^2. Since Δt2≤Δt t^2≤ t, define Cdisc:=C0[Lr+(R¯λ+LreLb)κLh+κ2Lh2C02λ].C_disc:=C_0 [L_r+ ( Rλ+L_re^L_b )κ L_h+ κ^2L_h^2C_02λ ]. Then JD(U~0)≤JC(u0)+CdiscΔt.J_D( U^0)≤ J_C(u^0)+C_disc t. Thus the discrete feasible set is nonempty for all N≥N0N≥ N_0. For fixed N, the map U↦xNUU x_N^U is continuous, and therefore the feasible set is closed. Moreover, JD(U)≥14λ‖U‖ℓΔt22−βr,J_D(U)≥ 14λ\|U\|_ ^2_ t^2- _r, so JDJ_D is coercive. Since the discrete control space is finite-dimensional, the discrete problem admits a global minimizer UD⋆(N)U_D (N). By optimality, JD(UD⋆(N))≤JD(U~0)≤JC(u0)+CdiscΔt.J_D (U_D (N) )≤ J_D( U^0)≤ J_C(u^0)+C_disc t. Using the coercivity estimate again, ‖UD⋆(N)‖ℓΔt22≤R⋆+4λCdiscΔt.\|U_D (N)\|_ ^2_ t^2≤ R_ +4λ C_disc t. Define N2:=⌈max⁡1,4λCdisc⌉,N1:=max⁡N0,N2.N_2:= \1,4λ C_disc \ , N_1:= \N_0,N_2\. Then, for every N≥N1N≥ N_1, ‖UD⋆(N)‖ℓΔt22≤R⋆+1=R¯.\|U_D (N)\|_ ^2_ t^2≤ R_ +1= R. Step 4: Two-sided comparison of the optimal values. We now apply the preceding constructions to the two optimal controls. Continuous optimum to discrete feasible control. Since ‖uC⋆‖L22≤R⋆<R¯,\|u_C \|_L^2^2≤ R_ < R, let U¯⋆ U denote the block average of uC⋆u_C . By Jensen’s inequality, ‖U¯⋆‖ℓΔt22≤R⋆.\| U \|_ ^2_ t^2≤ R_ . Step 2 gives ‖xNU¯⋆−x1uC⋆‖≤C0Δt.\|x_N U -x_1^u_C \|≤ C_0 t. Since uC⋆u_C is feasible, ‖[h⁡(xNU¯⋆)]+‖≤LhC0Δt.\|[h(x_N U )]_+\|≤ L_hC_0 t. For N≥N0N≥ N_0, terminal repairability provides a feasible discrete control U~⋆ U satisfying ‖U~⋆−U¯⋆‖ℓΔt2≤δN.\| U - U \|_ ^2_ t≤ _N. Repeating the preceding objective estimate gives JD(U~⋆)≤JC(uC⋆)+CdiscΔt.J_D( U )≤ J_C(u_C )+C_disc t. Since U~⋆ U is feasible for the discrete problem, JD(UD⋆(N))−JC(uC⋆)≤CdiscΔt.J_D (U_D (N) )-J_C(u_C )≤ C_disc t. Discrete optimum to continuous feasible control. For N≥N1N≥ N_1, Step 3 gives ‖UD⋆(N)‖ℓΔt22≤R¯.\|U_D (N)\|_ ^2_ t^2≤ R. Define its piecewise-constant interpolation by u^D⋆(t):=uD,k⋆,t∈[tk,tk+1). u_D (t):=u_D,k , t∈[t_k,t_k+1). The quadratic control energy is preserved exactly: ‖u^D⋆‖L22=‖UD⋆(N)‖ℓΔt22.\| u_D \|_L^2^2=\|U_D (N)\|_ ^2_ t^2. Let x^D⋆ x_D denote the corresponding exact continuous trajectory. Step 2 gives ‖x^D⋆(1)−xNUD⋆(N)‖≤C0Δt.\| x_D (1)-x_N^U_D (N)\|≤ C_0 t. Since UD⋆(N)U_D (N) is feasible for the discrete problem, ‖[h⁡(x^D⋆(1))]+‖≤LhC0Δt.\|[h( x_D (1))]_+\|≤ L_hC_0 t. For N≥N0N≥ N_0, Assumption 1 i provides a feasible continuous control u~D⋆ u_D satisfying ‖u~D⋆−u^D⋆‖L2≤δN.\| u_D - u_D \|_L^2≤ _N. The piecewise-constant interpolation preserves the quadratic running cost exactly, and the same endpoint and repair estimates therefore yield JC(u~D⋆)≤JD(UD⋆(N))+CdiscΔt.J_C( u_D )≤ J_D (U_D (N) )+C_disc t. Since u~D⋆ u_D is feasible for the continuous problem, optimality of uC⋆u_C implies JC(uC⋆)−JD(UD⋆(N))≤CdiscΔt.J_C(u_C )-J_D (U_D (N) )≤ C_disc t. Combining the two directions gives |JD(UD⋆(N))−JC(uC⋆)|≤CdiscΔt=CdiscN, |J_D (U_D (N) )-J_C(u_C ) |≤ C_disc t= C_discN, for every N≥N1N≥ N_1. This proves the claim. ∎ Remark 2 (What determines the discretization rate). The proof separates the error into two mechanisms. Regularity of the flow, together with block averaging of the additive control, gives the first-order Euler endpoint estimate ‖x1u−xNU¯‖=O⁡(Δt).\|x_1^u-x_N U\|=O( t). Local terminal repairability then converts the resulting O⁡(Δt)O( t) constraint violation into an O⁡(Δt)O( t) perturbation of the control. Consequently, enforcing exact terminal feasibility does not change the first-order discretization rate. B.2 Unconstrained Receding-Horizon Bounds For completeness, we record the unconstrained counterparts of the constrained receding-horizon bounds in the main text. These results are special cases in which the terminal feasibility constraint is absent. Lemma 5 (Unconstrained RH-H bound). Suppose the discretized dynamics are affine, xk+1=xk+Δt(Akxk+ck+uk),Δt=1N,x_k+1=x_k+ t(A_kx_k+c_k+u_k), t= 1N, with ‖Akx+ck‖≤M\|A_kx+c_k\|≤ M on the relevant reachable set, and let r be concave and LrL_r-Lipschitz. For 1≤H≤N1≤ H≤ N, let H,⋆ x^H, denote the exact unconstrained RH-H closed-loop trajectory with horizon H, and let ⋆ x denote the full-horizon optimum. Then 0≤J⁡(H,⋆)−J⁡(⋆)≤CHUN−HN,0≤ J( x^H, )-J( x )≤ C_H U N-HN, where CHU:=2max⁡LrM+λLr22,M22λ.C_H U:=2 \L_rM+ λ L_r^22, M^22λ \. In particular, CHUC_H U is independent of N and H. Proof. Recall that Jj,L⋆(x):=minuj:j+L−1∑k=j+L−1Δt2λ∥uk∥2−r(xj+L),J_j,L (x):= _u_j:j+L-1 \ _k=j^j+L-1 t2λ\|u_k\|^2-r(x_j+L) \, subject to xk+1=xk+Δt(Akxk+ck+uk),xj=x.x_k+1=x_k+ t(A_kx_k+c_k+u_k), x_j=x. We write Jj⋆(x):=Jj,N−j⋆(x)J_j (x):=J_j,N-j (x) for the full remaining-horizon value. Under the affine dynamics, the terminal state is affine in the control sequence. Since r is concave, −r-r is convex, while the quadratic control cost is coercive. Hence all finite-horizon unconstrained problems considered below attain their minima. Step 1: Local sensitivity to the terminal horizon. Fix j,L,qj,L,q such that j+L+q≤Nj+L+q≤ N, and set τ:=qΔt.τ:=q t. We first show that |Jj,L+q⋆(x)−Jj,L⋆(x)|≤Clocτ, |J_j,L+q (x)-J_j,L (x) |≤ C_locτ, where Cloc:=max⁡LrM+λLr22,M22λ.C_loc:= \L_rM+ λ L_r^22, M^22λ \. Let uj:j+L+q−1⋆u_j:j+L+q-1 be optimal for Jj,L+q⋆(x)J_j,L+q (x), with associated trajectory xk⋆x_k , and denote the energy of its last q controls by E:=∑k=j+Lj+L+q−1Δt‖uk⋆‖2.E:= _k=j+L^j+L+q-1 t\,\|u_k \|^2. Using only the first L controls gives an admissible candidate for Jj,L⋆(x)J_j,L (x). Therefore, Jj,L⋆(x)−Jj,L+q⋆(x) J_j,L (x)-J_j,L+q (x) ≤r⁡(xj+L+q⋆)−r⁡(xj+L⋆)−E2λ. ≤ r(x_j+L+q )-r(x_j+L )- E2λ. By the LrL_r-Lipschitz continuity of r, Jj,L⋆(x)−Jj,L+q⋆(x)≤Lr‖xj+L+q⋆−xj+L⋆‖−E2λ.J_j,L (x)-J_j,L+q (x)≤ L_r\|x_j+L+q -x_j+L \|- E2λ. Using the dynamics and the drift bound ‖Akxk+ck‖≤M\|A_kx_k+c_k\|≤ M from Assumption 2, ‖xj+L+q⋆−xj+L⋆‖ \|x_j+L+q -x_j+L \| ≤∑k=j+Lj+L+q−1Δt(M+‖uk⋆‖) ≤ _k=j+L^j+L+q-1 t (M+\|u_k \| ) ≤Mτ+τE, ≤ Mτ+ τ E, where the second inequality follows from Cauchy–Schwarz. Hence Jj,L⋆(x)−Jj,L+q⋆(x)≤LrMτ+LrτE−E2λ.J_j,L (x)-J_j,L+q (x)≤ L_rMτ+L_r τ E- E2λ. Completing the square in E E gives LrτE−E2λ≤λLr22τ,L_r τ E- E2λ≤ λ L_r^22τ, and therefore Jj,L⋆(x)−Jj,L+q⋆(x)≤(LrM+λLr22)τ.J_j,L (x)-J_j,L+q (x)≤ (L_rM+ λ L_r^22 )τ. For the reverse direction, let uj:j+L−1⋆u_j:j+L-1 be optimal for Jj,L⋆(x)J_j,L (x) and let x¯:=xj+L⋆ x:=x_j+L be its terminal state. Extend this trajectory over the next q steps by choosing uk=−(Akx¯+ck),k=j+L,…,j+L+q−1.u_k=- (A_k x+c_k ), k=j+L,…,j+L+q-1. Then xk+1=x¯,x_k+1= x, so the state, and hence the terminal reward, remains unchanged throughout the extension. Since ‖Akx¯+ck‖≤M\|A_k x+c_k\|≤ M, the additional control cost is at most ∑k=j+Lj+L+q−1Δt2λ‖uk‖2≤M22λτ. _k=j+L^j+L+q-1 t2λ\|u_k\|^2≤ M^22λτ. Thus, Jj,L+q⋆(x)−Jj,L⋆(x)≤M22λτ.J_j,L+q (x)-J_j,L (x)≤ M^22λτ. Combining the two directions yields |Jj,L+q⋆(x)−Jj,L⋆(x)|≤ClocqΔt. |J_j,L+q (x)-J_j,L (x) |≤ C_locq t. Step 2: Accumulation along the RH-H trajectory. Set m:=N−H,m:=N-H, so that mΔt=(N−H)/Nm t=(N-H)/N. For j=0,…,mj=0,…,m, define Wj(x):=Jj,H⋆(x).W_j(x):=J_j,H (x). At j=mj=m, the horizon H reaches the true terminal time, and therefore Wm(x)=Jm⋆(x).W_m(x)=J_m (x). Let xjH,⋆x_j^H, denote the RH-H closed-loop state and ujH,⋆u_j^H, the first control returned by the H-step problem at time j. By optimality of the truncated problem, Wj(xjH,⋆)=Δt2λ‖ujH,⋆‖2+Jj+1,H−1⋆(xj+1H,⋆),W_j(x_j^H, )= t2λ\|u_j^H, \|^2+J_j+1,H-1 (x_j+1^H, ), for j=0,…,m−1j=0,…,m-1. Therefore, Δt2λ‖ujH,⋆‖2 t2λ\|u_j^H, \|^2 =Wj(xjH,⋆)−Wj+1(xj+1H,⋆) =W_j(x_j^H, )-W_j+1(x_j+1^H, ) +[Wj+1(xj+1H,⋆)−Jj+1,H−1⋆(xj+1H,⋆)]. + [W_j+1(x_j+1^H, )-J_j+1,H-1 (x_j+1^H, ) ]. Applying the local horizon estimate with q=1q=1 gives Wj+1(xj+1H,⋆)−Jj+1,H−1⋆(xj+1H,⋆)≤ClocΔt.W_j+1(x_j+1^H, )-J_j+1,H-1 (x_j+1^H, )≤ C_loc t. Hence Δt2λ‖ujH,⋆‖2≤Wj(xjH,⋆)−Wj+1(xj+1H,⋆)+ClocΔt. t2λ\|u_j^H, \|^2≤ W_j(x_j^H, )-W_j+1(x_j+1^H, )+C_loc t. Summing from j=0j=0 to m−1m-1 yields ∑j=0m−1Δt2λ‖ujH,⋆‖2≤W0(x0)−Wm(xmH,⋆)+mClocΔt. _j=0^m-1 t2λ\|u_j^H, \|^2≤ W_0(x_0)-W_m(x_m^H, )+mC_loc t. From time m=N−Hm=N-H onward, the remaining horizon is at most H, so RH-H coincides with full remaining-horizon re-optimization. Consequently, its remaining cost from xmH,⋆x_m^H, is exactly Jm⋆(xmH,⋆)=Wm(xmH,⋆)J_m (x_m^H, )=W_m(x_m^H, ). Thus J⁡(H,⋆)≤W0(x0)+mClocΔt.J( x^H, )≤ W_0(x_0)+mC_loc t. Finally, W0(x0)=J0,H⋆(x0)W_0(x_0)=J_0,H (x_0), while J⁡(⋆)=J0,N⋆(x0)J( x )=J_0,N (x_0). Applying the local horizon estimate once more, now with L=HL=H and q=mq=m, gives W0(x0)−J⁡(⋆)≤ClocmΔt.W_0(x_0)-J( x )≤ C_locm t. Therefore, J⁡(H,⋆)−J⁡(⋆)≤2ClocmΔt.J( x^H, )-J( x )≤ 2C_locm t. Since ⋆ x is globally optimal, J⁡(H,⋆)−J⁡(⋆)≥0.J( x^H, )-J( x )≥ 0. Using m=N−Hm=N-H and Δt=1/N t=1/N, 0≤J⁡(H,⋆)−J⁡(⋆)≤2ClocN−HN.0≤ J( x^H, )-J( x )≤ 2C_loc N-HN. Defining CHU:=2Cloc=2max⁡LrM+λLr22,M22λC_H U:=2C_loc=2 \L_rM+ λ L_r^22, M^22λ \ proves the result. ∎ Lemma 6 (Unconstrained RH-T bound). Suppose the hypotheses of Lemma 5 hold and, in addition, κA:=sup0≤k≤ℓ≤N‖∏s=kℓ−1(I+ΔtAs)‖<∞. _A:= _0≤ k≤ ≤ N \| _s=k -1(I+ tA_s) \|<∞. For 1≤H≤N1≤ H≤ N, let ^H,⋆ x^H, denote the exact unconstrained RH-T closed-loop trajectory obtained by optimizing controls over an H-step window and propagating the local endpoint to time N under the uncontrolled pretrained flow. Then 0≤J⁡(^H,⋆)−J⁡(⋆)≤CTUN−HN,0≤ J( x^H, )-J( x )≤ C_T U N-HN, where CTU:=λ2κA2Lr2.C_T U:= λ2 _A^2L_r^2. In particular, CTUC_T U is independent of N and H. Proof. As in the proof of Lemma 5, the quadratic control penalty and Lipschitz terminal reward ensure that all finite-dimensional optimization problems below attain their minima. We therefore focus only on the additional structure introduced by terminal rollout. For j+L≤Nj+L≤ N, define the rollout value J^j,L⋆(x):=minuj:j+L−1∑k=j+L−1Δt2λ∥uk∥2−r(x^N), J_j,L (x):= _u_j:j+L-1 \ _k=j^j+L-1 t2λ\|u_k\|^2-r( x_N) \, where the state follows xk+1=xk+Δt(Akxk+ck+uk)x_k+1=x_k+ t(A_kx_k+c_k+u_k) for k=j,…,j+L−1k=j,…,j+L-1, and is subsequently propagated to time N under the uncontrolled flow, i.e., uk=0u_k=0 for k≥j+Lk≥ j+L. Step 1: Terminal rollout does not accumulate local horizon errors. Increasing the control horizon only enlarges the feasible set of control sequences. Hence Jj⋆(x)≤J^j,L+1⋆(x)≤J^j,L⋆(x).J_j (x)≤ J_j,L+1 (x)≤ J_j,L (x). Let m:=N−H,m:=N-H, and define W^j(x):=J^j,H⋆(x),j=0,…,m. W_j(x):= J_j,H (x), j=0,…,m. At j=mj=m, exactly H steps remain until the true terminal time, and therefore W^m(x)=Jm⋆(x). W_m(x)=J_m (x). Let x^jH,⋆ x_j^H, denote the RH-T closed-loop state and u^jH,⋆ u_j^H, the first control selected at time j. By optimality of the H-step rollout problem, W^j(x^jH,⋆)=Δt2λ‖u^jH,⋆‖2+J^j+1,H−1⋆(x^j+1H,⋆). W_j( x_j^H, )= t2λ\| u_j^H, \|^2+ J_j+1,H-1 ( x_j+1^H, ). Since allowing one additional control action cannot increase the optimal value, W^j+1(x^j+1H,⋆)=J^j+1,H⋆(x^j+1H,⋆)≤J^j+1,H−1⋆(x^j+1H,⋆). W_j+1( x_j+1^H, )= J_j+1,H ( x_j+1^H, )≤ J_j+1,H-1 ( x_j+1^H, ). Consequently, Δt2λ‖u^jH,⋆‖2≤W^j(x^jH,⋆)−W^j+1(x^j+1H,⋆). t2λ\| u_j^H, \|^2≤ W_j( x_j^H, )- W_j+1( x_j+1^H, ). Summing from j=0j=0 to m−1m-1 gives ∑j=0m−1Δt2λ‖u^jH,⋆‖2≤W^0(x0)−W^m(x^mH,⋆). _j=0^m-1 t2λ\| u_j^H, \|^2≤ W_0(x_0)- W_m( x_m^H, ). From time m onward, the control horizon reaches the true terminal time, so the remaining closed-loop cost is Jm⋆(x^mH,⋆)=W^m(x^mH,⋆)J_m ( x_m^H, )= W_m( x_m^H, ). Therefore, J⁡(^H,⋆)≤J^0,H⋆(x0).J( x^H, )≤ J_0,H (x_0). Since the full-horizon problem has a larger control set, J⁡(⋆)=J0⋆(x0)≤J^0,H⋆(x0).J( x )=J_0 (x_0)≤ J_0,H (x_0). It is therefore sufficient to bound the single initial rollout gap J^0,H⋆(x0)−J0⋆(x0) J_0,H (x_0)-J_0 (x_0). Step 2: Bounding the initial rollout gap. Let u0⋆,…,uN−1⋆u_0 ,…,u_N-1 denote the controls associated with the full-horizon optimum ⋆ x . Define τH:=(N−H)Δt=N−HN, _H:=(N-H) t= N-HN, and let EH:=∑k=HN−1Δt‖uk⋆‖2E_H:= _k=H^N-1 t\,\|u_k \|^2 be the control energy of its omitted tail. Construct a feasible candidate for the H-step rollout problem by keeping the first H optimal controls and setting all subsequent controls to zero: u~k=uk⋆,k=0,…,H−1,0,k=H,…,N−1. u_k= casesu_k ,&k=0,…,H-1,\\ 0,&k=H,…,N-1. cases Let x~k x_k denote the resulting trajectory. Since the two trajectories coincide at time H, define δk:=xk⋆−x~k,δH=0. _k:=x_k - x_k, _H=0. For k≥Hk≥ H, affine dynamics give δk+1=(I+ΔtAk)δk+Δtuk⋆. _k+1=(I+ tA_k) _k+ t\,u_k . Writing Φℓ,k:=∏s=kℓ−1(I+ΔtAs),Φk,k:=I, _ ,k:= _s=k -1(I+ tA_s), _k,k:=I, we obtain δN=∑k=HN−1ΔtΦN,k+1uk⋆. _N= _k=H^N-1 t\, _N,k+1u_k . By the definition of κA _A, ‖δN‖≤κA∑k=HN−1Δt‖uk⋆‖.\| _N\|≤ _A _k=H^N-1 t\,\|u_k \|. Cauchy–Schwarz then gives ‖δN‖≤κAτHEH.\| _N\|≤ _A _HE_H. Because u~ u is feasible for the rollout problem, J^0,H⋆(x0)−J0⋆(x0)≤J⁡(~)−J⁡(⋆). J_0,H (x_0)-J_0 (x_0)≤ J( x)-J( x ). The first H control costs cancel, while the omitted tail energy appears with a negative sign: J⁡(~)−J⁡(⋆)=−EH2λ+r⁡(xN⋆)−r⁡(x~N).J( x)-J( x )=- E_H2λ+r(x_N )-r( x_N). Using the LrL_r-Lipschitz continuity of r, J^0,H⋆(x0)−J0⋆(x0)≤LrκAτHEH−EH2λ. J_0,H (x_0)-J_0 (x_0)≤ L_r _A _HE_H- E_H2λ. Completing the square in EH E_H yields LrκAτHEH−EH2λ≤λ2κA2Lr2τH.L_r _A _HE_H- E_H2λ≤ λ2 _A^2L_r^2 _H. Therefore, 0≤J^0,H⋆(x0)−J0⋆(x0)≤λ2κA2Lr2N−HN.0≤ J_0,H (x_0)-J_0 (x_0)≤ λ2 _A^2L_r^2 N-HN. Combining this estimate with Step 1 and using the global optimality of ⋆ x gives 0≤J⁡(^H,⋆)−J⁡(⋆)≤λ2κA2Lr2N−HN.0≤ J( x^H, )-J( x )≤ λ2 _A^2L_r^2 N-HN. Thus, CTU=λ2κA2Lr2,C_T U= λ2 _A^2L_r^2, which proves the claim. ∎ B.3 Proof of Lemma 2 Proof. We follow the argument of Lemma 5, with one additional step: truncating a feasible long-horizon trajectory need not leave a feasible endpoint. The uniform Slater condition allows us to restore feasibility by convexly mixing the truncated trajectory with a strictly feasible one. For j+L≤Nj+L≤ N, let Jj,LC,⋆(x):=minuj:j+L−1∑k=j+L−1Δt2λ∥uk∥2−r(xj+L)J_j,L^C, (x):= _u_j:j+L-1 \ _k=j^j+L-1 t2λ\|u_k\|^2-r(x_j+L) \ subject to the affine dynamics, xj=x_j=x, and the terminal constraint h⁡(xj+L)≤0h(x_j+L)≤ 0. We write JjC,⋆(x):=Jj,N−jC,⋆(x).J_j^C, (x):=J_j,N-j^C, (x). Recall the constants χh:=Lh(J¯−J¯)γ,Gh:=Lr+χh, _h:= L_h( J- J)γ, G_h:=L_r+ _h, where γ is the uniform Slater margin, J¯ J uniformly bounds the objective of the strictly feasible Slater trajectories, and J¯ J is a uniform lower bound on the relevant finite-horizon objectives. The latter exists by the quadratic control penalty and the Lipschitz terminal reward, as in the coercivity argument used previously. Step 1: Local constrained horizon sensitivity. Fix j,L,qj,L,q with j+L+q≤Nj+L+q≤ N, and set τ:=qΔt.τ:=q t. We show that |Jj,L+qC,⋆(x)−Jj,LC,⋆(x)|≤ClocCτ, |J_j,L+q^C, (x)-J_j,L^C, (x) |≤ C_loc Cτ, where ClocC:=max⁡GhM+λGh22,M22λ.C_loc C:= \G_hM+ λ G_h^22, M^22λ \. Let uj:j+L+q−1⋆u_j:j+L+q-1 be optimal for Jj,L+qC,⋆(x)J_j,L+q^C, (x), with trajectory xk⋆x_k , and let E:=∑k=j+Lj+L+q−1Δt‖uk⋆‖2E:= _k=j+L^j+L+q-1 t\,\|u_k \|^2 denote the energy in the omitted tail. Define d:=Mτ+τE.d:=Mτ+ τ E. As in the proof of Lemma 5, ‖xj+L+q⋆−xj+L⋆‖≤d.\|x_j+L+q -x_j+L \|≤ d. Although xj+L+q⋆x_j+L+q is feasible, the truncated endpoint xj+L⋆x_j+L need not be. Since h⁡(xj+L+q⋆)≤0h(x_j+L+q )≤ 0 and h is LhL_h-Lipschitz, each component satisfies hi(xj+L⋆)≤Lhd.h_i(x_j+L )≤ L_hd. Set δ:=Lhd.δ:=L_hd. By Assumption 3, the L-step problem admits a strictly feasible Slater trajectory (xks,uks)(x_k s,u_k s) from the same initial state x, satisfying h⁡(xj+Ls)≤−γ,Jj,L(xs,us)≤J¯.h(x_j+L s)≤- 1, J_j,L(x s,u s)≤ J. Choose α:=δγ+δ.α:= δγ+δ. In particular, 0≤α≤δγ.0≤α≤ δγ. Consider the convexly mixed controls u~k:=(1−α)uk⋆+αuks,k=j,…,j+L−1. u_k:=(1-α)u_k +α u_k s, k=j,…,j+L-1. Because the dynamics are affine and both trajectories start from the same state, their corresponding states satisfy x~k=(1−α)xk⋆+αxks. x_k=(1-α)x_k +α x_k s. By componentwise convexity of h, hi(x~j+L) h_i( x_j+L) ≤(1−α)hi(xj+L⋆)+αhi(xj+Ls) ≤(1-α)h_i(x_j+L )+α h_i(x_j+L s) ≤(1−α)δ−αγ=0. ≤(1-α)δ-αγ=0. Thus the mixed trajectory is feasible for the L-step problem. Since the dynamics are affine, the quadratic running cost is convex in the controls, and −r-r is convex because r is concave. Hence the truncated objective is convex, so Jj,L(x~,u~)≤(1−α)Jj,L(x⋆,u⋆)+αJj,L(xs,us).J_j,L( x, u)≤(1-α)J_j,L(x ,u )+α J_j,L(x s,u s). Using Jj,L(x⋆,u⋆)≥J¯J_j,L(x ,u )≥ J and Jj,L(xs,us)≤J¯J_j,L(x s,u s)≤ J gives Jj,L(x~,u~)≤Jj,L(x⋆,u⋆)+α⁡(J¯−J¯).J_j,L( x, u)≤ J_j,L(x ,u )+α( J- J). Therefore, Jj,LC,⋆(x)≤Jj,L(x⋆,u⋆)+α⁡(J¯−J¯).J_j,L^C, (x)≤ J_j,L(x ,u )+α( J- J). The objective of the truncated prefix and that of the original (L+q)(L+q)-step optimum differ by Jj,L(x⋆,u⋆)−Jj,L+qC,⋆(x) J_j,L(x ,u )-J_j,L+q^C, (x) =r⁡(xj+L+q⋆)−r⁡(xj+L⋆)−E2λ =r(x_j+L+q )-r(x_j+L )- E2λ ≤Lrd−E2λ. ≤ L_rd- E2λ. Moreover, α⁡(J¯−J¯)≤Lh(J¯−J¯)γd=χhd.α( J- J)≤ L_h( J- J)γd= _hd. Combining the preceding inequalities, Jj,LC,⋆(x)−Jj,L+qC,⋆(x)≤Ghd−E2λ.J_j,L^C, (x)-J_j,L+q^C, (x)≤ G_hd- E2λ. Substituting d=Mτ+τEd=Mτ+ τ E, Jj,LC,⋆(x)−Jj,L+qC,⋆(x)≤GhMτ+GhτE−E2λ.J_j,L^C, (x)-J_j,L+q^C, (x)≤ G_hMτ+G_h τ E- E2λ. Completing the square gives GhτE−E2λ≤λGh22τ,G_h τ E- E2λ≤ λ G_h^22τ, and hence Jj,LC,⋆(x)−Jj,L+qC,⋆(x)≤(GhM+λGh22)τ.J_j,L^C, (x)-J_j,L+q^C, (x)≤ (G_hM+ λ G_h^22 )τ. For the reverse direction, let an optimal trajectory for Jj,LC,⋆(x)J_j,L^C, (x) terminate at x¯ x. Since h⁡(x¯)≤0h( x)≤ 0, extend it for the next q steps using uk=−(Akx¯+ck),k=j+L,…,j+L+q−1.u_k=- (A_k x+c_k ), k=j+L,…,j+L+q-1. This keeps the state fixed at x¯ x, so feasibility and the terminal reward are preserved. By the drift bound, ‖uk‖≤M,\|u_k\|≤ M, and therefore the extension costs at most M22λτ. M^22λτ. Thus Jj,L+qC,⋆(x)−Jj,LC,⋆(x)≤M22λτ.J_j,L+q^C, (x)-J_j,L^C, (x)≤ M^22λτ. Combining the two directions, |Jj,L+qC,⋆(x)−Jj,LC,⋆(x)|≤ClocCqΔt. |J_j,L+q^C, (x)-J_j,L^C, (x) |≤ C_loc Cq t. Step 2: Accumulation along the constrained RH-H trajectory. The remainder follows the telescoping argument in Lemma 5. We include the main steps for completeness. Let m:=N−H,WjC(x):=Jj,HC,⋆(x),j=0,…,m.m:=N-H, W_j C(x):=J_j,H^C, (x), j=0,…,m. At j=mj=m, the horizon reaches the true terminal time, so WmC(x)=JmC,⋆(x).W_m C(x)=J_m^C, (x). Let xjH,⋆x_j^H, and ujH,⋆u_j^H, denote the constrained RH-H closed-loop state and applied control. Bellman’s principle for the H-step constrained problem gives WjC(xjH,⋆)=Δt2λ‖ujH,⋆‖2+Jj+1,H−1C,⋆(xj+1H,⋆).W_j C(x_j^H, )= t2λ\|u_j^H, \|^2+J_j+1,H-1^C, (x_j+1^H, ). Applying the local horizon estimate with q=1q=1, Wj+1C(xj+1H,⋆)−Jj+1,H−1C,⋆(xj+1H,⋆)≤ClocCΔt.W_j+1 C(x_j+1^H, )-J_j+1,H-1^C, (x_j+1^H, )≤ C_loc C t. Consequently, Δt2λ‖ujH,⋆‖2≤WjC(xjH,⋆)−Wj+1C(xj+1H,⋆)+ClocCΔt. t2λ\|u_j^H, \|^2≤ W_j C(x_j^H, )-W_j+1 C(x_j+1^H, )+C_loc C t. Summing from j=0j=0 to m−1m-1 yields ∑j=0m−1Δt2λ‖ujH,⋆‖2≤W0C(x0)−WmC(xmH,⋆)+mClocCΔt. _j=0^m-1 t2λ\|u_j^H, \|^2≤ W_0 C(x_0)-W_m C(x_m^H, )+mC_loc C t. From time m=N−Hm=N-H onward, the remaining horizon is at most H, so RH-H coincides with full remaining-horizon re-optimization. Its remaining cost is therefore JmC,⋆(xmH,⋆)=WmC(xmH,⋆).J_m^C, (x_m^H, )=W_m C(x_m^H, ). Hence J⁡(H,⋆)≤W0C(x0)+mClocCΔt.J( x^H, )≤ W_0 C(x_0)+mC_loc C t. Finally, W0C(x0)=J0,HC,⋆(x0),J⁡(⋆)=J0,NC,⋆(x0).W_0 C(x_0)=J_0,H^C, (x_0), J( x )=J_0,N^C, (x_0). Applying the local constrained horizon estimate once more with L=HL=H and q=mq=m gives W0C(x0)−J⁡(⋆)≤ClocCmΔt.W_0 C(x_0)-J( x )≤ C_loc Cm t. Therefore, J⁡(H,⋆)−J⁡(⋆)≤2ClocCmΔt.J( x^H, )-J( x )≤ 2C_loc Cm t. Since ⋆ x is globally optimal among feasible trajectories, J⁡(H,⋆)−J⁡(⋆)≥0.J( x^H, )-J( x )≥ 0. Using m=N−Hm=N-H and Δt=1/N t=1/N, 0≤J⁡(H,⋆)−J⁡(⋆)≤CHCN−HN,0≤ J( x^H, )-J( x )≤ C_H C N-HN, where CHC:=2ClocC=2max⁡GhM+λGh22,M22λ.C_H C:=2C_loc C=2 \G_hM+ λ G_h^22, M^22λ \. This proves the claim. ∎ B.4 Proof of Lemma 3 Proof. The proof combines the no-accumulation argument of Lemma 6 with the Slater-mixing argument used in Lemma 2. The only new issue is that truncating the control horizon and rolling forward with zero control may violate the terminal constraint; the uniform Slater condition is used to repair this candidate without changing the order of the bound. For j+L≤Nj+L≤ N, define the constrained rollout value J^j,LC,⋆(x):=minuj:j+L−1∑k=j+L−1Δt2λ∥uk∥2−r(x^N), J_j,L^C, (x):= _u_j:j+L-1 \ _k=j^j+L-1 t2λ\|u_k\|^2-r( x_N) \, where the state follows the controlled affine dynamics for k=j,…,j+L−1k=j,…,j+L-1, is subsequently propagated with uk=0u_k=0 until time N, and satisfies h⁡(x^N)≤0.h( x_N)≤ 0. As before, JjC,⋆(x)J_j^C, (x) denotes the constrained full remaining-horizon value. Recall χh:=Lh(J¯−J¯)γ,Gh:=Lr+χh, _h:= L_h( J- J)γ, G_h:=L_r+ _h, where J¯ J is a uniform lower bound on the relevant rollout objectives and J¯ J bounds the objective of the corresponding strictly feasible Slater trajectories. Step 1: Terminal rollout does not accumulate horizon errors. The nesting argument from Lemma 6 remains valid in the constrained case. Indeed, any feasible L-step rollout control remains feasible for the (L+1)(L+1)-step problem by setting the additional control to zero; the resulting terminal state is unchanged. Hence JjC,⋆(x)≤J^j,L+1C,⋆(x)≤J^j,LC,⋆(x).J_j^C, (x)≤ J_j,L+1^C, (x)≤ J_j,L^C, (x). Let m:=N−H,W^jC(x):=J^j,HC,⋆(x),j=0,…,m.m:=N-H, W_j C(x):= J_j,H^C, (x), j=0,…,m. At j=mj=m, the control horizon reaches the true terminal time, so W^mC(x)=JmC,⋆(x). W_m C(x)=J_m^C, (x). Let x^jH,⋆ x_j^H, and u^jH,⋆ u_j^H, denote the constrained RH-T closed-loop state and applied control. Optimality gives W^jC(x^jH,⋆)=Δt2λ‖u^jH,⋆‖2+J^j+1,H−1C,⋆(x^j+1H,⋆). W_j C( x_j^H, )= t2λ\| u_j^H, \|^2+ J_j+1,H-1^C, ( x_j+1^H, ). By nestedness, W^j+1C(x^j+1H,⋆)≤J^j+1,H−1C,⋆(x^j+1H,⋆), W_j+1 C( x_j+1^H, )≤ J_j+1,H-1^C, ( x_j+1^H, ), and therefore Δt2λ‖u^jH,⋆‖2≤W^jC(x^jH,⋆)−W^j+1C(x^j+1H,⋆). t2λ\| u_j^H, \|^2≤ W_j C( x_j^H, )- W_j+1 C( x_j+1^H, ). Summing over j=0,…,m−1j=0,…,m-1 and using W^mC=JmC,⋆ W_m C=J_m^C, yields J⁡(^H,⋆)≤J^0,HC,⋆(x0).J( x^H, )≤ J_0,H^C, (x_0). Consequently, 0≤J⁡(^H,⋆)−J⁡(⋆)≤J^0,HC,⋆(x0)−J0C,⋆(x0).0≤ J( x^H, )-J( x )≤ J_0,H^C, (x_0)-J_0^C, (x_0). It therefore remains only to bound this initial constrained rollout gap. Step 2: Bounding the constrained rollout gap. Let u0⋆,…,uN−1⋆u_0 ,…,u_N-1 be the controls associated with the full constrained optimum ⋆ x , and define τH:=(N−H)Δt=N−HN, _H:=(N-H) t= N-HN, together with its omitted-tail energy EH:=∑k=HN−1Δt‖uk⋆‖2.E_H:= _k=H^N-1 t\,\|u_k \|^2. As in Lemma 6, retain the first H optimal controls and set the remaining controls to zero: u~k=uk⋆,k=0,…,H−1,0,k=H,…,N−1. u_k= casesu_k ,&k=0,…,H-1,\\ 0,&k=H,…,N-1. cases Let x~N x_N denote the resulting terminal state. The affine perturbation argument of Lemma 6 gives ‖xN⋆−x~N‖≤κAτHEH.\|x_N - x_N\|≤ _A _HE_H. Since the full-horizon optimum is feasible, h⁡(xN⋆)≤0h(x_N )≤ 0, and h is LhL_h-Lipschitz. Thus each terminal constraint violation of the truncated rollout is bounded by hi(x~N)≤δH,δH:=LhκAτHEH.h_i( x_N)≤ _H, _H:=L_h _A _HE_H. The zero-tail rollout need not itself be feasible. By Assumption 3, however, the H-step rollout problem admits a strictly feasible trajectory (xs,us)(x s,u s) from the same initial state satisfying h⁡(xNs)≤−γ,J^0,H(xs,us)≤J¯.h(x_N s)≤- 1, J_0,H(x s,u s)≤ J. Set α:=δHγ+δH≤δHγ,α:= _Hγ+ _H≤ _Hγ, and form the mixed rollout controls ukα:=(1−α)u~k+αuks,k=0,…,H−1,u_k^α:=(1-α) u_k+α u_k s, k=0,…,H-1, with zero control thereafter. Because the dynamics are affine, the resulting terminal state satisfies xNα=(1−α)x~N+αxNs.x_N^α=(1-α) x_N+α x_N s. By componentwise convexity of h, hi(xNα) h_i(x_N^α) ≤(1−α)hi(x~N)+αhi(xNs) ≤(1-α)h_i( x_N)+α h_i(x_N s) ≤(1−α)δH−αγ=0. ≤(1-α) _H-αγ=0. Hence the mixed rollout is feasible. The rollout objective is convex in the active controls, since the dynamics are affine, the running cost is quadratic, and −r-r is convex. Therefore, J^0,H(uα)≤(1−α)J^0,H(u~)+αJ^0,H(us). J_0,H(u^α)≤(1-α) J_0,H( u)+α J_0,H(u s). Using the uniform lower and upper bounds J^0,H(u~)≥J¯ J_0,H( u)≥ J and J^0,H(us)≤J¯ J_0,H(u s)≤ J, J^0,H(uα)≤J^0,H(u~)+α⁡(J¯−J¯). J_0,H(u^α)≤ J_0,H( u)+α( J- J). Moreover, α⁡(J¯−J¯)≤Lh(J¯−J¯)γκAτHEH=χhκAτHEH.α( J- J)≤ L_h( J- J)γ _A _HE_H= _h _A _HE_H. Since uαu^α is feasible for the constrained rollout problem, J^0,HC,⋆(x0)≤J^0,H(uα). J_0,H^C, (x_0)≤ J_0,H(u^α). On the other hand, comparing the zero-tail trajectory with the full optimum gives J^0,H(u~)−J0C,⋆(x0)=−EH2λ+r⁡(xN⋆)−r⁡(x~N). J_0,H( u)-J_0^C, (x_0)=- E_H2λ+r(x_N )-r( x_N). Using the LrL_r-Lipschitz continuity of r, J^0,H(u~)−J0C,⋆(x0)≤LrκAτHEH−EH2λ. J_0,H( u)-J_0^C, (x_0)≤ L_r _A _HE_H- E_H2λ. Combining the reward and feasibility-repair contributions, J^0,HC,⋆(x0)−J0C,⋆(x0)≤GhκAτHEH−EH2λ. J_0,H^C, (x_0)-J_0^C, (x_0)≤ G_h _A _HE_H- E_H2λ. Completing the square in EH E_H gives GhκAτHEH−EH2λ≤λ2κA2Gh2τH.G_h _A _HE_H- E_H2λ≤ λ2 _A^2G_h^2 _H. Therefore, 0≤J^0,HC,⋆(x0)−J0C,⋆(x0)≤λ2κA2Gh2N−HN.0≤ J_0,H^C, (x_0)-J_0^C, (x_0)≤ λ2 _A^2G_h^2 N-HN. Combining this estimate with Step 1 yields 0≤J⁡(^H,⋆)−J⁡(⋆)≤λ2κA2Gh2N−HN.0≤ J( x^H, )-J( x )≤ λ2 _A^2G_h^2 N-HN. Thus CTC=λ2κA2Gh2,C_T C= λ2 _A^2G_h^2, which proves the claim. ∎ B.5 Proof of Lemma 4 We first collect two properties of each frozen RH problem that are used in the proof. Fix a physical index j and current state x¯j x_j. Let Hj=min⁡H,N−jH_j= \H,N-j\ and write j=(xj+1,…,xj+Hj)x_j=(x_j+1,…,x_j+H_j). Under Assumption 2, the stacked residuals satisfy col⁡(ej,…,ej+Hj−1)=Bjj−aj−Sjx¯j,col(e_j,…,e_j+H_j-1)=B_jx_j-a_j-S_j x_j, where Bj=[Id0⋯0−(Id+ΔtAj+1)Id⋱0⋱00⋯−(Id+ΔtAj+Hj−1)Id],Sj=[Id+ΔtAj00],B_j= bmatrixI_d&0&·s&0\\ -(I_d+ tA_j+1)&I_d& & \\ 0& & &0\\ 0&·s&-(I_d+ tA_j+H_j-1)&I_d bmatrix, S_j= bmatrixI_d+ tA_j\\ 0\\ \\ 0 bmatrix, and aj=Δtcol⁡(cj,…,cj+Hj−1)a_j= t\,col(c_j,…,c_j+H_j-1). Recall σj:=σmin2(Bj)λΔt>0. _j:= _ ^2(B_j)λ t>0. For RH-T, let QjQ_j denote the linear part of the terminal map Φj _j; for RH-H, Qj=IQ_j=I. Thus Lg,j:=Lh‖Qj‖,Ld,j:=Lg,j2σj.L_g,j:=L_h\|Q_j\|, L_d,j:= L_g,j^2 _j. Lemma 7 (Strong convexity and well-posedness). Under Assumptions 2–3, φj(⋅,x¯j) _j(·; x_j) and ℒj(⋅,μ,x¯j)L_j(·,μ; x_j) are σj _j-strongly convex for every μ≥0μ≥ 0. Moreover, the primal response j(μ):=argminjℒj(j,μ,x¯j)x_j(μ):= _x_jL_j(x_j,μ; x_j) exists and is unique, and the constrained local problem has a unique optimizer j⋆x_j . Proof. Since BjB_j is block lower triangular with identity diagonal blocks, it is nonsingular. The Hessian of the quadratic path cost is 1λΔtBj⊤Bj⪰σjI. 1λ tB_j B_j _jI. Moreover, Φj _j is affine, r is concave, and each component of h is convex. Hence −r∘Φj-r _j and gj,i=hi∘Φjg_j,i=h_i _j are convex. Therefore φj _j is σj _j-strongly convex, and so is ℒj=φj+μ⊤gjL_j= _j+μ g_j for every μ≥0μ≥ 0. It remains to show attainment. Let X=jX=x_j, ζj=aj+Sjx¯j _j=a_j+S_j x_j, and write qj(X):=12λΔt‖BjX−ζj‖2,ψμ(X):=−r⁡(Φj(xj+Hj))+μ⊤gj(X).q_j(X):= 12λ t\|B_jX- _j\|^2, _μ(X):=-r( _j(x_j+H_j))+μ g_j(X). The finite-valued convex function ψμ _μ admits an affine lower bound: for any fixed X X and sμ∈∂ψμ(X^)s_μ∈∂ _μ( X), ψμ(X)≥ψμ(X^)+sμ⊤(X−X^). _μ(X)≥ _μ( X)+s_μ (X- X). Also, qj(X) q_j(X) ≥14λΔt‖BjX‖2−12λΔt‖ζj‖2 ≥ 14λ t\|B_jX\|^2- 12λ t\| _j\|^2 ≥σmin2(Bj)4λΔt‖X‖2−12λΔt‖ζj‖2. ≥ _ ^2(B_j)4λ t\|X\|^2- 12λ t\| _j\|^2. Thus ℒj(X,μ,x¯j)≥σmin2(Bj)4λΔt‖X‖2−‖sμ‖‖X‖−CL_j(X,μ; x_j)≥ _ ^2(B_j)4λ t\|X\|^2-\|s_μ\|\,\|X\|-C for a finite constant C, and hence the Lagrangian is coercive. Therefore its minimum is attained; strong convexity makes the minimizer unique. Finally, Assumption 3 gives a nonempty convex feasible set for the frozen RH problem. Since its objective is strongly convex and coercive, the constrained minimizer exists and is unique. ∎ Lemma 8 (Dual regularity). Under Assumptions 2–3, strong duality holds, ℳj⋆M_j is nonempty and compact, and djd_j is differentiable with ∇dj(μ)=gj(j(μ)).∇ d_j(μ)=g_j(x_j(μ)). Furthermore, ‖j(μ)−j(ν)‖≤Lg,jσj‖μ−ν‖,\|x_j(μ)-x_j(ν)\|≤ L_g,j _j\|μ-ν\|, and ‖∇dj(μ)−∇dj(ν)‖≤Ld,j‖μ−ν‖.\|∇ d_j(μ)-∇ d_j(ν)\|≤ L_d,j\|μ-ν\|. Proof. By Assumption 3, the local convex problem satisfies Slater’s condition. Standard convex duality therefore gives strong duality and existence of an optimal multiplier. To see boundedness, let ~j x_j be a Slater point with gj(~j)≤−γg_j( x_j)≤- 1. For any μ⋆∈ℳj⋆μ _j , pj⋆=dj(μ⋆) p_j =d_j(μ ) ≤ℒj(~j,μ⋆,x¯j) _j( x_j,μ ; x_j) ≤φj(~j,x¯j)−γ‖μ⋆‖1. ≤ _j( x_j; x_j)-γ\|μ \|_1. Hence ‖μ⋆‖1≤φj(~j,x¯j)−pj⋆γ.\|μ \|_1≤ _j( x_j; x_j)-p_j γ. Thus ℳj⋆M_j is bounded; closedness follows from continuity of djd_j, proved below, so it is compact. We next establish Lipschitzness of the primal response. Strong convexity and optimality of j(μ)x_j(μ) and j(ν)x_j(ν) give ℒj(j(ν),μ)≥ℒj(j(μ),μ)+σj2‖j(ν)−j(μ)‖2L_j(x_j(ν),μ) _j(x_j(μ),μ)+ _j2\|x_j(ν)-x_j(μ)\|^2 and ℒj(j(μ),ν)≥ℒj(j(ν),ν)+σj2‖j(μ)−j(ν)‖2.L_j(x_j(μ),ν) _j(x_j(ν),ν)+ _j2\|x_j(μ)-x_j(ν)\|^2. Adding these inequalities and cancelling the common objective terms yields σj‖j(μ)−j(ν)‖2≤(μ−ν)⊤(gj(j(ν))−gj(j(μ))). _j\|x_j(μ)-x_j(ν)\|^2≤(μ-ν) (g_j(x_j(ν))-g_j(x_j(μ)) ). Using Cauchy–Schwarz and the Lg,jL_g,j-Lipschitzness of gjg_j, σj‖j(μ)−j(ν)‖2≤Lg,j‖μ−ν‖‖j(μ)−j(ν)‖, _j\|x_j(μ)-x_j(ν)\|^2≤ L_g,j\|μ-ν\|\|x_j(μ)-x_j(ν)\|, which proves ‖j(μ)−j(ν)‖≤Lg,jσj‖μ−ν‖.\|x_j(μ)-x_j(ν)\|≤ L_g,j _j\|μ-ν\|. To identify the dual gradient, optimality of the two primal responses gives dj(ν)≤dj(μ)+(ν−μ)⊤gj(j(μ))d_j(ν)≤ d_j(μ)+(ν-μ) g_j(x_j(μ)) and dj(μ)≤dj(ν)+(μ−ν)⊤gj(j(ν)).d_j(μ)≤ d_j(ν)+(μ-ν) g_j(x_j(ν)). Therefore |dj(ν)−dj(μ)−(ν−μ)⊤gj(j(μ))| |d_j(ν)-d_j(μ)-(ν-μ) g_j(x_j(μ)) | ≤‖ν−μ‖‖gj(j(ν))−gj(j(μ))‖ ≤\|ν-μ\|\,\|g_j(x_j(ν))-g_j(x_j(μ))\| ≤Lg,j2σj‖ν−μ‖2. ≤ L_g,j^2 _j\|ν-μ\|^2. Hence ∇dj(μ)=gj(j(μ))∇ d_j(μ)=g_j(x_j(μ)). Finally, ‖∇dj(μ)−∇dj(ν)‖ \|∇ d_j(μ)-∇ d_j(ν)\| ≤Lg,j‖j(μ)−j(ν)‖ ≤ L_g,j\|x_j(μ)-x_j(ν)\| ≤Lg,j2σj‖μ−ν‖=Ld,j‖μ−ν‖. ≤ L_g,j^2 _j\|μ-ν\|=L_d,j\|μ-ν\|. ∎ Proof of Lemma 4. Fix j and suppress the dependence on x¯j x_j. From Lemma 8, ⟨∇dj(μ)−∇dj(ν),μ−ν⟩ ∇ d_j(μ)-∇ d_j(ν),μ-ν =(μ−ν)⊤(gj(j(μ))−gj(j(ν))) =(μ-ν) (g_j(x_j(μ))-g_j(x_j(ν)) ) ≤−σj‖j(μ)−j(ν)‖2 ≤- _j\|x_j(μ)-x_j(ν)\|^2 ≤−1Ld,j‖∇dj(μ)−∇dj(ν)‖2. ≤- 1L_d,j\|∇ d_j(μ)-∇ d_j(ν)\|^2. Let T(μ):=Projℝ+m(μ+ημ∇dj(μ)).T(μ):=Proj_R_+^m (μ+ _μ∇ d_j(μ) ). A point z maximizes djd_j over ℝ+mR_+^m if and only if T⁡(z)=zT(z)=z. For z∈ℳj⋆z _j , nonexpansiveness of the projection and the preceding inequality imply ‖μj(s+1)−z‖2 \| _j^(s+1)-z\|^2 ≤‖μj(s)−z+ημ(∇dj(μj(s))−∇dj(z))‖2 ≤\| _j^(s)-z+ _μ(∇ d_j( _j^(s))-∇ d_j(z))\|^2 ≤‖μj(s)−z‖2−ημ(2Ld,j−ημ)‖∇dj(μj(s))−∇dj(z)‖2. ≤\| _j^(s)-z\|^2- _μ ( 2L_d,j- _μ )\|∇ d_j( _j^(s))-∇ d_j(z)\|^2. Thus, whenever 0<ημ<2/Ld,j0< _μ<2/L_d,j, the sequence is Fejér monotone with respect to ℳj⋆M_j and therefore bounded. Using firm nonexpansiveness of the projection gives more precisely ‖μj(s+1)−z‖2+‖μj(s)−μj(s+1)+ημ(∇dj(μj(s))−∇dj(z))‖2 \| _j^(s+1)-z\|^2+\| _j^(s)- _j^(s+1)+ _μ(∇ d_j( _j^(s))-∇ d_j(z))\|^2 ≤‖μj(s)−z‖2−ημ(2Ld,j−ημ)‖∇dj(μj(s))−∇dj(z)‖2. ≤\| _j^(s)-z\|^2- _μ ( 2L_d,j- _μ )\|∇ d_j( _j^(s))-∇ d_j(z)\|^2. Summing over s shows both squared terms on the left are summable. Consequently, ‖μj(s+1)−μj(s)‖→0.\| _j^(s+1)- _j^(s)\|→ 0. Since the sequence is bounded, it admits a convergent subsequence μj(sk)→μ¯j _j^(s_k)→ μ_j. Passing to the limit in μj(s+1)=T⁡(μj(s)) _j^(s+1)=T( _j^(s)) gives T⁡(μ¯j)=μ¯jT( μ_j)= μ_j; hence μ¯j∈ℳj⋆ μ_j _j . Fejér monotonicity with respect to μ¯j μ_j then implies convergence of the full sequence: μj(s)→μ¯j. _j^(s)→ μ_j. Lemma 8 gives continuity of the primal response, and Lemma 7 gives uniqueness. Therefore j(s)=j(μj(s))⟶j(μ¯j)=j⋆.x_j^(s)=x_j( _j^(s)) _j( μ_j)=x_j . Finite-iteration bounds. Assume now 0<ημ≤1/Ld,j0< _μ≤ 1/L_d,j. The projection optimality condition gives, for any z≥0z≥ 0, ⟨∇dj(μj(k)),z−μj(k+1)⟩≤‖z−μj(k)‖2−‖z−μj(k+1)‖2−‖μj(k+1)−μj(k)‖22ημ. ∇ d_j( _j^(k)),z- _j^(k+1) ≤ \|z- _j^(k)\|^2-\|z- _j^(k+1)\|^2-\| _j^(k+1)- _j^(k)\|^22 _μ. By concavity and Ld,jL_d,j-smoothness, dj(z)−dj(μj(k+1)) d_j(z)-d_j( _j^(k+1)) ≤⟨∇dj(μj(k)),z−μj(k+1)⟩ ≤ ∇ d_j( _j^(k)),z- _j^(k+1) +Ld,j2‖μj(k+1)−μj(k)‖2. + L_d,j2\| _j^(k+1)- _j^(k)\|^2. Hence dj(z)−dj(μj(k+1)) d_j(z)-d_j( _j^(k+1)) ≤‖z−μj(k)‖2−‖z−μj(k+1)‖22ημ ≤ \|z- _j^(k)\|^2-\|z- _j^(k+1)\|^22 _μ −(12ημ−Ld,j2)‖μj(k+1)−μj(k)‖2. - ( 12 _μ- L_d,j2 )\| _j^(k+1)- _j^(k)\|^2. Taking z=μj(k)z= _j^(k) shows that dj(μj(k))d_j( _j^(k)) is nondecreasing. Let μj†∈ℳj⋆ _j _j satisfy ‖μj(0)−μj†‖=Rj:=dist⁡(μj(0),ℳj⋆).\| _j^(0)- _j \|=R_j:=dist( _j^(0),M_j ). Summing the preceding inequality with z=μj†z= _j gives ∑k=0s−1(pj⋆−dj(μj(k+1)))≤Rj22ημ. _k=0^s-1 (p_j -d_j( _j^(k+1)) )≤ R_j^22 _μ. Since these gaps are nonincreasing, 0≤pj⋆−dj(μj(s))≤Rj22ημs.0≤ p_j -d_j( _j^(s))≤ R_j^22 _μs. Finally, strong convexity of ℒj(⋅,μ)L_j(·,μ) and feasibility of j⋆x_j imply σj2‖j(μ)−j⋆‖2 _j2\|x_j(μ)-x_j \|^2 ≤ℒj(j⋆,μ)−dj(μ) _j(x_j ,μ)-d_j(μ) ≤pj⋆−dj(μ). ≤ p_j -d_j(μ). Taking μ=μj(s)μ= _j^(s) yields ‖j(μj(s))−j⋆‖≤Rjσjημs.\|x_j( _j^(s))-x_j \|≤ R_j _j _μs. Since [⋅]+[\,·\,]_+ is nonexpansive and gj(j⋆)≤0g_j(x_j )≤ 0, ‖[gj(j(μj(s)))]+‖ \|[g_j(x_j( _j^(s)))]_+\| ≤‖gj(j(μj(s)))−gj(j⋆)‖ ≤\|g_j(x_j( _j^(s)))-g_j(x_j )\| ≤Lg,j‖j(μj(s))−j⋆‖ ≤ L_g,j\|x_j( _j^(s))-x_j \| ≤Lg,jRjσjημs. ≤ L_g,jR_j _j _μs. This proves the lemma. ∎ B.6 Proof of Theorem 1 We next propagate the local optimization error through the receding-horizon closed loop. For a current state ξ, let j⋆(ξ)x_j (ξ) denote the exact optimizer of the jjth RH problem and define πj(ξ):=[j⋆(ξ)]1. _j(ξ):=[x_j (ξ)]_1. Thus, the exact RH trajectory satisfies x0RH=x0x_0 RH=x_0 and xj+1RH=πj(xjRH)x_j+1 RH= _j(x_j RH). Lemma 9 (Lipschitz exact RH policy). Under Assumptions 2–3, ‖j⋆(ξ)−j⋆(ζ)‖≤κj‖ξ−ζ‖,κj:=‖Bj⊤Sj‖λΔtσj.\|x_j (ξ)-x_j (ζ)\|≤ _j\|ξ-ζ\|, _j:= \|B_j S_j\|λ t\, _j. Consequently, ‖πj(ξ)−πj(ζ)‖≤κj‖ξ−ζ‖.\| _j(ξ)- _j(ζ)\|≤ _j\|ξ-ζ\|. Proof. The constraint gj(j)=h⁡(Φj(xj+Hj))g_j(x_j)=h( _j(x_j+H_j)) does not depend explicitly on the current state. Hence the feasible set is the same for initial states ξ and ζ. Since φj(⋅,ξ) _j(·;ξ) is σj _j-strongly convex and j⋆(ξ)x_j (ξ) minimizes it over this common convex feasible set, φj(j⋆(ζ),ξ)≥φj(j⋆(ξ),ξ)+σj2‖j⋆(ξ)−j⋆(ζ)‖2. _j(x_j (ζ);ξ)≥ _j(x_j (ξ);ξ)+ _j2\|x_j (ξ)-x_j (ζ)\|^2. Exchanging ξ and ζ gives φj(j⋆(ξ),ζ)≥φj(j⋆(ζ),ζ)+σj2‖j⋆(ξ)−j⋆(ζ)‖2. _j(x_j (ξ);ζ)≥ _j(x_j (ζ);ζ)+ _j2\|x_j (ξ)-x_j (ζ)\|^2. Adding the two inequalities, σj‖j⋆(ξ)−j⋆(ζ)‖2 _j\|x_j (ξ)-x_j (ζ)\|^2 ≤φj(j⋆(ξ),ζ)−φj(j⋆(ξ),ξ) ≤ _j(x_j (ξ);ζ)- _j(x_j (ξ);ξ) +φj(j⋆(ζ),ξ)−φj(j⋆(ζ),ζ). + _j(x_j (ζ);ξ)- _j(x_j (ζ);ζ). The terminal-reward terms cancel, and expanding the two quadratic residual costs gives σj‖j⋆(ξ)−j⋆(ζ)‖2 _j\|x_j (ξ)-x_j (ζ)\|^2 ≤1λΔt⟨Bj⊤Sj(ξ−ζ),j⋆(ξ)−j⋆(ζ)⟩ ≤ 1λ t B_j S_j(ξ-ζ),x_j (ξ)-x_j (ζ) ≤‖Bj⊤Sj‖λΔt‖ξ−ζ‖‖j⋆(ξ)−j⋆(ζ)‖. ≤ \|B_j S_j\|λ t\|ξ-ζ\|\|x_j (ξ)-x_j (ζ)\|. Cancelling the final factor proves ‖j⋆(ξ)−j⋆(ζ)‖≤κj‖ξ−ζ‖.\|x_j (ξ)-x_j (ζ)\|≤ _j\|ξ-ζ\|. Since πj _j extracts the first block of the optimized path and this projection has operator norm one, the same bound holds for πj _j. ∎ Lemma 10 (Propagation of local solve errors). Suppose the implemented trajectory satisfies ∥xj+1−πj(xj)∥≤ϵj,j=0,…,N−1.\|x_j+1- _j(x_j)\|≤ _j, j=0,…,N-1. Then, for every p=1,…,Np=1,…,N, ‖xp−xpRH‖≤∑j=0p−1(∏ℓ=j+1p−1κℓ)ϵj,\|x_p-x_p RH\|≤ _j=0^p-1 ( _ =j+1^p-1 _ ) _j, where an empty product equals one. Moreover, if h⁡(πN−1(xN−1))≤0h( _N-1(x_N-1))≤ 0, then ‖[h⁡(xN)]+‖≤LhϵN−1.\|[h(x_N)]_+\|≤ L_h _N-1. Proof. Since x0=x0RHx_0=x_0 RH, the claim is immediate for p=1p=1: ‖x1−x1RH‖=‖x1−π0(x0)‖≤ϵ0.\|x_1-x_1 RH\|=\|x_1- _0(x_0)\|≤ _0. Assume the bound holds at time p. Then ‖xp+1−xp+1RH‖ \|x_p+1-x_p+1 RH\| ≤‖xp+1−πp(xp)‖+‖πp(xp)−πp(xpRH)‖ ≤\|x_p+1- _p(x_p)\|+\| _p(x_p)- _p(x_p RH)\| ≤ϵp+κp‖xp−xpRH‖ ≤ _p+ _p\|x_p-x_p RH\| ≤ϵp+∑j=0p−1(∏ℓ=j+1pκℓ)ϵj ≤ _p+ _j=0^p-1 ( _ =j+1^p _ ) _j =∑j=0p(∏ℓ=j+1pκℓ)ϵj, = _j=0^p ( _ =j+1^p _ ) _j, which proves the result by induction. For terminal feasibility, since h⁡(πN−1(xN−1))≤0h( _N-1(x_N-1))≤ 0 and h is LhL_h-Lipschitz, ‖[h⁡(xN)]+‖ \|[h(x_N)]_+\| =‖[h⁡(xN)]+−[h⁡(πN−1(xN−1))]+‖ =\|[h(x_N)]_+-[h( _N-1(x_N-1))]_+\| ≤‖h⁡(xN)−h⁡(πN−1(xN−1))‖ ≤\|h(x_N)-h( _N-1(x_N-1))\| ≤Lh‖xN−πN−1(xN−1)‖ ≤ L_h\|x_N- _N-1(x_N-1)\| ≤LhϵN−1. ≤ L_h _N-1. ∎ Proof of Theorem 1. At physical window j, Lemma 4 with s=Ks=K gives ‖j(K)−j⋆(xj)‖≤RjσjημK,\|x_j^(K)-x_j (x_j)\|≤ R_j _j _μK, where Rj=dist⁡(μj(0),ℳj⋆(xj)).R_j=dist ( _j^(0),M_j (x_j) ). Since the implemented next state is the first block of j(K)x_j^(K), while πj(xj) _j(x_j) is the first block of j⋆(xj)x_j (x_j), ‖xj+1−πj(xj)‖≤RjσjημK=:ϵj.\|x_j+1- _j(x_j)\|≤ R_j _j _μK=: _j. Define R¯:=maxj⁡Rj,σ¯:=minj⁡σj,κ¯:=maxj⁡κj. R:= _jR_j, σ:= _j _j, κ:= _j _j. Lemma 10 therefore gives ‖xp−xpRH‖ \|x_p-x_p RH\| ≤∑j=0p−1(∏ℓ=j+1p−1κℓ)RjσjημK ≤ _j=0^p-1 ( _ =j+1^p-1 _ ) R_j _j _μK ≤R¯σ¯ημK∑j=0p−1κ¯p−1−j ≤ R σ\, _μK _j=0^p-1 κ^\,p-1-j =R¯Γp(κ¯)σ¯ημK, = R\, _p( κ) σ\, _μK, where Γp(κ¯):=∑q=0p−1κ¯q. _p( κ):= _q=0^p-1 κ^q. Finally, at the last physical window j=N−1j=N-1, we have HN−1=1H_N-1=1. Hence πN−1(xN−1) _N-1(x_N-1) is the terminal state of the exact feasible local solution, so h⁡(πN−1(xN−1))≤0.h( _N-1(x_N-1))≤ 0. Applying Lemma 10, ‖[h⁡(xN)]+‖ \|[h(x_N)]_+\| ≤LhϵN−1 ≤ L_h _N-1 ≤LhR¯σ¯ημK. ≤ L_h R σ\, _μK. This proves the theorem. ∎ Appendix C Additional Experimental Details and Results C.1 Implementation Details C.1.1 Low-dimensional experiment implementation details Model and sampler. The pretrained model is a flow map on the two-moons density, parameterised as Xθ(s,t,x)=x+(t−s)σdϕθ(s,t,x/σd)X_θ(s,t,x)=x+(t-s)\, _d\, _θ\! (s,\,t,\,x/ _d ) with ϕθ _θ an MLP of three hidden layers of 256256 GELU units, input (s,t,x/σd)∈ℝ4(s,t,x/ _d) ^4 and output in ℝ2R^2. We use the EMA weights at 1515k training steps. The data are sklearn’s make_moons with noise 0.050.05, standardised per coordinate and then scaled by 0.750.75; σd=0.750 _d=0.750 is the standard deviation of that training set, and the base law is p0=⁡(0,σd2I)p_0=N(0, _d^2I). Unlike the image experiment, time runs forward, tj=j/Nt_j=j/N with N=16N=16 uniform Euler steps and Δt=1/16 t=1/16; Every run uses the same batch of 512512 particles drawn from the same seed, shared across all methods and all λ, so each frontier is computed on identical initial conditions. Rewards. Two terminal rewards are used, both bounded in (0,1](0,1] so that reward gains are directly comparable across methods. The upper-moon reward is an RBF envelope over one mode, r(x)=exp(−minj∥x−x(j)∥2/2σ2)r(x)= \! (- _j x-x^(j) ^2/2σ^2 ) with σ=0.35σ=0.35 and 400400 reference points x(j)x^(j) drawn from the upper moon; it is ≈1≈ 1 on that mode and decays off it. The uncontrolled sampler already attains r¯=0.589 r=0.589, so the attainable headroom is small. The off-manifold reward is a single Gaussian bump r(x)=exp(−∥x−p∥2/2σ2)r(x)= \! (- x-p ^2/2σ^2 ) with p=(−0.1, 1.75)p=(-0.1,\,1.75) and σ=0.5σ=0.5. Its optimum lies off the data manifold: the upper moon’s apex (−0.36, 1.41)(-0.36,\,1.41) scores only 0.690.69, and the uncontrolled sampler attains r¯=0.062 r=0.062. Reaching high reward therefore requires leaving the support of the learned density, which is what separates the methods. Control methods. All three methods share the sampler, the initial batch and the reward, and no constraints are active, so ConicFlow again reduces to K=1K=1 in Algorithm 1. Two departures from the image configuration are specific to this setting. First, the window subproblem is not warm-started across windows: each window is initialised from a fresh uncontrolled rollout of its own H steps, which isolates the effect of the horizon from that of solver warm-starting. Second, we report two choices of the terminal map Φj _j in (RH-T): endpoint, Φj=Xj+H,N _j=X_j+H,N, which predicts t=1t=1 from the window’s terminal state through one flow-map call (the image-experiment choice), and horizon, Φj=id _j=id, which evaluates the reward directly at the window’s terminal state xj+Hx_j+H and makes no network call beyond the rollout. HardFlow is run in its unconstrained form, so its terminal subproblem reduces to proximal reward steering, x^N⋆=arg⁡minx−r⁡(x)+αt22Δtλ∥x−x¯N∥2 x_N = _x-r(x)+ _t^22 t\,λ x- x_N ^2, solved by the augmented-Lagrangian routine with ρ0=1 _0=1 and ρ×2ρ× 2 per outer round, warm-started across sampling steps. Its released code parameterises the proximal weight as w=λocαt2/(2Δt)w= _oc _t^2/(2 t), i.e. with the opposite sense: larger λoc _oc holds the sample closer to the pretrained path. We therefore set λoc=1/λ _oc=1/λ so that a single λ axis means the same thing for all three methods — larger λ is a weaker control penalty and more reward — which is what makes the frontiers superimposable. FMRG is run in both variants, J (differentiating through Xj,NX_j,N) and E (stop-gradient at the map), with the step z←z+(Δt/ninner)λ∇rz← z+( t/n_inner)λ\,∇ r. λ grids and per-λ selection. Because the state space is two-dimensional, the inner hyperparameters can be swept exhaustively rather than range-found. For every (method, variant, λ) we run the full inner grid of Table 2 and report the configuration minimising the control objective itself, J=12λ∑j∥Δxj∥2−r¯(xN),J\;=\; 12λ _j x_j ^2\;-\; r(x_N), (2) where Δxj=xj+1−x¯j+1 x_j=x_j+1- x_j+1 is the applied deviation from the uncontrolled Euler step. Each plotted marker is that single best run, so the curves are attainment frontiers under a fixed selection rule rather than non-dominated envelopes. The horizontal axis of both figures is the mean integrated control effort ∑j∥Δxj∥2 _j x_j ^2 averaged over the batch, which equals ΔtE t\,E for the energy E of (3), i.e. E/16E/16 at N=16N=16. For the upper-moon reward all three methods share λ∈0.1, 0.25, 0.5, 0.75, 1, 2, 3λ∈\0.1,\,0.25,\,0.5,\,0.75,\,1,\,2,\,3\; beyond λ≈3λ≈ 3 every method has saturated, since the reward is already near 11 on the mode the sampler visits. For the off-manifold reward HardFlow and FMRG saturate by λ≈5λ≈ 5 and are swept over the 1313 values 0.1,…,100\0.1,…,100\, whereas ConicFlow is still improving at λ=100λ=100 and is extended geometrically to λ=105λ=10^5 (2222 values in total) until its reward plateaus. Table 2: Hyperparameters for the low-dimensional experiments. Shared across methods Steps N 1616 Batch / seed 512512 particles / 4242, shared across all runs Base law p0p_0 ⁡(0,σd2I)N(0, _d^2I), σd=0.750 _d=0.750 Backbone MLP, 3×2563× 256 GELU; EMA weights at 1515k steps Per method ConicFlow horizon H 1,2,4,8,16\1,2,4,8,16\ terminal map Φj _j endpoint Xj+H,NX_j+H,N or horizon idid inner optimizer Adam inner iterations 25,50\25,50\, extended to 100,200\100,200\ at large λ learning rate 0.01,0.05\0.01,0.05\, extended with 0.001,0.005,0.2\0.001,0.005,0.2\ at large λ window initialisation fresh uncontrolled rollout (no cross-window warm start) HardFlow inner optimizer L-BFGS (strong Wolfe) or Adam outer × inner 3,5,10×10,20\3,5,10\×\10,20\ learning rate 0.1,0.5\0.1,0.5\ (L-BFGS), 0.01,0.05\0.01,0.05\ (Adam) penalty ρ0 _0 / growth 11 / ×2× 2 per outer round proximal weight αt2/(2Δtλ) _t^2/(2 t\,λ), i.e. λoc=1/λ _oc=1/λ warm start across sampling steps FMRG variant J (through Xj,NX_j,N) and E (stop-gradient) guidance steps 1,3,5,10\1,3,5,10\ per sampling step C.1.2 Unconstrained image experiment implementation details Model and sampler. We use FLUX.1-dev with the FlowMap LoRA adapter, so one forward pass evaluates the flow map bt(z,t→t′)b_t(z,\,t\!→\!t ): sampling advances with t′=t−Δt =t- t, while t′=0t =0 estimates the terminal latent in a single call. The N steps are uniform, tj=1−j/Nt_j=1-j/N, so Δt=1/N t=1/N for the sampler. Each run draws one latent from a seed, and the uncontrolled sample from that same seed is the reference throughout. Reward. CLIP (Radford et al., 2021) is the cosine similarity between normalised CLIP-ViT-H-14 image and text features. PickScore (Kirstain et al., 2023) is exp⁡(τ) (τ) times the same cosine under its own finetuned encoder, τ being the checkpoint’s learned logit scale. ImageReward (Xu et al., 2023) is the score its checkpoint already standardises to zero mean and unit variance. All three are evaluated on the decoded image, so every reward gradient is taken through the VAE decoder as well as the scorer. The three scales are unrelated, so a reward gain is only ever compared within a reward. No constraints are active in this experiment. Evaluation suite. Table 1 and the frontier figures are computed on twelve problems: three reward models × four (prompt, seed) instances. The instances are “a portrait of a man” → “a smiling face” at seeds 4242 and 77, “a photo of a sleeping cat” → “a photo of a cat with wide open eyes” at seed 4242, and “a field of green grass” → “a field of red poppies” at seed 4242. Within a problem all three methods share the initial latent and the λ grid, so each problem contributes one matched triple of frontiers, and the ± in Table 1 is the standard error over the four instances of that reward. Figure 2 (left) is the PickScore man → smiling instance at seed 4242; Figure 4 shows two of the others. The image rows in Figure 2 (right) are selected at matched achieved reward rather than at matched λ, and use three editing pairs: “a portrait of a man” → “a smiling face” (seed 77), “a portrait of a young man” → “a portrait of an old man” (seed 317317), and “a portrait of a woman with straight hair” → “a portrait of a woman with curly hair” (seed 12891289). The latter two come from a separate editing sweep and do not enter Table 1. Control methods. All three methods share the sampler, the initial latent and the reward. We run ConicFlow in its (RH-T) form with H=2H=2; since this experiment is unconstrained, Algorithm 1 reduces to K=1K=1 and its primal response is approximated by ninnern_inner Adam steps at learning rate 10−210^-2, warm-started from the previous window’s solution rather than from the pretrained path. We evaluate Φj=Xj+H,N _j=X_j+H,N through the recursion, i.e. N−j−HN-j-H pretrained steps. Within the window the drift is evaluated without gradient for ℓ>0 >0, dropping the Jb⊤J_b term from the inner gradient. HardFlow, ported here to FLUX latents, instead optimizes the predicted terminal latent x^=Xj,N(xj) x=X_j,N(x_j) obtained from a single call and never differentiates the transformer. FMRG-J takes ninnern_inner Jacobian-guidance steps per sampling step, z←z+(Δt/ninner)λ∇zr(Xj,N(z))z← z+( t/n_inner)\,λ\, _zr (X_j,N(z) ). Choosing λ. The useful range of the control-cost weight λ moves with the reward’s scale and with the problem, so we range-find it per (reward, prompt) pair at a cheap N=4N=4 proxy of the N=8N=8 sweep rather than fixing one grid. We step ×10× 10 down from a ‖∇zr‖\| _zr\|-scaled starting value until Δr r vanishes into the noise and ×10× 10 up until it saturates, then place six λ per method across the logistic rise fitted to Δr(λ) r(λ), densified with log-midpoints where a finer frontier was wanted. The same grid is reused across seeds. The three per-method grids are merged into one, so every arm runs at identical λ, at which all three apply the same first-order step λΔt∇rλ\, t\,∇ r per sampling step. The resulting grids differ by orders of magnitude across the suite; the PickScore man → smiling cell behind Figure 2 (left) runs at λ∈0.95,1.5,2.5,4.9,9.8,15,23,38,62,130,280λ∈\0.95,1.5,2.5,4.9,9.8,15,23,38,62,130,280\. Mean reward gain over a shared energy budget. Comparing methods at a common λ is fragile, since λ enters each method differently and the same value buys different amounts of control. We therefore summarise each method by its whole frontier instead. Define E=∑k=0N−1Δt∥uk∥2E= _k=0^N-1 t\, u_k ^2 to denote a run’s control energy. Fix a problem (reward, prompt, seed). For method m we first discard every run that another run of the same method beats on both axes (i.e., no less energy and no more reward), leaving a set of non-dominated runs (Ei,ri)\(E_i,r_i)\ along which r increases with E; write r^m r_m for its piecewise-linear interpolant in log10⁡E _10E. Because the methods span different energies, we evaluate all of them on the interval attained by all three, E1=maxm⁡miniEiE_1= _m _iE_i and E2=minm⁡maxiEiE_2= _m _iE_i, and report Gm=1log10⁡(E2/E1)∫log10⁡E1log10⁡E2(r^m(E)−runc)dlog10E,G_m\;=\; 1 _10(E_2/E_1) _ _10E_1 _10E_2 ( r_m(E)-r_unc )\;d _10E, (3) the mean height of the frontier over that budget range. GmG_m carries the units of the reward, and Gm=0G_m=0 means the method does no better than the uncontrolled sample anywhere in the range. Both endpoints lie inside every method’s frontier by construction, so no extrapolation is involved, and (3) is invariant to any uniform rescaling of the energy axis. Integrating in log10⁡E _10E rather than in E weights each decade of budget equally; on a linear axis the top decade would carry almost all the weight. ConicFlow is taken at 2424 inner iterations, as plotted in Fig. 2. Two properties of (3) are worth stating. It is conservative for the method that reaches furthest: the interval stops at E2E_2, the lowest of the three maximum energies, so any reward attained beyond it — typically ConicFlow’s — is discarded. And it scores the attainment frontier rather than the sweep itself, so a method whose reward is non-monotone in λ is represented by its upper envelope. Hardware and cost. One run occupies a single NVIDIA RTX A5000 (24 GB), one process per card, peaking at 16.316.3–17.717.7 GB. Median wall clock over the 11 values of λ: 317317 s and 626626 s for ConicFlow at 1212 and 2424 inner iterations, 4646 s for HardFlow, 4343 s for FMRG-J, against 33–44 s for the uncontrolled sample. Table 3: Hyperparameters for the unconstrained image experiment (Fig 2) Shared across methods Steps N / horizon H 88 / 22 Resolution / latent 256×256256× 256 / 16×32×3216× 32× 32 Guidance scale 3.53.5 Precision bf16 sampler, int8 transformer, fp32 reward Reward models CLIP, PickScore (both CLIP-ViT-H-14), ImageReward Seeds / batch 4242 and 77 / 11 λ grid range-found per (reward, prompt) Per method ConicFlow inner iterations 1212 and 2424 (reported as two configurations) Inner optimizer Adam Learning rate 10−210^-2 Terminal estimate rollout of the remaining N−j−HN-j-H steps drift term detached for ℓ>0 >0 HardFlow inner iterations 4040 Learning rate 2.52.5 (replaced by 0.5/reg0.5/reg where lr⋅reg≥1lr·reg≥ 1) terminal estimate single flow-map call to t=0t=0 effort penalty reg=N/(2λ)reg=N/(2λ), active from step 00 FMRG-J guidance steps 33 per sampling step C.1.3 Constrained image experiments implementation details The sampler, resolution, precision, guidance scale, PickScore scorer, N=8N=8, and H=2H=2 are those of Appendix C.1.2; only what differs is described here. Problem. There is no reward: r≡0r≡ 0, so each method’s objective is its control energy alone and the two methods are compared on what it costs them to reach the same feasible set. An instance is a base prompt c, one to three attribute prompts c1,…,cCc_1,…,c_C, and a seed. Constraint i is hi(z)=τi−s⁡(dec⁡(z),ci)≤ 0,h_i(z)\;=\; _i\;-\;s (dec(z),\,c_i )\;≤\;0, (4) where s⁡(⋅,ci)s(·,c_i) is the PickScore of the decoded image against cic_i, exactly the unnormalised score of Appendix C.1.2 with cic_i in place of the edit prompt. Every constraint is evaluated at the predicted terminal latent, and a run is feasible when maxi⁡hi≤0 _ih_i≤ 0 on the final decoded image. Multiple constraints are independent PickScore heads that share nothing but the image. The generation prompt stays c throughout. Thresholds. Absolute PickScore values are not comparable across prompts or seeds, so thresholds are set relative to each instance’s own uncontrolled sample. For each (prompt, seed) pair we first draw the uncontrolled sample from that seed and score it with every constraint to obtain baselines si0s_i^0 (19.219.2–22.522.5 across the instances here). A tightness multiplier m then sets every threshold of the instance at once, τi(m)=si0+(m−1)|si0|=msi0(si0>0), _i(m)\;=\;s_i^0+(m-1)\,|s_i^0|\;=\;m\,s_i^0 (s_i^0>0), (5) so m=1m=1 places the boundary exactly on the uncontrolled image and m=1.5m=1.5 demands a 50%50\% higher score on every attribute simultaneously. (The absolute-value form is what the code implements; it keeps the same grid meaningful for scorers with negative baselines, which does not arise for PickScore.) One (prompt, seed, m) triple is one comparison: both methods receive the same λ, the same thresholds and the same initial latent, and the two uncontrolled images are checked to be byte-identical before the figure is drawn. ConicFlow. The window problem is (RH-T) with the reward term replaced by the Lagrangian term ∑iμihi(e^(xH)) _i _i\,h_i ( e(x_H) ), μi≥0 _i≥ 0. After each window solve one projected dual-ascent step is taken on the constraint value the optimiser saw at its horizon endpoint, μi←max⁡(0,μi+ηhi(e^(xH))), _i\;←\; \! (0,\; _i+η\,h_i( e(x_H)) ), (6) and the multipliers persist across sampling steps. All multipliers start at μinit=0.5 _init=0.5 with η=0.5η=0.5; both were chosen on the pilot instance (“a portrait of a man” → “a smiling face”) and then frozen for every run reported. The remaining difference from the unconstrained runs is the terminal estimate inside the inner objective: rather than rolling the remaining N−j−HN-j-H steps out, e e is a single flow-map call from xHx_H to t=0t=0, which makes the cost of one sampling step linear rather than quadratic in N. Only the surrogate the optimiser differentiates changes; the applied trajectory, the final image and everything scored from it are computed as before. HardFlow. Constraints enter HardFlow’s inner problem as an augmented Lagrangian on the predicted terminal latent x x, as in the released code: the ascent objective is −∑iμihi(x^)−ρ2∑imax(0,hi(x^))2−regα2∥x^−x¯∥2- _i _ih_i( x)- ρ2 _i (0,h_i( x))^2-reg\,α^2\| x- x\|^2, climbed with 4040 explicit gradient steps; every 88 inner steps the multipliers take the projected step (6) with η=0.2η=0.2, and they are reset to 00 at the start of every sampling step (the released code’s semantics, retained). Two settings depart from that code and both were forced by stability. First, the released code doubles ρ whenever a violation exceeds 0.020.02 at a dual update; over 4040 inner steps that lets ρ grow from 55 to 160160 within a single sampling step, and at the tighter thresholds it produced control energies orders of magnitude above ConicFlow’s on the same instance. We hold ρ=5ρ=5 fixed. Second, the inner ascent is an explicit fixed-step scheme whose error contracts by 1−2lr⋅reg⋅α21-2\,lr·reg·α^2 per step, with α→1α→ 1 at the final sampling step; with the authors’ lr=2.5lr=2.5 and reg=N/(2λ)=0.2reg=N/(2λ)=0.2 that factor reaches 00 and the pilot runs diverged, so the inner learning rate is 0.50.5 (lr⋅reg=0.1lr·reg=0.1). The constraint margin and activation fraction are 00, as in Appendix C.1.2. Choosing λ. λ=20λ=20 for both methods throughout. With no reward, λ does not trade reward against effort; it only weights energy against the constraint term and therefore selects which feasible point each method settles on, not whether it reaches one. Instances and sweep. Table 5 lists the instances. Nine prompt pairs carry one or two or three constraints and are swept over a five-point grid of tightness multipliers, from m=1.1m=1.1 to m=1.5m=1.5 in steps of 0.10.1. Every instance is run at both seeds, giving 8686 comparisons, i.e. 172172 controlled runs plus one uncontrolled probe per (prompt, seed). The attribute prompts were written so that the base prompt does not already imply the attribute (e.g., “a young man with straight hair” → “curly hair”). Table 4: Hyperparameters for the constrained image experiments. Everything not listed (sampler, resolution, precision, guidance scale) is as in Table 3. Shared across methods Reward none (r≡0r≡ 0; objective is control energy) Constraint scorer PickScore, one head per attribute prompt Constraint form hi=τi−si≤0h_i= _i-s_i≤ 0, evaluated at the endpoint Thresholds τi=msi0 _i=m\,s_i^0, Eq. (5) Steps N / horizon H 88 / 22 λ 2020 (fixed) Seeds 310310, 698698 Per method ConicFlow inner iterations 2424 inner optimiser Adam, lr 10−210^-2, warm-started terminal estimate one flow-map call xH→t=0x_H→ t=0 (“jump”) drift term detached for ℓ>0 >0 multipliers μinit=0.5 _init=0.5, persist across steps dual step η 0.50.5, once per sampling step, Eq. (6) HardFlow inner iterations 4040 inner lr 0.50.5 (lr⋅reg=0.1lr·reg=0.1) terminal estimate single flow-map call to t=0t=0 effort penalty reg=N/(2λ)=0.2reg=N/(2λ)=0.2, active from step 00 multipliers μinit=0 _init=0, reset every sampling step dual step η 0.20.2, every 88 inner steps penalty weight ρ 55, held fixed (released code doubles it) margin / activation 00 / 00 Table 5: Constrained instances. Each attribute prompt is one PickScore constraint (4); semicolons separate the constraints of an instance. base prompt c attribute prompts c1,…,cCc_1,…,c_C a portrait of a man a smiling face a photo of a golden retriever a dog wearing sunglasses a city street at noon a city street at night a portrait of a woman a woman with bright red hair a portrait of a woman sad frowning woman portrait of a young man with straight hair man with curly hair; photo of very old man a portrait of a woman a woman wearing glasses; a woman with short blonde hair a house by a lake a wooden log cabin; a house in heavy fog a green garden on an overcast day a garden with red flowers; a garden in bright sunshine; a garden with a stone path C.2 Additional Results Unconstrained. We present additional visualizations of the unconstrained image editing comparisons in Figure 5, similar to Figure 2 in the main paper. We also present two additional pareto plots of instances in Figure 4 similar to the one in Figure 2 but with different reward/prompts which define a different problem. Table 6: Per-problem breakdown of Table 1. Instance (prompt → edit, seed) runcr_unc Budget [E1,E2][E_1,E_2] ConicFlow HardFlow FMRG CLIP man → smiling, 42 0.1490.149 [8.8, 217][8.8,\,217] 0.1510.151 0.1350.135 0.0830.083 man → smiling, 7 0.2100.210 [8.2, 136][8.2,\,136] 0.0940.094 0.0700.070 0.0460.046 sleeping cat → open eyes, 42 0.2470.247 [8.8, 127][8.8,\,127] 0.1060.106 0.1010.101 0.0650.065 green grass → red poppies, 42 0.1270.127 [7.6, 122][7.6,\,122] 0.1710.171 0.1660.166 0.0660.066 mean ± s.e. 0.130±0.0180.130± 0.018 0.118±0.0210.118± 0.021 0.065±0.0080.065± 0.008 ConicFlow wins – 4/44/4 4/44/4 ImageReward man → smiling, 42 −1.611-1.611 [9.3, 25.0][9.3,\,25.0] 2.6482.648 2.2032.203 2.5462.546 man → smiling, 7 0.7910.791 [8.0, 14.1][8.0,\,14.1] 0.9200.920 0.7380.738 0.9380.938 sleeping cat → open eyes, 42 −1.285-1.285 [9.9, 254][9.9,\,254] 1.7161.716 0.6100.610 0.1470.147 green grass → red poppies, 42 −2.282-2.282 [7.4, 1974][7.4,\,1974] 2.9082.908 1.5651.565 1.4831.483 mean ± s.e. 2.048±0.4552.048± 0.455 1.279±0.3741.279± 0.374 1.278±0.5041.278± 0.504 ConicFlow wins – 4/44/4 3/43/4 PickScore man → smiling, 42 19.42619.426 [11.0, 1723][11.0,\,1723] 6.0116.011 2.2642.264 2.5322.532 man → smiling, 7 20.17820.178 [9.9, 1746][9.9,\,1746] 4.9114.911 1.8301.830 2.5612.561 sleeping cat → open eyes, 42 20.30420.304 [9.1, 2934][9.1,\,2934] 5.4755.475 0.3770.377 2.8222.822 green grass → red poppies, 42 18.64618.646 [7.4, 2180][7.4,\,2180] 8.2368.236 2.4722.472 2.7662.766 mean ± s.e. 6.158±0.7286.158± 0.728 1.736±0.4721.736± 0.472 2.670±0.0732.670± 0.073 ConicFlow wins – 4/44/4 4/44/4 ConicFlow wins, all 12 problems – 12/1212/12 11/1211/12 (a) ImageReward, “a photo of a sleeping cat” → “a photo of a cat with wide open eyes”, seed 42. (b) CLIP, “a portrait of a man” → “a smiling face”, seed 7. Figure 4: Reward–effort trade-off on two further editing instances, N=8N=8, H=2H=2. The horizontal axis is control energy ∑k∥uk∥2Δt _k u_k ^2 t and the vertical axis the achieved reward r⁡(xN)r(x_N), so up and to the left is better. Each marker is one value of λ. (a) PickScore, seed 12891289: “a portrait of a young man” → “a portrait of an old man” (b) PickScore, seed 317317: “a woman with straight hair” → “…with curly hair” (c) PickScore, seed 317317: “a portrait of a smiling man” → “…an angry man” (d) PickScore, seed 12891289: “a portrait of a smiling man” → “…an angry man” (e) CLIP, seed 4242: “a portrait of a man” → “a smiling face” (f) ImageReward, seed 77: “a portrait of a man” → “a smiling face” (g) CLIP, seed 4242: “a field of green grass” → “a field of red poppies” (h) ImageReward, seed 4242: “a field of green grass” → “a field of red poppies” (i) CLIP, seed 4242: “a photo of a sleeping cat” → “…with wide open eyes” (j) ImageReward, seed 4242: “a photo of a sleeping cat” → “…with wide open eyes” Figure 5: More examples of unconstrained editing. Each panel is one (reward, prompt, seed) instance at a matched-reward level. Panels 5(g)–5(j) are the harder case: turning a green field red, or opening a sleeping cat’s eyes, asks for the image to be restructured rather than adjusted, and no method delivers it. Constrained We present additional visualizations of the constrained image editing comparisons in Figure 6, 7, and 8, similar to Figure 3 in the main paper. (a) “a portrait of a man” → “a smiling face” (b) “a photo of a golden retriever” → “a dog wearing sunglasses” (c) “a city street at noon” → “a city street at night” Figure 6: More examples of constrained editing. Each panel is one (prompt, seed) instance. Instance details are in Table 5. (a) “a portrait of a woman” → “a woman with bright red hair” (b) “a portrait of a woman” → “sad frowning woman” (c) “portrait of a young man with straight hair” → “man with curly hair; photo of very old man” Figure 7: More examples of constrained editing. Each panel is one (prompt, seed) instance. Instance details are in Table 5. (a) “a portrait of a woman” → “a woman wearing glasses; a woman with short blonde hair” (b) “a house by a lake” → “a wooden log cabin; a house in heavy fog” (c) “a green garden on an overcast day” → “a garden with red flowers; a garden in bright sunshine; a garden with a stone path” Figure 8: More examples of constrained editing. Each panel is one (prompt, seed) instance. Instance details are in Table 5. Appendix D Ablations Three questions are relevant to the comparison in Section 5: (1) whether ConicFlow’s advantage is an artifact of the planning horizon we happened to choose, (2) whether it survives once the baselines are given the same compute, and (3) how much of it comes from where the horizon objective is evaluated. We address all three in this section on the editing task at N=8N=8, on exactly the shared λ grids of Section 5. D.1 Planning horizon ConicFlow plans over H steps and applies the first. We re-ran it at H∈1,4,6,8H∈\1,4,6,8\ against the shipped H=2H=2 configuration: four editing prompts at seed 317317, ninner=24n_inner=24, lr=10−2lr=10^-2, terminal estimate by a single flow-map call (“jump”), six λ per cell, so 2424 operating points per H. HardFlow and FMRG-J have no planning horizon and are not re-run. Table 7: Planning horizon, editing, N=8N=8, seed 317317; Δr r, E and time are medians over the 2424 operating points in each H. Wall clock is measured on shared cards and is indicative only. H Δr r E time 11 3.983.98 78.578.5 207207 s 22 4.454.45 68.768.7 249249 s 44 4.394.39 70.670.6 288288 s 66 4.674.67 90.390.3 321321 s 88 4.954.95 101.4101.4 264264 s Moving from a one-step to a two-step horizon improves both quantities at once: Δr r rises 12%12\% and E falls 12%12\%. Beyond H=2H=2 the two diverge. Reward keeps rising, apart from a plateau at H=4H=4, and reaches its maximum at H=8H=8, 24%24\% above H=1H=1 for 24%24\% more NFE, where NFE=ninner(N⁡(N−1)/2+H)NFE=n_inner (N(N-1)/2+H ). Energy grows much faster over the same range: from H=2H=2 to H=8H=8, E rises 48%48\% while Δr r gains only 11%11\%. Longer horizons therefore buy reward with disproportionately more control effort, and on the objective H=2H=2 is the best of the five. D.2 Matched wall clock ConicFlow is the only arm that backpropagates through the transformer, so at equal λ it also spends more compute: 239239 s against 3636 s for HardFlow and 3333 s for FMRG-J in Section 5. We removed that advantage by retuning ninnern_inner to 1212, 167167 and 1313 respectively, which brings all three to ≈122≈\!122 s per run, and re-ran the four editing prompts at seeds 317317 and 12891289 on the same shared λ grids: 4848 operating points at equal λ and equal compute. Table 8: Equal λ and equal wall clock, editing, N=8N=8, seeds 317317 and 12891289; Δr r, E and time are medians. ninnern_inner time Δr r E ConicFlow 1212 124124 s 4.164.16 84.684.6 HardFlow 167167 123123 s 2.392.39 555.2555.2 FMRG-J 1313 128128 s 1.841.84 32.532.5 At matched compute ConicFlow still delivers the largest reward gain, a median Δr r of 4.164.16 against 2.392.39 for HardFlow and 1.841.84 for FMRG-J (Table 8). The two baselines fall short in opposite directions. HardFlow, given more than three times its shipped compute, reaches only 57%57\% of ConicFlow’s reward gain and pays 6.6×6.6× the control energy for it (E=555.2E=555.2 against 84.684.6). FMRG-J spends the least energy of any arm, E=32.5E=32.5, but only because it barely moves the sample: its reward gain is under half of ConicFlow’s. ConicFlow is the only arm that combines a large reward gain with moderate energy, and extra compute alone does not move either baseline there. D.3 Objective evaluation point The RH-T variant of ConicFlow scores the horizon objective at the flow-map prediction of the terminal latent rather than at the planning-window state xHx_H it actually reaches. We ran a comparison between RH-H and RH-T variants. We used four editing prompts at 2 seeds, ninner=24n_inner=24, lr=10−2lr=10^-2, six λ values per cell, so 4848 operating points. As seen in Table 9, evaluating at the terminal point is worth more than the compute. RH-H runs 2.6×2.6× faster but costs 1.541.54 in median Δr r while raising median energy by a quarter, so the RH-H is dominated by RH-T on both axes. Table 9: Where the horizon objective is evaluated, editing, N=8N=8, seeds 317317 and 12891289; Δr r, E and time are medians. objective at time Δr r E ConicFlow endpoint 238238 s 4.794.79 68.768.7 horizon 9090 s 3.253.25 86.986.9 HardFlow — 3636 s 2.002.00 553.9553.9 FMRG-J — 3333 s 2.392.39 87.087.0

Instant research discovery