metropolis算法

20世纪十大算法本世纪初，美国物理学会(American Institute of Physics)和IEEE计算机社团(IEEE Computer Society)的一本联合刊物《科学与工程中的计算》发表了由田纳西大学的Jack Dongarra和橡树岭国家实验室的Francis Sullivan联名撰写的“世纪十大算法”一文，该文“试图整理出在20世纪对科学和工程领域的发展产生最大影响力的十大算法”。

作者苦于“任何选择都将是充满争议的，因为实在是没有最好的算法”，他们只好用编年顺序依次列出了这十项算法领域人类智慧的巅峰之作——给出了一份没有排名的算法排行榜。

有趣的是，该期杂志还专门邀请了这些算法相关领域的“大拿”为这十大算法撰写十篇综述文章，实在是蔚为壮观。

本文的目的，便是要带领读者走马观花，一同回顾当年这一算法界的盛举。

1946蒙特卡洛方法在广场上画一个边长一米的正方形，在正方形内部随意用粉笔画一个不规则的形状，呃，能帮我算算这个不规则图形的面积么？蒙特卡洛(Monte Carlo)方法便是解决这个问题的巧妙方法:随机向该正方形内扔N（N是一个很大的自然数）个黄豆，随后数数有多少个黄豆在这个不规则几何形状内部，比如说有M个:那么，这个奇怪形状的面积便近似于M/N，N越大，算出来的值便越精确。

别小看这个数黄豆的笨办法，大到国家的民意测验，小到中子的移动轨迹，从金融市场的风险分析，到军事演习的沙盘推演，蒙特卡洛方法无处不在背后发挥着它的神奇威力。

蒙特卡洛方法由美国拉斯阿莫斯国家实验室的三位科学家John von Neumann（看清楚了，这位可是冯诺伊曼同志！）,Stan Ulam和Nick Metropolis共同发明。

就其本质而言，蒙特卡洛方法是用类似于物理实验的近似方法求解问题，它的魔力在于，对于那些规模极大的问题，求解难度随着问题的维数（自变量个数）的增加呈指数级别增长，出现所谓的“维数的灾难”（Course of Dimensionality）。

自适应采样多通路Metropolis算法贺怀清; 湛少胜; 刘浩翰【期刊名称】《《计算机工程与设计》》【年(卷),期】2019(040)010【总页数】9页(P2835-2842,2860)【关键词】模糊度; 多通路Metropolis; 八邻域自适应采样; 马尔可夫链; 采样权重【作者】贺怀清; 湛少胜; 刘浩翰【作者单位】中国民航大学计算机科学与技术学院天津300300【正文语种】中文【中图分类】TP391.410 引言采用蒙特卡罗基本框架且基于物理学原理的光照传播计算[1-3]旨在建立连接光源与视点之间的采样路径。

典型的算法如BDPT算法[4]不能完美地适应不同路径的采样，于是Eric Veach等将蕴涵“重用”思想的Metropolis采样引入蒙特卡罗光照计算领域[5]，称为Metropolis光照传播算法。

在此基础上Kelemen等提出了PSSMLT[6]算法，在基础样本空间内使用随机扰动的突变策略获取新的采样路径，以解决采样局部性问题。

Hachisuka等为改进PSSMLT算法路径探索开销过大问题，结合多重重要性采样思想提出了MMLT[7]算法。

曾笮等[8]通过对路径进行空间排序改进MMLT算法。

以上算法并未对路径样本对像素合成的影响做评价。

徐庆等[9]利用模糊不确定性评价像素内噪声水平；Sbert M等[10]结合像素内和像素间信息提出自适应采样标准；陈帅等[11]通过八邻域面片采样方式改进PSSMLT算法；贾洁等[12]采用方差过滤优化MMLT算法预采样阶段的初始样本。

本文借鉴文献[9]采用模糊不确定性评价像素内样本噪声水平的思想，同时借鉴文献[11]八邻域采样的思想，但本文只在八邻域位置中做筛选不进行同时采样。

在此基础上参考文献[10]兼顾像素内和像素间信息，引入模糊度作为像素噪声水平的评价标准，提出改进的MMLT自适应采样算法，达到渲染后期自适应地引导样本分布的效果，进而提升渲染效果。

模拟退⽕算法模拟退⽕（SA）物理过程由以下三个部分组成1.加温过程问题的初始解2.等温过程对应算法的Metropolis抽样的过程3.冷却过程控制参数的下降默认的模拟退⽕是⼀个求最⼩值的过程，其中Metropolis准则是SA算法收敛于全局最优解的关键所在，Metropolis准则以⼀定的概率接受恶化解，这样就使算法跳离局部最优的陷进1.模拟退⽕算法求解⼀元函数最值问题使⽤simulannealbnd - Simulated annealing algorithm⼯具箱求y=sin(10*pi*x)./x;在[1,2]的最值下图是⽤画图法求出最值的x=1:0.01:2;y=sin(10*pi*x)./x;figurehold onplot(x,y,'linewidth',1.5);ylim([-1.5,1.5]);xlabel('x');ylabel('y');title('y=sin(10*\pi*x)/x');[maxVal,maxIndex]=max(y);plot(x(maxIndex),maxVal,'r*');text(x(maxIndex),maxVal,{['x:' num2str(x(maxIndex))],['y:' num2str(maxVal)]});[minVal,minIndex]=min(y);plot(x(minIndex),minVal,'ro');text(x(minIndex),minVal,{['x:' num2str(x(minIndex))],['y:' num2str(minVal)]});hold off;⽤模拟退⽕⼯具箱来找最值求最⼩值function fitness=fitnessfun(x)fitness=sin(10*pi*x)./x;end求最⼤值function fitness=fitnessfun(x)fitness=-sin(10*pi*x)./x;endOptimization running.Objective function value: -0.9527670052175917Maximum number of iterations exceeded: increase options.MaxIterations.⽤⼯具箱求得的最⼤值为0.95276700521759172.⼆元函数优化[x,y]=meshgrid(-5:0.1:5,-5:0.1:5);z=x.^2+y.^2-10*cos(2*pi*x)-10*cos(2*pi*y)+20;figuremesh(x,y,z);hold onxlabel('x');ylabel('y');zlabel('z');title('z=x^2+y^2-10*cos(2*\pi*x)-10*cos(2*\pi*y)+20');maxVal=max(z(:));[maxIndexX,maxIndexY]=find(z==maxVal);%返回z==maxVal时，x和y的索引for i=1:length(maxIndexX)plot3(x(maxIndexX(i),maxIndexY(i)),y(maxIndexX(i),maxIndexY(i)),maxVal,'r*');text(x(maxIndexX(i),maxIndexY(i)),y(maxIndexX(i),maxIndexY(i)),maxVal,{['x:' num2str(x(maxIndexX(i)))] ['y:' num2str(y(maxIndexY(i)))] ['z:' num2str(maxVal)] }); endhold off;function fitness=fitnessfun(x)fitness=-(x(1).^2+x(2).^2-10*cos(2*pi*x(1))-10*cos(2*pi*x(2))+20);endOptimization running.Objective function value: -80.50038894455415Maximum number of iterations exceeded: increase options.MaxIterations.找到的最⼤值：80.500388944554153.解TSP问题（⽤的数据和前⼏天⽤遗传算法写TSP问题的数据⼀致，但是结果⽐遗传算法算出来效果差很多，不知道是不是我写错了，怀疑⼈⽣_(:з」∠)_中。

