ABSTRACT Structural dynamic response time histories required for seismic analysis are particularly susceptible to continuous data missing occurring simultaneously across all sensor channels during seismic hazard. Earthquake‐induced disruptions often induce power interruptions and communication failures, and thus may delay analytical procedures or render structural evaluation infeasible. Although recent deep learning imputation methods can capture spatiotemporal correlations, they require some completely observed measurements and therefore fail to recover fully missing time intervals across all channels. Conversely, matrix completion methods exploit intrinsic low‐rank characteristics of structural dynamic responses. However, the temporal evolution of structural dynamics is not explicitly encoded in their formulations, and convex relaxations tend to introduce over‐shrinkage and estimation bias. Therefore, this study proposes a smoothly clipped absolute deviation (SCAD) low‐rank informed generative adversarial imputation network with an enhanced temporal convolutional generator for the challenging problem of continuous and simultaneous missing data imputation in structural dynamic responses. The physics‐informed low‐rank constraint is embedded into the generative process to capture global inter‐sensor correlation structures. Meanwhile, the temporal convolutional network (TCN) generator explicitly models long‐range temporal dependencies, which supports the reconstruction of nonlinear vibration characteristics across extended missing intervals. Adversarial learning constrains the imputed responses to be statistically consistent with the observed data. The proposed method is validated through two real‐world case studies of a bridge and the Canton Tower under ambient vibration or in‐situ seismic excitation. Ablation studies are conducted to examine the individual contributions of the enhanced TCN generator and the SCAD low‐rank constraint. In addition, comparative studies with state‐of‐the‐art approaches demonstrate the reconstruction accuracy and computational efficiency, highlighting its practical potential for rapid post‐earthquake structural assessment.

Continuum damage mechanics is characterized by mesh-dependent results unless specific countermeasures are taken. The most popular remedies involve introducing either nonlocality via filtering or a gradient extension for the damage variable(s). Such approaches have their limitations, e.g., they are hard to integrate into conventional finite-element codes, involve parameters that are non-trivial to determine experimentally and are incompatible with a scale transition that is both physically and mathematically sensible. The work at hand considers an alternative route to obtain mesh-independent damage models, namely via convex relaxation. Such convex damage models were considered before, but they are usually not capable of representing softening behavior. Schwarz et al. (Continuum Mech. Thermodyn., 33, pp. 69–95, 2021) proposed such a strategy by considering the convex envelope of a rate-limited simple damage model, i.e., an isotropic damage model without tension-compression anisotropy at small strains. However, they were not able to compute the envelope explicitly and provided an approximation only. In the work at hand, we introduce a number of conditions on the damage-degradation function which permit us to compute the convex envelope analytically for a large class of damage-degradation functions used in small-strain isotropic damage models. Interestingly, the obtained models involve a one-dimensional damaged microstructure, i.e., damage distributions emerge naturally. The resulting model is structurally simple and purely local, i.e., gradient-free, thermodynamically consistent and readily integrated into standard finite-element codes via traditional user subroutines. We discuss the computational and solid mechanical aspects of the ensuing model and demonstrate its numerical robustness via dedicated computational experiments. We also show that the model permits to be homogenized by considering a representative volume element study for an industrial-scale fiber-reinforced composite.

The joint design of unimodular waveforms and receive filters is a pivotal technology in Multiple-Input Multiple-Output (MIMO) radar systems. However, most existing methods primarily focus on point target detection or ignore the impact of active jamming in extended target scenarios. To bridge this gap, this paper proposes an optimization framework for the joint design of unimodular waveforms and receive filters specifically for MIMO radar extended target detection in the presence of suppressive jamming. The problem is formulated to maximize the Signal-to-Interference-plus-Noise Ratio (SINR) while strictly satisfying the unimodular constraint and mitigating suppressive jamming. Due to the non-convexity of the unimodular constraint and the quadratic fractional nature of the SINR objective function, the optimization problem is highly challenging. Unlike conventional methods that rely on convex relaxation—which often leads to performance degradation—we exploit the geometric structure of the constraint set. Specifically, the unimodular constraints are modeled using complex circle manifolds, and the suppressive jamming suppression requirements are integrated into the objective function via a smooth penalty metric. Building on these characteristics, a Product Complex Circle Euclidean Manifold (PCCEM) method is developed. This approach transforms the constrained problem into an unconstrained optimization task on a product manifold, which is then efficiently solved using the limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) algorithm. Simulation results demonstrate that the proposed PCCEM method outperforms baseline algorithms in terms of computational efficiency, output SINR, and the depth of the formed jamming notches.

