Learning to Solve Differential Equation Constrained Optimization Problems

## Abstract

Abstract Differential equations (DE) constrained optimization plays a critical role in numerous scientific and engineering fields, including energy systems, aerospace engineering, ecology, and finance, where optimal configurations or control strategies must be determined for systems governed by ordinary or stochastic differential equations. Despite its significance, the computational challenges associated with these problems have limited their practical use. To address these limitations, this paper introduces
a learning-based approach to DE-constrained optimization
that combines techniques from proxy optimization [ 1 ] and neural differential equations [ 2 ] . The proposed approach uses a dual-network architecture, with one approximating the control strategies, focusing on steady-state constraints, and another solving the associated DEs. This combination enables the approximation of optimal strategies while accounting for dynamic constraints in near real-time.
Experiments across problems in energy optimization and finance modeling show that this method provides full compliance with dynamic constraints and it produces results up to 25 times more precise than other methods which do not explicitly model the system’s dynamic equations.

## 1 Introduction

In a wide array of scientific and engineering applications, differential equations (DEs) serve as a fundamental tool to model dynamic phenomena where precise predictions and optimal control are crucial. These applications range from energy optimization, where generator dynamics are required to assess system stability, to aerospace engineering, relying on trajectory optimization, and finance, where asset price prediction hinges on stochastic processes.
Central to these applications is the optimization of systems constrained by Ordinary (ODEs) or Stochastic (SDEs) Differential Equations, referred to as *DE-constrained optimization problems*. These problems entail not only solving the DEs but also optimizing decision variables subject to the dynamics dictated by these equations.

This dual requirement, however, poses significant computational challenges. Traditional approaches, such as shooting methods [], collocation methods [], and discretization techniques [], are known to struggle with scalability and precision, especially on high-dimensional and nonlinear systems, which are of interest in this paper.
To address these challenges, this paper introduces a novel learning-based DE-optimization proxy that integrates advancements from two key methodologies: *proxy optimizers* [] and *neural differential equations* (neural-DEs) []. In our approach, a neural network serves as a proxy optimizer, approximating solutions to the decision problem while simultaneously leveraging another neural network to solve the underlying DEs. This novel dual-network architecture exploits a primal-dual method to ensure that both the dynamics dictated by the DEs and the optimization objectives are concurrently learned and respected. Importantly, this integration allows for end-to-end differentiation
enabling efficient gradient-based optimization.

The proposed method is validated across two domains: energy systems and financial modeling and optimization. The experimental results show the ability to directly handle DE in the optimization surrogate, which allows our method to produce solutions that are up to 25 times more precise than
standard proxy optimization techniques, while also adhering to system dynamics. *This precision improvement is important and opens new avenues for research and application in fields that demand high-fidelity dynamic modeling and optimal decision-making at low computational budgets*.

Contributions.
The paper makes the following contributions:
(1) It introduces a novel learning-based method to efficiently approximate solutions of DE-constrained optimization problems. Our approach is a unique integration of neural-DE models to capture the system dynamics and proxy optimizers to approximate the problem
decision variables. These components are sinergistically integrated into model training via a primal-dual method.
(2) It empirically demonstrates the importance of incorporating the system dynamics into an optimization framework, by showing that proxy optimizers methods who neglect these dynamic behaviors systematically violate the DE requirements.
(3) It shows that capturing the system dynamics with neural differential surrogate models, leads up to 25 times higher solution quality compared to other learning-based approaches capturing the dynamics such as PINN [] or LSTM [].

## 2 Related works

A key technique in process control is Model Predictive Control (MPC) [], which predicts future states over a time horizon to select optimal control actions that minimize a cost function while satisfying constraints. However, the DE-optimization problems central to this work involve a parametrization of the state variable equations themselves, controlled by the optimization problem, rendering traditional MPC approaches unsuitable. Surrogate models have been developed for reduced-order adjoint solutions [], sensitivity PDEs [], and learning DE operators [] to address these challenges more effectively.

In recent years, there has been growing interest in leveraging neural network architectures to approximate solutions for challenging constrained optimization problems. Termed *proxy optimizers*, these methods create fast surrogate models by learning mappings from optimization problem parameters to optimal solutions []. Some approaches use supervised learning, requiring precomputed optimal decisions and problem parameters for training [], while others adopt self-supervised strategies, relying solely on the problem’s structure and parameter instances []. A key challenge in this setting is ensuring that these learned solutions satisfy optimization constraints, often addressed by the penalty-based methods [] or through implicit layers that incorporate constraints within the model architecture []. Residual constraint violations can also be corrected post-inference via efficient projection techniques [].

Additionally, various methods have been developed to approximate solutions of differential equations with neural networks. In particular, Physics-Informed Neural Networks (PINNs) [] encode and solve differential equations within their architecture, integrating physical laws as prior knowledge. However, they are notoriously hard to train [] and suffer from limited generalization capabilities []. Another approach is neural differential equations [], which parameterize the hidden state derivatives in a neural network to model system dynamics. These models often capture the underlying dynamics with high fidelity using only data observations.

Our work builds on these areas for surrogate modeling and introduces a learning-based approach to solve, for the first time to our knowledge, DE-constrained optimization problems in near real-time.

## 3 Settings and Goals

Consider an optimization problem constrained by a system of ordinary differential equations(^1^11To ease notation the paper focuses on ODEs, and refers the reader to Appendix A for an extension to SDE.):

$$ $\displaystyle\operatorname*{Minimize}_{{\color[rgb]{0,0,1}\definecolor[named]{ pgfstrokecolor}{rgb}{0,0,1}\bm{u}}}$ $\displaystyle\;\;\overbrace{L(\bm{{{\color[rgb]{0,0,1}\definecolor[named]{ pgfstrokecolor}{rgb}{0,0,1}u}}},\bm{{{\color[rgb]{.75,.5,.25}\definecolor[ named]{pgfstrokecolor}{rgb}{.75,.5,.25}y}}}(T))\,+\int_{t=0}^{\top}\Phi(\bm{{{ \color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}u}}},\bm{{{ \color[rgb]{.75,.5,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.75,.5,.25}y}} }(t),t)\,dt}^{{\cal J}(\bm{{{\color[rgb]{0,0,1}\definecolor[named]{ pgfstrokecolor}{rgb}{0,0,1}u}}},\bm{{{\color[rgb]{.75,.5,.25}\definecolor[ named]{pgfstrokecolor}{rgb}{.75,.5,.25}y}}}(t))}$ (1a) s.t. $\displaystyle\;\;d{\bm{{{\color[rgb]{.75,.5,.25}\definecolor[named]{ pgfstrokecolor}{rgb}{.75,.5,.25}y}}}}(t)=\bm{F}(\bm{{{\color[rgb]{0,0,1} \definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}u}}},\bm{{{\color[rgb]{ .75,.5,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.75,.5,.25}y}}}(t),t)dt$ (1b) $\displaystyle\;\;\bm{{{\color[rgb]{.75,.5,.25}\definecolor[named]{ pgfstrokecolor}{rgb}{.75,.5,.25}y}}}(0)=\bm{I}(\bm{{{\color[rgb]{0,0,1} \definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}u}}})$ (1c) $\displaystyle\;\;\bm{g}(\bm{{{\color[rgb]{0,0,1}\definecolor[named]{ pgfstrokecolor}{rgb}{0,0,1}u}}},\bm{{{\color[rgb]{.75,.5,.25}\definecolor[ named]{pgfstrokecolor}{rgb}{.75,.5,.25}y}}}(t))\leq 0;\;\;\bm{h}(\bm{{{\color[ rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}u}}},\bm{{{\color[ rgb]{.75,.5,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.75,.5,.25}y}}}(t))=0,$ (1d) $$

where $\bm{{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}u}}}=({
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}u}}_{1},
\ldots,{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}u}}_
{n})$ represents the vector of decision variables and $\bm{{{\color[rgb]{.75,.5,.25}\definecolor[named]{pgfstrokecolor}{rgb}{
.75,.5,.25}y}}}(t)=({{\color[rgb]{.75,.5,.25}\definecolor[named]{
pgfstrokecolor}{rgb}{.75,.5,.25}y}}_{1}(t),\ldots,{{\color[rgb]{.75,.5,.25}
\definecolor[named]{pgfstrokecolor}{rgb}{.75,.5,.25}y}}_{m}(t))$ denotes the state variables, each governed by a differential equation $d{{\color[rgb]{.75,.5,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.75,.5,.25}
y}}_{i}(t)=F_{i}(\bm{{{\color[rgb]{.75,.5,.25}\definecolor[named]{
pgfstrokecolor}{rgb}{.75,.5,.25}y}}}(t),\bm{{{\color[rgb]{0,0,1}\definecolor[
named]{pgfstrokecolor}{rgb}{0,0,1}u}}},t)dt$. Here each $F_{i}$ describes the dynamic behavior of the system through ODEs. The set of all such ODEs is captured by $\bm{F}$, as defined in Constraint (1b). Note that these DEs are parametrized by decision variables ${\bm{{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}u}}}}$, rendering the coupling between the control strategy and the system’s dynamic response highly interdependent. The objective function ${\cal J}$ (1a) aims to minimize a combination of the running cost $\Phi$, which varies with the state and decision variables over time, and the terminal cost $L$, which depends on the final state $\bm{{{\color[rgb]{.75,.5,.25}\definecolor[named]{pgfstrokecolor}{rgb}{
.75,.5,.25}y}}}(T)$ and the decision variables $\bm{{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}u}}}$. The time horizon $T$ defines the period over which the optimization takes place. Constraint (1c) sets the initial conditions for the state variables based on the decision variables $\bm{{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}u}}}$. Additional constraints (1d) enforce sets of inequality and equality conditions on the state and decision variables, ensuring that the system constraints are met throughout the decision process.

Figure: Figure 1: Decision variables ${{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}u}}$ represent generators outputs, which are influenced by state variables ${{\color[rgb]{.75,.5,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.75,.5,.25}y}}$ describing rotor angles and speed.
Refer to caption: x1.png

Energy system example.
For example, in the context of power system optimization, decision variables $\bm{{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}u}}}$ capture generators’ power outputs, and state variables $\bm{{{\color[rgb]{.75,.5,.25}\definecolor[named]{pgfstrokecolor}{rgb}{
.75,.5,.25}y}}}(t)$ describe generator rotor angles and speeds, which are key for system stability. The system dynamics in $\bm{F}$, capture the electro-mechanical interactions in the power network,
and their initial conditions, as determined by (1c), are set based on the decision variables. The objective function ${\cal J}$ aims to minimize immediate operational costs like fuel consumption ($\Phi$) and address long-term costs ($L$) over a specific time horizon $T$. Optimizing the generator outputs is finally subject to engineering operational limits and physics constraints (Constraints (1d) , e.g. Ohm’s and Kirkhoff’s laws). An illustration is provided in Figure 1 and the problem description in Appendix B.

### Challenges

While being fundamental for many applications, Problem (1) presents three key
challenges:

- 1.
Finding optimal solutions to Problem (1) is computationally intractable. Even without the differential equation constraints, the decision version of
the problem alone is NP-hard in general.
- 2.
Achieving high-quality approximations of the system dynamics (Equations 1b) and (1c) in near real-time, poses the second significant challenge. The high dimensionality and non-linearity of these dynamics further complicate the task.
- 3.
Finally, the integration of the system dynamics into the decision-making process for solving Problem (1) poses another challenge.
Indeed, including differential equations (1b) in the optimization framework renders traditional numerical methods impractical for real-time applications.

The next section focuses on providing a solution to each of these challenges.

## 4 DE-Optimization Proxy

To address the challenges outlined above, the paper introduces *DE-Optimization Proxy* (DE-OP): a fully differentiable DE-optimization surrogate. In a nutshell, DE-OP defines a dual-network architecture, where one neural network, named $\mathcal{F}_{\omega}$, approximates the *optimal* decision variables $\bm{u}^{*}$, and another, denoted as $\mathcal{N}_{\theta}$, approximates the associated state variables $\bm{y}(t)$, based on the concept of neural differential equations. Here $\omega$ and $\theta$ represents the models’ $\mathcal{F}_{\omega}$ and $\mathcal{N}_{\theta}$ trainable parameters, respectively. An illustration of the overall framework and the resulting interaction of the dual network is provided in Figure 2.
The subsequent discussion first describes these two components individually and then shows their integration by exploiting a primal-dual learning framework.

Optimizing over distribution of instances.
To enable a learnable mechanism for addressing DE-constrained optimization, DE-OP operates over a distribution $\Pi$ of problem instances induced by problem parameters $\bm{\zeta}$, and aims to train a model across this distribution. The learning framework takes problem parameters $\bm{\zeta}$ as inputs and generates outputs $\hat{\bm{u}}$, representing an approximation of the optimal decision variables while adhering to the constraint functions (1b)–(1d). With reference to (1a), the formal learning objective is:

$$ $\displaystyle\operatorname*{Minimize}_{\omega,\theta}\mathbb{E}_{\bm{\zeta} \sim\Pi}\left[{\cal J}\left(\mathcal{F}_{\omega}(\bm{\zeta}),\mathcal{N}_{ \theta}(\mathcal{F}_{\omega}(\bm{\zeta}),t);\bm{\zeta}\right)\right]$ (2a) s.t. ( 1b )–( 1d ), (2b) $$

where $\hat{\bm{u}}\!=\!\mathcal{F}_{\omega}(\bm{\zeta})$ and the state variables estimate is denoted as $\hat{\bm{y}}(t)\!=\!\mathcal{N}_{\theta}(\mathcal{F}_{\omega}(\bm{\zeta}),t)$. Here, $\bm{\zeta}$ parameterizes each problem instance, representing constants such as customer demands in the power system example. Although the problem structure remains consistent across instances, each one involves a distinct decision problem, leading to unique state variable trajectories.
Given the complexity of solving Problem (2), the goal is to develop a fast and accurate neural DE optimization surrogate. This approach uses the concept of *proxy optimizers*, which is detailed next.

Figure: Figure 2: DE-OP uses a dual network architecture consisting of a proxy optimization model $\mathcal{F}_{\omega}$ to estimate the decision variables $\hat{\bm{u}}$ and a neural-DE model $\mathcal{N}_{\theta}$ to estimate the state-variables $\hat{\bm{y}}(t)$, with the objective function $\mathcal{J}(\mathcal{F}_{\omega}(\bm{\zeta}),\mathcal{N}_{\theta}\left( \mathcal{F}_{\omega}(\bm{\zeta}),t\right);\bm{\zeta})$ capturing the overall loss.
Refer to caption: x2.png

### 4.1 Neural Optimization Surrogate

The first objective is to establish a neural network-based mapping $F:\Pi\to\mathcal{U}$
that transforms parameters $\bm{\zeta}\sim\Pi$ from a DE-constrained optimization problem (1) into optimal decisions
$\bm{u}^{\star}(\bm{\zeta})\in\mathcal{U}$, operating under the restriction that $T=0$, (i.e. the dynamics of the system are absent). Practically, this mapping is modeled as a neural network $\mathcal{F}_{\omega}$which learns to predict the optimal decision variables from the problem parameters.
The model training uses a dataset $\mathcal{D}=\{(\bm{\zeta}_{i},\bm{u}_{i}^{\star})\}_{i=1}^{N}$, of $N$ samples, with each sample $(\bm{\zeta}_{i},\bm{u}_{i}^{\star})$ including the observed problem parameters $\bm{\zeta}_{i}$ and the corresponding (steady-state) optimal decision variables $\bm{u}_{i}^{\star}$.
The training objective is to refine $\mathcal{F}_{\omega}$ such that it closely approximates the ideal mapping $F$.
Several approaches have been proposed to build such surrogate optimization solvers, many of which leverage mathematical optimization principles [, ] and implicit layers [], to encourage or ensure constraint satisfaction (see Appendix C.1 for an in-depth discussion of such methods).
Once trained, the model ${\cal F}_{\omega}$ can be used to generate near-optimal solutions at low inference times. We leverage this idea to learn the mapping $F$.
As shown in Figure 2, the estimated optimal decisions $\hat{\bm{u}}=\mathcal{F}_{\omega}(\bm{\zeta})$ are then feeded to a neural-DE model, which will be discussed next.

### 4.2 Neural Estimation of the State Variables

The second objective of DE-OP involves to efficiently capture the system dynamics of DE-constrained optimization problems. This is achieved by developing neural DE models $\mathcal{N}_{\theta}$ to learn solutions of a *parametric family* of differential equations []. *Since these DE-surrogates are fully differentiable, they are particularly suitable for integration with the optimization surrogate introduced in the previous section*, aligning with our goal defined in Equation (2).

The optimization proxy estimated solutions $\hat{\bm{u}}$ determine the initial state variables $\bm{y}(0)$ through the function $\bm{I}$: $\hat{\bm{y}}(0)=\bm{I}(\hat{\bm{u}})$. As shown in Figure 2, given a solution $\hat{\bm{u}}$,
a neural-DE model $\mathcal{N}_{\theta}$ generates an estimate $\hat{\bm{y}}$ of the state variables that satisfies:

$$ $\displaystyle d{\hat{\bm{y}}}(t)$ $\displaystyle=\mathcal{N}_{\theta}(\hat{\bm{u}},t)dt$ (3a) $\displaystyle\hat{\bm{y}}(0)$ $\displaystyle=\bm{I}(\hat{\bm{u}}).$ (3b) $$

The remainder of this section details the methods by which the state variables are precisely estimated using a Neural Differential Equation model.

Model initialization and training. Given that the neural-DE model $\mathcal{N}_{\theta}$ takes as input the estimated decision variables $\hat{\bm{u}}$ provided by the DE-OP’s optimization proxy $\mathcal{F}_{\omega}$, where $\hat{\bm{u}}=\mathcal{F}_{\omega}(\bm{\zeta})$, it is practical to initialize (or hot-start) the neural-DE model effectively. To achieve this, we construct a dataset $\mathcal{D}=\{(\bm{u}_{i}^{\prime},\bm{y}_{i}^{\star}(t))\}_{i=1}^{N}$, where each $\bm{u}_{i}^{\prime}$ is a near-optimal decision sampled within the bounds specified by constraints $\bm{g}$ (e.g., $u_{j}^{\prime}\sim\mathcal{U}(a_{j},b_{j})$ if $u_{j}$’s bound is defined by $a_{j}$ and $b_{j}$). The corresponding state trajectories $\bm{y}_{i}^{\star}(t)$ are obtained by numerically solving the differential equations with initial condition $\bm{y}(0)=\bm{I}(\bm{u}_{i}^{\prime})$.
The neural-DE model is trained by minimizing the loss:

$$ $\operatorname*{Minimize}_{\theta}\mathbb{E}_{(\bm{x},\bm{y})\sim\mathcal{D}} \left[\left\|\mathcal{N}_{\theta}(\bm{x},t)-\bm{y}(t)\right\|^{2}\right],$ (4) $$

where $\bm{x}=\bm{u}^{\prime}$ and $\hat{\bm{y}}(t)=\mathcal{N}_{\theta}(\bm{x},t)$.

Since $\mathcal{F}_{\omega}$ approximates $\bm{u}^{\star}$, it may introduce errors. To mitigate this, the dataset $\mathcal{D}$ is constructed using near-optimal decisions $\bm{u}^{\prime}$ sampled from the feasible bounds $\bm{g}$, ensuring that $\mathcal{N}_{\theta}$ is trained on a distribution $\Pi_{\bm{u}^{\prime}}\approx\Pi_{\hat{\bm{u}}}$. This approach does not require exact optimal decisions and assumes only small estimation errors from $\mathcal{F}_{\omega}$, which is typically valid in practice.

Once trained, the neural-DE model $\mathcal{N}_{\theta}$ can accurately estimate state variables $\hat{\bm{y}}(t)$ for decisions $\hat{\bm{u}}$ produced by $\mathcal{F}_{\omega}$. This integration enables end-to-end training of the DE-OP framework, ensuring that both decision and state variables are optimized cohesively.
While other learning-based methods, such as Physics-Informed Neural Networks (PINN) [] and Long Short-Term Memory (LSTM) models [], can be used to efficiently estimate the state variables, due to a generalization bias (PINN), and lack of dynamic equations modeling (LSTM), they produce substantially less precise predictions than the neural-DE models. For a comparison of neural-DE models with these alternative approaches across the experimental tasks described in Section 5, please refer to Appendices D.1 and E.1.

### 4.3 Handling Static and Dynamics Constraints Jointly

To integrate the neural-DE models within the decision process, this paper proposes a Lagrangian Dual (LD) learning approach, which is inspired by the generalized augmented Lagrangian relaxation technique [] adopted in classic optimization.
In Lagrangian relaxation, some or all the problem constraints are relaxed into the objective function using Lagrangian multipliers to capture the penalty induced by violating them.
The proposed formulation leverages Lagrangian duality to integrate trainable and weighted regularization terms that encapsulates both the static and dynamic constraints violations.
When all the constraints are relaxed, the violation-based Lagrangian relaxation of problem (1) defines as

$$ $\operatorname*{Minimize}_{\bm{u}}\mathcal{J}({\bm{u}},{\bm{y}}(t))+\bm{\lambda }^{\top}_{h^{\prime}}|\bm{h}^{\prime}({\bm{u}},{\bm{y}}(t))|+\bm{\lambda}^{ \top}_{g}\max(0,\bm{g}({\bm{u}},{\bm{y}}(t))),$ $$

where $\mathcal{J}$, $\bm{g}$ are defined in (1a), and (1d), respectively, and $\bm{h}^{\prime}$ is defined as follows,

$$ $\bm{h}^{\prime}(\bm{u},\bm{y}(t))=\begin{bmatrix}d\bm{y}(t)-\bm{F}(\bm{u},\bm{ y}(t),t)dt\\ \bm{y}(0)-\bm{I}(\bm{u})\\ \bm{h}(\bm{u},\bm{y}(t)).\end{bmatrix}$ (5) $$