基于Metropolis-Hastings算法的α稳定分布参数估计马洪斌;马岩;杨春梅;沈锋【摘要】To solve the problem of parameter estimation of α-stable distributions, a method based on Markov chain Monte Carlo (MCMC) is proposed in this paper. A hierarchical model in order to make inference on the parameters was constructed based on Bayesian theorem, and the Metropolis-Hastings (M-H) sampling algorithm was used to generate Markov Chain. Under Bayesian framework, all the unknown parameters were regarded as random variables that can be simultaneously estimated using posterior distribution. The procedure of parameters updating was shown and the formula of acceptance rate was derived. Theoretic analysis and computer simulations indicate that the method performs well in estimating all the four parameters of arbitrary a-stable distributions.%针对α稳定分布参数估计问题,提出了一种基于MCMC动态模拟的参数估计方法.该方法根据贝叶斯理论建立在α稳定分布层次模型的基础上,利用Metropolis-Hastings抽样方法生成Mark-OV链,在贝叶斯框架下将所有待估计参数视为随机变量,利用后验分布实现稳定分布参数的同时估计,给出了新方法的迭代更新过程,并推导了接受概率的计算公式.理论分析和仿真结果表明,该方法能准确地估计出α稳定分布的4个参数,实现了任意对称或非对称α稳定分布的参数估计.【期刊名称】《电机与控制学报》【年(卷),期】2012(016)012【总页数】5页(P94-98)【关键词】α稳定分布;参数估计;MCMC;Metropolis-Hastings算法;贝叶斯推断【作者】马洪斌;马岩;杨春梅;沈锋【作者单位】哈尔滨理工大学机械动力工程学院,黑龙江哈尔滨150080;东北林业大学研究生院,黑龙江哈尔滨150040;东北林业大学机电工程学院,黑龙江哈尔滨150040;东北林业大学机电工程学院,黑龙江哈尔滨150040;哈尔滨工程大学自动化学院,黑龙江哈尔滨150001【正文语种】中文【中图分类】TP2710 引言在传统的信号处理中，由于高斯模型可以很好地描述许多信号噪声并且算法结构简单，通常假定接收机接收的信号噪声为高斯分布白噪声[1]，然而，在实际的情况中所遇到的许多信号或噪声往往具有显著地尖峰脉冲特性，例如水声信号、低频大气噪声、生物医学信号和金融数据等[2－4]。

马尔可夫链蒙特卡洛方法简介蒙特卡洛方法是一种通过随机抽样来解决问题的数值计算方法。

而在蒙特卡洛方法中，马尔可夫链蒙特卡洛方法（Markov Chain Monte Carlo, MCMC）是一种重要的技术，它可以用于求解很多实际问题，比如概率分布的估计、贝叶斯统计推断等。

本文将对马尔可夫链蒙特卡洛方法进行简要介绍。

1. 马尔可夫链马尔可夫链是指一个具有马尔可夫性质的随机过程。

所谓马尔可夫性质是指一个系统在给定当前状态下，未来的状态只与当前状态有关，而与过去状态无关。

换句话说，马尔可夫链的未来状态只取决于当前状态，而与过去状态无关。

这种性质使得马尔可夫链在模拟复杂系统时非常有用。

2. 马尔可夫链蒙特卡洛方法在蒙特卡洛方法中，马尔可夫链蒙特卡洛方法是通过构造一个马尔可夫链，使得该链的平稳分布恰好是我们要求的概率分布。

通过对该马尔可夫链进行随机抽样，最终可以得到与平稳分布一致的样本，从而对概率分布进行估计。

3. Metropolis-Hastings算法Metropolis-Hastings算法是一种常用的马尔可夫链蒙特卡洛方法。

其基本思想是通过一系列状态转移来构造一个满足平稳分布的马尔可夫链。

具体而言，算法首先随机初始化一个状态，然后通过一定的转移规则来进行状态转移。

在每次状态转移后，我们都根据一定的准则来接受或者拒绝转移，以保证最终的样本满足平稳分布。

4. Gibbs采样Gibbs采样是一种特殊的Metropolis-Hastings算法。

它适用于高维参数的分布估计问题。

在Gibbs采样中，我们将多维参数分解为多个条件分布，然后通过依次对每个条件分布进行抽样来得到最终的样本。

Gibbs采样在贝叶斯统计推断等领域有着广泛的应用。

5. 贝叶斯统计推断马尔可夫链蒙特卡洛方法在贝叶斯统计推断中有着重要的应用。

在贝叶斯统计中，我们往往需要对参数的后验分布进行估计。

而马尔可夫链蒙特卡洛方法可以通过对后验分布进行抽样来进行估计，从而得到参数的后验分布的近似值。

metropolis准则和物理退化过程全文共四篇示例，供读者参考第一篇示例：《metropolis准则和物理退化过程》引言：从古至今，人类一直在探索自然界的奥秘，物理学作为自然科学的一个重要分支，帮助我们理解世界的运行规律。

而在物理学中，metropolis准则是一个重要的概念，在研究物理退化过程中起到了关键作用。

本文将分析metropolis准则和物理退化过程之间的关系，帮助我们更好地理解物质世界的变化。

一、metropolis准则的基本概念metropolis准则最早由意大利数学家Metropolis提出，是一个关于随机模拟算法中接受与拒绝的准则。

在物理模拟中，我们经常需要通过随机采样来模拟系统的运行情况，而metropolis准则则帮助我们判断是否接受新的状态。

其基本原理是：如果新状态的能量小于当前状态的能量，则接受新状态；如果新状态的能量大于当前状态的能量，则按照一定的概率接受新状态。

metropolis准则的应用在物理模拟中非常广泛，可帮助我们高效地模拟复杂的系统。

二、物理退化过程的定义物理退化过程是指物质系统在外部作用下从一个有序状态变为一个无序状态的过程。

热力学第二定律告诉我们，自然界中所有的系统都趋向于熵的增加，即系统的有序度会不断下降。

这种有序度的下降就是物理退化过程。

比如，水从高处向低处流动，热传导导致温度均匀化等过程都是物理退化过程的例子。

三、metropolis准则和物理退化过程的关系在物理学中，metropolis准则被广泛应用于模拟系统的状态变化。

而在物理退化过程中，系统的状态也会不断变化，有序度逐渐减小。

可以说，物理退化过程也可以看作是系统状态的一个转移过程，而metropolis准则则可以帮助我们理解系统在状态之间的转移规律。

当系统在有序状态和无序状态之间转移时，我们可以利用metropolis准则来判断系统是否接受新的状态，从而更好地理解物理退化过程。

例如，在研究热平衡下系统的状态变化时，可以通过metropolis准则来模拟系统的状态转移，帮助我们理解系统的行为。

马尔可夫链蒙特卡洛方法简介马尔可夫链蒙特卡洛方法是一种基于随机抽样的数值计算方法，适用于求解复杂的概率和统计问题。

它的核心思想是利用马尔可夫链的收敛性质，通过随机抽样来模拟目标分布，并利用大数定律得到概率和统计量的近似解。

