模拟退火

来源：互联网发布：怎么操作淘宝千人千面编辑：程序博客网时间：2024/06/10 08:57

问题：

已知函数，

前向神经网络结构为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）给定初始温度T₀，及初始点，计算该点的函数值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准则对应的转移概率P_i

当

其它

确定是否接受从当前解i到新解j的转移。式中的c∈R⁺表示控制参数。开始让c取较大的值（与固体的溶解温度相对应），在进行足够多的转移后，缓慢减小c的值，如此重复，直至满足某个停止准则时算法终止。因此，模拟退火算法可视为递减控制参数值时Metropolis算法的迭代。以上就是模拟退火算法的三函数两准则，即状态产生函数，状态接受函数，温度更新函数，内循环终止准则，外循环终止准则。

3 神经元的学习机理

每个神经元的学习机理如图4.1所示：

∑

θi

f [·]

X₁

X₂

X_n

w_i1

w_i2

w_in

y_i

图3.1

用文字描述其过程为：X₁，X₂，…，X_nn是神经元的输入,θ_i是i神经元的阀值,w_i1,w_i2…w_in分别是i神经元对X₁，X₂，…，X_n的权系数，y_i是i神经元的输出；f [·]是激发函数。

从神经元模型可以得到神经元的数学模型表达式：

令，则上式可写成。即对每个神经元来说

其输入为

输出为

4模拟退火算法过程

标准的模拟退火法的一般步骤可描述如下：

(1)给定初始温度t=t₀_，随机产生初始状态s_i=s₀ ,令k=0;

(2)Repeat;

(2.1) Repeat;

(2.1.1) 产生新状态s_j=Generate(s_i);

