Physics-Guided Surrogate Modeling for Cross-Architecture Normalized Power-Density Optimization in Perovskite Solar Cells

2.2 Physics and equations

The governing equations reported below, Eqs. (1)-(9), correspond to the semiconductor drift–diffusion framework implemented in COMSOL and are used to describe the electrical response of the PSC structures.

$\nabla \cdot(\varepsilon \nabla \psi)=-q\left(p-n+N_D^{+}-N_A^{-}\right)$                          (1)

$J_n=q \mu_n n E+q D_n \nabla_n$                          (2)

$J_p=q \mu_q p E-q D_p \nabla_p$                          (3)

$D=\frac{\mu k_B T}{q}$                          (4)

$E=-\nabla \Phi$                          (5)

$P=J \cdot E$                          (6)

$\eta=P_{\text {out }} / P_{\text {in }}$                          (7)

$\frac{\partial n}{\partial t}=\frac{1}{q} \nabla J_n-R+G$                          (8)

$\frac{\partial p}{\partial t}=\frac{1}{q} \nabla J_p-R+G$                          (9)

where, $J_n$ and $J_p$ are the electron and hole current densities, $n$ and $p$ are the electron and hole concentrations, $D_n$ and $D_p$ are the carrier diffusivities, $\mu_n$ and $\mu_p$ are the electron and hole mobilities, $k_B$ is the Boltzmann constant, $T$ is temperature, $q$ is the elementary charge, $\phi$ is the electrostatic potential, $\varepsilon$ is the permittivity, and $N_D^{+}$ and $N_A^{-}$ are the ionized donor and acceptor concentrations.

The equations reported above are used to describe the electrical response of the investigated PSC structures. They are not intended to represent a full wavelength-resolved optical-generation model. The present study therefore focuses on the influence of material transport parameters and layer thickness on the simulated current-density response, rather than on experimentally measured power conversion efficiency.

Accordingly, the reported target quantity denoted as $D_{\max}$. represents a simulated maximum normalized electrical powerdensity indicator. It is obtained from the bias-dependent current density calculated by the semiconductor model and the applied voltage according to the electrical power-density relation $P=V_{\text {app}}.J$. The resulting power-density response is normalized with respect to the standard one-sun irradiance of $1000 \mathrm{~W} /\mathrm{m}^2$, and $D_{\text {max}}$ is extracted as the maximum value over the applied voltage sweep. Thus, $D_{\text {max}}$ is interpreted as a simulation-based normalized power-density metric, not as an experimentally measured solar-cell power conversion efficiency.

The simulation employs the Semiconductor Module in COMSOL Multiphysics with a stationary study to investigate the bias-dependent electrical response of the layered perovskite solar-cell structures. Each device layer is modeled as a semiconductor domain using material parameters including relative permittivity, band gap, electron affinity, effective density of states, carrier mobility, doping concentration, and Shockley–Read–Hall carrier lifetime. The ETL is treated as an n-type layer with an assigned donor concentration, while the HTL is treated as a p-type layer with an assigned acceptor concentration. Non-ideal charge-carrier losses are represented through recombination models available in the Semiconductor interface, including Shockley–Read–Hall recombination and carrier-lifetime-based trap-assisted recombination.

A swept mesh is applied to resolve the layered structure, particularly near the ETL/perovskite and perovskite/HTL interfaces. A parametric voltage sweep from 0 to 1.2 V is conducted by applying the bias to the gold contact, while the bottom electrode is used as the reference contact. For each thickness combination in the 50–500 nm range, the stationary drift–diffusion problem is solved using the Newton–Raphson method. The simulated current density is then used to calculate the normalized power-density response, from which $D_{\max}$ is extracted.

For reproducibility, we clarify that a wavelength-resolved optical model based on material-specific complex refractive indices $n(\lambda)$ and $k(\lambda)$ was not included in the present implementation. Therefore, the study is regarded as a semiconductor-transport-based surrogate modeling framework for normalized electrical power-density optimization. Incorporating wavelength-dependent optical generation and experimentally validated optical constants is identified as a future extension of the framework.

The top contact, modeled as an ideal ohmic interface, consists of gold (Au), chosen for its high conductivity, chemical stability, and suitable work function alignment with the HTL, enhancing hole collection and minimizing energy loss. The bottom contact is assigned to fluorine-doped tin oxide (FTO), commonly used as a transparent conductive electrode in PSC structures and serves here as the bottom electrical contact.