本文将介绍马尔可夫链蒙特卡洛方法的基本原理、应用领域和一些典型算法。

基本原理马尔可夫链蒙特卡洛方法的基本原理是基于马尔可夫链的收敛性质。

马尔可夫链是一种具有马尔可夫性质的随机过程，即下一时刻的状态只依赖于当前时刻的状态，而与之前的状态无关。

这种特性使得马尔可夫链具有收敛到平稳分布的性质，即当经过足够长的时间后，链的状态会趋向于一个固定的分布。

马尔可夫链蒙特卡洛方法利用马尔可夫链的收敛性质，通过从某一初始状态出发，经过多次状态转移后，得到一个服从目标分布的样本。

然后利用这些样本来估计目标分布的统计特性，如均值、方差、分位数等。

当样本量足够大时，根据大数定律，这些估计值会逼近真实值。

应用领域马尔可夫链蒙特卡洛方法在概率和统计领域有着广泛的应用。

其中，最为典型的应用就是概率分布的抽样和统计推断。

在贝叶斯统计中，常常需要对后验分布进行抽样，而马尔可夫链蒙特卡洛方法正是一种有效的抽样工具。

此外，在金融工程、统计物理、机器学习等领域，马尔可夫链蒙特卡洛方法也得到了广泛的应用。

除了概率和统计领域，马尔可夫链蒙特卡洛方法还被应用于优化问题的求解。

例如，模拟退火算法和遗传算法就是基于马尔可夫链蒙特卡洛方法的一种优化算法。

这些算法通过模拟随机状态的转移，逐步搜索最优解，对于复杂的优化问题有着良好的表现。

典型算法马尔可夫链蒙特卡洛方法有许多典型的算法，其中最为著名的包括Metropolis-Hastings算法和Gibbs抽样算法。

Metropolis-Hastings算法是一种基础的马尔可夫链蒙特卡洛方法，通过接受-拒绝的原则，实现对目标分布的抽样。

Gibbs抽样算法则是一种特殊的Metropolis-Hastings算法，适用于多维分布的抽样问题，它利用条件概率的性质，实现对联合分布的抽样。

三种启发式算法在钢铁排产系统中的对比摘要：钢铁的生产工序繁多复杂，合理安排好生产过程将会极大地提高生产效率，排产系统在其中就发挥着十分重要的作用，而一个高效合理的排产系统的构建离不开合适算法的使用。

本文将会对目前应用于排产系统的一些算法进行简单介绍。

关键词：排产系统；启发式算法；逆推算法；模拟退火算法；遗传算法一、引言全球范围对智能制造技术不断创新和扩大应用,对制造业提出了进一步提升内部信息化、数字化、智能化管理水平的要求，国内制造业目前的发展趋势，旨在降低库存水平，以减少成本，同时仍然能够做回应，以更短的交货时间从而满足客户的需求。

这要求企业优化生产操作，减少或消除非增值性的作业，如设置或等待时间，这也正是排产系统要解决的问题。

就钢铁制造业而言，钢铁生产过程涉及多台转炉、多种精炼工序、多台连铸机以及天车，这些生产工序之间的协作就会影响到不同钢种的生产效率和质量，这是一个多工序、多平台、离散和连续相结合的生产过程，对生产的安排提出了很高的要求[1]。

一个高效合理的排产系统将会极大地提高钢铁的生产效率和质量，而算法作为排产系统的核心内容，也影响着整个排产系统的功能。

二、排产系统排产系统是计划安排生产的系统，而计划排产属于车间调度问题，车间调度问题通常定义如下：在一定的约束条件下，把有限的资源在时间上按照一定的顺序分配给若干个任务，以满足或优化一个或多个性能指标[2]。

钢铁企业信息系统主要分为四个等级：L1基础自动化级、L2过程控制级、L3制造执行系统（MES）、L4公司管理级（ERP）。

排产系统是制造执行系统的核心部分，也是钢铁企业信息系统的重要组成部分，统领信息化系统的稳定运行，使生产产能最大化和各工序负载实现均衡，其运行质量的高低会直接影响到整个炼钢过程的生产效率。

基于不锈钢的调度主要涉及AOD冶炼、LF/VOD以及连铸的生产作业计划，其主要目标就是在尽可能满足当前生产并保证连铸机连浇的条件下，在接收上级的浇次计划之后，利用算法计算出合理的生产流程，通过下达生产指令，安排各炉次在各工序上的加工设备、加工时间和加工顺序，以获得生产作业计划的总流程时间或生产成本最优化。

一、1946 蒙特卡洛方法[1946: John von Neumann, Stan Ulam, and Nick Metropolis, all at the Los Alamos Scientific Labor atory, cook up the Metropolis algorithm, also known as the Monte Carlo method.]1946年，美国拉斯阿莫斯国家实验室的三位科学家John von Neumann,Stan Ulam 和 Nick Metrop olis共同发明，被称为蒙特卡洛方法。

它的具体定义是：在广场上画一个边长一米的正方形，在正方形内部随意用粉笔画一个不规则的形状，现在要计算这个不规则图形的面积，怎么计算列?蒙特卡洛(Monte Carlo)方法告诉我们，均匀的向该正方形内撒N（N 是一个很大的自然数）个黄豆，随后数数有多少个黄豆在这个不规则几何形状内部，比如说有M个，那么，这个奇怪形状的面积便近似于M/N，N越大，算出来的值便越精确。

在这里我们要假定豆子都在一个平面上，相互之间没有重叠。

(撒黄豆只是一个比喻。

)蒙特卡洛方法可用于近似计算圆周率：让计算机每次随机生成两个0到1之间的数，看这两个实数是否在单位圆内。

生成一系列随机点，统计单位圆内的点数与总点数，内接圆面积和正方形面积之比为PI:4，PI为圆周率。

(多谢网友七里河蠢才指出:S内接圆：S正=PI：4。

具体，请看文下第99条评论。

十六日修正)，当随机点取得越多（但即使取10的9次方个随机点时，其结果也仅在前4位与圆周率吻合）时，其结果越接近于圆周率。

二、1947 单纯形法[1947: George Dantzig, at the RAND Corporation, creates the simplex method for linear program ming.]1947年，兰德公司的，Grorge Dantzig，发明了单纯形方法。

