The field of parametric programming focuses on constructing mappings from optimization problem parameters to their corresponding solutions (Fiacco, 1983; Pistikopoulos et al., 2002). This approach has been explored across various optimization domains, including linear programming (Gal, 2010), quadratic programming (Bemporad et al., 2002), and mixed-integer nonlinear programming (Dua and Pistikopoulos, 1999). A notable application of parametric programming appears in controller design through explicit Model Predictive Control (MPC) (Bemporad et al., 2002; Tøndel et al., 2003), which involves solving MPC problems offline and storing the resulting active constraint sets. During online computation, the current system parameters are used to identify the appropriate pre-computed solution. Bemporad et al. (2002) demonstrated that for constrained linear quadratic regulator problems, these active sets form polyhedrons. Dua and Pistikopoulos (2000) further proposed efficient parameter space partitioning algorithms to construct those polyhedrons.

For non-convex optimization problems, where active sets generally do not form polyhedrons, researchers have developed different approaches. Several studies (Hauser, 2016; Tang and Hauser, 2019; Zhu and Martius, 2020) employ non-parametric learning techniques such as K-nearest neighbor (KNN) to directly leverage stored solution data for generating warm-start solutions online. A more recent research has expanded these methods by incorporating neural networks to learn embeddings that map from larger parameter sets to optimization solutions (Cauligi et al., 2021). A significant advantage of neural network methods is their capability to generalize to out-of-distribution problem instances that lie considerably outside the training dataset (Power and Berenson, 2022).

Mixed-integer convex programs (MICPs) integrate both continuous and discrete decision variables within an optimization framework. Their formal mathematical definition hinges on the property that if all binary variables are relaxed (i.e., replaced by continuous variables between 0 and 1), the resulting problem must be convex (Conforti et al., 2014). Standard MICP solvers such as Gurobi leverage branch-and-bound methods (Boyd and Mattingley, 2007) and cutting plane techniques (Marchand et al., 2002). Although MICPs face exponential worst-case computational complexity with respect to the number of binary variables, practical experience reveals that many real-world problems can be solved by exploring only a small portion of the potential search tree (Williams, 2013). This has facilitated the deployment of mixed-integer approaches in time-sensitive robotics applications such as online motion planning (Tordesillas et al., 2019). Nevertheless, as problem complexity increases, solution times can still extend to minutes or even hours (Lin et al., 2021), presenting significant challenges for real-time implementation.

Mathematical programs with complementarity constraints (MPCCs) offer an alternative approach to handling discrete variables by reformulating them as continuous variables with additional complementarity conditions. These complementarity constraints enforce that when one variable takes a non-zero value, the other must be exactly zero. This formulation, while compact, presents significant computational challenges because it violates the majority of Constraint Qualifications established for standard nonlinear optimization (Luo et al., 1996). To address these difficulties, researchers have developed specialized algorithms including time-stepping methods (Anitescu and Potra, 1997), pivoting algorithms (Drumwright, 2015), and central path approaches (Kojima et al., 1991). In robotics applications, complementarity constraints have proven valuable for trajectory optimization across different locomotion gaits and implicit contact dynamics (Yunt and Glocker, 2006; Posa et al., 2014; Zhang et al., 2021; Le Cleac’h et al., 2024), where they capture the complementary relationship between contact forces and surface distances. Complementarity constraints can also be used to model binary decision variables (Park and Boyd, 2018).

Consider a set of optimization problems parameterized by Θ ∈ R p , drawn from a distribution D ( Θ ) . We aim to solve a parameterized mixed-integer bilinear program (MIBLP) (Fischetti and Monaci, 2020): minimize x , z ‖ x − x g Θ ‖ Q 2 (1a) s.t. A Θ x ≤ h Θ (1b) l i ≤ x i ≤ u i i = 1 , … , dim x (1c) x j ∈ 0,1 j ∈ J (1d) x r s = x p s x q s s = 1 , … , S (1e)

where x ∈ R dim ( x ) includes both continuous and binary variables, x g denotes the desired target values as a function of the parameters, and J specifies the index set for binary variables in x . Constraints (1b) are mixed-integer linear, meaning that if the binary variables are relaxed to continuous domains (i.e., replacing x [ j ] ∈ { 0,1 } with x [ j ] ∈ [ 0,1 ] , ∀ j ∈ J ), these constraints become linear. Note that equality constraints can be incorporated into (1b) by expressing them as pairs of opposing inequality constraints. Constraints (1c) provide simple bounds on the variables, while (1d) enforces the binary requirements. Finally, (1e) introduces bilinear constraints, the distinguishing feature of MIBLPs compared to mixed-integer linear or quadratic programs.

