Parallel Solution of Nonlinear Projection Equations in a Multi-Task Learning Framework Dawen Wu and Abdel Lisser

Abstract-Nonlinear projection equations (NPEs) provide a unified framework for addressing various constrained nonlinear optimization and engineering problems. However, when it comes to solving multiple NPEs, traditional numerical integration methods are not efficient enough. This is because traditional methods solve each NPE iteratively and independently. In this paper, we propose a novel approach based on Multi-Task Learning (MTL) for solving multiple NPEs. The solution procedure is outlined as follows. First, we model each NPE as a system of ordinary differential equations (ODEs) using neurodynamic optimization. Second, for each ODE system, we use a Physics-Informed Neural Network (PINN) as the solution. Third, we use a multi-branch MTL framework, where each branch corresponds to a PINN model. This allows us to solve multiple NPEs in parallel by training a single neural network model. Experimental results show that our approach has superior computational performance, especially when the number of NPEs to be solved is large.

Index Terms-Nonlinear projection equations, Multi-task learning, Ordinary differential equations, Physics-informed neural networks, Error analysis

I. INTRODUCTION

A nonlinear projection equation (NPE) consists of a nonlinear mapping function and a feasible set, and the objective is to find the fixed point of the NPE problem. As reported in the literature [START_REF] Harker | Finite-dimensional variational inequality and nonlinear complementarity problems: a survey of theory, algorithms and applications[END_REF], [START_REF] Robinson | Normal maps induced by linear transformations[END_REF], NPEs provide a unified framework for modeling a variety of nonlinear optimization problems, including complementarity problems, variational inequalities, and equilibrium point problems. Such problems have numerous real-world applications in fields such as economics, engineering, and computer science.

Typically, NPEs are addressed through neurodynamic optimization, which models the problem as a system of ordinary differential equations (ODEs) [START_REF] Xia | A new neural network for solving nonlinear projection equations[END_REF]- [START_REF] Che | A collaborative neurodynamic approach to global and combinatorial optimization[END_REF]. The constructed ODE system must exhibit the global convergence property, ensuring that the state solution of the system converges to the NPE solution as the time variable approaches infinity, irrespective of the initial point. Consequently, the NPE problem is transformed into solving the state solution of the ODE system. However, this ODE system is often highly nonlinear and lacks analytical solutions. As a result, numerical integration methods such as Runge-Kutta (RK) methods or backward differential formulation are frequently employed to solve for the state solution [START_REF] Burden | Numerical analysis[END_REF].

Motivation. Traditional methods for solving NPEs typically rely on neurodynamic approaches using numerical integration. While feasible, these methods are computationally inefficient for handling large or multiple instances of NPE problems. This inefficiency is particularly pronounced in applications such as energy markets and traffic management, where it is often necessary to solve multiple nonlinear optimization problems. It is worth noting that many of these nonlinear optimization problems can be reformulated as NPE problems [START_REF] Harker | Finite-dimensional variational inequality and nonlinear complementarity problems: a survey of theory, algorithms and applications[END_REF]. For example, pricing and demand management in electricity markets can be formulated as many NPE problems, and solving them iteratively is time-consuming [START_REF] Singh | Variational inequality formulation for the games with random payoffs[END_REF]. To overcome these challenges, we propose a Multi-Task Learning (MTL) framework. This framework solves these interrelated NPE problems in parallel, significantly reducing computational cost.

A. Related works

Neurodynamic optimization. Over the past few decades, a variety of neurodynamic models have been developed to address diverse constrained optimization problems, encompassing linear and quadratic programming [START_REF] Tank | Simple 'neural' optimization networks: An a/d converter, signal decision circuit, and a linear programming circuit[END_REF], [START_REF] Liu | A one-layer recurrent neural network with a discontinuous hard-limiting activation function for quadratic programming[END_REF], general convex programming [START_REF] Kennedy | Neural Networks for Nonlinear Programming[END_REF]- [START_REF] Guo | A one-layer recurrent neural network for pseudoconvex optimization subject to linear equality constraints[END_REF], biconvex optimization [START_REF] Che | A two-timescale duplex neurodynamic approach to biconvex optimization[END_REF], pseudoconvex optimization problems [START_REF] Liu | A one-layer recurrent neural network for nonsmooth pseudoconvex optimization with quasiconvex inequality and affine equality constraints[END_REF], time-varying optimization problems [START_REF] Zhang | An exponential-type antinoise varying-gain network for solving disturbed time-varying inversion systems[END_REF]- [START_REF] Zhang | A review on varying-parameter convergence differential neural network[END_REF]. Specialized applications have also been explored, such as supervised feature selection through fractional programming [START_REF] Wang | A neurodynamic optimization approach to supervised feature selection via fractional programming[END_REF], task assignment in multivehicle systems [START_REF] Wang | Multivehicle task assignment based on collaborative neurodynamic optimization with discrete hopfield networks[END_REF], and mitigated TDOA localization [START_REF] Xiong | Tdoa-based localization with nlos mitigation via robust model transformation and neurodynamic optimization[END_REF]. In particular, a projection neurodynamic model for solving NPEs was introduced and demonstrated global convergence to the exact solution under mild conditions [START_REF] Xia | A new neural network for solving nonlinear projection equations[END_REF]. This model also exhibited both asymptotic and exponential stability without requiring a smooth nonlinear mapping. To enhance performance, a biprojection neurodynamic model was devised to efficiently solve quadratic optimization problems [START_REF] Xia | A Bi-Projection Neural Network for Solving Constrained Quadratic Optimization Problems[END_REF]. Moreover, a collaborative approach combining the projection neurodynamic model with particle swarm optimization was introduced for global optimization problems [START_REF] Che | A collaborative neurodynamic approach to global and combinatorial optimization[END_REF].

Physics-informed neural networks.

