problem

Core optimization problem interface for trajectory optimization.

This module provides the Problem class, the main entry point for defining and solving trajectory optimization problems using Sequential Convex Programming (SCP).

Example

The prototypical flow is to define a problem, then initialize, solve, and post-process the results

problem = Problem(dynamics, constraints, states, controls, N, time)
problem.initialize()
result = problem.solve()
result = problem.post_process()

Problem

Source code in openscvx/problem.py

class Problem:
    def __init__(
        self,
        dynamics: dict,
        constraints: List[Union[Constraint, CTCS]],
        states: List[State],
        controls: List[Control],
        N: int,
        time: Time,
        *,
        dynamics_discrete: Optional[dict] = None,
        dynamics_prop: Optional[dict] = None,
        states_prop: Optional[List[State]] = None,
        algebraic_prop: Optional[dict] = None,
        licq_min: Union[float, Dict[int, float]] = 0.0,
        licq_max: Union[float, Dict[int, float]] = 1e-4,
        algorithm: Optional[Union[Algorithm, dict]] = None,
        discretizer: Optional[Union[Discretizer, dict]] = None,
        solver: Optional[Union[ConvexSolver, dict]] = None,
        byof: Optional[Union[ByofSpec, dict]] = None,
        float_dtype: str = "float32",
    ):
        """The primary class in charge of compiling and exporting the solvers.

        Args:
            dynamics (dict): Dictionary mapping state names to their dynamics expressions.
                Each key should be a state name, and each value should be an Expr
                representing the derivative of that state.
            constraints (List[Union[CTCSConstraint, NodalConstraint]]):
                List of constraints decorated with @ctcs or @nodal
            states (List[State]): List of State objects representing the state variables.
                May optionally include a State named "time" (see time parameter below).
            controls (List[Control]): List of Control objects representing the control variables
            N (int): Number of segments in the trajectory
            time (Time): Time configuration object with initial, final, min, max.
                Required. If including a "time" state in states, the Time object will be ignored
                and time properties should be set on the time State object instead.
            dynamics_prop (dict, optional): Dictionary mapping EXTRA state names to their
                dynamics expressions for propagation. Only specify additional states beyond
                optimization states (e.g., {"distance": speed}). Do NOT duplicate optimization
                state dynamics here.
            states_prop (List[State], optional): List of EXTRA State objects for propagation only.
                Only specify additional states beyond optimization states. Used with dynamics_prop.
            algebraic_prop (dict, optional): Dictionary mapping names to symbolic expressions
                for outputs evaluated (not integrated) during propagation.
            licq_min: Minimum LICQ constraint value. Defaults to 0.0.
                Either a scalar (applied to all CTCS groups) or a dict
                mapping CTCS group ``idx`` to per-group bounds.
            licq_max: Maximum LICQ constraint value. Defaults to 1e-4.
                Either a scalar (applied to all CTCS groups) or a dict
                mapping CTCS group ``idx`` to per-group bounds.
            algorithm: SCP algorithm configuration. Accepts:

                - ``None`` — uses ``PenalizedTrustRegion()`` with defaults.
                - An ``Algorithm`` instance — used directly.
                - A ``dict`` — passed as kwargs to ``PenalizedTrustRegion()``.
                  Supports a nested ``autotuner`` key in any of these forms:

                  - **string** — class name with default parameters, e.g.
                    ``"RampProximalWeight"``.
                  - **dict** — class name via ``"type"`` key plus parameter
                    overrides, e.g.
                    ``{"type": "RampProximalWeight", "ramp_factor": 1.04}``.
                  - **instance** — an already-constructed autotuner object,
                    e.g. ``ox.RampProximalWeight(ramp_factor=1.04)``.

                The ``lam_cost`` key accepts either a float (applied
                uniformly to all minimize/maximize states) or a dict
                mapping state names to per-state weights::

                    # Uniform cost weight
                    algorithm={"lam_cost": 5e-1}

                    # Per-state cost weights
                    algorithm={"lam_cost": {"velocity": 1e-1, "time": 1e0}}

                    # Per-component weights for vector states
                    algorithm={"lam_cost": {"position": [0, 0, 1e-6], "fuel": 1e0}}

                When a dict is provided, every state that has a
                minimize/maximize objective must have an entry.  Dict
                values may be scalars (broadcast to all components) or
                arrays matching the state's shape.  States without
                objectives are automatically assigned weight 0.  The
                dict is expanded to an array of shape ``(n_states,)``
                during ``Problem`` construction.

                Examples::

                    # Just tweak weights (default algorithm & autotuner)
                    algorithm={"lam_cost": 5e-1, "k_max": 50}

                    # Autotuner by name (default parameters)
                    algorithm={"autotuner": "RampProximalWeight"}

                    # Autotuner as dict with overrides
                    algorithm={
                        "lam_cost": 5e-1,
                        "autotuner": {"type": "RampProximalWeight", "ramp_factor": 1.04},
                    }

                    # Autotuner as instance
                    algorithm={"autotuner": ox.RampProximalWeight(ramp_factor=1.04)}
            discretizer: Discretization method configuration. Accepts:

                - ``None`` — uses ``LinearizeDiscretize()`` with defaults
                  (FOH, Tsit5). Uses sparse Jacobians and compact variational
                  integration when sparsity patterns exist on dynamics; otherwise
                  falls back to dense ``jax.jacfwd`` (same numerics).
                - A ``Discretizer`` instance — used directly.
                - A ``dict`` — resolved via :func:`~openscvx.discretization._resolve_discretizer`
                  (default class is ``LinearizeDiscretize`` unless ``"type"`` is set).

                Examples::

                    # Per-control hold + ODE solver on the discretizer
                    #   thrust = ox.Control("thrust", shape=(3,), parameterization="ZOH")
                    discretizer={"ode_solver": "Dopri8"}

                    # Pass integrator kwargs (forwarded to Diffrax / diffeqsolve)
                    discretizer={
                        "diffrax_kwargs": {
                            "max_steps": 20_000,
                            "num_substeps": 100,
                        }
                    }

                    # Instance
                    discretizer=ox.LinearizeDiscretize(dis_type="ZOH", ode_solver="Dopri8")
            solver: Convex subproblem solver configuration. Accepts:

                - ``None`` — uses ``PTRSolver()`` with defaults (QOCO backend).
                - A ``ConvexSolver`` instance — used directly.
                - A ``dict`` — passed as kwargs to ``PTRSolver()``.

                Examples::

                    # Change CVXPY backend solver and tolerances
                    solver={"cvx_solver": "CLARABEL", "solver_args": {"tol_gap_abs": 1e-7}}

                    # Just change solver_args
                    solver={"solver_args": {"abstol": 1e-6, "reltol": 1e-9}}

                    # Enable cvxpygen code generation
                    solver={"cvxpygen": True}

                    # Instance
                    solver=ox.PTRSolver(cvx_solver="CLARABEL")
            byof (ByofSpec, optional): Expert mode only. Raw JAX functions to
                bypass symbolic layer. See :class:`openscvx.expert.ByofSpec` for
                detailed documentation.
            float_dtype (str): Default floating-point dtype for JAX lowering.
                Must be ``\"float32\"`` or ``\"float64\"``. This sets JAX's
                ``jax_enable_x64`` config flag (True for float64, False for float32),
                which controls the dtype used in all lowered JAX functions (including
                both branches of ``jax.lax.cond``) to avoid dtype-mismatch errors
                during integration.

        Note:
            There are two approaches for handling time:
            1. Auto-create (simple): Don't include "time" in states, provide Time object
            2. User-provided (for time-dependent constraints): Include "time" State in states and
               in dynamics dict, don't provide Time object
        """

        # Configure JAX's default dtype (float32 vs float64) via jax_enable_x64.
        # This must happen before lowering so that all JAX-based lowerers
        # (including conditionals) produce tensors with a consistent dtype.
        # jax_enable_x64=True means float64, jax_enable_x64=False means float32.
        enable_x64 = float_dtype.lower() in ("float64", "f64", "double")
        jax.config.update("jax_enable_x64", enable_x64)

        # Also set the dtype in the JAX lowerer module so it's available during lowering
        # This ensures conditionals use the correct dtype even if JAX config doesn't take effect
        from openscvx.symbolic.lowerers.jax.logic import set_default_float_dtype

        set_default_float_dtype(float_dtype)

        # Persist so integration tests (and callers) can re-sync process-wide JAX config
        # before initialize()/solve() after other examples have been imported.
        self._float_dtype: str = float_dtype

        # Symbolic Preprocessing & Augmentation
        # Resolve byof: dict → ByofSpec (validates keys and nested specs)
        if byof is not None:
            byof = ByofSpec.model_validate(byof)

        self.symbolic: SymbolicProblem = preprocess_symbolic_problem(
            dynamics=dynamics,
            dynamics_discrete=dynamics_discrete,
            constraints=constraints,
            states=states,
            controls=controls,
            N=N,
            time=time,
            licq_min=licq_min,
            licq_max=licq_max,
            dynamics_prop_extra=dynamics_prop,
            states_prop_extra=states_prop,
            algebraic_prop=algebraic_prop,
            byof=byof,
        )

        # Validate byof early (after preprocessing, before lowering) to fail fast
        if byof is not None:
            from openscvx.expert.validation import validate_byof

            # Calculate unified state and control dimensions from preprocessed states/controls
            # These dimensions include symbolic augmentation (time, CTCS) but not byof CTCS
            # augmentation, which is exactly what user byof functions will see
            n_x = sum(
                state.shape[0] if len(state.shape) > 0 else 1 for state in self.symbolic.states
            )
            n_u = sum(
                control.shape[0] if len(control.shape) > 0 else 1
                for control in self.symbolic.controls
            )

            validate_byof(byof, self.symbolic.states, n_x, n_u, N, self.symbolic.parameters)

        # Store byof for cache hashing
        self._byof = byof

        # Resolve algorithm: instance → use directly, dict/None → validate & build
        if isinstance(algorithm, Algorithm):
            self._algorithm = algorithm
        else:
            config = PenalizedTrustRegionConfig.model_validate(algorithm or {})
            self._algorithm = config.to_algorithm(
                states=self.symbolic.states, controls=self.symbolic.controls
            )

        # Resolve discretizer: instance → use directly, dict/None → validate & build
        if isinstance(discretizer, Discretizer):
            self._discretizer = discretizer
        else:
            spec = resolve_discretizer_config(discretizer or {})
            self._discretizer = spec.build()

        # Resolve solver: instance → use directly, dict/None → validate & build
        if isinstance(solver, ConvexSolver):
            self._solver = solver
        else:
            spec = resolve_solver_config(solver or {})
            self._solver = spec.build()

        # Lower to JAX and CVXPy (byof handling happens inside lower_symbolic_problem)
        self._lowered: LoweredProblem = lower_symbolic_problem(
            self.symbolic, self._solver, byof=byof
        )

        # Store parameters in two forms:
        self._parameters = self.symbolic.parameters  # Plain dict for JAX functions
        # Wrapper dict for user access that auto-syncs
        self._parameter_wrapper = ParameterDict(self, self._parameters, self.symbolic.parameters)

        # Setup SCP Configuration
        self.settings = Config(
            sim=SimConfig(
                x=self._lowered.x_unified,
                x_prop=self._lowered.x_prop_unified,
                u=self._lowered.u_unified,
                total_time=self._lowered.x_unified.initial[self._lowered.x_unified.time_slice][0],
                n=N,
                n_states=self._lowered.x_unified.initial.shape[0],
                n_states_prop=self._lowered.x_prop_unified.initial.shape[0],
                ctcs_node_intervals=self.symbolic.node_intervals,
            ),
            dev=DevConfig(),
            prp=PropagationConfig(),
        )

        # Copy time grid setting from Time to sim config so the solver can
        # read it during constraint assembly.
        if isinstance(time, Time):
            self.settings.sim._uniform_time_grid = time.uniform_time_grid

        self._discretization_solver: callable = None
        self._discretization_solver_impulsive: callable = None

        # Set up emitter & queue (thread started in initialize() after columns are known)
        if self.settings.dev.printing:
            self.print_queue = queue.Queue()
            self.emitter_function = lambda data: self.print_queue.put(data)
            self.print_thread = None  # Started in initialize()
        else:
            # no-op emitter; nothing ever gets queued or printed
            self.print_queue = None
            self.emitter_function = lambda data: None
            self.print_thread = None

        # Columns for printing (set in initialize() based on algorithm + autotuner)
        self._columns = None

        self.timing_init = None
        self.timing_solve = None
        self.timing_post = None
        self._profiling_session = None

        # Compiled dynamics (vmapped versions, set in initialize())
        self._compiled_dynamics_prop: Optional[Dynamics] = None

        # Compiled constraints (JIT-compiled versions, set in initialize())
        self._compiled_constraints: Optional[LoweredJaxConstraints] = None

        # Solver state (created fresh for each solve)
        self._state: Optional[AlgorithmState] = None

        # Final solution state (saved after successful solve)
        self._solution: Optional[AlgorithmState] = None

        # SCP algorithm (resolved from `algorithm` parameter above)

    @property
    def solver(self) -> ConvexSolver:
        """Access the convex subproblem solver instance.

        Attributes such as ``cvx_solver``, ``solver_args``, ``cvxpygen``, and
        ``cvxpygen_override`` can be modified freely before ``initialize``
        is called::

            problem.solver.solver_args = {"abstol": 1e-6, "reltol": 1e-9}
            problem.solver.cvxpygen = True
            problem.initialize()

        !!! warning
            Solver settings are compiled into the solve function during
            ``initialize()``.  Changes made **after** ``initialize()``
            will have no effect on subsequent solves.

        Returns:
            The solver instance (e.g., PTRSolver).
        """
        return self._solver

    @property
    def algorithm(self) -> Algorithm:
        """Access the SCP algorithm instance.

        Returns:
            The algorithm instance (e.g., PenalizedTrustRegion).
        """
        return self._algorithm

    @property
    def discretizer(self) -> Discretizer:
        """Access the discretizer instance.

        Attributes such as `dis_type`, `ode_solver`, and `diffrax_kwargs`
        can be modified freely before `initialize` is called:

            problem.discretizer.dis_type = "ZOH"
            problem.discretizer.ode_solver = "Dopri8"
            problem.discretizer.diffrax_kwargs = {"num_substeps": 100}
            problem.initialize()

        !!! warning
            Discretizer settings are compiled into the JIT-cached solver
            during `initialize`.  Changes made **after**
            `initialize()` will have no effect on subsequent solves.

        Returns:
            The discretizer instance (e.g., ``LinearizeDiscretize``).
        """
        return self._discretizer

    @property
    def parameters(self) -> ParameterDict:
        """Get the parameters dictionary.

        The returned dictionary automatically syncs to CVXPy when modified:
            problem.parameters["obs_radius"] = 2.0  # Auto-syncs to CVXPy
            problem.parameters.update({"gate_0_center": center})  # Also syncs

        Returns:
            ParameterDict: Special dict that syncs to CVXPy on assignment.
        """
        return self._parameter_wrapper

    @parameters.setter
    def parameters(self, new_params: dict):
        """Replace the entire parameters dictionary and sync to CVXPy.

        Args:
            new_params: New parameters dictionary
        """
        self._parameters = dict(new_params)  # Create new plain dict
        self._parameter_wrapper = ParameterDict(self, self._parameters, new_params)
        self._sync_parameters()

    def _sync_parameters(self):
        """Sync all parameter values to CVXPy parameters."""
        if self._lowered is None:
            return

        if self._lowered.cvxpy_params is not None:
            for name, value in self._parameter_wrapper.items():
                if name in self._lowered.cvxpy_params:
                    self._lowered.cvxpy_params[name].value = value

    def _sync_boundary_conditions(self):
        """Sync boundary conditions from State objects to lowered representation.

        This method reads the current `.initial` and `.final` values and types
        from the original State objects and updates both the unified state
        representation and the CVXPy solver parameters. This enables workflows
        where initial conditions change between solves.

        Note:
            Safe to call before initialize() - it will simply do nothing.
        """
        if self._lowered is None:
            return

        # Sync initial/final values and types from State objects to optimization unified
        # representation
        for state in self.symbolic.states:
            if state.initial is not None:
                self._lowered.x_unified.initial[state._slice] = state.initial
                self._lowered.x_unified.initial_type[state._slice] = state.initial_type
            if state.final is not None:
                self._lowered.x_unified.final[state._slice] = state.final
                self._lowered.x_unified.final_type[state._slice] = state.final_type

        # Sync initial/final values and types to propagation unified representation
        # (states_prop includes optimization states, so we sync those too)
        for state in self.symbolic.states_prop:
            if state.initial is not None:
                self._lowered.x_prop_unified.initial[state._slice] = state.initial
                self._lowered.x_prop_unified.initial_type[state._slice] = state.initial_type
            if state.final is not None:
                self._lowered.x_prop_unified.final[state._slice] = state.final
                self._lowered.x_prop_unified.final_type[state._slice] = state.final_type

        # Update CVXPy solver parameters (only if solver is initialized)
        if self._solver._problem is not None:
            self._solver.update_boundary_conditions(
                x_init=self._lowered.x_unified.initial,
                x_term=self._lowered.x_unified.final,
            )

    def _sync_guesses(self):
        """Sync trajectory guesses from State/Control objects to lowered representation.

        This method reads the current `.guess` values from the original State and
        Control objects and updates the unified representations. This enables warm-starting
        workflows where the initial trajectory guess is updated between solves (e.g., shifting
        the previous solution).

        Note:
            This only updates the unified representation. The AlgorithmState is
            created from these values in reset() or initialize(), so this must
            be called before those methods to take effect.
        """
        if self._lowered is None:
            return

        # Sync optimization state guesses
        for state in self.symbolic.states:
            if state.guess is not None:
                self._lowered.x_unified.guess[:, state._slice] = state.guess

        # Sync propagation state guesses (includes optimization states)
        for state in self.symbolic.states_prop:
            if state.guess is not None:
                self._lowered.x_prop_unified.guess[:, state._slice] = state.guess

        # Sync control guesses
        for control in self.symbolic.controls:
            if control.guess is not None:
                self._lowered.u_unified.guess[:, control._slice] = control.guess

    def sync(self):
        """Sync parameters and boundary conditions to the solver.

        Call this after modifying State.initial/final or parameters when using
        step() without reset(). This allows warm-starting from the previous
        solution while updating problem data.

        Note:
            This is automatically called by solve() and reset(). Only needed
            when using step() directly with modified parameters or boundary
            conditions between iterations.

        Example:
            MPC with warm-starting::

                problem.initialize()
                while running:
                    # Update initial condition from measurement
                    pos.initial = measured_state
                    problem.sync()  # Sync without resetting algorithm state

                    # Continue from previous solution (warm-start)
                    for _ in range(max_iters):
                        if problem.step()["converged"]:
                            break
        """
        self._sync_parameters()
        self._sync_boundary_conditions()

    @property
    def state(self) -> Optional[AlgorithmState]:
        """Access the current solver state.

        The solver state contains all mutable state from the SCP iterations,
        including current guesses, costs, weights, and history.

        Returns:
            AlgorithmState if initialized, None otherwise

        Example:
            When using `Problem.step()` can use the state to check convergence _etc._

                problem.initialize()
                problem.step()
                print(f"Iteration {problem.state.k}, J_tr={problem.state.J_tr}")
        """
        return self._state

    @property
    def lowered(self) -> LoweredProblem:
        """Access the lowered problem containing JAX/CVXPy objects.

        Returns:
            LoweredProblem with dynamics, constraints, unified interfaces, and CVXPy vars
        """
        return self._lowered

    @property
    def x_unified(self):
        """Unified state interface (delegates to lowered.x_unified)."""
        return self._lowered.x_unified

    @property
    def u_unified(self):
        """Unified control interface (delegates to lowered.u_unified)."""
        return self._lowered.u_unified

    @property
    def slices(self) -> dict[str, slice]:
        """Get mapping of state and control names to their slices in unified vectors.

        This property returns a dictionary mapping each state and control variable name
        to its slice in the respective unified vector. This is particularly useful for
        expert users working with byof (bring-your-own functions) who need to manually
        index into the unified x and u vectors.

        Returns:
            Dictionary mapping variable names to slice objects.
                State variables map to slices in the x vector.
                Control variables map to slices in the u vector.

        Example:
            Usage with byof::

                problem = ox.Problem(dynamics, states, controls, ...)
                print(problem.slices)
                # {'position': slice(0, 3), 'velocity': slice(3, 6), 'theta': slice(0, 1)}

                # Use in byof functions
                byof = {
                    "nodal_constraints": [
                        lambda x, u, node, params: x[problem.slices["velocity"][0]] - 10.0,
                        lambda x, u, node, params: u[problem.slices["theta"][0]] - 1.57,
                    ]
                }
        """
        slices = {}
        slices.update({state.name: state.slice for state in self.symbolic.states})
        slices.update({control.name: control.slice for control in self.symbolic.controls})
        return slices

    def _format_result(self, state: AlgorithmState, converged: bool) -> OptimizationResults:
        """Format solver state as an OptimizationResults object.

        Converts the internal solver state into a user-facing results object,
        mapping state/control arrays to named fields based on symbolic metadata.

        Args:
            state: The AlgorithmState to extract results from.
            converged: Whether the optimization converged.

        Returns:
            OptimizationResults containing the solution data.
        """
        # Build nodes dictionary with all states and controls
        nodes_dict = {}
        has_impulsive_controls = any(
            control.parameterization == "impulsive" for control in self.symbolic.controls
        )

        # Add all states (user-defined and augmented)
        for sym_state in self.symbolic.states:
            state_nodes = state.x[:, sym_state._slice].copy()
            if has_impulsive_controls:
                state_nodes[0] = self.settings.sim.x.initial[sym_state._slice]
            nodes_dict[sym_state.name] = state_nodes

        # Add all controls (user-defined and augmented)
        for control in self.symbolic.controls:
            nodes_dict[control.name] = state.u[:, control._slice]

        return OptimizationResults(
            converged=converged,
            t_final=state.x[:, self.settings.sim.time_slice][-1],
            nodes=nodes_dict,
            trajectory={},  # Populated by post_process
            _states=self.symbolic.states_prop,  # Use propagation states for trajectory dict
            _controls=self.symbolic.controls,
            X=state.X,  # Single source of truth - x and u are properties
            U=state.U,
            discretization_history=state.V_history,
            J_tr_history=state.J_tr,
            J_vb_history=state.J_vb,
            J_vc_history=state.J_vc,
            TR_history=state.TR_history,
            VC_history=state.VC_history,
            lam_prox_history=state.lam_prox_history.copy(),
            actual_reduction_history=state.actual_reduction_history.copy(),
            pred_reduction_history=state.pred_reduction_history.copy(),
            acceptance_ratio_history=state.acceptance_ratio_history.copy(),
        )

    def initialize(self):
        """Compile dynamics, constraints, and solvers; prepare for optimization.

        This method vmaps dynamics, JIT-compiles constraints, builds the convex
        subproblem, and initializes the solver state. Must be called before solve().

        Example:
            Prior to calling the `.solve()` method it is necessary to initialize the problem

                problem = Problem(dynamics, constraints, states, controls, N, time)
                problem.initialize()  # Compile and prepare
                problem.solve()       # Run optimization
        """
        printing.intro()

        # Create a new profiling session (shared across initialize/solve/post_process)
        self._profiling_session = (
            profiling._create_session() if self.settings.dev.profiling else None
        )
        pr = profiling.profiling_start(self.settings.dev.profiling, self._profiling_session)

        t_0_while = time.time()
        # Ensure scaling matrices are correct
        self.settings.sim.__post_init__()

        # Create compiled (vmapped) propagation dynamics
        # This preserves the original un-vmapped versions in _lowered
        self._compiled_dynamics_prop = Dynamics(
            f=jax.vmap(self._lowered.dynamics_prop.f, in_axes=(0, 0, 0, None)),
        )

        # Create compiled (JIT-compiled) constraints as new instances
        # This preserves the original un-JIT'd versions in _lowered
        # TODO: (haynec) switch to AOT instead of JIT
        compiled_nodal = [
            LoweredNodalConstraint(
                func=jax.jit(c.func),
                grad_g_x=jax.jit(c.grad_g_x),
                grad_g_u=jax.jit(c.grad_g_u),
                nodes=c.nodes,
            )
            for c in self._lowered.jax_constraints.nodal
        ]

        compiled_cross_node = [
            LoweredCrossNodeConstraint(
                func=jax.jit(c.func),
                grad_g_X=jax.jit(c.grad_g_X),
                grad_g_U=jax.jit(c.grad_g_U),
            )
            for c in self._lowered.jax_constraints.cross_node
        ]

        self._compiled_constraints = LoweredJaxConstraints(
            nodal=compiled_nodal,
            cross_node=compiled_cross_node,
            ctcs=self._lowered.jax_constraints.ctcs,  # CTCS aren't JIT-compiled here
        )

        # Generate discretization solver via the discretizer (handles Jacobians + vmapping)
        self._discretization_solver = self._discretizer.get_solver(
            self._lowered.dynamics, self.settings
        )
        self._discretization_solver_impulsive = get_impulsive_discretization_solver(
            self._lowered.dynamics_discrete
        )
        self._propagation_solver = get_propagation_solver(
            self._compiled_dynamics_prop.f, self.settings, self._discretizer
        )

        # Build convex subproblem (solver was created in __init__, variables in lower)
        self._solver.initialize(self._lowered, self.settings)

        # Print problem summary (after solver is initialized so we can access problem stats)
        printing.print_problem_summary(
            self.settings,
            self._lowered,
            self._solver,
            self._algorithm,
            self._discretizer,
        )

        # Get cache file paths using symbolic AST hashing
        # This is more stable than hashing lowered JAX code
        dis_solver_file, prop_solver_file = get_solver_cache_paths(
            self.symbolic,
            dt=self.settings.prp.dt,
            total_time=self.settings.sim.total_time,
            byof=self._byof,
        )

        # Compile the discretization solver
        self._discretization_solver = load_or_compile_discretization_solver(
            self._discretization_solver,
            dis_solver_file,
            self._parameters,  # Plain dict for JAX
            self.settings.sim.n,
            self.settings.sim.n_states,
            self.settings.sim.n_controls,
            save_compiled=self.settings.sim.save_compiled,
            debug=self.settings.dev.debug,
        )

        # Compile the impulsive/discrete discretization solver with the same pipeline.
        # This solver is evaluated on node-wise inputs (x_nodes, u_nodes), shape (N, ...).
        # if has_impulsive and self._discretization_solver_impulsive is not None:
        dis_imp_solver_file = dis_solver_file.with_name(
            f"{dis_solver_file.stem}_impulsive{dis_solver_file.suffix}"
        )
        self._discretization_solver_impulsive = load_or_compile_discretization_solver(
            self._discretization_solver_impulsive,
            dis_imp_solver_file,
            self._parameters,  # Plain dict for JAX
            self.settings.sim.n,
            self.settings.sim.n_states,
            self.settings.sim.n_controls,
            save_compiled=self.settings.sim.save_compiled,
            debug=self.settings.dev.debug,
            name="discrete",
        )

        # Setup propagation solver parameters
        dtau = 1.0 / (self.settings.sim.n - 1)
        dt_max = self.settings.sim.u.max[self.settings.sim.time_dilation_slice][0] * dtau
        self.settings.prp.max_tau_len = int(dt_max / self.settings.prp.dt) + 2

        # Compile the propagation solver
        self._propagation_solver = load_or_compile_propagation_solver(
            self._propagation_solver,
            prop_solver_file,
            self._parameters,  # Plain dict for JAX
            self.settings.sim.n_states_prop,
            self.settings.sim.n_controls,
            self.settings.prp.max_tau_len,
            save_compiled=self.settings.sim.save_compiled,
            debug=self.settings.dev.debug,
        )

        # Build per-constraint lam_vb arrays from symbolic constraints.
        # Deferred to initialize() so that user-set lam_vb values
        # (assigned after Problem construction) are picked up.
        n_byof_nodal = len(self._byof.nodal_constraints) if self._byof else 0
        n_byof_cross = len(self._byof.cross_nodal_constraints) if self._byof else 0
        self._algorithm.weights.build_vb_arrays(
            N=self.symbolic.N,
            nodal_constraints=self.symbolic.constraints.nodal,
            cross_node_constraints=self.symbolic.constraints.cross_node,
            n_byof_nodal=n_byof_nodal,
            n_byof_cross=n_byof_cross,
        )

        # Initialize the SCP algorithm
        print("Initializing the SCvx Subproblem Solver...")
        self._algorithm.initialize(
            self._solver,
            self._discretization_solver,
            self._compiled_constraints,
            self.emitter_function,
            self._parameters,  # For warm-start only
            self.settings,  # For warm-start only
            discretization_solver_impulsive=self._discretization_solver_impulsive,
        )
        print("✓ SCvx Subproblem Solver initialized")

        # Get columns from algorithm (now that autotuner is set) and start print thread
        if self.settings.dev.printing:
            self._columns = self._algorithm.get_columns(self.settings.dev.verbosity)
            self.print_thread = threading.Thread(
                target=printing.intermediate,
                args=(self.print_queue, self.settings, self._columns),
                daemon=True,
            )
            self.print_thread.start()
        else:
            # Printing was disabled after __init__, disable emitter to avoid queue buildup
            self.emitter_function = lambda data: None

        # Create fresh solver state
        self._state = AlgorithmState.from_settings(self.settings, self._algorithm.weights)

        t_f_while = time.time()
        self.timing_init = t_f_while - t_0_while
        print("Total Initialization Time: ", self.timing_init)

        # Prime the propagation solver
        prime_propagation_solver(self._propagation_solver, self._parameters, self.settings)

        profiling.profiling_end(pr, "initialize")

    def reset(self):
        """Reset solver state to re-run optimization from initial conditions.

        Creates fresh AlgorithmState while preserving compiled dynamics and solvers.
        Use this to run multiple optimizations without re-initializing.

        This method automatically syncs:
            - Trajectory guesses from State/Control `.guess` attributes
            - Boundary conditions from State `.initial` and `.final` attributes

        Raises:
            ValueError: If initialize() has not been called yet.

        Example:
            After calling `.step()` it may be necessary to reset the problem back to the initial
            conditions

                problem.initialize()
                result1 = problem.step()
                problem.reset()
                result2 = problem.solve()  # Fresh run with same setup

            MPC with warm-starting from previous solution::

                for measured_state in measurements:
                    # Update initial condition
                    pos.initial = measured_state[:3]

                    # Warm-start: shift previous solution as new guess
                    pos.guess = np.roll(prev_result.nodes["pos"], -1, axis=0)

                    problem.reset()  # Syncs guesses and boundary conditions
                    result = problem.solve()
        """
        if self._compiled_dynamics_prop is None:
            raise ValueError("Problem has not been initialized. Call initialize() first")

        # Sync guesses from State/Control objects (must happen before AlgorithmState creation)
        self._sync_guesses()

        # Sync boundary conditions from State objects
        self._sync_boundary_conditions()

        # Create fresh solver state from settings
        self._state = AlgorithmState.from_settings(self.settings, self._algorithm.weights)

        # Reset solution
        self._solution = None

        # Reset timing
        self.timing_solve = None
        self.timing_post = None

    def step(self) -> dict:
        """Perform a single SCP iteration.

        Designed for real-time plotting and interactive optimization. Performs one
        iteration including subproblem solve, state update, and progress emission.

        Note:
            This method is NOT idempotent - it mutates internal state and advances
            the iteration counter. Use reset() to return to initial conditions.

        Returns:
            dict: Contains "converged" (bool) and current iteration state

        Example:
            Call `.step()` manually in a loop to control the algorithm directly

                problem.initialize()
                while not problem.step()["converged"]:
                    plot_trajectory(problem.state.trajs[-1])
        """
        if self._state is None:
            raise ValueError("Problem has not been initialized. Call initialize() first")

        converged = self._algorithm.step(
            self._state,
            self._parameters,  # May change between steps
            self.settings,  # May change between steps
        )

        # Return dict matching original API
        return {
            "converged": converged,
            "scp_k": self._state.k,
            "scp_J_tr": self._state.J_tr,
            "scp_J_vb": self._state.J_vb,
            "scp_J_vc": self._state.J_vc,
        }

    def solve(
        self,
        max_iters: Optional[int] = None,
        time_limit: Optional[float] = None,
        continuous: bool = False,
    ) -> OptimizationResults:
        """Run the SCP algorithm until convergence or iteration limit.

        Args:
            max_iters: Maximum iterations (default: algorithm.k_max)
            time_limit: Wall-clock time limit in seconds. Overrides
                ``algorithm.t_max`` when provided. ``None`` (default) falls
                back to ``algorithm.t_max``.
            continuous: If True, run all iterations regardless of convergence

        Returns:
            OptimizationResults with trajectory and convergence info
                (call post_process() for full propagation)
        """
        # Sync parameters and boundary conditions before solving
        self._sync_parameters()
        self._sync_boundary_conditions()

        required = [
            self._compiled_dynamics_prop,
            self._compiled_constraints,
            self._solver,
            self._discretization_solver,
            self._state,
        ]
        if any(r is None for r in required):
            raise ValueError("Problem has not been initialized. Call initialize() before solve()")

        # Enable the profiler (reuse session from initialize)
        pr = profiling.profiling_start(self.settings.dev.profiling, self._profiling_session)

        t_0_while = time.time()
        # Print top header for solver results
        if self.settings.dev.printing:
            printing.header(self._columns)

        k_max = max_iters if max_iters is not None else self._algorithm.k_max
        t_max = time_limit if time_limit is not None else self._algorithm.t_max

        while self._state.k <= k_max:
            result = self.step()
            if result["converged"] and not continuous:
                break
            if t_max is not None and (time.time() - t_0_while) >= t_max:
                break

        t_f_while = time.time()
        self.timing_solve = t_f_while - t_0_while

        # Wait for print queue to drain (only if thread is running)
        if self.print_thread is not None and self.print_thread.is_alive():
            while self.print_queue.qsize() > 0:
                time.sleep(0.1)

        # Print bottom footer for solver results
        if self.settings.dev.printing:
            printing.footer(self._columns)

        profiling.profiling_end(pr, "solve")

        # Store solution state
        self._solution = copy.deepcopy(self._state)

        timed_out = t_max is not None and self.timing_solve >= t_max
        return self._format_result(self._state, self._state.k <= k_max and not timed_out)

    def post_process(self) -> OptimizationResults:
        """Propagate solution through full nonlinear dynamics for high-fidelity trajectory.

        Integrates the converged SCP solution through the nonlinear dynamics to
        produce x_full, u_full, and t_full. Call after solve() for final results.

        Returns:
            OptimizationResults with propagated trajectory fields

        Raises:
            ValueError: If solve() has not been called yet.
        """
        if self._solution is None:
            raise ValueError("No solution available. Call solve() first.")

        # Enable the profiler (reuse session from initialize)
        pr = profiling.profiling_start(self.settings.dev.profiling, self._profiling_session)

        # Create result from stored solution state
        result = self._format_result(self._solution, self._solution.k <= self._algorithm.k_max)

        t_0_post = time.time()
        result = propagate_trajectory_results(
            self._parameters,
            self.settings,
            result,
            self._propagation_solver,
            dynamics_discrete=self._lowered.dynamics_discrete.f,
            algebraic_prop=self._lowered.algebraic_prop,
            discretizer=self._discretizer,
        )
        t_f_post = time.time()

        self.timing_post = t_f_post - t_0_post

        # Store the propagated result back into _solution for plotting
        # Store as a cached attribute on the _solution object
        self._solution._propagated_result = result

        # Print results summary
        printing.print_results_summary(
            result, self.timing_post, self.timing_init, self.timing_solve
        )

        profiling.profiling_end(pr, "postprocess")
        return result

    def citation(self) -> str:
        """Return BibTeX citations for all components used in this problem.

        Aggregates citations from the algorithm and other components (discretization,
        convex solver, etc.) Each section is prefixed with a comment indicating which component the
        citation is for.

        Returns:
            Formatted string containing all BibTeX citations with comments.

        Example:
            Print all citations for a problem::

                problem = Problem(dynamics, constraints, states, controls, N, time)
                print(problem.citation())
        """
        sections = []

        sections.append(r"% --- AUTO-GENERATED CITATIONS FOR OPENSCVX CONFIGURATION ---")

        # Algorithm citations
        algo_citations = self._algorithm.citation()
        if algo_citations:
            algo_name = type(self._algorithm).__name__
            header = f"% Algorithm: {algo_name}"
            citations = "\n".join(algo_citations)
            sections.append(f"{header}\n\n{citations}")

        # Solver citations
        solver_citations = self._solver.citation()
        if solver_citations:
            solver_name = type(self._solver).__name__
            header = f"% Convex Solver: {solver_name}"
            citations = "\n".join(solver_citations)
            sections.append(f"{header}\n\n{citations}")

        # Discretization citations
        dis_citations = self._discretizer.citation()
        if dis_citations:
            dis_name = type(self._discretizer).__name__
            header = f"% Discretization: {dis_name}"
            citations = "\n".join(dis_citations)
            sections.append(f"{header}\n\n{citations}")

        sections.append(r"% --- END AUTO-GENERATED CITATIONS")

        return "\n\n".join(sections)

