#  Class-incremental Learning using a Sequence of Partial Implicitly
Regularized Classifiers

1st Sobirdzhon Bobiev  Institute of Data Science and Artificial Intelligence
Innopolis University
Innopolis, Russia, 420500
<EMAIL_ADDRESS>2nd Adil Khan  Institute of Data Science and
Artificial Intelligence
Innopolis University
Innopolis, Russia, 420500
<EMAIL_ADDRESS>3rd Syed Muhammad Ahsan Raza Kazmi  Institute of Secure
and Cyber Physical System
Innopolis University
Innopolis, Russia, 420500
<EMAIL_ADDRESS>
###### Abstract

In class-incremental learning, the objective is to learn a number of classes
sequentially without having access to the whole training data. However, due to
a problem known as catastrophic forgetting, neural networks suffer substantial
performance drop in such settings. The problem is often approached by
experience replay, a method which stores a limited number of samples to be
replayed in future steps to reduce forgetting of the learned classes. When
using a pretrained network as a feature extractor, we show that instead of
training a single classifier incrementally, it is better to train a number of
specialized classifiers which do not interfere with each other yet can
cooperatively predict a single class. Our experiments on CIFAR100 dataset show
that the proposed method improves the performance over SOTA by a large margin.

###### Index Terms:

continual-learning, class-incremental learning, catastrophic forgetting

Figure 1: Model architecture. It consists of a feature extractor F (pale-blue
trapezoid) and a growing number of classifiers (orange trapezoids), followed
by the final layer $S$ which concatenates all classifier outputs together and
applies the softmax function. At each training session $i$, only the newly
added classifier $\textbf{C}_{i}$ is trained while the previous classifiers
and the feature extractor are kept frozen. The frozen modules are represented
by dashed borders.

##  I Introduction

Artificial neural networks(ANNs) have been at the top of the machine learning
landscape for a while. They have achieved impressive performances across
various applications, including object recognition[1], anomaly detection[2,
3], accident detection[4, 5], action recognition[6, 7, 8], scene
classification[9], hyperspectral image classification[10, 11], medical image
analysis[12, 13], machine translation[14, 15] etc. Their success is mainly
attributed to the availability of large amounts of data and sufficient
computational power. While ANNs are inspired by the biological brain, still
there are some shortcomings that makes them different. Specifically, they are
not designed to learn in incremental way, like humans do. Humans keep learning
new knowledge throughout their lives. However, experiments show that ANNs
almost completely forget their previous knowledge when they are trained on a
new task[16]. This is termed as “Catastrophic Forgetting”[17].

The common way to train ANNs is to provide them the whole dataset at once, and
let them iterate through it multiple times. Unfortunately, there are
situations where this is not a feasible option. There can be memory limits in
the learning device, making it impossible to store all data. Also, there can
be security concerns for storing the data if it contains sensitive
information. Catastrophic forgetting in ANNs was first addressed by
McCloskey[17] in 1989, but is still an unsolved problem hindering the progress
towards building AI agents that can learn continuously.

Three main continual learning scenarios are identified: task-incremental
learning, class-incremental learning, and domain-incremental learning [18]. In
task-incremental learning, the model is required to learn a sequence of tasks
sequentially, and during inference it will be provided with the task identity.
On the other hand, in class-incremental learning, no such information is
provided at inference time, and the model is required to predict both the
correct class and the task. A slightly different scenario is domain-
incremental learning, which, unlike class-incremental learning, does not
require predicting the task identity.

This paper is concerned about class-incremental learning, where at each
training session multiple new classes are to be learned while also maintaining
the knowledge previous classes. Despite many proposed solutions in the
literature, the baseline methods like “Experience Replay” (ER)[19] and
GDumb[20] are still shown to be as effective as the state of the art. Given
that each training session contains examples of only a few classes, the deep
networks are prone to overfitting. Therefore, it is plausible to use a frozen
feature extractor and smaller classifiers on top. ER trains only a single
network, extending its outputs units for the new classes. In each training
session, the network is trained on the examples of new classes as well as a
small number of examples from previous classes which have been stored in the
limited memory buffer. While retraining on these few stored examples helps
retain the knowledge of past classes, there is still no guarantee that the new
knowledge will not interfere with past. To solve this problem, we propose to
train a separate classifier for each group of new classes. These classifiers
do not share any weights with each other, implying that the newly acquired
knowledge will be stored separately without overriding others.

##  II Related Work

In the following, we present some of the continual learning method existing in
the literature. We divide these methods into three groups which mainly rely on
one of the three ideas: replay, regularization, and architectural techniques,
which are all described in the following subsections.

###  II-A Replay Methods

A straightforward solution to prevent catastrophic forgetting is to revisit
the previous tasks. Rehearsal methods accomplish this by storing a small
portion of the seen examples for later retraining. When the model is faced
with new training data, it augments it with memory samples to reduce the
forgetting of previous knowledge.

iCaRL[21] uses the nearest-mean-of-exemplars classification approach,
classifying items to the class with the nearest center. It uses a heuristic
approach for updating the memory buffer, prioritizing items based on their
proximity to their class mean. Together with network distillation, it has been
able to learn in class-incremental learning scenario. Gradient Episodic Memory
(GEM)[22] is designed for task-incremental learning from streaming data. In
this setting, the model receives a series of tuples $(x,t,y)$ consisting of
input, task label, and target, respectively. GEM uses the memory samples not
for replay but to serve the inequality constraints that prevents the loss for
the past tasks to decrease. Replay using Memory Indexing (REMIND) [23] applies
quantization to the extracted feature maps from deeper layers of a CNN and
stores their indices in memory. This results in a compressed representation
and allows for much more number of stored examples. REMIND’s architecture
consists of a frozen feature extractor followed by a trainable classifier.
REMIND was shown to outperform existing approaches in a streaming setup on the
datasets ImageNet, and CORe50[24].

##### Rethinking Experience Replay[19]

The authors emphasize that simple experience replay(or simple rehearsal)[25]
is still as effective as state of the art if certain issues are resolved. As
stated by the authors, the three important issues with this approach are:
overfitting to stored examples, biased prediction and accuracy towards the
later classes, and non-i.i.d stored data in cases where the buffer is small.
Their proposals include an additional bias correction layer (BiC)[26],
exponential decay of learning rate, and balanced sampling for the memory
buffer. In their reported results ER has outperformed SOTA sophisticated
replay methods such as iCaRL[21], GEM[22], A-GEM[27], and HAL[28], sometimes
by a very large margin.

##### GDumb [20]

is a recent paper that questioned our progress in continual learning. It
proposes a simple and most generic baseline for continual learning. It
basically maintains a balanced memory buffer and only trains on the samples
contained in it. The reported results have shown that it outperforms many SOTA
methods in their respective settings for which they were designed to. GDumb
serves a strong baseline for all continual settings including class-
incremental learning.

###  II-B Regularization Methods

These methods impose a type of regularization that helps retain the learned
knowledge of past tasks. They do it by adding additional loss terms to the
loss function.

Learning without forgetting (LwF)[29] presented a modification to the standard
fine-tuning. A new network is initialized as a copy of the old one with an
extension of the output layer(called multihead approach) for the new task. A
distillation loss is added between old and new task heads so the new network
is reminded of the old tasks indirectly. In backward pass, only the new
network is updated.

Elastic Weight Consolidation (EWC)[30] and Synaptic Intelligence (SI) [31]
both try to approximate the importance of each of the parameters for previous
tasks and selectively apply regularization to limit their change.
Regularization methods alone, as several papers demonstrate, are not
sufficient for proper class-incremental learning[16, 32].

###  II-C Architectural Methods

Architectural methods try to manipulate weights, neurons, layers, or
architecture of the network to protect the learned knowledge while acquiring a
new one. They either use fixed or dynamic architectures. These methods have
the advantage of completely eliminating interference between tasks while
allowing knowledge transfer between them. However, practically, they are
coupled with scalability issues as the number of tasks grow.

“Progressive Neural Networks”(PNN)[33] inherently target the task-incremental
learning. Their architecture grows laterally, each time adding a new neural
network, called “columns”, for a new task. The new columns have connections to
other layers of previous networks and thus, highly benefit from knowledge
transfer. Once a column is trained, it will be kept frozen making the PNNs
completely immune to catastrophic forgetting. A big issue with PNNs is that
they keep growing too large, limiting their use in practice to only a small
number of tasks. “Compacting, Picking and Growing” (CPG) [34] tries to
overcome this issue by dynamically controlling the architecture of the network
in an efficient way. The weights of the network are grouped into “compact” and
“free”, where the newly added weights are considered free until they are
trained and then compacted. For each new task, a learnable binary mask is
created and applied to the compact weights to select a subset of them. Then
this task is learnt using this subset of compact weights and the free weights.
The network is grown during the training if the performance does not reaching
a satisfactory level, providing more free weights. The training is followed by
pruning to compact the newly learned weights. This method is useful in
scenarios with less number but larger tasks.

####  II-C1 Similarity with other works

Here we specifically mention some of the similar works in the literature and
point out the differences with our ours. Our main difference is that we train
a set of specialized classifiers dedicated to each class group while also
making them able to detect if a sample is from previous classes.

Aljundi et al. [35] trained a specialized model for each new task and proposed
a method for choosing the relevant one at inference time. Specifically, by
training separate autoencoders describing each task’s data distribution,
during inference they are able to choose the most relevant model for which the
corresponding autoencoder has the least reconstruction error. After an
autoencoder for a task is trained, it is used to select the most relevant task
to the current task. Then the classifier model is trained based on the most
related tasks.

AR1[36] trains a linear classifier on top of a feature extractor. In each
phase, a new linear classifier is trained for the current group of classes.
Then a mean-normalization is applied to the weights of this classifier to
eliminate the prediction bias. In testing, the prediction of the model is the
softmax over outputs of all classifiers. The difference with our method is
that we use deeper classifiers, while they do not. On the other hand, we use
memory buffer, but AR1 doesn’t use one.

