Papers

Bilevel optimization provides a natural modeling language for hierarchical decision problems. However, applying existing numerical solvers usually requires substantial manual analysis and reformulation. In this paper, we introduce disciplined bilevel programming (DBLP), a symbolic framework that allows users to specify and solve optimistic bilevel problems in a high-level, human-readable way that is close to the mathematical formulation. For problems with a disciplined nonlinear upper problem and a convex lower problem satisfying the disciplined parameterized programming rules, DBLP automatically canonicalizes the lower problem into conic form and constructs an equivalent single-level reformulation using the conic Karush-Kuhn-Tucker conditions. We relax the resulting complementarity constraint and use a gap continuation procedure to approximately solve a sequence of smooth nonlinear problems. We implement DBLP in the open-source Python package BLVPY, an extension of CVXPY for bilevel programming. We demonstrate the modeling and solution capabilities of BLVPY on a range of bilevel optimization problems from several application domains. The proposed framework and implementation allow users to specify and solve bilevel optimization problems within a few lines of code, without prior expertise in bilevel modeling and numerical optimization.

Powell asked whether the smoothness assumptions underlying classical superlinear convergence force a fixed power law between adjacent iterates of exact-line-search variable-metric methods. We answer this question negatively for BFGS: within the smooth strongly convex setting, the smallest possible adjacent-iterate Q-order is one, and this boundary is attained by a single nonterminating run. In every finite dimension at least two, and for any prescribed radius and Hessian tolerance, we construct an infinitely differentiable, globally strongly convex objective that equals the standard quadratic outside the corresponding ball and whose Hessian remains within the prescribed tolerance of the identity in operator norm. The objective has its unique minimizer at the origin and identity Hessian there. Exact-line-search BFGS, initialized with the identity matrix and started inside that ball, converges Q-superlinearly, yet no fixed power greater than one controls all sufficiently late adjacent errors.

We study constrained smooth optimization problems on Hadamard manifolds with closed geodesically convex feasible sets. We analyze two projected gradient schemes: one with a constant stepsize and another with a backtracking line search. The constant-stepsize scheme is analyzed under the assumption that the objective function has a Lipschitz continuous Riemannian gradient, whereas the backtracking variant does not require this assumption to establish stationarity of accumulation points. For both schemes, we prove that every accumulation point of the generated sequence is first-order stationary under the respective assumptions, without requiring compactness of the feasible set; compactness is needed only to ensure the existence of accumulation points. When the objective function has a Lipschitz continuous Riemannian gradient, we derive iteration-complexity bounds of order \(O(1/\sqrt{N})\) for projection-based stationarity measures for both schemes, together with the corresponding \(\varepsilon\)-complexity estimates. For the backtracking scheme, the complexity analysis additionally requires the trial line-search stepsizes to be uniformly bounded away from zero. Under the same respective assumptions, the generated sequences are also asymptotically regular. Finally, we illustrate the practical performance of the methods by solving constrained Karcher mean problems on the manifold of symmetric positive definite matrices.

We establish a whole-space solution framework in augmented spectral Barron spaces for uniformly elliptic parabolic equations with state-dependent principal coefficients. A frozen-symbol parametrix yields bounded Green and terminal operators and a one-derivative smoothing estimate on arbitrary finite horizons, without requiring the spatial variation of the principal coefficient to be perturbatively small. We then apply this linear framework to a finite-horizon Hamilton–Jacobi–Bellman equation. A semi-explicit gradient iteration converges on sufficiently short horizons at a Gamma-factorial rate in a whole-space augmented Barron norm, and its limit is identified with the stochastic-control value function. Finally, temporal Jackson approximation and spatial spectral Barron approximation give joint shallow cosine-network approximations in space and time for the value function and the optimal feedback, with quantitative error and neuron-count bounds. Taken together, these results connect variable-coefficient parabolic well-posedness and smoothing in augmented spectral Barron spaces with nonlinear HJB solvability and quantitative neural-network approximation.

