Julien Legendre

Pierre-Olivier Chapuis

Operating conditions and thermodynamic bounds of dual radiative heat engines

A dual radiative heat engine is a device made of two facing optoelectronic components (diodes) and capable of generating electrical power from heat. It can operate in three different regimes depending on the component biases, namely in thermoradiative-negative electroluminescent (TRNEL), thermoradiative-photovoltaic (TRPV) or thermophotonic (TPX) regimes. The use of dual engines gives access to operating conditions which are unachievable by single radiative engines such as thermophotovoltaic systems: at the radiative limit, TRNEL devices systematically reach the Carnot efficiency, while TPX devices can achieve large power outputs by means of electroluminescent enhancement. Upper bounds of the maximum power output and related efficiency achieved by dual engines are derived analytically, and compared to usual bounds. Spectral filtering and nonradiative recombinations are also briefly considered. This work provides common framework and guidelines for the study of radiative engines, which represent a promising solution for reliable and scalable energy conversion.

The conversion of heat into electrical power by optoelectronic means has gathered sizeable attention in the past decades. The two most popular solutions are photovoltaics (PV) for solar application, and thermophotovoltaics (TPV) [START_REF] Burger | Present Efficiencies and Future Opportunities in Thermophotovoltaics[END_REF][START_REF] Lapotin | Thermophotovoltaic efficiency of 40%[END_REF][START_REF] Tervo | Efficient and scalable GaInAs thermophotovoltaic devices[END_REF][START_REF] Giteau | Thermodynamic performance bounds for radiative heat engines[END_REF]. In the latter case, the radiation comes from a hot emitter maintained at a high temperature by the input heat. This gives access to a broad range of applications, for instance in latent heat TPV batteries [START_REF] Datas | Latent heat thermophotovoltaic batteries[END_REF]. Apart from PV and TPV, other optoelectronic systems are able to convert heat to electricity. One is the thermoradiative (TR) cell: as opposed to the PV cell, it is able to generate electrical power by emitting negative electroluminescent radiation to the cold surroundings [START_REF] Strandberg | Theoretical efficiency limits for thermoradiative energy conversion[END_REF][START_REF] Pusch | Fundamental Efficiency Bounds for the Conversion of a Radiative Heat Engine's Own Emission into Work[END_REF] (e.g. the night sky, outer space) or towards a cold absorber [START_REF] Santhanam | Thermal-to-electrical energy conversion by diodes under negative illumination[END_REF]. TR cells, along with PV and TPV cells, are single radiative heat engines: they are heat engines in which one active component produces electrical power either by emitting or absorbing radiation [START_REF] Tervo | Nearfield radiative thermoelectric energy converters: a review[END_REF]. Their typical electrical characteristic is provided in Section I of Supp. Mat. [10]. Radiative engines are part of the larger group of solid-state heat engines [START_REF] Melnick | From thermoelectricity to phonoelectricity[END_REF], along with thermoelectric [START_REF] Shi | Advanced Thermoelectric Design: From Materials and Structures to Devices[END_REF] and thermionic generators [START_REF] Abdul Khalid | Review on Thermionic Energy Converters[END_REF], and thus compete with such systems in the field of reliable and scalable energy conversion.

It is also possible to couple two different optoelectronic components into a dual radiative heat engine (see Fig. 1). One such dual engine is the thermoradiative-photovoltaic (TRPV) device [START_REF] Liao | Thermoradiative-Photovoltaic Cells[END_REF][START_REF] Tervo | Solar Thermoradiative-Photovoltaic Energy Conversion[END_REF], in which a hot TR cell is coupled to a cold PV cell: both components are then able to generate electrical power, although production by the TR cell limits its emission due to negative electroluminescence, and therefore reduces the PV cell production. Nonetheless, having two components provides better control of the operating conditions. Thermophotonics (TPX) [16][START_REF] Mcsherry | Effects of narrowband transport on near-field and far-field thermophotonic conversion[END_REF][START_REF] Sadi | Prospects and requirements for thermophotonic waste heat energy harvesting[END_REF] has also gathered interest recently. In such a device, the hot emitter is a light-emitting diode (LED). While an LED consumes electrical power, it is

Heat source (T h ) Hot optoelectronic component (µ h , T h ) Cold optoelectronic component (µ c , T c ) Heat sink (T c ) (q h , Ṅh ) (q c , Ṅc ) q source q sink P P h P c
FIG. 1: Representation of the dual radiative engine (here in a thermophotonic configuration).

able to emit electroluminescent radiation towards the PV cell with an above-unity wall-plug efficiency, enabling a significant increase in power output. Recently, TPX devices have mostly been studied in near-field operation, as near-field radiation further increases the power output [START_REF] Zhao | High-performance near-field thermophotovoltaics for waste heat recovery[END_REF][START_REF] Legendre | GaAs-based near-field thermophotonic devices: Approaching the idealized case with one-dimensional PN junctions[END_REF] and limits the impact of non-radiative losses [START_REF] Legendre | Overcoming nonradiative losses with AlGaAs PIN junctions for near-field thermophotonic energy harvesting[END_REF].

Although the aforementioned systems are all radiative heat engines and share the same core physical laws, they are studied independently of one another in the literature. This makes the comparison of their performance less straightforward. Therefore, in the current Letter, we provide a complete overview of the performance of dual radiative heat engines achieved for various bandgaps, and highlight the interests of each operating regime. In particular, we focus on the maximum achievable power and efficiency, and derive analytical upper bounds for these two quantities.

In the following, we assume that the dual engine operates in the far field. To maximise the achievable power output, the emission of radiation is assumed to follow the modified Planck law. Furthermore, the bandgaps of the two components are assumed to be equal, and the radiation exchanged below the bandgap is neglected to obtain an upper bound of the efficiency. Thus, the emitted photon flux density Ṅ and the related heat flux density q can be expressed as

Ṅi = 1 4π 2 c 2 h3 ∞ Eg E 2 exp E-µi kBTi -1 dE, ( 1a
)
q i = 1 4π 2 c 2 h3 ∞ Eg E 3 exp E-µi kBTi -1 dE, ( 1b
)
where i relates to the emitting body (either "h" or "c") with bandgap energy E g and temperature T i . µ i is the chemical potential of the emitted radiation, which must remain strictly smaller than E g and is assumed to be related to the voltage U i applied to the component through µ i = eU i [START_REF] Callahan | Coupled Charge and Radiation Transport Processes in Thermophotovoltaic and Thermoradiative Cells[END_REF][START_REF] Legendre | Theoretical and numerical analysis of near-field thermophotonic energy harvesters[END_REF], e being the elementary charge. The above integrals can be expressed analytically using polylogarithms [START_REF] Green | Analytical treatment of Trivich-Flinn and Shockley-Queisser photovoltaic efficiency limits using polylogarithms[END_REF], leading to

Ṅi = (k B T i ) 3 4π 2 c 2 h3 e 2 g,i Li 1 (e -xi ) + 2e g,i Li 2 (e -xi ) + 2 Li 3 (e -xi ) , (2a)
q i = (k B T i ) 4 4π 2 c 2 h3 e 3 g,i Li 1 (e -xi ) + 3e 2 g,i Li 2 (e -xi ) + 6e g,i Li 3 (e -xi ) + 6 Li 4 (e -xi ) , (2b)
where Li n is the n-th order polylogarithm, e g,i = E g /k B T i and x i = (E g -µ i )/k B T i . Finally, to obtain an upper bound of the dual engine performance, we assume that any active component operates at the radiative limit (i.e. without any non-radiative losses). The electrical power generated by a component i facing a component j is then

P i = U i • e( Ṅj -Ṅi ) = µ i ( Ṅj -Ṅi ). (3)
The total power output is therefore

P = (µ c -µ h )( Ṅh -Ṅc ). (4)
Since q source = q h -q c + P h , the overall heat engine efficiency is

