A Double Nonmonotone Quasi-Newton Method for Nonlinear Complementarity Problem Based on Piecewise NCP Functions

In this paper, a double nonmonotone quasi-Newton method is proposed for the nonlinear complementarity problem. By using 31 piecewise and 4-1 piecewise nonlinear complementarity functions, the nonlinear complementarity problem is reformulated into a smooth equation. By a double nonmonotone line search, a smooth Broyden-like algorithm is proposed, where a single solution of a smooth equation at each iteration is required with the reduction in the scale of the calculation. Under suitable conditions, the global convergence of the algorithm is proved, and numerical results with some practical applications are given to show the efficiency of the algorithm.


Introduction
In this paper, we consider the following nonlinear complementarity problem (NCP): find x ∈ R n such that where F: R n ⟶ R n is continuously differentiable and the superscript T denotes the transpose operator. When F is linear, problem (1) reduces to a linear complementarity problem (LCP). roughout this paper, the solution set of problem (1), denoted by X * , is assumed to be nonempty. Nonlinear complementarity problems arisen in many practical applications, for example, the KKT systems of mathematical programming problem, the economic equilibrium, the engineering design problem, can be reformulated into the NCP [1][2][3].
During the past decades, various efficient numerical algorithms are proposed to solve the NCP. One of the most effective methods is to transform the NCP into the semismooth equations (based on nonlinear complementarity function, NCP function) so that the semismooth Newtontype method can be deployed. e most well-known NCP functions are the Fischer-Burmeister function [4] (FB NCP function) and the modified FB NCP function [5]. Sun and Qi [6] proposed several NCP functions, investigated their properties, and provided a numerical comparison between the behavior of different NCP functions. Based on NCP functions, some kinds of algorithm are designed, see, for example, [7][8][9][10][11].
Another well-known class of algorithm is the smoothing algorithm.
e main idea of smoothing algorithm is to reformulate the NCP to smooth equations by introducing the smoothing NCP functions. Some smoothing NCP functions and the corresponding algorithms can be found in [12][13][14][15].
Besides the NCP functions mentioned above, a 3-1 piecewise NCP function was proposed by Liu et al. [16], using it to solve the inequality-constrained nonlinear optimization. e advantage of the 3-1 piecewise lies in the absence of the smoothing parameter. Motivated by the 3-1 piecewise NCP function, Su and Yang [17,18] developed smooth-based Newton algorithms with nonmonotone line search for nonlinear complementarity and generalized nonlinear complementarity problems. Different from the previous methods, the authors introduced independent variable quantities to simplify the algorithm, reducing the amount of calculation without using the smoothing parameter.
In this paper, we will construct a 3-1 piecewise and 4-1 piecewise NCP functions and develop a double nonmonotone quasi Newton method to solve the nonlinear complementarity problems. Based on the piecewise NCP functions, the nonlinear complementarity problem is transformed into the smooth equation. Moreover, we only solve one smooth equation at each iteration. In order to get the better numerical results, a double nonmonotone line search is used by combining with the Broyden-like algorithm. Consequently, the omission of the parameter μ and the single calculation of the Jacobian matrix at each iteration have led to the simplicity and flexibility of this approach. Furthermore, let t � F(x) as an independent variable, which has no relationship with x, ensures the realization of our algorithm easier. Our algorithm is proved to be well-defined and globally convergent under suitable conditions. At the end of the paper, we give numerical results to prove the effectiveness of the algorithm. is paper is organized as follows: the piecewise linear NCP functions are introduced in Section 1. e double nonmonotone line search with quasi-Newton method is given in Section 2. In Section 3, the convergence properties of the algorithm are presented. We give some numerical results in Section 4, and the conclusion is drawn in Section 5.

Algorithm Analysis
To describe our algorithm, we first give the definitions of NCP function and P 0 function. We assume that F: R n ⟶ R n is a continuously differentiable P 0 -function; if, for all x, y ∈ R n with x ≠ y, there exists an index i such that and we regard a pair (a, b) ∈ R 2 as an NCP pair if a ≥ 0, b ≥ 0, and a T b � 0; a function Φ: R 2 ⟶ R is called an NCP function, and we have Φ(a, b) � 0 if and only if (a, b) is a NCP pair.
In what follows, we first introduce the 3-1 piecewise NCP function and then define a 4-1 piecewise NCP function: We define the 4-1 piecewise linear NCP function (k is any positive integer): where t is a sequence in the algorithm and t � F(x) holds at the optimal solution to NCP. Hence, the NCP can be written as the following minimization problem: To get the solution of (7), we introduce the notations as follows: Denote the Jacobian matrix of H(x k , t k ) by V(x k , t k ), we get e identity matrix of n × n, diagonal matrix whose ith diagonal element is α k i , and the diagonal matrix whose ith diagonal element is β k i are represented by I, diag(α k i ), and diag(β k i ), respectively. We use the nonmonotone line search to present Broyden-like method. e search directions d and λ are obtained by calculating a system of smooth equation, and the algorithm is described in detail in Algorithm1.

