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

Authors: Zhensheng Yu; Zilun Wang; Ke Su
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

In this paper, a double nonmonotone quasi-Newton method is proposed for the nonlinear complementarity problem. By using 3-1 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.


---

## Body

## 1. Introduction

In this paper, we consider the following nonlinear complementarity problem (NCP): findx∈Rn such that(1)x≥0,Fx≥0,xTFx=0.where F:Rn⟶Rn is continuously differentiable and the superscript T denotes the transpose operator. When F is linear, problem (1) reduces to a linear complementarity problem (LCP). Throughout 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–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 Newton-type method can be deployed. The 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–11].Another well-known class of algorithm is the smoothing algorithm. The 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–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. The 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.Smoothing procedure allows one to use successful quasi-Newton approaches, and there are many quasi Newton methods available for the nonlinear complementarity problems based on some smoothing functions [19–26].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=Fx 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. This paper is organized as follows: the piecewise linear NCP functions are introduced in Section 1. The 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.


## 2. Algorithm Analysis

To describe our algorithm, we first give the definitions of NCP function andP0 function. We assume that F:Rn⟶Rn is a continuously differentiable P0-function; if, for all x,y∈Rn with x≠y, there exists an index i such that(2)xi−yiTFix−Fiy≥0,xi≠yi,and we regard a pair a,b∈R2 as an NCP pair if a≥0, b≥0, and aTb=0; a function Φ:R2⟶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:(3)Φa,b=3a−a2b,b≥a>0or3b>−a≥0;3a−b2a,a>b>0or3a>−b≥0;9a+9b,else.Ifa,b≠0,0, then(4)∇Φa,b=3−2aba2b2,b≥a>0,or3b>−a≥0;b2a23−2ba,a>b>0or3a>−b≥0;99,else.We define the 4–1 piecewise linear NCP function (k is any positive integer):(5)Φa,b=k2a,ifb≥ka;2kb−b2a,ifa>bk;2k2a+2kb+b2a,ifa<−bk;k2a+4kb,ifb≤−ka<0.Ifa,b≠0,0, then(6)∇Φa,b=k20,ifb≥ka;b2a22k−2ba,ifa>bk;2k2−b2a22k+2ba,ifa<−bk;k24k,ifb≤−ka<0.DenoteH:R2n⟶R2n,Hx,t=t−FxΦx,t,where t is a sequence in the algorithm and t=Fx holds at the optimal solution to NCP.Hence, the NCP can be written as the following minimization problem:(7)minΨx,t=Hx,t.To get the solution of (7), we introduce the notations as follows: (8)αik,βik=1,1,x,t=0,0;∇Φx,t,otherwise.i=1,2,…,n.Obviously,αik>0andβik>0.Denote the Jacobian matrix ofHxk,tk by Vxk,tk, we get(9)Vxk,tk=−F′xkIdiagαikdiagβik.The identity matrix ofn×n, diagonal matrix whose ith diagonal element is αik, and the diagonal matrix whose ith diagonal element is βik are represented by I, diagαik, and diagβik, respectively.We use the nonmonotone line search to present Broyden-like method. The search directionsd and λ are obtained by calculating a system of smooth equation, and the algorithm is described in detail in Algorithm1.Algorithm 1:
Step 0: initialization.Given initial pointx0,t0∈R2n, μ∈0,1, ξ>0, ξ¯<1, B0=Vx0,t0, k=0.Step 1: ifΨxk,tk=0, then stop. Otherwise, calculate the search directionBkdλ=Fxk−tk−Φxk,tk. (a)By (a), we can obtaindk and λk.Step 2: modified linear search technique.Step 2.1 IfΨxk+dk,tk+λk≤ξΨxk,tk, (b)Ψxk+dk,tk+λk≤ξmax0≤r≤mk−1Ψxk−r,tk−r, (c)tk+λk−Fxk+dk≤ξmax0≤r≤mk−1tk−r−Fxk−r, (d)wherem0=0, 0≤mk≤minmk−1+1,M is a positive constant. Then, letxk+1=xk+dk,tk+1=tk+λk, (e)and go to Step 3; otherwise, go to Step 2.2.Step 2.2: forj=0,1,…,, check the following inequality with μj successivelyΦxk+μjdk,tk+μjλk≤ξmax0≤r≤mk−1Φxk−r,xk−r, (f)sk+μjλk−Fxk+μjdk≤ξmax0≤r≤mk−1tk−r−Fxk−r, (g)Letjk be the smallest nonnegative integer j such that (f) and (g) hold for μj. Set ρk:=μjk, andxk+1=xk+ρkdk,tk+1=tk+ρkλk, (h)and go to Step 3.Step 3: UpdateBk to get Bk+1,Bk+1=Bk+ξkzkTyk−Bkzk/zk2, (i)wherezk=xk+1tk+1−xktk,yk=Hxk+1,tk+1−Hxk,tk. (j)Selectξk to satisfy ξk−1≤ξ¯ and matrix Bk+1 is nonsingular.Step 4: Letk=k+1, go to Step 1.