The inclusion of two distinct material systems is critical for developing a unified machine-learning surrogate that captures relationships between material properties, layer geometry, and device performance across the investigated architectures. While a model trained on a single device architecture may capture specific trends, using electronically diverse systems allows for deeper learning across varying material classes and interface behaviors. The selected systems display significant differences in key optoelectronic parameters: TiO₂ has relatively low electron mobility (≈0.1 cm²/V·s) compared to ZnO (≈200 cm²/V.s); MAPbI₃ exhibits a direct bandgap of about 1.55 eV, while CsFAPbI₃ has a slightly narrower gap of ≈1.48 eV; Spiro-OMeTAD presents a hole mobility on the order of $2 \times 10^{-5}$ cm²/V·s, whereas PTAA shows ≈ $4 \times 10^{-5}$ cm²/V·s. These contrasts ensure that the ML model learns to correlate D_max not only with geometry but also with intrinsic material parameters such as mobility, energy alignment, and carrier lifetime. In addition, the different band alignments, TiO₂ with an electron affinity of ≈4.1 eV versus ZnO at ≈4.5 eV, and the contrasting hole extraction properties of Spiro-OMeTAD and PTAA allow the model to capture how interfacial effects influence performance. Training and validating on both systems reduce the likelihood of overfitting to a single architecture and enables the model to learn predictive relationships across the investigated dual-architecture dataset. Since both material systems are represented during training, the reported evaluation demonstrates predictive performance on previously unseen thickness combinations within the investigated design space rather than extrapolation to entirely unseen material chemistries. Although the primary focus of this work is on data-driven D_max prediction rather than detailed device physics analysis, a representative simulation output is included to illustrate the numerical modeling framework used for data generation. Figure 2 shows the electron concentration distribution across a ZnO (ETL)/CsFAPbI₃ (PVK)/PTAA (HTL) structure as obtained from the COMSOL Multiphysics semiconductor model.

The three-dimensional visualization was generated by revolving the two-dimensional axisymmetric semiconductor model about the symmetry axis. A full 360° rotation was applied using the COMSOL Multiphysics Revolution 2D plotting feature to illustrate the spatial distribution of electron concentration within the device structure. This example does not constitute a performance result but serves to demonstrate the physical consistency of the device structure, layer definitions, and carrier transport implementation, thereby illustrating the simulation environment from which the bias-dependent current-density response and $D_{\max}$ dataset were obtained.

In this figure, we can notice that the electron density is highest within the ETL, reflecting its role as an electron-selective contact that facilitates efficient carrier extraction, while it is strongly suppressed in the HTL, confirming effective electron blocking at the hole-collecting side. Within the perovskite absorber, the gradual variation in electron concentration indicates the interplay of carrier transport, drift–diffusion behavior, and recombination across the device thickness. The pronounced contrast in carrier density at the ETL/PVK and PVK/HTL interfaces demonstrates proper band alignment and junction formation, validating the physical consistency of the Multiphysics model used to generate simulated electrical-response data for the subsequent data-driven D_max analysis.

The perovskite solar-cell structure was implemented using a two-dimensional axisymmetric formulation in COMSOL Multiphysics. Although the investigated device stack is physically planar, the axisymmetric representation provides a computationally efficient framework for solving the semiconductor drift–diffusion equations across the layered structure while preserving the one-dimensional nature of carrier transport along the device thickness. Since the material properties, doping profiles, and electrical boundary conditions are defined uniformly in the lateral direction, the model response is governed primarily by the vertical transport across the ETL/perovskite/HTL stack. The axisymmetric formulation was therefore used as a practical numerical representation of the layered device and to enable three-dimensional visualization through the COMSOL Revolution 2D plotting feature shown in Figure 2. No lateral device effects are considered in the present model.

The next section examines the application of synthetic data within ML frameworks, outlining the procedures for data preparation, model assessment, and algorithm training used to predict the D_max normalized power-density indicator from structural parameters.

Mathematically, this smoothing bias induced by the L₂ penalty and dropout preferentially preserves lower-frequency components of the response surface, trading marginal high-frequency accuracy for improved generalization to unseen data, a fundamental and desirable compromise that enhances robustness against simulation noise and local numerical variability. We partitioned the validation dataset by material system and D_max regime, with the detailed performance metrics presented in Table 3. A clear D_max -dependence pattern emerged, as error magnitude correlates positively with the D_max level—a result anticipated because high-D_max devices correspond to optimal parameter combinations where the response surface exhibits steep gradients, and even small perturbations produce substantial changes in output. In contrast, low-D_max regimes display gentler gradients, yielding smaller absolute errors despite comparable relative errors. Cross-configuration analysis reveals that the ZnO/PTAA system exhibits approximately 1.8 times higher error than the TiO₂/SPIRO system, aligning with the hypothesis that material-specific recombination complexity introduces additional variance. Notably, when data are binned by D_max regime, the high-D_max domain shows elevated errors for both configurations, indicating that the response surface possesses genuine curvature in the high-performance region irrespective of the material system.

