Main Free Edition Commercial Edition Docs FAQ Forum About Us

Rocket: derivative-free MINLP example

This page demonstrates an example of using ALGLIB's MIVNS algorithm to solve a toy derivative-free simulation-based MINLP: designing a low cost rocket capable of delivering a tiny payload to low Earth orbit. The objective is to minimize rocket cost subject to the burnout velocity being at least 7.8 km/s. The solution involves multiple flight simulations with various designs and fuel loads.

Alongside MIVNS, which targets derivative-free problems, ALGLIB also provides BBSYNC, an NLP-based branch-and-bound method for models with analytic derivatives and relaxable integer variables. This page focuses on MIVNS only. BBSYNC and ALGLIB's broader support for MINLP are covered in more detail in the primary MINLP article.

Free (fully functional, no limits on problem size) and commercial cross-platform implementations are available for C++, C#, Java, Python, and Delphi.

Introduction

This toy MINLP designs an optimal rocket capable of launching a tiny payload to low Earth orbit. The rocket may use up to three stages: a large 9-engine stage, a smaller 3-engine stage, and a tiny single-engine upper stage.

The objective is to minimize rocket cost (stages + fuel) subject to the burnout velocity being at least 7.8 km/s. To showcase how MIVNS handles linear constraints on integer variables, we also limit the design to at most two stages (out of three available).

For plausibility we use Rocket Lab's Rutherford engine as a rough reference, and we account for Earth's gravity losses. However, this is a toy problem, not a detailed engineering study, so we make several simplifications:

• The rocket flies strictly vertically; we ignore the gravity turn needed for LEO.
• Reaching a vertical velocity of at least 7.8 km/s at burnout counts as success.
• We include fuel and engine mass but ignore auxiliary structures such as fuel tanks and supports. In our model, the rocket is held together by willpower.
• Differential equations are integrated using a very basic Euler method.

BBSYNC and other MINLP solvers are provided by the minlpsolvers subpackage of the MINLP package (the link opens the Reference Manual for the programming language of your choice). For readability, the article shows code snippets in C++. Complete code in all programming languages supported by ALGLIB is available in the ALGLIB Reference Manual entry for this example: minlp_expensive_dfo.

Solver instance initialization

int n = 6;
minlpsolverstate solver;
real_1d_array x0 = "[1,0,0,0,0,0]";
minlpsolvercreate(n, x0, solver);

Variables and scales

real_1d_array bndl = "[0,0,0,0,0,0]";
real_1d_array bndu = "[1,1,1,inf,inf,inf]";
real_1d_array s = "[1,1,1,1000,1000,1000]";
minlpsolversetintkth(solver, 0);
minlpsolversetintkth(solver, 1);
minlpsolversetintkth(solver, 2);
minlpsolversetbc(solver, bndl, bndu);
minlpsolversetscale(solver, s);

The problem uses six decision variables: three binary variables (include stage 1/2/3 or not) and three continuous ones (propellant mass per included stage). Integer variables are truly non-relaxable: you can not run the simulator with a stage partially included. Hence, BB-type solvers can not be used.

We set the following per-variable constraints:

• [0,1] plus integrality for x0, x1, x2
• non-negativity for x3, x4, x5

Derivative-free and surrogate model methods are highly sensitive to variable scaling. If variables differ by 10³ or more in magnitude, distance metrics and stopping criteria become skewed: a step that is large for one variable will be tiny for another.

Accordingly, we provide a vector of per-variable scales. Integer variables have unit scale; fuels are measured in kilograms but the optimizer is told that typical magnitudes are around thousands. The solver rescales the problem internally, transparently to the user.

Encoding structure with linear constraints

Fixing irrelevant variables

real_1d_array stage1bnd = "[15000,0,0,-1,0,0]";
real_1d_array stage2bnd = "[0,5000,0,0,-1,0]";
real_1d_array stage3bnd = "[0,0,1500,0,0,-1]";
minlpsolveraddlc2dense(solver, stage1bnd, 0.0, fp_posinf);
minlpsolveraddlc2dense(solver, stage2bnd, 0.0, fp_posinf);
minlpsolveraddlc2dense(solver, stage3bnd, 0.0, fp_posinf);

