# Deep energy method in topology optimization applications

Junyan He1 Shashank Kushwaha1 Charul Chadha1 Seid Koric1,2 Diab Abueidda2
Iwona Jasiuk1 1 Department of Mechanical Science and Engineering, University
of Illinois at Urbana-Champaign, Champaign, IL, USA
2 National Center for Supercomputing Applications, University of Illinois at
Urbana-Champaign, Champaign, IL, USA<EMAIL_ADDRESS>
###### Abstract

This paper explores the possibilities of applying physics-informed neural
networks (PINNs) in topology optimization (TO) by introducing a fully self-
supervised TO framework that is based on PINNs. This framework solves the
forward elasticity problem by the deep energy method (DEM). Instead of
training a separate neural network to update the density distribution, we
leverage the fact that the compliance minimization problem is self-adjoint to
express the element sensitivity directly in terms of the displacement field
from the DEM model, and thus no additional neural network is needed for the
inverse problem. The method of moving asymptotes is used as the optimizer for
updating density distribution. The implementation of Neumann, Dirichlet, and
periodic boundary conditions are described in the context of the DEM model.
Three numerical examples are presented to demonstrate framework capabilities:
(1) Compliance minimization in 2D under different geometries and loading, (2)
Compliance minimization in 3D, and (3) Maximization of homogenized shear
modulus to design 2D meta material unit cells. The results show that the
optimized designs from the DEM-based framework are very comparable to those
generated by the finite element method, and shed light on a new way of
integrating PINN-based simulation methods into classical computational
mechanics problems.

## Graphical abstract

###### keywords:

Deep energy method , Physics-informed neural networks , Topology optimization
, SIMP

††journal: arXiv

*[inlinelist,1]label=(),

##  1 Introduction

Neural networks (NNs) have seen wide applications in computational mechanics.
They have been used to learn the complex deformation process of binary
composites [1, 2, 3], thin-walled shell structures [4, 5], plasticity [6, 7,
8] and to learn constitutive laws [9, 10, 11] from a large collection of
finite element simulations training set. These models seek to approximate the
underlying relations embedded in finite element solutions, but do not seek to
solve the underlying partial differential equations (PDEs) representing the
governing physics principles. In contrast, physics-informed neural networks
(PINNs) [12, 13, 14, 15, 13] present an alternative to the finite element
method (FEM) to numerically solve the governing PDEs directly.

PINNs construct a NN approximation to the PDE solution at discrete points, and
the loss function is defined based on how well the underlying physics is
satisfied at these training points. One such example of PINN is the deep
collocation method [16, 17, 18, 13], which is based directly on the strong
form of the governing PDE, and seeks to minimize a loss function defined as
the $L_{2}$ norm of the residual in the PDE at discrete points. Since it is
based on the strong form, naturally it involves the computation of second-
order (e.g., in elasticity, heat transfer and the Navier-Stokes equations) or
higher-order (e.g., plate bending) spatial gradients of the outputs through
automatic differentiation of the NN, which can be quite expensive. However,
they are relatively straightforward to implement and is applicable to all
PDEs, as it mainly involves numerically evaluating the differential operator
in the PDE. Therefore, strong-form based PINNs have been applied to solve a
wide variety of problems such as linear elasticity, hyperelasticity,
plasticity, composite mechanics [19, 20], micromechanics [15], and fluid
mechanics [14, 21].

An alternative method to avoid calculating high-order gradient is to write the
governing equation into a variational form, by finding a functional whose
variation yields the governing PDE; this is termed the deep Ritz method [22,
23]. For some physical systems, the governing physics dictate that the
solution of the PDE yields the minimum of the energy functional, such as the
principle of minimum potential energy in structural mechanics [24, 25]. In
these cases, it is natural to define the energy functional as the loss
function, as the training of the NN is a process where the loss value is
minimized. This energy-based method is termed the deep energy method (DEM)
[26, 27, 8]. Although DEM avoids working direcly with the strong form of the
PDE, it is only limited to physical systems where the solution is given by a
minimum potential, which is not applicable in cases like fluid mechanics.
Additionally, researchers have proposed and implemented a mixed formulation,
where both the strong form and the energy method are used to capture the
mechanical response with high solution gradients and stress concentrations
[28, 29, 30].

NNs are also widely used in topology optimization (TO) applications. NNs have
been trained on optimized designs generated by finite-element TO simulations
to rapidly predict optimized structures for new loading and boundary
conditions [31, 32, 33, 34], and are used as a way to parameterize the TO
design space as an alternative to the element-based parameterization in FEM-
based TO [35, 36, 37]. Recently, the work by Zehnder et al. [38] outlined a
new way of applying PINNs to TO, where the forward elasticity problem is
solved via a DEM-based technique, and another neural network is trained for
the inverse design problem to generate density-based gradient and update the
design density. This is a fully self-supervised approach, where no inputs from
finite-element simulations are needed. Although the training time is longer
compared to FEM solution time, the NN-based TO framework generated designs
that have similar, if not better performance than FEM-based TO.

Inspired by the work of Zehnder et al. [38], we introduce a simpler, yet fully
self-supervised approach for TO based on the DEM method to shed light on how
PINN-based simulation techniques can be combined with traditional solution
techniques in computational mechanics to provide novel solutions to classical
problems. This paper is organized as follows: Section 2 presents an overview
of deep neural networks, deep energy method, and topology optimization.
Section 3 presents and discusses the results from three numerical examples.
Section 4 summarizes the outcomes and limitation of this current study and
highlights possible future works.

##  2 Methods

###  2.1 Deep neural networks

Deep neural networks (DNNs) consist of multiple layers of interconnected
neurons, and the layers themselves are also connected. The configuration and
function of DNNs bear resemblance to that of a human brain. Fig. 1 shows a
typical structure of a fully connected neural network model.

Figure 1: Typical architecture of a fully connected neural network.

The architecture of a DNN can be characterized by the number of hidden layers
and the number of neurons in each hidden layer. The neurons of consecutive
layers are connected by a set of weights $\bm{W}$ and biases $\bm{b}$. The
output $\bm{y}^{i}$ of layer $i$ can be related to its corresponding weights
and biases as:

$\bm{y}^{i}=f_{\rm{act}}^{i}(\bm{W}^{i}\bm{y}^{i-1}+\bm{b}^{i}),$ (1)

where $f_{\rm{act}}^{i}$ denote the activation function for this layer.
Training of the DNN refers to the iterative process where the weights and
biases of the model are updated by means of gradient descent [39] such that
the value of a user-defined loss function is minimized.

###  2.2 The deep energy method

In this section, we describe DEM in the context of the elasticity equation in
structural mechanics. We consider a homogeneous, isotropic, linear elastic
body undergoing small deformation. We remark that DEM is not limited to linear
elastic materials and has been applied in nonlinear hyperelastic material
models [27]. In the absence of any body and inertial forces, the equilibrium
equation in tensorial form reads:

$\displaystyle\nabla\cdot\bm{\sigma}=\bm{0},\;\;\forall\bm{X}\in\Omega.$ (2)

The system is subject to Dirichlet boundary condition:

$\displaystyle\bm{u}=\bar{\bm{u}},\;\;\forall\bm{X}\in\partial\Omega_{u},$ (3)

and/or Neumann boundary condition:

$\displaystyle\bm{\sigma}\cdot\bm{n}=\bar{\bm{t}},\;\;\forall\bm{X}\in\partial\Omega_{t}.$
(4)

In the small deformation setting, the strain tensor is given by:

$\bm{\epsilon}=\frac{1}{2}(\nabla\bm{u}+\nabla\bm{u}^{T}).$ (5)

Stress can be computed from the constitutive law as:

$\bm{\sigma}=\frac{E}{1+\nu}\bm{\epsilon}+\frac{E\nu}{(1+\nu)(1-2\nu)}tr(\bm{\epsilon})\bm{I},$
(6)

where $E$ and $\nu$ are the Young’s modulus and Poisson’s ratio. DEM seeks the
solution to the equilibrium equation via the principle of minimum potential
energy. For a body in static equilibrium with no applied body forces, the
potential energy of the system reads:

$\psi(\bm{u})=\frac{1}{2}\int_{\Omega}\bm{\sigma}:\bm{\epsilon}\,dV-\int_{\partial\Omega_{t}}\bar{\bm{t}}\cdot\bm{u}\,dA.$
(7)

The loss function ($\mathcal{L}$) in DEM is defined identically as the
potential of the system:

$\mathcal{L}(\bm{u})=\psi(\bm{u}).$ (8)

The training of the DEM model
$\mathcal{M}:\mathbb{R}^{n}\to\mathbb{R}^{n}/\bm{X}\Rightarrow\tilde{\bm{u}}(\bm{X})$
111Note that the immediate output of DEM model is $\tilde{\bm{u}}$, which is
different from $\bm{u}$ in Eq. (5) and Eq. (7).  can be viewed as minimizing
the loss function (system potential energy) $\mathcal{L}$ through an
optimizer. The solution to the elasticity problem, as given by the DEM model,
is defined as:

$\bm{u}^{*}={\rm{arg}}\min_{\bm{u}}\,\mathcal{L}(\bm{u}).$ (9)

In this work, the L-BFGS optimizer [40] was used to train the model. The
learning rate is another parameter that can be changed by the user to achieve
optimal framework performance. Similar to the early stopping schedule commonly
used in model training, we monitor the relative change of the loss function
and stop the training once the relative change in loss function value is less
than the user-specified tolerance $\epsilon_{tol}$.

The underlying architecture of the DEM model is a fully connected feed-forward
network, which is similar to the one depicted in Fig. 1 and was described in
the work of Nguyen et al. [27]. Fourier transform was applied using the
package Random Fourier Features Pytorch (RFF) [41] to the input features to
transform them to frequency domain, which was shown to increase the accuracy
of PINNs in the work of Wang et al. [42]. The NN architecture is fully
parametrizable, and is characterized by the number of layers, number of
neurons (per layer), the activation function, and two standard deviation
values used to initialize the weights and biases of the MLP and RFF.
Hyperparameter optimization can be performed to improve the performance of the
DEM model.

###  2.3 Enforcement of boundary conditions

