domain: info math

text:
Dawen Wu
email: dawen.wu@centralesupelec.fr

Abdel Lisser
email: abdel.lisser@l2s.centralesupelec.fr

A deep learning approach for solving linear programming problems

Keywords: Linear programming, ODE systems, Neural networks, Deep learning

Finding the optimal solution to a linear programming (LP) problem is a long-standing computational problem in Operations Research. This paper proposes a deep learning approach in the form of feed-forward neural networks to solve the LP problem. The latter is first modeled by an ordinary differential equations (ODE) system, the state solution of which globally converges to the optimal solution of the LP problem. A neural network model is constructed as an approximate state solution to the ODE system, such that the neural network model contains the prediction of the LP problem. Furthermore, we extend the capability of the neural network by taking the parameter of LP problems as an input variable so that one neural network can solve multiple LP instances in a one-shot manner. Finally, we validate the proposed method through a collection of specific LP examples and show concretely how the proposed method solves the example.

Introduction

Linear programs are a class of optimization problems where the goal is to minimize a linear objective subject to linear constraints. They play an important role in operations research and have a wide range of real-world applications, such as transportation scheduling, manufacturing planning, and resource allocation.

We consider the LP problem of the following general form

             min x D T x s.t. Ax -b ≤ 0, (1)
where A ∈ R m×n , b ∈ R m and D ∈ R n are problem data. x ∈ R n and u ∈ R m are decision variables and dual 5 variables, respectively. n and m represent the number of decision variables and the number of constraints, respectively. This paper proposes a deep learning approach to solve LP problems, which composes of two parts. 1) The path followed by the decision and dual variables of an LP problem toward the optimal solution is modeled by an ODE system. 2) A feed-forward neural network is used as an approximate state solution to the ODE system.

Literature review

The LP problems have many applications in areas such as economics, engineering, and computer science [START_REF] Dantzig | Linear programming and extensions[END_REF][START_REF] Gass | Linear programming: methods and applications[END_REF]. The simplex and interior-point methods are two classical and still widely used solution methods, that have been well developed and extensively studied in the last decades [START_REF] Karmarkar | A new polynomial-time algorithm for linear programming[END_REF][START_REF] Nocedal | Numerical optimization[END_REF][START_REF] Gondzio | Interior point methods 25 years later[END_REF]. Academic or commercial software implement these two methods for user convenience, such as Cvxpy, Matlab, and Gurobi [START_REF] Diamond | CVXPY: A Python-embedded modeling language for convex optimization[END_REF]Matlab, 2017;Gurobi Optimization, LLC, 2021).

The connection between ODE systems and optimization problems was first studied by [START_REF] Hopfield | neural" computation of decisions in optimization problems[END_REF], where the LP problem is modeled by the Hopfield network. Since then, ODE systems approaches have been applied to various optimization problems, including quadratic programming [START_REF] Xia | A new neural network for solving linear and quadratic programming problems[END_REF][START_REF] Effati | Neural network models and its application for solving linear and quadratic programming problems[END_REF], linear complementarity problems [START_REF] Li-Zhi | A neural network for the linear complementarity problem[END_REF][START_REF] Xia | A recurrent neural network for solving linear projection equations[END_REF], secondorder cone programming [START_REF] Ko | Recurrent neural networks for solving second-order cone programs[END_REF][START_REF] Nazemi | A new neural network framework for solving convex second-order cone constrained variational inequality problems with an application in multifinger robot hands[END_REF], nonlinear projection equations [START_REF] Xia | A new neural network for solving nonlinear projection equations[END_REF] and stochastic game problems [START_REF] Wu | A dynamical neural network approach for solving stochastic two-player zero-sum games[END_REF]. These articles mentioned above use an ODE system to model the optimization problem, and then the optimal solution is obtained by solving the ODE system. These ODE systems are shown to have global convergence properties, i.e., the state solutions converge to the optimal solution as the time variable goes to infinity.

The numerical integration methods for solving ODE systems can be divided into two categories: explicit and implicit methods [START_REF] Teschl | Ordinary differential equations and dynamical systems volume 140[END_REF]. The explicit method determines the later state based on the current state, e.g., RK45, RK23 [START_REF] Dormand | A family of embedded runge-kutta formulae[END_REF][START_REF] Bogacki | A 3 (2) pair of runge-kutta formulas[END_REF]. The implicit method integrate the ODE system by solving equations involving the current and later states, e.g., Radau and BDF [START_REF] Wanner | Solving ordinary differential equations II volume 375[END_REF][START_REF] Shampine | The matlab ode suite[END_REF]. Software has been develop to implement these methods [START_REF] Virtanen | & SciPy 1.0 Contributors[END_REF][START_REF] Rackauckas | Differentialequations.jl-a performant and feature-rich ecosystem for solving differential equations in julia[END_REF]. [START_REF] Dissanayake | Neural-network-based approximations for solving partial differential equations[END_REF] first uses a neural network to approximate the solution of differential equations, where the loss function contains two terms satisfying the initial/boundary condition and satisfying the differential equation. [START_REF] Lagaris | Artificial neural networks for solving ordinary and partial differential equations[END_REF] develop a trail solution method that can ensure initial conditions always be satisfied. [START_REF] Flamant | Solving differential equations using neural network solution bundles[END_REF] takes the parameter of ODE system models as an input variable to the neural network, such that a neural network can be the solution of multiple differential equations, namely solution bundles. The universal approximation theorem is foundational for this line of research, which states that a neural network can approximate any continuous function to arbitrary accuracy [START_REF] Cybenko | Approximation by superpositions of a sigmoidal function[END_REF][START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF]. Python packages are developed to implement these methods [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].

With the rapid growth of computing power and data, there has been a large amount of research on neural networks [START_REF] Goodfellow | Deep Learning[END_REF]. Neural networks have a large number of applications in computer vision, natural language processing, data mining, and other fields [START_REF] Voulodimos | Deep learning for computer vision: A brief review[END_REF][START_REF] Deng | Deep learning in natural language processing[END_REF][START_REF] Liu | A survey of deep neural network architectures and their applications[END_REF]. Most training methods for neural networks are based on gradient descent, and many variants have been proposed to speed up training or to improve the generalization ability [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF][START_REF] Nair | Rectified linear units improve restricted boltzmann machines[END_REF]. Extreme learning machine is a new paradigm for training neural networks that do not involve any gradient-based iterative process [START_REF] Huang | Extreme learning machine: theory and applications[END_REF][START_REF] Zhang | Non-iterative and fast deep learning: Multilayer extreme learning machines[END_REF][START_REF] Xiao | Class-specific cost regulation extreme learning machine for imbalanced classification[END_REF]. This method can apply to a wide range of deep learning tasks. One of our contributions is to link the LP problem to deep learning (DL) so that the LP problem can benefit from a wide range of research in the DL community, e.g., the extreme learning machine. In addition, our paper is expected to inspire new solution method to problems beyond optimization, e.g., optimal control problems [START_REF] Zhang | Asynchronous fault detection for interval type-2 fuzzy nonhomogeneous higher-level markov jump systems with uncertain transition probabilities[END_REF], or to be an integral part of existing methods [START_REF] Djordjevic | Data-driven control of hydraulic servo actuator based on adaptive dynamic programming[END_REF][START_REF] Zhuang | Iterative learning control for repetitive tasks with randomly varying trial lengths using successive projection[END_REF].

Contributions

The contributions of this paper are threefold and summarized as follows.

• We propose a deep learning approach to solve LP problems. The proposed approach uses neurodynamic optimization to model the LP problem by an ODE system. A neural network model is then trained to be an approximate state solution of the ODE system. In this way, the LP problem can be solved only by deep learning, involving neither optimization methods, e.g., simplex or interior point methods, nor numerical integration methods, e.g., RK45.

• We increase the model capacity by treating the parameter of LP problems as an input variable to the neural network. The neural network then becomes the solution to multiple LP instances. A novel algorithm is proposed to train such a neural network. Once trained, the model is able to solve multiple LP instances in a one-shot manner without any iterative process.

• In the experimental Section, we concretely show how to solve a specific LP instance by the proposed method. We explicitly present the LP instance, the corresponding ODE system, and the neural network model. We compare our proposed approach with the RK45 and convex solvers..

Organization

The remaining sections are organized as follows. Section 2 describes how to model an LP problem by an ODE system, and demonstrate the effectiveness of the ODE system from theoretical point of view. Section 3 sketches the case of using one neural network to solve one LP problem. Section 4 shows the case of using one neural network to solve multiple LP problems. Section 5 presents experimental results and discussions.

Section 6 summarizes this paper and gives future directions.

Preliminaries

Karush-Kuhn-Tucker conditions

We consider the LP problem as presented in (1). The corresponding Lagrangian is

L(x, u) = D T x + u T (Ax -b) (2)
The corresponding Karush-Kuhn-Tucker (KKT) conditions are

D + A T u = 0, u T (Ax -b) = 0, Ax -b ≤ 0, u ≥ 0.
(3)

For the LP problem, the KKT conditions are necessary and sufficient for optimality, i.e., the optimal solution

x * of the LP problem (1) can be obtained by finding a solution (x * , u * ) of the KKT condition (3).

ODE system for modeling LP

Let x(t) : R → R n and u(t) : R → R m be functions with respect to the time variable t, and y(t) = (x T (t), u T (t)) T . The following ODE system models the LP problem (1) via the KKT conditions (3).