## 3. 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)
SupposeF:Rn⟶Rnis P0-function and it is continuously differentiable.(b)
On the level set of(10)Lx0,t0=x,t∈R2n|Ψx,t≤Ψx0,t0,where F is Lipschitz continuously differentiable, namely, there exists a constant Lsuch that for all x1,x2∈Rn,(11)Fx1−Fx2≤Lx1−x2.Remark 1.
(see [27]). Fx is P0-function, then F′x is positive semidefinite.Lemma 1.
IfHx0,t0≠0, then B0=V0 is nonsingular.Proof.
AssumeHx0,t0≠0. If V0Tu,v=0 for some u,v∈R2n, where u=u1,u2,…,unT and v=v1,v2,…,vnT, then(12)−F′x0u+Iv=0,(13)diagα0u+diagβ0v=0.
By the definitions ofαi0 and βi0, for all i, αi0>0 and βi0>0. Therefore, diagβ0 is nonsingular. Then(14)v=−diagβ0−1diagα0u.
Substitutev in (12) by (14), and multiply by uT, we have(15)−uTFx0u−uTdiagβ0−1diagα0u=0.
According to the definition ofP0-function, all the principal minor determinants of F′x is nonnegative; hence, F′x is positive semidefinite. And matrix diagβ0−1diagα0 is positive definite. Therefore u=0. Together with (14), it holds that v=0, which implies B0 is nonsingular.Lemma 2.
Assume that Assumption1 holds. Then Φxk,tk⟶0, as k⟶∞.Proof.
For convenience, we defineΦlk=max0≤r≤mk−1Φk−r, where k−mk+1≤lk≤k. When mk+1≤mk+1, we have(16)Φlk+1=max0≤r≤mk+1−1Φk+1−r≤max0≤r≤mkΦk+1−r=maxΦlk,Φk+1=Φlk.
Which meansΦlk is decreasing monotonely, and hence, we have Φlk convergent. Based on (c) of Algorithm1, we have Φlk≤Φllk−1.
Byξ∈0,1, Φlk⟶0k⟶∞ holds, so according to Φk+1≤ξΦlk⟶0, the conclusion holds.Lemma 3.
Assume Assumption1 holds. Then tk−Fxk⟶0 as k⟶∞.Proof.
Definetlk−Fxlk=max0≤r≤mk−1tk−r−Fxk−r, where k−M≤lk≤k. For mk+1≤mk+1, we have(17)tlk+1−Fxlk+1=max0≤r≤mk+1−1tk+1−r−Fxk+1−r≤max0≤r≤mktk+1−r−Fxk+1−r=maxtlk−Fxlk,tk+1−Fxk+1=tlk−Fxlk.
From (17), tlk−Fxlk is decreasing in a monotone way; then tlk−Fxlkis convergent.
According to (g) of Algorithm1, tlk−Fxlk≤ξtllk−10−Fxllk−1. By ξ∈0,1, tlk−Fxlk⟶0k⟶∞ holds. That means tk+1−Fxk+1≤ξtlk−Fxlk⟶0 holds by Algorithm1, so the conclusion is as follows.Lemma 4.
Assume Assumption1 holds. Then dk⟶0, λk⟶0, and Hk⟶0, as k⟶∞.Proof.
We haveΦxk,tk⟶0, tk−Fxk⟶0, as k⟶∞ by Lemma 2 and Lemma 3.
So,Hxk,tk⟶0, as k⟶∞:(18)Bkdkλk=Fxk−tk−Φxk,tk=0.
We know thatBk is nonsingular by Algorithm1. So, dk⟶0, and λk⟶0, as k⟶∞.Theorem 1.
Under the same condition in Lemma4, equation (a) of Algorithm1 has solutions, and the definition of Algorithm1 is well.Proof.
On the one hand, we knowB0 is nonsingular by Lemma 1. And Bk 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 Φxk,tk⟶0 and tk−Fxk⟶0 as k⟶∞ by Lemma 2 and Lemma 3. So Hxk,tk⟶0 as k⟶∞.Lemma 5.
Assume Assumption1 holds, and let xk,tk be generated sequence by Algorithm1; then xk,tk⊂Lx0,t0.Proof.
By induction, fork=0, we have x0,t0∈Lx0,t0.
Assumexk,tk∈Lx0,t0; then we have Ψxk,tk≤Ψx0,t0. By (c) and (d) of Algorithm1, we get(19)Ψxk+1,tk+1=Φxk+1,tk+1+tk+1−Fxk+1≤ξmax0≤r≤mk−1Φk−r+tk−r−Fxk−r=ξmax0≤r≤mk−1Ψxk−r,tk−r≤Ψx0,t0.
So,xk+1,tk+1∈Lx0,t0. Based on the similar analysis, it is easy to see xk,tk⊂Lx0,t0 for all k.Theorem 2.
Assume Assumption1 holds, and xk,tk is generated by Algorithm1; then there exists an accumulation point x∗,t∗ of the sequence xk,tk which is solution of NCP(1).Proof.
From Lemma3 and Lemma 4, we know xk,tk⊂Lx0,t0. By Assumption 1(b), we see that Lx0,t0 is bounded. So, xk,tk has an accumulation point. Suppose there exists a subsequence xk,tkk∈K which has an accumulation point x∗,t∗. We should prove Hx∗,t∗=0.
Supposexk,tkk∈K be an infinite sequence generated by Algorithm1. By construction of the algorithm, we know there are two types of successive iteration. Let K1=k|xk+1=xk+dk,tk+1=tk+λk and K2=k|xk+1=xk+ρkdk,tk+1=tk+ρkλk. We need to prove the conclusion by the following two cases:Case I:K1 is an infinite index set. Let the sequence be xk,tkk∈K1, which satisfy (b) of Algorithm1. Therefore,(20)Ψk1≤ξΨk2≤ξ2Ψk3≤⋯≤ξm−1Ψkm.This suggests thatliminfk⟶∞Hxk,tk=0.Case II:K2 is an infinite index set. Let the sequence be xk,tkk∈K2, which satisfy (f) and (g) of Algorithm 1.
It is known thatΦlk is monotone decreasing and limk⟶∞Φk=0 by Lemma 2 and tlk−Fxlk is monotone decreasing and limk⟶∞tk−Fxk=0 by Lemma 3. Therefore, limk⟶∞Hxk,tk=limk⟶∞Φk+limk⟶∞tk−Fxk=0 as k∈K2.
Therefore, the conclusion is followed.