It denotes the set of *all* equality constraints of problem (1), thus extending the constraints $\bm{h}$ in (1d), with the system dynamics (1b) and the initial conditions equations (1c) written in an implicit form as above.
Therein, $\bm{\lambda}_{h^{\prime}}$ and $\bm{\lambda}_{g}$ are the vectors of Lagrange multipliers associated with functions $\bm{h}^{\prime}$ and $\bm{g}$, e.g. ${\lambda}^{i}_{h^{\prime}},{\lambda}^{j}_{g}$ are associated with the $i$-th equality $h^{\prime}_{i}$ in $\bm{h}^{\prime}$ and $j$-th inequality $g_{j}$ in $\bm{g}$, respectively. The key advantage of expressing the system dynamics (1b) and initial conditions (1c) in the same implicit form as the equality constraints $\bm{h}$, (as shown in (5)), is that the system dynamics can be treated in the same manner as the constraint functions $\bm{h}$. This enables us to satisfy the system dynamics and the static set of constraints ensuring that they are incorporated seamlessly into the optimization process.

Figure: Algorithm 1 Primal Dual Learning for DE-Constrained Optimization

The proposed primal-dual learning method uses an iterative approach to find good values of the *primal* $\omega$, $\theta$ and dual $\bm{\lambda}_{h^{\prime}},\bm{\lambda}_{g}$ variables; it uses an augmented modified Lagrangian as a loss function to train the prediction $\hat{\bm{u}}$, $\hat{\bm{y}}(t)$ as employed

$$ $\begin{split}\mathcal{L}^{\text{DE-OP}}(\hat{\bm{u}},\bm{u}^{\star},\hat{\bm{y }}(t))&=\|\hat{\bm{u}}-\bm{u}^{\star}\|^{2}+\bm{\lambda}^{\top}_{h^{\prime}}| \bm{h}^{\prime}(\hat{\bm{u}},\hat{\bm{y}}(t))|+\bm{\lambda}^{\top}_{g}\max(0, \bm{g}(\hat{\bm{u}},\hat{\bm{y}}(t))),\end{split}$ (6) $$

where $\|\hat{\bm{u}}-\bm{u}^{\star}\|^{2}$ represents the decision error, with respect to *steady-state* optimal decision $\bm{u}^{\star}$, while $\bm{\lambda}^{\top}_{h^{\prime}}|\bm{h}^{\prime}(\hat{\bm{u}},\hat{\bm{y}}(t))|$ and $\bm{\lambda}^{\top}_{g}\max(0,\bm{g}(\hat{\bm{u}},\hat{\bm{y}}(t)))$ measures the constraint violations incurred by prediction $\hat{\bm{u}}=\mathcal{F}_{\omega}(\bm{\zeta})$ and $\hat{\bm{y}}(t)=\mathcal{N}_{\theta}(\mathcal{F}_{\omega}(\bm{\zeta}),t)$. To ease notation, in Equation (6) the dependency from parameters $\bm{\zeta}$ is omitted.
At iteration $k+1$, finding the optimal parameters $\omega,\theta$ requires solving

$$ $\omega^{k+1},\theta^{k+1}=\arg\min_{\omega,\theta}\mathbb{E}_{(\bm{\zeta},\bm{ u}^{\star})\sim\mathcal{D}}\left[\mathcal{L}^{\text{DE-OP}}\left(\mathcal{F}_{ \omega}^{\bm{\lambda}^{k}}(\bm{\zeta}),\bm{u}^{\star},\mathcal{N}_{\theta}^{ \bm{\lambda}^{k}}\left(\mathcal{F}_{\omega}^{\bm{\lambda}^{k}}(\bm{\zeta}),t \right)\right)\right],$ $$

where $\mathcal{F}_{\omega}^{\lambda^{k}}$ and $\mathcal{N}_{\theta}^{\lambda^{k}}$ denote the DE-OP’s optimization and predictor models $\mathcal{F}_{\omega}$ and $\mathcal{N}_{\theta}$, at iteration $k$, with $\bm{\lambda}^{k}=[\bm{\lambda}^{k}_{h^{\prime}}\ \bm{\lambda}^{k}_{g}]^{\top}$. This step is approximated using a stochastic gradient descent method

$$ $\displaystyle\omega^{k+1}=\omega^{k}-\eta\nabla_{\omega}\mathcal{L}^{\text{DE- OP}}\left(\mathcal{F}_{\omega^{k}}^{\bm{\lambda}^{k}}(\bm{\zeta}),\bm{u}^{ \star},\mathcal{N}_{\theta^{k}}^{\bm{\lambda}^{k}}\left(\mathcal{F}_{\omega^{k }}^{\bm{\lambda}^{k}}(\bm{\zeta}),t\right)\right)$ $\displaystyle\theta^{k+1}=\theta^{k}-\eta\nabla_{\theta}\mathcal{L}^{\text{DE- OP}}\left(\mathcal{F}_{\omega^{k}}^{\bm{\lambda}^{k}}(\bm{\zeta}),\bm{u}^{ \star},\mathcal{N}_{\theta^{k}}^{\bm{\lambda}^{k}}\left(\mathcal{F}_{\omega^{k }}^{\bm{\lambda}^{k}}(\bm{\zeta}),t\right)\right),$ $$

where $\eta$ denotes the learning rate and $\nabla_{\omega}\mathcal{L}$ and $\nabla_{\theta}\mathcal{L}$ represent the gradients of the loss function $\mathcal{L}$ with respect to the parameters $\omega$ and $\theta$, respectively, at the current iteration $k$. Importantly, this step does not recomputes the training parameters from scratch, but updates the weights $\omega$, $\theta$ based on their value at the previous iteration.
Finally, the Lagrange multipliers are updated as

$$ $\displaystyle\bm{\lambda}^{k+1}_{h^{\prime}}=\bm{\lambda}^{k}_{h^{\prime}}+ \rho|\bm{h}^{\prime}(\hat{\bm{u}},\hat{\bm{y}}(t))|$ $\displaystyle\bm{\lambda}^{k+1}_{g}=\bm{\lambda}^{k}_{g}+\rho\max\left(0,\bm{g }(\hat{\bm{u}},\hat{\bm{y}}(t))\right),\$ $$

where $\rho$ denotes the Lagrange step size.
The overall training scheme is presented in Algorithm 1. It takes as input the training dataset $\mathcal{D}=\{(\bm{\zeta}_{i},\bm{u}_{i}^{\star})\}_{i=1}^{N}$, the learning rate $\eta>0$, and the Lagrange step size $\rho>0$. The Lagrange multipliers are initialized in line $2$.
As shown in Figure 2, for each sample in the dataset (line $4$), the DE-OP’s optimization model $\mathcal{F}_{\omega^{k}}$ computes the predicted decisions $\hat{\bm{u}}_{i}$, while $\mathcal{N}_{\theta}$ computes an estimate of the state variables $\hat{\bm{y}}_{i}(t)$ (line $5$).
The loss function $\mathcal{L}^{\text{DE-OP}}$ is computed (line $6$) incorporating both the objective and the constraints and using the predicted values $\hat{\bm{u}}$, $\hat{\bm{y}}(t)$ and the Lagrange multipliers $\bm{\lambda}^{k}_{h^{\prime}}$ and $\bm{\lambda}_{g}^{k}$.The weights $\omega,\theta$ of the DE-OP models $\mathcal{F}_{\omega},\mathcal{N}_{\theta}$ are then updated using stochastic gradient descent (SGD) (line $7$).
Finally, at the end of the epoch, the Lagrange multipliers are updated based on the respective constraint violations (line $8$).

While the DE-OP model training algorithm is described extending the Lagrangian Dual Learning approach [], the flexibility of DE-OP allows to leverage other proxy optimizer methods, such as the self-supervised Primal-Dual Learning [], which could similarly be extended to integrate the system dynamics via neural-DE modeling within the DE-OP framework, as the experiments will show.

## 5 Experimental setting

This section evaluates the DE-OP model financial modeling and energy optimization tasks. Given the absence of other methods capable of meeting the stringent time requirements for solving DE-constrained optimization problems, we compare several proxy optimizer methods as baselines. However, these baselines focus on the “steady-state” aspects of the problem by omitting the system dynamic components, such as the objective term $\Phi$ and the system dynamic constraints (1b) and (1c) from Problem (1). They aim to approximate the optimal decision variables $\hat{\bm{u}}$ to $\bm{u}^{*}$ that could be obtained if the system was at a steady-state. We evaluate the Lagrangian Dual approach (LD) of [], which uses a penalty-based method for constraint satisfaction, Deep Constraint Completion and Correction (DC3) from [] that enforces constraint satisfaction through a completion-correction technique, self-supervised learning (PDL) of [] using an augmented Lagrangian loss function, and a method (MSE) that minimizes the mean squared error between the predicted solutions $\hat{\bm{u}}$ and the pre-computed (steady-state) solutions $\bm{u}^{*}$. A comprehensive description of these methods is provided in Appendix C.1.

Furthermore, the comparison includes various learning-based DE-surrogate solvers in place of the network ${\cal N}_{\theta}$ in our framework, including neural-differential equations [], PINNs [], and LSTM networks. The experiments focus on two main aspects: (1) comparing DE-OP with proxy optimizers that capture only the steady-state problems, focusing on the system dynamics violations, and (2) assessing the effectiveness of the various surrogate DE-solver methods.

### 5.1 Dynamic Portfolio Optimization

The classical Markowitz Portfolio Optimization [], described by (9a)-(9c), consists of determining the investment allocations within a portfolio to maximize a balance of return and risk. The paper extends this task by incorporating the stochastic dynamic (9d) of the asset prices, based on a simplified Black-Scholes model []. This model represents a real-world scenario where asset prices fluctuates, and investment decisions are made in advance, such as at the market’s opening, based on the final asset prices forecast.
The task defines as:

$$ $\displaystyle\operatorname*{Minimize}_{\bm{u}}\quad$ $\displaystyle\mathbb{E}\left[-\bm{y}(T)^{\top}\bm{u}+\bm{u}^{\top}\bm{\Sigma} \bm{u}\right]$ (9a) s.t. $\displaystyle\bm{1}^{\top}\bm{u}=1$ (9b) $\displaystyle u_{i}\geq 0\ \;\;\forall i\in[n]$ (9c) $\displaystyle dy_{i}(t)=\mu_{i}y_{i}(t)dt+\sigma_{i}y_{i}(t)dW_{i}(t)\ \;\; \forall i\in[n]$ (9d) $\displaystyle y_{i}(0)={{\zeta}_{i}}\ \;\;\forall i\in[n],$ (9e) $$

where $\bm{y}(t)\in\mathbb{R}^{n}$ represents the asset prices trend and $\bm{y}(T)$ denotes the asset prices at time horizon $T$ in (9a). The asset price dynamics described by (9d) follow a stochastic differential equation with drift $\mu_{i}$, volatility $\sigma_{i}$, and Wiener process $W_{i}(t)$ []. Decisions $\bm{u}\in\mathbb{R}^{n}$ represent fractional portfolio allocations. The objective minimizes $J(\bm{u},\bm{y}(t))=L(\bm{u},\bm{y}(T))$ with $\Phi(\bm{u},\bm{y}(t),t)=0$, balancing risk via covariance matrix $\bm{\Sigma}$ and expected return $\bm{y}(T)^{\top}\bm{u}$.