In the context of elasticity problems, three different types of boundary
conditions (BCs) are common, namely the Neumann, Dirichlet, and periodic BCs.
Enforcement of Neumann BCs is straightforward. The value of the prescribed
traction $\bar{\bm{t}}$ enters directly in the boundary integral part of Eq.
(7), while the boundary displacements are provided by the DEM model. Traction-
free surfaces are satisfied trivially as the boundary integral vanishes on
these surfaces. A surface (in 3D) or edge (in 2D) integration rule needs to be
defined for numerically evaluating the surface integral in regions where
traction is prescribed.

Dirichlet BCs are typically enforced via the penalty method by adding a
penalty term to the loss function [27, 43, 17]:

$MSE_{\partial\Omega_{u}}=\frac{1}{N_{u}}\sum_{j=1}^{N_{u}}||\bm{u}-\bar{\bm{u}}||^{2},$
(10)

leading to a modified loss function definition:

$\mathcal{L}(\bm{u})=\psi(\bm{u})+w\,MSE_{\partial\Omega_{u}},$ (11)

where $w$ is a user-defined weight parameter. To avoid defining the additional
parameter $w$, we adopted a direct approach to enforce Dirichlet BCs via an
additive decomposition:

$\bm{u}(\bm{X})=\tilde{\bm{u}}(\bm{X})*\bm{m}(\bm{X})+\bm{u}_{0}(\bm{X}),$
(12)

where $*$ denotes element-wise multiplication between two vectors, the vectors
$\bm{m}$ and $\bm{u}_{0}$ are chosen such that:

$\bm{m}(\bm{X})=\bm{0}\,\,{\rm{and}}\,\,\bm{u}_{0}(\bm{X})=\bar{\bm{u}},\forall\bm{X}\in\partial\Omega_{u}$
(13)

By design, $\bm{u}(\bm{X})$ satisfies all Dirichlet boundary conditions for
all arbitrary $\tilde{\bm{u}}$ produced by the DEM model.

Periodic BCs are commonly used when simulating a representative volume element
that repeats infinitely in space. It is also used when simulating unit cells
of periodic meta materials. In the context of the FEM, periodic BCs are
typically stated in the form of constraint equations [44]:

$u_{i}^{k+}-u_{i}^{k-}=\epsilon^{0}_{ij}(x_{j}^{k+}-x_{j}^{k-}),$ (14)

where $k-$ and $k+$ denote the left and right instances of the $k^{th}$
periodic boundary pair, and $\epsilon^{0}_{ij}$ is the (known) applied strain
tensor. For easier convergence of the DEM training process, we adopted a mixed
implementation for periodic BCs, where we treat zero and non-zero applied
strain components differently. The boundary displacement induced by a non-zero
applied strain component is enforced as a pair of Dirichlet BCs, and zero
applied strain components are enforced by penalty. Without loss of generality,
let the non-zero component of $\epsilon^{0}_{ij}$ be $\epsilon_{a}$, which is
applied on the $k^{th}$ periodic boundary pair. This leads to two Dirichlet
boundary conditions:

$\displaystyle\bm{u}(\bm{X})=\bm{0},\;\;\forall\bm{X}\in\partial\Omega_{k-},$
(15) $\displaystyle\bm{u}(\bm{X})_{a}=\Delta
L_{k}\epsilon_{a},\;\;\forall\bm{X}\in\partial\Omega_{k+},$

where $\Delta L_{k}$ is the distance between the $k+$ and $k-$ boundary pair.
For all other boundary pairs that correspond to zero applied strain
components, the constraint equation Eq. (14) reduces to:

$u_{i}^{k+}=u_{i}^{k-}.$ (16)

Therefore, the following penalty term is added to the loss function:

$MSE_{PBC}=\frac{E}{N_{p}}\sum_{j=1}^{N_{p}}||u_{i}^{k+}-u_{i}^{k-}||^{2},$
(17)

where $E$ is the Young’s modulus of the material (a physics-informed penalty
weight) and $N_{p}$ is the number of points of the $k^{th}$ periodic boundary.

###  2.4 Topology optimization

In this section, we describe topology optimization (TO) in the context of DEM.
The mathematical formulation of the TO problem using DEM can be stated as:

$\displaystyle\min_{\bm{\rho}}f(\bm{u},\bm{\rho})$ (18)
$\displaystyle{\rm{s.t.}}\;\bm{u}={\rm{arg}}\min_{\bm{u}}\,\mathcal{L}(\bm{u},\bm{\rho}),$
$\displaystyle V_{e}\sum\bm{\rho}=\bar{V},$ $\displaystyle 0\leq\rho_{e}\leq
1,\,\forall e=1\cdots N,$

where $f(\bm{u},\bm{\rho})$, $\bm{\rho}$, $V_{e}$ and $\bar{V}$ denote the
objective function to minimize, element density vector, element volume and
target volume, respectively. A typical objective of TO is to minimize the
system compliance subject to a volume/mass constraint. For compliance
minimization, the objective function is:

$f_{\rm{comp}}(\bm{u},\bm{\rho})=\frac{1}{2}\int_{\Omega}\bm{\sigma}^{*}(\bm{\rho}):\bm{\epsilon}\,dV$
(19)

where $\bm{\sigma}^{*}(\bm{\rho})$ is the stress computed from the classical
solid isotropic material penalization (SIMP) scheme [45]:

$\bm{\sigma}^{*}(\bm{\rho})=\bm{\rho}^{p}*\bm{\sigma}.$ (20)

The penalty exponent $p$ is taken to be 3 in this work.

For designing periodic 2D meta structures, a commonly employed objective is to
maximize the homogenized shear modulus of the unit cell [31, 44, 46]:

$G_{\rm{homo}}=\frac{1}{A}\int_{\Omega}\bm{\sigma}^{p*}(\bm{\rho}):\bm{\epsilon}^{p}\,dA,$
(21)

where $A$ is the area of the 2D element, $\bm{\sigma}^{p}$ and
$\bm{\epsilon}^{p}$ denote the stress and strain fields computed from the
_periodic_ displacement field subjected to the following applied simple shear
loading ($\bm{\epsilon}^{0}$ in Eq. (14) ):

$\bm{\epsilon}^{0}=\begin{bmatrix}0&0.01\\\ 0.01&0\end{bmatrix}.$ (22)

To fit into a minimization framework, the objective function is defined as:

$f_{\rm{shear}}(\bm{u},\bm{\rho})=-G_{\rm{homo}}.$ (23)

Due to the introduction of the element density vector, we correspondingly
modify the definition of the potential energy and loss function as:

$\psi(\bm{u},\bm{\rho})=\mathcal{L}(\bm{u},\bm{\rho})=\frac{1}{2}\int_{\Omega}\bm{\sigma}^{*}(\bm{\rho}):\bm{\epsilon}\,dV-\int_{\partial\Omega_{t}}\bar{\bm{t}}\cdot\bm{u}\,dA.$
(24)

In this work, we chose to evaluate the volume integral in Eq. (24) by Gauss
quadrature integration. This necessitates the formation of quadrilateral (in
2D) or hexahedral (in 3D) isoparametric finite elements from the structured
nodes $\bm{X}$ in the DEM domain, which is very similar to a FEM-type
treatment. The need to form elements is a major disadvantage of the framework,
as compared to other meshless DEM methods like those developed in [17, 27]
which only require nodes but not elements to exist in the simulation domain.
The gradient operator in Eq. (5) is evaluated by the gradients of the finite
element shape functions, and not from automatic differentiation of the DEM
model. The density vector $\bm{\rho}$ is defined at the cell center, identical
in classical FEM-based TO. The method of moving asymptotes (MMA) [47] is used
to iteratively update the density based on the sensitivity.

###  2.5 Sensitivity analysis

In this section, we describe the sensitivity analysis in the context of DEM.
Comparing Eq. (19) and Eq. (23), we see that both are similar in form to a
classical compliance minimization problem in FEM-based TO. Since the
compliance minimization problem is self-adjoint, the adjoint vector is simply
the negative of the displacement, and therefore, no additional adjoint
analysis is needed to evaluate the sensitivity. In this case, one simply needs
to express the sensitivity in terms of the displacement field and its
gradients produced by DEM. For compliance minimization problem (Eq. (19)), the
sensitivity of element $e$ is given by:

$\frac{\partial
f_{\rm{comp}}}{\partial\rho_{e}}=-\frac{1}{2}p\rho_{e}^{p-1}\int_{\Omega_{e}}\bm{\sigma}:\bm{\epsilon}\,dV.$
(25)

Similarly, for maximization of shear modulus in 2D meta material unit cells
(Eq. (23)), the sensitivity of element $e$ is given by:

$\frac{\partial
f_{\rm{shear}}}{\partial\rho_{e}}=\frac{1}{A_{e}}p\rho_{e}^{p-1}\int_{\Omega_{e}}\bm{\sigma}:\bm{\epsilon}\,dA.$
(26)

To avoid checkerboarding and mesh-dependence, the density filtering technique
introduced by Bruns et al. [48] was employed. In this approach, the MMA solver
does not directly manipulate $\bm{\rho}$, but instead manipulates on a pseudo-
density (the design variables) $\bm{\xi}$ that is related to $\bm{\rho}$ as:

$\bm{\rho}=\bm{W}\bm{\xi},$ (27)

where $\bm{W}$ is the normalized filter matrix whose rows are given by the
normalized weights $\bar{q}_{ij}$:

$\displaystyle q_{ij}={\rm{max}}(r_{min}-||\bm{X}_{i}-\bm{X}_{j}||),$ (28)
$\displaystyle\bar{q}_{ij}=\frac{1}{\sum_{k=1}^{N}q_{ik}}q_{ij}.$

Algorithm 1 summarizes the process of TO using the DEM method.

Input: Network architecture, domain size, grid size, material properties,
$\epsilon_{tol}$ , $r_{min}$, $\bar{V}$, initial density distribution
$\bm{\xi}^{0}$, maximum TO iterations

1

Output: Optimized design defined by density array $\bm{\rho}$

2

/* Initialization */

3 Build filter matrix $\bm{W}$

4 Initialize weights and biases of the DEM model $\mathcal{M}$

5 $i\leftarrow 0$

/* Begin topology optimization */

6 while  _ $i<$ max TO iteration _ do

7 Apply filter: $\bm{\rho}^{i}=\bm{W}\bm{\xi}^{i}$

/* Begin DEM training */

8 while  _not converged_ do

9 Obtain $\tilde{\bm{u}}$ from $\mathcal{M}$, apply Dirichlet BCs to get
$\bm{u}$