##  III Methods

Our purpose is to sequentially train on a number of disjoint datasets
containing different classes. Relying on a pretrained deep feature extractor,
we only consider training small networks on top of it. The traditional way is
to construct and train a single head with expanding output units to facilitate
the prediction of new classes. However, despite the replay mechanism, the
single head is still not sufficiently equipped against catastrophic
forgetting. We improve the situation by instead training much smaller
classifiers, one at a time that are specific only to new classes in the
training session. For a descriptive diagram of our method, please see Fig. 1.
In the following paragraph, we further formalize the training setting and
details of our method.

In class-incremental learning we aim to train a network on a sequence of
datasets $\mathcal{D}_{i}={(x_{i}^{j},y_{i}^{j})}_{j=1}^{n_{i}}$ with inputs
$x_{i}^{j}$ and labels $y_{j}^{i}\in\mathcal{Y}_{i}$, where
$\mathcal{Y}_{i}\cap\mathcal{Y}_{k}=\emptyset$ for any $i\neq k$. Note that in
the training session $i$ only the dataset $\mathcal{D}_{i}$ is available.
Suppose we are allowed to store a small number of examples from the current
training session in a limited memory buffer. We denote the samples in the
memory buffer by $\mathcal{M}_{i}$ and its capacity by $B$, thus
$\mathcal{M}_{i}\leq B$. The memory can be updated with new examples at each
training session. Having access to past examples through the memory buffer,
the training data for the $i$-th session will then consist of
$\mathcal{D}_{i}\cup\mathcal{M}_{i}$.

Let $\mathbf{F}(\cdot)$ be the feature extraction network. Our approach is
that at $i$-th training session we create a new classifier $\mathbf{C}_{i}$
network with $|\mathcal{Y}_{i}|$ output units that classifies between the new
classes. The final classification is done through the final layer
$S(\cdot)=\text{Softmax}(\oplus(\cdot))$ which simply concatenates the outputs
of all classifiers and then applies the softmax. The objective is to minimize
the cross entropy loss between the outputs of the final layer and the true
target over all training samples:

$\mathcal{L}_{\theta_{i}}(x,y)=-\sum_{j}^{C}t_{j}\log(S(c_{0},\ldots,c_{i})_{j})\\\
$

where

$\displaystyle c_{k}$
$\displaystyle=\mathbf{C}_{k}(\mathbf{F}(x))\quad\text{for}\quad k=0,\ldots,j$
$\displaystyle\theta_{i}$ $\displaystyle-\text{ parameters of the last
classifier, }\mathbf{C}_{i}$ $\displaystyle t$ $\displaystyle-\text{the one-
hot encoded vector of target $y$}$ $\displaystyle C$
$\displaystyle=\sum_{k=1}^{i}|\mathcal{Y}_{k}|\quad\text{(i.e. the number of
classes seen so far)}$

Note that only the last classifier is trained, while the loss is a function of
the outputs of all classifiers. In this way, the last classifier is adjusting
itself to respect the prediction of previous classifiers. In other words, it
is learning to produce higher output values for the samples belonging to its
feature space (the new classes in the current training session) while
suppressing itself for the samples belonging to previous classifiers. Although
all of the classifiers are in a sense “partial”, i.e. they can only predict
the classes belonging to them, they also see previous classes during the
training as “negative” examples acting like a regularizer which makes them
even stronger and robust. We also emphasize that since the previous
classifiers are frozen, it eliminates the problem of “forgetting” for them.

###  III-A Memory buffer

As the memory should be kept updated to include new samples, a sampling
strategy has to decide which examples should be selected or removed. We adopt
the greedy sampling approach described in [20]. It randomly replaces some of
the old examples in the memory with the new ones, while trying to satisfy the
balancing constraint, that is, to maintain an equal number of examples for
each class.

###  III-B Training

Training consists of multiple sessions. At the beginning of each training
session, a classifier is constructed which has the output size equal to the
number of classes in the training data. Training objective is to learn the
parameters of this new classifier. We train it by mini-batch gradient descent
where each mini-batch contains an equal number of samples from the current
dataset and the memory buffer. At the end of each training session the memory
buffer is updated with new samples.

###  III-C Early Stopping

As the number of examples per class in the memory keeps decreasing over time,
the model becomes prone to overfitting. Therefore, an early stopping mechanism
is essential. We propose that a portion of the data should be held for
validation. Apart from splitting the incoming data into train and validation
parts, we also dedicate a part of the memory for validation samples. We make
sure that the part of the memory for training gets updated only with the
samples in the training set, and similarly the validation part of the memory
gets updated only with the samples in the validation set.

###  III-D Bias correction layer

Since the mini-batches contain nonequal number of samples from the old and new
classes, it poses a class imbalance problem. This problem has been first
addressed in [26] and the authors proposed a simple yet effective solution,
called Bias Correction (BiC). It is a layer containing only two parameters and
applies a linear transformation to the output logits belonging to new classes
$\mathcal{Y}_{i}$ as follows:

$\displaystyle q_{k}=\left\\{\begin{array}[]{lr}\alpha
o_{k}+\beta&k\in\mathcal{Y}_{i}\\\ o_{k}&otherwise\end{array}\right.$

where $\alpha$ and $\beta$ are the bias parameters and $o_{k}$ is the output
logits of the final layer. The BiC layer is trained separately at the end of
each training session with a small amount of data as it contains only two
parameters.

##  IV Experimental Results

We consider the class-incremental learning scenario on CIFAR100 dataset by
splitting it into 5, 10, and 20 disjoint parts, each part containing 20, 10,
and 5 classes respectively. In all of our experiments, we use the same class
ordering which is obtained by random shuffling. We also run experiments with
different memory sizes. We compare our method against two state-of-the-art
baselines: ER equipped with BiC, and GDumb.

###  IV-A Implementation details

For feature extraction, we use a EfficientNet-B0 network[37] that is
pretrained on ImageNet. We remove from it the final convolution layer, the
output layer and the final MBConv block. Since this feature extractor accepts
inputs of size $224\times 224$, we resize the CIFAR100 images which are
$32\times 32$ to match the input size. The feature extractor is kept frozen in
all experiments of our method and the baselines. For classifiers, we use a
single $1\times 1$ convolution layer followed by global average pooling, and a
dense layer. In all experiments the classifiers have the same architecture
except the number of filters in the convolution layer and the number of units
in the final dense layer. Overall, we have made sure that our method does not
use more trainable parameters in total than the baselines (See table II).

In all experiments, we stop a training session when the validation loss does
not improve for 10 consecutive epochs. We have dedicated 10% of the data for
validation. For all methods we start by learning rate of 0.01. For our method
we decrease the learning rate when the validation loss does not improve for 3
consecutive epochs. For GDumb and ER we apply exponential decay to learning
with rate 0.95. At the end of each training session, just before testing, we
train the BiC layer to remove the prediction bias towards later classes. Since
the whole purpose of BiC is to remove prediction bias, we will train on it
only on the validation part of the memory buffer(after it has been updated to
include all seen classes) which the model itself has never trained on.

$\displaystyle{0}$$\displaystyle{20}$$\displaystyle{40}$$\displaystyle{60}$$\displaystyle{80}$$\displaystyle{100}$Number
of
classes$\displaystyle{0}$$\displaystyle{20}$$\displaystyle{40}$$\displaystyle{60}$$\displaystyle{80}$$\displaystyle{100}$Accuracy
(%)Ours 2000ER w/ BiC 2000GDumb 2000 Figure 2: Accuracies at the end of each
training session. CIFAR100 split into 5 parts.

$\displaystyle{0}$$\displaystyle{10}$$\displaystyle{20}$$\displaystyle{30}$$\displaystyle{40}$$\displaystyle{50}$$\displaystyle{60}$$\displaystyle{70}$$\displaystyle{80}$$\displaystyle{90}$$\displaystyle{100}$Number
of
classes$\displaystyle{0}$$\displaystyle{20}$$\displaystyle{40}$$\displaystyle{60}$$\displaystyle{80}$$\displaystyle{100}$Accuracy
(%)Ours 500ER w/ BiC 500GDumb 500

$\displaystyle{0}$$\displaystyle{10}$$\displaystyle{20}$$\displaystyle{30}$$\displaystyle{40}$$\displaystyle{50}$$\displaystyle{60}$$\displaystyle{70}$$\displaystyle{80}$$\displaystyle{90}$$\displaystyle{100}$Number
of
classes$\displaystyle{0}$$\displaystyle{20}$$\displaystyle{40}$$\displaystyle{60}$$\displaystyle{80}$$\displaystyle{100}$Accuracy
(%)Ours 1000ER w/ BiC 1000GDumb 1000

$\displaystyle{0}$$\displaystyle{10}$$\displaystyle{20}$$\displaystyle{30}$$\displaystyle{40}$$\displaystyle{50}$$\displaystyle{60}$$\displaystyle{70}$$\displaystyle{80}$$\displaystyle{90}$$\displaystyle{100}$Number
of
classes$\displaystyle{0}$$\displaystyle{20}$$\displaystyle{40}$$\displaystyle{60}$$\displaystyle{80}$$\displaystyle{100}$Accuracy
(%)Ours 2000ER w/ BiC 2000GDumb 2000

Figure 3: Accuracies at the end of each training session. CIFAR100 split into
10.

$\displaystyle{0}$$\displaystyle{5}$$\displaystyle{10}$$\displaystyle{15}$$\displaystyle{20}$$\displaystyle{25}$$\displaystyle{30}$$\displaystyle{35}$$\displaystyle{40}$$\displaystyle{45}$$\displaystyle{50}$$\displaystyle{55}$$\displaystyle{60}$$\displaystyle{65}$$\displaystyle{70}$$\displaystyle{75}$$\displaystyle{80}$$\displaystyle{85}$$\displaystyle{90}$$\displaystyle{95}$$\displaystyle{100}$Number
of
classes$\displaystyle{0}$$\displaystyle{20}$$\displaystyle{40}$$\displaystyle{60}$$\displaystyle{80}$$\displaystyle{100}$Accuracy
(%)Ours 1000ER w/ BiC 1000GDumb 1000

$\displaystyle{0}$$\displaystyle{5}$$\displaystyle{10}$$\displaystyle{15}$$\displaystyle{20}$$\displaystyle{25}$$\displaystyle{30}$$\displaystyle{35}$$\displaystyle{40}$$\displaystyle{45}$$\displaystyle{50}$$\displaystyle{55}$$\displaystyle{60}$$\displaystyle{65}$$\displaystyle{70}$$\displaystyle{75}$$\displaystyle{80}$$\displaystyle{85}$$\displaystyle{90}$$\displaystyle{95}$$\displaystyle{100}$Number
of
classes$\displaystyle{0}$$\displaystyle{20}$$\displaystyle{40}$$\displaystyle{60}$$\displaystyle{80}$$\displaystyle{100}$Accuracy
(%)Ours 2000ER w/ BiC 2000GDumb 2000

Figure 4: Accuracies at the end of each training session. CIFAR100 split into 20 parts. TABLE I: Accuracies at the end of the training. CIFAR100 split into 5, 10, and 20 parts. Splits | 5 | 10 | 20
---|---|---|---
Memory | 2000 | 500 | 1000 | 2000 | 1000 | 2000
GDumb | 48.24 | 23.79 | 39.89 | 48.74 | 39.57 | 47.83
ER w/ BiC | 57.89 | 42.46 | 48.97 | 54.68 | 48.78 | 55.23
Ours | 67.83 | 56.45 | 63.48 | 65.79 | 61.60 | 63.25
TABLE II: Number of trainable parameters in each of the experiment settings, when training on CIFAR100 split into 5, 10, and 20 parts. Method | # trainable parameters
---|---
5 splits | 10 splits | 20 splits
GDumb | 118K | 118K | 118K
ER | 118K | 118K | 118K
Ours | 17K $\to$ 86K | 8K $\to$ 82K | 4K $\to$ 80K

###  IV-B Baselines

GDumb and ER use the same architecture, they only differ in training. Gdumb
only trains on memory samples. It updates the memory buffer at the beginning
of each training session and then trains only on the memory buffer discarding
the rest of the data. We also train a BiC layer at the end of each training
session of ER. BiC layer is not necessary for GDumb due to the fact that it
trains on balanced dataset, i.e., the memory buffer.

###  IV-C Analysis

We test the models at the end of each training session over all classes that
have been encountered so far. In table I we report the accuracy at the end of
the training in three different settings with varying memory buffer sizes
(500, 1000, and 2000) and dataset splitting into 5 (Fig. 2), 10 (Fig. 3), and
20 (Fig. 4) parts. In all settings we see a large gap in accuracy between our
method and the second best method, ER. ER always outperforms GDumbs, as
expected, because it trains on all available data and the early stopping
mechanism that we introduced here prevents it from unnecessary training which
causes more forgetting of past data.

We observe that the gap in accuracy tends to get larger when we shrink the
memory buffer size. Our method still reaches a remarkable performance of
56.45% when memory size is 500 leaving a gap of 14% relative to the next
method, ER. On the other hand, we see a large performance drop of GDumb. This
is mainly because it is highly dependent on the memory size as it trains only
on the memory samples. However, GDumb shows the least difference in
performance when training with different dataset splittings, meaning that it
might have an advantage in settings more close to online learning.

###  IV-D Ablation studies

We conduct experiments to see if all components of our method are indeed
important. The first component is that we freeze the classifiers after we
train them so that in this way we believe the problem of “forgetting” is
eliminated. Therefore we have conducted an experiment where we let our method
to keep updating all classifiers. The second component is the additional BiC
layer. Looking at the confusion matrices, we have observed a prediction bias
towards the last group of classes which was a signal to incorporate BiC layer.
We also present the results here where our method does not use a BiC layer.
These experiments are conducted in the case of CIFAR100 split into 10 classes,
and a memory buffer size of 2000. The reported results are in Table III. We
can see the benefit of using BiC layer, which confirms that our method would
have some prediction bias without it. On the other hand, we can see a large
drop in accuracy (65.79% $\to$ 57.04%) when we allow updating our classifiers.
Nevertheless, still the performance stays above the other two methods (looking
at Table I, 54.68% and 48.74% of ER and GDumb, respectively).

TABLE III: Ablation studies testing our method in two alternative forms: (1) without freezing the classifiers, and (2) without using BiC layer. Method | Accuracy %
---|---
Ours w/out freezing | 57.04
Ours w/out BiC | 62.86
Ours | 65.79

##  V Conclusion

In this paper, we proposed a new approach for class-incremental learning. We
tackled the problem by training a separate classifier for each new group of
classes. By freezing these classifiers after they have been trained, we have
limited the problem of “forgetting” and achieved big improvements over strong
SOTA baselines. Our method consistently achieved better performance when
tested on CIFAR100 learning different number of classes at a time.

## References

  * [1] A. Khan and K. Fraz, “Post-training iterative hierarchical data augmentation for deep networks,” _Advances in Neural Information Processing Systems_ , vol. 33, 2020.
  * [2] A. R. Rivera, A. Khan, I. E. I. Bekkouch, and T. S. Sheikh, “Anomaly detection based on zero-shot outlier synthesis and hierarchical feature distillation,” _IEEE Transactions on Neural Networks and Learning Systems_ , 2020.
  * [3] K. Yakovlev, I. E. I. Bekkouch, A. M. Khan, and A. M. Khattak, “Abstraction-based outlier detection for image data,” in _Proceedings of SAI Intelligent Systems Conference_. Springer, 2020, pp. 540–552.
  * [4] E. Batanina, I. E. I. Bekkouch, Y. Youssry, A. Khan, A. M. Khattak, and M. Bortnikov, “Domain adaptation for car accident detection in videos,” in _2019 Ninth International Conference on Image Processing Theory, Tools and Applications (IPTA)_. IEEE, 2019, pp. 1–6.
  * [5] M. Bortnikov, A. Khan, A. M. Khattak, and M. Ahmad, “Accident recognition via 3d cnns for automated traffic monitoring in smart cities,” in _Science and Information Conference_. Springer, 2019, pp. 256–264.
  * [6] Y. Gavrilin and A. Khan, “Across-sensor feature learning for energy-efficient activity recognition on mobile devices,” in _2019 International Joint Conference on Neural Networks (IJCNN)_. IEEE, 2019, pp. 1–7.
  * [7] K. Sozykin, S. Protasov, A. Khan, R. Hussain, and J. Lee, “Multi-label class-imbalanced action recognition in hockey videos via 3d convolutional neural networks,” in _2018 19th IEEE/ACIS International Conference on Software Engineering, Artificial Intelligence, Networking and Parallel/Distributed Computing (SNPD)_. IEEE, 2018, pp. 146–151.
  * [8] A. M. Khan, Y.-K. Lee, S. Lee, and T.-S. Kim, “Accelerometer’s position independent physical activity recognition system for long-term activity monitoring in the elderly,” _Medical & biological engineering & computing_, vol. 48, no. 12, pp. 1271–1279, 2010.
  * [9] S. Protasov, A. M. Khan, K. Sozykin, and M. Ahmad, “Using deep features for video scene detection and annotation,” _Signal, Image and Video Processing_ , vol. 12, no. 5, pp. 991–999, 2018.
  * [10] M. Ahmad, A. M. Khan, M. Mazzara, S. Distefano, M. Ali, and M. S. Sarfraz, “A fast and compact 3-d cnn for hyperspectral image classification,” _IEEE Geoscience and Remote Sensing Letters_ , 2020.
  * [11] M. Ahmad, A. M. Khan, M. Mazzara, and S. Distefano, “Multi-layer extreme learning machine-based autoencoder for hyperspectral image classification.” in _VISIGRAPP (4: VISAPP)_ , 2019, pp. 75–82.
  * [12] M. Gusarev, R. Kuleev, A. Khan, A. R. Rivera, and A. M. Khattak, “Deep learning models for bone suppression in chest radiographs,” in _2017 IEEE Conference on Computational Intelligence in Bioinformatics and Computational Biology (CIBCB)_. IEEE, 2017, pp. 1–7.
  * [13] A. Dobrenkii, R. Kuleev, A. Khan, A. R. Rivera, and A. M. Khattak, “Large residual multiple view 3d cnn for false positive reduction in pulmonary nodule detection,” in _2017 IEEE conference on computational intelligence in bioinformatics and computational biology (CIBCB)_. IEEE, 2017, pp. 1–6.
  * [14] A. Khusainova, A. Khan, and A. R. Rivera, “Sart-similarity, analogies, and relatedness for tatar language: New benchmark datasets for word embeddings evaluation,” _arXiv preprint arXiv:1904.00365_ , 2019.
  * [15] A. Valeev, I. Gibadullin, A. Khusainova, and A. Khan, “Application of low-resource machine translation techniques to russian-tatar language pair,” _arXiv preprint arXiv:1910.00368_ , 2019.
  * [16] R. Kemker, M. McClure, A. Abitino, T. Hayes, and C. Kanan, “Measuring catastrophic forgetting in neural networks,” _arXiv preprint arXiv:1708.02072_ , 2017.
  * [17] M. McCloskey and N. J. Cohen, “Catastrophic interference in connectionist networks: The sequential learning problem,” vol. 24, pp. 109–165, 1989.
  * [18] G. M. van de Ven and A. S. Tolias, “Three scenarios for continual learning,” 2019\.
  * [19] P. Buzzega, M. Boschini, A. Porrello, and S. Calderara, “Rethinking experience replay: a bag of tricks for continual learning,” _arXiv preprint arXiv:2010.05595_ , 2020.
  * [20] A. Prabhu, P. Torr, and P. Dokania, “Gdumb: A simple approach that questions our progress in continual learning,” in _The European Conference on Computer Vision (ECCV)_ , August 2020.
  * [21] S. Rebuffi, A. Kolesnikov, G. Sperl, and C. H. Lampert, “icarl: Incremental classifier and representation learning,” in _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2017, pp. 5533–5542.
  * [22] D. Lopez-Paz and M. Ranzato, “Gradient episodic memory for continual learning,” in _Advances in neural information processing systems_ , 2017, pp. 6467–6476.
  * [23] T. L. Hayes, K. Kafle, R. Shrestha, M. Acharya, and C. Kanan, “Remind your neural network to prevent catastrophic forgetting,” in _Proceedings of the European Conference on Computer Vision (ECCV)_ , 2020.
  * [24] V. Lomonaco and D. Maltoni, “Core50: a new dataset and benchmark for continuous object recognition,” 2017.
  * [25] R. Ratcliff, “Connectionist models of recognition memory: constraints imposed by learning and forgetting functions.” _Psychological review_ , vol. 97, no. 2, p. 285, 1990.
  * [26] Y. Wu, Y. Chen, L. Wang, Y. Ye, Z. Liu, Y. Guo, and Y. Fu, “Large scale incremental learning,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 374–382.
  * [27] A. Chaudhry, M. Ranzato, M. Rohrbach, and M. Elhoseiny, “Efficient lifelong learning with a-gem,” _arXiv preprint arXiv:1812.00420_ , 2018.
  * [28] A. Chaudhry, A. Gordo, P. K. Dokania, P. Torr, and D. Lopez-Paz, “Using hindsight to anchor past knowledge in continual learning,” _arXiv preprint arXiv:2002.08165_ , 2020.
  * [29] Z. Li and D. Hoiem, “Learning without forgetting,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 40, no. 12, pp. 2935–2947, 2017\.
  * [30] J. Kirkpatrick, R. Pascanu, N. Rabinowitz, J. Veness, G. Desjardins, A. A. Rusu, K. Milan, J. Quan, T. Ramalho, A. Grabska-Barwinska _et al._ , “Overcoming catastrophic forgetting in neural networks,” _Proceedings of the national academy of sciences_ , vol. 114, no. 13, pp. 3521–3526, 2017.
  * [31] F. Zenke, B. Poole, and S. Ganguli, “Continual learning through synaptic intelligence,” _Proceedings of machine learning research_ , vol. 70, p. 3987, 2017.
  * [32] T. Lesort, A. Stoian, and D. Filliat, “Regularization shortcomings for continual learning,” 2020.
  * [33] A. A. Rusu, N. C. Rabinowitz, G. Desjardins, H. Soyer, J. Kirkpatrick, K. Kavukcuoglu, R. Pascanu, and R. Hadsell, “Progressive neural networks,” _arXiv preprint arXiv:1606.04671_ , 2016.
  * [34] C.-Y. Hung, C.-H. Tu, C.-E. Wu, C.-H. Chen, Y.-M. Chan, and C.-S. Chen, “Compacting, picking and growing for unforgetting continual learning,” in _Advances in Neural Information Processing Systems_ , 2019, pp. 13 669–13 679.
  * [35] R. Aljundi, P. Chakravarty, and T. Tuytelaars, “Expert gate: Lifelong learning with a network of experts,” 2017.
  * [36] D. Maltoni and V. Lomonaco, “Continuous learning in single-incremental-task scenarios,” _Neural Networks_ , vol. 116, pp. 56–73, 2019.
  * [37] M. Tan and Q. Le, “Efficientnet: Rethinking model scaling for convolutional neural networks,” in _International Conference on Machine Learning_. PMLR, 2019, pp. 6105–6114.


# GM-TOuNN: Graded Multiscale Topology Optimization using Neural Networks

Aaditya Chandrasekhar , Saketh Sridhara ††footnotemark: , Krishnan Suresh
Department of Mechanical Engineering, University of Wisconsin-Madison
{achandrasek3, ssridhara<EMAIL_ADDRESS>Contributed equally

###### Abstract

Multiscale topology optimization (M-TO) entails generating an optimal global
topology, and an optimal set of microstructures at a smaller scale, for a
physics-constrained problem. With the advent of additive manufacturing, M-TO
has gained significant prominence. However, generating optimal microstructures
at various locations can be computationally very expensive. As an alternate,
_graded multiscale topology optimization_ (GM-TO) has been proposed where one
or more pre-selected and graded (parameterized) microstructural topologies are
used to fill the domain optimally. This leads to a significant reduction in
computation while retaining many of the benefits of M-TO.

A successful GM-TO framework must: (1) be capable of efficiently handling
numerous pre-selected microstructures, (2) be able to continuously switch
between these microstructures during optimization, (3) ensure that the
partition of unity is satisfied, and (4) discourage microstructure mixing at
termination.

In this paper, we propose to meet these requirements by exploiting the unique
classification capacity of neural networks. Specifically, we propose a _graded
multiscale topology optimization using neural-network_ (GM-TOuNN) framework
with the following features: (1) the number of design variables is only weakly
dependent on the number of pre-selected microstructures, (2) it guarantees
partition of unity while discouraging microstructure mixing, and (3) it
supports automatic differentiation, thereby eliminating manual sensitivity
analysis. The proposed framework is illustrated through several examples.

Overview of the proposed method: Given candidate microstructures and a
topology optimization problem, a neural network (NN) maps the spatial
coordinates to a volume fraction field and a microstructure field, that
produce a graded multiscale design.

_K_ eywords multiscale topology optimization, graded microstructure, neural
networks, automatic differentiation

##  1 Introduction

Topology optimization (TO) is a strategy for optimally distributing material
within a design domain ([1, 2]) to maximize a desired objective, while meeting
one or more constraints. Various TO methods have been proposed; these include
density-based methods([3]), homogenization ([4, 5]), level-set ([6, 7]),
evolutionary ([8, 9]), morphable components ([10]), and topological
sensitivity based methods ([11, 12]). These methods typically lead to a
single-scale design where all geometric features are of similar size. As an
example, for the cantilever problem posed in fig. 1(a), a single-scale TO that
minimizes the compliance, subject to a volume constraint, leads to the
topology in fig. 1(b).

Figure 1: TO problem and solutions: (a) Mid-cantilever problem. (b) Single-
scale solution. (c) Multiscale solution. (d) Graded multiscale solution using
a single X-shaped microstructure.

In contrast, in multiscale topology optimization (M-TO), geometric features
are computed at two or more length scales ([13],[14]). For the problem in fig.
1(a), fig. 1(c) illustrates a topology, based on a two-scale M-TO ([15]). Such
designs exhibit unique characteristics such as large surface-to-volume ratio,
excellent resilience, etc. This has led to a wide range of applications ([16,
17]) including energy absorption ([18]), heat exchangers ([19]), and
biomedical applications ([20]). Current M-TO methods include clustering ([21,
15, 22]), kriging ([23]), multi-material ([24]), conformal lattices ([25]),
etc. The advent of additive manufacturing has further spurned research in M-TO
([16]).

###  1.1 Graded Multiscale TO

However, one of the challenges in M-TO is the high computational cost since
one must generate and evaluate various microstructures (through
homogenization) during each step of the optimization process ([13, 21, 15,
26]). To reduce this cost, researchers have proposed the use of graded
variations of one or more pre-selected microstructural topologies ([14, 19,
27]). [14] expressed the mechanical properties as a density function, with a
B-spline-based interpolation. [19] designed graded cellular structures With
triply periodic shapes. [27] parameterized the radii of each member of a
lattice with linear and sinusoidal graded variations. These methods are
referred to as _graded multiscale TO_ , or GM-TO, leading to topologies such
as the one illustrated in fig. 1(d), where, as an example, graded variations
of a single X-shaped microstructure are used everywhere. This leads to a
significant reduction in computational cost ([28]) since homogenization is not
needed within the optimization loop. However, restricting the design to a
single microstructure can significantly reduce performance ([29]). To further
improve the performance of GM-TO, without substantially increasing the
computational cost, one can use _a finite number of graded microstructures_.

Density-based methods: Theoretically, using a finite number of microstructures
in GM-TO is analogous to using multiple materials in TO. Therefore, various
authors ([30, 31, 32, 33]) proposed the use of uniform multi-phase materials
interpolation (UMMI) models for GM-TO. Given the elasticity matrix of each
microstructure at a particular volume fraction, two different UMMI models have
been proposed for interpolation [34, 35]. The first type (UMMI-1) proposed by
[31, 33], simply adds the contributions from each microstructure, and does not
penalize microstructure mixing. The UMMI of the second type (UMMI-2) penalizes
microstructure mixing at the cost of non-linearity ([30]). To mitigate the
effect of non-linearity, [36] suggested a gradual increase in the penalty in
multi-materials, further explored as a multi-scale formulation by [37, 38].
Despite penalization, microstructure mixing is not entirely eliminated in
UMMI-2. Further, while multiple materials can often co-exist at a given point,
this is unacceptable in GM-TO. Finally, a handful (typically less than five)
materials can be sufficient in multi-material TO, while a large number of
microstructures are required in GM-TO.

Level-set methods: Level-set methods are also a popular choice for GM-TO. [39,
40] considered a level-set based method where the microstructural shape is
represented and parametrized by implicit functions thereby circumventing the
need for homogenization within every loop; instead one can rely on curve-
fitting. [41] utilized the level set function to evolve on topologies on both
the micro and macroscale, with connectivity ensured by the higher-order
continuity of the level-set function. [42] addressed the challenge of
connectivity of microstructures through shape metamorphosis to build graded
transition zones using a connectivity index. Finally, GM-TO with graded
variations of a single microstructure was implemented by [43] based on the
parametrized level-set functions using radial basis functions.

Data-driven methods: Recently, data-driven techniques have been proposed to
address GM-TO. For instance, [44, 45] utilize a latent variable Gaussian
process that embeds discrete microstructures in a continuous and
differentiable latent space. While the latent space is trained on the discrete
homogenized data, the procedure of snapping and re-optimization may lead to
sub-optimal results. Several authors ([46, 47, 48]) replaced the expensive
homogenization process with neural networks that map the design variables to
the homogenized elasticity matrix. [33] used neural networks to interpolate
the stiffness matrices from homogenized data. [49] used two neural networks
trained on topology optimized data: one for determining the microstructural
topology and the other to improve the connectivity of microstructures. One of
the challenges in data-driven methods is the computational expense in data
generation and training of the machine learning models. Another challenge with
the data-driven methods is that physical bounds on the elasticity matrices
(such as positive-definiteness) cannot be enforced explicitly, which can lead
to spurious stiffness matrices and non-convergence.

###  1.2 Paper Contributions

In this paper, we exploit the unique classification capability of neural
networks to address some of the limitations of current GM-TO methods.
Specifically, we propose a _graded multiscale topology optimization using
neural-networks_ (GM-TOuNN) (section 2) framework. GM-TOuNN extends the mesh-
independent neural-network (NN) based representation of the macro-scale
topology proposed by [50] to GM-TO. The mesh-independent representation allows
us to consider a large number of candidate microstructures without increasing
the number of design variables. Further, the partition of unity is implicitly
guaranteed by the NN construction. Finally, one can leverage automatic
differentiation (AD) of the NN computational framework to provide end-to-end
differentiability and automated sensitivity analysis. Numerical experiments
are presented in section 3. We conclude with limitations and future work in
section 4.

##  2 Proposed Method

###  2.1 Problem Specification

Consider a design domain with prescribed loads, restraints and a set of
microstructural topologies (fig. 2). In the current work, we seek to compute
an optimal multiscale design where we determine, at every location, the
appropriate microstructure and its size-parameter (gradation). For simplicity,
we will assume the compliance must be minimized subject to a volume
constraint.

Figure 2: A GM-TO problem.

To assist in optimization, we introduce the following design variables. The
presence or absence of a microstructure $m$ at any point $\bm{x}$ will be
denoted by the variable $\rho_{m}(\bm{x})$, where, ideally,
$\rho_{m}(\bm{x})\in\\{0,1\\}$. However, for continuous gradient-based
optimization, we will let $0\leq\rho_{m}\leq 1$ and drive it towards $0/1$
through penalization. Thus, at any point $\bm{x}$, one can define the vector
$\bm{\rho}(\bm{x})=\\{\rho_{1}(\bm{x}),\rho_{2}(\bm{x}),\ldots,\rho_{M}(\bm{x})\\}$
that captures the presence or absence of the $M$ microstructures with the
_partition of unity_ constraint $\sum\rho_{m}(\bm{x})=1$. Further, since the
microstructures are graded, we control their size by the scalar design
variable $0\leq v(\bm{x})\leq 1$ where 0 denotes a void and 1 denotes a
complete fill. Thus, in conclusion, $\bm{\rho}(\bm{x})$ dictates the type of
microstructure, while $v(\bm{x})$ dictates the size or volume fraction of the
microstructure at $\bm{x}$.

Consequently, one can pose the GM-TO problem in a discrete finite-element
setting as:

$\displaystyle\underset{\bm{\rho},\bm{v}}{\text{minimize}}$ $\displaystyle
J(\bm{\rho},\bm{v})=\bm{f}^{T}\bm{u}$ (1a) subject to
$\displaystyle\bm{K}(\bm{\rho},\bm{v})\bm{u}=\bm{f}$ (1b) $\displaystyle
g_{v}(\bm{v})\coloneqq\frac{\sum\limits_{e}v(\bm{x}_{e})A_{e}}{V_{f}^{*}\sum\limits_{e}A_{e}}-1\leq
0$ (1c)
$\displaystyle\sum\limits_{m=1}^{M}\rho_{m}(\bm{x})=1\;,\;\forall\bm{x}$ (1d)
$\displaystyle 0\leq\rho_{m}(\bm{x})\leq 1\;,\;\forall\bm{x}\;,\;\forall m$
(1e) $\displaystyle 0\leq v(\bm{x})\leq 1\;,\;\forall\bm{x}$ (1f)

where $J$ is the structural compliance, $\bm{K}$ is the finite element
stiffness matrix, $\bm{u}$ is the displacement field, $\bm{f}$ is the applied
force, $V_{f}^{*}$ the maximum allowed volume fraction, with $A_{e}$ is the
area (2D) of the element $e$.

###  2.2 Design Representation using Neural Networks

Typically, for such problems, the design variables are captured via the
underlying FE mesh ([51]), i.e., for the above problem, design variables
$\bm{\rho}_{e}$ and $v_{e}$ are defined at each element $e$. Thus the number
of design variables grows linearly with the mesh size. Further, observe that
the partition of unity constraint must be imposed over each element.

In this paper, we avoid this undesirable growth in complexity by indirectly
controlling the design variables via a coordinate-based neural-network ([50,
52]). The proposed neural-network (NN) architecture (see fig. 3) consists of
the following entities:

  1. 1.

Input Layer: The input to the NN are points $\bm{x}\in\mathbf{R}^{d}$ ($d=2$
for 2D) within the domain.

  2. 2.

Hidden Layers: The hidden layers consists of a series of Swish
($x.sigmoid(x)$, [53]) activated dense layers.

  3. 3.

Output Layer: The output layer consists of $M+1$ neurons corresponding to the
microstructure composition $\bm{\rho}=\\{\rho_{1},\rho_{2},\ldots\rho_{M}\\}$
and the volume fraction $v$. Further, the neurons associated with $\bm{\rho}$
are activated by a softmax function. This guarantees that the partition of
unity ($\sum\rho_{i}=1$) and physical validity $0\leq\rho_{i}\leq 1$ are
automatically satisfied. The output neuron associated with volume fraction is
activated via a Sigmoid function, ensuring that $0\leq v(\bm{x})\leq 1$. Thus,
no additional constraint are needed.

  4. 4.

Design Variables: The weights and bias, denoted by the $\bm{w}$, now become
the primary design variables, i.e., we have $\bm{\rho}(x,y;\bm{w})$ and
$v(x,y;\bm{w})$.

Figure 3: Neural network architecture used for GM-TOuNN.

Thus the strategy is to perform GM-TO via the NN weights $\bm{w}$, i.e., the
GM-TO problem in eq. 1 reduces to:

$\displaystyle\underset{\bm{w}}{\text{minimize}}$ $\displaystyle
J(\bm{w})=\bm{f}^{T}\bm{u}(\bm{w})$ (2a) subject to
$\displaystyle\bm{K}(\bm{w})\bm{u}=\bm{f}$ (2b) $\displaystyle
g_{v}(\bm{w})\coloneqq\frac{\sum\limits_{e}v(\bm{x}_{e};\bm{w})A_{e}}{V_{f}^{*}\sum\limits_{e}A_{e}}-1\leq
0$ (2c)

Observe that: _(1) no additional constraint is needed since they are
automatically satisfied by the NN, and (2) increasing the number of candidate
microstructures only increases the number of output neurons but not the number
of design variables._

###  2.3 Material Model

Given the NN-architecture, one can now proceed to construct the material model
for analysis. Towards this end, let $[C_{m}(v)]$ be the elasticity matrix of
microstructure $m$ at volume fraction $v$, where $[C_{m}(v)]$ consists of six
independent components (in 2D):

$[C_{m}(v)]=\begin{bmatrix}C_{11}(v)&C_{12}(v)&C_{13}(v)\\\
&C_{22}(v)&C_{23}(v)\\\
\lx@intercol\hfil\text{sym.}\hfil\lx@intercol&C_{33}(v)\end{bmatrix}$ (3)

To obtain $[C_{m}(v)]$, we adopt a simple constrained polynomial scheme
([54]). Specifically, the homogenized constitutive matrix [55] for
microstructure $m$ is evaluated at a few instances of volume fractions. Then,
a polynomial ensuring positive definiteness is interpolated to these instances
for each of the components. As an example, the polynomials for an X-type
microstructure are illustrated in fig. 4.

Figure 4: Polynomial fit of the components of constitutive matrix for the $X$
type microstructure (note: $C_{00}=C_{11}$)

Given $[C_{m}(v)]$ for each microstructure $m$, the effective elasticity
matrix at any location $\bm{x}$ is defined here as a weighted average ([34]):

$[C(\bm{\rho},v)]=\sum\limits_{m=1}^{M}\rho_{m}^{p}[C_{m}(v)]$ (4)

where $p$ is the solid isotropic material with penalization (SIMP) constant.
The penalization discourages intermediate volume fractions, i.e., discourages
microstructure mixing.

###  2.4 Finite Element Analysis

We will use conventional finite element analysis (FEA) as part of the
framework. Here, we use a structured quadrilateral mesh due to its simplicity.
The element stiffness matrix is defined as:

$[K_{e}]=\int\limits_{\Omega_{e}}[\nabla N_{e}]^{T}[C(\bm{x}_{e})][\nabla
N_{e}]d\Omega_{e}$ (5)

where $[\nabla N_{e}]$ is the gradient of the shape matrix, and
$[C(\bm{x}_{e})]$ is the elasticity tensor evaluated at the center of the
element.

In a straightforward, but naïve, implementation, one would evaluate
$\bm{\rho}$ and $v$ at the center of each element via the NN. Then, given the
pre-computed polynomials $[C_{m}(v)]$, one would then find the effective $[C]$
for that element via eq. 4. Finally the element stiffness matrices will be
computed via eq. 5. Thus the element stiffness matrices must be computed for
each element, and for each step of the optimization process. This can be
computationally very expensive. To significantly reduce the computation, we
exploit the concept of _template stiffness matrices_ proposed by [15].

Observe that one can express the element stiffness matrix as follows:

$\displaystyle=C_{11}[\hat{K}^{1}]+C_{22}[\hat{K}^{2}]+C_{33}[\hat{K}^{3}]$
(6)
$\displaystyle+C_{12}[\hat{K}^{4}]+C_{13}[\hat{K}^{5}]+C_{23}[\hat{K}^{6}]$

where,

$[\hat{K}^{i}]=\int\limits_{\Omega_{e}}[\nabla N_{e}]^{T}[\hat{C}^{i}][\nabla
N_{e}]d\Omega_{e},\quad i=1,2,\ldots 6$ (7)

and

$\displaystyle[\hat{C}^{1}]=\begin{bmatrix}1&0&0\\\ 0&0&0\\\ 0&0&0\\\
\end{bmatrix}\;[\hat{C}^{2}]=\begin{bmatrix}0&0&0\\\ 0&1&0\\\ 0&0&0\\\
\end{bmatrix}\quad\;[\hat{C}^{3}]=\begin{bmatrix}0&0&0\\\ 0&0&0\\\ 0&0&1\\\
\end{bmatrix}$ $\displaystyle[\hat{C}^{4}]=\begin{bmatrix}0&1&0\\\ 1&0&0\\\
0&0&0\\\ \end{bmatrix}\;[\hat{C}^{5}]=\begin{bmatrix}0&0&1\\\ 0&0&0\\\
1&0&0\\\ \end{bmatrix}\quad\;[\hat{C}^{6}]=\begin{bmatrix}0&0&0\\\ 0&0&1\\\
0&1&0\\\ \end{bmatrix}$ (8)

Thus the matrices $[\hat{K}^{i}]$ are computed at the start of the
optimization. Then, during optimization, the element stiffness matrices are
computed efficiently via eq. 6. In other words, one would evaluate $\bm{\rho}$
and $v$ at the center of each element via the NN. Then, given the pre-computed
polynomials $[C_{m}(v)]$, one would then find the effective $[C]$ for that
element (as before). But, to compute the element stiffness matrices, we will
rely on eq. 6 This will be followed by the assembly of the global stiffness
matrix $\bm{K}$.

###  2.5 Loss Function

We now consider solving the NN-based optimization problem in eq. 2. Neural
networks are designed to minimize a loss function using well-known
optimization techniques such as Adam procedure ([56]). We therefore convert
the constrained minimization problem into a loss function minimization by
employing a log-barrier scheme as proposed in [57]. Specifically, the loss
function is defined as

$L(\bm{w})=J+\psi(g_{v})$ (9)

where,

$\psi_{t}(g)=\begin{cases}-\frac{1}{t}\log(-g)\;,\quad
g\leq\frac{-1}{t^{2}}\\\
tg-\frac{1}{t}\log(\frac{1}{t^{2}})+\frac{1}{t}\;,\quad\text{otherwise}\end{cases}$
(10)

where the parameter $t$ is updated during each iteration, making the
enforcement of the constraint stricter as the optimization progresses.

###  2.6 Sensitivity Analysis

Since Adam is a gradient-based optimizer, it requires sensitivities, i.e.,
derivative of the loss function (eq. 9) with respect to the design variable
(weights of the NN $\bm{w}$). Fortunately, one can exploit modern automatic
differentiation frameworks ( [58]) to avoid manual sensitivity calculations.
In particular, expressing all our computation, including FEA within JAX
([59]), results in an end-to-end differentiable framework, as illustrated in
fig. 5.

Figure 5: Optimization loop of the proposed GM-TOuNN framework.

###  2.7 Algorithm

We summarize the proposed framework in algorithm 1. We will assume that a GM-
TO problem, a desired volume fraction, an NN configuration and a set of
microstructures with their interpolated coefficients are given.

The first step in the algorithm is to discretize the domain with a finite
element mesh, and sample the mesh (3) at the center of each element (these
serve as inputs to the NN). We compute the template stiffness matrices (4).
The log barrier penalty $t$ and SIMP penalty parameter are also initialized
(5).

In the main iteration, the element volume fraction $v$ and microstructure
field $\bm{\rho}$ are computed using the NN using the current values of
$\bm{w}$ (7). These fields are then used to construct the stiffness matrix and
to solve for the displacement (8 \- 12). Then we compute the objective (13)
and volume constraint (14), leading to the loss function (15). The weights
$\bm{w}$ are then updated using the Adam optimization scheme (17). The
optimizer requests the sensitivities which are computed in an automated
fashion (16). Finally the log-barrier penalty and SIMP penalty parameters are
updated (19 and 20). The process is repeated until termination, i.e., till the
relative change in loss is below a certain threshold or the iterations exceed
a maximum value. Upon termination, each element is replaced with an image of
the associated microstructure with the desired volume fraction. The framework
is schematically depicted in fig. 5.

Algorithm 1 Graded Mutliscale TO

1:procedure GMTOuNN($\Omega_{0}$, BC, $V_{f}^{*}$, …) $\triangleright$ Inputs

2: $\Omega^{0}\rightarrow\Omega^{0}_{h}$ $\triangleright$ Discretize domain
for FE

3: $\bm{x}=\\{x_{e},y_{e}\\}_{e\in\Omega^{0}_{h}}$ $\triangleright$ elem
centers; NN input

4: $\Omega^{0}_{h}\rightarrow[\hat{K}]$ $\triangleright$ Compute stiff.
templates eq. 7

5: k = 0; $t=t_{0}$; $p=1$ $\triangleright$ Initialization

6: repeat$\triangleright$ Optimization (Training)

7: $NN(\bm{x}_{e};\bm{w})\rightarrow\\{\bm{\rho},v\\}_{e}$ $\triangleright$
Fwd prop NN

8: $v\rightarrow[C_{m}]\;m=1,\ldots,M\;\forall e$$\triangleright$ elasticity
matr. eq. 3

9: $\\{[C_{m}],\rho_{m}\\}_{e}\rightarrow[C]_{e}$ $\triangleright$ net
elasticity matrix eq. 4

10: $\\{[C],[\hat{K}]\\}_{\forall e}\rightarrow[K]_{e}$ $\triangleright$ elem
stiff matrices eq. 6

11: $\bigcup\limits_{e}[K]_{e}\rightarrow[K]$ $\triangleright$ global stiff
matrix

12: $\\{[K],\bm{f}\\}\rightarrow\bm{u}$ $\triangleright$ FEA eq. 2b

13: $\\{\bm{f},\bm{u}\\}\rightarrow J$ $\triangleright$ Objective, Equation 2a

14: $\\{v_{e},V_{f}^{*}\\}\rightarrow g_{v}$ $\triangleright$ vol cons, eq. 2c

15: $\\{J,g_{v}\\}\rightarrow L$ $\triangleright$ Loss from Equation 9

16: $AD(L\leftarrow\bm{w})\rightarrow\nabla L$ $\triangleright$ sensitivity
analysis

17: $w+\Delta w(\nabla L)\rightarrow w$ $\triangleright$ Adam optimizer [56]
step

18: $\text{k}++$

19: $t=t_{0}\mu^{k}$ $\triangleright$ Increment $t$; log-barrier term

20: $p=p+\Delta p$ $\triangleright$ SIMP penalty update

21: until $||\Delta L||<\epsilon^{*}$ and $k<k_{max}$ $\triangleright$ convg.
criteria

22: return $\bm{w},\bm{\rho},\bm{v}$

23:end procedure

##  3 Numerical examples

In this section, we conduct several experiments to illustrate the method and
algorithm. The default settings are as follows:

  * •

Mesh: A mesh size of $60\times 30$ with elements of size $1\times 1$ is used
for all experiments, unless otherwise stated; the force is assumed to be 1
unit.

  * •

Material: The microstructures are assumed to be composed of an isotropic
material with a Poisson’s ratio $\nu=0.3$ and Young’s modulus $E=1$.

  * •

Neural Network: The NN comprises of 4 Swish-activated hidden layers with 20
neurons in each layer.

  * •

Candidate Microstructures: A set of 11 predefined microstructures (fig. 6) are
used in the experiments. Observe that connectivity between these
microstructures is guaranteed. A quintic polynomial is used to interpolate the
components of the constitutive matrix.

  * •

Material Penalization: The SIMP penalization constant $p$ in eq. 4 is
incremented every iteration by $0.02$ using the continuation scheme [60],
starting from a value of $1$, capped at $8$.

  * •

Loss Function: The constraint penalty values of $t_{0}=1$ and $\mu=1.01$ is
used.

  * •

Optimizer: Adam optimizer with a learning rate of $0.01$ is used. Further, the
gradients are clipped at a maximum norm of 1 to improve stability of
convergence ([61]). A maximum of 300 iterations were allowed with
$\epsilon^{*}=0.01$

  * •

Computing Environment: All experiments are conducted on a MacBook M1 Pro,
using the JAX ([59]) environment.

Figure 6: Microstructures considered in the experiments.

###  3.1 Validation

Consider the simple tensile bar problem in fig. 7(a). Although eleven
microstructures were provided (see fig. 6), the final topology in fig. 7(b)
consists of a single microtopology (number 1) with a complete fill for all
elements in the middle and void at the top and bottom, i.e., the result is
consistent with expectations. The optimization converged in 240 iterations,
taking 47 seconds.

(a) Loading of a tensile bar.

(b) GM-TO of a tensile bar with $V_{f}^{*}=0.4$ ($J^{*}=183$).

Figure 7: A tensile-bar problem.

Next, we consider an MBB beam in fig. 8(a). [54] reported compliance values
between $76$ and $88$ when considering only one microstructure at a time. In
comparison, we achieve a compliance of $68$ when all microstructures are
allowed; see fig. 8(b).

(a) Loading of an MBB beam.

(b) Optimized MBB beam.

Figure 8: GM-TO of an MBB beam.

Next, we consider an L-bracket (800 elements) illustrated in fig. 9. Using
only a single X-type microstructure, [54] achieved a compliance of 1130, 283,
171 for volume fractions of 0.1, 0.3 and 0.5 respectively. By considering all
microstructures (fig. 10), we achieved compliance of 1048, 263 and 158
respectively.

Figure 9: Loading of an L-bracket Figure 10: GM-TO of an L-bracket at various
volume fractions.

Finally, we consider an edge-loaded beam with a mesh of $32\times 20$ elements
illustrated in fig. 11(a). Here, the design is constrained to have the same
microstructure along the X-axis (but can vary along the Y-Axis), with a
desired volume fraction $V_{f}^{*}=0.6$. The topology and compliance obtained
(fig. 11(c)) is validated against the one reported by [15] (fig. 11(b)).

(a) Edge-loaded beam.

(b) M-TO design ($J=1.45e5$) obtained by [15].

(c) GM-TO design ($J=1.41e5$) obtained using the proposed method.

Figure 11: Validation for an edge-loaded beam.

###  3.2 Convergence Study

We illustrate the typical convergence of the proposed algorithm for a mid-
cantilever beam (fig. 1(a)) for $V_{f}^{*}=0.5$. The compliance, volume
fraction and the evolving topology at various instances are illustrated in
fig. 12. We observe that the log-barrier formulation leads to a stable
convergence. Similar convergence behavior was observed for all other examples.
The computation took $53$ secs.

Figure 12: Convergence of compliance and volume fraction for a mid-cantilever
beam. Topologies are illustrated at the $0^{\text{th}}$, $50^{\text{th}}$,
$100^{\text{th}}$, $150^{\text{th}}$, $200^{\text{th}}$ and
$250^{\text{th}}$(final) iterations.

###  3.3 Varying number of microstructures

Two central hypotheses of the current work is that one can achieve better
designs with larger number of candidate microstructures, and the framework is
(computationally) insensitive to the number of candidates. To validate these,
consider the problem in fig. 8(a). The topologies and compliance values when
allowing for varying number of microstructures, in the order in which they are
present in fig. 6, are illustrated in fig. 13, As expected, the compliance
improves as we allow for larger number of microstructures. As noted earlier,
each additional candidate microstructure requires only one additional output
neuron (2 additional design variables), enabling us to consider large set of
microstructures if necessary. The computational time was approximately $56$
seconds, independent of the number of microstructures.

Figure 13: Design and compliance using 2, 4, 7 and 11 microstructures.

###  3.4 Pareto Designs

Exploring the Pareto-front is critical in making design choices and
understanding the trade-off between the objective and constraint. Consider
once again the mid-cantilever beam in fig. 1(a). We illustrate the compliance
values and topologies for varying volume fractions $V_{f}^{*}$ in fig. 14.

Figure 14: Pareto designs of mid-cantilever beam.

###  3.5 Mesh and NN Dependency

Next we study the effect of NN size and FE mesh on the computed topology using
the tip cantilever (fig. 2) as an example; all other parameters are kept
constant. For varying NN size, the topologies are reported in fig. 15; no
appreciable difference in computational time was observed. Similarly, fig. 16
captures the topologies for varying mesh size; the computational times were
7.5, 24.7 and 114.3 seconds per 100 iterations for the $40\times 20$,
$60\times 30$ and $80\times 40$ mesh respectively.

Figure 15: GM-TO of a tip cantilever at $V_{f}^{*}=0.45$ with varying NN size.
Figure 16: GM-TO of a tip cantilever at $V_{f}^{*}=0.25$ with varying mesh
size.

###  3.6 High Resolution Design

In the proposed method, one can perform optimization on a coarse mesh and then
re-sample at a higher resolution to populate the microstructures, resulting in
smoother variation in the microstructure topologies. Note that, even after re-
sampling, the framework guarantees partition of unity. However, due to NN-
interpolation, microstructure mixing can occur on the finer mesh at the
interfaces. This was resolved by choosing the microstructure $m$ with the
largest $\rho_{m}$. This is illustrated in fig. 17 where we optimize on a
$40\times 20$ mesh and then re-sample on a $80\times 40$ mesh.

Figure 17: Extraction of high resolution design with smoother variation of
microstructures.

##  4 Conclusions

A framework for graded multi-scale topology optimization using neural network
was proposed and demonstrated. The salient features of this framework are: (1)
the number of design variables are only weakly dependent on the number of pre-
selected microstructures, (2) the partition of unity constraint is
automatically satisfied, and (3) manual sensitivity calculations is avoided.

The framework was limited to 2D compliance minimization problems, involving
microstructures governed by a single size parameter. Extension to 3D, non-
compliance problems (such as energy absorption [62], orthopedic implants [63],
resonant frequencies [64]) using more generic multi-parameter microstructures
need to be explored. Furthermore, while we relied on polynomials for directly
interpolating the elasticity components, it is desirable to consider an Eigen-
value decomposition ([54]) for increased robustness. The framework can
complement and might benefit from data driven approaches ([44, 45]). Finally,
while we relied on the simple Adam optimization procedure, second order
optimization methods such as LBFGS might result in better/faster convergence
([65]).

## Acknowledgements

The authors would like to thank the support of National Science Foundation
through grant CMMI 1561899, and the U. S. Office of Naval Research under
PANTHER award number N00014-21-1-2916 through Dr. Timothy Bentley.

## Replication of Results

The Python code pertinent to this paper is available at github.com/UW-
ERSL/GMTOuNN

## Compliance with ethical standards

The authors declare that they have no conflict of interest.

## References

  * [1] M. P. Bendsøe and O. Sigmund.  Topology optimization.  Springer Berlin Heidelberg, 2004.
  * [2] O. Sigmund and K. Maute.  Topology optimization approaches: A comparative review.  Structural and Multidisciplinary Optimization, 48(6):1031–1055, 2013.
  * [3] M. P. Bendsøe.  Optimal shape design as a material distribution problem.  Structural Optimization, 1(4):193–202, 1989.
  * [4] M. P. Bendsøe and N. Kikuchi.  Generating optimal topologies in structural design using a homogenization method.  Computer Methods in Applied Mechanics and Engineering, 71(2):197–224, nov 1988.
  * [5] B. Hassani and E. Hinton.  A review of homogenization and topology optimization I - Homogenization theory for media with periodic structure.  Computers and Structures, 69(6):707–717, 1998.
  * [6] J. A. Sethian and Andreas Wiegmann.  Structural Boundary Design via Level Set and Immersed Interface Methods.  Journal of Computational Physics, 163(2):489–528, 2000.
  * [7] M.Y. Wang, X. Wang, and D. Guo.  A level set method for structural topology optimization.  Computer Methods in Applied Mechanics and Engineering, 192(1-2):227–246, jan 2003.
  * [8] Y. M. Xie and G. P. Steven.  A simple evolutionary procedure for structural optimization.  Computers and Structures, 49(5):885–896, dec 1993.
  * [9] X. Y. Yang, Y. M. Xei, G. P. Steven, and O. M. Querin.  Bidirectional evolutionary method for stiffness optimization.  AIAA Journal, 37(11):1483–1488, nov 1999.
  * [10] W. Zhang, J. Yuan, J. Zhang, and X. Guo.  A new topology optimization approach based on Moving Morphable Components (MMC) and the ersatz material model.  Structural and Multidisciplinary Optimization, 53(6):1243–1260, 2016.
  * [11] K. Suresh.  A 199-line Matlab code for Pareto-optimal tracing in topology optimization.  Structural and Multidisciplinary Optimization, 42(5):665–679, 2010\.
  * [12] Shiguang Deng and Krishnan Suresh.  Multi-constrained topology optimization via the topological sensitivity.  Structural and Multidisciplinary Optimization, 51(5):987–1001, 2015\.
  * [13] H. C. Rodrigues, J. M. Guedes, and M. P. Bendsøe.  Hierarchical optimization of material and structure.  Structural and Multidisciplinary Optimization, 24(1):1–10, 2002\.
  * [14] Y. Wang, H. Xu, and D. Pasini.  Multiscale isogeometric topology optimization for lattice materials.  Computer Methods in Applied Mechanics and Engineering, 316:568–585, 2017.
  * [15] Tej Kumar and Krishnan Suresh.  A density-and-strain-based k-clustering approach to microstructural topology optimization.  Structural and Multidisciplinary Optimization, 61(4):1399–1415, 2020.
  * [16] Wenjin Tao and Ming C. Leu.  Design of lattice structure for additive manufacturing.  In 2016 International Symposium on Flexible Automation (ISFA), pages 325–332. IEEE, aug 2016.
  * [17] Y. Wang, L. Zhang, S. Daynes, H. Zhang, S. Feih, and M. Y. Wang.  Design of graded lattice structure with optimized mesostructures for additive manufacturing.  Materials and Design, 142:114–123, 2018.
  * [18] Lin Cheng, Xuan Liang, Jiaxi Bai, Qian Chen, John Lemon, and Albert To.  On utilizing topology optimization to design support structure to prevent residual stress induced build failure in laser powder bed metal additive manufacturing.  Additive Manufacturing, 27(March):290–304, 2019.
  * [19] Dawei Li, Ning Dai, Yunlong Tang, Guoying Dong, and Yaoyao Fiona Zhao.  Design and Optimization of Graded Cellular Structures with Triply Periodic Level Surface-Based Topological Shapes.  Journal of Mechanical Design, Transactions of the ASME, 141(7), 2019\.
  * [20] Fei Liu, Zhongfa Mao, Peng Zhang, David Z. Zhang, Junjie Jiang, and Zhibo Ma.  Functionally graded porous scaffolds in multiple patterns: New design method, physical and mechanical properties.  Materials & Design, 2018.
  * [21] R. Sivapuram, P. D. Dunning, and H. A. Kim.  Simultaneous material and structural optimization by multiscale topology optimization.  Structural and Multidisciplinary Optimization, 54(5):1267–1281, 2016.
  * [22] Jiao Jia, Daicong Da, Cha-Liang Loh, Haibin Zhao, Sha Yin, and Jun Xu.  Multiscale topology optimization for non-uniform microstructures with hybrid cellular automata.  Structural and Multidisciplinary Optimization, 62(2):757–770, 2020\.
  * [23] Yan Zhang, Hao Li, Mi Xiao, Liang Gao, Sheng Chu, and Jinhao Zhang.  Concurrent topology optimization for cellular structures with nonuniform microstructures based on the kriging metamodel.  Structural and Multidisciplinary Optimization, 59(4):1273–1299, apr 2019.
  * [24] Y. Zhang, M. Xiao, H. Li, L. Gao, and S. Chu.  Multiscale concurrent topology optimization for cellular structures with multiple microstructures based on ordered SIMP interpolation.  Computational Materials Science, 155:74–91, dec 2018.
  * [25] Tongyu Wu and Shu Li.  An efficient multiscale optimization method for conformal lattice materials.  Structural and Multidisciplinary Optimization, 63(3):1063–1083, 2021.
  * [26] Xuechen Gu, Shaoming He, Yihao Dong, and Tao Song.  An improved ordered simp approach for multiscale concurrent topology optimization with multiple microstructures.  Composite Structures, 287:115363, 2022.
  * [27] Dilaksan Thillaithevan, Paul Bruce, and Matthew Santer.  Stress-constrained optimization using graded lattice microstructures.  Structural and Multidisciplinary Optimization, 63(2):721–740, 2021\.
  * [28] C. Wang, J. H. Zhu, W. H. Zhang, S. Y. Li, and J. Kong.  Concurrent topology optimization design of structures and non-uniform parameterized lattice microstructures.  Structural and Multidisciplinary Optimization, 58(1):35–50, 2018\.
  * [29] Zijun Wu, Fei Fan, Renbin Xiao, and Lianqing Yu.  The substructuring‐based topology optimization for maximizing the first eigenvalue of hierarchical lattice structure.  International Journal for Numerical Methods in Engineering, page nme.6342, mar 2020.
  * [30] Zhen Liu, Liang Xia, Qi Xia, and Tielin Shi.  Data-driven design approach to hierarchical hybrid structures with multiple lattice configurations.  Structural and Multidisciplinary Optimization, 61(6):2227–2235, feb 2020.
  * [31] Liang Xu and Gengdong Cheng.  Two-scale concurrent topology optimization with multiple micro materials based on principal stress direction.  In World Congress of Structural and Multidisciplinary Optimisation, pages 1726–1737. Springer, 2017.
  * [32] Yaguang Wang and Zhan Kang.  Concurrent two-scale topological design of multiple unit cells and structure using combined velocity field level set and density model.  Computer methods in applied mechanics and engineering, 347:340–364, 2019.
  * [33] Han Zhou, Jihong Zhu, Chuang Wang, Yifei Zhang, Jiaqi Wang, and Weihong Zhang.  Hierarchical structure optimization with parameterized lattice and multiscale finite element method.  Structural and Multidisciplinary Optimization, 65(1):1–20, 2022\.
  * [34] J. Stegmann and E. Lund.  Discrete material optimization of general composite shell structures.  International Journal for Numerical Methods in Engineering, 62(14):2009–2027, apr 2005.
  * [35] Tong Gao and Weihong Zhang.  A mass constraint formulation for structural topology optimization with multiphase materials.  International Journal for Numerical Methods in Engineering, 88(8):774–796, nov 2011.
  * [36] Emily D. Sanders, Miguel A. Aguiló, and Glaucio H. Paulino.  Multi-material continuum topology optimization with arbitrary volume and mass constraints.  Computer Methods in Applied Mechanics and Engineering, 340:798–823, oct 2018.
  * [37] ED Sanders, A Pereira, and GH Paulino.  Optimal and continuous multilattice embedding.  Science Advances, 7(16):eabf4838, 2021.
  * [38] Fernando V Senhora, Emily D Sanders, and Glaucio H Paulino.  Optimally-tailored spinodal architected materials for multiscale design and manufacturing.  Advanced Materials, page 2109304, 2022.
  * [39] Cong Hong Phong Nguyen and Young Choi.  Multiscale design of functionally graded cellular structures for additive manufacturing using level-set descriptions.  Structural and Multidisciplinary Optimization, 64(4):1983–1995, 2021.
  * [40] Ruijie Zhao, Junpeng Zhao, and Chunjie Wang.  Stress-constrained multiscale topology optimization with connectable graded microstructures using the worst-case analysis.  International Journal for Numerical Methods in Engineering, 123(8):1882–1906, 2022.
  * [41] Chen Yu, Qifu Wang, Zhaohui Xia, Yingjun Wang, Chao Mei, and Yunhua Liu.  Multiscale topology optimization for graded cellular structures based on level set surface cutting.  Structural and Multidisciplinary Optimization, 65(1):1–17, 2022\.
  * [42] Xiao-Yi Zhou, Zongliang Du, and H Alicia Kim.  A level set shape metamorphosis with mechanical constraints for geometrically graded microstructures.  Structural and Multidisciplinary Optimization, 60(1):1–16, 2019\.
  * [43] Hui Liu, Hongming Zong, Ye Tian, Qingping Ma, and Michael Yu Wang.  A novel subdomain level set method for structural topology optimization and its application in graded cellular structure design.  Structural and Multidisciplinary Optimization, 60(6):2221–2247, 2019.
  * [44] Liwei Wang, Siyu Tao, Ping Zhu, and Wei Chen.  Data-driven topology optimization with multiclass microstructures using latent variable gaussian process.  Journal of Mechanical Design, 143(3), 2021.
  * [45] Liwei Wang, Anton van Beek, Daicong Da, Yu-Chin Chan, Ping Zhu, and Wei Chen.  Data-driven multiscale design of cellular composites with multiclass microstructures for natural frequency maximization.  Composite Structures, 280:114949, 2022.
  * [46] Li Zheng, Siddhant Kumar, and Dennis M Kochmann.  Data-driven topology optimization of spinodoid metamaterials with seamlessly tunable anisotropy.  Computer Methods in Applied Mechanics and Engineering, 383:113894, 2021.
  * [47] Seth Watts, William Arrighi, Jun Kudo, Daniel A. Tortorelli, and Daniel A. White.  Simple, accurate surrogate models of the elastic response of three-dimensional open truss micro-architectures with applications to multiscale topology design.  Structural and Multidisciplinary Optimization, nov 2019.
  * [48] Daniel A. White, William J. Arrighi, Jun Kudo, and Seth E. Watts.  Multiscale topology optimization using neural network surrogate models.  Computer Methods in Applied Mechanics and Engineering, 346:1118–1135, apr 2019.
  * [49] Darshil Patel, Dustin Bielecki, Rahul Rai, and Gary Dargush.  Improving connectivity and accelerating multiscale topology optimization using deep neural network techniques.  Structural and Multidisciplinary Optimization, 65(4):1–19, 2022\.
  * [50] Aaditya Chandrasekhar and Krishnan Suresh.  Tounn: topology optimization using neural networks.  Structural and Multidisciplinary Optimization, 63(3):1135–1149, 2021.
  * [51] Emily D. Sanders, Miguel A. Aguiló, and Glaucio H. Paulino.  Multi-material continuum topology optimization with arbitrary volume and mass constraints.  Computer Methods in Applied Mechanics and Engineering, 340:798–823, oct 2018.
  * [52] Aaditya Chandrasekhar and Krishnan Suresh.  Multi-material topology optimization using neural networks.  Computer-Aided Design, 136:103017, 2021.
  * [53] Prajit Ramachandran, Barret Zoph, and Quoc V Le.  Searching for activation functions.  arXiv preprint arXiv:1710.05941, 2017.
  * [54] Tej Kumar, Saketh Sridhara, Bhagyashree Prabhune, and Krishnan Suresh.  Spectral decomposition for graded multi-scale topology optimization.  Computer Methods in Applied Mechanics and Engineering, 377:113670, 2021.
  * [55] Erik Andreassen and Casper Schousboe Andreasen.  How to determine composite material properties using numerical homogenization.  Computational Materials Science, 83:488–495, 2014.
  * [56] Diederik P. Kingma and Jimmy Lei Ba.  Adam: A method for stochastic optimization.  In 3rd International Conference on Learning Representations, ICLR 2015 - Conference Track Proceedings. International Conference on Learning Representations, ICLR, dec 2015.
  * [57] Hoel Kervadec, Jose Dolz, Jing Yuan, Christian Desrosiers, Eric Granger, and I Ben Ayed.  Constrained deep networks: Lagrangian optimization via log-barrier extensions.  CoRR, abs/1904.04205, 2(3):4, 2019.
  * [58] Aaditya Chandrasekhar, Saketh Sridhara, and Krishnan Suresh.  Auto: a framework for automatic differentiation in topology optimization.  Structural and Multidisciplinary Optimization, 64(6):4355–4365, 2021.
  * [59] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang.  JAX: composable transformations of Python+NumPy programs, 2018\.
  * [60] O. Sigmund and J. Petersson.  Numerical instabilities in topology optimization: A survey on procedures dealing with checkerboards, mesh-dependencies and local minima.  Structural Optimization, 16(1):68–75, 1998.
  * [61] Razvan Pascanu, Tomas Mikolov, and Yoshua Bengio.  Understanding the exploding gradient problem.  CoRR, abs/1211.5063, 2012.
  * [62] Wen Zhang, TX Yu, and Jun Xu.  Uncover the underlying mechanisms of topology and structural hierarchy in energy absorption performances of bamboo-inspired tubular honeycomb.  Extreme Mechanics Letters, 52:101640, 2022.
  * [63] Nicola Ferro, Simona Perotto, Daniele Bianchi, Raffaele Ferrante, and Marco Mannisi.  Design of cellular materials for multiscale topology optimization: application to patient-specific orthopedic devices.  Structural and Multidisciplinary Optimization, 65(3):1–26, 2022\.
  * [64] Morgan Nightingale, Robert Hewson, and Matthew Santer.  Multiscale optimisation of resonant frequencies for lattice-based additive manufactured structures.  Structural and Multidisciplinary Optimization, 63(3):1187–1201, 2021.
  * [65] Jorge Nocedal and Stephen Wright.  Numerical optimization.  Springer Science & Business Media, 2006.