Datasets and methods.
The drift and volatility factors $\mu_{i}\sim\mathcal{U}(0.5,1)$ and $\sigma_{i}\sim\mathcal{U}(0.05,0.1)$ in (9d) are sampled from uniform distributions. Initial asset prices $\{\zeta_{i}\}_{i=1}^{n}$ are obtained from the Nasdaq database [] to form initial vectors $\{\bm{\zeta}^{j}\}_{j=1}^{10,000}$, split into 80% training, 10% validation, and 10% test sets. Asset price trends $\bm{y}(t)$ and final prices $\bm{y}(T)$ are generated using an SDE solver with Itô integration []. Given $\bm{y}(T)$, the convex solver cvxpy [] computes the optimal decision $\bm{u}^{\star}$ for supervision during training.

We evaluate the role of asset price predictors using three models: a neural-SDE model, an LSTM, and a 2-layer Feed Forward ReLU network. Each model estimates final asset prices $\hat{\bm{y}}(T)$, which inform the DE-OP model $\mathcal{F}_{\omega}$ to estimate optimal decision allocations $\hat{\bm{u}}=\mathcal{F}_{\omega}(\hat{\bm{y}}(T))$. A “static” baseline method uses only proxy optimizers (Lagrangian Dual, DC3, or PDL) to approximate decisions based on initial prices $\bm{y}(0)=\bm{\zeta}_{i}$. Detailed comparisons of these approaches are in Appendix E.2.

#### Results.

A comparison between the DE-OP and the baseline methods is provided in Figure 3. The figure reports experiments for $n=50$ variable, while additional experiments are relegated to the Appendix E.2.
The x-axis categorizes the methods based on the type of asset price predictor used, or the lack of thereof, and reports the average *optimality gap* (in percentage) on the test set, defined as $\frac{|L(\bm{u}^{\star}(\bm{y}(T)),\bm{y}(T))-L(\hat{\bm{u}}(\hat{\bm{y}}(T)),
\bm{y}(T))|}{|L(\bm{u}^{\star}(\bm{y}(T)),\bm{y}(T))|}\times 100$. It measures the sub-optimality of the predicted solutions $\hat{\bm{u}}$ with respect to ground truth final asset price $\bm{y}(T)$, across different methods.

Figure: Figure 3: Average Opt. gap with $n=50$ asset prices.
Refer to caption: x3.png

The figure highlights the substantial performance difference between the dynamic DE-OP models and the static-only proxy optimizer methods, marked in the last column, denoted by ${\cal N}_{\theta}=\emptyset$. These static methods (DC3, LD, and PDL) fail to incorporate asset price dynamics, resulting in notably higher optimality gaps (exceeding 100%). Notably, their predictions can be over twice as suboptimal as the optimal solutions derived from dynamic modeling.

In contrast, DE-OP models that incorporate SDE models significantly outperform static methods. In particular, using neural-SDE predictors to model asset price dynamics results in *much higher* decision quality compared to both LSTM and Feed Forward models. Specifically, DE-OP with neural-SDE achieves the lowest optimality gap at $9.11\%$, successfully capturing the dynamics through an explicit modeling of the asset prices’ governing equations. In contrast, the LSTM model results in a notably higher optimality gap of $21.17\%$, approximately 2.5 times greater than that of DE-OP with neural-SDE, attributable to its lack of dynamic modeling. The Feed Forward model performs significantly worse with an optimality gap of $102.45\%$ for LD, indicating its inability to capture the time-dependent nature of asset pricing data.

The dynamic forecasting results in Appendix E.1 display different levels of precision of the final asset prices predictions among the dynamic predictors considered, which ultimately led to different decision quality. *In particular, DE-OP with a neural-SDE model performs consistently better than the LSTM model and produces up to $25\times$ better decisions (measured in terms of optimality gap) than any static-only proxy optimizer method*. This stark contrast underscores the effectiveness of DE-OP models in leveraging dynamic asset price predictors to improve decision quality.

### 5.2 Stability-Constrained AC-Optimal Power Flow

We now turn on studying a real world problem in power systems, which arises from integration of Synchronous Generator Dynamics with the Alternating Current Optimal Power Flow (AC OPF) problem. The AC OPF, detailed in Appendix Model (1), is foundational in power systems for finding cost-effective generator dispatches that meet demand while complying with physical and engineering constraints. Traditionally addressed as a *steady-state* snapshot, the AC OPF problem requires frequent resolution (e.g., every 10-15 minutes) due to fluctuating loads, posing challenges in maintaining operational continuity and system stability [].
Given the non-convexity, high dimensionality, and computational demands of this problem, proxy optimizers have emerged as a viable alternative to traditional numerical solvers. However, as shown later in this section, existing approaches such as those in [] and [], which focus on the steady-state aspect, fail to address the dynamic system requirements adequately.

The integration of the generator dynamics and related stability constraints into the steady-state AC-OPF formulation leads to the *stability-constrained* AC-OPF problem, which is detailed in Appendix B.3. Here, the decision variables $\bm{u}$, comprising generator power outputs and bus voltages, influence the state variables $\bm{y}(t)$, representing generator rotor angles and speeds. This coupling renders the problem particularly challenging. The objective is to optimize power dispatch costs, while satisfying demand, network constraints, and ensuring system stability as described by (1d).
DE-OP enables, for the first time to our knowledge, the integration of generator dynamics within the optimization process. This integration is key for system stability. Our implementation uses neural-ODE [] models $\mathcal{N}_{\theta}$, assigned individually to each generator. A comparative analysis of neural-ODEs and PINNs for modeling these dynamics is available in Appendix D.1.

Datasets.
DE-OP is evaluated on two key power networks, the WSCC 9 and IEEE 57 bus-systems [], under various operating system conditions.
This assessment benchmarks our proposed method, DE-OP, against three leading proxy optimizer methods for AC-OPF, which operate under a “steady-state” assumption, and thus cannot cope with constraints (1b) and (1c) of the DE-constrained problem. The methods are LD [], DC3 [], and MSE [], introduced in details in Appendix C.1.

As outlined in Model 2 in Appendix B.3, the role of decision variables $\bm{u}$ in influencing both the initial conditions and the governing equations, along with their independence from the time variable $t$, makes Model Predictive Control methods unsuitable for addressing the stability-constrained AC OPF problem. Furthermore, discretizing the DE system through methods like direct collocation [] becomes highly impractical for real-time applications due to the high number of variables and nonlinear system dynamics associated with each generator in the system. In contrast, all methods used for comparison (including our DE-OP) produce estimates of the optimal decision variables within milliseconds, as shown in Appendix, Table 6. Additionally, we note that only fast inference times are of interest in the area of proxy optimizers, as once trained, these methods can be applied to various related but distinct problem instances.

For both the benchmark systems, DE-OP and each proxy optimizer model are trained on a dataset $\mathcal{D}=\{(\bm{\zeta}_{i},\bm{u}_{i}^{\star})\}_{i=1}^{10,000}$, where $\bm{\zeta}_{i}$ representing a load demand and $\bm{u}_{i}^{\star}$ the corresponding optimal, *steady-state* decision. The $j^{th}$ load of the $i$-th sample $\bm{\zeta}_{i}^{j}$ is generated by applying a uniform random perturbation of $\pm 20\%$ to the corresponding nominal load. The (steady-state) AC-OPF models are implemented in Julia and solved using IPOPT []. The dataset uses an 80/10/10 split.
By leveraging the knowledge of the decision variables’ bounds (see Appendix, Model 2, Constraints (11b)-(11c)), each neural-ODE model $\mathcal{N}_{\theta}$ is trained to learn the corresponding generator dynamics on a dataset of near-optimal decisions $\Pi_{u^{\prime}}$ as described in Section 4.2. We refer the reader to Appendix B.3 for further details on the neural-ODE models training.

Figure: Figure 4: WSCC $9$-bus - Percentage of unstable dynamics observed during training. Average of 40 experiments.
Refer to caption: x4.png

Table: Table 1: Average and standard deviation of constraint violations and (steady-state) optimality gap on the IEEE $57$-bus system for different approaches based on $40$ independent runs.

#### Results.

Figures 5 and 5 show the percentage of estimated decisions violating the stability constraints during the first $50$ epochs of training, across all methods and test cases.
On the WSCC 9-bus system (Figure 5) DE-OP learns rapidly to meet the dynamic constraints, which violations approaches zero level after epoch $10$ of training, whereas, all the baseline methods lacking dynamic modeling, consistently produce unstable dynamics. Notably, all the baseline methods tested, e.g., DC3, LD, and MSE, systematically fail to satisfy stability requirements.
In contrast, by integrating generator dynamics within its training model, DE-OP begins to satisfy these requirements early in training, as depicted in both figures (blue curves). DE-OP shows a rapid adjustment on both systems within the first few epochs, and bringing the violations to near zero. In contrast, all baseline methods continue to exhibit 4% to 8% unstable dynamics throughout the training, even for a much higher number of training epochs. As we will show later, these will reflect also in large stability constraints violations, when evaluated on the test set.

Table 1 displays the test-set results for DE-OP and the baseline methods on the IEEE 57-bus system. For a detailed discussion on the results of each method on the WSCC-9 bus system, please see Appendix D and Table 5. These tables reports the following metrics:

- •
Static and Stability Constraint Violations: These are quantified for each test instance. The $j$-th static equality and the $k$-th inequality violation are calculated as $\frac{1}{n_{\text{test}}}\sum_{i=1}^{n_{\text{test}}}|h^{\prime}_{j}(\hat{\bm{
u}}^{i},\hat{\bm{y}}^{i}(t))|$ and $\frac{1}{n_{\text{test}}}\sum_{i=1}^{n_{\text{test}}}\max(0,g_{k}(\hat{\bm{u}}
^{i},\hat{\bm{y}}^{i}(t)))$ respectively, where $n_{\text{test}}$ is the test-set size. Detailed descriptions of the problem constraints can be found in Appendices B.1 and B.3.
- •
Optimality gap (at steady-state): This metric is defined as
$\frac{|L(\bm{u}^{\star}(\zeta),\bm{y}^{\star}(T))-L(\hat{\bm{u}}(\zeta),\hat{
\bm{y}}(T))|}{|L(\bm{u}^{\star}(\zeta),\bm{y}^{\star}(T))|}\times 100$.
It measures the gap incurred by the predictions $\hat{\bm{u}},\hat{\bm{y}}(t)$ against the decisions $\bm{u}^{\star}$ which are computed under the assumption that the generators are in a steady-state condition. This assumption is crucial for evaluating how closely each solution approximates the AC-OPF optimal results, though it *does not necessarily reflect the results relative to the stability-constrained AC-OPF problem*, our main focus, but which is highly intractable. Given the non-linearity of both the dynamics and optimization in the stability-constrained AC-OPF, computing exact optimal decisions $\bm{u}^{\star}$ with traditional methods is not feasible. Consequently, while our method may show slightly higher steady-state optimality gaps, these should not be interpreted in the context of the dynamic problem.

