Complete data sets to reproduce any and all experiments are included in the Supplementary Materials. The algorithm will be publicly accessible after publication.

Benchmark Test Functions
In this section, we define the benchmark test functions, all of which are set to be maximized during our optimizations.
Contamination Problem
The contamination problem was introduced by [1]} and is used to test categorical variables with binary categories. The problem aims to maximize the reward function for applying a preventative measure to stop contamination in a food supply chain with \(D\)  stages. At each \(i^{th}\)  stage, where \(i\in [1,D]\) , decontamination efforts can be applied. However, this effort comes at a cost \(c\)  and will decrease the contamination by a random rate \(\Gamma _i\) . If no prevention effort is taken, the contamination spreads with a rate of \(\Omega _i\) . At each stage \(i\) , the fraction of contaminated food is given by the recursive relation:
\(Z_i=\Omega _i(1-w_i)(1-Z_{i-1}) + (1-\Sigma _i w_i) Z_{i-1}\) 

here \(w_i\in {0,1}\)  and is the decision variable to determine if preventative measures are taken at \(i^{th}\)  stage or not. The goal is to decide which stages \(i\)  action should be taken to make sure \(Z_i\)  does not exceed an upper limit \(U_i\) . \(\Omega _i\)  and \(\Sigma _i\)  are determined by a uniform distribution. We consider the problem setup with Langrangian relaxation [2]}:
\(f({\vec{w}})= -\sum _{i=1}^D\left(cw_i+\frac{\rho }{T}\sum _{k=1}^T1_{\lbrace Z_k>U_i\rbrace }\right) - \lambda ||\vec{w}||_1\) 

Here violation of \(Z_k<U_i\)  is penalized by \(\rho =1\)  and summing the contaminated stages if the limit is violated and our total stages or dimensions are \(D=21\) . The cost \(c\)  is set to be 0.2 and \(Z_1=0.01\) . As in the setup for [2]}, we use \(T=100\)  stages, \(U_i=0.1\) , \(\lambda =0.01\)  and \(\epsilon =0.05\) .

Encrypted Amalgamated
Analytic test functions generally cannot mimic mixed variables. To map the continuous output of a function into \(N\)  discrete ordinal or categorical variables, the continuous range of the test function's output is first discretized into \(N\)  discrete subranges by selecting \((N-1)\)  break points, often equally spaced, within the bounds of the range. Then, the continuous output variable is assigned the integer round-off value of the subrange defined by its surrounding pair of break points. If necessary, the domain of the test function's output is first mapped into a larger domain so that each subrange has a unique integer value. To mimic ordinal variables, we are done, but for categorical variables, a random vector for each categorical variable is then generated which scrambles or `encrypts' the indices of these values, thus creating random landscapes as is the case for categorical variables with a latent space. The optimization algorithm only sees the encrypted space and the random vector is only used when evaluating the black-box function.
We also define a new test function that we call the Amalgamated function, a piece-wise function formed from commonly used analytical test functions with different features (for more details on these functions we refer to [4]}). Amalgamated function is non-convex and anisotropic, unlike conventional test functions where isotropy can be exploited.
For \(i=1...n\) , \(k=\) mod\((i-1,7)\) :
\(f(\vec{w})=\sum _{i=1}^{D}g(w_i)\) 
where
\(\begin{aligned}g(w_i)=\left\lbrace \begin{array}{@{\:}l@{}l}sin(w_i) & \text{if}\: k=0, \ w_i\in (0,\pi )\\[1ex]-\frac{w_i^4-16w_i^2+5w_i}{2} & \text{if}\: k=1, \ w_i\in (-5,5)\\[1ex]-(w_i^2) & \text{if}\: k=2, \ w_i\in (-10,10)\\[1ex]-[10+w_i^2-10cos(2\pi w_i)] & \text{if}\: k=3, \ w_i\in (-5,5)\\[1ex]-[100(w_i-w_{i-1}^2)^2+(1-w_i)^2] \quad \qquad & \text{if}\: k=4, \ w_i\in (-2,2)\\[1ex]|cos(w_i)| & \text{if}\: k=5, \ w_i\in (-\pi /2,\pi /2)\\[1ex]- w_i & \text{if}\: k=6, \ w_i\in (-30,30)\end{array}\right.\end{aligned}\) 
To create the Encrypted Amalgamated function, for categorical and ordinal variables, equally spaced points are taken within the bounds defined above. For our current work, we use a \(D=13\)  with 8 categorical and 3 ordinal variables with 5 states each, and 2 continuous variables.

NK Landscapes
NK Landscapes were introduced by [5]} as a way of creating optimization problems with categorical variables. \(N\)  describes the number of genes or number of dimensions \(D\)  and \(K\)  is the number of epistatic links of each gene to other genes, which describes the `ruggedness' of the landscape. A large number of random landscapes can be created for given \(N\)  and \(K\)  values. The global optimum of a generated landscape for experimentation can only be computed through complete enumeration. The landscape cost for any vector is calculated as an average of each component cost. Each component cost is based on the random values generated for the categories, not only by its own alleles, but also by the alleles in the other genes connected through the random epistasis matrix, with \(K\)  probability or ruggedness. A \(K=1\)  ruggedness translates to a fully connected genome.
The \(NK\)  Landscapes from [5]} were formulated only for binary variables. They were extended by [7]} for multi-categorical problems, which is the formulation we use. Details of the \(NK\)  Landscape test-functions we use can be found in [7]}. For the current study, we use \(N=8\)  with 4 categories each and ruggedness \(K=0.2\) .