With the high penetration of distributed energy resources (DERs), which are characterized by stochasticity and intermittency, traditional centralized optimization methods face challenges such as communication packet loss, low reliability, and poor scalability in large-scale DC microgrids. Therefore, distributed optimization methods have attracted attention due to their robustness and scalability. This paper extends our previous conference work by proposing a convex-relaxation-based distributed control strategy for DC microgrids with constant power loads (CPLs) and maximum power point tracking (MPPT)-controlled distributed generations (MPPT-DGs). Furthermore, a control strategy based on distributed observers is designed to achieve global optimal control under sparse communication networks. First, an exact convex relaxation method is applied to transform the original non-convex optimal power flow (OPF) problem into a convex problem, with theoretical guarantees of exactness. Then, the Karush–Kuhn–Tucker (KKT) conditions are equivalently transformed into a consensus-based optimality condition and integrated into the distributed control framework. Next, small-signal stability analysis is performed to verify the system’s robustness. To reduce communication costs, a distributed observer-based control strategy is proposed, which can achieve optimal control under sparse communication networks. The impact of communication delays on system stability is also investigated. Finally, the simulation results verify the accuracy of convex relaxation, the effectiveness of the proposed control strategy, and its performance under communication delay.

The FLASH effect can significantly reduce radiation-induced normal tissue damage while maintaining tumour control, but requires ultra-high dose rates and high doses. This work proposes a single-field-uniform-dose-per-fraction simultaneous dose and dose rate optimization (SFUDPF-SDDRO) method for proton FLASH radiotherapy to ensure both dose rate and dose meet FLASH effect thresholds. The SFUDPF method focuses on delivering the prescription dose for each fraction from only a single field instead of multiple fields, which inherently supports the ultra-high dose rate and high dose necessary for the FLASH effect. We performed retrospective FLASH treatment planning utilizing SFUDPF-SDDRO on four clinical head-and-neck (HN) cases for this study. SFUDPF planning involves delivering each prescription fraction (8Gyx5 fx) in 1 beam angle as opposed to multiple beam angles per fraction for IMPT. For each beam delivery, we maximized the FLASH effect in a 1cm expansion of the HN CTV (CTV+1cm) by enforcing FLASH dose-rate and dose thresholds of 40Gy/s and 5Gy, respectively, in this region. The pencil-beam-scanning dose rate (PBSDR) was calculated voxel-wise by modeling the raster-scanning spot trajectory, while neglecting energy switching times under the assumption of a range modulator capable of expanding a single-energy beam into a spread-out Bragg peak (SOBP). Robust optimization at 3mm/3.5% was performed to address setup and range uncertainties. We employed iterative convex relaxation and alternating direction method of multipliers algorithms to solve the non-convex optimization problem posed by the SFUDPF-SDDRO model. The FLASH effect was modelled within this work by multiplying the proton dose with a constant 0.7 dose modification factor for voxels fulfilling the dose-rate and dose thresholds to obtain the FLASH effective dose (FED). Effects of FLASH sparing maximization via SFUDPF-SDDRO are verified by comparing with IMPT and VMAT on plan qualities such as (i) high-dose area sparing, (ii) conformity index (CI), and (iii) OAR doses. FLASH RT via SFUDPF-SDDRO compared with IMPT and VMAT was evaluated for four clinical HN cases with different tumor geometries. When compared with their VMAT counterparts, SFUD-SDDRO achieved a considerable reduction of FED for OAR directly adjacent to the CTV. Specifically in case 1, the brainstem D1% decreased from 87.57% to 62.26%, and the spinal cord D10% decreased from 87.36% to 60.74%; in case 2, the D10% of the carotid decreased from 102.46% to 63.30%; in case 3, the D10%of the oral cavity decreased from 94.72% to 62.66%, and the D10% of the oropharynx decreased from 102.5% to 69.09%; in case 4, the D10% of the oral cavity decreased from 88.56% to 59.81%. The SFUDPF-SDDRO achieved a satisfactory CI in terms of FED, indicating that conformity was not sacrificed to achieve the FLASH effect. The proposed SFUDPF-SDDRO method is feasible and shows potential clinical benefits for FLASH treatment planning. Maximizing the FLASH effect within a 1cm ring around the target substantially limits high-dose spillage and enhances OAR sparing compared with conventional approaches.

