{'fraction_non_alphanumeric': 0.048593660116867, 'fraction_numerical': 0.018317236368040207, 'mean_word_length': 4.32055412371134, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Heuristic construction of exact experimental designs under multiple resource constraints

The aim of this paper is twofold. First, we introduce"resource constraints"as a general concept that covers many practical restrictions on experimental design. Second, for computing efficient exact designs of experiments under any combination of resource constraints, we propose a tabu search heuristic that uses some ideas of the Detmax procedure. To illustrate the scope and performance of our heuristic, we computed D-efficient designs for 1) a block model with limits on the numbers of blocks and on the availability of experimental material; 2) a quadratic regression model with simultaneous marginal and cost constraints; 3) a non-linear regression model with simultaneous direct and cost constraints. As we show, the proposed heuristic generates comparable or better results than algorithms specialized for computing optimal designs under less general constraints.


Introduction
Optimal design of experiments is an approach to constructing experimental designs using a statistically motivated utility function called optimality criterion, see, e.g., [34], [15] and [1]. Construction of an optimal experimental design is generally a very challenging problem of theoretical mathematics and numerical optimization. The main aim of this paper is to propose a heuristic method for computing efficient exact experimental designs under a large class of restrictions encountered in practice.
Suppose that we intend to perform an experiment consisting of a set of trials (runs, measurements). For each trial, we must select a design point from a finite design space X. We assume that in general it is possible to select the same design point for more than one trial.
Let φ : [0, ∞) n → [0, ∞) be an optimality criterion that measures the quality of designs for statistical inference. Often, the goal of the experimenter is to estimate unknown parameters of an underlying statistical model, and the value φ(ξ) is a measure of the information about the parameters of interest obtained from the experiment ξ, see, e.g., [37]. In view of this interpretation, it is natural to adopt the following assumption of monotonicity: (M) Augmentation (extension) of an experiment by additional trials cannot decrease its quality for statistical inference, i.e., if designs ξ and ζ satisfy ξ ≤ ζ componentwise, then φ(ξ) ≤ φ(ζ).
The most classical example of φ is the criterion of Doptimality for linear regression models with uncorrelated homoscedastic errors, as we will briefly describe. Consider an experiment with an n-point design space X. Assume that for each trial in the design point x ∈ X, the real-valued random observation Y satisfies E(Y ) = f T ι(x) β and Var(Y ) = σ 2 < ∞, where f ι(x) ∈ R m is a regressor vector, ι(x) ∈ {1:n} is the index of x, β ∈ R m is a vector of unknown parameters of interest, and σ 2 is a constant variance. For different trials, the observations are assumed to be uncorrelated. Then, the criterion of D-optimality φ D : [0, ∞) n → [0, ∞) is defined by It is possible to show that φ D is continuous, concave and monotonic on [0, ∞) n in the sense of Assumption (M). Additionally, φ D is homogeneous, that is, φ D (cξ) = cφ D (ξ) for any design ξ and any c ≥ 0. Hence, the quality of two designs can be compared by their relative D-efficiency defined by eff D (ξ|ζ) = φ D (ξ)/φ D (ζ) for all ξ and ζ such that φ D (ζ) > 0. A design ξ * is called D-optimal, if it maximizes the value of φ D in a given set of competing designs. From the statistical point of view, the D-optimal design minimizes the generalized variance of the best linear unbiased estimator of β. For details, see [34], [15], [37], and [1]. Usually, the set of designs is only restricted by the number of trials, i.e., for each feasible exact design ξ, where N is a maximum "size" of the experiment. 2 This corresponds to the situation where each trial amounts to the same cost and the experimental budget allows performing at most N trials. Alternatively, condition (1) can represent the situation where the trials must be performed sequentially, each trial lasts the same amount of time, and the deadline permits performing at most N trials. However, there are practical situations where any feasible design must satisfy one or more constraints different than (1). In this paper, we propose to consider general "resource" constraints of the form where a ri represents the consumption of the r-th resource by a single trial in the i-th design point, and b r represents a limit on the r-th resource. 2 In fact, the usual requirement is that the number of trials is exactly equal to N , but assumption (M) implies that this requirement is equivalent to (1).
The interpretation of (2) leads to the following assumptions: (C1) Resource limits are positive and finite, i.e., b 1 , ..., b k ∈ (0, ∞). (C2) Augmenting designs can not decrease the consumption of any resource, i.e., a ri ≥ 0 for all r ∈ {1:k} and i ∈ {1:n}. (C3) No trial is completely free and its replication must eventually result in exceeding some resource limit, i.e., for all i ∈ {1:n} there is some r ∈ {1:k} such that a ri > 0.
Clearly, constraint (1) is a special case of (2) with k = 1, a 1i = 1 for all i ∈ {1:n}, and b 1 = N . Sometimes, however, the costs depend on design points and the budget of the experiment is limited by B financial units (e.g., Section 6 in [13], or [33]). This restriction can be formalized by one resource constraint (k = 1) such that b 1 = B, and a 11 , ..., a 1n represent possibly unequal costs of trials in individual design points. Another class of constraints corresponds to the socalled marginal restrictions (e.g., [8], [29]) or, more generally, strata restrictions ( [22]). In this case, design space X is partitioned into non-overlapping sets X 1 , ..., X k and any experimental design ξ must satisfy x∈Xr ξ ι(x) ≤ s r for all r ∈ {1:k}, where s 1 , ..., s k are given positive numbers. Here, we obtain (2) by setting a rι(x) = I[x ∈ X r ] (the indicator function) and b r = s r for all r ∈ {1:k}.
A natural type of restrictions are the so-called direct constraints (e.g., [14]), which correspond to performing at most l 1 , ..., l n trials in design points i −1 (1), ..., i −1 (n), respectively. Often, it is possible to perform at most one observation in any design point, i.e., l 1 = ... = l n = 1. This can be converted to the resource constraints (2) by setting k = n, a ri = δ ri (the Kronecker delta) and b r = l r for all r ∈ {1:k}.
The general resource constraints can also accommodate limits on the availability of treatment samples in block designs ( [3]), required time delay between consecutive trials (Section 5 in [39]), and other experimental constraints. See also the review paper [9] or Chapter 4 in [15].
However, we remark that there are some other possible restrictions that can not be represented by (2), for instance limits on the transitional costs (cf. [41]).
Let ξ (0) be either an exact design representing trials that have already been performed, or a required initial part of the experiment. We will assume that ξ (0) satisfies (2). In the situation without existing/required trials, the design ξ (0) is simply the zero vector 0 n ∈ R n .
In a matrix form, the system (2) can be written as Aξ ≤ b componentwise, where A is the k × n matrix of coefficients a ri , r ∈ {1:k}, i ∈ {1:n}, and b is the kdimensional vector of b 1 , ..., b k . The assumptions (C1)-(C3) guarantee that the set of all feasible exact experimental designs is non-empty and finite 3 . The assumptions also imply that the set of all feasible approximate designs is a non-empty, compact and convex polyhedron.
For small-size problems of type (3), it is possible to use an enumeration method, such as branch-andbound, that guarantees finding a globally optimal solution (e.g., [43]). Nevertheless, there is no practical hope of creating an algorithm that rapidly produces provably optimal solutions of large instances of (3). Often, the only possibility is to use a heuristic that usually leads to an efficient feasible experimental design.
A natural approach to solving (3) is to use a heuristic based on "excursions" within the set of designs, as in some early algorithms for computing D-optimal experimental designs under the standard constraint (1), see, e.g. [11], [47] and [31]. From these algorithms, the most relevant for our problem is the Detmax procedure proposed by Mitchell ([31]), which is related to the tabu search methods (e.g., [17], see also [20]).
Today, the most popular for solving the standard optimal design problems are the exchange heuristics (e.g., Chapter 12 in [1], and [19]). However, the basic principle of the exchange heuristics is difficult to modify to solve instances of the general problem (3), since the number of trials of the optimal resource constrained design is not always known in advance. Moreover, an exchange of two design points may render a feasible design non-feasible.
As far as general resource constraints are concerned, an idea similar to Detmax has already been used in reliability theory to efficiently solve redundancy optimization problems (see [25]). In the area of optimal design of experiments, another related procedure has been proposed, based on a sequential removal of a single design point and a subsequent augmentation of the resulting design by a greedy method ( [46]). Finally, the paper [3] describes a randomized heuristic similar to simulated annealing, that can be used to solve the general problem (3).
For computing efficient solutions of resource constrained problems, we will propose a simple and robust heuristic related to the Detmax procedure, employing a tabu search principle; see Section 2. We will explore the performance of the heuristic as a tool for computing D-efficient exact designs of regression experiments with uncorrelated homoscedastic errors and linear resource constraints of type (2) used in diverse application areas; see Section 3.