Rastringin
Rastringin function is a commonly used noon-convex optimization function [4]} with a large number of local optima. It is defined as:
\(f(\vec{w})= -[10+w_i^2-10cos(2\pi w_i)], \ w_i\in (-5,5)\) 
We use \(D=9\)  for testing with 6 ordinal with 5 discrete states and 3 continuous variables. The ordinal variables are equally spaced within the bounds.

Encrypted Syblinski-Tang
We use the Syblinski-Tang function [4]}, an isotropic non-convex function. The function is considered difficult to optimize because many search algorithms get `stuck' at a local optimum. For use with categorical variables, we encrypt it as described previously. The Syblinski-Tang function, in terms of input vector \({\vec{w}}\) , is defined as:
\(f(\vec{w})= -\frac{\sum _{i=1}^{D}w_i^4-16w_i^2+5w_i}{2}, \ w_i\in (-5,2.5)\) 
For the current study, this function was tested with \(D=10\)  categorical variables and 5 categories for each variable.

Encrypted ZDT6
\(ZDT\)  benchmarks are a suite of multi-objective problems, suggested by [11]}, and most commonly used for testing such problems. We use \(ZDT6\) , which is non-convex and non-uniform in its parameter space. We again modify the function by encrypting it to work with categorical problems. \(ZDT6\)  is defined as:
\(\ \begin{aligned}f_1(\vec{w})&=exp(-4w_1)sin^6(6 \pi w_1)-1 \\ f_2(\vec{w})&=-g(\vec{w})\left[1-(f_1(\vec{w})/g(\vec{w}))^{2}\right]\\ g(\vec{w})&=1+9\left[\left(\sum _{i=2}^Dw_i\right)/(n-1)\right]^{1/4}\end{aligned}\) 
 Here \( w_1 \in [0,1]\)  and \( w_i =0\)  for \(i = 2,\dots ,D\) . The function was tested for \(D=10\)  with 5 categories each. We note that to evaluate the performance of MixMOBO, we compared it against the NSGA-II variant [12]} that can deal with mixed variables (by running \(ZDT4\)  in a mixed variable setting and \(ZDT6\)  with categorical variables). No encryption is necessary for GAs. GAs required, on average, \(10^2\)  more function calls compared to MixMOBO.

Architected Materials - Manufacturing and Testing
The microlattice structures were fabricated with a hybrid organic-inorganic material Zr-DMAEMA (30 wt%). The composition of this material is 70 wt% zirconium propoxide and 10 wt% (2-dimethylaminoethyl) methacrylate (DMAEMA) (Sigma-Aldrich). Before the fabrication of the structures, the material was placed onto glass substrates and remained in vacuum for 24 h. Further information about the material preparation can be found elsewhere [13]}, [14]}.
<FIGURE>The test structures were fabricated by diffusion-assisted high-resolution direct femtosecond laser
writing, which employs MPL and the aforementioned photoresist for high-resolution fabrication. The
system consists of a FemtoFiber pro NIR laser with a wavelength of \(780 nm\) , pulse width of \(100 fs\) ,
and repetition rate of 80 MHz. Local photopolymerization of the photosensitive material was
accomplished by tight focusing the laser beam onto the material with a 100X microscope objective
lense (Plan-ApoChromat 100x/1.40 Oil M27, Zeiss).
Our experimental laser setup utilizes precise piezoelectric stages (Physik Instrumente M-110.1DG),
where their movements are controlled by a program (3dpoli, Femtika). The IGES file of each unit cell
was created in Solidworks 2019 and used to create the coordinates needed to translate the 3D model
into movements of the stages. The laser output energy for the fabrication was measured right before
the objective lens opening at 10 \(mW\)  and the scanning speed used was 80 \(\mu m/s\) . The fabrication of each scaffold was conducted by scanning the top layer of unit cells and then moving downwards to the interface on the glass substrate in a unit cell-by-unit cell manner. The stage movements were directed along the vectors of each unit cell from top to bottom avoiding slicing the scaffolds and
printing it in a layer-by-layer manner, thus ensuring the best mechanical stability and reproducibility of
the experiment. Only a single scan was made for each vector as it proved to be sufficient for the rigidity of the structures. A detailed description of the setup can be found in [15]}.
Uniaxial compression tests were performed in situ with a nanoindentation apparatus (PI 85 SEM PicoIndenter, Hysitron) mounted inside the chamber of a scanning electron microscope (FEI Quanta 3D FEG) to enable high-precision nanomechanical testing and real-time recording of the deformation of the test structures. We show the results of the compression tests in Figure REF . A flat molybdenum tip (model # 72SC-D3/035 (407A-M)) of \(90\mu m\)  diameter was used in all the compression tests. To prevent the displacement of the substrate during testing, the glass substrates on which the specimens were mounted were fixed onto an SEM pin stub mount (TED PELLA) with PELCO® Pro C100 Cyanoacrylate Glue, TED PELLA). All of the structures were aligned such that the front face to be visible during testing to facilitate the detection of the location and instant of buckling and fracture. To capture the structure collapse and track the proliferation of damage, each structure was deformed at a rate of \(500 nm/s\)  to a maximum compressive strain of \(55\%\) . To ensure repeatability, at least 5 tests were performed with structures of a given design. To capture the instigation of the first buckling instability and, consequently, compute the critical buckling load, the measured force-displacement curves were frame-by-frame juxtaposed with the recordings of the deformation to identify whether excessive deformation produced a singularity in the force-displacement response. HIM that provides high depth of focus and high resolution was used to track the evolution of damage in the interior of the tested structures [16]}.