10 Use shape function gradient operator to compute $\bm{\epsilon}$

11 Use $\bm{\rho}$ to compute $\bm{\sigma}^{*}$

12 Compute loss

13 Update the weights and biases of $\mathcal{M}$ through back-propagation

14

/* Begin sensitivity calculation */

15 Obtain converged $\tilde{\bm{u}}$ from $\mathcal{M}$, apply Dirichlet BCs

16 Compute $\bm{\epsilon}$ and $\bm{\sigma}^{*}$

17 Compute the element-wise filtered sensitivity

18 Invert filtering for sensitivity: $\frac{\partial
f}{\partial\bm{\xi}}=\bm{W}^{T}\frac{\partial f}{\partial\bm{\rho}}$

19 Pass sensitivity to MMA solver to obtain $\bm{\xi}^{i+1}$

20 $i\leftarrow i+1$

21

Algorithm 1 Topology optimization using the deep energy method

##  3 Results and discussion

In this section, we present three numerical examples that showcase the
capabilities of our DEM-based TO framework. In the first example, we performed
compliance minimization in 2D for two different geometries: 1. a long slender
bridge fixed on both ends subject to center downward load, and 2. a short
cantilever beam subject to downward load on its right free end. In the second
example, we extend the bridge example to 3D. In the last example, we maximize
the shear modulus of 2D periodic unit cells to design meta materials, and
compare the results at different initial conditions. In all cases, we compared
our DEM results with those obtained from FEM. All training of the DEM model
was done on Delta, an HPC cluster hosted at the National Center for
Supercomputing Applications (NCSA). A representative GPU compute node on Delta
has the following hardware specifications: one AMD Milan CPU, four NVIDIA A40
GPUs, and 256 GB RAM. The hyperparameters of the DEM model are provided in
Table 1, and were obtained following the procedure described in the work by
Chadha et al. [49].

Table 1: Hyper parameters of the DEM model | Layers | Neurons | Activation function | Max iteration | learning rate | $\sigma_{MLP}$ | $\sigma_{RFF}$
---|---|---|---|---|---|---|---
Value | 5 | 68 | RReLU | 100 | 1.735 | 0.0622 | 0.1192

The DEM model was implemented in PyTorch (version 1.11.0) [50]. Material
properties used in all three examples are: $E=200$ MPa, and $\nu=0.3$. Plane
stress condition and unit thickness were assumed for all 2D simulations. A
target volume fraction of 40% was used. A Python implementation of the MMA
solver by Chandrasekhar et al. [51] was used. For all DEM examples presented,
a total of 80 TO iterations were conducted using a filter radius of 0.25 m.

###  3.1 Compliance minimization in 2D

Two different cases are presented in this example. In the first case, a 2D
rectangular bridge of size 12-by-2 m2 was subjected to a downward surface
traction at the center of the top edge over a length of 0.5 m. The bridge was
fixed at its left and right edges and 3751 nodes were placed into the domain
forming a 121-by-31 grid. In the second case, a 2D rectangular beam of size
10-by-5 m2 was subjected to a downward concentrated load at the center of the
right edge, while the left edge was fixed. 4186 nodes were placed into the
domain, forming a 91-by-46 grid. The relative convergence tolerance
$\epsilon_{tol}$ was set to $5\times 10^{-6}$ for all DEM training. As
comparison, FEM-based TO was applied to solve the same problem using identical
node layout.

The designs from DEM- and FEM-based TO at different design iterations are
compared in Figures 2 and 3. The relative compliance reduction for both cases
as the TO optimizer progressed is presented in Fig. 4a. To give a quantitative
comparison of the similarity between DEM- and FEM-based designs, the dice
similarity coefficient (DSC) was computed for both cases and is presented in
Fig. 4b. The DSC is defined as:

${\rm{DSC}}=\frac{2|\bm{I}_{FEM}*\bm{I}_{DEM}|}{|\bm{I}_{FEM}|+|\bm{I}_{DEM}|},$
(29)

where $\bm{I}$ represents the binarized density using a (dynamic) threshold of
0.5 times the current maximum density. Finally, the computational cost for
running the two cases using DEM and FEM222The highly optimized Python TO code
found in [52] was used for time comparison, which was based on the MATLAB code
by Andreassen et al. [53] are summarized in Fig. 4c and are compared in Table
2.

(a) DEM, iteration 2 |  (b) DEM, iteration 15 |  (c) DEM, iteration 30 |  (d) DEM, iteration 80
---|---|---|---
(e) FEM, iteration 2 |  (f) FEM, iteration 15 |  (g) FEM, iteration 30 |  (h) FEM, iteration 80
Figure 2: Density plot of designs generated by DEM and FEM for the bridge example. Both left and right edges are fixed, and a downward load is applied on the center of the top edge. (a) DEM, iteration 2 |  (b) DEM, iteration 15 |  (c) DEM, iteration 30 |  (d) DEM, iteration 80
---|---|---|---
(e) FEM, iteration 2 |  (f) FEM, iteration 15 |  (g) FEM, iteration 30 |  (h) FEM, iteration 80
Figure 3: Density plot of designs generated by DEM and FEM for the beam
example. The left edge is fixed and a downward load is applied at the center
of the right edge.

(a)

(b)

(c)

Figure 4: Comparing DEM-based TO with FEM-based TO in 2D: LABEL:sub@fig:compliance_history Relative compliance evolution for the two cases. LABEL:sub@fig:DSC Dice similarity coefficient computed on the binarized designs. LABEL:sub@fig:train2d Training time for each iteration. Table 2: Computational time for the two cases in Example 1 | FEM-based | DEM-based
---|---|---
Bridge | 4.1s | 460.7s
Beam | 12.7s | 485.7s

From Fig. 2 and Fig. 3, we clearly see that the DEM- and FEM-based TO
frameworks produced almost identical simulation design evolution history and
final designs. The similarity in designs is further demonstrated in Fig. 4b,
where we see that the DSC is greater than 90% for the final designs in both
cases. For the bridge case, we noticed that the DSC remains greater than 98%
for all design iterations, showing close resemblance between the DEM- and FEM-
based approaches. The two frameworks also produce designs of similar
performance, as is evident in Fig. 4a, where we see that the compliance
reduction history is almost identical between the two methods due to the use
of identical MMA optimizer. Nonetheless, Table 2 shows that FEM still holds a
massive advantage over DEM in terms of computational time, which agrees with
the trends observed in the similar DEM-based TO work like [38]. However, a
very tight relative tolerance of $5\times 10^{-6}$ was used in training, which
could be loosened to gain computational efficiency. From Fig. 4c, we highlight
two important features of DEM-based TO. For both cases, we noticed that the
training time shows a decreasing trend as design progressed. This is due to
the fact that the weights and biases from the previous TO iteration are used
as the initial condition for the next iteration, a form of transfer learning.
This leads to gradual training time reduction. This is in sharp contrast of
FEM-based TO, where the computational time for each iteration is constant, and
can be very beneficial when many TO iterations are needed. Also, we note that
the DEM training time went up by merely 5.4% from case 1 to case 2, where the
FEM-based TO time went up over 200%. The insensitivity of training time to
domain size is due to the fact that the number of parameters of DEM-based TO
is controlled by the DNN architecture, _not_ by the number of nodes used to
discretize the domain, which is beneficial for problems that need to be
resolved by a fine mesh.

###  3.2 Compliance minimization in 3D

The second example concerns a 3D rectangular bridge of size 12-by-2-by-2 m3.
It was fixed at both end surfaces and was subjected to a downward load at the
center of the top surface in a square region of size 0.5-by-0.5 m2. 75625
nodes were placed into the domain, forming a 121-by-25-by-25 grid.
$\epsilon_{tol}$ was set to $5\times 10^{-5}$. As comparison, Abaqus/Standard
[54] was used to perform TO using the built-in TO functionalities and
optimizer.

(a) DEM, iteration 80 |  (b) FEM, iteration 80
---|---
(c) Overlapping two designs |  (d) Section view
Figure 5: Comparison of final designs after 80 design iterations. The FEM-
based design is shown in the form of a stl surface, which the DEM design is
shown as voxels.

The final designs from DEM and FEM are presented in Fig. 5a and Fig. 5b,
respectively. The STL file for the final FEM design was extracted using a
density threshold of 0.8, and the same threshold was used to obtain the
voxelated DEM final design. The two designs are overlaid and presented in Fig.
5c and Fig. 5d. The relative compliance reduction and training time was shown
in Fig. 6. The total TO time for DEM was 1050.4s.

(a)

(b)

Figure 6: DEM-based TO in 3D: LABEL:sub@fig:compliance_history Relative
compliance evolution for the 3D bridge design. LABEL:sub@fig:train3d Training
time for each iteration.

From the overlaid plot and the section view in Fig. 5, we see that the DEM-
based design largely resembles the FEM-based design, especially near both ends
of the bridge. Some design differences can be seen in the center, especially
from the section view. From Fig. 6a, we see that DEM was able to quickly
reduce the compliance in less than 10 iterations, but has a final compliance
that is larger than that of the FEM design. The evolution histories are
different, which roots from the fact that different optimizers were used to
optimize the material distribution. From Fig. 6b, we again see that training
time decreased as design progressed, a trend that is consistent with Section
3.1. Finally, we highlight that the same hyperparameters in Section 3.1 were
used directly in this 3D example without any additional hyperparameter
optimization, which shows the robustness of the DEM model with respect to
hyperparameter selection.

###  3.3 Shear modulus maximization in 2D

In this example, a 2D unit cell of size 10-by-10 m2 was subjected to a shear
load at the top edge. The unit cell was subjected to periodic BCs in the X-
and Y-direction. Further, 6561 nodes were placed into the domain forming a
81-by-81 grid. The objective is to maximize the homogenized shear modulus
while maintaining the volume fraction at 40%. It is well-known that the final
optimized design depends heavily on the initial material distribution [31, 46,
44]. Therefore, three different initial designs were considered with different
initial hole diameters, as shown in the first row of Fig. 7. $\epsilon_{tol}$
was set to $5\times 10^{-6}$. For comparison, the DEM-based TO results are
compared to their FEM-based counterparts produced by the MATLAB [55] code
developed by Xia et al. [44].

