Data-Driven Microstructure Property Relations

TL;DR

This paper develops a machine learning-based approach to predict the effective heat conductivity of heterogeneous microstructured materials from image data, emphasizing the use of 2-point correlation functions and incremental POD methods.

Contribution

It introduces a novel combination of snapshot POD and neural networks for efficient microstructure-property prediction from images.

Findings

Effective prediction accuracy demonstrated on synthetic microstructures

Proposed incremental POD methods reduce computational complexity

Open-source Python code available for implementation

Abstract

An image based prediction of the effective heat conductivity for highly heterogeneous microstructured materials is presented. The synthetic materials under consideration show different inclusion morphology, orientation, volume fraction and topology. The prediction of the effective property is made exclusively based on image data with the main emphasis being put on the 2-point spatial correlation function. This task is implemented using both unsupervised and supervised machine learning methods. First, a snapshot proper orthogonal decomposition (POD) is used to analyze big sets of random microstructures and thereafter compress significant characteristics of the microstructure into a low-dimensional feature vector. In order to manage the related amount of data and computations, three different incremental snapshot POD methods are proposed. In the second step, the obtained feature vector is used to predict the effective material property by using feed forward neural networks. Numerical examples regarding the incremental basis identification and the prediction accuracy of the approach are presented. A Python code illustrating the application of the surrogate is freely available.

Figures

Tables

Table 1. Table 1: Data of the unsupervised learning (incremental RB identification) for the nine considered scenarios; the parameters ε = 0.025 , n c = 100 ​  and  ​ n a = 75 formulae-sequence 𝜀 0.025 subscript 𝑛 c 100  and  subscript 𝑛 a 75 \varepsilon=0.025\,,n_{\rm c}=100\text{ and }n_{\rm a}=75 were used. Some numbers are rounded for easier readability.

Method	final basis size	snapshots with ${𝒫}_{\delta }\ge \epsilon$	snapshots with ${𝒫}_{\delta }\le \epsilon$	enrichment steps $n_{\mathrm{it}}$	time $[$s$]$	used microstructures
A	143	150	730	4	20
B	80	400	2400	7	70
C	96	800	7700	12	200
A	596	670	4500	11	150
B	294	2400	12700	34	500
C	312	2600	16500	37	550
A	464	560	2900	9	150
B	274	2000	16100	29	500
C	244	1540	8000	22	280

Table 2. Table 2: Percentage errors for κ ¯ 11 subscript ¯ 𝜅 11 \bar{\kappa}_{11} and κ ¯ 22 subscript ¯ 𝜅 22 \bar{\kappa}_{22} given for each of the best ANNs, validated on every test set.

			validated with
		circles		rectangles		mixed
trained with	error measures	${\overline{\kappa }}_{11}$	${\overline{\kappa }}_{22}$	${\overline{\kappa }}_{11}$	${\overline{\kappa }}_{22}$	${\overline{\kappa }}_{11}$	${\overline{\kappa }}_{22}$
circles	mean [%]	1.7	1.7	3.5	3.9	3.9	3.8
	max [%]	13.9	12.2	16.4	24.6	27.1	16.4
rectangles	mean [%]	3.4	3.6	1.9	1.9	4.2	3.7
	max [%]	19.5	19.9	8.9	10.9	24.2	21.0
mixed	mean [%]	3.4	3.6	2.6	2.6	2.7	2.4
	max [%]	13.7	19.9	12.8	15.6	15.6	14.7

Equations

$$f_a=1-f_b$$

$$c_{2}(\mbox{\boldmath$r$};a,b)=\bigl{\langle}\,\chi^{(a)}(\mbox{\boldmath$x$})\,\chi^{(b)}(\mbox{\boldmath$x$}+\mbox{\boldmath$r$})\bigl{\rangle}_{x}$$

$$\displaystyle c_{2}(\mbox{\boldmath$r$};a,b)$$

$$\displaystyle c_{2}(\mbox{\boldmath$r$};a,a)$$

$$\displaystyle c_{2}(\mbox{\boldmath$0$};a,a)$$

$$\displaystyle c_{2}(\mbox{\boldmath$0$};b,b)$$

$$\displaystyle c_{2}(\mbox{\boldmath$0$};a,b)$$

$$\displaystyle c_{2}(\mbox{\boldmath$r$};a,b)$$

$$\displaystyle\langle c_{2}(\mbox{\boldmath$x$};m,m)\rangle_{x}$$

$${\underline{s}}_i={\underline{c}}_{2_i}^0(b,b)-f_b^2\underline{1}={\underline{c}}_{2_i}^0(a,a)-f_a^2\underline{1}.$$

$${\mathcal{P}}_{\delta }=\frac{||\underline{\underline{S}}-\underline{\underline{B}}{\underline{\underline{B}}}^{\mathsf{T}}\underline{\underline{S}}|{|}_F}{||\underline{\underline{S}}|{|}_F}$$

$${\underline{\underline{C}}}_S={\underline{\underline{S}}}^{\mathsf{T}}\underline{\underline{S}}=\underline{\underline{V}}\underline{\underline{\Theta }}{\underline{\underline{V}}}^{\mathsf{T}}.$$

$${\underline{\underline{V}}}^{\mathsf{T}}\underline{\underline{V}}$$

$${\delta }_N=\sqrt{\frac{{\sum }_{j=N+1}^{n_s}{\theta }_j}{{\sum }_{i=1}^{n_s}{\theta }_i}}=\sqrt{\frac{{\sum }_{j=N+1}^{n_s}{\theta }_j}{||\underline{\underline{S}}|{|}_F^2}}=\sqrt{\Vert \underline{\underline{S}}{\Vert }_F^2-\sum _{j=1}^N{\theta }_j}\stackrel{!}{\le }\epsilon ,$$

$$\underline{\underline{B}}=\underline{\underline{S}}\widetilde{\underline{\underline{V}}}{\widetilde{\underline{\underline{\Theta }}}}^{-\frac{1}{2}}$$

$$\underline{\underline{S}}=\underline{\underline{U}}\underline{\underline{\Sigma }}{\underline{\underline{W}}}^{\mathsf{T}}$$

$$\underline{\underline{U}}\in {\mathbb{R}}^{n\times n_s}:{\underline{\underline{U}}}^{\mathsf{T}}\underline{\underline{U}}$$

$${\delta }_N$$

$$\underline{\underline{B}}=\widetilde{\underline{\underline{U}}}.$$

$${\mathcal{P}}_{\delta }=\frac{||{\underline{s}}_i-\underline{\underline{B}}{\underline{\underline{B}}}^{\mathsf{T}}{\underline{s}}_i|{|}_F}{||{\underline{s}}_i|{|}_F}.$$

$$\underline{\underline{\Delta \hat{S}}}$$

$$\underline{\underline{\Delta C}}$$

$$\underline{\underline{\Delta B}}$$

$$\underline{\underline{B}}$$

$$\frac{\Vert \underline{\underline{\Delta \hat{S}}}-\underline{\underline{\Delta B}}{\underline{\underline{\Delta B}}}^{\mathsf{T}}\underline{\underline{\Delta \hat{S}}}{\Vert }_F}{\Vert \underline{\underline{\Delta S}}{\Vert }_F}$$

$$\text{max}\frac{|{\underline{s}}_i-\underline{\underline{B}}{\underline{\underline{B}}}^{\mathsf{T}}{\underline{s}}_i|}{|s_i|}$$