algorithm: Algorithm property

Access the SCP algorithm instance.

Returns:

Type	Description
Algorithm	The algorithm instance (e.g., PenalizedTrustRegion).

discretizer: Discretizer property

Access the discretizer instance.

Attributes such as dis_type, ode_solver, and diffrax_kwargs can be modified freely before initialize is called:

problem.discretizer.dis_type = "ZOH"
problem.discretizer.ode_solver = "Dopri8"
problem.discretizer.diffrax_kwargs = {"num_substeps": 100}
problem.initialize()

Warning

Discretizer settings are compiled into the JIT-cached solver during initialize. Changes made after initialize() will have no effect on subsequent solves.

Returns:

Type	Description
Discretizer	The discretizer instance (e.g., LinearizeDiscretize).

lowered: LoweredProblem property

Access the lowered problem containing JAX/CVXPy objects.

Returns:

Type	Description
LoweredProblem	LoweredProblem with dynamics, constraints, unified interfaces, and CVXPy vars

parameters: ParameterDict property writable

Get the parameters dictionary.

The returned dictionary automatically syncs to CVXPy when modified

problem.parameters["obs_radius"] = 2.0 # Auto-syncs to CVXPy problem.parameters.update({"gate_0_center": center}) # Also syncs

Returns:

Name	Type	Description
ParameterDict	ParameterDict	Special dict that syncs to CVXPy on assignment.