The initial and final designs from DEM- and FEM-based TO are compared in Fig.
7. The relative change in shear modulus for both cases over all the iterations
are presented in Fig. 8a. Further, DSC coefficient is calculated for the
binarized final designs obtained from both cases using Eq. (29) and is
presented in Fig. 8b. The DEM training time for the three cases were 579.7,
934.4 and 1125.7s, respectively. Training time is summarized in Fig. 8c.

(a) Initial, r=H/4 |  (b) Initial, r=H/10 |  (c) Initial, r=H/20
---|---|---
(d) DEM, iteration 80, r=H/4 |  (e) DEM, iteration 80, r=H/10 |  (f) DEM, iteration 80, r=H/20
(g) FEM, iteration 80, r=H/4 |  (h) FEM, iteration 80, r=H/10 |  (i) FEM, iteration 80, r=H/20
Figure 7: Density plot of designs generated by DEM and FEM for 2D shear
modulus maximization. The first row shows three different initial
configurations, and the second and third rows show the corresponding final
designs.

(a)

(b)

(c)

Figure 8: Comparing DEM-based TO with FEM-based TO in 2D shear modulus
maximization: LABEL:sub@fig:compliance_meta Relative shear modulus evolution
for the two cases. LABEL:sub@fig:DSC_meta Dice similarity coefficient computed
on the binarized designs. LABEL:sub@fig:train_meta Training time for each
iteration.

From Fig. 7, we clearly see that the DEM- and FEM-based TO frameworks produced
very similar final designs for all three initial designs. However, comparing
Fig. 7f and Fig. 7i, we noticed an important deficiency of the DEM-based TO
method. The final FEM-based TO design (Fig. 7i) has four-fold symmetry, while
the DEM-based design (Fig. 7f) does not. In fact, the DEM-based final design
is not periodic at the lower left corner. This difference is due to that fact
that periodic BCs are enforced in FEM-based simulations by constraining the
corresponding degrees of freedom of the global stiffness matrix, hence
periodicity of the displacement field (hence the material distribution) is
also guaranteed. Also, if the initial design is symmetric, symmetry persists
in subsequent designs. This is not the case for DEM-based TO, as the periodic
BCs were enforced by penalty. The periodic penalty loss gradually creep up as
the design progressed, and the designs gradually lost periodicity due to minor
violations of periodicity in the displacement field and material distribution.
Nonetheless, it is observed that the DSC is greater than 90% for the final
designs in all three cases. The low DSC during the initial stage of evolution
is a result of the different optimizers used in the DEM- and FEM-based TO
code. Similar to the two examples presented above, we again see gradual
reduction of training time as TO iterations progressed.

##  4 Conclusions and future work

In this work, we present the use of DEM-based simulation technique in TO. The
proposed framework combines solution to PDEs using PINNs and classical TO,
where the momentum balance equation is solved by a deep neural network using
DEM by minimizing the system potential energy. Instead of training a separate
neural network to predict the updated density array as in the work of Zehnder
et al. [38], we leverage the fact that compliance minimization problem is
self-adjoint, and express element sensitivity directly using the displacement
field produced by the DEM model. Therefore, only one neural network needs to
be trained in each optimization iteration, making the framework much simpler.
The training of the DEM model also leverages transfer learning to reduce
training time, where the weights and biases from the previous iteration are
used as initial condition for the next iteration, which leads to significant
reduction of training time in later iterations. The proposed framework
combines two different types of optimizers, using a L-BFGS optimizer for the
DEM model training, and using a classical MMA optimizer for updating the
density array.

Three numerical examples are presented to demonstrate the capabilities of the
proposed framework. In the first example, we performed compliance minimization
in 2D on two different loading and boundary conditions, where we saw that the
DEM-based topology optimization framework generated almost identical designs
to its FEM counterpart. In the second example, we extended compliance
minimization to 3D, where we optimized a 3D bridge under transverse loading.
It was found that the hyperparameters used in 2D simulations were effective in
3D as well, generating a design that is quite similar to the one generated by
Abaqus. In the last example, we attempted the design of 2D meta materials by
maximizing the homogenized shear modulus, which necessitated the
implementation of periodic boundary conditions in DEM. For three different
initial designs, it was seen that our DEM framework was able to generate
designs very similar to FEM, albeit suffered from minor loss of symmetry due
to the inaccurate enforcement of periodic boundary conditions. In all examples
presented, the training time for DEM model was longer than the solution time
for FEM, but training time shows a gradual decrease with design iterations and
is relatively insensitive to the number of nodes in the simulation domain.
This makes the DEM-based TO method potentially beneficial when a fine mesh and
many TO iterations are required.

We conclude that the proposed DEM-based topology optimization framework can
produce optimized designs that are very similar to FEM. Although having a high
computational cost, we emphasize that the goal of this work is not to
outperform FEM in particular applications, but rather to highlight the diverse
application scenarios for the DEM method, and shed light on how neural-
network-based solutions for PDEs can be employed in engineering applications.

A primary limitation of this work is the need to form element discretization
of the domain, and using shape function gradients to compute gradients of
field variables produced by DEM, which is in contrast with many previous
studies [17, 38, 43] that used a meshless method to solve the underlying
governing equation. Extension to a meshless formulation using automatic
differentiation of the DEM model will be our future work. We will also explore
topology optimization for hyperelastic materials undergoing finite
deformation, where Newton-Raphson iterations need to be employed in FEM,
rendering it rather expensive. In addition, the effects of enhancing the loss
function with its gradient [56, 12] on the DEM model accuracy will be studied.

## Replication of results

The data and source code that support the findings of this study can be found
at: https://github.com/Jasiuk-Research-Group. Note to editor and reviewers:
the link above will be made public upon the publication of this manuscript.
During the review period, the data and source code can be made available upon
request to the corresponding author.

## Competing Interests

The authors declare that they have no conflict of interest.

## Acknowledgements

The authors would like to thank the National Center for Supercomputing
Applications (NCSA) for providing the computational resources on the Delta
cluster.

## CRediT author contributions

Junyan He: Conceptualization, Methodology, Software, Formal analysis,
Investigation, Data Curation, Writing - Original Draft. Shashank Kushwaha:
Software, Formal Analysis, Investigation, Writing - Original Draft. Charul
Chadha: Software, Formal Analysis, Investigation, Writing - Original Draft.
Seid Koric: Supervision, Resources, Writing - Review & Editing. Diab Abueidda:
Supervision, Resources, Writing - Review & Editing. Iwona Jasiuk: Supervision,
Resources, Writing - Review & Editing, Funding Acquisition.