$$\underline{\underline{C}}=\left[\begin{array}{cc}{\underline{\underline{S}}}^{\mathsf{T}}\underline{\underline{S}}& {\underline{\underline{S}}}^{\mathsf{T}}\underline{\underline{\Delta S}}\\ {\underline{\underline{\Delta S}}}^{\mathsf{T}}\underline{\underline{S}}& {\underline{\underline{\Delta S}}}^{\mathsf{T}}\underline{\underline{\Delta S}}\end{array}\right]=\left[\begin{array}{cc}{\underline{\underline{C}}}_0& {\underline{\underline{S}}}^{\mathsf{T}}\underline{\underline{\Delta S}}\\ {\underline{\underline{\Delta S}}}^{\mathsf{T}}\underline{\underline{S}}& {\underline{\underline{\Delta S}}}^{\mathsf{T}}\underline{\underline{\Delta S}}\end{array}\right].$$

$${\underline{\underline{S}}}^{\mathsf{T}}\underline{\underline{S}}$$

$$\underline{\underline{S}}\approx \underline{\underline{B}}{\widetilde{\underline{\underline{\Theta }}}}^{\frac{1}{2}}{\widetilde{\underline{\underline{V}}}}^{\mathsf{T}}.$$

$$\underline{\underline{C}}\approx \left[\begin{array}{cc}{\widetilde{\underline{\underline{C}}}}_0& \widetilde{\underline{\underline{V}}}{\widetilde{\underline{\underline{\Theta }}}}^{\frac{1}{2}}{\underline{\underline{B}}}^{\mathsf{T}}\underline{\underline{\Delta S}}\\ {\underline{\underline{\Delta S}}}^{\mathsf{T}}\underline{\underline{B}}{\widetilde{\underline{\underline{\Theta }}}}^{\frac{1}{2}}{\widetilde{\underline{\underline{V}}}}^{\mathsf{T}}& {\underline{\underline{\Delta S}}}^{\mathsf{T}}\underline{\underline{\Delta S}}\end{array}\right]=\underset{{\underline{\underline{V}}}_{\ast }}{\underbrace{\left[\begin{array}{cc}\widetilde{\underline{\underline{V}}}& \underline{\underline{0}}\\ \underline{\underline{0}}& \underline{\underline{I}}\end{array}\right]}}\underset{{\underline{\underline{C}}}_1}{\underbrace{\left[\begin{array}{cc}\widetilde{\underline{\underline{\Theta }}}& {\widetilde{\underline{\underline{\Theta }}}}^{\frac{1}{2}}{\underline{\underline{B}}}^{\mathsf{T}}\underline{\underline{\Delta S}}\\ \text{sym.}& {\underline{\underline{\Delta S}}}^{\mathsf{T}}\underline{\underline{\Delta S}}\end{array}\right]}}\underset{{\underline{\underline{V}}}_{\ast }^{\mathsf{T}}}{\underbrace{\left[\begin{array}{cc}{\widetilde{\underline{\underline{V}}}}^{\mathsf{T}}& \underline{\underline{0}}\\ \underline{\underline{0}}& \underline{\underline{I}}\end{array}\right]}}.$$

$${\underline{\underline{C}}}_1={\underline{\underline{V}}}_1{\underline{\underline{\Theta }}}_1{\underline{\underline{V}}}_1^{\mathsf{T}}.$$

$$\underline{\underline{B}}=\Big{[}\underline{\underline{S}}\quad\underline{\underline{\Delta S}}\Big{]}\,\widetilde{\underline{\underline{V}}}\,\widetilde{\underline{\underline{\Theta}}}^{-\frac{1}{2}}\approx\Big{[}\underline{\underline{B}}\,\widetilde{\underline{\underline{\Theta}}}^{\frac{1}{2}}\,{\widetilde{\underline{\underline{V}}}}^{\sf T}\quad\underline{\underline{\Delta S}}\Big{]}\,\underbrace{\underline{\underline{V}}_{\ast}\,\widetilde{\underline{\underline{V}}}_{1}}_{\widetilde{\underline{\underline{W}}}}\,\widetilde{\underline{\underline{\Theta}}}_{1}^{-\frac{1}{2}}.$$

Full text

\pubvolume

xx \issuenum1 \articlenumber5

\historyReceived: date; Accepted: date; Published: date

\TitleData-Driven Microstructure Property Relations \AuthorJulian Lißner1,†,‡ and Felix Fritzen2,†,‡\orcidB\AuthorNamesFirstname Lastname, Firstname Lastname and Firstname Lastname

\corresCorrespondence: fritzen@mechbau.uni-stuttgart.de; Tel.: +49-711-685 66283 (F.F.) \firstnoteUniversity of Stuttgart, Pfaffenwalding 7, 70569 Stuttgart \secondnoteThese authors contributed equally to this work. \MSC74-04, 74A40, 74E30, 74Q05, 74S30

1 Introduction

In material analysis and design of heterogeneous materials, multiscale modeling can be used for the discovery of microstructured materials with tuned properties for engineering applications. Thereby, it contributes to the improvement of the technical properties, reduces the amount of resources invested into the construction and enhances the reliability of "material modeling"/"the description of material behaviour". However, the discovery of materials with the desired material property, which is characterized by the microstructure of the solid, constitutes a highly challenging inverse problem.

The basis for all multiscale models and simulations is information on the microstructure and on the microscale material behaviour. If at hand, physical experiments can be replaced by – often costly – computations in order to determine the material properties by virtual testing Ghosh et al. (1995); Dhatt et al. (2012); Leuschner and Fritzen (2018). Separation of structural and microstructural length scales can often be assumed. This enables the use of the representative volume element (RVE) Torquato (2013) equipped with the preferable periodic fluctuation boundary conditions (Jiang et al., 2001). The RVE characterizes the highly heterogeneous material using a single frame (or image) and the (analytical or numerical) computation can be conducted on this frame.

The concurrent simulation of the underlying microstructure (e.g., through nested FE simulations, cf., e.g., Feyel (1999); Miehe (2002), or e.g. considering microstructure behaviour in the constitutive laws, e.g. Beyerlein and Tomé (2008)) and of the problem on the structural scale is computationally intractable. In view of the correlation between computational complexity and energy consumption, nested FE simulations should be limited in application. Therefore, efficient methods giving reliable prediction of the material property are an active field of research: POD-driven reduced order models with hyper-reduction (e.g., Ryckelynck (2009); Hernández et al. (2014)), with multiple reduced bases spanning also internal variable Fritzen and Hodapp (2016); Leuschner and Fritzen (2017) and for finite strains (e.g., Yvonnet and He (2007); Kunc and Fritzen (2019)). See also general review articles on the topic such as Kanouté et al. (2009); Matouš et al. (2017).

Supposing that two similar images representing microstructured materials are considered, it is natural to expect similar effective properties in many physically relevant problems such as elasticity, thermal and electric conduction to mention only two applications. The main task, thus, persists in finding low-dimensional parameterizations of the images that capture the relevant information, use these parameterizations to compress the image information into few numbers and build a surrogate model operating only on the reduced representation. A black-box approach exploiting precomputed data for the construction of the surrogate is the use of established machine learning methods which is also the topic of this paper.

As the no free lunch theorem Wolpert and Macready (1997) states, an algorithm can not be arbitrarily fast and arbitrarily accurate at the same time. Hence, there has to be a compromise either in accuracy, computational speed or in versatility. At the cost of generality, i.e. by focusing on subclasses of microstructures, fast and accurate models can be deployed while still allowing for considerable variations of the microstructures. This does not mean that these subclasses must be overly confined: for instance, inclusion volume fractions ranging from 20 up to 80% are considered in this work. Using a limited number of computations performed on relevant microstructure images, machine learned methods can be trained for the subclass under consideration. The sampling of the data, the feature extraction and the training of the ML algorithm constitute the offline phase in which the surrogate model is built. Typically, the evaluation of the surrogate can be realized almost in real-time (at least this is the aspired and ambitious objective), thereby enabling previously infeasible applications in microstructure tayloring, interactive user interfaces and computations on mobile devices.