η = P q source = (µ c -µ h )( Ṅh -Ṅc ) q h -q c -µ h ( Ṅh -Ṅc ) . (5)
Using these constitutive equations, we start by studying how the power output of the dual engine varies in the (µ h ,µ c ) plane. The results obtained considering T h = 600 K and E g = 5k B T h are illustrated in Fig. 2a. The dual system can also operate as a heat pump [START_REF] Santhanam | Thermoelectrically Pumped Light-Emitting Diodes Operating above Unity Efficiency[END_REF][START_REF] Radevici | Thermophotonic cooling in GaAs based light emitters[END_REF], and we indicate in a blue colour scale the cooling power. Note that such high cooling powers can only be achieved because below-bandgap radiation is neglected [START_REF] Châtelet | Performances of far and near-field thermophotonic refrigeration in the detailed-balance approach[END_REF]. For the bandgap selected, the dual system is able to operate both as a heat engine and a heat pump in three of the four quadrants: namely, the TPX quadrant (upperright), the TRPV quadrant (upper-left) but also the "TRNEL" quadrant, in which a cold negative electroluminescent (NEL) diode consumes electrical power to limit the radiation sent to the hot facing TR cell. To the best of our knowledge, the TRNEL device has never been mentioned in literature. Globally, both the power output and the cooling power increases with µ c , while P also increases with µ h : the maximum power point (MPP) is therefore located in the TPX quadrant, while the maximum cooling power is reached in TRPV operation for µ c → E g . Despite µ c → E g , the cooling power converges towards a constant value, as demonstrated in [START_REF] Zhao | Chemical potential of photons and its implications for controlling radiative heat transfer[END_REF]. Additionally, note that for low E g , the TPX device becomes unable to perform cooling (see Section II of Supp. Mat.

[10]).

Note how the dual system does not switch directly from heat engine to heat pump operation as µ varies. Indeed, there is a narrow region in-between where the system is neither able to generate electrical power nor to cool the cold source. This gap can especially be observed when E g -µ h becomes lower than k B T h , but globally narrows down as the bandgap increases: in the limit of infinite bandgaps, the two operating regions become adjacent as q h = q c and Ṅh = Ṅc are equivalent. This is schematised in Fig. 2b. In this situation, the term related to the lowest-order polylogarithm dominates in the expressions provided in Eq. [START_REF] Lapotin | Thermophotovoltaic efficiency of 40%[END_REF]. By linearising Li 1 , we obtain that the transition from one region to the other -which corresponds to open-circuit conditions -occurs approximately for

µ c = T c T h µ h + 1 - T c T h E g + k B T c ln T h T c , (6)
this linearised expression being a very good approximation as long as E g -µ i ≫ k B T i . The other side of the power production region is simply delimited by the condition µ h = µ c . These two expressions have already been mentioned for TPX devices in [START_REF] Legendre | GaAs-based near-field thermophotonic devices: Approaching the idealized case with one-dimensional PN junctions[END_REF].

To quantify and compare the performance of TPX, TRPV and TRNEL devices, we provide in Fig. 3 the η -P plots obtained by varying both µ h and µ c , considering three different bandgaps and a heat source temperature of 600 K. For each device, the shaded area corresponds to achievable operation, while the full line is the envelope of this area. Looking at the three cases considered, it is first interesting to notice that while the shape of the admissible η -P area changes significantly when considering each engine individually, it remains mostly unchanged with varying bandgaps for the full dual radiative engine. For any of the bandgaps considered, the efficiency at maximum power remains mostly constant (between 28% and 34% of the Carnot efficiency η C = 1 -T c /T h ), while the maximum efficiency is always η C and is reached at zero power.

If we now compare the different engines, TRNEL appears to be the optimal choice to maximise the efficiency, being always able to reach the Carnot efficiency. Note that this is true only because it operates at the radiative limit: in Section III of Supp. Mat. [10], we analyse the impact of non-radiative losses on performance by considering a spectrally flat quantum efficiency, and pinpoint that the interest of TRNEL operation vanishes when the quantum efficiency goes significantly below unity. Even at the radiative limit, the advantage of TRNEL operation shrinks for large bandgap: as E g → ∞, all radiative engines are able to approach the Carnot efficiency in opencircuit conditions. This has already been demonstrated for TR [START_REF] Strandberg | Theoretical efficiency limits for thermoradiative energy conversion[END_REF][START_REF] Pusch | Fundamental Efficiency Bounds for the Conversion of a Radiative Heat Engine's Own Emission into Work[END_REF] and TPV [START_REF] Roux | Main performance metrics of thermophotovoltaic devices: analyzing the state of the art[END_REF], but can in fact be shown for any radiative engine. In the limit of an infinite bandgap (q h -q c )/( Ṅh -Ṅc ) → E g , leading η to be expressed as (µ c -µ h )/(E g -µ h ). Using Eq. ( 6), and keeping in mind that E g -µ h ≫ k B T h , we obtain that the efficiency in FIG. 3: η -P plots of dual radiative engines at the radiative limit, for T h = 600 K and for various bandgaps.

open-circuit conditions equals η C . In contrast, to maximise the power output, TPX is almost always the best candidate, TRPV becoming optimal only for very low bandgap energies (here, for E g < k B T h ≈ 0.05 eV) which are not achievable physically. This remains true for lower or higher heat source temperatures (see Section IV of Supp. Mat. [10]). TPX devices generally outperform other radiative engines in terms of power output because the hot emitter operates as an LED: it is therefore able to largely enhance its emission by electroluminescence, increasing consequently the various energy flows in the system. For high enough bandgaps, this can even cause P max to exceed σT 4 h (see Fig. 3c), which is physically impossible for single radiative engines.

While TRNEL devices allow maximising efficiency and TPX devices power, TRPV systems can in some conditions provide interesting trade-offs between power and efficiency (see Fig. 3b). However, this can only be observed for low bandgaps (a few k B T h at most) which are hardly achievable for the heat source temperature considered.

The study of η -P plots directly highlights that the maximum efficiency of dual radiative engines is always the Carnot efficiency, reached for zero power output. Regarding the MPP, it was observed that P max increases significantly with E g while η MPP varies only slowly, but no analytical expressions of these quantities have yet been formulated. Therefore, in the following, we focus on the analytical derivation of P max and η MPP . To do so, we consider once again that E g → ∞, since this allows reaching the largest possible power output (see Section V of Supp. Mat. [10]). By doing so, the Li 1 term dominates the expressions provided in Eq. ( 2). Using that

µ c -µ h = x h k B T h -x c k B T c
, and defining X i as exp(-x i ), the power output can be expressed as:

P = E 2 g k 2 B T 2 h 4π 2 c 2 h3 ln(X h ) - T c T h ln(X c ) × ln(1 -X h ) - T c T h ln(1 -X c ) , (7)
since Li 1 (x) = ln(1 -x). By using the symmetry of this expression with respect to X h and X c , one can show that P reaches its maximum for X h = X c = 1/2, thus for µ i = E g -ln(2)k B T i (see Section VI of Supp. Mat.

[10]). Consequently, the maximum power output is

P max = 1 h ln(2)E g k B (T h -T c ) 2πch 2 , (8)
and varies quadratically with both the bandgap energy and the temperature difference. Such variations were already pointed out in [START_REF] Zhao | Chemical potential of photons and its implications for controlling radiative heat transfer[END_REF], although without a complete closed-form expression. Since E g ≫ k B T i , both chemical potentials are greater than 0 and the maximum is reached in TPX regime, consistently with the results from Fig. 3.

To derive a closed-form expression of the efficiency at maximum power, we use that x h,MPP = x c,MPP = ln(2) to obtain (µ c,MPP -µ h,MPP )/(E g -µ h,MPP ) = η C . Dividing both the numerator and the denominator of Eq. ( 5) by ( Ṅh -Ṅc )(E g -µ h ), one gets