Another research direction explored in this paper concerns the application of deep learning for solving differential equations. This concept was initially introduced in the 1990s, where neural networks were trained to minimize a loss function that incorporated both boundary conditions and differential equations [START_REF] Dissanayake | Neural-networkbased approximations for solving partial differential equations[END_REF]. Subsequent research demonstrated that network architecture could be specifically designed to satisfy boundary conditions [START_REF] Lagaris | Artificial neural networks for solving ordinary and partial differential equations[END_REF], [START_REF] Mcfall | Artificial neural network method for solution of boundary value problems with exact satisfaction of arbitrary boundary conditions[END_REF]. With the rise of deep learning, this approach has regained interest for solving high-dimensional nonlinear partial differential equations [START_REF] Han | Solving high-dimensional partial differential equations using deep learning[END_REF], [START_REF] Huang | Partial Differential Equations Meet Deep Neural Networks: A Survey[END_REF]. A significant contribution to this field is the development of physics-informed neural networks (PINNs) [START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF], which integrate physical laws and data errors into the loss function. The flexible network structure and efficient training algorithms of PINNs have led to numerous successes across a wide range of computational problems in physics and engineering [START_REF] Cai | Physicsinformed neural networks (PINNs) for fluid mechanics: A review[END_REF]. This expanding body of work is propelled by a collocation strategy that adapts PINNs to exploit the structural properties of the target problem. Consequently, several variations of PINNs have been developed to address distinct problem scenarios [START_REF] Lu | Physics-Informed Neural Networks with Hard Constraints for Inverse Design[END_REF], [START_REF] Zhang | Learning in Modal Space: Solving Time-Dependent Stochastic PDEs Using Physics-Informed Neural Networks[END_REF] and enhance computational performance [START_REF] Fang | A high-efficient hybrid physics-informed neural networks based on convolutional neural network[END_REF]- [START_REF] Rezaei | A mixed formulation for physics-informed neural networks as a potential solver for engineering problems in heterogeneous domains: Comparison with finite element method[END_REF]. To facilitate the use of deep learning for solving differential equations, several packages have been developed [START_REF] Lu | DeepXDE: A Deep Learning Library for Solving Differential Equations[END_REF], [START_REF] Chen | NeuroDiffEq: A Python package for solving differential equations with neural networks[END_REF].

Multi-task learning.

Multi-task learning (MTL) is a machine learning approach that simultaneously optimizes multiple related tasks, improving generalization and performance by sharing information and leveraging commonalities among tasks [START_REF] Ruder | An overview of multi-task learning in deep neural networks[END_REF]. MTL algorithms primarily focus on training a single model capable of solving various tasks, and this approach has been applied to diverse fields such as reinforcement learning, computer vision, natural language processing, and robotics [START_REF] Bakker | Task clustering and gating for bayesian multitask learning[END_REF]- [START_REF] Mayr | Skill-based multi-objective reinforcement learning of industrial robot tasks with planning and knowledge integration[END_REF]. Addressing the challenging learning problems in MTL involves various strategies. Architectural solutions include employing multiple modules or paths, attention-based architectures, or decomposing the problem into local problems corresponding to individual tasks, which are easier to learn [START_REF] Rosenbaum | Routing networks: Adaptive selection of non-linear functions for multi-task learning[END_REF]- [START_REF] Maninis | Attentive singletasking of multiple tasks[END_REF]. These local models are then integrated into a single multi-task policy using diverse distillation techniques [START_REF] Czarnecki | Distilling policy distillation[END_REF].

B. Contributions

The main contributions of this paper are as follows:

• We present an MTL framework designed to solve multiple NPE problems simultaneously, in contrast to traditional numerical methods that treat each NPE problem individually. The proposed MTL framework consists of two parts, the shared layers and the task-specific layers.

The shared layers are responsible for processing the common features across multiple NPE problems, while the task-specific layers are designed to handle each individual NPE. • We conduct an error analysis for the proposed neural network approach. Our study shows that the total error in our method consists of three components: Neurodynamic Error, Generalization Error, and Optimization Error. We investigate the factors that influence these errors, including the time range, the network architecture, and the model training. Finally, we show that even a single-layer neural network is capable of converging to the optimal solution of the NPE. • We apply the proposed MTL framework to a variety of NPE scenarios. Experimental results show that our MTL framework exhibits significant computational performance benefits, especially when a large number of instances need to be solved or when the required solution accuracy is relatively lax. In particular, when solving for 100 NPE instances with a target error reduction down to 1, our MTL framework is 3 times faster than traditional numerical methods. Moreover, we observe that the computational time for our MTL framework grows more slowly than that of traditional numerical methods as the number of instances to solve increases. These experimental results validate the efficiency of the proposed MTL framework.

C. Outline

The remainder of this paper is organized as follows. Section II provides the necessary background, including an introduction to the NPE problem and how a neurodynamic optimization approach models it. Section III reformulates an NPE problem into a neural network training problem using PINNs, accompanied by an error analysis. Section IV introduces the proposed MTL framework and demonstrates its ability to solve multiple NPE problems in parallel. Section V presents experimental results and compares them with various numerical solvers. Finally, Section VI summarizes the key findings of this paper and outlines directions for future research.

II. NEURODYNAMIC APPROACH FOR MODELLING NPE

A. NPE Definition 1 (Nonlinear projection equation). Consider a nonlinear mapping G : R d → R d and a feasible set Ω ⊂ R d . The projection function P Ω : R d → Ω maps a vector z ∈ R d onto Ω, such that:

P Ω (z) = arg min x∈Ω ∥z -x∥, (1)
where ∥ • ∥ denotes the Euclidean norm.

The NPE problem, denoted by N P E(Ω, G), is to find a vector x * ∈ Ω satisfying:

P Ω (x * -G(x * )) = x * . ( 2
)
Definition 2 (Nonlinear complementarity problem). Consider a nonlinear mapping G : R d → R d . The nonlinear complementarity problem, denoted by N CP (G), is to find a vector x * ∈ R d satisfying:

G(x * ) ≥ 0, x * ≥ 0, G(x * ) T x * = 0. (3)
Definition 3 (Variational inequality). Consider a nonlinear mapping G : R d → R d and a feasible set Ω ⊂ R d . The variational inequality problem, denoted by V I(Ω, G), is to find a vector x * ∈ Ω satisfying: According to the literature [START_REF] Harker | Finite-dimensional variational inequality and nonlinear complementarity problems: a survey of theory, algorithms and applications[END_REF], NPEs can be viewed as a unified framework for many nonlinear optimization problems. For example, the Karush-Kuhn-Tucker (KKT) conditions for linear and quadratic programming problems can be represented as linear complementarity problems, and the KKT conditions for convex constraint nonlinear optimization problems can be transformed into nonlinear complementarity problems [START_REF] Billups | Complementarity problems[END_REF]. Both are recast as NPEs according to Proposition 1. Many Nash equilibria in game theory can be represented by variational inequalities [START_REF] Facchinei | On generalized Nash games and variational inequalities[END_REF], [START_REF] Singh | Variational inequality formulation for the games with random payoffs[END_REF], which are then reformulated as an NPE via Proposition 2.