## 4. 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. Hx,t<10−6 was the stop criterion.The number of iterations, the CPU time in seconds, and the value ofxTFx at the final iteration are are listed in Table 1. x0 in Table 1 means the Initial point where ones i,1 means the i dimension of this problem.


### 4.1. Some Test Problems

Examples1–8 (NCP) are considered.Example 1.
Consider (1), where x∈Rnand Fx=Mx+q with(21)M=4−1−14−1⋱⋱⋱⋱⋱−1−14,q=−1−1⋮−1−1.
We usex0=1,1,…,1Tand t0=10−3,10−3,…,10−3Tas the starting points to text this problem.Example 2.
Consider (1), where x∈R3, and Fx:R3⟶R3 given by(22)Fx=x2x3−x2+x3+1.
Figure1 is the 3D diagram of Example 2. 0,λ,0 is an infinite solution of this problem, where λ∈0,1. The initial points x0, t0 are randomly generated, and these elements are in the interval 0,10.Figure 1
Diagram of Example2.Example 3.
Consider (1), where x∈R7, and Fx:R7⟶R7 given by(23)Fx=2x1−x3+x5+3x6−1x2+2x5+x6−x7−3−x1+2x3+x4+x5+2x6−4x7+1x3+x4+x5−x6−1−x1−2x2−x3−x4+5−3x1−x2−2x3+x4+4x2+4x3−1.5.Example 4.
Consider (1), where x∈R4 and Fx:R4⟶R4 given by(24)Fx=x13−8x2+x23−x3+3x2+x3+2x33−3x4+2x43.Example 5.
(Kojima–Shindo Problem). Consider (1), where x∈R4and Fx:R4⟶R4 given by(25)Fx=3x12+2x1x2+2x22+x3+3x4−62x12+x1+x22+10x3+2x4−23x12+x1x2+2x22+2x3+9x4−9x12+3x22+2x3+3x4−3.
6/2,0,0,0.5 is a degenerate solution, and 1,0,3,0 is a nondegenerate solution.Example 6.
(Modified Mathiesen Problem). Consider (1), where x∈R4 and Fx:R4⟶R4 given by(26)Fx=−x2+x3+x4x1−4.5x3+2.7x4x2+15−x1−0.5x3+0.3x4x3+13−x1.Example 7.
The functionfx is endowed with the component as follows:(27)Fx=f1x,f2x,…,fnxT,fix=exi−1,i=1,2,…,n−1,fnx=exn+xn−1.Example 8.
Consider (1), where x∈Rn and Fx=Mx+q with(28)M=diag1n,2n,…,1,q=−1,−1,…,−1T.
Table1 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 xTfx 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 xTfx value. If the curve is closer to 1, the better the ability to solve the problem (Figure 3).Table 1
Iterations, CPU time, andxTfx for NCP Examples 2–6 between Algorithm1 and FDA.
Problemx0Algorithm 1 with 3–1 piecewiseAlgorithm 1 with 4–1 piecewiseFeasible directions algorithmIterCPU timexTfxIterCPU timexTfxIterCPU timexTfx4.1ones100,120.0008145.49E – 09130.0017681.12E – 0850.0043016.58E – 0720.4738085.69E – 08131.5204581.21E – 0850.257754−8.81E – 14ones4069,1315.640096−1.58E – 111318.5517571.16E – 0858.5569597.20E – 14ones8138,1251.2044584.51E – 101370.1452721.44E – 08553.324527−1.27E – 121,1,1T30.000235−1.02E – 0780.000617-5.41E – 08100.0003544.66E – 074.210,10,10T40.0002885.25E – 08100.0006043.26E – 07130.0005293.76E – 071,5,9T70.000462−2.23E – 0850.0007424.64E – 07190.0014653.50E – 07ones7,1150.001014−2.23E – 08220.0017685.53E – 07260.0043016.58E – 074.310∗ones7,1150.001661−4.62E – 09240.0023661.05E – 08NaNNaNNaN10∗rand7,1160.001608−2.14E – 09200.0012041.20E – 08NaNNaNNaN4.41,1,1,1T140.007181−3.11E – 08230.001696-6.07E – 07>500∗InfNaN10,10,10,10T1450.0192681.80E – 08860.004996-2.26E – 07210.0071376.12E – 074.51,1,1,1T130.00103−5.62E – 09260.0016584.93E – 09360.0011216.43E – 071,2,3,4T230.002497−5.65E – 07300.003223-2.75E – 08650.0097649.61E – 074.61,1,1,1T100.0013511.52E – 0890.0014736.19E – 09160.007463.67E – 074.7ones100,1110.039089−9.29E – 11170.0017688.98E – 08460.0243018.81E – 07ones1024,1201.2301442.26E – 10201.745354-2.68E – 11521.1574437.57E – 07ones4069,11525.532848−2.61E – 122035.6157973.61E – 105550.8069548.14E – 07ones8138,115305.640096−4.09E – 1420370.1261464.62E – 0957415.7200277.56E – 07ones100,140.0020647.41E – 07570.1076851.32E – 07NaNNaNNaN4.8ones1024,160.6413495.40E – 0740722.0023663.10E – 06NaNNaNNaNones4069,1826.7277941.87E – 06>500∗InfNaNNaNNaNNaNones8138,110215.6400965.19E – 06>500∗InfNaNNaNNaNNaNFigure 2
Schematic diagram of the changes ofxTfx with iteration of the three algorithms (the same initial point of ones (n, 1)).Figure 3
Performance profile for Algorithm1 and feasible directions algorithm through Examples 1–8.
(a)(b)(c)