A currently active research for microstructure property linkage is the material knowledge system (MKS) framework Brough et al. (2017). Many branches thereof exist, all trying to attain low-dimensional microstructure descriptors from the truncation of selected $n$-point correlation functions. For instance, a PCA of the $n$-point correlation functions of the microstructure is performed and the principal scores are used to in a polynomial regression model in order to predict material properties. The MKS is actively researched for different material structures Paulson et al. (2017); Gupta et al. (2015); Kalidindi (2012). For instance, Paulson et al. (2017); Gupta et al. (2015) successfully predict the elastic strain and yield stress for the underlying microstructure using the MKS approach, however they confine their focus on either the topological features of the microstructure or a confined range of allowed volume fractions (0-20%), often held constant in individual studies.

The goal in the present study is to make accurate image based predictions for RVEs spanning large subsets, e.g. in terms of volume fraction, morphological and topological variations, of microstructure materials.

Similarly to key ideas of the MKS approach, a reduced basis is deployed to reduce the dimensionality of the microstructural features contained in the $n$-point correlation functions. With the sheer amount of samples required, conventional methods fail to capture the key features of the all the microstuctures and we propose three novel incremental reduced basis updates to make the computation possible. Using synthetic microstructures/RVE, the costly training of the reduced basis and of the artificial neural network (e.g., Basheer and Hajmeer (2000)) becomes feasible, allowing to build a surrogate model for the image-property linkage. The surrogate accepts binarized image representations of bi-phasic materials as inputs. The outputs constitute the effective heat conductivity of the considered material.

In Section 2 the microstructure classification and the three different incremental snapshot POD procedures used during feature extraction are presented (unsupervised learning). In Section 3 the use of feed forward artificial neural networks for the processing of the extracted features is discussed. Numerical examples are presented in Section 4 including different inclusion morphologies and an investigation of the relaxation of the microstructure subclass confinement of the procedure by using mixed data sets. A Python code illustrating the application of the surrogate is freely available Github.

2 Materials and Methods

2.1 Microstructure Classification

The microstructure is defined by the representative volume element (RVE) Torquato (2013), which is one periodic frame or image characterizing the heterogeneous material under consideration, see Figure 1 for examples of the microstructure and its 2-point spatial correlation function (see below for its definition). Due to their favorable properties regarding the size of the RVE, periodic fluctuation boundary conditions are used for the computations during the offline phase (Jiang et al., 2001).

The n-point spatial correlation functions represent a widely used mathematical framework for microstructural characterization Torquato and Stell (1982); Berryman (1985). Roughly described, the n-point correlation is obtained by placing a polyline consisting of $(n-1)$ nodes defined relative to the first point by vectors \mbox{\boldmathr}_{1},\mbox{\boldmathr}_{2},\dots.

By placing the first point uniformly randomly into the microstructure and computing the mean probability of finding a prescribed sequence of material phases at the nodes of the polyline (including the initial point) denotes the n-point correlation c_{n}(\mbox{\boldmathr}_{1},\mbox{\boldmathr}_{2},\dots,\mbox{\boldmathr}_{n-1};m_{1},m_{2},\dots,m_{n}), where $m_{k}$ is the material label expected to be found at the $k$th node.

For example, the 1-point spatial correlation function, i.e. the probability of finding phase $m$ $(m\in\{a,b,\dots\})$, yields the phase volume fraction $f_{m}$ of phase $m$. In the present study bi-phasic materials are considered. Here $m=a$ corresponds to the matrix material (drawn blue in Figure 1) and $m=b$ to the inclusion phase (drawn yellow in Figure 1). The trivial relation

[TABLE]

holds. The 2-point spatial correlation function (2PCF) c_{2}(\mbox{\boldmathr};a,b) places the vector $r$ in each pixel/voxel $x$ of the RVE and states the probability of starting in the matrix phase $a$ and ending in the inclusion phase $b$. Mathematically we have

[TABLE]

with $\chi^{(m)}$ being the indicator function of phase $m$, $r$ the point offset and $\langle\,\bullet\,\rangle_{x}$ denoting the averaging operator over the RVE. The 2PCF is efficiently computed in Fourier space by making use of the algorithmically sleek fast Fourier transform (FFT) Cooley and Tukey (1965); Frigo and Johnson (1998):

[TABLE]

where $\mathscr{F}$ and $\mathscr{F}^{-1}$ denote the forward and backward FFT, $\overline{\bullet}$ is the complex conjugate and $\odot$ denotes the point-wise multiplication, respectively. For bi-phasic materials the three different two-point functions $c_{2}(\bullet;a,b)$, $c_{2}(\bullet;a,a)$, $c_{2}(\bullet;b,b)$ are related via:

[TABLE]

In view of computational considerations this redundancy can be exploit. Some key characteristics of the non-negative 2PCF are

[TABLE]

In addition to that, a key property of the 2PCF is its invariance with respect to translations of the periodic microstructure. This property is of essential importance when it comes to the comparison of several images under consideration, i.e. during the evaluation of similarities within images.

Examples of c_{2}(\mbox{\boldmathr};b,b) (referred to also as auto-correlation of the inclusion phase) are depicted by the lower set of images in Figure 1. By the metric of vision, the following characteristics can be observed:

• •

the maximum of c_{2}(\mbox{\boldmathr};b,b) occurs at the corners of the domain (corresponding to \mbox{\boldmathr}=\mbox{\boldmath0});

• •

preferred directions of the inclusion placement and/or orientation correspond to laminate-like images (best seen in the third microstructure from the left);

• •

the domain around \mbox{\boldmathr}=\mbox{\boldmath0} partially reflects the average inclusion shape;

• •

some similarities are found, particularly with respect to shape of the 2PCF at the corners and in the center.

These observations hint at the existence of a low-dimensional parametrization of relevant microstructural features. In the following this property is exploited by using a snapshot proper orthogonal decomposition (snapshot POD) in order to capture reoccurent patterns of the 2PCF. By working on the two-point function the afore-mentioned elimination of possible translations of the images is an important feature.

The influence of higher order spatial correlation functions has been investigated in the literature (e.g., Torquato and Stell, 1982; Fast and Kalidindi, 2011). These considerations often yield minor gains relative to the additional computations and the increased dimensionality***For instance, the 3PCF takes to vectors \mbox{\boldmathr}_{1},\mbox{\boldmathr}_{2}\in\varOmega as inputs. Hence, the full 3PCF is basically inaccessible in practice but only after major truncation.. Although it has been demonstrated that the two point function does not suffice to uniquely describe the microstructure in periodic domains Fullwood et al. (2008), there is evidence that the level of microstructural ambiguity for identical 2PCF can be considered low. Therefor, only the n-point correlation functions up to second order are accounted for in the present study.

2.2 Unsupervised Learning via Snapshot Proper Orthogonal Decomposition