As introduced in Section 1, we compare the data-driven performance of two reformulation approaches: converting problem (1) into either an MICP or an MPCC, and solving these reformulations with data-driven methods. For the MICP reformulation, when a learning algorithm can predict binary variable assignments based on problem parameters Θ , the original problem (1) reduces to a convex optimization problem that can be efficiently solved with standard convex program solvers. Alternatively, for the MPCC reformulation (leading to an NLP with complementary constraints), when a learning algorithm provides high-quality initial guesses for all variables, it can significantly improve computational performance.

To implement these data-driven approaches, we first construct a dataset of size D consisting of tuples { ( Θ i , x i ) } for i = 1 , … , D . The parameters { Θ i } are drawn from D ( Θ ) , and each corresponding solution x i is obtained by solving the original problem using conventional non-data-driven methods. This offline computation, while potentially time-consuming, needs to be performed only once.

We then use this dataset to train models that can either predict binary variables for the MICP reformulation, or generate initial guesses for all variables for the MPCC reformulation. The MICP approach only requires binary variable prediction since fixing these variables yields a convex programming problem that can be solved efficiently without warm-starting. Formally, we seek to approximate the problem-solution map F :   Θ ↦ x using the collected data. For binary variable prediction, we focus on learning F b :   Θ ↦ { x [ j ] | j ∈ J } , while for initial guess generation, we approximate the full mapping F . In the following sections, we describe the implementation details for each reformulation approach.

In this section, we present two problems with two different robotic models that will serve as benchmarks for evaluating the data-driven MIBLP reformulation approaches described in Section 3.1. The first model is a linear inverted pendulum with soft contact walls. The second model is a single rigid body with mode switching and contact constraints. For each problem, we detail how they can be formulated as MIBLPs following the structure of problem (1).

For our first benchmark problem, we implemented a linearized cart-pole system with soft contact walls on either side. The system consists of a cart with mass m cart = 1.0 kg and an inverted pendulum with mass m pole = 0.4 kg and length l = 0.6 m. The system includes walls on both sides at distances d left = d right = 0.5 m from the origin. When the pendulum tip contacts these walls, a soft contact force is generated proportional to the penetration depth, with stiffness coefficients k 1 = k 2 = 50 N/m for the right and left walls, respectively. The control input is limited to | u | ≤ 20 N, and the pendulum angle is constrained to | θ | ≤ π / 2 . The control objective is to stabilize the pendulum in the upright position x g = [ 0,0,0,0 ] . We use a quadratic cost function with state cost matrix Q = diag ( 1,50,1,50 ) , control cost R = diag ( 1,0.1 , 0.1 ) , and terminal cost Q ter computed from the Discrete Algebraic Riccati equation.

To build our dataset for data-driven acceleration, we created a dataset with 50 and 500 samples distributed across the state space. We generated solutions for a grid of initial conditions, varying cart position: [ − 0.3 , 0.3 ] m, pendulum angle: [ − π / 4 , π / 4 ] rad, cart velocity: [ − 0.5 , 0.5 ] m/s, and angular velocity: [ − 0.5 , 0.5 ] rad/s. For each set of initial conditions, we solved the complete MIBLP formulation using Gurobi (Gurobi Optimization, LLC, 2024) to obtain warm-start solutions.

For testing, we generated 200 random initial conditions different from the training set but drawn from the same distribution. These solutions were stored and used to warm-start both reformulation approaches. For the MPCC reformulation, we use IPOPT (Wachter, 2002) as the NLP solver for the complete formulation with complementary constraints. For the MICP reformulation, we use Gurobi as the convex optimization solver after fixing the binary variables.

Results Table 1 shows the computational performance comparison between different solving approaches. The “No Initial Guess MPCC” column represents solving the MPCC formulation with default zero initialization. The “Manual Initial Guess MPCC” column represents solving MPCC with a straight-line trajectory initialization that linearly interpolates all state variables from the initial state to the final upright equilibrium, with all control inputs and binary variables initialized to zero. The “Enhanced Manual Initial Guess MPCC” column represents solving MPCC with a physics-informed warm-start: if the pole is angling toward a wall, the trajectory assumes the pole tip will first move toward the wall at constant velocity, remain in contact for a specified duration, then move away from the wall at constant velocity; otherwise, it uses the straight-line interpolation. Our data-driven approaches included “Data-driven MPCC” using KNN to generate initial guesses for the MPCC formulation, and “Data-driven MICP” using KNN to predict binary variables for the MICP reformulation.

