- Email: [email protected]
Neural network models and its application for solving linear and quadratic programming problems
Applied Mathematics and Computation 172 (2006) 305–331 www.elsevier.com/locate/amc
Neural network models and its application for solving linear and quadratic programming problems S. Effati *, A.R. Nazemi Department of Mathematics, Teacher Training University of Sabzevar, Sabzevar, Iran
Abstract In this paper we consider two recurrent neural network model for solving linear and quadratic programming problems. The first model is derived from an unconstraint minimization reformulation of the program. The second model directly is obtained of optimality condition for an optimization problem. By applying the energy function and the duality gap, we will compare the convergence these models. We also explore the existence and the convergence of the trajectory and stability properties for the neural networks models. Finally, in some numerical examples, the effectiveness of the methods is shown. 2005 Elsevier Inc. All rights reserved. Keywords: Neural network; Linear programming; Quadratic programming; Reformulation; Optimality condition; Convergence
*
Corresponding author. E-mail addresses: eff[email protected] (S. Effati), [email protected] (A.R. Nazemi). 0096-3003/$ - see front matter 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.02.005
306
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
1. Introduction Optimization problems arise in a wide variety of scientific and engineering applications including signal processing, system identification, filter design, function approximation, regression analysis, and so on. In many engineering and scientific applications, the real-time solution of optimization problems is widely required. However, traditional algorithms for digital computers may not be efficient since the computing time required for a solution is greatly dependent on the dimension and structure of the problems. One possible and very promising approach to real-time optimization is to apply artificial neural networks. Because of the inherent massive parallelism, the neural network approach can solve optimization problems in running time at the orders of magnitude much faster than those of the most popular optimization algorithms executed on general-purpose digital computers. The introduction of artificial neural networks in optimization was stated in 1980 . Since then, significant research results have been achieved for various optimization. For example, in 1984 Chua and Lin [2], in 1985 and 1986 Hopfield and Tank [11,23] proposed a neural network for solving linear programming problems. In 1987, Kennedy and Chua [12] proposed an improved model that always guaranteed convergence. However, their new model converges to only an approximation of the optimal solution. In 1990, Rodriguez-Vazquez et al. [21] proposed a class of neural networks for solving optimization problems. Based on dual and projection methods [10,13,16,17] Xia et al. [24– 27,29,30] presented several neural network for solving linear and quadratic programming problems. We consider a general convex nonlinear program of the following form: minimize f ðxÞ subject to gðxÞ 6 0;
ð1Þ ð2Þ
hðxÞ ¼ 0;
ð3Þ t
where f : Rn ! R, g ¼ ½g1 ; . . . ; gm : Rn ! Rm is an m-dimensional vector-valt ued continuous function of n variables, and h ¼ ½h1 ; . . . ; hl : Rn ! Rl is a l-dimensional vector valued continuous function of n variables and also the functions f and gks (k = 1,2, . . . , m) are convex on Rn and the functions hps (p = 1,2, . . . , l) are affine on Rn . A vector x is called a feasible solution to nonlinear programming (1)–(3) if and only if x satisfies in the m + l constraints of problem (1)–(3). The feasible solution x is said to be a regular point if the gradients $gi(x), i 2 I = {ijgi(x) = 0} and $hp(x), p = 1, . . . , l are linearly independent, also when f,g and h are continuously differentiable functions, and there is a regular point ^x which is feasible, from [15] it follows that ðx; u; vÞ 2 Rnþmþl is an optimal solution to (1)–(3) if and only if (x, u, v) satisfies the Kuhn–Tucker conditions below
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
u P 0;
gðxÞ 6 0; t
hðxÞ ¼ 0; t
rf ðxÞ þ rgðxÞ u þ rhðxÞ v ¼ 0;
ut gðxÞ ¼ 0;
307
ð4Þ
where $g(x) = ($g1(x), $g2(x), . . . , $gm(x))t, $h(x) = ($h1(x), $h2(x), . . . , $hl (x))t are m · n and l · n Jacobian matrixes, respectively. The vectors u and v are the lagrangian multiplier vectors. The essence of neural network approach for optimization is to establish an energy function (nonnegative) and a dynamic system which is a representation of an artificial neural network. The dynamic system is normally in the form of first order ordinary differential equations. It is expected that for an initial point, the dynamic system will approach its static state (or equilibrium point) which corresponds the solution of the underlying optimization problem. An important requirement is that the energy function decreases monotonically as the dynamic system approaches an equilibrium point. For example, note the following models: Utilizing the penalty method [15], a constrained optimization problem (1)– (3) can be approximated by the following unconstrained optimization problem: ( ) m l X s X 2 2 þ minimize E1 ðxÞ ¼ f ðxÞ þ ½g ðxÞ þ hp ðxÞ ð5Þ 2 k¼1 k p¼1 where s is a positive penalty parameter and gþ k ðxÞ ¼ maxf0; g k ðxÞg. Then gradient method transforms the minimization problem into an associated system of ordinary differential equations dx ¼ Cfrf ðxÞ þ s½rgðxÞgþ ðxÞ þ rhðxÞhðxÞg; ð6Þ dt x 2 Rn ;
ð7Þ
where $g(x) and $h(x) are defined as in (4), C = diag(c1, c2, . . . , cn), and ci > 0 which is to scale the convergence rate of (6) and (7), and t þ þ gþ ðxÞ ¼ ½gþ 1 ðxÞ; g2 ðxÞ; . . . ; gm ðxÞ . The system in (6) and (7) is referred to as Kennedy and Chuas neural networks model [12]. About the global convergence and stability of (6) and (7), by [28] we have that in convex problem (1)–(3) for any fixed s > 0, the neural network of (6) and (7) is Lyapunov stable and globally convergent to a minimizer of E1. Moreover, the minimizer is an approximate solution to the problem (1)–(3) when s is sufficiently large. For the convex optimization problem (1)–(3), by the penalty function methods we can convert it into an unconstrained optimization problem which approximates the optimization problem (1)–(3). For example, we consider the following unconstrained optimization problem: ( ) m l 1X 1X 3 2 þ minimize E2 ðxÞ ¼ f ðxÞ þ s ½g ðxÞ þ h ðxÞ ð8Þ 3 k¼1 k 2 p¼1 p
308
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
then we similarly get the gradient model as follows: dx ¼ frf ðxÞ þ s½rgðxÞ½gþ ðxÞ2 þ rhðxÞhðxÞg; dt
ð9Þ
x 2 Rn ;
ð10Þ 2
2
2
2 t
2
þ þ þ 2 where ½gþ ðxÞ ¼ ½gþ 1 ðxÞ ; g 2 ðxÞ ; . . . ; gm ðxÞ and gk ðxÞ ¼ maxf0; gk ðxÞg. For the convex optimization (1)–(3), about the related stability and convergence of (9) and (10), by [25] we know that for any fixed s > 0, the neural network of (9) and (10) is Lyapunov stable and is at least locally convergent to a minimizer of E2, and the minimizer is an approximate solution to the problem (1)– (3) when s is sufficiently large. From the Kuhn–Tucker conditions we can see that the following function:
E3 ðyÞ ¼
m 1X 1 3 t t 2 ½gþ ðxÞ þ krf ðxÞ þ rgðxÞ u þ rhðxÞ vk2 3 k¼1 k 2
1 1 þ ðut gðxÞÞ2 þ khðxÞk22 2 2
ð11Þ
guarantee that y* = (x*, u*, v*) is a minimizer of E3 if and only if y* is a Kuhn– Tucker point with respect to convexity the problem (1)–(3). Let the gradient model for solving the problem (1)–(3) be dy ¼ CrE3 ðyÞ; dt
ð12Þ
where y 2 H ¼ fðx; u; vÞ 2 Rnþmþl ju P 0g, C = diag(c1, c2, . . . , cn+m+l) and ci > 0. By the result concerning stability and convergence of the gradient model, we know that for the convex problem (1)–(3), the neural network described above is Lyapunov stable and is at least locally convergent to an optimal solution to the problem (1)–(3) (see [25]). This paper is organized as follows. In Section 2, we introduce some mathematical preliminaries on which the development and usage of the proposed method is based. In Section 3, we will state a quadratic programming problem, the corresponding dual and the necessary and sufficient optimality condition. In Section 4, using the Fischer–Burmeister NCP function [4–8], we reformulated the problem as an unconstrained minimization problem and derive a gradient model. So, we will only state results of the existence, the convergence of the trajectory, as well as, some stability results. In Section 5, we proposed neural networks directly by putting necessary and sufficient optimality condition as the equations of dynamic system. We analyze the convergence and stability of the stated model. In Section 6, operating characteristics of the proposed neural networks are demonstrated via some illustrative examples. Finally, Section 7 concludes this paper.
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
309
2. Preliminaries This section provides the necessary mathematical backgrounds which used to study the proposed method and its usage. 2.1. Functions and matrices Definition 1.2.1. A matrix M 2 Rnn is said to • A P–0 matrix if all of its principle minors are nonnegative; • A P matrix if all of its principle minors are positive; • Positive semidefinite if xtMx P 0 for all xð6¼ 0Þ 2 Rn , and positive definite if xtMx > 0 for all xð6¼ 0Þ 2 Rn . Definition 2.2.1. The function F : Rn ! Rn is said to be a • P–0 function if and only if the Jacobian matrix oF is a P–0 matrix for all ox x 2 Rn ; • P-function if and only if the Jacobian matrix oF is a P matrix for all x 2 Rn . ox Definition 3.2.1. A function F : Rn ! Rn is said to be • Lipschitz continuous with constant L on a set Rn if for each pair of points x; y 2 Rn kF ðxÞ F ðyÞk 6 Lkx yk; where k Æ k is the l2-norm in Rn ; • Locally Lipschitz continuous on Rn if each point of Rn has a neighborhood D0 Rn such that the above inequality holds for each pair of points x,y 2 D0.
2.2. Stability in differential equations Consider the following differential equation: x_ ðtÞ ¼ f ðxðtÞÞ;
xðt0 Þ ¼ x0 2 Rn :
ð13Þ
The following classical result on the existence and uniqueness of the solution to (13) hold. Theorem 1.2.2. Assume that f is a continuous mapping from Rn to Rn . Then for arbitrary t0 P 0 and x0 2 Rn there exists a local solution x(t), t 2 [t0, s) to (13)
310
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
for some s > t0. If f is locally Lipschitzian continuous at x0 then the solution is unique, and if f is Lipschitzian continuous in Rn then s can be extended to 1. If a local solution defined on [t0, s) that cannot be extended to a local solution on a larger interval [t0, s) is the maximal interval of existence. An arbitrary local solution has an existence to a maximal one. The maximal interval of existence associated with x0 is often denoted by [t0, s(x0)). Theorem 2.2.2. Assume that f is a continuous mapping from Rn to Rn , if x(t), t 2 [t0, s(x0)) is a maximal solution and s(x0) > 1 then lim kxðtÞk ¼ þ1:
t!sðx0 Þ
A point x 2 Rn is called an equilibrium point of (13) if f(x*) = 0. Now, we state some common definition on stability (see [14]). Definition 3.2.2 (Stability in the sense of Lyapunov). Let x(t) be a solution of (13). An isolated equilibrium point x* is Lyapunov stable if for any x0 = x(t0) and any scalar > 0 there exists a d > 0 so that if kx(t0) x*k < d then kx(t) x*k < for t P t0. Definition 4.2.2 (Asymptotic stability). An isolated equilibrium point x* is said to be asymptotic stable if addition to being a Lyapunov stable it has the property that x(t) ! x* as t ! 1, if kx(t0)x*k < d. Definition 5.2.2 (Lyapunov function). Let X Rn be an open neighborhood of x. A continuously differentiable function g : Rn ! R is said to be a Lyapunov function at the state x (over the set X) for Eq. (13) if ( gðxÞ ¼ 0; gðxÞ > 0 for x 2 X; x 6¼ x; ð14Þ t dgðxðtÞÞ ¼ ½r gðxðtÞÞ f ðxðtÞ 6 0Þ; 8x 2 X: xðtÞ dt A Lyapunov function is often called an energy function for (13).
. Theorem 6.2.2 (i) An isolated equilibrium point x* is Lyapunov stable if there exists a Lyapunov function over some neighborhood X of x*. (ii) An isolated equilibrium point x* is asymptotically stable if there exists a Lyapunov function over some neighborhood X of x* satisfying. dgðxðtÞÞ < 0; dt
8xðtÞ 2 X; xðtÞ 6¼ x :
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
311
3. Neural network models for constrained optimization problem We consider a quadratic program (QP) of the following form: 1 f ðxÞ ¼ xt Qx þ Dtx 2
ð15Þ
subject to gðxÞ ¼ Ax b 6 0;
ð16Þ
hðxÞ ¼ Ex f ¼ 0;
ð17Þ
minimize
where Q is an n · n symmetric positive definite matrix, A is an m · n matrix and b 2 Rn , E is an l · n matrix, x 2 Rn and rank(A, E) = m + l. We define the dual quadratic programming (DQP) problem as follows: 1 f ðu; vÞ ¼ xt Qx þ Dtx þ ut ðAx bÞ þ vt ðEx f Þ 2 subject to gðu; vÞ ¼ Qx þ D þ At u þ Et v ¼ 0;
maximize
ð18Þ ð19Þ ð20Þ
u P 0: nþmþl
By the Kuhn–Tucker conditions, ðx; u; vÞ 2 R is an optimal solution to quadratic programming (15)–(17) if and only if (x, u, v) satisfies u P 0;
Ax b 6 0;
Ex ¼ f ;
Qx þ D þ At u þ Et v ¼ 0;
ut ðAx bÞ ¼ 0:
ð21Þ ð22Þ
Thus the Kuhn–Tucker conditions are principle factor for construct neural network models for QP problems. Remark 1.3. Throughout this paper we assumed that the quadratic programming (15)–(17) has unique optimal solution.
4. The first neural network model The Fischer–Burmeister function / : R2 ! R is given by pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi /ða; bÞ ¼ a2 þ b2 a b:
ð23Þ
Here are some properties of this function. Proposition 1.4 (i) /(a, b) = 0 if and only if a P 0, b P 0, ab = 0 (/ is called a NCP function). (ii) the square of /(a, b) is continuously differentiable. (iii) / is twice continuously differentiable everywhere except at the origin, but it is strongly semismooth at the origin.
312
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
Proof. See [20]. h Proposition 2.4. Finding a solution for quadratic programming can be equivalently reformulated as finding a solution of the following equation: !3 2 n m l P P P uk ak1 þ vp ep1 7 6 / a1 ; qi1 xi þ d 1 þ 7 6 i¼1 k¼1 p¼1 6 !7 7 6 n m l 7 6 6 / a2 ; P q xi þ d 2 þ P uk ak2 þ P vp ep2 7 i2 7 6 7 6 i¼1 k¼1 p¼1 7 6 7 6 . 7 6 .. 6 !7 7 6 n m l 7 6 6 / an ; P q xn þ d n þ P uk akn þ P vp epn 7 in 7 6 i¼1 k¼1 p¼1 7 6 7 6 7 ¼ 0; 6 ð24Þ /ðx; u; vÞ ¼ 6 /ðu1 ; ðb AxÞ1 Þ 7 7 6 7 6 /ðu2 ; ðb AxÞ2 Þ 7 6 7 6 .. 7 6 . 7 6 7 6 7 6 /ðum ; ðb AxÞm Þ 7 6 7 6 /ðb1 ; ðf ExÞ1 Þ 7 6 7 6 7 6 /ðb ; ðf ExÞ Þ 2 2 7 6 7 6 .. 7 6 5 4 . /ðbl ; ðf ExÞl Þ where ðx; u; vÞ 2 Rnþmþl and aks (k = 1,2, . . . , n) and bps (p = 1,2, . . . , l) are positive constants that are not nearly to zero, and (b Ax)k (k = 1,2, . . . , m) and (f Ex)p (p = 1,2, . . . , l) are kth and pth rows of (b Ax) and (f Ex), respectively. Proof. By property (i) of the above proposition and Kuhn–Tucker conditions (21) and (22), it is clear that if there exists ðx; u; vÞ 2 Rnþmþl such that / (x, u, v) = 0, for each vector a > 0, b > 0, then it follows that (x, u, v) is an optimal solution for the problem QP (15)–(17) and DQP (18) and (19). Suppose a = (a1, a2, . . . , an)t, u = (u1, u2, . . . , um)t and b = (b1, b2, . . . , bl)t. Therefore (24) can be written in vector form as follows: 2 3 /ða; Qx þ D þ At u þ Et vÞ 6 7 ð25Þ /ðx; u; vÞ ¼ 4 /ðu; b AxÞ 5 ¼ 0: /ðb; f ExÞ
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
Now, we define F : Rnþmþl ! Rnþmþl as follows: 2 3 Qx þ D þ At u þ Et v 6 7 F ðx; u; vÞ ¼ 4 b Ax 5:
313
ð26Þ
f Ex Thus, solving the problems (15)–(17) and (18)–(20) is equivalent to find a solution for the following equation: /ðyÞ ¼ /ðz; F ðyÞÞ ¼ 0;
ð27Þ
where y ¼ ðx; u; vÞ 2 Rnþmþl , z ¼ ða; u; bÞ 2 Rnþmþl and F is defined in (26). We note that / is locally Lipschitz continuous on Rnþmþl , then by [3] Jacobian matrix o/ is well-defined at any point of Rnþmþl . Thus, F is continuously differenoy tiable for all y 2 Rnþmþl . h Lemma 3.4. The function F that has been defined in (26) is a P-function on Rnþmþl . Proof. We show that the Jacobian matrix oF of F is positive definite matrix on oy Rnþmþl . We note that 2 3 Q At E t oF 6 7 ¼ 4 A 0 0 5: oy E 0 0 Now, let K = (p, q, r)t is arbitrary vector in Rnþmþl that is not zero, such that pt = (p1, p2, . . . , pn), qt = (q1, q2, . . . , qm) and rt = (r1, r2, . . . , rl). Since Q is positive definite, we have Kt
oF K ¼ pt Qp > 0: oy
This proof is complete.
ð28Þ h
We define 1 EðyÞ ¼ k/ðyÞk2 : 2
ð29Þ
We note that E(y) P 0 for all y 2 Rnþmþl and y solves QP (15)–(17) and DQP (18)–(20) if and only if E(y) = 0. Hence solving QP (15)–(17) and DQP (18)– (20) is equivalent to finding the global minimizer of the following unconstrained minimization problem: minimize
ð30Þ
EðyÞ nþmþl
subject to y 2 R
:
ð31Þ
314
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
The first result in the next proposition follows directly from the semismoothness of /, see [18–20] and for proof of second result see [4]. The third result follows directly from Lemma 3.4 and [4]. Proposition 4.4 (i) For any y 2 Rnþmþl , we have k/ðy þ dÞ /ðyÞ Vdk ¼ oðkdkÞfor d ! 0
and
V 2
o/ : oy
ð32Þ
(ii) The function E is continuously differentiable with $E(y) = VT/(y) for an arbitrary element V 2 o/ . oy (iii) Any stationary point y of the optimization problem (30) and (31), that is $E(y) = 0, is a solution to QP (15)–(17) and DQP (18)–(20). Let x(Æ), u(Æ), v(Æ) depend on time variables. The objective function E in (30) is continuously differentiable for all y 2 Rnþmþl . Hence, to use the steepest descent, the new neural network model can be described by the following nonlinear dynamical system: dyðtÞ ¼ srEðyðtÞÞ; dt
s > 0;
yð0Þ ¼ y 0 2 Rnþmþl ;
ð33Þ ð34Þ
where s is to scale the convergence rate of (33) and (34). 4.1. Stability analysis of the first model In this section, first we state the uniqueness, the convergence, and other various properties of the trajectory. Then we focus on a particular case where the equilibrium point is isolated. Finally, we discuss stability to the neural network (33) and (34). 4.1.1. Existence of the trajectory and stability of an isolated equilibrium point The first property of the system (33) and (34) is stated in the following proposition. Proposition 1.4.1.1. y* is an optimal solution to QP (15)–(17) and DQP (18)– (20) if and only if y* is an equilibrium point of the neural network (33,34). Proof. Let y* be the optimal solution for the problems (15)–(20). From Proposition 2.4, it is easy to see that /(y*) = 0, also from Proposition 4.4 (ii),
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
315
$E(y*) = 0. Conversely, since an equilibrium point y* for (33, 34) is called a stationary point of (30) an (31), it is well known [4] that F(y) being a P-function (Lemma 3.4) is a sufficient condition for a stationary point to be solution to QP (15)–(17) and DQP (18)–(20). Let L(y0) denote the level set associated with initial point y0 and be given by Lðy 0 Þ ¼ fy 2 Rnþmþl : EðyÞ 6 Eðy 0 Þg: In the next proposition, we state the existence of solution trajectory of (33) and (34). h Proposition 2.4.1.1 (i) For an arbitrary initial state y0, there exists exactly one maximal solution y(t), t 2 [t0, s(y0)) of Eqs. (33) and (34). (ii) If the level set L(y0) is bounded, then s(y0) = 1. Proof (see [31]). Now, we state the solution trajectory convergence of (33) and (34) by following corollary. h Corollary 3.4.1.1 (i) Let y(t):t 2 [t0, s(y0)) be the unique maximal solution to (33) and (34). If s(y0) = 1 and {y(t)} is bounded, then lim rEðyðtÞÞ ¼ 0:
t!1
(ii) L(y0) is bounded and every accumulation point of the trajectory y(t) is a solution to QP (15)–(17) and DQP (18)–(20). Proof (see [31]). For considering asymptotic stability of the neural network (33) and (34), let y ¼ ðx ; u ; v Þ 2 Rnþmþl be a solution to QP (15)–(17) and DQP (18)–(20). Obviously y* is an equilibrium point of (33) and (34). We assume that there is a neighborhood X Rnþmþl of y* such that rEðy Þ ¼ 0;
rEðyÞ 6¼ 0;
8y 2 X and y 6¼ y :
Theorem 4.4.1.1. If y* is an isolated equilibrium point of (33) and (34), then y* is asymptotically stable for (33) and (34). Proof. First, it is clear that E(y) is nonnegative function over Rnþmþl . Since y* is a solution for QP (15)–(17) and DQP (18)–(20), E(y*) = 0. We claim that for any y 2 X*n{y*}, E(y) > 0. Otherwise if there is a y 2 X*n{y*}, such that E(y) = 0, then /(y) = 0, consequently $E(y) = 0. This contradicting with the isolatedness of y* in X*. Furthermore
316
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
dEðyðtÞÞ t dyðtÞ 2 ¼ ½rEðyðtÞÞ ¼ skrEðyðtÞÞk 6 0; dt dt
8 y 2 Rnþmþl
and dEðyðtÞÞ < 0; dt
8y 2 X
and
yðtÞ 6¼ y :
It follows from Theorem 6.2.2 that y* is asymptotically stable for (33) and (34). h Remark 5.4.1.1. Suppose in the problems (15)–(17) and (18)–(20), Q is positive semidefinite. Then, the function F is defined in (26) is a P–0 function. Similar to previous sections and by [31], we can prove that the neural network in (33) and (34) is asymptotic stable at y*.
5. The second neural network model For construct the second model, we consider the equivalent Kuhn–Tucker conditions as follows: Qx þ D þ At u þ Et v ¼ 0; þ
ð35Þ
ðu þ Ax bÞ ¼ u;
ð36Þ
Ex f ¼ 0:
ð37Þ
Solving problems (15)–(17) and (18)–(20) can be equivalent as finding a solution of the nonlinear equations (35)–(37). Let x(Æ), u(Æ), v(Æ) depend on time variables. Hence, the new neural network model corresponding to (35)–(37) can be described by the following nonlinear dynamical system: dx ¼ jðQx þ D þ At u þ Et vÞ; dt
ð38Þ
du þ ¼ jððu þ Ax bÞ uÞ; dt
ð39Þ
dv ¼ jðEx f Þ: dt
ð40Þ
where j is to scale the convergence rate of (38)–(40). Now define U : Rnþmþr ! Rnþmþr as follows:
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
3 ðQx þ D þ At u þ Et vÞ 6 7 U ðx; u; vÞ ¼ 4 ðu þ Ax bÞþ u 5:
317
2
ð41Þ
Ex f By choosing y ¼ ðx; u; vÞ 2 Rnþmþl , neural network models (38)–(40) can be written as follows: dy ¼ jU ðyÞ; dt
ð42Þ
j > 0;
y 2 Rnþmþl :
ð43Þ
5.1. Stability analysis The first, it is clear that U is continuously differentiable. Thus U is locally Lipschitz continuous in Rnþmþl with positive constant kJ(y)k where J(y) is Jacobian matrix for U(y). So, by Theorem 1.2.2, the solution y(t), t 2 [t0, s) to (42) and (43), for some s > t0 is unique as s ! 1. In this section, we show that the proposed neural network has a good stability performance. The proposed neural network in (42) and (43) has two basic properties below. Proposition 1.5.1. ðx ; u ; v Þ 2 Rnþmþl is an isolated equilibrium point of the neural network (42) and (43), if and only if (x*, u*, v*) is optimal solution to QP (15)–(17) and DQP (18)–(20). Proof. Note that, (x*, u*, v*) is an equilibrium point of (42) and (43) if and only if dxdt ¼ 0, dudt ¼ 0, and dvdt ¼ 0. Equivalently Qx þ D þ At u þ Et v ¼ 0; ðu þ Ax bÞþ u ¼ 0; Ex f ¼ 0: Furthermore, (u* + Ax* b)+ u* = 0 if and only if u P 0;
Ax b 6 0;
ut ðAx bÞ ¼ 0:
ð44Þ
Thus the optimal solution of the QP (15)–(17) and the corresponding dual in (18)–(20) and the equilibrium point of the (42) and (43) are equivalent.
318
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
We assume that there is a neighborhood D Rnþmþl of y* such that U(y*) = 0, U(y) 5 0 for each y 2 D* that y 5 y*, equivalently 8 t t > < Qx þ D þ A u þ E v ¼ 0; þ ð45Þ ðu þ Ax bÞ u ¼ 0; > : Ex f ¼ 0; and for each (x, u, v) 2 D* such that (x, u, v) 5 (x*, u*, v*) 8 t t > < Qx þ D þ A u þ E v 6¼ 0; þ ðu þ Ax bÞ u 6¼ 0; > : Ex f 6¼ 0: Now we state our main result in this section.
ð46Þ
h
Theorem 2.5.1. Let y* is an isolated equilibrium point of (42, 43). If the Jacobian matrix of U(y) is positive definite on Rnþmþl , then the neural network in (42) and (43) is asymptotically stable to y*. Proof. First we define a suitable Lyapunov function. Let E : Rmþnþl ! Rmþnþl is defined for (42) and (43) as follows: 1 EðyÞ ¼ U t U ; 2
ð47Þ
where U = (U1, U2, . . . , Un, Un+1, . . . , Un+m, Un+m+1, . . . , Un+m+l)t and U i ¼ ðQx þ D þ At u þ Et vÞi U nþk ¼ ðuk þ ðAx bÞk Þþ uk U nþmþp ¼ ðf ExÞp
ði ¼ 1; . . . ; nÞ; ðk ¼ 1; 2; . . . ; mÞ;
ðp ¼ 1; 2; . . . ; lÞ:
ð48Þ ð49Þ ð50Þ
We show that E(y) is a Lyapunov function over the set D* for Eqs. (42) and (43). We note that E(y) is nonnegative on Rnþmþl . Since y* is optimal solution for QP (15)–(17) and the DQP (18)–(20), obviously E(y*) = 0. For any y 2 D*n{y*}, we claim that E(y) > 0. Otherwise if there is a y 2 D*n{y*} satisfying E(y) = 0, then U(y) = 0. Hence y is an equilibrium point of (42) and (43) and this is contradicting with the isolatedness of y* in Rnþmþl . Now we check the second condition in Theorem 6.2.2. By (42) and (43), dU oU dy ¼ ¼ J ðyÞU ðyÞ; dt oy dt
ð51Þ
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
319
where J(y) is Jacobian matrix for U(y) and is as follows: 2
At 0
Q 6 J ðyÞ ¼ 4 A
0
E
3 Et 7 0 5: 0
and 2
2Q 6 J þ Jt ¼ 4 0 0
0 0
3 0 7 0 5:
0
0
where Jt is transpose of matrix J. Furthermore, dEðyðtÞÞ ¼ dt
t oU oU ¼ U t J t U þ U t JU ¼ U t ðJ t þ J ÞU : U þ Ut oy oy
ð52Þ
Since Q is positive definite, thus the matrix J + Jt is negative definite. thus dEðyðtÞÞ 6 0: dt
ð53Þ
Consequently, the function E(y) is a Lyapunov function for (42) and (43) over the set D*. Because the isolatedness of y*, we have from (53) that dEðyðtÞÞ < 0; dt
8yðtÞ 2 D and yðtÞ 6¼ y :
ð54Þ
Finally, it follows from Theorem 6.2.2 (ii) that y* is asymptotically stable for (42) and (43) and it completes the proof. h Remark 3.5.1. The neural network in (42) and (43) is simple and more intuitive than that of the first model, but similarly some analysis were given [9,22]), however one example in Section 4 shows that the stability of the system cannot be guaranteed in the case that Q is positive semidefinite.
6. Simulation results In this section, we discuss the simulation results through four examples. The simulation is conducted on Matlab 6.5, the ordinary differential equation solver engaged is ode15s, and the system of difference equations converting (33) and (34) and (42) and (43) are the Euler formula.
320
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
Example 1.6. Consider the following linear programming problem [25]: minimize
f ðxÞ ¼ x1 x2 5 35 x1 x2 6 ; 12 12 5 35 g2 ðxÞ ¼ x1 þ x2 6 ; 2 2 g3 ðxÞ ¼ x1 6 5;
subject to g1 ðxÞ ¼
g4 ðxÞ ¼ x2 6 5: The problem has a unique solution (x*, u*) = (5, 5, 0, 0.4, 0, 0.6)t. We obtain 2 5 3 2 35 3 1 12 12 6 5 6 35 7 0 0 1 1 7 6 2 7 627 ð55Þ Q¼ ; D¼ ; A¼6 7; b ¼ 6 7: 4 1 0 5 455 0 0 1 0 1 5 We use the neural networks in (33) and (34) and (42) and (43) to solve this problem. The first, from Lagrangian duality, one can see that (x*, u*) is an optimal solution to this problem, if and only if (x*, u*) satisfied D þ At u ¼ 0;
ð56Þ
Ax b 6 0;
ð57Þ
ut ðAx bÞ ¼ 0;
u P 0:
ð58Þ
We obtain U ðyÞ ¼
ðD þ At uÞ þ
ðu þ Ax bÞ u
;
0 oU ¼ oy A
At 0
;
ð59Þ
where y ¼ ðx; uÞ 2 R6 . The second model by using (42) and (43) is as follows: dx1 5 5 ¼ j 1 þ u1 þ u2 u3 ; 12 2 dt dx2 ¼ jð1 u1 þ u2 þ u4 Þ; dt þ du1 5 35 ¼ j u1 þ x1 x2 u1 ; 12 12 dt þ du2 5 35 ¼ j u2 þ x1 þ x2 u2 ; 2 2 dt
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
321
du3 ¼ jðu3 x1 5Þþ u3 ; dt du4 þ ¼ jðu4 þ x2 5Þ u4 : dt We consider in Table 1 the simulation results for the different scaling factor j and the different time parameter t: Table 1 shows that under same initial point, all solutions are not of higher accuracy. Fig. 1 shows that the second model is divergent in the continuoustime case. To make a comparison, we solve the above problem by using the neural network in (33) and (34). We define the function /ða; DÞ /ðyÞ ¼ ; ð60Þ /ðu; b AxÞ where a = (a1, a2), u 2 R4 . In this example, we suppose a1 = 3 and a2 = 1. Now by reformulation of the (60), we have
Table 1 Initial point y0 = (10, 9, 1, 3, 5, 5) The second model
Time parameter
x*
u*
j=1 j = 10 j = 102 j = 103
1500 1000 100 10
(4.947, 5.0128) (5.0671, 5.0352) (5.1398, 5.0628) (4.8932, 4.9534)
(0, 0.2652, 0, 0.5594) (0, 0.2675, 0, 0.582) (0, 0.4203, 0, 0.6038) (0, 0.3083, 0, 0.5862)
12
10
10 5 8 6
0
4 –5 2 0
0
(a)
2
4
6
8
10
–10 0
2
4
6
8
10
(b)
Fig. 1. Divergent behaviour of the continuous-time with the initial point x0 = (10, 9) and the rate convergent j = 10 and time parameter t = 10 for the second model.
322
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
minimize
EðyÞ ¼
6 1X /2 ðyÞ 2 i¼1 i
subject to y 2 R6 : Note that F and oF are as follows: oy 0 D þ At u oF ¼ F ðyÞ ¼ oy b Ax A
At 0
:
ð61Þ
The simulation results are listed in the following table. All simulation results listed in Table 2 show that the neural network in (33) and (34) is stable in the Lyapunov sense and is convergent to y*. Note that, in this example only the first model is asymptotically stable. The solutions are listed in Table 1, we can see that the second neural network can provide an approximate to exact solution. However, the accuracy is not the good. From the simulation result of both the models, we can conclude that the first model has the advantages of convergence and high accuracy of solutions. Fig. 2 represents that the first model is convergent in the continuous-time case.
Table 2 Initial point y0 = (10, 9, 1, 3, 5, 5) The first model
Time parameter
x*
u*
j=1 j = 10 j = 102 j = 103
3000 300 30 3
(5, 5) (5, 5) (5, 5) (5, 5)
(0, 0.4, 0, 0.6) (0, 0.4, 0, 0.6) (0, 0.4, 0,0.6) (0, 0.4, 0, 0.6)
10
10
9 5 8 7
0
6 –5 5 4
0
(a)
1
2
3
–10 0
1
2
3
(b)
Fig. 2. Convergent behaviour of the continuous-time with the initial point x0 = (10, 9) and the rate convergent j = 103 and time parameter t = 3 for the first model.
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
323
Example 2.6. Let us consider the convex quadratic optimization problem (15)–(17) where Q¼
2 1 1 2
;
D¼
1
ð62Þ
1
and A and b are the same as in Example 1.6. The primal problem and its dual have a unique solution (x*, u*) = (5, 5, 0, 6, 0, 9) (see [25]). We obtain /ðyÞ ¼
/ða; Qx þ DÞ ; /ðu; b AxÞ
F ðyÞ ¼
Qx þ D þ At u
b Ax
;
Q oF ¼ oy A
At
; 0 ð63Þ
where a = (2, 1) and F is P-function. In this example, all simulation results show that both of the neural network in (33) and (34) and (42) and (43) are stable in the Lyapunov sense and convergent to one solution. Under the same initial point and the same scaling factor, we use to solve this problem of the corresponding discrete models of (33) and (34) and (42) and (43) with step size h = 0.01. Since the actual value of the energy function enables us to estimate directly the quality of the solution, we compare the aspect of corresponding energy function. Hence, suppose the scaling factor j = 1 is chosen in our tests. The stopping criterion in our tests for the first model and the second model is E(y(t)) 6 1010. We summarize the simulation results in Tables 3 and 4. The tf in the following tables represent the final time when the stopping criterion is met.
Table 3 Initial point y0 = (10, 10, 0, 0, 0, 0) The first model The second model
x*
u*
Iteration
tf (s)
E(y(tf))
(5, 5) (5, 5)
(0, 6, 0, 9) (0, 6, 0, 9)
7591 2207
75.91 22.07
9.9692 · 1011 9.9087 · 1011
x*
u*
Iteration
tf (s)
E(y(tf))
(5, 5) (5, 5)
(0, 6, 0, 9) (0, 6, 0, 9)
6982 2012
69.82 20.12
9.9917 · 1011 9.7903 · 1011
Table 4 Initial point y0 = (10, 10, 0, 0, 0, 0) The first model The second model
324
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
It is clear that the second model is quicker than the first model, because the energy function for the second model approach to zero faster than that for the first model. Fig. 3(a) and (b) shows their transient behaviour of the energy function of the corresponding to discrete-time neural networks of (33) and (34) and (42) and (43) with initial point (10, 10) by taking the step size h = 0.01. Furthermore, in the discrete case, by taking the step size h = 0.01, the trajectories of the initial point (10, 10) are shown in Fig. 4(a) and (b), respectively. We can see that solutions are quite accurate.
4000
1000
3000
800 600
2000
400 1000
200 0
0 0
2000
4000
6000
0
500
1000
(a)
1500
2000
(b)
Fig. 3. (a) The transient behaviour of the energy function in the first model with initial point y0 = (10, 10, 0, 0, 0, 0). (b) The transient behaviour of the energy function in the second model with initial point y0 = (10, 10, 0, 0, 0, 0).
15
15
10
10
5
5
0
0
–5
–5
–10 0
20
40
(a)
60
80
–10 0
5
10
15
20
25
(b)
Fig. 4. (a) The transient behaviour of the discrete-time of the first model on t 2 [0, 69.82] with initial point x0 = (10, 10). (b) The transient behaviour of the discrete-time of the second model on t 2 [0,20.12] with initial point x0 = (10, 10).
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
325
Example 3.6. Consider the following quadratic optimization problem [1]: minimize f ðxÞ ¼ x21 þ x22 þ x23 subject to g1 ðxÞ ¼ 2x1 þ x2 5 6 0; g2 ðxÞ ¼ x1 þ x3 2 6 0; g3 ðxÞ ¼ x1 þ 1 6 0; g4 ðxÞ ¼ x2 þ 2 6 0; g5 ðxÞ ¼ x3 6 0; where
2
2
2
6 Q ¼ 40 0
0
0
3
2
7 0 5;
0
2
2 3 0 6 7 D ¼ 4 0 5; 0
2 6 6 1 6 A¼6 6 1 6 4 0 0
1 0 0 1 0
3 0 7 1 7 7 0 7 7; 7 0 5 1
3 5 7 6 6 2 7 7 6 7 b¼6 6 1 7: 7 6 4 2 5 0 2
ð64Þ The exact solution is (1, 2, 0). Since the actual value of the duality gap enables us to estimate directly the quality of the solution, we compare the models aspect of the corresponding squared duality gap. All simulation results in Tables 5 and 6 show that the solution trajectory always converges to the unique solution x* = (1, 2, 0) and the corresponding dual optimal solution is (0, 0, 2, 4, 0). For example, let the scaling factor j = 1 and the stopping criterion is the dual2 ity gap kut ðAx bÞk2 6 1010 and for the first model, we suppose a = (5, 2, 3). In the discrete-time case taking step size h = 0.01, the squared duality gap
Table 5 Initial point y0 = (1,1, 2, 1, 0, 1,0,0) kyðtf Þ y k22 The first model The second model
tf (s)
10
4.3631 · 10 3.5225 · 1011
134.95 17.98
kut ðAx bÞk22
Iteration
11
9.9901 · 10 9.9222 · 1011
13 495 1798
Table 6 Initial point y0 = (5, 5, 1, 1, 0, 1, 0, 2) kyðtf Þ y k22 The first model The second model
10
4.3608 · 10 9.5220 · 1011
tf (s) 142.6 17.98
kut ðAx bÞk22 11
9.9556 · 10 4.5409 · 1011
Iteration 14 260 1798
326
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331 25
15 20 15
10
10 5 5 0
0 0
5000
10000
15000
0
500
1000
1500
2000
(b)
(a)
Fig. 5. (a). The transient behaviour of duality gap of the first model with initial point y0 = (1, 1, 2, 1, 0, 1, 0, 0). (b) The transient behaviour of the duality gap of the second model with initial point y0 = (1, 1, 2, 1, 0, 1, 0, 0).
kut ðAx bÞk22 along the trajectory of the zero initial point is shown in Fig. 5(a) and (b). In comparison, the second model is quicker than the first model, because the duality gap for the second model approach to zero is faster than his the first model. Figs. 6 and 7 show trajectories of the neural network (33) and (34) and (42) and (43) with the six different initial points converges to x*, respectively.
2 x*=(1,2,0)
1 0 –1 –2 3
3
2
2
1
1
0
0 –1
–1
Fig. 6. The transient behaviour in the discrete-time of the neural network (33) and (34) with six different initial point, the step size h = 0.01 and the scaling factor j = 10.
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
327
x*=(1,2,0)
2 1 0 –1 –2 –2 –1
2 1
0
0
1
–1 2
–2
Fig. 7. The transient behaviour in the discrete-time of the neural network (42) and (43) with six different initial point, the step size h = 0.01 and the scaling factor j = 1.
Example 4.6. Let us consider the following optimization problem [15]: minimize f ðxÞ ¼ x21 þ x22 þ x23 2x1 3x2 subject to g1 ðxÞ ¼ x1 6 0; g2 ðxÞ ¼ x2 6 0; g3 ðxÞ ¼ x3 6 0; g4 ðxÞ ¼ x4 6 0; h1 ðxÞ ¼ 2x1 þ x2 þ x3 þ x4 ¼ 7; h2 ðxÞ ¼ x1 þ x2 þ 2x3 þ x4 ¼ 6: The optimal solution is x* = (1.1233, 0.6507, 1.8288, 0.5685) and u* = (0, 0, 0, 0, 1.0548, 2.3562). In this example, we consider the complete form of original quadratic problem (15)–(17). We only apply the corresponding discrete model of the neural network (42, 43) by taking step size h = 0.01 to solve this problem. Figs. 8 and 9 illustrate the values of energy function and squared gap over iterations along the trajectory of the recurrent neural network in this example, respectively. It shows that the squared duality gap decreasea approximately zig-zag while the energy function decreases monotonically. So, Fig. 10 depicts the convergence characteristics of the discrete-time recurrent neural network with three different values of the design parameter h, with respect with the scaling factor and the stopping criterion are j = 1 and both of the energy function 2 and the squared duality gap ðEðyðtÞÞ 6 1010 & kut ðAx bÞ þ vt ðEx f Þk2 6 10 10 Þ (see Table 7). Fig. 11 shows that the trajectories of the neural network (42) and (43) with initial point x0 = (5, 5, 1,1) converges to x*.
328
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331 600
400
200
0 –50 –200 0
500
1000
1500
2000
Fig. 8. The transient behaviour of the energy function in the neural network (42) and (43) with initial point y0 = (5, 5, 1, 1, 0, 0, 0, 2, 2, 1).
400
200
0 –50 –200 0
500
1000
1500
2000
Fig. 9. The transient behaviour of the duality gap in the neural network (42) and (43) with initial point y0 = (5, 5, 1, 1, 0, 0, 0, 2, 2, 1).
600 h=0.003 400
h=0.006 h=0.009
200
0 –50 –50
0
100
200
300
400
500
600
700
800
Fig. 10. The transient behaviour in the neural network (42) and (43) of the energy function with different size steps and initial point y0 = (5, 5, 1, 1, 0, 0, 0, 2, 2, 1).
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
329
Table 7 Initial point y0 = (5, 5, 1, 1, 0, 0, 0, 2, 2, 1) kyðtf Þ y k22 8
1.5847 · 10
tf (s)
kut ðAx bÞ þ vt ðEx f Þk22
E(y(tf))
Iteration
16.25
1.7216 · 1011
9.7216 · 1011
1625
5 4 3 2 1 0 –1 –2 –3 –4 –5 0
2
4
6
8
10
12
14
16
18
Fig. 11. The transient behaviour in the discrete-time of the neural network (42) and (43) with initial point x0 = (5, 5, 1, 1), the step size h = 0.01 and the scaling factor j = 1.
7. Concluding remarks In this paper we have proposed two types of convergent recurrent neural networks. For solving the linear programming problems, the same as in Example 1.6, it has been shown that the accuracy of the second model is lower than that of the first model. For strictly convex quadratic optimization, the same as in Examples 2.6, 3.6 and 4.6, all trajectories converge accurately to their corresponding static states, respectively. Furthermore, in this situation, both of the models are stable in the Lyapunov sense and convergent to unique solution of problems. To make a comparison, we applied the energy function and the squared duality gap. The numerical result has shown that the duality gap and the energy function for second model has faster approach to zero than that of the first model.
330
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
Note that the convergence of stated models is faster when a larger scaling factor is applied. In convergence, the trajectories of variables always ultimately converges to the optimal solution, regardless of whether or not we choose initial point located in the feasible region or not.
References [1] S.K. Agrawal, B.C. Fabien, Optimization of Dynamic System, Kluwer Academic publishers, Netherlands, 1999. [2] L.O. Chua, G.-N. Lin, Nonlinear programming without computation, IEEE Trans. Circuits Syst., CAS 31 (1984) 182–188. [3] L.O. Clarke, Optimization and Nonsmooth Analysis, John Wiley and Sons, New York, 1993. [4] T. De Luca, F. Facchinei, C. Kanzow, A smoothness equation approach to the solution of nonlinear complementarity problems, Math. Program. 75 (1996) 407–439. [5] F. Facchinei, C. Kanzow, A nonsmooth inexact Newton method for the solution of large-scale nonlinear complementarity problems, Math. Program. 76 (1997) 493–512. [6] F. Facchinei, J. Soares, A new merit function for nonlinear complementarity problems and a related algorithms, Siam J. Optim. 7 (1997) 225–247. [7] A. Fischer, A special Newton-type optimization method, Optimization 24 (1992) 269–284. [8] A. Fischer, An NCP-function and its use for the solution of complementarity problems, in: D. Du, L. Qi, R. Womersley (Eds.), Recent Advance in Nonsmooth Optimization, World Scientific Publishers, NJ, 1995, pp. 88–105. [9] T.L. Friesz, D.H. Bernstein, N.J. Mehta, R.L. Tobin, S. Ganjlizadeh, Day-to-day dynamic network disequilibria and idealized traveler information system, Oper.Res. 42 (1994) 1120– 1136. [10] E.M. Gafini, D.P. Bertsekas, Two metric projection methods for constraints optimization, SIAM J. Contr. Optim. 22 (1984) 936–964. [11] J.J. Hopfield, D.W. Tank, Neural computation of decisions in optimization problems, Biol. Cybernet. 52 (1985) 141–152. [12] M.P. Kennedy, L.O. Chua, Neural networks for nonlinear programming, IEEE Trans. Circuits Syst. 35 (1988) 554–562. [13] D. Kinderlehrer, G. Stampcchia, An Introduction to Variational Inequalities and Network Equilibrium Problems, Academic, New York, 1980. [14] J. LaSalle, S. Lefschetz, Stability by Liapunovs Direct Method with Applications, Academic Press, New York, 1961. [15] D.G. Luenberger, Linear and Nonlinear Programming, second ed., Addison-Wesley, Reading, MA, 1973. [16] P. Marcotte, Application of Khobotovs algorithms to variational inequalities and network equilibrium problems, Inform. Syst. Oper. Res. 29 (1991) 258–270. [17] J.S. Pang, A posteriori error bounds for the linearly-constrained variational inequality problem, Math. Oper. Res. 12 (1987) 474–484. [18] J.S. Pang, L. Qi, A nonsmooth equations: motivation and algorithms, SIAM J. Optim. 3 (1993) 443–465. [19] L. Qi, Convergence analysis of some algorithms for solving nonsmooth equations, Math. Oper. Res. 18 (1993) 227–244. [20] L. Qi, J. Sun, A nonsmooth version of Newtons method, Math Program. 58 (1993) 353–368. [21] A. Rodriguez-Vazquez, R. Dominguez-Castro, A. Rueda, J.L. Huertas, E. Sanchez-Sinencio, Nonlinear switched-capacitor neural networks for optimization problems, IEEE Trans. Circuits Syst. 37 (1990) 384–397.
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
331
[22] T.E. Smith, T.L. Friesz, D.H. Bernstein, Z. Suo, in: J.S. Pang, M.C. Ferries (Eds.), A Comparative Analysis of Two Minimum-norm Projective Dynamics and Their Relationship to Variational Inequalities, Complementary and Variational Problems, Siam Press, Philadelphia, 1997, pp. 405–424. [23] D.W. Tank, J.J. Hopfield, Simple neural optimization networks: an A/D converter, signal decision circuit, and a linear programming circuit, IEEE Trans. Circuits Syst. 33 (1986) 533– 541. [24] Y. Xia, J. wang, Neural network for solving linear programming problems with bounded variables, IEEE Trans. Neural Networks 6 (1995) 515–519. [25] Y. Wu, Y. Xia, J. Li, W. Chen, A high-performance neural network for solving linear and quadratic programming problems, IEEE Trans. Neural Networks 7 (1996) 643–651. [26] Y. Xia, A new neural network for solving linear programming problems and its applications, IEEE Trans. Neural Networks 7 (1996) 525–529. [27] Y. Xia, A new neural network for solving linear and quadratic programming problems, IEEE Trans. Neural Networks 7 (1996) 1544–1547. [28] Y. Xia, J. Wang, A general methodology for designing globally convergent optimization neural networks, IEEE Trans. Curcuits Syst. 9 (1998) 1331–1343. [29] Y. Xia, J. Wang, A recurrent neural network for solving linear projection equations, Neural Network 13 (2000) 337–350. [30] Y. Xia, H. Leng, J. Wang, A projection neural network and its application to constrained optimization problems, IEEE Trans. Circuits Syst. 49 (2002) 447–458. [31] L.-Z. Liao, Q. Houduo, Q. Liqun, Solving non-linear complementarity problems with neural networks:a reformulation method approach, This work is supported in part by grant FRG/9798/II-42 of Hong Kong baptist University and the Australian Research Council, 1999.
Neural network models and its application for solving linear and quadratic programming problems S. Effati *, A.R. Nazemi Department of Mathematics, Teacher Training University of Sabzevar, Sabzevar, Iran
Abstract In this paper we consider two recurrent neural network model for solving linear and quadratic programming problems. The first model is derived from an unconstraint minimization reformulation of the program. The second model directly is obtained of optimality condition for an optimization problem. By applying the energy function and the duality gap, we will compare the convergence these models. We also explore the existence and the convergence of the trajectory and stability properties for the neural networks models. Finally, in some numerical examples, the effectiveness of the methods is shown. 2005 Elsevier Inc. All rights reserved. Keywords: Neural network; Linear programming; Quadratic programming; Reformulation; Optimality condition; Convergence
*
Corresponding author. E-mail addresses: eff[email protected] (S. Effati), [email protected] (A.R. Nazemi). 0096-3003/$ - see front matter 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.02.005
306
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
1. Introduction Optimization problems arise in a wide variety of scientific and engineering applications including signal processing, system identification, filter design, function approximation, regression analysis, and so on. In many engineering and scientific applications, the real-time solution of optimization problems is widely required. However, traditional algorithms for digital computers may not be efficient since the computing time required for a solution is greatly dependent on the dimension and structure of the problems. One possible and very promising approach to real-time optimization is to apply artificial neural networks. Because of the inherent massive parallelism, the neural network approach can solve optimization problems in running time at the orders of magnitude much faster than those of the most popular optimization algorithms executed on general-purpose digital computers. The introduction of artificial neural networks in optimization was stated in 1980 . Since then, significant research results have been achieved for various optimization. For example, in 1984 Chua and Lin [2], in 1985 and 1986 Hopfield and Tank [11,23] proposed a neural network for solving linear programming problems. In 1987, Kennedy and Chua [12] proposed an improved model that always guaranteed convergence. However, their new model converges to only an approximation of the optimal solution. In 1990, Rodriguez-Vazquez et al. [21] proposed a class of neural networks for solving optimization problems. Based on dual and projection methods [10,13,16,17] Xia et al. [24– 27,29,30] presented several neural network for solving linear and quadratic programming problems. We consider a general convex nonlinear program of the following form: minimize f ðxÞ subject to gðxÞ 6 0;
ð1Þ ð2Þ
hðxÞ ¼ 0;
ð3Þ t
where f : Rn ! R, g ¼ ½g1 ; . . . ; gm : Rn ! Rm is an m-dimensional vector-valt ued continuous function of n variables, and h ¼ ½h1 ; . . . ; hl : Rn ! Rl is a l-dimensional vector valued continuous function of n variables and also the functions f and gks (k = 1,2, . . . , m) are convex on Rn and the functions hps (p = 1,2, . . . , l) are affine on Rn . A vector x is called a feasible solution to nonlinear programming (1)–(3) if and only if x satisfies in the m + l constraints of problem (1)–(3). The feasible solution x is said to be a regular point if the gradients $gi(x), i 2 I = {ijgi(x) = 0} and $hp(x), p = 1, . . . , l are linearly independent, also when f,g and h are continuously differentiable functions, and there is a regular point ^x which is feasible, from [15] it follows that ðx; u; vÞ 2 Rnþmþl is an optimal solution to (1)–(3) if and only if (x, u, v) satisfies the Kuhn–Tucker conditions below
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
u P 0;
gðxÞ 6 0; t
hðxÞ ¼ 0; t
rf ðxÞ þ rgðxÞ u þ rhðxÞ v ¼ 0;
ut gðxÞ ¼ 0;
307
ð4Þ
where $g(x) = ($g1(x), $g2(x), . . . , $gm(x))t, $h(x) = ($h1(x), $h2(x), . . . , $hl (x))t are m · n and l · n Jacobian matrixes, respectively. The vectors u and v are the lagrangian multiplier vectors. The essence of neural network approach for optimization is to establish an energy function (nonnegative) and a dynamic system which is a representation of an artificial neural network. The dynamic system is normally in the form of first order ordinary differential equations. It is expected that for an initial point, the dynamic system will approach its static state (or equilibrium point) which corresponds the solution of the underlying optimization problem. An important requirement is that the energy function decreases monotonically as the dynamic system approaches an equilibrium point. For example, note the following models: Utilizing the penalty method [15], a constrained optimization problem (1)– (3) can be approximated by the following unconstrained optimization problem: ( ) m l X s X 2 2 þ minimize E1 ðxÞ ¼ f ðxÞ þ ½g ðxÞ þ hp ðxÞ ð5Þ 2 k¼1 k p¼1 where s is a positive penalty parameter and gþ k ðxÞ ¼ maxf0; g k ðxÞg. Then gradient method transforms the minimization problem into an associated system of ordinary differential equations dx ¼ Cfrf ðxÞ þ s½rgðxÞgþ ðxÞ þ rhðxÞhðxÞg; ð6Þ dt x 2 Rn ;
ð7Þ
where $g(x) and $h(x) are defined as in (4), C = diag(c1, c2, . . . , cn), and ci > 0 which is to scale the convergence rate of (6) and (7), and t þ þ gþ ðxÞ ¼ ½gþ 1 ðxÞ; g2 ðxÞ; . . . ; gm ðxÞ . The system in (6) and (7) is referred to as Kennedy and Chuas neural networks model [12]. About the global convergence and stability of (6) and (7), by [28] we have that in convex problem (1)–(3) for any fixed s > 0, the neural network of (6) and (7) is Lyapunov stable and globally convergent to a minimizer of E1. Moreover, the minimizer is an approximate solution to the problem (1)–(3) when s is sufficiently large. For the convex optimization problem (1)–(3), by the penalty function methods we can convert it into an unconstrained optimization problem which approximates the optimization problem (1)–(3). For example, we consider the following unconstrained optimization problem: ( ) m l 1X 1X 3 2 þ minimize E2 ðxÞ ¼ f ðxÞ þ s ½g ðxÞ þ h ðxÞ ð8Þ 3 k¼1 k 2 p¼1 p
308
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
then we similarly get the gradient model as follows: dx ¼ frf ðxÞ þ s½rgðxÞ½gþ ðxÞ2 þ rhðxÞhðxÞg; dt
ð9Þ
x 2 Rn ;
ð10Þ 2
2
2
2 t
2
þ þ þ 2 where ½gþ ðxÞ ¼ ½gþ 1 ðxÞ ; g 2 ðxÞ ; . . . ; gm ðxÞ and gk ðxÞ ¼ maxf0; gk ðxÞg. For the convex optimization (1)–(3), about the related stability and convergence of (9) and (10), by [25] we know that for any fixed s > 0, the neural network of (9) and (10) is Lyapunov stable and is at least locally convergent to a minimizer of E2, and the minimizer is an approximate solution to the problem (1)– (3) when s is sufficiently large. From the Kuhn–Tucker conditions we can see that the following function:
E3 ðyÞ ¼
m 1X 1 3 t t 2 ½gþ ðxÞ þ krf ðxÞ þ rgðxÞ u þ rhðxÞ vk2 3 k¼1 k 2
1 1 þ ðut gðxÞÞ2 þ khðxÞk22 2 2
ð11Þ
guarantee that y* = (x*, u*, v*) is a minimizer of E3 if and only if y* is a Kuhn– Tucker point with respect to convexity the problem (1)–(3). Let the gradient model for solving the problem (1)–(3) be dy ¼ CrE3 ðyÞ; dt
ð12Þ
where y 2 H ¼ fðx; u; vÞ 2 Rnþmþl ju P 0g, C = diag(c1, c2, . . . , cn+m+l) and ci > 0. By the result concerning stability and convergence of the gradient model, we know that for the convex problem (1)–(3), the neural network described above is Lyapunov stable and is at least locally convergent to an optimal solution to the problem (1)–(3) (see [25]). This paper is organized as follows. In Section 2, we introduce some mathematical preliminaries on which the development and usage of the proposed method is based. In Section 3, we will state a quadratic programming problem, the corresponding dual and the necessary and sufficient optimality condition. In Section 4, using the Fischer–Burmeister NCP function [4–8], we reformulated the problem as an unconstrained minimization problem and derive a gradient model. So, we will only state results of the existence, the convergence of the trajectory, as well as, some stability results. In Section 5, we proposed neural networks directly by putting necessary and sufficient optimality condition as the equations of dynamic system. We analyze the convergence and stability of the stated model. In Section 6, operating characteristics of the proposed neural networks are demonstrated via some illustrative examples. Finally, Section 7 concludes this paper.
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
309
2. Preliminaries This section provides the necessary mathematical backgrounds which used to study the proposed method and its usage. 2.1. Functions and matrices Definition 1.2.1. A matrix M 2 Rnn is said to • A P–0 matrix if all of its principle minors are nonnegative; • A P matrix if all of its principle minors are positive; • Positive semidefinite if xtMx P 0 for all xð6¼ 0Þ 2 Rn , and positive definite if xtMx > 0 for all xð6¼ 0Þ 2 Rn . Definition 2.2.1. The function F : Rn ! Rn is said to be a • P–0 function if and only if the Jacobian matrix oF is a P–0 matrix for all ox x 2 Rn ; • P-function if and only if the Jacobian matrix oF is a P matrix for all x 2 Rn . ox Definition 3.2.1. A function F : Rn ! Rn is said to be • Lipschitz continuous with constant L on a set Rn if for each pair of points x; y 2 Rn kF ðxÞ F ðyÞk 6 Lkx yk; where k Æ k is the l2-norm in Rn ; • Locally Lipschitz continuous on Rn if each point of Rn has a neighborhood D0 Rn such that the above inequality holds for each pair of points x,y 2 D0.
2.2. Stability in differential equations Consider the following differential equation: x_ ðtÞ ¼ f ðxðtÞÞ;
xðt0 Þ ¼ x0 2 Rn :
ð13Þ
The following classical result on the existence and uniqueness of the solution to (13) hold. Theorem 1.2.2. Assume that f is a continuous mapping from Rn to Rn . Then for arbitrary t0 P 0 and x0 2 Rn there exists a local solution x(t), t 2 [t0, s) to (13)
310
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
for some s > t0. If f is locally Lipschitzian continuous at x0 then the solution is unique, and if f is Lipschitzian continuous in Rn then s can be extended to 1. If a local solution defined on [t0, s) that cannot be extended to a local solution on a larger interval [t0, s) is the maximal interval of existence. An arbitrary local solution has an existence to a maximal one. The maximal interval of existence associated with x0 is often denoted by [t0, s(x0)). Theorem 2.2.2. Assume that f is a continuous mapping from Rn to Rn , if x(t), t 2 [t0, s(x0)) is a maximal solution and s(x0) > 1 then lim kxðtÞk ¼ þ1:
t!sðx0 Þ
A point x 2 Rn is called an equilibrium point of (13) if f(x*) = 0. Now, we state some common definition on stability (see [14]). Definition 3.2.2 (Stability in the sense of Lyapunov). Let x(t) be a solution of (13). An isolated equilibrium point x* is Lyapunov stable if for any x0 = x(t0) and any scalar > 0 there exists a d > 0 so that if kx(t0) x*k < d then kx(t) x*k < for t P t0. Definition 4.2.2 (Asymptotic stability). An isolated equilibrium point x* is said to be asymptotic stable if addition to being a Lyapunov stable it has the property that x(t) ! x* as t ! 1, if kx(t0)x*k < d. Definition 5.2.2 (Lyapunov function). Let X Rn be an open neighborhood of x. A continuously differentiable function g : Rn ! R is said to be a Lyapunov function at the state x (over the set X) for Eq. (13) if ( gðxÞ ¼ 0; gðxÞ > 0 for x 2 X; x 6¼ x; ð14Þ t dgðxðtÞÞ ¼ ½r gðxðtÞÞ f ðxðtÞ 6 0Þ; 8x 2 X: xðtÞ dt A Lyapunov function is often called an energy function for (13).
. Theorem 6.2.2 (i) An isolated equilibrium point x* is Lyapunov stable if there exists a Lyapunov function over some neighborhood X of x*. (ii) An isolated equilibrium point x* is asymptotically stable if there exists a Lyapunov function over some neighborhood X of x* satisfying. dgðxðtÞÞ < 0; dt
8xðtÞ 2 X; xðtÞ 6¼ x :
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
311
3. Neural network models for constrained optimization problem We consider a quadratic program (QP) of the following form: 1 f ðxÞ ¼ xt Qx þ Dtx 2
ð15Þ
subject to gðxÞ ¼ Ax b 6 0;
ð16Þ
hðxÞ ¼ Ex f ¼ 0;
ð17Þ
minimize
where Q is an n · n symmetric positive definite matrix, A is an m · n matrix and b 2 Rn , E is an l · n matrix, x 2 Rn and rank(A, E) = m + l. We define the dual quadratic programming (DQP) problem as follows: 1 f ðu; vÞ ¼ xt Qx þ Dtx þ ut ðAx bÞ þ vt ðEx f Þ 2 subject to gðu; vÞ ¼ Qx þ D þ At u þ Et v ¼ 0;
maximize
ð18Þ ð19Þ ð20Þ
u P 0: nþmþl
By the Kuhn–Tucker conditions, ðx; u; vÞ 2 R is an optimal solution to quadratic programming (15)–(17) if and only if (x, u, v) satisfies u P 0;
Ax b 6 0;
Ex ¼ f ;
Qx þ D þ At u þ Et v ¼ 0;
ut ðAx bÞ ¼ 0:
ð21Þ ð22Þ
Thus the Kuhn–Tucker conditions are principle factor for construct neural network models for QP problems. Remark 1.3. Throughout this paper we assumed that the quadratic programming (15)–(17) has unique optimal solution.
4. The first neural network model The Fischer–Burmeister function / : R2 ! R is given by pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi /ða; bÞ ¼ a2 þ b2 a b:
ð23Þ
Here are some properties of this function. Proposition 1.4 (i) /(a, b) = 0 if and only if a P 0, b P 0, ab = 0 (/ is called a NCP function). (ii) the square of /(a, b) is continuously differentiable. (iii) / is twice continuously differentiable everywhere except at the origin, but it is strongly semismooth at the origin.
312
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
Proof. See [20]. h Proposition 2.4. Finding a solution for quadratic programming can be equivalently reformulated as finding a solution of the following equation: !3 2 n m l P P P uk ak1 þ vp ep1 7 6 / a1 ; qi1 xi þ d 1 þ 7 6 i¼1 k¼1 p¼1 6 !7 7 6 n m l 7 6 6 / a2 ; P q xi þ d 2 þ P uk ak2 þ P vp ep2 7 i2 7 6 7 6 i¼1 k¼1 p¼1 7 6 7 6 . 7 6 .. 6 !7 7 6 n m l 7 6 6 / an ; P q xn þ d n þ P uk akn þ P vp epn 7 in 7 6 i¼1 k¼1 p¼1 7 6 7 6 7 ¼ 0; 6 ð24Þ /ðx; u; vÞ ¼ 6 /ðu1 ; ðb AxÞ1 Þ 7 7 6 7 6 /ðu2 ; ðb AxÞ2 Þ 7 6 7 6 .. 7 6 . 7 6 7 6 7 6 /ðum ; ðb AxÞm Þ 7 6 7 6 /ðb1 ; ðf ExÞ1 Þ 7 6 7 6 7 6 /ðb ; ðf ExÞ Þ 2 2 7 6 7 6 .. 7 6 5 4 . /ðbl ; ðf ExÞl Þ where ðx; u; vÞ 2 Rnþmþl and aks (k = 1,2, . . . , n) and bps (p = 1,2, . . . , l) are positive constants that are not nearly to zero, and (b Ax)k (k = 1,2, . . . , m) and (f Ex)p (p = 1,2, . . . , l) are kth and pth rows of (b Ax) and (f Ex), respectively. Proof. By property (i) of the above proposition and Kuhn–Tucker conditions (21) and (22), it is clear that if there exists ðx; u; vÞ 2 Rnþmþl such that / (x, u, v) = 0, for each vector a > 0, b > 0, then it follows that (x, u, v) is an optimal solution for the problem QP (15)–(17) and DQP (18) and (19). Suppose a = (a1, a2, . . . , an)t, u = (u1, u2, . . . , um)t and b = (b1, b2, . . . , bl)t. Therefore (24) can be written in vector form as follows: 2 3 /ða; Qx þ D þ At u þ Et vÞ 6 7 ð25Þ /ðx; u; vÞ ¼ 4 /ðu; b AxÞ 5 ¼ 0: /ðb; f ExÞ
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
Now, we define F : Rnþmþl ! Rnþmþl as follows: 2 3 Qx þ D þ At u þ Et v 6 7 F ðx; u; vÞ ¼ 4 b Ax 5:
313
ð26Þ
f Ex Thus, solving the problems (15)–(17) and (18)–(20) is equivalent to find a solution for the following equation: /ðyÞ ¼ /ðz; F ðyÞÞ ¼ 0;
ð27Þ
where y ¼ ðx; u; vÞ 2 Rnþmþl , z ¼ ða; u; bÞ 2 Rnþmþl and F is defined in (26). We note that / is locally Lipschitz continuous on Rnþmþl , then by [3] Jacobian matrix o/ is well-defined at any point of Rnþmþl . Thus, F is continuously differenoy tiable for all y 2 Rnþmþl . h Lemma 3.4. The function F that has been defined in (26) is a P-function on Rnþmþl . Proof. We show that the Jacobian matrix oF of F is positive definite matrix on oy Rnþmþl . We note that 2 3 Q At E t oF 6 7 ¼ 4 A 0 0 5: oy E 0 0 Now, let K = (p, q, r)t is arbitrary vector in Rnþmþl that is not zero, such that pt = (p1, p2, . . . , pn), qt = (q1, q2, . . . , qm) and rt = (r1, r2, . . . , rl). Since Q is positive definite, we have Kt
oF K ¼ pt Qp > 0: oy
This proof is complete.
ð28Þ h
We define 1 EðyÞ ¼ k/ðyÞk2 : 2
ð29Þ
We note that E(y) P 0 for all y 2 Rnþmþl and y solves QP (15)–(17) and DQP (18)–(20) if and only if E(y) = 0. Hence solving QP (15)–(17) and DQP (18)– (20) is equivalent to finding the global minimizer of the following unconstrained minimization problem: minimize
ð30Þ
EðyÞ nþmþl
subject to y 2 R
:
ð31Þ
314
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
The first result in the next proposition follows directly from the semismoothness of /, see [18–20] and for proof of second result see [4]. The third result follows directly from Lemma 3.4 and [4]. Proposition 4.4 (i) For any y 2 Rnþmþl , we have k/ðy þ dÞ /ðyÞ Vdk ¼ oðkdkÞfor d ! 0
and
V 2
o/ : oy
ð32Þ
(ii) The function E is continuously differentiable with $E(y) = VT/(y) for an arbitrary element V 2 o/ . oy (iii) Any stationary point y of the optimization problem (30) and (31), that is $E(y) = 0, is a solution to QP (15)–(17) and DQP (18)–(20). Let x(Æ), u(Æ), v(Æ) depend on time variables. The objective function E in (30) is continuously differentiable for all y 2 Rnþmþl . Hence, to use the steepest descent, the new neural network model can be described by the following nonlinear dynamical system: dyðtÞ ¼ srEðyðtÞÞ; dt
s > 0;
yð0Þ ¼ y 0 2 Rnþmþl ;
ð33Þ ð34Þ
where s is to scale the convergence rate of (33) and (34). 4.1. Stability analysis of the first model In this section, first we state the uniqueness, the convergence, and other various properties of the trajectory. Then we focus on a particular case where the equilibrium point is isolated. Finally, we discuss stability to the neural network (33) and (34). 4.1.1. Existence of the trajectory and stability of an isolated equilibrium point The first property of the system (33) and (34) is stated in the following proposition. Proposition 1.4.1.1. y* is an optimal solution to QP (15)–(17) and DQP (18)– (20) if and only if y* is an equilibrium point of the neural network (33,34). Proof. Let y* be the optimal solution for the problems (15)–(20). From Proposition 2.4, it is easy to see that /(y*) = 0, also from Proposition 4.4 (ii),
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
315
$E(y*) = 0. Conversely, since an equilibrium point y* for (33, 34) is called a stationary point of (30) an (31), it is well known [4] that F(y) being a P-function (Lemma 3.4) is a sufficient condition for a stationary point to be solution to QP (15)–(17) and DQP (18)–(20). Let L(y0) denote the level set associated with initial point y0 and be given by Lðy 0 Þ ¼ fy 2 Rnþmþl : EðyÞ 6 Eðy 0 Þg: In the next proposition, we state the existence of solution trajectory of (33) and (34). h Proposition 2.4.1.1 (i) For an arbitrary initial state y0, there exists exactly one maximal solution y(t), t 2 [t0, s(y0)) of Eqs. (33) and (34). (ii) If the level set L(y0) is bounded, then s(y0) = 1. Proof (see [31]). Now, we state the solution trajectory convergence of (33) and (34) by following corollary. h Corollary 3.4.1.1 (i) Let y(t):t 2 [t0, s(y0)) be the unique maximal solution to (33) and (34). If s(y0) = 1 and {y(t)} is bounded, then lim rEðyðtÞÞ ¼ 0:
t!1
(ii) L(y0) is bounded and every accumulation point of the trajectory y(t) is a solution to QP (15)–(17) and DQP (18)–(20). Proof (see [31]). For considering asymptotic stability of the neural network (33) and (34), let y ¼ ðx ; u ; v Þ 2 Rnþmþl be a solution to QP (15)–(17) and DQP (18)–(20). Obviously y* is an equilibrium point of (33) and (34). We assume that there is a neighborhood X Rnþmþl of y* such that rEðy Þ ¼ 0;
rEðyÞ 6¼ 0;
8y 2 X and y 6¼ y :
Theorem 4.4.1.1. If y* is an isolated equilibrium point of (33) and (34), then y* is asymptotically stable for (33) and (34). Proof. First, it is clear that E(y) is nonnegative function over Rnþmþl . Since y* is a solution for QP (15)–(17) and DQP (18)–(20), E(y*) = 0. We claim that for any y 2 X*n{y*}, E(y) > 0. Otherwise if there is a y 2 X*n{y*}, such that E(y) = 0, then /(y) = 0, consequently $E(y) = 0. This contradicting with the isolatedness of y* in X*. Furthermore
316
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
dEðyðtÞÞ t dyðtÞ 2 ¼ ½rEðyðtÞÞ ¼ skrEðyðtÞÞk 6 0; dt dt
8 y 2 Rnþmþl
and dEðyðtÞÞ < 0; dt
8y 2 X
and
yðtÞ 6¼ y :
It follows from Theorem 6.2.2 that y* is asymptotically stable for (33) and (34). h Remark 5.4.1.1. Suppose in the problems (15)–(17) and (18)–(20), Q is positive semidefinite. Then, the function F is defined in (26) is a P–0 function. Similar to previous sections and by [31], we can prove that the neural network in (33) and (34) is asymptotic stable at y*.
5. The second neural network model For construct the second model, we consider the equivalent Kuhn–Tucker conditions as follows: Qx þ D þ At u þ Et v ¼ 0; þ
ð35Þ
ðu þ Ax bÞ ¼ u;
ð36Þ
Ex f ¼ 0:
ð37Þ
Solving problems (15)–(17) and (18)–(20) can be equivalent as finding a solution of the nonlinear equations (35)–(37). Let x(Æ), u(Æ), v(Æ) depend on time variables. Hence, the new neural network model corresponding to (35)–(37) can be described by the following nonlinear dynamical system: dx ¼ jðQx þ D þ At u þ Et vÞ; dt
ð38Þ
du þ ¼ jððu þ Ax bÞ uÞ; dt
ð39Þ
dv ¼ jðEx f Þ: dt
ð40Þ
where j is to scale the convergence rate of (38)–(40). Now define U : Rnþmþr ! Rnþmþr as follows:
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
3 ðQx þ D þ At u þ Et vÞ 6 7 U ðx; u; vÞ ¼ 4 ðu þ Ax bÞþ u 5:
317
2
ð41Þ
Ex f By choosing y ¼ ðx; u; vÞ 2 Rnþmþl , neural network models (38)–(40) can be written as follows: dy ¼ jU ðyÞ; dt
ð42Þ
j > 0;
y 2 Rnþmþl :
ð43Þ
5.1. Stability analysis The first, it is clear that U is continuously differentiable. Thus U is locally Lipschitz continuous in Rnþmþl with positive constant kJ(y)k where J(y) is Jacobian matrix for U(y). So, by Theorem 1.2.2, the solution y(t), t 2 [t0, s) to (42) and (43), for some s > t0 is unique as s ! 1. In this section, we show that the proposed neural network has a good stability performance. The proposed neural network in (42) and (43) has two basic properties below. Proposition 1.5.1. ðx ; u ; v Þ 2 Rnþmþl is an isolated equilibrium point of the neural network (42) and (43), if and only if (x*, u*, v*) is optimal solution to QP (15)–(17) and DQP (18)–(20). Proof. Note that, (x*, u*, v*) is an equilibrium point of (42) and (43) if and only if dxdt ¼ 0, dudt ¼ 0, and dvdt ¼ 0. Equivalently Qx þ D þ At u þ Et v ¼ 0; ðu þ Ax bÞþ u ¼ 0; Ex f ¼ 0: Furthermore, (u* + Ax* b)+ u* = 0 if and only if u P 0;
Ax b 6 0;
ut ðAx bÞ ¼ 0:
ð44Þ
Thus the optimal solution of the QP (15)–(17) and the corresponding dual in (18)–(20) and the equilibrium point of the (42) and (43) are equivalent.
318
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
We assume that there is a neighborhood D Rnþmþl of y* such that U(y*) = 0, U(y) 5 0 for each y 2 D* that y 5 y*, equivalently 8 t t > < Qx þ D þ A u þ E v ¼ 0; þ ð45Þ ðu þ Ax bÞ u ¼ 0; > : Ex f ¼ 0; and for each (x, u, v) 2 D* such that (x, u, v) 5 (x*, u*, v*) 8 t t > < Qx þ D þ A u þ E v 6¼ 0; þ ðu þ Ax bÞ u 6¼ 0; > : Ex f 6¼ 0: Now we state our main result in this section.
ð46Þ
h
Theorem 2.5.1. Let y* is an isolated equilibrium point of (42, 43). If the Jacobian matrix of U(y) is positive definite on Rnþmþl , then the neural network in (42) and (43) is asymptotically stable to y*. Proof. First we define a suitable Lyapunov function. Let E : Rmþnþl ! Rmþnþl is defined for (42) and (43) as follows: 1 EðyÞ ¼ U t U ; 2
ð47Þ
where U = (U1, U2, . . . , Un, Un+1, . . . , Un+m, Un+m+1, . . . , Un+m+l)t and U i ¼ ðQx þ D þ At u þ Et vÞi U nþk ¼ ðuk þ ðAx bÞk Þþ uk U nþmþp ¼ ðf ExÞp
ði ¼ 1; . . . ; nÞ; ðk ¼ 1; 2; . . . ; mÞ;
ðp ¼ 1; 2; . . . ; lÞ:
ð48Þ ð49Þ ð50Þ
We show that E(y) is a Lyapunov function over the set D* for Eqs. (42) and (43). We note that E(y) is nonnegative on Rnþmþl . Since y* is optimal solution for QP (15)–(17) and the DQP (18)–(20), obviously E(y*) = 0. For any y 2 D*n{y*}, we claim that E(y) > 0. Otherwise if there is a y 2 D*n{y*} satisfying E(y) = 0, then U(y) = 0. Hence y is an equilibrium point of (42) and (43) and this is contradicting with the isolatedness of y* in Rnþmþl . Now we check the second condition in Theorem 6.2.2. By (42) and (43), dU oU dy ¼ ¼ J ðyÞU ðyÞ; dt oy dt
ð51Þ
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
319
where J(y) is Jacobian matrix for U(y) and is as follows: 2
At 0
Q 6 J ðyÞ ¼ 4 A
0
E
3 Et 7 0 5: 0
and 2
2Q 6 J þ Jt ¼ 4 0 0
0 0
3 0 7 0 5:
0
0
where Jt is transpose of matrix J. Furthermore, dEðyðtÞÞ ¼ dt
t oU oU ¼ U t J t U þ U t JU ¼ U t ðJ t þ J ÞU : U þ Ut oy oy
ð52Þ
Since Q is positive definite, thus the matrix J + Jt is negative definite. thus dEðyðtÞÞ 6 0: dt
ð53Þ
Consequently, the function E(y) is a Lyapunov function for (42) and (43) over the set D*. Because the isolatedness of y*, we have from (53) that dEðyðtÞÞ < 0; dt
8yðtÞ 2 D and yðtÞ 6¼ y :
ð54Þ
Finally, it follows from Theorem 6.2.2 (ii) that y* is asymptotically stable for (42) and (43) and it completes the proof. h Remark 3.5.1. The neural network in (42) and (43) is simple and more intuitive than that of the first model, but similarly some analysis were given [9,22]), however one example in Section 4 shows that the stability of the system cannot be guaranteed in the case that Q is positive semidefinite.
6. Simulation results In this section, we discuss the simulation results through four examples. The simulation is conducted on Matlab 6.5, the ordinary differential equation solver engaged is ode15s, and the system of difference equations converting (33) and (34) and (42) and (43) are the Euler formula.
320
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
Example 1.6. Consider the following linear programming problem [25]: minimize
f ðxÞ ¼ x1 x2 5 35 x1 x2 6 ; 12 12 5 35 g2 ðxÞ ¼ x1 þ x2 6 ; 2 2 g3 ðxÞ ¼ x1 6 5;
subject to g1 ðxÞ ¼
g4 ðxÞ ¼ x2 6 5: The problem has a unique solution (x*, u*) = (5, 5, 0, 0.4, 0, 0.6)t. We obtain 2 5 3 2 35 3 1 12 12 6 5 6 35 7 0 0 1 1 7 6 2 7 627 ð55Þ Q¼ ; D¼ ; A¼6 7; b ¼ 6 7: 4 1 0 5 455 0 0 1 0 1 5 We use the neural networks in (33) and (34) and (42) and (43) to solve this problem. The first, from Lagrangian duality, one can see that (x*, u*) is an optimal solution to this problem, if and only if (x*, u*) satisfied D þ At u ¼ 0;
ð56Þ
Ax b 6 0;
ð57Þ
ut ðAx bÞ ¼ 0;
u P 0:
ð58Þ
We obtain U ðyÞ ¼
ðD þ At uÞ þ
ðu þ Ax bÞ u
;
0 oU ¼ oy A
At 0
;
ð59Þ
where y ¼ ðx; uÞ 2 R6 . The second model by using (42) and (43) is as follows: dx1 5 5 ¼ j 1 þ u1 þ u2 u3 ; 12 2 dt dx2 ¼ jð1 u1 þ u2 þ u4 Þ; dt þ du1 5 35 ¼ j u1 þ x1 x2 u1 ; 12 12 dt þ du2 5 35 ¼ j u2 þ x1 þ x2 u2 ; 2 2 dt
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
321
du3 ¼ jðu3 x1 5Þþ u3 ; dt du4 þ ¼ jðu4 þ x2 5Þ u4 : dt We consider in Table 1 the simulation results for the different scaling factor j and the different time parameter t: Table 1 shows that under same initial point, all solutions are not of higher accuracy. Fig. 1 shows that the second model is divergent in the continuoustime case. To make a comparison, we solve the above problem by using the neural network in (33) and (34). We define the function /ða; DÞ /ðyÞ ¼ ; ð60Þ /ðu; b AxÞ where a = (a1, a2), u 2 R4 . In this example, we suppose a1 = 3 and a2 = 1. Now by reformulation of the (60), we have
Table 1 Initial point y0 = (10, 9, 1, 3, 5, 5) The second model
Time parameter
x*
u*
j=1 j = 10 j = 102 j = 103
1500 1000 100 10
(4.947, 5.0128) (5.0671, 5.0352) (5.1398, 5.0628) (4.8932, 4.9534)
(0, 0.2652, 0, 0.5594) (0, 0.2675, 0, 0.582) (0, 0.4203, 0, 0.6038) (0, 0.3083, 0, 0.5862)
12
10
10 5 8 6
0
4 –5 2 0
0
(a)
2
4
6
8
10
–10 0
2
4
6
8
10
(b)
Fig. 1. Divergent behaviour of the continuous-time with the initial point x0 = (10, 9) and the rate convergent j = 10 and time parameter t = 10 for the second model.
322
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
minimize
EðyÞ ¼
6 1X /2 ðyÞ 2 i¼1 i
subject to y 2 R6 : Note that F and oF are as follows: oy 0 D þ At u oF ¼ F ðyÞ ¼ oy b Ax A
At 0
:
ð61Þ
The simulation results are listed in the following table. All simulation results listed in Table 2 show that the neural network in (33) and (34) is stable in the Lyapunov sense and is convergent to y*. Note that, in this example only the first model is asymptotically stable. The solutions are listed in Table 1, we can see that the second neural network can provide an approximate to exact solution. However, the accuracy is not the good. From the simulation result of both the models, we can conclude that the first model has the advantages of convergence and high accuracy of solutions. Fig. 2 represents that the first model is convergent in the continuous-time case.
Table 2 Initial point y0 = (10, 9, 1, 3, 5, 5) The first model
Time parameter
x*
u*
j=1 j = 10 j = 102 j = 103
3000 300 30 3
(5, 5) (5, 5) (5, 5) (5, 5)
(0, 0.4, 0, 0.6) (0, 0.4, 0, 0.6) (0, 0.4, 0,0.6) (0, 0.4, 0, 0.6)
10
10
9 5 8 7
0
6 –5 5 4
0
(a)
1
2
3
–10 0
1
2
3
(b)
Fig. 2. Convergent behaviour of the continuous-time with the initial point x0 = (10, 9) and the rate convergent j = 103 and time parameter t = 3 for the first model.
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
323
Example 2.6. Let us consider the convex quadratic optimization problem (15)–(17) where Q¼
2 1 1 2
;
D¼
1
ð62Þ
1
and A and b are the same as in Example 1.6. The primal problem and its dual have a unique solution (x*, u*) = (5, 5, 0, 6, 0, 9) (see [25]). We obtain /ðyÞ ¼
/ða; Qx þ DÞ ; /ðu; b AxÞ
F ðyÞ ¼
Qx þ D þ At u
b Ax
;
Q oF ¼ oy A
At
; 0 ð63Þ
where a = (2, 1) and F is P-function. In this example, all simulation results show that both of the neural network in (33) and (34) and (42) and (43) are stable in the Lyapunov sense and convergent to one solution. Under the same initial point and the same scaling factor, we use to solve this problem of the corresponding discrete models of (33) and (34) and (42) and (43) with step size h = 0.01. Since the actual value of the energy function enables us to estimate directly the quality of the solution, we compare the aspect of corresponding energy function. Hence, suppose the scaling factor j = 1 is chosen in our tests. The stopping criterion in our tests for the first model and the second model is E(y(t)) 6 1010. We summarize the simulation results in Tables 3 and 4. The tf in the following tables represent the final time when the stopping criterion is met.
Table 3 Initial point y0 = (10, 10, 0, 0, 0, 0) The first model The second model
x*
u*
Iteration
tf (s)
E(y(tf))
(5, 5) (5, 5)
(0, 6, 0, 9) (0, 6, 0, 9)
7591 2207
75.91 22.07
9.9692 · 1011 9.9087 · 1011
x*
u*
Iteration
tf (s)
E(y(tf))
(5, 5) (5, 5)
(0, 6, 0, 9) (0, 6, 0, 9)
6982 2012
69.82 20.12
9.9917 · 1011 9.7903 · 1011
Table 4 Initial point y0 = (10, 10, 0, 0, 0, 0) The first model The second model
324
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
It is clear that the second model is quicker than the first model, because the energy function for the second model approach to zero faster than that for the first model. Fig. 3(a) and (b) shows their transient behaviour of the energy function of the corresponding to discrete-time neural networks of (33) and (34) and (42) and (43) with initial point (10, 10) by taking the step size h = 0.01. Furthermore, in the discrete case, by taking the step size h = 0.01, the trajectories of the initial point (10, 10) are shown in Fig. 4(a) and (b), respectively. We can see that solutions are quite accurate.
4000
1000
3000
800 600
2000
400 1000
200 0
0 0
2000
4000
6000
0
500
1000
(a)
1500
2000
(b)
Fig. 3. (a) The transient behaviour of the energy function in the first model with initial point y0 = (10, 10, 0, 0, 0, 0). (b) The transient behaviour of the energy function in the second model with initial point y0 = (10, 10, 0, 0, 0, 0).
15
15
10
10
5
5
0
0
–5
–5
–10 0
20
40
(a)
60
80
–10 0
5
10
15
20
25
(b)
Fig. 4. (a) The transient behaviour of the discrete-time of the first model on t 2 [0, 69.82] with initial point x0 = (10, 10). (b) The transient behaviour of the discrete-time of the second model on t 2 [0,20.12] with initial point x0 = (10, 10).
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
325
Example 3.6. Consider the following quadratic optimization problem [1]: minimize f ðxÞ ¼ x21 þ x22 þ x23 subject to g1 ðxÞ ¼ 2x1 þ x2 5 6 0; g2 ðxÞ ¼ x1 þ x3 2 6 0; g3 ðxÞ ¼ x1 þ 1 6 0; g4 ðxÞ ¼ x2 þ 2 6 0; g5 ðxÞ ¼ x3 6 0; where
2
2
2
6 Q ¼ 40 0
0
0
3
2
7 0 5;
0
2
2 3 0 6 7 D ¼ 4 0 5; 0
2 6 6 1 6 A¼6 6 1 6 4 0 0
1 0 0 1 0
3 0 7 1 7 7 0 7 7; 7 0 5 1
3 5 7 6 6 2 7 7 6 7 b¼6 6 1 7: 7 6 4 2 5 0 2
ð64Þ The exact solution is (1, 2, 0). Since the actual value of the duality gap enables us to estimate directly the quality of the solution, we compare the models aspect of the corresponding squared duality gap. All simulation results in Tables 5 and 6 show that the solution trajectory always converges to the unique solution x* = (1, 2, 0) and the corresponding dual optimal solution is (0, 0, 2, 4, 0). For example, let the scaling factor j = 1 and the stopping criterion is the dual2 ity gap kut ðAx bÞk2 6 1010 and for the first model, we suppose a = (5, 2, 3). In the discrete-time case taking step size h = 0.01, the squared duality gap
Table 5 Initial point y0 = (1,1, 2, 1, 0, 1,0,0) kyðtf Þ y k22 The first model The second model
tf (s)
10
4.3631 · 10 3.5225 · 1011
134.95 17.98
kut ðAx bÞk22
Iteration
11
9.9901 · 10 9.9222 · 1011
13 495 1798
Table 6 Initial point y0 = (5, 5, 1, 1, 0, 1, 0, 2) kyðtf Þ y k22 The first model The second model
10
4.3608 · 10 9.5220 · 1011
tf (s) 142.6 17.98
kut ðAx bÞk22 11
9.9556 · 10 4.5409 · 1011
Iteration 14 260 1798
326
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331 25
15 20 15
10
10 5 5 0
0 0
5000
10000
15000
0
500
1000
1500
2000
(b)
(a)
Fig. 5. (a). The transient behaviour of duality gap of the first model with initial point y0 = (1, 1, 2, 1, 0, 1, 0, 0). (b) The transient behaviour of the duality gap of the second model with initial point y0 = (1, 1, 2, 1, 0, 1, 0, 0).
kut ðAx bÞk22 along the trajectory of the zero initial point is shown in Fig. 5(a) and (b). In comparison, the second model is quicker than the first model, because the duality gap for the second model approach to zero is faster than his the first model. Figs. 6 and 7 show trajectories of the neural network (33) and (34) and (42) and (43) with the six different initial points converges to x*, respectively.
2 x*=(1,2,0)
1 0 –1 –2 3
3
2
2
1
1
0
0 –1
–1
Fig. 6. The transient behaviour in the discrete-time of the neural network (33) and (34) with six different initial point, the step size h = 0.01 and the scaling factor j = 10.
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
327
x*=(1,2,0)
2 1 0 –1 –2 –2 –1
2 1
0
0
1
–1 2
–2
Fig. 7. The transient behaviour in the discrete-time of the neural network (42) and (43) with six different initial point, the step size h = 0.01 and the scaling factor j = 1.
Example 4.6. Let us consider the following optimization problem [15]: minimize f ðxÞ ¼ x21 þ x22 þ x23 2x1 3x2 subject to g1 ðxÞ ¼ x1 6 0; g2 ðxÞ ¼ x2 6 0; g3 ðxÞ ¼ x3 6 0; g4 ðxÞ ¼ x4 6 0; h1 ðxÞ ¼ 2x1 þ x2 þ x3 þ x4 ¼ 7; h2 ðxÞ ¼ x1 þ x2 þ 2x3 þ x4 ¼ 6: The optimal solution is x* = (1.1233, 0.6507, 1.8288, 0.5685) and u* = (0, 0, 0, 0, 1.0548, 2.3562). In this example, we consider the complete form of original quadratic problem (15)–(17). We only apply the corresponding discrete model of the neural network (42, 43) by taking step size h = 0.01 to solve this problem. Figs. 8 and 9 illustrate the values of energy function and squared gap over iterations along the trajectory of the recurrent neural network in this example, respectively. It shows that the squared duality gap decreasea approximately zig-zag while the energy function decreases monotonically. So, Fig. 10 depicts the convergence characteristics of the discrete-time recurrent neural network with three different values of the design parameter h, with respect with the scaling factor and the stopping criterion are j = 1 and both of the energy function 2 and the squared duality gap ðEðyðtÞÞ 6 1010 & kut ðAx bÞ þ vt ðEx f Þk2 6 10 10 Þ (see Table 7). Fig. 11 shows that the trajectories of the neural network (42) and (43) with initial point x0 = (5, 5, 1,1) converges to x*.
328
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331 600
400
200
0 –50 –200 0
500
1000
1500
2000
Fig. 8. The transient behaviour of the energy function in the neural network (42) and (43) with initial point y0 = (5, 5, 1, 1, 0, 0, 0, 2, 2, 1).
400
200
0 –50 –200 0
500
1000
1500
2000
Fig. 9. The transient behaviour of the duality gap in the neural network (42) and (43) with initial point y0 = (5, 5, 1, 1, 0, 0, 0, 2, 2, 1).
600 h=0.003 400
h=0.006 h=0.009
200
0 –50 –50
0
100
200
300
400
500
600
700
800
Fig. 10. The transient behaviour in the neural network (42) and (43) of the energy function with different size steps and initial point y0 = (5, 5, 1, 1, 0, 0, 0, 2, 2, 1).
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
329
Table 7 Initial point y0 = (5, 5, 1, 1, 0, 0, 0, 2, 2, 1) kyðtf Þ y k22 8
1.5847 · 10
tf (s)
kut ðAx bÞ þ vt ðEx f Þk22
E(y(tf))
Iteration
16.25
1.7216 · 1011
9.7216 · 1011
1625
5 4 3 2 1 0 –1 –2 –3 –4 –5 0
2
4
6
8
10
12
14
16
18
Fig. 11. The transient behaviour in the discrete-time of the neural network (42) and (43) with initial point x0 = (5, 5, 1, 1), the step size h = 0.01 and the scaling factor j = 1.
7. Concluding remarks In this paper we have proposed two types of convergent recurrent neural networks. For solving the linear programming problems, the same as in Example 1.6, it has been shown that the accuracy of the second model is lower than that of the first model. For strictly convex quadratic optimization, the same as in Examples 2.6, 3.6 and 4.6, all trajectories converge accurately to their corresponding static states, respectively. Furthermore, in this situation, both of the models are stable in the Lyapunov sense and convergent to unique solution of problems. To make a comparison, we applied the energy function and the squared duality gap. The numerical result has shown that the duality gap and the energy function for second model has faster approach to zero than that of the first model.
330
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
Note that the convergence of stated models is faster when a larger scaling factor is applied. In convergence, the trajectories of variables always ultimately converges to the optimal solution, regardless of whether or not we choose initial point located in the feasible region or not.
References [1] S.K. Agrawal, B.C. Fabien, Optimization of Dynamic System, Kluwer Academic publishers, Netherlands, 1999. [2] L.O. Chua, G.-N. Lin, Nonlinear programming without computation, IEEE Trans. Circuits Syst., CAS 31 (1984) 182–188. [3] L.O. Clarke, Optimization and Nonsmooth Analysis, John Wiley and Sons, New York, 1993. [4] T. De Luca, F. Facchinei, C. Kanzow, A smoothness equation approach to the solution of nonlinear complementarity problems, Math. Program. 75 (1996) 407–439. [5] F. Facchinei, C. Kanzow, A nonsmooth inexact Newton method for the solution of large-scale nonlinear complementarity problems, Math. Program. 76 (1997) 493–512. [6] F. Facchinei, J. Soares, A new merit function for nonlinear complementarity problems and a related algorithms, Siam J. Optim. 7 (1997) 225–247. [7] A. Fischer, A special Newton-type optimization method, Optimization 24 (1992) 269–284. [8] A. Fischer, An NCP-function and its use for the solution of complementarity problems, in: D. Du, L. Qi, R. Womersley (Eds.), Recent Advance in Nonsmooth Optimization, World Scientific Publishers, NJ, 1995, pp. 88–105. [9] T.L. Friesz, D.H. Bernstein, N.J. Mehta, R.L. Tobin, S. Ganjlizadeh, Day-to-day dynamic network disequilibria and idealized traveler information system, Oper.Res. 42 (1994) 1120– 1136. [10] E.M. Gafini, D.P. Bertsekas, Two metric projection methods for constraints optimization, SIAM J. Contr. Optim. 22 (1984) 936–964. [11] J.J. Hopfield, D.W. Tank, Neural computation of decisions in optimization problems, Biol. Cybernet. 52 (1985) 141–152. [12] M.P. Kennedy, L.O. Chua, Neural networks for nonlinear programming, IEEE Trans. Circuits Syst. 35 (1988) 554–562. [13] D. Kinderlehrer, G. Stampcchia, An Introduction to Variational Inequalities and Network Equilibrium Problems, Academic, New York, 1980. [14] J. LaSalle, S. Lefschetz, Stability by Liapunovs Direct Method with Applications, Academic Press, New York, 1961. [15] D.G. Luenberger, Linear and Nonlinear Programming, second ed., Addison-Wesley, Reading, MA, 1973. [16] P. Marcotte, Application of Khobotovs algorithms to variational inequalities and network equilibrium problems, Inform. Syst. Oper. Res. 29 (1991) 258–270. [17] J.S. Pang, A posteriori error bounds for the linearly-constrained variational inequality problem, Math. Oper. Res. 12 (1987) 474–484. [18] J.S. Pang, L. Qi, A nonsmooth equations: motivation and algorithms, SIAM J. Optim. 3 (1993) 443–465. [19] L. Qi, Convergence analysis of some algorithms for solving nonsmooth equations, Math. Oper. Res. 18 (1993) 227–244. [20] L. Qi, J. Sun, A nonsmooth version of Newtons method, Math Program. 58 (1993) 353–368. [21] A. Rodriguez-Vazquez, R. Dominguez-Castro, A. Rueda, J.L. Huertas, E. Sanchez-Sinencio, Nonlinear switched-capacitor neural networks for optimization problems, IEEE Trans. Circuits Syst. 37 (1990) 384–397.
S. Effati, A.R. Nazemi / Appl. Math. Comput. 172 (2006) 305–331
331
[22] T.E. Smith, T.L. Friesz, D.H. Bernstein, Z. Suo, in: J.S. Pang, M.C. Ferries (Eds.), A Comparative Analysis of Two Minimum-norm Projective Dynamics and Their Relationship to Variational Inequalities, Complementary and Variational Problems, Siam Press, Philadelphia, 1997, pp. 405–424. [23] D.W. Tank, J.J. Hopfield, Simple neural optimization networks: an A/D converter, signal decision circuit, and a linear programming circuit, IEEE Trans. Circuits Syst. 33 (1986) 533– 541. [24] Y. Xia, J. wang, Neural network for solving linear programming problems with bounded variables, IEEE Trans. Neural Networks 6 (1995) 515–519. [25] Y. Wu, Y. Xia, J. Li, W. Chen, A high-performance neural network for solving linear and quadratic programming problems, IEEE Trans. Neural Networks 7 (1996) 643–651. [26] Y. Xia, A new neural network for solving linear programming problems and its applications, IEEE Trans. Neural Networks 7 (1996) 525–529. [27] Y. Xia, A new neural network for solving linear and quadratic programming problems, IEEE Trans. Neural Networks 7 (1996) 1544–1547. [28] Y. Xia, J. Wang, A general methodology for designing globally convergent optimization neural networks, IEEE Trans. Curcuits Syst. 9 (1998) 1331–1343. [29] Y. Xia, J. Wang, A recurrent neural network for solving linear projection equations, Neural Network 13 (2000) 337–350. [30] Y. Xia, H. Leng, J. Wang, A projection neural network and its application to constrained optimization problems, IEEE Trans. Circuits Syst. 49 (2002) 447–458. [31] L.-Z. Liao, Q. Houduo, Q. Liqun, Solving non-linear complementarity problems with neural networks:a reformulation method approach, This work is supported in part by grant FRG/9798/II-42 of Hong Kong baptist University and the Australian Research Council, 1999.