Firstly, note that all methods report comparable static constraint violations and steady-state optimality gaps, with errors within the $10^{-3}$ to $10^{-4}$ range. Although DE-OP exhibits slightly higher steady-state optimality gaps, approximately 0.5% higher than the best performing baseline, it’s important to recall that this metric *does not* reflect the stability-constrained AC-OPF optimality gap, but rather that of the problem addressed by the baselines, placing DE-OP at a seeming disadvantage. The higher objective costs observed with DE-OP is intuitively attributed to a restricted feasible space resulting from the integration of generator stability constraints within the AC OPF problem.

Crucially, all baseline methods systematically fail to meet the stability requirements, often by margins exceeding those observed during training. This illustrates a typical scenario where prediction errors on decisions parametrizing system dynamics have cascading effects on the associated constraint violations. In stark contrast, DE-OP achieves full compliance with stability constraints, reporting zero violations in each instance analyzed.

These findings are important. They highlight the critical role of dynamic requirements in AC-OPF problems for achieving accurate and stable solutions. The results underscore DE-OP’s effectiveness in adjusting potentially unstable set points, as further detailed in Appendix D.2 and demonstrates DE-OP’s effectiveness in ensuring system stability compared to traditional methods that focus on optimality under steady-state assumptions.

## 6 Conclusion

This work was motivated by the efficiency requirements associated with solving differential equations (DE)-constrained optimization problems.
It introduced a novel learning-based framework, DE-OP, which incorporate differential equation constraints into optimization tasks for near real-time application. The approach uses a dual-network architecture, with one approximating the control strategies, focusing on steady-state constraints, and another solving the associated DEs. This architecture exploits a primal-dual method to ensure that both the dynamics dictated by the DEs and the optimization objectives are concurrently learned and respected. This integration allows for end-to-end differentiation enabling efficient gradient-based optimization, and, for the first time to our knowledge, solving DE-constrained optimization problems in near real-time.
Empirical evaluations across financial modeling and energy optimization tasks, illustrated DE-OP’s capability to adeptly address these complex challenges. The results demonstrate not only the effectiveness of our approach but also its broad potential applicability across various scientific and engineering domains where system dynamics are crucial to optimization or control processes.

## 7 Acknowledgments

This work was partially supported by NSF grants EPCN-2242931 and NSF CAREER-2143706. The view and conclusions are those of the authors only.

## Appendix A Stochastic Differential Equation Constrained Optimization

This section extends the problem description (1) presented in Section 3 of the main paper to a stochastic setting.
In presence of stochastic dynamics, the optimization problem constrained by differential equations (1) becomes

$$ $\displaystyle\operatorname*{Minimize}_{\bm{u}}$ $\displaystyle\;\;\mathbb{E}\left[L(\bm{u},\bm{y}(T))\,+\int_{t=0}^{T}\Phi(\bm{ u},\bm{y}(t),t)\,dt\right]$ (10a) s.t. $\displaystyle\;\;d\bm{y}(t)=\bm{F}(\bm{u},\bm{y}(t),t)dt+\bm{G}(\bm{u},\bm{y}( t),t)d\bm{W}(t)$ (10b) $\displaystyle\;\;\bm{y}(0)=\bm{I}(\bm{u})$ (10c) $\displaystyle\;\;\bm{g}(\bm{u},\bm{y}(t))\leq 0$ (10d) $\displaystyle\;\;\bm{h}(\bm{u},\bm{y}(t))=0.$ (10e) $$

The SDE constrained optimization problem (10) involves determining the optimal decision variables $\bm{u}=(u_{1},\ldots,u_{n})$ in a system where the state variables $\bm{y}(t)=(y_{1}(t),\ldots,y_{m}(t))$ evolve according to stochastic dynamics (10a) and initial conditions dictated by (10b). Each state variable $y_{i}(t)$ is governed by a stochastic differential equation $d\bm{y}_{i}(t)=F_{i}(\bm{y}(t),\bm{u},t)dt+G_{i}(\bm{y}(t),\bm{u},t)dW_{i}(t)$, where $F_{i}$ represents the deterministic part of the dynamics, and $G_{i}$ captures the stochastic component, where $W_{i}(t)$ is a Wiener process. The set of all such equations is described by $\bm{F}$ and $\bm{G}$, as defined by (10b). The initial condition for the state variables is set by constraints (10c), where $\bm{y}(0)=\bm{I}(\bm{u})$ defines the starting state based on the control variables $\bm{u}$. Constraints (10d) and (10e) enforce inequality and equality constraints, respectively, on the state and control variables, ensuring that the system behaves within specified bounds throughout the decision process.

The objective (10a) is to minimize the expected value of a combination of the running cost $\Phi(\bm{u},\bm{y}(t),t)$, which varies with the state and decision variables over time, and the terminal cost $L(\bm{u},\bm{y}(T))$, which depends on the final state $\bm{y}(T)$ and the decision variables $\bm{u}$. The optimization is performed over a time horizon $T$, which defines the period during which the decision-making process occurs.

## Appendix B Stability Constrained AC-Optimal Power Flow

This section describes the stability constrained AC-Optimal Power Flow problem; it first introduces the AC-Optimal Power Flow problem and the synchronous generator dynamics, to eventually integrate these two components to form the stability constrained AC-Optimal Power Flow problem.

### B.1 AC-Optimal Power Flow problem

The AC-Optimal Power Flow (OPF) problem determines the most cost-effective generator dispatch that satisfies demand within a power network subject to various physical and engineering power systems constraints. Typically, the OPF problem involves capturing a snapshot of the power network parameters and determine the bus voltages and generator set-points based on that fixed state.
A power network can be represented as a graph $(\mathcal{N,L})$ with the node set $\mathcal{N}$ consisting of $n$ buses, and the edge set $\mathcal{L}$ comprises $l$ lines. The set $\mathcal{L}$ is defined as a collection of directed arcs, with $\mathcal{L}^{R}$ indicating the arcs in $\mathcal{L}$ but in the opposite direction. $\mathcal{G}\subset\mathcal{N}$ represents the set of all synchronous generators in the system. The power generation and demand at a bus $i\in\mathcal{N}$ are represented by complex variables $S^{r}_{i}=p_{i}^{r}+jq_{i}^{r}$ and $S^{d}_{i}=p_{i}^{d}+jq_{i}^{d}$, respectively. The power flow across line $ij$ is denoted by $S_{ij}$, and $\theta_{i}$ symbolizes the phase angles at bus $i\in\mathcal{N}$.

Figure: Model 1 The AC Optimal Power Flow Problem (AC-OPF)

The AC power flow equations use complex numbers for current $I$, voltage $V$, admittance $Y$, and power $S$, interconnected through various constraints. Kirchhoff’s Current Law (KCL) is represented by $I^{r}_{i}-{I^{d}_{i}}=\sum_{(i,j)\in\mathcal{L}\cup\mathcal{L}^{R}}I_{ij}$, Ohm’s Law by $I_{ij}={Y}_{ij}(V_{i}-V_{j})$, and AC power flow is denoted as $S_{ij}=V_{i}I_{ij}^{\ast}$.
These principles form the AC Power Flow equations, described by (11f) and (11g), which formulation is described by Model 1. The goal is to minimize a function (11a) representing dispatch costs for each generator.
Constraints (11b)-(11c) represents voltage operational limits to bound voltage magnitudes and phase angle differences, while (11d)-(11e) set boundaries for generator output and line flow. Constraint (11h) sets the reference phase angle. Finally, constraints (11f) and (11g) enforce KCL and Ohm’s Law, respectively.
The classical, *steady-state* problem, described by Model 1, does *not* incorporate systems dynamics capturing the behavior of the synchronous generators, and as such, does *not* guarantee stable operations for a power system.
This paper extends this problem by introducing the *Stability-Constrained* AC-Optimal Power Flow Problem, which integrates the generator dynamics and related stability constraints within the AC-OPF problem (1).

### B.2 Generator dynamics

The generator dynamics are modeled using the "Classical machine model" (12), which is typically adopted to describe the dynamic behavior of synchronous generators []

$$ $\mkern-18.0mu\frac{d}{dt}\begin{bmatrix}\delta^{g}(t)\\ \omega^{g}(t)\end{bmatrix}\!=\!\begin{bmatrix}\omega_{s}(\omega^{g}(t)-\omega_ {s})\\ \frac{1}{M^{g}}\left(P_{m}^{g}-D^{g}(\omega^{g}(t)-\omega_{s})-\frac{E^{g}_{q, 0}v_{g}}{X^{g}_{d}}\sin(\delta^{g}(t)-\theta_{g})\right)\end{bmatrix}\!$ (12) $$

, where $\delta^{g}(t)$ and $\omega^{g}(t)$ represents the rotor angle and angular speed over time $t$ of generator $g\in\mathcal{G}$, $\omega_{s}$ the
synchronous angular frequency, $M^{g}$ the machine’s inertia constant, $D^{g}$ the damping coefficient, $P^{g}_{m}$ the mechanical power, $X_{d}^{g}$ the transient reactance and $E^{g}_{q,0}$ electromotive force.
The initial value of the rotor angle $\delta^{g}_{0}$, and electromotive force $E^{g}_{q,0}$ for each generator $g\in\mathcal{G}$ are derived from the
active and reactive power equations, assuming the generator dynamical system (12) being in a steady state condition at time instant $t=0$, $\frac{d}{dt}\begin{bmatrix}\delta^{g}(t)\ \omega^{g}(t)\end{bmatrix}_{t=0}^{T}
=\begin{bmatrix}0\ 0\end{bmatrix}^{T}$:

$$ $\displaystyle\frac{E^{g}_{q,0}v_{g}\sin(\delta^{g}_{0}-\theta_{g})}{X^{g}_{d}} -p_{g}^{r}=0,$ (13) $\displaystyle\frac{E^{g}_{q,0}v_{g}\cos(\delta^{g}_{0}-\theta_{g})-v_{g}^{2}}{ X^{g}_{d}}-q_{g}^{r}=0.$ (14) $$

Following the same assumptions, the initial rotor angular speed is set as

$$ $\omega^{g}_{0}=\omega_{s}.$ (15) $$

Figure: Model 2 The Stability Constrained AC-OPF Problem

#### Stability limit

To guarantee stability of a synchronous generator $g\in\mathcal{G}$, the rotor angle $\delta^{g}(t)$ is required to remain below an instability threshold $\delta^{\max}$, as defined by SIngle Machine Equivalent (SIME) model:

$$ $\displaystyle\delta^{g}(t)\leq\delta^{\max}\quad\forall t\geq 0.$ (17) $$

Unstable conditions arise when violating the inequality constraint (17), which is the principal binding constraint that necessitates re-dispatching.

### B.3 Stability-Constrained AC-Optimal Power Flow Problem

The generator dynamics (12) and its initial conditions equations (13)-(15), together with the associated stability constraints (17), are thus integrated within the steady-state AC-OPF Problem 1, giving rise to the Stability-Constrained AC-OPF problem, which is detailed in Model 2.
In this problem, parameters $\bm{\zeta}=\bm{S}^{d}$ represent customer demand, while decision variables $\bm{x}=(\bm{S}^{r},\bm{V})$ are the generator settings and bus voltages; the state variables $\bm{y}(t)=(\delta^{g}(t),\omega^{g}(t))\ \forall g\in\mathcal{G}$ represent the rotor angle and angular speed of the generators.