(2.1.2) if min{1,exp[-(C(s_j) –C(s_i)/t_k]}>=random[0,1] s_i=s_j;

(2.2.3)Until 抽样稳定准则满足；

(2.2)退温t_k_＋1=update(t_k)并令k=k+1;

(3)Until算法终止准则满足；

(4)输出算法搜索结果。

具体的描述如下：

①对权系数w_ij置初值。

对各层的全系数w_ij置一个较小的非零随机数

②输入一组样本

第一步，任选一初始状态S₀作为当前解，并设初始温度T₀令i=0;

第二步，令T=T_i,调用Metropolis抽样算法，返回其当前解S,令Si=S;

第三步，按一定方式降低温度，即T=T_i+1, （T_i+1<T_i ,i=I+1）;

第四步，检查退火过程是否结束，未结束转第二步，结束转第五步；

第五步，以当前解S_i作为最优解输出.

上述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.2 神经网络结构是1-3-4-3-1，

图5.2

5.3 权系数与阀值的选定

初始权系数与阀值的产生

double tempW = 2*((double)(rand()%10000))/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 结果分析

编号

输入x

0.1

0.2

0.3

0.4

0.5

原始输出y

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

编号

输入x

0.6

0.7

0.8

0.9

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

样本编号

输入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

编号

输入x

1.8

1.9

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

编号

输入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

编号

输入x

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

编号

输入x

3.6

3.7

3.8

3.9

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

编号

输入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

编号

输入x

4.8

4.9

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

编号

输入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

编号

输入x

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=%C9%F1%BE%AD%CD%F8%C2%E7&search_type=code

//程序源码：以下是主要代码，由于涉及界面，消息传递，多线程等附加部分，加上会有碍//视听。原始输入（original.dat）与目标输出(NeuroNetwork.dat)，还有每次迭代变化的//权系数与阀值(weight.dat)都存储在文件中。

#include "StdAfx.h"

#include "saa.h"

/************************************************************************/

/*函数名:SetOriginalX 功能:根据函数定义域和步长设置X */

/************************************************************************/

void SetOriginalX(double* OriginalX)

{

double dStep = 0;

dStep = (2*PI-0)/(NUMOFDEVIDE-1);

for (int i = 1;i < NUMOFDEVIDE;i++)

{

OriginalX[i] = (i-1)*dStep;

}

/************************************************************************/

/* 功能:根据f(x)=|sin(x)|*x^2 计算y */

/************************************************************************/

void ComputerY(double* X, double* Y)

{

for (int i = 1;i < NUMOFDEVIDE;i++)

{

Y[i] = fabs(sin(X[i]))*X[i]*X[i];

}

/************************************************************************/

/* 功能:将x归一化 */

/************************************************************************/

void TransformX(double* OriginalX, double* StandardX)

{

for (int i = 1;i < NUMOFDEVIDE;i++)

{

StandardX[i] = OriginalX[i] / (2*PI);

}

/************************************************************************/

/* 功能:将y归一化 */

/************************************************************************/

void TransformY(double* OriginalY, double* StandardY)

{

for (int i = 1;i < NUMOFDEVIDE;i++)

{

StandardY[i] = OriginalY[i] / (Y_MAX-Y_MIN);

}

/************************************************************************/

/* 初始化温度 */

/************************************************************************/

double InitTemperature(double T)

{

return T;

}

/************************************************************************/

/* 初始化状态 */

/************************************************************************/

void InitStatus(double* Status)

{

for (int i = 1; i < NUMOFW; i++)

{

double tempW = 2*((double)(rand()%10000))/10000.0-1;

Status[i] = tempW;

}

/************************************************************************/

/* 保存状态 */

/************************************************************************/

void StoreStatus(double* Status, double* DestStatus)

{

for (int i = 1; i < NUMOFW; i++)

{

DestStatus[i] = Status[i];

}

/************************************************************************/

/* 产生函数 */

/************************************************************************/

double GenerateFunction(double Cur1Status, double CurTemperature)

{

// double newValue = 0.0;

// double deltaW = exp(-Cur1Status*Cur1Status/CurTemperature);

// newValue = Cur1Status + deltaW;

// return newValue;

double newValue = 0.0;

double Random = ((double)(rand()%10000))/10000.0;

double deltaW = PI*tan(Random-0.5);//exp(-Cur1Status*Cur1Status/CurTemperature);

newValue = Cur1Status + deltaW;

return newValue;

}

/************************************************************************/

/* 产生新状态 */

/************************************************************************/

void GenerateNewStatus(double* CurStatus, double* NewStatus,double CurTemperature)

{

StoreStatus(CurStatus, NewStatus);

int nIndex;

while((nIndex = rand() % NUMOFW)==0);

NewStatus[nIndex] = GenerateFunction(CurStatus[nIndex],CurTemperature);

// int numToChange;

// StoreStatus(CurStatus, NewStatus);

// while((numToChange = rand() % NUMOFW)==0);

// for (int i = 1; i < numToChange;i++)

// {

// int nIndex;

// while((nIndex = rand() % NUMOFW)==0);

// NewStatus[nIndex] = GenerateFunction(CurStatus[i],CurTemperature);

// }

}

/************************************************************************/

/* 激发函数,计算下层神经网络的输出 */

/************************************************************************/

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];

}

/************************************************************************/

/* 计算新状态与当前状态总代价差 */

/************************************************************************/

double CalcTotalCost(double* NewStatus, double* CurStatus, double* StandardX, double* StandardY)

{

double NewCost = 0.0;

double CurCost = 0.0;

double CurY = 0.0;

double NewY = 0.0;

double deltaTotal = 0.0;

double CurDelta = 0.0;

double NewDelta = 0.0;

for (int i = 1; i < NUMOFDEVIDE; i++)

{

CurY = ComputeNetworkOutput(CurStatus, StandardX[i]); //计算出单个x经过神经网络的输出

CurCost = fabs(CurY-StandardY[i]);

NewY = ComputeNetworkOutput(NewStatus, StandardX[i]);

NewCost = fabs(NewY-StandardY[i]);

CurDelta += CurCost;

NewDelta += NewCost;

}

deltaTotal = NewDelta - CurDelta;

return deltaTotal;

}

/************************************************************************/

/* 计算与目标理想的差 */

/************************************************************************/

double CalcTotalDiffer(double *CurStatus,double* StandardX, double* StandardY)

{

double differ = 0.0;

double tempY;

double totalDiffer = 0.0;

for (int i = 1; i < NUMOFDEVIDE; i++)

{

tempY = ComputeNetworkOutput(CurStatus, StandardX[i]);

differ = fabs(tempY-StandardY[i]);

totalDiffer += differ;

}

return totalDiffer;

}

/************************************************************************/

/* 接受更新 */

/************************************************************************/

void AcceptNewStatus(double* NewStatus, double* CurStatus)

{

StoreStatus(NewStatus,CurStatus);

}

/************************************************************************/

/* 计算更新概率 */

/************************************************************************/

double UpdateProbability(double deltaTotalDiffer, double CurTemperature)

{

return min(1,1/(1+exp(deltaTotalDiffer/CurTemperature)));

}

/************************************************************************/

/* 内循环终止判断 */

/************************************************************************/

int InternalLoopOver(int Mode,int InternalIterNumber)

{

// switch(Mode)

// {

// case 0:

// if (InternalIterNumber > NUMOFITER)

// {

// return 1;

// }

// else

// {

// return 0;

// }

// break;

// default:

// return 0;

// }

if (InternalIterNumber > NUMOFITER)

{

return 1;

}

else

{

return 0;

}

/************************************************************************/

/* 外循环终止判断 */

/************************************************************************/

int ExternalLoopOver(int Mode,double CurTemperature, double TotalDiffer)

{

// switch(Mode)

// {

// case 0:

// if (CurTemperature <= 0.01)

// {

// return 1;

// }

// else

// {

// return 0;

// }

// //break;

// case 1:

// if ((deltaTotalDiffer < 0.05)||CurTemperature <= 0.01)

// {

// return 1;

// }

// else

// {

// return 0;

// }

// //break;

// default:

// return 0;

// }

if (TotalDiffer < PRECISION ||CurTemperature <= 0.01)

{

return 1;

}

else

{

return 0;

}

double UpdateTemperature(double CurTemperature,int Mode)

{

// switch(Mode)

// {

// case 0:

// return 0.95*CurTemperature;

// //break;

// default:

// return 0.01;

// }

return 0.95*CurTemperature;

// return TEMPERATURE/(1+log(CurTemperature));

}

void Output(int Mode, double* Status,double Temperature,int ExterIterNum)

{

cout<<"迭代次数为:"<<ExterIterNum<<endl;

switch(Mode)

{

case 0:

cout<<"初始权系数与阀值为:"<<endl;

break;

case 1:

cout<<"最终权系数与阀值为:"<<endl;

break;

default:

break;

}

for (int i = 1; i < NUMOFW; i++)

{

cout<<Status[i]<<endl;

}

switch(Mode)

{

case 0:

cout<<"初始温度为:"<<endl;

break;

case 1:

cout<<"最终温度为:"<<endl;

break;

default:

break;

}

cout<<Temperature<<endl;

}

/************************************************************************/

/* 存储权系数到文件 */

/************************************************************************/

void StoreWeight2File(double* Status,FILE *fp)

{

for(int i = 1;i<NUMOFW; i++)

{

fprintf(fp,"/t%f",Status[i]);

}

fprintf(fp,"/n");

}

void StoreWeight2Ctrl(double* Status,CString &str)

{

CString strtmp;

for(int i = 1;i<NUMOFW; i++)

{

strtmp.Format(",%f",Status[i]);

str+=strtmp;

}

double RandomFunc()

{

return ((double)(rand()%10000))/10000.0;

}

/************************************************************************/

/* 主函数:sa算法的主体 */

/************************************************************************/

UINT SAAMain(LPVOID pParam)

{

double* OriginalX0 = new double[NUMOFDEVIDE]; //保存原始的X

double* OriginalY0 = new double[NUMOFDEVIDE]; //保存原始的Y

double* StandardX0 = new double[NUMOFDEVIDE]; //保存归一化后的原始X

double* StandardY0 = new double[NUMOFDEVIDE]; //保存归一化后的原始Y

double T0 = 0; //保存初始温度

double CurrentTemperature = 0;

double Status0[NUMOFW]; //保存最初的状态(权系数与阀值)

double CurrentStatus[NUMOFW]; //当前状态(权系数与阀值)

double NewStatus[NUMOFW]; //新状态

double Probability = 0;

double Random = 0;

double TotalDiffer = 0.0;

int InternalIterNumber = 0;

int ExternalIterNumber = 0;

SetOriginalX(OriginalX0);

ComputerY(OriginalX0,OriginalY0);

TransformX(OriginalX0,StandardX0);

TransformY(OriginalY0,StandardY0);

// char* pszWeightFileName = "Weight.dat";

// CFile* pWeightFile(pszWeightFileName, CFile::modeCreate|CFile::modeWrite);

FILE* pf;

char* pFileName = "Original.dat";

FILE* pf1;

char* pFileName1 = "NeuroObject.dat";

FILE* pf_w;

char* pFileName_w = "Weight.dat";

if(((pf=fopen( pFileName, "wt")) != NULL))

{

cout <<"File open or create error!/n"<<endl;

}

if(((pf1=fopen( pFileName1, "wt")) != NULL))

{

cout <<"File open or create error!/n"<<endl;

}

if(((pf_w=fopen( pFileName_w, "wt")) != NULL))

{

cout <<"File open or create error!/n"<<endl;

}

CString strTemp;

T0 = InitTemperature(TEMPERATURE);

CurrentTemperature = T0;

SendMessage((HWND)pParam,MYWM_UPDATET,(WPARAM)(&CurrentTemperature),(LPARAM)0);

srand((unsigned)time(NULL));

InitStatus(CurrentStatus); //初始化最初状态

StoreStatus(CurrentStatus,Status0); //保存最初的状态

strTemp.Format("当前的温度为:%f/r/n迭代次数为:%d/r/n当前权系数为:/r/n",CurrentTemperature,

(InternalIterNumber+NUMOFITER*ExternalIterNumber));

// int nLength = strTemp.GetLength();

// char *sz = new char[nLength];

// sz = strTemp.GetBuffer(0);

// pWeightFile->Write(sz,nLength);

// CListBox *pListBox = (CListBox*)pParam;

// HWND hListBox = pListBox->GetSafeHwnd();

// fprintf(pf_w,"当前的温度为:%f/n迭代次数为:%d/n当前权系数为:/n",CurrentTemperature,

// (InternalIterNumber+NUMOFITER*ExternalIterNumber));

// StoreWeight2Ctrl(CurrentStatus,strTemp);

// ::SendMessage((HWND)pParam, MYWM_INFOEDIT, (WPARAM)(&strTemp), (LPARAM)0);

::SendMessage((HWND)pParam, MYWM_INFOEDIT, (WPARAM)(CurrentStatus), (LPARAM)0);

StoreWeight2File(CurrentStatus,pf_w);

TotalDiffer = CalcTotalDiffer(CurrentStatus,StandardX0,StandardY0);

if (TotalDiffer < PRECISION)

{

fclose(pf);

fclose(pf1);

fclose(pf_w);

return 0;

}

/* StoreWeight2File(CurrentStatus,pWeightFile);*/

// Output(0,Status0,T0,ExternalIterNumber);

/* T0 = InitTemperature();*/

while (1)

{

double deltaTotal = 0.0;

while (1)

{

GenerateNewStatus(CurrentStatus, NewStatus,CurrentTemperature);

deltaTotal = CalcTotalCost(NewStatus,CurrentStatus,StandardX0,StandardY0);

if (deltaTotal <= 0)

{

AcceptNewStatus(NewStatus,CurrentStatus); //接受修改

// strTemp.Format("当前的温度为:%f/r/n迭代次数为:%d/r/n当前权系数为:/n",CurrentTemperature,

// (InternalIterNumber+NUMOFITER*ExternalIterNumber));

// StoreWeight2Ctrl(CurrentStatus,strTemp);

// ::SendMessage((HWND)pParam, MYWM_INFOEDIT, (WPARAM)(&strTemp), (LPARAM)0);

::SendMessage((HWND)pParam, MYWM_INFOEDIT, (WPARAM)(CurrentStatus), (LPARAM)0);

fprintf(pf_w,"当前的温度为:%f/n迭代次数为:%d/n当前权系数为:/n",CurrentTemperature,

(InternalIterNumber+NUMOFITER*ExternalIterNumber));

StoreWeight2File(CurrentStatus,pf_w);

}

else

{

Probability = UpdateProbability(deltaTotal,CurrentTemperature);

Random = RandomFunc();

if (Probability > Random)

{

AcceptNewStatus(NewStatus,CurrentStatus); //接受修改

// strTemp.Format("当前的温度为:%f/r/n迭代次数为:%d/r/n当前权系数为:/n",CurrentTemperature,

// (InternalIterNumber+NUMOFITER*ExternalIterNumber));

// StoreWeight2Ctrl(CurrentStatus,strTemp);

// ::SendMessage((HWND)pParam, MYWM_INFOEDIT, (WPARAM)(&strTemp), (LPARAM)0);

::SendMessage((HWND)pParam, MYWM_INFOEDIT, (WPARAM)(CurrentStatus), (LPARAM)0);

fprintf(pf_w,"当前的温度为:%f/n迭代次数为:%d/n当前权系数为:/n",CurrentTemperature,

(InternalIterNumber+NUMOFITER*ExternalIterNumber));

StoreWeight2File(CurrentStatus,pf_w);

}

//else 当前状态不变

}

if (InternalLoopOver(0,InternalIterNumber))

{

break; //跳出内循环

}

else

{

InternalIterNumber++;

}

TotalDiffer = CalcTotalDiffer(CurrentStatus,StandardX0,StandardY0);

if (TotalDiffer < PRECISION)

{

break;

}

// TotalDiffer = CalcTotalDiffer(CurrentStatus,StandardX0,StandardY0);

if (ExternalLoopOver(1,CurrentTemperature,TotalDiffer))

{

break; //跳出外循环

}

else

{

CurrentTemperature=UpdateTemperature(CurrentTemperature,0);

SendMessage((HWND)pParam,MYWM_UPDATET,(WPARAM)(&CurrentTemperature),(LPARAM)0);

ExternalIterNumber++;

InternalIterNumber = 0;

}

// char* pszOriginalFileName = "Original.dat";

// CFile OriginalFile(pszOriginalFileName, CFile::modeCreate|CFile::modeWrite);

// char* pszObjectFileName = "NeuroObject.dat";

// CFile ObjectFile(pszObjectFileName, CFile::modeCreate|CFile::modeWrite);

// double dStep = 0.0;

// double dY = 0.0;

// for (;dStep < 2*PI; dStep += 0.1)

// {

// dY = fabs(sin(dStep))*dStep*dStep;

// strTemp.Format("/t%f",dY);

// int nLength = strTemp.GetLength();

// char *sz = new char[nLength];

// sz = strTemp.GetBuffer(0);

// OriginalFile.Write(sz,nLength);

// delete [] sz;

// }

// strTemp.Format("/n");

// nLength = strTemp.GetLength();

// sz = new char[nLength];

// sz = strTemp.GetBuffer(0);

// OriginalFile.Write(sz,nLength);

// delete [] sz;

// for (;dStep < 2*PI; dStep += 0.1)

// {

// dY = ComputeNetworkOutput(CurrentStatus, dStep/(2*PI))*(Y_MAX-Y_MIN);;

// strTemp.Format("/t%f",dY);

// int nLength = strTemp.GetLength();

// char *sz = new char[nLength];

// sz = strTemp.GetBuffer(0);

// OriginalFile.Write(sz,nLength);

// delete [] sz;

// }

// strTemp.Format("/n");

// nLength = strTemp.GetLength();

// sz = new char[nLength];

// sz = strTemp.GetBuffer(0);

// OriginalFile.Write(sz,nLength);

// delete [] sz;

// OriginalFile.Close();

// ObjectFile.Close();

// pWeightFile->Close();

/* Output(1,CurrentStatus,CurrentTemperature,ExternalIterNumber);*/

SendMessage((HWND)pParam,MYWM_DRAWPIC,(WPARAM)(CurrentStatus),(LPARAM)0);

double dStep = 0.0;

double dY = 0.0;

for (;dStep < 2*PI; dStep += 0.1)

{

dY = fabs(sin(dStep))*dStep*dStep;

fprintf(pf,"/t%f",dY);

}

for (dStep =0.0;dStep < 2*PI; dStep+=0.1)

{

dY = ComputeNetworkOutput(CurrentStatus, dStep/(2*PI))*(Y_MAX-Y_MIN);

fprintf(pf1,"/t%f",dY);

}

// cout <<"/nTotalDiffer:"

//<<CalcTotalDiffer(CurrentStatus,StandardX0,StandardY0)<<endl;

// double y = fabs(sin(PI/9))*(PI/9)*(PI/9)/(Y_MAX-Y_MIN);

// double y0 = ComputeNetworkOutput(CurrentStatus, 1.0/18);

// cout<< "Differ:"<<(y-y0)<<endl;

fclose(pf);

fclose(pf1);

fclose(pf_w);

SendMessage((HWND)pParam,MYWM_EXITBTN,(WPARAM)(&CurrentTemperature),(LPARAM)0);

delete [] OriginalX0;

delete [] OriginalY0;

delete [] StandardX0;

delete [] StandardY0;

return 0;

}

由于上传图片太麻烦，这里给出全文下载

http://download1.csdn.net/down3/20070601/01124436747.pdf