η MPP = η C 1 + 1 ρ -1 1 ln(2) Eg kBT h , (9)
where ρ = E g ( Ṅh -Ṅc )/(q h -q c ) corresponds to the fraction of radiative energy being useful to optoelectronic conversion. To express it, both Li 1 and Li 2 terms are necessary. One obtains

1 ρ -1 1 ln(2) E g k B T h ∼ Eg→∞ (2-η C ) 1 ln(2) Li 2 (1/2) Li 1 (1/2) . (10)
The polylogarithmic terms having closed-form expressions for x = 1/2, the efficiency at maximum power obtained as

E g → ∞ is η MPP,Eg→∞ = η C 1 + (2 -η C )χ , (11)
FIG. 4: Variation of the efficiency at maximum power with the heat source temperature. For simplicity, only the TPX quadrant has been considered.

χ being a constant equal to 1 2 1 6 π ln(2) 2 -1 ≈ 1.21.
The temperature variation of η MPP,Eg→∞ is provided in Fig. 4 (black line), and matches well the numerical results obtained for bandgaps larger than 100k B T h . It also gives a good estimate of the efficiency obtained for standard bandgaps, as long as E g ≫ k B T h : for T h = 600 K, η MPP,Eg→∞ = 17.7% while η MPP,1 eV =17.1% and η MPP,0.52 eV = 16.6%.

To better understand how efficient dual radiative engines are at maximum power, one can compare Eq. [START_REF] Melnick | From thermoelectricity to phonoelectricity[END_REF] with classical upper bounds for η MPP . We choose to perform the comparison with the Novikov-Curzon-Ahlborn [START_REF] Novikov | Efficiency of an atomic power generating installation[END_REF][START_REF] Curzon | Efficiency of a Carnot engine at maximum power output[END_REF] efficiency η NCA = 1 -T c /T h and the Schmiedl-Seifert [START_REF] Schmiedl | Efficiency at maximum power: An analytically solvable model for stochastic heat engines[END_REF] efficiency η SS = 2η C /(4 -η C ), which were found to be efficiency bounds for endoreversible and exoreversible thermoelectric generators respectively [START_REF] Apertet | Irreversibilities and efficiency at maximum power of heat engines: The illustrative case of a thermoelectric generator[END_REF]. For the temperatures previously considered, both efficiencies are close to 29%, hence 11 percent points higher than η MPP,Eg→∞ . This significant difference, which highlights the presence of additional losses in radiative engines, can be attributed to thermalisation losses -i.e., to the fraction of radiative energy exchanged which is useless to optoelectronic conversion. To prove this, we now consider that the radiation is quasi-monochromatic around E g , which makes the ratio ρ equal to 1. In this scenario, µ i,MPP is found to tend towards E g (see Section VII of Supp. Mat. [10]). Then, the Bose-Einstein distributions can be simplified using that [exp(x)-1] -1 ∼ x -1 around 0. Setting any of the two partial derivatives of P with respect to µ i to zero leads to

η MPP,hω=Eg = 1 - T c T h , (12)
exactly the Novikov-Curzon-Ahlborn efficiency. δE being the radiation spectral bandwidth, the maximum power FIG. 5: Variation of the maximum power and related efficiency for varying spectral bandwidth, considering T h = 600 K. For simplicity, only the TPX quadrant has been considered.

output is then

P max = 1 h E g √ k B ( √ T h - √ T c ) 2πch 2 δE, (13)
an expression similar to Eq. ( 8). Since P goes to zero as δE → 0, there is a trade-off between power and efficiency as the bandwidth varies, as illustrated in Fig. 5: to achieve non-zero output power, the efficiency must fall below the usual bounds. It is noteworthy that the efficiency starts to decrease for bandwidths as low as few meV (corresponding to a quality factor Q = E g /δE close to 100 for E g = 0.52 eV), while reaching the broadband limit for a bandwidth of few tenths of eV (i.e. for Q ≈ 1 considering E g = 0.52 eV). If the efficiency at maximum power is too low for a given application, two main leverages are thus available to increase it, although at the expense of power: decrease the radiation bandwidth, or change µ h and µ c to move in the broadband η -P plots provided in Fig. 3, which can allow exceeding the aforementioned bounds [START_REF] Giteau | Thermodynamic performance bounds for radiative heat engines[END_REF]. The interest of each of these leverages depends on the bandgap, and on how far the system operates from the radiative limit. In some cases, spectral filtering can allow extending the region of achievable operating conditions, limiting the power loss undergone when high efficiency is required (see Section VIII of Supp. Mat [10]).

In conclusion, we have studied the power output and efficiency achievable by dual radiative heat engines, especially when they operate at the radiative limit. A unified view allows shining light on their similarities and respective merits. In particular, TRNEL devices are found to systematically reach Carnot efficiency, while TPX devices are almost always the best dual engines in terms of power output, and offer the broadest range of operating conditions of all dual engines for bandgaps over a few k B T h . We have analytically derived an upper bound for the maximum power and related efficiency, the latter mostly lying between 30% and 40% of η C and being several percent points below usual efficiency bounds due to thermalisation losses (11 points below for T h = 600 K). Interestingly, spectral filtering can mitigate part of the power loss when high efficiencies are targeted. In practice, below-bandgap radiation, more rigorous nonradiative losses and resistive losses shall be included, as well as thermal resistance effects at the thermostats which reduce the operating temperature difference T h -T c [START_REF] Legendre | Theoretical and numerical analysis of near-field thermophotonic energy harvesters[END_REF]. In addition, it would be worth investigating how the performance of such engines changes in the near field, where radiative emission exceeds the modified Planck law [START_REF] Zhao | High-performance near-field thermophotovoltaics for waste heat recovery[END_REF][START_REF] Legendre | GaAs-based near-field thermophotonic devices: Approaching the idealized case with one-dimensional PN junctions[END_REF][START_REF] Legendre | Theoretical and numerical analysis of near-field thermophotonic energy harvesters[END_REF].

This work has received funding from the European Union's Horizon 2020 research and innovation programme under Grant Agreement No. 951976 (TPX-Power project). The authors thank T. Châtelet, P. Kivisaari, O. Merchiers, J. Oksanen and J. van Gastel.

FIG. 2 :

 2 FIG.2:(a) Performance of dual radiative systems operating as heat engines or heat pumps, at the radiative limit and for E g = 5k B T h = 0.26 eV. (b) Operating regions of dual radiative systems in the limit of an infinite bandgap.


Xiaohao Sun

Luxia Yu

Liang Yue

Kun Zhou

Frédéric Demoly

Renee Ruike

H Jerry Zhao

Qi

Ruike Renee Zhao

H Jerry Qi

Machine learning and sequential subdomain optimization for ultrafast inverse design of 4D-printed active composite structures

Keywords: active composites, morphing structures, 4D printing, inverse design, machine learning