(x -x * ) T G (x * ) ≥ 0, x ∈ Ω. (4

B. Neurodynamic Approach

Assumption 1.

• The function G(•) is locally Lipschitz continuous.

• The feasible set Ω is a box-constrained set, defined as

Ω = {x = (x 1 , x 2 , . . . , x n ) ∈ R n | l i ≤ x i ≤ h i },
where l i and h i denote the lower and upper bounds of x i respectively. In this case, the projection function in Eq. ( 1) is reduced to P Ω (x) = P 1 Ω (x 1 ), P 2 Ω (x 2 ), . . . , P d Ω (x d ) T , where P i Ω (x i ), i ∈ {1, 2, . . . , d} is defined as:

P i Ω (x i ) =      l i if x i < l i , x i if l i ≤ x i ≤ h i , h i if x i > h i . (5)
Consider a time-dependent function y : R → R n , where y(t) denotes the state at time t. The objective of neurodynamic optimization is to design a first-order ODE system to govern y(•). In this paper, we utilize the projection neurodynamic model proposed by [START_REF] Xia | A new neural network for solving nonlinear projection equations[END_REF] to model the NPE, wherein the ODE system is defined as follows:

dy dt = -G (P Ω (y)) + P Ω (y) -y. (6)
To simplify the discussion, we define:

Φ(y) = -G (P Ω (y)) + P Ω (y) -y. (7)
Hence, the ODE system (6) can be expressed as dy dt = Φ(y). Definition 4 (State solution). Given an ODE system dy dt = Φ(y), where Φ : R d → R d and an initial condition (t 0 ∈ R, y 0 ∈ R d ), a vector-valued function y : R → R d is deemed the state solution if it satisfies the ODE system dy dt = Φ(y) and the initial condition y(t 0 ) = y 0 . Given a time range [t 0 , T ], y(T ) is denoted as the end state of the interval.

Theorem 1 ( [3]

). Consider an NPE problem, N P E(Ω, G), and let Assumption 1 hold. For any initial condition, (t 0 , y 0 ), the state solution of the ODE system (6) converges to the optimal solution of N P E(Ω, G) as time t approaches infinity, i.e, lim t→∞

y(t) = x * , (8)
where x * is an optimal solution of N P E(Ω, G).

In particular, if N P E(Ω, G) contains only one optimal solution x * , then the ODE system is globally asymptotically stable at x * . Initial value problem (IVP) construction. In practice, to solve the ODE system of Eq. ( 6) requires the construction of an IVP, which consists of the following three:

• The ODE system, which is given in Eq. ( 6).

• An initial condition (t 0 , y 0 ), which is user-specified and represents the starting point of the state solution.

• A time range [t 0 , T ], which is user-specified and represents the domain of the state solution to be solved. After constructing the IVP, y(t), t ∈ [t 0 , T ] represents the state solution over the time range of [t 0 , T ]. According to Theorem 1, the end state, y(T ), serves as a predicted solution to the NPE, i.e., y(T ) ≈ x * . The state solution is typically determined by a numerical method, such as the Runge-Kutta method. In the next section, we show how to use the neural network to solve the IVP, thus avoiding the use of numerical methods.

III. REFORMULATION OF AN NPE AS A NEURAL NETWORK LEARNING PROBLEM

In this section, we illustrate the process of reformulating an NPE problem as a neural network training problem. Section III-A presents a modified PINN and shows how it solves for both the ODE system and the NPE. Section III-B presents the loss and objective functions used to train the neural network. Section III-C examines the error composition of the proposed method and discusses its convergence properties.

A. PINN Model

Model description. We use the following PINN model to solve the NPE problem, ŷ(t; w) = y 0 + 1 -e -(t-t0) N (t; w), t ∈ [t 0 , T ], [START_REF] Liu | A one-layer recurrent neural network with a discontinuous hard-limiting activation function for quadratic programming[END_REF] where N (t; w) denotes a fully connected neural network with trainable parameters w. The given time range is [t 0 , T ]. The auxiliary function 1 -e -(t-t0) ensures that the neural network always satisfies the initial condition ŷ(t = t 0 ; w) = y 0 , irrespective of the model parameters w. Fig. 1(A) describes the PINN model.

Construction technique. The designed PINN model, as described by Eq. ( 9), employs the construction technique proposed by Lagaris et al. [START_REF] Lagaris | Artificial neural networks for solving ordinary and partial differential equations[END_REF]. This technique aims to modify the output of the neural network N (t; w) so that it inherently satisfies initial or boundary conditions, independent of the values of the trainable parameters w. This construction technique has also been used in other PINN studies [START_REF] Lu | Physics-Informed Neural Networks with Hard Constraints for Inverse Design[END_REF], [START_REF] Pang | fpinns: Fractional physicsinformed neural networks[END_REF], [START_REF] Lagari | Systematic construction of neural forms for solving partial differential equations inside rectangular domains, subject to initial, boundary and interface conditions[END_REF]. In our work, we adapt this technique to the specific requirements of our problem setting, which considers only the initial conditions (t 0 , y 0 ).

Approximate state solution to the ODE system. As shown in Fig. 1 (B), the proposed model itself approximates the state solution to the ODE system (6) over the time range

[t 0 , T ], that is, ŷ(t; w) ≈ y(t), t ∈ [t 0 , T ], (10)
where y(•) represents the true state solution of the ODE system. Although the input time t of the model ŷ(t; w) can be any real number, we only regard ŷ(t; w) as the solution of the ODE within the time range [t 0 , T ]. Thus, we constrain the input to t ∈ [t 0 , T ]. The PINN model as an approximate state solution to the ODE system, where the end state marked with red stars is the prediction to the NPE.

(A) (B)

Predicted solution to the NPE. The end state of the proposed model, denoted by ŷ(t = T ; w), serves as the predicted solution to the NPE (2), as illustrated in Fig. 1 (B). The following equation demonstrates how ŷ(t = T ; w) approximates the optimal solution x * :

ŷ(t = T ; w) ≈ y(T ) ≈ x * , (11)
where ŷ(t = T ; w) ≈ y(T ) indicates that the end state of our model approximates the true end state, and y(T ) ≈ y * is derived from Theorem 1, indicating that y(T ) resolves the NPE.

B. Training Objective

Loss function. The loss function for the proposed neural network is defined as:

L(t, w) = ∂ ŷ(t; w) ∂t -Φ(ŷ(t; w)) , (12)
where Φ(•) refers to the ODE system [START_REF] Burden | Numerical analysis[END_REF], which corresponds to the NPE to be solved. Φ(ŷ(t; w)) represents the expected derivative according to the ODE system. ∂ ŷ(t;w) ∂t denotes the actual derivative of the model, which can be computed using automatic differentiation tools such as PyTorch or JAX [START_REF] Paszke | PyTorch: An imperative style, highperformance deep learning library[END_REF], [START_REF] Bradbury | JAX: composable transformations of Python+NumPy programs[END_REF]. L(t, w) expresses the difference between the two at time t and with network parameters w.

Incorporating the NPE into the loss function. First, the NPE is reformulated as an ODE system using neurodynamic optimization. This ODE system is then integrated into the loss computation process. Initially, a neural network is a versatile framework without a specific aim to solve a particular NPE. By incorporating the reformulated NPE as an ODE system into the loss function, the neural network is trained to solve both the ODE system and the NPE.

Objective function. Assuming a fixed neural network architecture, the objective function is defined as:

J(w) = T t0 ∂ ŷ(t; w) ∂t -Φ(ŷ(t; w)) dt, (13)
which is the integral of the loss function over the time range [t 0 , T ]. The loss value L(t, w) represents the error of the model at time t, while the objective function J(w) represents the total error of the model over the time range [t 0 , T ]. Batch loss. However, the objective function J(w) is computationally intractable due to its integral component. Hence, in practice, the model is trained by minimizing the following batch loss:

L(T, w) = 1 |T| t∈T L(t, w), (14)
where T is a set of randomly sampled times from the interval [t 0 , T ], and |T| denotes the size of the set. In this manner, the integral in the objective function J(w) can be approximated by a sum of loss values over the set of sampled times. By minimizing the batch loss, the model can be effectively trained to solve the NPE.

C. Error Analysis

Introduction.

In this subsection, we analyze the prediction error of the PINN solution, denoted by ŷ(T ; w), as given in Eq. [START_REF] Xia | A projection neural network and its application to constrained optimization problems[END_REF]. We discuss the various components that contribute to this error. We then examine how the choice of network architecture and time range [t 0 , T ], two important hyperparameters, affect the error composition. Finally, we show that there exists a neural network that can accurately solve the target NPE problem, supported by universal approximation theorems of neural networks and the global convergence theorem of neurodynamic optimization.

Notations Setup. Fig. 2-(A) illustrates the error decomposition for the PINN solution to an NPE problem. The mathematical notations used in the figure are explained below:

• H denotes a network architecture, represented by a set of neural networks with the same architecture. For example, for an architecture with a single hidden layer consisting of 100 neurons, H contains all the neural networks under that particular architecture. • x * is the optimal solution for the NPE being solved. • ŷ(T ; w * ) is the best achievable solution within the architecture H, where w * is the minimizer of the objective function [START_REF] Che | A two-timescale duplex neurodynamic approach to biconvex optimization[END_REF]. Error Decomposition. We define the total error between the predicted and the optimal solution as

ℓ total = |ŷ(T ; w) -x * |. (15)
As shown in Fig. 2-(A), the total error ℓ total consists of the following three components:

ℓ total = ℓ oe + ℓ ge + ℓ ne . (16)
Below, we explain each of the error components.

• Optimization Error is defined as ℓ oe = |ŷ(T ; w)ŷ(T ; w * )|. The optimization error captures the discrepancy between the actual PINN solution and the optimal solution achievable under the same neural network architecture H. This error is due to the inherent limitations of the training process in reaching a global optimum, often settling for a local one instead. This type of error is ubiquitous in various deep learning tasks and has been the subject of extensive research [START_REF] Chizat | On the global convergence of gradient descent for over-parameterized models using optimal transport[END_REF], [START_REF] Bonnans | Large-scale nonconvex optimization: randomization, gap estimation, and numerical resolution[END_REF]. • Generalization Error is given by ℓ ge = ∥ŷ(T ; w * )y(T )∥. The generalization error is due to the representational capacity of the chosen network architecture H. For example, consider a single-layer network with 100 neurons. Even if the global optimum w * is reached, there remains an unavoidable gap to y(T ).

• Neurodynamic Error is denoted by ℓ ne = ∥y(T ) -x * ∥.
This error is independent of the neural network and is due solely to the limited size of the chosen time range [t 0 , T ].

According to Theorem 1, as T approaches infinity, y(T ) converges to the optimal solution, i.e., lim T →∞ y(T ) = x * . In practice, however, one can only choose a finitely large [t 0 , T ], which leads to this error. Network Architecture and Time Range. As described above, the two hyperparameters, the network architecture H and the time range [t 0 , T ], play a critical role in the composition of the error. The implications of the choice of these hyperparameters are discussed below:

• The choice of the network architecture H affects the generalization error. As the complexity of the network architecture increases, characterized by an increase in the number of hidden layers and neurons, the upper bound on the representational capability of the network increases, thereby reducing the generalization error. However, a more complex architecture poses training challenges and increases the difficulty of minimizing the optimization error. • The choice of the time range [t 0 , T ] affects the neurodynamic error. For a fixed network architecture H, a larger time range tends to reduce the neurodynamic error, while potentially increasing the generalization error. The reduction of the neurodynamic error is due to the convergence properties outlined in Theorem 1. Conversely, the increased generalization error arises because a more complex network is required to accurately represent the true state solution y(t) for t ∈ [t 0 , T ]. Towards Zero ℓ ge and ℓ ne As shown in Fig. 2-(B), in an extreme scenario, where the network architecture tends towards infinite complexity and the time range [t 0 , T ] becomes infinitely large, both the generalization error ℓ ge and the neurodynamic error ℓ ne converge to zero. This focuses the objective solely on minimizing the optimization error, i.e., ℓ total = ℓ oe . In this context, the NPE problem transforms into a conventional neural network training task, similar to various deep learning challenges. Advanced research in neural network training can thus be leveraged to better address the NPE problem. Despite the theoretical appeal, training a neural network under these conditions introduces inherent challenges. In particular, although ℓ ge and ℓ ne tend to zero, ℓ oe will inevitably increase, complicating the optimization task.

Next, we demonstrate that within a single-layer fully connected network architecture, there exists a neural network that makes ŷ(T ; w) the optimal solution to the NPE problem. First, we give a formal definition of the network architecture under consideration, followed by a proposition establishing its existence, and finally, we conclude with some remarks. Definition 5. Consider a neural network that takes inputs over the interval [t 0 , T ] and contains a single hidden layer with k hidden units, d output units. We denote such a network architecture as C (k) (ψ), where ψ specifies the activation function of the hidden units. For a network with an arbitrary number of hidden units, the architecture is defined as:

C(ψ) = ∞ k=1 C (k) (ψ). ( 17
)
Assumption 2. The activation function of ( 17) is nonconstant, bounded, and continuous.

Proposition 3. Consider an NPE problem denoted as N P E(Ω, G), and x * is the optimal solution to this problem. Let Assumptions 1 and 2 hold. Consider any initial condition (t 0 , y 0 ) for the ODE system (6). Then, there exists a neural network N ∈ C(ψ) such that ŷ(T ; W ) is arbitrarily close to x * as T tends to infinity. Formally,

lim T →∞ ŷ(T ; w) = x * . (18)
Remark 1. The choice of the activation function significantly impacts the neural network's ability to globally approximate the optimal solution. As stated in Lemma 2, the activation function must be non-constant, bounded, and continuous. Consequently, the widely used rectified linear unit (ReLU) activation function is unsuitable due to its unbounded nature. Instead, the hyperbolic tangent (tanh) activation function is a more appropriate choice, as it fulfills all the three requirements.

Remark 2. Proposition 3 establishes the existence of a singlelayer neural network that allows the PINN model [START_REF] Liu | A one-layer recurrent neural network with a discontinuous hard-limiting activation function for quadratic programming[END_REF] to provide an optimal solution to the target NPE problem. This existence result is jointly derived from the global convergence theorem in neurodynamic optimization and the universal approximation theorem for neural networks. Specifically, the global convergence theorem ensures that the state solution y(t) converges to the optimal solution as time t tends to infinity. Subsequently, the universal approximation theorem guarantees that there exists a neural network capable of approximating the corresponding y(t).

IV. MULTI-TASK LEARNING FOR MULTIPLE NPES In this section, we introduce the multi-task learning (MTL) framework employed in our study, which allows us to address multiple NPE instances using a single neural network, thereby eliminating the need for repeated individual solutions. In Section IV-A, we present the MTL framework, the loss function, and the training objective. In Section IV-B, we present the training procedure of the MTL model.

A. Multi-Task Learning Framework

MTL framework description. We employ an MTL framework based on the model in ( 9) that facilitates the simultaneous training of a single model for solving multiple NPEs. As shown in Fig. 3, the model comprises two parts: shared layers and task-specific layers. The shared layers learn common features among NPEs, promoting knowledge transfer across different NPEs and improving generalization capabilities. The task-specific layers generate task-specific predictions for each NPE.

Model input and output. The MTL model can be expressed by the following equation:

M T L(t; w) = ŷ(1) (t; w), ŷ(2) (t; w), . . . , ŷ(n) (t; w) . ( 19
)
The model takes an input of t ∈ [t 0 , T ]. It is essential to note that the MTL model requires the same time range [t 0 , T ] for all NPEs, as they share the input. The model has n outputs, each of which corresponds to an NPE to be solved. The MTL solutions for the multiple NPEs can be directly generalized from Section III. Specifically, each output, ŷ(i) for t ∈ [t 0 , T ], is the approximated state solution for the i-th ODE system, which in turn is derived from the i-th NPE. When choosing the time to t = T , the output ŷ(i) (t = T ; w) is the prediction for the i-th NPE. For each NPE, we can provide a distinct initial condition (t (i) 0 , y (i) 0 ). MTL loss. We now define the loss function for the MTL model consisting of n NPEs. The MTL loss function is given by:

L MTL (t, w) = 1 n n i=1 ∂ ŷ(i) (t; w) ∂t -Φ (i) (ŷ (i) (t; w)) , (20)
where Φ (i) (•) represents the ODE system derived from the i-th NPE. The MTL loss in [START_REF] Xiong | Tdoa-based localization with nlos mitigation via robust model transformation and neurodynamic optimization[END_REF] essentially extends the loss function [START_REF] Guo | A one-layer recurrent neural network for pseudoconvex optimization subject to linear equality constraints[END_REF] to accommodate multiple NPEs. Similarly, we formulate the objective function and the batch loss for the MTL model as follows:

J MTL (w) = 1 n n i=1 T t0 ∂ ŷ(i) (t; w) ∂t -Φ (i) (ŷ (i) (t; w)) dt, (21)
L MTL (T, w) = 1 |T| t∈T L MTL (t, w). (22)
The MTL loss function allows the model to learn from multiple NPEs simultaneously, thereby enhancing its overall performance and generalization capabilities. It is important to note that in the loss functions [START_REF] Xiong | Tdoa-based localization with nlos mitigation via robust model transformation and neurodynamic optimization[END_REF], [START_REF] Dissanayake | Neural-networkbased approximations for solving partial differential equations[END_REF], and ( 22), we assume that the weights assigned to each NPE are equal. In practice, one has the flexibility to assign weights as needed.

Shared information between tasks. The shared layers of the MTL framework output a hidden feature vector in their final layer. This hidden feature vector serves as a representation of the shared information between multiple NPE tasks. However, due to the black-box nature of neural networks, it is challenging to interpret or explain the specific meaning encapsulated by this hidden feature vector. This is analogous to the difficulty of interpreting a single feature vector or feature

Input

Shared layers Task-specific layers ), ŷ(i) , and Φ(•) (i) represent the initial condition, the approximate state solution, and the corresponding ODE system, respectively. M.S.E. stands for mean square error. map in any neural network, as its semantics are generally not readily understandable.

