--- title: "Getting Started with OptOR" author: "Christian Palmes, Raluca Ilinca Schmitt, Adrian Funke" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true number_sections: false bibliography: references.bib link-citations: true vignette: > %\VignetteIndexEntry{Getting Started with OptOR} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} Sys.setenv(OMP_NUM_THREADS = "2") knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 5 ) ``` # Introduction In pharmaceutical development, a design space as described in ICH Q8(R2) [@ich2009q8] represents combinations of material attributes and process parameters that have been demonstrated to provide assurance of product quality. Such feasible regions may be curved, non-convex, or otherwise difficult to communicate and implement in routine manufacturing. An axis-aligned hyperrectangle provides a particularly simple operating region. It assigns an individual lower and upper operating limit to each material attribute and process parameter, independent of the settings of the other factors, and can therefore be implemented as a conventional set of parameter intervals. `OptOR`, short for *Optimal Operating Regions*, computes large axis-aligned hyperrectangles within multidimensional feasible regions. Two principal cases are supported: 1. a model-agnostic case in which feasibility is represented by a discrete binary array; and 2. a model-based case in which the feasible region is defined by quadratic response functions. First, the model-agnostic case is introduced. Second, the quadratic-response case is addressed by constructing conservative and optimistic grid approximations, deriving a guaranteed feasible and maximal, that is, non-expandable, continuous hyperrectangle, and quantifying the remaining gap to the unknown global optimum. All worked examples use two factors so that the feasible region, grid classification, and resulting operating regions can be visualized directly. The functions and algorithms apply analogously in higher dimensions. ```{r load-package} library(OptOR) ``` # Model-agnostic operating regions The most general use case starts from a \(d\)-dimensional binary array. Entries equal to `1` represent feasible grid points or cells, whereas entries equal to `0` represent infeasible locations. The binary classification may originate from experimental observations, process simulations, mechanistic models, machine-learning predictions, Monte Carlo evaluations, or any other external procedure. The optimization itself is therefore agnostic to the model or method used to establish feasibility. The following small example uses a two-dimensional binary array. ```{r discrete-array} X <- matrix( c( 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0 ), nrow = 3, byrow = TRUE ) X ``` The optimal axis-aligned rectangle in the discrete array is obtained with `optimal_grid_hr()`, short for *optimal grid hyperrectangle*. ```{r discrete-optimization} grid_result <- optimal_grid_hr( X, verbose = FALSE ) grid_result ``` The lower and upper index vectors define the selected rectangle, \([1,2] \times [1,2]\). Its discrete volume is the product of its side lengths and is therefore equal to \(2 \cdot 2 = 4\). ```{r discrete-volume} grid_widths <- grid_result$ul - grid_result$ll + 1L grid_widths prod(grid_widths) ``` This approach is useful whenever the feasible region is available as a classified grid but the underlying response model is unavailable, unsuitable for direct optimization, or intentionally kept separate from `OptOR`. # Operating regions defined by quadratic response functions A more specific case arises in response surface methodology, where critical quality attributes or other responses are represented by quadratic functions \[ q_j(x) = c_j + b_j^\top x + x^\top Q_j x, \qquad j = 1, \ldots, m, \] where \(x \in \mathbb{R}^d\) denotes the vector of factor settings, \(c_j \in \mathbb{R}\) is the intercept, \(b_j \in \mathbb{R}^d\) is the vector of linear coefficients, and \(Q_j \in \mathbb{R}^{d \times d}\) is the symmetric matrix of quadratic and interaction coefficients. Each response may be subject to a lower acceptance limit, an upper acceptance limit, or both. The feasible region is the set of factor combinations for which all response requirements are fulfilled simultaneously. In this setting, the objective is the global maximum-volume hyperrectangular operating region within the specified working region, rather than merely a locally optimal rectangle obtained from a particular starting point. The conservative and optimistic grid constructions introduced below provide a feasible inner solution and a global upper volume bound, respectively. Quadratic response models are commonly fitted to data from designed experiments. Their simultaneous response constraints may define a curved and potentially non-convex design space. A hyperrectangular operating region must therefore be validated over its complete continuous extent rather than only at its center or at a finite number of selected points. The following example defines the two-dimensional unit disk \[ x_1^2 + x_2^2 \leq 1. \] ```{r define-fct} disk_fct <- list( c = 0, b = c(0, 0), Q = diag(2), lim_ul = 1 ) fcts <- list(disk_fct) ``` # Grid classification `calc_X()` constructs and classifies a regular grid over a \(d\)-dimensional working region. The arguments `rg_ll` and `rg_ul` contain the lower and upper bounds of the factors, \[ \mathrm{rg\_ll} = (\ell_1, \ldots, \ell_d)^\top, \qquad \mathrm{rg\_ul} = (u_1, \ldots, u_d)^\top, \] and therefore define the axis-aligned working region, commonly referred to as the experimental region in design of experiments, \[ \mathcal{R} = [\ell_1, u_1] \times \cdots \times [\ell_d, u_d]. \] The argument `n` specifies the number of grid cells in each coordinate direction. Along factor \(k\), each cell has width \[ \Delta_k = \frac{u_k - \ell_k}{n}. \] The complete grid contains \(n^d\) cells. In point mode, each cell is represented by its midpoint, located halfway between its boundaries in every coordinate direction. The cell is classified as feasible if all response requirements are satisfied at this point. Point mode is useful for exploratory calculations, but feasibility at the midpoint does not guarantee feasibility throughout the entire cell. ```{r classify-grid} grid_point <- calc_X( fcts = fcts, n = 8, rg_ll = c(-1, -1), rg_ul = c(1, 1), gmode = "point" ) grid_point$X ``` In this example, `rg_ll = c(-1, -1)` and `rg_ul = c(1, 1)` define the two-dimensional working region \([-1, 1] \times [-1, 1]\). Each factor range is divided into \(8\) cells of width \(0.25\), resulting in an \(8 \times 8\) classification array. The displayed matrix contains the pointwise feasibility classification of the corresponding cell midpoints. The returned object also contains metadata describing the grid and the response functions. ```{r inspect-grid-result} str(grid_point, max.level = 1) ``` # Conservative grid construction To construct a continuously feasible operating region, the working region is first partitioned into regular grid cells. Each cell is then classified by examining the extrema of all quadratic response functions over the complete cell. For an upper response limit, the maximum response value within the cell must not exceed the specified limit. For a lower response limit, the minimum response value within the cell must not fall below the specified limit. The classification therefore accounts for the behavior of the response functions between grid points rather than evaluating only a single representative point. ```{r conservative-grid} grid_conservative <- calc_X( fcts = fcts, n = 8, rg_ll = c(-1, -1), rg_ul = c( 1, 1), gmode = "conservative" ) grid_conservative$X ``` Each cell is assigned one of three values: * `0` if every point in the cell is infeasible, * `1` if every point in the cell is feasible, and * `2` if the cell contains both feasible and infeasible points. The three-level cell classification itself is identical for `gmode = "conservative"` and `gmode = "optimistic"`. The selected mode determines how cells classified as `2` are interpreted when the maximum grid-aligned hyperrectangle is computed. In conservative mode, these cells are treated as infeasible and therefore handled in the same way as cells classified as `0`. In optimistic mode, they are treated as feasible and assigned the value `1`. As a result, only cells that are feasible throughout can be included in the conservative grid-aligned hyperrectangle. The resulting union of cells is guaranteed to lie within the continuous feasible region, and every hyperrectangle composed entirely of these cells is continuously feasible. The conservative grid therefore defines an inner continuous approximation to the feasible region, rather than merely providing a pointwise numerical discretization. # Continuous refinement by greedy expansion `optimal_cont_hr()`, short for *optimal continuous hyperrectangle*, first determines an optimal hyperrectangle in the conservative binary grid and then expands its boundaries in the continuous domain. ```{r continuous-optimization} cont_result <- optimal_cont_hr( fcts = fcts, n = 30, rg_ll = c(-1, -1), rg_ul = c( 1, 1), gmode = "conservative", verbose = FALSE ) cont_result ``` The conservative grid solution provides a guaranteed feasible starting rectangle. During the subsequent greedy expansion, lower and upper boundaries are moved outward while feasibility over the complete candidate hyperrectangle is repeatedly verified. Expansion directions that still permit an increase are retained. Directions that prevent further enlargement are successively removed. The final hyperrectangle is therefore not restricted to the original grid boundaries. ```{r continuous-volume} cont_widths <- cont_result$ul - cont_result$ll V_cons <- prod(cont_widths) cont_widths V_cons ``` Upon termination, the returned continuous hyperrectangle is maximal, that is, non-expandable: none of its boundaries can be moved further outward in any direction without violating at least one response constraint. A maximal hyperrectangle should not be confused with a global maximum-volume hyperrectangle. Another feasible hyperrectangle with a different location or aspect ratio may, in principle, have a larger volume. # Verifying continuous feasibility `find_extrema()` determines the global minimum and maximum of a quadratic function over a hyperrectangle. It can therefore be used to verify that the returned operating region satisfies a response constraint over the complete rectangle. ```{r verify-region} verification <- find_extrema( fct = disk_fct, hr_ll = cont_result$ll, hr_ul = cont_result$ul ) verification ``` For the unit-disk example, the maximum response over the complete rectangle must not exceed the upper response limit. ```{r verify-upper-limit} verification$max <= disk_fct$lim_ul + 1e-10 ``` # Optimistic grid construction The optimistic classification provides a complementary outer approximation. With `gmode = "optimistic"`, a cell is excluded only when the complete cell can be shown to be infeasible. Cells that may contain at least one feasible point are retained. Consider a feasible continuous hyperrectangle attaining the global maximum volume. Every grid cell intersecting this hyperrectangle contains at least one feasible point and therefore cannot be classified as `0`. Such cells are classified as either `1` or `2`, both of which are treated as feasible in optimistic mode. Hence, every global volume-maximizing feasible continuous hyperrectangle is covered by cells treated as feasible in the optimistic grid. It follows that the volume of the maximum grid-aligned hyperrectangle in the optimistic grid is an upper bound on the global maximum continuous volume. This does not imply that the identified optimistic grid-aligned hyperrectangle itself necessarily contains a global optimal continuous solution. The corresponding optimistic grid-aligned hyperrectangle and its upper volume bound are computed with `optimal_cont_hr()` using `gmode = "optimistic"`. ```{r optimistic-bound} optimistic_result <- optimal_cont_hr( fcts = fcts, n = 30, rg_ll = c(-1, -1), rg_ul = c( 1, 1), gmode = "optimistic", verbose = FALSE ) optimistic_widths <- optimistic_result$ul - optimistic_result$ll V_optimistic <- prod(optimistic_widths) optimistic_widths V_optimistic ``` The returned optimistic hyperrectangle must not subsequently be reduced or refined inward. Such a modification could remove parts of the outer approximation and would therefore invalidate its interpretation as providing a global upper bound on the maximum volume. The conservative solution and the optimistic outer hyperrectangle bracket the global optimum: \[ V_{\mathrm{cons}} \leq V_{\mathrm{continuous}}^* \leq V_{\mathrm{optimistic}}, \] where \(V_{\mathrm{cons}}\) is the volume of the guaranteed feasible continuous solution, \(V_{\mathrm{continuous}}^*\) is the global maximum continuous volume, and \(V_{\mathrm{optimistic}}\) is its global upper bound. As the grid resolution increases, the conservative and optimistic approximations are expected to become tighter under suitable regularity conditions. # A dimension-adjusted approximation metric The direct volume ratio \[ \frac{V_{\mathrm{cons}}}{V_{\mathrm{optimistic}}} \] becomes increasingly difficult to interpret as the dimension \(d\) increases, because differences between side lengths accumulate multiplicatively in the volume. A dimension-adjusted metric is therefore \[ R = \left( \frac{V_{\mathrm{cons}}} {V_{\mathrm{optimistic}}} \right)^{1/d}. \] ```{r approximation-metric} d <- 2L R <- (V_cons / V_optimistic)^(1 / d) R ``` Since the volume of a \(d\)-dimensional hyperrectangle is the product of its side lengths, \(V^{1/d}\) is the geometric mean of those side lengths. Consequently, \[ R = \frac{ V_{\mathrm{cons}}^{1/d} }{ V_{\mathrm{optimistic}}^{1/d} } \] is the ratio of the geometric-mean side lengths of the conservative and optimistic hyperrectangles. Because \(V_{\mathrm{optimistic}}\) is an upper bound on the unknown global maximum, \(R\) is a conservative lower bound on the fraction of the global optimal geometric-mean side length attained by the reported solution. For example, a value of \(R = 0.95\) means that the geometric-mean side length of the conservative solution is guaranteed to be at least \(95\%\) of the corresponding global optimum as bounded by the optimistic approximation. The following visual example compares this guaranteed ratio with the corresponding ratio to an analytically known true optimum. # Visualizing the inner and outer approximations The following two-dimensional example summarizes the complete construction for \[ x_1^2+x_2^2 \leq 1 \] over the coded working region \([-1,1]^2\). A deliberately coarse grid is used so that the distinction between conservative and optimistic cell classifications remains visible. The blue rectangle is the continuously valid conservative solution, the orange rectangle is the optimistic outer bound, and the red rectangle is the analytically known global optimum for this particular example. For the conservative \(8 \times 8\) grid, the chosen global optimal discrete rectangle is given by the index ranges \([3,6] \times [2,7]\). The subsequent continuous expansion enlarges the rectangle in the \(x_1\) direction. The resulting displacement of the blue rectangle away from the underlying grid-cell boundaries is visible in the figure and illustrates the refinement performed after the global optimal conservative grid rectangle has been identified. The extensive conversion of the classified arrays into plotting data is performed in a hidden chunk. The methodologically relevant calls to `calc_X()` and `optimal_cont_hr()` have already been shown above, while the figure remains fully reproducible when the vignette is rendered. ```{r prepare-circle-figure, include=FALSE} n_plot <- 8L plot_rg_ll <- c(-1, -1) plot_rg_ul <- c( 1, 1) plot_fct <- list( c = 0, b = c(0, 0), Q = diag(2), lim_ul = 1 ) plot_fcts <- list(plot_fct) plot_grid_conservative <- calc_X( fcts = plot_fcts, n = n_plot, rg_ll = plot_rg_ll, rg_ul = plot_rg_ul, gmode = "conservative" ) plot_grid_optimistic <- calc_X( fcts = plot_fcts, n = n_plot, rg_ll = plot_rg_ll, rg_ul = plot_rg_ul, gmode = "optimistic" ) plot_hr_conservative <- optimal_cont_hr( fcts = plot_fcts, n = n_plot, rg_ll = plot_rg_ll, rg_ul = plot_rg_ul, gmode = "conservative", verbose = FALSE ) # The conservative classification may contain boundary cells coded as 2. # Only cells coded as 1 are guaranteed feasible and may enter the discrete # conservative optimization. plot_X_conservative_binary <- array( as.integer(plot_grid_conservative$X == 1L), dim = dim(plot_grid_conservative$X) ) stopifnot(all(plot_X_conservative_binary %in% c(0L, 1L))) plot_grid_hr_conservative <- optimal_grid_hr( plot_X_conservative_binary, verbose = FALSE ) plot_hr_optimistic <- optimal_cont_hr( fcts = plot_fcts, n = n_plot, rg_ll = plot_rg_ll, rg_ul = plot_rg_ul, gmode = "optimistic", verbose = FALSE ) true_half_width <- 1 / sqrt(2) plot_hr_true <- list( ll = rep(-true_half_width, 2), ul = rep(true_half_width, 2) ) plot_cell_width <- (plot_rg_ul - plot_rg_ll) / n_plot plot_x_centres <- plot_rg_ll[1] + (seq_len(n_plot) - 0.5) * plot_cell_width[1] plot_y_centres <- plot_rg_ll[2] + (seq_len(n_plot) - 0.5) * plot_cell_width[2] plot_cells <- expand.grid( i = seq_len(n_plot), j = seq_len(n_plot) ) plot_cells$x <- plot_x_centres[plot_cells$i] plot_cells$y <- plot_y_centres[plot_cells$j] plot_cells$conservative_code <- mapply( function(i, j) plot_grid_conservative$X[i, j], plot_cells$i, plot_cells$j ) plot_cells$optimistic_code <- mapply( function(i, j) plot_grid_optimistic$X[i, j], plot_cells$i, plot_cells$j ) plot_cells$classification <- "Outside optimistic approximation" plot_cells$classification[plot_cells$optimistic_code != 0] <- "Optimistic approximation" plot_cells$classification[plot_cells$conservative_code == 1] <- "Conservative" plot_cells$classification <- factor( plot_cells$classification, levels = c( "Outside optimistic approximation", "Optimistic approximation", "Conservative" ) ) plot_theta <- seq(0, 2 * pi, length.out = 1000) plot_circle <- data.frame( x = cos(plot_theta), y = sin(plot_theta) ) plot_rectangles <- data.frame( region = factor( c("Conservative", "True optimal", "Optimistic"), levels = c("Conservative", "True optimal", "Optimistic") ), xmin = c( plot_hr_conservative$ll[1], plot_hr_true$ll[1], plot_hr_optimistic$ll[1] ), xmax = c( plot_hr_conservative$ul[1], plot_hr_true$ul[1], plot_hr_optimistic$ul[1] ), ymin = c( plot_hr_conservative$ll[2], plot_hr_true$ll[2], plot_hr_optimistic$ll[2] ), ymax = c( plot_hr_conservative$ul[2], plot_hr_true$ul[2], plot_hr_optimistic$ul[2] ) ) circle_plot <- ggplot2::ggplot() + ggplot2::geom_tile( data = plot_cells, ggplot2::aes(x = x, y = y, fill = classification), width = plot_cell_width[1], height = plot_cell_width[2], colour = "black", linewidth = 0.25 ) + ggplot2::geom_path( data = plot_circle, ggplot2::aes(x = x, y = y), linewidth = 1.2, colour = "black" ) + ggplot2::geom_rect( data = plot_rectangles, ggplot2::aes( xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax, colour = region ), fill = NA, linewidth = 1.5 ) + ggplot2::scale_fill_manual( values = c( "Outside optimistic approximation" = "white", "Optimistic approximation" = "grey75", "Conservative" = "grey45" ), name = "Grid classification" ) + ggplot2::scale_colour_manual( values = c( "Conservative" = "blue", "True optimal" = "red", "Optimistic" = "orange" ), name = "Operating region", guide = ggplot2::guide_legend( order = 1, override.aes = list(linewidth = 1.5, fill = NA) ) ) + ggplot2::guides( fill = ggplot2::guide_legend(order = 2) ) + ggplot2::coord_equal( xlim = c(plot_rg_ll[1], plot_rg_ul[1]), ylim = c(plot_rg_ll[2], plot_rg_ul[2]), expand = FALSE ) + ggplot2::labs( x = expression(x[1]), y = expression(x[2]), title = "Grid approximations and optimal rectangles", subtitle = expression(x[1]^2 + x[2]^2 <= 1) ) + ggplot2::theme_bw(base_size = 13) + ggplot2::theme( panel.border = ggplot2::element_blank(), legend.position = "right", plot.title.position = "plot", legend.key.width = grid::unit(1, "cm"), legend.key.height = grid::unit(1, "cm") ) V_cons_plot <- prod(plot_hr_conservative$ul - plot_hr_conservative$ll) V_true_plot <- prod(plot_hr_true$ul - plot_hr_true$ll) V_optimistic_plot <- prod(plot_hr_optimistic$ul - plot_hr_optimistic$ll) ``` ```{r circle-grid-figure, echo=FALSE, fig.width=10, fig.height=7.5, fig.align='center', out.width='100%', fig.cap='Conservative and optimistic grid approximations of the unit disk. The conservative operating region is guaranteed feasible, the optimistic region provides a global upper volume bound, and the analytically known maximum volume lies between them.'} circle_plot ``` For this validation example, both the analytically known optimum and the optimistic outer bound are available. The two corresponding dimension-adjusted ratios are therefore ```{r circle-approximation-metrics} d_plot <- 2L R_true <- (V_cons_plot / V_true_plot)^(1 / d_plot) R_outer <- (V_cons_plot / V_optimistic_plot)^(1 / d_plot) c( conservative_to_true = R_true, conservative_to_outer_bound = R_outer ) ``` The ratio based on the true optimum is very close to one. The conservative rectangle has a somewhat longer side in the \(x_2\) direction and a shorter side in the \(x_1\) direction than the true optimal square. These opposing differences largely compensate in the product of the side lengths, so that the geometric-mean side length of the conservative solution is close to that of the true optimum. By contrast, the ratio based on the optimistic outer rectangle is substantially smaller. In this deliberately coarse example, the outer rectangle is a rather loose global upper volume bound. The smaller value therefore reflects the conservatism of the bound rather than a comparably large actual distance between the conservative solution and the true global optimum. This distinction illustrates why the optimistic ratio is guaranteed and generally available, but can underestimate the practical quality of the reported conservative solution when the outer approximation is coarse. # Adding operating-region constraints The hyperrectangle with the largest volume is not necessarily the most useful operating region in practice. Its shape may be unbalanced, with wide intervals for some factors but impractically narrow intervals for others. Subject-matter considerations may also require particular factor settings or ranges to remain available, even at the cost of a smaller total volume. `OptOR` therefore supports factor-specific minimum-width and interval-containment constraints. A minimum-width constraint requires the operating interval of a selected factor to have at least a specified length. For example, the following call requires a minimum width of \(1.8\) for the first factor: ```{r width-constraint} width_result <- optimal_cont_hr( fcts = fcts, n = 30, rg_ll = c(-1, -1), rg_ul = c( 1, 1), ctype = c("width", "none"), cwidth = c(1.8, 0), gmode = "conservative", verbose = FALSE ) width_result$ul - width_result$ll ``` Minimum-width constraints prevent the optimization from selecting factor ranges that are too narrow for routine operation. Interval-containment constraints require the final operating interval to contain a prespecified range. If the lower and upper bounds of this range coincide, the constraint ensures that a particular factor setting, such as a preferred operating point, is included. Both constraint types are implemented directly in `optimal_grid_hr()` and `optimal_cont_hr()`. The grid-aligned hyperrectangle is first optimized globally subject to the response requirements and the additional factor-specific constraints, after which `optimal_cont_hr()` performs the continuous expansion described above. In the continuous case, these requirements could alternatively be represented by additional response functions. The dedicated implementation, however, exploits constraint-specific pruning rules that can substantially reduce the search effort and thereby improve computational efficiency. These constraints also support an interactive workflow with subject-matter experts. In a graphical user interface, an expert can modify minimum widths, required intervals, or preferred factor settings and immediately assess how the optimal hyperrectangle changes. This combines formal optimization criteria, such as volume, with process knowledge and operational requirements. For the intended low-dimensional applications with \(d \leq 5\), moderate grid resolutions can often be evaluated interactively on current hardware, although runtime depends strongly on the dimension, grid resolution, geometry of the feasible region, and imposed constraints. The final operating region can thus be selected through an iterative combination of global mathematical optimization and expert assessment rather than by maximizing volume alone. # Implementation notes The computational core of `OptOR` is implemented entirely in C according to the C99 standard. Where available, OpenMP is used to parallelize computationally intensive parts of the discrete search. The R package acts primarily as a convenient wrapper for data preparation, function calls, visualization, and reproducible analysis, while the standalone C99 implementation can also be integrated into other software environments independently of the R wrapper. The discrete optimization is based on the enumeration procedure developed by Palmes, Koch, and Schaudt [@palmes2026]. For a fixed dimension \(d\), the method enumerates maximal constant hyperrectangles in \(O(n^{2d-2})\) time. Several pruning techniques reduce the computational effort in practice. However, Palmes et al. also show that this worst-case complexity bound is tight, that is, it cannot be improved in general. For quadratic response functions, conservative and optimistic grid cells are classified by computing the exact global minimum and maximum of each quadratic function over the corresponding axis-aligned subbox. Subboxes that cannot yet be classified are recursively subdivided until their classification is resolved or the target grid resolution is reached. Repeated algebraic calculations involving the unchanged quadratic coefficient matrices are cached and reused. In the intended low-dimensional applications with \(d \leq 5\), grid construction is typically not the computational bottleneck, despite the known NP-hardness of globally optimizing a quadratic function over a box, and is therefore not parallelized in the current implementation.