### 4.2. 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. The 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.Considering the competitive economic model of production and investment, supposeH is a price system, in which there are N kinds of commodities, we use RN to express commodity space. For producer i, the set of production is Yi⊆RN. For consumer j, the set of consumption is Zj⊆RN. The number of producers and consumers in the system are l and k, respectively. The total production, total consumption, and initial commodity reserve are represented by Yi, Zj, 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 Zj, Yi, λj∈RN.To describe the model better, we assume the following definitions. In particular,Zj, Yi, and λj are independent of x.Definition 1.
Letzj∈Zj, yi∈Yi, x is the equilibrium price:(1)
For everyi, the maximum profit function is x⋅yi.(2)
For everyj, preference maximum element is zj=zj∈Zj|x⋅zj≤x⋅λj+∑j=1lϕji⋅x⋅yi.(3)
Economic equilibrium is defined as∑i=1lλi+∑i=1lx⋅yi−∑j=1kzj=0.
It can be seen from Definition1 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(29)F=∑i=1lλi+∑i=1lx⋅yi−∑j=1kzj,x≥0,x⋅F=0.
Equation (29) is not only the equilibrium state of free allocation, but also the model of linear complementarity problem. If Zj, Yi, 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(30)PQ=50001/γQ−1/γ,where Q is the total quantity produced, P is the market price, and γ is the elasticity of demand with respect to price. Let qi denote the output of firm i and let the total cost function for firm i be given by(31)fiqi=ciqi+βi1+βiL1/βiqiβi+1/βi,Fiq=fi′qi−p∑j=1nqj−qip′∑j=1nqj,i=1,2,…,n,F=F1q,F2q,…,Fiq.Example 9.
Data is given in Table2.Table 2
Data of Example9.
FirmiciLiβi11051.22851.136514450.95250.8Example 10.
Data is given in Table3.Table 3
Data of Example10.
FirmiciLiβi15101.2023101.0038100.9045100.6051101.5063101.0077100.7084101.1096100.95103100.75