The Area Under the Receiver Operating Characteristics Curve (AUC) is a widely used metric for evaluating model performance across all possible decision thresholds. Existing methods for AUC optimization typically assume a predefined parametric distribution of thresholds. However, the optimal decision threshold depends on the misclassification costs, which follow a non-parametric distribution.This motivates us to introduce a variant of AUC, termed Cost-aware AUC (CAUC), where the thresholds are conditioned on an empirically determined cost distribution. Unfortunately, as a bilevel problem, it is challenging to directly optimize the CAUC: 1) The inner problem of finding the optimal thresholds is non-convex, leading to potential issues with convergence; 2) The outer problem involves the derivative of False Positive Rate (FPR) w.r.t. the threshold, which is unavailable without an explicit formulation of threshold distribution. To address challenge 1), we utilize the convex relaxation technique to reshape the inner problem into a convex one. Facing challenge 2), we propose an adaptive kernel density estimation framework. Specifically, the derivative of FPR is considered an aggregation of various kernel functions. To avoid manually crafting the aggregation function, we propose a finite-difference-based stochastic algorithm to optimize the model without explicit aggregation function. Theoretically, the proposed algorithm enjoys a convergence rate of $\mathcal {O}(\epsilon ^{-4})$O(ε-4). Empirical studies across various datasets and cost distributions speak to the effectiveness and soundness of our framework.

Distribution network optimization is represented by non-convex power flow equations, where traditional convex relaxation methods may lead to inaccurate or infeasible solutions. To effectively solve the non-convex optimization problem, an eigenvalue-based hidden convexity (EBHC) algorithm is proposed for radial distribution network optimization with continuous variables and separable objectives. Specifically, the original nonconvex problem is reformulated within a standard ADMM framework by separating the convex and non-convex components of the quadratic equality constraints. The convex component is represented as a centralized second-order cone programming (SOCP) problem that can be efficiently solved, while the non-convex component is decomposed into multiple sub-problems, each taking the form of a quadratic program with a single quadratic constraint. Each sub-problem exhibits a special mathematical structure, enabling parallel solution via multiple generalized eigenvalue problems (GEPs). The theoretical guarantees of the algorithm rely on the assumption that Slater's condition holds and that the network is radial. Through iterative coordination between the convex and non-convex components, the non-convex quadratic equality constraints are strictly satisfied. Numerical experiment results demonstrate the effectiveness and accuracy of the proposed algorithm.

Our ability to structure materials at the nanoscale has, and continues to, enable key advances in optical control. In pursuit of optimal photonic designs, substantial progress has been made on two complementary fronts: bottom-up structural optimizations (inverse design) discover complex high-performing structures but offer no guarantees of optimality; top-down field optimizations (convex relaxations) reveal fundamental performance limits but offer no guarantees that structures meeting the limits exist. We bridge the gap between these two parallel paradigms by introducing a "verlan" initialization method that exploits the encoded local and global wave information in duality-based convex relaxations to guide inverse design towards better-performing structures. We first illustrate this technique via the challenging problem of Purcell enhancement, maximizing the power extracted from a small emitter in the vicinity of a photonic structure, where ill-conditioning and the presence of competing local maxima lead to sub-optimal designs for adjoint optimization. Structures discovered by our verlan method outperform standard (random) initializations by close to an order of magnitude and approach fundamental performance limits within a factor of two, highlighting the possibility of accessing significant untapped performance improvements. We further validate this method using a planewave absorption example with a lossy dielectric, showcasing how verlan initializations can mitigate getting trapped by sub-optimal local minima in photonic inverse design problems.

Abstract The generalized maximum-entropy sampling problem (GMESP) is to select an order- s principal submatrix from an order- n covariance matrix, to maximize the product of its t greatest eigenvalues, $$0<t\le s <n$$ 0 < t ≤ s < n . Introduced more than 25 years ago, GMESP is a natural generalization of two fundamental problems in statistical design theory: (i) maximum-entropy sampling problem (MESP); (ii) binary D-optimality (D-Opt). In the general case, it can be motivated by a selection problem in the context of principal component analysis (PCA). We introduce the first convex-optimization based relaxation for GMESP, study its behavior, compare it to an earlier spectral bound, and demonstrate its use in a branch-and-bound scheme. We find that such an approach is practical when $$s-t$$ s - t is very small.