General description of the heuristic
We will say that a design ζ ∈ Ξ ex is created from a design ξ ∈ Ξ ex by a forward step (or a backward step) if ζ = ξ+e i (or ζ = ξ−e i ) for some standard unit vector e i ∈ R n . A design ξ ∈ Ξ ex will be called "maximal" if it can not be augmented without violation of some of the resource constraints, that is, if all designs created from ξ by a forward step are non-feasible.
For a design ξ ∈ Ξ ex , an upper neighbour is any feasible design that can be obtained from ξ by a forward step, i.e., the set of all upper neighbours of ξ is Similarly, a lower neighbour of a feasible design ξ is any design that can be obtained from ξ by a backward step, that is, the set of all lower neighbours of ξ is Note that L(ξ (0) ) = ∅, and U(ξ) = ∅ if and only if ξ is maximal. We will also assume that ξ (0) is not maximal, which means that L(ξ)∪U(ξ) = ∅ for any feasible design ξ.
Clearly, properties (C1)-(C3) imply that any feasible exact design is reachable from any other feasible exact design by a sequence of forward and backward steps within Ξ ex . Moreover, an optimal solution of the problem (3) can be found among maximal designs, in view of assumption (M).
The proposed algorithm starts in a design ξ (1) ∈ Ξ ex (cf. Subsection 2.2) and builds an excursion in the set of feasible exact designs, guided by a "tabu" list V of characteristic attributes of the designs that have already been visited.
Let attr (ξ) be a characteristic attribute of ξ ∈ Ξ ex and let val (ξ) be a local heuristic evaluation of ξ, i.e., a real number that roughly estimates how promising is ξ as a part of an excursion leading to an efficient design (see Subsections 2.3 and 2.4 for a more detailed specification of attr and val ).
Let ξ represent the current design in the excursion. The algorithm first attempts a forward step (if attr (ξ) / ∈ V ) or a backward step (if attr(ξ) ∈ V ), moving to a neighbouring feasible exact design ζ. Design ζ is chosen such that it maximizes val among all designs satisfying attr(ζ) / ∈ V . If the algorithm attempts a forward step but there is no ζ ∈ U(ξ) such that attr(ζ) / ∈ V , or if it attempts a backward step but there is no ζ ∈ L(ξ) such that attr(ζ) / ∈ V , the algorithm tries to reverse the direction of the search. If all these attempts fail, i.e., if the attributes of all neighbouring designs of ξ are contained in the list V , the algorithm resolves this "blockage" by randomly selecting a design from L(ξ) ∪ U(ξ) for the next step.
Each time a maximal design is encountered, the algorithm checks whether it is better than the best feasible design ξ + found so far. If the number of backward steps of an excursion exceeds a constant back max , the excursion is declared to be a "failure" and the algorithm is restarted from the currently best design. Note that the list V is not cleared after the restart, that is, the new excursion will follow a different path. The algorithm is terminated once the computation time exceeds a user-supplied time limit t max . The idea of the algorithm is made more precise by its meta-heuristic scheme; see Algorithm 1.
Thus, the algorithm is similar to the Detmax procedure, because it attempts a forward or a backward step depending on an attribute of the current design. Since most of the designs encountered have not been previously visited, the excursions tend to move towards maximal designs. Note, however, that even under the standard constraint, Algorithm 1 differs from the Detmax procedure in several important aspects. For instance, unlike Algorithm 1, the Detmax procedure does not avoid the backward steps to the designs with attributes in the tabu list, which often leads to retracing the same excursions. Input : The model, matrix A, vector b, design ξ (0) to be augmented, initial feasible design ξ (1) , required time of computation t max , criterion φ, maximum number of backward steps back max . Output: A design ξ + as the best found feasible solution of (3).
Algorithm 1: A general scheme of the proposed heuristic for computing efficient experimental designs under constraints (3). For linear regression, the model at the input can be represented by an m × n matrix F = (f 1 , ..., f n ) of regressors corresponding to individual design points. The model at the input is used inside the functions φ, val , and can also be used inside attr. Matrix A and vector b define the resource constraints. They are implicitly used for computing the set U(ξ), and can also be used inside val . In our implementation, the function attr used an additional constant n round ; see Subsection 2.3. Design ξ (0) is implicitly used for computing the set L(ξ). The function time returns the time elapsed from the start of the computation.