B. Training Procedure

NPE error. To assess the effectiveness of a prediction x pred in solving the NPE (2), we introduce a performance metric called the NPE error. This metric is defined as follows:

NE(x pred ) = ∥P Ω (x pred -G(x pred )) -x pred ∥ ∞ . (23)
The NPE error quantifies the deviation of the prediction x pred from the true solution of the NPE, providing an accurate assessment of how well the model has learned to solve the given problem. By minimizing the NPE error during training, we ensure that the model's predictions closely adhere to the properties of the NPE and improve its overall performance.

Pipeline. Algorithm 1 outlines the process for utilizing the proposed MTL framework to tackle multiple NPEs concurrently. Let I = {1, 2, . . . , n} denote the index set for the NPEs. First, a time range [t 0 , T ] and a set of initial conditions {(t (i) 0 , y (i) 0 )} i∈I must be given to formulate the IVPs. Then, the MTL model is initialized with n task-specific heads, where each head ŷ(i) (t; w) corresponds to the initial condition (t

(i) 0 , y (i) 0
). The model is fine-tuned by applying gradient descent to the loss function [START_REF] Lagaris | Artificial neural networks for solving ordinary and partial differential equations[END_REF] to improve the approximation. Notably, our solver relies exclusively on deep learning infrastructure, eliminating the need for traditional optimization or numerical integration solvers.