4.3 Feature importance and physical interpretability via permutation and SHAP

Feature importance was quantified by evaluating the reduction in R² after randomly permuting each retained feature; the results are summarized in Table 4. A clear hierarchy emerges, with total thickness showing the highest relative importance (87.3%). This indicates that the surrogate relies strongly on layer-scaling information within the investigated design space. However, this relationship should be interpreted as a learned correlation from the COMSOL-generated dataset rather than as direct evidence of wavelength-resolved optical absorption or a verified causal mechanism. In particular, the association between increased thickness and higher predicted D_max should be treated cautiously, since thicker semiconducting layers may also increase carrier transport length, recombination probability, and bulk series resistance. Thus, total thickness is best understood here as a dominant geometric descriptor that captures combined layer-scaling, transport, and recombination-related effects in the simulated drift–diffusion response.

Table 4. Permutation-based feature importance

Aspect ratio exhibits secondary importance (53.9%), modulating the trade‑off between series resistance and interface recombination. Higher ratios correspond to thicker perovskite layers relative to the transport‑layer stack, which reduces bulk resistive losses but increases recombination at interfaces. Its moderate importance suggests this balance is a material‑independent, yet subordinate, design consideration. The material‑specific interaction term (Config‑PVK interaction) accounts for 23.1% of relative importance, indicating that the model successfully learned configuration‑dependent sensitivities—such as the influence of ZnO’s superior mobility on optimal thickness ratios—while confirming that thickness‑scaling laws retain a universal component.

To elucidate how individual features influence specific predictions, SHAP (SHapley Additive exPlanations [38]) analysis was applied to a sample of 50 test instances. SHAP decomposes each prediction into the sum of marginal contributions from each feature, quantifying deviations from the baseline expected output (Figure 5).

Key insights include: total thickness displays a monotonic positive impact across its range (150–1400 nm) with a diminishing slope of approximately 0.012 normalized % per nm, consistent with the sub-linear saturation expected from thickness-dependent transport and recombination trade-offs; aspect ratio shows a non-monotonic relationship, with negative SHAP values at low ratios (<1.5) and positive contributions at high ratios (>2.5), revealing an optimal regime near Γ ≈ 2.0 that balances series resistance against interface recombination; and the Config-PVK_interaction exhibits a switch-like behavior, wherein ZnO configurations produce systematically lower SHAP values than TiO₂ at equivalent perovskite fractions, capturing a configuration-dependent response in which ZnO-related transport properties modify the predicted sensitivity of D_max to perovskite thickness.

4.4 Training dynamics and convergence behavior

The model demonstrated rapid and stable convergence, reaching a validation‑loss plateau within 10 epochs (Figure 6). Initial losses were high (training MSE = 85.3, validation MSE = 78.4 at epoch 1), but by epoch 10 they had dropped to 0.62 MSE and 0.58 MSE, respectively, and ultimately settled at 0.51 MSE (training) and 0.52 MSE (validation) by epoch 38. The near‑perfect alignment between training and validation curves throughout training indicates that the combined regularization (dropout + L₂) effectively prevented overfitting, that the VIF‑selected feature set was appropriately balanced—neither excessively rich nor overly sparse—and that the model architecture (a five‑layer network with a maximum of 256 neurons) was well‑matched to the dataset size (1,600 training samples) [33]. Gradient‑flow analysis further confirmed optimization health: the gradient norm in the input layer stabilized around 0.008 and in the output layer around 0.042 within five epochs, yielding a stable output‑to‑input gradient ratio of approximately 5.25, well within the ideal range of 1–100. The complete absence of exploding (> 1) or vanishing (< 10⁻⁶) gradients verifies that the initialization and batch normalization successfully maintained stable signal propagation throughout the deep network.

Figure 7. 3D optimization point cloud derived from the surrogate model

Table 6 compares the surrogate-predicted D_max values with the corresponding COMSOL-simulated values. For the TiO₂/MAPbI₃/Spiro-OMeTAD architecture, strong agreement was obtained across all evaluated configurations. The absolute prediction error ranged from 0.36 to 0.55 normalized, corresponding to relative errors between 2.08% and 3.15%. This confirms that the surrogate accurately reproduces the local response surface near the TiO₂-based candidate optimum.