The snapshot POD Sirovich (1987) can be used to construct a reduced basis (RB) Liang et al. (2002); Camphouse et al. (2008); Quarteroni et al. (2016) that provides an optimal subspace for approximating a given snapshot matrix $\underline{\underline{S}}\in{\mathbb{R}}^{n\times n_{s}}$. The matrix $\underline{\underline{S}}$ consists of $n_{s}$ individual snapshots $\underline{s}_{i}\in{\mathbb{R}}^{n}$ with the size $n$ being the dimension of the discrete representation of the unreduced field information. In the case of the 2PCF $n$ denotes the total number of pixels within the RVE, i.e. the discrete two-dimensional 2PCF (representing image data) is recast into vector format for further processing ($\,\underline{c}^{0}_{2}(m,m)\in{\mathbb{R}}^{n}\,$). In the present study, the constructed RB is used for information compression, i.e. for the extraction of relevant microstructural features from the image data. The reduced basis $\underline{\underline{B}}\in{\mathbb{R}}^{n\times N}$ retains the $N$ most salient features of the data contained in $\underline{\underline{S}}$ in a few eigenmodes represented by the orthonormal columns of $\underline{\underline{B}}$.

The actual snapshot data contained in $\underline{\underline{S}}$ is thus constructed from the discrete 2-point function data $\underline{s}^{0}_{i}\in{\mathbb{R}}^{n}$ via

[TABLE]

where $\underline{1}\in{\mathbb{R}}^{n}$ is a vector containing ones at all entries. The reduced basis is computed under the premise to minimize the overall relative projection error

[TABLE]

with respect to the Frobenius norm $\|\bullet\|_{F}$. The RB can be constructed with multiple methods, e.g. with the snapshot correlation matrix $\underline{\underline{C}}_{S}$ and its eigenvalue decomposition, which is given by

[TABLE]

The following properties of the sorted eigenvalue decomposition hold:

[TABLE]

and $\delta_{ij}$ denotes the Kronecker delta. The dimension of the reduced basis is determined by the POD threshold, i.e. the truncation criterion is given by

[TABLE]

where $\varepsilon>0$ is a given tolerance denoting the admissible approximation error. Then, the reduced basis is computed via

[TABLE]

after truncation of the eigenvalue- and eigenvector matrices to reduced dimension $N$ represented with $\widetilde{\underline{\underline{\Theta}}}\in{\mathbb{R}}^{N\times N}$ and $\widetilde{\underline{\underline{V}}}\in{\mathbb{R}}^{n\times N}$, respectively. The sorting of the eigenvalues with their corresponding eigenvectors leads to the property that the least reoccurent information given in $\underline{\underline{S}}$ is omitted. Hence, the first eigenmode in $\underline{\underline{B}}$ has the most dominant pattern, the second eigenmode the second most etc.

The properties of the reduced basis computed with the snapshot correlation matrix remain the same as for the singular value decomposition (SVD) introduced below.

The SVD Klema and Laub (1980) of the snapshot matrix is given by

[TABLE]

with the following properties (asserting $n_{s}\geq n$)

[TABLE]

and the sorted non-negative singular values $\sigma_{i}$ such that $\sigma_{1}\geq\sigma_{2}\geq\cdot\cdot\cdot\geq\sigma_{n_{s}}\geq 0$. The criterion for determining the reduced dimension $N$ matching (14) takes the form

[TABLE]

Then, the reduced basis is given by truncation of the columns of $\underline{\underline{U}}\to\widetilde{\underline{\underline{U}}}\in{\mathbb{R}}^{n\times N}$

[TABLE]

More specifically, the left subspace associated with the leading singular values is the RB. Both introduced methods yield the exact same result for the same snapshot matrix $\underline{\underline{S}}$.

2.3 Incremental Generation of the Reduced Basis $\underline{\underline{B}}$

The RB is deployed in order to compress the information contained in $n_{s}$ snapshots into an $N$-dimensional set of eigenmodes stored in the columns of $\underline{\underline{B}}\in{\mathbb{R}}^{n\times N}$, where $N\ll n_{s}$ is asserted. Since the RB is computed with the snapshot matrix alone, the information contained in $\underline{\underline{S}}$ needs to contain data representing the relevant microstructure range, i.e. covering the parameter range used in the generation of the synthetic materials, in order for $\underline{\underline{B}}$ to be representative for the problem under consideration.

In the case of microstructural images containing $n$ pixels, a ludicrous amount of $2^{n}$ states could theoretically be considered when allowing for fully arbitrary microstructures. When limiting attention to certain microstructure classes, then less information is required. Still, thousands of snapshots are usually required, at least. In the following, attention is limited to synthetic materials generated using random sequential adsorption of morphological prototypes with variable size, orientation, aspect ratio, overlap and phase volume fraction. Due to the high variability of such microstructures (see, e.g., Figure 1), a large number of snapshots is required that can usually not be stored in memory, i.e. a monolithic snapshot matrix $\underline{\underline{S}}$ is not available. Although attention is limited to two-dimensional model problems in this study, the problem aggravates considerably for three-dimensional images which imply technical challenges of various sort (storage, processing time, data management, …).

In order to be able to generate a rich RB accounting for largely varying microstructural classes, the incremental basis generation represents a core concept within the present work. It enables the RB generation based on a sequence of input snapshots but without the need to store previously considered data except for the current RB. Three different methods are proposed, two of which rely on approximations of the snapshot correlation matrix $\underline{\underline{C}}_{S}$, and one of which relies on the SVD of an approximate snapshot matrix. The general incremental scheme depicted in Figure 2 remains the same for all the procedures, i.e. the only difference is found during the step labeled adjust.

The algorithm is initialized by a small sized set of initial snapshots of the shifted 2-point correlation function (see Section 2.2). Further, the algorithmic variables $n_{\delta}=0$ and $\underline{\underline{\Delta S}}=\underline{\underline{0}}$ are set. The initial RB is computed classically using either the correlation matrix or the SVD (see previous section for details). After computation of the RB, the snapshots are stored neither in memory nor on a hard drive. The algorithm then takes input snapshots in the order of appearance. For each new snapshot $\underline{s}_{i}$ the relative projection error with respect to the current RB is computed:

[TABLE]

If ${\cal P}_{\delta}$ is greater than the tolerance $\varepsilon>0$ the snapshot is considered as inappropriately represented by the existing RB. Consequently, $\underline{s}_{i}$ is appended to a buffer $\underline{\underline{\Delta S}}$ containing candidates for the next basis enrichment and the counter $n_{\delta}$ is incremented. Once the buffer contains a critical number of $n_{\rm a}$ elements the actual enrichment is triggered and the buffer is emptied thereafter. Thereby the computational overhead is reduced. The three different update procedures are described later on in detail. The procedure is continued until $n_{\rm c}>0$ consecutive snapshots were found to be approximated up to the relative tolerance $\varepsilon$. Then the basis is considered as converged for the microstructure class under consideration.

In the following three methods for the update procedure are described. Formally, the update of an existing basis $\underline{\underline{B}}$ with a block of snapshots contained in the buffer $\underline{\underline{\Delta S}}$ is sought-after. The new basis is required to remain orthonormal.

2.3.1 Method A: Append Eigenmodes to $\underline{\underline{B}}$

A trivial enrichment strategy is given in terms of appending new modes to the existing basis while preserving orthonormality of the basis. Therefore, the projection of $\underline{\underline{\Delta S}}$ onto the existing RB is subtracted in a first step:

[TABLE]

It is readily seen that $\underline{\underline{\Delta\widehat{S}}}$ is orthogonal to $\underline{\underline{B}}$. Then the correlation matrix of the additional data and its eigen-decomposition are computed according to

[TABLE]

Eventually, the enrichment is given through the truncated matrices $\widetilde{\underline{\underline{V}}}$ and $\widetilde{\underline{\underline{\Theta}}}$:

[TABLE]

The new basis is then obtained by

[TABLE]

Method A simply adds modes generated from the projection residual $\underline{\underline{\Delta\widehat{S}}}$ in a decoupled way, i.e. the existing basis is not modified. In order to compute the basis update, only the existing RB $\underline{\underline{B}}$ and the temporarily stored snapshots $\underline{\underline{\Delta S}}$ are required.