Best Prediction w.r.t NPE error. An important feature of Algorithm 1 is that it uses the NPE error [START_REF] Mcfall | Artificial neural network method for solution of boundary value problems with exact satisfaction of arbitrary boundary conditions[END_REF] as a metric to evaluate how well the prediction solves the NPE problem. At each iteration, the algorithm calculates the NPE error for the current iteration, denoted as NE (i) curr , and compares it to the smallest NPE error encountered so far, denoted as NE (i) best . The terms x (i) curr and x (i) best represent the prediction at the current iteration and the best prediction encountered so far, respectively. If NE (i) curr is less than NE

(i)
best , it indicates that the model has identified a better prediction during this iteration. As a result, the algorithm updates NE (i) best to equal NE (i) curr and retains the best prediction as x

(i) best = x (i)
curr . Such a mechanism ensures that the optimal prediction obtained by the model is maintained throughout the training process, ultimately improving the overall performance of the algorithm. Algorithm 1 presents the for-loop iteration over the index set I, which in practice can easily be replaced by a parallel approach using Numpy to enhance computational efficiency.

V. EXPERIMENTS

Experimental Setup. We use PyTorch 1.12.1 [START_REF] Paszke | PyTorch: An imperative style, highperformance deep learning library[END_REF] to implement the proposed MTL framework and JAX 0.4.1 [START_REF] Bradbury | JAX: composable transformations of Python+NumPy programs[END_REF] to implement the ODE system. Regarding the network architecture for the MTL framework, the shared layers consist of three fully connected layers of 50, 100, and 100 neurons, respectively. The task-specific layers contain two fully connected layers with 100 and 50 neurons, respectively. The activation function used is tanh. The optimizer used for training is Adam [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF], with a learning rate of 0.001 and a batch size of 128. All three subsections use t 0 = 0 and y 0 = 0 as the initial condition. Sections V-A and V-C choose [0, 10] as the time range, while Section V-B considers different choices of time ranges. The numerical solvers used for comparison are Runge-Kutta 45 (RK45), Runge-Kutta 23 (RK23), Dormand-Prince 853 (DOP853), Backward Differentiation Formula (BDF), Radau, and LSODA. All of these methods are available in the Scipy library [START_REF] Virtanen | [END_REF].