slices: dict[str, slice] property

Get mapping of state and control names to their slices in unified vectors.

This property returns a dictionary mapping each state and control variable name to its slice in the respective unified vector. This is particularly useful for expert users working with byof (bring-your-own functions) who need to manually index into the unified x and u vectors.

Returns:

Type	Description
dict[str, slice]	Dictionary mapping variable names to slice objects. State variables map to slices in the x vector. Control variables map to slices in the u vector.

Example

Usage with byof::

problem = ox.Problem(dynamics, states, controls, ...)
print(problem.slices)
# {'position': slice(0, 3), 'velocity': slice(3, 6), 'theta': slice(0, 1)}

# Use in byof functions
byof = {
    "nodal_constraints": [
        lambda x, u, node, params: x[problem.slices["velocity"][0]] - 10.0,
        lambda x, u, node, params: u[problem.slices["theta"][0]] - 1.57,
    ]
}

solver: ConvexSolver property

Access the convex subproblem solver instance.

Attributes such as cvx_solver, solver_args, cvxpygen, and cvxpygen_override can be modified freely before initialize is called::

problem.solver.solver_args = {"abstol": 1e-6, "reltol": 1e-9}
problem.solver.cvxpygen = True
problem.initialize()

Warning

Solver settings are compiled into the solve function during initialize(). Changes made after initialize() will have no effect on subsequent solves.