Remarks on Method A

1. A.1

The truncation parameter $\epsilon$ must be chosen carefully such that

[TABLE]

In particular, the normalization with respect to the original data prior to projection onto the existing RB must be taken. 2. A.2

By appending orthonormal modes to the existing basis it is a priori guaranteed that the accuracy of previously considered snapshots cannot worsen, i.e. an upper bound for the relative projection error of all snapshots considered until termination of the algorithm is given by the truncation parameter $\epsilon$ and $n_{\rm a}$:

[TABLE]

2.3.2 Method B: Approximate Reconstruction of the Snapshot Correlation Matrix

The goal of this iterative update scheme is the accurate approximation of the new correlation matrix

[TABLE]

Here $\underline{\underline{S}}$ denotes all snapshots considered in the RB so far and $\underline{\underline{\Delta S}}$ contains the candidate snapshots. However, the previously used snapshots formally written as $\underline{\underline{S}}$ are no longer available since they can not be stored due to storage limitations. Using the previously computed matrices $\underline{\underline{B}},\widetilde{\underline{\underline{V}}},\widetilde{\underline{\underline{\Theta}}}$ the following approximations are available:

[TABLE]

where the accuracy of the approximation is governed by the truncation threshold $\delta_{N}$. Using these approximations and using the property of the eigenvalue decomposition, the snapshot matrix $\underline{\underline{S}}$ can be approximated by

[TABLE]

Note that $\underline{\underline{B}}\in{\mathbb{R}}^{n\times N}$ is stored anyway, $\widetilde{\underline{\underline{\Theta}}}\in{\mathbb{R}}^{N\times N}$ is diagonal and $\widetilde{\underline{\underline{V}}}\in{\mathbb{R}}^{n\times N}$ is of manageable size. The snapshot correlation matrix $\underline{\underline{C}}$ that considers the additional snapshots can be approximated as

[TABLE]

In order to compute the updated basis, the inexpensive eigenvalue decomposition of $\underline{\underline{C}}_{1}\in{\mathbb{R}}^{(N+n_{\rm a})\times(N+n_{\rm a})}$ is computed

[TABLE]

Analogously to the previous RB computation in equation (15), the adjusted and enriched basis is computed by

[TABLE]

To update the RB the truncated eigenvector matrix ($\underline{\underline{B}}\,,\widetilde{\underline{\underline{V}}}\leftarrow\widetilde{\underline{\underline{W}}}\in{\mathbb{R}}^{n\times N}$) need to be stored as well as the diagonal eigenvalue matrix $\widetilde{\underline{\underline{\Theta}}}$.

Remarks on Method B

1. B.1

The existing RB is not preserved but it is updated using the newly available information. Thereby, the accuracy of the RB for the approximation of the previous snapshots is not guaranteed a priori. However, numerical experiments have shown no increase in the approximation errors of previously well-approximated snapshots. 2. B.2

In contrast to Method A the dimension of the RB can remain constant, i.e. a mere adjustment of existing modes is possible. The average number of added modes per enrichment is well below that of Method A. 3. B.3

The additional storage requirements are tolerable and the additional computations are of low algorithmic complexity.

2.3.3 Method C: Incremental SVD

Method C is closely related to Method B. However, instead of building on the use of the correlation matrix, it relies on the use of an updated SVD, i.e. an approximate truncated SVD is sought-after:

[TABLE]

Since the original snapshot matrix $\underline{\underline{S}}$ can not be stored, only an approximation of the actual truncated SVD in (33) can be computed. Methods to compute an incremental SVD are introduced in Gu and Eisenstat (1994); Fareed et al. (2018), with the latter referring to Brand’s incremental algorithm Horn and Johnson (1985) which is used in the present study with minor modifications. With the previously computed basis $\underline{\underline{B}}$ at hand, the approximation of $\underline{\underline{S}}$ is known

[TABLE]

Introducing the projection residual $\underline{\underline{\Delta\widehat{S}}}$ of the enrichment snapshots $\underline{\underline{\Delta S}}$ together with its SVD

[TABLE]

and using (33), (34) and (35) the full snapshot matrix can be approximated after some algebra by

[TABLE]

Considering the SVD of $\underline{\underline{\Gamma}}=\underline{\underline{U}}_{\Gamma}\,\underline{\underline{\Sigma}}_{\Gamma}\,{\underline{\underline{W}}}^{\sf T}_{\Gamma}\in{\mathbb{R}}^{(N+n_{\rm a})\times(N+n_{\rm a})}$, which is inexpensive due to the sparsity of $\underline{\underline{\Gamma}}$, approximation (36) is further rewritten as

[TABLE]

It is easily shown that the matrices $\underline{\underline{U}}_{\ast}$ and $\underline{\underline{W}}_{\ast}$ are column-orthogonal and that $\underline{\underline{\Sigma}}_{\ast}$ is diagonal and non-negative. Therefore, the three matrices constitute an approximate SVD of the enlarged snapshot matrix at low computational expense. This implies the following updates after each enrichment step

[TABLE]

after truncation of $\underline{\underline{B}}$, where the truncation criteria needs to ensure that $\underline{\underline{B}}$ does not decrease in size. To compute the enrichment of the RB, $\underline{\underline{B}}\in{\mathbb{R}}^{n\times N}$ and the sparse singular values $\underline{\underline{\Sigma}}\in{\mathbb{R}}^{N\times N}$ after truncation need to be stored.

Remarks on Method C

1. C.1

As highlighted for Method B in remark B.1, the existing RB is not preserved but adjusted by considering the newly added information. This ensures the accuracy of the approximation error of previously well-approximated snapshots, but a priori guarantees regarding the subset approximation accuracies cannot be made. 2. C.2

In contrast to Method A the dimension of the RB can remain constant, i.e. a mere adjustment of existing modes is possible. The average number of added modes per enrichment is well below that of Method A. 3. C.3

Each update step in equation (38) is computed separately and, consequently, storing $\underline{\underline{W}}$ is not required since only the RB $\underline{\underline{B}}$ is of interest. 4. C.4

The diagonal matrix $\underline{\underline{\Sigma}}$ has low storage requirements corresponding to that of a vector in ${\mathbb{R}}^{N}$.

3 Supervised learning using Feed Forward Neural Network

During the supervised learning phase, the machine is provided with data sets consisting of inputs and the related outputs. Hence, the supervised learning phase tries to learn a function relating inputs to outputs without or with limited prior knowledge of the structure of the unknown mapping. Artificial neural networks (ANN) are a powerful machine learning tool which have gained wide popularity in the recent decades due to the surge in computational power Widrow and Lehr (1990); Basheer and Hajmeer (2000).

The functionality of the ANN is inspired by the (human) brain, propagating a signal (input) through multiple neurons where it is lastly transformed into an action (output). Various types of neural networks have been invented, e.g. feed forward, recurrent, convolutional, being applicable to almost any field of interest Kimoto et al. (1990); Sundermeyer et al. (2012); Simonyan and Zisserman (2014); Angermueller et al. (2016).

In the present study a regression model from the input, i.e. the feature vector $\underline{\xi}$ which is derived with the converged basis $\underline{\underline{B}}$, to the output, i.e. the effective heat conduction tensor $\bar{\underline{\underline{\kappa}}}$, is deployed with a dense feed forward ANN.