## References

  * Chen and Gu [2019] Chun-Teh Chen and Grace X Gu.  Machine learning for composite materials.  _MRS Communications_ , 9(2):556–566, 2019.
  * Yang et al. [2020a] Charles Yang, Youngsoo Kim, Seunghwa Ryu, and Grace X Gu.  Prediction of composite microstructure stress-strain curves using convolutional neural networks.  _Materials & Design_, 189:108509, 2020a.
  * Luo et al. [2021] Ling Luo, Boming Zhang, Guowei Zhang, Xueqin Li, Xiaobin Fang, Weidong Li, and Zhenchong Zhang.  Rapid prediction and inverse design of distortion behaviors of composite materials using artificial neural networks.  _Polymers for Advanced Technologies_ , 32(3):1049–1060, 2021.
  * He et al. [2022] Junyan He, Shashank Kushwaha, Diab Abueidda, and Iwona Jasiuk.  Exploring the structure-property relations of thin-walled, 2d extruded lattices using neural networks, 2022.
  * Chen et al. [2019] Guorong Chen, Tiange Li, Qijun Chen, Shaofei Ren, Chao Wang, and Shaofan Li.  Application of deep learning neural network to identify collision load conditions based on permanent plastic deformation of shell structures.  _Computational Mechanics_ , 64(2):435–449, 2019\.
  * Stoffel et al. [2019] Marcus Stoffel, Franz Bamer, and Bernd Markert.  Neural network based constitutive modeling of nonlinear viscoplastic structural response.  _Mechanics Research Communications_ , 95:85–88, 2019.
  * Gorji et al. [2020] Maysam B Gorji, Mojtaba Mozaffar, Julian N Heidenreich, Jian Cao, and Dirk Mohr.  On the potential of recurrent neural networks for modeling path dependent plasticity.  _Journal of the Mechanics and Physics of Solids_ , 143:103972, 2020.
  * Abueidda et al. [2021a] Diab W Abueidda, Seid Koric, Nahil A Sobh, and Huseyin Sehitoglu.  Deep learning for plasticity and thermo-viscoplasticity.  _International Journal of Plasticity_ , 136:102852, 2021a.
  * Chen [2021] Guang Chen.  Recurrent neural networks (rnns) learn the constitutive law of viscoelasticity.  _Computational Mechanics_ , 67(3):1009–1019, 2021\.
  * Yang et al. [2020b] Hang Yang, Qian Xiang, Shan Tang, and Xu Guo.  Learning material law from displacement fields by artificial neural network.  _Theoretical and Applied Mechanics Letters_ , 10(3):202–206, 2020b.
  * Flaschel et al. [2021] Moritz Flaschel, Siddhant Kumar, and Laura De Lorenzis.  Unsupervised discovery of interpretable hyperelastic constitutive laws.  _Computer Methods in Applied Mechanics and Engineering_ , 381:113852, 2021.
  * Shukla et al. [2022] Khemraj Shukla, Mengjia Xu, Nathaniel Trask, and George E Karniadakis.  Scalable algorithms for physics-informed neural and graph networks.  _Data-Centric Engineering_ , 3, 2022.
  * Haghighat et al. [2021] Ehsan Haghighat, Maziar Raissi, Adrian Moure, Hector Gomez, and Ruben Juanes.  A physics-informed deep learning framework for inversion and surrogate modeling in solid mechanics.  _Computer Methods in Applied Mechanics and Engineering_ , 379:113741, 2021.
  * Cai et al. [2022] Shengze Cai, Zhiping Mao, Zhicheng Wang, Minglang Yin, and George Em Karniadakis.  Physics-informed neural networks (pinns) for fluid mechanics: A review.  _Acta Mechanica Sinica_ , pages 1–12, 2022.
  * Henkes et al. [2022] Alexander Henkes, Henning Wessels, and Rolf Mahnken.  Physics informed neural networks for continuum micromechanics.  _Computer Methods in Applied Mechanics and Engineering_ , 393:114790, 2022.
  * Raissi [2018] Maziar Raissi.  Deep hidden physics models: Deep learning of nonlinear partial differential equations.  _The Journal of Machine Learning Research_ , 19(1):932–955, 2018.
  * Abueidda et al. [2021b] Diab W Abueidda, Qiyue Lu, and Seid Koric.  Meshless physics-informed deep learning method for three-dimensional solid mechanics.  _International Journal for Numerical Methods in Engineering_ , 122(23):7182–7201, 2021b.
  * Guo et al. [2021] Hongwei Guo, Xiaoying Zhuang, and Timon Rabczuk.  A deep collocation method for the bending analysis of kirchhoff plate.  _arXiv preprint arXiv:2102.02617_ , 2021.
  * Yan et al. [2022] CA Yan, R Vescovini, and L Dozio.  A framework based on physics-informed neural networks and extreme learning for the analysis of composite structures.  _Computers & Structures_, 265:106761, 2022.
  * Niaki et al. [2021] Sina Amini Niaki, Ehsan Haghighat, Trevor Campbell, Anoush Poursartip, and Reza Vaziri.  Physics-informed neural network for modelling the thermochemical curing process of composite-tool systems during manufacture.  _Computer Methods in Applied Mechanics and Engineering_ , 384:113959, 2021.
  * Wessels et al. [2020] Henning Wessels, Christian Weißenfels, and Peter Wriggers.  The neural particle method–an updated lagrangian physics informed neural network for computational fluid dynamics.  _Computer Methods in Applied Mechanics and Engineering_ , 368:113127, 2020.
  * Yu et al. [2018] Bing Yu et al.  The deep ritz method: a deep learning-based numerical algorithm for solving variational problems.  _Communications in Mathematics and Statistics_ , 6(1):1–12, 2018.
  * Liao and Ming [2019] Yulei Liao and Pingbing Ming.  Deep nitsche method: Deep ritz method with essential boundary conditions.  _arXiv preprint arXiv:1912.01309_ , 2019.
  * Segerlind [1984] Larry J Segerlind.  Applied finite element analysis.  1984\.
  * Reddy [2014] Junuthula Narasimha Reddy.  _An Introduction to Nonlinear Finite Element Analysis Second Edition: with applications to heat transfer, fluid mechanics, and solid mechanics_.  OUP Oxford, 2014.
  * Samaniego et al. [2020] Esteban Samaniego, Cosmin Anitescu, Somdatta Goswami, Vien Minh Nguyen-Thanh, Hongwei Guo, Khader Hamdia, X Zhuang, and T Rabczuk.  An energy approach to the solution of partial differential equations in computational mechanics via machine learning: Concepts, implementation and applications.  _Computer Methods in Applied Mechanics and Engineering_ , 362:112790, 2020.
  * Nguyen-Thanh et al. [2020] Vien Minh Nguyen-Thanh, Xiaoying Zhuang, and Timon Rabczuk.  A deep energy method for finite deformation hyperelasticity.  _European Journal of Mechanics-A/Solids_ , 80:103874, 2020\.
  * Fuhg and Bouklas [2022] Jan N Fuhg and Nikolaos Bouklas.  The mixed deep energy method for resolving concentration features in finite strain hyperelasticity.  _Journal of Computational Physics_ , 451:110839, 2022.
  * Abueidda et al. [2022a] Diab W Abueidda, Seid Koric, Erman Guleryuz, and Nahil A Sobh.  Enhanced physics-informed neural networks for hyperelasticity.  _arXiv preprint arXiv:2205.14148_ , 2022a.
  * Rezaei et al. [2022] Shahed Rezaei, Ali Harandi, Ahmad Moeineddin, Bai-Xiang Xu, and Stefanie Reese.  A mixed formulation for physics-informed neural networks as a potential solver for engineering problems in heterogeneous domains: comparison with finite element method.  _arXiv preprint arXiv:2206.13103_ , 2022.
  * Kollmann et al. [2020] Hunter T Kollmann, Diab W Abueidda, Seid Koric, Erman Guleryuz, and Nahil A Sobh.  Deep learning for topology optimization of 2d metamaterials.  _Materials & Design_, 196:109098, 2020.
  * Abueidda et al. [2020] Diab W Abueidda, Seid Koric, and Nahil A Sobh.  Topology optimization of 2d structures with nonlinearities using deep learning.  _Computers & Structures_, 237:106283, 2020.
  * Sosnovik and Oseledets [2019] Ivan Sosnovik and Ivan Oseledets.  Neural networks for topology optimization.  _Russian Journal of Numerical Analysis and Mathematical Modelling_ , 34(4):215–223, 2019.
  * Banga et al. [2018] Saurabh Banga, Harsh Gehani, Sanket Bhilare, Sagar Patel, and Levent Kara.  3d topology optimization using convolutional neural networks.  _arXiv preprint arXiv:1808.07440_ , 2018.
  * Zhang et al. [2021] Zeyu Zhang, Yu Li, Weien Zhou, Xiaoqian Chen, Wen Yao, and Yong Zhao.  Tonr: An exploration for a novel way combining neural network with topology optimization.  _Computer Methods in Applied Mechanics and Engineering_ , 386:114083, 2021.
  * Hoyer et al. [2019] Stephan Hoyer, Jascha Sohl-Dickstein, and Sam Greydanus.  Neural reparameterization improves structural optimization.  _arXiv preprint arXiv:1909.04240_ , 2019.
  * Chandrasekhar and Suresh [2021] Aaditya Chandrasekhar and Krishnan Suresh.  Tounn: topology optimization using neural networks.  _Structural and Multidisciplinary Optimization_ , 63(3):1135–1149, 2021.
  * Zehnder et al. [2021] Jonas Zehnder, Yue Li, Stelian Coros, and Bernhard Thomaszewski.  Ntopo: Mesh-free topology optimization using implicit neural representations.  _Advances in Neural Information Processing Systems_ , 34:10368–10381, 2021.
  * Pattanayak et al. [2017] Santanu Pattanayak, John S Pattanayak, and Suresh John.  _Pro deep learning with tensorflow_.  Springer, 2017.
  * Zhu et al. [1997] Ciyou Zhu, Richard H Byrd, Peihuang Lu, and Jorge Nocedal.  Algorithm 778: L-bfgs-b: Fortran subroutines for large-scale bound-constrained optimization.  _ACM Transactions on mathematical software (TOMS)_ , 23(4):550–560, 1997.
  * Long [2021] Joshua M. Long.  Random fourier features pytorch.  _GitHub. Note: https://github.com/jmclong/random-fourier-features-pytorch_ , 2021.
  * Wang et al. [2021] Sifan Wang, Hanwen Wang, and Paris Perdikaris.  On the eigenvector bias of fourier feature networks: From regression to solving multi-scale pdes with physics-informed neural networks.  _Computer Methods in Applied Mechanics and Engineering_ , 384:113938, 2021.
  * Abueidda et al. [2022b] Diab W Abueidda, Seid Koric, Rashid Abu Al-Rub, Corey M Parrott, Kai A James, and Nahil A Sobh.  A deep learning energy method for hyperelasticity and viscoelasticity.  _European Journal of Mechanics-A/Solids_ , 95:104639, 2022b.
  * Xia and Breitkopf [2015] Liang Xia and Piotr Breitkopf.  Design of materials using topology optimization and energy-based homogenization approach in matlab.  _Structural and multidisciplinary optimization_ , 52(6):1229–1241, 2015.
  * Rozvany [2009] George IN Rozvany.  A critical review of established methods of structural topology optimization.  _Structural and multidisciplinary optimization_ , 37(3):217–237, 2009.
  * Zhang et al. [2019] Yan Zhang, Mi Xiao, Hao Li, and Liang Gao.  Topology optimization of material microstructures using energy-based homogenization method under specified initial material layout.  _Journal of Mechanical Science and Technology_ , 33(2):677–693, 2019.
  * Svanberg [1987] Krister Svanberg.  The method of moving asymptotes—a new method for structural optimization.  _International journal for numerical methods in engineering_ , 24(2):359–373, 1987.
  * Bruns and Tortorelli [2001] Tyler E Bruns and Daniel A Tortorelli.  Topology optimization of non-linear elastic structures and compliant mechanisms.  _Computer methods in applied mechanics and engineering_ , 190(26-27):3443–3459, 2001.
  * Chadha et al. [2022] Charul Chadha, Diab Abueidda, Seid Koric, Erman Guleryuz, and Iwona Jasiuk.  Optimizing hyperparameters and architecture of deep energy method.  2022\.
  * Paszke et al. [2019] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems 32_ , pages 8024–8035. Curran Associates, Inc., 2019.  URL http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf.
  * Chandrasekhar et al. [2021] 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.
  * [52] Topology optimization codes written in python.  https://www.topopt.mek.dtu.dk/Apps-and-software/Topology-optimization-codes-written-in-Python.  Accessed: 2022-06-25.
  * Andreassen et al. [2011] Erik Andreassen, Anders Clausen, Mattias Schevenels, Boyan S Lazarov, and Ole Sigmund.  Efficient topology optimization in matlab using 88 lines of code.  _Structural and Multidisciplinary Optimization_ , 43(1):1–16, 2011.
  * SIMULIA [2020] SIMULIA.  Abaqus, 2020.
  * MATLAB [2021] MATLAB.  _version R2021a_.  The MathWorks Inc., Natick, Massachusetts, 2021.
  * Yu et al. [2022] Jeremy Yu, Lu Lu, Xuhui Meng, and George Em Karniadakis.  Gradient-enhanced physics-informed neural networks for forward and inverse pde problems.  _Computer Methods in Applied Mechanics and Engineering_ , 393:114823, 2022.


#  From Observability to Observer Realization
A path via elementary block-diagram manipulations

Eder Baron-Prada, Renzo Caballero, Eric Feron