Returns:

Type	Description
ConvexSolver	The solver instance (e.g., PTRSolver).

state: Optional[AlgorithmState] property

Access the current solver state.

The solver state contains all mutable state from the SCP iterations, including current guesses, costs, weights, and history.

Returns:

Type	Description
Optional[AlgorithmState]	AlgorithmState if initialized, None otherwise

Example

When using Problem.step() can use the state to check convergence etc.

problem.initialize()
problem.step()
print(f"Iteration {problem.state.k}, J_tr={problem.state.J_tr}")

u_unified property

Unified control interface (delegates to lowered.u_unified).

x_unified property

Unified state interface (delegates to lowered.x_unified).

__init__(dynamics: dict, constraints: List[Union[Constraint, CTCS]], states: List[State], controls: List[Control], N: int, time: Time, *, dynamics_discrete: Optional[dict] = None, dynamics_prop: Optional[dict] = None, states_prop: Optional[List[State]] = None, algebraic_prop: Optional[dict] = None, licq_min: Union[float, Dict[int, float]] = 0.0, licq_max: Union[float, Dict[int, float]] = 0.0001, algorithm: Optional[Union[Algorithm, dict]] = None, discretizer: Optional[Union[Discretizer, dict]] = None, solver: Optional[Union[ConvexSolver, dict]] = None, byof: Optional[Union[ByofSpec, dict]] = None, float_dtype: str = 'float32')

