Wadih Sawaya
email: wadih.sawaya@telecom-lille1.eu

Stéphane Derrode
email: stephane.derrode@centrale-marseille.fr

Mustava Ould-Barikalla

Jacques Rivailler
email: jacquesrivaillier@orange.fr

Jacques Rivailler Detection

M Ould

Detection and Iterative Decoding of a 2D Alphabetic Barcode

An algorithm dedicated to the detection and decoding of an alphabetic 2D barcode, namely the Alphacode, from a captured image is developed. When the camera is freely handled, images are degraded by several artifacts (motion, acquisition, low resolution...), making recovery of the code a more difficult task. We propose a two-stage strategy. The code structure is first recovered using image processing steps including localization, geometric correction and code grid structure retrieval. The code is then recognized using an iterative structure based on a factor graph representation of the displaying/reading process and making use of the specific structure of the code. The proposed algorithm is evaluated in case of strongly blurred images and shows high improvements in performances compared to the basic processing which is unable to decode reliably.

I. INTRODUCTION

ith the improvement of the quality and the resolution of their equipped camera, mobile phone are no more restricted to radio communication. Used as a capture device in conjunction with barcodes, mobile phones are addressed to interact with the real-world [START_REF] Rohs | Real-World Interaction with Camera-Phones[END_REF]. Emerging applications, for large public as well as professional usages, require drastic increasing in coding capacity and robustness with respect to the usual 1D barcode. Almost all solutions converge now to 2D barcodes meeting consequently the markets of identification (e.g. electronic stamps), authentication (e.g. counterfeit medicines) and interacting (e.g. e-ticketing).

Many 2D barcodes are now available such as Datamatrix, Maxicode, QR-code and PDF417 (see [START_REF] Ouaviani | A common image processing framework for 2D barcode reading[END_REF] for an example of the different symbols). Beside all of these codes, the Alphacode, patented by the eponymous society [START_REF]Principal patents for the Dotem code[END_REF], is built on a very different principle. Whereas, the aforementioned 2D barcodes represent a monolithic block of symbols generated by dedicated software, Alphacode is based on a predetermined alphabet defined by a font a characters, the Dotem font, and a given code consists with a chain of characters. The number of characters is not limited allowing to code as much information as needed for a given application. An example of an Alphacode word with five characters (two columns each) is shown in Fig. 1.a. This font based implementation grants Alphacode with numerous advantages especially in flexibility and interoperability. Standard keyboards, word processors and printers are the only materials required to generate a code on a physical support (screen, paper …). Remarkably enough, the Dotem alphabet is full compatible with Unicode norm (8, 16 and 32).

In order to confirm this novel technology for the identification, authentication and mobile markets we propose robust solutions for detecting and decoding the Alphacode. Most of the barcode readers are designed with limited signal processing circuitry [START_REF] Ouaviani | A common image processing framework for 2D barcode reading[END_REF]. The response time satisfying the human operator wins out over quality, or false decoding and rejections. Besides, classical artifacts inherent to low-quality cameras (e.g. sensor noise, sensor resolution) or to motion caused by camera manipulation by end-user or to unfocused image and to geometric distortion…impact drastically reading performances [START_REF] Houni | 1D barcode reading: an information theoretic approach[END_REF]. Classical processing techniques have a rate of good detection in the context described above at best of about 90% before any error correction blocs. This is definitely not sufficient when aiming at developing robust, reliable and sustainable receivers. In digital transmission, error probability is a quality measure used to evaluate different receiver designs. A value of 10 -6 is a commonly acceptable error rate for data. We propose then cutting edge solutions that enable us to attain this low error probability when decoding the Alphacode with standard cameras.

We propose a two-steps strategy which consists first in localizing the word and recovering the grid structure, second, in decoding using iterative structure based on statistical processing and making use of the Dotem structure.

Detection and Iterative Decoding of a 2D Alphabetic Barcode

Wadih Sawaya, Stéphane Derrode, Mustava Ould-Barikalla and Jacques Rivaillier The remaining of the paper is as follows. Section 2 presents the structure of a Dotem character. Section 3 exposes the image processing steps to localize the code and to correct geometrical distortions due to acquisition. Section 4 studies some optimal decoding strategies, taking into account image blur. Finally, Section 5 concludes on the proposed algorithm and discusses possible improvements.

II. DOTEM CODE STRUCTURE

In the basic Dotem alphabet presented in Fig. 1.b. (namely truetype font DoteM02), a character is made of at least two columns and ten rows of dots. A dot is a basic black or white square element designing a binary symbol. Dots from rows 2 to 7 are called "data dots" and allow to code up to 2 12 =4096 characters in a 2-column structure, 2 18 =242164 characters for a 3-column structure… Rows 1 and 10 have a fixed pattern used to make reading easier, especially the half dot in row 1, column 2. The black dot in row 9, last column (see Section 2.2) is an indicator for the end of a character. Dots from row 7 equip the character with a parity check binary element used to enhance the decoding process (Section 4).

III. IMAGE RECTIFICATION

When acquiring a 2D barcode with a mobile camera such as a webcam or a camera in a mobile phone, the so-obtained image is degraded by several artifacts [START_REF] Houni | 1D barcode reading: an information theoretic approach[END_REF] (see figure 2 for a typical example):

• Complex background: the image may contain other objects lying in the scene which gives arise to some difficulties in detecting and localizing the code in the image;

• Distortions: geometric distortions can be expected from the acquisition angle between the sensor and the code resulting in a projective deformation. A step is required to rectify the image in order to recover the rectangularshaped code.

• Noise and blur: Apart from the inherent sensor noise, images are blurred by (i) motion during acquisition, and (ii) necessary interpolation for rectification. Hence dots intensity runs on neighboring ones, making the decision on the state of a dot more difficult to take. A typical image illustrating the various artifacts encountered in barcode reading from a camera-equipped mobile phone. In this section we develop an algorithm for the recovery of the underlying grid from an image as in Fig. 2. The latter will be used to illustrate the three algorithm steps.

A. Code localization (Fig. 3)

To localize the region of interest inside a complex scene, we opt for the image gradient which is able to exhibit code edges. Then we label all connected components based on local gradient density. Using a priori geometrical property of the desired code, we select the region of interest among possible candidates.

B. Image rectification (Fig. 4)

As explained before, images are distorted by some projective transformation. In this step, we intend to estimate the geometric parameters and apply the projection backward to get the code rectified. The first step consists in retrieving the four corners of the parallelepiped including the code, see Fig. 4. a. We use the points of the convex hull to set the left, right, up and down segments, which give by intersection the searched corners.

Then, the parallelepiped corners are registered to a rectangle one to estimate the projection matrix. The projection matrix we get for image in Fig. 2 (see [START_REF] Rohs | Real-World Interaction with Camera-Phones[END_REF]) is

          - - = 1

C. Grid recovery (Fig. 5)

To recover the number of characters and the number of columns by character we make use of the different tags present in each character: (i) the semi-dot in the first row and second column, (ii) the black dot in the ninth row and last column. Basic image processing based on column and line profiles allow us to recover the code structure. For the studied example, the mean width of a dot is 18 pixels, and the mean height is about 18.7 pixels.

The processing implemented to recover the grid structure takes about 1.5 seconds for a PC with an Intel Core 2, 2.0 GHz processor, using Matlab.

IV. CODE READING

According to the grid the observed image is now a matrix Y of observed dots with elements:

l k n n i j i p p j j l i k l k w h u y , , , , + = ∑ ∑ - = - = - - . ( 1
)
Elements u k,l represent black and white modules (dots) of the original code, w k,l are zero-mean additive white Gaussian noise samples and the [h i,j ] describe the blur affecting the image. From a communication point of view, expression (1) models the reading/displaying process as a linear 2D ISI channel. A factor graph representation of this model is derived. We propose then a robust decoder based on a message-passing algorithm running on this graph and involving interactions between adjacent modules or dots. We make use of the parity dot to enhance performances with an iterative structure. We aim then at proposing a simplified version of the studied message passing algorithm without attempting much on performances, motivated by the design of receivers with a humanly acceptable complexity. We finally show some simulation results and conclude on the necessity of using such structure when simple binarization is not able to satisfy a given quality of service.

The factor graph consists of variable nodes (circles) and factor nodes (squares) connected by edges (Fig. 6.). Variable nodes represent the hidden data dots {u k,l }, the resulting gray-levels from ISI {x k,l } and the observed variables {y k,l }. Factor nodes r k,l represent local functions connecting a . The sum-product algorithm (SPA) is a message-passing algorithm that computes the marginal a posteriori probabilities (APP) of any data dot P(u k,l /Y), by passing "messages" through the factor graph [START_REF] Kschischang | Factor Graphs and the Sum-product Algorithm[END_REF].

We apply SPA to detect data dots (detection step). The outputs of this detector constitute soft extrinsic information about each dot and enter the parity decoder (decoding step). We realize an iterative detection and decoding procedure as the detector and the decoder exchange extrinsic information between each other for several iterations. In the SPA, messages are resumed by the following rules [START_REF] Kschischang | Factor Graphs and the Sum-product Algorithm[END_REF] for

n i n ≤ ≤ - et p j p ≤ ≤ -
(Fig. 6.):