dy dt = Φ(y) =   dx dt du dt   =   -D + A T (u + Ax -b) + (u + Ax -b) + -u   . (4)
The initial value problem (IVP) is developed by composing of the ODE system (4), an initial condition y(t 0 ) = y 0 and a time range [0, T ], i.e.,

dy dt =   -D + A T (u + Ax -b) + (u + Ax -b) + -u   , y(t 0 ) = y 0 , t ∈ [0, T ] (5)
y(t), t ∈ [0, T ] is called the state solution of the initial value problem (5) if it satisfies both the ODE system and the initial condition. y * is called an equilibrium point of the ODE system (4) if it satisfies Φ(y * ) = 0.

The following theorems and definition show how the ODE system (4) solves the LP problem (1).

Theorem 1 [START_REF] Nazemi | An efficient dynamic model for solving the shortest path problem[END_REF]). y * = (x * , u * ) T is an equilibrium point of the ODE system (4) if and only if it is a satisfied point of the KKT conditions (3).

Definition 1. Let y(t) be a solution of the ODE system dy dt = Φ(y). Then the ODE system is said to be globally convergent to the solution set Θ * = {y * | y * = (x * , u * ) solves (3)} if any state solution with an arbitrary initial point satisfies

lim t→∞ dist (y(t), Θ * ) = 0, where dist (y, Θ * ) = inf y * ∈Θ * ∥y -y * ∥ and ∥ • ∥ denotes the l 2 -norm.
Theorem 2 [START_REF] Nazemi | An efficient dynamic model for solving the shortest path problem[END_REF]). The ODE system (4) is globally convergent to the solution set

Θ * = {y * | y * = (x * , u * ) solves (3)}.
According to Theorem 1, the optimal solution to the LP problem is equivalent to the equilibrium points of the corresponding ODE system. According to Theorem 2, for any given initial point, the state solution of the ODE system will always converge to an optimal solution of the LP problem as the time t goes to infinity, i.e., lim t→∞ y(t) = y * . We refer the reader to [START_REF] Nazemi | An efficient dynamic model for solving the shortest path problem[END_REF][START_REF] Nazemi | A high performance neural network model for solving chance constrained optimization problems[END_REF][START_REF] Xia | A recurrent neural network for solving linear projection equations[END_REF] for the proof of the theorems and other details.

One neural network for solving one LP problem

Before using the neural network, the ODE system requires two parameters to form an IVP, i.e., the initial point y 0 and the time range [0, T ]. Thanks to the global convergence property of the ODE system, the initial point can be simply fixed to zero without loss of convergence, i.e., y(0) = 0.

(6)

Let the neural network model be defined as

ŷ (t; w) = (1 -e -t )NN (t; w) , t ∈ [0, T ], (7)
where NN represents a fully-connected neural network with weights w. The multiplier (1 -e -t ) is used to guarantee the initial condition ŷ(0; w) = 0 must be satisfied regardless of weights w. The multiplier greatly affects the region around t = 0 and gradually converges to 1 as time goes on. The neural network model ŷ (t; w) is act as an approximate state solution to the IVP. The endpoint of the model (7), ŷ (T ; w) = (1 -e -T )NN (T ; w) ,

represents an approximate solution to the KKT conditions.

Figure 1 shows the workflow of solving an LP problem by the neural network model. Consider solving LP problem (1), and the solution procedure can be summarised as follows. 1) The KKT conditions (3) of the LP problem are derived; The objective becomes finding a feasible point to the KKT conditions. 2) The ODE system ( 4) is built to model the KKT conditions; By considering the initial point ( 6) and a time range [0, T ], the ODE system forms an initial value problem.

3) The model ŷ (t; w) as defined in ( 7) is used as an approximate state solution to this initial value problem; The endpoint ŷ (T ; w) is an approximate solution to the KKT conditions.

Since the initial condition ŷ(0; w) = 0 is always satisfied thanks to the multiplier construction, The training objective is to make ŷ(t; w) satisfying the corresponding ODE system by adjusting model's weights w. Thus, the loss function is defined as

E(w) = 1 |T| ti∈T ℓ ∂ ŷ (t i ; w) ∂t i , Φ(ŷ (t i ; w)) , (9)

The initial condition

The time range

The KKT conditions where Φ refers to the ODE system associated with the LP problem to be solved. T ⊂ R + refers to a set of collocation times to be trained, where each collocation time t i is sampled by a uniform distribution over the The training algorithm is unsupervised since there is no pre-given fixed training dataset. Instead, the ODE system is given as an instruction to guide the neural network model. Since each round of training uses fresh data, there is no over-fitting to a particular training data set. Training can go on forever to constantly optimize the neural network's weights w. However, it is realistic to set a stopping criteria to break the loop when the training is deemed good enough. Some possible stopping criteria could be: stop when the number of iterations exceeds a maximum; stop when the loss function or other metrics fall below a certain threshold.

interval [0, T ]. ∂ ŷ(ti;w) ∂ti ∈ R n+m represents the partial derivative of ŷ with respect to t i . Φ(ŷ (t j ; w)) ∈ R n+m
After the training, the neural network model ŷ (t; w) is expected to be well approximate the true state solutiong y(t) on interval [0, T ]. Furthermore, the endpoint ŷ(T ; w) = (x, û) is expected to solve the KKT conditions, where x is the final prediction to the LP problem.

One neural network for solving multiple LP problems

Section 3 presents the case of using one neural network model to solve one LP problem. This section demonstrates how to solve multiple LP problems using a single neural network model, significantly increasing computational efficiency.

Let the neural network model for solving multiple LP problems be defined as

ŷ (t; θ; w) = (1 -e -t )NN (t, θ; w) , (10)
where NN (t, θ; w) is a nerual network that receive the parameter θ of an LP problem as input, as shown in with a parameter θ as input and predict their solutions. As in (7), the initial point is fixed to y(0) = 0, and

(1 -e -t ) ensures the initial conditions are satisfied. When considering only one LP problem, the parameter θ is a constant and can be removed, and the model (10) degraded to (7).

The loss function of ( 10) is defined as

E(w) = 1 |Θ| * |T| θj ∈Θ ti∈T ℓ ∂ ŷ (t i , θ j ; w) ∂t i , Φ j (ŷ (t i , θ j ; w)) , (11)
where Θ refers to the set of parameters θ, and each θ j ∈ Θ relate to an ODE system Φ j . For a specific θ j , ∂ ŷ(ti,θj ;w) ∂tj ∈ R n+m and Φ j (ŷ (t i , θ j ; w)) ∈ R n+m represent the partial derivative of ŷ with respect to t i and the derivative from the ODE system Φ j , respectively. After trainings, the model ŷ (t, θ; w) is expected to be able to solve LP problems with parameter θ coming from the distribution P. Consider a LP problem with parameter θ; The model's final prediction to this LP problem is ŷ(T, θ; w).

Numerical results

Section 5.1 gives the experimental setup, including the environment, hyperparameters, and evaluation metrics. Section 5.2 shows the case of solving a specific LP instance, as introduced in Section 3 whilst Section 5.3 shows the case of solving multiple LP instances, as introduced in Section 4.

Experimental setup

The experiments are conducted on the Google Colab Pro+ platform, which provides 51GB RAM, eight Intel Xeon 2.2GHz CPUs, and one Nvidia Tesla P100 16GB GPU. The neural network model is built by Pytorch 1.9.1 and CUDA 10.2.

The hyperparameters used for training the neural network are as follows

• The optimizer is Adam with a learning rate of 0.001.

• The loss function ℓ(•, •) is the mean square error.

150

• The batch size is 128, and the maximum number of iterations is 500000.

We compare the proposed approach to the ODE solver and the convex solvers. For the ODE solver, we select the RK45 method implemented by Scipy [START_REF] Dormand | A family of embedded runge-kutta formulae[END_REF][START_REF] Virtanen | & SciPy 1.0 Contributors[END_REF]. For the convex solvers, we select the ECOS, OSQP, and Gurobi provided by CVXPY [START_REF] Domahidi | Ecos: An socp solver for embedded systems[END_REF][START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF]Gurobi Optimization, LLC, 2021;[START_REF] Diamond | CVXPY: A Python-embedded modeling language for convex optimization[END_REF].

155

We use the following mean square error (MSE) metric to measure the prediction accuracy, defined by

MSE = 1 B B i=1 (f i (x * i ) -f i (x i )) 2 , ( 12
)
where B refers to the number of instances in a test batch. Consider the i-th instance in the test batch. x * i refers to an optimal solution, and f i (x * i ) refers to the true objective value obtained by the LP solver. x i refers to a prediction, and f i (x i ) refers to the predicted objective value. The evaluation indicator measure on the objective value f i (x i ) rather than on the predicted point xi , because the true objective value is unique f i (x * i ), whereas there may be multiple optimal solutions x * i . We also provide the mean objective value (MOV) to roughly describe a test batch, defined as

MOV = 1 B B i=1 f i (x i ). ( 13
)

Case 1: Solving one LP problem

Consider the following specific LP problem min

x -9.54x 1 -8.16x 2 -4.26x 3 -11.43x 4 s.t. 3.18x 1 + 2.72x 2 + 1.42x 3 + 3.81x 4 ≤ 7.81. ( 14
)
According to (4), the corresponding ODE system for this LP problem is

dy dt = Φ(y) =   dx dt du dt   =            --
+ (u + 3.18x 1 + 2.72x 2 + 1.42x 3 + 3.81x 4 -7.81) + -u            , ( 15
)
where u is the dual variable.