#### Training setting of Neural Ordinary Differential Equation Models

As the generator dynamics are described by a system of ODEs, neural-ODE [] models, one for each synchronous generator $g\in\mathcal{G}$, are used to capture their dynamics.
Each neural-DE model $\mathcal{N}_{\theta}^{g}$ is trained in a supervised fashion, as described in Section 4.2, to obtain dynamic predictors that are capable of providing accurate estimate of the state variables $\bm{y}^{g}(t)$ across a family of instances of the generator model (12).
Specifically, for each generator $g\in\mathcal{G}$, the datasets $\mathcal{D}^{g}$ used for training the generator dynamic predictor $\mathcal{N}_{\theta}^{g}$, consists of pairs $(\bm{x}_{i},\bm{y}_{i}(t))\sim\mathcal{D}^{g}$, where $\bm{x}_{i}=(\delta_{0}^{g},\omega_{0}^{g},|V^{\prime}_{g}|,\theta^{\prime}_{g})$ is the input of the neural-ODE model, and $\bm{y}_{i}(t)=({\delta}^{g}(t),{\omega}^{g}(t),|V^{\prime}_{g}|(t),\theta^{
\prime}_{g}(t))$ the corresponding solution of (12) with initial conditions $\bm{y}_{i}(0)=\bm{I}(\bm{x}_{i})$, represented by (16e)-(16g) and computed using $Dopri5$, a numerical algorithm implementing an adaptive Runge-Kutta method. For each input $\bm{x}$, the OPF decision variables $|V^{\prime}_{g}|,\theta^{\prime}_{g}$ are sampled from a uniform distribution $\mathcal{U}(a,b)$, where $a$ and $b$ are given by the corresponding operational limits specified by Constraints (11b)-(11c).
Note that each of these variables influences the initial condition of the state variables $({\delta}^{g}(0),{\omega}^{g}(0))$ via Equation (13)-(15), as well as the governing equations of the generator.
As $|V^{\prime}_{g}|,\theta^{\prime}_{g}$ extend the actual generator state variables $\left(\delta^{g}(t),\omega^{g}(t)\right)$, we are implicitly augmenting the generator model with these two additional state variables that have no dynamics (e.g. $\frac{d|V^{\prime}_{g}|(t)}{dt}=0,\frac{d\theta_{g}(t)}{dt}=0$) and initial condition $|V^{\prime}_{g}|(0)=|V^{\prime}_{g}|,\theta^{\prime}_{g}(0)=\theta^{\prime}_{g}$. This trick allows us to explicitly inform the neural ODE model of the role played by the voltage magnitude $|V_{g}|$ and angle $\theta_{g}$ on the dynamics of each generator.
The generator characteristics parameters of model (12), such as the damping coefficient $D^{g}$, inertia constant $M^{g}$, and mechanical power $P^{g}_{m}$ are adopted from [].
Each dataset $\mathcal{D}^{g}$ contains approximately $50\%$ of unstable trajectories and $50\%$ of stable trajectories, generated as described in Section 4.2.
At training time, given a pair $(\bm{x}_{i},\bm{y}_{i}(t))\sim\mathcal{D}^{g}$, the target is constructed as $\bm{y}^{g}_{i}(t)=\bm{y}^{g}_{i}(0),\bm{y}^{g}_{i}(\Delta_{t}),\ldots,\bm{y}^{
g}_{i}(n\Delta_{t})$, with $\Delta_{t}=0.001$ and the number of points $n$, is set to $200$ at the beginning of the training, and gradually increases up to $1000$. This trick allows to avoid local minima during training []. At test time, we set $n=1000$.

### B.4 Dynamic Portfolio Optimization

DE-OP uses a neural-SDE [] model $\mathcal{N}_{\theta}$ to capture the asset price dynamics $\bm{y}(t)$. The neural-SDE model consists of $2$ separate neural network, $\mathcal{N}_{\theta}=(\mathcal{N}^{f}_{\theta},\mathcal{N}^{g}_{\theta})$ where $\mathcal{N}^{f}_{\theta}$ aims to capture the deterministic component of (9d), $\mu_{i}\bm{\zeta}(t)dt$ and $\mathcal{N}^{g}_{\theta}$ the stochastic component $\sigma_{i}\bm{\zeta}(t)dW_{i}(t)$.
$\mathcal{N}^{f}_{\theta}$ is a simple linear layer, while $\mathcal{N}^{g}_{\theta}$ is a $2$-layer ReLU neural network.
Given the initial asset price vector $\bm{y}(0)=\bm{\zeta}$, the neural-SDE model $\mathcal{N}_{\theta}$ generates an estimate $\hat{\bm{y}}(t)=\mathcal{N}_{\theta}(\bm{y}(0),t)$ of the asset prices trend, from which the final asset price $\hat{\bm{y}}(T)$ is obtained.
The LSTM and the Feed Forward model used to estimate the final asset price $\hat{\bm{y}}(T)$ as the dynamic component of the corresponding baseline method are both a $2$-layer ReLU neural network.
The final time instant $T=28,800$ seconds which corresponds to $8$ hours.
Given initial condition $\bm{y}^{j}(0)=\bm{\zeta}^{j}$, the asset price trend $\bm{y}^{j}(t)$ is obtained by Ito numerical integration of (9d).
The neural-SDE model is trained on a dataset $\{(\bm{\zeta}^{j},\bm{y}^{j}(t))\}_{j=1}^{N}$; the LSTM model is trained on a dataset $\{(\bm{y}^{j}(t),\bm{y}^{j}(T))\}_{j=1}^{N}$, where $\bm{y}^{j}(t)=\bm{y}^{j}(0),\bm{y}^{j}(\Delta_{t}),\ldots,\bm{y}^{j}(\Delta_{t
}K)$ is a time series, $\Delta_{t}=100$ seconds and $\Delta_{t}K=T-1$. The Feed Forward network is trained on a dataset $\{(\bm{y}^{j}(0),\bm{y}^{j}(T))\}_{j=1}^{N}$.

## Appendix C Additional Experimental Details

### C.1 Proxy optimizer methods

This subsection describes in brief the proxy optimizer methods adopted in the experiments to estimate the optimal decision variables $\bm{u}^{\star}$, within the operational setting described by (1) or (10).
Each description below assumes a DNN model $\mathcal{F}_{\omega}$ parameterized by $\omega$, which acts on problem parameters $\bm{\zeta}$ to produce an estimate of the decision variables $\hat{\bm{u}}\coloneqq\mathcal{F}_{\omega}(\bm{\zeta})$, so that $\hat{\bm{u}}\approx\bm{u}^{\star}(\bm{\zeta})$.
To ease notation, the dependency from problem parameters $\bm{\zeta}$ is omitted.

#### Lagrangian Dual Learning (LD).

Fioretto et. al.
[] constructs the following modified Lagrangian as a loss function:

$$ $\mathcal{L}^{{LD}}(\hat{\bm{u}},\bm{y}(t))=\|\hat{\bm{u}}-\bm{u}^{\star}(\bm{ \zeta})\|^{2}+\bm{\lambda}^{T}\left[\bm{g}(\hat{\bm{u}},\bm{y}(t))\right]_{+}+ \bm{\mu}^{T}\bm{h}(\hat{\bm{u}},\bm{y}(t)).$ $$

At each iteration of LD training, the model $\mathcal{F}_{\omega}$ is trained to minimize a balance of constraint violations and proximity to the precomputed target optima $\bm{u}^{\star}(\bm{\zeta})$. Updates to the multiplier vectors $\bm{\lambda}$ and $\bm{\mu}$ are calculated based on the average constraint violations incurred by the predictions $\hat{\bm{u}}$, mimicking a dual ascent method [].

#### Deep Constraint Completion and Correction (DC3).

Donti et. al. [] uses the loss function

$$ $\mathcal{L}^{\text{DC3}}(\hat{\bm{u}},\hat{\bm{y}}(t))=\mathcal{J}(\hat{\bm{u} },\hat{\bm{y}}(t))+\lambda\|\left[\bm{g}(\hat{\bm{u}},\hat{\bm{y}}(t))\right]_ {+}\|^{2}\\ +\mu\|\bm{h}(\hat{\bm{u}},\hat{\bm{y}}(t))\|^{2}$ $$

which relies on a completion-correction technique to enforce constraint satisfaction, while maximizing the empirical objective $\mathcal{J}$ in self-supervised training.

#### Self-Supervised Primal-Dual Learning (PDL).

Part et. al. [] uses an augmented Lagrangian loss function

$$ $\begin{split}\mathcal{L}^{\text{PDL}}(\hat{\bm{u}},\hat{\bm{y}}(t))=\mathcal{J }(\hat{\bm{u}},\bm{y}(t))+\hat{\bm{\lambda}}^{T}\bm{g}(\hat{\bm{u}},\hat{\bm{y }}(t))+\hat{\bm{\mu}}^{T}\bm{h}(\hat{\bm{u}},\hat{\bm{y}}(t))+\\ \frac{\rho}{2}\left(\sum_{j}[g_{j}(\hat{\bm{u}},\hat{\bm{y}}(t))]_{+}+\sum_{i} |h_{i}(\hat{\bm{u}},\hat{\bm{y}}(t))|\right),\end{split}$ $$

which consists of a primal network to approximate the decision variables, and a dual network to learn the Lagrangian multipliers update. The method is self-supervised, requiring no precomputation of target solutions for training.

#### MSE.

Zamzam et. al. [] uses the loss function:

$$ $\mathcal{L}^{\text{MSE}}(\hat{\bm{u}},\bm{u}^{\star})=\|\hat{\bm{u}}-\bm{u}^{ \star}\|^{2}$ $$

which minimizes the MSE between the predicted solution $\hat{\bm{u}}$ and the corresponding precomputed solution $\bm{u}^{\star}$.

### C.2 Stability Constrained AC-Optimal Power Flow experiment

#### Hyperparameters of the neural-ODE models

Each neural-ODE model is a fully connected feedforward ReLU neural network with 2 hidden layers, each with 200 units.
Each model is trained using Adam optimizer, with default learning rate $\eta=10^{-3}$ and default hyperparameters.

#### Hyperparameters of the DE-OP’s optimization model and the proxy optimizer methods

The DE-OP optimization component $\mathcal{F}_{\omega}$ and each baseline proxy optimizer model is trained with Adam optimizer and default learning rate $\eta=10^{-3}$.
Each proxy optimizer model is a fully connected FeedFoward ReLU neural network with $5$ hidden layers, each with $200$ units.
The DE-OP’s optimization model $\mathcal{F}_{\omega}$ and Lagrangian Dual proxy model are trained with a Lagrangian step size $\rho=10^{-1}$, while the Lagrangian multipliers $\bm{\lambda}_{h^{\prime}}$ and $\bm{\lambda}_{g}$ are updated at each epoch.
DC3’s proxy optimizer model is trained with the same set of hyperparameters for OPF, as reported in the original paper.