In this work, we study radically uncoupled learning in discounted general-sum Markov games. Assuming “$\mathsf{ETH}$ for $\mathsf{PPAD}$", we show that, for every fixed discount factor, there is no polynomial-time algorithm for computing inverse-polynomially accurate coarse correlated equilibria in discounted general-sum Markov games when players learn independently in decentralized settings. Complementing this hardness result, we provide what appears to be the first radically uncoupled algorithm with sub-exponential convergence guarantees to coarse correlated equilibria in discounted general-sum Markov games without imposing any structural restrictions on the game. Our algorithm is a layered variant of optimistic mirror descent with an increasing step-size schedule tailored to the multi-agent setting. Finally, we develop both full-feedback and partial feedback versions of the aforementioned algorithm and establish sub-exponential convergence guarantees for each case.

The history of extremum seeking is not merely the history of an algorithm; it is the history of an idea. Few ideas in control engineering have demonstrated the remarkable longevity of extremum seeking control (ESC). Invented more than one hundred years ago, ESC has continually reinvented itself while remaining faithful to its original objective: enabling systems to optimize their performance without relying on accurate mathematical models. This perspective article proposes that the evolution of ESC is best understood through three scientific revolutions. The first Engineering Revolution (1922-1999) established the engineering principles of model-free optimization; the second Mathematical Revolution (2000-2010) provided the rigorous mathematical foundations that transformed ESC into a mature discipline of nonlinear control; and the third ongoing Infinite-Dimensional and Cyber-Physical Revolution (2010-present) continues to expand its scope toward delays, partial differential equations, distributed optimization, event-triggered implementations, and increasingly complex cyber-physical systems. Beyond recounting this historical evolution, we offer a personal perspective on why ESC has remained relevant across successive technological eras. We argue that its enduring influence arises because the fundamental engineering challenge has never changed: "how can a dynamical system learn to improve its own performance when the optimum is unknown?" As optimization, learning, and feedback control become increasingly intertwined, ESC appears uniquely positioned to contribute to a new generation of intelligent autonomous systems, pointing toward what may become the field's fourth scientific revolution.

This paper introduces a structural taxonomy for constrained non-convex optimization based on the signature of Lagrange multipliers at KKT stationary points. Leveraging a unified game-theoretic interpretation of eight classical algorithm families–including block coordinate descent, ADMM, generalized Benders decomposition, successive convex approximation, interior-point methods, mirror descent, Frank-Wolfe, and Riemannian gradient descent–we show that the normalized multiplier vector carries an algorithm-independent structural fingerprint. Four scale-free shape features of this vector partition the dual space into five operational regimes: Unconstrained, Resource-Limited, Saturation, Strongly-Coupled, and Hybrid. We establish four structural theorems characterizing the partition: invariance under natural KKT symmetries, local stability under data perturbation with explicit Lipschitz margins from Robinson's strong regularity, codimension-one regime transitions, and the topological identification of the Hybrid regime as the Lebesgue-null boundary of the core regimes. A linear-time classifier is proposed with provable guarantees on correctness, iteration stabilization, sample complexity, and online tracking under data drift. Numerical experiments on 104 mixed-integer nonlinear programs and a downlink beamforming instance validate the theoretical predictions. The framework provides a foundational tool for regime-aware algorithm design and robustness analysis in non-convex optimization.

In power markets where the physics of the transmission grid is closely represented in market clearing models, certain cross-zonal power exchanges can take place only if other exchanges occur concurrently. This leads to challenges in cross-zonal balancing capacity markets, where the activation of reserves in real time remains uncertain. Ensuring reserve deliverability in all activation scenarios is naturally modeled as a stochastic programming (SP) problem, whose size grows exponentially with the number of locations in the network. This formulation scales poorly in real-world applications. We first show that activation scenarios can be expressed in an order-book-agnostic way, reducing the challenge to a network modeling problem and improving computational performance. We then introduce a general inner approximation principle that we use to derive two scalable inner approximations and one column generation algorithm for tackling one of the approximations. The first inner approximation is well known to practitioners, and relates to a basic result for describing boxes in a polytope, while the second model, its associated column generation algorithm and finite-dimensional reformulation are based on semi-infinite linear programming and robust linear optimization. We compare the inner approximations to the exact SP formulation, which we also solve via a Dantzig-Wolfe decomposition for comparison purposes. Numerical results show that these inner approximations are much more scalable than the stochastic programming formulation, while reaping most of the benefits of cross-zonal exchanges. The approaches are of particular interest for future pan-European cross-zonal balancing capacity markets, and can also accommodate co-optimization of energy and balancing capacity products.

