Generating OMARS Designs#

OMARS designs (orthogonal minimally aliased response surface designs) are three-level designs that sit between screening designs and full response surface designs. Every main effect is orthogonal to every other main effect and to all second-order terms (the pure quadratics and the two-factor interactions), so all aliasing is confined to the second-order block, where it is kept minimal.

generate_omars() constructs such designs on demand with an integer linear program (ILP). Unlike the minimal conference-foldover member produced by generate_design(design_type="omars") (which is saturated and leaves no error degrees of freedom), the designs from generate_omars are sized to support a full second-order analysis with analyze_omars().

Note

OMARS® is a registered trademark of KU Leuven Research & Development. process-improve uses the name only for the class of designs defined by Núñez Ares and Goos (2020), and is not affiliated with, or endorsed by, KU Leuven, Minitab or Effex.

Note

generate_omars solves its ILP with HiGHS, through scipy.optimize.milp, which is part of the core install, so no extra is needed. The ilp extra that used to supply the CBC solver is deprecated and now installs nothing.

Quick start#

from process_improve.experiments import Factor, generate_omars, analyze_omars

factors = [Factor(name=n, low=-1, high=1) for n in "ABCDE"]

# Smallest foldover OMARS design that still leaves error degrees of freedom.
result = generate_omars(factors)
print(result.metadata["n_runs_selected"], result.metadata["expected_error_df"])
print(result.metadata["omars_verified"])  # True

# The design is ready for the staged OMARS analysis. A simulated response
# stands in for the measurements here.
import numpy as np

design = result.design[result.factor_names]
rng = np.random.default_rng(1)
y = 3 * design["A"] - 2 * design["B"] + 1.5 * design["A"] ** 2 + rng.normal(0, 0.5, len(design))
analysis = analyze_omars(design, y)
print(analysis.success, analysis.active_main_effects)  # True ['A', 'B']

You can pin an exact run size or search a window:

result = generate_omars(factors, n_runs=33)
result = generate_omars(factors, n_runs_range=(31, 41))

The method#

Every design produced here is a foldover \([H; -H; 0]\): a half-design \(H\), its mirror image \(-H\), and a single centre run. The foldover structure makes three of the four OMARS-defining conditions hold automatically, which is what keeps the construction tractable.

For a design coded to \(\{-1, 0, +1\}\):

  • Balance is automatic: a run \(h\) and its mirror \(-h\) cancel, so every main-effect column sums to zero.

  • Main effects clear of the two-factor interactions is automatic: the term \(x_i x_a x_b\) is an odd function, so the contributions from \(h\) and \(-h\) cancel.

  • Main effects clear of the pure quadratics is automatic: \(x_i x_j^2\) is odd in \(x_i\), so those contributions cancel too.

  • Quadratics are estimable because the centre run makes each \(x_i^2\) column take the value 0 at least once (so it is not constant).

The condition that is not automatic is the mutual orthogonality of the main effects. That condition is linear in the binary “include half-run” variables \(s_r\): for every pair of factors \(i < j\),

\[\sum_r \left( x_{r,i}\, x_{r,j} \right) s_r = 0 .\]

Each factor must also reach an outer level in at least one half-run, \(\sum_r |x_{r,i}|\, s_r \ge 1\); a factor left at its middle level throughout has an all-zero column, which is orthogonal to everything but is not a three-level factor. The run count is \(N = 2\sum_r s_r + 1\). The ILP therefore selects a half-design from the \((3^k - 1)/2\) distinct non-mirror three-level runs subject to only \(k(k-1)/2\) equality constraints and \(k\) coverage constraints. Because the coefficients are integers, the constraints are exact. Every selection the solver returns is re-checked exactly, and a floating-point is_omars() re-check guards every accepted design as a sanity check.

The estimability frontier#

A foldover cannot fit a full second-order model at just any run size, and the threshold is higher than simply counting parameters.

Every second-order term is an even function of the factors, so the quadratic and interaction columns of \(H\) and \(-H\) are identical. The even block of the model matrix therefore has at most \(h + 1\) distinct rows, against \(1 + k(k+1)/2\) columns (an intercept, \(k\) quadratics and \(k(k-1)/2\) interactions). The main effects live in the odd block and contribute at most \(k\) more, so for every foldover

\[\mathrm{rank}(X) \le k + \min\!\left(h + 1,\; 1 + \tfrac{k(k+1)}{2}\right),\]

with equality for half-designs in general position. The full second-order model is therefore estimable only from

\[N = k^2 + k + 1\]

runs: 13, 21, 31, 43 and 57 runs for three to seven factors. This exceeds the parameter count \(1 + 2k + k(k-1)/2\) by exactly \(k(k-1)/2\), the number of two-factor interactions, so “more runs than parameters” is not a sufficient test at any factor count. Sizing starts at this frontier, and n_runs below it is refused.

Choosing model="main_quadratic" drops the interactions from the model being sized for, which lowers the frontier to the definitive screening design’s \(2k + 1\) runs.