In a dense feed forward ANN (Figure 3) a signal is propagated through the hidden layers where every output of the previous layer $\underline{a}^{l-1}$ affects the activation $\underline{z}^{l}$ of the current layer $l\,(l=1,\dots,L+1)$. The activation of each layer gets wrapped into an activation function $f$ where the output of each neuron in the layers is computed, i.e. $\underline{a}^{l}=f(\underline{z}^{l})$. Note that matrix/vector notation is used, where each entry in the vectors denotes one neuron in the respective layer.

The basic learning algorithm/optimizer usually employed for a feed forward ANN is the back propagation algorithm Hecht-Nielsen (1992) and modifications thereof. The learning of the network consists in the numerical identification of randomly initialized and unknown weights $\underline{\underline{W}}^{l}$ and biases $\underline{b}^{l}$ that minimizes a given cost function. The latter gives an indication of the quality of the ANN prediction. For instance, the gradient back propagation computes gradients of the cost function with respect to the weights and biases, in order to compute suitable corrections for these parameters.

The learning itself is an iterative procedure in which the training data is cycled multiple times through the ANN (one run called an epoch). In each epoch the internal parameters are updated with the aim of improving the mapping from the input to the output. The optimization problem itself is (usually) high-dimensional. In most situations it is not well-posed and local minima and maxima can hinder convergence. Therefor, multiple random instantations of the network parameters are usually required to assure that a good set of parameters is found, even if the network layout remains unaltered.

Notably, the training requires a substantial input data set as inputs. It is important to note that the training of the ANN usually results in a parameter set that is able to approximate the training data with high accuracy under the given meta-parameters describing the network architecture (number of layers, number of neurons per layer, type of activation function). However, the approximation quality of the ANN may be different for other queries. Thus, it is important to validate the generality of the discovered mapping for the underlying problem setting. Therefor an additional validation data set is introduced, where only the evaluation of the cost function is tracked during the training. Generally when overfitting†††Overfitting relates to the fact that a subset of the data is nicely matched but small variations in the inputs can lead to substantial loss in accuracy, similar to oscillating higher-order polynomial interpolation functions. occurs, the errors for the validation set increase whereas the errors of the training set decrease. The training should be halted if such a scenario is detected.

Since the choice of activation function as well as the number of hidden layers and the number of neurons within the individual layers are arbitrary (describing the ANN architecture), these meta-parameters should be tailored specifically for the desired mapping. Finding the best neural network architecture is not straight-forward and usually relies on intuition, experience and many numerical experiments. In order to find a well suited ANN, various (random) realizations of each tested ANN architecture need to be computed, before a decision regarding the optimal layout can be made.

In the present study the ANN training is performed using Tensorflow in Python Abadi et al. (2016). Tensorflow is an open source project by the Google team, providing highly efficient algorithms for ANN implementation. The Adam Kingma and Ba (2014) optimizer, which is a modification of the gradient back propagation, has been deployed for the learning.

4 Results

4.1 Generation of Synthetic Microstructures

All of the used synthetic microstructures have been generated by a random sequential adsortion algorithm with some examples shown in Figure 1. Two morphological prototypes were used: spheres and rectangles. The parameters used to instantiate the generation of a new microstructure were modeled as uniformly distributed variables:

1. M.1

the phase volume fraction $f_{b}$ of the inclusions (0.2-0.8); 2. M.2

the size of each inclusion (0.0–1.0); 3. M.3

for rectangles: the orientation (0–$\pi$) and the aspect ration (1.0–10.0); 4. M.4

the admissible relative overlap $\varrho$ for each inclusion (0.0–1.0).

For $\varrho=0$ and the spherical inclusion, a boolean model of hard spheres is obtained. Setting $\varrho=1$ induces a boolean model without placement restrictions, i.e. new inclusions can be placed independent of the existing ones. The generated microstructures were stored as images with resolution 400$\times$400. After the generation of the RVE, the 2-point spatial correlation function was computed for the RVE. This was then shifted, see Section 2.2, and used as a snapshot $\underline{s}_{i}$ for the identification of the reduced basis.

Additionally, a smaller random set of RVEs used for the supervised learning phase was simulated using the recent Fourier-based solver FANS Leuschner and Fritzen (2018) in order to compute the effective heat conduction tensor $\bar{\underline{\underline{\kappa}}}$. The heat conductivity of the matrix and of the inclusion phase are prescribed as

[TABLE]

Here $R>0$ denotes the material contrast. In the present study, $R=5$ was considered, i.e. the matrix of the microstructure has a five times higher conductivity than the inclusions. These values can be seen as typical values for metal ceramic composites (Figure 4).

An inverse phase contrast has exemplarily been studied, i.e. inclusions with $\kappa_{b}=1$ and $\kappa_{a}=\frac{\kappa_{b}}{5}$ (corresponding to $R=\tfrac{1}{5}$) has also been investigated. Qualitatively, the results for the inverse phase contrast did not show any new findings or qualitative differences. Therefor, the following results focus on $R=5$, corresponding to rather insulating inclusions.

The symmetric tensor $\bar{\underline{\underline{\kappa}}}$ can be represented as a three-dimensional vector $\bar{\underline{\kappa}}$ using the normalized Voigt notation

[TABLE]

For the supervised learning of the ANNs (see Section 3), multiple files each containing 1500 data sets for different inclusion morphologies were generated (circle only; rectangle only; mixed). Each data set contains the image of the microstructure, the respective autocorrelation of the inclusion phase $c_{2}(\bullet;b,b)$ and the effective heat conductivity $\bar{\underline{\kappa}}_{\mathrm{V}}$.

4.2 Unsupervised Learning

First, the reduced basis is identified using the iterative procedure presented in Section 2.3. All three proposed methods were considered and for each of these, three different sets of microstructures were considered: The first set of microstructures consisted of RVEs with only circular inclusions, the second set consisted of RVEs with only rectangular inclusions, and the third set was divided into equal parts, each part consisting of RVEs with either circular or rectangular inclusions‡‡‡(i.e. each structures contained exclusively one of the two morphological prototypes and the same number of realizations for each prototype was enforced), respectively. Each type of microstructure was processed using each of the three incremental RB schemes introduced in Section 2.3. Hence, a total of nine different trainings were conducted, each using different randomly generated snapshots.

For the iterative enrichment process, the initial RB was computed from 200 snapshots $\underline{\underline{S}}_{0}$. Thereafter, snapshots were randomly generated and processed by the enrichment algorithm sketched in Figure 2. The number of snapshots per enrichment step has been set to $n_{\rm a}=75$ and the number of consecutive snapshots with $\mathcal{P}_{\delta}<\varepsilon$ , used to indicate convergence, has been set to $n_{\rm c}=100$. The relative projection tolerance $\varepsilon=0.025$ was chosen. Note that this corresponds to the maximum value of the mean relative $\|\cdot\|_{L^{2}}$-error that is considered exact for the shifted snapshots. The actual accuracy in the reproduction of the 2PCF c_{2}(\mbox{\boldmathr};b,b) is significantly lower than this (results are given in Figure 7 below).

Key attributes for each of the nine trainings are provided in Table 1. There is an obvious discrepancy between Method A and the remaining methods in basically all outputs. While Method A claims the lowest computing times, it yields approximately twice the number of modes. However, the number of snapshots needed is substantially lower which can be relevant if the generation of the synthetic microstructures is computationally involved.

Note that methods B and C yield similar results, although for the rectangular and circular training method C needed significantly more snapshots, method B needed significantly more snapshots for the mixed training. The outliers in the number of snapshots needed are due to the randomness of the materials and the chosen convergence criterion. The resulting basis size of methods B and C indicate very similar results from these methods.