As transit agencies consider bus fleet electrification, understanding the efficiency and cost of replacing diesel buses (DBs) with battery electric buses (BEBs) is critical. To evaluate this, this study applies a mixed-fleet optimization model, integrating scheduling, charging, and fleet composition decisions, across five agencies: Santa Monica's Big Blue Bus (BBB), the Chicago Transit Authority (CTA), Knoxville Area Transit (KAT), the Metropolitan Atlanta Rapid Transit Authority (MARTA), and Manhattan's Metropolitan Transportation Authority (MTA) bus service. By calculating electric fleet share, the BEB/DB replacement ratio, transit-link density, and vehicle activity-time allocation, the study finds that while optimized fleets remain majority-electric, vehicle substitution is rarely one-to-one. Average replacement ratios range from 1.101 for CTA to 1.245 for KAT, with higher transit-link density networks like CTA and MTA requiring fewer replacement buses per diesel bus displaced than lower-density networks like MARTA and KAT. While these relationships are descriptive rather than causal, non-revenue vehicle activity may help explain the differences. By shifting the focus from simple electric fleet share to diesel replacement efficiency, this multi-agency comparison demonstrates that transit agencies should use the replacement ratio to accurately forecast additional fleet capacity requirements and avoid the costly assumption of strict one-to-one vehicle substitution.

We study the Blaschke–Santaló diagram associated with the volume, the diameter, and the Cheeger constant. In the class of bounded open subsets of $\mathbb{R}^m$, $m\ge 2$, we give a complete description of the diagram. In the class of convex bodies of $\mathbb R^m$, we prove that the diagram is closed and simply connected as it is given by the region between two continuous functions for which qualitative properties such as monotonicity, local behavior, and growth estimates are studied.

Augmented Lagrangian methods are among the most effective approaches for solving constrained convex optimization problems. However, classical complexity analyses of first-order methods applied within the augmented Lagrangian framework usually rely on the assumption that the objective function has a Lipschitz continuous gradient. This assumption excludes an important class of generalized smooth functions whose gradients may grow unboundedly. In this paper, we study an inexact augmented Lagrangian method for solving linearly constrained convex optimization problems with $(L_0,L_1)$-smooth objective functions. We show that the augmented Lagrangian subproblems preserve the $(L_0,L_1)$-smooth structure, with parameters depending on the penalty coefficient. This property allows us to employ recent accelerated first-order schemes designed for generalized smooth optimization instead of classical smooth optimization methods. In particular, we combine the inexact augmented Lagrangian framework with a two-stage acceleration procedure based on clipped gradient descent and accelerated optimization.

This article investigates the problem of assigning multiple identical interceptors over multiple identical aerial threats, where all interceptors are launched in a single salvo. For this problem, two strategies have been studied in the literature: (A) to spread all interceptors as evenly as possible over all threats, and (B) to randomly assign all interceptors over the threats. The main contributions of this article are as follows. First, the literature contains empirical evidence that Strategy A is more efficient than Strategy B in terms of the mean number of missed threats, when the number of interceptors is no less than the number of threats. This article gives a rigorous proof of this fact. Secondly, it demonstrates numerically that Strategy~A can be significantly more efficient than Strategy B. Thirdly, it shows that Strategy~A is not only superior to Strategy B, but also the bona fide optimal strategy. Finally, this article establishes that these conclusions also hold when the number of interceptors is less than the number of threats.