Evaluation Metrics. We evaluate the performance of the MTL framework using the following three metrics.

• NPE Error: As defined in Eq. ( 23), this metric measures how well a prediction solves the NPE problem. n number of NPEs: {N P E(Ω (1) , G (1) ), N P E(Ω (2) , G (2) ), . . . , N P E(Ω (n) , G (n) )}; n number of initial conditions: {(t

(1) 0 , y (1)
0 ), (t Derive ODE systems, {Φ (1) (•), Φ (2) (•), . . . , Φ (n) (•)}, corresponding to the given NPEs using (6). for i = 1, 2, . . . , n do 10:

(2) 0 , y (2)
Obtain the prediction from the i-th head of the MTL model:

x (i) curr = ŷ(i) (t = T ; w) 11:
Project

x (i)
curr onto the feasible set Ω using (1):

x (i) curr = P Ω (x (i) curr ) 12:
NE (i) curr = NE(x (i) curr ) ▷ Calculate the NPE error of x (i)
curr by [START_REF] Mcfall | Artificial neural network method for solution of boundary value problems with exact satisfaction of arbitrary boundary conditions[END_REF]. if

NE (i) curr < NE (i)
best then 14:

NE (i) best = NE (i) curr ▷ Update NE (i)
best .

15:

x

(i) best = x (i) curr
▷ Update the best prediction. return {x

best } 20: end function

A. Solving Multiple NPEs with the Proposed MTL Framework

Problem Definition. Consider the following NPE problem,

P Ω (x * -G r (x * )) = x * , (24)
where

r = [a, b, c, d] ∈ R 4 , G r (x) =     3x 1 -a x1+1 + 5x 2 -13 1.2x 1 + bx 2 cx 3 + 8x 4 1x 3 + 2x 4 -4 x4+2 -d     , Ω = {x ∈ R 4 | 1 ≤ x 1 ≤ 100, -3 ≤ x 2 ≤ 100, -10 ≤ x 3 ≤ 100, 1 ≤ x 4 ≤ 100}, (25)
and P Ω (•) is defined in [START_REF] Che | A collaborative neurodynamic approach to global and combinatorial optimization[END_REF]. r is the problem data, and selecting distinct values of r yields different NPE instances.

Problem Set Construction. We create a problem set by uniformly sampling r from the interval [1, 10] 4 . The set of the multiple NPEs is denoted as {N P E(G r1 , Ω), N P E(G r2 , Ω), . . . , N P E(G rn , Ω)}, where n is the number of instances. These NPE instances share the same feasible region Ω and a majority of the nonlinear function G r (•), distinguished only by the problem parameter r. Traditional numerical methods solve this set of instances individually, treating each instance as an independent problem. In contrast, our proposed MTL framework solves all these NPE instances simultaneously in a parallel manner.

Training of the MTL Framework. Fig. 4 shows the training of different numbers of NPE instances using our proposed MTL framework. Specifically, Fig. 4(A) displays the loss values, indicating the solution accuracy to the ODE systems, while Fig. 4(B) shows the NPE error, indicating the solution accuracy to the NPE instances. With only 1000 iterations, all MSE losses converge from about 150 to less than 1, and all NPE errors converge from about 13 to less than 1. Notably, the significant reduction in MSE loss occurs primarily between the 100th and 1000th iteration, while the reduction in NPE error occurs between the 10th and 100th iteration. The asynchronous decrease of these two metrics suggests that the NPE error can be reduced even if the MSE loss is not reduced or even increased. These results benefit from the use of the evlautation metric [START_REF] Mcfall | Artificial neural network method for solution of boundary value problems with exact satisfaction of arbitrary boundary conditions[END_REF] in Algorithm 1.

Performance of the MTL Framework. Fig. 5 demonstrates the computational efficiency of the MTL framework compared to traditional numerical solvers. Here, the number of NPE instances is fixed at 100. Experimental results show that our MTL framework outperforms the numerical solvers across various target accuracies. In particular, when the required target accuracy is relaxed, such as requiring an NPE error of 1, the computational efficiency advantage of the MTL framework is significant, at least 3 times faster than RK45 (the bestperforming numerical solver in this case). When the required target accuracy is stringent, such as requiring an NPE error of 0.05, the computational efficiency of the MTL framework still leads over various numerical solvers.

B. Application to Hopfield Network

The goal of this subsection is to apply the proposed MTL framework for solving equilibrium points in a Hopfield network. This subsection also serves as a hyperparameter study to

Average MSE loss

Average NPE error

where T ∈ R 10×10 is defined as

T =        3 3 . . . 3 3 3 . . . 3 . . . . . . . . . . . . 3 3 . . . 3 -3 -3 . . . -3        . ( 27
)
The feasible set Ω is defined as

Ω = {x ∈ R 10 | -1 ≤ x ≤ 1}. ( 28
)
Our goal is to find x * such that dx dt = 0. The equilibrium point computation of the Hopfield network has wide applications in various types of combination optimization problems [START_REF] Hopfield | Neural" computation of decisions in optimization problems[END_REF]- [START_REF] Paik | Image restoration using a modified hopfield network[END_REF].

Problem Set Construction. We construct a set of instances by selecting several different q ∈ R 10 in Eq. [START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF]. Specifically, we uniformly sample 50 different q from the interval [-1, 1] 10 to create a problem set. Next, we use the proposed MTL framework to solve these 50 instances in parallel.

Training of the MTL Framework. Selection of Time Range. As shown in Fig. 7, the choice of the time range significantly influences the convergence speed of the MTL framework. There is no a universally optimal choice for the time range. Instead, it should be determined based on the training resources allocated. For example, in Fig. 7, the time range [0, 5] gives the best results for 750 iterations, whereas [0, 10] gives the best results for 1500 iterations. Here are some guidelines on how to choose an appropriate time range based on the number of training iterations:

• If computational resources are limited and a quick result is desired, choose a smaller time range, such as [0, 5]. • If computational resources are sufficient, a larger time range should be chosen, as this will result in reduced neurodynamic errors, consistent with the discussion in Section III-C. • In most scenarios, the choice of time range should be customized to align with the available computational resources. For example, we find that the time range [0, 10] works best for 1500 training iterations.

C. Application to Variational Inequality

The goal of this subsection is to apply the proposed MTL framework for solving variational inequalities. This subsection also compares the proposed MTL framework with various iterative methods.

Problem Definition. Consider the following variational inequality:

(x -x * ) T G r (x * ) ≥ 0, x ∈ Ω. (29)
G r (x) is parameterized by r = [a, b, c, d] ∈ R 4 , and

G r (x) =      2x 1 e x 2 1 +(x2-1) 2 + x 1 -x 2 -x 3 + a 2(x 2 -1)e x 2 1 +(x2-1) 2 -x 1 + 2x 2 + 2x 3 + b -1 √ x1+x2+3 * x3+x4 + x 1 + x 2 + x 3 + x 4 + c x 3 + 2x 4 -4 x4+2 + d      . ( 30
) Ω is defined as

Ω = {x ∈ R 4 | x ≥ 0}. (31)
By Proposition 2, Eq. ( 29) is reformulated as an NPE problem. Problem Set Construction. We construct a problem set by choosing different vectors r in Eq. ( 29), and each r corresponds to a particular instance. Specifically, we sample r with a uniform distribution over the interval [-10, 10] 4 . In the following, we apply our proposed MTL framework to solve the problem set with different numbers of instances. Description of Fig. 8. Fig. 8 shows the total computation time required by different solvers to solve a varying number of instances and at different target accuracies. In the figure ,   • 'MTL' denotes the MTL framework proposed in this paper.

• 'NN' refers to solving multiple NPE problems iteratively and individually using the PINN model [START_REF] Liu | A one-layer recurrent neural network with a discontinuous hard-limiting activation function for quadratic programming[END_REF], where one PINN model corresponds to one instance. • 'RK45' and 'BDF' are explicit and implicit numerical solvers, respectively. • 'LSODA' is a numerical solver that flexibly switches between the above two solvers.

Note that only the proposed MTL framework is a parallel solution algorithm, the other four methods are iterative.

Comparative Efficiency Analysis of Solvers. Based on the analysis of our experimental results, we draw the following conclusions:

• For all NPE error thresholds, the computational time of our proposed MTL framework increases as the number of instances increases. However, its growth rate is significantly slower than that of other solvers, indicating that MTL is more efficient when dealing with large numbers of instances. • While the NN method performs comparably to the MTL framework with fewer instances, its computation time increases more rapidly as the number of instances increases, suggesting that it is less adept at handling large instance sets than MTL. • Conventional numerical solvers, RK45, BDF, and LSODA, show relatively consistent performance over different NPE error thresholds. However, their computation times are higher than those of the MTL and NN solvers. In particular, LSODA performs relatively better than the other two numerical methods for this problem.

• As the NPE error threshold decreases, i.e. as the solution accuracy requirement becomes more stringent, the computation time increases for all solvers. This indicates that higher accuracy requirements place greater demands on computational resources. In this scenario, our MTL solver continues to maintain a leading position.

VI. CONCLUSION In this paper, we have proposed a novel methodology for solving NPEs. The proposed approach employs three important tools: (1) neurodynamic optimization for modeling NPEs as ODE systems, (2) PINNs for solving the ODE systems, and (3) MTL framework for integrating multiple PINNs. In addition, we have conducted a comprehensive error analysis on the proposed method. Experimental results show that our approach is efficient and accurate, outperforming traditional numerical integration methods. The proposed method holds great promise for efficiently addressing various constrained nonlinear optimization problems.

Below we suggest some future directions for our work:

• Investigate advanced network architectures or training methods from the field of PINNs to further improve the performance of the proposed approach. • Investigate advanced MTL frameworks to improve or understand the shared knowledge across multiple NPEs.

• Adapt the proposed approach to solve other types of nonlinear optimization problems by collaborating with various neurodynamic techniques.

A. APPENDIX

In this Appendix, we provide a proof for Proposition 3. The structure of the appendix is organized as follows: A.1 and A.2 contain two lemmas concerning the universal approximation theorem for neural networks. Subsequently, in A.3, we employ these two lemmas to furnish the proof for Proposition 3.

The following metric is used to measure the closeness between two functions. Definition 6. f : [t 0 , T ] → R d and g : [t 0 , T ] → R d are two functions. The closeness between the two functions is defined as follows

ρ(f, g) = sup t∈[t0,T ] ∥f (t) -g(t)∥ ∞ , (32)
where ∥•∥ ∞ is the infinity norm.

A.1 Lemma 1 and its proof

Proof. Let C([t 0 , T ], R n ) be the space of all continuous functions from the interval [t 0 , T ] to R n , equipped with the metric ρ defined in [START_REF] Sharma | Accelerated Training of Physics Informed Neural Networks (PINNs) using Meshless Discretizations[END_REF]. Define the closure of C(ψ), denoted as Ĉ, as the set of all functions that can be uniformly approximated by functions from C(ψ) with arbitrary precision. In other words, if h ∈ Ĉ, then for any ϵ > 0, there exists a neural network N ∈ C(ψ) such that ρ(h, N ) < ϵ.

Since ψ(•) is non-constant, bounded, and continuous, Ĉ contains non-constant functions and is closed under addition and scalar multiplication. Furthermore, Ĉ separates points. By the Stone-Weierstrass theorem, Ĉ is dense in C([t 0 , T ], R n ). Now, fix y ∈ C([t 0 , T ], R n ) and ϵ > 0. Since Ĉ is dense in C([t 0 , T ], R n ), there exists a function h ∈ Ĉ such that ρ(h, y) < ϵ/2.

Since h ∈ Ĉ, there exists a neural network N ∈ C(ψ) that approximates h within an arbitrary degree of accuracy. Choose a neural network N such that ρ(h, N ) < ϵ/2. Therefore, we have: ρ(N, y) ≤ ρ(y, h) + ρ(h, N ) < ϵ/2 + ϵ/2 = ϵ This completes the proof, showing that a feedforward neural network with one hidden layer and a finite number of neurons can approximate the continuous function y with an error less than ϵ under the closeness metric ρ(•, •).

A.2 Lemma 2 and its proof Lemma 2. Let y : [t 0 , T ] → R d be a continuous function and choose a non-constant, bounded, and continuous activation function ψ(•). Then, for any ϵ > 0, there exists a neural network N ∈ C(ψ) such that ρ(y, ŷ) < ϵ, [START_REF] Lu | DeepXDE: A Deep Learning Library for Solving Differential Equations[END_REF] where ŷ : [t 0 , T ] → R n is the PINN model in Eq. (9).

Proof. First, we construct a function z(t) = y(t) -y 0 to be approximated. By Lemma 1, there exists a neural network N ′ ∈ C(ψ) that approximates z with an arbitrary degree of accuracy. That is, for any ϵ > 0, there exists N ′ ∈ C(ψ) such that

ρ(z, N ′ ) < ϵ 2 . ( 35
)
Now, let's define a new function N (t; w) = N ′ (t;w)

(1-e -(t-t 0 ) ) .

Note that since N ′ is a neural network with activation function ψ, it follows that N is also a neural network with activation function ψ, and thus N ∈ C(ψ).

We construct the neural network-based function ŷ(t; w) = y 0 + 1 -e -(t-t0) N (t; w). Now we need to show that ρ(y, ŷ) < ϵ.

We have:

ρ(y, ŷ) = ρ(y, y 0 + 1 -e -(t-t0) N ) = ρ(y 0 + z, y 0

+ N ′ ) = ρ(z, N ′ ) < ϵ 2 < ϵ. (36)
This completes the proof, showing that there exists a PINN model ŷ(t; w) which approximates the continuous function y with arbitrarily small ϵ under the closeness metric ρ(•, •).

A.3 Proof of Proposition 3

Proposition 3. Consider an NPE problem denoted as N P E(Ω, G), and x * is the optimal solution to this problem. Let Assumptions 1 and 2 hold. Consider any initial condition (t 0 , y 0 ) for the ODE system (6). Then, there exists a neural network N ∈ C(ψ) such that ŷ(T ; W ) is arbitrarily close to x * as T tends to infinity. Formally,

lim T →∞ ŷ(T ; w) = x * . ( 37
)
Proof. Consider the ODE system (6) that models N P E(Ω, G), and let the state solution of the ODE system be y : [t 0 , +∞) → R d . By Theorem 1, we have

lim t→∞ y(t) = x * , (38)
where x * is an optimal solution of N P E(Ω, G). Let ϵ > 0 be an arbitrary positive number representing the desired accuracy. Consider an increasingly large sequence of time ranges {[t 0 , T n ]} with T n → +∞ as n → ∞.

For each interval [t 0 , T n ], apply Lemma 2 to obtain the PINN model ŷn (t; w n ) that approximates y(t) on the interval [t 0 , T n ] with an error less than ϵ, i.e., ρ(ŷ n , y) < ϵ.

(

) 39
By the definition of ρ(•, •) in Eq. ( 32), at the end state t = T n , it holds that ∥ŷ n (t = T n ; w n ) -y(t = T n )∥ < ϵ (40)

As n → ∞, T n → +∞. Therefore, lim n→∞ ŷn (t = T n ; w n ) = lim t→∞ y(t).

Combining ( 41) with [START_REF] Sener | Multi-task learning as multi-objective optimization[END_REF], we conclude that lim n→∞ ŷn (t = T n ; w n ) = x * . [START_REF] Zhang | A survey of multi-task learning in natural language processing: Regarding task relatedness and training methods[END_REF]

Fig. 1 :

 1 Fig. 1: (A): The PINN model, where [N 1 , N 2 , . . . , N n ] is the output of the neural network, and [ŷ 1 , ŷ2 , . . . , ŷn ] is the output of the model. (B):The PINN model as an approximate state solution to the ODE system, where the end state marked with red stars is the prediction to the NPE.

Fig. 2 :

 2 Fig. 2: (A) Error decomposition showcasing optimization, generalization, and neurodynamic errors. (B) The extreme case of infinitely large network architecture and time range.

Fig. 3 :

 3 Fig.3: MTL model for solving multiple NPEs. This figure illustrates the process of solving NPEs using the MTL framework. For each NPE indexed by i = 1, 2, . . . , n, (t(i) 0 , y (i) 0), ŷ(i) , and Φ(•)(i) represent the initial condition, the approximate state solution, and the corresponding ODE system, respectively. M.S.E. stands for mean square error.

3 :

 3 Instantiate the MTL model based on the given initial conditions.

  best = NE(ŷ (i) (t = T ; w)).

5 : 7 :Forward propagation 8 :

 578 while iteration ≤ maximum iteration do 6: Sample T ∼ U (t 0 , T ) ▷ U (t 0 , T ) is the uniform distribution with interval [t 0 , T ] Compute the MTL batch loss: L MTL (T, w) ▷ Update weights: w ← ∇ w L MTL (T, w) ▷ Backward propagation 9:

  13:

Fig. 4 :Fig. 5 :

 45 Fig. 4: Training of the MTL framework. n denotes the number of NPE instances to be solved. (A) Average MSE loss versus training iteration. (B) Average NPE error versus training iteration.

Fig. 6 Fig. 6 :Fig. 7 :

 667 Fig. 6: NPE error versus training iteration for the 50 instances of Eq. (26). (A): the NPE errors from iteration 0 to 600. (B): the NPE errors from iteration 600 to 1500.Iteration 0~750

Fig. 8 :

 8 Fig. 8: Efficiency comparison of various solvers over different NPE error thresholds. The y-axis indicates the total computation time in seconds, while the x-axis indicates the number of instances.

Lemma 1 .

 1 Let y : [t 0 , T ] → R d be a continuous function, and choose a non-constant, bounded, and continuous activation function ψ(•). Then, for any ϵ > 0, there exists a neural network N ∈ C(ψ) such that ρ(y, N ) < ϵ.

  )

	Proposition 1 ( [1]). Let Ω ⊂ R d be a nonempty closed
	convex set. Then x * solves the problem N CP (G) if and only
	if x * solves N P E(R n + , G), where R d + = {x ∈ R d |x ≥ 0}
	represents the set of non-negative real vectors.
	Proposition 2 ( [1]). Let Ω ⊂ R d be a nonempty closed convex
	set. Then x * solves the problem V I(Ω, G) if and only if x *
	solves N P E(Ω, G).

  Training of the MTL model to solve multiple NPEs Input: A time range: [t 0 , T ];

• Mean Square Error (MSE) Loss: This metric measures how well a predicted state solution solves the ODE system.

• Computation Time: This metric measures the computational efficiency.

Algorithm 1

  Predictions for the NPEs, denoted by {x

	Output: (1) best , x (2) best , . . . , x	(n) best }.	(n) 0 , y 0 )}. (n)
	1: function NPES SOLVER
	2:

0 ), . . . , (t