By choosing the initial point to be all-zero and the time range to be [0, 10], the ODE system forms the initial value problem as follows We construct a model ŷ (t; w) as defined in (7) to solve this initial value problem and the LP problem,

dy dt =            --9.54 + 3.18 T (u + 3.18x 1 + 2.72x 2 + 1.42x 3 + 3.81x 4 -7.81) + --8.16 + 2.72 T (u + 3.18x 1 + 2.72x 2 + 1.42x 3 + 3.81x 4 -7.81) + --4.26 + 1.42 T (u + 3.18x 1 + 2.72x 2 + 1.42x 3 + 3.81x 4 -7.81) + --11.43 + 3.81 T (u + 3.18x 1 + 2.72x 2 + 1.42x 3 + 3.81x 4 -7.81) + (u + 3.18x 1 + 2.72x 2 + 1.42x 3 + 3.81x 4 -7.81) + -u            , y(0) =            0 0 0 0 0            , t ∈ [0, 10]. (16)
160
where Table 1 shows the structure of the neural network. As shown in Figure 3a, the loss value decreases from an initial value of 60 to less than 0.1, implying that the model is trained to satisfy the ODE system conditions(15). Figure 3b shows the true state solution obtained by the numerical integration method RK45, and Figure 3c Table 2 shows the final solutions to the LP problem by the three approaches. The solution from the numerical integration method is the state solution's endpoint, i.e., y(10), and the solution of our neural network method is ŷ(10; w). All three approaches achieve the same optimal solutions.

(u + 3.18x 1 + 2.72x 2 + 1.42x 3 + 3.81x 4 -b) + -u            , y(0) =            0 0 0 0 0            , t ∈ [0,
which is able to receive various b ∈ [0, 10] as input, and NN(t, b; w) is a fully connected network as shown in Table (3). Figure 4 shows the training of this model according to Algorithm 2. After 10,000 iterations, the loss value drops from 175 to less than 5, and after 500,000 iterations, it gradually converges to 0.01. Figure 175 5 shows the solutions of the neural network for three specific instances after the training is completed. network can solve all these instances in one shot without suffering from this linear growth. The computational advantage becomes even more significant when the number of instances is large. Additionally, both in the training and prediction stages, our approach is highly conducive to parallelization.

Numbers of instances

D as variable

Consider an LP problem with the following form min

computational time. Together, the figure and the table show how our neural network approach effectively and efficiently handles the situation where D is a variable.

Conclusion

In this paper, we present a deep learning approach in the form of feed-forward neural networks to solve LP problems. We show how a neural network can be used to solve a single LP problem and multiple LP problems, 200 respectively. We demonstrate the effectiveness of the proposed method through specific LP examples. This paper is a beginning of a series of future works. We relate the most basic optimization problem, linear programming, to the machine learning (ML) community so that a large number of cutting-edge research from ML can directly benefit the linear programming study, such as parallel computing and GPU computing.

Since the neural network only gives approximate solutions, to improve the predicted accuracy, we can 205 develop more sophisticated network architecture and training algorithms that take advantage of the problem structure. Another line of future direction could be extending the proposed approach to other optimization or operations research problems, such as quadratic or second-order conic programming.

Figure 1 :

 1 Figure1: The connection between the LP problem, ODE system and neural network.

  represents the partial derivative expected to be according to the ODE system. ℓ(•, •) is an error metric used 105 to measure the difference between these two vectors. Algorithm 1: Solving an LP problem by one neural network Require: A time range [0, T ] Input : An LP problem Output : Solution to the LP problem 1 Function Main: 2 Derive the ODE system Φ of the LP problem.3Construct a neural network model ŷ (t; w) = (1 -e -t )NN (t; w).

4 while True do 5 T

 5 ∼ U (0, T ) Uniformly sample a batch of t from the interval [0, T ] 6 Forward propagation: Compute E(w) based on T 7 Backward propagation: Train the model weight w by ∇E(w) the training of the neural network model to solve an LP problem. The algorithm requires a time range [0, T ], which determines the input space of the model. The LP problem in hand should determine the choice of T , e.g., a larger T for a complex problem. In each training round, the neural network trains on a batch of time T, which would is by the end of the iteration.

  110

Figure 2 :

 2 Figure 2: The neural network model that receives t and θ as input. t refers to the time variable. θ refers to the parameter of an LP problem, and w refers to the weights of the neural network.

Figure 2 .

 2 Figure 2. The parameter θ represents the problem data, such as D, A, b in (1) or a combination of them, and each different θ represents a different LP problem. The model ŷ (t, θ; w) receive different LP problem

Algorithm 2 :

 2 Solving multiple LP problems by one neural network Require: A time range: [0, T ] A probability distribution P for generating θ. 1 Function Main: 2 while True do 3 Θ ∼ P Generate a set of θ according to the given probability distribution P 4 T ∼ U (0, T ) Uniformly sample a batch of t from the interval [0, T ] 5 Forward propagation: Compute E(w) based on the Θ and T 6 Backward propagation: Update w by ∇E(w) the training of the neural network model involving multiple LP problems. The algorithm aims at reducing the loss function (11). At each training iteration, a batch of data Θ and T is sampled to train the model, where Θ is sampled from the given probability distribution P, and T is sampled from the given time range [0, T ].

Figure 3 :

 3 Figure 3: (a) The loss values in the first 5000 iterations (b) The true state solution y(t) to the initial value problem (16) (c) The neural network approximate state solution ŷ(t; w) after training

Figure 4 :

 4 Figure 4: The loss value versus iteration in the first 10,000 iterations

Figure 5 :

 5 Figure 5: Three instances of (17) solved by the RK45 method and the neural network Upper row: True state solutions obtained by the RK45 method Lower row: Approximate state solutions obtained by the neural network Left: For solving the LP problem of (17) with b = 8.16, the endpoints are y(t = 10) = [0.76, 0.65, 0.34, 0.91, 3.00], ŷ(t = 10, b = 8.16; w) = [0.79, 0.67, 0.47, 0.91, 3.03] Middle: For solving the LP problem of (17) with b = 6.72, y(t = 10) = [0.63, 0.53, 0.28, 0.75, 3.00], ŷ(t = 10, b = 6.72; w) = [0.67, 0.56, 0.36, 0.76, 3.04] Right: For solving the LP problem of (17) with b = 6.18, y(t = 10) = [0.58, 0.49, 0.26, 0.69, 3.00], ŷ(t = 10, b = 6.18; w) =[0.64, 0.53, 0.32, 0.71, 3.04]

  81x 4 -7.81) + -u

  A two-layer neural network to solve the LP problem of variable D (20) The input is t ∈ [0, 10] and D ∈ R 4

Figure 6 :

 6 Figure 6: The loss value versus iteration in the first 10,000 iterations

Linear programming ODE system Neural network An LP problem The ODE system Approximate state solution Approximate solution to the KKT conditions

  9.54 + 3.18 T (u + 3.18x 1 + 2.72x 2 + 1.42x 3 + 3.81x 4 -7.81) + --8.16 + 2.72 T (u + 3.18x 1 + 2.72x 2 + 1.42x 3 + 3.81x 4 -7.81) +

	--4.26 + 1.42 T (u + 3.18x 1 + 2.72x 2 + 1.42x 3 + 3.81x 4 -7.81) +
	--11.43 + 3.81 T (u + 3.18x 1 + 2.72x 2 + 1.42x 3 + 3.81x 4 -7.81)

Table 1 :

 1 A two-layer neural network to solve the LP problem (14). The input is t ∈ [0, 10], and the second layer of the neural network has no activation function.

	Layers	Input Size Output Size Activation function
	Fully connected	1	100	Tanh
	Fully connected	100	5

Table 2 :

 2 shows the approximate state solution , i.e., the model ŷ(t; w) after training. The approximate state solution ŷ(t; w) is almost identical to the true state solution y(t), which validates the effectiveness of Solutions of three different approaches to the LP problem(14)

	165
	the training.
	Approaches	x *	u *
	Convex optimization solver: ECOS	[0.7294, 0.6238, 0.3257, 0.8739] 2.9999
	Numerical integration method: RK45 [0.7294, 0.6238, 0.3257, 0.8739] 2.9999
	Neural network	[0.7294, 0.6238, 0.3257, 0.8739] 2.9999

Table 4 :

 4 Comparison between the neural network model, ODE solver and convex solvers on solving multiple instances of (17) The metric MOV is defined in (13). The error metric MSE is defined in (12), where the MSE for the ODE solver and convex solvers are omitted because they provide exact solutions. The CPU time is in seconds, and <= 0.01 means the CPU times are less than or equal to 0.01 seconds.

		Neural network	ODE solver: RK45		Convex solvers
		MOV MSE	CPU time	MOV	CPU time	MOV	ECOS CPU time	OSQP CPU time	GUROBI CPU time
	100	-16.90 0.089 <=0.01 -16.68	3.45	-16.69	0.37	0.61	0.46
	300	-14.83 0.099 <=0.01 -14.57	10.34	-14.58	1.10	1.74	1.29
	500	-14.52 0.101 <=0.01 -14.27	17.38	-14.28	1.89	2.93	2.25
	1000	-15.40 0.095 <=0.01 -15.16	34.50	-15.16	3.74	5.90	4.36
	3000	-15.22 0.095 <=0.01 -14.98	104.35	-14.99 10.97	17.45	13.26
	5000	-15.39 0.095 <=0.01 -15.15	173.35	-15.16 18.60	29.28	22.13
	10000	-15.25 0.095 <=0.01 -15.01	346.25	-15.02 37.28	57.99	44.73