### 4.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. The wheel-rail problem is a typical elastic contact problem.Figure4 [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. The 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.Figure 4
Schematic diagram of the two-dimensional contact problem [30].
(a)(b)(c)Regarding a pointx,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(32)∀x,y∈C,w=0,z≥0,∀x,y∈E,w>0,z=0.If the two-dimensional potential contact area with contact surface is discretized, usersmx×my grid is divided, and let n represent the total number of grids; then(33)u=Tz,z,u∈Rn,T∈Rn×n,and the problem can be changed into a linear complementarity problem LCPq,T; to find a pair w,z∈Rn, the following is satisfied(34)w=Tz+q≥0,z≥0,zTw=0,where the coefficient matrix [32] T is a Toeplitz matrix, satisfying(35)T=T0T−1⋯T2−nT1−nT1T0T−1⋯T2−n⋮T1T0⋱⋮Tn−2⋮⋱⋱T−1Tn−1Tn−2⋯T1T0.Example 11.
The diagonal elementTk of the coefficient matrix T is(36)Tk=21+k−1.2,k≠0;2,k=0.Example 12.
The diagonal elementTk of the coefficient matrix T is(37)Tk=2−k,k=0,1,…,n−1.Example 13.
The diagonal elementTk of the coefficient matrix T is(38)Tk=198+1n,k=0;−0.5,k=1;0.25,k=2;116,k=3;0,else.
Table4 shows the performance of Algorithm1 using different piecewise methods for practical application problems. From Figures 5 to 7, it can be seen that Algorithm1 using 3-1 piecewise has a stronger ability to solve all the above problems than Algorithm1 applying 4-1 piecewise.Table 4
Iterations, CPU time, andxTfx for problem 4.9–4.13.
Problemx0Algorithm 1 with 3–1 piecewiseAlgorithm 1 with 4–1 piecewiseIterCPU timeXTF(X)IterCPU timeXTF(X)4.920∗ones5,1200.005197−1.87E − 06230.00523−1.49E − 0630∗ones5,1180.010009−2.28E − 06210.019381−3.98E − 064.1020∗ones5,1450.0337593.61E − 06430.0503651.03E − 0630∗ones5,1550.050477.62E − 07490.055852−1.80E − 064.11ones100,1130.0523091.87E − 08150.0898981.65E − 08ones1024,1142.080973−3.66E − 09162.6163481.94E − 09ones4069,11460.3673151.05E − 091668.6597958.54E − 09ones8138,114574.2456273.05E − 0920682.960276−5.23E − 09ones100,1140.0585791.96E − 09140.072436−5.59E − 094.12ones1024,1142.376266−1.89E − 08142.348159−4.86E − 09ones4069,11253.195688−7.32E − 091461.3210382.96E − 12ones8138,117519.846528−2.89E − 0814582.9336573.25E − 11ones100,1170.0604491.99E − 08190.08139074.00E − 09ones1024,1172.8602665.62E − 08193.1458438.61E − 094.13ones4069,11774.5321557.27E − 081982.960276−2.45E − 09ones8138,121674.2268192.02E − 0823782.1556325.12E − 09Figure 5
Performance profile on CPU time for Algorithm1 with different piecewise functions.Figure 6
Performance profile on iterations for Algorithm1 with different piecewise functions.Figure 7
Performance profile onxTfx for Algorithm1 with different piecewise functions.


