Glass Box / Black Box Optimization

pounce.trf_minimize solves problems where part of the model is an equation and part is a program:

\[ \min_x f(x) \quad \text{s.t.} \quad h(x)=0,\; g(x)\le 0,\; y = d(w) \]

Here \(f, h, g\) are ordinary algebra with exact derivatives — the glass box — while \(d\) is a black box: a CFD solve, a converged unit model, a trained network. Something you can call but cannot hand to an NLP solver as equations. \(w\) and \(y\) are subvectors of \(x\).

This is common in process engineering, where a flowsheet is algebraic except for one unit that needs its own simulator.

Why not just fit a surrogate and optimize that?

Because it does not work, and it fails quietly. Consider

\[ \min\; x_1^2 + x_2^2 \quad \text{s.t.}\quad x_2 = x_1^3 + x_1^2 + 1 \]

whose solution is \((0, 1)\) with \(f = 1\). Replace the constraint with a linear surrogate \(x_2 = x_1 + b\), fit \(b\) so the surrogate matches the truth model at the current point, optimize, refit, repeat. Started exactly at the optimum, this iteration walks away and converges to \((-1, 1)\), where \(f = 2\) — a local maximum of the real problem (Biegler 2024, Fig. 2a).

The reason is that matching values gives feasibility but says nothing about optimality, which is a statement about gradients. The trust-region filter method fixes this by (a) correcting the surrogate to match the truth model in both value and slope, and (b) confining each step to a region where that correction is still valid.

Quick start

No derivatives are required. jac (the glass-box gradient) and truth_jac (the black-box Jacobian) are both optional; leave either out and it is approximated by finite differences. The two examples below solve the same problem both ways.

Without derivatives (finite differences)

import numpy as np
import pounce

# Minimize (z-1)^2 + x^2 subject to z = sin(x), treating sin as a black box.
# The variable vector is v = [x, z]: v[0] is the black-box input w,
# v[1] is its output y.
res = pounce.trf_minimize(
    fun=lambda v: (v[1] - 1.0) ** 2 + v[0] ** 2,
    x0=[0.5, 0.0],
    truth_model=lambda w: np.sin(w),
    w_index=[0],
    y_index=[1],
    basis="zero",
)
print(res.x, res.fun, res.n_truth_evals)
# [0.4787 0.4606] 0.52008 13

All you supply is the objective and something you can call. Finite differences are used in two places:

• The truth model. With no truth_jac, \(J_d(w_k)\) is finite-differenced at the sampling radius, which costs up to \(n_w\) extra truth-model calls per iteration. That is why this run needs 13 evaluations instead of 8.
• The glass box. With no jac, each subproblem goes to pounce.minimize without a gradient, so pounce.minimize finite-differences fun (and any constraints given without a jac). This costs only cheap glass-box calls, not truth-model calls. pounce.minimize emits a one-time UserWarning saying so; that warning is advice, not an error.

basis="zero" (the default, written out here for clarity) makes the surrogate a plain linearization of the truth model, so it needs no extra sampling of its own. That makes it the right partner for a derivative-free call: the only per-iteration truth-model cost is the base point plus the finite-difference stencil.

With derivatives

res = pounce.trf_minimize(
    fun=lambda v: (v[1] - 1.0) ** 2 + v[0] ** 2,
    x0=[0.5, 0.0],
    truth_model=lambda w: np.sin(w),
    w_index=[0],
    y_index=[1],
    basis="zero",
    jac=lambda v: np.array([2 * v[0], 2 * (v[1] - 1.0)]),
    truth_jac=lambda w: np.cos(w).reshape(1, 1),
)
print(res.x, res.fun, res.n_truth_evals)
# [0.4787 0.4606] 0.52008 8

Both runs take 7 iterations and reach the same point. The analytic run needs 8 truth-model calls against 13. You can also mix the two: supply jac for the algebra you wrote yourself and omit truth_jac for a simulator that has no sensitivities. See Supply truth_jac if you possibly can for when supplying it pays off.

The method still needs the black box to be smooth enough to have a derivative. Finite differences estimate the Jacobian; they do not remove the need for one. For a noisy truth model, see Limitations.

You do not write the relationship y = d(w) into constraints yourself — w_index and y_index tell the method which variables they are, and it installs the surrogate constraint into each subproblem.

A worked, runnable version of everything below is in python/notebooks/29_trust_region_filter.ipynb: the failure mode plotted against the objective contours, the ZOC/FOC identities demonstrated on a deliberately bad basis, a basis cost comparison, and the Eason & Biegler benchmark with its convergence traces.

How it works

At each iteration the method:

1. Builds a corrected surrogate at the current point \(w_k\): \[ r_k(w) = \bar r(w) + \big(d(w_k) - \bar r(w_k)\big) + \big(J_d(w_k) - J_{\bar r}(w_k)\big)(w - w_k) \] The last two terms are the zero- and first-order corrections (ZOC/FOC). They guarantee \(r_k(w_k) = d(w_k)\) and \(J_{r_k}(w_k) = J_d(w_k)\) for any differentiable basis \(\bar r\).

2. Solves an ordinary NLP with \(d\) replaced by \(r_k\), and the trust region imposed as bounds on the decision variables: \(\max(l, u_k - \Delta) \le u \le \min(b, u_k + \Delta)\).

3. Evaluates the truth model at the trial point and measures the mismatch \(\theta = \lVert y - d(w) \rVert\).

4. Accepts or rejects via a filter on \((\theta, f)\) — a Pareto front of trade-offs rather than a penalty function with a parameter to tune.

Because ZOC/FOC makes the surrogate κ-fully linear automatically, the compatibility check and criticality phase of the original 2016 algorithm are unnecessary and \(\Delta\) need not shrink to zero.

Choosing a basis

The basis \(\bar r\) affects efficiency, not correctness — the corrections hold regardless of how good it is.

You can also pass any object implementing the Basis protocol (fit, predict, jacobian) — a low-fidelity physical model, a symbolic regression fit, a Gaussian process. No adapter or registration is needed.

There is no "linear", and that is not an oversight

Only curvature survives the correction. Write the ZOC/FOC formula with an affine basis \(\bar r(w) = a + B(w - w_{\text{ref}})\). Its Jacobian is \(B\) everywhere, so the basis-dependent part is

\[ \bar r(w) - \bar r(w_k) - B(w - w_k) = B(w-w_k) - B(w-w_k) = 0 \]

leaving \(r_k(w) = d(w_k) + J_d(w_k)(w - w_k)\) — exactly the "zero" result. An affine basis is therefore provably incapable of changing anything, while costing \(n_w + 1\) truth-model calls per iteration to compute it. Before adding any basis, check whether it is affine in \(w\); if it is, it cannot help.

Two findings from the literature worth internalizing before reaching for something fancy:

• Fit quality does not predict optimization performance. Pedrozo et al. (2025) benchmarked five surrogate families on a CO₂ pooling problem. Radial basis functions had the best R² of any model and needed 8 TRF iterations; Kriging needed 2; global polynomials were worst at both.
• Simple often wins outright. On Williams-Otto, Eason & Biegler (2016) found linear interpolation beat Kriging by 91 truth-model calls to 3141.

Fit the basis once and freeze it

By default a string basis is re-fitted from fresh truth-model samples every iteration. That is the wrong trade when the truth model is expensive, and it is not how the literature uses surrogates: Pedrozo et al. fit an ALAMO model once from a designed dataset and then let ZOC/FOC re-anchor it at each new point.

Pass a pre-fitted object and trf_minimize freezes it — fit is never called, no per-iteration sampling happens, and each iteration costs one truth-model evaluation plus a gradient. On the sin example:

Freezing keeps the quadratic’s three-iteration convergence but drops the in-loop cost from 10 calls to 4.

Freezing beats refitting the same basis; it does not automatically beat "zero". At \(n_w = 3\) on a similar problem the numbers come out:

Freezing cuts the quadratic’s cost by more than half (41 → 17), but the free "zero" basis still wins outright. The upfront design is \((n_w{+}1)(n_w{+}2)/2\) calls — 10 at \(n_w=3\), 66 at \(n_w=10\) — paid whether or not the curvature turns out to help. It amortizes when the run is long or each truth call is genuinely expensive, and not otherwise.

Start with "zero". Reach for a frozen curved basis when you can see the iteration count is the bottleneck and you have samples to spare.

from pounce.trf import QuadraticBasis, quadratic_design

design = quadratic_design(w0, 0.05)
basis = QuadraticBasis().fit(design, np.vstack([truth_model(w) for w in design]))

res = pounce.trf_minimize(..., basis=basis)   # frozen automatically

The default is auto: string bases refit, user-supplied objects freeze — if you fitted a model yourself it will not be clobbered. Override either way with refit_basis=True|False.

This is also sound rather than merely cheap. Because ZOC/FOC forces \(r_k(w_k) = d(w_k)\) and \(J_{r_k}(w_k) = J_d(w_k)\) at every new base point, a basis fitted somewhere else entirely still converges to the truth model’s solution — it only contributes curvature.

Supply truth_jac if you possibly can