Choice of the initial design
An important part of Algorithm 1 is the choice of the initial design ξ (1) . Our experience shows that a reasonably efficient design can usually be obtained by choosing ξ (1) = ξ (0) , but for more complex problems it is better to use multiple restarts of the heuristic, with initial designs created by a sequence of random forward steps starting from ξ (0) .

Choice of the characteristic attributes of designs
The characteristic attribute should be chosen such that, loosely speaking, it assigns different values to substantially different designs and the same values to essentially same designs (for instance to algebraically isomorphic designs).
After some experimentation, we have decided to use the attribute attr(ζ) equal to the value φ(ζ) rounded to n round significant digits. Note that instead of storing complete designs, storing real-valued attributes in the list V not only makes the time and memory requirements much smaller, but makes the tabu principle itself more efficient in some optimum design problems.

Choice of the local heuristic evaluation of designs
In our implementation of Algorithm 1, the local heuristic evaluation val (ζ) of a design ζ is an estimate of the maximal value of φ on the set of all designs augmenting ζ, i.e., val (ζ) is an estimate of val * (ζ) = max{φ(η) : η ∈ Ξ ex , ζ ≤ η}.
The rationale behind this particular evaluation is that if we were able to use the exact values of val * , the initial greedy phase of the algorithm started from ξ (0) would directly lead to a globally optimal solution.
For ζ ∈ Ξ ex , let be the vector of residual amounts of resources. Note that after a forward or a backward step, it is possible to use the update formula r(ζ ± e i ) = r(ζ) ∓ (a 1i , ..., a ki ) T , i ∈ {1:n}.
Note that the heuristic evaluation of designs is chosen such that it depends only on the set Ξ ap of feasible designs, not on the choice of the algebraic definition of Ξ ap . Furthermore, the excursions only depend on the ordering of approximate designs determined by the criterion φ, not on the chosen "version" of the same criterion.
Clearly, there are many other methods of computing a local design evaluation in Algorithm 1. For instance, it is possible to use some version of the direct greedy method based on the relative change of φ with respect to a change in residual resources, similarly to [25]. These methods may allow for a more rapid construction of the excursion, nevertheless, they may also lose the above-mentioned invariance properties. Moreover, in our experiments with various modifications of the move selection rules, we did not observe a significant increase in the quality of results.