## 4.1. Some Test Problems

Examples1–8 (NCP) are considered.Example 1.
Consider (1), where x∈Rnand Fx=Mx+q with(21)M=4−1−14−1⋱⋱⋱⋱⋱−1−14,q=−1−1⋮−1−1.
We usex0=1,1,…,1Tand t0=10−3,10−3,…,10−3Tas the starting points to text this problem.Example 2.
Consider (1), where x∈R3, and Fx:R3⟶R3 given by(22)Fx=x2x3−x2+x3+1.
Figure1 is the 3D diagram of Example 2. 0,λ,0 is an infinite solution of this problem, where λ∈0,1. The initial points x0, t0 are randomly generated, and these elements are in the interval 0,10.Figure 1
Diagram of Example2.Example 3.
Consider (1), where x∈R7, and Fx:R7⟶R7 given by(23)Fx=2x1−x3+x5+3x6−1x2+2x5+x6−x7−3−x1+2x3+x4+x5+2x6−4x7+1x3+x4+x5−x6−1−x1−2x2−x3−x4+5−3x1−x2−2x3+x4+4x2+4x3−1.5.Example 4.
Consider (1), where x∈R4 and Fx:R4⟶R4 given by(24)Fx=x13−8x2+x23−x3+3x2+x3+2x33−3x4+2x43.Example 5.
(Kojima–Shindo Problem). Consider (1), where x∈R4and Fx:R4⟶R4 given by(25)Fx=3x12+2x1x2+2x22+x3+3x4−62x12+x1+x22+10x3+2x4−23x12+x1x2+2x22+2x3+9x4−9x12+3x22+2x3+3x4−3.
6/2,0,0,0.5 is a degenerate solution, and 1,0,3,0 is a nondegenerate solution.Example 6.
(Modified Mathiesen Problem). Consider (1), where x∈R4 and Fx:R4⟶R4 given by(26)Fx=−x2+x3+x4x1−4.5x3+2.7x4x2+15−x1−0.5x3+0.3x4x3+13−x1.Example 7.
The functionfx is endowed with the component as follows:(27)Fx=f1x,f2x,…,fnxT,fix=exi−1,i=1,2,…,n−1,fnx=exn+xn−1.Example 8.
Consider (1), where x∈Rn and Fx=Mx+q with(28)M=diag1n,2n,…,1,q=−1,−1,…,−1T.
Table1 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 xTfx 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 xTfx value. If the curve is closer to 1, the better the ability to solve the problem (Figure 3).Table 1
Iterations, CPU time, andxTfx for NCP Examples 2–6 between Algorithm1 and FDA.
Problemx0Algorithm 1 with 3–1 piecewiseAlgorithm 1 with 4–1 piecewiseFeasible directions algorithmIterCPU timexTfxIterCPU timexTfxIterCPU timexTfx4.1ones100,120.0008145.49E – 09130.0017681.12E – 0850.0043016.58E – 0720.4738085.69E – 08131.5204581.21E – 0850.257754−8.81E – 14ones4069,1315.640096−1.58E – 111318.5517571.16E – 0858.5569597.20E – 14ones8138,1251.2044584.51E – 101370.1452721.44E – 08553.324527−1.27E – 121,1,1T30.000235−1.02E – 0780.000617-5.41E – 08100.0003544.66E – 074.210,10,10T40.0002885.25E – 08100.0006043.26E – 07130.0005293.76E – 071,5,9T70.000462−2.23E – 0850.0007424.64E – 07190.0014653.50E – 07ones7,1150.001014−2.23E – 08220.0017685.53E – 07260.0043016.58E – 074.310∗ones7,1150.001661−4.62E – 09240.0023661.05E – 08NaNNaNNaN10∗rand7,1160.001608−2.14E – 09200.0012041.20E – 08NaNNaNNaN4.41,1,1,1T140.007181−3.11E – 08230.001696-6.07E – 07>500∗InfNaN10,10,10,10T1450.0192681.80E – 08860.004996-2.26E – 07210.0071376.12E – 074.51,1,1,1T130.00103−5.62E – 09260.0016584.93E – 09360.0011216.43E – 071,2,3,4T230.002497−5.65E – 07300.003223-2.75E – 08650.0097649.61E – 074.61,1,1,1T100.0013511.52E – 0890.0014736.19E – 09160.007463.67E – 074.7ones100,1110.039089−9.29E – 11170.0017688.98E – 08460.0243018.81E – 07ones1024,1201.2301442.26E – 10201.745354-2.68E – 11521.1574437.57E – 07ones4069,11525.532848−2.61E – 122035.6157973.61E – 105550.8069548.14E – 07ones8138,115305.640096−4.09E – 1420370.1261464.62E – 0957415.7200277.56E – 07ones100,140.0020647.41E – 07570.1076851.32E – 07NaNNaNNaN4.8ones1024,160.6413495.40E – 0740722.0023663.10E – 06NaNNaNNaNones4069,1826.7277941.87E – 06>500∗InfNaNNaNNaNNaNones8138,110215.6400965.19E – 06>500∗InfNaNNaNNaNNaNFigure 2
Schematic diagram of the changes ofxTfx with iteration of the three algorithms (the same initial point of ones (n, 1)).Figure 3
Performance profile for Algorithm1 and feasible directions algorithm through Examples 1–8.
(a)(b)(c)