Proton minibeam radiation therapy (pMBRT) employs spatially fractionated dose distributions to reduce normal tissue toxicity. A key component is the multi-slit collimator (MSC), which shapes the beam into narrow, spatially separated minibeams. Small lateral shifts of the MSC relative to the beam direction can substantially alter peak-valley dose patterns, target coverage, and organs-at-risk (OAR) sparing, making MSC positioning a critical planning parameter. We develop a novel collimator position optimization (CPO) algorithm for pMBRT that allows independent lateral shifts of the MSC at each beam angle to improve plan quality. The problem is formulated as a mixed-integer programming (MIP) model that jointly optimizes MSC positions and spot intensities. Binary variables select candidate lateral shifts per beam angle, while continuous variables represent spot intensities. The resulting non-convex problem is solved using an augmented Lagrangian framework with iterative convex relaxation and alternating direction method of multipliers (ADMM) decomposition. In three clinical cases, the proposed method achieved near-optimal solutions with substantially reduced computation time compared to exhaustive enumeration (e.g., 700 s vs. 15,000 s for an abdominal case). Allowing multiple MSC positions per beam angle led to consistent dosimetric improvements, particularly in OAR sparing; for example, mean oral cavity dose in a head-and-neck case decreased from 6.5 Gy to 4.6 Gy. MSC position optimization enhances pMBRT plan quality and can be efficiently integrated into clinical treatment planning.

Global Optimization for Generalized Linear Multiplicative Programs Through Simplicial Branch and Bound with Convex Relaxation

Quantum algorithms for estimating the ground state energy of a quantum system often operate by preparing a classically accessible quantum state and then applying quantum phase estimation. Whether this approach yields quantum advantage hinges on the state's energy spectrum, that is, the sequence of the state's overlaps with the energy eigenstates of the system Hamiltonian. We show that the energy spectrum of any entanglement-compressed quantum state must have large support if most energy eigenstates are highly entangled, an assumption supported by the eigenstate thermalization hypothesis. Furthermore, we show that if the compressed quantum state minimizes expected energy, then its energy spectrum decays with the inverse-squared energy eigenvalues under a convex relaxation of the compression constraint. This explains the main empirical finding of Silvester et al. [Phys. Rev. Lett. 134, 126503 (2025).PRLTAO0031-900710.1103/PhysRevLett.134.126503] that the energy spectra of matrix product states do not decay exponentially.

A key challenge limiting the industrial deployment of large-scale heavy-haul trainsets, such as 30,000-ton units, is ensuring stable maneuverability under complex railway and communication conditions. This study investigates the differences in communication range and maneuverability between traditional combined trains and emerging virtual coupled train platoons, focusing on the existing rail communication-signaling infrastructure. The goal is to ensure reliable communication and stable maneuverability in complex environments. To achieve this, a general optimisation model is developed for maneuvering heterogeneous heavy-haul trainsets with various configurations. The model incorporates both longitudinal train dynamics and communication factors, achieving convex relaxation by equivalently modelling three-dimensional railway lines, ensuring a globally efficient and stable solution. The model and method are validated through engineering case studies, demonstrating significant improvements in maneuverability safety and energy efficiency. The analysis shows that the primary safety concerns for 30,000-ton combined trains and train platoons are longitudinal impact and rear-end collisions, respectively. Based on these findings, it is recommended that the train to train(T2T) communication delay limits be set at 2.70 s for combined trains and 5.70 s for train platoons.

The growing adoption of hydrogen-powered transport demands scalable and robust control strategies for hydrogen refueling stations (HRSs) operating under uncertainty in supply, demand, and market conditions. This study presents a delay-aware predictive control framework for renewable-integrated HRSs equipped with heterogeneous hydrogen storage systems and dual interaction with electricity and hydrogen markets. The station architecture enables simultaneous refueling across multiple fuel cell electric vehicle (FCEV) classes, each served by a dedicated high-pressure storage unit, while auxiliary tanks function as buffers and market reserves. Inter-storage coordination is managed using receding horizon control across layered decision stages, allowing flexible hydrogen routing under dynamic operating conditions. A central challenge is hydrogen delivery delays, which introduce a temporal gap between procurement actions and actual availability. The proposed formulation incorporates these delays within the control horizon, classified as deterministic (fixed lead times), stochastic (modeled via discrete uncertainty sets), and logistics-based (dependent on route planning, fleet capacity, and congestion, captured through time-dependent concave functions). A mode-switching mechanism allows the operator to activate one of four control strategies: deterministic MPC (DMPC), scenario-based stochastic MPC (SMPC), convex relaxed MPC (RMPC), and scaled risk-averse SMPC (SRA-SMPC) with conditional value-at-risk and chance constraints. Convex relaxation techniques are applied to address combinatorial complexity from binary variables, nonlinear tank dynamics, and inter-market constraints, ensuring real-time tractability while preserving constraint feasibility and economic performance. Numerical simulations confirm the framework’s effectiveness in coordinating storage operations, meeting demand, and reducing costs under uncertainty, with substantial computational benefits compared to conventional approaches. • A unified control framework manages delivery delays in hydrogen refueling stations. • Delay types are categorized to enable adaptive controller selection. • Predictive strategies balance dispatch accuracy and computational efficiency. • Market interaction and renewable coordination improve cost and sustainability. • The framework supports future expansion to multi-station hydrogen networks.

