已知函数 ，
前向神经网络结构为 1-3-4 -3-1，
模拟退火算法，
实现对网络的权系数进行优化，对函数的模拟。
回答如下问题：
1、 给出神经网络的结构图及函数表达式
2、 简述对应算法
3、 给出学习过程
（1）、算法原理
（2）、学习误差e<0.005时结束
（3）、按学习的次数，按均匀的等分取10组权系数为结果，记录数据并给出这些数据值
（4）、选择5个典型权系数给出这5 个权系数的优化过程的变化曲线
（5）、比较学习结果的神经网络运行结果曲线与函数的输入输出曲线，指出其差异
（6）、讨论学习中存在的问题
（7）、提出来改进方法
（8）、给出源程序
（9）、给出参考文献
1 引言
模拟退火算法(simulated annealing，SA)作为局部搜索算法的扩展，在每一次修改模型的过程中，随机产生一个新的状态模型，然后以一定的概率选择邻域中能量值大的状态．这种接受新模型的方式使其成为一种全局最优算法，并得到理论证明和实际应用的验证．SA虽然在寻优能力上不容置疑，但它是以严密的退火计划为保证的，具体地讲，就是足够高的初始温度、缓慢的退火速度、大量的迭代次数及同一温度下足够的扰动次数[1].
本文只是通过前向神经网络结构，采用模拟退火算法，实现对网络的权系数进行优化，对函数 的图形模拟。第2部分介绍模拟退火算法（Simulated Annealing）的原理。第3部分分析神经元的学习机理。第4部分阐述模拟退火算法的基本过程以及。第5部分给出本题的流程图和解题结果，并对结果分析。第6部分给出改进方法。第7部分对本次设计的小结。
2 模拟退火算法的原理
模拟退火法(Simulated annealing, SA)是一种用于解决连续、有序离散和多模态优化问题的随机优化技术，它是模拟热力学中经典粒子系统的降温过程，来求解规划问题的极值。当孤立粒子系统的温度以足够慢的速度下降时，系统近似处于热力学平衡状态，最后系统将达到本身的最低能量状态，即基态，这相当于能量函数的全局极小点。由于模拟退火法可以有效地摆脱局部极小值，以任意接近于1地概率达到全局最小值点，因而不论在组合优化问题，还是在连续空间优化问题上都有广泛地应用。
模拟退火法的基本原理如下：
(1）给定初始温度T0，及初始点，计算该点的函数值f(x)。
(2）随机产生扰动 ，得到新点 ，计算新点函数值f(x′),及函数值差Δf=f(x′)-f(x)。
(3）若 ，则接受新点，作为下一次模拟的初始点；
(4）若 ，则计算新点接受概率： ，产生伪随机数r, ，如果 ≥r，则接受新点作为下一次模拟的初始点；否则放弃新点，仍取原来的点作为下一次模拟的初始点。
以上步骤称为Metropolis过程。按照一定的退火方案逐渐降低温度，重复Metropolis过程，就构成了模拟退火优化算法，简称模拟退火法。当系统温度足够低时，就认为达到了全局最优状态。按照热力学分子运动理论，粒子作无规则运动时，它具有的能量带有随机性，温度较高时，系统的内能较大，但是对某个粒子而言，它所具有的能量可能较小。因此算法要记录整个退火过程中出现的能量较小的点。 在模拟退火优化算法中，降温的方式对算法有很大影响。如果温度下降过快，可能会丢失极值点；如果温度下降过慢，算法的收敛速度又大大降低。为了提高模拟退火优化算法的有效性，许多学者提出了多种退火方案，有代表性的有：（1）经典退火方式：降温公式为： ；特点是温度下降很缓慢，因此，算法的收敛 速度也是很慢的。
（2）快速退火方式：降温公式为：tk+1=α/log(1+tk),式中α可以改善退火曲线的形态。这种退火方式的特点是在高温区，温度的下降是比较快的，而在低温区，降温的速率较小。这符合热力学分子运动理论中，某粒子在高温时，具有较低能量的概率要比在低温时小得多，因此寻优的重点应在低温区。
模拟退火算法的基本组成可以分为：
(1)产生机制：以一定的概率选择一种解，这个概率函数称为生成函数。
(2)接受准则：以一定的概率接受代价函数值的偶然上升，这个概率函数称为接受函数。
(3)控制参数：以一定的冷却方式退火，实际上用控制参数c进行退火控制，以确定随机引入的随机扰动的强度。
(4)代价函数：也成为目标函数，通过对代价函数的比较来决定当前解是否改变，同时，该函数也用来作为算法结束的判断条件。
模拟退火算法通过产生机制产生组合优化问题解的序列，并由与Metropolis准则对应的转移概率Pi
1
当
其它
=
确定是否接受从当前解i到新解j的转移。式中的c∈R+表示控制参数。开始让c取较大的值（与固体的溶解温度相对应），在进行足够多的转移后，缓慢减小c的值，如此重复，直至满足某个停止准则时算法终止。因此，模拟退火算法可视为递减控制参数值时Metropolis算法的迭代。以上就是模拟退火算法的三函数两准则，即状态产生函数，状态接受函数，温度更新函数，内循环终止准则，外循环终止准则。
3 神经元的学习机理
每个神经元的学习机理如图4.1所示：
∑
θi
f [·]
X1
X2
Xn
wi1
wi2
win
yi
图3.1
用文字描述其过程为：X1，X2，…，Xnn是神经元的输入,θi是i神经元的阀值,wi1,wi2…win分别是i神经元对X1，X2，…，Xn的权系数，yi是i神经元的输出；f [·]是激发函数。
从神经元模型可以得到神经元的数学模型表达式：
令 ，则上式可写成 。即对每个神经元来说
其输入为
输出为
4模拟退火算法过程
标准的模拟退火法的一般步骤可描述如下：
(1)给定初始温度t=t0，随机产生初始状态si=s0 ,令k=0;
(2)Repeat;
(2.1) Repeat;
( 2.1.1 ) 产生新状态sj=Generate(si);
( 2.1.2 ) if min{1,exp[-(C(sj) –C(si)/tk]}>=random[0,1] si=sj ;
( 2.2.3 )Until 抽样稳定准则满足；
(2.2)退温tk＋1=update(tk)并令k=k+1;
(3)Until算法终止准则满足；
(4)输出算法搜索结果。
具体的描述如下：
①对权系数wij置初值。
对各层的全系数wij置一个较小的非零随机数
②输入一组样本
第一步，任选一初始状态S0作为当前解，并设初始温度T0令i=0;
第二步，令T=Ti,调用Metropolis抽样算法，返回其当前解S,令Si=S;
第三步，按一定方式降低温度，即T=Ti+1 , （Ti+1<Ti ,i=I+1）;
第四步，检查退火过程是否结束，未结束转第二步，结束转第五步；
第五步，以当前解Si作为最优解输出.
上述Metropolis抽样算法如下：
（1） K=0,S(k)=S,在T下执行以下步骤：
（2） 根据当前解S(k),所处状态S，产生一个近邻子集N[S(k)] S, 由N[S(k)]随机地选一个新状态S，作为下一个候选解，计算ΔC=C(S’)-C(S(k))
（3） 若ΔC<0,则接受S’为下一个当前解；若ΔC〉0，则按概率exp(-ΔC/T)接受T’为下一个当前解。若S’被接受，则S(k+1)=S,否则S(k+1)=S（k）。
（4） K=K+1,由某个给定的收敛准则，判断算法是否应该中止。如果是，则转第5步，否则转第2步。
（5） 将当前解返回给调用它的退火过程算法。
③程序流程图
图4.1
5 解题过程与结果分析
5.1 函数 ,
如图5.1
如图5.1
5.2 神经网络结构是 1-3-4 -3-1，
图5.2
5.3 权系数与阀值的选定
初始权系数与阀值的产生
double tempW = 2*((double)(rand()000))/10000.0-1; //随机产生值
权系数与阀值的产生是通过伪随机数生成，范围在 之间。
这里随机给出一组数据：
0.449000 -0.507000 -0.859400 0.854400 -0.179600 0.628200 0.401600 0.625000 0.626600 -0.487400 0.528800 -0.624400 -0.080400 -0.342400 0.845800 -0.737400 -0.340600 -0.271200 0.752000 -0.018400 -0.969800 -0.045000 0.847000 0.135200 0.515600 -0.033200 0.699600 -0.072200 0.722000 -0.774200 -0.638000 -0.955400 0.577000 0.183800 -0.825400 -0.668400 0.755000 -0.449000 -0.235800 0.683400 0.864400
其中第4，5，6，19，20，21，22，35，36，37，41个数据是阀值。
5.4 通过网络计算输出
/************************************************************************/
/* 激发函数,用于计算下层神经网络的输出 */
/************************************************************************/
double ActivateFunction(double x)
{
return 1/(1+exp(-x));
}
/************************************************************************/
/*在设定状态下,计算输出 */
/************************************************************************/
double ComputeNetworkOutput(double* Status, double Xi)
{
double tempW[NUMOFW];
double tempX[NUMOFW];
tempW[4] = Status[1]*Xi-Status[4];
tempW[5] = Status[2]*Xi-Status[5];
tempW[6] = Status[3]*Xi-Status[6];
tempX[4] = ActivateFunction(tempW[4]);
tempX[5] = ActivateFunction(tempW[5]);
tempX[6] = ActivateFunction(tempW[6]);
tempW[19] = Status[7]*tempX[4]+Status[8]*tempX[5]+Status[9]*tempX[6]-Status[19];
tempW[20] = Status[10]*tempX[4]+Status[11]*tempX[5]+Status[12]*tempX[6]-Status[20];
tempW[21] = Status[13]*tempX[4]+Status[14]*tempX[5]+Status[15]*tempX[6]-Status[21];
tempW[22] = Status[16]*tempX[4]+Status[17]*tempX[5]+Status[18]*tempX[6]-Status[22];
tempX[19] = ActivateFunction(tempW[19]);
tempX[20] = ActivateFunction(tempW[20]);
tempX[21] = ActivateFunction(tempW[21]);
tempX[22] = ActivateFunction(tempW[22]);
tempW[35] = Status[23]*tempX[19]+Status[24]*tempX[20]+Status[25]*tempX[21]+Status[26]*tempX[22]-Status[35];
tempW[36] = Status[27]*tempX[19]+Status[28]*tempX[20]+Status[29]*tempX[21]+Status[30]*tempX[22]-Status[36];
tempW[37] = Status[31]*tempX[19]+Status[32]*tempX[20]+Status[33]*tempX[21]+Status[34]*tempX[22]-Status[37];
tempX[35] = ActivateFunction(tempW[35]);
tempX[36] = ActivateFunction(tempW[36]);
tempX[37] = ActivateFunction(tempW[37]);
tempW[41] = Status[38]*tempX[35]+Status[39]*tempX[36]+Status[40]*tempX[37]-Status[41];
tempX[41] = ActivateFunction(tempW[41]);
return tempX[41];
}
5.5 权系数与阀值的优化过程取样
这里分别取出第1，1001，2003，3000，4012，5002，6035，7006，7422次迭代之后权系数与阀值的图形。
图 5.5.1 图 5.5.2
图5.5.3 图5.5.4
图5.5.5 图5.5.6
图5.5.7 图5.5.8
以上各图对应得数据，请参看附本”weight.dat”
根据图形可以明显观察出，权系数在最初的变化比较大，但在后期，特别在温度变化很小的情况下，图形的变化很小。这就说明本程序表现的特征比较符合模拟退火算法。
5.6 结果分析
编号
1
2
3
4
5
6
输入x
0
0.1
0.2
0.3
0.4
0.5
原始输出y
0
0.000998
0.007947
0.026597
0.062307
0.119856
目标输出y0
0.00246
0.003731
0.006655
0.013767
0.031226
0.071413
差值(y0-y)
0.00246
0.002733
-0.00129
-0.01283
-0.03108
-0.04844
编号
7
8
9
10
11
12
输入x
0.6
0.7
0.8
0.9
1
1.1
原始输出y
0.203271
0.315667
0.459108
0.634495
0.841471
1.078361
目标输出y0
0.152279
0.288838
0.481395
0.712969
0.958491
1.195856
差值(y0-y)
-0.05099
-0.02683
0.022287
0.078474
0.11702
0.117495
样本编号
13
14
15
16
17
18
输入x
1.2
1.3
1.4
1.5
1.6
1.7
原始输出y
1.342136
1.628413
1.931481
2.244364
2.558908
2.865911
目标输出y0
1.411186
1.598567
1.757544
1.890646
2.001693
2.094845
差值(y0-y)
0.06905
-0.02985
-0.17394
-0.35372
-0.55722
-0.77107
编号
19
20
21
22
23
24
输入x
1.8
1.9
2
2.1
2.2
2.3
原始输出y
3.155266
3.416143
3.63719
3.806753
3.913123
3.944781
目标输出y0
2.174172
2.243509
2.30647
2.366546
2.427231
2.492185
差值(y0-y)
-0.98109
-1.17263
-1.3307
-1.44021
-1.48589
-1.4526
编号
25
26
27
28
29
30
输入x
2.4
2.5
2.6
2.7
2.8
2.9
原始输出y
3.890668
3.740451
3.484789
3.115599
2.626307
2.012087
目标输出y0
2.565404
2.65143
2.755588
2.884273
3.045295
3.248288
差值(y0-y)
-1.32526
-1.08902
-0.7292
-0.23133
0.418988
1.236201
编号
31
32
33
34
35
36
输入x
3
3.1
3.2
3.3
3.4
3.5
原始输出y
1.27008
0.39959
0.597751
1.717851
2.954055
4.297095
目标输出y0
3.50517
3.83059
4.242289
4.761086
5.410101
6.212669
差值(y0-y)
2.23509
3.431
3.644538
3.043235
2.456046
1.915574
编号
37
38
39
40
41
42
输入x
3.6
3.7
3.8
3.9
4
4.1
原始输出y
5.735065
7.253457
8.835228
10.460923
12.10884
13.755238
目标输出y0
7.188293
8.346381
9.678587
11.152284
12.70917
14.272196
差值(y0-y)
1.453228
1.092924
0.843359
0.691361
0.600329
0.516958
编号
43
44
45
46
47
48
输入x
4.2
4.3
4.4
4.5
4.6
4.7
原始输出y
15.374597
16.939908
18.423016
19.794985
21.026502
22.088305
目标输出y0
15.760202
17.104646
18.261258
19.212633
19.963027
20.529578
差值(y0-y)
0.385605
0.164738
-0.161758
-0.582352
-1.063475
-1.558727
编号
49
50
51
52
53
54
输入x
4.8
4.9
5
5.1
5.2
5.3
原始输出y
22.951633
23.588687
23.973107
24.08044
23.888614
23.378392
目标输出y0
20.933799
21.195243
21.327423
21.33524
21.213062
20.942788
差值(y0-y)
-2.017834
-2.393444
-2.645684
-2.7452
-2.675552
-2.435604
编号
55
56
57
58
59
60
输入x
5.4
5.5
5.6
5.7
5.8
5.9
原始输出y
22.533812
21.342595
19.796522
17.891773
15.629217
13.014647
目标输出y0
20.491986
19.813009
18.845831
17.529344
15.826149
13.759318
差值(y0-y)
-2.041826
-1.529586
-0.950691
-0.362429
0.196932
0.744671
编号
61
62
63
输入x
6
6.1
6.2
原始输出y
10.058958
6.778267
3.193957
目标输出y0
11.443003
9.074857
6.875808
差值(y0-y)
1.384045
2.29659
3.681851
以上数据对应的图形为图5-1
图5-1
说明一下图中给定的精度是总体样本差，因此当样本数大于2400是就可满足题中的0.005的单位样本差，但精度要求对于100的起始温度而言太高，所以在最后是通过外循环的“零度准则”（外循环终止准则）结束程序的运行。
另外，如图5.2，模拟图形与原图形相差就很小了。但数据没有保留，而且每次即使用同样的参数图形的输出也会相差比较大，学习的效果在多次试验中呈现明显的不确定性，而且数据一样则是不可预期的变化趋势，我想这对于计算机的算法而言是个亟待解决的问题。但我们知道，模拟退火算法的核心思想就是随机优化，也就是说在大量计算的情况下，满足真正的随机平衡过程时，结果才有可能往理想的变化趋势进展。经过同一组数据的多次试验，我们可以看出，该算法对图形的上凸学习效果要远远好于下凸的学习效果，而且即使下凸学习偏理想，也会导致下面图形的学习效果。最初，我以为是对 (函数的二阶导数)敏感，而对 则失去感觉，但在此函数的二阶导数只是相差一个符号，不应该会产生这么大的差距（具体什么原因，我也没弄明白）。
图5.2
还有一点，请参看图5.3。初始温度是5.0 内循环次数是5000，样本数是41。看上去图5.1的学习效果与图5.3的学习效果相差不大，图5.1用了7422次迭代，图5.3则用了10087次迭代。
在数据相差较大的情况下，学习效果却相差不大，这点虽比较符合随机优化中的随机思想，但我觉得与算法的实际应用而言，这是个需要该进之处。
图5.3
6模拟退火方法的不足与改进
6.1 不足之处
通过分析模拟退火法的基本原理，可知，模拟退火法在一系列递减温度下产生的点列，从理论上讲，可以看作是一系列的马尔可夫链（Markov Chains）。在某一温度下多次重复 Metropolis过程，目标函数值的分布规律将满足玻尔兹曼分布规律。如果系统温度以足够慢的速率下降，玻尔兹曼分布就趋向收敛于全局最小状态的均匀分布。也就是说，若按照一定的条件产生无限长的马尔可夫链，模拟退火法就能保证以概率1.0收敛于全局极小点，该结论已被文献［4］证明。文献［4］指出，在某一温度下，只要计算时间足够长，也就是马尔可夫链足够长，其起始点的函数值将以很高的概率低于终止点的函数值，即求得全局最小点。
通过上面的分析可以看出，模拟退火法主要存在以下不足：
（1） 尽管理论上只要计算时间足够长，模拟退火法就可以保证以概率1.0收敛于全局最优点。但是在实际算法的实现过程中，由于计算速度和时间的限制，在优化效果和计算时间二者之间存在矛盾，因而难以保证计算结果为全局最优点，优化效果不甚理想。
（2） 在每一温度下很难判定是否达到了平衡状态，即马尔可夫链的长度不易控制，反应到算法上，就是Metropolis过程的次数不易控制。
（3） 模拟退火算法中的两种退火方式，T始终按照优化前给定的规律变化而没有修正，这是不科学的。
为了解决模拟退火法存在的不足，台湾学者Feng-Tse Lin等人曾提出了用遗传算法的优点来改进模拟退火算法的模拟退火-遗传算法。本文将在退火过程中综合考虑计算结果及搜索路径而确定T的取值方法的自适应模拟退火法与遗传算法相结合，形成自适应模拟退火-遗传算法，从而改进模拟退火算法本身所存在的不足。
6.2 算法的改进
这里也主要引用相关文献的介绍，改进主要从两个方面进行。
6.2.1 改进退火过程
(1）给定初始温度To,随机产生初始状态s，令初始最优解S’=S,当前的状态为S(0)＝S, i,p=0;
(2)令T=Ti,令当前的状态S(i)=S’(k);
(3)判断C(S*<C(S*’)？若是，则令p=p+1;否则，令S*=S*’,p=0；
（4）退温Ti+1=update(Ti)，令i=i+1;
(5)判断p>step2?若是，则转到第6步，否则，返回第2步；
6.2.2 改进的抽样过程
(1) 令k=0时的初始当前状态S’(0)=S(i),初始S*’=S*,q=0
(2) 由状态i通过状态产生函数产生新状态S’,计算增量 .
(3) 若 ，则接受S’作为当前解，并判断C(S*’)>C(S’)?如果是，则令S*’=S’,q=0;
否则，令q=q+1；如果 ,则以概率exp(- )接受S’作为下一当前状态。若S’被接受，令S’(k+1)=S’, q=q+1;否则，令S’(k+1)=S’(k).
(4) 令k=k+1,判断q>step(1)?如果是，就转(5),否则返回(2).
(5) 将当前最优解S*’返回到改进退火过程。
7 小结
通过对模拟退火算法的课程设计，从中学到了一些重要的算法思想和方法
模拟退火算法的随机优化思想，随机在目标明确的问题优化中的应用，这实际是由一种不确定性，按一种确定的方向进行变化，由不确定走向确定的过程。模拟退火的关键是掌握问题的目标和优化的对象，在理清楚这两个方面后，模拟退火的核心就是在不断降温的逼近目标的过程中，通过一定的随机策略和控制策略，实现对对象的优化，降温与优化是个并行过程。
对于SA，三函数两准则就是算法的组成部分。这里主要对状态产生函数谈点感想。状态产生函数是应用一定的随机函数对当前的状态进行修改，本题是对权系数与阀值的修改，我认为这里体现了算法的精华，另一处就是 时的接受函数，当然降温函数也是个重要部分。比较而言，我认为状态产生函数最值得研究。由于有许多组数据需要优化，而每次只能对权系数与阀值在一组数据上进行优化，即对权系数与阀值的优化是对随机选择的一组数据优化。那如何对整个数据进行优化呢？整个数据，则会有每组数据都产生期望输出与实际输出（或称目标输出）的差异，因而每组数据都会生成代价，此代价对于单个数据而言是减小的，当对于取样的整个数据而言，可能是增加。因而，我认为应该计算整体的代价，即整体样本才是本题的优化对象。
在编程过程中，不断修改产生函数，最初选择的是 ,经试验，收敛太慢，后来改用 ，效果要好。最初的温度下降函数选用的t(i+1)=t0/(1+log(currentTemperature)),经比较也不如t(i+1)=t(i)*0.95,虽然前者在前期温度的变化较明显，但后来出于几乎不变的状态，这对不断降温明显是个设计问题。
另外，为了方便演示，我将整个程序设计成GUI，可以方便查看神经网络结构，修改精度，初始温度，内循环数以及样本数。同时，温度与权系数的变化也可以在界面上动态更新。该小程序还利用了matcom实现画图，将最终结果显示在画面上，由于权系数变化太快，不好随时更新，于是数据保留在文件中，最终只是以一定的间隔抽取，通过matcom画图。这次设计基本完成了题目给定的目标，还有许多地方要加以改进。
[1] 模拟退火算法机理研究 陈华根 吴健生 王家林 陈 冰 同济大学海洋地质教育部重点实验室，上海20009
[2] 问题求解的人工智能神经网络方法 王士同陈剑夫等编著 气象出版社
[3] 神经网络模糊逻辑控制 余永权著 电子工业出版社
[4] 智能优化算法及其应用 王凌 清华大学出版社，施普林格出版社，2001
[5] http://bbs.matwav.com/index.html
[6]http://www.programsalon.com/search_db.asp?keyword=