## 4.2. 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. The 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.Considering the competitive economic model of production and investment, supposeH is a price system, in which there are N kinds of commodities, we use RN to express commodity space. For producer i, the set of production is Yi⊆RN. For consumer j, the set of consumption is Zj⊆RN. The number of producers and consumers in the system are l and k, respectively. The total production, total consumption, and initial commodity reserve are represented by Yi, Zj, 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 Zj, Yi, λj∈RN.To describe the model better, we assume the following definitions. In particular,Zj, Yi, and λj are independent of x.Definition 1.
Letzj∈Zj, yi∈Yi, x is the equilibrium price:(1)
For everyi, the maximum profit function is x⋅yi.(2)
For everyj, preference maximum element is zj=zj∈Zj|x⋅zj≤x⋅λj+∑j=1lϕji⋅x⋅yi.(3)
Economic equilibrium is defined as∑i=1lλi+∑i=1lx⋅yi−∑j=1kzj=0.
It can be seen from Definition1 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(29)F=∑i=1lλi+∑i=1lx⋅yi−∑j=1kzj,x≥0,x⋅F=0.
Equation (29) is not only the equilibrium state of free allocation, but also the model of linear complementarity problem. If Zj, Yi, 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(30)PQ=50001/γQ−1/γ,where Q is the total quantity produced, P is the market price, and γ is the elasticity of demand with respect to price. Let qi denote the output of firm i and let the total cost function for firm i be given by(31)fiqi=ciqi+βi1+βiL1/βiqiβi+1/βi,Fiq=fi′qi−p∑j=1nqj−qip′∑j=1nqj,i=1,2,…,n,F=F1q,F2q,…,Fiq.Example 9.
Data is given in Table2.Table 2
Data of Example9.
FirmiciLiβi11051.22851.136514450.95250.8Example 10.
Data is given in Table3.Table 3
Data of Example10.
FirmiciLiβi15101.2023101.0038100.9045100.6051101.5063101.0077100.7084101.1096100.95103100.75