Examples
In this section, we will apply Algorithm 1 to the most common situation in optimal design of experimentscomputation of D-efficient experimental designs for linear regression models with uncorrelated homoscedastic errors,as described in the introduction. Although the chosen criterion is always the same, the selected optimization problems have very different sets of feasible designs. Our experience suggests that the feasible set has a more pronounced effect on the complexity of the optimization problem (3) than the choice of the criterion within the class of standard criteria used for optimal design.
To demonstrate the universality of Algorithm 1, we selected the same heuristic parameters in all examples (namely, back max = 16 steps and n round = 9 significant digits). We ran all computations for t max = 120 seconds 4 . To illustrate the statistical distribution of the quality of results and detect potentially difficult instances of the optimization problems, we always used a set of 10 independent initial designs generated by a random sequence of forward steps starting from ξ (0) . In Example 3.1, we used the R computing environment ( [38]), and in Examples 3.2 and 3.3 we used Matlab. The codes can be found at www.iam.fmph.uniba.sk/design/ All examples were computed on a 64 bit Windows 7 system running an Intel Core i5-2400 processor at 3.10 GHz with 4GB of RAM.

Designs for a block model with a constraint on the number of blocks
Consider a block model with N blocks of size two and v treatments. More precisely, assume that the uncorrelated observations Y 1 , ..., Y N satisfy where t 1 (j), t 2 (j) ∈ {1:v} are the treatments selected for the j-th block, with effects τ (t 1 (j)), τ (t 2 (j)), and Var(Y j ) = σ 2 < ∞, j ∈ {1:N }. An experimental design is given by a selection of treatments t 1 (j) and t 2 (j) to be compared in the j-th block, for all j ∈ {1:N }. Optimal designs for this model have been applied in two-channel microarray experiments (e.g., [44], [4]) and elsewhere.
In this setting, the design space can be viewed as the set of all possible pairs of treatments, i.e., X = {(1, 2), (1, 3), . . . , (v − 1, v)}, which can be indexed by The problem of the so-called Doptimal block designs is then equivalent to the standard D-optimal design problem as described in the introduction, with m = (v − 1)-dimensional regressors for all (t 1 , t 2 ) ∈ X, cf. [39].
In this example, the aim is to demonstrate that Algorithm 1 performs well under the standard constraint (1), i.e., if the only restriction is not to exceed the given number N of blocks.
We implemented Algorithm 1 in the environment R and used it to compute D-efficient designs for v = 16 treatments and N = 15, . . . , 120 blocks. We then compared the designs with the results of a simulated annealing procedure od implemented in R package "smida" (see [45]), with parameters criterion = "D", dye = FALSE, and the number n.iter of iterations chosen such that the computation time is approximately t max = 120 seconds 5 . Figure 1 shows that Algorithm 1 systematically produced either the same or better results than the simulated annealing method (with small exceptions for N = 28, 36, 53, 54, 55, 78). The numerical results suggest that the simulated annealing procedure has difficulties if N is a multiple of 8.
Any block design with blocks of size two can be represented by a "concurrence" graph with v vertices and N possibly multiple edges, such that the endpoints of edges correspond to the treatments used in the same blocks (e.g., [5]). Kirchhoff's matrix tree theorem implies that φ m D (ξ) is equal to the number of spanning trees of the concurrence graph of design ξ. Thus, the problem of D-optimal designs for this specific statistical model is equivalent to the problem of t-optimal graphs, that is, the concurrence graph of the D-optimal design maximizes the number of spanning trees in the class of graphs with fixed number of vertices and edges.
For some numbers N and v, the D-optimal designs (or t-optimal graphs) are known theoretically (see [6], [16], [35]). For instance, it is known that a complete almost-regular multipartite graph is t-optimal among all simple graphs with the same numbers of vertices 5 The main application area of Algorithm 1 is computing efficient designs under non-standard constraints. Therefore, we did not perform a detailed comparison of Algorithm 1 with the vast number of other known methods applicable to computing optimal design under the standard constraint. and edges. The number of spanning trees for a complete multipartite graphs with v vertices and p ≥ 2 partitions of sizes k 1 , ..., k p is given by ([2], [27]) With this formula, we can calculate the optimal value of the D-optimality criterion for v = 16 and N =64,  85, 96, 102, 106, 109, 112, 113, . . . , 120. It turns out that Algorithm 1 consistently finds the theoretically Doptimal designs for all of these values of N .
Compelling candidates for D−optimal designs are those that are represented by strongly regular graphs (srg), because of their high degree of symmetry. The imprimitive strongly regular graphs are either disconnected graphs or complete multipartite graphs with the partitions of the same size (D−optimal, as mentioned above). For v = 16 vertices there exist four primitive strongly regular graphs (see [10], Chapter VII.11). One of them, srg (16,10,6,6), that is, the Clebsh graph with N = 40 edges, was obtained by our heuristic, and we conjecture that it is D-optimal. However, the remaining three of the strongly regular graphs 6 are not Doptimal; their efficiencies compared to the designs found by Algorithm 1 are 98.65%, 99.68%, and 99.61%, respectively. In Figure 2, we depict the concurrence graph representation of the designs obtained by Algorithm 1 for v = 16 and N = 48, 72, 80. Interestingly, all these graphs contain a large number of complete bipartite subgraphs.