Choosing the run size and the design#

  • When n_runs is given it is used directly. n_runs - center_runs must be even, and n_runs must reach the estimability frontier described above; a smaller value is rejected rather than silently producing a design that cannot fit the model.

  • Otherwise the solver minimises the run count within a window to find the smallest feasible design that still leaves error degrees of freedom. A feasible size is not always a usable one: when no design found at that size can estimate the model, the search moves up the window one size at a time (omars_search.run_sizes_searched counts the sizes tried). A size that needs more half-runs than the \((3^k - 1)/2\) distinct ones can only be built by repeating half-runs, which the exhaustive search does for three and four factors.

  • Candidate designs are then collected at that run size, in one of two ways (metadata["search_mode"] says which):

    • "exhaustive": when the design class is small enough (three or four factors at moderate sizes), every feasible design is enumerated, so the selection below is exact.

    • "multistart": otherwise the ILP is solved repeatedly with random linear objectives (n_restarts of them). Each objective steers the solver towards a different feasible design. Each of these solves stops after a fixed number of branch-and-bound nodes (solver_options["node_limit"], default 100) and returns the best design found by then, which ends it at the same point on every run, so a fixed random_state reproduces the design. Designs that differ only in run order, factor order or factor signs are counted once. For the "dominance" and "d_efficiency" criteria, each design found is then improved by a local search that swaps one, two or three half-runs at a time while keeping the main effects orthogonal, moving to a swap whenever it raises the D-efficiency.

    A valid OMARS design can still leave the model it was sized for rank-deficient; such designs are set aside before the selection, and each one adds a constraint to later solves that excludes every design sharing its rank deficiency. The winner is chosen by selection_criterion:

    • "dominance" (default) keeps the Pareto front on D-efficiency (higher is better) and the maximum second-order correlation (lower is better), then prefers the smallest, most efficient design.

    • "d_efficiency" maximises the D-efficiency of the model the design is sized for.

    • "min_second_order_correlation" minimises the largest second-order correlation.

    • "a_optimal" minimises the summed coefficient variance \(\operatorname{tr}\left((X^\top X)^{-1}\right)\) of that model.

    Under every criterion, a design with a second-order column that never varies (a term it cannot estimate) ranks last.

  • Optionally, satisfice sets acceptability thresholds that are applied before the dominance/criterion step: a design is kept only if it clears every threshold. Supported keys are "d_efficiency" (a minimum) and "max_second_order_correlation" (a maximum), e.g. satisfice={"d_efficiency": 5.0, "max_second_order_correlation": 0.7}. Together these implement the satisficing-and-dominance multicriteria selection of Nunez Ares and Goos (2020): first discard designs that fail the minimum bars (satisficing), then drop dominated designs and choose from the Pareto front (dominance). This deliberately avoids collapsing several criteria into a single weighted score, which would hide the trade-offs. A ValueError is raised if no enumerated design meets the thresholds.

The returned DesignResult records the provenance and a search report under metadata (family, sparsity, expected_error_df, d_efficiency, max_second_order_correlation, solver_status and an omars_search report). The report counts the ILP solves and their time, the solves a node or time limit stopped, and the rank-deficient designs set aside, and it says whether the run size was proven to be the smallest in the window.

Performance: iterations and timing by factor count#

The table below reports, for the automatic smallest-size search at n_restarts=8, the size of the candidate half-pool, the run size of the design returned, the resulting error degrees of freedom, how the candidates were collected, the number of ILP solves performed (the iteration count: the minimise-size solve, plus a feasibility solve and the restarts on the multistart path) and the wall-clock time of the whole call. The run size is the estimability frontier \(k^2 + k + 1\) described above, and the error degrees of freedom follow as \(N - (1 + 2k + k(k-1)/2)\); both are exact. Times were measured on an x86_64 machine with CPython 3.11 and SciPy 1.17.1 (HiGHS 1.12), single-threaded, with the default solver_options; they are indicative, vary by machine, and scale with n_restarts.

Factors \(k\)

Half-pool size

Runs \(N\)

Error df

Search

ILP solves

Time (s)

3

13

13

3

exhaustive

1

< 0.1

4

40

21

6

exhaustive

1

11

5

121

31

10

multistart

10

2.8

6

364

43

15

multistart

10

13

7

1093

57

21

multistart

10

41

On the exhaustive path the time is the enumeration (about 250,000 four-factor designs), not the solver. On the multistart path the iteration count is set by n_restarts and can be lower when the search stops early because new solves only repeat designs already found. Each restart is bounded by solver_options["node_limit"] (default 100 branch-and-bound nodes), so its cost grows with the half-pool size \((3^k - 1)/2\) and the run size the frontier demands, but not without limit. The time beyond the solves is the local search. At the default n_restarts=50 the five-factor search takes about 20 s and the six-factor search about 70 s.

generate_design(factors, design_type="omars_ilp", budget=N) does not spend a fixed number of restarts: it stops at a run size once eight solves in a row have returned an estimable design that, after the local search, does not improve the Pareto front on D-efficiency and the maximum second-order correlation, with 50 restarts as the ceiling. For six factors and budget=17 that is about 20 solves, under 20 s.

solver_options["time_limit"] (default 60 s per solve) is a safety cap, not the budget: if it stops a solve, the result depends on the machine’s speed, and metadata["omars_search"].time_limited_solves counts those solves. Beyond seven factors the half-pool triples with each factor; pin n_runs or use model="main_quadratic", whose frontier is only \(2k + 1\).

Limitations#

  • Foldover family only. generate_omars builds the (dominant) foldover OMARS family. The rarer non-foldover members from the enumerated catalogue are a documented future extension.

  • Run-count parity. A foldover design has \(2h\) half-runs and mirrors plus its centre runs, so n_runs - center_runs must be even.

  • Second-order aliasing remains. Reaching the estimability frontier makes the full second-order model fittable; it does not make the second-order block orthogonal. The quadratics and interactions stay mutually aliased - that is the “minimally aliased” in OMARS, not a defect - which is why analyze_omars() resolves them in stages rather than fitting them all at once. Ask for more than the frontier run size to buy error degrees of freedom and lower second-order correlations.

References#

  • Nunez Ares, J. and Goos, P. (2020). “Enumeration and multicriteria selection of orthogonal minimally aliased response surface designs.” Technometrics, 62(1):21-36.

  • Nunez Ares, J. and Goos, P. (2019). “An integer linear programming approach to find trend-robust run orders of experimental designs.” Journal of Quality Technology.