Introductory state-space linear control courses focus on linear, time-
invariant systems and spend intense efforts by introducing system realizations
that allow the student to grasp fundamental concepts, among which
controllability, observability, and controller and observer design. Assuming
that the system is single-output, observability realizations allow the student
to easily perceive the limitations of the observed output in terms of what it
can estimate about the state of the system. Once observability is established,
the observer realization provides a rich framework allowing the educator to
establish the theorem stating the equivalence between observability and the
ability to choose arbitrary pole locations for the corresponding state
estimation error dynamics. The pole location is achieved by introducing a
linear gain which leads to a trivial observer pole placement process. While
the observer realizations are trivial to obtain from systems expressed as
transfer functions, such is not the case when starting from arbitrary
observable state-space representations and elementary but laborious matrix
manipulations transform the observability realization of a given system into
its observer realization. This note describes a graphical mechanism to
transform a system expressed in observability form into its observer form
based on elementary block diagram manipulations, a technique whose pedagogical
usefulness tends to be overlooked[1]. Given the well known duality principle
in linear systems [2][Section 5.2.4], [3][Section 4.3] with respect to
controllability and observability, the proposed graphical mechanism is
applicable to transform from controllability to controller realization.

An extensive bibliographic search reveals references that briefly elaborate on
observability and observer properties of the systems [4, 5, 6, 7, 8, 9, 10],
however these references do not explore these concepts in depth. Additional
references introduce the observability matrix from the original system which
is essential to construct the observability realization [11, 12, 13, 14, 15,
16, 17, 3, 18, 19, 20, 21, 22, 23, 24, 25, 2, 26]. Finally, other references
describe the block diagrams for the observability/observer and
controllability/controllable canonical realizations [27, 28, 18, 29, 30],
however they rely on the transformation theorem to reach the observer
realization from the observability form. The same search indicates that the
graphical approach, introduced here, to perform the state-space transformation
from observability realization to observer realization, is new. Further
investigations about the state-space transformation corresponding to these
block manipulations also reveal a close connection between the computation of
these state-space transformations and recent results on generalized Fibonacci
polynomials.

## Preliminaries

We consider the single-input, single-output system

$\begin{array}[]{rcl}\displaystyle\frac{\operatorname{d\\!}}{\operatorname{d\\!}t}\bm{x}&=&\bm{Ax}+\bm{B}u,\;\bm{x}(0)=\bm{x}_{0},\\\\[10.0pt]
y&=&\bm{Cx},\end{array}$ (1)

where $\bm{x}=\bm{x}(t)\in{\bf R}^{n}$ is the system’s state, $u=u(t)$ is the
system’s input, and $y=y(t)$ is the system’s output. As we are addressing the
observability characteristics of the system, for notation clarity purposes, we
assume a null input to the system and therefore omit the input matrix
$\bm{B}$. However, the reader must be warned that the following discussion
makes extensive use of the word “realization”, which is most often understood
as applying to systems with both inputs and outputs. Thus, although the input
matrix $\bm{B}$ is omitted below, it is always present in the discussions
below whenever the word “realization” is used.

All considerations are relative to the system’s observability properties and
observer design considerations. These considerations can be transposed to the
corresponding controllability and controller design considerations. A simple
criterion to establish the observability of the system (1) is established by
the classical observability theorem.

######  Theorem 1 (See [29]).

The system (1) is observable if and only if the square matrix

${\cal O}=\left[\begin{array}[]{l}\bm{C}\\\ \bm{CA}\\\ \bm{CA}^{2}\\\ \hskip
5.69054pt\vdots\\\ \bm{CA}^{n-1}\end{array}\right]$

is invertible.

The observability theorem exposed relates the linear independency of the $n$
rows of $\mathcal{O}$, with the ability to estimate each one of the $n$ states
of the system (1). The observability matrix, $\mathcal{O}$, is used to extract
the observability realization of the observable system (1).

######  Theorem 2.

Consider the basis formed by the vectors
$\bm{C}^{\top},\;\bm{A}^{\top}\bm{C}^{\top},\;\ldots,{\bm{A}^{n-1}}^{\top}\bm{C}^{\top}$.
In that basis, the system (1) is realized with the companion form

$\begin{array}[]{rcl}\displaystyle\frac{\operatorname{d\\!}}{\operatorname{d\\!}t}\tilde{\bm{x}}&=&\left[\begin{array}[]{ccccc}0&1&0&\ldots&0\\\
0&0&1&\ddots&\vdots\\\ \vdots&&\ddots&\ddots&\vdots\\\ 0&0&\ldots&0&1\\\
-a_{0}&-a_{1}&\ldots&-a_{n-2}&-a_{n-1}\end{array}\right]\tilde{\bm{x}}\\\\[40.0pt]
&=&\bm{A}_{\rm observability}\tilde{\bm{x}}\\\\[10.0pt]
y&=&\left[1\;0\;\ldots\;0\;0\right]\tilde{\bm{x}}\\\\[10.0pt] &=&\bm{C}_{\rm
observability}\tilde{\bm{x}},\end{array}$ (2)

with ${\bm{A}}_{\rm observability}={\cal O}\bm{A}{\cal O}^{-1}$,
${\bm{C}}_{\rm observability}=\bm{C}{\cal O}$, and $a_{0},\dots,\;a_{n-1}$ are
the coefficients of the characteristic polynomial $\lambda(s)$ of $\bm{A}$

$\lambda(s)=s^{n}+a_{n-1}s^{n-1}+\dots+a_{1}s+a_{0}.$ (3)

###### Proof.

Proof of Theorem 2 is well know in the literature. See for instance
[2][Section 5.2.2], [29][Section 6.2] and [13] [Section 3.2]. ∎

While the observability realization (2) and surrounding analyses provide
valuable and easy insight into the observability of (1), it fails to provide
an immediately constructive mechanism to specify observers in terms of
specifying arbitrary poles of observer error dynamics. A realization of (1)
that supports this purpose is the observer realization

$\begin{array}[]{rcl}\displaystyle\frac{\operatorname{d\\!}}{\operatorname{d\\!}t}\bar{\bm{x}}&=&\left[\begin{array}[]{ccccc}-a_{n-1}&1&0&\ldots&0\\\
-a_{n-2}&0&1&\ddots&\vdots\\\ \vdots&&\ddots&\ddots&\vdots\\\
-a_{1}&0&\ldots&0&1\\\
-a_{0}&0&0&\ldots&0\end{array}\right]\bar{\bm{x}}\\\\[40.0pt] &=&\bm{A}_{\rm
observer}\bar{\bm{x}}\\\\[10.0pt]
y&=&\left[1\;0\;\ldots\;0\;0\right]\bar{\bm{x}}\\\\[10.0pt] &=&\bm{C}_{\rm
observer}\bar{\bm{x}}.\end{array}$ (4)

This realization is convenient for the reader to understand the flexibility
offered by linear observers of the form

$\begin{array}[]{rcl}\displaystyle\frac{\operatorname{d\\!}}{\operatorname{d\\!}t}\hat{\bm{x}}&=&\bm{A}\hat{\bm{x}}+\bm{L}(y-\bm{C}\hat{\bm{x}})\\\
&=&(\bm{A}-\bm{LC})\hat{\bm{x}}+\bm{L}y.\end{array}$

This kind of observer is known as Luenberger observer, where $\bm{L}\in{\bf
R}^{n}$ is the observer gain. Indeed, in the transformed state-space that
defines $\bar{x}$, the observer becomes

$\begin{array}[]{rcl}\displaystyle\frac{\operatorname{d\\!}}{\operatorname{d\\!}t}\hat{\bar{\bm{x}}}&=&\bm{A}_{\rm
observer}\hat{\bar{\bm{x}}}+\bm{L}_{\rm observer}(y-\bm{C}_{\rm
observer}\hat{\bar{\bm{x}}})\\\\[10.0pt] &=&(\bm{A}_{\rm observer}-\bm{L}_{\rm
observer}\bm{C}_{\rm observer})\hat{\bar{\bm{x}}}+\bm{L}_{\rm
observer}y\\\\[10.0pt]
&=&\left[\begin{array}[]{ccccc}-a_{n-1}-l_{n-1}&1&0&\ldots&0\\\
-a_{n-2}-l_{n-2}&0&1&\ddots&\vdots\\\ \vdots&&\ddots&\ddots&\vdots\\\
-a_{1}-l_{1}&0&\ldots&0&1\\\
-a_{0}-l_{0}&0&0&\ldots&0\end{array}\right]\hat{\bar{\bm{x}}}+\left[\begin{array}[]{l}l_{n-1}\\\
l_{n-2}\\\ \hskip 2.84526pt\vdots\\\ l_{1}\\\
l_{0}\end{array}\right]y.\end{array}$

This realization clearly shows how the observer gains can be used to modify
the characteristic polynomial of the estimation error dynamics arbitrarily
since each gain directly and simply impacts each corresponding characteristic
polynomial coefficient. The derivation of the observer realization (4) of the
system (1) from its observability realization (2) is, however, not immediate
and existing derivations, which rely on a number of linear-algebraic
constructs, can constitute an obstacle to the senior undergraduate or first-
year graduate students who learns this topic.

We have found a way to go from Observability to Observer realization in the
current literature. It relies on the transfer function for the derivation of
the observer realization form. This is immediate but relies on having the
system’s transfer function. We propose two alternative derivations of the
observer realization via matrix manipulations and elementary block
manipulations that follow the spirit of those found in Kailath [29, pp.
37-45].

## Going from Observability to Observer realizations via Matrix Manipulations

Given the system (1), a linear state-space transformation is given by any
invertible matrix $\bm{P}\in{\bf R}^{n\times n}$ such that
$\hat{\bar{\bm{x}}}=\bm{Px}$, resulting in the new realization

$\displaystyle\frac{\operatorname{d\\!}}{\operatorname{d\\!}t}\hat{\bar{\bm{x}}}$
$\displaystyle=\bm{PAP}^{-1}\hat{\bar{\bm{x}}}$ (5) $\displaystyle y$
$\displaystyle=\bm{CP}^{-1}\hat{\bar{\bm{x}}}.$ (6)

Considering that
$\bm{A}_{\rm{observer}}=\bm{PA}_{\rm{observability}}\bm{P}^{-1}$, acceptable
linear state-space transformations may be found by solving for (7) and (8).

$\displaystyle\bm{A}_{\rm{observer}}\bm{P}-\bm{PA}_{\rm{observability}}=0,$
(7) $\displaystyle\bm{C}_{\rm{observer}}\bm{P}-\bm{C}_{\rm{observability}}=0.$
(8)

The complete transformation of the realization also involves the condition
$\bm{B}_{\rm observer}-\bm{PB}_{\rm observability}=0$. However, this condition
does not participate actively to the rest of the paper and is therefore
omitted.