Table 4

 4

reports our neural network model's CPU times and accuracy along with two comparable ap-

  = [D 1 , D 2 , D 3 , D 4 ] T ∈ R 4 is generated by a random number from the uniform distribution U (0, 10) multiplied by the vector [3.18, 2.72, 1.43, 3.81] T . Every different D represents a different LP problem. The--D 1 + 3.18 T (u + 3.18x 1 + 2.72x 2 + 1.42x 3 + 3.81x 4 -7.81) + --D 2 + 2.72 T (u + 3.18x 1 + 2.72x 2 + 1.42x 3 + 3.81x 4 -7.81) + --D 3 + 1.42 T (u + 3.18x 1 + 2.72x 2 + 1.42x 3 + 3.81x 4 -7.81) + --D 4 + 3.81 T (u + 3.18x 1 + 2.72x 2 + 1.42x 3 + 3.81x 4 -7.81) +(u + 3.18x 1 + 2.72x 2 + 1.42x 3 + 3.

			x	D T x
			s.t.	(20)
			3.18x 1 + 2.72x 2 + 1.42x 3 + 3.81x 4 ≤ 7.81,
	where D corresponding initial value problem is
			
	dy dt	=	         

which is able to receive various D as input, and NN(t, D; w) is a fully connect network as shown in Table 190 (5). Figure 6 shows the training of this model according to Algorithm 2. After 10,000 iterations, the loss value decreases from 600 to less than 10. Compared to the case of b as a variable, D being a variable is more challenging because of the higher dimensional of the input space.


Dawen Wu
email: dawen.wu@centralesupelec.fr

Abdel Lisser
email: abdel.lisser@l2s.centralesupelec.fr

Enhancing Neurodynamic Approach with Physics-Informed Neural Networks for Solving Non-Smooth Convex Optimization Problems

Keywords: Non-smooth convex optimization problem, Neurodynamic optimization, Physics-informed neural network, Numerical integration method, Ordinary differential equation

This paper proposes a deep learning approach for solving non-smooth convex optimization problems (NCOPs), which have broad applications in computer science, engineering, and physics. Our approach combines neurodynamic optimization with physics-informed neural networks (PINNs) to provide an efficient and accurate solution. We first use neurodynamic optimization to formulate an initial value problem (IVP) that involves a system of ordinary differential equations for the NCOP. We then introduce a modified PINN as an approximate state solution to the IVP. Finally, we develop a dedicated algorithm to train the model to solve the IVP and minimize the NCOP objective simultaneously. Unlike existing numerical integration methods, a key advantage of our approach is that it does not require the computation of a series of intermediate states to produce a prediction of the NCOP. Our experimental results show that this computational feature results in fewer iterations being required to produce more accurate prediction solutions. Furthermore, our approach is effective in finding feasible solutions that satisfy the NCOP constraint.

Introduction

Non-smooth convex optimization problems (NCOPs) are an important class of optimization problems that arise in various real-world applications. While traditional optimization methods have been successful in solving many smooth convex optimization problems, they are not well-suited to handle the nonsmoothness of the objective functions in NCOPs. This has led to the development of alternative approaches, such as neurodynamic optimization.

Neurodynamic optimization is a promising approach for solving NCOPs [START_REF] Qin | A two-layer recurrent neural network for nonsmooth convex optimization problems[END_REF]. It involves the implementation of a circuit-based neurodynamic model to solve the optimization problem in real time.

This approach has many potential applications, such as in resource allocation [START_REF] Jia | An adaptive continuous-time algorithm for nonsmooth convex resource allocation optimization[END_REF], feature selection [START_REF] Wang | Neurodynamics-driven supervised feature selection[END_REF], and coordination of multi-manipulator systems [START_REF] Hou | Distributed adaptive coordinated control of multimanipulator systems using neural networks[END_REF].

The use of neurodynamic optimization to solve NCOPs typically involves three steps. First, a first-order ordinary differential equation (ODE) system is used to model the NCOP. The ODE system must be shown to have a globally convergent property, meaning that the state solution of the system converges to the optimal solution of the NCOP as time approaches infinity, regardless of the starting point. Second, classical numerical integration methods such as Runge-Kutta methods or backward differentiation formulas are used to numerically obtain the state solution of the ODE system [START_REF] Burden | Numerical analysis[END_REF]. Third, the state at the end of a given time range, called the end state, is taken as to be the predicted solution of the NCOP. However, this approach has a drawback. That is, the numerical integration method has to compute all the intermediate states, from the initial state to the end state, to solve the optimization problem, which makes it inefficient.

To overcome this limitation, we propose a novel approach that combines neurodynamic optimization with physics-informed neural networks to provide an efficient and accurate solution. Unlike the traditional numerical integration method, our improved approach can directly predict the entire state solution, including the end state. This approach allows for efficient and accurate optimization of NCOPs without requiring the computation of all intermediate states.

Here we want to clarify the difference between our work and Neural ODEs. [START_REF] Chen | Neural ordinary differential equations[END_REF]. Neural ODEs are a class of neural networks that interpret transformations within the model as a continuous process and use ODEs to characterize this process. Training such a model involves solving an ODE using numerical integration methods, and recent research bypasses the use of numerical integration methods to improve computational efficiency [START_REF] Biloš | Neural flows: Efficient alternative to neural odes[END_REF]. However, it is important to emphasize that the primary goal of a neural ODE is not to solve a particular ODE per se, but rather to address traditional supervised learning challenges, such as time series modeling. In our case, the task at hand is to solve an ODE system derived from an NCOP. In this context, the neural network acts as a surrogate model that is specifically trained to approximate the solution of the ODE system.

Related Works

Neuraodynamic Optimization. Neurodynamic optimization is a class of methods that model constrained optimization problems using ODE systems. This approach was first introduced by [START_REF] Hopfield | Neural" computation of decisions in optimization problems[END_REF] to solve the traveling salesman problem. Since then, neurodynamic optimization has been applied to a wide range of optimization problems, including linear and quadratic programming problems [START_REF] Xia | A recurrent neural network for solving linear projection equations[END_REF], general convex programming problems [START_REF] Kennedy | Neural Networks for Nonlinear Programming[END_REF][START_REF] Xia | A novel recurrent neural network for solving nonlinear optimization problems with inequality constraints[END_REF], biconvex optimization problems [START_REF] Che | A two-timescale duplex neurodynamic approach to biconvex optimization[END_REF], global optimization problems [START_REF] Che | A collaborative neurodynamic approach to global and combinatorial optimization[END_REF], and pseudoconvex problems [START_REF] Guo | A one-layer recurrent neural network for pseudoconvex optimization subject to linear equality constraints[END_REF][START_REF] Bian | Neural network for nonsmooth pseudoconvex optimization with general convex constraints[END_REF][START_REF] Liu | A one-layer recurrent neural network for nonsmooth pseudoconvex optimization with quasiconvex inequality and affine equality constraints[END_REF]. These methods typically use the Lyapunov stability theorem to prove that the constructed ODE system has a global convergence property, which means that any state solution of the ODE system converges to an optimal solution of the target problem.

Physics-Informed Neural Networks (PINNs). Another line of research included in our work is the use of deep learning to solve differential equations. As a representative approach, 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] are a class of neural networks that integrate physical laws described by PDEs and boundary conditions into the training losses. Thus, the neural network is trained to become the solution of the PDE. PINNs have a wide range of successful applications, such as identifying material properties (Shukla et al., 2021a), optimizing sensor locations [START_REF] Jagtap | Deep kronecker neural networks: A general framework for neural networks with adaptive activation functions[END_REF], and investigating complex fluid dynamics, including supersonic and high-speed flows (Jagtap et al., 2022a;[START_REF] Mao | Physics-informed neural networks for high-speed flows[END_REF]. To facilitate the use and application of PINNs, several software packages have been developed and are available to the scientific community [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].

Numerous methods and theoretical explorations have been introduced to improve the computational performance of PINNs. From the theoretical side, research has been conducted to understand the error estimates associated with PINNs. [START_REF] Mishra | Estimates on the generalization error of physics-informed neural networks for approximating a class of inverse problems for pdes[END_REF] investigates the generalization error of PINNs in approximating a class of inverse problems (Mishra &[START_REF] Mishra | Estimates on the generalization error of physics-informed neural networks for approximating a class of inverse problems for pdes[END_REF][START_REF] De Ryck | Error estimates for physics informed neural networks approximating the navier-stokes equations[END_REF] provides error estimates for the PINN approximation of the Navier-Stokes equations. On the methodological side, a prominent strategy is the development of domain decomposition methods, which divide the problem domain into several subdomains, each of which is addressed by a separate PINN. This approach creates interconnected PINN solutions that allow independent treatment of each subdomain while preserving the interconnections across the partitioned domains (Jagtap et al., 2020c;[START_REF] Shukla | Parallel physics-informed neural networks via domain decomposition[END_REF]Hu et al., 2022a;[START_REF] Penwarden | A unified scalable framework for causal sweeping strategies for physics-informed neural networks (pinns) and their temporal decompositions[END_REF][START_REF] Hu | When do extended physics-informed neural networks (xpinns) improve generalization[END_REF]. To further improve the performance of PINNs, adaptive activation functions have been investigated, which add learnable parameters to the activation function [START_REF] Jagtap | Adaptive activation functions accelerate convergence in deep and physics-informed neural networks[END_REF][START_REF] Jagtap | Deep kronecker neural networks: A general framework for neural networks with adaptive activation functions[END_REF]Jagtap & Karniadakis, 2023). These innovations in PINNs are reshaping the landscape of solution methods for solving complex scientific and engineering problems.

Contributions

Our paper presents several key contributions:

• We propose a novel approach for solving NCOPs that enhances neurodynamic optimization with PINNs, avoiding the use of numerical integration methods.

• We improve the PINN to adapt to the NCOP by incorporating the initial condition into the neural network and using only one neural network to predict all system states.

• We design a specialized training algorithm that leverages the structure of NCOPs to optimize the performance of our proposed approach. The neural network is trained to simultaneously satisfy the ODE system and minimize the NCOP objective function.

• Our experimental results demonstrate several advantages of our approach over classical numerical integration methods and vanilla PINNs. Our approach requires fewer iterations to produce more accurate predicted solutions, and it finds a feasible solution to the problem more quickly.

Outline

The remaining sections are organized as follows: Section 2 provides the necessary background for understanding the paper, including an introduction to the NCOP and the neurodynamic optimization used to model it. In Section 3, we present our proposed neural network for solving the NCOP. Section 4 details the design of the loss function and the training algorithm for the neural network. Section 5 presents the results of our experiments on solving the NCOP using the proposed approach, and compares our approach with the existing methods. Finally, Section 6 summarizes the main results of the paper and outlines possible directions for future research.

Neurodynamic Approach for Modelling NCOP

Section 2.1 introduces the NCOP, which is the goal of this paper. In Section 2.2, we show how to model the NCOP using a neurodynamic approach. In Section 2.3, we introduce the PINN method for solving PDEs.

NCOP

We consider the following optimization problem:

                   min x f (x) s.t. g(x) ≤ 0 Ax = b, (1)
where x = (x 1 , x 2 , . . . , x n ) T ∈ R n is the decision variables, and f : R n → R is the objective function.

g(x) = (g 1 (x), g 2 (x), . . . , g m (x)) T : R n → R m represents the inequality constraints, and u = (u 1 , u 2 , . . . , u m ) represents the dual variables of the inequality constraints. Ax = b represents the equality constraints with A ∈ R k×n and b ∈ R k . n, m, and k denote the number of decision variables, inequality constraints, and equality constraints, respectively.

In this paper, we consider the case where f (x) and g(x) are convex but not necessarily smooth, and A is of full row rank. We denote x * and u * as the optimal primal and dual solutions, respectively.

Definition 1 (Subgradient and subdifferential).

A vector l ∈ R n is a subgradient of f : R n → R at x ∈ dom f if the following holds f (z) ≥ f (x) + l T (z -x), ∀z ∈ dom f. ( 2
)
The set of all subgradients of f at x is called the subdifferential of f at x and is denoted by ∂f (x).

Neurodynamic Approach

Now, let x : R → R n and u : R → R m be some time dependent functions. The aim of a neurodynamic approach is to construct a first-order ODE system to govern x(t) and u(t), such that they will settle down to the optimal primal and dual solutions of the NCOP (1). In this paper, the two-layer neurodynamic approach in Qin & Xue (2014) is adopted, where the ODE system is described as follows:

dx dt ∈ -(I -U ) ∂f (x) + ∂g(x) T (U + g(x)) + -A T ρ(Ax -b), du dt = 1 2 -U + (U + g(x)) + , (3)
where U = A T AA T -1 A, I ∈ R n×n is the identity matrix, ρ(s) = (ρ (s 1 ) , ρ (s 2 ) , . . . , ρ (s k ))

T , and for

i = 1, 2, . . . , k, ρ (s i ) =            1 if s i > 0, [-1, 1] if s i = 0, -1 if s i < 0. ( 4
)
To simplify the discussion, we denote y(t) = (x(t) T , u(t) T ) T and define:

Φ(y) =   -(I -U ) ∂f (x) + ∂g(x) T (U + g(x)) + -A T ρ(Ax -b) 1 2 (-U + (U + g(x)) + )   . (5)
Thus, the ODE system (3) can be written as dy dt ∈ Φ(y).

Definition 2 (State solution). Consider an ODE system dy dt ∈ Φ(y), where Φ : R n → R n . Given (t 0 ∈ R, y 0 ∈ R n ), a vector value function y : R → R n is called a state solution, if it satisfies the ODE system dy dt ∈ Φ(y) and the initial condition y(t 0 ) = y 0 .

In particular, we call y(t) the state at time t. Given a time interval [t 0 , T ], we call y(T ) the end state on that time interval.

Theorem 1 [START_REF] Qin | A two-layer recurrent neural network for nonsmooth convex optimization problems[END_REF]). Consider a NCOP (1) and its derived ODE system (3). Given any initial condition y(t 0 ) = y 0 , the state solution y(t) of the ODE system converges to an optimal solution y * as time t approaches infinity, i.e.

lim t→∞ y(t) = y * , (6)
where y * = (x * T , u * T ) T , x * and u * are the optimal primal and dual solutions to the NCOP.

In particular, if the NCOP contains only one optimal solution x * , then the ODE system is called globally asymptotically stable at y * .

Initial Value Problem (IVP) Construction. In practice, in order to use the neurodynamic approach to solve the NCOP, we need to construct an IVP consisting of three components: 1) the ODE system (3), 2)

