- Email: [email protected]
Solving degenerate optimization problems using networks of neural oscillators
Neural Networks, Vol. 5, pp, 949-959, 1992
0893-6080/92 $5.00 + .00 Copyright © 1992 Pergamon Press Ltd.
Printed in the USA. All rights reserved.
ORIGINAL C O N T R I B U T I O N
Solving Degenerate Optimization Problems Using Networks of Neural Oscillators D E R E K M . WELLS Florida Atlantic University
( Ra'eived 1 October 1990: revised and acct7~ted 27 Felwuao' 1992)
Abstract--Optimization problems containing multiple, identical solutions, in particular the travelling salesman problem (TSP), are soh,ed using networks of neural oscillators. Processing units are described by two state variables representing fast attd shin' membrane events, much like the Fitzhugh-Nagumo model neurons, and are passive oseillator~ in that at least two processing units with mutually inhihitoo' couplings are required to generate sustained oscillations. Tile emergent behavior of the network is to produce stable travelling waves (limit ~3'cles) o.[ftring rate impulses in the position coordinate of the network with the phase relationship between active (suprathreshold) processing units determining the TSP solution. This approach has the advantage of avoiding most h~cal energy minima with tile best or near-best sohttian being finmd eveo' time.6~r small net works ( < 10-~processing units). ,,Is network size is enlargened, the time required to reach a stable travelling wave can./igaration increases. To improve tile convergence rate, the strength of the network couplings encoding the cost./imction and tile strength q[the network couplings encoding tile syntactic constraints are altenlatel)' reduced slich that net work evohttion is dictated by only one sel o[eouplings at an)' given time. Good solutions are fimnd./hr net works containing up to 104 processing units without a large ~:vpenditure ~[computational t:[]brt.
Keywords--Neural oscillators. Optimization, TSP. 1. INTRODUCTION
act out-of-phase or at a different frequency) to the dynamics of the cooperative networks would be minimal. There is a large body of literature that deals with modelling coordinated movement in lamprey using nonlinear coupled oscillators (see Cohen, Holmes, & Rand, 1982; Kopell & Ermentrout, 1988) where the coupling strength between oscillators is proportional to a phase difference between body segments. The model neuron in these approaches is based upon the FitzhughNagumo (FN) neuron or variants of the FN neuron (Fitzhugh, 1961; Nagumo, Arimoto, & Yoshizawa, 1962). Abbot (1990) has studied some of the properties of FN networks and has shown that associative memory can be stored as a phase relationship between network elements, different memory patterns having different phase relationships. Perhaps the most fully understood neuronal circuitry that generates rhythmic motor patterns is the stomatograstric system of the lobster (see Selverston, 1988, for a review). The stomatogastric ganglion is composed of two networks which control the muscles of the gastric mill and pyloric regions of the lobster stomach. Selverston has shown that the pyloric network consists of neurons that can be classified as intrinsic oscillators (i.e., neurons that are capable of oscillating individually) whereas most of the neurons in the gastric network can be classified as passive os-
Many of the artificial neural networks that have been developed do not take into consideration temporal factors associated with real biological neural networks. In many biological systems, it is often the underlying neurodynamics that is the most important factor to understand with regard to coordination ( Kelso & Sch6ner, 1988; Sch6ner & Kelso, 1988) and movement in general. Based upon studies of the visual cortex in cats, Gray and Singer ( 1989 ) and Sporns, Gally, Reeke, and Edelman (1989) have suggested that network oscillatory responses (namely, phase and frequency locking) could act as a general mechanism for coordinating spatially distinct populations of cortical cells involved in sensorimotor tasks. The influence of additional inputs (that
Acknowledgements: I would like to thank Scott Kelso, Hirofumi Nagashino, Abhijit Pandya, Tim Kiemel, Gonzales DeGuzmann, and Pier Zanone for commenting on earlier versions of this manuscript and Ed Page for introducing me to the topic of network representations of optimization problems. This research was supported in part by a grant from the Otfice of Naval Research ( N00014-88-J- 1191 ) to JAS Kelso. Requests for reprints should be sent to Derek M. Wells. Center for Complex Systems. Florida Atlantic University. Boca Raton, FL 33431.
949
950
cillators (i.e., neurons that do not oscillate individually but only within a network structure). One approach to modeling passive oscillators has been described by Matsuoka ( 1985, 1987 ). A fast rise-time constant and a slow decay-time constant are associated with each neuron such that a nonmonotonic response curve to a step input develops. Sustained oscillations can be generated when two or more of these model neurons have mutually inhibitory connections. A useful method to study emergent behavior within networks of model neurons is to have the networks solve well defined yet computationally difficult problems. These problems often involve competing sets of constraints such as in optimization problems. The classic biological example of a behavior involving competing sets of constraints is the movement of the arm towards some target. For this motor task, a speed-accuracy tradeoff exists (which was parameterized by Fitts in 1954) wherein accuracy of the final hand position with respect to the target coordinates is compromised at high speeds, Most well defined problems however are distinctly nonbiological, although the network coupling architectures used to solve the given problem as well as the emergent network properties might have some relevance to information processing within real neural subsystems. In this paper, a network of passive oscillators is used to solve an optimization problem containing multiple, identical solutions in particular, the travelling salesman problem (TSP). This problem consists of finding the minimum tour distance through a set of nodes where each node is visited exactly once and where the tour returns to its starting location. For any tour through a set ofn nodes the solution is 2n -fold degenerate in that there are n possible starting positions and two possible tour directions that yield the same tour length. A neural network solution to this problem, which involves encoding a cost function (to minimize tour distance ) and a set of syntactic constraints (to uniquely assign each node a position in the tour) within the coupling matrix of the network, has been pursued by a number of researchers (Hopfield & Tank, 1985; Kamgar-Parsi & Kamgar-Parsi, 1990; Wilson & Pawley, 1988). Each group has reported that TSP scales poorly as the size of the network increases and that the success rate for finding valid tours through the nodes is at best modest. Here, the emergent behavior of the passive oscillator network is to develop travelling waves of firing rate impulses in the position coordinate of the network. The phase relationship between active (suprathreshold) processing unit impulses at a given instant in time dictates a solution to the optimization problem. Individual oscillators within the network exhibit frequency locking (i.e., the variability in individual oscillator frequencies tends toward zero as the system evolves) but do not exhibit phase locking (i.e., the midpoints of the firing rate impulses for each subset of oscillators that represent
D. M. Wells
a valid tour do not synchronize as the system evolves). A measure of the system stability is defined as the variation in this impulse phase offset over time. As network size is increased, the time (number of oscillator periods) required to reach a stable travelling wave configuration increases. To improve the convergence rate for large networks, the strength of the network couplings encoding the cost function and the strength of the network couplings encoding the syntactic constraints are alternately reduced such that network evolution is dictated by only one set of couplings at any given time. Good TSP solutions are found for networks containing up to 104 processing units without a large expenditure of computational effort. 2. T H E M O D E L N E U R O N AND N E T W O R K There are many different models of individual biological neurons which vary in their level of description of biophysical events. A mathematically tractable set of differential equations which are based on the HodgkinHuxley equations was developed independently by Fitzhugh ( 1961 ) and Nagumo et al. (1962). The differential equations describing the model neuron consist of two state variables ( U and V): U is the slow recovery component of the membrane current and V is the transmembrane potential. The Fitzhugh-Nagumo equations have the form: dV
C - ~ = - G ( V ) - U + ~] li,
(1)
i
L -~t dU = - R U + V - 0,
(2)
where C is the membrane capitance, G(V) is a function which mimics the N-shaped relationship between the fast membrane current and the cell potential V and usually is chosen to be a cubic polynomial: G ( V ) = a(I" + pV 2 + qV 3) where a, p, and q are constants that are model specific, R and L are the resistance and inductance associated with the slow current events, T I is the sum of all currents entering the cell both from other neurons in the network and from sources external to the network, and 0 is the membrane potential threshold. Electronically equivalent circuits of networks of these model neurons have been constructed using tunnel diodes (Hoppenstaedt, 1981 ) or sets of operational amplifiers (Keener, 1983). The FN equations have been used to model individual action potentials and at a different level of description (i.e., a slower time scale) to model impulses of neural activity where the output signal is considered to be an average firing rate. These model neurons exhibit rich dynamic properties and are able to mimic a number of the characteristics of real biological neurons including post-inhibitory rebound, plateau potentials, and frequency entrainment (Abbot, 1990).
Solving Degenerate Optimization Problems
951
The differential equations used to represent the model neurons in this paper are similar to the FN equations except that the fast current-voltage curve G ( V ) is some increasing monotonic function instead of being an N-shaped function. In the simulations which follow, G(V) was chosen to be a linear function, G(V) = a V, which greatly reduces the dynamic properties of an individual model neuron: d V / d t = d U / d t = 0 for large t. The neuron output response is modelled as an averaged firing rate which is related to transmembrane potential by a function F(V) that is usually chosen to be sigmoidal in shape. The equations for the state variables of the model neurons ( U i, Vi, i = 1. . . . . N) have the form: dVi dt
= - a I ~ - U, + I, + dUi dt
=
F(V, )
-{3(u,
-
y = 0.01
y
=
0,05
y=0.1
Z kL
y =0.2
y =0.5
N ~ woF(Vj), j=Lj,,i
(3)
time
v v , ),
(4)
FIGURE 1. Transient response of an isolated model neuron to a step input for a range of 3, values (with a/B = 10).
1 -
l +
y = 0.001
e -x(v'-°')
"
5)
where wo is the coupling strength from neuron j to neuron i, [3 = R L -t , 3' = R - ~and ~ is a proportionality constant which specifies the steepness of F(V); C was set equal to one which does not reduce the generality of the equations. The influence of the parameter 3' (with a/[3 = I0) on the response characteristics of a single model neuron to a step input ( I = 0 to I =/,,) is shown in Figure 1. In the simulations which follow, 3' was chosen such that the transient response of a model neuron is overdamped: 3' < ( a - [3)2/4[3. Matsuoka (1985, 1987) has used a similar set of equations as a model for the phenomenon of accommodation (the decrease in firing rate over time to a constant input) and as an approach to generating sustained oscillations in model networks. In the limit as 3' tends to zero, U = Uoe -~' also tends toward zero for large t hence, eqns ( 3 ) to ( 5 ) reduce to the Hopfield differential equations (Hopfield, 1982, 1984). Results presented in this paper were also reproduced using a piece-wise linear version of F(V) with no change in the qualitative features of the network behavior. The magnitude of the external input current, I~, and the strength and architecture of the processing unit couplings, ~'0, determines the network response.
F(I') =
To understand why sustained oscillations can be generated when two or more of these model neurons are coupled together, consider a network of two neurons with mutually inhibitory connections ( w~2 = w2~ = - 77) and constant external input current (I, = Iz = Io) as shown in Figure 2a. To facilitate analysis let F(V) be piecewise linear with slope 1/6 (0 = 0):
0< z'
1,
I ' > t3
O,
I'<0
(6)
and let a ,> /3 so that the system acts as a relaxation oscillator ( d V / d t ,> d U / d t ) . In this case, the time spent in the firing rate transition region (0 < F ( V ) < 1 ) is far less than the time spent in the saturation regions ( F ( V ) = O, F(V) = I) which permits the system dynamics to be approximated by only two of the four equations of motion at any given time. That is, if F( l~ ) equals 0 or 1 then Ui governs the overall system dynamics ( f, changes so rapidly with respect to Ui in this region that it can be considered constant for any value of Ui ). From eqn (3), d l ~ / d t = 0 gives, 1
I', = - ( - U ,
+ I , , - ~F(I,))),
(7)
and substituting eqn (7) into eqn (4) yields, dt
~ I +
U, + ~ ~7 (I,, - ~F(1))). -
(8)
Whereas. if0 < F(I~ ) < 1 then I~ governs the overall system dynamics ( Ui = U,,. = constant): dI',
2.1. Two Coupled Oscillators
I" y,
dt
-
c~I,} -
U,c +
1,,-
rlF(I~).
(9)
The usual approach for determining if a given system is capable of sustained oscillations is to calculate where the fixed points for the system ( d U / d t = d V / d t = O) lie with regard to the current state of the system; if stable fixed points exist in the same firing rate region as the current system state then sustained oscillations are not possible (see Grasman, 1987, for an overview
952
D. M. Wells I.
Io (w~: = w n < 0 )
(a)
u~ ',',
.:..÷.:,~.,-,.'. ,.:.
~,,.,,.:, ~.,.:.L.,,,.,. '-?,; ::.,
a
i ' _i . oi _~ o _
I.o-,~.,°-~,
Ui
(b) FIGURE 2. A two-oscillator network: (a) coupling architecture, (b) state-space trajectory (as denoted by arrows) projected onto the (U1, U2)-plane. There are nine regions corresponding to different values of F(V,) and F(V2): F(VI ) = O, F(V~ ) = 1, 0 < F(VI ) < 1 (i = 1, 2) where within each region limits for the system state variables were determined (e.g., in region I where F(V~) = 1 and F(V2) = 0, V 1 > ~, V 2 < 0, U1 < Io - / ; a and U2 > Io- ~). For the purpose of clarity, only regions 1, 3, and 8 (projected onto the (U~, U2)-plane) are shaded. Different fixed points for the system (dU~/dt = dVddt = 0) are associated with each region (e.g., the fixed point (projected onto the (U~, U2)-plane) for the system in region I is denoted by "b," whereas the fixed point associated with region 5 is denoted by "a"). Sustained oscillations were possible because the fixed points for the system in any region were either outside of the region the system state was in or were unstable.
of relaxation oscillators). Note that different fixed points are associated with different firing rate regions. For this two-neuron network, the firing rate space (F(Vt), F(V2)) can be partitioned into nine distinct regions ( F( V/ ) = 0, F ( V i ) = 1 , 0 < F ( V ~ ) < l , i = I, 2) where within each region the limits for each of the state variables can be established using eqns ( 6 ) to (9). The limits for the slowly changing variables Ut and U2 in each of the nine firing rate regions are shown in Figure 2b. For example, the parallelogram in the center of the (Uj, U2)-plane indicates the region where both F( Vt ) and F(1/"_,)are in transition. Under the conditions that ~ > 6a and 77 < Io < 6(a + ~), it can be shown that the fixed points associated with any firing rate region of the system either lie outside of the region the system state is in or are unstable. For this system, a
zone of bistability exists within the square that is centered in the parallelogram; the fixed points associated with regions 1 and 5 are stable but lie outside of this zone, whereas the fixed point associated with region 9 is located in this zone but is unstable (a saddle node point). A trajectory of the system state projected onto the ( U~, U2)-plane is denoted by the arrows in Figure 2b, in which sustained movement occurs between firing rate regions 1 and 5. 2.2. Network Representation of T S P In order to solve an optimization problem such as TSP using a neural network, an appropriate representation of the problem must be formulated. This involves choosing a suitable network topology as well as defining external input currents ( l i ) and coupling strengths between processing units (w,i). One such network mapping scheme for TSP has been described by Hopfield and Tank (1985) in which they represented a valid tour through a set o f n nodes by an n × n-permutation matrix; an element (i, j ) in the matrix represents the probability that node i will be in position j of the tour. An example of a valid tour through a set of five nodes is illustrated in Table I. The total tour distance can be decoded from this representation by taking the following sum: d = d.4c + dcB + dBE + dEo + doA, where the order of the n nodes in the tour are ACBEDA. A valid tour requires that exactly one element in each row and column of the matrix be set to "1," all other elements equaling zero hence, a total of n! valid tours are possible. However, only (n - l ) ! / 2 of the tours are distinct in that for any valid tour through the nodes there are n possible starting positions and two possible directions which yield the same tour length. Based on this representation scheme for valid TSP solutions, Hopfield and Tank (1985) developed a neural network in which each processing unit in the network corresponded to an element in the permutation matrix. For 3' = 0, Hopfield (1984) has shown that the equations of motion for networks with certain coupling architectures ( wij = ~t~,, Wig= 0) lead to stable fixed point configurations o f F ( V, ). A computational energy function that minimizes total tour length while simultaneously satisfying the row and column constraints can TABLE 1 The Permutation Matrix Representing the 5-node Tour ACBEDA
Position Node
1
2
3
4
5
A
I
0
B
0
0
0
0
0
I
0
C D E
0 0 0
I 0 0
0
0 0 0
0 0 I
0 I 0
Solving Degenerate Optimization Problems
953
be associated with the network. This function takes the general form: 1 N
E
= -
~
~
ill
Y÷X
dzF-J(z),
(10)
i=1
and is related to the equation of motion for a processing unit (in the special case where 3, = 0) by dV~/dt =
-OE/OF( V~ ). In order to more clearly show how the constraints and cost function for TSP are mapped into the network couplings, matrix notation will be used to specify the firing rate of an individual processing unit with the indices X and Y referring to nodes (rows of the network) and the indices i a n d j referring to positions in the tour (columns of the network); for example, F(Vxi) is the firing rate for the processing unit in row X and column i of the network. The coupling term, E,., of an energy function that is minimal when the row and column constraints of TSP are satisfied may now be specified: 1,I,"t
Z Z ~ F(Vx,)F(Vxj) X
i
j÷i
+-g~ ~ ~ ~ F(I'v,)F(I,'r,), --
i
(11)
X Y ÷X
where W~ is the coupling strength corresponding to the row and column constraints and all summations are from 1 to n (the number of nodes). Often, another term proportional to
#
E,. = T ~x ~i Z F(Vxi)F(I/xj) )÷i
W2
+ ~ E Z Z F(VxJF(VvJ +--~ Z Z Z 4rvF(Vx,) --
i
X
}'÷X
X"
Y+X
X(F(Vya+,)F(I/y,,_,)), +
(13)
where W2 is the coupling strength corresponding to the cost function, dvv is the distance between nodes X and Y and subscripts for the position coordinate of the network are expressed modulo n (i.e., F(Vr.,+j) = F(Vr.j). By taking the derivative - 0 E,./OF(Vxi), the coupling term in the equation of motion for each processing unit in the TSP network is determined:
(14)
Y÷X
The local neighborhood for a processing unit in a network representing a 5-node TSP is illustrated in Figure 3. Network solutions to TSP are possible when appropriate choices for the coupling strengths, Wt and W2, are made. The differential equations governing the network are summarized below: dl':v' -
dt
a V x , - Ux, + Io- l'Vt( ~ F(Vxj) + ~ F(Vr,) )
j,,,
r,,x
-dxvW2(~F(Vr.,+,)+ ~ F ( V r . H ) ) , Y÷X
(15)
Y÷X
dUxi _ dt
B( Ux; - 7 l';vi),
(16)
1
F(Vx,) - I + e -xv-v' "
(17)
In the case where 3' > 0, the processing units as defined by eqns ( 15 ) to ( 17 ) do not settle into a stable fixed point configuration in the F(V) space; instead they form limit cycles in the position coordinate of the network. Stability, in the context of this network of oscillators, is defined as the formation of travelling Position 1
2
3
,,2)
is included in E,. that specifies the total neural activity of the network. However, it can be shown (Tagliarini & Page, 1987) that this global coupling term is not required if the row and column coupling strengths (Wt) and the external input current (Ix~ = I,,) are suitably chosen. A more complete coupling term of an energy function that also favors short tours has the following form:
14 "1
~ F(Vri)) Y÷X
-dx,.l.t/'2(~_~F(Vr.i+,)+ ~F(Vy.i_l) ) .
woF(V~)F(Vj) - Z F(V, )I, +a
=T
j÷i
N
Z t=l= j = t
E,.
-H.'I(~F(I.~)+
o
4
QQ
5
@®OO ®0
00
FIGURE 3. The local coupling neighborhood of processing unit B2. Processing units that are symbolized by circles containing squares are negatively coupled to unit B2. Circles with open squares encode the row and column constraints whereas circles with filled squares encode the cost function which minimizes total tour distance. Note that columns 1 and 5 are considered adjacent when applying the cost function so that no edge effects occur. This particular coupling architecture allows travelling waves of firing rate impulses to form in the position coordinate of the network.
954
D. M. Wells
waves (limit cycles) of firing rate impulses that maintain a constant phase relationship. The set of processing units that are simultaneously active (F(V) > threshold) form a solution to TSP. Returning briefly to the 5-node tour illustrated in Table 1, a stable travelling wave of firing rate impulses results in the order of nodes for valid solutions changing from ACBEDA to CBEDAC to BEDACB, and so forth. For a given simulation, a travelling wave in the opposite direction (ACBEDA to DACBED to EDACBE, and so forth) is equally probable although, the direction of the travelling waves can be prespecified if the coupling term corresponding to the cost function in eqn (15) is modified: -dxrW2((l +~) ~ F(V,.,+,)+ ~ F(Vv,,-,)), (18) Y~X
Y~'X
where I~l ~ 1. It should be noted, however, that ~ does not eliminate the 2-fold degeneracy in tour direction, it simply improves the convergence rate to stable limit cycles. For the simulations which follow ~ was set equal to zero. A time period, T, can be associated with the travelling wave and is defined as the number of system iterations required to change the order of suprathreshold processing units from ACBEDA to CBEDAC to BEDACB, and so forth. In the following section T is used to scale the evolution of the network state. 3. S I M U L A T I O N RESULTS Simulations were performed for networks ranging in size from 102 processing units to 10 4 processing units. For small networks, the set of parameter values used in the simulations was: a -- 0.1,/3 = 0.01, 3, = 0.1, X = 0.035, Wt = 12.0, W2 = 8.0, and Io = 20.0. Distances between individual nodes (dxr) were scaled such that an average separation distance of 0.5 was achieved. For large networks, the ratio W~:W2 was systematically varied over the time course of the simulation (more details about this schedule are described in Section 3.2). The magnitude of the external input current, I,,, determines the number of system iterations a processing unit remains suprathreshold. That is, Io dictates the width of the firing rate impulse. I,, was set such that the impulse width associated with each processing unit at one-half the maximum firing rate (F(V) = 0.5 ), the so-called full width at half maximum value ( F W H M ) , was 1.5 times the magnitude of T. This value for I,, was required in order to generate stable travelling waves through the network and was determined empirically. The corresponding size of T was determined experimentally using a simple peak finding routine which measured the midpoint of the time a given processing unit was suprathreshold. Valid solutions to TSP were determined by detecting which processing units in the network were suprathreshold (F(V) > Fo) at time t, where Fo is a detection threshold level that is limited to the range of F(V). To ensure that solutions to TSP
are possible, the value for Fo must be chosen such that the width of the firing rate impulse at the corresponding level ofF,, is less than 2T. For Io equal to 20.0, a suitable choice for F,, is 0.75. If exactly one unit was active for each row and column of the network matrix then a valid solution to TSP was recorded. Equations ( 15 ) and ( 16 ) were solved using the Euler method for numerical integration in which each processing unit of the network was updated once (in a random sequence) per system iteration. The initial state of the system was determined as follows: F(Vxi ) = F W H M / n T = 1.5/n (this value was based, approximately, upon the amount of neural activity that is present at any given time after the system has settled into a travelling wave configuration), Vxi was calculated using eqn (5) and Uxi was obtained from the following relation:
Uxi = -c~[~r, + I,, - 2(n - I)(W I + dxrB~)F(VxA. (19) 3.1. Small Networks Emergent behavior within complex systems has been classified by Wolfram (1984) into four general types: (a) fixed point, (b) periodic, (c) chaotic, and (d) intermittent. For a network consisting of the TSP coupling architecture (in the appropriate parameter regime), the emergent behavior of the oscillators was to produce travelling waves with constant period T, type 2 behavior (Figure 4a). Initially, subsets of travelling waves formed in competing directions (waves moving upward versus downward) however, after an elapsed time of about n T, only travelling waves moving downward persisted. Global organization of the network oscillators is not possible for the case where 3' equals zero; fixed point solutions to TSP (which are highly dependent upon initial conditions) are formed at elapsed times of less than n T. This property of the network to produce stable travelling waves was dependent upon the values of W~ and W2. If, for instance, Wz was set equal to zero, oscillatory responses still formed but appeared to be aperiodic (Figure 4b). In addition, intermittent type behavior was possible for intermediary values of Wz. Throughout the rest of this paper, only TSP networks that generate stable travelling waves will be considered. Figure 5 shows typical solutions to 10-node and 20node problems. The best or near-best solution was found with 100% success rate. For some trials, a number of other valid solutions were found (usually with larger tour distances) before the network settled into a stable travelling wave configuration. Coupled oscillators are often characterized by phase and frequency locking. For the TSP network, oscillators representing a given node tended toward a common frequency ( T ) -~ as the system evolved ( Figure 6a) but the phase offsets between oscillator impulses for processing units that were simultaneously active did not
Solving Degenerate Optimization Problems
-~> . 2 ,
:
~,
955
.-),_2. •
~ - " - ~
",
~
~
~
(a)
"-~'-~---)..~
100: Time
(a) I
.
~
~ .
~
L -------r~
~ c __
~ - ~ s
--
~ ~._LA
A
Co) FIGURE 5. TSP network solutions: (a) 10-node, (b) 20-node.
, ~--~.~_
100
~
"
tend toward zero. Figure 6b shows the standard deviation in impulse phase, SD, for 6 trials of the same 10node problem where only the processing unit update schedules varied. SD was calculated as follows:
~
Time
SD (b) FIGURE 4. Time series for each processing unit firing rate, F(V), in a 100 unit TSP network where outputs that correspond to the 10 possible positions in the tour (columns in the network matrix) are displayed consecutively: (a) W2 = 8.0. Sustained firing rate impulses are generated and form travelling waves of constant period, T, in the position coordinate of the network; (b) W2 = 0.0. Sustained firing rate impulses are generated but appear to be aperiodic.
\/Z7-2 (~,,., - ~b,..,)2 n-I V
(20)
where 4~.,. is the midpoint of the m-th firing rate impulse for the active processing unit representing the ith node and (4~,,. - ¢~,,.) is the phase offset between the active processing units representing nodes i and 1. Midpoints corresponding to the m-th firing rate impulse were located within the range: ~b].m -- 0.5T < ~bi.,. < ~b~,,. + 0.5T. If no midpoint was found within this win-
956
D. M. H"ells A measure of system stability (S) was formed by calculating the change in phase offset as the system evolved: / ~ ['=2 (A~bi''' - A4~''-t)2 , n- 1
S = O4 1D 0 0-
0
0
t
I.~ 1
Time / 5 nT (a)
(21)
where A$~.,, = $~.,, - St ..... A plot of S for networks of various sizes is shown in Figure 7 (where/,, was constant: F W H M / T = 1.5). As network size was increased, the time to reach a stable travelling wave configuration increased and at this level of external input current, stability was not obtained for the network representing a 40-node problem. A stable solution might be reached by an elapsed time of 5n T by increasing the external input current but at some level of input, the width of the firing rate impulses (at the level ofF,,) will become greater than 2 T a n d solutions to TSP would no longer be detectable using the approach described in this section. 3.2. Large N e t w o r k s
I-O3
0
1 Time / 5 nT Co)
FIGURE 6. Emergent TSP network properties: (a) the time pedud (number of system iterations) between fidng rate impulses for processing units representing each node (row) of the network tended toward a common period, T; (b) the standard deviation of the impulse phase between active processing units did not tend toward zero as the system evolved. Over multiple trials (six of which are shown), where the resultant tour through the nodes was identical, the standard deviation in the impulse phase ranged from about 0.1T to 0.4T.
dow (which might occur prior to the formation of a stable travelling wave) then 4~,,,, was set equal to 4h .... + 0.5T. For each trial in Figure 6b, the resultant tour through the nodes was the same whereas the range in standard deviation for the impulse phase varied between 0 . I T and OAT, indicating that multiple phase configurations (at a time scale less than T) were stable. From a practical standpoint, the ability to detect valid tours becomes degraded for network oscillators that have nonzero phase offsets.
In order to improve the convergence rate, as well as to solve, large scale TSPs ( up to 100 node problems), the coupling strength ratio 14'j :I,V_, was alternated throughout the simulation such that either the couplings that encoded the row and column constraints or the couplings that encoded the cost function dominated the evolution of the network state, but never both sets of couplings simultaneously, as was the case for small networks where the ratio remained constant (H,'t :H~ = 3:2). Initially, the coupling strength for the row and column constraints ( W~ ) was set to a value of 4.0 ( H"t : I4 2 = 1:2). This lower value for H.'t caused the number of processing units that were simultaneously active in each row and column of the network matrix to increase 1
n=:O . . . . . .
n=20
n=30 n=4O
i i
:,.
,
'i
' .
0
, I
Time 15 nT
FIGURE 7. A plot of the stability measure, S, for networks ranging in s i z e from 100 units to 1,600 units (FWHM/T = 1.5). It should be noted that the time required to reduce S to about zero varied widely for different trials of the s a m e problem.
Soh, ing Degenerate Optimization Problems from 1 to about 5. L Active processing units tended to cluster into either one or two regions in the position coordinate of the network matrix hence, an approximate solution to TSP was formed. However, as the results in this section will show, there was no guarantee that the best solution to TSP was one of the permutations within the approximate solution. The time to reach a stable travelling wave configuration was similar to the time required for small networks: between n T and 2n T. A number of possible methods to solve the TSP given this stable travelling wave configuration can now be implemented. One approach to finding good TSP solutions would simply be to perform a sequential search through all possible combinations of neurons that comprised the approximate solution. If the number of neurons that were active in a given row or column was far less than n (the number of nodes) then the total number of possible tours would be greatly reduced from n! The approach taken here was to increase the coupling strength ratio H'~:HJ2 to 4:1 (W~ = 16.0, W_, = 4.0) which enabled the couplings that encoded the row and column constraints to subsequently dominate the system evolution instead of the couplings that encoded the cost function. In order to favor a final TSP solution that remained within the bounds of the initial approximate solution, the magnitude of the WTcouplings between processing units that were not both simultaneously active in the approximate solution was increased to a value of 8.0. Solutions were found with 100% success rate although the travelling waves became unstable if the system state was allowed to continue evolving given the high coupling strength ratio. By periodically alternating the coupling strength ratio (a time schedule of which is illustrated in Figure 8), multiple solutions to a given TSP were found. Results for the Lin-Kernighan tour (Lin & Kernighan, 1973) were consistently about 4.8 to 5.0 in tour length compared to the best value of 4.32.-' Network solutions to problems involving 50 and 100 nodes (where nodes were randomly distributed in the unit plane) were also readily found: Figure 9 shows typical results. In general, it was found that processing units which comprised the approximate TSP solution clustered into two regions in the position coordinate of the network (i.e., two travelling waves were simultaneously formed) hence, the approximate solution (and consequently, the resultant valid solution) involved two redundant loops through the nodes. The fact that the number of stable travelling waves that simultaneously formed within a network did not conSimultaneously active processingunits in each row and column of the matrix can also be achieved by increasing the external input current, 1,,, or decreasing ,X. 2Approximatecoordinates for the kin-Kernighan tour were provided by Ed Pesulima.
957 Wt : W, 4
3 2 l
0
I 2nT
I ~, t > T
Time
FIGURE 8. Time schedule for alternating the coupling strength ratio Wl :W=. For W1:14/2 = 0.5, stable traveling waves formed after an elapsed time of between n T and 2n T with multiple ( 4 6) active processing units in each row and column of the network matrix; an approximate TSP solution was thereby established. At time 2n T, W~ :14/2was increased to 4.0 so that the row and column couplings dominated the system evolution. In order to favor a valid TSP solution that was one of the permutations of the approximate solution, the coupling strength ratio between processing units that were not both a part of the approximate solution was decreased to 2.0. These higher coupling strength ratios were maintained until a valid solution was detected (as indicated by arrows), after which the ratio was reduced to its original value (0.5) and the process was then repeated. Multiple solutions to a given TSP were thereby found.
tinue to increase beyond two (as network size increased) might reflect some inherent limitation in this network approach to solving TSP. One such limitation is the fact that a stable traveling wave configuration of processing unit activity levels does not eliminate the 2fold degeneracy in tour direction. Clearly, more studies are required in order to better characterize the global organizational properties of these large passive oscillator networks. 4. S U M M A R Y Throughout this paper, some of the properties of networks of passive oscillators have been studied in the context of solving a well-known optimization problem that has multiple solutions of identical cost namely, the travelling salesman problem (TSP). The coupling architecture for the local neighborhood of each processing unit in the oscillator network allowed travelling waves of firing rate impulses to be formed along the position coordinate of the network matrix. For an n node problem in which stable travelling waves formed, a total of n solutions of identical tour length were detected for each pass of the travelling wave through the network. It is believed that by forcing the network to find multiple solutions to TSP better and more consistent results were realized. However, the best solution is not necessarily found: often the network would settle into near-best solutions that were stable. Through systematic variation of the coupling strengths WI and IV,, a variety of network responses were possible. In particular, by forcing the evolution of the network state to be dominated (at any given time) by only one set of couplings, good network solutions to 100-node problems were found without a large expenditure of computational effort. From
958
D. M. Wells
i
(a)
(b) FIGURE 9. TSP network solutions for nodes randomly distributed in the unit plane: (a) 50-node, (b) 100-node. Typically, valid tours were approximated by two loops through the nodes.
a technical stand-point, there might be certain applications whereby redundancy could be built into a solution for an optimization problem in order to take advantage of the nonlinear properties of the passive oscillators (namely, frequency locking and to some extent phase locking) that help to globally organize these networks. REFERENCES Abbot, L. F. (1990). A network of oscillators. Journal of Physics A (Mathematical and General), 23(16), 3835-3859. Cohen, A. H., Holmes, P. J., & Rand, R. H. (1982). The nature of the coupling between segmental oscillators of the lamprey spinal generator for locomotion: A mathematical model. Journal of Mathematical Biology 13, 345-369.
Fitts, R M. (1954). The information capacity of the human motor system in controlling the amplitude of movement. Journal of Experimental Psycholog): 47, 381-39 I. Fitzhugh, R. ( 1961 ). Impulses and physiological states in theoretical models of nerve membrane. Biophysical Journal. l, 445-466. Grasman, J. ( 1987 ). Asymptotic methods for relo_~ation oscillations and applications. Applied Mathematical Sciences (Vol. 63). New York: Springer-Verlag. Gray, C. M., & Singer, W. (1989). Stimulus-specific neuronal oscillations in orientation columns of cat visual cortex. Proceedings of National Academy of Sciences ( US.A. ), 86, 1698-1702. Hopfield, J. J. (1982). Neural networks and physical systems with emergent collective computational abilities. Proceedings of the National Academy of Sciences ( US.A. ), 79, 2554-2558. Hopfield, J. J. (1984). Neurons with graded response have collective computational properties like those of two-state neurons. Proceedings of the National Academy of Sciences ( US.A. ), 8 I, 30883092. Hopfield, J. J., & Tank, D. W. (1985). "'Neural" computation of decisions in optimization problems. Biological Cybernetics. 54, 141-152. Hoppensteadt, E C. ( 1981 ). Electrical models of neurons. Lectures in Applied Mathematics, 19, 327-344. Kamgar-Parsi, B., & Kamgar-Parsi, B. ([990). On problem solving with Hopfield neural networks. Biological Cybernetics, 62, 415423. Keener, J. R ( 1983 ). Analog circuitry for the van der Pol and FitzhughNagumo Equations. IEEE Transactions on Systems, Man and Cybernetics, 13(5), 1010-1014. Kelso, J. A. S., & Sch6ner, G. ( 1988 ). Self-organization of coordinative movement patterns. Human Movement Science, 7, 27--46. Kopell, N., & Ermentrout, G. B. ( 1988 ). Coupled oscillators and the design of central pattern generators. Mathematical Biosciences, 89, 14-23. Lin, S., & Kernighan, B. W. (1973). An effective heuristic algorithm for the traveling-salesman problem. Operations Research, 21,498516. Matsuoka, K. (1985). Sustained oscillations generated by mutually inhibiting neurons with adaptation. Biological Cybernetics, 52, 367-376. Matsuoka, K. ( 1987 ). Mechanisms of frequency and pattern control in the neural rhythm generators. Biological Cybernetics, 56, 345353. Nagumo, J. S., Arimoto, S., & Yoshizawa, S. (1962). An active pulse transmission line simulating nerve axon. Proceedings of the IRE, 50, 2061-2070. Sch6ner, G., & Kelso, J. A. S. (1988). Dynamic pattern generation in behavioral and neural systems. Science, 239, 1513-1520. Selverston, A. I. ( 1988 ). A consideration of invertebrate central pattern generators as computational data bases. Neural Netuvrks, 1, 109117. Sporns, O., Gaily, J. A., Reeke, G. N., Jr., & Edelman, G. M. (1989). Reentrant signaling among simulated neuronal groups leads to coherency in their oscillatory activity. Proceedings of the National Academy of Sciences ( US.A. ), 86, 7265-7269. Tagliarini, G. A., & Page, E. W. ( 1987 ). Solving constraint satisfaction problems with neural networks. 1EEE First International Conference on Neural Networks--Proceedings. 3, 741-747. Wilson, G. V., & Pawley, G. S. (1988). On the stability ofthe travelling salesman problem algorithm of Hopfield and Tank. Biological Cybernetics, 57, 63-70. Wolfram, S. (1984). Cellular automata as models of complexity. Nature, 311(4), 419-424. NOMENCLATURE
U
slow r e c o v e r y c o m p o n e n t membrane current
of the
Solving Degenerate Optimization Problems V F(V) I
L R C 0
G(V)
W ij
Wi
potential difference across membrane function that relates post-synaptic firing rate to membrane potential external input current inductance associated with slow current events resistance associated with slow current events membrane capacitance membrane potential threshold function that relates fast membrane current to cell potential positive proportionality constants coupling strength from processing unit j to processing unit i coupling strength for a 2-oscillator network coupling strength corresponding to TSP syntactic constraints
959
N d dxr E
E~ SD S T
/x4~ In
Fo FWHM
coupling strength corresponding to TSP cost function number of nodes in the optimization problem number of processing units in network total tour distance distance between node X and node Y computational energy function coupling term of the computational energy function standard deviation of impulse phase measure of system stability oscillator period phase offset index of firing rate impulses corresponding to a given node detection threshold level impulse width of F(V) = 0.5
0893-6080/92 $5.00 + .00 Copyright © 1992 Pergamon Press Ltd.
Printed in the USA. All rights reserved.
ORIGINAL C O N T R I B U T I O N
Solving Degenerate Optimization Problems Using Networks of Neural Oscillators D E R E K M . WELLS Florida Atlantic University
( Ra'eived 1 October 1990: revised and acct7~ted 27 Felwuao' 1992)
Abstract--Optimization problems containing multiple, identical solutions, in particular the travelling salesman problem (TSP), are soh,ed using networks of neural oscillators. Processing units are described by two state variables representing fast attd shin' membrane events, much like the Fitzhugh-Nagumo model neurons, and are passive oseillator~ in that at least two processing units with mutually inhihitoo' couplings are required to generate sustained oscillations. Tile emergent behavior of the network is to produce stable travelling waves (limit ~3'cles) o.[ftring rate impulses in the position coordinate of the network with the phase relationship between active (suprathreshold) processing units determining the TSP solution. This approach has the advantage of avoiding most h~cal energy minima with tile best or near-best sohttian being finmd eveo' time.6~r small net works ( < 10-~processing units). ,,Is network size is enlargened, the time required to reach a stable travelling wave can./igaration increases. To improve tile convergence rate, the strength of the network couplings encoding the cost./imction and tile strength q[the network couplings encoding tile syntactic constraints are altenlatel)' reduced slich that net work evohttion is dictated by only one sel o[eouplings at an)' given time. Good solutions are fimnd./hr net works containing up to 104 processing units without a large ~:vpenditure ~[computational t:[]brt.
Keywords--Neural oscillators. Optimization, TSP. 1. INTRODUCTION
act out-of-phase or at a different frequency) to the dynamics of the cooperative networks would be minimal. There is a large body of literature that deals with modelling coordinated movement in lamprey using nonlinear coupled oscillators (see Cohen, Holmes, & Rand, 1982; Kopell & Ermentrout, 1988) where the coupling strength between oscillators is proportional to a phase difference between body segments. The model neuron in these approaches is based upon the FitzhughNagumo (FN) neuron or variants of the FN neuron (Fitzhugh, 1961; Nagumo, Arimoto, & Yoshizawa, 1962). Abbot (1990) has studied some of the properties of FN networks and has shown that associative memory can be stored as a phase relationship between network elements, different memory patterns having different phase relationships. Perhaps the most fully understood neuronal circuitry that generates rhythmic motor patterns is the stomatograstric system of the lobster (see Selverston, 1988, for a review). The stomatogastric ganglion is composed of two networks which control the muscles of the gastric mill and pyloric regions of the lobster stomach. Selverston has shown that the pyloric network consists of neurons that can be classified as intrinsic oscillators (i.e., neurons that are capable of oscillating individually) whereas most of the neurons in the gastric network can be classified as passive os-
Many of the artificial neural networks that have been developed do not take into consideration temporal factors associated with real biological neural networks. In many biological systems, it is often the underlying neurodynamics that is the most important factor to understand with regard to coordination ( Kelso & Sch6ner, 1988; Sch6ner & Kelso, 1988) and movement in general. Based upon studies of the visual cortex in cats, Gray and Singer ( 1989 ) and Sporns, Gally, Reeke, and Edelman (1989) have suggested that network oscillatory responses (namely, phase and frequency locking) could act as a general mechanism for coordinating spatially distinct populations of cortical cells involved in sensorimotor tasks. The influence of additional inputs (that
Acknowledgements: I would like to thank Scott Kelso, Hirofumi Nagashino, Abhijit Pandya, Tim Kiemel, Gonzales DeGuzmann, and Pier Zanone for commenting on earlier versions of this manuscript and Ed Page for introducing me to the topic of network representations of optimization problems. This research was supported in part by a grant from the Otfice of Naval Research ( N00014-88-J- 1191 ) to JAS Kelso. Requests for reprints should be sent to Derek M. Wells. Center for Complex Systems. Florida Atlantic University. Boca Raton, FL 33431.
949
950
cillators (i.e., neurons that do not oscillate individually but only within a network structure). One approach to modeling passive oscillators has been described by Matsuoka ( 1985, 1987 ). A fast rise-time constant and a slow decay-time constant are associated with each neuron such that a nonmonotonic response curve to a step input develops. Sustained oscillations can be generated when two or more of these model neurons have mutually inhibitory connections. A useful method to study emergent behavior within networks of model neurons is to have the networks solve well defined yet computationally difficult problems. These problems often involve competing sets of constraints such as in optimization problems. The classic biological example of a behavior involving competing sets of constraints is the movement of the arm towards some target. For this motor task, a speed-accuracy tradeoff exists (which was parameterized by Fitts in 1954) wherein accuracy of the final hand position with respect to the target coordinates is compromised at high speeds, Most well defined problems however are distinctly nonbiological, although the network coupling architectures used to solve the given problem as well as the emergent network properties might have some relevance to information processing within real neural subsystems. In this paper, a network of passive oscillators is used to solve an optimization problem containing multiple, identical solutions in particular, the travelling salesman problem (TSP). This problem consists of finding the minimum tour distance through a set of nodes where each node is visited exactly once and where the tour returns to its starting location. For any tour through a set ofn nodes the solution is 2n -fold degenerate in that there are n possible starting positions and two possible tour directions that yield the same tour length. A neural network solution to this problem, which involves encoding a cost function (to minimize tour distance ) and a set of syntactic constraints (to uniquely assign each node a position in the tour) within the coupling matrix of the network, has been pursued by a number of researchers (Hopfield & Tank, 1985; Kamgar-Parsi & Kamgar-Parsi, 1990; Wilson & Pawley, 1988). Each group has reported that TSP scales poorly as the size of the network increases and that the success rate for finding valid tours through the nodes is at best modest. Here, the emergent behavior of the passive oscillator network is to develop travelling waves of firing rate impulses in the position coordinate of the network. The phase relationship between active (suprathreshold) processing unit impulses at a given instant in time dictates a solution to the optimization problem. Individual oscillators within the network exhibit frequency locking (i.e., the variability in individual oscillator frequencies tends toward zero as the system evolves) but do not exhibit phase locking (i.e., the midpoints of the firing rate impulses for each subset of oscillators that represent
D. M. Wells
a valid tour do not synchronize as the system evolves). A measure of the system stability is defined as the variation in this impulse phase offset over time. As network size is increased, the time (number of oscillator periods) required to reach a stable travelling wave configuration increases. To improve the convergence rate for large networks, the strength of the network couplings encoding the cost function and the strength of the network couplings encoding the syntactic constraints are alternately reduced such that network evolution is dictated by only one set of couplings at any given time. Good TSP solutions are found for networks containing up to 104 processing units without a large expenditure of computational effort. 2. T H E M O D E L N E U R O N AND N E T W O R K There are many different models of individual biological neurons which vary in their level of description of biophysical events. A mathematically tractable set of differential equations which are based on the HodgkinHuxley equations was developed independently by Fitzhugh ( 1961 ) and Nagumo et al. (1962). The differential equations describing the model neuron consist of two state variables ( U and V): U is the slow recovery component of the membrane current and V is the transmembrane potential. The Fitzhugh-Nagumo equations have the form: dV
C - ~ = - G ( V ) - U + ~] li,
(1)
i
L -~t dU = - R U + V - 0,
(2)
where C is the membrane capitance, G(V) is a function which mimics the N-shaped relationship between the fast membrane current and the cell potential V and usually is chosen to be a cubic polynomial: G ( V ) = a(I" + pV 2 + qV 3) where a, p, and q are constants that are model specific, R and L are the resistance and inductance associated with the slow current events, T I is the sum of all currents entering the cell both from other neurons in the network and from sources external to the network, and 0 is the membrane potential threshold. Electronically equivalent circuits of networks of these model neurons have been constructed using tunnel diodes (Hoppenstaedt, 1981 ) or sets of operational amplifiers (Keener, 1983). The FN equations have been used to model individual action potentials and at a different level of description (i.e., a slower time scale) to model impulses of neural activity where the output signal is considered to be an average firing rate. These model neurons exhibit rich dynamic properties and are able to mimic a number of the characteristics of real biological neurons including post-inhibitory rebound, plateau potentials, and frequency entrainment (Abbot, 1990).
Solving Degenerate Optimization Problems
951
The differential equations used to represent the model neurons in this paper are similar to the FN equations except that the fast current-voltage curve G ( V ) is some increasing monotonic function instead of being an N-shaped function. In the simulations which follow, G(V) was chosen to be a linear function, G(V) = a V, which greatly reduces the dynamic properties of an individual model neuron: d V / d t = d U / d t = 0 for large t. The neuron output response is modelled as an averaged firing rate which is related to transmembrane potential by a function F(V) that is usually chosen to be sigmoidal in shape. The equations for the state variables of the model neurons ( U i, Vi, i = 1. . . . . N) have the form: dVi dt
= - a I ~ - U, + I, + dUi dt
=
F(V, )
-{3(u,
-
y = 0.01
y
=
0,05
y=0.1
Z kL
y =0.2
y =0.5
N ~ woF(Vj), j=Lj,,i
(3)
time
v v , ),
(4)
FIGURE 1. Transient response of an isolated model neuron to a step input for a range of 3, values (with a/B = 10).
1 -
l +
y = 0.001
e -x(v'-°')
"
5)
where wo is the coupling strength from neuron j to neuron i, [3 = R L -t , 3' = R - ~and ~ is a proportionality constant which specifies the steepness of F(V); C was set equal to one which does not reduce the generality of the equations. The influence of the parameter 3' (with a/[3 = I0) on the response characteristics of a single model neuron to a step input ( I = 0 to I =/,,) is shown in Figure 1. In the simulations which follow, 3' was chosen such that the transient response of a model neuron is overdamped: 3' < ( a - [3)2/4[3. Matsuoka (1985, 1987) has used a similar set of equations as a model for the phenomenon of accommodation (the decrease in firing rate over time to a constant input) and as an approach to generating sustained oscillations in model networks. In the limit as 3' tends to zero, U = Uoe -~' also tends toward zero for large t hence, eqns ( 3 ) to ( 5 ) reduce to the Hopfield differential equations (Hopfield, 1982, 1984). Results presented in this paper were also reproduced using a piece-wise linear version of F(V) with no change in the qualitative features of the network behavior. The magnitude of the external input current, I~, and the strength and architecture of the processing unit couplings, ~'0, determines the network response.
F(I') =
To understand why sustained oscillations can be generated when two or more of these model neurons are coupled together, consider a network of two neurons with mutually inhibitory connections ( w~2 = w2~ = - 77) and constant external input current (I, = Iz = Io) as shown in Figure 2a. To facilitate analysis let F(V) be piecewise linear with slope 1/6 (0 = 0):
0< z'
1,
I ' > t3
O,
I'<0
(6)
and let a ,> /3 so that the system acts as a relaxation oscillator ( d V / d t ,> d U / d t ) . In this case, the time spent in the firing rate transition region (0 < F ( V ) < 1 ) is far less than the time spent in the saturation regions ( F ( V ) = O, F(V) = I) which permits the system dynamics to be approximated by only two of the four equations of motion at any given time. That is, if F( l~ ) equals 0 or 1 then Ui governs the overall system dynamics ( f, changes so rapidly with respect to Ui in this region that it can be considered constant for any value of Ui ). From eqn (3), d l ~ / d t = 0 gives, 1
I', = - ( - U ,
+ I , , - ~F(I,))),
(7)
and substituting eqn (7) into eqn (4) yields, dt
~ I +
U, + ~ ~7 (I,, - ~F(1))). -
(8)
Whereas. if0 < F(I~ ) < 1 then I~ governs the overall system dynamics ( Ui = U,,. = constant): dI',
2.1. Two Coupled Oscillators
I" y,
dt
-
c~I,} -
U,c +
1,,-
rlF(I~).
(9)
The usual approach for determining if a given system is capable of sustained oscillations is to calculate where the fixed points for the system ( d U / d t = d V / d t = O) lie with regard to the current state of the system; if stable fixed points exist in the same firing rate region as the current system state then sustained oscillations are not possible (see Grasman, 1987, for an overview
952
D. M. Wells I.
Io (w~: = w n < 0 )
(a)
u~ ',',
.:..÷.:,~.,-,.'. ,.:.
~,,.,,.:, ~.,.:.L.,,,.,. '-?,; ::.,
a
i ' _i . oi _~ o _
I.o-,~.,°-~,
Ui
(b) FIGURE 2. A two-oscillator network: (a) coupling architecture, (b) state-space trajectory (as denoted by arrows) projected onto the (U1, U2)-plane. There are nine regions corresponding to different values of F(V,) and F(V2): F(VI ) = O, F(V~ ) = 1, 0 < F(VI ) < 1 (i = 1, 2) where within each region limits for the system state variables were determined (e.g., in region I where F(V~) = 1 and F(V2) = 0, V 1 > ~, V 2 < 0, U1 < Io - / ; a and U2 > Io- ~). For the purpose of clarity, only regions 1, 3, and 8 (projected onto the (U~, U2)-plane) are shaded. Different fixed points for the system (dU~/dt = dVddt = 0) are associated with each region (e.g., the fixed point (projected onto the (U~, U2)-plane) for the system in region I is denoted by "b," whereas the fixed point associated with region 5 is denoted by "a"). Sustained oscillations were possible because the fixed points for the system in any region were either outside of the region the system state was in or were unstable.
of relaxation oscillators). Note that different fixed points are associated with different firing rate regions. For this two-neuron network, the firing rate space (F(Vt), F(V2)) can be partitioned into nine distinct regions ( F( V/ ) = 0, F ( V i ) = 1 , 0 < F ( V ~ ) < l , i = I, 2) where within each region the limits for each of the state variables can be established using eqns ( 6 ) to (9). The limits for the slowly changing variables Ut and U2 in each of the nine firing rate regions are shown in Figure 2b. For example, the parallelogram in the center of the (Uj, U2)-plane indicates the region where both F( Vt ) and F(1/"_,)are in transition. Under the conditions that ~ > 6a and 77 < Io < 6(a + ~), it can be shown that the fixed points associated with any firing rate region of the system either lie outside of the region the system state is in or are unstable. For this system, a
zone of bistability exists within the square that is centered in the parallelogram; the fixed points associated with regions 1 and 5 are stable but lie outside of this zone, whereas the fixed point associated with region 9 is located in this zone but is unstable (a saddle node point). A trajectory of the system state projected onto the ( U~, U2)-plane is denoted by the arrows in Figure 2b, in which sustained movement occurs between firing rate regions 1 and 5. 2.2. Network Representation of T S P In order to solve an optimization problem such as TSP using a neural network, an appropriate representation of the problem must be formulated. This involves choosing a suitable network topology as well as defining external input currents ( l i ) and coupling strengths between processing units (w,i). One such network mapping scheme for TSP has been described by Hopfield and Tank (1985) in which they represented a valid tour through a set o f n nodes by an n × n-permutation matrix; an element (i, j ) in the matrix represents the probability that node i will be in position j of the tour. An example of a valid tour through a set of five nodes is illustrated in Table I. The total tour distance can be decoded from this representation by taking the following sum: d = d.4c + dcB + dBE + dEo + doA, where the order of the n nodes in the tour are ACBEDA. A valid tour requires that exactly one element in each row and column of the matrix be set to "1," all other elements equaling zero hence, a total of n! valid tours are possible. However, only (n - l ) ! / 2 of the tours are distinct in that for any valid tour through the nodes there are n possible starting positions and two possible directions which yield the same tour length. Based on this representation scheme for valid TSP solutions, Hopfield and Tank (1985) developed a neural network in which each processing unit in the network corresponded to an element in the permutation matrix. For 3' = 0, Hopfield (1984) has shown that the equations of motion for networks with certain coupling architectures ( wij = ~t~,, Wig= 0) lead to stable fixed point configurations o f F ( V, ). A computational energy function that minimizes total tour length while simultaneously satisfying the row and column constraints can TABLE 1 The Permutation Matrix Representing the 5-node Tour ACBEDA
Position Node
1
2
3
4
5
A
I
0
B
0
0
0
0
0
I
0
C D E
0 0 0
I 0 0
0
0 0 0
0 0 I
0 I 0
Solving Degenerate Optimization Problems
953
be associated with the network. This function takes the general form: 1 N
E
= -
~
~
ill
Y÷X
dzF-J(z),
(10)
i=1
and is related to the equation of motion for a processing unit (in the special case where 3, = 0) by dV~/dt =
-OE/OF( V~ ). In order to more clearly show how the constraints and cost function for TSP are mapped into the network couplings, matrix notation will be used to specify the firing rate of an individual processing unit with the indices X and Y referring to nodes (rows of the network) and the indices i a n d j referring to positions in the tour (columns of the network); for example, F(Vxi) is the firing rate for the processing unit in row X and column i of the network. The coupling term, E,., of an energy function that is minimal when the row and column constraints of TSP are satisfied may now be specified: 1,I,"t
Z Z ~ F(Vx,)F(Vxj) X
i
j÷i
+-g~ ~ ~ ~ F(I'v,)F(I,'r,), --
i
(11)
X Y ÷X
where W~ is the coupling strength corresponding to the row and column constraints and all summations are from 1 to n (the number of nodes). Often, another term proportional to
#
E,. = T ~x ~i Z F(Vxi)F(I/xj) )÷i
W2
+ ~ E Z Z F(VxJF(VvJ +--~ Z Z Z 4rvF(Vx,) --
i
X
}'÷X
X"
Y+X
X(F(Vya+,)F(I/y,,_,)), +
(13)
where W2 is the coupling strength corresponding to the cost function, dvv is the distance between nodes X and Y and subscripts for the position coordinate of the network are expressed modulo n (i.e., F(Vr.,+j) = F(Vr.j). By taking the derivative - 0 E,./OF(Vxi), the coupling term in the equation of motion for each processing unit in the TSP network is determined:
(14)
Y÷X
The local neighborhood for a processing unit in a network representing a 5-node TSP is illustrated in Figure 3. Network solutions to TSP are possible when appropriate choices for the coupling strengths, Wt and W2, are made. The differential equations governing the network are summarized below: dl':v' -
dt
a V x , - Ux, + Io- l'Vt( ~ F(Vxj) + ~ F(Vr,) )
j,,,
r,,x
-dxvW2(~F(Vr.,+,)+ ~ F ( V r . H ) ) , Y÷X
(15)
Y÷X
dUxi _ dt
B( Ux; - 7 l';vi),
(16)
1
F(Vx,) - I + e -xv-v' "
(17)
In the case where 3' > 0, the processing units as defined by eqns ( 15 ) to ( 17 ) do not settle into a stable fixed point configuration in the F(V) space; instead they form limit cycles in the position coordinate of the network. Stability, in the context of this network of oscillators, is defined as the formation of travelling Position 1
2
3
,,2)
is included in E,. that specifies the total neural activity of the network. However, it can be shown (Tagliarini & Page, 1987) that this global coupling term is not required if the row and column coupling strengths (Wt) and the external input current (Ix~ = I,,) are suitably chosen. A more complete coupling term of an energy function that also favors short tours has the following form:
14 "1
~ F(Vri)) Y÷X
-dx,.l.t/'2(~_~F(Vr.i+,)+ ~F(Vy.i_l) ) .
woF(V~)F(Vj) - Z F(V, )I, +a
=T
j÷i
N
Z t=l= j = t
E,.
-H.'I(~F(I.~)+
o
4
5
@®OO ®0
00
FIGURE 3. The local coupling neighborhood of processing unit B2. Processing units that are symbolized by circles containing squares are negatively coupled to unit B2. Circles with open squares encode the row and column constraints whereas circles with filled squares encode the cost function which minimizes total tour distance. Note that columns 1 and 5 are considered adjacent when applying the cost function so that no edge effects occur. This particular coupling architecture allows travelling waves of firing rate impulses to form in the position coordinate of the network.
954
D. M. Wells
waves (limit cycles) of firing rate impulses that maintain a constant phase relationship. The set of processing units that are simultaneously active (F(V) > threshold) form a solution to TSP. Returning briefly to the 5-node tour illustrated in Table 1, a stable travelling wave of firing rate impulses results in the order of nodes for valid solutions changing from ACBEDA to CBEDAC to BEDACB, and so forth. For a given simulation, a travelling wave in the opposite direction (ACBEDA to DACBED to EDACBE, and so forth) is equally probable although, the direction of the travelling waves can be prespecified if the coupling term corresponding to the cost function in eqn (15) is modified: -dxrW2((l +~) ~ F(V,.,+,)+ ~ F(Vv,,-,)), (18) Y~X
Y~'X
where I~l ~ 1. It should be noted, however, that ~ does not eliminate the 2-fold degeneracy in tour direction, it simply improves the convergence rate to stable limit cycles. For the simulations which follow ~ was set equal to zero. A time period, T, can be associated with the travelling wave and is defined as the number of system iterations required to change the order of suprathreshold processing units from ACBEDA to CBEDAC to BEDACB, and so forth. In the following section T is used to scale the evolution of the network state. 3. S I M U L A T I O N RESULTS Simulations were performed for networks ranging in size from 102 processing units to 10 4 processing units. For small networks, the set of parameter values used in the simulations was: a -- 0.1,/3 = 0.01, 3, = 0.1, X = 0.035, Wt = 12.0, W2 = 8.0, and Io = 20.0. Distances between individual nodes (dxr) were scaled such that an average separation distance of 0.5 was achieved. For large networks, the ratio W~:W2 was systematically varied over the time course of the simulation (more details about this schedule are described in Section 3.2). The magnitude of the external input current, I,,, determines the number of system iterations a processing unit remains suprathreshold. That is, Io dictates the width of the firing rate impulse. I,, was set such that the impulse width associated with each processing unit at one-half the maximum firing rate (F(V) = 0.5 ), the so-called full width at half maximum value ( F W H M ) , was 1.5 times the magnitude of T. This value for I,, was required in order to generate stable travelling waves through the network and was determined empirically. The corresponding size of T was determined experimentally using a simple peak finding routine which measured the midpoint of the time a given processing unit was suprathreshold. Valid solutions to TSP were determined by detecting which processing units in the network were suprathreshold (F(V) > Fo) at time t, where Fo is a detection threshold level that is limited to the range of F(V). To ensure that solutions to TSP
are possible, the value for Fo must be chosen such that the width of the firing rate impulse at the corresponding level ofF,, is less than 2T. For Io equal to 20.0, a suitable choice for F,, is 0.75. If exactly one unit was active for each row and column of the network matrix then a valid solution to TSP was recorded. Equations ( 15 ) and ( 16 ) were solved using the Euler method for numerical integration in which each processing unit of the network was updated once (in a random sequence) per system iteration. The initial state of the system was determined as follows: F(Vxi ) = F W H M / n T = 1.5/n (this value was based, approximately, upon the amount of neural activity that is present at any given time after the system has settled into a travelling wave configuration), Vxi was calculated using eqn (5) and Uxi was obtained from the following relation:
Uxi = -c~[~r, + I,, - 2(n - I)(W I + dxrB~)F(VxA. (19) 3.1. Small Networks Emergent behavior within complex systems has been classified by Wolfram (1984) into four general types: (a) fixed point, (b) periodic, (c) chaotic, and (d) intermittent. For a network consisting of the TSP coupling architecture (in the appropriate parameter regime), the emergent behavior of the oscillators was to produce travelling waves with constant period T, type 2 behavior (Figure 4a). Initially, subsets of travelling waves formed in competing directions (waves moving upward versus downward) however, after an elapsed time of about n T, only travelling waves moving downward persisted. Global organization of the network oscillators is not possible for the case where 3' equals zero; fixed point solutions to TSP (which are highly dependent upon initial conditions) are formed at elapsed times of less than n T. This property of the network to produce stable travelling waves was dependent upon the values of W~ and W2. If, for instance, Wz was set equal to zero, oscillatory responses still formed but appeared to be aperiodic (Figure 4b). In addition, intermittent type behavior was possible for intermediary values of Wz. Throughout the rest of this paper, only TSP networks that generate stable travelling waves will be considered. Figure 5 shows typical solutions to 10-node and 20node problems. The best or near-best solution was found with 100% success rate. For some trials, a number of other valid solutions were found (usually with larger tour distances) before the network settled into a stable travelling wave configuration. Coupled oscillators are often characterized by phase and frequency locking. For the TSP network, oscillators representing a given node tended toward a common frequency ( T ) -~ as the system evolved ( Figure 6a) but the phase offsets between oscillator impulses for processing units that were simultaneously active did not
Solving Degenerate Optimization Problems
-~> . 2 ,
:
~,
955
.-),_2. •
~ - " - ~
",
~
~
~
(a)
"-~'-~---)..~
100: Time
(a) I
.
~
~ .
~
L -------r~
~ c __
~ - ~ s
--
~ ~._LA
A
Co) FIGURE 5. TSP network solutions: (a) 10-node, (b) 20-node.
, ~--~.~_
100
~
"
tend toward zero. Figure 6b shows the standard deviation in impulse phase, SD, for 6 trials of the same 10node problem where only the processing unit update schedules varied. SD was calculated as follows:
~
Time
SD (b) FIGURE 4. Time series for each processing unit firing rate, F(V), in a 100 unit TSP network where outputs that correspond to the 10 possible positions in the tour (columns in the network matrix) are displayed consecutively: (a) W2 = 8.0. Sustained firing rate impulses are generated and form travelling waves of constant period, T, in the position coordinate of the network; (b) W2 = 0.0. Sustained firing rate impulses are generated but appear to be aperiodic.
\/Z7-2 (~,,., - ~b,..,)2 n-I V
(20)
where 4~.,. is the midpoint of the m-th firing rate impulse for the active processing unit representing the ith node and (4~,,. - ¢~,,.) is the phase offset between the active processing units representing nodes i and 1. Midpoints corresponding to the m-th firing rate impulse were located within the range: ~b].m -- 0.5T < ~bi.,. < ~b~,,. + 0.5T. If no midpoint was found within this win-
956
D. M. H"ells A measure of system stability (S) was formed by calculating the change in phase offset as the system evolved: / ~ ['=2 (A~bi''' - A4~''-t)2 , n- 1
S = O4 1D 0 0-
0
0
t
I.~ 1
Time / 5 nT (a)
(21)
where A$~.,, = $~.,, - St ..... A plot of S for networks of various sizes is shown in Figure 7 (where/,, was constant: F W H M / T = 1.5). As network size was increased, the time to reach a stable travelling wave configuration increased and at this level of external input current, stability was not obtained for the network representing a 40-node problem. A stable solution might be reached by an elapsed time of 5n T by increasing the external input current but at some level of input, the width of the firing rate impulses (at the level ofF,,) will become greater than 2 T a n d solutions to TSP would no longer be detectable using the approach described in this section. 3.2. Large N e t w o r k s
I-O3
0
1 Time / 5 nT Co)
FIGURE 6. Emergent TSP network properties: (a) the time pedud (number of system iterations) between fidng rate impulses for processing units representing each node (row) of the network tended toward a common period, T; (b) the standard deviation of the impulse phase between active processing units did not tend toward zero as the system evolved. Over multiple trials (six of which are shown), where the resultant tour through the nodes was identical, the standard deviation in the impulse phase ranged from about 0.1T to 0.4T.
dow (which might occur prior to the formation of a stable travelling wave) then 4~,,,, was set equal to 4h .... + 0.5T. For each trial in Figure 6b, the resultant tour through the nodes was the same whereas the range in standard deviation for the impulse phase varied between 0 . I T and OAT, indicating that multiple phase configurations (at a time scale less than T) were stable. From a practical standpoint, the ability to detect valid tours becomes degraded for network oscillators that have nonzero phase offsets.
In order to improve the convergence rate, as well as to solve, large scale TSPs ( up to 100 node problems), the coupling strength ratio 14'j :I,V_, was alternated throughout the simulation such that either the couplings that encoded the row and column constraints or the couplings that encoded the cost function dominated the evolution of the network state, but never both sets of couplings simultaneously, as was the case for small networks where the ratio remained constant (H,'t :H~ = 3:2). Initially, the coupling strength for the row and column constraints ( W~ ) was set to a value of 4.0 ( H"t : I4 2 = 1:2). This lower value for H.'t caused the number of processing units that were simultaneously active in each row and column of the network matrix to increase 1
n=:O . . . . . .
n=20
n=30 n=4O
i i
:,.
,
'i
' .
0
, I
Time 15 nT
FIGURE 7. A plot of the stability measure, S, for networks ranging in s i z e from 100 units to 1,600 units (FWHM/T = 1.5). It should be noted that the time required to reduce S to about zero varied widely for different trials of the s a m e problem.
Soh, ing Degenerate Optimization Problems from 1 to about 5. L Active processing units tended to cluster into either one or two regions in the position coordinate of the network matrix hence, an approximate solution to TSP was formed. However, as the results in this section will show, there was no guarantee that the best solution to TSP was one of the permutations within the approximate solution. The time to reach a stable travelling wave configuration was similar to the time required for small networks: between n T and 2n T. A number of possible methods to solve the TSP given this stable travelling wave configuration can now be implemented. One approach to finding good TSP solutions would simply be to perform a sequential search through all possible combinations of neurons that comprised the approximate solution. If the number of neurons that were active in a given row or column was far less than n (the number of nodes) then the total number of possible tours would be greatly reduced from n! The approach taken here was to increase the coupling strength ratio H'~:HJ2 to 4:1 (W~ = 16.0, W_, = 4.0) which enabled the couplings that encoded the row and column constraints to subsequently dominate the system evolution instead of the couplings that encoded the cost function. In order to favor a final TSP solution that remained within the bounds of the initial approximate solution, the magnitude of the WTcouplings between processing units that were not both simultaneously active in the approximate solution was increased to a value of 8.0. Solutions were found with 100% success rate although the travelling waves became unstable if the system state was allowed to continue evolving given the high coupling strength ratio. By periodically alternating the coupling strength ratio (a time schedule of which is illustrated in Figure 8), multiple solutions to a given TSP were found. Results for the Lin-Kernighan tour (Lin & Kernighan, 1973) were consistently about 4.8 to 5.0 in tour length compared to the best value of 4.32.-' Network solutions to problems involving 50 and 100 nodes (where nodes were randomly distributed in the unit plane) were also readily found: Figure 9 shows typical results. In general, it was found that processing units which comprised the approximate TSP solution clustered into two regions in the position coordinate of the network (i.e., two travelling waves were simultaneously formed) hence, the approximate solution (and consequently, the resultant valid solution) involved two redundant loops through the nodes. The fact that the number of stable travelling waves that simultaneously formed within a network did not conSimultaneously active processingunits in each row and column of the matrix can also be achieved by increasing the external input current, 1,,, or decreasing ,X. 2Approximatecoordinates for the kin-Kernighan tour were provided by Ed Pesulima.
957 Wt : W, 4
3 2 l
0
I 2nT
I ~, t > T
Time
FIGURE 8. Time schedule for alternating the coupling strength ratio Wl :W=. For W1:14/2 = 0.5, stable traveling waves formed after an elapsed time of between n T and 2n T with multiple ( 4 6) active processing units in each row and column of the network matrix; an approximate TSP solution was thereby established. At time 2n T, W~ :14/2was increased to 4.0 so that the row and column couplings dominated the system evolution. In order to favor a valid TSP solution that was one of the permutations of the approximate solution, the coupling strength ratio between processing units that were not both a part of the approximate solution was decreased to 2.0. These higher coupling strength ratios were maintained until a valid solution was detected (as indicated by arrows), after which the ratio was reduced to its original value (0.5) and the process was then repeated. Multiple solutions to a given TSP were thereby found.
tinue to increase beyond two (as network size increased) might reflect some inherent limitation in this network approach to solving TSP. One such limitation is the fact that a stable traveling wave configuration of processing unit activity levels does not eliminate the 2fold degeneracy in tour direction. Clearly, more studies are required in order to better characterize the global organizational properties of these large passive oscillator networks. 4. S U M M A R Y Throughout this paper, some of the properties of networks of passive oscillators have been studied in the context of solving a well-known optimization problem that has multiple solutions of identical cost namely, the travelling salesman problem (TSP). The coupling architecture for the local neighborhood of each processing unit in the oscillator network allowed travelling waves of firing rate impulses to be formed along the position coordinate of the network matrix. For an n node problem in which stable travelling waves formed, a total of n solutions of identical tour length were detected for each pass of the travelling wave through the network. It is believed that by forcing the network to find multiple solutions to TSP better and more consistent results were realized. However, the best solution is not necessarily found: often the network would settle into near-best solutions that were stable. Through systematic variation of the coupling strengths WI and IV,, a variety of network responses were possible. In particular, by forcing the evolution of the network state to be dominated (at any given time) by only one set of couplings, good network solutions to 100-node problems were found without a large expenditure of computational effort. From
958
D. M. Wells
i
(a)
(b) FIGURE 9. TSP network solutions for nodes randomly distributed in the unit plane: (a) 50-node, (b) 100-node. Typically, valid tours were approximated by two loops through the nodes.
a technical stand-point, there might be certain applications whereby redundancy could be built into a solution for an optimization problem in order to take advantage of the nonlinear properties of the passive oscillators (namely, frequency locking and to some extent phase locking) that help to globally organize these networks. REFERENCES Abbot, L. F. (1990). A network of oscillators. Journal of Physics A (Mathematical and General), 23(16), 3835-3859. Cohen, A. H., Holmes, P. J., & Rand, R. H. (1982). The nature of the coupling between segmental oscillators of the lamprey spinal generator for locomotion: A mathematical model. Journal of Mathematical Biology 13, 345-369.
Fitts, R M. (1954). The information capacity of the human motor system in controlling the amplitude of movement. Journal of Experimental Psycholog): 47, 381-39 I. Fitzhugh, R. ( 1961 ). Impulses and physiological states in theoretical models of nerve membrane. Biophysical Journal. l, 445-466. Grasman, J. ( 1987 ). Asymptotic methods for relo_~ation oscillations and applications. Applied Mathematical Sciences (Vol. 63). New York: Springer-Verlag. Gray, C. M., & Singer, W. (1989). Stimulus-specific neuronal oscillations in orientation columns of cat visual cortex. Proceedings of National Academy of Sciences ( US.A. ), 86, 1698-1702. Hopfield, J. J. (1982). Neural networks and physical systems with emergent collective computational abilities. Proceedings of the National Academy of Sciences ( US.A. ), 79, 2554-2558. Hopfield, J. J. (1984). Neurons with graded response have collective computational properties like those of two-state neurons. Proceedings of the National Academy of Sciences ( US.A. ), 8 I, 30883092. Hopfield, J. J., & Tank, D. W. (1985). "'Neural" computation of decisions in optimization problems. Biological Cybernetics. 54, 141-152. Hoppensteadt, E C. ( 1981 ). Electrical models of neurons. Lectures in Applied Mathematics, 19, 327-344. Kamgar-Parsi, B., & Kamgar-Parsi, B. ([990). On problem solving with Hopfield neural networks. Biological Cybernetics, 62, 415423. Keener, J. R ( 1983 ). Analog circuitry for the van der Pol and FitzhughNagumo Equations. IEEE Transactions on Systems, Man and Cybernetics, 13(5), 1010-1014. Kelso, J. A. S., & Sch6ner, G. ( 1988 ). Self-organization of coordinative movement patterns. Human Movement Science, 7, 27--46. Kopell, N., & Ermentrout, G. B. ( 1988 ). Coupled oscillators and the design of central pattern generators. Mathematical Biosciences, 89, 14-23. Lin, S., & Kernighan, B. W. (1973). An effective heuristic algorithm for the traveling-salesman problem. Operations Research, 21,498516. Matsuoka, K. (1985). Sustained oscillations generated by mutually inhibiting neurons with adaptation. Biological Cybernetics, 52, 367-376. Matsuoka, K. ( 1987 ). Mechanisms of frequency and pattern control in the neural rhythm generators. Biological Cybernetics, 56, 345353. Nagumo, J. S., Arimoto, S., & Yoshizawa, S. (1962). An active pulse transmission line simulating nerve axon. Proceedings of the IRE, 50, 2061-2070. Sch6ner, G., & Kelso, J. A. S. (1988). Dynamic pattern generation in behavioral and neural systems. Science, 239, 1513-1520. Selverston, A. I. ( 1988 ). A consideration of invertebrate central pattern generators as computational data bases. Neural Netuvrks, 1, 109117. Sporns, O., Gaily, J. A., Reeke, G. N., Jr., & Edelman, G. M. (1989). Reentrant signaling among simulated neuronal groups leads to coherency in their oscillatory activity. Proceedings of the National Academy of Sciences ( US.A. ), 86, 7265-7269. Tagliarini, G. A., & Page, E. W. ( 1987 ). Solving constraint satisfaction problems with neural networks. 1EEE First International Conference on Neural Networks--Proceedings. 3, 741-747. Wilson, G. V., & Pawley, G. S. (1988). On the stability ofthe travelling salesman problem algorithm of Hopfield and Tank. Biological Cybernetics, 57, 63-70. Wolfram, S. (1984). Cellular automata as models of complexity. Nature, 311(4), 419-424. NOMENCLATURE
U
slow r e c o v e r y c o m p o n e n t membrane current
of the
Solving Degenerate Optimization Problems V F(V) I
L R C 0
G(V)
W ij
Wi
potential difference across membrane function that relates post-synaptic firing rate to membrane potential external input current inductance associated with slow current events resistance associated with slow current events membrane capacitance membrane potential threshold function that relates fast membrane current to cell potential positive proportionality constants coupling strength from processing unit j to processing unit i coupling strength for a 2-oscillator network coupling strength corresponding to TSP syntactic constraints
959
N d dxr E
E~ SD S T
/x4~ In
Fo FWHM
coupling strength corresponding to TSP cost function number of nodes in the optimization problem number of processing units in network total tour distance distance between node X and node Y computational energy function coupling term of the computational energy function standard deviation of impulse phase measure of system stability oscillator period phase offset index of firing rate impulses corresponding to a given node detection threshold level impulse width of F(V) = 0.5