The primary class in charge of compiling and exporting the solvers.

Parameters:

Name	Type	Description	Default
dynamics	dict	Dictionary mapping state names to their dynamics expressions. Each key should be a state name, and each value should be an Expr representing the derivative of that state.	required
constraints	List[Union[CTCSConstraint, NodalConstraint]]	List of constraints decorated with @ctcs or @nodal	required
states	List[State]	List of State objects representing the state variables. May optionally include a State named "time" (see time parameter below).	required
controls	List[Control]	List of Control objects representing the control variables	required
N	int	Number of segments in the trajectory	required
time	Time	Time configuration object with initial, final, min, max. Required. If including a "time" state in states, the Time object will be ignored and time properties should be set on the time State object instead.	required
dynamics_prop	dict	Dictionary mapping EXTRA state names to their dynamics expressions for propagation. Only specify additional states beyond optimization states (e.g., {"distance": speed}). Do NOT duplicate optimization state dynamics here.	None
states_prop	List[State]	List of EXTRA State objects for propagation only. Only specify additional states beyond optimization states. Used with dynamics_prop.	None
algebraic_prop	dict	Dictionary mapping names to symbolic expressions for outputs evaluated (not integrated) during propagation.	None
licq_min	Union[float, Dict[int, float]]	Minimum LICQ constraint value. Defaults to 0.0. Either a scalar (applied to all CTCS groups) or a dict mapping CTCS group idx to per-group bounds.	0.0
licq_max	Union[float, Dict[int, float]]	Maximum LICQ constraint value. Defaults to 1e-4. Either a scalar (applied to all CTCS groups) or a dict mapping CTCS group idx to per-group bounds.	0.0001
algorithm	Optional[Union[Algorithm, dict]]	SCP algorithm configuration. Accepts: None — uses PenalizedTrustRegion() with defaults. An Algorithm instance — used directly. A dict — passed as kwargs to PenalizedTrustRegion(). Supports a nested autotuner key in any of these forms: string — class name with default parameters, e.g. "RampProximalWeight". dict — class name via "type" key plus parameter overrides, e.g. {"type": "RampProximalWeight", "ramp_factor": 1.04}. instance — an already-constructed autotuner object, e.g. ox.RampProximalWeight(ramp_factor=1.04). The lam_cost key accepts either a float (applied uniformly to all minimize/maximize states) or a dict mapping state names to per-state weights:: # Uniform cost weight algorithm={"lam_cost": 5e-1} # Per-state cost weights algorithm={"lam_cost": {"velocity": 1e-1, "time": 1e0}} # Per-component weights for vector states algorithm={"lam_cost": {"position": [0, 0, 1e-6], "fuel": 1e0}} When a dict is provided, every state that has a minimize/maximize objective must have an entry. Dict values may be scalars (broadcast to all components) or arrays matching the state's shape. States without objectives are automatically assigned weight 0. The dict is expanded to an array of shape (n_states,) during Problem construction. Examples:: # Just tweak weights (default algorithm & autotuner) algorithm={"lam_cost": 5e-1, "k_max": 50} # Autotuner by name (default parameters) algorithm={"autotuner": "RampProximalWeight"} # Autotuner as dict with overrides algorithm={ "lam_cost": 5e-1, "autotuner": {"type": "RampProximalWeight", "ramp_factor": 1.04}, } # Autotuner as instance algorithm={"autotuner": ox.RampProximalWeight(ramp_factor=1.04)}	None
discretizer	Optional[Union[Discretizer, dict]]	Discretization method configuration. Accepts: None — uses LinearizeDiscretize() with defaults (FOH, Tsit5). Uses sparse Jacobians and compact variational integration when sparsity patterns exist on dynamics; otherwise falls back to dense jax.jacfwd (same numerics). A Discretizer instance — used directly. A dict — resolved via :func:~openscvx.discretization._resolve_discretizer (default class is LinearizeDiscretize unless "type" is set). Examples:: # Per-control hold + ODE solver on the discretizer # thrust = ox.Control("thrust", shape=(3,), parameterization="ZOH") discretizer={"ode_solver": "Dopri8"} # Pass integrator kwargs (forwarded to Diffrax / diffeqsolve) discretizer={ "diffrax_kwargs": { "max_steps": 20_000, "num_substeps": 100, } } # Instance discretizer=ox.LinearizeDiscretize(dis_type="ZOH", ode_solver="Dopri8")	None
solver	Optional[Union[ConvexSolver, dict]]	Convex subproblem solver configuration. Accepts: None — uses PTRSolver() with defaults (QOCO backend). A ConvexSolver instance — used directly. A dict — passed as kwargs to PTRSolver(). Examples:: # Change CVXPY backend solver and tolerances solver={"cvx_solver": "CLARABEL", "solver_args": {"tol_gap_abs": 1e-7}} # Just change solver_args solver={"solver_args": {"abstol": 1e-6, "reltol": 1e-9}} # Enable cvxpygen code generation solver={"cvxpygen": True} # Instance solver=ox.PTRSolver(cvx_solver="CLARABEL")	None
byof	ByofSpec	Expert mode only. Raw JAX functions to bypass symbolic layer. See :class:openscvx.expert.ByofSpec for detailed documentation.	None
float_dtype	str	Default floating-point dtype for JAX lowering. Must be "float32" or "float64". This sets JAX's jax_enable_x64 config flag (True for float64, False for float32), which controls the dtype used in all lowered JAX functions (including both branches of jax.lax.cond) to avoid dtype-mismatch errors during integration.	'float32'

Note

There are two approaches for handling time: 1. Auto-create (simple): Don't include "time" in states, provide Time object 2. User-provided (for time-dependent constraints): Include "time" State in states and in dynamics dict, don't provide Time object

Source code in openscvx/problem.py