de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Shape transformation of active composites (ACs) depend on the spatial distribution and active response of constituent materials. Multimaterial 3D/4D printing [START_REF] Ge | Active materials by four-dimension printing[END_REF][START_REF] Kuang | Grayscale digital light processing 3D printing for highly functionally graded materials[END_REF][START_REF] Cheng | Centrifugal multimaterial 3D printing of multifunctional heterogeneous objects[END_REF][START_REF] Yue | Single-vat single-cure grayscale digital light processing 3D printing of materials with large property difference and high stretchability[END_REF][START_REF] Yue | Cold-programmed shape-morphing structures based on grayscale digital light processing 4D printing[END_REF] allows for easy implementation of voxel-level distribution of active materials/properties, thus permitting a vast space of active shape changes that can be achieved. To fully exploit this space, inverse design (i.e., finding the optimal material distribution to achieve a target shape) is highly desired, which is yet challenging due to the tremendous design space [START_REF] Demoly | The status, barriers, challenges, and future in design for 4D printing[END_REF][START_REF] Sun | Perspective: Machine learning in design for 3D/4D printing[END_REF]. To tackle these challenges, gradient-based and gradient-free methods have been developed. For instance, the gradient-based topology optimization (TO) has been used for designing active shape changes of ACs [START_REF] Maute | Level Set Topology Optimization of Printed Active Composites[END_REF][START_REF] Geiss | Combined Level-Set-XFEM-Density Topology Optimization of Four-Dimensional Printed Structures Undergoing Large Deformation[END_REF][START_REF] Tanaka | Turing pattern-based design and fabrication of inflatable shapemorphing structures[END_REF] or optimizing structure compliances [START_REF] Zolfagharian | Topology-Optimized 4D Printing of a Soft Actuator[END_REF]. The TO, however, typically requires the complicated derivation of gradients and may encounter difficulties when the ACs involve geometric or material nonlinearity. Alternatively, the gradient-free evolutionary algorithm (EA) has also been adopted for designing active shape changes of ACs [START_REF] Hamel | Machine-learning based design of active composite structures for 4D printing[END_REF][START_REF] Wu | Evolutionary Algorithm-Guided Voxel-Encoding Printing of Functional Hard-Magnetic Soft Active Materials[END_REF][START_REF] Athinarayanarao | Computational design for 4D printing of topology optimized multi-material active composites[END_REF]. In these EA implementations, finite element (FE) methods are used for numerous shape predictions (referred to as FE-EA), leading to high computational cost.

Although reduced-order models have been developed to replace FE for faster EA, they suffer from either relatively low accuracy (against FE) [START_REF] Sossou | Design for 4D printing: Modeling and computation of smart materials distributions[END_REF] or limited design freedom (1D material distribution) [START_REF] Wang | Evolutionary design of magnetic soft continuum robots[END_REF]. Accurate and efficient inverse design strategies remain to be developed.

Machine learning (ML) [START_REF] Guo | Artificial intelligence and machine learning in design of mechanical materials[END_REF] opens up new avenues for developing fast, computationally affordable, and high-fidelity models for both forward prediction and inverse design of a variety of material responses. Existing works mainly focused on optimizing or predicting mechanical properties of materials, such as strength and toughness of composites [START_REF] Gu | Bioinspired hierarchical composite design using machine learning: simulation, additive manufacturing, and experiment[END_REF][START_REF] Chen | Effect of Constituent Materials on Composite Performance: Exploring Design Strategies via Machine Learning[END_REF], stress and strain fields of composites [START_REF] Yang | Deep learning model to predict complex stress and strain fields in hierarchical composites[END_REF], Poisson's ratio of auxetic metamaterials [START_REF] Wilt | Accelerating Auxetic Metamaterial Design with Deep Learning[END_REF], responses of soft pneumatic robots [START_REF] Zolfagharian | 4D printing soft robots guided by machine learning and finite element models[END_REF], among others [START_REF] Roach | Utilizing computer vision and artificial intelligence algorithms to predict and design the mechanical compression response of direct ink write 3D printed foam replacement structures[END_REF][START_REF] Rawat | A novel topology design approach using an integrated deep learning network architecture[END_REF]. However, there is limited work on exploiting ML for designing active shape-change responses of ACs [START_REF] Sun | Perspective: Machine learning in design for 3D/4D printing[END_REF]. For voxelized AC beams, Zhang et al. [START_REF] Zhang | Finite-Element-Based Deep-Learning Model for Deformation Behavior of Digital Materials[END_REF] explored performances of multiple ML models for the forward prediction problem. Our recent work [START_REF] Sun | Machine Learning-Evolutionary Algorithm Enabled Design for 4D-Printed Active Composite Structures[END_REF] demonstrated that the recurrent neural network (RNN) is particularly suited for the beam problem as it inherently preserves a sequential dependency, similar to that in beam deformation. The RNN-based ML model thus demonstrated exceptional accuracy in forward shape predictions based on material distributions. The ML was then integrated with EA (ML-EA), enabling fast inverse designs for complicated target 2,5D shapes, typically in 10-60 minutes. Although ML-EA has significantly outperformed conventional methods such as FE-EA, it also has two major limitations. First, there remains a strong demand for faster or even real-time inverse designs for AC beams. One key limiting factor for the design speed is that EA seeks solutions in a vast design space comprising all voxels, which is relatively inefficient. Second, an ML model, once trained, is specific to certain beam length and voxel configuration, making it inapplicable to the design of target shapes that have different lengths (hence different voxel numbers) or finer features (such as large curvature) beyond the design space of current voxels.

This work addresses these issues by leveraging the sequential dependency of beam deformations. We first present a design approach that seamlessly combines an RNN ML model and a sequential subdomain optimization (SSO) algorithm to realize ultrafast inverse designs for target shapes of various complexities, producing results in second(s), which is hundreds or thousands times faster than conventional ML-EA method. In addition, the sequential characteristic allows for accurate predictions for beams with fewer lengthwise voxels than that of the training data. It also permits the splicing of multiple predicted shapes, extending accuracy to beams with more lengthwise voxels. Moreover, this splicing strategy, when integrated with ML-SSO, also permits ultrafast inverse designs for a spectrum of target intricate shapes with different lengths, using lengthwise voxel configurations different from that for the ML training. Finally, the highly efficient ML-SSO is employed for the inverse design of active shape changes of 4D-printed lattice structures. Our approach thus facilitates the development of 4D printing towards intelligent and streamlined design and fabrication of various shape-morphing AC structures.

Results and discussions

Physical problem and ML model

The physical problem is illustrated in Figure 1a. Here we consider an AC beam whose left end is fixed. It is made of two materials: an active material '1' expands under the thermal actuation while a passive material '0' does not. The beam is composed of Nx (length-or xdirection) × Ny (thickness-or y-direction) voxels, each assigned with an active or passive material. The active shape change, which is induced by the expansion mismatch under the thermal actuation, depends on the material distribution (Figure 1a). With a dataset generated by FE simulations, an RNN-based ML model is trained to predict actuated shapes for given material distributions (Figure 1b). RNN can use the "past" information to predict the "current" response and thus is highly capable for learning sequential data. It is noted the deflection of the AC beam has a similar sequential characteristic, i.e., the displacement of any points only depends on the points to their left (or before them); the displacement of the points to their right (or after them) would not affect their displacement. This feature inspires us to use RNN. For each step, RNN receives a single column of voxel encodings and produces the coordinates of the column (sampled on mid-axis of the beam). By repeating this process, RNN sequentially processes all voxel columns to predict the shape of entire beam (Figure 1b). Details on the FE model, dataset, and the architecture and training of the ML model are provided in Sections 4.1 and 4.2.

As shown in our previous study [START_REF] Sun | Machine Learning-Evolutionary Algorithm Enabled Design for 4D-Printed Active Composite Structures[END_REF], the RNN ML model demonstrates high accuracy for cases with different numbers of voxels. In this work, we adopt the ML model trained on 24 (Nx) × 4 (Ny) voxels. Notably, the coordinates (x, y) are sampled from the mid-axis, rightmost mesh point of each voxel column to represent a shape, which is different from the middle point as used in our previous work (Figure 1b). This choice of sampling points has two merits. First, it achieves higher prediction accuracy with R 2 = 0.99987 for x and 0.99986 for y, as shown in regression (or density scatter) plots of the ground-truth versus ML-predicted coordinates of the testing dataset (Figure 1c). The increased accuracy is attributed to the better sequential characteristic of the column rightmost points. Second, the new sampling points enable the ultrafast SSO optimizations as will be introduced in Section 2.2.