Designs for a quadratic model with simultaneous marginal and cost constraints
Consider the D-optimal design problem for sintering uranium pellets that are to be used as a fuel in nuclear plants, as discussed in [29]. The explanatory variables represent the "initial density" (x 1 ) and the "percentage of additive U 3 O 8 " (x 2 ). The statistical model relating response and the explanatory variables is assumed to be the full quadratic linear regression model with uncorrelated homoscedastic errors determined by the regressors In (7), it is assumed that (x 1 , x 2 ) lies in the design space for all (x 1 , x 2 ) ∈ X such that x 1 ≥ 95.1. The nature of the experiment requires marginal constraints on variable x 1 representing available experimental material (uranium rods). If we denote the required marginal sums by b 1 , . . . , b 18 , the constraints on a feasible design ξ are (1,3,14,59,52,29,25,32,36,29,36,38,12,10,8,2,3,3).
Furthermore, we supposed that one percent of the additive costs one price unit (cf. [23], [39]) and the financial resources of the experimenter are limited. Therefore, we solved the problem with additional constraints of the form where B is a maximum possible cost of the experiment. We varied the maximum cost from 1100 to 3900 price units with a step 50 and, for each B, we used Algorithm 1 to compute 10 exact designs maximizing the criterion of D-optimality.
To express the quality of the resulting designs, we computed their D-efficiencies relative to the approximate D-optimal designs obtained by maxdet programming (see [42]). Note that the efficiencies relative to the approximate optimal designs represent lower bounds on the efficiencies relative to the (unknown) exact optimal designs. Figure 3 shows that in all cases the efficiencies were higher than 99.99%. Additionally, the results are very stable in spite of the completely random selection of initial designs.
There are two mathematical programming methods that can be applied to the constrained problem from this subsection. The heuristic method based on integer quadratic programming (IQP; [23]) is often fast and simple to use, but, for this particular problem, it tends to produce worse results than Algorithm 1. The approach based on mixed integer second order cone programming (MISOCP; [39]) gives more efficient designs than IQP, but its results are still slightly worse than the results of Algorithm 1.
To show a concrete example of a marginally and cost constrained design, we chose B = 1965; see Figure  4 for the result. The exact design obtained in [39] by MISOCP has D-efficiency of about 1 − 10 −5 relative to the exact design obtained by Algorithm 1. The IQP and the MISOCP methods are more complex than Algorithm 1, often provide worse designs, and require an advanced integer programming solver. Nevertheless, note that they can be applied under more general linear constraints than Algorithm 1. Moreover, the MISOCP method provides a non-trivial lower bound on the efficiency of the resulting design.