We derive a convex relaxation principle for a large class of non convex variational problems where the functional to be minimized involves a one homogeneous gradient energy. This applies directly to free boundary or multiphase problems in the case of the classical total variation or of some anisotropic variants. The underlying argument is an exclusion principle which states that any global minimizer avoids taking values in the intervals where the lower order potential is nonconvex. This allows using duality methods and deriving a saddle point characterization of the global minimizers. A numerical validation of our principle is presented in the case of several free boundary and multiphase problems that we treat through a primal-dual algorithm. The accuracy of the interfaces and the convergence of the algoritm benefit in a large way of a new epigraphical projection method that we introduced to tackle the non differentiability of the convexified Lagrangian.

In automotive radar systems, super-resolution angle-of- arrival estimation is required to distinguish neighboring targets for good imaging results. In order to adapt to the rapidly changing dynamic driving scenarios, the angle of arrival estimation must meet the real-time computational requirements and can be realized based on single snapshot data. Existing gridless methods based on atomic norm minimization are suitable for super-resolved angle estimation for single snapshot signals and avoid the grid mismatch problem. But the atomic norm acts as a convex relaxation of the atomic <italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">l₀</i> norm, and it leads to a degradation of the angle resolution. Moreover, it has high computational complexity and cannot meet the real-time requirement. In order to avoid the convex relaxation and satisfy the real-time solution requirements at the same time, we propose a single snapshot multi-objective DOA estimation model. It minimizes both the measurement error term and the atomic <italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">l₀</i> norm, and directly deals with the conflict between angle measurement accuracy and signal sparsity to jointly achieve target detection and angle estimation. To efficiently solve this NP-hard model, we design a cooperative coevolutionary crayfish optimization algorithm. The novelty is two-fold: one is the proposal of a new cooperative co-evolutionary decomposition strategy for multiple populations to efficiently decompose the multi-objective DOA estimation model into multiple single-objective problems. The crayfish optimization algorithm is applied to each population to solve the multi-objective problem with variable dimensions based on multiple swarm cooperative co-evolutionary techniques. The other innovation is the proposal of a variable-length neighboring orthogonal crossover operator for information exchange between populations. This mechanism effectively speeds up the convergence of the algorithm. Simulation results and hardware experiments show that, compared with the spatially smooth MUSIC algorithm, the atomic-norm-based gridless algorithm and its various variants, and the variational Bayesian-based gridless algorithm, our method exhibits superior performance in terms of angle measurements and target number detection. Moreover, it can realize angle measurement within 20ms, which is the only gridless algorithm among these gridless super-resolution algorithms that can meet the real-time computation requirements of automotive radar.

The pervasive threat of jamming attacks, particularly from adaptive jammers capable of optimizing their strategies, poses a significant challenge to the security and reliability of wireless communications. This paper addresses this issue by investigating anti-jamming communications empowered by an active reconfigurable intelligent surface. The strategic interaction between the legitimate system and the adaptive jammer is modeled as a Stackelberg game, where the legitimate user, acting as the leader, proactively designs its strategy while anticipating the jammer’s optimal response. We prove the existence of the Stackelberg equilibrium and derive it using a backward induction method. Particularly, the jammer’s optimal strategy is embedded into the leader’s problem, resulting in a bi-level optimization that jointly considers legitimate transmit power, transmit/receive beamformers, and active reflection. We tackle this complex, non-convex problem by using a block coordinate descent framework, wherein subproblems are iteratively solved via convex relaxation and successive convex approximation techniques. Simulation results demonstrate the significant superiority of the proposed active RIS-assisted scheme in enhancing legitimate transmissions and degrading jamming effects compared to baseline schemes across various scenarios. These findings highlight the effectiveness of combining active RIS technology with a strategic game-theoretic framework for anti-jamming communications.