Also note that the calculation of the relative projection error $\mathcal{P}_{\delta}$ grows linearly with the dimension of the RB, i.e. the faster offline time of method A can quickly be compensated by the costly online procedure induced by the high dimension of the RB in comparison to the competing techniques.

To compare the accuracy of the resulting basis as well as during the training, the relative projection error $\mathcal{P}_{\delta}$ of the snapshots used for the original basis construction $\underline{\underline{S}}_{0}$ are plotted in Figure 5.

Since each training conducted a different amount of incremental updates, the abscissa has been rescaled such that the relative progress of the basis generation is shown vs. the number of enrichments divided by $n_{\mathrm{it}}$ (Figure 5 and Table 1).

Method B and C yield again very similar results whereas method A achieves a lower projection error on convergence, but at the expense of a considerably improved dimension of the RB. Since there seems to be an obvious correlation between resulting accuracy and the final basis size for the initial snapshots $\underline{\underline{S}}_{0}$ (see Figure 5, Table 1), the general quality for arbitrary stochastic inputs must be investigated. In order to quantify the quality of the RB, the accuracy can be expressed in terms of the relative projection error of approximating additional, newly generated snapshot data $\underline{\underline{S}}$ as a function of the method (A, B, C) and the number of modes $N\geq 1$ via

[TABLE]

in Matlab notation.

This measure captures to what extend the first $N$ basis functions represent the 2PCF of the underlying microstructure class. In the current work sets of 1500 newly generated snapshots assure an unbiased validation, i.e. the data was used in neither of the three training procedures. The results are stated in Figure 6. Again, methods B and C yield similar results, achieving lower projection errors with fewer eigenmodes compared to Method A, i.e. the basis produced by method A cannot catch up with its two competitors. On a side note, the rectangular inclusions apparently lead to significantly richer microstructure information which can be seen by direct comparison of the left to the middle plot in Figure 6. For methods B and C and for circular inclusions the relative error of 5% is reached for approximately 15 modes while rectangular inclusions require more than 60 modes to attain a similar accuracy. This is supported also by the rightmost plot determined from a sort of blend of the two microstructural types.

Since all of the previous error measures are given on the shiftedsnapshot according to equation (10), the true relative projection error on the unshifted snapshot is also investigated as a function of the basis size. It describes the relative accuracy of the approximation of the 2PCF c_{2}(\mbox{\boldmathr};b,b) as a function of the basis size. The errors in the shifted data (Figure 7, left) and the corresponding reconstructed 2PCF (Figure 7, right) for five individual randomly selected snapshots show that the actual relative error in the 2PCF reconstruction is below 5% for 10 reduced coefficients even for the challenging rectangular inclusion morphology, while the error in the shifted and shifted snapshots is on the order of 50%. This highlights the statement made earlier regarding the choice of $\varepsilon$ which is not directly the accepted mean error in the 2PCF, but only after application of the shift. The high discrepancy in the two relative projection errors is due to the fact that the shifted snapshots fluctuate closely around 0, i.e. the homogeneous part of the 2PCF is obviously of high relevance.

The development, i.e. the stabilization of the mode shapes over the enrichment steps, of a few selected eigenmodes is shown in Figure 8 using RVEs with circular inclusions for training of method C. Similar results are expected for method B, whereas for method A the eigenmodes would remain unconditionally unchanged over the enrichment steps, i.e. a pure enlargement of the basis takes place. The faster stabilization of the leading eigenmodes indicates a quick stabilization of the lower order statistics of the microstructure ensemble, while the tracking of higher order fluctuations is more involved.

4.3 Supervised Learning

After the training of the RB, the input for the neural network, the feature vector $\underline{\xi}$ was derived using the 1- and 2-point spatial correlation functions of the $i$th RVE as

[TABLE]

The size of the feature vector is determined by the amount of reduced coefficients $1\leq h\leq N$, i.e. the snapshot is projected onto the leading $h$ eigenmodes of $\underline{\underline{B}}$.

Since the inputs and outputs have a highly varying magnitude, they need to be shifted such that they are equally representative. Therefore, each entry of the feature vector is separately shifted and scaled such that its distribution of all samples has zero mean and a standard deviation of one. The output is shifted combinedly such that the mean of $\underline{\bar{\kappa}}_{\rm V}$ is $\underline{0}$. The transformed inputs and outputs are then given to the ANN for the training phase. Thus, the outputs of the ANN need to undergo an inverse scaling in order to yield the sought-after vector representation of the heat conduction tensor. These shifts and scalings need to be extracted from the available training data. Hence, every data set used for training purposes has its own parameters.

The training for the neural network has been conducted for all of the three microstructure classes, i.e. using only RVEs with circular inclusions, only RVEs with rectangular inclusions and lastly using RVEs with either circular or rectangular inclusions (the split was up to 60:40, which was randomly assigned before the training). In order to derive the feature vector, the converged basis of method C has been used. Note that depending on the training set, only circular/rectangular or mixed inclusions contributed to the RB.

The training of the ANN was conducted with an early stop algorithm: up to 10000 epochs were considered and the best–not necessarily last–parameter set has been saved. The decision on the best ANN was taken on the basis of the cost function of the validation set (which came out of the same parameter range as the training set). Each data set consisted of 1500 samples. These were shuffled randomly and split into the training set ($n_{t}=1000$) and the validation set ($n_{v}=500$) before each ANN training.

In the following, the error measurements used and the term of unbiased testing refers to the prediction of 1500 unseen data points for each of the three microstructure classes named test sets. In order to find a good overall ANN, the network architecture has been intensely studied: the accuracy of the prediction after the training has been evaluated with various sizes of the feature vector, different network layouts and for different activation functions (Figure 9). The depicted error measure (Figure 9) is the mean/max absolute derivation, i.e. the absolute error (corresponding to the Euclidean norm of the vector-valued error) for the prediction and the actual heat conductivity of the microstructure of $\underline{\overline{\kappa}}_{\rm V}$ over each test set.

The conductivity $\overline{\kappa}_{12}$ fluctuates mildly around zero for all inputs. In order to accurately capture this fluctuation, only the specific training and RB dimensions of four or higher ($h\geq 4$) are required cf. Figure 10. However, the error of $\overline{\kappa}_{12}$ for other inclusion morphologies than that of the training increase with $h$, albeit the values can be considered small in comparison to the $\overline{\kappa}_{11}$ and $\overline{\kappa}_{22}$ errors.

Since the trend in Figure 9 of the rectangular and mixed training looks promising, a further study using more reduced coefficients has been conducted, however the results for retangular training showed a stronger tendency towards overfitting. Rectangular training with $17\leq h\leq 25$ had very similar results to circular training with $12\leq h\leq 17$ reduced coefficients. For the mixed training the overall results worsened with a higher amount of reduced coefficients and the ANN did not seem to find any mapping generalizing the property linkage for both types of deployed RVE with the given training data.

The increased overfitting with respect to the training microstructure class for increased dimension of the feature vector can be explained by the rather limited number of 1000 input samples. For instance, considering a six-dimensional feature vector induces that for the rather limited number of ten samples per independent direction a total of 106 data points would be needed. The dilemma is that each of the computations is expensive, particularly when considering three-dimensional simulations, i.e. millions of samples can not be realized in practice. Therefore, the number of samples must be rather low which could be a limiting factor in view of the number of features that can be accounted for.

Overall, a correlating trend with the accuracy for the RB (Figure 6) and the behaviour of the ANN training could be seen (Figure 9). Investigating the circular training, with more specific information about the 2PCF available, a slightly better mapping is found for RVEs with circular inclusions, whereas the prediction for the other two microstructure classes slightly worsen. Similarly with the rectangular training, more reduced coefficients increase the fit on the RVE with rectangular inclusions, though more reduced coefficients are required to deliver the same accuracies as the circular training, as more number of eigenmodes (of the basis) are required to yield approximately the same relative projection error.