Package‘metropolis’October13,2022Title The Metropolis AlgorithmVersion0.1.8Date2020-09-21Author Alexander Keil[aut,cre]Maintainer Alexander Keil<*************>Description Learning and using the Metropolis algorithm forBayesianfitting of a generalized linear model.The package vignetteincludes examples of hand-coding a logistic model using severalvariants of the Metropolis algorithm.The package also contains Rfunctions for simulating posterior distributions of Bayesiangeneralized linear model parameters using guided,adaptive,guided-adaptive and random walk Metropolis algorithms.The random walk Metropolis algorithm was originally described in Metropolis et al(1953);<doi:10.1063/1.1699114>.License GPL(>=2)Depends coda,R(>=3.5.0)Imports statsSuggests knitr,markdownVignetteBuilder knitrEncoding UTF-8Language en-USLazyData trueRoxygenNote7.1.0NeedsCompilation noRepository CRANDate/Publication2020-09-2121:20:07UTC12as.mcmc.metropolis.samples R topics documented:as.mcmc.metropolis.samples (2)expit (3)logistic_ll (4)magfields (4)metropolis.control (5)metropolis_glm (5)normal_ll (7)plot.metropolis.samples (8)print.metropolis.samples (9)summary.metropolis.samples (9)Index11 as.mcmc.metropolis.samplesConvert glm_metropolis output to mcmc object from package codaDescriptionAllows use of useful functions from coda packageUsage##S3method for class metropolis.samplesas.mcmc(x,...)Argumentsx an object from the function"metropolis"...not usedDetailsTBAValueAn object of type"mcmc"from the coda packageexpit3 Examples##Not run:library("coda")dat=data.frame(y=rbinom(100,1,0.5),x1=runif(100),x2=runif(100))res=metropolis_glm(y~x1+x2,data=dat,family=binomial(),iter=10000,burnin=3000, adapt=TRUE,guided=TRUE,block=FALSE)res2=as.mcmc(res)summary(res2)##End(Not run)expit Inverse logit transformDescriptionInverse logit transformUsageexpit(mu)Argumentsmu log-oddsValuereturns a scalar or vector the same length as mu with values that are the inverse logit transform of muExampleslogodds=rnorm(10)expit(logodds)logodds=log(1.0)expit(logodds)4magfields logistic_ll logistic log likelihoodDescriptionlogistic log likelihoodUsagelogistic_ll(y,X,par)Argumentsy binary outcomeX design matrixpar vector of model coefficientsValuea scalar quantity proportional to a binomial likelihood with logistic parameterization,given y,X,andparmagfields A case control study of childhood leukemia and mag-neticfields from Savitz,Wachtel,Barnes,et al(1998)doi:R hrefhttps:///10.1093/oxfordjournals.aje.a11494310.1093/oxfordjournals.aje.a114943.DescriptionA case control study of childhood leukemia and magneticfields from Savitz,Wachtel,Barnes,et al(1998)doi:10.1093/oxfordjournals.aje.a114943.UsagemagfieldsFormatA data frame with234rows and2variables:y childhood leukemiax exposure to magneticfieldmetropolis.control5 metropolis.control metropolis.controlDescriptionmetropolis.controlUsagemetropolis.control(adapt.start=25,adapt.window=200,adapt.update=25,min.sigma=0.001,prop.sigma.start=1,scale=2.4)Argumentsadapt.start start adapting after this many iterations;set to iter+1to turn off adaptation adapt.window base acceptance rate on maximum of this many iterationsadapt.update frequency of adaptationmin.sigma minimum of the proposal distribution standard deviation(if set to zero,posterior may get stuck)prop.sigma.startstarting value,orfixed value for proposal distribution s standard deviation scale scale value for adaptation(how much should the posterior variance estimate be scaled by?).Scale/sqrt(p)is used in metropolis_glm function,and Gelmanet al.(2014,ISBN:9781584883883)recommend a scale of2.4@return Alist of parameters used infitting with the following named objects adapt.start,adapt.window,adapt.update,min.sigma,prop.sigma.start,scalemetropolis_glm Use the Metropolis Hastings algorithm to estimate Bayesian glm pa-rametersDescriptionThis function carries out the Metropolis algorithm.6metropolis_glmUsagemetropolis_glm(f,data,family=binomial(),iter=100,burnin=round(iter/2),pm=NULL,pv=NULL,chain=1,prop.sigma.start=0.1,inits=NULL,adaptive=TRUE,guided=FALSE,block=TRUE,saveproposal=FALSE,control=metropolis.control())Argumentsf an R style formula(e.g.y~x1+x2)data an R data frame containing the variables in ffamily R glm style family that determines model form:gaussian()or binomial()iter number of iterations after burnin to keepburnin number of iterations at the beginning to throw out(also used for adaptive phase) pm vector of prior means for normal prior on log(scale)(if applicable)and regres-sion coefficients(set to NULL to use uniform priors)pv vector of prior variances for normal prior on log(scale)(if applicable)and re-gression coefficients(set to NULL to use uniform priors)chain chain id(plan to deprecate)prop.sigma.startproposal distribution standard deviation(starting point if adapt=TRUE) inits NULL,a vector with length equal to number of parameters(intercept+x+scale ;gaussian()family only model only),or"glm"to set priors based on an MLEfit adaptive logical,should proposal distribution be adaptive?(TRUE usually gives better answers)guided logical,should the"guided"algorithm be used(TRUE usually gives better an-swers)block logical or a vector that sums to total number of parameters(e.g.if there are 4random variables in the model,including intercept,then block=c(1,3)willupdate the intercept separately from the other three parameters.)If TRUE,thenupdates each ing guide=TRUE with block as a vector isnot advisedsaveproposal(logical,default=FALSE)save the rejected proposals(block=TRUE only)?control parameters that controlfitting algorithm.See metropolis.control()normal_ll7 DetailsImplements the Metropolis algorithm,which allows user specified proposal distributions or im-plements an adaptive algorithm as described by Gelman et al.(2014,ISBN:9781584883883).This function also allows the"Guided"Metropolis algorithm of Gustafson(1998)doi:10.1023/A:1008880707168.Note that by default all parameters are estimated simulataneously via"block"sampling,but this default behavior can be changed with the"block"parameter.When using guided=TRUE, block should be set to FALSE.ValueAn object of type"metropolis.samples"which is a named list containing posterior MCMC samplesas well as somefitting information.Examplesdat=data.frame(y=rbinom(100,1,0.5),x1=runif(100),x2=runif(100))res=metropolis_glm(y~x1+x2,data=dat,family=binomial(),iter=1000,burnin=3000,adapt=TRUE,guided=TRUE,block=FALSE)ressummary(res)apply(res$parms,2,mean)glm(y~x1+x2,family=binomial(),data=dat)dat=data.frame(y=rnorm(100,1,0.5),x1=runif(100),x2=runif(100),x3=rpois(100,.2))res=metropolis_glm(y~x1+x2+factor(x3),data=dat,family=gaussian(),inits="glm",iter=10000,burnin=3000,adapt=TRUE,guide=TRUE,block=FALSE)apply(res$parms,2,mean)glm(y~x1+x2+factor(x3),family=gaussian(),data=dat)normal_ll Gaussian log likelihoodDescriptionGaussian log likelihoodUsagenormal_ll(y,X,par)Argumentsy binary outcomeX design matrixpar vector of gaussian scale parameter followed by model coefficientsValuea scalar quantity proportional to a normal likelihood with linear parameterization,given y,X,andparplot.metropolis.samplesPlot the output from the metropolis functionDescriptionThis function allows you to summarize output from the metropolis function.Usage##S3method for class metropolis.samplesplot(x,keepburn=FALSE,parms=NULL,...)Argumentsx the outputted object from the"metropolis_glm"functionkeepburn keep the burnin iterations in calculations(if adapt=TRUE,keepburn=TRUE parms names of parameters to plot(plots thefirst by default,if TRUE,plots all)...other arguments to plotDetailsTBAValueNoneExamplesdat=data.frame(y=rbinom(100,1,0.5),x1=runif(100),x2=runif(100))res=metropolis_glm(y~x1+x2,data=dat,family=binomial(),iter=10000,burnin=3000, adapt=TRUE,guided=TRUE,block=FALSE)plot(res)print.metropolis.samplesPrint a metropolis.samples objectDescriptionThis function allows you to summarize output from the"metropolis_glm"function.Usage##S3method for class metropolis.samplesprint(x,...)Argumentsx a"metropolis.samples"object from the function"metropolis_glm"...not used.DetailsNoneValueAn unmodified"metropolis.samples"object(invisibly)summary.metropolis.samplesSummarize a probability distribution from a Markov ChainDescriptionThis function allows you to summarize output from the metropolis function.Usage##S3method for class metropolis.samplessummary(object,keepburn=FALSE,...)Argumentsobject an object from the function"metropolis"keepburn keep the burnin iterations in calculations(if adapt=TRUE,keepburn=TRUE will yield potentially invalid summaries)...not used10summary.metropolis.samplesDetailsTBAValuereturns a list with the followingfields:nsamples:number of simulated samples sd:standard devia-tion of parameter distributions se:standard deviation of parameter distribution means ESS_parms: effective sample size of parameter distribution means postmean:posterior means and normal based 95%credible intervals postmedian:posterior medians and percentile based95%credible intervals postmode:posterior modes and highest posterior density based95%credible intervals Examplesdat=data.frame(y=rbinom(100,1,0.5),x1=runif(100),x2=runif(100))res=metropolis_glm(y~x1+x2,data=dat,family=binomial(),iter=10000,burnin=3000, adapt=TRUE,guided=TRUE,block=FALSE)summary(res)Index∗datasetsmagfields,4as.mcmc.metropolis.samples,2expit,3logistic_ll,4magfields,4metropolis.control,5metropolis_glm,5normal_ll,7plot.metropolis.samples,8print.metropolis.samples,9summary.metropolis.samples,911。