Equation (7) is a homogeneous Sylvester equation in the unknown $\bm{P}$. The
equivalence of $\bm{A}_{\rm{observer}}$ and $\bm{A}_{\rm{observability}}$
implies that there exists an invertible matrix $\bm{P}$ that satisfies (7) and
(8). One such solution is

$\bm{P}=\left[\begin{array}[]{cccc}1&0&\cdots&0\\\ a_{n-1}&1&\ddots&\vdots\\\
\vdots&\ddots&\ddots&0\\\ a_{0}&\cdots&a_{n-1}&1\end{array}\right].$

$\bm{P}$ is a Toeplitz matrix of dimension $n\times n$. It is remarkable, that
as it was expected all the eigenvalues of the system are equal to one, given
that it means that the spectrum is not changed by doing this transformation.
The matrix $\bm{P}^{-1}$ is characterized by [31]

$\displaystyle\bm{P}^{-1}=\left[\begin{array}[]{cccc}1&0&\cdots&0\\\
F_{1}&1&\ddots&\vdots\\\ \vdots&\ddots&\ddots&0\\\
F_{n-1}&\cdots&F_{1}&1\end{array}\right],$ (13)

where $F_{k}$ is the generalized Fibonacci polynomial of order $k$ on the
variables $(a_{n-k},a_{n-k+1},a_{n-k+2},\dots,a_{n-1})$. The generalized
Fibonacci polynomials ($F_{k}$), with $k\in\\{0,\ldots,n-1\\}$, are defined
inductively by

$\displaystyle F_{0}$ $\displaystyle=1$ $\displaystyle F_{1}$
$\displaystyle=-a_{n-1}F_{0}$ $\displaystyle F_{2}$
$\displaystyle=-a_{n-1}F_{1}-a_{n-2}F_{0}$ $\displaystyle F_{3}$
$\displaystyle=-a_{n-1}F_{2}-a_{n-2}F_{1}-a_{n-3}F_{0}$ $\displaystyle\vdots$
$\displaystyle F_{k}$ $\displaystyle=\sum_{i=1}^{k}-a_{n-i}F_{k-i}$
$\displaystyle\vdots$ $\displaystyle F_{n-1}$
$\displaystyle=-a_{n-1}F_{n-2}-a_{n-2}F_{n-3}-\cdots-a_{2}F_{1}-a_{1}F_{0}.$

Each $F_{k}$ may also be calculated as the determinant of the lower Hessenberg
matrix $\bm{H}_{k}$ [31], with

$\bm{H}_{k}=\left[\begin{array}[]{ccccc}a_{n-1}&-1&0&\cdots&0\\\
a_{n-2}&a_{n-1}&-1&\cdots&0\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
a_{n-k+1}&a_{n-k+2}&a_{n-k+3}&\cdots&-1\\\
a_{n-k}&a_{n-k+1}&a_{k-2}&\cdots&a_{n-1}\end{array}\right].$

Although interesting in its own right, and apparently new to the control
community, the explicit computation of $\bm{P}$ and its inverse $\bm{P}^{-1}$,
the transformation matrix allowing to go from observability to observer
realizations is heavy with computations and not very friendly to the novice.
This is why we introduce a mechanical and graphical process to do the same by
means of a sequence of elementary block manipulations.

## Transforming the observability realization into the observer realization
via elementary block manipulations

We begin with the observability realization (2), which we represent
graphically in Fig. 1.

Figure 1: Block-diagram representation of the observability realization of the
system (1).
---
$y^{(n)}+y^{(n-1)}a_{n-1}+\dots+a_{0}y=0.$ (14)
Table 1: Block-diagram representation of the observability realization of the
system (1) with the output $y$ ODE.

The state-space transformation to the observer realization is obtained by
moving integrators around the block diagram. With this proposed method, the
system is transformed into several observable and non-observable realizations
in the described steps.

The first move is shown in Fig. 2, whereby the upper left integrator is moved
down, replacing $\operatorname{d\\!}^{n}y/\operatorname{d\\!}t^{n}$ by
$\operatorname{d\\!}^{n-1}y/\operatorname{d\\!}t^{n-1}$ on the way.

Figure 2: First move: The upper left integrator is moved down.
---
$y^{(n-1)}={\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\int}y^{(n-1)}(-a_{n-1}){\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}}+{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\int}\left[\sum_{i=1}^{n-1}(-a_{i-1})\int^{(n-i)}y^{(n-1)}\right]{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}}$
$y^{(n-1)}=y^{(n-2)}(-a_{n-1})+{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\int}\left[\sum_{i=1}^{n-1}(-a_{i-1})y^{(i-1)}\right]{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}}$
$y^{(n-1)}=y^{(n-2)}(-a_{n-1})+\sum_{i=1}^{n-1}(-a_{i-1})y^{(i-2)}$
$y^{(n-1)}=\sum_{i=1}^{n}(-a_{i-1})y^{(i-2)}$ So, we obtain an integrated
version of ODE 14 in TABLE 1:
$y^{(n-1)}+a_{n-1}y^{(n-2)}+\dots+a_{1}y+a_{0}\int y=0.$
Table 2: First move: The upper left integrator is moved down. We also observe
the associated ODE.

The movement of the same integrator behind the sum splits it due to linearity.
Also, the following upper left integrator is moved down, where a
differentiator is created next to the gain $a_{n-1}$. We observe these two
steps in Fig. 3.

Figure 3: Second and third move: The lower left integrator is moved before the
sum, where it is split due to linearity (second move). The next upper left
integrator is moved down, where a differentiator is created next to the gain
$a_{n-1}$ (third move).
---
$y^{(n-2)}={\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\int}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\int}{\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\frac{\operatorname{d\\!}}{\operatorname{d\\!}t}}\left((-a_{n-1})y^{(n-2)}\right){\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}}{\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}}\\\
+{\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\int}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\int}\left((-a_{n-2})y^{(n-2)}\right){\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}}{\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}}+{\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\int}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\int}\left(\sum_{i=3}^{n}\int^{(i-2)}(-a_{i-1})y^{(n-1)}\right){\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}}{\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}}$
$y^{(n-2)}=(-a_{n-1})y^{(n-3)}+(-a_{n-2})y^{(n-4)}+\sum_{i=3}^{n}(-a_{n-i})y^{(n-2-i)}$
$y^{(n-2)}=\sum_{i=1}^{n}(-a_{n-i})y^{(n-2-i)}.$ So, we obtain a two-times
integrated version of ODE 14: $y^{(n-2)}+a_{n-1}y^{(n-3)}+\dots+a_{1}\int
y+a_{0}\int\int y=0.$
Table 3: Second and third move: The lower left integrator is moved before the
sum, where it is split due to linearity (second move). The next upper left
integrator is moved down, where a differentiator is created next to the gain
$a_{n-1}$ (third move). We also observe the associated ODE.

Since blocks $s$ and $1/s$ next to gain $a_{n-1}$ represent inverse linear
operators, we cancel them and we reach the realization in Fig. 4.

Figure 4: Realization after canceling integrator and differentiator next to
gain $a_{n-1}$.
---
$y^{(n-2)}=\int(-a_{n-1})y^{(n-2)}+\int\int(-a_{n-2})y^{(n-2)}\\\
+\int\int\left(\sum_{i=3}^{n}\left(-a_{n-i}\int^{(i-2)}y^{(n-2)}\right)\right)$
$y^{(n-2)}=\sum_{i=1}^{n}\int^{(i)}(-a_{n-i})y^{(n-1)}.$ So, we obtain a two-
times integrated version of ODE 14:
$y^{(n-2)}+a_{n-1}y^{(n-3)}+\dots+a_{1}\int y+a_{0}\int\int y=0.$
Table 4: Realization and ODE for the output after canceling the integrator and
differentiator next to gain $a_{n-1}$.

We keep repeating the preciously described process a total of $m$ times, where
$m$ integrator are down, and $n-m$ are up. Also, there are no integrators nor
differentiators next to gains $a_{n-1},a_{n-2},\dots,a_{n-(m-1)}$ due to
cancelation. In Fig. 5 we observe the associated realization.

Figure 5: We observe the realization after moving down and aligning the first
$m$ integrators.
---
$y^{(n-m)}=\sum_{i=1}^{m}(-a_{n-i})\int^{(i)}y^{(n-m)}+\int^{(m)}\sum_{i=1}^{n-m}\int^{(i)}(-a_{n-m-i})y^{(n-m)}$
$y^{(n-m)}=\sum_{i=1}^{n}\int^{(i)}(-a_{n-i})y^{(n-m)}.$ So, we obtain a
$m$-times integrated version of ODE 14:
$y^{(n-m)}+a_{n-1}y^{(n-m-1)}+\dots+a_{1}\int^{(m-1)}y+a_{0}\int^{(m)}y=0.$
Table 5: Realization, state-space matrix, and output ODE for the system after
moving down $m$ integrators.

It is now possible to repeat the foregoing steps shown in the Fig. 2, until
Fig. 5 and progressively move all remaining upper integrators until none
remains, resulting in the observer realization shown in Fig. 6.

Figure 6: Observer realization.
---
$y+a_{n-1}\int y+a_{n-2}\int\int y+\dots+a_{0}\int^{(n)}y=0.$ (15)
Table 6: Observer realization.

In the next section we analyze the block manipulations as a succession of
elementary state-space transformations.

Before leaving this section, we would like to thank the reviewer for
suggesting to draw the equivalence between the suggested manipulations and the
progressive integration of the differential equation (14) to yield the
integral equation (15). We utilize $\int^{(n)}y$ to denote the process of
integrating $n$ times the function $y$ with respect to time.

## State-space interpretation of elementary block manipulations

We call $\bm{A}_{m}$ to the state-matrix associated to the $m$-th observable
realization during the block manipulation. We call block manipulation to the
process of going from Table 5 for step $m$ to Table 5 for step $m+1$.
Essentially, we start with $\bm{A}_{\rm{observability}}$, and we move through
the sequence $\bm{A}_{1},\dots,\bm{A}_{n-1}$, until we reach
$\bm{A}_{\rm{observer}}$. $\bm{A}_{m}$ is given by