Designs for a non-linear regression model with simultaneous direct a cost constraints
The third example is taken from [46]. Suppose that we wish to find the best sequence of sampling times for a model relating time and internal concentrations of fluoranthene in an organism. The mean internal con-   Fig. 4 The marginally constrained D-optimal approximate design and the best marginally constrained exact design found by Algorithm 1, with the additional cost constraint B = 1965. The relative weights of the D-optimal approximate design are denoted by the sizes of the grey discs. The integer numbers represent the exact design found by Algorithm 1. Note that the exact design is very close to the optimum since its D-efficiency is higher than 99.92%. See Subsection 3.2 for details.
centration at time t of the experiment is given by where θ 1 and θ 2 are parameters corresponding to the constant uptake and elimination rates. The experiment will be initiated at a starting time s and all observations need to be performed within the following 144 hours. Hence, an appropriate designs space is X = {0, 1, ..., 144} and the elements t ∈ X represent the time (in hours) elapsed from s. The model in consideration is non-linear, therefore we will compute the locally D-optimal designs (see, e.g., [36]). To this end, we need to linearise the model in some fixed parameters ϑ 1 and ϑ 2 . Since the model is linear in θ 1 , the choice of ϑ 1 is irrelevant. For θ 2 , we will select the nominal value ϑ 2 = 0.2381 suggested in [46] by earlier experiments.
Thus, we will consider the D-optimal design problem for the linear regression model with uncorrelated homoscedastic errors and two-dimensional regressors where t ∈ X is the design point corresponding to the time of the observation, ι(t) = t + 1 is the index of the design point, and ∇ denotes the gradient.
The experiment requires observations at t = 0, t = 72, and t = 144 hours of the experiment. Thus, the de-sign ξ (0) to be augmented satisfies ξ In accord with [46], we also assume that the experimental budget of the practitioner is limited by B = 13 price units. Moreover, the sampling costs vary throughout the week. For a starting time s, we can divide the design space X as where X s 1 denotes the sampling times with regular hourly wage on weekdays (8am -5pm, Monday to Friday), X s 2 denotes the sampling times with double wage on the weekend (7pm Friday -6am Monday), and X s 1.5 denotes the sampling times with 1.5 of regular wage (all other times). Hence, if c is the cost of taking a sample at a time with a regular hourly wage, the cost of taking one sample in X s 2 will be 2c and the cost of taking a sample in X s 1.5 will be 1.5c. We are interested in finding optimal designs that do not exceed the budget 13c. Additionally, we can perform at most one observation in each design point.
Formally, the constraints required by the experimental set-up can be expressed as: We have used Algorithm 1 as well as the heuristic from [46] 7 to compute D-efficient designs for starting times s = 0, . . . , 167. To express the quality of the obtained designs, we have evaluated their efficiencies relative to the locally D-optimal approximate designs computed by maxdet programming; see Figure 5. Similarly to the previous examples, the results of Algorithm 1 are very stable; all random restarts resulted in the same design except for s = 27, 34, 36, 108. Moreover, all 1680 results of Algorithm 1 were the same or better than the corresponding results from [46] with the following exceptions: 10 results for s = 12, 3 results for s = 27, 1 result for s = 34, and 10 results for s = 35. In the most problematic case of s = 35, the D-efficiency of all 10 results of Algorithm 1 is only 94.46% relative to the design found by the heuristic from [46]. Our computational experiments show that for s = 35 Algorithm 1 detects the optimal design after 200 to 300 seconds, depending on the initial design.
Concrete examples of the experimental designs (for the starting time s = 72) are depicted in Figure 6. For this case the relative efficiency of the design found by Algorithm 1 with respect to the approximate D-optimal 7 We used the Matlab code provided on the web page of the authors of [46]. design is 99.56%, whereas for the design obtained by [46], the D-efficiency is 94.09%.