1) Messages from variable u nodes to factor nodes r:

l k j i l k j l i k j l i k l k u D u r r u m m m , ) 0 , 0 ( ) ' , ' ( , ' , ' , , → → →         ∝ ∏ ≠ - - - - (2)
In the first formula we distinguish the message of the form m D→u coming from the decoder and calculated according to the extrinsic information that the latter deliver: )

= = = → → l k u D l k l k u D l k u m u P u m u P l k l k (3) 0 ( ) 0 ( ) 1 ( ) 1 ( , , , , , , = = =
Practically, all messages m r→u in (2) are initialized to 1 in order to avoid short cycles. We thus have

u D r u m m → → = .
2) Messages from factor nodes r to variable nodes u :

∏ ∑ ∑ ≠ - - - - → → → ⋅ ∝ ) , ( ) ' , ' ( , , , , , , , ' , ,
) , (

j i j i l k j l i k U l k x l k l k j l i k l k r u r x l k l k u r m m U x r m (4)
In expression (4) we define U ={u k-n,l-p ,…,u k+n,l+p }, -= U U ' {u k -i,l -j }, and l k r , is given by the following logical proposition:

        = = = ∑ ∑ - = - - - = n n i j l i k j i p p j l k l k l k u h x U x r , , , , , ) , ( (5)
such that r( x, U ) = [ x = = γ (U)] is defined by the following, r( x, U ) = 1 if x = γ (U ) and 0 if not. Finally, we have ) / ( , , , , ,

l k l k Y r x x y f m l k l k l k = →
. Assuming a white and Gaussian additive noise we rewrite (4) as:

( ) ( ) ∑ ∏ ∑ = → →         ⋅ - - ∝ ≠ - - - - 1 / ' 0 2 , , , ) , ( ) ' , ' ( ' , ' , , , exp l k j i j i j l i k l k j l i k l k r U u D x l k l k u r m N x y m (6)
The second sum in [START_REF] Morelos-Zaragoza | The Art of Error Correcting Coding[END_REF] is the probability P(v) that the intersymbol interference affecting the dot u k -i,l -j and observed from y k,l is equal to

j l i k j i l k u h x v - - - = , , ,
. The extrinsic information transmitted to decoder l k L , is then:

∑ ∑ - = - = → →             - = = - - - - n n i p p j l k u r l k u r l k u m u m L l k j l i k l k j l i k ) 1 ( ) 1 ( log , , , , , , , (7)
as there are (2n+1)x(2p+1) observed variables related to one hidden variable u k,l.

Complexity of the sum product algorithm grows exponentially with the memory size, which make its application practically unfeasible for some channels such for long memory channels. The objective of this work is primarily to propose a reduced complexity version of the sum product algorithm, aiming at proposing an algorithm with a linear complexity according to the memory size. In estimation theory, linearity is obtained when the joint density function of the hidden and the observed variables is Gaussian. Let us assume this hypothesis which could fairly be justified in case of long memory channel and yielding a Normal distribution for P(v) with mean and variance Ep/N0

Fig. 1 .

 1 Fig. 1. A 2D Alphacode illustration. (a) a word of five characters. (b) an isolated character with two columns and ten rows.

Fig. 2 .

 2 Fig. 2.A typical image illustrating the various artifacts encountered in barcode reading from a camera-equipped mobile phone.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Connected component based on image gradient intensity.

  image in Fig.4. b. We may expect from this transformation that dots are stretched according to one direction, giving rectangular-shaped dots.

Fig. 5 .

 5 Fig. 5. Result of grid structure recovery. The blue vertical lines denote character separation.

Fig. 6 .

 6 Fig. 6. Factor graph representation of the 2D ISI channel combined to the iterative equalizer making use of the parity code.

Fig. 7 .

 7 Fig. 7. A 2D Alphacode chain of character (60 data characters and 30 redundant characters) with a dot of 4x4 pixels and a blur matrix of 5x5. This image is corrupted with a noise with Ep/N0 = 0 dB, where Ep is the mean energy received by one dot and N0 the spectral density of the noise.

Fig. 8 .

 8 Fig. 8. A 2D Alphacode chain of character (60 data characters and 30 redundant characters) with a dot of 2x2 pixels and a blur matrix of 3x3. This image is corrupted with a noise with Ep/N0 = 9 dB, where Ep is the mean energy received by one dot and N0 the spectral density of the noise.

Fig. 9 .

 9 Fig.9. Performance in bit error rate of an Alphacode with the proposed algorithm (dashed lines) compared to a classical reader (solid line). The received image have a dot size of 4x4 pixels and a blur of size 5x5. Different RS encoders are used.

are the mean and the variance of the variable node ' , ' j l i k u --calculated according to the a priori information provided by the parity check decoder. This approximation enables us to conclude a simple expression for the extrinsic information in (7):

After several iterations this iterative structure converges and one makes the following decision:

which may be used by an error control decoder. We assume as in classical 2D barcodes, Reed Solomon (RS) encoder [START_REF] Morelos-Zaragoza | The Art of Error Correcting Coding[END_REF] protecting information from errors. The redundant symbols may be added at the end of the Alphacode chain of character. We consider an RS encoder constructed over Galois Field GF(256), and compare the data error probability obtained with the novel iterative equalizer and the classical binarisation based detector which equips almost all available readers in the market. We intentionally simulate poor quality images with high amount of blur. In these simulations the matrix modeling the blur is larger than a dot. Fig. 7 shows the bad quality of the received image with a dot of 4x4 pixels and a blur matrix size of 5x5, and Figure 8 represents an image with a dot of 2x2 pixels and a blur of 3x3. Fig. 9 and 10 show the performances obtained for the two latter cases. It is obvious that the proposed algorithm outperforms standard readers. These figures represent the data error probability with respect to the received signal to noise ratio expressed in terms of E p /N 0 in dB, where E p is the mean energy received by one dot and N 0 the spectral density of the noise.

V. CONCLUSION

In this work, we present a study of the 2D alphabetic barcode based on the Dotem alphabet in real situation, i.e. under difficult acquisition conditions (noise, blur, complex background and geometrical deformation). We have implemented a two-steps solution to decode a word from an image acquired using a hand-held camera. First, we defined some image processing tasks, including code localization in a complex scene and projective-based image rectification to recover the word structure. Then, we studied the decoding process using a modified sum-product algorithm running on a factor graph modeling the displaying/reading process. We make use of the parity dot in each column to derive an iterative structure of this equalizer. This algorithm has a linear complexity which makes it a realistic and attractive reading solution for camera-phones. Reed-Solomon encoders are then used to control errors. The resulting data error rate is greatly improved compared to basic processing techniques equipping most of the barcode readers for the identification and interacting markets.


Dawen Wu
email: dawen.wu@centralesupelec.fr

Abdel Lisser

Neuro-PINN: A Hybrid Framework for Efficient Nonlinear Projection Equation Solutions

Keywords: Nonlinear projection equation, Physics-informed neural network, Neurodynamic optimization, Ordinary differential equation

Nonlinear projection equations (NPEs) provide a unified framework for solving various constrained nonlinear optimization and engineering problems. This paper presents a deep learning approach for solving NPEs by incorporating neurodynamic optimization and physics-informed neural networks (PINNs). First, we model the NPE as a system of ordinary differential equations (ODEs) using neurodynamic optimization, and the objective becomes solving this ODE system. Second, we use a modified PINN to serve as the solution for the ODE system. Third, the neural network is trained using a dedicated algorithm to optimize both the ODE system and the NPE.

Unlike conventional numerical integration methods, the proposed approach predicts the end state without computing all the intermediate states, resulting in a more efficient solution. The effectiveness of the proposed framework is demonstrated on a collection of classical problems, such as variational inequalities and complementarity problems.

INTRODUCTION

Nonlinear optimization problems are prevalent across a wide range of fields, such as engineering, physics, and economics. Solving these problems requires finding an optimal solution that satisfies a defined set of constraints while optimizing the objective function [START_REF] Bertsekas | Nonlinear programming[END_REF] . Nonlinear projection equations (NPEs) are a useful tool for formulating these problems and serve as a unifying framework for treating various nonlinear optimization problems, including nonlinear complementarity problems, variational inequalities, and equilibrium point problems [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] .

The NPEs are typically addressed by 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] Xia | A Bi-Projection Neural Network for Solving Constrained Quadratic Optimization Problems[END_REF][START_REF] Che | A collaborative neurodynamic approach to global and combinatorial optimization[END_REF] . The constructed ODE system must be shown to have the global convergence property, i.e., the state solution of the system converges to the optimal solution of the NPE, regardless of the initial point. The NPE problem is then transformed into solving the state solution of the ODE system. However, the ODE system is typically highly nonlinear and has no analytical solutions. Therefore, numerical integration methods such as Runge-Kutta (RK) methods are commonly used to solve the state solution [START_REF] Burden | Numerical analysis[END_REF] .

Motivation. Despite the usefulness of numerical integration methods, they are not efficient enough for solving the NPE problems. This inefficiency is due to the fact that the state solution of the ODE system only provides a solution to the NPE at the end state. To reach the end state, numerical integration methods require the computation of all intermediate states starting from the initial state, making the process computationally intensive. Therefore, there is a need for a more efficient algorithm to solve the NPE problem.