In this work we are concerned with the controllability of time-varying lumped control systems governed by a semilinear differential equation. We consider the semilinear system as a perturbation of a linear system. Assuming the underlying linear system is controllable, and the nonlinear forcing function satisfies a boundedness condition which is adapted to the underlying linear system, we show that the semilinear system is also approximately controllable

We develop a data-driven framework for nonlinear port-Hamiltonian (pH) systems based on interpolation conditions to characterize consistency between observed data and structured dynamical models. Specifically, we derive necessary and sufficient conditions for the existence of a pH system with a smooth (convex) Hamiltonian instantly consistent with a given dataset, without requiring explicit parametrization. We further provide a semidefinite programming formulation to verify consistency with non-degenerate interconnection and dissipation structures. Our results provide a principled approach to assess instant data consistency with physical structure and pave the way for control design directly from data.

Traditional studies on evolutionary dynamics of cooperation have concentrated on an idealized game setup free of payoff-observation errors. However, in real-world scenarios, individuals frequently encounter errors when observing the payoffs of their opponents during game interactions, resulting from unintentional mistakes, such as data misrecording, overlooking critical details, and miscalculations. This gap between idealized error-free game models and the error-prone real-world game interactions leads to the lack of insight into the impact of payoff-observation errors on the evolutionary dynamics of population games, in particular the efficiency of cooperation incentives. In this paper, we construct a research framework for population games with payoff-observation errors, which enables us to investigate the effects of errors on the evolutionary dynamics of cooperation in the evolutionary Prisoner's Dilemma game with combined incentives. To quantify the implementation costs of incentives in the presence of errors, we devise an index function and employ optimal control theory to derive the optimal incentive protocols. Our theoretical and numerical results reveal that payoff-observation errors can lower the costs compared to error-free cases, and we also derive the theoretical conditions for these results. Finally, we formulate an optimization problem to explore the cost difference between the optimal incentive protocols with and without errors, and further design an algorithm to obtain the numerical solution that minimizes this difference.

Community detection methods must balance two competing objectives: identifying small, cohesive groups while avoiding excessive fragmentation. Modularity, the most widely adopted optimization criterion, typically merges small communities in large networks due to its resolution limit. In contrast, a persistence-based criterion promotes more granular partitions. We introduce Scaled Null-Adjusted Persistence (Scaled-NAP), a parametric family of quality functions that incorporates both these criteria. The definition exploits the exact identity between a cluster's modularity contribution and its Null-Adjusted Persistence (NAP) multiplied by its relative volume. Raising this volume factor to a parameter $\alpha\in[0,1]$ yields NAP at $\alpha=0$ and modularity at $\alpha=1$, while intermediate values control the scale of the detected partition. We derive conditions under which merging two communities improves the objective function and characterize the emergence of scale dependence, including resolution-limit behaviour on Caveman graphs. We develop the Milano algorithm, a multilevel Louvain-style heuristic for optimizing Scaled-NAP on large networks. Experiments on weighted and unweighted Lancichinetti-Fortunato-Radicchi benchmarks show that Scaled-NAP achieves the highest or tied-highest recovery wherever the ground truth structure is detectable, with its advantage increasing under community-size heterogeneity. Tests on three real networks with up to 1.1 million nodes confirm its capability to identify fine-grained ground-truth communities. The Milano algorithm also turned out to be the fastest method evaluated on large networks. These results show that Scaled-NAP provides an effective and scalable bridge between modularity-based and persistence-based community detection methodologies.

The dual active-set method of Goldfarb and Idnani solves the strictly convex quadratic program by adding and dropping one constraint at a time, requiring no phase one. Primal–dual active set and block principal pivoting skip the walk altogether: they guess an entire active set at once and repair it from returning signs. For bound constraints the guess is safe. The system solved on a candidate set is a principal submatrix of $G^{-1}$, positive definite regardless of the guess; this $P$-matrix property guarantees finite termination. For $C^\top x \ge b$ that hypothesis fails, which is our central focus. The working-set system is $C_{\mathcal{A}}^\top G^{-1} C_{\mathcal{A}}$, positive definite only when $C_{\mathcal{A}}$ has full column rank, a property of the guess, not the data. The $P$-matrix property is lost and the structural obstruction invalidates the guarantee. Strict convexity keeps the method usable: KKT conditions make candidate sets certified rather than trusted. Over 600 random instances across four constraint families this failure never occurs. However, duplicating columns of $C$ raises failure rates to 93%, dropping certified fraction from 100% to zero. On real long-only portfolio data, the active set reaches 442 of $494 + 1$ constraints.