Conclusions
We showed that the resource constraints (2) cover many types of experimental design restrictions, and that the optimal design problems associated with these restrictions can be efficiently solved by a common heuristic. For simplicity, we computed the numerical examples for the criterion of D-optimality, but it is straightforward to apply the heuristic to any monotonic criterion. Moreover, the algorithm can be as easily applied to statistical models different from the standard regression.
There are many variants of the proposed heuristic that could enhance its performance in specific situations. Except for alternative choices of initial designs, formulas for designs' characteristic attributes and local heuristic evaluations, it might also be possible to improve the efficiency of the heuristic by a different notion of a "failed" excursion, or variations in the manipulation with the list V . Since the heuristic is based on forward and backward steps, the speed of the execution could also be enhanced using the update formulas from [1], Chapter 12.
Clearly, multitude of nature-inspired optimization heuristics, such as physical, evolutionary and swarm algorithms are also applicable to solving problems of type (3), either directly or using a penalty approach to take the constraints into account (cf., e.g., [40], [30], [12], cf. also [18], [32], [28] for the applications of these methods to the standard optimum design problem). However, these methods usually require a large amount of programmer's experience, numerical experimentation and fine tuning of parameters to fit the specific properties of the optimization problem at hand.
Thus, except for introducing the general resource constrained problem in the area of experimental design, a secondary aim of this paper was to provide a simple, universal, yet reasonably efficient benchmark method for testing more advanced techniques that might be developed in the future.  Fig. 6 The D-optimal approximate design (upper panel) and two exact designs (middle and lower panels) for the model from Subsection 3.3 with starting time s = 72. The graph on the middle panel represents the exact design obtained in [46] and the graph on the lower panel represents the exact design obtained by Algorithm 1.