an initial condition y(t 0 ) = y 0 , and 3) a time range t ∈ [t 0 , T ]. y(t) for t ∈ [t 0 , T ] represents the state solution of this IVP problem over the time range [t 0 , T ], where the end state, y(T ), is considered to be the predicted solution to the NCOP. According to Theorem 1, the larger the time range [t 0 , T ], the closer y(T ) is to the optimal solution y * of the NCOP.

Physics-Informed Neural Networks (PINNs)

In this subsection, we introduce the powerful framework of PINNs for solving PDEs. For clarity and to make this subsection self-contained, please note that the mathematical notations used in this subsection are independent of those used in other sections. The mathematical notations should be interpreted only in the context of this subsection.

A typical PDE problem can be expressed in the following general form:

D x (u; λ) = f (x), x ∈ Ω ⊂ R d B k (u) = g k (x), x ∈ ∂Ω ⊂ R d , for k = 1, 2, . . . , n b , (7)
where D x (•) is the differential operator, and u : Ω ∪ ∂Ω → R is the solution to be found. λ denotes the PDE parameters. B k (•) denotes to the boundary conditions, which can be of the Dirichlet, Neumann, or mixed type. For problems involving temporal dynamics, time t is considered as part of x, and the initial conditions can be treated as a unique type of boundary condition.

Let {x

(i) b , u (i) b } N b i=1 and {x (i) d } N d i=1
be the sets of randomly selected training points and residual points, respectively. These points are usually drawn from an unknown distribution. Let u Θ be a surrogate function based on a neural network with model parameters Θ. The goal of the PINN approach is to learn a surrogate function u Θ to approximate the solution u for a given PDE.

The loss function for PINNs can be expressed as:

L(Θ) = W b MSE b Θ; {x (i) b , u (i) b } N b i=1 + W d MSE d Θ; {x (i) d } N d i=1 , (8)
where W b and W d represent the weights for the data and residual losses, respectively. MSE b Θ; {x

(i) b , u (i) b } N b i=1 + and MSE d Θ; {x (i) d } N d
i=1 are given by:

MSE b Θ; {x (i) b , u (i) b } N b i=1 = 1 N b N b i=1 u (i) -u Θ x (i) u 2 , MSE d Θ; {x (i) d } N d i=1 = 1 N d N d i=1 D x (u Θ (x (i) d ); λ) -f (x (i) d ) 2 , (9)
where MSE b measures the data mismatch term, which enforces the boundary conditions as constraints.

MSE d evaluates the PDE residual at a finite set of collocation points. The neural network parameters Θ are determined by minimizing the loss function in (8).

A distinct advantage of PINNs is their ability to integrate both experimental and synthetic training data into the loss function. One of the central mechanisms facilitating this integration is the computation of the PDE residual using automatic differentiation [START_REF] Baydin | Automatic differentiation in machine learning: a survey[END_REF], which provides an efficient and accurate evaluation of the PDE operator. Thus, PINNs transform the task of solving a PDE into an optimization problem, where the global minimum of the loss function signifies the solution to the PDE. This capability positions PINNs as a grid-free method, alleviating the often cumbersome process of mesh generation.

Remark 1. The ODE system (3) under consideration is a particular case of the PDE problem (7). The key differences are as follows:

• The solution of the ODE system is a function with only one input, while the PDE solution is a function with multiple inputs.

• The ODE system involves multiple differential equations, each corresponding to a state dynamic, whereas the PDE typically involves only a single equation with multiple variables.

Mathematically, PINNs typically aim to find a solution u : R d → R, where d is the number of variables.

In contrast, the solution sought in the ODE system is u : R → R n , where n is the number of differential equations in the system.

Modified PINN Model

NN as solution of ODE system NN provides prediction to NCOP

When

Fully connected network

When

Fully connected network

Figure 1: Neural network solution for the ODE system and the NCOP. Left: When t ∈ [t 0 , T ], the model ŷ(t; w) itself is considered to be an approximate state solution of the ODE system. Right: When t = T , the model returns the prediction, (x, û) = ŷ(t = T ; w), where x and û represent the primal and dual predictions for the NCOP, respectively.