We evaluated the success rate, solving time, and solution quality for each method. The “Amount of data” row indicates the number of offline-solved problem instances in the dataset (zero for non-data-driven methods). The “Success rate” row shows the percentage of test problems solved successfully: for manual initialization methods, this represents a single solving attempt; for data-driven methods, the format a / b / c indicates cumulative success rates when trying K = 1 / 2 / 3 nearest neighbors sequentially until a solution is found. The “Average solve time” represents the mean computational time across successful problem instances, including all KNN neighbor attempts until success. We also provide a histogram showing the distribution of solving times in subplot (A) of Figure 1, indicating how many problem instances are solved within each time interval for data-driven MPCC with 50 training samples. The distribution for 500 training samples exhibits a similar pattern. For data-driven MICP, all problem instances are solved within 1 m. The “Avg. # of neighbors” indicates the average number of KNN neighbors attempted before finding a successful solution. The remaining rows show the average objective value among successful solutions and the solver used.

Analysis First, MPCC with default zero initialization completely fails to find feasible solutions (0% success rate), while manual straight-line initialization achieves moderate success (42.5%). The enhanced manual initialization achieves a higher success rate of 58.0%. However, designing such heuristics requires problem-specific insight and manual effort that may not scale to more complex robotic systems with higher-dimensional state spaces. The data-driven MPCC approach significantly improves performance, achieving a 76% success rate with just 50 data samples. The data-driven MICP formulation achieves a lower success rate (56%) with the same amount of data. When increasing the dataset to 500 samples, both approaches improve further, with MPCC reaching a 90% success rate and MICP improving to 71.0%. The MICP approach offers nearly a 100 × speedup with 50 samples and 45 × with 500 samples, partly due to the commercial Gurobi solver compared to the open-source IPOPT solver used for the MPCC approach.

It is important to note that the relaxation parameter ϵ plays a crucial role in the MPCC formulation to trade off success rate with numerical accuracy. In this experiment, we selected ϵ = 0.005 for a good balance. All the resultant binary variables are within 0.01 of their true binary values, and the success rate is reasonable. If we were to reduce this value to ϵ = 0.001 , the accuracy of binary solutions is higher at the cost of decreasing the success rate.

For this experiment, we evaluate our MIBLP motion planner on a quadrupedal robot model. We use the Spine-enhanced Climbing Autonomous Legged Exploration Robot (SCALER) for demonstration, and implement an SRB model that matches the robot’s physical properties. The total mass is set to 10.0 k g , the moment of inertia is [ 0.15 , 0.30 , 0.35 ] k g ⋅ m 2 , and the body dimensions are 0.5 m × 0.3 m × 0.2 m . For detailed descriptions of this platform, we refer readers to Tanaka et al. (2022). The SCALER robot is equipped with 6-axis force-torque sensors at each limb, allowing us to implement admittance control that enables tracking of the planned contact forces (Schperberg et al., 2023).

To showcase the advantages of our MIBLP formulation, which can plan trajectories without pre-assigned gait sequences or small angle approximations, we design a challenging scenario. The environment contains a low-clearance obstacle positioned at height h obs between x = 1.5 m and x = 2.5 m that requires the robot to pass underneath. The target is an elevated platform located at height h target > h obs , requiring the robot to reach upward to achieve the goal position. The robot begins in bipedal configuration with an upright posture, then transition into quadrupedal mode to crawl beneath the low-clearance obstacle, and finally return to bipedal configuration to reach the elevated target position. The scenario is shown in Figure 2.

We sampled 30 and 80 initial conditions with horizontal positions spanning a 1-m range and initial orientation varying within [ − π / 6 , π / 6 ] from the vertical orientation. We solved these problems using a two-stage approach: first, we solve a low-accuracy MICP formulation using coarser McCormick envelopes indicated in Table 2 (column “data collection”, p sh , i denotes the shoulder position for leg i relative to the center of mass) to generate initial solutions, then refine them through the MPCC formulation using the MICP solution as a warm-start to ensure bilinear constraints are exactly satisfied. This approach can be viewed as a simplified single-iteration case of iterative schemes such as Alternating Direction Method of Multipliers (ADMM) and Generalized Benders Decomposition (GBD). For applications requiring higher accuracy or when the initial MICP solution fails to converge in the MPCC refinement step, these full iterative schemes can be utilized.

For testing, we generated 100 random initial conditions that were different from the training dataset and used them to evaluate both the MPCC and MICP formulations. Higher accuracy McCormick envelopes are used in these MICP tests to ensure accurate approximation of the bilinear constraints, with the number of segmentations given in Table 2 (column “testing”).