We have already declared variables #3, #4, and #5 (fuel masses) as non-negative. Their upper bounds (max tank capacities) depend linearly on on the corresponding stage-selection variables #0, #1, and #2:

x[3] ≤ 15000*x[0]
x[4] ≤ 5000*x[1]
x[5] ≤ 1500*x[2]

When a stage is not used, its fuel capacity is forced to zero. One could achieve a similar effect by charging for unused fuel in the objective so it converges to zero. However, encoding these relationships as explicit linear constraints (whose structure the solver understands) helps to prune irrelevant regions of continuous subspaces quickly before optimization begins.

Furthermore, fixing irrelevant variables (usually at lower bounds) greatly helps the solver to converge. Surrogate-based derivative-free solvers are particularly easy to confuse with irrelevant variables.

Forbidden designs

real_1d_array justtwostages = "[1,1,1,0,0,0]";
minlpsolveraddlc2dense(solver, justtwostages, fp_neginf, 2);

Adding a linear constraint that allows at most two stages prunes irrelevant points of the integer grid. The solver won't even try to simulate a design with three stages.

x[0]+x[1]+x[2] ≤ 2

Similar constraints involving integer variables can enforce common structures: one-hot encoded groups of binary variables, cardinality limits, permutation-like patterns, etc.

Objective and nonlinear constraints

Callback code

void minlp_rocket_fvec(const alglib::real_1d_array &x, alglib::real_1d_array &fi, void *ptr)
{
 //
 // This callback represents a toy MINLP problem of designing an optimal rocket capable
 // of launching a tiny payload to low Earth orbit. The rocket can include up to three
 // stages: a big 9-engine stage, a small 3-engine stage and a tiny single-engine one.
 //
 // The decision variables are:
 // * 3 indicator 0/1 variables for each of three stages (whether the stage is included
 // into the rocket or not)
 // * 3 variables controlling amount of fuel loaded into each stage.
 //
 // The problem is to minimize rocket cost (stages + fuel) subject to velocity at the
 // burnout being at least 7.8 km/s.
 //
 // The objective (rocket price) is straighforward to determine. The velocity at the burnout
 // (the nonlinear constraint) is determined by integrating differential equations describing
 // the rocket in order to minimize its velocity at the burnout. Thus, the entire task is an
 // interesting example of minimizing a bilinear objective subject to expensive simulation-based
 // constraint having no analytic derivatives.
 //
 const double g0 = 9.80665; // m/s^2
 const double payload = 120.0; // kg payload
 const int nstages = 3;
 double r_thrust = 24000.0; // engine thrust
 double r_price = 200000.0; // engine price
 const double thrust[3] = { 9*r_thrust, 3*r_thrust, r_thrust }; // N
 const double Isp[3] = { 311.0, 321.0, 331.0 }; // s (vacuum)
 const double engine_mass[3] = { 9*35.0, 3*35.0, 35.0 }; // kg (fixed, jettisoned after burn)
 const double engine_cost[3] = { 9*r_price, 3*r_price, r_price }; // USD
 const double fuel_cost_per_kg = 1.2; // USD/kg (LOX+RP1)
 //
 // compute velocity at burnout: integrate each stage burn with simple Euler (no drag)
 //
 double mass = payload;
 for(int i=0;i<nstages;i++)
 if( x[i]>0.5 )
 mass += engine_mass[i] + x[nstages+i];
 double v_burnout = 0;
 for(int i=0;i<nstages;i++)
 if( x[i]>0.5 )
 {
 const double T = thrust[i];
 const double isp = Isp[i];
 const double mdot = T/(isp*g0); // kg/s
 double fuel = x[nstages+i];
 const double burn_time = (fuel>0.0)?(fuel/mdot):0.0;
 double dt = 0.05;
 if(burn_time>0.0)
 {
 double target_steps = 2000.0;
 dt = burn_time/target_steps;
 dt = fmax(dt,0.01);
 dt = fmin(dt,0.10);
 }
 while(fuel>1e-12)
 {
 double step = (fuel/mdot<dt) ? fuel/mdot : dt;
 // Net acceleration (no drag): a = (T - mg)/m
 double a = (T - mass*g0)/mass;
 // Euler step
 v_burnout += a*step;
 // Propellant consumption
 double dm = mdot*step;
 fuel -= dm;
 mass -= dm;
 
 if(step<=0.0) break;
 }
 // Stage separation
 mass = fmax(mass-engine_mass[i],payload);
 }
 //
 // compute price
 //
 double price = 0.0;
 for(int i=0;i<nstages;i++)
 price += x[i]*(engine_cost[i]+fuel_cost_per_kg*x[nstages+i]);
 
 //
 // Output scaled values:
 // fi[0] = price measured in millions
 // fi[1] = velocity measures in km/s
 //
 // Having well-scaled price and velicity greatly helps the solver to converge.
 //
 fi[0] = price/1000000.0;
 fi[1] = v_burnout/1000.0;
}