Related Works

Neurodynamic optimization. Over the past few decades, numerous neurodynamic approaches have been developed to solve various constrained optimization problems, including linear and quadratic programming [START_REF] Tank | SIMPLE 'NEURAL' OPTIMIZATION NETWORKS: AN A/D CONVERTER, SIGNAL DE-CISION 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] Xia | A projection neural network and its application to constrained optimization problems[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] , non-smooth optimization [START_REF] Jia | An Adaptive Continuous-Time Algorithm for Nonsmooth Convex Resource Allocation Optimization[END_REF] , and 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] . In particular, a projection neurodynamic model for solving NPEs was introduced and found to have global convergence to the exact solution under mild conditions [START_REF] Xia | A new neural network for solving nonlinear projection equations[END_REF] . The model also showed both asymptotic and exponential stability without the need for a smooth nonlinear mapping. To further improve the performance, a bi-projection neurodynamic model was developed to efficiently solve quadratic optimization problems [START_REF] Xia | A Bi-Projection Neural Network for Solving Constrained Quadratic Optimization Problems[END_REF] . In addition, 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 avenue of research explored in this paper is the use of deep learning to solve differential equations. In the 1990s, Lagaris et al. used neural networks (NNs) to serve as solutions to both ordinary and partial differential equations (ODEs and PDEs), embedding boundary conditions directly into the network architecture [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] . The advent of deep learning has reinvigorated interest in using such methods to tackle high-dimensional, nonlinear PDEs [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 key contribution in this area 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 differential equations and data errors into the loss function. The versatile architecture and efficient training algorithms of PINNs have enabled numerous successes in a wide range of computational challenges in physics and engineering [START_REF] Samaniego | An energy approach to the solution of partial differential equations in computational mechanics via machine learning: Concepts, implementation and applications[END_REF][START_REF] Mao | Physics-informed neural networks for high-speed flows[END_REF][START_REF] Cai | Physics-informed neural networks (PINNs) for fluid mechanics: A review[END_REF][START_REF] Abueidda | Meshless physics-informed deep learning method for three-dimensional solid mechanics[END_REF][START_REF] Motlagh | Deep learning phase-field model for brittle fractures[END_REF] . This expanding research landscape is driven by a combination strategy that adapts PINNs to exploit the structural properties of the target problem. As a result, a variety of PINN variants have emerged to address different 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 improve computational efficiency [START_REF] Yu | Gradient-enhanced physics-informed neural networks for forward and inverse PDE problems[END_REF][START_REF] Sharma | Accelerated Training of Physics Informed Neural Networks (PINNs) using Meshless Discretizations[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][START_REF] Abueidda | Enhanced physics-informed neural networks for hyperelasticity[END_REF] . Apart from collocation-based PINN approaches, many studies use deep energy methods to solve PDEs, where the energy of the system is used as a loss function to train the NN model [START_REF] He | A deep learning energy-based method for classical elastoplasticity[END_REF][START_REF] Samaniego | An energy approach to the solution of partial differential equations in computational mechanics via machine learning: Concepts, implementation and applications[END_REF][START_REF] Yu | The deep Ritz method: a deep learning-based numerical algorithm for solving variational problems[END_REF] . Many packages have been developed to support the use of deep learning for solving PDEs [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] .

Key Contributions

This paper presents several key contributions:

• We propose a deep learning approach to NPEs that combines neurodynamic optimization with PINNs. Our approach reformulates an NPE problem as an NN training problem, thereby eliminating the need for numerical methods to solve NPEs.

• To improve the performance of our proposed approach, we design a specialized training algorithm that focuses on the end state of the ODE system. In each training iteration, the algorithm uses an evaluation metric based on the NPE error to assess the predictive accuracy of the NN model and retain the optimal result. Therefore, the NN model is not only trained to solve the ODE system, but also to solve the NPE problem.

• In our experiments, we demonstrate the effectiveness of the proposed approach in solving large-scale NPE problems. We compare our approach with the traditional RK method and PINN to validate its performance. In addition, we perform a hyperparameter sensitivity analysis to investigate the impact of different hyperparameter configurations on the performance of our proposed approach.

Paper Outline

The remaining sections are organized as follows: Section 2 provides the necessary background, including an introduction to the NPE problem, neurodynamic optimization, and the RK method. In Section 3, we present our NN model for the NPE problem. Section 4 details the design of the loss function and the training algorithm for the proposed NN model. Section 5 presents the experimental results, and we compare our approach with the RK method and PINN. Finally, Section 6 summarizes the main results of the paper and outlines possible directions for future research.

PRELIMINARIES

Section 2.1 introduces the NPEs that this paper aims to solve, and discusses their equivalence to other nonlinear optimization problems. Section 2.2 describes the neurodynamic approach to modeling NPEs as ODE systems. Section 2.3 introduces the RK method, a classical numerical integration method for solving the ODE system.

Nonlinear Projection Equation

Definition 1 (Nonlinear projection equation). Consider a nonlinear mapping 𝐺 ∶ ℝ 𝑛 → ℝ 𝑛 and a feasible set Ω ⊂ ℝ 𝑛 . The projection function 𝑃 Ω ∶ ℝ 𝑛 → Ω maps a point 𝑧 ∈ ℝ 𝑛 onto Ω, such that:

𝑃 Ω (𝑧) = arg min 𝑥∈Ω ‖𝑧 -𝑥‖ (1)
where ‖ ⋅ ‖ denotes the Euclidean norm.

The NPE problem, denoted by 𝑁𝑃 𝐸(Ω, 𝐺), is to find a vector 𝑥 * ∈ Ω satisfying:

𝑃 Ω (𝑥 * -𝐺(𝑥 * )) = 𝑥 * (2)
Definition 2 (Nonlinear complementarity problem). Consider a nonlinear mapping 𝐺 ∶ ℝ 𝑛 → ℝ 𝑛 . The nonlinear complementarity problem, denoted by 𝑁𝐶𝑃 (𝐺), is to find a vector 𝑥 * ∈ ℝ 𝑛 satisfying:

𝐺(𝑥 * ) ≥ 0, 𝑥 * ≥ 0, 𝐺(𝑥 * ) 𝑇 𝑥 * = 0. (3)
Definition 3 (Variational inequality). Consider a nonlinear mapping 𝐺 ∶ ℝ 𝑛 → ℝ 𝑛 and a feasible set Ω ⊂ ℝ 𝑛 . The variational inequality problem, denoted by 𝑉 𝐼(Ω, 𝐺), is to find a vector 𝑥 * ∈ Ω satisfying: According to the literature 2 , 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. The equilibrium point of a stochastic Nash game is expressed as a variational inequality [START_REF] Singh | Variational inequality formulation for the games with random payoffs[END_REF] , which is reformulated as an NPE by Proposition 2.

(𝑥 -𝑥 * ) 𝑇 𝐺 (𝑥 * ) ≥ 0, 𝑥 ∈ Ω. (4)

Neurodynamic Optimization

Assumption 1.

• The function 𝐺(⋅) is locally Lipschitz continuous.

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

Ω = {𝑥 = (𝑥 1 , 𝑥 2 , … , 𝑥 𝑛 ) ∈ ℝ 𝑛 | 𝑙 𝑖 ≤ 𝑥 𝑖 ≤ ℎ 𝑖 },
where 𝑙 𝑖 and ℎ 𝑖 denote the lower and upper bounds of 𝑥 𝑖 respectively. In this case, the projection function 𝑃 Ω (⋅) in equation ( 1) is reduced as follows:

𝑃

Ω (𝑥) = ( 𝑃 1 Ω (𝑥 1 ), 𝑃 2 Ω (𝑥 2 ), … , 𝑃 𝑛 Ω (𝑥 𝑛 ) ) 𝑇 , (5)
where 𝑃 𝑖 Ω (𝑥 𝑖 ), 𝑖 ∈ {1, 2, … , 𝑛} is defined as:

𝑃 𝑖 Ω (𝑥 𝑖 ) = ⎧ ⎪ ⎨ ⎪ ⎩ 𝑙 𝑖 if 𝑥 𝑖 < 𝑙 𝑖 , 𝑥 𝑖 if 𝑙 𝑖 ≤ 𝑥 𝑖 ≤ ℎ 𝑖 , ℎ 𝑖 if 𝑥 𝑖 > ℎ 𝑖 . ( 6
)
Let 𝑦 ∶ ℝ → ℝ 𝑛 be a time-dependent function, where 𝑦(𝑡) represents the state at time 𝑡. The aim of neurodynamic optimization is to design a first-order ODE system, 𝑑𝑦 𝑑𝑡 , to govern 𝑦(𝑡). In this paper, we employ the projection neurodynamic model introduced by [START_REF] Xia | A new neural network for solving nonlinear projection equations[END_REF] to model the NPE, where the ODE system is defined as follows:

𝑑𝑦 𝑑𝑡 = -𝐺 ( 𝑃 Ω (𝑦) ) + 𝑃 Ω (𝑦) -𝑦. ( 7
)
To simplify the discussion, we define:

Φ(𝑦) = -𝐺 ( 𝑃 Ω (𝑦) ) + 𝑃 Ω (𝑦) -𝑦. ( 8
)
Thus, the ODE system (7) can be written as

) 9
where 𝑥 * is a satisfied point of 𝑁𝑃 𝐸(Ω, 𝐺).