Model Description. We propose a modified PINN model to solve the NCOP. Our model can be expressed by the following equation:

ŷ(t; w) = y 0 + 1 -e -(t-t0) N(t; w), t ∈ [t 0 , T ], (10)
where N(t; w) is a fully connected neural network with trainable parameters w. y 0 is a given initial point for the ODE system. [t 0 , T ] is a given time range. The auxiliary function 1 -e -(t-t0) ensures that the neural network always satisfies the initial condition ŷ(t = t 0 ; w) = y 0 regardless of w.

Improvements to PINN. The proposed model ( 10) can be seen as a modified version of the PINN to adapt to the considered problem. In particular, we make three improvements. First, we use the Lagaris method to handle the initial conditions of the ODE system [START_REF] Lagaris | Artificial neural networks for solving ordinary and partial differential equations[END_REF]. Second, we use only one neural network to handle all the differential equations of the entire ODE system, whereas, in the original method, each neural network handles only one differential equation of the system. Finally, we use the exponential form for the auxiliary function, which has been shown to improve convergence, as demonstrated in previous studies [START_REF] Mattheakis | Hamiltonian neural networks for solving equations of motion[END_REF]. Our experiments confirm this improvement in convergence, which is likely due to the reduced influence of the neural network as time progresses further from the initial time.

Approximate State Solution to the ODE. As shown in Figure 1 (Left), the proposed model ( 10)

itself is an approximate state solution of the ODE system (3) on the time range [t 0 , T ], i.e,

ŷ(t; w) ≈ y(t), t ∈ [t 0 , T ], (11)
where y(t) is the true state solution of the ODE system. While the input time t of the model ŷ(t; w) can be take any real number, we specifically use ŷ(t; w) as the solution of the ODE over the time range [t 0 , T ]. As a result, we restrict the input of ŷ(t; w) to the time range t ∈ [t 0 , T ].

Predicted Solution to the NCOP. The end state of the proposed model, i.e., ŷ(t = T ; w), is used as the predicted solution to the NCOP (1), as shown in Figure 1 (Right). The following equation shows how ŷ(t = T ; w) approximates the optimal solution y * :

ŷ(t = T ; w) ≈ y(T ) ≈ y * . (12)
Here, ŷ(t = T ; w) ≈ y(T ) indicates that the end state of our model approximates the true end state, and y(T ) ≈ y * comes from Theorem 1, indicating that the true end state is the predicted solution of the NCOP.

Model training

Training Objective

Loss Function. We define the loss function of the proposed model (10) as follows:

L(t, w) = ∂ ŷ(t; w) ∂t -Φ(ŷ(t; w)) , (13)
where Φ(•) refers to the ODE system (3), which corresponds to the NCOP to be solved. ∥•∥ is the Euclidean norm. Φ(ŷ(t; w)) is the expected derivative according to the ODE system. ∂ ŷ(t;w) ∂t is 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, high-performance deep learning library[END_REF][START_REF] Bradbury | JAX: composable transformations of Python+NumPy programs[END_REF]. L(t, w) represents the difference between the two at time t and with network parameters w.

Embedding the NCOP into the Loss Function. The NCOP is integrated into the loss computation process through the ODE system rather than as a component of the neural network. A neural network is created as an empty framework without a specific goal to solve a particular NCOP. Instead, by reformulating the NCOP as an ODE system and embedding it into the loss function, the neural network is trained toward solving the NCOP.

Objective Function. The goal of training the proposed model is to minimize the following objective function: Batch Loss. However, the objective function J(w) is computationally intractable to compute due to its integral part. Therefore, in practice, we train the model by minimizing the following batch loss:

J(w) = T t0 L(t, w)dt, (14)
L(T, w) = 1 |T| t∈T L(t, w), (15)
where T is a set of randomly sampled time points from the interval [t 0 , T ], and |T| denotes the size of this set.

In this way, we can approximate the integral in the objective function J(w) by a sum of loss values over the set of sampled times. By minimizing the batch loss, we can effectively train the model to solve the NCOP.

Algorithm Design

Objective value under Constraints (OuC) Metric. We introduce an evaluation metric called OuC to measure how well a predicted solution, x pred , solves the NCOP:

OuC(x pred ) =      f (x pred ) if x pred ∈ Ω, +∞ otherwise , ( 16
)
where Ω is the feasible set defined as Ω = {x|x ≤ g(x), Ax = b}. The OuC metric evaluates a predicted solution by setting OuC to positive infinity or a very large real number when the predicted solution is not feasible, and to the objective value when it is feasible.

Projection Mapping onto Equality Constraints. To increase the likelihood of OuC(x pred ) being a real value instead of infinite, we employ the following projection function to map x pred to the set of equality constraints, By definition, the evaluation metric OuC(x pred ) attains a real value only when x pred lies within the feasible set Ω. As illustrated in Figure 2, there are two circumstances when P eq aids in projecting x pred onto the feasible set:

P eq (x pred ) = x pred -A T AA T -1 (Ax pred -b). (17)
• x pred satisfies the inequality constraints but fails to meet the equality constraints, i.e., x pred ∈ I -Q.

• x pred does not satisfy both the inequality and equality constraints, i.e., x pred / ∈ I ∪ Q.

In both of the above scenarios, P eq has a reasonable chance of mapping x pred onto the feasible set Ω, resulting in the value of OuC being a real number rather than infinite.

Algorithm 1: Deep learning solver for convex optimization problems Input : A COP; A time range [t 0 , T ]; An initial condition (t 0 , y 0 ) Output: Predicted primal and dual solutions xbest and ûbest . 1 Function Main: This mechanism ensures that the best prediction obtained by the model is maintained throughout the training process, improving the overall performance of the algorithm.

Neural Network Solution Procedure. Our approach employs gradient descent on the batch loss (15)

to improve the neural network prediction at each iteration. Assume that the maximum number of training iterations is M . The evolution of the neural network is represented by ŷ(t; w 1 ), ŷ(t; w 2 ), . . . , ŷ(t; w M ), where w i and ŷ(t; w i ) denote the network parameters and the approximate state solution at the i-th iteration, respectively. The predicted end states are ŷ(t = T ; w 1 ), ŷ(t = T ; w 2 ), . . . , ŷ(t = T ; w M ), where ŷ(t = T ; w i ) represents the prediction to the NCOP at the i-th iteration.

Experiments

Neural Network Setup. To implement our proposed model and the ODE system, we used PyTorch 1.12.1 with CUDA 11.2 [START_REF] Paszke | PyTorch: An imperative style, high-performance deep learning library[END_REF] and JAX 0.4.1 [START_REF] Bradbury | JAX: composable transformations of Python+NumPy programs[END_REF]. Our neural network architecture consisted of a single layer fully connected network with 100 neurons and a Tanh activation function.

For training, we used the Adam optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] with a learning rate of 0.001, and a batch size of 128.

Evaluation metrics. We used two metrics to measure the performance of the model: 1) The OuC metric, as defined in equation ( 16). 2) Computational time, which measures the time required to achieve a desired accuracy.

Comparisons with Numerical Integration Methods

In this subsection, we compare our proposed method with six classical numerical integration methods:

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 (Virtanen et al., 2020). The RK45, RK23, and DOP853 are explicit Runge-Kutta methods, while the BDF and Radau are implicit methods. LSODA is an adaptive method that automatically switches between explicit and implicit methods depending on the stiffness of the ODE system.

Example 1: We aim to solve the following NCOP:

min x f (x) = 10(x 1 + x 2 ) 2 + (x 1 -2) 2 + 20|x 3 -3| + e x3 s.t. g(x) = (x 1 + 3) 2 + x 2 -36 ≤ 0 h(x) = 2x 1 + 3x 2 + 5x 3 -7 = 0. ( 18
)
The feasible set of this problem is convex, and the objective function is convex but non-smooth due to its inclusion of absolute values. 1,1] Neural Network Solution. Figure 3 shows the solution process of our approach. Each horizontal row represents the evolution of a neural network trained by Algorithm 1. Each sub-figure shows ŷ(t; w) for t ∈ [0, 10], which represents the approximate state solution of the IVP at a given training iteration.

Construction of

Here, we must emphasize that ŷ(t; w) is implemented by one neural network with four output units, and ŷ(t; w) = (ŷ 1 (t; w), ŷ2 (t; w), ŷ3 (t; w), ŷ4 (t; w)). In particular, the end state has ŷ(t = 10; w) =   ŷ1 (t = 10; w), ŷ2 (t = 10; w), ŷ3 (t = 10; w) =x , ŷ4 (t = 10; w)

=û   , (19)
where x and û represent the predicted primal and dual solutions, respectively, of the problem (18). OuC Performance. Figure 4 presents a comparative analysis of the OuC drop rates for our proposed method and six traditional numerical integration methods over three different initial points. The following observations can be made:

• Our method outperforms traditional numerical integration methods in terms of OuC reduction, as evidenced by the lower OuC values achieved in fewer iterations. For example, for the initial point (IP) [0, 0, 0, 0], our method reduces OuC from 43.66 to 28 in just five iterations, whereas the best-performing numerical integration methods, namely RK45, DOP853, and Radau, require 20 iterations to achieve comparable results. Similar results are observed for the other two initial points.

