1. Introduction
Combinatorial optimization is a branch of mathematical optimization which has wide applications in artificial intelligence, theoretical computer science, applied mathematics, and so on [
1]. One of the most famous combinatorial optimization problems is the traveling salesman problem (TSP), which is ubiquitous in real-life scenarios such as supply chain management, aircraft routing, and job-shop scheduling. TSP aims to find the shortest route that visits every city exactly once and returns to the starting point, which is, however, an NP-hard problem that cannot be solved accurately. Swarm intelligence optimization algorithm [
2,
3] and artificial neural network [
4,
5] are the two mainstream methods to deal with TSP. Based on chaotic simulated annealing, Chen et al. proposed a transient chaotic neural network (TCNN) model to solve TSP in 1995 [
6]. Compared with the Hopfield neural network (HNN) and chaotic neural network, TCNN shows stronger search ability. However, as the annealing processing of TCNN is determined by a damping factor, TCNN is easy to fall into local optimum when solving large-size TSP. In 2004, by adding decaying stochastic noise into TCNN, a noisy chaotic neural network (NCNN) was constructed [
7]. The simulation results show that the NCNN achieved higher global optimal rate than TCNN when solving large-size TSP. Additionally, in order to provide the ability to characterize local features, Gauss wavelet self-feedback was introduced into TCNN [
8], resulting in the reduction in local minima rate. However, these works all aimed to improve the global optimum rate by modifying the mathematical expression of annealing processing, while paying little attention to the physical implementation of annealing function or the whole algorithm.
Recently, memristor has attracted widespread interest on combinatorial optimization. Pershin et al. [
9] showed that the memristive networks can efficiently solve shortest-path optimization problems thanks to its electrically programmable conductance, which shows potential as the hardware basis for implementing HNN or TCNN. For example, Fahimi et al. [
10] proposed a weight annealing approach for HNN on a fabricated 20 × 20 analog-grade TiO
memristive crossbar with the help of the designed mixed-signal accelerators. Liu et al. [
11] proposed to use the conductance values of the Hewlett-Packard (HP) memristor, which can be adjusted dynamically according to the annealing algorithm, as the self-feedback connection weights in TCNN. Based on this, an HP memristor-based transient chaotic neural network (HPM-TCNN) was developed. Simulation results demonstrate that HPM-TCNN could quickly converge toward the global minimum, which was attributed to the efficient nonlinear annealing process stemming from the strong nonlinearity of the HP memristor. Moreover, Yang et al. [
12] experimentally demonstrated the implementation of TCNN on a memristor crossbar, which realized high probability of global optimum and fast convergence by using the nonlinear long-term depression (LTD) processes of the diagonal memristors in the crossbar for annealing.
It is noteworthy that the memristors used so far for the implementation of TCNN are all filamentary memristors, which operate through the formation and rupture of conductive filaments. Although the filament formation/rupture is a nonlinear ion migration process, which can be well used to implement the annealing function, its intrinsic stochasticity may cause imprecise control of conductance, relatively poor endurance, and large device-to-device variation in filamentary memristors. This in turn limits the performance of a practical memristor-based TCNN. To address this issue, memristors with more controllable nonlinear switching dynamics are highly desirable.
A ferroelectric memristor emerges as a suitable candidate because its conductance change relies on the switching of lattice non-centrosymmetry-induced polarization, which is a purely electronic mechanism with strong nonlinearity and high controllability [
13]. On one hand, the nonlinearity of polarization switching is fundamentally different from that of filament formation/rupture, offering an opportunity for the ferroelectric memristor to realize a more favorable annealing function. On the other hand, the high controllability of polarization switching enables the ferroelectric memristor to exhibit precise conductance control [
14], high endurance [
15], and small device-to-device variation [
16]. These features make the ferroelectric memristor a promising building block for the hardware TCNN. However, it should be pointed out that although ferroelectric memristors have been widely investigated as synaptic devices for artificial neural networks [
15,
17] and convolutional neural networks [
18], their use in the hardware implementation of TCNN has never been attempted.
In this study, we theoretically establish a ferroelectric memristor-based TCNN (FM-TCNN) model and demonstrate its high performance in solving combinatorial optimization problems such as the TSP. We first establish a single-neuron model with an annealing function implemented by a prototype ferroelectric memristor, i.e., ferroelectric tunnel junction (FTJ). The conductance of the ferroelectric memristor is used as the self-feedback connection weight. In addition, the conductance is adjusted by applying a voltage pulse whose amplitude depends linearly on the neuron output, resulting in the dynamical modulation of the self-feedback connection weight. It is shown that the neuron can successfully converge to the stable state from the chaotic state, and the convergence speed is modulated by the slope between the pulse amplitude and the neuron output. Then, the FM-TCNN model is applied to solve the TSP. It achieves shorter average path distance, higher convergence speed, and higher global optimal compared with TCNN and HPM-TCNN, demonstrating its effectiveness in solving combinatorial optimization problems.
The remainder of this paper is organized as follows. In
Section 2, the ferroelectric memristor is modeled. The FM-TCNN model is described in
Section 3. In
Section 4, the application of FM-TCNN in solving TSP is carried out, and the experimental results are also discussed. The final section concludes this study.
2. Ferroelectric Memristor Model
A prototype ferroelectric memristor, i.e., FTJ, is used in this study. FTJ consists of an ultrathin ferroelectric film (i.e., the tunnel barrier) sandwiched between two different metal electrodes. The tunnel conductance depends on the average height of the tunnel barrier, which is modulated by the relative fraction of ferroelectric domains with polarization oriented toward one (or the other) electrode. More specifically, when all the domains are aligned downward, the barrier height at the ferroelectric/bottom (top) electrode interface is reduced (enhanced). The reverse scenario occurs when all the domains are aligned upward. Assuming that the bottom electrode has a larger screening length over dielectric constant than the top electrode, the polarization charge-induced barrier height modulation at the bottom interface is therefore more significant than that at the top interface. As a result, the FTJ exhibits the lowest (highest) average barrier height when all the domains aligned downward (upward), leading to the ON-state conductance (OFF-state conductance ).
As the domains are gradually switched from downward to upward (from upward to downward) through the application of voltage pulses, the device conductance (
G) gradually evolves from
to
(from
to
). To realize the gradual switching of domains, the number of domains underneath the electrode area is assumed to be sufficiently large. This can be achieved in two ways: (1) using sufficiently large electrodes and (2) using the domain engineering to reduce the domain size. For practical application, the second way is more favorable, because a small electrode area is always required for high integration density. Previous studies have demonstrated that the domain size can be reduced to be below 5 nm [
19]. For a 130 nm technology node, the number of domains in an individual device can reach >670, which may be sufficiently large for ensuring the gradual switching of domains.
Using a simple model of parallel resistances, the relationship between
G and the area fraction of downward domains
S can be described below [
20]:
where
and
have fixed values once a specific FTJ is used (here,
and
are used, which are extracted from the experimental data of a Pt/BaTiO
3/Nb:SrTiO
3 FTJ [
21]), while
S is a value between 0 and 1 depending on the amplitude and duration of the applied pulse.
To establish the relationship between
S and the amplitude and duration of the applied pulse, the polarization switching kinetics need to be understood. It is widely accepted that the polarization switching in FTJ can be described by a nucleation-limited-switching (NLS) model [
22]. In this model, the ferroelectric film is considered to consist of many different basic areas with independent switching kinetics. The switching time (
) in each basic area is governed by the domain nucleation, while the time for the domain growth can be neglected. Therefore, a distribution of
exists, and a Lorentzian distribution function is assumed here [
23]:
where
is the mean switching time, and
is the half-width at half-maximum. The dependence of
on the applied electric field
E (
, where
V is the voltage of the applied pulse,
d is the ferroelectric film thickness, and symbol ‘
’ means to take the absolute value.) follows the characteristic Merz’s law:
while the relationship between
and
E is given by
Simulations were performed using the experimental data given in ref. [
21], where
d = 2.8 nm and
V varies in a certain range. Some examples of the simulated Lorentzian distribution of
at different
V are shown in
Figure 1. It is seen that when
V becomes larger, the switching occurs earlier and in a narrower time window, in good agreement with the experimental results [
22,
23,
24].
Having obtained the Lorentzian distribution of
, the area fraction of downward domains
S can be approximated as a function of the amplitude
V and duration
t of the applied pulse
where ‘±’ relates to the sign of the applied pulse.
Using Equation (
5) and the same set of
and
, the evolutions of
S as a function of pulse duration
t at different pulse amplitudes
V are calculated and shown in
Figure 2. When the applied pulse voltage is positive, as shown in
Figure 2a, the larger the pulse duration, the larger the fraction of domains that can be switched downward (i.e., the closer
S is to 1). In addition, under the same pulse duration, a larger pulse amplitude produces a larger fraction of downward domains. In contrast, as shown in
Figure 2b, more domains are switched upward, i.e.,
S becomes closer to 0, with the increase in pulse amplitude and duration.
Then, a pulse train (see
Figure 3a) consisting of consecutive pulses whose amplitudes vary from 0.8 v to 2.5 v (in increments of 0.1 v), from 2.4 v to 0.8 v (in decrements of 0.1 v), and subsequently vary from −1.1 v to −3.4 v (in decrements of 0.1 v), from 3.3 v to −1.1 v (in increments of 0.1 v), and finally back to 0 v, is applied to the FTJ, and the resulting evolution of
S is investigated.
Figure 3b shows the evolution of
S as a function of pulse amplitude at different pulse durations. It is seen that
S progressively increases toward a saturated value with the application of positive pulses and then progressively decreases toward the initial value, i.e., 0, with the application of negative pulses. Such hysteretic evolution of
S agrees well with the hysteretic behavior of polarization switching. Additionally,
Figure 3b also shows that, as the pulse duration increases, the saturated value of
S after the application of positive pulses increases. This can be explained by the fact that positive pulses with longer duration can cause more domains to switch downward (see
Figure 2a).
Substituting the
S values into Equation (
1), we can obtain multiple
hysteresis loops, as shown in
Figure 3c. These
hysteresis loops are in good agreement with those observed experimentally in FTJs [
13,
23,
24], confirming the validity of our ferroelectric memristor model. Moreover, note that these loops are not fully closed, because the applied negative pulses are not sufficient to switch all the domains back to the upward direction. Therefore,
S does not return exactly to 0 (not visible in
Figure 3b), leading to the inconsistency between the initial and final
G values.
3. The FM-TCNN
By introducing transiently chaotic dynamics into neural networks, Chen et al. proposed TCNN [
6], which is governed by
where
are the output, the internal state, and the self-feedback connection weight of neuron
i, respectively.
is a steepness parameter of the output function.
is a damping factor of nerve membrane.
is a positive scaling parameter for inputs.
is the connection weight from neuron
j to neuron
i.
is the input bias of neuron
i.
is a positive parameter, and
is a damping factor of the time-dependent
, which is a critical factor for the annealing process. The structural diagram of TCNN is shown in
Figure 4.
Due to the simulated annealing mechanism, the TCNN is an effective method to solve combinatorial optimization problems. As shown in
Figure 5, the single neuron of TCNN can converge to the stable state from the chaotic state after a sufficient number of iterations. More importantly, the critical factor
determines the chaotic attenuation speed. The larger the
, the faster the chaos attenuates. In addition, with the different values of
, the stabilized neuronal outputs are different.
Then, the FTJ-based ferroelectric memristor model, as established in
Section 2, is introduced into TCNN. Specifically, the conductance of the ferroelectric memristor is used as self-feedback connection weight for the annealing process. Equation (
8) is therefore modified as:
where
is a coefficient. In addition, the conductance is adjusted by an applied voltage pulse whose amplitude depends linearly on the neuron output:
where
b is the slope and
c is the intercept. The resultant conductance variation under the applied voltage pulse follows the FTJ’s behavioral model, i.e., Equations (
1)–(
5). Therefore, the self-feedback connection weight can be dynamically modified according to the neuron state. Combining Equations (
1)–(
7), (
9), and (
10), a ferroelectric memristor-based transiently chaotic neural network (FM-TCNN) is thus established.
Note that Equation (
10) can be implemented by an operational amplifier, as shown in
Figure 6. Let R1 =
bR2, R4 =
bR3, V1 =
c, and V2 =
, thus Vout =
. Then, the output voltage Vout of operational amplifier is applied to the FTJ-based ferroelectric memristor M1.
When
, FM-TCNN degrades into a single-neuron model. The parameters are set as
, and the others are the same as the TCNN single-neuron model. The iterations of
are shown in
Figure 7. We can see from
Figure 7a that the neuron presents the chaotic state at the initial stage of iteration. Then, with the decrease in
, as shown in
Figure 7b, the neuron quickly enters the stable state. The greater the value of
b, the slower
decreases, and the slower the neuron stabilizes. What is more, the neuron stabilizes at the same value even when the values of
b are different, which is quite different from the behavior of TCNN.
4. FM-TCNN for TSP
The TSP is one of the typical combinatorial optimization problems. Taking a city as a starting point, the traveling salesman must pass all the other cities and finally return to the starting point. The limitation is that each city can only be visited once. The goal of the TSP is to obtain the shortest path. Based on it, the energy function is first constructed as [
6]:
where
n is the number of cities,
is the elements in matrix
,
,
are weight parameters, and
is the distance between city
i and city
k. Equation (
11) consists of two parts. The first part is the limitation for the effective path, while the last part contains information concerning the total distance. It should be noted that the minimum energy state of the network means the optimal travel route for TSP. Moreover, note that for purely memristive circuits, the Lyapunov function is exactly quadratic, as Equation (
11). Thus, searching the global minimization for the quadratic combinatorial optimization problems (e.g., the TSP) can be mapped directly at the circuit level [
25,
26].
Since the output of TSP is a
matrix, the transient chaotic neural network for the TSP should change to have a double subscript. Thus, the FM-TCNN for
n-city TSP can be obtain as
where the output
of neuron in the network represents to visit city
i in visiting order
j. According to Equations (
12)–(
15), the self-feedback connection weights of different neurons are different; therefore,
ferroelectric memristors are used for the 10-city TSP. The coordinates of the 10-city TSP and its optimal solution are shown in
Figure 8. The initial values of
are generated randomly in the range of
. After iterating Equations (
12)–(
15) for 1800 times, the steady-state outputs can be obtained. A case of the elements in matrix
X is that only
and
are equal to 1, as shown in
Figure 9, and the rest are close to 0, indicating an optimal solution for the 10-city TSP as
. We can see from
Figure 9 that the neurons are explored in the whole phase space at the beginning of the iteration. Then, the neurons enter the periodic state after about 190 iterations and finally enter the stable state after about 500 iterations. As shown in
Figure 10, the evolution of the energy function matches that of the neuron state, and it achieves the minimum after about 500 iterations. The tour length and the corresponding energy function of the optimal solution are approximately equal to 2.6907 and 2.9597, respectively.
As discussed for
Figure 7, the decreasing rate of
and the speed of the neuron entering the steady state are determined by the parameter
b. Thus, the influence of
b on the performance of solving the 10-city TSP is first evaluated here. For each value of
b, 1000 independent experiments are run. The iteration stops if the value of the energy function is less than 4 and the difference between the current energy function and the previous one is less than 0.001, or the number of iterations reaches 1800. The parameters are set as
, and the others are the same as the FM-TCNN single-neuron model. The experimental results are shown in
Table 1. It is seen that all the solutions are feasible in the investigated
b range of
. When
b is negative, a negative disturbance is applied. Accelerated convergence can be achieved at the price of global optimum rate with decreasing
b. By contrast, when
b is positive, a positive disturbance is applied. The convergence speed reduces with increasing
b. However, the global optimum rate reaches a peak value at
and then decreases with increasing
b. There is a clear trade-off between the convergence speed and global optimum rate. When
, the proposed FM-TCNN achieves relatively good performance for the 10-city TSP, demonstrating its effectiveness in solving combinatorial optimization problems.
In addition, the performance comparison between TCNN [
6], HPM-TCNN [
11], and the proposed FM-TCNN is shown in
Table 2. Note that the
value in TCNN and the
b and
c values in HPM-TCNN have already been optimized. The FM-TCNN achieves the highest global optimum rates, which are 2.44% and 13.25% higher than those of the TCNN and HPM-TCNN, respectively. Meanwhile, the FM-TCNN achieves the shortest average path distance, which is the closest to the global optimum value of 2.6907. In addition, the FM-TCNN achieves the fewest average iterations for convergence, and its convergence speed is 32.8% faster than the TCNN. Therefore, the advantages of the FM-TCNN over TCNN and HPM-TCNN, in terms of global optimal rate, average path distance, and convergence speed, are well-verified. We believe that this should be attributed to the good nonlinear performance of the FTJ-based ferroelectric memristor.
We further evaluate the performances of FM-TCNN for 30-city, 50-city, and 75-city TSPs, as shown in
Table 3. The iteration times is set to 30,000. Parameters are the same as above, except
b. Different values of
b are tested, and 50 independent experiments are run for each value of
b. The values of
b which lead to the top five global optimal rates for the 30-city TSP are selected to be reported in
Table 3. It is seen that it becomes more difficult for the FM-TCNN to find the global optimum when the number of cities increases. In addition, the feasible solution rate decreases with the number of cities. More specifically, for the 30-city TSP, the feasible solution rate is over 90%, and the best path length and average path length are close to the global optimum path length, while for the 50-city TSP and 75-city TSP, the gaps between the best path length, average path length, and global optimum path length are relatively large. Additionally, the feasible solution rates decrease dramatically and falls in the range of 60–80%. There is still room for improving the performance of the FM-TCNN on large-size TSPs. Nevertheless, for the same city scale, the best performance of the FM-TCNN can be achieved in a relatively wide range of
b values. This enables us to effectively reduce the time for parameter adjustment when using the FM-TCNN.