def __init__(
    self,
    dynamics: dict,
    constraints: List[Union[Constraint, CTCS]],
    states: List[State],
    controls: List[Control],
    N: int,
    time: Time,
    *,
    dynamics_discrete: Optional[dict] = None,
    dynamics_prop: Optional[dict] = None,
    states_prop: Optional[List[State]] = None,
    algebraic_prop: Optional[dict] = None,
    licq_min: Union[float, Dict[int, float]] = 0.0,
    licq_max: Union[float, Dict[int, float]] = 1e-4,
    algorithm: Optional[Union[Algorithm, dict]] = None,
    discretizer: Optional[Union[Discretizer, dict]] = None,
    solver: Optional[Union[ConvexSolver, dict]] = None,
    byof: Optional[Union[ByofSpec, dict]] = None,
    float_dtype: str = "float32",
):
    """The primary class in charge of compiling and exporting the solvers.

    Args:
        dynamics (dict): Dictionary mapping state names to their dynamics expressions.
            Each key should be a state name, and each value should be an Expr
            representing the derivative of that state.
        constraints (List[Union[CTCSConstraint, NodalConstraint]]):
            List of constraints decorated with @ctcs or @nodal
        states (List[State]): List of State objects representing the state variables.
            May optionally include a State named "time" (see time parameter below).
        controls (List[Control]): List of Control objects representing the control variables
        N (int): Number of segments in the trajectory
        time (Time): Time configuration object with initial, final, min, max.
            Required. If including a "time" state in states, the Time object will be ignored
            and time properties should be set on the time State object instead.
        dynamics_prop (dict, optional): Dictionary mapping EXTRA state names to their
            dynamics expressions for propagation. Only specify additional states beyond
            optimization states (e.g., {"distance": speed}). Do NOT duplicate optimization
            state dynamics here.
        states_prop (List[State], optional): List of EXTRA State objects for propagation only.
            Only specify additional states beyond optimization states. Used with dynamics_prop.
        algebraic_prop (dict, optional): Dictionary mapping names to symbolic expressions
            for outputs evaluated (not integrated) during propagation.
        licq_min: Minimum LICQ constraint value. Defaults to 0.0.
            Either a scalar (applied to all CTCS groups) or a dict
            mapping CTCS group ``idx`` to per-group bounds.
        licq_max: Maximum LICQ constraint value. Defaults to 1e-4.
            Either a scalar (applied to all CTCS groups) or a dict
            mapping CTCS group ``idx`` to per-group bounds.
        algorithm: SCP algorithm configuration. Accepts:

            - ``None`` — uses ``PenalizedTrustRegion()`` with defaults.
            - An ``Algorithm`` instance — used directly.
            - A ``dict`` — passed as kwargs to ``PenalizedTrustRegion()``.
              Supports a nested ``autotuner`` key in any of these forms:

              - **string** — class name with default parameters, e.g.
                ``"RampProximalWeight"``.
              - **dict** — class name via ``"type"`` key plus parameter
                overrides, e.g.
                ``{"type": "RampProximalWeight", "ramp_factor": 1.04}``.
              - **instance** — an already-constructed autotuner object,
                e.g. ``ox.RampProximalWeight(ramp_factor=1.04)``.

            The ``lam_cost`` key accepts either a float (applied
            uniformly to all minimize/maximize states) or a dict
            mapping state names to per-state weights::

                # Uniform cost weight
                algorithm={"lam_cost": 5e-1}

                # Per-state cost weights
                algorithm={"lam_cost": {"velocity": 1e-1, "time": 1e0}}

                # Per-component weights for vector states
                algorithm={"lam_cost": {"position": [0, 0, 1e-6], "fuel": 1e0}}

            When a dict is provided, every state that has a
            minimize/maximize objective must have an entry.  Dict
            values may be scalars (broadcast to all components) or
            arrays matching the state's shape.  States without
            objectives are automatically assigned weight 0.  The
            dict is expanded to an array of shape ``(n_states,)``
            during ``Problem`` construction.

            Examples::

                # Just tweak weights (default algorithm & autotuner)
                algorithm={"lam_cost": 5e-1, "k_max": 50}

                # Autotuner by name (default parameters)
                algorithm={"autotuner": "RampProximalWeight"}

                # Autotuner as dict with overrides
                algorithm={
                    "lam_cost": 5e-1,
                    "autotuner": {"type": "RampProximalWeight", "ramp_factor": 1.04},
                }

                # Autotuner as instance
                algorithm={"autotuner": ox.RampProximalWeight(ramp_factor=1.04)}
        discretizer: Discretization method configuration. Accepts:

            - ``None`` — uses ``LinearizeDiscretize()`` with defaults
              (FOH, Tsit5). Uses sparse Jacobians and compact variational
              integration when sparsity patterns exist on dynamics; otherwise
              falls back to dense ``jax.jacfwd`` (same numerics).
            - A ``Discretizer`` instance — used directly.
            - A ``dict`` — resolved via :func:`~openscvx.discretization._resolve_discretizer`
              (default class is ``LinearizeDiscretize`` unless ``"type"`` is set).

            Examples::

                # Per-control hold + ODE solver on the discretizer
                #   thrust = ox.Control("thrust", shape=(3,), parameterization="ZOH")
                discretizer={"ode_solver": "Dopri8"}

                # Pass integrator kwargs (forwarded to Diffrax / diffeqsolve)
                discretizer={
                    "diffrax_kwargs": {
                        "max_steps": 20_000,
                        "num_substeps": 100,
                    }
                }

                # Instance
                discretizer=ox.LinearizeDiscretize(dis_type="ZOH", ode_solver="Dopri8")
        solver: Convex subproblem solver configuration. Accepts:

            - ``None`` — uses ``PTRSolver()`` with defaults (QOCO backend).
            - A ``ConvexSolver`` instance — used directly.
            - A ``dict`` — passed as kwargs to ``PTRSolver()``.

            Examples::

                # Change CVXPY backend solver and tolerances
                solver={"cvx_solver": "CLARABEL", "solver_args": {"tol_gap_abs": 1e-7}}

                # Just change solver_args
                solver={"solver_args": {"abstol": 1e-6, "reltol": 1e-9}}

                # Enable cvxpygen code generation
                solver={"cvxpygen": True}

                # Instance
                solver=ox.PTRSolver(cvx_solver="CLARABEL")
        byof (ByofSpec, optional): Expert mode only. Raw JAX functions to
            bypass symbolic layer. See :class:`openscvx.expert.ByofSpec` for
            detailed documentation.
        float_dtype (str): Default floating-point dtype for JAX lowering.
            Must be ``\"float32\"`` or ``\"float64\"``. This sets JAX's
            ``jax_enable_x64`` config flag (True for float64, False for float32),
            which controls the dtype used in all lowered JAX functions (including
            both branches of ``jax.lax.cond``) to avoid dtype-mismatch errors
            during integration.

    Note:
        There are two approaches for handling time:
        1. Auto-create (simple): Don't include "time" in states, provide Time object
        2. User-provided (for time-dependent constraints): Include "time" State in states and
           in dynamics dict, don't provide Time object
    """

    # Configure JAX's default dtype (float32 vs float64) via jax_enable_x64.
    # This must happen before lowering so that all JAX-based lowerers
    # (including conditionals) produce tensors with a consistent dtype.
    # jax_enable_x64=True means float64, jax_enable_x64=False means float32.
    enable_x64 = float_dtype.lower() in ("float64", "f64", "double")
    jax.config.update("jax_enable_x64", enable_x64)

    # Also set the dtype in the JAX lowerer module so it's available during lowering
    # This ensures conditionals use the correct dtype even if JAX config doesn't take effect
    from openscvx.symbolic.lowerers.jax.logic import set_default_float_dtype

    set_default_float_dtype(float_dtype)

    # Persist so integration tests (and callers) can re-sync process-wide JAX config
    # before initialize()/solve() after other examples have been imported.
    self._float_dtype: str = float_dtype

    # Symbolic Preprocessing & Augmentation
    # Resolve byof: dict → ByofSpec (validates keys and nested specs)
    if byof is not None:
        byof = ByofSpec.model_validate(byof)

    self.symbolic: SymbolicProblem = preprocess_symbolic_problem(
        dynamics=dynamics,
        dynamics_discrete=dynamics_discrete,
        constraints=constraints,
        states=states,
        controls=controls,
        N=N,
        time=time,
        licq_min=licq_min,
        licq_max=licq_max,
        dynamics_prop_extra=dynamics_prop,
        states_prop_extra=states_prop,
        algebraic_prop=algebraic_prop,
        byof=byof,
    )

    # Validate byof early (after preprocessing, before lowering) to fail fast
    if byof is not None:
        from openscvx.expert.validation import validate_byof

        # Calculate unified state and control dimensions from preprocessed states/controls
        # These dimensions include symbolic augmentation (time, CTCS) but not byof CTCS
        # augmentation, which is exactly what user byof functions will see
        n_x = sum(
            state.shape[0] if len(state.shape) > 0 else 1 for state in self.symbolic.states
        )
        n_u = sum(
            control.shape[0] if len(control.shape) > 0 else 1
            for control in self.symbolic.controls
        )

        validate_byof(byof, self.symbolic.states, n_x, n_u, N, self.symbolic.parameters)

    # Store byof for cache hashing
    self._byof = byof

    # Resolve algorithm: instance → use directly, dict/None → validate & build
    if isinstance(algorithm, Algorithm):
        self._algorithm = algorithm
    else:
        config = PenalizedTrustRegionConfig.model_validate(algorithm or {})
        self._algorithm = config.to_algorithm(
            states=self.symbolic.states, controls=self.symbolic.controls
        )

    # Resolve discretizer: instance → use directly, dict/None → validate & build
    if isinstance(discretizer, Discretizer):
        self._discretizer = discretizer
    else:
        spec = resolve_discretizer_config(discretizer or {})
        self._discretizer = spec.build()

    # Resolve solver: instance → use directly, dict/None → validate & build
    if isinstance(solver, ConvexSolver):
        self._solver = solver
    else:
        spec = resolve_solver_config(solver or {})
        self._solver = spec.build()

    # Lower to JAX and CVXPy (byof handling happens inside lower_symbolic_problem)
    self._lowered: LoweredProblem = lower_symbolic_problem(
        self.symbolic, self._solver, byof=byof
    )

    # Store parameters in two forms:
    self._parameters = self.symbolic.parameters  # Plain dict for JAX functions
    # Wrapper dict for user access that auto-syncs
    self._parameter_wrapper = ParameterDict(self, self._parameters, self.symbolic.parameters)

    # Setup SCP Configuration
    self.settings = Config(
        sim=SimConfig(
            x=self._lowered.x_unified,
            x_prop=self._lowered.x_prop_unified,
            u=self._lowered.u_unified,
            total_time=self._lowered.x_unified.initial[self._lowered.x_unified.time_slice][0],
            n=N,
            n_states=self._lowered.x_unified.initial.shape[0],
            n_states_prop=self._lowered.x_prop_unified.initial.shape[0],
            ctcs_node_intervals=self.symbolic.node_intervals,
        ),
        dev=DevConfig(),
        prp=PropagationConfig(),
    )

    # Copy time grid setting from Time to sim config so the solver can
    # read it during constraint assembly.
    if isinstance(time, Time):
        self.settings.sim._uniform_time_grid = time.uniform_time_grid

    self._discretization_solver: callable = None
    self._discretization_solver_impulsive: callable = None

    # Set up emitter & queue (thread started in initialize() after columns are known)
    if self.settings.dev.printing:
        self.print_queue = queue.Queue()
        self.emitter_function = lambda data: self.print_queue.put(data)
        self.print_thread = None  # Started in initialize()
    else:
        # no-op emitter; nothing ever gets queued or printed
        self.print_queue = None
        self.emitter_function = lambda data: None
        self.print_thread = None

    # Columns for printing (set in initialize() based on algorithm + autotuner)
    self._columns = None

    self.timing_init = None
    self.timing_solve = None
    self.timing_post = None
    self._profiling_session = None

    # Compiled dynamics (vmapped versions, set in initialize())
    self._compiled_dynamics_prop: Optional[Dynamics] = None

    # Compiled constraints (JIT-compiled versions, set in initialize())
    self._compiled_constraints: Optional[LoweredJaxConstraints] = None

    # Solver state (created fresh for each solve)
    self._state: Optional[AlgorithmState] = None

    # Final solution state (saved after successful solve)
    self._solution: Optional[AlgorithmState] = None