• The speed of convergence varies for different initial points, with our method showing greater robustness to different initial settings compared to numerical integration methods. While most numerical integration methods converge faster for the IP [1, 0, -2, 3] and slower for [0, 0, 0, 0], our method converges at approximately the same rate for both.

• The starting OuC values differ for various initial points, with the IP [1, 0, -2, 3] having the highest starting OuC value of around 135, and [0, 0, 0, 0] having the lowest starting OuC of approximately 43. However, the results demonstrate that the OuC values do not significantly affect the speed of convergence.

Why Our Approach has Better OuC than RK45. It is important to emphasize that our approach has no advantage over classical numerical integrators such as RK45 when it comes to solving for the full solution on the IVP, i.e. the function y. As shown in Figure 4, our method outperforms these numerical integrators in the OuC metric because our approach focuses on improving the end point ŷ(t = T ; w), rather than the entirety of the function ŷ. It is also worth noting that the OuC performance metric only considers the endpoint ŷ(t = T ; w), not the entire function ŷ. The methodology we propose, described in sections 3 and 4, is deliberately designed with this specific goal in mind. Figure 5: Comparison of our method with the numerical integration methods in terms of computational efficiency. Time is measured in seconds.

Computational Time Performance. Figure 5 shows the time needed by different solution methods to obtain an acceptable solution (i.e., OuC ≤ 28) to Problem (18). Our proposed method outperforms all considered numerical integration methods in terms of computational efficiency. The most efficient numerical integration method, BDF, still requires 10 times more computational time than our method to obtain a satisfactory solution. Moreover, the computational time of our method is less affected by different IP settings, requiring about 1.78 seconds for all three initial points. In contrast, some numerical integration methods, such as Radau and LSODA, show significant variations in computation time at different initial points.

Comparisons with PINN

In this subsection, we perform an ablation study comparing the proposed method with 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] and Lagaris method [START_REF] Lagaris | Artificial neural networks for solving ordinary and partial differential equations[END_REF] to validate the effectiveness of our method.

Example 2:

min x f (x) = |2.3x 1 + x 3 -3.5| + |x 2 + 2x 3 -1.8| + |1.3x 1 + x 2 + x 3 + 3| s.t. g(x) = x 2 1 -x 2 + x 3 + 3 ≤ 0 h 1 (x) = x 1 + x 2 + x 3 = 0, h 2 (x) = 2x 2 + x 3 = 0. ( 20
)
We aim to solve the NCOP (20). The problem has a convex feasible set and a convex but non-smooth objective function due to its inclusion of absolute values.

Construction of IVPs. We model problem (20) as an ODE system (3) and set the time range as [t 0 , T ] = [0, 10]. To construct three IVPs, we choose three initial points, namely [1, -1, 0, 3], [2, 3, -2, 1], and

[2, -2, 1, -2].
Based on these initial points, we instantiate three proposed neural networks (10) as approximate state solutions.

Experimental Setup. We compare our proposed approach with two methods: vanilla PINN and PINN with Lagaris construction method (PINN+Lagaris). Our approach can be regarded as the PINN+Lagaris method enhanced by the ORR mechanism (PINN+Lagaris+ORR). The hyperparameters and training details are the same as those in Section 5.1. Table 1 shows the performance of the three methods over the first one hundred iterations, while Figure 6 shows their convergence behavior over the first one thousand iterations. The results reveal the following key observations:

• Our proposed approach yields an excellent predicted solution within the first 20 training iterations, consistent with the results presented in Section 5.1. In contrast, even after 1000 iterations, neither the PINN nor the PINN+Lagaris methods achieves a predicted solution that compares favorably to that of our approach.

• Our method has a higher probability of obtaining feasible solutions. As shown in method returns 'inf' 11 times for the three IP configurations, while the Lagaris method reduces this occurrence to 7. In contrast, our method returns 'inf' only twice, both times in the first round, indicating that it can reach a feasible solution more quickly.

• Neither the PINN method nor the PINN+Lagaris method maintains an optimal solution during the optimization process. As shown in Table 1, the PINN method and the PINN+Lagaris method achieve good OuC values of 12.33 and 12.04, respectively, for the IP [2, 3, -2, 1] at the 100th iteration. However, neither method maintains this level of performance, and their OuC values increase in subsequent iterations.

• The Lagaris method can improve the performance of vanilla PINN, as shown in Figure 6. However, this improvement is not substantial and varies depending on the IP configuration. For example, the improvement is significant for the first and third initial points but negligible for the second IP.

Hyperparameter Study

In this subsection, we perform a hyperparameter study on the following NCOP.

Example 3:

min x f (x) = ∥Cx -d∥ 1 s.t. g 1 (x) = x 2 1 -x 2 + x 3 + x 5 -x 8 -10 ≤ 0, g 2 (x) = |x1 -x3 + x4 + x7| -4.8 ≤ 0, h(x) = x 1 + x 3 + x 5 + x 7 -1 = 0, (21)
where ∥•∥ 1 denotes the L 1 norm, and

C =         1 4 2 2 1.3 4 2 1 2.8 2 1.6 3.2 0 2 1 1 1 4 2.3 2 2.5 0 5 1 1 1 1 3.1 2.3 0 0.8 1         , d =         1.5 -3.8 6.2 7.5         . ( 22
)
Example 3 involves a nonsmooth objective function and a nonsmooth inequality constraint g 2 (x).

To set up the algorithm, we choose the IP as an all-ones vector and the time range as [0, 10]. In the following, we discuss the computational performance for different neural network sizes and different learning rates.

Model Size. In Figures 7 (A) and (B), we investigate the computational performance of neural networks with various widths and depths, and the optimal result is obtained with a 700-neuron-wide, single-layer structure. In (A), with a maximum of 1000 training iterations, networks with fewer neurons (such as 100, 300, and 500) underperform due to their model capacity, limiting further improvements in OuC even with underscore the need to find the most appropriate network structure for a given NCOP problem. An overly complex network would require an excessive amount of computing resources for optimization that may not be necessary. Conversely, a network that is too simple would not find the appropriate solution, regardless of the amount of training. Thus, the size of the neural network should be determined by factors related to the NCOP being solved. These include the number of decision variables, the constraints, and the complexity of both the objective and constraint functions.

Learning Rate. Figure 7(C) shows the performance of neural networks trained with different learning rates. At iteration 1000, the optimal performance is observed at the learning rate of α = 0.0001. It is important to note that if we zoom into the first 100 iterations, a larger learning rate α = 0.001 is more effective. This suggests that the choice of learning rate should depend on the actual preferences of the user.

The advantage of a large learning rate is that it can find a better prediction for the NCOP in a short time, while the disadvantage is that it performs poorly in the long run. In contrast, a small learning rate finds better solutions in the long run. Overall, the choice of learning rate should be determined by the user's specific requirements for speed and accuracy.

L 1 Norm Minimization Problem

Consider the following NCOP problem:

Example 4: 2. In (A), the red circle indicates the first time a feasible solution is found. (A), (B), (C), and (D) show the results of the algorithm iterations 0∼100, 100∼1000, 1000∼3000, and 3000∼10000, respectively.

min x f (x) = ∥x∥ 1 s.t. g i (x) = x 2 10 * (i-1)+1 + x 2 10 * (i-1)+2 + • • • + x 2 10 * (i-1)+10 -20 ≤ 0, i = 1, 2, . . . , 300 h(x) = Ax -b = 0, (23)
where x ∈ R 3000 , A ∈ R 1×3000 , with the first half entries of A being 1 and the rest 3, and b = 16.

IPs Description. We examine five different IPs, which are listed in Based on the five IPs, the OuC values at different algorithm iterations are shown in Table 2 and Figure 8. We observe that:

• As shown in Figure 8-(A), our algorithm quickly finds a feasible solution that satisfies both the inequality and equality constraints. Moreover, once the first feasible solution is found, the subsequent solutions given by the algorithm are all within the feasible set.

• The final solutions given by the algorithm are acceptable. After going through the entire solution process, the OuC values associated with the IPs decrease significantly. For IP1, the OuC decreases from 3000 to 5.81 (100% → 0.2%), and similar results can be found for other IPs. Given the fact that problem (23) has a known lower bound 0 for the objective value. This indicates that the final solutions produced by the proposed algorithm are already very close to the optimal solution.

• The OuC decreasing speed or convergence rate varies under different IP configuration. The convergence rate and final result of IP1∼IP4 are significantly better than those of IP5. This may be because IP5 is the farthest from the optimal solution, and thus requires a larger time range and more model training. Nevertheless, the proposed algorithm still significantly improves the OuC performance of IP5 (100% → 1.6%) with the given experimental setup.

• Most of the decrease in OuC values occurs in the first 1000 iterations. In particular, IP1, IP2, and IP4

reduce the OuC values to about 10 within only 1000 iterations, which is already very close to the final result, demonstrating the efficiency of the algorithm. NCOPs to [0, 10] and the IP y 0 as an all-one vector.

Figure 9 shows the resolution of these 100 NCOPs using our neural network approach, and Table 3 shows the statistical information for these OuC results at different iterations. A clear trend can be seen is that all the OuC values decrease as training progresses. Starting from a mean of 372.79, the OuC value drops to 24.00 after 1000 iterations (a reduction from 100% to 6.4%). This trend of decreasing OuC values is not only observed at the mean, but also consistently observed at the 25%, 50%, 75%, and 90% quantiles.