In particular, if 𝑁𝑃 𝐸(Ω, 𝐺) contains only one satisfied point 𝑥 * , then the ODE system is globally asymptotically stable at 𝑥 * .

Initial value problem (IVP) construction.

In practice, in order to use the neurodynamic approach to solve the NPE, we need to construct an IVP consisting of three components: 1) the ODE system of Eq. ( 7), 2) an initial condition 𝑦(𝑡 0 ) = 𝑦 0 , and 3) a time range 𝑡 ∈ [𝑡 0 , 𝑇 ]. 𝑦(𝑡) for 𝑡 ∈ [𝑡 0 , 𝑇 ] represents the state solution of this IVP over the time range [𝑡 0 , 𝑇 ], where the end state, 𝑦(𝑇 ), is considered to be the predicted solution to the NPE. According to Theorem 1, the larger the time range [𝑡 0 , 𝑇 ], the closer 𝑦(𝑇 ) is to the optimal solution 𝑥 * of the NPE.

Runge-Kutta Method

Runge-Kutta (RK) methods are widely used numerical techniques for approximating the state solution of an ODE system. They are particularly useful when exact analytical solutions are difficult or infeasible to obtain. The methodology involves approximating the state solution over a series of discrete time points. In this paper, we concentrate on the fourth-order RK method, commonly known as the RK method, which is one of the most accurate and popular methods in the family of RK methods.

The RK method takes as input an ODE system 𝑑𝑦 𝑑𝑡 = Φ(𝑦), an initial condition 𝑦(𝑡 0 ) = 𝑦 0 , and a time range [𝑡 0 , 𝑇 ]. The method sets 𝑁 collocation points that are equally spaced on the time range [𝑡 0 , 𝑇 ]. The RK method returns an approximate state solution, denoted by ȳ(𝑡), where 𝑡 is a time point within the finite set {𝑡 0 , 𝑡 1 , … , 𝑡 𝑁 , 𝑡 𝑇 }. The RK algorithm is stated as follows:

• Step 1: Initialize the step size ℎ = (𝑇 -𝑡 0 ) 𝑁 + 1 , the initial time 𝑡 = 𝑡 0 , and set the initial state ȳ(𝑡 0 ) = 𝑦 0 .

• Step 2: For 𝑖 = 1, 2, … , 𝑁 + 1 do Steps 3 and 4.

• Step 3: Set

𝐾 1 = ℎΦ ( ȳ(𝑡)) , 𝐾 2 = ℎΦ ( ȳ(𝑡) + 𝐾 1 2
) , 𝐾 3 = ℎΦ ( ȳ(𝑡) + 𝐾 2 2
) ,

𝐾 4 = ℎΦ ( ȳ(𝑡) + 𝐾 3 ) . ( 10
)
• Step 4: Set

ȳ(𝑡 + ℎ) = ȳ(𝑡) + 𝐾 1 + 2𝐾 2 + 2𝐾 3 + 𝐾 4 6 , (11)
and 𝑡 = 𝑡 + ℎ. When applying the RK method to the NPE problem, the approximate end state ȳ(𝑇 ) is considered as the predicted solution of the NPE problem. Specifically, we have:

NN as

ȳ(𝑇 ) ≈ 𝑦(𝑇 ) ≈ 𝑥 * , ( 12
)
where ȳ(𝑇 ) ≈ 𝑦(𝑇 ) indicates that the end state obtained by the RK method is an approximation of the true end state, and 𝑦(𝑇 ) ≈ 𝑥 * indicates that the true end state is the predicted solution of the NPE problem according to Theorem 1.

NEURAL NETWORK MODEL

Model description. We propose a neural network (NN) model to solve the NPE problem. Our model can be expressed by the following equation:

ŷ(𝑡; 𝑤) = 𝑦 0 + ( 1 -𝑒 -(𝑡-𝑡 0 ) ) N(𝑡; 𝑤), 𝑡 ∈ [𝑡 0 , 𝑇 ], (13)
where N(𝑡; 𝑤) represents a fully connected NN with trainable parameters 𝑤, and [𝑡 0 , 𝑇 ] is a given time range. We use the Lagaris method to incorporate the initial conditions of the ODE system into the NN model [START_REF] Lagaris | Artificial neural networks for solving ordinary and partial differential equations[END_REF] . In particular, the auxiliary function

( 1 -𝑒 -(𝑡-𝑡 0 ) )
ensures that the NN model always satisfies the initial condition (𝑡 0 , 𝑦 0 ), i.e., ŷ(𝑡 = 𝑡 0 ; 𝑤) = 𝑦 0 , regardless of the model parameters 𝑤. The exponential form in the auxiliary function has been demonstrated to improve the convergence of the model [START_REF] Mattheakis | Hamiltonian neural networks for solving equations of motion[END_REF] .

Approximate state solution of the ODE. As shown in Fig. 1 (Left), the proposed model itself is an approximate state solution to the ODE system of Eq. ( 7) on the time range [𝑡 0 , 𝑇 ], i.e.,

ŷ(𝑡; 𝑤) ≈ 𝑦(𝑡), 𝑡 ∈ [𝑡 0 , 𝑇 ], (14)
where where ŷ(𝑡 = 𝑇 ; 𝑤) ≈ 𝑦(𝑇 ) indicates that the end state of our model approximates the true end state, and 𝑦(𝑇 ) ≈ 𝑦 * comes from Theorem 1, indicating that 𝑦(𝑇 ) solves the NPE. Unlike most PINN models, which aim to solve for the entire input space [𝑡 0 , 𝑇 ], our NN model focuses on the end state of the ODE system, since it represents the predicted solution to the NPE problem. In the following section, we will show how to train the NN model with an emphasis on improving the prediction accuracy of the end state.

MODEL TRAINING

Training Objective

Loss function. The loss function of the proposed NN model is defined as:

(𝑡, 𝑤) = ‖ ‖ ‖ ‖ 𝜕 ŷ(𝑡; 𝑤) 𝜕𝑡 -Φ( ŷ(𝑡; 𝑤)) ‖ ‖ ‖ ‖ , ( 16
)
where Φ(⋅) refers to the ODE system, corresponding to the NPE problem to be solved. Φ( ŷ(𝑡; 𝑤)) represents the expected derivative according to the ODE system, and 𝜕 ŷ(𝑡;𝑤) 𝜕𝑡 represents the actual derivative of the NN model. 𝜕 ŷ(𝑡;𝑤) 𝜕𝑡 can be computed using automatic differentiation tools such as PyTorch or JAX [START_REF] Paszke | PyTorch: An Imperative Style, High-Performance Deep Learning Library[END_REF][START_REF] Bradbury | JAX: composable transformations of Python+NumPy programs[END_REF] . (𝑡, 𝑤) represents the difference between the two at time 𝑡 and with network parameters 𝑤. As shown in Fig. 2, the NPE is first reformulated as an ODE system via neurodynamic optimization. The ODE system is then incorporated into the loss computational process.

Objective function. The objective function of the NN model is defined as: Batch loss. However, the objective function 𝐽 (𝑤) is computationally intractable due to its integral part. Therefore, in practice, we cannot directly use 𝐽 (𝑤) to train the NN model. Instead, we train the model by minimizing the following batch loss:

𝐽 (𝑤) = 𝑇 ∫ 𝑡 0 (𝑡, 𝑤)𝑑𝑡, ( 17
(𝕋 , 𝑤) = 1 |𝕋 | ∑ 𝑡∈𝕋 (𝑡, 𝑤), ( 18
)
where 𝕋 is a set of time points uniformly sampled from the interval [𝑡 0 , 𝑇 ], and |𝕋 | denotes the size of the set. The batch loss (𝕋 , 𝑤) approximates the objective function 𝐽 (𝑤) by a sum of loss values over a set of sampled time points. By minimizing the batch loss, we can effectively train the model to solve the NPE.

Algorithm Design

NPE error. We introduce an evaluation metric, called NPE error, to measure how well a prediction 𝑥 pred solves the NPE problem. The metric is defined as:

NE(𝑥 pred ) = ‖ ‖ ‖ 𝑃 Ω (𝑥 pred -𝐺(𝑥 pred )) -𝑥 pred ‖ ‖ ‖∞ , ( 19
)
where ‖⋅‖ ∞ represents the infinity norm.

Algorithm 1 Deep learning solver for NPE based on neurodynamic optimization

Input: 𝑁𝑃 𝐸(Ω, 𝐺) as defined in Eq. ( 2); Time range [𝑡 0 , 𝑇 ]; Initial condition (𝑡 0 , 𝑦 0 ). Output: 𝑥 best , the NN prediction to 𝑁𝑃 𝐸(Ω, 𝐺).

1: function SOLVER(𝑁 𝑃 𝐸(Ω, 𝐺), [𝑡 0 , 𝑇 ], 𝑦 0 ) 2:
Derive the ODE system, Φ(⋅), according to Eq. ( 7).

3:

Initialize a NN model ŷ(𝑡; 𝑤). return 𝑥 best 18: end function Pipeline. Alg. 1 summarizes how to use our proposed approach to solve an NPE problem. First, we need to specify an initial condition (𝑡 0 , 𝑦 0 ) and a time range [𝑡 0 , 𝑇 ] to construct an IVP for the NPE problem. Then, we initialize the proposed NN model (13), which serves as an approximate state solution for this IVP. The model is trained by performing gradient descent on the batch loss of Eq. ( 18) to improve the approximation. Note that our solver is completely based on the deep learning infrastructure and does not require any standard optimization solver or numerical integration solver.

Optimal result retention (ORR) mechanism. A key to Alg. 1 is that we use an ORR mechanism based on the evaluation metric of Eq. (19). Specifically, in each iteration, the algorithm compares the NPE error of the current iteration, denoted as NE curr , with the lowest NPE error found so far, denoted as NE best . Correspondingly, 𝑥 curr and 𝑥 best represent the current prediction and the best predictions found so far, respectively. If NE curr is less than NE best , it means that the model found a better prediction in the current iteration. The algorithm then updates NE best to equal NE curr and stores the best prediction as 𝑥 best = 𝑥 curr . This mechanism ensures that the best prediction obtained by the model is maintained throughout the training process, improving the overall performance of the algorithm.

Comparison with the RK Method

Fig. 3 compares the solution procedures between our proposed NN approach and the RK method, both of which solve the NPE problem by solving the IVP. Both approaches use the end state as the predicted solution for the NPE problem. However, they differ in how they enhance the accuracy of the end state.

NN solution procedure. The NN approach employs gradient descent on the batch loss of Eq. ( 18) to improve the NN prediction at each iteration. The evolution of the NN model is represented by ŷ(𝑡; 𝑤 1 ), ŷ(𝑡; 𝑤 2 ), … , ŷ(𝑡; 𝑤 𝑀 ), where 𝑤 𝑖 and ŷ(𝑡; 𝑤 𝑖 ) denote the model parameters and the approximate state solution at the 𝑖-th iteration, respectively. The predicted end states are ŷ(𝑡 = 𝑇 ; 𝑤 1 ), ŷ(𝑡 = 𝑇 ; 𝑤 2 ), … , ŷ(𝑡 = 𝑇 ; 𝑤 𝑀 ), where ŷ(𝑡 = 𝑇 ; 𝑤 𝑖 ) represents the NPE prediction at the 𝑖-th iteration.

RK solution procedure. In contrast, the RK method computes discrete collocation points iteratively. The method progresses by solving Eq. ( 10) and Eq. ( 11) to obtain ȳ𝑗 and 𝑡 𝑗 , which represent the solved state values and collocation point at the 𝑗-th iteration, respectively. The state value ȳ𝑗 incorporates all previously solved state values. At the end of the 𝑗-th iteration, ȳ𝑗 is used as the prediction for the NPE.

EXPERIMENTS

Section 5.1 delineates the application of the proposed NN approach for solving various types of NPE problems. Section 5.2 contrasts our approach with the PINN. Section 5.3 investigates the performance of our NN approach over different network architectures and hyperparameter configurations. Section 5.4 demonstrates the effectiveness of our NN approach in solving large-scale NPE problems. Finally, Section 5.5 discusses the distinctive features and limitations of our proposed NN approach, while also outlining possible avenues for future research.

Three Examples

Experimental setup of our NN approach. We used PyTorch 1.12.1 [START_REF] Paszke | PyTorch: An Imperative Style, High-Performance Deep Learning Library[END_REF] to implement the proposed NN model and JAX 0.4.1 [START_REF] Bradbury | JAX: composable transformations of Python+NumPy programs[END_REF] to implement the ODE system. The NN model consists of a single fully connected layer with 100 neurons, and the activation function is Tanh. For training, we used the Adam optimizer with a learning rate of 0.001, a batch size of 128.

Experimental setup of the RK method. We used the RK method [START_REF] Dormand | A family of embedded Runge-Kutta formulae[END_REF] for comparison and called it via the Scipy library [START_REF] Virtanen | SciPy 1.0: fundamental algorithms for scientific computing in Python[END_REF] . We set the number of collocation points to 50,000, evenly distributed over the time range.

Linear Complementarity Problem

Example 1: Consider the following linear complementarity problem:

𝑥 𝑇 (𝑀𝑥 + 𝑞) = 0, 𝑥 ≥ 𝟎, 𝑀𝑥 + 𝑞 ≥ 𝟎, (20)
where

𝑀 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣
50 -8 -6 -9 12 -8 33 -1 -25 3 -6 -1 38 10 -4 -9 -25 10 55 -24 12 3 -4 -24 20

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , 𝑞 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ -2 -20 -16 -12 -14 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 21
)
The goal is to find an optimal solution 𝑥 * ∈ ℝ 5 that solves Eq. ( 20). Example 1 is reformulated as 𝑁𝑃 𝐸(ℝ + , 𝐺) by Proposition 1, where 𝐺(𝑥) = 𝑀𝑥 + 𝑞. Then, 𝑁𝑃 𝐸(ℝ + , 𝐺) is modeled by the ODE system of Eq. (7). We establish an IVP by specifying the initial condition as 𝑦(0) = 𝟎 and the time range as [0, 10]. We use the NN model proposed in Eq. ( 13), denoted as ŷ(𝑡; 𝑤), to serve as an approximate state solution for this IVP. The end state, ŷ(𝑡 = 10; 𝑤), serves as the predicted solution for both 𝑁𝑃 𝐸(ℝ + , 𝐺) and Example 1. The NN model is trained using Alg. 1 to improve accuracy. ) shows the results of the RK method after accumulating 0, 5,000, 25,000, and 50,000 collocation points, and the predicted solutions are marked with yellow stars.

Table 1 shows the predictions of our NN approach and the RK method for Example 1 at different iterations. The results suggest that our approach is comparably accurate to the RK method given the same initial condition and time range. After 10,000 iterations, the NPE error of our prediction is reduced to less than 0.1. The final solution from our NN approach is [0.07, 2.75, 0.08, 3.32, 4.24], with an NPE error of 0.03.

𝑥 𝑇 𝐹 (𝑥) = 0, 𝑥 ≥ 𝟎, 𝐹 (𝑥) ≥ 𝟎, (22)
where

𝐹 (𝑥) = ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ 𝑥 1 𝑒 (𝑥 2 1 +(𝑥 2 -1) 2 ) + 𝑥 2 2 + 𝑥 3 -10 (𝑥 2 -1)𝑒 𝑥 2 1 +(𝑥 2 -1) 2 + 4𝑥 1 + 𝑥 2 𝑥 3 + 2𝑥 2 3 + 𝑒 𝑥 4 -2 𝑥 1 + 8𝑥 2 + 3𝑥 3 -3 𝑥 4 -4 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ . ( 23
)
We reformulate Example 2 as 𝑁𝑃 𝐸(ℝ + , 𝐹 ) by Proposition 1. Then, 𝑁𝑃 𝐸(ℝ + , 𝐹 ) is modeled by the ODE system of Eq. ( 7). We establish an IVP by specifying the initial point as 𝑦(0) = 𝟎 and the time range as [0, 10]. We use the NN model, ŷ(𝑡; 𝑤), to serve as an approximate state solution for this IVP, where the end state ŷ(𝑡 = 10; 𝑤) is the predicted solution for Example 2. Fig. 5(A) and Fig. 5(B) illustrate the decrease in loss and NPE error, respectively. Specifically, the loss decreases from 17.78 to 0.05, while the NPE error decreases from 5.04 to a value close to zero. Notably, the most substantial reduction in both loss and NPE error occurs within the first 10,000 iterations. In addition, the implementation of the ORR mechanism, as described in Alg. 1, ensures that the NPE error decreases consistently, despite occasional small increases in the loss values.

Fig. 5(C) shows the NN model ŷ(𝑡; 𝑤) at the 0th, 1,000th, 10,000th, and 50,000th training iterations, from left to right, and the predicted solution for Example 2 are marked with red stars. As shown in the figure, the NN model at the 1,000th iteration is already very close to the final result at the 50,000th iteration. Fig. 5(D) shows the results of the RK method after accumulating 0, 5,000, 25,000, and 50,000 collocation points, and the predicted solutions are marked with yellow stars.

Table 2 shows the predictions and their NPE errors of our NN approach and the RK method for Example 2 at different iterations. Thanks to the adoption of the projection function in Alg. 1, our NN approach has a lower initial NPE error compared to the RK method. By the 1,000th iteration, the NPE error associated with our NN approach has been reduced to 0.14. By the 5,000th iteration, the NPE error has further diminished to 0.01. Upon completion of 10,000 iterations, our NN approach yields an optimal solution of [1.08, 0.00, 0.64, 4.00], which accurately resolves Example 2.

Variational Inequality

Example 3: Consider the following variational inequality:

(𝑥 -𝑥 * ) 𝑇 𝐺 (𝑥 * ) ≥ 0, 𝑥 ∈ Ω, (24)
where