模拟退火算法一、模拟退火算法概念模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。

根据Metropolis准则，粒子在温度T 时趋于平衡的概率为e-ΔE/(kT)，其中E为温度T时的内能，ΔE为其改变量，k为Boltzmann 常数。

用固体退火模拟组合优化问题，将内能E模拟为目标函数值f，温度T演化成控制参数t，即得到解组合优化问题的模拟退火算法：由初始解i和控制参数初值t开始，对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代，并逐步衰减t值，算法终止时的当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。

退火过程由冷却进度表(Cooling Schedule)控制，包括控制参数的初值t及其衰减因子Δt、每个t值时的迭代次数L和停止条件S。

二、模拟退火算法的模型模拟退火算法可以分解为解空间、目标函数和初始解三部分。

模拟退火的基本思想:(1) 初始化：初始温度T(充分大)，初始解状态S(是算法迭代的起点)，每个T值的迭代次数L(2) 对k=1，……，L做第(3)至第6步：(3) 产生新解S′(4) 计算增量Δt′=C(S′)-C(S)，其中C(S)为评价函数(5) 若Δt′<0则接受S′作为新的当前解，否则以概率exp(-Δt′/T)接受S′作为新的当前解.(6) 如果满足终止条件则输出当前解作为最优解，结束程序。

终止条件通常取为连续若干个新解都没有被接受时终止算法。

(7) T逐渐减少，且T->0，然后转第2步。

算法对应动态演示图：模拟退火算法新解的产生和接受可分为如下四个步骤：第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。

联合Metropolis算法和MDL准则的密集多目标分辨算法李振兴;刘进忙;周政;郭相科;李延磊【摘要】针对单脉冲雷达处理同一分辨单元内密集多目标回波时收敛较慢的问题，提出了一种将 Metropolis算法和最小描述长度（MDL）准则相结合的密集多目标分辨算法。

构建 Metropolis 算法中的更新函数和迭代规则，促使待估参数的不同马尔可夫链间的融合；根据融合判定规则选取抽样样本，估计出对应的目标参数；利用 MDL 准则实现密集目标的准确分辨。

该算法不仅提升了收敛速度，而且具有较高的参数估计精度，提高了算法在多目标下的分辨性能。

仿真结果验证了所提算法的有效性和可行性。