Convergence Analysis
In this section, the global convergence properties of a Broyden-like algorithm with 3-1 piecewise NCP function are discussed. We give some assumptions to prove the convergence of the algorithm.

Assumption 1
(a) Suppose F: R n ⟶ R n is P 0 -function and it is continuously differentiable. (b) On the level set of where F is Lipschitz continuously differentiable, namely, there exists a constant Lsuch that for all Remark 1 (see [27]).
By the definitions of α 0 Substitute v in (12) by (14), and multiply by u T , we have According to the definition of P 0 -function, all the principal minor determinants of F Together with (14), it holds that v � 0, which implies B 0 is nonsingular. □ Lemma 2. Assume that Assumption 1 holds.
Which means ‖Φ l(k) ‖ is decreasing monotonely, and hence, we have
According to (g) of Algorithm 1, We know that B k is nonsingular by Algorithm1. So, d k ⟶ 0, and λ k ⟶ 0, as k ⟶ ∞. Proof. On the one hand, we know B 0 is nonsingular by Lemma 1. And B k produced by the Broyden-like iteration is nonsingular. Hence equation (a) of Algorithm1 has one and only one solution. On the other hand, we know Proof. From Lemma 3 and Lemma 4, we know (x k , t k ) ⊂ L(x 0 , t 0 ). By Assumption 1(b), we see that L(x 0 , t 0 ) is bounded. So, (x k , t k ) has an accumulation point. Suppose there exists a subsequence (x k , t k ) k∈K which has an accumulation point (x * , t * ). We should prove H(x * , t * ) � 0.
Suppose (x k , t k ) k∈K be an infinite sequence generated by Algorithm1. By construction of the algorithm, we know there are two types of successive iteration. Let We need to prove the conclusion by the following two cases: is suggests that liminf k⟶∞ H(x k , t k ) � 0. Case II: K 2 is an infinite index set. Let the sequence be (x k , t k ) k∈K 2 , which satisfy (f ) and (g) of Algorithm 1.