𝐺(𝑥) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 𝑥 1 -2 𝑥 1 +0.8 + 5𝑥 2 -13 1.2𝑥 1 + 7𝑥 2 3𝑥 3 + 8𝑥 4 1𝑥 3 + 2𝑥 4 -4 𝑥 4 +2 -12 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ , Ω = {𝑥 ∈ ℝ 4 |1 ≤ 𝑥 1 ≤ 100, -3 ≤ 𝑥 2 ≤ 100, -10 ≤ 𝑥 3 ≤ 100, 1 ≤ 𝑥 4 ≤ 100}. ( 25
)
We reformulate Example 3 as 𝑁𝑃 𝐸(Ω, 𝐺) by Proposition 2. Then, 𝑁𝑃 𝐸(Ω, 𝐺) is modeled by the ODE system of Eq. ( 7). We establish an IVP by specifying initial point as 𝑦(0) = 𝟎 and the time interval as [0, 10]. We use the NN model, ŷ(𝑡; 𝑤), to serve as an approximate state solution for this IVP, where the end state ŷ(𝑡 = 10; 𝑤) is the predicted solution for Example 3.

Fig. 6(A) and Fig. 6(B) illustrate the decrease in loss and NPE error, respectively. Specifically, the loss decreases from 126.99 to 0.10, while the NPE error decreases from 15.42 to 0.00. Note that there are small fluctuations in the loss value around the 18,000th iteration, but the NPE error remains unaffected because the training algorithm retains the best prediction from its training history.

Fig. 6(C) shows the NN model ŷ(𝑡; 𝑤) at the 0th, 1,000th, 10,000th, and 50,000th training iterations, from left to right, and the predicted solution for Example 3 are marked with red stars. In particular, at the 1,000th iteration, i.e., the first subplot on the left in Fig. 6(C), the NN model is already very close to the final result of the 50,000th iteration. Fig. 6(D) shows the results of the RK method after accumulating 0, 5,000, 25,000, and 50,000 collocation points, and the predicted solutions are marked with yellow stars. 3 shows the predictions for Example 3 provided by our NN approach and the RK method. Remarkably, after only 1,000 iterations, our NN approach yields an acceptable prediction with an NPE error of 0.46. After 5,000 iterations, our NN approach refines this prediction further to an NPE error of less than 0.1. After 10,000 iterations, our approach converges to the solution of [28.07, -3.00, -10.00, 11.15], which is an optimal solution of Example 3.

Comparision with PINN

In this subsection, we compare our approach with PINNs in terms of accuracy. Specifically, we have chosen the basic version of PINN [START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF] for comparison. Despite the existence of more advanced models such as CPINNs [START_REF] Jagtap | Conservative physics-informed neural networks on discrete domains for conservation laws: Applications to forward and inverse problems[END_REF] or XPINNs [START_REF] Jagtap | Extended Physics-informed Neural Networks (XPINNs): A Generalized Space-Time Domain Decomposition based Deep Learning Framework for Nonlinear Partial Differential Equations[END_REF] , however, we find that the basic version of PINN is sufficiently effective for the NPE problem, obviating the need for more complex variants. Regarding the setting for the PINN, we use the same network architecture and hyperparameters as in Section 5.1.

NPE error

Our PINN

Example 1 Example 2 Example 3

Figure 7 Comparison of NPE errors between our proposed approach and the PINN for the three examples given in Section 5.1.

Our approach can be viewed as a modification of the PINN, specifically designed for the NPE problem to improve computational performance. The core distinction between our approach and the PINN lies in our focus on the end state of the NN model. Building on this, we employ the ORR mechanism in Alg. 1, which continuously monitors the NPE error of the end state throughout the training process. Therefore, the NN model optimizes simultaneously for both the ODE system and the NPE problem at hand, consistently maintaining the best result throughout the solution process. This unique focus results in improved performance, as demonstrated below. Fig. 7 compares our approach with the PINN for solving the three examples given in Section 5.1. As shown in the figure, both approaches start with similar initial errors, but our approach is significantly outperforms the PINN as training progresses. Specifically, in Example 1, the NPE error converges to 0.03 with our approach, while the PINN converges to 0.7. In Example 2, our approach converges to less than 0.001, while the PINN could only converge to 0.4. In Example 3, our approach again converges to less than 0.001, while the PINN only converges to 0.1. These experimental results show that our approach provides superior solutions for solving these NPE problems. Moreover, this validates the effectiveness of the key designs in Alg. 1, such as the use of the ORR mechanism and the projection function. In this subsection, we explore the influence of various hyperparameters on the computational performance of our NN model. Specifically, we focus on Example 3, as given in Section 5.1.3. The hyperparameters under investigation include the initial point, the time range, the number of hidden layers, the number of neurons, and the activation function used. Initial point. Fig. 8(A) and Table 4 show the results of different initial point configurations with the same time range of [0, 10]. The results suggest the following:

Hyperparameter Study

• All initial points converge to same optimal solution, as supported by Theorem 1.

• The convergence rates of different initial points vary, with initial points closer to the optimal solution converging faster.

• The initial point 𝑦(0) = [1, 2, 3, 4] is closest to the optimal solution and achieves the fastest convergence with the smallest initial NPE error.

• The initial point 𝑦(0) = [-10, -15, -10, -14] is the furthest away from the optimal solution and still converges, but has a higher initial error and slower convergence rate.

Time range. Fig. 8(B) and Table 5 show the results of different time range configurations with the same initial points of 𝑦(0) = [0, 0, 0, 0]. The results suggest the following:

• Shorter time ranges lead to faster convergence but may result in less accurate predictions. As shown in the table, the shortest time range of [0, 5] converges very fast, but its NPE error does not decrease much after 3,000 iterations.

• Longer time ranges provide better predictions, but require more training iterations. The longest range of [0, 15] converges slowly, but with more training it can give better results than the other two ranges.

• The choice of time represents trade-off. ranges may enhance accuracy but require more training, whereas shorter ranges are easier to train but may yield less satisfactory predictions. As shown in Table 5, considering a fixed maximum number of iterations at 10,000, the time range of [0, 8] achieves the optimal performance. 6 shows the NPE errors of different model sizes, specifically different numbers of hidden layers and neurons. As shown in the table, all models of different sizes converge to solutions with NPE errors below 0.03 by the 10,000th iteration. However, the model size has a significant impact on the speed of convergence. Models that are too small or too large perform worse than the others. For example, the 100-neuron single-layer model has higher NPE errors in the first 1000 iterations. The 5-layer model with 500 neurons per layer has similar results.

Therefore, selecting an appropriately-sized model is importance to achieve optimal performance. Small model sizes have limited capacity, while large model sizes require a lot of training to converge, which may not be necessary. Among the model sizes considered in Table 6, the two-layer model with 500 neurons each achieves the best performance. It converges to an NPE error of 0.1 within 300 iterations and further reduces to 0.01 within 1000 iterations, outperforming other model architectures. 7 shows the NPE errors for different activation functions on the same NN model. The data show that regardless of the activation function chosen, all NN models converge to solutions with NPE errors less than 0.03. Notably, the influence of the activation function on model performance is relatively minor when compared to the effects of model size. Among the tested activation functions, ReLU and Leaky ReLU were found to be slightly less effective than the others. Specifically, the tanh activation function stands out as the most efficient, which corroborates its widespread adoption in research related to PINNs [START_REF] Jagtap | Adaptive activation functions accelerate convergence in deep and physicsinformed neural networks[END_REF][START_REF] De Ryck | On the approximation of functions by tanh neural networks[END_REF][START_REF] Mishra | Estimates on the generalization error of physics-informed neural networks for approximating a class of inverse problems for PDEs[END_REF] .

Large Scale NPE

Consider the following NPE problem: For 𝑖 = 1, 2, … , 1000,

( 𝑥 * 𝑖 - 1 2 √ (𝑀𝑥 * ) 𝑖 + 𝑞 𝑖 ) + = 𝑥 * 𝑖 . ( 26
)
The objective is to find an optimal solution 𝑥 * = [𝑥 * 1 , 𝑥 * 2 , … , 𝑥 * 1000 ] ∈ ℝ 1000 that solves Eq. ( 26). The problem data 𝑀 ∈ ℝ 1000×1000 is partitioned as

𝑀 = [ 𝑀 1 𝑀 2 𝑀 3 𝑀 4 ] , ( 27
)
where 𝑀 1 ∈ ℝ 500×500 , 𝑀 2 ∈ ℝ 500×500 , 𝑀 3 ∈ ℝ 500×500 , and 𝑀 4 ∈ ℝ 500×500 are given by