%When multiple targets' echoes fall into the same resolution cell of monopulse radar , the echoes are distinguished at a slow rate of convergence using the available processing algorithm . Aiming at this problem , a new processing algorithm is presented , which combines the Metropolis algorithm with minimum description length ( MDL) criteria . Constructing the update function and the rule of iteration can improve the convergence of the different markov chains of parameters estimation . Then , by using the fusion decision rule , valid samples are selected to calculate the parameters . Lastly , the numbers of multiple unresolved targets are calculated by using minimum description length criteria . The proposed algorithm has not only a higher rate of convergence , but also a better resolution capability under multiple targets with a more accurate parametric estimation . The efficacy and feasibility of the method are verified by simulation .【期刊名称】《西安电子科技大学学报（自然科学版）》【年(卷),期】2014(000)003【总页数】7页(P174-180)【关键词】密集多目标;分辨算法;Metropolis算法;MDL 准则;马尔可夫链【作者】李振兴;刘进忙;周政;郭相科;李延磊【作者单位】空军工程大学防空反导学院，陕西西安 710051;空军工程大学防空反导学院，陕西西安 710051;空军工程大学防空反导学院，陕西西安 710051;空军工程大学防空反导学院，陕西西安 710051;空军工程大学防空反导学院，陕西西安 710051【正文语种】中文【中图分类】TN95雷达目标分辨常指对雷达目标进行检测及定位，即估计雷达目标在探测空间中的位置参数和目标数目.当多个目标位于雷达的同一分辨单元时，使用常规探测手段难以对其进行有效分辨及跟踪.如果雷达对这种密集多目标情况缺乏特殊的处理能力，而只是当作单目标以传统方法来处理时，将会带来一系列的严重后果[1－2].因此，对密集多目标的处理是当前雷达技术领域所面临的前沿课题和紧迫任务，而其中的分辨技术则是解决其他问题的前提和基础.虽然目前复杂的阵列信号处理技术及多波束数字形成技术能较好地解决多目标的分辨及定位，但是利用单脉冲雷达进行相应处理仍占据着十分重要的地位，对此也进行了很多的研究[3－9].文献[3－5]提出了针对SwerlingⅡ 型目标的广义似然比检测法，并提出了对两点源进行角估计的矩估计法.文献[6]针对两个不可分辨的瑞利目标提出了最大似然角估计方法.上述文献都假设多个目标回波滤波输出均位于同一个采样点上，在相邻的采样点上没有目标回波能量.但在实际情况中，常会在相邻采样点上存在回波的能量.针对这个问题，文献[7－8]提出了对目标回波滤波输出取两个采样点，利用最大似然估计法（MLE）对目标参数进行估计.文献[9]针对文献[7]提出的MLE存在的收敛时间过长、易陷入局部最优等缺点，提出利用传统Gibbs抽样法进行参数估计，取得了估计性能上的改善.但是文献[9]采用的Gibbs抽样法存在冗余抽样，即burnin区问题[10]，降低了算法的适用性.针对这一问题，笔者引入MCMC方法中的Metropolis算法[11]，并将其与MDL准则相结合，提出了一种新的密集多目标分辨方法.该方法的优点是：不仅提升了收敛速度，而且具有更高的参数估计精度，提高了算法在多目标下的分辨性能.由于算法的参数估计性能与初始值设置无关，因而具有良好的稳定性，实际应用方便.1 问题描述1.1 信号的数学模型文中使用的信号参数模型是沿用文献[4，8]建立的参数模型，即假设探测目标为SwerlingⅡ型目标，此时处于同一分辨单元内的每个目标回波能量至多出现在相邻两个采样点上.设单目标情况下接收机同相及正交支路输出的目标信号幅度分别为x和y，即x＝βcosφ，y＝βsinφ，其中，目标回波幅度β服从瑞利分布，φ服从[0，2π]范围的均匀分布.基于雷达信号处理的相关知识，可知固定载频下矩形脉冲信号经过接收机匹配滤波器输出后，雷达的单目标视频回波具有三角形包络[2，8].因此，不失一般性，假设雷达发射脉冲宽度为T，匹配滤波输出的采样率为1/T.令目标位置参数α＝ΔT T，则对应两采样点上回波的幅度分别为（1－α）x和αx.单脉冲雷达主要通过3个位置参数来对不同目标进行区分，即距离参数α、方位参数ηh和俯仰参数ηv[8].值得注意的是，这里ηh和ηv只代表角信息，并不是目标的真实角度，而是和天线、差天线在目标真实角处的增益的比值（其定义参见文献[7－8]）.因此，可利用这3个参数对单脉冲雷达的同相及正交支路的和波束以及差波束的方位和仰角输出进行建模.假定在雷达同一分辨单元内同时存在N个SwerlingⅡ型目标，令sik（m）、dhik（m）和dvik（m）分别代表由接收机同相支路获得的第m（m＝1，2，…，M）个目标回波，在第k（k＝1，2）个采样点处的和波束、方位差波束和俯仰差波束的输出（j＝1，2，…，N）个目标的3个待估参数.接收机同相支路输出的具体公式参见文献[7].同理，也可以类似得到正交支路输出的6组观测数据.需要指出的是，由于正交支路与同相支路的相似性，这里的M个实数脉冲等价于M/2个复数脉冲.在探测目标假设为SwerlingⅡ 型目标的背景下，目标j对应的第m个子脉冲输出的目标信号幅度xj（m）服从零均值的高斯分布[8].设xj（m）的方差＝1，即服从于标准正态分布.因此，同相支路的观测量为1.2 现有分辨算法分析分辨密集多目标的处理思想主要是：首先依据观测值，估计出目标的空间位置参数；在此基础上，利用信息论准则对目标数目进行检测，从而实现分辨目的.现有分辨算法主要是采用参数估计算法同MDL准则相结合来进行处理.所采用的参数估计算法主要分为两种类型：最大似然法估计（MLE）法[7－8]和MCMC参数估计法[9].1.2.1 最大似然估计法最大似然估计法主要是求取观测量Z的最大似然函数，来估计出其中的分辨参数.它的对数似然函数为其中，协方差矩阵R的表达式参见文献[7].根据最大似然估计原理，有式（4）中，必须同时满足以下条件：0≤αj≤1，－1≤ηhj≤1，－1≤ηvj≤1，j＝1，2，…，N.文献[7－8]采用了梯度投影法来进行寻优，为了提高算法的效率，还使用了Levenberg－Marquardt法，将Hessian矩阵的对角线元素标记为最速下降方向，从而寻找出最优值.由于该方法需要完成包括矩阵求导等一系列计算，尤其在多个目标的情况下，协方差矩阵R包含的未知参数增加，导致其算法的运算量急剧增加，严重影响了算法的实用性.1.2.2 MCMC参数估计法文献[9]针对 M LE法中的问题，提出采用 M CMC中的Gibbs抽样法进行处理.假设Ψ ＝ [x，α，ηh，ηv]T，为要估计的参数集，假定xj、αj、ηhj、ηvj 的先验概率分布都服从标准正态分布，其中，j＝1，2，…，N.首先，计算参数的后验概率分布公式.由于观测值z服从高斯分布，因而各参数的后验概率分布也服从高斯分布，即满足：其次，进行迭代计算，这里采用直接判断的方法来确定收敛.最后，求取参数的后验平均估计，作为目标参数的估计结果.通过仿真发现，文献[9]利用Gibbs抽样法，通过迭代来获得全局最优解，避免了大量复杂的矩阵运算，取得了较好的效果.但是仍存在下面3个方面问题：（1）初始值设置问题.文献[9]在计算各参数的后验概率分布时，由于参数的估计误差一般很小，从而降低了预测协方差的修正作用，使得算法收敛于真值的速率仍较慢.一旦初始值选择的误差较大时，会严重影响参数的估计精度.（2）冗余抽样问题.原Gibbs估计法无法确定burn－in时间，因此常利用大量的迭代运算来进行参数估计.这大大降低了算法的收敛速度.（3）参数状态空间的选择问题.通过前面的信号数学模型可知，每个目标包括3个待估参数（αj，ηhj，ηvj），它们彼此之间构成了统一整体.如果仅仅是进行单个参数的简单迭代，这无疑忽略了同一目标参数在参数空间内的相关性.因此，在参数估计过程中应将同一目标的3个参数综合起来考虑.2 基于Metropolis算法的参数估计及目标数估计由于Metropolis算法具有许多Gibbs采样器的特征，而且还回避了从复杂分布采样的问题，因而具有更加普遍的应用[11].在Metropolis算法的基础上，为了增加算法的效率，文中还充分吸收了文献[10]的马尔可夫链融合思想，即当大量的马尔可夫链的初始点是状态空间中的任意点，随着马尔可夫链中状态的增加，所有马尔可夫链通过相应处理可在某一个时间融合在一起.通过采用融合判决来确定马尔可夫链的稳定时间，可以避免初始值设置和冗余抽样问题，提高了算法的估计效率. 2.1 更新函数选择更新函数的选择是一个核心问题.更新函数可将马尔可夫链的下一个状态表示为马尔可夫链当前状态和建议分布生成状态之间的函数，即其中，Θ（t＋1）为马尔可夫链的下一个状态，Θ（t）为马尔可夫链的当前状态，Φ（t＋1）为建议分布产生的状态，fupdate（·，·）为更新函数，R（t＋1）为[0，1]之间的随机数.更新函数为其中，Θ为状态空间的当前状态，Φ为状态空间的候选状态.并有其中，fMAP为状态选择函数.为了保证融合的有效性，Metropolis算法要求初始状态应该覆盖状态空间的所有状态，但这会造成Metropolis算法难以应用.如果能够找到一个状态空间上的具有单调保偏序性质的接受概率函数，那么可用状态空间中任意两个状态为两条马尔可夫链的初始状态，以这两条马尔可夫链的融合判决来确定马尔可夫链的稳定时间，最后由稳定后的马尔可夫链完成参数估计.因此，确定p（Θ）为由于在计算更新函数时，当Θ比Φ更接近于真实状态子空间时，能够满足Θ和Φ之间的序关系，即若p（Θ）≤p（Φ），则有Θ≺Φ；否则，Φ≺Θ.因而，可以满足算法单调保偏性的要求.2.2 Metropolis算法步骤Metropolis算法的步骤如下：Step 1 随机产生两条马尔可夫链的初始参数向j＝1，2.Step 2 对于每条马尔可夫链，在每个时刻将产生两个候选状态.第1个候选状态Φ1是利用Gibbs法对上一时刻的参数向量抽样得到.第2个候选状态Φ2是在Φ1的基础上，通过建议分布产生的，这里确定在Φ1的局部状态空间中抽样.具体步骤如下：FORi＝1：N .建立以Φ1中第i个目标的3个参数值为均值，δ2为方差的高斯分布，从中抽样得到第i个目标的子状态Φ2_i.随着迭代次数的增加，可逐渐适当减小δ值，以确保迭代初期较早收敛到真值附近，迭代后期在真值附近能得到精度较高的估计值. 由式（8）得到判决函数γ.产生R～U[0，1]，U[0，1]为[0，1]之间的均匀分布.由式（9）确定采用Φ1 或者Φ2_i.重复这一过程，直到所有N个目标循环结束.END.Step 3 进行融合判决.若 m a x｛－｝＜ε，则任取一个Θij（j＝1，2）作为初始状态，应用各参数的具体后验分布概率，即式（5）来形成平稳的含有q节点的马尔可夫链，并对其求期望，获得目标方位向量的（Markov）链的第i个目标的第m个元素.2.3 目标数目估计上述多目标参数估计是在目标数已知的情况下获得的，但实际中目标的个数是未知的.结合文中的应用背景，确定采用MDL准则来估计目标数目[8，12].MDL准则思想为：对于给定一组数据和一组关于待定参数的数据统计模型，最佳的模型应该对数据提供最小的描述长度.最小描述长度可定义为其中，n为待估计的目标个数，k为每个目标待估计的参数个数，Q为观测数据个数，LL（）n为估计参数的对数似然函数.最终确定的目标个数应该满足文献[8]经过严格的数学推导，确定两个采样点间有效分辨的目标数目至多为5个.因此，确定n的迭代区间为[1，5].3 仿真实验类似文献[7]的仿真场景设置，文中共设置了5个目标的仿真参数，包括分辨单元内距离参数α、方位参数ηh和俯仰参数ηv，设置的参数都是经过归一化处理的，具体请参考文献[7].这里定义信噪比RSN＝单目标回波幅度xj（m）的方差，为方便计算，这里设为1.3.1 参数估计仿真假设回波子脉冲积累数为20个，迭代次数为100次（经融合判决，算法在第50次时已经收敛），对收敛后的数据每隔4个进行抽取，对获得的数据进行遍历求均，求得参数估计.图1显示了信噪比为20dB，目标数为3时，文中算法对两条方位估计的马尔可夫链的融合效果.可以看到，经过有限次迭代后，两条马尔可夫链在信号真实值附近融合在一起，也验证了文中算法能有效减少迭代次数.图1 马尔可夫链融合的效果（N＝3）图2显示了信噪比为20dB时，进行50次蒙特卡罗实验后，Metropolis算法（N ＝3，M＝50，λ＝1×106，δ＝0.1，ε＝0.1）的角估计性能.从对比结果可知，Metropolis算法的收敛效果要优于Gibbs算法.图2 两种算法对3个目标的方位和俯仰参数的估计效果（M＝20）图3显示了在10～30dB下，两种方位对目标数分别为3和5的参数估计的均方差根（RMSE）比较.这里结合文中的仿真场景，首先分别计算不同目标的相同参数的均方差根值，然后对同一参数的多个均方差根值求取期望，从而得到对该参数的最终均方差根值.图3中实线为Metropolis算法估计结果（为了不遮挡图中结果显示，在图中表示为“PS－”），虚线为Gibbs算法估计结果.对比图3（a）和图3（b）可知，随着目标数目的增加，目标角估计误差的均方根也增大.这主要是由于随着目标数的增加，待估参数也随着增多，在估计过程中会对估计精度造成影响.此外，文中算法在角估计精度方面要高于文献[9]所提的Gibbs抽样估计法，尤其随着信噪比的增大，其提高的趋势更加明显.表1给出了信噪比为20dB时，两种算法在不同目标数情况下，平均运算总次数以及每次迭代运算的平均运算时间对比.文中使用的平台是3.19GHz Pentium Dual－Core CPU，matlab（R2010b）.从表1可知，当处理的积累脉冲增多，或真实目标数增多时，两种算法的每次迭代的平均运行时间都会增加.而当条件相同时，Gibbs估计法的每次迭代的平均运行时间要小于文中算法的，这是由于后者增加了判决函数γ的计算过程.但是，文中算法的平均运算总次数要远小于Gibbs估计法的，因而能够在较少的迭代次数下收敛，而Gibbs抽样收敛速率较慢，往往需要大量的迭代运算.因此，文中算法的整体收敛速度高于Gibbs估计法的，并且文中算法受初始值的影响比Gibbs估计法的更小，在整体上，文中算法估计性能优于Gibbs估计法的估计性能.图3 两种算法对参数估计的均方差根比较（M＝20）表1 两种算法平均运算总次数/单次平均运算时间对比 s真实目标数M＝10 M＝20 M＝30M＝10 M＝20 M＝30 1 23/0.004 9 24/0.012 6 26/0.024 897/0.001 7 103/0.004 8 109/0.009 1 2 25/0.009 4 29/0.022 8 32/0.046 2 229/0.003 6 247/0.010 2 283/0.021 2 3 49/0.014 1 52/0.034 3 54/0.068 1 316/0.006 1 348/0.014 9 361/0.030 9 4 60/0.019 7 64/0.046 2 69/0.091 3 481/0.008 3 503/0.019 2 523/0.038 5 5 73/0.029 8 79/0.062 4 83/0.117 4 554/0.010 7 582/0.024 3 624/0.051 9 Metropolis估计法Gibbs估计法3.2 目标数估计仿真下面进行目标数的检测仿真.在前述目标参数估计的基础上，回波子脉冲积累数为20个，通过MDL准则估计目标的数目.这里定义的正确检测率，是指通过仿真获得的正确目标数的次数与仿真次数的比值.文中进行蒙特卡罗仿真150次，获得仿真结果如图4所示.由图4可知，正确检测率与信噪比和真实目标数目是相关的.随着目标数目的增多，其相应的正确检测率降低，而当真实目标数目不变时，正确检测率会随着信噪比的增大而变大.由表2可看出，当真实目标数较多时，文中基于Metropolis算法与MDL准则结合获得的正确检测率要高于文献[8－9]算法所获得的检测率.表2 信噪比为25dB时3种算法的正确检测率对比?图4 文中算法的密集多目标的正确检测率4 结束语基于多目标分辨和参数估计信号模型，提出了一种基于Metropolis算法的目标参数估计方法，并利用MDL准则进行目标数目的估计，从而达到对多个目标源的分辨.通过仿真实验，验证了文中算法的估计精度更高，并且能有效减少算法的收敛时间，同时利用MDL准则可以在真实目标数较多时，得到更高的目标正确检测率.因此，文中算法对多个目标源具有更优的分辨性能.参考文献：[1] 赵锋，毕莉.一种空间不可分辨多目标存在性检测的新方法[J].电子学报，2010，38（10）：2258－2263.Zhao Feng，Bi Li.A New Method for Detecting the Presence of Multiple Unresolved Targets[J].Acta Electronica Sinica，2010，38（10）：2258－2263.[2] 张直中.雷达信号的选择与处理[M].北京：国防工业出版社，1979：190－191.[3] Blair W D，Brandt－Pearce M.Statistical Description of Monopulse Parameters for Tracking Rayleigh Targets[J].IEEE Transactions on Aerospace and Electronic Systems，1998，34（2）：597－611.[4] Blair W D，Brandt－Pearce M.Unresolved Rayleigh Target Detection Using Monopulse Measurements [J].IEEE Transactions on Aerospace and Electronic Systems，1998，34（2）：543－552.[5] Blair W D，Brandt－Pearce M.Monopulse DOA Estimation of Two Unresolved Rayleigh Targets[J].IEEE Transactions on Aerospace and Electronic Systems，2001，37（2）：452－469.[6] Sinha A，Kirubarajan T，Bar－Shalom Y.Maximum Likelihood Angle Extractor of Two Closely Spaced Targets[J].IEEE Transactions on Aerospace and Electronic Systems，2002，38（1）：183－203.[7] Zhang X，Willett P K，Bar－Shalom Y.Monopulse Radar Detection and Localization of Multiple Unresolved Targets Via Joint Bin Processing[J].IEEE Transactions on Signal Processing，2005，53（4）：1225－1236.[8] Zhang X，Willett P，Bar－Shalom Y.Detection and Localization of Multiple Unresolved Extended Targets Via Monopulse Radar Signal Processing[J].IEEE Transactions on Aerospace and Electronic Systems，2009，45（2）：455－470.[9] 李朝伟，王宏强，黎湘，等.一种基于Gibbs抽样的多信号源分辨方法[J].信号处理，2006，22（4）：449－453.Li Chaowei，Wang Hongqiang，Li Xiang，et al.A Method of Distinguishing Multi－sources Based on the Gibbs Sampling[J].Signal Processing，2006，22（4）：449－453.[10] Djuric P M，Huang Y，Ghirmai T.Perfect Sampling：a Review and Applications to Signal Processing[J].IEEE Transactions on Signal Processing，2002，50（2）：345－356.[11] 王进玲，曾声奎，马纪明.基于MCMC的模糊自适应重要性抽样法[J].系统工程与电子技术，2012，34（2）：317－322.Wang Jinling，Zeng Shengkui，Ma Jiming.Fuzzy Adaptive Importance Sampling Method Based onMCMC[J].Systems Engineering and Electronics，2012，34（2）：317－322.[12] 杨志伟，贺顺，张娟，等.利用数据最小冗余拟合的信源数估计方法[J].西安电子科技大学学报，2012，39（4）：52－56.Yang Zhiwei，He Shun，Zhang Juan，et al.Estimation of the Number of Sources Via Minimal Redundancy Fitting[J].Journal of Xidian University，2012，39（4）：52－56.。