As mentioned above, since the left end of the actuated beam is fixed, the coordinate of a voxel column is dependent on the voxels in the current column and the left of the current column (nearer the fixed end). Similarly, the architecture of RNN determines that the "current" output is dependent on the "current" and "past" inputs. To illustrate this feature, we display the gradients (or sensitivities) of RNN-predicted coordinates (x, y) of the 12th and 24th voxel columns with respect to all voxels (i.e., ∂xi/∂Mv and ∂yi/∂Mv for i=12 and 24, where v=1, 2, …, 4×24 represents the voxel number, Mv represents the material encoding of voxel v), based on a state where all voxels are initially passive materials (Mv=0) (Figure 1d). The gradients are calculated through automatic differentiation as described in Section 4.3. We see that (x12, y12) depend on the first 12 voxel columns only, while (x24, y24) depend on all the voxels. The results evidently showcase that the RNN inherently preserves the sequential data dependency of the beam shape.

ML-SSO approach for ultrafast inverse design

Although the ML model achieves high accuracy in the forward prediction, applying it to the inverse map (from shape to design) is still challenging due to the "one-to-many" characteristic (i.e., one actuated shape could be achieved by distinct designs) [START_REF] Sun | Perspective: Machine learning in design for 3D/4D printing[END_REF]. This is particularly the case considering the large design space and the high-dimensional data for shape representation. To address this challenge, a straightforward approach is to combine ML with optimization algorithms such as EA [START_REF] Sun | Machine Learning-Evolutionary Algorithm Enabled Design for 4D-Printed Active Composite Structures[END_REF]. However, as a stochastic search algorithm, EA seeks solutions in a vast design space comprising all voxels, which is inefficient. In this work, we present an ML-empowered SSO approach which can further leverage the inherent sequential characteristic. Instead of directly optimizing the entire domain, this approach sequentially optimizes subdomains, composed of Nsub columns, from left to right (Figure 2a). In each optimization step, SSO receives the coordinates of Nsub target points (highlighted red "o") and optimizes the material distribution in the corresponding subdomain (highlighted yellow domain) such that the RNN-predicted coordinates of actuated points (highlighted blue "+") match with the targets. Once a subdomain, e.g., (i-1), is optimized, the target points and design domain move one step forward for the next subdomain, i, as seen by the one-step sliding of the highlighting window in Figure 2a. Due to sequential dependency, the optimization of subdomain i would not affect the optimized preceding subdomains (those to the left of subdomain i) except the overlapping domain of (i-1) and i (see the rightmost panel of Figure 2a). Here, we choose Nsub to range from 1 to 3. When Nsub=1, there is no overlap between two consecutive subdomains. As will be detailed later, a larger Nsub can improve the optimization accuracy. For the optimization of individual subdomains, we employ a brute-force approach, calculating the subdomain shapes of all possible designs and retaining the optimal one. This ensures the local optimality of found solution while not requiring many shape evaluations as the chosen Nsub is small. For larger Nsub, the subdomain EA or other stochastic search algorithms may be used to reduce computational cost. The loss function Lsub for each subdomain is defined as [START_REF] Ge | Active materials by four-dimension printing[END_REF] where (xsub,j, ysub,j) are the ML-predicted sampled point coordinates for the j-th column of the considered subdomain, !𝑥 # !"#,% , 𝑦 # !"#,% & are the corresponding target coordinates.

Compared to the ML-EA on entire domain, ML-SSO reduces computational cost in two aspects. First, the design space of a single subdomain is much smaller than the whole domain, thus drastically reducing the number of potential design evaluations required. For

( ) ( ) 2 2 , , , , 1 1 ˆŝub
N sub sub j sub j sub j sub j j sub L x x y y N = é ù = - + - ê ú ë û å
instance, for Nsub = 2, the brute-force method assesses 256 (=2 Ny×Nsub =2 8 ) potential subdomain designs per step, totaling 5888 (=256×23, where 23 = (Nx-Nsub+1) is the number of subdomains to optimize) evaluations, which is much smaller than the number (typically 75,000-375,000 [START_REF] Sun | Machine Learning-Evolutionary Algorithm Enabled Design for 4D-Printed Active Composite Structures[END_REF]) of global ML-EA. Second, these are the numbers of potential designs that need to be evaluated. In ML-SSO, the evaluation of each candidate design (comprising Nsub columns only) is also more efficient. This is because during each optimization step, RNN permits only processing the current subdomain and predicting its coordinates (xsub,j, ysub,j), j ∈ [1,Nsub], based on the stored hidden state, thereby bypassing the need to evaluate the preceding or subsequent subdomains. This further significantly reduces computational costs compared to ML-EA which demands a full-domain shape evaluation for each candidate design.

Performance of ML-SSO on numerically generated target shapes

We now consider multiple numerically generated target shapes of different complexities to evaluate the performance of ML-SSO with different Nsub values. In this sub-section, we focus on the case with 24×4 voxels. The first two target shapes, the one-period and twoperiod sinusoidal shapes, are optimized using Nsub = 1. To provide a more intuitive understanding of ML-SSO, in Figure 2b and 2c, we show the optimized and target shapes with symbols of varying colors to visualize the optimization sequence. Each optimized point ("+") is the one closest to its the corresponding target point ("o") among all possible designs for the specific subdomain. For both two targets, excellent agreements between the optimized shapes and the targets are achieved. Although some discrepancies can be observed, e.g., in the middle of the two-period shape where the curvature is relatively large, the subsequent optimized points can still accurately capture the target shape (Figure 2c), demonstrating the high capability of our approach.

Next, two more complicated targets are considered: the three-period shape and the halfbutterfly shape. As shown in Figure 2d and 2e, the sensitivity of optimized shapes to the Nsub is studied. For both targets, Nsub=1 leads to notable discrepancies between the optimized and target shapes. With the increased Nsub, the optimization accuracy is significantly improved. This is because ML-SSO operates on a greedy strategy, selecting the best shape at each step. Such an approach might lead to a short-sighted problem for regions with sharp turns, missing the target of the subsequent subdomain in the favor of the best match for the current one. Increasing Nsub effectively broadens the optimization foresight and can, to some extent, mitigate this problem, although it concurrently escalates the computational cost. To gain more insights, we show in Figure 2d and 2e the optimization sequence for the cases with Nsub=2, where each voxel column (except the first and last ones) is involved in two consecutive subdomains and thus optimized twice (see, for example, the overlapping domain of (i-1) and i in the rightmost panel of Figure 2a).

The temporarily optimal point denotes the one achieved in the first optimization of its involved subdomain and is depicted in gray "+" symbols. The discrepancy between the temporarily optimal point and the optimal point often appears in regions with sharp turns, manifesting the short-sighted problem. We further quantitively compare the performance of ML-SSO with different Nsub values on these four targets in terms of accuracy (Figure 3a) and time cost (Figure 3b). The ML-EA optimization results are already available [START_REF] Sun | Machine Learning-Evolutionary Algorithm Enabled Design for 4D-Printed Active Composite Structures[END_REF] and also compared here. The accuracy is quantified by the root-mean-squared errors (RMSEs)

( ) ( ) 2 2 1 1 ˆRMSE 2 x N i i i i i x x x y y N = é ù = - + - ë û å
EA which takes 660 s. For the most complex three-period shape, ML-SSO with Nsub=3 takes 11.8s and achieves much better RMSEs than that of ML-EA, which takes 3240 s. In short, ML-SSO outperforms or matches ML-EA in optimization accuracy while being significantly (at least two orders of magnitude) faster. This demonstrates very high efficacy of our approach in solving inverse material design problems for complicated target shapes, which is attributed to the seamless integration of RNN-ML and SSO that fully leverages sequential characteristics for highly efficient optimal design searches.

Streamlined 4D printing design and fabrication based on hand-drawn lines