Objective

This problem has a simple objective: the rocket's cost is a bilinear function of binary and continuous decision variables. However, the solver has to infer this structure from sparsely sampled evaluations

Furthermore, it has to somehow infer how to compare changes in the objective against changes in the nonlinear constraints. Is a one-cent price reduction worth violating the velocity constraint by 1 mm/s? What about a thousand dollars versus 1 m/s? What about allowing a temporary violation at a trial point, to be corrected later?

In gradient-based optimization one can obtain sensitivity information by studying typical gradient magnitudes. In derivative-free settings there is no automatic, reliable analogue. So we help the solver by scaling the objective and nonlinear constraints so their typical magnitudes have roughly unit scale and thus easily comparable.

The price is computed in dollars, but typical rocket prices are measured in millions, so we rescale objective by dividing it by 10⁶.

Nonlinear constraint

//
// Add a nonlinear constraint: burnout velocity must be at least 7.8 km/s.
//
// This line only declares the constraint and its bounds; the actual computation is supplied
// by the callback.
//
minlpsolveraddnlc2(solver, 7.8, fp_posinf);

Our only nonlinear constraint is a long-running black-box numerical simulation with no analytic derivatives. We expect low level of numerical noise and reasonable degree of smoothness in the region of interest. These properties are essential for the surrogate-model optimization MIVNS uses in continuous subspaces.

The same scaling considerations that apply to the objective apply to the constraint as well. The solver assumes that constraint violations should be much smaller than 1 in magnitude. Internally we use SI units, so typical burnout velocities are between 1.000 m/s and 10.000 m/s in magnitude. If we keep the constraint in meters per second, the solver would try to satisfy them up to millimeters per second, wasting time chasing an unnecessarily tight tolerance.

We therefore rescale the constraint by dividing velocities by 1000, yielding the much better scaled constraint v_burnout ≥ 7.8.

Configuring the solver

Code snippet

int batchsize = 1;
int budget = 150;
int maxneighborhood = 0;
minlpsolversetalgomivns(solver, budget, maxneighborhood, batchsize);

MIVNS overview

MIVNS (Mixed-Integer Variable Neighborhood Search) is ALGLIB's derivative-free solver for mixed-integer nonlinear programs with non-relaxable integer variables and expensive objective/constraint evaluations.

It explores the integer grid by expanding and scanning neighborhoods around the current incumbent, while running short, RBF surrogate searches in the associated continuous subspaces. No convexity or derivative information is required; only basic smoothness with respect to the continuous variables.

MIVNS efficiently searches valid designs, allocating evaluations where improvement is most likely, making it a practical choice for low-to-medium accuracy solutions of costly black-box MINLPs like this rocket model. Its key properties are:

• Solves all-integer or mixed-integer problems with box, linear, and nonlinear (in)equalities.
• Respects variable integrality at trial points, respects linear constraints on integral variables when forming neighborhoods.
• Excellent budget discipline: provides good progress with limited evaluations and scales well with larger budgets.
• Well suited to "structural-switch" problems where integer choices rewire the model (e.g., the choice of rocket's stages), and to simulators with high per-call cost.

Computational budget

Budget (function evaluations): a soft cap on total work. MIVNS may exceed it slightly (due to batched requests), but it stops shortly after crossing the limit. How to choose the budget?

A 'neighborhood scan' is the minimum number of function evaluations needed to perform a basic analysis of the immediate neighborhood. For an N-dimensional problem with N_I integer variables and N_F continuous ones there are ≈N_I neighboring nodes, and each node needs ≈N_F evaluations to fit at least linear model of the objective/constraints. Thus, one MIVNS neighborhood scan costs about N_I·N_F=N_I·(N-N_I)=N_F·(N-N_F) function evaluations.

MIVNS does not share information between integer nodes because it assumes that objective landscape may change drastically when jumping between nodes. That's why we have N_I·N_F instead of N_I+N_F as an estimate.

When either N_F or N_I is small, the scan cost is O(N), and good progress is often achievable with 5N...100N objective values. However, when both N_I and N_F are near N/2, a scan needs ≈N² evaluations, leading to worse scaling.

In practice, significant improvement typically comes within 5-10 scans. Difficult problems may need more than 50 scans. For our toy problem we choose budget=150, which corresponds to roughly 15 neighborhood scans.

Stopping conditions

The simplest way to stop MIVNS is to exhaust the computational budget. This is a default behavior if no other condition is set. Alternatively, you can stop after (a) all nodes in the current neighborhood have been analyzed, and (b) the neighborhood size exceeds maxneighborhood.

Practical choices are often proportional to the number of integer variables (e.g., 5N to 10N as a starting point), but the optimal value is highly problem-dependent. If you can experiment, a useful approach is to set the maximum budget possible with no neighborhood size limit; treat the result as a reference point. Then re-run while increasing the neighborhood size until you can reproduce the reference. Increase neighborhood limit several-fold to add stability.

Batch size

MIVNS solver can simultaneously process up to batchsize nodes from the current neighborhood. Depending on parallelism settings, it can parallelize internal computations (neighbor search and surrogate RBF model optimization), objective/constraints evaluations (callback parallelism; needs thread-safe callbacks) or both. If a single node (callback+internals) takes longer than 1ms, at least one of these phases is long enough to pay off multithreading overhead. Large batches may allow you to parallelize even sub-millisecond tasks.

MIVNS generally scales well with the cores count. Up to 2N_int neighbors of a central node are already there, and the solver can expand the neighborhood proactively if some cores are idle.

In our example we used single-threaded execution. In the production code don't forget to compile ALGLIB with parallelism enabled and to tell ALGLIB that it can use parallelism by passing alglib::parallel|alglib::parallelcallbacks to minlpsolveroptimize().

User callbacks must be re-entrant in order to parallelize objective evaluations. Their re-entrancy is not assumed by default; one states it with the alglib::parallelcallbacks flag passed to the optimizer along with the callback. A good design pattern is to avoid writable global state, use thread-local buffers and preallocate temporaries to avoid allocator contention.

Run the solver

minlpsolverreport rep;
real_1d_array xf;
alglib::minlpsolveroptimize(solver, minlp_rocket_fvec);
minlpsolverresults(solver, xf, rep);
printf("%.6f\n", double(rep.f));

After the run, we have rocket cost roughly equal to 810.000 USD. The decision variables vector is [0,1,1,0.000,4982.0,1463.7], which corresponds to the medium and the tiny stages being used together, with both stages being fueled almost up to full capacity (5.000 kg and 1.500 kg correspondingly). The largest stage (9 engines) is not used.

As of 2025, Rocket Lab's launch price is several times higher than 0.8M USD. Perhaps, ignoring fuel tanks, structural supports, atmospheric drag and gravity turn resulted in some underestimation of actual costs.

This article is licensed for personal use only.

ALGLIB® - numerical analysis library, 1999-2025.
ALGLIB is a registered trademark of the ALGLIB Project.
Policies for this site: privacy policy, trademark policy.