This comment revisits Lemma 1 in [1], which plays a central role in the event-triggered stabilization analysis developed therein. We identify technical gaps in the proof of the lemma and provide a corrected argument. In particular, careful treatment of the exponentially decaying term shows that its decay rate must be retained in the resulting convergence estimate. The statement of the original lemma, with the exponential decay rate determined by the minimum of the characteristic decay rate and the decay rate of this term, remains valid.

The symmetry measure of a convex body $S\subset\mathbb{R}^n$ is given by: $\mathrm{sym}(S):=\max\{\alpha\ge0:\text{ there exists }x\in S\text{ such that }-\alpha(S-x)\subseteq S-x\}$, where such an $x$ is called a Minkowski center. We prove that every convex body $S$ admits a $\sqrt{\frac{n}{\mathrm{sym}(S)}}$-rounding of $S$, namely, there exists an origin-centered ellipsoid $E$ and a center $c$ such that $E\subseteq S-c\subseteq\sqrt{\frac{n}{\mathrm{sym}(S)}}\,E$. This result was conjectured in 2005 by Belloni and Freund. As special cases, this recovers an $n$-rounding of $S$ (since $\mathrm{sym}(S)\ge\frac{1}{n}$), and a $\sqrt{n}$-rounding when $\mathrm{sym}(S)=1$. In the case when $S$ is a polytope given as the convex hull of points, the desired rounding is produced by a regularized minimum-volume covering ellipsoid problem where the regularization is with respect to the Minkowski center. Similarly, when $S$ is a polytope given as the intersection of halfspaces, such a rounding is produced by a regularized maximum-volume inscribed ellipsoid problem. In both of these cases, the rounding can be computed by first solving a linear optimization problem (to compute $\mathrm{sym}(S)$ and a Minkowski center), and then solving a convex optimization problem with a logarithmic determinant objective, second-order cone constraints, and one semidefinite cone constraint. We also show that the factor $\sqrt{\frac{n}{\mathrm{sym}(S)}}$ is nearly tight in its dependence on dimension and symmetry. When $\frac{n+1}{1+\mathrm{sym}(S)}$ is an integer, we show by explicit construction that the factor $\sqrt{\frac{n}{\mathrm{sym}(S)}}$ is tight. In the more general case, for every dimension $n$ and every admissible symmetry value, we construct a polytope $S$ for which every rounding factor is at least $\sqrt{\frac{2}{3}}\sqrt{\frac{n}{\mathrm{sym}(S)}}$.

This paper considers decentralized contingency MPC for multi-agent control under a state-only information pattern, with particular focus on limited sensing and plug-and-play operation. The objective is to retain recursive feasibility, safety, and Lyapunov-type convergence while reducing conservatism in local interaction handling. The framework relies on agent-wise fallback regions (safe sets) in which a feasible contingency maneuver to a safe equilibrium is always available. A novel safe-set update mechanism is introduced that supports less conservative decentralized interaction while preserving the underlying guarantees. This, in turn, enables memory-free local interaction and finite sensing ranges without requiring agents to reconstruct the exact neighbor geometry. The resulting scheme remains fully decentralized and preserves the shared-first-input contingency MPC structure. Theoretical guarantees and simulation results illustrate the effectiveness of the approach in dense multi-agent scenarios.