### C.3 Dynamic Portfolio Experiment

#### Hyperparameters of the asset prices predictor models

The stochastic component of the neural-SDE, the LSTM and Feed Forward model are each 2-layers ReLU networks, each with $100$ units.
The neural-SDE, LSTM and Feed Forward models are all trained using Adam optimizer, with default learning rate $\eta=10^{-3}$ and hyperparameters.

#### Hyperparameters of the DE-OP’s optimization model and the proxy optimizer methods

The DE-OP optimization component $\mathcal{F}_{\omega}$ and each baseline proxy optimizer model is trained with Adam optimizer and default learning rate $\eta=10^{-3}$.
Each proxy optimizer model is a fully connected Feed Foward ReLU neural network with $2$ hidden layers, each with $50$ units.
The DE-OP’s optimization model $\mathcal{F}_{\omega}$ and Lagrangian Dual proxy model are trained with a Lagrangian step size $\rho=10^{-1}$, while the Lagrangian multipliers $\bm{\lambda}_{h^{\prime}}$ and $\bm{\lambda}_{g}$ are updated each $10$ epoch.
PDL’s and DC3’s proxy optimizer model uses the same hyperparameters for the Convex Quadratic task, as reported in the respective original paper.

## Appendix D Additional experimental results: Stability Constrained AC Optimal Power Flow

This section reports additional experimental results of DE-OP model and the baseline methods on the WSCC 9-bus system and the IEEE 57-bus system.
Specifically, we report:

- •
The inference time (measured in seconds) of neural-ODE which we compare to the computational time of a traditional numerical ODE solver and the precision of state variables’ estimate $\hat{\bm{y}}(t)$ (measured as the percentage $\ell^{2}$ error) of NODEs and PINNs [], a different learning-based approach for learning the system dynamics.
- •
The (steady-state) decision error (MSE) of the OPF variables of DE-OP and each proxy optimizer method, incurred by the respective approximation $\hat{\bm{u}}$, computed assuming that the generators are in steady-state.
- •
The (steady-state) optimality gaps incurred by DE-OP and the baselines proxy optimizers predictions’ $\hat{\bm{u}}$, and measured as $\frac{|L(\bm{u}^{\star}(\zeta),\bm{y}^{\star}(T))-L(\hat{\bm{u}}(\zeta),\hat{
\bm{y}}(T))|}{|L(\bm{u}^{\star}(\zeta),\bm{y}^{\star}(T))|}\times 100$,
where $L$ is objective function (equation 16a).
- •
The inference time (measured in seconds) of DE-OP and each proxy optimizer model to generate $\hat{\bm{u}}$.

**Table 2: Average and standard deviation of computational times for numerical solvers vs. neural-ODE inference time by method**
| Method | Numerical Solver | neural-ODE |
| --- | --- | --- |
| Dopri5 (default) | $0.135\pm 0.015$ (sec) | $0.008\pm 0$ (sec) |
| Bosh3 | $0.446\pm 0.039$ (sec) | $0.017\pm 0$ (sec) |

Figure: Figure 6: True and neural-ODE (NODE) solutions of the generator state variables in unstable conditions.
Refer to caption: x6.png

Figure: Figure 7: True, neural-ODE (NODE) and PINN estimate of the generator state variables in stable conditions.
Refer to caption: x7.png

### D.1 Learning the generator dynamics

#### Runtime comparison between neural-ODEs and a traditional ODE solver.

Here the goal is to evaluate the neural-ODE’ inference time to produce the generators’ state variables estimates and to compare it with the computational time of a traditional ODE solver.
Given the synchronous generator model described by
(12), a numerical ODE solver could be adopted to
determine the evolution in time of the state variables $\delta^{g}(t)$ and $\omega^{g}(t)$. However, in case of unstable conditions, the system response can be as rapid as, or even exceed, the time required for computing the ODE solution with a numerical solver.
This situation is depicted in Figure 6 where unstable conditions arise before a numerical solution to the system of differential equations (12) is computed. Conversely, the neural ODE model $\mathcal{N}^{g}_{\phi}$ is capable of detecting unstable conditions before the system transitions into an unstable state, while also providing a good approximation of the solution. This speed advantage arises from the neural-ODE’ vector field approximation of (12), which enables quicker computation of the forward pass of a numerical ODE solver []. Table 2
reports the average and standard deviation of computational time, for numerical solvers, and inference time, for neural-ODEs, given $2$ different numerical algorithms. *Neural-ODE models are, on average, about $20$ times faster than a numerical solver which uses the dynamic equations of (12). This aspect makes neural-ODE models natural candidates as dynamic predictors for the generator model in real-time applications.*

#### Comparison between neural-ODEs and PINNs.

Here the goal is to asses the precision of the neural-ODEs’ estimate of the generator state variables and to compare them with PINNs []. PINNs are ML-based models that incorporates known physical laws into the learning process. Instead of relying solely on data, PINNs use physics-based constraints to guide the training, ensuring that the model’s predictions are consistent with the underlying scientific principles.
Figure 7 shows the neural-ODE and PINNs’ state variables estimates in case of stable conditions.
While a neural-ODEs model produces highly accurate state variables’ predictions, a PINN model trained on the same dataset $\mathcal{D}^{g}$ but affected by a generalization bias, is incapable of capturing the generator dynamics across different instances of the generator model (12) and produces poor state variables estimates. Specifically, the percentage $\ell^{2}$ error between the numerical ODE solver solutions $\delta(t),\omega(t)$ and the neural-ODE (NODE) predictions $\delta_{\text{NODE}}(t),\omega_{\text{NODE}}(t)$ is $5.17\%$, while
for the PINN predictions $\delta_{\text{PINN}}(t),\omega_{\text{PINN}}(t)$ is significantly higher at $69.45\%$.

Table: Table 3: Average (steady-state) decision errors for the WSCC $9$-bus system across different approaches based on $40$ trials.

Table: Table 4: Average (steady-state) decision errors for the IEEE $57$-bus system across different approaches based on $40$ trials.

Table: Table 5: Average and standard deviation of constraint violations and (steady-state) optimality gap on the WSCC $9$-bus system for different approaches based on $40$ independent runs.

### D.2 Constraint violations and (steady-state) decision errors

Tables 3 and 4 report the (steady-state) decision error (MSE) at test time of DE-OP and the baseline PO methods on the WSCC 9-bus and IEEE 57-bus system, respectively. Specifically, the tables reports the MSE of estimated solutions with respect to the ground-truth solutions, with the assumption that the synchronous generators are in steady-state. The same considerations reported in Section 5.2 regarding the steady-state optimality gaps apply also for the following discussion. In other words, the ground truth variables used to compute the decision errosr, are obtained from solving the *steady-state* ACOPF problem, and as such, the decision errors here reported do not necessarily reflect that of each method for solving the Stability-Constrained AC OPF problem. Nonetheless, this metric provides valuable insights on the impact of the decision variables on the dynamics requirements of Problem 2.

Firstly, note that, for each test case, all the methods achieve similar decision error.
Despite that, as shown in Table 1 of the main paper and Table 5, and Figures 5, 5 of the main paper, DE-OP is the only method that satisfy exactly the dynamic requirements (17), while all the baseline methods systematically violate the stability constraint. These results suggest that DE-OP modifies potentially unstable set points, at a cost of a only slightly higher MSE than the baseline approaches. Note in particular the MSE error of the OPF variables $|V|$ and $\theta$; these variables directly affect the generator dynamics in (12), and thus their modification is necessary to satisfy the stability constraint.
This trade-off is crucial for practical applications, where the dynamic requirements must be addressed.
Table 5 shows the violation of the static (flow and boundaries), along with the optimality gap with the assumption that the generators are in steady-state, for each method on the WSCC-9 bus system. Similarly to the IEEE 57 test-case discussed in Section 5.2, DC3 is the only method which achieves steady-state constraint satisfactions. All methods except DC3 generate comparable violations of the flow balance constraints, which is the most difficult constraint to satisfy due to its non-linear nature defined by Constraint (11f) and (11g). LD satisfy the boundary constraint by projecting its output $\hat{\bm{u}}$ within the feasible set defined by Constraints (11b)-(11c). Empirically, we found that removing this projection operation within the DE-OP model $\mathcal{F}_{\omega}$, in some cases allows to satisfy the dynamic requirements. We did not thoroughly investigate this result, but our intuition is that in some cases the decision variables $|V|$, $\theta$ involved in the stability analysis must assume values close to their boundaries to satisfy stability constraints. This comes at a cost of minimal boundary constraint violations from DE-OP. MSE, lacking of a mechanism to encourage constraint satisfactions, produces solutions violating each constraint function.

### D.3 Steady-state optimality gaps

This subsection discusses the sub-optimality of the estimated solution $\bm{\hat{u}}$ with respect to the ground truth $\bm{u}^{\star}$, with the assumption that the generators are in steady-state conditions, given parameters $\bm{\zeta}$ and measured in terms of objective value (16a), of DE-OP and each baseline method. The same considerations reported in Section 5.2 regarding the steady-state optimality gaps and how this metric should be interpreted, apply also for the subsequent discussion.
Table 1 in the main paper and Table 5 report the optimality gaps on the WSCC-9 and IEEE-57 bus system, respectively.
The tables report that all the methods achieve comparable gaps on each test case. This is intuitive, since the optimality gap depends solely on the power generated $\bm{p}^{r}$ (see objective function (16a)), and all methods produce similar $\bm{p}^{r}$’s prediction error, as displayed in Tables 3 and 4.
For the WSCC bus system, DE-OP produces average optimality gap of $0.13\%$ while preserving system stability, that are comparable with the best optimality gap - LD with $0.11$ and DC3 with $0.11$ - which often violates stability constraints.

**Table 6: Average and standard deviation of inference times for different OPF learning approaches in the test cases.**
| Models | WSCC 9-bus | IEEE 57-bus |  |
| --- | --- | --- | --- |
| $\mathcal{F}_{\omega}$ | $\mathcal{N}_{\theta}$ | Inference Time (sec) |  |
| DE-OP (ours) | $0.001\pm 0.00$ | $0.009\pm 0.00$ |  |
| MSE | $\emptyset$ | $0.000\pm 0.00$ | $0.001\pm 0.00$ |
| LD | $\emptyset$ | $0.000\pm 0.00$ | $0.001\pm 0.00$ |
| DC3 | $\emptyset$ | $0.025\pm 0.00$ | $0.089\pm 0.00$ |

### D.4 Inference time

Finally, we evaluate the average inference time of DE-OP and each baseline proxy optimizer method. Table 6 shows the inference time of each proxy optimizer method on each test case.
On average, DE-OP produces near-optimal and stable solutions in $1$ (ms) and $9$ (ms) for the WSCC-9 bus and IEEE 57-bus, respectively, which is slightly higher but comparable with the MSE and LD approaches, and about $15\times$ faster than the DC3 method.
DC3 achieves the highest inference time, due to its correction and completion procedure, which requires solving a nonlinear system of equations and the Jacobian matrix computation. While DE-OP can be already used for near real-time applications, its efficiency could be improved by computing the state variables in parallel, since each dynamic predictor is independent. *This aspect makes DE-OP’ inference time independent from the size of number of state variables and dynamical systems, suggesting potential for large-scale and complex systems*.