For the ZnO/CsFAPbI₃/PTAA architecture, larger discrepancies were observed near the surrogate-predicted optimum. Absolute errors ranged from 1.72 to 3.65 normalized %, corresponding to relative errors between 11.70% and 24.78%. In all evaluated ZnO-based cases, the surrogate overestimated D_max relative to the COMSOL simulations. Therefore, the ZnO-predicted optimum should be interpreted as a candidate region requiring further refinement rather than as a confirmed optimum.

This behavior is consistent with a known limitation of surrogate-based optimization. Although the surrogate can provide accurate predictions across much of the investigated design space, optimization algorithms naturally search for regions of maximum predicted response, where small local approximation errors may be amplified. Consequently, the surrogate is useful for rapidly identifying promising candidate designs, but final optimization should be confirmed using high-fidelity physics-based simulations before experimental implementation.

Overall, the additional COMSOL validation demonstrates architecture-dependent surrogate reliability. The TiO₂/MAPbI₃/Spiro-OMeTAD candidate region shows strong surrogate–COMSOL consistency, whereas the ZnO/CsFAPbI₃/PTAA candidate region requires denser local sampling and further validation. These findings support the use of the surrogate model as a computational screening tool while emphasizing the importance of physics-based verification for surrogate-selected optima.

4.9 Limitations and scope for further validation

Several aspects of this study warrant explicit acknowledgment as current limitations, each identifying a concrete direction for follow‑up validation. Although grid-proximal virtual screening reduces the risk of off-grid extrapolation, targeted COMSOL validation showed that surrogate accuracy is not uniform across both architectures. The TiO₂/MAPbI₃/Spiro-OMeTAD candidate region exhibited relatively close agreement between simulated and predicted D_max, whereas the ZnO/CsFAPbI₃/PTAA candidate region showed larger local deviations. This indicates that the ZnO-predicted optimum should be interpreted as a candidate region requiring further refinement rather than a confirmed optimum. Future work should include denser COMSOL sampling around high-D_max  regions, uncertainty-aware surrogate modeling, and validation against illuminated J – V data. First, the reported R² and MAPE reflect performance on a random 80/20 hold‑out split; because material‑system identity (TiO₂/MAPbI₃/Spiro‑OMeTAD vs. ZnO/CsFAPbI₃/PTAA) is itself an input feature, this split evaluates interpolation within the combined dual‑architecture design space but does not by itself establish generalization to unseen chemistries. A leave‑one‑architecture‑out evaluation—training exclusively on one system and testing on the other—is required before cross‑architecture transferability can be claimed with confidence and has not yet been performed. Second, the VIF‑based feature selection used for model training retained only total thickness, aspect ratio, and the configuration–PVK interaction term; these three descriptors cannot uniquely resolve individual ETL and HTL thicknesses in the reduced feature space. The optimum thickness combinations reported in Section 4.5 (TiO₂/MAPbI₃/Spiro‑OMeTAD: ETL = 499 nm, HTL = 494 nm, PVK = 500 nm; ZnO/CsFAPbI₃/PTAA: ETL = 50 nm, HTL = 50 nm, PVK = 496 nm) should therefore be interpreted as geometry recommendations within the investigated design space rather than as proofs of uniquely optimized ETL and HTL thicknesses beyond this range. Third, although the one‑million‑point screening shows that surrogate predictions vary smoothly with distance from the original 50 nm thickness grid and the reported optima lie within approximately 4–6 nm of simulated thickness combinations, the near‑grid region has only recently been probed by targeted COMSOL drift–diffusion simulations at the corrected optima. These checks indicate architecture-dependent agreement: the TiO₂/MAPbI₃/Spiro-OMeTAD candidate region shows relatively close surrogate–COMSOL consistency, whereas the ZnO/CsFAPbI₃/PTAA candidate region shows larger local deviations. Extending validation to broader and denser off-grid thickness combinations remain an important follow-up task. Fourth, the DNN (~45,000 parameters trained on ~1,600 samples with 3 retained features) has now been benchmarked against simpler alternatives, including linear regression, random forest, and a compact multilayer perceptron, yet further work is needed to explore more compact architectures and uncertainty‑aware surrogates that can provide prediction intervals or calibrated error estimates. Fifth, model validation to date relies primarily on spatial carrier‑density profiles (e.g., the electron‑concentration map in Figure 2) rather than full current–voltage (J–V) characteristics; direct comparison of simulated J–V curves and conventional illuminated photovoltaic metrics such as Voc, Jsc, fill factor, and power‑conversion efficiency with representative experimental or literature data has not yet been carried out. Addressing these limitations will require additional COMSOL simulations and machine-learning experiments beyond the scope of the present study and is therefore identified as a priority for future work.