Metropolis Monte Carlo方法
个自旋位形中，任何一个位形S 在正则系模型中的宏观统计量
∑−a
a
S B )
)(1)(S H e
Z S w −=
∑a
a S H
Dr. Ernest Ising,
1900-1998
)
伊辛模型的蒙特卡罗模拟基本步骤
E new <E old E new ≥E old
不接受这次自旋翻转
建立晶格，并定义与所考察系综相对应的自旋和哈
密顿量，假定计数从n =1开始，同时设定n 0和n m a x
自旋翻转
接受这次自旋翻转计算变量A n 的值，并存储n >n 0的每一步A n 值
计算平均值，输出结果
计算p =e x p [-ΔH /(k B T )]
产生随机数，0<z <1
z>p z ≤p n=n+1MC 抽样下宏观统计量的表达
N H
E =∑=i
i
S N M 1
∑∑−−=),(,j i a
a
j i j i S B S S J H •伊辛模型系统的哈密顿量
•平均每个格点的磁场能•平均每个格点的磁化强度•平均定场比热容)
(1
223H H T
Nk C B B −=
A 2D lattice MC code
•Origin 7.5
•\samples\programming\Ising model.opj
结果举例—E-T曲线
磁畴结构。

第 22卷第 8期2023年 8月Vol.22 No.8Aug.2023软件导刊Software Guide融合反向学习与Metropolis准则求解TSP的遗传算法叶梓萌，张大斌（广东白云学院大数据与计算机学院，广东广州 510450）摘要：针对传统遗传算法（GA）在求解巡回旅行商（TSP）问题时，因种群多样性、后期局部搜索寻优能力不足所导致全局搜索范围受限、易陷入早熟现象与收敛速度较慢的问题，基于种群逆转遗传算法（EROGA）融合反向学习（OBL）思想，模拟退火算法（SA）的Metropolis准则与现实精英学习理念，提出一种改进的遗传算法OBLGSAA。

首先在生成初始种群环节采用反向学习方式，以提升最优解的精度与收敛速度；然后采用Metropolis准则改进交叉、变异算子，以提升算法的局部搜索能力；最后引入现实精英学习理念，通过贪心轮转学习机制进一步提升GA的局部搜索能力。

在多种巡回旅行商数据集的仿真实验结果表明，OBLGSAA能有效改善GA种群的多样性，使算法不易陷入早熟收敛，并在收敛性能与求解精度上相较于原始EROGA更优，验证了算法的有效性与优越性。

关键词：反向学习；Metropolis准则；遗传算法；旅行商问题DOI：10.11907/rjdk.231452开放科学（资源服务）标识码（OSID）：中图分类号：TP18 文献标识码：A文章编号：1672-7800（2023）008-0104-07Genetic Algorithm Combining Opposition-based Learning andMetropolis Criterion for TSPYE Zimeng， ZHANG Dabin（College of Big Data and Computer， Guangdong Baiyun University， Guangzhou 510450， China）Abstract：When traditional genetic algorithm （GA） is used to solve traveling salesman problem （TSP）， the global search scope is limited，and it is easy to fall into premature phenomenon and slow Rate of convergence due to the lack of population diversity and local search ability in the later stage. Based on the population inversion genetic algorithm （EROGA）， which combines the idea of reverse learning （OBL）， the Me⁃tropolis principle of Simulated annealing algorithm （SA） and the idea of real elite learning， an improved genetic algorithm OBLGSAA is pro⁃posed. First， reverse learning is adopted in the process of generating initial population to improve the accuracy and Rate of convergence of the optimal solution； Then， the Metropolis criterion is used to improve the crossover and mutation operators to enhance the local search ability of the algorithm； Finally， the concept of real-life elite learning is introduced to further enhance GA's local search ability through a greedy rota⁃tion learning mechanism. The simulation results on various traveling salesman datasets show that OBLGSAA can effectively improve the diver⁃sity of GA population， making the algorithm less prone to premature convergence， and outperforming the original EROGA in terms of conver⁃gence performance and solution accuracy， verifying the effectiveness and superiority of the algorithm.Key Words：opposition-based learning； Metropolis criterion； genetic algorithm； traveling salesman problem0 引言求解旅行商（Traveling Salesman Problem ，TSP）［1］本质是求解NP-hard［2］组合优化问题，映射到现实社会具有实际意义，已被广泛应用于物流运输［3］、旅行路径规划［4］、电网规划［5］、生产作业［6］等行业领域。