citation() -> str

Return BibTeX citations for all components used in this problem.

Aggregates citations from the algorithm and other components (discretization, convex solver, etc.) Each section is prefixed with a comment indicating which component the citation is for.

Returns:

Type	Description
str	Formatted string containing all BibTeX citations with comments.

Example

Print all citations for a problem::

problem = Problem(dynamics, constraints, states, controls, N, time)
print(problem.citation())

Source code in openscvx/problem.py

def citation(self) -> str:
    """Return BibTeX citations for all components used in this problem.

    Aggregates citations from the algorithm and other components (discretization,
    convex solver, etc.) Each section is prefixed with a comment indicating which component the
    citation is for.

    Returns:
        Formatted string containing all BibTeX citations with comments.

    Example:
        Print all citations for a problem::

            problem = Problem(dynamics, constraints, states, controls, N, time)
            print(problem.citation())
    """
    sections = []

    sections.append(r"% --- AUTO-GENERATED CITATIONS FOR OPENSCVX CONFIGURATION ---")

    # Algorithm citations
    algo_citations = self._algorithm.citation()
    if algo_citations:
        algo_name = type(self._algorithm).__name__
        header = f"% Algorithm: {algo_name}"
        citations = "\n".join(algo_citations)
        sections.append(f"{header}\n\n{citations}")

    # Solver citations
    solver_citations = self._solver.citation()
    if solver_citations:
        solver_name = type(self._solver).__name__
        header = f"% Convex Solver: {solver_name}"
        citations = "\n".join(solver_citations)
        sections.append(f"{header}\n\n{citations}")

    # Discretization citations
    dis_citations = self._discretizer.citation()
    if dis_citations:
        dis_name = type(self._discretizer).__name__
        header = f"% Discretization: {dis_name}"
        citations = "\n".join(dis_citations)
        sections.append(f"{header}\n\n{citations}")

    sections.append(r"% --- END AUTO-GENERATED CITATIONS")

    return "\n\n".join(sections)

initialize()

Compile dynamics, constraints, and solvers; prepare for optimization.

This method vmaps dynamics, JIT-compiles constraints, builds the convex subproblem, and initializes the solver state. Must be called before solve().

Example

Prior to calling the .solve() method it is necessary to initialize the problem

problem = Problem(dynamics, constraints, states, controls, N, time)
problem.initialize()  # Compile and prepare
problem.solve()       # Run optimization

Source code in openscvx/problem.py

def initialize(self):
    """Compile dynamics, constraints, and solvers; prepare for optimization.

    This method vmaps dynamics, JIT-compiles constraints, builds the convex
    subproblem, and initializes the solver state. Must be called before solve().

    Example:
        Prior to calling the `.solve()` method it is necessary to initialize the problem

            problem = Problem(dynamics, constraints, states, controls, N, time)
            problem.initialize()  # Compile and prepare
            problem.solve()       # Run optimization
    """
    printing.intro()

    # Create a new profiling session (shared across initialize/solve/post_process)
    self._profiling_session = (
        profiling._create_session() if self.settings.dev.profiling else None
    )
    pr = profiling.profiling_start(self.settings.dev.profiling, self._profiling_session)

    t_0_while = time.time()
    # Ensure scaling matrices are correct
    self.settings.sim.__post_init__()

    # Create compiled (vmapped) propagation dynamics
    # This preserves the original un-vmapped versions in _lowered
    self._compiled_dynamics_prop = Dynamics(
        f=jax.vmap(self._lowered.dynamics_prop.f, in_axes=(0, 0, 0, None)),
    )

    # Create compiled (JIT-compiled) constraints as new instances
    # This preserves the original un-JIT'd versions in _lowered
    # TODO: (haynec) switch to AOT instead of JIT
    compiled_nodal = [
        LoweredNodalConstraint(
            func=jax.jit(c.func),
            grad_g_x=jax.jit(c.grad_g_x),
            grad_g_u=jax.jit(c.grad_g_u),
            nodes=c.nodes,
        )
        for c in self._lowered.jax_constraints.nodal
    ]

    compiled_cross_node = [
        LoweredCrossNodeConstraint(
            func=jax.jit(c.func),
            grad_g_X=jax.jit(c.grad_g_X),
            grad_g_U=jax.jit(c.grad_g_U),
        )
        for c in self._lowered.jax_constraints.cross_node
    ]

    self._compiled_constraints = LoweredJaxConstraints(
        nodal=compiled_nodal,
        cross_node=compiled_cross_node,
        ctcs=self._lowered.jax_constraints.ctcs,  # CTCS aren't JIT-compiled here
    )

    # Generate discretization solver via the discretizer (handles Jacobians + vmapping)
    self._discretization_solver = self._discretizer.get_solver(
        self._lowered.dynamics, self.settings
    )
    self._discretization_solver_impulsive = get_impulsive_discretization_solver(
        self._lowered.dynamics_discrete
    )
    self._propagation_solver = get_propagation_solver(
        self._compiled_dynamics_prop.f, self.settings, self._discretizer
    )

    # Build convex subproblem (solver was created in __init__, variables in lower)
    self._solver.initialize(self._lowered, self.settings)

    # Print problem summary (after solver is initialized so we can access problem stats)
    printing.print_problem_summary(
        self.settings,
        self._lowered,
        self._solver,
        self._algorithm,
        self._discretizer,
    )

    # Get cache file paths using symbolic AST hashing
    # This is more stable than hashing lowered JAX code
    dis_solver_file, prop_solver_file = get_solver_cache_paths(
        self.symbolic,
        dt=self.settings.prp.dt,
        total_time=self.settings.sim.total_time,
        byof=self._byof,
    )

    # Compile the discretization solver
    self._discretization_solver = load_or_compile_discretization_solver(
        self._discretization_solver,
        dis_solver_file,
        self._parameters,  # Plain dict for JAX
        self.settings.sim.n,
        self.settings.sim.n_states,
        self.settings.sim.n_controls,
        save_compiled=self.settings.sim.save_compiled,
        debug=self.settings.dev.debug,
    )

    # Compile the impulsive/discrete discretization solver with the same pipeline.
    # This solver is evaluated on node-wise inputs (x_nodes, u_nodes), shape (N, ...).
    # if has_impulsive and self._discretization_solver_impulsive is not None:
    dis_imp_solver_file = dis_solver_file.with_name(
        f"{dis_solver_file.stem}_impulsive{dis_solver_file.suffix}"
    )
    self._discretization_solver_impulsive = load_or_compile_discretization_solver(
        self._discretization_solver_impulsive,
        dis_imp_solver_file,
        self._parameters,  # Plain dict for JAX
        self.settings.sim.n,
        self.settings.sim.n_states,
        self.settings.sim.n_controls,
        save_compiled=self.settings.sim.save_compiled,
        debug=self.settings.dev.debug,
        name="discrete",
    )

    # Setup propagation solver parameters
    dtau = 1.0 / (self.settings.sim.n - 1)
    dt_max = self.settings.sim.u.max[self.settings.sim.time_dilation_slice][0] * dtau
    self.settings.prp.max_tau_len = int(dt_max / self.settings.prp.dt) + 2

    # Compile the propagation solver
    self._propagation_solver = load_or_compile_propagation_solver(
        self._propagation_solver,
        prop_solver_file,
        self._parameters,  # Plain dict for JAX
        self.settings.sim.n_states_prop,
        self.settings.sim.n_controls,
        self.settings.prp.max_tau_len,
        save_compiled=self.settings.sim.save_compiled,
        debug=self.settings.dev.debug,
    )

    # Build per-constraint lam_vb arrays from symbolic constraints.
    # Deferred to initialize() so that user-set lam_vb values
    # (assigned after Problem construction) are picked up.
    n_byof_nodal = len(self._byof.nodal_constraints) if self._byof else 0
    n_byof_cross = len(self._byof.cross_nodal_constraints) if self._byof else 0
    self._algorithm.weights.build_vb_arrays(
        N=self.symbolic.N,
        nodal_constraints=self.symbolic.constraints.nodal,
        cross_node_constraints=self.symbolic.constraints.cross_node,
        n_byof_nodal=n_byof_nodal,
        n_byof_cross=n_byof_cross,
    )

    # Initialize the SCP algorithm
    print("Initializing the SCvx Subproblem Solver...")
    self._algorithm.initialize(
        self._solver,
        self._discretization_solver,
        self._compiled_constraints,
        self.emitter_function,
        self._parameters,  # For warm-start only
        self.settings,  # For warm-start only
        discretization_solver_impulsive=self._discretization_solver_impulsive,
    )
    print("✓ SCvx Subproblem Solver initialized")

    # Get columns from algorithm (now that autotuner is set) and start print thread
    if self.settings.dev.printing:
        self._columns = self._algorithm.get_columns(self.settings.dev.verbosity)
        self.print_thread = threading.Thread(
            target=printing.intermediate,
            args=(self.print_queue, self.settings, self._columns),
            daemon=True,
        )
        self.print_thread.start()
    else:
        # Printing was disabled after __init__, disable emitter to avoid queue buildup
        self.emitter_function = lambda data: None

    # Create fresh solver state
    self._state = AlgorithmState.from_settings(self.settings, self._algorithm.weights)

    t_f_while = time.time()
    self.timing_init = t_f_while - t_0_while
    print("Total Initialization Time: ", self.timing_init)

    # Prime the propagation solver
    prime_propagation_solver(self._propagation_solver, self._parameters, self.settings)

    profiling.profiling_end(pr, "initialize")

post_process() -> OptimizationResults

Propagate solution through full nonlinear dynamics for high-fidelity trajectory.