## Appendix E Additional experimental results: Dynamic Portfolio Optimization

Figure: Figure 8: Asset prices (blue), LSTM (red) and neural-SDE asset prices estimates.
Refer to caption: x8.png

### E.1 Learning the asset price dynamics

#### Comparison between neural-SDE and LSTM.

Fig. 8 illustrates the asset price trends given different initial asset prices $\zeta^{0}$, with estimates $\hat{\bm{y}}^{0}(t)$ provided by both a neural-SDE model and an LSTM model, alongside the true asset prices $\bm{y}^{0}(t)$ computed with a numerical SDE solver implementing Euler-Maruyama method [].
It is evident that, given different initial asset price $\zeta^{0}$, the neural-SDE model produces more accurate predictions than the LSTM model, by explicitly capturing the asset pricing dynamic equations. These accurate predictions lead to more informed and higher quality decision making, as discussed in Section 5.1.

Figure: Figure 9: Average percentage optimality gap with $n=20$ asset prices for each method across different proxy optimizer methods.
Refer to caption: x9.png

### E.2 Optimality gaps

This section report additional results of DE-OP and the baseline methods across different proxy optimizer methods and asset price predictors on the Dynamic Portfolio Optimization task with $n=20$ and $n=50$.
Figures 10 and 10 display the average, percentage *optimality gap* on the test set, across different methods, for $n=20$ and $n=50$, respectively.
In both figures, it is evident that for a given proxy optimizer method (e.g., Lagrangian Dual), by using neural-SDE predictors to capture the asset prices dynamics, DE-OP yields superior decision quality compared to the baseline methods.
For $n=20$, DE-OP achieves the lowest optimality gaps - $12.92\%$ for DC3, $5.23\%$ for LD, and $13.45\%$ for PDL - by capturing the asset prices dynamics via explicit modeling of the asset prices’ dynamics. Predicting the final asset price with LSTM leads to optimality gaps of $8.67\%$ for LD, and $19.00\%$ for PDL for PDL, performing consistently worse than DE-OP, due to its lack of explicit dynamic modeling. The Feed Forward model performs significantly worse, leading to significantly higher gaps - $121.56\%$ for DC3, $100.45\%$ for LD, and $110.33\%$ for PDL - highlighting its limitations in capturing the time-dependent nature of the data, similarly to the proxy optimizer methods which ignore the system dynamics, which achieve $103.98\%$ for DC3, $111.63\%$ for LD, and $115.11\%$ for PDL. Overall, these results follow the trend reported in Figure 3 and discussed in Section 5.1 of the main paper, concerning the optimization task with $n=50$. Among the proxy optimizer methods considered, Lagrangian Dual consistently outperforms DC3 and PDL, suggesting that precomputed solutions can enhance the accuracy and robustness of optimization proxies. The optimality gaps achieved by each method when $n=50$, increase with respect to the optimality gaps achieved by the corresponding method when $n=20$, likely due to a higher complexity of the optimization task.
These results highlight the importance of accurate dynamic predictions, which in turn enable, in a subsequent stage, generating high quality investment allocations.

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 8
presentation: 4
questions: - In the Related Works section, the authors do not mention collocation-based or differential algebraic programming (DAE), although these are mentioned later in the results/discussion section. How do these methods compare to proposed approach, why are they relevant?
- How does "multiparametric programming" (which is a field related to MPC, developed to identify "maps" of control actions under different parametrizations, for which authors have also previously used surrogate models) compare to proposed approach? 
- Finally, the authors mention that DE-optimization problems involve parametrization of the state variables, and thus traditional MPC is unsuitable. This statement must be made more clear because it is a little misleading. What is traditional MPC and why is it unsuitable? Due to implementation or computational cost concerns?
- How would the method steps (training for Fw, training for Nθ), using Langrangian approach scale for larger grid systems? Please discuss this in your conclusions. 
- How critical is it to the method to have near-optimal guesses to solve the optimization problem with IPOPT? Would this ever become an issue for larger systems? 
- For both surrogate models trained, what is the trade-off between training sample diversity (covering the entire space) vs. optimality (samples being close to optimal regions) for accuracy and constraint violation of the overall method?
soundness: 4
strengths: - Well-written paper, well organized and well-motivated.
- The selected case studies are important, challenging and sufficient to show the effectiveness of the proposed approach. 
- Authors have done a lot of interesting and relevant work, presenting many more comparisons and discussion in the Appendices.
summary: In this paper the authors propose a dual learning approach for solving differential equation-based optimization problems, relevant to dynamic optimization and systems control. The authors present a financial and a power-grid case study, showcasing the accuracy and speed of their proposed method. The method includes several novel aspects and its applicability is well-motivated by the authors.
weaknesses: Related works section: 
- This work is relevant to control applications, which is a really rich area with respect to literature. The Reviewer realizes that it would be impossible to review all of it, but some statements and key relevant areas could be better represented. Several questions on this have been added to the Questions section. 
- The motivation includes stochastic DE constrained optimization, which is not shown in results, rather briefly mentioned in Appendix A. The discussion provided in Appendix A is not sufficient to justify that this could extend to stochastic opt, which comes with more challenges. 
- The authors only solve up to a 57-bus system, which is a benchmark problem in power systems optimization. Real problems are much larger than this. The authors should discuss scalability of proposed method and challenges in each step as number or nodes, variables, constraints increases. 
- The authors use a local optimizer to generate data, ensuring near-optimal solutions are used to warm-start the solver. The sensitivity of the method to this is not discussed.

---

confidence: 3
contribution: 3
initial_rating: 8
presentation: 3
questions: 1. Does the dual-network architecture need to consider synchronization issues during training? When one network changes, can the other network quickly adjust to adapt to the new state?
2. The authors have not sufficiently addressed the practical challenges of applying their method. Key aspects such as computational resource requirements, scalability to larger systems, and performance in the presence of noisy or incomplete data require comprehensive discussion.
3. The paper does not specify the numerical solvers used for comparison. In equations (3a) and (3b), where N_\theta approximates the evolution of state variables, how does the performance of the neural network solver compare to traditional numerical solvers like Runge-Kutta or Euler methods? A comparison focusing on accuracy, efficiency, and convergence is needed for completeness.
4. The loss function integrates the optimization of state and control variables, but how are the contributions of these two networks balanced?
5. The proposed network should discuss its numerical stability during the solving process, especially in cases where F is a highly nonlinear or stiff equation. How does the network avoid issues of gradient explosion or vanishing?
6. The discussion on Neural ODE in the abstract should include a reference to the paper by Chen et al. (2018) titled "Neural Ordinary Differential Equations," presented at NeurIPS, to provide necessary background and support for the proposed method.
soundness: 2
strengths: 1.	The writing is clear and well-structured, making complex concepts accessible.
2.	The methodology is rigorous and well-executed, with robust experiments that provide reliable results.
3.	The authors provide thorough comparisons with existing optimization methods, demonstrating significant performance improvements.
summary: The paper proposes a novel method for solving optimization problems constrained by differential equations (DEs). This method combines neural networks with a proxy optimization model, employing a dual-network architecture: one network approximates control strategies, while the other solves the differential equations related to these decisions. Experiments have validated the effectiveness of this method in fields such as financial modeling and energy optimization. Compared to traditional proxy optimization methods, this approach shows higher accuracy when handling high-dimensional dynamic constraints.
weaknesses: 1.	The authors do not adequately address potential limitations or challenges faced when applying their method in practice.
2.	The absence of diagrams or flowcharts makes complex concepts challenging to follow. Adding a diagram illustrating the dual-network architecture, particularly the interaction between the proxy optimization model F_\omega  and the solver N_\theta significantly improves understanding.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: * Is the upper bound on the number of points a typo? It seems should be T rather than \top
* It is hard to follow the paragraph "Optimizing over the distribution of instances". Why are the parameters \zeta used as input of the neural network? 
* How is the dataset D acquired for the experiment?
* Why is DE-OP compared to LSTM models in the experiments? Recurrent neural networks are generally unsuitable for DE-constrained optimization, so the rationale behind this comparison is unclear. Provide a clear justification for including LSTM models in the comparison. If there's a specific reason for this choice, please explain it.
soundness: 2
strengths: * The paper is well-structured and presented in a clear academic style. 
* The experiment results demonstrate that the DE-OP can be 25 times faster than existing methods, which is a great improvement. 
* Benchmarks are chosen as some real-world scenarios (AC-Optimal power flow optimization and dynamic portfolio optimization), which are beneficial for actual applications.
* The constrained optimization with differential equations is quite important for real-world control problems.
summary: This paper introduces the Differential Equation Optimization Proxy (DE-OP) for solving differential equation (DE)-constrained optimization problems. The major contribution is a proposed primal-dual framework to ensure constraint satisfaction while allowing end-to-end differentiation. Experiments on several real-world background optimization problems show that the DE-OP achieves improvements in solution stability compared with existing static methods.
weaknesses: * Contribution 2 appears to be more of a methodological explanation of DE-OP rather than a clear, original contribution. Clarifying its novelty could help highlight its value within the paper. 
* The description of the challenges addressed by DE-OP could be more specific. Currently, DE-OP mentions a wide range of problems such as NP-hard, high dimensionality, and nonlinearity that are common to many “learning-to-optimize” frameworks. A detailed description of the unique challenges faced by DE-OP would make the problem more specific.
* The proposed DE-OP is similar to many existing optimal control networks, such as the actor-critic model. Provide more details and comparisons between them or at least discuss the differences between DE-OP with neural networks for solving differential equations. 
* For equation (6), when L^{DE-OP} is minimized, how does DE-OP ensure constraint satisfaction? This is not explained in the main text, and Appendix C.1 only deals with the proxy optimizer approach, leaving this key detail ambiguous.
* In the experimental setting, baselines focus on “steady-state” aspects, raising questions about whether these baselines can effectively address the problem as formulated in Equation (1). This focus also raises potential fairness concerns in the comparison. Authors could either justify their choice of baselines more thoroughly, or include additional baselines that are more directly comparable to DE-OP in terms of handling dynamic aspects of the problem.
* The results in Table 1 do not strongly support DE-OP's superiority over baselines, as only stability violations are minimized, while other metrics do not place DE-OP even as the second-best method.  Please provide more context or discussion around these results. For example, you could explain why minimizing stability violations is particularly important in their problem setting or discuss the trade-offs between different performance metrics.

Meta-review: This paper proposes a dual-learning approach to solving optimization problems with differential equation constraints. The approach uses one network to map a problem-specific vector to the control solution and another network to approximate the differential equations. The training algorithm leverages steady-state solutions as a warm-up and employs a primal-dual optimization framework to enforce constraints effectively. The method’s performance is demonstrated through two case studies: one in finance and the other in power-grid optimization, highlighting its accuracy and computational efficiency.

All three reviewers provided positive feedback on this work. The paper is well-structured, written in a clear academic style, and introduces several novel aspects. The authors also provide strong motivation for the method’s applicability. Based on these strengths, the AC recommends acceptance.

Decision: Accept (Spotlight)