When training with both types of RVEs (i.e. for the mixed input set), the training seemed more random than the others. Some resulting ANNs (which are not shown) had more of the property of a circular training, whereas some were more like the rectangular training. Although the ANN was trained with general RVEs, it overall failed to give good predictions for every type of RVE when using more than two or three reduced coefficients. After training with the mixed set, the prediction of $\bar{\kappa}_{11}$ and $\bar{\kappa}_{22}$ was generally better for RVEs with rectangular inclusions, whereas $\bar{\kappa}_{12}$ was better for RVEs with circular inclusions. During the training, the ANN found most likely some strange local minimum which fitted the training data quite well, however it was not a general mapping for all the microstructure classes, which hints at a too low number of input samples, as discussed earlier.

The ANN architecture did not seem to have a big impact on the quality of the prediction, there have been ANNs with a single hidden layer and six hidden neurons which delivered results comparable to an ANN comprising six hidden layers with more than 100 hidden neurons.

The used activation functions were the sigmoid, relu, tanh and softplus, where only some combinations delivered poor results. Not a clear trend of ANN architecture and quality of prediction could be seen and, consequently, the best ANN were randomly found based on the lowest error on the test set.

The prediction accuracies for each test set of three differently trained ANNs, which have been deemed the best, is given in Figure 11. As to be expected, the lowest errors are achieved on the diagonal, i.e. training set $\widehat{=}$ validation set.

The training and architecture of the best ANNs in Figure 11 had the following properties:

[TABLE]

For an easier readability, the percentage mean and max errors for each ANN training and prediction are given in Table 2, again using the same three ANNs which have been explicitly shown. Note that since the values of $\bar{\kappa}_{12}$ vary closely around 0 (Figure 4), relative errors are not sensible for the quantity of interest.

A GUI code is provided in Github, where the user can choose between the three proposed ANN, the input for the prediction is a $400\times 400$ image in matrix format written in a text file or a TIFF image and the output is the prediction for the heat conduction tensor as described above. In order to compile the code, Python3 with Tensorflow is required. Some exemplary RVE with their respective heat conductivity are uploaded in a subfolder.

5 Conclusion

5.1 Summary and Concluding Remarks

The computational homogenization of highly heterogeneous microstructures is a challenging procedure with massive computational requirements. In the present study a method to efficiently and accurately predict the heat conductivity for any RVE with the image and no further information is proposed. Key ideas of the Materials Knowledge System (MKS) Fast and Kalidindi (2011); Kalidindi (2012) have been adopted in the sense that a subset of the POD compressed 2-point correlation function is used to identify a low-dimensional microstructure description. In contrast to Fast and Kalidindi (2011) the 2PCF is not truncated to a small neighborhood, but the full field information is considered. Similar to other works related to the MKS Brough et al. (2017), a truncated PCA of the 2-point information is used to extract microstructural key features.

However, the classical truncated PCA used, e.g., in Brough et al. (2017) is not applicable to the considered rich class of microstructures due to the high number of needed samples and the related unmanageable computational resources. Therefor, our proposal is founded on a novel incremental procedure for the generation of the RB of the 2PCF. Similar techniques have not been considered in the literature to the best of the authors’ knowledge. The shifting of the images of 2PCF before entering the POD is another feature that can help in reducing the impact of the inclusion volume fraction, i.e. the shifted function has zero mean.

Other than in Fast and Kalidindi (2011) no higher-order statistics are used. This is by purpose as the selection of the relevant entries of the higher order PCF is ambiguous and a challenge in itself. Instead, the present study focus on the variability of the input images in terms phase volume fractions in a broad range (20-80%) alongside topological variations (impenetrable, partial overlap, unrestricted placement) and different morphologies (circles and rectangles). Generally speaking, a much higher microstructural variation is accounted for, than in previous studies. Therefore, the current study also investigates how the proposed technique and similar MKS related approach can possibly generalize towards truly arbitrary input images (e.g. millions of snapshots) in order to built a super-database.

In order to cope with the variability of the 2PCF, the truncated PCA or snapshot POD effected during unsupervised learning phase is replaced by novel incremental procedures for the construction of a small-sized reduced microstructure parameterization. Three increment POD methods are proposed and their results are compared regarding the computational effort, the projection accuracy of the snapshots and the quality of the basis in view of capturing random inputs.

The learned reduced bases are used to extract a low-dimensional feature vector which denotes the input of a fully connected feed forward artificial neural network. The ANN is used to predict the homogenized heat conductivity of the material defined by the microstructure. The mean relative error of the surrogate is lower than 2% for the majority of the considered test data. This is remarkable in view of the phase contrast $R=5$ and the particle volume fractions ranging from 0.2-0.8, as well as morphological and topological variations. Further, an immense speedup in computing time is achieved by the surrogate over FE or FFT simulations.

Importantly, the presented methodology can immediately be adopted to different physical settings such as thermo-elastic properties, fluid permeability, dielectricity constants etc. The same holds for three-dimensional problems. However, the limited number of samples in 3d could be problematic as more features are likely required to attain a sufficiently accurate RB.

5.2 Discussion and Outlook

A weakness of the current approach remains the computational complexity of the method: Although the feature vector is rather low-dimensional, it requires the evaluation of the 2PCF using the FFT which is of complexity ${\cal O}(n\,\text{log}(n))$ where $n$ is the number of pixels/voxels in the image. In order to extract the reduced coefficient vector from the 2PCF, the latter must be projected onto the RB. This operation scales with ${\cal O}(n\,N)$. Both operations are at least linear to the number of pixels/voxels of the image which can be critical, especially in three-dimensional settings. Therefor, future investigations aiming at a reduced effort for gathering the features are required in our opinion.

Another extension of the current scheme could account for variable phase contrast $R$ which was fixed as $R=5$ in this work. Thereby, the dimension of the feature vector will be incremented which can add to the already existing data scarcity dilemma observed when considering a sufficient number of reduced coefficients: the number of input samples for the supervised learning should grow exponentially with the dimension of the feature vector. In the authors opinion this dependence is the most pronounced short-coming of the method and future studies should focus on limiting the number of required input samples in order to fight the curse of dimensionality: more reduced coefficients require an exponential growth in the available data, making the offline procedure unaffordable, today.

Advantages of the current scheme comprise the independence of the underlying simulation scheme. This does allow for heterogeneous simulation environments, the use of commercial software, multi-fidelity input data and blended sources of information (e.g. in silico data supported by experimental results).

\authorcontributions

Conceptualization, J.L. and F.F.; Data curation, J.L.; Formal analysis, F.F.; Funding acquisition, F.F.; Investigation, J.L. and F.F.; Methodology, J.L. and F.F.; Project administration, F.F.; Resources, F.F.; Software, J.L.; Supervision, F.F.; Validation, J.L.; Visualization, J.L.; Writing – original draft, J.L.; Writing – review & editing, J.L..

\funding

This research was funded by Deutsche Forschungsgemeinschaft (DFG) within the Emmy-Noether programm under grant DFG-FR2702/6 (contributions of F.F.).

Acknowledgements.

Support from Mauricio Fernández on the implementation and layout for the training of the Artificial Neural Networks using Google’s TensorFlow is highly appreciated. Stimulating discussions within the Cluster of Excellence SimTech (DFG EXC2075) on machine learning and reduced basis methods are highly acknowledged. \conflictsofinterestThe authors declare no conflict of interest. \reftitleReferences