Integrates the converged SCP solution through the nonlinear dynamics to produce x_full, u_full, and t_full. Call after solve() for final results.

Returns:

Type	Description
OptimizationResults	OptimizationResults with propagated trajectory fields

Raises:

Type	Description
ValueError	If solve() has not been called yet.

Source code in openscvx/problem.py

def post_process(self) -> OptimizationResults:
    """Propagate solution through full nonlinear dynamics for high-fidelity trajectory.

    Integrates the converged SCP solution through the nonlinear dynamics to
    produce x_full, u_full, and t_full. Call after solve() for final results.

    Returns:
        OptimizationResults with propagated trajectory fields

    Raises:
        ValueError: If solve() has not been called yet.
    """
    if self._solution is None:
        raise ValueError("No solution available. Call solve() first.")

    # Enable the profiler (reuse session from initialize)
    pr = profiling.profiling_start(self.settings.dev.profiling, self._profiling_session)

    # Create result from stored solution state
    result = self._format_result(self._solution, self._solution.k <= self._algorithm.k_max)

    t_0_post = time.time()
    result = propagate_trajectory_results(
        self._parameters,
        self.settings,
        result,
        self._propagation_solver,
        dynamics_discrete=self._lowered.dynamics_discrete.f,
        algebraic_prop=self._lowered.algebraic_prop,
        discretizer=self._discretizer,
    )
    t_f_post = time.time()

    self.timing_post = t_f_post - t_0_post

    # Store the propagated result back into _solution for plotting
    # Store as a cached attribute on the _solution object
    self._solution._propagated_result = result

    # Print results summary
    printing.print_results_summary(
        result, self.timing_post, self.timing_init, self.timing_solve
    )

    profiling.profiling_end(pr, "postprocess")
    return result

reset()

Reset solver state to re-run optimization from initial conditions.

Creates fresh AlgorithmState while preserving compiled dynamics and solvers. Use this to run multiple optimizations without re-initializing.

This method automatically syncs

• Trajectory guesses from State/Control .guess attributes
• Boundary conditions from State .initial and .final attributes

Raises:

Type	Description
ValueError	If initialize() has not been called yet.

Example

After calling .step() it may be necessary to reset the problem back to the initial conditions

problem.initialize()
result1 = problem.step()
problem.reset()
result2 = problem.solve()  # Fresh run with same setup

MPC with warm-starting from previous solution::

for measured_state in measurements:
    # Update initial condition
    pos.initial = measured_state[:3]

    # Warm-start: shift previous solution as new guess
    pos.guess = np.roll(prev_result.nodes["pos"], -1, axis=0)

    problem.reset()  # Syncs guesses and boundary conditions
    result = problem.solve()

Source code in openscvx/problem.py

def reset(self):
    """Reset solver state to re-run optimization from initial conditions.

    Creates fresh AlgorithmState while preserving compiled dynamics and solvers.
    Use this to run multiple optimizations without re-initializing.

    This method automatically syncs:
        - Trajectory guesses from State/Control `.guess` attributes
        - Boundary conditions from State `.initial` and `.final` attributes

    Raises:
        ValueError: If initialize() has not been called yet.

    Example:
        After calling `.step()` it may be necessary to reset the problem back to the initial
        conditions

            problem.initialize()
            result1 = problem.step()
            problem.reset()
            result2 = problem.solve()  # Fresh run with same setup

        MPC with warm-starting from previous solution::

            for measured_state in measurements:
                # Update initial condition
                pos.initial = measured_state[:3]

                # Warm-start: shift previous solution as new guess
                pos.guess = np.roll(prev_result.nodes["pos"], -1, axis=0)

                problem.reset()  # Syncs guesses and boundary conditions
                result = problem.solve()
    """
    if self._compiled_dynamics_prop is None:
        raise ValueError("Problem has not been initialized. Call initialize() first")

    # Sync guesses from State/Control objects (must happen before AlgorithmState creation)
    self._sync_guesses()

    # Sync boundary conditions from State objects
    self._sync_boundary_conditions()

    # Create fresh solver state from settings
    self._state = AlgorithmState.from_settings(self.settings, self._algorithm.weights)

    # Reset solution
    self._solution = None

    # Reset timing
    self.timing_solve = None
    self.timing_post = None

solve(max_iters: Optional[int] = None, time_limit: Optional[float] = None, continuous: bool = False) -> OptimizationResults

Run the SCP algorithm until convergence or iteration limit.

Parameters:

Name	Type	Description	Default
max_iters	Optional[int]	Maximum iterations (default: algorithm.k_max)	None
time_limit	Optional[float]	Wall-clock time limit in seconds. Overrides algorithm.t_max when provided. None (default) falls back to algorithm.t_max.	None
continuous	bool	If True, run all iterations regardless of convergence	False

Returns:

Type	Description
OptimizationResults	OptimizationResults with trajectory and convergence info (call post_process() for full propagation)

Source code in openscvx/problem.py

def solve(
    self,
    max_iters: Optional[int] = None,
    time_limit: Optional[float] = None,
    continuous: bool = False,
) -> OptimizationResults:
    """Run the SCP algorithm until convergence or iteration limit.

    Args:
        max_iters: Maximum iterations (default: algorithm.k_max)
        time_limit: Wall-clock time limit in seconds. Overrides
            ``algorithm.t_max`` when provided. ``None`` (default) falls
            back to ``algorithm.t_max``.
        continuous: If True, run all iterations regardless of convergence

    Returns:
        OptimizationResults with trajectory and convergence info
            (call post_process() for full propagation)
    """
    # Sync parameters and boundary conditions before solving
    self._sync_parameters()
    self._sync_boundary_conditions()

    required = [
        self._compiled_dynamics_prop,
        self._compiled_constraints,
        self._solver,
        self._discretization_solver,
        self._state,
    ]
    if any(r is None for r in required):
        raise ValueError("Problem has not been initialized. Call initialize() before solve()")

    # Enable the profiler (reuse session from initialize)
    pr = profiling.profiling_start(self.settings.dev.profiling, self._profiling_session)

    t_0_while = time.time()
    # Print top header for solver results
    if self.settings.dev.printing:
        printing.header(self._columns)

    k_max = max_iters if max_iters is not None else self._algorithm.k_max
    t_max = time_limit if time_limit is not None else self._algorithm.t_max

    while self._state.k <= k_max:
        result = self.step()
        if result["converged"] and not continuous:
            break
        if t_max is not None and (time.time() - t_0_while) >= t_max:
            break

    t_f_while = time.time()
    self.timing_solve = t_f_while - t_0_while

    # Wait for print queue to drain (only if thread is running)
    if self.print_thread is not None and self.print_thread.is_alive():
        while self.print_queue.qsize() > 0:
            time.sleep(0.1)

    # Print bottom footer for solver results
    if self.settings.dev.printing:
        printing.footer(self._columns)

    profiling.profiling_end(pr, "solve")

    # Store solution state
    self._solution = copy.deepcopy(self._state)

    timed_out = t_max is not None and self.timing_solve >= t_max
    return self._format_result(self._state, self._state.k <= k_max and not timed_out)

step() -> dict

Perform a single SCP iteration.

Designed for real-time plotting and interactive optimization. Performs one iteration including subproblem solve, state update, and progress emission.

Note

This method is NOT idempotent - it mutates internal state and advances the iteration counter. Use reset() to return to initial conditions.

Returns:

Name	Type	Description
dict	dict	Contains "converged" (bool) and current iteration state

Example

Call .step() manually in a loop to control the algorithm directly

problem.initialize()
while not problem.step()["converged"]:
    plot_trajectory(problem.state.trajs[-1])

Source code in openscvx/problem.py

def step(self) -> dict:
    """Perform a single SCP iteration.

    Designed for real-time plotting and interactive optimization. Performs one
    iteration including subproblem solve, state update, and progress emission.

    Note:
        This method is NOT idempotent - it mutates internal state and advances
        the iteration counter. Use reset() to return to initial conditions.

    Returns:
        dict: Contains "converged" (bool) and current iteration state

    Example:
        Call `.step()` manually in a loop to control the algorithm directly

            problem.initialize()
            while not problem.step()["converged"]:
                plot_trajectory(problem.state.trajs[-1])
    """
    if self._state is None:
        raise ValueError("Problem has not been initialized. Call initialize() first")

    converged = self._algorithm.step(
        self._state,
        self._parameters,  # May change between steps
        self.settings,  # May change between steps
    )

    # Return dict matching original API
    return {
        "converged": converged,
        "scp_k": self._state.k,
        "scp_J_tr": self._state.J_tr,
        "scp_J_vb": self._state.J_vb,
        "scp_J_vc": self._state.J_vc,
    }

sync()

Sync parameters and boundary conditions to the solver.

Call this after modifying State.initial/final or parameters when using step() without reset(). This allows warm-starting from the previous solution while updating problem data.

Note

This is automatically called by solve() and reset(). Only needed when using step() directly with modified parameters or boundary conditions between iterations.

Example

MPC with warm-starting::

problem.initialize()
while running:
    # Update initial condition from measurement
    pos.initial = measured_state
    problem.sync()  # Sync without resetting algorithm state

    # Continue from previous solution (warm-start)
    for _ in range(max_iters):
        if problem.step()["converged"]:
            break

Source code in openscvx/problem.py

def sync(self):
    """Sync parameters and boundary conditions to the solver.

    Call this after modifying State.initial/final or parameters when using
    step() without reset(). This allows warm-starting from the previous
    solution while updating problem data.

    Note:
        This is automatically called by solve() and reset(). Only needed
        when using step() directly with modified parameters or boundary
        conditions between iterations.

    Example:
        MPC with warm-starting::

            problem.initialize()
            while running:
                # Update initial condition from measurement
                pos.initial = measured_state
                problem.sync()  # Sync without resetting algorithm state

                # Continue from previous solution (warm-start)
                for _ in range(max_iters):
                    if problem.step()["converged"]:
                        break
    """
    self._sync_parameters()
    self._sync_boundary_conditions()