The real stability radius problem asks for the smallest real perturbation that makes a given matrix lose rank, and it is a structured counterpart of the corresponding complex stability radius problem. In this paper, we consider an angular version of this question for sectorial matrices, where the size of a perturbation is measured by its largest phase instead of its largest singular value. More precisely, for a sectorial matrix $M$, we investigate the minimum required largest phase of a real and sectorial perturbation $\Delta$ for which $I+M\Delta$ loses rank. We show that for some matrices, the complex and the real angular stability radius are the same, but that they are in general different, such as when $M$ is diagonal or complex symmetric. For these cases, we derive a formula for the real angular stability radius, and this formula depends on the two largest phases of $M$. Moreover, in these cases we show that the result for $n \times n$ matrices can be reduced to that for $2 \times 2$ matrices by compressing an arbitrary real destabilizing perturbation to the real two-dimensional subspace spanned by the real and imaginary parts of a destabilizing vector.

We study equilibrium computation in simulation-based games motivated by competition in natural-hazard insurance markets. The motivating application considers a multi-region market in which each insurer selects regional pricing and reinsurance decisions to maximize expected financial performance subject to multi-period insolvency constraints. Competitors' decisions affect both profitability and insolvency feasibility through the insured portfolio. Consequently, the strategic interaction defines a generalized Nash game with discontinuous player objectives induced by discrete household purchase decisions and potentially nonconvex feasible sets generated by the insolvency constraints. Discontinuity, nonconvex feasibility, and simulation-defined objectives and constraints complicate the direct use of standard optimality-based equilibrium formulations. Moreover, repeated evaluation of the multi-period simulator makes equilibrium search computationally expensive. To address these computational challenges, we develop an evaluation-efficient framework for computing approximate local equilibria, which combines a response optimizer tailored to the distinct structures of pricing and reinsurance with damped better-response dynamics to mitigate cycling. In a case study of the North Carolina hurricane insurance market, the framework identifies multiple approximate local equilibrium candidates both with and without a household affordability cap and finds no profitable sampled unilateral deviations within the tested neighborhoods under the original simulation model. A representative equilibrium candidate is computed in approximately 5–10 minutes on a single GPU, making repeated equilibrium analyses computationally practical.

The Poisson–Nernst–Planck (PNP) equations strongly couple ionic transport and electrostatic interactions through the Poisson equation, posing substantial numerical challenges under small permittivity and complex potential boundary conditions. Underlying these equations is a natural Wasserstein gradient-flow structure, in which the Poisson equation serves as a local realization of the nonlocal electrostatic interaction energy. Exploiting this structure, we formulate each time step as a constrained convex minimization problem where the ionic continuity equations and the Poisson equation are incorporated as linear constraints, allowing the concentrations, fluxes, and electrostatic potential to be updated simultaneously. The variational structure of the scheme intrinsically guarantees the dissipation of the original free energy, mass conservation, and nonnegativity of ionic concentrations under general electrostatic boundary conditions. Moreover, the framework is structurally modular: extending from classical to modified PNP models with steric interactions and concentration-gradient corrections requires only modifying the energy functional, while all structure-preserving properties are automatically retained. To efficiently solve the resulting large-scale constrained problems, we develop preconditioned and transformed primal–dual algorithms equipped with tailored fast dual solvers, namely DCT-based direct and Schur-complement iterative methods, that exploit the coupled block structure of the PDE constraints. Numerical experiments on classical and modified PNP systems demonstrate the accuracy and structure-preserving properties of the scheme, and show that the proposed algorithms converge reliably in strongly coupled small-permittivity regimes without significant growth in computational cost.

The problem of global minimization of nonconvex quadratic functions subject to box constraints is studied. Applying Gershgorin theorem we reduce the main matrix A to its diagonal form which is used to obtain difference of convex representation of the quadratic function. Then {\epsilon}-subdifferentials of DC components are studied and used to design a method for minimizing the quadratic function globally.