This paper investigates multi-antenna beamforming for integrated sensing and covert communication (ISCC) in consumer-electronics-assisted relaying networks. Specifically, we consider an amplify-and-forward (AF) relay network where a dual-functional multi-antenna transmitter Alice and one multi-antenna AF relay jointly serve a consumer electronics receiver Bob, in the presence of one warden Willie, while simultaneously conducting environmental sensing. Recognizing the hardware constraints and dense deployment characteristics of consumer electronics, we exploit the high spatial degrees of freedom offered by multi-antenna arrays. The proposed system steers beams to deal with the interference at the warden while maintaining high-gain links for legitimate users, thereby turning the dense multipath environment into an advantage for covertness. The objective is to maximize Bob’s signal-to-interference-plus-noise ratio (SINR) under stringent covertness and power constraints while maintaining a minimum sensing quality. To this end, we jointly optimize Alice’s communication/sensing beamforming and the relay amplification factor through a fractional programming framework. The covertness constraints are reformulated into tractable power-domain inequalities by employing total variation distance, KL divergence, and S-procedure-based convex relaxations, where three eavesdropper channel state information (ECSI) of instantaneous, bounded-error, and statistical scenarios are considered. These scenarios are carefully selected to cover the full spectrum of warden capabilities, from active detection to passive surveillance. The resulting non-convex problem is solved via semidefinite relaxation combined with Dinkelbach’s method and alternating optimization, complemented by relay antenna selection to enhance consumer-device link reliability. Simulation results are provided to validate the effectiveness of the proposed method. In particular, for a 30 dBm transmit power, the proposed scheme with instantaneous ECSI achieves approximately 14 dB SINR at the consumer receiver, compared to 5.5 dB and 3.5 dB under bounded-error and statistical ECSI, respectively, demonstrating the proposed scheme’s robustness and effectiveness in covert, consumer-electronics-integrated wireless networks.

In this paper, a multi-UAV enabled integrated sensing, computing, and communication (ISCAC) system model is proposed, in which multiple UAVs sense ground Users and offload sensing data to a high altitude platform (HAP) for processing through mobile edge computing (MEC). To maximize sensing data acquisition while minimizing total energy consumption of the ISCAC system, we present a trade-off optimization problem between UAV radar sensing and total energy consumption of the system. We transform the established non-convex optimization problem into three subproblems: sensing scheduling optimization, UAV transmit power optimization, and UAV-HAP trajectory optimization. We solve these subproblems using successive convex approximation (SCA) and relaxation methods, and propose a three-layer iterative optimization algorithm to solve the original optimization problem. Simulation results demonstrate that, compared to the benchmark schemes, the proposed algorithm can significantly improve the system performance.

The large-scale deployment of hydrogen-powered transportation requires economically viable and operationally resilient control strategies for hydrogen refueling stations (HRSs) operating under renewable variability, dynamic hydrogen demand, and logistics-induced delivery delays associated with fueling hydrogen mobility. To address this, this study develops a relaxation-based model predictive control (RMPC) framework for a renewable-integrated HRS equipped with heterogeneous hydrogen storage units, battery energy storage systems (BESS), and dual participation in hydrogen and electricity markets. The system architecture incorporates high-pressure tanks dedicated to distinct fuel cell electric vehicle (FCEV) classes, along with auxiliary buffer and backup tanks that enable inter-storage routing and enhance flexibility in meeting demand. Hydrogen procurement from external markets is scheduled through the backup tank, providing temporal decoupling between uncertain deliveries and refueling requirements. Based on time-dependent concave delay functions, a logistics-aware delivery delay model is introduced to capture time-varying lead times arising from routing constraints, fleet availability, and congestion. This battery-hydrogen co-optimization strategy enhances operational flexibility and resilience, effectively decoupling on-site hydrogen supply from short-term market volatility and demand fluctuations. A convex relaxation approach addresses the complexity of binary tank-switching decisions, nonlinear dynamics, and inter-market coordination, while maintaining real-time scheduling capability. The proposed RMPC is assessed under idealized and delay-affected market conditions to establish bounds on service continuity and operating cost. In particular, coupling BESSs with hydrogen storage improves station resilience and economic feasibility. The proposed approach yields a significant economic advantage over conventional controllers, with simulations indicating up to a 57% increase in profit when using multi-day predictive optimization, while the HRS achieves > 98% refueling availability and reduces energy costs by approximately 15% relative to baseline operation. These results demonstrate the potential of the proposed framework to enable robust, cost-effective, and scalable operation of next-generation hydrogen refueling infrastructure.