Numerical Results
In this section, some numerical results are given. We used a personal computer with 4.0 GB memory and Intel(R) Core(TM)i5-5200U CPU @2.20 GHz to perform all experiments. We used Windows 10 as the operating system and Matlab R2018b to write the computer codes. In the whole experiment, the parameters used in Algorithm1 were ξ � 0.9, μ � 0.8, ξ k ≡ 1, M is an integer which is randomly selected from 2 to 5. ‖H(x, t)‖ < 10 − 6 was the stop criterion. e number of iterations, the CPU time in seconds, and the value of x (T) F(x) at the final iteration are are listed in Table 1. x 0 in Table 1 (1), where x ∈ R n and F( We use x 0 � (1, 1, . . . , 1) T and t 0 � (10 − 3 , 10 − 3 , . . . , 10 − 3 ) T as the starting points to text this problem. Figure 1 is the 3D diagram of Example 2. (0, λ, 0) is an infinite solution of this problem, where λ ∈ [0, 1]. e initial points x 0 , t 0 are randomly generated, and these elements are in the interval (0, 10).
Step 1: if Ψ(x k , t k ) � 0, then stop. Otherwise, calculate the search direction By (a), we can obtain d k and λ k .
Step 2: modified linear search technique.
Step 2.1 If and go to Step 3; otherwise, go to Step 2.2.
Step 2.2: for j � 0, 1, . . . ,, check the following inequality with μ j successively Let j k be the smallest nonnegative integer j such that (f ) and (g) hold for μ j . Set ρ k : � μ j k , and x k+1 � x k + ρ k d k , t k+1 � t k + ρ k λ k , (h) and go to Step 3.
Example 6 (Modified Mathiesen Problem). Consider (1), where x ∈ R 4 and F(x): R 4 ⟶ R 4 given by Example 7. e function f(x) is endowed with the component as follows: Table 1 shows the results of Examples 1-8 using 3-1 piecewise, 4-1 piecewise Algorithm1 and feasible direction method, respectively. It can be seen from the table that Algorithm1 applying 3-1 piecewise has a good solution to all the above problems. Algorithm1 applying 4-1 piecewise is slightly insufficient, and the feasible direction method has some difficulties in solving examples above, and some of the examples cannot be solved. Figure 2 shows how the x T f(x) value of the three algorithms decreases as the number of iterations increases in each specific example. We use performance profiles [28]-distribution functions for a performance metric-as a tool for comparing different algorithms. We consider the comprehensive performance of the above three algorithms in terms of CPU time, number of iterations, and x T f(x) value. If the curve is closer to 1, the better the ability to solve the problem (Figure 3).

Nash Equilibrium Problem.
General economic equilibrium [29] means that total supply and total demand are exactly equal in a price system. With the existing productivity and technical conditions, producers get the most profit, while consumers get the most utility when they meet the budget constraints. e theory of general economic equilibrium was first put forward by the French economist Walras. Walras believes that when the whole economy is in equilibrium, the prices of all consumer goods and factors of production will have a certain equilibrium value, and their output and supply will have a certain equilibrium quantity. It is assumed that the whole economic system is a large and complete trading market, and the equilibrium price system means that all commodities are traded in this market, and finally all commodities can be traded.  Mathematical Problems in Engineering Considering the competitive economic model of production and investment, suppose H is a price system, in which there are N kinds of commodities, we use R N to express commodity space. For producer i, the set of production is Y i ⊆R N . For consumer j, the set of consumption is Z j ⊆R N . e number of producers and consumers in the system are l and k, respectively. e total production, total consumption, and initial commodity reserve are represented by Y i , Z j , and λ j , respectively, and the proportion of consumer j in the profit of producer i is represented by ϕ ji . Specially, i � 1, . . . , l; j � 1, . . . , k; and Z j , To describe the model better, we assume the following definitions. In particular, Z j , Y i , and λ j are independent of x.
(1) For every i, the maximum profit function is x · y i .
(2) For every j, preference maximum element is z j � z j ∈ Z j |xt · nz j q ≤ hx · xλ j 7 + C l j�1 ϕ ji · x · y i . (3) Economic equilibrium is defined as l i�1 λ i + l i�1 x · y i − k j�1 z j � 0. It can be seen from Definition 1 that when price system H reaches economic equilibrium, the demands of both producers and consumers are satisfied and then all the commodities of price system H are sold, that is, the commodities are cleared. We define the conditions for clearing the goods as Equation (29) is not only the equilibrium state of free allocation, but also the model of linear complementarity problem. If Z j , Y i , and λ j are related to x, (29) will become a nonlinear complementarity problem (NCP).
Let the inverse demand function for the market be defined by where Q is the total quantity produced, P is the market price, and c is the elasticity of demand with respect to price. Let q i denote the output of firm i and let the total cost function for firm i be given by Example 9. Data is given in Table 2.
Example 10. Data is given in Table 3.

Two-Dimensional Contact Problem.
Under the conditions of nonpenetration and negligible attraction between objects, the elastic contact problem mainly requires the  contact surface and the pressure of the contact surface when two objects are pressed together. e wheel-rail problem is a typical elastic contact problem. Figure 4 [31] shows the geometric structure application of the wheel-rail contact phenomenon, where Figure 4(a) represents the overall geometric structure showing the forward speed V and angular velocity ω of the track when the wheel is rolling. e track is deformed by the wheel pressure Fw and the sleeper pressure Fs1 and Fs2. At the same time, the wheel deforms due to the wheel-rail pressure Fr, and Figures 4(b) and 4(c) represent the undeformed and deformed states, respectively.
Regarding a point (x, y) on the contact surface, if z represents the pressure on the point, u represents the displacement from the dashed line to the solid line along the normal direction, q represents the distance of the dashed line when the point is not deformed, and w represents its shape, the gap between the rear wheel and the track is w � u + q. Assume that C is the contact surface and E is the other external area, the geometric relationship shown in Figure 4 can be abbreviated as If the two-dimensional potential contact area with contact surface is discretized, users mx × my grid is divided, and let n represent the total number of grids; then u � Tz, z, u ∈ R n , T ∈ R n×n , and the problem can be changed into a linear complementarity problem LCP(q, T); to find a pair w, z ∈ R n , the following is satisfied where the coefficient matrix [32] T is a Toeplitz matrix, satisfying  T � Example 11. e diagonal element T k of the coefficient matrix T is Example 12. e diagonal element T k of the coefficient matrix T is  Table 4 shows the performance of Algorithm1 using different piecewise methods for practical application problems. From Figures 5 to 7, it can be seen that Algo-rithm1 using 3-1 piecewise has a stronger ability to solve all the above problems than Algorithm1 applying 4-1 piecewise.   Figure 4: Schematic diagram of the two-dimensional contact problem [30].

Conclusion
In this paper, by using 3-1 and 4-1 piecewise nonlinear complementarity problem functions, we reformulate the nonlinear complementarity problem into smooth equations. By using a new nonmonotone line search, a modified smooth Broyden-like algorithm is proposed and the global convergence of the proposed algorithm is obtained, and the numerical tests for some practical problems show the efficiency of the algorithm. How to get the local convergence under certain conditions is worth studying in the future.

Data Availability
e data used to support the findings of this study are included within the article.