$\bm{A}_{m}=\setcounter{MaxMatrixCols}{14}\begin{bmatrix}\textbf{0}&1&0&\dots&\dots&\dots&\dots&\dots&\dots&\dots&\dots&0&0&0\\\
0&\textbf{0}&1&0&\dots&\dots&\dots&\dots&\dots&\dots&\dots&0&0&0\\\
\dots&\dots&\dots&\dots&\dots&\dots&\dots&\dots&\dots&\dots&\dots&\dots&\dots&\dots\\\
0&0&0&\dots&\textbf{0}&1&0&\dots&\dots&\dots&\dots&0&0&0\\\
0&0&0&\dots&\dots&\textbf{0}&1&0&\dots&\dots&\dots&0&0&0\\\
0&0&0&\dots&\dots&\dots&\bm{-a_{n-1}}&1&0&\dots&\dots&0&0&0\\\
0&0&0&\dots&\dots&\dots&-a_{n-2}&\textbf{0}&1&0&\dots&0&0&0\\\
\dots&\dots&\dots&\dots&\dots&\dots&\dots&\dots&\dots&\dots&\dots&\dots&\dots&\dots\\\
0&0&0&\dots&\dots&\dots&-a_{n-(m+1)}&\dots&\dots&\dots&\dots&0&\textbf{0}&1\\\
-a_{0}&-a_{1}&-a_{2}&\dots&\dots&-a_{n-(m-1)}&-a_{n-m}&0&\dots&\dots&\dots&\dots&\dots&\textbf{0}\\\
\end{bmatrix}.$

Since all observable realizations have the same dynamics, i.e., they are
equivalent. There exists a change-of-basis matrix $\bm{P}_{m}$ such that
$\bm{A}_{m}=\bm{P}_{m}\bm{A}_{m-1}\bm{P}^{-1}_{m}$ leading to

$\bm{P}_{n-1}...\bm{P}_{1}\bm{A}_{\rm{observability}}\bm{P}^{-1}_{1}...\bm{P}^{-1}_{n-1}=\bm{A}_{\rm{observer}}.$
(16)

The change-of-basis matrix $\bm{P}_{i}$ is given by

$\bm{P}_{i}=\begin{bmatrix}1&0&0&\dots&\dots&\dots&\dots&0&0\\\
0&\ddots&\ddots&&&&&&0\\\ \vdots&\ddots&1&\ddots&&&&&\vdots\\\
\vdots&&0&1&\ddots&&&&\vdots\\\ \vdots&&\vdots&-a_{n-1}&1&\ddots&&&\vdots\\\
\vdots&&\vdots&-a_{n-2}&0&\ddots&\ddots&&\vdots\\\
\vdots&&\vdots&\vdots&\vdots&\ddots&\ddots&\ddots&0\\\
\vdots&&\vdots&-a_{n-(i-1)}&\vdots&&\ddots&\ddots&0\\\
0&\dots&0&-a_{n-i}&0&\dots&\dots&0&1\\\ \end{bmatrix},$

and its inverse $\bm{P}_{i}^{-1}$ by

$P_{i}^{-1}=\begin{bmatrix}1&0&0&\dots&\dots&\dots&\dots&0&0\\\
0&\ddots&\ddots&&&&&&0\\\ \vdots&\ddots&1&\ddots&&&&&\vdots\\\
\vdots&&0&1&\ddots&&&&\vdots\\\ \vdots&&\vdots&a_{n-1}&1&\ddots&&&\vdots\\\
\vdots&&\vdots&a_{n-2}&0&\ddots&\ddots&&\vdots\\\
\vdots&&\vdots&\vdots&\vdots&\ddots&\ddots&\ddots&0\\\
\vdots&&\vdots&a_{n-(i-1)}&\vdots&&\ddots&\ddots&0\\\
0&\dots&0&a_{n-i}&0&\dots&\dots&0&1\\\ \end{bmatrix},$

where the coefficients $a_{n-1},\ldots,a_{n-i}$ are placed on the column
$n-i$. Naturally, (16) leads to the matrices $\bm{P}$ and $\bm{P}^{-1}$ in
(Going from Observability to Observer realizations via Matrix Manipulations)
and (13) respectively, as

$\displaystyle\prod_{i=1}^{n-1}\bm{P}_{n-i}$ $\displaystyle=\bm{P}$
$\displaystyle\prod_{i=1}^{n-1}\bm{P}^{-1}_{i}$ $\displaystyle=\bm{P}^{-1}.$

Furthermore, given that $\bm{C}_{m}=\bm{C}_{m-1}\bm{P}^{-1}_{m}$, it is
possible to see that
$\bm{C}_{\rm{observability}}=\bm{C}_{m}=\bm{C}_{\rm{observer}}$.

Despite the intrinsic value of these elementary state-space transformations,
the authors hope that the reader will agree with them when they say the
equivalent block-diagram manipulations are far more accessible to the senior
undergraduate or a first-year graduate student interested in learning about
linear systems and navigating their canonical state-space representations.

## Conclusion

A visual and easy-to-understand transformation from a LTI system in
observability form to observer form is introduced. This transformation solely
relies on elementary block diagram manipulations, thus making it accessible to
students during their early exposure to state-space systems. The third author
author has made this transformation a part of his elementary linear systems
theory curriculum. It also applies to controllability and controller
realizations via duality.

## Author Biography

Eder Baron-Prada received a B.S. degree in Electrical Engineering and the M.S.
degree in Industrial Automation from the Universidad Nacional de Colombia,
Bogotá, Colombia, in 2016 and 2019 respectively. He is currently a Ph.D.
student at King Abdullah University of Science and Technology (KAUST) in
Electrical and Computer Engineering. His research interests include game
theory, optimization, and control of multi-agent systems applied to modern
power systems.

Renzo Caballero received an Electrical Engineering degree from Universidad de
la República, Uruguay, in 2016, and an M.S. in Applied Mathematics and
Computer Science degree from King Abdullah University of Science and
Technology (KAUST), KSA, in 2019. He is currently pursuing his PhD at KAUST.
His current research focuses on self-replicating robots, but he also is
interested in microelectronic design, synthesis of 2D materials, and
stochastic optimal control.

Eric Feron is a professor of Electrical, Computer, and Mechanical Engineering.
He is the director of the Robotics, Intelligent Systems, and Control (RISC)
Laboratory. He recently joined the KAUST CEMSE Division from the Georgia
Institute of Technology. Prior to his time at Georgia Tech, he was an active
faculty member in MIT’s Aeronautics and Astronautics department from 1993
until 2005. Feron’s career in academia began in Paris, France, where he
obtained the B.S. and M.S. from École Polytechnique and École Normale
Supérieure, respectively. He later completed the Ph.D. in aerospace
engineering at Stanford University, USA.

## References

  * [1] R. C. Dorf and R. H. Bishop, _Modern control systems_. Pearson Prentice Hall, 2008.
  * [2] A. Bacciotti, _Stability and Control of Linear Systems_. Springer, 2019.
  * [3] R. L. Williams, D. A. Lawrence _et al._ , _Linear state-space control systems_. Wiley Online Library, 2007.
  * [4] J. RUBIO, Ed., _The Theory of Linear Systems_. Academic Press, 1971. [Online]. Available: https://www.sciencedirect.com/science/article/pii/B9780126016505500021
  * [5] R. W. Brockett, _Finite dimensional linear systems_. SIAM, 2015.
  * [6] C. Heij, A. C. Ran, and F. van Schagen, _Introduction to mathematical systems theory: linear systems, identification and control_. Springer Science & Business Media, 2006.
  * [7] J. P. Hespanha, _Linear systems theory_. Princeton university press, 2018.
  * [8] R. Curfoin, “Infinite-dimensional linear systems theory,” 1978.
  * [9] L. Farina and S. Rinaldi, _Positive linear systems: theory and applications_. John Wiley & Sons, 2011, vol. 50.
  * [10] O. H. Bosgra, H. Kwakernaak, and G. Meinsma, “Design methods for control systems,” _Notes for a Course of the Dutch Institute of Systems and Control, Winter term_ , vol. 2002, 2001.
  * [11] H. L. Trentelman, A. A. Stoorvogel, and M. Hautus, _Control theory for linear systems_. Springer Science & Business Media, 2012.
  * [12] D. Thyagarajan, T.; Kalpana, _Linear and Non-Linear System Theory_. Taylor & Francis Group, 2020.
  * [13] P. J. Antsaklis and A. N. Michel, _Linear systems_. Springer Science & Business Media, 2006.
  * [14] A. Sinha, _Linear systems: optimal and robust control_. CRC press, 2007.
  * [15] G. Gu, _Discrete-time linear systems: theory and design with applications_. Springer Science & Business Media, 2012.
  * [16] J. O’Reilly, _Observers for linear systems_. Academic press, 1983, vol. 170.
  * [17] B. Friedland, _Control System Design: An Introduction to State-Space Methods_ , ser. Dover Books on Electrical Engineering. Dover Publications, 2005.
  * [18] L. Zadeh and C. Desoer, _Linear system theory: the state space approach_. Courier Dover Publications, 2008\.
  * [19] C.-T. Chen and B. Shafai, _Linear system theory and design_. Oxford university press New York, 1999, vol. 3.
  * [20] H. Bourlès and G. K. Kwan, _Linear systems_. John Wiley & Sons, 2013.
  * [21] E. Hendricks, O. Jannerup, and P. H. Sørensen, _Linear systems control: deterministic and stochastic methods_. Springer, 2008.
  * [22] J. Zabczyk, _Mathematical control theory_. Springer, 2020.
  * [23] K. Ogata, _Modern control engineering_. Prentice hall, 2010.
  * [24] W. L. Brogan, _Modern control theory_. Pearson education india, 1991.
  * [25] A. Isidori, _Nonlinear control systems_. Springer Science & Business Media, 2013.
  * [26] D. F. Delchamps, _State space and input-output linear systems_. Springer Science & Business Media, 2012.
  * [27] G. Basile and G. Marro, _Controlled and conditioned invariants in linear system theory_. Prentice Hall Englewood Cliffs, NJ, 1992.
  * [28] G. F. Franklin, J. D. Powell, and A. Emami-Naeni, _Feedback Control of Dynamic Systems_. Addison-Wesley, 1986\.
  * [29] T. Kailath, _Linear Systems_. New Jersey: Prentice-Hall, 1980.
  * [30] F. Szidarovszky, _Linear systems theory_. Routledge, 2018.
  * [31] A. Şahin, “Inverse and factorization of triangular toeplitz matrices,” _Miskolc Mathematical Notes_ , vol. 19, no. 1, pp. 527–536, 2018.