Given two quadratic functions \( f(x) = x^T Ax + 2a^T x + a_0 \) and \( g(x) = x^T Bx + 2b^T x + b_0 ,\) each associated with either the strict inequality ($<0$); non-strict inequality ($\leq 0$); or the equality ($=0$), it is a fundamental question to ask whether or not the joint system has a solution. For homogeneous quadratic systems ($a=b=0,~a_0=b_0=0$), starting from Finsler's lemma in 1936 until Yuan's alternative lemma in 1990, all combinations of the unsolvability for $\{x\in \mathbb{R} ^n\mid x^T Ax \mathbin{\star } 0\}\cap \{x\in \mathbb{R} ^n\mid x^T Bx \mathbin{\#} 0\}\subset \{0\}$, where $\star $ and $\#$ can be any of $\{<,\leq ,=\}$, have been shown to possess either a positive definite or a positive semi-definite matrix pencil of $A$ and $B.$ Extensions to nonhomogeneous quadratic systems $\{x\in \mathbb{R} ^n\mid f(x) \mathbin{\star } 0\}\cap \{x\in \mathbb{R} ^n\mid g(x) \mathbin{\#} 0\}=\emptyset $ have been done for several cases already. Two challenging cases remain open: the nonhomogeneous Calabi Theorem which determines when $\{f(x)=0\}\cap \{g(x)=0\}=\emptyset $; and the nonhomogeneous (strict) Finsler lemma to determine whether $\{f(x)\leq 0\}\cap \{g(x)=0\}=\emptyset .$ The paper provides the answers to both, in theorems and algorithms.

Automated material handling systems in semiconductor fabs are operated by a material control system (MCS) that must schedule a relay route for every transport command online, before execution. This is a data-driven scheduling problem in which route cost is dominated in the upper tail by queueing at heterogeneous, partially observable relay equipment, so route selection requires estimating both delivery time and congestion risk at the decision moment. This paper proposes a transport-network-aware dynamic congestion representation (TN-DCR). Built on a static directed transport graph induced by historically observed relay segments, TN-DCR combines structural route priors, multi-window network-wide congestion context, route-level bottleneck exposure, and an inductive graph-aware route embedding, all constructed under a prediction-time-safety invariant that admits only information observed strictly before the prediction moment. The representation feeds separate queue- and transfer-time regressors and an ordinal multi-label classifier producing calibrated multi-threshold exceedance scores, with an empirical-Bayes stock-key residual correction reducing systematic queue-time underprediction. The predictions serve as costs in a risk-constrained route-scheduling rule that minimizes predicted delivery time subject to a bound on extreme-congestion probability, embedding the learned predictors within a lightweight operations-research decision model. In a controlled closed-loop evaluation, mean delivery time falls by 16.4% and internal resource waiting time by 22.6% while throughput remains essentially unchanged.

Large-scale industrial robot fleets share constrained physical infrastructure, making vehicle travel times dependent on safety separation, intersection access, downstream blocking, and station contention. We study this problem in overhead hoist transport (OHT) systems, a representative ceiling-mounted material-handling system used in semiconductor fabs. Static shortest-path routing cannot account for these time-varying traffic costs, whereas tabular Q-routing adapts online but learns each destination–node–action value independently, limiting information sharing across sparsely visited routing contexts and making startup behavior sensitive to inaccurate value estimates. We propose Neural Double Q-routing, which replaces destination-indexed tables with a shared state–action value network. The network is warm-started through return-to-go regression on mixed simulator-generated routing trajectories and then refined online using Double-Q updates, local congestion correction, and event-stratified structured replay. Across nine matched fleet-size–arrival-rate settings with 100, 150, and 200 OHTs, the proposed framework reduces mean completion time relative to tabular Double Q-routing by $0.8\%$–$8.8\%$. It achieves the lowest mean completion time among all compared methods in the six 150- and 200-OHT settings, whereas Dijkstra remains best in the three 100-OHT settings. Completed-task counts remain within $1\%$ of tabular Double Q-routing in eight of nine settings, and 95th-percentile completion time decreases in eight settings. In two matched startup scenarios, offline initialization increases the number of completed tasks by up to $23\%$ and reduces tail completion time by up to $15\%$.

By focusing on algorithmic stability as a means of establishing out-of-sample bounds, we provide a system-theoretic interpretation of generalization in learning-enabled dynamical systems arising in data-driven optimization and feedback control approximation. Given two neighboring datasets, we specifically model sample replacement as an exogenous disturbance acting on a sensitivity system, while the incremental behavior of the data-dependent operator is encoded through an integral quadratic constraint. By relying on dissipativity arguments, we establish a matrix inequality-based certificate and a uniform stability bound that separates the one-sample sensitivity of the learned operator, and an algorithm-dependent dynamical gain. The latter can then be optimized, offering a tractable tool for certifying and comparing generalization capabilities of learning dynamics. We show that our results recover classical ones for gradient descent, apply naturally to momentum-based methods such as heavy-ball and Nesterov acceleration, and extend to data-driven control.

Healthcare workforce scheduling is an NP-hard optimization problem requiring simultaneous satisfaction of labor regulations, coverage requirements, employee preferences and cost objectives. Existing approaches (genetic algorithms, integer programming, constraint programming) model 6-12 constraints at shift-level granularity and cannot guarantee regulatory compliance. They also lack support for multi-role, multi-skill heterogeneity, mandatory break scheduling with midpoint control, acuity-weighted workload equity, sub-shift granularity, inter-week stability, and cross-midnight shifts. This paper presents CP-SAT: a Constraint Programming formulation for multi-role, multi-skill healthcare scheduling. CP-SAT enforces 14 hard constraints guaranteeing zero regulatory violations, while optimizing 15 soft objectives via a unified weighted penalty function. Contributions include a shift-window decomposition enabling break scheduling with centrality control, acuity-weighted workload equity, multi-granularity resolution from 15 minutes to 1 day, inter-week stability, and grid-offset preprocessing mapping cross-midnight shifts into a single scheduling day without solver changes. CP-SAT is evaluated on 18 instances: five synthetic hospital units (10-33 nurses), 10 INRC-II benchmarks (5-80 nurses, up to 8-week horizons) and 3 NRP-23 compatible instances (10-25 nurses) with cross-midnight Night shifts. Results: zero hard-constraint violations across all 18 instances by construction; proven optimality on INRC-II n005w4 (objective 118, gap 0.0%, 104s); feasible schedules scaling to 179,800 variables and 351,425 constraints (80 nurses); service quality improved 50-67% over MOGA; and model size scaling near-linearly at approximately 4,400 variables per employee. The formulation enforces 29 total constraints (14 hard, 15 soft), nearly three times the industry average.

Popular adaptive stochastic gradient descent (SGD) methods to train artificial intelligence (AI) systems include the RMSprop, the Adam, and the AdamW optimizers, where the adaptivity parts in Adam and AdamW basically just coincide with RMSprop. Such adaptive methods involve several hyperparameters including the regularization parameter $\epsilon$ (which ensures that one does not divide by 0 and is often chosen to be very close to zero such as $10^{-8}$ in PyTorch by default) and the second moment decay parameter $\beta$ (which is often chosen to be very close to $1$ such as 0.99 (RMSprop) and 0.999 (Adam and AdamW) in PyTorch by default). Despite the high relevance of such methods, it remains an open research problem to provide error estimates for such methods with the error constants being not exploding but uniformly bounded with the respect to the hyperparameters, even in the situation of convex stochastic optimization problems. It is the key contribution of this work to essentially solve this problem for RMSprop. Specifically, we bound the expectation of the stopped evaluation of the objective function at the RMSprop process from above by the sum of an initialization term that decays exponentially in the training time, a stochastic approximation remainder of order $\gamma_n$, and a memory error of order $( 1 - \beta)^2$ with the error constants being uniformly controlled over all admissible choices of the step sizes, the second moment decay parameter $\beta$ and the regularization parameter $\epsilon\in[0,1]$ (also covering $\epsilon=0$). Our non-asymptotic error estimates hold not just for all sufficiently large n but hold for every gradient step $n=1,2,3,...$ with all error constants being explicitly specified. The key innovative new feature in the proof of our analysis are suitable inverse moment bounds for the second moment process in RMSprop.