Next, we apply the ultrafast ML-SSO approach to hand-drawn target shapes to enable a streamlined design-fabrication paradigm for 4D printing. In general, these target shapes are more challenging to optimize, because the discrete target points are uniformly sampled from the identified shape, which could have inappropriate spacing and thus produce difficulties for the optimization. Following our previous work [START_REF] Sun | Machine Learning-Evolutionary Algorithm Enabled Design for 4D-Printed Active Composite Structures[END_REF], computer vision (CV) algorithms are employed for the automatic identification of target shapes from hand-drawn lines. The ML-SSO is then utilized for the inverse design, and the optimized material distributions can be readily 4D-printed using the grayscale digital light processing (g-DLP) [START_REF] Kuang | Grayscale digital light processing 3D printing for highly functionally graded materials[END_REF][START_REF] Yue | Single-vat single-cure grayscale digital light processing 3D printing of materials with large property difference and high stretchability[END_REF] technique. The material system is different from that of our previous work [START_REF] Sun | Machine Learning-Evolutionary Algorithm Enabled Design for 4D-Printed Active Composite Structures[END_REF].

More details on the printing method, material and actuation mechanism are provided in Section 4.4. Figure 4 shows the ML-SSO design and fabrication results for four different target shapes.

The ML-EA optimized designs are available [START_REF] Sun | Machine Learning-Evolutionary Algorithm Enabled Design for 4D-Printed Active Composite Structures[END_REF] and also presented here for the evaluation of ML-SSO performance. Shown in five columns (from left to right) are hand-drawn lines, optimized shapes by ML-SSO (with different Nsub values, solid lines) and ML-EA (dashed lines) compared against the targets (symbols), the corresponding performance in terms of RMSE and time cost, the optimal ML-SSO designs selected (with greater Nsub, and RMSE < 1 mm) out of those in column 3 (arrows), and experimentally 4D-printed, actuated shapes.

In column 5, all printed strips are initially flat and only the actuated shapes are shown. For all the targets, the optimized shapes by both ML-SSO and ML-EA agree well with the target shapes, whilst ML-SSO is much faster than ML-EA. The optimal ML-SSO designs are finally converted into grayscale slices for 4D printing through g-DLP, and the printed strips transform into target shapes upon actuation. These results demonstrate an ultrafast, streamlined design-fabrication paradigm for 4D printing.

Active beams with varying number of lengthwise voxels: forward predictions

The preceding sections demonstrate ultrafast inverse designs enabled by ML-SSO for both numerically generated and hand-drawn target shapes. However, these designs have been constrained to fixed beam length (80 mm) and voxel configuration (24×4 voxels) on which the ML model is trained, implying a relatively limited design space. For example, a configuration of Nx=24 lengthwise voxels might struggle to capture complicated target shapes, such as the six-period shape, which would be very challenging for Nx=24 as it would require four voxels in the length direction to capture one period. Therefore, more lengthwise voxels (Nx > 24) are needed to tackle this design problem.

To further enhance design capabilities, we aim for the ML-based methods to address both As shown in Figure 5a, an input sequence (referring to a material distribution) of any length Nx can be divided into multiple subsequences with overlaps, each with a length of 24, except for the last subsequence which can be shorter than 24. We use ML to predict the shape of all subsequences, then concatenate them in order to obtain the full-sequence shape.

The predicted shapes of all subsequences are initially based on a fixed boundary condition on the left end, which we refer to as being in their local spaces, while the shape of the first subsequence is in the real space. We sequentially transform each subsequence shape into the real space by matching the coordinates in the overlapping regions between two consecutive subsequences. Here, we use the length of overlapping region to be Noverlap=6.

Note that the correct local concatenation angle cannot be determined if there is no overlapping region.

Next, we perform forward predictions on multiple benchmark shapes of different lengths, including Nx=12, 24, 48, 72, 96 (Figure 5b). The input sequences with Nx≤24 are directly predicted using the RNN-ML without accuracy degradation. The input sequences with Nx>24 are obtained by concatenating multiple sequences of the three-period shape and predicted using the ML-splicing strategy. Directly applying ML to these sequences would yield less accurate results. By using the ML-splicing strategy, we can achieve very high prediction accuracy. Note that the ML prediction errors will inevitably accumulate with increasing Nx and splicing times. Nonetheless, the ML-splicing predictions remain accurate for Nx up to 72, three times Nx=24; a noticeable discrepancy only arises at Nx=96. Our strategy thus demonstrates excellent robustness significantly (three times) beyond the sequence length (Nx=24) of the ML training set.

Active beams with varying number of lengthwise voxels: inverse designs

Next, we combine ML-SSO with the splicing strategy for inverse design problems. The splicing strategy works in a slightly different way from the forward prediction. In the forward prediction, each subsequence shape is transformed into the real space to infer the correct shape. In the inverse design, instead of transforming the shapes of potential designs into the real space for comparisons against the target, we transform the target shape (which is initially in the real space) into the local space for the optimization. This enables the seamless use of ML-SSO for each subsequence, bypassing unnecessary shape transformations and improving the optimization speed. We use multiple complicated target shapes (i.e., involving relatively sharp turns when represented using Nx=24 sampling points)

to demonstrate the capability of our splicing-empowered ML-SSO. Figure 6 shows the optimization results for these target shapes obtained with different Nx values. The optimized shapes together with the corresponding RMSEs and time cost are presented. FE simulations for the optimized designs are also performed, and the FE predicted shapes and their RMSEs are provided. The first two targets are the most complex target shapes considered so far: the three-period shape (Figure 6a) and hand-drawn mountain (Figure 6b). With the original Nx=24, Nsub=2 is not sufficient to accurately capture the two targets (Figure 2d and4); using Nsub=3 significantly increases the accuracy (with RMSEs of 0.21 mm and 0.96 mm) but also takes longer optimization time (11.8 s) (Figure 6a and6b). By using greater Nx values (48 and 40 respectively), which offer more sampling points, these complex targets with sharp turns can be rapidly captured with Nsub=2, achieving similar accuracies (with RMSEs of 0.28 mm and 0.97 mm) in much shorter times (1.97 s and 1.62 s) (Figure 6a and6b). To put these results in perspective, for the three-period shape, FE-EA approach was estimated to need 227 days for the inverse design, the ML-EA approach obtained the design in 54 min, our new approach requires a merely 1.97 s.

To further evaluate the capability of our approach, we consider a more complex six-period target shape obtained by chaining two three-period shapes (Figure 6c). The original Nx=24 fails to capture this target (RMSE=3.33 mm) even with Nsub=3 (which takes 11.8 s), and therefore greater Nx values are used. We then use Nx=48 with Nsub=2 and achieve an improved accuracy (RMSE=2.01 mm) in a much shorter time (1.98 s), but an evident discrepancy can still be observed. By further increasing Nx=60, even a smaller Nsub=1 successfully yields a high-precision optimized shape (RMSE=1.38 mm) in only 0.53 s (Figure 6c). As the last example, we take the profile of a "Chinese dragon" art as the target shape and again achieve excellent optimization results (RMSE=0.29 mm) very rapidly (0.42 s with Nx=40 and Nsub=1) (Figure 6d).

These results show that although we train the ML model on a certain voxel configuration, the splicing strategy allows for the ultrafast and accurate forward predictions and inverse designs for structures with variable Nx values and lengths by leveraging the inherent sequential characteristics.

Applications to rapid design of 4D-printed lattice structures

The inverse designs demonstrated so far have been constrained to the AC beam structures.

However, the voxel-level, ultrafast inverse shape-change designs for more complex structures such as lattice structures are highly desired, which could have broader applications but is also more challenging due to increased number of involving voxels in an entire lattice. In this sub-section, we further show that by leveraging appropriate symmetries, the ML-SSO can be applied in designing the active shape transformation of 4D-printed lattice structures. We focus on a periodic square lattice and optimize the material distribution on its edges so that they achieve target shapes under actuation. The design strategy can be also applied to other lattice structures, such as triangular and hexagonal lattices. First, we consider an arbitrary target edge shape. The material distribution is designed to satisfy 4-fold rotational symmetry about each joint to prevent mismatches in structural deformation. Due to the size limit of the printing stage, we use half of the hand-drawn mountain as a target (Nx=12). We adapt its optimized design for the lattice structure, validate the design through FE simulations, and then 4D print the lattice structure (Figure 7a). The FE results show that all edges of the lattice transform to the desired shape. In experiments, we only print individual square lattices due to the size limit. These printed lattices, initially square, transform into the target shape upon actuation. When tessellated, the structure matches with the FE result. The tessellation is used solely for improved visualization and does not affect the actuated shape.