It is the single highest-value option here. It removes \(n_w\) truth-model calls per iteration and makes the surrogate exactly first-order accurate rather than accurate to the finite-difference step. Many simulators expose it — COMSOL’s sensitivity module, Aspen’s equation-oriented mode — and the ZOC/FOC variant is designed around the assumption that it is available.

Without it, the method finite-differences at the sampling radius \(\sigma_k\). That radius is deliberately not machine epsilon: the right perturbation is a property of the truth model. Eason & Biegler had to inflate it to \(\min(0.1, 0.8\Delta)\) for a boiler model “to compensate for greater numerical noise in the model outputs”, against \(10^{-5}\) for smooth steam tables.

Two radii, not one

trust_radius (\(\Delta\)) bounds the step. sampling_radius (\(\sigma \le \Delta\)) bounds where the surrogate is fit. Eason & Biegler (2018) introduced this separation and reported it more than doubled the number of problems their test set could solve: with one radius doing both jobs, the algorithm is forced into tiny steps near the solution purely to keep the model accurate.

Convergence

trf_minimize reports success when both:

• \(\theta \le\) feasibility_tol — the surrogate agrees with the truth model, so the point is feasible for the real problem; and
• the step in \(w\) is below criticality_tol — so the FOC gradient match is still valid at the solution.

The second is not optional. ZOC/FOC pins \(J_{r_k} = J_d\) only at the base point, so a large accepted step leaves the subproblem solving KKT conditions with a stale Jacobian. Testing feasibility alone will report convergence at points whose true gradient is of order one.

Limitations

Noise. The truth model must be deterministic and smooth. Eason & Biegler assume noise is negligible and list rigorous noise handling as open. For optimization against physical measurements, use a noise-aware method such as Bayesian optimization.

Local. Converges to a local KKT point. Pair it with find_minima for a multistart sweep.

Restoration is partial. The published algorithm calls a restoration procedure in two situations; this implementation handles one and detects the other.

Incompatible subproblem. The glass-box constraints may have no solution inside the current trust region. What trf_minimize does depends on the iterate:

• Infeasible iterate (θ > feasibility_tol). It takes a restoration step: the point inside the same box that minimizes the surrogate mismatch ½‖y − r(w)‖² subject to the glass-box constraints (always feasible — the current point satisfies them). The current pair enters the filter first, so the restored point must improve on it; the step is logged as restoration and the radius follows the θ-ratio test. Without it, a linear surrogate whose linearized constraints are nearly parallel can cycle: the step that reaches their intersection is rejected by the filter, the radius contracts to half of it, the subproblem is incompatible, and expanding the radius re-poses exactly the rejected subproblem (gh#1000, basis="zero" on a sintering model).
• Feasible iterate (θ ≈ 0). There is nothing to restore — the radius is simply too small for the constraints, common on the first iteration. Contracting would be exactly backwards, so trf_minimize expands the radius and retries, logging the iteration as incompatible. On Eason’s example 1 this fires three times before the first real step:

   0  incompatible (Infeasible_Problem_Detected); Delta -> 2.000e-01
   1  incompatible (Infeasible_Problem_Detected); Delta -> 4.000e-01
   2  incompatible (Infeasible_Problem_Detected); Delta -> 8.000e-01
   3  f= 0.2748077559  theta=4.513e-02  Delta=1.427e+00  |s|=7.134e-01  f-step

If the radius reaches trust_radius_max and the subproblem is still infeasible, the constraints are infeasible for reasons unrelated to the trust region, and you get a clear error rather than a wrong answer.

Blocked filter. When the filter rejects every candidate, the iteration can settle into an f-step / θ-step / rejected limit cycle in which θ creeps down but the step length never shrinks. Real restoration would find a point acceptable to the filter; without it, trf_minimize detects the stall and returns success=False with a Stalled: message rather than silently burning its iteration budget. If you hit it, try a larger trust_radius, a richer basis, or a looser feasibility_tol.

References

• Eason, J.P. & Biegler, L.T. A trust region filter method for glass box/black box optimization. AIChE J. 62, 3124–3136 (2016).
• Eason, J.P. & Biegler, L.T. Advanced trust region optimization strategies for glass box/black box models. AIChE J. 64, 3934–3943 (2018).
• Yoshio, N. & Biegler, L.T. Demand-based optimization of a chlorobenzene process… AIChE J. 67, e17054 (2021).
• Biegler, L.T. The trust region filter strategy. Digital Chemical Engineering 13, 100197 (2024).
• Pedrozo, H.A. et al. Surrogate model optimization: a comparison case study with pooling problems of CO₂ point sources. Comput. Chem. Eng. 200, 109199 (2025).