Of particular note is the impressive magnitude of this reduction. The significant reduction in OuC values indicates that our method efficiently navigates the solution space, making steady progress towards optimality.

This highlights the potential of our approach for tackling a wide variety of NCOPs. In summary, these results provide a strong indication of the effectiveness of our proposed method for solving NCOPs, demonstrating its robustness and efficiency over a wide range of problem datasets.

Conclusion

In this study, we present a deep learning-based methodology for solving NCOPs. The proposed methodology is a fruitful fusion of neurodynamic optimization and PINNs. Methodologically, we have extended the PINN approach to accommodate neurodynamic optimization. In addition, we have developed a novel training algorithm that increases computational efficiency by exploiting the problem structure of NCOPs.

Experimental results have demonstrated the effectiveness of the proposed method on a number of NCOPs.

The computational performance can be further improved by tuning the hyperparameters and refining the training details.

In addition, our results have identified several avenues for future research. Specifically, we recommend investigating better methods for selecting initial points and time ranges, exploring different network architectures and advanced neurodynamic optimization techniques. Further development in these areas will undoubtedly improve the effectiveness and robustness of our approach, making it an important tool for

  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 ].

Figure 2 :

 2 Figure 2: I and Q denote the feasible set of inequality constraints and equality constraints, respectively. Ω = I ∩ Q denotes the general feasible set of the problem. Peq is a projection function that maps x pred onto Q.

2

  ODE system Φ(•) ← the COP 3 Construct a NN model ŷ(t; w).

4

  (x best , ûbest ) ← ŷ(t = T ; w) 5 ϵ best = Evaluate(x best ) 6 while iter ≤ max iter do 7 T ∼ U (t 0 , T ) ; ▷ Sample a batch of times, T. 8 L(T, w) ; ▷ Forward propagation. 9 w = w -∇ w L(T, w) ; ▷ Update w through backward propagation. 10 (x curr , ûcurr ) = ŷ(t = T ; w) ; ▷ Extract the predicted solution from the model.

  IVPs. We model the problem (18) by the ODE system (3) and set the time range as [t 0 , T ] = [0, 10]. We choose three initial points to study, namely [0, 0, 0, 0], [1, 0, -2, 3], and [-1, 1, -1, 1], which result in three IVPs. Based on these three initial points, we construct each of the three proposed neural networks (10) as approximate state solutions to the IVPs.

Figure 3 :

 3 Figure 3: Neural network solutions to problem (18). The three neural networks are initialized with three different initial points (IPs). Each row shows a neural network as a function at selected training iterations.

Figure 4 :

 4 Figure4: Comparison of our proposed method with the numerical integration methods on the OuC metric. The OuC metric is defined in Equation (16).

Figure 6 :

 6 Figure 6: Comparison of our proposed method with PINN and the PINN+Lagaris method on the OuC metric. The experiment is conducted on Problem (20)

Figure 7 :

 7 Figure 7: Hyperparameter study. (A): OuC performance on different numbers of neurons in a single layer neural network with a learning rate of 0.01. (B): OuC performance of different layers in a multilayer neural network with 500 neurons per layer and a learning rate of 0.01. (C): OuC performance at different learning rates in a two-layer neural network with 500 neurons per layer.

Figure 8 :

 8 Figure 8: OuC performance with various IP configurations. The detailed descriptions for the five IPs are given in Table2. In (A), the red circle indicates the first time a feasible solution is found. (A), (B), (C), and (D) show the results of the algorithm iterations 0∼100, 100∼1000, 1000∼3000, and 3000∼10000, respectively.

5. 5 .

 5 NCOP Problem SetProblem Set Description. We construct a set of NCOPs based on Example 4 (23), where each NCOP problem takes the following form:min x f (x) = ∥x∥ 1 s.t. g i (x) = x 2 10 * (i-1)+1 + x 2 10 * (i-1)+2 + • • • + x 2 10 * (i-1)+10 -c (k) ≤ 0, i = 1, 2, . . . , 100 h(x) = A (k) x -b (k) = 0,(24)wherex ∈ R 1000 , A (k) ∈ R 1000 , b (k) ∈ R, and c (k) ∈ R. A (k) , b (k) , c(k) are sampled from uniform distributions U (1, 5), U (10, 20), U (20, 30), respectively. We randomly generate 100 different problem data{ A (k) , b (k) , c (k) } to form 100 different NCOPs. These problem datasets { A (k) , b(k) , c (k) } 100 k=1 can be accessed from the link 1 . Consistent with previous experimental subsections, we set the time range for all

Figure 9 :

 9 Figure 9: OuC values of the 100 NCOPs at different training iterations. Left and right show the results of iterations 0∼100 and 100∼1000, respectively.

  Algorithm 1 summarizes how to use our proposed method to solve the NCOP. First, we need to specify an initial condition y(t 0 ) = y 0 and a time range [t 0 , T ] to construct the IVP. Then, we instantiate Optimal Result Retention (ORR) Mechanism. A key to Algorithm 1 is that we use an ORR mechanism based on the OuC metric (16). Specifically, at each iteration, the algorithm compares the OuC value at the current iteration, denoted as OuC curr , with the best OuC value found so far, denoted as OuC best .(x curr , ûcurr ) and (x best , ûbest ) represent the current prediction and the best prediction found so far, respectively. If OuC curr is less than OuC best , it means that the model found a better prediction in this iteration. The algorithm updates OuC best to equal OuC curr and stores the best prediction as (x best , ûbest ) = (x curr , ûcurr ).

11 xcurr = P eq (x curr ) ; ▷ Project xcurr onto equality constraints. 12 OuC curr = OuC(x curr ) ; ▷ Calculate the OuC value of xcurr . 13 if OuC curr < OuC best then 14 OuC best = OuC curr ; ▷ Update the best OuC value. 15 (x best , ûbest ) = (x curr , ûcurr ) ; ▷ Update the best prediction. 16 end 17 return (x best , ûbest ) 18 end

Pipeline.

the proposed model (

10

), which serves as an approximate state solution for this IVP. The model is trained by performing gradient descent on the batch loss (15) 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.

.18, 0.19, 1.22 , 0. 00] OuC: 43.66 Iter: 10 Prediction: [2.95, -3.52, 2.33, 0.00] OuC: 27.39 Iter: 20 Prediction: [3.24, -3.42, 2.16, 1.22] OuC: 27.39 Iter: 20 Prediction: [2.88, -2.96, 2.03, 2.33] OuC: 27.41 Iter: 0 Prediction: [1.67, 1.08, 0.08, 1.94] OuC: 135.25 Iter: 10 Prediction: [1.20, -1.33, 1.72, 0.00] OuC: 32.06 Iter: 20 Prediction: [2.85, -3.05, 2.09, 2.27] OuC: 27.43 Iter: 0 Prediction: [-1.29, 2.22, 0.58, 1.32] OuC: 69.64 Iter: 10 Prediction: [2.13, -2.77, 2.21, 0.00] OuC: 28.94 NN 2 with IP: [1, 0, -2, 3] NN 3 with IP

		4							4							4
				y1(t;w)
		2		y2(t;w)				2							2
	y(t;w)	0		y3(t;w) y4(t;w)			y(t;w)	0						y(t;w)	0
		2			Iter: 0 Prediction: [0			2							2
		4	0	2	4	6	8	10	4	0	2	4	6	8	10	4	0	2	4	6	8	10
		4			t				4			t				4			t
				y1(t;w)
		2		y2(t;w)				2							2
	y(t;w)	0		y3(t;w) y4(t;w)			y(t;w)	0						y(t;w)	0
		2							2							2
		4	0	2	4	6	8	10	4	0	2	4	6	8	10	4	0	2	4	6	8	10
		4			t				4			t				4			t
				y1(t;w)
		2		y2(t;w)				2							2
	y(t;w)	0		y3(t;w) y4(t;w)			y(t;w)	0						y(t;w)	0
		2							2							2
		4	0	2	4	6	8	10	4	0	2	4	6	8	10	4	0	2	4	6	8	10
					t							t							t

  Evolution of the Predictions. As discussed in Section 4.2, the end state prediction ŷ(t = 10; w) is improved by training on the entire approximate state solution. At training iteration 0, the approximate state solution is far from the true state solution, resulting in a high OuC value of the endpoint prediction. After 20 training iterations, the approximate state solution gets closer to the true state solution, and the endpoint prediction improves significantly. Notably, as these networks are constructed with different initial conditions, they have varying initial OuC values. Nonetheless, our proposed algorithm ensures that all networks converge to a prediction with an OuC value of less than 28 after 20 training iterations.

Table 1

 1

	, the PINN

Table 2 :

 2 Description of IPs and their OuC performance at different algorithm iterations. Columns 2 to 6 describe the IPs and the their initial information, and columns 7 to 10 describe their OuC values at different training iterations.

Table 2 .

 2 The IPs are used to configure Algorithm 1 to solve problem (23). As shown in the table, all the five IPs have large initial objective values, with IP2, IP4, and IP5 failing to satisfy the inequality constraint, i.e., max

	i=1,2,...,300 (g i (x)) ≥ 0,

Table 3 :

 3 Statistical data of OuC values for 100 NCOPs at different training iterations. The table shows the mean, standard deviation (STD), and values at the 25%, 50% (median), 75%, and 90% quantiles of the OuC distribution at each iteration.

https://drive.google.com/drive/folders/1D_3HP-fBp9tew4IgDroIQtgVb8vU-vGO?usp=drive_link