Results Table 3 shows the computational performances of different solving methods. “Manual Initial Guess MPCC” uses a common straight-line trajectory to connect the initial and goal configurations as the initial guess for the MPCC solver. All binary variables are initialized to 0.5 to provide a neutral starting point without biasing toward any particular contact mode. “Two-stage MICP-MPCC” refers to the two-stage approach used for data collection described previously: first solving a coarse MICP formulation, then refining through MPCC with the MICP solution as warm-start. “Data-driven MPCC” uses our KNN approach to generate initial guesses for the MPCC formulation, while “Data-driven MICP” uses KNN to predict binary variables for the MICP reformulation.

Analysis We observe that both data-driven approaches benefit substantially from warm-starting with the collected dataset. Manual initialization achieves a 0% success rate, failing to find feasible solutions for any test cases. The two-stage MICP-MPCC approach used for data collection demonstrates feasibility but requires 13.25 s on average to solve. In contrast, with 30 training samples, MPCC achieves up to 84% success rate while MICP reaches 29%; with 80 training samples, MPCC improves to 86% while MICP achieves 47%.

Given the same amount of data, MPCC shows higher success rates, matching the trend observed in the cart-pole experiment where MPCC also outperformed MICP with limited training data. For computational speed, MICP demonstrates a significant advantage, solving problems in approximately 0.19 s compared to 6.91 − 8.13 s for MPCC, more than 35 × faster. This speedup makes MICP particularly suitable for real-time applications.

According to our experiments on both the linear inverted pendulum and single rigid body control problems, data-driven methods significantly improve the performance of both MPCC and MICP reformulations. As we increased the amount of data, we observed a general trend of improvement in success rates for both approaches. The key distinction between MPCC and MICP reformulations lies in the solving speeds and the amount of data used to achieve similar success rates.

For both the linear inverted pendulum problem and the single rigid body problem, the MPCC formulation achieves higher success rates than the MICP formulation with the same amount of data. We noticed that IPOPT rarely modifies the initial binary guesses for x [ j ] , j ∈ J . Therefore, this difference stems from how the bilinear constraints (1e) are handled. The MICP reformulation converts them into multiple linear envelopes with additional binary variables δ s . The learner directly provides δ s , which fixes the intervals [ x p s i , x p s i + 1 ] and [ x q s j , x q s j + 1 ] within which x [ p s ] and x [ q s ] must lie. The convex solver then finds x [ p s ] and x [ q s ] inside these intervals. If the learner predicts incorrect intervals, the resulting problem becomes infeasible, requiring the learner to try alternative neighbors. Thus, a sufficiently large dataset is necessary to train the learner to predict the correct intervals consistently. In contrast, for MPCC re-formulation, KNN provides initial guesses for x [ p s ] , x [ q s ] , then the NLP solver uses numerical schemes to satisfy (1e) and other constraints. Since NLP solvers are designed to satisfy nonlinear constraints, even restoring from infeasible initial guesses, the MPCC approach demonstrates greater robustness with limited data, resulting in higher success rates.

On the other hand, the MICP formulation demonstrates significantly faster solving speeds when provided with sufficient data. This advantage likely stems from the efficiency of modern convex optimization solvers. For MPCC problems solved with IPOPT, even when we relax the complementarity constraint (2) to values between 0.005 and 0.01 and provide quality warm-start solutions, the solver frequently encounters numerical difficulties and shows slow convergence. The complementarity constraint (2) remains inherently challenging to satisfy. Notably, for the linear inverted pendulum with contact, directly solving using MICP with Gurobi can achieve global optimality at 300 Hz without any data-driven acceleration (Lin, 2024). In addition, numerical accuracy remains a concern with the MPCC approach. Although we only relax constraint (2) to values between 0.005 and 0.01, less than 1% of the binary values, these small relaxations can still lead to accuracy issues in certain cases.

The amount of training data required for effective data-driven MIBLP acceleration depends on the dimensionality of Θ and the problem complexity, such as the number of constraints. Previous theoretical work by Hauser (2016) has established worst-case bounds for data requirements in parametric optimization problems, showing that the amount of data needed scales exponentially with the dimensionality of Θ . While Hauser (2016) assume globally optimal solutions, they provide useful guidance for understanding scaling behavior when using locally optimal warm-start solutions. In practice, the data requirements can be significantly better than these theoretical bounds, especially for the MPCC formulation since NLP solvers are often capable of converging from poor initial guesses to feasible solutions.

Therefore, despite the benefits of data-driven acceleration for solving optimization problems, different formulations exhibit distinct performance characteristics and dataset requirements. The choice of problem formulation is crucial for maximizing the effectiveness of data-driven optimization methods.