## 4.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. The wheel-rail problem is a typical elastic contact problem.Figure4 [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. The 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.Figure 4
Schematic diagram of the two-dimensional contact problem [30].
(a)(b)(c)Regarding a pointx,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(32)∀x,y∈C,w=0,z≥0,∀x,y∈E,w>0,z=0.If the two-dimensional potential contact area with contact surface is discretized, usersmx×my grid is divided, and let n represent the total number of grids; then(33)u=Tz,z,u∈Rn,T∈Rn×n,and the problem can be changed into a linear complementarity problem LCPq,T; to find a pair w,z∈Rn, the following is satisfied(34)w=Tz+q≥0,z≥0,zTw=0,where the coefficient matrix [32] T is a Toeplitz matrix, satisfying(35)T=T0T−1⋯T2−nT1−nT1T0T−1⋯T2−n⋮T1T0⋱⋮Tn−2⋮⋱⋱T−1Tn−1Tn−2⋯T1T0.Example 11.
The diagonal elementTk of the coefficient matrix T is(36)Tk=21+k−1.2,k≠0;2,k=0.Example 12.
The diagonal elementTk of the coefficient matrix T is(37)Tk=2−k,k=0,1,…,n−1.Example 13.
The diagonal elementTk of the coefficient matrix T is(38)Tk=198+1n,k=0;−0.5,k=1;0.25,k=2;116,k=3;0,else.
Table4 shows the performance of Algorithm1 using different piecewise methods for practical application problems. From Figures 5 to 7, it can be seen that Algorithm1 using 3-1 piecewise has a stronger ability to solve all the above problems than Algorithm1 applying 4-1 piecewise.Table 4
Iterations, CPU time, andxTfx for problem 4.9–4.13.
Problemx0Algorithm 1 with 3–1 piecewiseAlgorithm 1 with 4–1 piecewiseIterCPU timeXTF(X)IterCPU timeXTF(X)4.920∗ones5,1200.005197−1.87E − 06230.00523−1.49E − 0630∗ones5,1180.010009−2.28E − 06210.019381−3.98E − 064.1020∗ones5,1450.0337593.61E − 06430.0503651.03E − 0630∗ones5,1550.050477.62E − 07490.055852−1.80E − 064.11ones100,1130.0523091.87E − 08150.0898981.65E − 08ones1024,1142.080973−3.66E − 09162.6163481.94E − 09ones4069,11460.3673151.05E − 091668.6597958.54E − 09ones8138,114574.2456273.05E − 0920682.960276−5.23E − 09ones100,1140.0585791.96E − 09140.072436−5.59E − 094.12ones1024,1142.376266−1.89E − 08142.348159−4.86E − 09ones4069,11253.195688−7.32E − 091461.3210382.96E − 12ones8138,117519.846528−2.89E − 0814582.9336573.25E − 11ones100,1170.0604491.99E − 08190.08139074.00E − 09ones1024,1172.8602665.62E − 08193.1458438.61E − 094.13ones4069,11774.5321557.27E − 081982.960276−2.45E − 09ones8138,121674.2268192.02E − 0823782.1556325.12E − 09Figure 5
Performance profile on CPU time for Algorithm1 with different piecewise functions.Figure 6
Performance profile on iterations for Algorithm1 with different piecewise functions.Figure 7
Performance profile onxTfx for Algorithm1 with different piecewise functions.


## 5. 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.



---

*Source: 6642725-2020-12-16.xml*