Note that for an arbitrary target shape considered above, four (2×2) lattice cells form a periodic unit. Next, we consider anti-symmetric target edge shapes, which are common in lattice structures, such that a single lattice cell would be periodic. In this case, a translational symmetry in the material distribution is satisfied, ensuring the structural compatibility. Due to the anti-symmetry of the edge, the design is carried out on half of the domain. For a simple wavy target (Nx=6 on its half), we perform the inverse design and 4D

printing. The initially square lattices, when actuated and tessellated, agree excellently with the FE predictions (Figures 7b). Furthermore, we consider two more complex target shapes (Nx=12 on the half) and use them for multiple lattice structures with identical target shapes for all edges or distinct targets for horizontal and vertical edges. For all lattices, excellent agreement is achieved among the target, the optimized shape, and the actuated edge shapes of the lattice obtained from FE simulations (Figure 7c). These results demonstrate the high efficacy and flexibility of our design strategy for lattice structures.

It is worth noting that a translational symmetric material distribution can be used to ensure the structural compatibility as long as the target edge shape has equal slopes at two ends, without necessarily being antisymmetric. Additionally, by employing appropriate target shapes for horizontal and vertical edges, our approach can be readily used for designing lattices with target macroscopic deformations, such as horizontal and vertical normal strains, shear strains, and/or rigid rotations. Moreover, our design strategy is also readily applicable to other lattice structures, such as triangular and hexagonal lattices. These results and discussions highlight the high flexibility of our approach in designing the active shape change of lattice structures.

Conclusions

We present an approach for ultrafast inverse design of 4D-printed AC beams by combining ML and SSO. An RNN ML model is rapid and accurate in forward shape predictions based on the material distributions. We then integrate ML with SSO for the inverse design of material distributions based on the target shapes. For multiple target shapes of different complexities, ML-SSO outperforms or is comparable to ML-EA in optimization accuracy while being two to three orders of magnitude faster, delivering results in mere second(s).

The CV-integrated ML-SSO then demonstrates an ultrafast, streamlined design-fabrication paradigm based on hand-drawn targets. Furthermore, although the ML model is trained on a fixed voxel configuration (Nx=24, Ny=4), we present a splicing strategy that achieves accurate shape predictions for beams with varying numbers (Nx) of lengthwise voxels with no need for retraining the ML model. This strategy, when integrated with ML-SSO, enhances adaptability for ultrafast inverse designs, accommodating target shapes of diverse complexities or lengths. For highly complex target shapes, the splicing-empowered ML-SSO, simply with Nx>24, proves to be more robust and rapid than the regular ML-SSO (which requires Nx=24), achieving better or similar optimization results while reducing time cost from 11.8 s to 1.6-2 s. For example, as a comparison, for the benchmark three period shape, the FE-EA method was estimated to need 227 days for the inverse design;

the ML-EA achieved design in 57 min; the new ML-SSO with splicing strategy requires only 1.97 s. Finally, the highly efficient ML-SSO is employed for the inverse design of active shape changes of 4D-printed lattice structures. Our approach thus offers an intelligent design-fabrication paradigm for 4D printing of various shape-morphing AC structures.

Models and Methods

Finite element model for data generation

We perform FE simulations using the commercial software ABAQUS (version 2018, Simulia, Providence, RI). The FE model follows our previous work [START_REF] Sun | Machine Learning-Evolutionary Algorithm Enabled Design for 4D-Printed Active Composite Structures[END_REF] and is briefly summarized here. An active beam with the left end fixed and under the plane strain conditions is considered. It has a dimension of 80 mm long × 1 mm thick and is partitioned into Nx × Ny voxels for material assignment. Both active and passive materials are modeled as incompressible neo-Hookean solids with the same Young's modulus but different coefficients of thermal expansion (0.001 for active and 0 for passive materials). The shape change is induced by applying a 100 o C temperature increase to the entire beam, which results in a linear strain mismatch of 0.1 between the two materials. Mesh convergence study is performed and 960×12=11520 hybrid plane strain (CPE4H) elements are adopted.

The FE model is automatically generated and run through a Python script.

To generate a dataset for ML, we focus on the voxel configuration of 24 (Nx) × 4 (Ny), which has a large design space of 2 96 (≈7.92 × 10 28 ) possible material distributions. 8600 random material distributions are created, and FE simulations are used to obtain their ground-truth actuated shapes. The generated dataset is then split into training, validation, and testing datasets with fractions of 0.7, 0.15, and 0.15. Three datasets (training, validation, and testing) show similar statistical distributions [START_REF] Sun | Machine Learning-Evolutionary Algorithm Enabled Design for 4D-Printed Active Composite Structures[END_REF].

To use RNN, a similar data structure the beam deformation problem is utilized. A column (y-direction) of voxels is the input of a single step and the entire input of a material distribution is restructured into a sequence of Nx columns of voxels. Similarly, a sequence of outputs is formed by coordinates (x, y) sampled from all columns. In general, sampling points can be arbitrarily chosen from available mesh points. Here, to reduce computational cost and better leverage the sequential dependency, we choose the sampling points to be the mid-axis, rightmost mesh point of each column.

Construction and training of the RNN model

Our network architecture consists of a sequence input layer, an LSTM [START_REF] Hochreiter | Long Short-Term Memory[END_REF] Two ML models are trained to predict coordinate x and y separately to better identify their respective errors. Let Fx and Fy denote the built models for x and y, respectively, we have ,

where x and y are predicted coordinate vectors whose components, xi and yi, i ∈ [1,Nx], are coordinates of sampling respectively. M is the digital encoding of a material distribution (an Ny×Nx array composed of "1" and "0"). θx and θy denote the learnable parameter sets for Fx and Fy, respectively. The loss function is defined as the half-meansquared-error between the predicted coordinates (xi or yi) and true coordinates (xi true or yi true ) for a specific material distribution, i.e., ,

for coordinate x and y, respectively. The training of ML models is to find the optimal θx and θy that minimize the loss values on the training set, i.e., and .

The initial learning rate is set to 0.005, which decreases by multiplying a factor of 1 √2 ⁄ every 50 epochs. The training stops after the validation loss converges. The mini-batch size during training is set to 64 or 10% of the training set size, whichever is smaller. The adaptive moment estimation (Adam) [START_REF] Kingma | A method for stochastic optimization[END_REF] optimizer is used to train the network.

Moreover, to fully exploit the parallel computing capability of ML and improve the optimization speed, the shape evaluation process in ML-SSO is implemented in a vectorized form. In addition, we examine the time cost for ML prediction of massive (5000) data with different batch sizes, and the optimal batch size is found to be 512. Note that in ML-SSO, the number of all potential designs for each subdomain is 2^(Ny×Nsub). We thus choose the batch size to be 512 for designs with Nsub=3 and to be 2^(Ny×Nsub) for Nsub=1 ( ) ( )

; , ;

x x y y F F = = x θ M y θ M ( ) ( ) ( ) ( )

Method of gradient calculation for Figure 1d

