2500
IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. 26, NO. 10, OCTOBER 2015
Solving Nonlinear Equality Constrained Multiobjective Optimization Problems Using Neural Networks Mohammed Mestari, Mohammed Benzirar, Nadia Saber, and Meryem Khouil Abstract— This paper develops a neural network architecture and a new processing method for solving in real time, the nonlinear equality constrained multiobjective optimization problem (NECMOP), where several nonlinear objective functions must be optimized in a conflicting situation. In this processing method, the NECMOP is converted to an equivalent scalar optimization problem (SOP). The SOP is then decomposed into several-separable subproblems processable in parallel and in a reasonable time by multiplexing switched capacitor circuits. The approach which we propose makes use of a decomposition– coordination principle that allows nonlinearity to be treated at a local level and where coordination is achieved through the use of Lagrange multipliers. The modularity and the regularity of the neural networks architecture herein proposed make it suitable for very large scale integration implementation. An application to the resolution of a physical problem is given to show that the approach used here possesses some advantages of the point of algorithmic view, and provides processes of resolution often simpler than the usual techniques. Index Terms— Multiplexing switched capacitor circuits, neural networks architecture, nonlinear constrained multiobjective optimization problem and scalar optimization problem (SOP).
I. I NTRODUCTION
I
N ENGINEERING and management science optimization problems there often exist several criteria, which must be considered in a conflicting situation [20]. This situation is formulated as a nonlinear constrained multicriterion optimization (multiperformance, multiobjective, or vector optimization) problem [nonlinear equality constrained multiobjective optimization problem (NECMOP)] where not a single objective function but several functions are to be minimized or maximized simultaneously. We assume that the objective functions and the constraints are nonlinear functions of all or some of the variables under consideration. Direct application of classic optimization methods to NECMOPs is, in general, difficult to carry out. These methods Manuscript received March 29, 2014; revised July 16, 2014 and December 18, 2014; accepted December 27, 2014. Date of publication January 30, 2015; date of current version September 16, 2015. M. Mestari, N. Saber, and M. Khouil are with the Department of Mathematics and Computer Sciences, Laboratoire des Signaux, Systèmes Distribués et Intelligence Artificielle, École Normale Supérieure de l’Enseignement Technique, Mohammedia 20800, Morocco (e-mail: [email protected]; [email protected]; [email protected]). M. Benzirar is with the Department of Physics, Faculty of Sciences and Techniques, Mohammedia 20650, Morocco (e-mail: [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TNNLS.2015.2388511
often make use of very complex mathematical tools, and in the worst case, the amount of computation required may grow exponentially with the size of the problem. The search for new insights and effective solutions remains an active research endeavor. A number of new approaches using genetic algorithms have proposed in several manners to solve NECMOP [1]–[12]. Unfortunately, when the constraints of the problem considered become too difficult to satisfy or when the objective space is nonconvex, the multiobjective genetic algorithms converge with difficulty to optimal Pareto front [2] and [7]. In addition, these algorithms determine any bound optimal Pareto front of the problem. One of the most promising applications of artificial neural networks (ANNs) is probably in the area of different classes of optimization problems. The ability of analog neural networks to process simultaneously a large number of variables makes it possible to find solutions for complex multiobjective optimization problems in almost real time. In fact, with the advance of analog VLSI technologies and electrooptics, it is feasible today to design specialized analog hardware chips that can solve a specific optimization problem considerably faster than using a sequential algorithm on a general purpose digital computer and sometimes even faster than on a specialized digital hardware. The main objective of this paper is to present a new approach and new techniques of calculation in real time by multiplexing switched capacitor (SC) circuits to contribute to the resolution of the NECMOP. The procedure followed in this paper consists of modifying the original problem by addition of supplementary variables and constraints. At first sight, this appears to increase complexity. In fact, such complications will be largely compensated by the structural simplification it introduces. The approach which we propose makes use of a decomposition–coordination principle, that allows nonlinearity to be treated at a local level and where coordination is achieved through the use of Lagrange multipliers [31]. This localization of nonlinearities principle is, in effect, one of the directing principles of this paper. In this approach, the NECMOP is converted to an equivalent SOP with a single objective (cost) function. The SOP is then solved by mapping the differential equations into corresponding difference equations, which are simulated by discrete-time computing units (i.e., discretetime integrators). Such discrete-time computing units can be realized using SC techniques [31]. As will be shown later, employing SC techniques enables us to considerably reduce the number of basic elements.
2162-237X © 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
MESTARI et al.: SOLVING NECMOPs USING NEURAL NETWORKS
The construction and the calculation of the cost functions necessitate to use the weighting network (WN), that will be realized by incorporating adaptive nonlinear building blocks and a linear neuron, and the winner-take-all network (WTAN), that will be implemented on the basis of the order statistic filters adjustable order statistic networks (AOSNETs) [29]–[31]. Generally speaking, the function of the WTAN is to select the strongest (largest) objective function from a set of several objective functions considered in conflicting situations. On the other hand, the function of the WN is to transform a multiobjective problem in a problem into one objective while aggregating the different objectives under the form of a weighted sum. All neural networks herein proposed are hardware implementable with very simple circuit elements, such as SC circuits and linear operational amplifiers. The convergence of the analog algorithm is equivalent to the stability of the neural system, which is proven using the Lyapunov function approach. Unlike the digital parallel algorithms, the proposed neural networks solve NECMOP without using floating-point arithmetics, the central processing unit and the memory chips. An application to the resolution of a problem of the physics concerning the calculation of the coefficient of diffusion of a laser beam crossing a heated turbulent medium is given to show that the approach used here possesses some advantages of the point of algorithmic view, and provides processes of resolution often simpler than the usual techniques. This application consists in calculating the local diffusion coefficient Dμ (x) by solving a multiobjective optimal control problem, where the coefficient Dμ (x) designates the control action of a nonlinear dynamic system. This system is described by a mathematical model based on Markov random processes and on the Einstein–Fokker–Planck–Kolmogorov equations as a starting point. The organization of this paper is as follows. Section II is devoted to the formulation of the problem and to its transformation to an equivalent SOP with a single objective (cost) function. Section III is dedicated to the analysis and the resolution of the SOP problem. At the end of this section, an important study of the stability of the proposed method is given. In Section IV, we demonstrate how to develop in real time, an appropriate neural network architecture for implementing the method proposed in Section III. Section V is devoted, first, to the formulation of the physical problem and to its transformation in terms of nonlinear dynamic systems, in the form of an NECMOP problem. Second, we show how the resolution of such NECMOP problem makes calculation of the diffusion coefficient Dμ along the laser beam trajectory is possible. The conclusion is provided in Section VI. II. S TATEMENT OF THE P ROBLEM Consider the nonlinear discrete-time dynamical systems described by x k+1 = f (x k , u k ) (1) x 0 given
2501
where x k ∈ Rn and u k ∈ Rm . x k and u k represent, respectively, the state and the input of the system at time k. The problem we propose to solve consists of determining the control inputs, which simultaneously optimize (minimize or maximize) several objective functions considered in conflicting situations. For example, to determine the control inputs that enable a desired state x d to be achieved at time N, so doing with minimum efforts; i.e., to determine the sequence of control inputs u ∗k (k = 0, 1, . . . , N − 1), which simultaneously minimizes J1 (x, u) =
N−1 1 ||u k ||2 2
(2)
k=0
and 1 ||x N − x d ||2 (3) 2 in the presence of equality constraints (1), where x = [x 0T , x 1T , . . . , x NT ]T and u = [u 0T , u 1T , . . . , u TN−1 ]T . This problem can be stated as a multiobjective optimization problem min {J1 (x, u), J2 (x, u), . . . , J p (x, u)} {u ∗k | 0≤k≤N−1} (4) s.t. x k+1 = f (x k , u k ) with x 0 = x(0) J2 (x, u) =
where x
= [x 0T , x 1T , . . . , x NT ]T ;
u x(0)
= [u 0T , u 1T , . . . , u TN−1 ]T ; given initial condition.
For the sake of convenience, we will assume that all the objective functions (J1 , J2 , . . . , J p ) are to be minimized. The function to be maximized can be converted into the form, which allows their minimization by max Ji (x, u) = − min(−Ji (x, u)). u
u
(5)
In general, the objective functions Ji (x, u) are conflicting, so the final design must be a compromise among them. Unlike standard mathematical programming with a single objective function, an optimal solution in the sense that one minimizes all objective functions simultaneously does not necessarily exist in multiobjective optimization problems. Since, in general, there exists no complete optimal solution y ∗ , which satisfies J1 (y ∗ ) ≤ J1 (y),
J2 (y ∗ ) ≤ J2 (y), . . . , J p (y ∗ ) ≤ J p (y)
for all y ∈ Rn+m , we usually determine a Pareto optimal solution [13], [20], [24]. The Pareto optimum gives a set of noninferior solutions, i.e., solutions for which there is no way of improving any objective function without worsening at least one other objective function. For the multiobjective optimization problem (4), we will use the concept of Pareto or efficiency points [13], [20], [24], [28]. A point y ∗ is called a local efficient point of the set of objective functions {J1 , J2 , . . . , J p }, if there exists a neighborhood of y ∗ such that for any other point y in this neighborhood at least one of the functions increases in value relatively at the point y ∗ , i.e., Ji (y) > Ji (y ∗ ) for some i (i = 1, 2, . . . , p).
2502
IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. 26, NO. 10, OCTOBER 2015
The solution of nonlinear constrained multiobjective problems [as (4)] generally constitutes a difficult and often a frustrating task. Solving this problem requires determining a set of feasible points of the decision space, representing a compromise according to the different objectives of the problem. Both solving techniques we adopt thereafter to process the problem (4) are non-Pareto approaches. They transform a multiobjective optimization problem into a single objective problem. Whether in the form of a Chebyshev distance (minmax method) or as a weighted sum (weighting method). This transformation has the advantage of being easily implementable as we will see in Section IV using the decomposition–coordination method that will be introduced in Section III-A, whose principle is described by the diagram in Fig. 3. The principles of these two non-Pareto methods previously mentioned are given by the following methods. 1) Minimax Method [13], [20], [24], [37]: The desired solution is the one which gives the smallest value of the maximum values of all the objective functions Ji , i.e., the solution that minimizes the following scalar function: E(x, u) = max {ωi Ji (x, u)} 1≤i≤ p
where ωi ≥ 0 are the weighting coefficients representing the relative importance of the function Ji . We usually assume that the weighting coefficients satisfy the relation p
ωi = 1.
i=1
2) Weighting Method [13], [20], [24], [37]: This method utilizes as a scalar function the weighting function defined as E(x, u) =
p i=1
ωi Ji (x, u)
p where ωi ≥ 0 and i=1 ωi = 1. It is clear that solving a multiobjective optimization problem using these two methods for a fixed weight vector is used to calculate only a few Pareto optimal solutions. To obtain a set containing a large number of Pareto optimal solutions, we must change the values of ω every time. The technique of weighted sum cannot handle completely nonconvex problems and is therefore very sensitive to the form of the Pareto front. Fig. 1 shows this case in dimension 2. In fact, for a fixed weight vector ω = (ω1 , ω2 ), the attainable optimal value for the objective function created is d = ω1 J1 + ω2 J2 . Both Pareto optimal points found are A and C. By varying the vector ω, it is possible to find other optimal Pareto points. Only all these points will lie on the convex portions of the surface of compromise. There is no possible value for ω to find such point B. Indeed, this technique does not allow approaching the entire Pareto front when it is not
Fig. 1.
2-D example of Pareto optimality and the Pareto front.
convex. However, the Minimax method can treat nonconvex problems provided that the reference point is chosen wisely. Therefore, problem (4) can be formulated as ⎧ ⎨ ∗ min E(x, u) {u k | 0≤k≤N−1} (6) ⎩s.t. x k+1 = f (x k , u k ) with x 0 = x(0) where x
= [x 0T , x 1T , . . . , x NT ]T ;
u E(x, u)
= [u 0T , u 1T , . . . , u TN−1 ]T ; p = max1≤i≤ p {ωi Ji (x, u)} or i=1 ωi Ji (x, u) p with ωi ≥ 0 and i=1 ωi = 1;
x(0)
given initial condition.
The difficulty in numerically solving the problem (6) is that in the worst case, the amount of computation required may grow exponentially with the size of the problem. In this paper, we develop and implement new techniques and neural networks, which enable us to solve this problem in real time. The basic procedure employed in this paper consists of decomposing the nonlinear discrete-time dynamical system described by (1) into several subsystems, allowing the decomposition of the overall problem (6) into a number of smaller separable subsystems processable in parallel and in reasonable time by the ANNs. III. A NALYSIS OF THE P ROBLEM We now concentrate on the issue of how to solve the problem (6). The key step is to decompose the system described by (1) into several subsystems. The decomposition method employed allows passage from a dynamic nonlinear system to a group of N interconnected subsystems arranged according to a simple serial structure (Fig. 2), where z k is the output of subsystem k z k = f (x k , u k ), k = 0, 1, 2, . . . , N − 2
(7)
x k = z k−1 , k = 1, 2, . . . , N − 1.
(8)
and
MESTARI et al.: SOLVING NECMOPs USING NEURAL NETWORKS
Fig. 2.
2503
Overall system made up of N interconnected subsystems.
In terms of interconnected subsystems, problem (6) can be stated as follows:
satisfies the following: 1 ∂E ∂f T ∗ + μ + βk∗ = 0, for 1 ≤ k ≤ N − 1 N ∂ xk ∂ xk k (14) T 1 ∂E ∂f ∇u k L = + μ∗ = 0, for 0 ≤ k ≤ N − 1 (15) N ∂u k ∂u k k ∇xk L =
⎧ min E(x, u) ⎪ ⎨{u ∗ | 0≤k≤N−1} k
(9) s.t. z k = f (x k , u k ), k = 0, 1, . . . , N − 2 ⎪ ⎩ and x k = z k−1 , k = 1, 2, . . . , N − 1 with x 0 = x(0)
∗ = 0, for 0 ≤ k ≤ N − 2 ∇zk L = −μ∗k − βk+1
where
∇μk L =
f (x k∗ , u ∗k ) − z k∗ = 0, ∗ x k∗ − z k−1 = 0, for
for 0 ≤ k ≤ N − 1
(16) (17)
x
= [x 0T , x 1T , . . . , x NT ]T ;
u E(x, u)
= [u 0T , u 1T , . . . , u TN−1 ]T ; p = max1≤i≤ p {ωi Ji (x, u)} or i=1 ωi Ji (x, u) p with ωi ≥ 0 and i=1 ωi = 1;
Solving the equality-constrained minimization problem (9) is equivalent to solving the associated system of differential equations (14)–(18).
x(0)
given initial condition.
A. Decomposition–Coordination Method
Note that even if all functions Ji (x, u) are differentiable the cost function max1≤i≤ p {ωi Ji (x, u)} will have corners where differentiability fails; i.e., the cost function max1≤i≤ p {ωi Ji (x, u)} has discontinuous partial first-order derivatives at points where two or more functions Ji (x, u) are equal to max1≤i≤ p {ωi Ji (x, u)} even if the Ji (x, u) have continuous first-order partial derivatives. For this problem, we can construct the ordinary Lagrange function L=
N−1
Lk
(10)
k=0
where 1 E(x, u) + μ0T ( f (x 0 , u 0 ) − z 0 ) (11) N 1 L k = E(x, u) + μkT ( f (x k , u k ) − z k ) + βkT (x k − z k−1 ) N for 1 ≤ k ≤ N − 2 (12) 1 L N−1 = E(x, u) + μTN−1 ( f (x N−1 , u N−1 ) − x d ) N T (x N−1 − z N−2 ) (13) + β N−1 L0 =
μk (n components) and βk (n components) are the Lagrange multiplier vectors introduced to consider the equality constraints (7) and (8). The derivations of the ordinary Lagrange function (10) enable us to transform the equality-constrained minimization problem (9) into a set of differential equations on the basis of which we design ANN architectures. A stationary (equilibrium) point (x k∗ , u ∗k , μ∗k , βk∗ , z k∗ ), according to the Karush-Khun-Tucker conditions [26],
∇βk L =
1 ≤ k ≤ N − 1.
(18)
One way to solve the equality constrained minimization problem (9) is to decompose the treatment of the associated system of differential equations (14)–(18) between two levels, not arbitrarily but in such a way as to obtain a separable form at the inferior level, i.e., each subproblem k only brings about the intervention of unknown variables of indices k, for 0 ≤ k ≤ N − 1. In effect, the processing of the system of (14)–(18) is divided between two levels. The upper level handles (16) and (18) and fixes z k (0 ≤ k ≤ N − 2) and βk (1 ≤ k ≤ N − 1), which it proposes to the lower level. This allows the decomposition of problem (9) to the lower level; each subproblem becomes as follows. 1) Subproblem 0 ⎧ z being given by the upper level ⎪ ⎨ 0 min E(x, u) (19) ⎪ ⎩s.t. f (x , u ) = z . 0 0 0 2) Subproblem k ⎧ z (k = 1, 2, . . . , N − 2) and βk (k = 1, 2, . . . , N − 2) ⎪ ⎪ k ⎪ ⎨being given by the upper level (20) ⎪min E(x, u) + βkT x k − c, where c = βkT z k−1 ⎪ ⎪ ⎩ s.t. z k = f (x k , u k ). 3) Subproblem N − 1 ⎧ z N−2 and β N−1 being given by the upper level ⎪ ⎪ ⎪ ⎨min E(x, u) + β T x N−1 N−1 − c (21) T ⎪ where c = β N−1 z N−2 ⎪ ⎪ ⎩ s.t. f (x N−1 , u N−1 ) − x d = 0. The solution of each subproblem thus corresponds to the processing of (14), (15), and (17) for z k (k = 0, 1, . . . , N − 2)
2504
IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. 26, NO. 10, OCTOBER 2015
and βk (k = 1, . . . , N − 1) supplied by the upper level. Thus, applying a gradient method, we obtain the system of differential equations
system of difference equations (28)–(30) and determination of ( j) ( j) ( j) ( j) ( j) ( j) the variables x k∗ (z k , βk ), u ∗k (z k , βk ), and μ∗k (z k , βk ), which, respectively, satisfy (14), (15), and (17). The results
d xk = −λx ∇xk L (22) dt du k = −λu ∇u k L (23) dt dμk = λμ ∇μk L (24) dt where λx > 0, λu > 0, and λμ > 0. Equations (22)–(24) can be written in the following scalar form:
n d x kq 1 ∂E ∂ fi = −λx + μki + βkq dt N ∂ x kq ∂ x kq
x k∗ (z k , βk ) and μ∗k (z k , βk ) are supplied for the upper level, which verifies whether the previously supplied information was correct and corrects it if necessary, so that the lower level can recommence its work with information of increased validity. Such correction is necessary to bring about evolution toward the satisfaction of (16) and (18). The upper level proceeds progressively in making the necessary correction to the coordination parameters, coming closer and closer to the satisfaction of coordination equations (16) and (18). Thus, the coordination parameters at iteration j + 1 are improvements on the coordination parameters at iteration j ( j ) ( j ) ( j +1) ( j) ( j) = z kq − λz −μ∗iq z k , βk − βk+1,q z kq
i=1
du kq dt
( j)
k = 1, . . . , N − 1 and q = 1, . . . , n
(26)
k = 1, 2, . . . , N − 1 and q = 1, 2, . . . , n (32)
(27)
where λz > 0 and λβ > 0. The solution of the system of difference equations (28)–(30) is thus repeated until satisfactory coordination is obtained, i.e., satisfaction of coordination equations (16) and (18). The method thus defined is shown in Fig. 3.
( j +1)
βkq
i=1
i=1
k = 1, . . . , N − 1 and q = 1, . . . , n
n (l) 1 ∂ E (l) ∂ f i (l) (l) = u kq − λu + μ N ∂u kq ∂u kq ki
(28)
k = 1, . . . , N − 1 and q = 1, . . . , m (l) (l) ( j ) = μ(l) − z kq , kq + λμ f q x k , u k
(29)
k = 0, . . . , N − 1 and q = 1, . . . , n.
(30)
i=1
μ(l+1) kq
( j)
(25)
where f (x k , u k ) = ( f 1 (x k , u k ), f 2 (x k , u k ), . . . , fn (x k , u k ))T , x k = (x k1 , . . . , x kn )T , u k = (u k1 , . . . , u kn )T , μk = (μk1 , . . . , μkn )T , βk = (βk1 , . . . , βkn )T , and z k = (z k1 , . . . , z kn )T . To realize the discrete-time network, we can transform the system of differential equations (25)–(27) into an appropriate system of difference equations. For example, using the forward Euler rule, the system of differential equations (25)–(27) can be converted into the system of difference equations
n (l) 1 ∂ E (l) ∂ f i ( j) (l+1) (l) + μki + βkq x kq = x kq − λx N ∂ x kq ∂ x kq
(l+1)
( j)
k = 1, . . . , N − 1 and q = 1, . . . , n
n 1 ∂E ∂ fi = −λu + μki N ∂u kq ∂u kq
dμkq = λμ ( f q (x k , u k ) − z kq ) dt k = 0, . . . , N − 1 and q = 1, . . . , n
u kq
( j)
Equations (28)–(30) can only, then, be solved locally if it ( j) possesses the necessary information, βk (k = 1, . . . , N − 1) ( j) and z k (k = 0, . . . , N − 2) sent from the upper level. To transmit the information necessary to the functioning of the lower level with a view to overall optimization [i.e., satisfaction of all equations (14)–(18)], it is essential to coordinate the two levels. To ensure such coordination, the ( j) upper level works simultaneously on βk (k = 1, . . . , N − 1) ( j) and z k (k = 0, . . . , N − 2), which make up the coordination parameters. These coordination parameters are considered as known within the lower level, allowing local resolution of the
k = 0, 1, . . . , N − 2 and q = 1, 2, . . . , n (31) ( j) ( j) ( j ) ( j) ∗ z k , βk − z k−1,q = βkq + λβ x kq
B. Analysis Stability The aim of this section is to demonstrate that the convergence of the method described by the diagram in Fig. 3 may be reduced to that of the coordinating level. To facilitate this paper, we will introduce the following notations. x k and call v k∗ (z k∗ , βk∗ ), μ∗k (z k∗ , βk∗ ), Let us posit v k = uk ( j) ( j) z k∗ , and βk∗ the solution sought. Let us also call v k∗ (z k , βk ) ( j)
( j)
and μ∗k (z k , βk ), the variables to be determined at the lower level by the system of equations (28)–(30) so that (14), (15), ( j) ( j) and (17) may be locally satisfied, and z k and βk , the coordination variables to be modified at the upper level following algorithms of coordination (31) and (32). Let us put the problem, which interests us in a more condensed form ⎧ Gk (v k , μk , βk ) = 0 ⎪ ⎪ ⎪ lower level equations ⎪ ⎨ Pk (v k , z k ) = 0 (33) Rk (v k , z k−1 ) = 0 ⎪ ⎪ ⎪ ⎪ ⎩ Hk (μk , βk+1 ) = 0 upper level equations where Gk =
L xk , Pk = μk L, Rk = βk L, and u k L
Hk = zk L. For the solution sought, we may write ⎧ ∗ G =0 ⎪ ⎪ ⎨ ∗k Pk = 0 R∗ = 0 ⎪ ⎪ ⎩ k∗ Hk = 0
(34)
MESTARI et al.: SOLVING NECMOPs USING NEURAL NETWORKS
Fig. 3.
2505
Information transfer between lower and upper levels.
where Gk∗ = Gk (v k∗ (βk∗ , z k∗ ), μ∗k (βk∗ , z k∗ ), βk∗ ), Pk∗ = ∗ ), and H∗ = Pk (v k∗ (βk∗ , z k∗ ), z k∗ ), Rk∗ = Rk (v k∗ (βk∗ , z k∗ ), z k−1 k ∗ ∗ ∗ ∗ Hk (μk (βk , z k ), βk+1 ). At each iteration j of the coordination loop, we have ( j) Gk = 0 (35) ( j) Pk = 0 ( j)
( j)
( j)
( j)
( j)
( j)
= Gk (v k∗ (βk , z k ), μ∗k (βk , z k ), βk ), and where Gk ( j) ( j) ( j) ( j) Pk = Pk (v k∗ (βk , z k ), z k ). ( j) ( j) ( j) = v k∗ (βk , z k ) − v k∗ (βk∗ , z k∗ ), Let us posit ev k ( j) ( j ) ( j ) ( j) ( j) eμk = μ∗k (βk , z k ) − μ∗k (βk∗ , z k∗ ), ezk = z k − z k∗ , and ( j)
( j)
eβk = βk − βk∗ , which designate the errors calculated at iteration j of the coordination loop.
( j)
( j)
( j)
( j)
Theorem 1: Let ev k , eμk , ezk , and eβk be the errors computed at the iteration j of the coordination loop. Then, ( j) ( j) ( j) ( j) ev k −→ 0 and eμk −→ 0 if ezk −→ 0 and eβk −→ 0. Proof: Linearization of the lower level equations in the neighborhood of solution allows us to write ( j)
Gk∗ +
∂Gk∗ ( j ) ∂Gk∗ ( j ) ∂Gk∗ ( j ) e + e + e ∂v k v k ∂μk μk ∂βk βk
(36)
( j)
Pk∗ +
∂Pk∗ ( j ) ∂Pk∗ ( j ) e + ez ∂v k v k ∂z k k
(37)
Gk Pk
where considering (34) and (35), (36) and (37) become ∂Gk∗ ( j ) ∂Gk∗ ( j ) ∂Gk∗ ( j ) e + e + e =0 ∂v k v k ∂μk μk ∂βk βk
(38)
2506
IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. 26, NO. 10, OCTOBER 2015
and
Using (45) and (46) in (42), we have ∂Pk∗ ( j ) ∂Pk∗ ( j ) e + ez = 0. ∂v k v k ∂z k k
(39)
( j ) =
k=0
( j)
Equations (38) and (39) demonstrate that ev k −→ 0 and ( j) eμk
−→ 0 if
( j) ez k
−→ 0 and
( j) eβ k
−→ 0.
N−1
∗ ∂Rk∗ ( j ) ( j ) T ∂Rk ( j ) ev k + eβ k ez ∂v k ∂z k−1 k−1 1 ( j )T ( j ) ( j )T ( j ) + λ2 Hk Hk + Rk Rk 2 ( j )T
+ eβ k
♦
( j)
( j)
One need only study, then, the variation of ezk and eβk during the iterations of the coordination loop to determine the convergence of the method illustrated by the diagram in Fig. 3. Given the very major role played by coordination (the upper level’s task), a more in-depth study of the convergence (stability) of coordination algorithms (31) and (32) becomes necessary. Let us choose as Lyapunov function
∗ ∗ ( j ) T ∂Hk ( j ) ( j ) T ∂Hk ( j ) λ −ezk eμk − ezk e ∂μk ∂βk+1 βk+1
= A( j )λ2 + B( j )λ
(47)
where A( j ) =
N−1
( j )T
Hk
( j)
( j )T
H k + Rk
( j)
Rk ≥ 0
k=0
and ( j ) =
N−1
1 2
( j )T ( j ) ez k
ez k
( j )T ( j )
+ eβk eβk ≥ 0.
(40)
k=0
( j ) = ( j + 1) − ( j ) N−1 1 ( j +1)T ( j +1) ( j )T ( j ) ez k − ez k ez k = ez k 2 i=0
( j +1) T ( j +1) eβ k
+ eβ k
( j )T ( j )
− eβ k eβ k
.
(41)
The development of (41) allows us to write ( j ) =
( j )T
ez k
k=0
+
( j)
( j )T
( j)
ezk + eβk eβk
1 ( j )T ( j ) ( j )T ( j) ezk ezk + eβk eβk (42) 2
where ( j)
( j +1)
− ezk = −λHk
( j)
( j +1)
− eβk = λFk
ezk = ezk
eβk = eβk
( j) ( j)
( j)
( j)
( j)
Hk +
( j)
Rk∗ +
Rk
−ezk
∗ ∂Hk∗ ( j ) ( j ) T ∂Hk ( j ) eμk − ezk e ∂μk ∂βk+1 βk+1
( j )T + eβ k
∗ ∂Rk∗ ( j ) ( j ) T ∂Rk ( j ) . e + eβ k ez ∂v k v k ∂z k−1 k−1
( j )T
(∗)
∂Hk∗ ( j ) ∂Hk∗ ( j ) eμk + e ∂μk ∂βk+1 βk+1
∂Rk∗ ( j ) ∂Rk∗ ( j ) ev k + ez . ∂v k ∂z k−1 k−1
Considering (34), (43), and (44), we may write ∗ ∂Hk ( j ) ∂Hk∗ ( j ) ( j) ezk = −λ eμk + eβk+1 ∂μk ∂βk+1 ∗ ∂Rk ( j ) ∂Rk∗ ( j ) ( j) eβk = λ ev k + ezk−1 . ∂v k ∂z k−1
Theorem 2: Let λ be the adaptive coefficient for the coordinative algorithms. Then, the convergence is guaranteed if one of the matrices (∂Gk∗ /∂v k )T (k = 0, 1, . . . , N − 1) is positive definite and the others are only positive semidefinite, and furthermore, if A( j ) = 0, λ must be chosen as: 0 < λ 0, the equation ( j ) = 0 admits two real distinct roots 0 and −(B( j )/A( j )). Then, ( j ) < 0 if B( j ) < 0 (i.e., exactly one of the matrices
MESTARI et al.: SOLVING NECMOPs USING NEURAL NETWORKS
2519
(∂Gk∗ /∂v k )T (k = 0, 1, . . . , N − 1) is definite positive and the others are only positive semidefinite (78) and λ is chosen as B( j ) B( j ) . (79) 0
IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. 26, NO. 10, OCTOBER 2015
Solving Nonlinear Equality Constrained Multiobjective Optimization Problems Using Neural Networks Mohammed Mestari, Mohammed Benzirar, Nadia Saber, and Meryem Khouil Abstract— This paper develops a neural network architecture and a new processing method for solving in real time, the nonlinear equality constrained multiobjective optimization problem (NECMOP), where several nonlinear objective functions must be optimized in a conflicting situation. In this processing method, the NECMOP is converted to an equivalent scalar optimization problem (SOP). The SOP is then decomposed into several-separable subproblems processable in parallel and in a reasonable time by multiplexing switched capacitor circuits. The approach which we propose makes use of a decomposition– coordination principle that allows nonlinearity to be treated at a local level and where coordination is achieved through the use of Lagrange multipliers. The modularity and the regularity of the neural networks architecture herein proposed make it suitable for very large scale integration implementation. An application to the resolution of a physical problem is given to show that the approach used here possesses some advantages of the point of algorithmic view, and provides processes of resolution often simpler than the usual techniques. Index Terms— Multiplexing switched capacitor circuits, neural networks architecture, nonlinear constrained multiobjective optimization problem and scalar optimization problem (SOP).
I. I NTRODUCTION
I
N ENGINEERING and management science optimization problems there often exist several criteria, which must be considered in a conflicting situation [20]. This situation is formulated as a nonlinear constrained multicriterion optimization (multiperformance, multiobjective, or vector optimization) problem [nonlinear equality constrained multiobjective optimization problem (NECMOP)] where not a single objective function but several functions are to be minimized or maximized simultaneously. We assume that the objective functions and the constraints are nonlinear functions of all or some of the variables under consideration. Direct application of classic optimization methods to NECMOPs is, in general, difficult to carry out. These methods Manuscript received March 29, 2014; revised July 16, 2014 and December 18, 2014; accepted December 27, 2014. Date of publication January 30, 2015; date of current version September 16, 2015. M. Mestari, N. Saber, and M. Khouil are with the Department of Mathematics and Computer Sciences, Laboratoire des Signaux, Systèmes Distribués et Intelligence Artificielle, École Normale Supérieure de l’Enseignement Technique, Mohammedia 20800, Morocco (e-mail: [email protected]; [email protected]; [email protected]). M. Benzirar is with the Department of Physics, Faculty of Sciences and Techniques, Mohammedia 20650, Morocco (e-mail: [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TNNLS.2015.2388511
often make use of very complex mathematical tools, and in the worst case, the amount of computation required may grow exponentially with the size of the problem. The search for new insights and effective solutions remains an active research endeavor. A number of new approaches using genetic algorithms have proposed in several manners to solve NECMOP [1]–[12]. Unfortunately, when the constraints of the problem considered become too difficult to satisfy or when the objective space is nonconvex, the multiobjective genetic algorithms converge with difficulty to optimal Pareto front [2] and [7]. In addition, these algorithms determine any bound optimal Pareto front of the problem. One of the most promising applications of artificial neural networks (ANNs) is probably in the area of different classes of optimization problems. The ability of analog neural networks to process simultaneously a large number of variables makes it possible to find solutions for complex multiobjective optimization problems in almost real time. In fact, with the advance of analog VLSI technologies and electrooptics, it is feasible today to design specialized analog hardware chips that can solve a specific optimization problem considerably faster than using a sequential algorithm on a general purpose digital computer and sometimes even faster than on a specialized digital hardware. The main objective of this paper is to present a new approach and new techniques of calculation in real time by multiplexing switched capacitor (SC) circuits to contribute to the resolution of the NECMOP. The procedure followed in this paper consists of modifying the original problem by addition of supplementary variables and constraints. At first sight, this appears to increase complexity. In fact, such complications will be largely compensated by the structural simplification it introduces. The approach which we propose makes use of a decomposition–coordination principle, that allows nonlinearity to be treated at a local level and where coordination is achieved through the use of Lagrange multipliers [31]. This localization of nonlinearities principle is, in effect, one of the directing principles of this paper. In this approach, the NECMOP is converted to an equivalent SOP with a single objective (cost) function. The SOP is then solved by mapping the differential equations into corresponding difference equations, which are simulated by discrete-time computing units (i.e., discretetime integrators). Such discrete-time computing units can be realized using SC techniques [31]. As will be shown later, employing SC techniques enables us to considerably reduce the number of basic elements.
2162-237X © 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
MESTARI et al.: SOLVING NECMOPs USING NEURAL NETWORKS
The construction and the calculation of the cost functions necessitate to use the weighting network (WN), that will be realized by incorporating adaptive nonlinear building blocks and a linear neuron, and the winner-take-all network (WTAN), that will be implemented on the basis of the order statistic filters adjustable order statistic networks (AOSNETs) [29]–[31]. Generally speaking, the function of the WTAN is to select the strongest (largest) objective function from a set of several objective functions considered in conflicting situations. On the other hand, the function of the WN is to transform a multiobjective problem in a problem into one objective while aggregating the different objectives under the form of a weighted sum. All neural networks herein proposed are hardware implementable with very simple circuit elements, such as SC circuits and linear operational amplifiers. The convergence of the analog algorithm is equivalent to the stability of the neural system, which is proven using the Lyapunov function approach. Unlike the digital parallel algorithms, the proposed neural networks solve NECMOP without using floating-point arithmetics, the central processing unit and the memory chips. An application to the resolution of a problem of the physics concerning the calculation of the coefficient of diffusion of a laser beam crossing a heated turbulent medium is given to show that the approach used here possesses some advantages of the point of algorithmic view, and provides processes of resolution often simpler than the usual techniques. This application consists in calculating the local diffusion coefficient Dμ (x) by solving a multiobjective optimal control problem, where the coefficient Dμ (x) designates the control action of a nonlinear dynamic system. This system is described by a mathematical model based on Markov random processes and on the Einstein–Fokker–Planck–Kolmogorov equations as a starting point. The organization of this paper is as follows. Section II is devoted to the formulation of the problem and to its transformation to an equivalent SOP with a single objective (cost) function. Section III is dedicated to the analysis and the resolution of the SOP problem. At the end of this section, an important study of the stability of the proposed method is given. In Section IV, we demonstrate how to develop in real time, an appropriate neural network architecture for implementing the method proposed in Section III. Section V is devoted, first, to the formulation of the physical problem and to its transformation in terms of nonlinear dynamic systems, in the form of an NECMOP problem. Second, we show how the resolution of such NECMOP problem makes calculation of the diffusion coefficient Dμ along the laser beam trajectory is possible. The conclusion is provided in Section VI. II. S TATEMENT OF THE P ROBLEM Consider the nonlinear discrete-time dynamical systems described by x k+1 = f (x k , u k ) (1) x 0 given
2501
where x k ∈ Rn and u k ∈ Rm . x k and u k represent, respectively, the state and the input of the system at time k. The problem we propose to solve consists of determining the control inputs, which simultaneously optimize (minimize or maximize) several objective functions considered in conflicting situations. For example, to determine the control inputs that enable a desired state x d to be achieved at time N, so doing with minimum efforts; i.e., to determine the sequence of control inputs u ∗k (k = 0, 1, . . . , N − 1), which simultaneously minimizes J1 (x, u) =
N−1 1 ||u k ||2 2
(2)
k=0
and 1 ||x N − x d ||2 (3) 2 in the presence of equality constraints (1), where x = [x 0T , x 1T , . . . , x NT ]T and u = [u 0T , u 1T , . . . , u TN−1 ]T . This problem can be stated as a multiobjective optimization problem min {J1 (x, u), J2 (x, u), . . . , J p (x, u)} {u ∗k | 0≤k≤N−1} (4) s.t. x k+1 = f (x k , u k ) with x 0 = x(0) J2 (x, u) =
where x
= [x 0T , x 1T , . . . , x NT ]T ;
u x(0)
= [u 0T , u 1T , . . . , u TN−1 ]T ; given initial condition.
For the sake of convenience, we will assume that all the objective functions (J1 , J2 , . . . , J p ) are to be minimized. The function to be maximized can be converted into the form, which allows their minimization by max Ji (x, u) = − min(−Ji (x, u)). u
u
(5)
In general, the objective functions Ji (x, u) are conflicting, so the final design must be a compromise among them. Unlike standard mathematical programming with a single objective function, an optimal solution in the sense that one minimizes all objective functions simultaneously does not necessarily exist in multiobjective optimization problems. Since, in general, there exists no complete optimal solution y ∗ , which satisfies J1 (y ∗ ) ≤ J1 (y),
J2 (y ∗ ) ≤ J2 (y), . . . , J p (y ∗ ) ≤ J p (y)
for all y ∈ Rn+m , we usually determine a Pareto optimal solution [13], [20], [24]. The Pareto optimum gives a set of noninferior solutions, i.e., solutions for which there is no way of improving any objective function without worsening at least one other objective function. For the multiobjective optimization problem (4), we will use the concept of Pareto or efficiency points [13], [20], [24], [28]. A point y ∗ is called a local efficient point of the set of objective functions {J1 , J2 , . . . , J p }, if there exists a neighborhood of y ∗ such that for any other point y in this neighborhood at least one of the functions increases in value relatively at the point y ∗ , i.e., Ji (y) > Ji (y ∗ ) for some i (i = 1, 2, . . . , p).
2502
IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. 26, NO. 10, OCTOBER 2015
The solution of nonlinear constrained multiobjective problems [as (4)] generally constitutes a difficult and often a frustrating task. Solving this problem requires determining a set of feasible points of the decision space, representing a compromise according to the different objectives of the problem. Both solving techniques we adopt thereafter to process the problem (4) are non-Pareto approaches. They transform a multiobjective optimization problem into a single objective problem. Whether in the form of a Chebyshev distance (minmax method) or as a weighted sum (weighting method). This transformation has the advantage of being easily implementable as we will see in Section IV using the decomposition–coordination method that will be introduced in Section III-A, whose principle is described by the diagram in Fig. 3. The principles of these two non-Pareto methods previously mentioned are given by the following methods. 1) Minimax Method [13], [20], [24], [37]: The desired solution is the one which gives the smallest value of the maximum values of all the objective functions Ji , i.e., the solution that minimizes the following scalar function: E(x, u) = max {ωi Ji (x, u)} 1≤i≤ p
where ωi ≥ 0 are the weighting coefficients representing the relative importance of the function Ji . We usually assume that the weighting coefficients satisfy the relation p
ωi = 1.
i=1
2) Weighting Method [13], [20], [24], [37]: This method utilizes as a scalar function the weighting function defined as E(x, u) =
p i=1
ωi Ji (x, u)
p where ωi ≥ 0 and i=1 ωi = 1. It is clear that solving a multiobjective optimization problem using these two methods for a fixed weight vector is used to calculate only a few Pareto optimal solutions. To obtain a set containing a large number of Pareto optimal solutions, we must change the values of ω every time. The technique of weighted sum cannot handle completely nonconvex problems and is therefore very sensitive to the form of the Pareto front. Fig. 1 shows this case in dimension 2. In fact, for a fixed weight vector ω = (ω1 , ω2 ), the attainable optimal value for the objective function created is d = ω1 J1 + ω2 J2 . Both Pareto optimal points found are A and C. By varying the vector ω, it is possible to find other optimal Pareto points. Only all these points will lie on the convex portions of the surface of compromise. There is no possible value for ω to find such point B. Indeed, this technique does not allow approaching the entire Pareto front when it is not
Fig. 1.
2-D example of Pareto optimality and the Pareto front.
convex. However, the Minimax method can treat nonconvex problems provided that the reference point is chosen wisely. Therefore, problem (4) can be formulated as ⎧ ⎨ ∗ min E(x, u) {u k | 0≤k≤N−1} (6) ⎩s.t. x k+1 = f (x k , u k ) with x 0 = x(0) where x
= [x 0T , x 1T , . . . , x NT ]T ;
u E(x, u)
= [u 0T , u 1T , . . . , u TN−1 ]T ; p = max1≤i≤ p {ωi Ji (x, u)} or i=1 ωi Ji (x, u) p with ωi ≥ 0 and i=1 ωi = 1;
x(0)
given initial condition.
The difficulty in numerically solving the problem (6) is that in the worst case, the amount of computation required may grow exponentially with the size of the problem. In this paper, we develop and implement new techniques and neural networks, which enable us to solve this problem in real time. The basic procedure employed in this paper consists of decomposing the nonlinear discrete-time dynamical system described by (1) into several subsystems, allowing the decomposition of the overall problem (6) into a number of smaller separable subsystems processable in parallel and in reasonable time by the ANNs. III. A NALYSIS OF THE P ROBLEM We now concentrate on the issue of how to solve the problem (6). The key step is to decompose the system described by (1) into several subsystems. The decomposition method employed allows passage from a dynamic nonlinear system to a group of N interconnected subsystems arranged according to a simple serial structure (Fig. 2), where z k is the output of subsystem k z k = f (x k , u k ), k = 0, 1, 2, . . . , N − 2
(7)
x k = z k−1 , k = 1, 2, . . . , N − 1.
(8)
and
MESTARI et al.: SOLVING NECMOPs USING NEURAL NETWORKS
Fig. 2.
2503
Overall system made up of N interconnected subsystems.
In terms of interconnected subsystems, problem (6) can be stated as follows:
satisfies the following: 1 ∂E ∂f T ∗ + μ + βk∗ = 0, for 1 ≤ k ≤ N − 1 N ∂ xk ∂ xk k (14) T 1 ∂E ∂f ∇u k L = + μ∗ = 0, for 0 ≤ k ≤ N − 1 (15) N ∂u k ∂u k k ∇xk L =
⎧ min E(x, u) ⎪ ⎨{u ∗ | 0≤k≤N−1} k
(9) s.t. z k = f (x k , u k ), k = 0, 1, . . . , N − 2 ⎪ ⎩ and x k = z k−1 , k = 1, 2, . . . , N − 1 with x 0 = x(0)
∗ = 0, for 0 ≤ k ≤ N − 2 ∇zk L = −μ∗k − βk+1
where
∇μk L =
f (x k∗ , u ∗k ) − z k∗ = 0, ∗ x k∗ − z k−1 = 0, for
for 0 ≤ k ≤ N − 1
(16) (17)
x
= [x 0T , x 1T , . . . , x NT ]T ;
u E(x, u)
= [u 0T , u 1T , . . . , u TN−1 ]T ; p = max1≤i≤ p {ωi Ji (x, u)} or i=1 ωi Ji (x, u) p with ωi ≥ 0 and i=1 ωi = 1;
Solving the equality-constrained minimization problem (9) is equivalent to solving the associated system of differential equations (14)–(18).
x(0)
given initial condition.
A. Decomposition–Coordination Method
Note that even if all functions Ji (x, u) are differentiable the cost function max1≤i≤ p {ωi Ji (x, u)} will have corners where differentiability fails; i.e., the cost function max1≤i≤ p {ωi Ji (x, u)} has discontinuous partial first-order derivatives at points where two or more functions Ji (x, u) are equal to max1≤i≤ p {ωi Ji (x, u)} even if the Ji (x, u) have continuous first-order partial derivatives. For this problem, we can construct the ordinary Lagrange function L=
N−1
Lk
(10)
k=0
where 1 E(x, u) + μ0T ( f (x 0 , u 0 ) − z 0 ) (11) N 1 L k = E(x, u) + μkT ( f (x k , u k ) − z k ) + βkT (x k − z k−1 ) N for 1 ≤ k ≤ N − 2 (12) 1 L N−1 = E(x, u) + μTN−1 ( f (x N−1 , u N−1 ) − x d ) N T (x N−1 − z N−2 ) (13) + β N−1 L0 =
μk (n components) and βk (n components) are the Lagrange multiplier vectors introduced to consider the equality constraints (7) and (8). The derivations of the ordinary Lagrange function (10) enable us to transform the equality-constrained minimization problem (9) into a set of differential equations on the basis of which we design ANN architectures. A stationary (equilibrium) point (x k∗ , u ∗k , μ∗k , βk∗ , z k∗ ), according to the Karush-Khun-Tucker conditions [26],
∇βk L =
1 ≤ k ≤ N − 1.
(18)
One way to solve the equality constrained minimization problem (9) is to decompose the treatment of the associated system of differential equations (14)–(18) between two levels, not arbitrarily but in such a way as to obtain a separable form at the inferior level, i.e., each subproblem k only brings about the intervention of unknown variables of indices k, for 0 ≤ k ≤ N − 1. In effect, the processing of the system of (14)–(18) is divided between two levels. The upper level handles (16) and (18) and fixes z k (0 ≤ k ≤ N − 2) and βk (1 ≤ k ≤ N − 1), which it proposes to the lower level. This allows the decomposition of problem (9) to the lower level; each subproblem becomes as follows. 1) Subproblem 0 ⎧ z being given by the upper level ⎪ ⎨ 0 min E(x, u) (19) ⎪ ⎩s.t. f (x , u ) = z . 0 0 0 2) Subproblem k ⎧ z (k = 1, 2, . . . , N − 2) and βk (k = 1, 2, . . . , N − 2) ⎪ ⎪ k ⎪ ⎨being given by the upper level (20) ⎪min E(x, u) + βkT x k − c, where c = βkT z k−1 ⎪ ⎪ ⎩ s.t. z k = f (x k , u k ). 3) Subproblem N − 1 ⎧ z N−2 and β N−1 being given by the upper level ⎪ ⎪ ⎪ ⎨min E(x, u) + β T x N−1 N−1 − c (21) T ⎪ where c = β N−1 z N−2 ⎪ ⎪ ⎩ s.t. f (x N−1 , u N−1 ) − x d = 0. The solution of each subproblem thus corresponds to the processing of (14), (15), and (17) for z k (k = 0, 1, . . . , N − 2)
2504
IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. 26, NO. 10, OCTOBER 2015
and βk (k = 1, . . . , N − 1) supplied by the upper level. Thus, applying a gradient method, we obtain the system of differential equations
system of difference equations (28)–(30) and determination of ( j) ( j) ( j) ( j) ( j) ( j) the variables x k∗ (z k , βk ), u ∗k (z k , βk ), and μ∗k (z k , βk ), which, respectively, satisfy (14), (15), and (17). The results
d xk = −λx ∇xk L (22) dt du k = −λu ∇u k L (23) dt dμk = λμ ∇μk L (24) dt where λx > 0, λu > 0, and λμ > 0. Equations (22)–(24) can be written in the following scalar form:
n d x kq 1 ∂E ∂ fi = −λx + μki + βkq dt N ∂ x kq ∂ x kq
x k∗ (z k , βk ) and μ∗k (z k , βk ) are supplied for the upper level, which verifies whether the previously supplied information was correct and corrects it if necessary, so that the lower level can recommence its work with information of increased validity. Such correction is necessary to bring about evolution toward the satisfaction of (16) and (18). The upper level proceeds progressively in making the necessary correction to the coordination parameters, coming closer and closer to the satisfaction of coordination equations (16) and (18). Thus, the coordination parameters at iteration j + 1 are improvements on the coordination parameters at iteration j ( j ) ( j ) ( j +1) ( j) ( j) = z kq − λz −μ∗iq z k , βk − βk+1,q z kq
i=1
du kq dt
( j)
k = 1, . . . , N − 1 and q = 1, . . . , n
(26)
k = 1, 2, . . . , N − 1 and q = 1, 2, . . . , n (32)
(27)
where λz > 0 and λβ > 0. The solution of the system of difference equations (28)–(30) is thus repeated until satisfactory coordination is obtained, i.e., satisfaction of coordination equations (16) and (18). The method thus defined is shown in Fig. 3.
( j +1)
βkq
i=1
i=1
k = 1, . . . , N − 1 and q = 1, . . . , n
n (l) 1 ∂ E (l) ∂ f i (l) (l) = u kq − λu + μ N ∂u kq ∂u kq ki
(28)
k = 1, . . . , N − 1 and q = 1, . . . , m (l) (l) ( j ) = μ(l) − z kq , kq + λμ f q x k , u k
(29)
k = 0, . . . , N − 1 and q = 1, . . . , n.
(30)
i=1
μ(l+1) kq
( j)
(25)
where f (x k , u k ) = ( f 1 (x k , u k ), f 2 (x k , u k ), . . . , fn (x k , u k ))T , x k = (x k1 , . . . , x kn )T , u k = (u k1 , . . . , u kn )T , μk = (μk1 , . . . , μkn )T , βk = (βk1 , . . . , βkn )T , and z k = (z k1 , . . . , z kn )T . To realize the discrete-time network, we can transform the system of differential equations (25)–(27) into an appropriate system of difference equations. For example, using the forward Euler rule, the system of differential equations (25)–(27) can be converted into the system of difference equations
n (l) 1 ∂ E (l) ∂ f i ( j) (l+1) (l) + μki + βkq x kq = x kq − λx N ∂ x kq ∂ x kq
(l+1)
( j)
k = 1, . . . , N − 1 and q = 1, . . . , n
n 1 ∂E ∂ fi = −λu + μki N ∂u kq ∂u kq
dμkq = λμ ( f q (x k , u k ) − z kq ) dt k = 0, . . . , N − 1 and q = 1, . . . , n
u kq
( j)
Equations (28)–(30) can only, then, be solved locally if it ( j) possesses the necessary information, βk (k = 1, . . . , N − 1) ( j) and z k (k = 0, . . . , N − 2) sent from the upper level. To transmit the information necessary to the functioning of the lower level with a view to overall optimization [i.e., satisfaction of all equations (14)–(18)], it is essential to coordinate the two levels. To ensure such coordination, the ( j) upper level works simultaneously on βk (k = 1, . . . , N − 1) ( j) and z k (k = 0, . . . , N − 2), which make up the coordination parameters. These coordination parameters are considered as known within the lower level, allowing local resolution of the
k = 0, 1, . . . , N − 2 and q = 1, 2, . . . , n (31) ( j) ( j) ( j ) ( j) ∗ z k , βk − z k−1,q = βkq + λβ x kq
B. Analysis Stability The aim of this section is to demonstrate that the convergence of the method described by the diagram in Fig. 3 may be reduced to that of the coordinating level. To facilitate this paper, we will introduce the following notations. x k and call v k∗ (z k∗ , βk∗ ), μ∗k (z k∗ , βk∗ ), Let us posit v k = uk ( j) ( j) z k∗ , and βk∗ the solution sought. Let us also call v k∗ (z k , βk ) ( j)
( j)
and μ∗k (z k , βk ), the variables to be determined at the lower level by the system of equations (28)–(30) so that (14), (15), ( j) ( j) and (17) may be locally satisfied, and z k and βk , the coordination variables to be modified at the upper level following algorithms of coordination (31) and (32). Let us put the problem, which interests us in a more condensed form ⎧ Gk (v k , μk , βk ) = 0 ⎪ ⎪ ⎪ lower level equations ⎪ ⎨ Pk (v k , z k ) = 0 (33) Rk (v k , z k−1 ) = 0 ⎪ ⎪ ⎪ ⎪ ⎩ Hk (μk , βk+1 ) = 0 upper level equations where Gk =
L xk , Pk = μk L, Rk = βk L, and u k L
Hk = zk L. For the solution sought, we may write ⎧ ∗ G =0 ⎪ ⎪ ⎨ ∗k Pk = 0 R∗ = 0 ⎪ ⎪ ⎩ k∗ Hk = 0
(34)
MESTARI et al.: SOLVING NECMOPs USING NEURAL NETWORKS
Fig. 3.
2505
Information transfer between lower and upper levels.
where Gk∗ = Gk (v k∗ (βk∗ , z k∗ ), μ∗k (βk∗ , z k∗ ), βk∗ ), Pk∗ = ∗ ), and H∗ = Pk (v k∗ (βk∗ , z k∗ ), z k∗ ), Rk∗ = Rk (v k∗ (βk∗ , z k∗ ), z k−1 k ∗ ∗ ∗ ∗ Hk (μk (βk , z k ), βk+1 ). At each iteration j of the coordination loop, we have ( j) Gk = 0 (35) ( j) Pk = 0 ( j)
( j)
( j)
( j)
( j)
( j)
= Gk (v k∗ (βk , z k ), μ∗k (βk , z k ), βk ), and where Gk ( j) ( j) ( j) ( j) Pk = Pk (v k∗ (βk , z k ), z k ). ( j) ( j) ( j) = v k∗ (βk , z k ) − v k∗ (βk∗ , z k∗ ), Let us posit ev k ( j) ( j ) ( j ) ( j) ( j) eμk = μ∗k (βk , z k ) − μ∗k (βk∗ , z k∗ ), ezk = z k − z k∗ , and ( j)
( j)
eβk = βk − βk∗ , which designate the errors calculated at iteration j of the coordination loop.
( j)
( j)
( j)
( j)
Theorem 1: Let ev k , eμk , ezk , and eβk be the errors computed at the iteration j of the coordination loop. Then, ( j) ( j) ( j) ( j) ev k −→ 0 and eμk −→ 0 if ezk −→ 0 and eβk −→ 0. Proof: Linearization of the lower level equations in the neighborhood of solution allows us to write ( j)
Gk∗ +
∂Gk∗ ( j ) ∂Gk∗ ( j ) ∂Gk∗ ( j ) e + e + e ∂v k v k ∂μk μk ∂βk βk
(36)
( j)
Pk∗ +
∂Pk∗ ( j ) ∂Pk∗ ( j ) e + ez ∂v k v k ∂z k k
(37)
Gk Pk
where considering (34) and (35), (36) and (37) become ∂Gk∗ ( j ) ∂Gk∗ ( j ) ∂Gk∗ ( j ) e + e + e =0 ∂v k v k ∂μk μk ∂βk βk
(38)
2506
IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. 26, NO. 10, OCTOBER 2015
and
Using (45) and (46) in (42), we have ∂Pk∗ ( j ) ∂Pk∗ ( j ) e + ez = 0. ∂v k v k ∂z k k
(39)
( j ) =
k=0
( j)
Equations (38) and (39) demonstrate that ev k −→ 0 and ( j) eμk
−→ 0 if
( j) ez k
−→ 0 and
( j) eβ k
−→ 0.
N−1
∗ ∂Rk∗ ( j ) ( j ) T ∂Rk ( j ) ev k + eβ k ez ∂v k ∂z k−1 k−1 1 ( j )T ( j ) ( j )T ( j ) + λ2 Hk Hk + Rk Rk 2 ( j )T
+ eβ k
♦
( j)
( j)
One need only study, then, the variation of ezk and eβk during the iterations of the coordination loop to determine the convergence of the method illustrated by the diagram in Fig. 3. Given the very major role played by coordination (the upper level’s task), a more in-depth study of the convergence (stability) of coordination algorithms (31) and (32) becomes necessary. Let us choose as Lyapunov function
∗ ∗ ( j ) T ∂Hk ( j ) ( j ) T ∂Hk ( j ) λ −ezk eμk − ezk e ∂μk ∂βk+1 βk+1
= A( j )λ2 + B( j )λ
(47)
where A( j ) =
N−1
( j )T
Hk
( j)
( j )T
H k + Rk
( j)
Rk ≥ 0
k=0
and ( j ) =
N−1
1 2
( j )T ( j ) ez k
ez k
( j )T ( j )
+ eβk eβk ≥ 0.
(40)
k=0
( j ) = ( j + 1) − ( j ) N−1 1 ( j +1)T ( j +1) ( j )T ( j ) ez k − ez k ez k = ez k 2 i=0
( j +1) T ( j +1) eβ k
+ eβ k
( j )T ( j )
− eβ k eβ k
.
(41)
The development of (41) allows us to write ( j ) =
( j )T
ez k
k=0
+
( j)
( j )T
( j)
ezk + eβk eβk
1 ( j )T ( j ) ( j )T ( j) ezk ezk + eβk eβk (42) 2
where ( j)
( j +1)
− ezk = −λHk
( j)
( j +1)
− eβk = λFk
ezk = ezk
eβk = eβk
( j) ( j)
( j)
( j)
( j)
Hk +
( j)
Rk∗ +
Rk
−ezk
∗ ∂Hk∗ ( j ) ( j ) T ∂Hk ( j ) eμk − ezk e ∂μk ∂βk+1 βk+1
( j )T + eβ k
∗ ∂Rk∗ ( j ) ( j ) T ∂Rk ( j ) . e + eβ k ez ∂v k v k ∂z k−1 k−1
( j )T
(∗)
∂Hk∗ ( j ) ∂Hk∗ ( j ) eμk + e ∂μk ∂βk+1 βk+1
∂Rk∗ ( j ) ∂Rk∗ ( j ) ev k + ez . ∂v k ∂z k−1 k−1
Considering (34), (43), and (44), we may write ∗ ∂Hk ( j ) ∂Hk∗ ( j ) ( j) ezk = −λ eμk + eβk+1 ∂μk ∂βk+1 ∗ ∂Rk ( j ) ∂Rk∗ ( j ) ( j) eβk = λ ev k + ezk−1 . ∂v k ∂z k−1
Theorem 2: Let λ be the adaptive coefficient for the coordinative algorithms. Then, the convergence is guaranteed if one of the matrices (∂Gk∗ /∂v k )T (k = 0, 1, . . . , N − 1) is positive definite and the others are only positive semidefinite, and furthermore, if A( j ) = 0, λ must be chosen as: 0 < λ 0, the equation ( j ) = 0 admits two real distinct roots 0 and −(B( j )/A( j )). Then, ( j ) < 0 if B( j ) < 0 (i.e., exactly one of the matrices
MESTARI et al.: SOLVING NECMOPs USING NEURAL NETWORKS
2519
(∂Gk∗ /∂v k )T (k = 0, 1, . . . , N − 1) is definite positive and the others are only positive semidefinite (78) and λ is chosen as B( j ) B( j ) . (79) 0