𝑀 1 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 1.2 0.6 … 0.6 0.6 1.2 … 0.6 ⋮ ⋮ ⋱ ⋮ 0.6 0.6 … 1.2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ , 𝑀 2 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 1 0 … 0 0 1 … 0 ⋮ ⋮ ⋱ ⋮ 0 0 … 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ , 𝑀 3 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ -1 0 … 0 0 -1 … 0 ⋮ ⋮ ⋱ ⋮ 0 0 … -1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ , 𝑀 4 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 … 0 0 0 … 0 ⋮ ⋮ ⋱ ⋮ 0 0 … 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 28
)
For the problem data 𝑞 ∈ ℝ 1000 , the first half is sampled uniformly from the interval [-31, -28], while the second half is sampled uniformly from the interval [3, 10]. We create a problem set by generating ten different 𝑞, resulting in ten different NPE problems. All the sampled 𝑞 can be accessed through the link [START_REF] Bertsekas | Nonlinear programming[END_REF] . Solving these NPE problems is non-trivial due to their high dimensionality with 1,000 variables and the presence of multiple nonlinear operations. Therefore, the RK method may encounter computational inefficiencies when applied to these NPE problems. We use the proposed NN approach to solve the ten large-scale NPE problems. With respect to the experimental setup, the employed network architecture consisted of a three-layer fully connected network, each layer having 500 neurons and utilizing the tanh activation function. The time range chosen for the experiments is [0, 10], and the initial point is set to a zero vector. All 8 and9 present the experimental results of our NN approach compared to the RK method, focusing on both accuracy and computational time. In order to offer a comprehensive understanding of the experimental results, we provide statistical descriptors of the outcomes, including the mean, standard deviation (STD), and various percentiles.

In the following, we discuss the differences between our NN approach and the RK method in terms of computational efficiency, stability and convergence.

• Efficiency: As shown in the tables, our NN approach outperforms the RK method in terms of computational efficiency when solving large scale NPE problems. In particular, our NN approach requires less CPU time than the RK method for the same level of accuracy. For example, our approach achieves an NPE error of 0.01, requiring an average of 30,000 iterations and consuming only 670.38 seconds. On the contrary, the RK method requires spanning a time range of [0, 10], which requires a higher average time of 1542.07 seconds. Moreover, the efficiency gap becomes more significant as the required accuracy is relaxed. For example, for an NPE error threshold of less than 0.50, our NN approach takes on average only 500 iterations and 11.32 seconds, making it about 48 times faster than the RK method, which requires 547.93 seconds.

• Stability: Our NN approach demonstrates not only computational efficiency, but also greater stability, reflected in a much lower STD of CPU times -only 5.82 at the 30,000th iteration. This indicates consistent and reliable performance across different NPE problem instances. On the other hand, the RK method exhibits higher variability with an STD of 205.43 for the time range of [0, 10], signifying greater sensitivity to the specific NPE problem at hand.

• Convergence: While there are notable differences in computational efficiency and stability between our NN approach and the RK method, both techniques demonstrate comparable accuracy in the long run, each achieving an NPE error as low as 0.01. As discussed in Section 4.3, the convergence mechanisms of the two approaches are fundamentally different. Our NN approach operates within a fixed time range of [0, 10] and refines its predictive accuracy through iterative training. In contrast, the RK method improves its accuracy by progressively extending the time range.

Discussion

Below, we summarize the key features of our proposed NN approach:

• Our NN approach reliably converges to the optimal solutions of the NPE problems. This is supported by Theorem 1 in neurodynamic optimization and the universal approximation theorem of NNs. Experimentally, we have shown that the NN approach successfully found optimal solutions for the three types of NPE problems in Section 5.1, as well as ten large-scale NPE problems in Section 5.4.

• The proposed NN approach outperforms the PINN in solving NPE problems. This improvement is due to some modifications made to the basic PINN approach that allow it to exploit the problem structure of NPEs for better performance.

In particular, during each training iteration, the NN model evaluates the accuracy of its end state against the target NPE problem and retains the best performing solution. Additionally, we employ the projection function of Eq. ( 1) to further boost accuracy.

• The computational performance of our NN approach is significantly affected by the hyperparameter settings. As discussed in Section 5.3, choices regarding the time range and initial point have a strong impact on the convergence result. For example, smaller time ranges may prevent the NN model from converging to an optimal solution, regardless of the training time. Moreover, in our empirical observations, the network architecture and activation functions influence the model convergence rate.

• Our NN approach excels in solving large-scale NPE problems. As presented in Section 5.4, our NN approach outperforms the RK method in terms of computational efficiency. For the same level of accuracy, the CPU time required by our approach is less than that required by the RK method. In addition, our approach exhibits greater stability, with its solution time remaining consistent across different NPE problems.

However, we must acknowledge some limitations of our proposed NN approach. The most notable limitation is that, unlike traditional RK methods, the NN approach lacks rigorous theoretical underpinnings to guarantee convergence, primarily due to the black-box nature of NNs. Furthermore, as elaborated in Section 5.3, the performance of the model is highly sensitive to hyperparameter choices and architectural decisions, requiring careful tuning. Numerical integration methods, which are typically simpler, avoid these complexities.

To address these limitations, future research should focus on improving network design and training methods. This could include developing more effective hyperparameter search strategies or exploring advanced network architectures, such as attention-based models or transformers. In addition, investigating how to incorporate domain-specific knowledge into the network structure may improve solution accuracy and reduce training time.

CONCLUSION

In this paper, we presented an innovative deep learning-based approach for solving NPEs based on neurodynamic optimization and PINNs. We showed how our approach can efficiently solve NPEs and highlighted its advantages over PINNs and the RK method. The proposed approach transforms NPE problems into NN training problems, allowing the use of the latest advances in machine learning and deep learning to solve NPEs. We also identified areas for future research, such as exploring better methods for selecting initial points and time ranges, experimenting with different network architectures, and investigating advanced neurodynamic optimization techniques, among others.

In summary, our proposed framework shows considerable promise for improving computational efficiency in solving NPEs. With ongoing research and development, we expect to further strengthen the robustness of our approach and position it as a valuable computational tool for addressing a wide range of nonlinear optimization challenges in diverse applications.

Figure 2

 2 Figure 2 Integrating an NPE into loss computation for NN training.

Section 4 .

 4 1 provides a definition of the loss function and objective function for the proposed NN model. Section 4.2 presents a training algorithm for the NN model. Section 4.3 compares our proposed NN approach with the RK method.

  ) which is the integral of the loss function over the time range [𝑡 0 , 𝑇 ]. The loss value (𝑡, 𝑤) represents the error of the model at time 𝑡, while the objective function 𝐽 (𝑤) represents the total error of the model over the time range [𝑡 0 , 𝑇 ].

Figure 3

 3 Figure 3Comparison of the solution procedures between the NN approach and RK method.

Figure 4

 4 Figure 4 Solving Example 1: (A) Mean square error (MSE) loss versus training iterations. (B) NPE error versus training iterations, where the NPE error is defined in Eq. (2). (C) Evolution of the NN solution, ŷ(𝑡; 𝑤). (D) Evolution of the RK solution, ȳ(𝑡).