The gradients of actuated coordinates x or y with respect to the voxel encoding can be written as [START_REF] Yue | Cold-programmed shape-morphing structures based on grayscale digital light processing 4D printing[END_REF] Since x and y are functions of M as described by (3), which involves differentiable operations only, we use automatic differentiation to evaluate the gradients in Eq.( 5) through the dlgradient function in Matlab. As defined in Section 4.2, M is an array composed of "1" and "0". Here the component of M is denoted by Mv, which represents the material encoding of voxel v, where v∈ [START_REF] Ge | Active materials by four-dimension printing[END_REF]Ny×Nx] represents the voxel number. The gradient values, ∂xi/∂Mv and ∂yi/∂Mv, for i=12 and 24, based on an all-passive state (Mv=0 for any v∈ [START_REF] Ge | Active materials by four-dimension printing[END_REF]Ny×Nx]), are then displayed in Figure 1d.

Materials and 4D printing

The photocurable resin is prepared by mixing isobornyl acrylate (IOBA, Sigma-Aldrich) and aliphatic urethane diacrylate (AUD, Ebecryl 8402, Allnex, GA, USA) in a weight ratio of 1:1. Then, 1 wt% photoinitiator (Irgacure 819, Sigma-Aldrich), 0.08 wt% photo absorber (Sudan I, Sigma-Aldrich), and 0.04 wt% fluorescent dye (Solvent green 5, Orichem International Ltd., Hangzhou, Zhejiang, China) are added. The resin is thoroughly mixed before printing.

The grayscale digital light processing (g-DLP) printing technique [START_REF] Yue | Single-vat single-cure grayscale digital light processing 3D printing of materials with large property difference and high stretchability[END_REF] is used to print the designed structure, where the degree of curing (DoC) can be locally controlled by the or ¶ ¶ ¶ ¶ x M y M assigned light intensity. Our ML-SSO designs are transformed into grayscale printing slices, where the active ("1") and passive ("0") phases correspond to grayscale percentages of 0% (hence higher light intensity) and 60% (hence lower light intensity), respectively, which later spatially assign the high-DoC and low-DoC phases in the printed structure, respectively. The printed structure is then placed in an 80℃ oven for 8 hours to facilitate monomer volatilization. The low-DoC phase contains more residual monomers that can volatize and thus shows more volume shrinkage than the high-DoC phase at elevated temperatures. The shrinkage strain mismatch of the two phases induces the shape transformation.

Modification of optimal designs

Note that the material properties in experiments are different from those used in FE simulations. Experimental characterizations show that the printed two material phases show a modulus ratio of 0.06, while the ML-EA design assumes the identical modulus for two constituent phases. Additionally, the practical expansion mismatch is identified to be 0.05, which is also different from that used in ML-EA (i.e., 0.1). Such issues can be resolved by retraining the ML model based on the FE data with practical material properties (expansion mismatch and modulus difference) and re-running the ML-EA. Here, instead of retraining a new model, we adopt a design conversion strategy [START_REF] Sun | Machine Learning-Evolutionary Algorithm Enabled Design for 4D-Printed Active Composite Structures[END_REF] to approximately compensate effects of property difference of the two phases on the shape change, i.e., the optimal designs are converted using the analytical curvatures of multi-layer composite beams. The effectiveness of this strategy has been validated by FE simulations and experiments, as detailed in our previous work [START_REF] Sun | Machine Learning-Evolutionary Algorithm Enabled Design for 4D-Printed Active Composite Structures[END_REF].

Figure 1 .

 1 Figure 1. Schematic illustration of the active composite (AC) beam and the sequential dependency. (a) Actuation of an AC beam with voxel-level material distribution. The two materials, encoded as "1" (active) and "0" (passive), have a volume expansion mismatch. (b) Architecture of an RNN, which shows a similar structure with the AC beam. (c) Density scatter plots of the ground-truth versus ML-predicted values of the coordinate x (left) and y (right). The color indicates the relative point density. (d) Gradients (or sensitivities) of actuated coordinates (x, y) of the 12th and 24th voxel columns with respect to all voxels (i.e., ∂xi/∂Mv and ∂yi/∂Mv for i=12 and 24, where v=1, 2, …, 4×24 represents the voxel

Figure 2 .

 2 Figure 2. Demonstration of ML-SSO on inverse design problems. (a) Schematic of ML-SSO. (b) One-period target and ML-SSO optimized shape with Nsub=1. (c) Two-period target and ML-SSO optimized shape with Nsub=1. (d) Three-period target, ML-SSO optimized shape with different Nsub values, and illustration of optimization sequence with Nsub=2. (e) Half-butterfly target, ML-SSO optimized shape with different Nsub values, and illustration of optimization sequence with Nsub=2. All ML-SSO designs are based on Nx=24.

Figure 3 .

 3 Figure 3. (a) RMSEs of optimized shapes by ML-SSO (with different Nsub values) and ML-EA for the four target shapes considered in Figure 2. The RMSE is calculated using Eq.(2). (b) Time cost of the corresponding optimizations, which are run using the Intel Core i9-10900 CPU and NVIDIA Quadro P620 GPU.

( 2 )

 2 between the optimized coordinates (xi, yi) and target coordinates (𝑥 # & , 𝑦 # & ) for all Nx sampling points. As shown in Figure3, for the relatively simple one-period and two-period shapes, ML-SSO with Nsub=1 which only takes 0.23 s achieves similar RMSEs with that of ML-

Figure 4 .

 4 Figure 4. Ultrafast ML-SSO enabled streamlined design and fabrication for 4D printing based on hand-drawn lines. From top to bottom: "G", "T", half-butterfly, mountain. From left to right: hand-drawn lines, optimized shapes by ML-SSO (with different Nsub values) and ML-EA, the corresponding performance in terms of RMSE and time cost, the optimal ML-SSO designs selected (with greater Nsub, and RMSE < 1 mm) out of those in column 3 (arrows), and 4D-printed, actuated shapes. All scale bars: 10 mm.

  forward and inverse problems for beam structures with varying lengths or Nx values. Here, we fix the voxel number Ny=4 in thickness and the voxel size. Thus, the sequential (column-by-column) prediction ability of the RNN-ML allows it to be directly applied to structures with different Nx values and hence different lengths. For the case Nx>24, we propose an ML-based splicing method, with no need for retraining the model, for forward predictions with different Nx values.

Figure 5 .

 5 Figure 5. Demonstration of ML-splicing on forward prediction problems of beam structures with different Nx values. (a) Schematic of ML-splicing. (b) Prediction results for benchmark input sequences of different Nx values.

Figure 6 .

 6 Figure 6. Demonstration of the splicing-empowered ML-SSO on inverse design problems using different Nx (numbers lengthwise voxels) or numbers of sampling points. Multiple target shapes are considered, including (a) three-period shape, (b) hand-drawn mountain shape, (c) six-period shape, and (d) "Chinese dragon" profile. Image by brgfx on Freepik (https://www.freepik.com/free-vector/chinese-dragon-flying-clouds-isolated-white-

Figure 7 .

 7 Figure 7. Applications of ML-SSO to rapid design for active lattice structures. (a) Arbitrary

and 2 .

 2 These strategies exploit the ML capability to deliver ultrafast and massive predictions.

  are the raw and normalized feature values, respectively, and mean is the mean value and std is the standard deviation. The randomly generated raw

	inputs (numerous '1' and '0') show a mean value of 0.5 and a standard deviation of 0.5. As
	a result, the input state '0' and '1' become '-1' and '1' after normalization. Such normalization
	is found to improve the network performance. The hidden size (number of neurons of each

layer, a fully connected layer, and a regression layer. The LSTM is a special type of RNN that addresses the issues of vanishing or exploding gradients presented in long sequences. The implementation, training, and testing are conducted using Matlab (2020a, MathWorks, Natick, MA). Before the training, all the input and output data are normalized through x' = (x-mean(x))/std (x), where x and x' neural layer) of LSTM is set as 50, and the LSTM layer has 200×(50+Ny+1) learnable parameters. LSTM can process a sequence input of any lengths, and different time steps of LSTM share the same learnable parameters.

Acknowledgments H.J.Q. acknowledges the support of an AFOSR grant (FA9550-20-1-0306; Dr. B.-L. "Les"