Fig. 4 (

 4 Fig.4(A) and Fig.4(B)show the decreasing loss and NPE error, respectively, where the loss drops from 152.59 to 0.23 and the NPE error drops from 24.69 to 0.03. Fig.4(C) shows the NN model ŷ(𝑡; 𝑤) at the 0th, 1,000th, 10,000th, and 50,000th training iterations, from left to right, and the predicted solutions for Example 1 are marked with red stars. Fig.4(D)shows the results of the RK method after accumulating 0, 5,000, 25,000, and 50,000 collocation points, and the predicted solutions are marked with yellow stars.Table1shows the predictions of our NN approach and the RK method for Example 1 at different iterations. The results suggest that our approach is comparably accurate to the RK method given the same initial condition and time range. After 10,000 iterations, the NPE error of our prediction is reduced to less than 0.1. The final solution from our NN approach is [0.07, 2.75, 0.08, 3.32, 4.24], with an NPE error of 0.03.

Figure 5

 5 Figure 5 Solving Example 2: (A) MSE loss versus training iterations. (B) NPE error versus training iterations. (C) Evolution of the NN solution, ŷ(𝑡; 𝑤). (D) Evolution of the RK solution, ȳ(𝑡).

Figure 6

 6 Figure 6 Solving Example 3: (A) MSE loss versus training iterations. (B) NPE error versus training iterations. (C) Evolution of the NN solution, ŷ(𝑡; 𝑤). (D) Evolution of the the RK solution, ȳ(𝑡).

Figure 8 (

 8 Figure 8 (A) NPE errors versus training iterations for three different initial points. (B) NPE errors versus training iterations for three different time ranges.

  𝑛 → ℝ 𝑛 and an initial condition 𝑦(𝑡 0 ) = 𝑦 0 . A vector value function 𝑦 ∶ ℝ → ℝ 𝑛 is the state solution, if it satisfies the ODE system 𝑑𝑦 𝑑𝑡 = Φ(𝑦) and the initial condition𝑦(0) = 𝑦 0 .𝑦(𝑡) is called the state at time 𝑡. Given a time range [𝑡 0 , 𝑇 ], 𝑦(𝑇 ) is called the end state on the interval.

	𝑑𝑦 𝑑𝑡 = Φ(𝑦).
	Definition 4 (State solution). Consider an ODE system 𝑑𝑦 𝑑𝑡 = Φ(𝑦), where Φ ∶ ℝ Theorem 1 (Xia & Feng 4 ). Consider an NPE problem, 𝑁𝑃 𝐸(Ω, 𝐺), and let Assumption 1 hold. Given any initial condition,
	𝑦(𝑡 0 ) = 𝑦 0 , the state solution of the ODE system of Eq. (7) converges to the optimal solution of 𝑁𝑃 𝐸(Ω, 𝐺) as time 𝑡 goes to
	infinity, i.e,
	lim

𝑡→∞ 𝑦(𝑡) = 𝑥 * , (

approximate state solution to ODE system NN provides prediction to NPE When When Fully connected network Fully connected network End state Figure 1 Neural network (NN) solution for the ODE system and the NPE. Left: When

  𝑡 ∈ [𝑡 0 , 𝑇 ], the model ŷ(𝑡; 𝑤) itself is considered to be an approximate state solution of the ODE system. Right: When 𝑡 = 𝑇 , the end state of the model, ŷ(𝑡 = 𝑇 ; 𝑤), is a predicted solution to the NPE.

Table 1

 1 Comparison of predicted solutions for Example 1 between approach and the

		approach		The RK method
	Iteration Prediction	NPE error Collocation point Prediction	NPE error
	0	[0.26, 0.43, 0.37, 0.77, 0.97] 24.67	0	[0.00, 0.00, 0.00, 0.00, 0.00] 20.00
	100	[0.17, 3.07, 0.00, 3.64, 4.50] 2.32	100	[0.05, 0.34, 0.22, 0.24, 0.28] 14.60
	1000	[0.13, 2.65, 0.06, 3.13, 3.92] 1.62	1000	[0.15, 1.70, 0.26, 1.82, 2.20] 7.86
	5000	[0.08, 2.76, 0.13, 3.32, 4.25] 0.13	5000	[0.07, 2.71, 0.09, 3.25, 4.15] 0.37
	10000	[0.07, 2.76, 0.08, 3.33, 4.25] 0.08	10000	[0.07, 2.75, 0.08, 3.32, 4.24] 0.03
	30000	[0.07, 2.75, 0.08, 3.32, 4.24] 0.03	30000	[0.07, 2.75, 0.08, 3.32, 4.24] 0.03
	50000	[0.07, 2.75, 0.08, 3.32, 4.24] 0.03	50000	[0.07, 2.75, 0.08, 3.32, 4.24] 0.03
	5.1.2 Nonlinear Complementarity Problem
	Example 2: Consider the following nonlinear complementarity problem:

Table 2

 2 Comparison of predicted solutions for Example 2 between our approach and the RK method

	Our NN approach		The RK method
	Iteration Prediction	NPE error Collocation point Prediction	NPE error
	0	[0.86, 0.00, 0.00, 0.00] 5.04	0	[0.00, 0.00, 0.00, 0.00] 10.00
	100	[1.30, 0.94, 0.00, 2.01] 2.00	100	[0.19, 0.04, 0.05, 0.08] 9.44
	1000	[1.08, 0.00, 0.59, 3.99] 0.14	1000	[1.16, 0.19, 0.20, 0.73] 3.27
	5000	[1.08, 0.00, 0.64, 3.99] 0.01	5000	[1.22, 0.22, 0.00, 2.53] 1.47
	10000	[1.08, 0.00, 0.64, 4.00] 0.00	10000	[1.10, 0.01, 0.40, 3.46] 0.62
	30000	[1.08, 0.00, 0.64, 4.00] 0.00	30000	[1.08, 0.00, 0.64, 3.99] 0.08
	50000	[1.08, 0.00, 0.64, 4.00] 0.00	50000	[1.08, 0.00, 0.64, 4.00] 0.09

Table 3

 3 Comparison of predicted solutions for Example 3 between our approach and the RK method

	Our NN approach		The RK method
	Iteration Prediction	NPE error Collocation point Prediction	NPE error
	0	[1.00, -0.46, 0.24, 1.00]	15.42	0	[1.00, 0.00, 0.00, 1.00]	13.11
	100	[1.77, 0.12, -1.57, 1.66]	11.43	100	[1.00, -0.02, -0.16, 1.00]	13.22
	1000	[27.62, -3.00, -10.00, 10.92] 0.46	1000	[2.61, -0.22, -1.75, 2.28]	12.09
	5000	[28.15, -3.00, -10.00, 11.18] 0.08	5000	[11.38, -1.69, -10.00, 8.88]	10.26
	10000	[28.07, -3.00, -10.00, 11.15] 0.00	10000	[20.78, -3.00, -10.00, 10.87] 7.32
	30000	[28.07, -3.00, -10.00, 11.15] 0.00	30000	[27.96, -3.00, -10.00, 11.15] 0.11
	49999	[28.07, -3.00, -10.00, 11.15] 0.00	49999	[28.07, -3.00, -10.00, 11.15] 0.00
	Table

Table 6

 6 Comparison of NPE errors for NNs with different model sizes. The top half of the table presents the results for singlelayer networks with different numbers of neurons, while the bottom half presents the results for multi-layer networks with neurons per layer.

	Iteration 1 layer, 100 neurons 1 layer, 500 neurons 1 layer, 1000 neurons 1 layer, 1500 neurons 1 layer, 2000 neurons
	100	11.43	0.33	0.23	0.78	2.32
	300	7.08	0.33	0.23	0.75	0.71
	500	2.07	0.33	0.23	0.52	0.71
	1000	0.46	0.33	0.23	0.16	0.29
	3000	0.10	0.00	0.03	0.02	0.02
	5000	0.08	0.00	0.00	0.01	0.01
	10000	0.00	0.00	0.00	0.00	0.01
	Iteration	1 layer, each 500 neurons	2 layer, each 500 neurons	3 layer, each 500 neurons	4 layer, each 500 neurons	5 layer, each 500 neurons
	100	0.33	1.73	2.22	6.43	6.88
	300	0.33	0.10	0.76	1.82	6.70
	500	0.33	0.07	0.27	0.57	0.52
	1000	0.33	0.01	0.03	0.23	0.34
	3000	0.00	0.01	0.01	0.03	0.03
	5000	0.00	0.00	0.00	0.03	0.03
	10000	0.00	0.00	0.00	0.02	0.03
	Number of layers and neurons. Table

Table 7

 7 Comparison of NPE errors for different activation functions. The considered NN model is a two-layer network, with 500 neurons in each layer.

	Iteration	2 layers, each 500 neurons
		tanh sinx sigmoid relu leaky relu
	100	1.73 2.18 3.01	2.62 2.32
	300	0.10 0.73 0.28	1.13 0.71
	500	0.07 0.15 0.05	1.13 0.71
	1000	0.01 0.03 0.01	1.13 0.29
	3000	0.01 0.03 0.01	0.01 0.01
	5000	0.00 0.03 0.01	0.01 0.01
	10000	0.00 0.03 0.00	0.01 0.00

Table 8

 8 Performance of the NN approach for solving the set of NPE problems. STD stands for standard deviation. CPU time is measured in seconds.

	Iteration	M.S.E. loss (Mean ± STD)	CPU time (Mean ± STD)	NPE error (Mean ± STD)	NPE error (50%-quantile)	NPE error (75%-quantile)	NPE error (95%-quantile)
	0	470.29 ± 92.68 0.00 ± 0.00	4839.33 ± 673.10 4720.42	4911.21	6008.41
	100	45.85 ± 39.83	2.36 ± 0.60	15.34 ± 16.46	6.65	19.98	46.38
	500	32.10 ± 26.70	11.32 ± 0.68	0.67 ± 0.89	0.37	0.46	2.30
	1000	26.81 ± 11.53	22.57 ± 1.16	0.31 ± 0.18	0.27	0.36	0.60
	3000	15.20 ± 11.23	67.82 ± 1.40	0.19 ± 0.06	0.18	0.24	0.27
	5000	12.33 ± 10.44	112.91 ± 2.20	0.13 ± 0.07	0.11	0.17	0.24
	7000	8.16 ± 8.13	157.98 ± 2.86	0.10 ± 0.06	0.07	0.14	0.19
	10000	2.93 ± 3.38	225.47 ± 3.82	0.05 ± 0.04	0.04	0.07	0.12
	30000	0.73 ± 1.18	670.38 ± 5.82	0.01 ± 0.01	0.01	0.01	0.02

Table 9

 9 Performance of the RK method for solving the set of NPE problems STD stands for standard deviation. CPU time is measured in seconds.

	Time range	CPU time (Mean ± STD)	NPE error (Mean ± STD)	NPE error (50%-quantile)	NPE error (75%-quantile)	NPE error (95%-quantile)
	[0, 2]	369.33 ± 44.27	1.87 ± 0.16	1.86	2.01	2.08
	[0, 4]	547.93 ± 74.83	0.48 ± 0.06	0.50	0.51	0.55
	[0, 6]	739.74 ± 88.02	0.14 ± 0.02	0.15	0.15	0.16
	[0, 8]	1048.92 ± 88.20	0.04 ± 0.01	0.04	0.05	0.05
	[0, 10]	1542.07 ± 205.43 0.01 ± 0.01	0.01	0.01	0.02
	other hyperparameter settings are kept consistent with those described in Section 5.1. Tables

https://github.com/wuwudawen/IJNME_data_2023/blob/main/Q.npy

DATA AVAILABILITY STATEMENT

The problem data used in Section 5.4 are available in https://github.com/wuwudawen/IJNME_data_2023/blob/main/Q.npy.